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
  16 //_________________________________________________________________________
  17 // Class containing methods for the isolation cut.
  18 // An AOD candidate (AliAODPWG4ParticleCorrelation type)
  19 // is passed. Look in a cone around the candidate and study
  20 // the hadronic activity inside to decide if the candidate is isolated
  21 //
  22 //
  23 //*-- Author: Gustavo Conesa (LNF-INFN)
  24
  25 //-Yaxian Mao (add the possibility for different IC method with different pt range, 01/10/2010)
  26 //-Yaxian Mao (check the candidate particle is the leading particle or not at the same hemishere)
  27
  28 //////////////////////////////////////////////////////////////////////////////
  29
  30
  31 // --- ROOT system ---
  32 #include <TLorentzVector.h>
  33 #include <TObjArray.h>
  34
  35 // --- AliRoot system ---
  36 #include "AliIsolationCut.h"
  37 #include "AliAODPWG4ParticleCorrelation.h"
  38 #include "AliAODTrack.h"
  39 #include "AliVCluster.h"
  40 #include "AliCaloTrackReader.h"
  41 #include "AliMixedEvent.h"
  42 #include "AliCaloPID.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   //default ctor
  57
  58   //Initialize parameters
  59   InitParameters();
  60
  61 }
  62
  63 //____________________________________________
  64 TString AliIsolationCut::GetICParametersList()
  65 {
  66   //Put data member values in string to keep in output container
  67
  68   TString parList ; //this will be list of parameters used for this analysis.
  69   const Int_t buffersize = 255;
  70   char onePar[buffersize] ;
  71
  72   snprintf(onePar,buffersize,"--- AliIsolationCut ---\n") ;
  73   parList+=onePar ;
  74   snprintf(onePar,buffersize,"fConeSize: (isolation cone size) %1.2f\n",fConeSize) ;
  75   parList+=onePar ;
  76   snprintf(onePar,buffersize,"fPtThreshold =%1.2f (isolation pt threshold) \n",fPtThreshold) ;
  77   parList+=onePar ;
  78   snprintf(onePar,buffersize,"fPtFraction=%1.2f (isolation pt threshold fraction ) \n",fPtFraction) ;
  79   parList+=onePar ;
  80   snprintf(onePar,buffersize,"fICMethod=%d (isolation cut case) \n",fICMethod) ;
  81   parList+=onePar ;
  82   snprintf(onePar,buffersize,"fPartInCone=%d \n",fPartInCone) ;
  83   parList+=onePar ;
  84
  85   return parList;
  86 }
  87
  88 //____________________________________
  89 void AliIsolationCut::InitParameters()
  90 {
  91   //Initialize the parameters of the analysis.
  92
  93   fConeSize       = 0.4 ;
  94   fPtThreshold    = 1.  ;
  95   fSumPtThreshold = 0.5 ;
  96   fPtFraction     = 0.1 ;
  97   fPartInCone     = kOnlyCharged;
  98   fICMethod       = kSumPtFracIC; // 0 pt threshol method, 1 cone pt sum method
  99
 100 }
 101
 102 //________________________________________________________________________________
 103 void  AliIsolationCut::MakeIsolationCut(const TObjArray * plCTS,
 104                                         const TObjArray * plNe,
 105                                         const AliCaloTrackReader * reader,
 106                                         const AliCaloPID * pid,
 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
 136   //Check charged particles in cone.
 137   if(plCTS &&
 138      (fPartInCone==kOnlyCharged || fPartInCone==kNeutralAndCharged))
 139   {
 140     TVector3 p3;
 141     for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ )
 142     {
 143
 144       AliAODTrack* track = (AliAODTrack *)(plCTS->At(ipr)) ;
 145
 146       //Do not count the candidate (pion, conversion photon) or the daughters of the candidate
 147       if(track->GetID() == pCandidate->GetTrackLabel(0) || track->GetID() == pCandidate->GetTrackLabel(1) ||
 148          track->GetID() == pCandidate->GetTrackLabel(2) || track->GetID() == pCandidate->GetTrackLabel(3)   ) continue ;
 149
 150       p3.SetXYZ(track->Px(),track->Py(),track->Pz());
 151       pt   = p3.Pt();
 152       eta  = p3.Eta();
 153       phi  = p3.Phi() ;
 154       if(phi<0) phi+=TMath::TwoPi();
 155
 156       // Only loop the particle at the same side of candidate
 157       if(TMath::Abs(phi-phiC)>TMath::PiOver2()) continue ;
 158
 159       // If at the same side has particle larger than candidate,
 160       // then candidate can not be the leading, skip such events
 161       if(pt > ptC)
 162       {
 163         n         = -1;
 164         nfrac     = -1;
 165         coneptsum = -1;
 166         isolated  = kFALSE;
 167         if(bFillAOD && reftracks)
 168         {
 169           reftracks->Clear();
 170           delete reftracks;
 171         }
 172
 173         return ;
 174       }
 175
 176       //Check if there is any particle inside cone with pt larger than  fPtThreshold
 177
 178       rad = Radius(etaC, phiC, eta, phi);
 179
 180       if(rad < fConeSize)
 181       {
 182         if(bFillAOD)
 183         {
 184           ntrackrefs++;
 185           if(ntrackrefs == 1)
 186           {
 187             reftracks = new TObjArray(0);
 188             //reftracks->SetName(Form("Tracks%s",aodArrayRefName.Data()));
 189             TString tempo(aodArrayRefName)  ;
 190             tempo += "Tracks" ;
 191             reftracks->SetName(tempo);
 192             reftracks->SetOwner(kFALSE);
 193           }
 194           reftracks->Add(track);
 195         }
 196
 197         //printf("charged in isolation cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
 198         coneptsum+=pt;
 199         if(pt > fPtThreshold )    n++;
 200         if(pt > fPtFraction*ptC ) nfrac++;
 201
 202       } // Inside cone
 203     }// charged particle loop
 204   }//Tracks
 205
 206   //Check neutral particles in cone.
 207   if(plNe &&
 208      (fPartInCone==kOnlyNeutral || fPartInCone==kNeutralAndCharged))
 209   {
 210     TLorentzVector mom ;
 211
 212     for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ )
 213     {
 214       AliVCluster * calo = (AliVCluster *)(plNe->At(ipr)) ;
 215
 216       //Get the index where the cluster comes, to retrieve the corresponding vertex
 217       Int_t evtIndex = 0 ;
 218       if (reader->GetMixedEvent())
 219         evtIndex=reader->GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
 220
 221
 222       //Do not count the candidate (photon or pi0) or the daughters of the candidate
 223       if(calo->GetID() == pCandidate->GetCaloLabel(0) ||
 224          calo->GetID() == pCandidate->GetCaloLabel(1)   ) continue ;
 225
 226       //Skip matched clusters with tracks
 227       if( pid->IsTrackMatched(calo,reader->GetCaloUtils(),reader->GetInputEvent()) ) continue ;
 228
 229       //Assume that come from vertex in straight line
 230       calo->GetMomentum(mom,reader->GetVertex(evtIndex)) ;
 231
 232       pt   = mom.Pt();
 233       eta  = mom.Eta();
 234       phi  = mom.Phi() ;
 235       if(phi<0) phi+=TMath::TwoPi();
 236
 237       // Only loop the particle at the same side of candidate
 238       if(TMath::Abs(phi-phiC)>TMath::PiOver2()) continue ;
 239
 240       // If at the same side has particle larger than candidate,
 241       // then candidate can not be the leading, skip such events
 242       if(pt > ptC)
 243       {
 244         n         = -1;
 245         nfrac     = -1;
 246         coneptsum = -1;
 247         isolated  = kFALSE;
 248
 249         if(bFillAOD)
 250         {
 251           if(reftracks)
 252           {
 253             reftracks  ->Clear();
 254             delete reftracks;
 255           }
 256
 257           if(refclusters)
 258           {
 259             refclusters->Clear();
 260             delete refclusters;
 261           }
 262         }
 263         return ;
 264       }
 265
 266       //Check if there is any particle inside cone with pt larger than  fPtThreshold
 267
 268       rad = Radius(etaC, phiC, eta, phi);
 269
 270       if(rad < fConeSize)
 271       {
 272         if(bFillAOD)
 273         {
 274           nclusterrefs++;
 275           if(nclusterrefs==1)
 276           {
 277             refclusters = new TObjArray(0);
 278             //refclusters->SetName(Form("Clusters%s",aodArrayRefName.Data()));
 279             TString tempo(aodArrayRefName)  ;
 280             tempo += "Clusters" ;
 281             refclusters->SetName(tempo);
 282             refclusters->SetOwner(kFALSE);
 283           }
 284           refclusters->Add(calo);
 285         }
 286
 287         //printf("neutral in isolation cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
 288         coneptsum+=pt;
 289         if(pt > fPtThreshold )   n++;
 290         //if fPtFraction*ptC<fPtThreshold then consider the fPtThreshold directly
 291         if(fPtFraction*ptC<fPtThreshold)
 292         {
 293           if(pt>fPtThreshold)    nfrac++ ;
 294         }
 295         else
 296         {
 297           if(pt>fPtFraction*ptC) nfrac++;
 298         }
 299
 300       }//in cone
 301
 302     }// neutral particle loop
 303
 304   }//neutrals
 305
 306   //printf("Isolation Cut: in cone with: pT>pTthres %d, pT > pTfrac*pTcandidate %d \n",n,nfrac);
 307
 308   //Add reference arrays to AOD when filling AODs only
 309   if(bFillAOD)
 310   {
 311     if(refclusters)     pCandidate->AddObjArray(refclusters);
 312     if(reftracks)         pCandidate->AddObjArray(reftracks);
 313   }
 314
 315   //Check isolation, depending on selected isolation criteria
 316   if( fICMethod == kPtThresIC)
 317   {
 318     if(n==0) isolated = kTRUE ;
 319   }
 320   else if( fICMethod == kSumPtIC)
 321   {
 322     if(coneptsum < fSumPtThreshold)
 323       isolated  =  kTRUE ;
 324   }
 325   else if( fICMethod == kPtFracIC)
 326   {
 327     if(nfrac==0) isolated = kTRUE ;
 328   }
 329   else if( fICMethod == kSumPtFracIC)
 330   {
 331     //when the fPtFraction*ptC < fSumPtThreshold then consider the later case
 332     if(fPtFraction*ptC < fSumPtThreshold  && coneptsum < fSumPtThreshold) isolated  =  kTRUE ;
 333     if(fPtFraction*ptC > fSumPtThreshold  && coneptsum < fPtFraction*ptC) isolated  =  kTRUE ;
 334   }
 335
 336 }
 337
 338 //_____________________________________________________
 339 void AliIsolationCut::Print(const Option_t * opt) const
 340 {
 341
 342   //Print some relevant parameters set for the analysis
 343   if(! opt)
 344     return;
 345
 346   printf("**** Print %s %s **** \n", GetName(), GetTitle() ) ;
 347
 348   printf("IC method          =     %d\n",    fICMethod   ) ;
 349   printf("Cone Size          =     %1.2f\n", fConeSize   ) ;
 350   printf("pT threshold       =     %2.1f\n", fPtThreshold) ;
 351   printf("pT fraction        =     %3.1f\n", fPtFraction ) ;
 352   printf("particle type in cone =  %d\n",    fPartInCone ) ;
 353   printf("    \n") ;
 354
 355 }
 356
 357 //___________________________________________________________________________
 358 Float_t AliIsolationCut::Radius(const Float_t etaC, const Float_t phiC,
 359                                 const Float_t eta , const Float_t phi) const
 360 {
 361   // Calculate the distance to trigger from any particle
 362
 363   Float_t dEta = etaC-eta;
 364   Float_t dPhi = phiC-phi;
 365
 366   if(TMath::Abs(dPhi) >= TMath::Pi())
 367     dPhi = TMath::TwoPi()-TMath::Abs(dPhi);
 368
 369   return TMath::Sqrt( dEta*dEta + dPhi*dPhi );
 370
 371 }
 372
 373
 374