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