BOSS 7.0.3
BESIII Offline Software System
Loading...
Searching...
No Matches
EmcRecShowerPosLog.cxx
Go to the documentation of this file.
1//
2// Logarithmic weighted position attribute reconstruction
3//
4// Created by Zhe Wang 2004, 5, 16
5//
6#include "EmcRec/EmcRecShowerPosLog.h"
7
8#include <iostream>
9#include <fstream>
10
11#include "EmcRecGeoSvc/EmcRecGeoSvc.h"
12#include "GaudiKernel/Bootstrap.h"
13#include "GaudiKernel/ISvcLocator.h"
14
15#include "EmcRec/EmcRecParameter.h"
16
17
18void EmcRecShowerPosLog::Position(RecEmcShower& aShower)
19{
20 RecEmcFractionMap::const_iterator cit;
21 HepPoint3D pos(0,0,0);
22 HepPoint3D possum(0,0,0);
23 double w,wsum=0;
24
25 // cout<<"EmcRecShowerPosLog::Position Begin"<<endl;
26
27 IEmcRecGeoSvc* iGeoSvc;
28 ISvcLocator* svcLocator = Gaudi::svcLocator();
29 StatusCode sc = svcLocator->service("EmcRecGeoSvc",iGeoSvc);
30 if(sc!=StatusCode::SUCCESS) {
31 // cout<<"Error: Can't get EmcRecGeoSvc"<<endl;
32 }
33
34 EmcRecParameter& Para=EmcRecParameter::GetInstance();
35
36 double offset=Para.LogPosOffset();
37 //cout<<offset<<endl;
38
39 double e5x5 = -99999;
40 if(Para.PositionMode2()=="all") {
41 double etot=aShower.getEAll();
42 for(cit=aShower.Begin();cit!=aShower.End();cit++){
43 w=offset+log((cit->second.getEnergy()/etot)*cit->second.getFraction());
44 if(w>0){
45 wsum+=w;
46 pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
47 possum+=pos*w;
48 }
49 }
50 } else if(Para.PositionMode2()=="3x3") {
51 double e3x3=aShower.e3x3();
52 RecEmcFractionMap fracMap3x3=aShower.getFractionMap3x3();
53 for(cit=fracMap3x3.begin();cit!=fracMap3x3.end();cit++) {
54 w=offset+log((cit->second.getEnergy()/e3x3)*cit->second.getFraction());
55 if(w>0){
56 wsum+=w;
57 pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
58 possum+=pos*w;
59 }
60 }
61 } else {
62 e5x5=aShower.e5x5();
63 RecEmcFractionMap fracMap5x5=aShower.getFractionMap5x5();
64 for(cit=fracMap5x5.begin();cit!=fracMap5x5.end();cit++) {
65 w=offset+log((cit->second.getEnergy()/e5x5)*cit->second.getFraction());
66 if(w>0){
67 wsum+=w;
68 pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
69 possum+=pos*w;
70 }
71 }
72 }
73 if(wsum<=0) {
74 // cout<<"[[[[[[[[[[[[[[["<<wsum<<endl;
75 }
76 pos=possum/wsum;
77
78 // aShower.setPositionNoCor(pos);
79 //----------------------------------------------------------------------
80 //position correction
81 RecEmcID id = aShower.getShowerId();
82 const unsigned int module = EmcID::barrel_ec(id);
83 const unsigned int thetaModule = EmcID::theta_module(id);
84 const unsigned int phiModule = EmcID::phi_module(id);
85 HepPoint3D posCorr(pos);
86
87 if(module==1) { //barrel
88
89
90 double IRShift=0, parTheta[5],parPhi[5];
91
92 unsigned int theModule;
93 if(thetaModule>21) {
94 theModule = 43 - thetaModule;
95 id = EmcID::crystal_id(module,theModule,phiModule);
96 posCorr.setZ(-posCorr.z());
97 } else {
98 theModule = thetaModule;
99 }
100
101
102 if(theModule==21) {
103 IRShift = 0;
104 } else if(theModule==20) {
105 IRShift = 2.5;
106 } else {
107 IRShift = 5.;
108 }
109 HepPoint3D IR(0,0,IRShift);
110
111 //MethodMode=0 old parameters
112 //MethodMode=1 new parameters
113 if(Para.MethodMode()==0){ //using old parameter given by hemiao
114
115 for(int i=0;i<5;i++){
116 if(theModule==21) {
117 double par[5] = { 51.97, -0.1209, 6.175, -0.1431, -0.005497 };
118 parTheta[i] = par[i];
119 }else if(theModule==20){
120 double par[5] = { 56.31, -0.1135, 6.312, -1.216, -0.02978 };
121 parTheta[i] = par[i];
122 }else{
123 double par[5] = { 10.65, -0.2257, 2.216, -0.1562, -0.05552 };
124 parTheta[i] = par[i];
125 }
126 }
127
128 } else if(Para.MethodMode()==1){ //using new parameter given by liuchunxiu
129
130 for(int i=0;i<5;i++){
131
132 if(e5x5>1.0) parTheta[i] = Para.BarrLogThetaPara(theModule,i);
133 else if(e5x5<=1.0 && e5x5>0.5) parTheta[i] = Para.BarrLogThetaPara(theModule+22,i);
134 else if(e5x5<=0.5) parTheta[i] = Para.BarrLogThetaPara(theModule+44,i);
135
136 if(e5x5>1.0) parPhi[i] = Para.BarrLogPhiPara(theModule,i);
137 else if(e5x5<=1.0 && e5x5>0.5) parPhi[i] = Para.BarrLogPhiPara(theModule+22,i);
138 else if(e5x5<=0.5) parPhi[i] = Para.BarrLogPhiPara(theModule+44,i);
139
140
141 }
142
143 }
144
145 HepPoint3D center01, center23; //center of point0,1 and point2,3 in crystal
146 center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
147 center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
148
149 double theta01,theta23,length,d;
150 theta01 = (center01-IR).theta();
151 theta23 = (center23-IR).theta();
152 length = (center01-IR).mag();
153 d = length*tan(theta01-theta23); //length in x direction
154
155 double xIn,xOut,deltaTheta;
156 xIn = length*tan(theta01-(posCorr-IR).theta())-d/2;
157 xOut = xIn-( parTheta[0]*atan(xIn*parTheta[1]-parTheta[4])+ parTheta[2]*xIn + parTheta[3] );
158 deltaTheta = atan((xOut+d/2)/length)-atan((xIn+d/2)/length);
159
160 //obtained by photon, not used, just for comparison
161 // double parPhi[5] = { 51.1, -0.1499, 7.62, 0.1301, 0.005581 };
162 double yIn,deltaY,deltaPhi;
163
164 //1.843=1.5+0.343 degree
165 if(phiModule==0){
166 yIn = posCorr.phi()*180./CLHEP::pi -phiModule*3. - 1.843;
167 }else if (phiModule==119){
168 yIn = posCorr.phi()*180./CLHEP::pi + 360.-phiModule*3. -1.843;
169 }else {
170 yIn = posCorr.phi() < 0 ? posCorr.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
171 : posCorr.phi()*180./CLHEP::pi-phiModule*3.-1.843;
172 }
173
174 //deltaY=yIn-Ymc => Ymc=yIn-deltaY
175 deltaY = parPhi[0]*atan(yIn*parPhi[1]-parPhi[4])+parPhi[2]*yIn+parPhi[3];
176 deltaPhi = deltaY*CLHEP::pi/180.;
177
178 HepPoint3D rotateVector(0,1,0); //y axis
179 rotateVector.rotateZ(center01.phi());
180 posCorr.setZ(posCorr.z()-IRShift);
181 posCorr.rotate(-deltaTheta,rotateVector);
182 posCorr.setZ(posCorr.z()+IRShift);
183 //phi parameter is gotten by the average of e+ e- peak
184 if(Para.MethodMode()==0){
185 posCorr.rotateZ(-0.002994);
186 }else if(Para.MethodMode()==1){
187 posCorr.rotateZ(-deltaPhi);
188 }
189
190 if(thetaModule>21) {
191 posCorr.setZ(-posCorr.z());
192 }
193
194 double lengthemc;
195 lengthemc = abs(posCorr.z()/cos(posCorr.theta()));
196 // for Data
197 double posDataCorPar;
198 if(Para.DataMode()==1) {
199 posDataCorPar=Para.BarrPosDataCor(thetaModule,phiModule);
200 posCorr.setTheta(posCorr.theta()-posDataCorPar/lengthemc);
201 }
202
203 //for MC
204 double posMCCorPar;
205 if(Para.DataMode()==0) {
206 posMCCorPar=Para.BarrPosMCCor(thetaModule,phiModule);
207 posCorr.setTheta(posCorr.theta()-posMCCorPar/lengthemc);
208 }
209
210 } else { //endcap position corretion
211 if(Para.MethodMode()==1){
212
213 double IRShift=0, parTheta[5],parPhi[5];
214
215 if(module==0){ //east endcap
216
217 IRShift = 10;
218
219 for(int i=0;i<5;i++) {
220
221
222 if(e5x5>1.0) parTheta[i] = Para.EastLogThetaPara(thetaModule,i);
223 else if(e5x5<=1.0 && e5x5>0.5) parTheta[i] = Para.EastLogThetaPara(thetaModule+6,i);
224 else if (e5x5<=0.5) parTheta[i] = Para.EastLogThetaPara(thetaModule+6,i);
225
226 if(e5x5>1.0) parPhi[i] = Para.EastLogPhiPara(0,i);
227 else if(e5x5<=1.0 && e5x5>0.5) parPhi[i] = Para.EastLogPhiPara(1,i);
228 else if (e5x5<=0.5) parPhi[i] = Para.EastLogPhiPara(2,i);
229 }
230
231 }else if (module==2){ //west endcap
232 IRShift = -10;
233 for(int i=0;i<5;i++){
234
235 if(e5x5>1.0) parTheta[i] = Para.WestLogThetaPara(thetaModule,i);
236 else if(e5x5<=1.0 && e5x5>0.5) parTheta[i] = Para.WestLogThetaPara(thetaModule+6,i);
237 else if (e5x5<=0.5) parTheta[i] = Para.WestLogThetaPara(thetaModule+6,i);
238
239 if(e5x5>1.0) parPhi[i] = Para.WestLogPhiPara(0,i);
240 else if(e5x5<=1.0 && e5x5>0.5) parPhi[i] = Para.WestLogPhiPara(1,i);
241 else if (e5x5<=0.5) parPhi[i] = Para.WestLogPhiPara(2,i);
242
243 }
244
245 }
246
247 HepPoint3D IR(0,0,IRShift);
248
249 double theta12=-9999,theta03=-9999,theta23=-9999,theta14=-9999,delta=-9999;
250 double phi01=-9999,phi23=-9999,phi12=-9999,phi34=-9999,deltaphi=-9999;
251 double xIn,deltaX,deltaTheta,yIn,deltaY,deltaPhi;
252
253 double posphi=posCorr.phi();
254 if(phiModule==0){
255 posphi = posphi;
256 }else if( ((thetaModule==0||thetaModule==1)&&phiModule==63)
257 ||((thetaModule==2||thetaModule==3)&&phiModule==79)
258 ||((thetaModule==4||thetaModule==5)&&phiModule==95)
259 ){
260 posphi = posphi+2*CLHEP::pi;
261 }else {
262 posphi = posphi < 0 ? posphi+2*CLHEP::pi : posphi;
263 }
264
265 HepPoint3D center,center01, center23, center12,center03,center34,center14;
266
267
268 if( (thetaModule==1&&((phiModule+3)%4 == 0||(phiModule+2)%4 == 0))
269 ||(thetaModule==3&&((phiModule+4)%5 == 0||(phiModule+3)%5 == 0||(phiModule+2)%5 == 0))) //pentagon
270 {
271
272 center23=(iGeoSvc->GetCrystalPoint(id,2)+ iGeoSvc->GetCrystalPoint(id,3))/2.0;
273 center14=(iGeoSvc->GetCrystalPoint(id,1)+ iGeoSvc->GetCrystalPoint(id,4))/2.0;
274 center=center23;
275 theta23 = (center23 - IR).theta();
276 theta14 = (center14 - IR).theta();
277 delta = theta14 - theta23;
278 xIn = (((posCorr-IR).theta() - theta23) - delta*0.5)/delta;
279
280
281 center12=(iGeoSvc->GetCrystalPoint(id,1)+ iGeoSvc->GetCrystalPoint(id,2))/2.0;
282 center34=(iGeoSvc->GetCrystalPoint(id,3)+ iGeoSvc->GetCrystalPoint(id,4))/2.0;
283 phi12 = center12.phi();
284 phi34 = center34.phi();
285
286 phi12 = phi12 <0 ? phi12+2*CLHEP::pi : phi12;
287 phi34 = phi34 <0 ? phi34+2*CLHEP::pi : phi34;
288 deltaphi = phi34 - phi12;
289 yIn = ((posphi - phi12) - deltaphi*0.5)/deltaphi;
290
291 } else {
292
293 center12=(iGeoSvc->GetCrystalPoint(id,1)+ iGeoSvc->GetCrystalPoint(id,2))/2.0;
294 center03=(iGeoSvc->GetCrystalPoint(id,0)+ iGeoSvc->GetCrystalPoint(id,3))/2.0;
295 center=center12;
296
297 theta12 = (center12 - IR).theta();
298 theta03 = (center03 - IR).theta();
299 delta = theta03 - theta12;
300 xIn = (((posCorr-IR).theta() - theta12) - delta*0.5)/delta;
301
302
303 center01=(iGeoSvc->GetCrystalPoint(id,0)+ iGeoSvc->GetCrystalPoint(id,1))/2.0;
304 center23=(iGeoSvc->GetCrystalPoint(id,2)+ iGeoSvc->GetCrystalPoint(id,3))/2.0;
305 phi01 = center01.phi();
306 phi23 = center23.phi();
307 phi01 = phi01 <0 ? phi01+2*CLHEP::pi : phi01;
308 phi23 = phi23 <0 ? phi23+2*CLHEP::pi : phi23;
309 deltaphi = phi23 - phi01;
310 yIn = ((posphi - phi01) - deltaphi*0.5)/deltaphi;
311
312 }
313 //deltaX=xIn-Xext => Xext=xIn-deltaX
314 deltaX = parTheta[0]*atan(xIn*parTheta[1]-parTheta[4]) + parTheta[2]*xIn + parTheta[3] ;
315 deltaTheta = deltaX*delta;
316 //deltaY=yIn-Ymc => Ymc=yIn-deltaY
317 deltaY = parPhi[0]*atan(yIn*parPhi[1]-parPhi[4]) + parPhi[2]*yIn + parPhi[3] ;
318 deltaPhi = deltaY*deltaphi;
319
320 HepPoint3D rotateVector(0,1,0); //y axis
321 rotateVector.rotateZ(center.phi());
322 posCorr.setZ(posCorr.z()-IRShift);
323 posCorr.rotate(-deltaTheta,rotateVector);
324 posCorr.setZ(posCorr.z()+IRShift);
325 posCorr.rotateZ(-deltaPhi);
326
327 double lengthemc;
328 lengthemc = abs(posCorr.z()/cos(posCorr.theta()));
329 double posDataCorPar;
330 if(Para.DataMode()==1) {
331 if(module==0) posDataCorPar=Para.EastPosDataCor(thetaModule,phiModule);
332 if(module==2) posDataCorPar=Para.WestPosDataCor(thetaModule,phiModule);
333 posCorr.setTheta(posCorr.theta()-posDataCorPar/lengthemc);
334 }
335 }
336
337 }
338
339
340 //PosCorr=0 without position correction
341 //PosCorr=1 with position correction
342 if(Para.PosCorr()==0){ aShower.setPosition(pos);}
343
344 if(Para.PosCorr()==1){
345 aShower.setPosition(posCorr);
346 }
347 ////////////////////////////
348 if(aShower.energy()<1e-5) return;
349
350 double r,theta,phi;
351 double dtheta,dphi,dx,dy,dz;
352
353 r = posCorr.mag();
354 theta = posCorr.theta();
355 phi = posCorr.phi();
356 //cout<<"energy: "<<aShower.energy()<<" theta: "<<theta<<" phi: "<<phi<<endl;
357
358 if(aShower.energy()>0) {
359 dtheta = Para.SigTheta(0)/sqrt(aShower.energy())+Para.SigTheta(1);
360 dphi = Para.SigPhi(0)/sqrt(aShower.energy())+Para.SigPhi(1);
361 }
362 else {
363 dtheta = 0;
364 dphi = 0;
365 }
366 //cout<<"dtheta: "<<dtheta<<" dphi: "<<dphi<<endl;
367
368 dx = fabs(iGeoSvc->GetBarrelR()*sin(phi)*dphi);
369 dy = fabs(iGeoSvc->GetBarrelR()*cos(phi)*dphi);
370 if(theta>0)
371 dz = fabs(iGeoSvc->GetBarrelR()*dtheta/(sin(theta)*sin(theta)));
372 else
373 dz = 0;
374 //cout<<" dx: "<<dx<<" dy: "<<dy<<" dz: "<<dz<<endl;
375
376 aShower.setDtheta(dtheta);
377 aShower.setDphi(dphi);
378
379 //----------------------------------------------------------------------
380 // Error matrix
381 HepSymMatrix matrix(3); //error matrix of r, theta, phi
382 matrix[0][0]=0;
383 matrix[1][1]=dtheta*dtheta;
384 matrix[2][2]=dphi*dphi;
385 matrix[0][1]=0;
386 matrix[0][2]=0;
387 matrix[1][2]=0;
388 //cout<<"matrix:\n"<<matrix<<endl;
389
390 HepMatrix matrix1(3,3),matrix2(3,3);
391 matrix1[0][0]=sin(theta)*cos(phi);
392 matrix1[0][1]=r*cos(theta)*cos(phi);
393 matrix1[0][2]=-r*sin(theta)*sin(phi);
394 matrix1[1][0]=sin(theta)*sin(phi);
395 matrix1[1][1]=r*cos(theta)*sin(phi);
396 matrix1[1][2]=r*sin(theta)*cos(phi);
397 matrix1[2][0]=cos(theta);
398 matrix1[2][1]=-r*sin(theta);
399 matrix1[2][2]=0;
400 //cout<<"matrix1:\n"<<matrix1<<endl;
401
402 for(int i=0;i<3;i++) {
403 for(int j=0;j<3;j++) {
404 matrix2[i][j]=matrix1[j][i];
405 }
406 }
407 //cout<<"matrix2:\n"<<matrix2<<endl;
408
409 HepMatrix matrix4(3,3);
410 matrix4=matrix1*matrix*matrix2;
411 //cout<<"matrix4:\n"<<matrix4<<endl;
412
413 HepSymMatrix matrix5(3); //error matrix of x, y, z
414 for(int i=0;i<3;i++) {
415 for(int j=0;j<=i;j++) {
416 matrix5[i][j]=matrix4[i][j];
417 }
418 }
419 //cout<<"matrix5:\n"<<matrix5<<endl;
420
421 aShower.setErrorMatrix(matrix5);
422
423 if(matrix5[0][0]>0) dx=sqrt(matrix5[0][0]);
424 if(matrix5[1][1]>0) dy=sqrt(matrix5[1][1]);
425 if(matrix5[2][2]>0) dz=sqrt(matrix5[2][2]);
426}
427
etot
Double_t etot
Definition: DataBase/tau_mode.c:7
RecEmcFractionMap
map< RecEmcID, RecEmcFraction, less< RecEmcID > > RecEmcFractionMap
Definition: Emc/EmcRecEventModel/EmcRecEventModel-01-01-18/EmcRecEventModel/RecEmcFraction.h:56
w
double w
Definition: EvtPhokharaDef.hh:26
tan
double tan(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:216
sin
double sin(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:213
DstEmcShower::setErrorMatrix
void setErrorMatrix(const HepSymMatrix &error)
Definition: Event/DstEvent/DstEvent-00-02-51/DstEvent/DstEmcShower.h:75
EmcID::crystal_id
static Identifier crystal_id(const unsigned int barrel_ec, const unsigned int theta_module, const unsigned int phi_module)
For a single crystal.
Definition: EmcID.cxx:71
EmcID::barrel_ec
static unsigned int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0)
Definition: EmcID.cxx:38
EmcID::theta_module
static unsigned int theta_module(const Identifier &id)
Definition: EmcID.cxx:43
EmcID::phi_module
static unsigned int phi_module(const Identifier &id)
Definition: EmcID.cxx:48
EmcRecParameter::LogPosOffset
double LogPosOffset() const
Definition: EmcRecParameter.cxx:498
EmcRecParameter::BarrPosMCCor
double BarrPosMCCor(int ntheta, int nphi) const
Definition: EmcRecParameter.cxx:702
EmcRecParameter::GetInstance
static EmcRecParameter & GetInstance()
Definition: EmcRecParameter.cxx:461
EmcRecParameter::EastPosDataCor
double EastPosDataCor(int ntheta, int nphi) const
Definition: EmcRecParameter.cxx:697
EmcRecParameter::WestPosDataCor
double WestPosDataCor(int ntheta, int nphi) const
Definition: EmcRecParameter.cxx:692
EmcRecParameter::BarrLogPhiPara
double BarrLogPhiPara(int n, int m) const
Definition: EmcRecParameter.cxx:590
EmcRecParameter::PosCorr
double PosCorr() const
Definition: EmcRecParameter.cxx:519
EmcRecParameter::MethodMode
double MethodMode() const
Definition: EmcRecParameter.cxx:514
EmcRecParameter::SigTheta
double SigTheta(int n) const
Definition: EmcRecParameter.cxx:553
EmcRecParameter::WestLogThetaPara
double WestLogThetaPara(int n, int m) const
Definition: EmcRecParameter.cxx:646
EmcRecParameter::BarrPosDataCor
double BarrPosDataCor(int ntheta, int nphi) const
Definition: EmcRecParameter.cxx:687
EmcRecParameter::BarrLogThetaPara
double BarrLogThetaPara(int n, int m) const
Definition: EmcRecParameter.cxx:583
EmcRecParameter::WestLogPhiPara
double WestLogPhiPara(int n, int m) const
Definition: EmcRecParameter.cxx:656
EmcRecParameter::EastLogPhiPara
double EastLogPhiPara(int n, int m) const
Definition: EmcRecParameter.cxx:651
EmcRecParameter::EastLogThetaPara
double EastLogThetaPara(int n, int m) const
Definition: EmcRecParameter.cxx:641
EmcRecParameter::DataMode
double DataMode() const
Definition: EmcRecParameter.cxx:524
EmcRecParameter::SigPhi
double SigPhi(int n) const
Definition: EmcRecParameter.cxx:558
EmcRecShowerPosLog::Position
virtual void Position(RecEmcShower &aShower)
Definition: EmcRecShowerPosLog.cxx:18
IEmcRecGeoSvc::GetBarrelR
virtual double GetBarrelR() const =0
IEmcRecGeoSvc::GetCFrontCenter
virtual HepPoint3D GetCFrontCenter(const Identifier &id) const =0
IEmcRecGeoSvc::GetCrystalPoint
virtual HepPoint3D GetCrystalPoint(const Identifier &id, const int i) const =0
RecEmcShower::End
RecEmcFractionMap::const_iterator End() const
Definition: RecEmcShower.cxx:136
RecEmcShower::getFractionMap5x5
RecEmcFractionMap getFractionMap5x5() const
Definition: RecEmcShower.cxx:253
RecEmcShower::Begin
RecEmcFractionMap::const_iterator Begin() const
Definition: RecEmcShower.cxx:131
RecEmcShower::getFractionMap3x3
RecEmcFractionMap getFractionMap3x3() const
Definition: RecEmcShower.cxx:229