]>
Commit | Line | Data |
---|---|---|
1c5acb87 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
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 * | |
477d6cee | 12 | * about the suitability of this software for any purpose. It is * |
1c5acb87 | 13 | * provided "as is" without express or implied warranty. * |
14 | **************************************************************************/ | |
15 | /* $Id: $ */ | |
16 | ||
17 | //_________________________________________________________________________ | |
18 | // Class for reading data (Kinematics) in order to do prompt gamma | |
19 | // or other particle identification and correlations | |
20 | // | |
21 | //*-- Author: Gustavo Conesa (LNF-INFN) | |
22 | ////////////////////////////////////////////////////////////////////////////// | |
23 | ||
24 | ||
25 | // --- ROOT system --- | |
26 | ||
27 | #include <TParticle.h> | |
28 | #include <TRandom.h> | |
477d6cee | 29 | #include <TArrayI.h> |
1c5acb87 | 30 | //#include "Riostream.h" |
31 | ||
32 | //---- ANALYSIS system ---- | |
33 | #include "AliCaloTrackMCReader.h" | |
1c5acb87 | 34 | #include "AliGenEventHeader.h" |
35 | #include "AliStack.h" | |
36 | #include "AliAODCaloCluster.h" | |
37 | #include "AliAODTrack.h" | |
477d6cee | 38 | #include "AliFidutialCut.h" |
39 | #include "AliAODEvent.h" | |
40 | #include "TParticle.h" | |
1c5acb87 | 41 | |
477d6cee | 42 | ClassImp(AliCaloTrackMCReader) |
1c5acb87 | 43 | |
44 | //____________________________________________________________________________ | |
45 | AliCaloTrackMCReader::AliCaloTrackMCReader() : | |
46 | AliCaloTrackReader(), fDecayPi0(0), | |
477d6cee | 47 | fNeutralParticlesArray(0x0), fChargedParticlesArray(0x0), |
48 | fStatusArray(0x0), fKeepAllStatus(0), fCheckOverlap(0), | |
49 | fEMCALOverlapAngle(0),fPHOSOverlapAngle(0), fIndex2ndPhoton(0) | |
1c5acb87 | 50 | { |
51 | //Ctor | |
52 | ||
53 | //Initialize parameters | |
54 | InitParameters(); | |
55 | fDataType = kMC; | |
56 | ||
57 | } | |
58 | ||
59 | //____________________________________________________________________________ | |
60 | AliCaloTrackMCReader::AliCaloTrackMCReader(const AliCaloTrackMCReader & g) : | |
61 | AliCaloTrackReader(g), fDecayPi0(g.fDecayPi0), | |
62 | fNeutralParticlesArray(g.fNeutralParticlesArray?new TArrayI(*g.fNeutralParticlesArray):0x0), | |
63 | fChargedParticlesArray(g.fChargedParticlesArray?new TArrayI(*g.fChargedParticlesArray):0x0), | |
64 | fStatusArray(g.fStatusArray?new TArrayI(*g.fStatusArray):0x0), | |
477d6cee | 65 | fKeepAllStatus(g.fKeepAllStatus), fCheckOverlap(g.fCheckOverlap), |
691bdd02 | 66 | fEMCALOverlapAngle( g.fEMCALOverlapAngle), fPHOSOverlapAngle(g.fPHOSOverlapAngle), |
67 | fIndex2ndPhoton(g.fIndex2ndPhoton) | |
1c5acb87 | 68 | { |
69 | // cpy ctor | |
70 | } | |
71 | ||
72 | //_________________________________________________________________________ | |
73 | AliCaloTrackMCReader & AliCaloTrackMCReader::operator = (const AliCaloTrackMCReader & source) | |
74 | { | |
75 | // assignment operator | |
76 | ||
77 | if(&source == this) return *this; | |
78 | ||
79 | fDecayPi0 = source.fDecayPi0; | |
80 | ||
81 | delete fChargedParticlesArray; | |
82 | fChargedParticlesArray = source.fChargedParticlesArray?new TArrayI(*source.fChargedParticlesArray):0x0; | |
83 | ||
84 | delete fNeutralParticlesArray; | |
85 | fNeutralParticlesArray = source.fNeutralParticlesArray?new TArrayI(*source.fNeutralParticlesArray):0x0; | |
86 | ||
87 | delete fStatusArray; | |
88 | fStatusArray = source.fStatusArray?new TArrayI(*source.fStatusArray):0x0; | |
89 | ||
90 | fKeepAllStatus = source.fKeepAllStatus ; | |
1c5acb87 | 91 | |
92 | return *this; | |
93 | ||
94 | } | |
95 | ||
96 | //_________________________________ | |
97 | AliCaloTrackMCReader::~AliCaloTrackMCReader() { | |
98 | //Dtor | |
99 | ||
100 | if(fChargedParticlesArray) delete fChargedParticlesArray ; | |
101 | if(fNeutralParticlesArray) delete fNeutralParticlesArray ; | |
102 | if(fStatusArray) delete fStatusArray ; | |
103 | ||
104 | } | |
105 | ||
106 | //____________________________________________________________________________ | |
107 | void AliCaloTrackMCReader::GetVertex(Double_t v[3]) const { | |
108 | //Return vertex position | |
109 | ||
110 | TArrayF pv; | |
111 | GetGenEventHeader()->PrimaryVertex(pv); | |
112 | v[0]=pv.At(0); | |
113 | v[1]=pv.At(1); | |
114 | v[2]=pv.At(2); | |
115 | ||
116 | } | |
117 | ||
118 | ||
119 | //_______________________________________________________________ | |
120 | void AliCaloTrackMCReader::InitParameters() | |
121 | { | |
122 | ||
123 | //Initialize the parameters of the analysis. | |
124 | ||
125 | fDecayPi0 = kFALSE; | |
126 | ||
127 | fChargedParticlesArray = new TArrayI(1); | |
128 | fChargedParticlesArray->SetAt(11,0); | |
129 | //Int_t pdgarray[]={12,14,16};// skip neutrinos | |
130 | //fNeutralParticlesArray = new TArrayI(3, pdgarray); | |
131 | fNeutralParticlesArray = new TArrayI(3); | |
132 | fNeutralParticlesArray->SetAt(12,0); fNeutralParticlesArray->SetAt(14,1); fNeutralParticlesArray->SetAt(16,2); | |
133 | fStatusArray = new TArrayI(1); | |
134 | fStatusArray->SetAt(1,0); | |
135 | ||
136 | fKeepAllStatus = kTRUE; | |
1c5acb87 | 137 | |
691bdd02 | 138 | fCheckOverlap = kFALSE; |
139 | fEMCALOverlapAngle = 2.5 * TMath::DegToRad(); | |
140 | fPHOSOverlapAngle = 0.5 * TMath::DegToRad(); | |
141 | fIndex2ndPhoton = -1; | |
142 | } | |
143 | //____________________________________________________________________________ | |
144 | void AliCaloTrackMCReader::CheckOverlap(const Float_t anglethres, const Int_t imom, Int_t & iPrimary, Int_t & index, TLorentzVector & mom, Int_t & pdg) { | |
145 | //Check overlap of decay photons | |
146 | if( fIndex2ndPhoton==iPrimary ){ | |
147 | fIndex2ndPhoton=-1; | |
148 | return; | |
149 | } | |
150 | else fIndex2ndPhoton=-1; | |
151 | ||
152 | ||
153 | if(pdg!=22) return; | |
154 | ||
155 | TLorentzVector ph1, ph2; | |
156 | TParticle *meson = GetStack()->Particle(imom); | |
157 | Int_t mepdg = meson->GetPdgCode(); | |
158 | Int_t idaug1 = meson->GetFirstDaughter(); | |
159 | if((mepdg == 111 || mepdg == 221 ) && meson->GetNDaughters() == 2){ //Check only decay in 2 photons | |
160 | TParticle * d1 = GetStack()->Particle(idaug1); | |
161 | TParticle *d2 = GetStack()->Particle(idaug1+1); | |
162 | if(d1->GetPdgCode() == 22 && d2->GetPdgCode() == 22 ){ | |
163 | d1->Momentum(ph1); | |
164 | d2->Momentum(ph2); | |
165 | //printf("angle %2.2f\n",ph1.Angle(ph2.Vect())); | |
166 | ||
167 | if(anglethres > ph1.Angle(ph2.Vect())){ | |
168 | //Keep the meson | |
169 | pdg=mepdg; | |
170 | index=imom; | |
171 | meson->Momentum(mom); | |
172 | //printf("Overlap:: pt %2.2f, phi %2.2f, eta %2.2f\n",mom.Pt(),mom.Phi(),mom.Eta()); | |
173 | if(iPrimary == idaug1) iPrimary++; //skip next photon in list | |
174 | } | |
175 | else{ | |
176 | //Do not check overlapping for next decay photon from same meson | |
177 | if(iPrimary == idaug1) {fIndex2ndPhoton = idaug1+1; | |
178 | } | |
179 | ||
180 | } | |
181 | } | |
182 | }//Meson Decay with 2 photon daughters | |
1c5acb87 | 183 | } |
184 | ||
185 | //____________________________________________________________________________ | |
691bdd02 | 186 | void AliCaloTrackMCReader::FillCalorimeters(Int_t & iParticle, TParticle* particle, TLorentzVector momentum, |
477d6cee | 187 | Int_t &ncalo) { |
188 | //Fill AODCaloClusters or TParticles lists of PHOS or EMCAL | |
189 | //In PHOS | |
190 | if(fFillPHOS && fFidutialCut->IsInFidutialCut(momentum,"PHOS") && momentum.Pt() > fPHOSPtMin){ | |
191 | Int_t index = iParticle ; | |
192 | Int_t pdg = TMath::Abs(particle->GetPdgCode()); | |
193 | if(fCheckOverlap) | |
194 | CheckOverlap(fPHOSOverlapAngle,particle->GetFirstMother(),index, iParticle, momentum, pdg); | |
195 | ||
196 | Char_t ttype= AliAODCluster::kPHOSNeutral; | |
197 | Int_t labels[] = {index}; | |
198 | Float_t x[] = {momentum.X(), momentum.Y(), momentum.Z()}; | |
199 | //Create object and write it to file | |
200 | AliAODCaloCluster *calo = new((*(fOutputEvent->GetCaloClusters()))[ncalo++]) | |
201 | AliAODCaloCluster(index,1,labels,momentum.E(), x, NULL, ttype, 0); | |
202 | ||
203 | SetCaloClusterPID(pdg,calo) ; | |
204 | if(fDebug > 3 && momentum.Pt() > 0.2) | |
a3aebfff | 205 | printf("AliCaloTrackMCReader::FillCalorimeters() - PHOS : Selected cluster %s E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n", |
477d6cee | 206 | particle->GetName(),momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta()); |
207 | fAODPHOS->Add(calo);//reference the selected object to the list | |
208 | } | |
9415d854 | 209 | |
477d6cee | 210 | //In EMCAL |
9415d854 | 211 | if(fFillEMCAL && fFidutialCut->IsInFidutialCut(momentum,"EMCAL") && momentum.Pt() > fEMCALPtMin){ |
477d6cee | 212 | Int_t index = iParticle ; |
213 | Int_t pdg = TMath::Abs(particle->GetPdgCode()); | |
214 | //Int_t pdgorg=pdg; | |
215 | if(fCheckOverlap) | |
216 | CheckOverlap(fEMCALOverlapAngle,particle->GetFirstMother(),iParticle, index, momentum, pdg); | |
217 | ||
218 | Char_t ttype= AliAODCluster::kEMCALClusterv1; | |
219 | Int_t labels[] = {index}; | |
220 | Float_t x[] = {momentum.X(), momentum.Y(), momentum.Z()}; | |
221 | //Create object and write it to file | |
222 | AliAODCaloCluster *calo = new((*(fOutputEvent->GetCaloClusters()))[ncalo++]) | |
223 | AliAODCaloCluster(iParticle,1,labels,momentum.E(), x, NULL, ttype, 0); | |
224 | ||
225 | SetCaloClusterPID(pdg,calo) ; | |
226 | if(fDebug > 3 && momentum.Pt() > 0.2) | |
a3aebfff | 227 | printf("AliCaloTrackMCReader::FillCalorimeters() - EMCAL : Selected cluster %s E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n", |
477d6cee | 228 | particle->GetName(),momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta()); |
229 | fAODEMCAL->Add(calo);//reference the selected object to the list | |
230 | } | |
1c5acb87 | 231 | } |
232 | ||
233 | //____________________________________________________________________________ | |
6639984f | 234 | void AliCaloTrackMCReader::FillInputEvent(Int_t iEntry){ |
477d6cee | 235 | //Fill the event counter and input lists that are needed, called by the analysis maker. |
236 | ||
237 | fEventNumber = iEntry; | |
238 | Int_t iParticle = 0 ; | |
239 | Double_t charge = 0.; | |
240 | Int_t ncalo = (fOutputEvent->GetCaloClusters())->GetEntriesFast(); | |
241 | Int_t ntrack = (fOutputEvent->GetTracks())->GetEntriesFast(); | |
6639984f | 242 | |
477d6cee | 243 | for (iParticle = 0 ; iParticle < GetStack()->GetNtrack() ; iParticle++) { |
244 | TParticle * particle = GetStack()->Particle(iParticle); | |
245 | TLorentzVector momentum; | |
246 | Float_t p[3]; | |
247 | Float_t x[3]; | |
248 | Int_t pdg = particle->GetPdgCode(); | |
249 | ||
250 | //Keep particles with a given status | |
251 | if(KeepParticleWithStatus(particle->GetStatusCode()) && (particle->Pt() > 0) ){ | |
252 | ||
253 | //Skip bizarre particles, they crash when charge is calculated | |
254 | // if(TMath::Abs(pdg) == 3124 || TMath::Abs(pdg) > 10000000) continue ; | |
255 | ||
256 | charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge(); | |
257 | particle->Momentum(momentum); | |
258 | //---------- Charged particles ---------------------- | |
259 | if((charge != 0) && (momentum.Pt() > fCTSPtMin) && (fFidutialCut->IsInFidutialCut(momentum,"CTS"))){ | |
260 | if(fFillCTS){ | |
261 | //Particles in CTS acceptance | |
262 | if(fDebug > 3 && momentum.Pt() > 0.2) | |
a3aebfff | 263 | printf("AliCaloTrackMCReader::FillInputEvent() - CTS : Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n", |
477d6cee | 264 | momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta()); |
265 | ||
266 | x[0] = particle->Vx(); x[1] = particle->Vy(); x[2] = particle->Vz(); | |
267 | p[0] = particle->Px(); p[1] = particle->Py(); p[2] = particle->Pz(); | |
268 | //Create object and write it to file | |
269 | AliAODTrack *aodTrack = new((*(fOutputEvent->GetTracks()))[ntrack++]) | |
270 | AliAODTrack(0, iParticle, p, kTRUE, x, kFALSE,NULL, 0, 0, | |
271 | NULL, | |
272 | 0x0,//primary, | |
273 | kFALSE, // No fit performed | |
274 | kFALSE, // No fit performed | |
275 | AliAODTrack::kPrimary, | |
276 | 0); | |
277 | SetTrackChargeAndPID(pdg, aodTrack); | |
278 | fAODCTS->Add(aodTrack);//reference the selected object to the list | |
279 | } | |
280 | //Keep some charged particles in calorimeters lists | |
281 | if((fFillPHOS || fFillEMCAL) && KeepChargedParticles(pdg)) FillCalorimeters(iParticle, particle, momentum, ncalo); | |
1c5acb87 | 282 | |
477d6cee | 283 | }//Charged |
284 | ||
285 | //-------------Neutral particles ---------------------- | |
286 | else if(charge == 0 && (fFillPHOS || fFillEMCAL)){ | |
287 | //Skip neutrinos or other neutral particles | |
288 | //if(SkipNeutralParticles(pdg) || particle->IsPrimary()) continue ; // protection added (MG) | |
289 | if(SkipNeutralParticles(pdg)) continue ; | |
290 | //Fill particle/calocluster arrays | |
291 | if(!fDecayPi0) { | |
292 | FillCalorimeters(iParticle, particle, momentum, ncalo); | |
1c5acb87 | 293 | } |
477d6cee | 294 | else { |
295 | //Sometimes pi0 are stable for the generator, if needed decay it by hand | |
296 | if(pdg == 111 ){ | |
297 | if(momentum.Pt() > fPHOSPtMin || momentum.Pt() > fEMCALPtMin){ | |
298 | TLorentzVector lvGamma1, lvGamma2 ; | |
299 | //Double_t angle = 0; | |
300 | ||
301 | //Decay | |
302 | MakePi0Decay(momentum,lvGamma1,lvGamma2);//,angle); | |
303 | ||
304 | //Check if Pi0 is in the acceptance of the calorimeters, if aperture angle is small, keep it | |
305 | TParticle * pPhoton1 = new TParticle(22,1,iParticle,0,0,0,lvGamma1.Px(),lvGamma1.Py(), | |
306 | lvGamma1.Pz(),lvGamma1.E(),0,0,0,0); | |
307 | TParticle * pPhoton2 = new TParticle(22,1,iParticle,0,0,0,lvGamma2.Px(),lvGamma2.Py(), | |
308 | lvGamma2.Pz(),lvGamma2.E(),0,0,0,0); | |
309 | //Fill particle/calocluster arrays | |
310 | FillCalorimeters(iParticle,pPhoton1,lvGamma1, ncalo); | |
311 | FillCalorimeters(iParticle,pPhoton2,lvGamma2, ncalo); | |
312 | }//pt cut | |
313 | }//pi0 | |
314 | else FillCalorimeters(iParticle,particle, momentum, ncalo); //Add the rest | |
1c5acb87 | 315 | } |
477d6cee | 316 | }//neutral particles |
317 | } //particle with correct status | |
318 | }//particle loop | |
319 | ||
320 | fIndex2ndPhoton = -1; //In case of overlapping studies, reset for each event | |
1c5acb87 | 321 | } |
322 | ||
323 | //________________________________________________________________ | |
324 | Bool_t AliCaloTrackMCReader::KeepParticleWithStatus(Int_t status) const { | |
325 | //Check if status is equal to one of the list | |
326 | //These particles will be used in analysis. | |
327 | if(!fKeepAllStatus){ | |
328 | for(Int_t i= 0; i < fStatusArray->GetSize(); i++) | |
329 | if(status == fStatusArray->At(i)) return kTRUE ; | |
330 | ||
331 | return kFALSE; | |
332 | ||
333 | } | |
334 | else | |
335 | return kTRUE ; | |
336 | } | |
337 | ||
338 | //________________________________________________________________ | |
339 | Bool_t AliCaloTrackMCReader::KeepChargedParticles(Int_t pdg) const { | |
340 | //Check if pdg is equal to one of the charged particles list | |
341 | //These particles will be added to the calorimeters lists. | |
342 | ||
343 | for(Int_t i= 0; i < fChargedParticlesArray->GetSize(); i++) | |
344 | if(TMath::Abs(pdg) == fChargedParticlesArray->At(i)) return kTRUE ; | |
345 | ||
346 | return kFALSE; | |
347 | ||
348 | } | |
349 | ||
350 | //________________________________________________________________ | |
351 | void AliCaloTrackMCReader::Print(const Option_t * opt) const | |
352 | { | |
353 | //Print some relevant parameters set for the analysis | |
354 | if(! opt) | |
355 | return; | |
356 | ||
a3aebfff | 357 | printf("**** Print **** %s %s ****\n", GetName(), GetTitle() ) ; |
1c5acb87 | 358 | |
359 | printf("Decay Pi0? : %d\n", fDecayPi0) ; | |
691bdd02 | 360 | printf("Check Overlap in Calo? : %d\n", fCheckOverlap) ; |
1c5acb87 | 361 | printf("Keep all status? : %d\n", fKeepAllStatus) ; |
362 | ||
363 | if(!fKeepAllStatus) printf("Keep particles with status : "); | |
364 | for(Int_t i= 0; i < fStatusArray->GetSize(); i++) | |
365 | printf(" %d ; ", fStatusArray->At(i)); | |
366 | printf("\n"); | |
367 | ||
368 | printf("Skip neutral particles in calo : "); | |
369 | for(Int_t i= 0; i < fNeutralParticlesArray->GetSize(); i++) | |
370 | printf(" %d ; ", fNeutralParticlesArray->At(i)); | |
371 | printf("\n"); | |
372 | ||
373 | printf("Keep charged particles in calo : "); | |
374 | for(Int_t i= 0; i < fChargedParticlesArray->GetSize(); i++) | |
375 | printf(" %d ; ", fChargedParticlesArray->At(i)); | |
376 | printf("\n"); | |
377 | ||
378 | } | |
379 | ||
380 | //____________________________________________________________________________ | |
381 | void AliCaloTrackMCReader::MakePi0Decay(TLorentzVector &p0, TLorentzVector &p1, | |
382 | TLorentzVector &p2) const {//, Double_t &angle) { | |
383 | // Perform isotropic decay pi0 -> 2 photons | |
384 | // p0 is pi0 4-momentum (inut) | |
385 | // p1 and p2 are photon 4-momenta (output) | |
386 | // cout<<"Boost vector"<<endl; | |
387 | Double_t mPi0 = TDatabasePDG::Instance()->GetParticle(111)->Mass(); | |
388 | TVector3 b = p0.BoostVector(); | |
389 | //cout<<"Parameters"<<endl; | |
390 | //Double_t mPi0 = p0.M(); | |
391 | Double_t phi = TMath::TwoPi() * gRandom->Rndm(); | |
392 | Double_t cosThe = 2 * gRandom->Rndm() - 1; | |
393 | Double_t cosPhi = TMath::Cos(phi); | |
394 | Double_t sinPhi = TMath::Sin(phi); | |
395 | Double_t sinThe = TMath::Sqrt(1-cosThe*cosThe); | |
396 | Double_t ePi0 = mPi0/2.; | |
397 | //cout<<"ePi0 "<<ePi0<<endl; | |
398 | //cout<<"Components"<<endl; | |
399 | p1.SetPx(+ePi0*cosPhi*sinThe); | |
400 | p1.SetPy(+ePi0*sinPhi*sinThe); | |
401 | p1.SetPz(+ePi0*cosThe); | |
402 | p1.SetE(ePi0); | |
403 | //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl; | |
404 | //cout<<"p1 Mass: "<<p1.Px()*p1.Px()+p1.Py()*p1.Py()+p1.Pz()*p1.Pz()-p1.E()*p1.E()<<endl; | |
405 | p2.SetPx(-ePi0*cosPhi*sinThe); | |
406 | p2.SetPy(-ePi0*sinPhi*sinThe); | |
407 | p2.SetPz(-ePi0*cosThe); | |
408 | p2.SetE(ePi0); | |
409 | //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl; | |
410 | //cout<<"p2 Mass: "<<p2.Px()*p2.Px()+p2.Py()*p2.Py()+p2.Pz()*p2.Pz()-p2.E()*p2.E()<<endl; | |
411 | //cout<<"Boost "<<b.X()<<" "<<b.Y()<<" "<<b.Z()<<endl; | |
412 | p1.Boost(b); | |
413 | //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl; | |
414 | p2.Boost(b); | |
415 | //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl; | |
416 | //cout<<"angle"<<endl; | |
417 | //angle = p1.Angle(p2.Vect()); | |
418 | //cout<<angle<<endl; | |
419 | } | |
420 | ||
421 | //____________________________________________________________________________ | |
477d6cee | 422 | void AliCaloTrackMCReader::SetInputOutputMCEvent(AliVEvent* /*esd*/, AliAODEvent* aod, AliMCEvent* mc) { |
1c5acb87 | 423 | // Connect the data pointer |
477d6cee | 424 | SetMC(mc); |
425 | SetOutputEvent(aod); | |
1c5acb87 | 426 | } |
427 | ||
428 | ||
429 | //________________________________________________________________ | |
430 | Bool_t AliCaloTrackMCReader::SkipNeutralParticles(Int_t pdg) const { | |
431 | //Check if pdg is equal to one of the neutral particles list | |
432 | //These particles will be skipped from analysis. | |
433 | ||
434 | for(Int_t i= 0; i < fNeutralParticlesArray->GetSize(); i++) | |
435 | if(TMath::Abs(pdg) == fNeutralParticlesArray->At(i)) return kTRUE ; | |
436 | ||
437 | return kFALSE; | |
438 | ||
439 | } | |
440 | ||
441 | ||
442 | //____________________________________________________________________ | |
443 | void AliCaloTrackMCReader::SetTrackChargeAndPID(const Int_t pdgCode, AliAODTrack *track) const { | |
444 | //Give a PID weight for tracks equal to 1 depending on the particle type | |
445 | ||
446 | Float_t pid[10] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; | |
447 | ||
448 | switch (pdgCode) { | |
449 | ||
450 | case 22: // gamma | |
451 | track->SetCharge(0); | |
452 | pid[AliAODTrack::kUnknown] = 1.; | |
453 | track->SetPID(pid); | |
454 | break; | |
455 | ||
456 | case 11: // e- | |
457 | track->SetCharge(-1); | |
458 | pid[AliAODTrack::kElectron] = 1.; | |
459 | track->SetPID(pid); | |
460 | break; | |
461 | ||
462 | case -11: // e+ | |
463 | track->SetCharge(+1); | |
464 | pid[AliAODTrack::kElectron] = 1.; | |
465 | track->SetPID(pid); | |
466 | break; | |
467 | ||
468 | case 13: // mu- | |
469 | track->SetCharge(-1); | |
470 | pid[AliAODTrack::kMuon] = 1.; | |
471 | track->SetPID(pid); | |
472 | break; | |
473 | ||
474 | case -13: // mu+ | |
475 | track->SetCharge(+1); | |
476 | pid[AliAODTrack::kMuon] = 1.; | |
477 | track->SetPID(pid); | |
478 | break; | |
479 | ||
480 | case 111: // pi0 | |
481 | track->SetCharge(0); | |
482 | pid[AliAODTrack::kUnknown] = 1.; | |
483 | track->SetPID(pid); | |
484 | break; | |
485 | ||
486 | case 211: // pi+ | |
487 | track->SetCharge(+1); | |
488 | pid[AliAODTrack::kPion] = 1.; | |
489 | track->SetPID(pid); | |
490 | break; | |
491 | ||
492 | case -211: // pi- | |
493 | track->SetCharge(-1); | |
494 | pid[AliAODTrack::kPion] = 1.; | |
495 | track->SetPID(pid); | |
496 | break; | |
497 | ||
498 | case 130: // K0L | |
499 | track->SetCharge(0); | |
500 | pid[AliAODTrack::kUnknown] = 1.; | |
501 | track->SetPID(pid); | |
502 | break; | |
503 | ||
504 | case 321: // K+ | |
505 | track->SetCharge(+1); | |
506 | pid[AliAODTrack::kKaon] = 1.; | |
507 | track->SetPID(pid); | |
508 | break; | |
509 | ||
510 | case -321: // K- | |
511 | track->SetCharge(-1); | |
512 | pid[AliAODTrack::kKaon] = 1.; | |
513 | track->SetPID(pid); | |
514 | break; | |
515 | ||
516 | case 2112: // n | |
517 | track->SetCharge(0); | |
518 | pid[AliAODTrack::kUnknown] = 1.; | |
519 | track->SetPID(pid); | |
520 | break; | |
521 | ||
522 | case 2212: // p | |
523 | track->SetCharge(+1); | |
524 | pid[AliAODTrack::kProton] = 1.; | |
525 | track->SetPID(pid); | |
526 | break; | |
527 | ||
528 | case -2212: // anti-p | |
529 | track->SetCharge(-1); | |
530 | pid[AliAODTrack::kProton] = 1.; | |
531 | track->SetPID(pid); | |
532 | break; | |
533 | ||
534 | case 310: // K0S | |
535 | track->SetCharge(0); | |
536 | pid[AliAODTrack::kUnknown] = 1.; | |
537 | track->SetPID(pid); | |
538 | break; | |
539 | ||
540 | case 311: // K0 | |
541 | track->SetCharge(0); | |
542 | pid[AliAODTrack::kUnknown] = 1.; | |
543 | track->SetPID(pid); | |
544 | break; | |
545 | ||
546 | case -311: // anti-K0 | |
547 | track->SetCharge(0); | |
548 | pid[AliAODTrack::kUnknown] = 1.; | |
549 | track->SetPID(pid); | |
550 | break; | |
551 | ||
552 | case 221: // eta | |
553 | track->SetCharge(0); | |
554 | pid[AliAODTrack::kUnknown] = 1.; | |
555 | track->SetPID(pid); | |
556 | break; | |
557 | ||
558 | case 3122: // lambda | |
559 | track->SetCharge(0); | |
560 | pid[AliAODTrack::kUnknown] = 1.; | |
561 | track->SetPID(pid); | |
562 | break; | |
563 | ||
564 | case 3222: // Sigma+ | |
565 | track->SetCharge(+1); | |
566 | pid[AliAODTrack::kUnknown] = 1.; | |
567 | track->SetPID(pid); | |
568 | break; | |
569 | ||
570 | case 3212: // Sigma0 | |
571 | track->SetCharge(-1); | |
572 | pid[AliAODTrack::kUnknown] = 1.; | |
573 | track->SetPID(pid); | |
574 | break; | |
575 | ||
576 | case 3112: // Sigma- | |
577 | track->SetCharge(-1); | |
578 | pid[AliAODTrack::kUnknown] = 1.; | |
579 | track->SetPID(pid); | |
580 | break; | |
581 | ||
582 | case 3322: // Xi0 | |
583 | track->SetCharge(0); | |
584 | pid[AliAODTrack::kUnknown] = 1.; | |
585 | track->SetPID(pid); | |
586 | break; | |
587 | ||
588 | case 3312: // Xi- | |
589 | track->SetCharge(-1); | |
590 | pid[AliAODTrack::kUnknown] = 1.; | |
591 | track->SetPID(pid); | |
592 | break; | |
593 | ||
594 | case 3334: // Omega- | |
595 | track->SetCharge(-1); | |
596 | pid[AliAODTrack::kUnknown] = 1.; | |
597 | track->SetPID(pid); | |
598 | break; | |
599 | ||
600 | case -2112: // n-bar | |
601 | track->SetCharge(0); | |
602 | pid[AliAODTrack::kUnknown] = 1.; | |
603 | track->SetPID(pid); | |
604 | break; | |
605 | ||
606 | case -3122: // anti-Lambda | |
607 | track->SetCharge(0); | |
608 | pid[AliAODTrack::kUnknown] = 1.; | |
609 | track->SetPID(pid); | |
610 | break; | |
611 | ||
612 | case -3222: // anti-Sigma- | |
613 | track->SetCharge(-1); | |
614 | pid[AliAODTrack::kUnknown] = 1.; | |
615 | track->SetPID(pid); | |
616 | break; | |
617 | ||
618 | case -3212: // anti-Sigma0 | |
619 | track->SetCharge(0); | |
620 | pid[AliAODTrack::kUnknown] = 1.; | |
621 | track->SetPID(pid); | |
622 | break; | |
623 | ||
624 | case -3112: // anti-Sigma+ | |
625 | track->SetCharge(+1); | |
626 | pid[AliAODTrack::kUnknown] = 1.; | |
627 | track->SetPID(pid); | |
628 | break; | |
629 | ||
630 | case -3322: // anti-Xi0 | |
631 | track->SetCharge(0); | |
632 | pid[AliAODTrack::kUnknown] = 1.; | |
633 | track->SetPID(pid); | |
634 | break; | |
635 | ||
636 | case -3312: // anti-Xi+ | |
637 | track->SetCharge(+1); | |
638 | break; | |
639 | ||
640 | case -3334: // anti-Omega+ | |
641 | track->SetCharge(+1); | |
642 | pid[AliAODTrack::kUnknown] = 1.; | |
643 | track->SetPID(pid); | |
644 | break; | |
645 | ||
646 | case 411: // D+ | |
647 | track->SetCharge(+1); | |
648 | pid[AliAODTrack::kUnknown] = 1.; | |
649 | track->SetPID(pid); | |
650 | break; | |
651 | ||
652 | case -411: // D- | |
653 | track->SetCharge(-1); | |
654 | pid[AliAODTrack::kUnknown] = 1.; | |
655 | track->SetPID(pid); | |
656 | break; | |
657 | ||
658 | case 421: // D0 | |
659 | track->SetCharge(0); | |
660 | pid[AliAODTrack::kUnknown] = 1.; | |
661 | track->SetPID(pid); | |
662 | break; | |
663 | ||
664 | case -421: // anti-D0 | |
665 | track->SetCharge(0); | |
666 | pid[AliAODTrack::kUnknown] = 1.; | |
667 | track->SetPID(pid); | |
668 | break; | |
669 | ||
670 | default : // unknown | |
671 | track->SetCharge(-99); | |
672 | pid[AliAODTrack::kUnknown] = 1.; | |
673 | track->SetPID(pid); | |
674 | } | |
675 | ||
676 | track->SetPID(pid); | |
677 | ||
678 | return; | |
679 | } | |
680 | ||
681 | //____________________________________________________________________ | |
682 | void AliCaloTrackMCReader::SetCaloClusterPID(const Int_t pdgCode, AliAODCaloCluster *calo) const { | |
683 | //Give a PID weight for CaloClusters equal to 1 depending on the particle type | |
684 | ||
685 | Float_t pid[13] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; | |
686 | ||
687 | switch (pdgCode) { | |
688 | ||
689 | case 22: // gamma | |
690 | pid[AliAODCaloCluster::kPhoton] = 1.; | |
691 | calo->SetPID(pid); | |
692 | break; | |
693 | ||
694 | case 11: // e- | |
695 | pid[AliAODCaloCluster::kElectron] = 1.; | |
696 | calo->SetPID(pid); | |
697 | break; | |
698 | ||
699 | case -11: // e+ | |
700 | pid[AliAODCaloCluster::kElectron] = 1.; | |
701 | calo->SetPID(pid); | |
702 | break; | |
703 | ||
704 | case 13: // mu- | |
705 | pid[AliAODCaloCluster::kCharged] = 1.; | |
706 | calo->SetPID(pid); | |
707 | break; | |
708 | ||
709 | case -13: // mu+ | |
710 | pid[AliAODCaloCluster::kCharged] = 1.; | |
711 | calo->SetPID(pid); | |
712 | break; | |
713 | ||
714 | case 111: // pi0 | |
715 | pid[AliAODCaloCluster::kPi0] = 1.; | |
716 | calo->SetPID(pid); | |
717 | break; | |
718 | ||
719 | case 211: // pi+ | |
720 | pid[AliAODCaloCluster::kCharged] = 1.; | |
721 | calo->SetPID(pid); | |
722 | break; | |
723 | ||
724 | case -211: // pi- | |
725 | pid[AliAODCaloCluster::kCharged] = 1.; | |
726 | calo->SetPID(pid); | |
727 | break; | |
728 | ||
729 | case 130: // K0L | |
730 | pid[AliAODCaloCluster::kKaon0] = 1.; | |
731 | pid[AliAODCaloCluster::kNeutral] = 1; | |
732 | calo->SetPID(pid); | |
733 | break; | |
734 | ||
735 | case 321: // K+ | |
736 | pid[AliAODCaloCluster::kCharged] = 1.; | |
737 | calo->SetPID(pid); | |
738 | break; | |
739 | ||
740 | case -321: // K- | |
741 | pid[AliAODCaloCluster::kCharged] = 1.; | |
742 | calo->SetPID(pid); | |
743 | break; | |
744 | ||
745 | case 2112: // n | |
746 | pid[AliAODCaloCluster::kNeutron] = 1.; | |
747 | pid[AliAODCaloCluster::kNeutral] = 1.; | |
748 | calo->SetPID(pid); | |
749 | break; | |
750 | ||
751 | case 2212: // p | |
752 | pid[AliAODCaloCluster::kCharged] = 1.; | |
753 | calo->SetPID(pid); | |
754 | break; | |
755 | ||
756 | case -2212: // anti-p | |
757 | pid[AliAODCaloCluster::kCharged] = 1.; | |
758 | calo->SetPID(pid); | |
759 | break; | |
760 | ||
761 | case 310: // K0S | |
762 | pid[AliAODCaloCluster::kKaon0] = 1.; | |
763 | pid[AliAODCaloCluster::kNeutral] = 1.; | |
764 | calo->SetPID(pid); | |
765 | break; | |
766 | ||
767 | case 311: // K0 | |
768 | pid[AliAODCaloCluster::kKaon0] = 1.; | |
769 | pid[AliAODCaloCluster::kNeutral] = 1.; | |
770 | calo->SetPID(pid); | |
771 | break; | |
772 | ||
773 | case -311: // anti-K0 | |
774 | pid[AliAODCaloCluster::kKaon0] = 1.; | |
775 | pid[AliAODCaloCluster::kNeutral] = 1.; | |
776 | calo->SetPID(pid); | |
777 | break; | |
778 | ||
779 | case 221: // eta | |
780 | pid[AliAODCaloCluster::kNeutral] = 1.; | |
781 | calo->SetPID(pid); | |
782 | break; | |
783 | ||
784 | case 3122: // lambda | |
785 | pid[AliAODCaloCluster::kUnknown] = 1.; | |
786 | calo->SetPID(pid); | |
787 | break; | |
788 | ||
789 | case 3222: // Sigma+ | |
790 | pid[AliAODCaloCluster::kUnknown] = 1.; | |
791 | calo->SetPID(pid); | |
792 | break; | |
793 | ||
794 | case 3212: // Sigma0 | |
795 | pid[AliAODCaloCluster::kUnknown] = 1.; | |
796 | calo->SetPID(pid); | |
797 | break; | |
798 | ||
799 | case 3112: // Sigma- | |
800 | pid[AliAODCaloCluster::kUnknown] = 1.; | |
801 | calo->SetPID(pid); | |
802 | break; | |
803 | ||
804 | case 3322: // Xi0 | |
805 | pid[AliAODCaloCluster::kUnknown] = 1.; | |
806 | calo->SetPID(pid); | |
807 | break; | |
808 | ||
809 | case 3312: // Xi- | |
810 | pid[AliAODCaloCluster::kUnknown] = 1.; | |
811 | calo->SetPID(pid); | |
812 | break; | |
813 | ||
814 | case 3334: // Omega- | |
815 | pid[AliAODCaloCluster::kUnknown] = 1.; | |
816 | calo->SetPID(pid); | |
817 | break; | |
818 | ||
819 | case -2112: // n-bar | |
820 | pid[AliAODCaloCluster::kNeutron] = 1.; | |
821 | pid[AliAODCaloCluster::kNeutral] = 1.; | |
822 | calo->SetPID(pid); | |
823 | break; | |
824 | ||
825 | case -3122: // anti-Lambda | |
826 | pid[AliAODCaloCluster::kUnknown] = 1.; | |
827 | calo->SetPID(pid); | |
828 | break; | |
829 | ||
830 | case -3222: // anti-Sigma- | |
831 | pid[AliAODCaloCluster::kUnknown] = 1.; | |
832 | calo->SetPID(pid); | |
833 | break; | |
834 | ||
835 | case -3212: // anti-Sigma0 | |
836 | pid[AliAODCaloCluster::kUnknown] = 1.; | |
837 | calo->SetPID(pid); | |
838 | break; | |
839 | ||
840 | case -3112: // anti-Sigma+ | |
841 | pid[AliAODCaloCluster::kUnknown] = 1.; | |
842 | calo->SetPID(pid); | |
843 | break; | |
844 | ||
845 | case -3322: // anti-Xi0 | |
846 | pid[AliAODCaloCluster::kUnknown] = 1.; | |
847 | calo->SetPID(pid); | |
848 | break; | |
849 | ||
850 | case -3312: // anti-Xi+ | |
851 | pid[AliAODCaloCluster::kUnknown] = 1.; | |
852 | calo->SetPID(pid); | |
853 | break; | |
854 | ||
855 | case -3334: // anti-Omega+ | |
856 | pid[AliAODCaloCluster::kUnknown] = 1.; | |
857 | calo->SetPID(pid); | |
858 | break; | |
859 | ||
860 | case 411: // D+ | |
861 | pid[AliAODCaloCluster::kUnknown] = 1.; | |
862 | calo->SetPID(pid); | |
863 | break; | |
864 | ||
865 | case -411: // D- | |
866 | pid[AliAODCaloCluster::kUnknown] = 1.; | |
867 | calo->SetPID(pid); | |
868 | break; | |
869 | ||
870 | case 421: // D0 | |
871 | pid[AliAODCaloCluster::kUnknown] = 1.; | |
872 | calo->SetPID(pid); | |
873 | break; | |
874 | ||
875 | case -421: // anti-D0 | |
876 | pid[AliAODCaloCluster::kUnknown] = 1.; | |
877 | calo->SetPID(pid); | |
878 | break; | |
879 | ||
880 | default : // unknown | |
881 | pid[AliAODCaloCluster::kUnknown] = 1.; | |
882 | calo->SetPID(pid); | |
883 | } | |
884 | ||
885 | ||
886 | return; | |
887 | } |