PWG4/PartCorrBase/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 /* $Id:  $ */
  16
  17 //_________________________________________________________________________
  18 // Base class for reading data: MonteCarlo, ESD or AOD, of PHOS EMCAL and
  19 // Central Barrel Tracking detectors (CTS).
  20 // Not all MC particles/tracks/clusters are kept, some kinematical/fiducial restrictions are done.
  21 // Mother class of : AliCaloTrackESDReader: Fills ESD data in 3 TObjArrays (PHOS, EMCAL, CTS)
  22 //                 : AliCaloTrackMCReader: Fills Kinematics data in 3 TObjArrays (PHOS, EMCAL, CTS)
  23 //                 : AliCaloTrackReader: Fills AOD data in 3 TObjArrays (PHOS, EMCAL, CTS)
  24 //
  25 // This part is commented: Mixing analysis can be done, input AOD with events
  26 // is opened in the AliCaloTrackReader::Init()
  27
  28 //-- Author: Gustavo Conesa (LNF-INFN)
  29 //////////////////////////////////////////////////////////////////////////////
  30
  31
  32 // --- ROOT system ---
  33 #include "TFile.h"
  34
  35 //---- ANALYSIS system ----
  36 #include "AliCaloTrackReader.h"
  37 #include "AliMCEvent.h"
  38 #include "AliAODMCHeader.h"
  39 #include "AliGenPythiaEventHeader.h"
  40 #include "AliVEvent.h"
  41 #include "AliAODEvent.h"
  42 #include "AliVTrack.h"
  43 #include "AliVParticle.h"
  44 #include "AliMixedEvent.h"
  45 #include "AliESDtrack.h"
  46 #include "AliEMCALRecoUtils.h"
  47 #include "AliESDtrackCuts.h"
  48
  49 ClassImp(AliCaloTrackReader)
  50
  51
  52 //____________________________________________________________________________
  53   AliCaloTrackReader::AliCaloTrackReader() :
  54     TObject(), fEventNumber(-1), fCurrentFileName(""),fDataType(0), fDebug(0),
  55     fFiducialCut(0x0), fCheckFidCut(kFALSE), fComparePtHardAndJetPt(kFALSE), fPtHardAndJetPtFactor(7),
  56     fCTSPtMin(0), fEMCALPtMin(0),fPHOSPtMin(0), fAODBranchList(new TList ),
  57     fAODCTS(new TObjArray()), fAODEMCAL(new TObjArray()), fAODPHOS(new TObjArray()),
  58     fEMCALCells(0x0), fPHOSCells(0x0),
  59     fInputEvent(0x0), fOutputEvent(0x0),fMC(0x0),
  60     fFillCTS(0),fFillEMCAL(0),fFillPHOS(0),
  61     fFillEMCALCells(0),fFillPHOSCells(0),
  62 //    fSecondInputAODTree(0x0), fSecondInputAODEvent(0x0),
  63 //    fSecondInputFileName(""),fSecondInputFirstEvent(0),
  64 //    fAODCTSNormalInputEntries(0), fAODEMCALNormalInputEntries(0),
  65 //    fAODPHOSNormalInputEntries(0),
  66     fTrackStatus(0),   fESDtrackCuts(0), fTrackMult(0), fTrackMultEtaCut(0.8),
  67     fReadStack(kFALSE), fReadAODMCParticles(kFALSE),
  68     fDeltaAODFileName("deltaAODPartCorr.root"),fFiredTriggerClassName(""),
  69     fAnaLED(kFALSE),fTaskName(""),fCaloUtils(0x0),
  70     fMixedEvent(NULL), fNMixedEvent(1), fVertex(NULL),
  71     fWriteOutputDeltaAOD(kFALSE),fOldAOD(kFALSE),fCaloFilterPatch(kFALSE),
  72     fEMCALClustersListName("")
  73 {
  74   //Ctor
  75
  76   //Initialize parameters
  77   InitParameters();
  78 }
  79
  80 //_________________________________
  81 AliCaloTrackReader::~AliCaloTrackReader() {
  82   //Dtor
  83
  84   if(fFiducialCut) delete fFiducialCut ;
  85
  86   if(fAODBranchList){
  87     fAODBranchList->Delete();
  88     delete fAODBranchList ;
  89   }
  90
  91   if(fAODCTS){
  92     if(fDataType!=kMC)fAODCTS->Clear() ;
  93     else              fAODCTS->Delete() ;
  94     delete fAODCTS ;
  95   }
  96
  97   if(fAODEMCAL){
  98     if(fDataType!=kMC)fAODEMCAL->Clear("C") ;
  99     else              fAODEMCAL->Delete() ;
 100     delete fAODEMCAL ;
 101   }
 102
 103   if(fAODPHOS){
 104     if(fDataType!=kMC)fAODPHOS->Clear("C") ;
 105     else              fAODPHOS->Delete() ;
 106     delete fAODPHOS ;
 107   }
 108
 109 //  if(fEMCALCells){
 110 //    delete fEMCALCells ;
 111 //  }
 112 //
 113 //  if(fPHOSCells){
 114 //    delete fPHOSCells ;
 115 //  }
 116
 117   if(fVertex){
 118     for (Int_t i = 0; i < fNMixedEvent; i++) {
 119       delete [] fVertex[i] ;
 120
 121     }
 122     delete [] fVertex ;
 123         }
 124
 125   if(fESDtrackCuts)   delete fESDtrackCuts;
 126
 127 //  Pointers not owned, done by the analysis frame
 128 //  if(fInputEvent)  delete fInputEvent ;
 129 //  if(fOutputEvent) delete fOutputEvent ;
 130 //  if(fMC)          delete fMC ;
 131
 132 //  if(fSecondInputAODTree){
 133 //    fSecondInputAODTree->Clear();
 134 //    delete fSecondInputAODTree;
 135 //  }
 136 //
 137 //  if(fSecondInputAODEvent) delete fSecondInputAODEvent ;
 138
 139   //  Pointer not owned, deleted by maker
 140   //if (fCaloUtils) delete fCaloUtils ;
 141
 142 }
 143
 144 //_________________________________________________________________________
 145 Bool_t AliCaloTrackReader::ComparePtHardAndJetPt(){
 146   // Check the event, if the requested ptHard is much larger than the jet pT, then there is a problem.
 147   // Only for PYTHIA.
 148   if(!fReadStack) return kTRUE; //Information not filtered to AOD
 149
 150   if(!strcmp(GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader")){
 151     TParticle * jet =  0;
 152     AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetGenEventHeader();
 153     Int_t nTriggerJets =  pygeh->NTriggerJets();
 154     Float_t ptHard = pygeh->GetPtHard();
 155
 156     //if(fDebug > 1) printf("AliMCAnalysisUtils::PythiaEventHeader: Njets: %d, pT Hard %f\n",nTriggerJets, ptHard);
 157     Float_t tmpjet[]={0,0,0,0};
 158     for(Int_t ijet = 0; ijet< nTriggerJets; ijet++){
 159       pygeh->TriggerJet(ijet, tmpjet);
 160       jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
 161       //Compare jet pT and pt Hard
 162       //if(fDebug > 1) printf("AliMCAnalysisUtils:: %d pycell jet pT %f\n",ijet, jet->Pt());
 163       if(jet->Pt() > fPtHardAndJetPtFactor * ptHard) {
 164         printf("AliMCAnalysisUtils::PythiaEventHeader: Njets: %d, pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n",
 165                nTriggerJets, ptHard, jet->Pt(), fPtHardAndJetPtFactor);
 166         return kFALSE;
 167       }
 168     }
 169     if(jet) delete jet;
 170   }
 171
 172   return kTRUE ;
 173
 174 }
 175
 176 //____________________________________________________________________________
 177 AliStack* AliCaloTrackReader::GetStack() const {
 178   //Return pointer to stack
 179   if(fMC)
 180     return fMC->Stack();
 181   else{
 182     if(fDebug > 1) printf("AliCaloTrackReader::GetStack() - Stack is not available\n");
 183     return 0x0 ;
 184   }
 185 }
 186
 187 //____________________________________________________________________________
 188 AliHeader* AliCaloTrackReader::GetHeader() const {
 189   //Return pointer to header
 190   if(fMC)
 191     return fMC->Header();
 192   else{
 193     printf("AliCaloTrackReader::Header is not available\n");
 194     return 0x0 ;
 195   }
 196 }
 197 //____________________________________________________________________________
 198 AliGenEventHeader* AliCaloTrackReader::GetGenEventHeader() const {
 199   //Return pointer to Generated event header
 200   if(fMC)
 201     return fMC->GenEventHeader();
 202   else{
 203     printf("AliCaloTrackReader::GenEventHeader is not available\n");
 204     return 0x0 ;
 205   }
 206 }
 207
 208 //____________________________________________________________________________
 209 TClonesArray* AliCaloTrackReader::GetAODMCParticles(Int_t input) const {
 210   //Return list of particles in AOD. Do it for the corresponding input event.
 211
 212   TClonesArray * rv = NULL ;
 213   if(fDataType == kAOD){
 214
 215     if(input == 0){
 216       //Normal input AOD
 217       AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fInputEvent) ;
 218       if(evt)
 219         rv = (TClonesArray*)evt->FindListObject("mcparticles");
 220       else
 221         printf("AliCaloTrackReader::GetAODMCParticles() - wrong AOD input index? %d, or non existing tree? \n",input);
 222
 223     } //else if(input == 1 && fSecondInputAODEvent){ //Second input AOD
 224 //      rv = (TClonesArray*) fSecondInputAODEvent->FindListObject("mcparticles");
 225 //    }
 226
 227   } else {
 228       printf("AliCaloTrackReader::GetAODMCParticles() - Input are not AODs\n");
 229   }
 230
 231   return rv ;
 232 }
 233
 234 //____________________________________________________________________________
 235 AliAODMCHeader* AliCaloTrackReader::GetAODMCHeader(Int_t input) const {
 236   //Return MC header in AOD. Do it for the corresponding input event.
 237   AliAODMCHeader *mch = NULL;
 238   if(fDataType == kAOD){
 239     //Normal input AOD
 240     if(input == 0) {
 241       mch = (AliAODMCHeader*)((AliAODEvent*)fInputEvent)->FindListObject("mcheader");
 242     }
 243     //          //Second input AOD
 244     //          else if(input == 1){
 245     //       mch = (AliAODMCHeader*) fSecondInputAODEvent->FindListObject("mcheader");
 246     //  }
 247     else {
 248       printf("AliCaloTrackReader::GetAODMCHeader() - wrong AOD input index, %d\n",input);
 249     }
 250   }
 251   else {
 252     printf("AliCaloTrackReader::GetAODMCHeader() - Input are not AODs\n");
 253   }
 254
 255   return mch;
 256 }
 257
 258 //_______________________________________________________________
 259 void AliCaloTrackReader::Init()
 260 {
 261   //Init reader. Method to be called in AliAnaPartCorrMaker
 262
 263   //Get the file with second input events if the filename is given
 264   //Get the tree and connect the AODEvent. Only with AODs
 265
 266   if(fReadStack && fReadAODMCParticles){
 267     printf("AliCaloTrackReader::Init() - Cannot access stack and mcparticles at the same time, change them \n");
 268     fReadStack = kFALSE;
 269     fReadAODMCParticles = kFALSE;
 270   }
 271
 272 //      if(fSecondInputFileName!=""){
 273 //              if(fDataType == kAOD){
 274 //                      TFile * input2   = new TFile(fSecondInputFileName,"read");
 275 //                      printf("AliCaloTrackReader::Init() - Second input file opened: %s, size %d \n", input2->GetName(), (Int_t) input2->GetSize());
 276 //                      fSecondInputAODTree = (TTree*) input2->Get("aodTree");
 277 //                      if(fSecondInputAODTree) printf("AliCaloTrackReader::Init() - Second input tree opened: %s, entries %d \n",
 278 //                                                                                 fSecondInputAODTree->GetName(), (Int_t) fSecondInputAODTree->GetEntries());
 279 //                      else{
 280 //                       printf("AliCaloTrackReader::Init() - Second input tree not available, STOP \n");
 281 //                       abort();
 282 //                      }
 283 //                      fSecondInputAODEvent = new AliAODEvent;
 284 //                      fSecondInputAODEvent->ReadFromTree(fSecondInputAODTree);
 285 //                      if(fSecondInputFirstEvent >= fSecondInputAODTree->GetEntriesFast()){
 286 //                              printf("AliCaloTrackReader::Init() - Requested first event of second input %d, is larger than number of events %d, STOP\n",
 287 //                                         fSecondInputFirstEvent, (Int_t) fSecondInputAODTree->GetEntriesFast());
 288 //                              abort();
 289 //                      }
 290 //              }
 291 //              else printf("AliCaloTrackReader::Init() - Second input not added, reader is not AOD\n");
 292 //      }
 293 }
 294
 295 //_______________________________________________________________
 296 void AliCaloTrackReader::InitParameters()
 297 {
 298   //Initialize the parameters of the analysis.
 299   fDataType   = kESD ;
 300   fCTSPtMin   = 0.2 ;
 301   fEMCALPtMin = 0.2 ;
 302   fPHOSPtMin  = 0.2 ;
 303
 304   //Do not filter the detectors input by default.
 305   fFillEMCAL      = kFALSE;
 306   fFillPHOS       = kFALSE;
 307   fFillCTS        = kFALSE;
 308   fFillEMCALCells = kFALSE;
 309   fFillPHOSCells  = kFALSE;
 310
 311   //fSecondInputFileName   = "" ;
 312   //fSecondInputFirstEvent = 0 ;
 313   fReadStack             = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
 314   fReadAODMCParticles    = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
 315   fDeltaAODFileName      = "deltaAODPartCorr.root";
 316   fFiredTriggerClassName      = "";
 317
 318   fAnaLED = kFALSE;
 319
 320   //We want tracks fitted in the detectors:
 321   //fTrackStatus=AliESDtrack::kTPCrefit;
 322   //fTrackStatus|=AliESDtrack::kITSrefit;
 323
 324   fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
 325
 326   fV0ADC[0] = 0;   fV0ADC[1] = 0;
 327   fV0Mul[0] = 0;   fV0Mul[1] = 0;
 328
 329 }
 330
 331 //________________________________________________________________
 332 void AliCaloTrackReader::Print(const Option_t * opt) const
 333 {
 334
 335   //Print some relevant parameters set for the analysis
 336   if(! opt)
 337     return;
 338
 339   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
 340   printf("Task name      : %s\n", fTaskName.Data()) ;
 341   printf("Data type      : %d\n", fDataType) ;
 342   printf("CTS Min pT     : %2.1f GeV/c\n", fCTSPtMin) ;
 343   printf("EMCAL Min pT   : %2.1f GeV/c\n", fEMCALPtMin) ;
 344   printf("PHOS Min pT    : %2.1f GeV/c\n", fPHOSPtMin) ;
 345   printf("Use CTS         =     %d\n", fFillCTS) ;
 346   printf("Use EMCAL       =     %d\n", fFillEMCAL) ;
 347   printf("Use PHOS        =     %d\n", fFillPHOS) ;
 348   printf("Use EMCAL Cells =     %d\n", fFillEMCALCells) ;
 349   printf("Use PHOS  Cells =     %d\n", fFillPHOSCells) ;
 350   printf("Track status    =     %d\n", (Int_t) fTrackStatus) ;
 351   printf("Track Mult Eta Cut =  %d\n", (Int_t) fTrackMultEtaCut) ;
 352   printf("Write delta AOD =     %d\n", fWriteOutputDeltaAOD) ;
 353
 354   if(fComparePtHardAndJetPt)
 355           printf("Compare jet pt and pt hard to accept event, factor = %2.2f",fPtHardAndJetPtFactor);
 356
 357 //  if(fSecondInputFileName!="") {
 358 //        printf("Second Input File Name     =     %s\n", fSecondInputFileName.Data()) ;
 359 //        printf("Second Input First Event   =     %d\n", fSecondInputFirstEvent) ;
 360 //  }
 361
 362   printf("Read Kine from, stack? %d, AOD ? %d \n", fReadStack, fReadAODMCParticles) ;
 363   printf("Delta AOD File Name =     %s\n", fDeltaAODFileName.Data()) ;
 364   printf("    \n") ;
 365 }
 366
 367 //___________________________________________________
 368 Bool_t AliCaloTrackReader::FillInputEvent(const Int_t iEntry, const char * currentFileName) {
 369   //Fill the event counter and input lists that are needed, called by the analysis maker.
 370
 371   fEventNumber = iEntry;
 372   fCurrentFileName = TString(currentFileName);
 373   if(!fInputEvent) {
 374           if(fDebug >= 0) printf("AliCaloTrackReader::FillInputEvent() - Input event not available, skip event analysis\n");
 375           return kFALSE;
 376   }
 377   //Select events only fired by a certain trigger configuration if it is provided
 378   Int_t eventType = 0;
 379   if(fInputEvent->GetHeader())
 380           eventType = ((AliVHeader*)fInputEvent->GetHeader())->GetEventType();
 381   if( fFiredTriggerClassName  !="" && !fAnaLED){
 382     if(eventType!=7)
 383       return kFALSE; //Only physics event, do not use for simulated events!!!
 384     if(fDebug > 0)
 385       printf("AliCaloTrackReader::FillInputEvent() - FiredTriggerClass <%s>, selected class <%s>, compare name %d\n",
 386              GetFiredTriggerClasses().Data(),fFiredTriggerClassName.Data(), GetFiredTriggerClasses().Contains(fFiredTriggerClassName));
 387     if( !GetFiredTriggerClasses().Contains(fFiredTriggerClassName) ) return kFALSE;
 388   }
 389   else if(fAnaLED){
 390 //        kStartOfRun =       1,    // START_OF_RUN
 391 //        kEndOfRun =         2,    // END_OF_RUN
 392 //        kStartOfRunFiles =  3,    // START_OF_RUN_FILES
 393 //        kEndOfRunFiles =    4,    // END_OF_RUN_FILES
 394 //        kStartOfBurst =     5,    // START_OF_BURST
 395 //        kEndOfBurst =       6,    // END_OF_BURST
 396 //        kPhysicsEvent =     7,    // PHYSICS_EVENT
 397 //        kCalibrationEvent = 8,    // CALIBRATION_EVENT
 398 //        kFormatError =      9,    // EVENT_FORMAT_ERROR
 399 //        kStartOfData =      10,   // START_OF_DATA
 400 //        kEndOfData =        11,   // END_OF_DATA
 401 //        kSystemSoftwareTriggerEvent   = 12, // SYSTEM_SOFTWARE_TRIGGER_EVENT
 402 //        kDetectorSoftwareTriggerEvent = 13  // DETECTOR_SOFTWARE_TRIGGER_EVENT
 403
 404           if(eventType!=7 && fDebug > 1 )printf("AliCaloTrackReader::FillInputEvent() - DO LED, Event Type <%d>, 8 Calibration \n",  eventType);
 405           if(eventType!=8)return kFALSE;
 406   }
 407
 408   //In case of analysis of events with jets, skip those with jet pt > 5 pt hard
 409   if(fComparePtHardAndJetPt && GetStack()) {
 410     if(!ComparePtHardAndJetPt()) return kFALSE ;
 411   }
 412
 413   //In case of mixing events with other AOD file
 414  // if(fDataType == kAOD && fSecondInputAODTree){
 415 //
 416 //    if(fDebug > 1)
 417 //      printf("AliCaloTrackReader::FillInputEvent() - Get event %d from second input AOD file \n", iEntry+fSecondInputFirstEvent);
 418 //    if(fSecondInputAODTree->GetEntriesFast() <= iEntry+fSecondInputFirstEvent) {
 419 //      if(fSecondInputAODTree->GetEntriesFast() == iEntry+fSecondInputFirstEvent)
 420 //                       printf("AliCaloTrackReader::FillInputEvent() - Skip events from event %d, no more events in second AOD file \n", iEntry);
 421 //      return kFALSE;
 422 //    }
 423 //
 424 //    //Get the Event
 425 //    Int_t nbytes = fSecondInputAODTree->GetEvent(iEntry+fSecondInputFirstEvent);
 426 //    if ( nbytes == 0 ) {//If nothing in AOD
 427 //      printf("AliCaloTrackReader::FillInputEvent() - Nothing in Second AOD input, STOP\n");
 428 //      abort() ;
 429 //    }
 430 //
 431 //  }
 432
 433   //Fill Vertex array
 434
 435   FillVertexArray();
 436
 437   //Fill the arrays with cluster/tracks/cells data
 438    if(fFillEMCALCells)
 439     FillInputEMCALCells();
 440   if(fFillPHOSCells)
 441     FillInputPHOSCells();
 442
 443   if(fFillCTS){
 444     FillInputCTS();
 445     //Accept events with at least one track
 446     if(fTrackMult == 0) return kFALSE;
 447   }
 448
 449   //In case of data produced with calo filter, some information stored in non usual places
 450   if(IsCaloFilterPatchOn()){
 451     fTrackMult = (Int_t) ((AliAODHeader*)fInputEvent->GetHeader())->GetCentrality();
 452     //printf("Track multiplicity %d \n",fTrackMult);
 453     if(fTrackMult == 0) return kFALSE;
 454   }
 455
 456   if(fFillEMCAL)
 457     FillInputEMCAL();
 458   if(fFillPHOS)
 459     FillInputPHOS();
 460
 461   FillInputVZERO();
 462
 463   return kTRUE ;
 464 }
 465
 466 //__________________________________________________
 467 void AliCaloTrackReader::ResetLists() {
 468   //  Reset lists, called by the analysis maker
 469
 470   if(fAODCTS)     fAODCTS     -> Clear();
 471   if(fAODEMCAL)   fAODEMCAL   -> Clear("C");
 472   if(fAODPHOS)    fAODPHOS    -> Clear("C");
 473 //  if(fEMCALCells) fEMCALCells -> Clear("");
 474 //  if(fPHOSCells)  fPHOSCells  -> Clear("");
 475
 476   fV0ADC[0] = 0;   fV0ADC[1] = 0;
 477   fV0Mul[0] = 0;   fV0Mul[1] = 0;
 478
 479 }
 480
 481 //____________________________________________________________________________
 482 void AliCaloTrackReader::SetInputEvent(AliVEvent* const input)
 483 {
 484   fInputEvent  = input;
 485   fMixedEvent = dynamic_cast<AliMixedEvent*>(GetInputEvent()) ;
 486   if (fMixedEvent) {
 487     fNMixedEvent = fMixedEvent->GetNumberOfEvents() ;
 488   }
 489
 490   //Delete previous vertex
 491   if(fVertex){
 492     for (Int_t i = 0; i < fNMixedEvent; i++) {
 493       delete [] fVertex[i] ;
 494     }
 495     delete [] fVertex ;
 496   }
 497
 498   fVertex = new Double_t*[fNMixedEvent] ;
 499   for (Int_t i = 0; i < fNMixedEvent; i++) {
 500     fVertex[i] = new Double_t[3] ;
 501     fVertex[i][0] = 0.0 ;
 502     fVertex[i][1] = 0.0 ;
 503     fVertex[i][2] = 0.0 ;
 504   }
 505 }
 506
 507 //____________________________________________________________________________
 508 void AliCaloTrackReader::GetVertex(Double_t vertex[3]) const {
 509   //Return vertex position to be used for single event analysis
 510   vertex[0]=fVertex[0][0];
 511   vertex[1]=fVertex[0][1];
 512   vertex[2]=fVertex[0][2];
 513 }
 514
 515 //____________________________________________________________________________
 516 void AliCaloTrackReader::GetVertex(Double_t vertex[3], const Int_t evtIndex) const {
 517   //Return vertex position for mixed event, recover the vertex in a particular event.
 518
 519   //Int_t evtIndex = 0; // for single events only one vertex stored in position 0, default value
 520   //if (fMixedEvent && clusterID >=0) {
 521   //  evtIndex=GetMixedEvent()->EventIndexForCaloCluster(clusterID) ;
 522   //}
 523
 524   vertex[0]=fVertex[evtIndex][0];  vertex[1]=fVertex[evtIndex][1];  vertex[2]=fVertex[evtIndex][2];
 525
 526 }
 527 //
 528
 529
 530 //____________________________________________________________________________
 531 void AliCaloTrackReader::FillVertexArray() {
 532
 533   //Fill data member with vertex
 534   //In case of Mixed event, multiple vertices
 535
 536   //Delete previous vertex
 537   if(fVertex){
 538     for (Int_t i = 0; i < fNMixedEvent; i++) {
 539       delete [] fVertex[i] ;
 540     }
 541     delete [] fVertex ;
 542   }
 543
 544   fVertex = new Double_t*[fNMixedEvent] ;
 545   for (Int_t i = 0; i < fNMixedEvent; i++) {
 546     fVertex[i] = new Double_t[3] ;
 547     fVertex[i][0] = 0.0 ;
 548     fVertex[i][1] = 0.0 ;
 549     fVertex[i][2] = 0.0 ;
 550   }
 551
 552   if (!fMixedEvent) { //Single event analysis
 553     if(fDataType!=kMC){
 554
 555       if(fInputEvent->GetPrimaryVertex()){
 556         fInputEvent->GetPrimaryVertex()->GetXYZ(fVertex[0]);
 557       }
 558       else {
 559         printf("AliCaloTrackReader::FillVertexArray() - NULL primary vertex\n");
 560         fVertex[0][0]=0.;   fVertex[0][1]=0.;   fVertex[0][2]=0.;
 561       }//Primary vertex pointer do not exist
 562
 563     } else {//MC read event
 564       fVertex[0][0]=0.;   fVertex[0][1]=0.;   fVertex[0][2]=0.;
 565     }
 566
 567     if(fDebug > 1)
 568       printf("AliCaloTrackReader::FillVertexArray() - Single Event Vertex : %f,%f,%f\n",fVertex[0][0],fVertex[0][1],fVertex[0][2]);
 569
 570   } else { // MultiEvent analysis
 571     for (Int_t iev = 0; iev < fNMixedEvent; iev++) {
 572       if (fMixedEvent->GetVertexOfEvent(iev))
 573         fMixedEvent->GetVertexOfEvent(iev)->GetXYZ(fVertex[iev]);
 574       else { // no vertex found !!!!
 575         AliWarning("No vertex found");
 576       }
 577
 578       if(fDebug > 1)
 579         printf("AliCaloTrackReader::FillVertexArray() - Multi Event %d Vertex : %f,%f,%f\n",iev,fVertex[iev][0],fVertex[iev][1],fVertex[iev][2]);
 580
 581     }
 582   }
 583
 584 }
 585
 586 //____________________________________________________________________________
 587 void AliCaloTrackReader::FillInputCTS() {
 588   //Return array with Central Tracking System (CTS) tracks
 589
 590   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputCTS()\n");
 591
 592   Int_t nTracks   = fInputEvent->GetNumberOfTracks() ;
 593   Double_t p[3];
 594   fTrackMult = 0;
 595   Int_t nstatus = 0;
 596   for (Int_t itrack =  0; itrack <  nTracks; itrack++) {////////////// track loop
 597     AliVTrack * track = (AliVTrack*)fInputEvent->GetTrack(itrack) ; // retrieve track from esd
 598
 599     //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
 600     if (fTrackStatus && !((track->GetStatus() & fTrackStatus) == fTrackStatus))
 601       continue ;
 602
 603     nstatus++;
 604
 605     if(fDataType==kESD && !fESDtrackCuts->AcceptTrack((AliESDtrack*)track)) continue;
 606
 607     // Track filter selection
 608     //if (fTrackFilter) {
 609           //  selectInfo = fTrackFilter->IsSelected(esdTrack);
 610           //  if (!selectInfo && !(esd->GetPrimaryVertex())->UsesTrack(esdTrack->GetID())) continue;
 611    // }
 612
 613     //Count the tracks in eta < 0.9
 614     //printf("Eta %f cut  %f\n",TMath::Abs(track->Eta()),fTrackMultEtaCut);
 615     if(TMath::Abs(track->Eta())< fTrackMultEtaCut) fTrackMult++;
 616
 617     track->GetPxPyPz(p) ;
 618     TLorentzVector momentum(p[0],p[1],p[2],0);
 619
 620     if(fCTSPtMin < momentum.Pt()){
 621
 622       if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"CTS"))
 623         continue;
 624
 625       if(fDebug > 2 && momentum.Pt() > 0.1)
 626         printf("AliCaloTrackReader::FillInputCTS() - Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
 627                momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 628
 629       if (fMixedEvent) {
 630         track->SetID(itrack);
 631       }
 632
 633       fAODCTS->Add(track);
 634
 635     }//Pt and Fiducial cut passed.
 636   }// track loop
 637
 638   //fAODCTSNormalInputEntries = fAODCTS->GetEntriesFast();
 639   if(fDebug > 1)
 640     printf("AliCaloTrackReader::FillInputCTS()   - aod entries %d, input tracks %d, pass status %d, multipliticy %d\n", fAODCTS->GetEntriesFast(), nTracks, nstatus, fTrackMult);//fAODCTSNormalInputEntries);
 641
 642     //  //If second input event available, add the clusters.
 643     //  if(fSecondInputAODTree && fSecondInputAODEvent){
 644     //    nTracks   = fSecondInputAODEvent->GetNumberOfTracks() ;
 645     //    if(fDebug > 1) printf("AliCaloTrackReader::FillInputCTS()   - Add second input tracks, entries %d\n", nTracks) ;
 646     //    for (Int_t itrack =  0; itrack <  nTracks; itrack++) {////////////// track loop
 647     //            AliAODTrack * track = ((AliAODEvent*)fSecondInputAODEvent)->GetTrack(itrack) ; // retrieve track from esd
 648     //
 649     //            //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
 650     //            if (fTrackStatus && !((track->GetStatus() & fTrackStatus) == fTrackStatus)) continue ;
 651     //
 652     //            track->GetPxPyPz(p) ;
 653     //            TLorentzVector momentum(p[0],p[1],p[2],0);
 654     //
 655     //            if(fCTSPtMin < momentum.Pt()){
 656     //
 657     //                    if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"CTS")) continue;
 658     //
 659     //                    if(fDebug > 2 && momentum.Pt() > 0.1) printf("AliCaloTrackReader::FillInputCTS() - Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
 660     //                                                                 momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 661     //
 662     //                    fAODCTS->Add(track);
 663     //
 664     //            }//Pt and Fiducial cut passed.
 665     //    }// track loop
 666     //
 667     //    if(fDebug > 1) printf("AliCaloTrackReader::FillInputCTS()   - aod normal entries %d, after second input %d\n", fAODCTSNormalInputEntries, fAODCTS->GetEntriesFast());
 668     //  }       //second input loop
 669     //
 670 }
 671
 672 //____________________________________________________________________________
 673 void AliCaloTrackReader::FillInputEMCALAlgorithm(AliVCluster * clus, const Int_t iclus) {
 674   //Fill the EMCAL data in the array, do it
 675
 676   Int_t vindex = 0 ;
 677   if (fMixedEvent)
 678     vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
 679
 680   //Check if the cluster contains any bad channel and if close to calorimeter borders
 681   if(GetCaloUtils()->ClusterContainsBadChannel("EMCAL",clus->GetCellsAbsId(), clus->GetNCells()))
 682     return;
 683   if(!GetCaloUtils()->CheckCellFiducialRegion(clus, (AliVCaloCells*)fInputEvent->GetEMCALCells(), fInputEvent, vindex))
 684     return;
 685
 686   TLorentzVector momentum ;
 687
 688   clus->GetMomentum(momentum, fVertex[vindex]);
 689
 690   if(fEMCALPtMin < momentum.Pt()){
 691
 692     if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"EMCAL"))
 693       return;
 694
 695     if(fDebug > 2 && momentum.E() > 0.1)
 696       printf("AliCaloTrackReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
 697              momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 698
 699     //Float_t pos[3];
 700     //clus->GetPosition(pos);
 701     //printf("Before Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
 702
 703     //Recalibrate the cluster energy
 704     if(GetCaloUtils()->IsRecalibrationOn()) {
 705       Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, GetEMCALCells());
 706       clus->SetE(energy);
 707       //printf("Recalibrated Energy %f\n",clus->E());
 708       GetCaloUtils()->RecalculateClusterShowerShapeParameters(GetEMCALCells(),clus);
 709       GetCaloUtils()->RecalculateClusterPID(clus);
 710
 711     }
 712
 713     //Recalculate distance to bad channels, if new list of bad channels provided
 714     GetCaloUtils()->RecalculateClusterDistanceToBadChannel(GetEMCALCells(),clus);
 715
 716     //Recalculate cluster position
 717     if(GetCaloUtils()->IsRecalculationOfClusterPositionOn()){
 718       GetCaloUtils()->RecalculateClusterPosition(GetEMCALCells(),clus);
 719       //clus->GetPosition(pos);
 720       //printf("After  Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
 721     }
 722
 723     //Correct non linearity
 724     if(GetCaloUtils()->IsCorrectionOfClusterEnergyOn()){
 725       GetCaloUtils()->CorrectClusterEnergy(clus) ;
 726       //printf("Linearity Corrected Energy %f\n",clus->E());
 727     }
 728
 729     if (fMixedEvent)
 730       clus->SetID(iclus) ;
 731
 732     fAODEMCAL->Add(clus);
 733   }
 734 }
 735
 736 //____________________________________________________________________________
 737 void AliCaloTrackReader::FillInputEMCAL() {
 738   //Return array with EMCAL clusters in aod format
 739
 740   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputEMCAL()\n");
 741
 742   //Loop to select clusters in fiducial cut and fill container with aodClusters
 743   if(fEMCALClustersListName==""){
 744     Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
 745     for (Int_t iclus =  0; iclus <  nclusters; iclus++) {
 746       AliVCluster * clus = 0;
 747       if ( (clus = fInputEvent->GetCaloCluster(iclus)) ) {
 748         if (IsEMCALCluster(clus)){
 749           FillInputEMCALAlgorithm(clus, iclus);
 750         }//EMCAL cluster
 751       }// cluster exists
 752     }// cluster loop
 753   }//Get the clusters from the input event
 754   else {
 755     TClonesArray * clusterList = dynamic_cast<TClonesArray*> (fOutputEvent->FindListObject(fEMCALClustersListName));
 756     if(!clusterList){
 757       printf("AliCaloTrackReader::FillInputEMCAL() - Wrong name of list with clusters? <%s>\n",fEMCALClustersListName.Data());
 758       return;
 759     }
 760     Int_t nclusters = clusterList->GetEntriesFast();
 761     for (Int_t iclus =  0; iclus <  nclusters; iclus++) {
 762       AliVCluster * clus = dynamic_cast<AliVCluster*> (clusterList->At(iclus));
 763       //printf("E %f\n",clus->E());
 764       FillInputEMCALAlgorithm(clus, iclus);
 765     }// cluster loop
 766   }
 767
 768   //Recalculate track matching
 769   GetCaloUtils()->RecalculateClusterTrackMatching(fInputEvent);
 770
 771   //fAODEMCALNormalInputEntries = fAODEMCAL->GetEntriesFast();
 772   if(fDebug > 1) printf("AliCaloTrackReader::FillInputEMCAL() - aod entries %d\n",  fAODEMCAL->GetEntriesFast());//fAODEMCALNormalInputEntries);
 773
 774     //If second input event available, add the clusters.
 775     //  if(fSecondInputAODTree && fSecondInputAODEvent){
 776     //    GetSecondInputAODVertex(v);
 777     //    nclusters = ((AliAODEvent*)fSecondInputAODEvent)->GetNumberOfCaloClusters();
 778     //    if(fDebug > 1) printf("AliCaloTrackReader::FillInputEMCAL() - Add second input clusters, entries %d\n", nclusters) ;
 779     //          for (Int_t iclus =  0; iclus < nclusters; iclus++) {
 780     //                  AliAODCaloCluster * clus = 0;
 781     //                  if ( (clus = ((AliAODEvent*)fSecondInputAODEvent)->GetCaloCluster(iclus)) ) {
 782     //                          if (clus->IsEMCAL()){
 783     //                                  TLorentzVector momentum ;
 784     //                                  clus->GetMomentum(momentum, v);
 785     //
 786     //                                  if(fEMCALPtMin < momentum.Pt()){
 787     //
 788     //                                          if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"EMCAL")) continue;
 789     //
 790     //                                          if(fDebug > 2 && momentum.E() > 0.1) printf("AliCaloTrackReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
 791     //                                                                                                                                  momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 792     //                                    fAODEMCAL->Add(clus);
 793     //                                  }//Pt and Fiducial cut passed.
 794     //                          }//EMCAL cluster
 795     //                  }// cluster exists
 796     //          }// cluster loop
 797     //
 798     //    if(fDebug > 1) printf("AliCaloTrackReader::FillInputEMCAL() - aod normal entries %d, after second input %d\n", fAODEMCALNormalInputEntries, fAODEMCAL->GetEntriesFast());
 799     //
 800     //  } //second input loop
 801 }
 802
 803   //____________________________________________________________________________
 804 void AliCaloTrackReader::FillInputPHOS() {
 805   //Return array with PHOS clusters in aod format
 806
 807   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputPHOS()\n");
 808
 809   //Loop to select clusters in fiducial cut and fill container with aodClusters
 810   Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
 811   for (Int_t iclus = 0; iclus < nclusters; iclus++) {
 812     AliVCluster * clus = 0;
 813     if ( (clus = fInputEvent->GetCaloCluster(iclus)) ) {
 814       if (IsPHOSCluster(clus)){
 815         //Check if the cluster contains any bad channel and if close to calorimeter borders
 816         Int_t vindex = 0 ;
 817         if (fMixedEvent)
 818           vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
 819         if( GetCaloUtils()->ClusterContainsBadChannel("PHOS",clus->GetCellsAbsId(), clus->GetNCells()))
 820           continue;
 821         if(!GetCaloUtils()->CheckCellFiducialRegion(clus, fInputEvent->GetPHOSCells(), fInputEvent, vindex))
 822           continue;
 823
 824         TLorentzVector momentum ;
 825
 826         clus->GetMomentum(momentum, fVertex[vindex]);
 827
 828         if(fPHOSPtMin < momentum.Pt()){
 829
 830           if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"PHOS"))
 831             continue;
 832
 833           if(fDebug > 2 && momentum.E() > 0.1)
 834             printf("AliCaloTrackReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
 835                    momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 836
 837             //Recalibrate the cluster energy
 838           if(GetCaloUtils()->IsRecalibrationOn()) {
 839             Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliAODCaloCells*)GetPHOSCells());
 840             clus->SetE(energy);
 841           }
 842
 843           if (fMixedEvent) {
 844             clus->SetID(iclus) ;
 845           }
 846
 847           fAODPHOS->Add(clus);
 848
 849         }//Pt and Fiducial cut passed.
 850       }//PHOS cluster
 851     }//cluster exists
 852   }//esd cluster loop
 853
 854   //fAODPHOSNormalInputEntries = fAODPHOS->GetEntriesFast() ;
 855   if(fDebug > 1) printf("AliCaloTrackReader::FillInputPHOS()  - aod entries %d\n",  fAODPHOS->GetEntriesFast());//fAODPHOSNormalInputEntries);
 856
 857     //If second input event available, add the clusters.
 858     //  if(fSecondInputAODTree && fSecondInputAODEvent){
 859     //    GetSecondInputAODVertex(v);
 860     //    nclusters = ((AliAODEvent*)fSecondInputAODEvent)->GetNumberOfCaloClusters();
 861     //    if(fDebug > 1) printf("AliCaloTrackReader::FillInputPHOS()  - Add second input clusters, entries %d\n", nclusters);
 862     //          for (Int_t iclus =  0; iclus < nclusters; iclus++) {
 863     //                  AliAODCaloCluster * clus = 0;
 864     //                  if ( (clus = ((AliAODEvent*)fSecondInputAODEvent)->GetCaloCluster(iclus)) ) {
 865     //                          if (clus->IsPHOS()){
 866     //                                  TLorentzVector momentum ;
 867     //                                  clus->GetMomentum(momentum, v);
 868     //
 869     //                                  if(fPHOSPtMin < momentum.Pt()){
 870     //
 871     //                                          if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"PHOS")) continue;
 872     //
 873     //                                          if(fDebug > 2 && momentum.E() > 0.1) printf("AliCaloTrackReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
 874     //                                                                                                                                  momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
 875     //                                          fAODPHOS->Add(clus);
 876     //                                  }//Pt and Fiducial cut passed.
 877     //                          }//PHOS cluster
 878     //                  }// cluster exists
 879     //          }// cluster loop
 880     //          if(fDebug > 1) printf("AliCaloTrackReader::FillInputPHOS()  - aod normal entries %d, after second input %d\n", fAODPHOSNormalInputEntries, fAODPHOS->GetEntriesFast());
 881     //  }       //second input loop
 882
 883 }
 884
 885 //____________________________________________________________________________
 886 void AliCaloTrackReader::FillInputEMCALCells() {
 887     //Return array with EMCAL cells in aod format
 888
 889   fEMCALCells = fInputEvent->GetEMCALCells();
 890
 891 }
 892
 893 //____________________________________________________________________________
 894 void AliCaloTrackReader::FillInputPHOSCells() {
 895     //Return array with PHOS cells in aod format
 896
 897   fPHOSCells = fInputEvent->GetPHOSCells();
 898
 899 }
 900
 901
 902 //____________________________________________________________________________
 903 Bool_t AliCaloTrackReader::IsEMCALCluster(AliVCluster* cluster) const {
 904   // Check if it is a cluster from EMCAL. For old AODs cluster type has
 905   // different number and need to patch here
 906
 907   if(fDataType==kAOD && fOldAOD)
 908   {
 909     if (cluster->GetType() == 2) return kTRUE;
 910     else                         return kFALSE;
 911   }
 912   else
 913   {
 914     return cluster->IsEMCAL();
 915   }
 916
 917 }
 918
 919 //____________________________________________________________________________
 920 Bool_t AliCaloTrackReader::IsPHOSCluster(AliVCluster * cluster) const {
 921   //Check if it is a cluster from PHOS.For old AODs cluster type has
 922   // different number and need to patch here
 923
 924   if(fDataType==kAOD && fOldAOD)
 925   {
 926     Int_t type = cluster->GetType();
 927     if (type == 0 || type == 1) return kTRUE;
 928     else                        return kFALSE;
 929   }
 930   else
 931   {
 932     return cluster->IsPHOS();
 933   }
 934
 935 }
 936