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 "AliESDCaloCluster.h"
  39 #include "AliAODCaloCluster.h"
  40 #include "AliAODTrack.h"
  41 #include "AliAODPid.h"
  42 #include "AliAODEvent.h"
  43 #include "AliFiducialCut.h"
  44
  45 ClassImp(AliCaloTrackESDReader)
  46
  47 //____________________________________________________________________________
  48 AliCaloTrackESDReader::AliCaloTrackESDReader() :
  49 AliCaloTrackReader()
  50 {
  51         //Default Ctor
  52
  53         //Initialize parameters
  54         fDataType=kESD;
  55         fReadStack          = kTRUE;
  56         fReadAODMCParticles = kFALSE;
  57         //We want tracks fitted in the detectors:
  58         fTrackStatus=AliESDtrack::kTPCrefit;
  59     fTrackStatus|=AliESDtrack::kITSrefit;
  60
  61 }
  62
  63 //____________________________________________________________________________
  64 AliCaloTrackESDReader::AliCaloTrackESDReader(const AliCaloTrackESDReader & g) :
  65   AliCaloTrackReader(g)
  66 {
  67   // cpy ctor
  68 }
  69
  70 //_________________________________________________________________________
  71 //AliCaloTrackESDReader & AliCaloTrackESDReader::operator = (const AliCaloTrackESDReader & source)
  72 //{
  73 //      // assignment operator
  74 //
  75 //      if(&source == this) return *this;
  76 //
  77 //      return *this;
  78 //
  79 //}
  80
  81 //____________________________________________________________________________
  82 void AliCaloTrackESDReader::FillInputCTS() {
  83   //Return array with CTS tracks
  84   if(fDebug > 2 ) printf("AliCaloTrackESDReader::FillInputCTS()\n");
  85
  86   //TObjArray * fAODCTS = new TObjArray();
  87   Int_t nTracks   = fInputEvent->GetNumberOfTracks() ;
  88   Int_t naod = 0;
  89   Double_t pos[3];
  90   Double_t p[3];
  91   Double_t covTr[21];
  92   Double_t pid[10];
  93   Double_t bfield = ((AliESDEvent*)fInputEvent)->GetMagneticField();
  94
  95   Double_t      timezero        = 0;   //TO BE FIXED
  96
  97   //To be replaced by call to AliEMCALGeoUtils when the class becomes available
  98   Double_t radius = 441.0; //[cm] EMCAL radius +13cm
  99   if(fDebug > 1) printf("AliCaloTrackESDReader::FillInputCTS() - org entries %d\n", nTracks);
 100   for (Int_t itrack =  0; itrack <  nTracks; itrack++) {////////////// track loop
 101     AliESDtrack * track = (AliESDtrack*) ((AliESDEvent*)fInputEvent)->GetTrack(itrack) ; // retrieve track from esd
 102
 103     //We want tracks fitted in the detectors: TPCrefit, ITSrefit ... check the set bits.
 104         if (fTrackStatus && !((track->GetStatus() & fTrackStatus) == fTrackStatus)) continue ;
 105
 106         track->GetPxPyPz(p) ;
 107         TLorentzVector momentum(p[0],p[1],p[2],0);
 108
 109         if(fCTSPtMin < momentum.Pt() &&fFiducialCut->IsInFiducialCut(momentum,"CTS")){
 110
 111                 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",
 112                                                     momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 113
 114                 track->GetXYZ(pos);
 115                 track->GetCovarianceXYZPxPyPz(covTr);
 116                 track->GetESDpid(pid);
 117
 118                 Float_t impactXY, impactZ;
 119
 120                 track->GetImpactParameters(impactXY,impactZ);
 121
 122                 if (impactXY<3) {
 123                   // track inside the beam pipe
 124                   //Put new aod object in file in AOD tracks array
 125                   AliAODTrack *aodTrack = new((*(fOutputEvent->GetTracks()))[naod++])
 126                     AliAODTrack(track->GetID(), track->GetLabel(), p, kTRUE, pos, kFALSE,covTr, (Short_t)track->GetSign(), track->GetITSClusterMap(),
 127                                 pid,
 128                                 0x0,//primary,
 129                                 kTRUE, // check if this is right
 130                                 kTRUE, // check if this is right
 131                                 AliAODTrack::kPrimary,
 132                                 0);
 133
 134                   aodTrack->SetFlags(track->GetStatus());
 135                   aodTrack->ConvertAliPIDtoAODPID();
 136
 137                   //fill needed AliAODPid information, including
 138                   //Extrapolate track to EMCAL surface for AOD-level track-cluster matching
 139                   AliAODPid *aodpid = new AliAODPid;
 140                   aodpid->SetITSsignal(track->GetITSsignal());
 141                   aodpid->SetTPCsignal(track->GetTPCsignal());
 142                   //n TRD planes = 6
 143                   Int_t nslices = track->GetNumberOfTRDslices()*6;
 144                   Double_t *trdslices = new Double_t[nslices];
 145                   for(Int_t iSl =0; iSl < track->GetNumberOfTRDslices(); iSl++) {
 146                     for(Int_t iPl =0; iPl<6; iPl++) trdslices[iPl*track->GetNumberOfTRDslices()+iSl] = track->GetTRDslice(iPl,iSl);
 147                   }
 148                   aodpid->SetTRDsignal(track->GetNumberOfTRDslices()*6,trdslices);
 149                   Double_t times[AliAODPid::kSPECIES]; track->GetIntegratedTimes(times);
 150                   aodpid->SetIntegratedTimes(times);
 151
 152                   aodpid->SetTOFsignal(track->GetTOFsignal()-timezero); // to be fixed
 153                   aodpid->SetHMPIDsignal(track->GetHMPIDsignal());
 154
 155                   Double_t emcpos[3] = {0.,0.,0.};
 156                   Double_t emcmom[3] = {0.,0.,0.};
 157                   aodpid->SetEMCALPosition(emcpos);
 158                   aodpid->SetEMCALMomentum(emcmom);
 159
 160                   AliExternalTrackParam *outerparam = (AliExternalTrackParam*)track->GetOuterParam();
 161                   if(!outerparam) continue;
 162
 163                   Bool_t okpos = outerparam->GetXYZAt(radius,bfield,emcpos);
 164                   Bool_t okmom = outerparam->GetPxPyPzAt(radius,bfield,emcmom);
 165                   if(!(okpos && okmom)) continue;
 166
 167                   aodpid->SetEMCALPosition(emcpos);
 168                   aodpid->SetEMCALMomentum(emcmom);
 169
 170                   aodTrack->SetDetPID(aodpid);
 171                 }
 172                 else continue;   // outside the beam pipe: orphan track
 173         }//Pt and Fiducial cut passed.
 174   }// track loop
 175
 176   //Put references to selected tracks in array
 177   for(Int_t itrack = 0; itrack < (fOutputEvent->GetTracks())->GetEntriesFast(); itrack++){
 178     AliAODTrack * track =  (AliAODTrack*) (fOutputEvent->GetTracks())->At(itrack);
 179     fAODCTS->Add(track);
 180   }
 181
 182   if(fDebug > 1) printf("AliCaloTrackESDReader::FillInputCTS() - aod entries %d\n", fAODCTS->GetEntriesFast());
 183 }
 184
 185 //____________________________________________________________________________
 186 void AliCaloTrackESDReader::FillInputEMCAL() {
 187   //Return array with EMCAL clusters in aod format
 188   if(fDebug > 2 ) printf("AliCaloTrackESDReader::FillInputEMCAL()\n");
 189
 190   //Get vertex for momentum calculation
 191   Double_t v[3] ; //vertex ;
 192   GetVertex(v);
 193
 194   Float_t pos[3] ;
 195   Int_t naod      = (fOutputEvent->GetCaloClusters())->GetEntriesFast();
 196   Int_t nTracks   = fInputEvent->GetNumberOfTracks() ;
 197
 198   //Loop to select clusters in fiducial cut and fill container with aodClusters
 199   for (Int_t iclus = 0; iclus < ((AliESDEvent*)fInputEvent)->GetNumberOfCaloClusters(); iclus++) {
 200     AliESDCaloCluster * clus = 0;
 201     if ( (clus = ((AliESDEvent*)fInputEvent)->GetCaloCluster(iclus)) ) {
 202       if (clus->IsEMCAL()){
 203
 204         //Check if the cluster contains any bad channel and if close to calorimeter borders
 205         if( GetCaloUtils()->ClusterContainsBadChannel("EMCAL",clus->GetCellsAbsId(), clus->GetNCells())) continue;
 206         if(!GetCaloUtils()->CheckCellFiducialRegion(clus, ((AliESDEvent*)fInputEvent)->GetEMCALCells())) continue;
 207
 208         TLorentzVector momentum ;
 209         clus->GetMomentum(momentum, v);
 210         if(fDebug > 3 && momentum.E() > 0.1) printf("AliCaloTrackESDReader::FillInputEMCAL() - all clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
 211                                                    momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 212         if(fEMCALPtMin < momentum.Pt() && fFiducialCut->IsInFiducialCut(momentum,"EMCAL")){
 213
 214           if(fDebug > 2 && momentum.E() > 0.1) printf("AliCaloTrackESDReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
 215                                                      momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 216
 217           clus->GetPosition(pos) ;
 218           Int_t id = clus->GetID();
 219           Int_t nLabel = clus->GetNLabels();
 220           Int_t *labels=0x0;
 221           if(clus->GetLabels()) labels =  (clus->GetLabels())->GetArray();
 222
 223           Float_t energy = clus->E();
 224           Char_t ttype= AliAODCluster::kEMCALClusterv1;
 225
 226           //Recalibrate the cluster energy
 227                 if(GetCaloUtils()->IsRecalibrationOn()) energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliESDCaloCells*)GetEMCALCells());
 228
 229
 230           //Put new aod object in file in AOD calo clusters array
 231           AliAODCaloCluster *caloCluster = new((*(fOutputEvent->GetCaloClusters()))[naod++])
 232             AliAODCaloCluster(id,nLabel,labels,energy, pos, NULL,ttype,0);
 233
 234           //       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",
 235           //         pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
 236           //         pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]);
 237           caloCluster->SetPIDFromESD(clus->GetPid());
 238           caloCluster->SetCaloCluster(clus->GetDistanceToBadChannel(), clus->GetClusterDisp(),
 239                                       clus->GetM20(), clus->GetM02(),
 240                                       clus->GetEmcCpvDistance(),  clus->GetNExMax(), clus->GetTOF()) ;
 241
 242
 243           if(fDebug > 3 && momentum.E() > 0.1)
 244             printf("AliCaloTrackESDReader::FillInputEMCAL() - Selected clusters Distance BC %2.2f, dispersion %2.2f, M20 %f, M02 %3.2f, NexMax %d, TOF %e\n",
 245                    clus->GetDistanceToBadChannel(), clus->GetClusterDisp(),clus->GetM20(), clus->GetM02(),
 246                    clus->GetNExMax(), clus->GetTOF());
 247
 248           caloCluster->SetNCells(clus->GetNCells());
 249           caloCluster->SetCellsAbsId(clus->GetCellsAbsId());
 250           caloCluster->SetCellsAmplitudeFraction(clus->GetCellsAmplitudeFraction());
 251
 252           TArrayI* matchedT =   clus->GetTracksMatched();
 253           if (nTracks > 0 && matchedT && clus->GetTrackMatched() >= 0) {
 254             for (Int_t im = 0; im < matchedT->GetSize(); im++) {
 255               Int_t iESDtrack = matchedT->At(im);
 256               if(iESDtrack < nTracks && iESDtrack > -1)
 257                 caloCluster->AddTrackMatched((fInputEvent->GetTrack(iESDtrack)));
 258             }
 259           }
 260           //Fill reference array
 261         }//Pt and Fiducial cut passed.
 262       }//EMCAL cluster
 263     }//cluster exists
 264   }//esd cluster loop
 265
 266   //Put references to selected clusters in array
 267   for(Int_t iclus = 0; iclus < (fOutputEvent->GetCaloClusters())->GetEntriesFast(); iclus++){
 268     AliAODCaloCluster * clus =  (AliAODCaloCluster*) (fOutputEvent->GetCaloClusters())->At(iclus);
 269     fAODEMCAL->Add(clus);
 270   }
 271   if(fDebug > 1) printf("AliCaloTrackESDReader::FillInputEMCAL() - aod entries %d\n", fAODEMCAL->GetEntriesFast());
 272
 273 }
 274
 275 //____________________________________________________________________________
 276 void AliCaloTrackESDReader::FillInputPHOS() {
 277   //Return array with PHOS clusters in aod format
 278   if(fDebug > 2 ) printf("AliCaloTrackESDReader::FillInputPHOS()\n");
 279
 280   Int_t nEMCAL = (fOutputEvent->GetCaloClusters())->GetEntriesFast();
 281   //Get vertex for momentum calculation
 282   Double_t v[3] ; //vertex ;
 283   GetVertex(v);
 284
 285   Float_t pos[3] ;
 286   Double_t * pid = new Double_t[AliPID::kSPECIESN];
 287   Int_t naod      = (fOutputEvent->GetCaloClusters())->GetEntriesFast();
 288   Int_t nTracks   = fInputEvent->GetNumberOfTracks() ;
 289
 290   //Loop to select clusters in fiducial cut and fill container with aodClusters
 291   for (Int_t iclus = 0; iclus < ((AliESDEvent*)fInputEvent)->GetNumberOfCaloClusters(); iclus++) {
 292     AliESDCaloCluster * clus = 0;
 293     if ( (clus = ((AliESDEvent*)fInputEvent)->GetCaloCluster(iclus)) ) {
 294       if (clus->IsPHOS()){
 295
 296         //Check if the cluster contains any bad channel and if close to calorimeter borders
 297         if( GetCaloUtils()->ClusterContainsBadChannel("PHOS",clus->GetCellsAbsId(), clus->GetNCells())) continue;
 298         if(!GetCaloUtils()->CheckCellFiducialRegion(clus, ((AliESDEvent*)fInputEvent)->GetPHOSCells())) continue;
 299
 300         TLorentzVector momentum ;
 301         clus->GetMomentum(momentum, v);
 302         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",
 303                                                    momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 304         if(fPHOSPtMin < momentum.Pt() && fFiducialCut->IsInFiducialCut(momentum,"PHOS")){
 305
 306           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",
 307                                                      momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 308
 309           clus->GetPosition(pos) ;
 310           Int_t id = clus->GetID();
 311           Int_t nLabel = clus->GetNLabels();
 312           Int_t *labels=0x0;
 313           if(clus->GetLabels()) labels =  (clus->GetLabels())->GetArray();
 314           Float_t energy = clus->E();
 315
 316           //Phos cluster type
 317           pid = clus->GetPid();
 318           // 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",
 319           //     pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
 320           //     pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]);
 321           Char_t ttype= AliAODCluster::kPHOSNeutral;
 322           Float_t wNeutral = pid[AliPID::kNeutron]+ pid[AliPID::kKaon0]+pid[AliPID::kPhoton]+pid[AliPID::kPi0];
 323           Float_t wCharged = pid[AliPID::kMuon] + pid[AliPID::kElectron] + pid[AliPID::kEleCon]+
 324             pid[AliPID::kProton]+pid[AliPID::kKaon]+pid[AliPID::kPion];
 325           if( wCharged > wNeutral)  ttype= AliAODCluster::kPHOSCharged;
 326
 327           //Recalibrate the cluster energy
 328           if(GetCaloUtils()->IsRecalibrationOn()) energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliESDCaloCells*)GetPHOSCells());
 329
 330           //Put new aod object in file in AOD calo clusters array
 331           AliAODCaloCluster *caloCluster = new((*(fOutputEvent->GetCaloClusters()))[naod++])
 332             AliAODCaloCluster(id,nLabel,labels,energy, pos, NULL, ttype, 0);
 333
 334
 335           caloCluster->SetPIDFromESD(pid);
 336           caloCluster->SetCaloCluster(clus->GetDistanceToBadChannel(), clus->GetClusterDisp(),
 337                                       clus->GetM20(), clus->GetM02(),
 338                                       clus->GetEmcCpvDistance(),  clus->GetNExMax()) ;
 339
 340           if(fDebug > 3 && momentum.E() > 0.2)
 341             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",
 342                    clus->GetDistanceToBadChannel(), clus->GetClusterDisp(),clus->GetM20(), clus->GetM02(),
 343                    clus->GetEmcCpvDistance(),  clus->GetNExMax(), clus->GetTOF());
 344
 345           caloCluster->SetNCells(clus->GetNCells());
 346           caloCluster->SetCellsAbsId(clus->GetCellsAbsId());
 347           caloCluster->SetCellsAmplitudeFraction(clus->GetCellsAmplitudeFraction());
 348
 349           TArrayI* matchedT =   clus->GetTracksMatched();
 350           if (nTracks > 0 && matchedT && clus->GetTrackMatched() >= 0) {
 351             for (Int_t im = 0; im < matchedT->GetSize(); im++) {
 352               Int_t iESDtrack = matchedT->At(im);
 353               if(iESDtrack < nTracks && iESDtrack > -1)
 354                 caloCluster->AddTrackMatched((fInputEvent->GetTrack(iESDtrack)));
 355             }
 356           }
 357         }//Pt and Fiducial cut passed.
 358       }//PHOS cluster
 359     }//cluster exists
 360   }//esd cluster loop
 361
 362   //Put references to selected clusters in array
 363   for(Int_t iclus = nEMCAL; iclus < (fOutputEvent->GetCaloClusters())->GetEntriesFast(); iclus++){
 364     AliAODCaloCluster * clus =  (AliAODCaloCluster*) (fOutputEvent->GetCaloClusters())->At(iclus);
 365     fAODPHOS->Add(clus);
 366   }
 367   if(fDebug > 1) printf("AliCaloTrackESDReader::FillInputPHOS() - aod entries %d\n", fAODPHOS->GetEntriesFast());
 368
 369 }
 370
 371 //____________________________________________________________________________
 372 void AliCaloTrackESDReader::FillInputEMCALCells() {
 373   //Return array with EMCAL cells in esd format
 374
 375   fEMCALCells = (TNamed*) ((AliESDEvent*)fInputEvent)->GetEMCALCells();
 376
 377 }
 378
 379 //____________________________________________________________________________
 380 void AliCaloTrackESDReader::FillInputPHOSCells() {
 381   //Return array with PHOS cells in esd format
 382
 383   fPHOSCells = (TNamed*) ((AliESDEvent*)fInputEvent)->GetPHOSCells();
 384
 385 }
 386
 387
 388 //____________________________________________________________________________
 389 void AliCaloTrackESDReader::GetVertex(Double_t  v[3]) const {
 390   //Return vertex position
 391
 392   //((AliESDEvent*)fInputEvent)->GetVertex()->GetXYZ(v) ;//SPD
 393   ((AliESDEvent*)fInputEvent)->GetPrimaryVertex()->GetXYZ(v);
 394
 395 }
 396
 397
 398 //____________________________________________________________________________
 399 Double_t AliCaloTrackESDReader::GetBField() const {
 400   //Return magnetic field
 401
 402   Double_t bfield = ((AliESDEvent*)fInputEvent)->GetMagneticField();
 403
 404   return bfield;
 405 }
 406
 407
 408 //____________________________________________________________________________
 409 void AliCaloTrackESDReader::SetInputOutputMCEvent(AliVEvent* esd, AliAODEvent* aod, AliMCEvent* mc) {
 410   // Connect the data pointers
 411
 412   if(strcmp(esd->GetName(),"AliESDEvent")){
 413     printf("AliCaloTrackESDReader::SetInputOutputMCEvent() - STOP ::Wrong reader, here only ESDs. Input name: %s != AliESDEvent \n",esd->GetName());
 414     abort();
 415   }
 416
 417   SetInputEvent(esd);
 418   SetOutputEvent(aod);
 419   SetMC(mc);
 420
 421 }