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