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