StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ua1.h
1 
2 int DEBUG=0;
3 
4 TGraphErrors* ua1Data200();
5 
6 TF1* ua1Fit130(scale=1);
7 
8 void
9 scale(TGraphAsymmErrors* graph,double scale);
10 
11 void
12 scale(TGraphErrors* graph,double scale);
13 
14 // divide by a float (e.g. TAA)
15 void
16 divide(TGraphASymmErrors* graph, double div, double divErr=0);
17 
18 // sets default
19 // type==1 : lower limit of ua1
20 TGraphAsymmErrors*
21 makeUA1(const TGraphAsymmErrors* graph, int type=0);
22 
23 TGraphAsymmErrors*
24 divide(TGraphAsymmErrors* graphA,TGraphAsymmErrors* graphB);
25 
26 // scale is the overall scale of the ratio
27 // type is which ua1 spectrum to create
28 TGraphAsymmErrors*
29 makeUA1Ratio(TGraphAsymmErrors* gSTAR,
30  int type, double scale);
31 
32 // divide ua1 ratio by TAA
33 TGraphAsymmErrors*
34 makeUA1RatioDiv(TGraphAsymmErrors* gSTAR,
35  double TAA, double TAAErr=0);
36 
37 // error on TAA
38 TH1D*
39 makeUA1TAAHistError(TGraphAsymmErrors* gSTAR,
40  double TAA, double TAAError);
41 
42 // scaleLow is the lower limit of the scale
43 TGraphAsymmErrors*
44 makeUA1ScaleError(TGraphAsymmErrors* gSTAR,
45  float scale, float scaleLow,float scaleHigh);
46 
47 // manuel's h-
48 TGraphAsymmErrors*
49 makeHMinus();
50 
51 // get rid of x errors
52 TGraphAsymmErrors*
53 removeXErrors(TGraphAsymmErrors* g);
54 
55 // kludge. divide by background
56 void
57 kludgeBackground(TGraphAsymmErrors* graph,float val)
58 {
59  double* xValues = graph->GetX();
60  double* yValues = graph->GetY();
61  double* errXLow = graph->GetEXlow();
62  double* errXHigh = graph->GetEXhigh();
63 
64  double* errYLow = graph->GetEYlow();
65  double* errYHigh = graph->GetEYhigh();
66 
67  for(int i=0; i<graph->GetN(); i++){
68  double y=yValues[i]*(1-val);
69  double yErr=errYLow[i]*(1-val);
70 
71  graph->SetPoint(i,xValues[i],y);
72  graph->SetPointError(i,errXLow[i],errXHigh[i],
73  yErr,yErr);
74 
75  }
76 }
77 
78 // kludge systematic errors
79 void
80 kludgeSystematics(TGraphAsymmErrors* graph)
81 {
82  // add 8% up to 4 GeV. 10% for the second to last bin
83  // 115% for the last bin
84 
85  double* xValues = graph->GetX();
86  double* yValues = graph->GetY();
87  double* errXLow = graph->GetEXlow();
88  double* errXHigh = graph->GetEXhigh();
89 
90  double* errYLow = graph->GetEYlow();
91  double* errYHigh = graph->GetEYhigh();
92 
93  int n=graph->GetN();
94  for(int i=0; i<n; i++){
95  double yErr;
96 
97  if(i==(n-2)) yErr=yValues[i]*.10;
98  else if(i==(n-1)) yErr=yValues[i]*.15;
99  else yErr=yValues[i]*.08;
100 
101  yErr += errYHigh[i];
102 
103  graph->SetPointError(i,errXLow[i],errXHigh[i],
104  yErr,yErr);
105 
106  }
107 
108 
109 }
110 
111 //--------------------------------------
112 /*
113  just scales the TGraph by some value.
114  would be nice if root had this as a method fot TGraph like TH1.
115  */
116 void
117 scale(TGraphAsymmErrors* graph,double scale)
118 {
119  if(DEBUG)printf("SCALE=%.2e\n (1/SCALE=%.2e)\n",scale,1./scale);
120  double* xValues = graph->GetX();
121  double* yValues = graph->GetY();
122  double* errXLow = graph->GetEXlow();
123  double* errXHigh = graph->GetEXhigh();
124 
125  double* errYLow = graph->GetEYlow();
126  double* errYHigh = graph->GetEYhigh();
127 
128  for(int i=0; i<graph->GetN(); i++){
129  double y=yValues[i], eYLow=errYLow[i], eYHigh=errYHigh[i];
130  double yScaled=yValues[i]*scale;
131  double eYLowScaled = errYLow[i]*scale;
132  double eYHighScaled =errYHigh[i]*scale;
133 
134  graph->SetPoint(i,xValues[i],yScaled);
135  graph->SetPointError(i,errXLow[i],errXHigh[i],
136  eYLowScaled,eYHighScaled);
137  if(DEBUG)
138  printf("\tbin=%d : %.2f<%.2f<%.2f: y=%.2e, yScaled=%.2e\n\t\t yErrLow=%.2e, yErrLowScaled=%.2e\n",
139  i,xValues[i]-errXLow[i],xValues[i],xValues[i]+errXHigh[i],
140  y,yScaled,eYLow,eYLowScaled);
141 
142  }
143 }
144 
145 /*
146  divide TGraphAsymmErrors by a float (and hopefully propagate errors).
147  only symmetric y errors.
148 */
149 void
150 divide(TGraphAsymmErrors* graph, double div, double divErr)
151 {
152  if(DEBUG)printf("DIVIDE=%.2e, err=%.2e\n",div,divErr);
153  double* xValues = graph->GetX();
154  double* yValues = graph->GetY();
155  double* errXLow = graph->GetEXlow();
156  double* errXHigh = graph->GetEXhigh();
157 
158  double* errYLow = graph->GetEYlow();
159  double* errYHigh = graph->GetEYhigh();
160 
161  for(int i=0; i<graph->GetN(); i++){
162  double y=yValues[i], eYLow=errYLow[i], eYHigh=errYHigh[i];
163  double yDiv=yValues[i]/div;
164  double eY = errYLow[i];
165  double eYDiv = errYLow[i]/div;
166 
167  if(divErr>0){
168  eYDiv=yDiv*sqrt((eY/y)*(eY/y) + (divErr/div)*(divErr/div));
169  }
170 
171 
172  graph->SetPoint(i,xValues[i],yDiv);
173  graph->SetPointError(i,errXLow[i],errXHigh[i],
174  eYDiv,eYDiv);
175  if(DEBUG)
176  printf("\tbin=%d : %.2f<%.2f<%.2f: y=%.2e, yScaled=%.2e\n\t\t yErrLow=%.2e, yErrLowScaled=%.2e\n",
177  i,xValues[i]-errXLow[i],xValues[i],xValues[i]+errXHigh[i],
178  y,yDiv,eY,eYDiv);
179 
180  }
181 
182 }
183 
184 
185 /*
186  scale for TGraphErrors
187 */
188 void
189 scale(TGraphErrors* graph,double scale)
190 {
191  if(DEBUG)printf("SCALE=%.2e\n (1/SCALE=%.2e)\n",scale,1./scale);
192  double* xValues = graph->GetX();
193  double* yValues = graph->GetY();
194  double* errX = graph->GetEX();
195  double* errY = graph->GetEY();
196 
197  for(int i=0; i<graph->GetN(); i++){
198  double y=yValues[i],eY=errY[i];
199  double yScaled=yValues[i]*scale;
200 
201  double eYScaled = errY[i]*scale;
202 
203  graph->SetPoint(i,xValues[i],yScaled);
204  graph->SetPointError(i,errX[i],eYScaled);
205 
206  /*
207  if(DEBUG)
208  printf("\tbin=%d : %.2f<%.2f<%.2f: y=%.2e, yScaled=%.2e\n\t\t yErrLow=%.2e, yErrLowScaled=%.2e\n",
209  i,xValues[i]-errXLow[i],xValues[i],xValues[i]+errXHigh[i],
210  y,yScaled,eYLow,eYLowScaled);
211  */
212  }
213 }
214 // ua1 data
215 
216 TGraphErrors* ua1Data200()
217 {
218  TGraphErrors* g=new TGraphErrors;
219  g->SetName("gUA1Data200GeV");
220 
221  const int n=40;
222  double ex[n]={0};
223  double x[n],y[n],ey[n];
224 
225 
226  x[0] = 0.25; y[0] = 62.0; ey[0] = 3.1;
227  x[1] = 0.35; y[1] = 36.2; ey[1] = 1.8;
228  x[2] = 0.45; y[2] = 20.8; ey[2] = 1.0;
229  x[3] = 0.55; y[3] = 12.32; ey[3] = 0.62;
230  x[4] = 0.65; y[4] = 7.33; ey[4] = 0.37;
231  x[5] = 0.75; y[5] = 4.54; ey[5] = 0.23;
232  x[6] = 0.85; y[6] = 2.77; ey[6] = 0.14;
233  x[7] = 0.95; y[7] = 1.724; ey[7] = 0.088;
234  x[8] = 1.05; y[8] = 1.162; ey[8] = 0.060;
235  x[9] = 1.15; y[9] = 0.769; ey[9] = 0.040;
236  x[10] = 1.25; y[10] = 0.521; ey[10] = 0.028;
237  x[11] = 1.35; y[11] = 0.344; ey[11] = 0.019;
238  x[12] = 1.45; y[12] = 0.240; ey[12] = 0.013;
239  x[13] = 1.55; y[13] = 0.1680; ey[13] = 0.0095;
240  x[14] = 1.65; y[14] = 0.1180; ey[14] = 0.0069;
241  x[15] = 1.75; y[15] = 0.0861; ey[15] = 0.0053;
242  x[16] = 1.85; y[16] = 0.0579; ey[16] = 0.0038;
243  x[17] = 1.95; y[17] = 0.0440; ey[17] = 0.0030;
244  x[18] = 2.05; y[18] = 0.0308; ey[18] = 0.0023;
245  x[19] = 2.15; y[19] = 0.0223; ey[19] = 0.0018;
246  x[20] = 2.25; y[20] = 0.0167; ey[20] = 0.0014;
247  x[21] = 2.35; y[21] = 0.0143; ey[21] = 0.0013;
248  x[22] = 2.45; y[22] = 0.0092; ey[22] = 0.0010;
249  x[23] = 2.55; y[23] = 0.00575; ey[23] = 0.00071;
250  x[24] = 2.65; y[24] = 0.00558; ey[24] = 0.00069;
251  x[25] = 2.75; y[25] = 0.00376; ey[25] = 0.00054;
252  x[26] = 2.85; y[26] = 0.00362; ey[26] = 0.00052;
253  x[27] = 2.95; y[27] = 0.00322; ey[27] = 0.00048;
254  x[28] = 3.05; y[28] = 0.00191; ey[28] = 0.00035;
255  x[29] = 3.15; y[29] = 0.00179; ey[29] = 0.00034;
256  x[30] = 3.25; y[30] = 0.00126; ey[30] = 0.00028;
257  x[31] = 3.35; y[31] = 0.00155; ey[31] = 0.00030;
258  x[32] = 3.45; y[32] = 0.00072; ey[32] = 0.00020;
259  x[33] = 3.55; y[33] = 0.00060; ey[33] = 0.00018;
260  x[34] = 3.65; y[34] = 0.000286; ey[34] = 0.000086;
261  x[35] = 3.85; y[35] = 0.000389; ey[35] = 0.000098;
262  x[36] = 4.10; y[36] = 0.000172; ey[36] = 0.000051;
263  x[37] = 4.50; y[37] = 0.000114; ey[37] = 0.000035;
264  x[38] = 5.70; y[38] = 0.0000135; ey[38] = 0.0000067;
265  x[39] = 6.70; y[39] = 0.0000066; ey[39] = 0.0000043;
266 
267 
268  for(int i=0; i<n; i++){
269  g->SetPoint(i,x[i],y[i]);
270  g->SetPointError(i,ex[i],ey[i]);
271  }
272  return g;
273 }
274 
275 // UA1(130)
276 
277 TF1* ua1Fit130(double scale)
278 {
279  double A = 266.6; double p0 = 1.895; double n = 12.98;
280  const char* fcn = "[3]*[0]*(1+x/[1])^-[2]";
281 
282  fUA1 = new TF1("fUA1)130",fcn,lowPt,highPt);
283  fUA1->SetParameter(0,A);
284  fUA1->SetParameter(1,p0);
285  fUA1->SetParameter(2,n);
286  fUA1->SetParameter(3,scale);
287 
288  return fUA1;
289 }
290 
291 /*
292  Note: this is the 130 GeV version from thomas.
293  couple functions to make ua1 spectra (fitted).
294  'graph' is your spectrum.(needed for binning)
295  no UA1 errors set.
296  2pi*A(1+pt/p0)^n .
297 
298  type 0: best guess
299  1: lower limit (of ua1)
300  2: upper limit
301 
302 */
303 
304 
305 TGraphAsymmErrors*
306 makeUA1(const TGraphAsymmErrors* graph, int type)
307 {
308  if(DEBUG)printf("MAKEUA1 type=%d (2piA(1=pt/p0)^-n)\n",type);
309  //
310  // get the values from graph
311  //
312  const int nBin = graph->GetN();
313  double* xValues = graph->GetX();
314  double* yValues = graph->GetY();
315  double* errXLow = graph->GetEXlow();
316  double* errXHigh = graph->GetEXhigh();
317 
318  TGraphAsymmErrors* gUA1=new TGraphAsymmErrors;
319  gUA1->SetName("ua1");
320 
321  // ual parameters from thomas's fit
322  // no errors set
323  double A,p0,n;
324 
325  switch(type){
326  case 0: // best guess
327  A = 266.6; p0 = 1.895; n = 12.98; break;
328  case 1: // lower limit
329  A = 260.6; p0 = 2.0701; n = 13.9; break;
330  case 2: // upper
331  A = 270.5; p0 = 1.8053; n = 12.514; break;
332  default:
333  cout << "wrong type=" << type << endl; exit(1);
334  }
335 
336  double lowPt=0, highPt=6;
337 
338  // A p0 n
339  const char* fcn = "2*pi*[0]*(1+x/[1])^-[2]";
340  TF1* fUA1;
341  fUA1 = new TF1("fUA1",fcn,lowPt,highPt);
342  fUA1->SetParameter(0,A);
343  fUA1->SetParameter(1,p0);
344  fUA1->SetParameter(2,n);
345 
346 
347  for(int i=0;i<nBin; i++){
348  //
349  // find x value
350  //
351  double lowEdge = xValues[i]-errXLow[i];
352  double highEdge = xValues[i]+errXHigh[i];
353  double binWidth = highEdge-lowEdge;
354 
355  double value = fUA1->Integral(lowEdge,highEdge)/binWidth;
356 
357  if(DEBUG)
358  printf("\tbin=%d : %.2f<%.2f<%.2f : val=%.4e\n",
359  i,lowEdge,xValues[i],highEdge,value);
360 
361  gUA1->SetPoint(i,xValues[i],value);
362  gUA1->SetPointError(i,errXLow[i],errXHigh[i],0,0);
363 
364  }
365  return gUA1;
366 }
367 
368 /*
369  graphA/graphB.
370  only sets errors based on the high errors.
371 */
372 
373 TGraphAsymmErrors*
374 divide(TGraphAsymmErrors* graphA,TGraphAsymmErrors* graphB)
375 {
376  if(DEBUG)printf("DIVIDE\n");
377 
378  const int nBin = graphA->GetN();
379  if(nBin != graphB->GetN()){
380  cout << "divide(..) :different bins sizes" << endl; exit(1);
381  }
382  // A
383  double* xValuesA = graphA->GetX();
384  double* yValuesA = graphA->GetY();
385  double* errXLowA = graphA->GetEXlow();
386  double* errXHighA = graphA->GetEXhigh();
387  double* errYA = graphA->GetEYhigh();
388 
389  // B
390  double* xValuesB = graphB->GetX();
391  double* yValuesB = graphB->GetY();
392  double* errYB = graphB->GetEYhigh();
393 
394  TGraphAsymmErrors* gRatio = new TGraphAsymmErrors;
395  gRatio->SetName("gRatio"); gRatio->SetTitle("gRatio");
396 
397  for(int i=0; i<nBin; i++){
398  if(yValuesB[i]==static_cast<double>(0)) {
399  cout << "divide: cant divide by 0" << endl;
400  }
401  double yRatio = yValuesA[i]/yValuesB[i];
402 
403  double errY =
404  sqrt(yRatio*yRatio*(
405  errYA[i]*errYA[i]/(yValuesA[i]*yValuesA[i])+
406  errYB[i]*errYB[i]/(yValuesB[i]*yValuesB[i])
407  )
408  );
409  gRatio->SetPoint(i,xValuesA[i],yRatio);
410  gRatio->SetPointError(i,errXLowA[i],errXHighA[i],
411  errY,errY);
412  double low=xValuesA[i]-errXLowA[i];
413  double high=xValuesA[i]+errXHighA[i];
414  if(DEBUG)
415  printf("\tbin=%d : %.2f<%.2f<%.2f : \n\t\tvalA=%.2e, valB=%.2e, A/B=%.2e\n\t\t errA=%.2e, errB=%.2e, errA/B=%.2e\n",
416  i,low,xValuesA[i],high,yValuesA[i],yValuesB[i],
417  yRatio,errYA[i],errYB[i],errY);
418 
419  }
420  return gRatio;
421 }
422 
423 
424 /*
425  STAR/UA1
426  'type' - see makeUA1 (0 is best guess, 1 is lower limit, 2 is upper limit)
427  'scale' is the overall scale of the ratio.
428  e.g. for central, we want scale = 1./TAA.
429  */
430 
431 TGraphAsymmErrors*
432 makeUA1Ratio(TGraphAsymmErrors* gSTAR,
433  int type, double scale)
434 {
435  TGraphAsymmErrors* gUA1=makeUA1(gSTAR,type);
436  TGraphAsymmErrors* gRatio=divide(gSTAR,gUA1);
437  double *y=gRatio->GetY();
438  scale(gRatio,scale);
439  return gRatio;
440 
441 }
442 /*
443  STAR/UA1 - but divide by TAA. TAA errors optional
444 */
445 TGraphAsymmErrors*
446 makeUA1RatioDiv(TGraphAsymmErrors* gSTAR,
447  double TAA, double TAAErr)
448 {
449  TGraphAsymmErrors* gUA1=makeUA1(gSTAR,0);
450  TGraphAsymmErrors* gRatio=divide(gSTAR,gUA1);
451  divide(gRatio,TAA,TAAErr);
452  return gRatio;
453 }
454 
455 
456 TH1D*
457 makeUA1TAAHistError(TGraphAsymmErrors* gSTAR,
458  double TAA, double TAAErr)
459 {
460  if(DEBUG)printf("makeUA1TAAHistError\n");
461 
462 
463  // make a ratio with the TAA errors included
464  TGraphAsymmErrors* gRatioErr=makeUA1RatioDiv(gSTAR,TAA,TAAErr);
465 
466  double* xValues = gRatioErr->GetX();
467  double* yValues = gRatioErr->GetY();
468  double* errXLow = gRatioErr->GetEXlow();
469  double* errXHigh = gRatioErr->GetEXhigh();
470 
471  double* errYLow = gRatioErr->GetEYlow();
472  double* errYHigh = gRatioErr->GetEYhigh();
473 
474  // find the bins of the hist
475  const int nAry=gRatioErr->GetN()+1;
476  const int nBin=nAry-1;
477  double* xAry=new double[nAry];
478  double* yAry=new double[nBin];// by bin (0 is the first bin)
479  double* yErr=new double[nBin];// by bin ...
480 
481  for(int i=0;i<nBin;i++){
482  float lowEdge=xValues[i]-errXLow[i];
483  xAry[i]=lowEdge;
484  yAry[i]=yValues[i];
485  yErr[i]=errYLow[i]; // assume symmetric
486  }
487  // the upper edge
488  xAry[nAry-1]=xValues[gRatioErr->GetN()-1]+errXHigh[gRatioErr->GetN()-1];
489 
490  // now make the 1d hist
491  TH1D* h=new TH1D("h","h",nBin,xAry);
492  for(int i=0;i<nBin;i++){
493  // cout << h->GetXaxis()->GetBinLowEdge(i+1)
494  // << "--" << h->GetXaxis()->GetBinUpEdge(i+1)<<endl;
495  h->SetBinContent(i+1,yAry[i]);
496  h->SetBinError(i+1,yErr[i]);
497  }
498  if(DEBUG) dump(h);
499  h->SetFillColor(17);
500 
501 
502  return h;
503 
504 }
505 
506 /*
507  makes graph with total errors.
508  scaleLow,scaleHigh is the error on the scale in absolute terms.
509  i.e. TAA=26+/-2. scale=1./26., scaleLow=1./28 (will give a lower
510  value to the ratio)
511  */
512 
513 TGraphAsymmErrors*
514 makeUA1ScaleError(TGraphAsymmErrors* gSTAR,
515  float scale, float scaleLow,float scaleHigh)
516 {
517  if(DEBUG)printf("ERROR\n");
518  if(DEBUG)printf("BEST\n");
519  TGraphAsymmErrors* gBestRatio=makeUA1Ratio(gSTAR,0,scale);
520 
521  // dont get confused. gLowRatio is the smaller ratio
522  // or larger denominator.(upper limit of UA1, lower scale)
523  //
524  if(DEBUG)printf("LOW RATIO\n");
525  TGraphAsymmErrors* gLowRatio=makeUA1Ratio(gSTAR,2,scaleLow); //upper ua1
526  if(DEBUG)printf("HIGH RATIO\n");
527  TGraphAsymmErrors* gHighRatio=makeUA1Ratio(gSTAR,1,scaleHigh);
528 
529  // clone the best ratio
530  // all the points are the same except for the errors.
531  // the high errors will be:
532  // gHighRatio-gBestRatio + gBestRatio error.
533  // similar for low errors.
534 
535  TGraphAsymmErrors* gError=
536  (TGraphAsymmErrors*)gBestRatio->Clone();
537  gError->SetName("gError");gError->SetTitle("gError");
538 
539  // ey best (should just be star's errors)
540  double* eYBest=gBestRatio->GetEYlow();
541 
542  // yvalues
543  double* yValuesBest=gBestRatio->GetY();
544  double* yValuesLow=gLowRatio->GetY();
545  double* yValuesHigh=gHighRatio->GetY();
546 
547  // x values
548  double* xValues=gBestRatio->GetX();
549  double* eXHigh=gBestRatio->GetEXhigh();
550  double* eXLow=gBestRatio->GetEXlow();
551 
552  if(DEBUG)printf("SET ERRORS\n");
553 
554  for(int i=0; i<gBestRatio->GetN(); i++){
555 
556  double eYHigh=yValuesHigh[i]-yValuesBest[i]; // should be positive
557  double eYLow=yValuesBest[i]-yValuesLow[i];
558  gError->SetPointError(i,0,0,eYLow+eYBest[i],eYHigh+eYBest[i]);
559 
560  if(DEBUG)
561  printf("\tbin=%d : %.2f<%.2f<%.2f : yBest=%.2e\n\t\teYhigh=%.2e, eYlow=%.2e, eYbest=%.2e\n",
562  i,xValues[i]-eXLow[i],xValues[i],xValues[i]+eXHigh[i],
563  yValuesBest[i],eYHigh,eYLow,eYBest[i]);
564 
565  }
566  return gError;
567 
568 }
569 
570 
571 // manuel's h-
572 // frank's h-+h+/2
573 // 0 is minbias
574 // 1 is 60-80 (add 70-80 and 60-70 / 2)
575 // 2 is central
576 
577 TGraphAsymmErrors*
578 makeHMinus(int cent)
579 {
580  cout << "###makeHMinus : " << cent << endl;
581 
582  char* fname="~/afs/real/PJacobs.root";
583 
584  TFile* frank = new TFile(fname);
585  if(!frank ||!frank->IsOpen()) exit(1);
586 
587  TGraphAsymmErrors* gHMinus = new TGraphAsymmErrors;
588  gHMinus->SetName("hMinus");
589 
590  char buf[100]; char* basename="hPlusMinus";
591  TH1* hMinus, *hMinusb;
592 
593  // change cases later.
594  switch(cent){
595  case 0: // minbias
596  sprintf(buf,"%s_00",basename);
597  hMinus=(TH1*)frank->Get(buf);
598  break;
599  case 1: // 60-80
600  sprintf(buf,"%s_02",basename);
601  hMinus=(TH1*)frank->Get(buf);
602  sprintf(buf,"%s_03",basename);
603  hMinusb=(TH1*)frank->Get(buf);
604  hMinus->Add(hMinusb); hMinus->Scale(0.5);
605  break;
606  case 2: // central
607  sprintf(buf,"%s_10",basename);
608  hMinus=(TH1*)frank->Get(buf);
609  break;
610  case 3: // 10-20
611  sprintf(buf,"%s_08",basename);
612  hMinus=(TH1*)frank->Get(buf);
613  break;
614  case 4: // 20-30
615  sprintf(buf,"%s_07",basename);
616  hMinus=(TH1*)frank->Get(buf);
617  break;
618  case 5: // 30-40
619  sprintf(buf,"%s_06",basename);
620  hMinus=(TH1*)frank->Get(buf);
621  break;
622  case 6: // 40-60
623  sprintf(buf,"%s_05",basename);
624  hMinus=(TH1*)frank->Get(buf);
625  sprintf(buf,"%s_04",basename);
626  hMinusb=(TH1*)frank->Get(buf);
627  hMinus->Add(hMinusb); hMinus->Scale(0.5);
628  break;
629  default:
630  cout << "Unknown cent = " << cent << endl; exit(1);
631  }
632  float lastpt=2.0;
633  int firstbin=2;
634 
635  // dump(hMinus);
636 
637  for(int i=firstbin; i<=hMinus->GetNbinsX(); i++){
638  double y=hMinus->GetBinContent(i);
639  double yerr=hMinus->GetBinError(i);
640  double x =hMinus->GetBinCenter(i);
641  double xerr=hMinus->GetBinWidth(i)/2.;
642 
643  if(x>lastpt) break;
644  gHMinus->SetPoint(i-firstbin,x,y);
645  gHMinus->SetPointError(i-firstbin,xerr,xerr,yerr,yerr);
646 
647  }
648  /*
649  gHMinus->SetPoint(0,.15,2226); gHMinus->SetPointError(0,.05,.05,94.84,94.84);
650  gHMinus->SetPoint(1,.25,1729.07);gHMinus->SetPointError(1,.05,.05,36.92,36.92);
651  gHMinus->SetPoint(2,.35,1117.08);gHMinus->SetPointError(2,.05,.05,26.09,26.09);
652  gHMinus->SetPoint(3,.45,732.16); gHMinus->SetPointError(3,.05,.05,13.7,13.7);
653  gHMinus->SetPoint(4,.55,478.81); gHMinus->SetPointError(4,.05,.05,7.30,7.30);
654  gHMinus->SetPoint(5,.65,304.32); gHMinus->SetPointError(5,.05,.05,5.99,5.99);
655  gHMinus->SetPoint(6,.75,204.69); gHMinus->SetPointError(6,.05,.05,3.97,3.97);
656  gHMinus->SetPoint(7,.85,139.48); gHMinus->SetPointError(7,.05,.05,4.62,4.62);
657  gHMinus->SetPoint(8,.95,97.62); gHMinus->SetPointError(8,.05,.05,2.04,2.04);
658  gHMinus->SetPoint(9,1.05,68.36); gHMinus->SetPointError(9,.05,.05,1.86,1.86);
659  gHMinus->SetPoint(10,1.15,46.96);gHMinus->SetPointError(10,.05,.05,1.48,1.48);
660  gHMinus->SetPoint(11,1.25,34.15);gHMinus->SetPointError(11,.05,.05,1.15,1.15);
661  gHMinus->SetPoint(12,1.35,23.82);gHMinus->SetPointError(12,.05,.05,0.98,0.98);
662  gHMinus->SetPoint(13,1.45,16.61);gHMinus->SetPointError(13,.05,.05,0.93,0.93);
663  gHMinus->SetPoint(14,1.55,11.72);gHMinus->SetPointError(14,.05,.05,0.74,0.74);
664  gHMinus->SetPoint(15,1.65,8.69); gHMinus->SetPointError(15,.05,.05,0.46,0.46);
665  gHMinus->SetPoint(16,1.75,6.20); gHMinus->SetPointError(16,.05,.05,0.43,0.43);
666  gHMinus->SetPoint(17,1.85,4.44); gHMinus->SetPointError(17,.05,.05,0.39,0.39);
667  gHMinus->SetPoint(18,1.95,3.38); gHMinus->SetPointError(18,.05,.05,0.32,0.32);
668  */
669  return gHMinus;
670 }
671 
672 TGraphAsymmErrors*
673 removeXErrors(TGraphAsymmErrors* g)
674 {
675  TGraphAsymmErrors* gClone=(TGraphAsymmErrors*)g->Clone();
676 
677  double* eYHigh=gClone->GetEYhigh();
678  double* eYLow=gClone->GetEYlow();
679 
680  for(int i=0;i<g->GetN(); i++){
681  gClone->SetPointError(i,0,0,eYLow[i],eYHigh[i]);
682  }
683  return gClone;
684 }
685 
686 void drawAxisBins(TGraphAsymmErrors* g,float size,float yMax)
687 {
688 
689  TLine* li=new TLine; li->SetLineWidth(1);
690 
691  double* x=g->GetX();
692  double* eXLow=g->GetEXlow();
693  double* eXHigh=g->GetEXhigh();
694 
695  // get the bin edges from the tgraph
696  for(int i=0; i<g->GetN(); i++){
697  double low = x[i]-eXLow[i];
698  double high= x[i]+eXHigh[i];
699 
700  // cout << low << "<" << x[i] << "<" << high << endl;
701  // cout << yMax << endl;
702  li->DrawLine(low,yMax-size,low,yMax);
703  li->DrawLine(high,yMax-size,high,yMax);
704  }
705 
706 }
Generated by   doxygen 1.8.5