]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/vertexingHF/macros/CalculateAveragePt.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / macros / CalculateAveragePt.C
CommitLineData
ea0c35d7 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//
39////////////////////////////////////////////////////////////////////////////////////
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
18ceba37 55TH2F *hPtInvMass=0x0;
56Int_t nptbins=0;
57Double_t *ptbinlimits=0x0;
ea0c35d7 58Double_t *rawsignal=0x0,*rawback=0x0,*sigma=0x0,*meansignal=0x0;
59Double_t nsigmaSignal=3.,nsigmaSBstart=5.;
60Double_t mesonMass=1.8645;
18ceba37 61Int_t nbinsx;//=hPtInvMass->GetNbinsX();
62Int_t nbinsy;//=hPtInvMass->GetNbinsY();
63Double_t binwidthpt;//=hPtInvMass->GetYaxis()->GetBinWidth(1);// ASSUME ALL BINS ARE WITH THE SAME WIDTH
64Double_t binwidthInvMss;//=hPtInvMass->GetXaxis()->GetBinWidth(1);// ASSUME ALL BINS ARE WITH THE SAME WIDTH
65Double_t ptmin;//=hPtInvMass->GetYaxis()->GetBinLowEdge(1);
66Double_t ptmax;//=hPtInvMass->GetYaxis()->GetBinUpEdge(nbinsy);
67Bool_t useFitForSubtraction=kFALSE;
ea0c35d7 68Bool_t useParGenAccOverLimAcc=kFALSE;
18ceba37 69Bool_t subtractSB=kTRUE;
ea0c35d7 70Bool_t corrForEff=kFALSE;
18ceba37 71void SetSubtractSB(Bool_t subtract){subtractSB=subtract;}
ea0c35d7 72void SetCorrForEff(Bool_t correff){corrForEff=correff;}
18ceba37 73void SetUseFitForSubtraction(Bool_t useFit){useFitForSubtraction=useFit;}
ea0c35d7 74void SetNsigmaForSignal(Double_t nsigm){nsigmaSignal=nsigm;}
75void SetNsigmaStartSB(Double_t nsigm){nsigmaSBstart=nsigm;
76}
77void CalculateAveragePt(Int_t rebin=1,Int_t firstbin=0,Int_t lastbin=0);
78TF1* 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
84Int_t standrebin[8]={2,2,2,2,2,2,4,5};
18ceba37 85TCanvas **cPtDistrNoSubtr;
86TF1 **fitfunc;
ea0c35d7 87TH1D *hAvRawYieldSpectrum,*hSignal,*hMeanSignal,*hEfficNum,*hEfficDenum,*hBack,*hSigma;
88void SetHistosEfficiency(TH1D *hNum,TH1D *hDenum){
89 hEfficNum=(TH1D*)hNum;
90 hEfficDenum=(TH1D*)hDenum;
91
92}
18ceba37 93TGraphAsymmErrors *grAvRawYieldSpectrum;
94TGraphAsymmErrors *grAvPtVSPtmean;
95
ea0c35d7 96TGraphAsymmErrors *grBackAvPtVSPtmean;
97
98void SetPtBinLimits(const Int_t npt,Double_t *ptbinlim){
18ceba37 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}
ea0c35d7 110
18ceba37 111void 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
126void SetHistRawSignal(TH1D *hS){hSignal=hS;
127 rawsignal=&(hSignal->GetArray())[1];
128 return;}
ea0c35d7 129
130void SetHistMean(TH1D *hM){hMeanSignal=hM;
131 meansignal=&(hMeanSignal->GetArray())[1];
132 return;
133}
134
18ceba37 135void 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}
144void SetHistSigma(TH1D *hSig){hSigma=hSig;
145 sigma=&(hSigma->GetArray())[1];
146 return;
147}
148
149void 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);
ea0c35d7 157
158}
18ceba37 159
ea0c35d7 160TH1D* 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
18ceba37 235}
236
ea0c35d7 237Bool_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");
254 if(useParGenAccOverLimAcc){
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}
18ceba37 270
ea0c35d7 271TH1D *HistoPtShapeFromData(Int_t ptbin,Int_t rebin=1.){
18ceba37 272 if(!hPtInvMass){
273 printf("No histo set \n");
ea0c35d7 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);
18ceba37 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;
ea0c35d7 290 TH1D *hPtDistrSB=new TH1D(namehist.Data(),"Side band Pt distribution",nbinsy,ptmin,ptmax);
18ceba37 291 namehist="hPtDistrTotPeak_Bin";
292 namehist+=ptbin;
293 printf("After Setting SB hist \n");
ea0c35d7 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 }
18ceba37 301
ea0c35d7 302
303
18ceba37 304 printf("Before loop on xbins \n");
305 for(Int_t j=1;j<nbinsx;j++){
ea0c35d7 306 if((hPtInvMass->GetXaxis()->GetBinLowEdge(j)<=mean-nsigmaSBstart*sigma[ptbin])||(hPtInvMass->GetXaxis()->GetBinUpEdge(j)>=mean+nsigmaSBstart*sigma[ptbin])){
18ceba37 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));
ea0c35d7 310 hPtDistrSB->SetBinError(k,TMath::Sqrt(hPtDistrSB->GetBinError(k)*hPtDistrSB->GetBinError(k)+hPtInvMass->GetBinError(j,k)*hPtInvMass->GetBinError(j,k)));
18ceba37 311 }
312 }
313 }
ea0c35d7 314 else if((hPtInvMass->GetXaxis()->GetBinLowEdge(j)>=mean-nsigmaSignal*sigma[ptbin])&&(hPtInvMass->GetXaxis()->GetBinUpEdge(j)<=mean+nsigmaSignal*sigma[ptbin])){
18ceba37 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));
ea0c35d7 319 hPtDistrTotPeak->SetBinError(k,TMath::Sqrt(hPtDistrTotPeak->GetBinError(k)*hPtDistrTotPeak->GetBinError(k)+hPtInvMass->GetBinError(j,k)*hPtInvMass->GetBinError(j,k)));
18ceba37 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 }
ea0c35d7 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);
18ceba37 360
ea0c35d7 361 TH1D *hPtSignalFromSubtraction=new TH1D(*hPtDistrTotPeak);
18ceba37 362 namehist="hPtSignalFromSubtraction_Bin";
363 namehist+=ptbin;
364 hPtSignalFromSubtraction->SetName(namehist.Data());
365 hPtSignalFromSubtraction->SetTitle("Signal Pt distr from subtraction");
ea0c35d7 366
367
18ceba37 368 if(subtractSB){
ea0c35d7 369 if(!useFitForSubtraction)hPtSignalFromSubtraction->Add(hPtDistrSB,-rawback[ptbin]*(nsigmaSignal/3.)/hPtDistrSB->Integral());// ROUGH LINEAR SCALING OF BACKGOUND!!
18ceba37 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++){
ea0c35d7 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));
18ceba37 382 }
383 }
384 }
385
386
387 cPtDistrNoSubtr[ptbin]->Divide(1,2);
388 cPtDistrNoSubtr[ptbin]->cd(1);
389 hPtDistrSB->Draw();
390 cPtDistrNoSubtr[ptbin]->cd(2);
ea0c35d7 391 if(corrForEff){
392 if(!CorrectForEfficiency(hPtSignalFromSubtraction)){
393 printf("Correction for efficieny failed \n");
394 }
395 }
18ceba37 396 hPtDistrTotPeak->Draw();
397 hPtSignalFromSubtraction->Draw("sames");
398
399
400 return hPtSignalFromSubtraction;
401
402}
403
ea0c35d7 404void 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){
18ceba37 405
406 SetPtBinLimits(hB);
407 SetHistoMassPt((TH2F*)hMassPt);
408 SetHistRawBack(hB);
409 SetHistSigma(hSigm);
ea0c35d7 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 }
18ceba37 418 if(hS)SetHistRawSignal(hS);
ea0c35d7 419 if(hMean)SetHistMean(hMean);
18ceba37 420 printf("Histos set \n");
421 CalculateAveragePt(rebin,firstbin,lastbin);
422}
423
ea0c35d7 424void 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
18ceba37 438 grAvRawYieldSpectrum=new TGraphAsymmErrors();
439 grAvPtVSPtmean=new TGraphAsymmErrors();
ea0c35d7 440 grBackAvPtVSPtmean=new TGraphAsymmErrors();
441 grBackAvPtVSPtmean->SetName("grBackAvPtVSPtmean");
18ceba37 442
ea0c35d7 443
444 Double_t avptbin=0.,errA,errB,errC,errisq;
18ceba37 445 for(Int_t bin=firstbin;bin<=lastbin;bin++){
446 // printf("Before Getting Subtracted Histo \n");
447 avptbin=0.;
ea0c35d7 448 errA=0.;
449 errB=0.;
450 errC=0.;
451 errisq=0.;
452
453 TH1D *h=HistoPtShapeFromData(bin,rebin);
18ceba37 454 hPtSpectra->Add(h);
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);
ea0c35d7 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);
18ceba37 463 }
464 avptbin/=h->Integral("width");
ea0c35d7 465 errisq=TMath::Sqrt(errA-2.*avptbin*errB+avptbin*avptbin*errC)/h->Integral("width");
466
18ceba37 467 grAvRawYieldSpectrum->SetPoint(bin-firstbin,avptbin,h->Integral("width")/hAvRawYieldSpectrum->GetBinWidth(bin+1));
ea0c35d7 468 grAvRawYieldSpectrum->SetPointError(bin-firstbin,avptbin-ptbinlimits[bin],ptbinlimits[bin+1]-avptbin,errisq/2.,errisq/2.);
18ceba37 469
470 grAvPtVSPtmean->SetPoint(bin-firstbin,0.5*(ptbinlimits[bin+1]+ptbinlimits[bin]),avptbin);
ea0c35d7 471 grAvPtVSPtmean->SetPointError(bin-firstbin,0.5*(ptbinlimits[bin+1]-ptbinlimits[bin]),0.5*(ptbinlimits[bin+1]-ptbinlimits[bin]),errisq,errisq);
18ceba37 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));
476 // hPtSpectra->Add(HistoPtShapeFromData(bin,rebin),1./efficiency[bin]/nev*factorTheoryCompare*primaryCharmCorrect[bin]/binwidthpt);//(ptbinlimits[bin+1]-ptbinlimits[bin])
477 cPtDistrNoSubtr[bin]->Draw();
478 }
479
480 TString nameout="ptCorrection";
481 TString namegrRawAvPt="grAvRawVsAvPt",namegrAvPtPtmean="grAvPtMeanPt";
482 if(useFitForSubtraction){
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
ea0c35d7 509 // TFile *ftheory=TFile::Open("/Users/administrator/ALICE/CHARM/ppData_2010/PtSelInBin/predictions_D0.root");
510 // TH1D *histoThCompareCentral=(TH1D*)ftheory->Get("hD0Kpipred_central");
511 // histoThCompareCentral->Rebin(5);
18ceba37 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");
ea0c35d7 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");
539 leg->AddEntry(grAvPtVSPtmean,"signal","lp");
540 leg->AddEntry(grBackAvPtVSPtmean,"background","lp");
541 leg->Draw();
18ceba37 542
ea0c35d7 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";
18ceba37 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();
ea0c35d7 578 grBackAvPtVSPtmean->Write();
579 cPtSpectra->Write();
580 for(Int_t bin=firstbin;bin<=lastbin;bin++){
581 cPtDistrNoSubtr[bin]->Write();
582 }
583
584
585
18ceba37 586 fout->Close();
587
588 // histoThCompareCentral->Draw("sames");
589}
590
591
ea0c35d7 592void 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
595 TFile *fData=TFile::Open("/Users/administrator/ALICE/CHARM/PbPBdata_2011/TestTrain/2013June4TrainData92_MC61/Data/Standard/RAAvsNPart/FinalMassPlots_v2/RawYieldBoth_tight.root");
596 TFile *fCF=TFile::Open("/Users/administrator/ALICE/CHARM/PbPBdata_2011/2013Jun8MCptWeightDataMoreSplit/MergeWithMyCode/Standard/EffPrompt/fileEff_D0_from_c.root");//EffForHFspectrumLHC10f6NOTd4Corr.root");
597 TFile *fDataList=TFile::Open("/Users/administrator/ALICE/CHARM/PbPBdata_2011/TestTrain/2013June4TrainData92_MC61/Data/AnalysisResults.root");
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
611 /////////////////////////////////////////////////////////////////
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);
18ceba37 617
ea0c35d7 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");
634 /////////////////////////////////////////////////////////////////
18ceba37 635
ea0c35d7 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
18ceba37 643
644
ea0c35d7 645
646}
647
648
649void 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
652TFile *fData=TFile::Open("RawYieldBoth.root");
653TFile *fCF=TFile::Open("fileEff_Ds_CommonFramework_from_c_Enriched.root");
654
655TH1D *hSig=(TH1D*)fData->Get("hSigma");
656TH1D *hMean=(TH1D*)fData->Get("hMass");
657TH1D *hBackground=(TH1D*)fData->Get("hBackground");
658TH1D *hSign=(TH1D*)fData->Get("hSignal");
659
660TH1D *hEffNum=(TH1D*)fCF->Get("hRecoPIDpt");
661TH1D *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}
678void 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
683 TFile *fData=TFile::Open("/Users/administrator/ALICE/CHARM/ppData_2010/PtSelInBin/testWithEff/Dplus/RawYieldBoth.root");
684 TFile *fCF=TFile::Open("/Users/administrator/ALICE/CHARM/ppData_2010/PtSelInBin/testWithEff/Dplus/CFEfficiency.root");
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
696 /////////////////////////////////////////////////////////////////
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");
716 /////////////////////////////////////////////////////////////////
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
730TH1D *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
754void AverageD0DplusResults(TString fileD0="/Users/administrator/ALICE/CHARM/ppData_2010/2011_Jul_05/data/LHC10bcdeAOD057/AvPt/MassRegSel3PkMore5SBEffCorrMeanFit/ptCorrectionFitSB.root",TString fileDplus="/Users/administrator/ALICE/CHARM/ppData_2010/2011_Jul_05/Dplus/AvPt/2011Jul26Renu/average_ptNew.root"){
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}