]> 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
cleanup
[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"
13b0d3c1 12#include "AliESDtrack.h"
2fbf38ac 13#include "AliStack.h"
0651f6b4 14#include "AliVEvent.h"
2fbf38ac 15#include "AliMCEvent.h"
0651f6b4 16#include "AliESDEvent.h"
13b0d3c1 17#include "TH2F.h"
e2ee5727 18#include "TH3F.h"
cf6522d1 19#include "TParticle.h"
ce546038 20#include "AliGenHijingEventHeader.h"
21#include "AliGenPythiaEventHeader.h"
0651f6b4 22#include "TList.h"
23#include "AliESDCaloCluster.h"
0f6416f3 24#include "AliLog.h"
e2ee5727 25#include <iostream>
26#include <AliCentrality.h>
ce546038 27
e2ee5727 28using namespace std;
16abb579 29
30ClassImp(AliAnalysisEtMonteCarlo);
31
32
ce546038 33// ctor
0651f6b4 34AliAnalysisEtMonteCarlo::AliAnalysisEtMonteCarlo():AliAnalysisEt()
e2ee5727 35 ,fImpactParameter(0)
36 ,fNcoll(0)
37 ,fNpart(0)
e2ee5727 38
39 ,fHistDecayVertexNonRemovedCharged(0)
40 ,fHistDecayVertexRemovedCharged(0)
41 ,fHistDecayVertexNonRemovedNeutral(0)
42 ,fHistDecayVertexRemovedNeutral(0)
43
44 ,fHistRemovedOrNot(0)
476828ae 45 ,fHistEtNonRemovedProtons(0)
46 ,fHistEtNonRemovedAntiProtons(0)
47 ,fHistEtNonRemovedPiPlus(0)
48 ,fHistEtNonRemovedPiMinus(0)
49 ,fHistEtNonRemovedKaonPlus(0)
50 ,fHistEtNonRemovedKaonMinus(0)
51 ,fHistEtNonRemovedK0s(0)
52 ,fHistEtNonRemovedLambdas(0)
53 ,fHistEtNonRemovedElectrons(0)
54 ,fHistEtNonRemovedPositrons(0)
55 ,fHistEtNonRemovedMuPlus(0)
56 ,fHistEtNonRemovedMuMinus(0)
57 ,fHistEtNonRemovedNeutrons(0)
58 ,fHistEtNonRemovedAntiNeutrons(0)
59 ,fHistEtNonRemovedGammas(0)
60 ,fHistEtNonRemovedGammasFromPi0(0)
61 ,fHistEtRemovedGammas(0)
62 ,fHistEtRemovedNeutrons(0)
63 ,fHistEtRemovedAntiNeutrons(0)
64 ,fHistMultNonRemovedProtons(0)
65 ,fHistMultNonRemovedAntiProtons(0)
66 ,fHistMultNonRemovedPiPlus(0)
67 ,fHistMultNonRemovedPiMinus(0)
68 ,fHistMultNonRemovedKaonPlus(0)
69 ,fHistMultNonRemovedKaonMinus(0)
70 ,fHistMultNonRemovedK0s(0)
71 ,fHistMultNonRemovedLambdas(0)
72 ,fHistMultNonRemovedElectrons(0)
73 ,fHistMultNonRemovedPositrons(0)
74 ,fHistMultNonRemovedMuPlus(0)
75 ,fHistMultNonRemovedMuMinus(0)
76 ,fHistMultNonRemovedNeutrons(0)
77 ,fHistMultNonRemovedAntiNeutrons(0)
78 ,fHistMultNonRemovedGammas(0)
79 ,fHistMultRemovedGammas(0)
80 ,fHistMultRemovedNeutrons(0)
81 ,fHistMultRemovedAntiNeutrons(0)
82 ,fHistTrackMultvsNonRemovedCharged(0)
83 ,fHistTrackMultvsNonRemovedNeutral(0)
84 ,fHistTrackMultvsRemovedGamma(0)
85 ,fHistClusterMultvsNonRemovedCharged(0)
86 ,fHistClusterMultvsNonRemovedNeutral(0)
87 ,fHistClusterMultvsRemovedGamma(0)
88 ,fHistMultvsNonRemovedChargedE(0)
89 ,fHistMultvsNonRemovedNeutralE(0)
90 ,fHistMultvsRemovedGammaE(0)
91 ,fEtNonRemovedProtons(0)
92 ,fEtNonRemovedAntiProtons(0)
93 ,fEtNonRemovedPiPlus(0)
94 ,fEtNonRemovedPiMinus(0)
95 ,fEtNonRemovedKaonPlus(0)
96 ,fEtNonRemovedKaonMinus(0)
97 ,fEtNonRemovedK0s(0)
98 ,fEtNonRemovedLambdas(0)
99 ,fEtNonRemovedElectrons(0)
100 ,fEtNonRemovedPositrons(0)
101 ,fEtNonRemovedMuMinus(0)
102 ,fEtNonRemovedMuPlus(0)
103 ,fEtNonRemovedGammas(0)
104 ,fEtNonRemovedGammasFromPi0(0)
105 ,fEtNonRemovedNeutrons(0)
106 ,fEtNonRemovedAntiNeutrons(0)
107 ,fEtRemovedGammas(0)
108 ,fEtRemovedNeutrons(0)
109 ,fEtRemovedAntiNeutrons(0)
110 ,fMultNonRemovedProtons(0)
111 ,fMultNonRemovedAntiProtons(0)
112 ,fMultNonRemovedPiPlus(0)
113 ,fMultNonRemovedPiMinus(0)
114 ,fMultNonRemovedKaonPlus(0)
115 ,fMultNonRemovedKaonMinus(0)
116 ,fMultNonRemovedK0s(0)
117 ,fMultNonRemovedLambdas(0)
118 ,fMultNonRemovedElectrons(0)
119 ,fMultNonRemovedPositrons(0)
120 ,fMultNonRemovedMuMinus(0)
121 ,fMultNonRemovedMuPlus(0)
122 ,fMultNonRemovedGammas(0)
123 ,fMultNonRemovedNeutrons(0)
124 ,fMultNonRemovedAntiNeutrons(0)
125 ,fMultRemovedGammas(0)
126 ,fMultRemovedNeutrons(0)
127 ,fMultRemovedAntiNeutrons(0)
128 ,fTrackMultInAcc(0)
e2ee5727 129 ,fHistDxDzNonRemovedCharged(0)
476828ae 130 ,fHistDxDzRemovedCharged(0)
e2ee5727 131 ,fHistDxDzNonRemovedNeutral(0)
476828ae 132 ,fHistDxDzRemovedNeutral(0)
133 ,fHistPiPlusMult(0)
134 ,fHistPiMinusMult(0)
135 ,fHistPiZeroMult(0)
136 ,fHistPiPlusMultAcc(0)
137 ,fHistPiMinusMultAcc(0)
138 ,fHistPiZeroMultAcc(0)
139 ,fPiPlusMult(0)
140 ,fPiMinusMult(0)
141 ,fPiZeroMult(0)
142 ,fPiPlusMultAcc(0)
143 ,fPiMinusMultAcc(0)
144 ,fPiZeroMultAcc(0)
e2ee5727 145 ,fNeutralRemoved(0)
146 ,fChargedRemoved(0)
147 ,fChargedNotRemoved(0)
148 ,fNeutralNotRemoved(0)
149 ,fEnergyNeutralRemoved(0)
150 ,fEnergyChargedRemoved(0)
151 ,fEnergyChargedNotRemoved(0)
152 ,fEnergyNeutralNotRemoved(0)
153
ce546038 154{
e2ee5727 155 fTrackDistanceCut = 7.0;
ce546038 156}
157
158// dtor
e2ee5727 159AliAnalysisEtMonteCarlo::~AliAnalysisEtMonteCarlo()
d0c22dcc 160{//Destructor
161 delete fHistDecayVertexNonRemovedCharged; // Decay vertex for non-removed charged particles
162 delete fHistDecayVertexRemovedCharged; // Decay vertex for non-removed charged particles
163 delete fHistDecayVertexNonRemovedNeutral; // Decay vertex for non-removed charged particles
164 delete fHistDecayVertexRemovedNeutral; // Decay vertex for non-removed charged particles
165
166 delete fHistRemovedOrNot; // If charged/neutral particles were removed or not
167
168 delete fHistEtNonRemovedProtons; // enter comment here
169 delete fHistEtNonRemovedAntiProtons; // enter comment here
170 delete fHistEtNonRemovedPiPlus; // enter comment here
171 delete fHistEtNonRemovedPiMinus; // enter comment here
172 delete fHistEtNonRemovedKaonPlus; // enter comment here
173 delete fHistEtNonRemovedKaonMinus; // enter comment here
174 delete fHistEtNonRemovedK0s; // enter comment here
175 delete fHistEtNonRemovedLambdas; // enter comment here
176 delete fHistEtNonRemovedElectrons; // enter comment here
177 delete fHistEtNonRemovedPositrons; // enter comment here
178 delete fHistEtNonRemovedMuPlus; // enter comment here
179 delete fHistEtNonRemovedMuMinus; // enter comment here
180 delete fHistEtNonRemovedNeutrons; // enter comment here
181 delete fHistEtNonRemovedAntiNeutrons; // enter comment here
182 delete fHistEtNonRemovedGammas; // enter comment here
183 delete fHistEtNonRemovedGammasFromPi0; // enter comment here
184
185 delete fHistEtRemovedGammas; // enter comment here
186 delete fHistEtRemovedNeutrons; // enter comment here
187 delete fHistEtRemovedAntiNeutrons; // enter comment here
188
189
190 delete fHistMultNonRemovedProtons; // enter comment here
191 delete fHistMultNonRemovedAntiProtons; // enter comment here
192 delete fHistMultNonRemovedPiPlus; // enter comment here
193 delete fHistMultNonRemovedPiMinus; // enter comment here
194 delete fHistMultNonRemovedKaonPlus; // enter comment here
195 delete fHistMultNonRemovedKaonMinus; // enter comment here
196 delete fHistMultNonRemovedK0s; // enter comment here
197 delete fHistMultNonRemovedLambdas; // enter comment here
198 delete fHistMultNonRemovedElectrons; // enter comment here
199 delete fHistMultNonRemovedPositrons; // enter comment here
200 delete fHistMultNonRemovedMuPlus; // enter comment here
201 delete fHistMultNonRemovedMuMinus; // enter comment here
202 delete fHistMultNonRemovedNeutrons; // enter comment here
203 delete fHistMultNonRemovedAntiNeutrons; // enter comment here
204 delete fHistMultNonRemovedGammas; // enter comment here
205
206 delete fHistMultRemovedGammas; // enter comment here
207 delete fHistMultRemovedNeutrons; // enter comment here
208 delete fHistMultRemovedAntiNeutrons; // enter comment here
209
210 delete fHistTrackMultvsNonRemovedCharged; // enter comment here
211 delete fHistTrackMultvsNonRemovedNeutral; // enter comment here
212 delete fHistTrackMultvsRemovedGamma; // enter comment here
213
214 delete fHistClusterMultvsNonRemovedCharged; // enter comment here
215 delete fHistClusterMultvsNonRemovedNeutral; // enter comment here
216 delete fHistClusterMultvsRemovedGamma; // enter comment here
217
218 delete fHistMultvsNonRemovedChargedE; // enter comment here
219 delete fHistMultvsNonRemovedNeutralE; // enter comment here
220 delete fHistMultvsRemovedGammaE; // enter comment here
221 delete fHistDxDzNonRemovedCharged; // enter comment here
222 delete fHistDxDzRemovedCharged; // enter comment here
223 delete fHistDxDzNonRemovedNeutral; // enter comment here
224 delete fHistDxDzRemovedNeutral; // enter comment here
225
226 delete fHistPiPlusMult; // enter comment here
227 delete fHistPiMinusMult; // enter comment here
228 delete fHistPiZeroMult; // enter comment here
229
230 delete fHistPiPlusMultAcc; // enter comment here
231 delete fHistPiMinusMultAcc; // enter comment here
232 delete fHistPiZeroMultAcc; // enter comment here
ce546038 233}
2fbf38ac 234
235Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
cf6522d1 236{ // analyse MC event
e2ee5727 237 ResetEventValues();
238
239 fPiPlusMult = 0;
240 fPiMinusMult = 0;
241 fPiZeroMult = 0;
242 if (fCentrality)
243 {
244 fCentClass = fCentrality->GetCentralityClass10(fCentralityMethod);
245
246 }
247 fSparseTracks[0] = 0;
248 fSparseTracks[1] = 0;
249 fSparseTracks[2] = 0;
250 fSparseTracks[3] = 0;
251 fSparseTracks[4] = 0;
252 fSparseTracks[5] = 0;
253 fSparseTracks[6] = fCentClass;
254
255
256 // Get us an mc event
257 if (!ev) {
258 AliFatal("ERROR: Event does not exist");
259 return 0;
260 }
261 AliMCEvent *event = dynamic_cast<AliMCEvent*>(ev);
262 if (!event) {
263 AliFatal("ERROR: MC Event does not exist");
264 return 0;
265 }
266
267 Double_t protonMass =fgProtonMass;
268
269 // Hijing header
270 AliGenEventHeader* genHeader = event->GenEventHeader();
271 if (!genHeader) {
272 Printf("ERROR: Event generation header does not exist");
273 return 0;
274 }
275 AliGenHijingEventHeader* hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
276 if (hijingGenHeader) {
277 fImpactParameter = hijingGenHeader->ImpactParameter();
278 fNcoll = hijingGenHeader->HardScatters(); // or should this be some combination of NN() NNw() NwN() NwNw() ?
279 fNpart = hijingGenHeader->ProjectileParticipants() + hijingGenHeader->TargetParticipants();
280 /*
281 printf("Hijing: ImpactParameter %g ReactionPlaneAngle %g \n",
282 hijingGenHeader->ImpactParameter(), hijingGenHeader->ReactionPlaneAngle());
283 printf("HardScatters %d ProjecileParticipants %d TargetParticipants %d\n",
284 hijingGenHeader->HardScatters(), hijingGenHeader->ProjectileParticipants(), hijingGenHeader->TargetParticipants());
285 printf("ProjSpectatorsn %d ProjSpectatorsp %d TargSpectatorsn %d TargSpectatorsp %d\n",
286 hijingGenHeader->ProjSpectatorsn(), hijingGenHeader->ProjSpectatorsp(), hijingGenHeader->TargSpectatorsn(), hijingGenHeader->TargSpectatorsp());
287 printf("NN %d NNw %d NwN %d, NwNw %d\n",
288 hijingGenHeader->NN(), hijingGenHeader->NNw(), hijingGenHeader->NwN(), hijingGenHeader->NwNw());
289 */
290 }
291
292 /* // placeholder if we want to get some Pythia info later
293 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
294 if (pythiaGenHeader) { // not Hijing; try with Pythia
295 printf("Pythia: ProcessType %d GetPtHard %g \n",
296 pythiaGenHeader->ProcessType(), pythiaGenHeader->GetPtHard());
297 }
ce546038 298 */
e2ee5727 299
300 // Let's play with the stack!
301 AliStack *stack = event->Stack();
302
303 Int_t nPrim = stack->GetNtrack();
304
305 Int_t partCount = 0;
306 for (Int_t iPart = 0; iPart < nPrim; iPart++)
2fbf38ac 307 {
e2ee5727 308
309 TParticle *part = stack->Particle(iPart);
b9462451 310// TParticle *partMom = 0;
311// TParticle *partMomLast = 0;
e2ee5727 312
313 if (!part)
314 {
315 Printf("ERROR: Could not get particle %d", iPart);
316 continue;
317 }
318
b9462451 319// Int_t iPartMom = part->GetMother(0);
320// Int_t iPartLastMom = part->GetMother(1);
e2ee5727 321
322 TParticlePDG *pdg = part->GetPDG(0);
00efe469 323 //TParticlePDG *pdgMom = 0;
b9462451 324// TParticlePDG *pdgMomLast = 0;
e2ee5727 325
326 if (!pdg)
2fbf38ac 327 {
e2ee5727 328 //Printf("ERROR: Could not get particle PDG %d", iPart);
329 continue;
2fbf38ac 330 }
331
b9462451 332// if (iPartMom>0)
333// {
334// partMom = stack->Particle(iPartMom);
335// //pdgMom = partMom->GetPDG(0);
336// }
337
338// if (iPartLastMom>0)
339// {
340// partMomLast = stack->Particle(iPartLastMom);
341// pdgMomLast = partMomLast->GetPDG(0);
342// }
e2ee5727 343
344
345 Double_t particleMassPart = 0; //The mass part in the Et calculation for this particle
346
347 // Check if it is a primary particle
348 //if (!stack->IsPhysicalPrimary(iPart)) continue;
349
e2ee5727 350 //printf("MC: iPart %03d eta %4.3f phi %4.3f code %d charge %g \n", iPart, part->Eta(), part->Phi(), pdg->PdgCode(), pdg->Charge()); // tmp/debug printout
351
352 // Check for reasonable (for now neutral and singly charged) charge on the particle
353 //TODO:Maybe not only singly charged?
354 if ((stack->IsPhysicalPrimary(iPart))&&( TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloSingleChargedParticle())<1e-3 || TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloNeutralParticle())<1e-3))
355 {
356
357 fMultiplicity++;
358
359 if (TMath::Abs(part->Eta()) < fCuts->GetCommonEtaCut())
2fbf38ac 360 {
e2ee5727 361 // calculate E_T
362 if (
363 TMath::Abs(pdg->PdgCode()) == fgProtonCode ||
364 TMath::Abs(pdg->PdgCode()) == fgNeutronCode ||
365 TMath::Abs(pdg->PdgCode()) == fgLambdaCode ||
366 TMath::Abs(pdg->PdgCode()) == fgXiCode ||
367 TMath::Abs(pdg->PdgCode()) == fgXi0Code ||
368 TMath::Abs(pdg->PdgCode()) == fgOmegaCode
369 )
370 {
371 if (pdg->PdgCode() > 0) {
372 particleMassPart = - protonMass;
373 }
374 if (pdg->PdgCode() < 0) {
375 particleMassPart = protonMass;
376 }
377 }
378 Double_t et = part->Energy() * TMath::Sin(part->Theta()) + particleMassPart;
379
380 fSparseTracks[0] = pdg->PdgCode();
381 fSparseTracks[1] = pdg->Charge()*3;
382 fSparseTracks[2] = pdg->Mass();
383 fSparseTracks[3] = et;
384 fSparseTracks[4] = part->Pt();
385 fSparseTracks[5] = part->Eta();
386 fSparseHistTracks->Fill(fSparseTracks);
387
388 // Fill up total E_T counters for each particle species
389 if (pdg->PdgCode() == fgProtonCode || pdg->PdgCode() == fgAntiProtonCode)
390 {
391 fProtonEt += et;
392 }
393 if (pdg->PdgCode() == fgPiPlusCode || pdg->PdgCode() == fgPiMinusCode)
394 {
395 fPionEt += et;
396 if (pdg->PdgCode() == fgPiPlusCode)
397 {
398 fPiPlusMult++;
399 }
400 else
401 {
402 fPiMinusMult++;
403 }
404 }
405 if (pdg->PdgCode() == fgGammaCode)
406 {
407 fPiZeroMult++;
408 }
409 if (pdg->PdgCode() == fgKPlusCode || pdg->PdgCode() == fgKMinusCode)
410 {
411 fChargedKaonEt += et;
412 }
413 if (pdg->PdgCode() == fgMuPlusCode || pdg->PdgCode() == fgMuMinusCode)
414 {
415 fMuonEt += et;
416 }
417 if (pdg->PdgCode() == fgEPlusCode || pdg->PdgCode() == fgEMinusCode)
418 {
419 fElectronEt += et;
420 }
421
422 // some neutrals also
423 if (pdg->PdgCode() == fgNeutronCode)
424 {
425 fNeutronEt += et;
426 }
427 if (pdg->PdgCode() == fgAntiNeutronCode)
428 {
429 fAntiNeutronEt += et;
430 }
431 if (pdg->PdgCode() == fgGammaCode)
432 {
433 fGammaEt += et;
434 }
435
436 // Neutral particles
437 if (TMath::Abs(pdg->Charge() - fCuts->GetMonteCarloNeutralParticle()) <1e-3 )
438 {
439
440 if (et > fCuts->GetCommonClusterEnergyCut()) fTotNeutralEt += et;
441
442 // inside EMCal acceptance
443 if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
444 {
445 fNeutralMultiplicity++;
446 //std::cout << pdg->PdgCode() << ", " << iPart << ", " << part->GetMother(0) << ", " << stack->IsPhysicalPrimary(iPart) << ", " << part->GetNDaughters() << ", " << part->Energy() << ", " << part->Eta() << ", " << part->Phi() << std::endl;
447
448 //if(part->GetDaughter(0) > 0) std::cout << stack->Particle(part->GetDaughter(0))->GetPdgCode() << " " << stack->Particle(part->GetDaughter(1))->GetPdgCode() << std::endl;
449 if (et > fCuts->GetCommonClusterEnergyCut()) fTotNeutralEtAcc += et;
450 if (et > fCuts->GetCommonClusterEnergyCut()) fTotEtAcc += et;
451 if (part->Energy() > 0.05) partCount++;
e2ee5727 452 }
453 }
454 //Charged particles
455 else if (TMath::Abs( pdg->Charge() - fCuts->GetMonteCarloNeutralParticle())>1e-3 )
456 {
457 fChargedMultiplicity++;
458 fTotChargedEt += et;
459
460 // inside EMCal acceptance
461 if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
462 {
463 fTotChargedEtAcc += et;
464 fTotEtAcc += et;
465
466 if (pdg->PdgCode() == fgProtonCode || pdg->PdgCode() == fgAntiProtonCode)
467 {
468 fProtonEtAcc += et;
469 }
470 if (pdg->PdgCode() == fgPiPlusCode || pdg->PdgCode() == fgPiMinusCode)
471 {
472 fPionEtAcc += et;
473 }
474 if (pdg->PdgCode() == fgKPlusCode || pdg->PdgCode() == fgKMinusCode)
475 {
476 fChargedKaonEtAcc += et;
477 }
478 if (pdg->PdgCode() == fgMuPlusCode || pdg->PdgCode() == fgMuMinusCode)
479 {
480 fMuonEtAcc += et;
481 }
e1fa1966 482
e2ee5727 483 if (pdg->PdgCode() == fgEPlusCode || pdg->PdgCode() == fgEMinusCode)
484 {
e1fa1966 485 fElectronEtAcc += et;
486 }
e2ee5727 487 } // inside EMCal acceptance
488
489 // if (TrackHitsCalorimeter(part, event->GetMagneticField()))
490 if (TrackHitsCalorimeter(part)) // magnetic field info not filled?
491 {
492 if (pdg->Charge() > 0) fHistPhivsPtPos->Fill(part->Phi(),part->Pt());
493 else fHistPhivsPtNeg->Fill(part->Phi(), part->Pt());
494 }
495 }
496 }
497
498 }
2fbf38ac 499 }
e2ee5727 500
501 fTotEt = fTotChargedEt + fTotNeutralEt;
502 fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;
503
504
505 //std::cout << "Event done! # of particles: " << partCount << std::endl;
506 return 0;
2fbf38ac 507}
0651f6b4 508//Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliMCEvent* mcEvent,AliESDEvent* realEvent)
509Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
510{ // analyse MC and real event info
e2ee5727 511 //if(!mcEvent || !realEvent){
512 if (!ev || !ev2) {
513 AliFatal("ERROR: Event does not exist");
514 return 0;
515 }
516
517 fSparseClusters[0] = 0;
518 fSparseClusters[1] = 0;
519 fSparseClusters[2] = 0;
520 fSparseClusters[3] = 0;
521 fSparseClusters[4] = 0;
522 fSparseClusters[5] = 0;
523 fSparseClusters[6] = 0;
524 fSparseClusters[7] = 0;
525 fSparseClusters[8] = 0;
526 fSparseClusters[9] = fCentClass;
527 fSparseClusters[10] = 0;
528
529 //AnalyseEvent(mcEvent);
530 AnalyseEvent(ev);
531 AliMCEvent *mcEvent = dynamic_cast<AliMCEvent*>(ev);
532 AliESDEvent *realEvent = dynamic_cast<AliESDEvent*>(ev2);
533 if (!mcEvent || !realEvent) {
534 AliFatal("ERROR: mcEvent or realEvent does not exist");
535 return 0;
536 }
537
538 AliStack *stack = mcEvent->Stack();
539
540 // get all detector clusters
541 TRefArray* caloClusters = new TRefArray();
542 if (fDetector == fCuts->GetDetectorEmcal()) realEvent->GetEMCALClusters( caloClusters );
543 else if (fDetector == fCuts->GetDetectorPhos()) realEvent->GetPHOSClusters( caloClusters );
544 else {
545 AliFatal("Detector ID has not been specified");
546 return -1;
547 }
548
549 // Track loop to fill a pT spectrum
550 for (Int_t iTracks = 0; iTracks < realEvent->GetNumberOfTracks(); iTracks++)
0651f6b4 551 {
e2ee5727 552 AliESDtrack* track = realEvent->GetTrack(iTracks);
553 if (!track)
0651f6b4 554 {
e2ee5727 555 Printf("ERROR: Could not receive track %d", iTracks);
556 continue;
0651f6b4 557 }
e2ee5727 558
559 if ( track->Phi() < fCuts->GetGeometryPhosPhiAccMaxCut() * TMath::Pi()/180. && track->Phi() > fCuts->GetGeometryPhosPhiAccMinCut() * TMath::Pi()/180. && TMath::Abs(track->Eta()) < fCuts->GetGeometryPhosEtaAccCut())
560 {
561 fTrackMultInAcc++;
562 }
563 }
564
565 Int_t nCluster = caloClusters->GetEntries();
566
567 // loop the clusters
568 for (int iCluster = 0; iCluster < nCluster; iCluster++ )
569 {
570 AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
e1fa1966 571 //Float_t caloE = caloCluster->E();
e2ee5727 572
573 UInt_t iPart = (UInt_t)TMath::Abs(caloCluster->GetLabel());
574 TParticle *part = stack->Particle(iPart);
575 TParticle *partMom = 0;
576
577 if (!part)
0651f6b4 578 {
e2ee5727 579 Printf("No MC particle %d", iCluster);
580 continue;
581 }
582
583 // compare MC and Rec energies for all particles
e1fa1966 584 //fHistAllERecEMC->Fill(part->Energy(),caloE);
e2ee5727 585
586 TParticlePDG *pdg = part->GetPDG(0);
587 TParticlePDG *pdgMom = 0;
588
589 if (!pdg)
590 {
591 Printf("ERROR: Could not get particle PDG %d", iPart);
592 continue;
593 }
594
595 Int_t iPartMom = part->GetMother(0);
596
597 if (iPartMom>0)
598 {
599 partMom = stack->Particle(iPartMom);
600 pdgMom = partMom->GetPDG(0);
601 }
602
603 // Check if it is a primary particle
604 //if (!stack->IsPhysicalPrimary(iPart)) continue;
605 if (!stack->IsPhysicalPrimary(iPart)) // check whether particle is primary. we keep secondary electron and gamma for testing.
606 {
607 if (!((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode) || (pdg->PdgCode() == fgGammaCode)))
608 continue;
609 } // end of primary particle check
610
611 if (TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloSingleChargedParticle())<1e-3 && TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloNeutralParticle())<1e-3) continue;
612
e2ee5727 613 Double_t clEt = CalculateTransverseEnergy(caloCluster);
614 if (caloCluster->E() < fCuts->GetCommonClusterEnergyCut()) continue;
615 Float_t pos[3];
616
617 caloCluster->GetPosition(pos);
618 TVector3 cp(pos);
619
620 fSparseClusters[0] = pdg->PdgCode();
621 fSparseClusters[1] = pdg->Charge()/3;
622 fSparseClusters[2] = pdg->Mass();
623 fSparseClusters[3] = clEt;
624 fSparseClusters[4] = caloCluster->E();
625 fSparseClusters[5] = cp.Eta();
626 fSparseClusters[6] = part->Energy() * TMath::Sin(part->Theta());
627 fSparseClusters[7] = part->Pt();
628 fSparseClusters[8] = part->Eta();
629 fSparseClusters[9] = fCentClass;
630 fSparseClusters[10] = 0;
631 fSparseHistClusters->Fill(fSparseClusters);
632
633 //if(caloCluster->GetEmcCpvDistance() < fTrackDistanceCut)
634
635 if (TMath::Abs(caloCluster->GetTrackDx()) < fTrackDxCut && TMath::Abs(caloCluster->GetTrackDz()) < fTrackDzCut)
636 {
637
638 if (pdg->Charge() != 0)
639 {
640 //std::cout << "Removing charged: " << pdg->PdgCode() << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
641 fChargedRemoved++;
642 fEnergyChargedRemoved += caloCluster->E();
643 fHistRemovedOrNot->Fill(0.0, fCentClass);
644 //std::cout << part->Vx() << " " << part->Vy() << " " << part->Vz() << " " << std::endl;
645 //fHistDecayVertexRemovedCharged->Fill(part->Vy(), part->Vx(), part->Vz());
646 fHistDxDzRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
647
648
649 }
650 else
651 {
652 //std::cout << "Removing neutral: " << pdg->PdgCode() << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
653 //std::cout << "Mother is: " << stack->Particle(part->GetMother(0))->GetPdgCode() << std::endl;
654 fNeutralRemoved++;
655 fEnergyNeutralRemoved += caloCluster->E();
656 fHistRemovedOrNot->Fill(1.0, fCentClass);
657 //std::cout << part->Vx() << " " << part->Vy() << " " << part->Vz() << " " << std::endl;
658 //fHistDecayVertexRemovedNeutral->Fill(part->Vy(), part->Vx(), part->Vz());
659 fHistDxDzRemovedNeutral->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
660
661 if (pdg->PdgCode() == fgGammaCode)
662 {
663 fEtRemovedGammas += clEt;
664 fMultRemovedGammas++;
665 }
666 else if (pdg->PdgCode() == fgNeutronCode)
667 {
668 fEtRemovedNeutrons += clEt;
669 fMultRemovedNeutrons++;
670 }
671 else if (pdg->PdgCode() == fgAntiNeutronCode)
672 {
673 fEtRemovedAntiNeutrons += clEt;
674 fMultRemovedAntiNeutrons++;
675 }
676 //else std::cout << "Hmm, what is this (neutral: " << pdg->PdgCode() << std::endl;
677 }
678 }
679 else
680 {
681
682 if (pdg->Charge() != 0)
683 {
684 //std::cout << "Not removing charged: " << pdg->PdgCode() << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
685 //std::cout << "Mother is: " << stack->Particle(part->GetMother(0))->GetPdgCode() << std::endl;
686 fChargedNotRemoved++;
687 fEnergyChargedNotRemoved += caloCluster->E();
688 fHistRemovedOrNot->Fill(2.0, fCentClass);
689 //std::cout << fHistDecayVertexNonRemovedCharged << std::endl;
690 //std::cout << part->Vx() << " " << part->Vy() << " " << part->Vz() << " " << std::endl;
691 //fHistDecayVertexNonRemovedCharged->Fill(part->Vy(), part->Vx(), part->Vz());
692 fHistDxDzNonRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
693 if (pdg->PdgCode() == fgProtonCode)
694 {
695 //std::cout << clEt << std::endl;
696 fEtNonRemovedProtons += clEt;
697 fMultNonRemovedProtons++;
698 }
699 else if (pdg->PdgCode() == fgAntiProtonCode)
700 {
701 //std::cout << clEt << std::endl;
702 fEtNonRemovedAntiProtons += clEt;
703 fMultNonRemovedAntiProtons++;
704 }
705 else if (pdg->PdgCode() == fgPiPlusCode)
706 {
707 //std::cout << "PI+" << clEt << std::endl;
708 fEtNonRemovedPiPlus += clEt;
709 fMultNonRemovedPiPlus++;
710 }
711 else if (pdg->PdgCode() == fgPiMinusCode)
712 {
713 // std::cout << "PI-" << clEt << std::endl;
714 fEtNonRemovedPiMinus += clEt;
715 fMultNonRemovedPiMinus++;
716 }
717 else if (pdg->PdgCode() == fgKPlusCode)
718 {
719 //std::cout << clEt << std::endl;
720 fEtNonRemovedKaonPlus += clEt;
721 fMultNonRemovedKaonPlus++;
722 }
723 else if (pdg->PdgCode() == fgKMinusCode)
724 {
725 //std::cout << clEt << std::endl;
726 fEtNonRemovedKaonMinus += clEt;
727 fMultNonRemovedKaonMinus++;
728 }
729 else if (pdg->PdgCode() == fgEPlusCode)
730 {
731 //std::cout << clEt << std::endl;
732 if (TMath::Sqrt(part->Vx()*part->Vx() + part->Vy()*part->Vy()) < 430)
733 {
734 fEtNonRemovedPositrons += clEt;
735 fMultNonRemovedPositrons++;
736 }
737 }
738 else if (pdg->PdgCode() == fgEMinusCode)
739 {
740 //std::cout << clEt << std::endl;
741 if (TMath::Sqrt(part->Vx()*part->Vx() + part->Vy()*part->Vy()) < 430)
742 {
743 fEtNonRemovedElectrons += clEt;
744 fMultNonRemovedElectrons++;
745 }
746 }
747 else if (pdg->PdgCode() == fgMuPlusCode)
748 {
749 std::cout << clEt << std::endl;
750 fEtNonRemovedMuPlus += clEt;
751 fMultNonRemovedMuPlus++;
752 }
753 else if (pdg->PdgCode() == fgMuMinusCode)
754 {
755 std::cout << clEt << std::endl;
756 fEtNonRemovedMuMinus += clEt;
757 fMultNonRemovedMuMinus++;
758 }
759
760 }
761 else
762 {
763 //std::cout << "Accepted: " << pdg->PdgCode() << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
764 //std::cout << "Not removing charged: " << pdg->PdgCode() << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
476828ae 765 //std::cout << "Mother is: " << stack->Particle(part->GetMother(0))->GetPdgCode() << stack->Particle(part->GetMother(0))->GetPDG()->GetName() << ", E: " << part->Energy() << std::endl;
e2ee5727 766 fNeutralNotRemoved++;
767 fEnergyNeutralNotRemoved += caloCluster->E();
768 fHistRemovedOrNot->Fill(3.0, fCentClass);
769 //std::cout << part->Vx() << " " << part->Vy() << " " << part->Vz() << " " << std::endl;
770 //fHistDecayVertexNonRemovedNeutral->Fill(part->Vy(), part->Vx(), part->Vz());
771 fHistDxDzNonRemovedNeutral->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
772 if (pdg->PdgCode() == fgGammaCode)
773 {
774 fEtNonRemovedGammas += clEt;
775 fMultNonRemovedGammas++;
776 if (pdgMom)
777 {
778 if (TMath::Abs(pdgMom->PdgCode()) == fgPi0Code || TMath::Abs(pdgMom->PdgCode()) == fgEtaCode || TMath::Abs(pdgMom->PdgCode()) == 331)
779 {
780// std::cout << "Mother of gamma: " << pdgMom->PdgCode() << " " << pdgMom->GetName() << ", E: " << part->Energy() << std::endl;
781 fEtNonRemovedGammasFromPi0 += clEt;
782 }
783 }
784 }
785 else if (pdg->PdgCode() == fgNeutronCode)
786 {
787 fEtNonRemovedNeutrons += clEt;
788 fMultNonRemovedNeutrons++;
789 }
790 else if (pdg->PdgCode() == fgAntiNeutronCode)
791 {
792 fEtNonRemovedAntiNeutrons += clEt;
793 fMultNonRemovedAntiNeutrons++;
794 }
795 else if (pdg->PdgCode() == fgK0LCode || pdg->PdgCode() == fgK0SCode)
796 {
797 fEtNonRemovedK0s += clEt;
798 fMultNonRemovedK0s++;
799
800 }
801 else if (TMath::Abs(pdg->PdgCode()) == fgLambdaCode)
802 {
803 fEtNonRemovedLambdas += clEt;
804 fMultNonRemovedLambdas++;
805 }
806 else std::cout << "Hmm, what is this (neutral not removed): " << pdg->PdgCode() << " " << pdg->GetName() << ", ET: " << clEt << std::endl;
807 }
808 }
809 } // end of loop over clusters
810
811 delete caloClusters;
812 //std::cout << "Distance cut: " << fTrackDistanceCut << std::endl;
813 //std::cout << "Number of removed neutrals: " << fNeutralRemoved << std::endl;
814 //std::cout << "Number of removed charged: " << fChargedRemoved << std::endl;
815 //std::cout << "Number of non-removed charged: " << fChargedNotRemoved << std::endl;
816 //std::cout << "Number of non-removed neutral: " << fNeutralNotRemoved << std::endl;
817
818// std::cout << "Energy of removed neutrals: " << fEnergyNeutralRemoved << std::endl;
819// std::cout << "Energy of removed charged: " << fEnergyChargedRemoved << std::endl;
820// std::cout << "Energy of non-removed charged: " << fEnergyChargedNotRemoved << std::endl;
821// std::cout << "Energy of non-removed neutral: " << fEnergyNeutralNotRemoved << std::endl;
822
823 FillHistograms();
824 return 0;
825}
2fbf38ac 826
827void AliAnalysisEtMonteCarlo::Init()
ce546038 828{ // init
e2ee5727 829 AliAnalysisEt::Init();
2fbf38ac 830}
13b0d3c1 831
ce546038 832void AliAnalysisEtMonteCarlo::ResetEventValues()
833{ // reset event values
e2ee5727 834 AliAnalysisEt::ResetEventValues();
835
836 // collision geometry defaults for p+p:
837 fImpactParameter = 0;
838 fNcoll = 1;
839 fNpart = 2;
840
841 fEtNonRemovedProtons = 0;
842 fEtNonRemovedAntiProtons = 0;
843 fEtNonRemovedPiPlus = 0;
844 fEtNonRemovedPiMinus = 0;
845 fEtNonRemovedKaonPlus = 0;
846 fEtNonRemovedKaonMinus = 0;
847 fEtNonRemovedK0s = 0;
848 fEtNonRemovedLambdas = 0;
849 fEtNonRemovedElectrons = 0;
850 fEtNonRemovedPositrons = 0;
851 fEtNonRemovedMuPlus = 0;
852 fEtNonRemovedMuMinus = 0;
853 fEtNonRemovedNeutrons = 0;
854 fEtNonRemovedAntiNeutrons = 0;
855 fEtNonRemovedGammas = 0;
856 fEtNonRemovedGammasFromPi0 = 0;
857
858 fEtRemovedGammas = 0;
859 fEtRemovedNeutrons = 0;
860 fEtRemovedAntiNeutrons = 0;
861
862 fMultNonRemovedProtons = 0;
863 fMultNonRemovedAntiProtons = 0;
864 fMultNonRemovedPiPlus = 0;
865 fMultNonRemovedPiMinus = 0;
866 fMultNonRemovedKaonPlus = 0;
867 fMultNonRemovedKaonMinus = 0;
868 fMultNonRemovedK0s = 0;
869 fMultNonRemovedLambdas = 0;
870 fMultNonRemovedElectrons = 0;
871 fMultNonRemovedPositrons = 0;
872 fMultNonRemovedMuPlus = 0;
873 fMultNonRemovedMuMinus = 0;
874 fMultNonRemovedNeutrons = 0;
875 fMultNonRemovedAntiNeutrons = 0;
876 fMultNonRemovedGammas = 0;
877
878 fMultRemovedGammas = 0;
879 fMultRemovedNeutrons = 0;
880 fMultRemovedAntiNeutrons = 0;
881
882 fTrackMultInAcc = 0;
883
ce546038 884}
885
886void AliAnalysisEtMonteCarlo::CreateHistograms()
887{ // histogram related additions
e2ee5727 888 AliAnalysisEt::CreateHistograms();
889 if (fTree) {
890 fTree->Branch("fImpactParameter",&fImpactParameter,"fImpactParameter/D");
891 fTree->Branch("fNcoll",&fNcoll,"fNcoll/I");
892 fTree->Branch("fNpart",&fNpart,"fNpart/I");
893 }
e1fa1966 894
e2ee5727 895 //fHistDecayVertexNonRemovedCharged = new TH3F("fHistDecayVertexNonRemovedCharged","fHistDecayVertexNonRemovedCharged", 500, -470, 30, 500, -300, 300, 40, -20, 20);
896 //fHistDecayVertexRemovedCharged = new TH3F("fHistDecayVertexRemovedCharged","fHistDecayVertexRemovedCharged", 500, -470, 30, 500, -300, 300, 40, -20, 20);
897 //fHistDecayVertexNonRemovedNeutral = new TH3F("fHistDecayVertexNonRemovedNeutral","fHistDecayVertexNonRemovedNeutral", 500, -470, 30, 500, -300, 300, 40, -20, 20);
898 //fHistDecayVertexRemovedNeutral = new TH3F("fHistDecayVertexRemovedNeutral","fHistDecayVertexRemovedNeutral", 500, -470, 30, 500, -300, 300, 40, -20, 20);
899
900 fHistRemovedOrNot = new TH2F("fHistRemovedOrNot", "fHistRemovedOrNot", 4, -0.5, 3.5, 10, -0.5, 9.5);
901
902 fHistEtNonRemovedProtons = new TH2F("fHistEtNonRemovedProtons", "fHistEtNonRemovedProtons", 1500, 0, 30, 10, -0.5, 9.5);
903 fHistEtNonRemovedAntiProtons = new TH2F("fHistEtNonRemovedAntiProtons", "fHistEtNonRemovedAntiProtons", 1500, 0, 30, 10, -0.5, 9.5);
904 fHistEtNonRemovedPiPlus = new TH2F("fHistEtNonRemovedPiPlus", "fHistEtNonRemovedPiPlus", 1500, 0, 30, 10, -0.5, 9.5);
905 fHistEtNonRemovedPiMinus = new TH2F("fHistEtNonRemovedPiMinus", "fHistEtNonRemovedPiMinus", 1500, 0, 30, 10, -0.5, 9.5);
906 fHistEtNonRemovedKaonPlus = new TH2F("fHistEtNonRemovedKaonPlus", "fHistEtNonRemovedKaonPlus", 1500, 0, 30, 10, -0.5, 9.5);
907 fHistEtNonRemovedKaonMinus = new TH2F("fHistEtNonRemovedKaonMinus", "fHistEtNonRemovedKaonMinus", 1500, 0, 30, 10, -0.5, 9.5);
908 fHistEtNonRemovedK0s = new TH2F("fHistEtNonRemovedK0s", "fHistEtNonRemovedK0s", 1500, 0, 30, 10, -0.5, 9.5);
909 fHistEtNonRemovedLambdas = new TH2F("fHistEtNonRemovedLambdas", "fHistEtNonRemovedLambdas", 1500, 0, 30, 10, -0.5, 9.5);
910 fHistEtNonRemovedElectrons = new TH2F("fHistEtNonRemovedElectrons", "fHistEtNonRemovedElectrons", 1500, 0, 30, 10, -0.5, 9.5);
911 fHistEtNonRemovedPositrons = new TH2F("fHistEtNonRemovedPositrons", "fHistEtNonRemovedPositrons", 1500, 0, 30, 10, -0.5, 9.5);
912 fHistEtNonRemovedMuPlus = new TH2F("fHistEtNonRemovedMuPlus", "fHistEtNonRemovedMuPlus", 1500, 0, 30, 10, -0.5, 9.5);
913 fHistEtNonRemovedMuMinus = new TH2F("fHistEtNonRemovedMuMinus", "fHistEtNonRemovedMuMinus", 1500, 0, 30, 10, -0.5, 9.5);
914 fHistEtNonRemovedNeutrons = new TH2F("fHistEtNonRemovedNeutrons", "fHistEtNonRemovedNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
915 fHistEtNonRemovedAntiNeutrons = new TH2F("fHistEtNonRemovedAntiNeutrons", "fHistEtNonRemovedAntiNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
916
917 fHistEtNonRemovedGammas = new TH2F("fHistEtNonRemovedGammas", "fHistEtNonRemovedGammas", 1500, 0, 30, 10, -0.5, 9.5);
918 fHistEtNonRemovedGammasFromPi0 = new TH2F("fHistEtNonRemovedGammasFromPi0", "fHistEtNonRemovedGammasFromPi0", 1500, 0, 30, 10, -0.5, 9.5);
919
920 fHistEtRemovedGammas = new TH2F("fHistEtRemovedGammas", "fHistEtRemovedGammas", 1500, 0, 30, 10, -0.5, 9.5);
921 fHistEtRemovedNeutrons = new TH2F("fHistEtRemovedNeutrons", "fHistEtRemovedNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
922 fHistEtRemovedAntiNeutrons = new TH2F("fHistEtRemovedAntiNeutrons", "fHistEtRemovedAntiNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
923
924 fHistMultNonRemovedProtons = new TH2F("fHistMultNonRemovedProtons", "fHistMultNonRemovedProtons", 100, -0.5, 99.5, 10, -0.5, 9.5);
925 fHistMultNonRemovedAntiProtons = new TH2F("fHistMultNonRemovedAntiProtons", "fHistMultNonRemovedAntiProtons", 100, -0.5, 99.5, 10, -0.5, 9.5);
926 fHistMultNonRemovedPiPlus = new TH2F("fHistMultNonRemovedPiPlus", "fHistMultNonRemovedPiPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
927 fHistMultNonRemovedPiMinus = new TH2F("fHistMultNonRemovedPiMinus", "fHistMultNonRemovedPiMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
928 fHistMultNonRemovedKaonPlus = new TH2F("fHistMultNonRemovedKaonPlus", "fHistMultNonRemovedKaonPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
929 fHistMultNonRemovedKaonMinus = new TH2F("fHistMultNonRemovedKaonMinus", "fHistMultNonRemovedKaonMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
930 fHistMultNonRemovedK0s = new TH2F("fHistMultNonRemovedK0s", "fHistMultNonRemovedK0s", 100, -0.5, 99.5, 10, -0.5, 9.5);
931 fHistMultNonRemovedLambdas = new TH2F("fHistMultNonRemovedLambdas", "fHistMultNonRemovedLambdas", 100, -0.5, 99.5, 10, -0.5, 9.5);
932 fHistMultNonRemovedElectrons = new TH2F("fHistMultNonRemovedElectrons", "fHistMultNonRemovedElectrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
933 fHistMultNonRemovedPositrons = new TH2F("fHistMultNonRemovedPositrons", "fHistMultNonRemovedPositrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
934 fHistMultNonRemovedMuPlus = new TH2F("fHistMultNonRemovedMuPlus", "fHistMultNonRemovedMuPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
935 fHistMultNonRemovedMuMinus = new TH2F("fHistMultNonRemovedMuMinus", "fHistMultNonRemovedMuMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
936 fHistMultNonRemovedNeutrons = new TH2F("fHistMultNonRemovedNeutrons", "fHistMultNonRemovedNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
937 fHistMultNonRemovedAntiNeutrons = new TH2F("fHistMultNonRemovedAntiNeutrons", "fHistMultNonRemovedAntiNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
938
939 fHistMultNonRemovedGammas = new TH2F("fHistMultNonRemovedGammas", "fHistMultNonRemovedGammas", 100, -0.5, 99.5, 10, -0.5, 9.5);
940
941 fHistMultRemovedGammas = new TH2F("fHistMultRemovedGammas", "fHistMultRemovedGammas", 100, -0.5, 99.5, 10, -0.5, 9.5);
942 fHistMultRemovedNeutrons = new TH2F("fHistMultRemovedNeutrons", "fHistMultRemovedNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
943 fHistMultRemovedAntiNeutrons = new TH2F("fHistMultRemovedAntiNeutrons", "fHistMultRemovedAntiNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
944
945 fHistTrackMultvsNonRemovedCharged = new TH2F("fHistTrackMultvsNonRemovedCharged", "fHistTrackMultvsNonRemovedCharged", 1000, -0.5, 999.5, 100, -0.5, 99.5);
946 fHistTrackMultvsNonRemovedNeutral = new TH2F("fHistTrackMultvsNonRemovedNeutral", "fHistTrackMultvsNonRemovedNeutral", 1000, -0.5, 999.5, 100, -0.5, 99.5);
947 fHistTrackMultvsRemovedGamma = new TH2F("fHistTrackMultvsRemovedGamma", "fHistTrackMultvsRemovedGamma", 1000, -0.5, 999.5, 100, -0.5, 99.5);
948
949 fHistClusterMultvsNonRemovedCharged = new TH2F("fHistClusterMultvsNonRemovedCharged", "fHistClusterMultvsNonRemovedCharged", 1000, -0.5, 999.5, 100, -0.5, 99.5);
950 fHistClusterMultvsNonRemovedNeutral = new TH2F("fHistClusterMultvsNonRemovedNeutral", "fHistClusterMultvsNonRemovedNeutral", 1000, -0.5, 999.5, 100, -0.5, 99.5);
951 fHistClusterMultvsRemovedGamma = new TH2F("fHistClusterMultvsRemovedGamma", "fHistClusterMultvsRemovedGamma", 1000, -0.5, 999.5, 100, -0.5, 99.5);
952
953 fHistDxDzNonRemovedCharged = new TH2F("fHistDxDzNonRemovedCharged", "fHistDxDzNonRemovedCharged", 800, -200, 200, 800, -200, 200);
954 fHistDxDzRemovedCharged = new TH2F("fHistDxDzRemovedCharged", "fHistDxDzRemovedCharged", 800, -200, 200, 800, -200, 200);
955 fHistDxDzNonRemovedNeutral = new TH2F("fHistDxDzNonRemovedNeutral", "fHistDxDzNonRemovedNeutral", 800, -200, 200, 800, -200, 200);
956 fHistDxDzRemovedNeutral = new TH2F("fHistDxDzRemovedNeutral", "fHistDxDzRemovedNeutral", 800, -200, 200, 800, -200, 200);
957
958 fHistPiPlusMult = new TH1F("fHistPiPlusMult", "fHistPiPlusMult", 2000, -0.5, 1999.5);
959 fHistPiMinusMult = new TH1F("fHistPiMinusMult", "fHistPiMinusMult", 2000, -0.5, 1999.5);
960 fHistPiZeroMult = new TH1F("fHistPiZeroMult", "fHistPiZeroMult", 2000, -0.5, 1999.5);
961
962 fHistPiPlusMultAcc = new TH1F("fHistPiPlusMultAcc", "fHistPiPlusMultAcc", 2000, -0.5, 1999.5);
963 fHistPiMinusMultAcc = new TH1F("fHistPiMinusMultAcc", "fHistPiMinusMultAcc", 2000, -0.5, 1999.5);
964 fHistPiZeroMultAcc = new TH1F("fHistPiZeroMultAcc", "fHistPiZeroMultAcc", 2000, -0.5, 1999.5);
965
ce546038 966}
967
0651f6b4 968void AliAnalysisEtMonteCarlo::FillOutputList(TList *list)
969{//fill the output list
e2ee5727 970 AliAnalysisEt::FillOutputList(list);
971
e2ee5727 972 list->Add(fHistRemovedOrNot);
973
974 list->Add(fHistEtNonRemovedProtons);
975 list->Add(fHistEtNonRemovedAntiProtons);
976 list->Add(fHistEtNonRemovedPiPlus);
977 list->Add(fHistEtNonRemovedPiMinus);
978 list->Add(fHistEtNonRemovedKaonPlus);
979 list->Add(fHistEtNonRemovedKaonMinus);
980 list->Add(fHistEtNonRemovedK0s);
981 list->Add(fHistEtNonRemovedLambdas);
982 list->Add(fHistEtNonRemovedElectrons);
983 list->Add(fHistEtNonRemovedPositrons);
984 list->Add(fHistEtNonRemovedMuPlus);
985 list->Add(fHistEtNonRemovedMuMinus);
986 list->Add(fHistEtNonRemovedNeutrons);
987 list->Add(fHistEtNonRemovedAntiNeutrons);
988 list->Add(fHistEtNonRemovedGammas);
989 list->Add(fHistEtNonRemovedGammasFromPi0);
990
991 list->Add(fHistEtRemovedGammas);
992 list->Add(fHistEtRemovedNeutrons);
993 list->Add(fHistEtRemovedAntiNeutrons);
994
995
996 list->Add(fHistMultNonRemovedProtons);
997 list->Add(fHistMultNonRemovedAntiProtons);
998 list->Add(fHistMultNonRemovedPiPlus);
999 list->Add(fHistMultNonRemovedPiMinus);
1000 list->Add(fHistMultNonRemovedKaonPlus);
1001 list->Add(fHistMultNonRemovedKaonMinus);
1002 list->Add(fHistMultNonRemovedK0s);
1003 list->Add(fHistMultNonRemovedLambdas);
1004 list->Add(fHistMultNonRemovedElectrons);
1005 list->Add(fHistMultNonRemovedPositrons);
1006 list->Add(fHistMultNonRemovedMuPlus);
1007 list->Add(fHistMultNonRemovedMuMinus);
1008 list->Add(fHistMultNonRemovedNeutrons);
1009 list->Add(fHistMultNonRemovedAntiNeutrons);
1010 list->Add(fHistMultNonRemovedGammas);
1011
1012 list->Add(fHistMultRemovedGammas);
1013 list->Add(fHistMultRemovedNeutrons);
1014 list->Add(fHistMultRemovedAntiNeutrons);
1015
1016 list->Add(fHistTrackMultvsNonRemovedCharged);
1017 list->Add(fHistTrackMultvsNonRemovedNeutral);
1018 list->Add(fHistTrackMultvsRemovedGamma);
1019
1020 list->Add(fHistClusterMultvsNonRemovedCharged);
1021 list->Add(fHistClusterMultvsNonRemovedNeutral);
1022 list->Add(fHistClusterMultvsRemovedGamma);
1023
1024 //list->Add(fHistDecayVertexNonRemovedCharged);
1025 //list->Add(fHistDecayVertexNonRemovedNeutral);
1026 //list->Add(fHistDecayVertexRemovedCharged);
1027 //list->Add(fHistDecayVertexRemovedNeutral);
1028
1029 list->Add(fHistDxDzNonRemovedCharged);
1030 list->Add(fHistDxDzRemovedCharged);
1031 list->Add(fHistDxDzNonRemovedNeutral);
1032 list->Add(fHistDxDzRemovedNeutral);
1033
1034 list->Add(fHistPiPlusMult);
1035 list->Add(fHistPiMinusMult);
1036 list->Add(fHistPiZeroMult);
1037 list->Add(fHistPiPlusMultAcc);
1038 list->Add(fHistPiMinusMultAcc);
1039 list->Add(fHistPiZeroMultAcc);
1040
0651f6b4 1041}
1042
1043
13b0d3c1 1044bool AliAnalysisEtMonteCarlo::TrackHitsCalorimeter(TParticle* part, Double_t magField)
1045{
e2ee5727 1046 // printf(" TrackHitsCalorimeter - magField %f\n", magField);
1047 AliESDtrack *esdTrack = new AliESDtrack(part);
1048 // Printf("MC Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
1049
1050 Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
1051
1052 // if(prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
1053
1054 bool status = prop &&
1055 TMath::Abs(esdTrack->Eta()) < fEtaCutAcc &&
1056 esdTrack->Phi() > fPhiCutAccMin*TMath::Pi()/180. &&
1057 esdTrack->Phi() < fPhiCutAccMax*TMath::Pi()/180.;
1058 delete esdTrack;
1059
1060 return status;
1061}
1062
1063void AliAnalysisEtMonteCarlo::FillHistograms()
537e541d 1064{ // let base class fill its histograms, and us fill the local ones
e2ee5727 1065 AliAnalysisEt::FillHistograms();
1066 //std::cout << fEtNonRemovedPiPlus << " " << fCentClass << std::endl;
1067
1068 fHistEtNonRemovedProtons->Fill(fEtNonRemovedProtons, fCentClass);
1069 fHistEtNonRemovedAntiProtons->Fill(fEtNonRemovedAntiProtons, fCentClass);
1070 fHistEtNonRemovedKaonPlus->Fill(fEtNonRemovedKaonPlus, fCentClass);
1071 fHistEtNonRemovedKaonMinus->Fill(fEtNonRemovedKaonMinus, fCentClass);
1072 fHistEtNonRemovedK0s->Fill(fEtNonRemovedK0s, fCentClass);
1073 fHistEtNonRemovedLambdas->Fill(fEtNonRemovedLambdas, fCentClass);
1074 fHistEtNonRemovedPiPlus->Fill(fEtNonRemovedPiPlus, fCentClass);
1075 fHistEtNonRemovedPiMinus->Fill(fEtNonRemovedPiMinus, fCentClass);
1076 fHistEtNonRemovedElectrons->Fill(fEtNonRemovedElectrons, fCentClass);
1077 fHistEtNonRemovedPositrons->Fill(fEtNonRemovedPositrons, fCentClass);
1078 fHistEtNonRemovedMuPlus->Fill(fEtNonRemovedMuPlus, fCentClass);
1079 fHistEtNonRemovedMuMinus->Fill(fEtNonRemovedMuMinus, fCentClass);
1080 fHistEtNonRemovedNeutrons->Fill(fEtNonRemovedNeutrons, fCentClass);
1081 fHistEtNonRemovedAntiNeutrons->Fill(fEtNonRemovedAntiNeutrons, fCentClass);
1082 fHistEtNonRemovedGammas->Fill(fEtNonRemovedGammas, fCentClass);
1083 fHistEtNonRemovedGammasFromPi0->Fill(fEtNonRemovedGammasFromPi0, fCentClass);
1084
1085 fHistEtRemovedGammas->Fill(fEtRemovedGammas, fCentClass);
1086 fHistEtRemovedNeutrons->Fill(fEtRemovedNeutrons, fCentClass);
1087 fHistEtRemovedAntiNeutrons->Fill(fEtRemovedAntiNeutrons, fCentClass);
1088
1089 fHistMultNonRemovedProtons->Fill(fMultNonRemovedProtons, fCentClass);
1090 fHistMultNonRemovedAntiProtons->Fill(fMultNonRemovedAntiProtons, fCentClass);
1091 fHistMultNonRemovedKaonPlus->Fill(fMultNonRemovedKaonPlus, fCentClass);
1092 fHistMultNonRemovedKaonMinus->Fill(fMultNonRemovedKaonMinus, fCentClass);
1093 fHistMultNonRemovedK0s->Fill(fMultNonRemovedK0s, fCentClass);
1094 fHistMultNonRemovedLambdas->Fill(fMultNonRemovedLambdas, fCentClass);
1095 fHistMultNonRemovedPiPlus->Fill(fMultNonRemovedPiPlus, fCentClass);
1096 fHistMultNonRemovedPiMinus->Fill(fMultNonRemovedPiMinus, fCentClass);
1097 fHistMultNonRemovedElectrons->Fill(fMultNonRemovedElectrons, fCentClass);
1098 fHistMultNonRemovedPositrons->Fill(fMultNonRemovedPositrons, fCentClass);
1099 fHistMultNonRemovedMuPlus->Fill(fMultNonRemovedMuPlus, fCentClass);
1100 fHistMultNonRemovedMuMinus->Fill(fMultNonRemovedMuMinus, fCentClass);
1101 fHistMultNonRemovedNeutrons->Fill(fMultNonRemovedNeutrons, fCentClass);
1102 fHistMultNonRemovedAntiNeutrons->Fill(fMultNonRemovedAntiNeutrons, fCentClass);
1103 fHistMultNonRemovedGammas->Fill(fMultNonRemovedGammas, fCentClass);
1104
1105 fHistMultRemovedGammas->Fill(fMultRemovedGammas, fCentClass);
1106 fHistMultRemovedNeutrons->Fill(fMultRemovedNeutrons, fCentClass);
1107 fHistMultRemovedAntiNeutrons->Fill(fMultRemovedAntiNeutrons, fCentClass);
1108
1109 fHistTrackMultvsNonRemovedCharged->Fill(fTrackMultInAcc,
1110 fMultNonRemovedAntiProtons+fMultNonRemovedElectrons+fMultNonRemovedKaonMinus+fMultNonRemovedKaonPlus
1111 +fMultNonRemovedMuMinus+fMultNonRemovedMuPlus+fMultNonRemovedPiMinus+fMultNonRemovedPiPlus+fMultNonRemovedPositrons
1112 +fMultNonRemovedProtons);
1113
1114 fHistTrackMultvsNonRemovedNeutral->Fill(fTrackMultInAcc,
1115 fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedLambdas);
1116
1117 fHistTrackMultvsRemovedGamma->Fill(fTrackMultInAcc,
1118 fMultRemovedGammas);
1119
1120 fHistClusterMultvsNonRemovedCharged->Fill(fNeutralMultiplicity,
1121 fMultNonRemovedAntiProtons+fMultNonRemovedElectrons+fMultNonRemovedKaonMinus
1122 +fMultNonRemovedKaonPlus+fMultNonRemovedMuMinus+fMultNonRemovedMuPlus
1123 +fMultNonRemovedPiMinus+fMultNonRemovedPiPlus+fMultNonRemovedPositrons+fMultNonRemovedProtons);
1124
1125 fHistClusterMultvsNonRemovedNeutral->Fill(fTrackMultInAcc,
1126 fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedLambdas);
1127
1128 fHistClusterMultvsRemovedGamma->Fill(fTrackMultInAcc,
1129 fMultRemovedGammas);
1130
13b0d3c1 1131}
1132
e2ee5727 1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147