Apply misalignment also for the inactive detectors (as requested by R.Grosso)
[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),
30 fEsdTrackCuts(0)
31{
32 //
33 // Constructor. Initialization of pointers
34 //
35
36 for (Int_t i=0; i<4; ++i)
37 fdNdEtaCorrection[i] = 0;
38}
39
40AlidNdEtaSystematicsSelector::~AlidNdEtaSystematicsSelector()
41{
42 //
43 // Destructor
44 //
45}
46
47void AlidNdEtaSystematicsSelector::Begin(TTree* tree)
48{
49 // Begin function
50
51 AliSelectorRL::Begin(tree);
52
53 ReadUserObjects(tree);
54}
55
56void AlidNdEtaSystematicsSelector::ReadUserObjects(TTree* tree)
57{
58 // read the user objects, called from slavebegin and begin
59
60 if (!fEsdTrackCuts && fInput)
61 fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fInput->FindObject("AliESDtrackCuts"));
62
63 if (!fEsdTrackCuts && tree)
64 fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (tree->GetUserInfo()->FindObject("AliESDtrackCuts"));
65
66 if (!fEsdTrackCuts)
67 AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from input list.");
68}
69
70void AlidNdEtaSystematicsSelector::SlaveBegin(TTree* tree)
71{
72 // The SlaveBegin() function is called after the Begin() function.
73 // When running with PROOF SlaveBegin() is called on each slave server.
74 // The tree argument is deprecated (on PROOF 0 is passed).
75
76 AliSelector::SlaveBegin(tree);
77
78 ReadUserObjects(tree);
79
80 TString option(GetOption());
81
82 printf("Running AlidNdEtaSystematicsSelector with options %s\n", option.Data());
83
84 if (option.Contains("secondaries"))
85 {
86 fSecondaries = new TH3F("fSecondaries", "fSecondaries;NacceptedTracks;CratioSecondaries;p_{T}", 201, -0.5, 205.5, 16, 0.45, 2.05, 10, 0, 10);
87 }
88
89 if (option.Contains("particle-composition"))
90 {
91 for (Int_t i=0; i<4; ++i)
92 {
93 TString name;
94 name.Form("correction_%d", i);
95 fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
96 }
97 }
98}
99
100Bool_t AlidNdEtaSystematicsSelector::Process(Long64_t entry)
101{
102 // The Process() function is called for each entry in the tree (or possibly
103 // keyed object in the case of PROOF) to be processed. The entry argument
104 // specifies which entry in the currently loaded tree is to be processed.
105 // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
106 // to read either all or the required parts of the data. When processing
107 // keyed objects with PROOF, the object is already loaded and is available
108 // via the fObject pointer.
109 //
110 // This function should contain the "body" of the analysis. It can contain
111 // simple or elaborate selection criteria, run algorithms on the data
112 // of the event and typically fill histograms.
113
114 // WARNING when a selector is used with a TChain, you must use
115 // the pointer to the current TTree to call GetEntry(entry).
116 // The entry is always the local entry number in the current tree.
117 // Assuming that fTree is the pointer to the TChain being processed,
118 // use fTree->GetTree()->GetEntry(entry).
119
120 if (AliSelectorRL::Process(entry) == kFALSE)
121 return kFALSE;
122
123 // Check prerequisites
124 if (!fESD)
125 {
126 AliDebug(AliLog::kError, "ESD branch not available");
127 return kFALSE;
128 }
129
130 if (!fEsdTrackCuts)
131 {
132 AliDebug(AliLog::kError, "fESDTrackCuts not available");
133 return kFALSE;
134 }
135
136 AliStack* stack = GetStack();
137 if (!stack)
138 {
139 AliDebug(AliLog::kError, "Stack not available");
140 return kFALSE;
141 }
142
143 TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
144
145 if (fdNdEtaCorrection[0])
146 FillCorrectionMaps(list);
147
148 if (fSecondaries)
149 FillSecondaries(list);
150
151 delete list;
152 list = 0;
153
154 return kTRUE;
155}
156
157void AlidNdEtaSystematicsSelector::FillCorrectionMaps(TObjArray* listOfTracks)
158{
159 // fills the correction maps for different particle species
160
161 AliStack* stack = GetStack();
162 AliHeader* header = GetHeader();
163
164 Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD);
165 Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
166
167 for (Int_t i=0; i<4; ++i)
168 {
169 fdNdEtaCorrection[i]->IncreaseEventCount();
170 if (eventTriggered)
171 fdNdEtaCorrection[i]->IncreaseTriggeredEventCount();
172 }
173
174 // get the MC vertex
175 AliGenEventHeader* genHeader = header->GenEventHeader();
176
177 TArrayF vtxMC(3);
178 genHeader->PrimaryVertex(vtxMC);
179
180 // loop over mc particles
181 Int_t nPrim = stack->GetNprimary();
182
183 for (Int_t iMc = 0; iMc < nPrim; ++iMc)
184 {
185 TParticle* particle = stack->Particle(iMc);
186
187 if (!particle)
188 {
189 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
190 continue;
191 }
192
193 if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
194 continue;
195
196 Float_t eta = particle->Eta();
197 Float_t pt = particle->Pt();
198
199 Int_t id = -1;
200 switch (TMath::Abs(particle->GetPdgCode()))
201 {
202 case 211: id = 0; break;
203 case 321: id = 1; break;
204 case 2212: id = 2; break;
205 default: id = 3; break;
206 }
207
208 if (vertexReconstructed)
209 fdNdEtaCorrection[id]->FillParticle(vtxMC[2], eta, pt);
210
211 fdNdEtaCorrection[id]->FillParticleAllEvents(eta, pt);
212 if (eventTriggered)
213 fdNdEtaCorrection[id]->FillParticleWhenEventTriggered(eta, pt);
214 }// end of mc particle
215
216 // loop over esd tracks
217 TIterator* iter = listOfTracks->MakeIterator();
218 TObject* obj = 0;
219 while ((obj = iter->Next()) != 0)
220 {
221 AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj);
222 if (!esdTrack)
223 continue;
224
225 // using the properties of the mc particle
226 Int_t label = TMath::Abs(esdTrack->GetLabel());
227 if (label == 0)
228 {
229 AliDebug(AliLog::kWarning, Form("WARNING: cannot find corresponding mc part for track %d.", label));
230 continue;
231 }
232
233 TParticle* particle = stack->Particle(label);
234 if (!particle)
235 {
236 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
237 continue;
238 }
239
240 Int_t id = -1;
241 switch (TMath::Abs(particle->GetPdgCode()))
242 {
243 case 211: id = 0; break;
244 case 321: id = 1; break;
245 case 2212: id = 2; break;
246 default: id = 3; break;
247 }
248
249 if (vertexReconstructed)
250 fdNdEtaCorrection[id]->FillParticleWhenMeasuredTrack(vtxMC[2], particle->Eta(), particle->Pt());
251 } // end of track loop
252
253 delete iter;
254 iter = 0;
255}
256
257void AlidNdEtaSystematicsSelector::FillSecondaries(TObjArray* listOfTracks)
258{
259 // fills the secondary histograms
260
261 AliStack* stack = GetStack();
262
263 TH1* nPrimaries = new TH1F("secondaries_primaries", "secondaries_primaries", fSecondaries->GetZaxis()->GetNbins(), fSecondaries->GetZaxis()->GetXmin(), fSecondaries->GetZaxis()->GetXmax());
264 TH1* nSecondaries = new TH1F("secondaries_secondaries", "secondaries_secondaries", fSecondaries->GetZaxis()->GetNbins(), fSecondaries->GetZaxis()->GetXmin(), fSecondaries->GetZaxis()->GetXmax());
265
266 TIterator* iter = listOfTracks->MakeIterator();
267 TObject* obj = 0;
268 while ((obj = iter->Next()) != 0)
269 {
270 AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj);
271 if (!esdTrack)
272 continue;
273
274 Int_t label = TMath::Abs(esdTrack->GetLabel());
275 if (label == 0)
276 {
277 AliDebug(AliLog::kWarning, Form("WARNING: cannot find corresponding mc part for track %d.", label));
278 continue;
279 }
280
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)
305 fSecondaries->Fill((Double_t) primaries + (Double_t) secondaries * factor, factor, nPrimaries->GetBinCenter(i));
306 }
307 }
308
309 delete nPrimaries;
310 nPrimaries = 0;
311
312 delete nSecondaries;
313 nSecondaries = 0;
314
315 delete iter;
316 iter = 0;
317}
318
319void AlidNdEtaSystematicsSelector::SlaveTerminate()
320{
321 // The SlaveTerminate() function is called after all entries or objects
322 // have been processed. When running with PROOF SlaveTerminate() is called
323 // on each slave server.
324
325 AliSelector::SlaveTerminate();
326
327 // Add the histograms to the output on each slave server
328 if (!fOutput)
329 {
330 AliDebug(AliLog::kError, Form("ERROR: Output list not initialized."));
331 return;
332 }
333
334 if (fSecondaries)
335 fOutput->Add(fSecondaries);
336
337 for (Int_t i=0; i<4; ++i)
338 if (fdNdEtaCorrection[i])
339 fOutput->Add(fdNdEtaCorrection[i]);
340}
341
342void AlidNdEtaSystematicsSelector::Terminate()
343{
344 // The Terminate() function is the last function to be called during
345 // a query. It always runs on the client, it can be used to present
346 // the results graphically or save the results to file.
347
348 AliSelector::Terminate();
349
350 fSecondaries = dynamic_cast<TH3F*> (fOutput->FindObject("fSecondaries"));
351 for (Int_t i=0; i<4; ++i)
352 fdNdEtaCorrection[i] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
353
354 TFile* fout = TFile::Open("systematics.root", "RECREATE");
355
356 if (fEsdTrackCuts)
357 fEsdTrackCuts->SaveHistograms("esd_track_cuts");
358
359 if (fSecondaries)
360 fSecondaries->Write();
361
362 for (Int_t i=0; i<4; ++i)
363 if (fdNdEtaCorrection[i])
364 fdNdEtaCorrection[i]->SaveHistograms();
365
366 fout->Write();
367 fout->Close();
368
369 if (fSecondaries)
370 {
371 new TCanvas;
372 fSecondaries->Draw();
373 }
374}