In case of track matching, extrapolate only from inner track param
[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       trackParam =  const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam());
1176     }
1177     
1178     //If the input event is AOD, the starting point for extrapolation is at vertex
1179     //AOD tracks are selected according to its bit.
1180     else if(aodevent)
1181     {
1182       AliAODTrack *aodTrack = aodevent->GetTrack(itr);
1183       if(!aodTrack) continue;
1184       if(!aodTrack->TestFilterMask(fAODFilterMask)) continue; //Select AOD tracks that fulfill GetStandardITSTPCTrackCuts2010()
1185       if(aodTrack->Pt()<fCutMinTrackPt) continue;
1186       Double_t pos[3],mom[3];
1187       aodTrack->GetXYZ(pos);
1188       aodTrack->GetPxPyPz(mom);
1189       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()));
1190       trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
1191     }
1192     
1193     //Return if the input data is not "AOD" or "ESD"
1194     else
1195     {
1196       printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n");
1197       return;
1198     }
1199     
1200     if(!trackParam) continue;
1201     
1202     Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
1203     Int_t index = -1;
1204     if(!clusterArr){// get clusters from event
1205       for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1206       {
1207         AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1208         if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1209         AliExternalTrackParam trkPamTmp(*trackParam);//Retrieve the starting point every time before the extrapolation  
1210         Float_t tmpEta=-999, tmpPhi=-999;
1211         if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
1212         if(fCutEtaPhiSum)
1213         {
1214           Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1215           if(tmpR<dRMax)
1216           {
1217             dRMax=tmpR;
1218             dEtaMax=tmpEta;
1219             dPhiMax=tmpPhi;
1220             index=icl;
1221           }
1222         }
1223         else if(fCutEtaPhiSeparate)
1224         {
1225           if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
1226           {
1227             dEtaMax = tmpEta;
1228             dPhiMax = tmpPhi;
1229             index=icl;
1230           }
1231         }
1232         else
1233         {
1234           printf("Error: please specify your cut criteria\n");
1235           printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1236           printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1237           if(aodevent && trackParam) delete trackParam;
1238           return;
1239         }
1240       }//cluster loop
1241     } 
1242     else { // external cluster array, not from ESD event
1243       for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
1244       {
1245         AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
1246         if(!cluster){ 
1247           AliInfo("Cluster not found!!!");
1248           continue;
1249         }
1250         if(!cluster->IsEMCAL()) continue;
1251         AliExternalTrackParam trkPamTmp (*trackParam);//Retrieve the starting point every time before the extrapolation
1252         Float_t tmpEta=-999, tmpPhi=-999;
1253         if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
1254         if(fCutEtaPhiSum)
1255         {
1256           Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1257           if(tmpR<dRMax)
1258           {
1259             dRMax=tmpR;
1260             dEtaMax=tmpEta;
1261             dPhiMax=tmpPhi;
1262             index=icl;
1263           }
1264         }
1265         else if(fCutEtaPhiSeparate)
1266         {
1267           if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
1268           {
1269             dEtaMax = tmpEta;
1270             dPhiMax = tmpPhi;
1271             index=icl;
1272           }
1273         }
1274         else
1275         {
1276           printf("Error: please specify your cut criteria\n");
1277           printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1278           printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1279           if(aodevent && trackParam) delete trackParam;
1280           return;
1281         }
1282       }//cluster loop
1283     }// external list of clusters
1284     
1285     if(index>-1)
1286     {
1287       fMatchedTrackIndex  ->AddAt(itr,matched);
1288       fMatchedClusterIndex->AddAt(index,matched);
1289       fResidualEta          ->AddAt(dEtaMax,matched);
1290       fResidualPhi          ->AddAt(dPhiMax,matched);
1291       matched++;
1292     }
1293     if(aodevent && trackParam) delete trackParam;
1294   }//track loop
1295   
1296   AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1297   
1298   fMatchedTrackIndex  ->Set(matched);
1299   fMatchedClusterIndex->Set(matched);
1300   fResidualPhi          ->Set(matched);
1301   fResidualEta          ->Set(matched);
1302 }
1303
1304 //________________________________________________________________________________
1305 Int_t AliEMCALRecoUtils::FindMatchedCluster(AliESDtrack *track, AliVEvent *event, AliEMCALGeometry *geom)
1306 {
1307   //
1308   // This function returns the index of matched cluster to input track
1309   // Returns -1 if no match is found
1310   
1311   
1312   Float_t dRMax = fCutR, dEtaMax = fCutEta, dPhiMax = fCutPhi;
1313   Int_t index = -1;
1314   
1315   AliExternalTrackParam *trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
1316   
1317   if(!trackParam) return index;   
1318   for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1319   {
1320     AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1321     if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;   
1322     AliExternalTrackParam trkPamTmp (*trackParam);//Retrieve the starting point every time before the extrapolation
1323     Float_t tmpEta=-999, tmpPhi=-999;
1324     if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
1325     if(fCutEtaPhiSum)
1326     {
1327       Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1328       if(tmpR<dRMax)
1329             {
1330               dRMax=tmpR;
1331               dEtaMax=tmpEta;
1332               dPhiMax=tmpPhi;
1333               index=icl;
1334             }
1335     }
1336     else if(fCutEtaPhiSeparate)
1337     {
1338       if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
1339             {
1340               dEtaMax = tmpEta;
1341               dPhiMax = tmpPhi;
1342               index=icl;
1343             }
1344     }
1345     else
1346     {
1347       printf("Error: please specify your cut criteria\n");
1348       printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1349       printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1350       return -1;
1351     }
1352   }//cluster loop
1353   return index;
1354 }
1355
1356 //________________________________________________________________________________
1357 Bool_t  AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Float_t &tmpEta, Float_t &tmpPhi)
1358 {
1359   //
1360   //Return the residual by extrapolating a track to a cluster
1361   //
1362   if(!trkParam || !cluster) return kFALSE;
1363   Float_t clsPos[3];
1364   Double_t trkPos[3];
1365   cluster->GetPosition(clsPos); //Has been recalculated
1366   TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
1367   Double_t alpha =  ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
1368   vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
1369   trkParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
1370   if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), fMass, fStep,kFALSE, 0.8, -1)) return kFALSE; 
1371   trkParam->GetXYZ(trkPos); //Get the extrapolated global position
1372
1373   TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1374   TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
1375
1376   Float_t clsPhi = (Float_t)clsPosVec.Phi();
1377   if(clsPhi<0) clsPhi+=2*TMath::Pi();
1378   Float_t trkPhi = (Float_t)trkPosVec.Phi();
1379   if(trkPhi<0) trkPhi+=2*TMath::Pi();
1380   tmpPhi = clsPhi-trkPhi;  // track cluster matching
1381   tmpEta = clsPosVec.Eta()-trkPosVec.Eta();  // track cluster matching
1382
1383   return kTRUE;
1384 }
1385
1386 //________________________________________________________________________________
1387 void AliEMCALRecoUtils::GetMatchedResiduals(Int_t clsIndex, Float_t &dEta, Float_t &dPhi)
1388 {
1389   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1390   //Get the residuals dEta and dPhi for this cluster to the closest track
1391   //Works with ESDs and AODs
1392
1393   if( FindMatchedPosForCluster(clsIndex) >= 999 )
1394   {
1395     AliDebug(2,"No matched tracks found!\n");
1396     dEta=999.;
1397     dPhi=999.;
1398     return;
1399   }
1400   dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
1401   dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
1402 }
1403 //________________________________________________________________________________
1404 void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi)
1405 {
1406   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1407   //Get the residuals dEta and dPhi for this track to the closest cluster
1408   //Works with ESDs and AODs
1409
1410   if( FindMatchedPosForTrack(trkIndex) >= 999 )
1411   {
1412     AliDebug(2,"No matched cluster found!\n");
1413     dEta=999.;
1414     dPhi=999.;
1415     return;
1416   }
1417   dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
1418   dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
1419 }
1420
1421 //__________________________________________________________
1422 Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
1423 {
1424   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1425   //Get the index of matched track to this cluster
1426   //Works with ESDs and AODs
1427   
1428   if(IsClusterMatched(clsIndex))
1429     return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
1430   else 
1431     return -1; 
1432 }
1433
1434 //__________________________________________________________
1435 Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
1436 {
1437   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1438   //Get the index of matched cluster to this track
1439   //Works with ESDs and AODs
1440   
1441   if(IsTrackMatched(trkIndex))
1442     return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
1443   else 
1444     return -1; 
1445 }
1446
1447 //__________________________________________________
1448 Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex) const
1449 {
1450   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1451   //Returns if the cluster has a match
1452   if(FindMatchedPosForCluster(clsIndex) < 999) 
1453     return kTRUE;
1454   else
1455     return kFALSE;
1456 }
1457
1458 //__________________________________________________
1459 Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex) const 
1460 {
1461   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1462   //Returns if the track has a match
1463   if(FindMatchedPosForTrack(trkIndex) < 999) 
1464     return kTRUE;
1465   else
1466     return kFALSE;
1467 }
1468
1469 //__________________________________________________________
1470 UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
1471 {
1472   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1473   //Returns the position of the match in the fMatchedClusterIndex array
1474   Float_t tmpR = fCutR;
1475   UInt_t pos = 999;
1476   
1477   for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
1478   {
1479     if(fMatchedClusterIndex->At(i)==clsIndex)
1480       {
1481         Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
1482         if(r<tmpR)
1483           {
1484             pos=i;
1485             tmpR=r;
1486             AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
1487           }
1488       }
1489   }
1490   return pos;
1491 }
1492
1493 //__________________________________________________________
1494 UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
1495 {
1496   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1497   //Returns the position of the match in the fMatchedTrackIndex array
1498   Float_t tmpR = fCutR;
1499   UInt_t pos = 999;
1500   
1501   for(Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
1502   {
1503     if(fMatchedTrackIndex->At(i)==trkIndex)
1504       {
1505         Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
1506         if(r<tmpR)
1507           {
1508             pos=i;
1509             tmpR=r;
1510             AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
1511           }
1512       }
1513   }
1514   return pos;
1515 }
1516
1517 //__________________________________________________________
1518 Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster, AliEMCALGeometry *geom, AliVCaloCells* cells)
1519 {
1520   // check if the cluster survives some quality cut
1521   //
1522   //
1523   Bool_t isGood=kTRUE;
1524   if(!cluster || !cluster->IsEMCAL()) return kFALSE;
1525   if(ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
1526   if(!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
1527   if(fRejectExoticCluster && IsExoticCluster(cluster)) return kFALSE;
1528
1529   return isGood;
1530 }
1531
1532 //__________________________________________________________
1533 Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
1534 {
1535   // Given a esd track, return whether the track survive all the cuts
1536
1537   // The different quality parameter are first
1538   // retrieved from the track. then it is found out what cuts the
1539   // track did not survive and finally the cuts are imposed.
1540
1541   UInt_t status = esdTrack->GetStatus();
1542
1543   Int_t nClustersITS = esdTrack->GetITSclusters(0);
1544   Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
1545
1546   Float_t chi2PerClusterITS = -1;
1547   Float_t chi2PerClusterTPC = -1;
1548   if (nClustersITS!=0)
1549     chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1550   if (nClustersTPC!=0) 
1551     chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
1552
1553
1554   //DCA cuts
1555   Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1556   //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
1557   SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
1558
1559
1560   Float_t b[2];
1561   Float_t bCov[3];
1562   esdTrack->GetImpactParameters(b,bCov);
1563   if (bCov[0]<=0 || bCov[2]<=0) {
1564     AliDebug(1, "Estimated b resolution lower or equal zero!");
1565     bCov[0]=0; bCov[2]=0;
1566   }
1567
1568   Float_t dcaToVertexXY = b[0];
1569   Float_t dcaToVertexZ = b[1];
1570   Float_t dcaToVertex = -1;
1571
1572   if (fCutDCAToVertex2D)
1573     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1574   else
1575     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1576     
1577   // cut the track?
1578   
1579   Bool_t cuts[kNCuts];
1580   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1581   
1582   // track quality cuts
1583   if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1584     cuts[0]=kTRUE;
1585   if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1586     cuts[1]=kTRUE;
1587   if (nClustersTPC<fCutMinNClusterTPC)
1588     cuts[2]=kTRUE;
1589   if (nClustersITS<fCutMinNClusterITS) 
1590     cuts[3]=kTRUE;
1591   if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC) 
1592     cuts[4]=kTRUE; 
1593   if (chi2PerClusterITS>fCutMaxChi2PerClusterITS) 
1594     cuts[5]=kTRUE;  
1595   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1596     cuts[6]=kTRUE;
1597   if (fCutDCAToVertex2D && dcaToVertex > 1)
1598     cuts[7] = kTRUE;
1599   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1600     cuts[8] = kTRUE;
1601   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1602     cuts[9] = kTRUE;
1603
1604   //Require at least one SPD point + anything else in ITS
1605   if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1606     cuts[10] = kTRUE;
1607
1608   Bool_t cut=kFALSE;
1609   for (Int_t i=0; i<kNCuts; i++) 
1610     if (cuts[i]) {cut = kTRUE;}
1611
1612     // cut the track
1613   if (cut) 
1614     return kFALSE;
1615   else 
1616     return kTRUE;
1617 }
1618 //__________________________________________________
1619 void AliEMCALRecoUtils::InitTrackCuts()
1620 {
1621   //Intilize the track cut criteria
1622   //By default these cuts are set according to AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()
1623   //Also you can customize the cuts using the setters
1624   
1625   switch (fTrackCutsType)
1626     {
1627     case kTPCOnlyCut:
1628       {
1629         AliInfo(Form("Track cuts for matching: GetStandardTPCOnlyTrackCuts()"));
1630         //TPC
1631         SetMinNClustersTPC(70);
1632         SetMaxChi2PerClusterTPC(4);
1633         SetAcceptKinkDaughters(kFALSE);
1634         SetRequireTPCRefit(kFALSE);
1635
1636         //ITS
1637         SetRequireITSRefit(kFALSE);
1638         SetMaxDCAToVertexZ(3.2);
1639         SetMaxDCAToVertexXY(2.4);
1640         SetDCAToVertex2D(kTRUE);
1641
1642         break;
1643       }
1644     
1645     case kGlobalCut:
1646       {
1647         AliInfo(Form("Track cuts for matching: GetStandardITSTPCTrackCuts2010(kTURE)"));
1648         //TPC
1649         SetMinNClustersTPC(70);
1650         SetMaxChi2PerClusterTPC(4);
1651         SetAcceptKinkDaughters(kFALSE);
1652         SetRequireTPCRefit(kTRUE);
1653
1654         //ITS
1655         SetRequireITSRefit(kTRUE);
1656         SetMaxDCAToVertexZ(2);
1657         SetMaxDCAToVertexXY();
1658         SetDCAToVertex2D(kFALSE);
1659
1660         break;
1661       }
1662     }
1663 }
1664
1665 //___________________________________________________
1666 void AliEMCALRecoUtils::Print(const Option_t *) const 
1667 {
1668   // Print Parameters
1669   
1670   printf("AliEMCALRecoUtils Settings: \n");
1671   printf("Misalignment shifts\n");
1672   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, 
1673                                   fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1674                                   fMisalRotShift[i*3],  fMisalRotShift[i*3+1],  fMisalRotShift[i*3+2]   );
1675   printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1676   for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
1677   
1678   printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
1679
1680   printf("Matching criteria: ");
1681   if(fCutEtaPhiSum)
1682     {
1683       printf("sqrt(dEta^2+dPhi^2)<%2.2f\n",fCutR);
1684     }
1685   else if(fCutEtaPhiSeparate)
1686     {
1687       printf("dEta<%2.2f, dPhi<%2.2f\n",fCutEta,fCutPhi);
1688     }
1689   else
1690     {
1691       printf("Error\n");
1692       printf("please specify your cut criteria\n");
1693       printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1694       printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1695     }
1696
1697   printf("Mass hypothesis = %2.3f [GeV/c^2], extrapolation step = %2.2f[cm]\n",fMass,fStep);
1698
1699   printf("Track cuts: \n");
1700   printf("Minimum track pT: %1.2f\n",fCutMinTrackPt);
1701   printf("AOD track selection mask: %d\n",fAODFilterMask);
1702   printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1703   printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1704   printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1705   printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1706   printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
1707
1708 }
1709
1710 //_____________________________________________________________________
1711 void AliEMCALRecoUtils::SetTimeDependentCorrections(Int_t runnumber){
1712   //Get EMCAL time dependent corrections from file and put them in the recalibration histograms
1713   //Do it only once and only if it is requested
1714   
1715   if(!fUseTimeCorrectionFactors) return;
1716   if(fTimeCorrectionFactorsSet)  return;
1717   
1718   printf("AliEMCALRecoUtils::GetTimeDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber);
1719  
1720   AliEMCALCalibTimeDepCorrection  *corr =  new AliEMCALCalibTimeDepCorrection();
1721   corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber));
1722   
1723   SwitchOnRecalibration();
1724   for(Int_t ism = 0; ism < 4; ism++){
1725     for(Int_t icol = 0; icol < 48; icol++){
1726       for(Int_t irow = 0; irow < 24; irow++){
1727         Float_t orgRecalFactor = GetEMCALChannelRecalibrationFactors(ism)->GetBinContent(icol,irow);
1728         Float_t newRecalFactor = orgRecalFactor*corr->GetCorrection(ism, icol,irow,0);
1729         GetEMCALChannelRecalibrationFactors(ism)->SetBinContent(icol,irow,newRecalFactor);
1730         //printf("ism %d, icol %d, irow %d, corrections : org %f, time dep %f, final %f (org*time %f)\n",ism, icol, irow, 
1731         //        orgRecalFactor, corr->GetCorrection(ism, icol,irow,0),
1732         //       (GetEMCALChannelRecalibrationFactors(ism))->GetBinContent(icol,irow),newRecalFactor);
1733       }
1734     }
1735   }
1736    fTimeCorrectionFactorsSet = kTRUE;
1737 }
1738