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 fSumPtThresholdMax(10000.),
  57 fPtFraction(0.),
  58 fICMethod(0),
  59 fPartInCone(0),
  60 fDebug(-1),
  61 fFracIsThresh(1)
  62 {
  63   //default ctor
  64
  65   //Initialize parameters
  66   InitParameters();
  67
  68 }
  69
  70 //_________________________________________________________________________________________________________________________________
  71 void AliIsolationCut::CalculateUEBandClusterNormalization(AliCaloTrackReader * /*reader*/, Float_t   etaC, Float_t /*phiC*/,
  72                                                            Float_t   phiUEptsumCluster,     Float_t   etaUEptsumCluster,
  73                                                            Float_t & phiUEptsumClusterNorm, Float_t & etaUEptsumClusterNorm,
  74                                                            Float_t & excessFracEta,         Float_t & excessFracPhi         ) const
  75 {
  76   // Normalize cluster background band
  77
  78   Float_t coneA     = fConeSize*fConeSize*TMath::Pi(); // A = pi R^2, isolation cone area
  79
  80   //Careful here if EMCal limits changed .. 2010 (4 SM) to 2011-12 (10 SM), for the moment consider 100 deg in phi
  81   Float_t emcEtaSize = 0.7*2; // TO FIX
  82   Float_t emcPhiSize = TMath::DegToRad()*100.; // TO FIX
  83
  84   /* //Catherine code
  85    if(((((2*fConeSize*emcPhiSize)-coneA))*phiBandBadCellsCoeff)!=0)phiUEptsumClusterNorm = phiUEptsumCluster*(coneA*coneBadCellsCoeff / (((2*fConeSize*emcPhiSize)-coneA))*phiBandBadCellsCoeff); // pi * R^2 / (2 R * 2 100 deg) -  trigger cone
  86    if(((((2*(fConeSize-excess)*emcPhiSize)-(coneA-excessFracEta))*etaBandBadCellsCoeff))!=0)phiUEptsumClusterNorm = phiUEptsumCluster*(coneA *coneBadCellsCoeff/ (((2*(fConeSize-excess)*emcPhiSize)-(coneA/excessFracEta))*etaBandBadCellsCoeff));
  87    if(((2*(fConeSize-excess)*emcEtaSize)-(coneA-excessFracPhi))*phiBandBadCellsCoeff!=0) etaUEptsumClusterNorm = etaUEptsumCluster*(coneA*coneBadCellsCoeff / (((2*(fConeSize-excess)*emcEtaSize)-(coneA/excessFracPhi))*phiBandBadCellsCoeff));
  88    */
  89
  90   if((2*fConeSize*emcPhiSize-coneA)!=0) phiUEptsumClusterNorm = phiUEptsumCluster*(coneA / (((2*fConeSize*emcPhiSize)-coneA))); // pi * R^2 / (2 R * 2 100 deg) -  trigger cone
  91   if((2*fConeSize*emcEtaSize-coneA)!=0) etaUEptsumClusterNorm = etaUEptsumCluster*(coneA / (((2*fConeSize*emcEtaSize)-coneA))); // pi * R^2 / (2 R * 2*0.7)  -  trigger cone
  92
  93   //out of eta acceptance
  94   excessFracEta = 1;
  95   excessFracPhi = 1;
  96
  97   if(TMath::Abs(etaC)+fConeSize > emcEtaSize/2.)
  98   {
  99     Float_t excess = TMath::Abs(etaC) + fConeSize - emcEtaSize/2.;
 100     excessFracEta  = CalculateExcessAreaFraction(excess);
 101
 102     if ( excessFracEta != 0) coneA /=  excessFracEta;
 103
 104     //UE band is also out of acceptance, need to estimate corrected area
 105     if(((2*fConeSize-excess)*emcPhiSize-coneA) != 0 ) phiUEptsumClusterNorm = phiUEptsumCluster*(coneA / ((((2*fConeSize-excess)*emcPhiSize)-coneA)));
 106     if(( 2*fConeSize        *emcEtaSize-coneA) != 0 ) etaUEptsumClusterNorm = etaUEptsumCluster*(coneA / ((( 2*fConeSize        *emcEtaSize)-coneA)));
 107   }
 108
 109 }
 110
 111 //________________________________________________________________________________________________________________________________
 112 void AliIsolationCut::CalculateUEBandTrackNormalization  (AliCaloTrackReader * reader,    Float_t   etaC, Float_t /*phiC*/,
 113                                                           Float_t   phiUEptsumTrack,      Float_t   etaUEptsumTrack,
 114                                                           Float_t & phiUEptsumTrackNorm,  Float_t & etaUEptsumTrackNorm,
 115                                                           Float_t & excessFracEta,        Float_t & excessFracPhi          ) const
 116 {
 117   // Normalize track background band
 118
 119   Float_t coneA     = fConeSize*fConeSize*TMath::Pi(); // A = pi R^2, isolation cone area
 120
 121   // Get the cut used for the TPC tracks in the reader, +-0.8, +-0.9 ...
 122   // Only valid in simple fidutial cut case and if the cut is applied, careful!
 123   Float_t tpcEtaSize = reader->GetFiducialCut()->GetCTSFidCutMaxEtaArray()->At(0) -
 124   reader->GetFiducialCut()->GetCTSFidCutMinEtaArray()->At(0) ;
 125   Float_t tpcPhiSize = TMath::TwoPi();
 126
 127   /*//Catherine code
 128    //phiUEptsumTrackNorm = phiUEptsumTrack*(coneA*coneBadCellsCoeff / (((2*fConeSize*tpcPhiSize)-coneA))*phiBandBadCellsCoeff); // pi * R^2 / (2 R * 2 pi) -  trigger cone
 129    //etaUEptsumTrackNorm = etaUEptsumTrack*(coneA*coneBadCellsCoeff / (((2*fConeSize*tpcEtaSize)-coneA))*etaBandBadCellsCoeff); // pi * R^2 / (2 R * 1.6)  -  trigger cone
 130    if((2*fConeSize*tpcPhiSize-coneA)!=0)phiUEptsumTrackNorm = phiUEptsumTrack*(coneA / (((2*fConeSize*tpcPhiSize)-coneA))); // pi * R^2 / (2 R * 2 pi) -  trigger cone
 131    if((2*fConeSize*tpcEtaSize-coneA)!=0)etaUEptsumTrackNorm = etaUEptsumTrack*(coneA / (((2*fConeSize*tpcEtaSize)-coneA))); // pi * R^2 / (2 R * 1.6)  -  trigger cone
 132    if((2*(fConeSize-excess)*tpcPhiSize)-(coneA-excessFracEta)!=0)phiUEptsumTrackNorm = phiUEptsumTrack*(coneA / (((2*(fConeSize-excess)*tpcPhiSize)-(coneA/excessFracEta))));
 133    */ //end Catherine code
 134
 135   //correct out of eta acceptance
 136   excessFracEta = 1;
 137   excessFracPhi = 1;
 138
 139   if((2*fConeSize*tpcPhiSize-coneA)!=0) phiUEptsumTrackNorm = phiUEptsumTrack*(coneA / (((2*fConeSize*tpcPhiSize)-coneA))); // pi * R^2 / (2 R * 2 pi) -  trigger cone
 140   if((2*fConeSize*tpcEtaSize-coneA)!=0) etaUEptsumTrackNorm = etaUEptsumTrack*(coneA / (((2*fConeSize*tpcEtaSize)-coneA))); // pi * R^2 / (2 R * 1.6)  -  trigger cone
 141
 142   if(TMath::Abs(etaC)+fConeSize > tpcEtaSize/2.)
 143   {
 144     Float_t excess = TMath::Abs(etaC) + fConeSize - tpcEtaSize/2.;
 145     excessFracEta  = CalculateExcessAreaFraction(excess);
 146     if (excessFracEta != 0) coneA /=  excessFracEta;
 147
 148     //UE band is also out of acceptance, need to estimate corrected area
 149     if(((2*fConeSize-excess)*tpcPhiSize - coneA) !=0 ) phiUEptsumTrackNorm = phiUEptsumTrack*(coneA / ((((2*fConeSize-excess)*tpcPhiSize)-coneA)));
 150     if(( 2*fConeSize        *tpcEtaSize - coneA) !=0 ) etaUEptsumTrackNorm = etaUEptsumTrack*(coneA / ((( 2*fConeSize        *tpcEtaSize)-coneA)));
 151   }
 152
 153 }
 154
 155 //________________________________________________________________________
 156 Float_t AliIsolationCut::CalculateExcessAreaFraction(Float_t excess) const
 157 {
 158   // Area of a circunference segment segment 1/2 R^2 (angle-sin(angle)), angle = 2*ACos((R-excess)/R)
 159
 160
 161   Float_t angle   = 2*TMath::ACos( (fConeSize-excess) / fConeSize );
 162
 163   Float_t coneA   = fConeSize*fConeSize*TMath::Pi(); // A = pi R^2, isolation cone area
 164
 165   Float_t excessA = fConeSize*fConeSize / 2 * (angle-TMath::Sin(angle));
 166
 167   if(coneA > excessA) return coneA / (coneA-excessA);
 168   else
 169   {
 170     printf("AliIsolationCut::CalculateExcessAreaFraction() - Please Check : Excess Track %2.3f, coneA %2.2f,  excessA %2.2f, angle %2.2f,factor %2.2f\n",
 171            excess,coneA, excessA, angle*TMath::RadToDeg(), coneA / (coneA-excessA));
 172     return  1;
 173   }
 174 }
 175
 176 //_______________________________________________________________________________________
 177 Float_t AliIsolationCut::GetCellDensity(AliAODPWG4ParticleCorrelation * pCandidate,
 178                                         AliCaloTrackReader * reader) const
 179 {
 180   // Get good cell density (number of active cells over all cells in cone)
 181
 182   Double_t coneCells    = 0.; //number of cells in cone with radius fConeSize
 183   Double_t coneCellsBad = 0.; //number of bad cells in cone with radius fConeSize
 184   Double_t cellDensity  = 1.;
 185
 186   Float_t phiC  = pCandidate->Phi() ;
 187   if(phiC<0) phiC+=TMath::TwoPi();
 188   Float_t etaC  = pCandidate->Eta() ;
 189
 190   if(pCandidate->GetDetector()=="EMCAL")
 191   {
 192     AliEMCALGeometry* eGeom = AliEMCALGeometry::GetInstance();
 193     AliCalorimeterUtils *cu = reader->GetCaloUtils();
 194
 195     Int_t absId = -999;
 196     if (eGeom->GetAbsCellIdFromEtaPhi(etaC,phiC,absId))
 197     {
 198       //Get absolute (col,row) of candidate
 199       Int_t iEta=-1, iPhi=-1, iRCU = -1;
 200       Int_t nSupMod = cu->GetModuleNumberCellIndexes(absId, pCandidate->GetDetector(), iEta, iPhi, iRCU);
 201
 202       Int_t colC = iEta;
 203       if (nSupMod % 2) colC =  AliEMCALGeoParams::fgkEMCALCols + iEta ;
 204       Int_t rowC = iPhi + AliEMCALGeoParams::fgkEMCALRows*int(nSupMod/2);
 205
 206       Int_t sqrSize = int(fConeSize/0.0143) ; // Size of cell in radians
 207       //loop on cells in a square of side fConeSize to check cells in cone
 208       for(Int_t icol = colC-sqrSize; icol < colC+sqrSize;icol++)
 209       {
 210         for(Int_t irow = rowC-sqrSize; irow < rowC+sqrSize; irow++)
 211         {
 212           if (Radius(colC, rowC, icol, irow) < sqrSize)
 213           {
 214             coneCells += 1.;
 215
 216             Int_t cellSM  = -999;
 217             Int_t cellEta = -999;
 218             Int_t cellPhi = -999;
 219             if(icol > AliEMCALGeoParams::fgkEMCALCols-1)
 220             {
 221               cellSM = 0+int(irow/AliEMCALGeoParams::fgkEMCALRows)*2;
 222               cellEta = icol-AliEMCALGeoParams::fgkEMCALCols;
 223               cellPhi = irow-AliEMCALGeoParams::fgkEMCALRows*int(cellSM/2);
 224             }
 225             if(icol < AliEMCALGeoParams::fgkEMCALCols)
 226             {
 227               cellSM = 1+int(irow/AliEMCALGeoParams::fgkEMCALRows)*2;
 228               cellEta = icol;
 229               cellPhi = irow-AliEMCALGeoParams::fgkEMCALRows*int(cellSM/2);
 230             }
 231
 232             //Count as bad "cells" out of EMCAL acceptance
 233             if(icol < 0 || icol > AliEMCALGeoParams::fgkEMCALCols*2 ||
 234                irow < 0 || irow > AliEMCALGeoParams::fgkEMCALRows*16./3) //5*nRows+1/3*nRows
 235             {
 236               coneCellsBad += 1.;
 237             }
 238             //Count as bad "cells" marked as bad in the DataBase
 239             else if (cu->GetEMCALChannelStatus(cellSM,cellEta,cellPhi)==1)
 240             {
 241               coneCellsBad += 1. ;
 242             }
 243           }
 244         }
 245       }//end of cells loop
 246     }
 247
 248     else if(fDebug>0) printf("cluster with bad (eta,phi) in EMCal for energy density calculation\n");
 249
 250     if (coneCells > 0.)
 251     {
 252       cellDensity = (coneCells-coneCellsBad)/coneCells;
 253       //printf("Energy density = %f\n", cellDensity);
 254     }
 255   }
 256
 257   return cellDensity;
 258
 259 }
 260
 261 //__________________________________________________________________________________
 262 void AliIsolationCut::GetCoeffNormBadCell(AliAODPWG4ParticleCorrelation * pCandidate,
 263                                           AliCaloTrackReader * reader,
 264                                           Float_t &  coneBadCellsCoeff,
 265                                           Float_t &  etaBandBadCellsCoeff,
 266                                           Float_t & phiBandBadCellsCoeff)
 267 {
 268   // Get good cell density (number of active cells over all cells in cone)
 269
 270   Double_t coneCells    = 0.; //number of cells in cone with radius fConeSize
 271   Double_t phiBandCells = 0.; //number of cells in band phi
 272   Double_t etaBandCells = 0.; //number of cells in band eta
 273
 274   Float_t phiC  = pCandidate->Phi() ;
 275   if(phiC<0) phiC+=TMath::TwoPi();
 276   Float_t etaC  = pCandidate->Eta() ;
 277
 278   if(pCandidate->GetDetector()=="EMCAL")
 279   {
 280     AliEMCALGeometry* eGeom = AliEMCALGeometry::GetInstance();
 281     AliCalorimeterUtils *cu = reader->GetCaloUtils();
 282
 283     Int_t absId = -999;
 284     if (eGeom->GetAbsCellIdFromEtaPhi(etaC,phiC,absId))
 285     {
 286       //Get absolute (col,row) of candidate
 287       Int_t iEta=-1, iPhi=-1, iRCU = -1;
 288       Int_t nSupMod = cu->GetModuleNumberCellIndexes(absId, pCandidate->GetDetector(),
 289                                                      iEta, iPhi, iRCU);
 290
 291       Int_t colC = iEta;
 292       if (nSupMod % 2) colC =  AliEMCALGeoParams::fgkEMCALCols + iEta ;
 293       Int_t rowC = iPhi + AliEMCALGeoParams::fgkEMCALRows*int(nSupMod/2);
 294
 295       Int_t sqrSize = int(fConeSize/0.0143) ; // Size of cell in radians
 296       for(Int_t icol = 0; icol < 2*AliEMCALGeoParams::fgkEMCALCols-1;icol++)
 297       {
 298         for(Int_t irow = 0; irow < 5*AliEMCALGeoParams::fgkEMCALRows -1; irow++)
 299         {
 300           //loop on cells in a square of side fConeSize to check cells in cone
 301           if     ( Radius(colC, rowC, icol, irow) < sqrSize ) { coneCells    += 1.; }
 302           else if( icol>colC-sqrSize  &&  icol<colC+sqrSize ) { phiBandCells += 1 ; }
 303           else if( irow>rowC-sqrSize  &&  irow<rowC+sqrSize ) { etaBandCells += 1 ; }
 304
 305           Int_t cellSM  = -999;
 306           Int_t cellEta = -999;
 307           Int_t cellPhi = -999;
 308           if(icol > AliEMCALGeoParams::fgkEMCALCols-1)
 309           {
 310             cellSM = 0+int(irow/AliEMCALGeoParams::fgkEMCALRows)*2;
 311             cellEta = icol-AliEMCALGeoParams::fgkEMCALCols;
 312             cellPhi = irow-AliEMCALGeoParams::fgkEMCALRows*int(cellSM/2);
 313           }
 314           if(icol < AliEMCALGeoParams::fgkEMCALCols)
 315           {
 316             cellSM = 1+int(irow/AliEMCALGeoParams::fgkEMCALRows)*2;
 317             cellEta = icol;
 318             cellPhi = irow-AliEMCALGeoParams::fgkEMCALRows*int(cellSM/2);
 319           }
 320
 321           if( (icol < 0 || icol > AliEMCALGeoParams::fgkEMCALCols*2-1 ||
 322                irow < 0 || irow > AliEMCALGeoParams::fgkEMCALRows*5 - 1) //5*nRows+1/3*nRows //Count as bad "cells" out of EMCAL acceptance
 323              || (cu->GetEMCALChannelStatus(cellSM,cellEta,cellPhi)==1))  //Count as bad "cells" marked as bad in the DataBase
 324           {
 325             if     ( Radius(colC, rowC, icol, irow) < sqrSize ) coneBadCellsCoeff    += 1.;
 326             else if( icol>colC-sqrSize  &&  icol<colC+sqrSize ) phiBandBadCellsCoeff += 1 ;
 327                   else if( irow>rowC-sqrSize  &&  irow<rowC+sqrSize ) etaBandBadCellsCoeff += 1 ;
 328           }
 329         }
 330       }//end of cells loop
 331     }
 332
 333     else if(fDebug > 0) printf("cluster with bad (eta,phi) in EMCal for energy density coeff calculation\n");
 334
 335     if (coneCells > 0.)
 336     {
 337       //   printf("Energy density coneBadCellsCoeff= %.2f coneCells%.2f\n", coneBadCellsCoeff,coneCells);
 338       coneBadCellsCoeff = (coneCells-coneBadCellsCoeff)/coneCells;
 339       //  printf("coneBadCellsCoeff= %.2f\n", coneBadCellsCoeff);
 340     }
 341     if (phiBandCells > 0.)
 342     {
 343       // printf("Energy density phiBandBadCellsCoeff = %.2f phiBandCells%.2f\n", phiBandBadCellsCoeff,phiBandCells);
 344       phiBandBadCellsCoeff = (phiBandCells-phiBandBadCellsCoeff)/phiBandCells;
 345       // printf("phiBandBadCellsCoeff = %.2f\n", phiBandBadCellsCoeff);
 346     }
 347     if (etaBandCells > 0.)
 348     {
 349       //printf("Energy density etaBandBadCellsCoeff = %.2f etaBandCells%.2f\n", etaBandBadCellsCoeff,etaBandCells);
 350       etaBandBadCellsCoeff = (etaBandCells-etaBandBadCellsCoeff)/etaBandCells;
 351       // printf("etaBandBadCellsCoeff = %.2f\n",etaBandBadCellsCoeff);
 352     }
 353
 354   }
 355
 356 }
 357
 358 //____________________________________________
 359 TString AliIsolationCut::GetICParametersList()
 360 {
 361   //Put data member values in string to keep in output container
 362
 363   TString parList ; //this will be list of parameters used for this analysis.
 364   const Int_t buffersize = 255;
 365   char onePar[buffersize] ;
 366
 367   snprintf(onePar,buffersize,"--- AliIsolationCut ---\n") ;
 368   parList+=onePar ;
 369   snprintf(onePar,buffersize,"fConeSize: (isolation cone size) %1.2f\n",fConeSize) ;
 370   parList+=onePar ;
 371   snprintf(onePar,buffersize,"fPtThreshold >%2.2f;<%2.2f (isolation pt threshold) \n",fPtThreshold,fPtThresholdMax) ;
 372   parList+=onePar ;
 373   snprintf(onePar,buffersize,"fSumPtThreshold >%2.2f;<%2.2f (isolation sum pt threshold) \n",fSumPtThreshold,fSumPtThresholdMax) ;
 374   parList+=onePar ;
 375   snprintf(onePar,buffersize,"fPtFraction=%2.2f (isolation pt threshold fraction) \n",fPtFraction) ;
 376   parList+=onePar ;
 377   snprintf(onePar,buffersize,"fICMethod=%d (isolation cut case) \n",fICMethod) ;
 378   parList+=onePar ;
 379   snprintf(onePar,buffersize,"fPartInCone=%d \n",fPartInCone) ;
 380   parList+=onePar ;
 381   snprintf(onePar,buffersize,"fFracIsThresh=%i \n",fFracIsThresh) ;
 382   parList+=onePar ;
 383
 384   return parList;
 385 }
 386
 387 //____________________________________
 388 void AliIsolationCut::InitParameters()
 389 {
 390   //Initialize the parameters of the analysis.
 391
 392   fConeSize       = 0.4 ;
 393   fPtThreshold    = 0.5  ;
 394   fPtThresholdMax = 10000.  ;
 395   fSumPtThreshold    = 1.0 ;
 396   fSumPtThresholdMax = 10000. ;
 397   fPtFraction     = 0.1 ;
 398   fPartInCone     = kNeutralAndCharged;
 399   fICMethod       = kSumPtIC; // 0 pt threshol method, 1 cone pt sum method
 400   fFracIsThresh   = 1;
 401 }
 402
 403 //________________________________________________________________________________
 404 void  AliIsolationCut::MakeIsolationCut(TObjArray * plCTS,
 405                                         TObjArray * plNe,
 406                                         AliCaloTrackReader * reader,
 407                                         AliCaloPID * pid,
 408                                         Bool_t bFillAOD,
 409                                         AliAODPWG4ParticleCorrelation  *pCandidate,
 410                                         TString aodArrayRefName,
 411                                         Int_t   & n,
 412                                         Int_t   & nfrac,
 413                                         Float_t & coneptsum,
 414                                         Bool_t  & isolated) const
 415 {
 416   //Search in cone around a candidate particle if it is isolated
 417   Float_t ptC   = pCandidate->Pt() ;
 418   Float_t phiC  = pCandidate->Phi() ;
 419   if(phiC<0) phiC+=TMath::TwoPi();
 420   Float_t etaC  = pCandidate->Eta() ;
 421
 422   Float_t ptLead = -100. ;
 423   Float_t pt     = -100. ;
 424   Float_t eta    = -100. ;
 425   Float_t phi    = -100. ;
 426   Float_t rad    = -100. ;
 427
 428   Float_t coneptsumCluster = 0;
 429   Float_t coneptsumTrack   = 0;
 430
 431   Float_t  etaBandPtSumTrack   = 0;
 432   Float_t  phiBandPtSumTrack   = 0;
 433   Float_t  etaBandPtSumCluster = 0;
 434   Float_t  phiBandPtSumCluster = 0;
 435
 436   n         = 0 ;
 437   nfrac     = 0 ;
 438   isolated  = kFALSE;
 439
 440   if(fDebug>0)
 441   {
 442     printf("AliIsolationCut::MakeIsolationCut() - Cadidate pT %2.2f, eta %2.2f, phi %2.2f, cone %1.2f, thres %2.2f, Fill AOD? %d",
 443            pCandidate->Pt(), pCandidate->Eta(), pCandidate->Phi()*TMath::RadToDeg(), fConeSize,fPtThreshold,bFillAOD);
 444     if(plCTS) printf(", nTracks %d"  ,plCTS->GetEntriesFast());
 445     if(plNe)  printf(", nClusters %d",plNe ->GetEntriesFast());
 446
 447     printf("\n");
 448   }
 449
 450   //Initialize the array with refrences
 451   TObjArray * refclusters  = 0x0;
 452   TObjArray * reftracks    = 0x0;
 453   Int_t       ntrackrefs   = 0;
 454   Int_t       nclusterrefs = 0;
 455
 456   //Check charged particles in cone.
 457   if(plCTS &&
 458      (fPartInCone==kOnlyCharged || fPartInCone==kNeutralAndCharged))
 459   {
 460     TVector3 p3;
 461     for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ )
 462     {
 463       AliVTrack* track = dynamic_cast<AliVTrack*>(plCTS->At(ipr)) ;
 464
 465       if(track)
 466       {
 467         //Do not count the candidate (pion, conversion photon) or the daughters of the candidate
 468         if(track->GetID() == pCandidate->GetTrackLabel(0) || track->GetID() == pCandidate->GetTrackLabel(1) ||
 469            track->GetID() == pCandidate->GetTrackLabel(2) || track->GetID() == pCandidate->GetTrackLabel(3)   ) continue ;
 470
 471         p3.SetXYZ(track->Px(),track->Py(),track->Pz());
 472         pt  = p3.Pt();
 473         eta = p3.Eta();
 474         phi = p3.Phi() ;
 475       }
 476       else
 477       {// Mixed event stored in AliAODPWG4Particles
 478         AliAODPWG4Particle * trackmix = dynamic_cast<AliAODPWG4Particle*>(plCTS->At(ipr)) ;
 479         if(!trackmix)
 480         {
 481           printf("AliIsolationCut::MakeIsolationCut() - Wrong track data type, continue\n");
 482           continue;
 483         }
 484
 485         pt  = trackmix->Pt();
 486         eta = trackmix->Eta();
 487         phi = trackmix->Phi() ;
 488       }
 489
 490       if( phi < 0 ) phi+=TMath::TwoPi();
 491
 492       rad = Radius(etaC, phiC, eta, phi);
 493
 494       // ** For the background out of cone **
 495
 496       if(rad > fConeSize)
 497       {
 498         if(eta > (etaC-fConeSize) && eta < (etaC+fConeSize)) phiBandPtSumTrack += pt;
 499         if(phi > (phiC-fConeSize) && phi < (phiC+fConeSize)) etaBandPtSumTrack += pt;
 500       }
 501
 502       // ** For the isolated particle **
 503
 504       // Only loop the particle at the same side of candidate
 505       if(TMath::Abs(phi-phiC) > TMath::PiOver2()) continue ;
 506
 507       // If at the same side has particle larger than candidate,
 508       // then candidate can not be the leading, skip such events
 509       if(pt > ptC)
 510       {
 511         n         = -1;
 512         nfrac     = -1;
 513         coneptsumTrack = -1;
 514         isolated  = kFALSE;
 515
 516         pCandidate->SetLeadingParticle(kFALSE);
 517
 518         if(bFillAOD && reftracks)
 519         {
 520           reftracks->Clear();
 521           delete reftracks;
 522         }
 523
 524         return ;
 525       }
 526
 527       // // Check if there is any particle inside cone with pt larger than  fPtThreshold
 528       // Check if the leading particule inside the cone has a ptLead larger than fPtThreshold
 529
 530       if( fDebug > 0 )
 531         printf("\t track %d, pT %2.2f, eta %1.2f, phi %2.2f, R candidate %2.2f", ipr,pt,eta,phi,rad);
 532
 533       if(rad < fConeSize)
 534       {
 535         if(fDebug > 0)  printf(" -  inside candidate cone");
 536
 537         if(bFillAOD)
 538         {
 539           ntrackrefs++;
 540           if(ntrackrefs == 1)
 541           {
 542             reftracks = new TObjArray(0);
 543             //reftracks->SetName(Form("Tracks%s",aodArrayRefName.Data()));
 544             TString tempo(aodArrayRefName)  ;
 545             tempo += "Tracks" ;
 546             reftracks->SetName(tempo);
 547             reftracks->SetOwner(kFALSE);
 548           }
 549           reftracks->Add(track);
 550         }
 551
 552         coneptsumTrack+=pt;
 553
 554         if( ptLead < pt ) ptLead = pt;
 555
 556 //        // *Before*, count particles in cone
 557 //        if(pt > fPtThreshold && pt < fPtThresholdMax)  n++;
 558 //
 559 //        //if fPtFraction*ptC<fPtThreshold then consider the fPtThreshold directly
 560 //        if(fFracIsThresh)
 561 //        {
 562 //          if( fPtFraction*ptC < fPtThreshold )
 563 //          {
 564 //            if( pt > fPtThreshold )    nfrac++ ;
 565 //          }
 566 //          else
 567 //          {
 568 //            if( pt > fPtFraction*ptC ) nfrac++;
 569 //          }
 570 //        }
 571 //        else
 572 //        {
 573 //          if( pt > fPtFraction*ptC ) nfrac++;
 574 //        }
 575
 576       } // Inside cone
 577
 578       if( fDebug > 0 )  printf("\n");
 579
 580     }// charged particle loop
 581
 582   }//Tracks
 583
 584
 585   //Check neutral particles in cone.
 586   if(plNe &&
 587      (fPartInCone==kOnlyNeutral || fPartInCone==kNeutralAndCharged))
 588   {
 589     TLorentzVector mom ;
 590
 591     for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ )
 592     {
 593       AliVCluster * calo = dynamic_cast<AliVCluster *>(plNe->At(ipr)) ;
 594
 595       if(calo)
 596       {
 597         //Get the index where the cluster comes, to retrieve the corresponding vertex
 598         Int_t evtIndex = 0 ;
 599         if (reader->GetMixedEvent())
 600           evtIndex=reader->GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
 601
 602
 603         //Do not count the candidate (photon or pi0) or the daughters of the candidate
 604         if(calo->GetID() == pCandidate->GetCaloLabel(0) ||
 605            calo->GetID() == pCandidate->GetCaloLabel(1)   ) continue ;
 606
 607         //Skip matched clusters with tracks in case of neutral+charged analysis
 608         if( fPartInCone == kNeutralAndCharged &&
 609            pid->IsTrackMatched(calo,reader->GetCaloUtils(),reader->GetInputEvent()) ) continue ;
 610
 611         //Assume that come from vertex in straight line
 612         calo->GetMomentum(mom,reader->GetVertex(evtIndex)) ;
 613
 614         pt  = mom.Pt()  ;
 615         eta = mom.Eta() ;
 616         phi = mom.Phi() ;
 617       }
 618       else
 619       {// Mixed event stored in AliAODPWG4Particles
 620         AliAODPWG4Particle * calomix = dynamic_cast<AliAODPWG4Particle*>(plNe->At(ipr)) ;
 621         if(!calomix)
 622         {
 623           printf("AliIsolationCut::MakeIsolationCut() - Wrong calo data type, continue\n");
 624           continue;
 625         }
 626
 627         pt  = calomix->Pt();
 628         eta = calomix->Eta();
 629         phi = calomix->Phi() ;
 630       }
 631
 632       if( phi < 0 ) phi+=TMath::TwoPi();
 633
 634       rad = Radius(etaC, phiC, eta, phi);
 635
 636       // ** For the background out of cone **
 637
 638       if(rad > fConeSize)
 639       {
 640         if(eta > (etaC-fConeSize) && eta < (etaC+fConeSize)) phiBandPtSumCluster += pt;
 641         if(phi > (phiC-fConeSize) && phi < (phiC+fConeSize)) etaBandPtSumCluster += pt;
 642       }
 643
 644       // ** For the isolated particle **
 645
 646       // Only loop the particle at the same side of candidate
 647       if(TMath::Abs(phi-phiC)>TMath::PiOver2()) continue ;
 648
 649       // If at the same side has particle larger than candidate,
 650       // then candidate can not be the leading, skip such events
 651       if(pt > ptC)
 652       {
 653         n         = -1;
 654         nfrac     = -1;
 655         coneptsumCluster = -1;
 656         isolated  = kFALSE;
 657
 658         pCandidate->SetLeadingParticle(kFALSE);
 659
 660         if(bFillAOD)
 661         {
 662           if(reftracks)
 663           {
 664             reftracks  ->Clear();
 665             delete reftracks;
 666           }
 667
 668           if(refclusters)
 669           {
 670             refclusters->Clear();
 671             delete refclusters;
 672           }
 673         }
 674         return ;
 675       }
 676
 677       //Check if there is any particle inside cone with pt larger than  fPtThreshold
 678
 679       if(fDebug > 0 )
 680         printf("\t cluster %d, pT %2.2f, eta %1.2f, phi %2.2f, R candidate %2.2f", ipr,pt,eta,phi,rad);
 681
 682       if(rad < fConeSize)
 683       {
 684         if(fDebug > 0 )  printf(" - inside candidate cone");
 685
 686         if(bFillAOD)
 687         {
 688           nclusterrefs++;
 689           if(nclusterrefs==1)
 690           {
 691             refclusters = new TObjArray(0);
 692             //refclusters->SetName(Form("Clusters%s",aodArrayRefName.Data()));
 693             TString tempo(aodArrayRefName)  ;
 694             tempo += "Clusters" ;
 695             refclusters->SetName(tempo);
 696             refclusters->SetOwner(kFALSE);
 697           }
 698           refclusters->Add(calo);
 699         }
 700
 701         coneptsumCluster+=pt;
 702
 703         if( ptLead < pt ) ptLead = pt;
 704
 705 //        // *Before*, count particles in cone
 706 //        if(pt > fPtThreshold && pt < fPtThresholdMax)  n++;
 707 //
 708 //        //if fPtFraction*ptC<fPtThreshold then consider the fPtThreshold directly
 709 //        if(fFracIsThresh)
 710 //        {
 711 //          if( fPtFraction*ptC < fPtThreshold )
 712 //          {
 713 //            if( pt > fPtThreshold )    nfrac++ ;
 714 //          }
 715 //          else
 716 //          {
 717 //            if( pt > fPtFraction*ptC ) nfrac++;
 718 //          }
 719 //        }
 720 //        else
 721 //        {
 722 //          if( pt > fPtFraction*ptC ) nfrac++;
 723 //        }
 724
 725       }//in cone
 726
 727       if(fDebug > 0 )  printf("\n");
 728
 729     }// neutral particle loop
 730
 731   }//neutrals
 732
 733   //Add reference arrays to AOD when filling AODs only
 734   if(bFillAOD)
 735   {
 736     if(refclusters)     pCandidate->AddObjArray(refclusters);
 737     if(reftracks)         pCandidate->AddObjArray(reftracks);
 738   }
 739
 740   coneptsum = coneptsumCluster + coneptsumTrack;
 741
 742   // *Now*, just check the leading particle in the cone if the threshold is passed
 743   if(ptLead > fPtThreshold && ptLead < fPtThresholdMax)  n = 1;
 744
 745   //if fPtFraction*ptC<fPtThreshold then consider the fPtThreshold directly
 746   if(fFracIsThresh)
 747   {
 748     if( fPtFraction*ptC < fPtThreshold )
 749     {
 750       if( ptLead > fPtThreshold )    nfrac = 1 ;
 751     }
 752     else
 753     {
 754       if( ptLead > fPtFraction*ptC ) nfrac = 1;
 755     }
 756   }
 757   else
 758   {
 759     if( ptLead > fPtFraction*ptC ) nfrac = 1;
 760   }
 761
 762   //-------------------------------------------------------------------
 763   //Check isolation, depending on selected isolation criteria requested
 764
 765   if( fICMethod == kPtThresIC)
 766   {
 767     if( n == 0 ) isolated = kTRUE ;
 768   }
 769   else if( fICMethod == kSumPtIC )
 770   {
 771     if( coneptsum > fSumPtThreshold &&
 772         coneptsum < fSumPtThresholdMax )
 773       isolated  =  kFALSE ;
 774     else
 775       isolated  =  kTRUE  ;
 776   }
 777   else if( fICMethod == kPtFracIC )
 778   {
 779     if(nfrac == 0 ) isolated = kTRUE ;
 780   }
 781   else if( fICMethod == kSumPtFracIC )
 782   {
 783     //when the fPtFraction*ptC < fSumPtThreshold then consider the later case
 784     // printf("photon analysis IsDataMC() ?%i\n",IsDataMC());
 785     if( fFracIsThresh )
 786     {
 787       if( fPtFraction*ptC < fSumPtThreshold  && coneptsum < fSumPtThreshold ) isolated  =  kTRUE ;
 788       if( fPtFraction*ptC > fSumPtThreshold  && coneptsum < fPtFraction*ptC ) isolated  =  kTRUE ;
 789     }
 790     else
 791     {
 792       if( coneptsum < fPtFraction*ptC ) isolated  =  kTRUE ;
 793     }
 794   }
 795   else if( fICMethod == kSumDensityIC )
 796   {
 797     // Get good cell density (number of active cells over all cells in cone)
 798     // and correct energy in cone
 799
 800     Float_t cellDensity = GetCellDensity(pCandidate,reader);
 801
 802     if( coneptsum < fSumPtThreshold*cellDensity )
 803       isolated = kTRUE;
 804   }
 805   else if( fICMethod == kSumBkgSubIC )
 806   {
 807     Double_t coneptsumBkg = 0.;
 808     Float_t  etaBandPtSumTrackNorm   = 0;
 809     Float_t  phiBandPtSumTrackNorm   = 0;
 810     Float_t  etaBandPtSumClusterNorm = 0;
 811     Float_t  phiBandPtSumClusterNorm = 0;
 812
 813     Float_t  excessFracEtaTrack   = 1;
 814     Float_t  excessFracPhiTrack   = 1;
 815     Float_t  excessFracEtaCluster = 1;
 816     Float_t  excessFracPhiCluster = 1;
 817
 818     // Normalize background to cone area
 819     if     (fPartInCone != kOnlyCharged       )
 820       CalculateUEBandClusterNormalization(reader, etaC, phiC,
 821                                           phiBandPtSumCluster    , etaBandPtSumCluster,
 822                                           phiBandPtSumClusterNorm, etaBandPtSumClusterNorm,
 823                                           excessFracEtaCluster   , excessFracPhiCluster    );
 824     if     (fPartInCone != kOnlyNeutral       )
 825       CalculateUEBandTrackNormalization(reader, etaC, phiC,
 826                                           phiBandPtSumTrack    , etaBandPtSumTrack  ,
 827                                           phiBandPtSumTrackNorm, etaBandPtSumTrackNorm,
 828                                           excessFracEtaTrack   , excessFracPhiTrack    );
 829
 830     if     (fPartInCone == kOnlyCharged       ) coneptsumBkg = etaBandPtSumTrackNorm;
 831     else if(fPartInCone == kOnlyNeutral       ) coneptsumBkg = etaBandPtSumClusterNorm;
 832     else if(fPartInCone == kNeutralAndCharged ) coneptsumBkg = etaBandPtSumClusterNorm + etaBandPtSumTrackNorm;
 833
 834     //coneptsumCluster*=(coneBadCellsCoeff*excessFracEtaCluster*excessFracPhiCluster) ; // apply this correction earlier???
 835     // line commented out in last modif!!!
 836
 837     coneptsum = coneptsumCluster+coneptsumTrack;
 838
 839     coneptsum -= coneptsumBkg;
 840
 841     if( coneptsum > fSumPtThreshold && coneptsum < fSumPtThresholdMax )
 842       isolated  =  kFALSE ;
 843     else
 844       isolated  =  kTRUE  ;
 845
 846   }
 847
 848 }
 849
 850 //_____________________________________________________
 851 void AliIsolationCut::Print(const Option_t * opt) const
 852 {
 853
 854   //Print some relevant parameters set for the analysis
 855   if(! opt)
 856     return;
 857
 858   printf("**** Print %s %s **** \n", GetName(), GetTitle() ) ;
 859
 860   printf("IC method          =     %d\n",    fICMethod   ) ;
 861   printf("Cone Size          =     %1.2f\n", fConeSize   ) ;
 862   printf("pT threshold       =     >%2.1f;<%2.1f\n", fPtThreshold   ,   fPtThresholdMax) ;
 863   printf("Sum pT threshold   =     >%2.1f;<%2.1f\n", fSumPtThreshold,fSumPtThresholdMax) ;
 864   printf("pT fraction        =     %3.1f\n", fPtFraction ) ;
 865   printf("particle type in cone =  %d\n",    fPartInCone ) ;
 866   printf("using fraction for high pt leading instead of frac ? %i\n",fFracIsThresh);
 867   printf("    \n") ;
 868
 869 }
 870
 871 //___________________________________________________________________________
 872 Float_t AliIsolationCut::Radius(Float_t etaC, Float_t phiC,
 873                                 Float_t eta , Float_t phi) const
 874 {
 875   // Calculate the distance to trigger from any particle
 876
 877   Float_t dEta = etaC-eta;
 878   Float_t dPhi = phiC-phi;
 879
 880   if(TMath::Abs(dPhi) >= TMath::Pi())
 881     dPhi = TMath::TwoPi()-TMath::Abs(dPhi);
 882
 883   return TMath::Sqrt( dEta*dEta + dPhi*dPhi );
 884
 885 }
 886
 887
 888