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