Added Add method to the various correction classes and update AliSystematicSelection.
[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 <../PYTHIA6/AliGenPythiaEventHeader.h>
22 #include <../EVGEN/AliGenCocktailEventHeader.h>
23
24 #include "esdTrackCuts/AliESDtrackCuts.h"
25 #include "AliPWG0Helper.h"
26 #include "AliPWG0depHelper.h"
27 #include "dNdEta/AlidNdEtaCorrection.h"
28
29 ClassImp(AlidNdEtaSystematicsSelector)
30
31 AlidNdEtaSystematicsSelector::AlidNdEtaSystematicsSelector() :
32   AliSelectorRL(),
33   fSecondaries(0),
34   fSigmaVertex(0),
35   fEsdTrackCuts(0),
36   fPIDParticles(0),
37   fPIDTracks(0),
38   fSignMode(0),
39   fMultiplicityMode(0)
40 {
41   //
42   // Constructor. Initialization of pointers
43   //
44
45   for (Int_t i=0; i<4; ++i)
46     fdNdEtaCorrectionSpecies[i] = 0;
47
48   for (Int_t i=0; i<3; ++i) {
49     fdNdEtaCorrectionVertexReco[i] = 0;
50     fdNdEtaCorrectionTriggerBias[i] = 0;
51   }
52 }
53
54 AlidNdEtaSystematicsSelector::~AlidNdEtaSystematicsSelector()
55 {
56   //
57   // Destructor
58   //
59 }
60
61 void AlidNdEtaSystematicsSelector::Begin(TTree* tree)
62 {
63   // Begin function
64
65   AliSelectorRL::Begin(tree);
66
67   ReadUserObjects(tree);
68 }
69
70 void 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
84 void 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   {
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");
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);
117       fdNdEtaCorrectionSpecies[i] = new AlidNdEtaCorrection(name, name);
118     }
119   }
120
121   if (option.Contains("sigma-vertex"))
122   {
123     fSigmaVertex = new TH1F("fSigmaVertex", "fSigmaVertex;Nsigma2vertex;NacceptedTracks", 50, 0.05, 5.05);
124     printf("WARNING: sigma-vertex analysis enabled. This will produce weird results in the AliESDtrackCuts histograms\n");
125   }
126
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
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
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);
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   }
163 }
164
165 Bool_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   // --------------------------------------------------------------
209   // related to events:
210   //
211   // stuff for vertex reconstruction correction systematics
212   Bool_t vertexRecoStudy  = kFALSE;
213   Bool_t triggerBiasStudy = kFALSE;
214   if (fdNdEtaCorrectionVertexReco[0])  vertexRecoStudy = kTRUE;
215   if (fdNdEtaCorrectionTriggerBias[0]) triggerBiasStudy = kTRUE;
216
217   AliHeader* header = GetHeader();
218   if (!header) {
219     AliDebug(AliLog::kError, "Header not available");
220     return kFALSE;
221   }
222   
223   // getting process information NB: this only works for Pythia !!!
224   Int_t processtype = AliPWG0depHelper::GetPythiaEventProcessType(header);
225     
226   // can only read pythia headers, either directly or from cocktalil header
227   AliGenEventHeader* genHeader = (AliGenEventHeader*)(header->GenEventHeader());
228   
229   if (!genHeader) {
230     AliDebug(AliLog::kError, "Gen header not available");
231     return kFALSE;
232   }
233   
234   // get the MC vertex
235   TArrayF vtxMC(3);
236   genHeader->PrimaryVertex(vtxMC);
237   
238   TObjArray* listOfTracks    = fEsdTrackCuts->GetAcceptedTracks(fESD);
239   Int_t nGoodTracks          = listOfTracks->GetEntries();
240   
241   Bool_t eventTriggered      = AliPWG0Helper::IsEventTriggered(fESD);
242   Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD);
243   
244   // non diffractive
245   if (processtype!=92 && processtype!=93 && processtype!=94) { 
246     // NB: passing the wrong process type here (1), since the process type is defined by the index in the array (here non-diffractive)
247     if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[0]->FillEvent(vtxMC[2], nGoodTracks, eventTriggered, vertexReconstructed, 1);      
248     if (vertexRecoStudy)  fdNdEtaCorrectionVertexReco[0] ->FillEvent(vtxMC[2], nGoodTracks, eventTriggered, vertexReconstructed, 1);
249   }
250   
251   // single diffractive
252   if (processtype==92 || processtype==93) { 
253     if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[1]->FillEvent(vtxMC[2], nGoodTracks, eventTriggered, vertexReconstructed, 1);     
254     if (vertexRecoStudy)  fdNdEtaCorrectionVertexReco[1] ->FillEvent(vtxMC[2], nGoodTracks, eventTriggered, vertexReconstructed, 1);
255   }
256   
257   // double diffractive
258   if (processtype==94) { 
259     if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[2]->FillEvent(vtxMC[2], nGoodTracks, eventTriggered, vertexReconstructed, 1);     
260     if (vertexRecoStudy)  fdNdEtaCorrectionVertexReco[2] ->FillEvent(vtxMC[2], nGoodTracks, eventTriggered, vertexReconstructed, 1);
261   }
262
263   if (eventTriggered & vertexReconstructed) {
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   // --------------------------------------------------------------
271   // MC particle loop
272   //
273
274   Int_t nPrim  = stack->GetNprimary();
275
276   for (Int_t iMc = 0; iMc < nPrim; ++iMc)
277   {
278     TParticle* particle = stack->Particle(iMc);
279     
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;
288     
289     if (SignOK(particle->GetPDG()) == kFALSE)
290       continue;
291     
292     Float_t eta = particle->Eta();
293     Float_t pt = particle->Pt();
294     
295     Int_t id = -1;
296     switch (TMath::Abs(particle->GetPdgCode()))
297       {
298       case 211: id = 0; break;
299       case 321: id = 1; break;
300       case 2212: id = 2; break;
301       case 11:
302         {
303           /*if (pt < 0.1 && particle->GetMother(0) > -1)
304           {
305             TParticle* mother = stack->Particle(particle->GetMother(0));
306             printf("Mother pdg code is %d\n", mother->GetPdgCode());
307           } */
308           //particle->Dump();
309           //if (particle->GetUniqueID() == 1)
310
311           //printf("Found illegal particle mc id = %d file = %s event = %d\n", iMc,  fTree->GetCurrentFile()->GetName(), fTree->GetTree()->GetReadEntry());
312         }
313       default: id = 3; break;
314     }
315
316     if (eventTriggered & vertexReconstructed) {
317       if (fdNdEtaCorrectionSpecies[id])
318         fdNdEtaCorrectionSpecies[id]->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, vertexReconstructed, 1);
319       //if (pt < 0.1)
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);
338     }
339
340   }// end of mc particle
341
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
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)
362       continue;
363
364     Float_t eta = particle->Eta();
365     Float_t pt  = particle->Pt();
366     
367     // non diffractive
368     if (processtype!=92 && processtype!=93 && processtype!=94) { 
369       if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[0]->FillTrackedParticle(vtxMC[2], eta, pt);
370       if (vertexRecoStudy)  fdNdEtaCorrectionVertexReco[0] ->FillTrackedParticle(vtxMC[2], eta, pt);
371     }
372     
373     // single diffractive
374     if (processtype==92 || processtype==93) { 
375       if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[1]->FillTrackedParticle(vtxMC[2], eta, pt);
376       if (vertexRecoStudy)  fdNdEtaCorrectionVertexReco[1] ->FillTrackedParticle(vtxMC[2], eta, pt);
377     }
378
379     // double diffractive
380     if (processtype==94) { 
381       if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[2]->FillTrackedParticle(vtxMC[2], eta, pt);
382       if (vertexRecoStudy)  fdNdEtaCorrectionVertexReco[2] ->FillTrackedParticle(vtxMC[2], eta, pt);
383     }
384     
385
386     // using the properties of the mc particle
387     Int_t label = TMath::Abs(esdTrack->GetLabel());
388     TParticle* particle = stack->Particle(label);
389     if (!particle)
390     {
391       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
392       continue;
393     }
394
395     TParticle* mother = particle;
396     // find primary particle that created this particle
397     while (label >= nPrim)
398     {
399       //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
400
401       if (mother->GetMother(0) == -1)
402       {
403         AliDebug(AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
404         mother = 0;
405         break;
406       }
407
408       label = mother->GetMother(0);
409
410       mother = stack->Particle(label);
411       if (!mother)
412       {
413         AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (find mother loop).", label));
414         break;
415       }
416     }
417
418     if (!mother)
419       continue;
420
421     if (SignOK(mother->GetPDG()) == kFALSE)
422       continue;
423
424     //printf("We continue with particle %d (pdg %d)\n", label, mother->GetPdgCode());
425
426     Int_t id = -1;
427     switch (TMath::Abs(mother->GetPdgCode()))
428     {
429       case 211:  id = 0; break;
430       case 321:  id = 1; break;
431       case 2212: id = 2; break;
432       default:   id = 3; break;
433     }
434
435     if (vertexReconstructed && eventTriggered) {
436       if (fdNdEtaCorrectionSpecies[id])
437         fdNdEtaCorrectionSpecies[id]->FillTrackedParticle(vtxMC[2], eta, pt);
438       //if (pt < 0.1)
439       if (fPIDTracks)
440         fPIDTracks->Fill(particle->GetPdgCode());
441     }
442   } // end of track loop
443
444   
445   delete iter;
446   iter = 0;
447
448   if (fSecondaries)
449     FillSecondaries();
450
451   if (fSigmaVertex)
452     FillSigmaVertex();
453
454
455   delete listOfTracks;
456   listOfTracks = 0;
457
458   return kTRUE;
459 }
460
461 void AlidNdEtaSystematicsSelector::FillSecondaries()
462 {
463   // fills the secondary histograms
464   
465   AliStack* stack = GetStack();
466
467   Int_t particlesPrimaries = 0;
468   Int_t particlesSecondaries = 0;
469   Int_t particlesPrimariesPtCut = 0;
470   Int_t particlesSecondariesPtCut = 0;
471
472   for (Int_t iParticle = 0; iParticle < stack->GetNtrack(); iParticle++)
473   {
474     TParticle* particle = stack->Particle(iParticle);
475     if (!particle)
476     {
477       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (particle loop).", iParticle));
478       continue;
479     }
480
481     if (TMath::Abs(particle->Eta()) > 0.9)
482       continue;
483
484     Bool_t isPrimary = kFALSE;
485     if (iParticle < stack->GetNprimary())
486     {
487       if (AliPWG0Helper::IsPrimaryCharged(particle, stack->GetNprimary()) == kFALSE)
488         continue;
489
490       isPrimary = kTRUE;
491     }
492
493     if (isPrimary)
494       particlesPrimaries++;
495     else
496       particlesSecondaries++;
497
498     if (particle->Pt() < 0.3)
499       continue;
500
501     if (isPrimary)
502       particlesPrimariesPtCut++;
503     else
504       particlesSecondariesPtCut++;
505   }
506
507   fSecondaries->Fill(1, particlesPrimaries);
508   fSecondaries->Fill(2, particlesSecondaries);
509   fSecondaries->Fill(3, particlesPrimariesPtCut);
510   fSecondaries->Fill(4, particlesSecondariesPtCut);
511
512   Int_t allTracksPrimaries = 0;
513   Int_t allTracksSecondaries = 0;
514   Int_t acceptedTracksPrimaries = 0;
515   Int_t acceptedTracksSecondaries = 0;
516
517   for (Int_t iTrack = 0; iTrack < fESD->GetNumberOfTracks(); iTrack++)
518   {
519     AliESDtrack* esdTrack = fESD->GetTrack(iTrack);
520     if (!esdTrack)
521     {
522       AliDebug(AliLog::kError, Form("UNEXPECTED: Could not get track %d.", iTrack));
523       continue;
524     }
525
526     Int_t label = TMath::Abs(esdTrack->GetLabel());
527     TParticle* particle = stack->Particle(label);
528     if (!particle)
529     {
530       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
531       continue;
532     }
533
534    if (label < stack->GetNprimary())
535       allTracksPrimaries++;
536     else
537       allTracksSecondaries++;
538
539     if (!fEsdTrackCuts->AcceptTrack(esdTrack))
540       continue;
541
542     if (label < stack->GetNprimary())
543       acceptedTracksPrimaries++;
544     else
545       acceptedTracksSecondaries++;
546   }
547
548   fSecondaries->Fill(5, allTracksPrimaries);
549   fSecondaries->Fill(6, allTracksSecondaries);
550   fSecondaries->Fill(7, acceptedTracksPrimaries);
551   fSecondaries->Fill(8, acceptedTracksSecondaries);
552
553   //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);
554 }
555
556 void AlidNdEtaSystematicsSelector::FillSigmaVertex()
557 {
558   // fills the fSigmaVertex histogram
559
560   // save the old value
561   Float_t oldSigmaVertex = fEsdTrackCuts->GetMinNsigmaToVertex();
562
563   // set to maximum
564   fEsdTrackCuts->SetMinNsigmaToVertex(5);
565
566   TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
567
568   TIterator* iter = list->MakeIterator();
569   TObject* obj = 0;
570   while ((obj = iter->Next()) != 0)
571   {
572     AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj);
573     if (!esdTrack)
574       continue;
575
576     Float_t sigma2Vertex = fEsdTrackCuts->GetSigmaToVertex(esdTrack);
577
578     for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
579     {
580       if (sigma2Vertex < nSigma)
581         fSigmaVertex->Fill(nSigma);
582     }
583   }
584
585   delete iter;
586   iter = 0;
587
588   delete list;
589   list = 0;
590
591   // set back the old value
592   fEsdTrackCuts->SetMinNsigmaToVertex(oldSigmaVertex);
593 }
594
595 void AlidNdEtaSystematicsSelector::SlaveTerminate()
596 {
597   // The SlaveTerminate() function is called after all entries or objects
598   // have been processed. When running with PROOF SlaveTerminate() is called
599   // on each slave server.
600
601   AliSelector::SlaveTerminate();
602
603   // Add the histograms to the output on each slave server
604   if (!fOutput)
605   {
606     AliDebug(AliLog::kError, Form("ERROR: Output list not initialized."));
607     return;
608   }
609
610   if (fSecondaries)
611     fOutput->Add(fSecondaries);
612
613   for (Int_t i=0; i<4; ++i)
614     if (fdNdEtaCorrectionSpecies[i])
615       fOutput->Add(fdNdEtaCorrectionSpecies[i]);
616
617   if (fSigmaVertex)
618     fOutput->Add(fSigmaVertex);
619
620   for (Int_t i=0; i<3; ++i) {
621     if (fdNdEtaCorrectionVertexReco[i])
622       fOutput->Add(fdNdEtaCorrectionVertexReco[i]);
623     
624     if (fdNdEtaCorrectionTriggerBias[i])
625       fOutput->Add(fdNdEtaCorrectionTriggerBias[i]);
626   }
627   
628 }
629
630 void AlidNdEtaSystematicsSelector::Terminate()
631 {
632   // The Terminate() function is the last function to be called during
633   // a query. It always runs on the client, it can be used to present
634   // the results graphically or save the results to file.
635
636   AliSelector::Terminate();
637
638   fSecondaries = dynamic_cast<TH2F*> (fOutput->FindObject("fSecondaries"));
639   for (Int_t i=0; i<4; ++i)
640     fdNdEtaCorrectionSpecies[i] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
641   fSigmaVertex = dynamic_cast<TH1F*> (fOutput->FindObject("fSigmaVertex"));
642
643   fdNdEtaCorrectionVertexReco[0] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("vertexRecoND"));
644   fdNdEtaCorrectionVertexReco[1] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("vertexRecoSD"));
645   fdNdEtaCorrectionVertexReco[2] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("vertexRecoDD"));
646
647   fdNdEtaCorrectionTriggerBias[0] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("triggerBiasND"));
648   fdNdEtaCorrectionTriggerBias[1] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("triggerBiasSD"));
649   fdNdEtaCorrectionTriggerBias[2] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("triggerBiasDD"));
650
651   if (fPIDParticles)
652   {
653     TDatabasePDG* pdgDB = new TDatabasePDG;
654
655     for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i)
656       if (fPIDParticles->GetBinContent(i) > 0)
657         printf("PDG = %d (%s): generated: %d, reconstructed: %d, ratio: %f\n", 
658                (Int_t) fPIDParticles->GetBinCenter(i), 
659                pdgDB->GetParticle((Int_t) fPIDParticles->GetBinCenter(i))->GetName(), 
660                (Int_t) fPIDParticles->GetBinContent(i), 
661                (Int_t) fPIDTracks->GetBinContent(i), 
662                ((fPIDTracks->GetBinContent(i) > 0) ? fPIDParticles->GetBinContent(i) / fPIDTracks->GetBinContent(i) : -1));
663
664     delete pdgDB;
665     pdgDB = 0;
666   }
667
668   for (Int_t i=0; i<3; ++i) {
669     if (fdNdEtaCorrectionVertexReco[i])
670       fdNdEtaCorrectionVertexReco[i]->Finish();
671     
672     if (fdNdEtaCorrectionTriggerBias[i])
673       fdNdEtaCorrectionTriggerBias[i]->Finish();
674   }
675
676   TFile* fout = TFile::Open("systematics.root", "RECREATE");
677
678   if (fEsdTrackCuts)
679     fEsdTrackCuts->SaveHistograms("esd_track_cuts");
680
681   if (fSecondaries)
682     fSecondaries->Write();
683
684   if (fSigmaVertex)
685     fSigmaVertex->Write();
686   
687   for (Int_t i=0; i<4; ++i)
688     if (fdNdEtaCorrectionSpecies[i])
689       fdNdEtaCorrectionSpecies[i]->SaveHistograms();
690
691   for (Int_t i=0; i<3; ++i) {
692     if (fdNdEtaCorrectionVertexReco[i])
693       fdNdEtaCorrectionVertexReco[i]->SaveHistograms();
694
695     if (fdNdEtaCorrectionTriggerBias[i])
696       fdNdEtaCorrectionTriggerBias[i]->SaveHistograms();
697   }
698
699   fout->Write();
700   fout->Close();
701
702   if (fSecondaries)
703   {
704     new TCanvas;
705     fSecondaries->Draw("COLZ");
706   }
707
708   if (fSigmaVertex)
709   {
710     new TCanvas;
711     fSigmaVertex->Draw();
712   }
713 }
714
715 Bool_t AlidNdEtaSystematicsSelector::SignOK(TParticlePDG* particle)
716 {
717   // returns if a particle with this sign should be counted
718   // this is determined by the value of fSignMode, which should have the same sign
719   // as the charge
720   // if fSignMode is 0 all particles are counted
721
722   if (fSignMode == 0)
723     return kTRUE;
724
725   if (!particle)
726   {
727     printf("WARNING: not counting a particle that does not have a pdg particle\n");
728     return kFALSE;
729   }
730
731   Double_t charge = particle->Charge();
732
733   if (fSignMode > 0)
734     if (charge < 0)
735       return kFALSE;
736
737   if (fSignMode < 0)
738     if (charge > 0)
739       return kFALSE;
740
741   return kTRUE;
742 }