]> 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
Fixed permissions
[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"
e2ee5727 32using namespace std;
16abb579 33
34ClassImp(AliAnalysisEtMonteCarlo);
35
36
ce546038 37// ctor
0651f6b4 38AliAnalysisEtMonteCarlo::AliAnalysisEtMonteCarlo():AliAnalysisEt()
f61cec2f 39 ,fImpactParameter(0)
40 ,fNcoll(0)
41 ,fNpart(0)
42 ,fPrimaryTree(0)
43 ,fTotEtWithSecondaryRemoved(0)
44 ,fTotEtSecondaryFromEmEtPrimary(0)
45 ,fTotEtSecondary(0)
46 ,fPrimaryCode(0)
47 ,fPrimaryCharge(0)
48 ,fPrimaryE(0)
49 ,fPrimaryEt(0)
50 ,fPrimaryPx(0)
51 ,fPrimaryPy(0)
52 ,fPrimaryPz(0)
53 ,fPrimaryVx(0)
54 ,fPrimaryVy(0)
55 ,fPrimaryVz(0)
56 ,fPrimaryAccepted(0)
b2c10007 57 ,fPrimaryMatched(0)
f61cec2f 58 ,fDepositedCode(0)
b2c10007 59 ,fDepositedE(0)
f61cec2f 60 ,fDepositedEt(0)
61 ,fDepositedCharge(0)
62 ,fDepositedVx(0)
63 ,fDepositedVy(0)
64 ,fDepositedVz(0)
b2c10007 65 ,fSecondary(kFALSE)
66 ,fReconstructedE(0)
67 ,fReconstructedEt(0)
68 ,fTotPx(0)
69 ,fTotPy(0)
70 ,fTotPz(0)
71 ,fClusterMult(0)
f61cec2f 72 ,fHistDecayVertexNonRemovedCharged(0)
73 ,fHistDecayVertexRemovedCharged(0)
74 ,fHistDecayVertexNonRemovedNeutral(0)
75 ,fHistDecayVertexRemovedNeutral(0)
76
77 ,fHistRemovedOrNot(0)
78 ,fHistEtNonRemovedProtons(0)
79 ,fHistEtNonRemovedAntiProtons(0)
80 ,fHistEtNonRemovedPiPlus(0)
81 ,fHistEtNonRemovedPiMinus(0)
82 ,fHistEtNonRemovedKaonPlus(0)
83 ,fHistEtNonRemovedKaonMinus(0)
84 ,fHistEtNonRemovedK0s(0)
85 ,fHistEtNonRemovedK0L(0)
86 ,fHistEtNonRemovedLambdas(0)
87 ,fHistEtNonRemovedElectrons(0)
88 ,fHistEtNonRemovedPositrons(0)
89 ,fHistEtNonRemovedMuPlus(0)
90 ,fHistEtNonRemovedMuMinus(0)
91 ,fHistEtNonRemovedNeutrons(0)
92 ,fHistEtNonRemovedAntiNeutrons(0)
93 ,fHistEtNonRemovedGammas(0)
94 ,fHistEtNonRemovedGammasFromPi0(0)
95 ,fHistEtRemovedGammas(0)
96 ,fHistEtRemovedNeutrons(0)
97 ,fHistEtRemovedAntiNeutrons(0)
98 ,fHistEtRemovedCharged(0)
99 ,fHistEtRemovedNeutrals(0)
100 ,fHistEtNonRemovedCharged(0)
101 ,fHistEtNonRemovedNeutrals(0)
102 ,fHistMultNonRemovedProtons(0)
103 ,fHistMultNonRemovedAntiProtons(0)
104 ,fHistMultNonRemovedPiPlus(0)
105 ,fHistMultNonRemovedPiMinus(0)
106 ,fHistMultNonRemovedKaonPlus(0)
107 ,fHistMultNonRemovedKaonMinus(0)
108 ,fHistMultNonRemovedK0s(0)
109 ,fHistMultNonRemovedK0L(0)
110 ,fHistMultNonRemovedLambdas(0)
111 ,fHistMultNonRemovedElectrons(0)
112 ,fHistMultNonRemovedPositrons(0)
113 ,fHistMultNonRemovedMuPlus(0)
114 ,fHistMultNonRemovedMuMinus(0)
115 ,fHistMultNonRemovedNeutrons(0)
116 ,fHistMultNonRemovedAntiNeutrons(0)
117 ,fHistMultNonRemovedGammas(0)
118 ,fHistMultRemovedGammas(0)
119 ,fHistMultRemovedNeutrons(0)
120 ,fHistMultRemovedAntiNeutrons(0)
121 ,fHistMultRemovedCharged(0)
122 ,fHistMultRemovedNeutrals(0)
123 ,fHistMultNonRemovedCharged(0)
124 ,fHistMultNonRemovedNeutrals(0)
125 ,fHistTrackMultvsNonRemovedCharged(0)
126 ,fHistTrackMultvsNonRemovedNeutral(0)
127 ,fHistTrackMultvsRemovedGamma(0)
128 ,fHistClusterMultvsNonRemovedCharged(0)
129 ,fHistClusterMultvsNonRemovedNeutral(0)
130 ,fHistClusterMultvsRemovedGamma(0)
131 ,fHistMultvsNonRemovedChargedE(0)
132 ,fHistMultvsNonRemovedNeutralE(0)
03bfd268 133 ,fHistMultvsRemovedGammaE(0)
2aab9269 134 ,fHistK0EDepositsVsPtForEOver000MeV(0)
135 ,fHistK0EDepositsVsPtForEOver100MeV(0)
136 ,fHistK0EDepositsVsPtForEOver150MeV(0)
137 ,fHistK0EDepositsVsPtForEOver200MeV(0)
138 ,fHistK0EDepositsVsPtForEOver250MeV(0)
139 ,fHistK0EDepositsVsPtForEOver300MeV(0)
140 ,fHistK0EDepositsVsPtForEOver350MeV(0)
141 ,fHistK0EDepositsVsPtForEOver400MeV(0)
142 ,fHistK0EDepositsVsPtForEOver450MeV(0)
143 ,fHistK0EDepositsVsPtForEOver500MeV(0)
144 ,fHistK0EGammaVsPtForEOver000MeV(0)
145 ,fHistK0EGammaVsPtForEOver100MeV(0)
146 ,fHistK0EGammaVsPtForEOver150MeV(0)
147 ,fHistK0EGammaVsPtForEOver200MeV(0)
148 ,fHistK0EGammaVsPtForEOver250MeV(0)
149 ,fHistK0EGammaVsPtForEOver300MeV(0)
150 ,fHistK0EGammaVsPtForEOver350MeV(0)
151 ,fHistK0EGammaVsPtForEOver400MeV(0)
152 ,fHistK0EGammaVsPtForEOver450MeV(0)
153 ,fHistK0EGammaVsPtForEOver500MeV(0)
f61cec2f 154 ,fEtNonRemovedProtons(0)
155 ,fEtNonRemovedAntiProtons(0)
156 ,fEtNonRemovedPiPlus(0)
157 ,fEtNonRemovedPiMinus(0)
158 ,fEtNonRemovedKaonPlus(0)
159 ,fEtNonRemovedKaonMinus(0)
160 ,fEtNonRemovedK0S(0)
161 ,fEtNonRemovedK0L(0)
162 ,fEtNonRemovedLambdas(0)
163 ,fEtNonRemovedElectrons(0)
164 ,fEtNonRemovedPositrons(0)
165 ,fEtNonRemovedMuMinus(0)
166 ,fEtNonRemovedMuPlus(0)
167 ,fEtNonRemovedGammas(0)
168 ,fEtNonRemovedGammasFromPi0(0)
169 ,fEtNonRemovedNeutrons(0)
170 ,fEtNonRemovedAntiNeutrons(0)
171 ,fEtRemovedProtons(0)
172 ,fEtRemovedAntiProtons(0)
173 ,fEtRemovedPiPlus(0)
174 ,fEtRemovedPiMinus(0)
175 ,fEtRemovedKaonPlus(0)
176 ,fEtRemovedKaonMinus(0)
177 ,fEtRemovedK0s(0)
178 ,fEtRemovedK0L(0)
179 ,fEtRemovedLambdas(0)
180 ,fEtRemovedElectrons(0)
181 ,fEtRemovedPositrons(0)
182 ,fEtRemovedMuMinus(0)
183 ,fEtRemovedMuPlus(0)
184 ,fEtRemovedGammasFromPi0(0)
185 ,fEtRemovedGammas(0)
186 ,fEtRemovedNeutrons(0)
187 ,fEtRemovedAntiNeutrons(0)
188 ,fMultNonRemovedProtons(0)
189 ,fMultNonRemovedAntiProtons(0)
190 ,fMultNonRemovedPiPlus(0)
191 ,fMultNonRemovedPiMinus(0)
192 ,fMultNonRemovedKaonPlus(0)
193 ,fMultNonRemovedKaonMinus(0)
194 ,fMultNonRemovedK0s(0)
195 ,fMultNonRemovedK0L(0)
196 ,fMultNonRemovedLambdas(0)
197 ,fMultNonRemovedElectrons(0)
198 ,fMultNonRemovedPositrons(0)
199 ,fMultNonRemovedMuMinus(0)
200 ,fMultNonRemovedMuPlus(0)
201 ,fMultNonRemovedGammas(0)
202 ,fMultNonRemovedNeutrons(0)
203 ,fMultNonRemovedAntiNeutrons(0)
204 ,fMultRemovedProtons(0)
205 ,fMultRemovedAntiProtons(0)
206 ,fMultRemovedPiPlus(0)
207 ,fMultRemovedPiMinus(0)
208 ,fMultRemovedKaonPlus(0)
209 ,fMultRemovedKaonMinus(0)
210 ,fMultRemovedK0s(0)
211 ,fMultRemovedK0L(0)
212 ,fMultRemovedLambdas(0)
213 ,fMultRemovedElectrons(0)
214 ,fMultRemovedPositrons(0)
215 ,fMultRemovedMuMinus(0)
216 ,fMultRemovedMuPlus(0)
217 ,fMultRemovedGammas(0)
218 ,fMultRemovedNeutrons(0)
219 ,fMultRemovedAntiNeutrons(0)
220 ,fTrackMultInAcc(0)
221 ,fHistDxDzNonRemovedCharged(0)
222 ,fHistDxDzRemovedCharged(0)
223 ,fHistDxDzNonRemovedNeutral(0)
224 ,fHistDxDzRemovedNeutral(0)
225 ,fHistPiPlusMult(0)
226 ,fHistPiMinusMult(0)
227 ,fHistPiZeroMult(0)
228 ,fHistPiPlusMultAcc(0)
229 ,fHistPiMinusMultAcc(0)
230 ,fHistPiZeroMultAcc(0)
b2c10007 231// ,fPiPlusMult(0)
232// ,fPiMinusMult(0)
f61cec2f 233 ,fPiZeroMult(0)
234 ,fPiPlusMultAcc(0)
235 ,fPiMinusMultAcc(0)
236 ,fPiZeroMultAcc(0)
237 ,fNeutralRemoved(0)
238 ,fChargedRemoved(0)
239 ,fChargedNotRemoved(0)
240 ,fNeutralNotRemoved(0)
b2c10007 241 ,fGammaRemoved(0)
242 ,fSecondaryNotRemoved(0)
f61cec2f 243 ,fEnergyNeutralRemoved(0)
244 ,fEnergyChargedRemoved(0)
245 ,fEnergyChargedNotRemoved(0)
246 ,fEnergyNeutralNotRemoved(0)
b2c10007 247 ,fEnergyGammaRemoved(0)
f61cec2f 248 ,fNClusters(0)
249 ,fTotNeutralEtAfterMinEnergyCut(0)
b2c10007 250 ,fHistGammasFound(0)
251 ,fHistGammasGenerated(0)
ce546038 252{
253}
254
255// dtor
e2ee5727 256AliAnalysisEtMonteCarlo::~AliAnalysisEtMonteCarlo()
f61cec2f 257{ //Destructor
311c6540 258
259 if(fPrimaryTree){
260 fPrimaryTree->Clear();
261 delete fPrimaryTree;
262 }
f61cec2f 263 delete fHistDecayVertexNonRemovedCharged; // Decay vertex for non-removed charged particles
264 delete fHistDecayVertexRemovedCharged; // Decay vertex for non-removed charged particles
265 delete fHistDecayVertexNonRemovedNeutral; // Decay vertex for non-removed charged particles
266 delete fHistDecayVertexRemovedNeutral; // Decay vertex for non-removed charged particles
267
268 delete fHistRemovedOrNot; // If charged/neutral particles were removed or not
269
270 delete fHistEtNonRemovedProtons; // enter comment here
271 delete fHistEtNonRemovedAntiProtons; // enter comment here
272 delete fHistEtNonRemovedPiPlus; // enter comment here
273 delete fHistEtNonRemovedPiMinus; // enter comment here
274 delete fHistEtNonRemovedKaonPlus; // enter comment here
275 delete fHistEtNonRemovedKaonMinus; // enter comment here
276 delete fHistEtNonRemovedK0s; // enter comment here
277 delete fHistEtNonRemovedK0L; // enter comment here
278 delete fHistEtNonRemovedLambdas; // enter comment here
279 delete fHistEtNonRemovedElectrons; // enter comment here
280 delete fHistEtNonRemovedPositrons; // enter comment here
281 delete fHistEtNonRemovedMuPlus; // enter comment here
282 delete fHistEtNonRemovedMuMinus; // enter comment here
283 delete fHistEtNonRemovedNeutrons; // enter comment here
284 delete fHistEtNonRemovedAntiNeutrons; // enter comment here
285 delete fHistEtNonRemovedGammas; // enter comment here
286 delete fHistEtNonRemovedGammasFromPi0; // enter comment here
287
288 delete fHistEtRemovedGammas; // enter comment here
289 delete fHistEtRemovedNeutrons; // enter comment here
290 delete fHistEtRemovedAntiNeutrons; // enter comment here
291
292
293 delete fHistMultNonRemovedProtons; // enter comment here
294 delete fHistMultNonRemovedAntiProtons; // enter comment here
295 delete fHistMultNonRemovedPiPlus; // enter comment here
296 delete fHistMultNonRemovedPiMinus; // enter comment here
297 delete fHistMultNonRemovedKaonPlus; // enter comment here
298 delete fHistMultNonRemovedKaonMinus; // enter comment here
299 delete fHistMultNonRemovedK0s; // enter comment here
300 delete fHistMultNonRemovedK0L; // enter comment here
301 delete fHistMultNonRemovedLambdas; // enter comment here
302 delete fHistMultNonRemovedElectrons; // enter comment here
303 delete fHistMultNonRemovedPositrons; // enter comment here
304 delete fHistMultNonRemovedMuPlus; // enter comment here
305 delete fHistMultNonRemovedMuMinus; // enter comment here
306 delete fHistMultNonRemovedNeutrons; // enter comment here
307 delete fHistMultNonRemovedAntiNeutrons; // enter comment here
308 delete fHistMultNonRemovedGammas; // enter comment here
309
310 delete fHistMultRemovedGammas; // enter comment here
311 delete fHistMultRemovedNeutrons; // enter comment here
312 delete fHistMultRemovedAntiNeutrons; // enter comment here
313
314 delete fHistTrackMultvsNonRemovedCharged; // enter comment here
315 delete fHistTrackMultvsNonRemovedNeutral; // enter comment here
316 delete fHistTrackMultvsRemovedGamma; // enter comment here
317
318 delete fHistClusterMultvsNonRemovedCharged; // enter comment here
319 delete fHistClusterMultvsNonRemovedNeutral; // enter comment here
320 delete fHistClusterMultvsRemovedGamma; // enter comment here
321
322 delete fHistMultvsNonRemovedChargedE; // enter comment here
323 delete fHistMultvsNonRemovedNeutralE; // enter comment here
324 delete fHistMultvsRemovedGammaE; // enter comment here
2aab9269 325
326 delete fHistK0EDepositsVsPtForEOver000MeV; // enter comment here
327 delete fHistK0EDepositsVsPtForEOver100MeV; // enter comment here
328 delete fHistK0EDepositsVsPtForEOver150MeV; // enter comment here
329 delete fHistK0EDepositsVsPtForEOver200MeV; // enter comment here
330 delete fHistK0EDepositsVsPtForEOver250MeV; // enter comment here
331 delete fHistK0EDepositsVsPtForEOver300MeV; // enter comment here
332 delete fHistK0EDepositsVsPtForEOver350MeV; // enter comment here
333 delete fHistK0EDepositsVsPtForEOver400MeV; // enter comment here
334 delete fHistK0EDepositsVsPtForEOver450MeV; // enter comment here
335 delete fHistK0EDepositsVsPtForEOver500MeV; // enter comment here
336 delete fHistK0EGammaVsPtForEOver000MeV; // enter comment here
337 delete fHistK0EGammaVsPtForEOver100MeV; // enter comment here
338 delete fHistK0EGammaVsPtForEOver150MeV; // enter comment here
339 delete fHistK0EGammaVsPtForEOver200MeV; // enter comment here
340 delete fHistK0EGammaVsPtForEOver250MeV; // enter comment here
341 delete fHistK0EGammaVsPtForEOver300MeV; // enter comment here
342 delete fHistK0EGammaVsPtForEOver350MeV; // enter comment here
343 delete fHistK0EGammaVsPtForEOver400MeV; // enter comment here
344 delete fHistK0EGammaVsPtForEOver450MeV; // enter comment here
345 delete fHistK0EGammaVsPtForEOver500MeV; // enter comment here
346
f61cec2f 347 delete fHistDxDzNonRemovedCharged; // enter comment here
348 delete fHistDxDzRemovedCharged; // enter comment here
349 delete fHistDxDzNonRemovedNeutral; // enter comment here
350 delete fHistDxDzRemovedNeutral; // enter comment here
351
352 delete fHistPiPlusMult; // enter comment here
353 delete fHistPiMinusMult; // enter comment here
354 delete fHistPiZeroMult; // enter comment here
355
356 delete fHistPiPlusMultAcc; // enter comment here
357 delete fHistPiMinusMultAcc; // enter comment here
358 delete fHistPiZeroMultAcc; // enter comment here
4503e29d 359 delete fHistGammasFound; // enter comment here
360 delete fHistGammasGenerated; // enter comment here
ce546038 361}
2fbf38ac 362
363Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
f61cec2f 364{ // analyse MC event
e2ee5727 365 ResetEventValues();
366
b2c10007 367
e2ee5727 368 fPiPlusMult = 0;
369 fPiMinusMult = 0;
370 fPiZeroMult = 0;
371 if (fCentrality)
372 {
373 fCentClass = fCentrality->GetCentralityClass10(fCentralityMethod);
374
375 }
e2ee5727 376
377 // Get us an mc event
378 if (!ev) {
379 AliFatal("ERROR: Event does not exist");
380 return 0;
381 }
382 AliMCEvent *event = dynamic_cast<AliMCEvent*>(ev);
383 if (!event) {
384 AliFatal("ERROR: MC Event does not exist");
385 return 0;
386 }
387
388 Double_t protonMass =fgProtonMass;
389
390 // Hijing header
391 AliGenEventHeader* genHeader = event->GenEventHeader();
392 if (!genHeader) {
393 Printf("ERROR: Event generation header does not exist");
394 return 0;
395 }
396 AliGenHijingEventHeader* hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
397 if (hijingGenHeader) {
398 fImpactParameter = hijingGenHeader->ImpactParameter();
399 fNcoll = hijingGenHeader->HardScatters(); // or should this be some combination of NN() NNw() NwN() NwNw() ?
400 fNpart = hijingGenHeader->ProjectileParticipants() + hijingGenHeader->TargetParticipants();
f61cec2f 401 }
402
e2ee5727 403 // Let's play with the stack!
404 AliStack *stack = event->Stack();
405
406 Int_t nPrim = stack->GetNtrack();
407
408 Int_t partCount = 0;
409 for (Int_t iPart = 0; iPart < nPrim; iPart++)
2fbf38ac 410 {
e2ee5727 411
412 TParticle *part = stack->Particle(iPart);
b2c10007 413
414
e2ee5727 415
416 if (!part)
417 {
418 Printf("ERROR: Could not get particle %d", iPart);
419 continue;
420 }
e2ee5727 421 TParticlePDG *pdg = part->GetPDG(0);
e2ee5727 422
f61cec2f 423 if (!pdg)
2fbf38ac 424 {
f61cec2f 425 Printf("ERROR: Could not get particle PDG %d", iPart);
e2ee5727 426 continue;
2fbf38ac 427 }
428
e2ee5727 429 Double_t particleMassPart = 0; //The mass part in the Et calculation for this particle
f61cec2f 430 Int_t code = pdg->PdgCode();
b2c10007 431
432 if(stack->IsPhysicalPrimary(iPart))
433 {
434 fTotPx += part->Px();
435 fTotPy += part->Py();
436 fTotPz += part->Pz();
437 }
438
e2ee5727 439 // Check for reasonable (for now neutral and singly charged) charge on the particle
86e7d5db 440 if (fSelector->IsNeutralMcParticle(iPart,*stack,*pdg))
e2ee5727 441 {
442
443 fMultiplicity++;
b2c10007 444// PrintFamilyTree(iPart, stack);
f61cec2f 445//
e2ee5727 446 if (TMath::Abs(part->Eta()) < fCuts->GetCommonEtaCut())
2fbf38ac 447 {
f61cec2f 448 //Printf("Particle with eta: %f, pid: %d", part->Eta(), code);
e2ee5727 449 // calculate E_T
450 if (
ef647350 451 TMath::Abs(code) == fgProtonCode ||
452 TMath::Abs(code) == fgNeutronCode ||
453 TMath::Abs(code) == fgLambdaCode ||
454 TMath::Abs(code) == fgXiCode ||
455 TMath::Abs(code) == fgXi0Code ||
456 TMath::Abs(code) == fgOmegaCode
e2ee5727 457 )
458 {
ef647350 459 if (code > 0) {
e2ee5727 460 particleMassPart = - protonMass;
461 }
ef647350 462 if (code < 0) {
e2ee5727 463 particleMassPart = protonMass;
464 }
465 }
466 Double_t et = part->Energy() * TMath::Sin(part->Theta()) + particleMassPart;
467
e2ee5727 468
469 // Fill up total E_T counters for each particle species
ef647350 470 if (code == fgProtonCode || code == fgAntiProtonCode)
e2ee5727 471 {
e2ee5727 472 }
ef647350 473 if (code == fgPiPlusCode || code == fgPiMinusCode)
e2ee5727 474 {
ef647350 475 if (code == fgPiPlusCode)
e2ee5727 476 {
e2ee5727 477 }
478 else
479 {
e2ee5727 480 }
481 }
ef647350 482 if (code == fgGammaCode)
e2ee5727 483 {
e2ee5727 484 }
f61cec2f 485 if (code == fgKPlusCode)
486 {
487 }
488 if(code == fgKMinusCode)
e2ee5727 489 {
e2ee5727 490 }
ef647350 491 if (code == fgMuPlusCode || code == fgMuMinusCode)
e2ee5727 492 {
e2ee5727 493 }
ef647350 494 if (code == fgEPlusCode || code == fgEMinusCode)
e2ee5727 495 {
e2ee5727 496 }
e2ee5727 497 // some neutrals also
ef647350 498 if (code == fgNeutronCode)
e2ee5727 499 {
e2ee5727 500 }
ef647350 501 if (code == fgAntiNeutronCode)
e2ee5727 502 {
e2ee5727 503 }
ef647350 504 if (code == fgGammaCode)
e2ee5727 505 {
e2ee5727 506 }
507
508 // Neutral particles
ef647350 509 //if (TMath::Abs(pdg->Charge() - fCuts->GetMonteCarloNeutralParticle()) <1e-3 )
f61cec2f 510
ef647350 511 if(code == fgGammaCode || code == fgPi0Code || code == fgEtaCode)
e2ee5727 512 {
b2c10007 513 // PrintFamilyTree(iPart,stack);
f61cec2f 514 //Printf("Gamma, phi: %f, eta: %f, phi cut min: %f, phi cut max: %f, eta cut: %f", part->Phi(), part->Eta(), fPhiMinCutAcc, fPhiMaxCutAcc, fEtaCutAcc);
515 //if (et > fCuts->GetCommonClusterEnergyCut()) fTotNeutralEt += et;
e2ee5727 516
517 // inside EMCal acceptance
f61cec2f 518
519 //if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiMaxCutAcc && part->Phi() > fPhiMinCutAcc)
520
b2c10007 521 if(fSelector->CutGeometricalAcceptance(*part) )
e2ee5727 522 {
523 fNeutralMultiplicity++;
f61cec2f 524 fTotNeutralEt += et;
86e7d5db 525 if(fSelector->PassMinEnergyCut(*part))
f61cec2f 526 {
527 fTotNeutralEtAfterMinEnergyCut += et;
528 }
e2ee5727 529 if (part->Energy() > 0.05) partCount++;
e2ee5727 530 }
531 }
532 //Charged particles
533 else if (TMath::Abs( pdg->Charge() - fCuts->GetMonteCarloNeutralParticle())>1e-3 )
534 {
e2ee5727 535
536 // inside EMCal acceptance
f61cec2f 537 if (fSelector->CutGeometricalAcceptance(*part))
e2ee5727 538 {
f61cec2f 539
540 fChargedMultiplicity++;
541
542 fTotChargedEt += et;
e2ee5727 543
ef647350 544 if (code == fgProtonCode || code == fgAntiProtonCode)
e2ee5727 545 {
e2ee5727 546 }
ef647350 547 if (code == fgPiPlusCode || code == fgPiMinusCode)
e2ee5727 548 {
e2ee5727 549 }
ef647350 550 if (code == fgKPlusCode || code == fgKMinusCode)
e2ee5727 551 {
e2ee5727 552 }
ef647350 553 if (code == fgMuPlusCode || code == fgMuMinusCode)
e2ee5727 554 {
e2ee5727 555 }
f61cec2f 556
ef647350 557 if (code == fgEPlusCode || code == fgEMinusCode)
e2ee5727 558 {
f61cec2f 559 fTotNeutralEt += et; // calling electrons neutral
560 fTotChargedEt -= et;
e1fa1966 561 }
e2ee5727 562 } // inside EMCal acceptance
563
e2ee5727 564 if (TrackHitsCalorimeter(part)) // magnetic field info not filled?
565 {
566 if (pdg->Charge() > 0) fHistPhivsPtPos->Fill(part->Phi(),part->Pt());
567 else fHistPhivsPtNeg->Fill(part->Phi(), part->Pt());
568 }
569 }
570 }
571
572 }
2fbf38ac 573 }
b2c10007 574 // std::cout << "Total: p_x = " << fTotPx << ", p_y = " << fTotPy << ", p_z = " << fTotPz << std::endl;
e2ee5727 575 fTotEt = fTotChargedEt + fTotNeutralEt;
f61cec2f 576 //fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;//
e2ee5727 577 //std::cout << "Event done! # of particles: " << partCount << std::endl;
578 return 0;
2fbf38ac 579}
0651f6b4 580//Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliMCEvent* mcEvent,AliESDEvent* realEvent)
581Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
f61cec2f 582{ // analyse MC and real event info
e2ee5727 583 //if(!mcEvent || !realEvent){
584 if (!ev || !ev2) {
585 AliFatal("ERROR: Event does not exist");
586 return 0;
587 }
f61cec2f 588 AliAnalysisEt::AnalyseEvent(ev);
e2ee5727 589 AliMCEvent *mcEvent = dynamic_cast<AliMCEvent*>(ev);
590 AliESDEvent *realEvent = dynamic_cast<AliESDEvent*>(ev2);
591 if (!mcEvent || !realEvent) {
592 AliFatal("ERROR: mcEvent or realEvent does not exist");
b2c10007 593
e2ee5727 594 }
595
b2c10007 596 std::vector<Int_t> foundGammas;
597
f61cec2f 598 fSelector->SetEvent(realEvent);
599
600 AnalyseEvent(ev);
601
e2ee5727 602 AliStack *stack = mcEvent->Stack();
603
604 // get all detector clusters
f61cec2f 605 // TRefArray* caloClusters = new TRefArray();
606
607// if (fDetector == fCuts->GetDetectorEmcal()) realEvent->GetEMCALClusters( caloClusters );
608 //else if (fDetector == fCuts->GetDetectorPhos()) realEvent->GetPHOSClusters( caloClusters );
609 //else {
610 //AliFatal("Detector ID has not been specified");
611 //return -1;
612// }
613
4503e29d 614//Note that this only returns clusters for the selected detector. fSelector actually calls the right GetClusters... for the detector
f61cec2f 615 TRefArray *caloClusters = fSelector->GetClusters();
e2ee5727 616
e2ee5727 617 Int_t nCluster = caloClusters->GetEntries();
b2c10007 618 fClusterMult = nCluster;
ef647350 619 fNClusters = 0;
e2ee5727 620 // loop the clusters
621 for (int iCluster = 0; iCluster < nCluster; iCluster++ )
622 {
f61cec2f 623 Int_t cf = 0;
e2ee5727 624 AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
2aab9269 625 //Float_t caloE = caloCluster->E()
f61cec2f 626 fNClusters++;
b2c10007 627 const UInt_t iPart = (UInt_t)TMath::Abs(caloCluster->GetLabel());
ef647350 628 TParticle *part = stack->Particle(iPart);
e2ee5727 629
630 if (!part)
0651f6b4 631 {
e2ee5727 632 Printf("No MC particle %d", iCluster);
633 continue;
634 }
635
f61cec2f 636 int primIdx = iPart;
e2ee5727 637 if (!stack->IsPhysicalPrimary(iPart)) // check whether particle is primary. we keep secondary electron and gamma for testing.
638 {
f61cec2f 639 primIdx = GetPrimMother(iPart, stack);
b2c10007 640 if(primIdx != stack->GetPrimary(iPart))
641 {
642 //std::cout << primIdx << " = " << stack->GetPrimary(iPart) << std::endl;
643 //PrintFamilyTree(iPart, stack);
644 }
2aab9269 645 //if it is from a K0S
f61cec2f 646 if(primIdx < 0)
647 {
648 std::cout << "What!? No primary?" << std::endl;
649 continue;
650 }
651
e2ee5727 652 } // end of primary particle check
b2c10007 653
f61cec2f 654 const int primCode = stack->Particle(primIdx)->GetPdgCode();
ef647350 655 TParticlePDG *pdg = part->GetPDG();
f61cec2f 656 if (!pdg)
ef647350 657 {
658 Printf("ERROR: Could not get particle PDG %d", iPart);
659 continue;
f61cec2f 660 }
e2ee5727 661
f61cec2f 662 Int_t code = pdg->PdgCode();
b2c10007 663 if(primCode == fgGammaCode)
664 {
665
666
86e7d5db 667 if(fSelector->PassDistanceToBadChannelCut(*caloCluster))//&&fSelector->CutGeometricalAcceptance(*(stack->Particle(primIdx))))
b2c10007 668 {
669// std::cout << "Gamma primary: " << primIdx << std::endl;
670 foundGammas.push_back(primIdx);
671 }
672
673 }
674 fCutFlow->Fill(cf++);
86e7d5db 675 if(!fSelector->PassDistanceToBadChannelCut(*caloCluster)) continue;
b2c10007 676 Double_t clEt = CalculateTransverseEnergy(*caloCluster);
ef647350 677// if(code == fgK0SCode) std::cout << "K0 energy: " << caloCluster->E() << std::endl;
86e7d5db 678 if(!fSelector->PassMinEnergyCut(*caloCluster)) continue;
f61cec2f 679 fCutFlow->Fill(cf++);
e2ee5727 680 Float_t pos[3];
b2c10007 681 //PrintFamilyTree(
e2ee5727 682 caloCluster->GetPosition(pos);
683 TVector3 cp(pos);
684
f61cec2f 685 TParticle *primPart = stack->Particle(primIdx);
686 fPrimaryCode = primPart->GetPdgCode();
03bfd268 687 fPrimaryCharge = (Int_t) primPart->GetPDG()->Charge();
f61cec2f 688
689 fPrimaryE = primPart->Energy();
690 fPrimaryEt = primPart->Energy()*TMath::Sin(primPart->Theta());
691
692 fPrimaryPx = primPart->Px();
693 fPrimaryPy = primPart->Py();
694 fPrimaryPz = primPart->Pz();
695
696 fPrimaryVx = primPart->Vx();
697 fPrimaryVy = primPart->Vy();
698 fPrimaryVz = primPart->Vz();
699
700 fPrimaryAccepted = false;
b2c10007 701 fPrimaryMatched = false;
f61cec2f 702
703 fDepositedCode = part->GetPdgCode();
b2c10007 704 fDepositedE = part->Energy();
705 fDepositedEt = part->Energy()*TMath::Sin(part->Theta());
03bfd268 706 fDepositedCharge = (Int_t) part->GetPDG()->Charge();
f61cec2f 707
708 fDepositedVx = part->Vx();
709 fDepositedVy = part->Vy();
710 fDepositedVz = part->Vz();
b2c10007 711
712 //fSecondary = fSelector->FromSecondaryInteraction(*primPart, *stack);
713 fSecondary =fSelector->FromSecondaryInteraction(*part, *stack);
714 if(fSecondary)
715 {
716 //std::cout << "Have secondary!" << std::endl;
717 //PrintFamilyTree(iPart, stack);
718 }
719 fReconstructedE = caloCluster->E();
720 fReconstructedEt = caloCluster->E()*TMath::Sin(cp.Theta());
e2ee5727 721 //if(caloCluster->GetEmcCpvDistance() < fTrackDistanceCut)
b2c10007 722
723 pdg = primPart->GetPDG(0);
724 code = primPart->GetPdgCode();
725 //if (TMath::Abs(caloCluster->GetTrackDx()) < 5 && TMath::Abs(caloCluster->GetTrackDz()) < 5)
86e7d5db 726 if(!fSelector->PassTrackMatchingCut(*caloCluster))
e2ee5727 727 {
b2c10007 728 fPrimaryMatched = true;
e2ee5727 729 if (pdg->Charge() != 0)
730 {
e2ee5727 731 fChargedRemoved++;
f61cec2f 732 fEnergyChargedRemoved += clEt;
e2ee5727 733 fHistRemovedOrNot->Fill(0.0, fCentClass);
e2ee5727 734 fHistDxDzRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
b2c10007 735
e2ee5727 736 }
737 else
738 {
b2c10007 739 if(code==fgGammaCode)
740 {
741 fGammaRemoved++;
742 fGammaRemovedEt+=clEt;
743
744 }
745 else
746 {
747 fNeutralRemoved++;
748 fEnergyNeutralRemoved += clEt;
749 fHistRemovedOrNot->Fill(1.0, fCentClass);
750 fHistDxDzRemovedNeutral->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
751 }
e2ee5727 752 }
753 }
754 else
755 {
b2c10007 756 fPrimaryAccepted = true;
757 fPrimaryMatched = false;
f61cec2f 758 fCutFlow->Fill(cf++);
e2ee5727 759 if (pdg->Charge() != 0)
760 {
e2ee5727 761
b2c10007 762 if(fSecondary)
763 {
764 fSecondaryNotRemoved++;
765 }
766 else
767 {
768 fChargedNotRemoved++;
769 fEnergyChargedNotRemoved += clEt;
770 fHistRemovedOrNot->Fill(2.0, fCentClass);
771 fHistDxDzNonRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
2aab9269 772 //PrintFamilyTree(iPart, stack);
b2c10007 773 }
e2ee5727 774 }
775 else
776 {
b2c10007 777 if(!fSecondary)
778 {
779 fNeutralNotRemoved++;
780 fEnergyNeutralNotRemoved += clEt;
781 fHistRemovedOrNot->Fill(3.0, fCentClass);
782 fHistDxDzNonRemovedNeutral->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
783 }
784 else
785 {
2aab9269 786// //PrintFamilyTree(iPart, stack);
b2c10007 787 }
788 // if(TMath::Abs(part->Vx()) < 1.0 && TMath::Abs(part->Vy()) < 1.0 && TMath::Abs(part->Vz()) < 20 && fSelector->IsEmEtParticle(primCode))
789 if(fSecondary && fSelector->IsEmEtParticle(primCode))
f61cec2f 790 {
b2c10007 791 fTotEtSecondaryFromEmEtPrimary += clEt;
792 fSecondaryNotRemoved++;
f61cec2f 793 }
794 else if(fSelector->IsEmEtParticle(primCode))
795 {
b2c10007 796 fTotEtWithSecondaryRemoved += clEt;
797// PrintFamilyTree(iPart, stack);
798
f61cec2f 799 }
b2c10007 800 else
801 {
802 //fTotEtSecondary += clEt;
803// PrintFamilyTree(iPart,stack);
f61cec2f 804 }
805 //code = stack->Particle(primIdx)->GetPdgCode();
b2c10007 806 if (code == fgGammaCode && !fSecondary)
e2ee5727 807 {
808 fEtNonRemovedGammas += clEt;
809 fMultNonRemovedGammas++;
b2c10007 810 fNeutralNotRemoved--;
811 fEnergyNeutralNotRemoved -= clEt;
812 //PrintFamilyTree(iPart, stack);
ef647350 813// if (pdgMom)
f61cec2f 814 // {
815 // if (TMath::Abs(pdgMom->PdgCode()) == fgPi0Code || TMath::Abs(pdgMom->PdgCode()) == fgEtaCode || TMath::Abs(pdgMom->PdgCode()) == 331)
816 // {
e2ee5727 817// std::cout << "Mother of gamma: " << pdgMom->PdgCode() << " " << pdgMom->GetName() << ", E: " << part->Energy() << std::endl;
f61cec2f 818 // fEtNonRemovedGammasFromPi0 += clEt;
819 // }
820 // }
e2ee5727 821 }
e2ee5727 822 }
823 }
f61cec2f 824 fPrimaryTree->Fill();
e2ee5727 825 } // end of loop over clusters
826
b2c10007 827
828 std::sort(foundGammas.begin(), foundGammas.end());
829 for (Int_t iPart = 0; iPart < stack->GetNtrack(); iPart++)
830 {
831
832 if(!stack->IsPhysicalPrimary(iPart)) continue;
833
834 TParticle *part = stack->Particle(iPart);
835
836 if (!part)
837 {
838 Printf("ERROR: Could not get particle %d", iPart);
839 continue;
840 }
841 TParticlePDG *pdg = part->GetPDG(0);
842
843 if (!pdg)
844 {
845 Printf("ERROR: Could not get particle PDG %d", iPart);
846 continue;
847 }
848
849 if(pdg->PdgCode()==fgGammaCode && fSelector->CutGeometricalAcceptance(*part))// TMath::Abs(part->Eta()) < 0.12)
850 {
851 fHistGammasGenerated->Fill(part->Energy());
852// std::cout << "Searching for gammas..." << std::endl;
853 if(std::binary_search(foundGammas.begin(),foundGammas.end(),iPart))
854 {
855 fHistGammasFound->Fill(part->Energy());
856 }
857// for(int i = 0; i < foundGammas.size(); i++)
858// {
859// // std::cout << iPart << std::endl;
860// if(foundGammas[i] == iPart)
861// {
862// fHistGammasFound->Fill(part->Energy());
863// std::cout << "Match!" << std::endl;
864// break;
865// }
866// }
867 }
868
869 }
e2ee5727 870 //std::cout << "Distance cut: " << fTrackDistanceCut << std::endl;
871 //std::cout << "Number of removed neutrals: " << fNeutralRemoved << std::endl;
872 //std::cout << "Number of removed charged: " << fChargedRemoved << std::endl;
873 //std::cout << "Number of non-removed charged: " << fChargedNotRemoved << std::endl;
874 //std::cout << "Number of non-removed neutral: " << fNeutralNotRemoved << std::endl;
875
876// std::cout << "Energy of removed neutrals: " << fEnergyNeutralRemoved << std::endl;
877// std::cout << "Energy of removed charged: " << fEnergyChargedRemoved << std::endl;
878// std::cout << "Energy of non-removed charged: " << fEnergyChargedNotRemoved << std::endl;
879// std::cout << "Energy of non-removed neutral: " << fEnergyNeutralNotRemoved << std::endl;
2aab9269 880 Float_t etCuts[10] = {0.0,0.1,0.15,0.2,0.25, 0.3,0.35,0.4,0.45,0.5};
881 Int_t nEtCuts = 10;
882 //loop over simulated particles in order to find K0S
883 for (Int_t iPart = 0; iPart < stack->GetNtrack(); iPart++){
884 TParticle *part = stack->Particle(iPart);
885 if (!part){
886 //Printf("ERROR: Could not get particle %d", iPart);
887 continue;
888 }
889 TParticlePDG *pdg = part->GetPDG(0);
890 if (!pdg){
891 //Printf("ERROR: Could not get particle PDG %d", iPart);
892 continue;
893 }
894
895 if(stack->IsPhysicalPrimary(iPart)){//if it is a K0 it might have decayed into four pions
896 //fgK0SCode, fgGammaCode, fgPi0Code
897 Int_t code = pdg->PdgCode();
898 if(code == fgK0SCode){//this is a K0S
899 if(stack->Particle(iPart)->GetNDaughters()==2){//it has two daughters
900 Int_t daughterID = stack->Particle(iPart)->GetFirstDaughter();
901 if(stack->Particle(daughterID)->GetPDG(0)->PdgCode()==fgPi0Code){//and it decayed into a pi0
902 //cout<<"I am a gamma from a pi0 from a K0S"<<endl;
903 Float_t pTk0 = stack->Particle(iPart)->Pt();
904 Int_t gammaDaughterIDs[6] = {0,0,0,0,0,0};
905 Float_t gammaEts[4] = {0.0,0.0,0.0,0.0};
906 Float_t totalGammaEts[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
907 Float_t totalClusterEts[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
908 gammaDaughterIDs[4] = stack->Particle(iPart)->GetDaughter(0);
909 gammaDaughterIDs[5] = stack->Particle(iPart)->GetDaughter(1);
910 gammaDaughterIDs[0] = stack->Particle(gammaDaughterIDs[4])->GetDaughter(0);
911 gammaDaughterIDs[1] = stack->Particle(gammaDaughterIDs[4])->GetDaughter(1);
912 gammaDaughterIDs[2] = stack->Particle(gammaDaughterIDs[5])->GetDaughter(0);
913 gammaDaughterIDs[3] = stack->Particle(gammaDaughterIDs[5])->GetDaughter(1);
914 for(int k=0;k<4;k++){
915 gammaEts[k] = stack->Particle(gammaDaughterIDs[4])->Energy() * TMath::Sin(stack->Particle(gammaDaughterIDs[4])->Theta() );
916 }
917 //does not always have all of the daughters
918 if(gammaDaughterIDs[0]==gammaDaughterIDs[1]){
919 gammaDaughterIDs[1] = -1;
920 gammaEts[1] = 0.0;
921 //cout<<"Duplicates. This has "<<stack->Particle(gammaDaughterIDs[4])->GetNDaughters()<<" daughters"<<endl;
922 }
923 if(gammaDaughterIDs[2]==gammaDaughterIDs[3]){
924 gammaDaughterIDs[3] = -1;
925 gammaEts[3] = 0.0;
926 //cout<<"Duplicates. This has "<<stack->Particle(gammaDaughterIDs[5])->GetNDaughters()<<" daughters"<<endl;
927 }
928// cout<<"My daughters are ";
929// for(int j=0;j<6;j++){
930// cout <<gammaDaughterIDs[j]<<" ";
931// if(gammaDaughterIDs[j]>=0) cout<<stack->Particle(gammaDaughterIDs[j])->GetName();
932// cout<<",";
933// }
934// cout<<endl;
935 for(int l=0;l<nEtCuts;l++){//loop over cut values
936 for(int k=0;k<4;k++){//loop over gamma daughter energies
937 if(gammaEts[k]>=etCuts[l]){
938 totalGammaEts[l] += gammaEts[k];//if gamma et is above the cut off energy add it
939 }
940 }
941 }
942
2aab9269 943 for (int iCluster = 0; iCluster < nCluster; iCluster++ ){//if this cluster is from any of the decay daughters of our K0S
944 AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
945 const UInt_t myPart = (UInt_t)TMath::Abs(caloCluster->GetLabel());
946 for(int j=0;j<6;j++){
03bfd268 947 if( myPart==((UInt_t) gammaDaughterIDs[j]) ){
c31562f7 948 //cout<<"Found a matching cluster!";//<<endl;
2aab9269 949 Double_t clEt = CalculateTransverseEnergy(*caloCluster);
c31562f7 950 //cout<<" it has energy "<<clEt<<endl;
2aab9269 951 for(int l=0;l<nEtCuts;l++){//loop over cut values
952 if(clEt>=etCuts[l]){
953 totalClusterEts[l] += clEt;//if cluster et is above the cut off energy add it
954 }
955 }
956
957 }
958 }
959 }
960 //now we have the total energy measured for each case
961 fHistK0EDepositsVsPtForEOver000MeV->Fill(pTk0,totalClusterEts[0]);
962 fHistK0EDepositsVsPtForEOver100MeV->Fill(pTk0,totalClusterEts[1]);
963 fHistK0EDepositsVsPtForEOver150MeV->Fill(pTk0,totalClusterEts[2]);
964 fHistK0EDepositsVsPtForEOver200MeV->Fill(pTk0,totalClusterEts[3]);
965 fHistK0EDepositsVsPtForEOver250MeV->Fill(pTk0,totalClusterEts[4]);
966 fHistK0EDepositsVsPtForEOver300MeV->Fill(pTk0,totalClusterEts[5]);
967 fHistK0EDepositsVsPtForEOver350MeV->Fill(pTk0,totalClusterEts[6]);
968 fHistK0EDepositsVsPtForEOver400MeV->Fill(pTk0,totalClusterEts[7]);
969 fHistK0EDepositsVsPtForEOver450MeV->Fill(pTk0,totalClusterEts[8]);
970 fHistK0EDepositsVsPtForEOver500MeV->Fill(pTk0,totalClusterEts[9]);
971 fHistK0EGammaVsPtForEOver000MeV->Fill(pTk0,totalGammaEts[0]);
972 fHistK0EGammaVsPtForEOver100MeV->Fill(pTk0,totalGammaEts[1]);
973 fHistK0EGammaVsPtForEOver150MeV->Fill(pTk0,totalGammaEts[2]);
974 fHistK0EGammaVsPtForEOver200MeV->Fill(pTk0,totalGammaEts[3]);
975 fHistK0EGammaVsPtForEOver250MeV->Fill(pTk0,totalGammaEts[4]);
976 fHistK0EGammaVsPtForEOver300MeV->Fill(pTk0,totalGammaEts[5]);
977 fHistK0EGammaVsPtForEOver350MeV->Fill(pTk0,totalGammaEts[6]);
978 fHistK0EGammaVsPtForEOver400MeV->Fill(pTk0,totalGammaEts[7]);
979 fHistK0EGammaVsPtForEOver450MeV->Fill(pTk0,totalGammaEts[8]);
980 fHistK0EGammaVsPtForEOver500MeV->Fill(pTk0,totalGammaEts[9]);
981 }
982 }
983 }
984 }
985 }
986
e2ee5727 987 FillHistograms();
988 return 0;
989}
2fbf38ac 990
991void AliAnalysisEtMonteCarlo::Init()
f61cec2f 992{ // init
e2ee5727 993 AliAnalysisEt::Init();
2fbf38ac 994}
13b0d3c1 995
ce546038 996void AliAnalysisEtMonteCarlo::ResetEventValues()
f61cec2f 997{ // reset event values
e2ee5727 998 AliAnalysisEt::ResetEventValues();
999
f61cec2f 1000 fTotEtSecondary = 0;
1001 fTotEtSecondaryFromEmEtPrimary = 0;
1002 fTotEtWithSecondaryRemoved = 0;
1003
e2ee5727 1004 // collision geometry defaults for p+p:
1005 fImpactParameter = 0;
1006 fNcoll = 1;
1007 fNpart = 2;
1008
1009 fEtNonRemovedProtons = 0;
1010 fEtNonRemovedAntiProtons = 0;
1011 fEtNonRemovedPiPlus = 0;
1012 fEtNonRemovedPiMinus = 0;
1013 fEtNonRemovedKaonPlus = 0;
1014 fEtNonRemovedKaonMinus = 0;
ef647350 1015 fEtNonRemovedK0S = 0;
1016 fEtNonRemovedK0L = 0;
e2ee5727 1017 fEtNonRemovedLambdas = 0;
1018 fEtNonRemovedElectrons = 0;
1019 fEtNonRemovedPositrons = 0;
1020 fEtNonRemovedMuPlus = 0;
1021 fEtNonRemovedMuMinus = 0;
1022 fEtNonRemovedNeutrons = 0;
1023 fEtNonRemovedAntiNeutrons = 0;
1024 fEtNonRemovedGammas = 0;
1025 fEtNonRemovedGammasFromPi0 = 0;
1026
ef647350 1027 fEtRemovedProtons = 0;
1028 fEtRemovedAntiProtons = 0;
1029 fEtRemovedPiPlus = 0;
1030 fEtRemovedPiMinus = 0;
1031 fEtRemovedKaonPlus = 0;
1032 fEtRemovedKaonMinus = 0;
1033 fEtRemovedK0s = 0;
1034 fEtRemovedK0L = 0;
1035 fEtRemovedLambdas = 0;
1036 fEtRemovedElectrons = 0;
1037 fEtRemovedPositrons = 0;
1038 fEtRemovedMuPlus = 0;
1039 fEtRemovedMuMinus = 0;
1040 fEtRemovedNeutrons = 0;
1041
1042 fEtRemovedGammasFromPi0 = 0;
e2ee5727 1043 fEtRemovedGammas = 0;
1044 fEtRemovedNeutrons = 0;
1045 fEtRemovedAntiNeutrons = 0;
1046
1047 fMultNonRemovedProtons = 0;
1048 fMultNonRemovedAntiProtons = 0;
1049 fMultNonRemovedPiPlus = 0;
1050 fMultNonRemovedPiMinus = 0;
1051 fMultNonRemovedKaonPlus = 0;
1052 fMultNonRemovedKaonMinus = 0;
1053 fMultNonRemovedK0s = 0;
ef647350 1054 fMultNonRemovedK0L = 0;
e2ee5727 1055 fMultNonRemovedLambdas = 0;
1056 fMultNonRemovedElectrons = 0;
1057 fMultNonRemovedPositrons = 0;
1058 fMultNonRemovedMuPlus = 0;
1059 fMultNonRemovedMuMinus = 0;
1060 fMultNonRemovedNeutrons = 0;
1061 fMultNonRemovedAntiNeutrons = 0;
1062 fMultNonRemovedGammas = 0;
1063
ef647350 1064 fMultRemovedProtons = 0;
1065 fMultRemovedAntiProtons = 0;
1066 fMultRemovedPiPlus = 0;
1067 fMultRemovedPiMinus = 0;
1068 fMultRemovedKaonPlus = 0;
1069 fMultRemovedKaonMinus = 0;
1070 fMultRemovedK0s = 0;
1071 fMultRemovedK0L = 0;
1072 fMultRemovedLambdas = 0;
1073 fMultRemovedElectrons = 0;
1074 fMultRemovedPositrons = 0;
1075 fMultRemovedMuPlus = 0;
1076 fMultRemovedMuMinus = 0;
1077
e2ee5727 1078 fMultRemovedGammas = 0;
1079 fMultRemovedNeutrons = 0;
1080 fMultRemovedAntiNeutrons = 0;
1081
ef647350 1082 fEnergyChargedNotRemoved = 0;
1083 fEnergyChargedRemoved = 0;
1084 fEnergyNeutralNotRemoved = 0;
1085 fEnergyNeutralRemoved = 0;
f61cec2f 1086
ef647350 1087 fChargedNotRemoved = 0;
1088 fChargedRemoved = 0;
1089 fNeutralNotRemoved = 0;
1090 fNeutralRemoved = 0;
b2c10007 1091 fGammaRemoved = 0;
1092 fSecondaryNotRemoved = 0;
f61cec2f 1093
e2ee5727 1094 fTrackMultInAcc = 0;
1095
f61cec2f 1096 fTotNeutralEtAfterMinEnergyCut = 0;
b2c10007 1097
1098 fSecondaryNotRemoved = 0;
1099
1100 fTotPx = 0;
1101 fTotPy = 0;
1102 fTotPz = 0;
1103
1104
ce546038 1105}
1106
1107void AliAnalysisEtMonteCarlo::CreateHistograms()
f61cec2f 1108{ // histogram related additions
e2ee5727 1109 AliAnalysisEt::CreateHistograms();
f61cec2f 1110 if (fEventSummaryTree) {
1111 fEventSummaryTree->Branch("fImpactParameter",&fImpactParameter,"fImpactParameter/D");
1112 fEventSummaryTree->Branch("fNcoll",&fNcoll,"fNcoll/I");
1113 fEventSummaryTree->Branch("fNpart",&fNpart,"fNpart/I");
1114 fEventSummaryTree->Branch("fTotEtWithSecondaryRemoved", &fTotEtWithSecondaryRemoved, "fTotEtWithSecondaryRemoved/D");
1115 fEventSummaryTree->Branch("fTotEtSecondaryFromEmEtPrimary", &fTotEtSecondaryFromEmEtPrimary, "fTotEtSecondaryFromEmEtPrimary/D");
1116 fEventSummaryTree->Branch("fTotEtSecondary", &fTotEtSecondary, "fTotEtSecondary/D");
1117 fEventSummaryTree->Branch("fTotNeutralEtAfterMinEnergyCut", &fTotNeutralEtAfterMinEnergyCut, "fTotNeutralEtAfterMinEnergyCut/D");
b2c10007 1118 fEventSummaryTree->Branch("fSecondaryNotRemoved", &fSecondaryNotRemoved, "fSecondaryNotRemoved/I");
1119 fEventSummaryTree->Branch("fChargedNotRemoved", &fChargedNotRemoved, "fChargedNotRemoved/I");
1120 fEventSummaryTree->Branch("fNeutralNotRemoved", &fNeutralNotRemoved, "fNeutralNotRemoved/I");
1121 fEventSummaryTree->Branch("fChargedRemoved", &fChargedRemoved, "fChargedRemoved/I");
1122 fEventSummaryTree->Branch("fNeutralRemoved", &fNeutralRemoved, "fNeutralRemoved/I");
1123 fEventSummaryTree->Branch("fGammaRemoved", &fGammaRemoved, "fGammaRemoved/I");
1124 fEventSummaryTree->Branch("fTotPx", &fTotPx, "fTotPx/D");
1125 fEventSummaryTree->Branch("fTotPy", &fTotPy, "fTotPy/D");
1126 fEventSummaryTree->Branch("fTotPz", &fTotPz, "fTotPz/D");
1127// fEventSummaryTree->Branch("f
e2ee5727 1128 }
f61cec2f 1129
e2ee5727 1130 //fHistDecayVertexNonRemovedCharged = new TH3F("fHistDecayVertexNonRemovedCharged","fHistDecayVertexNonRemovedCharged", 500, -470, 30, 500, -300, 300, 40, -20, 20);
1131 //fHistDecayVertexRemovedCharged = new TH3F("fHistDecayVertexRemovedCharged","fHistDecayVertexRemovedCharged", 500, -470, 30, 500, -300, 300, 40, -20, 20);
1132 //fHistDecayVertexNonRemovedNeutral = new TH3F("fHistDecayVertexNonRemovedNeutral","fHistDecayVertexNonRemovedNeutral", 500, -470, 30, 500, -300, 300, 40, -20, 20);
1133 //fHistDecayVertexRemovedNeutral = new TH3F("fHistDecayVertexRemovedNeutral","fHistDecayVertexRemovedNeutral", 500, -470, 30, 500, -300, 300, 40, -20, 20);
1134
1135 fHistRemovedOrNot = new TH2F("fHistRemovedOrNot", "fHistRemovedOrNot", 4, -0.5, 3.5, 10, -0.5, 9.5);
1136
1137 fHistEtNonRemovedProtons = new TH2F("fHistEtNonRemovedProtons", "fHistEtNonRemovedProtons", 1500, 0, 30, 10, -0.5, 9.5);
1138 fHistEtNonRemovedAntiProtons = new TH2F("fHistEtNonRemovedAntiProtons", "fHistEtNonRemovedAntiProtons", 1500, 0, 30, 10, -0.5, 9.5);
1139 fHistEtNonRemovedPiPlus = new TH2F("fHistEtNonRemovedPiPlus", "fHistEtNonRemovedPiPlus", 1500, 0, 30, 10, -0.5, 9.5);
1140 fHistEtNonRemovedPiMinus = new TH2F("fHistEtNonRemovedPiMinus", "fHistEtNonRemovedPiMinus", 1500, 0, 30, 10, -0.5, 9.5);
1141 fHistEtNonRemovedKaonPlus = new TH2F("fHistEtNonRemovedKaonPlus", "fHistEtNonRemovedKaonPlus", 1500, 0, 30, 10, -0.5, 9.5);
1142 fHistEtNonRemovedKaonMinus = new TH2F("fHistEtNonRemovedKaonMinus", "fHistEtNonRemovedKaonMinus", 1500, 0, 30, 10, -0.5, 9.5);
1143 fHistEtNonRemovedK0s = new TH2F("fHistEtNonRemovedK0s", "fHistEtNonRemovedK0s", 1500, 0, 30, 10, -0.5, 9.5);
ef647350 1144 fHistEtNonRemovedK0L = new TH2F("fHistEtNonRemovedK0L", "fHistEtNonRemovedK0L", 1500, 0, 30, 10, -0.5, 9.5);
e2ee5727 1145 fHistEtNonRemovedLambdas = new TH2F("fHistEtNonRemovedLambdas", "fHistEtNonRemovedLambdas", 1500, 0, 30, 10, -0.5, 9.5);
1146 fHistEtNonRemovedElectrons = new TH2F("fHistEtNonRemovedElectrons", "fHistEtNonRemovedElectrons", 1500, 0, 30, 10, -0.5, 9.5);
1147 fHistEtNonRemovedPositrons = new TH2F("fHistEtNonRemovedPositrons", "fHistEtNonRemovedPositrons", 1500, 0, 30, 10, -0.5, 9.5);
1148 fHistEtNonRemovedMuPlus = new TH2F("fHistEtNonRemovedMuPlus", "fHistEtNonRemovedMuPlus", 1500, 0, 30, 10, -0.5, 9.5);
1149 fHistEtNonRemovedMuMinus = new TH2F("fHistEtNonRemovedMuMinus", "fHistEtNonRemovedMuMinus", 1500, 0, 30, 10, -0.5, 9.5);
1150 fHistEtNonRemovedNeutrons = new TH2F("fHistEtNonRemovedNeutrons", "fHistEtNonRemovedNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
1151 fHistEtNonRemovedAntiNeutrons = new TH2F("fHistEtNonRemovedAntiNeutrons", "fHistEtNonRemovedAntiNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
1152
1153 fHistEtNonRemovedGammas = new TH2F("fHistEtNonRemovedGammas", "fHistEtNonRemovedGammas", 1500, 0, 30, 10, -0.5, 9.5);
1154 fHistEtNonRemovedGammasFromPi0 = new TH2F("fHistEtNonRemovedGammasFromPi0", "fHistEtNonRemovedGammasFromPi0", 1500, 0, 30, 10, -0.5, 9.5);
1155
1156 fHistEtRemovedGammas = new TH2F("fHistEtRemovedGammas", "fHistEtRemovedGammas", 1500, 0, 30, 10, -0.5, 9.5);
1157 fHistEtRemovedNeutrons = new TH2F("fHistEtRemovedNeutrons", "fHistEtRemovedNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
1158 fHistEtRemovedAntiNeutrons = new TH2F("fHistEtRemovedAntiNeutrons", "fHistEtRemovedAntiNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
f61cec2f 1159
ef647350 1160 fHistEtRemovedCharged = new TH2F("fHistEtRemovedCharged", "fHistEtRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
1161 fHistEtRemovedNeutrals = new TH2F("fHistEtRemovedNeutrals", "fHistEtRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
e2ee5727 1162
ef647350 1163 fHistEtNonRemovedCharged = new TH2F("fHistEtNonRemovedCharged", "fHistEtNonRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
1164 fHistEtNonRemovedNeutrals = new TH2F("fHistEtNonRemovedNeutrals", "fHistEtNonRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
f61cec2f 1165
e2ee5727 1166 fHistMultNonRemovedProtons = new TH2F("fHistMultNonRemovedProtons", "fHistMultNonRemovedProtons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1167 fHistMultNonRemovedAntiProtons = new TH2F("fHistMultNonRemovedAntiProtons", "fHistMultNonRemovedAntiProtons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1168 fHistMultNonRemovedPiPlus = new TH2F("fHistMultNonRemovedPiPlus", "fHistMultNonRemovedPiPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1169 fHistMultNonRemovedPiMinus = new TH2F("fHistMultNonRemovedPiMinus", "fHistMultNonRemovedPiMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1170 fHistMultNonRemovedKaonPlus = new TH2F("fHistMultNonRemovedKaonPlus", "fHistMultNonRemovedKaonPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1171 fHistMultNonRemovedKaonMinus = new TH2F("fHistMultNonRemovedKaonMinus", "fHistMultNonRemovedKaonMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1172 fHistMultNonRemovedK0s = new TH2F("fHistMultNonRemovedK0s", "fHistMultNonRemovedK0s", 100, -0.5, 99.5, 10, -0.5, 9.5);
ef647350 1173 fHistMultNonRemovedK0L = new TH2F("fHistMultNonRemovedK0L", "fHistMultNonRemovedK0L", 100, -0.5, 99.5, 10, -0.5, 9.5);
e2ee5727 1174 fHistMultNonRemovedLambdas = new TH2F("fHistMultNonRemovedLambdas", "fHistMultNonRemovedLambdas", 100, -0.5, 99.5, 10, -0.5, 9.5);
1175 fHistMultNonRemovedElectrons = new TH2F("fHistMultNonRemovedElectrons", "fHistMultNonRemovedElectrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1176 fHistMultNonRemovedPositrons = new TH2F("fHistMultNonRemovedPositrons", "fHistMultNonRemovedPositrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1177 fHistMultNonRemovedMuPlus = new TH2F("fHistMultNonRemovedMuPlus", "fHistMultNonRemovedMuPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1178 fHistMultNonRemovedMuMinus = new TH2F("fHistMultNonRemovedMuMinus", "fHistMultNonRemovedMuMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1179 fHistMultNonRemovedNeutrons = new TH2F("fHistMultNonRemovedNeutrons", "fHistMultNonRemovedNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1180 fHistMultNonRemovedAntiNeutrons = new TH2F("fHistMultNonRemovedAntiNeutrons", "fHistMultNonRemovedAntiNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1181
ef647350 1182 fHistMultNonRemovedGammas = new TH2F("fHistMultNonRemovedGammas", "fHistMultNonRemovedGammas", 100, -0.5, 99.5, 100, -0.5, 99.5);
e2ee5727 1183
ef647350 1184 fHistMultRemovedGammas = new TH2F("fHistMultRemovedGammas", "fHistMultRemovedGammas", 100, -0.5, 99.5, 100, -0.5, 99.5);
e2ee5727 1185 fHistMultRemovedNeutrons = new TH2F("fHistMultRemovedNeutrons", "fHistMultRemovedNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1186 fHistMultRemovedAntiNeutrons = new TH2F("fHistMultRemovedAntiNeutrons", "fHistMultRemovedAntiNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
f61cec2f 1187 /*
1188 fHistMultRemovedCharged = new TH2F("fHistMultRemovedCharged", "fHistMultRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
1189 fHistMultRemovedNeutrals = new TH2F("fHistMultRemovedNeutrals", "fHistMultRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
1190
1191 fHistMultNonRemovedCharged = new TH2F("fHistMultNonRemovedCharged", "fHistMultNonRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
1192 fHistMultNonRemovedNeutrals = new TH2F("fHistMultNonRemovedNeutrals", "fHistMultNonRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);*/
ef647350 1193
ef647350 1194
ef647350 1195 fHistMultRemovedCharged = new TH2F("fHistMultRemovedCharged", "fHistMultRemovedCharged", 100, -0.5, 99.5, 100, -0.5, 99.5);
1196 fHistMultRemovedNeutrals = new TH2F("fHistMultRemovedNeutrals", "fHistMultRemovedNeutrals", 100, -0.5, 99.5, 100, -0.5, 99.5);
1197
1198 fHistMultNonRemovedCharged = new TH2F("fHistMultNonRemovedCharged", "fHistMultNonRemovedCharged", 100, -0.5, 99.5, 100, -0.5, 99.5);
1199 fHistMultNonRemovedNeutrals = new TH2F("fHistMultNonRemovedNeutrals", "fHistMultNonRemovedNeutrals", 100, -0.5, 99.5, 100, -0.5, 99.5);
e2ee5727 1200
1201 fHistTrackMultvsNonRemovedCharged = new TH2F("fHistTrackMultvsNonRemovedCharged", "fHistTrackMultvsNonRemovedCharged", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1202 fHistTrackMultvsNonRemovedNeutral = new TH2F("fHistTrackMultvsNonRemovedNeutral", "fHistTrackMultvsNonRemovedNeutral", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1203 fHistTrackMultvsRemovedGamma = new TH2F("fHistTrackMultvsRemovedGamma", "fHistTrackMultvsRemovedGamma", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1204
1205 fHistClusterMultvsNonRemovedCharged = new TH2F("fHistClusterMultvsNonRemovedCharged", "fHistClusterMultvsNonRemovedCharged", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1206 fHistClusterMultvsNonRemovedNeutral = new TH2F("fHistClusterMultvsNonRemovedNeutral", "fHistClusterMultvsNonRemovedNeutral", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1207 fHistClusterMultvsRemovedGamma = new TH2F("fHistClusterMultvsRemovedGamma", "fHistClusterMultvsRemovedGamma", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1208
1209 fHistDxDzNonRemovedCharged = new TH2F("fHistDxDzNonRemovedCharged", "fHistDxDzNonRemovedCharged", 800, -200, 200, 800, -200, 200);
1210 fHistDxDzRemovedCharged = new TH2F("fHistDxDzRemovedCharged", "fHistDxDzRemovedCharged", 800, -200, 200, 800, -200, 200);
1211 fHistDxDzNonRemovedNeutral = new TH2F("fHistDxDzNonRemovedNeutral", "fHistDxDzNonRemovedNeutral", 800, -200, 200, 800, -200, 200);
1212 fHistDxDzRemovedNeutral = new TH2F("fHistDxDzRemovedNeutral", "fHistDxDzRemovedNeutral", 800, -200, 200, 800, -200, 200);
1213
2aab9269 1214 fHistK0EDepositsVsPtForEOver000MeV = new TH2F("fHistK0EDepositsVsPtForEOver000MeV","K^{0}_{S} deposits over 0 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1215 fHistK0EDepositsVsPtForEOver100MeV = new TH2F("fHistK0EDepositsVsPtForEOver100MeV","K^{0}_{S} deposits over 100 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1216 fHistK0EDepositsVsPtForEOver150MeV = new TH2F("fHistK0EDepositsVsPtForEOver150MeV","K^{0}_{S} deposits over 150 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1217 fHistK0EDepositsVsPtForEOver200MeV = new TH2F("fHistK0EDepositsVsPtForEOver200MeV","K^{0}_{S} deposits over 200 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1218 fHistK0EDepositsVsPtForEOver250MeV = new TH2F("fHistK0EDepositsVsPtForEOver250MeV","K^{0}_{S} deposits over 250 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1219 fHistK0EDepositsVsPtForEOver300MeV = new TH2F("fHistK0EDepositsVsPtForEOver300MeV","K^{0}_{S} deposits over 300 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1220 fHistK0EDepositsVsPtForEOver350MeV = new TH2F("fHistK0EDepositsVsPtForEOver350MeV","K^{0}_{S} deposits over 350 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1221 fHistK0EDepositsVsPtForEOver400MeV = new TH2F("fHistK0EDepositsVsPtForEOver400MeV","K^{0}_{S} deposits over 400 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1222 fHistK0EDepositsVsPtForEOver450MeV = new TH2F("fHistK0EDepositsVsPtForEOver450MeV","K^{0}_{S} deposits over 450 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1223 fHistK0EDepositsVsPtForEOver500MeV = new TH2F("fHistK0EDepositsVsPtForEOver500MeV","K^{0}_{S} deposits over 500 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1224 fHistK0EGammaVsPtForEOver000MeV = new TH2F("fHistK0EGammaVsPtForEOver000MeV","K^{0}_{S} gamma over 0 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1225 fHistK0EGammaVsPtForEOver100MeV = new TH2F("fHistK0EGammaVsPtForEOver100MeV","K^{0}_{S} gamma over 100 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1226 fHistK0EGammaVsPtForEOver150MeV = new TH2F("fHistK0EGammaVsPtForEOver150MeV","K^{0}_{S} gamma over 150 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1227 fHistK0EGammaVsPtForEOver200MeV = new TH2F("fHistK0EGammaVsPtForEOver200MeV","K^{0}_{S} gamma over 200 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1228 fHistK0EGammaVsPtForEOver250MeV = new TH2F("fHistK0EGammaVsPtForEOver250MeV","K^{0}_{S} gamma over 250 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1229 fHistK0EGammaVsPtForEOver300MeV = new TH2F("fHistK0EGammaVsPtForEOver300MeV","K^{0}_{S} gamma over 300 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1230 fHistK0EGammaVsPtForEOver350MeV = new TH2F("fHistK0EGammaVsPtForEOver350MeV","K^{0}_{S} gamma over 350 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1231 fHistK0EGammaVsPtForEOver400MeV = new TH2F("fHistK0EGammaVsPtForEOver400MeV","K^{0}_{S} gamma over 400 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1232 fHistK0EGammaVsPtForEOver450MeV = new TH2F("fHistK0EGammaVsPtForEOver450MeV","K^{0}_{S} gamma over 450 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1233 fHistK0EGammaVsPtForEOver500MeV = new TH2F("fHistK0EGammaVsPtForEOver500MeV","K^{0}_{S} gamma over 500 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1234
e2ee5727 1235 fHistPiPlusMult = new TH1F("fHistPiPlusMult", "fHistPiPlusMult", 2000, -0.5, 1999.5);
1236 fHistPiMinusMult = new TH1F("fHistPiMinusMult", "fHistPiMinusMult", 2000, -0.5, 1999.5);
1237 fHistPiZeroMult = new TH1F("fHistPiZeroMult", "fHistPiZeroMult", 2000, -0.5, 1999.5);
1238
1239 fHistPiPlusMultAcc = new TH1F("fHistPiPlusMultAcc", "fHistPiPlusMultAcc", 2000, -0.5, 1999.5);
1240 fHistPiMinusMultAcc = new TH1F("fHistPiMinusMultAcc", "fHistPiMinusMultAcc", 2000, -0.5, 1999.5);
1241 fHistPiZeroMultAcc = new TH1F("fHistPiZeroMultAcc", "fHistPiZeroMultAcc", 2000, -0.5, 1999.5);
1242
f61cec2f 1243 if(fCuts->GetHistMakeTree())
1244 {
1245 TString treename = "fPrimaryTree" + fHistogramNameSuffix;
1246 fPrimaryTree = new TTree(treename, treename);
1247
1248 fPrimaryTree->Branch("fTotEt",&fTotEt,"fTotEt/D");
1249 fPrimaryTree->Branch("fNeutralMultiplicity", &fNeutralMultiplicity, "fNeutralMultiplicity/I");
1250 fPrimaryTree->Branch("fCentClass",&fCentClass,"fCentClass/I");
1251
1252 fPrimaryTree->Branch("fPrimaryCode", &fPrimaryCode, "fPrimaryCode/I");
1253 fPrimaryTree->Branch("fPrimaryCharge", &fPrimaryCharge, "fPrimaryCharge/I");
1254
1255 fPrimaryTree->Branch("fPrimaryE", &fPrimaryE, "fPrimaryE/D");
1256 fPrimaryTree->Branch("fPrimaryEt", &fPrimaryEt, "fPrimaryEt/D");
1257
1258 fPrimaryTree->Branch("fPrimaryPx", &fPrimaryPx, "fPrimaryPx/D");
1259 fPrimaryTree->Branch("fPrimaryPy", &fPrimaryPy, "fPrimaryPy/D");
1260 fPrimaryTree->Branch("fPrimaryPz", &fPrimaryPz, "fPrimaryPz/D");
1261
1262 fPrimaryTree->Branch("fPrimaryVx", &fPrimaryVx, "fPrimaryVx/D");
1263 fPrimaryTree->Branch("fPrimaryVy", &fPrimaryVy, "fPrimaryVy/D");
1264 fPrimaryTree->Branch("fPrimaryVz", &fPrimaryVz, "fPrimaryVz/D");
1265
1266 fPrimaryTree->Branch("fPrimaryAccepted", &fPrimaryAccepted, "fPrimaryAccepted/B");
b2c10007 1267 fPrimaryTree->Branch("fPrimaryMatched", &fPrimaryMatched, "fPrimaryMatched/B");
f61cec2f 1268
1269
1270 fPrimaryTree->Branch("fDepositedCode", &fDepositedCode, "fDepositedCode/I");
1271 fPrimaryTree->Branch("fDepositedCharge", &fDepositedCharge, "fDepositedCharge/I");
b2c10007 1272 fPrimaryTree->Branch("fDepositedE", &fDepositedE, "fDepositedE/D");
f61cec2f 1273 fPrimaryTree->Branch("fDepositedEt", &fDepositedEt, "fDepositedEt/D");
1274
1275 fPrimaryTree->Branch("fDepositedVx", &fDepositedVx, "fDepositedVx/D");
1276 fPrimaryTree->Branch("fDepositedVy", &fDepositedVy, "fDepositedVy/D");
1277 fPrimaryTree->Branch("fDepositedVz", &fDepositedVz, "fDepositedVz/D");
1278
b2c10007 1279 fPrimaryTree->Branch("fSecondary", &fSecondary, "fSecondary/I");
1280
1281
1282 fPrimaryTree->Branch("fReconstructedE", &fReconstructedE, "fReconstructedE/D");
1283 fPrimaryTree->Branch("fReconstructedEt", &fReconstructedEt, "fReconstructedEt/D");
1284
1285 fPrimaryTree->Branch("fClusterMult", &fClusterMult, "fClusterMult/I");
1286
1287
f61cec2f 1288 }
1289
b2c10007 1290 fHistGammasFound = new TH1F("fHistGammasFound", "fHistGammasFound",200, 0, 10);
1291 fHistGammasGenerated = new TH1F("fHistGammasGenerated", "fHistGammasGenerated",200, 0, 10);
ce546038 1292}
1293
0651f6b4 1294void AliAnalysisEtMonteCarlo::FillOutputList(TList *list)
f61cec2f 1295{ //fill the output list
e2ee5727 1296 AliAnalysisEt::FillOutputList(list);
1297
f61cec2f 1298 if(fCuts->GetHistMakeTree())
1299 {
1300 list->Add(fPrimaryTree);
1301 }
1302
e2ee5727 1303 list->Add(fHistRemovedOrNot);
1304
1305 list->Add(fHistEtNonRemovedProtons);
1306 list->Add(fHistEtNonRemovedAntiProtons);
1307 list->Add(fHistEtNonRemovedPiPlus);
1308 list->Add(fHistEtNonRemovedPiMinus);
1309 list->Add(fHistEtNonRemovedKaonPlus);
1310 list->Add(fHistEtNonRemovedKaonMinus);
1311 list->Add(fHistEtNonRemovedK0s);
ef647350 1312 list->Add(fHistEtNonRemovedK0L);
e2ee5727 1313 list->Add(fHistEtNonRemovedLambdas);
1314 list->Add(fHistEtNonRemovedElectrons);
1315 list->Add(fHistEtNonRemovedPositrons);
1316 list->Add(fHistEtNonRemovedMuPlus);
1317 list->Add(fHistEtNonRemovedMuMinus);
1318 list->Add(fHistEtNonRemovedNeutrons);
1319 list->Add(fHistEtNonRemovedAntiNeutrons);
1320 list->Add(fHistEtNonRemovedGammas);
1321 list->Add(fHistEtNonRemovedGammasFromPi0);
1322
1323 list->Add(fHistEtRemovedGammas);
1324 list->Add(fHistEtRemovedNeutrons);
1325 list->Add(fHistEtRemovedAntiNeutrons);
1326
ef647350 1327 list->Add(fHistEtRemovedCharged);
1328 list->Add(fHistEtRemovedNeutrals);
1329
1330 list->Add(fHistEtNonRemovedCharged);
1331 list->Add(fHistEtNonRemovedNeutrals);
e2ee5727 1332
1333 list->Add(fHistMultNonRemovedProtons);
1334 list->Add(fHistMultNonRemovedAntiProtons);
1335 list->Add(fHistMultNonRemovedPiPlus);
1336 list->Add(fHistMultNonRemovedPiMinus);
1337 list->Add(fHistMultNonRemovedKaonPlus);
1338 list->Add(fHistMultNonRemovedKaonMinus);
1339 list->Add(fHistMultNonRemovedK0s);
ef647350 1340 list->Add(fHistMultNonRemovedK0L);
e2ee5727 1341 list->Add(fHistMultNonRemovedLambdas);
1342 list->Add(fHistMultNonRemovedElectrons);
1343 list->Add(fHistMultNonRemovedPositrons);
1344 list->Add(fHistMultNonRemovedMuPlus);
1345 list->Add(fHistMultNonRemovedMuMinus);
1346 list->Add(fHistMultNonRemovedNeutrons);
1347 list->Add(fHistMultNonRemovedAntiNeutrons);
1348 list->Add(fHistMultNonRemovedGammas);
1349
1350 list->Add(fHistMultRemovedGammas);
1351 list->Add(fHistMultRemovedNeutrons);
1352 list->Add(fHistMultRemovedAntiNeutrons);
1353
ef647350 1354 list->Add(fHistMultRemovedCharged);
1355 list->Add(fHistMultRemovedNeutrals);
1356
1357 list->Add(fHistMultNonRemovedCharged);
1358 list->Add(fHistMultNonRemovedNeutrals);
1359
e2ee5727 1360 list->Add(fHistTrackMultvsNonRemovedCharged);
1361 list->Add(fHistTrackMultvsNonRemovedNeutral);
1362 list->Add(fHistTrackMultvsRemovedGamma);
1363
1364 list->Add(fHistClusterMultvsNonRemovedCharged);
1365 list->Add(fHistClusterMultvsNonRemovedNeutral);
1366 list->Add(fHistClusterMultvsRemovedGamma);
1367
1368 //list->Add(fHistDecayVertexNonRemovedCharged);
1369 //list->Add(fHistDecayVertexNonRemovedNeutral);
1370 //list->Add(fHistDecayVertexRemovedCharged);
1371 //list->Add(fHistDecayVertexRemovedNeutral);
1372
1373 list->Add(fHistDxDzNonRemovedCharged);
1374 list->Add(fHistDxDzRemovedCharged);
1375 list->Add(fHistDxDzNonRemovedNeutral);
1376 list->Add(fHistDxDzRemovedNeutral);
1377
2aab9269 1378 list->Add(fHistK0EDepositsVsPtForEOver000MeV);
1379 list->Add(fHistK0EDepositsVsPtForEOver100MeV);
1380 list->Add(fHistK0EDepositsVsPtForEOver150MeV);
1381 list->Add(fHistK0EDepositsVsPtForEOver200MeV);
1382 list->Add(fHistK0EDepositsVsPtForEOver250MeV);
1383 list->Add(fHistK0EDepositsVsPtForEOver300MeV);
1384 list->Add(fHistK0EDepositsVsPtForEOver350MeV);
1385 list->Add(fHistK0EDepositsVsPtForEOver400MeV);
1386 list->Add(fHistK0EDepositsVsPtForEOver450MeV);
1387 list->Add(fHistK0EDepositsVsPtForEOver500MeV);
1388 list->Add(fHistK0EGammaVsPtForEOver000MeV);
1389 list->Add(fHistK0EGammaVsPtForEOver100MeV);
1390 list->Add(fHistK0EGammaVsPtForEOver150MeV);
1391 list->Add(fHistK0EGammaVsPtForEOver200MeV);
1392 list->Add(fHistK0EGammaVsPtForEOver250MeV);
1393 list->Add(fHistK0EGammaVsPtForEOver300MeV);
1394 list->Add(fHistK0EGammaVsPtForEOver350MeV);
1395 list->Add(fHistK0EGammaVsPtForEOver400MeV);
1396 list->Add(fHistK0EGammaVsPtForEOver450MeV);
1397 list->Add(fHistK0EGammaVsPtForEOver500MeV);
1398
e2ee5727 1399 list->Add(fHistPiPlusMult);
1400 list->Add(fHistPiMinusMult);
1401 list->Add(fHistPiZeroMult);
1402 list->Add(fHistPiPlusMultAcc);
1403 list->Add(fHistPiMinusMultAcc);
1404 list->Add(fHistPiZeroMultAcc);
b2c10007 1405
1406 list->Add(fHistGammasFound);
1407 list->Add(fHistGammasGenerated);
e2ee5727 1408
0651f6b4 1409}
1410
1411
13b0d3c1 1412bool AliAnalysisEtMonteCarlo::TrackHitsCalorimeter(TParticle* part, Double_t magField)
1413{
e2ee5727 1414 // printf(" TrackHitsCalorimeter - magField %f\n", magField);
1415 AliESDtrack *esdTrack = new AliESDtrack(part);
1416 // Printf("MC Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
1417
1418 Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
1419
1420 // if(prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
1421
f61cec2f 1422 bool status = prop && fSelector->CutGeometricalAcceptance(*esdTrack);
e2ee5727 1423 delete esdTrack;
1424
1425 return status;
1426}
1427
1428void AliAnalysisEtMonteCarlo::FillHistograms()
f61cec2f 1429{ // let base class fill its histograms, and us fill the local ones
e2ee5727 1430 AliAnalysisEt::FillHistograms();
1431 //std::cout << fEtNonRemovedPiPlus << " " << fCentClass << std::endl;
1432
1433 fHistEtNonRemovedProtons->Fill(fEtNonRemovedProtons, fCentClass);
1434 fHistEtNonRemovedAntiProtons->Fill(fEtNonRemovedAntiProtons, fCentClass);
1435 fHistEtNonRemovedKaonPlus->Fill(fEtNonRemovedKaonPlus, fCentClass);
1436 fHistEtNonRemovedKaonMinus->Fill(fEtNonRemovedKaonMinus, fCentClass);
ef647350 1437 fHistEtNonRemovedK0s->Fill(fEtNonRemovedK0S, fCentClass);
1438 fHistEtNonRemovedK0L->Fill(fEtNonRemovedK0L, fCentClass);
e2ee5727 1439 fHistEtNonRemovedLambdas->Fill(fEtNonRemovedLambdas, fCentClass);
1440 fHistEtNonRemovedPiPlus->Fill(fEtNonRemovedPiPlus, fCentClass);
1441 fHistEtNonRemovedPiMinus->Fill(fEtNonRemovedPiMinus, fCentClass);
1442 fHistEtNonRemovedElectrons->Fill(fEtNonRemovedElectrons, fCentClass);
1443 fHistEtNonRemovedPositrons->Fill(fEtNonRemovedPositrons, fCentClass);
1444 fHistEtNonRemovedMuPlus->Fill(fEtNonRemovedMuPlus, fCentClass);
1445 fHistEtNonRemovedMuMinus->Fill(fEtNonRemovedMuMinus, fCentClass);
1446 fHistEtNonRemovedNeutrons->Fill(fEtNonRemovedNeutrons, fCentClass);
1447 fHistEtNonRemovedAntiNeutrons->Fill(fEtNonRemovedAntiNeutrons, fCentClass);
1448 fHistEtNonRemovedGammas->Fill(fEtNonRemovedGammas, fCentClass);
1449 fHistEtNonRemovedGammasFromPi0->Fill(fEtNonRemovedGammasFromPi0, fCentClass);
1450
ef647350 1451 fHistEtRemovedGammas->Fill(fEtRemovedGammas, fNClusters);
e2ee5727 1452 fHistEtRemovedNeutrons->Fill(fEtRemovedNeutrons, fCentClass);
1453 fHistEtRemovedAntiNeutrons->Fill(fEtRemovedAntiNeutrons, fCentClass);
1454
ef647350 1455// fHistEtRemovedCharged->Fill(fEtRemovedAntiProtons+fEtRemovedElectrons+fEtRemovedKaonMinus+fEtRemovedKaonPlus
1456// +fEtRemovedMuMinus+fEtRemovedMuPlus+fEtRemovedPiMinus+fEtRemovedPiPlus+fEtRemovedPositrons
1457// +fEtRemovedProtons.
1458// fCentClass);
1459// fHistEtRemovedNeutrals->Fill(fEtRemovedNeutrons+fEtRemovedAntiNeutrons, fCentClass);
f61cec2f 1460//
ef647350 1461// fHistEtNonRemovedCharged->Fill(fEtNonRemovedAntiProtons+fEtNonRemovedElectrons+fEtNonRemovedKaonMinus+fEtNonRemovedKaonPlus
1462// +fEtNonRemovedMuMinus+fEtNonRemovedMuPlus+fEtNonRemovedPiMinus+fEtNonRemovedPiPlus+fEtNonRemovedPositrons
1463// +fEtNonRemovedProtons,
1464// fCentClass);
1465// fHistEtRemovedNeutrals->Fill(fEtNonRemovedNeutrons+fEtNonRemovedAntiNeutrons, fCentClass);
1466
1467 fHistEtRemovedCharged->Fill(fEnergyChargedRemoved, fNClusters);
1468 fHistEtRemovedNeutrals->Fill(fEnergyNeutralRemoved, fNClusters);
1469 fHistEtNonRemovedCharged->Fill(fEnergyChargedNotRemoved, fNClusters);
1470 fHistEtNonRemovedNeutrals->Fill(fEnergyNeutralNotRemoved, fNClusters);
f61cec2f 1471
ef647350 1472 fHistMultRemovedCharged->Fill(fChargedRemoved, fNClusters);
1473 fHistMultRemovedNeutrals->Fill(fNeutralRemoved, fNClusters);
1474 fHistMultNonRemovedCharged->Fill(fChargedNotRemoved, fNClusters);
1475 fHistMultNonRemovedNeutrals->Fill(fNeutralNotRemoved, fNClusters);
f61cec2f 1476
1477
e2ee5727 1478 fHistMultNonRemovedProtons->Fill(fMultNonRemovedProtons, fCentClass);
1479 fHistMultNonRemovedAntiProtons->Fill(fMultNonRemovedAntiProtons, fCentClass);
1480 fHistMultNonRemovedKaonPlus->Fill(fMultNonRemovedKaonPlus, fCentClass);
1481 fHistMultNonRemovedKaonMinus->Fill(fMultNonRemovedKaonMinus, fCentClass);
1482 fHistMultNonRemovedK0s->Fill(fMultNonRemovedK0s, fCentClass);
ef647350 1483 fHistMultNonRemovedK0L->Fill(fMultNonRemovedK0L, fCentClass);
e2ee5727 1484 fHistMultNonRemovedLambdas->Fill(fMultNonRemovedLambdas, fCentClass);
1485 fHistMultNonRemovedPiPlus->Fill(fMultNonRemovedPiPlus, fCentClass);
1486 fHistMultNonRemovedPiMinus->Fill(fMultNonRemovedPiMinus, fCentClass);
1487 fHistMultNonRemovedElectrons->Fill(fMultNonRemovedElectrons, fCentClass);
1488 fHistMultNonRemovedPositrons->Fill(fMultNonRemovedPositrons, fCentClass);
1489 fHistMultNonRemovedMuPlus->Fill(fMultNonRemovedMuPlus, fCentClass);
1490 fHistMultNonRemovedMuMinus->Fill(fMultNonRemovedMuMinus, fCentClass);
1491 fHistMultNonRemovedNeutrons->Fill(fMultNonRemovedNeutrons, fCentClass);
1492 fHistMultNonRemovedAntiNeutrons->Fill(fMultNonRemovedAntiNeutrons, fCentClass);
1493 fHistMultNonRemovedGammas->Fill(fMultNonRemovedGammas, fCentClass);
1494
1495 fHistMultRemovedGammas->Fill(fMultRemovedGammas, fCentClass);
1496 fHistMultRemovedNeutrons->Fill(fMultRemovedNeutrons, fCentClass);
1497 fHistMultRemovedAntiNeutrons->Fill(fMultRemovedAntiNeutrons, fCentClass);
1498
1499 fHistTrackMultvsNonRemovedCharged->Fill(fTrackMultInAcc,
1500 fMultNonRemovedAntiProtons+fMultNonRemovedElectrons+fMultNonRemovedKaonMinus+fMultNonRemovedKaonPlus
1501 +fMultNonRemovedMuMinus+fMultNonRemovedMuPlus+fMultNonRemovedPiMinus+fMultNonRemovedPiPlus+fMultNonRemovedPositrons
1502 +fMultNonRemovedProtons);
1503
1504 fHistTrackMultvsNonRemovedNeutral->Fill(fTrackMultInAcc,
b2c10007 1505 fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedK0L+fMultNonRemovedLambdas+fK0sMult);
e2ee5727 1506
1507 fHistTrackMultvsRemovedGamma->Fill(fTrackMultInAcc,
1508 fMultRemovedGammas);
1509
ef647350 1510 fHistClusterMultvsNonRemovedCharged->Fill(fNClusters,
e2ee5727 1511 fMultNonRemovedAntiProtons+fMultNonRemovedElectrons+fMultNonRemovedKaonMinus
1512 +fMultNonRemovedKaonPlus+fMultNonRemovedMuMinus+fMultNonRemovedMuPlus
1513 +fMultNonRemovedPiMinus+fMultNonRemovedPiPlus+fMultNonRemovedPositrons+fMultNonRemovedProtons);
1514
ef647350 1515 fHistClusterMultvsNonRemovedNeutral->Fill(fNClusters,
b2c10007 1516 fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedK0L+fMultNonRemovedLambdas+fK0sMult);
e2ee5727 1517
ef647350 1518 fHistClusterMultvsRemovedGamma->Fill(fNClusters,
e2ee5727 1519 fMultRemovedGammas);
1520
13b0d3c1 1521}
1522
e2ee5727 1523
1524
1525
ef647350 1526Int_t AliAnalysisEtMonteCarlo::PrintFamilyTree(Int_t partIdx, AliStack* stack)
4d376d01 1527{ // print family tree
f61cec2f 1528 TParticle *part = stack->Particle(partIdx);
b2c10007 1529// if(part->GetPdgCode() == fgK0SCode)
f61cec2f 1530 {
1531 std::cout << "This is index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") , is it primary: " << stack->IsPhysicalPrimary(partIdx)<< std::endl;
1532 std::cout << "PID: " << part->GetPdgCode() << "/" << part->GetName() << std::endl;
1533 std::cout << "Energy: " << part->Energy() << std::endl;
1534 std::cout << "Vertex: " << part->Vx() << ", " << part->Vy() << ", " << part->Vz() << std::endl;
1535 }
1536 return PrintMothers(partIdx, stack, 1);
ef647350 1537}
1538
1539Int_t AliAnalysisEtMonteCarlo::PrintMothers(Int_t partIdx, AliStack* stack, Int_t gen)
4d376d01 1540{ // print mothers
f61cec2f 1541 char *tabs = new char[gen+1];
1542 for(Int_t i = 0; i < gen; ++i)
1543 {
1544 //std::cout << i << std::endl;
1545 tabs[i] = '\t';
1546 }
1547 tabs[gen] = '\0';
1548 Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
1549 if(mothIdx < 0)
1550 {
1551 delete [] tabs;
1552 return 0;
1553 }
1554 TParticle *mother = stack->Particle(mothIdx);
b2c10007 1555// if(mother->GetPdgCode() == fgK0SCode)
f61cec2f 1556 {
b2c10007 1557 //std::cout << tabs << "Mother of index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") is: " << mothIdx << ", is it primary: " << stack->IsPhysicalPrimary(mothIdx)<< std::endl;
1558 std::cout << tabs << "Index: " << mothIdx << std::endl;
1559 std::cout << tabs << "Primary: " << stack->IsPhysicalPrimary(mothIdx) << std::endl;
f61cec2f 1560 std::cout << tabs << "PID: " << mother->GetPdgCode() << "/" << mother->GetName() << std::endl;
1561 std::cout << tabs << "Energy: " << mother->Energy() << std::endl;
b2c10007 1562 if(mother->GetFirstMother() >= 0)
1563 {
1564 std::cout << tabs << "Mother(s): " << stack->Particle(mother->GetFirstMother())->GetPdgCode();
1565 if(mother->GetSecondMother() >= 0) std::cout << ", " << stack->Particle(mother->GetSecondMother())->GetPdgCode();
1566 std::cout << std::endl;
1567 }
f61cec2f 1568 std::cout << tabs << "Vertex: " << mother->Vx() << ", " << mother->Vy() << ", " << mother->Vz() << std::endl;
1569 }
1570 if(mother->GetPdgCode() == fgK0SCode)
1571 {
ef647350 1572// std::cout << "K0S!!!!!!!!!!!!!11111!!!!!" << std::endl;
f61cec2f 1573 }
ef647350 1574// std::cout << "Mother of index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") is: " << mothIdx << std::endl;
1575// std::cout << "PID: " << mother->GetPdgCode() << "/" << mother->GetName() << std::endl;
1576// std::cout << "Energy: " << mother->Energy() << std::endl;
1577// std::cout << "Vertex: " << mother->Vx() << ", " << mother->Vy() << ", " << mother->Vz() << std::endl;
f61cec2f 1578
1579 delete [] tabs;
1580 return PrintMothers(mothIdx, stack, gen+1) + 1;
ef647350 1581}
1582
1583Int_t AliAnalysisEtMonteCarlo::GetPrimMother(Int_t partIdx, AliStack *stack)
4d376d01 1584{ // get primary mother
f61cec2f 1585 if(partIdx >= 0)
ef647350 1586 {
b2c10007 1587 //return stack->GetPrimary(partIdx);
1588
f61cec2f 1589 Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
1590 if(mothIdx < 0) return -1;
1591 TParticle *mother = stack->Particle(mothIdx);
1592 if(mother)
1593 {
1594 // if(mother->GetPdgCode() == fgK0SCode)
1595 //{
ef647350 1596// std::cout << "!!!!!!!!!!!!!!!!! K0S !!!!!!!!!!!!!!!!!!" << std::endl;
f61cec2f 1597 //return mothIdx;
1598 // }
1599 //if(mother->GetPdgCode() == fgK0SCode&& stack->IsPhysicalPrimary(mothIdx))
1600 //{
ef647350 1601// std::cout << "!!!!!!!!!!!!!!!!! Primary K0S !!!!!!!!!!!!!!!!!!" << std::endl;
f61cec2f 1602 //return mothIdx;
1603 // }
1604 if(stack->IsPhysicalPrimary(mothIdx)) return mothIdx;
1605 else return GetPrimMother(mothIdx, stack);
1606 }
1607 else
1608 {
1609 return -1;
1610 }
ef647350 1611 }
f61cec2f 1612 return -1;
ef647350 1613}
1614
1615Int_t AliAnalysisEtMonteCarlo::GetK0InFamily(Int_t partIdx, AliStack* stack)
4d376d01 1616{ // get K0 in family
f61cec2f 1617 if(partIdx >= 0)
ef647350 1618 {
2aab9269 1619 if(stack->Particle(partIdx)->GetPdgCode() == fgK0SCode) return partIdx;
f61cec2f 1620 Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
1621 if(mothIdx < 0) return -1;
1622 TParticle *mother = stack->Particle(mothIdx);
1623 if(mother)
ef647350 1624 {
2aab9269 1625 if(mother->GetPdgCode() == fgK0SCode)
ef647350 1626 {
f61cec2f 1627 return mothIdx;
ef647350 1628 }
f61cec2f 1629 return GetK0InFamily(mothIdx, stack);
1630 }
1631 else
1632 {
1633 return -1;
ef647350 1634 }
1635 }
f61cec2f 1636 return -1;
ef647350 1637}
e2ee5727 1638