systematics can be studied for only positively or negatively charged particles
[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   TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
209
210   if (fMultiplicityMode == 1 && list->GetEntries() > 20 ||
211       fMultiplicityMode == 2 && list->GetEntries() < 40)
212   {
213     delete list;
214     list = 0;
215     return kTRUE;
216   }
217
218   if (fdNdEtaCorrectionSpecies[0])
219     FillCorrectionMaps(list);
220
221   if (fSecondaries)
222     FillSecondaries();
223
224   if (fSigmaVertex)
225     FillSigmaVertex();
226
227   // stuff for vertex reconstruction correction systematics
228   Bool_t vertexRecoStudy  = kFALSE;
229   Bool_t triggerBiasStudy = kFALSE;
230   if (fdNdEtaCorrectionVertexReco[0]) vertexRecoStudy = kTRUE;
231   if (fdNdEtaCorrectionTriggerBias[0]) triggerBiasStudy = kTRUE;
232
233   if (vertexRecoStudy || triggerBiasStudy) {
234     AliHeader* header = GetHeader();
235     if (!header) {
236       AliDebug(AliLog::kError, "Header not available");
237       return kFALSE;
238     }
239
240     // getting process information
241     Int_t processtype = AliPWG0depHelper::GetPythiaEventProcessType(header);
242
243     AliDebug(AliLog::kInfo, Form("Pythia process type %d.",processtype));
244
245     // can only read pythia headers, either directly or from cocktalil header
246     AliGenEventHeader* genHeader = (AliGenEventHeader*)(header->GenEventHeader());
247     
248     if (!genHeader) {
249       AliDebug(AliLog::kError, "Gen header not available");
250       return kFALSE;
251     }
252
253     // get the MC vertex
254     TArrayF vtxMC(3);
255     genHeader->PrimaryVertex(vtxMC);
256
257     Int_t nGoodTracks = list->GetEntries();
258
259     Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
260     Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD);
261
262     // non diffractive
263     if (processtype!=92 && processtype!=93 && processtype!=94) { 
264       if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[0]->FillEventAll(vtxMC[2], nGoodTracks);
265
266       if (eventTriggered) {
267         if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[0]->FillEventWithTrigger(vtxMC[2], nGoodTracks);
268         if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[0]->FillEventWithTrigger(vtxMC[2], nGoodTracks);
269
270         if (vertexReconstructed)
271           if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[0]->FillEventWithTriggerWithReconstructedVertex(vtxMC[2], nGoodTracks);
272       }
273     }
274
275     // single diffractive
276     if (processtype==92 || processtype==93) { 
277       if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[1]->FillEventAll(vtxMC[2], nGoodTracks);
278
279       if (eventTriggered) {
280         if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[1]->FillEventWithTrigger(vtxMC[2], nGoodTracks);
281         if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[1]->FillEventWithTrigger(vtxMC[2], nGoodTracks);
282
283         if (vertexReconstructed)
284           if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[1]->FillEventWithTriggerWithReconstructedVertex(vtxMC[2], nGoodTracks);
285       }
286     }
287
288     // double diffractive
289     if (processtype==94) { 
290       if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[2]->FillEventAll(vtxMC[2], nGoodTracks);
291
292       if (eventTriggered) {
293         if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[2]->FillEventWithTrigger(vtxMC[2], nGoodTracks);
294         if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[2]->FillEventWithTrigger(vtxMC[2], nGoodTracks);
295
296         if (vertexReconstructed)
297           if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[2]->FillEventWithTriggerWithReconstructedVertex(vtxMC[2], nGoodTracks);
298       }
299     }
300   }
301
302   delete list;
303   list = 0;
304
305   return kTRUE;
306 }
307
308 void AlidNdEtaSystematicsSelector::FillCorrectionMaps(TObjArray* listOfTracks)
309 {
310   // fills the correction maps for different particle species
311
312   // TODO fix the use of the FillParticle* functions
313
314   AliStack* stack = GetStack();
315   AliHeader* header = GetHeader();
316
317   Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD);
318   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
319
320   // get the MC vertex
321   AliGenEventHeader* genHeader = header->GenEventHeader();
322
323   TArrayF vtxMC(3);
324   genHeader->PrimaryVertex(vtxMC);
325
326   // loop over mc particles
327   Int_t nPrim  = stack->GetNprimary();
328
329   for (Int_t iMc = 0; iMc < nPrim; ++iMc)
330   {
331     TParticle* particle = stack->Particle(iMc);
332
333     if (!particle)
334     {
335       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
336       continue;
337     }
338
339     if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
340       continue;
341
342     if (SignOK(particle->GetPDG()) == kFALSE)
343       continue;
344
345     Float_t eta = particle->Eta();
346     Float_t pt = particle->Pt();
347
348     Int_t id = -1;
349     switch (TMath::Abs(particle->GetPdgCode()))
350     {
351       case 211: id = 0; break;
352       case 321: id = 1; break;
353       case 2212: id = 2; break;
354       case 11:
355         {
356           /*if (pt < 0.1 && particle->GetMother(0) > -1)
357           {
358             TParticle* mother = stack->Particle(particle->GetMother(0));
359             printf("Mother pdg code is %d\n", mother->GetPdgCode());
360           } */
361           //particle->Dump();
362           //if (particle->GetUniqueID() == 1)
363
364           //printf("Found illegal particle mc id = %d file = %s event = %d\n", iMc,  fTree->GetCurrentFile()->GetName(), fTree->GetTree()->GetReadEntry());
365         }
366       default: id = 3; break;
367     }
368
369     if (vertexReconstructed && eventTriggered)
370     {
371       fdNdEtaCorrectionSpecies[id]->FillParticleVertex(vtxMC[2], eta, pt);
372       //if (pt < 0.1)
373         fPIDParticles->Fill(particle->GetPdgCode());
374     }
375
376     //fdNdEtaCorrectionSpecies[id]->FillParticleAllEvents(eta, pt);
377     //if (eventTriggered)
378     // fdNdEtaCorrectionSpecies[id]->FillParticleWhenEventTriggered(eta, pt);
379   }// end of mc particle
380
381   // loop over esd tracks
382   TIterator* iter = listOfTracks->MakeIterator();
383   TObject* obj = 0;
384   while ((obj = iter->Next()) != 0)
385   {
386     AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj);
387     if (!esdTrack)
388       continue;
389
390     // using the properties of the mc particle
391     Int_t label = TMath::Abs(esdTrack->GetLabel());
392     TParticle* particle = stack->Particle(label);
393     if (!particle)
394     {
395       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
396       continue;
397     }
398
399     TParticle* mother = particle;
400     // find primary particle that created this particle
401     while (label >= nPrim)
402     {
403       //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
404
405       if (mother->GetMother(0) == -1)
406       {
407         AliDebug(AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
408         mother = 0;
409         break;
410       }
411
412       label = mother->GetMother(0);
413
414       mother = stack->Particle(label);
415       if (!mother)
416       {
417         AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (find mother loop).", label));
418         break;
419       }
420     }
421
422     if (!mother)
423       continue;
424
425     if (SignOK(mother->GetPDG()) == kFALSE)
426       continue;
427
428     //printf("We continue with particle %d (pdg %d)\n", label, mother->GetPdgCode());
429
430     Int_t id = -1;
431     switch (TMath::Abs(mother->GetPdgCode()))
432     {
433       case 211:  id = 0; break;
434       case 321:  id = 1; break;
435       case 2212: id = 2; break;
436       default:   id = 3; break;
437     }
438
439     if (vertexReconstructed && eventTriggered)
440     {
441       fdNdEtaCorrectionSpecies[id]->FillParticleTracked(vtxMC[2], particle->Eta(), particle->Pt());
442       //if (particle->Pt() < 0.1)
443         fPIDTracks->Fill(particle->GetPdgCode());
444     }
445   } // end of track loop
446
447   Int_t nGoodTracks = listOfTracks->GetEntries();
448
449   for (Int_t i=0; i<4; ++i)
450   {
451     fdNdEtaCorrectionSpecies[i]->FillEventAll(vtxMC[2], nGoodTracks);
452     if (eventTriggered)
453     {
454       fdNdEtaCorrectionSpecies[i]->FillEventWithTrigger(vtxMC[2], nGoodTracks);
455       if (vertexReconstructed)
456         fdNdEtaCorrectionSpecies[i]->FillEventWithTriggerWithReconstructedVertex(vtxMC[2], nGoodTracks);
457     }
458   }
459
460   delete iter;
461   iter = 0;
462 }
463
464 void AlidNdEtaSystematicsSelector::FillSecondaries()
465 {
466   // fills the secondary histograms
467
468   AliStack* stack = GetStack();
469
470   Int_t particlesPrimaries = 0;
471   Int_t particlesSecondaries = 0;
472   Int_t particlesPrimariesPtCut = 0;
473   Int_t particlesSecondariesPtCut = 0;
474
475   for (Int_t iParticle = 0; iParticle < stack->GetNtrack(); iParticle++)
476   {
477     TParticle* particle = stack->Particle(iParticle);
478     if (!particle)
479     {
480       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (particle loop).", iParticle));
481       continue;
482     }
483
484     if (TMath::Abs(particle->Eta()) > 0.9)
485       continue;
486
487     Bool_t isPrimary = kFALSE;
488     if (iParticle < stack->GetNprimary())
489     {
490       if (AliPWG0Helper::IsPrimaryCharged(particle, stack->GetNprimary()) == kFALSE)
491         continue;
492
493       isPrimary = kTRUE;
494     }
495
496     if (isPrimary)
497       particlesPrimaries++;
498     else
499       particlesSecondaries++;
500
501     if (particle->Pt() < 0.3)
502       continue;
503
504     if (isPrimary)
505       particlesPrimariesPtCut++;
506     else
507       particlesSecondariesPtCut++;
508   }
509
510   fSecondaries->Fill(1, particlesPrimaries);
511   fSecondaries->Fill(2, particlesSecondaries);
512   fSecondaries->Fill(3, particlesPrimariesPtCut);
513   fSecondaries->Fill(4, particlesSecondariesPtCut);
514
515   Int_t allTracksPrimaries = 0;
516   Int_t allTracksSecondaries = 0;
517   Int_t acceptedTracksPrimaries = 0;
518   Int_t acceptedTracksSecondaries = 0;
519
520   for (Int_t iTrack = 0; iTrack < fESD->GetNumberOfTracks(); iTrack++)
521   {
522     AliESDtrack* esdTrack = fESD->GetTrack(iTrack);
523     if (!esdTrack)
524     {
525       AliDebug(AliLog::kError, Form("UNEXPECTED: Could not get track %d.", iTrack));
526       continue;
527     }
528
529     Int_t label = TMath::Abs(esdTrack->GetLabel());
530     TParticle* particle = stack->Particle(label);
531     if (!particle)
532     {
533       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
534       continue;
535     }
536
537    if (label < stack->GetNprimary())
538       allTracksPrimaries++;
539     else
540       allTracksSecondaries++;
541
542     if (!fEsdTrackCuts->AcceptTrack(esdTrack))
543       continue;
544
545     if (label < stack->GetNprimary())
546       acceptedTracksPrimaries++;
547     else
548       acceptedTracksSecondaries++;
549   }
550
551   fSecondaries->Fill(5, allTracksPrimaries);
552   fSecondaries->Fill(6, allTracksSecondaries);
553   fSecondaries->Fill(7, acceptedTracksPrimaries);
554   fSecondaries->Fill(8, acceptedTracksSecondaries);
555
556   //printf("P = %d, S = %d, P_pt = %d, S_pt = %d, T_P = %d, T_S = %d, T_P_acc = %d, T_S_acc = %d\n", particlesPrimaries, particlesSecondaries, particlesPrimariesPtCut, particlesSecondariesPtCut, allTracksPrimaries, allTracksSecondaries, acceptedTracksPrimaries, acceptedTracksSecondaries);
557 }
558
559 void AlidNdEtaSystematicsSelector::FillSigmaVertex()
560 {
561   // fills the fSigmaVertex histogram
562
563   // save the old value
564   Float_t oldSigmaVertex = fEsdTrackCuts->GetMinNsigmaToVertex();
565
566   // set to maximum
567   fEsdTrackCuts->SetMinNsigmaToVertex(5);
568
569   TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
570
571   TIterator* iter = list->MakeIterator();
572   TObject* obj = 0;
573   while ((obj = iter->Next()) != 0)
574   {
575     AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj);
576     if (!esdTrack)
577       continue;
578
579     Float_t sigma2Vertex = fEsdTrackCuts->GetSigmaToVertex(esdTrack);
580
581     for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
582     {
583       if (sigma2Vertex < nSigma)
584         fSigmaVertex->Fill(nSigma);
585     }
586   }
587
588   delete iter;
589   iter = 0;
590
591   delete list;
592   list = 0;
593
594   // set back the old value
595   fEsdTrackCuts->SetMinNsigmaToVertex(oldSigmaVertex);
596 }
597
598 void AlidNdEtaSystematicsSelector::SlaveTerminate()
599 {
600   // The SlaveTerminate() function is called after all entries or objects
601   // have been processed. When running with PROOF SlaveTerminate() is called
602   // on each slave server.
603
604   AliSelector::SlaveTerminate();
605
606   // Add the histograms to the output on each slave server
607   if (!fOutput)
608   {
609     AliDebug(AliLog::kError, Form("ERROR: Output list not initialized."));
610     return;
611   }
612
613   if (fSecondaries)
614     fOutput->Add(fSecondaries);
615
616   for (Int_t i=0; i<4; ++i)
617     if (fdNdEtaCorrectionSpecies[i])
618       fOutput->Add(fdNdEtaCorrectionSpecies[i]);
619
620   if (fSigmaVertex)
621     fOutput->Add(fSigmaVertex);
622
623   for (Int_t i=0; i<3; ++i) {
624     if (fdNdEtaCorrectionVertexReco[i])
625       fOutput->Add(fdNdEtaCorrectionVertexReco[i]);
626     
627     if (fdNdEtaCorrectionTriggerBias[i])
628       fOutput->Add(fdNdEtaCorrectionTriggerBias[i]);
629   }
630   
631 }
632
633 void AlidNdEtaSystematicsSelector::Terminate()
634 {
635   // The Terminate() function is the last function to be called during
636   // a query. It always runs on the client, it can be used to present
637   // the results graphically or save the results to file.
638
639   AliSelector::Terminate();
640
641   fSecondaries = dynamic_cast<TH2F*> (fOutput->FindObject("fSecondaries"));
642   for (Int_t i=0; i<4; ++i)
643     fdNdEtaCorrectionSpecies[i] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
644   fSigmaVertex = dynamic_cast<TH1F*> (fOutput->FindObject("fSigmaVertex"));
645
646   fdNdEtaCorrectionVertexReco[0] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("vertexRecoND"));
647   fdNdEtaCorrectionVertexReco[1] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("vertexRecoSD"));
648   fdNdEtaCorrectionVertexReco[2] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("vertexRecoDD"));
649
650   fdNdEtaCorrectionTriggerBias[0] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("triggerBiasND"));
651   fdNdEtaCorrectionTriggerBias[1] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("triggerBiasSD"));
652   fdNdEtaCorrectionTriggerBias[2] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("triggerBiasDD"));
653
654   if (fPIDParticles)
655   {
656     TDatabasePDG* pdgDB = new TDatabasePDG;
657
658     for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i)
659       if (fPIDParticles->GetBinContent(i) > 0)
660         printf("PDG = %d (%s): generated: %d, reconstructed: %d, ratio: %f\n", 
661                (Int_t) fPIDParticles->GetBinCenter(i), 
662                pdgDB->GetParticle((Int_t) fPIDParticles->GetBinCenter(i))->GetName(), 
663                (Int_t) fPIDParticles->GetBinContent(i), 
664                (Int_t) fPIDTracks->GetBinContent(i), 
665                ((fPIDTracks->GetBinContent(i) > 0) ? fPIDParticles->GetBinContent(i) / fPIDTracks->GetBinContent(i) : -1));
666
667     delete pdgDB;
668     pdgDB = 0;
669   }
670
671   for (Int_t i=0; i<3; ++i) {
672     if (fdNdEtaCorrectionVertexReco[i])
673       fdNdEtaCorrectionVertexReco[i]->Finish();
674     
675     if (fdNdEtaCorrectionTriggerBias[i])
676       fdNdEtaCorrectionTriggerBias[i]->Finish();
677   }
678
679   TFile* fout = TFile::Open("systematics.root", "RECREATE");
680
681   if (fEsdTrackCuts)
682     fEsdTrackCuts->SaveHistograms("esd_track_cuts");
683
684   if (fSecondaries)
685     fSecondaries->Write();
686
687   if (fSigmaVertex)
688     fSigmaVertex->Write();
689   
690   for (Int_t i=0; i<4; ++i)
691     if (fdNdEtaCorrectionSpecies[i])
692       fdNdEtaCorrectionSpecies[i]->SaveHistograms();
693
694   for (Int_t i=0; i<3; ++i) {
695     if (fdNdEtaCorrectionVertexReco[i])
696       fdNdEtaCorrectionVertexReco[i]->SaveHistograms();
697
698     if (fdNdEtaCorrectionTriggerBias[i])
699       fdNdEtaCorrectionTriggerBias[i]->SaveHistograms();
700   }
701
702   fout->Write();
703   fout->Close();
704
705   if (fSecondaries)
706   {
707     new TCanvas;
708     fSecondaries->Draw("COLZ");
709   }
710
711   if (fSigmaVertex)
712   {
713     new TCanvas;
714     fSigmaVertex->Draw();
715   }
716 }
717
718 Bool_t AlidNdEtaSystematicsSelector::SignOK(TParticlePDG* particle)
719 {
720   // returns if a particle with this sign should be counted
721   // this is determined by the value of fSignMode, which should have the same sign
722   // as the charge
723   // if fSignMode is 0 all particles are counted
724
725   if (fSignMode == 0)
726     return kTRUE;
727
728   if (!particle)
729   {
730     printf("WARNING: not counting a particle that does not have a pdg particle\n");
731     return kFALSE;
732   }
733
734   Double_t charge = particle->Charge();
735
736   if (fSignMode > 0)
737     if (charge < 0)
738       return kFALSE;
739
740   if (fSignMode < 0)
741     if (charge > 0)
742       return kFALSE;
743
744   return kTRUE;
745 }