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