]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG4/PartCorrBase/AliIsolationCut.cxx
Correct and clean the vertex retrieval in case of SE or ME analysis
[u/mrichter/AliRoot.git] / PWG4 / PartCorrBase / AliIsolationCut.cxx
index 862e326a73ce66df1fe8cfe37ea23b437b352590..63408b27d2aa488c46363b9d7f28d0d0a4e3391d 100755 (executable)
 #include "AliIsolationCut.h" 
 #include "AliAODPWG4ParticleCorrelation.h"
 #include "AliAODTrack.h"
-#include "AliAODCaloCluster.h"
+#include "AliVCluster.h"
+#include "AliCaloTrackReader.h"
+#include "AliMixedEvent.h"
 
 ClassImp(AliIsolationCut)
   
 //____________________________________________________________________________
   AliIsolationCut::AliIsolationCut() : 
     TObject(),
-    fConeSize(0.),fPtThreshold(0.), fPtFraction(0.), fICMethod(0)
+    fConeSize(0.),fPtThreshold(0.), fPtFraction(0.), fICMethod(0),fPartInCone(0)
  
 {
   //default ctor
@@ -50,7 +52,7 @@ ClassImp(AliIsolationCut)
   InitParameters();
 
 }
-
+/*
 //____________________________________________________________________________
 AliIsolationCut::AliIsolationCut(const AliIsolationCut & g) : 
   TObject(g),
@@ -78,26 +80,29 @@ AliIsolationCut & AliIsolationCut::operator = (const AliIsolationCut & source)
   return *this;
   
 }
-
+*/
 //____________________________________________________________________________
 TString AliIsolationCut::GetICParametersList()
 {
   //Put data member values in string to keep in output container
   
   TString parList ; //this will be list of parameters used for this analysis.
-  char onePar[255] ;
+  const Int_t buffersize = 255;
+  char onePar[buffersize] ;
   
-  sprintf(onePar,"--- AliIsolationCut ---\n") ;
+  snprintf(onePar,buffersize,"--- AliIsolationCut ---\n") ;
   parList+=onePar ;    
-  sprintf(onePar,"fConeSize: (isolation cone size) %1.2f\n",fConeSize) ;
+  snprintf(onePar,buffersize,"fConeSize: (isolation cone size) %1.2f\n",fConeSize) ;
   parList+=onePar ;
-  sprintf(onePar,"fPtThreshold =%1.2f (isolation pt threshold) \n",fPtThreshold) ;
+  snprintf(onePar,buffersize,"fPtThreshold =%1.2f (isolation pt threshold) \n",fPtThreshold) ;
   parList+=onePar ;
-  sprintf(onePar,"fPtFraction=%1.2f (isolation pt threshold fraction ) \n",fPtFraction) ;
+  snprintf(onePar,buffersize,"fPtFraction=%1.2f (isolation pt threshold fraction ) \n",fPtFraction) ;
   parList+=onePar ;
-  sprintf(onePar,"fICMethod=%d (isolation cut case) \n",fICMethod) ;
+  snprintf(onePar,buffersize,"fICMethod=%d (isolation cut case) \n",fICMethod) ;
   parList+=onePar ;
-  
+  snprintf(onePar,buffersize,"fPartInCone=%d \n",fPartInCone) ;
+  parList+=onePar ;
+
   return parList; 
 }
 
@@ -106,43 +111,40 @@ void AliIsolationCut::InitParameters()
 {
   //Initialize the parameters of the analysis.
   
-  fConeSize             = 0.4 ; 
-  fPtThreshold         = 1. ; 
-  fPtFraction        = 0.1 ; 
-  
-  fICMethod = kPtThresIC; // 0 pt threshol method, 1 cone pt sum method
+  fConeSize    = 0.4 ; 
+  fPtThreshold = 1. ; 
+  fPtFraction  = 0.1 ; 
+  fPartInCone  = kNeutralAndCharged;
+  fICMethod    = kPtThresIC; // 0 pt threshol method, 1 cone pt sum method
   
 }
 
 //__________________________________________________________________
-void  AliIsolationCut::MakeIsolationCut(TObjArray * plCTS,  TObjArray * plNe, Double_t * vertex
+void  AliIsolationCut::MakeIsolationCut(TObjArray * const plCTS,  TObjArray * const plNe, AliCaloTrackReader * const reader
                                        const Bool_t fillAOD, AliAODPWG4ParticleCorrelation  *pCandidate, 
                                        const TString aodArrayRefName,
                                        Int_t & n, Int_t & nfrac, Float_t &coneptsum,  Bool_t  &isolated) const
 {  
   //Search in cone around a candidate particle if it is isolated 
   Float_t phiC  = pCandidate->Phi() ;
-  Float_t etaC = pCandidate->Eta() ;
-  Float_t ptC = pCandidate->Pt() ;
-  Float_t pt     = -100. ;
+  Float_t etaC  = pCandidate->Eta() ;
+  Float_t ptC   = pCandidate->Pt() ;
+  Float_t pt    = -100. ;
   Float_t eta   = -100.  ;
-  Float_t phi    = -100.  ;
+  Float_t phi   = -100.  ;
   Float_t rad   = -100 ;
   n = 0 ;
   coneptsum = 0.; 
   isolated = kFALSE;
-
+  
   //Initialize the array with refrences
   TObjArray * refclusters = 0x0;
-  TObjArray * reftracks    =0x0;
-
-  if(fillAOD) {
-    refclusters = new TObjArray;
-    reftracks    = new TObjArray;
-  }
-
+  TObjArray * reftracks   = 0x0;
+  Int_t ntrackrefs   = 0;
+  Int_t nclusterrefs = 0;
+  
   //Check charged particles in cone.
-  if(plCTS){
+  if(plCTS && (fPartInCone==kOnlyCharged || fPartInCone==kNeutralAndCharged)){
     TVector3 p3;
     for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ ){
       AliAODTrack* track = (AliAODTrack *)(plCTS->At(ipr)) ; 
@@ -158,29 +160,57 @@ void  AliIsolationCut::MakeIsolationCut(TObjArray * plCTS,  TObjArray * plNe, Do
       rad = TMath::Sqrt((eta-etaC)*(eta-etaC)+ (phi-phiC)*(phi-phiC));
       
       if(rad < fConeSize){
-       if(fillAOD) {
-         reftracks->Add(track);
-       }
-       //printf("charged in isolation cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
-       coneptsum+=pt;
-       if(pt > fPtThreshold ) n++;
-       if(pt > fPtFraction*ptC ) nfrac++;  
+        if(fillAOD) {
+          ntrackrefs++;
+          if(ntrackrefs == 1){
+            reftracks = new TObjArray(0);
+            reftracks->SetName(aodArrayRefName+"Tracks");
+            reftracks->SetOwner(kFALSE);
+          }
+          reftracks->Add(track);
+        }
+        //printf("charged in isolation cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
+        coneptsum+=pt;
+        if(pt > fPtThreshold ) n++;
+        if(pt > fPtFraction*ptC ) nfrac++;  
       }
     }// charged particle loop
   }//Tracks
   
   //Check neutral particles in cone.  
-  if(plNe){
+  if(plNe && (fPartInCone==kOnlyNeutral || fPartInCone==kNeutralAndCharged)){
+         
+    //Get vertex for photon momentum calculation
+    //Double_t vertex2[] = {0,0,0} ; //vertex second AOD input ;
+    //if(reader->GetDataType()!= AliCaloTrackReader::kMC) 
+    //{
+      //if(reader->GetSecondInputAODTree()) reader->GetSecondInputAODVertex(vertex2);
+    //}
     TLorentzVector mom ;
     for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ ){
-      AliAODCaloCluster * calo = (AliAODCaloCluster *)(plNe->At(ipr)) ;
+      AliVCluster * calo = (AliVCluster *)(plNe->At(ipr)) ;
+      
+      //Get the index where the cluster comes, to retrieve the corresponding vertex
+      Int_t evtIndex = 0 ; 
+      if (reader->GetMixedEvent()) {
+        evtIndex=reader->GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ; 
+      }
       
       //Do not count the candidate (photon or pi0) or the daughters of the candidate
       if(calo->GetID() == pCandidate->GetCaloLabel(0) || calo->GetID() == pCandidate->GetCaloLabel(1)) continue ;      //Skip matched clusters with tracks
       
       if(calo->GetNTracksMatched() > 0) continue ; 
       
-      calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
+      //Input from second AOD?
+      //Int_t input = 0;
+      //      if     (pCandidate->GetDetector() == "EMCAL" && reader->GetAODEMCALNormalInputEntries() <= ipr) input = 1 ;
+      //      else if(pCandidate->GetDetector() == "PHOS"  && reader->GetAODPHOSNormalInputEntries()  <= ipr) input = 1;
+      
+      //Get Momentum vector, 
+      //if     (input == 0) 
+      calo->GetMomentum(mom,reader->GetVertex(evtIndex)) ;//Assume that come from vertex in straight line
+      //else if(input == 1) calo->GetMomentum(mom,vertex2);//Assume that come from vertex in straight line  
+      
       pt   = mom.Pt();
       eta  = mom.Eta();
       phi  = mom.Phi() ;
@@ -189,13 +219,19 @@ void  AliIsolationCut::MakeIsolationCut(TObjArray * plCTS,  TObjArray * plNe, Do
       //Check if there is any particle inside cone with pt larger than  fPtThreshold
       rad = TMath::Sqrt((eta-etaC)*(eta-etaC)+ (phi-phiC)*(phi-phiC));
       if(rad < fConeSize){
-       if(fillAOD) {
-         refclusters->Add(calo);
-       }
-       //printf("neutral in isolation cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
-       coneptsum+=pt;
-       if(pt > fPtThreshold ) n++;
-       if(pt > fPtFraction*ptC ) nfrac++;
+        if(fillAOD) {
+          nclusterrefs++;
+          if(nclusterrefs==1){
+            refclusters = new TObjArray(0);
+            refclusters->SetName(aodArrayRefName+"Clusters");
+            refclusters->SetOwner(kFALSE);
+          }
+          refclusters->Add(calo);
+        }
+        //printf("neutral in isolation cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
+        coneptsum+=pt;
+        if(pt > fPtThreshold ) n++;
+        if(pt > fPtFraction*ptC ) nfrac++;
       }//in cone
     }// neutral particle loop
   }//neutrals
@@ -204,16 +240,10 @@ void  AliIsolationCut::MakeIsolationCut(TObjArray * plCTS,  TObjArray * plNe, Do
   
   //Add reference arrays to AOD when filling AODs only
   if(fillAOD) {
-       if(refclusters->GetEntriesFast() > 0){ 
-               refclusters->SetName(aodArrayRefName+"Clusters");
-               pCandidate->AddObjArray(refclusters);
-       }
-       if(reftracks->GetEntriesFast()   > 0){
-               reftracks->SetName(aodArrayRefName+"Tracks");
-               pCandidate->AddObjArray(reftracks);
-       } 
+    if(refclusters)    pCandidate->AddObjArray(refclusters);
+    if(reftracks)        pCandidate->AddObjArray(reftracks);
   }
-
+  
   //Check isolation, depending on method.
   if( fICMethod == kPtThresIC){
     if(n==0) isolated = kTRUE ;
@@ -229,6 +259,10 @@ void  AliIsolationCut::MakeIsolationCut(TObjArray * plCTS,  TObjArray * plNe, Do
     if(coneptsum < fPtFraction*ptC)
       isolated  =  kTRUE ;
   }
+  
+  //if(refclusters) delete refclusters;
+  //if(reftracks)   delete reftracks;
+  
 }
 
 //__________________________________________________________________
@@ -245,7 +279,7 @@ void AliIsolationCut::Print(const Option_t * opt) const
   printf("Cone Size          =     %1.2f\n", fConeSize) ; 
   printf("pT threshold       =     %2.1f\n", fPtThreshold) ;
   printf("pT fraction        =     %3.1f\n", fPtFraction) ;
-
+  printf("particle type in cone =  %d\n",fPartInCone);
   printf("    \n") ;
   
 }