Adding a reminder for coders
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALRecoUtils.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 /* $Id: AliEMCALRecoUtils.cxx 33808 2009-07-15 09:48:08Z gconesab $ */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 //
20 // Class AliEMCALRecoUtils
21 // Some utilities to recalculate the cluster position or energy linearity
22 //
23 //
24 // Author:  Gustavo Conesa (LPSC- Grenoble) 
25 //          Track matching part: Rongrong Ma (Yale)
26
27 ///////////////////////////////////////////////////////////////////////////////
28 // --- standard c ---
29
30 // standard C++ includes
31 //#include <Riostream.h>
32
33 // ROOT includes
34 #include <TGeoManager.h>
35 #include <TGeoMatrix.h>
36 #include <TGeoBBox.h>
37 #include <TH2F.h>
38 #include <TArrayI.h>
39 #include <TArrayF.h>
40 #include "TObjArray.h"
41
42 // STEER includes
43 #include "AliVCluster.h"
44 #include "AliVCaloCells.h"
45 #include "AliVEvent.h"
46 #include "AliLog.h"
47 #include "AliPID.h"
48 #include "AliESDEvent.h"
49 #include "AliAODEvent.h"
50 #include "AliESDtrack.h"
51 #include "AliAODTrack.h"
52 #include "AliExternalTrackParam.h"
53 #include "AliESDfriendTrack.h"
54 #include "AliTrackerBase.h"
55
56 // EMCAL includes
57 #include "AliEMCALRecoUtils.h"
58 #include "AliEMCALGeometry.h"
59 #include "AliEMCALTrack.h"
60 #include "AliEMCALCalibTimeDepCorrection.h"
61 #include "AliEMCALPIDUtils.h"
62
63 ClassImp(AliEMCALRecoUtils)
64   
65 //______________________________________________
66 AliEMCALRecoUtils::AliEMCALRecoUtils():
67   fNonLinearityFunction (kNoCorrection), fParticleType(kPhoton),
68   fPosAlgo(kUnchanged),fW0(4.), fNonLinearThreshold(30),
69   fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
70   fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(),
71   fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
72   fAODFilterMask(32),
73   fMatchedTrackIndex(0x0), fMatchedClusterIndex(0x0), 
74   fResidualEta(0x0), fResidualPhi(0x0), fCutEtaPhiSum(kTRUE), fCutEtaPhiSeparate(kFALSE), fCutR(0.1), fCutEta(0.02), fCutPhi(0.04), fMass(0.139), fStep(50),
75   fRejectExoticCluster(kFALSE),
76   fTrackCutsType(kTPCOnlyCut), fCutMinTrackPt(0), fCutMinNClusterTPC(-1), fCutMinNClusterITS(-1), fCutMaxChi2PerClusterTPC(1e10), fCutMaxChi2PerClusterITS(1e10),
77   fCutRequireTPCRefit(kFALSE), fCutRequireITSRefit(kFALSE), fCutAcceptKinkDaughters(kFALSE),
78   fCutMaxDCAToVertexXY(1e10), fCutMaxDCAToVertexZ(1e10),fCutDCAToVertex2D(kFALSE),fPIDUtils(),
79   fUseTimeCorrectionFactors(kFALSE),  fTimeCorrectionFactorsSet(kFALSE)
80 {
81 //
82   // Constructor.
83   // Initialize all constant values which have to be used
84   // during Reco algorithm execution
85   //
86   
87   //Misalignment matrices
88   for(Int_t i = 0; i < 15 ; i++) {
89       fMisalTransShift[i] = 0.; 
90       fMisalRotShift[i]   = 0.; 
91   }
92   
93   //Non linearity
94   for(Int_t i = 0; i < 7  ; i++) fNonLinearityParams[i] = 0.; 
95
96   //For kBeamTestCorrected case, but default is no correction
97   fNonLinearityParams[0] =  0.99078;
98   fNonLinearityParams[1] =  0.161499;
99   fNonLinearityParams[2] =  0.655166; 
100   fNonLinearityParams[3] =  0.134101;
101   fNonLinearityParams[4] =  163.282;
102   fNonLinearityParams[5] =  23.6904;
103   fNonLinearityParams[6] =  0.978;
104   
105   //For kPi0GammaGamma case
106   //fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
107   //fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
108   //fNonLinearityParams[2] = 1.046;
109
110   //Track matching
111   fMatchedTrackIndex     = new TArrayI();
112   fMatchedClusterIndex   = new TArrayI();
113   fResidualPhi           = new TArrayF();
114   fResidualEta           = new TArrayF();
115   
116   InitTrackCuts();
117   fPIDUtils              = new AliEMCALPIDUtils();
118 }
119
120 //______________________________________________________________________
121 AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco) 
122 : TNamed(reco), fNonLinearityFunction(reco.fNonLinearityFunction), 
123   fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo), fW0(reco.fW0), fNonLinearThreshold(reco.fNonLinearThreshold),
124   fRecalibration(reco.fRecalibration),fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
125   fRemoveBadChannels(reco.fRemoveBadChannels),fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
126   fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
127   fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
128   fAODFilterMask(reco.fAODFilterMask),
129   fMatchedTrackIndex(reco.fMatchedTrackIndex?new TArrayI(*reco.fMatchedTrackIndex):0x0),
130   fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
131   fResidualEta(reco.fResidualEta?new TArrayF(*reco.fResidualEta):0x0),
132   fResidualPhi(reco.fResidualPhi?new TArrayF(*reco.fResidualPhi):0x0),
133   fCutEtaPhiSum(reco.fCutEtaPhiSum), fCutEtaPhiSeparate(reco.fCutEtaPhiSeparate), fCutR(reco.fCutR), fCutEta(reco.fCutEta), fCutPhi(reco.fCutPhi),
134   fMass(reco.fMass), fStep(reco.fStep),
135   fRejectExoticCluster(reco.fRejectExoticCluster),
136   fTrackCutsType(reco.fTrackCutsType), fCutMinTrackPt(reco.fCutMinTrackPt), fCutMinNClusterTPC(reco.fCutMinNClusterTPC), fCutMinNClusterITS(reco.fCutMinNClusterITS), 
137   fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC), fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS),
138   fCutRequireTPCRefit(reco.fCutRequireTPCRefit), fCutRequireITSRefit(reco.fCutRequireITSRefit),
139   fCutAcceptKinkDaughters(reco.fCutAcceptKinkDaughters),
140   fCutMaxDCAToVertexXY(reco.fCutMaxDCAToVertexXY), fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ),fCutDCAToVertex2D(reco.fCutDCAToVertex2D),
141   fPIDUtils(reco.fPIDUtils), 
142   fUseTimeCorrectionFactors(reco.fUseTimeCorrectionFactors),  fTimeCorrectionFactorsSet(reco.fTimeCorrectionFactorsSet)
143 {
144   //Copy ctor
145   
146   for(Int_t i = 0; i < 15 ; i++) {
147       fMisalRotShift[i] = reco.fMisalRotShift[i]; 
148       fMisalTransShift[i] = reco.fMisalTransShift[i]; 
149   } 
150   for(Int_t i = 0; i < 7  ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i]; 
151
152 }
153
154
155 //______________________________________________________________________
156 AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco) 
157 {
158   //Assignment operator
159   
160   if(this == &reco)return *this;
161   ((TNamed *)this)->operator=(reco);
162
163   fNonLinearityFunction      = reco.fNonLinearityFunction;
164   fParticleType              = reco.fParticleType;
165   fPosAlgo                   = reco.fPosAlgo; 
166   fW0                        = reco.fW0;
167   fNonLinearThreshold        = reco.fNonLinearThreshold;
168   fRecalibration             = reco.fRecalibration;
169   fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
170   fRemoveBadChannels         = reco.fRemoveBadChannels;
171   fRecalDistToBadChannels    = reco.fRecalDistToBadChannels;
172   fEMCALBadChannelMap        = reco.fEMCALBadChannelMap;
173   fNCellsFromEMCALBorder     = reco.fNCellsFromEMCALBorder;
174   fNoEMCALBorderAtEta0       = reco.fNoEMCALBorderAtEta0;
175
176
177   for(Int_t i = 0; i < 15 ; i++) {fMisalTransShift[i] = reco.fMisalTransShift[i]; fMisalRotShift[i] = reco.fMisalRotShift[i];}
178   for(Int_t i = 0; i < 7  ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i]; 
179
180   fAODFilterMask              = reco.fAODFilterMask;
181   
182   fCutEtaPhiSum              = reco.fCutEtaPhiSum;
183   fCutEtaPhiSeparate         = reco.fCutEtaPhiSeparate;
184   fCutR                      = reco.fCutR;
185   fCutEta                    = reco.fCutEta;
186   fCutPhi                    = reco.fCutPhi;
187   fMass                      = reco.fMass;
188   fStep                      = reco.fStep;
189   fRejectExoticCluster       = reco.fRejectExoticCluster;
190
191   fTrackCutsType             = reco.fTrackCutsType;
192   fCutMinTrackPt             = reco.fCutMinTrackPt;
193   fCutMinNClusterTPC         = reco.fCutMinNClusterTPC;
194   fCutMinNClusterITS         = reco.fCutMinNClusterITS; 
195   fCutMaxChi2PerClusterTPC   = reco.fCutMaxChi2PerClusterTPC;
196   fCutMaxChi2PerClusterITS   = reco.fCutMaxChi2PerClusterITS;
197   fCutRequireTPCRefit        = reco.fCutRequireTPCRefit;
198   fCutRequireITSRefit        = reco.fCutRequireITSRefit;
199   fCutAcceptKinkDaughters    = reco.fCutAcceptKinkDaughters;
200   fCutMaxDCAToVertexXY       = reco.fCutMaxDCAToVertexXY;
201   fCutMaxDCAToVertexZ        = reco.fCutMaxDCAToVertexZ;
202   fCutDCAToVertex2D          = reco.fCutDCAToVertex2D;
203
204   fPIDUtils                  = reco.fPIDUtils;
205   
206   fUseTimeCorrectionFactors  = reco.fUseTimeCorrectionFactors;
207   fTimeCorrectionFactorsSet  = reco.fTimeCorrectionFactorsSet;
208
209   
210   if(reco.fResidualEta){
211     // assign or copy construct
212     if(fResidualEta){ 
213       *fResidualEta = *reco.fResidualEta;
214     }
215     else fResidualEta = new TArrayF(*reco.fResidualEta);
216   }
217   else{
218     if(fResidualEta)delete fResidualEta;
219     fResidualEta = 0;
220   }
221   
222   if(reco.fResidualPhi){
223     // assign or copy construct
224     if(fResidualPhi){ 
225       *fResidualPhi = *reco.fResidualPhi;
226     }
227     else fResidualPhi = new TArrayF(*reco.fResidualPhi);
228   }
229   else{
230     if(fResidualPhi)delete fResidualPhi;
231     fResidualPhi = 0;
232   }
233   
234   if(reco.fMatchedTrackIndex){
235     // assign or copy construct
236     if(fMatchedTrackIndex){ 
237       *fMatchedTrackIndex = *reco.fMatchedTrackIndex;
238     }
239     else fMatchedTrackIndex = new TArrayI(*reco.fMatchedTrackIndex);
240   }
241   else{
242     if(fMatchedTrackIndex)delete fMatchedTrackIndex;
243     fMatchedTrackIndex = 0;
244   }  
245   
246   if(reco.fMatchedClusterIndex){
247     // assign or copy construct
248     if(fMatchedClusterIndex){ 
249       *fMatchedClusterIndex = *reco.fMatchedClusterIndex;
250     }
251     else fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex);
252   }
253   else{
254     if(fMatchedClusterIndex)delete fMatchedClusterIndex;
255     fMatchedClusterIndex = 0;
256   }
257    
258   return *this;
259 }
260
261
262 //__________________________________________________
263 AliEMCALRecoUtils::~AliEMCALRecoUtils()
264 {
265   //Destructor.
266         
267         if(fEMCALRecalibrationFactors) { 
268                 fEMCALRecalibrationFactors->Clear();
269                 delete  fEMCALRecalibrationFactors;
270         }       
271   
272   if(fEMCALBadChannelMap) { 
273                 fEMCALBadChannelMap->Clear();
274                 delete  fEMCALBadChannelMap;
275         }
276  
277   delete fMatchedTrackIndex   ; 
278   delete fMatchedClusterIndex ; 
279   delete fResidualEta         ; 
280   delete fResidualPhi         ; 
281
282 }
283
284 //_______________________________________________________________
285 Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells) 
286 {
287         // Given the list of AbsId of the cluster, get the maximum cell and 
288         // check if there are fNCellsFromBorder from the calorimeter border
289         
290   if(!cluster){
291     AliInfo("Cluster pointer null!");
292     return kFALSE;
293   }
294   
295   //If the distance to the border is 0 or negative just exit accept all clusters
296         if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
297   
298   Int_t absIdMax        = -1, iSM =-1, ieta = -1, iphi = -1;
299   Bool_t shared = kFALSE;
300   GetMaxEnergyCell(geom, cells, cluster, absIdMax,  iSM, ieta, iphi, shared);
301
302   AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n", 
303            absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
304         
305         if(absIdMax==-1) return kFALSE;
306         
307         //Check if the cell is close to the borders:
308         Bool_t okrow = kFALSE;
309         Bool_t okcol = kFALSE;
310   
311   if(iSM < 0 || iphi < 0 || ieta < 0 ) {
312     AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
313                   iSM,ieta,iphi));
314   }
315   
316   //Check rows/phi
317   if(iSM < 10){
318     if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE; 
319   }
320   else{
321     if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE; 
322   }
323   
324   //Check columns/eta
325   if(!fNoEMCALBorderAtEta0){
326     if(ieta  > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE; 
327   }
328   else{
329     if(iSM%2==0){
330       if(ieta >= fNCellsFromEMCALBorder)     okcol = kTRUE;     
331     }
332     else {
333       if(ieta <  48-fNCellsFromEMCALBorder)  okcol = kTRUE;     
334     }
335   }//eta 0 not checked
336     
337   AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d:  column? %d, row? %d\nq",
338            fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
339         
340         if (okcol && okrow) {
341     //printf("Accept\n");
342     return kTRUE;
343   }
344         else  {
345     //printf("Reject\n");
346     AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
347     return kFALSE;
348   }
349         
350 }       
351
352
353 //_________________________________________________________________________________________________________
354 Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, const Int_t nCells){
355         // Check that in the cluster cells, there is no bad channel of those stored 
356         // in fEMCALBadChannelMap or fPHOSBadChannelMap
357         
358         if(!fRemoveBadChannels)  return kFALSE;
359         if(!fEMCALBadChannelMap) return kFALSE;
360         
361         Int_t icol = -1;
362         Int_t irow = -1;
363         Int_t imod = -1;
364         for(Int_t iCell = 0; iCell<nCells; iCell++){
365                 
366                 //Get the column and row
367     Int_t iTower = -1, iIphi = -1, iIeta = -1; 
368     geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta); 
369     if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
370     geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);                      
371     if(GetEMCALChannelStatus(imod, icol, irow)){
372       AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
373       return kTRUE;
374     }
375                 
376         }// cell cluster loop
377         
378         return kFALSE;
379         
380 }
381
382 //_________________________________________________
383 Bool_t AliEMCALRecoUtils::IsExoticCluster(AliVCluster *cluster) const {
384   // Check if the cluster has high energy  but small number of cells
385   // The criteria comes from Gustavo's study
386   //
387   
388   if(!cluster){
389     AliInfo("Cluster pointer null!");
390     return kFALSE;
391   }
392   
393   Int_t nc = cluster->GetNCells() ;
394   
395   if      ( nc > 8 )                   return kFALSE ; // Good cluster, needed for 3x3 clusterizer  
396   else if ( nc < 1 + cluster->E()/3. ) return kTRUE  ; // Bad cluster
397   else                                 return kFALSE ; // Good cluster
398   
399 }
400
401 //__________________________________________________
402 Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
403 // Correct cluster energy from non linearity functions
404   
405   if(!cluster){
406     AliInfo("Cluster pointer null!");
407     return 0;
408   }
409   
410   
411   Float_t energy = cluster->E();
412   
413   switch (fNonLinearityFunction) {
414       
415     case kPi0MC:
416     {
417       //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
418       //Double_t fNonLinearityParams[0] = 1.014;
419       //Double_t fNonLinearityParams[1] = -0.03329;
420       //Double_t fNonLinearityParams[2] = -0.3853;
421       //Double_t fNonLinearityParams[3] = 0.5423;
422       //Double_t fNonLinearityParams[4] = -0.4335;
423        energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
424                   ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
425                     exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
426       break;
427     }
428      
429     case kPi0GammaGamma:
430     {
431       //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
432       //Double_t fNonLinearityParams[0] = 1.04;
433       //Double_t fNonLinearityParams[1] = -0.1445;
434       //Double_t fNonLinearityParams[2] = 1.046;
435       energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
436       break;
437     }
438       
439     case kPi0GammaConversion:
440     {
441       //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
442       //fNonLinearityParams[0] = 0.139393/0.1349766;
443       //fNonLinearityParams[1] = 0.0566186;
444       //fNonLinearityParams[2] = 0.982133;
445       energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
446       
447       break;
448     }
449       
450     case kBeamTest:
451     {
452       //From beam test, Alexei's results, for different ZS thresholds
453       //                        th=30 MeV; th = 45 MeV; th = 75 MeV
454       //fNonLinearityParams[0] = 1.007;      1.003;      1.002 
455       //fNonLinearityParams[1] = 0.894;      0.719;      0.797 
456       //fNonLinearityParams[2] = 0.246;      0.334;      0.358 
457       //Rescale the param[0] with 1.03
458       energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
459       
460       break;
461     }
462       
463     case kBeamTestCorrected:
464     {
465       //From beam test, corrected for material between beam and EMCAL
466       //fNonLinearityParams[0] =  0.99078
467       //fNonLinearityParams[1] =  0.161499;
468       //fNonLinearityParams[2] =  0.655166; 
469       //fNonLinearityParams[3] =  0.134101;
470       //fNonLinearityParams[4] =  163.282;
471       //fNonLinearityParams[5] =  23.6904;
472       //fNonLinearityParams[6] =  0.978;
473         energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
474
475       break;
476     }
477       
478     case kNoCorrection:
479       AliDebug(2,"No correction on the energy\n");
480       break;
481       
482   }
483   
484   return energy;
485
486 }
487 //__________________________________________________
488 void AliEMCALRecoUtils::InitNonLinearityParam()
489 {
490         //Initialising Non Linearity Parameters
491         
492         if(fNonLinearityFunction == kPi0MC)
493                 {
494                         fNonLinearityParams[0] = 1.014;
495       fNonLinearityParams[1] = -0.03329;
496       fNonLinearityParams[2] = -0.3853;
497       fNonLinearityParams[3] = 0.5423;
498       fNonLinearityParams[4] = -0.4335;
499                 }
500
501         if(fNonLinearityFunction == kPi0GammaGamma)
502                 {
503                         fNonLinearityParams[0] = 1.04;
504                         fNonLinearityParams[1] = -0.1445;
505                         fNonLinearityParams[2] = 1.046;
506                 }       
507
508         if(fNonLinearityFunction == kPi0GammaConversion)
509                 {
510       fNonLinearityParams[0] = 0.139393;
511       fNonLinearityParams[1] = 0.0566186;
512       fNonLinearityParams[2] = 0.982133;
513                 }       
514
515         if(fNonLinearityFunction == kBeamTest)
516                 {
517                         if(fNonLinearThreshold == 30)
518                         {
519                                 fNonLinearityParams[0] = 1.007; 
520                                 fNonLinearityParams[1] = 0.894; 
521                                 fNonLinearityParams[2] = 0.246; 
522                         }
523                         if(fNonLinearThreshold == 45)
524                         {
525                                 fNonLinearityParams[0] = 1.003; 
526                                 fNonLinearityParams[1] = 0.719; 
527                                 fNonLinearityParams[2] = 0.334; 
528                         }
529                         if(fNonLinearThreshold == 75)
530                         {
531                                 fNonLinearityParams[0] = 1.002; 
532                                 fNonLinearityParams[1] = 0.797; 
533                                 fNonLinearityParams[2] = 0.358; 
534                         }
535                 }
536
537         if(fNonLinearityFunction == kBeamTestCorrected)
538                 {
539                         fNonLinearityParams[0] =  0.99078;
540                         fNonLinearityParams[1] =  0.161499;
541                         fNonLinearityParams[2] =  0.655166; 
542                         fNonLinearityParams[3] =  0.134101;
543                         fNonLinearityParams[4] =  163.282;
544                         fNonLinearityParams[5] =  23.6904;
545                         fNonLinearityParams[6] =  0.978;
546                 }
547 }
548
549 //__________________________________________________
550 Float_t  AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const 
551 {
552   //Calculate shower depth for a given cluster energy and particle type
553
554   // parameters 
555   Float_t x0    = 1.31;
556   Float_t ecr   = 8;
557   Float_t depth = 0;
558   
559   switch ( iParticle )
560   {
561     case kPhoton:
562       depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
563       break;
564       
565     case kElectron:
566       depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
567       break;
568       
569     case kHadron:
570       // hadron 
571       // boxes anc. here
572       if(gGeoManager){
573         gGeoManager->cd("ALIC_1/XEN1_1");
574         TGeoNode        *geoXEn1    = gGeoManager->GetCurrentNode();
575         TGeoNodeMatrix  *geoSM      = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
576         if(geoSM){
577           TGeoVolume      *geoSMVol   = geoSM->GetVolume(); 
578           TGeoShape       *geoSMShape = geoSMVol->GetShape();
579           TGeoBBox        *geoBox     = dynamic_cast<TGeoBBox *>(geoSMShape);
580           if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
581           else AliFatal("Null GEANT box");
582         }else AliFatal("NULL  GEANT node matrix");
583       }
584       else{//electron
585         depth = x0 * (TMath::Log(energy*1000 / ecr)  - 0.5); //Multiply energy by 1000 to transform to MeV
586       }
587         
588       break;
589       
590     default://photon
591       depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
592   }  
593   
594   return depth;
595   
596 }
597
598 //__________________________________________________
599 void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu, 
600                                          Int_t & absId,  Int_t& iSupMod, Int_t& ieta, Int_t& iphi, Bool_t &shared)
601 {
602   //For a given CaloCluster gets the absId of the cell 
603   //with maximum energy deposit.
604   
605   Double_t eMax        = -1.;
606   Double_t eCell       = -1.;
607   Float_t  fraction    = 1.;
608   Float_t  recalFactor = 1.;
609   Int_t    cellAbsId   = -1 ;
610
611   Int_t iTower  = -1;
612   Int_t iIphi   = -1;
613   Int_t iIeta   = -1;
614   Int_t iSupMod0= -1;
615
616   if(!clu){
617     AliInfo("Cluster pointer null!");
618     absId=-1; iSupMod0=-1, ieta = -1; iphi = -1; shared = -1;
619     return;
620   }
621   
622   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
623     cellAbsId = clu->GetCellAbsId(iDig);
624     fraction  = clu->GetCellAmplitudeFraction(iDig);
625     //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
626     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
627     geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta); 
628     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
629     if(iDig==0) iSupMod0=iSupMod;
630     else if(iSupMod0!=iSupMod) {
631       shared = kTRUE;
632       //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
633     }
634     if(IsRecalibrationOn()) {
635       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
636     }
637     eCell  = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
638     //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
639     if(eCell > eMax)  { 
640       eMax  = eCell; 
641       absId = cellAbsId;
642       //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
643     }
644   }// cell loop
645   
646   //Get from the absid the supermodule, tower and eta/phi numbers
647   geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta); 
648   //Gives SuperModule and Tower numbers
649   geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
650                                          iIphi, iIeta,iphi,ieta); 
651   //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
652   //printf("Max end---\n");
653   
654 }
655
656 //________________________________________________________________
657 void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
658         //Init EMCAL recalibration factors
659         AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
660         //In order to avoid rewriting the same histograms
661         Bool_t oldStatus = TH1::AddDirectoryStatus();
662         TH1::AddDirectory(kFALSE);
663   
664         fEMCALRecalibrationFactors = new TObjArray(10);
665         for (int i = 0; i < 10; i++) fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i),  48, 0, 48, 24, 0, 24));
666         //Init the histograms with 1
667         for (Int_t sm = 0; sm < 10; sm++) {
668                 for (Int_t i = 0; i < 48; i++) {
669                         for (Int_t j = 0; j < 24; j++) {
670                                 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
671                         }
672                 }
673         }
674         fEMCALRecalibrationFactors->SetOwner(kTRUE);
675         fEMCALRecalibrationFactors->Compress();
676         
677         //In order to avoid rewriting the same histograms
678         TH1::AddDirectory(oldStatus);           
679 }
680
681
682 //________________________________________________________________
683 void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
684         //Init EMCAL bad channels map
685         AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
686         //In order to avoid rewriting the same histograms
687         Bool_t oldStatus = TH1::AddDirectoryStatus();
688         TH1::AddDirectory(kFALSE);
689         
690         fEMCALBadChannelMap = new TObjArray(10);
691         //TH2F * hTemp = new  TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
692         for (int i = 0; i < 10; i++) {
693                 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
694         }
695         
696         //delete hTemp;
697         
698         fEMCALBadChannelMap->SetOwner(kTRUE);
699         fEMCALBadChannelMap->Compress();
700         
701         //In order to avoid rewriting the same histograms
702         TH1::AddDirectory(oldStatus);           
703 }
704
705 //________________________________________________________________
706 void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells){
707         // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
708         
709   if(!cluster){
710     AliInfo("Cluster pointer null!");
711     return;
712   }  
713   
714         //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
715         UShort_t * index    = cluster->GetCellsAbsId() ;
716         Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
717         Int_t ncells = cluster->GetNCells();
718         
719         //Initialize some used variables
720         Float_t energy = 0;
721         Int_t absId    = -1;
722   Int_t icol = -1, irow = -1, imod=1;
723         Float_t factor = 1, frac = 0;
724         
725         //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
726         for(Int_t icell = 0; icell < ncells; icell++){
727                 absId = index[icell];
728                 frac =  fraction[icell];
729                 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
730                 Int_t iTower = -1, iIphi = -1, iIeta = -1; 
731                 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta); 
732                 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
733                 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);                  
734                 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
735     AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
736              imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
737                 
738                 energy += cells->GetCellAmplitude(absId)*factor*frac;
739         }
740         
741         
742                 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
743         
744         cluster->SetE(energy);
745         
746 }
747
748
749 //__________________________________________________
750 void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
751 {
752   //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
753   
754   if(!clu){
755     AliInfo("Cluster pointer null!");
756     return;
757   }
758     
759   if     (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
760   else if(fPosAlgo==kPosTowerIndex)  RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
761   else   AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
762   
763 }  
764
765 //__________________________________________________
766 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
767 {
768   // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
769   // The algorithm is a copy of what is done in AliEMCALRecPoint
770   
771   Double_t eCell       = 0.;
772   Float_t  fraction    = 1.;
773   Float_t  recalFactor = 1.;
774   
775   Int_t    absId   = -1;
776   Int_t    iTower  = -1, iIphi  = -1, iIeta  = -1;
777   Int_t    iSupModMax = -1, iSM=-1, iphi   = -1, ieta   = -1;
778   Float_t  weight = 0.,  totalWeight=0.;
779   Float_t  newPos[3] = {0,0,0};
780   Double_t pLocal[3], pGlobal[3];
781   Bool_t shared = kFALSE;
782
783   Float_t  clEnergy = clu->E(); //Energy already recalibrated previously
784   GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi,shared);
785   Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
786   
787   //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
788   
789   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
790     absId = clu->GetCellAbsId(iDig);
791     fraction  = clu->GetCellAmplitudeFraction(iDig);
792     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
793     geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta); 
794     geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);                       
795     
796     if(IsRecalibrationOn()) {
797       recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
798     }
799     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
800     
801     weight = GetCellWeight(eCell,clEnergy);
802     //printf("cell energy %f, weight %f\n",eCell,weight);
803     totalWeight += weight;
804     geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
805     //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
806     geom->GetGlobal(pLocal,pGlobal,iSupModMax);
807     //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
808
809     for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
810     
811   }// cell loop
812   
813   if(totalWeight>0){
814     for(int i=0; i<3; i++ )    newPos[i] /= totalWeight;
815   }
816     
817   //Float_t pos[]={0,0,0};
818   //clu->GetPosition(pos);
819   //printf("OldPos  : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
820   //printf("NewPos  : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
821   
822         if(iSupModMax > 1) {//sector 1
823           newPos[0] +=fMisalTransShift[3];//-=3.093; 
824           newPos[1] +=fMisalTransShift[4];//+=6.82;
825           newPos[2] +=fMisalTransShift[5];//+=1.635;
826     //printf("   +    : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
827
828         }
829         else {//sector 0
830           newPos[0] +=fMisalTransShift[0];//+=1.134;
831           newPos[1] +=fMisalTransShift[1];//+=8.2;
832           newPos[2] +=fMisalTransShift[2];//+=1.197;
833     //printf("   +    : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
834
835         }
836   //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
837
838   clu->SetPosition(newPos);
839   
840 }  
841
842 //__________________________________________________
843 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
844 {
845   // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
846   // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
847   
848   Double_t eCell       = 1.;
849   Float_t  fraction    = 1.;
850   Float_t  recalFactor = 1.;
851   
852   Int_t absId   = -1;
853   Int_t iTower  = -1;
854   Int_t iIphi   = -1, iIeta   = -1;
855         Int_t iSupMod = -1, iSupModMax = -1;
856   Int_t iphi = -1, ieta =-1;
857   Bool_t shared = kFALSE;
858
859   Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
860   GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi,shared);
861   Float_t  depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
862
863   Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
864   Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
865   Int_t startingSM = -1;
866   
867   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
868     absId = clu->GetCellAbsId(iDig);
869     fraction  = clu->GetCellAmplitudeFraction(iDig);
870     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
871     geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta); 
872     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);                   
873     
874     if     (iDig==0)  startingSM = iSupMod;
875     else if(iSupMod != startingSM) areInSameSM = kFALSE;
876
877     eCell  = cells->GetCellAmplitude(absId);
878     
879     if(IsRecalibrationOn()) {
880       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
881     }
882     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
883     
884     weight = GetCellWeight(eCell,clEnergy);
885     if(weight < 0) weight = 0;
886     totalWeight += weight;
887     weightedCol += ieta*weight;
888     weightedRow += iphi*weight;
889     
890     //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
891     
892     }// cell loop
893     
894   Float_t xyzNew[]={0.,0.,0.};
895   if(areInSameSM == kTRUE) {
896     //printf("In Same SM\n");
897     weightedCol = weightedCol/totalWeight;
898     weightedRow = weightedRow/totalWeight;
899     geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew); 
900   }
901   else {
902     //printf("In Different SM\n");
903     geom->RecalculateTowerPosition(iphi,        ieta,        iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew); 
904   }
905   
906   clu->SetPosition(xyzNew);
907   
908 }
909
910 //____________________________________________________________________________
911 void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster){           
912         
913   //re-evaluate distance to bad channel with updated bad map
914   
915   if(!fRecalDistToBadChannels) return;
916   
917   if(!cluster){
918     AliInfo("Cluster pointer null!");
919     return;
920   }  
921   
922         //Get channels map of the supermodule where the cluster is.
923   Int_t absIdMax        = -1, iSupMod =-1, icolM = -1, irowM = -1;
924   Bool_t shared = kFALSE;
925   GetMaxEnergyCell(geom, cells, cluster, absIdMax,  iSupMod, icolM, irowM, shared);
926   TH2D* hMap  = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
927
928   Int_t dRrow, dRcol;   
929         Float_t  minDist = 10000.;
930         Float_t  dist    = 0.;
931   
932   //Loop on tower status map 
933         for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
934                 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
935                         //Check if tower is bad.
936                         if(hMap->GetBinContent(icol,irow)==0) continue;
937       //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
938       //       iSupMod,icol, irow, icolM,irowM);
939       
940       dRrow=TMath::Abs(irowM-irow);
941       dRcol=TMath::Abs(icolM-icol);
942       dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
943                         if(dist < minDist){
944         //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
945         minDist = dist;
946       }
947       
948                 }
949         }
950   
951         //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
952         if (shared) {
953                 TH2D* hMap2 = 0;
954                 Int_t iSupMod2 = -1;
955     
956                 //The only possible combinations are (0,1), (2,3) ... (8,9)
957                 if(iSupMod%2) iSupMod2 = iSupMod-1;
958                 else          iSupMod2 = iSupMod+1;
959                 hMap2  = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
960     
961                 //Loop on tower status map of second super module
962                 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
963                         for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
964                                 //Check if tower is bad.
965                                 if(hMap2->GetBinContent(icol,irow)==0) continue;
966                                 //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels(shared) - \n \t Bad channel in SM %d, col %d, row %d \n \t Cluster max in SM %d, col %d, row %d\n",
967           //     iSupMod2,icol, irow,iSupMod,icolM,irowM);
968         
969         dRrow=TMath::Abs(irow-irowM);
970         
971         if(iSupMod%2) {
972                                   dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
973                                 }
974         else {
975           dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
976                                 }                    
977         
978                                 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
979         if(dist < minDist) minDist = dist;        
980         
981                         }
982                 }
983     
984         }// shared cluster in 2 SuperModules
985   
986   AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
987   cluster->SetDistanceToBadChannel(minDist);
988   
989 }
990
991 //____________________________________________________________________________
992 void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){           
993         
994   //re-evaluate identification parameters with bayesian
995   
996   if(!cluster){
997     AliInfo("Cluster pointer null!");
998     return;
999   }
1000   
1001         if ( cluster->GetM02() != 0)
1002     fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
1003   
1004   Float_t pidlist[AliPID::kSPECIESN+1];
1005         for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
1006   
1007   cluster->SetPID(pidlist);
1008         
1009 }
1010
1011 //____________________________________________________________________________
1012 void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster)
1013 {
1014   // Calculates new center of gravity in the local EMCAL-module coordinates 
1015   // and tranfers into global ALICE coordinates
1016   // Calculates Dispersion and main axis
1017   
1018   if(!cluster){
1019     AliInfo("Cluster pointer null!");
1020     return;
1021   }
1022     
1023   Int_t nstat  = 0;
1024   Float_t wtot = 0. ;
1025   Double_t eCell       = 0.;
1026   Float_t  fraction    = 1.;
1027   Float_t  recalFactor = 1.;
1028
1029   Int_t iSupMod = -1;
1030   Int_t iTower  = -1;
1031   Int_t iIphi   = -1;
1032   Int_t iIeta   = -1;
1033   Int_t iphi    = -1;
1034   Int_t ieta    = -1;
1035   Double_t etai = -1.;
1036   Double_t phii = -1.;
1037   
1038   Double_t w     = 0.;
1039   Double_t d     = 0.;
1040   Double_t dxx   = 0.;
1041   Double_t dzz   = 0.;
1042   Double_t dxz   = 0.;  
1043   Double_t xmean = 0.;
1044   Double_t zmean = 0.;
1045     
1046   //Loop on cells
1047   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
1048     
1049     //Get from the absid the supermodule, tower and eta/phi numbers
1050     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta); 
1051     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);        
1052     
1053     //Get the cell energy, if recalibration is on, apply factors
1054     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
1055     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1056     if(IsRecalibrationOn()) {
1057       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1058     }
1059     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1060     
1061     if(cluster->E() > 0 && eCell > 0){
1062       
1063       w  = GetCellWeight(eCell,cluster->E());
1064       
1065       etai=(Double_t)ieta;
1066       phii=(Double_t)iphi;              
1067       if(w > 0.0) {
1068         wtot += w ;
1069         nstat++;                        
1070         //Shower shape
1071         dxx  += w * etai * etai ;
1072         xmean+= w * etai ;
1073         dzz  += w * phii * phii ;
1074         zmean+= w * phii ; 
1075         dxz  += w * etai * phii ; 
1076       }
1077     }
1078     else
1079       AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1080   }//cell loop
1081   
1082   //Normalize to the weight     
1083   if (wtot > 0) {
1084     xmean /= wtot ;
1085     zmean /= wtot ;
1086   }
1087   else
1088     AliError(Form("Wrong weight %f\n", wtot));
1089   
1090   //Calculate dispersion        
1091   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
1092     
1093     //Get from the absid the supermodule, tower and eta/phi numbers
1094     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta); 
1095     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1096     
1097     //Get the cell energy, if recalibration is on, apply factors
1098     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
1099     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1100     if(IsRecalibrationOn()) {
1101       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1102     }
1103     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1104     
1105     if(cluster->E() > 0 && eCell > 0){
1106       
1107       w  = GetCellWeight(eCell,cluster->E());
1108       
1109       etai=(Double_t)ieta;
1110       phii=(Double_t)iphi;              
1111       if(w > 0.0)  d +=  w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean)); 
1112     }
1113     else
1114       AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1115   }// cell loop
1116   
1117   //Normalize to the weigth and set shower shape parameters
1118   if (wtot > 0 && nstat > 1) {
1119     d /= wtot ;
1120     dxx /= wtot ;
1121     dzz /= wtot ;
1122     dxz /= wtot ;
1123     dxx -= xmean * xmean ;
1124     dzz -= zmean * zmean ;
1125     dxz -= xmean * zmean ;
1126     cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
1127     cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
1128   }
1129   else{
1130     d=0. ;
1131     cluster->SetM20(0.) ;
1132     cluster->SetM02(0.) ;
1133   }     
1134   
1135   if (d>=0)
1136     cluster->SetDispersion(TMath::Sqrt(d)) ;
1137   else    
1138     cluster->SetDispersion(0) ;
1139 }
1140
1141 //____________________________________________________________________________
1142 void AliEMCALRecoUtils::FindMatches(AliVEvent *event,TObjArray * clusterArr,  AliEMCALGeometry *geom)
1143 {
1144   //This function should be called before the cluster loop
1145   //Before call this function, please recalculate the cluster positions
1146   //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
1147   //Store matched cluster indexes and residuals
1148   
1149   fMatchedTrackIndex->Reset();
1150   fMatchedClusterIndex->Reset();
1151   fResidualPhi->Reset();
1152   fResidualEta->Reset();
1153   
1154   fMatchedTrackIndex->Set(500);
1155   fMatchedClusterIndex->Set(500);
1156   fResidualPhi->Set(500);
1157   fResidualEta->Set(500);
1158   
1159   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1160   AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
1161   
1162   Int_t    matched=0;
1163   Double_t cv[21];
1164   for (Int_t i=0; i<21;i++) cv[i]=0;
1165   for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
1166   {
1167     AliExternalTrackParam *trackParam = 0;
1168     
1169     //If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner 
1170     if(esdevent)
1171     {
1172       AliESDtrack *esdTrack = esdevent->GetTrack(itr);
1173       if(!esdTrack || !IsAccepted(esdTrack)) continue;
1174       if(esdTrack->Pt()<fCutMinTrackPt) continue;
1175       const AliESDfriendTrack*  friendTrack = esdTrack->GetFriendTrack();
1176       if(friendTrack && friendTrack->GetTPCOut())
1177       {
1178         //Use TPC Out as starting point if it is available
1179         trackParam=  const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
1180       }
1181       else
1182       {
1183         //Otherwise use TPC inner
1184         trackParam =  const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam());
1185       }
1186     }
1187     
1188     //If the input event is AOD, the starting point for extrapolation is at vertex
1189     //AOD tracks are selected according to its bit.
1190     else if(aodevent)
1191     {
1192       AliAODTrack *aodTrack = aodevent->GetTrack(itr);
1193       if(!aodTrack) continue;
1194       if(!aodTrack->TestFilterMask(fAODFilterMask)) continue; //Select AOD tracks that fulfill GetStandardITSTPCTrackCuts2010()
1195       if(aodTrack->Pt()<fCutMinTrackPt) continue;
1196       Double_t pos[3],mom[3];
1197       aodTrack->GetXYZ(pos);
1198       aodTrack->GetPxPyPz(mom);
1199       AliDebug(5,Form("aod track: i=%d | pos=(%5.4f,%5.4f,%5.4f) | mom=(%5.4f,%5.4f,%5.4f) | charge=%d\n",itr,pos[0],pos[1],pos[2],mom[0],mom[1],mom[2],aodTrack->Charge()));
1200       trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
1201     }
1202     
1203     //Return if the input data is not "AOD" or "ESD"
1204     else
1205     {
1206       printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n");
1207       return;
1208     }
1209     
1210     if(!trackParam) continue;
1211     
1212     Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
1213     Int_t index = -1;
1214     if(!clusterArr){// get clusters from event
1215       for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1216       {
1217         AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1218         if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1219         AliExternalTrackParam trkPamTmp(*trackParam);//Retrieve the starting point every time before the extrapolation  
1220         Float_t tmpEta=-999, tmpPhi=-999;
1221         if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
1222         if(fCutEtaPhiSum)
1223         {
1224           Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1225           if(tmpR<dRMax)
1226           {
1227             dRMax=tmpR;
1228             dEtaMax=tmpEta;
1229             dPhiMax=tmpPhi;
1230             index=icl;
1231           }
1232         }
1233         else if(fCutEtaPhiSeparate)
1234         {
1235           if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
1236           {
1237             dEtaMax = tmpEta;
1238             dPhiMax = tmpPhi;
1239             index=icl;
1240           }
1241         }
1242         else
1243         {
1244           printf("Error: please specify your cut criteria\n");
1245           printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1246           printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1247           if(aodevent && trackParam) delete trackParam;
1248           return;
1249         }
1250       }//cluster loop
1251     } 
1252     else { // external cluster array, not from ESD event
1253       for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
1254       {
1255         AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
1256         if(!cluster){ 
1257           AliInfo("Cluster not found!!!");
1258           continue;
1259         }
1260         if(!cluster->IsEMCAL()) continue;
1261         AliExternalTrackParam trkPamTmp (*trackParam);//Retrieve the starting point every time before the extrapolation
1262         Float_t tmpEta=-999, tmpPhi=-999;
1263         if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
1264         if(fCutEtaPhiSum)
1265         {
1266           Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1267           if(tmpR<dRMax)
1268           {
1269             dRMax=tmpR;
1270             dEtaMax=tmpEta;
1271             dPhiMax=tmpPhi;
1272             index=icl;
1273           }
1274         }
1275         else if(fCutEtaPhiSeparate)
1276         {
1277           if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
1278           {
1279             dEtaMax = tmpEta;
1280             dPhiMax = tmpPhi;
1281             index=icl;
1282           }
1283         }
1284         else
1285         {
1286           printf("Error: please specify your cut criteria\n");
1287           printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1288           printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1289           if(aodevent && trackParam) delete trackParam;
1290           return;
1291         }
1292       }//cluster loop
1293     }// external list of clusters
1294     
1295     if(index>-1)
1296     {
1297       fMatchedTrackIndex  ->AddAt(itr,matched);
1298       fMatchedClusterIndex->AddAt(index,matched);
1299       fResidualEta          ->AddAt(dEtaMax,matched);
1300       fResidualPhi          ->AddAt(dPhiMax,matched);
1301       matched++;
1302     }
1303     if(aodevent && trackParam) delete trackParam;
1304   }//track loop
1305   
1306   AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1307   
1308   fMatchedTrackIndex  ->Set(matched);
1309   fMatchedClusterIndex->Set(matched);
1310   fResidualPhi          ->Set(matched);
1311   fResidualEta          ->Set(matched);
1312 }
1313
1314 //________________________________________________________________________________
1315 Int_t AliEMCALRecoUtils::FindMatchedCluster(AliESDtrack *track, AliVEvent *event, AliEMCALGeometry *geom)
1316 {
1317   //
1318   // This function returns the index of matched cluster to input track
1319   // Returns -1 if no match is found
1320   
1321   
1322   Float_t dRMax = fCutR, dEtaMax = fCutEta, dPhiMax = fCutPhi;
1323   Int_t index = -1;
1324   
1325   AliExternalTrackParam *trackParam=0;
1326   const AliESDfriendTrack*  friendTrack = track->GetFriendTrack();
1327   if(friendTrack && friendTrack->GetTPCOut())
1328     trackParam= const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
1329   else
1330     trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
1331   
1332   if(!trackParam) return index;   
1333   for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1334   {
1335     AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1336     if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;   
1337     AliExternalTrackParam trkPamTmp (*trackParam);//Retrieve the starting point every time before the extrapolation
1338     Float_t tmpEta=-999, tmpPhi=-999;
1339     if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
1340     if(fCutEtaPhiSum)
1341     {
1342       Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1343       if(tmpR<dRMax)
1344             {
1345               dRMax=tmpR;
1346               dEtaMax=tmpEta;
1347               dPhiMax=tmpPhi;
1348               index=icl;
1349             }
1350     }
1351     else if(fCutEtaPhiSeparate)
1352     {
1353       if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
1354             {
1355               dEtaMax = tmpEta;
1356               dPhiMax = tmpPhi;
1357               index=icl;
1358             }
1359     }
1360     else
1361     {
1362       printf("Error: please specify your cut criteria\n");
1363       printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1364       printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1365       return -1;
1366     }
1367   }//cluster loop
1368   return index;
1369 }
1370
1371 //________________________________________________________________________________
1372 Bool_t  AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Float_t &tmpEta, Float_t &tmpPhi)
1373 {
1374   //
1375   //Return the residual by extrapolating a track to a cluster
1376   //
1377   if(!trkParam || !cluster) return kFALSE;
1378   Float_t clsPos[3];
1379   Double_t trkPos[3];
1380   cluster->GetPosition(clsPos); //Has been recalculated
1381   TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
1382   Double_t alpha =  ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
1383   vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
1384   trkParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
1385   if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), fMass, fStep,kFALSE, 0.8, -1)) return kFALSE; 
1386   trkParam->GetXYZ(trkPos); //Get the extrapolated global position
1387
1388   TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1389   TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
1390
1391   Float_t clsPhi = (Float_t)clsPosVec.Phi();
1392   if(clsPhi<0) clsPhi+=2*TMath::Pi();
1393   Float_t trkPhi = (Float_t)trkPosVec.Phi();
1394   if(trkPhi<0) trkPhi+=2*TMath::Pi();
1395   tmpPhi = clsPhi-trkPhi;  // track cluster matching
1396   tmpEta = clsPosVec.Eta()-trkPosVec.Eta();  // track cluster matching
1397
1398   return kTRUE;
1399 }
1400
1401 //________________________________________________________________________________
1402 void AliEMCALRecoUtils::GetMatchedResiduals(Int_t clsIndex, Float_t &dEta, Float_t &dPhi)
1403 {
1404   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1405   //Get the residuals dEta and dPhi for this cluster to the closest track
1406   //Works with ESDs and AODs
1407
1408   if( FindMatchedPosForCluster(clsIndex) >= 999 )
1409   {
1410     AliDebug(2,"No matched tracks found!\n");
1411     dEta=999.;
1412     dPhi=999.;
1413     return;
1414   }
1415   dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
1416   dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
1417 }
1418 //________________________________________________________________________________
1419 void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi)
1420 {
1421   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1422   //Get the residuals dEta and dPhi for this track to the closest cluster
1423   //Works with ESDs and AODs
1424
1425   if( FindMatchedPosForTrack(trkIndex) >= 999 )
1426   {
1427     AliDebug(2,"No matched cluster found!\n");
1428     dEta=999.;
1429     dPhi=999.;
1430     return;
1431   }
1432   dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
1433   dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
1434 }
1435
1436 //__________________________________________________________
1437 Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
1438 {
1439   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1440   //Get the index of matched track to this cluster
1441   //Works with ESDs and AODs
1442   
1443   if(IsClusterMatched(clsIndex))
1444     return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
1445   else 
1446     return -1; 
1447 }
1448
1449 //__________________________________________________________
1450 Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
1451 {
1452   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1453   //Get the index of matched cluster to this track
1454   //Works with ESDs and AODs
1455   
1456   if(IsTrackMatched(trkIndex))
1457     return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
1458   else 
1459     return -1; 
1460 }
1461
1462 //__________________________________________________
1463 Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex) const
1464 {
1465   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1466   //Returns if the cluster has a match
1467   if(FindMatchedPosForCluster(clsIndex) < 999) 
1468     return kTRUE;
1469   else
1470     return kFALSE;
1471 }
1472
1473 //__________________________________________________
1474 Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex) const 
1475 {
1476   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1477   //Returns if the track has a match
1478   if(FindMatchedPosForTrack(trkIndex) < 999) 
1479     return kTRUE;
1480   else
1481     return kFALSE;
1482 }
1483
1484 //__________________________________________________________
1485 UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
1486 {
1487   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1488   //Returns the position of the match in the fMatchedClusterIndex array
1489   Float_t tmpR = fCutR;
1490   UInt_t pos = 999;
1491   
1492   for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
1493   {
1494     if(fMatchedClusterIndex->At(i)==clsIndex)
1495       {
1496         Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
1497         if(r<tmpR)
1498           {
1499             pos=i;
1500             tmpR=r;
1501             AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
1502           }
1503       }
1504   }
1505   return pos;
1506 }
1507
1508 //__________________________________________________________
1509 UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
1510 {
1511   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1512   //Returns the position of the match in the fMatchedTrackIndex array
1513   Float_t tmpR = fCutR;
1514   UInt_t pos = 999;
1515   
1516   for(Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
1517   {
1518     if(fMatchedTrackIndex->At(i)==trkIndex)
1519       {
1520         Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
1521         if(r<tmpR)
1522           {
1523             pos=i;
1524             tmpR=r;
1525             AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
1526           }
1527       }
1528   }
1529   return pos;
1530 }
1531
1532 //__________________________________________________________
1533 Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster, AliEMCALGeometry *geom, AliVCaloCells* cells)
1534 {
1535   // check if the cluster survives some quality cut
1536   //
1537   //
1538   Bool_t isGood=kTRUE;
1539   if(!cluster || !cluster->IsEMCAL()) return kFALSE;
1540   if(ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
1541   if(!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
1542   if(fRejectExoticCluster && IsExoticCluster(cluster)) return kFALSE;
1543
1544   return isGood;
1545 }
1546
1547 //__________________________________________________________
1548 Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
1549 {
1550   // Given a esd track, return whether the track survive all the cuts
1551
1552   // The different quality parameter are first
1553   // retrieved from the track. then it is found out what cuts the
1554   // track did not survive and finally the cuts are imposed.
1555
1556   UInt_t status = esdTrack->GetStatus();
1557
1558   Int_t nClustersITS = esdTrack->GetITSclusters(0);
1559   Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
1560
1561   Float_t chi2PerClusterITS = -1;
1562   Float_t chi2PerClusterTPC = -1;
1563   if (nClustersITS!=0)
1564     chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1565   if (nClustersTPC!=0) 
1566     chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
1567
1568
1569   //DCA cuts
1570   Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1571   //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
1572   SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
1573
1574
1575   Float_t b[2];
1576   Float_t bCov[3];
1577   esdTrack->GetImpactParameters(b,bCov);
1578   if (bCov[0]<=0 || bCov[2]<=0) {
1579     AliDebug(1, "Estimated b resolution lower or equal zero!");
1580     bCov[0]=0; bCov[2]=0;
1581   }
1582
1583   Float_t dcaToVertexXY = b[0];
1584   Float_t dcaToVertexZ = b[1];
1585   Float_t dcaToVertex = -1;
1586
1587   if (fCutDCAToVertex2D)
1588     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1589   else
1590     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1591     
1592   // cut the track?
1593   
1594   Bool_t cuts[kNCuts];
1595   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1596   
1597   // track quality cuts
1598   if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1599     cuts[0]=kTRUE;
1600   if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1601     cuts[1]=kTRUE;
1602   if (nClustersTPC<fCutMinNClusterTPC)
1603     cuts[2]=kTRUE;
1604   if (nClustersITS<fCutMinNClusterITS) 
1605     cuts[3]=kTRUE;
1606   if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC) 
1607     cuts[4]=kTRUE; 
1608   if (chi2PerClusterITS>fCutMaxChi2PerClusterITS) 
1609     cuts[5]=kTRUE;  
1610   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1611     cuts[6]=kTRUE;
1612   if (fCutDCAToVertex2D && dcaToVertex > 1)
1613     cuts[7] = kTRUE;
1614   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1615     cuts[8] = kTRUE;
1616   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1617     cuts[9] = kTRUE;
1618
1619   //Require at least one SPD point + anything else in ITS
1620   if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1621     cuts[10] = kTRUE;
1622
1623   Bool_t cut=kFALSE;
1624   for (Int_t i=0; i<kNCuts; i++) 
1625     if (cuts[i]) {cut = kTRUE;}
1626
1627     // cut the track
1628   if (cut) 
1629     return kFALSE;
1630   else 
1631     return kTRUE;
1632 }
1633 //__________________________________________________
1634 void AliEMCALRecoUtils::InitTrackCuts()
1635 {
1636   //Intilize the track cut criteria
1637   //By default these cuts are set according to AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()
1638   //Also you can customize the cuts using the setters
1639   
1640   switch (fTrackCutsType)
1641     {
1642     case kTPCOnlyCut:
1643       {
1644         AliInfo(Form("Track cuts for matching: GetStandardTPCOnlyTrackCuts()"));
1645         //TPC
1646         SetMinNClustersTPC(70);
1647         SetMaxChi2PerClusterTPC(4);
1648         SetAcceptKinkDaughters(kFALSE);
1649         SetRequireTPCRefit(kFALSE);
1650
1651         //ITS
1652         SetRequireITSRefit(kFALSE);
1653         SetMaxDCAToVertexZ(3.2);
1654         SetMaxDCAToVertexXY(2.4);
1655         SetDCAToVertex2D(kTRUE);
1656
1657         break;
1658       }
1659     
1660     case kGlobalCut:
1661       {
1662         AliInfo(Form("Track cuts for matching: GetStandardITSTPCTrackCuts2010(kTURE)"));
1663         //TPC
1664         SetMinNClustersTPC(70);
1665         SetMaxChi2PerClusterTPC(4);
1666         SetAcceptKinkDaughters(kFALSE);
1667         SetRequireTPCRefit(kTRUE);
1668
1669         //ITS
1670         SetRequireITSRefit(kTRUE);
1671         SetMaxDCAToVertexZ(2);
1672         SetMaxDCAToVertexXY();
1673         SetDCAToVertex2D(kFALSE);
1674
1675         break;
1676       }
1677     }
1678 }
1679
1680 //___________________________________________________
1681 void AliEMCALRecoUtils::Print(const Option_t *) const 
1682 {
1683   // Print Parameters
1684   
1685   printf("AliEMCALRecoUtils Settings: \n");
1686   printf("Misalignment shifts\n");
1687   for(Int_t i=0; i<5; i++) printf("\t sector %d, traslation (x,y,z)=(%f,%f,%f), rotation (x,y,z)=(%f,%f,%f)\n",i, 
1688                                   fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1689                                   fMisalRotShift[i*3],  fMisalRotShift[i*3+1],  fMisalRotShift[i*3+2]   );
1690   printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1691   for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
1692   
1693   printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
1694
1695   printf("Matching criteria: ");
1696   if(fCutEtaPhiSum)
1697     {
1698       printf("sqrt(dEta^2+dPhi^2)<%2.2f\n",fCutR);
1699     }
1700   else if(fCutEtaPhiSeparate)
1701     {
1702       printf("dEta<%2.2f, dPhi<%2.2f\n",fCutEta,fCutPhi);
1703     }
1704   else
1705     {
1706       printf("Error\n");
1707       printf("please specify your cut criteria\n");
1708       printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1709       printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1710     }
1711
1712   printf("Mass hypothesis = %2.3f [GeV/c^2], extrapolation step = %2.2f[cm]\n",fMass,fStep);
1713
1714   printf("Track cuts: \n");
1715   printf("Minimum track pT: %1.2f\n",fCutMinTrackPt);
1716   printf("AOD track selection mask: %d\n",fAODFilterMask);
1717   printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1718   printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1719   printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1720   printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1721   printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
1722
1723 }
1724
1725 //_____________________________________________________________________
1726 void AliEMCALRecoUtils::SetTimeDependentCorrections(Int_t runnumber){
1727   //Get EMCAL time dependent corrections from file and put them in the recalibration histograms
1728   //Do it only once and only if it is requested
1729   
1730   if(!fUseTimeCorrectionFactors) return;
1731   if(fTimeCorrectionFactorsSet)  return;
1732   
1733   printf("AliEMCALRecoUtils::GetTimeDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber);
1734  
1735   AliEMCALCalibTimeDepCorrection  *corr =  new AliEMCALCalibTimeDepCorrection();
1736   corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber));
1737   
1738   SwitchOnRecalibration();
1739   for(Int_t ism = 0; ism < 4; ism++){
1740     for(Int_t icol = 0; icol < 48; icol++){
1741       for(Int_t irow = 0; irow < 24; irow++){
1742         Float_t orgRecalFactor = GetEMCALChannelRecalibrationFactors(ism)->GetBinContent(icol,irow);
1743         Float_t newRecalFactor = orgRecalFactor*corr->GetCorrection(ism, icol,irow,0);
1744         GetEMCALChannelRecalibrationFactors(ism)->SetBinContent(icol,irow,newRecalFactor);
1745         //printf("ism %d, icol %d, irow %d, corrections : org %f, time dep %f, final %f (org*time %f)\n",ism, icol, irow, 
1746         //        orgRecalFactor, corr->GetCorrection(ism, icol,irow,0),
1747         //       (GetEMCALChannelRecalibrationFactors(ism))->GetBinContent(icol,irow),newRecalFactor);
1748       }
1749     }
1750   }
1751    fTimeCorrectionFactorsSet = kTRUE;
1752 }
1753