Warning removed.
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AlidNdEtaAnalysisMCSelector.cxx
CommitLineData
dc740de4 1/* $Id$ */
2
4dd2ad81 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>
5af55649 11#include <TH1F.h>
dc740de4 12#include <TH3F.h>
7029240a 13#include <TTree.h>
16e24ca3 14#include <TFile.h>
4dd2ad81 15
16#include <AliLog.h>
17#include <AliGenEventHeader.h>
7029240a 18#include <AliHeader.h>
0bd1f8a0 19#include <AliStack.h>
9e952c39 20#include <AliESDtrack.h>
4dd2ad81 21
16e24ca3 22#include "dNdEta/dNdEtaAnalysis.h"
74fd10b3 23#include <AliPWG0Helper.h>
24#include <AliCorrection.h>
25#include <AliCorrectionMatrix2D.h>
9e952c39 26#include "esdTrackCuts/AliESDtrackCuts.h"
4dd2ad81 27
dc740de4 28
4dd2ad81 29ClassImp(AlidNdEtaAnalysisMCSelector)
30
dc740de4 31AlidNdEtaAnalysisMCSelector::AlidNdEtaAnalysisMCSelector() :
16e24ca3 32 AliSelectorRL(),
9e952c39 33 fEsdTrackCuts(0),
16e24ca3 34 fdNdEtaAnalysis(0),
c69449cc 35 fdNdEtaAnalysisTr(0),
36 fdNdEtaAnalysisTrVtx(0),
9e952c39 37 fdNdEtaAnalysisTracks(0),
5af55649 38 fVertex(0),
74fd10b3 39 fPartPt(0),
40 fEvents(0)
4dd2ad81 41{
42 //
43 // Constructor. Initialization of pointers
44 //
45}
46
47AlidNdEtaAnalysisMCSelector::~AlidNdEtaAnalysisMCSelector()
48{
49 //
50 // Destructor
51 //
52}
53
9e952c39 54void AlidNdEtaAnalysisMCSelector::Begin(TTree* tree)
55{
56 // Begin function
57
58 ReadUserObjects(tree);
59}
60
61void 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
16e24ca3 75void 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
9e952c39 83 ReadUserObjects(tree);
84
16e24ca3 85 fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
c69449cc 86 fdNdEtaAnalysisTr = new dNdEtaAnalysis("dndetaTr", "dndetaTr");
87 fdNdEtaAnalysisTrVtx = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
9e952c39 88 fdNdEtaAnalysisTracks = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
0bd1f8a0 89 fVertex = new TH3F("vertex_check", "vertex_check", 50, -50, 50, 50, -50, 50, 50, -50, 50);
c69449cc 90 for (Int_t i=0; i<3; ++i)
91 {
74fd10b3 92 fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 60, -6, 6);
c69449cc 93 fPartEta[i]->Sumw2();
94 }
4c351225 95 fPartPt = new TH1F("dndeta_check_pt", "dndeta_check_pt", 1000, 0, 10);
c69449cc 96 fPartPt->Sumw2();
74fd10b3 97 fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 40, -20, 20);
16e24ca3 98}
99
944f0536 100void AlidNdEtaAnalysisMCSelector::Init(TTree *tree)
101{
16e24ca3 102 AliSelectorRL::Init(tree);
944f0536 103
9e952c39 104 if (!tree)
105 {
106 AliDebug(AliLog::kError, "tree not available");
107 return;
108 }
109
4c351225 110 tree->SetBranchStatus("*", 0);
c69449cc 111 tree->SetBranchStatus("fTriggerMask", 1);
112 tree->SetBranchStatus("fSPDVertex*", 1);
9e952c39 113 tree->SetBranchStatus("fTracks.fLabel", 1);
114
115 AliESDtrackCuts::EnableNeededBranches(tree);
944f0536 116}
117
4dd2ad81 118Bool_t AlidNdEtaAnalysisMCSelector::Process(Long64_t entry)
119{
16e24ca3 120 // fill the dNdEtaAnalysis class from the monte carlo
4dd2ad81 121
16e24ca3 122 if (AliSelectorRL::Process(entry) == kFALSE)
4dd2ad81 123 return kFALSE;
124
c69449cc 125 // Check prerequisites
126 if (!fESD)
127 {
128 AliDebug(AliLog::kError, "ESD branch not available");
129 return kFALSE;
130 }
131
0bd1f8a0 132 AliStack* stack = GetStack();
133 if (!stack)
dc740de4 134 {
0bd1f8a0 135 AliDebug(AliLog::kError, "Stack not available");
4dd2ad81 136 return kFALSE;
dc740de4 137 }
138
139 AliHeader* header = GetHeader();
140 if (!header)
141 {
142 AliDebug(AliLog::kError, "Header not available");
143 return kFALSE;
144 }
4dd2ad81 145
c69449cc 146 Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD);
147 Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
148
4dd2ad81 149 // get the MC vertex
dc740de4 150 AliGenEventHeader* genHeader = header->GenEventHeader();
4dd2ad81 151
152 TArrayF vtxMC(3);
153 genHeader->PrimaryVertex(vtxMC);
154
0bd1f8a0 155 // loop over mc particles
156 Int_t nPrim = stack->GetNprimary();
4dd2ad81 157
74fd10b3 158 Int_t nAcceptedParticles = 0;
159
0bd1f8a0 160 for (Int_t iMc = 0; iMc < nPrim; ++iMc)
4dd2ad81 161 {
0bd1f8a0 162 TParticle* particle = stack->Particle(iMc);
4dd2ad81 163
164 if (!particle)
165 continue;
166
45e97e28 167 if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
4dd2ad81 168 continue;
169
0bd1f8a0 170 AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
74fd10b3 171 ++nAcceptedParticles;
dc740de4 172
74fd10b3 173 fdNdEtaAnalysis->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
dc740de4 174 fVertex->Fill(particle->Vx(), particle->Vy(), particle->Vz());
5af55649 175
c69449cc 176 if (eventTriggered)
177 {
74fd10b3 178 fdNdEtaAnalysisTr->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
c69449cc 179 if (vertexReconstructed)
74fd10b3 180 fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
c69449cc 181 }
182
74fd10b3 183 if (TMath::Abs(vtxMC[2]) < 20)
c69449cc 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 }
4c351225 192
193 if (TMath::Abs(particle->Eta()) < 0.8)
194 fPartPt->Fill(particle->Pt());
4dd2ad81 195 }
74fd10b3 196
197 fEvents->Fill(vtxMC[2]);
198
199 fdNdEtaAnalysis->FillEvent(vtxMC[2], nAcceptedParticles);
c69449cc 200 if (eventTriggered)
201 {
74fd10b3 202 fdNdEtaAnalysisTr->FillEvent(vtxMC[2], nAcceptedParticles);
c69449cc 203 if (vertexReconstructed)
74fd10b3 204 fdNdEtaAnalysisTrVtx->FillEvent(vtxMC[2], nAcceptedParticles);
c69449cc 205 }
5af55649 206
9e952c39 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
4dd2ad81 244 return kTRUE;
245}
dc740de4 246
16e24ca3 247void 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);
c69449cc 263 fOutput->Add(fdNdEtaAnalysisTr);
264 fOutput->Add(fdNdEtaAnalysisTrVtx);
9e952c39 265 fOutput->Add(fdNdEtaAnalysisTracks);
266
4c351225 267 fOutput->Add(fPartPt);
74fd10b3 268 fOutput->Add(fEvents);
c69449cc 269 for (Int_t i=0; i<3; ++i)
270 fOutput->Add(fPartEta[i]);
16e24ca3 271}
272
dc740de4 273void AlidNdEtaAnalysisMCSelector::Terminate()
274{
16e24ca3 275 //
276
277 AliSelectorRL::Terminate();
278
279 fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
c69449cc 280 fdNdEtaAnalysisTr = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTr"));
281 fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
9e952c39 282 fdNdEtaAnalysisTracks = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTracks"));
4c351225 283 fPartPt = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_pt"));
74fd10b3 284 fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
c69449cc 285 for (Int_t i=0; i<3; ++i)
286 fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
16e24ca3 287
74fd10b3 288 if (!fdNdEtaAnalysis || !fdNdEtaAnalysisTr || !fdNdEtaAnalysisTrVtx || !fPartPt || !fEvents)
16e24ca3 289 {
4c351225 290 AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p", (void*) fdNdEtaAnalysis, (void*) fPartPt));
16e24ca3 291 return;
292 }
293
b4b9cacc 294 fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone);
295 fdNdEtaAnalysisTr->Finish(0, -1, AlidNdEtaCorrection::kNone);
296 fdNdEtaAnalysisTrVtx->Finish(0, -1, AlidNdEtaCorrection::kNone);
9e952c39 297 fdNdEtaAnalysisTracks->Finish(0, -1, AlidNdEtaCorrection::kNone);
c69449cc 298
74fd10b3 299 Int_t events = (Int_t) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Integral();
c69449cc 300 fPartPt->Scale(1.0/events);
301 fPartPt->Scale(1.0/fPartPt->GetBinWidth(1));
302
303 if (fPartEta[0])
304 {
74fd10b3 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));
c69449cc 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 }
d09fb536 322
16e24ca3 323 TFile* fout = new TFile("analysis_mc.root","RECREATE");
324
325 fdNdEtaAnalysis->SaveHistograms();
c69449cc 326 fdNdEtaAnalysisTr->SaveHistograms();
327 fdNdEtaAnalysisTrVtx->SaveHistograms();
9e952c39 328 fdNdEtaAnalysisTracks->SaveHistograms();
c69449cc 329 fPartPt->Write();
16e24ca3 330
331 fout->Write();
332 fout->Close();
dc740de4 333
4c351225 334 if (fPartPt)
335 {
4c351225 336 new TCanvas("control2", "control2", 500, 500);
c69449cc 337 fPartPt->DrawCopy();
4c351225 338 }
74fd10b3 339
340 if (fEvents)
341 {
342 new TCanvas("control3", "control3", 500, 500);
343 fEvents->DrawCopy();
344 }
dc740de4 345}