Updated materials in geometry (M. Sitta)
[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
33
34// STEER includes
35#include "AliEMCALRecoUtils.h"
36#include "AliEMCALGeoUtils.h"
37#include "AliVCluster.h"
38#include "AliVCaloCells.h"
39#include "AliLog.h"
40
41ClassImp(AliEMCALRecoUtils)
42
43//______________________________________________
44AliEMCALRecoUtils::AliEMCALRecoUtils():
45 fNonLinearityFunction (kPi0GammaGamma)
46{
47//
48 // Constructor.
49 // Initialize all constant values which have to be used
50 // during Reco algorithm execution
51 //
52
2a71e873 53 for(Int_t i = 0; i < 15 ; i++) {fMisalTransShift[i] = 0.; fMisalRotShift[i] = 0.; }
d9b3567c 54 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = 0.;
55 //By default kPi0GammaGamma case
56 fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
57 fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
58 fNonLinearityParams[2] = 1.046;
59
60}
61
62//______________________________________________________________________
63AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
64: TNamed(reco), fNonLinearityFunction(reco.fNonLinearityFunction)
65{
66 //Copy ctor
67
2a71e873 68 for(Int_t i = 0; i < 15 ; i++) {fMisalRotShift[i] = reco.fMisalRotShift[i]; fMisalTransShift[i] = reco.fMisalTransShift[i]; }
d9b3567c 69 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
70}
71
72
73//______________________________________________________________________
74AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco)
75{
76 //Assignment operator
77
78 if(this == &reco)return *this;
79 ((TNamed *)this)->operator=(reco);
80
81 fNonLinearityFunction = reco.fNonLinearityFunction;
2a71e873 82 for(Int_t i = 0; i < 15 ; i++) {fMisalTransShift[i] = reco.fMisalTransShift[i]; fMisalRotShift[i] = reco.fMisalRotShift[i];}
d9b3567c 83 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
84
85 return *this;
86}
87
88
89//__________________________________________________
90Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
91// Correct cluster energy from non linearity functions
92 Float_t energy = cluster->E();
93
94 switch (fNonLinearityFunction) {
95
96 case kPi0MC:
97 //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
98 //Double_t par0 = 1.001;
99 //Double_t par1 = -0.01264;
100 //Double_t par2 = -0.03632;
101 //Double_t par3 = 0.1798;
102 //Double_t par4 = -0.522;
103 energy /= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
104 ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
105 exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
106 break;
107
108 case kPi0GammaGamma:
109
110 //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
111 //Double_t par0 = 0.1457;
112 //Double_t par1 = -0.02024;
113 //Double_t par2 = 1.046;
114 energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
115 break;
116
117 case kPi0GammaConversion:
118
119 //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
120 //Double_t C = 0.139393/0.1349766;
121 //Double_t a = 0.0566186;
122 //Double_t b = 0.982133;
123 energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
124
125 break;
126
127 case kNoCorrection:
128 AliDebug(2,"No correction on the energy\n");
129 break;
130
131 }
132
133 return energy;
134
135}
136
137//__________________________________________________
138void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeoUtils *geom, AliVCaloCells* cells, AliVCluster* clu, Int_t & absId, Int_t& iSupMod, Int_t& ieta, Int_t& iphi)
139{
140 //For a given CaloCluster gets the absId of the cell
141 //with maximum energy deposit.
142
143 Double_t eMax = -1.;
144 Double_t eCell = -1.;
145 Int_t cellAbsId = -1 ;
146 Int_t iTower = -1;
147 Int_t iIphi = -1;
148 Int_t iIeta = -1;
149
150 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
151 cellAbsId = clu->GetCellAbsId(iDig);
152 eCell = cells->GetCellAmplitude(cellAbsId);
153 if(eCell > eMax) {
154 eMax = eCell;
155 absId = cellAbsId;
156 //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
157 }
158 }// cell loop
159
160 //Get from the absid the supermodule, tower and eta/phi numbers
161 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
162 //Gives SuperModule and Tower numbers
163 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
164 iIphi, iIeta,iphi,ieta);
165
166}
167
168//__________________________________________________
169void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeoUtils *geom, AliVCaloCells* cells, AliVCluster* clu, Int_t iParticle)
170{
171 //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
172
173 Double_t eMax = -1.;
174 Double_t eCell = -1.;
175 Int_t cellAbsId = -1;
176
177 Int_t maxId = -1;
178 Int_t iTower = -1;
179 Int_t iIphi = -1;
180 Int_t iIeta = -1;
181 Int_t iSupMod = -1;
182 Int_t iphi = -1, ieta =-1;
183
184 Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
185 Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
186 Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
187 Int_t startingSM = -1;
188
189 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
190 cellAbsId = clu->GetCellAbsId(iDig);
191 geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
192 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
193 if (iDig==0) startingSM = iSupMod;
194 else if(iSupMod != startingSM) areInSameSM = kFALSE;
195
196 eCell = cells->GetCellAmplitude(cellAbsId);
197
198 weight = TMath::Log(eCell/clEnergy) + 4;
199 if(weight < 0) weight = 0;
200 totalWeight += weight;
201 weightedCol += ieta*weight;
202 weightedRow += iphi*weight;
203
204 //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
205
206 if(eCell > eMax) {
207 eMax = eCell;
208 maxId = cellAbsId;
209 //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
210 }
211 }// cell loop
212
213 //Get from the absid the supermodule, tower and eta/phi numbers
214 geom->GetCellIndex(maxId,iSupMod,iTower,iIphi,iIeta);
215 //Gives SuperModule and Tower numbers
216 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
217 iIphi, iIeta,iphi,ieta);
218
219 Float_t xyzNew[]={0.,0.,0.};
220 if(areInSameSM == kTRUE) {
221 //printf("In Same SM\n");
222 weightedCol = weightedCol/totalWeight;
223 weightedRow = weightedRow/totalWeight;
2a71e873 224 geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupMod, clEnergy, iParticle, fMisalTransShift, fMisalRotShift, xyzNew);
d9b3567c 225 }
226 else {
227 //printf("In Different SM\n");
2a71e873 228 geom->RecalculateTowerPosition(iphi, ieta, iSupMod, clEnergy, iParticle, fMisalTransShift, fMisalRotShift, xyzNew);
d9b3567c 229 }
230 clu->SetPosition(xyzNew);
231
232
233 //printf("\t Max : cell %d, iSupMod %d, ieta %d, iphi %d \n",maxId,iSupMod, ieta,iphi);
234
235}
236
237//__________________________________________________
238void AliEMCALRecoUtils::Print(const Option_t *) const
239{
240 // Print Parameters
241
242 printf("AliEMCALRecoUtils Settings: \n");
243 printf("Misalignment shifts\n");
2a71e873 244 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,
245 fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
246 fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
d9b3567c 247 printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
248 for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
249
250}