]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGUD/dNdPt/macros/GenerateCorrMatr_PbPb.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Transition PWG0 -> PWGUD
[u/mrichter/AliRoot.git] / PWGUD / dNdPt / macros / GenerateCorrMatr_PbPb.C
CommitLineData
62e3b4b6 1&
2Double_t GenerateCorrMatr_PbPb(const char* mcfile, const char* mctask, const char* idstring ,const char* outfile, const char* gifdir = 0)
3//Double_t GenerateCorrMatr_PbPb_TPCIT()
4{
5
6Int_t c_first = 1;
7Int_t c_last = 11;
8
9TString id = TString(idstring);
10if ( id.Contains("c0-5") ) { c_first = c_last = 1; }
11if ( id.Contains("c5-10") ) { c_first = c_last = 2; }
12if ( id.Contains("c10-20") ) { c_first = c_last = 3; }
13if ( id.Contains("c20-30") ) { c_first = c_last = 4; }
14if ( id.Contains("c30-40") ) { c_first = c_last = 5; }
15if ( id.Contains("c40-50") ) { c_first = c_last = 6; }
16if ( id.Contains("c50-60") ) { c_first = c_last = 7; }
17if ( id.Contains("c60-70") ) { c_first = c_last = 8; }
18if ( id.Contains("c70-80") ) { c_first = c_last = 9; }
19if ( id.Contains("c80-90") ) { c_first = c_last = 10; }
20if ( id.Contains("c90-100") ) { c_first = c_last = 11; }
21
22//tmp setting
23 //const char* mcfile = "/lustre/alice/train/V006.MC_PbPb/2011-03-15_0037.5917/mergedPeriods/MC_PbPb/2.76ATeV/LHC10h8/mknichel_dNdPtPbPb_TPCITS_VZERO1.root";
24
25 //const char* mcfile = "/lustre/alice/train/V006.MC_PbPb/2011-03-13_0236.5847/mergedRuns/MC_PbPb/2.76ATeV/LHC11a7/137161.ana/mknichel_dNdPtPbPb_TPCITS_VZERO1.root";
26 //const char* mctask = "mknichel_dNdPtPbPb_TPCITS_VZERO";
27 //const char* idstring = "c70";
28 //const char* outfile = "/u/mknichel/alice/dNdPt_PbPb/2011-03-15/corrMatr_LHC10h8_TPCITS_c70.root";
29 //const char* gifdir ="/u/mknichel/alice/dNdPt_PbPb/2011-03-15/LHC10h8";
30
31 // settings vor zVertex cut (event and track level)
32 Double_t zVert = 10.0;
33
34 // setting on eta cut (track level)
35 Double_t eta = 0.8;
36
37 // strangeness scaling factor (for secondaries from strange decays)
38 Double_t sscale = 2.0;
39
40 //load required libraries
41 //load required libraries
2bfe5463 42 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 43
44 gSystem->Load("libCore");
45 gSystem->Load("libPhysics");
46 gSystem->Load("libMinuit");
47 gSystem->Load("libGui");
48 gSystem->Load("libXMLParser");
49
50 gSystem->Load("libGeom");
51 gSystem->Load("libVMC");
52
53 gSystem->Load("libNet");
54
55
56 gSystem->Load("libSTEERBase");
57 gSystem->Load("libESD");
58 gSystem->Load("libCDB");
59 gSystem->Load("libRAWDatabase");
60 gSystem->Load("libRAWDatarec");
61 gSystem->Load("libANALYSIS");
62
63
64
65 gSystem->Load("libANALYSIS.so");
66 gSystem->Load("libANALYSISalice.so");
67 gSystem->Load("libTENDER.so");
68 gSystem->Load("libCORRFW.so");
69 gSystem->Load("libPWG0base.so");
70 gSystem->Load("libPWG0dep");
71 gSystem->Load("libPWG0selectors.so");
72
73
74
75
76
77 // make plots nicer
78 gROOT->SetStyle("Plain");
79 gStyle->SetPalette(1);
80
81 // array for all correction matrices
82 TObjArray* CorrMatr = new TObjArray();
83
84 // array for all control histograms
85 TObjArray* ContHist = new TObjArray();
86
87
88 // load mc information
89 TFile* fmc = TFile::Open(mcfile,"READ");
90 if (!fmc) return -1;
91 TList* lmc = dynamic_cast<TList*>(fmc->Get(mctask));
92 if (!lmc) return -1;
93 AlidNdPtAnalysisPbPb *obj = dynamic_cast<AlidNdPtAnalysisPbPb*>(lmc->FindObject("dNdPtAnalysisPbPb"));
94 if (!obj) return -1;
95
96//
97// create rebinned thnsparse
98//
99 //ptbins before rebinning: 73;
100 const Int_t etaNbins = 30;
101 const Int_t zvNbins = 12;
102 const Int_t centralityBins = 11;
103
104 //Double_t binsMult[multNbins+1] = {-0.5, 10000.5 }; // forPbPb --reduced!
105
106 // Set pt binning for corrections
107 const Int_t ptNbins = 30;
108 Double_t binsPt[31] = {0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.4, 1.6, 1.8, 2.0, 3.0, 4.0, 50.0, 100.0 };
109
110 Double_t binsEta[etaNbins+1] = {
111 -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6,
112 -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4,
113 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4,
114 1.5};
115
116 Double_t binsZv[zvNbins+1] = {-30.,-25.,-20.,-15.,-10.,-5.,0.,5.,10.,15.,20.,25.,30.};
117
118 Double_t binsCentrality[centralityBins+1] = {0., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.};
119
120 Int_t binsTrackEventCorrMatrix[4]={zvNbins,ptNbins,etaNbins,centralityBins};
121 //Int_t binsTrackEvent[4]={zvNbins,ptNbins,etaNbins,centralityBins};
122 //Int_t binsTrackPtCorrelationMatrix[4]={ptNbins,ptNbins,etaNbins,centralityBins};
123/*
124 THnSparseF* fGenPrimTrackMatrix = new THnSparseF("fGenPrimTrackMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
125 fGenPrimTrackMatrix->SetBinEdges(0,binsZv);
126 fGenPrimTrackMatrix->SetBinEdges(1,binsPt);
127 fGenPrimTrackMatrix->SetBinEdges(2,binsEta);
128 fGenPrimTrackMatrix->SetBinEdges(3,binsCentrality);
129 fGenPrimTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
130 fGenPrimTrackMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
131 fGenPrimTrackMatrix->GetAxis(2)->SetTitle("mcEta");
132 fGenPrimTrackMatrix->GetAxis(3)->SetTitle("Centrality");
133 fGenPrimTrackMatrix->Sumw2();
134
135 THnSparseF* fRecPrimTrackMatrix = new THnSparseF("fRecPrimTrackMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
136 fRecPrimTrackMatrix->SetBinEdges(0,binsZv);
137 fRecPrimTrackMatrix->SetBinEdges(1,binsPt);
138 fRecPrimTrackMatrix->SetBinEdges(2,binsEta);
139 fRecPrimTrackMatrix->SetBinEdges(3,binsCentrality);
140 fRecPrimTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
141 fRecPrimTrackMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
142 fRecPrimTrackMatrix->GetAxis(2)->SetTitle("mcEta");
143 fRecPrimTrackMatrix->GetAxis(3)->SetTitle("Centrality");
144 fRecPrimTrackMatrix->Sumw2();
145
146 THnSparseF* fRecTrackMatrix = new THnSparseF("fRecTrackMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
147 fRecTrackMatrix->SetBinEdges(0,binsZv);
148 fRecTrackMatrix->SetBinEdges(1,binsPt);
149 fRecTrackMatrix->SetBinEdges(2,binsEta);
150 fRecTrackMatrix->SetBinEdges(3,binsCentrality);
151 fRecTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
152 fRecTrackMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
153 fRecTrackMatrix->GetAxis(2)->SetTitle("mcEta");
154 fRecTrackMatrix->GetAxis(3)->SetTitle("Centrality");
155 fRecTrackMatrix->Sumw2();
156
157 THnSparseF* fRecSecTrackMatrix = new THnSparseF("fRecSecTrackMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
158 fRecSecTrackMatrix->SetBinEdges(0,binsZv);
159 fRecSecTrackMatrix->SetBinEdges(1,binsPt);
160 fRecSecTrackMatrix->SetBinEdges(2,binsEta);
161 fRecSecTrackMatrix->SetBinEdges(3,binsCentrality);
162 fRecSecTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
163 fRecSecTrackMatrix->GetAxis(1)->SetTitle("Pt (GeV/c)");
164 fRecSecTrackMatrix->GetAxis(2)->SetTitle("Eta");
165 fRecSecTrackMatrix->GetAxis(3)->SetTitle("Centrality");
166 fRecSecTrackMatrix->Sumw2();
167*/
168
169 //Event statistics
170 THnSparse *fRecEventMatrix = obj->GetRecEventMatrix(); //all reconstructed events
171 fRecEventMatrix->GetAxis(2)->SetRange(c_first,c_last); // select centrality
172 TH2D* h2RecEventAll = fRecEventMatrix->Projection(0,1)->Clone("h2RecEventAll");
173 fRecEventMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
174 TH2D* h2RecEvent = fRecEventMatrix->Projection(0,1)->Clone("h2RecEvent");
175 Double_t MCReconstructedEvents = h2RecEvent->Integral();
176 Double_t MCReconstructedEventsAll = h2RecEventAll->Integral();
177 ContHist->Add(h2RecEvent);
178 ContHist->Add(h2RecEventAll);
179
180 THnSparse *fTriggerEventMatrix = obj->GetTriggerEventMatrix(); //all triggered events
181 fTriggerEventMatrix->GetAxis(2)->SetRange(c_first,c_last); // select centrality
182 TH2D* h2TriggerEventAll = fTriggerEventMatrix->Projection(0,1)->Clone("h2TriggerEventAll");
183 fTriggerEventMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
184 TH2D* h2TriggerEvent = fTriggerEventMatrix->Projection(0,1)->Clone("h2TriggerEvent");
185 Double_t MCTriggeredEvents = h2TriggerEvent->Integral();
186 Double_t MCTriggeredEventsAll = h2TriggerEventAll->Integral();
187 ContHist->Add(h2TriggerEvent);
188 ContHist->Add(h2TriggerEventAll);
189
190 THnSparse *fGenEventMatrix = obj->GetGenEventMatrix(); //all generated events
191 fGenEventMatrix->GetAxis(2)->SetRange(c_first,c_last); // select centrality
192 TH2D* h2GenEventAll = fGenEventMatrix->Projection(0,1)->Clone("h2GenEventAll");
193 fGenEventMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
194 TH2D* h2GenEvent = fGenEventMatrix->Projection(0,1)->Clone("h2GenEvent");
195 Double_t MCGeneratedEvents = h2GenEvent->Integral();
196 Double_t MCGeneratedEventsAll = h2GenEventAll->Integral();
197 ContHist->Add(h2RecEvent);
198 ContHist->Add(h2RecEventAll);
199
200 printf("=== generated MC events for correction matrices %lf ===\n",MCGeneratedEvents);
201 printf("=== triggered MC events for correction matrices %lf ===\n",MCTriggeredEvents);
202 printf("=== recontructed MC events for correction matrices %lf ===\n",MCReconstructedEvents);
203 printf("\n");
204 printf("=== cut on the zVertex +- %lf ===\n",zVert);
205 printf("=== generated MC events (before zVertex cut) %lf ===\n",MCGeneratedEventsAll);
206 printf("=== triggered MC events (before zVertex cut) %lf ===\n",MCTriggeredEventsAll);
207 printf("=== recontructed MC events (before zVertex cut) %lf ===\n",MCReconstructedEventsAll);
208
209 TH1D* h1MCGeneratedEvents = new TH1D("h1MCGeneratedEvents","h1MCGeneratedEvents",1,0,1);
210 TH1D* h1MCTriggeredEvents = new TH1D("h1MCTriggeredEvents","h1MCTriggeredEvents",1,0,1);
211 TH1D* h1MCReconstructedEvents = new TH1D("h1MCReconstructedEvents","h1MCReconstructedEvents",1,0,1);
212 TH1D* h1MCGeneratedEventsAll = new TH1D("h1MCGeneratedEventsAll","h1MCGeneratedEventsAll",1,0,1);
213 TH1D* h1MCTriggeredEventsAll = new TH1D("h1MCTriggeredEventsAll","h1MCTriggeredEventsAll",1,0,1);
214 TH1D* h1MCReconstructedEventsAll = new TH1D("h1MCReconstructedEventsAll","h1MCReconstructedEventsAll",1,0,1);
215
216 h1MCGeneratedEvents->Fill(0,MCGeneratedEvents);
217 h1MCGeneratedEvents->SetEntries(MCGeneratedEvents);
218 h1MCTriggeredEvents->Fill(0,MCTriggeredEvents);
219 h1MCTriggeredEvents->SetEntries(MCTriggeredEvents);
220 h1MCReconstructedEvents->Fill(0,MCReconstructedEvents);
221 h1MCReconstructedEvents->SetEntries(MCReconstructedEvents);
222 h1MCGeneratedEventsAll->Fill(0,MCGeneratedEventsAll);
223 h1MCGeneratedEventsAll->SetEntries(MCGeneratedEventsAll);
224 h1MCTriggeredEventsAll->Fill(0,MCTriggeredEventsAll);
225 h1MCTriggeredEventsAll->SetEntries(MCTriggeredEventsAll);
226 h1MCReconstructedEventsAll->Fill(0,MCReconstructedEventsAll);
227 h1MCReconstructedEventsAll->SetEntries(MCReconstructedEventsAll);
228
229 ContHist->Add(h1MCGeneratedEvents);
230 ContHist->Add(h1MCTriggeredEvents);
231 ContHist->Add(h1MCReconstructedEvents);
232 ContHist->Add(h1MCGeneratedEventsAll);
233 ContHist->Add(h1MCTriggeredEventsAll);
234 ContHist->Add(h1MCReconstructedEventsAll);
235
236
237 // efficienfy and correction matrices for tigger and vertex efficiency
238 TH2D* h2EventTriggerEffAll = AlidNdPtHelper::GenerateCorrMatrix(h2TriggerEventAll,h2GenEventAll,"h2EventTriggerEffAll");
239 TH2D* h2EventTriggerCorrAll = AlidNdPtHelper::GenerateCorrMatrix(h2GenEventAll,h2TriggerEventAll,"h2EventTriggerCorrAll");
240 TH2D* h2EventTriggerEff = AlidNdPtHelper::GenerateCorrMatrix(h2TriggerEvent,h2GenEvent,"h2EventTriggerEff");
241 TH2D* h2EventTriggerCorr = AlidNdPtHelper::GenerateCorrMatrix(h2GenEvent,h2TriggerEvent,"h2EventTriggerCorr");
242
243 TH2D* h2EventRecEffAll = AlidNdPtHelper::GenerateCorrMatrix(h2RecEventAll,h2TriggerEventAll,"h2EventRecEffAll");
244 TH2D* h2EventRecCorrAll = AlidNdPtHelper::GenerateCorrMatrix(h2TriggerEventAll,h2RecEventAll,"h2EventRecCorrAll");
245 TH2D* h2EventRecEff = AlidNdPtHelper::GenerateCorrMatrix(h2RecEvent,h2TriggerEvent,"h2EventRecEff");
246 TH2D* h2EventRecCorr = AlidNdPtHelper::GenerateCorrMatrix(h2TriggerEvent,h2RecEvent,"h2EventRecCorr");
247
248 TH2D* h2EventEffAll = AlidNdPtHelper::GenerateCorrMatrix(h2RecEventAll,h2GenEventAll,"h2EventEffAll");
249 TH2D* h2EventCorrAll = AlidNdPtHelper::GenerateCorrMatrix(h2GenEventAll,h2RecEventAll,"h2EventCorrAll");
250 TH2D* h2EventEff = AlidNdPtHelper::GenerateCorrMatrix(h2RecEvent,h2GenEvent,"h2EventEff");
251 TH2D* h2EventCorr = AlidNdPtHelper::GenerateCorrMatrix(h2GenEvent,h2RecEvent,"h2EventCorr");
252
253 CorrMatr->Add(h2EventTriggerEffAll);
254 CorrMatr->Add(h2EventTriggerCorrAll);
255 CorrMatr->Add(h2EventTriggerEff);
256 CorrMatr->Add(h2EventTriggerCorr);
257 CorrMatr->Add(h2EventRecEffAll);
258 CorrMatr->Add(h2EventRecCorrAll);
259 CorrMatr->Add(h2EventRecEff);
260 CorrMatr->Add(h2EventRecCorr);
261 CorrMatr->Add(h2EventEffAll);
262 CorrMatr->Add(h2EventCorrAll);
263 CorrMatr->Add(h2EventEff);
264 CorrMatr->Add(h2EventCorr);
265
266
267
268 // all recontructed
269 THnSparse *fRecTrackMatrix = obj->GetRecTrackMatrix(); //all reconstructed tracks
270 fRecTrackMatrix->GetAxis(3)->SetRange(c_first,c_last); // select centrality
271 fRecTrackMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
272 fRecTrackMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
273
274 TH3D* h3RecTrack = fRecTrackMatrix->Projection(0,1,2)->Clone("h3RecTrack");
275 TH2D* h2RecTrack_zv_pt = fRecTrackMatrix->Projection(0,1)->Clone("h2RecTrack_zv_pt");
276 TH2D* h2RecTrack_zv_eta = fRecTrackMatrix->Projection(0,2)->Clone("h2RecTrack_zv_eta");
277 TH2D* h2RecTrack_pt_eta = fRecTrackMatrix->Projection(1,2)->Clone("h2RecTrack_pt_eta");
278 TH1D* h1RecTrack_zv = fRecTrackMatrix->Projection(0)->Clone("h1RecTrack_zv");
279 TH1D* h1RecTrack_pt = fRecTrackMatrix->Projection(1)->Clone("h1RecTrack_pt");
280 TH1D* h1RecTrack_eta = fRecTrackMatrix->Projection(2)->Clone("h1RecTrack_eta");
281 Double_t MCReconstructedTracks = h3RecTrack->Integral();
282
283 h1RecTrack_pt = (TH1D*) h1RecTrack_pt->Rebin(ptNbins,"",binsPt); // rebin hist
284
285 ContHist->Add(h3RecTrack);
286 ContHist->Add(h2RecTrack_zv_pt);
287 ContHist->Add(h2RecTrack_zv_eta);
288 ContHist->Add(h2RecTrack_pt_eta);
289 ContHist->Add(h1RecTrack_zv);
290 ContHist->Add(h1RecTrack_pt);
291 ContHist->Add(h1RecTrack_eta);
292
293 // recontructed primary tracks
294 THnSparse *fRecPrimTrackMatrix = obj->GetRecPrimTrackMatrix();
295 fRecPrimTrackMatrix->GetAxis(3)->SetRange(c_first,c_last); // select centrality
296 THnSparse *fRecTrackMatrixScaled = fRecPrimTrackMatrix->Clone("fRecTrackMatrixScaled"); //used later for secondaries scaling
297 fRecPrimTrackMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
298 fRecPrimTrackMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
299 TH3D* h3RecPrimTrack = fRecPrimTrackMatrix->Projection(0,1,2)->Clone("h3RecPrimTrack");
300 TH2D* h2RecPrimTrack_zv_pt = fRecPrimTrackMatrix->Projection(0,1)->Clone("h2RecPrimTrack_zv_pt");
301 TH2D* h2RecPrimTrack_zv_eta = fRecPrimTrackMatrix->Projection(0,2)->Clone("h2RecPrimTrack_zv_eta");
302 TH2D* h2RecPrimTrack_pt_eta = fRecPrimTrackMatrix->Projection(1,2)->Clone("h2RecPrimTrack_pt_eta");
303 TH1D* h1RecPrimTrack_zv = fRecPrimTrackMatrix->Projection(0)->Clone("h1RecPrimTrack_zv");
304 TH1D* h1RecPrimTrack_pt = fRecPrimTrackMatrix->Projection(1)->Clone("h1RecPrimTrack_pt");
305 TH1D* h1RecPrimTrack_eta = fRecPrimTrackMatrix->Projection(2)->Clone("h1RecPrimTrack_eta");
306 Double_t MCReconstructedPrimTracks = h3RecPrimTrack->Integral();
307
308 h1RecPrimTrack_pt = (TH1D*) h1RecPrimTrack_pt->Rebin(ptNbins,"",binsPt); // rebin hist
309
310 ContHist->Add(h3RecPrimTrack);
311 ContHist->Add(h2RecPrimTrack_zv_pt);
312 ContHist->Add(h2RecPrimTrack_zv_eta);
313 ContHist->Add(h2RecPrimTrack_pt_eta);
314 ContHist->Add(h1RecPrimTrack_zv);
315 ContHist->Add(h1RecPrimTrack_pt);
316 ContHist->Add(h1RecPrimTrack_eta);
317
318 // recontructed secondary tracks
319 THnSparse *fRecSecTrackMatrix = obj->GetRecSecTrackMatrix();
320 fRecSecTrackMatrix->GetAxis(3)->SetRange(c_first,c_last); // select centrality
321 THnSparse *fRecSecTrackMatrixScaled = fRecSecTrackMatrix->Clone("fRecSecTrackMatrixScaled"); //used later for secondaries scaling
322 fRecSecTrackMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
323 fRecSecTrackMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
324 TH3D* h3RecSecTrack = fRecSecTrackMatrix->Projection(0,1,2)->Clone("h3RecSecTrack");
325 TH2D* h2RecSecTrack_zv_pt = fRecSecTrackMatrix->Projection(0,1)->Clone("h2RecSecTrack_zv_pt");
326 TH2D* h2RecSecTrack_zv_eta = fRecSecTrackMatrix->Projection(0,2)->Clone("h2RecSecTrack_zv_eta");
327 TH2D* h2RecSecTrack_pt_eta = fRecSecTrackMatrix->Projection(1,2)->Clone("h2RecSecTrack_pt_eta");
328 TH1D* h1RecSecTrack_zv = fRecSecTrackMatrix->Projection(0)->Clone("h1RecSecTrack_zv");
329 TH1D* h1RecSecTrack_pt = fRecSecTrackMatrix->Projection(1)->Clone("h1RecSecTrack_pt");
330 TH1D* h1RecSecTrack_eta = fRecSecTrackMatrix->Projection(2)->Clone("h1RecSecTrack_eta");
331 Double_t MCReconstructedSecTracks = h3RecSecTrack->Integral();
332
333 h1RecSecTrack_pt = (TH1D*) h1RecSecTrack_pt->Rebin(ptNbins,"",binsPt); // rebin hist
334
335 ContHist->Add(h3RecSecTrack);
336 ContHist->Add(h2RecSecTrack_zv_pt);
337 ContHist->Add(h2RecSecTrack_zv_eta);
338 ContHist->Add(h2RecSecTrack_pt_eta);
339 ContHist->Add(h1RecSecTrack_zv);
340 ContHist->Add(h1RecSecTrack_pt);
341 ContHist->Add(h1RecSecTrack_eta);
342
343 // generated primary tracks
344 THnSparse *fGenPrimTrackMatrix =obj->GetGenPrimTrackMatrix();
345 fGenPrimTrackMatrix->GetAxis(3)->SetRange(c_first,c_last); // select centrality
346 fGenPrimTrackMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
347 fGenPrimTrackMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
348 TH3D* h3GenPrimTrack = fGenPrimTrackMatrix->Projection(0,1,2)->Clone("h3GenPrimTrack");
349 TH2D* h2GenPrimTrack_zv_Pt = fGenPrimTrackMatrix->Projection(0,1)->Clone("h2GenPrimTrack_zv_pt");
350 TH2D* h2GenPrimTrack_zv_eta = fGenPrimTrackMatrix->Projection(0,2)->Clone("h2GenPrimTrack_zv_eta");
351 TH2D* h2GenPrimTrack_pt_eta = fGenPrimTrackMatrix->Projection(1,2)->Clone("h2GenPrimTrack_pt_eta");
352 TH1D* h1GenPrimTrack_zv = fGenPrimTrackMatrix->Projection(0)->Clone("h1GenPrimTrack_zv");
353 TH1D* h1GenPrimTrack_pt = fGenPrimTrackMatrix->Projection(1)->Clone("h1GenPrimTrack_pt");
354 TH1D* h1GenPrimTrack_eta = fGenPrimTrackMatrix->Projection(2)->Clone("h1GenPrimTrack_eta");
355 Double_t MCGeneratedPrimTracks = h3GenPrimTrack->Integral();
356
357
358 // mc truth histogram (for self-consistency check)
359 TH1D* dNdPt_MC = h1GenPrimTrack_pt->Clone("dNdPt_MC");
360 // normalization and finalization
361 // 1/N_evt 1/(2 pi pt) 1/width 1/etarange
362 for (int ii=1; ii <= dNdPt_MC->GetNbinsX() ;ii++) {
363 Double_t pt = dNdPt_MC->GetBinCenter(ii);
364 Double_t width = dNdPt_MC->GetBinWidth(ii);
365 Double_t val = dNdPt_MC->GetBinContent(ii);
366 Double_t err = dNdPt_MC->GetBinError(ii);
367 Double_t cval = (val)/(width * 2.0 * TMath::Pi() * 1.6 * MCGeneratedEvents * pt);
368 Double_t cerr = (err)/(width * 2.0 * TMath::Pi() * 1.6 * MCGeneratedEvents * pt);
369 dNdPt_MC->SetBinContent(ii,cval);
370 dNdPt_MC->SetBinError(ii,cerr);
371 }
372 dNdPt_MC->SetMarkerStyle(21);
373 dNdPt_MC->SetTitle("; p_{T} (GeV/c) ; 1/N_{evt} 1/(2#pi p_{T}) (d^{2}N_{ch})/(d#eta dp_{T})^{-2}");
374 ContHist->Add(dNdPt_MC);
375
376
377
378 h1GenPrimTrack_pt = (TH1D*) h1GenPrimTrack_pt->Rebin(ptNbins,"",binsPt); // rebin hist
379
380 ContHist->Add(h3GenPrimTrack);
381 ContHist->Add(h2GenPrimTrack_zv_pt);
382 ContHist->Add(h2GenPrimTrack_zv_eta);
383 ContHist->Add(h2GenPrimTrack_pt_eta);
384 ContHist->Add(h1GenPrimTrack_zv);
385 ContHist->Add(h1GenPrimTrack_pt);
386 ContHist->Add(h1GenPrimTrack_eta);
387 printf("\n");
388 printf("==============================================================\n");
389 printf("=== recontructed MC tracks %lf ===\n",MCReconstructedTracks);
390 printf("=== recontructed MC secondary tracks %lf ===\n",MCReconstructedSecTracks);
391 printf("=== recontructed MC primary tracks %lf ===\n",MCReconstructedPrimTracks);
392 printf("=== generated MC primary track %lf ===\n",MCGeneratedPrimTracks);
393 printf("==============================================================\n");
394 printf("\n");
395
396// THnSparse *fSparseTriggerTrackEvent = obj->GetTriggerTrackEventMatrix();//Tracks from triggered events
397// THnSparse *fSparseVtxTrackEvent = obj->GetRecTrackEventMatrix();//Tracks from events with rec. vtx
398// THnSparse *fSparseGenTrackEvent = obj->GetGenTrackEventMatrix();//generated TrackEvent matrix
399
400 // tracking efficiencies + corrections
401 TH3D* h3TrackEff = AlidNdPtHelper::GenerateCorrMatrix(h3RecPrimTrack,h3GenPrimTrack,"h3TrackEff");
402 TH3D* h3TrackCorr = AlidNdPtHelper::GenerateCorrMatrix(h3GenPrimTrack,h3RecPrimTrack,"h3TrackCorr");
403
404 TH2D* h2TrackEff_zv_pt = AlidNdPtHelper::GenerateCorrMatrix(h2RecPrimTrack_zv_pt,h2GenPrimTrack_zv_pt,"h2TrackEff_zv_pt");
405 TH2D* h2TrackCorr_zv_pt = AlidNdPtHelper::GenerateCorrMatrix(h2GenPrimTrack_zv_pt,h2RecPrimTrack_zv_pt,"h2TrackCorr_zv_pt");
406 TH2D* h2TrackEff_zv_eta = AlidNdPtHelper::GenerateCorrMatrix(h2RecPrimTrack_zv_eta,h2GenPrimTrack_zv_eta,"h2TrackEff_zv_eta");
407 TH2D* h2TrackCorr_zv_eta = AlidNdPtHelper::GenerateCorrMatrix(h2GenPrimTrack_zv_eta,h2RecPrimTrack_zv_eta,"h2TrackCorr_zv_eta");
408 TH2D* h2TrackEff_pt_eta = AlidNdPtHelper::GenerateCorrMatrix(h2RecPrimTrack_pt_eta,h2GenPrimTrack_pt_eta,"h2TrackEff_pt_eta");
409 TH2D* h2TrackCorr_pt_eta = AlidNdPtHelper::GenerateCorrMatrix(h2GenPrimTrack_pt_eta,h2RecPrimTrack_pt_eta,"h2TrackCorr_pt_eta");
410
411 TH1D* h1TrackEff_zv = AlidNdPtHelper::GenerateCorrMatrix(h1RecPrimTrack_zv,h1GenPrimTrack_zv,"h1TrackEff_zv");
412 TH1D* h1TrackCorr_zv = AlidNdPtHelper::GenerateCorrMatrix(h1GenPrimTrack_zv,h1RecPrimTrack_zv,"h1TrackCorr_zv");
413 TH1D* h1TrackEff_pt = AlidNdPtHelper::GenerateCorrMatrix(h1RecPrimTrack_pt,h1GenPrimTrack_pt,"h1TrackEff_pt");
414 TH1D* h1TrackCorr_pt = AlidNdPtHelper::GenerateCorrMatrix(h1GenPrimTrack_pt,h1RecPrimTrack_pt,"h1TrackCorr_pt");
415 TH1D* h1TrackEff_eta = AlidNdPtHelper::GenerateCorrMatrix(h1RecPrimTrack_eta,h1GenPrimTrack_eta,"h1TrackEff_eta");
416 TH1D* h1TrackCorr_eta = AlidNdPtHelper::GenerateCorrMatrix(h1GenPrimTrack_eta,h1RecPrimTrack_eta,"h1TrackCorr_eta");
417
418 CorrMatr->Add(h3TrackEff);
419 CorrMatr->Add(h3TrackCorr);
420 CorrMatr->Add(h2TrackEff_zv_pt);
421 CorrMatr->Add(h2TrackCorr_zv_pt);
422 CorrMatr->Add(h2TrackEff_zv_eta);
423 CorrMatr->Add(h2TrackCorr_zv_eta);
424 CorrMatr->Add(h2TrackEff_pt_eta);
425 CorrMatr->Add(h2TrackCorr_pt_eta);
426 CorrMatr->Add(h1TrackEff_zv);
427 CorrMatr->Add(h1TrackCorr_zv);
428 CorrMatr->Add(h1TrackEff_pt);
429 CorrMatr->Add(h1TrackCorr_pt);
430 CorrMatr->Add(h1TrackEff_eta);
431 CorrMatr->Add(h1TrackCorr_eta);
432
433 // scale the secondaries before calculating correction matrices
434 for (Long64_t i = 0; i < fRecSecTrackMatrixScaled->GetNbins(); i++) {
435 Int_t c[3];
436 Double_t val = fRecSecTrackMatrixScaled->GetBinContent(i,c);
437 Double_t err = fRecSecTrackMatrixScaled->GetBinError(i);
438 Double_t pt = fRecSecTrackMatrixScaled->GetAxis(1)->GetBinCenter(c[1]);
439 Double_t scale = GetStrangenessCorrFactorPbPb(pt, sscale);
440 fRecSecTrackMatrixScaled->SetBinContent(c,val*scale);
441 fRecSecTrackMatrixScaled->SetBinError(c,err*scale);
442 }
443
444 // for correct determination of secondaries contamination, also the total total tracks have to be scaled
445 // this is done by taking primaries and adding the scaled secondaries
446 fRecTrackMatrixScaled->Add(fRecSecTrackMatrixScaled);
447
448 fRecSecTrackMatrixScaled->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
449 fRecSecTrackMatrixScaled->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
450 fRecSecTrackMatrixScaled->GetAxis(3)->SetRange(c_first,c_last); // select centrality
451
452 TH3D* h3RecSecTrackScaled = fRecSecTrackMatrixScaled->Projection(0,1,2)->Clone("h3RecSecTrackScaled");
453 TH2D* h2RecSecTrackScaled_zv_pt = fRecSecTrackMatrixScaled->Projection(0,1)->Clone("h2RecSecTrackScaled_zv_pt");
454 TH2D* h2RecSecTrackScaled_zv_eta = fRecSecTrackMatrixScaled->Projection(0,2)->Clone("h2RecSecTrackScaled_zv_eta");
455 TH2D* h2RecSecTrackScaled_pt_eta = fRecSecTrackMatrixScaled->Projection(1,2)->Clone("h2RecSecTrackScaled_pt_eta");
456 TH1D* h1RecSecTrackScaled_zv = fRecSecTrackMatrixScaled->Projection(0)->Clone("h1RecSecTrackScaled_zv");
457 TH1D* h1RecSecTrackScaled_pt = fRecSecTrackMatrixScaled->Projection(1)->Clone("h1RecSecTrackScaled_pt");
458 TH1D* h1RecSecTrackScaled_eta = fRecSecTrackMatrixScaled->Projection(2)->Clone("h1RecSecTrackScaled_eta");
459
460 h1RecSecTrackScaled_pt = (TH1D*) h1RecSecTrackScaled_pt->Rebin(ptNbins,"",binsPt); // rebin hist
461
462 ContHist->Add(h3RecSecTrackScaled);
463 ContHist->Add(h2RecSecTrackScaled_zv_pt);
464 ContHist->Add(h2RecSecTrackScaled_zv_eta);
465 ContHist->Add(h2RecSecTrackScaled_pt_eta);
466 ContHist->Add(h1RecSecTrackScaled_zv);
467 ContHist->Add(h1RecSecTrackScaled_pt);
468 ContHist->Add(h1RecSecTrackScaled_eta);
469
470 fRecTrackMatrixScaled->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
471 fRecTrackMatrixScaled->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
472
473 TH3D* h3RecTrackScaled = fRecTrackMatrixScaled->Projection(0,1,2)->Clone("h3RecTrackScaled");
474 TH2D* h2RecTrackScaled_zv_pt = fRecTrackMatrixScaled->Projection(0,1)->Clone("h2RecTrackScaled_zv_pt");
475 TH2D* h2RecTrackScaled_zv_eta = fRecTrackMatrixScaled->Projection(0,2)->Clone("h2RecTrackScaled_zv_eta");
476 TH2D* h2RecTrackScaled_pt_eta = fRecTrackMatrixScaled->Projection(1,2)->Clone("h2RecTrackScaled_pt_eta");
477 TH1D* h1RecTrackScaled_zv = fRecTrackMatrixScaled->Projection(0)->Clone("h1RecTrackScaled_zv");
478 TH1D* h1RecTrackScaled_pt = fRecTrackMatrixScaled->Projection(1)->Clone("h1RecTrackScaled_pt");
479 TH1D* h1RecTrackScaled_eta = fRecTrackMatrixScaled->Projection(2)->Clone("h1RecTrackScaled_eta");
480
481
482 h1RecTrackScaled_pt = (TH1D*) h1RecTrackScaled_pt->Rebin(ptNbins,"",binsPt); // rebin hist
483
484
485 ContHist->Add(h3RecTrackScaled);
486 ContHist->Add(h2RecTrackScaled_zv_pt);
487 ContHist->Add(h2RecTrackScaled_zv_eta);
488 ContHist->Add(h2RecTrackScaled_pt_eta);
489 ContHist->Add(h1RecTrackScaled_zv);
490 ContHist->Add(h1RecTrackScaled_pt);
491 ContHist->Add(h1RecTrackScaled_eta);
492
493 // create histograms for secondaries contamination and correction
494
495 TH3D* h3SecCont = AlidNdPtHelper::GenerateCorrMatrix(h3RecSecTrackScaled,h3RecTrackScaled,"h3SecCont");
496 TH3D* h3SecCorr = AlidNdPtHelper::GenerateContCorrMatrix(h3RecSecTrackScaled,h3RecTrackScaled,"h3SecCorr");
497 TH2D* h2SecCont_zv_pt = AlidNdPtHelper::GenerateCorrMatrix(h2RecSecTrackScaled_zv_pt,h2RecTrackScaled_zv_pt,"h2SecCont_zv_pt");
498 TH2D* h2SecCorr_zv_pt = AlidNdPtHelper::GenerateContCorrMatrix(h2RecSecTrackScaled_zv_pt,h2RecTrackScaled_zv_pt,"h2SecCorr_zv_pt");
499 TH2D* h2SecCont_zv_eta = AlidNdPtHelper::GenerateCorrMatrix(h2RecSecTrackScaled_zv_eta,h2RecTrackScaled_zv_eta,"h2SecCont_zv_eta");
500 TH2D* h2SecCorr_zv_eta = AlidNdPtHelper::GenerateContCorrMatrix(h2RecSecTrackScaled_zv_eta,h2RecTrackScaled_zv_eta,"h2SecCorr_zv_eta");
501 TH2D* h2SecCont_pt_eta = AlidNdPtHelper::GenerateCorrMatrix(h2RecSecTrackScaled_pt_eta,h2RecTrackScaled_pt_eta,"h2SecCont_pt_eta");
502 TH2D* h2SecCorr_pt_eta = AlidNdPtHelper::GenerateContCorrMatrix(h2RecSecTrackScaled_pt_eta,h2RecTrackScaled_pt_eta,"h2SecCorr_pt_eta");
503 TH1D* h1SecCont_zv = AlidNdPtHelper::GenerateCorrMatrix(h1RecSecTrackScaled_zv,h1RecTrackScaled_zv,"h1SecCont_zv");
504 TH1D* h1SecCorr_zv = AlidNdPtHelper::GenerateContCorrMatrix(h1RecSecTrackScaled_zv,h1RecTrackScaled_zv,"h1SecCorr_zv");
505 TH1D* h1SecCont_pt = AlidNdPtHelper::GenerateCorrMatrix(h1RecSecTrackScaled_pt,h1RecTrackScaled_pt,"h1SecCont_pt");
506 TH1D* h1SecCorr_pt = AlidNdPtHelper::GenerateContCorrMatrix(h1RecSecTrackScaled_pt,h1RecTrackScaled_pt,"h1SecCorr_pt");
507 TH1D* h1SecCont_eta = AlidNdPtHelper::GenerateCorrMatrix(h1RecSecTrackScaled_eta,h1RecTrackScaled_eta,"h1SecCont_eta");
508 TH1D* h1SecCorr_eta = AlidNdPtHelper::GenerateContCorrMatrix(h1RecSecTrackScaled_eta,h1RecTrackScaled_eta,"h1SecCorr_eta");
509
510 CorrMatr->Add(h3SecCont);
511 CorrMatr->Add(h3SecCorr);
512 CorrMatr->Add(h2SecCont_zv_pt);
513 CorrMatr->Add(h2SecCorr_zv_pt);
514 CorrMatr->Add(h2SecCont_zv_eta);
515 CorrMatr->Add(h2SecCorr_zv_eta);
516 CorrMatr->Add(h2SecCont_pt_eta);
517 CorrMatr->Add(h2SecCorr_pt_eta);
518 CorrMatr->Add(h1SecCont_zv);
519 CorrMatr->Add(h1SecCorr_zv);
520 CorrMatr->Add(h1SecCont_pt);
521 CorrMatr->Add(h1SecCorr_pt);
522 CorrMatr->Add(h1SecCont_eta);
523 CorrMatr->Add(h1SecCorr_eta);
524
525 // plot pictures and save to gifdir
526 for (i=0; i < CorrMatr->LastIndex(); i++) {
527 TCanvas* ctmp = PlotHist(CorrMatr->At(i),idstring);
528 if (gifdir && ctmp) {
529 TString gif(gifdir);
530 gif += '/';
531 gif += ctmp->GetName();
532 gif += ".gif";
533 ctmp->SaveAs(gif.Data(),"gif");
534 delete ctmp;
535 }
536 }
537 for (i=0; i < ContHist->LastIndex(); i++) {
538 TCanvas* ctmp = PlotHist(ContHist->At(i),idstring);
539 if (gifdir && ctmp) {
540 TString gif(gifdir);
541 gif += '/';
542 gif += ctmp->GetName();
543 gif += ".gif";
544 ctmp->SaveAs(gif.Data(),"gif");
545 delete ctmp;
546 }
547 }
548
549 // save all correction matrices and control histograms to file
550 if (!outfile) { return; }
551 TFile *out = TFile::Open(outfile,"RECREATE");
552 CorrMatr->Write();
553 ContHist->Write();
554 out->Close();
555
556 return MCReconstructedEvents;
557
558}
559
560
561//_____________________________________________________________________________
562Double_t GetStrangenessCorrFactorPbPb(Double_t pt, Double_t s)
563{
564 // data driven correction factor for secondaries (PbPb)
565
566 if (pt <= 0.25) return 1.0;
567 if (pt <= 0.5) return GetLinearInterpolationValue(0.25,1.0,0.5,0.60+0.40*s, pt);
568 if (pt <= 1.0) return GetLinearInterpolationValue(0.5,0.60+0.40*s,1.0,0.53+0.47*s, pt);
569 if (pt <= 2.0) return GetLinearInterpolationValue(1.0,0.53+0.47*s,2.0,0.44+0.56*s, pt);
570 if (pt <= 5.0) return GetLinearInterpolationValue(2.0,0.44+0.56*s,5.0,0.33+0.67*s, pt);
571 return 0.33+0.67*s;
572}
573
574
575//___________________________________________________________________________
576Double_t GetLinearInterpolationValue(const Double_t x1,const Double_t y1,const Double_t x2,const Double_t y2, const Double_t pt)
577{
578 //
579 // linear interpolation
580 //
581 return ((y2-y1)/(x2-x1))*pt+(y2-(((y2-y1)/(x2-x1))*x2));
582}
583
584//___________________________________________________________________________
585TCanvas* PlotHist(TObject* hobj, const char* label=0)
586{
587 TH1* h = dynamic_cast<TH1*>(hobj);
588 if (!h) return 0;
589 if (h->GetDimension() > 2) return 0;
590 h->SetStats(0);
591 if ( TString(h->GetName()).Contains("Events")) { h->SetStats(1); }
592 TString t(label);
593 if (label) t += "_";
594 t += h->GetName();
595 h->SetTitle(t.Data());
596 TCanvas* c = new TCanvas(t.Data(),t.Data());
597 if (h->GetDimension() >= 1) {
598 TString xlabel(h->GetXaxis()->GetTitle());
599 if (xlabel.Contains("Pt")) { c->SetLogx(); h->GetXaxis()->SetRangeUser(0.1 , 100.); }
600 }
601 if (h->GetDimension() == 1) {
602 if (xlabel.Contains("p_{T}")) { c->SetLogx(); c->SetLogy(); h->GetXaxis()->SetRangeUser(0.1 , 100.); }
603 }
604 if (h->GetDimension() == 2) {
605 TString ylabel(h->GetYaxis()->GetTitle());
606 if (ylabel.Contains("Pt")) { c->SetLogy(); h->GetYaxis()->SetRangeUser(0.1 , 100.); }
607 if (ylabel.Contains("p_{T}")) { c->SetLogy(); h->GetYaxis()->SetRangeUser(0.1 , 100.); }
608 h->Draw("COLZ");
609 }
610 if (h->GetDimension() == 1) {
611 h->Draw();
612 }
613 return c;
614
615}
616
617Int_t CheckLoadLibrary(const char* library)
618{
619 // checks if a library is already loaded, if not loads the library
620
621 if (strlen(gSystem->GetLibraries(Form("%s.so", library), "", kFALSE)) > 0)
622 return 1;
623
624 return gSystem->Load(library);
625}