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