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