2 /**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 /* History of cvs commits:
21 * Revision 1.2 2007/08/17 12:40:04 schutz
22 * New analysis classes by Gustavo Conesa
24 * Revision 1.1.2.1 2007/07/26 10:32:09 schutz
25 * new analysis classes in the the new analysis framework
30 //_________________________________________________________________________
31 // Class for reading data (Kinematics and ESDs) in order to do prompt gamma correlations
32 // Class created from old AliPHOSGammaJet
33 // (see AliRoot versions previous Release 4-09)
35 //*-- Author: Gustavo Conesa (LNF-INFN)
36 //////////////////////////////////////////////////////////////////////////////
39 // --- ROOT system ---
41 #include <TParticle.h>
43 //---- ANALYSIS system ----
44 #include "AliGammaMCDataReader.h"
45 #include "AliESDEvent.h"
46 #include "AliESDVertex.h"
47 #include "AliESDCaloCluster.h"
51 ClassImp(AliGammaMCDataReader)
53 //____________________________________________________________________________
54 AliGammaMCDataReader::AliGammaMCDataReader() :
59 //Initialize parameters
64 //____________________________________________________________________________
65 AliGammaMCDataReader::AliGammaMCDataReader(const AliGammaMCDataReader & g) :
71 //_________________________________________________________________________
72 AliGammaMCDataReader & AliGammaMCDataReader::operator = (const AliGammaMCDataReader & source)
74 // assignment operator
76 if(&source == this) return *this;
83 //____________________________________________________________________________
84 void AliGammaMCDataReader::CreateParticleList(TObject * data, TObject * kine,
86 TClonesArray * plEMCAL,
87 TClonesArray * plPHOS,
88 TClonesArray * plPrimCTS,
89 TClonesArray * plPrimEMCAL,
90 TClonesArray * plPrimPHOS){
92 //Create a list of particles from the ESD. These particles have been measured
93 //by the Central Tracking system (TPC+ITS+...), PHOS and EMCAL
94 //Also create particle list with mothers.
96 AliESDEvent* esd = (AliESDEvent*) data;
97 AliStack* stack = (AliStack*) kine;
100 Double_t *pid = new Double_t[AliPID::kSPECIESN];
101 AliDebug(3,"Fill particle lists");
103 //Get vertex for momentum calculation
104 Double_t v[3] ; //vertex ;
105 esd->GetVertex()->GetXYZ(v) ;
107 //########### CALORIMETERS ##############
108 Int_t nCaloCluster = esd->GetNumberOfCaloClusters() ;
109 Int_t indexPH = plPHOS->GetEntries() ;
110 Int_t indexEM = plEMCAL->GetEntries() ;
112 for (npar = 0; npar < nCaloCluster; npar++) {//////////////CaloCluster loop
113 AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve cluster from esd
114 Int_t type = clus->GetClusterType();
116 //########### PHOS ##############
117 if(fSwitchOnPHOS && type == AliESDCaloCluster::kPHOSCluster){
118 AliDebug(4,Form("PHOS clusters: E %f, match %d", clus->E(),clus->GetTrackMatched()));
120 if(clus->GetTrackMatched()==-1){
121 TLorentzVector momentum ;
122 clus->GetMomentum(momentum, v);
123 Double_t phi = momentum.Phi();
124 if(phi<0) phi+=TMath::TwoPi() ;
125 if(momentum.Pt() > fNeutralPtCut && TMath::Abs(momentum.Eta()) < fPHOSEtaCut &&
126 phi > fPhiPHOSCut[0] && phi < fPhiPHOSCut[1] ) {
132 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 ",
133 momentum.E(),pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
134 pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]));
136 Float_t wPhoton = fPHOSPhotonWeight;
137 Float_t wPi0 = fPHOSPi0Weight;
139 if(fPHOSWeightFormula){
140 wPhoton = fPHOSPhotonWeightFormula->Eval(momentum.E()) ;
141 wPi0 = fPHOSPi0WeightFormula->Eval(momentum.E());
144 if(pid[AliPID::kPhoton] > wPhoton)
146 else if(pid[AliPID::kPi0] > wPi0)
148 else if(pid[AliPID::kElectron] > fPHOSElectronWeight)
150 else if(pid[AliPID::kEleCon] > fPHOSElectronWeight)
152 else if(pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] > fPHOSChargeWeight)
153 pdg = kChargedHadron ;
154 else if(pid[AliPID::kKaon0]+pid[AliPID::kNeutron] > fPHOSNeutralWeight)
155 pdg = kNeutralHadron ;
157 else if(pid[AliPID::kElectron]+pid[AliPID::kEleCon]+pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] >
158 pid[AliPID::kPhoton] + pid[AliPID::kPi0]+pid[AliPID::kKaon0]+pid[AliPID::kNeutron])
159 pdg = kChargedUnknown ;
161 pdg = kNeutralUnknown ;
162 //neutral cluster, unidentifed.
165 if(pdg != kElectron && pdg != kEleCon && pdg !=kChargedHadron && pdg !=kChargedUnknown ){//keep only neutral particles in the array
170 Int_t label = clus->GetLabel();
171 TParticle * pmother = GetMotherParticle(label,stack, "PHOS",momentum);
173 TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
174 momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
176 AliDebug(4,Form("PHOS added: pdg %d, pt %f, phi %f, eta %f", pdg, particle->Pt(),particle->Phi(),particle->Eta()));
177 new((*plPHOS)[indexPH]) TParticle(*particle) ;
178 new((*plPrimPHOS)[indexPH]) TParticle(*pmother) ;
182 else AliDebug(4,Form("PHOS charged cluster NOT added: pdg %d, pt %f, phi %f, eta %f\n",
183 pdg, momentum.Pt(),momentum.Phi(),momentum.Eta()));
186 else AliDebug(4,"Particle not added");
190 //################ EMCAL ##############
191 else if(fSwitchOnEMCAL && type == AliESDCaloCluster::kEMCALClusterv1){
192 AliDebug(4,Form("EMCAL clusters: E %f, match %d", clus->E(),clus->GetTrackMatched()));
194 if(clus->GetTrackMatched()==-1 ){
195 TLorentzVector momentum ;
196 clus->GetMomentum(momentum, v);
197 Double_t phi = momentum.Phi();
198 if(phi<0) phi+=TMath::TwoPi() ;
199 if(momentum.Pt() > fNeutralPtCut && TMath::Abs(momentum.Eta()) < fEMCALEtaCut &&
200 phi > fPhiEMCALCut[0] && phi < fPhiEMCALCut[1] ) {
206 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 ",
207 momentum.E(),pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
208 pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]));
210 if(pid[AliPID::kPhoton] > fEMCALPhotonWeight)
212 else if(pid[AliPID::kPi0] > fEMCALPi0Weight)
214 else if(pid[AliPID::kElectron] > fEMCALElectronWeight)
216 else if(pid[AliPID::kEleCon] > fEMCALElectronWeight)
218 else if(pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] > fEMCALChargeWeight)
219 pdg = kChargedHadron ;
220 else if(pid[AliPID::kKaon0]+pid[AliPID::kNeutron] > fEMCALNeutralWeight)
221 pdg = kNeutralHadron ;
222 else if(pid[AliPID::kElectron]+pid[AliPID::kEleCon]+pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] >
223 pid[AliPID::kPhoton] + pid[AliPID::kPi0]+pid[AliPID::kKaon0]+pid[AliPID::kNeutron])
224 pdg = kChargedUnknown ;
226 pdg = kNeutralUnknown ;
229 if(pdg != kElectron && pdg != kEleCon && pdg !=kChargedHadron && pdg !=kChargedUnknown){//keep only neutral particles in the array
234 Int_t label = clus->GetLabel();
235 TParticle * pmother = GetMotherParticle(label,stack, "EMCAL",momentum);
237 TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
238 momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
239 AliDebug(4,Form("EMCAL cluster added: pdg %f, pt %f, phi %f, eta %f", pdg, particle->Pt(),particle->Phi(),particle->Eta()));
241 new((*plEMCAL)[indexEM]) TParticle(*particle) ;
242 new((*plPrimEMCAL)[indexEM]) TParticle(*pmother) ;
245 else AliDebug(4,Form("EMCAL charged cluster NOT added: pdg %d, pt %f, phi %f, eta %f",
246 pdg, momentum.Pt(),momentum.Phi(),momentum.Eta()));
249 else AliDebug(4,"Particle not added");
257 //########### CTS (TPC+ITS) #####################
258 Int_t nTracks = esd->GetNumberOfTracks() ;
259 Int_t indexCh = plCTS->GetEntries() ;
262 AliDebug(3,Form("Number of tracks %d",nTracks));
264 for (npar = 0; npar < nTracks; npar++) {////////////// track loop
265 AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd
267 //We want tracks fitted in the detectors:
268 ULong_t status=AliESDtrack::kTPCrefit;
269 status|=AliESDtrack::kITSrefit;
271 //We want tracks whose PID bit is set:
272 // ULong_t status =AliESDtrack::kITSpid;
273 // status|=AliESDtrack::kTPCpid;
275 if ( (track->GetStatus() & status) == status) {//Check if the bits we want are set
276 // Do something with the tracks which were successfully
278 Double_t en = 0; //track ->GetTPCsignal() ;
280 track->GetPxPyPz(mom) ;
281 Double_t px = mom[0];
282 Double_t py = mom[1];
283 Double_t pz = mom[2]; //Check with TPC people if this is correct.
284 Int_t pdg = 11; //Give any charged PDG code, in this case electron.
285 //I just want to tag the particle as charged
288 Int_t label = TMath::Abs(track->GetLabel());
289 TParticle * pmother = new TParticle();
290 pmother = stack->Particle(label);
292 TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
293 px, py, pz, en, v[0], v[1], v[2], 0);
295 if(particle->Pt() > fChargedPtCut && TMath::Abs(particle->Eta())<fCTSEtaCut){
296 new((*plCTS)[indexCh]) TParticle(*particle) ;
297 new((*plPrimCTS)[indexCh]) TParticle(*pmother) ;
300 }//kinematic selection
301 }//select track from refit
305 AliDebug(3,Form("Particle lists filled, tracks %d , clusters: EMCAL %d, PHOS %d", indexCh,indexEM,indexPH));
309 TParticle * AliGammaMCDataReader:: GetMotherParticle(Int_t label, AliStack *stack, TString calo, TLorentzVector momentum)
311 //Gets the primary particle and do some checks:
312 //Check if primary is inside calorimeter and look the mother outsie
313 //Check if mother is a decay photon, in which case check if decay was overlapped
315 Float_t minangle = 0;
317 TParticle * pmother = new TParticle();
320 ipdist= fPHOSIPDistance;
321 minangle = fPHOSMinAngle;
323 else if (calo == "EMCAL"){
324 ipdist = fEMCALIPDistance;
325 minangle = fEMCALMinAngle;
330 pmother = stack->Particle(label);
331 Int_t mostatus = pmother->GetStatusCode();
332 Int_t mopdg = TMath::Abs(pmother->GetPdgCode());
334 //Check if mother is a conversion inside the calorimeter
335 //In such case, return the mother before the calorimeter.
336 //First approximation.
337 Float_t vy = TMath::Abs(pmother->Vy()) ;
339 if( mostatus == 0 && vy >= ipdist){
341 //cout<<"Calo: "<<calo<<" vy "<<vy<<" E clus "<<momentum.E()<<" Emother "<<pmother->Energy()<<" "
342 // <<pmother->GetName()<<endl;
344 while( vy >= ipdist){//inside calorimeter
345 AliDebug(4,Form("Conversion inside calorimeter, mother vertex %0.2f, IP distance %0.2f", vy, ipdist));
346 pmother = stack->Particle(pmother->GetMother(0));
347 vy = TMath::Abs(pmother->Vy()) ;
348 //cout<<" label "<<label<<" Mother: "<<pmother->GetName()<<" E "<<pmother->Energy()<<" Status "<<pmother->GetStatusCode()<<" and vertex "<<vy<<endl;
349 mostatus = pmother->GetStatusCode();
350 mopdg = TMath::Abs(pmother->GetPdgCode());
351 }//while vertex is inside calorimeter
352 //cout<<"Calo: "<<calo<<" final vy "<<vy<<" E clus "<<momentum.E()<<" Emother "<<pmother->Energy()<<" "
353 // <<pmother->GetName()<<endl;
354 }//check status and vertex
356 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",
357 calo.Data(),momentum.E(),pmother->Energy(), label, pmother->GetPdgCode(),
358 pmother->GetStatusCode(), vy, pmother->GetMother(0), stack->GetPrimary(label)));
360 //Check Decay photons
363 //his mother was a pi0?
364 TParticle * pmotherpi0 = stack->Particle(pmother->GetMother(0));
365 if( pmotherpi0->GetPdgCode() == 111){
367 AliDebug(4,Form(" %s: E cluster %2.2f, E gamma %2.2f, Mother Pi0, E %0.2f, status %d, daughters %d\n",
368 calo.Data(), momentum.E(),pmother->Energy(),pmotherpi0->Energy(),
369 pmotherpi0->GetStatusCode(), pmotherpi0->GetNDaughters()));
371 if(pmotherpi0->GetNDaughters() == 1) mostatus = 2 ; //signal that this photon can only be decay, not overlapped.
372 else if(pmotherpi0->GetNDaughters() == 2){
374 TParticle * pd1 = stack->Particle(pmotherpi0->GetDaughter(0));
375 TParticle * pd2 = stack->Particle(pmotherpi0->GetDaughter(1));
376 //if(pmotherpi0->Energy()> 10 ){
377 // cout<<"Two "<<calo<<" daugthers, pi0 label "<<pmother->GetMother(0)<<" E :"<<pmotherpi0->Energy()<<endl;
378 // cout<<" 1) label "<<pmotherpi0->GetDaughter(0)<<" pdg "<<pd1->GetPdgCode()<<" E "<<pd1->Energy()
379 // <<" phi "<<pd1->Phi()*TMath::RadToDeg()<<" eta "<<pd1->Eta()
380 // <<" mother label "<<pd1->GetMother(0)<<" n daug "<<pd1->GetNDaughters() <<endl;
381 // cout<<" 2) label "<<pmotherpi0->GetDaughter(1)<<" pdg "<<pd2->GetPdgCode()<<" E "<<pd2->Energy()
382 // <<" phi "<<pd2->Phi()*TMath::RadToDeg()<<" eta "<<pd2->Eta()<<" mother label "
383 // <<pd2->GetMother(0)<<" n daug "<<pd2->GetNDaughters() <<endl;
385 if(pd1->GetPdgCode() == 22 && pd2->GetPdgCode() == 22){
386 TLorentzVector lv1 , lv2 ;
389 Double_t angle = lv1.Angle(lv2.Vect());
390 // if(pmotherpi0->Energy()> 10 )
391 // cout<<"angle*ipdist "<<angle*ipdist<<" mindist "<< minangle <<endl;
392 if (angle < minangle){
393 //There is overlapping, pass the pi0 as mother
395 AliDebug(4,Form(">>>> %s cluster with E %2.2f is a overlapped pi0, E %2.2f, angle %2.4f < anglemin %2.4f\n",
396 calo.Data(), momentum.E(), pmotherpi0->Energy(), angle, minangle));
398 pmother = pmotherpi0 ;
401 else mostatus = 2 ; // gamma decay not overlapped
402 }// daughters are photons
403 else mostatus = 2; // daughters are e-gamma or e-e, no overlapped, or charged cluster
405 else AliDebug(4,Form("pi0 has %d daughters",pmotherpi0->GetNDaughters()));
409 pmother->SetStatusCode(mostatus); // if status = 2, decay photon.
412 else AliInfo(Form("Negative Kinematic label of PHOS cluster: %d",label));