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