]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG4/PartCorrBase/AliIsolationCut.cxx
Fix Coverity
[u/mrichter/AliRoot.git] / PWG4 / PartCorrBase / AliIsolationCut.cxx
index 503293deac3ebb717133ae3948b1ec06a90bb9c3..a46fb9c67c7cb0f0cb9dc45f13058a31ad145e06 100755 (executable)
 //
 //
 //*-- Author: Gustavo Conesa (LNF-INFN) 
+
+//-Yaxian Mao (add the possibility for different IC method with different pt range, 01/10/2010)
+//-Yaxian Mao (check the candidate particle is the leading particle or not at the same hemishere)
+
 //////////////////////////////////////////////////////////////////////////////
   
   
 // --- ROOT system --- 
-//#include <Riostream.h>
 #include <TLorentzVector.h>
 #include <TObjArray.h>
 
 #include "AliIsolationCut.h" 
 #include "AliAODPWG4ParticleCorrelation.h"
 #include "AliAODTrack.h"
-#include "AliAODCaloCluster.h"
+#include "AliVCluster.h"
 #include "AliCaloTrackReader.h"
+#include "AliMixedEvent.h"
+#include "AliCaloPID.h"
 
 ClassImp(AliIsolationCut)
   
-//____________________________________________________________________________
-  AliIsolationCut::AliIsolationCut() : 
-    TObject(),
-    fConeSize(0.),fPtThreshold(0.), fPtFraction(0.), fICMethod(0),fPartInCone(0)
+//____________________________________
+AliIsolationCut::AliIsolationCut() : 
+TObject(),
+fConeSize(0.),
+fPtThreshold(0.), 
+fSumPtThreshold(0.), 
+fPtFraction(0.), 
+fICMethod(0),
+fPartInCone(0)
+
 {
   //default ctor
   
   //Initialize parameters
   InitParameters();
-
-}
-/*
-//____________________________________________________________________________
-AliIsolationCut::AliIsolationCut(const AliIsolationCut & g) : 
-  TObject(g),
-  fConeSize(g.fConeSize),
-  fPtThreshold(g.fPtThreshold),
-  fPtFraction(g.fPtFraction), 
-  fICMethod(g.fICMethod)
-{
-  // cpy ctor
   
 }
 
-//_________________________________________________________________________
-AliIsolationCut & AliIsolationCut::operator = (const AliIsolationCut & source)
-{
-  // assignment operator
-  
-  if(&source == this) return *this;
-  
-  fConeSize = source.fConeSize ;
-  fPtThreshold = source.fPtThreshold ; 
-  fICMethod = source.fICMethod ;
-  fPtFraction = source.fPtFraction ;
-
-  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 ;
-  sprintf(onePar,"fPartInCone=%d \n",fPartInCone) ;
+  snprintf(onePar,buffersize,"fPartInCone=%d \n",fPartInCone) ;
   parList+=onePar ;
-
+  
   return parList; 
 }
 
-//____________________________________________________________________________
+//____________________________________
 void AliIsolationCut::InitParameters()
 {
   //Initialize the parameters of the analysis.
   
-  fConeSize    = 0.4 ; 
-  fPtThreshold = 1. ; 
-  fPtFraction  = 0.1 ; 
-  fPartInCone  = kNeutralAndCharged;
-  fICMethod    = kPtThresIC; // 0 pt threshol method, 1 cone pt sum method
+  fConeSize       = 0.4 ; 
+  fPtThreshold    = 1.  ; 
+  fSumPtThreshold = 0.5 ; 
+  fPtFraction     = 0.1 ; 
+  fPartInCone     = kOnlyCharged;
+  fICMethod       = kSumPtFracIC; // 0 pt threshol method, 1 cone pt sum method
   
 }
 
-//__________________________________________________________________
-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
+//________________________________________________________________________________
+void  AliIsolationCut::MakeIsolationCut(const TObjArray * plCTS, 
+                                        const TObjArray * plNe, 
+                                        const AliCaloTrackReader * reader, 
+                                        const AliCaloPID * pid, 
+                                        const Bool_t bFillAOD, 
+                                        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() ;
+  if(phiC<0) phiC+=TMath::TwoPi();
   Float_t etaC  = pCandidate->Eta() ;
   Float_t ptC   = pCandidate->Pt() ;
   Float_t pt    = -100. ;
-  Float_t eta   = -100.  ;
-  Float_t phi   = -100.  ;
-  Float_t rad   = -100 ;
-  n = 0 ;
+  Float_t eta   = -100. ;
+  Float_t phi   = -100. ;
+  Float_t rad   = -100. ;
+  
+  n         = 0 ;
+  nfrac     = 0 ;
   coneptsum = 0.; 
-  isolated = kFALSE;
-
+  isolated  = kFALSE;
+  
   //Initialize the array with refrences
   TObjArray * refclusters = 0x0;
   TObjArray * reftracks   = 0x0;
@@ -147,83 +141,133 @@ void  AliIsolationCut::MakeIsolationCut(TObjArray * const plCTS,  TObjArray * co
     for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ ){
       AliAODTrack* track = (AliAODTrack *)(plCTS->At(ipr)) ; 
       //Do not count the candidate (pion, conversion photon) or the daughters of the candidate
-      if(track->GetID() == pCandidate->GetTrackLabel(0) || track->GetID() == pCandidate->GetTrackLabel(1)) continue ;
+      if(track->GetID() == pCandidate->GetTrackLabel(0) || track->GetID() == pCandidate->GetTrackLabel(1) 
+         || track->GetID() == pCandidate->GetTrackLabel(2) || track->GetID() == pCandidate->GetTrackLabel(3) 
+         ) continue ;
       p3.SetXYZ(track->Px(),track->Py(),track->Pz());
       pt   = p3.Pt();
       eta  = p3.Eta();
       phi  = p3.Phi() ;
       if(phi<0) phi+=TMath::TwoPi();
       
+      //only loop the particle at the same side of candidate
+      if(TMath::Abs(phi-phiC)>TMath::PiOver2()) continue ;
+      //if at the same side has particle larger than candidate, then candidate can not be the leading, skip such events
+      if(pt > ptC){
+        n         = -1;
+        nfrac     = -1;
+        coneptsum = -1;
+        isolated  = kFALSE;
+        if(bFillAOD && reftracks) {
+          reftracks->Clear(); 
+          delete reftracks;
+        }
+        return ;
+      }
+      
       //Check if there is any particle inside cone with pt larger than  fPtThreshold
-      rad = TMath::Sqrt((eta-etaC)*(eta-etaC)+ (phi-phiC)*(phi-phiC));
+
+      rad = Radius(etaC, phiC, eta, phi);
       
       if(rad < fConeSize){
-       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++;  
-      }
+        if(bFillAOD) {
+          ntrackrefs++;
+          if(ntrackrefs == 1){
+            reftracks = new TObjArray(0);
+            //reftracks->SetName(Form("Tracks%s",aodArrayRefName.Data()));
+            TString tempo(aodArrayRefName)  ; 
+            tempo += "Tracks" ; 
+            reftracks->SetName(tempo);
+            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++;  
+      } // Inside cone
     }// charged particle loop
   }//Tracks
   
   //Check neutral particles in cone.  
   if(plNe && (fPartInCone==kOnlyNeutral || fPartInCone==kNeutralAndCharged)){
          
-       //Get vertex for photon momentum calculation
-       Double_t vertex[]  = {0,0,0} ; //vertex ;
-       Double_t vertex2[] = {0,0,0} ; //vertex second AOD input ;
-       if(!reader->GetDataType()== AliCaloTrackReader::kMC) 
-       {
-               reader->GetVertex(vertex);
-               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)) ;
       
-      //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
+      //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()) ; 
+      }
       
-      if(calo->GetNTracksMatched() > 0) continue ; 
-      //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;
+      //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 ;      
       
-      //Get Momentum vector, 
-      if     (input == 0) calo->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
-      else if(input == 1) calo->GetMomentum(mom,vertex2);//Assume that come from vertex in straight line  
+      //Skip matched clusters with tracks
+      if( pid->IsTrackMatched(calo,reader->GetCaloUtils(),reader->GetInputEvent()) ) continue ;
+    
+      //Assume that come from vertex in straight line
+      calo->GetMomentum(mom,reader->GetVertex(evtIndex)) ;
       
       pt   = mom.Pt();
       eta  = mom.Eta();
       phi  = mom.Phi() ;
       if(phi<0) phi+=TMath::TwoPi();
+      //only loop the particle at the same side of candidate
+      
+      //if at the same side has particle larger than candidate, then candidate can not be the leading, skip such events
+      if(TMath::Abs(phi-phiC)>TMath::PiOver2()) continue ;
+      
+      if(pt > ptC){
+        n         = -1;
+        nfrac     = -1;
+        coneptsum = -1;
+        isolated  = kFALSE;
+        if(bFillAOD){
+          if(reftracks){  
+            reftracks  ->Clear();
+            delete reftracks;
+          }
+          if(refclusters){
+            refclusters->Clear(); 
+            delete refclusters;
+          }
+        }
+        return ;
+      }
       
       //Check if there is any particle inside cone with pt larger than  fPtThreshold
-      rad = TMath::Sqrt((eta-etaC)*(eta-etaC)+ (phi-phiC)*(phi-phiC));
+
+      rad = Radius(etaC, phiC, eta, phi);
+      
       if(rad < fConeSize){
-       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++;
+        if(bFillAOD) {
+          nclusterrefs++;
+          if(nclusterrefs==1){
+            refclusters = new TObjArray(0);
+            //refclusters->SetName(Form("Clusters%s",aodArrayRefName.Data()));
+            TString tempo(aodArrayRefName)  ; 
+            tempo += "Clusters" ; 
+            refclusters->SetName(tempo);
+            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 fPtFraction*ptC<fPtThreshold then consider the fPtThreshold directly
+        if(fPtFraction*ptC<fPtThreshold) {
+          if(pt>fPtThreshold)    nfrac++ ;
+        }
+        else {
+          if(pt>fPtFraction*ptC) nfrac++; 
+        }
       }//in cone
     }// neutral particle loop
   }//neutrals
@@ -231,33 +275,30 @@ void  AliIsolationCut::MakeIsolationCut(TObjArray * const plCTS,  TObjArray * co
   //printf("Isolation Cut: in cone with: pT>pTthres %d, pT > pTfrac*pTcandidate %d \n",n,nfrac);
   
   //Add reference arrays to AOD when filling AODs only
-  if(fillAOD) {
+  if(bFillAOD) {
     if(refclusters)    pCandidate->AddObjArray(refclusters);
-    if(reftracks)      pCandidate->AddObjArray(reftracks);
+    if(reftracks)        pCandidate->AddObjArray(reftracks);
   }
-
   //Check isolation, depending on method.
   if( fICMethod == kPtThresIC){
     if(n==0) isolated = kTRUE ;
   }
   else if( fICMethod == kSumPtIC){
-    if(coneptsum < fPtThreshold)
+    if(coneptsum < fSumPtThreshold)
       isolated  =  kTRUE ;
   }
   else if( fICMethod == kPtFracIC){
     if(nfrac==0) isolated = kTRUE ;
   }
   else if( fICMethod == kSumPtFracIC){
-    if(coneptsum < fPtFraction*ptC)
-      isolated  =  kTRUE ;
+    //when the fPtFraction*ptC < fSumPtThreshold then consider the later case
+    if(fPtFraction*ptC < fSumPtThreshold  && coneptsum < fSumPtThreshold) isolated  =  kTRUE ;
+    if(fPtFraction*ptC > fSumPtThreshold  && coneptsum < fPtFraction*ptC) isolated  =  kTRUE ;
   }
-
-  //if(refclusters) delete refclusters;
-  //if(reftracks)   delete reftracks;
-
+  
 }
 
-//__________________________________________________________________
+//_____________________________________________________
 void AliIsolationCut::Print(const Option_t * opt) const
 {
   
@@ -267,11 +308,30 @@ void AliIsolationCut::Print(const Option_t * opt) const
   
   printf("**** Print %s %s **** \n", GetName(), GetTitle() ) ;
   
-  printf("IC method          =     %d\n", fICMethod) ; 
-  printf("Cone Size          =     %1.2f\n", fConeSize) ; 
+  printf("IC method          =     %d\n",    fICMethod   ) ; 
+  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("pT fraction        =     %3.1f\n", fPtFraction ) ;
+  printf("particle type in cone =  %d\n",    fPartInCone ) ;
   printf("    \n") ;
   
 } 
+
+//___________________________________________________________________________
+Float_t AliIsolationCut::Radius(const Float_t etaC, const Float_t phiC, 
+                                const Float_t eta , const Float_t phi) const
+{
+  // Calculate the distance to trigger from any particle
+
+  Float_t dEta = etaC-eta;
+  Float_t dPhi = phiC-phi;
+  
+  if(TMath::Abs(dPhi) >= TMath::Pi()) 
+    dPhi = TMath::TwoPi()-TMath::Abs(dPhi);
+  
+  return TMath::Sqrt( dEta*dEta + dPhi*dPhi );
+  
+}
+
+
+