BOSS 7.1.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EmcRecShowerPosLog Class Reference

#include <EmcRecShowerPosLog.h>

+ Inheritance diagram for EmcRecShowerPosLog:

Public Member Functions

 EmcRecShowerPosLog ()
 
virtual ~EmcRecShowerPosLog ()
 
virtual void Position (RecEmcShower &aShower)
 
- Public Member Functions inherited from EmcRecShowerPosAbs
 EmcRecShowerPosAbs ()
 
virtual ~EmcRecShowerPosAbs ()
 
virtual void Position (RecEmcShower &aShower)=0
 

Detailed Description

Definition at line 17 of file EmcRecShowerPosLog.h.

Constructor & Destructor Documentation

◆ EmcRecShowerPosLog()

EmcRecShowerPosLog::EmcRecShowerPosLog ( )
inline

Definition at line 21 of file EmcRecShowerPosLog.h.

21{}

◆ ~EmcRecShowerPosLog()

virtual EmcRecShowerPosLog::~EmcRecShowerPosLog ( )
inlinevirtual

Definition at line 22 of file EmcRecShowerPosLog.h.

22{}

Member Function Documentation

◆ Position()

void EmcRecShowerPosLog::Position ( RecEmcShower aShower)
virtual

Implements EmcRecShowerPosAbs.

Definition at line 18 of file EmcRecShowerPosLog.cxx.

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 //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}
double tan(const BesAngle a)
Definition: BesAngle.h:216
double sin(const BesAngle a)
Definition: BesAngle.h:210
double cos(const BesAngle a)
Definition: BesAngle.h:213
const double delta
Double_t etot
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 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
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 BarrPosMCCor(int ntheta, int nphi) const
static EmcRecParameter & GetInstance()
double EastPosDataCor(int ntheta, int nphi) const
double WestPosDataCor(int ntheta, int nphi) const
double BarrLogPhiPara(int n, int m) const
double PosCorr() const
double MethodMode() const
std::string PositionMode2() const
double SigTheta(int n) const
double WestLogThetaPara(int n, int m) const
double BarrPosDataCor(int ntheta, int nphi) const
double BarrLogThetaPara(int n, int m) const
double WestLogPhiPara(int n, int m) const
double EastLogPhiPara(int n, int m) const
double EastLogThetaPara(int n, int m) const
double DataMode() const
double SigPhi(int n) const
virtual double GetBarrelR() const =0
virtual HepPoint3D GetCFrontCenter(const Identifier &id) const =0
virtual HepPoint3D GetCrystalPoint(const Identifier &id, const int i) 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

The documentation for this class was generated from the following files: