]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/totEt/AliAnalysisEtMonteCarlo.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
fix warnings
[u/mrichter/AliRoot.git] / PWGLF / totEt / AliAnalysisEtMonteCarlo.cxx
CommitLineData
cf6522d1 1//_________________________________________________________________________
2// Utility Class for transverse energy studies
3// Base class for MC analysis
4// - MC output
5// implementation file
6//
7//*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL)
e2ee5727 8//_________________________________________________________________________
cf6522d1 9
2fbf38ac 10#include "AliAnalysisEtMonteCarlo.h"
11#include "AliAnalysisEtCuts.h"
f61cec2f 12#include "AliAnalysisEtSelector.h"
13#include "AliAnalysisEtSelector.h"
13b0d3c1 14#include "AliESDtrack.h"
2fbf38ac 15#include "AliStack.h"
0651f6b4 16#include "AliVEvent.h"
2fbf38ac 17#include "AliMCEvent.h"
0651f6b4 18#include "AliESDEvent.h"
13b0d3c1 19#include "TH2F.h"
e2ee5727 20#include "TH3F.h"
cf6522d1 21#include "TParticle.h"
ce546038 22#include "AliGenHijingEventHeader.h"
23#include "AliGenPythiaEventHeader.h"
0651f6b4 24#include "TList.h"
25#include "AliESDCaloCluster.h"
0f6416f3 26#include "AliLog.h"
e2ee5727 27#include <iostream>
28#include <AliCentrality.h>
ef647350 29#include "AliPHOSGeoUtils.h"
30#include "AliPHOSGeometry.h"
31#include "TFile.h"
c79c36cc 32#include "AliESDtrackCuts.h"
e2ee5727 33using namespace std;
16abb579 34
35ClassImp(AliAnalysisEtMonteCarlo);
36
37
ce546038 38// ctor
0651f6b4 39AliAnalysisEtMonteCarlo::AliAnalysisEtMonteCarlo():AliAnalysisEt()
c79c36cc 40 ,fIsMC(kTRUE)
f61cec2f 41 ,fImpactParameter(0)
42 ,fNcoll(0)
43 ,fNpart(0)
44 ,fPrimaryTree(0)
45 ,fTotEtWithSecondaryRemoved(0)
46 ,fTotEtSecondaryFromEmEtPrimary(0)
47 ,fTotEtSecondary(0)
48 ,fPrimaryCode(0)
49 ,fPrimaryCharge(0)
50 ,fPrimaryE(0)
51 ,fPrimaryEt(0)
52 ,fPrimaryPx(0)
53 ,fPrimaryPy(0)
54 ,fPrimaryPz(0)
55 ,fPrimaryVx(0)
56 ,fPrimaryVy(0)
57 ,fPrimaryVz(0)
58 ,fPrimaryAccepted(0)
b2c10007 59 ,fPrimaryMatched(0)
f61cec2f 60 ,fDepositedCode(0)
b2c10007 61 ,fDepositedE(0)
f61cec2f 62 ,fDepositedEt(0)
63 ,fDepositedCharge(0)
64 ,fDepositedVx(0)
65 ,fDepositedVy(0)
66 ,fDepositedVz(0)
b2c10007 67 ,fSecondary(kFALSE)
68 ,fReconstructedE(0)
69 ,fReconstructedEt(0)
70 ,fTotPx(0)
71 ,fTotPy(0)
72 ,fTotPz(0)
73 ,fClusterMult(0)
f61cec2f 74 ,fHistDecayVertexNonRemovedCharged(0)
75 ,fHistDecayVertexRemovedCharged(0)
76 ,fHistDecayVertexNonRemovedNeutral(0)
77 ,fHistDecayVertexRemovedNeutral(0)
78
79 ,fHistRemovedOrNot(0)
80 ,fHistEtNonRemovedProtons(0)
81 ,fHistEtNonRemovedAntiProtons(0)
82 ,fHistEtNonRemovedPiPlus(0)
83 ,fHistEtNonRemovedPiMinus(0)
84 ,fHistEtNonRemovedKaonPlus(0)
85 ,fHistEtNonRemovedKaonMinus(0)
86 ,fHistEtNonRemovedK0s(0)
87 ,fHistEtNonRemovedK0L(0)
88 ,fHistEtNonRemovedLambdas(0)
89 ,fHistEtNonRemovedElectrons(0)
90 ,fHistEtNonRemovedPositrons(0)
91 ,fHistEtNonRemovedMuPlus(0)
92 ,fHistEtNonRemovedMuMinus(0)
93 ,fHistEtNonRemovedNeutrons(0)
94 ,fHistEtNonRemovedAntiNeutrons(0)
95 ,fHistEtNonRemovedGammas(0)
96 ,fHistEtNonRemovedGammasFromPi0(0)
97 ,fHistEtRemovedGammas(0)
98 ,fHistEtRemovedNeutrons(0)
99 ,fHistEtRemovedAntiNeutrons(0)
100 ,fHistEtRemovedCharged(0)
101 ,fHistEtRemovedNeutrals(0)
102 ,fHistEtNonRemovedCharged(0)
103 ,fHistEtNonRemovedNeutrals(0)
104 ,fHistMultNonRemovedProtons(0)
105 ,fHistMultNonRemovedAntiProtons(0)
106 ,fHistMultNonRemovedPiPlus(0)
107 ,fHistMultNonRemovedPiMinus(0)
108 ,fHistMultNonRemovedKaonPlus(0)
109 ,fHistMultNonRemovedKaonMinus(0)
110 ,fHistMultNonRemovedK0s(0)
111 ,fHistMultNonRemovedK0L(0)
112 ,fHistMultNonRemovedLambdas(0)
113 ,fHistMultNonRemovedElectrons(0)
114 ,fHistMultNonRemovedPositrons(0)
115 ,fHistMultNonRemovedMuPlus(0)
116 ,fHistMultNonRemovedMuMinus(0)
117 ,fHistMultNonRemovedNeutrons(0)
118 ,fHistMultNonRemovedAntiNeutrons(0)
119 ,fHistMultNonRemovedGammas(0)
120 ,fHistMultRemovedGammas(0)
121 ,fHistMultRemovedNeutrons(0)
122 ,fHistMultRemovedAntiNeutrons(0)
123 ,fHistMultRemovedCharged(0)
124 ,fHistMultRemovedNeutrals(0)
125 ,fHistMultNonRemovedCharged(0)
126 ,fHistMultNonRemovedNeutrals(0)
127 ,fHistTrackMultvsNonRemovedCharged(0)
128 ,fHistTrackMultvsNonRemovedNeutral(0)
129 ,fHistTrackMultvsRemovedGamma(0)
130 ,fHistClusterMultvsNonRemovedCharged(0)
131 ,fHistClusterMultvsNonRemovedNeutral(0)
132 ,fHistClusterMultvsRemovedGamma(0)
133 ,fHistMultvsNonRemovedChargedE(0)
134 ,fHistMultvsNonRemovedNeutralE(0)
03bfd268 135 ,fHistMultvsRemovedGammaE(0)
0861cc1f 136 ,fCalcForKaonCorrection(kFALSE)
137 ,fHistK0EDepositsVsPtInAcceptance(0)
138 ,fHistK0EGammaVsPtInAcceptance(0)
139 ,fHistK0EDepositsVsPtOutOfAcceptance(0)
140 ,fHistK0EGammaVsPtOutOfAcceptance(0)
141 ,fHistSimKaonsInAcceptance(0)
142 ,fHistSimK0SInAcceptance(0)
143 ,fHistSimKPlusInAcceptance(0)
144 ,fHistSimKMinusInAcceptance(0)
145 ,fHistSimK0LInAcceptance(0)
146 ,fHistSimKaonsOutOfAcceptance(0)
147 ,fHistSimKaonsInAcceptanceWithDepositsPrimaries(0)
148 ,fHistSimKaonsOutOfAcceptanceWithDepositsSecondaries(0)
149 ,fHistSimKaonsOutOfAcceptanceWithDepositsPrimaries(0)
f61cec2f 150 ,fEtNonRemovedProtons(0)
151 ,fEtNonRemovedAntiProtons(0)
152 ,fEtNonRemovedPiPlus(0)
153 ,fEtNonRemovedPiMinus(0)
154 ,fEtNonRemovedKaonPlus(0)
155 ,fEtNonRemovedKaonMinus(0)
156 ,fEtNonRemovedK0S(0)
157 ,fEtNonRemovedK0L(0)
158 ,fEtNonRemovedLambdas(0)
159 ,fEtNonRemovedElectrons(0)
160 ,fEtNonRemovedPositrons(0)
161 ,fEtNonRemovedMuMinus(0)
162 ,fEtNonRemovedMuPlus(0)
163 ,fEtNonRemovedGammas(0)
164 ,fEtNonRemovedGammasFromPi0(0)
165 ,fEtNonRemovedNeutrons(0)
166 ,fEtNonRemovedAntiNeutrons(0)
167 ,fEtRemovedProtons(0)
168 ,fEtRemovedAntiProtons(0)
169 ,fEtRemovedPiPlus(0)
170 ,fEtRemovedPiMinus(0)
171 ,fEtRemovedKaonPlus(0)
172 ,fEtRemovedKaonMinus(0)
173 ,fEtRemovedK0s(0)
174 ,fEtRemovedK0L(0)
175 ,fEtRemovedLambdas(0)
176 ,fEtRemovedElectrons(0)
177 ,fEtRemovedPositrons(0)
178 ,fEtRemovedMuMinus(0)
179 ,fEtRemovedMuPlus(0)
180 ,fEtRemovedGammasFromPi0(0)
181 ,fEtRemovedGammas(0)
182 ,fEtRemovedNeutrons(0)
183 ,fEtRemovedAntiNeutrons(0)
184 ,fMultNonRemovedProtons(0)
185 ,fMultNonRemovedAntiProtons(0)
186 ,fMultNonRemovedPiPlus(0)
187 ,fMultNonRemovedPiMinus(0)
188 ,fMultNonRemovedKaonPlus(0)
189 ,fMultNonRemovedKaonMinus(0)
190 ,fMultNonRemovedK0s(0)
191 ,fMultNonRemovedK0L(0)
192 ,fMultNonRemovedLambdas(0)
193 ,fMultNonRemovedElectrons(0)
194 ,fMultNonRemovedPositrons(0)
195 ,fMultNonRemovedMuMinus(0)
196 ,fMultNonRemovedMuPlus(0)
197 ,fMultNonRemovedGammas(0)
198 ,fMultNonRemovedNeutrons(0)
199 ,fMultNonRemovedAntiNeutrons(0)
200 ,fMultRemovedProtons(0)
201 ,fMultRemovedAntiProtons(0)
202 ,fMultRemovedPiPlus(0)
203 ,fMultRemovedPiMinus(0)
204 ,fMultRemovedKaonPlus(0)
205 ,fMultRemovedKaonMinus(0)
206 ,fMultRemovedK0s(0)
207 ,fMultRemovedK0L(0)
208 ,fMultRemovedLambdas(0)
209 ,fMultRemovedElectrons(0)
210 ,fMultRemovedPositrons(0)
211 ,fMultRemovedMuMinus(0)
212 ,fMultRemovedMuPlus(0)
213 ,fMultRemovedGammas(0)
214 ,fMultRemovedNeutrons(0)
215 ,fMultRemovedAntiNeutrons(0)
216 ,fTrackMultInAcc(0)
217 ,fHistDxDzNonRemovedCharged(0)
218 ,fHistDxDzRemovedCharged(0)
219 ,fHistDxDzNonRemovedNeutral(0)
220 ,fHistDxDzRemovedNeutral(0)
221 ,fHistPiPlusMult(0)
222 ,fHistPiMinusMult(0)
223 ,fHistPiZeroMult(0)
224 ,fHistPiPlusMultAcc(0)
225 ,fHistPiMinusMultAcc(0)
226 ,fHistPiZeroMultAcc(0)
b2c10007 227// ,fPiPlusMult(0)
228// ,fPiMinusMult(0)
f61cec2f 229 ,fPiZeroMult(0)
230 ,fPiPlusMultAcc(0)
231 ,fPiMinusMultAcc(0)
232 ,fPiZeroMultAcc(0)
233 ,fNeutralRemoved(0)
234 ,fChargedRemoved(0)
235 ,fChargedNotRemoved(0)
236 ,fNeutralNotRemoved(0)
b2c10007 237 ,fGammaRemoved(0)
238 ,fSecondaryNotRemoved(0)
f61cec2f 239 ,fEnergyNeutralRemoved(0)
240 ,fEnergyChargedRemoved(0)
241 ,fEnergyChargedNotRemoved(0)
242 ,fEnergyNeutralNotRemoved(0)
b2c10007 243 ,fEnergyGammaRemoved(0)
f61cec2f 244 ,fNClusters(0)
245 ,fTotNeutralEtAfterMinEnergyCut(0)
0861cc1f 246 ,fCalcTrackMatchVsMult(kFALSE)
247 ,fHistGammasFound(0)
9ef6f13f 248 ,fHistGammasGenerated(0)
6a152780 249 ,fHistGammasFoundMult(0)
250 ,fHistGammasGeneratedMult(0)
9ef6f13f 251 ,fHistChargedTracksCut(0)
252 ,fHistChargedTracksAccepted(0)
253 ,fHistGammasCut(0)
254 ,fHistGammasAccepted(0)
0861cc1f 255 ,fHistChargedTracksCutMult(0)
256 ,fHistChargedTracksAcceptedMult(0)
ac610b08 257 ,fHistChargedTracksAcceptedLowPtCent(0)
258 ,fHistChargedTracksAcceptedLowPtCentNoAntiProtons(0)
259 ,fHistChargedTracksAcceptedLowPtCentAntiProtons(0)
0861cc1f 260 ,fHistGammasCutMult(0)
261 ,fHistGammasAcceptedMult(0)
9ef6f13f 262 ,fHistBadTrackMatches(0)
263 ,fHistMatchedTracksEvspTBkgd(0)
264 ,fHistMatchedTracksEvspTSignal(0)
265 ,fHistMatchedTracksEvspTBkgdPeripheral(0)
266 ,fHistMatchedTracksEvspTSignalPeripheral(0)
ac610b08 267 ,fHistMatchedTracksEvspTBkgdvsCent(0)
268 ,fHistMatchedTracksEvspTSignalvsCent(0)
269 ,fHistMatchedTracksEvspTBkgdvsCentEffCorr(0)
270 ,fHistMatchedTracksEvspTSignalvsCentEffCorr(0)
0861cc1f 271
9ef6f13f 272 ,fHistChargedTracksCutPeripheral(0)
273 ,fHistChargedTracksAcceptedPeripheral(0)
274 ,fHistGammasCutPeripheral(0)
275 ,fHistGammasAcceptedPeripheral(0)
276 ,fHistBadTrackMatchesdPhidEta(0)
277 ,fHistGoodTrackMatchesdPhidEta(0)
c79c36cc 278 ,fHistHadronDepositsAll(0)
279 ,fHistHadronDepositsReco(0)
ac610b08 280 ,fHistHadronDepositsAllCent(0)
281 ,fHistHadronDepositsRecoCent(0)
282 ,fHistHadronsAllCent(0)
d3ce32b8 283 ,fHistMultChVsSignalVsMult(0)
ce546038 284{
285}
286
287// dtor
e2ee5727 288AliAnalysisEtMonteCarlo::~AliAnalysisEtMonteCarlo()
f61cec2f 289{ //Destructor
311c6540 290
291 if(fPrimaryTree){
292 fPrimaryTree->Clear();
293 delete fPrimaryTree;
294 }
f61cec2f 295 delete fHistDecayVertexNonRemovedCharged; // Decay vertex for non-removed charged particles
296 delete fHistDecayVertexRemovedCharged; // Decay vertex for non-removed charged particles
297 delete fHistDecayVertexNonRemovedNeutral; // Decay vertex for non-removed charged particles
298 delete fHistDecayVertexRemovedNeutral; // Decay vertex for non-removed charged particles
299
300 delete fHistRemovedOrNot; // If charged/neutral particles were removed or not
301
302 delete fHistEtNonRemovedProtons; // enter comment here
303 delete fHistEtNonRemovedAntiProtons; // enter comment here
304 delete fHistEtNonRemovedPiPlus; // enter comment here
305 delete fHistEtNonRemovedPiMinus; // enter comment here
306 delete fHistEtNonRemovedKaonPlus; // enter comment here
307 delete fHistEtNonRemovedKaonMinus; // enter comment here
308 delete fHistEtNonRemovedK0s; // enter comment here
309 delete fHistEtNonRemovedK0L; // enter comment here
310 delete fHistEtNonRemovedLambdas; // enter comment here
311 delete fHistEtNonRemovedElectrons; // enter comment here
312 delete fHistEtNonRemovedPositrons; // enter comment here
313 delete fHistEtNonRemovedMuPlus; // enter comment here
314 delete fHistEtNonRemovedMuMinus; // enter comment here
315 delete fHistEtNonRemovedNeutrons; // enter comment here
316 delete fHistEtNonRemovedAntiNeutrons; // enter comment here
317 delete fHistEtNonRemovedGammas; // enter comment here
318 delete fHistEtNonRemovedGammasFromPi0; // enter comment here
319
320 delete fHistEtRemovedGammas; // enter comment here
321 delete fHistEtRemovedNeutrons; // enter comment here
322 delete fHistEtRemovedAntiNeutrons; // enter comment here
323
324
325 delete fHistMultNonRemovedProtons; // enter comment here
326 delete fHistMultNonRemovedAntiProtons; // enter comment here
327 delete fHistMultNonRemovedPiPlus; // enter comment here
328 delete fHistMultNonRemovedPiMinus; // enter comment here
329 delete fHistMultNonRemovedKaonPlus; // enter comment here
330 delete fHistMultNonRemovedKaonMinus; // enter comment here
331 delete fHistMultNonRemovedK0s; // enter comment here
332 delete fHistMultNonRemovedK0L; // enter comment here
333 delete fHistMultNonRemovedLambdas; // enter comment here
334 delete fHistMultNonRemovedElectrons; // enter comment here
335 delete fHistMultNonRemovedPositrons; // enter comment here
336 delete fHistMultNonRemovedMuPlus; // enter comment here
337 delete fHistMultNonRemovedMuMinus; // enter comment here
338 delete fHistMultNonRemovedNeutrons; // enter comment here
339 delete fHistMultNonRemovedAntiNeutrons; // enter comment here
340 delete fHistMultNonRemovedGammas; // enter comment here
341
342 delete fHistMultRemovedGammas; // enter comment here
343 delete fHistMultRemovedNeutrons; // enter comment here
344 delete fHistMultRemovedAntiNeutrons; // enter comment here
345
346 delete fHistTrackMultvsNonRemovedCharged; // enter comment here
347 delete fHistTrackMultvsNonRemovedNeutral; // enter comment here
348 delete fHistTrackMultvsRemovedGamma; // enter comment here
349
350 delete fHistClusterMultvsNonRemovedCharged; // enter comment here
351 delete fHistClusterMultvsNonRemovedNeutral; // enter comment here
352 delete fHistClusterMultvsRemovedGamma; // enter comment here
353
354 delete fHistMultvsNonRemovedChargedE; // enter comment here
355 delete fHistMultvsNonRemovedNeutralE; // enter comment here
356 delete fHistMultvsRemovedGammaE; // enter comment here
0861cc1f 357 delete fHistK0EDepositsVsPtInAcceptance;//enter comment here
358 delete fHistK0EGammaVsPtInAcceptance;//enter comment here
359 delete fHistK0EDepositsVsPtOutOfAcceptance;
360 delete fHistK0EGammaVsPtOutOfAcceptance;
361
362 delete fHistSimKaonsInAcceptance;
363 delete fHistSimK0SInAcceptance;
364 delete fHistSimKPlusInAcceptance;
365 delete fHistSimKMinusInAcceptance;
366 delete fHistSimK0LInAcceptance;
367 delete fHistSimKaonsOutOfAcceptance;
368 delete fHistSimKaonsInAcceptanceWithDepositsPrimaries;
369 delete fHistSimKaonsOutOfAcceptanceWithDepositsSecondaries;
370 delete fHistSimKaonsOutOfAcceptanceWithDepositsPrimaries;
2aab9269 371
f61cec2f 372 delete fHistDxDzNonRemovedCharged; // enter comment here
373 delete fHistDxDzRemovedCharged; // enter comment here
374 delete fHistDxDzNonRemovedNeutral; // enter comment here
375 delete fHistDxDzRemovedNeutral; // enter comment here
376
377 delete fHistPiPlusMult; // enter comment here
378 delete fHistPiMinusMult; // enter comment here
379 delete fHistPiZeroMult; // enter comment here
380
381 delete fHistPiPlusMultAcc; // enter comment here
382 delete fHistPiMinusMultAcc; // enter comment here
383 delete fHistPiZeroMultAcc; // enter comment here
4503e29d 384 delete fHistGammasFound; // enter comment here
385 delete fHistGammasGenerated; // enter comment here
6a152780 386 delete fHistGammasFoundMult; // enter comment here
387 delete fHistGammasGeneratedMult; // enter comment here
9ef6f13f 388 delete fHistChargedTracksCut;
389 delete fHistChargedTracksAccepted;
390 delete fHistGammasCut;
391 delete fHistGammasAccepted;
0861cc1f 392 delete fHistChargedTracksCutMult;
393 delete fHistChargedTracksAcceptedMult;
ac610b08 394 delete fHistChargedTracksAcceptedLowPtCent;
395 delete fHistChargedTracksAcceptedLowPtCentNoAntiProtons;
396 delete fHistChargedTracksAcceptedLowPtCentAntiProtons;
0861cc1f 397 delete fHistGammasCutMult;
398 delete fHistGammasAcceptedMult;
9ef6f13f 399 delete fHistBadTrackMatches;
400 delete fHistMatchedTracksEvspTBkgd;
401 delete fHistMatchedTracksEvspTSignal;
402 delete fHistMatchedTracksEvspTBkgdPeripheral;
403 delete fHistMatchedTracksEvspTSignalPeripheral;
ac610b08 404 delete fHistMatchedTracksEvspTBkgdvsCent;
405 delete fHistMatchedTracksEvspTSignalvsCent;
406 delete fHistMatchedTracksEvspTBkgdvsCentEffCorr;
407 delete fHistMatchedTracksEvspTSignalvsCentEffCorr;
9ef6f13f 408 delete fHistChargedTracksCutPeripheral;
409 delete fHistChargedTracksAcceptedPeripheral;
410 delete fHistGammasCutPeripheral;
411 delete fHistGammasAcceptedPeripheral;
412 delete fHistBadTrackMatchesdPhidEta;
413 delete fHistGoodTrackMatchesdPhidEta;
c79c36cc 414 delete fHistHadronDepositsAll;
415 delete fHistHadronDepositsReco;
ac610b08 416 delete fHistHadronDepositsAllCent;
417 delete fHistHadronDepositsRecoCent;
418 delete fHistHadronsAllCent;
d3ce32b8 419 delete fHistMultChVsSignalVsMult;
ce546038 420}
2fbf38ac 421
422Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
f61cec2f 423{ // analyse MC event
c79c36cc 424 if(!fIsMC) return 0;
e2ee5727 425 ResetEventValues();
426
b2c10007 427
e2ee5727 428 fPiPlusMult = 0;
429 fPiMinusMult = 0;
430 fPiZeroMult = 0;
431 if (fCentrality)
432 {
ac610b08 433 fCentClass = fCentrality->GetCentralityClass5(fCentralityMethod);
e2ee5727 434
435 }
e2ee5727 436
437 // Get us an mc event
438 if (!ev) {
439 AliFatal("ERROR: Event does not exist");
440 return 0;
441 }
442 AliMCEvent *event = dynamic_cast<AliMCEvent*>(ev);
443 if (!event) {
444 AliFatal("ERROR: MC Event does not exist");
445 return 0;
446 }
447
448 Double_t protonMass =fgProtonMass;
449
450 // Hijing header
451 AliGenEventHeader* genHeader = event->GenEventHeader();
452 if (!genHeader) {
453 Printf("ERROR: Event generation header does not exist");
454 return 0;
455 }
456 AliGenHijingEventHeader* hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
457 if (hijingGenHeader) {
458 fImpactParameter = hijingGenHeader->ImpactParameter();
459 fNcoll = hijingGenHeader->HardScatters(); // or should this be some combination of NN() NNw() NwN() NwNw() ?
460 fNpart = hijingGenHeader->ProjectileParticipants() + hijingGenHeader->TargetParticipants();
f61cec2f 461 }
462
e2ee5727 463 // Let's play with the stack!
464 AliStack *stack = event->Stack();
465
466 Int_t nPrim = stack->GetNtrack();
467
468 Int_t partCount = 0;
469 for (Int_t iPart = 0; iPart < nPrim; iPart++)
2fbf38ac 470 {
e2ee5727 471
472 TParticle *part = stack->Particle(iPart);
b2c10007 473
474
e2ee5727 475
476 if (!part)
477 {
478 Printf("ERROR: Could not get particle %d", iPart);
479 continue;
480 }
e2ee5727 481 TParticlePDG *pdg = part->GetPDG(0);
e2ee5727 482
f61cec2f 483 if (!pdg)
2fbf38ac 484 {
f61cec2f 485 Printf("ERROR: Could not get particle PDG %d", iPart);
e2ee5727 486 continue;
2fbf38ac 487 }
488
e2ee5727 489 Double_t particleMassPart = 0; //The mass part in the Et calculation for this particle
f61cec2f 490 Int_t code = pdg->PdgCode();
b2c10007 491
492 if(stack->IsPhysicalPrimary(iPart))
493 {
494 fTotPx += part->Px();
495 fTotPy += part->Py();
496 fTotPz += part->Pz();
497 }
498
e2ee5727 499 // Check for reasonable (for now neutral and singly charged) charge on the particle
86e7d5db 500 if (fSelector->IsNeutralMcParticle(iPart,*stack,*pdg))
e2ee5727 501 {
502
503 fMultiplicity++;
b2c10007 504// PrintFamilyTree(iPart, stack);
f61cec2f 505//
e2ee5727 506 if (TMath::Abs(part->Eta()) < fCuts->GetCommonEtaCut())
2fbf38ac 507 {
f61cec2f 508 //Printf("Particle with eta: %f, pid: %d", part->Eta(), code);
e2ee5727 509 // calculate E_T
510 if (
ef647350 511 TMath::Abs(code) == fgProtonCode ||
512 TMath::Abs(code) == fgNeutronCode ||
513 TMath::Abs(code) == fgLambdaCode ||
514 TMath::Abs(code) == fgXiCode ||
515 TMath::Abs(code) == fgXi0Code ||
516 TMath::Abs(code) == fgOmegaCode
e2ee5727 517 )
518 {
ef647350 519 if (code > 0) {
e2ee5727 520 particleMassPart = - protonMass;
521 }
ef647350 522 if (code < 0) {
e2ee5727 523 particleMassPart = protonMass;
524 }
525 }
526 Double_t et = part->Energy() * TMath::Sin(part->Theta()) + particleMassPart;
527
e2ee5727 528
6a152780 529// // Fill up total E_T counters for each particle species
530// if (code == fgProtonCode || code == fgAntiProtonCode)
531// {
532// }
533// if (code == fgPiPlusCode || code == fgPiMinusCode)
534// {
535// if (code == fgPiPlusCode)
536// {
537// }
538// else
539// {
540// }
541// }
542// if (code == fgGammaCode)
543// {
544// }
545// if (code == fgKPlusCode)
546// {
547// }
548// if(code == fgKMinusCode)
549// {
550// }
551// if (code == fgMuPlusCode || code == fgMuMinusCode)
552// {
553// }
554// if (code == fgEPlusCode || code == fgEMinusCode)
555// {
556// }
557// // some neutrals also
558// if (code == fgNeutronCode)
559// {
560// }
561// if (code == fgAntiNeutronCode)
562// {
563// }
564// if (code == fgGammaCode)
565// {
566// }
e2ee5727 567
568 // Neutral particles
ef647350 569 //if (TMath::Abs(pdg->Charge() - fCuts->GetMonteCarloNeutralParticle()) <1e-3 )
f61cec2f 570
ef647350 571 if(code == fgGammaCode || code == fgPi0Code || code == fgEtaCode)
e2ee5727 572 {
b2c10007 573 // PrintFamilyTree(iPart,stack);
f61cec2f 574 //Printf("Gamma, phi: %f, eta: %f, phi cut min: %f, phi cut max: %f, eta cut: %f", part->Phi(), part->Eta(), fPhiMinCutAcc, fPhiMaxCutAcc, fEtaCutAcc);
575 //if (et > fCuts->GetCommonClusterEnergyCut()) fTotNeutralEt += et;
e2ee5727 576
577 // inside EMCal acceptance
f61cec2f 578
579 //if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiMaxCutAcc && part->Phi() > fPhiMinCutAcc)
580
b2c10007 581 if(fSelector->CutGeometricalAcceptance(*part) )
e2ee5727 582 {
583 fNeutralMultiplicity++;
f61cec2f 584 fTotNeutralEt += et;
86e7d5db 585 if(fSelector->PassMinEnergyCut(*part))
f61cec2f 586 {
587 fTotNeutralEtAfterMinEnergyCut += et;
588 }
e2ee5727 589 if (part->Energy() > 0.05) partCount++;
e2ee5727 590 }
591 }
592 //Charged particles
593 else if (TMath::Abs( pdg->Charge() - fCuts->GetMonteCarloNeutralParticle())>1e-3 )
594 {
e2ee5727 595
596 // inside EMCal acceptance
f61cec2f 597 if (fSelector->CutGeometricalAcceptance(*part))
e2ee5727 598 {
f61cec2f 599
600 fChargedMultiplicity++;
601
602 fTotChargedEt += et;
e2ee5727 603
ef647350 604 if (code == fgProtonCode || code == fgAntiProtonCode)
e2ee5727 605 {
e2ee5727 606 }
ef647350 607 if (code == fgPiPlusCode || code == fgPiMinusCode)
e2ee5727 608 {
e2ee5727 609 }
ef647350 610 if (code == fgKPlusCode || code == fgKMinusCode)
e2ee5727 611 {
e2ee5727 612 }
ef647350 613 if (code == fgMuPlusCode || code == fgMuMinusCode)
e2ee5727 614 {
e2ee5727 615 }
f61cec2f 616
ef647350 617 if (code == fgEPlusCode || code == fgEMinusCode)
e2ee5727 618 {
f61cec2f 619 fTotNeutralEt += et; // calling electrons neutral
620 fTotChargedEt -= et;
e1fa1966 621 }
e2ee5727 622 } // inside EMCal acceptance
623
e2ee5727 624 if (TrackHitsCalorimeter(part)) // magnetic field info not filled?
625 {
626 if (pdg->Charge() > 0) fHistPhivsPtPos->Fill(part->Phi(),part->Pt());
627 else fHistPhivsPtNeg->Fill(part->Phi(), part->Pt());
628 }
629 }
630 }
631
632 }
2fbf38ac 633 }
b2c10007 634 // std::cout << "Total: p_x = " << fTotPx << ", p_y = " << fTotPy << ", p_z = " << fTotPz << std::endl;
e2ee5727 635 fTotEt = fTotChargedEt + fTotNeutralEt;
f61cec2f 636 //fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;//
e2ee5727 637 //std::cout << "Event done! # of particles: " << partCount << std::endl;
638 return 0;
2fbf38ac 639}
0651f6b4 640//Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliMCEvent* mcEvent,AliESDEvent* realEvent)
641Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
f61cec2f 642{ // analyse MC and real event info
e2ee5727 643 //if(!mcEvent || !realEvent){
c79c36cc 644 if(!fIsMC) return 0;
e2ee5727 645 if (!ev || !ev2) {
646 AliFatal("ERROR: Event does not exist");
647 return 0;
648 }
f61cec2f 649 AliAnalysisEt::AnalyseEvent(ev);
e2ee5727 650 AliMCEvent *mcEvent = dynamic_cast<AliMCEvent*>(ev);
651 AliESDEvent *realEvent = dynamic_cast<AliESDEvent*>(ev2);
652 if (!mcEvent || !realEvent) {
653 AliFatal("ERROR: mcEvent or realEvent does not exist");
b2c10007 654
e2ee5727 655 }
656
b2c10007 657 std::vector<Int_t> foundGammas;
658
f61cec2f 659 fSelector->SetEvent(realEvent);
660
661 AnalyseEvent(ev);
662
e2ee5727 663 AliStack *stack = mcEvent->Stack();
664
9ef6f13f 665 // get all detector clusters
f61cec2f 666 // TRefArray* caloClusters = new TRefArray();
667
668// if (fDetector == fCuts->GetDetectorEmcal()) realEvent->GetEMCALClusters( caloClusters );
669 //else if (fDetector == fCuts->GetDetectorPhos()) realEvent->GetPHOSClusters( caloClusters );
670 //else {
671 //AliFatal("Detector ID has not been specified");
672 //return -1;
673// }
674
4503e29d 675//Note that this only returns clusters for the selected detector. fSelector actually calls the right GetClusters... for the detector
cfe23ff0 676//It does not apply any cuts on these clusters
f61cec2f 677 TRefArray *caloClusters = fSelector->GetClusters();
e2ee5727 678
e2ee5727 679 Int_t nCluster = caloClusters->GetEntries();
b2c10007 680 fClusterMult = nCluster;
ef647350 681 fNClusters = 0;
d3ce32b8 682 Int_t fClusterMultChargedTracks = 0;
683 Int_t fClusterMultGammas = 0;
e2ee5727 684 // loop the clusters
685 for (int iCluster = 0; iCluster < nCluster; iCluster++ )
686 {
f61cec2f 687 Int_t cf = 0;
e2ee5727 688 AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
2aab9269 689 //Float_t caloE = caloCluster->E()
f61cec2f 690 fNClusters++;
b2c10007 691 const UInt_t iPart = (UInt_t)TMath::Abs(caloCluster->GetLabel());
ef647350 692 TParticle *part = stack->Particle(iPart);
e2ee5727 693
694 if (!part)
0651f6b4 695 {
e2ee5727 696 Printf("No MC particle %d", iCluster);
697 continue;
698 }
9ef6f13f 699
f61cec2f 700 int primIdx = iPart;
e2ee5727 701 if (!stack->IsPhysicalPrimary(iPart)) // check whether particle is primary. we keep secondary electron and gamma for testing.
702 {
f61cec2f 703 primIdx = GetPrimMother(iPart, stack);
b2c10007 704 if(primIdx != stack->GetPrimary(iPart))
705 {
706 //std::cout << primIdx << " = " << stack->GetPrimary(iPart) << std::endl;
707 //PrintFamilyTree(iPart, stack);
708 }
2aab9269 709 //if it is from a K0S
f61cec2f 710 if(primIdx < 0)
711 {
ac610b08 712 //std::cout << "What!? No primary?" << std::endl;
713 //PrintFamilyTree(iPart, stack);
9ef6f13f 714 //continue;
715 //This is a work around to fix a bug. For the EMCal when you use the tender supply, the parent particle ID gets messed up.
716 primIdx = iPart;
f61cec2f 717 }
718
e2ee5727 719 } // end of primary particle check
a0dcb3c1 720 //const int primCode = stack->Particle(primIdx)->GetPdgCode();
ef647350 721 TParticlePDG *pdg = part->GetPDG();
f61cec2f 722 if (!pdg)
ef647350 723 {
724 Printf("ERROR: Could not get particle PDG %d", iPart);
725 continue;
f61cec2f 726 }
e2ee5727 727
a0dcb3c1 728 //Int_t code = pdg->PdgCode();
53ad5ec1 729// if(primCode == fgGammaCode)
730// {
731
d3ce32b8 732
733
734 Bool_t nottrackmatched = kTRUE;//default to no track matched
6a152780 735 Float_t matchedTrackp = 0.0;
736 Float_t matchedTrackpt = 0.0;
d3ce32b8 737 nottrackmatched = fSelector->PassTrackMatchingCut(*caloCluster);
738 //by default ALL matched tracks are accepted, whether or not the match is good. So we check to see if the track is good.
584f2478 739 if(!nottrackmatched){//if the track is trackmatched
d3ce32b8 740 Int_t trackMatchedIndex = caloCluster->GetTrackMatchedIndex();
741 if(trackMatchedIndex < 0) nottrackmatched=kTRUE;
742 AliESDtrack *track = realEvent->GetTrack(trackMatchedIndex);
6a152780 743 matchedTrackp = track->P();
744 matchedTrackpt = track->Pt();
d3ce32b8 745 //if this is a good track, accept track will return true. The track matched is good, so not track matched is false
584f2478 746 nottrackmatched = !(fEsdtrackCutsTPC->AcceptTrack(track));//if the track is bad, this is track matched
6a152780 747 //if(!nottrackmatched) cout<<"Matched track p: "<<matchedTrackp<<" sim "<<part->P()<<endl;
d3ce32b8 748 }
749
750
a0dcb3c1 751 for(UInt_t i = 0; i < caloCluster->GetNLabels(); i++)
b2c10007 752 {
53ad5ec1 753 Int_t pIdx = caloCluster->GetLabelAt(i);
0861cc1f 754 //TParticle *p = stack->Particle(pIdx);
755
53ad5ec1 756 if(!stack->IsPhysicalPrimary(pIdx))
757 {
0861cc1f 758// PrintFamilyTree(pIdx, stack);
53ad5ec1 759 pIdx = GetPrimMother(pIdx, stack);
760 }
0861cc1f 761 if(fSelector->PassDistanceToBadChannelCut(*caloCluster))//&&fSelector->CutGeometricalAcceptance(*(stack->Particle(primIdx))))
762 {
d3ce32b8 763
0861cc1f 764// std::cout << "Gamma primary: " << primIdx << std::endl;
765// foundGammas.push_back(primIdx);
d3ce32b8 766 if(nottrackmatched){
767 foundGammas.push_back(pIdx);
768 }
0861cc1f 769 }
b2c10007 770 }
53ad5ec1 771 fCutFlow->Fill(cf++);
86e7d5db 772 if(!fSelector->PassDistanceToBadChannelCut(*caloCluster)) continue;
47151f26 773 Double_t clEt = CorrectForReconstructionEfficiency(*caloCluster,fClusterMult);
ef647350 774// if(code == fgK0SCode) std::cout << "K0 energy: " << caloCluster->E() << std::endl;
86e7d5db 775 if(!fSelector->PassMinEnergyCut(*caloCluster)) continue;
53ad5ec1 776
777
f61cec2f 778 fCutFlow->Fill(cf++);
0861cc1f 779 Float_t pos[3];
780 //PrintFamilyTree(
781 caloCluster->GetPosition(pos);
782 TVector3 cp(pos);
783
f61cec2f 784 TParticle *primPart = stack->Particle(primIdx);
785 fPrimaryCode = primPart->GetPdgCode();
03bfd268 786 fPrimaryCharge = (Int_t) primPart->GetPDG()->Charge();
f61cec2f 787
788 fPrimaryE = primPart->Energy();
789 fPrimaryEt = primPart->Energy()*TMath::Sin(primPart->Theta());
f61cec2f 790 fPrimaryPx = primPart->Px();
791 fPrimaryPy = primPart->Py();
792 fPrimaryPz = primPart->Pz();
cfe23ff0 793 //cout<<"I have a cluster and it's good energy "<<caloCluster->E()<<" simulated "<<fPrimaryE<<endl;
f61cec2f 794 fPrimaryVx = primPart->Vx();
795 fPrimaryVy = primPart->Vy();
796 fPrimaryVz = primPart->Vz();
797
798 fPrimaryAccepted = false;
b2c10007 799 fPrimaryMatched = false;
f61cec2f 800
801 fDepositedCode = part->GetPdgCode();
b2c10007 802 fDepositedE = part->Energy();
803 fDepositedEt = part->Energy()*TMath::Sin(part->Theta());
03bfd268 804 fDepositedCharge = (Int_t) part->GetPDG()->Charge();
f61cec2f 805
806 fDepositedVx = part->Vx();
807 fDepositedVy = part->Vy();
808 fDepositedVz = part->Vz();
b2c10007 809
810 //fSecondary = fSelector->FromSecondaryInteraction(*primPart, *stack);
811 fSecondary =fSelector->FromSecondaryInteraction(*part, *stack);
9ef6f13f 812// if(fSecondary)
813// {
814// //std::cout << "Have secondary!" << std::endl;
815// //PrintFamilyTree(iPart, stack);
816// }
b2c10007 817 fReconstructedE = caloCluster->E();
818 fReconstructedEt = caloCluster->E()*TMath::Sin(cp.Theta());
b2c10007 819
820 pdg = primPart->GetPDG(0);
a0dcb3c1 821 //Int_t code = primPart->GetPdgCode();
e2ee5727 822
9ef6f13f 823 Bool_t written = kFALSE;
824
d3ce32b8 825// Bool_t nottrackmatched = kTRUE;//default to no track matched
826// nottrackmatched = fSelector->PassTrackMatchingCut(*caloCluster);
827// //by default ALL matched tracks are accepted, whether or not the match is good. So we check to see if the track is good.
828// if(!nottrackmatched){
829// Int_t trackMatchedIndex = caloCluster->GetTrackMatchedIndex();
830// if(trackMatchedIndex < 0) nottrackmatched=kTRUE;
831// AliESDtrack *track = realEvent->GetTrack(trackMatchedIndex);
832// //if this is a good track, accept track will return true. The track matched is good, so not track matched is false
833// nottrackmatched = !(fEsdtrackCutsTPC->AcceptTrack(track));
834// }
9ef6f13f 835
836 if(fSecondary){//all particles from secondary interactions
837 written = kTRUE;
838 if(nottrackmatched){//secondaries not removed
839 if (fDepositedCharge != 0){//charged track not removed
840 fChargedNotRemoved++;
841 fEnergyChargedNotRemoved += clEt;
842 fHistRemovedOrNot->Fill(2.0, fCentClass);
843 fHistDxDzNonRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
844 }
845 else{
846 fSecondaryNotRemoved++;
847 }
848 }
849 else{//secondaries removed
850 if (fDepositedCharge != 0){
851 fChargedRemoved++;
852 fEnergyChargedRemoved += clEt;
853 fHistRemovedOrNot->Fill(0.0, fCentClass);
854 fHistDxDzRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
855 fHistChargedTracksCut->Fill(fDepositedEt);
0861cc1f 856 if(fCalcTrackMatchVsMult){
857 fHistChargedTracksCutMult->Fill(fDepositedEt,fClusterMult);
858 if(fClusterMult<25){fHistChargedTracksCutPeripheral->Fill(fDepositedEt);}
859 }
6a152780 860 fHistMatchedTracksEvspTBkgd->Fill(matchedTrackp,fReconstructedE);
0861cc1f 861 if(fCalcTrackMatchVsMult){
6a152780 862 if(fClusterMult<25){fHistMatchedTracksEvspTBkgdPeripheral->Fill(matchedTrackp,fReconstructedEt);}
ac610b08 863 fHistMatchedTracksEvspTBkgdvsCent->Fill(matchedTrackp,fReconstructedEt, fCentClass);
864 fHistMatchedTracksEvspTBkgdvsCentEffCorr->Fill(matchedTrackp,clEt, fCentClass);//Fill with the efficiency corrected energy
0861cc1f 865 }
9ef6f13f 866 Int_t trackindex = (caloCluster->GetLabelsArray())->At(1);
867 if(caloCluster->GetLabel()!=trackindex){
868 fHistBadTrackMatches->Fill(part->Pt(),fReconstructedE);
869 fHistBadTrackMatchesdPhidEta->Fill(caloCluster->GetTrackDx(),caloCluster->GetTrackDz());
870 //cout<<"Track matched, label cluster "<<caloCluster->GetLabel()<<" track "<<trackindex<<endl;
871 }
872 else{
873 fHistGoodTrackMatchesdPhidEta->Fill(caloCluster->GetTrackDx(),caloCluster->GetTrackDz());
874 }
875 }
876 else{//neutral energy removed
877 fNeutralRemoved++;
878 fEnergyNeutralRemoved += clEt;
879 fHistRemovedOrNot->Fill(1.0, fCentClass);
880 fHistDxDzRemovedNeutral->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
881 }
882 }
883 }
884 else{
885
886 if (fDepositedCharge != 0 && fDepositedCode!=fgEMinusCode && fDepositedCode!=fgEPlusCode){//if the particle hitting the calorimeter is pi/k/p/mu
887 written = kTRUE;
d3ce32b8 888 fClusterMultChargedTracks++;
9ef6f13f 889 if(nottrackmatched){//not removed but should be
c79c36cc 890 fHistHadronDepositsAll->Fill(part->Pt());
ac610b08 891 fHistHadronDepositsAllCent->Fill(part->Pt(), fCentClass);
9ef6f13f 892 fChargedNotRemoved++;
893 fEnergyChargedNotRemoved += clEt;
894 fHistRemovedOrNot->Fill(2.0, fCentClass);
895 fHistDxDzNonRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
896 fHistChargedTracksAccepted->Fill(fDepositedEt);
0861cc1f 897 if(fCalcTrackMatchVsMult){
6a152780 898 if(matchedTrackpt<0.5){//if we could never have matched this because of its pt, how much energy did it deposit?
ac610b08 899 fHistChargedTracksAcceptedLowPtCent->Fill(fDepositedEt, fCentClass);
900 if(pdg->PdgCode()!=fgAntiProtonCode){
901 fHistChargedTracksAcceptedLowPtCentNoAntiProtons->Fill(fDepositedEt, fCentClass);
902 }
903 else{
904 fHistChargedTracksAcceptedLowPtCentAntiProtons->Fill(fDepositedEt, fCentClass);
905 }
6a152780 906 }
0861cc1f 907 fHistChargedTracksAcceptedMult->Fill(fDepositedEt,fClusterMult);
908 if(fClusterMult<25){fHistChargedTracksAcceptedPeripheral->Fill(fDepositedEt);}
909 }
9ef6f13f 910 }
911 else{//removed and should have been
912 Int_t trackindex = (caloCluster->GetLabelsArray())->At(0);
c79c36cc 913 fHistHadronDepositsReco->Fill(part->Pt());
ac610b08 914 fHistHadronDepositsRecoCent->Fill(part->Pt(), fCentClass);
c79c36cc 915 fHistHadronDepositsAll->Fill(part->Pt());
ac610b08 916 fHistHadronDepositsAllCent->Fill(part->Pt(), fCentClass);
9ef6f13f 917 if(caloCluster->GetLabel()!=trackindex){
918 fHistBadTrackMatches->Fill(part->Pt(),fReconstructedE);
919 fHistBadTrackMatchesdPhidEta->Fill(caloCluster->GetTrackDx(),caloCluster->GetTrackDz());
920 //cout<<"Track matched, label cluster "<<caloCluster->GetLabel()<<" track "<<trackindex<<endl;
921 }
922 else{
923 fHistGoodTrackMatchesdPhidEta->Fill(caloCluster->GetTrackDx(),caloCluster->GetTrackDz());
924 }
925 fChargedRemoved++;
926 fEnergyChargedRemoved += clEt;
927 fHistRemovedOrNot->Fill(0.0, fCentClass);
928 fHistDxDzRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
929 fHistChargedTracksCut->Fill(fDepositedEt);
0861cc1f 930 if(fCalcTrackMatchVsMult){
931 fHistChargedTracksCutMult->Fill(fDepositedEt,fClusterMult);
932 if(fClusterMult<25){fHistChargedTracksCutPeripheral->Fill(fDepositedEt);}
933 }
6a152780 934 fHistMatchedTracksEvspTBkgd->Fill(matchedTrackp,fReconstructedE);
0861cc1f 935 if(fCalcTrackMatchVsMult){
6a152780 936 if(fClusterMult<25){fHistMatchedTracksEvspTBkgdPeripheral->Fill(matchedTrackp,fReconstructedEt);}
ac610b08 937 fHistMatchedTracksEvspTBkgdvsCent->Fill(matchedTrackp,fReconstructedEt, fCentClass);
938 fHistMatchedTracksEvspTBkgdvsCentEffCorr->Fill(matchedTrackp,clEt, fCentClass);//fill with the efficiency corrected energy
0861cc1f 939 }
9ef6f13f 940 }
941 }
942 //K0L and any neutral particles from the decay of K+/- or K0S
943 if(!written && (fPrimaryCode==fgKPlusCode || fPrimaryCode==fgKMinusCode || fPrimaryCode==fgK0SCode ||fPrimaryCode==fgK0LCode)){
944 written = kTRUE;//At this point we are not tracking them but we don't count them as neutrals accidentally removed
cfe23ff0 945 }
9ef6f13f 946
947 if(!written && (fDepositedCode==fgGammaCode || fDepositedCode==fgEMinusCode || fDepositedCode ==fgEPlusCode)){//if the particle hitting the calorimeter is gamma, electron and not from a kaon
d3ce32b8 948 fClusterMultGammas++;
9ef6f13f 949 written = kTRUE;
950 if(nottrackmatched){//Not removed and not supposed to be removed - signal
951 fEtNonRemovedGammas += clEt;
952 fMultNonRemovedGammas++;
953 fNeutralNotRemoved--;
954 fEnergyNeutralNotRemoved -= clEt;
955 fHistGammasAccepted->Fill(fDepositedEt);
0861cc1f 956 if(fCalcTrackMatchVsMult){
957 fHistGammasAcceptedMult->Fill(fDepositedEt,fClusterMult);
958 if(fClusterMult<25){fHistGammasAcceptedPeripheral->Fill(fDepositedEt);}
959 }
9ef6f13f 960 }
961 else{//removed but shouldn't have been
962 Int_t trackindex = (caloCluster->GetLabelsArray())->At(1);
963 if(caloCluster->GetLabel()!=trackindex){
964 fHistBadTrackMatches->Fill(part->Pt(),fReconstructedE);
965 fHistBadTrackMatchesdPhidEta->Fill(caloCluster->GetTrackDx(),caloCluster->GetTrackDz());
966// cout<<"Track matched, label cluster "<<caloCluster->GetLabel()<<" track "<<trackindex<<endl;
967// PrintFamilyTree(trackindex, stack);
968// cout<<"Cluster"<<endl;
969 }
970 fGammaRemoved++;
971 fGammaRemovedEt+=clEt;
972 fHistGammasCut->Fill(fDepositedEt);
0861cc1f 973 if(fCalcTrackMatchVsMult){
974 fHistGammasCutMult->Fill(fDepositedEt,fClusterMult);
975 if(fClusterMult<25){fHistGammasCutPeripheral->Fill(fDepositedEt);}
976 }
6a152780 977 fHistMatchedTracksEvspTSignal->Fill(matchedTrackp,fReconstructedE);
0861cc1f 978 if(fCalcTrackMatchVsMult){
6a152780 979 if(fClusterMult<25){fHistMatchedTracksEvspTSignalPeripheral->Fill(matchedTrackp,fReconstructedEt);}
ac610b08 980 fHistMatchedTracksEvspTSignalvsCent->Fill(matchedTrackp,fReconstructedEt, fCentClass);
981 fHistMatchedTracksEvspTSignalvsCentEffCorr->Fill(matchedTrackp,clEt, fCentClass);
0861cc1f 982 }
9ef6f13f 983 }
984 }
985 //all other cases - neutron, anti-neutron, not aware of other cases
986 if(!written){
987 fNeutralNotRemoved++;
988 fEnergyNeutralNotRemoved += clEt;
989 fHistRemovedOrNot->Fill(3.0, fCentClass);
990 fHistDxDzNonRemovedNeutral->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
991 }
992 }
f61cec2f 993 fPrimaryTree->Fill();
9ef6f13f 994 } // end of loop over clusters
e2ee5727 995
b2c10007 996
997 std::sort(foundGammas.begin(), foundGammas.end());
998 for (Int_t iPart = 0; iPart < stack->GetNtrack(); iPart++)
999 {
1000
1001 if(!stack->IsPhysicalPrimary(iPart)) continue;
1002
1003 TParticle *part = stack->Particle(iPart);
1004
1005 if (!part)
1006 {
1007 Printf("ERROR: Could not get particle %d", iPart);
1008 continue;
1009 }
1010 TParticlePDG *pdg = part->GetPDG(0);
1011
1012 if (!pdg)
1013 {
1014 Printf("ERROR: Could not get particle PDG %d", iPart);
1015 continue;
1016 }
1017
1018 if(pdg->PdgCode()==fgGammaCode && fSelector->CutGeometricalAcceptance(*part))// TMath::Abs(part->Eta()) < 0.12)
1019 {
1020 fHistGammasGenerated->Fill(part->Energy());
6a152780 1021 fHistGammasGeneratedMult->Fill(part->Energy(),fClusterMult);
b2c10007 1022 if(std::binary_search(foundGammas.begin(),foundGammas.end(),iPart))
1023 {
1024 fHistGammasFound->Fill(part->Energy());
6a152780 1025 fHistGammasFoundMult->Fill(part->Energy(),fClusterMult);
b2c10007 1026 }
b2c10007 1027 }
6a152780 1028 if(pdg->PdgCode()==fgPiPlusCode || pdg->PdgCode()==fgPiMinusCode || pdg->PdgCode()==fgProtonCode || pdg->PdgCode()==fgAntiProtonCode){//section here for all hadrons generated
ac610b08 1029 fHistHadronsAllCent->Fill(part->Pt(), fCentClass);
6a152780 1030 }
1031
b2c10007 1032
1033 }
0861cc1f 1034 if(fCalcForKaonCorrection){
1035 Float_t etCuts[11] = {0.0,0.05,0.1,0.15,0.2,0.25, 0.3,0.35,0.4,0.45,0.5};
1036 Int_t nEtCuts = 11;
1037 //loop over simulated particles in order to find K0S
1038 for (Int_t iPart = 0; iPart < stack->GetNtrack(); iPart++){
1039 TParticle *part = stack->Particle(iPart);
1040 if (!part){
1041 //Printf("ERROR: Could not get particle %d", iPart);
1042 continue;
1043 }
1044 TParticlePDG *pdg = part->GetPDG(0);
1045 if (!pdg){
1046 //Printf("ERROR: Could not get particle PDG %d", iPart);
1047 continue;
1048 }
1049 //if(stack->IsPhysicalPrimary(iPart)){//if it is a K0 it might have decayed into four pions
2aab9269 1050 //fgK0SCode, fgGammaCode, fgPi0Code
1051 Int_t code = pdg->PdgCode();
0861cc1f 1052 if(code == fgK0SCode || code==fgKPlusCode || code==fgKMinusCode ||code==fgK0LCode || code==fgK0Code){//this is a kaon
1053 //cout<<"I am a kaon too! "<<stack->Particle(iPart)->GetName()<<" "<<code<<endl;
1054 Float_t pTk = stack->Particle(iPart)->Pt();
1055 if(TMath::Abs(stack->Particle(iPart)->Y())<0.5 && stack->IsPhysicalPrimary(iPart)){//these are particles which would be included in our spectra measurements
1056 fHistSimKaonsInAcceptance->Fill(pTk);
1057 if(code == fgK0SCode){fHistSimK0SInAcceptance->Fill(pTk);}
1058 if(code == fgK0LCode){fHistSimK0LInAcceptance->Fill(pTk);}
1059 if(code == fgKPlusCode){fHistSimKPlusInAcceptance->Fill(pTk);}
1060 if(code == fgKMinusCode){fHistSimKMinusInAcceptance->Fill(pTk);}
1061 if(code == fgK0Code){//Split K0's between the two
1062 fHistSimK0SInAcceptance->Fill(pTk,0.5);
1063 fHistSimK0LInAcceptance->Fill(pTk,0.5);
1064 }
1065 }
1066 else{
1067 fHistSimKaonsOutOfAcceptance->Fill(pTk);
1068 // if(!stack->IsPhysicalPrimary(iPart)){
1069 // PrintFamilyTree(iPart, stack);
1070 // }
1071 }
1072 Float_t totalGammaEts[11] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0};
1073 Float_t totalClusterEts[11] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0};
1074 for (int iCluster = 0; iCluster < nCluster; iCluster++ ){//if this cluster is from any of the decay daughters of any kaon... but there is no easy way to look at this so we loop over clusters...
1075 AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
1076 const UInt_t myPart = (UInt_t)TMath::Abs(caloCluster->GetLabel());
1077 //identify the primary particle which created this cluster
1078 int primIdx = myPart;
1079 if (!stack->IsPhysicalPrimary(myPart)){
1080 primIdx = GetPrimMother(iPart, stack);
1081 } // end of primary particle check
29b7ba7f 1082 TParticle *hitPart = stack->Particle(myPart);
1083 Bool_t hitsAsChargedKaon = kFALSE;
1084 if(hitPart->GetPdgCode()== fgKPlusCode || hitPart->GetPdgCode()== fgKPlusCode){
1085 if(myPart==primIdx){
1086 //The particle hits as a charged kaon and that kaon is a primary kaon - do not count because this is counted in the hadronic correction!
1087 hitsAsChargedKaon = kTRUE;
1088 //cout<<"Found primary charged kaon cluster!"<<endl;
1089 }
1090 }
1091 if(primIdx==iPart && primIdx>0 && !hitsAsChargedKaon){//This cluster is from our primary particle and our primary particle is a kaon
0861cc1f 1092 //cout<<"I have a particle match! prim code"<<code<<" id "<<primIdx <<endl;
1093 Float_t pos[3];
1094 caloCluster->GetPosition(pos);
1095 TVector3 cp(pos);
1096 Double_t clEt = caloCluster->E()*TMath::Sin(cp.Theta());
47151f26 1097 Double_t clEtCorr = CorrectForReconstructionEfficiency(*caloCluster,fClusterMult);
2aab9269 1098 for(int l=0;l<nEtCuts;l++){//loop over cut values
0861cc1f 1099 if(clEt>=etCuts[l]){
1100 //cout<<", "<<clEt<<">="<<etCuts[l];
1101 totalClusterEts[l] += clEtCorr;//if cluster et is above the cut off energy add it
1102 totalGammaEts[l] += clEt;//if cluster et is above the cut off energy add it
2aab9269 1103 }
1104 }
2aab9269 1105 }
1106 }
0861cc1f 1107 // cout<<"Deposits: pT: "<<pTk;
1108 // for(int l=0;l<nEtCuts;l++){//loop over cut values
1109 // cout<<" "<<totalClusterEts[l];
1110 // }
1111 // cout<<endl;
1112 if(TMath::Abs(stack->Particle(iPart)->Y())<0.5 && stack->IsPhysicalPrimary(iPart)){//within the acceptance of our spectra and is a primary particle
1113 if(totalClusterEts[0]>0.0){fHistSimKaonsInAcceptanceWithDepositsPrimaries->Fill(pTk);}
1114 //cout<<"I have a particle match! prim code"<<code<<" id "<<iPart <<endl;
1115 for(int l=0;l<nEtCuts;l++){
d3ce32b8 1116 fHistK0EDepositsVsPtInAcceptance->Fill(pTk,totalClusterEts[l],etCuts[l]+0.001);
1117 fHistK0EGammaVsPtInAcceptance->Fill(pTk,totalGammaEts[l],etCuts[l]+0.001);
0861cc1f 1118 }
1119 }
1120 else{//outside the acceptance of our spectra
1121 if(totalClusterEts[0]>0.0){
1122 if(stack->IsPhysicalPrimary(iPart)){fHistSimKaonsOutOfAcceptanceWithDepositsPrimaries->Fill(pTk);}
1123 else{fHistSimKaonsOutOfAcceptanceWithDepositsSecondaries->Fill(pTk);}
1124 }
1125 for(int l=0;l<nEtCuts;l++){
d3ce32b8 1126 fHistK0EDepositsVsPtOutOfAcceptance->Fill(pTk,totalClusterEts[l],etCuts[l]+0.001);
1127 fHistK0EGammaVsPtOutOfAcceptance->Fill(pTk,totalGammaEts[l],etCuts[l]+0.001);
0861cc1f 1128 }
1129 }
1130
2aab9269 1131 }
1132 }
1133 }
d3ce32b8 1134 fHistMultChVsSignalVsMult->Fill(fClusterMultChargedTracks,fClusterMultGammas,fClusterMult);
e2ee5727 1135 FillHistograms();
1136 return 0;
1137}
2fbf38ac 1138
1139void AliAnalysisEtMonteCarlo::Init()
f61cec2f 1140{ // init
e2ee5727 1141 AliAnalysisEt::Init();
2fbf38ac 1142}
13b0d3c1 1143
ce546038 1144void AliAnalysisEtMonteCarlo::ResetEventValues()
f61cec2f 1145{ // reset event values
e2ee5727 1146 AliAnalysisEt::ResetEventValues();
c79c36cc 1147 if(!fIsMC) return;
e2ee5727 1148
f61cec2f 1149 fTotEtSecondary = 0;
1150 fTotEtSecondaryFromEmEtPrimary = 0;
1151 fTotEtWithSecondaryRemoved = 0;
1152
e2ee5727 1153 // collision geometry defaults for p+p:
1154 fImpactParameter = 0;
1155 fNcoll = 1;
1156 fNpart = 2;
1157
1158 fEtNonRemovedProtons = 0;
1159 fEtNonRemovedAntiProtons = 0;
1160 fEtNonRemovedPiPlus = 0;
1161 fEtNonRemovedPiMinus = 0;
1162 fEtNonRemovedKaonPlus = 0;
1163 fEtNonRemovedKaonMinus = 0;
ef647350 1164 fEtNonRemovedK0S = 0;
1165 fEtNonRemovedK0L = 0;
e2ee5727 1166 fEtNonRemovedLambdas = 0;
1167 fEtNonRemovedElectrons = 0;
1168 fEtNonRemovedPositrons = 0;
1169 fEtNonRemovedMuPlus = 0;
1170 fEtNonRemovedMuMinus = 0;
1171 fEtNonRemovedNeutrons = 0;
1172 fEtNonRemovedAntiNeutrons = 0;
1173 fEtNonRemovedGammas = 0;
1174 fEtNonRemovedGammasFromPi0 = 0;
1175
ef647350 1176 fEtRemovedProtons = 0;
1177 fEtRemovedAntiProtons = 0;
1178 fEtRemovedPiPlus = 0;
1179 fEtRemovedPiMinus = 0;
1180 fEtRemovedKaonPlus = 0;
1181 fEtRemovedKaonMinus = 0;
1182 fEtRemovedK0s = 0;
1183 fEtRemovedK0L = 0;
1184 fEtRemovedLambdas = 0;
1185 fEtRemovedElectrons = 0;
1186 fEtRemovedPositrons = 0;
1187 fEtRemovedMuPlus = 0;
1188 fEtRemovedMuMinus = 0;
1189 fEtRemovedNeutrons = 0;
1190
1191 fEtRemovedGammasFromPi0 = 0;
e2ee5727 1192 fEtRemovedGammas = 0;
1193 fEtRemovedNeutrons = 0;
1194 fEtRemovedAntiNeutrons = 0;
1195
1196 fMultNonRemovedProtons = 0;
1197 fMultNonRemovedAntiProtons = 0;
1198 fMultNonRemovedPiPlus = 0;
1199 fMultNonRemovedPiMinus = 0;
1200 fMultNonRemovedKaonPlus = 0;
1201 fMultNonRemovedKaonMinus = 0;
1202 fMultNonRemovedK0s = 0;
ef647350 1203 fMultNonRemovedK0L = 0;
e2ee5727 1204 fMultNonRemovedLambdas = 0;
1205 fMultNonRemovedElectrons = 0;
1206 fMultNonRemovedPositrons = 0;
1207 fMultNonRemovedMuPlus = 0;
1208 fMultNonRemovedMuMinus = 0;
1209 fMultNonRemovedNeutrons = 0;
1210 fMultNonRemovedAntiNeutrons = 0;
1211 fMultNonRemovedGammas = 0;
1212
ef647350 1213 fMultRemovedProtons = 0;
1214 fMultRemovedAntiProtons = 0;
1215 fMultRemovedPiPlus = 0;
1216 fMultRemovedPiMinus = 0;
1217 fMultRemovedKaonPlus = 0;
1218 fMultRemovedKaonMinus = 0;
1219 fMultRemovedK0s = 0;
1220 fMultRemovedK0L = 0;
1221 fMultRemovedLambdas = 0;
1222 fMultRemovedElectrons = 0;
1223 fMultRemovedPositrons = 0;
1224 fMultRemovedMuPlus = 0;
1225 fMultRemovedMuMinus = 0;
1226
e2ee5727 1227 fMultRemovedGammas = 0;
1228 fMultRemovedNeutrons = 0;
1229 fMultRemovedAntiNeutrons = 0;
1230
ef647350 1231 fEnergyChargedNotRemoved = 0;
1232 fEnergyChargedRemoved = 0;
1233 fEnergyNeutralNotRemoved = 0;
1234 fEnergyNeutralRemoved = 0;
f61cec2f 1235
ef647350 1236 fChargedNotRemoved = 0;
1237 fChargedRemoved = 0;
1238 fNeutralNotRemoved = 0;
1239 fNeutralRemoved = 0;
b2c10007 1240 fGammaRemoved = 0;
1241 fSecondaryNotRemoved = 0;
f61cec2f 1242
e2ee5727 1243 fTrackMultInAcc = 0;
1244
f61cec2f 1245 fTotNeutralEtAfterMinEnergyCut = 0;
b2c10007 1246
1247 fSecondaryNotRemoved = 0;
1248
1249 fTotPx = 0;
1250 fTotPy = 0;
1251 fTotPz = 0;
1252
1253
ce546038 1254}
1255
1256void AliAnalysisEtMonteCarlo::CreateHistograms()
f61cec2f 1257{ // histogram related additions
e2ee5727 1258 AliAnalysisEt::CreateHistograms();
c79c36cc 1259 if(!fIsMC) return;
f61cec2f 1260 if (fEventSummaryTree) {
1261 fEventSummaryTree->Branch("fImpactParameter",&fImpactParameter,"fImpactParameter/D");
1262 fEventSummaryTree->Branch("fNcoll",&fNcoll,"fNcoll/I");
1263 fEventSummaryTree->Branch("fNpart",&fNpart,"fNpart/I");
1264 fEventSummaryTree->Branch("fTotEtWithSecondaryRemoved", &fTotEtWithSecondaryRemoved, "fTotEtWithSecondaryRemoved/D");
1265 fEventSummaryTree->Branch("fTotEtSecondaryFromEmEtPrimary", &fTotEtSecondaryFromEmEtPrimary, "fTotEtSecondaryFromEmEtPrimary/D");
1266 fEventSummaryTree->Branch("fTotEtSecondary", &fTotEtSecondary, "fTotEtSecondary/D");
1267 fEventSummaryTree->Branch("fTotNeutralEtAfterMinEnergyCut", &fTotNeutralEtAfterMinEnergyCut, "fTotNeutralEtAfterMinEnergyCut/D");
b2c10007 1268 fEventSummaryTree->Branch("fSecondaryNotRemoved", &fSecondaryNotRemoved, "fSecondaryNotRemoved/I");
1269 fEventSummaryTree->Branch("fChargedNotRemoved", &fChargedNotRemoved, "fChargedNotRemoved/I");
1270 fEventSummaryTree->Branch("fNeutralNotRemoved", &fNeutralNotRemoved, "fNeutralNotRemoved/I");
1271 fEventSummaryTree->Branch("fChargedRemoved", &fChargedRemoved, "fChargedRemoved/I");
1272 fEventSummaryTree->Branch("fNeutralRemoved", &fNeutralRemoved, "fNeutralRemoved/I");
1273 fEventSummaryTree->Branch("fGammaRemoved", &fGammaRemoved, "fGammaRemoved/I");
1274 fEventSummaryTree->Branch("fTotPx", &fTotPx, "fTotPx/D");
1275 fEventSummaryTree->Branch("fTotPy", &fTotPy, "fTotPy/D");
1276 fEventSummaryTree->Branch("fTotPz", &fTotPz, "fTotPz/D");
1277// fEventSummaryTree->Branch("f
e2ee5727 1278 }
f61cec2f 1279
e2ee5727 1280 //fHistDecayVertexNonRemovedCharged = new TH3F("fHistDecayVertexNonRemovedCharged","fHistDecayVertexNonRemovedCharged", 500, -470, 30, 500, -300, 300, 40, -20, 20);
1281 //fHistDecayVertexRemovedCharged = new TH3F("fHistDecayVertexRemovedCharged","fHistDecayVertexRemovedCharged", 500, -470, 30, 500, -300, 300, 40, -20, 20);
1282 //fHistDecayVertexNonRemovedNeutral = new TH3F("fHistDecayVertexNonRemovedNeutral","fHistDecayVertexNonRemovedNeutral", 500, -470, 30, 500, -300, 300, 40, -20, 20);
1283 //fHistDecayVertexRemovedNeutral = new TH3F("fHistDecayVertexRemovedNeutral","fHistDecayVertexRemovedNeutral", 500, -470, 30, 500, -300, 300, 40, -20, 20);
1284
1285 fHistRemovedOrNot = new TH2F("fHistRemovedOrNot", "fHistRemovedOrNot", 4, -0.5, 3.5, 10, -0.5, 9.5);
1286
1287 fHistEtNonRemovedProtons = new TH2F("fHistEtNonRemovedProtons", "fHistEtNonRemovedProtons", 1500, 0, 30, 10, -0.5, 9.5);
1288 fHistEtNonRemovedAntiProtons = new TH2F("fHistEtNonRemovedAntiProtons", "fHistEtNonRemovedAntiProtons", 1500, 0, 30, 10, -0.5, 9.5);
1289 fHistEtNonRemovedPiPlus = new TH2F("fHistEtNonRemovedPiPlus", "fHistEtNonRemovedPiPlus", 1500, 0, 30, 10, -0.5, 9.5);
1290 fHistEtNonRemovedPiMinus = new TH2F("fHistEtNonRemovedPiMinus", "fHistEtNonRemovedPiMinus", 1500, 0, 30, 10, -0.5, 9.5);
1291 fHistEtNonRemovedKaonPlus = new TH2F("fHistEtNonRemovedKaonPlus", "fHistEtNonRemovedKaonPlus", 1500, 0, 30, 10, -0.5, 9.5);
1292 fHistEtNonRemovedKaonMinus = new TH2F("fHistEtNonRemovedKaonMinus", "fHistEtNonRemovedKaonMinus", 1500, 0, 30, 10, -0.5, 9.5);
1293 fHistEtNonRemovedK0s = new TH2F("fHistEtNonRemovedK0s", "fHistEtNonRemovedK0s", 1500, 0, 30, 10, -0.5, 9.5);
ef647350 1294 fHistEtNonRemovedK0L = new TH2F("fHistEtNonRemovedK0L", "fHistEtNonRemovedK0L", 1500, 0, 30, 10, -0.5, 9.5);
e2ee5727 1295 fHistEtNonRemovedLambdas = new TH2F("fHistEtNonRemovedLambdas", "fHistEtNonRemovedLambdas", 1500, 0, 30, 10, -0.5, 9.5);
1296 fHistEtNonRemovedElectrons = new TH2F("fHistEtNonRemovedElectrons", "fHistEtNonRemovedElectrons", 1500, 0, 30, 10, -0.5, 9.5);
1297 fHistEtNonRemovedPositrons = new TH2F("fHistEtNonRemovedPositrons", "fHistEtNonRemovedPositrons", 1500, 0, 30, 10, -0.5, 9.5);
1298 fHistEtNonRemovedMuPlus = new TH2F("fHistEtNonRemovedMuPlus", "fHistEtNonRemovedMuPlus", 1500, 0, 30, 10, -0.5, 9.5);
1299 fHistEtNonRemovedMuMinus = new TH2F("fHistEtNonRemovedMuMinus", "fHistEtNonRemovedMuMinus", 1500, 0, 30, 10, -0.5, 9.5);
1300 fHistEtNonRemovedNeutrons = new TH2F("fHistEtNonRemovedNeutrons", "fHistEtNonRemovedNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
1301 fHistEtNonRemovedAntiNeutrons = new TH2F("fHistEtNonRemovedAntiNeutrons", "fHistEtNonRemovedAntiNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
1302
1303 fHistEtNonRemovedGammas = new TH2F("fHistEtNonRemovedGammas", "fHistEtNonRemovedGammas", 1500, 0, 30, 10, -0.5, 9.5);
1304 fHistEtNonRemovedGammasFromPi0 = new TH2F("fHistEtNonRemovedGammasFromPi0", "fHistEtNonRemovedGammasFromPi0", 1500, 0, 30, 10, -0.5, 9.5);
1305
1306 fHistEtRemovedGammas = new TH2F("fHistEtRemovedGammas", "fHistEtRemovedGammas", 1500, 0, 30, 10, -0.5, 9.5);
1307 fHistEtRemovedNeutrons = new TH2F("fHistEtRemovedNeutrons", "fHistEtRemovedNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
1308 fHistEtRemovedAntiNeutrons = new TH2F("fHistEtRemovedAntiNeutrons", "fHistEtRemovedAntiNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
f61cec2f 1309
ef647350 1310 fHistEtRemovedCharged = new TH2F("fHistEtRemovedCharged", "fHistEtRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
1311 fHistEtRemovedNeutrals = new TH2F("fHistEtRemovedNeutrals", "fHistEtRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
e2ee5727 1312
ef647350 1313 fHistEtNonRemovedCharged = new TH2F("fHistEtNonRemovedCharged", "fHistEtNonRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
1314 fHistEtNonRemovedNeutrals = new TH2F("fHistEtNonRemovedNeutrals", "fHistEtNonRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
f61cec2f 1315
e2ee5727 1316 fHistMultNonRemovedProtons = new TH2F("fHistMultNonRemovedProtons", "fHistMultNonRemovedProtons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1317 fHistMultNonRemovedAntiProtons = new TH2F("fHistMultNonRemovedAntiProtons", "fHistMultNonRemovedAntiProtons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1318 fHistMultNonRemovedPiPlus = new TH2F("fHistMultNonRemovedPiPlus", "fHistMultNonRemovedPiPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1319 fHistMultNonRemovedPiMinus = new TH2F("fHistMultNonRemovedPiMinus", "fHistMultNonRemovedPiMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1320 fHistMultNonRemovedKaonPlus = new TH2F("fHistMultNonRemovedKaonPlus", "fHistMultNonRemovedKaonPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1321 fHistMultNonRemovedKaonMinus = new TH2F("fHistMultNonRemovedKaonMinus", "fHistMultNonRemovedKaonMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1322 fHistMultNonRemovedK0s = new TH2F("fHistMultNonRemovedK0s", "fHistMultNonRemovedK0s", 100, -0.5, 99.5, 10, -0.5, 9.5);
ef647350 1323 fHistMultNonRemovedK0L = new TH2F("fHistMultNonRemovedK0L", "fHistMultNonRemovedK0L", 100, -0.5, 99.5, 10, -0.5, 9.5);
e2ee5727 1324 fHistMultNonRemovedLambdas = new TH2F("fHistMultNonRemovedLambdas", "fHistMultNonRemovedLambdas", 100, -0.5, 99.5, 10, -0.5, 9.5);
1325 fHistMultNonRemovedElectrons = new TH2F("fHistMultNonRemovedElectrons", "fHistMultNonRemovedElectrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1326 fHistMultNonRemovedPositrons = new TH2F("fHistMultNonRemovedPositrons", "fHistMultNonRemovedPositrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1327 fHistMultNonRemovedMuPlus = new TH2F("fHistMultNonRemovedMuPlus", "fHistMultNonRemovedMuPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1328 fHistMultNonRemovedMuMinus = new TH2F("fHistMultNonRemovedMuMinus", "fHistMultNonRemovedMuMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1329 fHistMultNonRemovedNeutrons = new TH2F("fHistMultNonRemovedNeutrons", "fHistMultNonRemovedNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1330 fHistMultNonRemovedAntiNeutrons = new TH2F("fHistMultNonRemovedAntiNeutrons", "fHistMultNonRemovedAntiNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1331
ef647350 1332 fHistMultNonRemovedGammas = new TH2F("fHistMultNonRemovedGammas", "fHistMultNonRemovedGammas", 100, -0.5, 99.5, 100, -0.5, 99.5);
e2ee5727 1333
ef647350 1334 fHistMultRemovedGammas = new TH2F("fHistMultRemovedGammas", "fHistMultRemovedGammas", 100, -0.5, 99.5, 100, -0.5, 99.5);
e2ee5727 1335 fHistMultRemovedNeutrons = new TH2F("fHistMultRemovedNeutrons", "fHistMultRemovedNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1336 fHistMultRemovedAntiNeutrons = new TH2F("fHistMultRemovedAntiNeutrons", "fHistMultRemovedAntiNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
f61cec2f 1337 /*
1338 fHistMultRemovedCharged = new TH2F("fHistMultRemovedCharged", "fHistMultRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
1339 fHistMultRemovedNeutrals = new TH2F("fHistMultRemovedNeutrals", "fHistMultRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
1340
1341 fHistMultNonRemovedCharged = new TH2F("fHistMultNonRemovedCharged", "fHistMultNonRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
1342 fHistMultNonRemovedNeutrals = new TH2F("fHistMultNonRemovedNeutrals", "fHistMultNonRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);*/
ef647350 1343
ef647350 1344
ef647350 1345 fHistMultRemovedCharged = new TH2F("fHistMultRemovedCharged", "fHistMultRemovedCharged", 100, -0.5, 99.5, 100, -0.5, 99.5);
1346 fHistMultRemovedNeutrals = new TH2F("fHistMultRemovedNeutrals", "fHistMultRemovedNeutrals", 100, -0.5, 99.5, 100, -0.5, 99.5);
1347
1348 fHistMultNonRemovedCharged = new TH2F("fHistMultNonRemovedCharged", "fHistMultNonRemovedCharged", 100, -0.5, 99.5, 100, -0.5, 99.5);
1349 fHistMultNonRemovedNeutrals = new TH2F("fHistMultNonRemovedNeutrals", "fHistMultNonRemovedNeutrals", 100, -0.5, 99.5, 100, -0.5, 99.5);
e2ee5727 1350
1351 fHistTrackMultvsNonRemovedCharged = new TH2F("fHistTrackMultvsNonRemovedCharged", "fHistTrackMultvsNonRemovedCharged", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1352 fHistTrackMultvsNonRemovedNeutral = new TH2F("fHistTrackMultvsNonRemovedNeutral", "fHistTrackMultvsNonRemovedNeutral", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1353 fHistTrackMultvsRemovedGamma = new TH2F("fHistTrackMultvsRemovedGamma", "fHistTrackMultvsRemovedGamma", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1354
1355 fHistClusterMultvsNonRemovedCharged = new TH2F("fHistClusterMultvsNonRemovedCharged", "fHistClusterMultvsNonRemovedCharged", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1356 fHistClusterMultvsNonRemovedNeutral = new TH2F("fHistClusterMultvsNonRemovedNeutral", "fHistClusterMultvsNonRemovedNeutral", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1357 fHistClusterMultvsRemovedGamma = new TH2F("fHistClusterMultvsRemovedGamma", "fHistClusterMultvsRemovedGamma", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1358
1359 fHistDxDzNonRemovedCharged = new TH2F("fHistDxDzNonRemovedCharged", "fHistDxDzNonRemovedCharged", 800, -200, 200, 800, -200, 200);
1360 fHistDxDzRemovedCharged = new TH2F("fHistDxDzRemovedCharged", "fHistDxDzRemovedCharged", 800, -200, 200, 800, -200, 200);
1361 fHistDxDzNonRemovedNeutral = new TH2F("fHistDxDzNonRemovedNeutral", "fHistDxDzNonRemovedNeutral", 800, -200, 200, 800, -200, 200);
1362 fHistDxDzRemovedNeutral = new TH2F("fHistDxDzRemovedNeutral", "fHistDxDzRemovedNeutral", 800, -200, 200, 800, -200, 200);
1363
0861cc1f 1364 if(fCalcForKaonCorrection){
1365 Int_t nEtCut = 11;
1366 Float_t etCutAxis[12] = {0.00,0.05,0.10,0.15,0.20, 0.25,0.30,0.35,0.40,0.45, 0.50,0.55};
1367 fHistK0EDepositsVsPtInAcceptance = new TH3F("fHistK0EDepositsVsPtInAcceptance","Kaon deposits with corrections for kaons with y<0.5",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis,nEtCut,etCutAxis);
1368 fHistK0EGammaVsPtInAcceptance = new TH3F("fHistK0EGammaVsPtInAcceptance","Kaon deposits without corrections for kaons with y<0.5",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis,nEtCut,etCutAxis);
1369 fHistK0EDepositsVsPtOutOfAcceptance = new TH3F("fHistK0EDepositsVsPtOutOfAcceptance","Kaon deposits with corrections for kaons with y>0.5",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis,nEtCut,etCutAxis);
1370 fHistK0EGammaVsPtOutOfAcceptance = new TH3F("fHistK0EGammaVsPtOutOfAcceptance","Kaon deposits without corrections for kaons with y>0.5",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis,nEtCut,etCutAxis);
1371 fHistK0EDepositsVsPtInAcceptance->GetXaxis()->SetTitle("Kaon p_{T}");
1372 fHistK0EGammaVsPtInAcceptance->GetXaxis()->SetTitle("Kaon p_{T}");
1373 fHistK0EDepositsVsPtOutOfAcceptance->GetXaxis()->SetTitle("Kaon p_{T}");
1374 fHistK0EGammaVsPtOutOfAcceptance->GetXaxis()->SetTitle("Kaon p_{T}");
1375 fHistK0EDepositsVsPtInAcceptance->GetYaxis()->SetTitle("Deposited E_{T}");
1376 fHistK0EGammaVsPtInAcceptance->GetYaxis()->SetTitle("Deposited E_{T}");
1377 fHistK0EDepositsVsPtOutOfAcceptance->GetYaxis()->SetTitle("Deposited E_{T}");
1378 fHistK0EGammaVsPtOutOfAcceptance->GetYaxis()->SetTitle("Deposited E_{T}");
1379 fHistK0EDepositsVsPtInAcceptance->GetZaxis()->SetTitle("E_{T} cut");
1380 fHistK0EGammaVsPtInAcceptance->GetZaxis()->SetTitle("E_{T} cut");
1381 fHistK0EDepositsVsPtOutOfAcceptance->GetZaxis()->SetTitle("E_{T} cut");
1382 fHistK0EGammaVsPtOutOfAcceptance->GetZaxis()->SetTitle("E_{T} cut");
1383
1384 fHistSimKaonsInAcceptance = new TH1F("fHistSimKaonsInAcceptance","Kaons with y<0.5",fgNumOfPtBins,fgPtAxis);
1385 fHistSimK0SInAcceptance = new TH1F("fHistSimK0SInAcceptance","Kaons with y<0.5",fgNumOfPtBins,fgPtAxis);
1386 fHistSimK0LInAcceptance = new TH1F("fHistSimK0LInAcceptance","Kaons with y<0.5",fgNumOfPtBins,fgPtAxis);
1387 fHistSimKPlusInAcceptance = new TH1F("fHistSimKPlusInAcceptance","Kaons with y<0.5",fgNumOfPtBins,fgPtAxis);
1388 fHistSimKMinusInAcceptance = new TH1F("fHistSimKMinusInAcceptance","Kaons with y<0.5",fgNumOfPtBins,fgPtAxis);
1389 fHistSimKaonsOutOfAcceptance = new TH1F("fHistSimKaonsOutOfAcceptance","Kaons with y>0.5",fgNumOfPtBins,fgPtAxis);
1390 fHistSimKaonsInAcceptanceWithDepositsPrimaries = new TH1F("fHistSimKaonsInAcceptanceWithDepositsPrimaries","Primary Kaons which deposited energy in calorimeter with y>0.5",fgNumOfPtBins,fgPtAxis);
1391 fHistSimKaonsOutOfAcceptanceWithDepositsSecondaries = new TH1F("fHistSimKaonsOutOfAcceptanceWithDepositsSecondaries","Secondary Kaons which deposited energy in calorimeter with y>0.5",fgNumOfPtBins,fgPtAxis);
1392 fHistSimKaonsOutOfAcceptanceWithDepositsPrimaries = new TH1F("fHistSimKaonsOutOfAcceptanceWithDepositsPrimaries","Primary Kaons which deposited energy in calorimeter with y>0.5",fgNumOfPtBins,fgPtAxis);
1393 }
2aab9269 1394
e2ee5727 1395 fHistPiPlusMult = new TH1F("fHistPiPlusMult", "fHistPiPlusMult", 2000, -0.5, 1999.5);
1396 fHistPiMinusMult = new TH1F("fHistPiMinusMult", "fHistPiMinusMult", 2000, -0.5, 1999.5);
1397 fHistPiZeroMult = new TH1F("fHistPiZeroMult", "fHistPiZeroMult", 2000, -0.5, 1999.5);
1398
1399 fHistPiPlusMultAcc = new TH1F("fHistPiPlusMultAcc", "fHistPiPlusMultAcc", 2000, -0.5, 1999.5);
1400 fHistPiMinusMultAcc = new TH1F("fHistPiMinusMultAcc", "fHistPiMinusMultAcc", 2000, -0.5, 1999.5);
1401 fHistPiZeroMultAcc = new TH1F("fHistPiZeroMultAcc", "fHistPiZeroMultAcc", 2000, -0.5, 1999.5);
1402
f61cec2f 1403 if(fCuts->GetHistMakeTree())
1404 {
1405 TString treename = "fPrimaryTree" + fHistogramNameSuffix;
1406 fPrimaryTree = new TTree(treename, treename);
1407
1408 fPrimaryTree->Branch("fTotEt",&fTotEt,"fTotEt/D");
1409 fPrimaryTree->Branch("fNeutralMultiplicity", &fNeutralMultiplicity, "fNeutralMultiplicity/I");
1410 fPrimaryTree->Branch("fCentClass",&fCentClass,"fCentClass/I");
1411
1412 fPrimaryTree->Branch("fPrimaryCode", &fPrimaryCode, "fPrimaryCode/I");
1413 fPrimaryTree->Branch("fPrimaryCharge", &fPrimaryCharge, "fPrimaryCharge/I");
1414
1415 fPrimaryTree->Branch("fPrimaryE", &fPrimaryE, "fPrimaryE/D");
1416 fPrimaryTree->Branch("fPrimaryEt", &fPrimaryEt, "fPrimaryEt/D");
1417
1418 fPrimaryTree->Branch("fPrimaryPx", &fPrimaryPx, "fPrimaryPx/D");
1419 fPrimaryTree->Branch("fPrimaryPy", &fPrimaryPy, "fPrimaryPy/D");
1420 fPrimaryTree->Branch("fPrimaryPz", &fPrimaryPz, "fPrimaryPz/D");
1421
1422 fPrimaryTree->Branch("fPrimaryVx", &fPrimaryVx, "fPrimaryVx/D");
1423 fPrimaryTree->Branch("fPrimaryVy", &fPrimaryVy, "fPrimaryVy/D");
1424 fPrimaryTree->Branch("fPrimaryVz", &fPrimaryVz, "fPrimaryVz/D");
1425
1426 fPrimaryTree->Branch("fPrimaryAccepted", &fPrimaryAccepted, "fPrimaryAccepted/B");
b2c10007 1427 fPrimaryTree->Branch("fPrimaryMatched", &fPrimaryMatched, "fPrimaryMatched/B");
f61cec2f 1428
1429
1430 fPrimaryTree->Branch("fDepositedCode", &fDepositedCode, "fDepositedCode/I");
1431 fPrimaryTree->Branch("fDepositedCharge", &fDepositedCharge, "fDepositedCharge/I");
b2c10007 1432 fPrimaryTree->Branch("fDepositedE", &fDepositedE, "fDepositedE/D");
f61cec2f 1433 fPrimaryTree->Branch("fDepositedEt", &fDepositedEt, "fDepositedEt/D");
1434
1435 fPrimaryTree->Branch("fDepositedVx", &fDepositedVx, "fDepositedVx/D");
1436 fPrimaryTree->Branch("fDepositedVy", &fDepositedVy, "fDepositedVy/D");
1437 fPrimaryTree->Branch("fDepositedVz", &fDepositedVz, "fDepositedVz/D");
1438
b2c10007 1439 fPrimaryTree->Branch("fSecondary", &fSecondary, "fSecondary/I");
1440
1441
1442 fPrimaryTree->Branch("fReconstructedE", &fReconstructedE, "fReconstructedE/D");
1443 fPrimaryTree->Branch("fReconstructedEt", &fReconstructedEt, "fReconstructedEt/D");
1444
1445 fPrimaryTree->Branch("fClusterMult", &fClusterMult, "fClusterMult/I");
1446
1447
f61cec2f 1448 }
1449
9ef6f13f 1450 fHistGammasFound = new TH1F("fHistGammasFound", "fHistGammasFound",200, 0, 10);
1451 fHistGammasGenerated = new TH1F("fHistGammasGenerated", "fHistGammasGenerated",200, 0, 10);
6a152780 1452 fHistGammasFoundMult = new TH2F("fHistGammasFoundMult", "fHistGammasFoundMult",200, 0, 10,10,0,100);
1453 fHistGammasGeneratedMult = new TH2F("fHistGammasGeneratedMult", "fHistGammasGeneratedMult",200, 0, 10,10,0,100);
9ef6f13f 1454 fHistChargedTracksCut = new TH1F("fHistChargedTracksCut", "fHistChargedTracksCut",100, 0, 5);
1455 fHistChargedTracksAccepted = new TH1F("fHistChargedTracksAccepted", "fHistChargedTracksAccepted",100, 0, 5);
1456 fHistGammasCut = new TH1F("fHistGammasTracksCut", "fHistGammasTracksCut",100, 0, 5);
1457 fHistGammasAccepted = new TH1F("fHistGammasTracksAccepted", "fHistGammasTracksAccepted",100, 0, 5);
0861cc1f 1458
1459 if(fCalcTrackMatchVsMult){
1460 fHistChargedTracksCutMult = new TH2F("fHistChargedTracksCutMult", "fHistChargedTracksCutMult",100, 0, 5,10,0,100);
1461 fHistChargedTracksAcceptedMult = new TH2F("fHistChargedTracksAcceptedMult", "fHistChargedTracksAcceptedMult",100, 0, 5,10,0,100);
ac610b08 1462 fHistChargedTracksAcceptedLowPtCent = new TH2F("fHistChargedTracksAcceptedLowPtCent", "fHistChargedTracksAcceptedLowPtCent",100, 0, 5,20,0,20);
1463 fHistChargedTracksAcceptedLowPtCentNoAntiProtons = new TH2F("fHistChargedTracksAcceptedLowPtCentNoAntiProtons", "fHistChargedTracksAcceptedLowPtCentNoAntiProtons",100, 0, 5,20,0,20);
1464 fHistChargedTracksAcceptedLowPtCentAntiProtons = new TH2F("fHistChargedTracksAcceptedLowPtCentAntiProtons", "fHistChargedTracksAcceptedLowPtCentAntiProtons",100, 0, 5,20,0,20);
0861cc1f 1465 fHistGammasCutMult = new TH2F("fHistGammasTracksCutMult", "fHistGammasTracksCutMult",100, 0, 5,10,0,100);
1466 fHistGammasAcceptedMult = new TH2F("fHistGammasTracksAcceptedMult", "fHistGammasTracksAcceptedMult",100, 0, 5,10,0,100);
1467 }
1468
9ef6f13f 1469 fHistBadTrackMatches = new TH1F("fHistBadTrackMatches", "fHistBadTrackMatches",100, 0, 5);
1470 fHistMatchedTracksEvspTBkgd = new TH2F("fHistMatchedTracksEvspTBkgd", "fHistMatchedTracksEvspTBkgd",100, 0, 3,100,0,3);
1471 fHistMatchedTracksEvspTSignal = new TH2F("fHistMatchedTracksEvspTSignal", "fHistMatchedTracksEvspTSignal",100, 0, 3,100,0,3);
0861cc1f 1472 if(fCalcTrackMatchVsMult){
1473 fHistMatchedTracksEvspTBkgdPeripheral = new TH2F("fHistMatchedTracksEvspTBkgdPeripheral", "fHistMatchedTracksEvspTBkgd",100, 0, 3,100,0,3);
1474 fHistMatchedTracksEvspTSignalPeripheral = new TH2F("fHistMatchedTracksEvspTSignalPeripheral", "fHistMatchedTracksEvspTSignal",100, 0, 3,100,0,3);
1475
ac610b08 1476 fHistMatchedTracksEvspTBkgdvsCent = new TH3F("fHistMatchedTracksEvspTBkgdvsCent", "fHistMatchedTracksEvspTBkgdvsCent",100, 0, 3,100,0,3,20,0,20);
1477 fHistMatchedTracksEvspTSignalvsCent = new TH3F("fHistMatchedTracksEvspTSignalvsCent", "fHistMatchedTracksEvspTSignalvsCent",100, 0, 3,100,0,3,20,0,20);
1478 fHistMatchedTracksEvspTBkgdvsCentEffCorr = new TH3F("fHistMatchedTracksEvspTBkgdvsCentEffCorr", "fHistMatchedTracksEvspTBkgdvsCent",100, 0, 3,100,0,3,20,0,20);
1479 fHistMatchedTracksEvspTSignalvsCentEffCorr = new TH3F("fHistMatchedTracksEvspTSignalvsCentEffCorr", "fHistMatchedTracksEvspTSignalvsCent",100, 0, 3,100,0,3,20,0,20);
0861cc1f 1480
1481
1482 fHistChargedTracksCutPeripheral = new TH1F("fHistChargedTracksCutPeripheral", "fHistChargedTracksCut",100, 0, 5);
1483 fHistChargedTracksAcceptedPeripheral = new TH1F("fHistChargedTracksAcceptedPeripheral", "fHistChargedTracksAccepted",100, 0, 5);
1484 fHistGammasCutPeripheral = new TH1F("fHistGammasTracksCutPeripheral", "fHistGammasTracksCut",100, 0, 5);
1485 fHistGammasAcceptedPeripheral = new TH1F("fHistGammasTracksAcceptedPeripheral", "fHistGammasTracksAccepted",100, 0, 5);
1486 }
9ef6f13f 1487 fHistBadTrackMatchesdPhidEta = new TH2F("fHistBadTrackMatchesdPhidEta", "fHistBadTrackMatchesdPhidEta",20, -0.1, 0.1,20,-.1,0.1);
1488 fHistGoodTrackMatchesdPhidEta = new TH2F("fHistGoodTrackMatchesdPhidEta", "fHistGoodTrackMatchesdPhidEta",20, -0.1, 0.1,20,-.1,0.1);
c79c36cc 1489
1490 fHistHadronDepositsAll = new TH1F("fHistHadronDepositsAll","All Hadrons which deposited energy in calorimeter",fgNumOfPtBins,fgPtAxis);
1491 fHistHadronDepositsReco = new TH1F("fHistHadronDepositsReco","Reconstructed Hadrons which deposited energy in calorimeter",fgNumOfPtBins,fgPtAxis);
1492 //,10,0,100
1493 Int_t nMult = 20;
1494 Float_t nMultCuts[21] = { 0, 5,10,15,20, 25,30,35,40,45,
1495 50,55,60,65,70, 75,80,85,90,95,
1496 100};
ac610b08 1497 Int_t nCent = 20;
1498 Float_t nCentCuts[21] = { 0, 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20};
c79c36cc 1499
ac610b08 1500 fHistHadronDepositsAllCent = new TH2F("fHistHadronDepositsAllCent","All Hadrons which deposited energy in calorimeter",fgNumOfPtBins,fgPtAxis,nCent,nCentCuts);
1501 fHistHadronDepositsRecoCent = new TH2F("fHistHadronDepositsRecoCent","Reconstructed Hadrons which deposited energy in calorimeter",fgNumOfPtBins,fgPtAxis,nCent,nCentCuts);
d3ce32b8 1502
ac610b08 1503 fHistHadronsAllCent = new TH2F("fHistHadronsAllCent","All Hadrons vs cluster mult",fgNumOfPtBins,fgPtAxis,nCent,nCentCuts);
d3ce32b8 1504
1505 fHistMultChVsSignalVsMult = new TH3F("fHistMultChVsSignalVsMult","Charged particle Multiplicity vs Signal particle multiplicity vs Cluster Mult",nMult,nMultCuts,nMult,nMultCuts,nMult,nMultCuts);
ce546038 1506}
1507
0651f6b4 1508void AliAnalysisEtMonteCarlo::FillOutputList(TList *list)
f61cec2f 1509{ //fill the output list
e2ee5727 1510 AliAnalysisEt::FillOutputList(list);
1511
c79c36cc 1512 if(!fIsMC) return;
f61cec2f 1513 if(fCuts->GetHistMakeTree())
1514 {
1515 list->Add(fPrimaryTree);
1516 }
1517
e2ee5727 1518 list->Add(fHistRemovedOrNot);
1519
1520 list->Add(fHistEtNonRemovedProtons);
1521 list->Add(fHistEtNonRemovedAntiProtons);
1522 list->Add(fHistEtNonRemovedPiPlus);
1523 list->Add(fHistEtNonRemovedPiMinus);
1524 list->Add(fHistEtNonRemovedKaonPlus);
1525 list->Add(fHistEtNonRemovedKaonMinus);
1526 list->Add(fHistEtNonRemovedK0s);
ef647350 1527 list->Add(fHistEtNonRemovedK0L);
e2ee5727 1528 list->Add(fHistEtNonRemovedLambdas);
1529 list->Add(fHistEtNonRemovedElectrons);
1530 list->Add(fHistEtNonRemovedPositrons);
1531 list->Add(fHistEtNonRemovedMuPlus);
1532 list->Add(fHistEtNonRemovedMuMinus);
1533 list->Add(fHistEtNonRemovedNeutrons);
1534 list->Add(fHistEtNonRemovedAntiNeutrons);
1535 list->Add(fHistEtNonRemovedGammas);
1536 list->Add(fHistEtNonRemovedGammasFromPi0);
1537
1538 list->Add(fHistEtRemovedGammas);
1539 list->Add(fHistEtRemovedNeutrons);
1540 list->Add(fHistEtRemovedAntiNeutrons);
1541
ef647350 1542 list->Add(fHistEtRemovedCharged);
1543 list->Add(fHistEtRemovedNeutrals);
1544
1545 list->Add(fHistEtNonRemovedCharged);
1546 list->Add(fHistEtNonRemovedNeutrals);
e2ee5727 1547
1548 list->Add(fHistMultNonRemovedProtons);
1549 list->Add(fHistMultNonRemovedAntiProtons);
1550 list->Add(fHistMultNonRemovedPiPlus);
1551 list->Add(fHistMultNonRemovedPiMinus);
1552 list->Add(fHistMultNonRemovedKaonPlus);
1553 list->Add(fHistMultNonRemovedKaonMinus);
1554 list->Add(fHistMultNonRemovedK0s);
ef647350 1555 list->Add(fHistMultNonRemovedK0L);
e2ee5727 1556 list->Add(fHistMultNonRemovedLambdas);
1557 list->Add(fHistMultNonRemovedElectrons);
1558 list->Add(fHistMultNonRemovedPositrons);
1559 list->Add(fHistMultNonRemovedMuPlus);
1560 list->Add(fHistMultNonRemovedMuMinus);
1561 list->Add(fHistMultNonRemovedNeutrons);
1562 list->Add(fHistMultNonRemovedAntiNeutrons);
1563 list->Add(fHistMultNonRemovedGammas);
1564
1565 list->Add(fHistMultRemovedGammas);
1566 list->Add(fHistMultRemovedNeutrons);
1567 list->Add(fHistMultRemovedAntiNeutrons);
1568
ef647350 1569 list->Add(fHistMultRemovedCharged);
1570 list->Add(fHistMultRemovedNeutrals);
1571
1572 list->Add(fHistMultNonRemovedCharged);
1573 list->Add(fHistMultNonRemovedNeutrals);
1574
e2ee5727 1575 list->Add(fHistTrackMultvsNonRemovedCharged);
1576 list->Add(fHistTrackMultvsNonRemovedNeutral);
1577 list->Add(fHistTrackMultvsRemovedGamma);
1578
1579 list->Add(fHistClusterMultvsNonRemovedCharged);
1580 list->Add(fHistClusterMultvsNonRemovedNeutral);
1581 list->Add(fHistClusterMultvsRemovedGamma);
1582
1583 //list->Add(fHistDecayVertexNonRemovedCharged);
1584 //list->Add(fHistDecayVertexNonRemovedNeutral);
1585 //list->Add(fHistDecayVertexRemovedCharged);
1586 //list->Add(fHistDecayVertexRemovedNeutral);
1587
1588 list->Add(fHistDxDzNonRemovedCharged);
1589 list->Add(fHistDxDzRemovedCharged);
1590 list->Add(fHistDxDzNonRemovedNeutral);
1591 list->Add(fHistDxDzRemovedNeutral);
1592
0861cc1f 1593 if(fCalcForKaonCorrection){
1594 list->Add(fHistK0EDepositsVsPtInAcceptance);
1595 list->Add(fHistK0EGammaVsPtInAcceptance);
1596 list->Add(fHistK0EDepositsVsPtOutOfAcceptance);
1597 list->Add(fHistK0EGammaVsPtOutOfAcceptance);
1598 list->Add(fHistSimKaonsInAcceptance);
1599 list->Add(fHistSimK0SInAcceptance);
1600 list->Add(fHistSimK0LInAcceptance);
1601 list->Add(fHistSimKPlusInAcceptance);
1602 list->Add(fHistSimKMinusInAcceptance);
1603 list->Add(fHistSimKaonsOutOfAcceptance);
1604 list->Add(fHistSimKaonsInAcceptanceWithDepositsPrimaries);
1605 list->Add(fHistSimKaonsOutOfAcceptanceWithDepositsSecondaries);
1606 list->Add(fHistSimKaonsOutOfAcceptanceWithDepositsPrimaries);
1607 }
2aab9269 1608
e2ee5727 1609 list->Add(fHistPiPlusMult);
1610 list->Add(fHistPiMinusMult);
1611 list->Add(fHistPiZeroMult);
1612 list->Add(fHistPiPlusMultAcc);
1613 list->Add(fHistPiMinusMultAcc);
1614 list->Add(fHistPiZeroMultAcc);
b2c10007 1615
1616 list->Add(fHistGammasFound);
1617 list->Add(fHistGammasGenerated);
6a152780 1618 list->Add(fHistGammasFoundMult);
1619 list->Add(fHistGammasGeneratedMult);
9ef6f13f 1620 list->Add(fHistChargedTracksCut);
1621 list->Add(fHistChargedTracksAccepted);
1622 list->Add(fHistGammasCut);
1623 list->Add(fHistGammasAccepted);
0861cc1f 1624 if(fCalcTrackMatchVsMult){
1625 list->Add(fHistChargedTracksCutMult);
1626 list->Add(fHistChargedTracksAcceptedMult);
ac610b08 1627 list->Add(fHistChargedTracksAcceptedLowPtCent);
1628 list->Add(fHistChargedTracksAcceptedLowPtCentNoAntiProtons);
1629 list->Add(fHistChargedTracksAcceptedLowPtCentAntiProtons);
0861cc1f 1630 list->Add(fHistGammasCutMult);
1631 list->Add(fHistGammasAcceptedMult);
1632 }
9ef6f13f 1633 list->Add(fHistBadTrackMatches);
1634 list->Add(fHistMatchedTracksEvspTBkgd);
1635 list->Add(fHistMatchedTracksEvspTSignal);
0861cc1f 1636 if(fCalcTrackMatchVsMult){
1637 list->Add(fHistMatchedTracksEvspTBkgdPeripheral);
1638 list->Add(fHistMatchedTracksEvspTSignalPeripheral);
ac610b08 1639 list->Add(fHistMatchedTracksEvspTBkgdvsCent);
1640 list->Add(fHistMatchedTracksEvspTSignalvsCent);
1641 list->Add(fHistMatchedTracksEvspTBkgdvsCentEffCorr);
1642 list->Add(fHistMatchedTracksEvspTSignalvsCentEffCorr);
0861cc1f 1643 list->Add(fHistChargedTracksCutPeripheral);
1644 list->Add(fHistChargedTracksAcceptedPeripheral);
1645 list->Add(fHistGammasCutPeripheral);
1646 list->Add(fHistGammasAcceptedPeripheral);
1647 }
9ef6f13f 1648 list->Add(fHistBadTrackMatchesdPhidEta);
1649 list->Add(fHistGoodTrackMatchesdPhidEta);
c79c36cc 1650 list->Add(fHistHadronDepositsAll);
1651 list->Add(fHistHadronDepositsReco);
ac610b08 1652 list->Add(fHistHadronDepositsAllCent);
1653 list->Add(fHistHadronDepositsRecoCent);
1654 list->Add(fHistHadronsAllCent);
d3ce32b8 1655 list->Add(fHistMultChVsSignalVsMult);
e2ee5727 1656
0651f6b4 1657}
1658
1659
13b0d3c1 1660bool AliAnalysisEtMonteCarlo::TrackHitsCalorimeter(TParticle* part, Double_t magField)
1661{
e2ee5727 1662 // printf(" TrackHitsCalorimeter - magField %f\n", magField);
1663 AliESDtrack *esdTrack = new AliESDtrack(part);
1664 // Printf("MC Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
1665
1666 Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
1667
1668 // if(prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
1669
f61cec2f 1670 bool status = prop && fSelector->CutGeometricalAcceptance(*esdTrack);
e2ee5727 1671 delete esdTrack;
1672
1673 return status;
1674}
1675
1676void AliAnalysisEtMonteCarlo::FillHistograms()
f61cec2f 1677{ // let base class fill its histograms, and us fill the local ones
e2ee5727 1678 AliAnalysisEt::FillHistograms();
c79c36cc 1679 if(!fIsMC) return;
e2ee5727 1680 //std::cout << fEtNonRemovedPiPlus << " " << fCentClass << std::endl;
1681
1682 fHistEtNonRemovedProtons->Fill(fEtNonRemovedProtons, fCentClass);
1683 fHistEtNonRemovedAntiProtons->Fill(fEtNonRemovedAntiProtons, fCentClass);
1684 fHistEtNonRemovedKaonPlus->Fill(fEtNonRemovedKaonPlus, fCentClass);
1685 fHistEtNonRemovedKaonMinus->Fill(fEtNonRemovedKaonMinus, fCentClass);
ef647350 1686 fHistEtNonRemovedK0s->Fill(fEtNonRemovedK0S, fCentClass);
1687 fHistEtNonRemovedK0L->Fill(fEtNonRemovedK0L, fCentClass);
e2ee5727 1688 fHistEtNonRemovedLambdas->Fill(fEtNonRemovedLambdas, fCentClass);
1689 fHistEtNonRemovedPiPlus->Fill(fEtNonRemovedPiPlus, fCentClass);
1690 fHistEtNonRemovedPiMinus->Fill(fEtNonRemovedPiMinus, fCentClass);
1691 fHistEtNonRemovedElectrons->Fill(fEtNonRemovedElectrons, fCentClass);
1692 fHistEtNonRemovedPositrons->Fill(fEtNonRemovedPositrons, fCentClass);
1693 fHistEtNonRemovedMuPlus->Fill(fEtNonRemovedMuPlus, fCentClass);
1694 fHistEtNonRemovedMuMinus->Fill(fEtNonRemovedMuMinus, fCentClass);
1695 fHistEtNonRemovedNeutrons->Fill(fEtNonRemovedNeutrons, fCentClass);
1696 fHistEtNonRemovedAntiNeutrons->Fill(fEtNonRemovedAntiNeutrons, fCentClass);
1697 fHistEtNonRemovedGammas->Fill(fEtNonRemovedGammas, fCentClass);
1698 fHistEtNonRemovedGammasFromPi0->Fill(fEtNonRemovedGammasFromPi0, fCentClass);
1699
ef647350 1700 fHistEtRemovedGammas->Fill(fEtRemovedGammas, fNClusters);
e2ee5727 1701 fHistEtRemovedNeutrons->Fill(fEtRemovedNeutrons, fCentClass);
1702 fHistEtRemovedAntiNeutrons->Fill(fEtRemovedAntiNeutrons, fCentClass);
1703
ef647350 1704// fHistEtRemovedCharged->Fill(fEtRemovedAntiProtons+fEtRemovedElectrons+fEtRemovedKaonMinus+fEtRemovedKaonPlus
1705// +fEtRemovedMuMinus+fEtRemovedMuPlus+fEtRemovedPiMinus+fEtRemovedPiPlus+fEtRemovedPositrons
1706// +fEtRemovedProtons.
1707// fCentClass);
1708// fHistEtRemovedNeutrals->Fill(fEtRemovedNeutrons+fEtRemovedAntiNeutrons, fCentClass);
f61cec2f 1709//
ef647350 1710// fHistEtNonRemovedCharged->Fill(fEtNonRemovedAntiProtons+fEtNonRemovedElectrons+fEtNonRemovedKaonMinus+fEtNonRemovedKaonPlus
1711// +fEtNonRemovedMuMinus+fEtNonRemovedMuPlus+fEtNonRemovedPiMinus+fEtNonRemovedPiPlus+fEtNonRemovedPositrons
1712// +fEtNonRemovedProtons,
1713// fCentClass);
1714// fHistEtRemovedNeutrals->Fill(fEtNonRemovedNeutrons+fEtNonRemovedAntiNeutrons, fCentClass);
1715
1716 fHistEtRemovedCharged->Fill(fEnergyChargedRemoved, fNClusters);
1717 fHistEtRemovedNeutrals->Fill(fEnergyNeutralRemoved, fNClusters);
1718 fHistEtNonRemovedCharged->Fill(fEnergyChargedNotRemoved, fNClusters);
1719 fHistEtNonRemovedNeutrals->Fill(fEnergyNeutralNotRemoved, fNClusters);
f61cec2f 1720
ef647350 1721 fHistMultRemovedCharged->Fill(fChargedRemoved, fNClusters);
1722 fHistMultRemovedNeutrals->Fill(fNeutralRemoved, fNClusters);
1723 fHistMultNonRemovedCharged->Fill(fChargedNotRemoved, fNClusters);
1724 fHistMultNonRemovedNeutrals->Fill(fNeutralNotRemoved, fNClusters);
f61cec2f 1725
1726
e2ee5727 1727 fHistMultNonRemovedProtons->Fill(fMultNonRemovedProtons, fCentClass);
1728 fHistMultNonRemovedAntiProtons->Fill(fMultNonRemovedAntiProtons, fCentClass);
1729 fHistMultNonRemovedKaonPlus->Fill(fMultNonRemovedKaonPlus, fCentClass);
1730 fHistMultNonRemovedKaonMinus->Fill(fMultNonRemovedKaonMinus, fCentClass);
1731 fHistMultNonRemovedK0s->Fill(fMultNonRemovedK0s, fCentClass);
ef647350 1732 fHistMultNonRemovedK0L->Fill(fMultNonRemovedK0L, fCentClass);
e2ee5727 1733 fHistMultNonRemovedLambdas->Fill(fMultNonRemovedLambdas, fCentClass);
1734 fHistMultNonRemovedPiPlus->Fill(fMultNonRemovedPiPlus, fCentClass);
1735 fHistMultNonRemovedPiMinus->Fill(fMultNonRemovedPiMinus, fCentClass);
1736 fHistMultNonRemovedElectrons->Fill(fMultNonRemovedElectrons, fCentClass);
1737 fHistMultNonRemovedPositrons->Fill(fMultNonRemovedPositrons, fCentClass);
1738 fHistMultNonRemovedMuPlus->Fill(fMultNonRemovedMuPlus, fCentClass);
1739 fHistMultNonRemovedMuMinus->Fill(fMultNonRemovedMuMinus, fCentClass);
1740 fHistMultNonRemovedNeutrons->Fill(fMultNonRemovedNeutrons, fCentClass);
1741 fHistMultNonRemovedAntiNeutrons->Fill(fMultNonRemovedAntiNeutrons, fCentClass);
1742 fHistMultNonRemovedGammas->Fill(fMultNonRemovedGammas, fCentClass);
1743
1744 fHistMultRemovedGammas->Fill(fMultRemovedGammas, fCentClass);
1745 fHistMultRemovedNeutrons->Fill(fMultRemovedNeutrons, fCentClass);
1746 fHistMultRemovedAntiNeutrons->Fill(fMultRemovedAntiNeutrons, fCentClass);
1747
1748 fHistTrackMultvsNonRemovedCharged->Fill(fTrackMultInAcc,
1749 fMultNonRemovedAntiProtons+fMultNonRemovedElectrons+fMultNonRemovedKaonMinus+fMultNonRemovedKaonPlus
1750 +fMultNonRemovedMuMinus+fMultNonRemovedMuPlus+fMultNonRemovedPiMinus+fMultNonRemovedPiPlus+fMultNonRemovedPositrons
1751 +fMultNonRemovedProtons);
1752
1753 fHistTrackMultvsNonRemovedNeutral->Fill(fTrackMultInAcc,
b2c10007 1754 fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedK0L+fMultNonRemovedLambdas+fK0sMult);
e2ee5727 1755
1756 fHistTrackMultvsRemovedGamma->Fill(fTrackMultInAcc,
1757 fMultRemovedGammas);
1758
ef647350 1759 fHistClusterMultvsNonRemovedCharged->Fill(fNClusters,
e2ee5727 1760 fMultNonRemovedAntiProtons+fMultNonRemovedElectrons+fMultNonRemovedKaonMinus
1761 +fMultNonRemovedKaonPlus+fMultNonRemovedMuMinus+fMultNonRemovedMuPlus
1762 +fMultNonRemovedPiMinus+fMultNonRemovedPiPlus+fMultNonRemovedPositrons+fMultNonRemovedProtons);
1763
ef647350 1764 fHistClusterMultvsNonRemovedNeutral->Fill(fNClusters,
b2c10007 1765 fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedK0L+fMultNonRemovedLambdas+fK0sMult);
e2ee5727 1766
ef647350 1767 fHistClusterMultvsRemovedGamma->Fill(fNClusters,
e2ee5727 1768 fMultRemovedGammas);
1769
13b0d3c1 1770}
1771
e2ee5727 1772
1773
1774
ef647350 1775Int_t AliAnalysisEtMonteCarlo::PrintFamilyTree(Int_t partIdx, AliStack* stack)
4d376d01 1776{ // print family tree
f61cec2f 1777 TParticle *part = stack->Particle(partIdx);
b2c10007 1778// if(part->GetPdgCode() == fgK0SCode)
f61cec2f 1779 {
1780 std::cout << "This is index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") , is it primary: " << stack->IsPhysicalPrimary(partIdx)<< std::endl;
1781 std::cout << "PID: " << part->GetPdgCode() << "/" << part->GetName() << std::endl;
1782 std::cout << "Energy: " << part->Energy() << std::endl;
1783 std::cout << "Vertex: " << part->Vx() << ", " << part->Vy() << ", " << part->Vz() << std::endl;
1784 }
1785 return PrintMothers(partIdx, stack, 1);
ef647350 1786}
1787
1788Int_t AliAnalysisEtMonteCarlo::PrintMothers(Int_t partIdx, AliStack* stack, Int_t gen)
4d376d01 1789{ // print mothers
f61cec2f 1790 char *tabs = new char[gen+1];
1791 for(Int_t i = 0; i < gen; ++i)
1792 {
1793 //std::cout << i << std::endl;
1794 tabs[i] = '\t';
1795 }
1796 tabs[gen] = '\0';
1797 Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
1798 if(mothIdx < 0)
1799 {
1800 delete [] tabs;
1801 return 0;
1802 }
1803 TParticle *mother = stack->Particle(mothIdx);
b2c10007 1804// if(mother->GetPdgCode() == fgK0SCode)
f61cec2f 1805 {
b2c10007 1806 //std::cout << tabs << "Mother of index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") is: " << mothIdx << ", is it primary: " << stack->IsPhysicalPrimary(mothIdx)<< std::endl;
1807 std::cout << tabs << "Index: " << mothIdx << std::endl;
1808 std::cout << tabs << "Primary: " << stack->IsPhysicalPrimary(mothIdx) << std::endl;
f61cec2f 1809 std::cout << tabs << "PID: " << mother->GetPdgCode() << "/" << mother->GetName() << std::endl;
1810 std::cout << tabs << "Energy: " << mother->Energy() << std::endl;
b2c10007 1811 if(mother->GetFirstMother() >= 0)
1812 {
1813 std::cout << tabs << "Mother(s): " << stack->Particle(mother->GetFirstMother())->GetPdgCode();
1814 if(mother->GetSecondMother() >= 0) std::cout << ", " << stack->Particle(mother->GetSecondMother())->GetPdgCode();
1815 std::cout << std::endl;
1816 }
f61cec2f 1817 std::cout << tabs << "Vertex: " << mother->Vx() << ", " << mother->Vy() << ", " << mother->Vz() << std::endl;
1818 }
1819 if(mother->GetPdgCode() == fgK0SCode)
1820 {
ef647350 1821// std::cout << "K0S!!!!!!!!!!!!!11111!!!!!" << std::endl;
f61cec2f 1822 }
ef647350 1823// std::cout << "Mother of index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") is: " << mothIdx << std::endl;
1824// std::cout << "PID: " << mother->GetPdgCode() << "/" << mother->GetName() << std::endl;
1825// std::cout << "Energy: " << mother->Energy() << std::endl;
1826// std::cout << "Vertex: " << mother->Vx() << ", " << mother->Vy() << ", " << mother->Vz() << std::endl;
f61cec2f 1827
1828 delete [] tabs;
1829 return PrintMothers(mothIdx, stack, gen+1) + 1;
ef647350 1830}
1831
1832Int_t AliAnalysisEtMonteCarlo::GetPrimMother(Int_t partIdx, AliStack *stack)
4d376d01 1833{ // get primary mother
f61cec2f 1834 if(partIdx >= 0)
ef647350 1835 {
b2c10007 1836 //return stack->GetPrimary(partIdx);
1837
f61cec2f 1838 Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
1839 if(mothIdx < 0) return -1;
1840 TParticle *mother = stack->Particle(mothIdx);
1841 if(mother)
1842 {
f61cec2f 1843 if(stack->IsPhysicalPrimary(mothIdx)) return mothIdx;
1844 else return GetPrimMother(mothIdx, stack);
1845 }
1846 else
1847 {
1848 return -1;
1849 }
ef647350 1850 }
f61cec2f 1851 return -1;
ef647350 1852}
1853
1854Int_t AliAnalysisEtMonteCarlo::GetK0InFamily(Int_t partIdx, AliStack* stack)
4d376d01 1855{ // get K0 in family
f61cec2f 1856 if(partIdx >= 0)
ef647350 1857 {
2aab9269 1858 if(stack->Particle(partIdx)->GetPdgCode() == fgK0SCode) return partIdx;
f61cec2f 1859 Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
1860 if(mothIdx < 0) return -1;
1861 TParticle *mother = stack->Particle(mothIdx);
1862 if(mother)
ef647350 1863 {
2aab9269 1864 if(mother->GetPdgCode() == fgK0SCode)
ef647350 1865 {
f61cec2f 1866 return mothIdx;
ef647350 1867 }
f61cec2f 1868 return GetK0InFamily(mothIdx, stack);
1869 }
1870 else
1871 {
1872 return -1;
ef647350 1873 }
1874 }
f61cec2f 1875 return -1;
ef647350 1876}
e2ee5727 1877