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/fiducial 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 // This part is commented: Mixing analysis can be done, input AOD with events
  26 // is opened in the AliCaloTrackReader::Init()
  27
  28 //-- Author: Gustavo Conesa (LNF-INFN)
  29 //////////////////////////////////////////////////////////////////////////////
  30
  31
  32 // --- ROOT system ---
  33 #include "TFile.h"
  34
  35 //---- ANALYSIS system ----
  36 #include "AliCaloTrackReader.h"
  37 #include "AliMCEvent.h"
  38 #include "AliAODMCHeader.h"
  39 #include "AliGenPythiaEventHeader.h"
  40 #include "AliVEvent.h"
  41 #include "AliAODEvent.h"
  42 #include "AliVTrack.h"
  43 #include "AliVParticle.h"
  44 #include "AliMixedEvent.h"
  45
  46 ClassImp(AliCaloTrackReader)
  47
  48
  49 //____________________________________________________________________________
  50   AliCaloTrackReader::AliCaloTrackReader() :
  51     TObject(), fEventNumber(-1), fCurrentFileName(""),fDataType(0), fDebug(0),
  52     fFiducialCut(0x0), fCheckFidCut(kFALSE), fComparePtHardAndJetPt(kFALSE), fPtHardAndJetPtFactor(7),
  53     fCTSPtMin(0), fEMCALPtMin(0),fPHOSPtMin(0),
  54     fAODCTS(new TObjArray()), fAODEMCAL(new TObjArray()), fAODPHOS(new TObjArray()),
  55     fEMCALCells(0x0), fPHOSCells(0x0),
  56     fInputEvent(0x0), fOutputEvent(0x0),fMC(0x0),
  57     fFillCTS(0),fFillEMCAL(0),fFillPHOS(0),
  58     fFillEMCALCells(0),fFillPHOSCells(0),
  59 //    fSecondInputAODTree(0x0), fSecondInputAODEvent(0x0),
  60 //    fSecondInputFileName(""),fSecondInputFirstEvent(0),
  61 //    fAODCTSNormalInputEntries(0), fAODEMCALNormalInputEntries(0),
  62 //    fAODPHOSNormalInputEntries(0),
  63     fTrackStatus(0),
  64     fReadStack(kFALSE), fReadAODMCParticles(kFALSE),
  65     fDeltaAODFileName("deltaAODPartCorr.root"),fFiredTriggerClassName(""),
  66     fAnaLED(kFALSE),fTaskName(""),fCaloUtils(0x0),
  67     fMixedEvent(NULL), fNMixedEvent(1), fVertex(NULL), fWriteOutputStdAOD(kFALSE)
  68 {
  69   //Ctor
  70
  71   //Initialize parameters
  72   InitParameters();
  73 }
  74 /*
  75 //____________________________________________________________________________
  76 AliCaloTrackReader::AliCaloTrackReader(const AliCaloTrackReader & reader) :
  77   TObject(reader), fEventNumber(reader.fEventNumber), fCurrentFileName(reader.fCurrentFileName),
  78   fDataType(reader.fDataType), fDebug(reader.fDebug),
  79   fFiducialCut(reader.fFiducialCut),
  80   fComparePtHardAndJetPt(reader.fComparePtHardAndJetPt),
  81   fPtHardAndJetPtFactor(reader.fPtHardAndJetPtFactor),
  82   fCTSPtMin(reader.fCTSPtMin), fEMCALPtMin(reader.fEMCALPtMin),fPHOSPtMin(reader.fPHOSPtMin),
  83   fAODCTS(new TObjArray(*reader.fAODCTS)),
  84   fAODEMCAL(new TObjArray(*reader.fAODEMCAL)),
  85   fAODPHOS(new TObjArray(*reader.fAODPHOS)),
  86   fEMCALCells(new TNamed(*reader.fEMCALCells)),
  87   fPHOSCells(new TNamed(*reader.fPHOSCells)),
  88   fInputEvent(reader.fInputEvent), fOutputEvent(reader.fOutputEvent), fMC(reader.fMC),
  89   fFillCTS(reader.fFillCTS),fFillEMCAL(reader.fFillEMCAL),fFillPHOS(reader.fFillPHOS),
  90   fFillEMCALCells(reader.fFillEMCALCells),fFillPHOSCells(reader.fFillPHOSCells),
  91   fSecondInputAODTree(reader.fSecondInputAODTree),
  92   fSecondInputAODEvent(reader.fSecondInputAODEvent),
  93   fSecondInputFileName(reader.fSecondInputFileName),
  94   fSecondInputFirstEvent(reader.fSecondInputFirstEvent),
  95   fAODCTSNormalInputEntries(reader.fAODCTSNormalInputEntries),
  96   fAODEMCALNormalInputEntries(reader.fAODEMCALNormalInputEntries),
  97   fAODPHOSNormalInputEntries(reader.fAODPHOSNormalInputEntries),
  98   fTrackStatus(reader.fTrackStatus),
  99   fReadStack(reader.fReadStack), fReadAODMCParticles(reader.fReadAODMCParticles),
 100   fFiredTriggerClassName(reader.fFiredTriggerClassName),
 101   fAnaLED(reader.fAnaLED),
 102   fTaskName(reader.fTaskName),
 103   fCaloUtils(new AliCalorimeterUtils(*reader.fCaloUtils))
 104 {
 105   // cpy ctor
 106 }
 107 */
 108 //_________________________________________________________________________
 109 //AliCaloTrackReader & AliCaloTrackReader::operator = (const AliCaloTrackReader & source)
 110 //{
 111 //  // assignment operator
 112 //
 113 //  if(&source == this) return *this;
 114 //
 115 //  fDataType    = source.fDataType ;
 116 //  fDebug       = source.fDebug ;
 117 //  fEventNumber = source.fEventNumber ;
 118 //  fCurrentFileName = source.fCurrentFileName ;
 119 //  fFiducialCut = source.fFiducialCut;
 120 //
 121 //  fComparePtHardAndJetPt = source.fComparePtHardAndJetPt;
 122 //  fPtHardAndJetPtFactor  = source.fPtHardAndJetPtFactor;
 123 //
 124 //  fCTSPtMin    = source.fCTSPtMin ;
 125 //  fEMCALPtMin  = source.fEMCALPtMin ;
 126 //  fPHOSPtMin   = source.fPHOSPtMin ;
 127 //
 128 //  fAODCTS     = new TObjArray(*source.fAODCTS) ;
 129 //  fAODEMCAL   = new TObjArray(*source.fAODEMCAL) ;
 130 //  fAODPHOS    = new TObjArray(*source.fAODPHOS) ;
 131 //  fEMCALCells = new TNamed(*source.fEMCALCells) ;
 132 //  fPHOSCells  = new TNamed(*source.fPHOSCells) ;
 133 //
 134 //  fInputEvent  = source.fInputEvent;
 135 //  fOutputEvent = source.fOutputEvent;
 136 //  fMC          = source.fMC;
 137 //
 138 //  fFillCTS        = source.fFillCTS;
 139 //  fFillEMCAL      = source.fFillEMCAL;
 140 //  fFillPHOS       = source.fFillPHOS;
 141 //  fFillEMCALCells = source.fFillEMCALCells;
 142 //  fFillPHOSCells  = source.fFillPHOSCells;
 143 //
 144 //  fSecondInputAODTree    = source.fSecondInputAODTree;
 145 //  fSecondInputAODEvent   = source.fSecondInputAODEvent;
 146 //  fSecondInputFileName   = source.fSecondInputFileName;
 147 //  fSecondInputFirstEvent = source.fSecondInputFirstEvent;
 148 //
 149 //  fAODCTSNormalInputEntries   = source.fAODCTSNormalInputEntries;
 150 //  fAODEMCALNormalInputEntries = source.fAODEMCALNormalInputEntries;
 151 //  fAODPHOSNormalInputEntries  = source.fAODPHOSNormalInputEntries;
 152 //
 153 //  fTrackStatus        = source.fTrackStatus;
 154 //  fReadStack          = source.fReadStack;
 155 //  fReadAODMCParticles = source.fReadAODMCParticles;
 156 //
 157 //  fDeltaAODFileName   = source.fDeltaAODFileName;
 158 //
 159 //  fFiredTriggerClassName = source.fFiredTriggerClassName  ;
 160 //
 161 //  return *this;
 162 //
 163 //}
 164
 165 //_________________________________
 166 AliCaloTrackReader::~AliCaloTrackReader() {
 167   //Dtor
 168
 169   if(fFiducialCut) delete fFiducialCut ;
 170
 171   if(fAODCTS){
 172     fAODCTS->Clear() ;
 173     delete fAODCTS ;
 174   }
 175
 176   if(fAODEMCAL){
 177     fAODEMCAL->Clear() ;
 178     delete fAODEMCAL ;
 179   }
 180
 181   if(fAODPHOS){
 182     fAODPHOS->Clear() ;
 183     delete fAODPHOS ;
 184   }
 185
 186   if(fEMCALCells){
 187     fEMCALCells->Clear() ;
 188     delete fEMCALCells ;
 189   }
 190
 191   if(fPHOSCells){
 192     fPHOSCells->Clear() ;
 193     delete fPHOSCells ;
 194   }
 195   if (fMixedEvent) {
 196     for (Int_t i = 0; i < fNMixedEvent; i++) {
 197       delete [] fVertex[i] ;
 198     }
 199   }
 200   delete [] fVertex ;
 201
 202 //  Pointers not owned, done by the analysis frame
 203 //  if(fInputEvent)  delete fInputEvent ;
 204 //  if(fOutputEvent) delete fOutputEvent ;
 205 //  if(fMC)          delete fMC ;
 206
 207 //  if(fSecondInputAODTree){
 208 //    fSecondInputAODTree->Clear();
 209 //    delete fSecondInputAODTree;
 210 //  }
 211 //
 212 //  if(fSecondInputAODEvent) delete fSecondInputAODEvent ;
 213
 214   //  Pointer not owned, deleted by maker
 215   //if (fCaloUtils) delete fCaloUtils ;
 216
 217 }
 218
 219 //_________________________________________________________________________
 220 Bool_t AliCaloTrackReader::ComparePtHardAndJetPt(){
 221         // Check the event, if the requested ptHard is much larger than the jet pT, then there is a problem.
 222         // Only for PYTHIA.
 223         if(!fReadStack) return kTRUE; //Information not filtered to AOD
 224
 225         if(!strcmp(GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader")){
 226                 TParticle * jet =  new TParticle;
 227                 AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetGenEventHeader();
 228                 Int_t nTriggerJets =  pygeh->NTriggerJets();
 229                 Float_t ptHard = pygeh->GetPtHard();
 230
 231                 //if(fDebug > 1) printf("AliMCAnalysisUtils::PythiaEventHeader: Njets: %d, pT Hard %f\n",nTriggerJets, ptHard);
 232             Float_t tmpjet[]={0,0,0,0};
 233                 for(Int_t ijet = 0; ijet< nTriggerJets; ijet++){
 234                         pygeh->TriggerJet(ijet, tmpjet);
 235                         jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
 236                         //Compare jet pT and pt Hard
 237                         //if(fDebug > 1) printf("AliMCAnalysisUtils:: %d pycell jet pT %f\n",ijet, jet->Pt());
 238                         if(jet->Pt() > fPtHardAndJetPtFactor * ptHard) {
 239                                 printf("AliMCAnalysisUtils::PythiaEventHeader: Njets: %d, pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n",
 240                                            nTriggerJets, ptHard, jet->Pt(), fPtHardAndJetPtFactor);
 241                                 return kFALSE;
 242                         }
 243                 }
 244         }
 245
 246         return kTRUE ;
 247
 248 }
 249
 250 //____________________________________________________________________________
 251 AliStack* AliCaloTrackReader::GetStack() const {
 252   //Return pointer to stack
 253   if(fMC)
 254     return fMC->Stack();
 255   else{
 256     if(fDebug > 1) printf("AliCaloTrackReader::GetStack() - Stack is not available\n");
 257     return 0x0 ;
 258   }
 259 }
 260
 261 //____________________________________________________________________________
 262 AliHeader* AliCaloTrackReader::GetHeader() const {
 263   //Return pointer to header
 264   if(fMC)
 265     return fMC->Header();
 266   else{
 267     printf("AliCaloTrackReader::Header is not available\n");
 268     return 0x0 ;
 269   }
 270 }
 271 //____________________________________________________________________________
 272 AliGenEventHeader* AliCaloTrackReader::GetGenEventHeader() const {
 273   //Return pointer to Generated event header
 274   if(fMC)
 275     return fMC->GenEventHeader();
 276   else{
 277     printf("AliCaloTrackReader::GenEventHeader is not available\n");
 278     return 0x0 ;
 279   }
 280 }
 281
 282 //____________________________________________________________________________
 283 TClonesArray* AliCaloTrackReader::GetAODMCParticles(Int_t input) const {
 284   //Return list of particles in AOD. Do it for the corresponding input event.
 285
 286   TClonesArray * rv = NULL ;
 287   if(fDataType == kAOD){
 288
 289     if(input == 0){
 290       //Normal input AOD
 291       AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fInputEvent) ;
 292       if(evt)
 293         rv = (TClonesArray*)evt->FindListObject("mcparticles");
 294       else
 295         printf("AliCaloTrackReader::GetAODMCParticles() - wrong AOD input index? %d, or non existing tree? \n",input);
 296
 297     } //else if(input == 1 && fSecondInputAODEvent){ //Second input AOD
 298 //      rv = (TClonesArray*) fSecondInputAODEvent->FindListObject("mcparticles");
 299 //    }
 300
 301   } else {
 302       printf("AliCaloTrackReader::GetAODMCParticles() - Input are not AODs\n");
 303   }
 304
 305   return rv ;
 306 }
 307
 308 //____________________________________________________________________________
 309 AliAODMCHeader* AliCaloTrackReader::GetAODMCHeader(Int_t input) const {
 310         //Return MC header in AOD. Do it for the corresponding input event.
 311   AliAODMCHeader *mch = NULL;
 312         if(fDataType == kAOD){
 313                 //Normal input AOD
 314                 if(input == 0) {
 315       mch = (AliAODMCHeader*)((AliAODEvent*)fInputEvent)->FindListObject("mcheader");
 316     }
 317 //              //Second input AOD
 318 //              else if(input == 1){
 319 //       mch = (AliAODMCHeader*) fSecondInputAODEvent->FindListObject("mcheader");
 320 //  }
 321                 else {
 322                         printf("AliCaloTrackReader::GetAODMCHeader() - wrong AOD input index, %d\n",input);
 323                 }
 324         }
 325         else {
 326                 printf("AliCaloTrackReader::GetAODMCHeader() - Input are not AODs\n");
 327         }
 328
 329   return mch;
 330 }
 331
 332 //_______________________________________________________________
 333 void AliCaloTrackReader::Init()
 334 {
 335         //Init reader. Method to be called in AliAnaPartCorrMaker
 336
 337         //Get the file with second input events if the filename is given
 338         //Get the tree and connect the AODEvent. Only with AODs
 339
 340         if(fReadStack && fReadAODMCParticles){
 341                 printf("AliCaloTrackReader::Init() - Cannot access stack and mcparticles at the same time, change them \n");
 342                 fReadStack = kFALSE;
 343                 fReadAODMCParticles = kFALSE;
 344         }
 345
 346 //      if(fSecondInputFileName!=""){
 347 //              if(fDataType == kAOD){
 348 //                      TFile * input2   = new TFile(fSecondInputFileName,"read");
 349 //                      printf("AliCaloTrackReader::Init() - Second input file opened: %s, size %d \n", input2->GetName(), (Int_t) input2->GetSize());
 350 //                      fSecondInputAODTree = (TTree*) input2->Get("aodTree");
 351 //                      if(fSecondInputAODTree) printf("AliCaloTrackReader::Init() - Second input tree opened: %s, entries %d \n",
 352 //                                                                                 fSecondInputAODTree->GetName(), (Int_t) fSecondInputAODTree->GetEntries());
 353 //                      else{
 354 //                       printf("AliCaloTrackReader::Init() - Second input tree not available, STOP \n");
 355 //                       abort();
 356 //                      }
 357 //                      fSecondInputAODEvent = new AliAODEvent;
 358 //                      fSecondInputAODEvent->ReadFromTree(fSecondInputAODTree);
 359 //                      if(fSecondInputFirstEvent >= fSecondInputAODTree->GetEntriesFast()){
 360 //                              printf("AliCaloTrackReader::Init() - Requested first event of second input %d, is larger than number of events %d, STOP\n",
 361 //                                         fSecondInputFirstEvent, (Int_t) fSecondInputAODTree->GetEntriesFast());
 362 //                              abort();
 363 //                      }
 364 //              }
 365 //              else printf("AliCaloTrackReader::Init() - Second input not added, reader is not AOD\n");
 366 //      }
 367 }
 368
 369 //_______________________________________________________________
 370 void AliCaloTrackReader::InitParameters()
 371 {
 372   //Initialize the parameters of the analysis.
 373   fDataType   = kESD ;
 374   fCTSPtMin   = 0.2 ;
 375   fEMCALPtMin = 0.2 ;
 376   fPHOSPtMin  = 0.2 ;
 377
 378   //Do not filter the detectors input by default.
 379   fFillEMCAL      = kFALSE;
 380   fFillPHOS       = kFALSE;
 381   fFillCTS        = kFALSE;
 382   fFillEMCALCells = kFALSE;
 383   fFillPHOSCells  = kFALSE;
 384
 385   //fSecondInputFileName   = "" ;
 386   //fSecondInputFirstEvent = 0 ;
 387   fReadStack             = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
 388   fReadAODMCParticles    = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
 389   fDeltaAODFileName      = "deltaAODPartCorr.root";
 390   fFiredTriggerClassName      = "";
 391
 392   fAnaLED = kFALSE;
 393 }
 394
 395 //________________________________________________________________
 396 void AliCaloTrackReader::Print(const Option_t * opt) const
 397 {
 398
 399   //Print some relevant parameters set for the analysis
 400   if(! opt)
 401     return;
 402
 403   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
 404   printf("Task name      : %s\n", fTaskName.Data()) ;
 405   printf("Data type      : %d\n", fDataType) ;
 406   printf("CTS Min pT     : %2.1f GeV/c\n", fCTSPtMin) ;
 407   printf("EMCAL Min pT   : %2.1f GeV/c\n", fEMCALPtMin) ;
 408   printf("PHOS Min pT    : %2.1f GeV/c\n", fPHOSPtMin) ;
 409   printf("Use CTS         =     %d\n", fFillCTS) ;
 410   printf("Use EMCAL       =     %d\n", fFillEMCAL) ;
 411   printf("Use PHOS        =     %d\n", fFillPHOS) ;
 412   printf("Use EMCAL Cells =     %d\n", fFillEMCALCells) ;
 413   printf("Use PHOS  Cells =     %d\n", fFillPHOSCells) ;
 414   printf("Track status    =     %d\n", (Int_t) fTrackStatus) ;
 415   printf("Write std AOD   =     %d\n", fWriteOutputStdAOD) ;
 416
 417   if(fComparePtHardAndJetPt)
 418           printf("Compare jet pt and pt hard to accept event, factor = %2.2f",fPtHardAndJetPtFactor);
 419
 420 //  if(fSecondInputFileName!="") {
 421 //        printf("Second Input File Name     =     %s\n", fSecondInputFileName.Data()) ;
 422 //        printf("Second Input First Event   =     %d\n", fSecondInputFirstEvent) ;
 423 //  }
 424
 425   printf("Read Kine from, stack? %d, AOD ? %d \n", fReadStack, fReadAODMCParticles) ;
 426   printf("Delta AOD File Name =     %s\n", fDeltaAODFileName.Data()) ;
 427   printf("    \n") ;
 428 }
 429
 430 //___________________________________________________
 431 Bool_t AliCaloTrackReader::FillInputEvent(const Int_t iEntry, const char * currentFileName) {
 432   //Fill the event counter and input lists that are needed, called by the analysis maker.
 433
 434   fEventNumber = iEntry;
 435   fCurrentFileName = TString(currentFileName);
 436   if(!fInputEvent) {
 437           if(fDebug >= 0) printf("AliCaloTrackReader::FillInputEvent() - Input event not available, skip event analysis\n");
 438           return kFALSE;
 439   }
 440   //Select events only fired by a certain trigger configuration if it is provided
 441   Int_t eventType = 0;
 442   if(fInputEvent->GetHeader())
 443           eventType = ((AliVHeader*)fInputEvent->GetHeader())->GetEventType();
 444   if( fFiredTriggerClassName  !="" && !fAnaLED){
 445     if(eventType!=7)
 446       return kFALSE; //Only physics event, do not use for simulated events!!!
 447     if(fDebug > 0)
 448       printf("AliCaloTrackReader::FillInputEvent() - FiredTriggerClass <%s>, selected class <%s>, compare name %d\n",
 449              GetFiredTriggerClasses().Data(),fFiredTriggerClassName.Data(), GetFiredTriggerClasses().Contains(fFiredTriggerClassName));
 450     if( !GetFiredTriggerClasses().Contains(fFiredTriggerClassName) ) return kFALSE;
 451   }
 452   else if(fAnaLED){
 453 //        kStartOfRun =       1,    // START_OF_RUN
 454 //        kEndOfRun =         2,    // END_OF_RUN
 455 //        kStartOfRunFiles =  3,    // START_OF_RUN_FILES
 456 //        kEndOfRunFiles =    4,    // END_OF_RUN_FILES
 457 //        kStartOfBurst =     5,    // START_OF_BURST
 458 //        kEndOfBurst =       6,    // END_OF_BURST
 459 //        kPhysicsEvent =     7,    // PHYSICS_EVENT
 460 //        kCalibrationEvent = 8,    // CALIBRATION_EVENT
 461 //        kFormatError =      9,    // EVENT_FORMAT_ERROR
 462 //        kStartOfData =      10,   // START_OF_DATA
 463 //        kEndOfData =        11,   // END_OF_DATA
 464 //        kSystemSoftwareTriggerEvent   = 12, // SYSTEM_SOFTWARE_TRIGGER_EVENT
 465 //        kDetectorSoftwareTriggerEvent = 13  // DETECTOR_SOFTWARE_TRIGGER_EVENT
 466
 467           if(eventType!=7 && fDebug > 1 )printf("AliCaloTrackReader::FillInputEvent() - DO LED, Event Type <%d>, 8 Calibration \n",  eventType);
 468           if(eventType!=8)return kFALSE;
 469   }
 470
 471   if(fOutputEvent && (fDataType != kAOD) && ((fOutputEvent->GetCaloClusters())->GetEntriesFast()!=0 ||(fOutputEvent->GetTracks())->GetEntriesFast()!=0)){
 472     if (fFillCTS || fFillEMCAL || fFillPHOS) {
 473       printf("AliCaloTrackReader::AODCaloClusters or AODTracks already filled by the filter, do not use the ESD reader, use the AOD reader, STOP\n");
 474                   abort();
 475     }
 476   }
 477
 478   //In case of analysis of events with jets, skip those with jet pt > 5 pt hard
 479   if(fComparePtHardAndJetPt && GetStack()) {
 480     if(!ComparePtHardAndJetPt()) return kFALSE ;
 481   }
 482
 483   //In case of mixing events with other AOD file
 484  // if(fDataType == kAOD && fSecondInputAODTree){
 485 //
 486 //    if(fDebug > 1)
 487 //      printf("AliCaloTrackReader::FillInputEvent() - Get event %d from second input AOD file \n", iEntry+fSecondInputFirstEvent);
 488 //    if(fSecondInputAODTree->GetEntriesFast() <= iEntry+fSecondInputFirstEvent) {
 489 //      if(fSecondInputAODTree->GetEntriesFast() == iEntry+fSecondInputFirstEvent)
 490 //                       printf("AliCaloTrackReader::FillInputEvent() - Skip events from event %d, no more events in second AOD file \n", iEntry);
 491 //      return kFALSE;
 492 //    }
 493 //
 494 //    //Get the Event
 495 //    Int_t nbytes = fSecondInputAODTree->GetEvent(iEntry+fSecondInputFirstEvent);
 496 //    if ( nbytes == 0 ) {//If nothing in AOD
 497 //      printf("AliCaloTrackReader::FillInputEvent() - Nothing in Second AOD input, STOP\n");
 498 //      abort() ;
 499 //    }
 500 //
 501 //  }
 502
 503
 504   for (Int_t iev = 0; iev < fNMixedEvent; iev++) {
 505     if (!fMixedEvent) {
 506       GetVertex() ;
 507     }
 508     else
 509       fMixedEvent->GetVertexOfEvent(iev)->GetXYZ(fVertex[iev]);
 510   }
 511
 512   if(fFillEMCALCells)
 513     FillInputEMCALCells();
 514   if(fFillPHOSCells)
 515     FillInputPHOSCells();
 516
 517   if(fFillCTS)
 518     FillInputCTS();
 519   if(fFillEMCAL)
 520     FillInputEMCAL();
 521   if(fFillPHOS)
 522     FillInputPHOS();
 523
 524
 525   return kTRUE ;
 526 }
 527
 528 //__________________________________________________
 529 void AliCaloTrackReader::ResetLists() {
 530   //  Reset lists, called by the analysis maker
 531
 532   if(fAODCTS)   fAODCTS -> Clear();
 533   if(fAODEMCAL) fAODEMCAL -> Clear();
 534   if(fAODPHOS)  fAODPHOS -> Clear();
 535   if(fEMCALCells) fEMCALCells -> Clear();
 536   if(fPHOSCells)  fPHOSCells -> Clear();
 537
 538 }
 539
 540 //____________________________________________________________________________
 541 void AliCaloTrackReader::SetInputEvent(AliVEvent* const input)
 542 {
 543   fInputEvent  = input;
 544   fMixedEvent = dynamic_cast<AliMixedEvent*>(GetInputEvent()) ;
 545   if (fMixedEvent) {
 546     fNMixedEvent = fMixedEvent->GetNumberOfEvents() ;
 547   }
 548   fVertex = new Double_t*[fNMixedEvent] ;
 549   for (Int_t i = 0; i < fNMixedEvent; i++) {
 550     fVertex[i] = new Double_t[3] ;
 551     fVertex[i][0] = 0.0 ;
 552     fVertex[i][1] = 0.0 ;
 553     fVertex[i][2] = 0.0 ;
 554   }
 555 }
 556
 557 //____________________________________________________________________________
 558 Double_t * AliCaloTrackReader::GetVertex() {
 559     //Return vertex position
 560   if (fMixedEvent)
 561     return NULL ;
 562   if(fInputEvent->GetPrimaryVertex())
 563     fInputEvent->GetPrimaryVertex()->GetXYZ(fVertex[0]);
 564   else {
 565     printf("AliCaloTrackReader::GetVertex() - No vertex available, InputEvent()->GetPrimaryVertex() = 0, return (0,0,0)\n");
 566     fVertex[0][0]=0.0; fVertex[0][1]=0.0; fVertex[0][2]=0.0;
 567   }
 568   return fVertex[0] ;
 569 }
 570
 571   //____________________________________________________________________________
 572 void AliCaloTrackReader::FillInputCTS() {
 573     //Return array with Central Tracking System (CTS) tracks
 574   if(fDebug > 2 ) printf("AliCaloTrackAODReader::FillInputCTS()\n");
 575   Int_t nTracks   = fInputEvent->GetNumberOfTracks() ;
 576   Int_t naod = 0;
 577   Double_t p[3];
 578   for (Int_t itrack =  0; itrack <  nTracks; itrack++) {////////////// track loop
 579     AliVTrack * track = (AliVTrack*)fInputEvent->GetTrack(itrack) ; // retrieve track from esd
 580
 581       //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
 582     if (fTrackStatus && !((track->GetStatus() & fTrackStatus) == fTrackStatus))
 583       continue ;
 584
 585     track->GetPxPyPz(p) ;
 586     TLorentzVector momentum(p[0],p[1],p[2],0);
 587
 588     if(fCTSPtMin < momentum.Pt()){
 589
 590       if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"CTS"))
 591         continue;
 592
 593       if(fDebug > 2 && momentum.Pt() > 0.1)
 594         printf("AliCaloTrackAODReader::FillInputCTS() - Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
 595                momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 596
 597       if (fMixedEvent) {
 598         track->SetID(itrack);
 599       }
 600       if(fWriteOutputStdAOD){
 601         AliAODTrack * aodTrack = dynamic_cast<AliAODTrack*> (track) ;
 602         if (aodTrack) {
 603           AliAODTrack* newtrack = new((*(fOutputEvent->GetTracks()))[naod++]) AliAODTrack(*aodTrack);
 604           fAODCTS->Add(newtrack); //Use AOD stored in output for references.
 605         } else {
 606           AliError("Can only write AOD tracks to output AOD !") ;
 607         }
 608       }
 609       else {
 610         fAODCTS->Add(track);
 611       }
 612     }//Pt and Fiducial cut passed.
 613   }// track loop
 614
 615   //fAODCTSNormalInputEntries = fAODCTS->GetEntriesFast();
 616   if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputCTS()   - aod entries %d\n", fAODCTS->GetEntriesFast());//fAODCTSNormalInputEntries);
 617
 618     //  //If second input event available, add the clusters.
 619     //  if(fSecondInputAODTree && fSecondInputAODEvent){
 620     //    nTracks   = fSecondInputAODEvent->GetNumberOfTracks() ;
 621     //    if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputCTS()   - Add second input tracks, entries %d\n", nTracks) ;
 622     //    for (Int_t itrack =  0; itrack <  nTracks; itrack++) {////////////// track loop
 623     //            AliAODTrack * track = ((AliAODEvent*)fSecondInputAODEvent)->GetTrack(itrack) ; // retrieve track from esd
 624     //
 625     //            //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
 626     //            if (fTrackStatus && !((track->GetStatus() & fTrackStatus) == fTrackStatus)) continue ;
 627     //
 628     //            track->GetPxPyPz(p) ;
 629     //            TLorentzVector momentum(p[0],p[1],p[2],0);
 630     //
 631     //            if(fCTSPtMin < momentum.Pt()){
 632     //
 633     //                    if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"CTS")) continue;
 634     //
 635     //                    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",
 636     //                                                                 momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 637     //
 638     //                    if(fWriteOutputStdAOD){
 639     //                            AliAODTrack* newtrack = new((*(fOutputEvent->GetTracks()))[naod++]) AliAODTrack(*track);
 640     //                            fAODCTS->Add(newtrack); //Use AOD stored in output for references.
 641     //                    }
 642     //                    else fAODCTS->Add(track);
 643     //
 644     //            }//Pt and Fiducial cut passed.
 645     //    }// track loop
 646     //
 647     //    if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputCTS()   - aod normal entries %d, after second input %d\n", fAODCTSNormalInputEntries, fAODCTS->GetEntriesFast());
 648     //  }       //second input loop
 649     //
 650 }
 651
 652   //____________________________________________________________________________
 653 void AliCaloTrackReader::FillInputEMCAL() {
 654     //Return array with EMCAL clusters in aod format
 655   if(fDebug > 2 ) printf("AliCaloTrackAODReader::FillInputEMCAL()\n");
 656
 657   Int_t naod =  (fOutputEvent->GetCaloClusters())->GetEntriesFast();
 658     //Loop to select clusters in fiducial cut and fill container with aodClusters
 659   Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
 660   for (Int_t iclus =  0; iclus <  nclusters; iclus++) {
 661     AliVCluster * clus = 0;
 662     if ( (clus = fInputEvent->GetCaloCluster(iclus)) ) {
 663       if (clus->IsEMCAL()){
 664
 665         //Check if the cluster contains any bad channel and if close to calorimeter borders
 666         Int_t vindex = 0 ;
 667         if (fMixedEvent)
 668           vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
 669
 670         if(GetCaloUtils()->ClusterContainsBadChannel("EMCAL",clus->GetCellsAbsId(), clus->GetNCells()))
 671           continue;
 672         if(!GetCaloUtils()->CheckCellFiducialRegion(clus, (AliVCaloCells*)fInputEvent->GetEMCALCells(), fInputEvent, vindex))
 673           continue;
 674
 675         TLorentzVector momentum ;
 676
 677         clus->GetMomentum(momentum, fVertex[vindex]);
 678
 679         if(fEMCALPtMin < momentum.Pt()){
 680
 681           if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"EMCAL"))
 682             continue;
 683
 684           if(fDebug > 2 && momentum.E() > 0.1)
 685             printf("AliCaloTrackAODReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
 686                    momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 687
 688           //Recalibrate the cluster energy
 689           if(GetCaloUtils()->IsRecalibrationOn()) {
 690             Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliAODCaloCells*)GetEMCALCells());
 691             clus->SetE(energy);
 692           }
 693
 694           if (fMixedEvent) {
 695             clus->SetID(iclus) ;
 696           }
 697           if(fWriteOutputStdAOD){
 698             AliAODCaloCluster * aodClus = dynamic_cast<AliAODCaloCluster*> (clus) ;
 699             if (aodClus) {
 700               AliAODCaloCluster * newclus = new((*(fOutputEvent->GetCaloClusters()))[naod++])AliAODCaloCluster(*aodClus);
 701               fAODEMCAL->Add(newclus);
 702             } else {
 703               AliError("Can only write AOD clusters to output AOD !") ;
 704             }
 705           } else {
 706             fAODEMCAL->Add(clus);
 707           }
 708         }//Pt and Fiducial cut passed.
 709       }//EMCAL cluster
 710     }// cluster exists
 711   }// cluster loop
 712   //fAODEMCALNormalInputEntries = fAODEMCAL->GetEntriesFast();
 713   if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputEMCAL() - aod entries %d\n",  fAODEMCAL->GetEntriesFast());//fAODEMCALNormalInputEntries);
 714
 715     //If second input event available, add the clusters.
 716     //  if(fSecondInputAODTree && fSecondInputAODEvent){
 717     //    GetSecondInputAODVertex(v);
 718     //    nclusters = ((AliAODEvent*)fSecondInputAODEvent)->GetNumberOfCaloClusters();
 719     //    if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputEMCAL() - Add second input clusters, entries %d\n", nclusters) ;
 720     //          for (Int_t iclus =  0; iclus < nclusters; iclus++) {
 721     //                  AliAODCaloCluster * clus = 0;
 722     //                  if ( (clus = ((AliAODEvent*)fSecondInputAODEvent)->GetCaloCluster(iclus)) ) {
 723     //                          if (clus->IsEMCAL()){
 724     //                                  TLorentzVector momentum ;
 725     //                                  clus->GetMomentum(momentum, v);
 726     //
 727     //                                  if(fEMCALPtMin < momentum.Pt()){
 728     //
 729     //                                          if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"EMCAL")) continue;
 730     //
 731     //                                          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",
 732     //                                                                                                                                  momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 733     //                                          if(fWriteOutputStdAOD){
 734     //                                                  AliAODCaloCluster * newclus = new((*(fOutputEvent->GetCaloClusters()))[naod++])AliAODCaloCluster(*clus);
 735     //                                              fAODEMCAL->Add(newclus);
 736     //                                          }
 737     //                                          else fAODEMCAL->Add(clus);
 738     //                                  }//Pt and Fiducial cut passed.
 739     //                          }//EMCAL cluster
 740     //                  }// cluster exists
 741     //          }// cluster loop
 742     //
 743     //    if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputEMCAL() - aod normal entries %d, after second input %d\n", fAODEMCALNormalInputEntries, fAODEMCAL->GetEntriesFast());
 744     //
 745     //  } //second input loop
 746 }
 747
 748   //____________________________________________________________________________
 749 void AliCaloTrackReader::FillInputPHOS() {
 750     //Return array with PHOS clusters in aod format
 751   if(fDebug > 2 ) printf("AliCaloTrackAODReader::FillInputPHOS()\n");
 752
 753     //Get vertex for momentum calculation
 754
 755   Int_t naod =  (fOutputEvent->GetCaloClusters())->GetEntriesFast();
 756     //Loop to select clusters in fiducial cut and fill container with aodClusters
 757   Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
 758   for (Int_t iclus = 0; iclus < nclusters; iclus++) {
 759     AliVCluster * clus = 0;
 760     if ( (clus = fInputEvent->GetCaloCluster(iclus)) ) {
 761       if (clus->IsPHOS()){
 762           //Check if the cluster contains any bad channel and if close to calorimeter borders
 763         Int_t vindex = 0 ;
 764         if (fMixedEvent)
 765           vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
 766         if( GetCaloUtils()->ClusterContainsBadChannel("PHOS",clus->GetCellsAbsId(), clus->GetNCells()))
 767           continue;
 768         if(!GetCaloUtils()->CheckCellFiducialRegion(clus, fInputEvent->GetPHOSCells(), fInputEvent, vindex))
 769           continue;
 770
 771         TLorentzVector momentum ;
 772
 773         clus->GetMomentum(momentum, fVertex[vindex]);
 774
 775         if(fPHOSPtMin < momentum.Pt()){
 776
 777           if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"PHOS"))
 778             continue;
 779
 780           if(fDebug > 2 && momentum.E() > 0.1)
 781             printf("AliCaloTrackAODReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
 782                    momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 783
 784             //Recalibrate the cluster energy
 785           if(GetCaloUtils()->IsRecalibrationOn()) {
 786             Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliAODCaloCells*)GetPHOSCells());
 787             clus->SetE(energy);
 788           }
 789
 790           if (fMixedEvent) {
 791             clus->SetID(iclus) ;
 792           }
 793
 794           if(fWriteOutputStdAOD){
 795             AliAODCaloCluster * aodClus = dynamic_cast<AliAODCaloCluster*> (clus) ;
 796             if (aodClus) {
 797               AliAODCaloCluster * newclus = new((*(fOutputEvent->GetCaloClusters()))[naod++])AliAODCaloCluster(*aodClus);
 798               fAODPHOS->Add(newclus);
 799             } else {
 800               AliError("Can only write AOD clusters to output AOD !") ;
 801             }
 802           } else {
 803             fAODPHOS->Add(clus);
 804           }
 805         }//Pt and Fiducial cut passed.
 806       }//PHOS cluster
 807     }//cluster exists
 808   }//esd cluster loop
 809
 810   //fAODPHOSNormalInputEntries = fAODPHOS->GetEntriesFast() ;
 811   if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputPHOS()  - aod entries %d\n",  fAODPHOS->GetEntriesFast());//fAODPHOSNormalInputEntries);
 812
 813     //If second input event available, add the clusters.
 814     //  if(fSecondInputAODTree && fSecondInputAODEvent){
 815     //    GetSecondInputAODVertex(v);
 816     //    nclusters = ((AliAODEvent*)fSecondInputAODEvent)->GetNumberOfCaloClusters();
 817     //    if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputPHOS()  - Add second input clusters, entries %d\n", nclusters);
 818     //          for (Int_t iclus =  0; iclus < nclusters; iclus++) {
 819     //                  AliAODCaloCluster * clus = 0;
 820     //                  if ( (clus = ((AliAODEvent*)fSecondInputAODEvent)->GetCaloCluster(iclus)) ) {
 821     //                          if (clus->IsPHOS()){
 822     //                                  TLorentzVector momentum ;
 823     //                                  clus->GetMomentum(momentum, v);
 824     //
 825     //                                  if(fPHOSPtMin < momentum.Pt()){
 826     //
 827     //                                          if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"PHOS")) continue;
 828     //
 829     //                                          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",
 830     //                                                                                                                                  momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 831     //                                          if(fWriteOutputStdAOD){
 832     //                                                  AliAODCaloCluster * newclus = new((*(fOutputEvent->GetCaloClusters()))[naod++])AliAODCaloCluster(*clus);
 833     //                                              fAODPHOS->Add(newclus);
 834     //                                          }
 835     //                                          else fAODPHOS->Add(clus);
 836     //                                  }//Pt and Fiducial cut passed.
 837     //                          }//PHOS cluster
 838     //                  }// cluster exists
 839     //          }// cluster loop
 840     //          if(fDebug > 1) printf("AliCaloTrackAODReader::FillInputPHOS()  - aod normal entries %d, after second input %d\n", fAODPHOSNormalInputEntries, fAODPHOS->GetEntriesFast());
 841     //  }       //second input loop
 842
 843 }
 844
 845   //____________________________________________________________________________
 846 void AliCaloTrackReader::FillInputEMCALCells() {
 847     //Return array with EMCAL cells in aod format
 848
 849   fEMCALCells = fInputEvent->GetEMCALCells();
 850
 851 }
 852
 853   //____________________________________________________________________________
 854 void AliCaloTrackReader::FillInputPHOSCells() {
 855     //Return array with PHOS cells in aod format
 856
 857   fPHOSCells = fInputEvent->GetPHOSCells();
 858
 859 }
 860