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