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