Methods for shower shape parameters and PID recalculations added
[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 "AliLog.h"
43 #include "AliEMCALPIDUtils.h"
44 #include "AliPID.h"
45
46 ClassImp(AliEMCALRecoUtils)
47   
48 //______________________________________________
49 AliEMCALRecoUtils::AliEMCALRecoUtils():
50   fNonLinearityFunction (kNoCorrection), fParticleType(kPhoton),
51   fPosAlgo(kUnchanged),fW0(4.),
52   fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
53   fRemoveBadChannels(kFALSE),fEMCALBadChannelMap(),
54   fNCellsFromEMCALBorder(0),fNoEMCALBorderAtEta0(kTRUE), fPIDUtils()
55 {
56 //
57   // Constructor.
58   // Initialize all constant values which have to be used
59   // during Reco algorithm execution
60   //
61   
62   for(Int_t i = 0; i < 15 ; i++) {
63       fMisalTransShift[i] = 0.; 
64       fMisalRotShift[i] = 0.; 
65   }
66   for(Int_t i = 0; i < 6  ; i++) fNonLinearityParams[i] = 0.; 
67   //For kPi0GammaGamma case, but default is no correction
68   fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
69   fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
70   fNonLinearityParams[2] = 1.046;
71   
72   fPIDUtils = new AliEMCALPIDUtils();
73
74 }
75
76 //______________________________________________________________________
77 AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco) 
78 : TNamed(reco), fNonLinearityFunction(reco.fNonLinearityFunction), 
79   fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo), fW0(reco.fW0), 
80   fRecalibration(reco.fRecalibration),fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
81   fRemoveBadChannels(reco.fRemoveBadChannels),fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
82   fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
83   fPIDUtils(reco.fPIDUtils)
84
85 {
86   //Copy ctor
87   
88   for(Int_t i = 0; i < 15 ; i++) {
89       fMisalRotShift[i] = reco.fMisalRotShift[i]; 
90       fMisalTransShift[i] = reco.fMisalTransShift[i]; 
91   } 
92   for(Int_t i = 0; i < 6  ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i]; 
93 }
94
95
96 //______________________________________________________________________
97 AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco) 
98 {
99   //Assignment operator
100   
101   if(this == &reco)return *this;
102   ((TNamed *)this)->operator=(reco);
103
104   fNonLinearityFunction  = reco.fNonLinearityFunction;
105   fParticleType          = reco.fParticleType;
106   fPosAlgo               = reco.fPosAlgo; 
107   fW0                    = reco.fW0;
108   fRecalibration         = reco.fRecalibration;
109   fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
110   fRemoveBadChannels     = reco.fRemoveBadChannels;
111   fEMCALBadChannelMap    = reco.fEMCALBadChannelMap;
112   fNCellsFromEMCALBorder = reco.fNCellsFromEMCALBorder;
113   fNoEMCALBorderAtEta0   = reco.fNoEMCALBorderAtEta0;
114   fPIDUtils              = reco.fPIDUtils;
115
116   for(Int_t i = 0; i < 15 ; i++) {fMisalTransShift[i] = reco.fMisalTransShift[i]; fMisalRotShift[i] = reco.fMisalRotShift[i];}
117   for(Int_t i = 0; i < 6  ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i]; 
118   
119   return *this;
120 }
121
122
123 //__________________________________________________
124 AliEMCALRecoUtils::~AliEMCALRecoUtils()
125 {
126   //Destructor.
127         
128         if(fEMCALRecalibrationFactors) { 
129                 fEMCALRecalibrationFactors->Clear();
130                 delete  fEMCALRecalibrationFactors;
131         }       
132   
133   if(fEMCALBadChannelMap) { 
134                 fEMCALBadChannelMap->Clear();
135                 delete  fEMCALBadChannelMap;
136         }
137   
138 }
139
140 //_______________________________________________________________
141 Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells) 
142 {
143         // Given the list of AbsId of the cluster, get the maximum cell and 
144         // check if there are fNCellsFromBorder from the calorimeter border
145         
146   //If the distance to the border is 0 or negative just exit accept all clusters
147         if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
148   
149   Int_t absIdMax        = -1, iSM =-1, ieta = -1, iphi = -1;  
150   GetMaxEnergyCell(geom, cells, cluster, absIdMax,  iSM, ieta, iphi);
151
152   AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n", 
153            absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
154         
155         if(absIdMax==-1) return kFALSE;
156         
157         //Check if the cell is close to the borders:
158         Bool_t okrow = kFALSE;
159         Bool_t okcol = kFALSE;
160   
161   if(iSM < 0 || iphi < 0 || ieta < 0 ) {
162     AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
163                   iSM,ieta,iphi));
164   }
165   
166   //Check rows/phi
167   if(iSM < 10){
168     if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE; 
169   }
170   else{
171     if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE; 
172   }
173   
174   //Check columns/eta
175   if(!fNoEMCALBorderAtEta0){
176     if(ieta  > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE; 
177   }
178   else{
179     if(iSM%2==0){
180       if(ieta >= fNCellsFromEMCALBorder)     okcol = kTRUE;     
181     }
182     else {
183       if(ieta <  48-fNCellsFromEMCALBorder)  okcol = kTRUE;     
184     }
185   }//eta 0 not checked
186     
187   AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d:  column? %d, row? %d\nq",
188            fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
189         
190         if (okcol && okrow) {
191     //printf("Accept\n");
192     return kTRUE;
193   }
194         else  {
195     //printf("Reject\n");
196     AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
197     return kFALSE;
198   }
199         
200 }       
201
202
203 //_________________________________________________________________________________________________________
204 Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, Int_t nCells){
205         // Check that in the cluster cells, there is no bad channel of those stored 
206         // in fEMCALBadChannelMap or fPHOSBadChannelMap
207         
208         if(!fRemoveBadChannels)  return kFALSE;
209         if(!fEMCALBadChannelMap) return kFALSE;
210         
211         Int_t icol = -1;
212         Int_t irow = -1;
213         Int_t imod = -1;
214         for(Int_t iCell = 0; iCell<nCells; iCell++){
215                 
216                 //Get the column and row
217     Int_t iTower = -1, iIphi = -1, iIeta = -1; 
218     geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta); 
219     if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
220     geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);                      
221     if(GetEMCALChannelStatus(imod, icol, irow)){
222       AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
223       return kTRUE;
224     }
225                 
226         }// cell cluster loop
227         
228         return kFALSE;
229         
230 }
231
232 //__________________________________________________
233 Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
234 // Correct cluster energy from non linearity functions
235   Float_t energy = cluster->E();
236   
237   switch (fNonLinearityFunction) {
238       
239     case kPi0MC:
240       //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
241       //Double_t par0 = 1.001;
242       //Double_t par1 = -0.01264;
243       //Double_t par2 = -0.03632;
244       //Double_t par3 = 0.1798;
245       //Double_t par4 = -0.522;
246        energy /= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
247                   ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
248                     exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
249       break;
250       
251     case kPi0GammaGamma:
252
253       //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
254       //Double_t par0 = 0.1457;
255       //Double_t par1 = -0.02024;
256       //Double_t par2 = 1.046;
257       energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
258       break;
259       
260     case kPi0GammaConversion:
261       
262       //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
263       //Double_t C = 0.139393/0.1349766;
264       //Double_t a = 0.0566186;
265       //Double_t b = 0.982133;
266       energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
267       
268       break;
269       
270     case kNoCorrection:
271       AliDebug(2,"No correction on the energy\n");
272       break;
273       
274   }
275   
276   return energy;
277
278 }
279
280 //__________________________________________________
281 Float_t  AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const 
282 {
283   //Calculate shower depth for a given cluster energy and particle type
284
285   // parameters 
286   Float_t x0    = 1.23;
287   Float_t ecr   = 8;
288   Float_t depth = 0;
289   
290   switch ( iParticle )
291   {
292     case kPhoton:
293       depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
294       break;
295       
296     case kElectron:
297       depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
298       break;
299       
300     case kHadron:
301       // hadron 
302       // boxes anc. here
303       if(gGeoManager){
304         gGeoManager->cd("ALIC_1/XEN1_1");
305         TGeoNode        *geoXEn1    = gGeoManager->GetCurrentNode();
306         TGeoNodeMatrix  *geoSM      = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
307         if(geoSM){
308           TGeoVolume      *geoSMVol   = geoSM->GetVolume(); 
309           TGeoShape       *geoSMShape = geoSMVol->GetShape();
310           TGeoBBox        *geoBox     = dynamic_cast<TGeoBBox *>(geoSMShape);
311           if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
312           else AliFatal("Null GEANT box");
313         }else AliFatal("NULL  GEANT node matrix");
314       }
315       else{//electron
316         depth = x0 * (TMath::Log(energy*1000 / ecr)  - 0.5); //Multiply energy by 1000 to transform to MeV
317       }
318         
319       break;
320       
321     default://photon
322       depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
323   }  
324   
325   return depth;
326   
327 }
328
329 //__________________________________________________
330 void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu, Int_t & absId,  Int_t& iSupMod, Int_t& ieta, Int_t& iphi)
331 {
332   //For a given CaloCluster gets the absId of the cell 
333   //with maximum energy deposit.
334   
335   Double_t eMax        = -1.;
336   Double_t eCell       = -1.;
337   Float_t  fraction    = 1.;
338   Float_t  recalFactor = 1.;
339   Int_t    cellAbsId   = -1 ;
340
341   Int_t iTower  = -1;
342   Int_t iIphi   = -1;
343   Int_t iIeta   = -1;
344         //printf("---Max?\n");
345   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
346     cellAbsId = clu->GetCellAbsId(iDig);
347     fraction  = clu->GetCellAmplitudeFraction(iDig);
348     //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
349     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
350     if(IsRecalibrationOn()) {
351       geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta); 
352       geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
353       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
354     }
355     eCell  = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
356     //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
357     if(eCell > eMax)  { 
358       eMax  = eCell; 
359       absId = cellAbsId;
360       //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
361     }
362   }// cell loop
363   
364   //Get from the absid the supermodule, tower and eta/phi numbers
365   geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta); 
366   //Gives SuperModule and Tower numbers
367   geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
368                                          iIphi, iIeta,iphi,ieta); 
369   //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
370   //printf("Max end---\n");
371   
372 }
373
374 //________________________________________________________________
375 void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
376         //Init EMCAL recalibration factors
377         AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
378         //In order to avoid rewriting the same histograms
379         Bool_t oldStatus = TH1::AddDirectoryStatus();
380         TH1::AddDirectory(kFALSE);
381   
382         fEMCALRecalibrationFactors = new TObjArray(12);
383         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));
384         //Init the histograms with 1
385         for (Int_t sm = 0; sm < 12; sm++) {
386                 for (Int_t i = 0; i < 48; i++) {
387                         for (Int_t j = 0; j < 24; j++) {
388                                 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
389                         }
390                 }
391         }
392         fEMCALRecalibrationFactors->SetOwner(kTRUE);
393         fEMCALRecalibrationFactors->Compress();
394         
395         //In order to avoid rewriting the same histograms
396         TH1::AddDirectory(oldStatus);           
397 }
398
399
400 //________________________________________________________________
401 void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
402         //Init EMCAL bad channels map
403         AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
404         //In order to avoid rewriting the same histograms
405         Bool_t oldStatus = TH1::AddDirectoryStatus();
406         TH1::AddDirectory(kFALSE);
407         
408         fEMCALBadChannelMap = new TObjArray(12);
409         //TH2F * hTemp = new  TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
410         for (int i = 0; i < 12; i++) {
411                 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
412         }
413         
414         //delete hTemp;
415         
416         fEMCALBadChannelMap->SetOwner(kTRUE);
417         fEMCALBadChannelMap->Compress();
418         
419         //In order to avoid rewriting the same histograms
420         TH1::AddDirectory(oldStatus);           
421 }
422
423 //________________________________________________________________
424 void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells){
425         // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
426         
427         //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
428         UShort_t * index    = cluster->GetCellsAbsId() ;
429         Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
430         Int_t ncells = cluster->GetNCells();
431         
432         //Initialize some used variables
433         Float_t energy = 0;
434         Int_t absId    = -1;
435   Int_t icol = -1, irow = -1, imod=1;
436         Float_t factor = 1, frac = 0;
437         
438         //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
439         for(Int_t icell = 0; icell < ncells; icell++){
440                 absId = index[icell];
441                 frac =  fraction[icell];
442                 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
443                 Int_t iTower = -1, iIphi = -1, iIeta = -1; 
444                 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta); 
445                 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
446                 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);                  
447                 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
448     AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
449              imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
450                 
451                 energy += cells->GetCellAmplitude(absId)*factor*frac;
452         }
453         
454         
455                 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
456         
457         cluster->SetE(energy);
458         
459 }
460
461
462 //__________________________________________________
463 void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
464 {
465   //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
466   
467   if     (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
468   else if(fPosAlgo==kPosTowerIndex)  RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
469   else   AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
470   
471 }  
472
473 //__________________________________________________
474 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
475 {
476   // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
477   // The algorithm is a copy of what is done in AliEMCALRecPoint
478   
479   Double_t eCell       = 0.;
480   Float_t  fraction    = 1.;
481   Float_t  recalFactor = 1.;
482   
483   Int_t    absId   = -1;
484   Int_t    iTower  = -1, iIphi  = -1, iIeta  = -1;
485   Int_t    iSupModMax = -1, iSM=-1, iphi   = -1, ieta   = -1;
486   Float_t  weight = 0.,  totalWeight=0.;
487   Float_t  newPos[3] = {0,0,0};
488   Double_t pLocal[3], pGlobal[3];
489   
490   Float_t  clEnergy = clu->E(); //Energy already recalibrated previously
491   GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi);
492   Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
493   
494   //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
495   
496   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
497     absId = clu->GetCellAbsId(iDig);
498     fraction  = clu->GetCellAmplitudeFraction(iDig);
499     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
500     geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta); 
501     geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);                       
502     
503     if(IsRecalibrationOn()) {
504       recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
505     }
506     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
507     
508     weight = GetCellWeight(eCell,clEnergy);
509     //printf("cell energy %f, weight %f\n",eCell,weight);
510     totalWeight += weight;
511     geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
512     //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
513     geom->GetGlobal(pLocal,pGlobal,iSupModMax);
514     //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
515
516     for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
517     
518   }// cell loop
519   
520   if(totalWeight>0){
521     for(int i=0; i<3; i++ )    newPos[i] /= totalWeight;
522   }
523     
524   //Float_t pos[]={0,0,0};
525   //clu->GetPosition(pos);
526   //printf("OldPos  : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
527   //printf("NewPos  : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
528   
529         if(iSupModMax > 1) {//sector 1
530           newPos[0] +=fMisalTransShift[3];//-=3.093; 
531           newPos[1] +=fMisalTransShift[4];//+=6.82;
532           newPos[2] +=fMisalTransShift[5];//+=1.635;
533     //printf("   +    : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
534
535         }
536         else {//sector 0
537           newPos[0] +=fMisalTransShift[0];//+=1.134;
538           newPos[1] +=fMisalTransShift[1];//+=8.2;
539           newPos[2] +=fMisalTransShift[2];//+=1.197;
540     //printf("   +    : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
541
542         }
543   //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
544
545   clu->SetPosition(newPos);
546   
547 }  
548
549 //__________________________________________________
550 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
551 {
552   // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
553   // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
554   
555   Double_t eCell       = 1.;
556   Float_t  fraction    = 1.;
557   Float_t  recalFactor = 1.;
558   
559   Int_t absId   = -1;
560   Int_t iTower  = -1;
561   Int_t iIphi   = -1, iIeta   = -1;
562         Int_t iSupMod = -1, iSupModMax = -1;
563   Int_t iphi = -1, ieta =-1;
564   
565   Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
566   GetMaxEnergyCell(geom, cells, clu, absId,  iSupModMax, ieta, iphi);
567   Float_t  depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
568
569   Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
570   Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
571   Int_t startingSM = -1;
572   
573   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
574     absId = clu->GetCellAbsId(iDig);
575     fraction  = clu->GetCellAmplitudeFraction(iDig);
576     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
577     geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta); 
578     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);                   
579     
580     if     (iDig==0)  startingSM = iSupMod;
581     else if(iSupMod != startingSM) areInSameSM = kFALSE;
582
583     eCell  = cells->GetCellAmplitude(absId);
584     
585     if(IsRecalibrationOn()) {
586       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
587     }
588     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
589     
590     weight = GetCellWeight(eCell,clEnergy);
591     if(weight < 0) weight = 0;
592     totalWeight += weight;
593     weightedCol += ieta*weight;
594     weightedRow += iphi*weight;
595     
596     //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
597     
598     }// cell loop
599     
600   Float_t xyzNew[]={0.,0.,0.};
601   if(areInSameSM == kTRUE) {
602     //printf("In Same SM\n");
603     weightedCol = weightedCol/totalWeight;
604     weightedRow = weightedRow/totalWeight;
605     geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew); 
606   }
607   else {
608     //printf("In Different SM\n");
609     geom->RecalculateTowerPosition(iphi,        ieta,        iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew); 
610   }
611   
612   clu->SetPosition(xyzNew);
613   
614 }
615
616 //____________________________________________________________________________
617 void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){           
618         
619   //re-evaluate identification parameters with bayesian
620
621         if ( cluster->GetM02() != 0)
622     fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
623   
624   Float_t pidlist[AliPID::kSPECIESN+1];
625         for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
626   
627   cluster->SetPID(pidlist);
628         
629 }
630
631 //____________________________________________________________________________
632 void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster)
633 {
634   // Calculates new center of gravity in the local EMCAL-module coordinates 
635   // and tranfers into global ALICE coordinates
636   // Calculates Dispersion and main axis
637   
638   Int_t nstat  = 0;
639   Float_t wtot = 0. ;
640   Double_t eCell       = 0.;
641   Float_t  fraction    = 1.;
642   Float_t  recalFactor = 1.;
643
644   Int_t iSupMod = -1;
645   Int_t iTower  = -1;
646   Int_t iIphi   = -1;
647   Int_t iIeta   = -1;
648   Int_t iphi    = -1;
649   Int_t ieta    = -1;
650   Double_t etai = -1.;
651   Double_t phii = -1.;
652   
653   Double_t w     = 0.;
654   Double_t d     = 0.;
655   Double_t dxx   = 0.;
656   Double_t dzz   = 0.;
657   Double_t dxz   = 0.;  
658   Double_t xmean = 0.;
659   Double_t zmean = 0.;
660     
661   //Loop on cells
662   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
663     
664     //Get from the absid the supermodule, tower and eta/phi numbers
665     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta); 
666     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);        
667     
668     //Get the cell energy, if recalibration is on, apply factors
669     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
670     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
671     if(IsRecalibrationOn()) {
672       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
673     }
674     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
675     
676     if(cluster->E() > 0 && eCell > 0){
677       
678       w  = GetCellWeight(eCell,cluster->E());
679       
680       etai=(Double_t)ieta;
681       phii=(Double_t)iphi;              
682       if(w > 0.0) {
683         wtot += w ;
684         nstat++;                        
685         //Shower shape
686         dxx  += w * etai * etai ;
687         xmean+= w * etai ;
688         dzz  += w * phii * phii ;
689         zmean+= w * phii ; 
690         dxz  += w * etai * phii ; 
691       }
692     }
693     else
694       AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
695   }//cell loop
696   
697   //Normalize to the weight     
698   if (wtot > 0) {
699     xmean /= wtot ;
700     zmean /= wtot ;
701   }
702   else
703     AliError(Form("Wrong weight %f\n", wtot));
704   
705   //Calculate dispersion        
706   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
707     
708     //Get from the absid the supermodule, tower and eta/phi numbers
709     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta); 
710     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
711     
712     //Get the cell energy, if recalibration is on, apply factors
713     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
714     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
715     if(IsRecalibrationOn()) {
716       recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
717     }
718     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
719     
720     if(cluster->E() > 0 && eCell > 0){
721       
722       w  = GetCellWeight(eCell,cluster->E());
723       
724       etai=(Double_t)ieta;
725       phii=(Double_t)iphi;              
726       if(w > 0.0)  d +=  w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean)); 
727     }
728     else
729       AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
730   }// cell loop
731   
732   //Normalize to the weigth and set shower shape parameters
733   if (wtot > 0 && nstat > 1) {
734     d /= wtot ;
735     dxx /= wtot ;
736     dzz /= wtot ;
737     dxz /= wtot ;
738     dxx -= xmean * xmean ;
739     dzz -= zmean * zmean ;
740     dxz -= xmean * zmean ;
741     cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
742     cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
743   }
744   else{
745     d=0. ;
746     cluster->SetM20(0.) ;
747     cluster->SetM02(0.) ;
748   }     
749   
750   if (d>=0)
751     cluster->SetDispersion(TMath::Sqrt(d)) ;
752   else    
753     cluster->SetDispersion(0) ;
754   
755 }
756
757
758 //__________________________________________________
759 void AliEMCALRecoUtils::Print(const Option_t *) const 
760 {
761   // Print Parameters
762   
763   printf("AliEMCALRecoUtils Settings: \n");
764   printf("Misalignment shifts\n");
765   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, 
766                                   fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
767                                   fMisalRotShift[i*3],  fMisalRotShift[i*3+1],  fMisalRotShift[i*3+2]   );
768   printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
769   for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
770   
771   printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
772     
773 }