BOSS 7.0.6
BESIII Offline Software System
Loading...
Searching...
No Matches
EmcRecShowerPosLogShMax.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//
7
8#include <iostream>
9#include <fstream>
10
12#include "GaudiKernel/Bootstrap.h"
13#include "GaudiKernel/ISvcLocator.h"
14
16
17
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
35
36 double offset=Para.LogPosOffset();
37 //int NX0=Para.DepthOfShowerMax(); //unit X_0=1.86
38 //cout<<offset<<endl;
39
40 double e5x5 = -99999;
41 if(Para.PositionMode2()=="all") {
42 double etot=aShower.getEAll();
43 for(cit=aShower.Begin();cit!=aShower.End();cit++){
44 w=offset+log((cit->second.getEnergy()/etot)*cit->second.getFraction());
45 if(w>0){
46 wsum+=w;
47 pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
48 possum+=pos*w;
49 }
50 }
51 } else if(Para.PositionMode2()=="3x3") {
52 double e3x3=aShower.e3x3();
53 RecEmcFractionMap fracMap3x3=aShower.getFractionMap3x3();
54 for(cit=fracMap3x3.begin();cit!=fracMap3x3.end();cit++) {
55 w=offset+log((cit->second.getEnergy()/e3x3)*cit->second.getFraction());
56 if(w>0){
57 wsum+=w;
58 pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
59 possum+=pos*w;
60 }
61 }
62 } else {
63 e5x5=aShower.e5x5();
64 RecEmcFractionMap fracMap5x5=aShower.getFractionMap5x5();
65 for(cit=fracMap5x5.begin();cit!=fracMap5x5.end();cit++) {
66 w=offset+log((cit->second.getEnergy()/e5x5)*cit->second.getFraction());
67 if(w>0){
68 wsum+=w;
69 //pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
70
71 HepPoint3D Rc,R;
72 R=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
73 Rc=iGeoSvc->GetCCenter(cit->second.getCellId());
74 int NX0=6;
75 pos = R + (Rc-R)/(Rc-R).mag()*NX0*1.86; //X_CsI=1.86cm
76
77 possum+=pos*w;
78 }
79 }
80 }
81 if(wsum<=0) {
82 // cout<<"[[[[[[[[[[[[[[["<<wsum<<endl;
83 }
84 pos=possum/wsum;
85
86
87 int NX0=6;
88 HepPoint3D posFront=pos-pos/pos.mag()*NX0*1.86; //X_CsI=1.86cm
89
90
91 // aShower.setPositionNoCor(pos);
92 //----------------------------------------------------------------------
93 //position correction
94 RecEmcID id = aShower.getShowerId();
95 const unsigned int module = EmcID::barrel_ec(id);
96 const unsigned int thetaModule = EmcID::theta_module(id);
97 const unsigned int phiModule = EmcID::phi_module(id);
98 //HepPoint3D posCorr(pos);
99
100 HepPoint3D posCorr(posFront);
101
102 if(module==1) { //barrel
103
104
105 double IRShift, parTheta[5],parPhi[5];
106
107
108 for(int i=0;i<5;i++){
109
110 if(e5x5>1.0) parTheta[i] = Para.BarrLogShMaxThetaPara(thetaModule,i);
111 else if(e5x5<=1.0 && e5x5>0.5) parTheta[i] = Para.BarrLogShMaxThetaPara(thetaModule+44,i);
112 else if(e5x5<=0.5) parTheta[i] = Para.BarrLogShMaxThetaPara(thetaModule+88,i);
113
114 if(e5x5>1.0) parPhi[i] = Para.BarrLogShMaxPhiPara(thetaModule,i);
115 else if(e5x5<=1.0 && e5x5>0.5) parPhi[i] = Para.BarrLogShMaxPhiPara(thetaModule+44,i);
116 else if(e5x5<=0.5) parPhi[i] = Para.BarrLogShMaxPhiPara(thetaModule+88,i);
117
118
119 }
120
121
122 HepPoint3D center01, center23, center12,center03; //center of point0,1 and point2,3 in crystal
123 center01 = (iGeoSvc->GetCrystalPoint(id,0)+ iGeoSvc->GetCrystalPoint(id,1))/2;
124 center23 = (iGeoSvc->GetCrystalPoint(id,2)+ iGeoSvc->GetCrystalPoint(id,3))/2;
125 center03 = (iGeoSvc->GetCrystalPoint(id,0)+ iGeoSvc->GetCrystalPoint(id,3))/2.0;
126 center12 = (iGeoSvc->GetCrystalPoint(id,1)+ iGeoSvc->GetCrystalPoint(id,2))/2.0;
127
128 double theta01,theta23,dtheta;
129
130 theta01 = (center01).theta();
131 theta23 = (center23).theta();
132
133 dtheta = theta01-theta23; //length in x direction
134 double xIn,xOut,deltaTheta,deltaX;
135 xIn = (((posCorr).theta()-theta23)-dtheta/2)/dtheta;
136 //emcX-extX=deltaX ==> extX=emcX-deltaX,Xout=xin-deltaX
137 deltaX = parTheta[0]*atan(xIn*parTheta[1]-parTheta[4]) + parTheta[2]*xIn + parTheta[3] ;
138 deltaTheta = deltaX*dtheta;
139
140 //obtained by photon, not used, just for comparison
141 // double parPhi[5] = { 51.1, -0.1499, 7.62, 0.1301, 0.005581 };
142
143 double phi12=-9999,phi03=-9999,deltaphi=-9999;
144
145 phi03 = center03.phi();
146 phi12 = center12.phi();
147
148
149 double posPhi;
150 posPhi=posCorr.phi();
151
152
153 if( phiModule==0){
154 posPhi = posPhi;
155 phi03 = phi03;
156 phi12 = phi12;
157
158 }else if( phiModule==119){
159 posPhi = posPhi+2*CLHEP::pi;
160 phi03 = phi03+2*CLHEP::pi;
161 phi12 = phi12+2*CLHEP::pi;
162
163 }else {
164 posPhi = posPhi <0 ? posPhi+2*CLHEP::pi : posPhi;
165 phi03 = phi03 <0 ? phi03+2*CLHEP::pi : phi03;
166 phi12 = phi12 <0 ? phi12+2*CLHEP::pi : phi12;
167 }
168
169 deltaphi = phi12 - phi03;
170
171
172 double yIn,deltaY,deltaPhi;
173 yIn = ((posPhi - phi03) - deltaphi*0.5)/deltaphi;
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*deltaphi;
177
178 // deltaPhi = deltaY*CLHEP::pi/180.;
179 /*
180 HepPoint3D rotateVector(0,1,0); //y axis
181 rotateVector.rotateZ(center01.phi());
182 posCorr.setZ(posCorr.z()-IRShift);
183 posCorr.rotate(-deltaTheta,rotateVector);
184 posCorr.setZ(posCorr.z()+IRShift);
185 */
186 /*
187 HepPoint3D possh;
188 possh=posCorr;//-IRShift;
189 //possh.setTheta(possh.theta()-deltaTheta);
190 // posCorr=possh+IRShift;
191 posCorr.setTheta(possh.theta()-deltaTheta);
192 */
193
194 HepPoint3D rotateVector(0,1,0); //y axis
195 //rotateVector.rotateZ(center01.phi());
196 rotateVector.rotateZ(center23.phi());
197
198 posCorr.rotate(-deltaTheta,rotateVector);
199 posCorr.rotateZ(-deltaPhi);
200
201 /*
202 double lengthemc;
203 lengthemc = abs(posCorr.z()/cos(posCorr.theta()));
204 // for Data
205 double posDataCorPar;
206 if(Para.DataMode()==1) {
207 posDataCorPar=Para.BarrPosDataCor(thetaModule,phiModule);
208 posCorr.setTheta(posCorr.theta()-posDataCorPar/lengthemc);
209 }
210 */
211 /* //for MC
212 double posMCCorPar;
213 if(Para.DataMode()==0) {
214 posMCCorPar=Para.BarrPosMCCor(thetaModule,phiModule);
215 posCorr.setTheta(posCorr.theta()-posMCCorPar/lengthemc);
216 }
217 */
218
219 } else { //endcap position corretion
220 // if(Para.MethodMode()==1){
221
222 double IRShift=0, parTheta[5],parPhi[5];
223
224 if(module==0){ //east endcap
225
226 IRShift = 10;
227
228 for(int i=0;i<5;i++) {
229
230
231 if(e5x5>1.0) parTheta[i] = Para.EastLogShMaxThetaPara(thetaModule,i);
232 else if(e5x5<=1.0 && e5x5>0.5) parTheta[i] = Para.EastLogShMaxThetaPara(thetaModule+6,i);
233 else if (e5x5<=0.5) parTheta[i] = Para.EastLogShMaxThetaPara(thetaModule+6,i);
234
235 if(e5x5>1.0) parPhi[i] = Para.EastLogShMaxPhiPara(0,i);
236 else if(e5x5<=1.0 && e5x5>0.5) parPhi[i] = Para.EastLogShMaxPhiPara(1,i);
237 else if (e5x5<=0.5) parPhi[i] = Para.EastLogShMaxPhiPara(2,i);
238 }
239
240 }else if (module==2){ //west endcap
241 IRShift = -10;
242 for(int i=0;i<5;i++){
243
244 if(e5x5>1.0) parTheta[i] = Para.WestLogShMaxThetaPara(thetaModule,i);
245 else if(e5x5<=1.0 && e5x5>0.5) parTheta[i] = Para.WestLogShMaxThetaPara(thetaModule+6,i);
246 else if (e5x5<=0.5) parTheta[i] = Para.WestLogShMaxThetaPara(thetaModule+6,i);
247
248 if(e5x5>1.0) parPhi[i] = Para.WestLogShMaxPhiPara(0,i);
249 else if(e5x5<=1.0 && e5x5>0.5) parPhi[i] = Para.WestLogShMaxPhiPara(1,i);
250 else if (e5x5<=0.5) parPhi[i] = Para.WestLogShMaxPhiPara(2,i);
251
252 }
253
254 }
255
256 HepPoint3D IR(0,0,IRShift);
257
258 double theta12=-9999,theta03=-9999,theta23=-9999,theta14=-9999,delta=-9999;
259 double phi01=-9999,phi23=-9999,phi12=-9999,phi34=-9999,deltaphi=-9999;
260 double xIn,deltaX,deltaTheta,yIn,deltaY,deltaPhi;
261
262 double posphi=posCorr.phi();
263 if(phiModule==0){
264 posphi = posphi;
265 }else if( ((thetaModule==0||thetaModule==1)&&phiModule==63)
266 ||((thetaModule==2||thetaModule==3)&&phiModule==79)
267 ||((thetaModule==4||thetaModule==5)&&phiModule==95)
268 ){
269 posphi = posphi+2*CLHEP::pi;
270 }else {
271 posphi = posphi < 0 ? posphi+2*CLHEP::pi : posphi;
272 }
273
274 HepPoint3D center,center01, center23, center12,center03,center34,center14;
275
276
277 if( (thetaModule==1&&((phiModule+3)%4 == 0||(phiModule+2)%4 == 0))
278 ||(thetaModule==3&&((phiModule+4)%5 == 0||(phiModule+3)%5 == 0||(phiModule+2)%5 == 0))) //pentagon
279 {
280
281 center23=(iGeoSvc->GetCrystalPoint(id,2)+ iGeoSvc->GetCrystalPoint(id,3))/2.0;
282 center14=(iGeoSvc->GetCrystalPoint(id,1)+ iGeoSvc->GetCrystalPoint(id,4))/2.0;
283 center=center23;
284 theta23 = (center23 - IR).theta();
285 theta14 = (center14 - IR).theta();
286 delta = theta14 - theta23;
287 xIn = (((posCorr-IR).theta() - theta23) - delta*0.5)/delta;
288
289
290 center12=(iGeoSvc->GetCrystalPoint(id,1)+ iGeoSvc->GetCrystalPoint(id,2))/2.0;
291 center34=(iGeoSvc->GetCrystalPoint(id,3)+ iGeoSvc->GetCrystalPoint(id,4))/2.0;
292 phi12 = center12.phi();
293 phi34 = center34.phi();
294
295 if(phiModule==0){
296 phi12 = phi12;
297 phi34 = phi34;
298 }else if( ((thetaModule==0||thetaModule==1)&&phiModule==63)
299 ||((thetaModule==2||thetaModule==3)&&phiModule==79)
300 ||((thetaModule==4||thetaModule==5)&&phiModule==95)
301 ){
302 phi12 = phi12+2*CLHEP::pi;
303 phi34 = phi34+2*CLHEP::pi;
304 }else {
305 phi12 = phi12 <0 ? phi12+2*CLHEP::pi : phi12;
306 phi34 = phi34 <0 ? phi34+2*CLHEP::pi : phi34;
307 }
308
309 deltaphi = phi34 - phi12;
310 yIn = ((posphi - phi12) - deltaphi*0.5)/deltaphi;
311
312 } else {
313
314 center12=(iGeoSvc->GetCrystalPoint(id,1)+ iGeoSvc->GetCrystalPoint(id,2))/2.0;
315 center03=(iGeoSvc->GetCrystalPoint(id,0)+ iGeoSvc->GetCrystalPoint(id,3))/2.0;
316 center=center12;
317
318 theta12 = (center12 - IR).theta();
319 theta03 = (center03 - IR).theta();
320 delta = theta03 - theta12;
321 xIn = (((posCorr-IR).theta() - theta12) - delta*0.5)/delta;
322
323
324 center01=(iGeoSvc->GetCrystalPoint(id,0)+ iGeoSvc->GetCrystalPoint(id,1))/2.0;
325 center23=(iGeoSvc->GetCrystalPoint(id,2)+ iGeoSvc->GetCrystalPoint(id,3))/2.0;
326 phi01 = center01.phi();
327 phi23 = center23.phi();
328
329
330 if(phiModule==0){
331 phi01 = phi01;
332 phi23 = phi23;
333 }else if( ((thetaModule==0||thetaModule==1)&&phiModule==63)
334 ||((thetaModule==2||thetaModule==3)&&phiModule==79)
335 ||((thetaModule==4||thetaModule==5)&&phiModule==95)
336 ){
337 phi01 = phi01+2*CLHEP::pi;
338 phi23 = phi23+2*CLHEP::pi;
339 }else {
340 phi01 = phi01 <0 ? phi01+2*CLHEP::pi : phi01;
341 phi23 = phi23 <0 ? phi23+2*CLHEP::pi : phi23;
342 }
343
344 deltaphi = phi23 - phi01;
345 yIn = ((posphi - phi01) - deltaphi*0.5)/deltaphi;
346
347 }
348 //deltaX=xIn-Xext => Xext=xIn-deltaX
349 deltaX = parTheta[0]*atan(xIn*parTheta[1]-parTheta[4]) + parTheta[2]*xIn + parTheta[3] ;
350 deltaTheta = deltaX*delta;
351 //deltaY=yIn-Ymc => Ymc=yIn-deltaY
352 deltaY = parPhi[0]*atan(yIn*parPhi[1]-parPhi[4]) + parPhi[2]*yIn + parPhi[3] ;
353 deltaPhi = deltaY*deltaphi;
354
355 HepPoint3D rotateVector(0,1,0); //y axis
356 rotateVector.rotateZ(center.phi());
357 posCorr.setZ(posCorr.z()-IRShift);
358 posCorr.rotate(-deltaTheta,rotateVector);
359 posCorr.setZ(posCorr.z()+IRShift);
360 posCorr.rotateZ(-deltaPhi);
361 /*
362 double lengthemc;
363 lengthemc = abs(posCorr.z()/cos(posCorr.theta()));
364 double posDataCorPar;
365 if(Para.DataMode()==1) {
366 if(module==0) posDataCorPar=Para.EastPosDataCor(thetaModule,phiModule);
367 if(module==2) posDataCorPar=Para.WestPosDataCor(thetaModule,phiModule);
368 posCorr.setTheta(posCorr.theta()-posDataCorPar/lengthemc);
369 }
370 */
371 //}
372
373 }
374
375
376 //PosCorr=0 without position correction
377 //PosCorr=1 with position correction
378 if(Para.PosCorr()==0){
379 //int NX0=6;
380 //HepPoint3D posFront=pos-pos/pos.mag()*NX0*1.86; //X_CsI=1.86cm
381 aShower.setPosition(posFront);
382 // aShower.setPosition(pos);
383
384 }
385
386 if(Para.PosCorr()==1){
387
388 //pos from the shmax to the front endface of crystal
389 // int NX0=6;
390 //HepPoint3D posCorrFront=posCorr-posCorr/posCorr.mag()*NX0*1.86; //X_CsI=1.86cm
391 //aShower.setPosition(posCorrFront);
392 ////////////
393
394 aShower.setPosition(posCorr);
395 }
396 ////////////////////////////
397 if(aShower.energy()<1e-5) return;
398
399 double r,theta,phi;
400 double dtheta,dphi,dx,dy,dz;
401
402 r = posCorr.mag();
403 theta = posCorr.theta();
404 phi = posCorr.phi();
405 //cout<<"energy: "<<aShower.energy()<<" theta: "<<theta<<" phi: "<<phi<<endl;
406
407 if(aShower.energy()>0) {
408 dtheta = Para.SigTheta(0)/sqrt(aShower.energy())+Para.SigTheta(1);
409 dphi = Para.SigPhi(0)/sqrt(aShower.energy())+Para.SigPhi(1);
410 }
411 else {
412 dtheta = 0;
413 dphi = 0;
414 }
415 //cout<<"dtheta: "<<dtheta<<" dphi: "<<dphi<<endl;
416
417 dx = fabs(iGeoSvc->GetBarrelR()*sin(phi)*dphi);
418 dy = fabs(iGeoSvc->GetBarrelR()*cos(phi)*dphi);
419 if(theta>0)
420 dz = fabs(iGeoSvc->GetBarrelR()*dtheta/(sin(theta)*sin(theta)));
421 else
422 dz = 0;
423 //cout<<" dx: "<<dx<<" dy: "<<dy<<" dz: "<<dz<<endl;
424
425 aShower.setDtheta(dtheta);
426 aShower.setDphi(dphi);
427
428 //----------------------------------------------------------------------
429 // Error matrix
430 HepSymMatrix matrix(3); //error matrix of r, theta, phi
431 matrix[0][0]=0;
432 matrix[1][1]=dtheta*dtheta;
433 matrix[2][2]=dphi*dphi;
434 matrix[0][1]=0;
435 matrix[0][2]=0;
436 matrix[1][2]=0;
437 //cout<<"matrix:\n"<<matrix<<endl;
438
439 HepMatrix matrix1(3,3),matrix2(3,3);
440 matrix1[0][0]=sin(theta)*cos(phi);
441 matrix1[0][1]=r*cos(theta)*cos(phi);
442 matrix1[0][2]=-r*sin(theta)*sin(phi);
443 matrix1[1][0]=sin(theta)*sin(phi);
444 matrix1[1][1]=r*cos(theta)*sin(phi);
445 matrix1[1][2]=r*sin(theta)*cos(phi);
446 matrix1[2][0]=cos(theta);
447 matrix1[2][1]=-r*sin(theta);
448 matrix1[2][2]=0;
449 //cout<<"matrix1:\n"<<matrix1<<endl;
450
451 for(int i=0;i<3;i++) {
452 for(int j=0;j<3;j++) {
453 matrix2[i][j]=matrix1[j][i];
454 }
455 }
456 //cout<<"matrix2:\n"<<matrix2<<endl;
457
458 HepMatrix matrix4(3,3);
459 matrix4=matrix1*matrix*matrix2;
460 //cout<<"matrix4:\n"<<matrix4<<endl;
461
462 HepSymMatrix matrix5(3); //error matrix of x, y, z
463 for(int i=0;i<3;i++) {
464 for(int j=0;j<=i;j++) {
465 matrix5[i][j]=matrix4[i][j];
466 }
467 }
468 //cout<<"matrix5:\n"<<matrix5<<endl;
469
470 aShower.setErrorMatrix(matrix5);
471
472 if(matrix5[0][0]>0) dx=sqrt(matrix5[0][0]);
473 if(matrix5[1][1]>0) dy=sqrt(matrix5[1][1]);
474 if(matrix5[2][2]>0) dz=sqrt(matrix5[2][2]);
475}
476
double sin(const BesAngle a)
Definition: BesAngle.h:210
double cos(const BesAngle a)
Definition: BesAngle.h:213
const double delta
double w
map< RecEmcID, RecEmcFraction, less< RecEmcID > > RecEmcFractionMap
void setDphi(double dpi)
Definition: DstEmcShower.h:66
double e3x3() const
Definition: DstEmcShower.h:48
void setDtheta(double dt)
Definition: DstEmcShower.h:65
double e5x5() const
Definition: DstEmcShower.h:49
void setPosition(const HepPoint3D &pos)
Definition: DstEmcShower.h:62
double energy() const
Definition: DstEmcShower.h:45
void setErrorMatrix(const HepSymMatrix &error)
Definition: DstEmcShower.h:75
static unsigned int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0)
Definition: EmcID.cxx:38
static unsigned int theta_module(const Identifier &id)
Definition: EmcID.cxx:43
static unsigned int phi_module(const Identifier &id)
Definition: EmcID.cxx:48
double LogPosOffset() const
double WestLogShMaxPhiPara(int n, int m) const
double BarrLogShMaxPhiPara(int n, int m) const
static EmcRecParameter & GetInstance()
double EastLogShMaxPhiPara(int n, int m) const
double EastLogShMaxThetaPara(int n, int m) const
double PosCorr() const
std::string PositionMode2() const
double SigTheta(int n) const
double WestLogShMaxThetaPara(int n, int m) const
double SigPhi(int n) const
double BarrLogShMaxThetaPara(int n, int m) const
virtual void Position(RecEmcShower &aShower)
virtual double GetBarrelR() const =0
virtual HepPoint3D GetCFrontCenter(const Identifier &id) const =0
virtual HepPoint3D GetCrystalPoint(const Identifier &id, const int i) const =0
virtual HepPoint3D GetCCenter(const Identifier &id) const =0
RecEmcFractionMap::const_iterator End() const
RecEmcFractionMap getFractionMap5x5() const
RecEmcID getShowerId() const
Definition: RecEmcShower.h:55
RecEmcFractionMap::const_iterator Begin() const
RecEmcEnergy getEAll() const
Definition: RecEmcShower.h:92
RecEmcFractionMap getFractionMap3x3() const