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 fComparePtHardAndClusterPt(0),fPtHardAndClusterPtFactor(0),
  67 fCTSPtMin(0),                fEMCALPtMin(0),                  fPHOSPtMin(0),
  68 fCTSPtMax(0),                fEMCALPtMax(0),                  fPHOSPtMax(0),
  69 fUseEMCALTimeCut(1),         fUseParamTimeCut(0),             fUseTrackTimeCut(0),
  70 fEMCALTimeCutMin(-10000),    fEMCALTimeCutMax(10000),
  71 fEMCALParamTimeCutMin(),     fEMCALParamTimeCutMax(),
  72 fTrackTimeCutMin(-10000),    fTrackTimeCutMax(10000),
  73 fUseTrackDCACut(0),
  74 fAODBranchList(0x0),
  75 fCTSTracks(0x0),             fEMCALClusters(0x0),             fPHOSClusters(0x0),
  76 fEMCALCells(0x0),            fPHOSCells(0x0),
  77 fInputEvent(0x0),            fOutputEvent(0x0),fMC(0x0),
  78 fFillCTS(0),                 fFillEMCAL(0),                   fFillPHOS(0),
  79 fFillEMCALCells(0),          fFillPHOSCells(0),
  80 fRecalculateClusters(kFALSE),fSelectEmbeddedClusters(kFALSE),
  81 fTrackStatus(0),             fTrackFilterMask(0),
  82 fESDtrackCuts(0),            fESDtrackComplementaryCuts(0),   fConstrainTrack(kFALSE),
  83 fSelectHybridTracks(0),      fSelectSPDHitTracks(kFALSE),
  84 fTrackMult(0),               fTrackMultEtaCut(0.9),
  85 fReadStack(kFALSE),          fReadAODMCParticles(kFALSE),
  86 fDeltaAODFileName(""),       fFiredTriggerClassName(""),
  87 fEventTriggerMask(0),        fMixEventTriggerMask(0),         fEventTriggerAtSE(0),
  88 fAnaLED(kFALSE),
  89 fTaskName(""),               fCaloUtils(0x0),
  90 fMixedEvent(NULL),           fNMixedEvent(0),                 fVertex(NULL),
  91 fListMixedTracksEvents(),    fListMixedCaloEvents(),
  92 fLastMixedTracksEvent(-1),   fLastMixedCaloEvent(-1),
  93 fWriteOutputDeltaAOD(kFALSE),fOldAOD(kFALSE),                 fCaloFilterPatch(kFALSE),
  94 fEMCALClustersListName(""),  fZvtxCut(0.),
  95 fAcceptFastCluster(kFALSE),  fRemoveLEDEvents(kTRUE),
  96 fRemoveExoticEvents(kFALSE), fExoticTrigger(3),               fIsExoticEvent(kFALSE),
  97 fDoEventSelection(kFALSE),   fDoV0ANDEventSelection(kFALSE),
  98 fDoVertexBCEventSelection(kFALSE),
  99 fDoRejectNoTrackEvents(kFALSE),
 100 fUseEventsWithPrimaryVertex(kFALSE),
 101 fTriggerAnalysis (0x0),      fTimeStampEventSelect(0),
 102 fTimeStampEventFracMin(0),   fTimeStampEventFracMax(0),
 103 fTimeStampRunMin(0),         fTimeStampRunMax(0),
 104 fNPileUpClusters(-1),        fNNonPileUpClusters(-1),         fNPileUpClustersCut(3),
 105 fVertexBC(-200),             fRecalculateVertexBC(0),
 106 fCentralityClass(""),            fCentralityOpt(0),
 107 fEventPlaneMethod(""),       fImportGeometryFromFile(kFALSE), fImportGeometryFilePath("")
 108 {
 109   //Ctor
 110
 111   //Initialize parameters
 112   InitParameters();
 113 }
 114
 115 //_______________________________________
 116 AliCaloTrackReader::~AliCaloTrackReader()
 117 {
 118   //Dtor
 119
 120   delete fFiducialCut ;
 121
 122   if(fAODBranchList)
 123   {
 124     fAODBranchList->Delete();
 125     delete fAODBranchList ;
 126   }
 127
 128   if(fCTSTracks)
 129   {
 130     if(fDataType!=kMC)fCTSTracks->Clear() ;
 131     else              fCTSTracks->Delete() ;
 132     delete fCTSTracks ;
 133   }
 134
 135   if(fEMCALClusters)
 136   {
 137     if(fDataType!=kMC)fEMCALClusters->Clear("C") ;
 138     else              fEMCALClusters->Delete() ;
 139     delete fEMCALClusters ;
 140   }
 141
 142   if(fPHOSClusters)
 143   {
 144     if(fDataType!=kMC)fPHOSClusters->Clear("C") ;
 145     else              fPHOSClusters->Delete() ;
 146     delete fPHOSClusters ;
 147   }
 148
 149   if(fVertex)
 150   {
 151     for (Int_t i = 0; i < fNMixedEvent; i++)
 152     {
 153       delete [] fVertex[i] ;
 154
 155     }
 156     delete [] fVertex ;
 157         }
 158
 159   delete fESDtrackCuts;
 160   delete fESDtrackComplementaryCuts;
 161   delete fTriggerAnalysis;
 162
 163   //  Pointers not owned, done by the analysis frame
 164   //  if(fInputEvent)  delete fInputEvent ;
 165   //  if(fOutputEvent) delete fOutputEvent ;
 166   //  if(fMC)          delete fMC ;
 167   //  Pointer not owned, deleted by maker
 168   //  if (fCaloUtils) delete fCaloUtils ;
 169
 170 }
 171
 172 //________________________________________________________________________
 173 Bool_t  AliCaloTrackReader::AcceptDCA(const Float_t pt, const Float_t dca)
 174 {
 175   // Accept track if DCA is smaller than function
 176
 177   Float_t cut = fTrackDCACut[0]+fTrackDCACut[1]/TMath::Power(pt,fTrackDCACut[2]);
 178
 179   if(TMath::Abs(dca) < cut)
 180     return kTRUE;
 181   else
 182     return kFALSE;
 183
 184 }
 185
 186 //________________________________________________
 187 Bool_t AliCaloTrackReader::ComparePtHardAndJetPt()
 188 {
 189   // Check the event, if the requested ptHard is much smaller than the jet pT, then there is a problem.
 190   // Only for PYTHIA.
 191
 192   //printf("AliCaloTrackReader::ComparePtHardAndJetPt() - GenHeaderName : %s\n",GetGenEventHeader()->ClassName());
 193
 194   if(!strcmp(GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader"))
 195   {
 196     TParticle * jet =  0;
 197     AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetGenEventHeader();
 198     Int_t nTriggerJets =  pygeh->NTriggerJets();
 199     Float_t ptHard = pygeh->GetPtHard();
 200
 201     if(fDebug > 1)
 202       printf("AliCaloTrackReader::ComparePtHardAndJetPt() - Njets: %d, pT Hard %f\n",nTriggerJets, ptHard);
 203
 204     Float_t tmpjet[]={0,0,0,0};
 205     for(Int_t ijet = 0; ijet< nTriggerJets; ijet++)
 206     {
 207       pygeh->TriggerJet(ijet, tmpjet);
 208       jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
 209
 210       if(fDebug > 1)
 211         printf("AliCaloTrackReader::ComparePtHardAndJetPt() - jet %d; pycell jet pT %f\n",ijet, jet->Pt());
 212
 213       //Compare jet pT and pt Hard
 214       if(jet->Pt() > fPtHardAndJetPtFactor * ptHard)
 215       {
 216         printf("AliCaloTrackReader::ComparePtHardAndJetPt() - Reject jet event with : pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n",
 217                 ptHard, jet->Pt(), fPtHardAndJetPtFactor);
 218         return kFALSE;
 219       }
 220     }
 221
 222     if(jet) delete jet;
 223   }
 224
 225   return kTRUE ;
 226
 227 }
 228
 229 //____________________________________________________________________
 230 Bool_t AliCaloTrackReader::ComparePtHardAndClusterPt()
 231 {
 232   // Check the event, if the requested ptHard is smaller than the calorimeter cluster E, then there is a problem.
 233   // Only for PYTHIA.
 234
 235   if(!strcmp(GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader"))
 236   {
 237     AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetGenEventHeader();
 238     Float_t ptHard = pygeh->GetPtHard();
 239
 240     Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
 241     for (Int_t iclus =  0; iclus <  nclusters; iclus++)
 242     {
 243       AliVCluster * clus = fInputEvent->GetCaloCluster(iclus) ;
 244       Float_t ecluster = clus->E();
 245
 246       if(ecluster > fPtHardAndClusterPtFactor*ptHard)
 247       {
 248         printf("AliCaloTrackReader::ComparePtHardAndClusterPt() - Reject : ecluster %2.2f, calo %d, factor %2.2f, ptHard %f\n",ecluster,clus->GetType(),fPtHardAndClusterPtFactor,ptHard);
 249
 250         return kFALSE;
 251       }
 252     }
 253
 254   }
 255
 256   return kTRUE ;
 257
 258 }
 259
 260 //____________________________________________
 261 AliStack* AliCaloTrackReader::GetStack() const
 262 {
 263   //Return pointer to stack
 264   if(fMC)
 265     return fMC->Stack();
 266   else
 267   {
 268     if(fDebug > 1) printf("AliCaloTrackReader::GetStack() - Stack is not available\n");
 269     return 0x0 ;
 270   }
 271 }
 272
 273 //__________________________________________________
 274 TString AliCaloTrackReader::GetFiredTriggerClasses()
 275 {
 276   // List of triggered classes in a TString
 277
 278   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (GetInputEvent());
 279   AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (GetInputEvent());
 280
 281   if     (esdevent) return esdevent->GetFiredTriggerClasses();
 282   else if(aodevent) return aodevent->GetFiredTriggerClasses();
 283   else              return ""; // Mixed Event, MC event, does not have this trigger info
 284
 285 }
 286
 287 //______________________________________________
 288 AliHeader* AliCaloTrackReader::GetHeader() const
 289 {
 290   //Return pointer to header
 291   if(fMC)
 292   {
 293     return fMC->Header();
 294   }
 295   else
 296   {
 297     printf("AliCaloTrackReader::Header is not available\n");
 298     return 0x0 ;
 299   }
 300 }
 301
 302 //______________________________________________________________
 303 AliGenEventHeader* AliCaloTrackReader::GetGenEventHeader() const
 304 {
 305   //Return pointer to Generated event header
 306   if     (ReadStack() && fMC)
 307   {
 308     return fMC->GenEventHeader();
 309   }
 310   else if(ReadAODMCParticles() && GetAODMCHeader())
 311   {
 312     //printf("AliCaloTrackReader::GetGenEventHeader() - N headers %d\n",GetAODMCHeader()->GetNCocktailHeaders());
 313     if( GetAODMCHeader()->GetNCocktailHeaders() > 0)
 314       return GetAODMCHeader()->GetCocktailHeader(0) ;
 315     else
 316       return 0x0;
 317   }
 318   else
 319   {
 320     //printf("AliCaloTrackReader::GetGenEventHeader() - MC header not available! \n");
 321     return 0;
 322   }
 323 }
 324
 325 //____________________________________________________________________
 326 TClonesArray* AliCaloTrackReader::GetAODMCParticles() const
 327 {
 328   //Return list of particles in AOD. Do it for the corresponding input event.
 329
 330   TClonesArray * rv = NULL ;
 331   if(fDataType == kAOD)
 332   {
 333     //Normal input AOD
 334     AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fInputEvent) ;
 335     if(evt)
 336       rv = (TClonesArray*)evt->FindListObject("mcparticles");
 337     else
 338       printf("AliCaloTrackReader::GetAODMCParticles() - Null AOD event \n");
 339   }
 340   else
 341   {
 342     printf("AliCaloTrackReader::GetAODMCParticles() - Input are not AODs\n");
 343   }
 344
 345   return rv ;
 346 }
 347
 348 //________________________________________________________
 349 AliAODMCHeader* AliCaloTrackReader::GetAODMCHeader() const
 350 {
 351   //Return MC header in AOD. Do it for the corresponding input event.
 352
 353   AliAODMCHeader *mch = NULL;
 354
 355   if(fDataType == kAOD)
 356   {
 357     AliAODEvent * aod = dynamic_cast<AliAODEvent*> (fInputEvent);
 358     if(aod) mch = dynamic_cast<AliAODMCHeader*>(aod->FindListObject("mcHeader"));
 359   }
 360   else
 361   {
 362     printf("AliCaloTrackReader::GetAODMCHeader() - Input are not AODs\n");
 363   }
 364
 365   return mch;
 366 }
 367
 368 //___________________________________________________________
 369 Int_t AliCaloTrackReader::GetVertexBC(const AliVVertex * vtx)
 370 {
 371   // Get the vertex BC
 372
 373   Int_t vertexBC=vtx->GetBC();
 374   if(!fRecalculateVertexBC) return vertexBC;
 375
 376   // In old AODs BC not stored, recalculate it
 377   // loop over the global track and select those which have small DCA to primary vertex (e.g. primary).
 378   // If at least one of these primaries has valid BC != 0, then this vertex is a pile-up candidate.
 379   // Execute after CTS
 380   Double_t bz  = fInputEvent->GetMagneticField();
 381   Bool_t   bc0 = kFALSE;
 382   Int_t    ntr = GetCTSTracks()->GetEntriesFast();
 383   //printf("N Tracks %d\n",ntr);
 384
 385   for(Int_t i = 0 ; i < ntr ; i++)
 386   {
 387     AliVTrack * track =  (AliVTrack*) (GetCTSTracks()->At(i));
 388
 389     //Check if has TOF info, if not skip
 390     ULong_t status  = track->GetStatus();
 391     Bool_t  okTOF   = (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ;
 392     vertexBC        = track->GetTOFBunchCrossing(bz);
 393     Float_t pt      = track->Pt();
 394
 395     if(!okTOF) continue;
 396
 397     // Get DCA x, y
 398     Double_t dca[2]   = {1e6,1e6};
 399     Double_t covar[3] = {1e6,1e6,1e6};
 400     track->PropagateToDCA(vtx,bz,100.,dca,covar);
 401
 402     if(AcceptDCA(pt,dca[0]))
 403     {
 404       if     (vertexBC !=0 && fVertexBC != AliVTrack::kTOFBCNA) return vertexBC;
 405       else if(vertexBC == 0)                                    bc0 = kTRUE;
 406     }
 407   }
 408
 409   if( bc0 ) vertexBC = 0 ;
 410   else      vertexBC = AliVTrack::kTOFBCNA ;
 411
 412   return vertexBC;
 413
 414 }
 415
 416 //_____________________________
 417 void AliCaloTrackReader::Init()
 418 {
 419   //Init reader. Method to be called in AliAnaPartCorrMaker
 420
 421   //printf(" AliCaloTrackReader::Init() %p \n",gGeoManager);
 422
 423   if(fReadStack && fReadAODMCParticles)
 424   {
 425     printf("AliCaloTrackReader::Init() - Cannot access stack and mcparticles at the same time, change them \n");
 426     fReadStack          = kFALSE;
 427     fReadAODMCParticles = kFALSE;
 428   }
 429
 430   // Init geometry, I do not like much to do it like this ...
 431   if(fImportGeometryFromFile && !gGeoManager)
 432   {
 433     if(fImportGeometryFilePath=="") // If not specified, set a default location
 434     fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2011.root"; // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
 435
 436     printf("AliCaloTrackReader::Init() - Import %s\n",fImportGeometryFilePath.Data());
 437     TGeoManager::Import(fImportGeometryFilePath) ; // default need file "geometry.root" in local dir!!!!
 438   }
 439
 440   if(!fESDtrackCuts)
 441     fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts(); //initialize with TPC only tracks
 442
 443 }
 444
 445 //_______________________________________
 446 void AliCaloTrackReader::InitParameters()
 447 {
 448   //Initialize the parameters of the analysis.
 449   fDataType   = kESD ;
 450   fCTSPtMin   = 0.1 ;
 451   fEMCALPtMin = 0.1 ;
 452   fPHOSPtMin  = 0.1 ;
 453   fCTSPtMax   = 1000. ;
 454   fEMCALPtMax = 1000. ;
 455   fPHOSPtMax  = 1000. ;
 456
 457   //Track DCA cuts
 458   // dca_xy cut = 0.0105+0.0350/TMath::Power(pt,1.1);
 459   fTrackDCACut[0] = 0.0105;
 460   fTrackDCACut[1] = 0.0350;
 461   fTrackDCACut[2] = 1.1;
 462
 463   //Do not filter the detectors input by default.
 464   fFillEMCAL      = kFALSE;
 465   fFillPHOS       = kFALSE;
 466   fFillCTS        = kFALSE;
 467   fFillEMCALCells = kFALSE;
 468   fFillPHOSCells  = kFALSE;
 469
 470   fReadStack             = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
 471   fReadAODMCParticles    = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
 472   fDeltaAODFileName      = "deltaAODPartCorr.root";
 473   fFiredTriggerClassName = "";
 474   fEventTriggerMask      = AliVEvent::kAny;
 475   fMixEventTriggerMask   = AliVEvent::kAnyINT;
 476   fEventTriggerAtSE      = kTRUE; // Use only events that pass event selection at SE base class
 477
 478   fAcceptFastCluster = kTRUE;
 479   fAnaLED            = kFALSE;
 480
 481   //We want tracks fitted in the detectors:
 482   //fTrackStatus=AliESDtrack::kTPCrefit;
 483   //fTrackStatus|=AliESDtrack::kITSrefit;
 484   fTrackStatus     = 0;
 485   fTrackFilterMask = 128; //For AODs, but what is the difference between fTrackStatus and fTrackFilterMask?
 486
 487   fESDtrackCuts = 0;
 488   fESDtrackComplementaryCuts = 0;
 489
 490   fConstrainTrack = kFALSE ; // constrain tracks to vertex
 491
 492   fV0ADC[0] = 0;   fV0ADC[1] = 0;
 493   fV0Mul[0] = 0;   fV0Mul[1] = 0;
 494
 495   fZvtxCut   = 10.;
 496
 497   fNMixedEvent = 1;
 498
 499   fPtHardAndJetPtFactor     = 7.;
 500   fPtHardAndClusterPtFactor = 1.;
 501
 502   //Centrality
 503   fCentralityClass  = "V0M";
 504   fCentralityOpt    = 10;
 505   fCentralityBin[0] = fCentralityBin[1]=-1;
 506
 507   fEventPlaneMethod = "V0";
 508
 509   // Allocate memory (not sure this is the right place)
 510   fCTSTracks       = new TObjArray();
 511   fEMCALClusters   = new TObjArray();
 512   fPHOSClusters    = new TObjArray();
 513   fTriggerAnalysis = new AliTriggerAnalysis;
 514   fAODBranchList   = new TList ;
 515
 516   fImportGeometryFromFile = kFALSE;
 517
 518   fPileUpParamSPD[0] = 3   ; fPileUpParamSPD[1] = 0.8 ;
 519   fPileUpParamSPD[2] = 3.0 ; fPileUpParamSPD[3] = 2.0 ; fPileUpParamSPD[4] = 5.0;
 520
 521   // Parametrized time cut (LHC11d)
 522   fEMCALParamTimeCutMin[0] =-5; fEMCALParamTimeCutMin[1] =-1 ; fEMCALParamTimeCutMin[2] = 3.5 ; fEMCALParamTimeCutMin[3] = 1.  ;
 523   fEMCALParamTimeCutMax[0] = 5; fEMCALParamTimeCutMax[1] = 50; fEMCALParamTimeCutMax[2] = 0.45; fEMCALParamTimeCutMax[3] = 1.25;
 524
 525   // Parametrized time cut (LHC11c)
 526   //fEMCALParamTimeCutMin[0] =-5;   fEMCALParamTimeCutMin[1] =-1 ; fEMCALParamTimeCutMin[2] = 1.87; fEMCALParamTimeCutMin[3] = 0.4;
 527   //fEMCALParamTimeCutMax[0] = 3.5; fEMCALParamTimeCutMax[1] = 50; fEMCALParamTimeCutMax[2] = 0.15; fEMCALParamTimeCutMax[3] = 1.6;
 528
 529   fTimeStampRunMin = -1;
 530   fTimeStampRunMax = 1e12;
 531   fTimeStampEventFracMin = -1;
 532   fTimeStampEventFracMax = 2;
 533
 534   for(Int_t i = 0; i < 19; i++)
 535   {
 536     fEMCalBCEvent   [i] = 0;
 537     fEMCalBCEventCut[i] = 0;
 538     fTrackBCEvent   [i] = 0;
 539     fTrackBCEventCut[i] = 0;
 540   }
 541
 542 }
 543
 544 //___________________________________________________________
 545 Bool_t AliCaloTrackReader::IsInTimeWindow(const Double_t tof, const Float_t energy) const
 546 {
 547   // Cluster time selection window
 548
 549   // Parametrized cut depending on E
 550   if(fUseParamTimeCut)
 551   {
 552     Float_t minCut= fEMCALParamTimeCutMin[0]+fEMCALParamTimeCutMin[1]*TMath::Exp(-(energy-fEMCALParamTimeCutMin[2])/fEMCALParamTimeCutMin[3]);
 553     Float_t maxCut= fEMCALParamTimeCutMax[0]+fEMCALParamTimeCutMax[1]*TMath::Exp(-(energy-fEMCALParamTimeCutMax[2])/fEMCALParamTimeCutMax[3]);
 554     //printf("tof %f, minCut %f, maxCut %f\n",tof,minCut,maxCut);
 555     if( tof < minCut || tof > maxCut )  return kFALSE ;
 556   }
 557
 558   //In any case, the time should to be larger than the fixed window ...
 559   if( tof < fEMCALTimeCutMin  || tof > fEMCALTimeCutMax )  return kFALSE ;
 560
 561   return kTRUE ;
 562 }
 563
 564 //________________________________________________
 565 Bool_t AliCaloTrackReader::IsPileUpFromSPD() const
 566 {
 567   // Check if event is from pile-up determined by SPD
 568   // Default values: (3, 0.8, 3., 2., 5.)
 569   return fInputEvent->IsPileupFromSPD((Int_t) fPileUpParamSPD[0] , fPileUpParamSPD[1] ,
 570                                               fPileUpParamSPD[2] , fPileUpParamSPD[3] , fPileUpParamSPD[4] );
 571   //printf("Param : %d, %2.2f, %2.2f, %2.2f, %2.2f\n",(Int_t) fPileUpParamSPD[0], fPileUpParamSPD[1], fPileUpParamSPD[2], fPileUpParamSPD[3], fPileUpParamSPD[4]);
 572
 573 }
 574
 575 //__________________________________________________
 576 Bool_t AliCaloTrackReader::IsPileUpFromEMCal() const
 577 {
 578   // Check if event is from pile-up determined by EMCal
 579   if(fNPileUpClusters > fNPileUpClustersCut) return kTRUE ;
 580   else                                       return kFALSE;
 581 }
 582
 583 //________________________________________________________
 584 Bool_t AliCaloTrackReader::IsPileUpFromSPDAndEMCal() const
 585 {
 586   // Check if event is from pile-up determined by SPD and EMCal
 587   if( IsPileUpFromSPD() && IsPileUpFromEMCal()) return kTRUE ;
 588   else                                          return kFALSE;
 589 }
 590
 591 //_______________________________________________________
 592 Bool_t AliCaloTrackReader::IsPileUpFromSPDOrEMCal() const
 593 {
 594   // Check if event is from pile-up determined by SPD or EMCal
 595   if( IsPileUpFromSPD() || IsPileUpFromEMCal()) return kTRUE ;
 596   else                                          return kFALSE;
 597 }
 598
 599 //___________________________________________________________
 600 Bool_t AliCaloTrackReader::IsPileUpFromSPDAndNotEMCal() const
 601 {
 602   // Check if event is from pile-up determined by SPD and not by EMCal
 603   if( IsPileUpFromSPD() && !IsPileUpFromEMCal()) return kTRUE ;
 604   else                                          return kFALSE;
 605 }
 606
 607 //___________________________________________________________
 608 Bool_t AliCaloTrackReader::IsPileUpFromEMCalAndNotSPD() const
 609 {
 610   // Check if event is from pile-up determined by EMCal, not by SPD
 611   if( !IsPileUpFromSPD() && IsPileUpFromEMCal()) return kTRUE ;
 612   else                                           return kFALSE;
 613 }
 614
 615 //______________________________________________________________
 616 Bool_t AliCaloTrackReader::IsPileUpFromNotSPDAndNotEMCal() const
 617 {
 618   // Check if event not from pile-up determined neither by SPD nor by EMCal
 619   if( !IsPileUpFromSPD() && !IsPileUpFromEMCal()) return kTRUE ;
 620   else                                            return kFALSE;
 621 }
 622
 623 //________________________________________________________
 624 void AliCaloTrackReader::Print(const Option_t * opt) const
 625 {
 626
 627   //Print some relevant parameters set for the analysis
 628   if(! opt)
 629     return;
 630
 631   printf("***** Print: %s %s ******\n",    GetName(), GetTitle() ) ;
 632   printf("Task name      : %s\n",          fTaskName.Data()) ;
 633   printf("Data type      : %d\n",          fDataType) ;
 634   printf("CTS Min pT     : %2.1f GeV/c\n", fCTSPtMin) ;
 635   printf("EMCAL Min pT   : %2.1f GeV/c\n", fEMCALPtMin) ;
 636   printf("PHOS Min pT    : %2.1f GeV/c\n", fPHOSPtMin) ;
 637   printf("CTS Max pT     : %2.1f GeV/c\n", fCTSPtMax) ;
 638   printf("EMCAL Max pT   : %2.1f GeV/c\n", fEMCALPtMax) ;
 639   printf("PHOS Max pT    : %2.1f GeV/c\n", fPHOSPtMax) ;
 640   printf("EMCAL Time Cut: %3.1f < TOF  < %3.1f\n", fEMCALTimeCutMin, fEMCALTimeCutMax);
 641   printf("Use CTS         =     %d\n",     fFillCTS) ;
 642   printf("Use EMCAL       =     %d\n",     fFillEMCAL) ;
 643   printf("Use PHOS        =     %d\n",     fFillPHOS) ;
 644   printf("Use EMCAL Cells =     %d\n",     fFillEMCALCells) ;
 645   printf("Use PHOS  Cells =     %d\n",     fFillPHOSCells) ;
 646   printf("Track status    =     %d\n", (Int_t) fTrackStatus) ;
 647   printf("AODs Track filter mask  =  %d or hybrid %d, SPD hit %d\n", (Int_t) fTrackFilterMask,fSelectHybridTracks,fSelectSPDHitTracks) ;
 648   printf("Track Mult Eta Cut =  %d\n", (Int_t) fTrackMultEtaCut) ;
 649   printf("Write delta AOD =     %d\n",     fWriteOutputDeltaAOD) ;
 650   printf("Recalculate Clusters = %d\n",    fRecalculateClusters) ;
 651
 652   printf("Use Triggers selected in SE base class %d; If not what trigger Mask? %d; Trigger max for mixed %d \n",
 653          fEventTriggerAtSE, fEventTriggerMask,fMixEventTriggerMask);
 654
 655   if(fComparePtHardAndClusterPt)
 656           printf("Compare jet pt and pt hard to accept event, factor = %2.2f",fPtHardAndJetPtFactor);
 657
 658   if(fComparePtHardAndClusterPt)
 659           printf("Compare cluster pt and pt hard to accept event, factor = %2.2f",fPtHardAndClusterPtFactor);
 660
 661   printf("Read Kine from, stack? %d, AOD ? %d \n", fReadStack, fReadAODMCParticles) ;
 662   printf("Delta AOD File Name =     %s\n", fDeltaAODFileName.Data()) ;
 663   printf("Centrality: Class %s, Option %d, Bin [%d,%d] \n", fCentralityClass.Data(),fCentralityOpt,fCentralityBin[0], fCentralityBin[1]) ;
 664
 665   printf("    \n") ;
 666
 667 }
 668
 669 //_________________________________________________________________________
 670 Bool_t AliCaloTrackReader::FillInputEvent(const Int_t iEntry,
 671                                           const char * /*currentFileName*/)
 672 {
 673   //Fill the event counter and input lists that are needed, called by the analysis maker.
 674
 675   fEventNumber = iEntry;
 676   //fCurrentFileName = TString(currentFileName);
 677   if(!fInputEvent)
 678   {
 679           if(fDebug >= 0) printf("AliCaloTrackReader::FillInputEvent() - Input event not available, skip event analysis\n");
 680           return kFALSE;
 681   }
 682
 683   //Select events only fired by a certain trigger configuration if it is provided
 684   Int_t eventType = 0;
 685   if(fInputEvent->GetHeader())
 686           eventType = ((AliVHeader*)fInputEvent->GetHeader())->GetEventType();
 687
 688   if (GetFiredTriggerClasses().Contains("FAST")  && !GetFiredTriggerClasses().Contains("ALL") && !fAcceptFastCluster)
 689   {
 690     if(fDebug > 0)  printf("AliCaloTrackReader::FillInputEvent - Do not count events from fast cluster, trigger name %s\n",fFiredTriggerClassName.Data());
 691     return kFALSE;
 692   }
 693
 694
 695   //---------------------------------------------------------------------
 696   // Do not count events that where likely triggered by an exotic cluster
 697   //---------------------------------------------------------------------
 698   if(fRemoveExoticEvents && GetFiredTriggerClasses().Contains("EMC"))
 699   {
 700     RejectExoticEvents();
 701     //printf("Is this event exotic? %d\n",fIsExoticEvent);
 702     if(fIsExoticEvent) return kFALSE;
 703   }
 704
 705   //-------------------------------------------------------------------------------------
 706   // Reject event if large clusters with large energy
 707   // Use only for LHC11a data for the moment, and if input is clusterizer V1 or V1+unfolding
 708   // If clusterzer NxN or V2 it does not help
 709   //-------------------------------------------------------------------------------------
 710   Int_t run = fInputEvent->GetRunNumber();
 711   if( fRemoveLEDEvents && run > 146857  && run < 146861 )
 712   {
 713     Bool_t reject = RejectLEDEvents();
 714     if(reject) return kFALSE;
 715   }// Remove LED events
 716
 717   //------------------------
 718   // Reject pure LED events?
 719   //-------------------------
 720   if( fFiredTriggerClassName  !="" && !fAnaLED)
 721   {
 722     //printf("Event type %d\n",eventType);
 723     if(eventType!=7)
 724       return kFALSE; //Only physics event, do not use for simulated events!!!
 725
 726     if(fDebug > 0)
 727       printf("AliCaloTrackReader::FillInputEvent() - FiredTriggerClass <%s>, selected class <%s>, compare name %d\n",
 728              GetFiredTriggerClasses().Data(),fFiredTriggerClassName.Data(), GetFiredTriggerClasses().Contains(fFiredTriggerClassName));
 729
 730     if( !GetFiredTriggerClasses().Contains(fFiredTriggerClassName) ) return kFALSE;
 731     else if(fDebug > 0) printf("AliCaloTrackReader::FillInputEvent() - Accepted triggered event\n");
 732   }
 733   else if(fAnaLED)
 734   {
 735     //    kStartOfRun =       1,    // START_OF_RUN
 736     //    kEndOfRun =         2,    // END_OF_RUN
 737     //    kStartOfRunFiles =  3,    // START_OF_RUN_FILES
 738     //    kEndOfRunFiles =    4,    // END_OF_RUN_FILES
 739     //    kStartOfBurst =     5,    // START_OF_BURST
 740     //    kEndOfBurst =       6,    // END_OF_BURST
 741     //    kPhysicsEvent =     7,    // PHYSICS_EVENT
 742     //    kCalibrationEvent = 8,    // CALIBRATION_EVENT
 743     //    kFormatError =      9,    // EVENT_FORMAT_ERROR
 744     //    kStartOfData =      10,   // START_OF_DATA
 745     //    kEndOfData =        11,   // END_OF_DATA
 746     //    kSystemSoftwareTriggerEvent   = 12, // SYSTEM_SOFTWARE_TRIGGER_EVENT
 747     //    kDetectorSoftwareTriggerEvent = 13  // DETECTOR_SOFTWARE_TRIGGER_EVENT
 748
 749           if(eventType!=7 && fDebug > 1 )printf("AliCaloTrackReader::FillInputEvent() - DO LED, Event Type <%d>, 8 Calibration \n",  eventType);
 750           if(eventType!=8)return kFALSE;
 751   }
 752
 753   //In case of analysis of events with jets, skip those with jet pt > 5 pt hard
 754   if(fComparePtHardAndJetPt)
 755   {
 756     if(!ComparePtHardAndJetPt()) return kFALSE ;
 757   }
 758
 759   if(fComparePtHardAndClusterPt)
 760   {
 761     if(!ComparePtHardAndClusterPt()) return kFALSE ;
 762   }
 763
 764   //Fill Vertex array
 765   FillVertexArray();
 766   //Reject events with Z vertex too large, only for SE analysis, if not, cut on the analysis code
 767   if(!GetMixedEvent() && TMath::Abs(fVertex[0][2]) > fZvtxCut) return kFALSE;
 768
 769   //------------------------------------------------------
 770   //Event rejection depending on vertex, pileup, v0and
 771   //------------------------------------------------------
 772   if(fDataType==kESD && fTimeStampEventSelect)
 773   {
 774     AliESDEvent* esd = dynamic_cast<AliESDEvent*> (fInputEvent);
 775     if(esd)
 776     {
 777       Int_t timeStamp = esd->GetTimeStamp();
 778       Float_t timeStampFrac = 1.*(timeStamp-fTimeStampRunMin) / (fTimeStampRunMax-fTimeStampRunMin);
 779
 780       //printf("stamp0 %d, max0 %d, frac %f\n", timeStamp-fTimeStampRunMin,fTimeStampRunMax-fTimeStampRunMin, timeStampFrac);
 781
 782       if(timeStampFrac < fTimeStampEventFracMin || timeStampFrac > fTimeStampEventFracMax) return kFALSE;
 783     }
 784     //printf("\t accept time stamp\n");
 785   }
 786
 787
 788   //------------------------------------------------------
 789   //Event rejection depending on vertex, pileup, v0and
 790   //------------------------------------------------------
 791
 792   if(fUseEventsWithPrimaryVertex)
 793   {
 794     if( !CheckForPrimaryVertex() )              return kFALSE;
 795     if( TMath::Abs(fVertex[0][0] ) < 1.e-6 &&
 796         TMath::Abs(fVertex[0][1] ) < 1.e-6 &&
 797         TMath::Abs(fVertex[0][2] ) < 1.e-6    ) return kFALSE;
 798   }
 799
 800   //printf("Reader : IsPileUp %d, Multi %d\n",IsPileUpFromSPD(),fInputEvent->IsPileupFromSPDInMultBins());
 801
 802   if(fDoEventSelection)
 803   {
 804     if(!fCaloFilterPatch)
 805     {
 806       // Do not analyze events with pileup
 807       Bool_t bPileup = IsPileUpFromSPD();
 808       //IsPileupFromSPDInMultBins() // method to try
 809       //printf("pile-up %d, %d, %2.2f, %2.2f, %2.2f, %2.2f\n",bPileup, (Int_t) fPileUpParamSPD[0], fPileUpParamSPD[1], fPileUpParamSPD[2], fPileUpParamSPD[3], fPileUpParamSPD[4]);
 810       if(bPileup) return kFALSE;
 811
 812       if(fDoV0ANDEventSelection)
 813       {
 814         Bool_t bV0AND = kTRUE;
 815         AliESDEvent* esd = dynamic_cast<AliESDEvent*> (fInputEvent);
 816         if(esd)
 817           bV0AND = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0AND);
 818         //else bV0AND = //FIXME FOR AODs
 819         if(!bV0AND) return kFALSE;
 820       }
 821     }//CaloFilter patch
 822     else
 823     {
 824       if(fInputEvent->GetNumberOfCaloClusters() > 0)
 825       {
 826         AliVCluster * calo = fInputEvent->GetCaloCluster(0);
 827         if(calo->GetNLabels() == 4)
 828         {
 829           Int_t * selection = calo->GetLabels();
 830           Bool_t bPileup = selection[0];
 831           if(bPileup) return kFALSE;
 832
 833           Bool_t bGoodV = selection[1];
 834           if(fUseEventsWithPrimaryVertex && !bGoodV) return kFALSE;
 835
 836           if(fDoV0ANDEventSelection)
 837           {
 838             Bool_t bV0AND = selection[2];
 839             if(!bV0AND) return kFALSE;
 840           }
 841
 842           fTrackMult = selection[3];
 843           if(fTrackMult == 0) return kFALSE;
 844         }
 845         else
 846         {
 847           //First filtered AODs, track multiplicity stored there.
 848           fTrackMult = (Int_t) ((AliAODHeader*)fInputEvent->GetHeader())->GetCentrality();
 849           if(fTrackMult == 0) return kFALSE;
 850         }
 851       }//at least one cluster
 852       else
 853       {
 854         //printf("AliCaloTrackReader::FillInputEvent() - No clusters in event\n");
 855         //Remove events with  vertex (0,0,0), bad vertex reconstruction
 856         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;
 857
 858         //First filtered AODs, track multiplicity stored there.
 859         fTrackMult = (Int_t) ((AliAODHeader*)fInputEvent->GetHeader())->GetCentrality();
 860         if(fTrackMult == 0) return kFALSE;
 861       }// no cluster
 862     }// CaloFileter patch
 863   }// Event selection/AliceSoft/AliRoot/trunk/PWG/CaloTrackCorrBase/AliCaloTrackReader.h
 864   //------------------------------------------------------
 865
 866   //Check if there is a centrality value, PbPb analysis, and if a centrality bin selection is requested
 867   //If we need a centrality bin, we select only those events in the corresponding bin.
 868   if(GetCentrality() && fCentralityBin[0]>=0 && fCentralityBin[1]>=0 && fCentralityOpt==100)
 869   {
 870     Int_t cen = GetEventCentrality();
 871     if(cen > fCentralityBin[1] || cen < fCentralityBin[0]) return kFALSE; //reject events out of bin.
 872   }
 873
 874   //Fill the arrays with cluster/tracks/cells data
 875
 876   if(!fEventTriggerAtSE)
 877   {
 878     // In case of mixing analysis, accept MB events, not only Trigger
 879     // Track and cluster arrays filled for MB in order to create the pool in the corresponding analysis
 880     // via de method in the base class FillMixedEventPool()
 881
 882     AliAnalysisManager *manager = AliAnalysisManager::GetAnalysisManager();
 883     AliInputEventHandler *inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
 884
 885     if(!inputHandler) return kFALSE ;  // to content coverity
 886
 887     UInt_t isTrigger = inputHandler->IsEventSelected() & fEventTriggerMask;
 888     UInt_t isMB      = inputHandler->IsEventSelected() & fMixEventTriggerMask;
 889
 890     if(!isTrigger && !isMB) return kFALSE;
 891
 892     //printf("Selected triggered event : %s\n",GetFiredTriggerClasses().Data());
 893   }
 894
 895   // Get the main vertex BC, in case not available
 896   // it is calculated in FillCTS checking the BC of tracks
 897   // with DCA small (if cut applied, if open)
 898   fVertexBC=fInputEvent->GetPrimaryVertex()->GetBC();
 899
 900   if(fFillCTS)
 901   {
 902     FillInputCTS();
 903     //Accept events with at least one track
 904     if(fTrackMult == 0 && fDoRejectNoTrackEvents) return kFALSE ;
 905   }
 906
 907   if(fDoVertexBCEventSelection)
 908   {
 909     if(fVertexBC!=0 && fVertexBC!=AliVTrack::kTOFBCNA) return kFALSE ;
 910   }
 911
 912   if(fFillEMCALCells)
 913     FillInputEMCALCells();
 914
 915   if(fFillPHOSCells)
 916     FillInputPHOSCells();
 917
 918   if(fFillEMCAL)
 919     FillInputEMCAL();
 920
 921   if(fFillPHOS)
 922     FillInputPHOS();
 923
 924   FillInputVZERO();
 925
 926
 927   return kTRUE ;
 928 }
 929
 930 //___________________________________
 931 void AliCaloTrackReader::ResetLists()
 932 {
 933   //  Reset lists, called by the analysis maker
 934
 935   if(fCTSTracks)       fCTSTracks     -> Clear();
 936   if(fEMCALClusters)   fEMCALClusters -> Clear("C");
 937   if(fPHOSClusters)    fPHOSClusters  -> Clear("C");
 938
 939   fV0ADC[0] = 0;   fV0ADC[1] = 0;
 940   fV0Mul[0] = 0;   fV0Mul[1] = 0;
 941
 942 }
 943
 944 //____________________________________________________________
 945 void AliCaloTrackReader::SetInputEvent(AliVEvent* const input)
 946 {
 947   fInputEvent  = input;
 948   fMixedEvent = dynamic_cast<AliMixedEvent*>(GetInputEvent()) ;
 949   if (fMixedEvent)
 950     fNMixedEvent = fMixedEvent->GetNumberOfEvents() ;
 951
 952   //Delete previous vertex
 953   if(fVertex)
 954   {
 955     for (Int_t i = 0; i < fNMixedEvent; i++)
 956     {
 957       delete [] fVertex[i] ;
 958     }
 959     delete [] fVertex ;
 960   }
 961
 962   fVertex = new Double_t*[fNMixedEvent] ;
 963   for (Int_t i = 0; i < fNMixedEvent; i++)
 964   {
 965     fVertex[i] = new Double_t[3] ;
 966     fVertex[i][0] = 0.0 ;
 967     fVertex[i][1] = 0.0 ;
 968     fVertex[i][2] = 0.0 ;
 969   }
 970 }
 971
 972 //__________________________________________________
 973 Int_t AliCaloTrackReader::GetEventCentrality() const
 974 {
 975   //Return current event centrality
 976
 977   if(GetCentrality())
 978   {
 979     if     (fCentralityOpt==100) return (Int_t) GetCentrality()->GetCentralityPercentile(fCentralityClass); // 100 bins max
 980     else if(fCentralityOpt==10)  return GetCentrality()->GetCentralityClass10(fCentralityClass);// 10 bins max
 981     else if(fCentralityOpt==20)  return GetCentrality()->GetCentralityClass5(fCentralityClass); // 20 bins max
 982     else
 983     {
 984       printf("AliCaloTrackReader::GetEventCentrality() - Unknown centrality option %d, use 10, 20 or 100\n",fCentralityOpt);
 985       return -1;
 986     }
 987   }
 988   else return -1;
 989
 990 }
 991
 992 //_____________________________________________________
 993 Double_t AliCaloTrackReader::GetEventPlaneAngle() const
 994 {
 995   //Return current event centrality
 996
 997   if(GetEventPlane())
 998   {
 999     Float_t ep =  GetEventPlane()->GetEventplane(GetEventPlaneMethod(), GetInputEvent());
1000
1001     if(GetEventPlaneMethod()=="Q" && (ep < 0 || ep > TMath::Pi()))
1002     {
1003       if(fDebug > 0 ) printf("AliCaloTrackReader::GetEventPlaneAngle() -  Bad EP for <Q> method : %f\n",ep);
1004       return -1000;
1005     }
1006     else if(GetEventPlaneMethod().Contains("V0")  )
1007     {
1008       if((ep > TMath::Pi()/2 || ep < -TMath::Pi()/2))
1009       {
1010         if(fDebug > 0 ) printf("AliCaloTrackReader::GetEventPlaneAngle() -  Bad EP for <%s> method : %f\n",GetEventPlaneMethod().Data(), ep);
1011         return -1000;
1012       }
1013
1014       ep+=TMath::Pi()/2; // put same range as for <Q> method
1015
1016     }
1017
1018     //printf("AliCaloTrackReader::GetEventPlaneAngle() = %f\n",ep);
1019     if(fDebug > 0 )
1020     {
1021       if     (ep > TMath::Pi()) printf("AliCaloTrackReader::GetEventPlaneAngle() - Too large angle = %f\n",ep);
1022       else if(ep < 0          ) printf("AliCaloTrackReader::GetEventPlaneAngle() - Negative angle = %f\n" ,ep);
1023     }
1024
1025     return ep;
1026   }
1027   else
1028   {
1029     if(fDataType!=kMC && fDebug > 0) printf("AliCaloTrackReader::GetEventPlaneAngle() -  No EP pointer\n");
1030     return -1000;
1031   }
1032
1033 }
1034
1035 //__________________________________________________________
1036 void AliCaloTrackReader::GetVertex(Double_t vertex[3]) const
1037 {
1038   //Return vertex position to be used for single event analysis
1039   vertex[0]=fVertex[0][0];
1040   vertex[1]=fVertex[0][1];
1041   vertex[2]=fVertex[0][2];
1042 }
1043
1044 //____________________________________________________________
1045 void AliCaloTrackReader::GetVertex(Double_t vertex[3],
1046                                    const Int_t evtIndex) const
1047 {
1048   //Return vertex position for mixed event, recover the vertex in a particular event.
1049
1050   vertex[0]=fVertex[evtIndex][0];  vertex[1]=fVertex[evtIndex][1];  vertex[2]=fVertex[evtIndex][2];
1051
1052 }
1053
1054 //________________________________________
1055 void AliCaloTrackReader::FillVertexArray()
1056 {
1057
1058   //Fill data member with vertex
1059   //In case of Mixed event, multiple vertices
1060
1061   //Delete previous vertex
1062   if(fVertex)
1063   {
1064     for (Int_t i = 0; i < fNMixedEvent; i++)
1065     {
1066       delete [] fVertex[i] ;
1067     }
1068     delete [] fVertex ;
1069   }
1070
1071   fVertex = new Double_t*[fNMixedEvent] ;
1072   for (Int_t i = 0; i < fNMixedEvent; i++)
1073   {
1074     fVertex[i] = new Double_t[3] ;
1075     fVertex[i][0] = 0.0 ;
1076     fVertex[i][1] = 0.0 ;
1077     fVertex[i][2] = 0.0 ;
1078   }
1079
1080   if (!fMixedEvent)
1081   { //Single event analysis
1082     if(fDataType!=kMC)
1083     {
1084
1085       if(fInputEvent->GetPrimaryVertex())
1086       {
1087         fInputEvent->GetPrimaryVertex()->GetXYZ(fVertex[0]);
1088       }
1089       else
1090       {
1091         printf("AliCaloTrackReader::FillVertexArray() - NULL primary vertex\n");
1092         fVertex[0][0]=0.;   fVertex[0][1]=0.;   fVertex[0][2]=0.;
1093       }//Primary vertex pointer do not exist
1094
1095     } else
1096     {//MC read event
1097       fVertex[0][0]=0.;   fVertex[0][1]=0.;   fVertex[0][2]=0.;
1098     }
1099
1100     if(fDebug > 1)
1101       printf("AliCaloTrackReader::FillVertexArray() - Single Event Vertex : %f,%f,%f\n",fVertex[0][0],fVertex[0][1],fVertex[0][2]);
1102
1103   } else
1104   { // MultiEvent analysis
1105     for (Int_t iev = 0; iev < fNMixedEvent; iev++)
1106     {
1107       if (fMixedEvent->GetVertexOfEvent(iev))
1108         fMixedEvent->GetVertexOfEvent(iev)->GetXYZ(fVertex[iev]);
1109       else
1110       { // no vertex found !!!!
1111         AliWarning("No vertex found");
1112       }
1113
1114       if(fDebug > 1)
1115         printf("AliCaloTrackReader::FillVertexArray() - Multi Event %d Vertex : %f,%f,%f\n",iev,fVertex[iev][0],fVertex[iev][1],fVertex[iev][2]);
1116
1117     }
1118   }
1119
1120 }
1121
1122 //_____________________________________
1123 void AliCaloTrackReader::FillInputCTS()
1124 {
1125   //Return array with Central Tracking System (CTS) tracks
1126
1127   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputCTS()\n");
1128
1129   Double_t pTrack[3] = {0,0,0};
1130
1131   Int_t nTracks = fInputEvent->GetNumberOfTracks() ;
1132   fTrackMult    = 0;
1133   Int_t nstatus = 0;
1134   Double_t bz   = GetInputEvent()->GetMagneticField();
1135
1136   for(Int_t i = 0; i < 19; i++)
1137   {
1138     fTrackBCEvent   [i] = 0;
1139     fTrackBCEventCut[i] = 0;
1140   }
1141
1142   Bool_t   bc0  = kFALSE;
1143   if(fRecalculateVertexBC) fVertexBC=AliVTrack::kTOFBCNA;
1144
1145   for (Int_t itrack =  0; itrack <  nTracks; itrack++)
1146   {////////////// track loop
1147     AliVTrack * track = (AliVTrack*)fInputEvent->GetTrack(itrack) ; // retrieve track from esd
1148
1149     //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
1150     ULong_t status = track->GetStatus();
1151
1152     if (fTrackStatus && !((status & fTrackStatus) == fTrackStatus))
1153       continue ;
1154
1155     nstatus++;
1156
1157     Float_t dcaTPC =-999;
1158
1159     if     (fDataType==kESD)
1160     {
1161       AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (track);
1162
1163       if(esdTrack)
1164       {
1165         if(fESDtrackCuts->AcceptTrack(esdTrack))
1166         {
1167           track->GetPxPyPz(pTrack) ;
1168
1169           if(fConstrainTrack)
1170           {
1171             if(esdTrack->GetConstrainedParam())
1172             {
1173               const AliExternalTrackParam* constrainParam = esdTrack->GetConstrainedParam();
1174               esdTrack->Set(constrainParam->GetX(),constrainParam->GetAlpha(),constrainParam->GetParameter(),constrainParam->GetCovariance());
1175               esdTrack->GetConstrainedPxPyPz(pTrack);
1176             }
1177             else continue;
1178
1179           } // use constrained tracks
1180
1181           if(fSelectSPDHitTracks)
1182           {//Not much sense to use with TPC only or Hybrid tracks
1183             if(!esdTrack->HasPointOnITSLayer(0) && !esdTrack->HasPointOnITSLayer(1)) continue ;
1184           }
1185         }
1186         // Complementary track to global : Hybrids (make sure that the previous selection is for Global)
1187         else  if(fESDtrackComplementaryCuts && fESDtrackComplementaryCuts->AcceptTrack(esdTrack))
1188         {
1189           // constrain the track
1190           if(esdTrack->GetConstrainedParam())
1191           {
1192             esdTrack->Set(esdTrack->GetConstrainedParam()->GetX(),esdTrack->GetConstrainedParam()->GetAlpha(),esdTrack->GetConstrainedParam()->GetParameter(),esdTrack->GetConstrainedParam()->GetCovariance());
1193
1194             track->GetPxPyPz(pTrack) ;
1195
1196           }
1197           else continue;
1198         }
1199         else continue;
1200       }
1201     } // ESD
1202     else if(fDataType==kAOD)
1203     {
1204       AliAODTrack *aodtrack = dynamic_cast <AliAODTrack*>(track);
1205
1206       if(aodtrack)
1207       {
1208        if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputCTS():AOD track type: %d (primary %d), hybrid? %d \n",
1209                               aodtrack->GetType(),AliAODTrack::kPrimary,
1210                               aodtrack->IsHybridGlobalConstrainedGlobal());
1211
1212
1213         if (fSelectHybridTracks)
1214         {
1215           if (!aodtrack->IsHybridGlobalConstrainedGlobal())       continue ;
1216         }
1217         else
1218         {
1219           if ( aodtrack->TestFilterBit(fTrackFilterMask)==kFALSE) continue ;
1220         }
1221
1222         if(fSelectSPDHitTracks)
1223         {//Not much sense to use with TPC only or Hybrid tracks
1224           if(!aodtrack->HasPointOnITSLayer(0) && !aodtrack->HasPointOnITSLayer(1)) continue ;
1225         }
1226
1227         if (aodtrack->GetType()!= AliAODTrack::kPrimary)          continue ;
1228
1229         if (fDebug > 2 ) printf("AliCaloTrackReader::FillInputCTS(): \t accepted track! \n");
1230
1231         //In case of AODs, TPC tracks cannot be propagated back to primary vertex,
1232         // info stored here
1233         dcaTPC = aodtrack->DCA();
1234
1235         track->GetPxPyPz(pTrack) ;
1236
1237       } // aod track exists
1238       else continue ;
1239
1240     } // AOD
1241
1242     TLorentzVector momentum(pTrack[0],pTrack[1],pTrack[2],0);
1243
1244     Bool_t okTOF  = ( (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ) ;
1245     Double_t tof  = -1000;
1246     Int_t trackBC = -1000 ;
1247
1248     if(okTOF)
1249     {
1250       trackBC = track->GetTOFBunchCrossing(bz);
1251       SetTrackEventBC(trackBC+9);
1252
1253       tof = track->GetTOFsignal()*1e-3;
1254     }
1255
1256     if(fUseTrackDCACut)
1257     {
1258       //normal way to get the dca, cut on dca_xy
1259       if(dcaTPC==-999)
1260       {
1261         Double_t dca[2]   = {1e6,1e6};
1262         Double_t covar[3] = {1e6,1e6,1e6};
1263         Bool_t okDCA = track->PropagateToDCA(fInputEvent->GetPrimaryVertex(),bz,100.,dca,covar);
1264         if( okDCA) okDCA = AcceptDCA(momentum.Pt(),dca[0]);
1265         if(!okDCA)
1266         {
1267           //printf("AliCaloTrackReader::FillInputCTS() - Reject track pt %2.2f, dca_xy %2.4f, BC %d\n",momentum.Pt(),dca[0],trackBC);
1268           continue ;
1269         }
1270       }
1271     }// DCA cuts
1272
1273     if(okTOF)
1274     {
1275       //SetTrackEventBCcut(bc);
1276       SetTrackEventBCcut(trackBC+9);
1277
1278       //After selecting tracks with small DCA, pointing to vertex, set vertex BC depeding on tracks BC
1279       if(fRecalculateVertexBC)
1280       {
1281         if     (trackBC !=0 && trackBC != AliVTrack::kTOFBCNA) fVertexBC = trackBC;
1282         else if(trackBC == 0)                                  bc0 = kTRUE;
1283       }
1284
1285       //In any case, the time should to be larger than the fixed window ...
1286       if( fUseTrackTimeCut && (trackBC!=0 || tof < fTrackTimeCutMin  || tof > fTrackTimeCutMax) )
1287       {
1288         //printf("Remove track time %f and bc = %d\n",tof,trackBC);
1289         continue ;
1290       }
1291       //else printf("Accept track time %f and bc = %d\n",tof,trackBC);
1292
1293     }
1294
1295     //Count the tracks in eta < 0.9
1296     //printf("Eta %f cut  %f\n",TMath::Abs(track->Eta()),fTrackMultEtaCut);
1297     if(TMath::Abs(track->Eta())< fTrackMultEtaCut) fTrackMult++;
1298
1299     if(fCTSPtMin > momentum.Pt() || fCTSPtMax < momentum.Pt()) continue ;
1300
1301     if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"CTS")) continue;
1302
1303     if(fDebug > 2 && momentum.Pt() > 0.1)
1304       printf("AliCaloTrackReader::FillInputCTS() - Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
1305              momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
1306
1307     if (fMixedEvent)  track->SetID(itrack);
1308
1309     fCTSTracks->Add(track);
1310
1311   }// track loop
1312
1313   if(fVertexBC ==0 || fVertexBC == AliVTrack::kTOFBCNA)
1314   {
1315     if( bc0 ) fVertexBC = 0 ;
1316     else      fVertexBC = AliVTrack::kTOFBCNA ;
1317   }
1318
1319   if(fDebug > 1)
1320     printf("AliCaloTrackReader::FillInputCTS()   - aod entries %d, input tracks %d, pass status %d, multipliticy %d\n", fCTSTracks->GetEntriesFast(), nTracks, nstatus, fTrackMult);//fCTSTracksNormalInputEntries);
1321
1322 }
1323
1324 //__________________________________________________________________
1325 void AliCaloTrackReader::FillInputEMCALAlgorithm(AliVCluster * clus,
1326                                                  const Int_t iclus)
1327 {
1328   //Fill the EMCAL data in the array, do it
1329
1330   Int_t vindex = 0 ;
1331   if (fMixedEvent)
1332     vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
1333
1334   //Reject clusters with bad channels, close to borders and exotic;
1335   if(!GetCaloUtils()->GetEMCALRecoUtils()->IsGoodCluster(clus,GetCaloUtils()->GetEMCALGeometry(),GetEMCALCells(),fInputEvent->GetBunchCrossNumber())) return;
1336
1337   //Mask all cells in collumns facing ALICE thick material if requested
1338   if(GetCaloUtils()->GetNMaskCellColumns())
1339   {
1340     Int_t absId   = -1;
1341     Int_t iSupMod = -1;
1342     Int_t iphi    = -1;
1343     Int_t ieta    = -1;
1344     Bool_t shared = kFALSE;
1345     GetCaloUtils()->GetEMCALRecoUtils()->GetMaxEnergyCell(GetCaloUtils()->GetEMCALGeometry(), GetEMCALCells(),clus,absId,iSupMod,ieta,iphi,shared);
1346     if(GetCaloUtils()->MaskFrameCluster(iSupMod, ieta)) return;
1347   }
1348
1349   if(fSelectEmbeddedClusters)
1350   {
1351     if(clus->GetNLabels()==0 || clus->GetLabel() < 0) return;
1352     //else printf("Embedded cluster,  %d, n label %d label %d  \n",iclus,clus->GetNLabels(),clus->GetLabel());
1353   }
1354
1355   //Float_t pos[3];
1356   //clus->GetPosition(pos);
1357   //printf("Before Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
1358
1359   if(fRecalculateClusters)
1360   {
1361     //Recalibrate the cluster energy
1362     if(GetCaloUtils()->IsRecalibrationOn())
1363     {
1364       Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, GetEMCALCells());
1365
1366       clus->SetE(energy);
1367       //printf("Recalibrated Energy %f\n",clus->E());
1368
1369       GetCaloUtils()->RecalculateClusterShowerShapeParameters(GetEMCALCells(),clus);
1370       GetCaloUtils()->RecalculateClusterPID(clus);
1371
1372     } // recalculate E
1373
1374     //Recalculate distance to bad channels, if new list of bad channels provided
1375     GetCaloUtils()->RecalculateClusterDistanceToBadChannel(GetEMCALCells(),clus);
1376
1377     //Recalculate cluster position
1378     if(GetCaloUtils()->IsRecalculationOfClusterPositionOn())
1379     {
1380       GetCaloUtils()->RecalculateClusterPosition(GetEMCALCells(),clus);
1381       //clus->GetPosition(pos);
1382       //printf("After  Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
1383     }
1384
1385     // Recalculate TOF
1386     if(GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn())
1387     {
1388       Double_t tof      = clus->GetTOF();
1389       Float_t  frac     =-1;
1390       Int_t    absIdMax = GetCaloUtils()->GetMaxEnergyCell(fEMCALCells, clus,frac);
1391
1392       if(fDataType==AliCaloTrackReader::kESD)
1393       {
1394         tof = fEMCALCells->GetCellTime(absIdMax);
1395       }
1396
1397       GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
1398
1399       clus->SetTOF(tof);
1400
1401     }// Time recalibration
1402   }
1403
1404   //Correct non linearity
1405   if(GetCaloUtils()->IsCorrectionOfClusterEnergyOn())
1406   {
1407     GetCaloUtils()->CorrectClusterEnergy(clus) ;
1408     //printf("Linearity Corrected Energy %f\n",clus->E());
1409
1410     //In case of MC analysis, to match resolution/calibration in real data
1411     Float_t rdmEnergy = GetCaloUtils()->GetEMCALRecoUtils()->SmearClusterEnergy(clus);
1412     // printf("\t Energy %f, smeared %f\n", clus->E(),rdmEnergy);
1413     clus->SetE(rdmEnergy);
1414   }
1415
1416   Double_t tof = clus->GetTOF()*1e9;
1417
1418   Int_t bc = TMath::Nint(tof/50) + 9;
1419   //printf("tof %2.2f, bc+5=%d\n",tof,bc);
1420
1421   SetEMCalEventBC(bc);
1422
1423   if(fEMCALPtMin > clus->E() || fEMCALPtMax < clus->E()) return ;
1424
1425   TLorentzVector momentum ;
1426
1427   clus->GetMomentum(momentum, fVertex[vindex]);
1428
1429   if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"EMCAL")) return ;
1430
1431   SetEMCalEventBCcut(bc);
1432
1433   if(!IsInTimeWindow(tof,clus->E()))
1434   {
1435     fNPileUpClusters++ ;
1436     if(fUseEMCALTimeCut) return ;
1437   }
1438   else
1439     fNNonPileUpClusters++;
1440
1441   if(fDebug > 2 && momentum.E() > 0.1)
1442     printf("AliCaloTrackReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
1443            momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
1444
1445   if (fMixedEvent)
1446     clus->SetID(iclus) ;
1447
1448   //Correct MC label for AODs
1449
1450   fEMCALClusters->Add(clus);
1451
1452 }
1453
1454 //_______________________________________
1455 void AliCaloTrackReader::FillInputEMCAL()
1456 {
1457   //Return array with EMCAL clusters in aod format
1458
1459   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputEMCAL()\n");
1460
1461   // First recalibrate cells, time or energy
1462   //  if(GetCaloUtils()->IsRecalibrationOn())
1463   //    GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCells(GetCaloUtils()->GetEMCALGeometry(),
1464   //                                                          GetEMCALCells(),
1465   //                                                          fInputEvent->GetBunchCrossNumber());
1466
1467   fNPileUpClusters    = 0; // Init counter
1468   fNNonPileUpClusters = 0; // Init counter
1469   for(Int_t i = 0; i < 19; i++)
1470   {
1471     fEMCalBCEvent   [i] = 0;
1472     fEMCalBCEventCut[i] = 0;
1473   }
1474
1475   //Loop to select clusters in fiducial cut and fill container with aodClusters
1476   if(fEMCALClustersListName=="")
1477   {
1478     Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
1479     for (Int_t iclus =  0; iclus <  nclusters; iclus++)
1480     {
1481       AliVCluster * clus = 0;
1482       if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
1483       {
1484         if (IsEMCALCluster(clus))
1485         {
1486           FillInputEMCALAlgorithm(clus, iclus);
1487         }//EMCAL cluster
1488       }// cluster exists
1489     }// cluster loop
1490
1491     //Recalculate track matching
1492     GetCaloUtils()->RecalculateClusterTrackMatching(fInputEvent);
1493
1494   }//Get the clusters from the input event
1495   else
1496   {
1497     TClonesArray * clusterList = 0x0;
1498
1499     if      (fInputEvent->FindListObject(fEMCALClustersListName))
1500     {
1501       clusterList = dynamic_cast<TClonesArray*> (fInputEvent->FindListObject(fEMCALClustersListName));
1502     }
1503     else if(fOutputEvent)
1504     {
1505       clusterList = dynamic_cast<TClonesArray*> (fOutputEvent->FindListObject(fEMCALClustersListName));
1506     }
1507
1508     if(!clusterList)
1509     {
1510         printf("AliCaloTrackReader::FillInputEMCAL() - Wrong name of list with clusters?  <%s>\n",fEMCALClustersListName.Data());
1511         return;
1512     }
1513
1514     Int_t nclusters = clusterList->GetEntriesFast();
1515     for (Int_t iclus =  0; iclus <  nclusters; iclus++)
1516     {
1517       AliVCluster * clus = dynamic_cast<AliVCluster*> (clusterList->At(iclus));
1518       //printf("E %f\n",clus->E());
1519       if (clus) FillInputEMCALAlgorithm(clus, iclus);
1520       else printf("AliCaloTrackReader::FillInputEMCAL() - Null cluster in list!\n");
1521     }// cluster loop
1522
1523     // Recalculate the pile-up time, in case long time clusters removed during clusterization
1524     //printf("Input event INIT : Pile-up clusters %d, NO pile-up %d\n",fNPileUpClusters,fNNonPileUpClusters);
1525
1526     fNPileUpClusters    = 0; // Init counter
1527     fNNonPileUpClusters = 0; // Init counter
1528     for(Int_t i = 0; i < 19; i++)
1529     {
1530       fEMCalBCEvent   [i] = 0;
1531       fEMCalBCEventCut[i] = 0;
1532     }
1533
1534     for (Int_t iclus =  0; iclus < fInputEvent->GetNumberOfCaloClusters(); iclus++)
1535     {
1536       AliVCluster * clus = 0;
1537
1538       if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
1539       {
1540         if (IsEMCALCluster(clus))
1541         {
1542
1543           Float_t  frac     =-1;
1544           Int_t    absIdMax = GetCaloUtils()->GetMaxEnergyCell(fEMCALCells, clus,frac);
1545           Double_t tof = clus->GetTOF();
1546           GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
1547           tof*=1e9;
1548
1549           //printf("Input event cluster : AbsIdMax %d, E %2.2f, time %2.2f \n", absIdMax,clus->E(),tof);
1550
1551           //Reject clusters with bad channels, close to borders and exotic;
1552           if(!GetCaloUtils()->GetEMCALRecoUtils()->IsGoodCluster(clus,GetCaloUtils()->GetEMCALGeometry(),GetEMCALCells(),fInputEvent->GetBunchCrossNumber()))  continue;
1553
1554           Int_t bc = TMath::Nint(tof/50) + 9;
1555           SetEMCalEventBC(bc);
1556
1557           if(fEMCALPtMin > clus->E() || fEMCALPtMax < clus->E()) continue ;
1558
1559           TLorentzVector momentum ;
1560
1561           clus->GetMomentum(momentum, fVertex[0]);
1562
1563           if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"EMCAL")) return ;
1564
1565           SetEMCalEventBCcut(bc);
1566
1567           if(!IsInTimeWindow(tof,clus->E()))
1568             fNPileUpClusters++ ;
1569           else
1570             fNNonPileUpClusters++;
1571
1572         }
1573       }
1574     }
1575
1576     //printf("Input event : Pile-up clusters %d, NO pile-up %d\n",fNPileUpClusters,fNNonPileUpClusters);
1577
1578     // Recalculate track matching, not necessary if already done in the reclusterization task.
1579     // in case it was not done ...
1580     GetCaloUtils()->RecalculateClusterTrackMatching(fInputEvent,clusterList);
1581
1582   }
1583
1584   if(fDebug > 1) printf("AliCaloTrackReader::FillInputEMCAL() - aod entries %d, n pile-up clusters %d, n non pile-up %d \n",  fEMCALClusters->GetEntriesFast(),fNPileUpClusters,fNNonPileUpClusters);
1585
1586 }
1587
1588 //______________________________________
1589 void AliCaloTrackReader::FillInputPHOS()
1590 {
1591   //Return array with PHOS clusters in aod format
1592
1593   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputPHOS()\n");
1594
1595   //Loop to select clusters in fiducial cut and fill container with aodClusters
1596   Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
1597   for (Int_t iclus = 0; iclus < nclusters; iclus++)
1598   {
1599     AliVCluster * clus = 0;
1600     if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
1601     {
1602       if (IsPHOSCluster(clus))
1603       {
1604         //Check if the cluster contains any bad channel and if close to calorimeter borders
1605         Int_t vindex = 0 ;
1606         if (fMixedEvent)
1607           vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
1608         if( GetCaloUtils()->ClusterContainsBadChannel("PHOS",clus->GetCellsAbsId(), clus->GetNCells()))
1609           continue;
1610         if(!GetCaloUtils()->CheckCellFiducialRegion(clus, fInputEvent->GetPHOSCells(), fInputEvent, vindex))
1611           continue;
1612
1613         if(fRecalculateClusters)
1614         {
1615           //Recalibrate the cluster energy
1616           if(GetCaloUtils()->IsRecalibrationOn())
1617           {
1618             Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliAODCaloCells*)GetPHOSCells());
1619             clus->SetE(energy);
1620           }
1621         }
1622
1623         TLorentzVector momentum ;
1624
1625         clus->GetMomentum(momentum, fVertex[vindex]);
1626
1627         if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"PHOS")) continue;
1628
1629         if(fPHOSPtMin > momentum.E() || fPHOSPtMax < momentum.E())          continue;
1630
1631         if(fDebug > 2 && momentum.E() > 0.1)
1632           printf("AliCaloTrackReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
1633                  momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
1634
1635
1636         if (fMixedEvent)
1637         {
1638           clus->SetID(iclus) ;
1639         }
1640
1641         fPHOSClusters->Add(clus);
1642
1643       }//PHOS cluster
1644     }//cluster exists
1645   }//esd cluster loop
1646
1647   if(fDebug > 1) printf("AliCaloTrackReader::FillInputPHOS()  - aod entries %d\n",  fPHOSClusters->GetEntriesFast());
1648
1649 }
1650
1651 //____________________________________________
1652 void AliCaloTrackReader::FillInputEMCALCells()
1653 {
1654   //Return array with EMCAL cells in aod format
1655
1656   fEMCALCells = fInputEvent->GetEMCALCells();
1657
1658 }
1659
1660 //___________________________________________
1661 void AliCaloTrackReader::FillInputPHOSCells()
1662 {
1663   //Return array with PHOS cells in aod format
1664
1665   fPHOSCells = fInputEvent->GetPHOSCells();
1666
1667 }
1668
1669 //_______________________________________
1670 void AliCaloTrackReader::FillInputVZERO()
1671 {
1672   //Fill VZERO information in data member, add all the channels information.
1673   AliVVZERO* v0 = fInputEvent->GetVZEROData();
1674   //printf("Init V0: ADC (%d,%d), Multiplicity (%d,%d) \n",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]);
1675
1676   if (v0)
1677   {
1678     AliESDVZERO* esdV0 = dynamic_cast<AliESDVZERO*> (v0);
1679     for (Int_t i = 0; i < 32; i++)
1680     {
1681       if(esdV0)
1682       {//Only available in ESDs
1683         fV0ADC[0] += (Int_t)esdV0->GetAdcV0C(i);
1684         fV0ADC[1] += (Int_t)esdV0->GetAdcV0A(i);
1685       }
1686
1687       fV0Mul[0] += (Int_t)v0->GetMultiplicityV0C(i);
1688       fV0Mul[1] += (Int_t)v0->GetMultiplicityV0A(i);
1689     }
1690     if(fDebug > 0)
1691       printf("V0: ADC (%d,%d), Multiplicity (%d,%d) \n",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]);
1692   }
1693   else
1694   {
1695     if(fDebug > 0)
1696       printf("Cannot retrieve V0 ESD! Run w/ null V0 charges\n ");
1697   }
1698 }
1699
1700
1701 //___________________________________________________________________
1702 Bool_t AliCaloTrackReader::IsEMCALCluster(AliVCluster* cluster) const
1703 {
1704   // Check if it is a cluster from EMCAL. For old AODs cluster type has
1705   // different number and need to patch here
1706
1707   if(fDataType==kAOD && fOldAOD)
1708   {
1709     if (cluster->GetType() == 2) return kTRUE;
1710     else                         return kFALSE;
1711   }
1712   else
1713   {
1714     return cluster->IsEMCAL();
1715   }
1716
1717 }
1718
1719 //___________________________________________________________________
1720 Bool_t AliCaloTrackReader::IsPHOSCluster(AliVCluster * cluster) const
1721 {
1722   //Check if it is a cluster from PHOS.For old AODs cluster type has
1723   // different number and need to patch here
1724
1725   if(fDataType==kAOD && fOldAOD)
1726   {
1727     Int_t type = cluster->GetType();
1728     if (type == 0 || type == 1) return kTRUE;
1729     else                        return kFALSE;
1730   }
1731   else
1732   {
1733     return cluster->IsPHOS();
1734   }
1735
1736 }
1737
1738 //________________________________________________
1739 Bool_t AliCaloTrackReader::CheckForPrimaryVertex()
1740 {
1741   //Check if the vertex was well reconstructed, copy from V0Reader of conversion group
1742   //Only for ESDs ...
1743
1744   AliESDEvent * event = dynamic_cast<AliESDEvent*> (fInputEvent);
1745   if(!event) return kTRUE;
1746
1747   if(event->GetPrimaryVertexTracks()->GetNContributors() > 0)
1748   {
1749     return kTRUE;
1750   }
1751
1752   if(event->GetPrimaryVertexTracks()->GetNContributors() < 1)
1753   {
1754     // SPD vertex
1755     if(event->GetPrimaryVertexSPD()->GetNContributors() > 0)
1756     {
1757       //cout<<"spd vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
1758       return kTRUE;
1759
1760     }
1761     if(event->GetPrimaryVertexSPD()->GetNContributors() < 1)
1762     {
1763       //      cout<<"bad vertex type::"<< event->GetPrimaryVertex()->GetName() << endl;
1764       return kFALSE;
1765     }
1766   }
1767
1768   return kFALSE;
1769
1770 }
1771
1772 //____________________________________________
1773 void  AliCaloTrackReader::RejectExoticEvents()
1774 {
1775   // Reject events triggered by exotic cells
1776   // simple method
1777   Int_t nOfHighECl = 0 ;
1778   Int_t nExo       = 0 ;
1779   Int_t nclusters  = fInputEvent->GetNumberOfCaloClusters();
1780   for (Int_t iclus =  0; iclus <  nclusters; iclus++)
1781   {
1782     AliVCluster * clus = 0;
1783     if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
1784     {
1785       if (IsEMCALCluster(clus) && clus->E() > fExoticTrigger)
1786       {
1787         Bool_t exotic = GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCluster(clus, fInputEvent->GetEMCALCells());
1788         if(exotic) nExo++;
1789         else       nOfHighECl++;
1790       }//EMCAL cluster
1791     }// cluster exists
1792   }// cluster loop
1793
1794   //printf("trigger %2.1f, n Exotic %d, n high energy %d\n", fExoticTrigger,nExo,nOfHighECl);
1795
1796   if ((nExo>0)&&(nOfHighECl==0)) fIsExoticEvent = kTRUE ;
1797   else                           fIsExoticEvent = kFALSE;
1798
1799 }
1800
1801 //__________________________________________
1802 Bool_t  AliCaloTrackReader::RejectLEDEvents()
1803 {
1804   // LED Events in period LHC11a contaminated sample, simple method
1805   // to reject such events
1806
1807   // Count number of cells with energy larger than 0.1 in SM3, cut on this number
1808   Int_t ncellsSM3 = 0;
1809   for(Int_t icell = 0; icell < fInputEvent->GetEMCALCells()->GetNumberOfCells(); icell++)
1810   {
1811     Int_t absID = fInputEvent->GetEMCALCells()->GetCellNumber(icell);
1812     Int_t sm    = GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absID);
1813     if(fInputEvent->GetEMCALCells()->GetAmplitude(icell) > 0.1 && sm==3) ncellsSM3++;
1814   }
1815
1816   Int_t ncellcut = 21;
1817   if(GetFiredTriggerClasses().Contains("EMC")) ncellcut = 35;
1818
1819   if(ncellsSM3 >= ncellcut)
1820   {
1821     if(fDebug > 0)
1822       printf(" AliCaloTrackReader::FillInputEvent() - reject event with ncells in SM3 %d, cut %d, trig %s\n",
1823              ncellsSM3,ncellcut,GetFiredTriggerClasses().Data());
1824     return kTRUE;
1825   }
1826
1827   return kFALSE;
1828
1829 }
1830
1831 //____________________________________________________________
1832 void  AliCaloTrackReader::SetTrackCuts(AliESDtrackCuts * cuts)
1833 {
1834   // Set Track cuts
1835
1836   if(fESDtrackCuts) delete fESDtrackCuts ;
1837
1838   fESDtrackCuts = cuts ;
1839
1840 }
1841
1842 //_________________________________________________________________________
1843 void  AliCaloTrackReader::SetTrackComplementaryCuts(AliESDtrackCuts * cuts)
1844 {
1845   // Set Track cuts for complementary tracks (hybrids)
1846
1847   if(fESDtrackComplementaryCuts) delete fESDtrackComplementaryCuts ;
1848
1849   fESDtrackComplementaryCuts = cuts ;
1850
1851 }
1852
1853