CGEM BOSS 6.6.5.h
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
17{
18 verboseLevel = 0;
19}
20
23
25{
26 // Read EMC parameters from database
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
115{
116 // Print EMC parameters
117 ;
118}
119
121{
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
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
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}
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
void ComputeEMCParameters()
G4double GetThetaPosition(G4int NbCrystal)
G4double GetYPosition(G4int NbCrystal)
G4double GetZHalfLength(G4int NbCrystal)
G4double GetXHalfLength2(G4int NbCrystal)
G4double GetTanAlpha2(G4int NbCrystal)
G4double GetTanAlpha1(G4int NbCrystal)
G4double GetYHalfLength1(G4int NbCrystal)
G4double GetThetaAxis(G4int NbCrystal)
G4double GetXHalfLength4(G4int NbCrystal)
G4double GetYHalfLength2(G4int NbCrystal)
G4VPhysicalVolume * GetPhysiBSCCrystal(G4int NbCrystal)
G4double GetXPosition(G4int NbCrystal)
G4double GetZPosition(G4int NbCrystal)
void SetCasingThickness(G4ThreeVector)
G4double GetXHalfLength1(G4int NbCrystal)
G4double GetPhiAxis(G4int NbCrystal)
G4double GetXHalfLength3(G4int NbCrystal)
G4double GetOrgGlassLengthY()
G4double GetBSCYFront0()
G4double GetBSCPosition1()
G4double GetEndRingDr()
G4double GetCrystalLength()
G4double GetBSCRmax1()
G4double GetRearCasingThickness()
G4double GetOrgGlassLengthZ()
G4double GetOCGirderAngle()
G4double GetBSCRmin2()
G4double GetPDLengthX()
G4double GetSPBarThickness()
G4double GetBSCYFront1()
G4double GetBSCPosition0()
G4double GetTaperRingInnerLength()
G4double GetTaperRingThickness1()
static BesEmcParameter & GetInstance()
G4double GetBSCYFront()
G4double GetTyvekThickness()
G4double GetOrgGlassLengthX()
G4double GetTaperRingDz()
G4double GetRearBoxDz()
G4double GetSPBarwidth()
G4double GetTaperRingThickness3()
G4double GetBSCRmin1()
G4double GetMylarThickness()
G4double GetBSCRmin()
G4double GetWaterPipeThickness()
G4double GetPDLengthY()
G4double GetPDLengthZ()
G4double GetAlPlateDz()
G4double GetWaterPipeDr()
G4double GetTaperRingTheta()
G4double GetHangingPlateDz()
G4double GetRearBoxLength()
G4double GetEndRingRmin()
G4double GetCrystalLength1()
G4double GetEndRingDz()
G4double GetSPBarThickness1()
G4double GetTaperRingThickness2()
G4double GetTaperRingOuterLength()
G4double GetAlThickness()
G4double GetPABoxThickness()
G4double GetBSCAngleRotat()