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