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