PWG/CaloTrackCorrBase/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
  16 //_________________________________________________________________________
  17 // Base class for reading data: MonteCarlo, ESD or AOD, of PHOS EMCAL and
  18 // Central Barrel Tracking detectors (CTS).
  19 // Not all MC particles/tracks/clusters are kept, some kinematical/fiducial restrictions are done.
  20 // Mother class of : AliCaloTrackESDReader: Fills ESD data in 3 TObjArrays (PHOS, EMCAL, CTS)
  21 //                 : AliCaloTrackMCReader : Fills Kinematics data in 3 TObjArrays (PHOS, EMCAL, CTS)
  22 //                 : AliCaloTrackAODReader: Fills AOD data in 3 TObjArrays (PHOS, EMCAL, CTS)
  23 //-- Author: Gustavo Conesa (LNF-INFN)
  24 //////////////////////////////////////////////////////////////////////////////
  25
  26
  27 // --- ROOT system ---
  28 #include <TFile.h>
  29 #include <TGeoManager.h>
  30
  31 // ---- ANALYSIS system ----
  32 #include "AliMCEvent.h"
  33 #include "AliAODMCHeader.h"
  34 #include "AliGenPythiaEventHeader.h"
  35 #include "AliESDEvent.h"
  36 #include "AliAODEvent.h"
  37 #include "AliVTrack.h"
  38 #include "AliVParticle.h"
  39 #include "AliMixedEvent.h"
  40 #include "AliESDtrack.h"
  41 #include "AliESDtrackCuts.h"
  42 #include "AliTriggerAnalysis.h"
  43 #include "AliESDVZERO.h"
  44 #include "AliVCaloCells.h"
  45 #include "AliAnalysisManager.h"
  46 #include "AliInputEventHandler.h"
  47
  48 // ---- Detectors ----
  49 #include "AliPHOSGeoUtils.h"
  50 #include "AliEMCALGeometry.h"
  51 #include "AliEMCALRecoUtils.h"
  52
  53 // ---- CaloTrackCorr ---
  54 #include "AliCalorimeterUtils.h"
  55 #include "AliCaloTrackReader.h"
  56
  57 ClassImp(AliCaloTrackReader)
  58
  59
  60 //________________________________________
  61 AliCaloTrackReader::AliCaloTrackReader() :
  62 TObject(),                   fEventNumber(-1), //fCurrentFileName(""),
  63 fDataType(0),                fDebug(0),
  64 fFiducialCut(0x0),           fCheckFidCut(kFALSE),
  65 fComparePtHardAndJetPt(0),   fPtHardAndJetPtFactor(0),
  66 fCTSPtMin(0),                fEMCALPtMin(0),                  fPHOSPtMin(0),
  67 fCTSPtMax(0),                fEMCALPtMax(0),                  fPHOSPtMax(0),
  68 fEMCALTimeCutMin(-10000),    fEMCALTimeCutMax(10000),
  69 fAODBranchList(0x0),
  70 fCTSTracks(0x0),             fEMCALClusters(0x0),             fPHOSClusters(0x0),
  71 fEMCALCells(0x0),            fPHOSCells(0x0),
  72 fInputEvent(0x0),            fOutputEvent(0x0),fMC(0x0),
  73 fFillCTS(0),                 fFillEMCAL(0),                   fFillPHOS(0),
  74 fFillEMCALCells(0),          fFillPHOSCells(0),
  75 fRecalculateClusters(kFALSE),fSelectEmbeddedClusters(kFALSE),
  76 fTrackStatus(0),             fTrackFilterMask(0),
  77 fESDtrackCuts(0),            fConstrainTrack(kFALSE),         fSelectHybridTracks(0),
  78 fTrackMult(0),               fTrackMultEtaCut(0.8),
  79 fReadStack(kFALSE),          fReadAODMCParticles(kFALSE),
  80 fDeltaAODFileName(""),       fFiredTriggerClassName(""),
  81 fEventTriggerMask(0),        fMixEventTriggerMask(0),         fEventTriggerAtSE(0),
  82 fAnaLED(kFALSE),
  83 fTaskName(""),               fCaloUtils(0x0),
  84 fMixedEvent(NULL),           fNMixedEvent(0),                 fVertex(NULL),
  85 fListMixedTracksEvents(),    fListMixedCaloEvents(),
  86 fLastMixedTracksEvent(-1),   fLastMixedCaloEvent(-1),
  87 fWriteOutputDeltaAOD(kFALSE),fOldAOD(kFALSE),                 fCaloFilterPatch(kFALSE),
  88 fEMCALClustersListName(""),  fZvtxCut(0.),
  89 fAcceptFastCluster(kFALSE),  fRemoveLEDEvents(kTRUE),
  90 fDoEventSelection(kFALSE),   fDoV0ANDEventSelection(kFALSE),  fUseEventsWithPrimaryVertex(kFALSE),
  91 fTriggerAnalysis (0x0),
  92 fCentralityClass(""),        fCentralityOpt(0),
  93 fEventPlaneMethod(""),       fImportGeometryFromFile(kFALSE), fImportGeometryFilePath("")
  94 {
  95   //Ctor
  96
  97   //Initialize parameters
  98   InitParameters();
  99 }
 100
 101 //_______________________________________
 102 AliCaloTrackReader::~AliCaloTrackReader()
 103 {
 104   //Dtor
 105
 106   delete fFiducialCut ;
 107
 108   if(fAODBranchList)
 109   {
 110     fAODBranchList->Delete();
 111     delete fAODBranchList ;
 112   }
 113
 114   if(fCTSTracks)
 115   {
 116     if(fDataType!=kMC)fCTSTracks->Clear() ;
 117     else              fCTSTracks->Delete() ;
 118     delete fCTSTracks ;
 119   }
 120
 121   if(fEMCALClusters)
 122   {
 123     if(fDataType!=kMC)fEMCALClusters->Clear("C") ;
 124     else              fEMCALClusters->Delete() ;
 125     delete fEMCALClusters ;
 126   }
 127
 128   if(fPHOSClusters)
 129   {
 130     if(fDataType!=kMC)fPHOSClusters->Clear("C") ;
 131     else              fPHOSClusters->Delete() ;
 132     delete fPHOSClusters ;
 133   }
 134
 135   if(fVertex)
 136   {
 137     for (Int_t i = 0; i < fNMixedEvent; i++)
 138     {
 139       delete [] fVertex[i] ;
 140
 141     }
 142     delete [] fVertex ;
 143         }
 144
 145   delete fESDtrackCuts;
 146   delete fTriggerAnalysis;
 147
 148   //  Pointers not owned, done by the analysis frame
 149   //  if(fInputEvent)  delete fInputEvent ;
 150   //  if(fOutputEvent) delete fOutputEvent ;
 151   //  if(fMC)          delete fMC ;
 152   //  Pointer not owned, deleted by maker
 153   //  if (fCaloUtils) delete fCaloUtils ;
 154
 155 }
 156
 157 //________________________________________________
 158 Bool_t AliCaloTrackReader::ComparePtHardAndJetPt()
 159 {
 160   // Check the event, if the requested ptHard is much larger than the jet pT, then there is a problem.
 161   // Only for PYTHIA.
 162   if(!fReadStack) return kTRUE; //Information not filtered to AOD
 163
 164   if(!strcmp(GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader"))
 165   {
 166     TParticle * jet =  0;
 167     AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetGenEventHeader();
 168     Int_t nTriggerJets =  pygeh->NTriggerJets();
 169     Float_t ptHard = pygeh->GetPtHard();
 170
 171     //if(fDebug > 1) printf("AliMCAnalysisUtils::PythiaEventHeader: Njets: %d, pT Hard %f\n",nTriggerJets, ptHard);
 172     Float_t tmpjet[]={0,0,0,0};
 173     for(Int_t ijet = 0; ijet< nTriggerJets; ijet++)
 174     {
 175       pygeh->TriggerJet(ijet, tmpjet);
 176       jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
 177       //Compare jet pT and pt Hard
 178       //if(fDebug > 1) printf("AliMCAnalysisUtils:: %d pycell jet pT %f\n",ijet, jet->Pt());
 179       if(jet->Pt() > fPtHardAndJetPtFactor * ptHard)
 180       {
 181         printf("AliMCAnalysisUtils::PythiaEventHeader: Njets: %d, pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n",
 182                nTriggerJets, ptHard, jet->Pt(), fPtHardAndJetPtFactor);
 183         return kFALSE;
 184       }
 185     }
 186     if(jet) delete jet;
 187   }
 188
 189   return kTRUE ;
 190
 191 }
 192
 193 //____________________________________________
 194 AliStack* AliCaloTrackReader::GetStack() const
 195 {
 196   //Return pointer to stack
 197   if(fMC)
 198     return fMC->Stack();
 199   else
 200   {
 201     if(fDebug > 1) printf("AliCaloTrackReader::GetStack() - Stack is not available\n");
 202     return 0x0 ;
 203   }
 204 }
 205
 206 //__________________________________________________
 207 TString AliCaloTrackReader::GetFiredTriggerClasses()
 208 {
 209   // List of triggered classes in a TString
 210
 211   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (GetInputEvent());
 212   AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (GetInputEvent());
 213
 214   if     (esdevent) return esdevent->GetFiredTriggerClasses();
 215   else if(aodevent) return aodevent->GetFiredTriggerClasses();
 216   else              return ""; // Mixed Event, MC event, does not have this trigger info
 217
 218 }
 219
 220 //______________________________________________
 221 AliHeader* AliCaloTrackReader::GetHeader() const
 222 {
 223   //Return pointer to header
 224   if(fMC)
 225   {
 226     return fMC->Header();
 227   }
 228   else
 229   {
 230     printf("AliCaloTrackReader::Header is not available\n");
 231     return 0x0 ;
 232   }
 233 }
 234
 235 //______________________________________________________________
 236 AliGenEventHeader* AliCaloTrackReader::GetGenEventHeader() const
 237 {
 238   //Return pointer to Generated event header
 239   if(fMC)
 240   {
 241     return fMC->GenEventHeader();
 242   }
 243   else
 244   {
 245     printf("AliCaloTrackReader::GenEventHeader is not available\n");
 246     return 0x0 ;
 247   }
 248 }
 249
 250 //____________________________________________________________________
 251 TClonesArray* AliCaloTrackReader::GetAODMCParticles(Int_t input) const
 252 {
 253   //Return list of particles in AOD. Do it for the corresponding input event.
 254
 255   TClonesArray * rv = NULL ;
 256   if(fDataType == kAOD)
 257   {
 258     if(input == 0)
 259     {
 260       //Normal input AOD
 261       AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fInputEvent) ;
 262       if(evt)
 263         rv = (TClonesArray*)evt->FindListObject("mcparticles");
 264       else
 265         printf("AliCaloTrackReader::GetAODMCParticles() - wrong AOD input index? %d, or non existing tree? \n",input);
 266     }
 267
 268   }
 269   else
 270   {
 271     printf("AliCaloTrackReader::GetAODMCParticles() - Input are not AODs\n");
 272   }
 273
 274   return rv ;
 275 }
 276
 277 //___________________________________________________________________
 278 AliAODMCHeader* AliCaloTrackReader::GetAODMCHeader(Int_t input) const
 279 {
 280   //Return MC header in AOD. Do it for the corresponding input event.
 281
 282   AliAODMCHeader *mch = NULL;
 283
 284   if(fDataType == kAOD)
 285   {
 286     //Normal input AOD
 287     if(input == 0)
 288     {
 289       mch = (AliAODMCHeader*)((AliAODEvent*)fInputEvent)->FindListObject("mcheader");
 290     }
 291     else
 292     {
 293       printf("AliCaloTrackReader::GetAODMCHeader() - wrong AOD input index, %d\n",input);
 294     }
 295   }
 296   else
 297   {
 298     printf("AliCaloTrackReader::GetAODMCHeader() - Input are not AODs\n");
 299   }
 300
 301   return mch;
 302 }
 303
 304 //_____________________________
 305 void AliCaloTrackReader::Init()
 306 {
 307   //Init reader. Method to be called in AliAnaPartCorrMaker
 308
 309   //printf(" AliCaloTrackReader::Init() %p \n",gGeoManager);
 310
 311   if(fReadStack && fReadAODMCParticles)
 312   {
 313     printf("AliCaloTrackReader::Init() - Cannot access stack and mcparticles at the same time, change them \n");
 314     fReadStack          = kFALSE;
 315     fReadAODMCParticles = kFALSE;
 316   }
 317
 318   // Init geometry, I do not like much to do it like this ...
 319   if(fImportGeometryFromFile && !gGeoManager)
 320   {
 321     if(fImportGeometryFilePath=="") // If not specified, set a default location
 322     fImportGeometryFilePath = "$ALICE_ROOT/PWGGA/EMCALTasks/macros/geometry.root"; // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
 323
 324     printf("AliCaloTrackReader::Init() - Import %s\n",fImportGeometryFilePath.Data());
 325     TGeoManager::Import(fImportGeometryFilePath) ; // default need file "geometry.root" in local dir!!!!
 326   }
 327
 328 }
 329
 330 //_______________________________________
 331 void AliCaloTrackReader::InitParameters()
 332 {
 333   //Initialize the parameters of the analysis.
 334   fDataType   = kESD ;
 335   fCTSPtMin   = 0.1 ;
 336   fEMCALPtMin = 0.1 ;
 337   fPHOSPtMin  = 0.1 ;
 338   fCTSPtMax   = 1000. ;
 339   fEMCALPtMax = 1000. ;
 340   fPHOSPtMax  = 1000. ;
 341
 342   //Do not filter the detectors input by default.
 343   fFillEMCAL      = kFALSE;
 344   fFillPHOS       = kFALSE;
 345   fFillCTS        = kFALSE;
 346   fFillEMCALCells = kFALSE;
 347   fFillPHOSCells  = kFALSE;
 348
 349   fReadStack             = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
 350   fReadAODMCParticles    = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
 351   fDeltaAODFileName      = "deltaAODPartCorr.root";
 352   fFiredTriggerClassName = "";
 353   fEventTriggerMask      = AliVEvent::kAny;
 354   fMixEventTriggerMask   = AliVEvent::kAnyINT;
 355   fEventTriggerAtSE      = kTRUE; // Use only events that pass event selection at SE base class
 356
 357   fAcceptFastCluster = kTRUE;
 358   fAnaLED            = kFALSE;
 359
 360   //We want tracks fitted in the detectors:
 361   //fTrackStatus=AliESDtrack::kTPCrefit;
 362   //fTrackStatus|=AliESDtrack::kITSrefit;
 363   fTrackStatus     = 0;
 364   fTrackFilterMask = 128; //For AODs, but what is the difference between fTrackStatus and fTrackFilterMask?
 365
 366   fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts(); //initialize with TPC only tracks
 367
 368   fConstrainTrack = kFALSE ; // constrain tracks to vertex
 369
 370   fV0ADC[0] = 0;   fV0ADC[1] = 0;
 371   fV0Mul[0] = 0;   fV0Mul[1] = 0;
 372
 373   fZvtxCut   = 10.;
 374
 375   fNMixedEvent = 1;
 376
 377   fPtHardAndJetPtFactor = 7;
 378
 379   //Centrality
 380   fCentralityClass  = "V0M";
 381   fCentralityOpt    = 10;
 382   fCentralityBin[0] = fCentralityBin[1]=-1;
 383
 384   fEventPlaneMethod = "V0";
 385
 386   // Allocate memory (not sure this is the right place)
 387   fCTSTracks       = new TObjArray();
 388   fEMCALClusters   = new TObjArray();
 389   fPHOSClusters    = new TObjArray();
 390   fTriggerAnalysis = new AliTriggerAnalysis;
 391   fAODBranchList   = new TList ;
 392
 393   fImportGeometryFromFile = kFALSE;
 394
 395 }
 396
 397 //________________________________________________________
 398 void AliCaloTrackReader::Print(const Option_t * opt) const
 399 {
 400
 401   //Print some relevant parameters set for the analysis
 402   if(! opt)
 403     return;
 404
 405   printf("***** Print: %s %s ******\n",    GetName(), GetTitle() ) ;
 406   printf("Task name      : %s\n",          fTaskName.Data()) ;
 407   printf("Data type      : %d\n",          fDataType) ;
 408   printf("CTS Min pT     : %2.1f GeV/c\n", fCTSPtMin) ;
 409   printf("EMCAL Min pT   : %2.1f GeV/c\n", fEMCALPtMin) ;
 410   printf("PHOS Min pT    : %2.1f GeV/c\n", fPHOSPtMin) ;
 411   printf("CTS Max pT     : %2.1f GeV/c\n", fCTSPtMax) ;
 412   printf("EMCAL Max pT   : %2.1f GeV/c\n", fEMCALPtMax) ;
 413   printf("PHOS Max pT    : %2.1f GeV/c\n", fPHOSPtMax) ;
 414   printf("EMCAL Time Cut: %3.1f < TOF  < %3.1f\n", fEMCALTimeCutMin, fEMCALTimeCutMax);
 415   printf("Use CTS         =     %d\n",     fFillCTS) ;
 416   printf("Use EMCAL       =     %d\n",     fFillEMCAL) ;
 417   printf("Use PHOS        =     %d\n",     fFillPHOS) ;
 418   printf("Use EMCAL Cells =     %d\n",     fFillEMCALCells) ;
 419   printf("Use PHOS  Cells =     %d\n",     fFillPHOSCells) ;
 420   printf("Track status    =     %d\n", (Int_t) fTrackStatus) ;
 421   printf("AODs Track filter mask  =  %d or hybrid %d\n", (Int_t) fTrackFilterMask,fSelectHybridTracks) ;
 422   printf("Track Mult Eta Cut =  %d\n", (Int_t) fTrackMultEtaCut) ;
 423   printf("Write delta AOD =     %d\n",     fWriteOutputDeltaAOD) ;
 424   printf("Recalculate Clusters = %d\n",    fRecalculateClusters) ;
 425
 426   printf("Use Triggers selected in SE base class %d; If not what trigger Mask? %d; Trigger max for mixed %d \n",
 427          fEventTriggerAtSE, fEventTriggerMask,fMixEventTriggerMask);
 428
 429   if(fComparePtHardAndJetPt)
 430           printf("Compare jet pt and pt hard to accept event, factor = %2.2f",fPtHardAndJetPtFactor);
 431
 432   printf("Read Kine from, stack? %d, AOD ? %d \n", fReadStack, fReadAODMCParticles) ;
 433   printf("Delta AOD File Name =     %s\n", fDeltaAODFileName.Data()) ;
 434   printf("Centrality: Class %s, Option %d, Bin [%d,%d] \n", fCentralityClass.Data(),fCentralityOpt,fCentralityBin[0], fCentralityBin[1]) ;
 435
 436   printf("    \n") ;
 437
 438 }
 439
 440 //_________________________________________________________________________
 441 Bool_t AliCaloTrackReader::FillInputEvent(const Int_t iEntry,
 442                                           const char * /*currentFileName*/)
 443 {
 444   //Fill the event counter and input lists that are needed, called by the analysis maker.
 445
 446   fEventNumber = iEntry;
 447   //fCurrentFileName = TString(currentFileName);
 448   if(!fInputEvent)
 449   {
 450           if(fDebug >= 0) printf("AliCaloTrackReader::FillInputEvent() - Input event not available, skip event analysis\n");
 451           return kFALSE;
 452   }
 453
 454   //Select events only fired by a certain trigger configuration if it is provided
 455   Int_t eventType = 0;
 456   if(fInputEvent->GetHeader())
 457           eventType = ((AliVHeader*)fInputEvent->GetHeader())->GetEventType();
 458
 459   if (GetFiredTriggerClasses().Contains("FAST")  && !GetFiredTriggerClasses().Contains("ALL") && !fAcceptFastCluster)
 460   {
 461     if(fDebug > 0)  printf("AliCaloTrackReader::FillInputEvent - Do not count events from fast cluster, trigger name %s\n",fFiredTriggerClassName.Data());
 462     return kFALSE;
 463   }
 464
 465   //-------------------------------------------------------------------------------------
 466   // Reject event if large clusters with large energy
 467   // Use only for LHC11a data for the moment, and if input is clusterizer V1 or V1+unfolding
 468   // If clusterzer NxN or V2 it does not help
 469   //-------------------------------------------------------------------------------------
 470   Int_t run = fInputEvent->GetRunNumber();
 471   if( fRemoveLEDEvents && run > 146857  && run < 146861 )
 472   {
 473     //printf("Event %d\n",GetEventNumber());
 474
 475     // Count number of cells with energy larger than 0.1 in SM3, cut on this number
 476     Int_t ncellsSM3 = 0;
 477     for(Int_t icell = 0; icell < fInputEvent->GetEMCALCells()->GetNumberOfCells(); icell++)
 478     {
 479       Int_t absID = fInputEvent->GetEMCALCells()->GetCellNumber(icell);
 480       Int_t sm    = GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absID);
 481       if(fInputEvent->GetEMCALCells()->GetAmplitude(icell) > 0.1 && sm==3) ncellsSM3++;
 482     }
 483
 484     Int_t ncellcut = 21;
 485     if(fFiredTriggerClassName.Contains("EMC")) ncellcut = 35;
 486
 487     if(ncellsSM3 >= ncellcut)
 488     {
 489       if(fDebug > 0) printf(" AliCaloTrackReader::FillInputEvent() - reject event with ncells in SM3 %d\n",ncellsSM3);
 490       return kFALSE;
 491     }
 492   }// Remove LED events
 493
 494   // Reject pure LED events?
 495   if( fFiredTriggerClassName  !="" && !fAnaLED)
 496   {
 497     //printf("Event type %d\n",eventType);
 498     if(eventType!=7)
 499       return kFALSE; //Only physics event, do not use for simulated events!!!
 500
 501     if(fDebug > 0)
 502       printf("AliCaloTrackReader::FillInputEvent() - FiredTriggerClass <%s>, selected class <%s>, compare name %d\n",
 503              GetFiredTriggerClasses().Data(),fFiredTriggerClassName.Data(), GetFiredTriggerClasses().Contains(fFiredTriggerClassName));
 504
 505     if( !GetFiredTriggerClasses().Contains(fFiredTriggerClassName) ) return kFALSE;
 506     else if(fDebug > 0) printf("AliCaloTrackReader::FillInputEvent() - Accepted triggered event\n");
 507   }
 508   else if(fAnaLED)
 509   {
 510     //    kStartOfRun =       1,    // START_OF_RUN
 511     //    kEndOfRun =         2,    // END_OF_RUN
 512     //    kStartOfRunFiles =  3,    // START_OF_RUN_FILES
 513     //    kEndOfRunFiles =    4,    // END_OF_RUN_FILES
 514     //    kStartOfBurst =     5,    // START_OF_BURST
 515     //    kEndOfBurst =       6,    // END_OF_BURST
 516     //    kPhysicsEvent =     7,    // PHYSICS_EVENT
 517     //    kCalibrationEvent = 8,    // CALIBRATION_EVENT
 518     //    kFormatError =      9,    // EVENT_FORMAT_ERROR
 519     //    kStartOfData =      10,   // START_OF_DATA
 520     //    kEndOfData =        11,   // END_OF_DATA
 521     //    kSystemSoftwareTriggerEvent   = 12, // SYSTEM_SOFTWARE_TRIGGER_EVENT
 522     //    kDetectorSoftwareTriggerEvent = 13  // DETECTOR_SOFTWARE_TRIGGER_EVENT
 523
 524           if(eventType!=7 && fDebug > 1 )printf("AliCaloTrackReader::FillInputEvent() - DO LED, Event Type <%d>, 8 Calibration \n",  eventType);
 525           if(eventType!=8)return kFALSE;
 526   }
 527
 528   //In case of analysis of events with jets, skip those with jet pt > 5 pt hard
 529   if(fComparePtHardAndJetPt && GetStack())
 530   {
 531     if(!ComparePtHardAndJetPt()) return kFALSE ;
 532   }
 533
 534   //Fill Vertex array
 535   FillVertexArray();
 536   //Reject events with Z vertex too large, only for SE analysis, if not, cut on the analysis code
 537   if(!GetMixedEvent() && TMath::Abs(fVertex[0][2]) > fZvtxCut) return kFALSE;
 538
 539   //------------------------------------------------------
 540   //Event rejection depending on vertex, pileup, v0and
 541   //------------------------------------------------------
 542
 543   if(fUseEventsWithPrimaryVertex)
 544   {
 545     if( !CheckForPrimaryVertex() )              return kFALSE;
 546     if( TMath::Abs(fVertex[0][0] ) < 1.e-6 &&
 547         TMath::Abs(fVertex[0][1] ) < 1.e-6 &&
 548         TMath::Abs(fVertex[0][2] ) < 1.e-6    ) return kFALSE;
 549   }
 550
 551   if(fDoEventSelection)
 552   {
 553     if(!fCaloFilterPatch)
 554     {
 555       //Do not analyze events with pileup
 556       Bool_t bPileup = fInputEvent->IsPileupFromSPD(3, 0.8, 3., 2., 5.); //Default values, if not it does not compile
 557       //Bool_t bPileup = event->IsPileupFromSPD();
 558       if(bPileup) return kFALSE;
 559
 560       if(fDoV0ANDEventSelection)
 561       {
 562         Bool_t bV0AND = kTRUE;
 563         AliESDEvent* esd = dynamic_cast<AliESDEvent*> (fInputEvent);
 564         if(esd)
 565           bV0AND = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0AND);
 566         //else bV0AND = //FIXME FOR AODs
 567         if(!bV0AND) return kFALSE;
 568       }
 569     }//CaloFilter patch
 570     else
 571     {
 572       if(fInputEvent->GetNumberOfCaloClusters() > 0)
 573       {
 574         AliVCluster * calo = fInputEvent->GetCaloCluster(0);
 575         if(calo->GetNLabels() == 4)
 576         {
 577           Int_t * selection = calo->GetLabels();
 578           Bool_t bPileup = selection[0];
 579           if(bPileup) return kFALSE;
 580
 581           Bool_t bGoodV = selection[1];
 582           if(fUseEventsWithPrimaryVertex && !bGoodV) return kFALSE;
 583
 584           if(fDoV0ANDEventSelection)
 585           {
 586             Bool_t bV0AND = selection[2];
 587             if(!bV0AND) return kFALSE;
 588           }
 589
 590           fTrackMult = selection[3];
 591           if(fTrackMult == 0) return kFALSE;
 592         }
 593         else
 594         {
 595           //First filtered AODs, track multiplicity stored there.
 596           fTrackMult = (Int_t) ((AliAODHeader*)fInputEvent->GetHeader())->GetCentrality();
 597           if(fTrackMult == 0) return kFALSE;
 598         }
 599       }//at least one cluster
 600       else
 601       {
 602         //printf("AliCaloTrackReader::FillInputEvent() - No clusters in event\n");
 603         //Remove events with  vertex (0,0,0), bad vertex reconstruction
 604         if(fUseEventsWithPrimaryVertex && TMath::Abs(fVertex[0][0]) < 1.e-6 && TMath::Abs(fVertex[0][1]) < 1.e-6 && TMath::Abs(fVertex[0][2]) < 1.e-6) return kFALSE;
 605
 606         //First filtered AODs, track multiplicity stored there.
 607         fTrackMult = (Int_t) ((AliAODHeader*)fInputEvent->GetHeader())->GetCentrality();
 608         if(fTrackMult == 0) return kFALSE;
 609       }// no cluster
 610     }// CaloFileter patch
 611   }// Event selection
 612   //------------------------------------------------------
 613
 614   //Check if there is a centrality value, PbPb analysis, and if a centrality bin selection is requested
 615   //If we need a centrality bin, we select only those events in the corresponding bin.
 616   if(GetCentrality() && fCentralityBin[0]>=0 && fCentralityBin[1]>=0 && fCentralityOpt==100)
 617   {
 618     Int_t cen = GetEventCentrality();
 619     if(cen > fCentralityBin[1] || cen < fCentralityBin[0]) return kFALSE; //reject events out of bin.
 620   }
 621
 622   //Fill the arrays with cluster/tracks/cells data
 623
 624   if(fFillEMCALCells)
 625     FillInputEMCALCells();
 626
 627   if(fFillPHOSCells)
 628     FillInputPHOSCells();
 629
 630   FillInputVZERO();
 631
 632   if(fEventTriggerAtSE)
 633   {
 634     if(fFillCTS)
 635     {
 636       FillInputCTS();
 637       //Accept events with at least one track
 638       if(fTrackMult == 0 && fDoEventSelection) return kFALSE;
 639     }
 640     if(fFillEMCAL)
 641       FillInputEMCAL();
 642     if(fFillPHOS)
 643       FillInputPHOS();
 644   }
 645   else
 646   {
 647     // In case of mixing analysis, all triggers accepted, but trigger particles selected
 648     // only for the specific trigered events selected here. Mixing done only for MB events,
 649     // tracks array filled also for those events and not the others.
 650
 651     AliAnalysisManager *manager = AliAnalysisManager::GetAnalysisManager();
 652     AliInputEventHandler *inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
 653
 654     if(!inputHandler) return kFALSE ;  // to content coverity
 655
 656     UInt_t isTrigger = inputHandler->IsEventSelected() & fEventTriggerMask;
 657     UInt_t isMB      = inputHandler->IsEventSelected() & fMixEventTriggerMask;
 658
 659     if(fFillCTS && (isMB || isTrigger))
 660     {
 661       FillInputCTS();
 662       //Accept events with at least one track
 663       if(fTrackMult == 0 && fDoEventSelection) return kFALSE;
 664     }
 665
 666     if(isTrigger)
 667     {
 668       if(fFillEMCAL)
 669         FillInputEMCAL();
 670       if(fFillPHOS)
 671         FillInputPHOS();
 672     }
 673   }
 674
 675   return kTRUE ;
 676 }
 677
 678 //___________________________________
 679 void AliCaloTrackReader::ResetLists()
 680 {
 681   //  Reset lists, called by the analysis maker
 682
 683   if(fCTSTracks)       fCTSTracks     -> Clear();
 684   if(fEMCALClusters)   fEMCALClusters -> Clear("C");
 685   if(fPHOSClusters)    fPHOSClusters  -> Clear("C");
 686
 687   fV0ADC[0] = 0;   fV0ADC[1] = 0;
 688   fV0Mul[0] = 0;   fV0Mul[1] = 0;
 689
 690 }
 691
 692 //____________________________________________________________
 693 void AliCaloTrackReader::SetInputEvent(AliVEvent* const input)
 694 {
 695   fInputEvent  = input;
 696   fMixedEvent = dynamic_cast<AliMixedEvent*>(GetInputEvent()) ;
 697   if (fMixedEvent)
 698     fNMixedEvent = fMixedEvent->GetNumberOfEvents() ;
 699
 700   //Delete previous vertex
 701   if(fVertex)
 702   {
 703     for (Int_t i = 0; i < fNMixedEvent; i++)
 704     {
 705       delete [] fVertex[i] ;
 706     }
 707     delete [] fVertex ;
 708   }
 709
 710   fVertex = new Double_t*[fNMixedEvent] ;
 711   for (Int_t i = 0; i < fNMixedEvent; i++)
 712   {
 713     fVertex[i] = new Double_t[3] ;
 714     fVertex[i][0] = 0.0 ;
 715     fVertex[i][1] = 0.0 ;
 716     fVertex[i][2] = 0.0 ;
 717   }
 718 }
 719
 720 //__________________________________________________
 721 Int_t AliCaloTrackReader::GetEventCentrality() const
 722 {
 723   //Return current event centrality
 724
 725   if(GetCentrality())
 726   {
 727     if     (fCentralityOpt==100) return (Int_t) GetCentrality()->GetCentralityPercentile(fCentralityClass); // 100 bins max
 728     else if(fCentralityOpt==10)  return GetCentrality()->GetCentralityClass10(fCentralityClass);// 10 bins max
 729     else if(fCentralityOpt==20)  return GetCentrality()->GetCentralityClass5(fCentralityClass); // 20 bins max
 730     else
 731     {
 732       printf("AliCaloTrackReader::GetEventCentrality() - Unknown centrality option %d, use 10, 20 or 100\n",fCentralityOpt);
 733       return -1;
 734     }
 735   }
 736   else return -1;
 737
 738 }
 739
 740 //_____________________________________________________
 741 Double_t AliCaloTrackReader::GetEventPlaneAngle() const
 742 {
 743   //Return current event centrality
 744
 745   if(GetEventPlane())
 746   {
 747     Float_t ep =  GetEventPlane()->GetEventplane(GetEventPlaneMethod(), GetInputEvent());
 748
 749     if(GetEventPlaneMethod()=="Q" && (ep < 0 || ep > TMath::Pi()))
 750     {
 751       printf("AliCaloTrackReader::GetEventPlaneAngle() -  Bad EP for <Q> method : %f\n",ep);
 752       return -1000;
 753     }
 754     else if(GetEventPlaneMethod().Contains("V0")  )
 755     {
 756       if((ep > TMath::Pi()/2 || ep < -TMath::Pi()/2))
 757       {
 758         printf("AliCaloTrackReader::GetEventPlaneAngle() -  Bad EP for <%s> method : %f\n",GetEventPlaneMethod().Data(), ep);
 759         return -1000;
 760       }
 761
 762       ep+=TMath::Pi()/2; // put same range as for <Q> method
 763
 764     }
 765
 766     //printf("AliCaloTrackReader::GetEventPlaneAngle() = %f\n",ep);
 767     if     (ep > TMath::Pi()) printf("AliCaloTrackReader::GetEventPlaneAngle() - Too large angle = %f\n",ep);
 768     else if(ep < 0          ) printf("AliCaloTrackReader::GetEventPlaneAngle() - Negative angle = %f\n" ,ep);
 769
 770     return ep;
 771   }
 772   else
 773   {
 774     printf("AliCaloTrackReader::GetEventPlaneAngle() -  No EP pointer\n");
 775     return -1000;
 776   }
 777
 778 }
 779
 780 //__________________________________________________________
 781 void AliCaloTrackReader::GetVertex(Double_t vertex[3]) const
 782 {
 783   //Return vertex position to be used for single event analysis
 784   vertex[0]=fVertex[0][0];
 785   vertex[1]=fVertex[0][1];
 786   vertex[2]=fVertex[0][2];
 787 }
 788
 789 //____________________________________________________________
 790 void AliCaloTrackReader::GetVertex(Double_t vertex[3],
 791                                    const Int_t evtIndex) const
 792 {
 793   //Return vertex position for mixed event, recover the vertex in a particular event.
 794
 795   vertex[0]=fVertex[evtIndex][0];  vertex[1]=fVertex[evtIndex][1];  vertex[2]=fVertex[evtIndex][2];
 796
 797 }
 798
 799 //________________________________________
 800 void AliCaloTrackReader::FillVertexArray()
 801 {
 802
 803   //Fill data member with vertex
 804   //In case of Mixed event, multiple vertices
 805
 806   //Delete previous vertex
 807   if(fVertex)
 808   {
 809     for (Int_t i = 0; i < fNMixedEvent; i++)
 810     {
 811       delete [] fVertex[i] ;
 812     }
 813     delete [] fVertex ;
 814   }
 815
 816   fVertex = new Double_t*[fNMixedEvent] ;
 817   for (Int_t i = 0; i < fNMixedEvent; i++)
 818   {
 819     fVertex[i] = new Double_t[3] ;
 820     fVertex[i][0] = 0.0 ;
 821     fVertex[i][1] = 0.0 ;
 822     fVertex[i][2] = 0.0 ;
 823   }
 824
 825   if (!fMixedEvent)
 826   { //Single event analysis
 827     if(fDataType!=kMC)
 828     {
 829
 830       if(fInputEvent->GetPrimaryVertex())
 831       {
 832         fInputEvent->GetPrimaryVertex()->GetXYZ(fVertex[0]);
 833       }
 834       else
 835       {
 836         printf("AliCaloTrackReader::FillVertexArray() - NULL primary vertex\n");
 837         fVertex[0][0]=0.;   fVertex[0][1]=0.;   fVertex[0][2]=0.;
 838       }//Primary vertex pointer do not exist
 839
 840     } else
 841     {//MC read event
 842       fVertex[0][0]=0.;   fVertex[0][1]=0.;   fVertex[0][2]=0.;
 843     }
 844
 845     if(fDebug > 1)
 846       printf("AliCaloTrackReader::FillVertexArray() - Single Event Vertex : %f,%f,%f\n",fVertex[0][0],fVertex[0][1],fVertex[0][2]);
 847
 848   } else
 849   { // MultiEvent analysis
 850     for (Int_t iev = 0; iev < fNMixedEvent; iev++)
 851     {
 852       if (fMixedEvent->GetVertexOfEvent(iev))
 853         fMixedEvent->GetVertexOfEvent(iev)->GetXYZ(fVertex[iev]);
 854       else
 855       { // no vertex found !!!!
 856         AliWarning("No vertex found");
 857       }
 858
 859       if(fDebug > 1)
 860         printf("AliCaloTrackReader::FillVertexArray() - Multi Event %d Vertex : %f,%f,%f\n",iev,fVertex[iev][0],fVertex[iev][1],fVertex[iev][2]);
 861
 862     }
 863   }
 864
 865 }
 866
 867 //_____________________________________
 868 void AliCaloTrackReader::FillInputCTS()
 869 {
 870   //Return array with Central Tracking System (CTS) tracks
 871
 872   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputCTS()\n");
 873
 874   Double_t pTrack[3] = {0,0,0};
 875
 876   Int_t nTracks = fInputEvent->GetNumberOfTracks() ;
 877   fTrackMult    = 0;
 878   Int_t nstatus = 0;
 879
 880   for (Int_t itrack =  0; itrack <  nTracks; itrack++)
 881   {////////////// track loop
 882     AliVTrack * track = (AliVTrack*)fInputEvent->GetTrack(itrack) ; // retrieve track from esd
 883
 884     //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
 885     if (fTrackStatus && !((track->GetStatus() & fTrackStatus) == fTrackStatus))
 886       continue ;
 887
 888     nstatus++;
 889
 890     if     (fDataType==kESD)
 891     {
 892       AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (track);
 893
 894       if(esdTrack && fESDtrackCuts->AcceptTrack(esdTrack))
 895       {
 896         track->GetPxPyPz(pTrack) ;
 897
 898         if(fConstrainTrack)
 899         {
 900           if(esdTrack->GetConstrainedParam())
 901           {
 902             const AliExternalTrackParam* constrainParam = esdTrack->GetConstrainedParam();
 903             esdTrack->Set(constrainParam->GetX(),constrainParam->GetAlpha(),constrainParam->GetParameter(),constrainParam->GetCovariance());
 904             esdTrack->GetConstrainedPxPyPz(pTrack);
 905           }
 906           else continue;
 907
 908         } // use constrained tracks
 909
 910       }
 911       else continue;
 912
 913     } // ESD
 914     else if(fDataType==kAOD)
 915     {
 916       AliAODTrack *aodtrack = dynamic_cast <AliAODTrack*>(track);
 917
 918       if(aodtrack)
 919       {
 920        if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputCTS():AOD track type: %d (primary %d), hybrid? %d \n",
 921                               aodtrack->GetType(),AliAODTrack::kPrimary,
 922                               aodtrack->IsHybridGlobalConstrainedGlobal());
 923
 924
 925         if (fSelectHybridTracks)
 926         {
 927           if (!aodtrack->IsHybridGlobalConstrainedGlobal())       continue ;
 928         }
 929         else
 930         {
 931           if ( aodtrack->TestFilterBit(fTrackFilterMask)==kFALSE) continue ;
 932         }
 933
 934         if (aodtrack->GetType()!= AliAODTrack::kPrimary)          continue ;
 935
 936         if (fDebug > 2 ) printf("AliCaloTrackReader::FillInputCTS(): \t accepted track! \n");
 937
 938         track->GetPxPyPz(pTrack) ;
 939
 940       } // aod track exists
 941       else continue ;
 942
 943     } // AOD
 944
 945     //Count the tracks in eta < 0.9
 946     //printf("Eta %f cut  %f\n",TMath::Abs(track->Eta()),fTrackMultEtaCut);
 947     if(TMath::Abs(track->Eta())< fTrackMultEtaCut) fTrackMult++;
 948
 949     TLorentzVector momentum(pTrack[0],pTrack[1],pTrack[2],0);
 950
 951     if(fCTSPtMin < momentum.Pt() && fCTSPtMax > momentum.Pt())
 952     {
 953       if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"CTS"))
 954         continue;
 955
 956       if(fDebug > 2 && momentum.Pt() > 0.1)
 957         printf("AliCaloTrackReader::FillInputCTS() - Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
 958                momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 959
 960       if (fMixedEvent)
 961       {
 962         track->SetID(itrack);
 963       }
 964
 965       fCTSTracks->Add(track);
 966
 967     }//Pt and Fiducial cut passed.
 968   }// track loop
 969
 970   if(fDebug > 1)
 971     printf("AliCaloTrackReader::FillInputCTS()   - aod entries %d, input tracks %d, pass status %d, multipliticy %d\n", fCTSTracks->GetEntriesFast(), nTracks, nstatus, fTrackMult);//fCTSTracksNormalInputEntries);
 972
 973 }
 974
 975 //__________________________________________________________________
 976 void AliCaloTrackReader::FillInputEMCALAlgorithm(AliVCluster * clus,
 977                                                  const Int_t iclus)
 978 {
 979   //Fill the EMCAL data in the array, do it
 980
 981   Int_t vindex = 0 ;
 982   if (fMixedEvent)
 983     vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
 984
 985   //Reject clusters with bad channels, close to borders and exotic;
 986   if(!GetCaloUtils()->GetEMCALRecoUtils()->IsGoodCluster(clus,GetCaloUtils()->GetEMCALGeometry(),GetEMCALCells(),fInputEvent->GetBunchCrossNumber())) return;
 987
 988   //Mask all cells in collumns facing ALICE thick material if requested
 989   if(GetCaloUtils()->GetNMaskCellColumns())
 990   {
 991     Int_t absId   = -1;
 992     Int_t iSupMod = -1;
 993     Int_t iphi    = -1;
 994     Int_t ieta    = -1;
 995     Bool_t shared = kFALSE;
 996     GetCaloUtils()->GetEMCALRecoUtils()->GetMaxEnergyCell(GetCaloUtils()->GetEMCALGeometry(), GetEMCALCells(),clus,absId,iSupMod,ieta,iphi,shared);
 997     if(GetCaloUtils()->MaskFrameCluster(iSupMod, ieta)) return;
 998   }
 999
1000   if(fSelectEmbeddedClusters)
1001   {
1002     if(clus->GetNLabels()==0 || clus->GetLabel() < 0) return;
1003     //else printf("Embedded cluster,  %d, n label %d label %d  \n",iclus,clus->GetNLabels(),clus->GetLabel());
1004   }
1005
1006   //Float_t pos[3];
1007   //clus->GetPosition(pos);
1008   //printf("Before Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
1009
1010   if(fRecalculateClusters)
1011   {
1012     //Recalibrate the cluster energy
1013     if(GetCaloUtils()->IsRecalibrationOn())
1014     {
1015       Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, GetEMCALCells());
1016
1017       clus->SetE(energy);
1018       //printf("Recalibrated Energy %f\n",clus->E());
1019
1020       GetCaloUtils()->RecalculateClusterShowerShapeParameters(GetEMCALCells(),clus);
1021       GetCaloUtils()->RecalculateClusterPID(clus);
1022
1023     } // recalculate E
1024
1025     //Recalculate distance to bad channels, if new list of bad channels provided
1026     GetCaloUtils()->RecalculateClusterDistanceToBadChannel(GetEMCALCells(),clus);
1027
1028     //Recalculate cluster position
1029     if(GetCaloUtils()->IsRecalculationOfClusterPositionOn())
1030     {
1031       GetCaloUtils()->RecalculateClusterPosition(GetEMCALCells(),clus);
1032       //clus->GetPosition(pos);
1033       //printf("After  Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
1034     }
1035
1036     // Recalculate TOF
1037     if(GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn())
1038     {
1039       Double_t tof      = clus->GetTOF();
1040       Float_t  frac     =-1;
1041       Int_t    absIdMax = GetCaloUtils()->GetMaxEnergyCell(fEMCALCells, clus,frac);
1042
1043       if(fDataType==AliCaloTrackReader::kESD)
1044       {
1045         tof = fEMCALCells->GetCellTime(absIdMax);
1046       }
1047
1048       GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
1049
1050       clus->SetTOF(tof);
1051
1052     }// Time recalibration
1053   }
1054
1055   //Correct non linearity
1056   if(GetCaloUtils()->IsCorrectionOfClusterEnergyOn())
1057   {
1058     GetCaloUtils()->CorrectClusterEnergy(clus) ;
1059     //printf("Linearity Corrected Energy %f\n",clus->E());
1060
1061     //In case of MC analysis, to match resolution/calibration in real data
1062     Float_t rdmEnergy = GetCaloUtils()->GetEMCALRecoUtils()->SmearClusterEnergy(clus);
1063     // printf("\t Energy %f, smeared %f\n", clus->E(),rdmEnergy);
1064     clus->SetE(rdmEnergy);
1065   }
1066
1067   TLorentzVector momentum ;
1068
1069   clus->GetMomentum(momentum, fVertex[vindex]);
1070
1071   if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"EMCAL")) return;
1072
1073   if(fEMCALPtMin > momentum.E() || fEMCALPtMax < momentum.E())         return;
1074
1075   Double_t tof = clus->GetTOF()*1e9;
1076   if(tof < fEMCALTimeCutMin     || tof > fEMCALTimeCutMax)             return;
1077
1078   if(fDebug > 2 && momentum.E() > 0.1)
1079     printf("AliCaloTrackReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
1080            momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
1081
1082   if (fMixedEvent)
1083     clus->SetID(iclus) ;
1084
1085   fEMCALClusters->Add(clus);
1086
1087 }
1088
1089 //_______________________________________
1090 void AliCaloTrackReader::FillInputEMCAL()
1091 {
1092   //Return array with EMCAL clusters in aod format
1093
1094   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputEMCAL()\n");
1095
1096   // First recalibrate cells, time or energy
1097   //  if(GetCaloUtils()->IsRecalibrationOn())
1098   //    GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCells(GetCaloUtils()->GetEMCALGeometry(),
1099   //                                                          GetEMCALCells(),
1100   //                                                          fInputEvent->GetBunchCrossNumber());
1101
1102   //Loop to select clusters in fiducial cut and fill container with aodClusters
1103   if(fEMCALClustersListName=="")
1104   {
1105     Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
1106     for (Int_t iclus =  0; iclus <  nclusters; iclus++)
1107     {
1108       AliVCluster * clus = 0;
1109       if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
1110       {
1111         if (IsEMCALCluster(clus))
1112         {
1113           FillInputEMCALAlgorithm(clus, iclus);
1114         }//EMCAL cluster
1115       }// cluster exists
1116     }// cluster loop
1117
1118     //Recalculate track matching
1119     GetCaloUtils()->RecalculateClusterTrackMatching(fInputEvent);
1120
1121   }//Get the clusters from the input event
1122   else
1123   {
1124     TClonesArray * clusterList = 0x0;
1125
1126     if      (fInputEvent->FindListObject(fEMCALClustersListName))
1127     {
1128       clusterList = dynamic_cast<TClonesArray*> (fInputEvent ->FindListObject(fEMCALClustersListName));
1129     }
1130     else if(fOutputEvent)
1131     {
1132       clusterList = dynamic_cast<TClonesArray*> (fOutputEvent->FindListObject(fEMCALClustersListName));
1133     }
1134
1135     if(!clusterList)
1136     {
1137       //printf("AliCaloTrackReader::FillInputEMCAL() - Wrong name of list with clusters? Try input event <%s>\n",fEMCALClustersListName.Data());
1138       //List not in output event, try input event
1139       clusterList = dynamic_cast<TClonesArray*> (fInputEvent->FindListObject(fEMCALClustersListName));
1140       if(!clusterList)
1141       {
1142         printf("AliCaloTrackReader::FillInputEMCAL() - Wrong name of list with clusters?  <%s>\n",fEMCALClustersListName.Data());
1143         return;
1144       }
1145     }
1146
1147     Int_t nclusters = clusterList->GetEntriesFast();
1148     for (Int_t iclus =  0; iclus <  nclusters; iclus++)
1149     {
1150       AliVCluster * clus = dynamic_cast<AliVCluster*> (clusterList->At(iclus));
1151       //printf("E %f\n",clus->E());
1152       if (clus) FillInputEMCALAlgorithm(clus, iclus);
1153       else printf("AliCaloTrackReader::FillInputEMCAL() - Null cluster in list!\n");
1154     }// cluster loop
1155
1156     // Recalculate track matching, not necessary if already done in the reclusterization task.
1157     // in case it was not done ...
1158     GetCaloUtils()->RecalculateClusterTrackMatching(fInputEvent,clusterList);
1159
1160   }
1161
1162   if(fDebug > 1) printf("AliCaloTrackReader::FillInputEMCAL() - aod entries %d\n",  fEMCALClusters->GetEntriesFast());
1163
1164 }
1165
1166 //______________________________________
1167 void AliCaloTrackReader::FillInputPHOS()
1168 {
1169   //Return array with PHOS clusters in aod format
1170
1171   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputPHOS()\n");
1172
1173   //Loop to select clusters in fiducial cut and fill container with aodClusters
1174   Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
1175   for (Int_t iclus = 0; iclus < nclusters; iclus++)
1176   {
1177     AliVCluster * clus = 0;
1178     if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
1179     {
1180       if (IsPHOSCluster(clus))
1181       {
1182         //Check if the cluster contains any bad channel and if close to calorimeter borders
1183         Int_t vindex = 0 ;
1184         if (fMixedEvent)
1185           vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
1186         if( GetCaloUtils()->ClusterContainsBadChannel("PHOS",clus->GetCellsAbsId(), clus->GetNCells()))
1187           continue;
1188         if(!GetCaloUtils()->CheckCellFiducialRegion(clus, fInputEvent->GetPHOSCells(), fInputEvent, vindex))
1189           continue;
1190
1191         if(fRecalculateClusters)
1192         {
1193           //Recalibrate the cluster energy
1194           if(GetCaloUtils()->IsRecalibrationOn())
1195           {
1196             Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliAODCaloCells*)GetPHOSCells());
1197             clus->SetE(energy);
1198           }
1199         }
1200
1201         TLorentzVector momentum ;
1202
1203         clus->GetMomentum(momentum, fVertex[vindex]);
1204
1205         if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"PHOS")) continue;
1206
1207         if(fPHOSPtMin > momentum.E() || fPHOSPtMax < momentum.E())          continue;
1208
1209         if(fDebug > 2 && momentum.E() > 0.1)
1210           printf("AliCaloTrackReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
1211                  momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
1212
1213
1214         if (fMixedEvent)
1215         {
1216           clus->SetID(iclus) ;
1217         }
1218
1219         fPHOSClusters->Add(clus);
1220
1221       }//PHOS cluster
1222     }//cluster exists
1223   }//esd cluster loop
1224
1225   if(fDebug > 1) printf("AliCaloTrackReader::FillInputPHOS()  - aod entries %d\n",  fPHOSClusters->GetEntriesFast());
1226
1227 }
1228
1229 //____________________________________________
1230 void AliCaloTrackReader::FillInputEMCALCells()
1231 {
1232   //Return array with EMCAL cells in aod format
1233
1234   fEMCALCells = fInputEvent->GetEMCALCells();
1235
1236 }
1237
1238 //___________________________________________
1239 void AliCaloTrackReader::FillInputPHOSCells()
1240 {
1241   //Return array with PHOS cells in aod format
1242
1243   fPHOSCells = fInputEvent->GetPHOSCells();
1244
1245 }
1246
1247 //_______________________________________
1248 void AliCaloTrackReader::FillInputVZERO()
1249 {
1250   //Fill VZERO information in data member, add all the channels information.
1251   AliVVZERO* v0 = fInputEvent->GetVZEROData();
1252   //printf("Init V0: ADC (%d,%d), Multiplicity (%d,%d) \n",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]);
1253
1254   if (v0)
1255   {
1256     AliESDVZERO* esdV0 = dynamic_cast<AliESDVZERO*> (v0);
1257     for (Int_t i = 0; i < 32; i++)
1258     {
1259       if(esdV0)
1260       {//Only available in ESDs
1261         fV0ADC[0] += (Int_t)esdV0->GetAdcV0C(i);
1262         fV0ADC[1] += (Int_t)esdV0->GetAdcV0A(i);
1263       }
1264
1265       fV0Mul[0] += (Int_t)v0->GetMultiplicityV0C(i);
1266       fV0Mul[1] += (Int_t)v0->GetMultiplicityV0A(i);
1267     }
1268     if(fDebug > 0)
1269       printf("V0: ADC (%d,%d), Multiplicity (%d,%d) \n",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]);
1270   }
1271   else
1272   {
1273     if(fDebug > 0)
1274       printf("Cannot retrieve V0 ESD! Run w/ null V0 charges\n ");
1275   }
1276 }
1277
1278
1279 //___________________________________________________________________
1280 Bool_t AliCaloTrackReader::IsEMCALCluster(AliVCluster* cluster) const
1281 {
1282   // Check if it is a cluster from EMCAL. For old AODs cluster type has
1283   // different number and need to patch here
1284
1285   if(fDataType==kAOD && fOldAOD)
1286   {
1287     if (cluster->GetType() == 2) return kTRUE;
1288     else                         return kFALSE;
1289   }
1290   else
1291   {
1292     return cluster->IsEMCAL();
1293   }
1294
1295 }
1296
1297 //___________________________________________________________________
1298 Bool_t AliCaloTrackReader::IsPHOSCluster(AliVCluster * cluster) const
1299 {
1300   //Check if it is a cluster from PHOS.For old AODs cluster type has
1301   // different number and need to patch here
1302
1303   if(fDataType==kAOD && fOldAOD)
1304   {
1305     Int_t type = cluster->GetType();
1306     if (type == 0 || type == 1) return kTRUE;
1307     else                        return kFALSE;
1308   }
1309   else
1310   {
1311     return cluster->IsPHOS();
1312   }
1313
1314 }
1315
1316 //________________________________________________
1317 Bool_t AliCaloTrackReader::CheckForPrimaryVertex()
1318 {
1319   //Check if the vertex was well reconstructed, copy from V0Reader of conversion group
1320   //Only for ESDs ...
1321
1322   AliESDEvent * event = dynamic_cast<AliESDEvent*> (fInputEvent);
1323   if(!event) return kTRUE;
1324
1325   if(event->GetPrimaryVertexTracks()->GetNContributors() > 0)
1326   {
1327     return kTRUE;
1328   }
1329
1330   if(event->GetPrimaryVertexTracks()->GetNContributors() < 1)
1331   {
1332     // SPD vertex
1333     if(event->GetPrimaryVertexSPD()->GetNContributors() > 0)
1334     {
1335       //cout<<"spd vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
1336       return kTRUE;
1337
1338     }
1339     if(event->GetPrimaryVertexSPD()->GetNContributors() < 1)
1340     {
1341       //      cout<<"bad vertex type::"<< event->GetPrimaryVertex()->GetName() << endl;
1342       return kFALSE;
1343     }
1344   }
1345
1346   return kFALSE;
1347
1348 }
1349
1350 //____________________________________________________________
1351 void  AliCaloTrackReader::SetTrackCuts(AliESDtrackCuts * cuts)
1352 {
1353   // Set Track cuts
1354
1355   if(fESDtrackCuts) delete fESDtrackCuts ;
1356
1357   fESDtrackCuts = cuts ;
1358
1359 }
1360