PWG/CaloTrackCorrBase/AliCaloTrackReader.cxx

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