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