]>
Commit | Line | Data |
---|---|---|
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 | ||
34 | enum {kD0toKpi, kDplusKpipi}; | |
35 | enum {kBoth, kParticleOnly, kAntiParticleOnly}; | |
36 | enum {kExpo=0, kLinear, kPol2}; | |
37 | enum {kGaus=0, kDoubleGaus}; | |
38 | enum {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 }; | |
46 | const Int_t nPtBins=5; | |
47 | Double_t ptlims[nPtBins+1]={1.,2.,4.,8.,12.,24.}; | |
48 | Int_t rebin[nPtBins]={4,4,6,8,8}; | |
49 | Double_t sigmapt[nPtBins]={ 0.008, 0.014, 0.019, 0.027, 0.033 }; | |
50 | Bool_t fixPeakSigma = kTRUE; | |
51 | // | |
52 | const Int_t nMultbins=7; | |
53 | Double_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 | // | |
57 | Int_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 | // | |
59 | Bool_t isMC=kFALSE; | |
60 | Int_t typeb=kExpo; | |
61 | Int_t types=kGaus; | |
62 | Int_t optPartAntiPart=kBoth; | |
63 | Int_t factor4refl=0; | |
64 | Float_t massRangeForCounting=0.05; // GeV --> it is 3 sigmapt[binpt] | |
8f1140d0 | 65 | Float_t nSigmaRangeForCounting=3.0; // 3 sigmapt[binpt] |
0987c810 | 66 | TH2F* hPtMass=0x0; |
67 | TString suffix="StdPid"; | |
68 | ||
69 | ||
70 | // Functions | |
71 | Bool_t LoadDplusHistos(TObjArray* listFiles, TH3F** hPtMassMult, TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, const char *CutsType, Int_t Option); | |
72 | Bool_t LoadD0toKpiHistos(TObjArray* listFiles, TH3F** hPtMassMult, TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, AliNormalizationCounter *counter, const char *CutsType, Int_t Option); | |
73 | Bool_t CheckNtrVsZvtx(TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, Int_t nFiles); | |
74 | TH1F* RebinHisto(TH1F* hOrig, Int_t reb, Int_t firstUse=-1); | |
75 | ||
76 | ||
77 | void 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 | //_____________________________________________________________________________________________ | |
541 | Bool_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 | //_____________________________________________________________________________________________ | |
665 | Bool_t LoadDplusHistos(TObjArray* listFiles, TH3F** hPtMassMult, TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, const char *CutsType, Int_t Option) | |
666 | { | |
667 | Int_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 | //_____________________________________________________________________________________________ | |
755 | TH1F* 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 | //_____________________________________________________________________________________________ | |
795 | Bool_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 | } |