correct minor coding violation, load root libraries in case of par file compilation...
[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   GetMaxEnergyCell(geom, cells, cluster, absIdMax,  iSM, ieta, iphi);
235
236   AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n", 
237            absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
238         
239         if(absIdMax==-1) return kFALSE;
240         
241         //Check if the cell is close to the borders:
242         Bool_t okrow = kFALSE;
243         Bool_t okcol = kFALSE;
244   
245   if(iSM < 0 || iphi < 0 || ieta < 0 ) {
246     AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
247                   iSM,ieta,iphi));
248   }
249   
250   //Check rows/phi
251   if(iSM < 10){
252     if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE; 
253   }
254   else{
255     if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE; 
256   }
257   
258   //Check columns/eta
259   if(!fNoEMCALBorderAtEta0){
260     if(ieta  > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE; 
261   }
262   else{
263     if(iSM%2==0){
264       if(ieta >= fNCellsFromEMCALBorder)     okcol = kTRUE;     
265     }
266     else {
267       if(ieta <  48-fNCellsFromEMCALBorder)  okcol = kTRUE;     
268     }
269   }//eta 0 not checked
270     
271   AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d:  column? %d, row? %d\nq",
272            fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
273         
274         if (okcol && okrow) {
275     //printf("Accept\n");
276     return kTRUE;
277   }
278         else  {
279     //printf("Reject\n");
280     AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
281     return kFALSE;
282   }
283         
284 }       
285
286
287 //_________________________________________________________________________________________________________
288 Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, Int_t nCells){
289         // Check that in the cluster cells, there is no bad channel of those stored 
290         // in fEMCALBadChannelMap or fPHOSBadChannelMap
291         
292         if(!fRemoveBadChannels)  return kFALSE;
293         if(!fEMCALBadChannelMap) return kFALSE;
294         
295         Int_t icol = -1;
296         Int_t irow = -1;
297         Int_t imod = -1;
298         for(Int_t iCell = 0; iCell<nCells; iCell++){
299                 
300                 //Get the column and row
301     Int_t iTower = -1, iIphi = -1, iIeta = -1; 
302     geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta); 
303     if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
304     geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);                      
305     if(GetEMCALChannelStatus(imod, icol, irow)){
306       AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
307       return kTRUE;
308     }
309                 
310         }// cell cluster loop
311         
312         return kFALSE;
313         
314 }
315
316 //__________________________________________________
317 Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
318 // Correct cluster energy from non linearity functions
319   Float_t energy = cluster->E();
320   
321   switch (fNonLinearityFunction) {
322       
323     case kPi0MC:
324       //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
325       //Double_t par0 = 1.001;
326       //Double_t par1 = -0.01264;
327       //Double_t par2 = -0.03632;
328       //Double_t par3 = 0.1798;
329       //Double_t par4 = -0.522;
330        energy /= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
331                   ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
332                     exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
333       break;
334       
335     case kPi0GammaGamma:
336
337       //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
338       //Double_t par0 = 0.1457;
339       //Double_t par1 = -0.02024;
340       //Double_t par2 = 1.046;
341       energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
342       break;
343       
344     case kPi0GammaConversion:
345       
346       //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
347       //Double_t C = 0.139393/0.1349766;
348       //Double_t a = 0.0566186;
349       //Double_t b = 0.982133;
350       energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
351       
352       break;
353       
354     case kNoCorrection:
355       AliDebug(2,"No correction on the energy\n");
356       break;
357       
358   }
359   
360   return energy;
361
362 }
363
364 //__________________________________________________
365 Float_t  AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const 
366 {
367   //Calculate shower depth for a given cluster energy and particle type
368
369   // parameters 
370   Float_t x0    = 1.23;
371   Float_t ecr   = 8;
372   Float_t depth = 0;
373   
374   switch ( iParticle )
375   {
376     case kPhoton:
377       depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
378       break;
379       
380     case kElectron:
381       depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
382       break;
383       
384     case kHadron:
385       // hadron 
386       // boxes anc. here
387       if(gGeoManager){
388         gGeoManager->cd("ALIC_1/XEN1_1");
389         TGeoNode        *geoXEn1    = gGeoManager->GetCurrentNode();
390         TGeoNodeMatrix  *geoSM      = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
391         if(geoSM){
392           TGeoVolume      *geoSMVol   = geoSM->GetVolume(); 
393           TGeoShape       *geoSMShape = geoSMVol->GetShape();
394           TGeoBBox        *geoBox     = dynamic_cast<TGeoBBox *>(geoSMShape);
395           if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
396           else AliFatal("Null GEANT box");
397         }else AliFatal("NULL  GEANT node matrix");
398       }
399       else{//electron
400         depth = x0 * (TMath::Log(energy*1000 / ecr)  - 0.5); //Multiply energy by 1000 to transform to MeV
401       }
402         
403       break;
404       
405     default://photon
406       depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
407   }  
408   
409   return depth;
410   
411 }
412
413 //__________________________________________________
414 void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu, Int_t & absId,  Int_t& iSupMod, Int_t& ieta, Int_t& iphi)
415 {
416   //For a given CaloCluster gets the absId of the cell 
417   //with maximum energy deposit.
418   
419   Double_t eMax        = -1.;
420   Double_t eCell       = -1.;
421   Float_t  fraction    = 1.;
422   Float_t  recalFactor = 1.;
423   Int_t    cellAbsId   = -1 ;
424
425   Int_t iTower  = -1;
426   Int_t iIphi   = -1;
427   Int_t iIeta   = -1;
428         //printf("---Max?\n");
429   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
430     cellAbsId = clu->GetCellAbsId(iDig);
431     fraction  = clu->GetCellAmplitudeFraction(iDig);
432     //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
433     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
434     if(IsRecalibrationOn()) {
435       geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta); 
436       geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
437       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
438     }
439     eCell  = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
440     //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
441     if(eCell > eMax)  { 
442       eMax  = eCell; 
443       absId = cellAbsId;
444       //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
445     }
446   }// cell loop
447   
448   //Get from the absid the supermodule, tower and eta/phi numbers
449   geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta); 
450   //Gives SuperModule and Tower numbers
451   geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
452                                          iIphi, iIeta,iphi,ieta); 
453   //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
454   //printf("Max end---\n");
455   
456 }
457
458 //________________________________________________________________
459 void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
460         //Init EMCAL recalibration factors
461         AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
462         //In order to avoid rewriting the same histograms
463         Bool_t oldStatus = TH1::AddDirectoryStatus();
464         TH1::AddDirectory(kFALSE);
465   
466         fEMCALRecalibrationFactors = new TObjArray(12);
467         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));
468         //Init the histograms with 1
469         for (Int_t sm = 0; sm < 12; sm++) {
470                 for (Int_t i = 0; i < 48; i++) {
471                         for (Int_t j = 0; j < 24; j++) {
472                                 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
473                         }
474                 }
475         }
476         fEMCALRecalibrationFactors->SetOwner(kTRUE);
477         fEMCALRecalibrationFactors->Compress();
478         
479         //In order to avoid rewriting the same histograms
480         TH1::AddDirectory(oldStatus);           
481 }
482
483
484 //________________________________________________________________
485 void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
486         //Init EMCAL bad channels map
487         AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
488         //In order to avoid rewriting the same histograms
489         Bool_t oldStatus = TH1::AddDirectoryStatus();
490         TH1::AddDirectory(kFALSE);
491         
492         fEMCALBadChannelMap = new TObjArray(12);
493         //TH2F * hTemp = new  TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
494         for (int i = 0; i < 12; i++) {
495                 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
496         }
497         
498         //delete hTemp;
499         
500         fEMCALBadChannelMap->SetOwner(kTRUE);
501         fEMCALBadChannelMap->Compress();
502         
503         //In order to avoid rewriting the same histograms
504         TH1::AddDirectory(oldStatus);           
505 }
506
507 //________________________________________________________________
508 void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells){
509         // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
510         
511         //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
512         UShort_t * index    = cluster->GetCellsAbsId() ;
513         Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
514         Int_t ncells = cluster->GetNCells();
515         
516         //Initialize some used variables
517         Float_t energy = 0;
518         Int_t absId    = -1;
519   Int_t icol = -1, irow = -1, imod=1;
520         Float_t factor = 1, frac = 0;
521         
522         //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
523         for(Int_t icell = 0; icell < ncells; icell++){
524                 absId = index[icell];
525                 frac =  fraction[icell];
526                 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
527                 Int_t iTower = -1, iIphi = -1, iIeta = -1; 
528                 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta); 
529                 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
530                 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);                  
531                 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
532     AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
533              imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
534                 
535                 energy += cells->GetCellAmplitude(absId)*factor*frac;
536         }
537         
538         
539                 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
540         
541         cluster->SetE(energy);
542         
543 }
544
545
546 //__________________________________________________
547 void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
548 {
549   //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
550   
551   if     (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
552   else if(fPosAlgo==kPosTowerIndex)  RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
553   else   AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
554   
555 }  
556
557 //__________________________________________________
558 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
559 {
560   // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
561   // The algorithm is a copy of what is done in AliEMCALRecPoint
562   
563   Double_t eCell       = 0.;
564   Float_t  fraction    = 1.;
565   Float_t  recalFactor = 1.;
566   
567   Int_t    absId   = -1;
568   Int_t    iTower  = -1, iIphi  = -1, iIeta  = -1;
569   Int_t    iSupModMax = -1, iSM=-1, iphi   = -1, ieta   = -1;
570   Float_t  weight = 0.,  totalWeight=0.;
571   Float_t  newPos[3] = {0,0,0};
572   Double_t pLocal[3], pGlobal[3];
573   
574   Float_t  clEnergy = clu->E(); //Energy already recalibrated previously
575   GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi);
576   Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
577   
578   //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
579   
580   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
581     absId = clu->GetCellAbsId(iDig);
582     fraction  = clu->GetCellAmplitudeFraction(iDig);
583     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
584     geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta); 
585     geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);                       
586     
587     if(IsRecalibrationOn()) {
588       recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
589     }
590     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
591     
592     weight = GetCellWeight(eCell,clEnergy);
593     //printf("cell energy %f, weight %f\n",eCell,weight);
594     totalWeight += weight;
595     geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
596     //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
597     geom->GetGlobal(pLocal,pGlobal,iSupModMax);
598     //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
599
600     for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
601     
602   }// cell loop
603   
604   if(totalWeight>0){
605     for(int i=0; i<3; i++ )    newPos[i] /= totalWeight;
606   }
607     
608   //Float_t pos[]={0,0,0};
609   //clu->GetPosition(pos);
610   //printf("OldPos  : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
611   //printf("NewPos  : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
612   
613         if(iSupModMax > 1) {//sector 1
614           newPos[0] +=fMisalTransShift[3];//-=3.093; 
615           newPos[1] +=fMisalTransShift[4];//+=6.82;
616           newPos[2] +=fMisalTransShift[5];//+=1.635;
617     //printf("   +    : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
618
619         }
620         else {//sector 0
621           newPos[0] +=fMisalTransShift[0];//+=1.134;
622           newPos[1] +=fMisalTransShift[1];//+=8.2;
623           newPos[2] +=fMisalTransShift[2];//+=1.197;
624     //printf("   +    : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
625
626         }
627   //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
628
629   clu->SetPosition(newPos);
630   
631 }  
632
633 //__________________________________________________
634 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
635 {
636   // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
637   // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
638   
639   Double_t eCell       = 1.;
640   Float_t  fraction    = 1.;
641   Float_t  recalFactor = 1.;
642   
643   Int_t absId   = -1;
644   Int_t iTower  = -1;
645   Int_t iIphi   = -1, iIeta   = -1;
646         Int_t iSupMod = -1, iSupModMax = -1;
647   Int_t iphi = -1, ieta =-1;
648   
649   Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
650   GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi);
651   Float_t  depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
652
653   Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
654   Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
655   Int_t startingSM = -1;
656   
657   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
658     absId = clu->GetCellAbsId(iDig);
659     fraction  = clu->GetCellAmplitudeFraction(iDig);
660     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
661     geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta); 
662     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);                   
663     
664     if     (iDig==0)  startingSM = iSupMod;
665     else if(iSupMod != startingSM) areInSameSM = kFALSE;
666
667     eCell  = cells->GetCellAmplitude(absId);
668     
669     if(IsRecalibrationOn()) {
670       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
671     }
672     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
673     
674     weight = GetCellWeight(eCell,clEnergy);
675     if(weight < 0) weight = 0;
676     totalWeight += weight;
677     weightedCol += ieta*weight;
678     weightedRow += iphi*weight;
679     
680     //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
681     
682     }// cell loop
683     
684   Float_t xyzNew[]={0.,0.,0.};
685   if(areInSameSM == kTRUE) {
686     //printf("In Same SM\n");
687     weightedCol = weightedCol/totalWeight;
688     weightedRow = weightedRow/totalWeight;
689     geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew); 
690   }
691   else {
692     //printf("In Different SM\n");
693     geom->RecalculateTowerPosition(iphi,        ieta,        iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew); 
694   }
695   
696   clu->SetPosition(xyzNew);
697   
698 }
699
700 //____________________________________________________________________________
701 void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){           
702         
703   //re-evaluate identification parameters with bayesian
704
705         if ( cluster->GetM02() != 0)
706     fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
707   
708   Float_t pidlist[AliPID::kSPECIESN+1];
709         for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
710   
711   cluster->SetPID(pidlist);
712         
713 }
714
715 //____________________________________________________________________________
716 void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster)
717 {
718   // Calculates new center of gravity in the local EMCAL-module coordinates 
719   // and tranfers into global ALICE coordinates
720   // Calculates Dispersion and main axis
721   
722   Int_t nstat  = 0;
723   Float_t wtot = 0. ;
724   Double_t eCell       = 0.;
725   Float_t  fraction    = 1.;
726   Float_t  recalFactor = 1.;
727
728   Int_t iSupMod = -1;
729   Int_t iTower  = -1;
730   Int_t iIphi   = -1;
731   Int_t iIeta   = -1;
732   Int_t iphi    = -1;
733   Int_t ieta    = -1;
734   Double_t etai = -1.;
735   Double_t phii = -1.;
736   
737   Double_t w     = 0.;
738   Double_t d     = 0.;
739   Double_t dxx   = 0.;
740   Double_t dzz   = 0.;
741   Double_t dxz   = 0.;  
742   Double_t xmean = 0.;
743   Double_t zmean = 0.;
744     
745   //Loop on cells
746   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
747     
748     //Get from the absid the supermodule, tower and eta/phi numbers
749     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta); 
750     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);        
751     
752     //Get the cell energy, if recalibration is on, apply factors
753     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
754     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
755     if(IsRecalibrationOn()) {
756       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
757     }
758     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
759     
760     if(cluster->E() > 0 && eCell > 0){
761       
762       w  = GetCellWeight(eCell,cluster->E());
763       
764       etai=(Double_t)ieta;
765       phii=(Double_t)iphi;              
766       if(w > 0.0) {
767         wtot += w ;
768         nstat++;                        
769         //Shower shape
770         dxx  += w * etai * etai ;
771         xmean+= w * etai ;
772         dzz  += w * phii * phii ;
773         zmean+= w * phii ; 
774         dxz  += w * etai * phii ; 
775       }
776     }
777     else
778       AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
779   }//cell loop
780   
781   //Normalize to the weight     
782   if (wtot > 0) {
783     xmean /= wtot ;
784     zmean /= wtot ;
785   }
786   else
787     AliError(Form("Wrong weight %f\n", wtot));
788   
789   //Calculate dispersion        
790   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
791     
792     //Get from the absid the supermodule, tower and eta/phi numbers
793     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta); 
794     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
795     
796     //Get the cell energy, if recalibration is on, apply factors
797     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
798     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
799     if(IsRecalibrationOn()) {
800       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
801     }
802     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
803     
804     if(cluster->E() > 0 && eCell > 0){
805       
806       w  = GetCellWeight(eCell,cluster->E());
807       
808       etai=(Double_t)ieta;
809       phii=(Double_t)iphi;              
810       if(w > 0.0)  d +=  w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean)); 
811     }
812     else
813       AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
814   }// cell loop
815   
816   //Normalize to the weigth and set shower shape parameters
817   if (wtot > 0 && nstat > 1) {
818     d /= wtot ;
819     dxx /= wtot ;
820     dzz /= wtot ;
821     dxz /= wtot ;
822     dxx -= xmean * xmean ;
823     dzz -= zmean * zmean ;
824     dxz -= xmean * zmean ;
825     cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
826     cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
827   }
828   else{
829     d=0. ;
830     cluster->SetM20(0.) ;
831     cluster->SetM02(0.) ;
832   }     
833   
834   if (d>=0)
835     cluster->SetDispersion(TMath::Sqrt(d)) ;
836   else    
837     cluster->SetDispersion(0) ;
838   
839 }
840
841 //__________________________________________________
842 void AliEMCALRecoUtils::FindMatches(AliVEvent *event)
843 {
844   //This function should be called before the cluster loop
845   //Before call this function, please recalculate the cluster positions
846   //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
847   //Store matched cluster indexes and residuals
848   //It only works with ESDs, not AODs
849   
850   fMatchedClusterIndex->Reset();
851   fResidualZ->Reset();
852   fResidualR->Reset();
853   
854   fMatchedClusterIndex->Set(100);
855   fResidualZ->Set(100);
856   fResidualR->Set(100);
857   
858   Int_t    matched=0;
859   Float_t  clsPos[3];
860   Double_t trkPos[3];
861   for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
862   {
863     AliESDtrack *track = ((AliESDEvent*)event)->GetTrack(itr);
864     if(!track || !IsAccepted(track)) continue;
865     
866     Float_t dRMax = fCutR, dZMax = fCutZ;
867     Int_t index = -1;
868     AliEMCALTrack *emctrack = new AliEMCALTrack(*track);
869     for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
870     {
871       AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
872       if(!cluster->IsEMCAL()) continue;
873       cluster->GetPosition(clsPos); //Has been recalculated
874       if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) )  continue;
875       emctrack->GetXYZ(trkPos);
876       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) );
877       Float_t tmpZ = TMath::Abs(clsPos[2]-trkPos[2]);
878       
879       if(tmpR<dRMax)
880             {
881               dRMax=tmpR;
882               dZMax=tmpZ;
883               index=icl;
884             }
885       
886     }//cluser loop
887     
888     if(index>-1)
889     {
890       fMatchedClusterIndex->AddAt(index,matched);
891       fResidualZ->AddAt(dZMax,matched);
892       fResidualR->AddAt(dRMax,matched);
893       matched++;
894     }
895     delete emctrack;
896   }//track loop
897   fMatchedClusterIndex->Set(matched);
898   fResidualZ->Set(matched);
899   fResidualR->Set(matched);
900   
901   //printf("Number of matched pairs: %d\n",matched);
902 }
903
904 //__________________________________________________
905 void AliEMCALRecoUtils::GetMatchedResiduals(Int_t index, Float_t &dR, Float_t &dZ)
906 {
907   //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
908   //Get the residuals dR and dZ for this cluster
909   //It only works with ESDs, not AODs
910
911   if( FindMatchedPos(index)==-1 )
912   {
913     AliDebug(2,"No matched tracks found!\n");
914     dR=999.;
915     dZ=999.;
916     return;
917   }
918   dR = fResidualR->At(FindMatchedPos(index));
919   dZ = fResidualZ->At(FindMatchedPos(index));
920 }
921
922 //__________________________________________________
923 Bool_t AliEMCALRecoUtils::IsMatched(Int_t index)
924 {
925   //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
926   //Returns if cluster has a match
927   if(FindMatchedPos(index)>-1) return kTRUE;
928   else
929     return kFALSE;
930 }
931 //__________________________________________________
932 Int_t AliEMCALRecoUtils::FindMatchedPos(Int_t index) const
933 {
934   //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
935   //Returns the position of the match in the fMatchedClusterIndex array
936   Float_t tmpR = fCutR;
937   Int_t pos=-1;
938
939   for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
940     {
941       if(fMatchedClusterIndex->At(i)==index && fResidualR->At(i)<tmpR)
942         {
943           pos=i;
944           tmpR=fResidualR->At(i);
945         }
946       //printf("Matched cluster pos: %d, index: %d, dR: %2.4f, dZ: %2.4f.\n",i,fMatchedClusterIndex->At(i),fResidualR->At(i),fResidualZ->At(i));
947     }
948   return pos;
949 }
950
951 Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
952 {
953   // Given a esd track, return whether the track survive all the cuts
954
955   // The different quality parameter are first
956   // retrieved from the track. then it is found out what cuts the
957   // track did not survive and finally the cuts are imposed.
958
959   UInt_t status = esdTrack->GetStatus();
960
961   Int_t nClustersITS = esdTrack->GetITSclusters(0);
962   Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
963
964   Float_t chi2PerClusterITS = -1;
965   Float_t chi2PerClusterTPC = -1;
966   if (nClustersITS!=0)
967     chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
968   if (nClustersTPC!=0) 
969     chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
970       
971   Float_t b[2];
972   Float_t bCov[3];
973   esdTrack->GetImpactParameters(b,bCov);
974   if (bCov[0]<=0 || bCov[2]<=0) {
975     AliDebug(1, "Estimated b resolution lower or equal zero!");
976     bCov[0]=0; bCov[2]=0;
977   }
978
979   Float_t dcaToVertexXY = b[0];
980   Float_t dcaToVertexZ = b[1];
981   Float_t dcaToVertex = -1;
982
983   if (fCutDCAToVertex2D)
984     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
985   else
986     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
987     
988   // cut the track?
989   
990   Bool_t cuts[kNCuts];
991   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
992   
993   // track quality cuts
994   if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
995     cuts[0]=kTRUE;
996   if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
997     cuts[1]=kTRUE;
998   if (nClustersTPC<fCutMinNClusterTPC)
999     cuts[2]=kTRUE;
1000   if (nClustersITS<fCutMinNClusterITS) 
1001     cuts[3]=kTRUE;
1002   if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC) 
1003     cuts[4]=kTRUE; 
1004   if (chi2PerClusterITS>fCutMaxChi2PerClusterITS) 
1005     cuts[5]=kTRUE;  
1006   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1007     cuts[6]=kTRUE;
1008   if (fCutDCAToVertex2D && dcaToVertex > 1)
1009     cuts[7] = kTRUE;
1010   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1011     cuts[8] = kTRUE;
1012   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1013     cuts[9] = kTRUE;
1014
1015   Bool_t cut=kFALSE;
1016   for (Int_t i=0; i<kNCuts; i++) 
1017     if (cuts[i]) {cut = kTRUE;}
1018
1019     // cut the track
1020   if (cut) 
1021     return kFALSE;
1022   else 
1023     return kTRUE;
1024 }
1025 //__________________________________________________
1026 void AliEMCALRecoUtils::InitTrackCuts()
1027 {
1028   //Intilize the track cut criteria
1029   //By default these cuts are set according to AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()
1030   //Also you can customize the cuts using the setters
1031
1032   SetMinNClustersTPC(50);
1033   SetMaxChi2PerClusterTPC(4);
1034   SetAcceptKinkDaughters(kFALSE);
1035   SetMaxDCAToVertexZ(3.2);
1036   SetMaxDCAToVertexXY(2.4);
1037   SetDCAToVertex2D(kTRUE);
1038 }
1039
1040 //__________________________________________________
1041 void AliEMCALRecoUtils::Print(const Option_t *) const 
1042 {
1043   // Print Parameters
1044   
1045   printf("AliEMCALRecoUtils Settings: \n");
1046   printf("Misalignment shifts\n");
1047   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, 
1048                                   fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1049                                   fMisalRotShift[i*3],  fMisalRotShift[i*3+1],  fMisalRotShift[i*3+2]   );
1050   printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1051   for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
1052   
1053   printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
1054
1055   printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
1056
1057   printf("Track cuts: \n");
1058   printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1059   printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1060   printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1061   printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1062   printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
1063
1064     
1065 }