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