1 //_________________________________________________________________________
2 // Utility Class for transverse energy studies
3 // Base class for MC analysis
7 //*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL)
8 //_________________________________________________________________________
10 #include "AliAnalysisEtMonteCarlo.h"
11 #include "AliAnalysisEtCuts.h"
12 #include "AliAnalysisEtSelector.h"
13 #include "AliAnalysisEtSelector.h"
14 #include "AliESDtrack.h"
16 #include "AliVEvent.h"
17 #include "AliMCEvent.h"
18 #include "AliESDEvent.h"
21 #include "TParticle.h"
22 #include "AliGenHijingEventHeader.h"
23 #include "AliGenPythiaEventHeader.h"
25 #include "AliESDCaloCluster.h"
28 #include <AliCentrality.h>
29 #include "AliPHOSGeoUtils.h"
30 #include "AliPHOSGeometry.h"
34 ClassImp(AliAnalysisEtMonteCarlo);
38 AliAnalysisEtMonteCarlo::AliAnalysisEtMonteCarlo():AliAnalysisEt()
43 ,fTotEtWithSecondaryRemoved(0)
44 ,fTotEtSecondaryFromEmEtPrimary(0)
63 ,fHistDecayVertexNonRemovedCharged(0)
64 ,fHistDecayVertexRemovedCharged(0)
65 ,fHistDecayVertexNonRemovedNeutral(0)
66 ,fHistDecayVertexRemovedNeutral(0)
69 ,fHistEtNonRemovedProtons(0)
70 ,fHistEtNonRemovedAntiProtons(0)
71 ,fHistEtNonRemovedPiPlus(0)
72 ,fHistEtNonRemovedPiMinus(0)
73 ,fHistEtNonRemovedKaonPlus(0)
74 ,fHistEtNonRemovedKaonMinus(0)
75 ,fHistEtNonRemovedK0s(0)
76 ,fHistEtNonRemovedK0L(0)
77 ,fHistEtNonRemovedLambdas(0)
78 ,fHistEtNonRemovedElectrons(0)
79 ,fHistEtNonRemovedPositrons(0)
80 ,fHistEtNonRemovedMuPlus(0)
81 ,fHistEtNonRemovedMuMinus(0)
82 ,fHistEtNonRemovedNeutrons(0)
83 ,fHistEtNonRemovedAntiNeutrons(0)
84 ,fHistEtNonRemovedGammas(0)
85 ,fHistEtNonRemovedGammasFromPi0(0)
86 ,fHistEtRemovedGammas(0)
87 ,fHistEtRemovedNeutrons(0)
88 ,fHistEtRemovedAntiNeutrons(0)
89 ,fHistEtRemovedCharged(0)
90 ,fHistEtRemovedNeutrals(0)
91 ,fHistEtNonRemovedCharged(0)
92 ,fHistEtNonRemovedNeutrals(0)
93 ,fHistMultNonRemovedProtons(0)
94 ,fHistMultNonRemovedAntiProtons(0)
95 ,fHistMultNonRemovedPiPlus(0)
96 ,fHistMultNonRemovedPiMinus(0)
97 ,fHistMultNonRemovedKaonPlus(0)
98 ,fHistMultNonRemovedKaonMinus(0)
99 ,fHistMultNonRemovedK0s(0)
100 ,fHistMultNonRemovedK0L(0)
101 ,fHistMultNonRemovedLambdas(0)
102 ,fHistMultNonRemovedElectrons(0)
103 ,fHistMultNonRemovedPositrons(0)
104 ,fHistMultNonRemovedMuPlus(0)
105 ,fHistMultNonRemovedMuMinus(0)
106 ,fHistMultNonRemovedNeutrons(0)
107 ,fHistMultNonRemovedAntiNeutrons(0)
108 ,fHistMultNonRemovedGammas(0)
109 ,fHistMultRemovedGammas(0)
110 ,fHistMultRemovedNeutrons(0)
111 ,fHistMultRemovedAntiNeutrons(0)
112 ,fHistMultRemovedCharged(0)
113 ,fHistMultRemovedNeutrals(0)
114 ,fHistMultNonRemovedCharged(0)
115 ,fHistMultNonRemovedNeutrals(0)
116 ,fHistTrackMultvsNonRemovedCharged(0)
117 ,fHistTrackMultvsNonRemovedNeutral(0)
118 ,fHistTrackMultvsRemovedGamma(0)
119 ,fHistClusterMultvsNonRemovedCharged(0)
120 ,fHistClusterMultvsNonRemovedNeutral(0)
121 ,fHistClusterMultvsRemovedGamma(0)
122 ,fHistMultvsNonRemovedChargedE(0)
123 ,fHistMultvsNonRemovedNeutralE(0)
124 ,fHistMultvsRemovedGammaE(0)
125 ,fEtNonRemovedProtons(0)
126 ,fEtNonRemovedAntiProtons(0)
127 ,fEtNonRemovedPiPlus(0)
128 ,fEtNonRemovedPiMinus(0)
129 ,fEtNonRemovedKaonPlus(0)
130 ,fEtNonRemovedKaonMinus(0)
133 ,fEtNonRemovedLambdas(0)
134 ,fEtNonRemovedElectrons(0)
135 ,fEtNonRemovedPositrons(0)
136 ,fEtNonRemovedMuMinus(0)
137 ,fEtNonRemovedMuPlus(0)
138 ,fEtNonRemovedGammas(0)
139 ,fEtNonRemovedGammasFromPi0(0)
140 ,fEtNonRemovedNeutrons(0)
141 ,fEtNonRemovedAntiNeutrons(0)
142 ,fEtRemovedProtons(0)
143 ,fEtRemovedAntiProtons(0)
145 ,fEtRemovedPiMinus(0)
146 ,fEtRemovedKaonPlus(0)
147 ,fEtRemovedKaonMinus(0)
150 ,fEtRemovedLambdas(0)
151 ,fEtRemovedElectrons(0)
152 ,fEtRemovedPositrons(0)
153 ,fEtRemovedMuMinus(0)
155 ,fEtRemovedGammasFromPi0(0)
157 ,fEtRemovedNeutrons(0)
158 ,fEtRemovedAntiNeutrons(0)
159 ,fMultNonRemovedProtons(0)
160 ,fMultNonRemovedAntiProtons(0)
161 ,fMultNonRemovedPiPlus(0)
162 ,fMultNonRemovedPiMinus(0)
163 ,fMultNonRemovedKaonPlus(0)
164 ,fMultNonRemovedKaonMinus(0)
165 ,fMultNonRemovedK0s(0)
166 ,fMultNonRemovedK0L(0)
167 ,fMultNonRemovedLambdas(0)
168 ,fMultNonRemovedElectrons(0)
169 ,fMultNonRemovedPositrons(0)
170 ,fMultNonRemovedMuMinus(0)
171 ,fMultNonRemovedMuPlus(0)
172 ,fMultNonRemovedGammas(0)
173 ,fMultNonRemovedNeutrons(0)
174 ,fMultNonRemovedAntiNeutrons(0)
175 ,fMultRemovedProtons(0)
176 ,fMultRemovedAntiProtons(0)
177 ,fMultRemovedPiPlus(0)
178 ,fMultRemovedPiMinus(0)
179 ,fMultRemovedKaonPlus(0)
180 ,fMultRemovedKaonMinus(0)
183 ,fMultRemovedLambdas(0)
184 ,fMultRemovedElectrons(0)
185 ,fMultRemovedPositrons(0)
186 ,fMultRemovedMuMinus(0)
187 ,fMultRemovedMuPlus(0)
188 ,fMultRemovedGammas(0)
189 ,fMultRemovedNeutrons(0)
190 ,fMultRemovedAntiNeutrons(0)
192 ,fHistDxDzNonRemovedCharged(0)
193 ,fHistDxDzRemovedCharged(0)
194 ,fHistDxDzNonRemovedNeutral(0)
195 ,fHistDxDzRemovedNeutral(0)
199 ,fHistPiPlusMultAcc(0)
200 ,fHistPiMinusMultAcc(0)
201 ,fHistPiZeroMultAcc(0)
210 ,fChargedNotRemoved(0)
211 ,fNeutralNotRemoved(0)
212 ,fEnergyNeutralRemoved(0)
213 ,fEnergyChargedRemoved(0)
214 ,fEnergyChargedNotRemoved(0)
215 ,fEnergyNeutralNotRemoved(0)
217 ,fTotNeutralEtAfterMinEnergyCut(0)
223 AliAnalysisEtMonteCarlo::~AliAnalysisEtMonteCarlo()
227 fPrimaryTree->Clear();
230 delete fHistDecayVertexNonRemovedCharged; // Decay vertex for non-removed charged particles
231 delete fHistDecayVertexRemovedCharged; // Decay vertex for non-removed charged particles
232 delete fHistDecayVertexNonRemovedNeutral; // Decay vertex for non-removed charged particles
233 delete fHistDecayVertexRemovedNeutral; // Decay vertex for non-removed charged particles
235 delete fHistRemovedOrNot; // If charged/neutral particles were removed or not
237 delete fHistEtNonRemovedProtons; // enter comment here
238 delete fHistEtNonRemovedAntiProtons; // enter comment here
239 delete fHistEtNonRemovedPiPlus; // enter comment here
240 delete fHistEtNonRemovedPiMinus; // enter comment here
241 delete fHistEtNonRemovedKaonPlus; // enter comment here
242 delete fHistEtNonRemovedKaonMinus; // enter comment here
243 delete fHistEtNonRemovedK0s; // enter comment here
244 delete fHistEtNonRemovedK0L; // enter comment here
245 delete fHistEtNonRemovedLambdas; // enter comment here
246 delete fHistEtNonRemovedElectrons; // enter comment here
247 delete fHistEtNonRemovedPositrons; // enter comment here
248 delete fHistEtNonRemovedMuPlus; // enter comment here
249 delete fHistEtNonRemovedMuMinus; // enter comment here
250 delete fHistEtNonRemovedNeutrons; // enter comment here
251 delete fHistEtNonRemovedAntiNeutrons; // enter comment here
252 delete fHistEtNonRemovedGammas; // enter comment here
253 delete fHistEtNonRemovedGammasFromPi0; // enter comment here
255 delete fHistEtRemovedGammas; // enter comment here
256 delete fHistEtRemovedNeutrons; // enter comment here
257 delete fHistEtRemovedAntiNeutrons; // enter comment here
260 delete fHistMultNonRemovedProtons; // enter comment here
261 delete fHistMultNonRemovedAntiProtons; // enter comment here
262 delete fHistMultNonRemovedPiPlus; // enter comment here
263 delete fHistMultNonRemovedPiMinus; // enter comment here
264 delete fHistMultNonRemovedKaonPlus; // enter comment here
265 delete fHistMultNonRemovedKaonMinus; // enter comment here
266 delete fHistMultNonRemovedK0s; // enter comment here
267 delete fHistMultNonRemovedK0L; // enter comment here
268 delete fHistMultNonRemovedLambdas; // enter comment here
269 delete fHistMultNonRemovedElectrons; // enter comment here
270 delete fHistMultNonRemovedPositrons; // enter comment here
271 delete fHistMultNonRemovedMuPlus; // enter comment here
272 delete fHistMultNonRemovedMuMinus; // enter comment here
273 delete fHistMultNonRemovedNeutrons; // enter comment here
274 delete fHistMultNonRemovedAntiNeutrons; // enter comment here
275 delete fHistMultNonRemovedGammas; // enter comment here
277 delete fHistMultRemovedGammas; // enter comment here
278 delete fHistMultRemovedNeutrons; // enter comment here
279 delete fHistMultRemovedAntiNeutrons; // enter comment here
281 delete fHistTrackMultvsNonRemovedCharged; // enter comment here
282 delete fHistTrackMultvsNonRemovedNeutral; // enter comment here
283 delete fHistTrackMultvsRemovedGamma; // enter comment here
285 delete fHistClusterMultvsNonRemovedCharged; // enter comment here
286 delete fHistClusterMultvsNonRemovedNeutral; // enter comment here
287 delete fHistClusterMultvsRemovedGamma; // enter comment here
289 delete fHistMultvsNonRemovedChargedE; // enter comment here
290 delete fHistMultvsNonRemovedNeutralE; // enter comment here
291 delete fHistMultvsRemovedGammaE; // enter comment here
292 delete fHistDxDzNonRemovedCharged; // enter comment here
293 delete fHistDxDzRemovedCharged; // enter comment here
294 delete fHistDxDzNonRemovedNeutral; // enter comment here
295 delete fHistDxDzRemovedNeutral; // enter comment here
297 delete fHistPiPlusMult; // enter comment here
298 delete fHistPiMinusMult; // enter comment here
299 delete fHistPiZeroMult; // enter comment here
301 delete fHistPiPlusMultAcc; // enter comment here
302 delete fHistPiMinusMultAcc; // enter comment here
303 delete fHistPiZeroMultAcc; // enter comment here
306 Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
307 { // analyse MC event
315 fCentClass = fCentrality->GetCentralityClass10(fCentralityMethod);
319 // Get us an mc event
321 AliFatal("ERROR: Event does not exist");
324 AliMCEvent *event = dynamic_cast<AliMCEvent*>(ev);
326 AliFatal("ERROR: MC Event does not exist");
330 Double_t protonMass =fgProtonMass;
333 AliGenEventHeader* genHeader = event->GenEventHeader();
335 Printf("ERROR: Event generation header does not exist");
338 AliGenHijingEventHeader* hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
339 if (hijingGenHeader) {
340 fImpactParameter = hijingGenHeader->ImpactParameter();
341 fNcoll = hijingGenHeader->HardScatters(); // or should this be some combination of NN() NNw() NwN() NwNw() ?
342 fNpart = hijingGenHeader->ProjectileParticipants() + hijingGenHeader->TargetParticipants();
345 // Let's play with the stack!
346 AliStack *stack = event->Stack();
348 Int_t nPrim = stack->GetNtrack();
351 for (Int_t iPart = 0; iPart < nPrim; iPart++)
354 TParticle *part = stack->Particle(iPart);
358 Printf("ERROR: Could not get particle %d", iPart);
362 TParticlePDG *pdg = part->GetPDG(0);
366 Printf("ERROR: Could not get particle PDG %d", iPart);
370 Double_t particleMassPart = 0; //The mass part in the Et calculation for this particle
371 Int_t code = pdg->PdgCode();
373 // Check for reasonable (for now neutral and singly charged) charge on the particle
374 //TODO:Maybe not only singly charged?
375 if (fSelector->CutNeutralMcParticle(iPart,*stack,*pdg))
379 //PrintFamilyTree(iPart, stack);
381 if (TMath::Abs(part->Eta()) < fCuts->GetCommonEtaCut())
383 //Printf("Particle with eta: %f, pid: %d", part->Eta(), code);
386 TMath::Abs(code) == fgProtonCode ||
387 TMath::Abs(code) == fgNeutronCode ||
388 TMath::Abs(code) == fgLambdaCode ||
389 TMath::Abs(code) == fgXiCode ||
390 TMath::Abs(code) == fgXi0Code ||
391 TMath::Abs(code) == fgOmegaCode
395 particleMassPart = - protonMass;
398 particleMassPart = protonMass;
401 Double_t et = part->Energy() * TMath::Sin(part->Theta()) + particleMassPart;
404 // Fill up total E_T counters for each particle species
405 if (code == fgProtonCode || code == fgAntiProtonCode)
408 if (code == fgPiPlusCode || code == fgPiMinusCode)
410 if (code == fgPiPlusCode)
417 if (code == fgGammaCode)
420 if (code == fgKPlusCode)
423 if(code == fgKMinusCode)
426 if (code == fgMuPlusCode || code == fgMuMinusCode)
429 if (code == fgEPlusCode || code == fgEMinusCode)
432 // some neutrals also
433 if (code == fgNeutronCode)
436 if (code == fgAntiNeutronCode)
439 if (code == fgGammaCode)
444 //if (TMath::Abs(pdg->Charge() - fCuts->GetMonteCarloNeutralParticle()) <1e-3 )
446 if(code == fgGammaCode || code == fgPi0Code || code == fgEtaCode)
449 //Printf("Gamma, phi: %f, eta: %f, phi cut min: %f, phi cut max: %f, eta cut: %f", part->Phi(), part->Eta(), fPhiMinCutAcc, fPhiMaxCutAcc, fEtaCutAcc);
450 //if (et > fCuts->GetCommonClusterEnergyCut()) fTotNeutralEt += et;
452 // inside EMCal acceptance
454 //if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiMaxCutAcc && part->Phi() > fPhiMinCutAcc)
456 if(fSelector->CutGeometricalAcceptance(*part))
458 fNeutralMultiplicity++;
460 if(fSelector->CutMinEnergy(*part))
462 fTotNeutralEtAfterMinEnergyCut += et;
464 if (part->Energy() > 0.05) partCount++;
468 else if (TMath::Abs( pdg->Charge() - fCuts->GetMonteCarloNeutralParticle())>1e-3 )
471 // inside EMCal acceptance
472 if (fSelector->CutGeometricalAcceptance(*part))
475 fChargedMultiplicity++;
479 if (code == fgProtonCode || code == fgAntiProtonCode)
482 if (code == fgPiPlusCode || code == fgPiMinusCode)
485 if (code == fgKPlusCode || code == fgKMinusCode)
488 if (code == fgMuPlusCode || code == fgMuMinusCode)
492 if (code == fgEPlusCode || code == fgEMinusCode)
494 fTotNeutralEt += et; // calling electrons neutral
497 } // inside EMCal acceptance
499 if (TrackHitsCalorimeter(part)) // magnetic field info not filled?
501 if (pdg->Charge() > 0) fHistPhivsPtPos->Fill(part->Phi(),part->Pt());
502 else fHistPhivsPtNeg->Fill(part->Phi(), part->Pt());
510 fTotEt = fTotChargedEt + fTotNeutralEt;
511 //fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;//
512 //std::cout << "Event done! # of particles: " << partCount << std::endl;
515 //Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliMCEvent* mcEvent,AliESDEvent* realEvent)
516 Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
517 { // analyse MC and real event info
518 //if(!mcEvent || !realEvent){
520 AliFatal("ERROR: Event does not exist");
523 AliAnalysisEt::AnalyseEvent(ev);
524 AliMCEvent *mcEvent = dynamic_cast<AliMCEvent*>(ev);
525 AliESDEvent *realEvent = dynamic_cast<AliESDEvent*>(ev2);
526 if (!mcEvent || !realEvent) {
527 AliFatal("ERROR: mcEvent or realEvent does not exist");
531 fSelector->SetEvent(realEvent);
535 AliStack *stack = mcEvent->Stack();
537 // get all detector clusters
538 // TRefArray* caloClusters = new TRefArray();
540 // if (fDetector == fCuts->GetDetectorEmcal()) realEvent->GetEMCALClusters( caloClusters );
541 //else if (fDetector == fCuts->GetDetectorPhos()) realEvent->GetPHOSClusters( caloClusters );
543 //AliFatal("Detector ID has not been specified");
547 TRefArray *caloClusters = fSelector->GetClusters();
549 Int_t nCluster = caloClusters->GetEntries();
552 for (int iCluster = 0; iCluster < nCluster; iCluster++ )
555 AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
556 //Float_t caloE = caloCluster->E();
557 fCutFlow->Fill(cf++);
558 if(!fSelector->CutDistanceToBadChannel(*caloCluster)) continue;
561 UInt_t iPart = (UInt_t)TMath::Abs(caloCluster->GetLabel());
562 TParticle *part = stack->Particle(iPart);
566 Printf("No MC particle %d", iCluster);
571 if (!stack->IsPhysicalPrimary(iPart)) // check whether particle is primary. we keep secondary electron and gamma for testing.
573 primIdx = GetPrimMother(iPart, stack);
576 std::cout << "What!? No primary?" << std::endl;
580 } // end of primary particle check
581 const int primCode = stack->Particle(primIdx)->GetPdgCode();
582 TParticlePDG *pdg = part->GetPDG();
585 Printf("ERROR: Could not get particle PDG %d", iPart);
589 Int_t code = pdg->PdgCode();
590 Double_t clEt = CalculateTransverseEnergy(caloCluster);
591 // if(code == fgK0SCode) std::cout << "K0 energy: " << caloCluster->E() << std::endl;
592 if(!fSelector->CutMinEnergy(*caloCluster)) continue;
593 fCutFlow->Fill(cf++);
596 caloCluster->GetPosition(pos);
599 TParticle *primPart = stack->Particle(primIdx);
600 fPrimaryCode = primPart->GetPdgCode();
601 fPrimaryCharge = primPart->GetPDG()->Charge();
603 fPrimaryE = primPart->Energy();
604 fPrimaryEt = primPart->Energy()*TMath::Sin(primPart->Theta());
606 fPrimaryPx = primPart->Px();
607 fPrimaryPy = primPart->Py();
608 fPrimaryPz = primPart->Pz();
610 fPrimaryVx = primPart->Vx();
611 fPrimaryVy = primPart->Vy();
612 fPrimaryVz = primPart->Vz();
614 fPrimaryAccepted = false;
616 fDepositedCode = part->GetPdgCode();
617 fDepositedEt = caloCluster->E() * TMath::Sin(cp.Eta());
618 fDepositedCharge = part->GetPDG()->Charge();
620 fDepositedVx = part->Vx();
621 fDepositedVy = part->Vy();
622 fDepositedVz = part->Vz();
623 //if(caloCluster->GetEmcCpvDistance() < fTrackDistanceCut)
624 pdg = part->GetPDG(0);
625 //if (TMath::Abs(caloCluster->GetTrackDx()) < fTrackDxCut && TMath::Abs(caloCluster->GetTrackDz()) < fTrackDzCut)
626 if(!fSelector->CutTrackMatching(*caloCluster))
627 //if(caloCluster->GetTrackMatchedIndex() > -1 && (fCuts->GetPhosTrackRCut() > TestCPV(caloCluster->GetTrackDx(), caloCluster->GetTrackDz(), ev->GetTrack(caloCluster->GetTrackMatchedIndex())->Pt(), ev->GetTrack(caloCluster->GetTrackMatchedIndex())->Charge(), ev)))
630 if (pdg->Charge() != 0)
632 //std::cout << "Removing charged: " << code << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
634 //fEnergyChargedRemoved += caloCluster->E();
635 fEnergyChargedRemoved += clEt;
636 fHistRemovedOrNot->Fill(0.0, fCentClass);
637 //std::cout << part->Vx() << " " << part->Vy() << " " << part->Vz() << " " << std::endl;
638 //fHistDecayVertexRemovedCharged->Fill(part->Vy(), part->Vx(), part->Vz());
639 fHistDxDzRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
640 if(code == fgProtonCode)
642 fProtonRemovedEt += clEt;
643 fProtonRemovedMult++;
646 else if(code == fgAntiProtonCode)
648 fAntiProtonRemovedEt += clEt;
649 fAntiProtonRemovedMult++;
651 else if(code == fgPiPlusCode)
653 fPiPlusRemovedEt += clEt;
654 fPiPlusRemovedMult++;
656 else if(code == fgPiMinusCode)
658 fPiMinusRemovedEt += clEt;
659 fPiMinusRemovedMult++;
661 else if(code == fgKPlusCode)
663 fKPlusRemovedEt += clEt;
666 else if(code == fgKMinusCode)
668 fKMinusRemovedEt += clEt;
669 fKMinusRemovedMult++;
671 else if(code == fgMuMinusCode)
673 fMuMinusRemovedEt += clEt;
674 fMuMinusRemovedMult++;
676 else if(code == fgMuPlusCode)
678 fMuPlusRemovedEt += clEt;
679 fMuPlusRemovedMult++;
681 else if(code == fgEMinusCode)
683 fEMinusRemovedEt += clEt;
684 fEMinusRemovedMult++;
686 else if(code == fgEPlusCode)
688 fEPlusRemovedEt += clEt;
693 Printf("Removed: %d, with E_T: %f", code, clEt);
695 if (primCode == fgGammaCode)
697 fGammaRemovedEt += clEt;
700 else if (primCode == fgPi0Code)
702 fPi0RemovedEt += clEt;
705 else if (primCode == fgEtaCode)
707 fPi0RemovedEt += clEt;
710 else if(primCode == fgK0LCode)
712 fK0lRemovedEt += clEt;
715 else if(primCode == fgK0SCode)
717 fK0sRemovedEt += clEt;
724 //std::cout << "Removing neutral: " << code << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
725 //std::cout << "Mother is: " << stack->Particle(part->GetMother(0))->GetPdgCode() << std::endl;
727 //fEnergyNeutralRemoved += caloCluster->E();
728 fEnergyNeutralRemoved += clEt;
729 fHistRemovedOrNot->Fill(1.0, fCentClass);
730 //std::cout << part->Vx() << " " << part->Vy() << " " << part->Vz() << " " << std::endl;
731 //fHistDecayVertexRemovedNeutral->Fill(part->Vy(), part->Vx(), part->Vz());
732 fHistDxDzRemovedNeutral->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
734 if (code == fgGammaCode || code == fgPi0Code || code == fgEtaCode)
736 fEtRemovedGammas += clEt;
737 fMultRemovedGammas++;
739 else if (code == fgNeutronCode)
741 fEtRemovedNeutrons += clEt;
742 fMultRemovedNeutrons++;
743 fNeutronRemovedEt += clEt;
744 fNeutronRemovedMult++;
746 else if (code == fgAntiNeutronCode)
748 fEtRemovedAntiNeutrons += clEt;
749 fMultRemovedAntiNeutrons++;
750 fAntiNeutronRemovedEt += clEt;
751 fAntiNeutronRemovedMult++;
753 if (primCode == fgGammaCode)
755 fGammaRemovedEt += clEt;
758 else if(primCode == fgPi0Code)
760 fPi0RemovedEt += clEt;
763 else if(primCode == fgEtaCode)
765 fPi0RemovedEt += clEt;
768 else if(primCode == fgK0LCode)
770 fK0lRemovedEt += clEt;
773 else if(primCode == fgK0SCode)
775 fK0sRemovedEt += clEt;
782 fCutFlow->Fill(cf++);
783 if (pdg->Charge() != 0)
785 //std::cout << "Not removing charged: " << code << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
786 //std::cout << "Mother is: " << stack->Particle(part->GetMother(0))->GetPdgCode() << std::endl;
787 fChargedNotRemoved++;
788 //fEnergyChargedNotRemoved += caloCluster->E();
789 fEnergyChargedNotRemoved += clEt;
790 fHistRemovedOrNot->Fill(2.0, fCentClass);
791 //std::cout << fHistDecayVertexNonRemovedCharged << std::endl;
792 //std::cout << part->Vx() << " " << part->Vy() << " " << part->Vz() << " " << std::endl;
793 //fHistDecayVertexNonRemovedCharged->Fill(part->Vy(), part->Vx(), part->Vz());
794 fHistDxDzNonRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
795 if (code == fgProtonCode)
797 //std::cout << clEt << std::endl;
798 fEtNonRemovedProtons += clEt;
799 fMultNonRemovedProtons++;
803 else if (code == fgAntiProtonCode)
805 //std::cout << clEt << std::endl;
806 fEtNonRemovedAntiProtons += clEt;
807 fMultNonRemovedAntiProtons++;
808 fAntiProtonEt += clEt;
811 else if (code == fgPiPlusCode)
813 //std::cout << "PI+" << clEt << std::endl;
814 fEtNonRemovedPiPlus += clEt;
815 fMultNonRemovedPiPlus++;
819 else if (code == fgPiMinusCode)
821 // std::cout << "PI-" << clEt << std::endl;
822 fEtNonRemovedPiMinus += clEt;
823 fMultNonRemovedPiMinus++;
827 else if (code == fgKPlusCode)
829 //std::cout << clEt << std::endl;
830 fEtNonRemovedKaonPlus += clEt;
831 fMultNonRemovedKaonPlus++;
835 else if (code == fgKMinusCode)
837 //std::cout << clEt << std::endl;
838 fEtNonRemovedKaonMinus += clEt;
839 fMultNonRemovedKaonMinus++;
843 else if (code == fgEPlusCode)
845 //std::cout << clEt << std::endl;
846 if (TMath::Sqrt(part->Vx()*part->Vx() + part->Vy()*part->Vy()) < 440)
848 fEtNonRemovedPositrons += clEt;
849 fMultNonRemovedPositrons++;
854 else if (code == fgEMinusCode)
856 //std::cout << clEt << std::endl;
857 if (TMath::Sqrt(part->Vx()*part->Vx() + part->Vy()*part->Vy()) < 440)
859 fEtNonRemovedElectrons += clEt;
860 fMultNonRemovedElectrons++;
865 else if (code == fgMuPlusCode)
867 fEtNonRemovedMuPlus += clEt;
868 fMultNonRemovedMuPlus++;
872 else if (code == fgMuMinusCode)
874 fEtNonRemovedMuMinus += clEt;
875 fMultNonRemovedMuMinus++;
879 if (primCode == fgGammaCode)
884 else if(primCode == fgPi0Code)
889 else if(primCode == fgEtaCode)
894 else if(primCode == fgK0LCode)
899 else if(primCode == fgK0SCode)
908 fPrimaryAccepted = true;
909 //std::cout << "Not removing charged: " << code << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
910 //std::cout << "Mother is: " << stack->Particle(part->GetMother(0))->GetPdgCode() << stack->Particle(part->GetMother(0))->GetPDG()->GetName() << ", E: " << part->Energy() << std::endl;
911 fNeutralNotRemoved++;
912 fEnergyNeutralNotRemoved += clEt;
913 fHistRemovedOrNot->Fill(3.0, fCentClass);
914 fHistDxDzNonRemovedNeutral->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
916 if(TMath::Abs(part->Vx()) < 1.0 && TMath::Abs(part->Vy()) < 1.0 && TMath::Abs(part->Vz()) < 20 && fSelector->IsEmEtParticle(primCode))
918 fTotEtWithSecondaryRemoved += clEt;
920 else if(fSelector->IsEmEtParticle(primCode))
922 fTotEtSecondaryFromEmEtPrimary += clEt;
925 fTotEtSecondary += clEt;
927 //code = stack->Particle(primIdx)->GetPdgCode();
928 if (code == fgGammaCode)
930 fEtNonRemovedGammas += clEt;
931 fMultNonRemovedGammas++;
934 // if (TMath::Abs(pdgMom->PdgCode()) == fgPi0Code || TMath::Abs(pdgMom->PdgCode()) == fgEtaCode || TMath::Abs(pdgMom->PdgCode()) == 331)
936 // std::cout << "Mother of gamma: " << pdgMom->PdgCode() << " " << pdgMom->GetName() << ", E: " << part->Energy() << std::endl;
937 // fEtNonRemovedGammasFromPi0 += clEt;
941 else if(TMath::Abs(code) == fgPi0Code)
943 fEtNonRemovedGammasFromPi0 += clEt;
944 fMultNonRemovedGammas++;
946 else if(TMath::Abs(code) == fgEtaCode)
948 fEtNonRemovedGammasFromPi0 += clEt;
949 fMultNonRemovedGammas++;
951 else if(TMath::Abs(code) == 331)
953 fEtNonRemovedGammasFromPi0 += clEt;
954 fMultNonRemovedGammas++;
956 else if (code == fgNeutronCode)
958 fEtNonRemovedNeutrons += clEt;
959 fMultNonRemovedNeutrons++;
961 else if (code == fgAntiNeutronCode)
963 fEtNonRemovedAntiNeutrons += clEt;
964 fMultNonRemovedAntiNeutrons++;
966 //else if (code == fgK0LCode || pdgMom->PdgCode() == fgK0SCode)
967 else if (code == fgK0SCode)
969 //std::cout << "K0 with energy: " << clEt << std::endl;
970 fEtNonRemovedK0S += clEt;
971 fMultNonRemovedK0s++;
973 else if(TMath::Abs(code) == fgK0LCode)
976 fEtNonRemovedK0L += clEt;
977 fMultNonRemovedK0L++;
980 else if (TMath::Abs(code) == fgLambdaCode)
982 fEtNonRemovedLambdas += clEt;
983 fMultNonRemovedLambdas++;
985 else std::cout << "Hmm, what is this (neutral not removed): " << code << " " << pdg->GetName() << ", ET: " << clEt << std::endl;
986 if (primCode == fgGammaCode)
991 else if(primCode == fgPi0Code)
996 else if(primCode == fgEtaCode)
1001 else if(primCode == fgK0LCode)
1006 else if(primCode == fgK0SCode)
1013 fPrimaryTree->Fill();
1014 } // end of loop over clusters
1016 //std::cout << "Distance cut: " << fTrackDistanceCut << std::endl;
1017 //std::cout << "Number of removed neutrals: " << fNeutralRemoved << std::endl;
1018 //std::cout << "Number of removed charged: " << fChargedRemoved << std::endl;
1019 //std::cout << "Number of non-removed charged: " << fChargedNotRemoved << std::endl;
1020 //std::cout << "Number of non-removed neutral: " << fNeutralNotRemoved << std::endl;
1022 // std::cout << "Energy of removed neutrals: " << fEnergyNeutralRemoved << std::endl;
1023 // std::cout << "Energy of removed charged: " << fEnergyChargedRemoved << std::endl;
1024 // std::cout << "Energy of non-removed charged: " << fEnergyChargedNotRemoved << std::endl;
1025 // std::cout << "Energy of non-removed neutral: " << fEnergyNeutralNotRemoved << std::endl;
1031 void AliAnalysisEtMonteCarlo::Init()
1033 AliAnalysisEt::Init();
1036 void AliAnalysisEtMonteCarlo::ResetEventValues()
1037 { // reset event values
1038 AliAnalysisEt::ResetEventValues();
1040 fTotEtSecondary = 0;
1041 fTotEtSecondaryFromEmEtPrimary = 0;
1042 fTotEtWithSecondaryRemoved = 0;
1044 // collision geometry defaults for p+p:
1045 fImpactParameter = 0;
1049 fEtNonRemovedProtons = 0;
1050 fEtNonRemovedAntiProtons = 0;
1051 fEtNonRemovedPiPlus = 0;
1052 fEtNonRemovedPiMinus = 0;
1053 fEtNonRemovedKaonPlus = 0;
1054 fEtNonRemovedKaonMinus = 0;
1055 fEtNonRemovedK0S = 0;
1056 fEtNonRemovedK0L = 0;
1057 fEtNonRemovedLambdas = 0;
1058 fEtNonRemovedElectrons = 0;
1059 fEtNonRemovedPositrons = 0;
1060 fEtNonRemovedMuPlus = 0;
1061 fEtNonRemovedMuMinus = 0;
1062 fEtNonRemovedNeutrons = 0;
1063 fEtNonRemovedAntiNeutrons = 0;
1064 fEtNonRemovedGammas = 0;
1065 fEtNonRemovedGammasFromPi0 = 0;
1067 fEtRemovedProtons = 0;
1068 fEtRemovedAntiProtons = 0;
1069 fEtRemovedPiPlus = 0;
1070 fEtRemovedPiMinus = 0;
1071 fEtRemovedKaonPlus = 0;
1072 fEtRemovedKaonMinus = 0;
1075 fEtRemovedLambdas = 0;
1076 fEtRemovedElectrons = 0;
1077 fEtRemovedPositrons = 0;
1078 fEtRemovedMuPlus = 0;
1079 fEtRemovedMuMinus = 0;
1080 fEtRemovedNeutrons = 0;
1082 fEtRemovedGammasFromPi0 = 0;
1083 fEtRemovedGammas = 0;
1084 fEtRemovedNeutrons = 0;
1085 fEtRemovedAntiNeutrons = 0;
1087 fMultNonRemovedProtons = 0;
1088 fMultNonRemovedAntiProtons = 0;
1089 fMultNonRemovedPiPlus = 0;
1090 fMultNonRemovedPiMinus = 0;
1091 fMultNonRemovedKaonPlus = 0;
1092 fMultNonRemovedKaonMinus = 0;
1093 fMultNonRemovedK0s = 0;
1094 fMultNonRemovedK0L = 0;
1095 fMultNonRemovedLambdas = 0;
1096 fMultNonRemovedElectrons = 0;
1097 fMultNonRemovedPositrons = 0;
1098 fMultNonRemovedMuPlus = 0;
1099 fMultNonRemovedMuMinus = 0;
1100 fMultNonRemovedNeutrons = 0;
1101 fMultNonRemovedAntiNeutrons = 0;
1102 fMultNonRemovedGammas = 0;
1104 fMultRemovedProtons = 0;
1105 fMultRemovedAntiProtons = 0;
1106 fMultRemovedPiPlus = 0;
1107 fMultRemovedPiMinus = 0;
1108 fMultRemovedKaonPlus = 0;
1109 fMultRemovedKaonMinus = 0;
1110 fMultRemovedK0s = 0;
1111 fMultRemovedK0L = 0;
1112 fMultRemovedLambdas = 0;
1113 fMultRemovedElectrons = 0;
1114 fMultRemovedPositrons = 0;
1115 fMultRemovedMuPlus = 0;
1116 fMultRemovedMuMinus = 0;
1118 fMultRemovedGammas = 0;
1119 fMultRemovedNeutrons = 0;
1120 fMultRemovedAntiNeutrons = 0;
1122 fEnergyChargedNotRemoved = 0;
1123 fEnergyChargedRemoved = 0;
1124 fEnergyNeutralNotRemoved = 0;
1125 fEnergyNeutralRemoved = 0;
1127 fChargedNotRemoved = 0;
1128 fChargedRemoved = 0;
1129 fNeutralNotRemoved = 0;
1130 fNeutralRemoved = 0;
1133 fTrackMultInAcc = 0;
1135 fTotNeutralEtAfterMinEnergyCut = 0;
1139 void AliAnalysisEtMonteCarlo::CreateHistograms()
1140 { // histogram related additions
1141 AliAnalysisEt::CreateHistograms();
1142 if (fEventSummaryTree) {
1143 fEventSummaryTree->Branch("fImpactParameter",&fImpactParameter,"fImpactParameter/D");
1144 fEventSummaryTree->Branch("fNcoll",&fNcoll,"fNcoll/I");
1145 fEventSummaryTree->Branch("fNpart",&fNpart,"fNpart/I");
1146 fEventSummaryTree->Branch("fTotEtWithSecondaryRemoved", &fTotEtWithSecondaryRemoved, "fTotEtWithSecondaryRemoved/D");
1147 fEventSummaryTree->Branch("fTotEtSecondaryFromEmEtPrimary", &fTotEtSecondaryFromEmEtPrimary, "fTotEtSecondaryFromEmEtPrimary/D");
1148 fEventSummaryTree->Branch("fTotEtSecondary", &fTotEtSecondary, "fTotEtSecondary/D");
1149 fEventSummaryTree->Branch("fTotNeutralEtAfterMinEnergyCut", &fTotNeutralEtAfterMinEnergyCut, "fTotNeutralEtAfterMinEnergyCut/D");
1152 //fHistDecayVertexNonRemovedCharged = new TH3F("fHistDecayVertexNonRemovedCharged","fHistDecayVertexNonRemovedCharged", 500, -470, 30, 500, -300, 300, 40, -20, 20);
1153 //fHistDecayVertexRemovedCharged = new TH3F("fHistDecayVertexRemovedCharged","fHistDecayVertexRemovedCharged", 500, -470, 30, 500, -300, 300, 40, -20, 20);
1154 //fHistDecayVertexNonRemovedNeutral = new TH3F("fHistDecayVertexNonRemovedNeutral","fHistDecayVertexNonRemovedNeutral", 500, -470, 30, 500, -300, 300, 40, -20, 20);
1155 //fHistDecayVertexRemovedNeutral = new TH3F("fHistDecayVertexRemovedNeutral","fHistDecayVertexRemovedNeutral", 500, -470, 30, 500, -300, 300, 40, -20, 20);
1157 fHistRemovedOrNot = new TH2F("fHistRemovedOrNot", "fHistRemovedOrNot", 4, -0.5, 3.5, 10, -0.5, 9.5);
1159 fHistEtNonRemovedProtons = new TH2F("fHistEtNonRemovedProtons", "fHistEtNonRemovedProtons", 1500, 0, 30, 10, -0.5, 9.5);
1160 fHistEtNonRemovedAntiProtons = new TH2F("fHistEtNonRemovedAntiProtons", "fHistEtNonRemovedAntiProtons", 1500, 0, 30, 10, -0.5, 9.5);
1161 fHistEtNonRemovedPiPlus = new TH2F("fHistEtNonRemovedPiPlus", "fHistEtNonRemovedPiPlus", 1500, 0, 30, 10, -0.5, 9.5);
1162 fHistEtNonRemovedPiMinus = new TH2F("fHistEtNonRemovedPiMinus", "fHistEtNonRemovedPiMinus", 1500, 0, 30, 10, -0.5, 9.5);
1163 fHistEtNonRemovedKaonPlus = new TH2F("fHistEtNonRemovedKaonPlus", "fHistEtNonRemovedKaonPlus", 1500, 0, 30, 10, -0.5, 9.5);
1164 fHistEtNonRemovedKaonMinus = new TH2F("fHistEtNonRemovedKaonMinus", "fHistEtNonRemovedKaonMinus", 1500, 0, 30, 10, -0.5, 9.5);
1165 fHistEtNonRemovedK0s = new TH2F("fHistEtNonRemovedK0s", "fHistEtNonRemovedK0s", 1500, 0, 30, 10, -0.5, 9.5);
1166 fHistEtNonRemovedK0L = new TH2F("fHistEtNonRemovedK0L", "fHistEtNonRemovedK0L", 1500, 0, 30, 10, -0.5, 9.5);
1167 fHistEtNonRemovedLambdas = new TH2F("fHistEtNonRemovedLambdas", "fHistEtNonRemovedLambdas", 1500, 0, 30, 10, -0.5, 9.5);
1168 fHistEtNonRemovedElectrons = new TH2F("fHistEtNonRemovedElectrons", "fHistEtNonRemovedElectrons", 1500, 0, 30, 10, -0.5, 9.5);
1169 fHistEtNonRemovedPositrons = new TH2F("fHistEtNonRemovedPositrons", "fHistEtNonRemovedPositrons", 1500, 0, 30, 10, -0.5, 9.5);
1170 fHistEtNonRemovedMuPlus = new TH2F("fHistEtNonRemovedMuPlus", "fHistEtNonRemovedMuPlus", 1500, 0, 30, 10, -0.5, 9.5);
1171 fHistEtNonRemovedMuMinus = new TH2F("fHistEtNonRemovedMuMinus", "fHistEtNonRemovedMuMinus", 1500, 0, 30, 10, -0.5, 9.5);
1172 fHistEtNonRemovedNeutrons = new TH2F("fHistEtNonRemovedNeutrons", "fHistEtNonRemovedNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
1173 fHistEtNonRemovedAntiNeutrons = new TH2F("fHistEtNonRemovedAntiNeutrons", "fHistEtNonRemovedAntiNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
1175 fHistEtNonRemovedGammas = new TH2F("fHistEtNonRemovedGammas", "fHistEtNonRemovedGammas", 1500, 0, 30, 10, -0.5, 9.5);
1176 fHistEtNonRemovedGammasFromPi0 = new TH2F("fHistEtNonRemovedGammasFromPi0", "fHistEtNonRemovedGammasFromPi0", 1500, 0, 30, 10, -0.5, 9.5);
1178 fHistEtRemovedGammas = new TH2F("fHistEtRemovedGammas", "fHistEtRemovedGammas", 1500, 0, 30, 10, -0.5, 9.5);
1179 fHistEtRemovedNeutrons = new TH2F("fHistEtRemovedNeutrons", "fHistEtRemovedNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
1180 fHistEtRemovedAntiNeutrons = new TH2F("fHistEtRemovedAntiNeutrons", "fHistEtRemovedAntiNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
1182 fHistEtRemovedCharged = new TH2F("fHistEtRemovedCharged", "fHistEtRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
1183 fHistEtRemovedNeutrals = new TH2F("fHistEtRemovedNeutrals", "fHistEtRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
1185 fHistEtNonRemovedCharged = new TH2F("fHistEtNonRemovedCharged", "fHistEtNonRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
1186 fHistEtNonRemovedNeutrals = new TH2F("fHistEtNonRemovedNeutrals", "fHistEtNonRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
1188 fHistMultNonRemovedProtons = new TH2F("fHistMultNonRemovedProtons", "fHistMultNonRemovedProtons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1189 fHistMultNonRemovedAntiProtons = new TH2F("fHistMultNonRemovedAntiProtons", "fHistMultNonRemovedAntiProtons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1190 fHistMultNonRemovedPiPlus = new TH2F("fHistMultNonRemovedPiPlus", "fHistMultNonRemovedPiPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1191 fHistMultNonRemovedPiMinus = new TH2F("fHistMultNonRemovedPiMinus", "fHistMultNonRemovedPiMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1192 fHistMultNonRemovedKaonPlus = new TH2F("fHistMultNonRemovedKaonPlus", "fHistMultNonRemovedKaonPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1193 fHistMultNonRemovedKaonMinus = new TH2F("fHistMultNonRemovedKaonMinus", "fHistMultNonRemovedKaonMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1194 fHistMultNonRemovedK0s = new TH2F("fHistMultNonRemovedK0s", "fHistMultNonRemovedK0s", 100, -0.5, 99.5, 10, -0.5, 9.5);
1195 fHistMultNonRemovedK0L = new TH2F("fHistMultNonRemovedK0L", "fHistMultNonRemovedK0L", 100, -0.5, 99.5, 10, -0.5, 9.5);
1196 fHistMultNonRemovedLambdas = new TH2F("fHistMultNonRemovedLambdas", "fHistMultNonRemovedLambdas", 100, -0.5, 99.5, 10, -0.5, 9.5);
1197 fHistMultNonRemovedElectrons = new TH2F("fHistMultNonRemovedElectrons", "fHistMultNonRemovedElectrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1198 fHistMultNonRemovedPositrons = new TH2F("fHistMultNonRemovedPositrons", "fHistMultNonRemovedPositrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1199 fHistMultNonRemovedMuPlus = new TH2F("fHistMultNonRemovedMuPlus", "fHistMultNonRemovedMuPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1200 fHistMultNonRemovedMuMinus = new TH2F("fHistMultNonRemovedMuMinus", "fHistMultNonRemovedMuMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
1201 fHistMultNonRemovedNeutrons = new TH2F("fHistMultNonRemovedNeutrons", "fHistMultNonRemovedNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1202 fHistMultNonRemovedAntiNeutrons = new TH2F("fHistMultNonRemovedAntiNeutrons", "fHistMultNonRemovedAntiNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1204 fHistMultNonRemovedGammas = new TH2F("fHistMultNonRemovedGammas", "fHistMultNonRemovedGammas", 100, -0.5, 99.5, 100, -0.5, 99.5);
1206 fHistMultRemovedGammas = new TH2F("fHistMultRemovedGammas", "fHistMultRemovedGammas", 100, -0.5, 99.5, 100, -0.5, 99.5);
1207 fHistMultRemovedNeutrons = new TH2F("fHistMultRemovedNeutrons", "fHistMultRemovedNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1208 fHistMultRemovedAntiNeutrons = new TH2F("fHistMultRemovedAntiNeutrons", "fHistMultRemovedAntiNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
1210 fHistMultRemovedCharged = new TH2F("fHistMultRemovedCharged", "fHistMultRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
1211 fHistMultRemovedNeutrals = new TH2F("fHistMultRemovedNeutrals", "fHistMultRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
1213 fHistMultNonRemovedCharged = new TH2F("fHistMultNonRemovedCharged", "fHistMultNonRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
1214 fHistMultNonRemovedNeutrals = new TH2F("fHistMultNonRemovedNeutrals", "fHistMultNonRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);*/
1217 fHistMultRemovedCharged = new TH2F("fHistMultRemovedCharged", "fHistMultRemovedCharged", 100, -0.5, 99.5, 100, -0.5, 99.5);
1218 fHistMultRemovedNeutrals = new TH2F("fHistMultRemovedNeutrals", "fHistMultRemovedNeutrals", 100, -0.5, 99.5, 100, -0.5, 99.5);
1220 fHistMultNonRemovedCharged = new TH2F("fHistMultNonRemovedCharged", "fHistMultNonRemovedCharged", 100, -0.5, 99.5, 100, -0.5, 99.5);
1221 fHistMultNonRemovedNeutrals = new TH2F("fHistMultNonRemovedNeutrals", "fHistMultNonRemovedNeutrals", 100, -0.5, 99.5, 100, -0.5, 99.5);
1223 fHistTrackMultvsNonRemovedCharged = new TH2F("fHistTrackMultvsNonRemovedCharged", "fHistTrackMultvsNonRemovedCharged", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1224 fHistTrackMultvsNonRemovedNeutral = new TH2F("fHistTrackMultvsNonRemovedNeutral", "fHistTrackMultvsNonRemovedNeutral", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1225 fHistTrackMultvsRemovedGamma = new TH2F("fHistTrackMultvsRemovedGamma", "fHistTrackMultvsRemovedGamma", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1227 fHistClusterMultvsNonRemovedCharged = new TH2F("fHistClusterMultvsNonRemovedCharged", "fHistClusterMultvsNonRemovedCharged", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1228 fHistClusterMultvsNonRemovedNeutral = new TH2F("fHistClusterMultvsNonRemovedNeutral", "fHistClusterMultvsNonRemovedNeutral", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1229 fHistClusterMultvsRemovedGamma = new TH2F("fHistClusterMultvsRemovedGamma", "fHistClusterMultvsRemovedGamma", 1000, -0.5, 999.5, 100, -0.5, 99.5);
1231 fHistDxDzNonRemovedCharged = new TH2F("fHistDxDzNonRemovedCharged", "fHistDxDzNonRemovedCharged", 800, -200, 200, 800, -200, 200);
1232 fHistDxDzRemovedCharged = new TH2F("fHistDxDzRemovedCharged", "fHistDxDzRemovedCharged", 800, -200, 200, 800, -200, 200);
1233 fHistDxDzNonRemovedNeutral = new TH2F("fHistDxDzNonRemovedNeutral", "fHistDxDzNonRemovedNeutral", 800, -200, 200, 800, -200, 200);
1234 fHistDxDzRemovedNeutral = new TH2F("fHistDxDzRemovedNeutral", "fHistDxDzRemovedNeutral", 800, -200, 200, 800, -200, 200);
1236 fHistPiPlusMult = new TH1F("fHistPiPlusMult", "fHistPiPlusMult", 2000, -0.5, 1999.5);
1237 fHistPiMinusMult = new TH1F("fHistPiMinusMult", "fHistPiMinusMult", 2000, -0.5, 1999.5);
1238 fHistPiZeroMult = new TH1F("fHistPiZeroMult", "fHistPiZeroMult", 2000, -0.5, 1999.5);
1240 fHistPiPlusMultAcc = new TH1F("fHistPiPlusMultAcc", "fHistPiPlusMultAcc", 2000, -0.5, 1999.5);
1241 fHistPiMinusMultAcc = new TH1F("fHistPiMinusMultAcc", "fHistPiMinusMultAcc", 2000, -0.5, 1999.5);
1242 fHistPiZeroMultAcc = new TH1F("fHistPiZeroMultAcc", "fHistPiZeroMultAcc", 2000, -0.5, 1999.5);
1244 if(fCuts->GetHistMakeTree())
1246 TString treename = "fPrimaryTree" + fHistogramNameSuffix;
1247 fPrimaryTree = new TTree(treename, treename);
1249 fPrimaryTree->Branch("fTotEt",&fTotEt,"fTotEt/D");
1250 fPrimaryTree->Branch("fNeutralMultiplicity", &fNeutralMultiplicity, "fNeutralMultiplicity/I");
1251 fPrimaryTree->Branch("fCentClass",&fCentClass,"fCentClass/I");
1253 fPrimaryTree->Branch("fPrimaryCode", &fPrimaryCode, "fPrimaryCode/I");
1254 fPrimaryTree->Branch("fPrimaryCharge", &fPrimaryCharge, "fPrimaryCharge/I");
1256 fPrimaryTree->Branch("fPrimaryE", &fPrimaryE, "fPrimaryE/D");
1257 fPrimaryTree->Branch("fPrimaryEt", &fPrimaryEt, "fPrimaryEt/D");
1259 fPrimaryTree->Branch("fPrimaryPx", &fPrimaryPx, "fPrimaryPx/D");
1260 fPrimaryTree->Branch("fPrimaryPy", &fPrimaryPy, "fPrimaryPy/D");
1261 fPrimaryTree->Branch("fPrimaryPz", &fPrimaryPz, "fPrimaryPz/D");
1263 fPrimaryTree->Branch("fPrimaryVx", &fPrimaryVx, "fPrimaryVx/D");
1264 fPrimaryTree->Branch("fPrimaryVy", &fPrimaryVy, "fPrimaryVy/D");
1265 fPrimaryTree->Branch("fPrimaryVz", &fPrimaryVz, "fPrimaryVz/D");
1267 fPrimaryTree->Branch("fPrimaryAccepted", &fPrimaryAccepted, "fPrimaryAccepted/B");
1270 fPrimaryTree->Branch("fDepositedCode", &fDepositedCode, "fDepositedCode/I");
1271 fPrimaryTree->Branch("fDepositedCharge", &fDepositedCharge, "fDepositedCharge/I");
1272 fPrimaryTree->Branch("fDepositedEt", &fDepositedEt, "fDepositedEt/D");
1274 fPrimaryTree->Branch("fDepositedVx", &fDepositedVx, "fDepositedVx/D");
1275 fPrimaryTree->Branch("fDepositedVy", &fDepositedVy, "fDepositedVy/D");
1276 fPrimaryTree->Branch("fDepositedVz", &fDepositedVz, "fDepositedVz/D");
1282 void AliAnalysisEtMonteCarlo::FillOutputList(TList *list)
1283 { //fill the output list
1284 AliAnalysisEt::FillOutputList(list);
1286 if(fCuts->GetHistMakeTree())
1288 list->Add(fPrimaryTree);
1291 list->Add(fHistRemovedOrNot);
1293 list->Add(fHistEtNonRemovedProtons);
1294 list->Add(fHistEtNonRemovedAntiProtons);
1295 list->Add(fHistEtNonRemovedPiPlus);
1296 list->Add(fHistEtNonRemovedPiMinus);
1297 list->Add(fHistEtNonRemovedKaonPlus);
1298 list->Add(fHistEtNonRemovedKaonMinus);
1299 list->Add(fHistEtNonRemovedK0s);
1300 list->Add(fHistEtNonRemovedK0L);
1301 list->Add(fHistEtNonRemovedLambdas);
1302 list->Add(fHistEtNonRemovedElectrons);
1303 list->Add(fHistEtNonRemovedPositrons);
1304 list->Add(fHistEtNonRemovedMuPlus);
1305 list->Add(fHistEtNonRemovedMuMinus);
1306 list->Add(fHistEtNonRemovedNeutrons);
1307 list->Add(fHistEtNonRemovedAntiNeutrons);
1308 list->Add(fHistEtNonRemovedGammas);
1309 list->Add(fHistEtNonRemovedGammasFromPi0);
1311 list->Add(fHistEtRemovedGammas);
1312 list->Add(fHistEtRemovedNeutrons);
1313 list->Add(fHistEtRemovedAntiNeutrons);
1315 list->Add(fHistEtRemovedCharged);
1316 list->Add(fHistEtRemovedNeutrals);
1318 list->Add(fHistEtNonRemovedCharged);
1319 list->Add(fHistEtNonRemovedNeutrals);
1321 list->Add(fHistMultNonRemovedProtons);
1322 list->Add(fHistMultNonRemovedAntiProtons);
1323 list->Add(fHistMultNonRemovedPiPlus);
1324 list->Add(fHistMultNonRemovedPiMinus);
1325 list->Add(fHistMultNonRemovedKaonPlus);
1326 list->Add(fHistMultNonRemovedKaonMinus);
1327 list->Add(fHistMultNonRemovedK0s);
1328 list->Add(fHistMultNonRemovedK0L);
1329 list->Add(fHistMultNonRemovedLambdas);
1330 list->Add(fHistMultNonRemovedElectrons);
1331 list->Add(fHistMultNonRemovedPositrons);
1332 list->Add(fHistMultNonRemovedMuPlus);
1333 list->Add(fHistMultNonRemovedMuMinus);
1334 list->Add(fHistMultNonRemovedNeutrons);
1335 list->Add(fHistMultNonRemovedAntiNeutrons);
1336 list->Add(fHistMultNonRemovedGammas);
1338 list->Add(fHistMultRemovedGammas);
1339 list->Add(fHistMultRemovedNeutrons);
1340 list->Add(fHistMultRemovedAntiNeutrons);
1342 list->Add(fHistMultRemovedCharged);
1343 list->Add(fHistMultRemovedNeutrals);
1345 list->Add(fHistMultNonRemovedCharged);
1346 list->Add(fHistMultNonRemovedNeutrals);
1348 list->Add(fHistTrackMultvsNonRemovedCharged);
1349 list->Add(fHistTrackMultvsNonRemovedNeutral);
1350 list->Add(fHistTrackMultvsRemovedGamma);
1352 list->Add(fHistClusterMultvsNonRemovedCharged);
1353 list->Add(fHistClusterMultvsNonRemovedNeutral);
1354 list->Add(fHistClusterMultvsRemovedGamma);
1356 //list->Add(fHistDecayVertexNonRemovedCharged);
1357 //list->Add(fHistDecayVertexNonRemovedNeutral);
1358 //list->Add(fHistDecayVertexRemovedCharged);
1359 //list->Add(fHistDecayVertexRemovedNeutral);
1361 list->Add(fHistDxDzNonRemovedCharged);
1362 list->Add(fHistDxDzRemovedCharged);
1363 list->Add(fHistDxDzNonRemovedNeutral);
1364 list->Add(fHistDxDzRemovedNeutral);
1366 list->Add(fHistPiPlusMult);
1367 list->Add(fHistPiMinusMult);
1368 list->Add(fHistPiZeroMult);
1369 list->Add(fHistPiPlusMultAcc);
1370 list->Add(fHistPiMinusMultAcc);
1371 list->Add(fHistPiZeroMultAcc);
1376 bool AliAnalysisEtMonteCarlo::TrackHitsCalorimeter(TParticle* part, Double_t magField)
1378 // printf(" TrackHitsCalorimeter - magField %f\n", magField);
1379 AliESDtrack *esdTrack = new AliESDtrack(part);
1380 // Printf("MC Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
1382 Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
1384 // if(prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
1386 bool status = prop && fSelector->CutGeometricalAcceptance(*esdTrack);
1392 void AliAnalysisEtMonteCarlo::FillHistograms()
1393 { // let base class fill its histograms, and us fill the local ones
1394 AliAnalysisEt::FillHistograms();
1395 //std::cout << fEtNonRemovedPiPlus << " " << fCentClass << std::endl;
1397 fHistEtNonRemovedProtons->Fill(fEtNonRemovedProtons, fCentClass);
1398 fHistEtNonRemovedAntiProtons->Fill(fEtNonRemovedAntiProtons, fCentClass);
1399 fHistEtNonRemovedKaonPlus->Fill(fEtNonRemovedKaonPlus, fCentClass);
1400 fHistEtNonRemovedKaonMinus->Fill(fEtNonRemovedKaonMinus, fCentClass);
1401 fHistEtNonRemovedK0s->Fill(fEtNonRemovedK0S, fCentClass);
1402 fHistEtNonRemovedK0L->Fill(fEtNonRemovedK0L, fCentClass);
1403 fHistEtNonRemovedLambdas->Fill(fEtNonRemovedLambdas, fCentClass);
1404 fHistEtNonRemovedPiPlus->Fill(fEtNonRemovedPiPlus, fCentClass);
1405 fHistEtNonRemovedPiMinus->Fill(fEtNonRemovedPiMinus, fCentClass);
1406 fHistEtNonRemovedElectrons->Fill(fEtNonRemovedElectrons, fCentClass);
1407 fHistEtNonRemovedPositrons->Fill(fEtNonRemovedPositrons, fCentClass);
1408 fHistEtNonRemovedMuPlus->Fill(fEtNonRemovedMuPlus, fCentClass);
1409 fHistEtNonRemovedMuMinus->Fill(fEtNonRemovedMuMinus, fCentClass);
1410 fHistEtNonRemovedNeutrons->Fill(fEtNonRemovedNeutrons, fCentClass);
1411 fHistEtNonRemovedAntiNeutrons->Fill(fEtNonRemovedAntiNeutrons, fCentClass);
1412 fHistEtNonRemovedGammas->Fill(fEtNonRemovedGammas, fCentClass);
1413 fHistEtNonRemovedGammasFromPi0->Fill(fEtNonRemovedGammasFromPi0, fCentClass);
1415 fHistEtRemovedGammas->Fill(fEtRemovedGammas, fNClusters);
1416 fHistEtRemovedNeutrons->Fill(fEtRemovedNeutrons, fCentClass);
1417 fHistEtRemovedAntiNeutrons->Fill(fEtRemovedAntiNeutrons, fCentClass);
1419 // fHistEtRemovedCharged->Fill(fEtRemovedAntiProtons+fEtRemovedElectrons+fEtRemovedKaonMinus+fEtRemovedKaonPlus
1420 // +fEtRemovedMuMinus+fEtRemovedMuPlus+fEtRemovedPiMinus+fEtRemovedPiPlus+fEtRemovedPositrons
1421 // +fEtRemovedProtons.
1423 // fHistEtRemovedNeutrals->Fill(fEtRemovedNeutrons+fEtRemovedAntiNeutrons, fCentClass);
1425 // fHistEtNonRemovedCharged->Fill(fEtNonRemovedAntiProtons+fEtNonRemovedElectrons+fEtNonRemovedKaonMinus+fEtNonRemovedKaonPlus
1426 // +fEtNonRemovedMuMinus+fEtNonRemovedMuPlus+fEtNonRemovedPiMinus+fEtNonRemovedPiPlus+fEtNonRemovedPositrons
1427 // +fEtNonRemovedProtons,
1429 // fHistEtRemovedNeutrals->Fill(fEtNonRemovedNeutrons+fEtNonRemovedAntiNeutrons, fCentClass);
1431 fHistEtRemovedCharged->Fill(fEnergyChargedRemoved, fNClusters);
1432 fHistEtRemovedNeutrals->Fill(fEnergyNeutralRemoved, fNClusters);
1433 fHistEtNonRemovedCharged->Fill(fEnergyChargedNotRemoved, fNClusters);
1434 fHistEtNonRemovedNeutrals->Fill(fEnergyNeutralNotRemoved, fNClusters);
1436 fHistMultRemovedCharged->Fill(fChargedRemoved, fNClusters);
1437 fHistMultRemovedNeutrals->Fill(fNeutralRemoved, fNClusters);
1438 fHistMultNonRemovedCharged->Fill(fChargedNotRemoved, fNClusters);
1439 fHistMultNonRemovedNeutrals->Fill(fNeutralNotRemoved, fNClusters);
1442 fHistMultNonRemovedProtons->Fill(fMultNonRemovedProtons, fCentClass);
1443 fHistMultNonRemovedAntiProtons->Fill(fMultNonRemovedAntiProtons, fCentClass);
1444 fHistMultNonRemovedKaonPlus->Fill(fMultNonRemovedKaonPlus, fCentClass);
1445 fHistMultNonRemovedKaonMinus->Fill(fMultNonRemovedKaonMinus, fCentClass);
1446 fHistMultNonRemovedK0s->Fill(fMultNonRemovedK0s, fCentClass);
1447 fHistMultNonRemovedK0L->Fill(fMultNonRemovedK0L, fCentClass);
1448 fHistMultNonRemovedLambdas->Fill(fMultNonRemovedLambdas, fCentClass);
1449 fHistMultNonRemovedPiPlus->Fill(fMultNonRemovedPiPlus, fCentClass);
1450 fHistMultNonRemovedPiMinus->Fill(fMultNonRemovedPiMinus, fCentClass);
1451 fHistMultNonRemovedElectrons->Fill(fMultNonRemovedElectrons, fCentClass);
1452 fHistMultNonRemovedPositrons->Fill(fMultNonRemovedPositrons, fCentClass);
1453 fHistMultNonRemovedMuPlus->Fill(fMultNonRemovedMuPlus, fCentClass);
1454 fHistMultNonRemovedMuMinus->Fill(fMultNonRemovedMuMinus, fCentClass);
1455 fHistMultNonRemovedNeutrons->Fill(fMultNonRemovedNeutrons, fCentClass);
1456 fHistMultNonRemovedAntiNeutrons->Fill(fMultNonRemovedAntiNeutrons, fCentClass);
1457 fHistMultNonRemovedGammas->Fill(fMultNonRemovedGammas, fCentClass);
1459 fHistMultRemovedGammas->Fill(fMultRemovedGammas, fCentClass);
1460 fHistMultRemovedNeutrons->Fill(fMultRemovedNeutrons, fCentClass);
1461 fHistMultRemovedAntiNeutrons->Fill(fMultRemovedAntiNeutrons, fCentClass);
1463 fHistTrackMultvsNonRemovedCharged->Fill(fTrackMultInAcc,
1464 fMultNonRemovedAntiProtons+fMultNonRemovedElectrons+fMultNonRemovedKaonMinus+fMultNonRemovedKaonPlus
1465 +fMultNonRemovedMuMinus+fMultNonRemovedMuPlus+fMultNonRemovedPiMinus+fMultNonRemovedPiPlus+fMultNonRemovedPositrons
1466 +fMultNonRemovedProtons);
1468 fHistTrackMultvsNonRemovedNeutral->Fill(fTrackMultInAcc,
1469 fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedK0L+fMultNonRemovedLambdas);
1471 fHistTrackMultvsRemovedGamma->Fill(fTrackMultInAcc,
1472 fMultRemovedGammas);
1474 fHistClusterMultvsNonRemovedCharged->Fill(fNClusters,
1475 fMultNonRemovedAntiProtons+fMultNonRemovedElectrons+fMultNonRemovedKaonMinus
1476 +fMultNonRemovedKaonPlus+fMultNonRemovedMuMinus+fMultNonRemovedMuPlus
1477 +fMultNonRemovedPiMinus+fMultNonRemovedPiPlus+fMultNonRemovedPositrons+fMultNonRemovedProtons);
1479 fHistClusterMultvsNonRemovedNeutral->Fill(fNClusters,
1480 fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedK0L+fMultNonRemovedLambdas);
1482 fHistClusterMultvsRemovedGamma->Fill(fNClusters,
1483 fMultRemovedGammas);
1490 Int_t AliAnalysisEtMonteCarlo::PrintFamilyTree(Int_t partIdx, AliStack* stack)
1491 { // print family tree
1492 TParticle *part = stack->Particle(partIdx);
1493 if(part->GetPdgCode() == fgK0SCode)
1495 std::cout << "This is index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") , is it primary: " << stack->IsPhysicalPrimary(partIdx)<< std::endl;
1496 std::cout << "PID: " << part->GetPdgCode() << "/" << part->GetName() << std::endl;
1497 std::cout << "Energy: " << part->Energy() << std::endl;
1498 std::cout << "Vertex: " << part->Vx() << ", " << part->Vy() << ", " << part->Vz() << std::endl;
1500 return PrintMothers(partIdx, stack, 1);
1503 Int_t AliAnalysisEtMonteCarlo::PrintMothers(Int_t partIdx, AliStack* stack, Int_t gen)
1505 char *tabs = new char[gen+1];
1506 for(Int_t i = 0; i < gen; ++i)
1508 //std::cout << i << std::endl;
1512 Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
1518 TParticle *mother = stack->Particle(mothIdx);
1519 if(mother->GetPdgCode() == fgK0SCode)
1521 std::cout << tabs << "Mother of index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") is: " << mothIdx << ", is it primary: " << stack->IsPhysicalPrimary(mothIdx)<< std::endl;
1522 std::cout << tabs << "PID: " << mother->GetPdgCode() << "/" << mother->GetName() << std::endl;
1523 std::cout << tabs << "Energy: " << mother->Energy() << std::endl;
1524 std::cout << tabs << "Vertex: " << mother->Vx() << ", " << mother->Vy() << ", " << mother->Vz() << std::endl;
1526 if(mother->GetPdgCode() == fgK0SCode)
1528 // std::cout << "K0S!!!!!!!!!!!!!11111!!!!!" << std::endl;
1530 // std::cout << "Mother of index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") is: " << mothIdx << std::endl;
1531 // std::cout << "PID: " << mother->GetPdgCode() << "/" << mother->GetName() << std::endl;
1532 // std::cout << "Energy: " << mother->Energy() << std::endl;
1533 // std::cout << "Vertex: " << mother->Vx() << ", " << mother->Vy() << ", " << mother->Vz() << std::endl;
1536 return PrintMothers(mothIdx, stack, gen+1) + 1;
1539 Int_t AliAnalysisEtMonteCarlo::GetPrimMother(Int_t partIdx, AliStack *stack)
1540 { // get primary mother
1543 Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
1544 if(mothIdx < 0) return -1;
1545 TParticle *mother = stack->Particle(mothIdx);
1548 // if(mother->GetPdgCode() == fgK0SCode)
1550 // std::cout << "!!!!!!!!!!!!!!!!! K0S !!!!!!!!!!!!!!!!!!" << std::endl;
1553 //if(mother->GetPdgCode() == fgK0SCode&& stack->IsPhysicalPrimary(mothIdx))
1555 // std::cout << "!!!!!!!!!!!!!!!!! Primary K0S !!!!!!!!!!!!!!!!!!" << std::endl;
1558 if(stack->IsPhysicalPrimary(mothIdx)) return mothIdx;
1559 else return GetPrimMother(mothIdx, stack);
1569 Int_t AliAnalysisEtMonteCarlo::GetK0InFamily(Int_t partIdx, AliStack* stack)
1570 { // get K0 in family
1573 if(stack->Particle(partIdx)->GetPdgCode() == fgK0SCode) return partIdx;
1574 Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
1575 if(mothIdx < 0) return -1;
1576 TParticle *mother = stack->Particle(mothIdx);
1579 if(mother->GetPdgCode() == fgK0SCode)
1583 return GetK0InFamily(mothIdx, stack);