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