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