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