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