StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
finish.C
1 /*
2  add all the corrections, backgrounds, etc, to the spectra
3 */
4 
5 #include <vector>
6 #include <utility>
7 #include "finish.h" // histograms defined here
8 #include "common/Name.cc"
9 #include "commonmacro/histutil.h"
10 
11 char* inrootname="~/wrk/real/links/P01hi.minbias.2000.hist/hispectra_cut2100.hist.root";
12 char* correctname="~/wrk/assoc/links/P01hi.central.HighpT_piminus_101.hist/momcorrection_cut2100.hist.root";
13 char* bkgname = "~/afs/real/background_test.hist.root";
14 char* outrootname="~/wrk/real/links/temp.hist/finish_test.root";
15 char* textdir = "TEXT";
16 
17 void finish(const char* inRootName=inrootname,
18  const char* correctionName=correctname,
19  const char* backgroundName=bkgname,
20  const char* outRootName=outrootname,
21  const char* textDir=textdir,
22  int cut = 2100,
23  int iter = 0)
24 {
25 
26  cout << "In file : " << inRootName << endl
27  << "mom correct : " << correctionName << endl
28  << "bkg : " << backgroundName << endl
29  << "out file : " << outRootName << endl
30  << "text dir : " << textDir << endl
31  << "cut : " << cut << endl
32  << "iter : " << iter << endl
33  << "offset type : " << offSetType << endl
34  << "fit high : " << fitHighPt << endl
35  << "cor high : " << corHighPt << endl;
36 
37  Int_t stat=0;
38  int doMomCorrection = (iter>0);
39 
40  gSystem->Load("StHiMicroAnalysis");
41  Cut::SetCut(cut); Cut::ShowCuts();
42 
43  char *trigger = "test";
44  //
45  // make output text files
46  //
47  initTextFiles(textDir,cut,iter);
48 
49  *ofVal << ">>>>>>>>>>" << endl;
50  *ofVal << "Doing momentum correction? " << doMomCorrection << endl;
51  cout << ">>>>>>>>>>>" << endl;
52  cout << "Doing momentum correction? " << doMomCorrection << endl;
53 
54  //
55  // read the momentum correction histograms
56  //
57  if(doMomCorrection) getCorrectionHistograms(correctionName);
58 
59 
60  //
61  // get the background stuff
62  //
63  getBackgroundHistograms(backgroundName);
64 
65 
66  //
67  // read the histograms from the root file
68  //
69  getHistograms(inRootName);
70 
71  //
72  // add the raw yields
73  //
74  doBothChargeSigns();
75 
76  //
77  // output root file
78  //
79  outRoot = new TFile(outRootName,"RECREATE");
80 
81  //
82  // compute all the corrections (momentum and background)
83  //
84  doCorrections();
85 
86 
87 
88  //
89  // two ways to get the offset of the pt bin.
90  // 1. find the weighted mean according to a fit of the histograms
91  // 2. use the histogram of the weighted mean.
92  // which is better?
93  //
94 
95  initFit();
96 
97  switch(offSetType){
98  case 1:
99  fitHistograms(); // fit to obtain the weighted mean
100  doOffSetGraphs(1); // from fit
101  break;
102  case 2:
103  doOffSetGraphs(2); // from histograms
104  break;
105  default:
106  cout << "wrong type"; return;
107  }
108 
109  //
110  // fit the graphs. these parameters will be used for mom correction
111  //
112  fitGraphs();
113 
114  //
115  // make spectra (just devide by 1/pt)
116  //
117  doSpectraGraphs();
118 
119  //
120  // make some plots
121  //
122  //draw(psDir,trigger,cut);
123 
124 
125 
126  if(1){
127  writeHistograms(outRootName);
128  }
129 
130 
131 
132 }
133 //____________________
134 
135 void initTextFiles(const char* textDir, int cut, int iter)
136 {
137 
138  sprintf(name,"%s/corValues_cut%d_iter%d.txt",textDir,cut,iter);
139  cout << "correction values : " << name << endl;
140  ofCor = new ofstream(name);
141  if(!ofCor) { cout << "cannot create " << name << endl;}
142 
143  // sprintf(name,"%s/ptBinValues_cut%d_iter%d.txt",textDir,cut,iter);
144  //ofPtBin = new ofstream(name);
145 
146  sprintf(name,"%s/values_cut%d_iter%d.txt",textDir,cut,iter);
147  cout << "values : " << name << endl;
148  ofVal = new ofstream(name);
149 
150  //sprintf(name,"%s/ua1_%s_cut%d.txt",textDir,trigger,cut);
151  // ofRatio = new ofstream(name);
152 }
153 
154 
155 //____________________
156 
157 void getCorrectionHistograms(const char* correctionName)
158 {
159  cout << "\t**getCorrectionHistograms**" << endl;
160 
161  Int_t stat = 0;;
162  //
163  // momentum correction
164  //
165  TFile* corRoot = new TFile(correctionName);
166 
167  if(!corRoot) return 1;
168 
169  for(Int_t iBin=0; iBin<nVarBin; iBin++){
170  // deltra crippled
171  setName(name,"h1PtDCrRatio",iBin);
172  correction[iBin].crip.h1PtDCrRatio = (TH1D*) corRoot->Get(name);
173  if(!correction[iBin].crip.h1PtDCrRatio) {
174  cout << "no " << name << endl; exit(1);
175  }
176  showHistogramValues(correction[iBin].crip.h1PtDCrRatio,
177  "MOMENTUM CORRECTIONS");
178 
179  // weight used (reality check)
180  setName(name,"fCorrected",iBin);
181  correction[iBin].crip.fCorrected = (TF1*)corRoot->Get(name);
182  if(!correction[iBin].crip.fCorrected){
183  cout << "no " << name << endl; exit(1);
184  }
185  for(Int_t iCharge=0; iCharge<2; iCharge++){
186 
187  setName(name,"h1PtDCrRatio",iBin,sPM[iCharge].Data());
188  correction[iBin].cripPM[iCharge].h1PtDCrRatio = (TH1D*)corRoot->Get(name);
189  showHistogramValues(correction[iBin].cripPM[iCharge].h1PtDCrRatio,
190  "MOMENTUM CORRECTIONS");
191 
192  if(!correction[iBin].cripPM[iCharge].h1PtDCrRatio) {
193  cout << "no " << name << endl; exit(1);
194  }
195  setName(name,"fCorrected",iBin,sPM[iCharge].Data());
196  correction[iBin].cripPM[iCharge].fCorrected = (TF1*)corRoot->Get(name);
197  if(!correction[iBin].cripPM[iCharge].fCorrected){
198  cout << "no " << name << endl; exit(1);
199  stat++;
200  }
201  }
202 
203 
204  } // varbin
205 
206  //
207  // reality check
208  //
209 
210  cout << "Momentum correction used weights ..." << endl;
211 
212  for(Int_t iBin=0; iBin<nVarBin; iBin++){
213  cout << "Bin " << iBin << endl;
214  cout << "\tBoth charges"
215  << "\tpt0="
216  <<correction[iBin].crip.fCorrected->GetParameter(0)
217  << "\t\tptPow="
218  <<correction[iBin].crip.fCorrected->GetParameter(1)
219  <<endl;
220 
221  for(Int_t iCharge=0; iCharge<2; iCharge++){
222  cout << "\t"<< sPM[iCharge].Data()
223  << "\tpt0="
224  <<correction[iBin].cripPM[iCharge].fCorrected->GetParameter(0)
225  << "\t\tptPow="
226  <<correction[iBin].cripPM[iCharge].fCorrected->GetParameter(1)
227  <<endl;
228 
229  }
230  }
231 
232 }
233 //________________
234 //
235 // need to fix this
236 //
237 
238 void getBackgroundHistograms(const char* backgroundName)
239 {
240  cout << "\t**getBackgroundHistograms()**" << endl;
241 
242  TFile* bkgRoot = new TFile(backgroundName);
243  Int_t stat = 0;
244 
245  // for now, set the +&- backgrounds to be the same
246 
247  strcpy(name,"h1BackGround");
248 
249  char* dcaType = "DcaGl";
250 
251  for(int iBin=0; iBin<nVarBin; iBin++){
252  setName(name,"h1Background",iBin);
253  sprintf(name,"%s%s",name,dcaType);
254  correction[iBin].bkg.h1Background=(TH1D*) bkgRoot->Get(name);
255  //background[iBin].h1Background = (TH1D*) bkgRoot->Get(name);
256  if(!correction[iBin].bkg.h1Background){
257  cout << "Cannot find " << name << " in "
258  << backgroundName << endl; exit(1);
259  }
260  showHistogramValues(correction[iBin].bkg.h1Background,"");
261 
262  for(int ic=0; ic<2; ic++){
263  setName(name,"h1Background",iBin,sPM[ic].Data());
264  sprintf(name,"%s%s",name,dcaType);
265  correction[iBin].bkgPM[ic].h1Background=(TH1D*) bkgRoot->Get(name);
266 
267  if(!correction[iBin].bkg.h1Background){
268  cout << "Cannot find " << name << " in "
269  << backgroundName << endl; exit(1);
270  }
271  showHistogramValues(correction[iBin].bkgPM[ic].h1Background,"");
272  }
273  }
274 
275 }
276 //____________________
277 
278 void getHistograms(const char* inRootName)
279 {
280  cout <<"\t**getHistograms()**" << endl;
281 
282 
283  TFile* inRoot = new TFile(inRootName);
284  if(!inRoot || !inRoot->IsOpen()){
285  cout << "cannot open " << inRootName << endl; exit(1);
286  }
287 
288  // vertex z
289  //
290  h1VertexZCut = (TH1D*)inRoot->Get("h1VertexZCut");
291 
292  // number of events and the eta cuts
293  //
294  h1NEvent = (TH1D*) inRoot->Get("h1NEvent");
295  h1EtaCut = (TH1D*) inRoot->Get("h1EtaCut");
296 
297 
298  nEvent = h1NEvent->GetBinContent(1);
299  etaCut[0] = h1EtaCut->GetBinContent(1);
300  etaCut[1] = h1EtaCut->GetBinContent(2);
301  scale = (1./nEvent)*(1./(etaCut[1]-etaCut[0]));
302 
303  cout << "n event : " << nEvent << endl;
304  cout << "eta cut : " << etaCut[0] << ", " << etaCut[1] << endl;
305  cout << "scale : " << scale << endl;
306 
307  *ofVal << "n event : " << nEvent << endl
308  << "eta cut : " << etaCut[0] << ", " << etaCut[1] << endl
309  << "scale : " << scale << endl;
310 
311  for(Int_t iBin=0; iBin<nVarBin; iBin++){
312 
313  setName(name,"h1OneOverPt",iBin);
314  varBin[iBin].mean.h1OneOverPt = (TH1D*)inRoot->Get(name);
315 
316  setName(name,"h1WeightedMean",iBin);
317  varBin[iBin].mean.h1WeightedMean = (TH1D*)inRoot->Get(name);
318 
319  // raw
320  //
321  setName(name,"h1Raw",iBin);
322  varBin[iBin].spec.h1Raw = (TH1D*) inRoot->Get(name);
323 
324  showHistogramValues(varBin[iBin].spec.h1Raw,"YIELDS");
325 
326  //
327  // efficiency corrected
328  //
329  setName(name,"h1EffCorrected",iBin);
330  varBin[iBin].spec.h1EffCorrected = (TH1D*) inRoot->Get(name);
331  showHistogramValues(varBin[iBin].spec.h1EffCorrected,"YIELDS");
332 
333  //
334  // corrected (background + momentum resolution)
335  //
336  setName(name,"h1Corrected",iBin);
337  varBin[iBin].spec.h1Corrected = (TH1D*) inRoot->Get(name);
338  showHistogramValues(varBin[iBin].spec.h1Corrected,"YIELDS");
339 
340 
341 
342 
343  for(Int_t iCharge=0; iCharge<2; iCharge++){
344 
345  // mean stuff
346  //
347  setName(name,"h1OneOverPt",iBin,sPM[iCharge].Data());
348  varBin[iBin].meanPM[iCharge].h1OneOverPt = (TH1D*)inRoot->Get(name);
349 
350  setName(name,"h1WeightedMean",iBin,sPM[iCharge].Data());
351  varBin[iBin].meanPM[iCharge].h1WeightedMean = (TH1D*)inRoot->Get(name);
352 
353  //
354  // raw
355  //
356  setName(name,"h1Raw",iBin,sPM[iCharge].Data());
357  varBin[iBin].specPM[iCharge].h1Raw =
358  (TH1D*) inRoot->Get(name);
359  showHistogramValues(varBin[iBin].specPM[iCharge].h1Raw,"YIELDS");
360  //
361  // efficiency corrected
362  //
363  setName(name,"h1EffCorrected",iBin,sPM[iCharge].Data());
364  varBin[iBin].specPM[iCharge].h1EffCorrected =
365  (TH1D*) inRoot->Get(name);
366  showHistogramValues(varBin[iBin].specPM[iCharge].h1EffCorrected,"YIELDS");
367 
368  //
369  // corrected (background + momentum resolution)
370  //
371  setName(name,"h1Corrected",iBin,sPM[iCharge].Data());
372  varBin[iBin].specPM[iCharge].h1Corrected =
373  (TH1D*) inRoot->Get(name);
374  showHistogramValues(varBin[iBin].specPM[iCharge].h1Corrected,"YIELDS");
375  }
376 
377  //***** ratios
378 
379  setName(name,"h1RawRatio",iBin);
380  varBin[iBin].h1RawRatio = (TH1D*) inRoot->Get(name);
381 
382  setName(name,"h1CorrectedRatio",iBin);
383  varBin[iBin].h1CorrectedRatio = (TH1D*) inRoot->Get(name);
384 
385  } // varBin
386  // cout << "done " <<endl;
387 }
388 //__________________
389 
390 void doBothChargeSigns()
391 {
392  cout << "\t**doBothChargsSigns()**" << endl;
393 
394  //
395  // just scale by a half
396  //
397 
398  for(Int_t iBin=0; iBin<nVarBin; iBin++){
399  varBin[iBin].spec.h1Raw->Scale(0.5);
400  //showHistogramValues(varBin[iBin].spec.h1Raw,"SCALED 1/2");
401 
402  varBin[iBin].spec.h1EffCorrected->Scale(0.5);
403  //showHistogramValues(varBin[iBin].spec.h1EffCorrected,"SCALED 1/2");
404 
405  varBin[iBin].spec.h1Corrected->Scale(0.5);
406  //showHistogramValues(varBin[iBin].spec.h1Corrected,"SCALED 1/2");
407  }
408 
409 }
410 
411 //____________________
412 
413 void doCorrections()
414 {
415  cout << "\t**doCorrections()**" << endl;
416 
417  // bkg and momentum resolution
418 
419  for(Int_t iBin=0; iBin<nVarBin; iBin++){
420 
421  computeCorrection(varBin[iBin].spec.h1EffCorrected,
422  correction[iBin].crip.h1PtDCrRatio,
423  correction[iBin].bkg.h1Background,
424  varBin[iBin].spec.h1Corrected);
425 
426  showHistogramValues(varBin[iBin].spec.h1Corrected,"CORRECTED");
427 
428  for(Int_t iCharge=0; iCharge<2; iCharge++){
429 
430  computeCorrection(varBin[iBin].specPM[iCharge].h1EffCorrected,
431  correction[iBin].cripPM[iCharge].h1PtDCrRatio,
432  correction[iBin].bkgPM[iCharge].h1Background,
433  varBin[iBin].specPM[iCharge].h1Corrected);
434  showHistogramValues(varBin[iBin].specPM[iCharge].h1Corrected,
435  "CORRECTED");
436  }
437 
438 
439  }
440 
441  // scale the eff corrected separately for safety
442  //
443  for(Int_t iBin=0; iBin<nVarBin; iBin++){
444  varBin[iBin].spec.h1Raw->Scale(scale);
445  varBin[iBin].spec.h1Corrected->Scale(scale);
446  varBin[iBin].spec.h1EffCorrected->Scale(scale);
447  for(Int_t iCharge=0; iCharge<2; iCharge++){
448  varBin[iBin].specPM[iCharge].h1Raw->Scale(scale);
449  varBin[iBin].specPM[iCharge].h1Corrected->Scale(scale);
450  varBin[iBin].specPM[iCharge].h1EffCorrected->Scale(scale);
451  }
452 
453  }
454 
455 
456  // compute the pt bin separately
457  //
458  for(Int_t iBin=0; iBin<nVarBin; iBin++){
459  computePtBin(varBin[iBin].spec.h1Raw);
460  computePtBin(varBin[iBin].spec.h1EffCorrected);
461  computePtBin(varBin[iBin].spec.h1Corrected);
462 
463  for(Int_t iCharge=0; iCharge<2; iCharge++){
464  computePtBin(varBin[iBin].specPM[iCharge].h1Raw);
465  computePtBin(varBin[iBin].specPM[iCharge].h1EffCorrected);
466  computePtBin(varBin[iBin].specPM[iCharge].h1Corrected);
467 
468  }
469  }
470 }
471 
472 //____________________
473 
474 void initFit()
475 {
476  cout << "\tinitFit()" << endl;
477 
478 
479  fDist = new TF1("fDist",dist,fitLowPt,fitHighPt);
480  fMean = new TF1("fMean",mean,fitLowPt,fitHighPt);
481 
482  Double_t p2 = 1e4;
483 
484  fDist->SetParameters(2.0,13.,p2);
485  fMean->SetParameters(2.0,13.,p2);
486  const int nParam = 3;
487 
488  Double_t param[nParam];
489  fDist->GetParameters(param);
490 
491  //cout << "initial param : " ;
492  //for(int i=0; i<nParam; i++){
493  // cout << "\t" << param[i] ;
494  // }
495  cout << endl;
496 }
497 
498 //____________________
499 //
500 // fit the spectra. use the weighted mean of each bin
501 // for the 1/pt scale
502 
503 void fitHistograms()
504 {
505  cout << "\t**fit histograms**" << endl;
506 
507  //
508  // dummy canvas
509  //
510  TCanvas* dummy = new TCanvas("dummy","dummy",100,100,500,500);
511  gPad->SetLogy();
512 
513 
514  for(Int_t iBin=0; iBin<nVarBin; iBin++){
515  varBin[iBin].spec.fRawTmp =
516  fit(varBin[iBin].spec.h1Raw);
517 
518  varBin[iBin].spec.fEffCorrectedTmp =
519  fit(varBin[iBin].spec.h1EffCorrected);
520 
521  varBin[iBin].spec.fCorrectedTmp =
522  fit(varBin[iBin].spec.h1Corrected);
523 
524  for(Int_t iCharge=0; iCharge<2; iCharge++){
525  varBin[iBin].specPM[iCharge].fRawTmp =
526  fit(varBin[iBin].specPM[iCharge].h1Raw);
527 
528  varBin[iBin].specPM[iCharge].fEffCorrectedTmp =
529  fit(varBin[iBin].specPM[iCharge].h1EffCorrected);
530 
531  varBin[iBin].specPM[iCharge].fCorrectedTmp =
532  fit(varBin[iBin].specPM[iCharge].h1Corrected);
533  }
534  }
535 
536  delete dummy;
537 }
538 
539 
540 //____________________
541 //
542 // only used for getting the weighted mean from the fit
543 //
544 
545 TF1* fit(TH1* h1)
546 {
547  cout << "\tfitting " << h1->GetName() << endl;
548 
549  double param[3];
550  TF1* f=0;
551 
552  // first get a reasonable initial parameters
553  TH1* hClone=(TH1*)h1->Clone();
554  hClone->Fit("fDist","QR");
555  f = hClone->GetFunction("fDist");
556 
557  f->GetParameters(param);
558  fDist->SetParameters(param);
559 
560  // fit for real
561  hClone->Fit("fDist","QR");
562  f = hClone->GetFunction("fDist");
563  /*
564  cout << "\t\tparams : "
565  << "p0=" << f->GetParameter(0)
566  << ",p1=" << f->GetParameter(1)
567  << ",p2=" << f->GetParameter(2)
568  << endl;
569  */
570  TString sName = h1->GetName();
571 
572  Ssiz_t last = findLast(sName,".");
573 
574  sName.Replace(last,2,"fTmp");
575 
576  f->SetName(sName.Data());
577  f->SetTitle(sName.Data());
578 
579  return f;
580 }
581 
582 //_____________________
583 //
584 // fit tgraphs to obtain fit parameters for mom correction
585 //
586 
587 TF1* fit(TGraphAsymmErrors* g)
588 {
589  cout << "\tfitting " << g->GetName() << endl;
590 
591  double param[3];
592  TF1* f=0;
593 
594  // first get a reasonable initial parameters
595  TGraphAsymmErrors* gClone=(TGraphAsymmErrors*)g->Clone();
596  gClone->Fit("fDist","QR");
597  f = gClone->GetFunction("fDist");
598 
599  f->GetParameters(param);
600  fDist->SetParameters(param);
601 
602  // fit for real
603  gClone->Fit("fDist","QR");
604  f = gClone->GetFunction("fDist");
605 
606  cout << "\t\tparams : "
607  << "p0=" << f->GetParameter(0)
608  << ",p1=" << f->GetParameter(1)
609  << ",p2=" << f->GetParameter(2)
610  << endl;
611 
612  TString sName = g->GetName();
613 
614  Ssiz_t last = findLast(sName,".");
615 
616  sName.Replace(last,1,"f");
617  f->SetName(sName.Data());
618  f->SetTitle(sName.Data());
619 
620  cout << "\t\t fcn name " << f->GetName() << endl;
621 
622  return f;
623 }
624 
625 //_____________________
626 
627 
628 Double_t getWeightedMean(TF1* f,Double_t lowEdge, Double_t upEdge)
629 {
630  Double_t param[3];
631 
632  f->GetParameters(param);
633 
634  //
635  // set the parameters for the mean function
636  //
637  fMean->SetParameters(param); // defined in the header file
638 
639  return fMean->Integral(lowEdge,upEdge)/f->Integral(lowEdge,upEdge);
640 
641 }
642 
643 //_________________
644 //
645 // make the off set tgraphs
646 
647 
648 TGraphAsymmErrors*
649 graphOffset(TH1D* h1,TF1* f,TH1* hMean,Int_t type)
650 {
651 
652  Int_t lowBin = h1->GetXaxis()->FindBin(fitLowPt);
653  Int_t highBin = h1->GetXaxis()->FindBin(fitHighPt-.0001);
654  Int_t nBin = highBin-lowBin+1;
655 
656  TArrayD yErrorAry(nBin), yValueAry(nBin), xValueAry(nBin),
657  xLowErrorAry(nBin), xHighErrorAry(nBin);
658 
659  // dont fill the first and last bin.
660  // histograms bin numbered from 1 to n.
661  // array number from 0 to n-1.
662 
663  Int_t binAry(0);
664 
665  for(Int_t iBin=lowBin; iBin<=highBin; iBin++,binAry++){
666 
667  //
668  // x error
669  //
670  Double_t lowEdge = h1->GetXaxis()->GetBinLowEdge(iBin);
671  Double_t upEdge = h1->GetXaxis()->GetBinUpEdge(iBin);
672  Double_t mean = 0;
673 
674  switch(type){
675  case 1:
676  mean = getWeightedMean(f,lowEdge,upEdge); break;
677  case 2:
678  mean = h1Mean->GetBinContent(h1Mean->GetXaxis()->FindBin((lowEdge+upEdge)/2)); break;
679  default:
680  mean = getWeightedMean(f,lowEdge,upEdge);
681  }
682 
683  xLowErrorAry.AddAt(fabs(mean-lowEdge),binAry);
684  xHighErrorAry.AddAt(fabs(mean-upEdge),binAry);
685 
686  xValueAry.AddAt(mean,binAry); // x
687 
688  yErrorAry.AddAt(h1->GetBinError(iBin),binAry); // y error
689 
690  yValueAry.AddAt(h1->GetBinContent(iBin),binAry); // y
691 
692  }
693 
694  //
695  // now actually make the graph
696  //
697 
698  TGraphAsymmErrors* g = new TGraphAsymmErrors(nBin,
699  xValueAry.GetArray(),
700  yValueAry.GetArray(),
701  xLowErrorAry.GetArray(),
702  xHighErrorAry.GetArray(),
703  yErrorAry.GetArray(),
704  yErrorAry.GetArray());
705 
706  TString sName = h1->GetName();
707  Ssiz_t last = findLast(sName,".");
708 
709  sName.Replace(last,2,"g");
710 
711  g->SetName(sName.Data());
712  g->SetTitle(sName.Data());
713 
714  return g;
715 }
716 //__________________
717 
718 TGraphAsymmErrors*
719 graphSpectra(TGraphAsymmErrors* g)
720 {
721 
722  const Int_t nBin = g->GetN();
723 
724  Double_t* x = g->GetX();
725  Double_t* y = g->GetY();
726  Double_t* eXlow = g->GetEXlow();
727  Double_t* eXhigh = g->GetEXhigh();
728  Double_t* eYlow = g->GetEYlow();
729 
730  TArrayD yErrorAry(nBin), yValueAry(nBin);
731 
732  // divide each bin by the value of that bin.
733  // also divide the errors.
734  //
735 
736  for(Int_t i=0; i<nBin; i++){
737  Double_t mean = x[i];
738 
739  yErrorAry.AddAt(eYlow[i]/mean,i);
740  yValueAry.AddAt(y[i]/mean,i);
741  }
742 
743 
744  TGraphAsymmErrors* gSpec
745  = new TGraphAsymmErrors(nBin,x,yValueAry.GetArray(),eXlow,eXhigh,
746  yErrorAry.GetArray(),yErrorAry.GetArray());
747 
748  TString sName = g->GetName();
749  Ssiz_t last = findLast(sName,".");
750 
751  sName.Replace(last,1,"gSpec");
752 
753  gSpec->SetName(sName.Data());
754  gSpec->SetTitle(sName.Data());
755 
756  return gSpec;
757 
758 
759 }
760 
761 //____________________
762 //
763 //
764 
765 void doOffSetGraphs(Int_t type)
766 {
767 
768  for(Int_t iBin=0; iBin<nVarBin; iBin++){
769 
770  varBin[iBin].spec.gRaw =
771  graphOffset(varBin[iBin].spec.h1Raw,
772  varBin[iBin].spec.fRawTmp,
773  varBin[iBin].mean.h1WeightedMean,type);
774 
775  varBin[iBin].spec.gEffCorrected =
776  graphOffset(varBin[iBin].spec.h1EffCorrected,
777  varBin[iBin].spec.fEffCorrectedTmp,
778  varBin[iBin].mean.h1WeightedMean,type);
779 
780  varBin[iBin].spec.gCorrected =
781  graphOffset(varBin[iBin].spec.h1Corrected,
782  varBin[iBin].spec.fCorrectedTmp,
783  varBin[iBin].mean.h1WeightedMean,type);
784 
785 
786  for(Int_t iCharge=0; iCharge<2; iCharge++){
787 
788  varBin[iBin].specPM[iCharge].gRaw =
789  graphOffset(varBin[iBin].specPM[iCharge].h1Raw,
790  varBin[iBin].specPM[iCharge].fRawTmp,
791  varBin[iBin].mean.h1WeightedMean,type);
792 
793 
794  varBin[iBin].specPM[iCharge].gEffCorrected =
795  graphOffset(varBin[iBin].specPM[iCharge].h1EffCorrected,
796  varBin[iBin].specPM[iCharge].fEffCorrectedTmp,
797  varBin[iBin].meanPM[iCharge].h1WeightedMean,type);
798 
799  varBin[iBin].specPM[iCharge].gCorrected =
800  graphOffset(varBin[iBin].specPM[iCharge].h1Corrected,
801  varBin[iBin].specPM[iCharge].fCorrectedTmp,
802  varBin[iBin].meanPM[iCharge].h1WeightedMean,type);
803 
804  }
805 
806  }
807 }
808 
809 //____________________
810 //
811 // fit the offset graphs.
812 // these parameters are used for the momentum corrections
813 //
814 
815 void fitGraphs()
816 {
817  cout << "\t**fitGraphs()**" << endl;
818 
819  //
820  // dummy canvas
821  //
822  TCanvas* dummy = new TCanvas("dummy","dummy",100,100,500,500);
823  gPad->SetLogy();
824 
825  for(Int_t iBin=0; iBin<nVarBin; iBin++){
826 
827  varBin[iBin].spec.fEffCorrected =
828  fit(varBin[iBin].spec.gEffCorrected);
829 
830  varBin[iBin].spec.fCorrected =
831  fit(varBin[iBin].spec.gCorrected);
832 
833  for(Int_t iCharge=0; iCharge<2; iCharge++){
834 
835  varBin[iBin].specPM[iCharge].fEffCorrected =
836  fit(varBin[iBin].specPM[iCharge].gEffCorrected);
837 
838  varBin[iBin].specPM[iCharge].fCorrected =
839  fit(varBin[iBin].specPM[iCharge].gCorrected);
840  }
841  }
842 
843  delete dummy;
844 
845 }
846 //____________________
847 
848 void doSpectraGraphs()
849 {
850  cout << "\t**doSpectraGraphs()**" << endl;
851 
852  for(Int_t iBin=0; iBin<nVarBin; iBin++){
853  varBin[iBin].spec.gSpecCorrected =
854  graphSpectra(varBin[iBin].spec.gCorrected);
855 
856  for(Int_t iCharge=0; iCharge<2; iCharge++){
857 
858  varBin[iBin].specPM[iCharge].gSpecCorrected =
859  graphSpectra(varBin[iBin].specPM[iCharge].gCorrected);
860  }
861 
862  }
863 
864 }
865 //____________________
866 
867 
868 
869 //____________________
870 
871 void draw(const char* psDir, const char* trigger, const char* cut)
872 {
873 
874 }
875 
876 
877 //____________________
878 //
879 // set the background error at 100%
880 //
881 
882 void computeCorrection(TH1D* hRc, TH1D* hDCrOverRc,
883  TH1D* hBkg, TH1D* hCorrected)
884 {
885 
886  Stat_t rc, dCr, rcError, dCrError, error, corrected, bkgCorrected,
887  bkgCorrectedError, center, bkg, bkgError;
888 
889  Int_t bin;
890  Int_t nBin = hRc->GetNbinsX();
891  Int_t maxBkgBin = hBkg->GetNbinsX();
892 
893  *ofCor << "******************************************************" << endl;
894  *ofCor << "## computeCorrection ##" << endl;
895  *ofCor << "## corrected=" << hCorrected->GetName() << endl
896  << " rc=" << hRc->GetName() << endl
897  << " bkg=" << hBkg->GetName() << endl;
898 
899  int firstCorBin=hRc->GetXaxis()->FindBin(corLowPt);
900  *ofCor << "First correction bin is " << firstCorBin << endl;
901  for(Int_t iBin=1; iBin<=nBin; iBin++){
902 
903  rc = hRc->GetBinContent(iBin); // real data
904 
905  // dCr weighted w/flat input
906  dCr = (hDCrOverRc) ? hDCrOverRc->GetBinContent(iBin) :0;
907 
908  // background
909  center = hRc->GetXaxis()->GetBinCenter(iBin);
910  bin = hBkg->GetXaxis()->FindBin(center);
911  if(bin==0 || bin>maxBkgBin){
912  bkg = 0;
913  }
914  else{
915  bkg = hBkg->GetBinContent(bin);
916  }
917 
918  // changed the definition a little
919 
920  bkgCorrected = rc*(1-bkg);
921 
922  rcError = hRc->GetBinError(iBin);
923 
924  // set the background error as the same as the bkg
925  bkgError = bkg;
926 
927  bkgCorrectedError
928  = TMath::Sqrt(rcError*rcError*(1+bkg*bkg) + rc*rc*bkgError*bkgError);
929 
930  // now do the dCr stuff
931  //
932  if(iBin>=firstCorBin){
933  corrected = bkgCorrected*(1-dCr);
934  }
935  else{
936  *ofCor << "\tskipping momentum correction for this bin" << endl;
937  corrected = bkgCorrected;
938  }
939  // forget about the crippled error
940  //
941  error = bkgCorrectedError;
942 
943  hCorrected->SetBinContent(iBin,corrected);
944  hCorrected->SetBinError(iBin,error);
945 
946  float fracError = (rc>0) ? ((error/rc)*100) : 0;
947  float fracRcError = (rc>0) ? ((rcError/rc)*100) : 0;
948 
949 
950  *ofCor << " rc : " << hRc->GetXaxis()->GetBinLowEdge(iBin) << "-"
951  << hRc->GetXaxis()->GetBinUpEdge(iBin)
952  << "\trc bin content : " << rc
953  << ", crippled : " << dCr
954  << ", background : " << bkg << endl
955  << "\tbkg corrected : " << bkgCorrected << endl
956  << "\tcorrected : " << corrected << endl;
957  *ofCor << "\t rc error : " << rcError
958  << "(" << fracRcError << " %)"
959  << ", dCrError : " << dCrError
960  << ", bkg error: " << bkgError
961  << endl
962  << "\tbkg cor error : " << bkgCorrectedError << endl
963  << "\ttotal error : " << error << endl
964  << "\ttotal err % : " << fracError << endl;
965 
966  }
967 }
968 //____________________
969 // just divide histograms by bin center
970 
971 
972 void ptDivide(TH1* h1)
973 {
974  Int_t nBin=h1->GetNbinsX();
975  for(int iBin=1; iBin<=nBin; iBin++){
976  double content= h1->GetBinContent(iBin);
977  double center = h1->GetBinCenter(iBin);
978  double error = h1->GetBinError(iBin);
979  if(content) h1->SetBinContent(iBin,content/center);
980  if(error)h1->SetBinError(iBin,error/center);
981  }
982 
983 }
984 
985 
986 void computePtBin(TH1* h1)
987 {
988 
989  double binWidth;
990  double binContent;
991  double binError;
992 
993  Int_t nBin = h1->GetNbinsX();
994  for(Int_t iBin=1; iBin<=nBin; iBin++){
995  binWidth = h1->GetBinWidth(iBin);
996  binContent= h1->GetBinContent(iBin);
997  if(binContent>0) binContent /= binWidth;
998  binError = h1->GetBinError(iBin);
999  if(binError!=0) binError /= binWidth;
1000 
1001  h1->SetBinContent(iBin,binContent);
1002  h1->SetBinError(iBin,binError);
1003 
1004  }
1005 }
1006 //_______________
1007 //
1008 // also saves the TF1's and TGraphs
1009 
1010 
1011 void writeHistograms(const char* outRootName)
1012 {
1013  cout << "writeHistograms() : " << outRootName << endl;
1014 
1015  // TFile* outRoot = new TFile(outRootName,"RECREATE");
1016 
1017  for(Int_t iBin=0; iBin<nVarBin; iBin++){
1018 
1019  varBin[iBin].spec.h1Raw->Write();
1020  varBin[iBin].spec.h1EffCorrected->Write();
1021  varBin[iBin].spec.h1Corrected->Write();
1022  varBin[iBin].spec.gRaw->Write();
1023  varBin[iBin].spec.gEffCorrected->Write();
1024  varBin[iBin].spec.gCorrected->Write();
1025  varBin[iBin].spec.gSpecCorrected->Write();
1026  varBin[iBin].spec.fCorrected->Write();
1027 
1028  showHistogramValues(varBin[iBin].spec.h1Raw,"FINAL");
1029  showHistogramValues(varBin[iBin].spec.h1EffCorrected,"FINAL");
1030  showHistogramValues(varBin[iBin].spec.h1Corrected,"FINAL");
1031 
1032  showTGraphValues(varBin[iBin].spec.gRaw,"FINAL");
1033  showTGraphValues(varBin[iBin].spec.gEffCorrected,"FINAL");
1034  showTGraphValues(varBin[iBin].spec.gCorrected,"FINAL");
1035  showTGraphValues(varBin[iBin].spec.gSpecCorrected,"FINAL");
1036 
1037  for(Int_t iCharge=0; iCharge<2; iCharge++){
1038 
1039  // ua1 stuff
1040  //ua1[iBin].ratioPM[iCharge].gRatioWeight->Write();
1041  //ua1[iBin].ratioPM[iCharge].gRatio->Write();
1042 
1043  varBin[iBin].specPM[iCharge].h1Raw->Write();
1044  varBin[iBin].specPM[iCharge].h1EffCorrected->Write();
1045  varBin[iBin].specPM[iCharge].h1Corrected->Write();
1046  varBin[iBin].specPM[iCharge].gEffCorrected->Write();
1047  varBin[iBin].specPM[iCharge].gCorrected->Write();
1048  varBin[iBin].specPM[iCharge].gSpecCorrected->Write();
1049  varBin[iBin].specPM[iCharge].fCorrected->Write();
1050 
1051  showHistogramValues(varBin[iBin].specPM[iCharge].h1Raw,"FINAL");
1052  showHistogramValues(varBin[iBin].specPM[iCharge].h1EffCorrected,"FINAL");
1053  showHistogramValues(varBin[iBin].specPM[iCharge].h1Corrected,"FINAL");
1054 
1055  showTGraphValues(varBin[iBin].specPM[iCharge].gRaw,"FINAL");
1056  showTGraphValues(varBin[iBin].specPM[iCharge].gEffCorrected,"FINAL");
1057  showTGraphValues(varBin[iBin].specPM[iCharge].gCorrected,"FINAL");
1058  showTGraphValues(varBin[iBin].specPM[iCharge].gSpecCorrected,"FINAL");
1059 
1060  }
1061 
1062  }
1063 
1064 }
1065 
1066 
1067 
1068 
1069 //________________________
1070 
1071 Ssiz_t findLast(TString temp, const char* a)
1072 {
1073  Ssiz_t last = 0;
1074 
1075  while(1){
1076  Ssiz_t pos = temp.First(a);
1077  if(pos<0) break;
1078  last += ++pos;
1079  temp.Remove(0,pos);
1080  }
1081 
1082  return last;
1083 
1084 }
Generated by   doxygen 1.8.5