1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 //_________________________________________________________________________
18 // Class for reading data (Kinematics) in order to do prompt gamma
19 // or other particle identification and correlations
21 //*-- Author: Gustavo Conesa (LNF-INFN)
22 //////////////////////////////////////////////////////////////////////////////
25 // --- ROOT system ---
27 #include <TParticle.h>
29 #include <TClonesArray.h>
30 //#include "Riostream.h"
32 //---- ANALYSIS system ----
33 #include "AliCaloTrackMCReader.h"
35 #include "AliGenEventHeader.h"
37 #include "AliAODCaloCluster.h"
38 #include "AliAODTrack.h"
40 ClassImp(AliCaloTrackMCReader)
42 //____________________________________________________________________________
43 AliCaloTrackMCReader::AliCaloTrackMCReader() :
44 AliCaloTrackReader(), fDecayPi0(0),
45 fNeutralParticlesArray(0x0), fChargedParticlesArray(0x0),
46 fStatusArray(0x0), fKeepAllStatus(0), fClonesArrayType(0)
50 //Initialize parameters
56 //____________________________________________________________________________
57 AliCaloTrackMCReader::AliCaloTrackMCReader(const AliCaloTrackMCReader & g) :
58 AliCaloTrackReader(g), fDecayPi0(g.fDecayPi0),
59 fNeutralParticlesArray(g.fNeutralParticlesArray?new TArrayI(*g.fNeutralParticlesArray):0x0),
60 fChargedParticlesArray(g.fChargedParticlesArray?new TArrayI(*g.fChargedParticlesArray):0x0),
61 fStatusArray(g.fStatusArray?new TArrayI(*g.fStatusArray):0x0),
62 fKeepAllStatus(g.fKeepAllStatus), fClonesArrayType(g.fClonesArrayType)
67 //_________________________________________________________________________
68 AliCaloTrackMCReader & AliCaloTrackMCReader::operator = (const AliCaloTrackMCReader & source)
70 // assignment operator
72 if(&source == this) return *this;
74 fDecayPi0 = source.fDecayPi0;
76 delete fChargedParticlesArray;
77 fChargedParticlesArray = source.fChargedParticlesArray?new TArrayI(*source.fChargedParticlesArray):0x0;
79 delete fNeutralParticlesArray;
80 fNeutralParticlesArray = source.fNeutralParticlesArray?new TArrayI(*source.fNeutralParticlesArray):0x0;
83 fStatusArray = source.fStatusArray?new TArrayI(*source.fStatusArray):0x0;
85 fKeepAllStatus = source.fKeepAllStatus ;
86 fClonesArrayType = source.fClonesArrayType ;
92 //_________________________________
93 AliCaloTrackMCReader::~AliCaloTrackMCReader() {
96 if(fChargedParticlesArray) delete fChargedParticlesArray ;
97 if(fNeutralParticlesArray) delete fNeutralParticlesArray ;
98 if(fStatusArray) delete fStatusArray ;
102 //____________________________________________________________________________
103 void AliCaloTrackMCReader::GetVertex(Double_t v[3]) const {
104 //Return vertex position
107 GetGenEventHeader()->PrimaryVertex(pv);
115 //_______________________________________________________________
116 void AliCaloTrackMCReader::InitParameters()
119 //Initialize the parameters of the analysis.
123 fChargedParticlesArray = new TArrayI(1);
124 fChargedParticlesArray->SetAt(11,0);
125 //Int_t pdgarray[]={12,14,16};// skip neutrinos
126 //fNeutralParticlesArray = new TArrayI(3, pdgarray);
127 fNeutralParticlesArray = new TArrayI(3);
128 fNeutralParticlesArray->SetAt(12,0); fNeutralParticlesArray->SetAt(14,1); fNeutralParticlesArray->SetAt(16,2);
129 fStatusArray = new TArrayI(1);
130 fStatusArray->SetAt(1,0);
132 fKeepAllStatus = kTRUE;
133 fClonesArrayType = kAliAOD ;
137 //____________________________________________________________________________
138 void AliCaloTrackMCReader::FillCalorimeters(const Int_t iParticle, TParticle* particle, TLorentzVector momentum,
139 Int_t &indexPHOS, Int_t &indexEMCAL) {
140 //Fill AODCaloClusters or TParticles lists of PHOS or EMCAL
142 if(fFillPHOS && fFidutialCut->IsInFidutialCut(momentum,"PHOS") && momentum.Pt() > fPHOSPtMin){
144 if(fClonesArrayType == kTParticle) new((*fAODPHOS)[indexPHOS++]) TParticle(*particle) ;
147 Char_t ttype= AliAODCluster::kPHOSNeutral;
148 Int_t labels[] = {iParticle};
149 Float_t x[] = {momentum.X(), momentum.Y(), momentum.Z()};
150 AliAODCaloCluster *calo = new((*fAODPHOS)[indexPHOS++])
151 AliAODCaloCluster(iParticle,1,labels,momentum.E(), x, NULL, ttype, 0);
152 SetCaloClusterPID(particle->GetPdgCode(),calo) ;
153 if(fDebug > 3 && momentum.Pt() > 0.2)
154 printf("Fill MC PHOS :: Selected tracks %s E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
155 particle->GetName(),momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
159 else if(fFillEMCAL && fFidutialCut->IsInFidutialCut(momentum,"EMCAL") && momentum.Pt() > fEMCALPtMin){
160 if(fClonesArrayType == kTParticle) new((*fAODEMCAL)[indexEMCAL++]) TParticle(*particle) ;
162 Char_t ttype= AliAODCluster::kEMCALClusterv1;
163 Int_t labels[] = {iParticle};
164 Float_t x[] = {momentum.X(), momentum.Y(), momentum.Z()};
165 AliAODCaloCluster *calo = new((*fAODEMCAL)[indexEMCAL++])
166 AliAODCaloCluster(iParticle,1,labels,momentum.E(), x, NULL, ttype, 0);
167 SetCaloClusterPID(particle->GetPdgCode(),calo) ;
168 if(fDebug > 3 && momentum.Pt() > 0.2)
169 printf("Fill MC EMCAL :: Selected tracks %s E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
170 particle->GetName(),momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
175 //____________________________________________________________________________
176 void AliCaloTrackMCReader::FillInputEvent()
178 //Create list of particles from EMCAL, PHOS and CTS.
180 if(fClonesArrayType == kTParticle){
181 fAODCTS = new TClonesArray("TParticle",0);
182 fAODEMCAL = new TClonesArray("TParticle",0);
183 fAODPHOS = new TClonesArray("TParticle",0);
185 else if(fClonesArrayType == kAliAOD){
186 fAODCTS = new TClonesArray("AliAODTrack",0);
187 fAODPHOS = new TClonesArray("AliAODCaloCluster",0);
188 fAODEMCAL = new TClonesArray("AliAODCaloCluster",0);
190 else {AliFatal("Wrong clones type");}
194 Int_t indexEMCAL = 0 ;
195 Int_t indexPHOS = 0 ;
197 Int_t iParticle = 0 ;
198 Double_t charge = 0.;
200 for (iParticle=0 ; iParticle < GetStack()->GetNtrack() ; iParticle++) {
201 TParticle * particle = GetStack()->Particle(iParticle);
202 TLorentzVector momentum;
205 Int_t pdg = particle->GetPdgCode();
207 //Keep particles with a given status
208 if(KeepParticleWithStatus(particle->GetStatusCode()) && (particle->Pt() > 0) ){
210 //Skip bizarre particles, they crash when charge is calculated
211 // if(TMath::Abs(pdg) == 3124 || TMath::Abs(pdg) > 10000000) continue ;
213 charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
214 particle->Momentum(momentum);
215 //---------- Charged particles ----------------------
216 if((charge != 0) && (momentum.Pt() > fCTSPtMin) && (fFidutialCut->IsInFidutialCut(momentum,"CTS"))){
218 //Particles in CTS acceptance
219 if(fDebug > 3 && momentum.Pt() > 0.2)
220 printf("Fill MC CTS :: Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
221 momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
223 if(fClonesArrayType == kTParticle) new((*fAODCTS)[indexCh++]) TParticle(*particle) ;
225 x[0] = particle->Vx(); x[1] = particle->Vy(); x[2] = particle->Vz();
226 p[0] = particle->Px(); p[1] = particle->Py(); p[2] = particle->Pz();
227 AliAODTrack *aodTrack = new((*fAODCTS)[indexCh++])
228 AliAODTrack(0, iParticle, p, kTRUE, x, kFALSE,NULL, 0, 0,
231 kFALSE, // No fit performed
232 kFALSE, // No fit performed
233 AliAODTrack::kPrimary,
235 SetTrackChargeAndPID(pdg, aodTrack);
238 //Keep some charged particles in calorimeters lists
239 if((fFillPHOS || fFillEMCAL) && KeepChargedParticles(pdg)) FillCalorimeters(iParticle, particle, momentum, indexPHOS, indexEMCAL);
243 //-------------Neutral particles ----------------------
244 else if(charge == 0 && (fFillPHOS || fFillEMCAL)){
245 //Skip neutrinos or other neutral particles
246 //if(SkipNeutralParticles(pdg) || particle->IsPrimary()) continue ; // protection added (MG)
247 if(SkipNeutralParticles(pdg)) continue ;
249 //Fill particle/calocluster arrays
251 FillCalorimeters(iParticle, particle, momentum, indexPHOS, indexEMCAL);
254 //Sometimes pi0 are stable for the generator, if needed decay it by hand
256 if(momentum.Pt() > fPHOSPtMin || momentum.Pt() > fEMCALPtMin){
257 TLorentzVector lvGamma1, lvGamma2 ;
258 //Double_t angle = 0;
261 MakePi0Decay(momentum,lvGamma1,lvGamma2);//,angle);
263 //Check if Pi0 is in the acceptance of the calorimeters, if aperture angle is small, keep it
264 TParticle * pPhoton1 = new TParticle(22,1,iParticle,0,0,0,lvGamma1.Px(),lvGamma1.Py(),
265 lvGamma1.Pz(),lvGamma1.E(),0,0,0,0);
266 TParticle * pPhoton2 = new TParticle(22,1,iParticle,0,0,0,lvGamma2.Px(),lvGamma2.Py(),
267 lvGamma2.Pz(),lvGamma2.E(),0,0,0,0);
268 //Fill particle/calocluster arrays
269 FillCalorimeters(iParticle,pPhoton1,lvGamma1, indexPHOS, indexEMCAL);
270 FillCalorimeters(iParticle,pPhoton2,lvGamma2, indexPHOS, indexEMCAL);
273 else FillCalorimeters(iParticle,particle, momentum, indexPHOS, indexEMCAL); //Add the rest
276 } //particle with correct status
281 //________________________________________________________________
282 Bool_t AliCaloTrackMCReader::KeepParticleWithStatus(Int_t status) const {
283 //Check if status is equal to one of the list
284 //These particles will be used in analysis.
286 for(Int_t i= 0; i < fStatusArray->GetSize(); i++)
287 if(status == fStatusArray->At(i)) return kTRUE ;
296 //________________________________________________________________
297 Bool_t AliCaloTrackMCReader::KeepChargedParticles(Int_t pdg) const {
298 //Check if pdg is equal to one of the charged particles list
299 //These particles will be added to the calorimeters lists.
301 for(Int_t i= 0; i < fChargedParticlesArray->GetSize(); i++)
302 if(TMath::Abs(pdg) == fChargedParticlesArray->At(i)) return kTRUE ;
308 //________________________________________________________________
309 void AliCaloTrackMCReader::Print(const Option_t * opt) const
311 //Print some relevant parameters set for the analysis
315 Info("**** Print **** ", "%s %s", GetName(), GetTitle() ) ;
317 printf("Decay Pi0? : %d\n", fDecayPi0) ;
318 printf("TClonesArray type : %d\n", fClonesArrayType) ;
319 printf("Keep all status? : %d\n", fKeepAllStatus) ;
321 if(!fKeepAllStatus) printf("Keep particles with status : ");
322 for(Int_t i= 0; i < fStatusArray->GetSize(); i++)
323 printf(" %d ; ", fStatusArray->At(i));
326 printf("Skip neutral particles in calo : ");
327 for(Int_t i= 0; i < fNeutralParticlesArray->GetSize(); i++)
328 printf(" %d ; ", fNeutralParticlesArray->At(i));
331 printf("Keep charged particles in calo : ");
332 for(Int_t i= 0; i < fChargedParticlesArray->GetSize(); i++)
333 printf(" %d ; ", fChargedParticlesArray->At(i));
338 //____________________________________________________________________________
339 void AliCaloTrackMCReader::MakePi0Decay(TLorentzVector &p0, TLorentzVector &p1,
340 TLorentzVector &p2) const {//, Double_t &angle) {
341 // Perform isotropic decay pi0 -> 2 photons
342 // p0 is pi0 4-momentum (inut)
343 // p1 and p2 are photon 4-momenta (output)
344 // cout<<"Boost vector"<<endl;
345 Double_t mPi0 = TDatabasePDG::Instance()->GetParticle(111)->Mass();
346 TVector3 b = p0.BoostVector();
347 //cout<<"Parameters"<<endl;
348 //Double_t mPi0 = p0.M();
349 Double_t phi = TMath::TwoPi() * gRandom->Rndm();
350 Double_t cosThe = 2 * gRandom->Rndm() - 1;
351 Double_t cosPhi = TMath::Cos(phi);
352 Double_t sinPhi = TMath::Sin(phi);
353 Double_t sinThe = TMath::Sqrt(1-cosThe*cosThe);
354 Double_t ePi0 = mPi0/2.;
355 //cout<<"ePi0 "<<ePi0<<endl;
356 //cout<<"Components"<<endl;
357 p1.SetPx(+ePi0*cosPhi*sinThe);
358 p1.SetPy(+ePi0*sinPhi*sinThe);
359 p1.SetPz(+ePi0*cosThe);
361 //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
362 //cout<<"p1 Mass: "<<p1.Px()*p1.Px()+p1.Py()*p1.Py()+p1.Pz()*p1.Pz()-p1.E()*p1.E()<<endl;
363 p2.SetPx(-ePi0*cosPhi*sinThe);
364 p2.SetPy(-ePi0*sinPhi*sinThe);
365 p2.SetPz(-ePi0*cosThe);
367 //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
368 //cout<<"p2 Mass: "<<p2.Px()*p2.Px()+p2.Py()*p2.Py()+p2.Pz()*p2.Pz()-p2.E()*p2.E()<<endl;
369 //cout<<"Boost "<<b.X()<<" "<<b.Y()<<" "<<b.Z()<<endl;
371 //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
373 //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
374 //cout<<"angle"<<endl;
375 //angle = p1.Angle(p2.Vect());
379 //____________________________________________________________________________
380 void AliCaloTrackMCReader::SetInputEvent(TObject* /*esd*/, TObject* aod, TObject* mc) {
381 // Connect the data pointer
382 SetMC((AliMCEvent*) mc);
383 SetAOD((AliAODEvent*) aod);
387 //________________________________________________________________
388 Bool_t AliCaloTrackMCReader::SkipNeutralParticles(Int_t pdg) const {
389 //Check if pdg is equal to one of the neutral particles list
390 //These particles will be skipped from analysis.
392 for(Int_t i= 0; i < fNeutralParticlesArray->GetSize(); i++)
393 if(TMath::Abs(pdg) == fNeutralParticlesArray->At(i)) return kTRUE ;
400 //____________________________________________________________________
401 void AliCaloTrackMCReader::SetTrackChargeAndPID(const Int_t pdgCode, AliAODTrack *track) const {
402 //Give a PID weight for tracks equal to 1 depending on the particle type
404 Float_t pid[10] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
410 pid[AliAODTrack::kUnknown] = 1.;
415 track->SetCharge(-1);
416 pid[AliAODTrack::kElectron] = 1.;
421 track->SetCharge(+1);
422 pid[AliAODTrack::kElectron] = 1.;
427 track->SetCharge(-1);
428 pid[AliAODTrack::kMuon] = 1.;
433 track->SetCharge(+1);
434 pid[AliAODTrack::kMuon] = 1.;
440 pid[AliAODTrack::kUnknown] = 1.;
445 track->SetCharge(+1);
446 pid[AliAODTrack::kPion] = 1.;
451 track->SetCharge(-1);
452 pid[AliAODTrack::kPion] = 1.;
458 pid[AliAODTrack::kUnknown] = 1.;
463 track->SetCharge(+1);
464 pid[AliAODTrack::kKaon] = 1.;
469 track->SetCharge(-1);
470 pid[AliAODTrack::kKaon] = 1.;
476 pid[AliAODTrack::kUnknown] = 1.;
481 track->SetCharge(+1);
482 pid[AliAODTrack::kProton] = 1.;
486 case -2212: // anti-p
487 track->SetCharge(-1);
488 pid[AliAODTrack::kProton] = 1.;
494 pid[AliAODTrack::kUnknown] = 1.;
500 pid[AliAODTrack::kUnknown] = 1.;
504 case -311: // anti-K0
506 pid[AliAODTrack::kUnknown] = 1.;
512 pid[AliAODTrack::kUnknown] = 1.;
518 pid[AliAODTrack::kUnknown] = 1.;
523 track->SetCharge(+1);
524 pid[AliAODTrack::kUnknown] = 1.;
529 track->SetCharge(-1);
530 pid[AliAODTrack::kUnknown] = 1.;
535 track->SetCharge(-1);
536 pid[AliAODTrack::kUnknown] = 1.;
542 pid[AliAODTrack::kUnknown] = 1.;
547 track->SetCharge(-1);
548 pid[AliAODTrack::kUnknown] = 1.;
553 track->SetCharge(-1);
554 pid[AliAODTrack::kUnknown] = 1.;
560 pid[AliAODTrack::kUnknown] = 1.;
564 case -3122: // anti-Lambda
566 pid[AliAODTrack::kUnknown] = 1.;
570 case -3222: // anti-Sigma-
571 track->SetCharge(-1);
572 pid[AliAODTrack::kUnknown] = 1.;
576 case -3212: // anti-Sigma0
578 pid[AliAODTrack::kUnknown] = 1.;
582 case -3112: // anti-Sigma+
583 track->SetCharge(+1);
584 pid[AliAODTrack::kUnknown] = 1.;
588 case -3322: // anti-Xi0
590 pid[AliAODTrack::kUnknown] = 1.;
594 case -3312: // anti-Xi+
595 track->SetCharge(+1);
598 case -3334: // anti-Omega+
599 track->SetCharge(+1);
600 pid[AliAODTrack::kUnknown] = 1.;
605 track->SetCharge(+1);
606 pid[AliAODTrack::kUnknown] = 1.;
611 track->SetCharge(-1);
612 pid[AliAODTrack::kUnknown] = 1.;
618 pid[AliAODTrack::kUnknown] = 1.;
622 case -421: // anti-D0
624 pid[AliAODTrack::kUnknown] = 1.;
629 track->SetCharge(-99);
630 pid[AliAODTrack::kUnknown] = 1.;
639 //____________________________________________________________________
640 void AliCaloTrackMCReader::SetCaloClusterPID(const Int_t pdgCode, AliAODCaloCluster *calo) const {
641 //Give a PID weight for CaloClusters equal to 1 depending on the particle type
643 Float_t pid[13] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
648 pid[AliAODCaloCluster::kPhoton] = 1.;
653 pid[AliAODCaloCluster::kElectron] = 1.;
658 pid[AliAODCaloCluster::kElectron] = 1.;
663 pid[AliAODCaloCluster::kCharged] = 1.;
668 pid[AliAODCaloCluster::kCharged] = 1.;
673 pid[AliAODCaloCluster::kPi0] = 1.;
678 pid[AliAODCaloCluster::kCharged] = 1.;
683 pid[AliAODCaloCluster::kCharged] = 1.;
688 pid[AliAODCaloCluster::kKaon0] = 1.;
689 pid[AliAODCaloCluster::kNeutral] = 1;
694 pid[AliAODCaloCluster::kCharged] = 1.;
699 pid[AliAODCaloCluster::kCharged] = 1.;
704 pid[AliAODCaloCluster::kNeutron] = 1.;
705 pid[AliAODCaloCluster::kNeutral] = 1.;
710 pid[AliAODCaloCluster::kCharged] = 1.;
714 case -2212: // anti-p
715 pid[AliAODCaloCluster::kCharged] = 1.;
720 pid[AliAODCaloCluster::kKaon0] = 1.;
721 pid[AliAODCaloCluster::kNeutral] = 1.;
726 pid[AliAODCaloCluster::kKaon0] = 1.;
727 pid[AliAODCaloCluster::kNeutral] = 1.;
731 case -311: // anti-K0
732 pid[AliAODCaloCluster::kKaon0] = 1.;
733 pid[AliAODCaloCluster::kNeutral] = 1.;
738 pid[AliAODCaloCluster::kNeutral] = 1.;
743 pid[AliAODCaloCluster::kUnknown] = 1.;
748 pid[AliAODCaloCluster::kUnknown] = 1.;
753 pid[AliAODCaloCluster::kUnknown] = 1.;
758 pid[AliAODCaloCluster::kUnknown] = 1.;
763 pid[AliAODCaloCluster::kUnknown] = 1.;
768 pid[AliAODCaloCluster::kUnknown] = 1.;
773 pid[AliAODCaloCluster::kUnknown] = 1.;
778 pid[AliAODCaloCluster::kNeutron] = 1.;
779 pid[AliAODCaloCluster::kNeutral] = 1.;
783 case -3122: // anti-Lambda
784 pid[AliAODCaloCluster::kUnknown] = 1.;
788 case -3222: // anti-Sigma-
789 pid[AliAODCaloCluster::kUnknown] = 1.;
793 case -3212: // anti-Sigma0
794 pid[AliAODCaloCluster::kUnknown] = 1.;
798 case -3112: // anti-Sigma+
799 pid[AliAODCaloCluster::kUnknown] = 1.;
803 case -3322: // anti-Xi0
804 pid[AliAODCaloCluster::kUnknown] = 1.;
808 case -3312: // anti-Xi+
809 pid[AliAODCaloCluster::kUnknown] = 1.;
813 case -3334: // anti-Omega+
814 pid[AliAODCaloCluster::kUnknown] = 1.;
819 pid[AliAODCaloCluster::kUnknown] = 1.;
824 pid[AliAODCaloCluster::kUnknown] = 1.;
829 pid[AliAODCaloCluster::kUnknown] = 1.;
833 case -421: // anti-D0
834 pid[AliAODCaloCluster::kUnknown] = 1.;
839 pid[AliAODCaloCluster::kUnknown] = 1.;