PWG/CaloTrackCorrBase/AliCaloTrackReader.cxx

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