]> 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
undid removing directories where par files are extracted
[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"
f8416b14 14
15using namespace std;
a23f7c97 16
17AliAnalysisMultPbTrackHistoManager * hManData = 0;
18AliAnalysisMultPbTrackHistoManager * hManCorr = 0;
19TH2D * hEvStatData = 0;
20TH2D * hEvStatCorr = 0;
21
a28a49f6 22const Int_t kHistoFitCompoments = 3;
f8416b14 23TH1D * gHistoCompoments[kHistoFitCompoments];
24
a23f7c97 25void LoadLibs();
26void LoadData(TString dataFolder, TString correctionFolder);
27void SetStyle();
e0376287 28void CheckSecondaries(Double_t & fracWeak, Double_t &fracMaterial);
29void CheckVz();
a23f7c97 30void ShowAcceptanceInVzSlices() ;
f8416b14 31TH1D * GetRatioIntegratedFractions (TH1 * hNum, TH1 * hDenum) ;
e0376287 32TH1D * GetCumulativeHisto (TH1 * h) ;
f8416b14 33static Double_t HistoSum(const double * x, const double* p);
34TF1 * GetFunctionHistoSum() ;
35TF1 * GetMTExp(Float_t norm=68, Float_t t=25) ;
36TF1 * GetHagedorn(Float_t norm=68, Float_t pt0=25, Float_t n=13) ;
37TF1 * GetLevy(Double_t temp=0.1, Double_t n=7, Double_t norm=10, const char * name="fLevy") ;
e0376287 38void PrintCanvas(TCanvas* c,const TString formats) ;
f8416b14 39
e0376287 40// global switches
41Bool_t doPrint=kFALSE;// disable PrintCanvas
f8416b14 42
43
e0376287 44void correct(TString dataFolder = "./output/LHC10g2d_130844_V0M_bin_10/", TString correctionFolder = "./output/LHC10g2a_130844_V0M_bin_10/") {
a23f7c97 45
46 // Load stuff and set some styles
47 LoadLibs();
48 LoadData(dataFolder,correctionFolder);
49 SetStyle();
e0376287 50 ShowAcceptanceInVzSlices();
51 return;
a23f7c97 52
53 // TODO add some cool printout for cuts and centrality selection
54
e0376287 55 CheckVz();
56
57 Double_t fractionWeak = 1, fractionMaterial=1;
58 // CheckSecondaries(fractionWeak, fractionMaterial);
59 cout << "Rescaling secondaries correction, weak: " << fractionWeak << ", material: " << fractionMaterial <<endl;
a23f7c97 60
61
62 // Some shorthands
63 // FIXME: Gen should be projected including overflow in z?
e0376287 64 TH1D * hDataPt = (TH1D*) hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRec, -0.5,0.5,-10,10)->Clone("hDataPt");
65 TH1D * hMCPtGen = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen, -0.5,0.5,-10,10);
66 TH1D * hMCPtRec = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRec, -0.5,0.5,-10,10);
67 TH1D * hMCPtPri = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim, -0.5,0.5,-10,10);
68 TH1D * hMCPtSeM = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat, -0.5,0.5,-10,10);
69 TH1D * hMCPtSeW = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak, -0.5,0.5,-10,10);
70 TH1D * hMCPtFak = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecFake, -0.5,0.5,-10,10);
a23f7c97 71
72 TCanvas * cdata = new TCanvas ("cData", "Data");
e0376287 73 cdata->SetLogy();
a23f7c97 74 hDataPt->Draw();
75 // hMCPtRec->Draw("same");
76
77 TCanvas * cMC = new TCanvas ("cMC", "Monte Carlo");
e0376287 78 cMC->SetLogy();
f8416b14 79 cMC->cd();
a23f7c97 80 hMCPtGen ->Draw();
81 hMCPtRec ->Draw("same");
82 hMCPtPri ->Draw("same");
83 hMCPtSeM ->Draw("same");
84 hMCPtSeW ->Draw("same");
85 hMCPtFak ->Draw("same");
86
87 cout << "Fake/All Rec = " << hMCPtFak->Integral()/hMCPtRec->Integral() << endl;
88 cout << "SM/All Rec = " << hMCPtSeM->Integral()/hMCPtRec->Integral() << endl;
89 cout << "SW/All Rec = " << hMCPtSeW->Integral()/hMCPtRec->Integral() << endl;
90
91
92 // Compute efficiency and subtract secondaries and fakes after rescaling to data
93 // PRIM_DATA = ALL_DATA - SEC_MC/ALL_MC*ALL_DATA - FAK_MC/ALL_MC*ALL_DATA
94 // TRUE_DATA = PRIM_DATA * GEN_MC/PRIM_MC
95
96 TH1D * hEffPt = (TH1D*) hMCPtPri->Clone("hEffPt");
97 hEffPt->Divide(hMCPtPri,hMCPtGen,1,1,"B");
98
99 TH1D * hCorSeM = (TH1D*) hMCPtSeM->Clone("hEffPt");
100 hCorSeM->Divide(hMCPtSeM,hMCPtRec,1,1,"B");
e0376287 101 hCorSeM->Scale(fractionMaterial);// rescale material correction
a23f7c97 102 hCorSeM->Multiply(hDataPt);
103
104 TH1D * hCorSeW = (TH1D*) hMCPtSeW->Clone("hEffPt");
105 hCorSeW->Divide(hMCPtSeW,hMCPtRec,1,1,"B");
106 hCorSeW->Scale(fractionWeak);// rescale weak correction
107 hCorSeW->Multiply(hDataPt);
108
109 TH1D * hCorFak = (TH1D*) hMCPtFak->Clone("hEffPt");
110 hCorFak->Divide(hMCPtFak,hMCPtRec,1,1,"B");
111 hCorFak->Multiply(hDataPt);
112
113 TH1D * hCorrected = (TH1D*) hDataPt->Clone("hCorrected");
114 hCorrected->Add(hCorSeM,-1);
115 hCorrected->Add(hCorSeW,-1);
116 hCorrected->Add(hCorFak,-1);
117 hCorrected->Divide(hEffPt);
118 hCorrected->SetMarkerStyle(kOpenStar);
119
120 cdata->cd();
121 hDataPt->Draw();
122 hCorrected->SetLineColor(kBlack);
123 hCorrected->SetMarkerColor(kBlack);
124 hCorrected->Draw("same");
125 hMCPtGen->DrawClone("same");
126 TF1 * f = GetLevy();
127 hCorrected->Fit(f,"", "same");
128 hCorrected->Fit(f,"IME", "same");
129 cout << "dN/deta = " << f->Integral(0,100) << " +- " << f->IntegralError(0,100) << endl;
130 cout << "Generated dN/deta (correction) = " << hMCPtGen->Integral() << endl;
131 // FIXME: comment this out
132 TH1D * hDataGen = hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen, -0.5,0.5,-10,10);
133 cout << "Generated dN/deta (data) = " << hDataGen->Integral() << endl;
134 hDataGen->Draw("same");
135}
136
e0376287 137void CheckSecondaries(Double_t &fracWeak, Double_t &fracMaterial) {
a23f7c97 138 // Returns the fraction you need to rescale the secondaries from weak decays for.
139
140 // Some shorthands
141 TH1D * hDataDCA = hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec );
f8416b14 142 // TH1D * hMCDCAGen = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen );
a23f7c97 143 TH1D * hMCDCARec = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec );
144 TH1D * hMCDCAPri = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim );
145 TH1D * hMCDCASW = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak );
146 TH1D * hMCDCASM = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat );
147 TH1D * hMCDCAFak = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake );
148
149
150 TCanvas * cCumulative = new TCanvas("cDCAculumative","DCA cumulative distributions");
f8416b14 151 cCumulative->cd();
a23f7c97 152 GetCumulativeHisto(hMCDCAPri )->Draw();
153 GetRatioIntegratedFractions(hMCDCASW, hMCDCARec )->Draw("same");
154 GetRatioIntegratedFractions(hMCDCASM, hMCDCARec )->Draw("same");
155 GetRatioIntegratedFractions(hMCDCAPri,hMCDCARec )->Draw("same");
156
157
158 TCanvas * c1 = new TCanvas("cDCAFit", "Fit to the DCA distributions");
159 c1->SetLogy();
160 // Draw all
e0376287 161 // hDataDCA->Draw();
162 // hMCDCARec ->Draw("same");
163 // hMCDCAPri ->Draw("same");
164 // hMCDCASW ->Draw("same");
165 // hMCDCASM ->Draw("same");
166 // hMCDCAFak ->Draw("same");
167 // return;
a23f7c97 168
e0376287 169 // Fit the DCA distribution. Uses a TF1 made by summing histograms
a23f7c97 170 TH1D * hMCPrimSMFak = (TH1D*) hMCDCAPri->Clone("hMCPrimSMFak");
a28a49f6 171 // hMCPrimSMFak->Add(hMCDCASM);
a23f7c97 172 hMCPrimSMFak->Add(hMCDCAFak);
f8416b14 173
e0376287 174 // Set the components which are used in HistoSum, the static
175 // function for GetFunctionHistoSum
f8416b14 176 gHistoCompoments[0] = (TH1D*) hMCPrimSMFak->Clone();
177 gHistoCompoments[1] = (TH1D*) hMCDCASW->Clone();
a28a49f6 178 gHistoCompoments[2] = (TH1D*) hMCDCASM->Clone();
f8416b14 179 TF1 * fHistos = GetFunctionHistoSum();
e0376287 180 fHistos->SetParameters(1,1,1);
181 fHistos->SetLineColor(kRed);
182 // Fit!
a28a49f6 183 hDataDCA->Fit(fHistos,"","",0,200);
e0376287 184 // Rescale the components and draw to see how it looks like
f8416b14 185 hMCPrimSMFak->Scale(fHistos->GetParameter(0));
186 hMCDCASW ->Scale(fHistos->GetParameter(1));
a28a49f6 187 hMCDCASM ->Scale(fHistos->GetParameter(2));
f8416b14 188 hMCPrimSMFak->Draw("same");
189 hMCDCASW ->Draw("same");
a28a49f6 190 hMCDCASM ->Draw("same");
e0376287 191 // compute scaling factors
192 fracWeak = hMCDCASW->Integral()/(hMCPrimSMFak->Integral()+hMCDCASW->Integral()+hMCDCASM->Integral());
193 fracMaterial = hMCDCASM->Integral()/(hMCPrimSMFak->Integral()+hMCDCASW->Integral()+hMCDCASM->Integral());
f8416b14 194
a23f7c97 195
196}
197
e0376287 198void CheckVz() {
199 // compares the Vz distribution in data and in MC
200 TCanvas * c1 = new TCanvas("cVz", "Vertex Z distribution");
201 c1->cd();
202 TH1D * hData = hManData->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoRec );
203 TH1D * hCorr = hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoRec );
204 hCorr->Draw("");
205 hData->Draw("same");
206
207}
208
a23f7c97 209void LoadLibs() {
210
211 gSystem->Load("libVMC");
212 gSystem->Load("libTree");
213 gSystem->Load("libSTEERBase");
214 gSystem->Load("libESD");
215 gSystem->Load("libAOD");
216 gSystem->Load("libANALYSIS");
217 gSystem->Load("libANALYSISalice");
218 gSystem->Load("libCORRFW");
219 gSystem->Load("libMinuit");
220 gSystem->Load("libPWG2Spectra");
221 gSystem->Load("libPWG0base");
222
a28a49f6 223 gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWG0/multPbPb"));
f8416b14 224 gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWG1/background"));
a23f7c97 225 // Load helper classes
226 // TODO: replace this by a list of TOBJStrings
a28a49f6 227 TString taskName("$ALICE_ROOT/PWG0/multPbPb/AliAnalysisTaskMultPbTracks.cxx+");
228 TString histoManName("$ALICE_ROOT/PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.cxx+");
229 TString centrName("$ALICE_ROOT/PWG0/multPbPb/AliAnalysisMultPbCentralitySelector.cxx+");
a23f7c97 230 TString listName("$ALICE_ROOT/PWG1/background/AliHistoListWrapper.cxx+");
231
a28a49f6 232 gSystem->ExpandPathName(taskName);
233 gSystem->ExpandPathName(histoManName);
234 gSystem->ExpandPathName(centrName);
235 gSystem->ExpandPathName(listName);
236
a23f7c97 237 Bool_t debug=0;
238 gROOT->LoadMacro(listName +(debug?"+g":""));
239 gROOT->LoadMacro(histoManName+(debug?"+g":""));
240 gROOT->LoadMacro(centrName +(debug?"+g":""));
241 gROOT->LoadMacro(taskName +(debug?"+g":""));
242
243 // Histo fitter
244 gROOT->LoadMacro("/Users/mfloris/Work/ALICE/ANALYSIS/HistoFitter/fcn.cxx+g");
f8416b14 245 gROOT->LoadMacro("/Users/mfloris/Work/ALICE/ANALYSIS/HistoFitter/AliHistoFitter.cxx+g");
a23f7c97 246
247
248}
249
250
251void SetStyle() {
252
253 hManData->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec) ->SetLineColor(kBlack);
254 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen) ->SetLineColor(kRed );
255 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec) ->SetLineColor(kRed );
256 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetLineColor(kGreen);
257 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetLineColor(kBlue );
258 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetLineColor(kBlue );
259 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetLineColor(kCyan );
260
261 hManData->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec) ->SetMarkerColor(kBlack);
262 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen) ->SetMarkerColor(kRed );
263 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec) ->SetMarkerColor(kRed );
264 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerColor(kGreen);
265 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerColor(kBlue );
266 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerColor(kBlue );
267 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerColor(kCyan );
268
269 hManData->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec) ->SetMarkerStyle(kFullCircle);
270 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen) ->SetMarkerStyle(kFullSquare);
271 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec) ->SetMarkerStyle(kOpenSquare);
272 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerStyle(kOpenSquare);
273 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerStyle(kOpenSquare);
274 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerStyle(kOpenCircle);
275 hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerStyle(kOpenSquare);
276
277 hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec) ->SetLineColor(kBlack);
278 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen) ->SetLineColor(kRed );
279 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec) ->SetLineColor(kRed );
280 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetLineColor(kGreen);
281 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetLineColor(kBlue );
a28a49f6 282 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetLineColor(kBlue-7 );
a23f7c97 283 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetLineColor(kCyan );
284
285 hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec) ->SetMarkerColor(kBlack);
286 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen) ->SetMarkerColor(kRed );
287 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec) ->SetMarkerColor(kRed );
288 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerColor(kGreen);
289 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerColor(kBlue );
a28a49f6 290 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerColor(kBlue-7 );
a23f7c97 291 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerColor(kCyan );
292
293 hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec) ->SetMarkerStyle(kFullCircle);
294 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen) ->SetMarkerStyle(kFullSquare);
295 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec) ->SetMarkerStyle(kOpenSquare);
296 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerStyle(kOpenSquare);
297 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerStyle(kOpenSquare);
298 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerStyle(kOpenSquare);
299 hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerStyle(kOpenSquare);
300
301
302}
303
304void LoadData(TString dataFolder, TString correctionFolder){
305
306 // Get histo manager for data and MC + stat histos
307 TFile * fData = new TFile(dataFolder+"multPbPbtracks.root");
308 TFile * fCorr = new TFile(correctionFolder+"multPbPbtracks.root");
309 TFile * fStatData = new TFile(dataFolder+"event_stat.root");
310 TFile * fStatCorr = new TFile(correctionFolder+"event_stat.root");
311
312 hManData = (AliAnalysisMultPbTrackHistoManager*) fData->Get("histoManager");
313 hManCorr = (AliAnalysisMultPbTrackHistoManager*) fCorr->Get("histoManager");
314 AliESDtrackCuts * cutsData = (AliESDtrackCuts*) fData->Get("AliESDtrackCuts");
315 AliESDtrackCuts * cutsCorr = (AliESDtrackCuts*) fCorr->Get("AliESDtrackCuts");
316 if (cutsData != cutsCorr) {
317 cout << "ERROR: MC and data do not have the same cuts" << endl;
318 // FIXME: exit here
319 }
320 cutsData->Print();
321 hEvStatData = (TH2D*) fStatData->Get("fHistStatistics");
322 hEvStatCorr = (TH2D*) fStatCorr->Get("fHistStatistics");
323
2bbfd8c6 324 AliAnalysisMultPbCentralitySelector * centrData = (AliAnalysisMultPbCentralitySelector*) fData->Get("");
325
a23f7c97 326 // Normalize
327 Int_t irowGoodTrigger = 1;
328 if (hEvStatCorr && hEvStatData) {
329 // 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 330 // hManData->ScaleHistos(hEvStatData->GetBinContent(AliPhysicsSelection::kStatAccepted,irowGoodTrigger));
331 // hManCorr->ScaleHistos(hEvStatCorr->GetBinContent(AliPhysicsSelection::kStatAccepted,irowGoodTrigger));
332 hManData->ScaleHistos(hManData->GetHistoStats()->GetBinContent(2));
333 hManCorr->ScaleHistos(hManCorr->GetHistoStats()->GetBinContent(2));
a23f7c97 334 } else {
335 cout << "WARNING!!! ARBITRARY SCALING" << endl;
336 hManData->ScaleHistos(1000);
337 hManCorr->ScaleHistos(1000);
338 }
339}
340
f8416b14 341TF1 * GetHagedorn(Float_t norm, Float_t pt0, Float_t n) {
a23f7c97 342
343 TF1 * f =0;
344 Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
345
346 f=new TF1("fHagedorn",Form("(x/sqrt(x*x+%f*%f))*x*[0]*( 1 + x/[1] )^(-[2])",mass,mass),0,10);
347 f->SetParameters(norm, pt0, n);
348 f->SetParLimits(1, 0.01, 10);
349 f->SetParNames("norm", "pt0", "n");
350 f->SetLineWidth(1);
351 return f;
352
353
354}
355
f8416b14 356TF1 * GetMTExp(Float_t norm, Float_t t) {
a23f7c97 357
358 TF1 * f =0;
359 Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
360
361 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);
362 f->SetParameters(norm, t);
363 // f->SetParLimits(1, 0.01);
364 f->SetParNames("norm", "T");
365 f->SetLineWidth(1);
366 return f;
367
368
369}
370
f8416b14 371TF1 * GetLevy(Double_t temp, Double_t n, Double_t norm, const char * name) {
a23f7c97 372
373 char formula[500];
374 Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
375
376
377 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])");
378 TF1* f=new TF1(name,formula,0,10);
379 f->SetParameters(norm, n, temp,mass);
380 f->SetParLimits(2, 0.01, 10);
381 f->SetParNames("norm (dN/dy)", "n", "T", "mass");
382 f->FixParameter(3,mass);
383 f->SetLineWidth(1);
384 return f;
385
386
387}
388
389
390TH1D * GetCumulativeHisto (TH1 * h) {
391 // Returns a cumulative histogram
392 // FIXME: put this method in a tools class
f8416b14 393 TH1D * hInt = (TH1D*) h->Clone(TString(h->GetName())+"_Int");
a23f7c97 394 hInt->Reset();
395 Float_t integral = h->Integral(-1,-1); // Consider under/over flow!
396 Int_t nbin = h->GetNbinsX();
397 for(Int_t ibin = 1; ibin <= nbin; ibin++){
398 Double_t content = h->Integral(-1,ibin) / integral;
399 hInt->SetBinContent(ibin, content);
400 }
401 return hInt;
402}
403
404TH1D * GetRatioIntegratedFractions (TH1 * hNum, TH1 * hDenum) {
405 // Returns the the ratio of integrated histograms
406 // FIXME: put this method in a tools class
f8416b14 407 TH1D * hRatio = (TH1D*) hNum->Clone(TString(hNum->GetName())+hDenum->GetName()+"_RatioIntegrated");
a23f7c97 408 hRatio->Reset();
409 Int_t nbin = hNum->GetNbinsX();
410 for(Int_t ibin = 1; ibin <= nbin; ibin++){
411 Double_t content = hNum->Integral(-1,ibin) / hDenum->Integral(-1,ibin);// consider underflow
412 hRatio->SetBinContent(ibin, content);
413 }
414 return hRatio;
415}
416
417void ShowAcceptanceInVzSlices() {
418 TCanvas * cvz = new TCanvas("cvz","acc #times eff vs vz");
e0376287 419 for(Int_t ivz = -10; ivz < -6; ivz+=2){
420 ivz=0;//FIXME
421 Bool_t first = kTRUE;
422 TH1D * hMCPtPri = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim , -0.5,0.5,ivz,ivz+2);
423 TH1D * hMCPtGen = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen, -0.5,0.5,ivz,ivz+2);
424 TH1D * hMCPtPriD = hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim , -0.5,0.5,ivz,ivz+2);
425 TH1D * hMCPtGenD = hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen, -0.5,0.5,ivz,ivz+2);
a23f7c97 426 // hEff= hMCPtGen;
e0376287 427 TH1D * hEff = (TH1D*)hMCPtPri->Clone(Form("hEff_vz_%d_%d",ivz,ivz+2));
a23f7c97 428 hEff->Divide(hMCPtPri,hMCPtGen,1,1,"B");
e0376287 429 TH1D * hEffD = (TH1D*)hMCPtPriD->Clone(Form("hEffD_vz_%d_%d",ivz,ivz+2));
430 hEffD->Divide(hMCPtPriD,hMCPtGenD,1,1,"B");
431 hEffD->SetLineColor(kRed);
a23f7c97 432 cout << "ivz " << ivz << endl;
e0376287 433 if(first) {
434 first=kFALSE;
a23f7c97 435 cout << "First" << endl;
436 hEff->Draw();
e0376287 437 hEffD->Draw("same");
a23f7c97 438 // hMCPtGen->Draw();
439 // hMCPtPri->Draw("same");
440 }
441 else {
442 cout << "Same" << endl;
443 hEff->Draw("same");
e0376287 444 hEffD->Draw("same");
a23f7c97 445 // hMCPtGen->Draw("");
446 // hMCPtPri->Draw("same");
447 }
448 cvz->Update();
449 // cvz->WaitPrimitive();
450 }
451
452 //hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim )->Draw();
e0376287 453 // hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoRec )->Draw("");
454 // hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoGen )->Draw("same");
a23f7c97 455
456
457}
f8416b14 458
459TF1 * GetFunctionHistoSum() {
460
461 TF1 * f = new TF1 ("fHistoSum",HistoSum,0,50,kHistoFitCompoments);
a28a49f6 462 f->SetParNames("Primaries+Fakes", "Sec. Weak decays","Sec. Material");
f8416b14 463 f->SetNpx(1000);
464 return f;
465
466}
467
468static Double_t HistoSum(const double * x, const double* p){
469 // This function uses a global array of histograms, rescaled by some
470 // parameters, to define a function
471 // The array is called gHistoCompoments
472 // The size of the arrays is given by the global constant kHistoFitCompoments
473
474 Double_t xx = x[0];
475 Double_t value = 0;
476 for(Int_t icomp = 0; icomp < kHistoFitCompoments; icomp++){
477 // value += gHistoCompoments[icomp]-Interpolate(xx) * p[icomp];
478 Int_t ibin = gHistoCompoments[icomp]->FindBin(xx);
479 value += gHistoCompoments[icomp]->GetBinContent(ibin) * p[icomp];
480 }
481
482 return value;
483
484}
e0376287 485
486void PrintCanvas(TCanvas* c,const TString formats) {
487 // print a canvas in every of the given comma-separated formats
488 // ensure the canvas is updated
489 if(!doPrint) return;
490 c->Update();
491 gSystem->ProcessEvents();
492 TObjArray * arr = formats.Tokenize(",");
493 TIterator * iter = arr->MakeIterator();
494 TObjString * element = 0;
495 TString name =c ->GetName();
496 name.ReplaceAll(" ","_");
497 name.ReplaceAll("+","Plus");
498 name.ReplaceAll("-","");
499 while ((element = (TObjString*) iter->Next())) {
500 c->Print(name+ "."+element->GetString());
501 }
502}