effc++ changes
[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 <TH3F.h>
12 #include <TH1F.h>
13 #include <TParticle.h>
14
15 #include <AliLog.h>
16 #include <AliESD.h>
17 #include <AliGenEventHeader.h>
18 #include <AliStack.h>
19 #include <AliHeader.h>
20
21 #include "esdTrackCuts/AliESDtrackCuts.h"
22 #include "AliPWG0Helper.h"
23 #include "AlidNdEtaCorrection.h"
24
25 ClassImp(AlidNdEtaSystematicsSelector)
26
27 AlidNdEtaSystematicsSelector::AlidNdEtaSystematicsSelector() :
28   AliSelectorRL(),
29   fSecondaries(0),
30   fSigmaVertex(0),
31   fEsdTrackCuts(0),
32   fOverallPrimaries(0),
33   fOverallSecondaries(0),
34   fPIDParticles(0),
35   fPIDTracks(0)
36 {
37   //
38   // Constructor. Initialization of pointers
39   //
40
41   for (Int_t i=0; i<4; ++i)
42     fdNdEtaCorrection[i] = 0;
43 }
44
45 AlidNdEtaSystematicsSelector::~AlidNdEtaSystematicsSelector()
46 {
47   //
48   // Destructor
49   //
50 }
51
52 void AlidNdEtaSystematicsSelector::Begin(TTree* tree)
53 {
54   // Begin function
55
56   AliSelectorRL::Begin(tree);
57
58   ReadUserObjects(tree);
59 }
60
61 void AlidNdEtaSystematicsSelector::ReadUserObjects(TTree* tree)
62 {
63   // read the user objects, called from slavebegin and begin
64
65   if (!fEsdTrackCuts && fInput)
66     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fInput->FindObject("AliESDtrackCuts"));
67
68   if (!fEsdTrackCuts && tree)
69     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (tree->GetUserInfo()->FindObject("AliESDtrackCuts"));
70
71   if (!fEsdTrackCuts)
72      AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from input list.");
73 }
74
75 void AlidNdEtaSystematicsSelector::SlaveBegin(TTree* tree)
76 {
77   // The SlaveBegin() function is called after the Begin() function.
78   // When running with PROOF SlaveBegin() is called on each slave server.
79   // The tree argument is deprecated (on PROOF 0 is passed).
80
81   AliSelector::SlaveBegin(tree);
82
83   ReadUserObjects(tree);
84
85   TString option(GetOption());
86
87   printf("Running AlidNdEtaSystematicsSelector with options %s\n", option.Data());
88
89   if (option.Contains("secondaries"))
90   {
91     fSecondaries = new TH3F("fSecondaries", "fSecondaries;NacceptedTracks;CratioSecondaries;p_{T}", 2000, -0.5, 205.5, 16, 0.45, 2.05, 10, 0, 10);
92   }
93
94   if (option.Contains("particle-composition"))
95   {
96     for (Int_t i=0; i<4; ++i)
97     {
98       TString name;
99       name.Form("correction_%d", i);
100       fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
101     }
102   }
103
104   if (option.Contains("sigma-vertex"))
105   {
106     fSigmaVertex = new TH1F("fSigmaVertex", "fSigmaVertex;Nsigma2vertex;NacceptedTracks", 10, 0.25, 5.25);
107     printf("WARNING: sigma-vertex analysis enabled. This will produce weird results in the AliESDtrackCuts histograms\n");
108   }
109
110   fPIDParticles = new TH1F("pid_particles", "PID of generated primary particles", 10001, -5000.5, 5000.5);
111   fPIDTracks = new TH1F("pid_tracks", "MC PID of reconstructed tracks", 10001, -5000.5, 5000.5);
112 }
113
114 Bool_t AlidNdEtaSystematicsSelector::Process(Long64_t entry)
115 {
116   // The Process() function is called for each entry in the tree (or possibly
117   // keyed object in the case of PROOF) to be processed. The entry argument
118   // specifies which entry in the currently loaded tree is to be processed.
119   // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
120   // to read either all or the required parts of the data. When processing
121   // keyed objects with PROOF, the object is already loaded and is available
122   // via the fObject pointer.
123   //
124   // This function should contain the "body" of the analysis. It can contain
125   // simple or elaborate selection criteria, run algorithms on the data
126   // of the event and typically fill histograms.
127
128   // WARNING when a selector is used with a TChain, you must use
129   //  the pointer to the current TTree to call GetEntry(entry).
130   //  The entry is always the local entry number in the current tree.
131   //  Assuming that fTree is the pointer to the TChain being processed,
132   //  use fTree->GetTree()->GetEntry(entry).
133
134   if (AliSelectorRL::Process(entry) == kFALSE)
135     return kFALSE;
136
137   // Check prerequisites
138   if (!fESD)
139   {
140     AliDebug(AliLog::kError, "ESD branch not available");
141     return kFALSE;
142   }
143
144   if (!fEsdTrackCuts)
145   {
146     AliDebug(AliLog::kError, "fESDTrackCuts not available");
147     return kFALSE;
148   }
149
150   AliStack* stack = GetStack();
151   if (!stack)
152   {
153     AliDebug(AliLog::kError, "Stack not available");
154     return kFALSE;
155   }
156
157   TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
158
159   if (fdNdEtaCorrection[0])
160     FillCorrectionMaps(list);
161
162   if (fSecondaries)
163     FillSecondaries(list);
164
165   if (fSigmaVertex)
166     FillSigmaVertex();
167
168   delete list;
169   list = 0;
170
171   return kTRUE;
172 }
173
174 void AlidNdEtaSystematicsSelector::FillCorrectionMaps(TObjArray* listOfTracks)
175 {
176   // fills the correction maps for different particle species
177
178   AliStack* stack = GetStack();
179   AliHeader* header = GetHeader();
180
181   Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD);
182   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
183
184   for (Int_t i=0; i<4; ++i)
185   {
186     fdNdEtaCorrection[i]->IncreaseEventCount();
187     if (eventTriggered)
188       fdNdEtaCorrection[i]->IncreaseTriggeredEventCount();
189   }
190
191   // get the MC vertex
192   AliGenEventHeader* genHeader = header->GenEventHeader();
193
194   TArrayF vtxMC(3);
195   genHeader->PrimaryVertex(vtxMC);
196
197   // loop over mc particles
198   Int_t nPrim  = stack->GetNprimary();
199
200   for (Int_t iMc = 0; iMc < nPrim; ++iMc)
201   {
202     TParticle* particle = stack->Particle(iMc);
203
204     if (!particle)
205     {
206       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
207       continue;
208     }
209
210     if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
211       continue;
212
213     Float_t eta = particle->Eta();
214     Float_t pt = particle->Pt();
215
216     Int_t id = -1;
217     switch (TMath::Abs(particle->GetPdgCode()))
218     {
219       case 211: id = 0; break;
220       case 321: id = 1; break;
221       case 2212: id = 2; break;
222       case 11: 
223         {
224           if (pt < 0.1 && particle->GetMother(0) > -1)
225           {
226             TParticle* mother = stack->Particle(particle->GetMother(0));
227             printf("Mother pdg code is %d\n", mother->GetPdgCode());
228           }   
229           //particle->Dump();
230           //if (particle->GetUniqueID() == 1)
231
232           //printf("Found illegal particle mc id = %d file = %s event = %d\n", iMc,  fTree->GetCurrentFile()->GetName(), fTree->GetTree()->GetReadEntry()); 
233         }
234       default: id = 3; break;
235     }
236
237     if (vertexReconstructed)
238     {
239       fdNdEtaCorrection[id]->FillParticle(vtxMC[2], eta, pt);
240       if (pt < 0.1)
241         fPIDParticles->Fill(particle->GetPdgCode());
242     }
243
244     fdNdEtaCorrection[id]->FillParticleAllEvents(eta, pt);
245     if (eventTriggered)
246       fdNdEtaCorrection[id]->FillParticleWhenEventTriggered(eta, pt);
247   }// end of mc particle
248
249   // loop over esd tracks
250   TIterator* iter = listOfTracks->MakeIterator();
251   TObject* obj = 0;
252   while ((obj = iter->Next()) != 0)
253   {
254     AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj);
255     if (!esdTrack)
256       continue;
257
258     // using the properties of the mc particle
259     Int_t label = TMath::Abs(esdTrack->GetLabel());
260     TParticle* particle = stack->Particle(label);
261     if (!particle)
262     {
263       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
264       continue;
265     }
266
267     Int_t id = -1;
268     switch (TMath::Abs(particle->GetPdgCode()))
269     {
270       case 211:  id = 0; break;
271       case 321:  id = 1; break;
272       case 2212: id = 2; break;
273       default:   id = 3; break;
274     }
275
276     if (vertexReconstructed)
277     {
278       fdNdEtaCorrection[id]->FillParticleWhenMeasuredTrack(vtxMC[2], particle->Eta(), particle->Pt());
279       if (particle->Pt() < 0.1)
280         fPIDTracks->Fill(particle->GetPdgCode());
281     }
282   } // end of track loop
283
284   delete iter;
285   iter = 0;
286 }
287
288 void AlidNdEtaSystematicsSelector::FillSecondaries(TObjArray* listOfTracks)
289 {
290   // fills the secondary histograms
291
292   AliStack* stack = GetStack();
293
294   TH1* nPrimaries = new TH1F("secondaries_primaries", "secondaries_primaries", fSecondaries->GetZaxis()->GetNbins(), fSecondaries->GetZaxis()->GetXmin(), fSecondaries->GetZaxis()->GetXmax());
295   TH1* nSecondaries = new TH1F("secondaries_secondaries", "secondaries_secondaries", fSecondaries->GetZaxis()->GetNbins(), fSecondaries->GetZaxis()->GetXmin(), fSecondaries->GetZaxis()->GetXmax());
296
297   TIterator* iter = listOfTracks->MakeIterator();
298   TObject* obj = 0;
299   while ((obj = iter->Next()) != 0)
300   {
301     AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj);
302     if (!esdTrack)
303       continue;
304
305     Int_t label = TMath::Abs(esdTrack->GetLabel());
306     TParticle* particle = stack->Particle(label);
307     if (!particle)
308     {
309       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
310       continue;
311     }
312
313     Int_t nPrim  = stack->GetNprimary();
314     if (label < nPrim)
315       nPrimaries->Fill(particle->Pt());
316     else
317       nSecondaries->Fill(particle->Pt());
318   }
319
320   for (Int_t i=1; i<=nPrimaries->GetNbinsX(); ++i)
321   {
322     Int_t primaries = (Int_t) nPrimaries->GetBinContent(i);
323     Int_t secondaries = (Int_t) nSecondaries->GetBinContent(i);
324
325     if (primaries + secondaries > 0)
326     {
327       AliDebug(AliLog::kDebug, Form("The ratio between primaries and secondaries is %d:%d = %f", primaries, secondaries, ((secondaries > 0) ? (Double_t) primaries / secondaries : -1)));
328
329       for (Double_t factor = 0.5; factor < 2.01; factor += 0.1)
330       {
331         Double_t nTracks = (Double_t) primaries + (Double_t) secondaries * factor;
332         fSecondaries->Fill(nTracks, factor, nPrimaries->GetBinCenter(i));
333         //if (secondaries > 0) printf("We fill: %f %f %f\n", nTracks, factor, nPrimaries->GetBinCenter(i));
334       }
335     }
336   }
337
338   fOverallPrimaries += (Int_t) nPrimaries->Integral();
339   fOverallSecondaries += (Int_t) nSecondaries->Integral();
340
341   delete nPrimaries;
342   nPrimaries = 0;
343
344   delete nSecondaries;
345   nSecondaries = 0;
346
347   delete iter;
348   iter = 0;
349 }
350
351 void AlidNdEtaSystematicsSelector::FillSigmaVertex()
352 {
353   // fills the fSigmaVertex histogram
354
355   // save the old value
356   Float_t oldSigmaVertex = fEsdTrackCuts->GetMinNsigmaToVertex();
357
358   // set to maximum
359   fEsdTrackCuts->SetMinNsigmaToVertex(5);
360
361   TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
362
363   TIterator* iter = list->MakeIterator();
364   TObject* obj = 0;
365   while ((obj = iter->Next()) != 0)
366   {
367     AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj);
368     if (!esdTrack)
369       continue;
370
371     Float_t sigma2Vertex = fEsdTrackCuts->GetSigmaToVertex(esdTrack);
372
373     for (Double_t nSigma = 0.5; nSigma < 5.1; nSigma += 0.5)
374     {
375       if (sigma2Vertex < nSigma)
376         fSigmaVertex->Fill(nSigma);
377     }
378   }
379
380   delete iter;
381   iter = 0;
382
383   delete list;
384   list = 0;
385
386   // set back the old value
387   fEsdTrackCuts->SetMinNsigmaToVertex(oldSigmaVertex);
388 }
389
390 void AlidNdEtaSystematicsSelector::SlaveTerminate()
391 {
392   // The SlaveTerminate() function is called after all entries or objects
393   // have been processed. When running with PROOF SlaveTerminate() is called
394   // on each slave server.
395
396   AliSelector::SlaveTerminate();
397
398   // Add the histograms to the output on each slave server
399   if (!fOutput)
400   {
401     AliDebug(AliLog::kError, Form("ERROR: Output list not initialized."));
402     return;
403   }
404
405   if (fSecondaries)
406     fOutput->Add(fSecondaries);
407
408   for (Int_t i=0; i<4; ++i)
409     if (fdNdEtaCorrection[i])
410       fOutput->Add(fdNdEtaCorrection[i]);
411
412   if (fSigmaVertex)
413     fOutput->Add(fSigmaVertex);
414 }
415
416 void AlidNdEtaSystematicsSelector::Terminate()
417 {
418   // The Terminate() function is the last function to be called during
419   // a query. It always runs on the client, it can be used to present
420   // the results graphically or save the results to file.
421
422   AliSelector::Terminate();
423
424   fSecondaries = dynamic_cast<TH3F*> (fOutput->FindObject("fSecondaries"));
425   for (Int_t i=0; i<4; ++i)
426     fdNdEtaCorrection[i] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
427   fSigmaVertex = dynamic_cast<TH1F*> (fOutput->FindObject("fSigmaVertex"));
428
429   TDatabasePDG* pdgDB = new TDatabasePDG;
430
431   for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i)
432     if (fPIDParticles->GetBinContent(i) > 0)
433       printf("PDG = %d (%s): generated: %d, reconstructed: %d, ratio: %f\n", (Int_t) fPIDParticles->GetBinCenter(i), pdgDB->GetParticle((Int_t) fPIDParticles->GetBinCenter(i))->GetName(), (Int_t) fPIDParticles->GetBinContent(i), (Int_t) fPIDTracks->GetBinContent(i), ((fPIDTracks->GetBinContent(i) > 0) ? fPIDParticles->GetBinContent(i) / fPIDTracks->GetBinContent(i) : -1));
434
435   delete pdgDB;
436   pdgDB = 0;
437
438   TFile* fout = TFile::Open("systematics.root", "RECREATE");
439
440   if (fEsdTrackCuts)
441     fEsdTrackCuts->SaveHistograms("esd_track_cuts");
442
443   if (fSecondaries)
444   {
445     fSecondaries->Write();
446     printf("We had %d primaries and %d secondaries.\n", (Int_t) fOverallPrimaries, (Int_t) fOverallSecondaries);
447   }
448
449   if (fSigmaVertex)
450     fSigmaVertex->Write();
451
452   for (Int_t i=0; i<4; ++i)
453     if (fdNdEtaCorrection[i])
454       fdNdEtaCorrection[i]->SaveHistograms();
455
456   fout->Write();
457   fout->Close();
458
459   if (fSecondaries)
460   {
461     new TCanvas;
462     fSecondaries->Draw();
463   }
464
465   if (fSigmaVertex)
466   {
467     new TCanvas;
468     fSigmaVertex->Draw();
469   }
470 }