459d70c3905979f36a239c99aca58be10d6a2e2e
[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 "AliLog.h"
46 #include "AliPID.h"
47 #include "AliESDEvent.h"
48 #include "AliAODEvent.h"
49 #include "AliESDtrack.h"
50 #include "AliAODTrack.h"
51 #include "AliExternalTrackParam.h"
52 #include "AliESDfriendTrack.h"
53 #include "AliTrackerBase.h"
54
55 // EMCAL includes
56 #include "AliEMCALRecoUtils.h"
57 #include "AliEMCALGeometry.h"
58 #include "AliTrackerBase.h"
59 #include "AliEMCALCalibTimeDepCorrection.h" // Run dependent
60 #include "AliEMCALPIDUtils.h"
61
62
63 ClassImp(AliEMCALRecoUtils)
64   
65 //_____________________________________
66 AliEMCALRecoUtils::AliEMCALRecoUtils():
67   fParticleType(0),                       fPosAlgo(0),                            fW0(0), 
68   fNonLinearityFunction(0),               fNonLinearThreshold(0),
69   fSmearClusterEnergy(kFALSE),            fRandom(),
70   fCellsRecalibrated(kFALSE),             fRecalibration(kFALSE),                 fEMCALRecalibrationFactors(),
71   fTimeRecalibration(kFALSE),             fEMCALTimeRecalibrationFactors(),
72   fUseRunCorrectionFactors(kFALSE),       fRunCorrectionFactorsSet(kFALSE),
73   fRemoveBadChannels(kFALSE),             fRecalDistToBadChannels(kFALSE),        fEMCALBadChannelMap(),
74   fNCellsFromEMCALBorder(0),              fNoEMCALBorderAtEta0(kTRUE),
75   fRejectExoticCluster(kFALSE),           fRejectExoticCells(kFALSE), 
76   fExoticCellFraction(0),                 fExoticCellDiffTime(0),                 fExoticCellMinAmplitude(0),
77   fPIDUtils(),                            fAODFilterMask(0),
78   fMatchedTrackIndex(0x0),                fMatchedClusterIndex(0x0), 
79   fResidualEta(0x0), fResidualPhi(0x0),   fCutEtaPhiSum(kFALSE),                  fCutEtaPhiSeparate(kFALSE), 
80   fCutR(0),                               fCutEta(0),                             fCutPhi(0),
81   fClusterWindow(0),                      fMass(0),                           
82   fStepSurface(0),                        fStepCluster(0),
83   fTrackCutsType(0),                      fCutMinTrackPt(0),                      fCutMinNClusterTPC(0), 
84   fCutMinNClusterITS(0),                  fCutMaxChi2PerClusterTPC(0),            fCutMaxChi2PerClusterITS(0),
85   fCutRequireTPCRefit(kFALSE),            fCutRequireITSRefit(kFALSE),            fCutAcceptKinkDaughters(kFALSE),
86   fCutMaxDCAToVertexXY(0),                fCutMaxDCAToVertexZ(0),                 fCutDCAToVertex2D(kFALSE)
87 {
88 //
89   // Constructor.
90   // Initialize all constant values which have to be used
91   // during Reco algorithm execution
92   //
93   
94   // Init parameters
95   InitParameters();
96   
97   //Track matching
98   fMatchedTrackIndex     = new TArrayI();
99   fMatchedClusterIndex   = new TArrayI();
100   fResidualPhi           = new TArrayF();
101   fResidualEta           = new TArrayF();
102   fPIDUtils              = new AliEMCALPIDUtils();
103
104   InitTrackCuts();
105 }
106
107 //______________________________________________________________________
108 AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco) 
109 : TNamed(reco), 
110   fParticleType(reco.fParticleType),                         fPosAlgo(reco.fPosAlgo),     fW0(reco.fW0),
111   fNonLinearityFunction(reco.fNonLinearityFunction),         fNonLinearThreshold(reco.fNonLinearThreshold),
112   fSmearClusterEnergy(reco.fSmearClusterEnergy),             fRandom(),
113   fCellsRecalibrated(reco.fCellsRecalibrated),
114   fRecalibration(reco.fRecalibration),                       fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
115   fTimeRecalibration(reco.fTimeRecalibration),               fEMCALTimeRecalibrationFactors(reco.fEMCALTimeRecalibrationFactors),
116   fUseRunCorrectionFactors(reco.fUseRunCorrectionFactors),   fRunCorrectionFactorsSet(reco.fRunCorrectionFactorsSet),
117   fRemoveBadChannels(reco.fRemoveBadChannels),               fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
118   fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
119   fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),       fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
120   fRejectExoticCluster(reco.fRejectExoticCluster),           fRejectExoticCells(reco.fRejectExoticCells), 
121   fExoticCellFraction(reco.fExoticCellFraction),             fExoticCellDiffTime(reco.fExoticCellDiffTime),               
122   fExoticCellMinAmplitude(reco.fExoticCellMinAmplitude),
123   fPIDUtils(reco.fPIDUtils),                                 fAODFilterMask(reco.fAODFilterMask),
124   fMatchedTrackIndex(  reco.fMatchedTrackIndex?  new TArrayI(*reco.fMatchedTrackIndex):0x0),
125   fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
126   fResidualEta(        reco.fResidualEta?        new TArrayF(*reco.fResidualEta):0x0),
127   fResidualPhi(        reco.fResidualPhi?        new TArrayF(*reco.fResidualPhi):0x0),
128   fCutEtaPhiSum(reco.fCutEtaPhiSum),                         fCutEtaPhiSeparate(reco.fCutEtaPhiSeparate), 
129   fCutR(reco.fCutR),        fCutEta(reco.fCutEta),           fCutPhi(reco.fCutPhi),
130   fClusterWindow(reco.fClusterWindow),
131   fMass(reco.fMass),        fStepSurface(reco.fStepSurface), fStepCluster(reco.fStepCluster),
132   fTrackCutsType(reco.fTrackCutsType),                       fCutMinTrackPt(reco.fCutMinTrackPt), 
133   fCutMinNClusterTPC(reco.fCutMinNClusterTPC),               fCutMinNClusterITS(reco.fCutMinNClusterITS), 
134   fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC),   fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS),
135   fCutRequireTPCRefit(reco.fCutRequireTPCRefit),             fCutRequireITSRefit(reco.fCutRequireITSRefit),
136   fCutAcceptKinkDaughters(reco.fCutAcceptKinkDaughters),     fCutMaxDCAToVertexXY(reco.fCutMaxDCAToVertexXY),    
137   fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ),             fCutDCAToVertex2D(reco.fCutDCAToVertex2D)
138 {
139   //Copy ctor
140   
141   for(Int_t i = 0; i < 15 ; i++) { fMisalRotShift[i]      = reco.fMisalRotShift[i]      ; 
142                                    fMisalTransShift[i]    = reco.fMisalTransShift[i]    ; } 
143   for(Int_t i = 0; i < 7  ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
144   for(Int_t i = 0; i < 3  ; i++) { fSmearClusterParam[i]  = reco.fSmearClusterParam[i]  ; }
145
146 }
147
148
149 //______________________________________________________________________
150 AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco) 
151 {
152   //Assignment operator
153   
154   if(this == &reco)return *this;
155   ((TNamed *)this)->operator=(reco);
156
157   for(Int_t i = 0; i < 15 ; i++) { fMisalTransShift[i]    = reco.fMisalTransShift[i]    ; 
158                                    fMisalRotShift[i]      = reco.fMisalRotShift[i]      ; }
159   for(Int_t i = 0; i < 7  ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
160   for(Int_t i = 0; i < 3  ; i++) { fSmearClusterParam[i]  = reco.fSmearClusterParam[i]  ; }   
161   
162   fParticleType              = reco.fParticleType;
163   fPosAlgo                   = reco.fPosAlgo; 
164   fW0                        = reco.fW0;
165   
166   fNonLinearityFunction      = reco.fNonLinearityFunction;
167   fNonLinearThreshold        = reco.fNonLinearThreshold;
168   fSmearClusterEnergy        = reco.fSmearClusterEnergy;
169
170   fCellsRecalibrated         = reco.fCellsRecalibrated;
171   fRecalibration             = reco.fRecalibration;
172   fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
173
174   fTimeRecalibration             = reco.fTimeRecalibration;
175   fEMCALTimeRecalibrationFactors = reco.fEMCALTimeRecalibrationFactors;
176
177   fUseRunCorrectionFactors   = reco.fUseRunCorrectionFactors;
178   fRunCorrectionFactorsSet   = reco.fRunCorrectionFactorsSet;
179   
180   fRemoveBadChannels         = reco.fRemoveBadChannels;
181   fRecalDistToBadChannels    = reco.fRecalDistToBadChannels;
182   fEMCALBadChannelMap        = reco.fEMCALBadChannelMap;
183   
184   fNCellsFromEMCALBorder     = reco.fNCellsFromEMCALBorder;
185   fNoEMCALBorderAtEta0       = reco.fNoEMCALBorderAtEta0;
186   
187   fRejectExoticCluster       = reco.fRejectExoticCluster;           
188   fRejectExoticCells         = reco.fRejectExoticCells; 
189   fExoticCellFraction        = reco.fExoticCellFraction;
190   fExoticCellDiffTime        = reco.fExoticCellDiffTime;              
191   fExoticCellMinAmplitude    = reco.fExoticCellMinAmplitude;
192   
193   fPIDUtils                  = reco.fPIDUtils;
194
195   fAODFilterMask             = reco.fAODFilterMask;
196   
197   fCutEtaPhiSum              = reco.fCutEtaPhiSum;
198   fCutEtaPhiSeparate         = reco.fCutEtaPhiSeparate;
199   fCutR                      = reco.fCutR;
200   fCutEta                    = reco.fCutEta;
201   fCutPhi                    = reco.fCutPhi;
202   fClusterWindow             = reco.fClusterWindow;
203   fMass                      = reco.fMass;
204   fStepSurface               = reco.fStepSurface;
205   fStepCluster               = reco.fStepCluster;
206
207   fTrackCutsType             = reco.fTrackCutsType;
208   fCutMinTrackPt             = reco.fCutMinTrackPt;
209   fCutMinNClusterTPC         = reco.fCutMinNClusterTPC;
210   fCutMinNClusterITS         = reco.fCutMinNClusterITS; 
211   fCutMaxChi2PerClusterTPC   = reco.fCutMaxChi2PerClusterTPC;
212   fCutMaxChi2PerClusterITS   = reco.fCutMaxChi2PerClusterITS;
213   fCutRequireTPCRefit        = reco.fCutRequireTPCRefit;
214   fCutRequireITSRefit        = reco.fCutRequireITSRefit;
215   fCutAcceptKinkDaughters    = reco.fCutAcceptKinkDaughters;
216   fCutMaxDCAToVertexXY       = reco.fCutMaxDCAToVertexXY;
217   fCutMaxDCAToVertexZ        = reco.fCutMaxDCAToVertexZ;
218   fCutDCAToVertex2D          = reco.fCutDCAToVertex2D;
219   
220   if(reco.fResidualEta)
221   {
222     // assign or copy construct
223     if(fResidualEta)
224     { 
225       *fResidualEta = *reco.fResidualEta;
226     }
227     else 
228     {
229       fResidualEta = new TArrayF(*reco.fResidualEta);
230     }
231   }
232   else
233   {
234     if(fResidualEta)delete fResidualEta;
235     fResidualEta = 0;
236   }
237   
238   if(reco.fResidualPhi)
239   {
240     // assign or copy construct
241     if(fResidualPhi)
242     { 
243       *fResidualPhi = *reco.fResidualPhi;
244     }
245     else 
246     {
247       fResidualPhi = new TArrayF(*reco.fResidualPhi);
248     }
249   }
250   else
251   {
252     if(fResidualPhi)delete fResidualPhi;
253     fResidualPhi = 0;
254   }
255   
256   if(reco.fMatchedTrackIndex)
257   {
258     // assign or copy construct
259     if(fMatchedTrackIndex)
260     { 
261       *fMatchedTrackIndex = *reco.fMatchedTrackIndex;
262     }
263     else 
264     { 
265       fMatchedTrackIndex = new TArrayI(*reco.fMatchedTrackIndex);
266     }
267   }
268   else
269   {
270     if(fMatchedTrackIndex)delete fMatchedTrackIndex;
271     fMatchedTrackIndex = 0;
272   }  
273   
274   if(reco.fMatchedClusterIndex)
275   {
276     // assign or copy construct
277     if(fMatchedClusterIndex)
278     { 
279       *fMatchedClusterIndex = *reco.fMatchedClusterIndex;
280     }
281     else 
282     {
283       fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex);
284     }
285   }
286   else
287   {
288     if(fMatchedClusterIndex)delete fMatchedClusterIndex;
289     fMatchedClusterIndex = 0;
290   }
291    
292   return *this;
293 }
294
295
296 //_____________________________________
297 AliEMCALRecoUtils::~AliEMCALRecoUtils()
298 {
299   //Destructor.
300   
301   if(fEMCALRecalibrationFactors) 
302   { 
303     fEMCALRecalibrationFactors->Clear();
304     delete fEMCALRecalibrationFactors;
305   }  
306   
307   if(fEMCALTimeRecalibrationFactors) 
308   { 
309     fEMCALTimeRecalibrationFactors->Clear();
310     delete fEMCALTimeRecalibrationFactors;
311   }  
312   
313   if(fEMCALBadChannelMap) 
314   { 
315     fEMCALBadChannelMap->Clear();
316     delete fEMCALBadChannelMap;
317   }
318  
319   delete fMatchedTrackIndex   ; 
320   delete fMatchedClusterIndex ; 
321   delete fResidualEta         ; 
322   delete fResidualPhi         ; 
323   delete fPIDUtils            ;
324
325   InitTrackCuts();
326 }
327
328 //_______________________________________________________________________________
329 Bool_t AliEMCALRecoUtils::AcceptCalibrateCell(const Int_t absID, const Int_t bc,
330                                               Float_t  & amp,    Double_t & time, 
331                                               AliVCaloCells* cells) 
332 {
333   // Reject cell if criteria not passed and calibrate it
334   
335   AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
336   
337   if(absID < 0 || absID >= 24*48*geom->GetNumberOfSuperModules()) return kFALSE;
338   
339   Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
340   
341   if(!geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta)) 
342   {
343     // cell absID does not exist
344     amp=0; time = 1.e9;
345     return kFALSE; 
346   }
347   
348   geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);  
349
350   // Do not include bad channels found in analysis,
351   if( IsBadChannelsRemovalSwitchedOn() && GetEMCALChannelStatus(imod, ieta, iphi)) 
352   {
353     return kFALSE;
354   }
355   
356   //Recalibrate energy
357   amp  = cells->GetCellAmplitude(absID);
358   if(!fCellsRecalibrated && IsRecalibrationOn())
359     amp *= GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
360   
361   
362   // Recalibrate time
363   time = cells->GetCellTime(absID);
364   
365   RecalibrateCellTime(absID,bc,time);
366   
367   return kTRUE;
368 }
369
370 //_____________________________________________________________________________
371 Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(const AliEMCALGeometry* geom, 
372                                                   const AliVCluster* cluster, 
373                                                   AliVCaloCells* cells) 
374 {
375   // Given the list of AbsId of the cluster, get the maximum cell and 
376   // check if there are fNCellsFromBorder from the calorimeter border
377   
378   if(!cluster)
379   {
380     AliInfo("Cluster pointer null!");
381     return kFALSE;
382   }
383   
384   //If the distance to the border is 0 or negative just exit accept all clusters
385   if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
386   
387   Int_t absIdMax  = -1, iSM =-1, ieta = -1, iphi = -1;
388   Bool_t shared = kFALSE;
389   GetMaxEnergyCell(geom, cells, cluster, absIdMax,  iSM, ieta, iphi, shared);
390   
391   AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n", 
392                   absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
393   
394   if(absIdMax==-1) return kFALSE;
395   
396   //Check if the cell is close to the borders:
397   Bool_t okrow = kFALSE;
398   Bool_t okcol = kFALSE;
399   
400   if(iSM < 0 || iphi < 0 || ieta < 0 ) 
401   {
402     AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
403                   iSM,ieta,iphi));
404   }
405   
406   //Check rows/phi
407   if(iSM < 10)
408   {
409     if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE; 
410   }
411   else if (iSM >=10 && ( ( geom->GetEMCGeometry()->GetGeoName()).Contains("12SMV1"))) 
412   {
413     if(iphi >= fNCellsFromEMCALBorder && iphi < 8-fNCellsFromEMCALBorder) okrow =kTRUE; //1/3 sm case
414   }
415   else 
416   {
417     if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE; // half SM case
418   }
419   
420   //Check columns/eta
421   if(!fNoEMCALBorderAtEta0)
422   {
423     if(ieta  > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE; 
424   }
425   else
426   {
427     if(iSM%2==0)
428     {
429       if(ieta >= fNCellsFromEMCALBorder)     okcol = kTRUE;  
430     }
431     else 
432     {
433       if(ieta <  48-fNCellsFromEMCALBorder)  okcol = kTRUE;  
434     }
435   }//eta 0 not checked
436   
437   AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d:  column? %d, row? %d\nq",
438                   fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
439   
440   if (okcol && okrow) 
441   {
442     //printf("Accept\n");
443     return kTRUE;
444   }
445   else  
446   {
447     //printf("Reject\n");
448     AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
449     return kFALSE;
450   }
451   
452 }  
453
454
455 //_______________________________________________________________________________
456 Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(const AliEMCALGeometry* geom, 
457                                                     const UShort_t* cellList, 
458                                                     const Int_t nCells)
459 {
460   // Check that in the cluster cells, there is no bad channel of those stored 
461   // in fEMCALBadChannelMap or fPHOSBadChannelMap
462   
463   if(!fRemoveBadChannels)  return kFALSE;
464   if(!fEMCALBadChannelMap) return kFALSE;
465   
466   Int_t icol = -1;
467   Int_t irow = -1;
468   Int_t imod = -1;
469   for(Int_t iCell = 0; iCell<nCells; iCell++)
470   {
471     //Get the column and row
472     Int_t iTower = -1, iIphi = -1, iIeta = -1; 
473     geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta); 
474     if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
475     geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);      
476     if(GetEMCALChannelStatus(imod, icol, irow))
477     {
478       AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
479       return kTRUE;
480     }
481     
482   }// cell cluster loop
483   
484   return kFALSE;
485 }
486
487 //_____________________________________________________________________________________________
488 Bool_t AliEMCALRecoUtils::IsExoticCell(const Int_t absID, AliVCaloCells* cells, const Int_t bc)
489 {
490   // Look to cell neighbourhood and reject if it seems exotic
491   // Do before recalibrating the cells
492
493   if(!fRejectExoticCells) return kFALSE;
494   
495   AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
496   
497   Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
498   geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta); 
499   geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);  
500   
501   //Get close cells index, energy and time, not in corners
502
503   Int_t absID1 = -1;
504   Int_t absID2 = -1;
505   
506   if( iphi < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
507   if( iphi > 0 )                                absID2 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
508   
509   // In case of cell in eta = 0 border, depending on SM shift the cross cell index
510
511   Int_t absID3 = -1;
512   Int_t absID4 = -1;
513   
514   
515   if     ( ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2) )
516   {
517     absID3 = geom-> GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
518     absID4 = geom-> GetAbsCellIdFromCellIndexes(imod,   iphi, ieta-1); 
519   }
520   else if( ieta == 0 && imod%2 )
521   {
522     absID3 = geom-> GetAbsCellIdFromCellIndexes(imod,   iphi, ieta+1);
523     absID4 = geom-> GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1); 
524   }
525   else
526   {
527     if( ieta < AliEMCALGeoParams::fgkEMCALCols-1 ) 
528       absID3 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
529     if( ieta > 0 )                                 
530       absID4 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1); 
531   }
532
533   //printf("IMOD %d, AbsId %d, a %d, b %d, c %d e %d \n",imod,absID,absID1,absID2,absID3,absID4);
534
535   
536   Float_t  ecell  = 0, ecell1  = 0, ecell2  = 0, ecell3  = 0, ecell4  = 0;
537   Double_t tcell  = 0, tcell1  = 0, tcell2  = 0, tcell3  = 0, tcell4  = 0;
538   Bool_t   accept = 0, accept1 = 0, accept2 = 0, accept3 = 0, accept4 = 0;
539   
540   accept  = AcceptCalibrateCell(absID, bc, ecell ,tcell ,cells); 
541     
542   if(!accept) return kTRUE; // reject this cell
543   
544   if(ecell < fExoticCellMinAmplitude) return kFALSE; // do not reject low energy cells
545   
546   accept1 = AcceptCalibrateCell(absID1,bc, ecell1,tcell1,cells); 
547   accept2 = AcceptCalibrateCell(absID2,bc, ecell2,tcell2,cells); 
548   accept3 = AcceptCalibrateCell(absID3,bc, ecell3,tcell3,cells); 
549   accept4 = AcceptCalibrateCell(absID4,bc, ecell4,tcell4,cells); 
550   
551   /*
552     printf("Cell absID %d \n",absID);
553     printf("\t  accept1 %d, accept2 %d, accept3 %d, accept4 %d\n",
554            accept1,accept2,accept3,accept4);
555     printf("\t id %d: id1 %d, id2 %d, id3 %d, id4 %d\n",
556            absID,absID1,absID2,absID3,absID4);
557     printf("\t e %f: e1 %f, e2 %f, e3 %f, e4 %f\n",
558            ecell,ecell1,ecell2,ecell3,ecell4);
559     printf("\t t %f: t1 %f, t2 %f, t3 %f, t4 %f;\n dt1 %f, dt2 %f, dt3 %f, dt4 %f\n",
560            tcell*1.e9,tcell1*1.e9,tcell2*1.e9,tcell3*1.e9,tcell4*1.e9,
561            TMath::Abs(tcell-tcell1)*1.e9, TMath::Abs(tcell-tcell2)*1.e9, TMath::Abs(tcell-tcell3)*1.e9, TMath::Abs(tcell-tcell4)*1.e9);
562   */
563   
564   if(TMath::Abs(tcell-tcell1)*1.e9 > fExoticCellDiffTime) ecell1 = 0 ;
565   if(TMath::Abs(tcell-tcell2)*1.e9 > fExoticCellDiffTime) ecell2 = 0 ;
566   if(TMath::Abs(tcell-tcell3)*1.e9 > fExoticCellDiffTime) ecell3 = 0 ;
567   if(TMath::Abs(tcell-tcell4)*1.e9 > fExoticCellDiffTime) ecell4 = 0 ;
568
569   Float_t eCross = ecell1+ecell2+ecell3+ecell4;
570
571   //printf("\t eCell %f, eCross %f, 1-eCross/eCell %f\n",ecell,eCross,1-eCross/ecell);
572   
573   if(1-eCross/ecell > fExoticCellFraction) 
574   {
575     AliDebug(2,Form("AliEMCALRecoUtils::IsExoticCell() - EXOTIC CELL id %d, eCell %f, eCross %f, 1-eCross/eCell %f\n",
576                     absID,ecell,eCross,1-eCross/ecell));
577     return kTRUE;
578   }
579
580   return kFALSE;
581 }
582
583 //___________________________________________________________________
584 Bool_t AliEMCALRecoUtils::IsExoticCluster(const AliVCluster *cluster, 
585                                           AliVCaloCells *cells, 
586                                           const Int_t bc) 
587 {
588   // Check if the cluster highest energy tower is exotic
589   
590   if(!cluster)
591   {
592     AliInfo("Cluster pointer null!");
593     return kFALSE;
594   }
595   
596   if(!fRejectExoticCluster) return kFALSE;
597   
598   // Get highest energy tower
599   AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
600   Int_t iSupMod = -1, absId = -1, ieta = -1, iphi = -1;
601   Bool_t shared = kFALSE;
602   GetMaxEnergyCell(geom, cells, cluster, absId, iSupMod, ieta, iphi, shared);
603   
604   return IsExoticCell(absId,cells,bc);
605   
606 }
607
608 //_______________________________________________________________________
609 Float_t AliEMCALRecoUtils::SmearClusterEnergy(const AliVCluster* cluster) 
610 {
611   //In case of MC analysis, smear energy to match resolution/calibration in real data
612   
613   if(!cluster)
614   {
615     AliInfo("Cluster pointer null!");
616     return 0;
617   }
618   
619   Float_t energy    = cluster->E() ;
620   Float_t rdmEnergy = energy ;
621   if(fSmearClusterEnergy)
622   {
623     rdmEnergy = fRandom.Gaus(energy,fSmearClusterParam[0] * TMath::Sqrt(energy) +
624                                     fSmearClusterParam[1] * energy +
625                                     fSmearClusterParam[2] );
626     AliDebug(2, Form("Energy: original %f, smeared %f\n", energy, rdmEnergy));
627   }
628   
629   return rdmEnergy;
630 }
631
632 //____________________________________________________________________________
633 Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster)
634 {
635   // Correct cluster energy from non linearity functions
636   
637   if(!cluster)
638   {
639     AliInfo("Cluster pointer null!");
640     return 0;
641   }
642   
643   Float_t energy = cluster->E();
644
645   switch (fNonLinearityFunction) 
646   {
647       
648     case kPi0MC:
649     {
650       //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
651       //Double_t fNonLinearityParams[0] = 1.014;
652       //Double_t fNonLinearityParams[1] = -0.03329;
653       //Double_t fNonLinearityParams[2] = -0.3853;
654       //Double_t fNonLinearityParams[3] = 0.5423;
655       //Double_t fNonLinearityParams[4] = -0.4335;
656        energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
657                   ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
658                     exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
659       break;
660     }
661      
662     case kPi0GammaGamma:
663     {
664       //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
665       //Double_t fNonLinearityParams[0] = 1.04;
666       //Double_t fNonLinearityParams[1] = -0.1445;
667       //Double_t fNonLinearityParams[2] = 1.046;
668       energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
669       break;
670     }
671       
672     case kPi0GammaConversion:
673     {
674       //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
675       //fNonLinearityParams[0] = 0.139393/0.1349766;
676       //fNonLinearityParams[1] = 0.0566186;
677       //fNonLinearityParams[2] = 0.982133;
678       energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
679       
680       break;
681     }
682       
683     case kBeamTest:
684     {
685       //From beam test, Alexei's results, for different ZS thresholds
686       //                        th=30 MeV; th = 45 MeV; th = 75 MeV
687       //fNonLinearityParams[0] = 1.007;      1.003;      1.002 
688       //fNonLinearityParams[1] = 0.894;      0.719;      0.797 
689       //fNonLinearityParams[2] = 0.246;      0.334;      0.358 
690       //Rescale the param[0] with 1.03
691       energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
692       
693       break;
694     }
695       
696     case kBeamTestCorrected:
697     {
698       //From beam test, corrected for material between beam and EMCAL
699       //fNonLinearityParams[0] =  0.99078
700       //fNonLinearityParams[1] =  0.161499;
701       //fNonLinearityParams[2] =  0.655166; 
702       //fNonLinearityParams[3] =  0.134101;
703       //fNonLinearityParams[4] =  163.282;
704       //fNonLinearityParams[5] =  23.6904;
705       //fNonLinearityParams[6] =  0.978;
706         energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
707
708       break;
709     }
710       
711     case kNoCorrection:
712       AliDebug(2,"No correction on the energy\n");
713       break;
714       
715   }
716
717   return energy;
718 }
719
720 //__________________________________________________
721 void AliEMCALRecoUtils::InitNonLinearityParam()
722 {
723   //Initialising Non Linearity Parameters
724   
725   if(fNonLinearityFunction == kPi0MC) 
726   {
727     fNonLinearityParams[0] = 1.014;
728     fNonLinearityParams[1] = -0.03329;
729     fNonLinearityParams[2] = -0.3853;
730     fNonLinearityParams[3] = 0.5423;
731     fNonLinearityParams[4] = -0.4335;
732   }
733   
734   if(fNonLinearityFunction == kPi0GammaGamma) 
735   {
736     fNonLinearityParams[0] = 1.04;
737     fNonLinearityParams[1] = -0.1445;
738     fNonLinearityParams[2] = 1.046;
739   }  
740
741   if(fNonLinearityFunction == kPi0GammaConversion) 
742   {
743     fNonLinearityParams[0] = 0.139393;
744     fNonLinearityParams[1] = 0.0566186;
745     fNonLinearityParams[2] = 0.982133;
746   }  
747
748   if(fNonLinearityFunction == kBeamTest) 
749   {
750     if(fNonLinearThreshold == 30) 
751     {
752       fNonLinearityParams[0] = 1.007; 
753       fNonLinearityParams[1] = 0.894; 
754       fNonLinearityParams[2] = 0.246; 
755     }
756     if(fNonLinearThreshold == 45) 
757     {
758       fNonLinearityParams[0] = 1.003; 
759       fNonLinearityParams[1] = 0.719; 
760       fNonLinearityParams[2] = 0.334; 
761     }
762     if(fNonLinearThreshold == 75) 
763     {
764       fNonLinearityParams[0] = 1.002; 
765       fNonLinearityParams[1] = 0.797; 
766       fNonLinearityParams[2] = 0.358; 
767     }
768   }
769
770   if(fNonLinearityFunction == kBeamTestCorrected) 
771   {
772     fNonLinearityParams[0] =  0.99078;
773     fNonLinearityParams[1] =  0.161499;
774     fNonLinearityParams[2] =  0.655166; 
775     fNonLinearityParams[3] =  0.134101;
776     fNonLinearityParams[4] =  163.282;
777     fNonLinearityParams[5] =  23.6904;
778     fNonLinearityParams[6] =  0.978;
779   }
780 }
781
782 //_________________________________________________________
783 Float_t  AliEMCALRecoUtils::GetDepth(const Float_t energy, 
784                                      const Int_t iParticle, 
785                                      const Int_t iSM) const 
786 {
787   //Calculate shower depth for a given cluster energy and particle type
788
789   // parameters 
790   Float_t x0    = 1.31;
791   Float_t ecr   = 8;
792   Float_t depth = 0;
793   
794   switch ( iParticle )
795   {
796     case kPhoton:
797       depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
798       break;
799       
800     case kElectron:
801       depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
802       break;
803       
804     case kHadron:
805       // hadron 
806       // boxes anc. here
807       if(gGeoManager)
808       {
809         gGeoManager->cd("ALIC_1/XEN1_1");
810         TGeoNode        *geoXEn1    = gGeoManager->GetCurrentNode();
811         TGeoNodeMatrix  *geoSM      = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
812         if(geoSM)
813         {
814           TGeoVolume      *geoSMVol   = geoSM->GetVolume(); 
815           TGeoShape       *geoSMShape = geoSMVol->GetShape();
816           TGeoBBox        *geoBox     = dynamic_cast<TGeoBBox *>(geoSMShape);
817           if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
818           else AliFatal("Null GEANT box");
819         }
820         else AliFatal("NULL  GEANT node matrix");
821       }
822       else
823       {//electron
824         depth = x0 * (TMath::Log(energy*1000 / ecr)  - 0.5); //Multiply energy by 1000 to transform to MeV
825       }
826         
827       break;
828       
829     default://photon
830       depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
831   }  
832   
833   return depth;
834 }
835
836 //____________________________________________________________________
837 void AliEMCALRecoUtils::GetMaxEnergyCell(const AliEMCALGeometry *geom, 
838                                          AliVCaloCells* cells, 
839                                          const AliVCluster* clu, 
840                                          Int_t  & absId,  
841                                          Int_t  & iSupMod, 
842                                          Int_t  & ieta, 
843                                          Int_t  & iphi, 
844                                          Bool_t & shared)
845 {
846   //For a given CaloCluster gets the absId of the cell 
847   //with maximum energy deposit.
848   
849   Double_t eMax        = -1.;
850   Double_t eCell       = -1.;
851   Float_t  fraction    = 1.;
852   Float_t  recalFactor = 1.;
853   Int_t    cellAbsId   = -1 ;
854
855   Int_t iTower  = -1;
856   Int_t iIphi   = -1;
857   Int_t iIeta   = -1;
858   Int_t iSupMod0= -1;
859
860   if(!clu)
861   {
862     AliInfo("Cluster pointer null!");
863     absId=-1; iSupMod0=-1, ieta = -1; iphi = -1; shared = -1;
864     return;
865   }
866   
867   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) 
868   {
869     cellAbsId = clu->GetCellAbsId(iDig);
870     fraction  = clu->GetCellAmplitudeFraction(iDig);
871     //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
872     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
873     geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta); 
874     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
875     if     (iDig==0) 
876     {
877       iSupMod0=iSupMod;
878     }
879     else if(iSupMod0!=iSupMod) 
880     {
881       shared = kTRUE;
882       //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
883     }
884     if(!fCellsRecalibrated && IsRecalibrationOn()) 
885     {
886       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
887     }
888     eCell  = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
889     //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
890     if(eCell > eMax)  
891     { 
892       eMax  = eCell; 
893       absId = cellAbsId;
894       //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
895     }
896   }// cell loop
897   
898   //Get from the absid the supermodule, tower and eta/phi numbers
899   geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta); 
900   //Gives SuperModule and Tower numbers
901   geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
902                                          iIphi, iIeta,iphi,ieta); 
903   //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
904   //printf("Max end---\n");
905 }
906
907 //______________________________________
908 void AliEMCALRecoUtils::InitParameters()
909 {
910   // Initialize data members with default values
911   
912   fParticleType = kPhoton;
913   fPosAlgo      = kUnchanged;
914   fW0           = 4.5;
915   
916   fNonLinearityFunction = kNoCorrection;
917   fNonLinearThreshold   = 30;
918   
919   fExoticCellFraction     = 0.97;
920   fExoticCellDiffTime     = 1e6;
921   fExoticCellMinAmplitude = 0.5;
922   
923   fAODFilterMask = 32;
924   
925   fCutEtaPhiSum      = kTRUE;
926   fCutEtaPhiSeparate = kFALSE;
927   
928   fCutR   = 0.05; 
929   fCutEta = 0.025; 
930   fCutPhi = 0.05;
931   
932   fClusterWindow = 100;
933   fMass          = 0.139;
934   
935   fStepSurface   = 20.;                      
936   fStepCluster   = 5.;
937   fTrackCutsType = kLooseCut;
938   
939   fCutMinTrackPt     = 0;
940   fCutMinNClusterTPC = -1;
941   fCutMinNClusterITS = -1;
942   
943   fCutMaxChi2PerClusterTPC  = 1e10;
944   fCutMaxChi2PerClusterITS  = 1e10;
945   
946   fCutRequireTPCRefit     = kFALSE;
947   fCutRequireITSRefit     = kFALSE;
948   fCutAcceptKinkDaughters = kFALSE;
949   
950   fCutMaxDCAToVertexXY = 1e10;             
951   fCutMaxDCAToVertexZ  = 1e10;              
952   fCutDCAToVertex2D    = kFALSE;
953   
954   
955   //Misalignment matrices
956   for(Int_t i = 0; i < 15 ; i++) 
957   {
958     fMisalTransShift[i] = 0.; 
959     fMisalRotShift[i]   = 0.; 
960   }
961   
962   //Non linearity
963   for(Int_t i = 0; i < 7  ; i++) fNonLinearityParams[i] = 0.; 
964   
965   //For kBeamTestCorrected case, but default is no correction
966   fNonLinearityParams[0] =  0.99078;
967   fNonLinearityParams[1] =  0.161499;
968   fNonLinearityParams[2] =  0.655166; 
969   fNonLinearityParams[3] =  0.134101;
970   fNonLinearityParams[4] =  163.282;
971   fNonLinearityParams[5] =  23.6904;
972   fNonLinearityParams[6] =  0.978;
973   
974   //For kPi0GammaGamma case
975   //fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
976   //fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
977   //fNonLinearityParams[2] = 1.046;
978   
979   //Cluster energy smearing
980   fSmearClusterEnergy   = kFALSE;
981   fSmearClusterParam[0] = 0.07; // * sqrt E term
982   fSmearClusterParam[1] = 0.00; // * E term
983   fSmearClusterParam[2] = 0.00; // constant
984 }
985
986 //_____________________________________________________
987 void AliEMCALRecoUtils::InitEMCALRecalibrationFactors()
988 {
989   //Init EMCAL recalibration factors
990   AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
991   //In order to avoid rewriting the same histograms
992   Bool_t oldStatus = TH1::AddDirectoryStatus();
993   TH1::AddDirectory(kFALSE);
994   
995   fEMCALRecalibrationFactors = new TObjArray(12);
996   for (int i = 0; i < 12; i++) 
997     fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),
998                                              Form("EMCALRecalFactors_SM%d",i),  48, 0, 48, 24, 0, 24));
999   //Init the histograms with 1
1000   for (Int_t sm = 0; sm < 12; sm++) 
1001   {
1002     for (Int_t i = 0; i < 48; i++) 
1003     {
1004       for (Int_t j = 0; j < 24; j++) 
1005       {
1006         SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
1007       }
1008     }
1009   }
1010   
1011   fEMCALRecalibrationFactors->SetOwner(kTRUE);
1012   fEMCALRecalibrationFactors->Compress();
1013   
1014   //In order to avoid rewriting the same histograms
1015   TH1::AddDirectory(oldStatus);    
1016 }
1017
1018 //_________________________________________________________
1019 void AliEMCALRecoUtils::InitEMCALTimeRecalibrationFactors()
1020 {
1021   //Init EMCAL recalibration factors
1022   AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
1023   //In order to avoid rewriting the same histograms
1024   Bool_t oldStatus = TH1::AddDirectoryStatus();
1025   TH1::AddDirectory(kFALSE);
1026   
1027   fEMCALTimeRecalibrationFactors = new TObjArray(4);
1028   for (int i = 0; i < 4; i++) 
1029     fEMCALTimeRecalibrationFactors->Add(new TH1F(Form("hAllTimeAvBC%d",i),
1030                                                  Form("hAllTimeAvBC%d",i),  
1031                                                  48*24*12,0.,48*24*12)          );
1032   //Init the histograms with 1
1033   for (Int_t bc = 0; bc < 4; bc++) 
1034   {
1035     for (Int_t i = 0; i < 48*24*12; i++) 
1036       SetEMCALChannelTimeRecalibrationFactor(bc,i,0.);
1037   }
1038   
1039   fEMCALTimeRecalibrationFactors->SetOwner(kTRUE);
1040   fEMCALTimeRecalibrationFactors->Compress();
1041   
1042   //In order to avoid rewriting the same histograms
1043   TH1::AddDirectory(oldStatus);    
1044 }
1045
1046 //____________________________________________________
1047 void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()
1048 {
1049   //Init EMCAL bad channels map
1050   AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
1051   //In order to avoid rewriting the same histograms
1052   Bool_t oldStatus = TH1::AddDirectoryStatus();
1053   TH1::AddDirectory(kFALSE);
1054   
1055   fEMCALBadChannelMap = new TObjArray(12);
1056   //TH2F * hTemp = new  TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
1057   for (int i = 0; i < 12; i++) 
1058   {
1059     fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
1060   }
1061   
1062   fEMCALBadChannelMap->SetOwner(kTRUE);
1063   fEMCALBadChannelMap->Compress();
1064   
1065   //In order to avoid rewriting the same histograms
1066   TH1::AddDirectory(oldStatus);    
1067 }
1068
1069 //____________________________________________________________________________
1070 void AliEMCALRecoUtils::RecalibrateClusterEnergy(const AliEMCALGeometry* geom, 
1071                                                  AliVCluster * cluster, 
1072                                                  AliVCaloCells * cells, 
1073                                                  const Int_t bc)
1074 {
1075   // Recalibrate the cluster energy and Time, considering the recalibration map 
1076   // and the energy of the cells and time that compose the cluster.
1077   // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
1078   
1079   if(!cluster)
1080   {
1081     AliInfo("Cluster pointer null!");
1082     return;
1083   }  
1084   
1085   //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
1086   UShort_t * index    = cluster->GetCellsAbsId() ;
1087   Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
1088   Int_t ncells = cluster->GetNCells();
1089   
1090   //Initialize some used variables
1091   Float_t energy = 0;
1092   Int_t   absId  =-1;
1093   Int_t   icol   =-1, irow =-1, imod=1;
1094   Float_t factor = 1, frac = 0;
1095   Int_t   absIdMax = -1;
1096   Float_t emax     = 0;
1097   
1098   //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
1099   for(Int_t icell = 0; icell < ncells; icell++)
1100   {
1101     absId = index[icell];
1102     frac =  fraction[icell];
1103     if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
1104     
1105     if(!fCellsRecalibrated && IsRecalibrationOn()) 
1106     {
1107       // Energy  
1108       Int_t iTower = -1, iIphi = -1, iIeta = -1; 
1109       geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta); 
1110       if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
1111       geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);      
1112       factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
1113       
1114       AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
1115                       imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
1116       
1117     } 
1118     
1119     energy += cells->GetCellAmplitude(absId)*factor*frac;
1120     
1121     if(emax < cells->GetCellAmplitude(absId)*factor*frac)
1122     {
1123       emax     = cells->GetCellAmplitude(absId)*factor*frac;
1124       absIdMax = absId;
1125     }
1126   }
1127   
1128   AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f \n",cluster->E(),energy));
1129
1130   cluster->SetE(energy);
1131
1132   // Recalculate time of cluster   
1133   Double_t timeorg = cluster->GetTOF();
1134   if(!fCellsRecalibrated && IsTimeRecalibrationOn())
1135   {
1136     Double_t time = timeorg;
1137     RecalibrateCellTime(absIdMax,bc,time);
1138     cluster->SetTOF(time);
1139   } 
1140
1141   AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Time before %f, after %f \n",timeorg,cluster->GetTOF()));
1142
1143 }
1144
1145 //_____________________________________________________________
1146 void AliEMCALRecoUtils::RecalibrateCells(AliVCaloCells * cells,
1147                                          const Int_t bc)
1148 {
1149   // Recalibrate the cells time and energy, considering the recalibration map and the energy 
1150   // of the cells that compose the cluster.
1151   // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
1152
1153   if(!IsRecalibrationOn() && !IsTimeRecalibrationOn()) return;
1154   
1155   if(!cells)
1156   {
1157     AliInfo("Cells pointer null!");
1158     return;
1159   }  
1160   
1161   Short_t  absId  =-1;
1162   Bool_t   accept = kFALSE;
1163   Float_t  ecell  = 0;
1164   Double_t tcell  = 0;
1165   Double_t ecellin = 0;
1166   Double_t tcellin = 0;
1167   Short_t  mclabel = -1;
1168   Double_t efrac = 0;
1169   
1170   Int_t nEMcell  = cells->GetNumberOfCells() ;  
1171   for (Int_t iCell = 0; iCell < nEMcell; iCell++) 
1172   { 
1173     cells->GetCell( iCell, absId, ecellin, tcellin, mclabel, efrac );
1174     
1175     accept = AcceptCalibrateCell(absId, bc, ecell ,tcell ,cells); 
1176     if(!accept) 
1177     {
1178       ecell = 0;
1179       tcell = -1;
1180     }
1181     
1182     //Set new values
1183     cells->SetCell(iCell,absId,ecell, tcell, mclabel, efrac);
1184   }
1185
1186   fCellsRecalibrated = kTRUE;
1187 }
1188
1189 //_______________________________________________________________________________________________________
1190 void AliEMCALRecoUtils::RecalibrateCellTime(const Int_t absId, const Int_t bc, Double_t & celltime) const
1191 {
1192   // Recalibrate time of cell with absID  considering the recalibration map 
1193   // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
1194   
1195   if(!fCellsRecalibrated && IsTimeRecalibrationOn() && bc >= 0)
1196   {
1197     celltime -= GetEMCALChannelTimeRecalibrationFactor(bc%4,absId)*1.e-9;    ;  
1198   }
1199 }
1200   
1201 //______________________________________________________________________________
1202 void AliEMCALRecoUtils::RecalculateClusterPosition(const AliEMCALGeometry *geom, 
1203                                                    AliVCaloCells* cells, 
1204                                                    AliVCluster* clu)
1205 {
1206   //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
1207   
1208   if(!clu)
1209   {
1210     AliInfo("Cluster pointer null!");
1211     return;
1212   }
1213     
1214   if     (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
1215   else if(fPosAlgo==kPosTowerIndex)  RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
1216   else   AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
1217 }  
1218
1219 //_____________________________________________________________________________________________
1220 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(const AliEMCALGeometry *geom, 
1221                                                                   AliVCaloCells* cells, 
1222                                                                   AliVCluster* clu)
1223 {
1224   // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
1225   // The algorithm is a copy of what is done in AliEMCALRecPoint
1226   
1227   Double_t eCell       = 0.;
1228   Float_t  fraction    = 1.;
1229   Float_t  recalFactor = 1.;
1230   
1231   Int_t    absId   = -1;
1232   Int_t    iTower  = -1, iIphi  = -1, iIeta  = -1;
1233   Int_t    iSupModMax = -1, iSM=-1, iphi   = -1, ieta   = -1;
1234   Float_t  weight = 0.,  totalWeight=0.;
1235   Float_t  newPos[3] = {0,0,0};
1236   Double_t pLocal[3], pGlobal[3];
1237   Bool_t shared = kFALSE;
1238
1239   Float_t  clEnergy = clu->E(); //Energy already recalibrated previously
1240   GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi,shared);
1241   Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
1242   
1243   //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
1244   
1245   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) 
1246   {
1247     absId = clu->GetCellAbsId(iDig);
1248     fraction  = clu->GetCellAmplitudeFraction(iDig);
1249     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1250     
1251     if (!fCellsRecalibrated)
1252     {
1253       geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta); 
1254       geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);      
1255       
1256       if(IsRecalibrationOn()) 
1257       {
1258         recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
1259       }
1260     }
1261     
1262     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
1263     
1264     weight = GetCellWeight(eCell,clEnergy);
1265     totalWeight += weight;
1266     
1267     geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
1268     //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
1269     geom->GetGlobal(pLocal,pGlobal,iSupModMax);
1270     //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
1271
1272     for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
1273   }// cell loop
1274   
1275   if(totalWeight>0)
1276   {
1277     for(int i=0; i<3; i++ )    newPos[i] /= totalWeight;
1278   }
1279     
1280   //Float_t pos[]={0,0,0};
1281   //clu->GetPosition(pos);
1282   //printf("OldPos  : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
1283   //printf("NewPos  : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
1284   
1285   if(iSupModMax > 1) //sector 1
1286   {
1287     newPos[0] +=fMisalTransShift[3];//-=3.093; 
1288     newPos[1] +=fMisalTransShift[4];//+=6.82;
1289     newPos[2] +=fMisalTransShift[5];//+=1.635;
1290     //printf("   +    : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
1291   } else //sector 0
1292   {
1293     newPos[0] +=fMisalTransShift[0];//+=1.134;
1294     newPos[1] +=fMisalTransShift[1];//+=8.2;
1295     newPos[2] +=fMisalTransShift[2];//+=1.197;
1296     //printf("   +    : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
1297   }
1298   //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
1299
1300   clu->SetPosition(newPos);
1301 }  
1302
1303 //____________________________________________________________________________________________
1304 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(const AliEMCALGeometry *geom, 
1305                                                                  AliVCaloCells* cells, 
1306                                                                  AliVCluster* clu)
1307 {
1308   // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
1309   // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
1310   
1311   Double_t eCell       = 1.;
1312   Float_t  fraction    = 1.;
1313   Float_t  recalFactor = 1.;
1314   
1315   Int_t absId   = -1;
1316   Int_t iTower  = -1;
1317   Int_t iIphi   = -1, iIeta   = -1;
1318   Int_t iSupMod = -1, iSupModMax = -1;
1319   Int_t iphi = -1, ieta =-1;
1320   Bool_t shared = kFALSE;
1321
1322   Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
1323   GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi,shared);
1324   Float_t  depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
1325
1326   Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
1327   Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
1328   Int_t startingSM = -1;
1329   
1330   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) 
1331   {
1332     absId = clu->GetCellAbsId(iDig);
1333     fraction  = clu->GetCellAmplitudeFraction(iDig);
1334     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1335
1336     if     (iDig==0)  startingSM = iSupMod;
1337     else if(iSupMod != startingSM) areInSameSM = kFALSE;
1338
1339     eCell  = cells->GetCellAmplitude(absId);
1340     
1341     geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta); 
1342     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);    
1343     
1344     if (!fCellsRecalibrated)
1345     {
1346       if(IsRecalibrationOn()) 
1347       {
1348         recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1349       }
1350     }
1351     
1352     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
1353     
1354     weight = GetCellWeight(eCell,clEnergy);
1355     if(weight < 0) weight = 0;
1356     totalWeight += weight;
1357     weightedCol += ieta*weight;
1358     weightedRow += iphi*weight;
1359     
1360     //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
1361   }// cell loop
1362     
1363   Float_t xyzNew[]={0.,0.,0.};
1364   if(areInSameSM == kTRUE) 
1365   {
1366     //printf("In Same SM\n");
1367     weightedCol = weightedCol/totalWeight;
1368     weightedRow = weightedRow/totalWeight;
1369     geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew); 
1370   } 
1371   else 
1372   {
1373     //printf("In Different SM\n");
1374     geom->RecalculateTowerPosition(iphi,        ieta,        iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew); 
1375   }
1376   
1377   clu->SetPosition(xyzNew);
1378 }
1379
1380 //___________________________________________________________________________________________
1381 void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(const AliEMCALGeometry * geom, 
1382                                                                AliVCaloCells* cells, 
1383                                                                AliVCluster * cluster)
1384 {           
1385   //re-evaluate distance to bad channel with updated bad map
1386   
1387   if(!fRecalDistToBadChannels) return;
1388   
1389   if(!cluster)
1390   {
1391     AliInfo("Cluster pointer null!");
1392     return;
1393   }  
1394   
1395   //Get channels map of the supermodule where the cluster is.
1396   Int_t absIdMax  = -1, iSupMod =-1, icolM = -1, irowM = -1;
1397   Bool_t shared = kFALSE;
1398   GetMaxEnergyCell(geom, cells, cluster, absIdMax,  iSupMod, icolM, irowM, shared);
1399   TH2D* hMap  = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
1400
1401   Int_t dRrow, dRcol;  
1402   Float_t  minDist = 10000.;
1403   Float_t  dist    = 0.;
1404   
1405   //Loop on tower status map 
1406   for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
1407   {
1408     for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
1409     {
1410       //Check if tower is bad.
1411       if(hMap->GetBinContent(icol,irow)==0) continue;
1412       //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
1413       //       iSupMod,icol, irow, icolM,irowM);
1414       
1415       dRrow=TMath::Abs(irowM-irow);
1416       dRcol=TMath::Abs(icolM-icol);
1417       dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
1418       if(dist < minDist)
1419       {
1420         //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
1421         minDist = dist;
1422       }
1423     }
1424   }
1425   
1426   //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
1427   if (shared) 
1428   {
1429     TH2D* hMap2 = 0;
1430     Int_t iSupMod2 = -1;
1431     
1432     //The only possible combinations are (0,1), (2,3) ... (8,9)
1433     if(iSupMod%2) iSupMod2 = iSupMod-1;
1434     else          iSupMod2 = iSupMod+1;
1435     hMap2  = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
1436     
1437     //Loop on tower status map of second super module
1438     for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
1439     {
1440       for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
1441       {
1442         //Check if tower is bad.
1443         if(hMap2->GetBinContent(icol,irow)==0) continue;
1444         //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",
1445         //     iSupMod2,icol, irow,iSupMod,icolM,irowM);
1446         dRrow=TMath::Abs(irow-irowM);
1447         
1448         if(iSupMod%2) 
1449         {
1450           dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
1451         } else 
1452         {
1453           dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
1454         }                    
1455         
1456         dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
1457         if(dist < minDist) minDist = dist;        
1458       }
1459     }
1460   }// shared cluster in 2 SuperModules
1461   
1462   AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
1463   cluster->SetDistanceToBadChannel(minDist);
1464 }
1465
1466 //__________________________________________________________________
1467 void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster)
1468 {           
1469   //re-evaluate identification parameters with bayesian
1470   
1471   if(!cluster)
1472   {
1473     AliInfo("Cluster pointer null!");
1474     return;
1475   }
1476   
1477   if ( cluster->GetM02() != 0)
1478     fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
1479   
1480   Float_t pidlist[AliPID::kSPECIESN+1];
1481   for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
1482         
1483   cluster->SetPID(pidlist);
1484 }
1485
1486 //___________________________________________________________________________________________________________________
1487 void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(const AliEMCALGeometry * geom, 
1488                                                                 AliVCaloCells* cells, 
1489                                                                 AliVCluster * cluster,
1490                                                                 Float_t & l0,   Float_t & l1,   
1491                                                                 Float_t & disp, Float_t & dEta, Float_t & dPhi,
1492                                                                 Float_t & sEta, Float_t & sPhi, Float_t & sEtaPhi)
1493 {
1494   // Calculates new center of gravity in the local EMCAL-module coordinates 
1495   // and tranfers into global ALICE coordinates
1496   // Calculates Dispersion and main axis
1497   
1498   if(!cluster)
1499   {
1500     AliInfo("Cluster pointer null!");
1501     return;
1502   }
1503     
1504   Double_t eCell       = 0.;
1505   Float_t  fraction    = 1.;
1506   Float_t  recalFactor = 1.;
1507
1508   Int_t    iSupMod = -1;
1509   Int_t    iTower  = -1;
1510   Int_t    iIphi   = -1;
1511   Int_t    iIeta   = -1;
1512   Int_t    iphi    = -1;
1513   Int_t    ieta    = -1;
1514   Double_t etai    = -1.;
1515   Double_t phii    = -1.;
1516   
1517   Int_t    nstat   = 0 ;
1518   Float_t  wtot    = 0.;
1519   Double_t w       = 0.;
1520   Double_t etaMean = 0.;
1521   Double_t phiMean = 0.;
1522     
1523   //Loop on cells
1524   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) 
1525   {
1526     //Get from the absid the supermodule, tower and eta/phi numbers
1527     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta); 
1528     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);        
1529     
1530     //Get the cell energy, if recalibration is on, apply factors
1531     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
1532     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1533     
1534     if (!fCellsRecalibrated)
1535     {
1536       if(IsRecalibrationOn()) 
1537       {
1538         recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1539       }
1540     }
1541     
1542     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1543     
1544     if(cluster->E() > 0 && eCell > 0)
1545     {
1546       w  = GetCellWeight(eCell,cluster->E());
1547       
1548       etai=(Double_t)ieta;
1549       phii=(Double_t)iphi;  
1550       
1551       if(w > 0.0) 
1552       {
1553         wtot += w ;
1554         nstat++;            
1555         //Shower shape
1556         sEta     += w * etai * etai ;
1557         etaMean  += w * etai ;
1558         sPhi     += w * phii * phii ;
1559         phiMean  += w * phii ; 
1560         sEtaPhi  += w * etai * phii ; 
1561       }
1562     } 
1563     else
1564       AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1565   }//cell loop
1566   
1567   //Normalize to the weight  
1568   if (wtot > 0) 
1569   {
1570     etaMean /= wtot ;
1571     phiMean /= wtot ;
1572   }
1573   else
1574     AliError(Form("Wrong weight %f\n", wtot));
1575   
1576   //Calculate dispersion  
1577   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) 
1578   {
1579     //Get from the absid the supermodule, tower and eta/phi numbers
1580     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta); 
1581     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1582     
1583     //Get the cell energy, if recalibration is on, apply factors
1584     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
1585     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1586     if (IsRecalibrationOn()) 
1587     {
1588       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1589     }
1590     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1591     
1592     if(cluster->E() > 0 && eCell > 0)
1593     {
1594       w  = GetCellWeight(eCell,cluster->E());
1595       
1596       etai=(Double_t)ieta;
1597       phii=(Double_t)iphi;    
1598       if(w > 0.0) 
1599       { 
1600         disp +=  w *((etai-etaMean)*(etai-etaMean)+(phii-phiMean)*(phii-phiMean)); 
1601         dEta +=  w * (etai-etaMean)*(etai-etaMean) ; 
1602         dPhi +=  w * (phii-phiMean)*(phii-phiMean) ; 
1603       }
1604     }
1605     else
1606       AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1607   }// cell loop
1608   
1609   //Normalize to the weigth and set shower shape parameters
1610   if (wtot > 0 && nstat > 1) 
1611   {
1612     disp    /= wtot ;
1613     dEta    /= wtot ;
1614     dPhi    /= wtot ;
1615     sEta    /= wtot ;
1616     sPhi    /= wtot ;
1617     sEtaPhi /= wtot ;
1618     
1619     sEta    -= etaMean * etaMean ;
1620     sPhi    -= phiMean * phiMean ;
1621     sEtaPhi -= etaMean * phiMean ;
1622     
1623     l0 = (0.5 * (sEta + sPhi) + TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
1624     l1 = (0.5 * (sEta + sPhi) - TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
1625   }
1626   else
1627   {
1628     l0   = 0. ;
1629     l1   = 0. ;
1630     dEta = 0. ; dPhi = 0. ; disp    = 0. ;
1631     sEta = 0. ; sPhi = 0. ; sEtaPhi = 0. ;
1632   }  
1633   
1634 }
1635
1636 //____________________________________________________________________________________________
1637 void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(const AliEMCALGeometry * geom, 
1638                                                                 AliVCaloCells* cells, 
1639                                                                 AliVCluster * cluster)
1640 {
1641   // Calculates new center of gravity in the local EMCAL-module coordinates 
1642   // and tranfers into global ALICE coordinates
1643   // Calculates Dispersion and main axis and puts them into the cluster
1644   
1645   Float_t l0   = 0., l1   = 0.;
1646   Float_t disp = 0., dEta = 0., dPhi    = 0.; 
1647   Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
1648   
1649   AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(geom,cells,cluster,l0,l1,disp,
1650                                                              dEta, dPhi, sEta, sPhi, sEtaPhi);
1651   
1652   cluster->SetM02(l0);
1653   cluster->SetM20(l1);
1654   if(disp > 0. ) cluster->SetDispersion(TMath::Sqrt(disp)) ;
1655   
1656
1657
1658 //____________________________________________________________________________
1659 void AliEMCALRecoUtils::FindMatches(AliVEvent *event,
1660                                     TObjArray * clusterArr,  
1661                                     const AliEMCALGeometry *geom)
1662 {
1663   //This function should be called before the cluster loop
1664   //Before call this function, please recalculate the cluster positions
1665   //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
1666   //Store matched cluster indexes and residuals
1667   
1668   fMatchedTrackIndex  ->Reset();
1669   fMatchedClusterIndex->Reset();
1670   fResidualPhi->Reset();
1671   fResidualEta->Reset();
1672   
1673   fMatchedTrackIndex  ->Set(1000);
1674   fMatchedClusterIndex->Set(1000);
1675   fResidualPhi->Set(1000);
1676   fResidualEta->Set(1000);
1677   
1678   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1679   AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
1680   
1681   // init the magnetic field if not already on
1682   if(!TGeoGlobalMagField::Instance()->GetField())
1683   {
1684     AliInfo("Init the magnetic field\n");
1685     if     (esdevent) 
1686     {
1687       esdevent->InitMagneticField();
1688     }
1689     else if(aodevent)
1690     {
1691       Double_t curSol = 30000*aodevent->GetMagneticField()/5.00668;
1692       Double_t curDip = 6000 *aodevent->GetMuonMagFieldScale();
1693       AliMagF *field  = AliMagF::CreateFieldMap(curSol,curDip);
1694       TGeoGlobalMagField::Instance()->SetField(field);
1695     }
1696     else
1697     {
1698       AliInfo("Mag Field not initialized, null esd/aod evetn pointers");
1699     }
1700     
1701   } // Init mag field
1702   
1703   TObjArray *clusterArray = 0x0;
1704   if(!clusterArr)
1705     {
1706       clusterArray = new TObjArray(event->GetNumberOfCaloClusters());
1707       for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1708   {
1709     AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1710     if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1711     clusterArray->AddAt(cluster,icl);
1712   }
1713     }
1714   
1715   Int_t    matched=0;
1716   Double_t cv[21];
1717   for (Int_t i=0; i<21;i++) cv[i]=0;
1718   for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
1719   {
1720     AliExternalTrackParam *trackParam = 0;
1721     
1722     //If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner 
1723     AliESDtrack *esdTrack = 0;
1724     AliAODTrack *aodTrack = 0;
1725     if(esdevent)
1726     {
1727       esdTrack = esdevent->GetTrack(itr);
1728       if(!esdTrack) continue;
1729       if(!IsAccepted(esdTrack)) continue;
1730       if(esdTrack->Pt()<fCutMinTrackPt) continue;
1731       Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
1732       if(TMath::Abs(esdTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue;
1733       trackParam =  const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam());
1734     }
1735     
1736     //If the input event is AOD, the starting point for extrapolation is at vertex
1737     //AOD tracks are selected according to its filterbit.
1738     else if(aodevent)
1739     {
1740       aodTrack = aodevent->GetTrack(itr);
1741       if(!aodTrack) continue;
1742       if(!aodTrack->TestFilterMask(fAODFilterMask)) continue; //Select AOD tracks that fulfill GetStandardITSTPCTrackCuts2010()
1743       if(aodTrack->Pt()<fCutMinTrackPt) continue;
1744       Double_t phi = aodTrack->Phi()*TMath::RadToDeg();
1745       if(TMath::Abs(aodTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue;
1746       Double_t pos[3],mom[3];
1747       aodTrack->GetXYZ(pos);
1748       aodTrack->GetPxPyPz(mom);
1749       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()));
1750       trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
1751     }
1752     
1753     //Return if the input data is not "AOD" or "ESD"
1754     else
1755     {
1756       printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n");
1757       if(clusterArray)
1758   {
1759     clusterArray->Clear();
1760     delete clusterArray;
1761   }
1762       return;
1763     }
1764     
1765     if(!trackParam) continue;
1766
1767     //Extrapolate the track to EMCal surface
1768     AliExternalTrackParam emcalParam(*trackParam);
1769     Float_t eta, phi;
1770     if(!ExtrapolateTrackToEMCalSurface(&emcalParam, 430., fMass, fStepSurface, eta, phi)) 
1771       {
1772   if(aodevent && trackParam) delete trackParam;
1773   continue;
1774       }
1775
1776 //    if(esdevent)
1777 //      {
1778 //  esdTrack->SetOuterParam(&emcalParam,AliExternalTrackParam::kMultSec);
1779 //      }
1780
1781     if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad())
1782       {
1783   if(aodevent && trackParam) delete trackParam;
1784   continue;
1785       }
1786
1787
1788     //Find matched clusters
1789     Int_t index = -1;
1790     Float_t dEta = -999, dPhi = -999;
1791     if(!clusterArr)
1792       {
1793   index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArray, dEta, dPhi);  
1794       }
1795     else
1796       {
1797   index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);  
1798       }  
1799     
1800     if(index>-1)
1801     {
1802       fMatchedTrackIndex   ->AddAt(itr,matched);
1803       fMatchedClusterIndex ->AddAt(index,matched);
1804       fResidualEta         ->AddAt(dEta,matched);
1805       fResidualPhi         ->AddAt(dPhi,matched);
1806       matched++;
1807     }
1808     if(aodevent && trackParam) delete trackParam;
1809   }//track loop
1810
1811   if(clusterArray)
1812     {
1813       clusterArray->Clear();
1814       delete clusterArray;
1815     }
1816   
1817   AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1818   
1819   fMatchedTrackIndex   ->Set(matched);
1820   fMatchedClusterIndex ->Set(matched);
1821   fResidualPhi         ->Set(matched);
1822   fResidualEta         ->Set(matched);
1823 }
1824
1825 //________________________________________________________________________________
1826 Int_t AliEMCALRecoUtils::FindMatchedClusterInEvent(const AliESDtrack *track, 
1827                                                    const AliVEvent *event, 
1828                                                    const AliEMCALGeometry *geom, 
1829                                                    Float_t &dEta, Float_t &dPhi)
1830 {
1831   //
1832   // This function returns the index of matched cluster to input track
1833   // Returns -1 if no match is found
1834   Int_t index = -1;
1835   Double_t phiV = track->Phi()*TMath::RadToDeg();
1836   if(TMath::Abs(track->Eta())>0.8 || phiV <= 20 || phiV >= 240 ) return index;
1837   AliExternalTrackParam *trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
1838   if(!trackParam) return index;
1839   AliExternalTrackParam emcalParam(*trackParam);
1840   Float_t eta, phi;
1841   if(!ExtrapolateTrackToEMCalSurface(&emcalParam, 430., fMass, fStepSurface, eta, phi)) return index;
1842   if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()) return index;
1843
1844   TObjArray *clusterArr = new TObjArray(event->GetNumberOfCaloClusters());
1845
1846   for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1847   {
1848     AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1849     if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1850     clusterArr->AddAt(cluster,icl);
1851   }
1852
1853   index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);  
1854   clusterArr->Clear();
1855   delete clusterArr;
1856   
1857   return index;
1858 }
1859
1860 //_______________________________________________________________________________________________
1861 Int_t  AliEMCALRecoUtils::FindMatchedClusterInClusterArr(const AliExternalTrackParam *emcalParam, 
1862                                                          AliExternalTrackParam *trkParam, 
1863                                                          const TObjArray * clusterArr, 
1864                                                          Float_t &dEta, Float_t &dPhi)
1865 {
1866   // Find matched cluster in array
1867   
1868   dEta=-999, dPhi=-999;
1869   Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
1870   Int_t index = -1;
1871   Float_t tmpEta=-999, tmpPhi=-999;
1872
1873   Double_t exPos[3] = {0.,0.,0.};
1874   if(!emcalParam->GetXYZ(exPos)) return index;
1875
1876   Float_t clsPos[3] = {0.,0.,0.};
1877   for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
1878     {
1879       AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
1880       if(!cluster || !cluster->IsEMCAL()) continue;
1881       cluster->GetPosition(clsPos);
1882       Double_t dR = TMath::Sqrt(TMath::Power(exPos[0]-clsPos[0],2)+TMath::Power(exPos[1]-clsPos[1],2)+TMath::Power(exPos[2]-clsPos[2],2));
1883       if(dR > fClusterWindow) continue;
1884
1885       AliExternalTrackParam trkPamTmp (*trkParam);//Retrieve the starting point every time before the extrapolation
1886       if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, fMass, fStepCluster, tmpEta, tmpPhi)) continue;
1887       if(fCutEtaPhiSum)
1888         {
1889           Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1890           if(tmpR<dRMax)
1891       {
1892         dRMax=tmpR;
1893         dEtaMax=tmpEta;
1894         dPhiMax=tmpPhi;
1895         index=icl;
1896       }
1897         }
1898       else if(fCutEtaPhiSeparate)
1899         {
1900           if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
1901       {
1902         dEtaMax = tmpEta;
1903         dPhiMax = tmpPhi;
1904         index=icl;
1905       }
1906         }
1907       else
1908         {
1909           printf("Error: please specify your cut criteria\n");
1910           printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1911           printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1912           return index;
1913         }
1914     }
1915
1916   dEta=dEtaMax;
1917   dPhi=dPhiMax;
1918
1919   return index;
1920 }
1921
1922 //------------------------------------------------------------------------------------
1923 Bool_t AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(AliExternalTrackParam *trkParam, 
1924                                                          const Double_t emcalR,
1925                                                          const Double_t mass, 
1926                                                          const Double_t step, 
1927                                                          Float_t &eta, 
1928                                                          Float_t &phi)
1929 {
1930   //Extrapolate track to EMCAL surface
1931   
1932   eta = -999, phi = -999;
1933   if(!trkParam) return kFALSE;
1934   if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, emcalR, mass, step, kTRUE, 0.8, -1)) return kFALSE;
1935   Double_t trkPos[3] = {0.,0.,0.};
1936   if(!trkParam->GetXYZ(trkPos)) return kFALSE;
1937   TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
1938   eta = trkPosVec.Eta();
1939   phi = trkPosVec.Phi();
1940   if(phi<0)
1941     phi += 2*TMath::Pi();
1942
1943   return kTRUE;
1944 }
1945
1946 //-----------------------------------------------------------------------------------
1947 Bool_t AliEMCALRecoUtils::ExtrapolateTrackToPosition(AliExternalTrackParam *trkParam, 
1948                                                      const Float_t *clsPos, 
1949                                                      Double_t mass, 
1950                                                      Double_t step, 
1951                                                      Float_t &tmpEta, 
1952                                                      Float_t &tmpPhi)
1953 {
1954   //
1955   //Return the residual by extrapolating a track param to a global position
1956   //
1957   tmpEta = -999;
1958   tmpPhi = -999;
1959   if(!trkParam) return kFALSE;
1960   Double_t trkPos[3] = {0.,0.,0.};
1961   TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
1962   Double_t alpha =  ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
1963   vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
1964   if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), mass, step,kTRUE, 0.8, -1)) return kFALSE;
1965   if(!trkParam->GetXYZ(trkPos)) return kFALSE; //Get the extrapolated global position
1966
1967   TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1968   TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
1969
1970   // track cluster matching
1971   tmpPhi = clsPosVec.DeltaPhi(trkPosVec);    // tmpPhi is between -pi and pi
1972   tmpEta = clsPosVec.Eta()-trkPosVec.Eta();
1973
1974   return kTRUE;
1975 }
1976
1977 //----------------------------------------------------------------------------------
1978 Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, 
1979                                                     const AliVCluster *cluster, 
1980                                                     const Double_t mass, 
1981                                                     const Double_t step, 
1982                                                     Float_t &tmpEta, 
1983                                                     Float_t &tmpPhi)
1984 {
1985   //
1986   //Return the residual by extrapolating a track param to a cluster
1987   //
1988   tmpEta = -999;
1989   tmpPhi = -999;
1990   if(!cluster || !trkParam) return kFALSE;
1991
1992   Float_t clsPos[3] = {0.,0.,0.};
1993   cluster->GetPosition(clsPos);
1994
1995   return ExtrapolateTrackToPosition(trkParam, clsPos, mass, step, tmpEta, tmpPhi);
1996 }
1997
1998 //---------------------------------------------------------------------------------
1999 Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, 
2000                                                     const AliVCluster *cluster, 
2001                                                     Float_t &tmpEta, 
2002                                                     Float_t &tmpPhi)
2003 {
2004   //
2005   //Return the residual by extrapolating a track param to a clusterfStepCluster
2006   //
2007
2008   return ExtrapolateTrackToCluster(trkParam, cluster, fMass, fStepCluster, tmpEta, tmpPhi);
2009 }
2010
2011 //_______________________________________________________________________
2012 void AliEMCALRecoUtils::GetMatchedResiduals(const Int_t clsIndex, 
2013                                             Float_t &dEta, Float_t &dPhi)
2014 {
2015   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
2016   //Get the residuals dEta and dPhi for this cluster to the closest track
2017   //Works with ESDs and AODs
2018
2019   if( FindMatchedPosForCluster(clsIndex) >= 999 )
2020   {
2021     AliDebug(2,"No matched tracks found!\n");
2022     dEta=999.;
2023     dPhi=999.;
2024     return;
2025   }
2026   dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
2027   dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
2028 }
2029
2030 //______________________________________________________________________________________________
2031 void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi)
2032 {
2033   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
2034   //Get the residuals dEta and dPhi for this track to the closest cluster
2035   //Works with ESDs and AODs
2036
2037   if( FindMatchedPosForTrack(trkIndex) >= 999 )
2038   {
2039     AliDebug(2,"No matched cluster found!\n");
2040     dEta=999.;
2041     dPhi=999.;
2042     return;
2043   }
2044   dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
2045   dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
2046 }
2047
2048 //__________________________________________________________
2049 Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
2050 {
2051   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
2052   //Get the index of matched track to this cluster
2053   //Works with ESDs and AODs
2054   
2055   if(IsClusterMatched(clsIndex))
2056     return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
2057   else 
2058     return -1; 
2059 }
2060
2061 //__________________________________________________________
2062 Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
2063 {
2064   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
2065   //Get the index of matched cluster to this track
2066   //Works with ESDs and AODs
2067   
2068   if(IsTrackMatched(trkIndex))
2069     return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
2070   else 
2071     return -1; 
2072 }
2073
2074 //______________________________________________________________
2075 Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex) const
2076 {
2077   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
2078   //Returns if the cluster has a match
2079   if(FindMatchedPosForCluster(clsIndex) < 999) 
2080     return kTRUE;
2081   else
2082     return kFALSE;
2083 }
2084
2085 //____________________________________________________________
2086 Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex) const 
2087 {
2088   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
2089   //Returns if the track has a match
2090   if(FindMatchedPosForTrack(trkIndex) < 999) 
2091     return kTRUE;
2092   else
2093     return kFALSE;
2094 }
2095
2096 //______________________________________________________________________
2097 UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
2098 {
2099   //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
2100   //Returns the position of the match in the fMatchedClusterIndex array
2101   Float_t tmpR = fCutR;
2102   UInt_t pos = 999;
2103   
2104   for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++) 
2105   {
2106     if(fMatchedClusterIndex->At(i)==clsIndex) 
2107     {
2108       Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
2109       if(r<tmpR) 
2110       {
2111         pos=i;
2112         tmpR=r;
2113         AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
2114                         fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
2115       }
2116     }
2117   }
2118   return pos;
2119 }
2120
2121 //____________________________________________________________________
2122 UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
2123 {
2124   //Given a track index as in AliESDEvent::GetTrack(trkIndex)
2125   //Returns the position of the match in the fMatchedTrackIndex array
2126   Float_t tmpR = fCutR;
2127   UInt_t pos = 999;
2128   
2129   for(Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++) 
2130   {
2131     if(fMatchedTrackIndex->At(i)==trkIndex) 
2132     {
2133       Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
2134       if(r<tmpR) 
2135       {
2136         pos=i;
2137         tmpR=r;
2138         AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
2139                         fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
2140       }
2141     }
2142   }
2143   return pos;
2144 }
2145
2146 //__________________________________________________________________________
2147 Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster, 
2148                                         const AliEMCALGeometry *geom, 
2149                                         AliVCaloCells* cells,const Int_t bc)
2150 {
2151   // check if the cluster survives some quality cut
2152   //
2153   //
2154   Bool_t isGood=kTRUE;
2155
2156   if(!cluster || !cluster->IsEMCAL())              return kFALSE;
2157   
2158   if(ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
2159   
2160   if(!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
2161   
2162   if(IsExoticCluster(cluster, cells,bc))           return kFALSE;
2163
2164   return isGood;
2165 }
2166
2167 //__________________________________________________________
2168 Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
2169 {
2170   // Given a esd track, return whether the track survive all the cuts
2171
2172   // The different quality parameter are first
2173   // retrieved from the track. then it is found out what cuts the
2174   // track did not survive and finally the cuts are imposed.
2175
2176   UInt_t status = esdTrack->GetStatus();
2177
2178   Int_t nClustersITS = esdTrack->GetITSclusters(0);
2179   Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
2180
2181   Float_t chi2PerClusterITS = -1;
2182   Float_t chi2PerClusterTPC = -1;
2183   if (nClustersITS!=0)
2184     chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
2185   if (nClustersTPC!=0) 
2186     chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
2187
2188
2189   //DCA cuts
2190   if(fTrackCutsType==kGlobalCut)
2191     {
2192       Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
2193       //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
2194       SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
2195     }
2196
2197
2198   Float_t b[2];
2199   Float_t bCov[3];
2200   esdTrack->GetImpactParameters(b,bCov);
2201   if (bCov[0]<=0 || bCov[2]<=0) 
2202   {
2203     AliDebug(1, "Estimated b resolution lower or equal zero!");
2204     bCov[0]=0; bCov[2]=0;
2205   }
2206
2207   Float_t dcaToVertexXY = b[0];
2208   Float_t dcaToVertexZ = b[1];
2209   Float_t dcaToVertex = -1;
2210
2211   if (fCutDCAToVertex2D)
2212     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
2213   else
2214     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
2215     
2216   // cut the track?
2217   
2218   Bool_t cuts[kNCuts];
2219   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
2220   
2221   // track quality cuts
2222   if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
2223     cuts[0]=kTRUE;
2224   if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
2225     cuts[1]=kTRUE;
2226   if (nClustersTPC<fCutMinNClusterTPC)
2227     cuts[2]=kTRUE;
2228   if (nClustersITS<fCutMinNClusterITS) 
2229     cuts[3]=kTRUE;
2230   if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC) 
2231     cuts[4]=kTRUE; 
2232   if (chi2PerClusterITS>fCutMaxChi2PerClusterITS) 
2233     cuts[5]=kTRUE;  
2234   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
2235     cuts[6]=kTRUE;
2236   if (fCutDCAToVertex2D && dcaToVertex > 1)
2237     cuts[7] = kTRUE;
2238   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
2239     cuts[8] = kTRUE;
2240   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
2241     cuts[9] = kTRUE;
2242
2243   if(fTrackCutsType==kGlobalCut)
2244     {
2245       //Require at least one SPD point + anything else in ITS
2246       if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
2247   cuts[10] = kTRUE;
2248     }
2249
2250   Bool_t cut=kFALSE;
2251   for (Int_t i=0; i<kNCuts; i++)
2252     if (cuts[i]) { cut = kTRUE ; }
2253
2254     // cut the track
2255   if (cut) 
2256     return kFALSE;
2257   else 
2258     return kTRUE;
2259 }
2260
2261 //_____________________________________
2262 void AliEMCALRecoUtils::InitTrackCuts()
2263 {
2264   //Intilize the track cut criteria
2265   //By default these cuts are set according to AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()
2266   //Also you can customize the cuts using the setters
2267   
2268   switch (fTrackCutsType)
2269   {
2270     case kTPCOnlyCut:
2271     {
2272       AliInfo(Form("Track cuts for matching: GetStandardTPCOnlyTrackCuts()"));
2273       //TPC
2274       SetMinNClustersTPC(70);
2275       SetMaxChi2PerClusterTPC(4);
2276       SetAcceptKinkDaughters(kFALSE);
2277       SetRequireTPCRefit(kFALSE);
2278       
2279       //ITS
2280       SetRequireITSRefit(kFALSE);
2281       SetMaxDCAToVertexZ(3.2);
2282       SetMaxDCAToVertexXY(2.4);
2283       SetDCAToVertex2D(kTRUE);
2284       
2285       break;
2286     }
2287       
2288     case kGlobalCut:
2289     {
2290       AliInfo(Form("Track cuts for matching: GetStandardITSTPCTrackCuts2010(kTURE)"));
2291       //TPC
2292       SetMinNClustersTPC(70);
2293       SetMaxChi2PerClusterTPC(4);
2294       SetAcceptKinkDaughters(kFALSE);
2295       SetRequireTPCRefit(kTRUE);
2296       
2297       //ITS
2298       SetRequireITSRefit(kTRUE);
2299       SetMaxDCAToVertexZ(2);
2300       SetMaxDCAToVertexXY();
2301       SetDCAToVertex2D(kFALSE);
2302       
2303       break;
2304     }
2305       
2306     case kLooseCut:
2307     {
2308       AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut"));
2309       SetMinNClustersTPC(50);
2310       SetAcceptKinkDaughters(kTRUE);
2311       
2312       break;
2313     }
2314   }
2315 }
2316
2317
2318 //________________________________________________________________________
2319 void AliEMCALRecoUtils::SetClusterMatchedToTrack(const AliVEvent *event)
2320 {
2321   // Checks if tracks are matched to EMC clusters and set the matched EMCAL cluster index to ESD track. 
2322
2323   Int_t nTracks = event->GetNumberOfTracks();
2324   for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack) 
2325   {
2326     AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrack));
2327     if (!track) 
2328     {
2329       AliWarning(Form("Could not receive track %d", iTrack));
2330       continue;
2331     }
2332     
2333     Int_t matchClusIndex = GetMatchedClusterIndex(iTrack);       
2334     track->SetEMCALcluster(matchClusIndex); //sets -1 if track not matched within residual
2335     /*the following can be done better if AliVTrack::SetStatus will be there. Patch pending with Andreas/Peter*/
2336     AliESDtrack* esdtrack = dynamic_cast<AliESDtrack*>(track);
2337     if (esdtrack) { 
2338       if(matchClusIndex != -1) 
2339         esdtrack->SetStatus(AliESDtrack::kEMCALmatch);
2340       else
2341         esdtrack->ResetStatus(AliESDtrack::kEMCALmatch);
2342     } else {
2343       AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(track);
2344       if(matchClusIndex != -1) 
2345         aodtrack->SetStatus(AliESDtrack::kEMCALmatch);
2346       else
2347         aodtrack->ResetStatus(AliESDtrack::kEMCALmatch);
2348     }
2349
2350   }
2351   AliDebug(2,"Track matched to closest cluster");  
2352 }
2353
2354 //_________________________________________________________________________
2355 void AliEMCALRecoUtils::SetTracksMatchedToCluster(const AliVEvent *event)
2356 {
2357   // Checks if EMC clusters are matched to ESD track.
2358   // Adds track indexes of all the tracks matched to a cluster withing residuals in ESDCalocluster.
2359   
2360   for (Int_t iClus=0; iClus < event->GetNumberOfCaloClusters(); ++iClus) 
2361   {
2362     AliVCluster *cluster = event->GetCaloCluster(iClus);
2363     if (!cluster->IsEMCAL()) 
2364       continue;
2365     
2366     Int_t nTracks = event->GetNumberOfTracks();
2367     TArrayI arrayTrackMatched(nTracks);
2368     
2369     // Get the closest track matched to the cluster
2370     Int_t nMatched = 0;
2371     Int_t matchTrackIndex = GetMatchedTrackIndex(iClus);
2372     if (matchTrackIndex != -1) 
2373     {
2374       arrayTrackMatched[nMatched] = matchTrackIndex;
2375       nMatched++;
2376     }
2377     
2378     // Get all other tracks matched to the cluster
2379     for(Int_t iTrk=0; iTrk<nTracks; ++iTrk) 
2380     {
2381       AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrk));
2382       if(iTrk == matchTrackIndex) continue;
2383       if(track->GetEMCALcluster() == iClus)
2384       {
2385         arrayTrackMatched[nMatched] = iTrk;
2386         ++nMatched;
2387       }
2388     }
2389     
2390     //printf("Tender::SetTracksMatchedToCluster - cluster E %f, N matches %d, first match %d\n",cluster->E(),nMatched,arrayTrackMatched[0]);
2391     
2392     arrayTrackMatched.Set(nMatched);
2393     AliESDCaloCluster *esdcluster = dynamic_cast<AliESDCaloCluster*>(cluster);
2394     if (esdcluster) 
2395       esdcluster->AddTracksMatched(arrayTrackMatched);
2396     else if (nMatched>0) {
2397       AliAODCaloCluster *aodcluster = dynamic_cast<AliAODCaloCluster*>(cluster);
2398       if (aodcluster)
2399         aodcluster->AddTrackMatched(event->GetTrack(arrayTrackMatched.At(0)));
2400     }
2401     
2402     Float_t eta= -999, phi = -999;
2403     if (matchTrackIndex != -1) 
2404       GetMatchedResiduals(iClus, eta, phi);
2405     cluster->SetTrackDistance(phi, eta);
2406   }
2407   
2408   AliDebug(2,"Cluster matched to tracks");  
2409 }
2410
2411
2412 //___________________________________________________
2413 void AliEMCALRecoUtils::Print(const Option_t *) const 
2414 {
2415   // Print Parameters
2416   
2417   printf("AliEMCALRecoUtils Settings: \n");
2418   printf("Misalignment shifts\n");
2419   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, 
2420                                   fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
2421                                   fMisalRotShift[i*3],  fMisalRotShift[i*3+1],  fMisalRotShift[i*3+2]   );
2422   printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
2423   for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
2424   
2425   printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
2426
2427   printf("Matching criteria: ");
2428   if(fCutEtaPhiSum)
2429     {
2430       printf("sqrt(dEta^2+dPhi^2)<%4.3f\n",fCutR);
2431     }
2432   else if(fCutEtaPhiSeparate)
2433     {
2434       printf("dEta<%4.3f, dPhi<%4.3f\n",fCutEta,fCutPhi);
2435     }
2436   else
2437     {
2438       printf("Error\n");
2439       printf("please specify your cut criteria\n");
2440       printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
2441       printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
2442     }
2443
2444   printf("Mass hypothesis = %2.3f [GeV/c^2], extrapolation step to surface = %2.2f[cm], step to cluster = %2.2f[cm]\n",fMass,fStepSurface, fStepCluster);
2445   printf("Cluster selection window: dR < %2.0f\n",fClusterWindow);
2446
2447   printf("Track cuts: \n");
2448   printf("Minimum track pT: %1.2f\n",fCutMinTrackPt);
2449   printf("AOD track selection mask: %d\n",fAODFilterMask);
2450   printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
2451   printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
2452   printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
2453   printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
2454   printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
2455 }
2456
2457 //_________________________________________________________________
2458 void AliEMCALRecoUtils::SetRunDependentCorrections(Int_t runnumber)
2459 {
2460   //Get EMCAL time dependent corrections from file and put them in the recalibration histograms
2461   //Do it only once and only if it is requested
2462   
2463   if(!fUseRunCorrectionFactors) return;
2464   if(fRunCorrectionFactorsSet)  return;
2465   
2466   AliInfo(Form("AliEMCALRecoUtils::GetRunDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber));
2467  
2468   AliEMCALCalibTimeDepCorrection  *corr =  new AliEMCALCalibTimeDepCorrection();
2469   corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber));
2470   
2471   SwitchOnRecalibration();
2472   
2473   AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
2474   
2475   for(Int_t ism = 0; ism < geom->GetNumberOfSuperModules(); ism++)
2476   {
2477     for(Int_t icol = 0; icol < 48; icol++)
2478     {
2479       for(Int_t irow = 0; irow < 24; irow++)
2480       {
2481         Float_t orgRecalFactor = GetEMCALChannelRecalibrationFactors(ism)->GetBinContent(icol,irow);
2482         Float_t newRecalFactor = orgRecalFactor*corr->GetCorrection(ism, icol,irow,0);
2483         GetEMCALChannelRecalibrationFactors(ism)->SetBinContent(icol,irow,newRecalFactor);
2484         //printf("ism %d, icol %d, irow %d, corrections : org %f, time dep %f, final %f (org*time %f)\n",ism, icol, irow, 
2485         //        orgRecalFactor, corr->GetCorrection(ism, icol,irow,0),
2486         //       (GetEMCALChannelRecalibrationFactors(ism))->GetBinContent(icol,irow),newRecalFactor);
2487       }
2488     }
2489   }
2490   
2491   fRunCorrectionFactorsSet = kTRUE;
2492   
2493 }