AliPhysics  2797316 (2797316)
CalculateAveragePt.C
Go to the documentation of this file.
1 // MACRO TO COMPUTE THE <pt> from data for D mesons
2 // works for D0, Ds and D+; D* case should be tested since the side bands definition is less clear
3 //
4 // Requires
5 // a 2D (pt,mass) histogram, with fine pt bins in pt (finer than those used for the spectrum, e.g. 100 or 200 MeV/c)
6 // a histogram of background yield vs pt (with the pt bins in which the <pt> has to be computed)
7 // a histogram of the invariant mass sigma (with the pt bins in which the <pt> has to be computed)
8 // optionally a histogram of signal yield (not fully needed) (with the pt bins in which the <pt> has to be computed)
9 // optionally a histogram of the mean of the signal peak obtained from the fit (with the pt bins in which the <pt> has to be computed)
10 // The last 4 histograms are among the output of the standard macro used for fitting the mass distributions
11 // optionally (but strongly recommended) an histograms with prompt D reco efficiencies (in finer bins than those in which the <pt> is wanted, otherwise they will have no effect).
12 // The binning of the efficiency histo should match that of the pt axis of the 2D (pt,mass) plots x rebin (where rebin is used to rebin the data histogram with the method Rebin)
13 // e.g. if you have data in 200 MeV/c and efficiency in 500 MeV/c bins than the minimum binning for which data and efficiencis match is 1 GeV/c
14 // -> you have to use rebin=5 (or multiples of 5) to have the macro working (it stops otherwise)
15 // Options:
16 // rebin (see above)
17 // usefitforsubtraction: with this switched on the pt shape of the backgrond is fitted with an exponential and the fit function is used for the background subtraction in the signal region
18 // (useful if there are stat fluctuations in the background pt shape distribution, but take care and check the fit quality in the control plots)
19 // correct for efficiency: if switched on, the pt distribution of the signal (i.e. that in the signal peak after back subtraction) is weighted by the efficiency
20 // it accounts for the pt dependence of the efficiency -> compulsory if a realistic (corrected) <pt> has to be computed
21 //
22 // useParGenAccOverLimAcc: uses a parametrization for correcting the efficiencies by GenAcc/LimAcc. It is useful since the GenAcc/LimAcc factor might
23 // not be available in very fine bins, especially at high pt
24 //
25 //
26 // There are standard methods hard-coded for D0 (DoStandardForD0), D+ and Ds as an example (of course you have to change the paths of the files and the name of the histos!!).
27 // To run them (e.g. for D0):
28 // -1 (my advice): copy the macro in your working directory and change the name of the paths and histo there (it helps to retrieve the results in the future).
29 // Or copy the ingredient files in the working directory,
30 // 0- set the path of the files
31 // 1- .L CalculateAveragePt.C++ (compile the macro, not really needed)
32 // 2- DoStandardForD0(5,kFALSE,kTRUE,kTRUE,0,3)
33 // ... it should be easy :)
34
35 //
36 //
37 // For problems write to: andrea.rossi@cern.ch
38 //
40
41 #include <TH2F.h>
42 #include <TH1F.h>
43 #include <TH2D.h>
44 #include <TH1D.h>
45 #include <TCanvas.h>
46 #include <TF1.h>
47 #include <TFile.h>
48 #include <TDirectory.h>
49 #include <TMath.h>
50 #include <TString.h>
51 #include <TGraphAsymmErrors.h>
52 #include <TAxis.h>
53 #include <TLegend.h>
54
55 TH2F *hPtInvMass=0x0;
56 Int_t nptbins=0;
57 Double_t *ptbinlimits=0x0;
58 Double_t *rawsignal=0x0,*rawback=0x0,*sigma=0x0,*meansignal=0x0;
59 Double_t nsigmaSignal=3.,nsigmaSBstart=5.;
60 Double_t mesonMass=1.8645;
61 Int_t nbinsx;//=hPtInvMass->GetNbinsX();
62 Int_t nbinsy;//=hPtInvMass->GetNbinsY();
63 Double_t binwidthpt;//=hPtInvMass->GetYaxis()->GetBinWidth(1);// ASSUME ALL BINS ARE WITH THE SAME WIDTH
64 Double_t binwidthInvMss;//=hPtInvMass->GetXaxis()->GetBinWidth(1);// ASSUME ALL BINS ARE WITH THE SAME WIDTH
65 Double_t ptmin;//=hPtInvMass->GetYaxis()->GetBinLowEdge(1);
66 Double_t ptmax;//=hPtInvMass->GetYaxis()->GetBinUpEdge(nbinsy);
67 Bool_t useFitForSubtraction=kFALSE;
68 Bool_t useParGenAccOverLimAcc=kFALSE;
69 Bool_t subtractSB=kTRUE;
70 Bool_t corrForEff=kFALSE;
71 void SetSubtractSB(Bool_t subtract){subtractSB=subtract;}
72 void SetCorrForEff(Bool_t correff){corrForEff=correff;}
73 void SetUseFitForSubtraction(Bool_t useFit){useFitForSubtraction=useFit;}
74 void SetNsigmaForSignal(Double_t nsigm){nsigmaSignal=nsigm;}
75 void SetNsigmaStartSB(Double_t nsigm){nsigmaSBstart=nsigm;
76 }
77 void CalculateAveragePt(Int_t rebin=1,Int_t firstbin=0,Int_t lastbin=0);
78 TF1* ParametricGenAccOverLimAccCorr(){// parametric form of GenAcc/LimAcc correction for D0->Kpi
79  TF1 *fGenAccOverLimAcc=new TF1("fat7","[2]+([0]*TMath::ATan([1]*x+[3])+[4]*TMath::Log(1.+x/[5]))",0,24);
80  fGenAccOverLimAcc->SetParameters(2.76153e-01,6.69658e-01,8.45865e-01,-1.87709e+00,2.89845e+03,1.39651e+05);
81  return fGenAccOverLimAcc;
82 }
83
84 Int_t standrebin[8]={2,2,2,2,2,2,4,5};
85 TCanvas **cPtDistrNoSubtr;
86 TF1 **fitfunc;
88 void SetHistosEfficiency(TH1D *hNum,TH1D *hDenum){
89  hEfficNum=(TH1D*)hNum;
90  hEfficDenum=(TH1D*)hDenum;
91
92 }
93 TGraphAsymmErrors *grAvRawYieldSpectrum;
94 TGraphAsymmErrors *grAvPtVSPtmean;
95
96 TGraphAsymmErrors *grBackAvPtVSPtmean;
97
98 void SetPtBinLimits(const Int_t npt,Double_t *ptbinlim){
99  if(nptbins>0){
100  delete ptbinlimits;ptbinlimits=0x0;
101  }
102  nptbins=npt;
103  ptbinlimits=new Double_t[nptbins+1];
104  for(Int_t bin=0;bin<=nptbins;bin++){
105  ptbinlimits[bin]=ptbinlim[bin];
106  }
107  cPtDistrNoSubtr=new TCanvas*[nptbins];
108  fitfunc=new TF1*[nptbins];
109 }
110
111 void SetPtBinLimits(TH1 *histo){
112
113  if(nptbins>0){
114  delete ptbinlimits;ptbinlimits=0x0;
115  }
116  TAxis *xax=histo->GetXaxis();
117  nptbins=xax->GetNbins();
118  ptbinlimits=new Double_t[nptbins+1];
119  for(Int_t bin=0;bin<=nptbins;bin++){
120  ptbinlimits[bin]=xax->GetBinLowEdge(bin+1);
121  }
122  cPtDistrNoSubtr=new TCanvas*[nptbins];
123  fitfunc=new TF1*[nptbins];
124 }
125
126 void SetHistRawSignal(TH1D *hS){hSignal=hS;
127  rawsignal=&(hSignal->GetArray())[1];
128  return;}
129
130 void SetHistMean(TH1D *hM){hMeanSignal=hM;
131  meansignal=&(hMeanSignal->GetArray())[1];
132  return;
133 }
134
135 void SetHistRawBack(TH1D *hB){hBack=hB;
136  rawback=&(hBack->GetArray())[1];
137  hAvRawYieldSpectrum=new TH1D(*hBack);
138  hAvRawYieldSpectrum->SetName("hAvRawYieldDistrTotPeak");
139  hAvRawYieldSpectrum->SetTitle("Average Raw Yield in bins after subtraction from bin counting");
140  hAvRawYieldSpectrum->Reset(0);
141  hAvRawYieldSpectrum->SetLineColor(kRed);
142  return;
143 }
144 void SetHistSigma(TH1D *hSig){hSigma=hSig;
145  sigma=&(hSigma->GetArray())[1];
146  return;
147 }
148
149 void SetHistoMassPt(TH2F *h2){
150  hPtInvMass=new TH2F(*h2);
151  nbinsx=hPtInvMass->GetNbinsX();
152  nbinsy=hPtInvMass->GetNbinsY();
153  binwidthpt=hPtInvMass->GetYaxis()->GetBinWidth(1);// ASSUME ALL BINS ARE WITH THE SAME WIDTH
154  binwidthInvMss=hPtInvMass->GetXaxis()->GetBinWidth(1);// ASSUME ALL BINS ARE WITH THE SAME WIDTH
155  ptmin=hPtInvMass->GetYaxis()->GetBinLowEdge(1);
156  ptmax=hPtInvMass->GetYaxis()->GetBinUpEdge(nbinsy);
157
158 }
159
160 TH1D* CheckBinningAndMerge(TH1D *hA,TH1D *hB,Double_t precision=0.001,Double_t minX=-9999.,Double_t maxX=-9999.){// Checks that histo hB binning is finer and fits hA binning
161  // and return an Histo with hB bins integrated to match than hA bins
162  // TO BE ADDED A CHECK ON THE MIN AND MAX VALUES OF THE AXES IN THE 2 PLOTS
163  // NO PARTICULAR TREATMENT OF THE ERROR IN THE INTEGRATION (SHOULD BE OK BUT FOR EFFIENCIES)
164
165  TH1D *hWork=new TH1D(*hA);
166  hWork->SetName(Form("%sNewBins",hB->GetName()));
167  hWork->SetTitle(hB->GetTitle());
168
169  Double_t lebinA,uebinA;
170  Int_t binBleA=0,binBueA=0;
171  Int_t binBleALast=0,binBueALast=0;
172  Int_t maxbinA=hA->GetNbinsX();
173  Int_t minbinA=1;
174  if(maxX!=-9999.){
175  maxbinA=hA->FindBin(maxX);
176  }
177  if(minX!=-9999.){
178  minbinA=hA->FindBin(minX);
179  }
180
181  for(Int_t binA=minbinA;binA<=maxbinA;binA++){
182  uebinA=hA->GetXaxis()->GetBinUpEdge(binA);
183  lebinA=hA->GetXaxis()->GetBinLowEdge(binA);
184  // Check that Low Edge of binA is close to the low or upper edge of a bin in hB
185  // within precition
186  binBleA=hB->FindBin(lebinA);
187
188  if(TMath::Abs(hB->GetXaxis()->GetBinLowEdge(binBleA)-lebinA)>precision&&TMath::Abs(hB->GetXaxis()->GetBinUpEdge(binBleA)-lebinA)>precision){
189  printf("The bins of the efficiency histo do not match those of the pt spectrum histo (edges 1):\n Cannot correct for efficiensies \n");
190  delete hWork;
191  return 0x0;
192  }
193  if(TMath::Abs(hB->GetXaxis()->GetBinLowEdge(binBleA)-lebinA)>TMath::Abs(hB->GetXaxis()->GetBinUpEdge(binBleA)-lebinA))binBleA++;
194  if(binBleA<=binBleALast){
195  printf("Problems with rebinning Hist B, non compatibility with different bin of Hist A (1)\n");
196  delete hWork;
197  return 0x0;
198  }
199  if(hB->GetBinWidth(binBleA)-hA->GetBinWidth(binA)>precision){
200  printf("The bins of the efficiency histo do not match those of the pt spectrum histo (bin width 1):\n Cannot correct for efficiensies \n");
201  delete hWork;
202  return 0x0;
203  }
204  // Check that Upper Edge of binA is close to the low or upper edge of a bin in hB
205  // within precition
206  binBueA=hB->FindBin(uebinA);
207  if(TMath::Abs(hB->GetXaxis()->GetBinLowEdge(binBueA)-uebinA)>precision&&TMath::Abs(hB->GetXaxis()->GetBinUpEdge(binBueA)-uebinA)>precision){
208  printf("The bins of the efficiency histo do not match those of the pt spectrum histo (edges 2) (bin %d Vs %d):\n Cannot correct for efficiensies \n",binA,binBueA);
209  printf("%f Vs %f or %f \n",uebinA,hB->GetXaxis()->GetBinLowEdge(binBueA),hB->GetXaxis()->GetBinUpEdge(binBueA));
210  delete hWork;
211  return 0x0;
212  }
213  if(TMath::Abs(hB->GetXaxis()->GetBinLowEdge(binBueA)-uebinA)<TMath::Abs(hB->GetXaxis()->GetBinUpEdge(binBueA)-uebinA))binBueA--;
214  if(binBueA<=binBueALast){
215  printf("Problems with rebinning Hist B, non compatibility with different bin of Hist A (2) \n");
216  delete hWork;
217  return 0x0;
218  }
219  if(hB->GetBinWidth(binBueA)-hA->GetBinWidth(binA)>precision){
220  printf("The bins of the efficiency histo do not match those of the pt spectrum histo (bin width 2):\n Cannot correct for efficiensies \n");
221  delete hWork;
222  return 0x0;
223  }
224  if(binBleA>binBueA){
225  printf("Mathcing failure (probably a bug in the code \n");
226  delete hWork;
227  return 0x0;
228  }
229
230  hWork->SetBinContent(binA,hB->Integral(binBleA,binBueA));
231  }
232
233  return hWork;
234
235 }
236
237 Bool_t CorrectForEfficiency(TH1D *hPtHisto){
238  TH1D *hWorkNum=0x0,*hWorkDenum=0x0;
239  // hWorkNum=CheckBinningAndMerge(hPtHisto,hEfficNum,0.01,1.,23.9);
240  // printf("ptmin and pt max: %f, %f \n",ptbinlimits[0],ptbinlimits[nptbins]);
241  hWorkNum=CheckBinningAndMerge(hPtHisto,hEfficNum,0.01,ptbinlimits[0]*1.0001,ptbinlimits[nptbins]*0.9999);
242  hWorkNum->Sumw2();
243
244  if(!hWorkNum)return kFALSE;
245  hWorkDenum=CheckBinningAndMerge(hPtHisto,hEfficDenum,0.01,ptbinlimits[0]*1.0001,ptbinlimits[nptbins]*0.9999);
246  // hWorkDenum=CheckBinningAndMerge(hPtHisto,hEfficDenum,0.01,1.,23.9);;
247  hWorkDenum->Sumw2();
248  if(!hWorkDenum){
249  delete hWorkNum;
250  return kFALSE;
251  }
252  TH1D *hWorkNumCl=(TH1D*)hWorkDenum->Clone("hWorkNumCl");
253  hWorkNumCl->Divide(hWorkNum,hWorkDenum,1,1,"B");
255  TF1 *fGenAccoOverLimAcc=ParametricGenAccOverLimAccCorr();
256  for(Int_t j=1;j<=hWorkNumCl->GetNbinsX();j++){
257  printf(" EFFF for %f: %f * %f \n",hWorkNumCl->GetBinCenter(j), hWorkNumCl->GetBinContent(j),fGenAccoOverLimAcc->Eval(hWorkNumCl->GetBinCenter(j)));
258  hWorkNumCl->SetBinContent(j,hWorkNumCl->GetBinContent(j)*fGenAccoOverLimAcc->Eval(hWorkNumCl->GetBinCenter(j)));
259  hWorkNumCl->SetBinError(j,hWorkNumCl->GetBinError(j)*fGenAccoOverLimAcc->Eval(hWorkNumCl->GetBinCenter(j)));
260
261  }
262  }
263  hPtHisto->Divide(hWorkNumCl);
264
265  delete hWorkNum;
266  delete hWorkNumCl;
267  delete hWorkDenum;
268  return kTRUE;
269 }
270
271 TH1D *HistoPtShapeFromData(Int_t ptbin,Int_t rebin=1.){
272  if(!hPtInvMass){
273  printf("No histo set \n");
274  return 0x0;
275  }
276  Double_t mean=mesonMass;
277  if(meansignal!=0x0){
278  mean=meansignal[ptbin];
279  printf("USING MEAN FROM MASS FIT: %f\n",mean);
280  }
281  TString namehist="hPtDistrSB_Bin",funcname;
282  namehist+=ptbin;
283
284  TString cname="cPtDistrNoSubBin";
285  cname+=ptbin;
286  cPtDistrNoSubtr[ptbin]=new TCanvas(cname.Data(),"CPtDistrNoSub",700,700);
287
288  printf("After Setting Canva \n");
289  Int_t binleftpt=-1,binrightpt=-1;
290  TH1D *hPtDistrSB=new TH1D(namehist.Data(),"Side band Pt distribution",nbinsy,ptmin,ptmax);
291  namehist="hPtDistrTotPeak_Bin";
292  namehist+=ptbin;
293  printf("After Setting SB hist \n");
294  TH1D *hPtDistrTotPeak=new TH1D(namehist.Data(),"Total Pt distribution in signal mass reagion",nbinsy,ptmin,ptmax);
295  for(Int_t jbi=1;jbi<hPtDistrTotPeak->GetNbinsX();jbi++){
296  hPtDistrTotPeak->SetBinContent(jbi,0);
297  hPtDistrTotPeak->SetBinError(jbi,0);
298  hPtDistrSB->SetBinContent(jbi,0);
299  hPtDistrSB->SetBinError(jbi,0);
300  }
301
302
303
304  printf("Before loop on xbins \n");
305  for(Int_t j=1;j<nbinsx;j++){
306  if((hPtInvMass->GetXaxis()->GetBinLowEdge(j)<=mean-nsigmaSBstart*sigma[ptbin])||(hPtInvMass->GetXaxis()->GetBinUpEdge(j)>=mean+nsigmaSBstart*sigma[ptbin])){
307  for(Int_t k=1;k<nbinsy;k++){
308  if(hPtInvMass->GetYaxis()->GetBinLowEdge(k)>=ptbinlimits[ptbin]&&hPtInvMass->GetYaxis()->GetBinUpEdge(k)<=ptbinlimits[ptbin+1]){
309  hPtDistrSB->SetBinContent(k,hPtDistrSB->GetBinContent(k)+hPtInvMass->GetBinContent(j,k));
310  hPtDistrSB->SetBinError(k,TMath::Sqrt(hPtDistrSB->GetBinError(k)*hPtDistrSB->GetBinError(k)+hPtInvMass->GetBinError(j,k)*hPtInvMass->GetBinError(j,k)));
311  }
312  }
313  }
314  else if((hPtInvMass->GetXaxis()->GetBinLowEdge(j)>=mean-nsigmaSignal*sigma[ptbin])&&(hPtInvMass->GetXaxis()->GetBinUpEdge(j)<=mean+nsigmaSignal*sigma[ptbin])){
315  for(Int_t k=1;k<nbinsy;k++){
316  if(hPtInvMass->GetYaxis()->GetBinLowEdge(k)>=ptbinlimits[ptbin]&&hPtInvMass->GetYaxis()->GetBinUpEdge(k)<=ptbinlimits[ptbin+1]){
317  if(binleftpt==-1)binleftpt=k;
318  hPtDistrTotPeak->SetBinContent(k,hPtDistrTotPeak->GetBinContent(k)+hPtInvMass->GetBinContent(j,k));
319  hPtDistrTotPeak->SetBinError(k,TMath::Sqrt(hPtDistrTotPeak->GetBinError(k)*hPtDistrTotPeak->GetBinError(k)+hPtInvMass->GetBinError(j,k)*hPtInvMass->GetBinError(j,k)));
320  binrightpt=k;
321  }
322  }
323  }
324  }
325  hPtDistrSB->Sumw2();
326  hPtDistrTotPeak->Sumw2();
327  if(rebin!=1){//Rebin and recalculate needed info
328  hPtDistrSB->Rebin(rebin);
329  hPtDistrTotPeak->Rebin(rebin);
330  binleftpt=-1;
331  binleftpt=-1;
332  for(Int_t k=1;k<=hPtDistrTotPeak->GetNbinsX();k++){
333  if(hPtDistrTotPeak->GetBinLowEdge(k)>=ptbinlimits[ptbin]&&hPtDistrTotPeak->GetBinLowEdge(k+1)<=ptbinlimits[ptbin+1]){
334  if(binleftpt==-1)binleftpt=k;
335  binrightpt=k;
336  }
337  }
338  binwidthpt=hPtDistrTotPeak->GetBinWidth(1);
339  }
340
341  Double_t avptbin=0.,errA=0.,errB=0.,errC=0.,errisq=0.;
342  // CALCULATE BACKGROUND AVERAGE PT
343  for(Int_t bbin=1;bbin<=hPtDistrSB->GetNbinsX();bbin++){
344  avptbin+=hPtDistrSB->GetBinCenter(bbin)*hPtDistrSB->GetBinContent(bbin)*hPtDistrSB->GetBinWidth(bbin);
345  errisq=hPtDistrSB->GetBinError(bbin)*hPtDistrSB->GetBinError(bbin);
346  errA+=hPtDistrSB->GetBinCenter(bbin)*hPtDistrSB->GetBinCenter(bbin)*errisq;
347  errB+=hPtDistrSB->GetBinCenter(bbin)*errisq;
348  errC+=errisq;
349  // printf("ERROR VALUES: %f, %f \n",errA,errB);
350  }
351  avptbin/=hPtDistrSB->Integral("width");
352  errisq=TMath::Sqrt(errA-2.*avptbin*errB+avptbin*avptbin*errC)/hPtDistrSB->Integral("width");
353  if(!grBackAvPtVSPtmean){
354  grBackAvPtVSPtmean=new TGraphAsymmErrors();
355  grBackAvPtVSPtmean->SetName("grBackAvPtVSPtmean");
356  }
357  Int_t grbin=grBackAvPtVSPtmean->GetN();
358  grBackAvPtVSPtmean->SetPoint(grbin,0.5*(ptbinlimits[ptbin+1]+ptbinlimits[ptbin]),avptbin);
359  grBackAvPtVSPtmean->SetPointError(grbin,0.5*(ptbinlimits[ptbin+1]-ptbinlimits[ptbin]),0.5*(ptbinlimits[ptbin+1]-ptbinlimits[ptbin]),errisq,errisq);
360
361  TH1D *hPtSignalFromSubtraction=new TH1D(*hPtDistrTotPeak);
362  namehist="hPtSignalFromSubtraction_Bin";
363  namehist+=ptbin;
364  hPtSignalFromSubtraction->SetName(namehist.Data());
365  hPtSignalFromSubtraction->SetTitle("Signal Pt distr from subtraction");
366
367
368  if(subtractSB){
369  if(!useFitForSubtraction)hPtSignalFromSubtraction->Add(hPtDistrSB,-rawback[ptbin]*(nsigmaSignal/3.)/hPtDistrSB->Integral());// ROUGH LINEAR SCALING OF BACKGOUND!!
370  else{
371  funcname="expofit";
372  funcname+=ptbin;
373  fitfunc[ptbin]=new TF1(funcname.Data(),"expo",ptbinlimits[ptbin],ptbinlimits[ptbin+1]);
374  TCanvas *cTempt=new TCanvas();
375  cTempt->cd();
376  hPtDistrSB->Fit(fitfunc[ptbin],"RLEV","",ptbinlimits[ptbin],ptbinlimits[ptbin+1]);
377  delete cTempt;
378  for(Int_t jb=binleftpt;jb<=binrightpt;jb++){
379  printf("Err before sub: %f \n",hPtSignalFromSubtraction->GetBinError(jb));
380  hPtSignalFromSubtraction->SetBinContent(jb,hPtSignalFromSubtraction->GetBinContent(jb)-fitfunc[ptbin]->Eval(hPtSignalFromSubtraction->GetBinCenter(jb))/fitfunc[ptbin]->Integral(ptbinlimits[ptbin],ptbinlimits[ptbin+1])*binwidthpt*rawback[ptbin]*(nsigmaSignal/3.));// ROUGH LINEAR SCALING OF BACKGOUND!!
381  printf("Err after sub: %f \n\n",hPtSignalFromSubtraction->GetBinError(jb));
382  }
383  }
384  }
385
386
387  cPtDistrNoSubtr[ptbin]->Divide(1,2);
388  cPtDistrNoSubtr[ptbin]->cd(1);
389  hPtDistrSB->Draw();
390  cPtDistrNoSubtr[ptbin]->cd(2);
391  if(corrForEff){
392  if(!CorrectForEfficiency(hPtSignalFromSubtraction)){
393  printf("Correction for efficieny failed \n");
394  }
395  }
396  hPtDistrTotPeak->Draw();
397  hPtSignalFromSubtraction->Draw("sames");
398
399
400  return hPtSignalFromSubtraction;
401
402 }
403
404 void CalculateAveragePt(TH2* hMassPt,TH1D *hB,TH1D *hSigm,TH1D *hEffNum=0x0,TH1D *hEffDenum=0x0,TH1D *hS=0x0,TH1D *hMean=0x0,Int_t rebin=1,Int_t firstbin=3,Int_t lastbin=8){
405
406  SetPtBinLimits(hB);
407  SetHistoMassPt((TH2F*)hMassPt);
408  SetHistRawBack(hB);
409  SetHistSigma(hSigm);
410  if(((!hEffNum)&&(hEffDenum))||((hEffNum)&&(!hEffDenum))){
411  printf("Two histos are needed for the efficiency: numerator and denumerator (for rebinning) \n");
412  }
413  else if(hEffNum){
414  SetHistosEfficiency(hEffNum,hEffDenum);
415  // SetHistEffNum(hEffNum);
416  // SetHistEffDenum(hEffDenum);
417  }
418  if(hS)SetHistRawSignal(hS);
419  if(hMean)SetHistMean(hMean);
420  printf("Histos set \n");
421  CalculateAveragePt(rebin,firstbin,lastbin);
422 }
423
424 void CalculateAveragePt(Int_t rebin,Int_t firstbin,Int_t lastbin){
425
426  TH1D *hPtSpectra;
427  if(rebin>0){
428  hPtSpectra=new TH1D("hPtDistrTotPeak","Total Pt distribution in signal mass reagion",nbinsy,ptmin,ptmax);
429  if(rebin!=1)hPtSpectra->Rebin(rebin);
430  }
431  else if(rebin==-1){
432  printf("Standard rebinning not implemented yet \n");return;
433  }
434  else {
435  printf("Wrong rebin numnber \n");return;
436  }
437
438  grAvRawYieldSpectrum=new TGraphAsymmErrors();
439  grAvPtVSPtmean=new TGraphAsymmErrors();
440  grBackAvPtVSPtmean=new TGraphAsymmErrors();
441  grBackAvPtVSPtmean->SetName("grBackAvPtVSPtmean");
442
443
444  Double_t avptbin=0.,errA,errB,errC,errisq;
445  for(Int_t bin=firstbin;bin<=lastbin;bin++){
446  // printf("Before Getting Subtracted Histo \n");
447  avptbin=0.;
448  errA=0.;
449  errB=0.;
450  errC=0.;
451  errisq=0.;
452
453  TH1D *h=HistoPtShapeFromData(bin,rebin);
455  hAvRawYieldSpectrum->SetBinContent(bin+1,h->Integral("width")/hAvRawYieldSpectrum->GetBinWidth(bin+1));
456  for(Int_t bbin=1;bbin<=h->GetNbinsX();bbin++){
457  avptbin+=h->GetBinCenter(bbin)*h->GetBinContent(bbin)*h->GetBinWidth(bbin);
458  errisq=h->GetBinError(bbin)*h->GetBinError(bbin);
459  errA+=h->GetBinCenter(bbin)*h->GetBinCenter(bbin)*errisq;
460  errB+=h->GetBinCenter(bbin)*errisq;
461  errC+=errisq;
462  // printf("ERROR VALUES: %f, %f \n",errA,errB);
463  }
464  avptbin/=h->Integral("width");
465  errisq=TMath::Sqrt(errA-2.*avptbin*errB+avptbin*avptbin*errC)/h->Integral("width");
466
467  grAvRawYieldSpectrum->SetPoint(bin-firstbin,avptbin,h->Integral("width")/hAvRawYieldSpectrum->GetBinWidth(bin+1));
468  grAvRawYieldSpectrum->SetPointError(bin-firstbin,avptbin-ptbinlimits[bin],ptbinlimits[bin+1]-avptbin,errisq/2.,errisq/2.);
469
470  grAvPtVSPtmean->SetPoint(bin-firstbin,0.5*(ptbinlimits[bin+1]+ptbinlimits[bin]),avptbin);
471  grAvPtVSPtmean->SetPointError(bin-firstbin,0.5*(ptbinlimits[bin+1]-ptbinlimits[bin]),0.5*(ptbinlimits[bin+1]-ptbinlimits[bin]),errisq,errisq);
472
473
474  hSignal->SetBinContent(bin+1,hSignal->GetBinContent(bin+1)*hPtSpectra->GetBinWidth(1)/hSignal->GetBinWidth(bin+1));
475  hSignal->SetBinError(bin+1,hSignal->GetBinError(bin+1)*hPtSpectra->GetBinWidth(1)/hSignal->GetBinWidth(bin+1));
477  cPtDistrNoSubtr[bin]->Draw();
478  }
479
480  TString nameout="ptCorrection";
481  TString namegrRawAvPt="grAvRawVsAvPt",namegrAvPtPtmean="grAvPtMeanPt";
483  nameout.Append("FitSB");
484  namegrRawAvPt.Append("FitSB");
485  namegrAvPtPtmean.Append("FitSB");
486  }
487  if(!subtractSB){
488  namegrRawAvPt.Append("NoSBsub.root");
489  namegrAvPtPtmean.Append("NoSBsub.root");
490  nameout.Append("NoSBsubtract.root");
491  }
492  else nameout.Append(".root");
493
494
495  grAvRawYieldSpectrum->SetName(namegrRawAvPt.Data());
496  grAvRawYieldSpectrum->SetMarkerStyle(21);
497  grAvRawYieldSpectrum->SetMarkerSize(1.2);
498  grAvRawYieldSpectrum->SetMarkerColor(kBlue);
499  grAvRawYieldSpectrum->SetLineColor(kBlue);
500
501  grAvPtVSPtmean->SetName(namegrAvPtPtmean.Data());
502  grAvPtVSPtmean->SetMarkerStyle(21);
503  grAvPtVSPtmean->SetMarkerSize(1.2);
504  grAvPtVSPtmean->SetMarkerColor(kBlue);
505  grAvPtVSPtmean->SetLineColor(kBlue);
506
507
508
510  // TH1D *histoThCompareCentral=(TH1D*)ftheory->Get("hD0Kpipred_central");
511  // histoThCompareCentral->Rebin(5);
512
513
514
515
516  TCanvas *cPtSpectra=new TCanvas("cPtSpectra","cPtSpectra",700,700);
517  cPtSpectra->cd();
518  hPtSpectra->Draw();
519  //hAvRawYieldSpectrum->Draw("same");
520  hSignal->SetLineColor(kGreen);
521  hSignal->Draw("same");
522  grAvRawYieldSpectrum->Draw("p");
523
524  TCanvas *cAvPtVSPtmean=new TCanvas("cAvPtVSPtmean","cAvPtVSPtmean",700,700);
525  cAvPtVSPtmean->cd();
526  grAvPtVSPtmean->Draw("ap");
527  grAvPtVSPtmean->GetYaxis()->SetTitle("<p_{t}> (GeV/c)");
528  grAvPtVSPtmean->GetXaxis()->SetTitle("bin centre p_{t} (GeV/c)");
529
530  grBackAvPtVSPtmean->SetMarkerColor(kRed);
531  grBackAvPtVSPtmean->SetMarkerStyle(24);
532  grBackAvPtVSPtmean->SetLineColor(kRed);
533  grBackAvPtVSPtmean->GetYaxis()->SetTitle("<p_{t}> (GeV/c)");
534  grBackAvPtVSPtmean->GetXaxis()->SetTitle("bin centre p_{t} (GeV/c)");
535  grBackAvPtVSPtmean->Draw("p");
536
537
538  TLegend *leg=new TLegend(0.7,0.5,0.95,0.2,"","NDC");
541  leg->Draw();
542
543  // PRINTING RESULTS
544  Double_t *ygrSignal,*ygrBack,*eygrSignal,*eygrBack;
545  ygrSignal=grAvPtVSPtmean->GetY();
546  eygrSignal=grAvPtVSPtmean->GetEYhigh();
547  ygrBack=grBackAvPtVSPtmean->GetY();
548  eygrBack=grBackAvPtVSPtmean->GetEYhigh();
549  printf("Av Pt for signal: \n");
550  for(Int_t in=0;in<grAvPtVSPtmean->GetN();in++){
551  cout<<ygrSignal[in]<<endl;
552  }
553  printf("Error on Av Pt for signal: \n");
554  for(Int_t in=0;in<grAvPtVSPtmean->GetN();in++){
555  cout<<eygrSignal[in]<<endl;
556  }
557
558  printf("Av Pt for back: \n");
559  for(Int_t in=0;in<grAvPtVSPtmean->GetN();in++){
560  cout<<ygrBack[in]<<endl;
561  }
562  printf("Error on Av Pt for back: \n");
563  for(Int_t in=0;in<grAvPtVSPtmean->GetN();in++){
564  cout<<eygrBack[in]<<endl;
565  }
566
567
568
569
570  nameout="ptCorrection";
571  if(useFitForSubtraction)nameout.Append("FitSB");
572  if(!subtractSB)nameout.Append("NoSBsubtract.root");
573  else nameout.Append(".root");
574  TFile *fout=new TFile(nameout.Data(),"RECREATE");
575  fout->cd();
576  grAvPtVSPtmean->Write();
577  grAvRawYieldSpectrum->Write();
578  grBackAvPtVSPtmean->Write();
579  cPtSpectra->Write();
580  for(Int_t bin=firstbin;bin<=lastbin;bin++){
581  cPtDistrNoSubtr[bin]->Write();
582  }
583
584
585
586  fout->Close();
587
588  // histoThCompareCentral->Draw("sames");
589 }
590
591
592 void DoStandardForD0(Int_t rebin,Bool_t usefit,Bool_t corrforeff=kTRUE,Bool_t useParGenAccLimacc=kTRUE,Int_t firstbin=0,Int_t lastbin=3){// firstbin and lastbin are w.r.t. the signal histo plot, starting from 0
593
594  //N.B. ADAPT THE PART BELOW WITH YOUR PATHS AND HISTO NAMES
598  TDirectory *fdir=(TDirectory*)fDataList->Get("PWG3_D2H_d0D0");
599  TList *lslist=(TList*)fdir->Get("clistTGHCsign_d0D0");
600  TH2F *hMassPt=(TH2F*)lslist->FindObject("hInvMassPtTGHCsign");
601  TH1D *hSig=(TH1D*)fData->Get("hSigma");
602  TH1D *hBackground=(TH1D*)fData->Get("hBackground");
603  TH1D *hSign=(TH1D*)fData->Get("hSignal");
604  TH1D *hMean=(TH1D*)fData->Get("hMass");
605
606  mesonMass=1.86484;
607
608  TH1D *hEffNum=(TH1D*)fCF->Get("hRecoPIDpt");
609  TH1D *hEffDeNum=(TH1D*)fCF->Get("hMCAccpt");
610
612  // FOR DEBUGGING: CHECK EFFICIENCY PLOTS REBINNING
613  // TCanvas *cEff=new TCanvas();
614  // cEff->cd();
615  // TH1D *hRatio=new TH1D(*(TH1D*)hEffNum);
616  // hRatio->Divide(hEffDeNum);
617
618  // TH1D *hPt=(TH1D*)(hMassPt->ProjectionY());
619
620  // hPt->Rebin(rebin);cout<< hPt->GetBinWidth(1)<<endl;
621  // cout<< hEffNum->GetBinWidth(1)<<endl;
622  // TH1D *hWorkNum=0x0;
623  // hWorkNum=CheckBinningAndMerge(hPt,(TH1D*)hEffNum,0.01,23.9);
624  // TH1D *hWorkDenum=0x0;
625  // hWorkDenum=CheckBinningAndMerge(hPt,(TH1D*)hEffDeNum,0.01.,23.9);
626  // TH1D *hWorkRatio=new TH1D(*hWorkNum);
627  // hWorkRatio->Divide(hWorkDenum);
628
629
630  // hWorkRatio->Draw();
631  // hPt->Draw("sames");
632
633  // hRatio->Draw("sames");
635
636  SetNsigmaForSignal(3.);
637  SetNsigmaStartSB(5.);
638  SetUseFitForSubtraction(usefit);
639  corrForEff=corrforeff;
640  useParGenAccOverLimAcc=useParGenAccLimacc;
641  CalculateAveragePt(hMassPt,hBackground,hSig,hEffNum,hEffDeNum,hSign,hMean,rebin,firstbin,lastbin);
642
643
644
645
646 }
647
648
649 void DoStandardForDs(Int_t rebin,Bool_t usefit,Bool_t corrforeff=kTRUE,Int_t firstbin=0,Int_t lastbin=3){// firstbin and lastbin are w.r.t. the signal histo plot, starting from 0
650
651  //N.B. ADAPT THE PART BELOW WITH YOUR PATHS AND HISTO NAMES
652 TFile *fData=TFile::Open("RawYieldBoth.root");
653 TFile *fCF=TFile::Open("fileEff_Ds_CommonFramework_from_c_Enriched.root");
654
655 TH1D *hSig=(TH1D*)fData->Get("hSigma");
656 TH1D *hMean=(TH1D*)fData->Get("hMass");
657 TH1D *hBackground=(TH1D*)fData->Get("hBackground");
658 TH1D *hSign=(TH1D*)fData->Get("hSignal");
659
660 TH1D *hEffNum=(TH1D*)fCF->Get("hRecoPIDpt");
661 TH1D *hEffDeNum=(TH1D*)fCF->Get("hMCAccpt");
662  if(!hEffNum)return;
663  if(!hEffDeNum)return;
664  TFile *file=new TFile("AnalysisResults.root");
665  TDirectory *dirFile=(TDirectory*)file->Get("PWG3_D2H_InvMassDs");
666  TList *cOutput = (TList*)dirFile->Get("coutputDs");
667  TH2F *hMassPt=(TH2F*)cOutput->FindObject("hPtVsMass");
668  //TH2F *hMassPt=(TH2F*)fData->Get("hPtVsMass");
669
670  mesonMass=1.96847;
671
672  SetUseFitForSubtraction(usefit);
673  corrForEff=corrforeff;
674  CalculateAveragePt(hMassPt,hBackground,hSig,hEffNum,hEffDeNum,hSign,hMean,rebin,firstbin,lastbin);
675
676
677 }
678 void DoStandardForDplus(Int_t rebin,Bool_t usefit,Bool_t corrforeff=kTRUE,Int_t firstbin=0,Int_t lastbin=2){
679 // firstbin and lastbin are w.r.t. the signal histo plot, starting from 0
680
681  //N.B. ADAPT THE PART BELOW WITH YOUR PATHS AND HISTO NAMES
682
685  TH2F *hMassPt=(TH2F*)fData->Get("hPtVsMassTC");
686  TH1D *hSig=(TH1D*)fData->Get("hSigma");
687  TH1D *hMean=(TH1D*)fData->Get("hMean");
688  TH1D *hBackground=(TH1D*)fData->Get("hBackground");
689  TH1D *hSign=(TH1D*)fData->Get("hSignal");
690
691  TH1D *hEffNum=(TH1D*)fCF->Get("RecoPID");
692  TH1D *hEffDeNum=(TH1D*)fCF->Get("GenAcc");
693
694  mesonMass=1.86962;
695
697  // FOR DEBUGGING: CHECK EFFICIENCY PLOTS REBINNING
698  // TCanvas *cEff=new TCanvas();
699  // cEff->cd();
700  // TH1D *hRatio=new TH1D(*(TH1D*)hEffNum);
701  // hRatio->Divide(hEffDeNum);
702
703  // TH1D *hPt=(TH1D*)(hMassPt->ProjectionY());
704
705  // hPt->Rebin(rebin);cout<< hPt->GetBinWidth(1)<<endl;
706  // cout<< hEffNum->GetBinWidth(1)<<endl;
707  // TH1D *hWorkNum=0x0;
708  // hWorkNum=CheckBinningAndMerge(hPt,(TH1D*)hEffNum,0.01,23.9);
709  // TH1D *hWorkDenum=0x0;
710  // hWorkDenum=CheckBinningAndMerge(hPt,(TH1D*)hEffDeNum,0.01.,23.9);
711  // TH1D *hWorkRatio=new TH1D(*hWorkNum);
712  // hWorkRatio->Divide(hWorkDenum);
713  // hWorkRatio->Draw();
714  // hPt->Draw("sames");
715  // hRatio->Draw("sames");
717
718  SetUseFitForSubtraction(usefit);
719  corrForEff=corrforeff;
720  CalculateAveragePt(hMassPt,hBackground,hSig,hEffNum,hEffDeNum,hSign,hMean,rebin,firstbin,lastbin);
721
722
723
724
725 }
726
727
728
729
730 TH1D *SmearEffHisto(TH1D *hInput,TString name="hEffNum",Double_t maxPt=40.,Double_t step=0.1){
731  Int_t nbins=(Int_t)(maxPt/step);
732  Int_t maxBinA=hInput->GetNbinsX();
733  Int_t binBl,binBu;
734  Double_t binlea,binuea,contInMin,contInMax;
735
736  TH1D *h=new TH1D(name.Data(),name.Data(),nbins,0.,maxPt);
737  for(Int_t bin=1;bin<=maxBinA;bin++){
738  contInMin=hInput->GetBinContent(bin);
739  if(bin==maxBinA) contInMax=contInMin*1.2;
740  else contInMax=hInput->GetBinContent(bin+1);
741  binlea=hInput->GetBinLowEdge(bin);
742  binuea=hInput->GetXaxis()->GetBinUpEdge(bin);
743  binBl=h->FindBin(binlea);
744  binBu=h->FindBin(binuea);
745  for(Int_t bb=binBl;bb<binBu;bb++){
746  h->SetBinContent(bb,contInMin+(bb-binBl)*(contInMax-contInMin)/(binBu-binBl));
747  }
748  }
749  return h;
750
751 }
752
753
755
756  Double_t *yD0,*eyD0,*yDplus,*eyDplus,*x,*exLow,*exHigh;
757  Double_t error=0.,average=0.;
758
759  TFile *fD0=TFile::Open(fileD0.Data());
760  // TCanvas *cD0=(TCanvas*)fD0->Get("cAvPtVSPtmean");
761  // TGraphAsymmErrors *grD0=(TGraphAsymmErrors*)cD0->FindObject("grAvPtMeanPtFitSB");
762  TGraphAsymmErrors *grD0=(TGraphAsymmErrors*)fD0->Get("grAvPtMeanPtFitSB");
763
764  grD0->SetName("grAvPtMeanPtFitSB_D0");
765  yD0=grD0->GetY();
766  eyD0=grD0->GetEYhigh();
767
768
769  TFile *fDplus=TFile::Open(fileDplus.Data());
770  TCanvas *cDplus=(TCanvas*)fDplus->Get("cAvPtVSPtmean");
771  TGraphAsymmErrors *grDplus=(TGraphAsymmErrors*)cDplus->FindObject("grAvPtMeanPtFitSB");
772  grDplus->SetName("grAvPtMeanPtFitSB_Dplus");
773  yDplus=grDplus->GetY();
774  eyDplus=grDplus->GetEYhigh();
775
776
777  Int_t nlmax=grDplus->GetN();
778  Int_t nlD0=grD0->GetN();
779  Int_t nlDplus=grDplus->GetN();
780
781  if(grD0->GetN()>nlmax){
782  nlmax=grD0->GetN();
783  x=grD0->GetX();
784  exLow=grD0->GetEXlow();
785  exHigh=grD0->GetEXhigh();
786  }
787  else {
788  x=grDplus->GetX();
789  exLow=grDplus->GetEXlow();
790  exHigh=grDplus->GetEXhigh();
791  }
792
793
794  TGraphAsymmErrors *grAverage=new TGraphAsymmErrors();
795  grAverage->SetName("grAverageD0Dplus");
796
797  for(Int_t j=0;j<nlmax;j++){
798  if(j<nlD0&&j<nlDplus){
799  average=(yD0[j]/(eyD0[j]*eyD0[j])+yDplus[j]/(eyDplus[j]*eyDplus[j]))/(1./(eyD0[j]*eyD0[j])+1./(eyDplus[j]*eyDplus[j]));
800  error=TMath::Sqrt(1./(1./(eyD0[j]*eyD0[j])+1./(eyDplus[j]*eyDplus[j])));
801  }
802  else if(j>=nlD0){
803  average=yDplus[j];
804  error=eyDplus[j];
805  }
806  else if(j>=nlDplus){
807  average=yD0[j];
808  error=eyD0[j];
809  }
810  grAverage->SetPoint(j,x[j],average);
811  grAverage->SetPointError(j,exLow[j],exHigh[j],error,error);
812  }
813
814  // Printing Info: D0
815  printf("Average Value D0: \n");
816  for(Int_t j=0;j<nlD0;j++){
817  cout<<yD0[j]<<endl;
818
819  }
820
821  printf("Average Errors D0: \n");
822  for(Int_t j=0;j<nlD0;j++){
823  cout<<eyD0[j]<<endl;
824
825  }
826
827  // Printing Info: Dplus
828  printf("Average Value Dplus: \n");
829  for(Int_t j=0;j<nlDplus;j++){
830  cout<<yDplus[j]<<endl;
831
832  }
833
834  printf("Average Errors Dplus: \n");
835  for(Int_t j=0;j<nlDplus;j++){
836  cout<<eyDplus[j]<<endl;
837
838  }
839
840
841  // Printing Info: Dplus D0 average
842  yDplus=grAverage->GetY();
843  eyDplus=grAverage->GetEYhigh();
844  printf("Average Value D0 + Dplus: \n");
845  for(Int_t j=0;j<nlDplus;j++){
846  cout<<yDplus[j]<<endl;
847
848  }
849
850  printf("Average Errors D0+Dplus: \n");
851  for(Int_t j=0;j<nlDplus;j++){
852  cout<<eyDplus[j]<<endl;
853
854  }
855
856
857
858
859  TCanvas *cCompare=new TCanvas("cCompareD0Dplus","cCompareD0Dplus",700,700);
860  cCompare->cd();
861
862  grDplus->SetMarkerStyle(21);
863  grDplus->SetMarkerSize(1.2);
864  grDplus->SetMarkerColor(kBlue);
865  grDplus->SetLineColor(kBlue);
866  grDplus->Draw("ap");
867
868  grD0->SetMarkerStyle(22);
869  grD0->SetMarkerSize(1.2);
870  grD0->SetMarkerColor(kRed);
871  grD0->SetLineColor(kRed);
872  grD0->Draw("p");
873
874  grAverage->SetMarkerStyle(20);
875  grAverage->SetMarkerSize(1.2);
876  grAverage->SetMarkerColor(kBlack);
877  grAverage->SetLineColor(kBlack);
878  grAverage->Draw("p");
879
880
881
882 }
TH1D * hSignal
void SetPtBinLimits(const Int_t npt, Double_t *ptbinlim)
Int_t standrebin[8]
TH1D * hEfficDenum
TH1D * hSigma
Bool_t useFitForSubtraction
void SetHistoMassPt(TH2F *h2)
TGraphAsymmErrors * grBackAvPtVSPtmean
Int_t nbinsy
Double_t nsigmaSBstart
void SetSubtractSB(Bool_t subtract)
TCanvas ** cPtDistrNoSubtr
void CalculateAveragePt(Int_t rebin=1, Int_t firstbin=0, Int_t lastbin=0)
Double_t binwidthInvMss
Double_t * meansignal
Double_t nsigmaSignal
Double_t * ptbinlimits
Bool_t useParGenAccOverLimAcc
void SetNsigmaStartSB(Double_t nsigm)
TH1D * hAvRawYieldSpectrum
void SetHistRawBack(TH1D *hB)
Double_t mesonMass
void SetUseFitForSubtraction(Bool_t useFit)
TH1D * CheckBinningAndMerge(TH1D *hA, TH1D *hB, Double_t precision=0.001, Double_t minX=-9999., Double_t maxX=-9999.)
Bool_t subtractSB
Double_t ptmax
Double_t * sigma
TH1D * hEfficNum
void SetNsigmaForSignal(Double_t nsigm)
void DoStandardForDplus(Int_t rebin, Bool_t usefit, Bool_t corrforeff=kTRUE, Int_t firstbin=0, Int_t lastbin=2)
TF1 * ParametricGenAccOverLimAccCorr()
void DoStandardForDs(Int_t rebin, Bool_t usefit, Bool_t corrforeff=kTRUE, Int_t firstbin=0, Int_t lastbin=3)
TH1D * hMeanSignal
TH2F * hPtInvMass
Double_t * rawback
void SetHistRawSignal(TH1D *hS)
Double_t binwidthpt
Bool_t CorrectForEfficiency(TH1D *hPtHisto)
TH1D * SmearEffHisto(TH1D *hInput, TString name="hEffNum", Double_t maxPt=40., Double_t step=0.1)
Int_t nbinsx
void SetHistMean(TH1D *hM)
void SetHistSigma(TH1D *hSig)
Double_t ptmin
TH1D * HistoPtShapeFromData(Int_t ptbin, Int_t rebin=1.)
TGraphAsymmErrors * grAvPtVSPtmean
Double_t * rawsignal
void SetCorrForEff(Bool_t correff)
TFile * file
Int_t rebin
TF1 ** fitfunc
const Int_t nbins
TH1D * hBack
Bool_t corrForEff
Int_t nptbins
void SetHistosEfficiency(TH1D *hNum, TH1D *hDenum)
TGraphAsymmErrors * grAvRawYieldSpectrum
void DoStandardForD0(Int_t rebin, Bool_t usefit, Bool_t corrforeff=kTRUE, Bool_t useParGenAccLimacc=kTRUE, Int_t firstbin=0, Int_t lastbin=3)