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