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