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
 205         if(ClusterContainsBadChannel("EMCAL",clus->GetCellsAbsId(), clus->GetNCells())) continue;
 206
 207         TLorentzVector momentum ;
 208         clus->GetMomentum(momentum, v);
 209         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",
 210                                                    momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 211         if(fEMCALPtMin < momentum.Pt() && fFiducialCut->IsInFiducialCut(momentum,"EMCAL")){
 212
 213           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",
 214                                                      momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 215
 216           clus->GetPosition(pos) ;
 217           Int_t id = clus->GetID();
 218           Int_t nLabel = clus->GetNLabels();
 219           Int_t *labels=0x0;
 220           if(clus->GetLabels()) labels =  (clus->GetLabels())->GetArray();
 221
 222           Float_t energy = clus->E();
 223           Char_t ttype= AliAODCluster::kEMCALClusterv1;
 224
 225           //Put new aod object in file in AOD calo clusters array
 226           AliAODCaloCluster *caloCluster = new((*(fOutputEvent->GetCaloClusters()))[naod++])
 227             AliAODCaloCluster(id,nLabel,labels,energy, pos, NULL,ttype,0);
 228
 229           //       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",
 230           //         pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
 231           //         pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]);
 232           caloCluster->SetPIDFromESD(clus->GetPid());
 233           caloCluster->SetCaloCluster(clus->GetDistanceToBadChannel(), clus->GetClusterDisp(),
 234                                       clus->GetM20(), clus->GetM02(),
 235                                       clus->GetEmcCpvDistance(),  clus->GetNExMax(), clus->GetTOF()) ;
 236
 237
 238           if(fDebug > 3 && momentum.E() > 0.1)
 239             printf("AliCaloTrackESDReader::FillInputEMCAL() - Selected clusters Distance BC %2.2f, dispersion %2.2f, M20 %f, M02 %3.2f, NexMax %d, TOF %e\n",
 240                    clus->GetDistanceToBadChannel(), clus->GetClusterDisp(),clus->GetM20(), clus->GetM02(),
 241                    clus->GetNExMax(), clus->GetTOF());
 242
 243           caloCluster->SetNCells(clus->GetNCells());
 244           caloCluster->SetCellsAbsId(clus->GetCellsAbsId());
 245           caloCluster->SetCellsAmplitudeFraction(clus->GetCellsAmplitudeFraction());
 246
 247           TArrayI* matchedT =   clus->GetTracksMatched();
 248           if (nTracks > 0 && matchedT && clus->GetTrackMatched() >= 0) {
 249             for (Int_t im = 0; im < matchedT->GetSize(); im++) {
 250               Int_t iESDtrack = matchedT->At(im);
 251               if(iESDtrack < nTracks && iESDtrack > -1)
 252                 caloCluster->AddTrackMatched((fInputEvent->GetTrack(iESDtrack)));
 253             }
 254           }
 255           //Fill reference array
 256         }//Pt and Fiducial cut passed.
 257       }//EMCAL cluster
 258     }//cluster exists
 259   }//esd cluster loop
 260
 261   //Put references to selected clusters in array
 262   for(Int_t iclus = 0; iclus < (fOutputEvent->GetCaloClusters())->GetEntriesFast(); iclus++){
 263     AliAODCaloCluster * clus =  (AliAODCaloCluster*) (fOutputEvent->GetCaloClusters())->At(iclus);
 264     fAODEMCAL->Add(clus);
 265   }
 266   if(fDebug > 1) printf("AliCaloTrackESDReader::FillInputEMCAL() - aod entries %d\n", fAODEMCAL->GetEntriesFast());
 267
 268 }
 269
 270 //____________________________________________________________________________
 271 void AliCaloTrackESDReader::FillInputPHOS() {
 272   //Return array with PHOS clusters in aod format
 273   if(fDebug > 2 ) printf("AliCaloTrackESDReader::FillInputPHOS()\n");
 274
 275   Int_t nEMCAL = (fOutputEvent->GetCaloClusters())->GetEntriesFast();
 276   //Get vertex for momentum calculation
 277   Double_t v[3] ; //vertex ;
 278   GetVertex(v);
 279
 280   Float_t pos[3] ;
 281   Double_t * pid = new Double_t[AliPID::kSPECIESN];
 282   Int_t naod      = (fOutputEvent->GetCaloClusters())->GetEntriesFast();
 283   Int_t nTracks   = fInputEvent->GetNumberOfTracks() ;
 284
 285   //Loop to select clusters in fiducial cut and fill container with aodClusters
 286   for (Int_t iclus = 0; iclus < ((AliESDEvent*)fInputEvent)->GetNumberOfCaloClusters(); iclus++) {
 287     AliESDCaloCluster * clus = 0;
 288     if ( (clus = ((AliESDEvent*)fInputEvent)->GetCaloCluster(iclus)) ) {
 289       if (clus->IsPHOS()){
 290
 291         //Check if the cluster contains any bad channel
 292         if(ClusterContainsBadChannel("PHOS",clus->GetCellsAbsId(), clus->GetNCells())) continue;
 293
 294         TLorentzVector momentum ;
 295         clus->GetMomentum(momentum, v);
 296         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",
 297                                                    momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 298         if(fPHOSPtMin < momentum.Pt() && fFiducialCut->IsInFiducialCut(momentum,"PHOS")){
 299
 300           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",
 301                                                      momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 302
 303           clus->GetPosition(pos) ;
 304           Int_t id = clus->GetID();
 305           Int_t nLabel = clus->GetNLabels();
 306           Int_t *labels=0x0;
 307           if(clus->GetLabels()) labels =  (clus->GetLabels())->GetArray();
 308           Float_t energy = clus->E();
 309
 310           //Phos cluster type
 311           pid = clus->GetPid();
 312           // 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",
 313           //     pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
 314           //     pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]);
 315           Char_t ttype= AliAODCluster::kPHOSNeutral;
 316           Float_t wNeutral = pid[AliPID::kNeutron]+ pid[AliPID::kKaon0]+pid[AliPID::kPhoton]+pid[AliPID::kPi0];
 317           Float_t wCharged = pid[AliPID::kMuon] + pid[AliPID::kElectron] + pid[AliPID::kEleCon]+
 318             pid[AliPID::kProton]+pid[AliPID::kKaon]+pid[AliPID::kPion];
 319           if( wCharged > wNeutral)  ttype= AliAODCluster::kPHOSCharged;
 320
 321           //Put new aod object in file in AOD calo clusters array
 322           AliAODCaloCluster *caloCluster = new((*(fOutputEvent->GetCaloClusters()))[naod++])
 323             AliAODCaloCluster(id,nLabel,labels,energy, pos, NULL, ttype, 0);
 324
 325
 326           caloCluster->SetPIDFromESD(pid);
 327           caloCluster->SetCaloCluster(clus->GetDistanceToBadChannel(), clus->GetClusterDisp(),
 328                                       clus->GetM20(), clus->GetM02(),
 329                                       clus->GetEmcCpvDistance(),  clus->GetNExMax()) ;
 330
 331           if(fDebug > 3 && momentum.E() > 0.2)
 332             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",
 333                    clus->GetDistanceToBadChannel(), clus->GetClusterDisp(),clus->GetM20(), clus->GetM02(),
 334                    clus->GetEmcCpvDistance(),  clus->GetNExMax(), clus->GetTOF());
 335
 336           caloCluster->SetNCells(clus->GetNCells());
 337           caloCluster->SetCellsAbsId(clus->GetCellsAbsId());
 338           caloCluster->SetCellsAmplitudeFraction(clus->GetCellsAmplitudeFraction());
 339
 340           TArrayI* matchedT =   clus->GetTracksMatched();
 341           if (nTracks > 0 && matchedT && clus->GetTrackMatched() >= 0) {
 342             for (Int_t im = 0; im < matchedT->GetSize(); im++) {
 343               Int_t iESDtrack = matchedT->At(im);
 344               if(iESDtrack < nTracks && iESDtrack > -1)
 345                 caloCluster->AddTrackMatched((fInputEvent->GetTrack(iESDtrack)));
 346             }
 347           }
 348         }//Pt and Fiducial cut passed.
 349       }//PHOS cluster
 350     }//cluster exists
 351   }//esd cluster loop
 352
 353   //Put references to selected clusters in array
 354   for(Int_t iclus = nEMCAL; iclus < (fOutputEvent->GetCaloClusters())->GetEntriesFast(); iclus++){
 355     AliAODCaloCluster * clus =  (AliAODCaloCluster*) (fOutputEvent->GetCaloClusters())->At(iclus);
 356     fAODPHOS->Add(clus);
 357   }
 358   if(fDebug > 1) printf("AliCaloTrackESDReader::FillInputPHOS() - aod entries %d\n", fAODPHOS->GetEntriesFast());
 359
 360 }
 361
 362 //____________________________________________________________________________
 363 void AliCaloTrackESDReader::FillInputEMCALCells() {
 364   //Return array with EMCAL cells in esd format
 365
 366   fEMCALCells = (TNamed*) ((AliESDEvent*)fInputEvent)->GetEMCALCells();
 367
 368 }
 369
 370 //____________________________________________________________________________
 371 void AliCaloTrackESDReader::FillInputPHOSCells() {
 372   //Return array with PHOS cells in esd format
 373
 374   fPHOSCells = (TNamed*) ((AliESDEvent*)fInputEvent)->GetPHOSCells();
 375
 376 }
 377
 378 //____________________________________________________________________________
 379 void AliCaloTrackESDReader::GetVertex(Double_t  v[3]) const {
 380   //Return vertex position
 381
 382   //((AliESDEvent*)fInputEvent)->GetVertex()->GetXYZ(v) ;//SPD
 383   ((AliESDEvent*)fInputEvent)->GetPrimaryVertex()->GetXYZ(v);
 384
 385 }
 386
 387
 388 //____________________________________________________________________________
 389 Double_t AliCaloTrackESDReader::GetBField() const {
 390   //Return magnetic field
 391
 392   Double_t bfield = ((AliESDEvent*)fInputEvent)->GetMagneticField();
 393
 394   return bfield;
 395 }
 396
 397
 398 //____________________________________________________________________________
 399 void AliCaloTrackESDReader::SetInputOutputMCEvent(AliVEvent* esd, AliAODEvent* aod, AliMCEvent* mc) {
 400   // Connect the data pointers
 401
 402   if(strcmp(esd->GetName(),"AliESDEvent")){
 403     printf("AliCaloTrackESDReader::SetInputOutputMCEvent() - STOP ::Wrong reader, here only ESDs. Input name: %s != AliESDEvent \n",esd->GetName());
 404     abort();
 405   }
 406
 407   SetInputEvent(esd);
 408   SetOutputEvent(aod);
 409   SetMC(mc);
 410
 411 }