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