Update
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AlidNdEtaSystematicsSelector.cxx
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 <TH2F.h>
12 #include <TH1F.h>
13 #include <TParticle.h>
14 #include <TParticlePDG.h>
15
16 #include <AliLog.h>
17 #include <AliESD.h>
18 #include <AliStack.h>
19 #include <AliHeader.h>
20 #include <AliGenEventHeader.h>
21 #include <../STEER/AliGenPythiaEventHeader.h>
22 #include <../STEER/AliGenCocktailEventHeader.h>
23
24 #include "AliESDtrackCuts.h"
25 #include "AliPWG0Helper.h"
26 #include "dNdEta/AlidNdEtaCorrection.h"
27
28 ClassImp(AlidNdEtaSystematicsSelector)
29
30 AlidNdEtaSystematicsSelector::AlidNdEtaSystematicsSelector() :
31   AliSelectorRL(),
32   fSecondaries(0),
33   fSigmaVertex(0),
34   fEsdTrackCuts(0),
35   fPIDParticles(0),
36   fPIDTracks(0),
37   fSignMode(0),
38   fMultiplicityMode(0)
39 {
40   //
41   // Constructor. Initialization of pointers
42   //
43
44   for (Int_t i=0; i<4; ++i)
45     fdNdEtaCorrectionSpecies[i] = 0;
46
47   for (Int_t i=0; i<3; ++i) {
48     fdNdEtaCorrectionVertexReco[i] = 0;
49     fdNdEtaCorrectionTriggerBias[i] = 0;
50   }
51 }
52
53 AlidNdEtaSystematicsSelector::~AlidNdEtaSystematicsSelector()
54 {
55   //
56   // Destructor
57   //
58 }
59
60 void AlidNdEtaSystematicsSelector::Begin(TTree* tree)
61 {
62   // Begin function
63
64   AliSelectorRL::Begin(tree);
65
66   ReadUserObjects(tree);
67 }
68
69 void 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
83 void 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
97   // Options: secondaries particle-composition sigma-vertex vertexreco triggerbias
98
99   if (option.Contains("secondaries"))
100   {
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");
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);
118       fdNdEtaCorrectionSpecies[i] = new AlidNdEtaCorrection(name, name);
119     }
120   }
121
122   if (option.Contains("sigma-vertex"))
123   {
124     fSigmaVertex = new TH1F("fSigmaVertex", "fSigmaVertex;Nsigma2vertex;NacceptedTracks", 50, 0.05, 5.05);
125     printf("WARNING: sigma-vertex analysis enabled. This will produce weird results in the AliESDtrackCuts histograms\n");
126   }
127
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
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
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);
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   }
164 }
165
166 Bool_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
209   // --------------------------------------------------------------
210   // related to events:
211   //
212   // stuff for vertex reconstruction correction systematics
213   Bool_t vertexRecoStudy  = kFALSE;
214   Bool_t triggerBiasStudy = kFALSE;
215   if (fdNdEtaCorrectionVertexReco[0])  vertexRecoStudy = kTRUE;
216   if (fdNdEtaCorrectionTriggerBias[0]) triggerBiasStudy = kTRUE;
217
218   AliHeader* header = GetHeader();
219   if (!header) {
220     AliDebug(AliLog::kError, "Header not available");
221     return kFALSE;
222   }
223   
224   // getting process information
225   Int_t processtype = AliPWG0Helper::GetEventProcessType(header);
226
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;
233   }
234   
235   // get the MC vertex
236   TArrayF vtxMC(3);
237   genHeader->PrimaryVertex(vtxMC);
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
246   if (processtype == AliPWG0Helper::kND) {
247     // NB: passing the wrong process type here (1), since the process type is defined by the index in the array (here non-diffractive) // ???? CKB ???
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 == AliPWG0Helper::kSD) { 
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==AliPWG0Helper::kDD) { 
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
264   for (Int_t i=0; i<4; ++i) {
265     if (fdNdEtaCorrectionSpecies[i])
266       fdNdEtaCorrectionSpecies[i]->FillEvent(vtxMC[2], nGoodTracks, eventTriggered, vertexReconstructed, 1);
267   }
268
269   // --------------------------------------------------------------
270   // MC particle loop
271   //
272
273   Int_t nPrim  = stack->GetNprimary();
274
275   for (Int_t iMc = 0; iMc < nPrim; ++iMc)
276   {
277     TParticle* particle = stack->Particle(iMc);
278     
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;
287     
288     if (SignOK(particle->GetPDG()) == kFALSE)
289       continue;
290     
291     Float_t eta = particle->Eta();
292     Float_t pt = particle->Pt();
293     
294     Int_t id = -1;
295     switch (TMath::Abs(particle->GetPdgCode()))
296       {
297       case 211: id = 0; break;
298       case 321: id = 1; break;
299       case 2212: id = 2; break;
300       case 11:
301         {
302           /*if (pt < 0.1 && particle->GetMother(0) > -1)
303           {
304             TParticle* mother = stack->Particle(particle->GetMother(0));
305             printf("Mother pdg code is %d\n", mother->GetPdgCode());
306           } */
307           //particle->Dump();
308           //if (particle->GetUniqueID() == 1)
309
310           //printf("Found illegal particle mc id = %d file = %s event = %d\n", iMc,  fTree->GetCurrentFile()->GetName(), fTree->GetTree()->GetReadEntry());
311         }
312       default: id = 3; break;
313     }
314
315     if (eventTriggered & vertexReconstructed) {
316       if (fdNdEtaCorrectionSpecies[id])
317         fdNdEtaCorrectionSpecies[id]->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, vertexReconstructed, 1);
318       //if (pt < 0.1)
319       if (fPIDParticles)
320         fPIDParticles->Fill(particle->GetPdgCode());
321     }
322     
323     // non diffractive
324     if (processtype==AliPWG0Helper::kND) { 
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==AliPWG0Helper::kSD) { 
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==AliPWG0Helper::kDD) { 
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);
337     }
338
339   }// end of mc particle
340
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
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)
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));
369       continue;
370     }
371
372     Float_t eta = particle->Eta();
373     Float_t pt  = particle->Pt();
374     
375     // non diffractive
376     if (processtype==AliPWG0Helper::kND) { 
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==AliPWG0Helper::kSD) { 
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==AliPWG0Helper::kDD) { 
389       if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[2]->FillTrackedParticle(vtxMC[2], eta, pt);
390       if (vertexRecoStudy)  fdNdEtaCorrectionVertexReco[2] ->FillTrackedParticle(vtxMC[2], eta, pt);
391     }
392
393     // find primary particle that created this particle
394     TParticle* mother = AliPWG0Helper::FindPrimaryMother(stack, label);
395     if (!mother)
396       continue;
397
398     if (SignOK(mother->GetPDG()) == kFALSE)
399       continue;
400
401     //printf("We continue with particle %d (pdg %d)\n", label, mother->GetPdgCode());
402
403     Int_t id = -1;
404     switch (TMath::Abs(mother->GetPdgCode()))
405     {
406       case 211:  id = 0; break;
407       case 321:  id = 1; break;
408       case 2212: id = 2; break;
409       default:   id = 3; break;
410     }
411
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());
418     }
419   } // end of track loop
420
421   
422   delete iter;
423   iter = 0;
424
425   if (fSecondaries)
426     FillSecondaries();
427
428   if (fSigmaVertex)
429     FillSigmaVertex();
430
431
432   delete listOfTracks;
433   listOfTracks = 0;
434
435   return kTRUE;
436 }
437
438 void AlidNdEtaSystematicsSelector::FillSecondaries()
439 {
440   // fills the secondary histograms
441   
442   AliStack* stack = GetStack();
443
444   Int_t particlesPrimaries = 0;
445   Int_t particlesSecondaries = 0;
446   Int_t particlesPrimariesPtCut = 0;
447   Int_t particlesSecondariesPtCut = 0;
448
449   for (Int_t iParticle = 0; iParticle < stack->GetNtrack(); iParticle++)
450   {
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);
497     if (!esdTrack)
498     {
499       AliDebug(AliLog::kError, Form("UNEXPECTED: Could not get track %d.", iTrack));
500       continue;
501     }
502
503     Int_t label = TMath::Abs(esdTrack->GetLabel());
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
511    if (label < stack->GetNprimary())
512       allTracksPrimaries++;
513     else
514       allTracksSecondaries++;
515
516     if (!fEsdTrackCuts->AcceptTrack(esdTrack))
517       continue;
518
519     if (label < stack->GetNprimary())
520       acceptedTracksPrimaries++;
521     else
522       acceptedTracksSecondaries++;
523   }
524
525   fSecondaries->Fill(5, allTracksPrimaries);
526   fSecondaries->Fill(6, allTracksSecondaries);
527   fSecondaries->Fill(7, acceptedTracksPrimaries);
528   fSecondaries->Fill(8, acceptedTracksSecondaries);
529
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);
531 }
532
533 void 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
555     for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
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
572 void 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)
591     if (fdNdEtaCorrectionSpecies[i])
592       fOutput->Add(fdNdEtaCorrectionSpecies[i]);
593
594   if (fSigmaVertex)
595     fOutput->Add(fSigmaVertex);
596
597   for (Int_t i=0; i<3; ++i) {
598     if (fdNdEtaCorrectionVertexReco[i])
599       fOutput->Add(fdNdEtaCorrectionVertexReco[i]);
600     
601     if (fdNdEtaCorrectionTriggerBias[i])
602       fOutput->Add(fdNdEtaCorrectionTriggerBias[i]);
603   }
604   
605 }
606
607 void 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
615   fSecondaries = dynamic_cast<TH2F*> (fOutput->FindObject("fSecondaries"));
616   for (Int_t i=0; i<4; ++i)
617     fdNdEtaCorrectionSpecies[i] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
618   fSigmaVertex = dynamic_cast<TH1F*> (fOutput->FindObject("fSigmaVertex"));
619
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
628   if (fPIDParticles)
629   {
630     TDatabasePDG* pdgDB = new TDatabasePDG;
631
632     for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i)
633       if (fPIDParticles->GetBinContent(i) > 0)
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));
640
641     delete pdgDB;
642     pdgDB = 0;
643   }
644
645   for (Int_t i=0; i<3; ++i) {
646     if (fdNdEtaCorrectionVertexReco[i])
647       fdNdEtaCorrectionVertexReco[i]->Finish();
648     
649     if (fdNdEtaCorrectionTriggerBias[i])
650       fdNdEtaCorrectionTriggerBias[i]->Finish();
651   }
652
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
661   if (fSigmaVertex)
662     fSigmaVertex->Write();
663   
664   for (Int_t i=0; i<4; ++i)
665     if (fdNdEtaCorrectionSpecies[i])
666       fdNdEtaCorrectionSpecies[i]->SaveHistograms();
667
668   for (Int_t i=0; i<3; ++i) {
669     if (fdNdEtaCorrectionVertexReco[i])
670       fdNdEtaCorrectionVertexReco[i]->SaveHistograms();
671
672     if (fdNdEtaCorrectionTriggerBias[i])
673       fdNdEtaCorrectionTriggerBias[i]->SaveHistograms();
674   }
675
676   fout->Write();
677   fout->Close();
678
679   if (fSecondaries)
680   {
681     new TCanvas;
682     fSecondaries->Draw("COLZ");
683   }
684
685   if (fSigmaVertex)
686   {
687     new TCanvas;
688     fSigmaVertex->Draw();
689   }
690 }
691
692 Bool_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 }