Store and recover the index of the matched track, possibility to do matching with...
[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
29 // --- standard c ---
30
31 // standard C++ includes
32 //#include <Riostream.h>
33
34 // ROOT includes
35 #include <TGeoManager.h>
36 #include <TGeoMatrix.h>
37 #include <TGeoBBox.h>
38
39 // STEER includes
40 #include "AliVCluster.h"
41 #include "AliVCaloCells.h"
42 #include "AliVEvent.h"
43 #include "AliLog.h"
44 #include "AliPID.h"
45 #include "AliESDEvent.h"
46 #include "AliESDtrack.h"
47
48 // EMCAL includes
49 #include "AliEMCALRecoUtils.h"
50 #include "AliEMCALGeometry.h"
51 #include "AliEMCALTrack.h"
52 #include "AliEMCALCalibTimeDepCorrection.h"
53 #include "AliEMCALPIDUtils.h"
54
55 ClassImp(AliEMCALRecoUtils)
56   
57 //______________________________________________
58 AliEMCALRecoUtils::AliEMCALRecoUtils():
59   fNonLinearityFunction (kNoCorrection), fParticleType(kPhoton),
60   fPosAlgo(kUnchanged),fW0(4.),
61   fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
62   fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(),
63   fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
64   fMatchedTrackIndex(0x0), fMatchedClusterIndex(0x0), 
65   fResidualZ(0x0), fResidualR(0x0), fCutR(20), fCutZ(20),
66   fCutMinNClusterTPC(0), fCutMinNClusterITS(0), fCutMaxChi2PerClusterTPC(0), fCutMaxChi2PerClusterITS(0),
67   fCutRequireTPCRefit(0), fCutRequireITSRefit(0), fCutAcceptKinkDaughters(0),
68   fCutMaxDCAToVertexXY(0), fCutMaxDCAToVertexZ(0),fCutDCAToVertex2D(0),fPIDUtils(),
69   fUseTimeCorrectionFactors(kFALSE),  fTimeCorrectionFactorsSet(kFALSE)
70 {
71 //
72   // Constructor.
73   // Initialize all constant values which have to be used
74   // during Reco algorithm execution
75   //
76   
77   //Misalignment matrices
78   for(Int_t i = 0; i < 15 ; i++) {
79       fMisalTransShift[i] = 0.; 
80       fMisalRotShift[i]   = 0.; 
81   }
82   
83   //Non linearity
84   for(Int_t i = 0; i < 6  ; i++) fNonLinearityParams[i] = 0.; 
85   //For kPi0GammaGamma case, but default is no correction
86   fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
87   fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
88   fNonLinearityParams[2] = 1.046;
89
90   //Track matching
91   fMatchedTrackIndex   = new TArrayI();
92   fMatchedClusterIndex = new TArrayI();
93   fResidualZ           = new TArrayF();
94   fResidualR           = new TArrayF();
95   
96   InitTrackCuts();
97   
98   fPIDUtils            = new AliEMCALPIDUtils();
99
100
101 }
102
103 //______________________________________________________________________
104 AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco) 
105 : TNamed(reco), fNonLinearityFunction(reco.fNonLinearityFunction), 
106   fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo), fW0(reco.fW0), 
107   fRecalibration(reco.fRecalibration),fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
108   fRemoveBadChannels(reco.fRemoveBadChannels),fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
109   fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
110   fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
111   fMatchedTrackIndex(reco.fMatchedTrackIndex?new TArrayI(*reco.fMatchedTrackIndex):0x0),
112   fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
113   fResidualZ(reco.fResidualZ?new TArrayF(*reco.fResidualZ):0x0),
114   fResidualR(reco.fResidualR?new TArrayF(*reco.fResidualR):0x0),
115   fCutR(reco.fCutR),fCutZ(reco.fCutZ),
116   fCutMinNClusterTPC(reco.fCutMinNClusterTPC), fCutMinNClusterITS(reco.fCutMinNClusterITS), 
117   fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC), fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS),
118   fCutRequireTPCRefit(reco.fCutRequireTPCRefit), fCutRequireITSRefit(reco.fCutRequireITSRefit),
119   fCutAcceptKinkDaughters(reco.fCutAcceptKinkDaughters),
120   fCutMaxDCAToVertexXY(reco.fCutMaxDCAToVertexXY), fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ),fCutDCAToVertex2D(reco.fCutDCAToVertex2D),
121   fPIDUtils(reco.fPIDUtils), 
122   fUseTimeCorrectionFactors(reco.fUseTimeCorrectionFactors),  fTimeCorrectionFactorsSet(reco.fTimeCorrectionFactorsSet)
123 {
124   //Copy ctor
125   
126   for(Int_t i = 0; i < 15 ; i++) {
127       fMisalRotShift[i] = reco.fMisalRotShift[i]; 
128       fMisalTransShift[i] = reco.fMisalTransShift[i]; 
129   } 
130   for(Int_t i = 0; i < 6  ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i]; 
131
132 }
133
134
135 //______________________________________________________________________
136 AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco) 
137 {
138   //Assignment operator
139   
140   if(this == &reco)return *this;
141   ((TNamed *)this)->operator=(reco);
142
143   fNonLinearityFunction      = reco.fNonLinearityFunction;
144   fParticleType              = reco.fParticleType;
145   fPosAlgo                   = reco.fPosAlgo; 
146   fW0                        = reco.fW0;
147   fRecalibration             = reco.fRecalibration;
148   fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
149   fRemoveBadChannels         = reco.fRemoveBadChannels;
150   fRecalDistToBadChannels    = reco.fRecalDistToBadChannels;
151   fEMCALBadChannelMap        = reco.fEMCALBadChannelMap;
152   fNCellsFromEMCALBorder     = reco.fNCellsFromEMCALBorder;
153   fNoEMCALBorderAtEta0       = reco.fNoEMCALBorderAtEta0;
154
155
156   for(Int_t i = 0; i < 15 ; i++) {fMisalTransShift[i] = reco.fMisalTransShift[i]; fMisalRotShift[i] = reco.fMisalRotShift[i];}
157   for(Int_t i = 0; i < 6  ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i]; 
158   
159   fCutR                      = reco.fCutR;
160   fCutZ                      = reco.fCutZ;
161
162   fCutMinNClusterTPC         = reco.fCutMinNClusterTPC;
163   fCutMinNClusterITS         = reco.fCutMinNClusterITS; 
164   fCutMaxChi2PerClusterTPC   = reco.fCutMaxChi2PerClusterTPC;
165   fCutMaxChi2PerClusterITS   = reco.fCutMaxChi2PerClusterITS;
166   fCutRequireTPCRefit        = reco.fCutRequireTPCRefit;
167   fCutRequireITSRefit        = reco.fCutRequireITSRefit;
168   fCutAcceptKinkDaughters    = reco.fCutAcceptKinkDaughters;
169   fCutMaxDCAToVertexXY       = reco.fCutMaxDCAToVertexXY;
170   fCutMaxDCAToVertexZ        = reco.fCutMaxDCAToVertexZ;
171   fCutDCAToVertex2D          = reco.fCutDCAToVertex2D;
172
173   fPIDUtils                  = reco.fPIDUtils;
174   
175   fUseTimeCorrectionFactors  = reco.fUseTimeCorrectionFactors;
176   fTimeCorrectionFactorsSet  = reco.fTimeCorrectionFactorsSet;
177
178   
179   if(reco.fResidualR){
180     // assign or copy construct
181     if(fResidualR){ 
182       *fResidualR = *reco.fResidualR;
183     }
184     else fResidualR = new TArrayF(*reco.fResidualR);
185   }
186   else{
187     if(fResidualR)delete fResidualR;
188     fResidualR = 0;
189   }
190   
191   if(reco.fResidualZ){
192     // assign or copy construct
193     if(fResidualZ){ 
194       *fResidualZ = *reco.fResidualZ;
195     }
196     else fResidualZ = new TArrayF(*reco.fResidualZ);
197   }
198   else{
199     if(fResidualZ)delete fResidualZ;
200     fResidualZ = 0;
201   }
202   
203   if(reco.fMatchedTrackIndex){
204     // assign or copy construct
205     if(fMatchedTrackIndex){ 
206       *fMatchedTrackIndex = *reco.fMatchedTrackIndex;
207     }
208     else fMatchedTrackIndex = new TArrayI(*reco.fMatchedTrackIndex);
209   }
210   else{
211     if(fMatchedTrackIndex)delete fMatchedTrackIndex;
212     fMatchedTrackIndex = 0;
213   }  
214   
215   if(reco.fMatchedClusterIndex){
216     // assign or copy construct
217     if(fMatchedClusterIndex){ 
218       *fMatchedClusterIndex = *reco.fMatchedClusterIndex;
219     }
220     else fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex);
221   }
222   else{
223     if(fMatchedClusterIndex)delete fMatchedClusterIndex;
224     fMatchedClusterIndex = 0;
225   }
226   
227   
228   return *this;
229 }
230
231
232 //__________________________________________________
233 AliEMCALRecoUtils::~AliEMCALRecoUtils()
234 {
235   //Destructor.
236         
237         if(fEMCALRecalibrationFactors) { 
238                 fEMCALRecalibrationFactors->Clear();
239                 delete  fEMCALRecalibrationFactors;
240         }       
241   
242   if(fEMCALBadChannelMap) { 
243                 fEMCALBadChannelMap->Clear();
244                 delete  fEMCALBadChannelMap;
245         }
246  
247   if(fMatchedTrackIndex)   {delete fMatchedTrackIndex;   fMatchedTrackIndex=0;}
248   if(fMatchedClusterIndex) {delete fMatchedClusterIndex; fMatchedClusterIndex=0;}
249   if(fResidualR)           {delete fResidualR;           fResidualR=0;}
250   if(fResidualZ)           {delete fResidualZ;           fResidualZ=0;}
251
252 }
253
254 //_______________________________________________________________
255 Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells) 
256 {
257         // Given the list of AbsId of the cluster, get the maximum cell and 
258         // check if there are fNCellsFromBorder from the calorimeter border
259         
260   //If the distance to the border is 0 or negative just exit accept all clusters
261         if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
262   
263   Int_t absIdMax        = -1, iSM =-1, ieta = -1, iphi = -1;
264   Bool_t shared = kFALSE;
265   GetMaxEnergyCell(geom, cells, cluster, absIdMax,  iSM, ieta, iphi, shared);
266
267   AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n", 
268            absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
269         
270         if(absIdMax==-1) return kFALSE;
271         
272         //Check if the cell is close to the borders:
273         Bool_t okrow = kFALSE;
274         Bool_t okcol = kFALSE;
275   
276   if(iSM < 0 || iphi < 0 || ieta < 0 ) {
277     AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
278                   iSM,ieta,iphi));
279   }
280   
281   //Check rows/phi
282   if(iSM < 10){
283     if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE; 
284   }
285   else{
286     if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE; 
287   }
288   
289   //Check columns/eta
290   if(!fNoEMCALBorderAtEta0){
291     if(ieta  > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE; 
292   }
293   else{
294     if(iSM%2==0){
295       if(ieta >= fNCellsFromEMCALBorder)     okcol = kTRUE;     
296     }
297     else {
298       if(ieta <  48-fNCellsFromEMCALBorder)  okcol = kTRUE;     
299     }
300   }//eta 0 not checked
301     
302   AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d:  column? %d, row? %d\nq",
303            fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
304         
305         if (okcol && okrow) {
306     //printf("Accept\n");
307     return kTRUE;
308   }
309         else  {
310     //printf("Reject\n");
311     AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
312     return kFALSE;
313   }
314         
315 }       
316
317
318 //_________________________________________________________________________________________________________
319 Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, Int_t nCells){
320         // Check that in the cluster cells, there is no bad channel of those stored 
321         // in fEMCALBadChannelMap or fPHOSBadChannelMap
322         
323         if(!fRemoveBadChannels)  return kFALSE;
324         if(!fEMCALBadChannelMap) return kFALSE;
325         
326         Int_t icol = -1;
327         Int_t irow = -1;
328         Int_t imod = -1;
329         for(Int_t iCell = 0; iCell<nCells; iCell++){
330                 
331                 //Get the column and row
332     Int_t iTower = -1, iIphi = -1, iIeta = -1; 
333     geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta); 
334     if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
335     geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);                      
336     if(GetEMCALChannelStatus(imod, icol, irow)){
337       AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
338       return kTRUE;
339     }
340                 
341         }// cell cluster loop
342         
343         return kFALSE;
344         
345 }
346
347 //__________________________________________________
348 Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
349 // Correct cluster energy from non linearity functions
350   Float_t energy = cluster->E();
351   
352   switch (fNonLinearityFunction) {
353       
354     case kPi0MC:
355     {
356       //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
357       //Double_t fNonLinearityParams[0] = 1.001;
358       //Double_t fNonLinearityParams[1] = -0.01264;
359       //Double_t fNonLinearityParams[2] = -0.03632;
360       //Double_t fNonLinearityParams[3] = 0.1798;
361       //Double_t fNonLinearityParams[4] = -0.522;
362        energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
363                   ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
364                     exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
365       break;
366     }
367       
368     case kPi0GammaGamma:
369     {
370       //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
371       //Double_t fNonLinearityParams[0] = 1.04;
372       //Double_t fNonLinearityParams[1] = -0.1445;
373       //Double_t fNonLinearityParams[2] = 1.046;
374       energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
375       break;
376     }
377       
378     case kPi0GammaConversion:
379     {
380       //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
381       //fNonLinearityParams[0] = 0.139393/0.1349766;
382       //fNonLinearityParams[1] = 0.0566186;
383       //fNonLinearityParams[2] = 0.982133;
384       energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
385       
386       break;
387     }
388       
389     case kBeamTest:
390     {
391       //From beam test, Alexei's results, for different ZS thresholds
392       //                        th=30 MeV; th = 45 MeV; th = 75 MeV
393       //fNonLinearityParams[0] = 1.007;      1.003;      1.002 
394       //fNonLinearityParams[1] = 0.894;      0.719;      0.797 
395       //fNonLinearityParams[2] = 0.246;      0.334;      0.358 
396       //Rescale the param[0] with 1.03
397       energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
398       
399       break;
400     }
401       
402     case kNoCorrection:
403       AliDebug(2,"No correction on the energy\n");
404       break;
405       
406   }
407   
408   return energy;
409
410 }
411
412 //__________________________________________________
413 Float_t  AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const 
414 {
415   //Calculate shower depth for a given cluster energy and particle type
416
417   // parameters 
418   Float_t x0    = 1.31;
419   Float_t ecr   = 8;
420   Float_t depth = 0;
421   
422   switch ( iParticle )
423   {
424     case kPhoton:
425       depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
426       break;
427       
428     case kElectron:
429       depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
430       break;
431       
432     case kHadron:
433       // hadron 
434       // boxes anc. here
435       if(gGeoManager){
436         gGeoManager->cd("ALIC_1/XEN1_1");
437         TGeoNode        *geoXEn1    = gGeoManager->GetCurrentNode();
438         TGeoNodeMatrix  *geoSM      = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
439         if(geoSM){
440           TGeoVolume      *geoSMVol   = geoSM->GetVolume(); 
441           TGeoShape       *geoSMShape = geoSMVol->GetShape();
442           TGeoBBox        *geoBox     = dynamic_cast<TGeoBBox *>(geoSMShape);
443           if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
444           else AliFatal("Null GEANT box");
445         }else AliFatal("NULL  GEANT node matrix");
446       }
447       else{//electron
448         depth = x0 * (TMath::Log(energy*1000 / ecr)  - 0.5); //Multiply energy by 1000 to transform to MeV
449       }
450         
451       break;
452       
453     default://photon
454       depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
455   }  
456   
457   return depth;
458   
459 }
460
461 //__________________________________________________
462 void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu, 
463                                          Int_t & absId,  Int_t& iSupMod, Int_t& ieta, Int_t& iphi, Bool_t &shared)
464 {
465   //For a given CaloCluster gets the absId of the cell 
466   //with maximum energy deposit.
467   
468   Double_t eMax        = -1.;
469   Double_t eCell       = -1.;
470   Float_t  fraction    = 1.;
471   Float_t  recalFactor = 1.;
472   Int_t    cellAbsId   = -1 ;
473
474   Int_t iTower  = -1;
475   Int_t iIphi   = -1;
476   Int_t iIeta   = -1;
477   Int_t iSupMod0= -1;
478         //printf("---Max?\n");
479   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
480     cellAbsId = clu->GetCellAbsId(iDig);
481     fraction  = clu->GetCellAmplitudeFraction(iDig);
482     //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
483     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
484     geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta); 
485     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
486     if(iDig==0) iSupMod0=iSupMod;
487     else if(iSupMod0!=iSupMod) {
488       shared = kTRUE;
489       //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
490     }
491     if(IsRecalibrationOn()) {
492       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
493     }
494     eCell  = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
495     //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
496     if(eCell > eMax)  { 
497       eMax  = eCell; 
498       absId = cellAbsId;
499       //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
500     }
501   }// cell loop
502   
503   //Get from the absid the supermodule, tower and eta/phi numbers
504   geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta); 
505   //Gives SuperModule and Tower numbers
506   geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
507                                          iIphi, iIeta,iphi,ieta); 
508   //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
509   //printf("Max end---\n");
510   
511 }
512
513 //________________________________________________________________
514 void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
515         //Init EMCAL recalibration factors
516         AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
517         //In order to avoid rewriting the same histograms
518         Bool_t oldStatus = TH1::AddDirectoryStatus();
519         TH1::AddDirectory(kFALSE);
520   
521         fEMCALRecalibrationFactors = new TObjArray(10);
522         for (int i = 0; i < 12; i++) fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i),  48, 0, 48, 24, 0, 24));
523         //Init the histograms with 1
524         for (Int_t sm = 0; sm < 12; sm++) {
525                 for (Int_t i = 0; i < 48; i++) {
526                         for (Int_t j = 0; j < 24; j++) {
527                                 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
528                         }
529                 }
530         }
531         fEMCALRecalibrationFactors->SetOwner(kTRUE);
532         fEMCALRecalibrationFactors->Compress();
533         
534         //In order to avoid rewriting the same histograms
535         TH1::AddDirectory(oldStatus);           
536 }
537
538
539 //________________________________________________________________
540 void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
541         //Init EMCAL bad channels map
542         AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
543         //In order to avoid rewriting the same histograms
544         Bool_t oldStatus = TH1::AddDirectoryStatus();
545         TH1::AddDirectory(kFALSE);
546         
547         fEMCALBadChannelMap = new TObjArray(10);
548         //TH2F * hTemp = new  TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
549         for (int i = 0; i < 10; i++) {
550                 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
551         }
552         
553         //delete hTemp;
554         
555         fEMCALBadChannelMap->SetOwner(kTRUE);
556         fEMCALBadChannelMap->Compress();
557         
558         //In order to avoid rewriting the same histograms
559         TH1::AddDirectory(oldStatus);           
560 }
561
562 //________________________________________________________________
563 void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells){
564         // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
565         
566         //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
567         UShort_t * index    = cluster->GetCellsAbsId() ;
568         Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
569         Int_t ncells = cluster->GetNCells();
570         
571         //Initialize some used variables
572         Float_t energy = 0;
573         Int_t absId    = -1;
574   Int_t icol = -1, irow = -1, imod=1;
575         Float_t factor = 1, frac = 0;
576         
577         //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
578         for(Int_t icell = 0; icell < ncells; icell++){
579                 absId = index[icell];
580                 frac =  fraction[icell];
581                 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
582                 Int_t iTower = -1, iIphi = -1, iIeta = -1; 
583                 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta); 
584                 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
585                 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);                  
586                 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
587     AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
588              imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
589                 
590                 energy += cells->GetCellAmplitude(absId)*factor*frac;
591         }
592         
593         
594                 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
595         
596         cluster->SetE(energy);
597         
598 }
599
600
601 //__________________________________________________
602 void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
603 {
604   //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
605   
606   if     (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
607   else if(fPosAlgo==kPosTowerIndex)  RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
608   else   AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
609   
610 }  
611
612 //__________________________________________________
613 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
614 {
615   // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
616   // The algorithm is a copy of what is done in AliEMCALRecPoint
617   
618   Double_t eCell       = 0.;
619   Float_t  fraction    = 1.;
620   Float_t  recalFactor = 1.;
621   
622   Int_t    absId   = -1;
623   Int_t    iTower  = -1, iIphi  = -1, iIeta  = -1;
624   Int_t    iSupModMax = -1, iSM=-1, iphi   = -1, ieta   = -1;
625   Float_t  weight = 0.,  totalWeight=0.;
626   Float_t  newPos[3] = {0,0,0};
627   Double_t pLocal[3], pGlobal[3];
628   Bool_t shared = kFALSE;
629
630   Float_t  clEnergy = clu->E(); //Energy already recalibrated previously
631   GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi,shared);
632   Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
633   
634   //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
635   
636   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
637     absId = clu->GetCellAbsId(iDig);
638     fraction  = clu->GetCellAmplitudeFraction(iDig);
639     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
640     geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta); 
641     geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);                       
642     
643     if(IsRecalibrationOn()) {
644       recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
645     }
646     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
647     
648     weight = GetCellWeight(eCell,clEnergy);
649     //printf("cell energy %f, weight %f\n",eCell,weight);
650     totalWeight += weight;
651     geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
652     //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
653     geom->GetGlobal(pLocal,pGlobal,iSupModMax);
654     //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
655
656     for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
657     
658   }// cell loop
659   
660   if(totalWeight>0){
661     for(int i=0; i<3; i++ )    newPos[i] /= totalWeight;
662   }
663     
664   //Float_t pos[]={0,0,0};
665   //clu->GetPosition(pos);
666   //printf("OldPos  : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
667   //printf("NewPos  : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
668   
669         if(iSupModMax > 1) {//sector 1
670           newPos[0] +=fMisalTransShift[3];//-=3.093; 
671           newPos[1] +=fMisalTransShift[4];//+=6.82;
672           newPos[2] +=fMisalTransShift[5];//+=1.635;
673     //printf("   +    : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
674
675         }
676         else {//sector 0
677           newPos[0] +=fMisalTransShift[0];//+=1.134;
678           newPos[1] +=fMisalTransShift[1];//+=8.2;
679           newPos[2] +=fMisalTransShift[2];//+=1.197;
680     //printf("   +    : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
681
682         }
683   //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
684
685   clu->SetPosition(newPos);
686   
687 }  
688
689 //__________________________________________________
690 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
691 {
692   // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
693   // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
694   
695   Double_t eCell       = 1.;
696   Float_t  fraction    = 1.;
697   Float_t  recalFactor = 1.;
698   
699   Int_t absId   = -1;
700   Int_t iTower  = -1;
701   Int_t iIphi   = -1, iIeta   = -1;
702         Int_t iSupMod = -1, iSupModMax = -1;
703   Int_t iphi = -1, ieta =-1;
704   Bool_t shared = kFALSE;
705
706   Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
707   GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi,shared);
708   Float_t  depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
709
710   Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
711   Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
712   Int_t startingSM = -1;
713   
714   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
715     absId = clu->GetCellAbsId(iDig);
716     fraction  = clu->GetCellAmplitudeFraction(iDig);
717     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
718     geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta); 
719     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);                   
720     
721     if     (iDig==0)  startingSM = iSupMod;
722     else if(iSupMod != startingSM) areInSameSM = kFALSE;
723
724     eCell  = cells->GetCellAmplitude(absId);
725     
726     if(IsRecalibrationOn()) {
727       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
728     }
729     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
730     
731     weight = GetCellWeight(eCell,clEnergy);
732     if(weight < 0) weight = 0;
733     totalWeight += weight;
734     weightedCol += ieta*weight;
735     weightedRow += iphi*weight;
736     
737     //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
738     
739     }// cell loop
740     
741   Float_t xyzNew[]={0.,0.,0.};
742   if(areInSameSM == kTRUE) {
743     //printf("In Same SM\n");
744     weightedCol = weightedCol/totalWeight;
745     weightedRow = weightedRow/totalWeight;
746     geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew); 
747   }
748   else {
749     //printf("In Different SM\n");
750     geom->RecalculateTowerPosition(iphi,        ieta,        iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew); 
751   }
752   
753   clu->SetPosition(xyzNew);
754   
755 }
756
757 //____________________________________________________________________________
758 void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster){           
759         
760   //re-evaluate distance to bad channel with updated bad map
761   
762   if(!fRecalDistToBadChannels) return;
763   
764         //Get channels map of the supermodule where the cluster is.
765   Int_t absIdMax        = -1, iSupMod =-1, icolM = -1, irowM = -1;
766   Bool_t shared = kFALSE;
767   GetMaxEnergyCell(geom, cells, cluster, absIdMax,  iSupMod, icolM, irowM, shared);
768   TH2D* hMap  = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
769
770   Int_t dRrow, dRcol;   
771         Float_t  minDist = 10000.;
772         Float_t  dist    = 0.;
773   
774   //Loop on tower status map 
775         for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
776                 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
777                         //Check if tower is bad.
778                         if(hMap->GetBinContent(icol,irow)==0) continue;
779       //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
780       //       iSupMod,icol, irow, icolM,irowM);
781       
782       dRrow=TMath::Abs(irowM-irow);
783       dRcol=TMath::Abs(icolM-icol);
784       dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
785                         if(dist < minDist){
786         //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
787         minDist = dist;
788       }
789       
790                 }
791         }
792   
793         //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
794         if (shared) {
795                 TH2D* hMap2 = 0;
796                 Int_t iSupMod2 = -1;
797     
798                 //The only possible combinations are (0,1), (2,3) ... (8,9)
799                 if(iSupMod%2) iSupMod2 = iSupMod-1;
800                 else          iSupMod2 = iSupMod+1;
801                 hMap2  = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
802     
803                 //Loop on tower status map of second super module
804                 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
805                         for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
806                                 //Check if tower is bad.
807                                 if(hMap2->GetBinContent(icol,irow)==0) continue;
808                                 //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",
809           //     iSupMod2,icol, irow,iSupMod,icolM,irowM);
810         
811         dRrow=TMath::Abs(irow-irowM);
812         
813         if(iSupMod%2) {
814                                   dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
815                                 }
816         else {
817           dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
818                                 }                    
819         
820                                 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
821         if(dist < minDist) minDist = dist;        
822         
823                         }
824                 }
825     
826         }// shared cluster in 2 SuperModules
827   
828   AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
829   cluster->SetDistanceToBadChannel(minDist);
830   
831 }
832
833 //____________________________________________________________________________
834 void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){           
835         
836   //re-evaluate identification parameters with bayesian
837
838         if ( cluster->GetM02() != 0)
839     fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
840   
841   Float_t pidlist[AliPID::kSPECIESN+1];
842         for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
843   
844   cluster->SetPID(pidlist);
845         
846 }
847
848 //____________________________________________________________________________
849 void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster)
850 {
851   // Calculates new center of gravity in the local EMCAL-module coordinates 
852   // and tranfers into global ALICE coordinates
853   // Calculates Dispersion and main axis
854   
855   Int_t nstat  = 0;
856   Float_t wtot = 0. ;
857   Double_t eCell       = 0.;
858   Float_t  fraction    = 1.;
859   Float_t  recalFactor = 1.;
860
861   Int_t iSupMod = -1;
862   Int_t iTower  = -1;
863   Int_t iIphi   = -1;
864   Int_t iIeta   = -1;
865   Int_t iphi    = -1;
866   Int_t ieta    = -1;
867   Double_t etai = -1.;
868   Double_t phii = -1.;
869   
870   Double_t w     = 0.;
871   Double_t d     = 0.;
872   Double_t dxx   = 0.;
873   Double_t dzz   = 0.;
874   Double_t dxz   = 0.;  
875   Double_t xmean = 0.;
876   Double_t zmean = 0.;
877     
878   //Loop on cells
879   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
880     
881     //Get from the absid the supermodule, tower and eta/phi numbers
882     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta); 
883     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);        
884     
885     //Get the cell energy, if recalibration is on, apply factors
886     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
887     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
888     if(IsRecalibrationOn()) {
889       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
890     }
891     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
892     
893     if(cluster->E() > 0 && eCell > 0){
894       
895       w  = GetCellWeight(eCell,cluster->E());
896       
897       etai=(Double_t)ieta;
898       phii=(Double_t)iphi;              
899       if(w > 0.0) {
900         wtot += w ;
901         nstat++;                        
902         //Shower shape
903         dxx  += w * etai * etai ;
904         xmean+= w * etai ;
905         dzz  += w * phii * phii ;
906         zmean+= w * phii ; 
907         dxz  += w * etai * phii ; 
908       }
909     }
910     else
911       AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
912   }//cell loop
913   
914   //Normalize to the weight     
915   if (wtot > 0) {
916     xmean /= wtot ;
917     zmean /= wtot ;
918   }
919   else
920     AliError(Form("Wrong weight %f\n", wtot));
921   
922   //Calculate dispersion        
923   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
924     
925     //Get from the absid the supermodule, tower and eta/phi numbers
926     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta); 
927     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
928     
929     //Get the cell energy, if recalibration is on, apply factors
930     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
931     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
932     if(IsRecalibrationOn()) {
933       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
934     }
935     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
936     
937     if(cluster->E() > 0 && eCell > 0){
938       
939       w  = GetCellWeight(eCell,cluster->E());
940       
941       etai=(Double_t)ieta;
942       phii=(Double_t)iphi;              
943       if(w > 0.0)  d +=  w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean)); 
944     }
945     else
946       AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
947   }// cell loop
948   
949   //Normalize to the weigth and set shower shape parameters
950   if (wtot > 0 && nstat > 1) {
951     d /= wtot ;
952     dxx /= wtot ;
953     dzz /= wtot ;
954     dxz /= wtot ;
955     dxx -= xmean * xmean ;
956     dzz -= zmean * zmean ;
957     dxz -= xmean * zmean ;
958     cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
959     cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
960   }
961   else{
962     d=0. ;
963     cluster->SetM20(0.) ;
964     cluster->SetM02(0.) ;
965   }     
966   
967   if (d>=0)
968     cluster->SetDispersion(TMath::Sqrt(d)) ;
969   else    
970     cluster->SetDispersion(0) ;
971   
972 }
973
974 //____________________________________________________________________________
975 void AliEMCALRecoUtils::FindMatches(AliVEvent *event, TObjArray * clusterArr)
976 {
977   //This function should be called before the cluster loop
978   //Before call this function, please recalculate the cluster positions
979   //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
980   //Store matched cluster indexes and residuals
981   //It only works with ESDs, not AODs
982   
983   fMatchedTrackIndex  ->Reset();
984   fMatchedClusterIndex->Reset();
985   fResidualZ          ->Reset();
986   fResidualR          ->Reset();
987   
988   fMatchedTrackIndex  ->Set(100);
989   fMatchedClusterIndex->Set(100);
990   fResidualZ          ->Set(100);
991   fResidualR          ->Set(100);
992   
993   Int_t    matched=0;
994   Float_t  clsPos[3];
995   Double_t trkPos[3];
996   for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
997   {
998     AliESDtrack *track = ((AliESDEvent*)event)->GetTrack(itr);
999     if(!track || !IsAccepted(track)) continue;
1000     
1001     Float_t dRMax = fCutR, dZMax = fCutZ;
1002     Int_t index = -1;
1003     AliEMCALTrack *emctrack = new AliEMCALTrack(*track);
1004     if(!clusterArr){// get clusters from event
1005       for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1006       {
1007         AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1008         if(!cluster->IsEMCAL()) continue;
1009         cluster->GetPosition(clsPos); //Has been recalculated
1010         if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) )  continue;
1011         emctrack->GetXYZ(trkPos);
1012         Float_t tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
1013         Float_t tmpZ = TMath::Abs(clsPos[2]-trkPos[2]);
1014       
1015         if(tmpR<dRMax)
1016         {
1017           dRMax=tmpR;
1018           dZMax=tmpZ;
1019           index=icl;
1020         }
1021       }//cluster loop
1022     } else { // external cluster array, not from ESD event
1023       for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
1024       {
1025         AliVCluster *cluster = (AliVCluster*) clusterArr->At(icl);
1026         if(!cluster->IsEMCAL()) continue;
1027         cluster->GetPosition(clsPos); //Has been recalculated
1028         if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) )  continue;
1029         emctrack->GetXYZ(trkPos);
1030         Float_t tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
1031         Float_t tmpZ = TMath::Abs(clsPos[2]-trkPos[2]);
1032         
1033         if(tmpR<dRMax)
1034         {
1035           dRMax=tmpR;
1036           dZMax=tmpZ;
1037           index=icl;
1038         }
1039       }//cluster loop
1040     }// external list of clusters
1041
1042     if(index>-1)
1043     {
1044       fMatchedTrackIndex  ->AddAt(itr,matched);
1045       fMatchedClusterIndex->AddAt(index,matched);
1046       fResidualZ          ->AddAt(dZMax,matched);
1047       fResidualR          ->AddAt(dRMax,matched);
1048       matched++;
1049     }
1050     delete emctrack;
1051   }//track loop
1052   
1053   AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1054   
1055   fMatchedTrackIndex  ->Set(matched);
1056   fMatchedClusterIndex->Set(matched);
1057   fResidualZ          ->Set(matched);
1058   fResidualR          ->Set(matched);
1059   
1060   //printf("Number of matched pairs: %d\n",matched);
1061 }
1062
1063 //________________________________________________________________________________
1064 void AliEMCALRecoUtils::GetMatchedResiduals(Int_t index, Float_t &dR, Float_t &dZ)
1065 {
1066   //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1067   //Get the residuals dR and dZ for this cluster
1068   //It only works with ESDs, not AODs
1069
1070   if( FindMatchedPos(index) >= 999 )
1071   {
1072     AliDebug(2,"No matched tracks found!\n");
1073     dR=999.;
1074     dZ=999.;
1075     return;
1076   }
1077   dR = fResidualR->At(FindMatchedPos(index));
1078   dZ = fResidualZ->At(FindMatchedPos(index));
1079   //printf("dR %f, dZ %f\n",dR, dZ);
1080 }
1081
1082 //__________________________________________________________
1083 Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t index)
1084 {
1085   //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1086   //Get the index of matched track for this cluster
1087   //It only works with ESDs, not AODs
1088   
1089   if(IsMatched(index))
1090     return fMatchedTrackIndex->At(FindMatchedPos(index));
1091   else 
1092     return -1; 
1093 }
1094
1095
1096 //__________________________________________________
1097 Bool_t AliEMCALRecoUtils::IsMatched(Int_t index)
1098 {
1099   //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1100   //Returns if cluster has a match
1101   if(FindMatchedPos(index) < 999) 
1102     return kTRUE;
1103   else
1104     return kFALSE;
1105 }
1106 //__________________________________________________________
1107 UInt_t AliEMCALRecoUtils::FindMatchedPos(Int_t index) const
1108 {
1109   //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1110   //Returns the position of the match in the fMatchedClusterIndex array
1111   Float_t tmpR = fCutR;
1112   UInt_t pos = 999;
1113   
1114   for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
1115   {
1116     if(fMatchedClusterIndex->At(i)==index && fResidualR->At(i)<tmpR)
1117     {
1118       pos=i;
1119       tmpR=fResidualR->At(i);
1120       AliDebug(3,Form("Matched cluster pos: %d, index: %d, dR: %2.4f, dZ: %2.4f.\n",i,fMatchedClusterIndex->At(i),fResidualR->At(i),fResidualZ->At(i)));
1121     }
1122   }
1123   return pos;
1124 }
1125
1126 //__________________________________________________________
1127 Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
1128 {
1129   // Given a esd track, return whether the track survive all the cuts
1130
1131   // The different quality parameter are first
1132   // retrieved from the track. then it is found out what cuts the
1133   // track did not survive and finally the cuts are imposed.
1134
1135   UInt_t status = esdTrack->GetStatus();
1136
1137   Int_t nClustersITS = esdTrack->GetITSclusters(0);
1138   Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
1139
1140   Float_t chi2PerClusterITS = -1;
1141   Float_t chi2PerClusterTPC = -1;
1142   if (nClustersITS!=0)
1143     chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1144   if (nClustersTPC!=0) 
1145     chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
1146
1147
1148   //DCA cuts
1149   Float_t MaxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1150   //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
1151   SetMaxDCAToVertexXY(MaxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
1152
1153
1154   Float_t b[2];
1155   Float_t bCov[3];
1156   esdTrack->GetImpactParameters(b,bCov);
1157   if (bCov[0]<=0 || bCov[2]<=0) {
1158     AliDebug(1, "Estimated b resolution lower or equal zero!");
1159     bCov[0]=0; bCov[2]=0;
1160   }
1161
1162   Float_t dcaToVertexXY = b[0];
1163   Float_t dcaToVertexZ = b[1];
1164   Float_t dcaToVertex = -1;
1165
1166   if (fCutDCAToVertex2D)
1167     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1168   else
1169     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1170     
1171   // cut the track?
1172   
1173   Bool_t cuts[kNCuts];
1174   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1175   
1176   // track quality cuts
1177   if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1178     cuts[0]=kTRUE;
1179   if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1180     cuts[1]=kTRUE;
1181   if (nClustersTPC<fCutMinNClusterTPC)
1182     cuts[2]=kTRUE;
1183   if (nClustersITS<fCutMinNClusterITS) 
1184     cuts[3]=kTRUE;
1185   if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC) 
1186     cuts[4]=kTRUE; 
1187   if (chi2PerClusterITS>fCutMaxChi2PerClusterITS) 
1188     cuts[5]=kTRUE;  
1189   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1190     cuts[6]=kTRUE;
1191   if (fCutDCAToVertex2D && dcaToVertex > 1)
1192     cuts[7] = kTRUE;
1193   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1194     cuts[8] = kTRUE;
1195   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1196     cuts[9] = kTRUE;
1197
1198   //Require at least one SPD point + anything else in ITS
1199   if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1200     cuts[10] = kTRUE;
1201
1202   Bool_t cut=kFALSE;
1203   for (Int_t i=0; i<kNCuts; i++) 
1204     if (cuts[i]) {cut = kTRUE;}
1205
1206     // cut the track
1207   if (cut) 
1208     return kFALSE;
1209   else 
1210     return kTRUE;
1211 }
1212 //__________________________________________________
1213 void AliEMCALRecoUtils::InitTrackCuts()
1214 {
1215   //Intilize the track cut criteria
1216   //By default these cuts are set according to AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1217   //Also you can customize the cuts using the setters
1218   
1219   //TPC
1220   SetMinNClustersTPC(70);
1221   SetMaxChi2PerClusterTPC(4);
1222   SetAcceptKinkDaughters(kFALSE);
1223   SetRequireTPCRefit(kTRUE);
1224   
1225   //ITS
1226   SetRequireITSRefit(kTRUE);
1227   SetMaxDCAToVertexZ(2);
1228   SetDCAToVertex2D(kFALSE);
1229   SetMaxChi2PerClusterITS(); //which by default sets the value to 1e10.
1230   SetMinNClustersITS();
1231 }
1232
1233 //___________________________________________________
1234 void AliEMCALRecoUtils::Print(const Option_t *) const 
1235 {
1236   // Print Parameters
1237   
1238   printf("AliEMCALRecoUtils Settings: \n");
1239   printf("Misalignment shifts\n");
1240   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, 
1241                                   fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1242                                   fMisalRotShift[i*3],  fMisalRotShift[i*3+1],  fMisalRotShift[i*3+2]   );
1243   printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1244   for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
1245   
1246   printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
1247
1248   printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
1249
1250   printf("Track cuts: \n");
1251   printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1252   printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1253   printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1254   printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1255   printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
1256
1257 }
1258
1259 //_____________________________________________________________________
1260 void AliEMCALRecoUtils::SetTimeDependentCorrections(Int_t runnumber){
1261   //Get EMCAL time dependent corrections from file and put them in the recalibration histograms
1262   //Do it only once and only if it is requested
1263   
1264   if(!fUseTimeCorrectionFactors) return;
1265   if(fTimeCorrectionFactorsSet)  return;
1266   
1267   printf("AliEMCALRecoUtils::GetTimeDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber);
1268  
1269   AliEMCALCalibTimeDepCorrection  *corr =  new AliEMCALCalibTimeDepCorrection();
1270   corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber));
1271   
1272   SwitchOnRecalibration();
1273   for(Int_t ism = 0; ism < 4; ism++){
1274     for(Int_t icol = 0; icol < 48; icol++){
1275       for(Int_t irow = 0; irow < 24; irow++){
1276         Float_t orgRecalFactor = GetEMCALChannelRecalibrationFactors(ism)->GetBinContent(icol,irow);
1277         Float_t newRecalFactor = orgRecalFactor*corr->GetCorrection(ism, icol,irow,0);
1278         GetEMCALChannelRecalibrationFactors(ism)->SetBinContent(icol,irow,newRecalFactor);
1279         //printf("ism %d, icol %d, irow %d, corrections : org %f, time dep %f, final %f (org*time %f)\n",ism, icol, irow, 
1280         //        orgRecalFactor, corr->GetCorrection(ism, icol,irow,0),
1281         //       (GetEMCALChannelRecalibrationFactors(ism))->GetBinContent(icol,irow),newRecalFactor);
1282       }
1283     }
1284   }
1285    fTimeCorrectionFactorsSet = kTRUE;
1286 }
1287