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