Warning removed.
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AlidNdEtaAnalysisMCSelector.cxx
1 /* $Id$ */
2
3 #include "AlidNdEtaAnalysisMCSelector.h"
4
5 #include <TStyle.h>
6 #include <TSystem.h>
7 #include <TCanvas.h>
8 #include <TParticle.h>
9 #include <TParticlePDG.h>
10 #include <TVector3.h>
11 #include <TH1F.h>
12 #include <TH3F.h>
13 #include <TTree.h>
14 #include <TFile.h>
15
16 #include <AliLog.h>
17 #include <AliGenEventHeader.h>
18 #include <AliHeader.h>
19 #include <AliStack.h>
20 #include <AliESDtrack.h>
21
22 #include "dNdEta/dNdEtaAnalysis.h"
23 #include <AliPWG0Helper.h>
24 #include <AliCorrection.h>
25 #include <AliCorrectionMatrix2D.h>
26 #include "esdTrackCuts/AliESDtrackCuts.h"
27
28
29 ClassImp(AlidNdEtaAnalysisMCSelector)
30
31 AlidNdEtaAnalysisMCSelector::AlidNdEtaAnalysisMCSelector() :
32   AliSelectorRL(),
33   fEsdTrackCuts(0),
34   fdNdEtaAnalysis(0),
35   fdNdEtaAnalysisTr(0),
36   fdNdEtaAnalysisTrVtx(0),
37   fdNdEtaAnalysisTracks(0),
38   fVertex(0),
39   fPartPt(0),
40   fEvents(0)
41 {
42   //
43   // Constructor. Initialization of pointers
44   //
45 }
46
47 AlidNdEtaAnalysisMCSelector::~AlidNdEtaAnalysisMCSelector()
48 {
49   //
50   // Destructor
51   //
52 }
53
54 void AlidNdEtaAnalysisMCSelector::Begin(TTree* tree)
55 {
56   // Begin function
57
58   ReadUserObjects(tree);
59 }
60
61 void AlidNdEtaAnalysisMCSelector::ReadUserObjects(TTree* tree)
62 {
63   // read the user objects, called from slavebegin and begin
64
65   if (!fEsdTrackCuts && fInput)
66     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fInput->FindObject("AliESDtrackCuts"));
67
68   if (!fEsdTrackCuts && tree)
69     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (tree->GetUserInfo()->FindObject("AliESDtrackCuts"));
70
71   if (!fEsdTrackCuts)
72      AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from input list.");
73 }
74
75 void AlidNdEtaAnalysisMCSelector::SlaveBegin(TTree * tree)
76 {
77   // The SlaveBegin() function is called after the Begin() function.
78   // When running with PROOF SlaveBegin() is called on each slave server.
79   // The tree argument is deprecated (on PROOF 0 is passed).
80
81   AliSelectorRL::SlaveBegin(tree);
82
83   ReadUserObjects(tree);
84
85   fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
86   fdNdEtaAnalysisTr = new dNdEtaAnalysis("dndetaTr", "dndetaTr");
87   fdNdEtaAnalysisTrVtx = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
88   fdNdEtaAnalysisTracks = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
89   fVertex = new TH3F("vertex_check", "vertex_check", 50, -50, 50, 50, -50, 50, 50, -50, 50);
90   for (Int_t i=0; i<3; ++i)
91   {
92     fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 60, -6, 6);
93     fPartEta[i]->Sumw2();
94   }
95   fPartPt =  new TH1F("dndeta_check_pt", "dndeta_check_pt", 1000, 0, 10);
96   fPartPt->Sumw2();
97   fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 40, -20, 20);
98 }
99
100 void AlidNdEtaAnalysisMCSelector::Init(TTree *tree)
101 {
102   AliSelectorRL::Init(tree);
103
104   if (!tree)
105   {
106     AliDebug(AliLog::kError, "tree not available");
107     return;
108   }
109
110   tree->SetBranchStatus("*", 0);
111   tree->SetBranchStatus("fTriggerMask", 1);
112   tree->SetBranchStatus("fSPDVertex*", 1);
113   tree->SetBranchStatus("fTracks.fLabel", 1);
114
115   AliESDtrackCuts::EnableNeededBranches(tree);
116 }
117
118 Bool_t AlidNdEtaAnalysisMCSelector::Process(Long64_t entry)
119 {
120   // fill the dNdEtaAnalysis class from the monte carlo
121
122   if (AliSelectorRL::Process(entry) == kFALSE)
123     return kFALSE;
124
125   // Check prerequisites
126   if (!fESD)
127   {
128     AliDebug(AliLog::kError, "ESD branch not available");
129     return kFALSE;
130   }
131
132   AliStack* stack = GetStack();
133   if (!stack)
134   {
135     AliDebug(AliLog::kError, "Stack not available");
136     return kFALSE;
137   }
138
139   AliHeader* header = GetHeader();
140   if (!header)
141   {
142     AliDebug(AliLog::kError, "Header not available");
143     return kFALSE;
144   }
145
146   Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD);
147   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
148
149   // get the MC vertex
150   AliGenEventHeader* genHeader = header->GenEventHeader();
151
152   TArrayF vtxMC(3);
153   genHeader->PrimaryVertex(vtxMC);
154
155   // loop over mc particles
156   Int_t nPrim  = stack->GetNprimary();
157
158   Int_t nAcceptedParticles = 0;
159
160   for (Int_t iMc = 0; iMc < nPrim; ++iMc)
161   {
162     TParticle* particle = stack->Particle(iMc);
163
164     if (!particle)
165       continue;
166
167     if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
168       continue;
169
170     AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
171     ++nAcceptedParticles;
172
173     fdNdEtaAnalysis->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
174     fVertex->Fill(particle->Vx(), particle->Vy(), particle->Vz());
175
176     if (eventTriggered)
177     {
178       fdNdEtaAnalysisTr->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
179       if (vertexReconstructed)
180         fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
181     }
182
183     if (TMath::Abs(vtxMC[2]) < 20)
184     {
185       fPartEta[0]->Fill(particle->Eta());
186
187       if (vtxMC[2] < 0)
188         fPartEta[1]->Fill(particle->Eta());
189       else
190         fPartEta[2]->Fill(particle->Eta());
191     }
192
193     if (TMath::Abs(particle->Eta()) < 0.8)
194       fPartPt->Fill(particle->Pt());
195   }
196
197   fEvents->Fill(vtxMC[2]);
198
199   fdNdEtaAnalysis->FillEvent(vtxMC[2], nAcceptedParticles);
200   if (eventTriggered)
201   {
202     fdNdEtaAnalysisTr->FillEvent(vtxMC[2], nAcceptedParticles);
203     if (vertexReconstructed)
204       fdNdEtaAnalysisTrVtx->FillEvent(vtxMC[2], nAcceptedParticles);
205   }
206
207   if (!eventTriggered || !vertexReconstructed)
208     return kTRUE;
209
210   // from tracks is only done for triggered and vertex reconstructed events
211
212   // get number of "good" tracks
213   TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
214   Int_t nGoodTracks = list->GetEntries();
215
216   // loop over esd tracks
217   for (Int_t t=0; t<nGoodTracks; t++)
218   {
219     AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(t));
220     if (!esdTrack)
221     {
222       AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", t));
223       continue;
224     }
225
226     Int_t label = TMath::Abs(esdTrack->GetLabel());
227
228     TParticle* particle = stack->Particle(label);
229     if (!particle)
230     {
231       AliDebug(AliLog::kError, Form("ERROR: Could not retrieve particle %d.", esdTrack->GetLabel()));
232       continue;
233     }
234
235     fdNdEtaAnalysisTracks->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
236   } // end of track loop
237
238   delete list;
239   list = 0;
240
241   // for event count per vertex
242   fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], nGoodTracks);
243
244   return kTRUE;
245 }
246
247 void AlidNdEtaAnalysisMCSelector::SlaveTerminate()
248 {
249   // The SlaveTerminate() function is called after all entries or objects
250   // have been processed. When running with PROOF SlaveTerminate() is called
251   // on each slave server.
252
253   AliSelectorRL::SlaveTerminate();
254
255   // Add the histograms to the output on each slave server
256   if (!fOutput)
257   {
258     AliDebug(AliLog::kError, Form("ERROR: Output list not initialized."));
259     return;
260   }
261
262   fOutput->Add(fdNdEtaAnalysis);
263   fOutput->Add(fdNdEtaAnalysisTr);
264   fOutput->Add(fdNdEtaAnalysisTrVtx);
265   fOutput->Add(fdNdEtaAnalysisTracks);
266
267   fOutput->Add(fPartPt);
268   fOutput->Add(fEvents);
269   for (Int_t i=0; i<3; ++i)
270     fOutput->Add(fPartEta[i]);
271 }
272
273 void AlidNdEtaAnalysisMCSelector::Terminate()
274 {
275   //
276
277   AliSelectorRL::Terminate();
278
279   fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
280   fdNdEtaAnalysisTr = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTr"));
281   fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
282   fdNdEtaAnalysisTracks = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTracks"));
283   fPartPt = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_pt"));
284   fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
285   for (Int_t i=0; i<3; ++i)
286     fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
287
288   if (!fdNdEtaAnalysis || !fdNdEtaAnalysisTr || !fdNdEtaAnalysisTrVtx || !fPartPt || !fEvents)
289   {
290     AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p", (void*) fdNdEtaAnalysis, (void*) fPartPt));
291     return;
292   }
293
294   fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone);
295   fdNdEtaAnalysisTr->Finish(0, -1, AlidNdEtaCorrection::kNone);
296   fdNdEtaAnalysisTrVtx->Finish(0, -1, AlidNdEtaCorrection::kNone);
297   fdNdEtaAnalysisTracks->Finish(0, -1, AlidNdEtaCorrection::kNone);
298
299   Int_t events = (Int_t) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Integral();
300   fPartPt->Scale(1.0/events);
301   fPartPt->Scale(1.0/fPartPt->GetBinWidth(1));
302
303   if (fPartEta[0])
304   {
305     Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-19.9), fEvents->GetXaxis()->FindBin(-0.001));
306     Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(19.9));
307
308     fPartEta[0]->Scale(1.0 / (events1 + events2));
309     fPartEta[1]->Scale(1.0 / events1);
310     fPartEta[2]->Scale(1.0 / events2);
311
312     for (Int_t i=0; i<3; ++i)
313       fPartEta[i]->Scale(1.0 / fPartEta[i]->GetBinWidth(1));
314
315     new TCanvas("control", "control", 500, 500);
316     for (Int_t i=0; i<3; ++i)
317     {
318       fPartEta[i]->SetLineColor(i+1);
319       fPartEta[i]->Draw((i==0) ? "" : "SAME");
320     }
321   }
322
323   TFile* fout = new TFile("analysis_mc.root","RECREATE");
324
325   fdNdEtaAnalysis->SaveHistograms();
326   fdNdEtaAnalysisTr->SaveHistograms();
327   fdNdEtaAnalysisTrVtx->SaveHistograms();
328   fdNdEtaAnalysisTracks->SaveHistograms();
329   fPartPt->Write();
330
331   fout->Write();
332   fout->Close();
333
334   if (fPartPt)
335   {
336     new TCanvas("control2", "control2", 500, 500);
337     fPartPt->DrawCopy();
338   }
339
340   if (fEvents)
341   {
342     new TCanvas("control3", "control3", 500, 500);
343     fEvents->DrawCopy();
344   }
345 }