]>
Commit | Line | Data |
---|---|---|
10ebe68d | 1 | /* $Id$ */ |
2 | ||
3 | #include "AlidNdEtaSystematicsSelector.h" | |
4 | ||
5 | #include <TStyle.h> | |
6 | #include <TSystem.h> | |
7 | #include <TCanvas.h> | |
8 | #include <TVector3.h> | |
9 | #include <TChain.h> | |
10 | #include <TFile.h> | |
11 | #include <TH3F.h> | |
12 | #include <TH1F.h> | |
13 | #include <TParticle.h> | |
14 | ||
15 | #include <AliLog.h> | |
16 | #include <AliESD.h> | |
17 | #include <AliGenEventHeader.h> | |
18 | #include <AliStack.h> | |
19 | #include <AliHeader.h> | |
20 | ||
21 | #include "esdTrackCuts/AliESDtrackCuts.h" | |
22 | #include "AliPWG0Helper.h" | |
23 | #include "AlidNdEtaCorrection.h" | |
24 | ||
25 | ClassImp(AlidNdEtaSystematicsSelector) | |
26 | ||
27 | AlidNdEtaSystematicsSelector::AlidNdEtaSystematicsSelector() : | |
28 | AliSelectorRL(), | |
29 | fSecondaries(0), | |
9f469bf5 | 30 | fSigmaVertex(0), |
f31d5d49 | 31 | fEsdTrackCuts(0), |
32 | fOverallPrimaries(0), | |
61385583 | 33 | fOverallSecondaries(0), |
34 | fPIDParticles(0), | |
35 | fPIDTracks(0) | |
10ebe68d | 36 | { |
37 | // | |
38 | // Constructor. Initialization of pointers | |
39 | // | |
40 | ||
41 | for (Int_t i=0; i<4; ++i) | |
42 | fdNdEtaCorrection[i] = 0; | |
43 | } | |
44 | ||
45 | AlidNdEtaSystematicsSelector::~AlidNdEtaSystematicsSelector() | |
46 | { | |
47 | // | |
48 | // Destructor | |
49 | // | |
50 | } | |
51 | ||
52 | void AlidNdEtaSystematicsSelector::Begin(TTree* tree) | |
53 | { | |
54 | // Begin function | |
55 | ||
56 | AliSelectorRL::Begin(tree); | |
57 | ||
58 | ReadUserObjects(tree); | |
59 | } | |
60 | ||
61 | void AlidNdEtaSystematicsSelector::ReadUserObjects(TTree* tree) | |
62 | { | |
63 | // read the user objects, called from slavebegin and begin | |
64 | ||
65 | if (!fEsdTrackCuts && fInput) | |
66 | fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fInput->FindObject("AliESDtrackCuts")); | |
67 | ||
68 | if (!fEsdTrackCuts && tree) | |
69 | fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (tree->GetUserInfo()->FindObject("AliESDtrackCuts")); | |
70 | ||
71 | if (!fEsdTrackCuts) | |
72 | AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from input list."); | |
73 | } | |
74 | ||
75 | void AlidNdEtaSystematicsSelector::SlaveBegin(TTree* tree) | |
76 | { | |
77 | // The SlaveBegin() function is called after the Begin() function. | |
78 | // When running with PROOF SlaveBegin() is called on each slave server. | |
79 | // The tree argument is deprecated (on PROOF 0 is passed). | |
80 | ||
81 | AliSelector::SlaveBegin(tree); | |
82 | ||
83 | ReadUserObjects(tree); | |
84 | ||
85 | TString option(GetOption()); | |
86 | ||
87 | printf("Running AlidNdEtaSystematicsSelector with options %s\n", option.Data()); | |
88 | ||
89 | if (option.Contains("secondaries")) | |
90 | { | |
9f469bf5 | 91 | fSecondaries = new TH3F("fSecondaries", "fSecondaries;NacceptedTracks;CratioSecondaries;p_{T}", 2000, -0.5, 205.5, 16, 0.45, 2.05, 10, 0, 10); |
10ebe68d | 92 | } |
93 | ||
94 | if (option.Contains("particle-composition")) | |
95 | { | |
96 | for (Int_t i=0; i<4; ++i) | |
97 | { | |
98 | TString name; | |
99 | name.Form("correction_%d", i); | |
100 | fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name); | |
101 | } | |
102 | } | |
9f469bf5 | 103 | |
104 | if (option.Contains("sigma-vertex")) | |
105 | { | |
106 | fSigmaVertex = new TH1F("fSigmaVertex", "fSigmaVertex;Nsigma2vertex;NacceptedTracks", 10, 0.25, 5.25); | |
107 | printf("WARNING: sigma-vertex analysis enabled. This will produce weird results in the AliESDtrackCuts histograms\n"); | |
108 | } | |
61385583 | 109 | |
110 | fPIDParticles = new TH1F("pid_particles", "PID of generated primary particles", 10001, -5000.5, 5000.5); | |
111 | fPIDTracks = new TH1F("pid_tracks", "MC PID of reconstructed tracks", 10001, -5000.5, 5000.5); | |
10ebe68d | 112 | } |
113 | ||
114 | Bool_t AlidNdEtaSystematicsSelector::Process(Long64_t entry) | |
115 | { | |
116 | // The Process() function is called for each entry in the tree (or possibly | |
117 | // keyed object in the case of PROOF) to be processed. The entry argument | |
118 | // specifies which entry in the currently loaded tree is to be processed. | |
119 | // It can be passed to either TTree::GetEntry() or TBranch::GetEntry() | |
120 | // to read either all or the required parts of the data. When processing | |
121 | // keyed objects with PROOF, the object is already loaded and is available | |
122 | // via the fObject pointer. | |
123 | // | |
124 | // This function should contain the "body" of the analysis. It can contain | |
125 | // simple or elaborate selection criteria, run algorithms on the data | |
126 | // of the event and typically fill histograms. | |
127 | ||
128 | // WARNING when a selector is used with a TChain, you must use | |
129 | // the pointer to the current TTree to call GetEntry(entry). | |
130 | // The entry is always the local entry number in the current tree. | |
131 | // Assuming that fTree is the pointer to the TChain being processed, | |
132 | // use fTree->GetTree()->GetEntry(entry). | |
133 | ||
134 | if (AliSelectorRL::Process(entry) == kFALSE) | |
135 | return kFALSE; | |
136 | ||
137 | // Check prerequisites | |
138 | if (!fESD) | |
139 | { | |
140 | AliDebug(AliLog::kError, "ESD branch not available"); | |
141 | return kFALSE; | |
142 | } | |
143 | ||
144 | if (!fEsdTrackCuts) | |
145 | { | |
146 | AliDebug(AliLog::kError, "fESDTrackCuts not available"); | |
147 | return kFALSE; | |
148 | } | |
149 | ||
150 | AliStack* stack = GetStack(); | |
151 | if (!stack) | |
152 | { | |
153 | AliDebug(AliLog::kError, "Stack not available"); | |
154 | return kFALSE; | |
155 | } | |
156 | ||
157 | TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD); | |
158 | ||
159 | if (fdNdEtaCorrection[0]) | |
160 | FillCorrectionMaps(list); | |
161 | ||
162 | if (fSecondaries) | |
163 | FillSecondaries(list); | |
164 | ||
9f469bf5 | 165 | if (fSigmaVertex) |
166 | FillSigmaVertex(); | |
167 | ||
10ebe68d | 168 | delete list; |
169 | list = 0; | |
170 | ||
171 | return kTRUE; | |
172 | } | |
173 | ||
174 | void AlidNdEtaSystematicsSelector::FillCorrectionMaps(TObjArray* listOfTracks) | |
175 | { | |
176 | // fills the correction maps for different particle species | |
177 | ||
178 | AliStack* stack = GetStack(); | |
179 | AliHeader* header = GetHeader(); | |
180 | ||
181 | Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD); | |
182 | Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD); | |
183 | ||
184 | for (Int_t i=0; i<4; ++i) | |
185 | { | |
186 | fdNdEtaCorrection[i]->IncreaseEventCount(); | |
187 | if (eventTriggered) | |
188 | fdNdEtaCorrection[i]->IncreaseTriggeredEventCount(); | |
189 | } | |
190 | ||
191 | // get the MC vertex | |
192 | AliGenEventHeader* genHeader = header->GenEventHeader(); | |
193 | ||
194 | TArrayF vtxMC(3); | |
195 | genHeader->PrimaryVertex(vtxMC); | |
196 | ||
197 | // loop over mc particles | |
198 | Int_t nPrim = stack->GetNprimary(); | |
199 | ||
200 | for (Int_t iMc = 0; iMc < nPrim; ++iMc) | |
201 | { | |
202 | TParticle* particle = stack->Particle(iMc); | |
203 | ||
204 | if (!particle) | |
205 | { | |
206 | AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc)); | |
207 | continue; | |
208 | } | |
209 | ||
210 | if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE) | |
211 | continue; | |
212 | ||
213 | Float_t eta = particle->Eta(); | |
214 | Float_t pt = particle->Pt(); | |
215 | ||
216 | Int_t id = -1; | |
217 | switch (TMath::Abs(particle->GetPdgCode())) | |
218 | { | |
219 | case 211: id = 0; break; | |
220 | case 321: id = 1; break; | |
221 | case 2212: id = 2; break; | |
61385583 | 222 | case 11: |
223 | { | |
224 | if (pt < 0.1 && particle->GetMother(0) > -1) | |
225 | { | |
226 | TParticle* mother = stack->Particle(particle->GetMother(0)); | |
227 | printf("Mother pdg code is %d\n", mother->GetPdgCode()); | |
228 | } | |
229 | //particle->Dump(); | |
230 | //if (particle->GetUniqueID() == 1) | |
231 | ||
232 | //printf("Found illegal particle mc id = %d file = %s event = %d\n", iMc, fTree->GetCurrentFile()->GetName(), fTree->GetTree()->GetReadEntry()); | |
233 | } | |
10ebe68d | 234 | default: id = 3; break; |
235 | } | |
236 | ||
237 | if (vertexReconstructed) | |
61385583 | 238 | { |
10ebe68d | 239 | fdNdEtaCorrection[id]->FillParticle(vtxMC[2], eta, pt); |
61385583 | 240 | if (pt < 0.1) |
241 | fPIDParticles->Fill(particle->GetPdgCode()); | |
242 | } | |
10ebe68d | 243 | |
244 | fdNdEtaCorrection[id]->FillParticleAllEvents(eta, pt); | |
245 | if (eventTriggered) | |
246 | fdNdEtaCorrection[id]->FillParticleWhenEventTriggered(eta, pt); | |
247 | }// end of mc particle | |
248 | ||
249 | // loop over esd tracks | |
250 | TIterator* iter = listOfTracks->MakeIterator(); | |
251 | TObject* obj = 0; | |
252 | while ((obj = iter->Next()) != 0) | |
253 | { | |
254 | AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj); | |
255 | if (!esdTrack) | |
256 | continue; | |
257 | ||
258 | // using the properties of the mc particle | |
259 | Int_t label = TMath::Abs(esdTrack->GetLabel()); | |
10ebe68d | 260 | TParticle* particle = stack->Particle(label); |
261 | if (!particle) | |
262 | { | |
263 | AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label)); | |
264 | continue; | |
265 | } | |
266 | ||
267 | Int_t id = -1; | |
268 | switch (TMath::Abs(particle->GetPdgCode())) | |
269 | { | |
9f469bf5 | 270 | case 211: id = 0; break; |
271 | case 321: id = 1; break; | |
10ebe68d | 272 | case 2212: id = 2; break; |
9f469bf5 | 273 | default: id = 3; break; |
10ebe68d | 274 | } |
275 | ||
276 | if (vertexReconstructed) | |
61385583 | 277 | { |
10ebe68d | 278 | fdNdEtaCorrection[id]->FillParticleWhenMeasuredTrack(vtxMC[2], particle->Eta(), particle->Pt()); |
61385583 | 279 | if (particle->Pt() < 0.1) |
280 | fPIDTracks->Fill(particle->GetPdgCode()); | |
281 | } | |
10ebe68d | 282 | } // end of track loop |
283 | ||
284 | delete iter; | |
285 | iter = 0; | |
286 | } | |
287 | ||
288 | void AlidNdEtaSystematicsSelector::FillSecondaries(TObjArray* listOfTracks) | |
289 | { | |
290 | // fills the secondary histograms | |
291 | ||
292 | AliStack* stack = GetStack(); | |
293 | ||
294 | TH1* nPrimaries = new TH1F("secondaries_primaries", "secondaries_primaries", fSecondaries->GetZaxis()->GetNbins(), fSecondaries->GetZaxis()->GetXmin(), fSecondaries->GetZaxis()->GetXmax()); | |
295 | TH1* nSecondaries = new TH1F("secondaries_secondaries", "secondaries_secondaries", fSecondaries->GetZaxis()->GetNbins(), fSecondaries->GetZaxis()->GetXmin(), fSecondaries->GetZaxis()->GetXmax()); | |
296 | ||
297 | TIterator* iter = listOfTracks->MakeIterator(); | |
298 | TObject* obj = 0; | |
299 | while ((obj = iter->Next()) != 0) | |
300 | { | |
301 | AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj); | |
302 | if (!esdTrack) | |
303 | continue; | |
304 | ||
305 | Int_t label = TMath::Abs(esdTrack->GetLabel()); | |
10ebe68d | 306 | TParticle* particle = stack->Particle(label); |
307 | if (!particle) | |
308 | { | |
309 | AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label)); | |
310 | continue; | |
311 | } | |
312 | ||
313 | Int_t nPrim = stack->GetNprimary(); | |
314 | if (label < nPrim) | |
315 | nPrimaries->Fill(particle->Pt()); | |
316 | else | |
317 | nSecondaries->Fill(particle->Pt()); | |
318 | } | |
319 | ||
320 | for (Int_t i=1; i<=nPrimaries->GetNbinsX(); ++i) | |
321 | { | |
322 | Int_t primaries = (Int_t) nPrimaries->GetBinContent(i); | |
323 | Int_t secondaries = (Int_t) nSecondaries->GetBinContent(i); | |
324 | ||
325 | if (primaries + secondaries > 0) | |
326 | { | |
327 | AliDebug(AliLog::kDebug, Form("The ratio between primaries and secondaries is %d:%d = %f", primaries, secondaries, ((secondaries > 0) ? (Double_t) primaries / secondaries : -1))); | |
328 | ||
329 | for (Double_t factor = 0.5; factor < 2.01; factor += 0.1) | |
9f469bf5 | 330 | { |
331 | Double_t nTracks = (Double_t) primaries + (Double_t) secondaries * factor; | |
332 | fSecondaries->Fill(nTracks, factor, nPrimaries->GetBinCenter(i)); | |
333 | //if (secondaries > 0) printf("We fill: %f %f %f\n", nTracks, factor, nPrimaries->GetBinCenter(i)); | |
334 | } | |
10ebe68d | 335 | } |
336 | } | |
337 | ||
f31d5d49 | 338 | fOverallPrimaries += (Int_t) nPrimaries->Integral(); |
339 | fOverallSecondaries += (Int_t) nSecondaries->Integral(); | |
340 | ||
10ebe68d | 341 | delete nPrimaries; |
342 | nPrimaries = 0; | |
343 | ||
344 | delete nSecondaries; | |
345 | nSecondaries = 0; | |
346 | ||
347 | delete iter; | |
348 | iter = 0; | |
349 | } | |
350 | ||
9f469bf5 | 351 | void AlidNdEtaSystematicsSelector::FillSigmaVertex() |
352 | { | |
353 | // fills the fSigmaVertex histogram | |
354 | ||
355 | // save the old value | |
356 | Float_t oldSigmaVertex = fEsdTrackCuts->GetMinNsigmaToVertex(); | |
357 | ||
358 | // set to maximum | |
359 | fEsdTrackCuts->SetMinNsigmaToVertex(5); | |
360 | ||
361 | TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD); | |
362 | ||
363 | TIterator* iter = list->MakeIterator(); | |
364 | TObject* obj = 0; | |
365 | while ((obj = iter->Next()) != 0) | |
366 | { | |
367 | AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj); | |
368 | if (!esdTrack) | |
369 | continue; | |
370 | ||
371 | Float_t sigma2Vertex = fEsdTrackCuts->GetSigmaToVertex(esdTrack); | |
372 | ||
373 | for (Double_t nSigma = 0.5; nSigma < 5.1; nSigma += 0.5) | |
374 | { | |
375 | if (sigma2Vertex < nSigma) | |
376 | fSigmaVertex->Fill(nSigma); | |
377 | } | |
378 | } | |
379 | ||
380 | delete iter; | |
381 | iter = 0; | |
382 | ||
383 | delete list; | |
384 | list = 0; | |
385 | ||
386 | // set back the old value | |
387 | fEsdTrackCuts->SetMinNsigmaToVertex(oldSigmaVertex); | |
388 | } | |
389 | ||
10ebe68d | 390 | void AlidNdEtaSystematicsSelector::SlaveTerminate() |
391 | { | |
392 | // The SlaveTerminate() function is called after all entries or objects | |
393 | // have been processed. When running with PROOF SlaveTerminate() is called | |
394 | // on each slave server. | |
395 | ||
396 | AliSelector::SlaveTerminate(); | |
397 | ||
398 | // Add the histograms to the output on each slave server | |
399 | if (!fOutput) | |
400 | { | |
401 | AliDebug(AliLog::kError, Form("ERROR: Output list not initialized.")); | |
402 | return; | |
403 | } | |
404 | ||
405 | if (fSecondaries) | |
406 | fOutput->Add(fSecondaries); | |
407 | ||
408 | for (Int_t i=0; i<4; ++i) | |
409 | if (fdNdEtaCorrection[i]) | |
410 | fOutput->Add(fdNdEtaCorrection[i]); | |
9f469bf5 | 411 | |
412 | if (fSigmaVertex) | |
413 | fOutput->Add(fSigmaVertex); | |
10ebe68d | 414 | } |
415 | ||
416 | void AlidNdEtaSystematicsSelector::Terminate() | |
417 | { | |
418 | // The Terminate() function is the last function to be called during | |
419 | // a query. It always runs on the client, it can be used to present | |
420 | // the results graphically or save the results to file. | |
421 | ||
422 | AliSelector::Terminate(); | |
423 | ||
424 | fSecondaries = dynamic_cast<TH3F*> (fOutput->FindObject("fSecondaries")); | |
425 | for (Int_t i=0; i<4; ++i) | |
426 | fdNdEtaCorrection[i] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject(Form("correction_%d", i))); | |
9f469bf5 | 427 | fSigmaVertex = dynamic_cast<TH1F*> (fOutput->FindObject("fSigmaVertex")); |
10ebe68d | 428 | |
61385583 | 429 | TDatabasePDG* pdgDB = new TDatabasePDG; |
430 | ||
431 | for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i) | |
432 | if (fPIDParticles->GetBinContent(i) > 0) | |
433 | printf("PDG = %d (%s): generated: %d, reconstructed: %d, ratio: %f\n", (Int_t) fPIDParticles->GetBinCenter(i), pdgDB->GetParticle((Int_t) fPIDParticles->GetBinCenter(i))->GetName(), (Int_t) fPIDParticles->GetBinContent(i), (Int_t) fPIDTracks->GetBinContent(i), ((fPIDTracks->GetBinContent(i) > 0) ? fPIDParticles->GetBinContent(i) / fPIDTracks->GetBinContent(i) : -1)); | |
434 | ||
435 | delete pdgDB; | |
436 | pdgDB = 0; | |
437 | ||
10ebe68d | 438 | TFile* fout = TFile::Open("systematics.root", "RECREATE"); |
439 | ||
440 | if (fEsdTrackCuts) | |
441 | fEsdTrackCuts->SaveHistograms("esd_track_cuts"); | |
442 | ||
443 | if (fSecondaries) | |
f31d5d49 | 444 | { |
10ebe68d | 445 | fSecondaries->Write(); |
f31d5d49 | 446 | printf("We had %d primaries and %d secondaries.\n", (Int_t) fOverallPrimaries, (Int_t) fOverallSecondaries); |
447 | } | |
10ebe68d | 448 | |
9f469bf5 | 449 | if (fSigmaVertex) |
450 | fSigmaVertex->Write(); | |
451 | ||
10ebe68d | 452 | for (Int_t i=0; i<4; ++i) |
453 | if (fdNdEtaCorrection[i]) | |
454 | fdNdEtaCorrection[i]->SaveHistograms(); | |
455 | ||
456 | fout->Write(); | |
457 | fout->Close(); | |
458 | ||
459 | if (fSecondaries) | |
460 | { | |
461 | new TCanvas; | |
462 | fSecondaries->Draw(); | |
463 | } | |
9f469bf5 | 464 | |
465 | if (fSigmaVertex) | |
466 | { | |
467 | new TCanvas; | |
468 | fSigmaVertex->Draw(); | |
469 | } | |
10ebe68d | 470 | } |