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 //_________________________________________________________________________
19 // Class for reading data (ESDs) in order to do prompt gamma
20 // or other particle identification and correlations
22 // It is a filtering class, transforms ESD tracks or CaloClusters
23 // into AOD tracks and calocluters, which are the basic input of the analysis
24 // classes in this frame.
25 // It is recommended to use the official filter AliAnalysisTaskESDfilter, and
26 // then the reader for AODs AliCaloTrackAODReader.
29 //*-- Author: Gustavo Conesa (LNF-INFN)
30 //////////////////////////////////////////////////////////////////////////////
33 // --- ROOT system ---
34 //#include "Riostream.h"
36 //---- ANALYSIS system ----
37 #include "AliCaloTrackESDReader.h"
38 #include "AliESDEvent.h"
39 #include "AliESDCaloCluster.h"
40 #include "AliAODCaloCluster.h"
41 #include "AliAODTrack.h"
42 #include "AliAODPid.h"
43 #include "AliAODEvent.h"
44 #include "AliFidutialCut.h"
46 ClassImp(AliCaloTrackESDReader)
48 //____________________________________________________________________________
49 AliCaloTrackESDReader::AliCaloTrackESDReader() :
54 //Initialize parameters
57 fReadAODMCParticles = kFALSE;
58 //We want tracks fitted in the detectors:
59 fTrackStatus=AliESDtrack::kTPCrefit;
60 fTrackStatus|=AliESDtrack::kITSrefit;
64 //____________________________________________________________________________
65 AliCaloTrackESDReader::AliCaloTrackESDReader(const AliCaloTrackESDReader & g) :
71 //_________________________________________________________________________
72 //AliCaloTrackESDReader & AliCaloTrackESDReader::operator = (const AliCaloTrackESDReader & source)
74 // // assignment operator
76 // if(&source == this) return *this;
82 //____________________________________________________________________________
83 void AliCaloTrackESDReader::FillInputCTS() {
84 //Return array with CTS tracks
85 //TObjArray * fAODCTS = new TObjArray();
86 Int_t nTracks = fInputEvent->GetNumberOfTracks() ;
92 Double_t bfield = ((AliESDEvent*)fInputEvent)->GetMagneticField();
94 Double_t timezero = 0; //TO BE FIXED
96 //To be replaced by call to AliEMCALGeoUtils when the class becomes available
97 Double_t radius = 441.0; //[cm] EMCAL radius +13cm
99 for (Int_t itrack = 0; itrack < nTracks; itrack++) {////////////// track loop
100 AliESDtrack * track = (AliESDtrack*) ((AliESDEvent*)fInputEvent)->GetTrack(itrack) ; // retrieve track from esd
102 //We want tracks fitted in the detectors: TPCrefit, ITSrefit ... check the set bits.
103 if (fTrackStatus && !((track->GetStatus() & fTrackStatus) == fTrackStatus)) continue ;
105 track->GetPxPyPz(p) ;
106 TLorentzVector momentum(p[0],p[1],p[2],0);
108 if(fCTSPtMin < momentum.Pt() &&fFidutialCut->IsInFidutialCut(momentum,"CTS")){
110 if(fDebug > 3 && momentum.Pt() > 0.2) printf("AliCaloTrackESDReader::FillInputCTS() - Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
111 momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
114 track->GetCovarianceXYZPxPyPz(covTr);
115 track->GetESDpid(pid);
117 Float_t impactXY, impactZ;
119 track->GetImpactParameters(impactXY,impactZ);
122 // track inside the beam pipe
123 //Put new aod object in file in AOD tracks array
124 AliAODTrack *aodTrack = new((*(fOutputEvent->GetTracks()))[naod++])
125 AliAODTrack(track->GetID(), track->GetLabel(), p, kTRUE, pos, kFALSE,covTr, (Short_t)track->GetSign(), track->GetITSClusterMap(),
128 kTRUE, // check if this is right
129 kTRUE, // check if this is right
130 AliAODTrack::kPrimary,
133 aodTrack->SetFlags(track->GetStatus());
134 aodTrack->ConvertAliPIDtoAODPID();
136 //fill needed AliAODPid information, including
137 //Extrapolate track to EMCAL surface for AOD-level track-cluster matching
138 AliAODPid *aodpid = new AliAODPid;
139 aodpid->SetITSsignal(track->GetITSsignal());
140 aodpid->SetTPCsignal(track->GetTPCsignal());
142 Int_t nslices = track->GetNumberOfTRDslices()*6;
143 Double_t *trdslices = new Double_t[nslices];
144 for(Int_t iSl =0; iSl < track->GetNumberOfTRDslices(); iSl++) {
145 for(Int_t iPl =0; iPl<6; iPl++) trdslices[iPl*track->GetNumberOfTRDslices()+iSl] = track->GetTRDslice(iPl,iSl);
147 aodpid->SetTRDsignal(track->GetNumberOfTRDslices()*6,trdslices);
148 Double_t times[AliAODPid::kSPECIES]; track->GetIntegratedTimes(times);
149 aodpid->SetIntegratedTimes(times);
151 aodpid->SetTOFsignal(track->GetTOFsignal()-timezero); // to be fixed
152 aodpid->SetHMPIDsignal(track->GetHMPIDsignal());
154 Double_t emcpos[3] = {0.,0.,0.};
155 Double_t emcmom[3] = {0.,0.,0.};
156 aodpid->SetEMCALPosition(emcpos);
157 aodpid->SetEMCALMomentum(emcmom);
159 AliExternalTrackParam *outerparam = (AliExternalTrackParam*)track->GetOuterParam();
160 if(!outerparam) return;
162 Bool_t okpos = outerparam->GetXYZAt(radius,bfield,emcpos);
163 Bool_t okmom = outerparam->GetPxPyPzAt(radius,bfield,emcmom);
164 if(!(okpos && okmom)) return;
166 aodpid->SetEMCALPosition(emcpos);
167 aodpid->SetEMCALMomentum(emcmom);
169 aodTrack->SetDetPID(aodpid);
171 else continue; // outside the beam pipe: orphan track
172 }//Pt and Fidutial cut passed.
175 //Put references to selected tracks in array
176 for(Int_t itrack = 0; itrack < (fOutputEvent->GetTracks())->GetEntriesFast(); itrack++){
177 AliAODTrack * track = (AliAODTrack*) (fOutputEvent->GetTracks())->At(itrack);
181 if(fDebug > 1) printf("AliCaloTrackESDReader::FillInputCTS() - aod entries %d\n", fAODCTS->GetEntriesFast());
184 //____________________________________________________________________________
185 void AliCaloTrackESDReader::FillInputEMCAL() {
186 //Return array with EMCAL clusters in aod format
188 //Get vertex for momentum calculation
189 Double_t v[3] ; //vertex ;
193 Int_t naod = (fOutputEvent->GetCaloClusters())->GetEntriesFast();
194 Int_t nTracks = fInputEvent->GetNumberOfTracks() ;
196 //Loop to select clusters in fidutial cut and fill container with aodClusters
197 for (Int_t iclus = 0; iclus < ((AliESDEvent*)fInputEvent)->GetNumberOfCaloClusters(); iclus++) {
198 AliESDCaloCluster * clus = 0;
199 if ( (clus = ((AliESDEvent*)fInputEvent)->GetCaloCluster(iclus)) ) {
200 if (clus->IsEMCAL()){
201 TLorentzVector momentum ;
202 clus->GetMomentum(momentum, v);
203 if(fDebug > 3 && momentum.E() > 0.2) printf("AliCaloTrackESDReader::FillInputEMCAL() - all clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
204 momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
205 if(fEMCALPtMin < momentum.Pt() && fFidutialCut->IsInFidutialCut(momentum,"EMCAL")){
207 if(fDebug > 2 && momentum.E() > 0.2) printf("AliCaloTrackESDReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
208 momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
210 clus->GetPosition(pos) ;
211 Int_t id = clus->GetID();
212 Int_t nLabel = clus->GetNLabels();
214 if(clus->GetLabels()) labels = (clus->GetLabels())->GetArray();
216 Float_t energy = clus->E();
217 Char_t ttype= AliAODCluster::kEMCALClusterv1;
219 //Put new aod object in file in AOD calo clusters array
220 AliAODCaloCluster *caloCluster = new((*(fOutputEvent->GetCaloClusters()))[naod++])
221 AliAODCaloCluster(id,nLabel,labels,energy, pos, NULL,ttype,0);
223 // printf("Reader PID ESD: 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 \n",
224 // pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
225 // pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]);
226 caloCluster->SetPIDFromESD(clus->GetPid());
227 caloCluster->SetCaloCluster(clus->GetDistanceToBadChannel(), clus->GetClusterDisp(),
228 clus->GetM20(), clus->GetM02(),
229 clus->GetEmcCpvDistance(), clus->GetNExMax(), clus->GetTOF()) ;
232 if(fDebug > 3 && momentum.E() > 0.2)
233 printf("AliCaloTrackESDReader::FillInputEMCAL() - Selected clusters Distance BC %2.2f, dispersion %2.2f, M20 %f, M02 %3.2f, NexMax %d, TOF %e\n",
234 clus->GetDistanceToBadChannel(), clus->GetClusterDisp(),clus->GetM20(), clus->GetM02(),
235 clus->GetNExMax(), clus->GetTOF());
237 caloCluster->SetNCells(clus->GetNCells());
238 caloCluster->SetCellsAbsId(clus->GetCellsAbsId());
239 caloCluster->SetCellsAmplitudeFraction(clus->GetCellsAmplitudeFraction());
241 TArrayI* matchedT = clus->GetTracksMatched();
242 if (nTracks > 0 && matchedT && clus->GetTrackMatched() >= 0) {
243 for (Int_t im = 0; im < matchedT->GetSize(); im++) {
244 Int_t iESDtrack = matchedT->At(im);
245 if(iESDtrack < nTracks && iESDtrack > -1)
246 caloCluster->AddTrackMatched((fInputEvent->GetTrack(iESDtrack)));
249 //Fill reference array
250 }//Pt and Fidutial cut passed.
255 //Put references to selected clusters in array
256 for(Int_t iclus = 0; iclus < (fOutputEvent->GetCaloClusters())->GetEntriesFast(); iclus++){
257 AliAODCaloCluster * clus = (AliAODCaloCluster*) (fOutputEvent->GetCaloClusters())->At(iclus);
258 fAODEMCAL->Add(clus);
260 if(fDebug > 1) printf("AliCaloTrackESDReader::FillInputEMCAL() - aod entries %d\n", fAODEMCAL->GetEntriesFast());
264 //____________________________________________________________________________
265 void AliCaloTrackESDReader::FillInputPHOS() {
266 //Return array with PHOS clusters in aod format
267 Int_t nEMCAL = (fOutputEvent->GetCaloClusters())->GetEntriesFast();
268 //Get vertex for momentum calculation
269 Double_t v[3] ; //vertex ;
273 Double_t * pid = new Double_t[AliPID::kSPECIESN];
274 Int_t naod = (fOutputEvent->GetCaloClusters())->GetEntriesFast();
275 Int_t nTracks = fInputEvent->GetNumberOfTracks() ;
277 //Loop to select clusters in fidutial cut and fill container with aodClusters
278 for (Int_t iclus = 0; iclus < ((AliESDEvent*)fInputEvent)->GetNumberOfCaloClusters(); iclus++) {
279 AliESDCaloCluster * clus = 0;
280 if ( (clus = ((AliESDEvent*)fInputEvent)->GetCaloCluster(iclus)) ) {
282 TLorentzVector momentum ;
283 clus->GetMomentum(momentum, v);
284 if(fDebug > 3 && momentum.E() > 0.1)printf("AliCaloTrackESDReader::FillInputPHOS() - all clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
285 momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
286 if(fPHOSPtMin < momentum.Pt() && fFidutialCut->IsInFidutialCut(momentum,"PHOS")){
288 if(fDebug > 2 && momentum.E() > 0.1)printf("AliCaloTrackESDReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
289 momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
291 clus->GetPosition(pos) ;
292 Int_t id = clus->GetID();
293 Int_t nLabel = clus->GetNLabels();
295 if(clus->GetLabels()) labels = (clus->GetLabels())->GetArray();
296 Float_t energy = clus->E();
299 pid = clus->GetPid();
300 // printf("Reader PID ESD: 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 \n",
301 // pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
302 // pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]);
303 Char_t ttype= AliAODCluster::kPHOSNeutral;
304 Float_t wNeutral = pid[AliPID::kNeutron]+ pid[AliPID::kKaon0]+pid[AliPID::kPhoton]+pid[AliPID::kPi0];
305 Float_t wCharged = pid[AliPID::kMuon] + pid[AliPID::kElectron] + pid[AliPID::kEleCon]+
306 pid[AliPID::kProton]+pid[AliPID::kKaon]+pid[AliPID::kPion];
307 if( wCharged > wNeutral) ttype= AliAODCluster::kPHOSCharged;
309 //Put new aod object in file in AOD calo clusters array
310 AliAODCaloCluster *caloCluster = new((*(fOutputEvent->GetCaloClusters()))[naod++])
311 AliAODCaloCluster(id,nLabel,labels,energy, pos, NULL, ttype, 0);
314 caloCluster->SetPIDFromESD(pid);
315 caloCluster->SetCaloCluster(clus->GetDistanceToBadChannel(), clus->GetClusterDisp(),
316 clus->GetM20(), clus->GetM02(),
317 clus->GetEmcCpvDistance(), clus->GetNExMax()) ;
319 if(fDebug > 3 && momentum.E() > 0.2)
320 printf("AliCaloTrackESDReader::FillInputPHOS() - Selected clusters Distance BC %2.2f, dispersion %2.2f, M20 %f, M02 %3.2f, EmcCpvDist %3.3f, NexMax %d, TOF %e\n",
321 clus->GetDistanceToBadChannel(), clus->GetClusterDisp(),clus->GetM20(), clus->GetM02(),
322 clus->GetEmcCpvDistance(), clus->GetNExMax(), clus->GetTOF());
324 caloCluster->SetNCells(clus->GetNCells());
325 caloCluster->SetCellsAbsId(clus->GetCellsAbsId());
326 caloCluster->SetCellsAmplitudeFraction(clus->GetCellsAmplitudeFraction());
328 TArrayI* matchedT = clus->GetTracksMatched();
329 if (nTracks > 0 && matchedT && clus->GetTrackMatched() >= 0) {
330 for (Int_t im = 0; im < matchedT->GetSize(); im++) {
331 Int_t iESDtrack = matchedT->At(im);
332 if(iESDtrack < nTracks && iESDtrack > -1)
333 caloCluster->AddTrackMatched((fInputEvent->GetTrack(iESDtrack)));
336 }//Pt and Fidutial cut passed.
341 //Put references to selected clusters in array
342 for(Int_t iclus = nEMCAL; iclus < (fOutputEvent->GetCaloClusters())->GetEntriesFast(); iclus++){
343 AliAODCaloCluster * clus = (AliAODCaloCluster*) (fOutputEvent->GetCaloClusters())->At(iclus);
346 if(fDebug > 1) printf("AliCaloTrackESDReader::FillInputPHOS() - aod entries %d\n", fAODPHOS->GetEntriesFast());
350 //____________________________________________________________________________
351 void AliCaloTrackESDReader::FillInputEMCALCells() {
352 //Return array with EMCAL cells in esd format
354 fEMCALCells = (TNamed*) ((AliESDEvent*)fInputEvent)->GetEMCALCells();
358 //____________________________________________________________________________
359 void AliCaloTrackESDReader::FillInputPHOSCells() {
360 //Return array with PHOS cells in esd format
362 fPHOSCells = (TNamed*) ((AliESDEvent*)fInputEvent)->GetPHOSCells();
366 //____________________________________________________________________________
367 void AliCaloTrackESDReader::GetVertex(Double_t v[3]) const {
368 //Return vertex position
370 //((AliESDEvent*)fInputEvent)->GetVertex()->GetXYZ(v) ;//SPD
371 ((AliESDEvent*)fInputEvent)->GetPrimaryVertex()->GetXYZ(v);
376 //____________________________________________________________________________
377 Double_t AliCaloTrackESDReader::GetBField() const {
378 //Return magnetic field
380 Double_t bfield = ((AliESDEvent*)fInputEvent)->GetMagneticField();
386 //____________________________________________________________________________
387 void AliCaloTrackESDReader::SetInputOutputMCEvent(AliVEvent* esd, AliAODEvent* aod, AliMCEvent* mc) {
388 // Connect the data pointers
390 if(strcmp(esd->GetName(),"AliESDEvent")){
391 printf("AliCaloTrackESDReader::SetInputOutputMCEvent() - STOP ::Wrong reader, here only ESDs. Input name: %s != AliESDEvent \n",esd->GetName());