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
  31 // ---- ANALYSIS system ----
  32 #include "AliMCEvent.h"
  33 #include "AliAODMCHeader.h"
  34 #include "AliGenPythiaEventHeader.h"
  35 #include "AliESDEvent.h"
  36 #include "AliAODEvent.h"
  37 #include "AliVTrack.h"
  38 #include "AliVParticle.h"
  39 #include "AliMixedEvent.h"
  40 #include "AliESDtrack.h"
  41 #include "AliESDtrackCuts.h"
  42 #include "AliTriggerAnalysis.h"
  43 #include "AliESDVZERO.h"
  44 #include "AliVCaloCells.h"
  45 #include "AliAnalysisManager.h"
  46 #include "AliInputEventHandler.h"
  47 #include "AliAODMCParticle.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 fEventTriggerMask(0),        fMixEventTriggerMask(0),         fEventTriggerAtSE(0),
  89 fAnaLED(kFALSE),
  90 fTaskName(""),               fCaloUtils(0x0),
  91 fMixedEvent(NULL),           fNMixedEvent(0),                 fVertex(NULL),
  92 fListMixedTracksEvents(),    fListMixedCaloEvents(),
  93 fLastMixedTracksEvent(-1),   fLastMixedCaloEvent(-1),
  94 fWriteOutputDeltaAOD(kFALSE),fOldAOD(kFALSE),                 fCaloFilterPatch(kFALSE),
  95 fEMCALClustersListName(""),  fZvtxCut(0.),
  96 fAcceptFastCluster(kFALSE),  fRemoveLEDEvents(kTRUE),
  97 fDoEventSelection(kFALSE),   fDoV0ANDEventSelection(kFALSE),
  98 fDoVertexBCEventSelection(kFALSE),
  99 fDoRejectNoTrackEvents(kFALSE),
 100 fUseEventsWithPrimaryVertex(kFALSE),
 101 fTriggerAnalysis (0x0),      fTimeStampEventSelect(0),
 102 fTimeStampEventFracMin(0),   fTimeStampEventFracMax(0),
 103 fTimeStampRunMin(0),         fTimeStampRunMax(0),
 104 fNPileUpClusters(-1),        fNNonPileUpClusters(-1),         fNPileUpClustersCut(3),
 105 fVertexBC(-200),             fRecalculateVertexBC(0),
 106 fCentralityClass(""),            fCentralityOpt(0),
 107 fEventPlaneMethod(""),       fImportGeometryFromFile(kFALSE), fImportGeometryFilePath("")
 108 {
 109   //Ctor
 110
 111   //Initialize parameters
 112   InitParameters();
 113 }
 114
 115 //_______________________________________
 116 AliCaloTrackReader::~AliCaloTrackReader()
 117 {
 118   //Dtor
 119
 120   delete fFiducialCut ;
 121
 122   if(fAODBranchList)
 123   {
 124     fAODBranchList->Delete();
 125     delete fAODBranchList ;
 126   }
 127
 128   if(fCTSTracks)
 129   {
 130     if(fDataType!=kMC)fCTSTracks->Clear() ;
 131     else              fCTSTracks->Delete() ;
 132     delete fCTSTracks ;
 133   }
 134
 135   if(fEMCALClusters)
 136   {
 137     if(fDataType!=kMC)fEMCALClusters->Clear("C") ;
 138     else              fEMCALClusters->Delete() ;
 139     delete fEMCALClusters ;
 140   }
 141
 142   if(fPHOSClusters)
 143   {
 144     if(fDataType!=kMC)fPHOSClusters->Clear("C") ;
 145     else              fPHOSClusters->Delete() ;
 146     delete fPHOSClusters ;
 147   }
 148
 149   if(fVertex)
 150   {
 151     for (Int_t i = 0; i < fNMixedEvent; i++)
 152     {
 153       delete [] fVertex[i] ;
 154
 155     }
 156     delete [] fVertex ;
 157         }
 158
 159   delete fESDtrackCuts;
 160   delete fESDtrackComplementaryCuts;
 161   delete fTriggerAnalysis;
 162
 163   //  Pointers not owned, done by the analysis frame
 164   //  if(fInputEvent)  delete fInputEvent ;
 165   //  if(fOutputEvent) delete fOutputEvent ;
 166   //  if(fMC)          delete fMC ;
 167   //  Pointer not owned, deleted by maker
 168   //  if (fCaloUtils) delete fCaloUtils ;
 169
 170 }
 171
 172 //________________________________________________________________________
 173 Bool_t  AliCaloTrackReader::AcceptDCA(const Float_t pt, const Float_t dca)
 174 {
 175   // Accept track if DCA is smaller than function
 176
 177   Float_t cut = fTrackDCACut[0]+fTrackDCACut[1]/TMath::Power(pt,fTrackDCACut[2]);
 178
 179   if(TMath::Abs(dca) < cut)
 180     return kTRUE;
 181   else
 182     return kFALSE;
 183
 184 }
 185
 186 //________________________________________________
 187 Bool_t AliCaloTrackReader::ComparePtHardAndJetPt()
 188 {
 189   // Check the event, if the requested ptHard is much smaller than the jet pT, then there is a problem.
 190   // Only for PYTHIA.
 191
 192   //printf("AliCaloTrackReader::ComparePtHardAndJetPt() - GenHeaderName : %s\n",GetGenEventHeader()->ClassName());
 193
 194   if(!strcmp(GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader"))
 195   {
 196     TParticle * jet =  0;
 197     AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetGenEventHeader();
 198     Int_t nTriggerJets =  pygeh->NTriggerJets();
 199     Float_t ptHard = pygeh->GetPtHard();
 200
 201     if(fDebug > 1)
 202       printf("AliCaloTrackReader::ComparePtHardAndJetPt() - Njets: %d, pT Hard %f\n",nTriggerJets, ptHard);
 203
 204     Float_t tmpjet[]={0,0,0,0};
 205     for(Int_t ijet = 0; ijet< nTriggerJets; ijet++)
 206     {
 207       pygeh->TriggerJet(ijet, tmpjet);
 208       jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
 209
 210       if(fDebug > 1)
 211         printf("AliCaloTrackReader::ComparePtHardAndJetPt() - jet %d; pycell jet pT %f\n",ijet, jet->Pt());
 212
 213       //Compare jet pT and pt Hard
 214       if(jet->Pt() > fPtHardAndJetPtFactor * ptHard)
 215       {
 216         printf("AliCaloTrackReader::ComparePtHardAndJetPt() - Reject jet event with : pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n",
 217                 ptHard, jet->Pt(), fPtHardAndJetPtFactor);
 218         return kFALSE;
 219       }
 220     }
 221
 222     if(jet) delete jet;
 223   }
 224
 225   return kTRUE ;
 226
 227 }
 228
 229 //____________________________________________________________________
 230 Bool_t AliCaloTrackReader::ComparePtHardAndClusterPt()
 231 {
 232   // Check the event, if the requested ptHard is smaller than the calorimeter cluster E, then there is a problem.
 233   // Only for PYTHIA.
 234
 235   if(!strcmp(GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader"))
 236   {
 237     AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetGenEventHeader();
 238     Float_t ptHard = pygeh->GetPtHard();
 239
 240     Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
 241     for (Int_t iclus =  0; iclus <  nclusters; iclus++)
 242     {
 243       AliVCluster * clus = fInputEvent->GetCaloCluster(iclus) ;
 244       Float_t ecluster = clus->E();
 245
 246       if(ecluster > fPtHardAndClusterPtFactor*ptHard)
 247       {
 248         printf("AliCaloTrackReader::ComparePtHardAndClusterPt() - Reject : ecluster %2.2f, calo %d, factor %2.2f, ptHard %f\n",ecluster,clus->GetType(),fPtHardAndClusterPtFactor,ptHard);
 249
 250         return kFALSE;
 251       }
 252     }
 253
 254   }
 255
 256   return kTRUE ;
 257
 258 }
 259
 260 //_________________________________________________________________
 261 void AliCaloTrackReader::CorrectMCLabelForAODs(AliVCluster * clus)
 262 {
 263  // AODs filter particles, not anymore correspondance with MC position in array
 264  // Check if label is correct and if not, change it
 265
 266   Int_t * labels = clus->GetLabels();
 267
 268   for(UInt_t ilabel = 0; ilabel < clus->GetNLabels(); ilabel++)
 269   {
 270     Int_t orgLabel = labels[ilabel];
 271
 272     TClonesArray * arr = GetAODMCParticles()  ;
 273
 274     if(!arr)
 275     {
 276       printf("AliCaloTrackReader::CorrectMCLabelForAODs() - Input array not available\n");
 277       return ;
 278     }
 279
 280     AliAODMCParticle *  particle = (AliAODMCParticle *)arr->At(orgLabel);
 281
 282     if(orgLabel != particle->Label())
 283     {
 284       // loop on the particles list and check if there is one with the same label
 285       for(Int_t ind = 0; ind < arr->GetEntriesFast(); ind++ )
 286       {
 287         particle = (AliAODMCParticle *) arr->At(ind);
 288
 289         if(orgLabel == particle->Label()) labels[ilabel] = ind;
 290       }
 291     }
 292
 293     //if(orgLabel!=labels[ilabel]) printf("\t Label in %d - out %d \n",orgLabel, clus->GetLabels()[ilabel]);
 294
 295   }
 296 }
 297
 298 //____________________________________________
 299 AliStack* AliCaloTrackReader::GetStack() const
 300 {
 301   //Return pointer to stack
 302   if(fMC)
 303     return fMC->Stack();
 304   else
 305   {
 306     if(fDebug > 1) printf("AliCaloTrackReader::GetStack() - Stack is not available\n");
 307     return 0x0 ;
 308   }
 309 }
 310
 311 //__________________________________________________
 312 TString AliCaloTrackReader::GetFiredTriggerClasses()
 313 {
 314   // List of triggered classes in a TString
 315
 316   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (GetInputEvent());
 317   AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (GetInputEvent());
 318
 319   if     (esdevent) return esdevent->GetFiredTriggerClasses();
 320   else if(aodevent) return aodevent->GetFiredTriggerClasses();
 321   else              return ""; // Mixed Event, MC event, does not have this trigger info
 322
 323 }
 324
 325 //______________________________________________
 326 AliHeader* AliCaloTrackReader::GetHeader() const
 327 {
 328   //Return pointer to header
 329   if(fMC)
 330   {
 331     return fMC->Header();
 332   }
 333   else
 334   {
 335     printf("AliCaloTrackReader::Header is not available\n");
 336     return 0x0 ;
 337   }
 338 }
 339
 340 //______________________________________________________________
 341 AliGenEventHeader* AliCaloTrackReader::GetGenEventHeader() const
 342 {
 343   //Return pointer to Generated event header
 344   if     (ReadStack() && fMC)
 345   {
 346     return fMC->GenEventHeader();
 347   }
 348   else if(ReadAODMCParticles() && GetAODMCHeader())
 349   {
 350     //printf("AliCaloTrackReader::GetGenEventHeader() - N headers %d\n",GetAODMCHeader()->GetNCocktailHeaders());
 351     if( GetAODMCHeader()->GetNCocktailHeaders() > 0)
 352       return GetAODMCHeader()->GetCocktailHeader(0) ;
 353     else
 354       return 0x0;
 355   }
 356   else
 357   {
 358     //printf("AliCaloTrackReader::GetGenEventHeader() - MC header not available! \n");
 359     return 0;
 360   }
 361 }
 362
 363 //____________________________________________________________________
 364 TClonesArray* AliCaloTrackReader::GetAODMCParticles() const
 365 {
 366   //Return list of particles in AOD. Do it for the corresponding input event.
 367
 368   TClonesArray * rv = NULL ;
 369   if(fDataType == kAOD)
 370   {
 371     //Normal input AOD
 372     AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fInputEvent) ;
 373     if(evt)
 374       rv = (TClonesArray*)evt->FindListObject("mcparticles");
 375     else
 376       printf("AliCaloTrackReader::GetAODMCParticles() - Null AOD event \n");
 377   }
 378   else
 379   {
 380     printf("AliCaloTrackReader::GetAODMCParticles() - Input are not AODs\n");
 381   }
 382
 383   return rv ;
 384 }
 385
 386 //________________________________________________________
 387 AliAODMCHeader* AliCaloTrackReader::GetAODMCHeader() const
 388 {
 389   //Return MC header in AOD. Do it for the corresponding input event.
 390
 391   AliAODMCHeader *mch = NULL;
 392
 393   if(fDataType == kAOD)
 394   {
 395     AliAODEvent * aod = dynamic_cast<AliAODEvent*> (fInputEvent);
 396     if(aod) mch = dynamic_cast<AliAODMCHeader*>(aod->FindListObject("mcHeader"));
 397   }
 398   else
 399   {
 400     printf("AliCaloTrackReader::GetAODMCHeader() - Input are not AODs\n");
 401   }
 402
 403   return mch;
 404 }
 405
 406 //___________________________________________________________
 407 Int_t AliCaloTrackReader::GetVertexBC(const AliVVertex * vtx)
 408 {
 409   // Get the vertex BC
 410
 411   Int_t vertexBC=vtx->GetBC();
 412   if(!fRecalculateVertexBC) return vertexBC;
 413
 414   // In old AODs BC not stored, recalculate it
 415   // loop over the global track and select those which have small DCA to primary vertex (e.g. primary).
 416   // If at least one of these primaries has valid BC != 0, then this vertex is a pile-up candidate.
 417   // Execute after CTS
 418   Double_t bz  = fInputEvent->GetMagneticField();
 419   Bool_t   bc0 = kFALSE;
 420   Int_t    ntr = GetCTSTracks()->GetEntriesFast();
 421   //printf("N Tracks %d\n",ntr);
 422
 423   for(Int_t i = 0 ; i < ntr ; i++)
 424   {
 425     AliVTrack * track =  (AliVTrack*) (GetCTSTracks()->At(i));
 426
 427     //Check if has TOF info, if not skip
 428     ULong_t status  = track->GetStatus();
 429     Bool_t  okTOF   = (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ;
 430     vertexBC        = track->GetTOFBunchCrossing(bz);
 431     Float_t pt      = track->Pt();
 432
 433     if(!okTOF) continue;
 434
 435     // Get DCA x, y
 436     Double_t dca[2]   = {1e6,1e6};
 437     Double_t covar[3] = {1e6,1e6,1e6};
 438     track->PropagateToDCA(vtx,bz,100.,dca,covar);
 439
 440     if(AcceptDCA(pt,dca[0]))
 441     {
 442       if     (vertexBC !=0 && fVertexBC != AliVTrack::kTOFBCNA) return vertexBC;
 443       else if(vertexBC == 0)                                    bc0 = kTRUE;
 444     }
 445   }
 446
 447   if( bc0 ) vertexBC = 0 ;
 448   else      vertexBC = AliVTrack::kTOFBCNA ;
 449
 450   return vertexBC;
 451
 452 }
 453
 454 //_____________________________
 455 void AliCaloTrackReader::Init()
 456 {
 457   //Init reader. Method to be called in AliAnaPartCorrMaker
 458
 459   //printf(" AliCaloTrackReader::Init() %p \n",gGeoManager);
 460
 461   if(fReadStack && fReadAODMCParticles)
 462   {
 463     printf("AliCaloTrackReader::Init() - Cannot access stack and mcparticles at the same time, change them \n");
 464     fReadStack          = kFALSE;
 465     fReadAODMCParticles = kFALSE;
 466   }
 467
 468   // Init geometry, I do not like much to do it like this ...
 469   if(fImportGeometryFromFile && !gGeoManager)
 470   {
 471     if(fImportGeometryFilePath=="") // If not specified, set a default location
 472     fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2011.root"; // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
 473
 474     printf("AliCaloTrackReader::Init() - Import %s\n",fImportGeometryFilePath.Data());
 475     TGeoManager::Import(fImportGeometryFilePath) ; // default need file "geometry.root" in local dir!!!!
 476   }
 477
 478   if(!fESDtrackCuts)
 479     fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts(); //initialize with TPC only tracks
 480
 481 }
 482
 483 //_______________________________________
 484 void AliCaloTrackReader::InitParameters()
 485 {
 486   //Initialize the parameters of the analysis.
 487   fDataType   = kESD ;
 488   fCTSPtMin   = 0.1 ;
 489   fEMCALPtMin = 0.1 ;
 490   fPHOSPtMin  = 0.1 ;
 491   fCTSPtMax   = 1000. ;
 492   fEMCALPtMax = 1000. ;
 493   fPHOSPtMax  = 1000. ;
 494
 495   //Track DCA cuts
 496   // dca_xy cut = 0.0105+0.0350/TMath::Power(pt,1.1);
 497   fTrackDCACut[0] = 0.0105;
 498   fTrackDCACut[1] = 0.0350;
 499   fTrackDCACut[2] = 1.1;
 500
 501   //Do not filter the detectors input by default.
 502   fFillEMCAL      = kFALSE;
 503   fFillPHOS       = kFALSE;
 504   fFillCTS        = kFALSE;
 505   fFillEMCALCells = kFALSE;
 506   fFillPHOSCells  = kFALSE;
 507
 508   fReadStack             = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
 509   fReadAODMCParticles    = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
 510   fDeltaAODFileName      = "deltaAODPartCorr.root";
 511   fFiredTriggerClassName = "";
 512   fEventTriggerMask      = AliVEvent::kAny;
 513   fMixEventTriggerMask   = AliVEvent::kAnyINT;
 514   fEventTriggerAtSE      = kTRUE; // Use only events that pass event selection at SE base class
 515
 516   fAcceptFastCluster = kTRUE;
 517   fAnaLED            = kFALSE;
 518
 519   //We want tracks fitted in the detectors:
 520   //fTrackStatus=AliESDtrack::kTPCrefit;
 521   //fTrackStatus|=AliESDtrack::kITSrefit;
 522   fTrackStatus     = 0;
 523   fTrackFilterMask = 128; //For AODs, but what is the difference between fTrackStatus and fTrackFilterMask?
 524
 525   fESDtrackCuts = 0;
 526   fESDtrackComplementaryCuts = 0;
 527
 528   fConstrainTrack = kFALSE ; // constrain tracks to vertex
 529
 530   fV0ADC[0] = 0;   fV0ADC[1] = 0;
 531   fV0Mul[0] = 0;   fV0Mul[1] = 0;
 532
 533   fZvtxCut   = 10.;
 534
 535   fNMixedEvent = 1;
 536
 537   fPtHardAndJetPtFactor     = 7.;
 538   fPtHardAndClusterPtFactor = 1.;
 539
 540   //Centrality
 541   fCentralityClass  = "V0M";
 542   fCentralityOpt    = 10;
 543   fCentralityBin[0] = fCentralityBin[1]=-1;
 544
 545   fEventPlaneMethod = "V0";
 546
 547   // Allocate memory (not sure this is the right place)
 548   fCTSTracks       = new TObjArray();
 549   fEMCALClusters   = new TObjArray();
 550   fPHOSClusters    = new TObjArray();
 551   fTriggerAnalysis = new AliTriggerAnalysis;
 552   fAODBranchList   = new TList ;
 553
 554   fImportGeometryFromFile = kFALSE;
 555
 556   fPileUpParamSPD[0] = 3   ; fPileUpParamSPD[1] = 0.8 ;
 557   fPileUpParamSPD[2] = 3.0 ; fPileUpParamSPD[3] = 2.0 ; fPileUpParamSPD[4] = 5.0;
 558
 559   // Parametrized time cut (LHC11d)
 560   fEMCALParamTimeCutMin[0] =-5; fEMCALParamTimeCutMin[1] =-1 ; fEMCALParamTimeCutMin[2] = 3.5 ; fEMCALParamTimeCutMin[3] = 1.  ;
 561   fEMCALParamTimeCutMax[0] = 5; fEMCALParamTimeCutMax[1] = 50; fEMCALParamTimeCutMax[2] = 0.45; fEMCALParamTimeCutMax[3] = 1.25;
 562
 563   // Parametrized time cut (LHC11c)
 564   //fEMCALParamTimeCutMin[0] =-5;   fEMCALParamTimeCutMin[1] =-1 ; fEMCALParamTimeCutMin[2] = 1.87; fEMCALParamTimeCutMin[3] = 0.4;
 565   //fEMCALParamTimeCutMax[0] = 3.5; fEMCALParamTimeCutMax[1] = 50; fEMCALParamTimeCutMax[2] = 0.15; fEMCALParamTimeCutMax[3] = 1.6;
 566
 567   fTimeStampRunMin = -1;
 568   fTimeStampRunMax = 1e12;
 569   fTimeStampEventFracMin = -1;
 570   fTimeStampEventFracMax = 2;
 571
 572   for(Int_t i = 0; i < 19; i++)
 573   {
 574     fEMCalBCEvent   [i] = 0;
 575     fEMCalBCEventCut[i] = 0;
 576     fTrackBCEvent   [i] = 0;
 577     fTrackBCEventCut[i] = 0;
 578   }
 579
 580 }
 581
 582 //___________________________________________________________
 583 Bool_t AliCaloTrackReader::IsInTimeWindow(const Double_t tof, const Float_t energy) const
 584 {
 585   // Cluster time selection window
 586
 587   // Parametrized cut depending on E
 588   if(fUseParamTimeCut)
 589   {
 590     Float_t minCut= fEMCALParamTimeCutMin[0]+fEMCALParamTimeCutMin[1]*TMath::Exp(-(energy-fEMCALParamTimeCutMin[2])/fEMCALParamTimeCutMin[3]);
 591     Float_t maxCut= fEMCALParamTimeCutMax[0]+fEMCALParamTimeCutMax[1]*TMath::Exp(-(energy-fEMCALParamTimeCutMax[2])/fEMCALParamTimeCutMax[3]);
 592     //printf("tof %f, minCut %f, maxCut %f\n",tof,minCut,maxCut);
 593     if( tof < minCut || tof > maxCut )  return kFALSE ;
 594   }
 595
 596   //In any case, the time should to be larger than the fixed window ...
 597   if( tof < fEMCALTimeCutMin  || tof > fEMCALTimeCutMax )  return kFALSE ;
 598
 599   return kTRUE ;
 600 }
 601
 602 //________________________________________________
 603 Bool_t AliCaloTrackReader::IsPileUpFromSPD() const
 604 {
 605   // Check if event is from pile-up determined by SPD
 606   // Default values: (3, 0.8, 3., 2., 5.)
 607   return fInputEvent->IsPileupFromSPD((Int_t) fPileUpParamSPD[0] , fPileUpParamSPD[1] ,
 608                                               fPileUpParamSPD[2] , fPileUpParamSPD[3] , fPileUpParamSPD[4] );
 609   //printf("Param : %d, %2.2f, %2.2f, %2.2f, %2.2f\n",(Int_t) fPileUpParamSPD[0], fPileUpParamSPD[1], fPileUpParamSPD[2], fPileUpParamSPD[3], fPileUpParamSPD[4]);
 610
 611 }
 612
 613 //__________________________________________________
 614 Bool_t AliCaloTrackReader::IsPileUpFromEMCal() const
 615 {
 616   // Check if event is from pile-up determined by EMCal
 617   if(fNPileUpClusters > fNPileUpClustersCut) return kTRUE ;
 618   else                                       return kFALSE;
 619 }
 620
 621 //________________________________________________________
 622 Bool_t AliCaloTrackReader::IsPileUpFromSPDAndEMCal() const
 623 {
 624   // Check if event is from pile-up determined by SPD and EMCal
 625   if( IsPileUpFromSPD() && IsPileUpFromEMCal()) return kTRUE ;
 626   else                                          return kFALSE;
 627 }
 628
 629 //_______________________________________________________
 630 Bool_t AliCaloTrackReader::IsPileUpFromSPDOrEMCal() const
 631 {
 632   // Check if event is from pile-up determined by SPD or EMCal
 633   if( IsPileUpFromSPD() || IsPileUpFromEMCal()) return kTRUE ;
 634   else                                          return kFALSE;
 635 }
 636
 637 //___________________________________________________________
 638 Bool_t AliCaloTrackReader::IsPileUpFromSPDAndNotEMCal() const
 639 {
 640   // Check if event is from pile-up determined by SPD and not by EMCal
 641   if( IsPileUpFromSPD() && !IsPileUpFromEMCal()) return kTRUE ;
 642   else                                          return kFALSE;
 643 }
 644
 645 //___________________________________________________________
 646 Bool_t AliCaloTrackReader::IsPileUpFromEMCalAndNotSPD() const
 647 {
 648   // Check if event is from pile-up determined by EMCal, not by SPD
 649   if( !IsPileUpFromSPD() && IsPileUpFromEMCal()) return kTRUE ;
 650   else                                           return kFALSE;
 651 }
 652
 653 //______________________________________________________________
 654 Bool_t AliCaloTrackReader::IsPileUpFromNotSPDAndNotEMCal() const
 655 {
 656   // Check if event not from pile-up determined neither by SPD nor by EMCal
 657   if( !IsPileUpFromSPD() && !IsPileUpFromEMCal()) return kTRUE ;
 658   else                                            return kFALSE;
 659 }
 660
 661 //________________________________________________________
 662 void AliCaloTrackReader::Print(const Option_t * opt) const
 663 {
 664
 665   //Print some relevant parameters set for the analysis
 666   if(! opt)
 667     return;
 668
 669   printf("***** Print: %s %s ******\n",    GetName(), GetTitle() ) ;
 670   printf("Task name      : %s\n",          fTaskName.Data()) ;
 671   printf("Data type      : %d\n",          fDataType) ;
 672   printf("CTS Min pT     : %2.1f GeV/c\n", fCTSPtMin) ;
 673   printf("EMCAL Min pT   : %2.1f GeV/c\n", fEMCALPtMin) ;
 674   printf("PHOS Min pT    : %2.1f GeV/c\n", fPHOSPtMin) ;
 675   printf("CTS Max pT     : %2.1f GeV/c\n", fCTSPtMax) ;
 676   printf("EMCAL Max pT   : %2.1f GeV/c\n", fEMCALPtMax) ;
 677   printf("PHOS Max pT    : %2.1f GeV/c\n", fPHOSPtMax) ;
 678   printf("EMCAL Time Cut: %3.1f < TOF  < %3.1f\n", fEMCALTimeCutMin, fEMCALTimeCutMax);
 679   printf("Use CTS         =     %d\n",     fFillCTS) ;
 680   printf("Use EMCAL       =     %d\n",     fFillEMCAL) ;
 681   printf("Use PHOS        =     %d\n",     fFillPHOS) ;
 682   printf("Use EMCAL Cells =     %d\n",     fFillEMCALCells) ;
 683   printf("Use PHOS  Cells =     %d\n",     fFillPHOSCells) ;
 684   printf("Track status    =     %d\n", (Int_t) fTrackStatus) ;
 685   printf("AODs Track filter mask  =  %d or hybrid %d, SPD hit %d\n", (Int_t) fTrackFilterMask,fSelectHybridTracks,fSelectSPDHitTracks) ;
 686   printf("Track Mult Eta Cut =  %d\n", (Int_t) fTrackMultEtaCut) ;
 687   printf("Write delta AOD =     %d\n",     fWriteOutputDeltaAOD) ;
 688   printf("Recalculate Clusters = %d\n",    fRecalculateClusters) ;
 689
 690   printf("Use Triggers selected in SE base class %d; If not what trigger Mask? %d; Trigger max for mixed %d \n",
 691          fEventTriggerAtSE, fEventTriggerMask,fMixEventTriggerMask);
 692
 693   if(fComparePtHardAndClusterPt)
 694           printf("Compare jet pt and pt hard to accept event, factor = %2.2f",fPtHardAndJetPtFactor);
 695
 696   if(fComparePtHardAndClusterPt)
 697           printf("Compare cluster pt and pt hard to accept event, factor = %2.2f",fPtHardAndClusterPtFactor);
 698
 699   printf("Read Kine from, stack? %d, AOD ? %d \n", fReadStack, fReadAODMCParticles) ;
 700   printf("Delta AOD File Name =     %s\n", fDeltaAODFileName.Data()) ;
 701   printf("Centrality: Class %s, Option %d, Bin [%d,%d] \n", fCentralityClass.Data(),fCentralityOpt,fCentralityBin[0], fCentralityBin[1]) ;
 702
 703   printf("    \n") ;
 704
 705 }
 706
 707 //_________________________________________________________________________
 708 Bool_t AliCaloTrackReader::FillInputEvent(const Int_t iEntry,
 709                                           const char * /*currentFileName*/)
 710 {
 711   //Fill the event counter and input lists that are needed, called by the analysis maker.
 712
 713   fEventNumber = iEntry;
 714   //fCurrentFileName = TString(currentFileName);
 715   if(!fInputEvent)
 716   {
 717           if(fDebug >= 0) printf("AliCaloTrackReader::FillInputEvent() - Input event not available, skip event analysis\n");
 718           return kFALSE;
 719   }
 720
 721   //Select events only fired by a certain trigger configuration if it is provided
 722   Int_t eventType = 0;
 723   if(fInputEvent->GetHeader())
 724           eventType = ((AliVHeader*)fInputEvent->GetHeader())->GetEventType();
 725
 726   if (GetFiredTriggerClasses().Contains("FAST")  && !GetFiredTriggerClasses().Contains("ALL") && !fAcceptFastCluster)
 727   {
 728     if(fDebug > 0)  printf("AliCaloTrackReader::FillInputEvent - Do not count events from fast cluster, trigger name %s\n",fFiredTriggerClassName.Data());
 729     return kFALSE;
 730   }
 731
 732   //-------------------------------------------------------------------------------------
 733   // Reject event if large clusters with large energy
 734   // Use only for LHC11a data for the moment, and if input is clusterizer V1 or V1+unfolding
 735   // If clusterzer NxN or V2 it does not help
 736   //-------------------------------------------------------------------------------------
 737   Int_t run = fInputEvent->GetRunNumber();
 738   if( fRemoveLEDEvents && run > 146857  && run < 146861 )
 739   {
 740     //printf("Event %d\n",GetEventNumber());
 741
 742     // Count number of cells with energy larger than 0.1 in SM3, cut on this number
 743     Int_t ncellsSM3 = 0;
 744     for(Int_t icell = 0; icell < fInputEvent->GetEMCALCells()->GetNumberOfCells(); icell++)
 745     {
 746       Int_t absID = fInputEvent->GetEMCALCells()->GetCellNumber(icell);
 747       Int_t sm    = GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absID);
 748       if(fInputEvent->GetEMCALCells()->GetAmplitude(icell) > 0.1 && sm==3) ncellsSM3++;
 749     }
 750
 751     Int_t ncellcut = 21;
 752     if(fFiredTriggerClassName.Contains("EMC")) ncellcut = 35;
 753
 754     if(ncellsSM3 >= ncellcut)
 755     {
 756       if(fDebug > 0) printf(" AliCaloTrackReader::FillInputEvent() - reject event with ncells in SM3 %d\n",ncellsSM3);
 757       return kFALSE;
 758     }
 759   }// Remove LED events
 760
 761   // Reject pure LED events?
 762   if( fFiredTriggerClassName  !="" && !fAnaLED)
 763   {
 764     //printf("Event type %d\n",eventType);
 765     if(eventType!=7)
 766       return kFALSE; //Only physics event, do not use for simulated events!!!
 767
 768     if(fDebug > 0)
 769       printf("AliCaloTrackReader::FillInputEvent() - FiredTriggerClass <%s>, selected class <%s>, compare name %d\n",
 770              GetFiredTriggerClasses().Data(),fFiredTriggerClassName.Data(), GetFiredTriggerClasses().Contains(fFiredTriggerClassName));
 771
 772     if( !GetFiredTriggerClasses().Contains(fFiredTriggerClassName) ) return kFALSE;
 773     else if(fDebug > 0) printf("AliCaloTrackReader::FillInputEvent() - Accepted triggered event\n");
 774   }
 775   else if(fAnaLED)
 776   {
 777     //    kStartOfRun =       1,    // START_OF_RUN
 778     //    kEndOfRun =         2,    // END_OF_RUN
 779     //    kStartOfRunFiles =  3,    // START_OF_RUN_FILES
 780     //    kEndOfRunFiles =    4,    // END_OF_RUN_FILES
 781     //    kStartOfBurst =     5,    // START_OF_BURST
 782     //    kEndOfBurst =       6,    // END_OF_BURST
 783     //    kPhysicsEvent =     7,    // PHYSICS_EVENT
 784     //    kCalibrationEvent = 8,    // CALIBRATION_EVENT
 785     //    kFormatError =      9,    // EVENT_FORMAT_ERROR
 786     //    kStartOfData =      10,   // START_OF_DATA
 787     //    kEndOfData =        11,   // END_OF_DATA
 788     //    kSystemSoftwareTriggerEvent   = 12, // SYSTEM_SOFTWARE_TRIGGER_EVENT
 789     //    kDetectorSoftwareTriggerEvent = 13  // DETECTOR_SOFTWARE_TRIGGER_EVENT
 790
 791           if(eventType!=7 && fDebug > 1 )printf("AliCaloTrackReader::FillInputEvent() - DO LED, Event Type <%d>, 8 Calibration \n",  eventType);
 792           if(eventType!=8)return kFALSE;
 793   }
 794
 795   //In case of analysis of events with jets, skip those with jet pt > 5 pt hard
 796   if(fComparePtHardAndJetPt)
 797   {
 798     if(!ComparePtHardAndJetPt()) return kFALSE ;
 799   }
 800
 801   if(fComparePtHardAndClusterPt)
 802   {
 803     if(!ComparePtHardAndClusterPt()) return kFALSE ;
 804   }
 805
 806   //Fill Vertex array
 807   FillVertexArray();
 808   //Reject events with Z vertex too large, only for SE analysis, if not, cut on the analysis code
 809   if(!GetMixedEvent() && TMath::Abs(fVertex[0][2]) > fZvtxCut) return kFALSE;
 810
 811   //------------------------------------------------------
 812   //Event rejection depending on vertex, pileup, v0and
 813   //------------------------------------------------------
 814   if(fDataType==kESD && fTimeStampEventSelect)
 815   {
 816     AliESDEvent* esd = dynamic_cast<AliESDEvent*> (fInputEvent);
 817     if(esd)
 818     {
 819       Int_t timeStamp = esd->GetTimeStamp();
 820       Float_t timeStampFrac = 1.*(timeStamp-fTimeStampRunMin) / (fTimeStampRunMax-fTimeStampRunMin);
 821
 822       //printf("stamp0 %d, max0 %d, frac %f\n", timeStamp-fTimeStampRunMin,fTimeStampRunMax-fTimeStampRunMin, timeStampFrac);
 823
 824       if(timeStampFrac < fTimeStampEventFracMin || timeStampFrac > fTimeStampEventFracMax) return kFALSE;
 825     }
 826     //printf("\t accept time stamp\n");
 827   }
 828
 829
 830   //------------------------------------------------------
 831   //Event rejection depending on vertex, pileup, v0and
 832   //------------------------------------------------------
 833
 834   if(fUseEventsWithPrimaryVertex)
 835   {
 836     if( !CheckForPrimaryVertex() )              return kFALSE;
 837     if( TMath::Abs(fVertex[0][0] ) < 1.e-6 &&
 838         TMath::Abs(fVertex[0][1] ) < 1.e-6 &&
 839         TMath::Abs(fVertex[0][2] ) < 1.e-6    ) return kFALSE;
 840   }
 841
 842   //printf("Reader : IsPileUp %d, Multi %d\n",IsPileUpFromSPD(),fInputEvent->IsPileupFromSPDInMultBins());
 843
 844   if(fDoEventSelection)
 845   {
 846     if(!fCaloFilterPatch)
 847     {
 848       // Do not analyze events with pileup
 849       Bool_t bPileup = IsPileUpFromSPD();
 850       //IsPileupFromSPDInMultBins() // method to try
 851       //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]);
 852       if(bPileup) return kFALSE;
 853
 854       if(fDoV0ANDEventSelection)
 855       {
 856         Bool_t bV0AND = kTRUE;
 857         AliESDEvent* esd = dynamic_cast<AliESDEvent*> (fInputEvent);
 858         if(esd)
 859           bV0AND = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0AND);
 860         //else bV0AND = //FIXME FOR AODs
 861         if(!bV0AND) return kFALSE;
 862       }
 863     }//CaloFilter patch
 864     else
 865     {
 866       if(fInputEvent->GetNumberOfCaloClusters() > 0)
 867       {
 868         AliVCluster * calo = fInputEvent->GetCaloCluster(0);
 869         if(calo->GetNLabels() == 4)
 870         {
 871           Int_t * selection = calo->GetLabels();
 872           Bool_t bPileup = selection[0];
 873           if(bPileup) return kFALSE;
 874
 875           Bool_t bGoodV = selection[1];
 876           if(fUseEventsWithPrimaryVertex && !bGoodV) return kFALSE;
 877
 878           if(fDoV0ANDEventSelection)
 879           {
 880             Bool_t bV0AND = selection[2];
 881             if(!bV0AND) return kFALSE;
 882           }
 883
 884           fTrackMult = selection[3];
 885           if(fTrackMult == 0) return kFALSE;
 886         }
 887         else
 888         {
 889           //First filtered AODs, track multiplicity stored there.
 890           fTrackMult = (Int_t) ((AliAODHeader*)fInputEvent->GetHeader())->GetCentrality();
 891           if(fTrackMult == 0) return kFALSE;
 892         }
 893       }//at least one cluster
 894       else
 895       {
 896         //printf("AliCaloTrackReader::FillInputEvent() - No clusters in event\n");
 897         //Remove events with  vertex (0,0,0), bad vertex reconstruction
 898         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;
 899
 900         //First filtered AODs, track multiplicity stored there.
 901         fTrackMult = (Int_t) ((AliAODHeader*)fInputEvent->GetHeader())->GetCentrality();
 902         if(fTrackMult == 0) return kFALSE;
 903       }// no cluster
 904     }// CaloFileter patch
 905   }// Event selection/AliceSoft/AliRoot/trunk/PWG/CaloTrackCorrBase/AliCaloTrackReader.h
 906   //------------------------------------------------------
 907
 908   //Check if there is a centrality value, PbPb analysis, and if a centrality bin selection is requested
 909   //If we need a centrality bin, we select only those events in the corresponding bin.
 910   if(GetCentrality() && fCentralityBin[0]>=0 && fCentralityBin[1]>=0 && fCentralityOpt==100)
 911   {
 912     Int_t cen = GetEventCentrality();
 913     if(cen > fCentralityBin[1] || cen < fCentralityBin[0]) return kFALSE; //reject events out of bin.
 914   }
 915
 916   //Fill the arrays with cluster/tracks/cells data
 917
 918   if(!fEventTriggerAtSE)
 919   {
 920     // In case of mixing analysis, accept MB events, not only Trigger
 921     // Track and cluster arrays filled for MB in order to create the pool in the corresponding analysis
 922     // via de method in the base class FillMixedEventPool()
 923
 924     AliAnalysisManager *manager = AliAnalysisManager::GetAnalysisManager();
 925     AliInputEventHandler *inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
 926
 927     if(!inputHandler) return kFALSE ;  // to content coverity
 928
 929     UInt_t isTrigger = inputHandler->IsEventSelected() & fEventTriggerMask;
 930     UInt_t isMB      = inputHandler->IsEventSelected() & fMixEventTriggerMask;
 931
 932     if(!isTrigger && !isMB) return kFALSE;
 933
 934     //printf("Selected triggered event : %s\n",GetFiredTriggerClasses().Data());
 935   }
 936
 937   // Get the main vertex BC, in case not available
 938   // it is calculated in FillCTS checking the BC of tracks
 939   // with DCA small (if cut applied, if open)
 940   fVertexBC=fInputEvent->GetPrimaryVertex()->GetBC();
 941
 942   if(fFillCTS)
 943   {
 944     FillInputCTS();
 945     //Accept events with at least one track
 946     if(fTrackMult == 0 && fDoRejectNoTrackEvents) return kFALSE ;
 947   }
 948
 949   if(fDoVertexBCEventSelection)
 950   {
 951     if(fVertexBC!=0 && fVertexBC!=AliVTrack::kTOFBCNA) return kFALSE ;
 952   }
 953
 954   if(fFillEMCALCells)
 955     FillInputEMCALCells();
 956
 957   if(fFillPHOSCells)
 958     FillInputPHOSCells();
 959
 960   if(fFillEMCAL)
 961     FillInputEMCAL();
 962
 963   if(fFillPHOS)
 964     FillInputPHOS();
 965
 966   FillInputVZERO();
 967
 968
 969   return kTRUE ;
 970 }
 971
 972 //___________________________________
 973 void AliCaloTrackReader::ResetLists()
 974 {
 975   //  Reset lists, called by the analysis maker
 976
 977   if(fCTSTracks)       fCTSTracks     -> Clear();
 978   if(fEMCALClusters)   fEMCALClusters -> Clear("C");
 979   if(fPHOSClusters)    fPHOSClusters  -> Clear("C");
 980
 981   fV0ADC[0] = 0;   fV0ADC[1] = 0;
 982   fV0Mul[0] = 0;   fV0Mul[1] = 0;
 983
 984 }
 985
 986 //____________________________________________________________
 987 void AliCaloTrackReader::SetInputEvent(AliVEvent* const input)
 988 {
 989   fInputEvent  = input;
 990   fMixedEvent = dynamic_cast<AliMixedEvent*>(GetInputEvent()) ;
 991   if (fMixedEvent)
 992     fNMixedEvent = fMixedEvent->GetNumberOfEvents() ;
 993
 994   //Delete previous vertex
 995   if(fVertex)
 996   {
 997     for (Int_t i = 0; i < fNMixedEvent; i++)
 998     {
 999       delete [] fVertex[i] ;
1000     }
1001     delete [] fVertex ;
1002   }
1003
1004   fVertex = new Double_t*[fNMixedEvent] ;
1005   for (Int_t i = 0; i < fNMixedEvent; i++)
1006   {
1007     fVertex[i] = new Double_t[3] ;
1008     fVertex[i][0] = 0.0 ;
1009     fVertex[i][1] = 0.0 ;
1010     fVertex[i][2] = 0.0 ;
1011   }
1012 }
1013
1014 //__________________________________________________
1015 Int_t AliCaloTrackReader::GetEventCentrality() const
1016 {
1017   //Return current event centrality
1018
1019   if(GetCentrality())
1020   {
1021     if     (fCentralityOpt==100) return (Int_t) GetCentrality()->GetCentralityPercentile(fCentralityClass); // 100 bins max
1022     else if(fCentralityOpt==10)  return GetCentrality()->GetCentralityClass10(fCentralityClass);// 10 bins max
1023     else if(fCentralityOpt==20)  return GetCentrality()->GetCentralityClass5(fCentralityClass); // 20 bins max
1024     else
1025     {
1026       printf("AliCaloTrackReader::GetEventCentrality() - Unknown centrality option %d, use 10, 20 or 100\n",fCentralityOpt);
1027       return -1;
1028     }
1029   }
1030   else return -1;
1031
1032 }
1033
1034 //_____________________________________________________
1035 Double_t AliCaloTrackReader::GetEventPlaneAngle() const
1036 {
1037   //Return current event centrality
1038
1039   if(GetEventPlane())
1040   {
1041     Float_t ep =  GetEventPlane()->GetEventplane(GetEventPlaneMethod(), GetInputEvent());
1042
1043     if(GetEventPlaneMethod()=="Q" && (ep < 0 || ep > TMath::Pi()))
1044     {
1045       if(fDebug > 0 ) printf("AliCaloTrackReader::GetEventPlaneAngle() -  Bad EP for <Q> method : %f\n",ep);
1046       return -1000;
1047     }
1048     else if(GetEventPlaneMethod().Contains("V0")  )
1049     {
1050       if((ep > TMath::Pi()/2 || ep < -TMath::Pi()/2))
1051       {
1052         if(fDebug > 0 ) printf("AliCaloTrackReader::GetEventPlaneAngle() -  Bad EP for <%s> method : %f\n",GetEventPlaneMethod().Data(), ep);
1053         return -1000;
1054       }
1055
1056       ep+=TMath::Pi()/2; // put same range as for <Q> method
1057
1058     }
1059
1060     //printf("AliCaloTrackReader::GetEventPlaneAngle() = %f\n",ep);
1061     if(fDebug > 0 )
1062     {
1063       if     (ep > TMath::Pi()) printf("AliCaloTrackReader::GetEventPlaneAngle() - Too large angle = %f\n",ep);
1064       else if(ep < 0          ) printf("AliCaloTrackReader::GetEventPlaneAngle() - Negative angle = %f\n" ,ep);
1065     }
1066
1067     return ep;
1068   }
1069   else
1070   {
1071     if(fDataType!=kMC && fDebug > 0) printf("AliCaloTrackReader::GetEventPlaneAngle() -  No EP pointer\n");
1072     return -1000;
1073   }
1074
1075 }
1076
1077 //__________________________________________________________
1078 void AliCaloTrackReader::GetVertex(Double_t vertex[3]) const
1079 {
1080   //Return vertex position to be used for single event analysis
1081   vertex[0]=fVertex[0][0];
1082   vertex[1]=fVertex[0][1];
1083   vertex[2]=fVertex[0][2];
1084 }
1085
1086 //____________________________________________________________
1087 void AliCaloTrackReader::GetVertex(Double_t vertex[3],
1088                                    const Int_t evtIndex) const
1089 {
1090   //Return vertex position for mixed event, recover the vertex in a particular event.
1091
1092   vertex[0]=fVertex[evtIndex][0];  vertex[1]=fVertex[evtIndex][1];  vertex[2]=fVertex[evtIndex][2];
1093
1094 }
1095
1096 //________________________________________
1097 void AliCaloTrackReader::FillVertexArray()
1098 {
1099
1100   //Fill data member with vertex
1101   //In case of Mixed event, multiple vertices
1102
1103   //Delete previous vertex
1104   if(fVertex)
1105   {
1106     for (Int_t i = 0; i < fNMixedEvent; i++)
1107     {
1108       delete [] fVertex[i] ;
1109     }
1110     delete [] fVertex ;
1111   }
1112
1113   fVertex = new Double_t*[fNMixedEvent] ;
1114   for (Int_t i = 0; i < fNMixedEvent; i++)
1115   {
1116     fVertex[i] = new Double_t[3] ;
1117     fVertex[i][0] = 0.0 ;
1118     fVertex[i][1] = 0.0 ;
1119     fVertex[i][2] = 0.0 ;
1120   }
1121
1122   if (!fMixedEvent)
1123   { //Single event analysis
1124     if(fDataType!=kMC)
1125     {
1126
1127       if(fInputEvent->GetPrimaryVertex())
1128       {
1129         fInputEvent->GetPrimaryVertex()->GetXYZ(fVertex[0]);
1130       }
1131       else
1132       {
1133         printf("AliCaloTrackReader::FillVertexArray() - NULL primary vertex\n");
1134         fVertex[0][0]=0.;   fVertex[0][1]=0.;   fVertex[0][2]=0.;
1135       }//Primary vertex pointer do not exist
1136
1137     } else
1138     {//MC read event
1139       fVertex[0][0]=0.;   fVertex[0][1]=0.;   fVertex[0][2]=0.;
1140     }
1141
1142     if(fDebug > 1)
1143       printf("AliCaloTrackReader::FillVertexArray() - Single Event Vertex : %f,%f,%f\n",fVertex[0][0],fVertex[0][1],fVertex[0][2]);
1144
1145   } else
1146   { // MultiEvent analysis
1147     for (Int_t iev = 0; iev < fNMixedEvent; iev++)
1148     {
1149       if (fMixedEvent->GetVertexOfEvent(iev))
1150         fMixedEvent->GetVertexOfEvent(iev)->GetXYZ(fVertex[iev]);
1151       else
1152       { // no vertex found !!!!
1153         AliWarning("No vertex found");
1154       }
1155
1156       if(fDebug > 1)
1157         printf("AliCaloTrackReader::FillVertexArray() - Multi Event %d Vertex : %f,%f,%f\n",iev,fVertex[iev][0],fVertex[iev][1],fVertex[iev][2]);
1158
1159     }
1160   }
1161
1162 }
1163
1164 //_____________________________________
1165 void AliCaloTrackReader::FillInputCTS()
1166 {
1167   //Return array with Central Tracking System (CTS) tracks
1168
1169   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputCTS()\n");
1170
1171   Double_t pTrack[3] = {0,0,0};
1172
1173   Int_t nTracks = fInputEvent->GetNumberOfTracks() ;
1174   fTrackMult    = 0;
1175   Int_t nstatus = 0;
1176   Double_t bz   = GetInputEvent()->GetMagneticField();
1177
1178   for(Int_t i = 0; i < 19; i++)
1179   {
1180     fTrackBCEvent   [i] = 0;
1181     fTrackBCEventCut[i] = 0;
1182   }
1183
1184   Bool_t   bc0  = kFALSE;
1185   if(fRecalculateVertexBC) fVertexBC=AliVTrack::kTOFBCNA;
1186
1187   for (Int_t itrack =  0; itrack <  nTracks; itrack++)
1188   {////////////// track loop
1189     AliVTrack * track = (AliVTrack*)fInputEvent->GetTrack(itrack) ; // retrieve track from esd
1190
1191     //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
1192     ULong_t status = track->GetStatus();
1193
1194     if (fTrackStatus && !((status & fTrackStatus) == fTrackStatus))
1195       continue ;
1196
1197     nstatus++;
1198
1199     Float_t dcaTPC =-999;
1200
1201     if     (fDataType==kESD)
1202     {
1203       AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (track);
1204
1205       if(esdTrack)
1206       {
1207         if(fESDtrackCuts->AcceptTrack(esdTrack))
1208         {
1209           track->GetPxPyPz(pTrack) ;
1210
1211           if(fConstrainTrack)
1212           {
1213             if(esdTrack->GetConstrainedParam())
1214             {
1215               const AliExternalTrackParam* constrainParam = esdTrack->GetConstrainedParam();
1216               esdTrack->Set(constrainParam->GetX(),constrainParam->GetAlpha(),constrainParam->GetParameter(),constrainParam->GetCovariance());
1217               esdTrack->GetConstrainedPxPyPz(pTrack);
1218             }
1219             else continue;
1220
1221           } // use constrained tracks
1222
1223           if(fSelectSPDHitTracks)
1224           {//Not much sense to use with TPC only or Hybrid tracks
1225             if(!esdTrack->HasPointOnITSLayer(0) && !esdTrack->HasPointOnITSLayer(1)) continue ;
1226           }
1227         }
1228         // Complementary track to global : Hybrids (make sure that the previous selection is for Global)
1229         else  if(fESDtrackComplementaryCuts && fESDtrackComplementaryCuts->AcceptTrack(esdTrack))
1230         {
1231           // constrain the track
1232           if(esdTrack->GetConstrainedParam())
1233           {
1234             esdTrack->Set(esdTrack->GetConstrainedParam()->GetX(),esdTrack->GetConstrainedParam()->GetAlpha(),esdTrack->GetConstrainedParam()->GetParameter(),esdTrack->GetConstrainedParam()->GetCovariance());
1235
1236             track->GetPxPyPz(pTrack) ;
1237
1238           }
1239           else continue;
1240         }
1241         else continue;
1242       }
1243     } // ESD
1244     else if(fDataType==kAOD)
1245     {
1246       AliAODTrack *aodtrack = dynamic_cast <AliAODTrack*>(track);
1247
1248       if(aodtrack)
1249       {
1250        if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputCTS():AOD track type: %d (primary %d), hybrid? %d \n",
1251                               aodtrack->GetType(),AliAODTrack::kPrimary,
1252                               aodtrack->IsHybridGlobalConstrainedGlobal());
1253
1254
1255         if (fSelectHybridTracks)
1256         {
1257           if (!aodtrack->IsHybridGlobalConstrainedGlobal())       continue ;
1258         }
1259         else
1260         {
1261           if ( aodtrack->TestFilterBit(fTrackFilterMask)==kFALSE) continue ;
1262         }
1263
1264         if(fSelectSPDHitTracks)
1265         {//Not much sense to use with TPC only or Hybrid tracks
1266           if(!aodtrack->HasPointOnITSLayer(0) && !aodtrack->HasPointOnITSLayer(1)) continue ;
1267         }
1268
1269         if (aodtrack->GetType()!= AliAODTrack::kPrimary)          continue ;
1270
1271         if (fDebug > 2 ) printf("AliCaloTrackReader::FillInputCTS(): \t accepted track! \n");
1272
1273         //In case of AODs, TPC tracks cannot be propagated back to primary vertex,
1274         // info stored here
1275         dcaTPC = aodtrack->DCA();
1276
1277         track->GetPxPyPz(pTrack) ;
1278
1279       } // aod track exists
1280       else continue ;
1281
1282     } // AOD
1283
1284     TLorentzVector momentum(pTrack[0],pTrack[1],pTrack[2],0);
1285
1286     Bool_t okTOF  = ( (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ) ;
1287     Double_t tof  = -1000;
1288     Int_t trackBC = -1000 ;
1289
1290     if(okTOF)
1291     {
1292       trackBC = track->GetTOFBunchCrossing(bz);
1293       SetTrackEventBC(trackBC+9);
1294
1295       tof = track->GetTOFsignal()*1e-3;
1296     }
1297
1298     if(fUseTrackDCACut)
1299     {
1300       //normal way to get the dca, cut on dca_xy
1301       if(dcaTPC==-999)
1302       {
1303         Double_t dca[2]   = {1e6,1e6};
1304         Double_t covar[3] = {1e6,1e6,1e6};
1305         Bool_t okDCA = track->PropagateToDCA(fInputEvent->GetPrimaryVertex(),bz,100.,dca,covar);
1306         if( okDCA) okDCA = AcceptDCA(momentum.Pt(),dca[0]);
1307         if(!okDCA)
1308         {
1309           //printf("AliCaloTrackReader::FillInputCTS() - Reject track pt %2.2f, dca_xy %2.4f, BC %d\n",momentum.Pt(),dca[0],trackBC);
1310           continue ;
1311         }
1312       }
1313     }// DCA cuts
1314
1315     if(okTOF)
1316     {
1317       //SetTrackEventBCcut(bc);
1318       SetTrackEventBCcut(trackBC+9);
1319
1320       //After selecting tracks with small DCA, pointing to vertex, set vertex BC depeding on tracks BC
1321       if(fRecalculateVertexBC)
1322       {
1323         if     (trackBC !=0 && trackBC != AliVTrack::kTOFBCNA) fVertexBC = trackBC;
1324         else if(trackBC == 0)                                  bc0 = kTRUE;
1325       }
1326
1327       //In any case, the time should to be larger than the fixed window ...
1328       if( fUseTrackTimeCut && (trackBC!=0 || tof < fTrackTimeCutMin  || tof > fTrackTimeCutMax) )
1329       {
1330         //printf("Remove track time %f and bc = %d\n",tof,trackBC);
1331         continue ;
1332       }
1333       //else printf("Accept track time %f and bc = %d\n",tof,trackBC);
1334
1335     }
1336
1337     //Count the tracks in eta < 0.9
1338     //printf("Eta %f cut  %f\n",TMath::Abs(track->Eta()),fTrackMultEtaCut);
1339     if(TMath::Abs(track->Eta())< fTrackMultEtaCut) fTrackMult++;
1340
1341     if(fCTSPtMin > momentum.Pt() || fCTSPtMax < momentum.Pt()) continue ;
1342
1343     if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"CTS")) continue;
1344
1345     if(fDebug > 2 && momentum.Pt() > 0.1)
1346       printf("AliCaloTrackReader::FillInputCTS() - Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
1347              momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
1348
1349     if (fMixedEvent)  track->SetID(itrack);
1350
1351     fCTSTracks->Add(track);
1352
1353   }// track loop
1354
1355   if(fVertexBC ==0 || fVertexBC == AliVTrack::kTOFBCNA)
1356   {
1357     if( bc0 ) fVertexBC = 0 ;
1358     else      fVertexBC = AliVTrack::kTOFBCNA ;
1359   }
1360
1361   if(fDebug > 1)
1362     printf("AliCaloTrackReader::FillInputCTS()   - aod entries %d, input tracks %d, pass status %d, multipliticy %d\n", fCTSTracks->GetEntriesFast(), nTracks, nstatus, fTrackMult);//fCTSTracksNormalInputEntries);
1363
1364 }
1365
1366 //__________________________________________________________________
1367 void AliCaloTrackReader::FillInputEMCALAlgorithm(AliVCluster * clus,
1368                                                  const Int_t iclus)
1369 {
1370   //Fill the EMCAL data in the array, do it
1371
1372   Int_t vindex = 0 ;
1373   if (fMixedEvent)
1374     vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
1375
1376   //Reject clusters with bad channels, close to borders and exotic;
1377   if(!GetCaloUtils()->GetEMCALRecoUtils()->IsGoodCluster(clus,GetCaloUtils()->GetEMCALGeometry(),GetEMCALCells(),fInputEvent->GetBunchCrossNumber())) return;
1378
1379   //Mask all cells in collumns facing ALICE thick material if requested
1380   if(GetCaloUtils()->GetNMaskCellColumns())
1381   {
1382     Int_t absId   = -1;
1383     Int_t iSupMod = -1;
1384     Int_t iphi    = -1;
1385     Int_t ieta    = -1;
1386     Bool_t shared = kFALSE;
1387     GetCaloUtils()->GetEMCALRecoUtils()->GetMaxEnergyCell(GetCaloUtils()->GetEMCALGeometry(), GetEMCALCells(),clus,absId,iSupMod,ieta,iphi,shared);
1388     if(GetCaloUtils()->MaskFrameCluster(iSupMod, ieta)) return;
1389   }
1390
1391   if(fSelectEmbeddedClusters)
1392   {
1393     if(clus->GetNLabels()==0 || clus->GetLabel() < 0) return;
1394     //else printf("Embedded cluster,  %d, n label %d label %d  \n",iclus,clus->GetNLabels(),clus->GetLabel());
1395   }
1396
1397   //Float_t pos[3];
1398   //clus->GetPosition(pos);
1399   //printf("Before Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
1400
1401   if(fRecalculateClusters)
1402   {
1403     //Recalibrate the cluster energy
1404     if(GetCaloUtils()->IsRecalibrationOn())
1405     {
1406       Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, GetEMCALCells());
1407
1408       clus->SetE(energy);
1409       //printf("Recalibrated Energy %f\n",clus->E());
1410
1411       GetCaloUtils()->RecalculateClusterShowerShapeParameters(GetEMCALCells(),clus);
1412       GetCaloUtils()->RecalculateClusterPID(clus);
1413
1414     } // recalculate E
1415
1416     //Recalculate distance to bad channels, if new list of bad channels provided
1417     GetCaloUtils()->RecalculateClusterDistanceToBadChannel(GetEMCALCells(),clus);
1418
1419     //Recalculate cluster position
1420     if(GetCaloUtils()->IsRecalculationOfClusterPositionOn())
1421     {
1422       GetCaloUtils()->RecalculateClusterPosition(GetEMCALCells(),clus);
1423       //clus->GetPosition(pos);
1424       //printf("After  Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
1425     }
1426
1427     // Recalculate TOF
1428     if(GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn())
1429     {
1430       Double_t tof      = clus->GetTOF();
1431       Float_t  frac     =-1;
1432       Int_t    absIdMax = GetCaloUtils()->GetMaxEnergyCell(fEMCALCells, clus,frac);
1433
1434       if(fDataType==AliCaloTrackReader::kESD)
1435       {
1436         tof = fEMCALCells->GetCellTime(absIdMax);
1437       }
1438
1439       GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
1440
1441       clus->SetTOF(tof);
1442
1443     }// Time recalibration
1444   }
1445
1446   //Correct non linearity
1447   if(GetCaloUtils()->IsCorrectionOfClusterEnergyOn())
1448   {
1449     GetCaloUtils()->CorrectClusterEnergy(clus) ;
1450     //printf("Linearity Corrected Energy %f\n",clus->E());
1451
1452     //In case of MC analysis, to match resolution/calibration in real data
1453     Float_t rdmEnergy = GetCaloUtils()->GetEMCALRecoUtils()->SmearClusterEnergy(clus);
1454     // printf("\t Energy %f, smeared %f\n", clus->E(),rdmEnergy);
1455     clus->SetE(rdmEnergy);
1456   }
1457
1458   Double_t tof = clus->GetTOF()*1e9;
1459
1460   Int_t bc = TMath::Nint(tof/50) + 9;
1461   //printf("tof %2.2f, bc+5=%d\n",tof,bc);
1462
1463   SetEMCalEventBC(bc);
1464
1465   if(fEMCALPtMin > clus->E() || fEMCALPtMax < clus->E()) return ;
1466
1467   TLorentzVector momentum ;
1468
1469   clus->GetMomentum(momentum, fVertex[vindex]);
1470
1471   if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"EMCAL")) return ;
1472
1473   SetEMCalEventBCcut(bc);
1474
1475   if(!IsInTimeWindow(tof,clus->E()))
1476   {
1477     fNPileUpClusters++ ;
1478     if(fUseEMCALTimeCut) return ;
1479   }
1480   else
1481     fNNonPileUpClusters++;
1482
1483   if(fDebug > 2 && momentum.E() > 0.1)
1484     printf("AliCaloTrackReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
1485            momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
1486
1487   if (fMixedEvent)
1488     clus->SetID(iclus) ;
1489
1490   //Correct MC label for AODs
1491
1492   if(ReadAODMCParticles())
1493     CorrectMCLabelForAODs(clus);
1494
1495   fEMCALClusters->Add(clus);
1496
1497 }
1498
1499 //_______________________________________
1500 void AliCaloTrackReader::FillInputEMCAL()
1501 {
1502   //Return array with EMCAL clusters in aod format
1503
1504   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputEMCAL()\n");
1505
1506   // First recalibrate cells, time or energy
1507   //  if(GetCaloUtils()->IsRecalibrationOn())
1508   //    GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCells(GetCaloUtils()->GetEMCALGeometry(),
1509   //                                                          GetEMCALCells(),
1510   //                                                          fInputEvent->GetBunchCrossNumber());
1511
1512   fNPileUpClusters    = 0; // Init counter
1513   fNNonPileUpClusters = 0; // Init counter
1514   for(Int_t i = 0; i < 19; i++)
1515   {
1516     fEMCalBCEvent   [i] = 0;
1517     fEMCalBCEventCut[i] = 0;
1518   }
1519
1520   //Loop to select clusters in fiducial cut and fill container with aodClusters
1521   if(fEMCALClustersListName=="")
1522   {
1523     Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
1524     for (Int_t iclus =  0; iclus <  nclusters; iclus++)
1525     {
1526       AliVCluster * clus = 0;
1527       if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
1528       {
1529         if (IsEMCALCluster(clus))
1530         {
1531           FillInputEMCALAlgorithm(clus, iclus);
1532         }//EMCAL cluster
1533       }// cluster exists
1534     }// cluster loop
1535
1536     //Recalculate track matching
1537     GetCaloUtils()->RecalculateClusterTrackMatching(fInputEvent);
1538
1539   }//Get the clusters from the input event
1540   else
1541   {
1542     TClonesArray * clusterList = 0x0;
1543
1544     if      (fInputEvent->FindListObject(fEMCALClustersListName))
1545     {
1546       clusterList = dynamic_cast<TClonesArray*> (fInputEvent->FindListObject(fEMCALClustersListName));
1547     }
1548     else if(fOutputEvent)
1549     {
1550       clusterList = dynamic_cast<TClonesArray*> (fOutputEvent->FindListObject(fEMCALClustersListName));
1551     }
1552
1553     if(!clusterList)
1554     {
1555         printf("AliCaloTrackReader::FillInputEMCAL() - Wrong name of list with clusters?  <%s>\n",fEMCALClustersListName.Data());
1556         return;
1557     }
1558
1559     Int_t nclusters = clusterList->GetEntriesFast();
1560     for (Int_t iclus =  0; iclus <  nclusters; iclus++)
1561     {
1562       AliVCluster * clus = dynamic_cast<AliVCluster*> (clusterList->At(iclus));
1563       //printf("E %f\n",clus->E());
1564       if (clus) FillInputEMCALAlgorithm(clus, iclus);
1565       else printf("AliCaloTrackReader::FillInputEMCAL() - Null cluster in list!\n");
1566     }// cluster loop
1567
1568     // Recalculate the pile-up time, in case long time clusters removed during clusterization
1569     //printf("Input event INIT : Pile-up clusters %d, NO pile-up %d\n",fNPileUpClusters,fNNonPileUpClusters);
1570
1571     fNPileUpClusters    = 0; // Init counter
1572     fNNonPileUpClusters = 0; // Init counter
1573     for(Int_t i = 0; i < 19; i++)
1574     {
1575       fEMCalBCEvent   [i] = 0;
1576       fEMCalBCEventCut[i] = 0;
1577     }
1578
1579     for (Int_t iclus =  0; iclus < fInputEvent->GetNumberOfCaloClusters(); iclus++)
1580     {
1581       AliVCluster * clus = 0;
1582
1583       if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
1584       {
1585         if (IsEMCALCluster(clus))
1586         {
1587
1588           Float_t  frac     =-1;
1589           Int_t    absIdMax = GetCaloUtils()->GetMaxEnergyCell(fEMCALCells, clus,frac);
1590           Double_t tof = clus->GetTOF();
1591           GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
1592           tof*=1e9;
1593
1594           //printf("Input event cluster : AbsIdMax %d, E %2.2f, time %2.2f \n", absIdMax,clus->E(),tof);
1595
1596           //Reject clusters with bad channels, close to borders and exotic;
1597           if(!GetCaloUtils()->GetEMCALRecoUtils()->IsGoodCluster(clus,GetCaloUtils()->GetEMCALGeometry(),GetEMCALCells(),fInputEvent->GetBunchCrossNumber()))  continue;
1598
1599           Int_t bc = TMath::Nint(tof/50) + 9;
1600           SetEMCalEventBC(bc);
1601
1602           if(fEMCALPtMin > clus->E() || fEMCALPtMax < clus->E()) continue ;
1603
1604           TLorentzVector momentum ;
1605
1606           clus->GetMomentum(momentum, fVertex[0]);
1607
1608           if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"EMCAL")) return ;
1609
1610           SetEMCalEventBCcut(bc);
1611
1612           if(!IsInTimeWindow(tof,clus->E()))
1613             fNPileUpClusters++ ;
1614           else
1615             fNNonPileUpClusters++;
1616
1617         }
1618       }
1619     }
1620
1621     //printf("Input event : Pile-up clusters %d, NO pile-up %d\n",fNPileUpClusters,fNNonPileUpClusters);
1622
1623     // Recalculate track matching, not necessary if already done in the reclusterization task.
1624     // in case it was not done ...
1625     GetCaloUtils()->RecalculateClusterTrackMatching(fInputEvent,clusterList);
1626
1627   }
1628
1629   if(fDebug > 1) printf("AliCaloTrackReader::FillInputEMCAL() - aod entries %d, n pile-up clusters %d, n non pile-up %d \n",  fEMCALClusters->GetEntriesFast(),fNPileUpClusters,fNNonPileUpClusters);
1630
1631 }
1632
1633 //______________________________________
1634 void AliCaloTrackReader::FillInputPHOS()
1635 {
1636   //Return array with PHOS clusters in aod format
1637
1638   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputPHOS()\n");
1639
1640   //Loop to select clusters in fiducial cut and fill container with aodClusters
1641   Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
1642   for (Int_t iclus = 0; iclus < nclusters; iclus++)
1643   {
1644     AliVCluster * clus = 0;
1645     if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
1646     {
1647       if (IsPHOSCluster(clus))
1648       {
1649         //Check if the cluster contains any bad channel and if close to calorimeter borders
1650         Int_t vindex = 0 ;
1651         if (fMixedEvent)
1652           vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
1653         if( GetCaloUtils()->ClusterContainsBadChannel("PHOS",clus->GetCellsAbsId(), clus->GetNCells()))
1654           continue;
1655         if(!GetCaloUtils()->CheckCellFiducialRegion(clus, fInputEvent->GetPHOSCells(), fInputEvent, vindex))
1656           continue;
1657
1658         if(fRecalculateClusters)
1659         {
1660           //Recalibrate the cluster energy
1661           if(GetCaloUtils()->IsRecalibrationOn())
1662           {
1663             Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliAODCaloCells*)GetPHOSCells());
1664             clus->SetE(energy);
1665           }
1666         }
1667
1668         TLorentzVector momentum ;
1669
1670         clus->GetMomentum(momentum, fVertex[vindex]);
1671
1672         if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"PHOS")) continue;
1673
1674         if(fPHOSPtMin > momentum.E() || fPHOSPtMax < momentum.E())          continue;
1675
1676         if(fDebug > 2 && momentum.E() > 0.1)
1677           printf("AliCaloTrackReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
1678                  momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
1679
1680
1681         if (fMixedEvent)
1682         {
1683           clus->SetID(iclus) ;
1684         }
1685
1686         if(ReadAODMCParticles())
1687           CorrectMCLabelForAODs(clus);
1688
1689         fPHOSClusters->Add(clus);
1690
1691       }//PHOS cluster
1692     }//cluster exists
1693   }//esd cluster loop
1694
1695   if(fDebug > 1) printf("AliCaloTrackReader::FillInputPHOS()  - aod entries %d\n",  fPHOSClusters->GetEntriesFast());
1696
1697 }
1698
1699 //____________________________________________
1700 void AliCaloTrackReader::FillInputEMCALCells()
1701 {
1702   //Return array with EMCAL cells in aod format
1703
1704   fEMCALCells = fInputEvent->GetEMCALCells();
1705
1706 }
1707
1708 //___________________________________________
1709 void AliCaloTrackReader::FillInputPHOSCells()
1710 {
1711   //Return array with PHOS cells in aod format
1712
1713   fPHOSCells = fInputEvent->GetPHOSCells();
1714
1715 }
1716
1717 //_______________________________________
1718 void AliCaloTrackReader::FillInputVZERO()
1719 {
1720   //Fill VZERO information in data member, add all the channels information.
1721   AliVVZERO* v0 = fInputEvent->GetVZEROData();
1722   //printf("Init V0: ADC (%d,%d), Multiplicity (%d,%d) \n",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]);
1723
1724   if (v0)
1725   {
1726     AliESDVZERO* esdV0 = dynamic_cast<AliESDVZERO*> (v0);
1727     for (Int_t i = 0; i < 32; i++)
1728     {
1729       if(esdV0)
1730       {//Only available in ESDs
1731         fV0ADC[0] += (Int_t)esdV0->GetAdcV0C(i);
1732         fV0ADC[1] += (Int_t)esdV0->GetAdcV0A(i);
1733       }
1734
1735       fV0Mul[0] += (Int_t)v0->GetMultiplicityV0C(i);
1736       fV0Mul[1] += (Int_t)v0->GetMultiplicityV0A(i);
1737     }
1738     if(fDebug > 0)
1739       printf("V0: ADC (%d,%d), Multiplicity (%d,%d) \n",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]);
1740   }
1741   else
1742   {
1743     if(fDebug > 0)
1744       printf("Cannot retrieve V0 ESD! Run w/ null V0 charges\n ");
1745   }
1746 }
1747
1748
1749 //___________________________________________________________________
1750 Bool_t AliCaloTrackReader::IsEMCALCluster(AliVCluster* cluster) const
1751 {
1752   // Check if it is a cluster from EMCAL. For old AODs cluster type has
1753   // different number and need to patch here
1754
1755   if(fDataType==kAOD && fOldAOD)
1756   {
1757     if (cluster->GetType() == 2) return kTRUE;
1758     else                         return kFALSE;
1759   }
1760   else
1761   {
1762     return cluster->IsEMCAL();
1763   }
1764
1765 }
1766
1767 //___________________________________________________________________
1768 Bool_t AliCaloTrackReader::IsPHOSCluster(AliVCluster * cluster) const
1769 {
1770   //Check if it is a cluster from PHOS.For old AODs cluster type has
1771   // different number and need to patch here
1772
1773   if(fDataType==kAOD && fOldAOD)
1774   {
1775     Int_t type = cluster->GetType();
1776     if (type == 0 || type == 1) return kTRUE;
1777     else                        return kFALSE;
1778   }
1779   else
1780   {
1781     return cluster->IsPHOS();
1782   }
1783
1784 }
1785
1786 //________________________________________________
1787 Bool_t AliCaloTrackReader::CheckForPrimaryVertex()
1788 {
1789   //Check if the vertex was well reconstructed, copy from V0Reader of conversion group
1790   //Only for ESDs ...
1791
1792   AliESDEvent * event = dynamic_cast<AliESDEvent*> (fInputEvent);
1793   if(!event) return kTRUE;
1794
1795   if(event->GetPrimaryVertexTracks()->GetNContributors() > 0)
1796   {
1797     return kTRUE;
1798   }
1799
1800   if(event->GetPrimaryVertexTracks()->GetNContributors() < 1)
1801   {
1802     // SPD vertex
1803     if(event->GetPrimaryVertexSPD()->GetNContributors() > 0)
1804     {
1805       //cout<<"spd vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
1806       return kTRUE;
1807
1808     }
1809     if(event->GetPrimaryVertexSPD()->GetNContributors() < 1)
1810     {
1811       //      cout<<"bad vertex type::"<< event->GetPrimaryVertex()->GetName() << endl;
1812       return kFALSE;
1813     }
1814   }
1815
1816   return kFALSE;
1817
1818 }
1819
1820 //____________________________________________________________
1821 void  AliCaloTrackReader::SetTrackCuts(AliESDtrackCuts * cuts)
1822 {
1823   // Set Track cuts
1824
1825   if(fESDtrackCuts) delete fESDtrackCuts ;
1826
1827   fESDtrackCuts = cuts ;
1828
1829 }
1830
1831 //_________________________________________________________________________
1832 void  AliCaloTrackReader::SetTrackComplementaryCuts(AliESDtrackCuts * cuts)
1833 {
1834   // Set Track cuts for complementary tracks (hybrids)
1835
1836   if(fESDtrackComplementaryCuts) delete fESDtrackComplementaryCuts ;
1837
1838   fESDtrackComplementaryCuts = cuts ;
1839
1840 }
1841
1842