HMPID friends (Giacomo)
[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():
094786cc 48 fNonLinearityFunction (kPi0GammaGamma), fParticleType(kPhoton),
49 fPosAlgo(kPosTowerIndex),fW0(4.),fRecalibration(kFALSE), fEMCALRecalibrationFactors()
d9b3567c 50{
51//
52 // Constructor.
53 // Initialize all constant values which have to be used
54 // during Reco algorithm execution
55 //
56
2a71e873 57 for(Int_t i = 0; i < 15 ; i++) {fMisalTransShift[i] = 0.; fMisalRotShift[i] = 0.; }
d9b3567c 58 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = 0.;
59 //By default kPi0GammaGamma case
60 fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
61 fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
62 fNonLinearityParams[2] = 1.046;
63
64}
65
66//______________________________________________________________________
67AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
094786cc 68: TNamed(reco), fNonLinearityFunction(reco.fNonLinearityFunction),
69 fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo),
70 fW0(reco.fW0), fRecalibration(reco.fRecalibration),fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors)
d9b3567c 71{
72 //Copy ctor
73
2a71e873 74 for(Int_t i = 0; i < 15 ; i++) {fMisalRotShift[i] = reco.fMisalRotShift[i]; fMisalTransShift[i] = reco.fMisalTransShift[i]; }
d9b3567c 75 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
76}
77
78
79//______________________________________________________________________
80AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco)
81{
82 //Assignment operator
83
84 if(this == &reco)return *this;
85 ((TNamed *)this)->operator=(reco);
86
87 fNonLinearityFunction = reco.fNonLinearityFunction;
094786cc 88 fParticleType = reco.fParticleType;
89 fPosAlgo = reco.fPosAlgo;
90 fW0 = reco.fW0;
91 fRecalibration = reco.fRecalibration;
92 fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
93
2a71e873 94 for(Int_t i = 0; i < 15 ; i++) {fMisalTransShift[i] = reco.fMisalTransShift[i]; fMisalRotShift[i] = reco.fMisalRotShift[i];}
d9b3567c 95 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
96
97 return *this;
98}
99
100
101//__________________________________________________
094786cc 102AliEMCALRecoUtils::~AliEMCALRecoUtils()
103{
104 //Destructor.
105
106 if(fEMCALRecalibrationFactors) {
107 fEMCALRecalibrationFactors->Clear();
108 delete fEMCALRecalibrationFactors;
109 }
110}
111
112
113//__________________________________________________
d9b3567c 114Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
115// Correct cluster energy from non linearity functions
116 Float_t energy = cluster->E();
117
118 switch (fNonLinearityFunction) {
119
120 case kPi0MC:
121 //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
122 //Double_t par0 = 1.001;
123 //Double_t par1 = -0.01264;
124 //Double_t par2 = -0.03632;
125 //Double_t par3 = 0.1798;
126 //Double_t par4 = -0.522;
127 energy /= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
128 ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
129 exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
130 break;
131
132 case kPi0GammaGamma:
133
134 //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
135 //Double_t par0 = 0.1457;
136 //Double_t par1 = -0.02024;
137 //Double_t par2 = 1.046;
138 energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
139 break;
140
141 case kPi0GammaConversion:
142
143 //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
144 //Double_t C = 0.139393/0.1349766;
145 //Double_t a = 0.0566186;
146 //Double_t b = 0.982133;
147 energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
148
149 break;
150
151 case kNoCorrection:
152 AliDebug(2,"No correction on the energy\n");
153 break;
154
155 }
156
157 return energy;
158
159}
160
161//__________________________________________________
094786cc 162Float_t AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const
163{
164 //Calculate shower depth for a given cluster energy and particle type
165
166 // parameters
167 Float_t x0 = 1.23;
168 Float_t ecr = 8;
169 Float_t depth = 0;
170
171 switch ( iParticle )
172 {
173 case kPhoton:
174 depth = x0 * TMath::Log( (energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
175 break;
176
177 case kElectron:
178 depth = x0 * TMath::Log( (energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
179 break;
180
181 case kHadron:
182 // hadron
183 // boxes anc. here
184 if(gGeoManager){
185 gGeoManager->cd("ALIC_1/XEN1_1");
186 TGeoNode *geoXEn1 = gGeoManager->GetCurrentNode();
187 TGeoNodeMatrix *geoSM = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
188 TGeoVolume *geoSMVol = geoSM->GetVolume();
189 TGeoShape *geoSMShape = geoSMVol->GetShape();
190 TGeoBBox *geoBox = dynamic_cast<TGeoBBox *>(geoSMShape);
191 Float_t l = geoBox->GetDX()*2 ;
192 depth = 0.5 * l;
193 }
194 else{//electron
195 depth = x0 * TMath::Log( (energy*1000 / ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
196 }
197
198 break;
199
200 default://photon
201 depth = x0 * TMath::Log( (energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
202 }
203
204 return depth;
205
206}
207
208//__________________________________________________
209void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu, Int_t & absId, Int_t& iSupMod, Int_t& ieta, Int_t& iphi)
d9b3567c 210{
211 //For a given CaloCluster gets the absId of the cell
212 //with maximum energy deposit.
213
214 Double_t eMax = -1.;
215 Double_t eCell = -1.;
094786cc 216 Float_t fraction = 1.;
217 Float_t recalFactor = 1.;
d9b3567c 218 Int_t cellAbsId = -1 ;
094786cc 219
d9b3567c 220 Int_t iTower = -1;
221 Int_t iIphi = -1;
222 Int_t iIeta = -1;
223
224 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
094786cc 225 cellAbsId = clu->GetCellAbsId(iDig);
226 fraction = clu->GetCellAmplitudeFraction(iDig);
227 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
228 if(IsRecalibrationOn()) {
229 geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
230 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
231 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
232 }
233 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
234
235 if(eCell > eMax) {
d9b3567c 236 eMax = eCell;
237 absId = cellAbsId;
238 //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
239 }
240 }// cell loop
241
242 //Get from the absid the supermodule, tower and eta/phi numbers
243 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
244 //Gives SuperModule and Tower numbers
245 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
246 iIphi, iIeta,iphi,ieta);
247
248}
249
094786cc 250//________________________________________________________________
251void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
252 //Init EMCAL recalibration factors
253 AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
254 //In order to avoid rewriting the same histograms
255 Bool_t oldStatus = TH1::AddDirectoryStatus();
256 TH1::AddDirectory(kFALSE);
257
258 fEMCALRecalibrationFactors = new TObjArray(12);
259 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));
260 //Init the histograms with 1
261 for (Int_t sm = 0; sm < 12; sm++) {
262 for (Int_t i = 0; i < 48; i++) {
263 for (Int_t j = 0; j < 24; j++) {
264 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
265 }
266 }
267 }
268 fEMCALRecalibrationFactors->SetOwner(kTRUE);
269 fEMCALRecalibrationFactors->Compress();
270
271 //In order to avoid rewriting the same histograms
272 TH1::AddDirectory(oldStatus);
273}
274
275
276//________________________________________________________________
277void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells){
278 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
279
280 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
281 UShort_t * index = cluster->GetCellsAbsId() ;
282 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
283 Int_t ncells = cluster->GetNCells();
284
285 //Initialize some used variables
286 Float_t energy = 0;
287 Int_t absId = -1;
288 Int_t icol = -1, irow = -1, imod=1;
289 Float_t factor = 1, frac = 0;
290
291 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
292 for(Int_t icell = 0; icell < ncells; icell++){
293 absId = index[icell];
294 frac = fraction[icell];
295 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
296 Int_t iTower = -1, iIphi = -1, iIeta = -1;
297 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
298 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
299 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
300 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
301 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
302 imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
303
304 energy += cells->GetCellAmplitude(absId)*factor*frac;
305 }
306
307
308 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
309
310 cluster->SetE(energy);
311
312}
313
314
d9b3567c 315//__________________________________________________
094786cc 316void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
d9b3567c 317{
318 //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
319
094786cc 320 if (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
321 else if(fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
322
323}
324
325//__________________________________________________
326void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
327{
328 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
329 // The algorithm is a copy of what is done in AliEMCALRecPoint
330
331 Double_t eCell = 0.;
332 Float_t fraction = 1.;
333 Float_t recalFactor = 1.;
334
335 Int_t absId = -1;
336 Int_t iTower = -1, iIphi = -1, iIeta = -1;
337 Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
338 Float_t weight = 0., totalWeight=0.;
339 Float_t newPos[3] = {0,0,0};
340 Double_t pLocal[3], pGlobal[3];
341
342 Float_t clEnergy = clu->E(); //Energy already recalibrated previously
343
344 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi);
345 Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
346
347 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
348 absId = clu->GetCellAbsId(iDig);
349 fraction = clu->GetCellAmplitudeFraction(iDig);
350 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
351 geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
352 geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
353
354 if(IsRecalibrationOn()) {
355 recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
356 }
357 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
358
359 weight = GetCellWeight(eCell,clEnergy);
360 totalWeight += weight;
361 geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
362 geom->GetGlobal(pLocal,pGlobal,iSupModMax);
363
364 for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
365
366 }// cell loop
367
368 if(totalWeight>0){
369 for(int i=0; i<3; i++ ) newPos[i] /= totalWeight;
370 }
371
372 //printf("iSM %d \n",iSupMod);
373 //Float_t pos[]={0,0,0};
374 //clu->GetPosition(pos);
375 //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
376
377
378 //printf("NewPos a: %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
379
380 if(iSupModMax > 1) {//sector 1
381 newPos[0] +=fMisalTransShift[3];//-=3.093;
382 newPos[1] +=fMisalTransShift[4];//+=6.82;
383 newPos[2] +=fMisalTransShift[5];//+=1.635;
384 }
385 else {//sector 0
386 newPos[0] +=fMisalTransShift[0];//+=1.134;
387 newPos[1] +=fMisalTransShift[1];//+=8.2;
388 newPos[2] +=fMisalTransShift[2];//+=1.197;
389 }
390
391 clu->SetPosition(newPos);
392
393
394
395}
396
397//__________________________________________________
398void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
399{
400 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
401 // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
402
403 Double_t eCell = 1.;
404 Float_t fraction = 1.;
405 Float_t recalFactor = 1.;
406
407 Int_t absId = -1;
d9b3567c 408 Int_t iTower = -1;
094786cc 409 Int_t iIphi = -1, iIeta = -1;
410 Int_t iSupMod = -1, iSupModMax = -1;
d9b3567c 411 Int_t iphi = -1, ieta =-1;
412
413 Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
094786cc 414 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi);
415 Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
416
d9b3567c 417 Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
094786cc 418 Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
419 Int_t startingSM = -1;
d9b3567c 420
421 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
094786cc 422 absId = clu->GetCellAbsId(iDig);
423 fraction = clu->GetCellAmplitudeFraction(iDig);
424 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
425 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
d9b3567c 426 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
094786cc 427
d9b3567c 428 if (iDig==0) startingSM = iSupMod;
429 else if(iSupMod != startingSM) areInSameSM = kFALSE;
094786cc 430
431 eCell = cells->GetCellAmplitude(absId);
d9b3567c 432
094786cc 433 if(IsRecalibrationOn()) {
434 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
435 }
436 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
d9b3567c 437
094786cc 438 weight = GetCellWeight(eCell,clEnergy);
d9b3567c 439 if(weight < 0) weight = 0;
440 totalWeight += weight;
441 weightedCol += ieta*weight;
442 weightedRow += iphi*weight;
443
444 //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
445
094786cc 446 }// cell loop
447
d9b3567c 448 Float_t xyzNew[]={0.,0.,0.};
449 if(areInSameSM == kTRUE) {
450 //printf("In Same SM\n");
451 weightedCol = weightedCol/totalWeight;
452 weightedRow = weightedRow/totalWeight;
094786cc 453 geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
d9b3567c 454 }
455 else {
456 //printf("In Different SM\n");
094786cc 457 geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
d9b3567c 458 }
d9b3567c 459
094786cc 460 clu->SetPosition(xyzNew);
d9b3567c 461
462}
463
464//__________________________________________________
465void AliEMCALRecoUtils::Print(const Option_t *) const
466{
467 // Print Parameters
468
469 printf("AliEMCALRecoUtils Settings: \n");
470 printf("Misalignment shifts\n");
2a71e873 471 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,
472 fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
473 fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
d9b3567c 474 printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
475 for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
094786cc 476
477 printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
d9b3567c 478
479}