056acd5baa7a45fdcd04a7a5aa4bc31323984d61
[u/mrichter/AliRoot.git] / PWGLF / totEt / AliAnalysisEtMonteCarlo.cxx
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)
8 //_________________________________________________________________________
9
10 #include "AliAnalysisEtMonteCarlo.h"
11 #include "AliAnalysisEtCuts.h"
12 #include "AliAnalysisEtSelector.h"
13 #include "AliAnalysisEtSelector.h"
14 #include "AliESDtrack.h"
15 #include "AliStack.h"
16 #include "AliVEvent.h"
17 #include "AliMCEvent.h"
18 #include "AliESDEvent.h"
19 #include "TH2F.h"
20 #include "TH3F.h"
21 #include "TParticle.h"
22 #include "AliGenHijingEventHeader.h"
23 #include "AliGenPythiaEventHeader.h"
24 #include "TList.h"
25 #include "AliESDCaloCluster.h"
26 #include "AliLog.h"
27 #include <iostream>
28 #include <AliCentrality.h>
29 #include "AliPHOSGeoUtils.h"
30 #include "AliPHOSGeometry.h"
31 #include "TFile.h"
32 using namespace std;
33
34 ClassImp(AliAnalysisEtMonteCarlo);
35
36
37 // ctor
38 AliAnalysisEtMonteCarlo::AliAnalysisEtMonteCarlo():AliAnalysisEt()
39     ,fImpactParameter(0)
40     ,fNcoll(0)
41     ,fNpart(0)
42     ,fPrimaryTree(0)
43     ,fTotEtWithSecondaryRemoved(0)
44     ,fTotEtSecondaryFromEmEtPrimary(0)
45     ,fTotEtSecondary(0)
46     ,fPrimaryCode(0)
47     ,fPrimaryCharge(0)
48     ,fPrimaryE(0)
49     ,fPrimaryEt(0)
50     ,fPrimaryPx(0)
51     ,fPrimaryPy(0)
52     ,fPrimaryPz(0)
53     ,fPrimaryVx(0)
54     ,fPrimaryVy(0)
55     ,fPrimaryVz(0)
56     ,fPrimaryAccepted(0)
57     ,fPrimaryMatched(0)
58     ,fDepositedCode(0)
59     ,fDepositedE(0)
60     ,fDepositedEt(0)
61     ,fDepositedCharge(0)
62     ,fDepositedVx(0)
63     ,fDepositedVy(0)
64     ,fDepositedVz(0)
65     ,fSecondary(kFALSE)
66     ,fReconstructedE(0)
67     ,fReconstructedEt(0)
68     ,fTotPx(0)
69     ,fTotPy(0)
70     ,fTotPz(0)
71     ,fClusterMult(0)
72     ,fHistDecayVertexNonRemovedCharged(0)
73     ,fHistDecayVertexRemovedCharged(0)
74     ,fHistDecayVertexNonRemovedNeutral(0)
75     ,fHistDecayVertexRemovedNeutral(0)
76
77     ,fHistRemovedOrNot(0)
78     ,fHistEtNonRemovedProtons(0)
79     ,fHistEtNonRemovedAntiProtons(0)
80     ,fHistEtNonRemovedPiPlus(0)
81     ,fHistEtNonRemovedPiMinus(0)
82     ,fHistEtNonRemovedKaonPlus(0)
83     ,fHistEtNonRemovedKaonMinus(0)
84     ,fHistEtNonRemovedK0s(0)
85     ,fHistEtNonRemovedK0L(0)
86     ,fHistEtNonRemovedLambdas(0)
87     ,fHistEtNonRemovedElectrons(0)
88     ,fHistEtNonRemovedPositrons(0)
89     ,fHistEtNonRemovedMuPlus(0)
90     ,fHistEtNonRemovedMuMinus(0)
91     ,fHistEtNonRemovedNeutrons(0)
92     ,fHistEtNonRemovedAntiNeutrons(0)
93     ,fHistEtNonRemovedGammas(0)
94     ,fHistEtNonRemovedGammasFromPi0(0)
95     ,fHistEtRemovedGammas(0)
96     ,fHistEtRemovedNeutrons(0)
97     ,fHistEtRemovedAntiNeutrons(0)
98     ,fHistEtRemovedCharged(0)
99     ,fHistEtRemovedNeutrals(0)
100     ,fHistEtNonRemovedCharged(0)
101     ,fHistEtNonRemovedNeutrals(0)
102     ,fHistMultNonRemovedProtons(0)
103     ,fHistMultNonRemovedAntiProtons(0)
104     ,fHistMultNonRemovedPiPlus(0)
105     ,fHistMultNonRemovedPiMinus(0)
106     ,fHistMultNonRemovedKaonPlus(0)
107     ,fHistMultNonRemovedKaonMinus(0)
108     ,fHistMultNonRemovedK0s(0)
109     ,fHistMultNonRemovedK0L(0)
110     ,fHistMultNonRemovedLambdas(0)
111     ,fHistMultNonRemovedElectrons(0)
112     ,fHistMultNonRemovedPositrons(0)
113     ,fHistMultNonRemovedMuPlus(0)
114     ,fHistMultNonRemovedMuMinus(0)
115     ,fHistMultNonRemovedNeutrons(0)
116     ,fHistMultNonRemovedAntiNeutrons(0)
117     ,fHistMultNonRemovedGammas(0)
118     ,fHistMultRemovedGammas(0)
119     ,fHistMultRemovedNeutrons(0)
120     ,fHistMultRemovedAntiNeutrons(0)
121     ,fHistMultRemovedCharged(0)
122     ,fHistMultRemovedNeutrals(0)
123     ,fHistMultNonRemovedCharged(0)
124     ,fHistMultNonRemovedNeutrals(0)
125     ,fHistTrackMultvsNonRemovedCharged(0)
126     ,fHistTrackMultvsNonRemovedNeutral(0)
127     ,fHistTrackMultvsRemovedGamma(0)
128     ,fHistClusterMultvsNonRemovedCharged(0)
129     ,fHistClusterMultvsNonRemovedNeutral(0)
130     ,fHistClusterMultvsRemovedGamma(0)
131     ,fHistMultvsNonRemovedChargedE(0)
132     ,fHistMultvsNonRemovedNeutralE(0)
133     ,fHistMultvsRemovedGammaE(0)
134                                                   ,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)
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)
231 //     ,fPiPlusMult(0)
232 //     ,fPiMinusMult(0)
233     ,fPiZeroMult(0)
234     ,fPiPlusMultAcc(0)
235     ,fPiMinusMultAcc(0)
236     ,fPiZeroMultAcc(0)
237     ,fNeutralRemoved(0)
238     ,fChargedRemoved(0)
239     ,fChargedNotRemoved(0)
240     ,fNeutralNotRemoved(0)
241     ,fGammaRemoved(0)
242     ,fSecondaryNotRemoved(0)
243     ,fEnergyNeutralRemoved(0)
244     ,fEnergyChargedRemoved(0)
245     ,fEnergyChargedNotRemoved(0)
246     ,fEnergyNeutralNotRemoved(0)
247     ,fEnergyGammaRemoved(0)
248     ,fNClusters(0)
249     ,fTotNeutralEtAfterMinEnergyCut(0)
250     ,fHistGammasFound(0)
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)
267 {
268 }
269
270 // dtor
271 AliAnalysisEtMonteCarlo::~AliAnalysisEtMonteCarlo()
272 {   //Destructor
273
274   if(fPrimaryTree){
275     fPrimaryTree->Clear();
276     delete fPrimaryTree;
277   }
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
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
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
374     delete fHistGammasFound; // enter comment here
375     delete fHistGammasGenerated; // enter comment here
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;
391 }
392
393 Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
394 {   // analyse MC event
395     ResetEventValues();
396
397     
398     fPiPlusMult = 0;
399     fPiMinusMult = 0;
400     fPiZeroMult = 0;
401     if (fCentrality)
402     {
403         fCentClass = fCentrality->GetCentralityClass10(fCentralityMethod);
404
405     }
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();
431     }
432
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++)
440     {
441
442         TParticle *part = stack->Particle(iPart);
443         
444         
445
446         if (!part)
447         {
448             Printf("ERROR: Could not get particle %d", iPart);
449             continue;
450         }
451         TParticlePDG *pdg = part->GetPDG(0);
452
453         if (!pdg)
454         {
455             Printf("ERROR: Could not get particle PDG %d", iPart);
456             continue;
457         }
458
459         Double_t particleMassPart = 0; //The mass part in the Et calculation for this particle
460         Int_t code = pdg->PdgCode();
461         
462         if(stack->IsPhysicalPrimary(iPart))
463         {
464           fTotPx += part->Px();
465           fTotPy += part->Py();
466           fTotPz += part->Pz();
467         }
468         
469         // Check for reasonable (for now neutral and singly charged) charge on the particle
470         if (fSelector->IsNeutralMcParticle(iPart,*stack,*pdg))
471         {
472
473             fMultiplicity++;
474 //             PrintFamilyTree(iPart, stack);
475 //
476             if (TMath::Abs(part->Eta()) < fCuts->GetCommonEtaCut())
477             {
478                 //Printf("Particle with eta: %f, pid: %d", part->Eta(), code);
479                 // calculate E_T
480                 if (
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
487                 )
488                 {
489                     if (code > 0) {
490                         particleMassPart = - protonMass;
491                     }
492                     if (code < 0) {
493                         particleMassPart = protonMass;
494                     }
495                 }
496                 Double_t et = part->Energy() * TMath::Sin(part->Theta()) + particleMassPart;
497
498
499                 // Fill up total E_T counters for each particle species
500                 if (code == fgProtonCode || code == fgAntiProtonCode)
501                 {
502                 }
503                 if (code == fgPiPlusCode || code == fgPiMinusCode)
504                 {
505                     if (code == fgPiPlusCode)
506                     {
507                     }
508                     else
509                     {
510                     }
511                 }
512                 if (code == fgGammaCode)
513                 {
514                 }
515                 if (code == fgKPlusCode)
516                 {
517                 }
518                 if(code == fgKMinusCode)
519                 {
520                 }
521                 if (code == fgMuPlusCode || code == fgMuMinusCode)
522                 {
523                 }
524                 if (code == fgEPlusCode || code == fgEMinusCode)
525                 {
526                 }
527                 // some neutrals also
528                 if (code == fgNeutronCode)
529                 {
530                 }
531                 if (code == fgAntiNeutronCode)
532                 {
533                 }
534                 if (code == fgGammaCode)
535                 {
536                 }
537
538                 // Neutral particles
539                 //if (TMath::Abs(pdg->Charge() - fCuts->GetMonteCarloNeutralParticle()) <1e-3 )
540
541                 if(code == fgGammaCode || code == fgPi0Code || code == fgEtaCode)
542                 {
543                 //  PrintFamilyTree(iPart,stack);
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;
546
547                     // inside EMCal acceptance
548
549                     //if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiMaxCutAcc && part->Phi() > fPhiMinCutAcc)
550
551                     if(fSelector->CutGeometricalAcceptance(*part) )
552                     {
553                         fNeutralMultiplicity++;
554                         fTotNeutralEt += et;
555                         if(fSelector->PassMinEnergyCut(*part))
556                         {
557                             fTotNeutralEtAfterMinEnergyCut += et;
558                         }
559                         if (part->Energy() > 0.05) partCount++;
560                     }
561                 }
562                 //Charged particles
563                 else if (TMath::Abs( pdg->Charge() - fCuts->GetMonteCarloNeutralParticle())>1e-3 )
564                 {
565
566                     // inside EMCal acceptance
567                     if (fSelector->CutGeometricalAcceptance(*part))
568                     {
569
570                         fChargedMultiplicity++;
571
572                         fTotChargedEt += et;
573
574                         if (code == fgProtonCode || code == fgAntiProtonCode)
575                         {
576                         }
577                         if (code == fgPiPlusCode || code == fgPiMinusCode)
578                         {
579                         }
580                         if (code == fgKPlusCode || code == fgKMinusCode)
581                         {
582                         }
583                         if (code == fgMuPlusCode || code == fgMuMinusCode)
584                         {
585                         }
586
587                         if (code == fgEPlusCode || code == fgEMinusCode)
588                         {
589                             fTotNeutralEt += et; // calling electrons neutral
590                             fTotChargedEt -= et;
591                         }
592                     } // inside EMCal acceptance
593
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         }
603     }
604    // std::cout << "Total: p_x = " << fTotPx << ", p_y = " << fTotPy << ", p_z = " << fTotPz << std::endl;
605     fTotEt = fTotChargedEt + fTotNeutralEt;
606     //fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;//
607     //std::cout << "Event done! # of particles: " << partCount << std::endl;
608     return 0;
609 }
610 //Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliMCEvent* mcEvent,AliESDEvent* realEvent)
611 Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
612 {   // analyse MC and real event info
613     //if(!mcEvent || !realEvent){
614     if (!ev || !ev2) {
615         AliFatal("ERROR: Event does not exist");
616         return 0;
617     }
618     AliAnalysisEt::AnalyseEvent(ev);
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");
623        
624     }
625
626     std::vector<Int_t> foundGammas;
627     
628     fSelector->SetEvent(realEvent);
629
630     AnalyseEvent(ev);
631
632     AliStack *stack = mcEvent->Stack();
633
634     // get all detector clusters
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
644 //Note that this only returns clusters for the selected detector.  fSelector actually calls the right GetClusters... for the detector
645 //It does not apply any cuts on these clusters
646     TRefArray *caloClusters = fSelector->GetClusters();
647
648     Int_t nCluster = caloClusters->GetEntries();
649     fClusterMult = nCluster;
650     fNClusters = 0;
651     // loop the clusters
652     for (int iCluster = 0; iCluster < nCluster; iCluster++ )
653     {
654         Int_t cf = 0;
655         AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
656         //Float_t caloE = caloCluster->E()
657         fNClusters++;
658         const UInt_t iPart = (UInt_t)TMath::Abs(caloCluster->GetLabel());
659         TParticle *part  =  stack->Particle(iPart);
660
661         if (!part)
662         {
663             Printf("No MC particle %d", iCluster);
664             continue;
665         }
666
667         int primIdx = iPart;
668         if (!stack->IsPhysicalPrimary(iPart)) // check whether particle is primary. we keep secondary electron and gamma for testing.
669         {
670             primIdx = GetPrimMother(iPart, stack);
671             if(primIdx != stack->GetPrimary(iPart))
672             {
673               //std::cout << primIdx << " = " << stack->GetPrimary(iPart) << std::endl;
674               //PrintFamilyTree(iPart, stack);
675             }
676             //if it is from a K0S
677             if(primIdx < 0)
678             {
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;
684             }
685
686         } // end of primary particle check
687         
688         //const int primCode = stack->Particle(primIdx)->GetPdgCode();
689         TParticlePDG *pdg = part->GetPDG();
690         if (!pdg)
691         {
692             Printf("ERROR: Could not get particle PDG %d", iPart);
693             continue;
694         }
695
696         //Int_t code = pdg->PdgCode();
697 //      if(primCode == fgGammaCode) 
698 //      {
699           
700         for(UInt_t i = 0; i < caloCluster->GetNLabels(); i++)
701         {
702           Int_t pIdx = caloCluster->GetLabelAt(i);
703           //TParticle *p = stack->Particle(pIdx);
704           
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))))
711           {
712 //          std::cout << "Gamma primary: " << primIdx << std::endl;
713 //          foundGammas.push_back(primIdx); 
714             foundGammas.push_back(pIdx); 
715           }
716         }
717         fCutFlow->Fill(cf++);
718         if(!fSelector->PassDistanceToBadChannelCut(*caloCluster)) continue;
719         Double_t clEt = CorrectForReconstructionEfficiency(*caloCluster);
720 //      if(code == fgK0SCode) std::cout << "K0 energy: " << caloCluster->E() << std::endl;
721         if(!fSelector->PassMinEnergyCut(*caloCluster)) continue;
722
723         
724         fCutFlow->Fill(cf++);
725         Float_t pos[3];
726         //PrintFamilyTree(
727         caloCluster->GetPosition(pos);
728         TVector3 cp(pos);
729
730         TParticle *primPart = stack->Particle(primIdx);
731         fPrimaryCode = primPart->GetPdgCode();
732         fPrimaryCharge = (Int_t) primPart->GetPDG()->Charge();
733
734         fPrimaryE = primPart->Energy();
735         fPrimaryEt = primPart->Energy()*TMath::Sin(primPart->Theta());
736         fPrimaryPx = primPart->Px();
737         fPrimaryPy = primPart->Py();
738         fPrimaryPz = primPart->Pz();
739         //cout<<"I have a cluster and it's good energy "<<caloCluster->E()<<" simulated "<<fPrimaryE<<endl;
740         fPrimaryVx = primPart->Vx();
741         fPrimaryVy = primPart->Vy();
742         fPrimaryVz = primPart->Vz();
743
744         fPrimaryAccepted = false;
745         fPrimaryMatched = false;
746
747         fDepositedCode = part->GetPdgCode();
748         fDepositedE = part->Energy();
749         fDepositedEt = part->Energy()*TMath::Sin(part->Theta());
750         fDepositedCharge = (Int_t) part->GetPDG()->Charge();
751
752         fDepositedVx = part->Vx();
753         fDepositedVy = part->Vy();
754         fDepositedVz = part->Vz();
755
756         //fSecondary = fSelector->FromSecondaryInteraction(*primPart, *stack);
757         fSecondary =fSelector->FromSecondaryInteraction(*part, *stack);
758 //      if(fSecondary) 
759 //      {
760 //        //std::cout << "Have secondary!" << std::endl;
761 //        //PrintFamilyTree(iPart, stack);
762 //      }
763         fReconstructedE = caloCluster->E();
764         fReconstructedEt = caloCluster->E()*TMath::Sin(cp.Theta());
765         
766         pdg = primPart->GetPDG(0);
767         //Int_t code = primPart->GetPdgCode();
768
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
849           }
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         }
886         fPrimaryTree->Fill();
887     } // end of loop over clusters     
888
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());
914           if(std::binary_search(foundGammas.begin(),foundGammas.end(),iPart))
915           {
916             fHistGammasFound->Fill(part->Energy());
917           }
918         }
919         
920     }
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);
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               }
959               for(int k=0;k<4;k++){
960                 //if( TMath::Abs(stack->Particle(gammaDaughterIDs[k])->Eta()) <= fCuts->GetGeometryEmcalEtaAccCut() ){//only add the energy if it's within the detector acceptance
961                 if(gammaDaughterIDs[k]==-1) continue;
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                 }
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               }
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
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++){
993                   if( myPart==((UInt_t) gammaDaughterIDs[j]) ){
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;
998                     for(int l=0;l<nEtCuts;l++){//loop over cut values
999                       //cout<<", cut = "<<etCuts[l];
1000                       if(clEt>=etCuts[l]){
1001                         //cout<<", "<<clEt<<">="<<etCuts[l];
1002                         totalClusterEts[l] += clEt;//if cluster et is above the cut off energy add it
1003                       }
1004                     }
1005                     //cout<<endl;
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     
1037     FillHistograms();
1038     return 0;
1039 }
1040
1041 void AliAnalysisEtMonteCarlo::Init()
1042 {   // init
1043     AliAnalysisEt::Init();
1044 }
1045
1046 void AliAnalysisEtMonteCarlo::ResetEventValues()
1047 {   // reset event values
1048     AliAnalysisEt::ResetEventValues();
1049
1050     fTotEtSecondary = 0;
1051     fTotEtSecondaryFromEmEtPrimary = 0;
1052     fTotEtWithSecondaryRemoved = 0;
1053
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;
1065     fEtNonRemovedK0S = 0;
1066     fEtNonRemovedK0L = 0;
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
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;
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;
1104     fMultNonRemovedK0L = 0;
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
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
1128     fMultRemovedGammas = 0;
1129     fMultRemovedNeutrons = 0;
1130     fMultRemovedAntiNeutrons = 0;
1131
1132     fEnergyChargedNotRemoved = 0;
1133     fEnergyChargedRemoved = 0;
1134     fEnergyNeutralNotRemoved = 0;
1135     fEnergyNeutralRemoved = 0;
1136
1137     fChargedNotRemoved = 0;
1138     fChargedRemoved = 0;
1139     fNeutralNotRemoved = 0;
1140     fNeutralRemoved = 0;
1141     fGammaRemoved = 0;
1142     fSecondaryNotRemoved = 0;
1143
1144     fTrackMultInAcc = 0;
1145
1146     fTotNeutralEtAfterMinEnergyCut = 0;
1147   
1148     fSecondaryNotRemoved = 0;
1149     
1150     fTotPx = 0;
1151     fTotPy = 0;
1152     fTotPz = 0;
1153     
1154     
1155 }
1156
1157 void AliAnalysisEtMonteCarlo::CreateHistograms()
1158 {   // histogram related additions
1159     AliAnalysisEt::CreateHistograms();
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");
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
1178     }
1179
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);
1194     fHistEtNonRemovedK0L = new TH2F("fHistEtNonRemovedK0L", "fHistEtNonRemovedK0L", 1500, 0, 30, 10, -0.5, 9.5);
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);
1209
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);
1212
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);
1215
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);
1223     fHistMultNonRemovedK0L = new TH2F("fHistMultNonRemovedK0L", "fHistMultNonRemovedK0L", 100, -0.5, 99.5, 10, -0.5, 9.5);
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
1232     fHistMultNonRemovedGammas = new  TH2F("fHistMultNonRemovedGammas", "fHistMultNonRemovedGammas", 100, -0.5, 99.5, 100, -0.5, 99.5);
1233
1234     fHistMultRemovedGammas = new  TH2F("fHistMultRemovedGammas", "fHistMultRemovedGammas", 100, -0.5, 99.5, 100, -0.5, 99.5);
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);
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);*/
1243
1244
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);
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
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
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
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");
1317         fPrimaryTree->Branch("fPrimaryMatched", &fPrimaryMatched, "fPrimaryMatched/B");
1318
1319
1320         fPrimaryTree->Branch("fDepositedCode", &fDepositedCode, "fDepositedCode/I");
1321         fPrimaryTree->Branch("fDepositedCharge", &fDepositedCharge, "fDepositedCharge/I");
1322         fPrimaryTree->Branch("fDepositedE", &fDepositedE, "fDepositedE/D");
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
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         
1338     }
1339
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);
1357 }
1358
1359 void AliAnalysisEtMonteCarlo::FillOutputList(TList *list)
1360 {   //fill the output list
1361     AliAnalysisEt::FillOutputList(list);
1362
1363     if(fCuts->GetHistMakeTree())
1364     {
1365         list->Add(fPrimaryTree);
1366     }
1367
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);
1377     list->Add(fHistEtNonRemovedK0L);
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
1392     list->Add(fHistEtRemovedCharged);
1393     list->Add(fHistEtRemovedNeutrals);
1394
1395     list->Add(fHistEtNonRemovedCharged);
1396     list->Add(fHistEtNonRemovedNeutrals);
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);
1405     list->Add(fHistMultNonRemovedK0L);
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
1419     list->Add(fHistMultRemovedCharged);
1420     list->Add(fHistMultRemovedNeutrals);
1421
1422     list->Add(fHistMultNonRemovedCharged);
1423     list->Add(fHistMultNonRemovedNeutrals);
1424
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
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
1464     list->Add(fHistPiPlusMult);
1465     list->Add(fHistPiMinusMult);
1466     list->Add(fHistPiZeroMult);
1467     list->Add(fHistPiPlusMultAcc);
1468     list->Add(fHistPiMinusMultAcc);
1469     list->Add(fHistPiZeroMultAcc);
1470     
1471     list->Add(fHistGammasFound);
1472     list->Add(fHistGammasGenerated);
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);
1488
1489 }
1490
1491
1492 bool AliAnalysisEtMonteCarlo::TrackHitsCalorimeter(TParticle* part, Double_t magField)
1493 {
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
1502     bool status = prop && fSelector->CutGeometricalAcceptance(*esdTrack);
1503     delete esdTrack;
1504
1505     return status;
1506 }
1507
1508 void AliAnalysisEtMonteCarlo::FillHistograms()
1509 {   // let base class fill its histograms, and us fill the local ones
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);
1517     fHistEtNonRemovedK0s->Fill(fEtNonRemovedK0S, fCentClass);
1518     fHistEtNonRemovedK0L->Fill(fEtNonRemovedK0L, fCentClass);
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
1531     fHistEtRemovedGammas->Fill(fEtRemovedGammas, fNClusters);
1532     fHistEtRemovedNeutrons->Fill(fEtRemovedNeutrons, fCentClass);
1533     fHistEtRemovedAntiNeutrons->Fill(fEtRemovedAntiNeutrons, fCentClass);
1534
1535 //     fHistEtRemovedCharged->Fill(fEtRemovedAntiProtons+fEtRemovedElectrons+fEtRemovedKaonMinus+fEtRemovedKaonPlus
1536 //                                             +fEtRemovedMuMinus+fEtRemovedMuPlus+fEtRemovedPiMinus+fEtRemovedPiPlus+fEtRemovedPositrons
1537 //                                             +fEtRemovedProtons.
1538 //                              fCentClass);
1539 //     fHistEtRemovedNeutrals->Fill(fEtRemovedNeutrons+fEtRemovedAntiNeutrons, fCentClass);
1540 //
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);
1551
1552     fHistMultRemovedCharged->Fill(fChargedRemoved, fNClusters);
1553     fHistMultRemovedNeutrals->Fill(fNeutralRemoved, fNClusters);
1554     fHistMultNonRemovedCharged->Fill(fChargedNotRemoved, fNClusters);
1555     fHistMultNonRemovedNeutrals->Fill(fNeutralNotRemoved, fNClusters);
1556
1557
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);
1563     fHistMultNonRemovedK0L->Fill(fMultNonRemovedK0L, fCentClass);
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,
1585                                             fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedK0L+fMultNonRemovedLambdas+fK0sMult);
1586
1587     fHistTrackMultvsRemovedGamma->Fill(fTrackMultInAcc,
1588                                        fMultRemovedGammas);
1589
1590     fHistClusterMultvsNonRemovedCharged->Fill(fNClusters,
1591             fMultNonRemovedAntiProtons+fMultNonRemovedElectrons+fMultNonRemovedKaonMinus
1592             +fMultNonRemovedKaonPlus+fMultNonRemovedMuMinus+fMultNonRemovedMuPlus
1593             +fMultNonRemovedPiMinus+fMultNonRemovedPiPlus+fMultNonRemovedPositrons+fMultNonRemovedProtons);
1594
1595     fHistClusterMultvsNonRemovedNeutral->Fill(fNClusters,
1596             fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedK0L+fMultNonRemovedLambdas+fK0sMult);
1597
1598     fHistClusterMultvsRemovedGamma->Fill(fNClusters,
1599                                          fMultRemovedGammas);
1600
1601 }
1602
1603
1604
1605
1606 Int_t AliAnalysisEtMonteCarlo::PrintFamilyTree(Int_t partIdx, AliStack* stack)
1607 { // print family tree
1608     TParticle *part = stack->Particle(partIdx);
1609 //     if(part->GetPdgCode() == fgK0SCode)
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);
1617 }
1618
1619 Int_t AliAnalysisEtMonteCarlo::PrintMothers(Int_t partIdx, AliStack* stack, Int_t gen)
1620 { // print mothers
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);
1635 //     if(mother->GetPdgCode() == fgK0SCode)
1636     {
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;
1640         std::cout << tabs << "PID: " << mother->GetPdgCode() << "/" << mother->GetName() << std::endl;
1641         std::cout << tabs << "Energy: " << mother->Energy() << std::endl;
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         }
1648         std::cout << tabs << "Vertex: " << mother->Vx() << ", " << mother->Vy() << ", " << mother->Vz() << std::endl;
1649     }
1650     if(mother->GetPdgCode() == fgK0SCode)
1651     {
1652 //      std::cout << "K0S!!!!!!!!!!!!!11111!!!!!" << std::endl;
1653     }
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;
1658
1659     delete [] tabs;
1660     return PrintMothers(mothIdx, stack, gen+1) + 1;
1661 }
1662
1663 Int_t AliAnalysisEtMonteCarlo::GetPrimMother(Int_t partIdx, AliStack *stack)
1664 { // get primary mother
1665     if(partIdx >= 0)
1666     {
1667         //return stack->GetPrimary(partIdx);
1668         
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         {
1674             if(stack->IsPhysicalPrimary(mothIdx)) return mothIdx;
1675             else return GetPrimMother(mothIdx, stack);
1676         }
1677         else
1678         {
1679             return -1;
1680         }
1681     }
1682     return -1;
1683 }
1684
1685 Int_t AliAnalysisEtMonteCarlo::GetK0InFamily(Int_t partIdx, AliStack* stack)
1686 { // get K0 in family
1687     if(partIdx >= 0)
1688     {
1689         if(stack->Particle(partIdx)->GetPdgCode() == fgK0SCode) return partIdx;
1690         Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
1691         if(mothIdx < 0) return -1;
1692         TParticle *mother = stack->Particle(mothIdx);
1693         if(mother)
1694         {
1695           if(mother->GetPdgCode() == fgK0SCode)
1696             {
1697                 return mothIdx;
1698             }
1699             return GetK0InFamily(mothIdx, stack);
1700         }
1701         else
1702         {
1703             return -1;
1704         }
1705     }
1706     return -1;
1707 }
1708