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