PWG/CaloTrackCorrBase/AliCaloTrackReader.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //_________________________________________________________________________
  17 // Base class for reading data: MonteCarlo, ESD or AOD, of PHOS EMCAL and
  18 // Central Barrel Tracking detectors (CTS).
  19 // Not all MC particles/tracks/clusters are kept, some kinematical/fiducial restrictions are done.
  20 // Mother class of : AliCaloTrackESDReader: Fills ESD data in 3 TObjArrays (PHOS, EMCAL, CTS)
  21 //                 : AliCaloTrackMCReader : Fills Kinematics data in 3 TObjArrays (PHOS, EMCAL, CTS)
  22 //                 : AliCaloTrackAODReader: Fills AOD data in 3 TObjArrays (PHOS, EMCAL, CTS)
  23 //-- Author: Gustavo Conesa (LNF-INFN)
  24 //////////////////////////////////////////////////////////////////////////////
  25
  26
  27 // --- ROOT system ---
  28 #include <TFile.h>
  29 #include <TGeoManager.h>
  30 #include <TGeoGlobalMagField.h>
  31
  32 // ---- ANALYSIS system ----
  33 #include "AliMCEvent.h"
  34 #include "AliAODMCHeader.h"
  35 #include "AliGenPythiaEventHeader.h"
  36 #include "AliESDEvent.h"
  37 #include "AliAODEvent.h"
  38 #include "AliVTrack.h"
  39 #include "AliVParticle.h"
  40 #include "AliMixedEvent.h"
  41 #include "AliESDtrack.h"
  42 #include "AliESDtrackCuts.h"
  43 #include "AliTriggerAnalysis.h"
  44 #include "AliESDVZERO.h"
  45 #include "AliVCaloCells.h"
  46 #include "AliMagF.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 fCTSPtMin(0),                fEMCALPtMin(0),                  fPHOSPtMin(0),
  67 fCTSPtMax(0),                fEMCALPtMax(0),                  fPHOSPtMax(0),
  68 fAODBranchList(0x0),
  69 fCTSTracks(0x0),             fEMCALClusters(0x0),             fPHOSClusters(0x0),
  70 fEMCALCells(0x0),            fPHOSCells(0x0),
  71 fInputEvent(0x0),            fOutputEvent(0x0),fMC(0x0),
  72 fFillCTS(0),                 fFillEMCAL(0),                   fFillPHOS(0),
  73 fFillEMCALCells(0),          fFillPHOSCells(0),
  74 fRecalculateClusters(kFALSE),fSelectEmbeddedClusters(kFALSE),
  75 fTrackStatus(0),             fTrackFilterMask(0),             fESDtrackCuts(0),
  76 fTrackMult(0),               fTrackMultEtaCut(0.8),
  77 fReadStack(kFALSE),          fReadAODMCParticles(kFALSE),
  78 fDeltaAODFileName(""),       fFiredTriggerClassName(""),      fAnaLED(kFALSE),
  79 fTaskName(""),               fCaloUtils(0x0),
  80 fMixedEvent(NULL),           fNMixedEvent(0),                 fVertex(NULL),
  81 fWriteOutputDeltaAOD(kFALSE),fOldAOD(kFALSE),                 fCaloFilterPatch(kFALSE),
  82 fEMCALClustersListName(""),  fZvtxCut(0.),
  83 fAcceptFastCluster(kFALSE),  fRemoveLEDEvents(kFALSE),
  84 fDoEventSelection(kFALSE),   fDoV0ANDEventSelection(kFALSE),  fUseEventsWithPrimaryVertex(kFALSE),
  85 fTriggerAnalysis (0x0),
  86 fCentralityClass(""),        fCentralityOpt(0),
  87 fEventPlaneMethod(""),       fImportGeometryFromFile(kFALSE), fImportGeometryFilePath("")
  88
  89 {
  90   //Ctor
  91
  92   //Initialize parameters
  93   InitParameters();
  94 }
  95
  96 //_______________________________________
  97 AliCaloTrackReader::~AliCaloTrackReader()
  98 {
  99   //Dtor
 100
 101   delete fFiducialCut ;
 102
 103   if(fAODBranchList){
 104     fAODBranchList->Delete();
 105     delete fAODBranchList ;
 106   }
 107
 108   if(fCTSTracks){
 109     if(fDataType!=kMC)fCTSTracks->Clear() ;
 110     else              fCTSTracks->Delete() ;
 111     delete fCTSTracks ;
 112   }
 113
 114   if(fEMCALClusters){
 115     if(fDataType!=kMC)fEMCALClusters->Clear("C") ;
 116     else              fEMCALClusters->Delete() ;
 117     delete fEMCALClusters ;
 118   }
 119
 120   if(fPHOSClusters){
 121     if(fDataType!=kMC)fPHOSClusters->Clear("C") ;
 122     else              fPHOSClusters->Delete() ;
 123     delete fPHOSClusters ;
 124   }
 125
 126   if(fVertex){
 127     for (Int_t i = 0; i < fNMixedEvent; i++) {
 128       delete [] fVertex[i] ;
 129
 130     }
 131     delete [] fVertex ;
 132         }
 133
 134   delete fESDtrackCuts;
 135   delete fTriggerAnalysis;
 136
 137   //  Pointers not owned, done by the analysis frame
 138   //  if(fInputEvent)  delete fInputEvent ;
 139   //  if(fOutputEvent) delete fOutputEvent ;
 140   //  if(fMC)          delete fMC ;
 141   //  Pointer not owned, deleted by maker
 142   //  if (fCaloUtils) delete fCaloUtils ;
 143
 144 }
 145
 146 //________________________________________________
 147 Bool_t AliCaloTrackReader::ComparePtHardAndJetPt()
 148 {
 149   // Check the event, if the requested ptHard is much larger than the jet pT, then there is a problem.
 150   // Only for PYTHIA.
 151   if(!fReadStack) return kTRUE; //Information not filtered to AOD
 152
 153   if(!strcmp(GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader")){
 154     TParticle * jet =  0;
 155     AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetGenEventHeader();
 156     Int_t nTriggerJets =  pygeh->NTriggerJets();
 157     Float_t ptHard = pygeh->GetPtHard();
 158
 159     //if(fDebug > 1) printf("AliMCAnalysisUtils::PythiaEventHeader: Njets: %d, pT Hard %f\n",nTriggerJets, ptHard);
 160     Float_t tmpjet[]={0,0,0,0};
 161     for(Int_t ijet = 0; ijet< nTriggerJets; ijet++){
 162       pygeh->TriggerJet(ijet, tmpjet);
 163       jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
 164       //Compare jet pT and pt Hard
 165       //if(fDebug > 1) printf("AliMCAnalysisUtils:: %d pycell jet pT %f\n",ijet, jet->Pt());
 166       if(jet->Pt() > fPtHardAndJetPtFactor * ptHard) {
 167         printf("AliMCAnalysisUtils::PythiaEventHeader: Njets: %d, pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n",
 168                nTriggerJets, ptHard, jet->Pt(), fPtHardAndJetPtFactor);
 169         return kFALSE;
 170       }
 171     }
 172     if(jet) delete jet;
 173   }
 174
 175   return kTRUE ;
 176
 177 }
 178
 179 //____________________________________________
 180 AliStack* AliCaloTrackReader::GetStack() const
 181 {
 182   //Return pointer to stack
 183   if(fMC)
 184     return fMC->Stack();
 185   else{
 186     if(fDebug > 1) printf("AliCaloTrackReader::GetStack() - Stack is not available\n");
 187     return 0x0 ;
 188   }
 189 }
 190
 191 //__________________________________________________
 192 TString AliCaloTrackReader::GetFiredTriggerClasses()
 193 {
 194   // List of triggered classes in a TString
 195
 196   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (GetInputEvent());
 197   AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (GetInputEvent());
 198
 199   if     (esdevent) return esdevent->GetFiredTriggerClasses();
 200   else if(aodevent) return aodevent->GetFiredTriggerClasses();
 201   else              return ""; // Mixed Event, MC event, does not have this trigger info
 202
 203 }
 204
 205 //______________________________________________
 206 AliHeader* AliCaloTrackReader::GetHeader() const
 207 {
 208   //Return pointer to header
 209   if(fMC)
 210     return fMC->Header();
 211   else{
 212     printf("AliCaloTrackReader::Header is not available\n");
 213     return 0x0 ;
 214   }
 215 }
 216
 217 //______________________________________________________________
 218 AliGenEventHeader* AliCaloTrackReader::GetGenEventHeader() const
 219 {
 220   //Return pointer to Generated event header
 221   if(fMC)
 222     return fMC->GenEventHeader();
 223   else{
 224     printf("AliCaloTrackReader::GenEventHeader is not available\n");
 225     return 0x0 ;
 226   }
 227 }
 228
 229 //____________________________________________________________________
 230 TClonesArray* AliCaloTrackReader::GetAODMCParticles(Int_t input) const
 231 {
 232   //Return list of particles in AOD. Do it for the corresponding input event.
 233
 234   TClonesArray * rv = NULL ;
 235   if(fDataType == kAOD){
 236
 237     if(input == 0){
 238       //Normal input AOD
 239       AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fInputEvent) ;
 240       if(evt)
 241         rv = (TClonesArray*)evt->FindListObject("mcparticles");
 242       else
 243         printf("AliCaloTrackReader::GetAODMCParticles() - wrong AOD input index? %d, or non existing tree? \n",input);
 244
 245     }
 246
 247   } else {
 248     printf("AliCaloTrackReader::GetAODMCParticles() - Input are not AODs\n");
 249   }
 250
 251   return rv ;
 252 }
 253
 254 //___________________________________________________________________
 255 AliAODMCHeader* AliCaloTrackReader::GetAODMCHeader(Int_t input) const
 256 {
 257   //Return MC header in AOD. Do it for the corresponding input event.
 258   AliAODMCHeader *mch = NULL;
 259   if(fDataType == kAOD){
 260     //Normal input AOD
 261     if(input == 0) {
 262       mch = (AliAODMCHeader*)((AliAODEvent*)fInputEvent)->FindListObject("mcheader");
 263     }
 264     //          //Second input AOD
 265     //          else if(input == 1){
 266     //       mch = (AliAODMCHeader*) fSecondInputAODEvent->FindListObject("mcheader");
 267     //  }
 268     else {
 269       printf("AliCaloTrackReader::GetAODMCHeader() - wrong AOD input index, %d\n",input);
 270     }
 271   }
 272   else {
 273     printf("AliCaloTrackReader::GetAODMCHeader() - Input are not AODs\n");
 274   }
 275
 276   return mch;
 277 }
 278
 279 //_____________________________
 280 void AliCaloTrackReader::Init()
 281 {
 282   //Init reader. Method to be called in AliAnaPartCorrMaker
 283
 284   //printf(" AliCaloTrackReader::Init() %p \n",gGeoManager);
 285
 286   if(fReadStack && fReadAODMCParticles){
 287     printf("AliCaloTrackReader::Init() - Cannot access stack and mcparticles at the same time, change them \n");
 288     fReadStack          = kFALSE;
 289     fReadAODMCParticles = kFALSE;
 290   }
 291
 292   // Init geometry, I do not like much to do it like this ...
 293   if(fImportGeometryFromFile && !gGeoManager) {
 294     printf("AliCaloTrackReader::Init() - Import geometry.root file\n");
 295     TGeoManager::Import(Form("%s/geometry.root", fImportGeometryFilePath.Data())) ; // default need file "geometry.root" in local dir!!!!
 296   }
 297
 298 }
 299
 300 //_______________________________________
 301 void AliCaloTrackReader::InitParameters()
 302 {
 303   //Initialize the parameters of the analysis.
 304   fDataType   = kESD ;
 305   fCTSPtMin   = 0.1 ;
 306   fEMCALPtMin = 0.1 ;
 307   fPHOSPtMin  = 0.1 ;
 308   fCTSPtMax   = 1000. ;
 309   fEMCALPtMax = 1000. ;
 310   fPHOSPtMax  = 1000. ;
 311
 312   //Do not filter the detectors input by default.
 313   fFillEMCAL      = kFALSE;
 314   fFillPHOS       = kFALSE;
 315   fFillCTS        = kFALSE;
 316   fFillEMCALCells = kFALSE;
 317   fFillPHOSCells  = kFALSE;
 318
 319   fReadStack             = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
 320   fReadAODMCParticles    = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
 321   fDeltaAODFileName      = "deltaAODPartCorr.root";
 322   fFiredTriggerClassName = "";
 323
 324   fAcceptFastCluster = kTRUE;
 325   fAnaLED            = kFALSE;
 326
 327   //We want tracks fitted in the detectors:
 328   //fTrackStatus=AliESDtrack::kTPCrefit;
 329   //fTrackStatus|=AliESDtrack::kITSrefit;
 330   fTrackStatus     = 0;
 331   fTrackFilterMask = 128; //For AODs, but what is the difference between fTrackStatus and fTrackFilterMask?
 332
 333   fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts(); //initialize with TPC only tracks
 334
 335   fV0ADC[0] = 0;   fV0ADC[1] = 0;
 336   fV0Mul[0] = 0;   fV0Mul[1] = 0;
 337
 338   fZvtxCut   = 10.;
 339
 340   fNMixedEvent = 1;
 341
 342   fPtHardAndJetPtFactor = 7;
 343
 344   //Centrality
 345   fCentralityClass  = "V0M";
 346   fCentralityOpt    = 10;
 347   fCentralityBin[0] = fCentralityBin[1]=-1;
 348
 349   fEventPlaneMethod = "Q";
 350
 351   // Allocate memory (not sure this is the right place)
 352   fCTSTracks       = new TObjArray();
 353   fEMCALClusters   = new TObjArray();
 354   fPHOSClusters    = new TObjArray();
 355   fTriggerAnalysis = new AliTriggerAnalysis;
 356   fAODBranchList   = new TList ;
 357
 358   fImportGeometryFromFile = kFALSE;
 359   fImportGeometryFilePath = ".";
 360
 361 }
 362
 363 //________________________________________________________
 364 void AliCaloTrackReader::Print(const Option_t * opt) const
 365 {
 366
 367   //Print some relevant parameters set for the analysis
 368   if(! opt)
 369     return;
 370
 371   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
 372   printf("Task name      : %s\n", fTaskName.Data()) ;
 373   printf("Data type      : %d\n", fDataType) ;
 374   printf("CTS Min pT     : %2.1f GeV/c\n", fCTSPtMin) ;
 375   printf("EMCAL Min pT   : %2.1f GeV/c\n", fEMCALPtMin) ;
 376   printf("PHOS Min pT    : %2.1f GeV/c\n", fPHOSPtMin) ;
 377   printf("CTS Max pT     : %2.1f GeV/c\n", fCTSPtMax) ;
 378   printf("EMCAL Max pT   : %2.1f GeV/c\n", fEMCALPtMax) ;
 379   printf("PHOS Max pT    : %2.1f GeV/c\n", fPHOSPtMax) ;
 380   printf("Use CTS         =     %d\n", fFillCTS) ;
 381   printf("Use EMCAL       =     %d\n", fFillEMCAL) ;
 382   printf("Use PHOS        =     %d\n", fFillPHOS) ;
 383   printf("Use EMCAL Cells =     %d\n", fFillEMCALCells) ;
 384   printf("Use PHOS  Cells =     %d\n", fFillPHOSCells) ;
 385   printf("Track status    =     %d\n", (Int_t) fTrackStatus) ;
 386   printf("Track filter mask (AODs) =  %d\n", (Int_t) fTrackFilterMask) ;
 387   printf("Track Mult Eta Cut =  %d\n", (Int_t) fTrackMultEtaCut) ;
 388   printf("Write delta AOD =     %d\n", fWriteOutputDeltaAOD) ;
 389   printf("Recalculate Clusters = %d\n", fRecalculateClusters) ;
 390
 391   if(fComparePtHardAndJetPt)
 392           printf("Compare jet pt and pt hard to accept event, factor = %2.2f",fPtHardAndJetPtFactor);
 393
 394   printf("Read Kine from, stack? %d, AOD ? %d \n", fReadStack, fReadAODMCParticles) ;
 395   printf("Delta AOD File Name =     %s\n", fDeltaAODFileName.Data()) ;
 396   printf("Centrality: Class %s, Option %d, Bin [%d,%d] \n", fCentralityClass.Data(),fCentralityOpt,fCentralityBin[0], fCentralityBin[1]) ;
 397
 398   printf("    \n") ;
 399
 400 }
 401
 402 //_________________________________________________________________________
 403 Bool_t AliCaloTrackReader::FillInputEvent(const Int_t iEntry,
 404                                           const char * /*currentFileName*/)
 405 {
 406   //Fill the event counter and input lists that are needed, called by the analysis maker.
 407
 408   fEventNumber = iEntry;
 409   //fCurrentFileName = TString(currentFileName);
 410   if(!fInputEvent) {
 411           if(fDebug >= 0) printf("AliCaloTrackReader::FillInputEvent() - Input event not available, skip event analysis\n");
 412           return kFALSE;
 413   }
 414   //Select events only fired by a certain trigger configuration if it is provided
 415   Int_t eventType = 0;
 416   if(fInputEvent->GetHeader())
 417           eventType = ((AliVHeader*)fInputEvent->GetHeader())->GetEventType();
 418
 419   if (GetFiredTriggerClasses().Contains("FAST")  && !GetFiredTriggerClasses().Contains("ALL") && !fAcceptFastCluster) {
 420     if(fDebug > 0)  printf("AliCaloTrackReader::FillInputEvent - Do not count events from fast cluster, trigger name %s\n",fFiredTriggerClassName.Data());
 421     return kFALSE;
 422   }
 423
 424   //-------------------------------------------------------------------------------------
 425   // Reject event if large clusters with large energy
 426   // Use only for LHC11a data for the moment, and if input is clusterizer V1 or V1+unfolding
 427   // If clusterzer NxN or V2 it does not help
 428   //-------------------------------------------------------------------------------------
 429   if(fRemoveLEDEvents){
 430
 431     //printf("Event %d\n",GetEventNumber());
 432     for (Int_t i = 0; i < fInputEvent->GetNumberOfCaloClusters(); i++)
 433     {
 434       AliVCluster *clus = fInputEvent->GetCaloCluster(i);
 435       if(clus->IsEMCAL()){
 436         if ((clus->E() > 500 && clus->GetNCells() > 200 ) || clus->GetNCells() > 200) {
 437           Int_t absID = clus->GetCellsAbsId()[0];
 438           Int_t sm = GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absID);
 439           if(fDebug > 0) printf("AliCaloTrackReader::FillInputEvent - reject event %d with cluster : E %f, ncells %d, absId(0) %d, SM %d\n",GetEventNumber(),clus->E(),  clus->GetNCells(),absID, sm);
 440           return kFALSE;
 441         }
 442       }
 443     }
 444
 445     // Count number of cells with energy larger than 0.1 in SM3, cut on this number
 446     Int_t ncellsSM3 = 0;
 447     Int_t ncellsSM4 = 0;
 448     for(Int_t icell = 0; icell < fInputEvent->GetEMCALCells()->GetNumberOfCells(); icell++){
 449       Int_t absID = fInputEvent->GetEMCALCells()->GetCellNumber(icell);
 450       Int_t sm = GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absID);
 451       if(fInputEvent->GetEMCALCells()->GetAmplitude(icell) > 0.1 && sm==3) ncellsSM3++;
 452       if(fInputEvent->GetEMCALCells()->GetAmplitude(icell) > 0.1 && sm==4) ncellsSM4++;
 453     }
 454
 455     Int_t ncellcut = 21;
 456     if(fFiredTriggerClassName.Contains("EMC")) ncellcut = 35;
 457
 458     if(ncellsSM3 >= ncellcut || ncellsSM4 >= 100) {
 459       if(fDebug > 0) printf(" AliCaloTrackReader::FillInputEvent() - reject event with ncells in SM3 %d and SM4 %d\n",ncellsSM3, ncellsSM4);
 460       return kFALSE;
 461     }
 462   }// Remove LED events
 463
 464   // Reject pure LED events?
 465   if( fFiredTriggerClassName  !="" && !fAnaLED){
 466     if(eventType!=7)
 467       return kFALSE; //Only physics event, do not use for simulated events!!!
 468     if(fDebug > 0)
 469       printf("AliCaloTrackReader::FillInputEvent() - FiredTriggerClass <%s>, selected class <%s>, compare name %d\n",
 470              GetFiredTriggerClasses().Data(),fFiredTriggerClassName.Data(), GetFiredTriggerClasses().Contains(fFiredTriggerClassName));
 471     if( !GetFiredTriggerClasses().Contains(fFiredTriggerClassName) ) return kFALSE;
 472     else if(fDebug > 0) printf("AliCaloTrackReader::FillInputEvent() - Accepted triggered event\n");
 473   }
 474   else if(fAnaLED){
 475     //    kStartOfRun =       1,    // START_OF_RUN
 476     //    kEndOfRun =         2,    // END_OF_RUN
 477     //    kStartOfRunFiles =  3,    // START_OF_RUN_FILES
 478     //    kEndOfRunFiles =    4,    // END_OF_RUN_FILES
 479     //    kStartOfBurst =     5,    // START_OF_BURST
 480     //    kEndOfBurst =       6,    // END_OF_BURST
 481     //    kPhysicsEvent =     7,    // PHYSICS_EVENT
 482     //    kCalibrationEvent = 8,    // CALIBRATION_EVENT
 483     //    kFormatError =      9,    // EVENT_FORMAT_ERROR
 484     //    kStartOfData =      10,   // START_OF_DATA
 485     //    kEndOfData =        11,   // END_OF_DATA
 486     //    kSystemSoftwareTriggerEvent   = 12, // SYSTEM_SOFTWARE_TRIGGER_EVENT
 487     //    kDetectorSoftwareTriggerEvent = 13  // DETECTOR_SOFTWARE_TRIGGER_EVENT
 488
 489           if(eventType!=7 && fDebug > 1 )printf("AliCaloTrackReader::FillInputEvent() - DO LED, Event Type <%d>, 8 Calibration \n",  eventType);
 490           if(eventType!=8)return kFALSE;
 491   }
 492
 493   //In case of analysis of events with jets, skip those with jet pt > 5 pt hard
 494   if(fComparePtHardAndJetPt && GetStack()) {
 495     if(!ComparePtHardAndJetPt()) return kFALSE ;
 496   }
 497
 498   //Fill Vertex array
 499   FillVertexArray();
 500   //Reject events with Z vertex too large, only for SE analysis, if not, cut on the analysis code
 501   if(!GetMixedEvent() && TMath::Abs(fVertex[0][2]) > fZvtxCut) return kFALSE;
 502
 503   //------------------------------------------------------
 504   //Event rejection depending on vertex, pileup, v0and
 505   //------------------------------------------------------
 506   if(fDoEventSelection){
 507     if(!fCaloFilterPatch){
 508       //Do not analyze events with pileup
 509       Bool_t bPileup = fInputEvent->IsPileupFromSPD(3, 0.8, 3., 2., 5.); //Default values, if not it does not compile
 510       //Bool_t bPileup = event->IsPileupFromSPD();
 511       if(bPileup) return kFALSE;
 512
 513       if(fDoV0ANDEventSelection){
 514         Bool_t bV0AND = kTRUE;
 515         AliESDEvent* esd = dynamic_cast<AliESDEvent*> (fInputEvent);
 516         if(esd)
 517           bV0AND = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0AND);
 518         //else bV0AND = //FIXME FOR AODs
 519         if(!bV0AND) return kFALSE;
 520       }
 521
 522       if(fUseEventsWithPrimaryVertex && !CheckForPrimaryVertex()) return kFALSE;
 523
 524     }//CaloFilter patch
 525     else{
 526       if(fInputEvent->GetNumberOfCaloClusters() > 0) {
 527         AliVCluster * calo = fInputEvent->GetCaloCluster(0);
 528         if(calo->GetNLabels() == 4){
 529           Int_t * selection = calo->GetLabels();
 530           Bool_t bPileup = selection[0];
 531           if(bPileup) return kFALSE;
 532
 533           Bool_t bGoodV = selection[1];
 534           if(fUseEventsWithPrimaryVertex && !bGoodV) return kFALSE;
 535
 536           if(fDoV0ANDEventSelection){
 537             Bool_t bV0AND = selection[2];
 538             if(!bV0AND) return kFALSE;
 539           }
 540
 541           fTrackMult = selection[3];
 542           if(fTrackMult == 0) return kFALSE;
 543         } else {
 544           //First filtered AODs, track multiplicity stored there.
 545           fTrackMult = (Int_t) ((AliAODHeader*)fInputEvent->GetHeader())->GetCentrality();
 546           if(fTrackMult == 0) return kFALSE;
 547         }
 548       }//at least one cluster
 549       else {
 550         //printf("AliCaloTrackReader::FillInputEvent() - No clusters in event\n");
 551         //Remove events with  vertex (0,0,0), bad vertex reconstruction
 552         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;
 553
 554         //First filtered AODs, track multiplicity stored there.
 555         fTrackMult = (Int_t) ((AliAODHeader*)fInputEvent->GetHeader())->GetCentrality();
 556         if(fTrackMult == 0) return kFALSE;
 557       }// no cluster
 558     }// CaloFileter patch
 559   }// Event selection
 560   //------------------------------------------------------
 561
 562   //Check if there is a centrality value, PbPb analysis, and if a centrality bin selection is requested
 563   //If we need a centrality bin, we select only those events in the corresponding bin.
 564   if(GetCentrality() && fCentralityBin[0]>=0 && fCentralityBin[1]>=0 && fCentralityOpt==100){
 565     Int_t cen = GetEventCentrality();
 566     if(cen > fCentralityBin[1] || cen < fCentralityBin[0]) return kFALSE; //reject events out of bin.
 567   }
 568
 569   //Fill the arrays with cluster/tracks/cells data
 570   if(fFillEMCALCells)
 571     FillInputEMCALCells();
 572   if(fFillPHOSCells)
 573     FillInputPHOSCells();
 574
 575   if(fFillCTS){
 576     FillInputCTS();
 577     //Accept events with at least one track
 578     if(fTrackMult == 0 && fDoEventSelection) return kFALSE;
 579   }
 580
 581   if(fFillEMCAL)
 582     FillInputEMCAL();
 583   if(fFillPHOS)
 584     FillInputPHOS();
 585
 586   FillInputVZERO();
 587
 588   return kTRUE ;
 589 }
 590
 591 //___________________________________
 592 void AliCaloTrackReader::ResetLists()
 593 {
 594   //  Reset lists, called by the analysis maker
 595
 596   if(fCTSTracks)       fCTSTracks     -> Clear();
 597   if(fEMCALClusters)   fEMCALClusters -> Clear("C");
 598   if(fPHOSClusters)    fPHOSClusters  -> Clear("C");
 599
 600   fV0ADC[0] = 0;   fV0ADC[1] = 0;
 601   fV0Mul[0] = 0;   fV0Mul[1] = 0;
 602
 603 }
 604
 605 //____________________________________________________________
 606 void AliCaloTrackReader::SetInputEvent(AliVEvent* const input)
 607 {
 608   fInputEvent  = input;
 609   fMixedEvent = dynamic_cast<AliMixedEvent*>(GetInputEvent()) ;
 610   if (fMixedEvent) {
 611     fNMixedEvent = fMixedEvent->GetNumberOfEvents() ;
 612   }
 613
 614   //Delete previous vertex
 615   if(fVertex){
 616     for (Int_t i = 0; i < fNMixedEvent; i++) {
 617       delete [] fVertex[i] ;
 618     }
 619     delete [] fVertex ;
 620   }
 621
 622   fVertex = new Double_t*[fNMixedEvent] ;
 623   for (Int_t i = 0; i < fNMixedEvent; i++) {
 624     fVertex[i] = new Double_t[3] ;
 625     fVertex[i][0] = 0.0 ;
 626     fVertex[i][1] = 0.0 ;
 627     fVertex[i][2] = 0.0 ;
 628   }
 629 }
 630
 631 //__________________________________________________
 632 Int_t AliCaloTrackReader::GetEventCentrality() const
 633 {
 634   //Return current event centrality
 635
 636   if(GetCentrality()){
 637     if(fCentralityOpt==100)     return (Int_t) GetCentrality()->GetCentralityPercentile(fCentralityClass); // 100 bins max
 638     else if(fCentralityOpt==10) return GetCentrality()->GetCentralityClass10(fCentralityClass);// 10 bins max
 639     else if(fCentralityOpt==20) return GetCentrality()->GetCentralityClass5(fCentralityClass); // 20 bins max
 640     else {
 641       printf("AliAnaPartCorrBaseClass::Unknown centrality option %d, use 10, 20 or 100\n",fCentralityOpt);
 642       return 0;
 643     }
 644   }
 645   else return 0;
 646
 647 }
 648
 649 //__________________________________________________________
 650 void AliCaloTrackReader::GetVertex(Double_t vertex[3]) const
 651 {
 652   //Return vertex position to be used for single event analysis
 653   vertex[0]=fVertex[0][0];
 654   vertex[1]=fVertex[0][1];
 655   vertex[2]=fVertex[0][2];
 656 }
 657
 658 //____________________________________________________________
 659 void AliCaloTrackReader::GetVertex(Double_t vertex[3],
 660                                    const Int_t evtIndex) const
 661 {
 662   //Return vertex position for mixed event, recover the vertex in a particular event.
 663
 664   vertex[0]=fVertex[evtIndex][0];  vertex[1]=fVertex[evtIndex][1];  vertex[2]=fVertex[evtIndex][2];
 665
 666 }
 667
 668 //________________________________________
 669 void AliCaloTrackReader::FillVertexArray()
 670 {
 671
 672   //Fill data member with vertex
 673   //In case of Mixed event, multiple vertices
 674
 675   //Delete previous vertex
 676   if(fVertex){
 677     for (Int_t i = 0; i < fNMixedEvent; i++) {
 678       delete [] fVertex[i] ;
 679     }
 680     delete [] fVertex ;
 681   }
 682
 683   fVertex = new Double_t*[fNMixedEvent] ;
 684   for (Int_t i = 0; i < fNMixedEvent; i++) {
 685     fVertex[i] = new Double_t[3] ;
 686     fVertex[i][0] = 0.0 ;
 687     fVertex[i][1] = 0.0 ;
 688     fVertex[i][2] = 0.0 ;
 689   }
 690
 691   if (!fMixedEvent) { //Single event analysis
 692     if(fDataType!=kMC){
 693
 694       if(fInputEvent->GetPrimaryVertex()){
 695         fInputEvent->GetPrimaryVertex()->GetXYZ(fVertex[0]);
 696       }
 697       else {
 698         printf("AliCaloTrackReader::FillVertexArray() - NULL primary vertex\n");
 699         fVertex[0][0]=0.;   fVertex[0][1]=0.;   fVertex[0][2]=0.;
 700       }//Primary vertex pointer do not exist
 701
 702     } else {//MC read event
 703       fVertex[0][0]=0.;   fVertex[0][1]=0.;   fVertex[0][2]=0.;
 704     }
 705
 706     if(fDebug > 1)
 707       printf("AliCaloTrackReader::FillVertexArray() - Single Event Vertex : %f,%f,%f\n",fVertex[0][0],fVertex[0][1],fVertex[0][2]);
 708
 709   } else { // MultiEvent analysis
 710     for (Int_t iev = 0; iev < fNMixedEvent; iev++) {
 711       if (fMixedEvent->GetVertexOfEvent(iev))
 712         fMixedEvent->GetVertexOfEvent(iev)->GetXYZ(fVertex[iev]);
 713       else { // no vertex found !!!!
 714         AliWarning("No vertex found");
 715       }
 716
 717       if(fDebug > 1)
 718         printf("AliCaloTrackReader::FillVertexArray() - Multi Event %d Vertex : %f,%f,%f\n",iev,fVertex[iev][0],fVertex[iev][1],fVertex[iev][2]);
 719
 720     }
 721   }
 722
 723 }
 724
 725 //_____________________________________
 726 void AliCaloTrackReader::FillInputCTS()
 727 {
 728   //Return array with Central Tracking System (CTS) tracks
 729
 730   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputCTS()\n");
 731
 732   Int_t nTracks   = fInputEvent->GetNumberOfTracks() ;
 733   Double_t p[3];
 734   fTrackMult = 0;
 735   Int_t nstatus = 0;
 736   for (Int_t itrack =  0; itrack <  nTracks; itrack++) {////////////// track loop
 737     AliVTrack * track = (AliVTrack*)fInputEvent->GetTrack(itrack) ; // retrieve track from esd
 738
 739     //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
 740     if (fTrackStatus && !((track->GetStatus() & fTrackStatus) == fTrackStatus))
 741       continue ;
 742
 743     nstatus++;
 744
 745     if     (fDataType==kESD && !fESDtrackCuts->AcceptTrack((AliESDtrack*)track))
 746     {
 747       continue;
 748     }
 749     else if(fDataType==kAOD)
 750     {
 751       AliAODTrack *aodtrack = dynamic_cast <AliAODTrack*>(track);
 752       if(aodtrack){
 753         if(fDebug > 2 )
 754           printf("AliCaloTrackReader::FillInputCTS():AOD track type: %c \n", aodtrack->GetType());
 755         if (fDataType!=kMC && aodtrack->TestFilterBit(fTrackFilterMask)==kFALSE) continue;
 756         if(aodtrack->GetType()!=AliAODTrack::kPrimary) continue;
 757       }
 758     }
 759
 760     //Count the tracks in eta < 0.9
 761     //printf("Eta %f cut  %f\n",TMath::Abs(track->Eta()),fTrackMultEtaCut);
 762     if(TMath::Abs(track->Eta())< fTrackMultEtaCut) fTrackMult++;
 763
 764     track->GetPxPyPz(p) ;
 765     TLorentzVector momentum(p[0],p[1],p[2],0);
 766
 767     if(fCTSPtMin < momentum.Pt() && fCTSPtMax > momentum.Pt()){
 768
 769       if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"CTS"))
 770         continue;
 771
 772       if(fDebug > 2 && momentum.Pt() > 0.1)
 773         printf("AliCaloTrackReader::FillInputCTS() - Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
 774                momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 775
 776       if (fMixedEvent) {
 777         track->SetID(itrack);
 778       }
 779
 780       fCTSTracks->Add(track);
 781
 782     }//Pt and Fiducial cut passed.
 783   }// track loop
 784
 785   //fCTSTracksNormalInputEntries = fCTSTracks->GetEntriesFast();
 786   if(fDebug > 1)
 787     printf("AliCaloTrackReader::FillInputCTS()   - aod entries %d, input tracks %d, pass status %d, multipliticy %d\n", fCTSTracks->GetEntriesFast(), nTracks, nstatus, fTrackMult);//fCTSTracksNormalInputEntries);
 788
 789 }
 790
 791 //__________________________________________________________________
 792 void AliCaloTrackReader::FillInputEMCALAlgorithm(AliVCluster * clus,
 793                                                  const Int_t iclus)
 794 {
 795   //Fill the EMCAL data in the array, do it
 796
 797   Int_t vindex = 0 ;
 798   if (fMixedEvent)
 799     vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
 800
 801   //Reject clusters with bad channels, close to borders and exotic;
 802   if(!GetCaloUtils()->GetEMCALRecoUtils()->IsGoodCluster(clus,GetCaloUtils()->GetEMCALGeometry(),GetEMCALCells(),fInputEvent->GetBunchCrossNumber())) return;
 803
 804   //Mask all cells in collumns facing ALICE thick material if requested
 805   if(GetCaloUtils()->GetNMaskCellColumns()){
 806     Int_t absId   = -1;
 807     Int_t iSupMod = -1;
 808     Int_t iphi    = -1;
 809     Int_t ieta    = -1;
 810     Bool_t shared = kFALSE;
 811     GetCaloUtils()->GetEMCALRecoUtils()->GetMaxEnergyCell(GetCaloUtils()->GetEMCALGeometry(), GetEMCALCells(),clus,absId,iSupMod,ieta,iphi,shared);
 812     if(GetCaloUtils()->MaskFrameCluster(iSupMod, ieta)) return;
 813   }
 814
 815   if(fSelectEmbeddedClusters){
 816     if(clus->GetNLabels()==0 || clus->GetLabel() < 0) return;
 817     //else printf("Embedded cluster,  %d, n label %d label %d  \n",iclus,clus->GetNLabels(),clus->GetLabel());
 818   }
 819
 820   //Float_t pos[3];
 821   //clus->GetPosition(pos);
 822   //printf("Before Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
 823
 824   if(fRecalculateClusters){
 825     //Recalibrate the cluster energy
 826     if(GetCaloUtils()->IsRecalibrationOn()) {
 827
 828       Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, GetEMCALCells());
 829
 830       clus->SetE(energy);
 831       //printf("Recalibrated Energy %f\n",clus->E());
 832
 833       GetCaloUtils()->RecalculateClusterShowerShapeParameters(GetEMCALCells(),clus);
 834       GetCaloUtils()->RecalculateClusterPID(clus);
 835
 836     } // recalculate E
 837
 838     //Recalculate distance to bad channels, if new list of bad channels provided
 839     GetCaloUtils()->RecalculateClusterDistanceToBadChannel(GetEMCALCells(),clus);
 840
 841     //Recalculate cluster position
 842     if(GetCaloUtils()->IsRecalculationOfClusterPositionOn()){
 843       GetCaloUtils()->RecalculateClusterPosition(GetEMCALCells(),clus);
 844       //clus->GetPosition(pos);
 845       //printf("After  Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
 846     }
 847   }
 848
 849   // Recalculate TOF
 850   if(GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn()) {
 851
 852     Double_t tof      = clus->GetTOF();
 853     Float_t  frac     =-1;
 854     Int_t    absIdMax = GetCaloUtils()->GetMaxEnergyCell(fEMCALCells, clus,frac);
 855
 856     if(fDataType==AliCaloTrackReader::kESD){
 857       tof = fEMCALCells->GetCellTime(absIdMax);
 858     }
 859
 860     GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
 861
 862     clus->SetTOF(tof);
 863
 864   }// Time recalibration
 865
 866   //Correct non linearity
 867   if(GetCaloUtils()->IsCorrectionOfClusterEnergyOn()){
 868     GetCaloUtils()->CorrectClusterEnergy(clus) ;
 869     //printf("Linearity Corrected Energy %f\n",clus->E());
 870
 871     //In case of MC analysis, to match resolution/calibration in real data
 872     Float_t rdmEnergy = GetCaloUtils()->GetEMCALRecoUtils()->SmearClusterEnergy(clus);
 873     // printf("\t Energy %f, smeared %f\n", clus->E(),rdmEnergy);
 874     clus->SetE(rdmEnergy);
 875   }
 876
 877   TLorentzVector momentum ;
 878
 879   clus->GetMomentum(momentum, fVertex[vindex]);
 880
 881   if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"EMCAL")) return;
 882
 883   if(fEMCALPtMin > momentum.E() || fEMCALPtMax < momentum.E())         return;
 884
 885   if(fDebug > 2 && momentum.E() > 0.1)
 886     printf("AliCaloTrackReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
 887            momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 888
 889   if (fMixedEvent)
 890     clus->SetID(iclus) ;
 891
 892   fEMCALClusters->Add(clus);
 893
 894 }
 895
 896 //_______________________________________
 897 void AliCaloTrackReader::FillInputEMCAL()
 898 {
 899   //Return array with EMCAL clusters in aod format
 900
 901   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputEMCAL()\n");
 902
 903   // First recalibrate cells, time or energy
 904   //  if(GetCaloUtils()->IsRecalibrationOn())
 905   //    GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCells(GetCaloUtils()->GetEMCALGeometry(),
 906   //                                                          GetEMCALCells(),
 907   //                                                          fInputEvent->GetBunchCrossNumber());
 908
 909   //Loop to select clusters in fiducial cut and fill container with aodClusters
 910   if(fEMCALClustersListName==""){
 911     Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
 912     for (Int_t iclus =  0; iclus <  nclusters; iclus++) {
 913       AliVCluster * clus = 0;
 914       if ( (clus = fInputEvent->GetCaloCluster(iclus)) ) {
 915         if (IsEMCALCluster(clus)){
 916           FillInputEMCALAlgorithm(clus, iclus);
 917         }//EMCAL cluster
 918       }// cluster exists
 919     }// cluster loop
 920
 921     //Recalculate track matching
 922     GetCaloUtils()->RecalculateClusterTrackMatching(fInputEvent);
 923
 924   }//Get the clusters from the input event
 925   else {
 926     TClonesArray * clusterList = 0x0;
 927     if(fOutputEvent) clusterList = dynamic_cast<TClonesArray*> (fOutputEvent->FindListObject(fEMCALClustersListName));
 928     if(!clusterList){
 929       //printf("AliCaloTrackReader::FillInputEMCAL() - Wrong name of list with clusters? Try input event <%s>\n",fEMCALClustersListName.Data());
 930       //List not in output event, try input event
 931       clusterList = dynamic_cast<TClonesArray*> (fInputEvent->FindListObject(fEMCALClustersListName));
 932       if(!clusterList){
 933         printf("AliCaloTrackReader::FillInputEMCAL() - Wrong name of list with clusters?  <%s>\n",fEMCALClustersListName.Data());
 934         return;
 935       }
 936     }
 937
 938     Int_t nclusters = clusterList->GetEntriesFast();
 939     for (Int_t iclus =  0; iclus <  nclusters; iclus++) {
 940       AliVCluster * clus = dynamic_cast<AliVCluster*> (clusterList->At(iclus));
 941       //printf("E %f\n",clus->E());
 942       if (clus) FillInputEMCALAlgorithm(clus, iclus);
 943       else printf("AliCaloTrackReader::FillInputEMCAL() - Null cluster in list!\n");
 944
 945     }// cluster loop
 946
 947     // Recalculate track matching, not necessary if already done in the reclusterization task.
 948     // in case it was not done ...
 949     GetCaloUtils()->RecalculateClusterTrackMatching(fInputEvent,clusterList);
 950
 951   }
 952
 953   if(fDebug > 1) printf("AliCaloTrackReader::FillInputEMCAL() - aod entries %d\n",  fEMCALClusters->GetEntriesFast());
 954
 955 }
 956
 957 //______________________________________
 958 void AliCaloTrackReader::FillInputPHOS()
 959 {
 960   //Return array with PHOS clusters in aod format
 961
 962   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputPHOS()\n");
 963
 964   //Loop to select clusters in fiducial cut and fill container with aodClusters
 965   Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
 966   for (Int_t iclus = 0; iclus < nclusters; iclus++) {
 967     AliVCluster * clus = 0;
 968     if ( (clus = fInputEvent->GetCaloCluster(iclus)) ) {
 969       if (IsPHOSCluster(clus)){
 970         //Check if the cluster contains any bad channel and if close to calorimeter borders
 971         Int_t vindex = 0 ;
 972         if (fMixedEvent)
 973           vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
 974         if( GetCaloUtils()->ClusterContainsBadChannel("PHOS",clus->GetCellsAbsId(), clus->GetNCells()))
 975           continue;
 976         if(!GetCaloUtils()->CheckCellFiducialRegion(clus, fInputEvent->GetPHOSCells(), fInputEvent, vindex))
 977           continue;
 978
 979         if(fRecalculateClusters){
 980
 981           //Recalibrate the cluster energy
 982           if(GetCaloUtils()->IsRecalibrationOn()) {
 983             Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliAODCaloCells*)GetPHOSCells());
 984             clus->SetE(energy);
 985           }
 986
 987         }
 988
 989         TLorentzVector momentum ;
 990
 991         clus->GetMomentum(momentum, fVertex[vindex]);
 992
 993         if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"PHOS")) continue;
 994
 995         if(fPHOSPtMin > momentum.E() || fPHOSPtMax < momentum.E())          continue;
 996
 997         if(fDebug > 2 && momentum.E() > 0.1)
 998           printf("AliCaloTrackReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
 999                  momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
1000
1001
1002         if (fMixedEvent) {
1003           clus->SetID(iclus) ;
1004         }
1005
1006         fPHOSClusters->Add(clus);
1007
1008       }//PHOS cluster
1009     }//cluster exists
1010   }//esd cluster loop
1011
1012   if(fDebug > 1) printf("AliCaloTrackReader::FillInputPHOS()  - aod entries %d\n",  fPHOSClusters->GetEntriesFast());
1013
1014 }
1015
1016 //____________________________________________
1017 void AliCaloTrackReader::FillInputEMCALCells()
1018 {
1019   //Return array with EMCAL cells in aod format
1020
1021   fEMCALCells = fInputEvent->GetEMCALCells();
1022
1023 }
1024
1025 //___________________________________________
1026 void AliCaloTrackReader::FillInputPHOSCells()
1027 {
1028   //Return array with PHOS cells in aod format
1029
1030   fPHOSCells = fInputEvent->GetPHOSCells();
1031
1032 }
1033
1034 //_______________________________________
1035 void AliCaloTrackReader::FillInputVZERO()
1036 {
1037   //Fill VZERO information in data member, add all the channels information.
1038   AliVVZERO* v0 = fInputEvent->GetVZEROData();
1039   //printf("Init V0: ADC (%d,%d), Multiplicity (%d,%d) \n",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]);
1040
1041   if (v0)
1042   {
1043     AliESDVZERO* esdV0 = dynamic_cast<AliESDVZERO*> (v0);
1044     for (Int_t i = 0; i < 32; i++)
1045     {
1046       if(esdV0){//Only available in ESDs
1047         fV0ADC[0] += (Int_t)esdV0->GetAdcV0C(i);
1048         fV0ADC[1] += (Int_t)esdV0->GetAdcV0A(i);
1049       }
1050       fV0Mul[0] += (Int_t)v0->GetMultiplicityV0C(i);
1051       fV0Mul[1] += (Int_t)v0->GetMultiplicityV0A(i);
1052     }
1053     if(fDebug > 0)
1054       printf("V0: ADC (%d,%d), Multiplicity (%d,%d) \n",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]);
1055   }
1056   else
1057   {
1058     if(fDebug > 0)
1059       printf("Cannot retrieve V0 ESD! Run w/ null V0 charges\n ");
1060   }
1061 }
1062
1063
1064 //___________________________________________________________________
1065 Bool_t AliCaloTrackReader::IsEMCALCluster(AliVCluster* cluster) const
1066 {
1067   // Check if it is a cluster from EMCAL. For old AODs cluster type has
1068   // different number and need to patch here
1069
1070   if(fDataType==kAOD && fOldAOD)
1071   {
1072     if (cluster->GetType() == 2) return kTRUE;
1073     else                         return kFALSE;
1074   }
1075   else
1076   {
1077     return cluster->IsEMCAL();
1078   }
1079
1080 }
1081
1082 //___________________________________________________________________
1083 Bool_t AliCaloTrackReader::IsPHOSCluster(AliVCluster * cluster) const
1084 {
1085   //Check if it is a cluster from PHOS.For old AODs cluster type has
1086   // different number and need to patch here
1087
1088   if(fDataType==kAOD && fOldAOD)
1089   {
1090     Int_t type = cluster->GetType();
1091     if (type == 0 || type == 1) return kTRUE;
1092     else                        return kFALSE;
1093   }
1094   else
1095   {
1096     return cluster->IsPHOS();
1097   }
1098
1099 }
1100
1101 //________________________________________________
1102 Bool_t AliCaloTrackReader::CheckForPrimaryVertex()
1103 {
1104   //Check if the vertex was well reconstructed, copy from V0Reader of conversion group
1105   //Only for ESDs ...
1106
1107   AliESDEvent * event = dynamic_cast<AliESDEvent*> (fInputEvent);
1108   if(!event) return kTRUE;
1109
1110   if(event->GetPrimaryVertexTracks()->GetNContributors() > 0) {
1111     return kTRUE;
1112   }
1113
1114   if(event->GetPrimaryVertexTracks()->GetNContributors() < 1) {
1115     // SPD vertex
1116     if(event->GetPrimaryVertexSPD()->GetNContributors() > 0) {
1117       //cout<<"spd vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
1118       return kTRUE;
1119
1120     }
1121     if(event->GetPrimaryVertexSPD()->GetNContributors() < 1) {
1122       //      cout<<"bad vertex type::"<< event->GetPrimaryVertex()->GetName() << endl;
1123       return kFALSE;
1124     }
1125   }
1126
1127   return kFALSE;
1128
1129 }
1130
1131 //____________________________________________________________
1132 void  AliCaloTrackReader::SetTrackCuts(AliESDtrackCuts * cuts)
1133 {
1134   // Set Track cuts
1135
1136   if(fESDtrackCuts) delete fESDtrackCuts ;
1137
1138   fESDtrackCuts = cuts ;
1139
1140 }
1141