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