BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
Public Member Functions | List of all members
EmcRecShowerPosLoglin Class Reference

#include <EmcRecShowerPosLoglin.h>

+ Inheritance diagram for EmcRecShowerPosLoglin:

Public Member Functions

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

Detailed Description

Definition at line 17 of file EmcRecShowerPosLoglin.h.

Constructor & Destructor Documentation

◆ EmcRecShowerPosLoglin()

EmcRecShowerPosLoglin::EmcRecShowerPosLoglin ( )
inline

Definition at line 21 of file EmcRecShowerPosLoglin.h.

21{}

◆ ~EmcRecShowerPosLoglin()

virtual EmcRecShowerPosLoglin::~EmcRecShowerPosLoglin ( )
inlinevirtual

Definition at line 22 of file EmcRecShowerPosLoglin.h.

22{}

Member Function Documentation

◆ Position()

void EmcRecShowerPosLoglin::Position ( RecEmcShower & aShower)
virtual

Implements EmcRecShowerPosAbs.

Definition at line 19 of file EmcRecShowerPosLoglin.cxx.

20{
21 RecEmcFractionMap::const_iterator cit;
22 HepPoint3D pos(0,0,0);
23 HepPoint3D possum(0,0,0);
24 double w,w1,w2,wsum=0;
25 double num,num0;
26
27 // cout<<"EmcRecShowerPosLoglin::Position Begin"<<endl;
28
29 IEmcRecGeoSvc* iGeoSvc;
30 ISvcLocator* svcLocator = Gaudi::svcLocator();
31 StatusCode sc = svcLocator->service("EmcRecGeoSvc",iGeoSvc);
32 if(sc!=StatusCode::SUCCESS) {
33 // cout<<"Error: Can't get EmcRecGeoSvc"<<endl;
34 }
35
36 EmcRecParameter& Para=EmcRecParameter::GetInstance();
37
38
39 // double offset=Para.LogPosOffset();
40 double a0=4;
41 double offset;
42 //cout<<offset<<endl;
43 double e5x5 = -99999;
44 if(Para.PositionMode2()=="all") {
45 double etot=aShower.getEAll();
46 for(cit=aShower.Begin();cit!=aShower.End();cit++){
47 w=offset+log((cit->second.getEnergy()/etot)*cit->second.getFraction());
48 if(w>0){
49 wsum+=w;
50 pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
51 possum+=pos*w;
52 }
53 }
54 } else if(Para.PositionMode2()=="3x3") {
55 double e3x3=aShower.e3x3();
56 RecEmcFractionMap fracMap3x3=aShower.getFractionMap3x3();
57 for(cit=fracMap3x3.begin();cit!=fracMap3x3.end();cit++) {
58 w=offset+log((cit->second.getEnergy()/e3x3)*cit->second.getFraction());
59 if(w>0){
60 wsum+=w;
61 pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
62 possum+=pos*w;
63 }
64 }
65 } else {
66 e5x5=aShower.e5x5();
67 offset=a0-1.594*exp(-2.543*e5x5);
68 RecEmcFractionMap fracMap5x5=aShower.getFractionMap5x5();
69
70 num0=0;
71 num=0;
72 for(cit=fracMap5x5.begin();cit!=fracMap5x5.end();cit++) {
73 w1=offset+log((cit->second.getEnergy()/e5x5)*cit->second.getFraction());
74 w2=(cit->second.getEnergy()/e5x5)*cit->second.getFraction();
75
76 num0++;
77
78 if(w1>0){
79 num++;
80 }
81
82 }
83
84 // cout<<"num0=="<<num0<<endl;
85 // cout<<"num==="<<num<<endl;
86 for(cit=fracMap5x5.begin();cit!=fracMap5x5.end();cit++) {
87 w1=offset+log((cit->second.getEnergy()/e5x5)*cit->second.getFraction());
88 w2=(cit->second.getEnergy()/e5x5)*cit->second.getFraction();
89
90 num0++;
91
92 if(w1>0){
93 if(num>1){
94 // cout<<"w1 be used"<<endl;
95 wsum+=w1;
96 pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
97 possum+=pos*w1;
98 } else{
99 // cout<<"w2 be used"<<endl;
100 wsum+=w2;
101 pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
102 possum+=pos*w2;
103 }
104 }
105 }
106
107
108 }
109 if(wsum<=0) {
110 // cout<<"[[[[[[[[[[[[[[["<<wsum<<endl;
111 }
112 pos=possum/wsum;
113
114 //PosCorr=0 without position correction
115 //PosCorr=1 with position correction
116 if(Para.PosCorr()==0){
117
118 aShower.setPosition(pos);
119 }
120
121 //----------------------------------------------------------------------
122 //position correction
123 RecEmcID id = aShower.getShowerId();
124 const unsigned int module = EmcID::barrel_ec(id);;
125 const unsigned int thetaModule = EmcID::theta_module(id);
126 const unsigned int phiModule = EmcID::phi_module(id);
127 HepPoint3D posCorr(pos);
128
129 if(module==1) { //barrel
130 unsigned int theModule;
131 if(thetaModule>21) {
132 theModule = 43 - thetaModule;
133 id = EmcID::crystal_id(module,theModule,phiModule);
134 posCorr.setZ(-posCorr.z());
135 } else {
136 theModule = thetaModule;
137 }
138
139 double IRShift=0, parTheta[5],parPhi[5];
140 if(theModule==21) {
141 IRShift = 0;
142 } else if(theModule==20) {
143 IRShift = 2.5;
144 } else {
145 IRShift = 5.;
146 }
147
148
149 EmcRecParameter& Para=EmcRecParameter::GetInstance();
150 for(int i=0;i<5;i++){
151
152 parTheta[i] = Para.BarrLoglinThetaPara(theModule,i);
153
154 parPhi[i] = Para.BarrLoglinPhiPara(0,i);
155
156// cout<<"Para.BarrLoglinThetaPara"<<i<<"===="<<Para.BarrLoglinThetaPara(theModule,i)<<endl;
157// cout<<"Para.BarrLoglinPhiPara"<<i<<"===="<<Para.BarrLoglinPhiPara(0,i)<<endl;
158
159// cout<<"theModule===="<<theModule<<endl;
160// cout<<"partheta"<<i<<"===="<<parTheta[i]<<endl;
161// cout<<"parphi"<<i<<"===="<<parPhi[i]<<endl;
162
163
164 }
165
166 HepPoint3D IR(0,0,IRShift);
167 HepPoint3D center01, center23; //center of point0,1 and point2,3 in crystal
168 center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
169 center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
170
171 double theta01,theta23,length,d;
172 theta01 = (center01-IR).theta();
173 theta23 = (center23-IR).theta();
174 length = (center01-IR).mag();
175 d = length*tan(theta01-theta23); //length in x direction
176
177 double xIn,xOut,deltaTheta;
178 xIn = length*tan(theta01-(posCorr-IR).theta())-d/2;
179 xOut = xIn-( parTheta[0]*atan(xIn*parTheta[1]-parTheta[4])+ parTheta[2]*xIn + parTheta[3] );
180 deltaTheta = atan((xOut+d/2)/length)-atan((xIn+d/2)/length);
181 //cout<<"xIn: "<<xIn<<" xOut: "<<xOut<<" deltaTheta: "<<deltaTheta<<endl;
182 //cout<<"dx: "<<parTheta[0]*atan(xIn*parTheta[1]-parTheta[4])
183 // + parTheta[2]*xIn + parTheta[3] <<endl;
184
185 //obtained by photon, not used, just for comparison
186 // double parPhi[5] = { 51.1, -0.1499, 7.62, 0.1301, 0.005581 };
187 double yIn,yOut,deltaPhi;
188 // yIn = posCorr.phi() < 0 ? posCorr.phi()*180./CLHEP::pi+360.-phiModule*3.-1.1537
189 // : posCorr.phi()*180./CLHEP::pi-phiModule*3.-1.1537;
190 // yOut = parPhi[0]*atan(yIn*parPhi[1]-parPhi[4])+parPhi[2]*yIn+parPhi[3];
191
192 yIn = posCorr.phi() < 0 ? posCorr.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
193 : posCorr.phi()*180./CLHEP::pi-phiModule*3.-1.843;
194
195 if(phiModule==0&&posCorr.phi()<0){
196 yIn = posCorr.phi()*180./CLHEP::pi+360.-119*3-1.843;
197 }else if (phiModule==119&&posCorr.phi()>0){
198 yIn = posCorr.phi()*180./CLHEP::pi-1.843;
199 }else {
200 yIn = posCorr.phi() < 0 ? posCorr.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843 : posCorr.phi()*180./CLHEP::pi-phiModule*3.-1.843;
201 }
202 yOut = parPhi[0]*atan(yIn*parPhi[1]-parPhi[4])+parPhi[2]*yIn+parPhi[3];
203 deltaPhi = yOut*CLHEP::pi/180.;
204
205 HepPoint3D rotateVector(0,1,0); //y axis
206 rotateVector.rotateZ(center01.phi());
207 posCorr.setZ(posCorr.z()-IRShift);
208 posCorr.rotate(-deltaTheta,rotateVector);
209 posCorr.setZ(posCorr.z()+IRShift);
210 //phi parameter is gotten by the average of e+ e- peak
211
212 posCorr.rotateZ(-deltaPhi);
213
214 if(thetaModule>21) {
215 posCorr.setZ(-posCorr.z());
216 }
217 }
218 if(Para.PosCorr()==1){
219
220 aShower.setPosition(posCorr);
221 }
222 ////////////////////////////
223 if(aShower.energy()<1e-5) return;
224
225 double r,theta,phi;
226 double dtheta,dphi,dx,dy,dz;
227
228 r = posCorr.mag();
229 theta = posCorr.theta();
230 phi = posCorr.phi();
231 //cout<<"energy: "<<aShower.energy()<<" theta: "<<theta<<" phi: "<<phi<<endl;
232
233 if(aShower.energy()>0) {
234 dtheta = Para.SigTheta(0)/sqrt(aShower.energy())+Para.SigTheta(1);
235 dphi = Para.SigPhi(0)/sqrt(aShower.energy())+Para.SigPhi(1);
236 }
237 else {
238 dtheta = 0;
239 dphi = 0;
240 }
241 //cout<<"dtheta: "<<dtheta<<" dphi: "<<dphi<<endl;
242
243 dx = fabs(iGeoSvc->GetBarrelR()*sin(phi)*dphi);
244 dy = fabs(iGeoSvc->GetBarrelR()*cos(phi)*dphi);
245 if(theta>0)
246 dz = fabs(iGeoSvc->GetBarrelR()*dtheta/(sin(theta)*sin(theta)));
247 else
248 dz = 0;
249 //cout<<" dx: "<<dx<<" dy: "<<dy<<" dz: "<<dz<<endl;
250
251 aShower.setDtheta(dtheta);
252 aShower.setDphi(dphi);
253
254 //----------------------------------------------------------------------
255 // Error matrix
256 HepSymMatrix matrix(3); //error matrix of r, theta, phi
257 matrix[0][0]=0;
258 matrix[1][1]=dtheta*dtheta;
259 matrix[2][2]=dphi*dphi;
260 matrix[0][1]=0;
261 matrix[0][2]=0;
262 matrix[1][2]=0;
263 //cout<<"matrix:\n"<<matrix<<endl;
264
265 HepMatrix matrix1(3,3),matrix2(3,3);
266 matrix1[0][0]=sin(theta)*cos(phi);
267 matrix1[0][1]=r*cos(theta)*cos(phi);
268 matrix1[0][2]=-r*sin(theta)*sin(phi);
269 matrix1[1][0]=sin(theta)*sin(phi);
270 matrix1[1][1]=r*cos(theta)*sin(phi);
271 matrix1[1][2]=r*sin(theta)*cos(phi);
272 matrix1[2][0]=cos(theta);
273 matrix1[2][1]=-r*sin(theta);
274 matrix1[2][2]=0;
275 //cout<<"matrix1:\n"<<matrix1<<endl;
276
277 for(int i=0;i<3;i++) {
278 for(int j=0;j<3;j++) {
279 matrix2[i][j]=matrix1[j][i];
280 }
281 }
282 //cout<<"matrix2:\n"<<matrix2<<endl;
283
284 HepMatrix matrix4(3,3);
285 matrix4=matrix1*matrix*matrix2;
286 //cout<<"matrix4:\n"<<matrix4<<endl;
287
288 HepSymMatrix matrix5(3); //error matrix of x, y, z
289 for(int i=0;i<3;i++) {
290 for(int j=0;j<=i;j++) {
291 matrix5[i][j]=matrix4[i][j];
292 }
293 }
294 //cout<<"matrix5:\n"<<matrix5<<endl;
295
296 aShower.setErrorMatrix(matrix5);
297
298 if(matrix5[0][0]>0) dx=sqrt(matrix5[0][0]);
299 if(matrix5[1][1]>0) dy=sqrt(matrix5[1][1]);
300 if(matrix5[2][2]>0) dz=sqrt(matrix5[2][2]);
301}
tan
double tan(const BesAngle a)
Definition BesAngle.h:216
sin
double sin(const BesAngle a)
Definition BesAngle.h:210
cos
double cos(const BesAngle a)
Definition BesAngle.h:213
etot
Double_t etot
Definition DataBase/tau_mode.c:7
a0
character *LEPTONflag integer iresonances real zeta5 real a0
Definition Ekhara/Ekhara-00-01-02/src/ekhara/common.h:4
exp
EvtComplex exp(const EvtComplex &c)
Definition EvtComplex.hh:252
w
double w
Definition EvtPhokharaDef.hh:29
RecEmcFractionMap
map< RecEmcID, RecEmcFraction, less< RecEmcID > > RecEmcFractionMap
Definition RecEmcFraction.h:56
DstEmcShower::setDphi
void setDphi(double dpi)
Definition DstEmcShower.h:66
DstEmcShower::e3x3
double e3x3() const
Definition DstEmcShower.h:48
DstEmcShower::setDtheta
void setDtheta(double dt)
Definition DstEmcShower.h:65
DstEmcShower::e5x5
double e5x5() const
Definition DstEmcShower.h:49
DstEmcShower::setPosition
void setPosition(const HepPoint3D &pos)
Definition DstEmcShower.h:62
DstEmcShower::energy
double energy() const
Definition DstEmcShower.h:45
DstEmcShower::setErrorMatrix
void setErrorMatrix(const HepSymMatrix &error)
Definition 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::GetInstance
static EmcRecParameter & GetInstance()
Definition EmcRecParameter.cxx:537
EmcRecParameter::PosCorr
double PosCorr() const
Definition EmcRecParameter.cxx:595
EmcRecParameter::PositionMode2
std::string PositionMode2() const
Definition EmcRecParameter.h:194
EmcRecParameter::SigTheta
double SigTheta(int n) const
Definition EmcRecParameter.cxx:629
EmcRecParameter::BarrLoglinThetaPara
double BarrLoglinThetaPara(int n, int m) const
Definition EmcRecParameter.cxx:672
EmcRecParameter::BarrLoglinPhiPara
double BarrLoglinPhiPara(int n, int m) const
Definition EmcRecParameter.cxx:679
EmcRecParameter::SigPhi
double SigPhi(int n) const
Definition EmcRecParameter.cxx:634
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::getShowerId
RecEmcID getShowerId() const
Definition RecEmcShower.h:55
RecEmcShower::Begin
RecEmcFractionMap::const_iterator Begin() const
Definition RecEmcShower.cxx:131
RecEmcShower::getEAll
RecEmcEnergy getEAll() const
Definition RecEmcShower.h:92
RecEmcShower::getFractionMap3x3
RecEmcFractionMap getFractionMap3x3() const
Definition RecEmcShower.cxx:229
The documentation for this class was generated from the following files: