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