PWG4/PartCorrBase/AliCaloTrackReader.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 /* $Id:  $ */
  16
  17 //_________________________________________________________________________
  18 // Base class for reading data: MonteCarlo, ESD or AOD, of PHOS EMCAL and
  19 // Central Barrel Tracking detectors (CTS).
  20 // Not all MC particles/tracks/clusters are kept, some kinematical/fidutial restrictions are done.
  21 // Mother class of : AliCaloTrackESDReader: Fills ESD data in 3 TObjArrays (PHOS, EMCAL, CTS)
  22 //                 : AliCaloTrackMCReader: Fills Kinematics data in 3 TObjArrays (PHOS, EMCAL, CTS)
  23 //                 : AliCaloTrackAODReader: Fills AOD data in 3 TObjArrays (PHOS, EMCAL, CTS)
  24 //
  25 //-- Author: Gustavo Conesa (LNF-INFN)
  26 //////////////////////////////////////////////////////////////////////////////
  27
  28
  29 // --- ROOT system ---
  30 #include "TFile.h"
  31
  32 //---- ANALYSIS system ----
  33 #include "AliCaloTrackReader.h"
  34 #include "AliFidutialCut.h"
  35 #include "AliMCEvent.h"
  36 #include "AliAODMCHeader.h"
  37 #include "AliGenPythiaEventHeader.h"
  38 #include "AliAODEvent.h"
  39
  40 ClassImp(AliCaloTrackReader)
  41
  42
  43 //____________________________________________________________________________
  44   AliCaloTrackReader::AliCaloTrackReader() :
  45     TObject(), fEventNumber(-1), fCurrentFileName(""),fDataType(0), fDebug(0),
  46     fFidutialCut(0x0), fComparePtHardAndJetPt(kFALSE), fPtHardAndJetPtFactor(7),
  47     fCTSPtMin(0), fEMCALPtMin(0),fPHOSPtMin(0),
  48     fAODCTS(new TObjArray()), fAODEMCAL(new TObjArray()), fAODPHOS(new TObjArray()),
  49     fEMCALCells(0x0), fPHOSCells(0x0),
  50     fInputEvent(0x0), fOutputEvent(0x0),fMC(0x0),
  51     fFillCTS(0),fFillEMCAL(0),fFillPHOS(0),
  52     fFillEMCALCells(0),fFillPHOSCells(0),
  53         fSecondInputAODTree(0x0), fSecondInputAODEvent(0x0),
  54         fSecondInputFileName(""),fSecondInputFirstEvent(0),
  55         fAODCTSNormalInputEntries(0), fAODEMCALNormalInputEntries(0),
  56     fAODPHOSNormalInputEntries(0), fTrackStatus(0),
  57         fReadStack(kFALSE), fReadAODMCParticles(kFALSE),
  58         fWriteOutputStdAOD(kFALSE), fDeltaAODFileName("deltaAODPartCorr.root")
  59 {
  60   //Ctor
  61
  62   //Initialize parameters
  63   InitParameters();
  64 }
  65
  66 //____________________________________________________________________________
  67 AliCaloTrackReader::AliCaloTrackReader(const AliCaloTrackReader & g) :
  68   TObject(g), fEventNumber(g.fEventNumber), fCurrentFileName(g.fCurrentFileName),
  69   fDataType(g.fDataType), fDebug(g.fDebug),
  70   fFidutialCut(g.fFidutialCut),
  71   fComparePtHardAndJetPt(g.fComparePtHardAndJetPt),
  72   fPtHardAndJetPtFactor(g.fPtHardAndJetPtFactor),
  73   fCTSPtMin(g.fCTSPtMin), fEMCALPtMin(g.fEMCALPtMin),fPHOSPtMin(g.fPHOSPtMin),
  74   fAODCTS(new TObjArray(*g.fAODCTS)),
  75   fAODEMCAL(new TObjArray(*g.fAODEMCAL)),
  76   fAODPHOS(new TObjArray(*g.fAODPHOS)),
  77   fEMCALCells(new TNamed(*g.fEMCALCells)),
  78   fPHOSCells(new TNamed(*g.fPHOSCells)),
  79   fInputEvent(g.fInputEvent), fOutputEvent(g.fOutputEvent), fMC(g.fMC),
  80   fFillCTS(g.fFillCTS),fFillEMCAL(g.fFillEMCAL),fFillPHOS(g.fFillPHOS),
  81   fFillEMCALCells(g.fFillEMCALCells),fFillPHOSCells(g.fFillPHOSCells),
  82   fSecondInputAODTree(g.fSecondInputAODTree),
  83   fSecondInputAODEvent(g.fSecondInputAODEvent),
  84   fSecondInputFileName(g.fSecondInputFileName),
  85   fSecondInputFirstEvent(g.fSecondInputFirstEvent),
  86   fAODCTSNormalInputEntries(g.fAODCTSNormalInputEntries),
  87   fAODEMCALNormalInputEntries(g.fAODEMCALNormalInputEntries),
  88   fAODPHOSNormalInputEntries(g.fAODPHOSNormalInputEntries),
  89   fTrackStatus(g.fTrackStatus),
  90   fReadStack(g.fReadStack), fReadAODMCParticles(g.fReadAODMCParticles),
  91   fWriteOutputStdAOD(g.fWriteOutputStdAOD), fDeltaAODFileName(g.fDeltaAODFileName)
  92 {
  93   // cpy ctor
  94
  95 }
  96
  97 //_________________________________________________________________________
  98 AliCaloTrackReader & AliCaloTrackReader::operator = (const AliCaloTrackReader & source)
  99 {
 100   // assignment operator
 101
 102   if(&source == this) return *this;
 103
 104   fDataType    = source.fDataType ;
 105   fDebug       = source.fDebug ;
 106   fEventNumber = source.fEventNumber ;
 107   fCurrentFileName = source.fCurrentFileName ;
 108   fFidutialCut = source.fFidutialCut;
 109
 110   fComparePtHardAndJetPt = source.fComparePtHardAndJetPt;
 111   fPtHardAndJetPtFactor  = source.fPtHardAndJetPtFactor;
 112
 113   fCTSPtMin    = source.fCTSPtMin ;
 114   fEMCALPtMin  = source.fEMCALPtMin ;
 115   fPHOSPtMin   = source.fPHOSPtMin ;
 116
 117   fAODCTS     = new TObjArray(*source.fAODCTS) ;
 118   fAODEMCAL   = new TObjArray(*source.fAODEMCAL) ;
 119   fAODPHOS    = new TObjArray(*source.fAODPHOS) ;
 120   fEMCALCells = new TNamed(*source.fEMCALCells) ;
 121   fPHOSCells  = new TNamed(*source.fPHOSCells) ;
 122
 123   fInputEvent  = source.fInputEvent;
 124   fOutputEvent = source.fOutputEvent;
 125   fMC          = source.fMC;
 126
 127   fFillCTS        = source.fFillCTS;
 128   fFillEMCAL      = source.fFillEMCAL;
 129   fFillPHOS       = source.fFillPHOS;
 130   fFillEMCALCells = source.fFillEMCALCells;
 131   fFillPHOSCells  = source.fFillPHOSCells;
 132
 133   fSecondInputAODTree    = source.fSecondInputAODTree;
 134   fSecondInputAODEvent   = source.fSecondInputAODEvent;
 135   fSecondInputFileName   = source.fSecondInputFileName;
 136   fSecondInputFirstEvent = source.fSecondInputFirstEvent;
 137
 138   fAODCTSNormalInputEntries   = source.fAODCTSNormalInputEntries;
 139   fAODEMCALNormalInputEntries = source.fAODEMCALNormalInputEntries;
 140   fAODPHOSNormalInputEntries  = source.fAODPHOSNormalInputEntries;
 141
 142   fTrackStatus        = source.fTrackStatus;
 143   fReadStack          = source.fReadStack;
 144   fReadAODMCParticles = source.fReadAODMCParticles;
 145
 146   fWriteOutputStdAOD  = source.fWriteOutputStdAOD;
 147   fDeltaAODFileName   = source.fDeltaAODFileName;
 148
 149   return *this;
 150
 151 }
 152
 153 //_________________________________
 154 AliCaloTrackReader::~AliCaloTrackReader() {
 155   //Dtor
 156
 157   if(fFidutialCut) delete fFidutialCut ;
 158
 159   if(fAODCTS){
 160     fAODCTS->Clear() ;
 161     delete fAODCTS ;
 162   }
 163
 164   if(fAODEMCAL){
 165     fAODEMCAL->Clear() ;
 166     delete fAODEMCAL ;
 167   }
 168
 169   if(fAODPHOS){
 170     fAODPHOS->Clear() ;
 171     delete fAODPHOS ;
 172   }
 173
 174   if(fEMCALCells){
 175     fEMCALCells->Clear() ;
 176     delete fEMCALCells ;
 177   }
 178
 179   if(fPHOSCells){
 180     fPHOSCells->Clear() ;
 181     delete fPHOSCells ;
 182   }
 183
 184   if(fInputEvent)  delete fInputEvent ;
 185   if(fOutputEvent) delete fOutputEvent ;
 186   if(fMC)          delete fMC ;
 187
 188   if(fSecondInputAODTree){
 189         fSecondInputAODTree->Clear();
 190         delete fSecondInputAODTree;
 191   }
 192
 193   if(fSecondInputAODEvent) delete fSecondInputAODEvent ;
 194
 195 }
 196
 197
 198 //_________________________________________________________________________
 199 Bool_t AliCaloTrackReader::ComparePtHardAndJetPt(){
 200         // Check the event, if the requested ptHard is much larger than the jet pT, then there is a problem.
 201         // Only for PYTHIA.
 202         if(!fReadStack) return kTRUE; //Information not filtered to AOD
 203
 204         if(!strcmp(GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader")){
 205                 TParticle * jet =  new TParticle;
 206                 AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetGenEventHeader();
 207                 Int_t nTriggerJets =  pygeh->NTriggerJets();
 208                 Float_t ptHard = pygeh->GetPtHard();
 209
 210                 //if(fDebug > 1) printf("AliMCAnalysisUtils::PythiaEventHeader: Njets: %d, pT Hard %f\n",nTriggerJets, ptHard);
 211             Float_t tmpjet[]={0,0,0,0};
 212                 for(Int_t ijet = 0; ijet< nTriggerJets; ijet++){
 213                         pygeh->TriggerJet(ijet, tmpjet);
 214                         jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
 215                         //Compare jet pT and pt Hard
 216                         //if(fDebug > 1) printf("AliMCAnalysisUtils:: %d pycell jet pT %f\n",ijet, jet->Pt());
 217                         if(jet->Pt() > fPtHardAndJetPtFactor * ptHard) {
 218                                 printf("AliMCAnalysisUtils::PythiaEventHeader: Njets: %d, pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n",
 219                                            nTriggerJets, ptHard, jet->Pt(), fPtHardAndJetPtFactor);
 220                                 return kFALSE;
 221                         }
 222                 }
 223         }
 224
 225         return kTRUE ;
 226
 227 }
 228
 229 //____________________________________________________________________________
 230 AliStack* AliCaloTrackReader::GetStack() const {
 231   //Return pointer to stack
 232   if(fMC)
 233     return fMC->Stack();
 234   else{
 235     if(fDebug > 1) printf("AliCaloTrackReader::GetStack() - Stack is not available\n");
 236     return 0x0 ;
 237   }
 238 }
 239
 240 //____________________________________________________________________________
 241 AliHeader* AliCaloTrackReader::GetHeader() const {
 242   //Return pointer to header
 243   if(fMC)
 244     return fMC->Header();
 245   else{
 246     printf("AliCaloTrackReader::Header is not available\n");
 247     return 0x0 ;
 248   }
 249 }
 250 //____________________________________________________________________________
 251 AliGenEventHeader* AliCaloTrackReader::GetGenEventHeader() const {
 252   //Return pointer to Generated event header
 253   if(fMC)
 254     return fMC->GenEventHeader();
 255   else{
 256     printf("AliCaloTrackReader::GenEventHeader is not available\n");
 257     return 0x0 ;
 258   }
 259 }
 260
 261 //____________________________________________________________________________
 262 TClonesArray* AliCaloTrackReader::GetAODMCParticles(Int_t input) const {
 263         //Return list of particles in AOD. Do it for the corresponding input event.
 264         if(fDataType == kAOD){
 265          //Normal input AOD
 266          if(input == 0) return (TClonesArray*)((AliAODEvent*)fInputEvent)->FindListObject("mcparticles");
 267           //Second input AOD
 268          else if(input == 1 && fSecondInputAODEvent) return (TClonesArray*) fSecondInputAODEvent->FindListObject("mcparticles");
 269          else {
 270              printf("AliCaloTrackReader::GetAODMCParticles() - wrong AOD input index? %d, or non existing tree? \n",input);
 271                  return 0x0;
 272          }
 273         }
 274         else {
 275                 printf("AliCaloTrackReader::GetAODMCParticles() - Input are not AODs\n");
 276                 return 0x0;
 277         }
 278 }
 279
 280 //____________________________________________________________________________
 281 AliAODMCHeader* AliCaloTrackReader::GetAODMCHeader(Int_t input) const {
 282         //Return MC header in AOD. Do it for the corresponding input event.
 283         if(fDataType == kAOD){
 284                 //Normal input AOD
 285                 if(input == 0) return (AliAODMCHeader*)((AliAODEvent*)fInputEvent)->FindListObject("mcheader");
 286                 //Second input AOD
 287                 else if(input == 1) return  (AliAODMCHeader*) fSecondInputAODEvent->FindListObject("mcheader");
 288                 else {
 289                         printf("AliCaloTrackReader::GetAODMCHeader() - wrong AOD input index, %d\n",input);
 290                         return 0x0;
 291                 }
 292         }
 293         else {
 294                 printf("AliCaloTrackReader::GetAODMCHeader() - Input are not AODs\n");
 295                 return 0x0;
 296         }
 297 }
 298
 299 //_______________________________________________________________
 300 void AliCaloTrackReader::Init()
 301 {
 302         //Init reader. Method to be called in AliAnaPartCorrMaker
 303
 304         //Get the file with second input events if the filename is given
 305         //Get the tree and connect the AODEvent. Only with AODs
 306
 307         if(fReadStack && fReadAODMCParticles){
 308                 printf("AliCaloTrackReader::Init() - Cannot access stack and mcparticles at the same time, change them \n");
 309                 fReadStack = kFALSE;
 310                 fReadAODMCParticles = kFALSE;
 311         }
 312
 313         if(fSecondInputFileName!=""){
 314                 if(fDataType == kAOD){
 315                         TFile * input2   = new TFile(fSecondInputFileName,"read");
 316                         printf("AliCaloTrackReader::Init() - Second input file opened: %s, size %d \n", input2->GetName(), (Int_t) input2->GetSize());
 317                         fSecondInputAODTree = (TTree*) input2->Get("aodTree");
 318                         if(fSecondInputAODTree) printf("AliCaloTrackReader::Init() - Second input tree opened: %s, entries %d \n",
 319                                                                                    fSecondInputAODTree->GetName(), (Int_t) fSecondInputAODTree->GetEntries());
 320                         else{
 321                          printf("AliCaloTrackReader::Init() - Second input tree not available, STOP \n");
 322                          abort();
 323                         }
 324                         fSecondInputAODEvent = new AliAODEvent;
 325                         fSecondInputAODEvent->ReadFromTree(fSecondInputAODTree);
 326                         if(fSecondInputFirstEvent >= fSecondInputAODTree->GetEntriesFast()){
 327                                 printf("AliCaloTrackReader::Init() - Requested first event of second input %d, is larger than number of events %d, STOP\n",
 328                                            fSecondInputFirstEvent, (Int_t) fSecondInputAODTree->GetEntriesFast());
 329                                 abort();
 330                         }
 331                 }
 332                 else printf("AliCaloTrackReader::Init() - Second input not added, reader is not AOD\n");
 333         }
 334
 335
 336 }
 337 //_______________________________________________________________
 338 void AliCaloTrackReader::InitParameters()
 339 {
 340
 341   //Initialize the parameters of the analysis.
 342   fDataType   = kESD ;
 343   fCTSPtMin   = 0.2 ;
 344   fEMCALPtMin = 0.2 ;
 345   fPHOSPtMin  = 0.2 ;
 346
 347   fFillEMCAL      = kTRUE;
 348   fFillPHOS       = kTRUE;
 349   fFillCTS        = kTRUE;
 350   fFillEMCALCells = kFALSE;
 351   fFillPHOSCells  = kFALSE;
 352
 353   fFidutialCut           = new AliFidutialCut();
 354   fSecondInputFileName   = "" ;
 355   fSecondInputFirstEvent = 0 ;
 356   fReadStack             = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
 357   fReadAODMCParticles    = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
 358   fWriteOutputStdAOD     = kTRUE;  // Leave it kTRUE in case we run different tasks at the same time after the filter, when running alone better switch it off.
 359   fDeltaAODFileName      = "deltaAODPartCorr.root";
 360 }
 361
 362
 363 //________________________________________________________________
 364 void AliCaloTrackReader::Print(const Option_t * opt) const
 365 {
 366
 367   //Print some relevant parameters set for the analysis
 368   if(! opt)
 369     return;
 370
 371   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
 372   printf("Data type      : %d\n", fDataType) ;
 373   printf("CTS Min pT     : %2.1f GeV/c\n", fCTSPtMin) ;
 374   printf("EMCAL Min pT   : %2.1f GeV/c\n", fEMCALPtMin) ;
 375   printf("PHOS Min pT    : %2.1f GeV/c\n", fPHOSPtMin) ;
 376   printf("Use CTS         =     %d\n", fFillCTS) ;
 377   printf("Use EMCAL       =     %d\n", fFillEMCAL) ;
 378   printf("Use PHOS        =     %d\n", fFillPHOS) ;
 379   printf("Use EMCAL Cells =     %d\n", fFillEMCALCells) ;
 380   printf("Use PHOS  Cells =     %d\n", fFillPHOSCells) ;
 381   printf("Track status    =     %d\n", (Int_t) fTrackStatus) ;
 382   if(fComparePtHardAndJetPt)
 383           printf("Compare jet pt and pt hard to accept event, factor = %2.2f",fPtHardAndJetPtFactor);
 384
 385   if(fSecondInputFileName!="") {
 386           printf("Second Input File Name     =     %s\n", fSecondInputFileName.Data()) ;
 387           printf("Second Input First Event   =     %d\n", fSecondInputFirstEvent) ;
 388   }
 389
 390   printf("Read Kine from, stack? %d, AOD ? %d \n", fReadStack, fReadAODMCParticles) ;
 391   printf("Write in std AOD    =     %d\n", fWriteOutputStdAOD) ;
 392   printf("Delta AOD File Name =     %s\n", fDeltaAODFileName.Data()) ;
 393
 394   printf("    \n") ;
 395 }
 396
 397 //___________________________________________________
 398 Bool_t AliCaloTrackReader::FillInputEvent(const Int_t iEntry, const char * currentFileName) {
 399   //Fill the event counter and input lists that are needed, called by the analysis maker.
 400
 401   fEventNumber = iEntry;
 402   fCurrentFileName = TString(currentFileName);
 403
 404   if((fDataType != kAOD) && ((fOutputEvent->GetCaloClusters())->GetEntriesFast()!=0 ||(fOutputEvent->GetTracks())->GetEntriesFast()!=0)){
 405     printf("AliCaloTrackReader::AODCaloClusters or AODTracks already filled by the filter, do not use the ESD reader, use the AOD reader, STOP\n");
 406     abort();
 407   }
 408
 409   //In case of analysis of events with jets, skip those with jet pt > 5 pt hard
 410   if(fComparePtHardAndJetPt && GetStack()) {
 411                 if(!ComparePtHardAndJetPt()) return kFALSE ;
 412   }
 413
 414   //In case of mixing events with other AOD file
 415   if(fDataType == kAOD && fSecondInputAODTree){
 416
 417          if(fDebug > 1)
 418                  printf("AliCaloTrackReader::FillInputEvent() - Get event %d from second input AOD file \n", iEntry+fSecondInputFirstEvent);
 419          if(fSecondInputAODTree->GetEntriesFast() <= iEntry+fSecondInputFirstEvent) {
 420                  if(fSecondInputAODTree->GetEntriesFast() == iEntry+fSecondInputFirstEvent)
 421                          printf("AliCaloTrackReader::FillInputEvent() - Skip events from event %d, no more events in second AOD file \n", iEntry);
 422                  return kFALSE;
 423          }
 424
 425          //Get the Event
 426          Int_t nbytes = fSecondInputAODTree->GetEvent(iEntry+fSecondInputFirstEvent);
 427          if ( nbytes == 0 ) {//If nothing in AOD
 428                  printf("AliCaloTrackReader::FillInputEvent() - Nothing in Second AOD input, STOP\n");
 429                  abort() ;
 430          }
 431
 432   }
 433   //if(iEntry > 10) return kFALSE;
 434   if(fFillCTS)   FillInputCTS();
 435   if(fFillEMCAL) FillInputEMCAL();
 436   if(fFillPHOS)  FillInputPHOS();
 437   if(fFillEMCALCells) FillInputEMCALCells();
 438   if(fFillPHOSCells)  FillInputPHOSCells();
 439
 440   return kTRUE ;
 441 }
 442
 443 //__________________________________________________
 444 void AliCaloTrackReader::ResetLists() {
 445   //  Reset lists, called by the analysis maker
 446
 447   if(fAODCTS)   fAODCTS -> Clear();
 448   if(fAODEMCAL) fAODEMCAL -> Clear();
 449   if(fAODPHOS)  fAODPHOS -> Clear();
 450   if(fEMCALCells) fEMCALCells -> Clear();
 451   if(fPHOSCells)  fPHOSCells -> Clear();
 452   if(!fWriteOutputStdAOD){
 453           //Only keep copied tracks and clusters if requested
 454           fOutputEvent->GetTracks()      ->Clear();
 455           fOutputEvent->GetCaloClusters()->Clear();
 456   }
 457 }