EffC++ warnings corrected.
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AlidNdEtaSystematicsSelector.cxx
CommitLineData
10ebe68d 1/* $Id$ */
2
3#include "AlidNdEtaSystematicsSelector.h"
4
5#include <TStyle.h>
6#include <TSystem.h>
7#include <TCanvas.h>
8#include <TVector3.h>
9#include <TChain.h>
10#include <TFile.h>
11#include <TH3F.h>
12#include <TH1F.h>
13#include <TParticle.h>
14
15#include <AliLog.h>
16#include <AliESD.h>
17#include <AliGenEventHeader.h>
18#include <AliStack.h>
19#include <AliHeader.h>
20
21#include "esdTrackCuts/AliESDtrackCuts.h"
22#include "AliPWG0Helper.h"
23#include "AlidNdEtaCorrection.h"
24
25ClassImp(AlidNdEtaSystematicsSelector)
26
27AlidNdEtaSystematicsSelector::AlidNdEtaSystematicsSelector() :
28 AliSelectorRL(),
29 fSecondaries(0),
9f469bf5 30 fSigmaVertex(0),
f31d5d49 31 fEsdTrackCuts(0),
32 fOverallPrimaries(0),
33 fOverallSecondaries(0)
10ebe68d 34{
35 //
36 // Constructor. Initialization of pointers
37 //
38
39 for (Int_t i=0; i<4; ++i)
40 fdNdEtaCorrection[i] = 0;
41}
42
43AlidNdEtaSystematicsSelector::~AlidNdEtaSystematicsSelector()
44{
45 //
46 // Destructor
47 //
48}
49
50void AlidNdEtaSystematicsSelector::Begin(TTree* tree)
51{
52 // Begin function
53
54 AliSelectorRL::Begin(tree);
55
56 ReadUserObjects(tree);
57}
58
59void AlidNdEtaSystematicsSelector::ReadUserObjects(TTree* tree)
60{
61 // read the user objects, called from slavebegin and begin
62
63 if (!fEsdTrackCuts && fInput)
64 fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fInput->FindObject("AliESDtrackCuts"));
65
66 if (!fEsdTrackCuts && tree)
67 fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (tree->GetUserInfo()->FindObject("AliESDtrackCuts"));
68
69 if (!fEsdTrackCuts)
70 AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from input list.");
71}
72
73void AlidNdEtaSystematicsSelector::SlaveBegin(TTree* tree)
74{
75 // The SlaveBegin() function is called after the Begin() function.
76 // When running with PROOF SlaveBegin() is called on each slave server.
77 // The tree argument is deprecated (on PROOF 0 is passed).
78
79 AliSelector::SlaveBegin(tree);
80
81 ReadUserObjects(tree);
82
83 TString option(GetOption());
84
85 printf("Running AlidNdEtaSystematicsSelector with options %s\n", option.Data());
86
87 if (option.Contains("secondaries"))
88 {
9f469bf5 89 fSecondaries = new TH3F("fSecondaries", "fSecondaries;NacceptedTracks;CratioSecondaries;p_{T}", 2000, -0.5, 205.5, 16, 0.45, 2.05, 10, 0, 10);
10ebe68d 90 }
91
92 if (option.Contains("particle-composition"))
93 {
94 for (Int_t i=0; i<4; ++i)
95 {
96 TString name;
97 name.Form("correction_%d", i);
98 fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
99 }
100 }
9f469bf5 101
102 if (option.Contains("sigma-vertex"))
103 {
104 fSigmaVertex = new TH1F("fSigmaVertex", "fSigmaVertex;Nsigma2vertex;NacceptedTracks", 10, 0.25, 5.25);
105 printf("WARNING: sigma-vertex analysis enabled. This will produce weird results in the AliESDtrackCuts histograms\n");
106 }
10ebe68d 107}
108
109Bool_t AlidNdEtaSystematicsSelector::Process(Long64_t entry)
110{
111 // The Process() function is called for each entry in the tree (or possibly
112 // keyed object in the case of PROOF) to be processed. The entry argument
113 // specifies which entry in the currently loaded tree is to be processed.
114 // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
115 // to read either all or the required parts of the data. When processing
116 // keyed objects with PROOF, the object is already loaded and is available
117 // via the fObject pointer.
118 //
119 // This function should contain the "body" of the analysis. It can contain
120 // simple or elaborate selection criteria, run algorithms on the data
121 // of the event and typically fill histograms.
122
123 // WARNING when a selector is used with a TChain, you must use
124 // the pointer to the current TTree to call GetEntry(entry).
125 // The entry is always the local entry number in the current tree.
126 // Assuming that fTree is the pointer to the TChain being processed,
127 // use fTree->GetTree()->GetEntry(entry).
128
129 if (AliSelectorRL::Process(entry) == kFALSE)
130 return kFALSE;
131
132 // Check prerequisites
133 if (!fESD)
134 {
135 AliDebug(AliLog::kError, "ESD branch not available");
136 return kFALSE;
137 }
138
139 if (!fEsdTrackCuts)
140 {
141 AliDebug(AliLog::kError, "fESDTrackCuts not available");
142 return kFALSE;
143 }
144
145 AliStack* stack = GetStack();
146 if (!stack)
147 {
148 AliDebug(AliLog::kError, "Stack not available");
149 return kFALSE;
150 }
151
152 TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
153
154 if (fdNdEtaCorrection[0])
155 FillCorrectionMaps(list);
156
157 if (fSecondaries)
158 FillSecondaries(list);
159
9f469bf5 160 if (fSigmaVertex)
161 FillSigmaVertex();
162
10ebe68d 163 delete list;
164 list = 0;
165
166 return kTRUE;
167}
168
169void AlidNdEtaSystematicsSelector::FillCorrectionMaps(TObjArray* listOfTracks)
170{
171 // fills the correction maps for different particle species
172
173 AliStack* stack = GetStack();
174 AliHeader* header = GetHeader();
175
176 Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD);
177 Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
178
179 for (Int_t i=0; i<4; ++i)
180 {
181 fdNdEtaCorrection[i]->IncreaseEventCount();
182 if (eventTriggered)
183 fdNdEtaCorrection[i]->IncreaseTriggeredEventCount();
184 }
185
186 // get the MC vertex
187 AliGenEventHeader* genHeader = header->GenEventHeader();
188
189 TArrayF vtxMC(3);
190 genHeader->PrimaryVertex(vtxMC);
191
192 // loop over mc particles
193 Int_t nPrim = stack->GetNprimary();
194
195 for (Int_t iMc = 0; iMc < nPrim; ++iMc)
196 {
197 TParticle* particle = stack->Particle(iMc);
198
199 if (!particle)
200 {
201 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
202 continue;
203 }
204
205 if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
206 continue;
207
208 Float_t eta = particle->Eta();
209 Float_t pt = particle->Pt();
210
211 Int_t id = -1;
212 switch (TMath::Abs(particle->GetPdgCode()))
213 {
214 case 211: id = 0; break;
215 case 321: id = 1; break;
216 case 2212: id = 2; break;
217 default: id = 3; break;
218 }
219
220 if (vertexReconstructed)
221 fdNdEtaCorrection[id]->FillParticle(vtxMC[2], eta, pt);
222
223 fdNdEtaCorrection[id]->FillParticleAllEvents(eta, pt);
224 if (eventTriggered)
225 fdNdEtaCorrection[id]->FillParticleWhenEventTriggered(eta, pt);
226 }// end of mc particle
227
228 // loop over esd tracks
229 TIterator* iter = listOfTracks->MakeIterator();
230 TObject* obj = 0;
231 while ((obj = iter->Next()) != 0)
232 {
233 AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj);
234 if (!esdTrack)
235 continue;
236
237 // using the properties of the mc particle
238 Int_t label = TMath::Abs(esdTrack->GetLabel());
10ebe68d 239 TParticle* particle = stack->Particle(label);
240 if (!particle)
241 {
242 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
243 continue;
244 }
245
246 Int_t id = -1;
247 switch (TMath::Abs(particle->GetPdgCode()))
248 {
9f469bf5 249 case 211: id = 0; break;
250 case 321: id = 1; break;
10ebe68d 251 case 2212: id = 2; break;
9f469bf5 252 default: id = 3; break;
10ebe68d 253 }
254
255 if (vertexReconstructed)
256 fdNdEtaCorrection[id]->FillParticleWhenMeasuredTrack(vtxMC[2], particle->Eta(), particle->Pt());
257 } // end of track loop
258
259 delete iter;
260 iter = 0;
261}
262
263void AlidNdEtaSystematicsSelector::FillSecondaries(TObjArray* listOfTracks)
264{
265 // fills the secondary histograms
266
267 AliStack* stack = GetStack();
268
269 TH1* nPrimaries = new TH1F("secondaries_primaries", "secondaries_primaries", fSecondaries->GetZaxis()->GetNbins(), fSecondaries->GetZaxis()->GetXmin(), fSecondaries->GetZaxis()->GetXmax());
270 TH1* nSecondaries = new TH1F("secondaries_secondaries", "secondaries_secondaries", fSecondaries->GetZaxis()->GetNbins(), fSecondaries->GetZaxis()->GetXmin(), fSecondaries->GetZaxis()->GetXmax());
271
272 TIterator* iter = listOfTracks->MakeIterator();
273 TObject* obj = 0;
274 while ((obj = iter->Next()) != 0)
275 {
276 AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj);
277 if (!esdTrack)
278 continue;
279
280 Int_t label = TMath::Abs(esdTrack->GetLabel());
10ebe68d 281 TParticle* particle = stack->Particle(label);
282 if (!particle)
283 {
284 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
285 continue;
286 }
287
288 Int_t nPrim = stack->GetNprimary();
289 if (label < nPrim)
290 nPrimaries->Fill(particle->Pt());
291 else
292 nSecondaries->Fill(particle->Pt());
293 }
294
295 for (Int_t i=1; i<=nPrimaries->GetNbinsX(); ++i)
296 {
297 Int_t primaries = (Int_t) nPrimaries->GetBinContent(i);
298 Int_t secondaries = (Int_t) nSecondaries->GetBinContent(i);
299
300 if (primaries + secondaries > 0)
301 {
302 AliDebug(AliLog::kDebug, Form("The ratio between primaries and secondaries is %d:%d = %f", primaries, secondaries, ((secondaries > 0) ? (Double_t) primaries / secondaries : -1)));
303
304 for (Double_t factor = 0.5; factor < 2.01; factor += 0.1)
9f469bf5 305 {
306 Double_t nTracks = (Double_t) primaries + (Double_t) secondaries * factor;
307 fSecondaries->Fill(nTracks, factor, nPrimaries->GetBinCenter(i));
308 //if (secondaries > 0) printf("We fill: %f %f %f\n", nTracks, factor, nPrimaries->GetBinCenter(i));
309 }
10ebe68d 310 }
311 }
312
f31d5d49 313 fOverallPrimaries += (Int_t) nPrimaries->Integral();
314 fOverallSecondaries += (Int_t) nSecondaries->Integral();
315
10ebe68d 316 delete nPrimaries;
317 nPrimaries = 0;
318
319 delete nSecondaries;
320 nSecondaries = 0;
321
322 delete iter;
323 iter = 0;
324}
325
9f469bf5 326void AlidNdEtaSystematicsSelector::FillSigmaVertex()
327{
328 // fills the fSigmaVertex histogram
329
330 // save the old value
331 Float_t oldSigmaVertex = fEsdTrackCuts->GetMinNsigmaToVertex();
332
333 // set to maximum
334 fEsdTrackCuts->SetMinNsigmaToVertex(5);
335
336 TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
337
338 TIterator* iter = list->MakeIterator();
339 TObject* obj = 0;
340 while ((obj = iter->Next()) != 0)
341 {
342 AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj);
343 if (!esdTrack)
344 continue;
345
346 Float_t sigma2Vertex = fEsdTrackCuts->GetSigmaToVertex(esdTrack);
347
348 for (Double_t nSigma = 0.5; nSigma < 5.1; nSigma += 0.5)
349 {
350 if (sigma2Vertex < nSigma)
351 fSigmaVertex->Fill(nSigma);
352 }
353 }
354
355 delete iter;
356 iter = 0;
357
358 delete list;
359 list = 0;
360
361 // set back the old value
362 fEsdTrackCuts->SetMinNsigmaToVertex(oldSigmaVertex);
363}
364
10ebe68d 365void AlidNdEtaSystematicsSelector::SlaveTerminate()
366{
367 // The SlaveTerminate() function is called after all entries or objects
368 // have been processed. When running with PROOF SlaveTerminate() is called
369 // on each slave server.
370
371 AliSelector::SlaveTerminate();
372
373 // Add the histograms to the output on each slave server
374 if (!fOutput)
375 {
376 AliDebug(AliLog::kError, Form("ERROR: Output list not initialized."));
377 return;
378 }
379
380 if (fSecondaries)
381 fOutput->Add(fSecondaries);
382
383 for (Int_t i=0; i<4; ++i)
384 if (fdNdEtaCorrection[i])
385 fOutput->Add(fdNdEtaCorrection[i]);
9f469bf5 386
387 if (fSigmaVertex)
388 fOutput->Add(fSigmaVertex);
10ebe68d 389}
390
391void AlidNdEtaSystematicsSelector::Terminate()
392{
393 // The Terminate() function is the last function to be called during
394 // a query. It always runs on the client, it can be used to present
395 // the results graphically or save the results to file.
396
397 AliSelector::Terminate();
398
399 fSecondaries = dynamic_cast<TH3F*> (fOutput->FindObject("fSecondaries"));
400 for (Int_t i=0; i<4; ++i)
401 fdNdEtaCorrection[i] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
9f469bf5 402 fSigmaVertex = dynamic_cast<TH1F*> (fOutput->FindObject("fSigmaVertex"));
10ebe68d 403
404 TFile* fout = TFile::Open("systematics.root", "RECREATE");
405
406 if (fEsdTrackCuts)
407 fEsdTrackCuts->SaveHistograms("esd_track_cuts");
408
409 if (fSecondaries)
f31d5d49 410 {
10ebe68d 411 fSecondaries->Write();
f31d5d49 412 printf("We had %d primaries and %d secondaries.\n", (Int_t) fOverallPrimaries, (Int_t) fOverallSecondaries);
413 }
10ebe68d 414
9f469bf5 415 if (fSigmaVertex)
416 fSigmaVertex->Write();
417
10ebe68d 418 for (Int_t i=0; i<4; ++i)
419 if (fdNdEtaCorrection[i])
420 fdNdEtaCorrection[i]->SaveHistograms();
421
422 fout->Write();
423 fout->Close();
424
425 if (fSecondaries)
426 {
427 new TCanvas;
428 fSecondaries->Draw();
429 }
9f469bf5 430
431 if (fSigmaVertex)
432 {
433 new TCanvas;
434 fSigmaVertex->Draw();
435 }
10ebe68d 436}