Record changes.
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AlidNdEtaSystematicsSelector.cxx
CommitLineData
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>
7af955da 11#include <TH2F.h>
10ebe68d 12#include <TH1F.h>
13#include <TParticle.h>
38faa1b9 14#include <TParticlePDG.h>
10ebe68d 15
16#include <AliLog.h>
17#include <AliESD.h>
10ebe68d 18#include <AliStack.h>
19#include <AliHeader.h>
6de7db91 20#include <AliGenEventHeader.h>
4c351225 21#include <../PYTHIA6/AliGenPythiaEventHeader.h>
22#include <../EVGEN/AliGenCocktailEventHeader.h>
10ebe68d 23
24#include "esdTrackCuts/AliESDtrackCuts.h"
25#include "AliPWG0Helper.h"
4c351225 26#include "AliPWG0depHelper.h"
0bd1f8a0 27#include "dNdEta/AlidNdEtaCorrection.h"
10ebe68d 28
29ClassImp(AlidNdEtaSystematicsSelector)
30
31AlidNdEtaSystematicsSelector::AlidNdEtaSystematicsSelector() :
32 AliSelectorRL(),
33 fSecondaries(0),
9f469bf5 34 fSigmaVertex(0),
f31d5d49 35 fEsdTrackCuts(0),
61385583 36 fPIDParticles(0),
38faa1b9 37 fPIDTracks(0),
38 fSignMode(0),
39 fMultiplicityMode(0)
10ebe68d 40{
41 //
42 // Constructor. Initialization of pointers
43 //
44
45 for (Int_t i=0; i<4; ++i)
6de7db91 46 fdNdEtaCorrectionSpecies[i] = 0;
47
7b0956f3 48 for (Int_t i=0; i<3; ++i) {
6de7db91 49 fdNdEtaCorrectionVertexReco[i] = 0;
7b0956f3 50 fdNdEtaCorrectionTriggerBias[i] = 0;
51 }
10ebe68d 52}
53
54AlidNdEtaSystematicsSelector::~AlidNdEtaSystematicsSelector()
55{
56 //
57 // Destructor
58 //
59}
60
61void AlidNdEtaSystematicsSelector::Begin(TTree* tree)
62{
63 // Begin function
64
65 AliSelectorRL::Begin(tree);
66
67 ReadUserObjects(tree);
68}
69
70void AlidNdEtaSystematicsSelector::ReadUserObjects(TTree* tree)
71{
72 // read the user objects, called from slavebegin and begin
73
74 if (!fEsdTrackCuts && fInput)
75 fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fInput->FindObject("AliESDtrackCuts"));
76
77 if (!fEsdTrackCuts && tree)
78 fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (tree->GetUserInfo()->FindObject("AliESDtrackCuts"));
79
80 if (!fEsdTrackCuts)
81 AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from input list.");
82}
83
84void AlidNdEtaSystematicsSelector::SlaveBegin(TTree* tree)
85{
86 // The SlaveBegin() function is called after the Begin() function.
87 // When running with PROOF SlaveBegin() is called on each slave server.
88 // The tree argument is deprecated (on PROOF 0 is passed).
89
90 AliSelector::SlaveBegin(tree);
91
92 ReadUserObjects(tree);
93
94 TString option(GetOption());
95
96 printf("Running AlidNdEtaSystematicsSelector with options %s\n", option.Data());
97
98 if (option.Contains("secondaries"))
99 {
7af955da 100 fSecondaries = new TH2F("fSecondaries", "fSecondaries;Case;N", 8, 0.5, 8.5, 2001, -0.5, 2000.5);
101 fSecondaries->GetXaxis()->SetBinLabel(1, "All Primaries");
102 fSecondaries->GetXaxis()->SetBinLabel(2, "All Secondaries");
103 fSecondaries->GetXaxis()->SetBinLabel(3, "Primaries pT > 0.3 GeV/c");
104 fSecondaries->GetXaxis()->SetBinLabel(4, "Secondaries pT > 0.3 GeV/c");
105 fSecondaries->GetXaxis()->SetBinLabel(5, "Tracks from Primaries");
106 fSecondaries->GetXaxis()->SetBinLabel(6, "Tracks from Secondaries");
107 fSecondaries->GetXaxis()->SetBinLabel(7, "Accepted Tracks from Primaries");
108 fSecondaries->GetXaxis()->SetBinLabel(8, "Accepted Tracks from Secondaries");
10ebe68d 109 }
110
111 if (option.Contains("particle-composition"))
112 {
113 for (Int_t i=0; i<4; ++i)
114 {
115 TString name;
116 name.Form("correction_%d", i);
6de7db91 117 fdNdEtaCorrectionSpecies[i] = new AlidNdEtaCorrection(name, name);
10ebe68d 118 }
119 }
9f469bf5 120
121 if (option.Contains("sigma-vertex"))
122 {
72e597d7 123 fSigmaVertex = new TH1F("fSigmaVertex", "fSigmaVertex;Nsigma2vertex;NacceptedTracks", 50, 0.05, 5.05);
9f469bf5 124 printf("WARNING: sigma-vertex analysis enabled. This will produce weird results in the AliESDtrackCuts histograms\n");
125 }
61385583 126
6de7db91 127 if (option.Contains("vertexreco")) {
128 fdNdEtaCorrectionVertexReco[0] = new AlidNdEtaCorrection("vertexRecoND", "vertexRecoND");
129 fdNdEtaCorrectionVertexReco[1] = new AlidNdEtaCorrection("vertexRecoSD", "vertexRecoSD");
130 fdNdEtaCorrectionVertexReco[2] = new AlidNdEtaCorrection("vertexRecoDD", "vertexRecoDD");
131 }
132
7b0956f3 133 if (option.Contains("triggerbias")) {
134 fdNdEtaCorrectionTriggerBias[0] = new AlidNdEtaCorrection("triggerBiasND", "triggerBiasND");
135 fdNdEtaCorrectionTriggerBias[1] = new AlidNdEtaCorrection("triggerBiasSD", "triggerBiasSD");
136 fdNdEtaCorrectionTriggerBias[2] = new AlidNdEtaCorrection("triggerBiasDD", "triggerBiasDD");
137 }
138
61385583 139 fPIDParticles = new TH1F("pid_particles", "PID of generated primary particles", 10001, -5000.5, 5000.5);
140 fPIDTracks = new TH1F("pid_tracks", "MC PID of reconstructed tracks", 10001, -5000.5, 5000.5);
38faa1b9 141
142 if (option.Contains("only-positive"))
143 {
144 AliInfo("Processing only positive particles.");
145 fSignMode = 1;
146 }
147 else if (option.Contains("only-negative"))
148 {
149 AliInfo("Processing only negative particles.");
150 fSignMode = -1;
151 }
152
153 if (option.Contains("low-multiplicity"))
154 {
155 AliInfo("Processing only events with low multiplicity.");
156 fMultiplicityMode = 1;
157 }
158 else if (option.Contains("high-multiplicity"))
159 {
160 AliInfo("Processing only events with high multiplicity.");
161 fMultiplicityMode = 2;
162 }
10ebe68d 163}
164
165Bool_t AlidNdEtaSystematicsSelector::Process(Long64_t entry)
166{
167 // The Process() function is called for each entry in the tree (or possibly
168 // keyed object in the case of PROOF) to be processed. The entry argument
169 // specifies which entry in the currently loaded tree is to be processed.
170 // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
171 // to read either all or the required parts of the data. When processing
172 // keyed objects with PROOF, the object is already loaded and is available
173 // via the fObject pointer.
174 //
175 // This function should contain the "body" of the analysis. It can contain
176 // simple or elaborate selection criteria, run algorithms on the data
177 // of the event and typically fill histograms.
178
179 // WARNING when a selector is used with a TChain, you must use
180 // the pointer to the current TTree to call GetEntry(entry).
181 // The entry is always the local entry number in the current tree.
182 // Assuming that fTree is the pointer to the TChain being processed,
183 // use fTree->GetTree()->GetEntry(entry).
184
185 if (AliSelectorRL::Process(entry) == kFALSE)
186 return kFALSE;
187
188 // Check prerequisites
189 if (!fESD)
190 {
191 AliDebug(AliLog::kError, "ESD branch not available");
192 return kFALSE;
193 }
194
195 if (!fEsdTrackCuts)
196 {
197 AliDebug(AliLog::kError, "fESDTrackCuts not available");
198 return kFALSE;
199 }
200
201 AliStack* stack = GetStack();
202 if (!stack)
203 {
204 AliDebug(AliLog::kError, "Stack not available");
205 return kFALSE;
206 }
207
208 TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
209
38faa1b9 210 if (fMultiplicityMode == 1 && list->GetEntries() > 20 ||
211 fMultiplicityMode == 2 && list->GetEntries() < 40)
212 {
213 delete list;
214 list = 0;
215 return kTRUE;
216 }
217
6de7db91 218 if (fdNdEtaCorrectionSpecies[0])
10ebe68d 219 FillCorrectionMaps(list);
220
221 if (fSecondaries)
7af955da 222 FillSecondaries();
10ebe68d 223
9f469bf5 224 if (fSigmaVertex)
225 FillSigmaVertex();
226
6de7db91 227 // stuff for vertex reconstruction correction systematics
7b0956f3 228 Bool_t vertexRecoStudy = kFALSE;
229 Bool_t triggerBiasStudy = kFALSE;
230 if (fdNdEtaCorrectionVertexReco[0]) vertexRecoStudy = kTRUE;
231 if (fdNdEtaCorrectionTriggerBias[0]) triggerBiasStudy = kTRUE;
232
233 if (vertexRecoStudy || triggerBiasStudy) {
6de7db91 234 AliHeader* header = GetHeader();
235 if (!header) {
236 AliDebug(AliLog::kError, "Header not available");
237 return kFALSE;
238 }
239
7b0956f3 240 // getting process information
4c351225 241 Int_t processtype = AliPWG0depHelper::GetPythiaEventProcessType(header);
6de7db91 242
7b0956f3 243 AliDebug(AliLog::kInfo, Form("Pythia process type %d.",processtype));
244
245 // can only read pythia headers, either directly or from cocktalil header
246 AliGenEventHeader* genHeader = (AliGenEventHeader*)(header->GenEventHeader());
6de7db91 247
7b0956f3 248 if (!genHeader) {
249 AliDebug(AliLog::kError, "Gen header not available");
250 return kFALSE;
251 }
6de7db91 252
253 // get the MC vertex
254 TArrayF vtxMC(3);
255 genHeader->PrimaryVertex(vtxMC);
256
257 Int_t nGoodTracks = list->GetEntries();
258
259 Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
260 Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD);
261
6de7db91 262 // non diffractive
263 if (processtype!=92 && processtype!=93 && processtype!=94) {
7b0956f3 264 if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[0]->FillEventAll(vtxMC[2], nGoodTracks);
265
6de7db91 266 if (eventTriggered) {
4c351225 267 if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[0]->FillEventWithTrigger(vtxMC[2], nGoodTracks);
268 if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[0]->FillEventWithTrigger(vtxMC[2], nGoodTracks);
7b0956f3 269
4c351225 270 if (vertexReconstructed)
271 if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[0]->FillEventWithTriggerWithReconstructedVertex(vtxMC[2], nGoodTracks);
6de7db91 272 }
273 }
274
275 // single diffractive
276 if (processtype==92 || processtype==93) {
7b0956f3 277 if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[1]->FillEventAll(vtxMC[2], nGoodTracks);
278
6de7db91 279 if (eventTriggered) {
4c351225 280 if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[1]->FillEventWithTrigger(vtxMC[2], nGoodTracks);
281 if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[1]->FillEventWithTrigger(vtxMC[2], nGoodTracks);
7b0956f3 282
4c351225 283 if (vertexReconstructed)
284 if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[1]->FillEventWithTriggerWithReconstructedVertex(vtxMC[2], nGoodTracks);
6de7db91 285 }
286 }
287
288 // double diffractive
289 if (processtype==94) {
7b0956f3 290 if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[2]->FillEventAll(vtxMC[2], nGoodTracks);
291
6de7db91 292 if (eventTriggered) {
4c351225 293 if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[2]->FillEventWithTrigger(vtxMC[2], nGoodTracks);
294 if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[2]->FillEventWithTrigger(vtxMC[2], nGoodTracks);
7b0956f3 295
4c351225 296 if (vertexReconstructed)
297 if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[2]->FillEventWithTriggerWithReconstructedVertex(vtxMC[2], nGoodTracks);
6de7db91 298 }
299 }
300 }
301
10ebe68d 302 delete list;
303 list = 0;
304
305 return kTRUE;
306}
307
308void AlidNdEtaSystematicsSelector::FillCorrectionMaps(TObjArray* listOfTracks)
309{
310 // fills the correction maps for different particle species
311
38faa1b9 312 // TODO fix the use of the FillParticle* functions
313
10ebe68d 314 AliStack* stack = GetStack();
315 AliHeader* header = GetHeader();
316
317 Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD);
318 Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
319
10ebe68d 320 // get the MC vertex
321 AliGenEventHeader* genHeader = header->GenEventHeader();
322
323 TArrayF vtxMC(3);
324 genHeader->PrimaryVertex(vtxMC);
325
326 // loop over mc particles
327 Int_t nPrim = stack->GetNprimary();
328
329 for (Int_t iMc = 0; iMc < nPrim; ++iMc)
330 {
331 TParticle* particle = stack->Particle(iMc);
332
333 if (!particle)
334 {
335 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
336 continue;
337 }
338
339 if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
340 continue;
341
38faa1b9 342 if (SignOK(particle->GetPDG()) == kFALSE)
343 continue;
344
10ebe68d 345 Float_t eta = particle->Eta();
346 Float_t pt = particle->Pt();
347
348 Int_t id = -1;
349 switch (TMath::Abs(particle->GetPdgCode()))
350 {
351 case 211: id = 0; break;
352 case 321: id = 1; break;
353 case 2212: id = 2; break;
7af955da 354 case 11:
61385583 355 {
7af955da 356 /*if (pt < 0.1 && particle->GetMother(0) > -1)
61385583 357 {
358 TParticle* mother = stack->Particle(particle->GetMother(0));
359 printf("Mother pdg code is %d\n", mother->GetPdgCode());
7af955da 360 } */
61385583 361 //particle->Dump();
362 //if (particle->GetUniqueID() == 1)
363
7af955da 364 //printf("Found illegal particle mc id = %d file = %s event = %d\n", iMc, fTree->GetCurrentFile()->GetName(), fTree->GetTree()->GetReadEntry());
61385583 365 }
10ebe68d 366 default: id = 3; break;
367 }
368
38faa1b9 369 if (vertexReconstructed && eventTriggered)
61385583 370 {
38faa1b9 371 fdNdEtaCorrectionSpecies[id]->FillParticleVertex(vtxMC[2], eta, pt);
7af955da 372 //if (pt < 0.1)
61385583 373 fPIDParticles->Fill(particle->GetPdgCode());
374 }
10ebe68d 375
7b0956f3 376 //fdNdEtaCorrectionSpecies[id]->FillParticleAllEvents(eta, pt);
377 //if (eventTriggered)
378 // fdNdEtaCorrectionSpecies[id]->FillParticleWhenEventTriggered(eta, pt);
10ebe68d 379 }// end of mc particle
380
381 // loop over esd tracks
382 TIterator* iter = listOfTracks->MakeIterator();
383 TObject* obj = 0;
384 while ((obj = iter->Next()) != 0)
385 {
386 AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj);
387 if (!esdTrack)
388 continue;
389
390 // using the properties of the mc particle
391 Int_t label = TMath::Abs(esdTrack->GetLabel());
10ebe68d 392 TParticle* particle = stack->Particle(label);
393 if (!particle)
394 {
395 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
396 continue;
397 }
398
7af955da 399 TParticle* mother = particle;
400 // find primary particle that created this particle
401 while (label >= nPrim)
402 {
403 //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
404
405 if (mother->GetMother(0) == -1)
406 {
407 AliDebug(AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
408 mother = 0;
409 break;
410 }
411
412 label = mother->GetMother(0);
413
414 mother = stack->Particle(label);
415 if (!mother)
416 {
417 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (find mother loop).", label));
418 break;
419 }
420 }
421
422 if (!mother)
423 continue;
424
38faa1b9 425 if (SignOK(mother->GetPDG()) == kFALSE)
426 continue;
427
7af955da 428 //printf("We continue with particle %d (pdg %d)\n", label, mother->GetPdgCode());
429
10ebe68d 430 Int_t id = -1;
7af955da 431 switch (TMath::Abs(mother->GetPdgCode()))
10ebe68d 432 {
9f469bf5 433 case 211: id = 0; break;
434 case 321: id = 1; break;
10ebe68d 435 case 2212: id = 2; break;
9f469bf5 436 default: id = 3; break;
10ebe68d 437 }
438
38faa1b9 439 if (vertexReconstructed && eventTriggered)
61385583 440 {
38faa1b9 441 fdNdEtaCorrectionSpecies[id]->FillParticleTracked(vtxMC[2], particle->Eta(), particle->Pt());
7af955da 442 //if (particle->Pt() < 0.1)
61385583 443 fPIDTracks->Fill(particle->GetPdgCode());
444 }
10ebe68d 445 } // end of track loop
446
7af955da 447 Int_t nGoodTracks = listOfTracks->GetEntries();
448
449 for (Int_t i=0; i<4; ++i)
450 {
7b0956f3 451 fdNdEtaCorrectionSpecies[i]->FillEventAll(vtxMC[2], nGoodTracks);
7af955da 452 if (eventTriggered)
453 {
6de7db91 454 fdNdEtaCorrectionSpecies[i]->FillEventWithTrigger(vtxMC[2], nGoodTracks);
7af955da 455 if (vertexReconstructed)
6de7db91 456 fdNdEtaCorrectionSpecies[i]->FillEventWithTriggerWithReconstructedVertex(vtxMC[2], nGoodTracks);
7af955da 457 }
458 }
459
10ebe68d 460 delete iter;
461 iter = 0;
462}
463
7af955da 464void AlidNdEtaSystematicsSelector::FillSecondaries()
10ebe68d 465{
466 // fills the secondary histograms
467
468 AliStack* stack = GetStack();
469
7af955da 470 Int_t particlesPrimaries = 0;
471 Int_t particlesSecondaries = 0;
472 Int_t particlesPrimariesPtCut = 0;
473 Int_t particlesSecondariesPtCut = 0;
10ebe68d 474
7af955da 475 for (Int_t iParticle = 0; iParticle < stack->GetNtrack(); iParticle++)
10ebe68d 476 {
7af955da 477 TParticle* particle = stack->Particle(iParticle);
478 if (!particle)
479 {
480 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (particle loop).", iParticle));
481 continue;
482 }
483
484 if (TMath::Abs(particle->Eta()) > 0.9)
485 continue;
486
487 Bool_t isPrimary = kFALSE;
488 if (iParticle < stack->GetNprimary())
489 {
490 if (AliPWG0Helper::IsPrimaryCharged(particle, stack->GetNprimary()) == kFALSE)
491 continue;
492
493 isPrimary = kTRUE;
494 }
495
496 if (isPrimary)
497 particlesPrimaries++;
498 else
499 particlesSecondaries++;
500
501 if (particle->Pt() < 0.3)
502 continue;
503
504 if (isPrimary)
505 particlesPrimariesPtCut++;
506 else
507 particlesSecondariesPtCut++;
508 }
509
510 fSecondaries->Fill(1, particlesPrimaries);
511 fSecondaries->Fill(2, particlesSecondaries);
512 fSecondaries->Fill(3, particlesPrimariesPtCut);
513 fSecondaries->Fill(4, particlesSecondariesPtCut);
514
515 Int_t allTracksPrimaries = 0;
516 Int_t allTracksSecondaries = 0;
517 Int_t acceptedTracksPrimaries = 0;
518 Int_t acceptedTracksSecondaries = 0;
519
520 for (Int_t iTrack = 0; iTrack < fESD->GetNumberOfTracks(); iTrack++)
521 {
522 AliESDtrack* esdTrack = fESD->GetTrack(iTrack);
10ebe68d 523 if (!esdTrack)
7af955da 524 {
525 AliDebug(AliLog::kError, Form("UNEXPECTED: Could not get track %d.", iTrack));
10ebe68d 526 continue;
7af955da 527 }
10ebe68d 528
529 Int_t label = TMath::Abs(esdTrack->GetLabel());
10ebe68d 530 TParticle* particle = stack->Particle(label);
531 if (!particle)
532 {
533 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
534 continue;
535 }
536
7af955da 537 if (label < stack->GetNprimary())
538 allTracksPrimaries++;
10ebe68d 539 else
7af955da 540 allTracksSecondaries++;
10ebe68d 541
7af955da 542 if (!fEsdTrackCuts->AcceptTrack(esdTrack))
543 continue;
10ebe68d 544
7af955da 545 if (label < stack->GetNprimary())
546 acceptedTracksPrimaries++;
547 else
548 acceptedTracksSecondaries++;
10ebe68d 549 }
550
7af955da 551 fSecondaries->Fill(5, allTracksPrimaries);
552 fSecondaries->Fill(6, allTracksSecondaries);
553 fSecondaries->Fill(7, acceptedTracksPrimaries);
554 fSecondaries->Fill(8, acceptedTracksSecondaries);
10ebe68d 555
7af955da 556 //printf("P = %d, S = %d, P_pt = %d, S_pt = %d, T_P = %d, T_S = %d, T_P_acc = %d, T_S_acc = %d\n", particlesPrimaries, particlesSecondaries, particlesPrimariesPtCut, particlesSecondariesPtCut, allTracksPrimaries, allTracksSecondaries, acceptedTracksPrimaries, acceptedTracksSecondaries);
10ebe68d 557}
558
9f469bf5 559void AlidNdEtaSystematicsSelector::FillSigmaVertex()
560{
561 // fills the fSigmaVertex histogram
562
563 // save the old value
564 Float_t oldSigmaVertex = fEsdTrackCuts->GetMinNsigmaToVertex();
565
566 // set to maximum
567 fEsdTrackCuts->SetMinNsigmaToVertex(5);
568
569 TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
570
571 TIterator* iter = list->MakeIterator();
572 TObject* obj = 0;
573 while ((obj = iter->Next()) != 0)
574 {
575 AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj);
576 if (!esdTrack)
577 continue;
578
579 Float_t sigma2Vertex = fEsdTrackCuts->GetSigmaToVertex(esdTrack);
580
72e597d7 581 for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
9f469bf5 582 {
583 if (sigma2Vertex < nSigma)
584 fSigmaVertex->Fill(nSigma);
585 }
586 }
587
588 delete iter;
589 iter = 0;
590
591 delete list;
592 list = 0;
593
594 // set back the old value
595 fEsdTrackCuts->SetMinNsigmaToVertex(oldSigmaVertex);
596}
597
10ebe68d 598void AlidNdEtaSystematicsSelector::SlaveTerminate()
599{
600 // The SlaveTerminate() function is called after all entries or objects
601 // have been processed. When running with PROOF SlaveTerminate() is called
602 // on each slave server.
603
604 AliSelector::SlaveTerminate();
605
606 // Add the histograms to the output on each slave server
607 if (!fOutput)
608 {
609 AliDebug(AliLog::kError, Form("ERROR: Output list not initialized."));
610 return;
611 }
612
613 if (fSecondaries)
614 fOutput->Add(fSecondaries);
615
616 for (Int_t i=0; i<4; ++i)
6de7db91 617 if (fdNdEtaCorrectionSpecies[i])
618 fOutput->Add(fdNdEtaCorrectionSpecies[i]);
9f469bf5 619
620 if (fSigmaVertex)
621 fOutput->Add(fSigmaVertex);
6de7db91 622
7b0956f3 623 for (Int_t i=0; i<3; ++i) {
6de7db91 624 if (fdNdEtaCorrectionVertexReco[i])
625 fOutput->Add(fdNdEtaCorrectionVertexReco[i]);
7b0956f3 626
627 if (fdNdEtaCorrectionTriggerBias[i])
628 fOutput->Add(fdNdEtaCorrectionTriggerBias[i]);
629 }
630
10ebe68d 631}
632
633void AlidNdEtaSystematicsSelector::Terminate()
634{
635 // The Terminate() function is the last function to be called during
636 // a query. It always runs on the client, it can be used to present
637 // the results graphically or save the results to file.
638
639 AliSelector::Terminate();
640
7af955da 641 fSecondaries = dynamic_cast<TH2F*> (fOutput->FindObject("fSecondaries"));
10ebe68d 642 for (Int_t i=0; i<4; ++i)
6de7db91 643 fdNdEtaCorrectionSpecies[i] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
9f469bf5 644 fSigmaVertex = dynamic_cast<TH1F*> (fOutput->FindObject("fSigmaVertex"));
10ebe68d 645
4c351225 646 fdNdEtaCorrectionVertexReco[0] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("vertexRecoND"));
647 fdNdEtaCorrectionVertexReco[1] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("vertexRecoSD"));
648 fdNdEtaCorrectionVertexReco[2] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("vertexRecoDD"));
649
650 fdNdEtaCorrectionTriggerBias[0] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("triggerBiasND"));
651 fdNdEtaCorrectionTriggerBias[1] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("triggerBiasSD"));
652 fdNdEtaCorrectionTriggerBias[2] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("triggerBiasDD"));
653
0bd1f8a0 654 if (fPIDParticles)
655 {
656 TDatabasePDG* pdgDB = new TDatabasePDG;
61385583 657
0bd1f8a0 658 for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i)
659 if (fPIDParticles->GetBinContent(i) > 0)
6de7db91 660 printf("PDG = %d (%s): generated: %d, reconstructed: %d, ratio: %f\n",
661 (Int_t) fPIDParticles->GetBinCenter(i),
662 pdgDB->GetParticle((Int_t) fPIDParticles->GetBinCenter(i))->GetName(),
663 (Int_t) fPIDParticles->GetBinContent(i),
664 (Int_t) fPIDTracks->GetBinContent(i),
665 ((fPIDTracks->GetBinContent(i) > 0) ? fPIDParticles->GetBinContent(i) / fPIDTracks->GetBinContent(i) : -1));
61385583 666
0bd1f8a0 667 delete pdgDB;
668 pdgDB = 0;
669 }
61385583 670
7b0956f3 671 for (Int_t i=0; i<3; ++i) {
6de7db91 672 if (fdNdEtaCorrectionVertexReco[i])
673 fdNdEtaCorrectionVertexReco[i]->Finish();
7b0956f3 674
675 if (fdNdEtaCorrectionTriggerBias[i])
676 fdNdEtaCorrectionTriggerBias[i]->Finish();
677 }
6de7db91 678
10ebe68d 679 TFile* fout = TFile::Open("systematics.root", "RECREATE");
680
681 if (fEsdTrackCuts)
682 fEsdTrackCuts->SaveHistograms("esd_track_cuts");
683
684 if (fSecondaries)
685 fSecondaries->Write();
686
9f469bf5 687 if (fSigmaVertex)
688 fSigmaVertex->Write();
7b0956f3 689
10ebe68d 690 for (Int_t i=0; i<4; ++i)
6de7db91 691 if (fdNdEtaCorrectionSpecies[i])
692 fdNdEtaCorrectionSpecies[i]->SaveHistograms();
693
7b0956f3 694 for (Int_t i=0; i<3; ++i) {
6de7db91 695 if (fdNdEtaCorrectionVertexReco[i])
696 fdNdEtaCorrectionVertexReco[i]->SaveHistograms();
10ebe68d 697
7b0956f3 698 if (fdNdEtaCorrectionTriggerBias[i])
699 fdNdEtaCorrectionTriggerBias[i]->SaveHistograms();
700 }
701
10ebe68d 702 fout->Write();
703 fout->Close();
704
705 if (fSecondaries)
706 {
707 new TCanvas;
7af955da 708 fSecondaries->Draw("COLZ");
10ebe68d 709 }
9f469bf5 710
711 if (fSigmaVertex)
712 {
713 new TCanvas;
714 fSigmaVertex->Draw();
715 }
10ebe68d 716}
38faa1b9 717
718Bool_t AlidNdEtaSystematicsSelector::SignOK(TParticlePDG* particle)
719{
720 // returns if a particle with this sign should be counted
721 // this is determined by the value of fSignMode, which should have the same sign
722 // as the charge
723 // if fSignMode is 0 all particles are counted
724
725 if (fSignMode == 0)
726 return kTRUE;
727
728 if (!particle)
729 {
730 printf("WARNING: not counting a particle that does not have a pdg particle\n");
731 return kFALSE;
732 }
733
734 Double_t charge = particle->Charge();
735
736 if (fSignMode > 0)
737 if (charge < 0)
738 return kFALSE;
739
740 if (fSignMode < 0)
741 if (charge > 0)
742 return kFALSE;
743
744 return kTRUE;
745}