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