]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG0/multPbPb/correct.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Changes for event-by-event and cosmetics
[u/mrichter/AliRoot.git] / PWG0 / multPbPb / correct.C
CommitLineData
f8416b14 1#include "AliAnalysisMultPbTrackHistoManager.h"
2#include "TH2D.h"
3#include "TH1D.h"
4#include "TF1.h"
5#include "TCanvas.h"
6#include "TFile.h"
7#include "TSystem.h"
8#include "TROOT.h"
9#include <iostream>
10#include "TDatabasePDG.h"
11#include "AliPhysicsSelection.h"
12#include "AliESDtrackCuts.h"
2bbfd8c6 13#include "AliAnalysisMultPbCentralitySelector.h"
4d0aa70f 14#include "TLegend.h"
abd808b9 15#include "AliBWTools.h"
f8416b14 16
17using namespace std;
a23f7c97 18
19AliAnalysisMultPbTrackHistoManager * hManData = 0;
20AliAnalysisMultPbTrackHistoManager * hManCorr = 0;
21TH2D * hEvStatData = 0;
22TH2D * hEvStatCorr = 0;
23
bcc49144 24const Int_t kHistoFitCompoments = 2;
f8416b14 25TH1D * gHistoCompoments[kHistoFitCompoments];
26
aa6151eb 27void LoadLibs( Bool_t debug=0);
a23f7c97 28void LoadData(TString dataFolder, TString correctionFolder);
29void SetStyle();
e0376287 30void CheckSecondaries(Double_t & fracWeak, Double_t &fracMaterial);
72491d7c 31void CheckSanity();
9441cdaa 32void CheckCorrections();
a23f7c97 33void ShowAcceptanceInVzSlices() ;
f8416b14 34TH1D * GetRatioIntegratedFractions (TH1 * hNum, TH1 * hDenum) ;
e0376287 35TH1D * GetCumulativeHisto (TH1 * h) ;
f8416b14 36static Double_t HistoSum(const double * x, const double* p);
37TF1 * GetFunctionHistoSum() ;
38TF1 * GetMTExp(Float_t norm=68, Float_t t=25) ;
39TF1 * GetHagedorn(Float_t norm=68, Float_t pt0=25, Float_t n=13) ;
40TF1 * GetLevy(Double_t temp=0.1, Double_t n=7, Double_t norm=10, const char * name="fLevy") ;
e0376287 41void PrintCanvas(TCanvas* c,const TString formats) ;
f8416b14 42
e0376287 43// global switches
44Bool_t doPrint=kFALSE;// disable PrintCanvas
55e8e56a 45Float_t zmin = -10;
46Float_t zmax = 10;
72491d7c 47Float_t etaMin = -0.5;
48Float_t etaMax = 0.5;
bfae2924 49Int_t gCentralityBin = -1;
f8416b14 50
aa6151eb 51#define CORRECT_2D
f8416b14 52
72491d7c 53void correct(TString dataFolder = "./output/LHC10g2d_130844_V0M_bin_10/", TString correctionFolder = "./output/LHC10g2a_130844_V0M_bin_10/",
bfae2924 54 Float_t zminl=-10, Float_t zmaxl=10, Float_t etaMinl = -0.5, Float_t etaMaxl=0.5, Float_t npart=381.3, Float_t scaleWeakByHand = -1, Int_t centrBin = -1) {
72491d7c 55
bfae2924 56 // Set globals
72491d7c 57 zmin = zminl;
58 zmax = zmaxl;
59 etaMin = etaMinl;
60 etaMax = etaMaxl;
bfae2924 61 gCentralityBin = centrBin;
a23f7c97 62
63 // Load stuff and set some styles
64 LoadLibs();
65 LoadData(dataFolder,correctionFolder);
66 SetStyle();
3b8cbf2d 67 // ShowAcceptanceInVzSlices();
68 // return;
a23f7c97 69
70 // TODO add some cool printout for cuts and centrality selection
71
72491d7c 72 CheckSanity();
e0376287 73
74 Double_t fractionWeak = 1, fractionMaterial=1;
bcc49144 75 CheckSecondaries(fractionWeak, fractionMaterial);
e0376287 76 cout << "Rescaling secondaries correction, weak: " << fractionWeak << ", material: " << fractionMaterial <<endl;
9441cdaa 77 if(scaleWeakByHand >= 0) {
78 fractionWeak = scaleWeakByHand;
79 cout << "Scaling Strangeness by hand (factor " << fractionWeak <<")" << endl;
80
81 }
a23f7c97 82 // Some shorthands
9441cdaa 83
aa6151eb 84#if defined CORRECT_1D
72491d7c 85 TH1D * hDataPt = (TH1D*) hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRec, etaMin,etaMax,zmin,zmax)->Clone("hDataPt");
9441cdaa 86 TH1D * hMCPtGen = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen, etaMin,etaMax,zmin,zmax);
72491d7c 87 //zTH1D * hMCPtGen = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen, etaMin,etaMax,zmin,zmax);
88 TH1D * hMCPtRec = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRec, etaMin,etaMax,zmin,zmax);
89 TH1D * hMCPtPri = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim, etaMin,etaMax,zmin,zmax);
90 TH1D * hMCPtSeM = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat, etaMin,etaMax,zmin,zmax);
91 TH1D * hMCPtSeW = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak, etaMin,etaMax,zmin,zmax);
92 TH1D * hMCPtFak = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecFake, etaMin,etaMax,zmin,zmax);
aa6151eb 93#elif defined CORRECT_2D
72491d7c 94 TH1 * hDataPt = (TH2D*) hManData->GetHistoPtVz(AliAnalysisMultPbTrackHistoManager::kHistoRec, etaMin,etaMax)->Clone("hDataPt");
95 TH1 * hMCPtGen = (TH2D*) hManCorr->GetHistoPtVz(AliAnalysisMultPbTrackHistoManager::kHistoGen, etaMin,etaMax);
96 TH1 * hMCPtRec = (TH2D*) hManCorr->GetHistoPtVz(AliAnalysisMultPbTrackHistoManager::kHistoRec, etaMin,etaMax);
97 TH1 * hMCPtPri = (TH2D*) hManCorr->GetHistoPtVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim, etaMin,etaMax);
98 TH1 * hMCPtSeM = (TH2D*) hManCorr->GetHistoPtVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat, etaMin,etaMax);
99 TH1 * hMCPtSeW = (TH2D*) hManCorr->GetHistoPtVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak, etaMin,etaMax);
100 TH1 * hMCPtFak = (TH2D*) hManCorr->GetHistoPtVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake, etaMin,etaMax);
aa6151eb 101#endif
a23f7c97 102
a23f7c97 103 // hMCPtRec->Draw("same");
104
105 TCanvas * cMC = new TCanvas ("cMC", "Monte Carlo");
e0376287 106 cMC->SetLogy();
f8416b14 107 cMC->cd();
72491d7c 108 // hMCPtSeW->Add(hMCPtSeM);
109 // hMCPtSeW->Divide(hMCPtRec);
110
111
a23f7c97 112 hMCPtGen ->Draw();
113 hMCPtRec ->Draw("same");
114 hMCPtPri ->Draw("same");
115 hMCPtSeM ->Draw("same");
116 hMCPtSeW ->Draw("same");
117 hMCPtFak ->Draw("same");
72491d7c 118#ifdef CORRECT_1D
4d0aa70f 119 hMCPtGen ->GetXaxis()->SetRangeUser(0,4.5);
120 hMCPtGen ->GetYaxis()->SetRangeUser(0.1,1e4);
72491d7c 121#endif
4d0aa70f 122 TLegend * lMC = new TLegend(0.505034, 0.59965, 0.877517, 0.926573,"Monte Carlo");
123 lMC->AddEntry( hMCPtGen, "Generated");
124 lMC->AddEntry( hMCPtRec, "All Rec");
125 lMC->AddEntry( hMCPtPri, "Rec Primaries");
126 lMC->AddEntry( hMCPtSeM, "Rec Sec. Material");
127 lMC->AddEntry( hMCPtSeW, "Rec Sec. Weak");
128 lMC->AddEntry( hMCPtFak, "Rec Fakes");
129 lMC->Draw();
130
a23f7c97 131
132 cout << "Fake/All Rec = " << hMCPtFak->Integral()/hMCPtRec->Integral() << endl;
133 cout << "SM/All Rec = " << hMCPtSeM->Integral()/hMCPtRec->Integral() << endl;
134 cout << "SW/All Rec = " << hMCPtSeW->Integral()/hMCPtRec->Integral() << endl;
135
136
137 // Compute efficiency and subtract secondaries and fakes after rescaling to data
138 // PRIM_DATA = ALL_DATA - SEC_MC/ALL_MC*ALL_DATA - FAK_MC/ALL_MC*ALL_DATA
139 // TRUE_DATA = PRIM_DATA * GEN_MC/PRIM_MC
140
141 TH1D * hEffPt = (TH1D*) hMCPtPri->Clone("hEffPt");
142 hEffPt->Divide(hMCPtPri,hMCPtGen,1,1,"B");
143
4d0aa70f 144 TH1D * hCorSeM = (TH1D*) hMCPtSeM->Clone("hCorSeM");
a23f7c97 145 hCorSeM->Divide(hMCPtSeM,hMCPtRec,1,1,"B");
e0376287 146 hCorSeM->Scale(fractionMaterial);// rescale material correction
4d0aa70f 147 hCorSeM->Multiply(hDataPt);
a23f7c97 148
4d0aa70f 149 TH1D * hCorSeW = (TH1D*) hMCPtSeW->Clone("hCorSeW");
a23f7c97 150 hCorSeW->Divide(hMCPtSeW,hMCPtRec,1,1,"B");
151 hCorSeW->Scale(fractionWeak);// rescale weak correction
152 hCorSeW->Multiply(hDataPt);
153
4d0aa70f 154 TH1D * hCorFak = (TH1D*) hMCPtFak->Clone("hCorFak");
a23f7c97 155 hCorFak->Divide(hMCPtFak,hMCPtRec,1,1,"B");
156 hCorFak->Multiply(hDataPt);
157
158 TH1D * hCorrected = (TH1D*) hDataPt->Clone("hCorrected");
159 hCorrected->Add(hCorSeM,-1);
160 hCorrected->Add(hCorSeW,-1);
161 hCorrected->Add(hCorFak,-1);
162 hCorrected->Divide(hEffPt);
163 hCorrected->SetMarkerStyle(kOpenStar);
164
4d0aa70f 165
166 TCanvas * cCorrections = new TCanvas("cCorrections", "cCorrections");
167
168 hEffPt->Draw();
169 // hCorSeM->Draw();
170 // hCorSeM->SetLineColor(kRed);
171 // hCorSeM->SetMarkerColor(kRed);
172 // hMCPtSeM->Draw("same");
173 // hCorSeW->Draw("same");
174 // hCorSeW->SetLineColor(kRed);
175 // hCorSeW->SetMarkerColor(kRed);
176 // hMCPtSeW->Draw("same");
177
178
aa6151eb 179 // hDataPt->Draw();
180 // return;
72491d7c 181 TCanvas * cdata = new TCanvas ("cData", "Data");
182 cdata->SetLogy();
aa6151eb 183 cdata->cd();
184#ifdef CORRECT_2D
185 Int_t minProj = hDataPt->GetYaxis()->FindBin(zmin);
186 Int_t maxProj = hDataPt->GetYaxis()->FindBin(zmax-0.0001);
187
188 cout << minProj << "-" << maxProj << endl;
189
190 // This correction accounts for the vertex cut;
191
192 hDataPt = ((TH2D*)hDataPt )->ProjectionX("_px",minProj,maxProj,"E");
193 hCorrected = ((TH2D*)hCorrected)->ProjectionX("_px",minProj,maxProj,"E");
194 hMCPtGen = ((TH2D*)hMCPtGen )->ProjectionX("_px",minProj,maxProj,"E");
195
196 Float_t vertexCutCorrection =
197 hManCorr->GetHistoVzEvent(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Integral(-1,-1) /
198 hManCorr->GetHistoVzEvent(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Integral(minProj,maxProj) ;
199 // cout << vertexCutCorrection << " " << hMCPtGen->Integral(-1,-1) << " " << hMCPtPri->Integral() << endl;
200 // vertexCutCorrection /= hMCPtGen->Integral(-1,-1);
55e8e56a 201 vertexCutCorrection = 1; // FIXME
aa6151eb 202 cout << "Vertex cut correction " << vertexCutCorrection << " (Efficiency " << 1./vertexCutCorrection << ")" << endl;
203
204 hDataPt ->Scale(1.,"width");
205 hCorrected ->Scale(vertexCutCorrection,"width");
206 hMCPtGen ->Scale(1.,"width");
207#endif
208
a23f7c97 209 hDataPt->Draw();
4d0aa70f 210 hDataPt ->GetXaxis()->SetRangeUser(0,4.5);
211 hDataPt ->GetYaxis()->SetRangeUser(0.1,1e4);
a23f7c97 212 hCorrected->SetLineColor(kBlack);
213 hCorrected->SetMarkerColor(kBlack);
214 hCorrected->Draw("same");
215 hMCPtGen->DrawClone("same");
bcc49144 216 // TF1 * f = GetLevy();
217 TF1 * f = GetHagedorn();
a23f7c97 218 hCorrected->Fit(f,"", "same");
aa6151eb 219 hCorrected->Fit(f,"IME", "same",0,2);
4d0aa70f 220 cout << "dN/deta (function) = " << f->Integral(0,100) << " +- " << f->IntegralError(0,100) << endl;
9441cdaa 221 cout << "dN/deta (data>0.15) = " << hCorrected->Integral(4,-1,"width") << endl;//
4d0aa70f 222 cout << "dN/deta (func+data) = " << f->Integral(0,0.1) + hCorrected->Integral(3,-1,"width") << endl;//
72491d7c 223 cout << "dN/deta (func 0.15+data) = " << f->Integral(0,0.15) + hCorrected->Integral(4,-1,"width") << endl;//
eef42d18 224 cout << "Generated dN/deta (correction) = " << hMCPtGen->Integral("width") << endl;
72491d7c 225 cout << "Generated dN/deta (correction, <0.13) = " << hMCPtGen->Integral(1,2,"width") << endl;
226 cout << "-------" << endl;
9441cdaa 227 cout << "Rescaling secondaries correction, weak: " << fractionWeak << ", material: " << fractionMaterial <<endl;
72491d7c 228 Double_t errorData = 0;
9441cdaa 229 Float_t dNdeta = (f->Integral(0,0.1) + hCorrected->IntegralAndError(1,-1,errorData, "width")) / (etaMax-etaMin);
72491d7c 230 Float_t dNdetaE = TMath::Sqrt(f->IntegralError(0,0.1)*f->IntegralError(0,0.1) + errorData*errorData) / (etaMax-etaMin);
bcc49144 231 cout << "dN/deta " << dNdeta << " +- " << dNdetaE << endl;
232 cout << "(dN/deta)/0.5 Npart " << dNdeta/npart*2 << " +- " << dNdetaE/npart*2 << endl;
abd808b9 233 Double_t mean, meanerr;
234 AliBWTools::GetMeanDataAndExtrapolation(hCorrected, f, mean, meanerr,0,1000);
235 cout << "Mean Pt " << mean << "+-" << meanerr << endl;
72491d7c 236
abd808b9 237
72491d7c 238 TH1D * hDataGen = hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen, etaMin,etaMax,zmin,zmax);
eef42d18 239 cout << "Generated dN/deta (data) = " << hDataGen->Integral("width") << endl;
4d0aa70f 240 hDataGen->Draw("same");
241 TLegend * l = new TLegend(0.520134, 0.676573, 0.885906, 0.923077,"137161, p1+++");
242 l->AddEntry(hDataPt, "Raw data");
243 l->AddEntry(hCorrected, "Corrected data");
244 l->AddEntry(hMCPtGen, "Monte Carlo (generated)");
bcc49144 245 l->AddEntry(f, "Hagedorn Fit");
4d0aa70f 246 l->Draw();
abd808b9 247
248
249 Double_t errXCheck = 0, yieldXCheck=0;
250 AliBWTools::GetYield(hCorrected,f,yieldXCheck,errXCheck);
251 cout << "BWTools crosscheck: " << yieldXCheck << "+-" << errXCheck << endl;
252
253
a23f7c97 254}
255
e0376287 256void CheckSecondaries(Double_t &fracWeak, Double_t &fracMaterial) {
a23f7c97 257 // Returns the fraction you need to rescale the secondaries from weak decays for.
258
259 // Some shorthands
bcc49144 260 TH2D * hDataDCAPt = hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec );
f8416b14 261 // TH1D * hMCDCAGen = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen );
bcc49144 262 TH2D * hMCDCARecPt = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec );
263 TH2D * hMCDCAPriPt = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim );
264 TH2D * hMCDCASWPt = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak );
265 TH2D * hMCDCASMPt = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat );
266 TH2D * hMCDCAFakPt = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake );
a23f7c97 267
268
269 TCanvas * cCumulative = new TCanvas("cDCAculumative","DCA cumulative distributions");
f8416b14 270 cCumulative->cd();
bcc49144 271 GetCumulativeHisto(hMCDCAPriPt )->Draw();
272 GetRatioIntegratedFractions(hMCDCASWPt, hMCDCARecPt )->Draw("same");
273 GetRatioIntegratedFractions(hMCDCASMPt, hMCDCARecPt )->Draw("same");
274 GetRatioIntegratedFractions(hMCDCAPriPt,hMCDCARecPt )->Draw("same");
a23f7c97 275
276
277 TCanvas * c1 = new TCanvas("cDCAFit", "Fit to the DCA distributions");
bfae2924 278 // c1->Divide(4,2);
a23f7c97 279 c1->SetLogy();
280 // Draw all
e0376287 281 // hDataDCA->Draw();
282 // hMCDCARec ->Draw("same");
283 // hMCDCAPri ->Draw("same");
284 // hMCDCASW ->Draw("same");
285 // hMCDCASM ->Draw("same");
286 // hMCDCAFak ->Draw("same");
287 // return;
a23f7c97 288
e0376287 289 // Fit the DCA distribution. Uses a TF1 made by summing histograms
bcc49144 290 TH2D * hMCPrimSMFakPt = (TH2D*) hMCDCAPriPt->Clone("hMCPrimSMFak");
291 hMCPrimSMFakPt->Add(hMCDCASMPt);
292 hMCPrimSMFakPt->Add(hMCDCAFakPt);
f8416b14 293
e0376287 294 // Set the components which are used in HistoSum, the static
295 // function for GetFunctionHistoSum
bcc49144 296 // Project onti DCA axis
abd808b9 297 //const Int_t ptBinsFit[] = {3,5,7,9,11,15,19,23,31,-1};
bfae2924 298 const Int_t ptBinsFit[] = {3,20,-1};
bcc49144 299 Int_t ibinPt = -1;
300 while(ptBinsFit[(++ibinPt)+1]!=-1){
301 c1->cd(ibinPt+1);
302 gPad->SetLogy();
303 Int_t minPtBin=ptBinsFit[ibinPt];
304 Int_t maxPtBin=ptBinsFit[ibinPt+1]-1;
305
306 TH1D * hDataDCA = (TH1D*) hDataDCAPt ->ProjectionX(TString(hDataDCAPt ->GetName())+long(ibinPt)+"_px",minPtBin,maxPtBin,"E")->Clone();
307 TH1D * hMCPrimSMFak = (TH1D*) hMCPrimSMFakPt->ProjectionX(TString(hMCPrimSMFakPt->GetName())+long(ibinPt)+"_px",minPtBin,maxPtBin,"E")->Clone();
308 TH1D * hMCDCASW = (TH1D*) hMCDCASWPt ->ProjectionX(TString(hMCDCASWPt ->GetName())+long(ibinPt)+"_px",minPtBin,maxPtBin,"E")->Clone();
309 TH1D * hMCDCASM = (TH1D*) hMCDCASMPt ->ProjectionX(TString(hMCDCASMPt ->GetName())+long(ibinPt)+"_px",minPtBin,maxPtBin,"E")->Clone();
310 gHistoCompoments[0] = (TH1D*) hMCPrimSMFakPt->ProjectionX(TString(hMCPrimSMFakPt->GetName())+long(ibinPt)+"_px",minPtBin,maxPtBin,"E")->Clone();
311 gHistoCompoments[1] = (TH1D*) hMCDCASWPt ->ProjectionX(TString(hMCDCASWPt ->GetName())+long(ibinPt)+"_px",minPtBin,maxPtBin,"E")->Clone();
312 // gHistoCompoments[2] = (TH1D*) hMCDCASMPt ->ProjectionX("_px",minPtBin,maxPtBin,"E")->Clone();
313 TF1 * fHistos = GetFunctionHistoSum();
314 fHistos->SetParameters(1,1);
315 // fHistos->FixParameter(2,1);
316 fHistos->SetLineColor(kRed);
317 hDataDCA= (TH1D*) hDataDCAPt->ProjectionX(TString(hDataDCAPt->GetName())+long(ibinPt)+"_px",minPtBin,maxPtBin,"E")->Clone();
318 hDataDCA->Draw();
319 // Fit!
320 hDataDCA->Fit(fHistos,"0Q","",-2,2);
321 // Rescale the components and draw to see how it looks like
322 hMCPrimSMFak->Scale(fHistos->GetParameter(0));
323 hMCDCASW ->Scale(fHistos->GetParameter(1));
324 hMCDCASM ->Scale(fHistos->GetParameter(2));
325 hMCPrimSMFak->Draw("same");
326 hMCDCASW ->Draw("same");
327 hMCDCASM ->Draw("same");
328 TH1D* hMCSum = (TH1D*) hMCPrimSMFak->Clone();
329 hMCSum->Add(hMCDCASW);
330 // hMCSum->Add(hMCDCASM);
331 hMCSum->SetLineColor(kRed);
332 hMCSum->SetMarkerStyle(1);
333 hMCSum->Draw("lhist,same");
334 // fHistos->Draw("same");
335 // compute scaling factors
336 // fracWeak = fHistos->GetParameter(1)/(fHistos->GetParameter(0)+fHistos->GetParameter(1));
337 cout << hDataDCAPt->GetYaxis()->GetBinLowEdge(minPtBin) << " - ";
338 cout << hDataDCAPt->GetYaxis()->GetBinLowEdge(maxPtBin+1) << " = ";
339 cout << 2*fHistos->GetParameter(1)/(fHistos->GetParameter(0)+fHistos->GetParameter(1));
340 cout << " ("<< fHistos->GetChisquare() <<")";
341 cout << " " << fHistos->GetParameter(0) << "," << fHistos->GetParameter(1) << endl;
342
343
344 fracWeak = 2*fHistos->GetParameter(1)/(fHistos->GetParameter(0)+fHistos->GetParameter(1));
345 // fracMaterial = fHistos->GetParameter(2)/(fHistos->GetParameter(0);
346 }
a23f7c97 347
348}
349
72491d7c 350void CheckSanity() {
351 // compares various distributions in data and in MC
e0376287 352 TCanvas * c1 = new TCanvas("cVz", "Vertex Z distribution");
353 c1->cd();
9441cdaa 354 TH1 * hData = hManData->GetHistoVzEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec );
355 TH1 * hCorr = hManCorr->GetHistoVzEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec );
72491d7c 356 hData->SetLineColor (kBlack);
357 hData->SetMarkerColor(kBlack);
358 hData->SetMarkerStyle(kFullCircle);
359 hCorr->SetLineColor (kRed);
360 hCorr->SetMarkerColor(kRed);
361 hCorr->SetMarkerStyle(kOpenSquare);
362 TLegend * l = new TLegend(0.575503, 0.70979, 0.86745, 0.917832);
363 l->AddEntry(hData, "Data");
364 l->AddEntry(hCorr, "MC");
e0376287 365 hCorr->Draw("");
366 hData->Draw("same");
72491d7c 367 l->Draw();
368
369
370 TCanvas * c2 = new TCanvas("cEta", "Eta distribution");
371 c2->cd();
372 hData = hManData->GetHistoEta(AliAnalysisMultPbTrackHistoManager::kHistoRec );
373 hCorr = hManCorr->GetHistoEta(AliAnalysisMultPbTrackHistoManager::kHistoRec );
374 hData->SetLineColor (kBlack);
375 hData->SetMarkerColor(kBlack);
376 hData->SetMarkerStyle(kFullCircle);
377 hCorr->SetLineColor (kRed);
378 hCorr->SetMarkerColor(kRed);
379 hCorr->SetMarkerStyle(kOpenSquare);
380 l = new TLegend(0.575503, 0.70979, 0.86745, 0.917832);
381 l->AddEntry(hData, "Data");
382 l->AddEntry(hCorr, "MC");
383 hData->Draw("");
384 hCorr->Draw("same");
385 l->Draw();
386
387 TCanvas * c3 = new TCanvas("cMult", "Multiplicity distribution");
388 c3->cd();
389 hData = hManData->GetHistoMult(AliAnalysisMultPbTrackHistoManager::kHistoRec );
390 hCorr = hManCorr->GetHistoMult(AliAnalysisMultPbTrackHistoManager::kHistoRec );
391 hData->SetLineColor (kBlack);
392 hData->SetMarkerColor(kBlack);
393 hData->SetMarkerStyle(kFullCircle);
394 hCorr->SetLineColor (kRed);
395 hCorr->SetMarkerColor(kRed);
396 hCorr->SetMarkerStyle(kOpenSquare);
397 l = new TLegend(0.575503, 0.70979, 0.86745, 0.917832);
398 l->AddEntry(hData, "Data");
399 l->AddEntry(hCorr, "MC");
400 hData->Draw("");
401 hCorr->Draw("same");
402 l->Draw();
e0376287 403
9441cdaa 404 TCanvas * c4 = new TCanvas("cStatistics", "cStats");
405 c4->cd();
406 hData = hManData->GetHistoStats();
407 hCorr = hManCorr->GetHistoStats();
408 hData->SetLineColor (kBlack);
409 hData->SetMarkerColor(kBlack);
410 // hData->SetMarkerStyle(kFullCircle);
411 hCorr->SetLineColor (kRed);
412 hCorr->SetMarkerColor(kRed);
413 // hCorr->SetMarkerStyle(kOpenSquare);
414 l = new TLegend(0.575503, 0.70979, 0.86745, 0.917832);
415 l->AddEntry(hData, "Data");
416 l->AddEntry(hCorr, "MC");
417 hData->Draw("");
418 hCorr->Draw("same");
419 l->Draw();
420
421
422}
423
424void CheckCorrections() {
425
bcc49144 426 // TH1D * hDataPt = (TH1D*) hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRec, etaMin,etaMax,zmin,zmax)->Clone("hDataPt");
9441cdaa 427 TH1D * hMCPtGen = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen, etaMin,etaMax,zmin,zmax);
428 TH1D * hMCPtRec = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRec, etaMin,etaMax,zmin,zmax);
429 TH1D * hMCPtPri = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim, etaMin,etaMax,zmin,zmax);
430 TH1D * hMCPtSeM = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat, etaMin,etaMax,zmin,zmax);
431 TH1D * hMCPtSeW = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak, etaMin,etaMax,zmin,zmax);
432 TH1D * hMCPtFak = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecFake, etaMin,etaMax,zmin,zmax);
433
434 hMCPtSeM->Divide(hMCPtRec);
435 hMCPtSeW->Divide(hMCPtRec);
436 hMCPtPri->Divide(hMCPtGen);
437 TCanvas * c1 = new TCanvas("cSecondaries", "Secondaries");
438 hMCPtSeM->Draw();
439 hMCPtSeW->Draw("same");
440
441 TCanvas * c2 = new TCanvas("cEfficiency", "Efficiency");
442 hMCPtPri->Draw();
443
e0376287 444}
445
aa6151eb 446void LoadLibs( Bool_t debug) {
a23f7c97 447
448 gSystem->Load("libVMC");
449 gSystem->Load("libTree");
450 gSystem->Load("libSTEERBase");
451 gSystem->Load("libESD");
452 gSystem->Load("libAOD");
453 gSystem->Load("libANALYSIS");
454 gSystem->Load("libANALYSISalice");
455 gSystem->Load("libCORRFW");
456 gSystem->Load("libMinuit");
457 gSystem->Load("libPWG2Spectra");
458 gSystem->Load("libPWG0base");
459
a28a49f6 460 gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWG0/multPbPb"));
f8416b14 461 gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWG1/background"));
abd808b9 462 gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWG2/SPECTRA/Fit"));
463
a23f7c97 464 // Load helper classes
465 // TODO: replace this by a list of TOBJStrings
a28a49f6 466 TString taskName("$ALICE_ROOT/PWG0/multPbPb/AliAnalysisTaskMultPbTracks.cxx+");
467 TString histoManName("$ALICE_ROOT/PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.cxx+");
468 TString centrName("$ALICE_ROOT/PWG0/multPbPb/AliAnalysisMultPbCentralitySelector.cxx+");
a23f7c97 469 TString listName("$ALICE_ROOT/PWG1/background/AliHistoListWrapper.cxx+");
470
a28a49f6 471 gSystem->ExpandPathName(taskName);
472 gSystem->ExpandPathName(histoManName);
473 gSystem->ExpandPathName(centrName);
474 gSystem->ExpandPathName(listName);
475
aa6151eb 476
a23f7c97 477 gROOT->LoadMacro(listName +(debug?"+g":""));
478 gROOT->LoadMacro(histoManName+(debug?"+g":""));
479 gROOT->LoadMacro(centrName +(debug?"+g":""));
480 gROOT->LoadMacro(taskName +(debug?"+g":""));
481
482 // Histo fitter
483 gROOT->LoadMacro("/Users/mfloris/Work/ALICE/ANALYSIS/HistoFitter/fcn.cxx+g");
f8416b14 484 gROOT->LoadMacro("/Users/mfloris/Work/ALICE/ANALYSIS/HistoFitter/AliHistoFitter.cxx+g");
a23f7c97 485
486
487}
488
489
490void SetStyle() {
491
492 hManData->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec) ->SetLineColor(kBlack);
493 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen) ->SetLineColor(kRed );
494 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec) ->SetLineColor(kRed );
495 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetLineColor(kGreen);
496 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetLineColor(kBlue );
497 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetLineColor(kBlue );
498 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetLineColor(kCyan );
499
500 hManData->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec) ->SetMarkerColor(kBlack);
501 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen) ->SetMarkerColor(kRed );
502 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec) ->SetMarkerColor(kRed );
503 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerColor(kGreen);
504 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerColor(kBlue );
bcc49144 505 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerColor(kCyan );
506 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerColor(kYellow );
a23f7c97 507
508 hManData->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec) ->SetMarkerStyle(kFullCircle);
509 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen) ->SetMarkerStyle(kFullSquare);
510 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec) ->SetMarkerStyle(kOpenSquare);
511 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerStyle(kOpenSquare);
512 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerStyle(kOpenSquare);
513 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerStyle(kOpenCircle);
514 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerStyle(kOpenSquare);
515
516 hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec) ->SetLineColor(kBlack);
517 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen) ->SetLineColor(kRed );
518 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec) ->SetLineColor(kRed );
519 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetLineColor(kGreen);
520 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetLineColor(kBlue );
bcc49144 521 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetLineColor(kCyan);
522 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetLineColor(kYellow );
a23f7c97 523
524 hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec) ->SetMarkerColor(kBlack);
525 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen) ->SetMarkerColor(kRed );
526 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec) ->SetMarkerColor(kRed );
527 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerColor(kGreen);
528 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerColor(kBlue );
bcc49144 529 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerColor(kCyan );
530 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerColor(kYellow );
a23f7c97 531
532 hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec) ->SetMarkerStyle(kFullCircle);
533 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen) ->SetMarkerStyle(kFullSquare);
534 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec) ->SetMarkerStyle(kOpenSquare);
535 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerStyle(kOpenSquare);
536 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerStyle(kOpenSquare);
537 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerStyle(kOpenSquare);
538 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerStyle(kOpenSquare);
539
540
541}
542
543void LoadData(TString dataFolder, TString correctionFolder){
544
bfae2924 545 TString file = "/multPbPbtracks.root";
546 if (gCentralityBin != -1 ) file.ReplaceAll(".root", Form("_%2.2d.root", gCentralityBin));
a23f7c97 547 // Get histo manager for data and MC + stat histos
bfae2924 548 TFile * fData = new TFile(dataFolder+file);
549 TFile * fCorr = new TFile(correctionFolder+file);
72491d7c 550 TFile * fStatData = new TFile(dataFolder+"/event_stat.root");
551 TFile * fStatCorr = new TFile(correctionFolder+"/event_stat.root");
a23f7c97 552
553 hManData = (AliAnalysisMultPbTrackHistoManager*) fData->Get("histoManager");
554 hManCorr = (AliAnalysisMultPbTrackHistoManager*) fCorr->Get("histoManager");
555 AliESDtrackCuts * cutsData = (AliESDtrackCuts*) fData->Get("AliESDtrackCuts");
556 AliESDtrackCuts * cutsCorr = (AliESDtrackCuts*) fCorr->Get("AliESDtrackCuts");
9441cdaa 557
558 // FIXME: implement operator= and Print in AliESDtrackCuts
559 // if (cutsData != cutsCorr) {
560 // cout << "ERROR: MC and data do not have the same cuts" << endl;
561 // // }
a23f7c97 562 cutsData->Print();
563 hEvStatData = (TH2D*) fStatData->Get("fHistStatistics");
564 hEvStatCorr = (TH2D*) fStatCorr->Get("fHistStatistics");
565
aa6151eb 566 AliAnalysisMultPbCentralitySelector * centrData = (AliAnalysisMultPbCentralitySelector*) fData->Get("Cuts");
567 if(!centrData) {
568 cout << "ERROR: cannot read centrality data" << endl;
569 }
570 centrData->Print();
a23f7c97 571 // Normalize
572 Int_t irowGoodTrigger = 1;
573 if (hEvStatCorr && hEvStatData) {
574 // hManData->ScaleHistos(75351.36/1.015);// Nint for run 104892 estimated correcting for the trigger efficiency, multiplied for the physics selection efficiency which I'm not correcting for the time being
e0376287 575 // hManData->ScaleHistos(hEvStatData->GetBinContent(AliPhysicsSelection::kStatAccepted,irowGoodTrigger));
576 // hManCorr->ScaleHistos(hEvStatCorr->GetBinContent(AliPhysicsSelection::kStatAccepted,irowGoodTrigger));
55e8e56a 577 // hManData->ScaleHistos(hManData->GetHistoStats()->GetBinContent(AliAnalysisMultPbTrackHistoManager::kStatVtx));
578 // hManCorr->ScaleHistos(hManCorr->GetHistoStats()->GetBinContent(AliAnalysisMultPbTrackHistoManager::kStatVtx));
579 TH1D* hvzData = hManData->GetHistoVzEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec);
580 TH1D* hvzCorr = hManCorr->GetHistoVzEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec);
581 hManData->ScaleHistos(hvzData->Integral(hvzData->FindBin(zmin),hvzData->FindBin(zmax-0.001)));
582 hManCorr->ScaleHistos(hvzCorr->Integral(hvzCorr->FindBin(zmin),hvzCorr->FindBin(zmax-0.001)));
a23f7c97 583 } else {
584 cout << "WARNING!!! ARBITRARY SCALING" << endl;
585 hManData->ScaleHistos(1000);
586 hManCorr->ScaleHistos(1000);
587 }
588}
589
f8416b14 590TF1 * GetHagedorn(Float_t norm, Float_t pt0, Float_t n) {
a23f7c97 591
592 TF1 * f =0;
593 Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
594
595 f=new TF1("fHagedorn",Form("(x/sqrt(x*x+%f*%f))*x*[0]*( 1 + x/[1] )^(-[2])",mass,mass),0,10);
596 f->SetParameters(norm, pt0, n);
597 f->SetParLimits(1, 0.01, 10);
598 f->SetParNames("norm", "pt0", "n");
599 f->SetLineWidth(1);
600 return f;
601
602
603}
604
f8416b14 605TF1 * GetMTExp(Float_t norm, Float_t t) {
a23f7c97 606
607 TF1 * f =0;
608 Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
609
610 f=new TF1("fMTExp",Form("(x/sqrt(x*x+%f*%f))*x*[0]*exp(-sqrt(x*x+%f*%f)/[1])",mass,mass,mass,mass),0,10);
611 f->SetParameters(norm, t);
612 // f->SetParLimits(1, 0.01);
613 f->SetParNames("norm", "T");
614 f->SetLineWidth(1);
615 return f;
616
617
618}
619
f8416b14 620TF1 * GetLevy(Double_t temp, Double_t n, Double_t norm, const char * name) {
a23f7c97 621
622 char formula[500];
623 Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
624
625
626 sprintf(formula,"(x/sqrt(x*x+[3]*[3]))*x*( [0]*([1]-1)*([1]-2) )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) ) ) * ( 1 + (sqrt([3]*[3]+x*x) -[3])/([1]*[2]) )^(-[1])");
627 TF1* f=new TF1(name,formula,0,10);
628 f->SetParameters(norm, n, temp,mass);
629 f->SetParLimits(2, 0.01, 10);
630 f->SetParNames("norm (dN/dy)", "n", "T", "mass");
631 f->FixParameter(3,mass);
632 f->SetLineWidth(1);
633 return f;
634
635
636}
637
638
639TH1D * GetCumulativeHisto (TH1 * h) {
640 // Returns a cumulative histogram
9441cdaa 641 // TODO: put this method in a tools class
f8416b14 642 TH1D * hInt = (TH1D*) h->Clone(TString(h->GetName())+"_Int");
a23f7c97 643 hInt->Reset();
644 Float_t integral = h->Integral(-1,-1); // Consider under/over flow!
645 Int_t nbin = h->GetNbinsX();
646 for(Int_t ibin = 1; ibin <= nbin; ibin++){
647 Double_t content = h->Integral(-1,ibin) / integral;
648 hInt->SetBinContent(ibin, content);
649 }
650 return hInt;
651}
652
653TH1D * GetRatioIntegratedFractions (TH1 * hNum, TH1 * hDenum) {
654 // Returns the the ratio of integrated histograms
9441cdaa 655 // TODO: put this method in a tools class
f8416b14 656 TH1D * hRatio = (TH1D*) hNum->Clone(TString(hNum->GetName())+hDenum->GetName()+"_RatioIntegrated");
a23f7c97 657 hRatio->Reset();
658 Int_t nbin = hNum->GetNbinsX();
659 for(Int_t ibin = 1; ibin <= nbin; ibin++){
660 Double_t content = hNum->Integral(-1,ibin) / hDenum->Integral(-1,ibin);// consider underflow
661 hRatio->SetBinContent(ibin, content);
662 }
663 return hRatio;
664}
665
666void ShowAcceptanceInVzSlices() {
667 TCanvas * cvz = new TCanvas("cvz","acc #times eff vs vz");
e0376287 668 for(Int_t ivz = -10; ivz < -6; ivz+=2){
669 ivz=0;//FIXME
670 Bool_t first = kTRUE;
72491d7c 671 TH1D * hMCPtPri = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim , etaMin,etaMax,ivz,ivz+2);
672 TH1D * hMCPtGen = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen, etaMin,etaMax,ivz,ivz+2);
673 TH1D * hMCPtPriD = hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim , etaMin,etaMax,ivz,ivz+2);
674 TH1D * hMCPtGenD = hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen, etaMin,etaMax,ivz,ivz+2);
a23f7c97 675 // hEff= hMCPtGen;
e0376287 676 TH1D * hEff = (TH1D*)hMCPtPri->Clone(Form("hEff_vz_%d_%d",ivz,ivz+2));
a23f7c97 677 hEff->Divide(hMCPtPri,hMCPtGen,1,1,"B");
e0376287 678 TH1D * hEffD = (TH1D*)hMCPtPriD->Clone(Form("hEffD_vz_%d_%d",ivz,ivz+2));
679 hEffD->Divide(hMCPtPriD,hMCPtGenD,1,1,"B");
680 hEffD->SetLineColor(kRed);
a23f7c97 681 cout << "ivz " << ivz << endl;
e0376287 682 if(first) {
683 first=kFALSE;
a23f7c97 684 cout << "First" << endl;
685 hEff->Draw();
e0376287 686 hEffD->Draw("same");
a23f7c97 687 // hMCPtGen->Draw();
688 // hMCPtPri->Draw("same");
689 }
690 else {
691 cout << "Same" << endl;
692 hEff->Draw("same");
e0376287 693 hEffD->Draw("same");
a23f7c97 694 // hMCPtGen->Draw("");
695 // hMCPtPri->Draw("same");
696 }
697 cvz->Update();
698 // cvz->WaitPrimitive();
699 }
700
701 //hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim )->Draw();
e0376287 702 // hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoRec )->Draw("");
703 // hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoGen )->Draw("same");
a23f7c97 704
705
706}
f8416b14 707
708TF1 * GetFunctionHistoSum() {
709
bcc49144 710 TF1 * f = new TF1 ("fHistoSum",HistoSum,-1,1,kHistoFitCompoments);
abd808b9 711 // f->SetParNames("Primaries+Fakes", "Sec. Weak decays","Sec. Material");
712 f->SetParNames("Primaries+Sec. Material", "Sec. Weak decays");
f8416b14 713 f->SetNpx(1000);
714 return f;
715
716}
717
718static Double_t HistoSum(const double * x, const double* p){
719 // This function uses a global array of histograms, rescaled by some
720 // parameters, to define a function
721 // The array is called gHistoCompoments
722 // The size of the arrays is given by the global constant kHistoFitCompoments
723
724 Double_t xx = x[0];
725 Double_t value = 0;
726 for(Int_t icomp = 0; icomp < kHistoFitCompoments; icomp++){
727 // value += gHistoCompoments[icomp]-Interpolate(xx) * p[icomp];
728 Int_t ibin = gHistoCompoments[icomp]->FindBin(xx);
729 value += gHistoCompoments[icomp]->GetBinContent(ibin) * p[icomp];
730 }
731
732 return value;
733
734}
e0376287 735
736void PrintCanvas(TCanvas* c,const TString formats) {
737 // print a canvas in every of the given comma-separated formats
738 // ensure the canvas is updated
739 if(!doPrint) return;
740 c->Update();
741 gSystem->ProcessEvents();
742 TObjArray * arr = formats.Tokenize(",");
743 TIterator * iter = arr->MakeIterator();
744 TObjString * element = 0;
745 TString name =c ->GetName();
746 name.ReplaceAll(" ","_");
747 name.ReplaceAll("+","Plus");
748 name.ReplaceAll("-","");
749 while ((element = (TObjString*) iter->Next())) {
750 c->Print(name+ "."+element->GetString());
751 }
752}
9441cdaa 753