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