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),
 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
 698   //fCurrentFileName = TString(currentFileName);
 699   if(!fInputEvent)
 700   {
 701           if(fDebug >= 0) printf("AliCaloTrackReader::FillInputEvent() - Input event not available, skip event analysis\n");
 702           return kFALSE;
 703   }
 704
 705   //Select events only fired by a certain trigger configuration if it is provided
 706   Int_t eventType = 0;
 707   if(fInputEvent->GetHeader())
 708           eventType = ((AliVHeader*)fInputEvent->GetHeader())->GetEventType();
 709
 710   if (GetFiredTriggerClasses().Contains("FAST")  && !GetFiredTriggerClasses().Contains("ALL") && !fAcceptFastCluster)
 711   {
 712     if(fDebug > 0)  printf("AliCaloTrackReader::FillInputEvent - Do not count events from fast cluster, trigger name %s\n",fFiredTriggerClassName.Data());
 713     return kFALSE;
 714   }
 715
 716
 717   //-------------------------------------------------------------------------------------
 718   // Reject event if large clusters with large energy
 719   // Use only for LHC11a data for the moment, and if input is clusterizer V1 or V1+unfolding
 720   // If clusterzer NxN or V2 it does not help
 721   //-------------------------------------------------------------------------------------
 722   Int_t run = fInputEvent->GetRunNumber();
 723   if( fRemoveLEDEvents && run > 146857  && run < 146861 )
 724   {
 725     Bool_t reject = RejectLEDEvents();
 726     if(reject) return kFALSE;
 727   }// Remove LED events
 728
 729   //------------------------
 730   // Reject pure LED events?
 731   //-------------------------
 732   if( fFiredTriggerClassName  !="" && !fAnaLED)
 733   {
 734     //printf("Event type %d\n",eventType);
 735     if(eventType!=7)
 736       return kFALSE; //Only physics event, do not use for simulated events!!!
 737
 738     if(fDebug > 0)
 739       printf("AliCaloTrackReader::FillInputEvent() - FiredTriggerClass <%s>, selected class <%s>, compare name %d\n",
 740              GetFiredTriggerClasses().Data(),fFiredTriggerClassName.Data(), GetFiredTriggerClasses().Contains(fFiredTriggerClassName));
 741
 742     if( !GetFiredTriggerClasses().Contains(fFiredTriggerClassName) ) return kFALSE;
 743     else if(fDebug > 0) printf("AliCaloTrackReader::FillInputEvent() - Accepted triggered event\n");
 744   }
 745   else if(fAnaLED)
 746   {
 747     //    kStartOfRun =       1,    // START_OF_RUN
 748     //    kEndOfRun =         2,    // END_OF_RUN
 749     //    kStartOfRunFiles =  3,    // START_OF_RUN_FILES
 750     //    kEndOfRunFiles =    4,    // END_OF_RUN_FILES
 751     //    kStartOfBurst =     5,    // START_OF_BURST
 752     //    kEndOfBurst =       6,    // END_OF_BURST
 753     //    kPhysicsEvent =     7,    // PHYSICS_EVENT
 754     //    kCalibrationEvent = 8,    // CALIBRATION_EVENT
 755     //    kFormatError =      9,    // EVENT_FORMAT_ERROR
 756     //    kStartOfData =      10,   // START_OF_DATA
 757     //    kEndOfData =        11,   // END_OF_DATA
 758     //    kSystemSoftwareTriggerEvent   = 12, // SYSTEM_SOFTWARE_TRIGGER_EVENT
 759     //    kDetectorSoftwareTriggerEvent = 13  // DETECTOR_SOFTWARE_TRIGGER_EVENT
 760
 761           if(eventType!=7 && fDebug > 1 )printf("AliCaloTrackReader::FillInputEvent() - DO LED, Event Type <%d>, 8 Calibration \n",  eventType);
 762           if(eventType!=8)return kFALSE;
 763   }
 764
 765   //In case of analysis of events with jets, skip those with jet pt > 5 pt hard
 766   if(fComparePtHardAndJetPt)
 767   {
 768     if(!ComparePtHardAndJetPt()) return kFALSE ;
 769   }
 770
 771   if(fComparePtHardAndClusterPt)
 772   {
 773     if(!ComparePtHardAndClusterPt()) return kFALSE ;
 774   }
 775
 776   //Fill Vertex array
 777   FillVertexArray();
 778   //Reject events with Z vertex too large, only for SE analysis, if not, cut on the analysis code
 779   if(!GetMixedEvent() && TMath::Abs(fVertex[0][2]) > fZvtxCut) return kFALSE;
 780
 781   //------------------------------------------------------
 782   //Event rejection depending on vertex, pileup, v0and
 783   //------------------------------------------------------
 784   if(fDataType==kESD && fTimeStampEventSelect)
 785   {
 786     AliESDEvent* esd = dynamic_cast<AliESDEvent*> (fInputEvent);
 787     if(esd)
 788     {
 789       Int_t timeStamp = esd->GetTimeStamp();
 790       Float_t timeStampFrac = 1.*(timeStamp-fTimeStampRunMin) / (fTimeStampRunMax-fTimeStampRunMin);
 791
 792       //printf("stamp0 %d, max0 %d, frac %f\n", timeStamp-fTimeStampRunMin,fTimeStampRunMax-fTimeStampRunMin, timeStampFrac);
 793
 794       if(timeStampFrac < fTimeStampEventFracMin || timeStampFrac > fTimeStampEventFracMax) return kFALSE;
 795     }
 796     //printf("\t accept time stamp\n");
 797   }
 798
 799
 800   //------------------------------------------------------
 801   //Event rejection depending on vertex, pileup, v0and
 802   //------------------------------------------------------
 803
 804   if(fUseEventsWithPrimaryVertex)
 805   {
 806     if( !CheckForPrimaryVertex() )              return kFALSE;
 807     if( TMath::Abs(fVertex[0][0] ) < 1.e-6 &&
 808         TMath::Abs(fVertex[0][1] ) < 1.e-6 &&
 809         TMath::Abs(fVertex[0][2] ) < 1.e-6    ) return kFALSE;
 810   }
 811
 812   //printf("Reader : IsPileUp %d, Multi %d\n",IsPileUpFromSPD(),fInputEvent->IsPileupFromSPDInMultBins());
 813
 814   if(fDoEventSelection)
 815   {
 816     if(!fCaloFilterPatch)
 817     {
 818       // Do not analyze events with pileup
 819       Bool_t bPileup = IsPileUpFromSPD();
 820       //IsPileupFromSPDInMultBins() // method to try
 821       //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]);
 822       if(bPileup) return kFALSE;
 823
 824       if(fDoV0ANDEventSelection)
 825       {
 826         Bool_t bV0AND = kTRUE;
 827         AliESDEvent* esd = dynamic_cast<AliESDEvent*> (fInputEvent);
 828         if(esd)
 829           bV0AND = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0AND);
 830         //else bV0AND = //FIXME FOR AODs
 831         if(!bV0AND) return kFALSE;
 832       }
 833     }//CaloFilter patch
 834     else
 835     {
 836       if(fInputEvent->GetNumberOfCaloClusters() > 0)
 837       {
 838         AliVCluster * calo = fInputEvent->GetCaloCluster(0);
 839         if(calo->GetNLabels() == 4)
 840         {
 841           Int_t * selection = calo->GetLabels();
 842           Bool_t bPileup = selection[0];
 843           if(bPileup) return kFALSE;
 844
 845           Bool_t bGoodV = selection[1];
 846           if(fUseEventsWithPrimaryVertex && !bGoodV) return kFALSE;
 847
 848           if(fDoV0ANDEventSelection)
 849           {
 850             Bool_t bV0AND = selection[2];
 851             if(!bV0AND) return kFALSE;
 852           }
 853
 854           fTrackMult = selection[3];
 855           if(fTrackMult == 0) return kFALSE;
 856         }
 857         else
 858         {
 859           //First filtered AODs, track multiplicity stored there.
 860           fTrackMult = (Int_t) ((AliAODHeader*)fInputEvent->GetHeader())->GetCentrality();
 861           if(fTrackMult == 0) return kFALSE;
 862         }
 863       }//at least one cluster
 864       else
 865       {
 866         //printf("AliCaloTrackReader::FillInputEvent() - No clusters in event\n");
 867         //Remove events with  vertex (0,0,0), bad vertex reconstruction
 868         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;
 869
 870         //First filtered AODs, track multiplicity stored there.
 871         fTrackMult = (Int_t) ((AliAODHeader*)fInputEvent->GetHeader())->GetCentrality();
 872         if(fTrackMult == 0) return kFALSE;
 873       }// no cluster
 874     }// CaloFileter patch
 875   }// Event selection/AliceSoft/AliRoot/trunk/PWG/CaloTrackCorrBase/AliCaloTrackReader.h
 876
 877   //------------------------------------------------------
 878
 879   //Check if there is a centrality value, PbPb analysis, and if a centrality bin selection is requested
 880   //If we need a centrality bin, we select only those events in the corresponding bin.
 881   if(GetCentrality() && fCentralityBin[0]>=0 && fCentralityBin[1]>=0 && fCentralityOpt==100)
 882   {
 883     Int_t cen = GetEventCentrality();
 884     if(cen > fCentralityBin[1] || cen < fCentralityBin[0]) return kFALSE; //reject events out of bin.
 885   }
 886
 887   //Fill the arrays with cluster/tracks/cells data
 888
 889   if(!fEventTriggerAtSE)
 890   {
 891     // In case of mixing analysis, accept MB events, not only Trigger
 892     // Track and cluster arrays filled for MB in order to create the pool in the corresponding analysis
 893     // via de method in the base class FillMixedEventPool()
 894
 895     AliAnalysisManager *manager = AliAnalysisManager::GetAnalysisManager();
 896     AliInputEventHandler *inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
 897
 898     if(!inputHandler) return kFALSE ;  // to content coverity
 899
 900     UInt_t isTrigger = inputHandler->IsEventSelected() & fEventTriggerMask;
 901     UInt_t isMB      = inputHandler->IsEventSelected() & fMixEventTriggerMask;
 902
 903     if(!isTrigger && !isMB) return kFALSE;
 904
 905     //printf("Selected triggered event : %s\n",GetFiredTriggerClasses().Data());
 906   }
 907
 908   //----------------------------------------------------------------------
 909   // Do not count events that where likely triggered by an exotic cluster
 910   // or out BC cluster
 911   //----------------------------------------------------------------------
 912
 913   //Get Patches that triggered
 914   TArrayI patches = GetL0TriggerPatches();
 915 /*
 916   if(fRemoveExoticEvents)
 917   {
 918     RejectExoticEvents(patches);
 919     if(fIsExoticEvent)
 920     {
 921       //printf("AliCaloTrackReader::FillInputEvent() - REJECT exotic triggered event \n");
 922       return kFALSE;
 923     }
 924   }
 925
 926   RejectTriggeredEventsByPileUp(patches);
 927   //printf("AliCaloTrackReader::FillInputEvent(), Trigger BC = %d\n",fTriggerClusterBC);
 928
 929   if(fRemoveTriggerOutBCEvents)
 930   {
 931     if(fTriggerClusterBC != 0 && fTriggerClusterBC != 6)
 932     {
 933       //printf("\t REJECT, bad trigger cluster BC\n");
 934       return kFALSE;
 935     }
 936   }
 937 */
 938
 939   MatchTriggerCluster(patches);
 940
 941   if(fRemoveBadTriggerEvents)
 942   {
 943     printf("ACCEPT triggered event? - exotic? %d - bad cell %d - BC %d  - Matched %d\n",fIsExoticEvent,fIsBadCellEvent,fTriggerClusterBC,fIsTriggerMatch);
 944     if     (fIsExoticEvent)         return kFALSE;
 945     else if(fIsBadCellEvent)        return kFALSE;
 946     else if(fTriggerClusterBC == 0) return kFALSE;
 947     printf("\t *** YES\n");
 948   }
 949
 950   patches.Reset();
 951
 952   // Get the main vertex BC, in case not available
 953   // it is calculated in FillCTS checking the BC of tracks
 954   // with DCA small (if cut applied, if open)
 955   fVertexBC=fInputEvent->GetPrimaryVertex()->GetBC();
 956
 957   if(fFillCTS)
 958   {
 959     FillInputCTS();
 960     //Accept events with at least one track
 961     if(fTrackMult == 0 && fDoRejectNoTrackEvents) return kFALSE ;
 962   }
 963
 964   if(fDoVertexBCEventSelection)
 965   {
 966     if(fVertexBC!=0 && fVertexBC!=AliVTrack::kTOFBCNA) return kFALSE ;
 967   }
 968
 969   if(fFillEMCALCells)
 970     FillInputEMCALCells();
 971
 972   if(fFillPHOSCells)
 973     FillInputPHOSCells();
 974
 975   if(fFillEMCAL)
 976     FillInputEMCAL();
 977
 978   if(fFillPHOS)
 979     FillInputPHOS();
 980
 981   FillInputVZERO();
 982
 983
 984   return kTRUE ;
 985 }
 986
 987 //___________________________________
 988 void AliCaloTrackReader::ResetLists()
 989 {
 990   //  Reset lists, called by the analysis maker
 991
 992   if(fCTSTracks)       fCTSTracks     -> Clear();
 993   if(fEMCALClusters)   fEMCALClusters -> Clear("C");
 994   if(fPHOSClusters)    fPHOSClusters  -> Clear("C");
 995
 996   fV0ADC[0] = 0;   fV0ADC[1] = 0;
 997   fV0Mul[0] = 0;   fV0Mul[1] = 0;
 998
 999 }
1000
1001 //____________________________________________________________
1002 void AliCaloTrackReader::SetInputEvent(AliVEvent* const input)
1003 {
1004   fInputEvent  = input;
1005   fMixedEvent = dynamic_cast<AliMixedEvent*>(GetInputEvent()) ;
1006   if (fMixedEvent)
1007     fNMixedEvent = fMixedEvent->GetNumberOfEvents() ;
1008
1009   //Delete previous vertex
1010   if(fVertex)
1011   {
1012     for (Int_t i = 0; i < fNMixedEvent; i++)
1013     {
1014       delete [] fVertex[i] ;
1015     }
1016     delete [] fVertex ;
1017   }
1018
1019   fVertex = new Double_t*[fNMixedEvent] ;
1020   for (Int_t i = 0; i < fNMixedEvent; i++)
1021   {
1022     fVertex[i] = new Double_t[3] ;
1023     fVertex[i][0] = 0.0 ;
1024     fVertex[i][1] = 0.0 ;
1025     fVertex[i][2] = 0.0 ;
1026   }
1027 }
1028
1029 //__________________________________________________
1030 Int_t AliCaloTrackReader::GetEventCentrality() const
1031 {
1032   //Return current event centrality
1033
1034   if(GetCentrality())
1035   {
1036     if     (fCentralityOpt==100) return (Int_t) GetCentrality()->GetCentralityPercentile(fCentralityClass); // 100 bins max
1037     else if(fCentralityOpt==10)  return GetCentrality()->GetCentralityClass10(fCentralityClass);// 10 bins max
1038     else if(fCentralityOpt==20)  return GetCentrality()->GetCentralityClass5(fCentralityClass); // 20 bins max
1039     else
1040     {
1041       printf("AliCaloTrackReader::GetEventCentrality() - Unknown centrality option %d, use 10, 20 or 100\n",fCentralityOpt);
1042       return -1;
1043     }
1044   }
1045   else return -1;
1046
1047 }
1048
1049 //_____________________________________________________
1050 Double_t AliCaloTrackReader::GetEventPlaneAngle() const
1051 {
1052   //Return current event centrality
1053
1054   if(GetEventPlane())
1055   {
1056     Float_t ep =  GetEventPlane()->GetEventplane(GetEventPlaneMethod(), GetInputEvent());
1057
1058     if(GetEventPlaneMethod()=="Q" && (ep < 0 || ep > TMath::Pi()))
1059     {
1060       if(fDebug > 0 ) printf("AliCaloTrackReader::GetEventPlaneAngle() -  Bad EP for <Q> method : %f\n",ep);
1061       return -1000;
1062     }
1063     else if(GetEventPlaneMethod().Contains("V0")  )
1064     {
1065       if((ep > TMath::Pi()/2 || ep < -TMath::Pi()/2))
1066       {
1067         if(fDebug > 0 ) printf("AliCaloTrackReader::GetEventPlaneAngle() -  Bad EP for <%s> method : %f\n",GetEventPlaneMethod().Data(), ep);
1068         return -1000;
1069       }
1070
1071       ep+=TMath::Pi()/2; // put same range as for <Q> method
1072
1073     }
1074
1075     //printf("AliCaloTrackReader::GetEventPlaneAngle() = %f\n",ep);
1076     if(fDebug > 0 )
1077     {
1078       if     (ep > TMath::Pi()) printf("AliCaloTrackReader::GetEventPlaneAngle() - Too large angle = %f\n",ep);
1079       else if(ep < 0          ) printf("AliCaloTrackReader::GetEventPlaneAngle() - Negative angle = %f\n" ,ep);
1080     }
1081
1082     return ep;
1083   }
1084   else
1085   {
1086     if(fDataType!=kMC && fDebug > 0) printf("AliCaloTrackReader::GetEventPlaneAngle() -  No EP pointer\n");
1087     return -1000;
1088   }
1089
1090 }
1091
1092 //__________________________________________________________
1093 void AliCaloTrackReader::GetVertex(Double_t vertex[3]) const
1094 {
1095   //Return vertex position to be used for single event analysis
1096   vertex[0]=fVertex[0][0];
1097   vertex[1]=fVertex[0][1];
1098   vertex[2]=fVertex[0][2];
1099 }
1100
1101 //____________________________________________________________
1102 void AliCaloTrackReader::GetVertex(Double_t vertex[3],
1103                                    const Int_t evtIndex) const
1104 {
1105   //Return vertex position for mixed event, recover the vertex in a particular event.
1106
1107   vertex[0]=fVertex[evtIndex][0];  vertex[1]=fVertex[evtIndex][1];  vertex[2]=fVertex[evtIndex][2];
1108
1109 }
1110
1111 //________________________________________
1112 void AliCaloTrackReader::FillVertexArray()
1113 {
1114
1115   //Fill data member with vertex
1116   //In case of Mixed event, multiple vertices
1117
1118   //Delete previous vertex
1119   if(fVertex)
1120   {
1121     for (Int_t i = 0; i < fNMixedEvent; i++)
1122     {
1123       delete [] fVertex[i] ;
1124     }
1125     delete [] fVertex ;
1126   }
1127
1128   fVertex = new Double_t*[fNMixedEvent] ;
1129   for (Int_t i = 0; i < fNMixedEvent; i++)
1130   {
1131     fVertex[i] = new Double_t[3] ;
1132     fVertex[i][0] = 0.0 ;
1133     fVertex[i][1] = 0.0 ;
1134     fVertex[i][2] = 0.0 ;
1135   }
1136
1137   if (!fMixedEvent)
1138   { //Single event analysis
1139     if(fDataType!=kMC)
1140     {
1141
1142       if(fInputEvent->GetPrimaryVertex())
1143       {
1144         fInputEvent->GetPrimaryVertex()->GetXYZ(fVertex[0]);
1145       }
1146       else
1147       {
1148         printf("AliCaloTrackReader::FillVertexArray() - NULL primary vertex\n");
1149         fVertex[0][0]=0.;   fVertex[0][1]=0.;   fVertex[0][2]=0.;
1150       }//Primary vertex pointer do not exist
1151
1152     } else
1153     {//MC read event
1154       fVertex[0][0]=0.;   fVertex[0][1]=0.;   fVertex[0][2]=0.;
1155     }
1156
1157     if(fDebug > 1)
1158       printf("AliCaloTrackReader::FillVertexArray() - Single Event Vertex : %f,%f,%f\n",fVertex[0][0],fVertex[0][1],fVertex[0][2]);
1159
1160   } else
1161   { // MultiEvent analysis
1162     for (Int_t iev = 0; iev < fNMixedEvent; iev++)
1163     {
1164       if (fMixedEvent->GetVertexOfEvent(iev))
1165         fMixedEvent->GetVertexOfEvent(iev)->GetXYZ(fVertex[iev]);
1166       else
1167       { // no vertex found !!!!
1168         AliWarning("No vertex found");
1169       }
1170
1171       if(fDebug > 1)
1172         printf("AliCaloTrackReader::FillVertexArray() - Multi Event %d Vertex : %f,%f,%f\n",iev,fVertex[iev][0],fVertex[iev][1],fVertex[iev][2]);
1173
1174     }
1175   }
1176
1177 }
1178
1179 //_____________________________________
1180 void AliCaloTrackReader::FillInputCTS()
1181 {
1182   //Return array with Central Tracking System (CTS) tracks
1183
1184   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputCTS()\n");
1185
1186   Double_t pTrack[3] = {0,0,0};
1187
1188   Int_t nTracks = fInputEvent->GetNumberOfTracks() ;
1189   fTrackMult    = 0;
1190   Int_t nstatus = 0;
1191   Double_t bz   = GetInputEvent()->GetMagneticField();
1192
1193   for(Int_t i = 0; i < 19; i++)
1194   {
1195     fTrackBCEvent   [i] = 0;
1196     fTrackBCEventCut[i] = 0;
1197   }
1198
1199   Bool_t   bc0  = kFALSE;
1200   if(fRecalculateVertexBC) fVertexBC=AliVTrack::kTOFBCNA;
1201
1202   for (Int_t itrack =  0; itrack <  nTracks; itrack++)
1203   {////////////// track loop
1204     AliVTrack * track = (AliVTrack*)fInputEvent->GetTrack(itrack) ; // retrieve track from esd
1205
1206     //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
1207     ULong_t status = track->GetStatus();
1208
1209     if (fTrackStatus && !((status & fTrackStatus) == fTrackStatus))
1210       continue ;
1211
1212     nstatus++;
1213
1214     Float_t dcaTPC =-999;
1215
1216     if     (fDataType==kESD)
1217     {
1218       AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (track);
1219
1220       if(esdTrack)
1221       {
1222         if(fESDtrackCuts->AcceptTrack(esdTrack))
1223         {
1224           track->GetPxPyPz(pTrack) ;
1225
1226           if(fConstrainTrack)
1227           {
1228             if(esdTrack->GetConstrainedParam())
1229             {
1230               const AliExternalTrackParam* constrainParam = esdTrack->GetConstrainedParam();
1231               esdTrack->Set(constrainParam->GetX(),constrainParam->GetAlpha(),constrainParam->GetParameter(),constrainParam->GetCovariance());
1232               esdTrack->GetConstrainedPxPyPz(pTrack);
1233             }
1234             else continue;
1235
1236           } // use constrained tracks
1237
1238           if(fSelectSPDHitTracks)
1239           {//Not much sense to use with TPC only or Hybrid tracks
1240             if(!esdTrack->HasPointOnITSLayer(0) && !esdTrack->HasPointOnITSLayer(1)) continue ;
1241           }
1242         }
1243         // Complementary track to global : Hybrids (make sure that the previous selection is for Global)
1244         else  if(fESDtrackComplementaryCuts && fESDtrackComplementaryCuts->AcceptTrack(esdTrack))
1245         {
1246           // constrain the track
1247           if(esdTrack->GetConstrainedParam())
1248           {
1249             esdTrack->Set(esdTrack->GetConstrainedParam()->GetX(),esdTrack->GetConstrainedParam()->GetAlpha(),esdTrack->GetConstrainedParam()->GetParameter(),esdTrack->GetConstrainedParam()->GetCovariance());
1250
1251             track->GetPxPyPz(pTrack) ;
1252
1253           }
1254           else continue;
1255         }
1256         else continue;
1257       }
1258     } // ESD
1259     else if(fDataType==kAOD)
1260     {
1261       AliAODTrack *aodtrack = dynamic_cast <AliAODTrack*>(track);
1262
1263       if(aodtrack)
1264       {
1265        if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputCTS():AOD track type: %d (primary %d), hybrid? %d \n",
1266                               aodtrack->GetType(),AliAODTrack::kPrimary,
1267                               aodtrack->IsHybridGlobalConstrainedGlobal());
1268
1269
1270         if (fSelectHybridTracks)
1271         {
1272           if (!aodtrack->IsHybridGlobalConstrainedGlobal())       continue ;
1273         }
1274         else
1275         {
1276           if ( aodtrack->TestFilterBit(fTrackFilterMask)==kFALSE) continue ;
1277         }
1278
1279         if(fSelectSPDHitTracks)
1280         {//Not much sense to use with TPC only or Hybrid tracks
1281           if(!aodtrack->HasPointOnITSLayer(0) && !aodtrack->HasPointOnITSLayer(1)) continue ;
1282         }
1283
1284         if (aodtrack->GetType()!= AliAODTrack::kPrimary)          continue ;
1285
1286         if (fDebug > 2 ) printf("AliCaloTrackReader::FillInputCTS(): \t accepted track! \n");
1287
1288         //In case of AODs, TPC tracks cannot be propagated back to primary vertex,
1289         // info stored here
1290         dcaTPC = aodtrack->DCA();
1291
1292         track->GetPxPyPz(pTrack) ;
1293
1294       } // aod track exists
1295       else continue ;
1296
1297     } // AOD
1298
1299     TLorentzVector momentum(pTrack[0],pTrack[1],pTrack[2],0);
1300
1301     Bool_t okTOF  = ( (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ) ;
1302     Double_t tof  = -1000;
1303     Int_t trackBC = -1000 ;
1304
1305     if(okTOF)
1306     {
1307       trackBC = track->GetTOFBunchCrossing(bz);
1308       SetTrackEventBC(trackBC+9);
1309
1310       tof = track->GetTOFsignal()*1e-3;
1311     }
1312
1313     if(fUseTrackDCACut)
1314     {
1315       //normal way to get the dca, cut on dca_xy
1316       if(dcaTPC==-999)
1317       {
1318         Double_t dca[2]   = {1e6,1e6};
1319         Double_t covar[3] = {1e6,1e6,1e6};
1320         Bool_t okDCA = track->PropagateToDCA(fInputEvent->GetPrimaryVertex(),bz,100.,dca,covar);
1321         if( okDCA) okDCA = AcceptDCA(momentum.Pt(),dca[0]);
1322         if(!okDCA)
1323         {
1324           //printf("AliCaloTrackReader::FillInputCTS() - Reject track pt %2.2f, dca_xy %2.4f, BC %d\n",momentum.Pt(),dca[0],trackBC);
1325           continue ;
1326         }
1327       }
1328     }// DCA cuts
1329
1330     if(okTOF)
1331     {
1332       //SetTrackEventBCcut(bc);
1333       SetTrackEventBCcut(trackBC+9);
1334
1335       //After selecting tracks with small DCA, pointing to vertex, set vertex BC depeding on tracks BC
1336       if(fRecalculateVertexBC)
1337       {
1338         if     (trackBC !=0 && trackBC != AliVTrack::kTOFBCNA) fVertexBC = trackBC;
1339         else if(trackBC == 0)                                  bc0 = kTRUE;
1340       }
1341
1342       //In any case, the time should to be larger than the fixed window ...
1343       if( fUseTrackTimeCut && (trackBC!=0 || tof < fTrackTimeCutMin  || tof > fTrackTimeCutMax) )
1344       {
1345         //printf("Remove track time %f and bc = %d\n",tof,trackBC);
1346         continue ;
1347       }
1348       //else printf("Accept track time %f and bc = %d\n",tof,trackBC);
1349
1350     }
1351
1352     //Count the tracks in eta < 0.9
1353     //printf("Eta %f cut  %f\n",TMath::Abs(track->Eta()),fTrackMultEtaCut);
1354     if(TMath::Abs(track->Eta())< fTrackMultEtaCut) fTrackMult++;
1355
1356     if(fCTSPtMin > momentum.Pt() || fCTSPtMax < momentum.Pt()) continue ;
1357
1358     if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"CTS")) continue;
1359
1360     if(fDebug > 2 && momentum.Pt() > 0.1)
1361       printf("AliCaloTrackReader::FillInputCTS() - Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
1362              momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
1363
1364     if (fMixedEvent)  track->SetID(itrack);
1365
1366     fCTSTracks->Add(track);
1367
1368   }// track loop
1369
1370   if(fVertexBC ==0 || fVertexBC == AliVTrack::kTOFBCNA)
1371   {
1372     if( bc0 ) fVertexBC = 0 ;
1373     else      fVertexBC = AliVTrack::kTOFBCNA ;
1374   }
1375
1376   if(fDebug > 1)
1377     printf("AliCaloTrackReader::FillInputCTS()   - aod entries %d, input tracks %d, pass status %d, multipliticy %d\n", fCTSTracks->GetEntriesFast(), nTracks, nstatus, fTrackMult);//fCTSTracksNormalInputEntries);
1378
1379 }
1380
1381 //__________________________________________________________________
1382 void AliCaloTrackReader::FillInputEMCALAlgorithm(AliVCluster * clus,
1383                                                  const Int_t iclus)
1384 {
1385   //Fill the EMCAL data in the array, do it
1386
1387   Int_t vindex = 0 ;
1388   if (fMixedEvent)
1389     vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
1390
1391   //Reject clusters with bad channels, close to borders and exotic;
1392   if(!GetCaloUtils()->GetEMCALRecoUtils()->IsGoodCluster(clus,GetCaloUtils()->GetEMCALGeometry(),GetEMCALCells(),fInputEvent->GetBunchCrossNumber())) return;
1393
1394   //Mask all cells in collumns facing ALICE thick material if requested
1395   if(GetCaloUtils()->GetNMaskCellColumns())
1396   {
1397     Int_t absId   = -1;
1398     Int_t iSupMod = -1;
1399     Int_t iphi    = -1;
1400     Int_t ieta    = -1;
1401     Bool_t shared = kFALSE;
1402     GetCaloUtils()->GetEMCALRecoUtils()->GetMaxEnergyCell(GetCaloUtils()->GetEMCALGeometry(), GetEMCALCells(),clus,absId,iSupMod,ieta,iphi,shared);
1403     if(GetCaloUtils()->MaskFrameCluster(iSupMod, ieta)) return;
1404   }
1405
1406   if(fSelectEmbeddedClusters)
1407   {
1408     if(clus->GetNLabels()==0 || clus->GetLabel() < 0) return;
1409     //else printf("Embedded cluster,  %d, n label %d label %d  \n",iclus,clus->GetNLabels(),clus->GetLabel());
1410   }
1411
1412   //Float_t pos[3];
1413   //clus->GetPosition(pos);
1414   //printf("Before Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
1415
1416   if(fRecalculateClusters)
1417   {
1418     //Recalibrate the cluster energy
1419     if(GetCaloUtils()->IsRecalibrationOn())
1420     {
1421       Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, GetEMCALCells());
1422
1423       clus->SetE(energy);
1424       //printf("Recalibrated Energy %f\n",clus->E());
1425
1426       GetCaloUtils()->RecalculateClusterShowerShapeParameters(GetEMCALCells(),clus);
1427       GetCaloUtils()->RecalculateClusterPID(clus);
1428
1429     } // recalculate E
1430
1431     //Recalculate distance to bad channels, if new list of bad channels provided
1432     GetCaloUtils()->RecalculateClusterDistanceToBadChannel(GetEMCALCells(),clus);
1433
1434     //Recalculate cluster position
1435     if(GetCaloUtils()->IsRecalculationOfClusterPositionOn())
1436     {
1437       GetCaloUtils()->RecalculateClusterPosition(GetEMCALCells(),clus);
1438       //clus->GetPosition(pos);
1439       //printf("After  Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
1440     }
1441
1442     // Recalculate TOF
1443     if(GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn())
1444     {
1445       Double_t tof      = clus->GetTOF();
1446       Float_t  frac     =-1;
1447       Int_t    absIdMax = GetCaloUtils()->GetMaxEnergyCell(fEMCALCells, clus,frac);
1448
1449       if(fDataType==AliCaloTrackReader::kESD)
1450       {
1451         tof = fEMCALCells->GetCellTime(absIdMax);
1452       }
1453
1454       GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
1455
1456       clus->SetTOF(tof);
1457
1458     }// Time recalibration
1459   }
1460
1461   //Correct non linearity
1462   if(GetCaloUtils()->IsCorrectionOfClusterEnergyOn())
1463   {
1464     GetCaloUtils()->CorrectClusterEnergy(clus) ;
1465     //printf("Linearity Corrected Energy %f\n",clus->E());
1466
1467     //In case of MC analysis, to match resolution/calibration in real data
1468     Float_t rdmEnergy = GetCaloUtils()->GetEMCALRecoUtils()->SmearClusterEnergy(clus);
1469     // printf("\t Energy %f, smeared %f\n", clus->E(),rdmEnergy);
1470     clus->SetE(rdmEnergy);
1471   }
1472
1473   Double_t tof = clus->GetTOF()*1e9;
1474
1475   Int_t bc = TMath::Nint(tof/50) + 9;
1476   //printf("tof %2.2f, bc+5=%d\n",tof,bc);
1477
1478   SetEMCalEventBC(bc);
1479
1480   if(fEMCALPtMin > clus->E() || fEMCALPtMax < clus->E()) return ;
1481
1482   TLorentzVector momentum ;
1483
1484   clus->GetMomentum(momentum, fVertex[vindex]);
1485
1486   if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"EMCAL")) return ;
1487
1488   SetEMCalEventBCcut(bc);
1489
1490   if(!IsInTimeWindow(tof,clus->E()))
1491   {
1492     fNPileUpClusters++ ;
1493     if(fUseEMCALTimeCut) return ;
1494   }
1495   else
1496     fNNonPileUpClusters++;
1497
1498   if(fDebug > 2 && momentum.E() > 0.1)
1499     printf("AliCaloTrackReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
1500            momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
1501
1502   if (fMixedEvent)
1503     clus->SetID(iclus) ;
1504
1505   //Correct MC label for AODs
1506
1507   fEMCALClusters->Add(clus);
1508
1509 }
1510
1511 //_______________________________________
1512 void AliCaloTrackReader::FillInputEMCAL()
1513 {
1514   //Return array with EMCAL clusters in aod format
1515
1516   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputEMCAL()\n");
1517
1518   // First recalibrate cells, time or energy
1519   //  if(GetCaloUtils()->IsRecalibrationOn())
1520   //    GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCells(GetCaloUtils()->GetEMCALGeometry(),
1521   //                                                          GetEMCALCells(),
1522   //                                                          fInputEvent->GetBunchCrossNumber());
1523
1524   fNPileUpClusters    = 0; // Init counter
1525   fNNonPileUpClusters = 0; // Init counter
1526   for(Int_t i = 0; i < 19; i++)
1527   {
1528     fEMCalBCEvent   [i] = 0;
1529     fEMCalBCEventCut[i] = 0;
1530   }
1531
1532   //Loop to select clusters in fiducial cut and fill container with aodClusters
1533   if(fEMCALClustersListName=="")
1534   {
1535     Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
1536     for (Int_t iclus =  0; iclus <  nclusters; iclus++)
1537     {
1538       AliVCluster * clus = 0;
1539       if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
1540       {
1541         if (IsEMCALCluster(clus))
1542         {
1543           FillInputEMCALAlgorithm(clus, iclus);
1544         }//EMCAL cluster
1545       }// cluster exists
1546     }// cluster loop
1547
1548     //Recalculate track matching
1549     GetCaloUtils()->RecalculateClusterTrackMatching(fInputEvent);
1550
1551   }//Get the clusters from the input event
1552   else
1553   {
1554     TClonesArray * clusterList = 0x0;
1555
1556     if      (fInputEvent->FindListObject(fEMCALClustersListName))
1557     {
1558       clusterList = dynamic_cast<TClonesArray*> (fInputEvent->FindListObject(fEMCALClustersListName));
1559     }
1560     else if(fOutputEvent)
1561     {
1562       clusterList = dynamic_cast<TClonesArray*> (fOutputEvent->FindListObject(fEMCALClustersListName));
1563     }
1564
1565     if(!clusterList)
1566     {
1567         printf("AliCaloTrackReader::FillInputEMCAL() - Wrong name of list with clusters?  <%s>\n",fEMCALClustersListName.Data());
1568         return;
1569     }
1570
1571     Int_t nclusters = clusterList->GetEntriesFast();
1572     for (Int_t iclus =  0; iclus <  nclusters; iclus++)
1573     {
1574       AliVCluster * clus = dynamic_cast<AliVCluster*> (clusterList->At(iclus));
1575       //printf("E %f\n",clus->E());
1576       if (clus) FillInputEMCALAlgorithm(clus, iclus);
1577       else printf("AliCaloTrackReader::FillInputEMCAL() - Null cluster in list!\n");
1578     }// cluster loop
1579
1580     // Recalculate the pile-up time, in case long time clusters removed during clusterization
1581     //printf("Input event INIT : Pile-up clusters %d, NO pile-up %d\n",fNPileUpClusters,fNNonPileUpClusters);
1582
1583     fNPileUpClusters    = 0; // Init counter
1584     fNNonPileUpClusters = 0; // Init counter
1585     for(Int_t i = 0; i < 19; i++)
1586     {
1587       fEMCalBCEvent   [i] = 0;
1588       fEMCalBCEventCut[i] = 0;
1589     }
1590
1591     for (Int_t iclus =  0; iclus < fInputEvent->GetNumberOfCaloClusters(); iclus++)
1592     {
1593       AliVCluster * clus = 0;
1594
1595       if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
1596       {
1597         if (IsEMCALCluster(clus))
1598         {
1599
1600           Float_t  frac     =-1;
1601           Int_t    absIdMax = GetCaloUtils()->GetMaxEnergyCell(fEMCALCells, clus,frac);
1602           Double_t tof = clus->GetTOF();
1603           GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
1604           tof*=1e9;
1605
1606           //printf("Input event cluster : AbsIdMax %d, E %2.2f, time %2.2f \n", absIdMax,clus->E(),tof);
1607
1608           //Reject clusters with bad channels, close to borders and exotic;
1609           if(!GetCaloUtils()->GetEMCALRecoUtils()->IsGoodCluster(clus,GetCaloUtils()->GetEMCALGeometry(),GetEMCALCells(),fInputEvent->GetBunchCrossNumber()))  continue;
1610
1611           Int_t bc = TMath::Nint(tof/50) + 9;
1612           SetEMCalEventBC(bc);
1613
1614           if(fEMCALPtMin > clus->E() || fEMCALPtMax < clus->E()) continue ;
1615
1616           TLorentzVector momentum ;
1617
1618           clus->GetMomentum(momentum, fVertex[0]);
1619
1620           if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"EMCAL")) return ;
1621
1622           SetEMCalEventBCcut(bc);
1623
1624           if(!IsInTimeWindow(tof,clus->E()))
1625             fNPileUpClusters++ ;
1626           else
1627             fNNonPileUpClusters++;
1628
1629         }
1630       }
1631     }
1632
1633     //printf("Input event : Pile-up clusters %d, NO pile-up %d\n",fNPileUpClusters,fNNonPileUpClusters);
1634
1635     // Recalculate track matching, not necessary if already done in the reclusterization task.
1636     // in case it was not done ...
1637     GetCaloUtils()->RecalculateClusterTrackMatching(fInputEvent,clusterList);
1638
1639   }
1640
1641   if(fDebug > 1) printf("AliCaloTrackReader::FillInputEMCAL() - aod entries %d, n pile-up clusters %d, n non pile-up %d \n",  fEMCALClusters->GetEntriesFast(),fNPileUpClusters,fNNonPileUpClusters);
1642
1643 }
1644
1645 //______________________________________
1646 void AliCaloTrackReader::FillInputPHOS()
1647 {
1648   //Return array with PHOS clusters in aod format
1649
1650   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputPHOS()\n");
1651
1652   //Loop to select clusters in fiducial cut and fill container with aodClusters
1653   Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
1654   for (Int_t iclus = 0; iclus < nclusters; iclus++)
1655   {
1656     AliVCluster * clus = 0;
1657     if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
1658     {
1659       if (IsPHOSCluster(clus))
1660       {
1661         //Check if the cluster contains any bad channel and if close to calorimeter borders
1662         Int_t vindex = 0 ;
1663         if (fMixedEvent)
1664           vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
1665         if( GetCaloUtils()->ClusterContainsBadChannel("PHOS",clus->GetCellsAbsId(), clus->GetNCells()))
1666           continue;
1667         if(!GetCaloUtils()->CheckCellFiducialRegion(clus, fInputEvent->GetPHOSCells(), fInputEvent, vindex))
1668           continue;
1669
1670         if(fRecalculateClusters)
1671         {
1672           //Recalibrate the cluster energy
1673           if(GetCaloUtils()->IsRecalibrationOn())
1674           {
1675             Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliAODCaloCells*)GetPHOSCells());
1676             clus->SetE(energy);
1677           }
1678         }
1679
1680         TLorentzVector momentum ;
1681
1682         clus->GetMomentum(momentum, fVertex[vindex]);
1683
1684         if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"PHOS")) continue;
1685
1686         if(fPHOSPtMin > momentum.E() || fPHOSPtMax < momentum.E())          continue;
1687
1688         if(fDebug > 2 && momentum.E() > 0.1)
1689           printf("AliCaloTrackReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
1690                  momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
1691
1692
1693         if (fMixedEvent)
1694         {
1695           clus->SetID(iclus) ;
1696         }
1697
1698         fPHOSClusters->Add(clus);
1699
1700       }//PHOS cluster
1701     }//cluster exists
1702   }//esd cluster loop
1703
1704   if(fDebug > 1) printf("AliCaloTrackReader::FillInputPHOS()  - aod entries %d\n",  fPHOSClusters->GetEntriesFast());
1705
1706 }
1707
1708 //____________________________________________
1709 void AliCaloTrackReader::FillInputEMCALCells()
1710 {
1711   //Return array with EMCAL cells in aod format
1712
1713   fEMCALCells = fInputEvent->GetEMCALCells();
1714
1715 }
1716
1717 //___________________________________________
1718 void AliCaloTrackReader::FillInputPHOSCells()
1719 {
1720   //Return array with PHOS cells in aod format
1721
1722   fPHOSCells = fInputEvent->GetPHOSCells();
1723
1724 }
1725
1726 //_______________________________________
1727 void AliCaloTrackReader::FillInputVZERO()
1728 {
1729   //Fill VZERO information in data member, add all the channels information.
1730   AliVVZERO* v0 = fInputEvent->GetVZEROData();
1731   //printf("Init V0: ADC (%d,%d), Multiplicity (%d,%d) \n",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]);
1732
1733   if (v0)
1734   {
1735     AliESDVZERO* esdV0 = dynamic_cast<AliESDVZERO*> (v0);
1736     for (Int_t i = 0; i < 32; i++)
1737     {
1738       if(esdV0)
1739       {//Only available in ESDs
1740         fV0ADC[0] += (Int_t)esdV0->GetAdcV0C(i);
1741         fV0ADC[1] += (Int_t)esdV0->GetAdcV0A(i);
1742       }
1743
1744       fV0Mul[0] += (Int_t)v0->GetMultiplicityV0C(i);
1745       fV0Mul[1] += (Int_t)v0->GetMultiplicityV0A(i);
1746     }
1747     if(fDebug > 0)
1748       printf("V0: ADC (%d,%d), Multiplicity (%d,%d) \n",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]);
1749   }
1750   else
1751   {
1752     if(fDebug > 0)
1753       printf("Cannot retrieve V0 ESD! Run w/ null V0 charges\n ");
1754   }
1755 }
1756
1757
1758 //___________________________________________________________________
1759 Bool_t AliCaloTrackReader::IsEMCALCluster(AliVCluster* cluster) const
1760 {
1761   // Check if it is a cluster from EMCAL. For old AODs cluster type has
1762   // different number and need to patch here
1763
1764   if(fDataType==kAOD && fOldAOD)
1765   {
1766     if (cluster->GetType() == 2) return kTRUE;
1767     else                         return kFALSE;
1768   }
1769   else
1770   {
1771     return cluster->IsEMCAL();
1772   }
1773
1774 }
1775
1776 //___________________________________________________________________
1777 Bool_t AliCaloTrackReader::IsPHOSCluster(AliVCluster * cluster) const
1778 {
1779   //Check if it is a cluster from PHOS.For old AODs cluster type has
1780   // different number and need to patch here
1781
1782   if(fDataType==kAOD && fOldAOD)
1783   {
1784     Int_t type = cluster->GetType();
1785     if (type == 0 || type == 1) return kTRUE;
1786     else                        return kFALSE;
1787   }
1788   else
1789   {
1790     return cluster->IsPHOS();
1791   }
1792
1793 }
1794
1795 //________________________________________________
1796 Bool_t AliCaloTrackReader::CheckForPrimaryVertex()
1797 {
1798   //Check if the vertex was well reconstructed, copy from V0Reader of conversion group
1799   //Only for ESDs ...
1800
1801   AliESDEvent * event = dynamic_cast<AliESDEvent*> (fInputEvent);
1802   if(!event) return kTRUE;
1803
1804   if(event->GetPrimaryVertexTracks()->GetNContributors() > 0)
1805   {
1806     return kTRUE;
1807   }
1808
1809   if(event->GetPrimaryVertexTracks()->GetNContributors() < 1)
1810   {
1811     // SPD vertex
1812     if(event->GetPrimaryVertexSPD()->GetNContributors() > 0)
1813     {
1814       //cout<<"spd vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
1815       return kTRUE;
1816
1817     }
1818     if(event->GetPrimaryVertexSPD()->GetNContributors() < 1)
1819     {
1820       //      cout<<"bad vertex type::"<< event->GetPrimaryVertex()->GetName() << endl;
1821       return kFALSE;
1822     }
1823   }
1824
1825   return kFALSE;
1826
1827 }
1828
1829 //_______________________________________________
1830 TArrayI AliCaloTrackReader::GetL0TriggerPatches()
1831 {
1832   //Select the patches that triggered
1833
1834   // some variables
1835   Int_t  trigtimes[30], globCol, globRow,ntimes, i;
1836   Int_t  absId  = -1; //[100];
1837   Int_t  nPatch = 0;
1838
1839   TArrayI patches(0);
1840
1841   // get object pointer
1842   AliVCaloTrigger *caloTrigger = GetInputEvent()->GetCaloTrigger( "EMCAL" );
1843
1844   // class is not empty
1845   if( caloTrigger->GetEntries() > 0 )
1846   {
1847     // must reset before usage, or the class will fail
1848     caloTrigger->Reset();
1849
1850     // go throuth the trigger channels
1851     while( caloTrigger->Next() )
1852     {
1853       // get position in global 2x2 tower coordinates
1854       caloTrigger->GetPosition( globCol, globRow );
1855
1856       // get dimension of time arrays
1857       caloTrigger->GetNL0Times( ntimes );
1858
1859       // no L0s in this channel
1860       // presence of the channel in the iterator still does not guarantee that L0 was produced!!
1861       if( ntimes < 1 )
1862         continue;
1863
1864       // get timing array
1865       caloTrigger->GetL0Times( trigtimes );
1866
1867       //printf("trigger time window %d - %d\n",fTriggerPatchTimeWindow[0],fTriggerPatchTimeWindow[1]);
1868       // go through the array
1869       for( i = 0; i < ntimes; i++ )
1870       {
1871         // check if in cut - 8,9 shall be accepted in 2011
1872         if( trigtimes[i] >= fTriggerPatchTimeWindow[0] && trigtimes[i] <= fTriggerPatchTimeWindow[1] )
1873         {
1874           //printf("Accepted trigger time %d \n",trigtimes[i]);
1875           //if(nTrig > 99) continue;
1876           GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(globCol,globRow, absId);//absIDTrig[nTrig]) ;
1877           //printf("pass the time cut globCol %d, globRow %d absId %d\n",globCol,globRow, absIDTrig[nTrig]);
1878           patches.Set(nPatch+1);
1879           patches.AddAt(absId,nPatch++);
1880         }
1881       } // trigger time array
1882
1883     } // trigger iterator
1884   } // go thorough triggers
1885
1886   //printf("N patches %d, test %d,first %d, last %d\n",patches.GetSize(), nTrig, patches.At(0), patches.At(patches.GetSize()-1));
1887
1888   return patches;
1889 }
1890
1891
1892 /*
1893 //___________________________________________________________
1894 void  AliCaloTrackReader::RejectExoticEvents(TArrayI patches)
1895 {
1896   // Reject events triggered by exotic cells
1897
1898   if(!GetFiredTriggerClasses().Contains("EMC") && !fForceExoticRejection)
1899   {
1900     fIsExoticEvent = kFALSE;
1901     return;
1902   }
1903
1904   Int_t nclusters  = fInputEvent->GetNumberOfCaloClusters();
1905
1906   // simple method, just count how many exotics and how many high energy clusters
1907   // over the trigger threshold there are
1908
1909   if(!fTriggerPatchExoticRejection)
1910   {
1911     Int_t nOfHighECl = 0 ;
1912     Int_t nExo       = 0 ;
1913     for (Int_t iclus = 0 ; iclus <  nclusters; iclus++)
1914     {
1915       AliVCluster * clus = 0;
1916       if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
1917       {
1918         if (IsEMCALCluster(clus) && clus->E() > fTriggerEventThreshold)
1919         {
1920           Bool_t exotic = GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCluster(clus, fInputEvent->GetEMCALCells());
1921           Bool_t bad    = GetCaloUtils()->GetEMCALRecoUtils()->ClusterContainsBadChannel(GetCaloUtils()->GetEMCALGeometry(),
1922                                                                                          clus->GetCellsAbsId(),clus->GetNCells());
1923
1924           if     ( exotic)
1925           {
1926             //printf("- Exotic E %f, tof %f\n",clus->E(),clus->GetTOF()*1e9);
1927             nExo++;
1928           }
1929           else if(!bad  && !exotic) nOfHighECl++;
1930
1931         }//EMCAL cluster
1932       }// cluster exists
1933     }// cluster loop
1934
1935     //printf("- trigger %2.1f, n Exotic %d, n high energy %d\n", fTriggerEventThreshold,nExo,nOfHighECl);
1936
1937     if ( ( nExo > 0 ) && ( nOfHighECl == 0 ) )
1938     {
1939       fIsExoticEvent = kTRUE ;
1940       //printf("*** Exotic Event\n");
1941     }
1942     else
1943       fIsExoticEvent = kFALSE;
1944
1945   }
1946   //Check if there is any trigger patch that has an associated exotic cluster
1947   else
1948   {
1949     //printf("Trigger Exotic?\n");
1950
1951
1952     Int_t nPatch = patches.GetSize();
1953
1954 //    for(Int_t iabsId =0; iabsId < nTrig; iabsId++)
1955 //    {
1956 //      Int_t absIDCell[4];
1957 //      GetCaloUtils()->GetEMCALGeometry()->GetCellIndexFromFastORIndex(absIDTrig[iabsId], absIDCell);
1958 //      printf("Patch %d absIdTrigger %d: AbsIdCells: %d - %d - %d - %d \n",iabsId,absIDTrig[iabsId],absIDCell[0],absIDCell[1],absIDCell[2],absIDCell[3]);
1959 //    }
1960
1961     // Loop on the clusters, check if there is any that falls into one of the patches
1962     for (Int_t iclus =  0; iclus <  nclusters; iclus++)
1963     {
1964       AliVCluster * clus = 0;
1965       fIsExoticEvent = kFALSE;
1966       if ( (clus = fInputEvent->GetCaloCluster(iclus)) && IsEMCALCluster(clus))
1967       {
1968         Bool_t exotic = GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCluster(clus, fInputEvent->GetEMCALCells());
1969
1970 //        if(nPatch > 4)
1971 //        {
1972 //          Float_t frac = -1;
1973 //          Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(fInputEvent->GetEMCALCells(), clus,frac);
1974 //          Double_t tof    = clus->GetTOF();
1975 //          Bool_t bad    = GetCaloUtils()->GetEMCALRecoUtils()->ClusterContainsBadChannel(GetCaloUtils()->GetEMCALGeometry(),
1976 //                                                                                         clus->GetCellsAbsId(),clus->GetNCells());
1977 //          GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
1978 //          printf("cluster %d, E %2.2f, tof %2.2f, AbsId max %d, exo %d, bad %d\n",iclus,clus->E(),tof*1e9,absIdMax,exotic,bad);
1979 //        }
1980
1981         if(exotic)
1982         {
1983           Float_t frac = -1;
1984           Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(fInputEvent->GetEMCALCells(), clus,frac);
1985           //Double_t tof    = clus->GetTOF();
1986           //GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
1987           //printf("cluster %d, E %2.2f, tof %2.2f, AbsId max %d\n",iclus,clus->E(),tof*1e9,absIdMax);
1988
1989           for(Int_t iabsId =0; iabsId < nPatch; iabsId++)
1990           {
1991             Int_t absIDCell[4];
1992             GetCaloUtils()->GetEMCALGeometry()->GetCellIndexFromFastORIndex(patches.At(iabsId), absIDCell);
1993             for(Int_t ipatch = 0; ipatch < 4; ipatch++)
1994             {
1995               //Loop on the cluster cells and check if any is in patch, ideally
1996               // max id should suffice, but not always
1997               //              for(Int_t iCell = 0; iCell<clus->GetNCells(); iCell++)
1998               //              {
1999               //                Int_t id = clus->GetCellsAbsId()[iCell];
2000               //                Double_t tofCell =  fInputEvent->GetEMCALCells()->GetCellTime(id);
2001               //                GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(id,fInputEvent->GetBunchCrossNumber(),tofCell);
2002               //                //printf(" id %d - E %2.2f - tof %2.2f; ",clus->GetCellsAbsId()[iCell],
2003               //                //       fInputEvent->GetEMCALCells()->GetCellAmplitude(id), tofCell*1e9);
2004               //                if(clus->GetCellsAbsId()[iCell] == absIDCell[ipatch])
2005               //                {
2006               //                  //printf("\n *** Exotic trigger : absId %d, E %2.1f \n",clus->GetCellsAbsId()[iCell],clus->E());
2007               //                  fIsExoticEvent = kTRUE;
2008               //                  //return;
2009               //                }
2010               //                //else printf("\n");
2011               //              }
2012
2013               if(absIdMax == absIDCell[ipatch])
2014               {
2015                 Double_t tof =  clus->GetTOF();
2016                 GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
2017                 //printf("*** Exotic trigger : absId %d, E %2.1f, tof %f \n",absIdMax,clus->E(), clus->GetTOF()*1e9);
2018
2019                 fIsExoticEvent = kTRUE;
2020                 return;
2021               }
2022             }// cell patch loop
2023           }// trigger patch loop
2024         }// exotic cluster
2025       }// EMCal cluster
2026     }// Cluster loop
2027   }// exotic and trigger patch event flag
2028
2029 }
2030 */
2031
2032 /*
2033 //______________________________________________________________________
2034 void  AliCaloTrackReader::RejectTriggeredEventsByPileUp(TArrayI patches)
2035 {
2036   // Reject events triggered by exotic cells
2037
2038   if(!GetFiredTriggerClasses().Contains("EMC") && !fForceExoticRejection)
2039   {
2040     fTriggerClusterBC = 0;
2041     return;
2042   }
2043
2044   TClonesArray * clusterList = 0;
2045   if      (fInputEvent->FindListObject(fEMCALClustersListName))
2046   {
2047     clusterList = dynamic_cast<TClonesArray*> (fInputEvent->FindListObject(fEMCALClustersListName));
2048   }
2049   else if(fOutputEvent)
2050   {
2051     clusterList = dynamic_cast<TClonesArray*> (fOutputEvent->FindListObject(fEMCALClustersListName));
2052   }
2053
2054   Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
2055   if(clusterList) nclusters = clusterList->GetEntriesFast();
2056
2057   // simple method, just count how many exotics and how many high energy clusters
2058   // over the trigger threshold there are
2059
2060   if(!fTriggerPatchExoticRejection)
2061   {
2062     Int_t nOfHighECl   = 0 ;
2063     Int_t nOutBC       = 0 ;
2064     Float_t tofcluster = 1000;
2065     Float_t eBCN       =-1;
2066
2067     for (Int_t iclus = 0 ; iclus <  nclusters; iclus++)
2068     {
2069       AliVCluster * clus = 0;
2070
2071       if(clusterList) clus = (AliVCluster*) clusterList->At(iclus);
2072       else            clus = fInputEvent->GetCaloCluster(iclus);
2073
2074       if ( clus )
2075       {
2076         if (IsEMCALCluster(clus) && clus->E() > fTriggerEventThreshold)
2077         {
2078           Bool_t bad    = GetCaloUtils()->GetEMCALRecoUtils()->ClusterContainsBadChannel(GetCaloUtils()->GetEMCALGeometry(),
2079                                                                                          clus->GetCellsAbsId(),clus->GetNCells());
2080
2081           Bool_t exotic = GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCluster(clus, fInputEvent->GetEMCALCells());
2082
2083           if(bad || exotic || clus->E() < 1) continue;
2084
2085           Float_t frac = -1;
2086           Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(fInputEvent->GetEMCALCells(), clus,frac);
2087
2088           Double_t tof   = clus->GetTOF();
2089           if(GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn())
2090             GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
2091           tof *=1.e9;
2092
2093           if ( TMath::Abs(tof) > 25 )
2094           {
2095             nOutBC++;
2096             if(clus->E() > eBCN)
2097             {
2098               tofcluster = tof;
2099               eBCN       = clus->E();
2100             }
2101
2102             //printf("Simple: Large time:  E %f, tof %f, absId %d\n",clus->E(),tofcluster, absIdMax);
2103
2104           }
2105           else
2106           {
2107             //printf("Simple: OK cluster E %f, tof %f\n",clus->E(),tofcluster);
2108             nOfHighECl++;
2109           }
2110
2111         }//EMCAL cluster
2112       }// cluster exists
2113     }// cluster loop
2114
2115     //printf("- n OutBC %d, n high energy %d\n", nOutBC,nOfHighECl);
2116
2117     Double_t tofclusterUS = TMath::Abs(tofcluster);
2118     if ( ( nOutBC > 0 ) && ( nOfHighECl == 0 ) )
2119     {
2120       if     (tofclusterUS < 75 ) fTriggerClusterBC = 1 ;
2121       else if(tofclusterUS < 125) fTriggerClusterBC = 2 ;
2122       else if(tofclusterUS < 175) fTriggerClusterBC = 3 ;
2123       else if(tofclusterUS < 225) fTriggerClusterBC = 4 ;
2124       else if(tofclusterUS < 275) fTriggerClusterBC = 5 ;
2125       else                        fTriggerClusterBC = 6 ;
2126
2127       if(tofcluster < 0) fTriggerClusterBC*=-1;
2128
2129     }
2130     else if(( nOutBC == 0 ) && ( nOfHighECl == 0 ) )
2131     {
2132       fTriggerClusterBC = 7 ;
2133     }
2134     else
2135     {
2136       fTriggerClusterBC = 0;
2137     }
2138
2139     //printf("*** Simple: Trigger tag BC %d, tof %2.2f, E %2.2f\n",fTriggerClusterBC, tofcluster, eBCN);
2140
2141   }
2142   //Check if there is any trigger patch that has an associated exotic cluster
2143   else
2144   {
2145     //printf("Trigger Exotic?\n");
2146
2147     Int_t nPatch = patches.GetSize();
2148
2149     // Loop on the clusters, check if there is any that falls into one of the patches
2150
2151     Float_t tofcluster = 1000;
2152     Float_t eBCN       =-1;
2153
2154     Bool_t ok   = kFALSE;
2155     Bool_t ok2  = kFALSE;
2156
2157     for (Int_t iclus =  0; iclus <  nclusters; iclus++)
2158     {
2159       AliVCluster * clus = 0;
2160       if(clusterList) clus = (AliVCluster*) clusterList->At(iclus);
2161       else            clus = fInputEvent->GetCaloCluster(iclus);
2162
2163       if ( clus )
2164       {
2165         Bool_t bad    = GetCaloUtils()->GetEMCALRecoUtils()->ClusterContainsBadChannel(GetCaloUtils()->GetEMCALGeometry(),
2166                                                                                        clus->GetCellsAbsId(),clus->GetNCells());
2167
2168         Bool_t exotic = GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCluster(clus, fInputEvent->GetEMCALCells());
2169
2170         if(bad || exotic || clus->E() < 1) continue;
2171
2172         Float_t frac   = -1;
2173         Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(fInputEvent->GetEMCALCells(), clus,frac);
2174
2175         Double_t tof   = clus->GetTOF();
2176
2177         if(GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn())
2178           GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
2179         tof *=1.e9;
2180
2181         // in case no final match with the trigger, check if there was a high energy cluster
2182         // with the timing of BC=0
2183         if(clus->E() > fTriggerEventThreshold)
2184         {
2185           if(TMath::Abs(tof) < 25) ok  = kTRUE;
2186           else                     ok2 = kTRUE;
2187         }
2188         //printf("cluster %d, E %2.2f, tof %2.2f, AbsId max %d, exotic %d, bad %d\n",iclus,clus->E(),tof,absIdMax, exotic, bad);
2189
2190  //       if ( TMath::Abs(tof) > 25 )
2191  //       {
2192           //printf("cluster %d, E %2.2f, tof %2.2f, AbsId max %d\n",iclus,clus->E(),tof*1e9,absIdMax);
2193
2194           for(Int_t iabsId =0; iabsId < nPatch; iabsId++)
2195           {
2196             Int_t absIDCell[4];
2197             GetCaloUtils()->GetEMCALGeometry()->GetCellIndexFromFastORIndex(patches.At(iabsId), absIDCell);
2198             //if(tof > 75 ) printf("E %2.2f TOF %2.2f Trigger patch %d, cells : %d, %d, %d, %d\n",
2199             //                     clus->E(),tof,patches.At(iabsId), absIDCell[0],absIDCell[1],absIDCell[2],absIDCell[3]);
2200
2201
2202             for(Int_t ipatch = 0; ipatch < 4; ipatch++)
2203             {
2204               if(absIdMax == absIDCell[ipatch])
2205               {
2206                 //printf("*** Patches : absId %d, E %2.1f, tof %f \n",absIdMax,clus->E(), tof);
2207                 if(clus->E() > eBCN)
2208                 {
2209                   tofcluster = tof;
2210                   eBCN       = clus->E();
2211                 }
2212               }
2213             }// cell patch loop
2214           }// trigger patch loop
2215 //        }// out BC
2216       }// EMCal cluster
2217     }// Cluster loop
2218
2219     Double_t tofclusterUS = TMath::Abs(tofcluster);
2220
2221     if     (tofcluster == 1000)
2222     {
2223       if(ok) fTriggerClusterBC = 6 ; // no trigger match but high energy cluster with time at BC=0
2224       else   fTriggerClusterBC = 7 ; // no trigger match and no likely good cluster
2225     }
2226     else if(tofclusterUS < 25 ) fTriggerClusterBC = 0 ;
2227     else if(tofclusterUS < 75 ) fTriggerClusterBC = 1 ;
2228     else if(tofclusterUS < 125) fTriggerClusterBC = 2 ;
2229     else if(tofclusterUS < 175) fTriggerClusterBC = 3 ;
2230     else if(tofclusterUS < 225) fTriggerClusterBC = 4 ;
2231     else                        fTriggerClusterBC = 5 ;
2232
2233     //printf(" selected tof %f\n",tofcluster);
2234
2235     if(tofcluster < 0) fTriggerClusterBC*=-1;
2236
2237     //if(fTriggerClusterBC != 0) printf("*** Patches: Trigger out of BC %d, tof %2.2f, E %2.2f\n",fTriggerClusterBC,tofcluster,eBCN);
2238
2239     //if(fTriggerClusterBC==7) printf("*** No trigger match, high energy cluster? %d\n",ok2);
2240     //if(fTriggerClusterBC==6) printf("*** No trigger match, but high energy cluster in BC0\n");
2241
2242   }// exotic and trigger patch event flag
2243
2244 }
2245 */
2246
2247 //______________________________________________________________________
2248 void  AliCaloTrackReader::MatchTriggerCluster(TArrayI patches)
2249 {
2250   // Finds the cluster that triggered
2251
2252   // Init info from previous event
2253   fTriggerClusterIndex = -1;
2254   fTriggerClusterId    = -1;
2255   fIsTriggerMatch      = kFALSE;
2256   fTriggerClusterBC    = -10000;
2257   fIsExoticEvent       = kFALSE;
2258   fIsBadCellEvent      = kFALSE;
2259
2260   // Do only analysis for triggered events
2261   if(!GetFiredTriggerClasses().Contains("EMC"))
2262   {
2263     fTriggerClusterBC    = 0;
2264     return;
2265   }
2266
2267   //Recover the list of clusters
2268   TClonesArray * clusterList = 0;
2269   if      (fInputEvent->FindListObject(fEMCALClustersListName))
2270   {
2271     clusterList = dynamic_cast<TClonesArray*> (fInputEvent->FindListObject(fEMCALClustersListName));
2272   }
2273   else if(fOutputEvent)
2274   {
2275     clusterList = dynamic_cast<TClonesArray*> (fOutputEvent->FindListObject(fEMCALClustersListName));
2276   }
2277
2278   // Get number of clusters and of trigger patches
2279   Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
2280   if(clusterList) nclusters = clusterList->GetEntriesFast();
2281
2282   Int_t nPatch = patches.GetSize();
2283
2284
2285   //Init some variables used in the cluster loop
2286   Float_t tofPatchMax = 100000;
2287   Float_t ePatchMax   =-1;
2288
2289   Float_t tofMax      = 100000;
2290   Float_t eMax        =-1;
2291
2292   Int_t   clusMax     =-1;
2293   Int_t   idclusMax   =-1;
2294   Bool_t  badMax      = kFALSE;
2295   Bool_t  exoMax      = kFALSE;
2296
2297   Int_t   nOfHighECl  = 0 ;
2298
2299   // Loop on the clusters, check if there is any that falls into one of the patches
2300   for (Int_t iclus =  0; iclus <  nclusters; iclus++)
2301   {
2302     AliVCluster * clus = 0;
2303     if(clusterList) clus = (AliVCluster*) clusterList->At(iclus);
2304     else            clus = fInputEvent->GetCaloCluster(iclus);
2305
2306     if ( !clus )                continue ;
2307
2308     if ( !IsEMCALCluster(clus)) continue ;
2309
2310     if ( clus->E() < 1 )        continue ;
2311
2312     Bool_t bad    = GetCaloUtils()->GetEMCALRecoUtils()->ClusterContainsBadChannel(GetCaloUtils()->GetEMCALGeometry(),
2313                                                                                    clus->GetCellsAbsId(),clus->GetNCells());
2314
2315     Bool_t exotic = GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCluster(clus, fInputEvent->GetEMCALCells());
2316
2317     Float_t  energy = clus->E();
2318     Int_t    idclus = clus->GetID();
2319
2320     Float_t frac   = -1;
2321     Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(fInputEvent->GetEMCALCells(), clus,frac);
2322
2323     Double_t tof   = clus->GetTOF();
2324     if(GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn())
2325       GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
2326     tof *=1.e9;
2327
2328     //printf("cluster %d, ID %d, E %2.2f, tof %2.2f, AbsId max %d, exotic %d, bad %d\n",iclus,idclus, energy,tof,absIdMax, exotic, bad);
2329
2330     // Find the highest energy cluster, avobe trigger threshold
2331     // in the event in case no match to trigger is found
2332     if( energy > eMax )// && energy > fTriggerEventThreshold )
2333     {
2334       tofMax    = tof;
2335       eMax      = energy;
2336       badMax    = bad;
2337       exoMax    = exotic;
2338       clusMax   = iclus;
2339       idclusMax = idclus;
2340     }
2341
2342     // count the good clusters in the event avobe the trigger threshold
2343     // to check the exotic events
2344     if(!bad && !exotic)// && energy > fTriggerEventThreshold)
2345       nOfHighECl++;
2346
2347     // Find match to trigger
2348     if(fTriggerPatchClusterMatch)
2349     {
2350       for(Int_t iabsId =0; iabsId < nPatch; iabsId++)
2351       {
2352         Int_t absIDCell[4];
2353         GetCaloUtils()->GetEMCALGeometry()->GetCellIndexFromFastORIndex(patches.At(iabsId), absIDCell);
2354         //if(tof > 75 ) printf("E %2.2f TOF %2.2f Trigger patch %d, cells : %d, %d, %d, %d\n",
2355         //                     clus->E(),tof,patches.At(iabsId), absIDCell[0],absIDCell[1],absIDCell[2],absIDCell[3]);
2356
2357         for(Int_t ipatch = 0; ipatch < 4; ipatch++)
2358         {
2359           if(absIdMax == absIDCell[ipatch])
2360           {
2361             //printf("*** Patches : absId %d, E %2.1f, tof %f \n",absIdMax,clus->E(), tof);
2362             if(energy > ePatchMax)
2363             {
2364               tofPatchMax          = tof;
2365               ePatchMax            = energy;
2366               fIsBadCellEvent      = bad;
2367               fIsExoticEvent       = exotic;
2368               fTriggerClusterIndex = iclus;
2369               fTriggerClusterId    = idclus;
2370               fIsTriggerMatch      = kTRUE;
2371             }
2372           }
2373         }// cell patch loop
2374       }// trigger patch loop
2375     } // Do trigger patch matching
2376
2377   }// Cluster loop
2378
2379   // If there was no match, assign as trigger
2380   // the highest energy cluster in the event
2381   if(!fIsTriggerMatch)
2382   {
2383     tofPatchMax          = tofMax;
2384     ePatchMax            = eMax;
2385     fIsBadCellEvent      = badMax;
2386     fIsExoticEvent       = exoMax;
2387     fTriggerClusterIndex = clusMax;
2388     fTriggerClusterId    = idclusMax;
2389   }
2390
2391   Double_t tofPatchMaxUS = TMath::Abs(tofPatchMax);
2392
2393   if     (tofPatchMaxUS < 25 ) fTriggerClusterBC = 0 ;
2394   else if(tofPatchMaxUS < 75 ) fTriggerClusterBC = 1 ;
2395   else if(tofPatchMaxUS < 125) fTriggerClusterBC = 2 ;
2396   else if(tofPatchMaxUS < 175) fTriggerClusterBC = 3 ;
2397   else if(tofPatchMaxUS < 225) fTriggerClusterBC = 4 ;
2398   else if(tofPatchMaxUS < 275) fTriggerClusterBC = 5 ;
2399   else
2400   {
2401     //printf("AliCaloTrackReader::MatchTriggerCluster() - Large BC - tof %2.3f - Index %d\n",tofPatchMaxUS,fTriggerClusterIndex);
2402     if(fTriggerClusterIndex >= 0) fTriggerClusterBC = 6 ;
2403     else
2404     {
2405       fTriggerClusterIndex = -2;
2406       fTriggerClusterId    = -2;
2407     }
2408   }
2409
2410   if(tofPatchMax < 0) fTriggerClusterBC*=-1;
2411
2412 //  printf("AliCaloTrackReader::MatchTriggerCluster(TArrayI patches) - Trigger cluster: index %d, ID %d, E = %2.2f, tof = %2.2f (BC = %d), bad cell? %d, exotic? %d, patch match? %d, n High E cluster %d\n",
2413 //         fTriggerClusterIndex, fTriggerClusterId,ePatchMax, tofPatchMax,
2414 //         fTriggerClusterBC, fIsBadCellEvent,fIsExoticEvent, fIsTriggerMatch, nOfHighECl);
2415 //
2416 //  if(!fIsTriggerMatch)  printf("\t highest energy cluster:  index %d, ID %d, E = %2.2f, tof = %2.2f, bad cell? %d, exotic? %d\n",clusMax, idclusMax, eMax,tofMax, badMax,exoMax);
2417
2418   if(fIsBadCellEvent) fIsExoticEvent = kFALSE;
2419
2420 }
2421
2422 //__________________________________________
2423 Bool_t  AliCaloTrackReader::RejectLEDEvents()
2424 {
2425   // LED Events in period LHC11a contaminated sample, simple method
2426   // to reject such events
2427
2428   // Count number of cells with energy larger than 0.1 in SM3, cut on this number
2429   Int_t ncellsSM3 = 0;
2430   for(Int_t icell = 0; icell < fInputEvent->GetEMCALCells()->GetNumberOfCells(); icell++)
2431   {
2432     Int_t absID = fInputEvent->GetEMCALCells()->GetCellNumber(icell);
2433     Int_t sm    = GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absID);
2434     if(fInputEvent->GetEMCALCells()->GetAmplitude(icell) > 0.1 && sm==3) ncellsSM3++;
2435   }
2436
2437   Int_t ncellcut = 21;
2438   if(GetFiredTriggerClasses().Contains("EMC")) ncellcut = 35;
2439
2440   if(ncellsSM3 >= ncellcut)
2441   {
2442     if(fDebug > 0)
2443       printf(" AliCaloTrackReader::FillInputEvent() - reject event with ncells in SM3 %d, cut %d, trig %s\n",
2444              ncellsSM3,ncellcut,GetFiredTriggerClasses().Data());
2445     return kTRUE;
2446   }
2447
2448   return kFALSE;
2449
2450 }
2451
2452 //____________________________________________________________
2453 void  AliCaloTrackReader::SetTrackCuts(AliESDtrackCuts * cuts)
2454 {
2455   // Set Track cuts
2456
2457   if(fESDtrackCuts) delete fESDtrackCuts ;
2458
2459   fESDtrackCuts = cuts ;
2460
2461 }
2462
2463 //_________________________________________________________________________
2464 void  AliCaloTrackReader::SetTrackComplementaryCuts(AliESDtrackCuts * cuts)
2465 {
2466   // Set Track cuts for complementary tracks (hybrids)
2467
2468   if(fESDtrackComplementaryCuts) delete fESDtrackComplementaryCuts ;
2469
2470   fESDtrackComplementaryCuts = cuts ;
2471
2472 }
2473
2474