PWG4/PartCorrBase/AliCaloTrackAODReader.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 (AODs) in order to do prompt gamma
  20 //  or other particle identification and correlations.
  21 // Mixing analysis can be done, input AOD with events
  22 // is opened in the AliCaloTrackReader::Init()
  23 //
  24 //
  25 //*-- Author: Gustavo Conesa (LNF-INFN)
  26 //////////////////////////////////////////////////////////////////////////////
  27
  28
  29 // --- ROOT system ---
  30 //#include "Riostream.h"
  31
  32 //---- ANALYSIS system ----
  33 #include "AliCaloTrackAODReader.h"
  34 #include "AliAODEvent.h"
  35 #include "AliAODCaloCluster.h"
  36 #include "AliAODTrack.h"
  37 #include "AliESDtrack.h"
  38 #include "AliFidutialCut.h"
  39
  40 ClassImp(AliCaloTrackAODReader)
  41
  42 //____________________________________________________________________________
  43 AliCaloTrackAODReader::AliCaloTrackAODReader() :
  44   AliCaloTrackReader(), fWriteOutputAOD(kFALSE)
  45 {
  46   //Default Ctor
  47
  48   //Initialize parameters
  49   fDataType=kAOD;
  50   fReadStack          = kTRUE;
  51   fReadAODMCParticles = kFALSE;
  52   //We want tracks fitted in the detectors:
  53   fTrackStatus=AliESDtrack::kTPCrefit;
  54   fTrackStatus|=AliESDtrack::kITSrefit;
  55
  56 }
  57
  58 //____________________________________________________________________________
  59 AliCaloTrackAODReader::AliCaloTrackAODReader(const AliCaloTrackAODReader & aodr) :
  60   AliCaloTrackReader(aodr),fWriteOutputAOD(aodr.fWriteOutputAOD)
  61 {
  62   // cpy ctor
  63 }
  64
  65 //_________________________________________________________________________
  66 //AliCaloTrackAODReader & AliCaloTrackAODReader::operator = (const AliCaloTrackAODReader & source)
  67 //{
  68 //  // assignment operator
  69 //
  70 //  if(&source == this) return *this;
  71 //
  72 //  return *this;
  73 //
  74 //}
  75
  76 //____________________________________________________________________________
  77 void AliCaloTrackAODReader::FillInputCTS() {
  78   //Return array with Central Tracking System (CTS) tracks
  79
  80   Int_t nTracks   = fInputEvent->GetNumberOfTracks() ;
  81   Int_t naod = 0;
  82   Double_t p[3];
  83   for (Int_t itrack =  0; itrack <  nTracks; itrack++) {////////////// track loop
  84     AliAODTrack * track = ((AliAODEvent*)fInputEvent)->GetTrack(itrack) ; // retrieve track from esd
  85
  86     //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
  87         if (fTrackStatus && !((track->GetStatus() & fTrackStatus) == fTrackStatus)) continue ;
  88
  89     track->GetPxPyPz(p) ;
  90     TLorentzVector momentum(p[0],p[1],p[2],0);
  91
  92     if(fCTSPtMin < momentum.Pt() && fFidutialCut->IsInFidutialCut(momentum,"CTS")){
  93
  94       if(fDebug > 2 && momentum.Pt() > 0.1) printf("AliCaloTrackAODReader::FillInputCTS() - Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
  95                                                   momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
  96
  97           if(fWriteOutputAOD){
  98                    AliAODTrack* newtrack = new((*(fOutputEvent->GetTracks()))[naod++]) AliAODTrack(*track);
  99                    fAODCTS->Add(newtrack); //Use AOD stored in output for references.
 100           }
 101           else fAODCTS->Add(track);
 102     }//Pt and Fidutial cut passed.
 103   }// track loop
 104
 105   fAODCTSNormalInputEntries = fAODCTS->GetEntriesFast();
 106   if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputCTS()   - aod entries %d\n", fAODCTSNormalInputEntries);
 107
 108   //If second input event available, add the clusters.
 109   if(fSecondInputAODTree && fSecondInputAODEvent){
 110           nTracks   = fSecondInputAODEvent->GetNumberOfTracks() ;
 111           if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputCTS()   - Add second input tracks, entries %d\n", nTracks) ;
 112           for (Int_t itrack =  0; itrack <  nTracks; itrack++) {////////////// track loop
 113                   AliAODTrack * track = ((AliAODEvent*)fSecondInputAODEvent)->GetTrack(itrack) ; // retrieve track from esd
 114
 115                   //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
 116                   if (fTrackStatus && !((track->GetStatus() & fTrackStatus) == fTrackStatus)) continue ;
 117
 118                   track->GetPxPyPz(p) ;
 119                   TLorentzVector momentum(p[0],p[1],p[2],0);
 120
 121                   if(fCTSPtMin < momentum.Pt() && fFidutialCut->IsInFidutialCut(momentum,"CTS")){
 122
 123                           if(fDebug > 2 && momentum.Pt() > 0.1) printf("AliCaloTrackAODReader::FillInputCTS() - Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
 124                                                                                                                    momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 125
 126                           if(fWriteOutputAOD){
 127                                   AliAODTrack* newtrack = new((*(fOutputEvent->GetTracks()))[naod++]) AliAODTrack(*track);
 128                                   fAODCTS->Add(newtrack); //Use AOD stored in output for references.
 129                           }
 130                           else fAODCTS->Add(track);
 131
 132                   }//Pt and Fidutial cut passed.
 133           }// track loop
 134
 135           if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputCTS()   - aod normal entries %d, after second input %d\n", fAODCTSNormalInputEntries, fAODCTS->GetEntriesFast());
 136   }     //second input loop
 137
 138 }
 139
 140 //____________________________________________________________________________
 141 void AliCaloTrackAODReader::FillInputEMCAL() {
 142   //Return array with EMCAL clusters in aod format
 143
 144   //Get vertex for momentum calculation
 145   Double_t v[3] ; //vertex ;
 146   GetVertex(v);
 147
 148   Int_t naod =  (fOutputEvent->GetCaloClusters())->GetEntriesFast();
 149   //Loop to select clusters in fidutial cut and fill container with aodClusters
 150   Int_t nclusters = ((AliAODEvent*)fInputEvent)->GetNCaloClusters();
 151   for (Int_t iclus =  0; iclus <  nclusters; iclus++) {
 152     AliAODCaloCluster * clus = 0;
 153     if ( (clus = ((AliAODEvent*)fInputEvent)->GetCaloCluster(iclus)) ) {
 154       if (clus->IsEMCALCluster()){
 155         TLorentzVector momentum ;
 156         clus->GetMomentum(momentum, v);
 157
 158         if(fEMCALPtMin < momentum.Pt() && fFidutialCut->IsInFidutialCut(momentum,"EMCAL")){
 159
 160           if(fDebug > 2 && momentum.E() > 0.1) printf("AliCaloTrackAODReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
 161                                                      momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 162
 163           if(fWriteOutputAOD){
 164                 AliAODCaloCluster * newclus = new((*(fOutputEvent->GetCaloClusters()))[naod++])AliAODCaloCluster(*clus);
 165                 fAODEMCAL->Add(newclus);
 166           }
 167           else fAODEMCAL->Add(clus);
 168         }//Pt and Fidutial cut passed.
 169       }//EMCAL cluster
 170     }// cluster exists
 171   }// cluster loop
 172
 173   fAODEMCALNormalInputEntries = fAODEMCAL->GetEntriesFast();
 174   if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputEMCAL() - aod entries %d\n", fAODEMCALNormalInputEntries);
 175
 176   //If second input event available, add the clusters.
 177   if(fSecondInputAODTree && fSecondInputAODEvent){
 178           GetSecondInputAODVertex(v);
 179           nclusters = ((AliAODEvent*)fSecondInputAODEvent)->GetNCaloClusters();
 180           if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputEMCAL() - Add second input clusters, entries %d\n", nclusters) ;
 181                 for (Int_t iclus =  0; iclus < nclusters; iclus++) {
 182                         AliAODCaloCluster * clus = 0;
 183                         if ( (clus = ((AliAODEvent*)fSecondInputAODEvent)->GetCaloCluster(iclus)) ) {
 184                                 if (clus->IsEMCALCluster()){
 185                                         TLorentzVector momentum ;
 186                                         clus->GetMomentum(momentum, v);
 187
 188                                         if(fEMCALPtMin < momentum.Pt() && fFidutialCut->IsInFidutialCut(momentum,"EMCAL")){
 189
 190                                                 if(fDebug > 2 && momentum.E() > 0.1) printf("AliCaloTrackAODReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
 191                                                                                                                                         momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 192                                                 if(fWriteOutputAOD){
 193                                                         AliAODCaloCluster * newclus = new((*(fOutputEvent->GetCaloClusters()))[naod++])AliAODCaloCluster(*clus);
 194                                                     fAODEMCAL->Add(newclus);
 195                                                 }
 196                                                 else fAODEMCAL->Add(clus);
 197                                         }//Pt and Fidutial cut passed.
 198                                 }//EMCAL cluster
 199                         }// cluster exists
 200                 }// cluster loop
 201
 202           if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputEMCAL() - aod normal entries %d, after second input %d\n", fAODEMCALNormalInputEntries, fAODEMCAL->GetEntriesFast());
 203
 204         } //second input loop
 205 }
 206
 207 //____________________________________________________________________________
 208 void AliCaloTrackAODReader::FillInputPHOS() {
 209   //Return array with PHOS clusters in aod format
 210
 211   //Get vertex for momentum calculation
 212   Double_t v[3] ; //vertex ;
 213   GetVertex(v);
 214
 215   Int_t naod =  (fOutputEvent->GetCaloClusters())->GetEntriesFast();
 216   //Loop to select clusters in fidutial cut and fill container with aodClusters
 217   Int_t nclusters = ((AliAODEvent*)fInputEvent)->GetNCaloClusters();
 218   for (Int_t iclus =  0; iclus < nclusters; iclus++) {
 219     AliAODCaloCluster * clus = 0;
 220     if ( (clus = ((AliAODEvent*)fInputEvent)->GetCaloCluster(iclus)) ) {
 221       if (clus->IsPHOSCluster()){
 222         TLorentzVector momentum ;
 223         clus->GetMomentum(momentum, v);
 224
 225         if(fPHOSPtMin < momentum.Pt() && fFidutialCut->IsInFidutialCut(momentum,"PHOS")){
 226
 227           if(fDebug > 2 && momentum.E() > 0.1) printf("AliCaloTrackAODReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
 228                                                      momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 229
 230                 if(fWriteOutputAOD){
 231                         AliAODCaloCluster * newclus = new((*(fOutputEvent->GetCaloClusters()))[naod++])AliAODCaloCluster(*clus);
 232                         fAODPHOS->Add(newclus);
 233                 }
 234                 else fAODPHOS->Add(clus);
 235         }//Pt and Fidutial cut passed.
 236       }//PHOS cluster
 237     }//cluster exists
 238   }//esd cluster loop
 239
 240   fAODPHOSNormalInputEntries = fAODPHOS->GetEntriesFast() ;
 241   if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputPHOS()  - aod entries %d\n", fAODPHOSNormalInputEntries);
 242
 243   //If second input event available, add the clusters.
 244   if(fSecondInputAODTree && fSecondInputAODEvent){
 245           GetSecondInputAODVertex(v);
 246           nclusters = ((AliAODEvent*)fSecondInputAODEvent)->GetNCaloClusters();
 247           if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputPHOS()  - Add second input clusters, entries %d\n", nclusters);
 248                 for (Int_t iclus =  0; iclus < nclusters; iclus++) {
 249                         AliAODCaloCluster * clus = 0;
 250                         if ( (clus = ((AliAODEvent*)fSecondInputAODEvent)->GetCaloCluster(iclus)) ) {
 251                                 if (clus->IsPHOSCluster()){
 252                                         TLorentzVector momentum ;
 253                                         clus->GetMomentum(momentum, v);
 254
 255                                         if(fPHOSPtMin < momentum.Pt() && fFidutialCut->IsInFidutialCut(momentum,"PHOS")){
 256
 257                                                 if(fDebug > 2 && momentum.E() > 0.1) printf("AliCaloTrackAODReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
 258                                                                                                                                         momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 259                                                 if(fWriteOutputAOD){
 260                                                         AliAODCaloCluster * newclus = new((*(fOutputEvent->GetCaloClusters()))[naod++])AliAODCaloCluster(*clus);
 261                                                     fAODPHOS->Add(newclus);
 262                                                 }
 263                                                 else fAODPHOS->Add(clus);
 264                                         }//Pt and Fidutial cut passed.
 265                                 }//PHOS cluster
 266                         }// cluster exists
 267                 }// cluster loop
 268                 if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputPHOS()  - aod normal entries %d, after second input %d\n", fAODPHOSNormalInputEntries, fAODPHOS->GetEntriesFast());
 269   }     //second input loop
 270
 271 }
 272
 273 //____________________________________________________________________________
 274 void AliCaloTrackAODReader::FillInputEMCALCells() {
 275   //Return array with EMCAL cells in aod format
 276
 277   fEMCALCells = (TNamed*) ((AliAODEvent*)fInputEvent)->GetEMCALCells();
 278
 279 }
 280
 281 //____________________________________________________________________________
 282 void AliCaloTrackAODReader::FillInputPHOSCells() {
 283   //Return array with PHOS cells in aod format
 284
 285   fPHOSCells = (TNamed*) ((AliAODEvent*)fInputEvent)->GetPHOSCells();
 286
 287 }
 288
 289 //____________________________________________________________________________
 290 void AliCaloTrackAODReader::GetVertex(Double_t  v[3]) const {
 291   //Return vertex position
 292
 293   v[0] = ((AliAODEvent*)fInputEvent)->GetVertex(0)->GetX() ;//CHECK!!!
 294   v[1] = ((AliAODEvent*)fInputEvent)->GetVertex(0)->GetY() ;//CHECK!!!
 295   v[2] = ((AliAODEvent*)fInputEvent)->GetVertex(0)->GetZ() ;//CHECK!!!
 296 }
 297
 298 //____________________________________________________________________________
 299 void AliCaloTrackAODReader::GetSecondInputAODVertex(Double_t  v[3]) const {
 300         //Return vertex position of second AOD input
 301
 302         v[0] = ((AliAODEvent*)fSecondInputAODEvent)->GetVertex(0)->GetX() ;//CHECK!!!
 303         v[1] = ((AliAODEvent*)fSecondInputAODEvent)->GetVertex(0)->GetY() ;//CHECK!!!
 304         v[2] = ((AliAODEvent*)fSecondInputAODEvent)->GetVertex(0)->GetZ() ;//CHECK!!!
 305 }
 306
 307 //____________________________________________________________________________
 308 Double_t AliCaloTrackAODReader::GetBField() const {
 309   //Return magnetic field
 310
 311   Double_t bfield =  ((AliAODEvent*)fInputEvent)->GetMagneticField();
 312
 313   return bfield;
 314
 315 }
 316
 317 //____________________________________________________________________________
 318 void AliCaloTrackAODReader::SetInputOutputMCEvent(AliVEvent* input, AliAODEvent* aod, AliMCEvent* mc) {
 319   // Connect the data pointers
 320   // If input is AOD, do analysis with input, if not, do analysis with the output aod.
 321   if(!strcmp(input->GetName(),"AliESDEvent"))   {
 322     SetInputEvent(aod);
 323     SetOutputEvent(aod);
 324   }
 325   else if(!strcmp(input->GetName(),"AliAODEvent")){
 326     SetInputEvent(input);
 327     SetOutputEvent(aod);
 328   }
 329   else{
 330     printf("AliCaloTrackAODReader::SetInputOutputMCEvent() - STOP : Wrong data format: %s\n",input->GetName());
 331     abort();
 332   }
 333
 334   SetMC(mc);
 335
 336 }