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