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

#include <EmcRecShowerPosLinShMax.h>

+ Inheritance diagram for EmcRecShowerPosLinShMax:

Public Member Functions

 EmcRecShowerPosLinShMax ()
 
virtual ~EmcRecShowerPosLinShMax ()
 
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 15 of file EmcRecShowerPosLinShMax.h.

Constructor & Destructor Documentation

◆ EmcRecShowerPosLinShMax()

EmcRecShowerPosLinShMax::EmcRecShowerPosLinShMax ( )
inline

Definition at line 18 of file EmcRecShowerPosLinShMax.h.

18{}

◆ ~EmcRecShowerPosLinShMax()

virtual EmcRecShowerPosLinShMax::~EmcRecShowerPosLinShMax ( )
inlinevirtual

Definition at line 19 of file EmcRecShowerPosLinShMax.h.

19{}

Member Function Documentation

◆ Position()

void EmcRecShowerPosLinShMax::Position ( RecEmcShower &  aShower)
virtual

Implements EmcRecShowerPosAbs.

Definition at line 18 of file EmcRecShowerPosLinShMax.cxx.

19{
20 RecEmcFractionMap::const_iterator cit;
21 HepPoint3D pos(0,0,0);
22 HepPoint3D posMax(0,0,0);
23 HepPoint3D possum(0,0,0);
24 double w,wsum=0;
25 double etot;
26
27//cout<<"EmcRecShowerPosLinShMax::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 for(cit=aShower.Begin();cit!=aShower.End();cit++){
39 etot+=(cit->second.getEnergy()*cit->second.getFraction());
40 //cout<<etot<<endl;
41 }
42 wsum=0;
43
44 double etota=aShower.getEAll();
45 double ltot=(log(etota/0.0145)+0.5)*1.86;
46 double e55=aShower.e5x5();
47 double len55=(log(etota/0.0145)+0.5)*1.86;
48 Identifier seedID = aShower.getShowerId();
49 HepPoint3D seedPoint(iGeoSvc->GetCFrontCenter(seedID));
50 //Rs is the distance from orgin to showerCentroid, unit cm
51// HepDouble Rs=seedPoint.mag();
52 for(cit=aShower.Begin();cit!=aShower.End();cit++){
53 w=(cit->second.getEnergy()*cit->second.getFraction());
54
55 HepPoint3D hitPoint(iGeoSvc->GetCFrontCenter(cit->second.getCellId()));
56// HepDouble Ang=seedPoint.angle(hitPoint);
57 // theDistance is the distance from orgin to showerMax
58 double theDistance=0;
59 if(cit->second.getCellId()==seedID) {
60// theDistance=(Rs+ltot);
61 theDistance=(seedPoint.mag()+ltot);
62 } else {
63 theDistance=(seedPoint.mag()+ltot)/cos(seedPoint.angle(hitPoint));
64// theDistance=(Rs+ltot)/cos(Ang);
65 }
66
67 //cout<<w<<endl;
68 wsum+=w;
69 pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
70
71 posMax = (theDistance/pos.mag())*pos;
72
73 possum+=posMax*w;
74 }
75
76 if(wsum<=0) {
77 cout<<"[[[[[[[[[[[[[[["<<wsum;
78 }
79 pos=possum/wsum;
80// pos=((pos.mag()-len55)/pos.mag())*pos;
81// aShower.setPosition(pos);
82 if(Para.PosCorr()==0){ aShower.setPosition(pos);}
83
84 //----------------------------------------------------------------------
85 //position correction
86 RecEmcID id = aShower.getShowerId();
87 const unsigned int module = EmcID::barrel_ec(id);;
88 const unsigned int thetaModule = EmcID::theta_module(id);
89 const unsigned int phiModule = EmcID::phi_module(id);
90 HepPoint3D posCorr(pos);
91
92 if(module==1) { //barrel
93 unsigned int theModule;
94 if(thetaModule>21) {
95 theModule = 43 - thetaModule;
96 id = EmcID::crystal_id(module,theModule,phiModule);
97 posCorr.setZ(-posCorr.z());
98 } else {
99 theModule = thetaModule;
100 }
101
102 double IRShift, parTheta[5],parPhi[5];
103 if(theModule==21) {
104 IRShift = 0;
105 } else if(theModule==20) {
106 IRShift = 2.5;
107 } else {
108 IRShift = 5.;
109 }
110
111 EmcRecParameter& Para=EmcRecParameter::GetInstance();
112
113 for(int i=0;i<5;i++){
114
115 if(e55>1.0) parTheta[i] = Para.BarrShLinThetaPara(theModule,i);
116 else if(e55<=1.0 && e55>0.5) parTheta[i] = Para.BarrShLinThetaPara(theModule+22,i);
117 else if(e55<=0.5) parTheta[i] = Para.BarrShLinThetaPara(theModule+44,i);
118
119
120 if(e55>1.0) parPhi[i] = Para.BarrShLinPhiPara(0,i);
121 else if(e55<=1.0 && e55>0.5) parPhi[i] = Para.BarrShLinPhiPara(1,i);
122 else if(e55<=0.5) parPhi[i] = Para.BarrShLinPhiPara(2,i);
123
124 }
125
126 HepPoint3D IR(0,0,IRShift);
127 HepPoint3D center01, center23; //center of point0,1 and point2,3 in crystal
128 center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
129 center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
130
131 double theta01,theta23,length,d;
132 theta01 = (center01-IR).theta();
133 theta23 = (center23-IR).theta();
134 length = (center01-IR).mag();
135 d = length*tan(theta01-theta23); //length in x direction
136
137 double xIn,xOut,deltaTheta;
138 xIn = length*tan(theta01-(posCorr-IR).theta())-d/2;
139 xOut = xIn-( parTheta[0]*atan(xIn*parTheta[1]-parTheta[4])
140 + parTheta[2]*xIn + parTheta[3] );
141 deltaTheta = atan((xOut+d/2)/length)-atan((xIn+d/2)/length);
142 //cout<<"xIn: "<<xIn<<" xOut: "<<xOut<<" deltaTheta: "<<deltaTheta<<endl;
143 //cout<<"dx: "<<parTheta[0]*atan(xIn*parTheta[1]-parTheta[4])
144 // + parTheta[2]*xIn + parTheta[3] <<endl;
145/////////////////////////////////////////////////////////////////////////////////////////////////
146 double yIn,yOut,deltaPhi;
147
148 if(phiModule==0&&posCorr.phi()<0){
149
150 yIn = posCorr.phi()*180./CLHEP::pi+360.-119*3-1.843;
151
152 }else if (phiModule==119&&posCorr.phi()>0){
153
154 yIn = posCorr.phi()*180./CLHEP::pi-1.843;
155
156 }else {
157
158 yIn = posCorr.phi() < 0 ? posCorr.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
159 : posCorr.phi()*180./CLHEP::pi-phiModule*3.-1.843;
160
161 }
162 yOut = parPhi[0]*atan(yIn*parPhi[1]-parPhi[4])+parPhi[2]*yIn+parPhi[3];
163
164// 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;
165
166 deltaPhi = yOut*CLHEP::pi/180.;
167
168 HepPoint3D rotateVector(0,1,0); //y axis
169 rotateVector.rotateZ(center01.phi());
170 posCorr.setZ(posCorr.z()-IRShift);
171 posCorr.rotate(-deltaTheta,rotateVector);
172 posCorr.setZ(posCorr.z()+IRShift);
173 //phi parameter is gotten by the average of e+ e- peak
174// posCorr.rotateZ(-0.002994);
175 posCorr.rotateZ(-deltaPhi);
176
177 if(thetaModule>21) {
178 posCorr.setZ(-posCorr.z());
179 }
180 }
181 posCorr=((posCorr.mag()-len55)/posCorr.mag())*posCorr;
182 // aShower.setPosition(posCorr);
183 if(Para.PosCorr()==1){ aShower.setPosition(posCorr);}
184
185 ////////////////////////////
186 if(aShower.energy()<1e-5) return;
187
188 double r,theta,phi;
189 double dtheta,dphi,dx,dy,dz;
190
191 r = posCorr.mag();
192 theta = posCorr.theta();
193 phi = posCorr.phi();
194 //cout<<"energy: "<<aShower.energy()<<" theta: "<<theta<<" phi: "<<phi<<endl;
195
196 if(aShower.energy()>0) {
197 dtheta = Para.SigTheta(0)/sqrt(aShower.energy())+Para.SigTheta(1);
198 dphi = Para.SigPhi(0)/sqrt(aShower.energy())+Para.SigPhi(1);
199 }
200 else {
201 dtheta = 0;
202 dphi = 0;
203 }
204 //cout<<"dtheta: "<<dtheta<<" dphi: "<<dphi<<endl;
205
206 dx = fabs(iGeoSvc->GetBarrelR()*sin(phi)*dphi);
207 dy = fabs(iGeoSvc->GetBarrelR()*cos(phi)*dphi);
208 if(theta>0)
209 dz = fabs(iGeoSvc->GetBarrelR()*dtheta/(sin(theta)*sin(theta)));
210 else
211 dz = 0;
212 //cout<<" dx: "<<dx<<" dy: "<<dy<<" dz: "<<dz<<endl;
213
214 aShower.setDtheta(dtheta);
215 aShower.setDphi(dphi);
216
217 //----------------------------------------------------------------------
218 // Error matrix
219 HepSymMatrix matrix(3); //error matrix of r, theta, phi
220 matrix[0][0]=0;
221 matrix[1][1]=dtheta*dtheta;
222 matrix[2][2]=dphi*dphi;
223 matrix[0][1]=0;
224 matrix[0][2]=0;
225 matrix[1][2]=0;
226 //cout<<"matrix:\n"<<matrix<<endl;
227
228 HepMatrix matrix1(3,3),matrix2(3,3);
229 matrix1[0][0]=sin(theta)*cos(phi);
230 matrix1[0][1]=r*cos(theta)*cos(phi);
231 matrix1[0][2]=-r*sin(theta)*sin(phi);
232 matrix1[1][0]=sin(theta)*sin(phi);
233 matrix1[1][1]=r*cos(theta)*sin(phi);
234 matrix1[1][2]=r*sin(theta)*cos(phi);
235 matrix1[2][0]=cos(theta);
236 matrix1[2][1]=-r*sin(theta);
237 matrix1[2][2]=0;
238 //cout<<"matrix1:\n"<<matrix1<<endl;
239
240 for(int i=0;i<3;i++) {
241 for(int j=0;j<3;j++) {
242 matrix2[i][j]=matrix1[j][i];
243 }
244 }
245 //cout<<"matrix2:\n"<<matrix2<<endl;
246
247 HepMatrix matrix4(3,3);
248 matrix4=matrix1*matrix*matrix2;
249 //cout<<"matrix4:\n"<<matrix4<<endl;
250
251 HepSymMatrix matrix5(3); //error matrix of x, y, z
252 for(int i=0;i<3;i++) {
253 for(int j=0;j<=i;j++) {
254 matrix5[i][j]=matrix4[i][j];
255 }
256 }
257 //cout<<"matrix5:\n"<<matrix5<<endl;
258
259 aShower.setErrorMatrix(matrix5);
260
261 if(matrix5[0][0]>0) dx=sqrt(matrix5[0][0]);
262 if(matrix5[1][1]>0) dy=sqrt(matrix5[1][1]);
263 if(matrix5[2][2]>0) dz=sqrt(matrix5[2][2]);
264}
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
w
double w
Definition: EvtPhokharaDef.hh:26
DstEmcShower::setDphi
void setDphi(double dpi)
Definition: DstEmcShower.h:66
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::BarrShLinPhiPara
double BarrShLinPhiPara(int n, int m) const
Definition: EmcRecParameter.cxx:711
EmcRecParameter::SigTheta
double SigTheta(int n) const
Definition: EmcRecParameter.cxx:629
EmcRecParameter::BarrShLinThetaPara
double BarrShLinThetaPara(int n, int m) const
Definition: EmcRecParameter.cxx:706
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::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
The documentation for this class was generated from the following files: