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