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