]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/vertexingHF/macros/ReadDvsMultiplicity.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Adding Tab values for QpPb analysis
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / macros / ReadDvsMultiplicity.C
CommitLineData
0987c810 1#if !defined(__CINT__) || defined(__MAKECINT__)
2#include <TInterpreter.h>
3#include <TString.h>
4#include <TObjString.h>
5#include <TObjArray.h>
6#include <TMath.h>
7#include <TFile.h>
8#include <TCanvas.h>
9#include <TH1.h>
10#include <TH1F.h>
11#include <TH2F.h>
12#include <TH3.h>
13#include <TH3F.h>
14#include <TH1D.h>
15#include <TF1.h>
16#include <TSystem.h>
17#include <TStyle.h>
18#include <TLegend.h>
19#include <TList.h>
20#include <TLegendEntry.h>
21#include <TDatabasePDG.h>
22#include <TGraph.h>
23
24#include "AliAODRecoDecayHF.h"
25#include "AliRDHFCuts.h"
26#include "AliRDHFCutsDplustoKpipi.h"
27#include "AliRDHFCutsD0toKpi.h"
28#include "AliHFMassFitter.h"
29#include "AliNormalizationCounter.h"
30#endif
31
32/* $Id$ */
33
34enum {kD0toKpi, kDplusKpipi};
35enum {kBoth, kParticleOnly, kAntiParticleOnly};
36enum {kExpo=0, kLinear, kPol2};
37enum {kGaus=0, kDoubleGaus};
38enum {kCorr=0, kUnCorr, kNoPid};
39
40
41// Common variables: to be configured by the user
42// const Int_t nPtBins=6;
43// Double_t ptlims[nPtBins+1]={1., 2.,4.,6.,8.,12.,24.};
44// Int_t rebin[nPtBins]={4,4,6,6,8,8};
45// Double_t sigmapt[nPtBins]={ 0.008, 0.010, 0.012, 0.016, 0.018, 0.020 };
46const Int_t nPtBins=5;
47Double_t ptlims[nPtBins+1]={1.,2.,4.,8.,12.,24.};
48Int_t rebin[nPtBins]={4,4,6,8,8};
49Double_t sigmapt[nPtBins]={ 0.008, 0.014, 0.019, 0.027, 0.033 };
50Bool_t fixPeakSigma = kTRUE;
51//
52const Int_t nMultbins=7;
53Double_t multlims[nMultbins+1]={1.,9.,14.,20.,31.,50.,81.,100.};
54// const Int_t nMultbins=1;
55// Double_t multlims[nMultbins+1]={0.,500.};
56//
57Int_t firstUsedBin[nPtBins]={-1,-1,-1,-1,-1}; // -1 uses all bins, >=1 to set the lower bin to be accepted from original histo
58//
59Bool_t isMC=kFALSE;
60Int_t typeb=kExpo;
61Int_t types=kGaus;
62Int_t optPartAntiPart=kBoth;
63Int_t factor4refl=0;
64Float_t massRangeForCounting=0.05; // GeV --> it is 3 sigmapt[binpt]
8f1140d0 65Float_t nSigmaRangeForCounting=3.0; // 3 sigmapt[binpt]
0987c810 66TH2F* hPtMass=0x0;
67TString suffix="StdPid";
68
69
70// Functions
71Bool_t LoadDplusHistos(TObjArray* listFiles, TH3F** hPtMassMult, TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, const char *CutsType, Int_t Option);
72Bool_t LoadD0toKpiHistos(TObjArray* listFiles, TH3F** hPtMassMult, TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, AliNormalizationCounter *counter, const char *CutsType, Int_t Option);
73Bool_t CheckNtrVsZvtx(TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, Int_t nFiles);
74TH1F* RebinHisto(TH1F* hOrig, Int_t reb, Int_t firstUse=-1);
75
76
77void ReadDvsMultiplicity(Int_t analysisType=kD0toKpi,
78 TString fileNameb="AnalysisResults.root",
79 TString fileNamec="",
80 TString fileNamed="",
81 TString fileNamee="",
82 const char *CutsType="",
83 Int_t Option=kCorr)
84{
85 // gSystem->SetIncludePath("-I. -I$ROOTSYS/include -I$ALICE_ROOT -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -I$ALICE_ROOT/TPC -I$ALICE_ROOT/CONTAINERS -I$ALICE_ROOT/STEER/STEER -I$ALICE_ROOT/STEER/STEERBase -I$ALICE_ROOT/STEER/ESD -I$ALICE_ROOT/STEER/AOD -I$ALICE_ROOT/TRD -I$ALICE_ROOT/macros -I$ALICE_ROOT/ANALYSIS -I$ALICE_ROOT/OADB -I$ALICE_ROOT/PWGHF -I$ALICE_ROOT/PWGHF/base -I$ALICE_ROOT/PWGHF/vertexingHF -I$ALICE_ROOT/PWG/FLOW/Case -I$ALICE_ROOT/PWG/FLOW/Tasks -g");
86
87 // gInterpreter->ExecuteMacro("$ALICE_ROOT/PWGHF/vertexingHF/macros/LoadLibraries.C");
88 gStyle->SetOptTitle(1);
89
90 // for(int j=0; j<=nMultbins; j++) multlims[j] *= (68./8.8);
91
92 Int_t nFiles=0;
93 TObjArray* listFiles=new TObjArray();
94 if(fileNameb!="") { listFiles->AddLast(new TObjString(fileNameb.Data())); nFiles++; }
95 if(fileNamec!="") { listFiles->AddLast(new TObjString(fileNamec.Data())); nFiles++; }
96 if(fileNamed!="") { listFiles->AddLast(new TObjString(fileNamed.Data())); nFiles++; }
97 if(fileNamee!="") { listFiles->AddLast(new TObjString(fileNamee.Data())); nFiles++; }
98 if(listFiles->GetEntries()==0){
99 printf("Missing file names in input\n");
100 return;
101 }
102 TH3F** hPtMassMult=new TH3F*[1];
103 hPtMassMult[0]=0x0;
104 TH1F** hmass=new TH1F*[nPtBins*nMultbins];
105 for(Int_t i=0; i<nPtBins*nMultbins; i++) hmass[i]=0x0;
106 TH2F** hNtrZvtx=new TH2F*[4];
107 for(Int_t i=0; i<4; i++) hNtrZvtx[i]=0x0;
108 TH2F** hNtrZvtxCorr=new TH2F*[4];
109 for(Int_t i=0; i<4; i++) hNtrZvtxCorr[i]=0x0;
110 AliNormalizationCounter *counter=0x0;
111 TH1F * hNormalization = new TH1F("hNormalization","Events in the norm counter",nMultbins,multlims);
112
113 Float_t massD;
114 Bool_t retCode;
115 if(analysisType==kD0toKpi) {
116 massD=1.86484003067016602e+00;//(Float_t)(TDatabasePDG::Instance()->GetParticle(421)->Mass());
117 retCode=LoadD0toKpiHistos(listFiles,hPtMassMult,hNtrZvtx,hNtrZvtxCorr,counter,CutsType,Option);
118 }
119 else if(analysisType==kDplusKpipi) {
120 massD=1.86961996555328369e+00;//(Float_t)(TDatabasePDG::Instance()->GetParticle(411)->Mass());
121 retCode=LoadDplusHistos(listFiles,hPtMassMult,hNtrZvtx,hNtrZvtxCorr,CutsType,Option);
122 }
123 else{
124 printf("Wrong analysis type parameter\n");
125 return;
126 }
127 if(!retCode){
128 printf("ERROR in reading input files\n");
129 return;
130 }
131 //
132 // Check the multiplicity variables
133 //
134 Bool_t isMult = CheckNtrVsZvtx(hNtrZvtx,hNtrZvtxCorr,nFiles);
135 if(!isMult) return;
136 //
137 //
138 printf(" Preparing the mass fits");
139 TH1F* hCntSig1[nMultbins];
140 TH1F* hCntSig2[nMultbins];
141 TH1F* hNDiffCntSig1[nMultbins];
142 TH1F* hNDiffCntSig2[nMultbins];
143 TH1F* hSignal[nMultbins];
144 TH1F* hRelErrSig[nMultbins];
145 TH1F* hInvSignif[nMultbins];
146 TH1F* hBackground[nMultbins];
147 TH1F* hBackgroundNormSigma[nMultbins];
148 TH1F* hSignificance[nMultbins];
149 TH1F* hMass[nMultbins];
150 TH1F* hSigma[nMultbins];
151 for(Int_t i=0; i<nMultbins; i++){
152 hCntSig1[i]=new TH1F(Form("hCntSig1_%d",i),Form("hCntSig1_%d",i),nPtBins,ptlims);
153 hCntSig2[i]=new TH1F(Form("hCntSig2_%d",i),Form("hCntSig2_%d",i),nPtBins,ptlims);
154 hNDiffCntSig1[i]=new TH1F(Form("hNDiffCntSig1_%d",i),Form("hNDiffCntSig1_%d",i),nPtBins,ptlims);
155 hNDiffCntSig2[i]=new TH1F(Form("hNDiffCntSig2_%d",i),Form("hNDiffCntSig2_%d",i),nPtBins,ptlims);
156 hSignal[i]=new TH1F(Form("hSignal_%d",i),Form("hSignal_%d",i),nPtBins,ptlims);
157 hRelErrSig[i]=new TH1F(Form("hRelErrSig_%d",i),Form("hRelErrSig_%d",i),nPtBins,ptlims);
158 hInvSignif[i]=new TH1F(Form("hInvSignif_%d",i),Form("hInvSignif_%d",i),nPtBins,ptlims);
159 hBackground[i]=new TH1F(Form("hBackground_%d",i),Form("hBackground_%d",i),nPtBins,ptlims);
160 hBackgroundNormSigma[i]=new TH1F(Form("hBackgroundNormSigma_%d",i),Form("hBackgroundNormSigma_%d",i),nPtBins,ptlims);
161 hSignificance[i]=new TH1F(Form("hSignificance_%d",i),Form("hSignificance_%d",i),nPtBins,ptlims);
162 hMass[i]=new TH1F(Form("hMass_%d",i),Form("hMass_%d",i),nPtBins,ptlims);
163 hSigma[i]=new TH1F(Form("hSigma_%d",i),Form("hSigma_%d",i),nPtBins,ptlims);
164 }
165 printf(", defined...\n");
166
167 // std::cout << " htemp :"<<hPtMassMult[0]<<std::endl;
168 TH1F* hptaxis = (TH1F*)hPtMassMult[0]->ProjectionZ("hptaxis");
169 TH1F* hmassaxis = (TH1F*)hPtMassMult[0]->ProjectionY("hmassaxis");
170 TH1F* hmultaxis = (TH1F*)hPtMassMult[0]->ProjectionX("hmultaxis");
171 Int_t nMassBins=hmassaxis->GetNbinsX();
172 Double_t hmin=hmassaxis->GetBinLowEdge(3);
173 Double_t hmax=hmassaxis->GetBinLowEdge(nMassBins-2) + hmassaxis->GetBinWidth(nMassBins-2);
174 Int_t iPad=1;
175
176 printf("Now initializing the fit functions\n");
177 // TF1* funBckStore1=0x0;
178 // TF1* funBckStore2=0x0;
179 // TF1* funBckStore3=0x0;
180
181 Int_t nPtMultbins = nPtBins*nMultbins;
182 AliHFMassFitter** fitter=new AliHFMassFitter*[nPtMultbins];
183 Double_t arrchisquare0[nPtBins];
184 Double_t arrchisquare1[nPtBins];
185 Double_t arrchisquare2[nPtBins];
186 Double_t arrchisquare3[nPtBins];
187 Double_t arrchisquare4[nPtBins];
188 Double_t arrchisquare5[nPtBins];
189 for(Int_t i=0; i<nPtBins; i++){
190 arrchisquare0[i]=0.;
191 arrchisquare1[i]=0.;
192 arrchisquare2[i]=0.;
193 arrchisquare3[i]=0.;
194 arrchisquare4[i]=0.;
195 arrchisquare5[i]=0.;
196 }
197
198 TCanvas* canvas[nMultbins];
199 Int_t nx=2, ny=2;
200 if(nPtBins>4){
201 nx=3;
202 ny=2;
203 }
204 if(nPtBins>6){
205 nx=4;
206 ny=3;
207 }
208 if(nPtBins>12){
209 nx=5;
210 ny=4;
211 }
212 for(Int_t i=0; i<nMultbins; i++ ){
213 canvas[i] = new TCanvas(Form("c%d",i),Form("summary canvas for mult bin %d",i));
214 canvas[i]->Divide(nx,ny);
215 }
216 TCanvas *myCanvas[nPtMultbins];
217
218 //
219 //
220 // Loop on multiplicity bins
221 //
222 Int_t massBin=0;
223 for(Int_t j=0; j<nMultbins; j++){
224 Double_t sig,errsig,s,errs,b,errb;
225 // printf(" Studying multiplicity bin %d\n",j);
226 Int_t multbinlow = hmultaxis->FindBin(multlims[j]);
227 Int_t multbinhigh = hmultaxis->FindBin(multlims[j+1])-1;
228 Float_t val = multbinlow + (multbinhigh-multbinlow)/2.;
229 Int_t hnbin = hNormalization->FindBin(val);
230 Float_t nevents = 0.;
231 if(counter) { nevents = counter->GetNEventsForNorm(multbinlow,multbinhigh); cout << endl<<endl<<" Nevents ("<<multbinlow<<","<<multbinhigh<<") ="<<nevents<<endl<<endl<<endl;}
232 hNormalization->SetBinContent(hnbin,nevents);
233 //
234 // Loop on pt bins
235 //
236 iPad=1;
237 for(Int_t iBin=0; iBin<nPtBins; iBin++){
238 canvas[j]->cd(iPad++);
239 // printf(" projecting to the mass histogram\n");
240 Int_t ptbinlow = hptaxis->FindBin(ptlims[iBin]);
241 Int_t ptbinhigh = hptaxis->FindBin(ptlims[iBin+1])-1;
242 hmass[massBin] = (TH1F*)hPtMassMult[0]->ProjectionY(Form("hmass_%d_%d",j,iBin),multbinlow,multbinhigh,ptbinlow,ptbinhigh);
243 if( hmass[massBin]->GetEntries() < 60 ) {
244 massBin++;
245 continue;
246 }
247 Int_t origNbins=hmass[massBin]->GetNbinsX();
248 // printf(" rebinning the mass histogram\n");
249 TH1F* hRebinned=RebinHisto(hmass[massBin],rebin[iBin],firstUsedBin[iBin]);
250 hmin=hRebinned->GetBinLowEdge(2);
251 hmax=hRebinned->GetBinLowEdge(hRebinned->GetNbinsX());
252 // printf(" define the mass fitter and options \n");
253 fitter[massBin] = new AliHFMassFitter(hRebinned,hmin, hmax,1,typeb,types);
254 fitter[massBin]->SetRangeFit(1.65,2.10);
255 Int_t rebinItem = origNbins/fitter[massBin]->GetBinN();
256 fitter[massBin]->SetReflectionSigmaFactor(factor4refl);
257 fitter[massBin]->SetInitialGaussianMean(massD);
258 fitter[massBin]->SetInitialGaussianSigma(sigmapt[iBin]);
259 if(fixPeakSigma) {
260 fitter[massBin]->SetFixGaussianMean(massD);
261 fitter[massBin]->SetFixGaussianSigma(sigmapt[iBin],kTRUE);
262 }
263 Bool_t out=fitter[massBin]->MassFitter(0);
264 if(!out) {
265 fitter[massBin]->GetHistoClone()->Draw();
266 massBin++;
267 continue;
268 }
269 // printf(" getting the fit parameters\n");
270 Double_t mass=fitter[massBin]->GetMean();
271 Double_t massUnc=fitter[massBin]->GetMeanUncertainty();
272 Double_t sigma=fitter[massBin]->GetSigma();
273 Double_t sigmaUnc=fitter[massBin]->GetSigmaUncertainty();
274 if(j==0) arrchisquare0[iBin]=fitter[massBin]->GetReducedChiSquare();
275 else if(j==1) arrchisquare1[iBin]=fitter[massBin]->GetReducedChiSquare();
276 else if(j==2) arrchisquare2[iBin]=fitter[massBin]->GetReducedChiSquare();
277 else if(j==3) arrchisquare3[iBin]=fitter[massBin]->GetReducedChiSquare();
278 else if(j==4) arrchisquare4[iBin]=fitter[massBin]->GetReducedChiSquare();
279 else if(j==5) arrchisquare5[iBin]=fitter[massBin]->GetReducedChiSquare();
280 TF1* fB1=fitter[massBin]->GetBackgroundFullRangeFunc();
281 TF1* fB2=fitter[massBin]->GetBackgroundRecalcFunc();
282 TF1* fM=fitter[massBin]->GetMassFunc();
283 // if(iBin==0 && fB1) funBckStore1=new TF1(*fB1);
284 // if(iBin==0 && fB2) funBckStore2=new TF1(*fB2);
285 // if(iBin==0 && fM) funBckStore3=new TF1(*fM);
286
287 fitter[massBin]->DrawHere(gPad);
288 fitter[massBin]->Signal(3,s,errs);
289 fitter[massBin]->Background(3,b,errb);
290 fitter[massBin]->Significance(3,sig,errsig);
291 Double_t ry=fitter[massBin]->GetRawYield();
292 Double_t ery=fitter[massBin]->GetRawYieldError();
293 myCanvas[massBin] = new TCanvas(Form("myCanvas_%d_%d",j,iBin),Form("Invariant mass mult bin %d, pt bin %d",j,iBin));
294 fitter[massBin]->DrawHere(gPad);
295
296 Float_t cntSig1=0.;
297 Float_t cntSig2=0.;
298 Float_t cntErr=0.;
8f1140d0 299 massRangeForCounting = nSigmaRangeForCounting*sigmapt[iBin];
0987c810 300 // cout << " pt bin "<< iBin << " mass range = "<< massRangeForCounting<<endl;
301 Float_t minBinSum=hmassaxis->FindBin(massD-massRangeForCounting);
302 Float_t maxBinSum=hmassaxis->FindBin(massD+massRangeForCounting);
303 for(Int_t iMB=minBinSum; iMB<=maxBinSum; iMB++){
304 Float_t bkg1=fB1 ? fB1->Eval(hmass[massBin]->GetBinCenter(iMB))/rebinItem : 0;
305 Float_t bkg2=fB2 ? fB2->Eval(hmass[massBin]->GetBinCenter(iMB))/rebinItem : 0;
306 cntSig1+=(hmass[massBin]->GetBinContent(iMB)-bkg1);
307 cntSig2+=(hmass[massBin]->GetBinContent(iMB)-bkg2);
308 cntErr+=(hmass[massBin]->GetBinContent(iMB));
309 }
310 hCntSig1[j]->SetBinContent(iBin+1,cntSig1);
311 hCntSig1[j]->SetBinError(iBin+1,TMath::Sqrt(cntErr));
312 hNDiffCntSig1[j]->SetBinContent(iBin+1,(s-cntSig1)/s);
313 hNDiffCntSig1[j]->SetBinError(iBin+1,TMath::Sqrt(cntErr)/s);
314 hCntSig2[j]->SetBinContent(iBin+1,cntSig2);
315 hNDiffCntSig2[j]->SetBinContent(iBin+1,(s-cntSig2)/s);
316 hNDiffCntSig2[j]->SetBinError(iBin+1,TMath::Sqrt(cntErr)/s);
317 hCntSig2[j]->SetBinError(iBin+1,TMath::Sqrt(cntErr));
318 hSignal[j]->SetBinContent(iBin+1,ry);
319 hSignal[j]->SetBinError(iBin+1,ery);
320 hRelErrSig[j]->SetBinContent(iBin+1,errs/s);
321 hInvSignif[j]->SetBinContent(iBin+1,1/sig);
322 hInvSignif[j]->SetBinError(iBin+1,errsig/(sig*sig));
323 hBackground[j]->SetBinContent(iBin+1,b); //consider sigma
324 hBackground[j]->SetBinError(iBin+1,errb);
325 hBackgroundNormSigma[j]->SetBinContent(iBin+1,b/(3*fitter[massBin]->GetSigma())*(3*0.012)); //consider sigma
326 hBackgroundNormSigma[j]->SetBinError(iBin+1,errb);
327 hSignificance[j]->SetBinContent(iBin+1,sig);
328 hSignificance[j]->SetBinError(iBin+1,errsig);
329 hMass[j]->SetBinContent(iBin+1,mass);
330 hMass[j]->SetBinError(iBin+1,massUnc);
331 hSigma[j]->SetBinContent(iBin+1,sigma);
332 hSigma[j]->SetBinError(iBin+1,sigmaUnc);
333
334 massBin++;
335 delete hRebinned;
336 }// end loop on pt bins
337
338 canvas[j]->Update();
339 canvas[j]->SaveAs(Form("hMass%s_%d_%d.eps",CutsType,typeb,j));
340 // canvas[j]->SaveAs(Form("hMass%s_%d_%d_MultInt.eps",CutsType,typeb,j));
341
342 }// end loop on multiplicity bins
343
344
345 TCanvas *cpar=new TCanvas("cpar","Fit params",1200,600);
346 cpar->Divide(2,1);
347 cpar->cd(1);
348 for(Int_t imult=0; imult<nMultbins; imult++) {
349 hMass[imult]->SetMarkerStyle(20);
350 hMass[imult]->GetXaxis()->SetTitle("Pt (GeV/c)");
351 hMass[imult]->GetYaxis()->SetTitle("Mass (GeV/c^{2})");
352 hMass[imult]->SetMarkerColor(2*imult);
353 if(imult==5) hMass[imult]->SetMarkerColor(2*imult-3);
354 if(imult==0) {
355 hMass[imult]->SetMarkerColor(kBlack);
356 hMass[imult]->Draw("PE");
357 }
358 else hMass[imult]->Draw("PEsame");
359 }
360 cpar->cd(2);
361 for(Int_t imult=0; imult<nMultbins; imult++) {
362 hSigma[imult]->SetMarkerStyle(20);
363 // hSigma[0]->Draw("PE");
364 hSigma[imult]->GetXaxis()->SetTitle("Pt (GeV/c)");
365 hSigma[imult]->GetYaxis()->SetTitle("Sigma (GeV/c^{2})");
366 hSigma[imult]->SetMarkerColor(2*imult);
367 if(imult==5) hSigma[imult]->SetMarkerColor(2*imult-3);
368 if(imult==0) {
369 hSigma[imult]->SetMarkerColor(kBlack);
370 hSigma[imult]->Draw("PE");
371 }
372 else hSigma[imult]->Draw("PEsame");
373 }
374 cpar->Update();
375
376 /*
377 TCanvas** csig;//= new TCanvas*[nMultbins];
378 TCanvas** cDiffS;//=new TCanvas*[nMultbins];
379 for(Int_t i=0; i<nMultbins; i++){
380 csig[i] =new TCanvas(Form("csig_%d",i),Form("Results, multiplicity bin %d",i),1200,600);
381 csig[i]->Divide(3,1);
382 csig[i]->cd(1);
383 hSignal[i]->SetMarkerStyle(20);
384 hSignal[i]->SetMarkerColor(4);
385 hSignal[i]->SetLineColor(4);
386 hSignal[i]->GetXaxis()->SetTitle("Pt (GeV/c)");
387 hSignal[i]->GetYaxis()->SetTitle("Signal");
388 hSignal[i]->Draw("P");
389 hCntSig1[i]->SetMarkerStyle(26);
390 hCntSig1[i]->SetMarkerColor(2);
391 hCntSig1[i]->SetLineColor(2);
392 hCntSig1[i]->Draw("PSAME");
393 hCntSig2[i]->SetMarkerStyle(29);
394 hCntSig2[i]->SetMarkerColor(kGray+1);
395 hCntSig2[i]->SetLineColor(kGray+1);
396 hCntSig2[i]->Draw("PSAME");
397 TLegend* leg=new TLegend(0.4,0.7,0.89,0.89);
398 leg->SetFillColor(0);
399 TLegendEntry* ent=leg->AddEntry(hSignal[i],"From Fit","PL");
400 ent->SetTextColor(hSignal[i]->GetMarkerColor());
401 ent=leg->AddEntry(hCntSig1[i],"From Counting1","PL");
402 ent->SetTextColor(hCntSig1[i]->GetMarkerColor());
403 ent=leg->AddEntry(hCntSig2[i],"From Counting2","PL");
404 ent->SetTextColor(hCntSig2[i]->GetMarkerColor());
405 leg->Draw();
406 csig[i]->cd(2);
407 hBackground[i]->SetMarkerStyle(20);
408 hBackground[i]->Draw("P");
409 hBackground[i]->GetXaxis()->SetTitle("Pt (GeV/c)");
410 hBackground[i]->GetYaxis()->SetTitle("Background");
411 csig[i]->cd(3);
412 hSignificance[i]->SetMarkerStyle(20);
413 hSignificance[i]->Draw("P");
414 hSignificance[i]->GetXaxis()->SetTitle("Pt (GeV/c)");
415 hSignificance[i]->GetYaxis()->SetTitle("Significance");
416
417 cDiffS[i]=new TCanvas(Form("cDiffS_%d",i),Form("Difference, multiplicity bin %d",i),1200,600);
418 cDiffS[i]->Divide(2,1);
419 cDiffS[i]->cd(1);
420 hRelErrSig[i]->SetMarkerStyle(20); //fullcircle
421 hRelErrSig[i]->SetTitleOffset(1.2);
422 hRelErrSig[i]->SetTitle("Rel Error from Fit;p_{t} (GeV/c);Signal Relative Error");
423 hRelErrSig[i]->Draw("P");
424 hInvSignif[i]->SetMarkerStyle(21); //fullsquare
425 hInvSignif[i]->SetMarkerColor(kMagenta+1);
426 hInvSignif[i]->SetLineColor(kMagenta+1);
427 hInvSignif[i]->Draw("PSAMES");
428 TLegend* leg2=new TLegend(0.4,0.7,0.89,0.89);
429 leg2->SetFillColor(0);
430 TLegendEntry* ent2=leg2->AddEntry(hRelErrSig[i],"From Fit","P");
431 ent2->SetTextColor(hRelErrSig[i]->GetMarkerColor());
432 ent2=leg2->AddEntry(hInvSignif[i],"1/Significance","PL");
433 ent2->SetTextColor(hInvSignif[i]->GetMarkerColor());
434 leg2->Draw();
435
436 cDiffS[i]->cd(2);
437 hNDiffCntSig1[i]->SetMarkerStyle(26);
438 hNDiffCntSig1[i]->SetMarkerColor(2);
439 hNDiffCntSig1[i]->SetLineColor(2);
440 hNDiffCntSig1[i]->SetTitle("Cmp Fit-Count;p_{t} (GeV/c);(S_{fit}-S_{count})/S_{fit}");
441 hNDiffCntSig1[i]->Draw("P");
442 hNDiffCntSig2[i]->SetMarkerStyle(29);
443 hNDiffCntSig2[i]->SetMarkerColor(kGray+1);
444 hNDiffCntSig2[i]->SetLineColor(kGray+1);
445 hNDiffCntSig2[i]->Draw("PSAME");
446 TLegend* leg1=new TLegend(0.4,0.7,0.89,0.89);
447 leg1->SetFillColor(0);
448 TLegendEntry* ent1=leg1->AddEntry(hNDiffCntSig1[i],"From Counting1","PL");
449 ent1->SetTextColor(hNDiffCntSig1[i]->GetMarkerColor());
450 ent1=leg1->AddEntry(hNDiffCntSig2[i],"From Counting2","PL");
451 ent1->SetTextColor(hNDiffCntSig2[i]->GetMarkerColor());
452 leg1->Draw();
453 }
454*/
455
456 TGraph* grReducedChiSquare0=new TGraph(nPtBins,ptlims,arrchisquare0);
457 grReducedChiSquare0->SetName("grReducedChiSquare0");
458 grReducedChiSquare0->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");
459 TGraph* grReducedChiSquare1=new TGraph(nPtBins,ptlims,arrchisquare1);
460 grReducedChiSquare1->SetName("grReducedChiSquare1");
461 grReducedChiSquare1->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");
462 TGraph* grReducedChiSquare2=new TGraph(nPtBins,ptlims,arrchisquare2);
463 grReducedChiSquare2->SetName("grReducedChiSquare2");
464 grReducedChiSquare2->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");
465 TGraph* grReducedChiSquare3=new TGraph(nPtBins,ptlims,arrchisquare3);
466 grReducedChiSquare3->SetName("grReducedChiSquare3");
467 grReducedChiSquare3->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");
468 TCanvas *cChi2=new TCanvas("cChi2","reduced chi square",600,600);
469 cChi2->cd();
470 grReducedChiSquare0->SetMarkerStyle(21);
471 grReducedChiSquare0->Draw("AP");
472 grReducedChiSquare1->SetMarkerStyle(22);
473 grReducedChiSquare1->Draw("Psame");
474 grReducedChiSquare2->SetMarkerStyle(23);
475 grReducedChiSquare2->Draw("Psame");
476 grReducedChiSquare3->SetMarkerStyle(24);
477 grReducedChiSquare3->Draw("Psame");
478
479 TCanvas* cbkgNormSigma=new TCanvas("cbkgNormSigma","Background normalized to sigma",400,600);
480 cbkgNormSigma->cd();
481 for(Int_t i=0; i<nMultbins; i++){
482 hBackgroundNormSigma[i]->SetMarkerStyle(20);
483 hBackgroundNormSigma[i]->GetXaxis()->SetTitle("Pt (GeV/c)");
484 hBackgroundNormSigma[i]->GetYaxis()->SetTitle("Background #times 3 #times 0.012/ (3 #times #sigma)");
485 hBackgroundNormSigma[i]->SetMarkerColor(2*i);
486 if(i==5) hBackgroundNormSigma[i]->SetMarkerColor(2*i-3);
487 if(i==0) {
488 hBackgroundNormSigma[i]->SetMarkerColor(kBlack);
489 hBackgroundNormSigma[i]->Draw("PE");
490 }
491 else hBackgroundNormSigma[i]->Draw("Psame");
492 }
493
494
495 TString partname="Both";
496 if(optPartAntiPart==kParticleOnly) {
497 if(analysisType==kD0toKpi) partname="D0";
498 if(analysisType==kDplusKpipi) partname="Dplus";
499 }
500 if(optPartAntiPart==kAntiParticleOnly) {
501 if(analysisType==kD0toKpi) partname="D0bar";
502 if(analysisType==kDplusKpipi) partname="Dminus";
503 }
504
505 TString outfilename = Form("RawYield_Mult_%s_%s",partname.Data(),CutsType);
506 // outfilename += "_MultInt";
507 if(fixPeakSigma) outfilename += "_SigmaFixed";
8f1140d0 508 outfilename += Form("_BCin%1.1fSigma",nSigmaRangeForCounting);
0987c810 509 if(typeb==0) outfilename += "_Expo.root";
510 else if(typeb==1) outfilename += "_Linear.root";
511 else if(typeb==2) outfilename += "_Pol2.root";
512
513 TFile* outf=new TFile(outfilename,"recreate");
514 outf->cd();
515 hNormalization->Write();
516 for(Int_t j=0; j<massBin; j++) hmass[j]->Write();
517 for(Int_t j=0; j<nMultbins; j++){
518 hMass[j]->Write();
519 hSigma[j]->Write();
520 hCntSig1[j]->Write();
521 hCntSig2[j]->Write();
522 hNDiffCntSig1[j]->Write();
523 hNDiffCntSig2[j]->Write();
524 hSignal[j]->Write();
525 hRelErrSig[j]->Write();
526 hInvSignif[j]->Write();
527 hBackground[j]->Write();
528 hBackgroundNormSigma[j]->Write();
529 hSignificance[j]->Write();
530 }
531 grReducedChiSquare0->Write();
532 grReducedChiSquare1->Write();
533 grReducedChiSquare2->Write();
534 grReducedChiSquare3->Write();
535 // hPtMass->Write();
536 outf->Close();
537
538}
539
540//_____________________________________________________________________________________________
541Bool_t LoadD0toKpiHistos(TObjArray* listFiles, TH3F** hPtMassMult, TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, AliNormalizationCounter *counter, const char *CutsType, Int_t Option)
542{
543 Int_t nFiles=listFiles->GetEntries();
544 TList **hlist=new TList*[nFiles];
545 TList **hlistNorm=new TList*[nFiles];
546 AliRDHFCutsD0toKpi** dcuts=new AliRDHFCutsD0toKpi*[nFiles];
547
548 Int_t nReadFiles=0;
549 for(Int_t iFile=0; iFile<nFiles; iFile++){
550 TString fName=((TObjString*)listFiles->At(iFile))->GetString();
551 TFile *f=TFile::Open(fName.Data());
552 if(!f){
553 printf("ERROR: file %s does not exist\n",fName.Data());
554 continue;
555 }
556 printf("Open File %s\n",f->GetName());
557
558 TString dirname="PWG3_D2H_DMult_D0";
559 if(optPartAntiPart==kParticleOnly) dirname+="D0";
560 else if(optPartAntiPart==kAntiParticleOnly) dirname+="D0bar";
561 dirname += CutsType;
562 TDirectory *dir = (TDirectory*)f->Get(dirname);
563 if(!dir){
564 printf("ERROR: directory %s not found in %s\n",dirname.Data(),fName.Data());
565 continue;
566 }
567
568 TString listmassname="coutputD0";
569 if(optPartAntiPart==kParticleOnly) listmassname+="D0";
570 else if(optPartAntiPart==kAntiParticleOnly) listmassname+="D0bar";
571 listmassname += CutsType;
572 printf("List mass name %s\n",listmassname.Data());
573 hlist[nReadFiles]=(TList*)dir->Get(listmassname);
574
575 TString listnorm="coutputNormD0";
576 if(optPartAntiPart==kParticleOnly) listnorm+="D0";
577 else if(optPartAntiPart==kAntiParticleOnly) listnorm+="D0bar";
578 listnorm += CutsType;
579 printf("List norm name %s\n",listnorm.Data());
580 hlistNorm[nReadFiles]=(TList*)dir->Get(listnorm);
581 // AliNormalizationCounter *tmpcounter = (AliNormalizationCounter*)hlistNorm[nReadFiles]->FindObject("NormCounterCorrMult");
582 // counter->Add(tmpcounter);
583 // delete tmpcounter;
584 // counter = tmpcounter;
585
586 TString cutsobjname="coutputCutsD0";
587 if(optPartAntiPart==kParticleOnly) cutsobjname+="D0";
588 else if(optPartAntiPart==kAntiParticleOnly) cutsobjname+="D0bar";
589 cutsobjname += CutsType;
590 printf("Cuts name %s\n",cutsobjname.Data());
591 dcuts[nReadFiles]=(AliRDHFCutsD0toKpi*)dir->Get(cutsobjname);
592 if(!dcuts[nReadFiles]) {
593 printf("ERROR: Cut objects do not match\n");
594 return kFALSE;
595 }
596 /*
597 if(nReadFiles>0){
598 Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);
599 if(!sameCuts){
600 printf("ERROR: Cut objects do not match\n");
601 return kFALSE;
602 }
603 }
604 */
605 nReadFiles++;
606 }
607 if(nReadFiles<nFiles){
608 printf("WARNING: not all requested files have been found\n");
609 if (nReadFiles==0) {
610 printf("ERROR: Any file/dir found\n");
611 return kFALSE;
612 }
613 }
614 // printf("Cuts type %s, particle/antipart %d\n",CutsType,optPartAntiPart);
615
616 /*
617 Int_t nPtBinsCuts=dcuts[0]->GetNPtBins();
618 printf("Number of pt bins for cut object = %d\n",nPtBins);
619 Float_t *ptlimsCuts=dcuts[0]->GetPtBinLimits();
620 ptlimsCuts[nPtBinsCuts]=ptlimsCuts[nPtBinsCuts-1]+4.;
621 */
622
623 printf("Get the 3D histogram \n");
624 const char *histoname="";
625 if(optPartAntiPart==kParticleOnly) histoname= "hPtVsMassvsMultPart";
626 else if(optPartAntiPart==kAntiParticleOnly) histoname="hPtVsMassvsMultAntiPart";
627 else if(optPartAntiPart==kBoth) histoname="hPtVsMassvsMult";
628 if(Option==kUnCorr) histoname="hPtVsMassvsMultUncorr";
629 if(Option==kNoPid) histoname="hPtVsMassvsMultNoPid";
630
631 TH3F * htemp;
632 for(Int_t iFile=0; iFile<nReadFiles; iFile++){
633 printf(" Looking for histo histMass %s for file %d\n",histoname,iFile);
634 htemp=(TH3F*)hlist[iFile]->FindObject(Form("%s",histoname));
635 // cout << " htemp :"<<htemp<<endl;
636 if(!hPtMassMult[0]){
637 hPtMassMult[0]=new TH3F(*htemp);
638 }else{
639 hPtMassMult[0]->Add(htemp);
640 }
641 hNtrZvtx[iFile] = (TH2F*)hlist[iFile]->FindObject("hNtrVsZvtx");
642 hNtrZvtxCorr[iFile] = (TH2F*)hlist[iFile]->FindObject("hNtrCorrVsZvtx");
643 }
644
645 // cout<<" hPtMassMult:"<<hPtMassMult[0]<<endl;
646
647 TString partname="Both";
648 if(optPartAntiPart==kParticleOnly) partname="D0";
649 if(optPartAntiPart==kAntiParticleOnly) partname="D0bar";
650
651 TString outfilename = Form("RawYield%s_%s",partname.Data(),CutsType);
652 if(fixPeakSigma) outfilename += "_SigmaFixed";
653 if(typeb==0) outfilename += "_Expo.root";
654 else if(typeb==1) outfilename += "_Linear.root";
655 else if(typeb==2) outfilename += "_Pol2.root";
656 TFile* outf=new TFile(outfilename,"recreate");
657 outf->cd();
658 dcuts[0]->Write();
659 outf->Close();
660 return kTRUE;
661
662}
663
664//_____________________________________________________________________________________________
665Bool_t LoadDplusHistos(TObjArray* listFiles, TH3F** hPtMassMult, TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, const char *CutsType, Int_t Option)
666{
667Int_t nFiles=listFiles->GetEntries();
668 TList **hlist=new TList*[nFiles];
669 TList **hlistNorm=new TList*[nFiles];
670 AliRDHFCutsDplustoKpipi** dcuts=new AliRDHFCutsDplustoKpipi*[nFiles];
671
672 Int_t nReadFiles=0;
673 for(Int_t iFile=0; iFile<nFiles; iFile++){
674 TString fName=((TObjString*)listFiles->At(iFile))->GetString();
675 TFile *f=TFile::Open(fName.Data());
676 if(!f){
677 printf("ERROR: file %s does not exist\n",fName.Data());
678 continue;
679 }
680 printf("Open File %s\n",f->GetName());
681 TDirectory *dir = (TDirectory*)f->Get(Form("PWG3_D2H_DMult_Dplus%s",suffix.Data()));
682 // TDirectory *dir = (TDirectory*)f->Get("PWG3_D2H_DMult");
683 if(!dir){
684 printf("ERROR: directory PWG3_D2H_DMult not found in %s\n",fName.Data());
685 continue;
686 }
687 hlist[nReadFiles]=(TList*)dir->Get(Form("coutputDplus%s",suffix.Data()));
688 TList *listcut = (TList*)dir->Get(Form("coutputCutsDplus%s",suffix.Data()));
689 TList *listNorm = (TList*)dir->Get(Form("coutputNormDplus%s",suffix.Data()));
690 dcuts[nReadFiles]=(AliRDHFCutsDplustoKpipi*)listcut->At(0);
691 if(nReadFiles>0){
692 Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);
693 if(!sameCuts){
694
695 printf("ERROR: Cut objects do not match\n");
696 return kFALSE;
697 }
698 }
699
700
701
702 nReadFiles++;
703 }
704 if(nReadFiles<nFiles){
705 printf("WARNING: not all requested files have been found\n");
706 if (nReadFiles==0) {
707 printf("ERROR: Any file/dir found\n");
708 return kFALSE;
709 }
710 }
711
712
713 printf("Get the 3D histogram \n");
714 const char *histoname="";
715 if(optPartAntiPart==kParticleOnly) histoname= "hPtVsMassvsMultPart";
716 else if(optPartAntiPart==kAntiParticleOnly) histoname="hPtVsMassvsMultAntiPart";
717 else if(optPartAntiPart==kBoth) histoname="hPtVsMassvsMult";
718 if(Option==kUnCorr) histoname="hPtVsMassvsMultUncorr";
719 if(Option==kNoPid) histoname="hPtVsMassvsMultNoPid";
720
721 TH3F * htemp;
722 for(Int_t iFile=0; iFile<nReadFiles; iFile++){
723 printf(" Looking for histo histMass %s for file %d\n",histoname,iFile);
724 htemp=(TH3F*)hlist[iFile]->FindObject(Form("%s",histoname));
725 // cout << " htemp :"<<htemp<<endl;
726 if(!hPtMassMult[0]){
727 hPtMassMult[0]=new TH3F(*htemp);
728 }else{
729 hPtMassMult[0]->Add(htemp);
730 }
731 // (TH2F*)hNtrZvtx[iFile] = (TH2F*)hlist[iFile]->FindObject("hNtrVsZvtx");
732 //(TH2F*)hNtrZvtxCorr[iFile] = (TH2F*)hlist[iFile]->FindObject("hNtrVsZvtxCorr");
733 }
734
735 // cout<<" hPtMassMult:"<<hPtMassMult[0]<<endl;
736
737 TString partname="Both";
738 if(optPartAntiPart==kParticleOnly) partname="Dplus";
739 if(optPartAntiPart==kAntiParticleOnly) partname="Dminus";
740
741 TString outfilename = Form("RawYield%s_%s",partname.Data(),CutsType);
742 if(fixPeakSigma) outfilename += "_SigmaFixed";
743 if(typeb==0) outfilename += "_Expo.root";
744 else if(typeb==1) outfilename += "_Linear.root";
745 else if(typeb==2) outfilename += "_Pol2.root";
746 TFile* outf=new TFile(outfilename,"recreate");
747 outf->cd();
748 dcuts[0]->Write();
749 outf->Close();
750 return kTRUE;
751
752}
753
754//_____________________________________________________________________________________________
755TH1F* RebinHisto(TH1F* hOrig, Int_t reb, Int_t firstUse){
756 // Rebin histogram, from bin firstUse to lastUse
757 // Use all bins if firstUse=-1
758
759 Int_t nBinOrig=hOrig->GetNbinsX();
760 Int_t firstBinOrig=1;
761 Int_t lastBinOrig=nBinOrig;
762 Int_t nBinOrigUsed=nBinOrig;
763 Int_t nBinFinal=nBinOrig/reb;
764 if(firstUse>=1){
765 firstBinOrig=firstUse;
766 nBinFinal=(nBinOrig-firstUse+1)/reb;
767 nBinOrigUsed=nBinFinal*reb;
768 lastBinOrig=firstBinOrig+nBinOrigUsed-1;
769 }else{
770 Int_t exc=nBinOrigUsed%reb;
771 if(exc!=0){
772 nBinOrigUsed-=exc;
773 firstBinOrig+=exc/2;
774 lastBinOrig=firstBinOrig+nBinOrigUsed-1;
775 }
776 }
777
778 printf("Rebin from %d bins to %d bins -- Used bins=%d in range %d-%d\n",nBinOrig,nBinFinal,nBinOrigUsed,firstBinOrig,lastBinOrig);
779 Float_t lowLim=hOrig->GetXaxis()->GetBinLowEdge(firstBinOrig);
780 Float_t hiLim=hOrig->GetXaxis()->GetBinUpEdge(lastBinOrig);
781 TH1F* hRebin=new TH1F(Form("%s-rebin",hOrig->GetName()),hOrig->GetTitle(),nBinFinal,lowLim,hiLim);
782 Int_t lastSummed=firstBinOrig-1;
783 for(Int_t iBin=1;iBin<=nBinFinal; iBin++){
784 Float_t sum=0.;
785 for(Int_t iOrigBin=0;iOrigBin<reb;iOrigBin++){
786 sum+=hOrig->GetBinContent(lastSummed+1);
787 lastSummed++;
788 }
789 hRebin->SetBinContent(iBin,sum);
790 }
791 return hRebin;
792}
793
794//_____________________________________________________________________________________________
795Bool_t CheckNtrVsZvtx(TH2F** hNtrackVsVtxZ, TH2F** hNtrackVsVtxZCorr, Int_t nFiles)
796{
797
798 TCanvas *cNtrVsZvtx = new TCanvas("cNtrVsZvtx","Ntr Vs Zvtx");
799 cNtrVsZvtx->Divide(2,2);
800 for(Int_t i=0; i<nFiles; i++){
801 cNtrVsZvtx->cd(i+1);
802 // hNtrackVsVtxZ[i]->Fit("pol4");
803 hNtrackVsVtxZ[i]->Draw("colz");
804 cNtrVsZvtx->Update();
805 }
806
807 TCanvas *cNtrVsZvtxCorr = new TCanvas("cNtrVsZvtxCorr","Ntr Vs Zvtx Corr");
808 cNtrVsZvtxCorr->Divide(2,2);
809 for(Int_t i=0; i<nFiles; i++){
810 cNtrVsZvtxCorr->cd(i+1);
811 // hNtrackVsVtxZCorr[i]->Fit("pol4");
812 hNtrackVsVtxZCorr[i]->Draw("colz");
813 }
814
815 TH1F *hNtrAxis = (TH1F*)hNtrackVsVtxZ[0]->ProjectionY("hNtrAxis");
816 TH1F *hZvtx[nMultbins];
817 Int_t firstbin=0, lastbin=0;
818 TCanvas *cZvtx = new TCanvas("cZvtx","Zvtx projections");
819 cZvtx->Divide(2,2);
820 for(Int_t i=0; i<nFiles; i++){
821 cZvtx->cd(i+1);
822 firstbin = hNtrAxis->FindBin( multlims[i] );
823 lastbin = hNtrAxis->FindBin( multlims[i+1] ) -1;
824 hZvtx[i] = (TH1F*)hNtrackVsVtxZ[i]->ProjectionX(Form("hZvtx_%d",i),firstbin,lastbin);
825 hZvtx[i]->Draw();
826 }
827 TH1F *hZvtxCorr[nMultbins];
828 TCanvas *cZvtxCorr = new TCanvas("cZvtxCorr","Zvtx projections Corr");
829 cZvtxCorr->Divide(2,2);
830 for(Int_t i=0; i<nFiles; i++){
831 cZvtxCorr->cd(i+1);
832 firstbin = hNtrAxis->FindBin( multlims[i] );
833 lastbin = hNtrAxis->FindBin( multlims[i+1] ) -1;
834 hZvtxCorr[i] = (TH1F*)hNtrackVsVtxZCorr[i]->ProjectionX(Form("hZvtxCorr_%d",i),firstbin,lastbin);
835 hZvtxCorr[i]->Draw();
836 }
837
838 return kTRUE;
839}