PWG4/PartCorrBase/AliIsolationCut.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 // Class containing methods for the isolation cut.
  19 // An AOD candidate (AliAODPWG4ParticleCorrelation type)
  20 // is passed. Look in a cone around the candidate and study
  21 // the hadronic activity inside to decide if the candidate is isolated
  22 //
  23 //
  24 //*-- Author: Gustavo Conesa (LNF-INFN)
  25
  26 //-Yaxian Mao (add the possibility for different IC method with different pt range, 01/10/2010)
  27 //-Yaxian Mao (check the candidate particle is the leading particle or not at the same hemishere)
  28
  29 //////////////////////////////////////////////////////////////////////////////
  30
  31
  32 // --- ROOT system ---
  33 #include <TLorentzVector.h>
  34 #include <TObjArray.h>
  35
  36 // --- AliRoot system ---
  37 #include "AliIsolationCut.h"
  38 #include "AliAODPWG4ParticleCorrelation.h"
  39 #include "AliAODTrack.h"
  40 #include "AliVCluster.h"
  41 #include "AliCaloTrackReader.h"
  42 #include "AliMixedEvent.h"
  43
  44 ClassImp(AliIsolationCut)
  45
  46 //____________________________________
  47 AliIsolationCut::AliIsolationCut() :
  48 TObject(),
  49 fConeSize(0.),
  50 fPtThreshold(0.),
  51 fSumPtThreshold(0.),
  52 fPtFraction(0.),
  53 fICMethod(0),
  54 fPartInCone(0)
  55
  56 {
  57   //default ctor
  58
  59   //Initialize parameters
  60   InitParameters();
  61
  62 }
  63
  64 //____________________________________________
  65 TString AliIsolationCut::GetICParametersList()
  66 {
  67   //Put data member values in string to keep in output container
  68
  69   TString parList ; //this will be list of parameters used for this analysis.
  70   const Int_t buffersize = 255;
  71   char onePar[buffersize] ;
  72
  73   snprintf(onePar,buffersize,"--- AliIsolationCut ---\n") ;
  74   parList+=onePar ;
  75   snprintf(onePar,buffersize,"fConeSize: (isolation cone size) %1.2f\n",fConeSize) ;
  76   parList+=onePar ;
  77   snprintf(onePar,buffersize,"fPtThreshold =%1.2f (isolation pt threshold) \n",fPtThreshold) ;
  78   parList+=onePar ;
  79   snprintf(onePar,buffersize,"fPtFraction=%1.2f (isolation pt threshold fraction ) \n",fPtFraction) ;
  80   parList+=onePar ;
  81   snprintf(onePar,buffersize,"fICMethod=%d (isolation cut case) \n",fICMethod) ;
  82   parList+=onePar ;
  83   snprintf(onePar,buffersize,"fPartInCone=%d \n",fPartInCone) ;
  84   parList+=onePar ;
  85
  86   return parList;
  87 }
  88
  89 //____________________________________________________________________________
  90 void AliIsolationCut::InitParameters()
  91 {
  92   //Initialize the parameters of the analysis.
  93
  94   fConeSize       = 0.4 ;
  95   fPtThreshold    = 1.  ;
  96   fSumPtThreshold = 0.5 ;
  97   fPtFraction     = 0.1 ;
  98   fPartInCone     = kOnlyCharged;
  99   fICMethod       = kSumPtFracIC; // 0 pt threshol method, 1 cone pt sum method
 100
 101 }
 102
 103 //________________________________________________________________________________
 104 void  AliIsolationCut::MakeIsolationCut(const TObjArray * plCTS,
 105                                         const TObjArray * plNe,
 106                                         const AliCaloTrackReader * reader,
 107                                         const Bool_t bFillAOD,
 108                                         AliAODPWG4ParticleCorrelation  *pCandidate,
 109                                         const TString & aodArrayRefName,
 110                                         Int_t & n,
 111                                         Int_t & nfrac,
 112                                         Float_t &coneptsum,
 113                                         Bool_t  &isolated) const
 114 {
 115   //Search in cone around a candidate particle if it is isolated
 116   Float_t phiC  = pCandidate->Phi() ;
 117   if(phiC<0) phiC+=TMath::TwoPi();
 118   Float_t etaC  = pCandidate->Eta() ;
 119   Float_t ptC   = pCandidate->Pt() ;
 120   Float_t pt    = -100. ;
 121   Float_t eta   = -100. ;
 122   Float_t phi   = -100. ;
 123   Float_t rad   = -100. ;
 124
 125   n         = 0 ;
 126   nfrac     = 0 ;
 127   coneptsum = 0.;
 128   isolated  = kFALSE;
 129
 130   //Initialize the array with refrences
 131   TObjArray * refclusters = 0x0;
 132   TObjArray * reftracks   = 0x0;
 133   Int_t ntrackrefs   = 0;
 134   Int_t nclusterrefs = 0;
 135   //Check charged particles in cone.
 136   if(plCTS && (fPartInCone==kOnlyCharged || fPartInCone==kNeutralAndCharged)){
 137     TVector3 p3;
 138     for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ ){
 139       AliAODTrack* track = (AliAODTrack *)(plCTS->At(ipr)) ;
 140       //Do not count the candidate (pion, conversion photon) or the daughters of the candidate
 141       if(track->GetID() == pCandidate->GetTrackLabel(0) || track->GetID() == pCandidate->GetTrackLabel(1)
 142          || track->GetID() == pCandidate->GetTrackLabel(2) || track->GetID() == pCandidate->GetTrackLabel(3)
 143          ) continue ;
 144       p3.SetXYZ(track->Px(),track->Py(),track->Pz());
 145       pt   = p3.Pt();
 146       eta  = p3.Eta();
 147       phi  = p3.Phi() ;
 148       if(phi<0) phi+=TMath::TwoPi();
 149
 150       //only loop the particle at the same side of candidate
 151       if(TMath::Abs(phi-phiC)>TMath::PiOver2()) continue ;
 152       //if at the same side has particle larger than candidate, then candidate can not be the leading, skip such events
 153       if(pt > ptC){
 154         n         = -1;
 155         nfrac     = -1;
 156         coneptsum = -1;
 157         isolated  = kFALSE;
 158         if(bFillAOD && reftracks) {
 159           reftracks->Clear();
 160           delete reftracks;
 161         }
 162         return ;
 163       }
 164       //Check if there is any particle inside cone with pt larger than  fPtThreshold
 165       rad = TMath::Sqrt((eta-etaC)*(eta-etaC)+ (phi-phiC)*(phi-phiC));
 166
 167       if(rad < fConeSize){
 168         if(bFillAOD) {
 169           ntrackrefs++;
 170           if(ntrackrefs == 1){
 171             reftracks = new TObjArray(0);
 172             //reftracks->SetName(Form("Tracks%s",aodArrayRefName.Data()));
 173             TString tempo(aodArrayRefName)  ;
 174             tempo += "Tracks" ;
 175             reftracks->SetName(tempo);
 176             reftracks->SetOwner(kFALSE);
 177           }
 178           reftracks->Add(track);
 179         }
 180         //printf("charged in isolation cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
 181         coneptsum+=pt;
 182         if(pt > fPtThreshold )    n++;
 183         if(pt > fPtFraction*ptC ) nfrac++;
 184       } // Inside cone
 185     }// charged particle loop
 186   }//Tracks
 187
 188   //Check neutral particles in cone.
 189   if(plNe && (fPartInCone==kOnlyNeutral || fPartInCone==kNeutralAndCharged)){
 190
 191
 192     TLorentzVector mom ;
 193     for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ ){
 194       AliVCluster * calo = (AliVCluster *)(plNe->At(ipr)) ;
 195
 196       //Get the index where the cluster comes, to retrieve the corresponding vertex
 197       Int_t evtIndex = 0 ;
 198       if (reader->GetMixedEvent()) {
 199         evtIndex=reader->GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
 200       }
 201
 202       //Do not count the candidate (photon or pi0) or the daughters of the candidate
 203       if(calo->GetID() == pCandidate->GetCaloLabel(0) || calo->GetID() == pCandidate->GetCaloLabel(1)) continue ;      //Skip matched clusters with tracks
 204
 205       if(calo->GetNTracksMatched() > 0) continue ;
 206
 207       calo->GetMomentum(mom,reader->GetVertex(evtIndex)) ;//Assume that come from vertex in straight line
 208
 209       pt   = mom.Pt();
 210       eta  = mom.Eta();
 211       phi  = mom.Phi() ;
 212       if(phi<0) phi+=TMath::TwoPi();
 213       //only loop the particle at the same side of candidate
 214
 215       if(TMath::Abs(phi-phiC)>TMath::PiOver2()) continue ;
 216       //if at the same side has particle larger than candidate, then candidate can not be the leading, skip such events
 217       if(pt > ptC){
 218         n         = -1;
 219         nfrac     = -1;
 220         coneptsum = -1;
 221         isolated  = kFALSE;
 222         if(bFillAOD){
 223           if(reftracks){
 224             reftracks  ->Clear();
 225             delete reftracks;
 226           }
 227           if(refclusters){
 228             refclusters->Clear();
 229             delete refclusters;
 230           }
 231         }
 232         return ;
 233       }
 234
 235       //Check if there is any particle inside cone with pt larger than  fPtThreshold
 236       rad = TMath::Sqrt((eta-etaC)*(eta-etaC)+ (phi-phiC)*(phi-phiC));
 237       if(rad < fConeSize){
 238         if(bFillAOD) {
 239           nclusterrefs++;
 240           if(nclusterrefs==1){
 241             refclusters = new TObjArray(0);
 242             //refclusters->SetName(Form("Clusters%s",aodArrayRefName.Data()));
 243             TString tempo(aodArrayRefName)  ;
 244             tempo += "Clusters" ;
 245             refclusters->SetName(tempo);
 246             refclusters->SetOwner(kFALSE);
 247           }
 248           refclusters->Add(calo);
 249         }
 250         //printf("neutral in isolation cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
 251         coneptsum+=pt;
 252         if(pt > fPtThreshold )     n++;
 253         //if fPtFraction*ptC<fPtThreshold then consider the fPtThreshold directly
 254         if(fPtFraction*ptC<fPtThreshold) {
 255           if(pt>fPtThreshold)    nfrac++ ;
 256         }
 257         else {
 258           if(pt>fPtFraction*ptC) nfrac++;
 259         }
 260       }//in cone
 261     }// neutral particle loop
 262   }//neutrals
 263
 264   //printf("Isolation Cut: in cone with: pT>pTthres %d, pT > pTfrac*pTcandidate %d \n",n,nfrac);
 265
 266   //Add reference arrays to AOD when filling AODs only
 267   if(bFillAOD) {
 268     if(refclusters)     pCandidate->AddObjArray(refclusters);
 269     if(reftracks)         pCandidate->AddObjArray(reftracks);
 270   }
 271   //Check isolation, depending on method.
 272   if( fICMethod == kPtThresIC){
 273     if(n==0) isolated = kTRUE ;
 274   }
 275   else if( fICMethod == kSumPtIC){
 276     if(coneptsum < fSumPtThreshold)
 277       isolated  =  kTRUE ;
 278   }
 279   else if( fICMethod == kPtFracIC){
 280     if(nfrac==0) isolated = kTRUE ;
 281   }
 282   else if( fICMethod == kSumPtFracIC){
 283     //when the fPtFraction*ptC < fSumPtThreshold then consider the later case
 284     if(fPtFraction*ptC < fSumPtThreshold  && coneptsum < fSumPtThreshold) isolated  =  kTRUE ;
 285     if(fPtFraction*ptC > fSumPtThreshold  && coneptsum < fPtFraction*ptC) isolated  =  kTRUE ;
 286   }
 287
 288 }
 289
 290 //_____________________________________________________
 291 void AliIsolationCut::Print(const Option_t * opt) const
 292 {
 293
 294   //Print some relevant parameters set for the analysis
 295   if(! opt)
 296     return;
 297
 298   printf("**** Print %s %s **** \n", GetName(), GetTitle() ) ;
 299
 300   printf("IC method          =     %d\n",    fICMethod   ) ;
 301   printf("Cone Size          =     %1.2f\n", fConeSize   ) ;
 302   printf("pT threshold       =     %2.1f\n", fPtThreshold) ;
 303   printf("pT fraction        =     %3.1f\n", fPtFraction ) ;
 304   printf("particle type in cone =  %d\n",    fPartInCone ) ;
 305   printf("    \n") ;
 306
 307 }