CGEM BOSS 6.6.5.f
BESIII Offline Software System
Loading...
Searching...
No Matches
BesEmcGeometry.cc
Go to the documentation of this file.
1//---------------------------------------------------------------------------//
2// BESIII Object_Oreiented Simulation and Reconstruction Tool //
3//---------------------------------------------------------------------------//
4//Descpirtion: Geometry of EMC detector
5//Author: Fu Chengdong
6//Created: Oct 23, 2003
7//Comment:
8//---------------------------------------------------------------------------//
9//
10#include "G4ThreeVector.hh"
11
12#include "BesEmcGeometry.hh"
13#include "BesEmcParameter.hh"
14#include "G4ios.hh"
15
16BesEmcGeometry::BesEmcGeometry()
17{
18 verboseLevel = 0;
19}
20
21BesEmcGeometry::~BesEmcGeometry()
22{;}
23
24void BesEmcGeometry::ReadEMCParameters()
25{
26 // Read EMC parameters from database
27 BesEmcParameter& emcPara=BesEmcParameter::GetInstance();
28
29 TaperRingDz = emcPara.GetTaperRingDz();
30 TaperRingThickness1 = emcPara.GetTaperRingThickness1();
31 TaperRingThickness2 = emcPara.GetTaperRingThickness2();
32 TaperRingThickness3 = emcPara.GetTaperRingThickness3();
33 TaperRingTheta = emcPara.GetTaperRingTheta()*deg;
34 TaperRingInnerLength = emcPara.GetTaperRingInnerLength();
35 TaperRingOuterLength = emcPara.GetTaperRingOuterLength();
36
37 BSCRmin = emcPara.GetBSCRmin();
38 BSCDz = emcPara.GetBSCDz()-TaperRingThickness3;
39
40 BSCRmin1 = emcPara.GetBSCRmin1();
41 BSCRmax1 = emcPara.GetBSCRmax1();
42 BSCRmin2 = emcPara.GetBSCRmin2();
43 //BSCRmax2 = emcPara.GetBSCRmax2();
44 BSCDz1 = emcPara.GetBSCDz1()-TaperRingThickness3;
45
46 BSCAngleRotat = emcPara.GetBSCAngleRotat()*deg;
47 BSCNbPhi = emcPara.GetBSCNbPhi();
48 BSCNbTheta = emcPara.GetBSCNbTheta();
49
50 BSCCryLength = emcPara.GetCrystalLength();
51 BSCCryLength1 = emcPara.GetCrystalLength1();
52 BSCYFront0 = emcPara.GetBSCYFront0();
53 BSCYFront = emcPara.GetBSCYFront();
54 BSCYFront1 = emcPara.GetBSCYFront1();
55 BSCPosition0 = emcPara.GetBSCPosition0();
56 BSCPosition1 = emcPara.GetBSCPosition1();
57
58 fTyvekThickness = emcPara.GetTyvekThickness();
59 fAlThickness = emcPara.GetAlThickness();
60 fMylarThickness = emcPara.GetMylarThickness();
61
62 rearBoxLength = emcPara.GetRearBoxLength();
63 rearBoxDz = emcPara.GetRearBoxDz();
64
65 HangingPlateDz = emcPara.GetHangingPlateDz();
66 OCGirderAngle = emcPara.GetOCGirderAngle()*deg;
67
68 rearCasingThickness = emcPara.GetRearCasingThickness();
69
70 orgGlassLengthX = emcPara.GetOrgGlassLengthX();
71 orgGlassLengthY = emcPara.GetOrgGlassLengthY();
72 orgGlassLengthZ = emcPara.GetOrgGlassLengthZ();
73
74 PDLengthX = emcPara.GetPDLengthX();
75 PDLengthY = emcPara.GetPDLengthY();
76 PDLengthZ = emcPara.GetPDLengthZ();
77
78 AlPlateDz = emcPara.GetAlPlateDz();
79 PABoxDz = emcPara.GetPABoxDz();
80 PABoxThickness = emcPara.GetPABoxThickness();
81
82 cableDr = emcPara.GetCableDr();
83 waterPipeDr = emcPara.GetWaterPipeDr();
84 waterPipeThickness= emcPara.GetWaterPipeThickness();
85
86 SPBarThickness = emcPara.GetSPBarThickness();
87 SPBarThickness1 = emcPara.GetSPBarThickness1();
88 SPBarwidth = emcPara.GetSPBarwidth();
89
90 EndRingDz = emcPara.GetEndRingDz();
91 EndRingDr = emcPara.GetEndRingDr();
92 EndRingRmin = emcPara.GetEndRingRmin();
93
94 for(G4int i=0;i<BSCNbTheta*6;i++)
95 {
96 zHalfLength[i] =0;
97 thetaAxis[i] =0;
98 phiAxis [i] =0;
99 yHalfLength1[i] =0;
100 xHalfLength1[i] =0;
101 xHalfLength2[i] =0;
102 tanAlpha1[i] =0;
103 yHalfLength2[i] =0;
104 xHalfLength3[i] =0;
105 xHalfLength4[i] =0;
106 tanAlpha2[i] =0;
107 thetaPosition[i]=0;
108 xPosition[i] =0;
109 yPosition[i] =0;
110 zPosition[i] =0;
111 }
112}
113
114void BesEmcGeometry::PrintEMCParameters()
115{
116 // Print EMC parameters
117 ;
118}
119
120void BesEmcGeometry::ComputeEMCParameters()
121{
122 ReadEMCParameters();
123
124 // Compute derived parameters of the calorimeter
125 G4double theta0;
126 G4int i;
127 //for debug
128 //G4Exception("BesEmcGeometry::ComputeEMCParameters() starting.......");
129 //
130 //support frame
131 //
132 const G4double delta=1.*mm; //for opening-cut girder
133 TaperRingRmin1 = BSCRmax1-TaperRingThickness1;
134 TaperRingRmin2 = TaperRingRmin1+TaperRingDz*tan(TaperRingTheta);
135 TaperRingDr = TaperRingThickness2/cos(TaperRingTheta);
136 TaperRingOuterLength1 = TaperRingOuterLength+TaperRingThickness3*tan(TaperRingTheta);
137
138 BSCRmax2 = TaperRingRmin2+TaperRingDr-TaperRingThickness3*tan(TaperRingTheta);
139 BSCRmax = BSCRmin+33.8*cm;
140 BSCPhiDphi= 360.*deg/BSCNbPhi;
141 BSCPhiSphi= -BSCPhiDphi/2.;
142 BSCPhiDz = BSCDz;
143 BSCPhiRmax= BSCRmax-0.5*cm;
144 BSCPhiRmin= BSCRmin/sin(90.*deg+BSCPhiDphi/2.)
145 *sin(90.*deg+BSCPhiDphi/2.-BSCAngleRotat); //Rmin after rotate
146 SPBarDphi=2*atan(SPBarwidth/(2*BSCRmax));
147
148 //crystal No.1(from z=0 to big)
149 zHalfLength[0] = BSCCryLength/2.;
150 yHalfLength1[0]= BSCYFront0/2.;
151 xHalfLength1[0]= BSCPhiRmin*tan(BSCPhiDphi)/2.;
152 xHalfLength2[0]= xHalfLength1[0];
153
154 //there are 3 choices of rminprojected
155 G4double rminprojected=BSCRmin*cos(BSCAngleRotat);
156 rminprojected=BSCRmin;//according to hardware design
157 //rminprojected=BSCPhiRmin;
158 yHalfLength2[0]= yHalfLength1[0]
159 +(BSCYFront0-BSCPosition0)/rminprojected*BSCCryLength/2.;
160 xHalfLength3[0]= xHalfLength1[0]*(BSCPhiRmin+BSCCryLength)/BSCPhiRmin;
161 xHalfLength4[0]= xHalfLength2[0]*(BSCPhiRmin+BSCCryLength)/BSCPhiRmin;
162
163 thetaPosition[0]= 90.*deg;
164 xPosition[0] = BSCPhiRmin+BSCCryLength/2.;
165 yPosition[0] =
166 -(xHalfLength1[0]+xHalfLength3[0]+xHalfLength2[0]+xHalfLength4[0])/4.;
167 zPosition[0] = (yHalfLength1[0]+yHalfLength2[0])/2.;
168
169 theta0= 90.*deg-atan((BSCYFront0-BSCPosition0)/rminprojected);
170 if(verboseLevel>1)
171 {
172 G4cout << "------------------------------------>1"<< G4endl
173 << "The direction of crystal is:"
174 << theta0/deg <<"(deg)."<< G4endl;
175 }
176
177 rearBoxPosX[0] = xPosition[0]+BSCCryLength/2.+rearBoxDz/2.;
178 rearBoxPosY[0] = -(xHalfLength3[0]+xHalfLength4[0])/2.;
179 rearBoxPosZ[0] = yHalfLength2[0];
180
181 OCGirderRmin1[0] = rearBoxPosX[0]+rearBoxDz/2.+delta;
182 OCGirderRmin2[0] = rearBoxPosX[0]+rearBoxDz/2.+delta;
183 OCGirderDz[0] = rearBoxLength;
184 OCGirderPosZ[0] = rearBoxLength/2;
185
186 cableLength[0] = BSCDz-rearBoxPosZ[0];
187 cablePosX[0] = BSCPhiRmax-cableDr-2.*mm;
188 cablePosY[0] = rearBoxPosY[0]-3*cableDr;
189 cablePosZ[0] = BSCDz-cableLength[0]/2;
190
191 //crystal No.2
192 zHalfLength[1]= BSCCryLength/2.;
193 yHalfLength1[1]= BSCYFront0/2.;
194 G4double deltaR= BSCYFront0*cos(theta0);
195 xHalfLength1[1]= xHalfLength1[0];
196 xHalfLength2[1]= xHalfLength1[1]*(BSCPhiRmin+deltaR)/BSCPhiRmin;
197 yHalfLength2[1]= yHalfLength1[1]
198 +tan(theta0-atan(rminprojected/(BSCYFront0*(1.+1./sin(theta0))-
199 BSCPosition1)))*BSCCryLength/2.;
200 deltaR=deltaR+BSCCryLength*sin(theta0);
201 xHalfLength4[1]= xHalfLength1[1]*(BSCPhiRmin+deltaR)/BSCPhiRmin;
202 deltaR=deltaR-yHalfLength2[1]*2.*cos(theta0);
203 xHalfLength3[1]= xHalfLength1[1]*(BSCPhiRmin+deltaR)/BSCPhiRmin;
204
205 thetaPosition[1]= theta0;
206 xPosition[1] = BSCPhiRmin+BSCCryLength/2.*sin(theta0)
207 + (3.*yHalfLength1[1]-yHalfLength2[1])/2.*cos(theta0);
208 yPosition[1] =
209 -(xHalfLength1[1]+xHalfLength3[1]+xHalfLength2[1]+xHalfLength4[1])/4.;
210 zPosition[1] = yHalfLength1[0]*2.
211 + (yHalfLength1[1]*2./tan(theta0)+BSCCryLength/2.)*cos(theta0)
212 + (yHalfLength1[1]+yHalfLength2[1])/2.*sin(theta0);
213
214 rearBoxPosX[1] = xPosition[1]+(BSCCryLength/2.+rearBoxDz/2.)*sin(theta0);
215 rearBoxPosY[1] = -(xHalfLength3[1]+xHalfLength4[1])/2.;
216 rearBoxPosZ[1] = zPosition[1]+(BSCCryLength/2.+rearBoxDz/2.)*cos(theta0);
217
218 OCGirderRmin1[1] = rearBoxPosX[1]+rearBoxDz*sin(theta0)/2.+rearBoxLength*cos(theta0)/2+delta;
219 OCGirderRmin2[1] = rearBoxPosX[1]+rearBoxDz*sin(theta0)/2.-rearBoxLength*cos(theta0)/2+delta;
220 OCGirderDz[1] = rearBoxLength*sin(theta0);
221 OCGirderPosZ[1] = rearBoxPosZ[1]+rearBoxDz*cos(theta0)/2.;
222
223 cableLength[1] = BSCDz-rearBoxPosZ[1];
224 cablePosX[1] = cablePosX[0];
225 cablePosY[1] = cablePosY[0]+2*cableDr;
226 cablePosZ[1] = BSCDz-cableLength[1]/2;
227
228 theta0= theta0-atan((yHalfLength2[1]-yHalfLength1[1])*2./BSCCryLength);
229
230 for(i=2;i<BSCNbTheta;i++)
231 {
232 if(verboseLevel>1)
233 {
234 G4cout << "------------------------------------>"<<i<< G4endl
235 << "The direction of crystal is:"
236 << theta0/deg <<"(deg)." << G4endl;
237 }
238 G4double CryLength;
239 if(i==BSCNbTheta-1) //the rightest crystal
240 {
241 CryLength=BSCCryLength1; //24cm
242 yHalfLength1[i]=BSCYFront1/2.;
243 } else {
244 //crystal No.i+1
245 CryLength=BSCCryLength; //28cm
246 yHalfLength1[i]=BSCYFront/2.;
247 }
248 zHalfLength[i]=CryLength/2;
249
250 deltaR=yHalfLength1[i]*2.*cos(theta0);
251 xHalfLength1[i]=xHalfLength1[0];
252 xHalfLength2[i]=xHalfLength1[i]/BSCPhiRmin*(BSCPhiRmin+deltaR);
253 yHalfLength2[i]=yHalfLength1[i]
254 *(1.+CryLength/(rminprojected/sin(theta0)
255 +yHalfLength1[i]*2./tan(theta0)));
256 deltaR=deltaR+CryLength*sin(theta0);
257 xHalfLength4[i]=xHalfLength1[i]/BSCPhiRmin*(BSCPhiRmin+deltaR);
258 deltaR=deltaR-yHalfLength2[i]*2.*cos(theta0);
259 xHalfLength3[i]=xHalfLength1[i]/BSCPhiRmin*(BSCPhiRmin+deltaR);
260
261 thetaPosition[i]=theta0;
262 xPosition[i]=BSCPhiRmin+zHalfLength[i]*sin(theta0)
263 + (3.*yHalfLength1[i]-yHalfLength2[i])/2.*cos(theta0);
264 yPosition[i]=
265 -(xHalfLength1[i]+xHalfLength3[i]+xHalfLength2[i]+xHalfLength4[i])/4.;
266 zPosition[i]=BSCPosition1+rminprojected/tan(theta0)
267 +(2.*yHalfLength1[i]/tan(theta0)+zHalfLength[i])*cos(theta0)
268 +(yHalfLength1[i]+yHalfLength2[i])/2.*sin(theta0);
269
270 rearBoxPosX[i] = xPosition[i]+(CryLength/2.+rearBoxDz/2.)*sin(theta0);
271 rearBoxPosY[i] = -(xHalfLength3[i]+xHalfLength4[i])/2.;
272 rearBoxPosZ[i] = zPosition[i]+(CryLength/2.+rearBoxDz/2.)*cos(theta0);
273
274 OCGirderRmin1[i] = rearBoxPosX[i]+rearBoxDz*sin(theta0)/2.+rearBoxLength*cos(theta0)/2+delta;
275 OCGirderRmin2[i] = rearBoxPosX[i]+rearBoxDz*sin(theta0)/2.-rearBoxLength*cos(theta0)/2+delta;
276 OCGirderDz[i] = rearBoxLength*sin(theta0);
277 OCGirderPosZ[i] = rearBoxPosZ[i]+rearBoxDz*cos(theta0)/2.;
278
279 G4int xCable,yCable;
280 xCable = i/4;
281 yCable = i-4*(G4int)(i/4);
282 cableLength[i] = BSCDz-(rearBoxPosZ[i]+rearBoxDz/2.*cos(theta0));
283 cablePosX[i] = cablePosX[0]-xCable*2*cableDr;
284 cablePosY[i] = cablePosY[0]+yCable*2*cableDr;
285 cablePosZ[i] = BSCDz-cableLength[i]/2;
286
287 theta0=theta0-atan(2.*yHalfLength1[i]/(rminprojected/sin(theta0)
288 +2.*yHalfLength1[i]/tan(theta0)));
289
290 }
291 thetaPosition[BSCNbTheta]=theta0;
292 if(verboseLevel>1)
293 {
294 G4cout << "------------------------------------>"<<i<< G4endl
295 << "The direction of crystal is:"
296 << theta0/deg <<"(deg)." << G4endl;
297 }
298
299 G4double oop;
300 G4double dx=0.001*mm; //zoom out of crystals to avoid overlap
301 for(i=0;i<BSCNbTheta;i++)
302 {
303 xHalfLength1[i]-=dx;
304 xHalfLength2[i]-=dx;
305 xHalfLength3[i]-=dx;
306 xHalfLength4[i]-=dx;
307 yHalfLength1[i]-=dx;
308 yHalfLength2[i]-=dx;
309
310 G4double CryLength;
311 if(i==BSCNbTheta-1) //the rightest crystal
312 {
313 CryLength=BSCCryLength1;
314 } else {
315 CryLength=BSCCryLength;
316 }
317 oop=sqrt((yHalfLength2[i]-yHalfLength1[i])*(yHalfLength2[i]
318 -yHalfLength1[i])
319 +(xHalfLength3[i]+xHalfLength4[i]-xHalfLength1[i]
320 -xHalfLength2[i])*(xHalfLength3[i]+xHalfLength4[i]
321 -xHalfLength1[i]-xHalfLength2[i])/4);
322 thetaAxis[i]=atan(oop/CryLength);
323 phiAxis[i] =180.*deg+atan((yHalfLength2[i]-yHalfLength1[i])
324 /(xHalfLength3[i]+xHalfLength4[i]
325 -xHalfLength1[i]-xHalfLength2[i])*2.);
326 tanAlpha2[i]=-(xHalfLength4[i]-xHalfLength3[i])/yHalfLength2[i]/2.;
327 tanAlpha1[i]=-(xHalfLength2[i]-xHalfLength1[i])/yHalfLength1[i]/2.;
328 //G4cout <<thetaAxis[i]<<", "
329 // <<phiAxis[i]<<", "
330 // <<tanAlpha1[i]<<", "
331 // <<tanAlpha2[i]<<G4endl;
332 }
333
334}
335
336void BesEmcGeometry::ModifyForCasing()
337{
338 // Compute the sizes of the naked crystals and the casing
339 // Casing size
340 // BSCNbTheta---->2*BSCNbTheta-1 : before crystal
341 // 2*BSCNbTheta-->3*BSCNbTheta-1 : a side of crystal
342 // 3*BSCNbTheta-->4*BSCNbTheta-1 : b side of crystal
343 // 4*BSCNbTheta-->5*BSCNbTheta-1 : c side of crystal
344 // 5*BSCNbTheta-->6*BSCNbTheta-1 : d side of crystal
345 // d
346 // ----------
347 // | |
348 // | |
349 // c | | a
350 // | |
351 // | |
352 // | /
353 // | /
354 // | /
355 // b
356 G4double totalThickness=fTyvekThickness+fAlThickness+fMylarThickness;//
357 G4double delta=0.,angle1=0.*deg,angle2=0.*deg;
358 G4double oop; //the distance of the centers of the two parallel plane
359
360 G4double rminprojected=BSCRmin*cos(BSCAngleRotat);
361 rminprojected=BSCRmin;//accord with hardware design
362 //rminprojected=BSCPhiRmin;
363
364 G4int i;
365 for(i=0;i<BSCNbTheta;i++)
366 {
367 G4double CryLength;
368 if(i==BSCNbTheta-1) //the rightest crystal
369 {
370 CryLength=BSCCryLength1;
371 } else {
372 CryLength=BSCCryLength;
373 }
374 zHalfLength[BSCNbTheta+i]=totalThickness/2.;
375 yHalfLength1[BSCNbTheta+i]=yHalfLength1[i];
376 yHalfLength2[BSCNbTheta+i]=yHalfLength1[i]
377 +(yHalfLength2[i]-yHalfLength1[i])*totalThickness/CryLength;
378 xHalfLength1[BSCNbTheta+i]=xHalfLength1[i];
379 xHalfLength2[BSCNbTheta+i]=xHalfLength2[i];
380 xHalfLength1[BSCNbTheta*2+i]=xHalfLength3[BSCNbTheta+i]=
381 xHalfLength1[i]*(CryLength-totalThickness)/CryLength
382 +xHalfLength3[i]*totalThickness/CryLength;
383 xHalfLength2[BSCNbTheta*4+i]=xHalfLength4[BSCNbTheta+i]=
384 xHalfLength2[i]*(CryLength-totalThickness)/CryLength
385 +xHalfLength4[i]*totalThickness/CryLength;
386
387 zHalfLength[BSCNbTheta*5+i]=zHalfLength[BSCNbTheta*4+i]=
388 zHalfLength[BSCNbTheta*3+i]=zHalfLength[BSCNbTheta*2+i]=
389 zHalfLength[i]-totalThickness/2.;
390
391 yHalfLength2[BSCNbTheta*2+i]=yHalfLength1[BSCNbTheta*2+i]=
392 totalThickness/cos(thetaPosition[i+1]-thetaPosition[i])/2.;
393 xHalfLength3[BSCNbTheta*2+i]=xHalfLength3[i];
394 xHalfLength4[BSCNbTheta*2+i]=xHalfLength3[i]
395 +(xHalfLength4[i]-xHalfLength3[i])*yHalfLength2[BSCNbTheta*2+i]
396 /yHalfLength2[i];
397 xHalfLength2[BSCNbTheta*2+i]=xHalfLength3[BSCNbTheta+i]
398 +(xHalfLength4[BSCNbTheta+i]-xHalfLength3[BSCNbTheta+i])
399 *yHalfLength1[BSCNbTheta*2+i]/yHalfLength2[BSCNbTheta*1+i];
400
401 yHalfLength2[BSCNbTheta*4+i]=yHalfLength1[BSCNbTheta*4+i]=
402 totalThickness/2.;
403 xHalfLength4[BSCNbTheta*4+i]=xHalfLength4[i];
404 xHalfLength3[BSCNbTheta*4+i]=xHalfLength4[i]
405 -(xHalfLength4[i]-xHalfLength3[i])*yHalfLength2[BSCNbTheta*4+i]
406 /yHalfLength2[i];
407 xHalfLength1[BSCNbTheta*4+i]=xHalfLength4[BSCNbTheta+i]
408 -(xHalfLength4[BSCNbTheta+i]-xHalfLength3[BSCNbTheta+i])
409 *yHalfLength1[BSCNbTheta*4+i]/yHalfLength2[BSCNbTheta*1+i];
410
411 delta=totalThickness/2.+yHalfLength1[BSCNbTheta*2+i];
412 angle1=atan(yHalfLength2[i]/(xHalfLength4[i]-xHalfLength3[i]));
413 angle2=atan(2.*(xHalfLength4[i]-xHalfLength2[i])*sin(angle1)
414 /CryLength);
415
416 yHalfLength1[BSCNbTheta*5+i]=yHalfLength1[BSCNbTheta*3+i]=
417 yHalfLength1[i]-delta;
418 yHalfLength2[BSCNbTheta*5+i]=yHalfLength2[BSCNbTheta*3+i]=
419 yHalfLength2[i]-delta;
420 xHalfLength4[BSCNbTheta*3+i]=xHalfLength3[BSCNbTheta*3+i]=
421 xHalfLength2[BSCNbTheta*3+i]=xHalfLength1[BSCNbTheta*3+i]=
422 totalThickness/cos(angle2)/sin(angle1)/2.;
423 xHalfLength4[BSCNbTheta*5+i]=xHalfLength3[BSCNbTheta*5+i]=
424 xHalfLength2[BSCNbTheta*5+i]=xHalfLength1[BSCNbTheta*5+i]=
425 totalThickness/2.;
426
427 zHalfLength[i]=zHalfLength[i]-totalThickness/2.;
428 yHalfLength1[i]=yHalfLength1[i]-delta;
429 yHalfLength2[i]=yHalfLength2[i]-delta;
430 delta=totalThickness*(1.+1./cos(angle2)/sin(angle1))/2.;
431 xHalfLength1[i]=xHalfLength1[i]-delta;
432 xHalfLength2[i]=xHalfLength2[i]-delta;
433 xHalfLength3[i]=xHalfLength3[i]-delta;
434 xHalfLength4[i]=xHalfLength4[i]-delta;
435
436 oop=sqrt((yHalfLength2[i]-yHalfLength1[i])*(yHalfLength2[i]
437 -yHalfLength1[i])
438 +(xHalfLength3[i]+xHalfLength4[i]-xHalfLength1[i]
439 -xHalfLength2[i])*(xHalfLength3[i]+xHalfLength4[i]
440 -xHalfLength1[i]-xHalfLength2[i])/4);
441 thetaAxis[i]=atan(oop/CryLength);
442 // -phi --->180+phi
443 phiAxis[i] =180.*deg+atan((yHalfLength2[i]-yHalfLength1[i])
444 /(xHalfLength3[i]+xHalfLength4[i]
445 -xHalfLength1[i]-xHalfLength2[i])*2.);
446
447 oop=sqrt((yHalfLength2[BSCNbTheta+i]-yHalfLength1[BSCNbTheta+i])
448 *(yHalfLength2[BSCNbTheta+i]-yHalfLength1[BSCNbTheta+i])
449 +(xHalfLength3[BSCNbTheta+i]+xHalfLength4[BSCNbTheta+i]
450 -xHalfLength1[BSCNbTheta+i]-xHalfLength2[BSCNbTheta+i])
451 *(xHalfLength3[BSCNbTheta+i]+xHalfLength4[BSCNbTheta+i]
452 -xHalfLength1[BSCNbTheta+i]-xHalfLength2[BSCNbTheta+i])/4);
453 thetaAxis[BSCNbTheta+i]=atan(oop/totalThickness);
454 phiAxis [BSCNbTheta+i]=
455 -atan((yHalfLength2[BSCNbTheta+i]-yHalfLength1[BSCNbTheta+i])
456 /(xHalfLength3[BSCNbTheta+i]+xHalfLength4[BSCNbTheta+i]
457 -xHalfLength1[BSCNbTheta+i]-xHalfLength2[BSCNbTheta+i])*2.);
458
459 oop=sqrt((yHalfLength2[i]-yHalfLength1[i])*(yHalfLength2[i]
460 -yHalfLength1[i])*4
461 +(xHalfLength3[BSCNbTheta*2+i]+xHalfLength4[BSCNbTheta*2+i]
462 -xHalfLength1[BSCNbTheta*2+i]-xHalfLength2[BSCNbTheta*2+i])
463 *(xHalfLength3[BSCNbTheta*2+i]+xHalfLength4[BSCNbTheta*2+i]
464 -xHalfLength1[BSCNbTheta*2+i]-xHalfLength2[BSCNbTheta*2+i])
465 /4);
466 thetaAxis[BSCNbTheta*2+i]=atan(oop/(zHalfLength[BSCNbTheta*2+i]*2));
467 phiAxis [BSCNbTheta*2+i]=
468 -atan((yHalfLength2[i]-yHalfLength1[i])
469 /(xHalfLength3[BSCNbTheta*2+i]+xHalfLength4[BSCNbTheta*2+i]
470 -xHalfLength1[BSCNbTheta*2+i]-xHalfLength2[BSCNbTheta*2+i])*4);
471
472 oop=sqrt((yHalfLength2[i]-yHalfLength1[i])*(yHalfLength2[i]
473 -yHalfLength1[i])*4
474 +(xHalfLength4[i]-xHalfLength2[i])
475 *(xHalfLength4[i]-xHalfLength2[i])*4);
476 thetaAxis[BSCNbTheta*3+i]=atan(oop/(zHalfLength[BSCNbTheta*3+i]*2));
477 phiAxis [BSCNbTheta*3+i]=-atan((yHalfLength2[i]-yHalfLength1[i])
478 /(xHalfLength4[i]-xHalfLength2[i]));
479
480 thetaAxis[BSCNbTheta*4+i]=
481 atan((xHalfLength4[BSCNbTheta*4+i]+xHalfLength3[BSCNbTheta*4+i]
482 -xHalfLength2[BSCNbTheta*4+i]-xHalfLength1[BSCNbTheta*4+i])/2.
483 /(zHalfLength[BSCNbTheta*4+i]*2));
484 phiAxis [BSCNbTheta*4+i]=0;
485
486 thetaAxis[BSCNbTheta*5+i]=atan((xHalfLength3[BSCNbTheta*5+i]
487 -xHalfLength1[BSCNbTheta*5+i])
488 /(zHalfLength[BSCNbTheta*5+i]*2));
489 phiAxis [BSCNbTheta*5+i]=-90.*deg;
490
491 tanAlpha2[BSCNbTheta+i]=tanAlpha1[BSCNbTheta+i]=tanAlpha1[i]=
492 -(xHalfLength2[i]-xHalfLength1[i])/yHalfLength1[i]/2.;
493 tanAlpha1[BSCNbTheta*2+i]=(xHalfLength2[BSCNbTheta*2+i]
494 -xHalfLength1[BSCNbTheta*2+i])
495 /yHalfLength1[BSCNbTheta*2+i]/2.;
496 tanAlpha1[BSCNbTheta*3+i]=tanAlpha1[i]*2.;
497 tanAlpha1[BSCNbTheta*4+i]=(xHalfLength2[BSCNbTheta*4+i]
498 -xHalfLength1[BSCNbTheta*4+i])
499 /yHalfLength1[BSCNbTheta*4+i]/2.;
500 tanAlpha1[BSCNbTheta*5+i]=(xHalfLength2[BSCNbTheta*5+i]
501 -xHalfLength1[BSCNbTheta*5+i])
502 /yHalfLength1[BSCNbTheta*5+i]/2.;
503
504 tanAlpha2[i]=-(xHalfLength4[i]-xHalfLength3[i])/yHalfLength2[i]/2.;
505 //tanAlpha2[BSCNbTheta+i]=(xHalfLength4[BSCNbTheta+i]-xHalfLength3[BSCNbTheta+i])/yHalfLength2[BSCNbTheta+i]/2.;
506 tanAlpha2[BSCNbTheta*2+i]=(xHalfLength4[BSCNbTheta*2+i]
507 -xHalfLength3[BSCNbTheta*2+i])
508 /yHalfLength2[BSCNbTheta*2+i]/2.;
509 tanAlpha2[BSCNbTheta*3+i]=tanAlpha2[i]*2.;
510 tanAlpha2[BSCNbTheta*4+i]=(xHalfLength4[BSCNbTheta*4+i]
511 -xHalfLength3[BSCNbTheta*4+i])
512 /yHalfLength2[BSCNbTheta*4+i]/2.;
513 tanAlpha2[BSCNbTheta*5+i]=(xHalfLength4[BSCNbTheta*5+i]
514 -xHalfLength3[BSCNbTheta*5+i])
515 /yHalfLength2[BSCNbTheta*5+i]/2.;
516
517 zPosition[BSCNbTheta*5+i]=zPosition[BSCNbTheta*3+i]=zPosition[i]=
518 zPosition[i]+totalThickness/2.*cos(thetaPosition[i])
519 -yHalfLength1[BSCNbTheta*2+i]*sin(thetaPosition[i]);
520 zPosition[i]=totalThickness/2.;//for the newest method
521 xPosition[BSCNbTheta*5+i]=xPosition[BSCNbTheta*3+i]=xPosition[i]=
522 xPosition[i]+totalThickness/2.*sin(thetaPosition[i])
523 +totalThickness*(1./cos(thetaPosition[i+1]-thetaPosition[i])-1)/2.
524 *cos(thetaPosition[i]);
525 xPosition[i]=totalThickness*(1.-1./cos(angle2)/sin(angle1))/2.;
526 //for the newest method
527 yPosition[i]=yPosition[i]
528 +totalThickness*(1.-1./cos(angle2)/sin(angle1))/2.;
529 yPosition[i]=yHalfLength1[BSCNbTheta*2+i]-totalThickness/2.;//for the newest method
530 yPosition[BSCNbTheta*3+i]=yPosition[i]*2.+xHalfLength1[BSCNbTheta*3+i];
531 yPosition[BSCNbTheta*5+i]=xHalfLength1[BSCNbTheta*5+i];
532
533 xPosition[BSCNbTheta+i]=BSCPhiRmin
534 +zHalfLength[BSCNbTheta+i]*sin(thetaPosition[i])
535 +(3.*yHalfLength1[BSCNbTheta+i]-yHalfLength2[BSCNbTheta+i])/2.
536 *cos(thetaPosition[i]);
537 yPosition[BSCNbTheta+i]=(xHalfLength1[BSCNbTheta+i]
538 +xHalfLength3[BSCNbTheta+i]
539 +xHalfLength2[BSCNbTheta+i]
540 +xHalfLength4[BSCNbTheta+i])/4.;
541 zPosition[BSCNbTheta+i]=BSCPosition1+rminprojected/tan(thetaPosition[i])
542 +(2.*yHalfLength1[BSCNbTheta+i]/tan(thetaPosition[i])
543 +zHalfLength[BSCNbTheta+i])*cos(thetaPosition[i])
544 +(yHalfLength1[BSCNbTheta+i]+yHalfLength2[BSCNbTheta+i])/2.
545 *sin(thetaPosition[i]);
546
547 xPosition[BSCNbTheta*2+i]=xPosition[i]
548 +((yHalfLength1[i]+yHalfLength2[i])/2.+yHalfLength1[BSCNbTheta*2+i])
549 *cos(thetaPosition[i]);
550 zPosition[BSCNbTheta*2+i]=zPosition[i]
551 -((yHalfLength1[i]+yHalfLength2[i])/2.+yHalfLength1[BSCNbTheta*2+i])
552 *sin(thetaPosition[i]);
553 yPosition[BSCNbTheta*2+i]=(xHalfLength1[BSCNbTheta*2+i]
554 +xHalfLength3[BSCNbTheta*2+i]
555 +xHalfLength2[BSCNbTheta*2+i]
556 +xHalfLength4[BSCNbTheta*2+i])/4.;
557
558 xPosition[BSCNbTheta*4+i]=xPosition[i]
559 -((yHalfLength1[i]+yHalfLength2[i])/2.+yHalfLength1[BSCNbTheta*4+i])
560 *cos(thetaPosition[i]);
561 zPosition[BSCNbTheta*4+i]=zPosition[i]
562 -((yHalfLength1[i]+yHalfLength2[i])/2.+yHalfLength1[BSCNbTheta*4+i])
563 *sin(thetaPosition[i]);
564 yPosition[BSCNbTheta*4+i]=(xHalfLength1[BSCNbTheta*4+i]
565 +xHalfLength3[BSCNbTheta*4+i]
566 +xHalfLength2[BSCNbTheta*4+i]
567 +xHalfLength4[BSCNbTheta*4+i])/4.;
568
569 }
570
571 if(verboseLevel>1)
572 for(i=0;i<BSCNbTheta*6;i++)
573 {
574 G4cout << "The sizes of the " << i+1 << " crystal are:" << G4endl
575 << "zHalfLength =" << zHalfLength[i]/cm << "(cm)," << G4endl
576 << "thetaAxis =" << thetaAxis[i]/deg << "(deg),"<< G4endl
577 << "phiAxis =" << phiAxis[i]/deg << "(deg),"<< G4endl
578 << "yHalfLength1=" << yHalfLength1[i]/cm << "(cm)," << G4endl
579 << "xHalfLength1=" << xHalfLength1[i]/cm << "(cm)," << G4endl
580 << "xHalfLength2=" << xHalfLength2[i]/cm << "(cm)," << G4endl
581 << "tanAlpha1 =" << tanAlpha1[i] << G4endl
582 << "yHalfLength2=" << yHalfLength2[i]/cm << "(cm)," << G4endl
583 << "xHalfLength3=" << xHalfLength3[i]/cm << "(cm)," << G4endl
584 << "xHalfLength4=" << xHalfLength4[i]/cm << "(cm)," << G4endl
585 << "tanAlpha2 =" << tanAlpha2[i] << "." << G4endl
586 << "The position of the " << i+1 << " crystal is:" << G4endl
587 << "(" << xPosition[i]/cm << ","
588 << yPosition[i]/cm << ","
589 << zPosition[i]/cm << ")cm" << G4endl;
590 }
591
592}
593
594void BesEmcGeometry::SetCasingThickness(G4ThreeVector val)
595{
596 // change Gap thickness and recompute the calorimeter parameters
597 fTyvekThickness = val('X');
598 fAlThickness = val('Y');
599 fMylarThickness = val('Z');
600}
601
602G4double BesEmcGeometry::GetXPosition(G4int NbCrystal)
603{
604 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
605 {
606 return xPosition[NbCrystal];
607 }
608 else
609 {
610 return 0.;
611 }
612}
613
614G4double BesEmcGeometry::GetYPosition(G4int NbCrystal)
615{
616 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
617 {
618 return yPosition[NbCrystal];
619 }
620 else
621 {
622 return 0.;
623 }
624}
625
626G4double BesEmcGeometry::GetZPosition(G4int NbCrystal)
627{
628 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
629 {
630 return zPosition[NbCrystal];
631 }
632 else
633 {
634 return 0.;
635 }
636}
637
638G4double BesEmcGeometry::GetThetaPosition(G4int NbCrystal)
639{
640 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
641 {
642 return thetaPosition[NbCrystal];
643 }
644 else
645 {
646 return 0.;
647 }
648}
649
650G4double BesEmcGeometry::GetZHalfLength(G4int NbCrystal)
651{
652 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
653 {
654 return zHalfLength[NbCrystal];
655 }
656 else
657 {
658 return 0.;
659 }
660}
661
662G4double BesEmcGeometry::GetThetaAxis(G4int NbCrystal)
663{
664 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
665 {
666 return thetaAxis[NbCrystal];
667 }
668 else
669 {
670 return 0.;
671 }
672}
673
674G4double BesEmcGeometry::GetPhiAxis(G4int NbCrystal)
675{
676 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
677 {
678 return phiAxis[NbCrystal];
679 }
680 else
681 {
682 return 0.;
683 }
684}
685
686G4double BesEmcGeometry::GetYHalfLength1(G4int NbCrystal)
687{
688 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
689 {
690 return yHalfLength1[NbCrystal];
691 }
692 else
693 {
694 return 0.;
695 }
696}
697
698G4double BesEmcGeometry::GetXHalfLength1(G4int NbCrystal)
699{
700 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
701 {
702 return xHalfLength1[NbCrystal];
703 }
704 else
705 {
706 return 0.;
707 }
708}
709
710G4double BesEmcGeometry::GetXHalfLength2(G4int NbCrystal)
711{
712 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
713 {
714 return xHalfLength2[NbCrystal];
715 }
716 else
717 {
718 return 0.;
719 }
720}
721
722G4double BesEmcGeometry::GetTanAlpha1(G4int NbCrystal)
723{
724 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
725 {
726 return tanAlpha1[NbCrystal];
727 }
728 else
729 {
730 return 0.;
731 }
732}
733
734G4double BesEmcGeometry::GetYHalfLength2(G4int NbCrystal)
735{
736 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
737 {
738 return yHalfLength2[NbCrystal];
739 }
740 else
741 {
742 return 0.;
743 }
744}
745
746G4double BesEmcGeometry::GetXHalfLength3(G4int NbCrystal)
747{
748 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
749 {
750 return xHalfLength3[NbCrystal];
751 }
752 else
753 {
754 return 0.;
755 }
756}
757
758G4double BesEmcGeometry::GetXHalfLength4(G4int NbCrystal)
759{
760 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
761 {
762 return xHalfLength4[NbCrystal];
763 }
764 else
765 {
766 return 0.;
767 }
768}
769
770G4double BesEmcGeometry::GetTanAlpha2(G4int NbCrystal)
771{
772 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
773 {
774 return tanAlpha2[NbCrystal];
775 }
776 else
777 {
778 return 0.;
779 }
780}
781
782G4VPhysicalVolume* BesEmcGeometry::GetPhysiBSCCrystal(G4int NbCrystal)
783{
784 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
785 {
786 return physiBSCCrystal[NbCrystal];
787 }
788 else
789 {
790 return NULL;
791 }
792}
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
BesEmcGeometry::ComputeEMCParameters
void ComputeEMCParameters()
Definition: BesEmcGeometry.cc:120
BesEmcGeometry::GetThetaPosition
G4double GetThetaPosition(G4int NbCrystal)
Definition: BesEmcGeometry.cc:638
BesEmcGeometry::ReadEMCParameters
void ReadEMCParameters()
Definition: BesEmcGeometry.cc:24
BesEmcGeometry::GetYPosition
G4double GetYPosition(G4int NbCrystal)
Definition: BesEmcGeometry.cc:614
BesEmcGeometry::GetZHalfLength
G4double GetZHalfLength(G4int NbCrystal)
Definition: BesEmcGeometry.cc:650
BesEmcGeometry::GetXHalfLength2
G4double GetXHalfLength2(G4int NbCrystal)
Definition: BesEmcGeometry.cc:710
BesEmcGeometry::GetTanAlpha2
G4double GetTanAlpha2(G4int NbCrystal)
Definition: BesEmcGeometry.cc:770
BesEmcGeometry::GetTanAlpha1
G4double GetTanAlpha1(G4int NbCrystal)
Definition: BesEmcGeometry.cc:722
BesEmcGeometry::GetYHalfLength1
G4double GetYHalfLength1(G4int NbCrystal)
Definition: BesEmcGeometry.cc:686
BesEmcGeometry::GetThetaAxis
G4double GetThetaAxis(G4int NbCrystal)
Definition: BesEmcGeometry.cc:662
BesEmcGeometry::GetXHalfLength4
G4double GetXHalfLength4(G4int NbCrystal)
Definition: BesEmcGeometry.cc:758
BesEmcGeometry::GetYHalfLength2
G4double GetYHalfLength2(G4int NbCrystal)
Definition: BesEmcGeometry.cc:734
BesEmcGeometry::GetPhysiBSCCrystal
G4VPhysicalVolume * GetPhysiBSCCrystal(G4int NbCrystal)
Definition: BesEmcGeometry.cc:782
BesEmcGeometry::GetXPosition
G4double GetXPosition(G4int NbCrystal)
Definition: BesEmcGeometry.cc:602
BesEmcGeometry::GetZPosition
G4double GetZPosition(G4int NbCrystal)
Definition: BesEmcGeometry.cc:626
BesEmcGeometry::SetCasingThickness
void SetCasingThickness(G4ThreeVector)
Definition: BesEmcGeometry.cc:594
BesEmcGeometry::GetXHalfLength1
G4double GetXHalfLength1(G4int NbCrystal)
Definition: BesEmcGeometry.cc:698
BesEmcGeometry::GetPhiAxis
G4double GetPhiAxis(G4int NbCrystal)
Definition: BesEmcGeometry.cc:674
BesEmcGeometry::PrintEMCParameters
void PrintEMCParameters()
Definition: BesEmcGeometry.cc:114
BesEmcGeometry::GetXHalfLength3
G4double GetXHalfLength3(G4int NbCrystal)
Definition: BesEmcGeometry.cc:746
BesEmcParameter::GetInstance
static BesEmcParameter & GetInstance()
Definition: BesEmcParameter.cc:33