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