adding selector that creates histograms needed for systematic uncertainty
[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     Float_t eta = particle->Eta();
219     Float_t pt = particle->Pt();
220
221     if (vertexReconstructed)
222     {
223       fdNdEtaCorrection->FillParticle(vtxMC[2], eta, pt);
224       if (pt > 0.1 && pt < 0.2)
225         fPIDParticles->Fill(particle->GetPdgCode());
226     }
227
228     fdNdEtaCorrection->FillParticleAllEvents(eta, pt);
229     if (eventTriggered)
230       fdNdEtaCorrection->FillParticleWhenEventTriggered(eta, pt);
231   }// end of mc particle
232
233   // ########################################################
234   // loop over esd tracks
235   Int_t nTracks = fESD->GetNumberOfTracks();
236
237   // count the number of "good" tracks as parameter for vertex reconstruction efficiency
238   Int_t nGoodTracks = 0;
239   for (Int_t t=0; t<nTracks; t++)
240   {
241     AliDebug(AliLog::kDebug+1, Form("ESD Loop: Processing track %d.", t));
242
243     AliESDtrack* esdTrack = fESD->GetTrack(t);
244
245     // cut the esd track?
246     if (!fEsdTrackCuts->AcceptTrack(esdTrack))
247       continue;
248
249     nGoodTracks++;
250
251     // using the properties of the mc particle
252     Int_t label = TMath::Abs(esdTrack->GetLabel());
253     if (label == 0)
254     {
255       AliDebug(AliLog::kWarning, Form("WARNING: cannot find corresponding mc part for track %d.", t));
256       continue;
257     }
258
259     TParticle* particle = stack->Particle(label);
260     if (!particle)
261     {
262       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
263       continue;
264     }
265
266     if (SignOK(particle->GetPDG()) == kFALSE)
267         continue;
268
269     if (vertexReconstructed)
270     {
271       fdNdEtaCorrection->FillParticleWhenMeasuredTrack(vtxMC[2], particle->Eta(), particle->Pt());
272       if (particle->Pt() > 0.1 && particle->Pt() < 0.2)
273       {
274         fPIDTracks->Fill(particle->GetPdgCode());
275         if (particle->GetPDG()->Charge() > 0)
276         {
277           fClustersITSPos->Fill(esdTrack->GetITSclusters(0));
278           fClustersTPCPos->Fill(esdTrack->GetTPCclusters(0));
279         }
280         else
281         {
282           fClustersITSNeg->Fill(esdTrack->GetITSclusters(0));
283           fClustersTPCNeg->Fill(esdTrack->GetTPCclusters(0));
284         }
285       }
286     }
287   } // end of track loop
288
289   fdNdEtaCorrection->FillEvent(vtxMC[2], nGoodTracks);
290   if (eventTriggered)
291   {
292     fdNdEtaCorrection->FillEventWithTrigger(vtxMC[2], nGoodTracks);
293     if (vertexReconstructed)
294       fdNdEtaCorrection->FillEventWithTriggerWithReconstructedVertex(vtxMC[2], nGoodTracks);
295   }
296
297   return kTRUE;
298 }
299
300 void AlidNdEtaCorrectionSelector::SlaveTerminate()
301 {
302   // The SlaveTerminate() function is called after all entries or objects
303   // have been processed. When running with PROOF SlaveTerminate() is called
304   // on each slave server.
305
306   AliSelectorRL::SlaveTerminate();
307
308   // Add the histograms to the output on each slave server
309   if (!fOutput)
310   {
311     AliDebug(AliLog::kError, "ERROR: Output list not initialized");
312     return;
313   }
314
315   fOutput->Add(fdNdEtaCorrection);
316 }
317
318 void AlidNdEtaCorrectionSelector::Terminate()
319 {
320   // The Terminate() function is the last function to be called during
321   // a query. It always runs on the client, it can be used to present
322   // the results graphically or save the results to file.
323
324   AliSelectorRL::Terminate();
325
326   fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction"));
327   if (!fdNdEtaCorrection)
328   {
329     AliDebug(AliLog::kError, "Could not read object from output list");
330     return;
331   }
332
333   fdNdEtaCorrection->Finish();
334
335   TFile* fout = new TFile(Form("correction_map%s.root", GetOption()), "RECREATE");
336
337   fEsdTrackCuts->SaveHistograms("esd_track_cuts");
338   fdNdEtaCorrection->SaveHistograms();
339
340   fout->Write();
341   fout->Close();
342
343   fdNdEtaCorrection->DrawHistograms();
344
345   new TCanvas("pidcanvas", "pidcanvas", 500, 500);
346
347   fPIDParticles->Draw();
348   fPIDTracks->SetLineColor(2);
349   fPIDTracks->Draw("SAME");
350
351   TDatabasePDG* pdgDB = new TDatabasePDG;
352
353   for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i)
354     if (fPIDParticles->GetBinContent(i) > 0)
355       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));
356
357   delete pdgDB;
358   pdgDB = 0;
359
360   TCanvas* canvas = new TCanvas("clusters", "clusters", 1000, 500);
361   canvas->Divide(2, 1);
362
363   canvas->cd(1);
364   fClustersITSPos->Draw();
365   fClustersITSNeg->SetLineColor(kRed);
366   fClustersITSNeg->Draw("SAME");
367
368   canvas->cd(2);
369   fClustersTPCPos->Draw();
370   fClustersTPCNeg->SetLineColor(kRed);
371   fClustersTPCNeg->Draw("SAME");
372 }