]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EMCAL/AliEMCALRecoUtils.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Move option to reclusterizer with calibrated amplitudes from NxN to mother Clusterize...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALRecoUtils.cxx
CommitLineData
d9b3567c 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
094786cc 33#include <TGeoManager.h>
34#include <TGeoMatrix.h>
35#include <TGeoBBox.h>
d9b3567c 36
37// STEER includes
38#include "AliEMCALRecoUtils.h"
094786cc 39#include "AliEMCALGeometry.h"
d9b3567c 40#include "AliVCluster.h"
41#include "AliVCaloCells.h"
42#include "AliLog.h"
83bfd77a 43#include "AliEMCALPIDUtils.h"
44#include "AliPID.h"
d9b3567c 45
46ClassImp(AliEMCALRecoUtils)
47
48//______________________________________________
49AliEMCALRecoUtils::AliEMCALRecoUtils():
fd6df01c 50 fNonLinearityFunction (kNoCorrection), fParticleType(kPhoton),
51 fPosAlgo(kUnchanged),fW0(4.),
52 fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
53 fRemoveBadChannels(kFALSE),fEMCALBadChannelMap(),
83bfd77a 54 fNCellsFromEMCALBorder(0),fNoEMCALBorderAtEta0(kTRUE), fPIDUtils()
d9b3567c 55{
56//
57 // Constructor.
58 // Initialize all constant values which have to be used
59 // during Reco algorithm execution
60 //
61
fd6df01c 62 for(Int_t i = 0; i < 15 ; i++) {
63 fMisalTransShift[i] = 0.;
64 fMisalRotShift[i] = 0.;
65 }
d9b3567c 66 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = 0.;
fd6df01c 67 //For kPi0GammaGamma case, but default is no correction
d9b3567c 68 fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
69 fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
70 fNonLinearityParams[2] = 1.046;
71
83bfd77a 72 fPIDUtils = new AliEMCALPIDUtils();
73
d9b3567c 74}
75
76//______________________________________________________________________
77AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
094786cc 78: TNamed(reco), fNonLinearityFunction(reco.fNonLinearityFunction),
fd6df01c 79 fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo), fW0(reco.fW0),
80 fRecalibration(reco.fRecalibration),fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
81 fRemoveBadChannels(reco.fRemoveBadChannels),fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
83bfd77a 82 fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
83 fPIDUtils(reco.fPIDUtils)
fd6df01c 84
d9b3567c 85{
86 //Copy ctor
87
fd6df01c 88 for(Int_t i = 0; i < 15 ; i++) {
89 fMisalRotShift[i] = reco.fMisalRotShift[i];
90 fMisalTransShift[i] = reco.fMisalTransShift[i];
91 }
d9b3567c 92 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
93}
94
95
96//______________________________________________________________________
97AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco)
98{
99 //Assignment operator
100
101 if(this == &reco)return *this;
102 ((TNamed *)this)->operator=(reco);
103
fd6df01c 104 fNonLinearityFunction = reco.fNonLinearityFunction;
105 fParticleType = reco.fParticleType;
106 fPosAlgo = reco.fPosAlgo;
107 fW0 = reco.fW0;
108 fRecalibration = reco.fRecalibration;
094786cc 109 fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
fd6df01c 110 fRemoveBadChannels = reco.fRemoveBadChannels;
111 fEMCALBadChannelMap = reco.fEMCALBadChannelMap;
112 fNCellsFromEMCALBorder = reco.fNCellsFromEMCALBorder;
113 fNoEMCALBorderAtEta0 = reco.fNoEMCALBorderAtEta0;
83bfd77a 114 fPIDUtils = reco.fPIDUtils;
115
2a71e873 116 for(Int_t i = 0; i < 15 ; i++) {fMisalTransShift[i] = reco.fMisalTransShift[i]; fMisalRotShift[i] = reco.fMisalRotShift[i];}
d9b3567c 117 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
118
119 return *this;
120}
121
122
094786cc 123//__________________________________________________
124AliEMCALRecoUtils::~AliEMCALRecoUtils()
125{
126 //Destructor.
127
128 if(fEMCALRecalibrationFactors) {
129 fEMCALRecalibrationFactors->Clear();
130 delete fEMCALRecalibrationFactors;
131 }
fd6df01c 132
133 if(fEMCALBadChannelMap) {
134 fEMCALBadChannelMap->Clear();
135 delete fEMCALBadChannelMap;
136 }
137
094786cc 138}
139
fd6df01c 140//_______________________________________________________________
141Bool_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
83bfd77a 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));
fd6df01c 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
83bfd77a 187 AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d: column? %d, row? %d\nq",
fd6df01c 188 fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
189
83bfd77a 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 }
fd6df01c 199
200}
201
202
203//_________________________________________________________________________________________________________
204Bool_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);
83bfd77a 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 }
fd6df01c 225
226 }// cell cluster loop
227
228 return kFALSE;
229
230}
094786cc 231
d9b3567c 232//__________________________________________________
233Float_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//__________________________________________________
094786cc 281Float_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:
fd6df01c 293 depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
094786cc 294 break;
295
296 case kElectron:
fd6df01c 297 depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
094786cc 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));
fd6df01c 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");
094786cc 314 }
315 else{//electron
fd6df01c 316 depth = x0 * (TMath::Log(energy*1000 / ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
094786cc 317 }
318
319 break;
320
321 default://photon
fd6df01c 322 depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
094786cc 323 }
324
325 return depth;
326
327}
328
329//__________________________________________________
330void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu, Int_t & absId, Int_t& iSupMod, Int_t& ieta, Int_t& iphi)
d9b3567c 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.;
094786cc 337 Float_t fraction = 1.;
338 Float_t recalFactor = 1.;
d9b3567c 339 Int_t cellAbsId = -1 ;
094786cc 340
d9b3567c 341 Int_t iTower = -1;
342 Int_t iIphi = -1;
343 Int_t iIeta = -1;
83bfd77a 344 //printf("---Max?\n");
d9b3567c 345 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
094786cc 346 cellAbsId = clu->GetCellAbsId(iDig);
347 fraction = clu->GetCellAmplitudeFraction(iDig);
83bfd77a 348 //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
094786cc 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;
83bfd77a 356 //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
094786cc 357 if(eCell > eMax) {
d9b3567c 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,
83bfd77a 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");
d9b3567c 371
372}
373
094786cc 374//________________________________________________________________
375void 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
fd6df01c 400//________________________________________________________________
401void 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
094786cc 423//________________________________________________________________
424void 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
d9b3567c 462//__________________________________________________
094786cc 463void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
d9b3567c 464{
465 //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
466
094786cc 467 if (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
468 else if(fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
fd6df01c 469 else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
094786cc 470
471}
472
473//__________________________________________________
474void 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
094786cc 491 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi);
492 Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
493
83bfd77a 494 //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
495
094786cc 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);
83bfd77a 509 //printf("cell energy %f, weight %f\n",eCell,weight);
094786cc 510 totalWeight += weight;
511 geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
83bfd77a 512 //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
094786cc 513 geom->GetGlobal(pLocal,pGlobal,iSupModMax);
83bfd77a 514 //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
515
094786cc 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
094786cc 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]);
83bfd77a 527 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
094786cc 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;
83bfd77a 533 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
534
094786cc 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;
83bfd77a 540 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
541
094786cc 542 }
83bfd77a 543 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
544
094786cc 545 clu->SetPosition(newPos);
546
094786cc 547}
548
549//__________________________________________________
550void 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;
d9b3567c 560 Int_t iTower = -1;
094786cc 561 Int_t iIphi = -1, iIeta = -1;
562 Int_t iSupMod = -1, iSupModMax = -1;
d9b3567c 563 Int_t iphi = -1, ieta =-1;
564
565 Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
094786cc 566 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi);
567 Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
568
d9b3567c 569 Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
094786cc 570 Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
571 Int_t startingSM = -1;
d9b3567c 572
573 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
094786cc 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);
d9b3567c 578 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
094786cc 579
d9b3567c 580 if (iDig==0) startingSM = iSupMod;
581 else if(iSupMod != startingSM) areInSameSM = kFALSE;
094786cc 582
583 eCell = cells->GetCellAmplitude(absId);
d9b3567c 584
094786cc 585 if(IsRecalibrationOn()) {
586 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
587 }
588 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
d9b3567c 589
094786cc 590 weight = GetCellWeight(eCell,clEnergy);
d9b3567c 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
094786cc 598 }// cell loop
599
d9b3567c 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;
094786cc 605 geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
d9b3567c 606 }
607 else {
608 //printf("In Different SM\n");
094786cc 609 geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
d9b3567c 610 }
d9b3567c 611
094786cc 612 clu->SetPosition(xyzNew);
d9b3567c 613
614}
615
83bfd77a 616//____________________________________________________________________________
617void 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//____________________________________________________________________________
632void 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
d9b3567c 758//__________________________________________________
759void AliEMCALRecoUtils::Print(const Option_t *) const
760{
761 // Print Parameters
762
763 printf("AliEMCALRecoUtils Settings: \n");
764 printf("Misalignment shifts\n");
2a71e873 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] );
d9b3567c 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]);
094786cc 770
771 printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
d9b3567c 772
773}