PWGUD/dNdPt -> PWGLF/SPECTRA/ChargedHadrons/dNdPt
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / ChargedHadrons / dNdPt / macros / GenerateCorrMatr.C
CommitLineData
7c96ea29 1Double_t GenerateCorrMatr(const char* mcfile, const char* mctask, const char* idstring ,const char* outfile, const char* gifdir = 0)
2//Double_t GenerateCorrMatr()
3{
4
5//tmp setting
6// const char* mcfile = "/lustre/alice/train/V005.PbPb_MC/2010-11-29_0108.4251/mergedPeriods/MC/PbPb/LHC10h8/mknichel_dNdPtPbPb_gt_v0_c0_syst4.root";
7// const char* mctask = "mknichel_dNdPtPbPb_gt_v0_c0_syst4";
8// const char* idstring = "gt_v0_c0_syst4";
9// const char* outfile = "/u/mknichel/alice/dNdPt/2010-11-29_0318/corrMatr_gt_v0_c0_syst4.root";
10// const char* gifdir = "/u/mknichel/alice/dNdPt/2010-11-29";
11
12 // settings vor zVertex cut (event and track level)
13 Double_t zVert = 10.0;
14
15 // setting on eta cut (track level)
16 Double_t eta = 0.8;
17
18 // strangeness scaling factor (for secondaries from strange decays)
19 Double_t sscale = 1.4;
20
21 //load required libraries
22 //load required libraries
2bfe5463 23 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");
7c96ea29 24
25 gSystem->Load("libCore");
26 gSystem->Load("libPhysics");
27 gSystem->Load("libMinuit");
28 gSystem->Load("libGui");
29 gSystem->Load("libXMLParser");
30
31 gSystem->Load("libGeom");
32 gSystem->Load("libVMC");
33
34 gSystem->Load("libNet");
35
36
37 gSystem->Load("libSTEERBase");
38 gSystem->Load("libESD");
39 gSystem->Load("libCDB");
40 gSystem->Load("libRAWDatabase");
41 gSystem->Load("libRAWDatarec");
42 gSystem->Load("libANALYSIS");
43
44
45
46 gSystem->Load("libANALYSIS.so");
47 gSystem->Load("libANALYSISalice.so");
48 gSystem->Load("libTENDER.so");
49 gSystem->Load("libCORRFW.so");
50 gSystem->Load("libPWG0base.so");
51 gSystem->Load("libPWG0dep");
52 gSystem->Load("libPWG0selectors.so");
53
54
55
56
57
58 // make plots nicer
59 gROOT->SetStyle("Plain");
60 gStyle->SetPalette(1);
61
62 // array for all correction matrices
63 TObjArray* CorrMatr = new TObjArray();
64
65 // array for all control histograms
66 TObjArray* ContHist = new TObjArray();
67
68
69 // load mc information
70 TFile* fmc = TFile::Open(mcfile,"READ");
71 if (!fmc) return -1;
72 TList* lmc = dynamic_cast<TList*>(fmc->Get(mctask));
73 if (!lmc) return -1;
74 AlidNdPtAnalysisPbPb *obj = dynamic_cast<AlidNdPtAnalysisPbPb*>(lmc->FindObject("dNdPtAnalysisPbPb"));
75 if (!obj) return -1;
76
77 //Event statistics
78 THnSparse *fRecEventMatrix = obj->GetRecEventMatrix(); //all reconstructed events
79 TH2D* h2RecEventAll = fRecEventMatrix->Projection(0,1)->Clone("h2RecEventAll");
80 fRecEventMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
81 TH2D* h2RecEvent = fRecEventMatrix->Projection(0,1)->Clone("h2RecEvent");
82 Double_t MCReconstructedEvents = h2RecEvent->Integral();
83 Double_t MCReconstructedEventsAll = h2RecEventAll->Integral();
84 ContHist->Add(h2RecEvent);
85 ContHist->Add(h2RecEventAll);
86
87 THnSparse *fTriggerEventMatrix = obj->GetTriggerEventMatrix(); //all triggered events
88 TH2D* h2TriggerEventAll = fTriggerEventMatrix->Projection(0,1)->Clone("h2TriggerEventAll");
89 fTriggerEventMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
90 TH2D* h2TriggerEvent = fTriggerEventMatrix->Projection(0,1)->Clone("h2TriggerEvent");
91 Double_t MCTriggeredEvents = h2TriggerEvent->Integral();
92 Double_t MCTriggeredEventsAll = h2TriggerEventAll->Integral();
93 ContHist->Add(h2TriggerEvent);
94 ContHist->Add(h2TriggerEventAll);
95
96 THnSparse *fGenEventMatrix = obj->GetGenEventMatrix(); //all generated events
97 TH2D* h2GenEventAll = fGenEventMatrix->Projection(0,1)->Clone("h2GenEventAll");
98 fGenEventMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
99 TH2D* h2GenEvent = fGenEventMatrix->Projection(0,1)->Clone("h2GenEvent");
100 Double_t MCGeneratedEvents = h2GenEvent->Integral();
101 Double_t MCGeneratedEventsAll = h2GenEventAll->Integral();
102 ContHist->Add(h2RecEvent);
103 ContHist->Add(h2RecEventAll);
104
105 printf("=== generated MC events for correction matrices %lf ===\n",MCGeneratedEvents);
106 printf("=== triggered MC events for correction matrices %lf ===\n",MCTriggeredEvents);
107 printf("=== recontructed MC events for correction matrices %lf ===\n",MCReconstructedEvents);
108 printf("\n");
109 printf("=== cut on the zVertex +- %lf ===\n",zVert);
110 printf("=== generated MC events (before zVertex cut) %lf ===\n",MCGeneratedEventsAll);
111 printf("=== triggered MC events (before zVertex cut) %lf ===\n",MCTriggeredEventsAll);
112 printf("=== recontructed MC events (before zVertex cut) %lf ===\n",MCReconstructedEventsAll);
113
114 TH1D* h1MCGeneratedEvents = new TH1D("h1MCGeneratedEvents","h1MCGeneratedEvents",1,0,1);
115 TH1D* h1MCTriggeredEvents = new TH1D("h1MCTriggeredEvents","h1MCTriggeredEvents",1,0,1);
116 TH1D* h1MCReconstructedEvents = new TH1D("h1MCReconstructedEvents","h1MCReconstructedEvents",1,0,1);
117 TH1D* h1MCGeneratedEventsAll = new TH1D("h1MCGeneratedEventsAll","h1MCGeneratedEventsAll",1,0,1);
118 TH1D* h1MCTriggeredEventsAll = new TH1D("h1MCTriggeredEventsAll","h1MCTriggeredEventsAll",1,0,1);
119 TH1D* h1MCReconstructedEventsAll = new TH1D("h1MCReconstructedEventsAll","h1MCReconstructedEventsAll",1,0,1);
120
121 h1MCGeneratedEvents->Fill(0,MCGeneratedEvents);
122 h1MCGeneratedEvents->SetEntries(MCGeneratedEvents);
123 h1MCTriggeredEvents->Fill(0,MCTriggeredEvents);
124 h1MCTriggeredEvents->SetEntries(MCTriggeredEvents);
125 h1MCReconstructedEvents->Fill(0,MCReconstructedEvents);
126 h1MCReconstructedEvents->SetEntries(MCReconstructedEvents);
127 h1MCGeneratedEventsAll->Fill(0,MCGeneratedEventsAll);
128 h1MCGeneratedEventsAll->SetEntries(MCGeneratedEventsAll);
129 h1MCTriggeredEventsAll->Fill(0,MCTriggeredEventsAll);
130 h1MCTriggeredEventsAll->SetEntries(MCTriggeredEventsAll);
131 h1MCReconstructedEventsAll->Fill(0,MCReconstructedEventsAll);
132 h1MCReconstructedEventsAll->SetEntries(MCReconstructedEventsAll);
133
134 ContHist->Add(h1MCGeneratedEvents);
135 ContHist->Add(h1MCTriggeredEvents);
136 ContHist->Add(h1MCReconstructedEvents);
137 ContHist->Add(h1MCGeneratedEventsAll);
138 ContHist->Add(h1MCTriggeredEventsAll);
139 ContHist->Add(h1MCReconstructedEventsAll);
140
141 // efficienfy and correction matrices for tigger and vertex efficiency
142 TH2D* h2EventTriggerEffAll = AlidNdPtHelper::GenerateCorrMatrix(h2TriggerEventAll,h2GenEventAll,"h2EventTriggerEffAll");
143 TH2D* h2EventTriggerCorrAll = AlidNdPtHelper::GenerateCorrMatrix(h2GenEventAll,h2TriggerEventAll,"h2EventTriggerCorrAll");
144 TH2D* h2EventTriggerEff = AlidNdPtHelper::GenerateCorrMatrix(h2TriggerEvent,h2GenEvent,"h2EventTriggerEff");
145 TH2D* h2EventTriggerCorr = AlidNdPtHelper::GenerateCorrMatrix(h2GenEvent,h2TriggerEvent,"h2EventTriggerCorr");
146
147 TH2D* h2EventRecEffAll = AlidNdPtHelper::GenerateCorrMatrix(h2RecEventAll,h2TriggerEventAll,"h2EventRecEffAll");
148 TH2D* h2EventRecCorrAll = AlidNdPtHelper::GenerateCorrMatrix(h2TriggerEventAll,h2RecEventAll,"h2EventRecCorrAll");
149 TH2D* h2EventRecEff = AlidNdPtHelper::GenerateCorrMatrix(h2RecEvent,h2TriggerEvent,"h2EventRecEff");
150 TH2D* h2EventRecCorr = AlidNdPtHelper::GenerateCorrMatrix(h2TriggerEvent,h2RecEvent,"h2EventRecCorr");
151
152 TH2D* h2EventEffAll = AlidNdPtHelper::GenerateCorrMatrix(h2RecEventAll,h2GenEventAll,"h2EventEffAll");
153 TH2D* h2EventCorrAll = AlidNdPtHelper::GenerateCorrMatrix(h2GenEventAll,h2RecEventAll,"h2EventCorrAll");
154 TH2D* h2EventEff = AlidNdPtHelper::GenerateCorrMatrix(h2RecEvent,h2GenEvent,"h2EventEff");
155 TH2D* h2EventCorr = AlidNdPtHelper::GenerateCorrMatrix(h2GenEvent,h2RecEvent,"h2EventCorr");
156
157 CorrMatr->Add(h2EventTriggerEffAll);
158 CorrMatr->Add(h2EventTriggerCorrAll);
159 CorrMatr->Add(h2EventTriggerEff);
160 CorrMatr->Add(h2EventTriggerCorr);
161 CorrMatr->Add(h2EventRecEffAll);
162 CorrMatr->Add(h2EventRecCorrAll);
163 CorrMatr->Add(h2EventRecEff);
164 CorrMatr->Add(h2EventRecCorr);
165 CorrMatr->Add(h2EventEffAll);
166 CorrMatr->Add(h2EventCorrAll);
167 CorrMatr->Add(h2EventEff);
168 CorrMatr->Add(h2EventCorr);
169
170
171
172 // all recontructed
173 THnSparse *fRecTrackMatrix = obj->GetRecTrackMatrix();
174 fRecTrackMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
175 fRecTrackMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
176 TH3D* h3RecTrack = fRecTrackMatrix->Projection(0,1,2)->Clone("h3RecTrack");
177 TH2D* h2RecTrack_zv_pt = fRecTrackMatrix->Projection(0,1)->Clone("h2RecTrack_zv_pt");
178 TH2D* h2RecTrack_zv_eta = fRecTrackMatrix->Projection(0,2)->Clone("h2RecTrack_zv_eta");
179 TH2D* h2RecTrack_pt_eta = fRecTrackMatrix->Projection(1,2)->Clone("h2RecTrack_pt_eta");
180 TH1D* h1RecTrack_zv = fRecTrackMatrix->Projection(0)->Clone("h1RecTrack_zv");
181 TH1D* h1RecTrack_pt = fRecTrackMatrix->Projection(1)->Clone("h1RecTrack_pt");
182 TH1D* h1RecTrack_eta = fRecTrackMatrix->Projection(2)->Clone("h1RecTrack_eta");
183 Double_t MCReconstructedTracks = h3RecTrack->Integral();
184
185 ContHist->Add(h3RecTrack);
186 ContHist->Add(h2RecTrack_zv_pt);
187 ContHist->Add(h2RecTrack_zv_eta);
188 ContHist->Add(h2RecTrack_pt_eta);
189 ContHist->Add(h1RecTrack_zv);
190 ContHist->Add(h1RecTrack_pt);
191 ContHist->Add(h1RecTrack_eta);
192
193 // recontructed primary tracks
194 THnSparse *fRecPrimTrackMatrix = obj->GetRecPrimTrackMatrix();
195 THnSparse *fRecTrackMatrixScaled = fRecPrimTrackMatrix->Clone("fRecTrackMatrixScaled"); //used later for secondaries scaling
196 fRecPrimTrackMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
197 fRecPrimTrackMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
198 TH3D* h3RecPrimTrack = fRecPrimTrackMatrix->Projection(0,1,2)->Clone("h3RecPrimTrack");
199 TH2D* h2RecPrimTrack_zv_pt = fRecPrimTrackMatrix->Projection(0,1)->Clone("h2RecPrimTrack_zv_pt");
200 TH2D* h2RecPrimTrack_zv_eta = fRecPrimTrackMatrix->Projection(0,2)->Clone("h2RecPrimTrack_zv_eta");
201 TH2D* h2RecPrimTrack_pt_eta = fRecPrimTrackMatrix->Projection(1,2)->Clone("h2RecPrimTrack_pt_eta");
202 TH1D* h1RecPrimTrack_zv = fRecPrimTrackMatrix->Projection(0)->Clone("h1RecPrimTrack_zv");
203 TH1D* h1RecPrimTrack_pt = fRecPrimTrackMatrix->Projection(1)->Clone("h1RecPrimTrack_pt");
204 TH1D* h1RecPrimTrack_eta = fRecPrimTrackMatrix->Projection(2)->Clone("h1RecPrimTrack_eta");
205 Double_t MCReconstructedPrimTracks = h3RecPrimTrack->Integral();
206
207 ContHist->Add(h3RecPrimTrack);
208 ContHist->Add(h2RecPrimTrack_zv_pt);
209 ContHist->Add(h2RecPrimTrack_zv_eta);
210 ContHist->Add(h2RecPrimTrack_pt_eta);
211 ContHist->Add(h1RecPrimTrack_zv);
212 ContHist->Add(h1RecPrimTrack_pt);
213 ContHist->Add(h1RecPrimTrack_eta);
214
215 // recontructed secondary tracks
216 THnSparse *fRecSecTrackMatrix = obj->GetRecSecTrackMatrix();
217 THnSparse *fRecSecTrackMatrixScaled = fRecSecTrackMatrix->Clone("fRecSecTrackMatrixScaled"); //used later for secondaries scaling
218 fRecSecTrackMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
219 fRecSecTrackMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
220 TH3D* h3RecSecTrack = fRecSecTrackMatrix->Projection(0,1,2)->Clone("h3RecSecTrack");
221 TH2D* h2RecSecTrack_zv_pt = fRecSecTrackMatrix->Projection(0,1)->Clone("h2RecSecTrack_zv_pt");
222 TH2D* h2RecSecTrack_zv_eta = fRecSecTrackMatrix->Projection(0,2)->Clone("h2RecSecTrack_zv_eta");
223 TH2D* h2RecSecTrack_pt_eta = fRecSecTrackMatrix->Projection(1,2)->Clone("h2RecSecTrack_pt_eta");
224 TH1D* h1RecSecTrack_zv = fRecSecTrackMatrix->Projection(0)->Clone("h1RecSecTrack_zv");
225 TH1D* h1RecSecTrack_pt = fRecSecTrackMatrix->Projection(1)->Clone("h1RecSecTrack_pt");
226 TH1D* h1RecSecTrack_eta = fRecSecTrackMatrix->Projection(2)->Clone("h1RecSecTrack_eta");
227 Double_t MCReconstructedSecTracks = h3RecSecTrack->Integral();
228
229 ContHist->Add(h3RecSecTrack);
230 ContHist->Add(h2RecSecTrack_zv_pt);
231 ContHist->Add(h2RecSecTrack_zv_eta);
232 ContHist->Add(h2RecSecTrack_pt_eta);
233 ContHist->Add(h1RecSecTrack_zv);
234 ContHist->Add(h1RecSecTrack_pt);
235 ContHist->Add(h1RecSecTrack_eta);
236
237 // generated primary tracks
238 THnSparse *fGenPrimTrackMatrix = obj->GetGenPrimTrackMatrix();
239 fGenPrimTrackMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
240 fGenPrimTrackMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
241 TH3D* h3GenPrimTrack = fGenPrimTrackMatrix->Projection(0,1,2)->Clone("h3GenPrimTrack");
242 TH2D* h2GenPrimTrack_zv_Pt = fGenPrimTrackMatrix->Projection(0,1)->Clone("h2GenPrimTrack_zv_pt");
243 TH2D* h2GenPrimTrack_zv_eta = fGenPrimTrackMatrix->Projection(0,2)->Clone("h2GenPrimTrack_zv_eta");
244 TH2D* h2GenPrimTrack_pt_eta = fGenPrimTrackMatrix->Projection(1,2)->Clone("h2GenPrimTrack_pt_eta");
245 TH1D* h1GenPrimTrack_zv = fGenPrimTrackMatrix->Projection(0)->Clone("h1GenPrimTrack_zv");
246 TH1D* h1GenPrimTrack_pt = fGenPrimTrackMatrix->Projection(1)->Clone("h1GenPrimTrack_pt");
247 TH1D* h1GenPrimTrack_eta = fGenPrimTrackMatrix->Projection(2)->Clone("h1GenPrimTrack_eta");
248 Double_t MCGeneratedPrimTracks = h3GenPrimTrack->Integral();
249
250 ContHist->Add(h3GenPrimTrack);
251 ContHist->Add(h2GenPrimTrack_zv_pt);
252 ContHist->Add(h2GenPrimTrack_zv_eta);
253 ContHist->Add(h2GenPrimTrack_pt_eta);
254 ContHist->Add(h1GenPrimTrack_zv);
255 ContHist->Add(h1GenPrimTrack_pt);
256 ContHist->Add(h1GenPrimTrack_eta);
257 printf("\n");
258 printf("==============================================================\n");
259 printf("=== recontructed MC tracks %lf ===\n",MCReconstructedTracks);
260 printf("=== recontructed MC secondary tracks %lf ===\n",MCReconstructedSecTracks);
261 printf("=== recontructed MC primary tracks %lf ===\n",MCReconstructedPrimTracks);
262 printf("=== generated MC primary track %lf ===\n",MCGeneratedPrimTracks);
263 printf("==============================================================\n");
264 printf("\n");
265
266// THnSparse *fSparseTriggerTrackEvent = obj->GetTriggerTrackEventMatrix();//Tracks from triggered events
267// THnSparse *fSparseVtxTrackEvent = obj->GetRecTrackEventMatrix();//Tracks from events with rec. vtx
268// THnSparse *fSparseGenTrackEvent = obj->GetGenTrackEventMatrix();//generated TrackEvent matrix
269
270 // tracking efficiencies + corrections
271 TH3D* h3TrackEff = AlidNdPtHelper::GenerateCorrMatrix(h3RecPrimTrack,h3GenPrimTrack,"h3TrackEff");
272 TH3D* h3TrackCorr = AlidNdPtHelper::GenerateCorrMatrix(h3GenPrimTrack,h3RecPrimTrack,"h3TrackCorr");
273
274 TH2D* h2TrackEff_zv_pt = AlidNdPtHelper::GenerateCorrMatrix(h2RecPrimTrack_zv_pt,h2GenPrimTrack_zv_pt,"h2TrackEff_zv_pt");
275 TH2D* h2TrackCorr_zv_pt = AlidNdPtHelper::GenerateCorrMatrix(h2GenPrimTrack_zv_pt,h2RecPrimTrack_zv_pt,"h2TrackCorr_zv_pt");
276 TH2D* h2TrackEff_zv_eta = AlidNdPtHelper::GenerateCorrMatrix(h2RecPrimTrack_zv_eta,h2GenPrimTrack_zv_eta,"h2TrackEff_zv_eta");
277 TH2D* h2TrackCorr_zv_eta = AlidNdPtHelper::GenerateCorrMatrix(h2GenPrimTrack_zv_eta,h2RecPrimTrack_zv_eta,"h2TrackCorr_zv_eta");
278 TH2D* h2TrackEff_pt_eta = AlidNdPtHelper::GenerateCorrMatrix(h2RecPrimTrack_pt_eta,h2GenPrimTrack_pt_eta,"h2TrackEff_pt_eta");
279 TH2D* h2TrackCorr_pt_eta = AlidNdPtHelper::GenerateCorrMatrix(h2GenPrimTrack_pt_eta,h2RecPrimTrack_pt_eta,"h2TrackCorr_pt_eta");
280
281 TH1D* h1TrackEff_zv = AlidNdPtHelper::GenerateCorrMatrix(h1RecPrimTrack_zv,h1GenPrimTrack_zv,"h1TrackEff_zv");
282 TH1D* h1TrackCorr_zv = AlidNdPtHelper::GenerateCorrMatrix(h1GenPrimTrack_zv,h1RecPrimTrack_zv,"h1TrackCorr_zv");
283 TH1D* h1TrackEff_pt = AlidNdPtHelper::GenerateCorrMatrix(h1RecPrimTrack_pt,h1GenPrimTrack_pt,"h1TrackEff_pt");
284 TH1D* h1TrackCorr_pt = AlidNdPtHelper::GenerateCorrMatrix(h1GenPrimTrack_pt,h1RecPrimTrack_pt,"h1TrackCorr_pt");
285 TH1D* h1TrackEff_eta = AlidNdPtHelper::GenerateCorrMatrix(h1RecPrimTrack_eta,h1GenPrimTrack_eta,"h1TrackEff_eta");
286 TH1D* h1TrackCorr_eta = AlidNdPtHelper::GenerateCorrMatrix(h1GenPrimTrack_eta,h1RecPrimTrack_eta,"h1TrackCorr_eta");
287
288 CorrMatr->Add(h3TrackEff);
289 CorrMatr->Add(h3TrackCorr);
290 CorrMatr->Add(h2TrackEff_zv_pt);
291 CorrMatr->Add(h2TrackCorr_zv_pt);
292 CorrMatr->Add(h2TrackEff_zv_eta);
293 CorrMatr->Add(h2TrackCorr_zv_eta);
294 CorrMatr->Add(h2TrackEff_pt_eta);
295 CorrMatr->Add(h2TrackCorr_pt_eta);
296 CorrMatr->Add(h1TrackEff_zv);
297 CorrMatr->Add(h1TrackCorr_zv);
298 CorrMatr->Add(h1TrackEff_pt);
299 CorrMatr->Add(h1TrackCorr_pt);
300 CorrMatr->Add(h1TrackEff_eta);
301 CorrMatr->Add(h1TrackCorr_eta);
302
303 // scale the secondaries before calculating correction matrices
304 for (Long64_t i = 0; i < fRecSecTrackMatrixScaled->GetNbins(); i++) {
305 Int_t c[3];
306 Double_t val = fRecSecTrackMatrixScaled->GetBinContent(i,c);
307 Double_t err = fRecSecTrackMatrixScaled->GetBinError(i);
308 Double_t pt = fRecSecTrackMatrixScaled->GetAxis(1)->GetBinCenter(c[1]);
309 Double_t scale = GetStrangenessCorrFactorPbPb(pt, sscale);
310 fRecSecTrackMatrixScaled->SetBinContent(c,val*scale);
311 fRecSecTrackMatrixScaled->SetBinError(c,err*scale);
312 }
313
314 // for correct determination of secondaries contamination, also the total total tracks have to be scaled
315 // this is done by taking primaries and adding the scaled secondaries
316 fRecTrackMatrixScaled->Add(fRecSecTrackMatrixScaled);
317
318 fRecSecTrackMatrixScaled->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
319 fRecSecTrackMatrixScaled->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
320
321 TH3D* h3RecSecTrackScaled = fRecSecTrackMatrixScaled->Projection(0,1,2)->Clone("h3RecSecTrackScaled");
322 TH2D* h2RecSecTrackScaled_zv_pt = fRecSecTrackMatrixScaled->Projection(0,1)->Clone("h2RecSecTrackScaled_zv_pt");
323 TH2D* h2RecSecTrackScaled_zv_eta = fRecSecTrackMatrixScaled->Projection(0,2)->Clone("h2RecSecTrackScaled_zv_eta");
324 TH2D* h2RecSecTrackScaled_pt_eta = fRecSecTrackMatrixScaled->Projection(1,2)->Clone("h2RecSecTrackScaled_pt_eta");
325 TH1D* h1RecSecTrackScaled_zv = fRecSecTrackMatrixScaled->Projection(0)->Clone("h1RecSecTrackScaled_zv");
326 TH1D* h1RecSecTrackScaled_pt = fRecSecTrackMatrixScaled->Projection(1)->Clone("h1RecSecTrackScaled_pt");
327 TH1D* h1RecSecTrackScaled_eta = fRecSecTrackMatrixScaled->Projection(2)->Clone("h1RecSecTrackScaled_eta");
328
329 ContHist->Add(h3RecSecTrackScaled);
330 ContHist->Add(h2RecSecTrackScaled_zv_pt);
331 ContHist->Add(h2RecSecTrackScaled_zv_eta);
332 ContHist->Add(h2RecSecTrackScaled_pt_eta);
333 ContHist->Add(h1RecSecTrackScaled_zv);
334 ContHist->Add(h1RecSecTrackScaled_pt);
335 ContHist->Add(h1RecSecTrackScaled_eta);
336
337 fRecTrackMatrixScaled->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
338 fRecTrackMatrixScaled->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
339
340 TH3D* h3RecTrackScaled = fRecTrackMatrixScaled->Projection(0,1,2)->Clone("h3RecTrackScaled");
341 TH2D* h2RecTrackScaled_zv_pt = fRecTrackMatrixScaled->Projection(0,1)->Clone("h2RecTrackScaled_zv_pt");
342 TH2D* h2RecTrackScaled_zv_eta = fRecTrackMatrixScaled->Projection(0,2)->Clone("h2RecTrackScaled_zv_eta");
343 TH2D* h2RecTrackScaled_pt_eta = fRecTrackMatrixScaled->Projection(1,2)->Clone("h2RecTrackScaled_pt_eta");
344 TH1D* h1RecTrackScaled_zv = fRecTrackMatrixScaled->Projection(0)->Clone("h1RecTrackScaled_zv");
345 TH1D* h1RecTrackScaled_pt = fRecTrackMatrixScaled->Projection(1)->Clone("h1RecTrackScaled_pt");
346 TH1D* h1RecTrackScaled_eta = fRecTrackMatrixScaled->Projection(2)->Clone("h1RecTrackScaled_eta");
347
348 ContHist->Add(h3RecTrackScaled);
349 ContHist->Add(h2RecTrackScaled_zv_pt);
350 ContHist->Add(h2RecTrackScaled_zv_eta);
351 ContHist->Add(h2RecTrackScaled_pt_eta);
352 ContHist->Add(h1RecTrackScaled_zv);
353 ContHist->Add(h1RecTrackScaled_pt);
354 ContHist->Add(h1RecTrackScaled_eta);
355
356 // create histograms for secondaries contamination and correction
357
358 TH3D* h3SecCont = AlidNdPtHelper::GenerateCorrMatrix(h3RecSecTrackScaled,h3RecTrackScaled,"h3SecCont");
359 TH3D* h3SecCorr = AlidNdPtHelper::GenerateContCorrMatrix(h3RecSecTrackScaled,h3RecTrackScaled,"h3SecCorr");
360 TH2D* h2SecCont_zv_pt = AlidNdPtHelper::GenerateCorrMatrix(h2RecSecTrackScaled_zv_pt,h2RecTrackScaled_zv_pt,"h2SecCont_zv_pt");
361 TH2D* h2SecCorr_zv_pt = AlidNdPtHelper::GenerateContCorrMatrix(h2RecSecTrackScaled_zv_pt,h2RecTrackScaled_zv_pt,"h2SecCorr_zv_pt");
362 TH2D* h2SecCont_zv_eta = AlidNdPtHelper::GenerateCorrMatrix(h2RecSecTrackScaled_zv_eta,h2RecTrackScaled_zv_eta,"h2SecCont_zv_eta");
363 TH2D* h2SecCorr_zv_eta = AlidNdPtHelper::GenerateContCorrMatrix(h2RecSecTrackScaled_zv_eta,h2RecTrackScaled_zv_eta,"h2SecCorr_zv_eta");
364 TH2D* h2SecCont_pt_eta = AlidNdPtHelper::GenerateCorrMatrix(h2RecSecTrackScaled_pt_eta,h2RecTrackScaled_pt_eta,"h2SecCont_pt_eta");
365 TH2D* h2SecCorr_pt_eta = AlidNdPtHelper::GenerateContCorrMatrix(h2RecSecTrackScaled_pt_eta,h2RecTrackScaled_pt_eta,"h2SecCorr_pt_eta");
366 TH1D* h1SecCont_zv = AlidNdPtHelper::GenerateCorrMatrix(h1RecSecTrackScaled_zv,h1RecTrackScaled_zv,"h1SecCont_zv");
367 TH1D* h1SecCorr_zv = AlidNdPtHelper::GenerateContCorrMatrix(h1RecSecTrackScaled_zv,h1RecTrackScaled_zv,"h1SecCorr_zv");
368 TH1D* h1SecCont_pt = AlidNdPtHelper::GenerateCorrMatrix(h1RecSecTrackScaled_pt,h1RecTrackScaled_pt,"h1SecCont_pt");
369 TH1D* h1SecCorr_pt = AlidNdPtHelper::GenerateContCorrMatrix(h1RecSecTrackScaled_pt,h1RecTrackScaled_pt,"h1SecCorr_pt");
370 TH1D* h1SecCont_eta = AlidNdPtHelper::GenerateCorrMatrix(h1RecSecTrackScaled_eta,h1RecTrackScaled_eta,"h1SecCont_eta");
371 TH1D* h1SecCorr_eta = AlidNdPtHelper::GenerateContCorrMatrix(h1RecSecTrackScaled_eta,h1RecTrackScaled_eta,"h1SecCorr_eta");
372
373 CorrMatr->Add(h3SecCont);
374 CorrMatr->Add(h3SecCorr);
375 CorrMatr->Add(h2SecCont_zv_pt);
376 CorrMatr->Add(h2SecCorr_zv_pt);
377 CorrMatr->Add(h2SecCont_zv_eta);
378 CorrMatr->Add(h2SecCorr_zv_eta);
379 CorrMatr->Add(h2SecCont_pt_eta);
380 CorrMatr->Add(h2SecCorr_pt_eta);
381 CorrMatr->Add(h1SecCont_zv);
382 CorrMatr->Add(h1SecCorr_zv);
383 CorrMatr->Add(h1SecCont_pt);
384 CorrMatr->Add(h1SecCorr_pt);
385 CorrMatr->Add(h1SecCont_eta);
386 CorrMatr->Add(h1SecCorr_eta);
387
388 // plot pictures and save to gifdir
389 for (i=0; i < CorrMatr->LastIndex(); i++) {
390 TCanvas* ctmp = PlotHist(CorrMatr->At(i),idstring);
391 if (gifdir && ctmp) {
392 TString gif(gifdir);
393 gif += '/';
394 gif += ctmp->GetName();
395 gif += ".gif";
396 ctmp->SaveAs(gif.Data(),"gif");
397 delete ctmp;
398 }
399 }
400 for (i=0; i < ContHist->LastIndex(); i++) {
401 TCanvas* ctmp = PlotHist(ContHist->At(i),idstring);
402 if (gifdir && ctmp) {
403 TString gif(gifdir);
404 gif += '/';
405 gif += ctmp->GetName();
406 gif += ".gif";
407 ctmp->SaveAs(gif.Data(),"gif");
408 delete ctmp;
409 }
410 }
411
412 // save all correction matrices and control histograms to file
413 if (!outfile) { return; }
414 TFile *out = TFile::Open(outfile,"RECREATE");
415 CorrMatr->Write();
416 ContHist->Write();
417 out->Close();
418
419 return MCReconstructedEvents;
420
421}
422
423
424//_____________________________________________________________________________
425Double_t GetStrangenessCorrFactorPbPb(Double_t pt, Double_t s)
426{
427 // data driven correction factor for secondaries (PbPb)
428
429 if (pt <= 0.25) return 1.0;
430 if (pt <= 0.5) return GetLinearInterpolationValue(0.25,1.0,0.5,0.60+0.40*s, pt);
431 if (pt <= 1.0) return GetLinearInterpolationValue(0.5,0.60+0.40*s,1.0,0.53+0.47*s, pt);
432 if (pt <= 2.0) return GetLinearInterpolationValue(1.0,0.53+0.47*s,2.0,0.44+0.56*s, pt);
433 if (pt <= 5.0) return GetLinearInterpolationValue(2.0,0.44+0.56*s,5.0,0.33+0.67*s, pt);
434 return 0.33+0.67*s;
435}
436
437
438//___________________________________________________________________________
439Double_t GetLinearInterpolationValue(const Double_t x1,const Double_t y1,const Double_t x2,const Double_t y2, const Double_t pt)
440{
441 //
442 // linear interpolation
443 //
444 return ((y2-y1)/(x2-x1))*pt+(y2-(((y2-y1)/(x2-x1))*x2));
445}
446
447//___________________________________________________________________________
448TCanvas* PlotHist(TObject* hobj, const char* label=0)
449{
450 TH1* h = dynamic_cast<TH1*>(hobj);
451 if (!h) return 0;
452 if (h->GetDimension() > 2) return 0;
453 h->SetStats(0);
454 if ( TString(h->GetName()).Contains("Events")) { h->SetStats(1); }
455 TString t(label);
456 if (label) t += "_";
457 t += h->GetName();
458 h->SetTitle(t.Data());
459 TCanvas* c = new TCanvas(t.Data(),t.Data());
460 if (h->GetDimension() >= 1) {
461 TString xlabel(h->GetXaxis()->GetTitle());
462 if (xlabel.Contains("Pt")) { c->SetLogx(); h->GetXaxis()->SetRangeUser(0.1 , 50.); }
463 }
464 if (h->GetDimension() == 1) {
465 if (xlabel.Contains("p_{T}")) { c->SetLogx(); c->SetLogy(); h->GetXaxis()->SetRangeUser(0.1 , 50.); }
466 }
467 if (h->GetDimension() == 2) {
468 TString ylabel(h->GetYaxis()->GetTitle());
469 if (ylabel.Contains("Pt")) { c->SetLogy(); h->GetYaxis()->SetRangeUser(0.1 , 50.); }
470 if (ylabel.Contains("p_{T}")) { c->SetLogy(); h->GetYaxis()->SetRangeUser(0.1 , 50.); }
471 h->Draw("COLZ");
472 }
473 if (h->GetDimension() == 1) {
474 h->Draw();
475 }
476 return c;
477
478}
479
480Int_t CheckLoadLibrary(const char* library)
481{
482 // checks if a library is already loaded, if not loads the library
483
484 if (strlen(gSystem->GetLibraries(Form("%s.so", library), "", kFALSE)) > 0)
485 return 1;
486
487 return gSystem->Load(library);
488}