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