PWGUD/dNdPt -> PWGLF/SPECTRA/ChargedHadrons/dNdPt
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / ChargedHadrons / dNdPt / macros / ApplyCorrections_PbPb.C
CommitLineData
62e3b4b6 1Double_t ApplyCorrections_PbPb(const char* datafile, const char* datatask, const char* corrfile, const char* idstring ,const char* outfile, const char* gifdir = 0)
2//Double_t ApplyCorrections_PbPb()
3{
4
5// tmp setting
6// const char* datafile = "/lustre/alice/train/V006.PbPb/2011-03-15_0009.5917/mergedPeriods/PbPb/2.76ATeV/LHC10h.pass1/mknichel_dNdPtPbPb_TPCITS_VZERO1.root";
7// const char* datatask = "mknichel_dNdPtPbPb_TPCITS_VZERO";
8// const char* corrfile = "/u/mknichel/alice/dNdPt_PbPb/2011-03-15/corrMatr_LHC10h8_TPCITS.root";
9// const char* idstring = "c0";
10// const char* outfile = "/u/mknichel/alice/dNdPt_PbPb/2011-03-15/finalSpectra_LHC10h.pass1_TPCITS.root";
11// const char* gifdir = "/u/mknichel/alice/dNdPt_PbPb/2011-03-15/LHC10h.pass1";
12
13
14Int_t c_first = 1;
15Int_t c_last = 11;
16
17TString id = TString(idstring);
18if ( id.Contains("c0-5") ) { c_first = c_last = 1; }
19if ( id.Contains("c5-10") ) { c_first = c_last = 2; }
20if ( id.Contains("c10-20") ) { c_first = c_last = 3; }
21if ( id.Contains("c20-30") ) { c_first = c_last = 4; }
22if ( id.Contains("c30-40") ) { c_first = c_last = 5; }
23if ( id.Contains("c40-50") ) { c_first = c_last = 6; }
24if ( id.Contains("c50-60") ) { c_first = c_last = 7; }
25if ( id.Contains("c60-70") ) { c_first = c_last = 8; }
26if ( id.Contains("c70-80") ) { c_first = c_last = 9; }
27if ( id.Contains("c80-90") ) { c_first = c_last = 10; }
28if ( id.Contains("c90-100") ) { c_first = c_last = 11; }
29
30 // settings vor zVertex cut (event and track level)
31 Double_t zVert = 10.0;
32
33 // setting on eta cut (track level)
34 Double_t eta = 0.8;
35
36 //load required libraries
37 //load required libraries
2bfe5463 38 gSystem->AddIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/ -I$ALICE_ROOT/include -I$ALICE_ROOT/STEER -I$ALICE_ROOT/ANALYSIS -I$ALICE_ROOT/PWG0 -I$ALICE_ROOT/PWGPP -I$ALICE_ROOT/PWG2 -I$ALICE_ROOT/PWG3 -I$ALICE_ROOT/PWG3/vertexingHF -I$ALICE_ROOT/PWG4 -I$ALICE_ROOT/CORRFW -I$ALICE_ROOT/TPC -I$ALICE_ROOT/TRD -I$ALICE_ROOT/PWG3/muon -I$ALICE_ROOT/JETAN -I$ALICE_ROOT/ANALYSIS/Tender");
62e3b4b6 39
40 gSystem->Load("libCore");
41 gSystem->Load("libPhysics");
42 gSystem->Load("libMinuit");
43 gSystem->Load("libGui");
44 gSystem->Load("libXMLParser");
45
46 gSystem->Load("libGeom");
47 gSystem->Load("libVMC");
48
49 gSystem->Load("libNet");
50
51 gSystem->Load("libSTEERBase");
52 gSystem->Load("libESD");
53 gSystem->Load("libCDB");
54 gSystem->Load("libRAWDatabase");
55 gSystem->Load("libRAWDatarec");
56 gSystem->Load("libANALYSIS");
57
58
59
60 gSystem->Load("libANALYSIS.so");
61 gSystem->Load("libANALYSISalice.so");
62 gSystem->Load("libTENDER.so");
63 gSystem->Load("libCORRFW.so");
64 gSystem->Load("libPWG0base.so");
65 gSystem->Load("libPWG0dep");
66 gSystem->Load("libPWG0selectors.so");
67
68
69 // make plots nicer
70 gROOT->SetStyle("Plain");
71 gStyle->SetPalette(1);
72
73 // array for all histograms to be saved
74 TObjArray* Hists = new TObjArray();
75
76 // open file with correction matrices
77 TFile *fcorr = TFile::Open(corrfile,"READ");
78 if (!fcorr) return -2;
79
80 // load data
81 TFile* fdata = TFile::Open(datafile,"READ");
82 if (!fdata) return -1;
83 TList* ldata = dynamic_cast<TList*>(fdata->Get(datatask));
84 if (!ldata) return -1;
85 AlidNdPtAnalysisPbPb *obj = dynamic_cast<AlidNdPtAnalysisPbPb*>(ldata->FindObject("dNdPtAnalysisPbPb"));
86 if (!obj) return -1;
87
88 //Event statistics
89 THnSparse *fRecEventHist2 = obj->GetRecEventHist2(); //reconstructed events
90 fRecEventHist2->GetAxis(2)->SetRange(c_first,c_last); // select centrality
91 TH2D* h2RecEvent2All = (TH2D*) fRecEventHist2->Projection(0,1)->Clone("h2RecEvent2All");
92 fRecEventHist2->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
93 TH2D* h2RecEvent2 = (TH2D*) fRecEventHist2->Projection(0,1)->Clone("h2RecEvent2");
94
95 Double_t ReconstructedEvents = h2RecEvent2->Integral();
96 Double_t ReconstructedEventsAll = h2RecEvent2All->Integral();
97
98 Hists->Add(h2RecEvent2);
99 Hists->Add(h2RecEvent2All);
100
101 printf("=== Number of events from DATA %lf ===\n",ReconstructedEvents);
102 printf("=== Number of events from DATA (before zVertex cut) %lf ===\n",ReconstructedEventsAll);
103
104 TH1D* h1ReconstructedEvents = new TH1D("h1ReconstructedEvents","h1ReconstructedEvents",1,0,1);
105 TH1D* h1ReconstructedEventsAll = new TH1D("h1ReconstructedEventsAll","h1ReconstructedEventsAll",1,0,1);
106
107 h1ReconstructedEvents->Fill(0,ReconstructedEvents);
108 h1ReconstructedEvents->SetEntries(ReconstructedEvents);
109 h1ReconstructedEventsAll->Fill(0,ReconstructedEventsAll);
110 h1ReconstructedEventsAll->SetEntries(ReconstructedEventsAll);
111
112 Hists->Add(h1ReconstructedEvents);
113 Hists->Add(h1ReconstructedEventsAll);
114
115 // retrieve tracks
116 THnSparse* fRecTrackHist2 = obj->GetRecTrackHist2(2); //after all cuts (2)
117 fRecTrackHist2->GetAxis(3)->SetRange(c_first,c_last); // select centrality
118 fRecTrackHist2->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01); //zVertex
119 fRecTrackHist2->GetAxis(2)->SetRangeUser(-eta, eta-0.01); // eta
120
121 TH1D* h1RecTrackHist2_zv = fRecTrackHist2->Projection(0)->Clone("h1RecTrackHist2_zv");
122 TH1D* h1RecTrackHist2_pt = fRecTrackHist2->Projection(1)->Clone("h1RecTrackHist2_pt");
123 TH1D* h1RecTrackHist2_eta = fRecTrackHist2->Projection(2)->Clone("h1RecTrackHist2_eta");
124
125 Hists->Add(h1RecTrackHist2_zv);
126 Hists->Add(h1RecTrackHist2_pt);
127 Hists->Add(h1RecTrackHist2_eta);
128
129 // retrieve correction matrices for tracking efficiency (note different binning!)
130 TH1D* h1TrackCorr_pt = (TH1D*)fcorr->Get("h1TrackCorr_pt");
131 TH1D* h1TrackCorr_eta = (TH1D*)fcorr->Get("h1TrackCorr_eta");
132
133 // retrieve correction matrices for secondaries (note different binning!)
134 TH1D* h1SecCorr_pt = (TH1D*)fcorr->Get("h1SecCorr_pt");
135 TH1D* h1SecCorr_eta = (TH1D*)fcorr->Get("h1SecCorr_eta");
136
137 // create corrected spectra (as clone of raw data)
138 TH1D* h1Corrected_pt = h1RecTrackHist2_pt->Clone("h1Corrected_pt");
139 TH1D* h1Corrected_eta = h1RecTrackHist2_eta->Clone("h1Corrected_eta");
140
141 // secondaries correction for pt spectrum
142 for (int i=1; i <= h1Corrected_pt->GetNbinsX() ; i++) {
143 Double_t pt = h1Corrected_pt->GetBinCenter(i);
144 Double_t val = h1Corrected_pt->GetBinContent(i);
145 Double_t err = h1Corrected_pt->GetBinError(i);
146 if (pt >= 50) { pt = 49.5; } // above 50 GeV corr matr have low statistics
147 Double_t secCorr = h1SecCorr_pt->GetBinContent(h1SecCorr_pt->FindBin(pt));
148 Double_t secCorrErr = h1SecCorr_pt->GetBinError(h1SecCorr_pt->FindBin(pt));
149 Double_t effCorr = h1TrackCorr_pt->GetBinContent(h1TrackCorr_pt->FindBin(pt));
150 Double_t effCorrErr = h1TrackCorr_pt->GetBinError(h1TrackCorr_pt->FindBin(pt));
151 Double_t corr = effCorr*secCorr;
152 Double_t corrErr = effCorr*secCorrErr + secCorr*effCorrErr; // errors are correlated
153 Double_t cval = val*corr;
154 Double_t cerr = TMath::Sqrt(val*val*corrErr*corrErr + err*err*corr*corr);
155 h1Corrected_pt->SetBinContent(i,cval);
156 h1Corrected_pt->SetBinError(i,cerr);
157 }
158
159
160 // for eta the correction is simpler because of same binning
161 h1Corrected_eta->Multiply(h1SecCorr_eta);
162 h1Corrected_eta->Multiply(h1TrackCorr_eta);
163
164 Hists->Add(h1Corrected_pt);
165 Hists->Add(h1Corrected_eta);
166
167 // create final spectra (as clone of corrected data)
168 TH1D* dNdPt = h1Corrected_pt->Clone("dNdPt");
169 TH1D* dNdEta = h1Corrected_eta->Clone("dNdEta");
170
171 // also uncorrected spectra (as clone of raw data)
172 TH1D* dNdPt_raw = h1RecTrackHist2_pt->Clone("dNdPt_raw");
173 TH1D* dNdEta_raw = h1RecTrackHist2_eta->Clone("dNdEta_raw");
174
175 TH1D* dNdPt_nores = h1Corrected_pt->Clone("dNdPt_nores");
176
177
178//TF1 *fperi = new TF1("fperi","1.00343-0.000608425*x-6.7038e-05*x*x",5.,40.);
179//TF1 *fcent = new TF1("cent","1.01074e+00-1.98127e-03*x-1.19903e-04*x*x",5.,40.);
180TFile* fptcorr = TFile::Open("ptcorrPbPb_150511.root");
181
182TF1 * fun = 0;
183
184 TF1 *fcent = new TF1("fcent","1.01074e+00-1.98127e-03*x-1.19903e-04*x*x",5.,50.);
185 TF1 *fperi = new TF1("fperi","1.00343-0.000608425*x-6.7038e-05*x*x",5.,50.);
186Double_t downscale = 1.;
187 if (c_first != c_last) {
188 cout << "++++++++++++++++++++++++++++++++++++++++" << endl;
189 cout << "WARNING: pt resolution correction error!" << endl;
190 cout << " (works only for single centraliy bins) " << endl;
191 cout << "++++++++++++++++++++++++++++++++++++++++" << endl;
192 }
193 if (c_first == 1) {
194 fun = (TF1*) fptcorr->Get("ptcorr_c0");
195 downscale = 0.8;
196 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
197 cout << "0-5% central pt-resolution correction used!" << endl;
198 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
199 } else if (c_first == 2) {
200 downscale = 0.8;
201 fun =(TF1*) fptcorr->Get("ptcorr_c5");
202 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
203 cout << "5-10% central pt-resolution correction used!" << endl;
204 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
205 } else if (c_first == 3) {
206 downscale = 0.8;
207 fun =(TF1*) fptcorr->Get("ptcorr_c10");
208 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
209 cout << "10-20% central pt-resolution correction used!" << endl;
210 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
211 } else if (c_first == 4) {
212 downscale = 0.9;
213 fun =(TF1*) fptcorr->Get("ptcorr_c20");
214 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
215 cout << "20-30% central pt-resolution correction used!" << endl;
216 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
217 } else if (c_first == 5) {
218 fun =(TF1*) fptcorr->Get("ptcorr_c30");
219 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
220 cout << "30-40% central pt-resolution correction used!" << endl;
221 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
222 } else if (c_first == 6) {
223 fun =(TF1*) fptcorr->Get("ptcorr_c40");
224 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
225 cout << "40-50% central pt-resolution correction used!" << endl;
226 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
227 } else if (c_first == 7) {
228 fun =(TF1*) fptcorr->Get("ptcorr_c50");
229 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
230 cout << "50-60% central pt-resolution correction used!" << endl;
231 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
232 } else if (c_first == 8) {
233 fun =(TF1*) fptcorr->Get("ptcorr_c60");
234 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
235 cout << "60-70% central pt-resolution correction used!" << endl;
236 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
237 } else {
238 fun =(TF1*) fptcorr->Get("ptcorr_c70");
239 cout << "+++++++++++++++++++++++++++++++++++++++++++++" << endl;
240 cout << "70-80% central pt-resolution correction used!" << endl;
241 cout << "+++++++++++++++++++++++++++++++++++++++++++++" << endl;
242 }
243
244
245 // normalization and finalization
246 // 1/N_evt 1/(2 pi pt) 1/width 1/etarange
247
248 for (int i=1; i <= dNdPt->GetNbinsX() ;i++) {
249 Double_t pt = dNdPt->GetBinCenter(i);
250 Double_t width = dNdPt->GetBinWidth(i);
251 Double_t val = dNdPt->GetBinContent(i);
252 Double_t err = dNdPt->GetBinError(i);
253 Double_t corrPtResol = 1.0;
254 corrPtResol = 1.-((1.-fun->Eval(pt))*downscale );
255 if (pt < 10.) { corrPtResol = 1.0; }
256 Double_t cval = (val * corrPtResol)/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * pt);
257 Double_t cerr = (err * corrPtResol)/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * pt);
258 dNdPt->SetBinContent(i,cval);
259 dNdPt->SetBinError(i,cerr);
260 }
261 // for dndeta again simpler
262 dNdEta->Scale(1,"width");
263 dNdEta->Scale(1./ReconstructedEvents);
264
265 // normalization and finalization without resolution correction
266 // 1/N_evt 1/(2 pi pt) 1/width 1/etarange
267 for (int i=1; i <= dNdPt_nores->GetNbinsX() ;i++) {
268 Double_t pt = dNdPt_nores->GetBinCenter(i);
269 Double_t width = dNdPt_nores->GetBinWidth(i);
270 Double_t val = dNdPt_nores->GetBinContent(i);
271 Double_t err = dNdPt_nores->GetBinError(i);
272 Double_t cval = (val)/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * pt);
273 Double_t cerr = (err)/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * pt);
274 dNdPt_nores->SetBinContent(i,cval);
275 dNdPt_nores->SetBinError(i,cerr);
276 }
277
278 // normalization and finalization
279 // 1/N_evt 1/(2 pi pt) 1/width 1/etarange
280 for (int i=1; i <= dNdPt_raw->GetNbinsX() ;i++) {
281 Double_t pt = dNdPt_raw->GetBinCenter(i);
282 Double_t width = dNdPt_raw->GetBinWidth(i);
283 Double_t val = dNdPt_raw->GetBinContent(i);
284 Double_t err = dNdPt_raw->GetBinError(i);
285 Double_t cval = val/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * pt);
286 Double_t cerr = err/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * pt);
287 dNdPt_raw->SetBinContent(i,cval);
288 dNdPt_raw->SetBinError(i,cerr);
289 }
290 // for dndeta again simpler
291 dNdEta_raw->Scale(1,"width");
292 dNdEta_raw->Scale(1./ReconstructedEvents);
293
294 dNdEta->SetMarkerStyle(21);
295 dNdPt->SetMarkerStyle(21);
296 dNdPt->SetTitle("; p_{T} (GeV/c) ; 1/N_{evt} 1/(2#pi p_{T}) (d^{2}N_{ch})/(d#eta dp_{T}) (GeV/c)^{-2}");
297 dNdEta->SetTitle("; #eta ; 1/N_{evt} (d^{2}N_{ch})/(d#eta)");
298
299 dNdPt_raw->SetTitle("; p_{T} (GeV/c) ; 1/N_{evt} 1/(2#pi p_{T}) (d^{2}N_{ch})/(d#eta dp_{T}) (GeV/c)^{-2}");
300 dNdEta_raw->SetTitle("; #eta ; 1/N_{evt} (d^{2}N_{ch})/(d#eta)");
301
302 Hists->Add(dNdEta_raw);
303 Hists->Add(dNdPt_raw);
304 Hists->Add(dNdEta);
305 Hists->Add(dNdPt);
306 Hists->Add(dNdPt_nores);
307
308
309 // plot pictures and save to gifdir
310 for (i=0; i < Hists->LastIndex(); i++) {
311 TCanvas* ctmp = PlotHist(Hists->At(i),idstring);
312 if (gifdir && ctmp) {
313 TString gif(gifdir);
314 gif += '/';
315 gif += ctmp->GetName();
316 gif += ".gif";
317 ctmp->SaveAs(gif.Data(),"gif");
318 delete ctmp;
319 }
320 }
321
322 // save all correction matrices and control histograms to file
323 if (!outfile) { return; }
324 TFile *out = TFile::Open(outfile,"RECREATE");
325 Hists->Write();
326 out->Close();
327
328 return ReconstructedEvents;
329
330}
331
332//___________________________________________________________________________
333TCanvas* PlotHist(TObject* hobj, const char* label=0)
334{
335 TH1* h = dynamic_cast<TH1*>(hobj);
336 if (!h) return 0;
337 if (h->GetDimension() > 2) return 0;
338 h->SetStats(0);
339 if ( TString(h->GetName()).Contains("Events")) { h->SetStats(1); }
340 TString t(label);
341 if (label) t += "_";
342 t += h->GetName();
343 h->SetTitle(t.Data());
344 TCanvas* c = new TCanvas(t.Data(),t.Data());
345 if (h->GetDimension() >= 1) {
346 TString xlabel(h->GetXaxis()->GetTitle());
347 if (xlabel.Contains("Pt")) { c->SetLogx(); c->SetLogy(); h->GetXaxis()->SetRangeUser(0.1 , 100.); }
348 if (xlabel.Contains("p_{T}")) { c->SetLogx(); c->SetLogy(); h->GetXaxis()->SetRangeUser(0.1 , 100.); }
349 }
350 if (h->GetDimension() == 2) {
351 TString ylabel(h->GetYaxis()->GetTitle());
352 if (ylabel.Contains("Pt")) { c->SetLogy(); h->GetYaxis()->SetRangeUser(0.1 , 100.); }
353 if (ylabel.Contains("p_{T}")) { c->SetLogy(); h->GetYaxis()->SetRangeUser(0.1 , 100.); }
354 h->Draw("COLZ");
355 }
356 if (h->GetDimension() == 1) {
357 h->Draw();
358 }
359 return c;
360
361}