]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/vertexingHF/macros/FitMassSpectra.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
temporary hijack for the systematic uncertainties of pPb vs centrality (assigned...
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / macros / FitMassSpectra.C
CommitLineData
b92c0195 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 <TH1F.h>
0587ab27 10#include <TH2F.h>
11#include <TH1D.h>
b92c0195 12#include <TF1.h>
cb2336fe 13#include <TStyle.h>
14#include <TLegend.h>
15#include <TLegendEntry.h>
b92c0195 16#include <TDatabasePDG.h>
5f1ad8a6 17#include <TGraph.h>
b92c0195 18
19#include "AliAODRecoDecayHF.h"
20#include "AliRDHFCuts.h"
21#include "AliRDHFCutsDplustoKpipi.h"
470a7b55 22#include "AliRDHFCutsDStartoKpipi.h"
23851d0b 23#include "AliRDHFCutsDstoKKpi.h"
cb2336fe 24#include "AliRDHFCutsD0toKpi.h"
b92c0195 25#include "AliHFMassFitter.h"
642de1ee 26#include "AliNormalizationCounter.h"
b92c0195 27#endif
28
29
30// MACRO to perform fits to D meson invariant mass spectra
31// and store raw yields and cut object into a root output file
32// Origin: F. Prino (prino@to.infn.it)
cb2336fe 33// D0: C. Bianchin (cbianchi@pd.infn.it)
b92c0195 34
35
36
37//
470a7b55 38enum {kD0toKpi, kDplusKpipi, kDStarD0pi, kDsKKpi};
cb2336fe 39enum {kBoth, kParticleOnly, kAntiParticleOnly};
470a7b55 40enum {kExpo=0, kThrExpo=5, kLinear, kPol2};
b92c0195 41enum {kGaus=0, kDoubleGaus};
42
43
44// Common variables: to be configured by the user
0587ab27 45const Int_t nPtBins=6;
46Double_t ptlims[nPtBins+1]={2.,3.,4.,5.,6.,8.,12.};
1918f6ea 47Int_t rebin[nPtBins]={2,4,4,4,4,4};
48Int_t firstUsedBin[nPtBins]={-1,-1,-1,-1,-1,-1}; // -1 uses all bins, >=1 to set the lower bin to be accepted from original histo
36f9ecb7 49
642de1ee 50TString suffix="Loose_SecondSet1236_ForCF08";
51
52
0587ab27 53//const Int_t nPtBins=7;//6;
54//Double_t ptlims[nPtBins+1]={1.,2.,3.,4.,5.,6.,8.,12.};
5f1ad8a6 55//Int_t rebin[nPtBins+1]={8,6,10,10,10,10,10,10}; //for looser cuts
0587ab27 56//Int_t rebin[nPtBins+1]={10,10,10,10,10,10,10,10}; //for Chiara's cuts
5f1ad8a6 57Int_t typeb=kExpo;
b92c0195 58Int_t types=kGaus;
cb2336fe 59Int_t optPartAntiPart=kBoth;
b92c0195 60Int_t factor4refl=0;
e654a7da 61Double_t minMassForFit=1.7;
62Double_t maxMassForFit=2.1;
b92c0195 63Float_t massRangeForCounting=0.05; // GeV
0587ab27 64TH2F* hPtMass=0x0;
642de1ee 65Double_t nEventsForNorm=0.;
b92c0195 66
cb2336fe 67//for D0only
36f9ecb7 68Bool_t cutsappliedondistr=kTRUE;
5f1ad8a6 69const Int_t nsamples=2;//3;
70//Int_t nevents[nsamples]={8.5635859e+07 /*LHC10b+cpart*/,8.9700624e+07/*LHC10dpart*/};
71//Int_t nevents[nsamples]={9.0374946e+07 /*LHC10b+c*/,9.7593785e+07/*LHC10d*/};
72//Int_t nevents[nsamples]={1.1777812e+08 /*LHC10dnewTPCpid*/,9.7593785e+07/*LHC10d*/};
73//Int_t nevents[nsamples]={1.1777812e+08 /*LHC10dnewTPCpid*/,9.7593785e+07/*LHC10d*/,9.0374946e+07 /*LHC10b+c*/};
74Int_t nevents[nsamples]={1.18860695e+08 /*LHC10dnewTPCpid*/,9.0374946e+07 /*LHC10b+c*/};
75//Int_t nevents[nsamples]={1. /*LHC10dnewTPCpid*/,1 /*LHC10dnewTPCpidrescale*/};
b92c0195 76// Functions
77
78Bool_t LoadDplusHistos(TObjArray* listFiles, TH1F** hMass);
23851d0b 79Bool_t LoadDsHistos(TObjArray* listFiles, TH1F** hMass);
b92c0195 80Bool_t LoadD0toKpiHistos(TObjArray* listFiles, TH1F** hMass);
470a7b55 81Bool_t LoadDstarD0piHistos(TObjArray* listFiles, TH1F** hMass);
82
1918f6ea 83TH1F* RebinHisto(TH1F* hOrig, Int_t reb, Int_t firstUse=-1);
b92c0195 84
85
86
87
88void FitMassSpectra(Int_t analysisType=kDplusKpipi,
23851d0b 89 TString fileNameb="",
90 TString fileNamec="",
91 TString fileNamed="",
92 TString fileNamee=""
b92c0195 93 ){
94 //
95
1918f6ea 96 gInterpreter->ExecuteMacro("$ALICE_ROOT/PWGHF/vertexingHF/macros/LoadLibraries.C");
b92c0195 97 gStyle->SetOptTitle(0);
98
99 TObjArray* listFiles=new TObjArray();
100 if(fileNameb!="") listFiles->AddLast(new TObjString(fileNameb.Data()));
101 if(fileNamec!="") listFiles->AddLast(new TObjString(fileNamec.Data()));
102 if(fileNamed!="") listFiles->AddLast(new TObjString(fileNamed.Data()));
23851d0b 103 if(fileNamee!="") listFiles->AddLast(new TObjString(fileNamee.Data()));
b92c0195 104 if(listFiles->GetEntries()==0){
105 printf("Missing file names in input\n");
106 return;
107 }
108
b92c0195 109 TH1F** hmass=new TH1F*[nPtBins];
110 for(Int_t i=0;i<nPtBins;i++) hmass[i]=0x0;
cb2336fe 111
470a7b55 112 Float_t massD, massD0_fDstar;
b92c0195 113 Bool_t retCode;
114 if(analysisType==kD0toKpi){
115 retCode=LoadD0toKpiHistos(listFiles,hmass);
116 massD=TDatabasePDG::Instance()->GetParticle(421)->Mass();
117 }
118 else if(analysisType==kDplusKpipi){
119 retCode=LoadDplusHistos(listFiles,hmass);
120 massD=TDatabasePDG::Instance()->GetParticle(411)->Mass();
121 }
470a7b55 122 else if(analysisType==kDStarD0pi){
123 retCode=LoadDstarD0piHistos(listFiles,hmass);
124 massD=TDatabasePDG::Instance()->GetParticle(413)->Mass();
125 massD0_fDstar=TDatabasePDG::Instance()->GetParticle(421)->Mass();
126 massD =massD-massD0_fDstar;
127 }
23851d0b 128 else if(analysisType==kDsKKpi){
129 retCode=LoadDsHistos(listFiles,hmass);
130 massD=TDatabasePDG::Instance()->GetParticle(431)->Mass();
131 }
b92c0195 132 else{
23851d0b 133 printf("Wrong analysis type parameter\n");
b92c0195 134 return;
135 }
136 if(!retCode){
137 printf("ERROR in reading input files\n");
138 return;
139 }
140
141
142
0587ab27 143 TH1D* hCntSig1=new TH1D("hCntSig1","hCntSig1",nPtBins,ptlims);
144 TH1D* hCntSig2=new TH1D("hCntSig2","hCntSig2",nPtBins,ptlims);
145 TH1D* hNDiffCntSig1=new TH1D("hNDiffCntSig1","hNDiffCntSig1",nPtBins,ptlims);
146 TH1D* hNDiffCntSig2=new TH1D("hNDiffCntSig2","hNDiffCntSig2",nPtBins,ptlims);
147 TH1D* hSignal=new TH1D("hSignal","hSignal",nPtBins,ptlims);
148 TH1D* hRelErrSig=new TH1D("hRelErrSig","hRelErrSig",nPtBins,ptlims);
149 TH1D* hInvSignif=new TH1D("hInvSignif","hInvSignif",nPtBins,ptlims);
150 TH1D* hBackground=new TH1D("hBackground","hBackground",nPtBins,ptlims);
151 TH1D* hBackgroundNormSigma=new TH1D("hBackgroundNormSigma","hBackgroundNormSigma",nPtBins,ptlims);
152 TH1D* hSignificance=new TH1D("hSignificance","hSignificance",nPtBins,ptlims);
153 TH1D* hMass=new TH1D("hMass","hMass",nPtBins,ptlims);
154 TH1D* hSigma=new TH1D("hSigma","hSigma",nPtBins,ptlims);
155
b92c0195 156
1918f6ea 157 Int_t nMassBins=hmass[0]->GetNbinsX();
e654a7da 158 Double_t hmin=TMath::Max(minMassForFit,hmass[0]->GetBinLowEdge(2));
159 Double_t hmax=TMath::Min(maxMassForFit,hmass[0]->GetBinLowEdge(nMassBins-2)+hmass[0]->GetBinWidth(nMassBins-2));
1918f6ea 160 Float_t minBinSum=hmass[0]->FindBin(massD-massRangeForCounting);
161 Float_t maxBinSum=hmass[0]->FindBin(massD+massRangeForCounting);
b92c0195 162 Int_t iPad=1;
163
cb2336fe 164 TF1* funBckStore1=0x0;
165 TF1* funBckStore2=0x0;
166 TF1* funBckStore3=0x0;
b92c0195 167
168 AliHFMassFitter** fitter=new AliHFMassFitter*[nPtBins];
5f1ad8a6 169 Double_t arrchisquare[nPtBins];
b92c0195 170 TCanvas* c1= new TCanvas("c1","MassSpectra",1000,800);
171 Int_t nx=3, ny=2;
172 if(nPtBins>6){
173 nx=4;
174 ny=3;
175 }
176 if(nPtBins>12){
177 nx=5;
178 ny=4;
179 }
180 c1->Divide(nx,ny);
181
182 Double_t sig,errsig,s,errs,b,errb;
183 for(Int_t iBin=0; iBin<nPtBins; iBin++){
184 c1->cd(iPad++);
53d1574e 185 Int_t origNbins=hmass[iBin]->GetNbinsX();
1918f6ea 186 TH1F* hRebinned=RebinHisto(hmass[iBin],rebin[iBin],firstUsedBin[iBin]);
e654a7da 187 hmin=TMath::Max(minMassForFit,hRebinned->GetBinLowEdge(2));
188 hmax=TMath::Min(maxMassForFit,hRebinned->GetBinLowEdge(hRebinned->GetNbinsX()));
1918f6ea 189 fitter[iBin]=new AliHFMassFitter(hRebinned,hmin, hmax,1,typeb,types);
53d1574e 190 rebin[iBin]=origNbins/fitter[iBin]->GetBinN();
b92c0195 191 fitter[iBin]->SetReflectionSigmaFactor(factor4refl);
192 fitter[iBin]->SetInitialGaussianMean(massD);
470a7b55 193 if(analysisType==kDStarD0pi) fitter[iBin]->SetInitialGaussianSigma(0.00065);
b92c0195 194 Bool_t out=fitter[iBin]->MassFitter(0);
5f1ad8a6 195 if(!out) {
196 fitter[iBin]->GetHistoClone()->Draw();
197 continue;
198 }
0587ab27 199 Double_t mass=fitter[iBin]->GetMean();
200 Double_t sigma=fitter[iBin]->GetSigma();
5f1ad8a6 201 arrchisquare[iBin]=fitter[iBin]->GetReducedChiSquare();
b92c0195 202 TF1* fB1=fitter[iBin]->GetBackgroundFullRangeFunc();
203 TF1* fB2=fitter[iBin]->GetBackgroundRecalcFunc();
204 TF1* fM=fitter[iBin]->GetMassFunc();
205 if(iBin==0 && fB1) funBckStore1=new TF1(*fB1);
206 if(iBin==0 && fB2) funBckStore2=new TF1(*fB2);
207 if(iBin==0 && fM) funBckStore3=new TF1(*fM);
208
209 fitter[iBin]->DrawHere(gPad);
210 fitter[iBin]->Signal(3,s,errs);
211 fitter[iBin]->Background(3,b,errb);
212 fitter[iBin]->Significance(3,sig,errsig);
e654a7da 213 Double_t ry=fitter[iBin]->GetRawYield();
214 Double_t ery=fitter[iBin]->GetRawYieldError();
b92c0195 215 Float_t cntSig1=0.;
216 Float_t cntSig2=0.;
217 Float_t cntErr=0.;
218 for(Int_t iMB=minBinSum; iMB<=maxBinSum; iMB++){
cb2336fe 219 Float_t bkg1=fB1 ? fB1->Eval(hmass[iBin]->GetBinCenter(iMB))/rebin[iBin] : 0;
220 Float_t bkg2=fB2 ? fB2->Eval(hmass[iBin]->GetBinCenter(iMB))/rebin[iBin] : 0;
b92c0195 221 cntSig1+=(hmass[iBin]->GetBinContent(iMB)-bkg1);
222 cntSig2+=(hmass[iBin]->GetBinContent(iMB)-bkg2);
223 cntErr+=(hmass[iBin]->GetBinContent(iMB));
224 }
225 hCntSig1->SetBinContent(iBin+1,cntSig1);
226 hCntSig1->SetBinError(iBin+1,TMath::Sqrt(cntErr));
36f9ecb7 227 hNDiffCntSig1->SetBinContent(iBin+1,(s-cntSig1)/s);
228 hNDiffCntSig1->SetBinError(iBin+1,TMath::Sqrt(cntErr)/s);
b92c0195 229 hCntSig2->SetBinContent(iBin+1,cntSig2);
36f9ecb7 230 hNDiffCntSig2->SetBinContent(iBin+1,(s-cntSig2)/s);
231 hNDiffCntSig2->SetBinError(iBin+1,TMath::Sqrt(cntErr)/s);
b92c0195 232 hCntSig2->SetBinError(iBin+1,TMath::Sqrt(cntErr));
e654a7da 233 hSignal->SetBinContent(iBin+1,ry);
234 hSignal->SetBinError(iBin+1,ery);
36f9ecb7 235 hRelErrSig->SetBinContent(iBin+1,errs/s);
236 hInvSignif->SetBinContent(iBin+1,1/sig);
237 hInvSignif->SetBinError(iBin+1,errsig/(sig*sig));
5f1ad8a6 238 hBackground->SetBinContent(iBin+1,b); //consider sigma
b92c0195 239 hBackground->SetBinError(iBin+1,errb);
5f1ad8a6 240 hBackgroundNormSigma->SetBinContent(iBin+1,b/(3*fitter[iBin]->GetSigma())*(3*0.012)); //consider sigma
241 hBackgroundNormSigma->SetBinError(iBin+1,errb);
b92c0195 242 hSignificance->SetBinContent(iBin+1,sig);
243 hSignificance->SetBinError(iBin+1,errsig);
0587ab27 244 hMass->SetBinContent(iBin+1,mass);
245 hMass->SetBinError(iBin+1,0.0001);
246 hSigma->SetBinContent(iBin+1,sigma);
247 hSigma->SetBinError(iBin+1,0.0001);
b92c0195 248
249 }
0587ab27 250
5f1ad8a6 251 /*
b92c0195 252 c1->cd(1); // is some cases the fitting function of 1st bin get lost
253 funBckStore1->Draw("same");
254 funBckStore2->Draw("same");
255 funBckStore3->Draw("same");
5f1ad8a6 256 */
0587ab27 257
258 TCanvas *cpar=new TCanvas("cpar","Fit params",1200,600);
259 cpar->Divide(2,1);
260 cpar->cd(1);
261 hMass->SetMarkerStyle(20);
262 hMass->Draw("PE");
263 hMass->GetXaxis()->SetTitle("Pt (GeV/c)");
264 hMass->GetXaxis()->SetTitle("Mass (GeV/c^{2})");
265 cpar->cd(2);
266 hSigma->SetMarkerStyle(20);
267 hSigma->Draw("PE");
268 hSigma->GetXaxis()->SetTitle("Pt (GeV/c)");
269 hSigma->GetXaxis()->SetTitle("Sigma (GeV/c^{2})");
270
b92c0195 271 TCanvas* csig=new TCanvas("csig","Results",1200,600);
272 csig->Divide(3,1);
273 csig->cd(1);
274 hSignal->SetMarkerStyle(20);
275 hSignal->SetMarkerColor(4);
276 hSignal->SetLineColor(4);
277 hSignal->GetXaxis()->SetTitle("Pt (GeV/c)");
278 hSignal->GetYaxis()->SetTitle("Signal");
279 hSignal->Draw("P");
280 hCntSig1->SetMarkerStyle(26);
281 hCntSig1->SetMarkerColor(2);
282 hCntSig1->SetLineColor(2);
283 hCntSig1->Draw("PSAME");
284 hCntSig2->SetMarkerStyle(29);
285 hCntSig2->SetMarkerColor(kGray+1);
286 hCntSig2->SetLineColor(kGray+1);
287 hCntSig2->Draw("PSAME");
288 TLegend* leg=new TLegend(0.4,0.7,0.89,0.89);
289 leg->SetFillColor(0);
290 TLegendEntry* ent=leg->AddEntry(hSignal,"From Fit","PL");
291 ent->SetTextColor(hSignal->GetMarkerColor());
292 ent=leg->AddEntry(hCntSig1,"From Counting1","PL");
293 ent->SetTextColor(hCntSig1->GetMarkerColor());
294 ent=leg->AddEntry(hCntSig2,"From Counting2","PL");
295 ent->SetTextColor(hCntSig2->GetMarkerColor());
296 leg->Draw();
297 csig->cd(2);
298 hBackground->SetMarkerStyle(20);
299 hBackground->Draw("P");
300 hBackground->GetXaxis()->SetTitle("Pt (GeV/c)");
301 hBackground->GetYaxis()->SetTitle("Background");
302 csig->cd(3);
303 hSignificance->SetMarkerStyle(20);
304 hSignificance->Draw("P");
305 hSignificance->GetXaxis()->SetTitle("Pt (GeV/c)");
306 hSignificance->GetYaxis()->SetTitle("Significance");
307
36f9ecb7 308 TCanvas* cDiffS=new TCanvas("cDiffS","Difference",1200,600);
309 cDiffS->Divide(2,1);
310 cDiffS->cd(1);
311 hRelErrSig->SetMarkerStyle(20); //fullcircle
312 hRelErrSig->SetTitleOffset(1.2);
5f1ad8a6 313 hRelErrSig->SetTitle("Rel Error from Fit;p_{t} (GeV/c);Signal Relative Error");
36f9ecb7 314 hRelErrSig->Draw("P");
315 hInvSignif->SetMarkerStyle(21); //fullsquare
316 hInvSignif->SetMarkerColor(kMagenta+1);
317 hInvSignif->SetLineColor(kMagenta+1);
318 hInvSignif->Draw("PSAMES");
319 TLegend* leg2=new TLegend(0.4,0.7,0.89,0.89);
320 leg2->SetFillColor(0);
321 TLegendEntry* ent2=leg2->AddEntry(hRelErrSig,"From Fit","P");
322 ent2->SetTextColor(hRelErrSig->GetMarkerColor());
323 ent2=leg2->AddEntry(hInvSignif,"1/Significance","PL");
324 ent2->SetTextColor(hInvSignif->GetMarkerColor());
325 leg2->Draw();
326
327 cDiffS->cd(2);
328 hNDiffCntSig1->SetMarkerStyle(26);
329 hNDiffCntSig1->SetMarkerColor(2);
330 hNDiffCntSig1->SetLineColor(2);
5f1ad8a6 331 hNDiffCntSig1->SetTitle("Cmp Fit-Count;p_{t} (GeV/c);(S_{fit}-S_{count})/S_{fit}");
36f9ecb7 332 hNDiffCntSig1->Draw("P");
333 hNDiffCntSig2->SetMarkerStyle(29);
334 hNDiffCntSig2->SetMarkerColor(kGray+1);
335 hNDiffCntSig2->SetLineColor(kGray+1);
336 hNDiffCntSig2->Draw("PSAME");
337 TLegend* leg1=new TLegend(0.4,0.7,0.89,0.89);
338 leg1->SetFillColor(0);
339 TLegendEntry* ent1=leg1->AddEntry(hNDiffCntSig1,"From Counting1","PL");
340 ent1->SetTextColor(hNDiffCntSig1->GetMarkerColor());
341 ent1=leg1->AddEntry(hNDiffCntSig2,"From Counting2","PL");
342 ent1->SetTextColor(hNDiffCntSig2->GetMarkerColor());
343 leg1->Draw();
344
5f1ad8a6 345 TGraph* grReducedChiSquare=new TGraph(nPtBins,ptlims,arrchisquare);
346 grReducedChiSquare->SetName("grReducedChiSquare");
347 grReducedChiSquare->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");
348 TCanvas *cChi2=new TCanvas("cChi2","reduced chi square",600,600);
349 cChi2->cd();
350 grReducedChiSquare->SetMarkerStyle(21);
351 grReducedChiSquare->Draw("AP");
352
353 TCanvas* cbkgNormSigma=new TCanvas("cbkgNormSigma","Background normalized to sigma",400,600);
354 cbkgNormSigma->cd();
355 hBackgroundNormSigma->SetMarkerStyle(20);
356 hBackgroundNormSigma->Draw("P");
357 hBackgroundNormSigma->GetXaxis()->SetTitle("Pt (GeV/c)");
358 hBackgroundNormSigma->GetYaxis()->SetTitle("Background #times 3 #times 0.012/ (3 #times #sigma)");
359 hBackgroundNormSigma->Draw();
360
361
362 TString partname="Both";
363 if(optPartAntiPart==kParticleOnly) {
364 if(analysisType==kD0toKpi) partname="D0";
365 if(analysisType==kDplusKpipi) partname="Dplus";
23851d0b 366 if(analysisType==kDsKKpi) partname="Dsplus";
5f1ad8a6 367 }
368 if(optPartAntiPart==kAntiParticleOnly) {
369 if(analysisType==kD0toKpi) partname="D0bar";
370 if(analysisType==kDplusKpipi) partname="Dminus";
23851d0b 371 if(analysisType==kDsKKpi) partname="Dsminus";
5f1ad8a6 372 }
36f9ecb7 373
642de1ee 374 printf("Events for norm = %f\n",nEventsForNorm);
375 TH1F* hNEvents=new TH1F("hNEvents","",1,0.,1.);
376 hNEvents->SetBinContent(1,nEventsForNorm);
377
5f1ad8a6 378 TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"update");
b92c0195 379 outf->cd();
642de1ee 380 hNEvents->Write();
0587ab27 381 hMass->Write();
382 hSigma->Write();
b92c0195 383 hCntSig1->Write();
384 hCntSig2->Write();
36f9ecb7 385 hNDiffCntSig1->Write();
386 hNDiffCntSig2->Write();
b92c0195 387 hSignal->Write();
36f9ecb7 388 hRelErrSig->Write();
389 hInvSignif->Write();
b92c0195 390 hBackground->Write();
5f1ad8a6 391 hBackgroundNormSigma->Write();
b92c0195 392 hSignificance->Write();
5f1ad8a6 393 grReducedChiSquare->Write();
0587ab27 394 hPtMass->Write();
b92c0195 395 outf->Close();
396}
397
398
399Bool_t LoadDplusHistos(TObjArray* listFiles, TH1F** hMass){
400
401 Int_t nFiles=listFiles->GetEntries();
402 TList **hlist=new TList*[nFiles];
403 AliRDHFCutsDplustoKpipi** dcuts=new AliRDHFCutsDplustoKpipi*[nFiles];
404
405 Int_t nReadFiles=0;
406 for(Int_t iFile=0; iFile<nFiles; iFile++){
407 TString fName=((TObjString*)listFiles->At(iFile))->GetString();
408 TFile *f=TFile::Open(fName.Data());
409 if(!f){
410 printf("ERROR: file %s does not exist\n",fName.Data());
411 continue;
412 }
413 printf("Open File %s\n",f->GetName());
414 TDirectory *dir = (TDirectory*)f->Get("PWG3_D2H_InvMassDplus");
415 if(!dir){
cb2336fe 416 printf("ERROR: directory PWG3_D2H_InvMassDplus not found in %s\n",fName.Data());
b92c0195 417 continue;
418 }
642de1ee 419 hlist[nReadFiles]=(TList*)dir->Get(Form("coutputDplus%s",suffix.Data()));
420 TList *listcut = (TList*)dir->Get(Form("coutputDplusCuts%s",suffix.Data()));
421 dcuts[nReadFiles]=(AliRDHFCutsDplustoKpipi*)listcut->At(0);
b92c0195 422 if(nReadFiles>0){
423 Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);
424 if(!sameCuts){
425 printf("ERROR: Cut objects do not match\n");
426 return kFALSE;
427 }
428 }
642de1ee 429 AliNormalizationCounter* c=(AliNormalizationCounter*)dir->Get(Form("coutputDplusNorm%s",suffix.Data()));
430 printf("Events for normalization = %f\n",c->GetNEventsForNorm());
431 nEventsForNorm+=c->GetNEventsForNorm();
b92c0195 432 nReadFiles++;
433 }
434 if(nReadFiles<nFiles){
435 printf("WARNING: not all requested files have been found\n");
436 }
437
438 Int_t nPtBinsCuts=dcuts[0]->GetNPtBins();
439 printf("Number of pt bins for cut object = %d\n",nPtBins);
440 Float_t *ptlimsCuts=dcuts[0]->GetPtBinLimits();
441 ptlimsCuts[nPtBinsCuts]=ptlimsCuts[nPtBinsCuts-1]+4.;
442
443 Int_t iFinBin=0;
444 for(Int_t i=0;i<nPtBinsCuts;i++){
445 if(ptlimsCuts[i]>=ptlims[iFinBin+1]) iFinBin+=1;
446 if(iFinBin>nPtBins) break;
447 if(ptlimsCuts[i]>=ptlims[iFinBin] &&
448 ptlimsCuts[i+1]<=ptlims[iFinBin+1]){
449 for(Int_t iFile=0; iFile<nReadFiles; iFile++){
cb2336fe 450 TString histoName;
451 if(optPartAntiPart==kBoth) histoName.Form("hMassPt%dTC",i);
452 else if(optPartAntiPart==kParticleOnly) histoName.Form("hMassPt%dTCPlus",i);
453 else if(optPartAntiPart==kAntiParticleOnly) histoName.Form("hMassPt%dTCMinus",i);
454 TH1F* htemp=(TH1F*)hlist[iFile]->FindObject(histoName.Data());
455 if(!htemp){
456 printf("ERROR: Histogram %s not found\n",histoName.Data());
457 return kFALSE;
458 }
b92c0195 459 if(!hMass[iFinBin]){
460 hMass[iFinBin]=new TH1F(*htemp);
461 }else{
462 hMass[iFinBin]->Add(htemp);
463 }
464 }
465 }
466 }
5f1ad8a6 467 TString partname="Both";
468 if(optPartAntiPart==kParticleOnly) partname="Dplus";
469 if(optPartAntiPart==kAntiParticleOnly) partname="Dminus";
b92c0195 470
0587ab27 471 for(Int_t iFile=0; iFile<nReadFiles; iFile++){
472 TH2F* htemp2=(TH2F*)hlist[iFile]->FindObject("hPtVsMassTC");
473 if(iFile==0){
474 hPtMass=new TH2F(*htemp2);
475 }else{
476 hPtMass->Add(htemp2);
477 }
478 }
479
5f1ad8a6 480 TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"recreate");
b92c0195 481 outf->cd();
482 dcuts[0]->Write();
483 outf->Close();
484
485 return kTRUE;
486
487}
488
23851d0b 489
490Bool_t LoadDsHistos(TObjArray* listFiles, TH1F** hMass){
491
492 Int_t nFiles=listFiles->GetEntries();
493 TList **hlist=new TList*[nFiles];
494 AliRDHFCutsDstoKKpi** dcuts=new AliRDHFCutsDstoKKpi*[nFiles];
495
496 Int_t nReadFiles=0;
497 for(Int_t iFile=0; iFile<nFiles; iFile++){
498 TString fName=((TObjString*)listFiles->At(iFile))->GetString();
499 TFile *f=TFile::Open(fName.Data());
500 if(!f){
501 printf("ERROR: file %s does not exist\n",fName.Data());
502 continue;
503 }
504 printf("Open File %s\n",f->GetName());
505 TDirectory *dir = (TDirectory*)f->Get("PWG3_D2H_InvMassDs");
506 if(!dir){
507 printf("ERROR: directory PWG3_D2H_InvMassDs not found in %s\n",fName.Data());
508 continue;
509 }
510 hlist[nReadFiles]=(TList*)dir->Get("coutputDs");
511 TList *listcut = (TList*)dir->Get("coutputDsCuts");
512 dcuts[nReadFiles]=(AliRDHFCutsDstoKKpi*)listcut->At(0);
513 cout<< dcuts[nReadFiles]<<endl;
514 if(nReadFiles>0){
515 Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);
516 if(!sameCuts){
517 printf("ERROR: Cut objects do not match\n");
518 return kFALSE;
519 }
520 }
521 nReadFiles++;
522 }
523 if(nReadFiles<nFiles){
524 printf("WARNING: not all requested files have been found\n");
525 }
526
527 Int_t nPtBinsCuts=dcuts[0]->GetNPtBins();
528 printf("Number of pt bins for cut object = %d\n",nPtBins);
529 Float_t *ptlimsCuts=dcuts[0]->GetPtBinLimits();
530 ptlimsCuts[nPtBinsCuts]=ptlimsCuts[nPtBinsCuts-1]+4.;
531
532 Int_t iFinBin=0;
533 for(Int_t i=0;i<nPtBinsCuts;i++){
534 if(ptlimsCuts[i]>=ptlims[iFinBin+1]) iFinBin+=1;
535 if(iFinBin>nPtBins) break;
536 if(ptlimsCuts[i]>=ptlims[iFinBin] &&
537 ptlimsCuts[i+1]<=ptlims[iFinBin+1]){
538 for(Int_t iFile=0; iFile<nReadFiles; iFile++){
539 TString histoName;
540 if(optPartAntiPart==kBoth) histoName.Form("hMassAllPt%dphi",i);
541 else if(optPartAntiPart==kParticleOnly){
542 printf("Particle/Antiparticle not yet enabled for Ds");
543 histoName.Form("hMassAllPt%dphi",i);
544 }
545 else if(optPartAntiPart==kAntiParticleOnly){
546 printf("Particle/Antiparticle not yet enabled for Ds");
547 histoName.Form("hMassAllPt%dphi",i);
548 }
549 TH1F* htemp=(TH1F*)hlist[iFile]->FindObject(histoName.Data());
550 if(!htemp){
551 printf("ERROR: Histogram %s not found\n",histoName.Data());
552 return kFALSE;
553 }
554 if(!hMass[iFinBin]){
555 hMass[iFinBin]=new TH1F(*htemp);
556 }else{
557 hMass[iFinBin]->Add(htemp);
558 }
559 }
560 }
561 }
562 TString partname="Both";
563 if(optPartAntiPart==kParticleOnly) partname="Both";
564 if(optPartAntiPart==kAntiParticleOnly) partname="Both";
565
566 for(Int_t iFile=0; iFile<nReadFiles; iFile++){
567 TH2F* htemp2=(TH2F*)hlist[iFile]->FindObject("hPtVsMassPhi");
568 if(iFile==0){
569 hPtMass=new TH2F(*htemp2);
570 }else{
571 hPtMass->Add(htemp2);
572 }
573 }
574
575 TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"recreate");
576 outf->cd();
577 dcuts[0]->Write();
578 outf->Close();
579
580 return kTRUE;
581
582}
583
584
b92c0195 585Bool_t LoadD0toKpiHistos(TObjArray* listFiles, TH1F** hMass){
cb2336fe 586 //
587 Int_t nFiles=listFiles->GetEntries();
588 TList **hlist=new TList*[nFiles];
589 AliRDHFCutsD0toKpi** dcuts=new AliRDHFCutsD0toKpi*[nFiles];
590
591 Int_t nReadFiles=0;
592 for(Int_t iFile=0; iFile<nFiles; iFile++){
593 TString fName=((TObjString*)listFiles->At(iFile))->GetString();
594 TFile *f=TFile::Open(fName.Data());
595 if(!f){
596 printf("ERROR: file %s does not exist\n",fName.Data());
597 continue;
598 }
599 printf("Open File %s\n",f->GetName());
600
601 TString dirname="PWG3_D2H_D0InvMass";
602 if(optPartAntiPart==kParticleOnly) dirname+="D0";
603 else if(optPartAntiPart==kAntiParticleOnly) dirname+="D0bar";
604 if(cutsappliedondistr) dirname+="C";
605 TDirectory *dir = (TDirectory*)f->Get(dirname);
606 if(!dir){
607 printf("ERROR: directory %s not found in %s\n",dirname.Data(),fName.Data());
608 continue;
609 }
610 TString listmassname="coutputmassD0Mass";
5f1ad8a6 611 if(optPartAntiPart==kParticleOnly) listmassname+="D0";
612 else if(optPartAntiPart==kAntiParticleOnly) listmassname+="D0bar";
cb2336fe 613 if(cutsappliedondistr) listmassname+="C";
614
615 hlist[nReadFiles]=(TList*)dir->Get(listmassname);
616
617 TString cutsobjname="cutsD0";
5f1ad8a6 618 if(optPartAntiPart==kParticleOnly) cutsobjname+="D0";
619 else if(optPartAntiPart==kAntiParticleOnly) cutsobjname+="D0bar";
cb2336fe 620 if(cutsappliedondistr) cutsobjname+="C";
621
622 dcuts[nReadFiles]=(AliRDHFCutsD0toKpi*)dir->Get(cutsobjname);
623 if(nReadFiles>0){
624 Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);
625 if(!sameCuts){
626 printf("ERROR: Cut objects do not match\n");
627 return kFALSE;
628 }
629 }
630 nReadFiles++;
631 }
632 if(nReadFiles<nFiles){
633 printf("WARNING: not all requested files have been found\n");
36f9ecb7 634 if (nReadFiles==0) {
635 printf("ERROR: Any file/dir found\n");
636 return kFALSE;
637 }
cb2336fe 638 }
639
640 Int_t nPtBinsCuts=dcuts[0]->GetNPtBins();
641 printf("Number of pt bins for cut object = %d\n",nPtBins);
642 Float_t *ptlimsCuts=dcuts[0]->GetPtBinLimits();
643 ptlimsCuts[nPtBinsCuts]=ptlimsCuts[nPtBinsCuts-1]+4.;
644
645 Int_t iFinBin=0;
646 for(Int_t i=0;i<nPtBinsCuts;i++){
647 if(ptlimsCuts[i]>=ptlims[iFinBin+1]) iFinBin+=1;
648 if(iFinBin>nPtBins) break;
649 if(ptlimsCuts[i]>=ptlims[iFinBin] &&
650 ptlimsCuts[i+1]<=ptlims[iFinBin+1]){
651 for(Int_t iFile=0; iFile<nReadFiles; iFile++){
652 TH1F* htemp=(TH1F*)hlist[iFile]->FindObject(Form("histMass_%d",i));
653 if(!hMass[iFinBin]){
654 hMass[iFinBin]=new TH1F(*htemp);
655 }else{
656 hMass[iFinBin]->Add(htemp);
657 }
658 }
659 }
660 }
5f1ad8a6 661 TString partname="Both";
662 if(optPartAntiPart==kParticleOnly) partname="D0";
663 if(optPartAntiPart==kAntiParticleOnly) partname="D0bar";
cb2336fe 664
5f1ad8a6 665 TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"recreate");
cb2336fe 666 outf->cd();
667 dcuts[0]->Write();
668 outf->Close();
669 return kTRUE;
b92c0195 670}
5f1ad8a6 671
470a7b55 672Bool_t LoadDstarD0piHistos(TObjArray* listFiles, TH1F** hMass){
673 //
674 Int_t nFiles=listFiles->GetEntries();
675 TList **hlist=new TList*[nFiles];
676 AliRDHFCutsDStartoKpipi** dcuts=new AliRDHFCutsDStartoKpipi*[nFiles];
677 Int_t nReadFiles=0;
678 for(Int_t iFile=0; iFile<nFiles; iFile++){
679 TString fName=((TObjString*)listFiles->At(iFile))->GetString();
680 TFile *f=TFile::Open(fName.Data());
681 if(!f){
682 printf("ERROR: file %s does not exist\n",fName.Data());
683 continue;
684 }
685 printf("Open File %s\n",f->GetName());
686 TString dirname="PWG3_D2H_DStarSpectra";
687 TDirectory *dir = (TDirectory*)f->Get(dirname);
688 if(!dir){
689 printf("ERROR: directory %s not found in %s\n",dirname.Data(),fName.Data());
690 continue;
691 }
692 TString listmassname="DStarPID00";
693
694 hlist[nReadFiles]=(TList*)dir->Get(listmassname);
695 TString cutsobjname="DStartoKpipiCuts";
696 dcuts[nReadFiles]=(AliRDHFCutsDStartoKpipi*)dir->Get(cutsobjname);
697 if(nReadFiles>0){
698 Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);
699 if(!sameCuts){
700 printf("ERROR: Cut objects do not match\n");
701 return kFALSE;
702 }
703 }
704 nReadFiles++;
705 }
706 if(nReadFiles<nFiles){
707 printf("WARNING: not all requested files have been found\n");
708 if (nReadFiles==0) {
709 printf("ERROR: Any file/dir found\n");
710 return kFALSE;
711 }
712 }
713 Int_t nPtBinsCuts=dcuts[0]->GetNPtBins();
714 printf("Number of pt bins for cut object = %d\n",nPtBins);
715 Float_t *ptlimsCuts=dcuts[0]->GetPtBinLimits();
716 ptlimsCuts[nPtBinsCuts]=ptlimsCuts[nPtBinsCuts-1]+4.;
717
718
719
720
721 Int_t iFinBin=0;
722 for(Int_t i=0;i<nPtBinsCuts;i++){
723 if(ptlimsCuts[i]>=ptlims[iFinBin+1]) iFinBin+=1;
724 if(iFinBin>nPtBins) break;
725 if(ptlimsCuts[i]>=ptlims[iFinBin] &&
726 ptlimsCuts[i+1]<=ptlims[iFinBin+1]){
727 for(Int_t iFile=0; iFile<nReadFiles; iFile++){
728 TH1F* htemp=(TH1F*)hlist[iFile]->FindObject(Form("histDeltaMass_%d",i));
729 if(!hMass[iFinBin]){
730 hMass[iFinBin]=new TH1F(*htemp);
731 }else{
732 hMass[iFinBin]->Add(htemp);
733 }
734 }
735 }
736 }
737 TString partname="Both";
738 if(optPartAntiPart==kParticleOnly) partname="DStar";
739 if(optPartAntiPart==kAntiParticleOnly) partname="DStarbar";
740
741 TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"recreate");
742 outf->cd();
743 dcuts[0]->Write();
744 outf->Close();
745 return kTRUE;
746}
747
5f1ad8a6 748void CompareFitTypes(TString* paths, TString* legtext,Int_t ncmp=3,TString* filenameYield=0x0){
749 //read ncmp RawYield.roots and draw them together
750 //set the global variable nevents before running
751 //arguments:
752 // - paths= vector of ncmp dimension with the paths of the file RawYield.root
753 // - legtext= vector of ncmp dimension with the label for the legend
754 // - ncmp= number of files to compare (default is 3)
755 // - filenameYield= optional ncmp-dimensional array with the difference between the names of the files to be compared (useful if the 2 files are in the same directory but have different names)
756
757 gStyle->SetOptStat(0);
758 gStyle->SetOptTitle(0);
759 gStyle->SetCanvasColor(0);
760 gStyle->SetFrameFillColor(0);
761 gStyle->SetTitleFillColor(0);
762 gStyle->SetFrameBorderMode(0);
763
764 if(!filenameYield) filenameYield=new TString[ncmp];
765
766 for(Int_t k=0;k<ncmp;k++){
767 if(!filenameYield) filenameYield[k]="RawYield.root";
768 filenameYield[k].Prepend(paths[k]);
769 }
770
771 TCanvas* cSig=new TCanvas("cSig","Results",1200,600);
772 cSig->Divide(3,1);
773 TCanvas* cBkgN=new TCanvas("cBkgN","Background normalized to sigma",400,600);
774 TCanvas* cDiffS=new TCanvas("cDiffS","Difference",1200,600);
775 cDiffS->Divide(2,1);
776 TCanvas *cChi2=new TCanvas("cChi2","reduced chi square",600,600);
777
778 TLegend* leg1=new TLegend(0.4,0.7,0.89,0.89);
779 leg1->SetFillColor(0);
780 TLegend* leg2=(TLegend*)leg1->Clone();
781 TLegend* leg3=(TLegend*)leg1->Clone();
782 TLegend* leg4=new TLegend(0.4,0.6,0.8,0.89);
783 leg4->SetFillColor(0);
784
785 for(Int_t iTypes=0;iTypes<ncmp;iTypes++){
786 TFile* fin=new TFile(filenameYield[iTypes]);
787 if(!fin){
788 printf("WARNING: %s not found",filenameYield[iTypes].Data());
789 continue;
790 }
791
792 TH1F* hSignal=(TH1F*)fin->Get("hSignal");
793 TH1F* hBackground=(TH1F*)fin->Get("hBackground");
794 TH1F* hBackgroundNormSigma=(TH1F*)fin->Get("hBackgroundNormSigma");
795 TH1F* hSignificance=(TH1F*)fin->Get("hSignificance");
796 hSignal->SetName(Form("%s%d",hSignal->GetName(),iTypes));
797 hBackground->SetName(Form("%s%d",hBackground->GetName(),iTypes));
798 hBackgroundNormSigma->SetName(Form("%s%d",hBackgroundNormSigma->GetName(),iTypes));
799 hSignificance->SetName(Form("%s%d",hSignificance->GetName(),iTypes));
800
801 hSignal->SetMarkerColor(iTypes+2);
802 hSignal->SetLineColor(iTypes+2);
803 hBackground->SetMarkerColor(iTypes+2);
804 hBackground->SetLineColor(iTypes+2);
805 hBackgroundNormSigma->SetMarkerColor(iTypes+2);
806 hBackgroundNormSigma->SetLineColor(iTypes+2);
807 hSignificance->SetMarkerColor(iTypes+2);
808 hSignificance->SetLineColor(iTypes+2);
809
810 TLegendEntry* ent4=leg4->AddEntry(hSignal,Form("%s",legtext[iTypes].Data()),"PL");
811 ent4->SetTextColor(hSignal->GetMarkerColor());
812
813 if(ncmp==nsamples){
814 printf("Info: Normalizing signal, background and significance to the number of events\n");
815 hSignal->Scale(1./nevents[iTypes]);
816 hBackground->Scale(1./nevents[iTypes]);
817 hBackgroundNormSigma->Scale(1./nevents[iTypes]);
818 hSignificance->Scale(1./TMath::Sqrt(nevents[iTypes]));
819 }
820
821 if(iTypes==0){
822 cSig->cd(1);
823 hSignal->DrawClone("P");
824 cSig->cd(2);
825 hBackground->DrawClone("P");
826 cSig->cd(3);
827 hSignificance->DrawClone("P");
828 cBkgN->cd();
829 hBackgroundNormSigma->DrawClone("P");
830 } else{
831 cSig->cd(1);
832 hSignal->DrawClone("Psames");
833 cSig->cd(2);
834 hBackground->DrawClone("Psames");
835 cSig->cd(3);
836 hSignificance->DrawClone("Psames");
837 cBkgN->cd();
838 hBackgroundNormSigma->DrawClone("Psames");
839 }
840
841 TH1F* hRelErrSig=(TH1F*)fin->Get("hRelErrSig");
842 TH1F* hInvSignif=(TH1F*)fin->Get("hInvSignif");
843 hRelErrSig->SetName(Form("%s%d",hRelErrSig->GetName(),iTypes));
844 hInvSignif->SetName(Form("%s%d",hInvSignif->GetName(),iTypes));
845
846 hRelErrSig->SetMarkerColor(iTypes+2);
847 hRelErrSig->SetLineColor(iTypes+2);
848 hInvSignif->SetMarkerColor(iTypes+2);
849 hInvSignif->SetLineColor(iTypes+2);
850
851 TLegendEntry* ent1=leg1->AddEntry(hRelErrSig,Form("From Fit (%s)",legtext[iTypes].Data()),"P");
852 ent1->SetTextColor(hRelErrSig->GetMarkerColor());
853 ent1=leg1->AddEntry(hInvSignif,Form("1/Significance (%s)",legtext[iTypes].Data()),"PL");
854 ent1->SetTextColor(hInvSignif->GetMarkerColor());
855
856 cDiffS->cd(1);
857 if(iTypes==0){
858 hRelErrSig->DrawClone("P");
859 hInvSignif->DrawClone();
860 } else{
861 hRelErrSig->DrawClone("Psames");
862 hInvSignif->DrawClone("sames");
863 }
864
865 TH1F* hNDiffCntSig1=(TH1F*)fin->Get("hNDiffCntSig1");
866 TH1F* hNDiffCntSig2=(TH1F*)fin->Get("hNDiffCntSig2");
867 hNDiffCntSig1->SetName(Form("%s%d",hNDiffCntSig1->GetName(),iTypes));
868 hNDiffCntSig2->SetName(Form("%s%d",hNDiffCntSig2->GetName(),iTypes));
869
870 hNDiffCntSig1->SetMarkerColor(iTypes+2);
871 hNDiffCntSig1->SetLineColor(iTypes+2);
872 hNDiffCntSig2->SetMarkerColor(iTypes+2);
873 hNDiffCntSig2->SetLineColor(iTypes+2);
874 TLegendEntry* ent2=leg2->AddEntry(hNDiffCntSig1,Form("From Counting1 (%s)",legtext[iTypes].Data()),"PL");
875 ent2->SetTextColor(hNDiffCntSig1->GetMarkerColor());
876 ent2=leg2->AddEntry(hNDiffCntSig2,Form("From Counting2 (%s)",legtext[iTypes].Data()),"PL");
877 ent2->SetTextColor(hNDiffCntSig2->GetMarkerColor());
878
879 cDiffS->cd(2);
880 if(iTypes==0){
881 hNDiffCntSig1->DrawClone();
882 hNDiffCntSig2->DrawClone();
883 }else{
884 hNDiffCntSig1->DrawClone("sames");
885 hNDiffCntSig2->DrawClone("sames");
886 }
887
888 TGraph* grReducedChiSquare=(TGraph*)fin->Get("grReducedChiSquare");
889 grReducedChiSquare->SetName(Form("%s%d",grReducedChiSquare->GetName(),iTypes));
890
891 grReducedChiSquare->SetMarkerColor(iTypes+2);
892 grReducedChiSquare->SetLineColor(iTypes+2);
893 TLegendEntry* ent3=leg3->AddEntry(grReducedChiSquare,Form("%s",legtext[iTypes].Data()),"PL");
894 ent3->SetTextColor(grReducedChiSquare->GetMarkerColor());
895
896 cChi2->cd();
897 if(iTypes==0) grReducedChiSquare->DrawClone("AP");
898 else grReducedChiSquare->DrawClone("P");
899 }
900
901 cSig->cd(1);
902 leg4->Draw();
903
904 cDiffS->cd(1);
905 leg1->Draw();
906
907 cDiffS->cd(2);
908 leg2->Draw();
909
910 cChi2->cd();
911 leg3->Draw();
912
913 TFile* fout=new TFile("ComparisonRawYield.root","RECREATE");
914 fout->cd();
915 cDiffS->Write();
916 cChi2->Write();
917 fout->Close();
918}
1918f6ea 919
920
921TH1F* RebinHisto(TH1F* hOrig, Int_t reb, Int_t firstUse){
922 // Rebin histogram, from bin firstUse to lastUse
923 // Use all bins if firstUse=-1
924
925 Int_t nBinOrig=hOrig->GetNbinsX();
926 Int_t firstBinOrig=1;
927 Int_t lastBinOrig=nBinOrig;
928 Int_t nBinOrigUsed=nBinOrig;
929 Int_t nBinFinal=nBinOrig/reb;
930 if(firstUse>=1){
931 firstBinOrig=firstUse;
932 nBinFinal=(nBinOrig-firstUse+1)/reb;
933 nBinOrigUsed=nBinFinal*reb;
934 lastBinOrig=firstBinOrig+nBinOrigUsed-1;
935 }else{
936 Int_t exc=nBinOrigUsed%reb;
937 if(exc!=0){
938 nBinOrigUsed-=exc;
939 firstBinOrig+=exc/2;
940 lastBinOrig=firstBinOrig+nBinOrigUsed-1;
941 }
942 }
943
944 printf("Rebin from %d bins to %d bins -- Used bins=%d in range %d-%d\n",nBinOrig,nBinFinal,nBinOrigUsed,firstBinOrig,lastBinOrig);
945 Float_t lowLim=hOrig->GetXaxis()->GetBinLowEdge(firstBinOrig);
946 Float_t hiLim=hOrig->GetXaxis()->GetBinUpEdge(lastBinOrig);
947 TH1F* hRebin=new TH1F(Form("%s-rebin",hOrig->GetName()),hOrig->GetTitle(),nBinFinal,lowLim,hiLim);
948 Int_t lastSummed=firstBinOrig-1;
949 for(Int_t iBin=1;iBin<=nBinFinal; iBin++){
950 Float_t sum=0.;
951 for(Int_t iOrigBin=0;iOrigBin<reb;iOrigBin++){
952 sum+=hOrig->GetBinContent(lastSummed+1);
953 lastSummed++;
954 }
955 hRebin->SetBinContent(iBin,sum);
956 }
957 return hRebin;
958}