BOSS 7.0.4
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//
6#include "EmcRec/EmcRecShowerPosLogShMax.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
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_t etot
map< RecEmcID, RecEmcFraction, less< RecEmcID > > RecEmcFractionMap
double w
double sin(const BesAngle a)
double cos(const BesAngle a)
void setErrorMatrix(const HepSymMatrix &error)
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
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
RecEmcFractionMap::const_iterator Begin() const
RecEmcFractionMap getFractionMap3x3() const