]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG4/PartCorrBase/AliCaloTrackReader.cxx
Many fixes:
[u/mrichter/AliRoot.git] / PWG4 / PartCorrBase / AliCaloTrackReader.cxx
index 076381235f22799c374a23b8d367ae0a5120ebde..a331d006e13a43ebf99ba07334a47dfe03ee9d9f 100755 (executable)
@@ -42,6 +42,9 @@
 #include "AliVTrack.h"
 #include "AliVParticle.h"
 #include "AliMixedEvent.h"
+#include "AliESDtrack.h"
+#include "AliEMCALRecoUtils.h"
+#include "AliESDtrackCuts.h"
 
 ClassImp(AliCaloTrackReader)
   
@@ -60,109 +63,17 @@ ClassImp(AliCaloTrackReader)
 //    fSecondInputFileName(""),fSecondInputFirstEvent(0), 
 //    fAODCTSNormalInputEntries(0), fAODEMCALNormalInputEntries(0), 
 //    fAODPHOSNormalInputEntries(0), 
-    fTrackStatus(0), 
+    fTrackStatus(0),   fESDtrackCuts(0), fTrackMult(0), fTrackMultEtaCut(0.9),
     fReadStack(kFALSE), fReadAODMCParticles(kFALSE), 
     fDeltaAODFileName("deltaAODPartCorr.root"),fFiredTriggerClassName(""),
     fAnaLED(kFALSE),fTaskName(""),fCaloUtils(0x0), 
     fMixedEvent(NULL), fNMixedEvent(1), fVertex(NULL), 
-    fWriteOutputDeltaAOD(kFALSE),fOldAOD(kFALSE)
-{
+    fWriteOutputDeltaAOD(kFALSE),fOldAOD(kFALSE){
   //Ctor
   
   //Initialize parameters
   InitParameters();
 }
-/*
-//____________________________________________________________________________
-AliCaloTrackReader::AliCaloTrackReader(const AliCaloTrackReader & reader) :   
-  TObject(reader), fEventNumber(reader.fEventNumber), fCurrentFileName(reader.fCurrentFileName), 
-  fDataType(reader.fDataType), fDebug(reader.fDebug),
-  fFiducialCut(reader.fFiducialCut),
-  fComparePtHardAndJetPt(reader.fComparePtHardAndJetPt),
-  fPtHardAndJetPtFactor(reader.fPtHardAndJetPtFactor),
-  fCTSPtMin(reader.fCTSPtMin), fEMCALPtMin(reader.fEMCALPtMin),fPHOSPtMin(reader.fPHOSPtMin), 
-  fAODBranchList(new TList()),
-  fAODCTS(new TObjArray(*reader.fAODCTS)),  
-  fAODEMCAL(new TObjArray(*reader.fAODEMCAL)),
-  fAODPHOS(new TObjArray(*reader.fAODPHOS)),
-  fEMCALCells(new TNamed(*reader.fEMCALCells)),
-  fPHOSCells(new TNamed(*reader.fPHOSCells)),
-  fInputEvent(reader.fInputEvent), fOutputEvent(reader.fOutputEvent), fMC(reader.fMC),
-  fFillCTS(reader.fFillCTS),fFillEMCAL(reader.fFillEMCAL),fFillPHOS(reader.fFillPHOS),
-  fFillEMCALCells(reader.fFillEMCALCells),fFillPHOSCells(reader.fFillPHOSCells),
-  fSecondInputAODTree(reader.fSecondInputAODTree), 
-  fSecondInputAODEvent(reader.fSecondInputAODEvent),
-  fSecondInputFileName(reader.fSecondInputFileName), 
-  fSecondInputFirstEvent(reader.fSecondInputFirstEvent),
-  fAODCTSNormalInputEntries(reader.fAODCTSNormalInputEntries), 
-  fAODEMCALNormalInputEntries(reader.fAODEMCALNormalInputEntries), 
-  fAODPHOSNormalInputEntries(reader.fAODPHOSNormalInputEntries),
-  fTrackStatus(reader.fTrackStatus),
-  fReadStack(reader.fReadStack), fReadAODMCParticles(reader.fReadAODMCParticles),
-  fFiredTriggerClassName(reader.fFiredTriggerClassName),
-  fAnaLED(reader.fAnaLED), 
-  fTaskName(reader.fTaskName),
-  fCaloUtils(new AliCalorimeterUtils(*reader.fCaloUtils))
-{
-  // cpy ctor  
-}
-*/
-//_________________________________________________________________________
-//AliCaloTrackReader & AliCaloTrackReader::operator = (const AliCaloTrackReader & source)
-//{
-//  // assignment operator
-//  
-//  if(&source == this) return *this;
-//  
-//  fDataType    = source.fDataType ;
-//  fDebug       = source.fDebug ;
-//  fEventNumber = source.fEventNumber ;
-//  fCurrentFileName = source.fCurrentFileName ;
-//  fFiducialCut = source.fFiducialCut;
-//     
-//  fComparePtHardAndJetPt = source.fComparePtHardAndJetPt;
-//  fPtHardAndJetPtFactor  = source.fPtHardAndJetPtFactor;
-//     
-//  fCTSPtMin    = source.fCTSPtMin ;
-//  fEMCALPtMin  = source.fEMCALPtMin ;
-//  fPHOSPtMin   = source.fPHOSPtMin ; 
-//  
-//  fAODCTS     = new TObjArray(*source.fAODCTS) ;
-//  fAODEMCAL   = new TObjArray(*source.fAODEMCAL) ;
-//  fAODPHOS    = new TObjArray(*source.fAODPHOS) ;
-//  fEMCALCells = new TNamed(*source.fEMCALCells) ;
-//  fPHOSCells  = new TNamed(*source.fPHOSCells) ;
-//
-//  fInputEvent  = source.fInputEvent;
-//  fOutputEvent = source.fOutputEvent;
-//  fMC          = source.fMC;
-//  
-//  fFillCTS        = source.fFillCTS;
-//  fFillEMCAL      = source.fFillEMCAL;
-//  fFillPHOS       = source.fFillPHOS;
-//  fFillEMCALCells = source.fFillEMCALCells;
-//  fFillPHOSCells  = source.fFillPHOSCells;
-//
-//  fSecondInputAODTree    = source.fSecondInputAODTree;
-//  fSecondInputAODEvent   = source.fSecondInputAODEvent;
-//  fSecondInputFileName   = source.fSecondInputFileName;
-//  fSecondInputFirstEvent = source.fSecondInputFirstEvent;
-//
-//  fAODCTSNormalInputEntries   = source.fAODCTSNormalInputEntries; 
-//  fAODEMCALNormalInputEntries = source.fAODEMCALNormalInputEntries; 
-//  fAODPHOSNormalInputEntries  = source.fAODPHOSNormalInputEntries;
-//     
-//  fTrackStatus        = source.fTrackStatus;
-//  fReadStack          = source.fReadStack;
-//  fReadAODMCParticles = source.fReadAODMCParticles;  
-//     
-//  fDeltaAODFileName   = source.fDeltaAODFileName;
-//     
-//  fFiredTriggerClassName = source.fFiredTriggerClassName  ;
-//     
-//  return *this;
-//  
-//}
 
 //_________________________________
 AliCaloTrackReader::~AliCaloTrackReader() {
@@ -176,17 +87,20 @@ AliCaloTrackReader::~AliCaloTrackReader() {
   }  
   
   if(fAODCTS){
-    fAODCTS->Clear() ; 
+    if(fDataType!=kMC)fAODCTS->Clear() ; 
+    else              fAODCTS->Delete() ; 
     delete fAODCTS ;
   }
   
   if(fAODEMCAL){
-    fAODEMCAL->Clear() ; 
+    if(fDataType!=kMC)fAODEMCAL->Clear() ; 
+    else              fAODEMCAL->Delete() ; 
     delete fAODEMCAL ;
   }
   
   if(fAODPHOS){
-    fAODPHOS->Clear() ; 
+    if(fDataType!=kMC)fAODPHOS->Clear() ; 
+    else              fAODPHOS->Delete() ; 
     delete fAODPHOS ;
   }
   
@@ -199,13 +113,17 @@ AliCaloTrackReader::~AliCaloTrackReader() {
     fPHOSCells->Clear() ; 
     delete fPHOSCells ;
   }
-  if (fMixedEvent) {
+
+  if(fVertex){
     for (Int_t i = 0; i < fNMixedEvent; i++) {
-      delete [] fVertex[i] ; 
+      delete [] fVertex[i] ;
+
     }
-  }
-  delete [] fVertex ;
-       
+    delete [] fVertex ;
+       }
+
+  if(fESDtrackCuts)   delete fESDtrackCuts;
+  
 //  Pointers not owned, done by the analysis frame
 //  if(fInputEvent)  delete fInputEvent ;
 //  if(fOutputEvent) delete fOutputEvent ;
@@ -230,13 +148,13 @@ Bool_t AliCaloTrackReader::ComparePtHardAndJetPt(){
        if(!fReadStack) return kTRUE; //Information not filtered to AOD
        
        if(!strcmp(GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader")){
-               TParticle * jet =  new TParticle;
+               TParticle * jet =  0;
                AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetGenEventHeader();
                Int_t nTriggerJets =  pygeh->NTriggerJets();
                Float_t ptHard = pygeh->GetPtHard();
                
                //if(fDebug > 1) printf("AliMCAnalysisUtils::PythiaEventHeader: Njets: %d, pT Hard %f\n",nTriggerJets, ptHard);
-           Float_t tmpjet[]={0,0,0,0};
+    Float_t tmpjet[]={0,0,0,0};
                for(Int_t ijet = 0; ijet< nTriggerJets; ijet++){
                        pygeh->TriggerJet(ijet, tmpjet);
                        jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
@@ -248,8 +166,9 @@ Bool_t AliCaloTrackReader::ComparePtHardAndJetPt(){
                                return kFALSE;
                        }
                }
+    if(jet) delete jet; 
        }
-       
+         
        return kTRUE ;
        
 }
@@ -397,6 +316,14 @@ void AliCaloTrackReader::InitParameters()
   fFiredTriggerClassName      = "";
                
   fAnaLED = kFALSE;
+  
+  //We want tracks fitted in the detectors:
+  //fTrackStatus=AliESDtrack::kTPCrefit;
+  //fTrackStatus|=AliESDtrack::kITSrefit;  
+  
+  fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
+
+  
 }
 
 //________________________________________________________________
@@ -419,6 +346,7 @@ void AliCaloTrackReader::Print(const Option_t * opt) const
   printf("Use EMCAL Cells =     %d\n", fFillEMCALCells) ;
   printf("Use PHOS  Cells =     %d\n", fFillPHOSCells) ;
   printf("Track status    =     %d\n", (Int_t) fTrackStatus) ;
+  printf("Track Mult Eta Cut =  %d\n", (Int_t) fTrackMultEtaCut) ;
   printf("Write delta AOD =     %d\n", fWriteOutputDeltaAOD) ;
 
   if(fComparePtHardAndJetPt)
@@ -500,16 +428,12 @@ Bool_t AliCaloTrackReader::FillInputEvent(const Int_t iEntry, const char * curre
 //    
 //  }
        
-
-  for (Int_t iev = 0; iev < fNMixedEvent; iev++) {
-    if (!fMixedEvent) {
-      GetVertex() ;
-    }      
-    else 
-      fMixedEvent->GetVertexOfEvent(iev)->GetXYZ(fVertex[iev]);     
-  }
+  //Fill Vertex array
   
-  if(fFillEMCALCells) 
+  FillVertexArray();
+  
+  //Fill the arrays with cluster/tracks/cells data
+   if(fFillEMCALCells) 
     FillInputEMCALCells();
   if(fFillPHOSCells)  
     FillInputPHOSCells();
@@ -529,12 +453,11 @@ Bool_t AliCaloTrackReader::FillInputEvent(const Int_t iEntry, const char * curre
 void AliCaloTrackReader::ResetLists() {
   //  Reset lists, called by the analysis maker 
 
-  if(fAODCTS)   fAODCTS -> Clear();
-  if(fAODEMCAL) fAODEMCAL -> Clear();
-  if(fAODPHOS)  fAODPHOS -> Clear();
+  if(fAODCTS)     fAODCTS     -> Clear();
+  if(fAODEMCAL)   fAODEMCAL   -> Clear();
+  if(fAODPHOS)    fAODPHOS    -> Clear();
   if(fEMCALCells) fEMCALCells -> Clear();
-  if(fPHOSCells)  fPHOSCells -> Clear();
-
+  if(fPHOSCells)  fPHOSCells  -> Clear();
 }
 
 //____________________________________________________________________________
@@ -545,6 +468,15 @@ void AliCaloTrackReader::SetInputEvent(AliVEvent* const input)
   if (fMixedEvent) {
     fNMixedEvent = fMixedEvent->GetNumberOfEvents() ; 
   }
+
+  //Delete previous vertex
+  if(fVertex){
+    for (Int_t i = 0; i < fNMixedEvent; i++) {
+      delete [] fVertex[i] ; 
+    }
+    delete [] fVertex ;
+  }
+  
   fVertex = new Double_t*[fNMixedEvent] ; 
   for (Int_t i = 0; i < fNMixedEvent; i++) {
     fVertex[i] = new Double_t[3] ; 
@@ -555,19 +487,72 @@ void AliCaloTrackReader::SetInputEvent(AliVEvent* const input)
 }
 
 //____________________________________________________________________________
-Double_t * AliCaloTrackReader::GetVertex() {
-    //Return vertex position
-  if (fMixedEvent)
-    return NULL ; 
-  if(fInputEvent->GetPrimaryVertex())
-    fInputEvent->GetPrimaryVertex()->GetXYZ(fVertex[0]); 
-  else {
-    printf("AliCaloTrackReader::GetVertex() - No vertex available, InputEvent()->GetPrimaryVertex() = 0, return (0,0,0)\n");
-    fVertex[0][0]=0.0; fVertex[0][1]=0.0; fVertex[0][2]=0.0;
+void AliCaloTrackReader::GetVertex(Double_t vertex[3]) const {
+  //Return vertex position to be used for single event analysis
+  vertex[0]=fVertex[0][0];  
+  vertex[1]=fVertex[0][1];  
+  vertex[2]=fVertex[0][2];
+}
+
+//____________________________________________________________________________
+void AliCaloTrackReader::GetVertex(Double_t vertex[3], const Int_t evtIndex) const {
+  //Return vertex position for mixed event, recover the vertex in a particular event.
+  
+  //Int_t evtIndex = 0; // for single events only one vertex stored in position 0, default value
+  //if (fMixedEvent && clusterID >=0) {
+  //  evtIndex=GetMixedEvent()->EventIndexForCaloCluster(clusterID) ; 
+  //}
+  
+  vertex[0]=fVertex[evtIndex][0];  vertex[1]=fVertex[evtIndex][1];  vertex[2]=fVertex[evtIndex][2];
+  
+}
+//
+
+
+//____________________________________________________________________________
+void AliCaloTrackReader::FillVertexArray() {
+  
+  //Fill data member with vertex
+  //In case of Mixed event, multiple vertices
+  
+  //Delete previous vertex
+  if(fVertex){
+    for (Int_t i = 0; i < fNMixedEvent; i++) {
+      delete [] fVertex[i] ; 
+    }
+    delete [] fVertex ;  
+  }
+  
+  fVertex = new Double_t*[fNMixedEvent] ; 
+  for (Int_t i = 0; i < fNMixedEvent; i++) {
+    fVertex[i] = new Double_t[3] ; 
+    fVertex[i][0] = 0.0 ; 
+    fVertex[i][1] = 0.0 ; 
+    fVertex[i][2] = 0.0 ; 
+  }          
+  
+  if (!fMixedEvent) { //Single event analysis
+    
+    if(fDataType!=kMC)fInputEvent->GetPrimaryVertex()->GetXYZ(fVertex[0]); 
+    else { 
+      fVertex[0][0]=0.;   fVertex[0][1]=0.;   fVertex[0][2]=0.;
+    }
+      
+    if(fDebug > 1)
+      printf("AliCaloTrackReader::FillVertexArray() - Single Event Vertex : %f,%f,%f\n",fVertex[0][0],fVertex[0][1],fVertex[0][2]);
+
+  } else { // MultiEvent analysis
+    for (Int_t iev = 0; iev < fNMixedEvent; iev++) {
+      fMixedEvent->GetVertexOfEvent(iev)->GetXYZ(fVertex[iev]);
+      if(fDebug > 1)
+        printf("AliCaloTrackReader::FillVertexArray() - Multi Event %d Vertex : %f,%f,%f\n",iev,fVertex[iev][0],fVertex[iev][1],fVertex[iev][2]);
+
+    }
   }
-  return fVertex[0] ; 
+  
 }
 
+
 //____________________________________________________________________________
 void AliCaloTrackReader::FillInputCTS() {
   //Return array with Central Tracking System (CTS) tracks
@@ -576,13 +561,23 @@ void AliCaloTrackReader::FillInputCTS() {
   
   Int_t nTracks   = fInputEvent->GetNumberOfTracks() ;
   Double_t p[3];
+  fTrackMult = 0;
+  Int_t nstatus = 0;
   for (Int_t itrack =  0; itrack <  nTracks; itrack++) {////////////// track loop
     AliVTrack * track = (AliVTrack*)fInputEvent->GetTrack(itrack) ; // retrieve track from esd
-    
-      //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
+
+    //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
     if (fTrackStatus && !((track->GetStatus() & fTrackStatus) == fTrackStatus)) 
       continue ;
     
+    nstatus++;
+    
+    if(fDataType==kESD && !fESDtrackCuts->AcceptTrack((AliESDtrack*)track)) continue;
+    
+    //Count the tracks in eta < 0.9
+    //printf("Eta %f cut  %f\n",TMath::Abs(track->Eta()),fTrackMultEtaCut);
+    if(TMath::Abs(track->Eta())< fTrackMultEtaCut) fTrackMult++;
+    
     track->GetPxPyPz(p) ;
     TLorentzVector momentum(p[0],p[1],p[2],0);
     
@@ -605,7 +600,8 @@ void AliCaloTrackReader::FillInputCTS() {
   }// track loop
        
   //fAODCTSNormalInputEntries = fAODCTS->GetEntriesFast();
-  if(fDebug > 1) printf("AliCaloTrackReader::FillInputCTS()   - aod entries %d\n", fAODCTS->GetEntriesFast());//fAODCTSNormalInputEntries);
+  if(fDebug > 1) 
+    printf("AliCaloTrackReader::FillInputCTS()   - aod entries %d, input tracks %d, pass status %d, multipliticy %d\n", fAODCTS->GetEntriesFast(), nTracks, nstatus, fTrackMult);//fAODCTSNormalInputEntries);
   
     //  //If second input event available, add the clusters.
     //  if(fSecondInputAODTree && fSecondInputAODEvent){
@@ -673,10 +669,31 @@ void AliCaloTrackReader::FillInputEMCAL() {
             printf("AliCaloTrackReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
                    momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
           
+          //Float_t pos[3];
+          //clus->GetPosition(pos);
+          //printf("Before Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
+          
           //Recalibrate the cluster energy 
           if(GetCaloUtils()->IsRecalibrationOn()) {
-            Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliAODCaloCells*)GetEMCALCells());
+            Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, GetEMCALCells());
             clus->SetE(energy);
+            //printf("Recalibrated Energy %f\n",clus->E());  
+            GetCaloUtils()->RecalculateClusterShowerShapeParameters(GetEMCALCells(),clus);
+            GetCaloUtils()->RecalculateClusterPID(clus);
+
+          }
+          
+          //Correct non linearity
+          if(GetCaloUtils()->IsCorrectionOfClusterEnergyOn()){
+            GetCaloUtils()->CorrectClusterEnergy(clus) ;
+            //printf("Linearity Corrected Energy %f\n",clus->E());  
+          }
+          
+          //Recalculate cluster position
+          if(GetCaloUtils()->IsRecalculationOfClusterPositionOn()){
+            GetCaloUtils()->RecalculateClusterPosition(GetEMCALCells(),clus); 
+            //clus->GetPosition(pos);
+            //printf("After  Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
           }
           
           if (fMixedEvent) {
@@ -854,5 +871,3 @@ Bool_t AliCaloTrackReader::IsPHOSCluster(AliVCluster * cluster) const {
   
 }
 
-
-