Compilation on Windows/Cygwin
[u/mrichter/AliRoot.git] / PWG4 / AliGammaMCDataReader.cxx
CommitLineData
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
57ClassImp(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//____________________________________________________________________________
71AliGammaMCDataReader::AliGammaMCDataReader(const AliGammaMCDataReader & g) :
4b707925 72 AliGammaReader(g)
bdcfac30 73{
74 // cpy ctor
75}
76
77//_________________________________________________________________________
78AliGammaMCDataReader & 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//____________________________________________________________________________
90void 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
320TParticle * 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}