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