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