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