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 "AliEMCALGeometry.h"
  39 #include "AliEMCALGeoParams.h"
  40 #include "AliCalorimeterUtils.h"
  41 #include "AliAODTrack.h"
  42 #include "AliVCluster.h"
  43 #include "AliCaloTrackReader.h"
  44 #include "AliMixedEvent.h"
  45 #include "AliCaloPID.h"
  46
  47 ClassImp(AliIsolationCut)
  48
  49 //____________________________________
  50 AliIsolationCut::AliIsolationCut() :
  51 TObject(),
  52 fConeSize(0.),
  53 fPtThreshold(0.),
  54 fSumPtThreshold(0.),
  55 fPtFraction(0.),
  56 fICMethod(0),
  57 fPartInCone(0),
  58 fDebug(-1),
  59 fFracIsThresh(1)
  60 {
  61   //default ctor
  62
  63   //Initialize parameters
  64   InitParameters();
  65
  66 }
  67
  68
  69 //__________________________________________________________________________________
  70 Float_t AliIsolationCut::GetCellDensity(const AliAODPWG4ParticleCorrelation * pCandidate,
  71                                         const AliCaloTrackReader * reader) const
  72 {
  73   // Get good cell density (number of active cells over all cells in cone)
  74
  75   Double_t coneCells    = 0.; //number of cells in cone with radius fConeSize
  76   Double_t coneCellsBad = 0.; //number of bad cells in cone with radius fConeSize
  77   Double_t cellDensity  = 1.;
  78
  79   Float_t phiC  = pCandidate->Phi() ;
  80   if(phiC<0) phiC+=TMath::TwoPi();
  81   Float_t etaC  = pCandidate->Eta() ;
  82
  83   if(pCandidate->GetDetector()=="EMCAL")
  84   {
  85     AliEMCALGeometry* eGeom = AliEMCALGeometry::GetInstance();
  86     AliCalorimeterUtils *cu = reader->GetCaloUtils();
  87
  88     Int_t absId = -999;
  89     if (eGeom->GetAbsCellIdFromEtaPhi(etaC,phiC,absId))
  90     {
  91       //Get absolute (col,row) of candidate
  92       Int_t iEta=-1, iPhi=-1, iRCU = -1;
  93       Int_t nSupMod = cu->GetModuleNumberCellIndexes(absId, pCandidate->GetDetector(), iEta, iPhi, iRCU);
  94
  95       Int_t colC = iEta;
  96       if (nSupMod % 2) colC =  AliEMCALGeoParams::fgkEMCALCols + iEta ;
  97       Int_t rowC = iPhi + AliEMCALGeoParams::fgkEMCALRows*int(nSupMod/2);
  98
  99       Int_t sqrSize = int(fConeSize/0.0143) ; // Size of cell in radians
 100       //loop on cells in a square of side fConeSize to check cells in cone
 101       for(Int_t icol = colC-sqrSize; icol < colC+sqrSize;icol++)
 102       {
 103         for(Int_t irow = rowC-sqrSize; irow < rowC+sqrSize; irow++)
 104         {
 105           if (Radius(colC, rowC, icol, irow) < sqrSize)
 106           {
 107             coneCells += 1.;
 108
 109             Int_t cellSM  = -999;
 110             Int_t cellEta = -999;
 111             Int_t cellPhi = -999;
 112             if(icol > AliEMCALGeoParams::fgkEMCALCols-1)
 113             {
 114               cellSM = 0+int(irow/AliEMCALGeoParams::fgkEMCALRows)*2;
 115               cellEta = icol-AliEMCALGeoParams::fgkEMCALCols;
 116               cellPhi = irow-AliEMCALGeoParams::fgkEMCALRows*int(cellSM/2);
 117             }
 118             if(icol < AliEMCALGeoParams::fgkEMCALCols)
 119             {
 120               cellSM = 1+int(irow/AliEMCALGeoParams::fgkEMCALRows)*2;
 121               cellEta = icol;
 122               cellPhi = irow-AliEMCALGeoParams::fgkEMCALRows*int(cellSM/2);
 123             }
 124
 125             //Count as bad "cells" out of EMCAL acceptance
 126             if(icol < 0 || icol > AliEMCALGeoParams::fgkEMCALCols*2 ||
 127                irow < 0 || irow > AliEMCALGeoParams::fgkEMCALRows*16./3) //5*nRows+1/3*nRows
 128             {
 129               coneCellsBad += 1.;
 130             }
 131             //Count as bad "cells" marked as bad in the DataBase
 132             else if (cu->GetEMCALChannelStatus(cellSM,cellEta,cellPhi)==1)
 133             {
 134               coneCellsBad += 1. ;
 135             }
 136           }
 137         }
 138       }//end of cells loop
 139     }
 140
 141     else if(fDebug>0) printf("cluster with bad (eta,phi) in EMCal for energy density calculation\n");
 142
 143     if (coneCells > 0.)
 144     {
 145       cellDensity = (coneCells-coneCellsBad)/coneCells;
 146       //printf("Energy density = %f\n", cellDensity);
 147     }
 148   }
 149
 150   return cellDensity;
 151
 152 }
 153
 154 //____________________________________________
 155 TString AliIsolationCut::GetICParametersList()
 156 {
 157   //Put data member values in string to keep in output container
 158
 159   TString parList ; //this will be list of parameters used for this analysis.
 160   const Int_t buffersize = 255;
 161   char onePar[buffersize] ;
 162
 163   snprintf(onePar,buffersize,"--- AliIsolationCut ---\n") ;
 164   parList+=onePar ;
 165   snprintf(onePar,buffersize,"fConeSize: (isolation cone size) %1.2f\n",fConeSize) ;
 166   parList+=onePar ;
 167   snprintf(onePar,buffersize,"fPtThreshold =%1.2f (isolation pt threshold) \n",fPtThreshold) ;
 168   parList+=onePar ;
 169   snprintf(onePar,buffersize,"fPtFraction=%1.2f (isolation pt threshold fraction ) \n",fPtFraction) ;
 170   parList+=onePar ;
 171   snprintf(onePar,buffersize,"fICMethod=%d (isolation cut case) \n",fICMethod) ;
 172   parList+=onePar ;
 173   snprintf(onePar,buffersize,"fPartInCone=%d \n",fPartInCone) ;
 174   parList+=onePar ;
 175  snprintf(onePar,buffersize,"fFracIsThresh=%i \n",fFracIsThresh) ;
 176   parList+=onePar ;
 177
 178   return parList;
 179 }
 180
 181 //____________________________________
 182 void AliIsolationCut::InitParameters()
 183 {
 184   //Initialize the parameters of the analysis.
 185
 186   fConeSize       = 0.4 ;
 187   fPtThreshold    = 1.  ;
 188   fSumPtThreshold = 0.5 ;
 189   fPtFraction     = 0.1 ;
 190   fPartInCone     = kOnlyCharged;
 191   fICMethod       = kSumPtFracIC; // 0 pt threshol method, 1 cone pt sum method
 192   fFracIsThresh   = 1;
 193 }
 194
 195 //________________________________________________________________________________
 196 void  AliIsolationCut::MakeIsolationCut(const TObjArray * plCTS,
 197                                         const TObjArray * plNe,
 198                                         const AliCaloTrackReader * reader,
 199                                         const AliCaloPID * pid,
 200                                         const Bool_t bFillAOD,
 201                                         AliAODPWG4ParticleCorrelation  *pCandidate,
 202                                         const TString & aodArrayRefName,
 203                                         Int_t   & n,
 204                                         Int_t   & nfrac,
 205                                         Float_t & coneptsum,
 206                                         Bool_t  & isolated) const
 207 {
 208   //Search in cone around a candidate particle if it is isolated
 209   Float_t ptC   = pCandidate->Pt() ;
 210   Float_t phiC  = pCandidate->Phi() ;
 211   if(phiC<0) phiC+=TMath::TwoPi();
 212   Float_t etaC  = pCandidate->Eta() ;
 213   Float_t pt    = -100. ;
 214   Float_t eta   = -100. ;
 215   Float_t phi   = -100. ;
 216   Float_t rad   = -100. ;
 217
 218   n         = 0 ;
 219   nfrac     = 0 ;
 220   coneptsum = 0.;
 221   isolated  = kFALSE;
 222
 223   if(fDebug>0)
 224   {
 225     printf("AliIsolationCut::MakeIsolationCut() - Cadidate pT %2.2f, eta %2.2f, phi %2.2f, cone %1.2f, thres %2.2f, Fill AOD? %d",
 226            pCandidate->Pt(), pCandidate->Eta(), pCandidate->Phi()*TMath::RadToDeg(), fConeSize,fPtThreshold,bFillAOD);
 227     if(plCTS) printf(", nTracks %d"  ,plCTS->GetEntriesFast());
 228     if(plNe)  printf(", nClusters %d",plNe ->GetEntriesFast());
 229
 230     printf("\n");
 231   }
 232
 233   //Initialize the array with refrences
 234   TObjArray * refclusters = 0x0;
 235   TObjArray * reftracks   = 0x0;
 236   Int_t ntrackrefs   = 0;
 237   Int_t nclusterrefs = 0;
 238
 239   //Check charged particles in cone.
 240   if(plCTS &&
 241      (fPartInCone==kOnlyCharged || fPartInCone==kNeutralAndCharged))
 242   {
 243     TVector3 p3;
 244     for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ )
 245     {
 246
 247       AliAODTrack* track = (AliAODTrack *)(plCTS->At(ipr)) ;
 248
 249       //Do not count the candidate (pion, conversion photon) or the daughters of the candidate
 250       if(track->GetID() == pCandidate->GetTrackLabel(0) || track->GetID() == pCandidate->GetTrackLabel(1) ||
 251          track->GetID() == pCandidate->GetTrackLabel(2) || track->GetID() == pCandidate->GetTrackLabel(3)   ) continue ;
 252
 253       p3.SetXYZ(track->Px(),track->Py(),track->Pz());
 254       pt   = p3.Pt();
 255       eta  = p3.Eta();
 256       phi  = p3.Phi() ;
 257       if(phi<0) phi+=TMath::TwoPi();
 258
 259       // Only loop the particle at the same side of candidate
 260       if(TMath::Abs(phi-phiC)>TMath::PiOver2()) continue ;
 261
 262       // If at the same side has particle larger than candidate,
 263       // then candidate can not be the leading, skip such events
 264       if(pt > ptC)
 265       {
 266         n         = -1;
 267         nfrac     = -1;
 268         coneptsum = -1;
 269         isolated  = kFALSE;
 270
 271         pCandidate->SetLeadingParticle(kFALSE);
 272
 273         if(bFillAOD && reftracks)
 274         {
 275           reftracks->Clear();
 276           delete reftracks;
 277         }
 278
 279         return ;
 280       }
 281
 282       //Check if there is any particle inside cone with pt larger than  fPtThreshold
 283
 284       rad = Radius(etaC, phiC, eta, phi);
 285
 286       if(fDebug>0)
 287         printf("\t track %d, pT %2.2f, eta %1.2f, phi %2.2f, R candidate %2.2f", ipr,pt,eta,phi,rad);
 288
 289       if(rad < fConeSize)
 290       {
 291         if(fDebug>0)  printf(" -  inside candidate cone");
 292
 293         if(bFillAOD)
 294         {
 295           ntrackrefs++;
 296           if(ntrackrefs == 1)
 297           {
 298             reftracks = new TObjArray(0);
 299             //reftracks->SetName(Form("Tracks%s",aodArrayRefName.Data()));
 300             TString tempo(aodArrayRefName)  ;
 301             tempo += "Tracks" ;
 302             reftracks->SetName(tempo);
 303             reftracks->SetOwner(kFALSE);
 304           }
 305           reftracks->Add(track);
 306         }
 307
 308
 309         coneptsum+=pt;
 310         if(pt > fPtThreshold )    n++;
 311         if(pt > fPtFraction*ptC ) nfrac++;
 312
 313       } // Inside cone
 314
 315       if(fDebug>0)  printf("\n");
 316
 317     }// charged particle loop
 318
 319
 320   }//Tracks
 321
 322
 323   //Check neutral particles in cone.
 324   if(plNe &&
 325      (fPartInCone==kOnlyNeutral || fPartInCone==kNeutralAndCharged))
 326   {
 327     TLorentzVector mom ;
 328
 329     for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ )
 330     {
 331       AliVCluster * calo = (AliVCluster *)(plNe->At(ipr)) ;
 332
 333       //Get the index where the cluster comes, to retrieve the corresponding vertex
 334       Int_t evtIndex = 0 ;
 335       if (reader->GetMixedEvent())
 336         evtIndex=reader->GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
 337
 338
 339       //Do not count the candidate (photon or pi0) or the daughters of the candidate
 340       if(calo->GetID() == pCandidate->GetCaloLabel(0) ||
 341          calo->GetID() == pCandidate->GetCaloLabel(1)   ) continue ;
 342
 343       //Skip matched clusters with tracks in case of neutral+charged analysis
 344       if( fPartInCone == kNeutralAndCharged &&
 345           pid->IsTrackMatched(calo,reader->GetCaloUtils(),reader->GetInputEvent()) ) continue ;
 346
 347       //Assume that come from vertex in straight line
 348       calo->GetMomentum(mom,reader->GetVertex(evtIndex)) ;
 349
 350       pt   = mom.Pt();
 351       eta  = mom.Eta();
 352       phi  = mom.Phi() ;
 353       if(phi<0) phi+=TMath::TwoPi();
 354
 355       // Only loop the particle at the same side of candidate
 356       if(TMath::Abs(phi-phiC)>TMath::PiOver2()) continue ;
 357
 358       // If at the same side has particle larger than candidate,
 359       // then candidate can not be the leading, skip such events
 360       if(pt > ptC)
 361       {
 362         n         = -1;
 363         nfrac     = -1;
 364         coneptsum = -1;
 365         isolated  = kFALSE;
 366
 367         pCandidate->SetLeadingParticle(kFALSE);
 368
 369         if(bFillAOD)
 370         {
 371           if(reftracks)
 372           {
 373             reftracks  ->Clear();
 374             delete reftracks;
 375           }
 376
 377           if(refclusters)
 378           {
 379             refclusters->Clear();
 380             delete refclusters;
 381           }
 382         }
 383         return ;
 384       }
 385
 386       //Check if there is any particle inside cone with pt larger than  fPtThreshold
 387
 388       rad = Radius(etaC, phiC, eta, phi);
 389
 390       if(fDebug>0)
 391         printf("\t cluster %d, pT %2.2f, eta %1.2f, phi %2.2f, R candidate %2.2f", ipr,pt,eta,phi,rad);
 392
 393       if(rad < fConeSize)
 394       {
 395         if(fDebug>0)  printf(" - inside candidate cone");
 396
 397         if(bFillAOD)
 398         {
 399           nclusterrefs++;
 400           if(nclusterrefs==1)
 401           {
 402             refclusters = new TObjArray(0);
 403             //refclusters->SetName(Form("Clusters%s",aodArrayRefName.Data()));
 404             TString tempo(aodArrayRefName)  ;
 405             tempo += "Clusters" ;
 406             refclusters->SetName(tempo);
 407             refclusters->SetOwner(kFALSE);
 408           }
 409           refclusters->Add(calo);
 410         }
 411
 412         coneptsum+=pt;
 413         if(pt > fPtThreshold )   n++;
 414         //if fPtFraction*ptC<fPtThreshold then consider the fPtThreshold directly
 415         if(fFracIsThresh){
 416           if( fPtFraction*ptC<fPtThreshold)
 417             {
 418               if(pt>fPtThreshold)    nfrac++ ;
 419             }
 420           else
 421             {
 422              if(pt>fPtFraction*ptC) nfrac++;
 423             }
 424         }
 425         else {
 426           if(pt>fPtFraction*ptC) nfrac++;
 427         }
 428
 429       }//in cone
 430
 431       if(fDebug>0)  printf("\n");
 432
 433     }// neutral particle loop
 434
 435   }//neutrals
 436
 437
 438   //Add reference arrays to AOD when filling AODs only
 439   if(bFillAOD)
 440   {
 441     if(refclusters)     pCandidate->AddObjArray(refclusters);
 442     if(reftracks)         pCandidate->AddObjArray(reftracks);
 443   }
 444
 445   //Check isolation, depending on selected isolation criteria
 446   if( fICMethod == kPtThresIC)
 447   {
 448     if(n==0) isolated = kTRUE ;
 449   }
 450   else if( fICMethod == kSumPtIC)
 451   {
 452     if(coneptsum < fSumPtThreshold)
 453       isolated  =  kTRUE ;
 454   }
 455   else if( fICMethod == kPtFracIC)
 456   {
 457     if(nfrac==0) isolated = kTRUE ;
 458   }
 459   else if( fICMethod == kSumPtFracIC)
 460   {
 461     //when the fPtFraction*ptC < fSumPtThreshold then consider the later case
 462     if(fFracIsThresh ){
 463       if(fPtFraction*ptC < fSumPtThreshold  && coneptsum < fSumPtThreshold) isolated  =  kTRUE ;
 464       if( fPtFraction*ptC > fSumPtThreshold  && coneptsum < fPtFraction*ptC) isolated  =  kTRUE ;
 465     }
 466     else
 467       {
 468         if(coneptsum < fPtFraction*ptC) isolated  =  kTRUE ;
 469       }
 470   }
 471  else if( fICMethod == kSumDensityIC)
 472   {
 473     // Get good cell density (number of active cells over all cells in cone)
 474     // and correct energy in cone
 475     Float_t cellDensity = GetCellDensity(pCandidate,reader);
 476     if(coneptsum < fSumPtThreshold*cellDensity)
 477       isolated = kTRUE;
 478   }
 479
 480 }
 481
 482 //_____________________________________________________
 483 void AliIsolationCut::Print(const Option_t * opt) const
 484 {
 485
 486   //Print some relevant parameters set for the analysis
 487   if(! opt)
 488     return;
 489
 490   printf("**** Print %s %s **** \n", GetName(), GetTitle() ) ;
 491
 492   printf("IC method          =     %d\n",    fICMethod   ) ;
 493   printf("Cone Size          =     %1.2f\n", fConeSize   ) ;
 494   printf("pT threshold       =     %2.1f\n", fPtThreshold) ;
 495   printf("pT fraction        =     %3.1f\n", fPtFraction ) ;
 496   printf("particle type in cone =  %d\n",    fPartInCone ) ;
 497   printf("using fraction for high pt leading instead of frac ? %i\n",fFracIsThresh);
 498   printf("    \n") ;
 499
 500 }
 501
 502 //___________________________________________________________________________
 503 Float_t AliIsolationCut::Radius(const Float_t etaC, const Float_t phiC,
 504                                 const Float_t eta , const Float_t phi) const
 505 {
 506   // Calculate the distance to trigger from any particle
 507
 508   Float_t dEta = etaC-eta;
 509   Float_t dPhi = phiC-phi;
 510
 511   if(TMath::Abs(dPhi) >= TMath::Pi())
 512     dPhi = TMath::TwoPi()-TMath::Abs(dPhi);
 513
 514   return TMath::Sqrt( dEta*dEta + dPhi*dPhi );
 515
 516 }
 517
 518
 519