84d96ad27cf27e83f9d4928d567930a212ff5458
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AlidNdEtaCorrectionSelector.cxx
1 /* $Id$ */
2
3 #include "AlidNdEtaCorrectionSelector.h"
4
5 #include <TStyle.h>
6 #include <TSystem.h>
7 #include <TCanvas.h>
8 #include <TParticle.h>
9 #include <TParticlePDG.h>
10 #include <TH1F.h>
11
12 #include <TChain.h>
13 #include <TSelector.h>
14 #include <TFile.h>
15
16 #include <AliLog.h>
17 #include <AliGenEventHeader.h>
18 #include <AliTracker.h>
19 #include <AliHeader.h>
20 #include <AliESDVertex.h>
21 #include <AliESD.h>
22 #include <AliESDtrack.h>
23 #include <AliRunLoader.h>
24 #include <AliStack.h>
25
26 #include "esdTrackCuts/AliESDtrackCuts.h"
27 #include "dNdEta/AlidNdEtaCorrection.h"
28 #include "AliPWG0Helper.h"
29
30 ClassImp(AlidNdEtaCorrectionSelector)
31
32 AlidNdEtaCorrectionSelector::AlidNdEtaCorrectionSelector() :
33   AliSelectorRL(),
34   fEsdTrackCuts(0),
35   fdNdEtaCorrection(0),
36   fPIDParticles(0),
37   fPIDTracks(0),
38   fSignMode(0)
39 {
40   //
41   // Constructor. Initialization of pointers
42   //
43 }
44
45 AlidNdEtaCorrectionSelector::~AlidNdEtaCorrectionSelector()
46 {
47   //
48   // Destructor
49   //
50
51   // histograms are in the output list and deleted when the output
52   // list is deleted by the TSelector dtor
53 }
54
55 Bool_t AlidNdEtaCorrectionSelector::SignOK(TParticlePDG* particle)
56 {
57   // returns if a particle with this sign should be counted
58   // this is determined by the value of fSignMode, which should have the same sign
59   // as the charge
60   // if fSignMode is 0 all particles are counted
61
62   if (fSignMode == 0)
63     return kTRUE;
64
65   if (!particle)
66   {
67     printf("WARNING: not counting a particle that does not have a pdg particle\n");
68     return kFALSE;
69   }
70
71   Double_t charge = particle->Charge();
72
73   if (fSignMode > 0)
74     if (charge < 0)
75       return kFALSE;
76
77   if (fSignMode < 0)
78     if (charge > 0)
79       return kFALSE;
80
81   return kTRUE;
82 }
83
84 void AlidNdEtaCorrectionSelector::Begin(TTree * tree)
85 {
86   // The Begin() function is called at the start of the query.
87   // When running with PROOF Begin() is only called on the client.
88   // The tree argument is deprecated (on PROOF 0 is passed).
89
90   AliSelectorRL::Begin(tree);
91
92   TString option = GetOption();
93   AliInfo(Form("Called with option %s.", option.Data()));
94
95   if (option.Contains("only-positive"))
96   {
97     AliInfo("Processing only positive particles.");
98     fSignMode = 1;
99   }
100   else if (option.Contains("only-negative"))
101   {
102     AliInfo("Processing only negative particles.");
103     fSignMode = -1;
104   }
105 }
106
107 void AlidNdEtaCorrectionSelector::SlaveBegin(TTree * tree)
108 {
109   // The SlaveBegin() function is called after the Begin() function.
110   // When running with PROOF SlaveBegin() is called on each slave server.
111   // The tree argument is deprecated (on PROOF 0 is passed).
112
113   AliSelectorRL::SlaveBegin(tree);
114
115   fdNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
116
117   if (fTree)
118     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fTree->GetUserInfo()->FindObject("AliESDtrackCuts"));
119
120   if (!fEsdTrackCuts)
121     AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from user info");
122
123   fPIDParticles = new TH1F("pid_particles", "PID of generated primary particles", 10001, -5000.5, 5000.5);
124   fPIDTracks = new TH1F("pid_tracks", "MC PID of reconstructed tracks", 10001, -5000.5, 5000.5);
125
126   fClustersITSPos = new TH1F("clusters_its_pos", "clusters_its_pos", 7, -0.5, 6.5);
127   fClustersTPCPos = new TH1F("clusters_tpc_pos", "clusters_tpc_pos", 160, -0.5, 159.5);
128
129   fClustersITSNeg = new TH1F("clusters_its_neg", "clusters_its_neg", 7, -0.5, 6.5);
130   fClustersTPCNeg = new TH1F("clusters_tpc_neg", "clusters_tpc_neg", 160, -0.5, 159.5);
131 }
132
133 Bool_t AlidNdEtaCorrectionSelector::Process(Long64_t entry)
134 {
135   // The Process() function is called for each entry in the tree (or possibly
136   // keyed object in the case of PROOF) to be processed. The entry argument
137   // specifies which entry in the currently loaded tree is to be processed.
138   // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
139   // to read either all or the required parts of the data. When processing
140   // keyed objects with PROOF, the object is already loaded and is available
141   // via the fObject pointer.
142   //
143   // This function should contain the "body" of the analysis. It can contain
144   // simple or elaborate selection criteria, run algorithms on the data
145   // of the event and typically fill histograms.
146
147   // WARNING when a selector is used with a TChain, you must use
148   //  the pointer to the current TTree to call GetEntry(entry).
149   //  The entry is always the local entry number in the current tree.
150   //  Assuming that fTree is the pointer to the TChain being processed,
151   //  use fTree->GetTree()->GetEntry(entry).
152
153   if (AliSelectorRL::Process(entry) == kFALSE)
154     return kFALSE;
155
156   // check prerequesites
157   if (!fESD)
158   {
159     AliDebug(AliLog::kError, "ESD branch not available");
160     return kFALSE;
161   }
162
163   AliHeader* header = GetHeader();
164   if (!header)
165   {
166     AliDebug(AliLog::kError, "Header not available");
167     return kFALSE;
168   }
169
170   AliStack* stack = GetStack();
171   if (!stack)
172   {
173     AliDebug(AliLog::kError, "Stack not available");
174     return kFALSE;
175   }
176
177   if (!fEsdTrackCuts)
178   {
179     AliDebug(AliLog::kError, "fESDTrackCuts not available");
180     return kFALSE;
181   }
182
183   Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD);
184
185   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
186
187   fdNdEtaCorrection->IncreaseEventCount();
188   if (eventTriggered)
189     fdNdEtaCorrection->IncreaseTriggeredEventCount();
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     AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
203
204     TParticle* particle = stack->Particle(iMc);
205
206     if (!particle)
207     {
208       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
209       continue;
210     }
211
212     if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
213       continue;
214
215     if (SignOK(particle->GetPDG()) == kFALSE)
216         continue;
217
218     //if (TMath::Abs(particle->GetPdgCode()) == 2212)      continue;
219
220     Float_t eta = particle->Eta();
221     Float_t pt = particle->Pt();
222
223     if (vertexReconstructed)
224     {
225       fdNdEtaCorrection->FillParticle(vtxMC[2], eta, pt);
226       if (pt > 0.1 && pt < 0.2)
227         fPIDParticles->Fill(particle->GetPdgCode());
228     }
229
230     fdNdEtaCorrection->FillParticleAllEvents(eta, pt);
231     if (eventTriggered)
232       fdNdEtaCorrection->FillParticleWhenEventTriggered(eta, pt);
233   }// end of mc particle
234
235   // ########################################################
236   // loop over esd tracks
237   Int_t nTracks = fESD->GetNumberOfTracks();
238
239   // count the number of "good" tracks as parameter for vertex reconstruction efficiency
240   Int_t nGoodTracks = 0;
241   for (Int_t t=0; t<nTracks; t++)
242   {
243     AliDebug(AliLog::kDebug+1, Form("ESD Loop: Processing track %d.", t));
244
245     AliESDtrack* esdTrack = fESD->GetTrack(t);
246
247     // cut the esd track?
248     if (!fEsdTrackCuts->AcceptTrack(esdTrack))
249       continue;
250
251     nGoodTracks++;
252
253     // using the properties of the mc particle
254     Int_t label = TMath::Abs(esdTrack->GetLabel());
255     if (label == 0)
256     {
257       AliDebug(AliLog::kWarning, Form("WARNING: cannot find corresponding mc part for track %d.", t));
258       continue;
259     }
260
261     TParticle* particle = stack->Particle(label);
262     if (!particle)
263     {
264       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
265       continue;
266     }
267
268     if (SignOK(particle->GetPDG()) == kFALSE)
269         continue;
270
271     //if (TMath::Abs(particle->GetPdgCode()) == 2212)      continue;
272
273     if (vertexReconstructed)
274     {
275       fdNdEtaCorrection->FillParticleWhenMeasuredTrack(vtxMC[2], particle->Eta(), particle->Pt());
276       if (particle->Pt() > 0.1 && particle->Pt() < 0.2)
277       {
278         fPIDTracks->Fill(particle->GetPdgCode());
279         if (particle->GetPDG()->Charge() > 0)
280         {
281           fClustersITSPos->Fill(esdTrack->GetITSclusters(0));
282           fClustersTPCPos->Fill(esdTrack->GetTPCclusters(0));
283         }
284         else
285         {
286           fClustersITSNeg->Fill(esdTrack->GetITSclusters(0));
287           fClustersTPCNeg->Fill(esdTrack->GetTPCclusters(0));
288         }
289       }
290     }
291   } // end of track loop
292
293   fdNdEtaCorrection->FillEvent(vtxMC[2], nGoodTracks);
294   if (eventTriggered)
295   {
296     fdNdEtaCorrection->FillEventWithTrigger(vtxMC[2], nGoodTracks);
297     if (vertexReconstructed)
298       fdNdEtaCorrection->FillEventWithTriggerWithReconstructedVertex(vtxMC[2], nGoodTracks);
299   }
300
301   return kTRUE;
302 }
303
304 void AlidNdEtaCorrectionSelector::SlaveTerminate()
305 {
306   // The SlaveTerminate() function is called after all entries or objects
307   // have been processed. When running with PROOF SlaveTerminate() is called
308   // on each slave server.
309
310   AliSelectorRL::SlaveTerminate();
311
312   // Add the histograms to the output on each slave server
313   if (!fOutput)
314   {
315     AliDebug(AliLog::kError, "ERROR: Output list not initialized");
316     return;
317   }
318
319   fOutput->Add(fdNdEtaCorrection);
320 }
321
322 void AlidNdEtaCorrectionSelector::Terminate()
323 {
324   // The Terminate() function is the last function to be called during
325   // a query. It always runs on the client, it can be used to present
326   // the results graphically or save the results to file.
327
328   AliSelectorRL::Terminate();
329
330   fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction"));
331   if (!fdNdEtaCorrection)
332   {
333     AliDebug(AliLog::kError, "Could not read object from output list");
334     return;
335   }
336
337   fdNdEtaCorrection->Finish();
338
339   TFile* fout = new TFile(Form("correction_map%s.root", GetOption()), "RECREATE");
340
341   fEsdTrackCuts->SaveHistograms("esd_track_cuts");
342   fdNdEtaCorrection->SaveHistograms();
343
344   fout->Write();
345   fout->Close();
346
347   fdNdEtaCorrection->DrawHistograms();
348
349   new TCanvas("pidcanvas", "pidcanvas", 500, 500);
350
351   fPIDParticles->Draw();
352   fPIDTracks->SetLineColor(2);
353   fPIDTracks->Draw("SAME");
354
355   TDatabasePDG* pdgDB = new TDatabasePDG;
356
357   for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i)
358     if (fPIDParticles->GetBinContent(i) > 0)
359       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));
360
361   delete pdgDB;
362   pdgDB = 0;
363
364   TCanvas* canvas = new TCanvas("clusters", "clusters", 1000, 500);
365   canvas->Divide(2, 1);
366
367   canvas->cd(1);
368   fClustersITSPos->Draw();
369   fClustersITSNeg->SetLineColor(kRed);
370   fClustersITSNeg->Draw("SAME");
371
372   canvas->cd(2);
373   fClustersTPCPos->Draw();
374   fClustersTPCNeg->SetLineColor(kRed);
375   fClustersTPCNeg->Draw("SAME");
376 }