]>
Commit | Line | Data |
---|---|---|
bdcfac30 | 1 | |
2 | /************************************************************************** | |
3 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
4 | * * | |
5 | * Author: The ALICE Off-line Project. * | |
6 | * Contributors are mentioned in the code where appropriate. * | |
7 | * * | |
8 | * Permission to use, copy, modify and distribute this software and its * | |
9 | * documentation strictly for non-commercial purposes is hereby granted * | |
10 | * without fee, provided that the above copyright notice appears in all * | |
11 | * copies and that both the copyright notice and this permission notice * | |
12 | * appear in the supporting documentation. The authors make no claims * | |
13 | * about the suitability of this software for any purpose. It is * | |
14 | * provided "as is" without express or implied warranty. * | |
15 | **************************************************************************/ | |
16 | /* $Id$ */ | |
17 | ||
18 | /* History of cvs commits: | |
19 | * | |
20 | * $Log$ | |
307cf3c6 | 21 | * Revision 1.5 2007/12/14 10:54:15 gustavo |
22 | * Protection against cluster labels larger than kinematics entries | |
23 | * | |
b398a6cb | 24 | * Revision 1.4 2007/10/29 13:48:42 gustavo |
25 | * Corrected coding violations | |
26 | * | |
4b707925 | 27 | * Revision 1.2 2007/08/17 12:40:04 schutz |
28 | * New analysis classes by Gustavo Conesa | |
29 | * | |
bdcfac30 | 30 | * Revision 1.1.2.1 2007/07/26 10:32:09 schutz |
31 | * new analysis classes in the the new analysis framework | |
32 | * | |
33 | * | |
34 | */ | |
35 | ||
36 | //_________________________________________________________________________ | |
37 | // Class for reading data (Kinematics and ESDs) in order to do prompt gamma correlations | |
38 | // Class created from old AliPHOSGammaJet | |
39 | // (see AliRoot versions previous Release 4-09) | |
40 | // | |
41 | //*-- Author: Gustavo Conesa (LNF-INFN) | |
42 | ////////////////////////////////////////////////////////////////////////////// | |
43 | ||
44 | ||
45 | // --- ROOT system --- | |
3bb2c538 | 46 | #include <TFormula.h> |
bdcfac30 | 47 | #include <TParticle.h> |
3bb2c538 | 48 | |
bdcfac30 | 49 | //---- ANALYSIS system ---- |
50 | #include "AliGammaMCDataReader.h" | |
bdcfac30 | 51 | #include "AliESDEvent.h" |
52 | #include "AliESDVertex.h" | |
53 | #include "AliESDCaloCluster.h" | |
54 | #include "AliStack.h" | |
3bb2c538 | 55 | #include "AliLog.h" |
bdcfac30 | 56 | |
57 | ClassImp(AliGammaMCDataReader) | |
58 | ||
59 | //____________________________________________________________________________ | |
60 | AliGammaMCDataReader::AliGammaMCDataReader() : | |
4b707925 | 61 | AliGammaReader() |
bdcfac30 | 62 | { |
63 | //Default Ctor | |
64 | ||
65 | //Initialize parameters | |
66 | fDataType=kMCData; | |
bdcfac30 | 67 | |
68 | } | |
69 | ||
70 | //____________________________________________________________________________ | |
71 | AliGammaMCDataReader::AliGammaMCDataReader(const AliGammaMCDataReader & g) : | |
4b707925 | 72 | AliGammaReader(g) |
bdcfac30 | 73 | { |
74 | // cpy ctor | |
75 | } | |
76 | ||
77 | //_________________________________________________________________________ | |
78 | AliGammaMCDataReader & AliGammaMCDataReader::operator = (const AliGammaMCDataReader & source) | |
79 | { | |
80 | // assignment operator | |
81 | ||
82 | if(&source == this) return *this; | |
83 | ||
4b707925 | 84 | |
bdcfac30 | 85 | return *this; |
86 | ||
87 | } | |
88 | ||
89 | //____________________________________________________________________________ | |
90 | void AliGammaMCDataReader::CreateParticleList(TObject * data, TObject * kine, | |
91 | TClonesArray * plCTS, | |
92 | TClonesArray * plEMCAL, | |
4b707925 | 93 | TClonesArray * plPHOS, |
94 | TClonesArray * plPrimCTS, | |
95 | TClonesArray * plPrimEMCAL, | |
96 | TClonesArray * plPrimPHOS){ | |
bdcfac30 | 97 | |
98 | //Create a list of particles from the ESD. These particles have been measured | |
3bb2c538 | 99 | //by the Central Tracking system (TPC+ITS+...), PHOS and EMCAL |
100 | //Also create particle list with mothers. | |
bdcfac30 | 101 | |
102 | AliESDEvent* esd = (AliESDEvent*) data; | |
103 | AliStack* stack = (AliStack*) kine; | |
104 | ||
105 | Int_t npar = 0 ; | |
4b707925 | 106 | Double_t *pid = new Double_t[AliPID::kSPECIESN]; |
bdcfac30 | 107 | AliDebug(3,"Fill particle lists"); |
108 | ||
109 | //Get vertex for momentum calculation | |
110 | Double_t v[3] ; //vertex ; | |
111 | esd->GetVertex()->GetXYZ(v) ; | |
112 | ||
3bb2c538 | 113 | //########### CALORIMETERS ############## |
114 | Int_t nCaloCluster = esd->GetNumberOfCaloClusters() ; | |
bdcfac30 | 115 | Int_t indexPH = plPHOS->GetEntries() ; |
3bb2c538 | 116 | Int_t indexEM = plEMCAL->GetEntries() ; |
bdcfac30 | 117 | |
3bb2c538 | 118 | for (npar = 0; npar < nCaloCluster; npar++) {//////////////CaloCluster loop |
119 | AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve cluster from esd | |
120 | Int_t type = clus->GetClusterType(); | |
121 | ||
122 | //########### PHOS ############## | |
123 | if(fSwitchOnPHOS && type == AliESDCaloCluster::kPHOSCluster){ | |
124 | AliDebug(4,Form("PHOS clusters: E %f, match %d", clus->E(),clus->GetTrackMatched())); | |
125 | ||
126 | if(clus->GetTrackMatched()==-1){ | |
127 | TLorentzVector momentum ; | |
128 | clus->GetMomentum(momentum, v); | |
129 | Double_t phi = momentum.Phi(); | |
130 | if(phi<0) phi+=TMath::TwoPi() ; | |
131 | if(momentum.Pt() > fNeutralPtCut && TMath::Abs(momentum.Eta()) < fPHOSEtaCut && | |
132 | phi > fPhiPHOSCut[0] && phi < fPhiPHOSCut[1] ) { | |
4b707925 | 133 | |
3bb2c538 | 134 | pid=clus->GetPid(); |
135 | Int_t pdg = 22; | |
4b707925 | 136 | |
3bb2c538 | 137 | if(IsPHOSPIDOn()){ |
138 | AliDebug(5,Form("E %1.2f; PID: ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f,pi %0.2f, k %0.2f, p %0.2f, k0 %0.2f, n %0.2f, mu %0.2f ", | |
139 | momentum.E(),pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion], | |
140 | pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon])); | |
141 | ||
142 | Float_t wPhoton = fPHOSPhotonWeight; | |
143 | Float_t wPi0 = fPHOSPi0Weight; | |
144 | ||
145 | if(fPHOSWeightFormula){ | |
146 | wPhoton = fPHOSPhotonWeightFormula->Eval(momentum.E()) ; | |
147 | wPi0 = fPHOSPi0WeightFormula->Eval(momentum.E()); | |
148 | } | |
149 | ||
150 | if(pid[AliPID::kPhoton] > wPhoton) | |
151 | pdg = kPhoton ; | |
152 | else if(pid[AliPID::kPi0] > wPi0) | |
153 | pdg = kPi0 ; | |
154 | else if(pid[AliPID::kElectron] > fPHOSElectronWeight) | |
155 | pdg = kElectron ; | |
156 | else if(pid[AliPID::kEleCon] > fPHOSElectronWeight) | |
157 | pdg = kEleCon ; | |
158 | else if(pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] > fPHOSChargeWeight) | |
159 | pdg = kChargedHadron ; | |
160 | else if(pid[AliPID::kKaon0]+pid[AliPID::kNeutron] > fPHOSNeutralWeight) | |
161 | pdg = kNeutralHadron ; | |
162 | ||
163 | else if(pid[AliPID::kElectron]+pid[AliPID::kEleCon]+pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] > | |
164 | pid[AliPID::kPhoton] + pid[AliPID::kPi0]+pid[AliPID::kKaon0]+pid[AliPID::kNeutron]) | |
165 | pdg = kChargedUnknown ; | |
166 | else | |
167 | pdg = kNeutralUnknown ; | |
168 | //neutral cluster, unidentifed. | |
4b707925 | 169 | } |
170 | ||
3bb2c538 | 171 | if(pdg != kElectron && pdg != kEleCon && pdg !=kChargedHadron && pdg !=kChargedUnknown ){//keep only neutral particles in the array |
b398a6cb | 172 | |
173 | TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0); | |
174 | new((*plPHOS)[indexPH]) TParticle(*particle) ; | |
175 | AliDebug(4,Form("PHOS added: pdg %d, pt %f, phi %f, eta %f", pdg, particle->Pt(),particle->Phi(),particle->Eta())); | |
176 | ||
3bb2c538 | 177 | //############### |
178 | //Check kinematics | |
179 | //############### | |
b398a6cb | 180 | TParticle * pmother = new TParticle(); |
3bb2c538 | 181 | Int_t label = clus->GetLabel(); |
307cf3c6 | 182 | if(label < stack->GetNtrack()) |
b398a6cb | 183 | pmother = GetMotherParticle(label,stack, "PHOS",momentum); |
184 | else | |
307cf3c6 | 185 | AliInfo(Form("PHOS: Bad label %d, too large, NPrimaries %d",label,stack->GetNtrack())); |
3bb2c538 | 186 | new((*plPrimPHOS)[indexPH]) TParticle(*pmother) ; |
3bb2c538 | 187 | |
b398a6cb | 188 | indexPH++; |
3bb2c538 | 189 | } |
190 | else AliDebug(4,Form("PHOS charged cluster NOT added: pdg %d, pt %f, phi %f, eta %f\n", | |
191 | pdg, momentum.Pt(),momentum.Phi(),momentum.Eta())); | |
4b707925 | 192 | |
3bb2c538 | 193 | }//pt, eta, phi cut |
194 | else AliDebug(4,"Particle not added"); | |
195 | }//track-match? | |
196 | }//PHOS cluster | |
4b707925 | 197 | |
3bb2c538 | 198 | //################ EMCAL ############## |
199 | else if(fSwitchOnEMCAL && type == AliESDCaloCluster::kEMCALClusterv1){ | |
200 | AliDebug(4,Form("EMCAL clusters: E %f, match %d", clus->E(),clus->GetTrackMatched())); | |
201 | ||
202 | if(clus->GetTrackMatched()==-1 ){ | |
203 | TLorentzVector momentum ; | |
204 | clus->GetMomentum(momentum, v); | |
205 | Double_t phi = momentum.Phi(); | |
206 | if(phi<0) phi+=TMath::TwoPi() ; | |
207 | if(momentum.Pt() > fNeutralPtCut && TMath::Abs(momentum.Eta()) < fEMCALEtaCut && | |
208 | phi > fPhiEMCALCut[0] && phi < fPhiEMCALCut[1] ) { | |
4b707925 | 209 | |
3bb2c538 | 210 | pid=clus->GetPid(); |
211 | Int_t pdg = 22; | |
4b707925 | 212 | |
3bb2c538 | 213 | if(IsEMCALPIDOn()){ |
214 | AliDebug(5,Form("E %1.2f; PID: ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f,pi %0.2f, k %0.2f, p %0.2f, k0 %0.2f, n %0.2f, mu %0.2f ", | |
215 | momentum.E(),pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion], | |
216 | pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon])); | |
217 | ||
218 | if(pid[AliPID::kPhoton] > fEMCALPhotonWeight) | |
219 | pdg = kPhoton ; | |
220 | else if(pid[AliPID::kPi0] > fEMCALPi0Weight) | |
221 | pdg = kPi0 ; | |
222 | else if(pid[AliPID::kElectron] > fEMCALElectronWeight) | |
223 | pdg = kElectron ; | |
224 | else if(pid[AliPID::kEleCon] > fEMCALElectronWeight) | |
225 | pdg = kEleCon ; | |
226 | else if(pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] > fEMCALChargeWeight) | |
227 | pdg = kChargedHadron ; | |
228 | else if(pid[AliPID::kKaon0]+pid[AliPID::kNeutron] > fEMCALNeutralWeight) | |
229 | pdg = kNeutralHadron ; | |
230 | else if(pid[AliPID::kElectron]+pid[AliPID::kEleCon]+pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] > | |
231 | pid[AliPID::kPhoton] + pid[AliPID::kPi0]+pid[AliPID::kKaon0]+pid[AliPID::kNeutron]) | |
232 | pdg = kChargedUnknown ; | |
233 | else | |
234 | pdg = kNeutralUnknown ; | |
235 | } | |
236 | ||
237 | if(pdg != kElectron && pdg != kEleCon && pdg !=kChargedHadron && pdg !=kChargedUnknown){//keep only neutral particles in the array | |
b398a6cb | 238 | |
239 | TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0); | |
240 | new((*plEMCAL)[indexEM]) TParticle(*particle) ; | |
241 | AliDebug(4,Form("EMCAL cluster added: pdg %f, pt %f, phi %f, eta %f", pdg, particle->Pt(),particle->Phi(),particle->Eta())); | |
3bb2c538 | 242 | |
243 | //############### | |
244 | //Check kinematics | |
245 | //############### | |
b398a6cb | 246 | TParticle * pmother = new TParticle(); |
3bb2c538 | 247 | Int_t label = clus->GetLabel(); |
307cf3c6 | 248 | if(label < stack->GetNtrack()) |
b398a6cb | 249 | pmother = GetMotherParticle(label,stack, "EMCAL",momentum); |
250 | else | |
307cf3c6 | 251 | AliInfo(Form("EMCAL: Bad label %d, too large, NPrimaries %d",label,stack->GetNtrack())); |
3bb2c538 | 252 | new((*plPrimEMCAL)[indexEM]) TParticle(*pmother) ; |
b398a6cb | 253 | |
3bb2c538 | 254 | indexEM++; |
255 | } | |
256 | else AliDebug(4,Form("EMCAL charged cluster NOT added: pdg %d, pt %f, phi %f, eta %f", | |
257 | pdg, momentum.Pt(),momentum.Phi(),momentum.Eta())); | |
258 | ||
259 | }//pt, phi, eta cut | |
260 | else AliDebug(4,"Particle not added"); | |
261 | }//track-matched | |
262 | }//EMCAL cluster | |
4b707925 | 263 | |
bdcfac30 | 264 | }//cluster loop |
265 | ||
3bb2c538 | 266 | |
267 | ||
bdcfac30 | 268 | //########### CTS (TPC+ITS) ##################### |
3bb2c538 | 269 | Int_t nTracks = esd->GetNumberOfTracks() ; |
bdcfac30 | 270 | Int_t indexCh = plCTS->GetEntries() ; |
4b707925 | 271 | |
3bb2c538 | 272 | if(fSwitchOnCTS){ |
273 | AliDebug(3,Form("Number of tracks %d",nTracks)); | |
bdcfac30 | 274 | |
3bb2c538 | 275 | for (npar = 0; npar < nTracks; npar++) {////////////// track loop |
276 | AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd | |
bdcfac30 | 277 | |
3bb2c538 | 278 | //We want tracks fitted in the detectors: |
279 | ULong_t status=AliESDtrack::kTPCrefit; | |
280 | status|=AliESDtrack::kITSrefit; | |
281 | ||
282 | //We want tracks whose PID bit is set: | |
283 | // ULong_t status =AliESDtrack::kITSpid; | |
284 | // status|=AliESDtrack::kTPCpid; | |
bdcfac30 | 285 | |
3bb2c538 | 286 | if ( (track->GetStatus() & status) == status) {//Check if the bits we want are set |
287 | // Do something with the tracks which were successfully | |
288 | // re-fitted | |
289 | Double_t en = 0; //track ->GetTPCsignal() ; | |
290 | Double_t mom[3]; | |
291 | track->GetPxPyPz(mom) ; | |
292 | Double_t px = mom[0]; | |
293 | Double_t py = mom[1]; | |
294 | Double_t pz = mom[2]; //Check with TPC people if this is correct. | |
295 | Int_t pdg = 11; //Give any charged PDG code, in this case electron. | |
296 | //I just want to tag the particle as charged | |
4b707925 | 297 | |
4b707925 | 298 | //Check kinematics |
3bb2c538 | 299 | Int_t label = TMath::Abs(track->GetLabel()); |
300 | TParticle * pmother = new TParticle(); | |
301 | pmother = stack->Particle(label); | |
4b707925 | 302 | |
303 | TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, | |
3bb2c538 | 304 | px, py, pz, en, v[0], v[1], v[2], 0); |
4b707925 | 305 | |
3bb2c538 | 306 | if(particle->Pt() > fChargedPtCut && TMath::Abs(particle->Eta())<fCTSEtaCut){ |
307 | new((*plCTS)[indexCh]) TParticle(*particle) ; | |
308 | new((*plPrimCTS)[indexCh]) TParticle(*pmother) ; | |
309 | indexCh++; | |
310 | ||
311 | }//kinematic selection | |
312 | }//select track from refit | |
313 | }//track loop | |
314 | }//CTS | |
bdcfac30 | 315 | |
3bb2c538 | 316 | AliDebug(3,Form("Particle lists filled, tracks %d , clusters: EMCAL %d, PHOS %d", indexCh,indexEM,indexPH)); |
317 | ||
318 | } | |
4b707925 | 319 | |
320 | TParticle * AliGammaMCDataReader:: GetMotherParticle(Int_t label, AliStack *stack, TString calo, TLorentzVector momentum) | |
bdcfac30 | 321 | { |
4b707925 | 322 | //Gets the primary particle and do some checks: |
323 | //Check if primary is inside calorimeter and look the mother outsie | |
324 | //Check if mother is a decay photon, in which case check if decay was overlapped | |
325 | ||
326 | Float_t minangle = 0; | |
327 | Float_t ipdist = 0; | |
328 | TParticle * pmother = new TParticle(); | |
bdcfac30 | 329 | |
4b707925 | 330 | if(calo == "PHOS"){ |
331 | ipdist= fPHOSIPDistance; | |
332 | minangle = fPHOSMinAngle; | |
333 | } | |
334 | else if (calo == "EMCAL"){ | |
335 | ipdist = fEMCALIPDistance; | |
336 | minangle = fEMCALMinAngle; | |
337 | } | |
bdcfac30 | 338 | |
bdcfac30 | 339 | |
4b707925 | 340 | if(label>=0){ |
341 | pmother = stack->Particle(label); | |
342 | Int_t mostatus = pmother->GetStatusCode(); | |
343 | Int_t mopdg = TMath::Abs(pmother->GetPdgCode()); | |
344 | ||
345 | //Check if mother is a conversion inside the calorimeter | |
346 | //In such case, return the mother before the calorimeter. | |
347 | //First approximation. | |
348 | Float_t vy = TMath::Abs(pmother->Vy()) ; | |
bdcfac30 | 349 | |
4b707925 | 350 | if( mostatus == 0 && vy >= ipdist){ |
bdcfac30 | 351 | |
4b707925 | 352 | //cout<<"Calo: "<<calo<<" vy "<<vy<<" E clus "<<momentum.E()<<" Emother "<<pmother->Energy()<<" " |
353 | // <<pmother->GetName()<<endl; | |
354 | ||
355 | while( vy >= ipdist){//inside calorimeter | |
356 | AliDebug(4,Form("Conversion inside calorimeter, mother vertex %0.2f, IP distance %0.2f", vy, ipdist)); | |
357 | pmother = stack->Particle(pmother->GetMother(0)); | |
358 | vy = TMath::Abs(pmother->Vy()) ; | |
359 | //cout<<" label "<<label<<" Mother: "<<pmother->GetName()<<" E "<<pmother->Energy()<<" Status "<<pmother->GetStatusCode()<<" and vertex "<<vy<<endl; | |
360 | mostatus = pmother->GetStatusCode(); | |
361 | mopdg = TMath::Abs(pmother->GetPdgCode()); | |
362 | }//while vertex is inside calorimeter | |
363 | //cout<<"Calo: "<<calo<<" final vy "<<vy<<" E clus "<<momentum.E()<<" Emother "<<pmother->Energy()<<" " | |
364 | // <<pmother->GetName()<<endl; | |
365 | }//check status and vertex | |
366 | ||
367 | AliDebug(4,Form("%s, ESD E %2.2f, PID %d, mother: E %2.2f, label %d, status %d, vertex %3.2f, mother 2 %d, grandmother %d \n", | |
368 | calo.Data(),momentum.E(),pmother->Energy(), label, pmother->GetPdgCode(), | |
369 | pmother->GetStatusCode(), vy, pmother->GetMother(0), stack->GetPrimary(label))); | |
370 | ||
371 | //Check Decay photons | |
372 | if(mopdg == 22){ | |
373 | ||
374 | //his mother was a pi0? | |
375 | TParticle * pmotherpi0 = stack->Particle(pmother->GetMother(0)); | |
376 | if( pmotherpi0->GetPdgCode() == 111){ | |
377 | ||
378 | AliDebug(4,Form(" %s: E cluster %2.2f, E gamma %2.2f, Mother Pi0, E %0.2f, status %d, daughters %d\n", | |
379 | calo.Data(), momentum.E(),pmother->Energy(),pmotherpi0->Energy(), | |
380 | pmotherpi0->GetStatusCode(), pmotherpi0->GetNDaughters())); | |
381 | ||
382 | if(pmotherpi0->GetNDaughters() == 1) mostatus = 2 ; //signal that this photon can only be decay, not overlapped. | |
383 | else if(pmotherpi0->GetNDaughters() == 2){ | |
384 | ||
385 | TParticle * pd1 = stack->Particle(pmotherpi0->GetDaughter(0)); | |
386 | TParticle * pd2 = stack->Particle(pmotherpi0->GetDaughter(1)); | |
387 | //if(pmotherpi0->Energy()> 10 ){ | |
388 | // cout<<"Two "<<calo<<" daugthers, pi0 label "<<pmother->GetMother(0)<<" E :"<<pmotherpi0->Energy()<<endl; | |
389 | // cout<<" 1) label "<<pmotherpi0->GetDaughter(0)<<" pdg "<<pd1->GetPdgCode()<<" E "<<pd1->Energy() | |
390 | // <<" phi "<<pd1->Phi()*TMath::RadToDeg()<<" eta "<<pd1->Eta() | |
391 | // <<" mother label "<<pd1->GetMother(0)<<" n daug "<<pd1->GetNDaughters() <<endl; | |
392 | // cout<<" 2) label "<<pmotherpi0->GetDaughter(1)<<" pdg "<<pd2->GetPdgCode()<<" E "<<pd2->Energy() | |
393 | // <<" phi "<<pd2->Phi()*TMath::RadToDeg()<<" eta "<<pd2->Eta()<<" mother label " | |
394 | // <<pd2->GetMother(0)<<" n daug "<<pd2->GetNDaughters() <<endl; | |
395 | //} | |
396 | if(pd1->GetPdgCode() == 22 && pd2->GetPdgCode() == 22){ | |
397 | TLorentzVector lv1 , lv2 ; | |
398 | pd1->Momentum(lv1); | |
399 | pd2->Momentum(lv2); | |
400 | Double_t angle = lv1.Angle(lv2.Vect()); | |
401 | // if(pmotherpi0->Energy()> 10 ) | |
402 | // cout<<"angle*ipdist "<<angle*ipdist<<" mindist "<< minangle <<endl; | |
403 | if (angle < minangle){ | |
404 | //There is overlapping, pass the pi0 as mother | |
405 | ||
406 | AliDebug(4,Form(">>>> %s cluster with E %2.2f is a overlapped pi0, E %2.2f, angle %2.4f < anglemin %2.4f\n", | |
407 | calo.Data(), momentum.E(), pmotherpi0->Energy(), angle, minangle)); | |
408 | ||
409 | pmother = pmotherpi0 ; | |
410 | ||
411 | } | |
412 | else mostatus = 2 ; // gamma decay not overlapped | |
413 | }// daughters are photons | |
414 | else mostatus = 2; // daughters are e-gamma or e-e, no overlapped, or charged cluster | |
415 | }//2 daughters | |
416 | else AliDebug(4,Form("pi0 has %d daughters",pmotherpi0->GetNDaughters())); | |
417 | }//pi0 decay photon? | |
418 | }//photon | |
419 | ||
420 | pmother->SetStatusCode(mostatus); // if status = 2, decay photon. | |
421 | ||
422 | }//label >=0 | |
423 | else AliInfo(Form("Negative Kinematic label of PHOS cluster: %d",label)); | |
424 | ||
425 | return pmother ; | |
bdcfac30 | 426 | |
427 | } |