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