]>
Commit | Line | Data |
---|---|---|
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 | ||
25 | ClassImp(AlidNdEtaSystematicsSelector) | |
26 | ||
27 | AlidNdEtaSystematicsSelector::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 | ||
43 | AlidNdEtaSystematicsSelector::~AlidNdEtaSystematicsSelector() | |
44 | { | |
45 | // | |
46 | // Destructor | |
47 | // | |
48 | } | |
49 | ||
50 | void AlidNdEtaSystematicsSelector::Begin(TTree* tree) | |
51 | { | |
52 | // Begin function | |
53 | ||
54 | AliSelectorRL::Begin(tree); | |
55 | ||
56 | ReadUserObjects(tree); | |
57 | } | |
58 | ||
59 | void 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 | ||
73 | void 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 | ||
109 | Bool_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 | ||
169 | void 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 | ||
263 | void 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 | 326 | void 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 | 365 | void 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 | ||
391 | void 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 | } |