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