]>
Commit | Line | Data |
---|---|---|
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 | |
17 | using namespace std; | |
a23f7c97 | 18 | |
19 | AliAnalysisMultPbTrackHistoManager * hManData = 0; | |
20 | AliAnalysisMultPbTrackHistoManager * hManCorr = 0; | |
21 | TH2D * hEvStatData = 0; | |
22 | TH2D * hEvStatCorr = 0; | |
23 | ||
bcc49144 | 24 | const Int_t kHistoFitCompoments = 2; |
f8416b14 | 25 | TH1D * gHistoCompoments[kHistoFitCompoments]; |
26 | ||
aa6151eb | 27 | void LoadLibs( Bool_t debug=0); |
a23f7c97 | 28 | void LoadData(TString dataFolder, TString correctionFolder); |
29 | void SetStyle(); | |
e0376287 | 30 | void CheckSecondaries(Double_t & fracWeak, Double_t &fracMaterial); |
72491d7c | 31 | void CheckSanity(); |
9441cdaa | 32 | void CheckCorrections(); |
a23f7c97 | 33 | void ShowAcceptanceInVzSlices() ; |
f8416b14 | 34 | TH1D * GetRatioIntegratedFractions (TH1 * hNum, TH1 * hDenum) ; |
e0376287 | 35 | TH1D * GetCumulativeHisto (TH1 * h) ; |
f8416b14 | 36 | static Double_t HistoSum(const double * x, const double* p); |
37 | TF1 * GetFunctionHistoSum() ; | |
38 | TF1 * GetMTExp(Float_t norm=68, Float_t t=25) ; | |
39 | TF1 * GetHagedorn(Float_t norm=68, Float_t pt0=25, Float_t n=13) ; | |
40 | TF1 * GetLevy(Double_t temp=0.1, Double_t n=7, Double_t norm=10, const char * name="fLevy") ; | |
e0376287 | 41 | void PrintCanvas(TCanvas* c,const TString formats) ; |
f8416b14 | 42 | |
e0376287 | 43 | // global switches |
44 | Bool_t doPrint=kFALSE;// disable PrintCanvas | |
55e8e56a | 45 | Float_t zmin = -10; |
46 | Float_t zmax = 10; | |
72491d7c | 47 | Float_t etaMin = -0.5; |
48 | Float_t etaMax = 0.5; | |
bfae2924 | 49 | Int_t gCentralityBin = -1; |
f8416b14 | 50 | |
aa6151eb | 51 | #define CORRECT_2D |
f8416b14 | 52 | |
72491d7c | 53 | void 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 | 256 | void 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 |
2b42cee2 | 297 | const Int_t ptBinsFit[] = {3,5,7,9,11,15,19,23,31,-1}; |
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 | 350 | void 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 | ||
424 | void 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 | 446 | void 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")); |
2bfe5463 | 461 | gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWGPP/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+"); | |
2bfe5463 | 469 | TString listName("$ALICE_ROOT/PWGPP/background/AliHistoListWrapper.cxx+"); |
a23f7c97 | 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 | ||
490 | void 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 | ||
543 | void 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 | 590 | TF1 * 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 | 605 | TF1 * 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 | 620 | TF1 * 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 | ||
639 | TH1D * 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 | ||
653 | TH1D * 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 | ||
666 | void 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 | |
708 | TF1 * 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 | ||
718 | static 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 | |
736 | void 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 |