]> 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
fixing tree saving in the output
[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)
2aab9269 133 ,fHistK0EDepositsVsPtForEOver000MeV(0)
134 ,fHistK0EDepositsVsPtForEOver100MeV(0)
135 ,fHistK0EDepositsVsPtForEOver150MeV(0)
136 ,fHistK0EDepositsVsPtForEOver200MeV(0)
137 ,fHistK0EDepositsVsPtForEOver250MeV(0)
138 ,fHistK0EDepositsVsPtForEOver300MeV(0)
139 ,fHistK0EDepositsVsPtForEOver350MeV(0)
140 ,fHistK0EDepositsVsPtForEOver400MeV(0)
141 ,fHistK0EDepositsVsPtForEOver450MeV(0)
142 ,fHistK0EDepositsVsPtForEOver500MeV(0)
143 ,fHistK0EGammaVsPtForEOver000MeV(0)
144 ,fHistK0EGammaVsPtForEOver100MeV(0)
145 ,fHistK0EGammaVsPtForEOver150MeV(0)
146 ,fHistK0EGammaVsPtForEOver200MeV(0)
147 ,fHistK0EGammaVsPtForEOver250MeV(0)
148 ,fHistK0EGammaVsPtForEOver300MeV(0)
149 ,fHistK0EGammaVsPtForEOver350MeV(0)
150 ,fHistK0EGammaVsPtForEOver400MeV(0)
151 ,fHistK0EGammaVsPtForEOver450MeV(0)
152 ,fHistK0EGammaVsPtForEOver500MeV(0)
f61cec2f 153 ,fHistMultvsRemovedGammaE(0)
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();
687 fPrimaryCharge = primPart->GetPDG()->Charge();
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());
f61cec2f 706 fDepositedCharge = part->GetPDG()->Charge();
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
943 float energyMeasured = 0.0;
944 for (int iCluster = 0; iCluster < nCluster; iCluster++ ){//if this cluster is from any of the decay daughters of our K0S
945 AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
946 const UInt_t myPart = (UInt_t)TMath::Abs(caloCluster->GetLabel());
947 for(int j=0;j<6;j++){
948 if(myPart==gammaDaughterIDs[j]){
949 cout<<"Found a matching cluster!";//<<endl;
950 Double_t clEt = CalculateTransverseEnergy(*caloCluster);
951 cout<<" it has energy "<<clEt<<endl;
952 for(int l=0;l<nEtCuts;l++){//loop over cut values
953 if(clEt>=etCuts[l]){
954 totalClusterEts[l] += clEt;//if cluster et is above the cut off energy add it
955 }
956 }
957
958 }
959 }
960 }
961 //now we have the total energy measured for each case
962 fHistK0EDepositsVsPtForEOver000MeV->Fill(pTk0,totalClusterEts[0]);
963 fHistK0EDepositsVsPtForEOver100MeV->Fill(pTk0,totalClusterEts[1]);
964 fHistK0EDepositsVsPtForEOver150MeV->Fill(pTk0,totalClusterEts[2]);
965 fHistK0EDepositsVsPtForEOver200MeV->Fill(pTk0,totalClusterEts[3]);
966 fHistK0EDepositsVsPtForEOver250MeV->Fill(pTk0,totalClusterEts[4]);
967 fHistK0EDepositsVsPtForEOver300MeV->Fill(pTk0,totalClusterEts[5]);
968 fHistK0EDepositsVsPtForEOver350MeV->Fill(pTk0,totalClusterEts[6]);
969 fHistK0EDepositsVsPtForEOver400MeV->Fill(pTk0,totalClusterEts[7]);
970 fHistK0EDepositsVsPtForEOver450MeV->Fill(pTk0,totalClusterEts[8]);
971 fHistK0EDepositsVsPtForEOver500MeV->Fill(pTk0,totalClusterEts[9]);
972 fHistK0EGammaVsPtForEOver000MeV->Fill(pTk0,totalGammaEts[0]);
973 fHistK0EGammaVsPtForEOver100MeV->Fill(pTk0,totalGammaEts[1]);
974 fHistK0EGammaVsPtForEOver150MeV->Fill(pTk0,totalGammaEts[2]);
975 fHistK0EGammaVsPtForEOver200MeV->Fill(pTk0,totalGammaEts[3]);
976 fHistK0EGammaVsPtForEOver250MeV->Fill(pTk0,totalGammaEts[4]);
977 fHistK0EGammaVsPtForEOver300MeV->Fill(pTk0,totalGammaEts[5]);
978 fHistK0EGammaVsPtForEOver350MeV->Fill(pTk0,totalGammaEts[6]);
979 fHistK0EGammaVsPtForEOver400MeV->Fill(pTk0,totalGammaEts[7]);
980 fHistK0EGammaVsPtForEOver450MeV->Fill(pTk0,totalGammaEts[8]);
981 fHistK0EGammaVsPtForEOver500MeV->Fill(pTk0,totalGammaEts[9]);
982 }
983 }
984 }
985 }
986 }
987
e2ee5727 988 FillHistograms();
989 return 0;
990}
2fbf38ac 991
992void AliAnalysisEtMonteCarlo::Init()
f61cec2f 993{ // init
e2ee5727 994 AliAnalysisEt::Init();
2fbf38ac 995}
13b0d3c1 996
ce546038 997void AliAnalysisEtMonteCarlo::ResetEventValues()
f61cec2f 998{ // reset event values
e2ee5727 999 AliAnalysisEt::ResetEventValues();
1000
f61cec2f 1001 fTotEtSecondary = 0;
1002 fTotEtSecondaryFromEmEtPrimary = 0;
1003 fTotEtWithSecondaryRemoved = 0;
1004
e2ee5727 1005 // collision geometry defaults for p+p:
1006 fImpactParameter = 0;
1007 fNcoll = 1;
1008 fNpart = 2;
1009
1010 fEtNonRemovedProtons = 0;
1011 fEtNonRemovedAntiProtons = 0;
1012 fEtNonRemovedPiPlus = 0;
1013 fEtNonRemovedPiMinus = 0;
1014 fEtNonRemovedKaonPlus = 0;
1015 fEtNonRemovedKaonMinus = 0;
ef647350 1016 fEtNonRemovedK0S = 0;
1017 fEtNonRemovedK0L = 0;
e2ee5727 1018 fEtNonRemovedLambdas = 0;
1019 fEtNonRemovedElectrons = 0;
1020 fEtNonRemovedPositrons = 0;
1021 fEtNonRemovedMuPlus = 0;
1022 fEtNonRemovedMuMinus = 0;
1023 fEtNonRemovedNeutrons = 0;
1024 fEtNonRemovedAntiNeutrons = 0;
1025 fEtNonRemovedGammas = 0;
1026 fEtNonRemovedGammasFromPi0 = 0;
1027
ef647350 1028 fEtRemovedProtons = 0;
1029 fEtRemovedAntiProtons = 0;
1030 fEtRemovedPiPlus = 0;
1031 fEtRemovedPiMinus = 0;
1032 fEtRemovedKaonPlus = 0;
1033 fEtRemovedKaonMinus = 0;
1034 fEtRemovedK0s = 0;
1035 fEtRemovedK0L = 0;
1036 fEtRemovedLambdas = 0;
1037 fEtRemovedElectrons = 0;
1038 fEtRemovedPositrons = 0;
1039 fEtRemovedMuPlus = 0;
1040 fEtRemovedMuMinus = 0;
1041 fEtRemovedNeutrons = 0;
1042
1043 fEtRemovedGammasFromPi0 = 0;
e2ee5727 1044 fEtRemovedGammas = 0;
1045 fEtRemovedNeutrons = 0;
1046 fEtRemovedAntiNeutrons = 0;
1047
1048 fMultNonRemovedProtons = 0;
1049 fMultNonRemovedAntiProtons = 0;
1050 fMultNonRemovedPiPlus = 0;
1051 fMultNonRemovedPiMinus = 0;
1052 fMultNonRemovedKaonPlus = 0;
1053 fMultNonRemovedKaonMinus = 0;
1054 fMultNonRemovedK0s = 0;
ef647350 1055 fMultNonRemovedK0L = 0;
e2ee5727 1056 fMultNonRemovedLambdas = 0;
1057 fMultNonRemovedElectrons = 0;
1058 fMultNonRemovedPositrons = 0;
1059 fMultNonRemovedMuPlus = 0;
1060 fMultNonRemovedMuMinus = 0;
1061 fMultNonRemovedNeutrons = 0;
1062 fMultNonRemovedAntiNeutrons = 0;
1063 fMultNonRemovedGammas = 0;
1064
ef647350 1065 fMultRemovedProtons = 0;
1066 fMultRemovedAntiProtons = 0;
1067 fMultRemovedPiPlus = 0;
1068 fMultRemovedPiMinus = 0;
1069 fMultRemovedKaonPlus = 0;
1070 fMultRemovedKaonMinus = 0;
1071 fMultRemovedK0s = 0;
1072 fMultRemovedK0L = 0;
1073 fMultRemovedLambdas = 0;
1074 fMultRemovedElectrons = 0;
1075 fMultRemovedPositrons = 0;
1076 fMultRemovedMuPlus = 0;
1077 fMultRemovedMuMinus = 0;
1078
e2ee5727 1079 fMultRemovedGammas = 0;
1080 fMultRemovedNeutrons = 0;
1081 fMultRemovedAntiNeutrons = 0;
1082
ef647350 1083 fEnergyChargedNotRemoved = 0;
1084 fEnergyChargedRemoved = 0;
1085 fEnergyNeutralNotRemoved = 0;
1086 fEnergyNeutralRemoved = 0;
f61cec2f 1087
ef647350 1088 fChargedNotRemoved = 0;
1089 fChargedRemoved = 0;
1090 fNeutralNotRemoved = 0;
1091 fNeutralRemoved = 0;
b2c10007 1092 fGammaRemoved = 0;
1093 fSecondaryNotRemoved = 0;
f61cec2f 1094
e2ee5727 1095 fTrackMultInAcc = 0;
1096
f61cec2f 1097 fTotNeutralEtAfterMinEnergyCut = 0;
b2c10007 1098
1099 fSecondaryNotRemoved = 0;
1100
1101 fTotPx = 0;
1102 fTotPy = 0;
1103 fTotPz = 0;
1104
1105
ce546038 1106}
1107
1108void AliAnalysisEtMonteCarlo::CreateHistograms()
f61cec2f 1109{ // histogram related additions
e2ee5727 1110 AliAnalysisEt::CreateHistograms();
f61cec2f 1111 if (fEventSummaryTree) {
1112 fEventSummaryTree->Branch("fImpactParameter",&fImpactParameter,"fImpactParameter/D");
1113 fEventSummaryTree->Branch("fNcoll",&fNcoll,"fNcoll/I");
1114 fEventSummaryTree->Branch("fNpart",&fNpart,"fNpart/I");
1115 fEventSummaryTree->Branch("fTotEtWithSecondaryRemoved", &fTotEtWithSecondaryRemoved, "fTotEtWithSecondaryRemoved/D");
1116 fEventSummaryTree->Branch("fTotEtSecondaryFromEmEtPrimary", &fTotEtSecondaryFromEmEtPrimary, "fTotEtSecondaryFromEmEtPrimary/D");
1117 fEventSummaryTree->Branch("fTotEtSecondary", &fTotEtSecondary, "fTotEtSecondary/D");
1118 fEventSummaryTree->Branch("fTotNeutralEtAfterMinEnergyCut", &fTotNeutralEtAfterMinEnergyCut, "fTotNeutralEtAfterMinEnergyCut/D");
b2c10007 1119 fEventSummaryTree->Branch("fSecondaryNotRemoved", &fSecondaryNotRemoved, "fSecondaryNotRemoved/I");
1120 fEventSummaryTree->Branch("fChargedNotRemoved", &fChargedNotRemoved, "fChargedNotRemoved/I");
1121 fEventSummaryTree->Branch("fNeutralNotRemoved", &fNeutralNotRemoved, "fNeutralNotRemoved/I");
1122 fEventSummaryTree->Branch("fChargedRemoved", &fChargedRemoved, "fChargedRemoved/I");
1123 fEventSummaryTree->Branch("fNeutralRemoved", &fNeutralRemoved, "fNeutralRemoved/I");
1124 fEventSummaryTree->Branch("fGammaRemoved", &fGammaRemoved, "fGammaRemoved/I");
1125 fEventSummaryTree->Branch("fTotPx", &fTotPx, "fTotPx/D");
1126 fEventSummaryTree->Branch("fTotPy", &fTotPy, "fTotPy/D");
1127 fEventSummaryTree->Branch("fTotPz", &fTotPz, "fTotPz/D");
1128// fEventSummaryTree->Branch("f
e2ee5727 1129 }
f61cec2f 1130
e2ee5727 1131 //fHistDecayVertexNonRemovedCharged = new TH3F("fHistDecayVertexNonRemovedCharged","fHistDecayVertexNonRemovedCharged", 500, -470, 30, 500, -300, 300, 40, -20, 20);
1132 //fHistDecayVertexRemovedCharged = new TH3F("fHistDecayVertexRemovedCharged","fHistDecayVertexRemovedCharged", 500, -470, 30, 500, -300, 300, 40, -20, 20);
1133 //fHistDecayVertexNonRemovedNeutral = new TH3F("fHistDecayVertexNonRemovedNeutral","fHistDecayVertexNonRemovedNeutral", 500, -470, 30, 500, -300, 300, 40, -20, 20);
1134 //fHistDecayVertexRemovedNeutral = new TH3F("fHistDecayVertexRemovedNeutral","fHistDecayVertexRemovedNeutral", 500, -470, 30, 500, -300, 300, 40, -20, 20);
1135
1136 fHistRemovedOrNot = new TH2F("fHistRemovedOrNot", "fHistRemovedOrNot", 4, -0.5, 3.5, 10, -0.5, 9.5);
1137
1138 fHistEtNonRemovedProtons = new TH2F("fHistEtNonRemovedProtons", "fHistEtNonRemovedProtons", 1500, 0, 30, 10, -0.5, 9.5);
1139 fHistEtNonRemovedAntiProtons = new TH2F("fHistEtNonRemovedAntiProtons", "fHistEtNonRemovedAntiProtons", 1500, 0, 30, 10, -0.5, 9.5);
1140 fHistEtNonRemovedPiPlus = new TH2F("fHistEtNonRemovedPiPlus", "fHistEtNonRemovedPiPlus", 1500, 0, 30, 10, -0.5, 9.5);
1141 fHistEtNonRemovedPiMinus = new TH2F("fHistEtNonRemovedPiMinus", "fHistEtNonRemovedPiMinus", 1500, 0, 30, 10, -0.5, 9.5);
1142 fHistEtNonRemovedKaonPlus = new TH2F("fHistEtNonRemovedKaonPlus", "fHistEtNonRemovedKaonPlus", 1500, 0, 30, 10, -0.5, 9.5);
1143 fHistEtNonRemovedKaonMinus = new TH2F("fHistEtNonRemovedKaonMinus", "fHistEtNonRemovedKaonMinus", 1500, 0, 30, 10, -0.5, 9.5);
1144 fHistEtNonRemovedK0s = new TH2F("fHistEtNonRemovedK0s", "fHistEtNonRemovedK0s", 1500, 0, 30, 10, -0.5, 9.5);
ef647350 1145 fHistEtNonRemovedK0L = new TH2F("fHistEtNonRemovedK0L", "fHistEtNonRemovedK0L", 1500, 0, 30, 10, -0.5, 9.5);
e2ee5727 1146 fHistEtNonRemovedLambdas = new TH2F("fHistEtNonRemovedLambdas", "fHistEtNonRemovedLambdas", 1500, 0, 30, 10, -0.5, 9.5);
1147 fHistEtNonRemovedElectrons = new TH2F("fHistEtNonRemovedElectrons", "fHistEtNonRemovedElectrons", 1500, 0, 30, 10, -0.5, 9.5);
1148 fHistEtNonRemovedPositrons = new TH2F("fHistEtNonRemovedPositrons", "fHistEtNonRemovedPositrons", 1500, 0, 30, 10, -0.5, 9.5);
1149 fHistEtNonRemovedMuPlus = new TH2F("fHistEtNonRemovedMuPlus", "fHistEtNonRemovedMuPlus", 1500, 0, 30, 10, -0.5, 9.5);
1150 fHistEtNonRemovedMuMinus = new TH2F("fHistEtNonRemovedMuMinus", "fHistEtNonRemovedMuMinus", 1500, 0, 30, 10, -0.5, 9.5);
1151 fHistEtNonRemovedNeutrons = new TH2F("fHistEtNonRemovedNeutrons", "fHistEtNonRemovedNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
1152 fHistEtNonRemovedAntiNeutrons = new TH2F("fHistEtNonRemovedAntiNeutrons", "fHistEtNonRemovedAntiNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
1153
1154 fHistEtNonRemovedGammas = new TH2F("fHistEtNonRemovedGammas", "fHistEtNonRemovedGammas", 1500, 0, 30, 10, -0.5, 9.5);
1155 fHistEtNonRemovedGammasFromPi0 = new TH2F("fHistEtNonRemovedGammasFromPi0", "fHistEtNonRemovedGammasFromPi0", 1500, 0, 30, 10, -0.5, 9.5);
1156
1157 fHistEtRemovedGammas = new TH2F("fHistEtRemovedGammas", "fHistEtRemovedGammas", 1500, 0, 30, 10, -0.5, 9.5);
1158 fHistEtRemovedNeutrons = new TH2F("fHistEtRemovedNeutrons", "fHistEtRemovedNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
1159 fHistEtRemovedAntiNeutrons = new TH2F("fHistEtRemovedAntiNeutrons", "fHistEtRemovedAntiNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
f61cec2f 1160
ef647350 1161 fHistEtRemovedCharged = new TH2F("fHistEtRemovedCharged", "fHistEtRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
1162 fHistEtRemovedNeutrals = new TH2F("fHistEtRemovedNeutrals", "fHistEtRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
e2ee5727 1163
ef647350 1164 fHistEtNonRemovedCharged = new TH2F("fHistEtNonRemovedCharged", "fHistEtNonRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
1165 fHistEtNonRemovedNeutrals = new TH2F("fHistEtNonRemovedNeutrals", "fHistEtNonRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
f61cec2f 1166
e2ee5727 1167 fHistMultNonRemovedProtons = new TH2F("fHistMultNonRemovedProtons", "fHistMultNonRemovedProtons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1168 fHistMultNonRemovedAntiProtons = new TH2F("fHistMultNonRemovedAntiProtons", "fHistMultNonRemovedAntiProtons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1169 fHistMultNonRemovedPiPlus = new TH2F("fHistMultNonRemovedPiPlus", "fHistMultNonRemovedPiPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1170 fHistMultNonRemovedPiMinus = new TH2F("fHistMultNonRemovedPiMinus", "fHistMultNonRemovedPiMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1171 fHistMultNonRemovedKaonPlus = new TH2F("fHistMultNonRemovedKaonPlus", "fHistMultNonRemovedKaonPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1172 fHistMultNonRemovedKaonMinus = new TH2F("fHistMultNonRemovedKaonMinus", "fHistMultNonRemovedKaonMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1173 fHistMultNonRemovedK0s = new TH2F("fHistMultNonRemovedK0s", "fHistMultNonRemovedK0s", 100, -0.5, 99.5, 10, -0.5, 9.5);
ef647350 1174 fHistMultNonRemovedK0L = new TH2F("fHistMultNonRemovedK0L", "fHistMultNonRemovedK0L", 100, -0.5, 99.5, 10, -0.5, 9.5);
e2ee5727 1175 fHistMultNonRemovedLambdas = new TH2F("fHistMultNonRemovedLambdas", "fHistMultNonRemovedLambdas", 100, -0.5, 99.5, 10, -0.5, 9.5);
1176 fHistMultNonRemovedElectrons = new TH2F("fHistMultNonRemovedElectrons", "fHistMultNonRemovedElectrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1177 fHistMultNonRemovedPositrons = new TH2F("fHistMultNonRemovedPositrons", "fHistMultNonRemovedPositrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1178 fHistMultNonRemovedMuPlus = new TH2F("fHistMultNonRemovedMuPlus", "fHistMultNonRemovedMuPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1179 fHistMultNonRemovedMuMinus = new TH2F("fHistMultNonRemovedMuMinus", "fHistMultNonRemovedMuMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1180 fHistMultNonRemovedNeutrons = new TH2F("fHistMultNonRemovedNeutrons", "fHistMultNonRemovedNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1181 fHistMultNonRemovedAntiNeutrons = new TH2F("fHistMultNonRemovedAntiNeutrons", "fHistMultNonRemovedAntiNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1182
ef647350 1183 fHistMultNonRemovedGammas = new TH2F("fHistMultNonRemovedGammas", "fHistMultNonRemovedGammas", 100, -0.5, 99.5, 100, -0.5, 99.5);
e2ee5727 1184
ef647350 1185 fHistMultRemovedGammas = new TH2F("fHistMultRemovedGammas", "fHistMultRemovedGammas", 100, -0.5, 99.5, 100, -0.5, 99.5);
e2ee5727 1186 fHistMultRemovedNeutrons = new TH2F("fHistMultRemovedNeutrons", "fHistMultRemovedNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1187 fHistMultRemovedAntiNeutrons = new TH2F("fHistMultRemovedAntiNeutrons", "fHistMultRemovedAntiNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
f61cec2f 1188 /*
1189 fHistMultRemovedCharged = new TH2F("fHistMultRemovedCharged", "fHistMultRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
1190 fHistMultRemovedNeutrals = new TH2F("fHistMultRemovedNeutrals", "fHistMultRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
1191
1192 fHistMultNonRemovedCharged = new TH2F("fHistMultNonRemovedCharged", "fHistMultNonRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
1193 fHistMultNonRemovedNeutrals = new TH2F("fHistMultNonRemovedNeutrals", "fHistMultNonRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);*/
ef647350 1194
ef647350 1195
ef647350 1196 fHistMultRemovedCharged = new TH2F("fHistMultRemovedCharged", "fHistMultRemovedCharged", 100, -0.5, 99.5, 100, -0.5, 99.5);
1197 fHistMultRemovedNeutrals = new TH2F("fHistMultRemovedNeutrals", "fHistMultRemovedNeutrals", 100, -0.5, 99.5, 100, -0.5, 99.5);
1198
1199 fHistMultNonRemovedCharged = new TH2F("fHistMultNonRemovedCharged", "fHistMultNonRemovedCharged", 100, -0.5, 99.5, 100, -0.5, 99.5);
1200 fHistMultNonRemovedNeutrals = new TH2F("fHistMultNonRemovedNeutrals", "fHistMultNonRemovedNeutrals", 100, -0.5, 99.5, 100, -0.5, 99.5);
e2ee5727 1201
1202 fHistTrackMultvsNonRemovedCharged = new TH2F("fHistTrackMultvsNonRemovedCharged", "fHistTrackMultvsNonRemovedCharged", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1203 fHistTrackMultvsNonRemovedNeutral = new TH2F("fHistTrackMultvsNonRemovedNeutral", "fHistTrackMultvsNonRemovedNeutral", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1204 fHistTrackMultvsRemovedGamma = new TH2F("fHistTrackMultvsRemovedGamma", "fHistTrackMultvsRemovedGamma", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1205
1206 fHistClusterMultvsNonRemovedCharged = new TH2F("fHistClusterMultvsNonRemovedCharged", "fHistClusterMultvsNonRemovedCharged", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1207 fHistClusterMultvsNonRemovedNeutral = new TH2F("fHistClusterMultvsNonRemovedNeutral", "fHistClusterMultvsNonRemovedNeutral", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1208 fHistClusterMultvsRemovedGamma = new TH2F("fHistClusterMultvsRemovedGamma", "fHistClusterMultvsRemovedGamma", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1209
1210 fHistDxDzNonRemovedCharged = new TH2F("fHistDxDzNonRemovedCharged", "fHistDxDzNonRemovedCharged", 800, -200, 200, 800, -200, 200);
1211 fHistDxDzRemovedCharged = new TH2F("fHistDxDzRemovedCharged", "fHistDxDzRemovedCharged", 800, -200, 200, 800, -200, 200);
1212 fHistDxDzNonRemovedNeutral = new TH2F("fHistDxDzNonRemovedNeutral", "fHistDxDzNonRemovedNeutral", 800, -200, 200, 800, -200, 200);
1213 fHistDxDzRemovedNeutral = new TH2F("fHistDxDzRemovedNeutral", "fHistDxDzRemovedNeutral", 800, -200, 200, 800, -200, 200);
1214
2aab9269 1215 fHistK0EDepositsVsPtForEOver000MeV = new TH2F("fHistK0EDepositsVsPtForEOver000MeV","K^{0}_{S} deposits over 0 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1216 fHistK0EDepositsVsPtForEOver100MeV = new TH2F("fHistK0EDepositsVsPtForEOver100MeV","K^{0}_{S} deposits over 100 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1217 fHistK0EDepositsVsPtForEOver150MeV = new TH2F("fHistK0EDepositsVsPtForEOver150MeV","K^{0}_{S} deposits over 150 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1218 fHistK0EDepositsVsPtForEOver200MeV = new TH2F("fHistK0EDepositsVsPtForEOver200MeV","K^{0}_{S} deposits over 200 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1219 fHistK0EDepositsVsPtForEOver250MeV = new TH2F("fHistK0EDepositsVsPtForEOver250MeV","K^{0}_{S} deposits over 250 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1220 fHistK0EDepositsVsPtForEOver300MeV = new TH2F("fHistK0EDepositsVsPtForEOver300MeV","K^{0}_{S} deposits over 300 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1221 fHistK0EDepositsVsPtForEOver350MeV = new TH2F("fHistK0EDepositsVsPtForEOver350MeV","K^{0}_{S} deposits over 350 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1222 fHistK0EDepositsVsPtForEOver400MeV = new TH2F("fHistK0EDepositsVsPtForEOver400MeV","K^{0}_{S} deposits over 400 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1223 fHistK0EDepositsVsPtForEOver450MeV = new TH2F("fHistK0EDepositsVsPtForEOver450MeV","K^{0}_{S} deposits over 450 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1224 fHistK0EDepositsVsPtForEOver500MeV = new TH2F("fHistK0EDepositsVsPtForEOver500MeV","K^{0}_{S} deposits over 500 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1225 fHistK0EGammaVsPtForEOver000MeV = new TH2F("fHistK0EGammaVsPtForEOver000MeV","K^{0}_{S} gamma over 0 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1226 fHistK0EGammaVsPtForEOver100MeV = new TH2F("fHistK0EGammaVsPtForEOver100MeV","K^{0}_{S} gamma over 100 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1227 fHistK0EGammaVsPtForEOver150MeV = new TH2F("fHistK0EGammaVsPtForEOver150MeV","K^{0}_{S} gamma over 150 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1228 fHistK0EGammaVsPtForEOver200MeV = new TH2F("fHistK0EGammaVsPtForEOver200MeV","K^{0}_{S} gamma over 200 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1229 fHistK0EGammaVsPtForEOver250MeV = new TH2F("fHistK0EGammaVsPtForEOver250MeV","K^{0}_{S} gamma over 250 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1230 fHistK0EGammaVsPtForEOver300MeV = new TH2F("fHistK0EGammaVsPtForEOver300MeV","K^{0}_{S} gamma over 300 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1231 fHistK0EGammaVsPtForEOver350MeV = new TH2F("fHistK0EGammaVsPtForEOver350MeV","K^{0}_{S} gamma over 350 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1232 fHistK0EGammaVsPtForEOver400MeV = new TH2F("fHistK0EGammaVsPtForEOver400MeV","K^{0}_{S} gamma over 400 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1233 fHistK0EGammaVsPtForEOver450MeV = new TH2F("fHistK0EGammaVsPtForEOver450MeV","K^{0}_{S} gamma over 450 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1234 fHistK0EGammaVsPtForEOver500MeV = new TH2F("fHistK0EGammaVsPtForEOver500MeV","K^{0}_{S} gamma over 500 MeV",fgNumOfPtBins,fgPtAxis,fgNumOfPtBins,fgPtAxis);
1235
e2ee5727 1236 fHistPiPlusMult = new TH1F("fHistPiPlusMult", "fHistPiPlusMult", 2000, -0.5, 1999.5);
1237 fHistPiMinusMult = new TH1F("fHistPiMinusMult", "fHistPiMinusMult", 2000, -0.5, 1999.5);
1238 fHistPiZeroMult = new TH1F("fHistPiZeroMult", "fHistPiZeroMult", 2000, -0.5, 1999.5);
1239
1240 fHistPiPlusMultAcc = new TH1F("fHistPiPlusMultAcc", "fHistPiPlusMultAcc", 2000, -0.5, 1999.5);
1241 fHistPiMinusMultAcc = new TH1F("fHistPiMinusMultAcc", "fHistPiMinusMultAcc", 2000, -0.5, 1999.5);
1242 fHistPiZeroMultAcc = new TH1F("fHistPiZeroMultAcc", "fHistPiZeroMultAcc", 2000, -0.5, 1999.5);
1243
f61cec2f 1244 if(fCuts->GetHistMakeTree())
1245 {
1246 TString treename = "fPrimaryTree" + fHistogramNameSuffix;
1247 fPrimaryTree = new TTree(treename, treename);
1248
1249 fPrimaryTree->Branch("fTotEt",&fTotEt,"fTotEt/D");
1250 fPrimaryTree->Branch("fNeutralMultiplicity", &fNeutralMultiplicity, "fNeutralMultiplicity/I");
1251 fPrimaryTree->Branch("fCentClass",&fCentClass,"fCentClass/I");
1252
1253 fPrimaryTree->Branch("fPrimaryCode", &fPrimaryCode, "fPrimaryCode/I");
1254 fPrimaryTree->Branch("fPrimaryCharge", &fPrimaryCharge, "fPrimaryCharge/I");
1255
1256 fPrimaryTree->Branch("fPrimaryE", &fPrimaryE, "fPrimaryE/D");
1257 fPrimaryTree->Branch("fPrimaryEt", &fPrimaryEt, "fPrimaryEt/D");
1258
1259 fPrimaryTree->Branch("fPrimaryPx", &fPrimaryPx, "fPrimaryPx/D");
1260 fPrimaryTree->Branch("fPrimaryPy", &fPrimaryPy, "fPrimaryPy/D");
1261 fPrimaryTree->Branch("fPrimaryPz", &fPrimaryPz, "fPrimaryPz/D");
1262
1263 fPrimaryTree->Branch("fPrimaryVx", &fPrimaryVx, "fPrimaryVx/D");
1264 fPrimaryTree->Branch("fPrimaryVy", &fPrimaryVy, "fPrimaryVy/D");
1265 fPrimaryTree->Branch("fPrimaryVz", &fPrimaryVz, "fPrimaryVz/D");
1266
1267 fPrimaryTree->Branch("fPrimaryAccepted", &fPrimaryAccepted, "fPrimaryAccepted/B");
b2c10007 1268 fPrimaryTree->Branch("fPrimaryMatched", &fPrimaryMatched, "fPrimaryMatched/B");
f61cec2f 1269
1270
1271 fPrimaryTree->Branch("fDepositedCode", &fDepositedCode, "fDepositedCode/I");
1272 fPrimaryTree->Branch("fDepositedCharge", &fDepositedCharge, "fDepositedCharge/I");
b2c10007 1273 fPrimaryTree->Branch("fDepositedE", &fDepositedE, "fDepositedE/D");
f61cec2f 1274 fPrimaryTree->Branch("fDepositedEt", &fDepositedEt, "fDepositedEt/D");
1275
1276 fPrimaryTree->Branch("fDepositedVx", &fDepositedVx, "fDepositedVx/D");
1277 fPrimaryTree->Branch("fDepositedVy", &fDepositedVy, "fDepositedVy/D");
1278 fPrimaryTree->Branch("fDepositedVz", &fDepositedVz, "fDepositedVz/D");
1279
b2c10007 1280 fPrimaryTree->Branch("fSecondary", &fSecondary, "fSecondary/I");
1281
1282
1283 fPrimaryTree->Branch("fReconstructedE", &fReconstructedE, "fReconstructedE/D");
1284 fPrimaryTree->Branch("fReconstructedEt", &fReconstructedEt, "fReconstructedEt/D");
1285
1286 fPrimaryTree->Branch("fClusterMult", &fClusterMult, "fClusterMult/I");
1287
1288
f61cec2f 1289 }
1290
b2c10007 1291 fHistGammasFound = new TH1F("fHistGammasFound", "fHistGammasFound",200, 0, 10);
1292 fHistGammasGenerated = new TH1F("fHistGammasGenerated", "fHistGammasGenerated",200, 0, 10);
ce546038 1293}
1294
0651f6b4 1295void AliAnalysisEtMonteCarlo::FillOutputList(TList *list)
f61cec2f 1296{ //fill the output list
e2ee5727 1297 AliAnalysisEt::FillOutputList(list);
1298
f61cec2f 1299 if(fCuts->GetHistMakeTree())
1300 {
1301 list->Add(fPrimaryTree);
1302 }
1303
e2ee5727 1304 list->Add(fHistRemovedOrNot);
1305
1306 list->Add(fHistEtNonRemovedProtons);
1307 list->Add(fHistEtNonRemovedAntiProtons);
1308 list->Add(fHistEtNonRemovedPiPlus);
1309 list->Add(fHistEtNonRemovedPiMinus);
1310 list->Add(fHistEtNonRemovedKaonPlus);
1311 list->Add(fHistEtNonRemovedKaonMinus);
1312 list->Add(fHistEtNonRemovedK0s);
ef647350 1313 list->Add(fHistEtNonRemovedK0L);
e2ee5727 1314 list->Add(fHistEtNonRemovedLambdas);
1315 list->Add(fHistEtNonRemovedElectrons);
1316 list->Add(fHistEtNonRemovedPositrons);
1317 list->Add(fHistEtNonRemovedMuPlus);
1318 list->Add(fHistEtNonRemovedMuMinus);
1319 list->Add(fHistEtNonRemovedNeutrons);
1320 list->Add(fHistEtNonRemovedAntiNeutrons);
1321 list->Add(fHistEtNonRemovedGammas);
1322 list->Add(fHistEtNonRemovedGammasFromPi0);
1323
1324 list->Add(fHistEtRemovedGammas);
1325 list->Add(fHistEtRemovedNeutrons);
1326 list->Add(fHistEtRemovedAntiNeutrons);
1327
ef647350 1328 list->Add(fHistEtRemovedCharged);
1329 list->Add(fHistEtRemovedNeutrals);
1330
1331 list->Add(fHistEtNonRemovedCharged);
1332 list->Add(fHistEtNonRemovedNeutrals);
e2ee5727 1333
1334 list->Add(fHistMultNonRemovedProtons);
1335 list->Add(fHistMultNonRemovedAntiProtons);
1336 list->Add(fHistMultNonRemovedPiPlus);
1337 list->Add(fHistMultNonRemovedPiMinus);
1338 list->Add(fHistMultNonRemovedKaonPlus);
1339 list->Add(fHistMultNonRemovedKaonMinus);
1340 list->Add(fHistMultNonRemovedK0s);
ef647350 1341 list->Add(fHistMultNonRemovedK0L);
e2ee5727 1342 list->Add(fHistMultNonRemovedLambdas);
1343 list->Add(fHistMultNonRemovedElectrons);
1344 list->Add(fHistMultNonRemovedPositrons);
1345 list->Add(fHistMultNonRemovedMuPlus);
1346 list->Add(fHistMultNonRemovedMuMinus);
1347 list->Add(fHistMultNonRemovedNeutrons);
1348 list->Add(fHistMultNonRemovedAntiNeutrons);
1349 list->Add(fHistMultNonRemovedGammas);
1350
1351 list->Add(fHistMultRemovedGammas);
1352 list->Add(fHistMultRemovedNeutrons);
1353 list->Add(fHistMultRemovedAntiNeutrons);
1354
ef647350 1355 list->Add(fHistMultRemovedCharged);
1356 list->Add(fHistMultRemovedNeutrals);
1357
1358 list->Add(fHistMultNonRemovedCharged);
1359 list->Add(fHistMultNonRemovedNeutrals);
1360
e2ee5727 1361 list->Add(fHistTrackMultvsNonRemovedCharged);
1362 list->Add(fHistTrackMultvsNonRemovedNeutral);
1363 list->Add(fHistTrackMultvsRemovedGamma);
1364
1365 list->Add(fHistClusterMultvsNonRemovedCharged);
1366 list->Add(fHistClusterMultvsNonRemovedNeutral);
1367 list->Add(fHistClusterMultvsRemovedGamma);
1368
1369 //list->Add(fHistDecayVertexNonRemovedCharged);
1370 //list->Add(fHistDecayVertexNonRemovedNeutral);
1371 //list->Add(fHistDecayVertexRemovedCharged);
1372 //list->Add(fHistDecayVertexRemovedNeutral);
1373
1374 list->Add(fHistDxDzNonRemovedCharged);
1375 list->Add(fHistDxDzRemovedCharged);
1376 list->Add(fHistDxDzNonRemovedNeutral);
1377 list->Add(fHistDxDzRemovedNeutral);
1378
2aab9269 1379 list->Add(fHistK0EDepositsVsPtForEOver000MeV);
1380 list->Add(fHistK0EDepositsVsPtForEOver100MeV);
1381 list->Add(fHistK0EDepositsVsPtForEOver150MeV);
1382 list->Add(fHistK0EDepositsVsPtForEOver200MeV);
1383 list->Add(fHistK0EDepositsVsPtForEOver250MeV);
1384 list->Add(fHistK0EDepositsVsPtForEOver300MeV);
1385 list->Add(fHistK0EDepositsVsPtForEOver350MeV);
1386 list->Add(fHistK0EDepositsVsPtForEOver400MeV);
1387 list->Add(fHistK0EDepositsVsPtForEOver450MeV);
1388 list->Add(fHistK0EDepositsVsPtForEOver500MeV);
1389 list->Add(fHistK0EGammaVsPtForEOver000MeV);
1390 list->Add(fHistK0EGammaVsPtForEOver100MeV);
1391 list->Add(fHistK0EGammaVsPtForEOver150MeV);
1392 list->Add(fHistK0EGammaVsPtForEOver200MeV);
1393 list->Add(fHistK0EGammaVsPtForEOver250MeV);
1394 list->Add(fHistK0EGammaVsPtForEOver300MeV);
1395 list->Add(fHistK0EGammaVsPtForEOver350MeV);
1396 list->Add(fHistK0EGammaVsPtForEOver400MeV);
1397 list->Add(fHistK0EGammaVsPtForEOver450MeV);
1398 list->Add(fHistK0EGammaVsPtForEOver500MeV);
1399
e2ee5727 1400 list->Add(fHistPiPlusMult);
1401 list->Add(fHistPiMinusMult);
1402 list->Add(fHistPiZeroMult);
1403 list->Add(fHistPiPlusMultAcc);
1404 list->Add(fHistPiMinusMultAcc);
1405 list->Add(fHistPiZeroMultAcc);
b2c10007 1406
1407 list->Add(fHistGammasFound);
1408 list->Add(fHistGammasGenerated);
e2ee5727 1409
0651f6b4 1410}
1411
1412
13b0d3c1 1413bool AliAnalysisEtMonteCarlo::TrackHitsCalorimeter(TParticle* part, Double_t magField)
1414{
e2ee5727 1415 // printf(" TrackHitsCalorimeter - magField %f\n", magField);
1416 AliESDtrack *esdTrack = new AliESDtrack(part);
1417 // Printf("MC Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
1418
1419 Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
1420
1421 // if(prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
1422
f61cec2f 1423 bool status = prop && fSelector->CutGeometricalAcceptance(*esdTrack);
e2ee5727 1424 delete esdTrack;
1425
1426 return status;
1427}
1428
1429void AliAnalysisEtMonteCarlo::FillHistograms()
f61cec2f 1430{ // let base class fill its histograms, and us fill the local ones
e2ee5727 1431 AliAnalysisEt::FillHistograms();
1432 //std::cout << fEtNonRemovedPiPlus << " " << fCentClass << std::endl;
1433
1434 fHistEtNonRemovedProtons->Fill(fEtNonRemovedProtons, fCentClass);
1435 fHistEtNonRemovedAntiProtons->Fill(fEtNonRemovedAntiProtons, fCentClass);
1436 fHistEtNonRemovedKaonPlus->Fill(fEtNonRemovedKaonPlus, fCentClass);
1437 fHistEtNonRemovedKaonMinus->Fill(fEtNonRemovedKaonMinus, fCentClass);
ef647350 1438 fHistEtNonRemovedK0s->Fill(fEtNonRemovedK0S, fCentClass);
1439 fHistEtNonRemovedK0L->Fill(fEtNonRemovedK0L, fCentClass);
e2ee5727 1440 fHistEtNonRemovedLambdas->Fill(fEtNonRemovedLambdas, fCentClass);
1441 fHistEtNonRemovedPiPlus->Fill(fEtNonRemovedPiPlus, fCentClass);
1442 fHistEtNonRemovedPiMinus->Fill(fEtNonRemovedPiMinus, fCentClass);
1443 fHistEtNonRemovedElectrons->Fill(fEtNonRemovedElectrons, fCentClass);
1444 fHistEtNonRemovedPositrons->Fill(fEtNonRemovedPositrons, fCentClass);
1445 fHistEtNonRemovedMuPlus->Fill(fEtNonRemovedMuPlus, fCentClass);
1446 fHistEtNonRemovedMuMinus->Fill(fEtNonRemovedMuMinus, fCentClass);
1447 fHistEtNonRemovedNeutrons->Fill(fEtNonRemovedNeutrons, fCentClass);
1448 fHistEtNonRemovedAntiNeutrons->Fill(fEtNonRemovedAntiNeutrons, fCentClass);
1449 fHistEtNonRemovedGammas->Fill(fEtNonRemovedGammas, fCentClass);
1450 fHistEtNonRemovedGammasFromPi0->Fill(fEtNonRemovedGammasFromPi0, fCentClass);
1451
ef647350 1452 fHistEtRemovedGammas->Fill(fEtRemovedGammas, fNClusters);
e2ee5727 1453 fHistEtRemovedNeutrons->Fill(fEtRemovedNeutrons, fCentClass);
1454 fHistEtRemovedAntiNeutrons->Fill(fEtRemovedAntiNeutrons, fCentClass);
1455
ef647350 1456// fHistEtRemovedCharged->Fill(fEtRemovedAntiProtons+fEtRemovedElectrons+fEtRemovedKaonMinus+fEtRemovedKaonPlus
1457// +fEtRemovedMuMinus+fEtRemovedMuPlus+fEtRemovedPiMinus+fEtRemovedPiPlus+fEtRemovedPositrons
1458// +fEtRemovedProtons.
1459// fCentClass);
1460// fHistEtRemovedNeutrals->Fill(fEtRemovedNeutrons+fEtRemovedAntiNeutrons, fCentClass);
f61cec2f 1461//
ef647350 1462// fHistEtNonRemovedCharged->Fill(fEtNonRemovedAntiProtons+fEtNonRemovedElectrons+fEtNonRemovedKaonMinus+fEtNonRemovedKaonPlus
1463// +fEtNonRemovedMuMinus+fEtNonRemovedMuPlus+fEtNonRemovedPiMinus+fEtNonRemovedPiPlus+fEtNonRemovedPositrons
1464// +fEtNonRemovedProtons,
1465// fCentClass);
1466// fHistEtRemovedNeutrals->Fill(fEtNonRemovedNeutrons+fEtNonRemovedAntiNeutrons, fCentClass);
1467
1468 fHistEtRemovedCharged->Fill(fEnergyChargedRemoved, fNClusters);
1469 fHistEtRemovedNeutrals->Fill(fEnergyNeutralRemoved, fNClusters);
1470 fHistEtNonRemovedCharged->Fill(fEnergyChargedNotRemoved, fNClusters);
1471 fHistEtNonRemovedNeutrals->Fill(fEnergyNeutralNotRemoved, fNClusters);
f61cec2f 1472
ef647350 1473 fHistMultRemovedCharged->Fill(fChargedRemoved, fNClusters);
1474 fHistMultRemovedNeutrals->Fill(fNeutralRemoved, fNClusters);
1475 fHistMultNonRemovedCharged->Fill(fChargedNotRemoved, fNClusters);
1476 fHistMultNonRemovedNeutrals->Fill(fNeutralNotRemoved, fNClusters);
f61cec2f 1477
1478
e2ee5727 1479 fHistMultNonRemovedProtons->Fill(fMultNonRemovedProtons, fCentClass);
1480 fHistMultNonRemovedAntiProtons->Fill(fMultNonRemovedAntiProtons, fCentClass);
1481 fHistMultNonRemovedKaonPlus->Fill(fMultNonRemovedKaonPlus, fCentClass);
1482 fHistMultNonRemovedKaonMinus->Fill(fMultNonRemovedKaonMinus, fCentClass);
1483 fHistMultNonRemovedK0s->Fill(fMultNonRemovedK0s, fCentClass);
ef647350 1484 fHistMultNonRemovedK0L->Fill(fMultNonRemovedK0L, fCentClass);
e2ee5727 1485 fHistMultNonRemovedLambdas->Fill(fMultNonRemovedLambdas, fCentClass);
1486 fHistMultNonRemovedPiPlus->Fill(fMultNonRemovedPiPlus, fCentClass);
1487 fHistMultNonRemovedPiMinus->Fill(fMultNonRemovedPiMinus, fCentClass);
1488 fHistMultNonRemovedElectrons->Fill(fMultNonRemovedElectrons, fCentClass);
1489 fHistMultNonRemovedPositrons->Fill(fMultNonRemovedPositrons, fCentClass);
1490 fHistMultNonRemovedMuPlus->Fill(fMultNonRemovedMuPlus, fCentClass);
1491 fHistMultNonRemovedMuMinus->Fill(fMultNonRemovedMuMinus, fCentClass);
1492 fHistMultNonRemovedNeutrons->Fill(fMultNonRemovedNeutrons, fCentClass);
1493 fHistMultNonRemovedAntiNeutrons->Fill(fMultNonRemovedAntiNeutrons, fCentClass);
1494 fHistMultNonRemovedGammas->Fill(fMultNonRemovedGammas, fCentClass);
1495
1496 fHistMultRemovedGammas->Fill(fMultRemovedGammas, fCentClass);
1497 fHistMultRemovedNeutrons->Fill(fMultRemovedNeutrons, fCentClass);
1498 fHistMultRemovedAntiNeutrons->Fill(fMultRemovedAntiNeutrons, fCentClass);
1499
1500 fHistTrackMultvsNonRemovedCharged->Fill(fTrackMultInAcc,
1501 fMultNonRemovedAntiProtons+fMultNonRemovedElectrons+fMultNonRemovedKaonMinus+fMultNonRemovedKaonPlus
1502 +fMultNonRemovedMuMinus+fMultNonRemovedMuPlus+fMultNonRemovedPiMinus+fMultNonRemovedPiPlus+fMultNonRemovedPositrons
1503 +fMultNonRemovedProtons);
1504
1505 fHistTrackMultvsNonRemovedNeutral->Fill(fTrackMultInAcc,
b2c10007 1506 fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedK0L+fMultNonRemovedLambdas+fK0sMult);
e2ee5727 1507
1508 fHistTrackMultvsRemovedGamma->Fill(fTrackMultInAcc,
1509 fMultRemovedGammas);
1510
ef647350 1511 fHistClusterMultvsNonRemovedCharged->Fill(fNClusters,
e2ee5727 1512 fMultNonRemovedAntiProtons+fMultNonRemovedElectrons+fMultNonRemovedKaonMinus
1513 +fMultNonRemovedKaonPlus+fMultNonRemovedMuMinus+fMultNonRemovedMuPlus
1514 +fMultNonRemovedPiMinus+fMultNonRemovedPiPlus+fMultNonRemovedPositrons+fMultNonRemovedProtons);
1515
ef647350 1516 fHistClusterMultvsNonRemovedNeutral->Fill(fNClusters,
b2c10007 1517 fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedK0L+fMultNonRemovedLambdas+fK0sMult);
e2ee5727 1518
ef647350 1519 fHistClusterMultvsRemovedGamma->Fill(fNClusters,
e2ee5727 1520 fMultRemovedGammas);
1521
13b0d3c1 1522}
1523
e2ee5727 1524
1525
1526
ef647350 1527Int_t AliAnalysisEtMonteCarlo::PrintFamilyTree(Int_t partIdx, AliStack* stack)
4d376d01 1528{ // print family tree
f61cec2f 1529 TParticle *part = stack->Particle(partIdx);
b2c10007 1530// if(part->GetPdgCode() == fgK0SCode)
f61cec2f 1531 {
1532 std::cout << "This is index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") , is it primary: " << stack->IsPhysicalPrimary(partIdx)<< std::endl;
1533 std::cout << "PID: " << part->GetPdgCode() << "/" << part->GetName() << std::endl;
1534 std::cout << "Energy: " << part->Energy() << std::endl;
1535 std::cout << "Vertex: " << part->Vx() << ", " << part->Vy() << ", " << part->Vz() << std::endl;
1536 }
1537 return PrintMothers(partIdx, stack, 1);
ef647350 1538}
1539
1540Int_t AliAnalysisEtMonteCarlo::PrintMothers(Int_t partIdx, AliStack* stack, Int_t gen)
4d376d01 1541{ // print mothers
f61cec2f 1542 char *tabs = new char[gen+1];
1543 for(Int_t i = 0; i < gen; ++i)
1544 {
1545 //std::cout << i << std::endl;
1546 tabs[i] = '\t';
1547 }
1548 tabs[gen] = '\0';
1549 Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
1550 if(mothIdx < 0)
1551 {
1552 delete [] tabs;
1553 return 0;
1554 }
1555 TParticle *mother = stack->Particle(mothIdx);
b2c10007 1556// if(mother->GetPdgCode() == fgK0SCode)
f61cec2f 1557 {
b2c10007 1558 //std::cout << tabs << "Mother of index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") is: " << mothIdx << ", is it primary: " << stack->IsPhysicalPrimary(mothIdx)<< std::endl;
1559 std::cout << tabs << "Index: " << mothIdx << std::endl;
1560 std::cout << tabs << "Primary: " << stack->IsPhysicalPrimary(mothIdx) << std::endl;
f61cec2f 1561 std::cout << tabs << "PID: " << mother->GetPdgCode() << "/" << mother->GetName() << std::endl;
1562 std::cout << tabs << "Energy: " << mother->Energy() << std::endl;
b2c10007 1563 if(mother->GetFirstMother() >= 0)
1564 {
1565 std::cout << tabs << "Mother(s): " << stack->Particle(mother->GetFirstMother())->GetPdgCode();
1566 if(mother->GetSecondMother() >= 0) std::cout << ", " << stack->Particle(mother->GetSecondMother())->GetPdgCode();
1567 std::cout << std::endl;
1568 }
f61cec2f 1569 std::cout << tabs << "Vertex: " << mother->Vx() << ", " << mother->Vy() << ", " << mother->Vz() << std::endl;
1570 }
1571 if(mother->GetPdgCode() == fgK0SCode)
1572 {
ef647350 1573// std::cout << "K0S!!!!!!!!!!!!!11111!!!!!" << std::endl;
f61cec2f 1574 }
ef647350 1575// std::cout << "Mother of index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") is: " << mothIdx << std::endl;
1576// std::cout << "PID: " << mother->GetPdgCode() << "/" << mother->GetName() << std::endl;
1577// std::cout << "Energy: " << mother->Energy() << std::endl;
1578// std::cout << "Vertex: " << mother->Vx() << ", " << mother->Vy() << ", " << mother->Vz() << std::endl;
f61cec2f 1579
1580 delete [] tabs;
1581 return PrintMothers(mothIdx, stack, gen+1) + 1;
ef647350 1582}
1583
1584Int_t AliAnalysisEtMonteCarlo::GetPrimMother(Int_t partIdx, AliStack *stack)
4d376d01 1585{ // get primary mother
f61cec2f 1586 if(partIdx >= 0)
ef647350 1587 {
b2c10007 1588 //return stack->GetPrimary(partIdx);
1589
f61cec2f 1590 Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
1591 if(mothIdx < 0) return -1;
1592 TParticle *mother = stack->Particle(mothIdx);
1593 if(mother)
1594 {
1595 // if(mother->GetPdgCode() == fgK0SCode)
1596 //{
ef647350 1597// std::cout << "!!!!!!!!!!!!!!!!! K0S !!!!!!!!!!!!!!!!!!" << std::endl;
f61cec2f 1598 //return mothIdx;
1599 // }
1600 //if(mother->GetPdgCode() == fgK0SCode&& stack->IsPhysicalPrimary(mothIdx))
1601 //{
ef647350 1602// std::cout << "!!!!!!!!!!!!!!!!! Primary K0S !!!!!!!!!!!!!!!!!!" << std::endl;
f61cec2f 1603 //return mothIdx;
1604 // }
1605 if(stack->IsPhysicalPrimary(mothIdx)) return mothIdx;
1606 else return GetPrimMother(mothIdx, stack);
1607 }
1608 else
1609 {
1610 return -1;
1611 }
ef647350 1612 }
f61cec2f 1613 return -1;
ef647350 1614}
1615
1616Int_t AliAnalysisEtMonteCarlo::GetK0InFamily(Int_t partIdx, AliStack* stack)
4d376d01 1617{ // get K0 in family
f61cec2f 1618 if(partIdx >= 0)
ef647350 1619 {
2aab9269 1620 if(stack->Particle(partIdx)->GetPdgCode() == fgK0SCode) return partIdx;
f61cec2f 1621 Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
1622 if(mothIdx < 0) return -1;
1623 TParticle *mother = stack->Particle(mothIdx);
1624 if(mother)
ef647350 1625 {
2aab9269 1626 if(mother->GetPdgCode() == fgK0SCode)
ef647350 1627 {
f61cec2f 1628 return mothIdx;
ef647350 1629 }
f61cec2f 1630 return GetK0InFamily(mothIdx, stack);
1631 }
1632 else
1633 {
1634 return -1;
ef647350 1635 }
1636 }
f61cec2f 1637 return -1;
ef647350 1638}
e2ee5727 1639