PWG4/PartCorrBase/AliCaloTrackESDReader.cxx

   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 //_________________________________________________________________________
  19 // Class for reading data (ESDs) in order to do prompt gamma
  20 // or other particle identification and correlations
  21 //
  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.
  27 //
  28 //
  29 //*-- Author: Gustavo Conesa (LNF-INFN)
  30 //////////////////////////////////////////////////////////////////////////////
  31
  32
  33 // --- ROOT system ---
  34 //#include "Riostream.h"
  35
  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"
  45
  46 ClassImp(AliCaloTrackESDReader)
  47
  48 //____________________________________________________________________________
  49 AliCaloTrackESDReader::AliCaloTrackESDReader() :
  50 AliCaloTrackReader()
  51 {
  52         //Default Ctor
  53
  54         //Initialize parameters
  55         fDataType=kESD;
  56         fReadStack          = kTRUE;
  57         fReadAODMCParticles = kFALSE;
  58         //We want tracks fitted in the detectors:
  59         fTrackStatus=AliESDtrack::kTPCrefit;
  60     fTrackStatus|=AliESDtrack::kITSrefit;
  61
  62 }
  63
  64 //____________________________________________________________________________
  65 AliCaloTrackESDReader::AliCaloTrackESDReader(const AliCaloTrackESDReader & g) :
  66   AliCaloTrackReader(g)
  67 {
  68   // cpy ctor
  69 }
  70
  71 //_________________________________________________________________________
  72 //AliCaloTrackESDReader & AliCaloTrackESDReader::operator = (const AliCaloTrackESDReader & source)
  73 //{
  74 //      // assignment operator
  75 //
  76 //      if(&source == this) return *this;
  77 //
  78 //      return *this;
  79 //
  80 //}
  81
  82 //____________________________________________________________________________
  83 void AliCaloTrackESDReader::FillInputCTS() {
  84   //Return array with CTS tracks
  85   //TObjArray * fAODCTS = new TObjArray();
  86   Int_t nTracks   = fInputEvent->GetNumberOfTracks() ;
  87   Int_t naod = 0;
  88   Double_t pos[3];
  89   Double_t p[3];
  90   Double_t covTr[21];
  91   Double_t pid[10];
  92   Double_t bfield = ((AliESDEvent*)fInputEvent)->GetMagneticField();
  93
  94   Double_t      timezero        = 0;   //TO BE FIXED
  95
  96   //To be replaced by call to AliEMCALGeoUtils when the class becomes available
  97   Double_t radius = 441.0; //[cm] EMCAL radius +13cm
  98   if(fDebug > 1) printf("AliCaloTrackESDReader::FillInputCTS() - org entries %d\n", nTracks);
  99   for (Int_t itrack =  0; itrack <  nTracks; itrack++) {////////////// track loop
 100     AliESDtrack * track = (AliESDtrack*) ((AliESDEvent*)fInputEvent)->GetTrack(itrack) ; // retrieve track from esd
 101
 102     //We want tracks fitted in the detectors: TPCrefit, ITSrefit ... check the set bits.
 103         if (fTrackStatus && !((track->GetStatus() & fTrackStatus) == fTrackStatus)) continue ;
 104
 105         track->GetPxPyPz(p) ;
 106         TLorentzVector momentum(p[0],p[1],p[2],0);
 107
 108         if(fCTSPtMin < momentum.Pt() &&fFidutialCut->IsInFidutialCut(momentum,"CTS")){
 109
 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());
 112
 113                 track->GetXYZ(pos);
 114                 track->GetCovarianceXYZPxPyPz(covTr);
 115                 track->GetESDpid(pid);
 116
 117                 Float_t impactXY, impactZ;
 118
 119                 track->GetImpactParameters(impactXY,impactZ);
 120
 121                 if (impactXY<3) {
 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(),
 126                                 pid,
 127                                 0x0,//primary,
 128                                 kTRUE, // check if this is right
 129                                 kTRUE, // check if this is right
 130                                 AliAODTrack::kPrimary,
 131                                 0);
 132
 133                   aodTrack->SetFlags(track->GetStatus());
 134                   aodTrack->ConvertAliPIDtoAODPID();
 135
 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());
 141                   //n TRD planes = 6
 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);
 146                   }
 147                   aodpid->SetTRDsignal(track->GetNumberOfTRDslices()*6,trdslices);
 148                   Double_t times[AliAODPid::kSPECIES]; track->GetIntegratedTimes(times);
 149                   aodpid->SetIntegratedTimes(times);
 150
 151                   aodpid->SetTOFsignal(track->GetTOFsignal()-timezero); // to be fixed
 152                   aodpid->SetHMPIDsignal(track->GetHMPIDsignal());
 153
 154                   Double_t emcpos[3] = {0.,0.,0.};
 155                   Double_t emcmom[3] = {0.,0.,0.};
 156                   aodpid->SetEMCALPosition(emcpos);
 157                   aodpid->SetEMCALMomentum(emcmom);
 158
 159                   AliExternalTrackParam *outerparam = (AliExternalTrackParam*)track->GetOuterParam();
 160                   if(!outerparam) continue;
 161
 162                   Bool_t okpos = outerparam->GetXYZAt(radius,bfield,emcpos);
 163                   Bool_t okmom = outerparam->GetPxPyPzAt(radius,bfield,emcmom);
 164                   if(!(okpos && okmom)) continue;
 165
 166                   aodpid->SetEMCALPosition(emcpos);
 167                   aodpid->SetEMCALMomentum(emcmom);
 168
 169                   aodTrack->SetDetPID(aodpid);
 170                 }
 171                 else continue;   // outside the beam pipe: orphan track
 172         }//Pt and Fidutial cut passed.
 173   }// track loop
 174
 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);
 178     fAODCTS->Add(track);
 179   }
 180
 181   if(fDebug > 1) printf("AliCaloTrackESDReader::FillInputCTS() - aod entries %d\n", fAODCTS->GetEntriesFast());
 182 }
 183
 184 //____________________________________________________________________________
 185 void AliCaloTrackESDReader::FillInputEMCAL() {
 186   //Return array with EMCAL clusters in aod format
 187
 188   //Get vertex for momentum calculation
 189   Double_t v[3] ; //vertex ;
 190   GetVertex(v);
 191
 192   Float_t pos[3] ;
 193   Int_t naod      = (fOutputEvent->GetCaloClusters())->GetEntriesFast();
 194   Int_t nTracks   = fInputEvent->GetNumberOfTracks() ;
 195
 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")){
 206
 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());
 209
 210           clus->GetPosition(pos) ;
 211           Int_t id = clus->GetID();
 212           Int_t nLabel = clus->GetNLabels();
 213           Int_t *labels=0x0;
 214           if(clus->GetLabels()) labels =  (clus->GetLabels())->GetArray();
 215
 216           Float_t energy = clus->E();
 217           Char_t ttype= AliAODCluster::kEMCALClusterv1;
 218
 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);
 222
 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()) ;
 230
 231
 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());
 236
 237           caloCluster->SetNCells(clus->GetNCells());
 238           caloCluster->SetCellsAbsId(clus->GetCellsAbsId());
 239           caloCluster->SetCellsAmplitudeFraction(clus->GetCellsAmplitudeFraction());
 240
 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)));
 247             }
 248           }
 249           //Fill reference array
 250         }//Pt and Fidutial cut passed.
 251       }//EMCAL cluster
 252     }//cluster exists
 253   }//esd cluster loop
 254
 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);
 259   }
 260   if(fDebug > 1) printf("AliCaloTrackESDReader::FillInputEMCAL() - aod entries %d\n", fAODEMCAL->GetEntriesFast());
 261
 262 }
 263
 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 ;
 270   GetVertex(v);
 271
 272   Float_t pos[3] ;
 273   Double_t * pid = new Double_t[AliPID::kSPECIESN];
 274   Int_t naod      = (fOutputEvent->GetCaloClusters())->GetEntriesFast();
 275   Int_t nTracks   = fInputEvent->GetNumberOfTracks() ;
 276
 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)) ) {
 281       if (clus->IsPHOS()){
 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")){
 287
 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());
 290
 291           clus->GetPosition(pos) ;
 292           Int_t id = clus->GetID();
 293           Int_t nLabel = clus->GetNLabels();
 294           Int_t *labels=0x0;
 295           if(clus->GetLabels()) labels =  (clus->GetLabels())->GetArray();
 296           Float_t energy = clus->E();
 297
 298           //Phos cluster type
 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;
 308
 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);
 312
 313
 314           caloCluster->SetPIDFromESD(pid);
 315           caloCluster->SetCaloCluster(clus->GetDistanceToBadChannel(), clus->GetClusterDisp(),
 316                                       clus->GetM20(), clus->GetM02(),
 317                                       clus->GetEmcCpvDistance(),  clus->GetNExMax()) ;
 318
 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());
 323
 324           caloCluster->SetNCells(clus->GetNCells());
 325           caloCluster->SetCellsAbsId(clus->GetCellsAbsId());
 326           caloCluster->SetCellsAmplitudeFraction(clus->GetCellsAmplitudeFraction());
 327
 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)));
 334             }
 335           }
 336         }//Pt and Fidutial cut passed.
 337       }//PHOS cluster
 338     }//cluster exists
 339   }//esd cluster loop
 340
 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);
 344     fAODPHOS->Add(clus);
 345   }
 346   if(fDebug > 1) printf("AliCaloTrackESDReader::FillInputPHOS() - aod entries %d\n", fAODPHOS->GetEntriesFast());
 347
 348 }
 349
 350 //____________________________________________________________________________
 351 void AliCaloTrackESDReader::FillInputEMCALCells() {
 352   //Return array with EMCAL cells in esd format
 353
 354   fEMCALCells = (TNamed*) ((AliESDEvent*)fInputEvent)->GetEMCALCells();
 355
 356 }
 357
 358 //____________________________________________________________________________
 359 void AliCaloTrackESDReader::FillInputPHOSCells() {
 360   //Return array with PHOS cells in esd format
 361
 362   fPHOSCells = (TNamed*) ((AliESDEvent*)fInputEvent)->GetPHOSCells();
 363
 364 }
 365
 366 //____________________________________________________________________________
 367 void AliCaloTrackESDReader::GetVertex(Double_t  v[3]) const {
 368   //Return vertex position
 369
 370   //((AliESDEvent*)fInputEvent)->GetVertex()->GetXYZ(v) ;//SPD
 371   ((AliESDEvent*)fInputEvent)->GetPrimaryVertex()->GetXYZ(v);
 372
 373 }
 374
 375
 376 //____________________________________________________________________________
 377 Double_t AliCaloTrackESDReader::GetBField() const {
 378   //Return magnetic field
 379
 380   Double_t bfield = ((AliESDEvent*)fInputEvent)->GetMagneticField();
 381
 382   return bfield;
 383 }
 384
 385
 386 //____________________________________________________________________________
 387 void AliCaloTrackESDReader::SetInputOutputMCEvent(AliVEvent* esd, AliAODEvent* aod, AliMCEvent* mc) {
 388   // Connect the data pointers
 389
 390   if(strcmp(esd->GetName(),"AliESDEvent")){
 391     printf("AliCaloTrackESDReader::SetInputOutputMCEvent() - STOP ::Wrong reader, here only ESDs. Input name: %s != AliESDEvent \n",esd->GetName());
 392     abort();
 393   }
 394
 395   SetInputEvent(esd);
 396   SetOutputEvent(aod);
 397   SetMC(mc);
 398
 399 }