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