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