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

#include <EmcRecShowerPosLogShMax.h>

+ Inheritance diagram for EmcRecShowerPosLogShMax:

Public Member Functions

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

Detailed Description

Definition at line 15 of file InstallArea/include/EmcRec/EmcRec/EmcRecShowerPosLogShMax.h.

Constructor & Destructor Documentation

◆ EmcRecShowerPosLogShMax() [1/2]

EmcRecShowerPosLogShMax::EmcRecShowerPosLogShMax ( )
inline

Definition at line 18 of file InstallArea/include/EmcRec/EmcRec/EmcRecShowerPosLogShMax.h.

18{}

◆ ~EmcRecShowerPosLogShMax() [1/2]

virtual EmcRecShowerPosLogShMax::~EmcRecShowerPosLogShMax ( )
inlinevirtual

Definition at line 19 of file InstallArea/include/EmcRec/EmcRec/EmcRecShowerPosLogShMax.h.

19{}

◆ EmcRecShowerPosLogShMax() [2/2]

EmcRecShowerPosLogShMax::EmcRecShowerPosLogShMax ( )
inline

Definition at line 18 of file Reconstruction/EmcRec/EmcRec-01-02-57/EmcRec/EmcRecShowerPosLogShMax.h.

18{}

◆ ~EmcRecShowerPosLogShMax() [2/2]

virtual EmcRecShowerPosLogShMax::~EmcRecShowerPosLogShMax ( )
inlinevirtual

Definition at line 19 of file Reconstruction/EmcRec/EmcRec-01-02-57/EmcRec/EmcRecShowerPosLogShMax.h.

19{}

Member Function Documentation

◆ Position() [1/2]

void EmcRecShowerPosLogShMax::Position ( RecEmcShower &  aShower)
virtual

Implements EmcRecShowerPosAbs.

Definition at line 17 of file EmcRecShowerPosLogShMax.cxx.

18{
19 RecEmcFractionMap::const_iterator cit;
20 HepPoint3D pos(0,0,0);
21 HepPoint3D posMax(0,0,0);
22 HepPoint3D possum(0,0,0);
23 double w,wsum=0;
24
25// cout<<"EmcRecShowerPosLogShMax::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 double offset=Para.LogPosOffset();
36 double e55 = aShower.e5x5();
37 double len55=(log(e55/0.0145)+0.5)*1.86;
38
39 if(Para.PositionMode2()=="all") {
40 double etot=aShower.getEAll();
41 double ltot=(log(etot/0.0145)+0.5)*1.86;
42 Identifier seedID = aShower.getShowerId();
43 HepPoint3D seedPoint(iGeoSvc->GetCFrontCenter(seedID));
44 //Rs is the distance from orgin to showerCentroid, unit cm
45// HepDouble Rs=seedPoint.mag();
46
47 for(cit=aShower.Begin();cit!=aShower.End();cit++){
48 w=offset+log((cit->second.getEnergy()/etot)*cit->second.getFraction());
49
50 HepPoint3D hitPoint(iGeoSvc->GetCFrontCenter(cit->second.getCellId()));
51// HepDouble Ang=seedPoint.angle(hitPoint);
52 // theDistance is the distance from orgin to showerMax
53 double theDistance;
54 if(cit->second.getCellId()==seedID) {
55// theDistance=(seedPoint.mag()+ltot);
56 theDistance=(seedPoint.mag()+ltot);
57
58 } else {
59// theDistance=(seedPoint.mag()+ltot)/cos(seedPoint.angle(hitPoint));
60 theDistance=(seedPoint.mag()+ltot)/cos(seedPoint.angle(hitPoint));
61 }
62 if(w>0){
63 wsum+=w;
64 pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
65 posMax = (theDistance/pos.mag())*pos;
66 possum+=posMax*w;
67 }
68 }
69 } else if(Para.PositionMode2()=="3x3") {
70 double e3x3=aShower.e3x3();
71 double l3x3=(log(e3x3/0.0145)+0.5)*1.86;
72 Identifier seedID = aShower.getShowerId();
73 HepPoint3D seedPoint(iGeoSvc->GetCFrontCenter(seedID));
74 //Rs is the distance from orgin to showerCentroid, unit cm
75// HepDouble Rs=seedPoint.mag();
76
77 RecEmcFractionMap fracMap3x3=aShower.getFractionMap3x3();
78 for(cit=fracMap3x3.begin();
79 cit!=fracMap3x3.end();
80 cit++) {
81 w=offset+log((cit->second.getEnergy()/e3x3)*cit->second.getFraction());
82 HepPoint3D hitPoint(iGeoSvc->GetCFrontCenter(cit->second.getCellId()));
83// HepDouble Ang=seedPoint.angle(hitPoint);
84 // theDistance is the distance from orgin to showerMax
85 double theDistance;
86 if(cit->second.getCellId()==seedID) {
87 theDistance=(seedPoint.mag()+l3x3);
88 } else {
89 theDistance=(seedPoint.mag()+l3x3)/cos(seedPoint.angle(hitPoint));
90 }
91 if(w>0){
92 wsum+=w;
93 pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
94 posMax = (theDistance/pos.mag())*pos;
95 possum+=posMax*w;
96 }
97 }
98 } else {
99 double e5x5=aShower.e5x5();
100 double l5x5=(log(e5x5/0.0145)+0.5)*1.86;
101 Identifier seedID = aShower.getShowerId();
102 HepPoint3D seedPoint(iGeoSvc->GetCFrontCenter(seedID));
103 //Rs is the distance from orgin to showerCentroid, unit cm
104// HepDouble Rs=seedPoint.mag();
105
106 RecEmcFractionMap fracMap5x5=aShower.getFractionMap5x5();
107 for(cit=fracMap5x5.begin();
108 cit!=fracMap5x5.end();
109 cit++) {
110 w=offset+log((cit->second.getEnergy()/e5x5)*cit->second.getFraction());
111
112 HepPoint3D hitPoint(iGeoSvc->GetCFrontCenter(cit->second.getCellId()));
113// HepDouble Ang=seedPoint.angle(hitPoint);
114 // theDistance is the distance from orgin to showerMax
115 double theDistance;
116 if(cit->second.getCellId()==seedID) {
117 theDistance=(seedPoint.mag()+l5x5);
118 } else {
119 theDistance=(seedPoint.mag()+l5x5)/cos(seedPoint.angle(hitPoint));
120 }
121
122 if(w>0){
123 wsum+=w;
124 pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
125 posMax = (theDistance/pos.mag())*pos;
126 possum+=posMax*w;
127 //cout<<"posMax="<<posMax<<endl;
128 }
129 }
130 }
131
132 if(wsum<=0) {
133 cout<<"[[[[[[[[[[[[[[["<<wsum;
134 }
135 pos=possum/wsum;
136
137// pos=((pos.mag()-len55)/pos.mag())*pos;
138// aShower.setPosition(pos);
139
140 if(Para.PosCorr()==0){ aShower.setPosition(pos);}
141
142 //----------------------------------------------------------------------
143 //position correction
144 RecEmcID id = aShower.getShowerId();
145 const unsigned int module = EmcID::barrel_ec(id);;
146 const unsigned int thetaModule = EmcID::theta_module(id);
147 const unsigned int phiModule = EmcID::phi_module(id);
148 HepPoint3D posCorr(pos);
149
150 if(module==1) { //barrel
151 unsigned int theModule;
152 if(thetaModule>21) {
153 theModule = 43 - thetaModule;
154 id = EmcID::crystal_id(module,theModule,phiModule);
155 posCorr.setZ(-posCorr.z());
156 } else {
157 theModule = thetaModule;
158 }
159
160 double IRShift, parTheta[5],parPhi[5];
161 if(theModule==21) {
162 IRShift = 0;
163 } else if(theModule==20) {
164 IRShift = 2.5;
165 } else {
166 IRShift = 5.;
167 }
168
169 EmcRecParameter& Para=EmcRecParameter::GetInstance();
170
171 for(int i=0;i<5;i++){
172
173 if(e55>1.0) parTheta[i] = Para.BarrShLogThetaPara(theModule,i);
174 else if(e55<=1.0 && e55>0.5) parTheta[i] = Para.BarrShLogThetaPara(theModule+22,i);
175 else if(e55<=0.5) parTheta[i] = Para.BarrShLogThetaPara(theModule+44,i);
176
177 if(e55>1.0) parPhi[i] = Para.BarrShLogPhiPara(0,i);
178 else if(e55<=1.0 && e55>0.5) parPhi[i] = Para.BarrShLogPhiPara(1,i);
179 else if(e55<=0.5) parPhi[i] = Para.BarrShLogPhiPara(2,i);
180 }
181
182 HepPoint3D IR(0,0,IRShift);
183 HepPoint3D center01, center23; //center of point0,1 and point2,3 in crystal
184 center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
185 center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
186
187 double theta01,theta23,length,d;
188 theta01 = (center01-IR).theta();
189 theta23 = (center23-IR).theta();
190 length = (center01-IR).mag();
191 d = length*tan(theta01-theta23); //length in x direction
192
193 double xIn,xOut,deltaTheta;
194 xIn = length*tan(theta01-(posCorr-IR).theta())-d/2;
195 xOut = xIn-( parTheta[0]*atan(xIn*parTheta[1]-parTheta[4])
196 + parTheta[2]*xIn + parTheta[3] );
197 deltaTheta = atan((xOut+d/2)/length)-atan((xIn+d/2)/length);
198 //cout<<"xIn: "<<xIn<<" xOut: "<<xOut<<" deltaTheta: "<<deltaTheta<<endl;
199 //cout<<"dx: "<<parTheta[0]*atan(xIn*parTheta[1]-parTheta[4])
200 // + parTheta[2]*xIn + parTheta[3] <<endl;
201/////////////////////////////////////////////////////////////////////////////////////////////////
202 double yIn,yOut,deltaPhi;
203
204 if(phiModule==0&&posCorr.phi()<0){
205
206 yIn = posCorr.phi()*180./CLHEP::pi+360.-119*3-1.843;
207
208 }else if (phiModule==119&&posCorr.phi()>0){
209
210 yIn = posCorr.phi()*180./CLHEP::pi-1.843;
211
212 }else {
213
214 yIn = posCorr.phi() < 0 ? posCorr.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
215 : posCorr.phi()*180./CLHEP::pi-phiModule*3.-1.843;
216
217 }
218 yOut = parPhi[0]*atan(yIn*parPhi[1]-parPhi[4])+parPhi[2]*yIn+parPhi[3];
219
220// yOut = parPhi[0]+parPhi[1]*yIn+parPhi[2]*yIn*yIn+parPhi[3]*yIn*yIn*yIn+parPhi[4]*yIn*yIn*yIn*yIn+parPhi[5]*yIn*yIn*yIn*yIn*yIn;
221 deltaPhi = yOut*CLHEP::pi/180.;
222
223 HepPoint3D rotateVector(0,1,0); //y axis
224 rotateVector.rotateZ(center01.phi());
225 posCorr.setZ(posCorr.z()-IRShift);
226 posCorr.rotate(-deltaTheta,rotateVector);
227 posCorr.setZ(posCorr.z()+IRShift);
228// posCorr.rotateZ(-0.002994);
229 posCorr.rotateZ(-deltaPhi);
230
231 if(thetaModule>21) {
232 posCorr.setZ(-posCorr.z());
233 }
234 }
235 posCorr=((posCorr.mag()-len55)/posCorr.mag())*posCorr;
236 // aShower.setPosition(posCorr);
237 if(Para.PosCorr()==1){ aShower.setPosition(posCorr);}
238
239 ////////////////////////////
240 if(aShower.energy()<1e-5) return;
241
242 double r,theta,phi;
243 double dtheta,dphi,dx,dy,dz;
244
245 r = posCorr.mag();
246 theta = posCorr.theta();
247 phi = posCorr.phi();
248 //cout<<"energy: "<<aShower.energy()<<" theta: "<<theta<<" phi: "<<phi<<endl;
249
250 if(aShower.energy()>0) {
251 dtheta = Para.SigTheta(0)/sqrt(aShower.energy())+Para.SigTheta(1);
252 dphi = Para.SigPhi(0)/sqrt(aShower.energy())+Para.SigPhi(1);
253 }
254 else {
255 dtheta = 0;
256 dphi = 0;
257 }
258 //cout<<"dtheta: "<<dtheta<<" dphi: "<<dphi<<endl;
259
260 dx = fabs(iGeoSvc->GetBarrelR()*sin(phi)*dphi);
261 dy = fabs(iGeoSvc->GetBarrelR()*cos(phi)*dphi);
262 if(theta>0)
263 dz = fabs(iGeoSvc->GetBarrelR()*dtheta/(sin(theta)*sin(theta)));
264 else
265 dz = 0;
266 //cout<<" dx: "<<dx<<" dy: "<<dy<<" dz: "<<dz<<endl;
267
268 aShower.setDtheta(dtheta);
269 aShower.setDphi(dphi);
270
271 //----------------------------------------------------------------------
272 // Error matrix
273 HepSymMatrix matrix(3); //error matrix of r, theta, phi
274 matrix[0][0]=0;
275 matrix[1][1]=dtheta*dtheta;
276 matrix[2][2]=dphi*dphi;
277 matrix[0][1]=0;
278 matrix[0][2]=0;
279 matrix[1][2]=0;
280 //cout<<"matrix:\n"<<matrix<<endl;
281
282 HepMatrix matrix1(3,3),matrix2(3,3);
283 matrix1[0][0]=sin(theta)*cos(phi);
284 matrix1[0][1]=r*cos(theta)*cos(phi);
285 matrix1[0][2]=-r*sin(theta)*sin(phi);
286 matrix1[1][0]=sin(theta)*sin(phi);
287 matrix1[1][1]=r*cos(theta)*sin(phi);
288 matrix1[1][2]=r*sin(theta)*cos(phi);
289 matrix1[2][0]=cos(theta);
290 matrix1[2][1]=-r*sin(theta);
291 matrix1[2][2]=0;
292 //cout<<"matrix1:\n"<<matrix1<<endl;
293
294 for(int i=0;i<3;i++) {
295 for(int j=0;j<3;j++) {
296 matrix2[i][j]=matrix1[j][i];
297 }
298 }
299 //cout<<"matrix2:\n"<<matrix2<<endl;
300
301 HepMatrix matrix4(3,3);
302 matrix4=matrix1*matrix*matrix2;
303 //cout<<"matrix4:\n"<<matrix4<<endl;
304
305 HepSymMatrix matrix5(3); //error matrix of x, y, z
306 for(int i=0;i<3;i++) {
307 for(int j=0;j<=i;j++) {
308 matrix5[i][j]=matrix4[i][j];
309 }
310 }
311 //cout<<"matrix5:\n"<<matrix5<<endl;
312
313 aShower.setErrorMatrix(matrix5);
314
315 if(matrix5[0][0]>0) dx=sqrt(matrix5[0][0]);
316 if(matrix5[1][1]>0) dy=sqrt(matrix5[1][1]);
317 if(matrix5[2][2]>0) dz=sqrt(matrix5[2][2]);
318}
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
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::GetInstance
static EmcRecParameter & GetInstance()
Definition: EmcRecParameter.cxx:461
EmcRecParameter::PosCorr
double PosCorr() const
Definition: EmcRecParameter.cxx:519
EmcRecParameter::SigTheta
double SigTheta(int n) const
Definition: EmcRecParameter.cxx:553
EmcRecParameter::BarrShLogThetaPara
double BarrShLogThetaPara(int n, int m) const
Definition: EmcRecParameter.cxx:620
EmcRecParameter::BarrShLogPhiPara
double BarrShLogPhiPara(int n, int m) const
Definition: EmcRecParameter.cxx:625
EmcRecParameter::SigPhi
double SigPhi(int n) const
Definition: EmcRecParameter.cxx:558
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

◆ Position() [2/2]

virtual void EmcRecShowerPosLogShMax::Position ( RecEmcShower &  aShower)
virtual

Implements EmcRecShowerPosAbs.

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