403e64d81fdcc05d9414ab8cd80842c5bdf92860
[u/mrichter/AliRoot.git] / 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     if(fDebug > 0 )
770       printf("pT Cand %2.2f, pT Lead %2.2f, %2.2f<pT Lead< %2.2f, isolated %d\n",
771              ptC,ptLead,fPtThreshold,fPtThresholdMax,isolated);
772   }
773   else if( fICMethod == kSumPtIC )
774   {
775     if( coneptsum > fSumPtThreshold &&
776         coneptsum < fSumPtThresholdMax )
777       isolated  =  kFALSE ;
778     else
779       isolated  =  kTRUE  ;
780     
781     if(fDebug > 0 )
782       printf("pT Cand %2.2f, SumPt %2.2f, %2.2f<Sum pT< %2.2f, isolated %d\n",
783              ptC,ptLead,fSumPtThreshold,fSumPtThresholdMax,isolated);
784   }
785   else if( fICMethod == kPtFracIC )
786   {
787     if(nfrac == 0 ) isolated = kTRUE ;
788   }
789   else if( fICMethod == kSumPtFracIC )
790   {
791     //when the fPtFraction*ptC < fSumPtThreshold then consider the later case
792     // printf("photon analysis IsDataMC() ?%i\n",IsDataMC());
793     if( fFracIsThresh )
794     {
795       if( fPtFraction*ptC < fSumPtThreshold  && coneptsum < fSumPtThreshold ) isolated  =  kTRUE ;
796       if( fPtFraction*ptC > fSumPtThreshold  && coneptsum < fPtFraction*ptC ) isolated  =  kTRUE ;
797     }
798     else 
799     {
800       if( coneptsum < fPtFraction*ptC ) isolated  =  kTRUE ;
801     }
802   }
803   else if( fICMethod == kSumDensityIC )
804   {    
805     // Get good cell density (number of active cells over all cells in cone)
806     // and correct energy in cone
807     
808     Float_t cellDensity = GetCellDensity(pCandidate,reader);
809     
810     if( coneptsum < fSumPtThreshold*cellDensity )
811       isolated = kTRUE;
812   }
813   else if( fICMethod == kSumBkgSubIC )
814   {
815     Double_t coneptsumBkg = 0.;
816     Float_t  etaBandPtSumTrackNorm   = 0;
817     Float_t  phiBandPtSumTrackNorm   = 0;
818     Float_t  etaBandPtSumClusterNorm = 0;
819     Float_t  phiBandPtSumClusterNorm = 0;
820     
821     Float_t  excessFracEtaTrack   = 1;
822     Float_t  excessFracPhiTrack   = 1;
823     Float_t  excessFracEtaCluster = 1;
824     Float_t  excessFracPhiCluster = 1;
825     
826     // Normalize background to cone area
827     if     (fPartInCone != kOnlyCharged       )
828       CalculateUEBandClusterNormalization(reader, etaC, phiC,
829                                           phiBandPtSumCluster    , etaBandPtSumCluster,
830                                           phiBandPtSumClusterNorm, etaBandPtSumClusterNorm,
831                                           excessFracEtaCluster   , excessFracPhiCluster    );
832     if     (fPartInCone != kOnlyNeutral       )
833       CalculateUEBandTrackNormalization(reader, etaC, phiC,
834                                           phiBandPtSumTrack    , etaBandPtSumTrack  ,
835                                           phiBandPtSumTrackNorm, etaBandPtSumTrackNorm,
836                                           excessFracEtaTrack   , excessFracPhiTrack    );
837     
838     if     (fPartInCone == kOnlyCharged       ) coneptsumBkg = etaBandPtSumTrackNorm;
839     else if(fPartInCone == kOnlyNeutral       ) coneptsumBkg = etaBandPtSumClusterNorm;
840     else if(fPartInCone == kNeutralAndCharged ) coneptsumBkg = etaBandPtSumClusterNorm + etaBandPtSumTrackNorm;
841     
842     //coneptsumCluster*=(coneBadCellsCoeff*excessFracEtaCluster*excessFracPhiCluster) ; // apply this correction earlier???
843     // line commented out in last modif!!!
844     
845     coneptsum = coneptsumCluster+coneptsumTrack;
846     
847     coneptsum -= coneptsumBkg;
848     
849     if( coneptsum > fSumPtThreshold && coneptsum < fSumPtThresholdMax )
850       isolated  =  kFALSE ;
851     else
852       isolated  =  kTRUE  ;
853
854   }
855   
856 }
857
858 //_____________________________________________________
859 void AliIsolationCut::Print(const Option_t * opt) const
860 {
861   
862   //Print some relevant parameters set for the analysis
863   if(! opt)
864     return;
865   
866   printf("**** Print %s %s **** \n", GetName(), GetTitle() ) ;
867   
868   printf("IC method          =     %d\n",    fICMethod   ) ; 
869   printf("Cone Size          =     %1.2f\n", fConeSize   ) ; 
870   printf("pT threshold       =     >%2.1f;<%2.1f\n", fPtThreshold   ,   fPtThresholdMax) ;
871   printf("Sum pT threshold   =     >%2.1f;<%2.1f\n", fSumPtThreshold,fSumPtThresholdMax) ;
872   printf("pT fraction        =     %3.1f\n", fPtFraction ) ;
873   printf("particle type in cone =  %d\n",    fPartInCone ) ;
874   printf("using fraction for high pt leading instead of frac ? %i\n",fFracIsThresh);
875   printf("    \n") ;
876   
877
878
879 //___________________________________________________________________________
880 Float_t AliIsolationCut::Radius(Float_t etaC, Float_t phiC,
881                                 Float_t eta , Float_t phi) const
882 {
883   // Calculate the distance to trigger from any particle
884
885   Float_t dEta = etaC-eta;
886   Float_t dPhi = phiC-phi;
887   
888   if(TMath::Abs(dPhi) >= TMath::Pi()) 
889     dPhi = TMath::TwoPi()-TMath::Abs(dPhi);
890   
891   return TMath::Sqrt( dEta*dEta + dPhi*dPhi );
892   
893 }
894
895
896