]>
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> | |
7af955da | 11 | #include <TH2F.h> |
10ebe68d | 12 | #include <TH1F.h> |
13 | #include <TParticle.h> | |
38faa1b9 | 14 | #include <TParticlePDG.h> |
10ebe68d | 15 | |
16 | #include <AliLog.h> | |
17 | #include <AliESD.h> | |
10ebe68d | 18 | #include <AliStack.h> |
19 | #include <AliHeader.h> | |
6de7db91 | 20 | #include <AliGenEventHeader.h> |
4c351225 | 21 | #include <../PYTHIA6/AliGenPythiaEventHeader.h> |
22 | #include <../EVGEN/AliGenCocktailEventHeader.h> | |
10ebe68d | 23 | |
24 | #include "esdTrackCuts/AliESDtrackCuts.h" | |
25 | #include "AliPWG0Helper.h" | |
4c351225 | 26 | #include "AliPWG0depHelper.h" |
0bd1f8a0 | 27 | #include "dNdEta/AlidNdEtaCorrection.h" |
10ebe68d | 28 | |
29 | ClassImp(AlidNdEtaSystematicsSelector) | |
30 | ||
31 | AlidNdEtaSystematicsSelector::AlidNdEtaSystematicsSelector() : | |
32 | AliSelectorRL(), | |
33 | fSecondaries(0), | |
9f469bf5 | 34 | fSigmaVertex(0), |
f31d5d49 | 35 | fEsdTrackCuts(0), |
61385583 | 36 | fPIDParticles(0), |
38faa1b9 | 37 | fPIDTracks(0), |
38 | fSignMode(0), | |
39 | fMultiplicityMode(0) | |
10ebe68d | 40 | { |
41 | // | |
42 | // Constructor. Initialization of pointers | |
43 | // | |
44 | ||
45 | for (Int_t i=0; i<4; ++i) | |
6de7db91 | 46 | fdNdEtaCorrectionSpecies[i] = 0; |
47 | ||
7b0956f3 | 48 | for (Int_t i=0; i<3; ++i) { |
6de7db91 | 49 | fdNdEtaCorrectionVertexReco[i] = 0; |
7b0956f3 | 50 | fdNdEtaCorrectionTriggerBias[i] = 0; |
51 | } | |
10ebe68d | 52 | } |
53 | ||
54 | AlidNdEtaSystematicsSelector::~AlidNdEtaSystematicsSelector() | |
55 | { | |
56 | // | |
57 | // Destructor | |
58 | // | |
59 | } | |
60 | ||
61 | void AlidNdEtaSystematicsSelector::Begin(TTree* tree) | |
62 | { | |
63 | // Begin function | |
64 | ||
65 | AliSelectorRL::Begin(tree); | |
66 | ||
67 | ReadUserObjects(tree); | |
68 | } | |
69 | ||
70 | void AlidNdEtaSystematicsSelector::ReadUserObjects(TTree* tree) | |
71 | { | |
72 | // read the user objects, called from slavebegin and begin | |
73 | ||
74 | if (!fEsdTrackCuts && fInput) | |
75 | fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fInput->FindObject("AliESDtrackCuts")); | |
76 | ||
77 | if (!fEsdTrackCuts && tree) | |
78 | fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (tree->GetUserInfo()->FindObject("AliESDtrackCuts")); | |
79 | ||
80 | if (!fEsdTrackCuts) | |
81 | AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from input list."); | |
82 | } | |
83 | ||
84 | void AlidNdEtaSystematicsSelector::SlaveBegin(TTree* tree) | |
85 | { | |
86 | // The SlaveBegin() function is called after the Begin() function. | |
87 | // When running with PROOF SlaveBegin() is called on each slave server. | |
88 | // The tree argument is deprecated (on PROOF 0 is passed). | |
89 | ||
90 | AliSelector::SlaveBegin(tree); | |
91 | ||
92 | ReadUserObjects(tree); | |
93 | ||
94 | TString option(GetOption()); | |
95 | ||
96 | printf("Running AlidNdEtaSystematicsSelector with options %s\n", option.Data()); | |
97 | ||
98 | if (option.Contains("secondaries")) | |
99 | { | |
7af955da | 100 | fSecondaries = new TH2F("fSecondaries", "fSecondaries;Case;N", 8, 0.5, 8.5, 2001, -0.5, 2000.5); |
101 | fSecondaries->GetXaxis()->SetBinLabel(1, "All Primaries"); | |
102 | fSecondaries->GetXaxis()->SetBinLabel(2, "All Secondaries"); | |
103 | fSecondaries->GetXaxis()->SetBinLabel(3, "Primaries pT > 0.3 GeV/c"); | |
104 | fSecondaries->GetXaxis()->SetBinLabel(4, "Secondaries pT > 0.3 GeV/c"); | |
105 | fSecondaries->GetXaxis()->SetBinLabel(5, "Tracks from Primaries"); | |
106 | fSecondaries->GetXaxis()->SetBinLabel(6, "Tracks from Secondaries"); | |
107 | fSecondaries->GetXaxis()->SetBinLabel(7, "Accepted Tracks from Primaries"); | |
108 | fSecondaries->GetXaxis()->SetBinLabel(8, "Accepted Tracks from Secondaries"); | |
10ebe68d | 109 | } |
110 | ||
111 | if (option.Contains("particle-composition")) | |
112 | { | |
113 | for (Int_t i=0; i<4; ++i) | |
114 | { | |
115 | TString name; | |
116 | name.Form("correction_%d", i); | |
6de7db91 | 117 | fdNdEtaCorrectionSpecies[i] = new AlidNdEtaCorrection(name, name); |
10ebe68d | 118 | } |
119 | } | |
9f469bf5 | 120 | |
121 | if (option.Contains("sigma-vertex")) | |
122 | { | |
72e597d7 | 123 | fSigmaVertex = new TH1F("fSigmaVertex", "fSigmaVertex;Nsigma2vertex;NacceptedTracks", 50, 0.05, 5.05); |
9f469bf5 | 124 | printf("WARNING: sigma-vertex analysis enabled. This will produce weird results in the AliESDtrackCuts histograms\n"); |
125 | } | |
61385583 | 126 | |
6de7db91 | 127 | if (option.Contains("vertexreco")) { |
128 | fdNdEtaCorrectionVertexReco[0] = new AlidNdEtaCorrection("vertexRecoND", "vertexRecoND"); | |
129 | fdNdEtaCorrectionVertexReco[1] = new AlidNdEtaCorrection("vertexRecoSD", "vertexRecoSD"); | |
130 | fdNdEtaCorrectionVertexReco[2] = new AlidNdEtaCorrection("vertexRecoDD", "vertexRecoDD"); | |
131 | } | |
132 | ||
7b0956f3 | 133 | if (option.Contains("triggerbias")) { |
134 | fdNdEtaCorrectionTriggerBias[0] = new AlidNdEtaCorrection("triggerBiasND", "triggerBiasND"); | |
135 | fdNdEtaCorrectionTriggerBias[1] = new AlidNdEtaCorrection("triggerBiasSD", "triggerBiasSD"); | |
136 | fdNdEtaCorrectionTriggerBias[2] = new AlidNdEtaCorrection("triggerBiasDD", "triggerBiasDD"); | |
137 | } | |
138 | ||
61385583 | 139 | fPIDParticles = new TH1F("pid_particles", "PID of generated primary particles", 10001, -5000.5, 5000.5); |
140 | fPIDTracks = new TH1F("pid_tracks", "MC PID of reconstructed tracks", 10001, -5000.5, 5000.5); | |
38faa1b9 | 141 | |
142 | if (option.Contains("only-positive")) | |
143 | { | |
144 | AliInfo("Processing only positive particles."); | |
145 | fSignMode = 1; | |
146 | } | |
147 | else if (option.Contains("only-negative")) | |
148 | { | |
149 | AliInfo("Processing only negative particles."); | |
150 | fSignMode = -1; | |
151 | } | |
152 | ||
153 | if (option.Contains("low-multiplicity")) | |
154 | { | |
155 | AliInfo("Processing only events with low multiplicity."); | |
156 | fMultiplicityMode = 1; | |
157 | } | |
158 | else if (option.Contains("high-multiplicity")) | |
159 | { | |
160 | AliInfo("Processing only events with high multiplicity."); | |
161 | fMultiplicityMode = 2; | |
162 | } | |
10ebe68d | 163 | } |
164 | ||
165 | Bool_t AlidNdEtaSystematicsSelector::Process(Long64_t entry) | |
166 | { | |
167 | // The Process() function is called for each entry in the tree (or possibly | |
168 | // keyed object in the case of PROOF) to be processed. The entry argument | |
169 | // specifies which entry in the currently loaded tree is to be processed. | |
170 | // It can be passed to either TTree::GetEntry() or TBranch::GetEntry() | |
171 | // to read either all or the required parts of the data. When processing | |
172 | // keyed objects with PROOF, the object is already loaded and is available | |
173 | // via the fObject pointer. | |
174 | // | |
175 | // This function should contain the "body" of the analysis. It can contain | |
176 | // simple or elaborate selection criteria, run algorithms on the data | |
177 | // of the event and typically fill histograms. | |
178 | ||
179 | // WARNING when a selector is used with a TChain, you must use | |
180 | // the pointer to the current TTree to call GetEntry(entry). | |
181 | // The entry is always the local entry number in the current tree. | |
182 | // Assuming that fTree is the pointer to the TChain being processed, | |
183 | // use fTree->GetTree()->GetEntry(entry). | |
184 | ||
185 | if (AliSelectorRL::Process(entry) == kFALSE) | |
186 | return kFALSE; | |
187 | ||
188 | // Check prerequisites | |
189 | if (!fESD) | |
190 | { | |
191 | AliDebug(AliLog::kError, "ESD branch not available"); | |
192 | return kFALSE; | |
193 | } | |
194 | ||
195 | if (!fEsdTrackCuts) | |
196 | { | |
197 | AliDebug(AliLog::kError, "fESDTrackCuts not available"); | |
198 | return kFALSE; | |
199 | } | |
200 | ||
201 | AliStack* stack = GetStack(); | |
202 | if (!stack) | |
203 | { | |
204 | AliDebug(AliLog::kError, "Stack not available"); | |
205 | return kFALSE; | |
206 | } | |
207 | ||
85ae901e | 208 | // -------------------------------------------------------------- |
209 | // related to events: | |
210 | // | |
6de7db91 | 211 | // stuff for vertex reconstruction correction systematics |
7b0956f3 | 212 | Bool_t vertexRecoStudy = kFALSE; |
213 | Bool_t triggerBiasStudy = kFALSE; | |
85ae901e | 214 | if (fdNdEtaCorrectionVertexReco[0]) vertexRecoStudy = kTRUE; |
7b0956f3 | 215 | if (fdNdEtaCorrectionTriggerBias[0]) triggerBiasStudy = kTRUE; |
216 | ||
85ae901e | 217 | AliHeader* header = GetHeader(); |
218 | if (!header) { | |
219 | AliDebug(AliLog::kError, "Header not available"); | |
220 | return kFALSE; | |
221 | } | |
222 | ||
223 | // getting process information NB: this only works for Pythia !!! | |
224 | Int_t processtype = AliPWG0depHelper::GetPythiaEventProcessType(header); | |
6de7db91 | 225 | |
85ae901e | 226 | // can only read pythia headers, either directly or from cocktalil header |
227 | AliGenEventHeader* genHeader = (AliGenEventHeader*)(header->GenEventHeader()); | |
228 | ||
229 | if (!genHeader) { | |
230 | AliDebug(AliLog::kError, "Gen header not available"); | |
231 | return kFALSE; | |
6de7db91 | 232 | } |
85ae901e | 233 | |
10ebe68d | 234 | // get the MC vertex |
10ebe68d | 235 | TArrayF vtxMC(3); |
236 | genHeader->PrimaryVertex(vtxMC); | |
85ae901e | 237 | |
238 | TObjArray* listOfTracks = fEsdTrackCuts->GetAcceptedTracks(fESD); | |
239 | Int_t nGoodTracks = listOfTracks->GetEntries(); | |
240 | ||
241 | Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD); | |
242 | Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD); | |
243 | ||
244 | // non diffractive | |
245 | if (processtype!=92 && processtype!=93 && processtype!=94) { | |
246 | // NB: passing the wrong process type here (1), since the process type is defined by the index in the array (here non-diffractive) | |
247 | if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[0]->FillEvent(vtxMC[2], nGoodTracks, eventTriggered, vertexReconstructed, 1); | |
248 | if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[0] ->FillEvent(vtxMC[2], nGoodTracks, eventTriggered, vertexReconstructed, 1); | |
249 | } | |
250 | ||
251 | // single diffractive | |
252 | if (processtype==92 || processtype==93) { | |
253 | if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[1]->FillEvent(vtxMC[2], nGoodTracks, eventTriggered, vertexReconstructed, 1); | |
254 | if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[1] ->FillEvent(vtxMC[2], nGoodTracks, eventTriggered, vertexReconstructed, 1); | |
255 | } | |
256 | ||
257 | // double diffractive | |
258 | if (processtype==94) { | |
259 | if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[2]->FillEvent(vtxMC[2], nGoodTracks, eventTriggered, vertexReconstructed, 1); | |
260 | if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[2] ->FillEvent(vtxMC[2], nGoodTracks, eventTriggered, vertexReconstructed, 1); | |
261 | } | |
262 | ||
263 | if (eventTriggered & vertexReconstructed) { | |
264 | for (Int_t i=0; i<4; ++i) { | |
265 | if (fdNdEtaCorrectionSpecies[i]) | |
266 | fdNdEtaCorrectionSpecies[i]->FillEvent(vtxMC[2], nGoodTracks, eventTriggered, vertexReconstructed, 1); | |
267 | } | |
268 | } | |
269 | ||
270 | // -------------------------------------------------------------- | |
271 | // MC particle loop | |
272 | // | |
10ebe68d | 273 | |
10ebe68d | 274 | Int_t nPrim = stack->GetNprimary(); |
275 | ||
276 | for (Int_t iMc = 0; iMc < nPrim; ++iMc) | |
277 | { | |
278 | TParticle* particle = stack->Particle(iMc); | |
85ae901e | 279 | |
10ebe68d | 280 | if (!particle) |
281 | { | |
282 | AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc)); | |
283 | continue; | |
284 | } | |
285 | ||
286 | if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE) | |
287 | continue; | |
85ae901e | 288 | |
38faa1b9 | 289 | if (SignOK(particle->GetPDG()) == kFALSE) |
290 | continue; | |
85ae901e | 291 | |
10ebe68d | 292 | Float_t eta = particle->Eta(); |
293 | Float_t pt = particle->Pt(); | |
85ae901e | 294 | |
10ebe68d | 295 | Int_t id = -1; |
296 | switch (TMath::Abs(particle->GetPdgCode())) | |
85ae901e | 297 | { |
10ebe68d | 298 | case 211: id = 0; break; |
299 | case 321: id = 1; break; | |
300 | case 2212: id = 2; break; | |
7af955da | 301 | case 11: |
61385583 | 302 | { |
7af955da | 303 | /*if (pt < 0.1 && particle->GetMother(0) > -1) |
61385583 | 304 | { |
305 | TParticle* mother = stack->Particle(particle->GetMother(0)); | |
306 | printf("Mother pdg code is %d\n", mother->GetPdgCode()); | |
7af955da | 307 | } */ |
61385583 | 308 | //particle->Dump(); |
309 | //if (particle->GetUniqueID() == 1) | |
310 | ||
7af955da | 311 | //printf("Found illegal particle mc id = %d file = %s event = %d\n", iMc, fTree->GetCurrentFile()->GetName(), fTree->GetTree()->GetReadEntry()); |
61385583 | 312 | } |
10ebe68d | 313 | default: id = 3; break; |
314 | } | |
315 | ||
85ae901e | 316 | if (eventTriggered & vertexReconstructed) { |
317 | if (fdNdEtaCorrectionSpecies[id]) | |
318 | fdNdEtaCorrectionSpecies[id]->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, vertexReconstructed, 1); | |
7af955da | 319 | //if (pt < 0.1) |
85ae901e | 320 | if (fPIDParticles) |
321 | fPIDParticles->Fill(particle->GetPdgCode()); | |
322 | } | |
323 | ||
324 | // non diffractive | |
325 | if (processtype!=92 && processtype!=93 && processtype!=94) { | |
326 | if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[0]->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, vertexReconstructed, 1); | |
327 | if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[0] ->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, vertexReconstructed, 1); | |
328 | } | |
329 | // single diffractive | |
330 | if (processtype==92 || processtype==93) { | |
331 | if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[1]->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, vertexReconstructed, 1); | |
332 | if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[1] ->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, vertexReconstructed, 1); | |
333 | } | |
334 | // double diffractive | |
335 | if (processtype==94) { | |
336 | if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[2]->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, vertexReconstructed, 1); | |
337 | if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[2] ->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, vertexReconstructed, 1); | |
61385583 | 338 | } |
10ebe68d | 339 | |
10ebe68d | 340 | }// end of mc particle |
341 | ||
85ae901e | 342 | |
343 | // -------------------------------------------------------------- | |
344 | // ESD track loop | |
345 | // | |
346 | ||
347 | if (fMultiplicityMode == 1 && listOfTracks->GetEntries() > 20 || | |
348 | fMultiplicityMode == 2 && listOfTracks->GetEntries() < 40) | |
349 | { | |
350 | delete listOfTracks; | |
351 | listOfTracks = 0; | |
352 | return kTRUE; | |
353 | } | |
354 | ||
10ebe68d | 355 | // loop over esd tracks |
356 | TIterator* iter = listOfTracks->MakeIterator(); | |
357 | TObject* obj = 0; | |
358 | while ((obj = iter->Next()) != 0) | |
359 | { | |
360 | AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj); | |
361 | if (!esdTrack) | |
362 | continue; | |
363 | ||
364 | // using the properties of the mc particle | |
365 | Int_t label = TMath::Abs(esdTrack->GetLabel()); | |
10ebe68d | 366 | TParticle* particle = stack->Particle(label); |
367 | if (!particle) | |
368 | { | |
369 | AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label)); | |
370 | continue; | |
371 | } | |
372 | ||
7af955da | 373 | TParticle* mother = particle; |
374 | // find primary particle that created this particle | |
375 | while (label >= nPrim) | |
376 | { | |
377 | //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0)); | |
378 | ||
379 | if (mother->GetMother(0) == -1) | |
380 | { | |
381 | AliDebug(AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label)); | |
382 | mother = 0; | |
383 | break; | |
384 | } | |
385 | ||
386 | label = mother->GetMother(0); | |
387 | ||
388 | mother = stack->Particle(label); | |
389 | if (!mother) | |
390 | { | |
391 | AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (find mother loop).", label)); | |
392 | break; | |
393 | } | |
394 | } | |
395 | ||
396 | if (!mother) | |
397 | continue; | |
398 | ||
38faa1b9 | 399 | if (SignOK(mother->GetPDG()) == kFALSE) |
400 | continue; | |
401 | ||
7af955da | 402 | //printf("We continue with particle %d (pdg %d)\n", label, mother->GetPdgCode()); |
403 | ||
10ebe68d | 404 | Int_t id = -1; |
7af955da | 405 | switch (TMath::Abs(mother->GetPdgCode())) |
10ebe68d | 406 | { |
9f469bf5 | 407 | case 211: id = 0; break; |
408 | case 321: id = 1; break; | |
10ebe68d | 409 | case 2212: id = 2; break; |
9f469bf5 | 410 | default: id = 3; break; |
10ebe68d | 411 | } |
85ae901e | 412 | Float_t eta = particle->Eta(); |
413 | Float_t pt = particle->Pt(); | |
10ebe68d | 414 | |
85ae901e | 415 | if (vertexReconstructed && eventTriggered) { |
416 | if (fdNdEtaCorrectionSpecies[id]) | |
417 | fdNdEtaCorrectionSpecies[id]->FillTrackedParticle(vtxMC[2], eta, pt); | |
418 | //if (pt < 0.1) | |
419 | if (fPIDTracks) | |
420 | fPIDTracks->Fill(particle->GetPdgCode()); | |
61385583 | 421 | } |
10ebe68d | 422 | |
85ae901e | 423 | // non diffractive |
424 | if (processtype!=92 && processtype!=93 && processtype!=94) { | |
425 | if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[0]->FillTrackedParticle(vtxMC[2], eta, pt); | |
426 | if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[0] ->FillTrackedParticle(vtxMC[2], eta, pt); | |
427 | } | |
428 | ||
429 | // single diffractive | |
430 | if (processtype==92 || processtype==93) { | |
431 | if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[1]->FillTrackedParticle(vtxMC[2], eta, pt); | |
432 | if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[1] ->FillTrackedParticle(vtxMC[2], eta, pt); | |
433 | } | |
7af955da | 434 | |
85ae901e | 435 | // double diffractive |
436 | if (processtype==94) { | |
437 | if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[2]->FillTrackedParticle(vtxMC[2], eta, pt); | |
438 | if (vertexRecoStudy) fdNdEtaCorrectionVertexReco[2] ->FillTrackedParticle(vtxMC[2], eta, pt); | |
7af955da | 439 | } |
7af955da | 440 | |
85ae901e | 441 | } // end of track loop |
442 | ||
443 | ||
10ebe68d | 444 | delete iter; |
445 | iter = 0; | |
85ae901e | 446 | |
447 | if (fSecondaries) | |
448 | FillSecondaries(); | |
449 | ||
450 | if (fSigmaVertex) | |
451 | FillSigmaVertex(); | |
452 | ||
453 | ||
454 | delete listOfTracks; | |
455 | listOfTracks = 0; | |
456 | ||
457 | return kTRUE; | |
10ebe68d | 458 | } |
459 | ||
7af955da | 460 | void AlidNdEtaSystematicsSelector::FillSecondaries() |
10ebe68d | 461 | { |
462 | // fills the secondary histograms | |
85ae901e | 463 | |
10ebe68d | 464 | AliStack* stack = GetStack(); |
465 | ||
7af955da | 466 | Int_t particlesPrimaries = 0; |
467 | Int_t particlesSecondaries = 0; | |
468 | Int_t particlesPrimariesPtCut = 0; | |
469 | Int_t particlesSecondariesPtCut = 0; | |
10ebe68d | 470 | |
7af955da | 471 | for (Int_t iParticle = 0; iParticle < stack->GetNtrack(); iParticle++) |
10ebe68d | 472 | { |
7af955da | 473 | TParticle* particle = stack->Particle(iParticle); |
474 | if (!particle) | |
475 | { | |
476 | AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (particle loop).", iParticle)); | |
477 | continue; | |
478 | } | |
479 | ||
480 | if (TMath::Abs(particle->Eta()) > 0.9) | |
481 | continue; | |
482 | ||
483 | Bool_t isPrimary = kFALSE; | |
484 | if (iParticle < stack->GetNprimary()) | |
485 | { | |
486 | if (AliPWG0Helper::IsPrimaryCharged(particle, stack->GetNprimary()) == kFALSE) | |
487 | continue; | |
488 | ||
489 | isPrimary = kTRUE; | |
490 | } | |
491 | ||
492 | if (isPrimary) | |
493 | particlesPrimaries++; | |
494 | else | |
495 | particlesSecondaries++; | |
496 | ||
497 | if (particle->Pt() < 0.3) | |
498 | continue; | |
499 | ||
500 | if (isPrimary) | |
501 | particlesPrimariesPtCut++; | |
502 | else | |
503 | particlesSecondariesPtCut++; | |
504 | } | |
505 | ||
506 | fSecondaries->Fill(1, particlesPrimaries); | |
507 | fSecondaries->Fill(2, particlesSecondaries); | |
508 | fSecondaries->Fill(3, particlesPrimariesPtCut); | |
509 | fSecondaries->Fill(4, particlesSecondariesPtCut); | |
510 | ||
511 | Int_t allTracksPrimaries = 0; | |
512 | Int_t allTracksSecondaries = 0; | |
513 | Int_t acceptedTracksPrimaries = 0; | |
514 | Int_t acceptedTracksSecondaries = 0; | |
515 | ||
516 | for (Int_t iTrack = 0; iTrack < fESD->GetNumberOfTracks(); iTrack++) | |
517 | { | |
518 | AliESDtrack* esdTrack = fESD->GetTrack(iTrack); | |
10ebe68d | 519 | if (!esdTrack) |
7af955da | 520 | { |
521 | AliDebug(AliLog::kError, Form("UNEXPECTED: Could not get track %d.", iTrack)); | |
10ebe68d | 522 | continue; |
7af955da | 523 | } |
10ebe68d | 524 | |
525 | Int_t label = TMath::Abs(esdTrack->GetLabel()); | |
10ebe68d | 526 | TParticle* particle = stack->Particle(label); |
527 | if (!particle) | |
528 | { | |
529 | AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label)); | |
530 | continue; | |
531 | } | |
532 | ||
7af955da | 533 | if (label < stack->GetNprimary()) |
534 | allTracksPrimaries++; | |
10ebe68d | 535 | else |
7af955da | 536 | allTracksSecondaries++; |
10ebe68d | 537 | |
7af955da | 538 | if (!fEsdTrackCuts->AcceptTrack(esdTrack)) |
539 | continue; | |
10ebe68d | 540 | |
7af955da | 541 | if (label < stack->GetNprimary()) |
542 | acceptedTracksPrimaries++; | |
543 | else | |
544 | acceptedTracksSecondaries++; | |
10ebe68d | 545 | } |
546 | ||
7af955da | 547 | fSecondaries->Fill(5, allTracksPrimaries); |
548 | fSecondaries->Fill(6, allTracksSecondaries); | |
549 | fSecondaries->Fill(7, acceptedTracksPrimaries); | |
550 | fSecondaries->Fill(8, acceptedTracksSecondaries); | |
10ebe68d | 551 | |
7af955da | 552 | //printf("P = %d, S = %d, P_pt = %d, S_pt = %d, T_P = %d, T_S = %d, T_P_acc = %d, T_S_acc = %d\n", particlesPrimaries, particlesSecondaries, particlesPrimariesPtCut, particlesSecondariesPtCut, allTracksPrimaries, allTracksSecondaries, acceptedTracksPrimaries, acceptedTracksSecondaries); |
10ebe68d | 553 | } |
554 | ||
9f469bf5 | 555 | void AlidNdEtaSystematicsSelector::FillSigmaVertex() |
556 | { | |
557 | // fills the fSigmaVertex histogram | |
558 | ||
559 | // save the old value | |
560 | Float_t oldSigmaVertex = fEsdTrackCuts->GetMinNsigmaToVertex(); | |
561 | ||
562 | // set to maximum | |
563 | fEsdTrackCuts->SetMinNsigmaToVertex(5); | |
564 | ||
565 | TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD); | |
566 | ||
567 | TIterator* iter = list->MakeIterator(); | |
568 | TObject* obj = 0; | |
569 | while ((obj = iter->Next()) != 0) | |
570 | { | |
571 | AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (obj); | |
572 | if (!esdTrack) | |
573 | continue; | |
574 | ||
575 | Float_t sigma2Vertex = fEsdTrackCuts->GetSigmaToVertex(esdTrack); | |
576 | ||
72e597d7 | 577 | for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1) |
9f469bf5 | 578 | { |
579 | if (sigma2Vertex < nSigma) | |
580 | fSigmaVertex->Fill(nSigma); | |
581 | } | |
582 | } | |
583 | ||
584 | delete iter; | |
585 | iter = 0; | |
586 | ||
587 | delete list; | |
588 | list = 0; | |
589 | ||
590 | // set back the old value | |
591 | fEsdTrackCuts->SetMinNsigmaToVertex(oldSigmaVertex); | |
592 | } | |
593 | ||
10ebe68d | 594 | void AlidNdEtaSystematicsSelector::SlaveTerminate() |
595 | { | |
596 | // The SlaveTerminate() function is called after all entries or objects | |
597 | // have been processed. When running with PROOF SlaveTerminate() is called | |
598 | // on each slave server. | |
599 | ||
600 | AliSelector::SlaveTerminate(); | |
601 | ||
602 | // Add the histograms to the output on each slave server | |
603 | if (!fOutput) | |
604 | { | |
605 | AliDebug(AliLog::kError, Form("ERROR: Output list not initialized.")); | |
606 | return; | |
607 | } | |
608 | ||
609 | if (fSecondaries) | |
610 | fOutput->Add(fSecondaries); | |
611 | ||
612 | for (Int_t i=0; i<4; ++i) | |
6de7db91 | 613 | if (fdNdEtaCorrectionSpecies[i]) |
614 | fOutput->Add(fdNdEtaCorrectionSpecies[i]); | |
9f469bf5 | 615 | |
616 | if (fSigmaVertex) | |
617 | fOutput->Add(fSigmaVertex); | |
6de7db91 | 618 | |
7b0956f3 | 619 | for (Int_t i=0; i<3; ++i) { |
6de7db91 | 620 | if (fdNdEtaCorrectionVertexReco[i]) |
621 | fOutput->Add(fdNdEtaCorrectionVertexReco[i]); | |
7b0956f3 | 622 | |
623 | if (fdNdEtaCorrectionTriggerBias[i]) | |
624 | fOutput->Add(fdNdEtaCorrectionTriggerBias[i]); | |
625 | } | |
626 | ||
10ebe68d | 627 | } |
628 | ||
629 | void AlidNdEtaSystematicsSelector::Terminate() | |
630 | { | |
631 | // The Terminate() function is the last function to be called during | |
632 | // a query. It always runs on the client, it can be used to present | |
633 | // the results graphically or save the results to file. | |
634 | ||
635 | AliSelector::Terminate(); | |
636 | ||
7af955da | 637 | fSecondaries = dynamic_cast<TH2F*> (fOutput->FindObject("fSecondaries")); |
10ebe68d | 638 | for (Int_t i=0; i<4; ++i) |
6de7db91 | 639 | fdNdEtaCorrectionSpecies[i] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject(Form("correction_%d", i))); |
9f469bf5 | 640 | fSigmaVertex = dynamic_cast<TH1F*> (fOutput->FindObject("fSigmaVertex")); |
10ebe68d | 641 | |
4c351225 | 642 | fdNdEtaCorrectionVertexReco[0] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("vertexRecoND")); |
643 | fdNdEtaCorrectionVertexReco[1] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("vertexRecoSD")); | |
644 | fdNdEtaCorrectionVertexReco[2] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("vertexRecoDD")); | |
645 | ||
646 | fdNdEtaCorrectionTriggerBias[0] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("triggerBiasND")); | |
647 | fdNdEtaCorrectionTriggerBias[1] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("triggerBiasSD")); | |
648 | fdNdEtaCorrectionTriggerBias[2] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("triggerBiasDD")); | |
649 | ||
0bd1f8a0 | 650 | if (fPIDParticles) |
651 | { | |
652 | TDatabasePDG* pdgDB = new TDatabasePDG; | |
61385583 | 653 | |
0bd1f8a0 | 654 | for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i) |
655 | if (fPIDParticles->GetBinContent(i) > 0) | |
6de7db91 | 656 | printf("PDG = %d (%s): generated: %d, reconstructed: %d, ratio: %f\n", |
657 | (Int_t) fPIDParticles->GetBinCenter(i), | |
658 | pdgDB->GetParticle((Int_t) fPIDParticles->GetBinCenter(i))->GetName(), | |
659 | (Int_t) fPIDParticles->GetBinContent(i), | |
660 | (Int_t) fPIDTracks->GetBinContent(i), | |
661 | ((fPIDTracks->GetBinContent(i) > 0) ? fPIDParticles->GetBinContent(i) / fPIDTracks->GetBinContent(i) : -1)); | |
61385583 | 662 | |
0bd1f8a0 | 663 | delete pdgDB; |
664 | pdgDB = 0; | |
665 | } | |
61385583 | 666 | |
7b0956f3 | 667 | for (Int_t i=0; i<3; ++i) { |
6de7db91 | 668 | if (fdNdEtaCorrectionVertexReco[i]) |
669 | fdNdEtaCorrectionVertexReco[i]->Finish(); | |
7b0956f3 | 670 | |
671 | if (fdNdEtaCorrectionTriggerBias[i]) | |
672 | fdNdEtaCorrectionTriggerBias[i]->Finish(); | |
673 | } | |
6de7db91 | 674 | |
10ebe68d | 675 | TFile* fout = TFile::Open("systematics.root", "RECREATE"); |
676 | ||
677 | if (fEsdTrackCuts) | |
678 | fEsdTrackCuts->SaveHistograms("esd_track_cuts"); | |
679 | ||
680 | if (fSecondaries) | |
681 | fSecondaries->Write(); | |
682 | ||
9f469bf5 | 683 | if (fSigmaVertex) |
684 | fSigmaVertex->Write(); | |
7b0956f3 | 685 | |
10ebe68d | 686 | for (Int_t i=0; i<4; ++i) |
6de7db91 | 687 | if (fdNdEtaCorrectionSpecies[i]) |
688 | fdNdEtaCorrectionSpecies[i]->SaveHistograms(); | |
689 | ||
7b0956f3 | 690 | for (Int_t i=0; i<3; ++i) { |
6de7db91 | 691 | if (fdNdEtaCorrectionVertexReco[i]) |
692 | fdNdEtaCorrectionVertexReco[i]->SaveHistograms(); | |
10ebe68d | 693 | |
7b0956f3 | 694 | if (fdNdEtaCorrectionTriggerBias[i]) |
695 | fdNdEtaCorrectionTriggerBias[i]->SaveHistograms(); | |
696 | } | |
697 | ||
10ebe68d | 698 | fout->Write(); |
699 | fout->Close(); | |
700 | ||
701 | if (fSecondaries) | |
702 | { | |
703 | new TCanvas; | |
7af955da | 704 | fSecondaries->Draw("COLZ"); |
10ebe68d | 705 | } |
9f469bf5 | 706 | |
707 | if (fSigmaVertex) | |
708 | { | |
709 | new TCanvas; | |
710 | fSigmaVertex->Draw(); | |
711 | } | |
10ebe68d | 712 | } |
38faa1b9 | 713 | |
714 | Bool_t AlidNdEtaSystematicsSelector::SignOK(TParticlePDG* particle) | |
715 | { | |
716 | // returns if a particle with this sign should be counted | |
717 | // this is determined by the value of fSignMode, which should have the same sign | |
718 | // as the charge | |
719 | // if fSignMode is 0 all particles are counted | |
720 | ||
721 | if (fSignMode == 0) | |
722 | return kTRUE; | |
723 | ||
724 | if (!particle) | |
725 | { | |
726 | printf("WARNING: not counting a particle that does not have a pdg particle\n"); | |
727 | return kFALSE; | |
728 | } | |
729 | ||
730 | Double_t charge = particle->Charge(); | |
731 | ||
732 | if (fSignMode > 0) | |
733 | if (charge < 0) | |
734 | return kFALSE; | |
735 | ||
736 | if (fSignMode < 0) | |
737 | if (charge > 0) | |
738 | return kFALSE; | |
739 | ||
740 | return kTRUE; | |
741 | } |