BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
ExtBesEmcEndGeometry.cxx
Go to the documentation of this file.
1#include "G4ThreeVector.hh"
2
3#include "ExtBesEmcEndGeometry.h"
4#include "ExtBesEmcParameter.h"
5#include "G4ios.hh"
6
7#include <assert.h>
8#include <fstream>
9#include <iomanip>
10#include <iostream>
11#include <sstream>
12// #include "ReadBoostRoot.hh"
13#include "G4SystemOfUnits.hh"
14
15using namespace std;
16
17ExtBesEmcEndGeometry::ExtBesEmcEndGeometry() { ; }
18
19ExtBesEmcEndGeometry::~ExtBesEmcEndGeometry() { ; }
20
21void ExtBesEmcEndGeometry::ReadParameters() {
22 ExtBesEmcParameter emcPara;
23 emcPara.ReadData();
24
25 WorldRmin1 = emcPara.GetWorldRmin1();
26 WorldRmax1 = emcPara.GetWorldRmax1();
27 WorldRmin2 = emcPara.GetWorldRmin2();
28 WorldRmax2 = emcPara.GetWorldRmax2();
29 WorldDz = emcPara.GetWorldDz();
30 WorldZPosition = emcPara.GetWorldZPosition();
31 CrystalLength = emcPara.GetCrystalLength();
32
33 for ( G4int i = 0; i < 6; i++ ) cryNumInOneLayer[i] = emcPara.GetCryInOneLayer( i );
34 for ( G4int i = 0; i < 5; i++ ) pentaInOneSector[i] = emcPara.GetPentaInOneSector( i );
35
36 fTyvekThickness = emcPara.GetTyvekThickness();
37 fAlThickness = emcPara.GetAlThickness();
38 fMylarThickness = emcPara.GetMylarThickness();
39 totalThickness = fTyvekThickness + fAlThickness + fMylarThickness;
40
41 // G4String ParaPath = ReadBoostRoot::GetBoostRoot();
42 G4String ParaPath = getenv( "SIMUTILDATAROOT" );
43 if ( !ParaPath )
44 {
45 G4cout << "BOOST environment not set!" << G4endl;
46 exit( -1 );
47 }
48 ParaPath += "/dat/EmcEndGeometry.dat";
49
50 ifstream fin;
51 fin.open( ParaPath );
52 assert( fin );
53 for ( G4int i = 0; i < 30; i++ )
54 for ( G4int j = 0; j < 24; j++ ) fin >> param[i][j];
55 for ( G4int i = 0; i < 5; i++ )
56 for ( G4int j = 0; j < 6; j++ ) fin >> penta[i][j];
57 fin.close();
58}
59
60void ExtBesEmcEndGeometry::ComputeParameters() {
61
62 ////////////////////////////////////////////////////////////////////////
63 // emc endcap crystal //
64 //////////////////////////////////////////////////////////////////////////
65 // 1 //
66 // 0 __ -- //
67 // -- \ //
68 // | \ \ //
69 // | \ __ - 2 //
70 // | \ -- | //
71 // | | 3 | //
72 // | | | //
73 // | | | //
74 // | | | //
75 // | | 5 | //
76 // 4 \ | | //
77 // \ | __ - 6 //
78 // \| -- //
79 // 7 //
80 /////////////////////////////////////////////////////////////////////////
81
82 ReadParameters();
83
84 G4int pentaNb = 0;
85 for ( G4int i = 0; i < 30; i++ )
86 {
87 for ( G4int j = 0; j < 8; j++ )
88 {
89 // use 24 parameters to construct 8 point of one crystal
90 G4ThreeVector* pPnt =
91 new G4ThreeVector( param[i][j], param[i][j + 8], param[i][j + 16] - WorldZPosition );
92 fPnt[i][j] = *pPnt;
93 }
94
95 if ( i == 5 || i == 6 || i == 14 || i == 15 || i == 16 )
96 {
97 for ( G4int j = 0; j < 8; j++ ) fPnt[30 + pentaNb][j] = fPnt[i][j];
98
99 // compute the 5th point of the pentagonal crystal
100 G4double y0, y1, y4, y5;
101 G4ThreeVector v0, v1, v4, v5;
102 y0 = penta[pentaNb][0];
103 y1 = penta[pentaNb][1];
104 y4 = penta[pentaNb][2];
105 y5 = penta[pentaNb][3];
106 v0 = ( y0 - fPnt[i][3].getY() ) * ( fPnt[i][0] - fPnt[i][3] ) /
107 ( fPnt[i][0].getY() - fPnt[i][3].getY() ) +
108 fPnt[i][3];
109 v1 = ( y1 - fPnt[i][2].getY() ) * ( fPnt[i][1] - fPnt[i][2] ) /
110 ( fPnt[i][1].getY() - fPnt[i][2].getY() ) +
111 fPnt[i][2];
112 v4 = ( y4 - fPnt[i][7].getY() ) * ( fPnt[i][4] - fPnt[i][7] ) /
113 ( fPnt[i][4].getY() - fPnt[i][7].getY() ) +
114 fPnt[i][7];
115 v5 = ( y5 - fPnt[i][6].getY() ) * ( fPnt[i][5] - fPnt[i][6] ) /
116 ( fPnt[i][5].getY() - fPnt[i][6].getY() ) +
117 fPnt[i][6];
118
119 G4double x1, x5;
120 x1 = penta[pentaNb][4];
121 x5 = penta[pentaNb][5];
122 v1 = ( x1 - v0.getX() ) * ( v1 - v0 ) / ( v1.getX() - v0.getX() ) +
123 v0; // v1', the fifth point
124 v5 = ( x5 - v4.getX() ) * ( v5 - v4 ) / ( v5.getX() - v4.getX() ) + v4; // v5'
125
126 fPnt[i][0] = v0;
127 fPnt[i][1] = v1;
128 fPnt[i][4] = v4;
129 fPnt[i][5] = v5;
130
131 fPnt[30 + pentaNb][2] = v1;
132 fPnt[30 + pentaNb][3] = v0;
133 fPnt[30 + pentaNb][6] = v5;
134 fPnt[30 + pentaNb][7] = v4;
135
136 pentaNb++;
137 }
138 }
139
140 // reflect point in one sector according to new design
141 G4ThreeVector temp[35][8];
142 for ( G4int i = 0; i < 35; i++ )
143 {
144 for ( G4int j = 0; j < 8; j++ )
145 {
146 temp[i][j] = fPnt[i][j];
147 fPnt[i][j].rotateZ( 157.5 * deg );
148 fPnt[i][j].setX( -fPnt[i][j].getX() );
149 }
150
151 // point 0<-->3, 1<-->2, 4<-->7, 5<-->6
152 for ( G4int j = 0; j < 8; j++ )
153 {
154 if ( j < 2 )
155 {
156 G4ThreeVector v = fPnt[i][j];
157 fPnt[i][j] = fPnt[i][3 - j];
158 fPnt[i][3 - j] = v;
159 }
160 else if ( j >= 4 && j < 6 )
161 {
162 G4ThreeVector v = fPnt[i][j];
163 fPnt[i][j] = fPnt[i][11 - j];
164 fPnt[i][11 - j] = v;
165 }
166 }
167 }
168
169 // exchange sequence in the same layer
170 Exchange( 0, 3 );
171 Exchange( 1, 2 );
172 Exchange( 4, 7 );
173 Exchange( 5, 31 );
174 Exchange( 6, 30 );
175 Exchange( 8, 12 );
176 Exchange( 9, 11 );
177 Exchange( 13, 17 );
178 Exchange( 14, 34 );
179 Exchange( 15, 33 );
180 Exchange( 16, 32 );
181 Exchange( 18, 23 );
182 Exchange( 19, 22 );
183 Exchange( 20, 21 );
184 Exchange( 24, 29 );
185 Exchange( 25, 28 );
186 Exchange( 26, 27 );
187
188 /*for(G4int i=0;i<35;i++)
189 {
190 G4cout<<"crystal number: "<<i<<G4endl;
191 for(G4int j=0;j<8;j++)
192 G4cout<<fPnt[i][j]<<G4endl;
193 }*/
194
195 // compute the 6 crystal beside the 20mm gap
196 for ( G4int i = 0; i < 35; i++ )
197 {
198 for ( G4int j = 0; j < 8; j++ )
199 {
200 G4ThreeVector pPnt1 = temp[i][j];
201 fPnt1[i][j] = pPnt1.rotateZ( 67.5 * deg );
202 }
203 if ( ( i == 3 ) || ( i == 7 ) || ( i == 12 ) || ( i == 17 ) || ( i == 23 ) || ( i == 29 ) )
204 {
205 fPnt1[i][0].setX( 10 );
206 fPnt1[i][1].setX( 10 );
207 fPnt1[i][4].setX( 10 );
208 fPnt1[i][5].setX( 10 );
209
210 G4double y0 = fPnt1[i][0].getY() +
211 10 * ( fPnt1[i][3].getY() - fPnt1[i][0].getY() ) / fPnt1[i][3].getX();
212 G4double y1 = fPnt1[i][1].getY() +
213 10 * ( fPnt1[i][2].getY() - fPnt1[i][1].getY() ) / fPnt1[i][2].getX();
214 G4double y4 = fPnt1[i][4].getY() +
215 10 * ( fPnt1[i][7].getY() - fPnt1[i][4].getY() ) / fPnt1[i][7].getX();
216 G4double y5 = fPnt1[i][5].getY() +
217 10 * ( fPnt1[i][6].getY() - fPnt1[i][5].getY() ) / fPnt1[i][6].getX();
218
219 G4double z0 = fPnt1[i][0].getZ() +
220 10 * ( fPnt1[i][3].getZ() - fPnt1[i][0].getZ() ) / fPnt1[i][3].getX();
221 G4double z1 = fPnt1[i][1].getZ() +
222 10 * ( fPnt1[i][2].getZ() - fPnt1[i][1].getZ() ) / fPnt1[i][2].getX();
223 G4double z4 = fPnt1[i][4].getZ() +
224 10 * ( fPnt1[i][7].getZ() - fPnt1[i][4].getZ() ) / fPnt1[i][7].getX();
225 G4double z5 = fPnt1[i][5].getZ() +
226 10 * ( fPnt1[i][6].getZ() - fPnt1[i][5].getZ() ) / fPnt1[i][6].getX();
227
228 fPnt1[i][0].setY( y0 );
229 fPnt1[i][1].setY( y1 );
230 fPnt1[i][4].setY( y4 );
231 fPnt1[i][5].setY( y5 );
232
233 fPnt1[i][0].setZ( z0 );
234 fPnt1[i][1].setZ( z1 );
235 fPnt1[i][4].setZ( z4 );
236 fPnt1[i][5].setZ( z5 );
237 }
238 }
239}
240
241void ExtBesEmcEndGeometry::Exchange( G4int cry1, G4int cry2 ) {
242 G4ThreeVector v;
243 for ( G4int i = 0; i < 8; i++ )
244 {
245 v = fPnt[cry1][i];
246 fPnt[cry1][i] = fPnt[cry2][i];
247 fPnt[cry2][i] = v;
248 }
249}
250
251// reflect x axis to construct sectors on the left
252void ExtBesEmcEndGeometry::ReflectX() {
253 for ( G4int i = 0; i < 35; i++ )
254 {
255 for ( G4int j = 0; j < 8; j++ ) { fPnt1[i][j].setX( -fPnt1[i][j].getX() ); }
256
257 // change the sequence of eight points according to the requirment of IrregBox
258 // point 0<-->1, 2<-->3, 4<-->5, 6<-->7
259 for ( G4int j = 0; j < 8; j++ )
260 {
261 if ( ( j % 2 ) == 0 )
262 {
263 G4ThreeVector v2 = fPnt1[i][j];
264 fPnt1[i][j] = fPnt1[i][j + 1];
265 fPnt1[i][j + 1] = v2;
266 }
267 }
268 }
269}
270
271// subtract the casing
272void ExtBesEmcEndGeometry::ModifyForCasing( G4ThreeVector pos[8], G4int CryNb ) {
273 G4ThreeVector center1( 0, 0, 0 ); // the center of large surface
274 G4ThreeVector center2( 0, 0, 0 ); // the center of small surface
275
276 if ( CryNb == 5 || CryNb == 6 || CryNb == 14 || CryNb == 15 || CryNb == 16 )
277 {
278 center1 = ( pos[0] + pos[1] ) / 2;
279 center2 = ( pos[4] + pos[5] ) / 2;
280 }
281 else if ( CryNb >= 30 && CryNb < 35 )
282 {
283 center1 = ( pos[2] + pos[3] ) / 2;
284 center2 = ( pos[6] + pos[7] ) / 2;
285 }
286 else
287 {
288 center1 = ( pos[0] + pos[1] + pos[2] + pos[3] ) / 4;
289 center2 = ( pos[4] + pos[5] + pos[6] + pos[7] ) / 4;
290 }
291
292 G4double r1 = ( pos[1] - center1 ).r();
293 G4double r2 = ( pos[2] - center1 ).r();
294 G4double r12 = ( pos[1] - pos[2] ).r();
295 G4double theta = acos( ( r2 * r2 + r12 * r12 - r1 * r1 ) / ( 2 * r2 * r12 ) );
296 G4double h = r2 * sin( theta ); // distance from the center to the vertical side
297 G4double t1 = totalThickness / h;
298
299 r1 = ( pos[5] - center2 ).r();
300 r2 = ( pos[6] - center2 ).r();
301 r12 = ( pos[5] - pos[6] ).r();
302 theta = acos( ( r2 * r2 + r12 * r12 - r1 * r1 ) / ( 2 * r2 * r12 ) );
303 h = r2 * sin( theta );
304 G4double t2 = totalThickness / h;
305
306 for ( G4int i = 0; i < 8; i++ )
307 {
308 if ( i < 4 ) { cryPoint[i] = ( 1 - t1 ) * pos[i] + t1 * center1; }
309 else
310 {
311 G4ThreeVector temp = ( 1 - t2 ) * pos[i] + t2 * center2;
312 cryPoint[i] = ( 1 - totalThickness / CrystalLength ) * temp +
313 ( totalThickness / CrystalLength ) * pos[i - 4];
314 }
315 // G4cout<<"casing="<<pos[i]<<"\t"<<"crystal="<<cryPoint[i]<<G4endl;
316 }
317}
318
319G4ThreeVector ExtBesEmcEndGeometry::ComputeDimAndPos( const G4int partId, const G4int numTheta,
320 const G4int numPhi ) {
321 // ComputeParameters();
322 G4int sector = -1, cryNb = -1;
323 G4int leftFlag = 0;
324 G4int downFlag = 0;
325 G4int pentaFlag = 0;
326 G4int flag = 0;
327 G4double A1 = 0, a1 = 0, B1 = 0, b1 = 0, C1 = 0, c1 = 0, D1 = 0, d1 = 0, E1 = 0,
328 e1 = 0; // dimension of pentagonal crystal
329 G4int m_numPhi = 0;
330 G4ThreeVector position( 0, 0, 0 );
331 G4int cryInOneSector =
332 cryNumInOneLayer[numTheta] / 16; // number of crystal in one layer in one sector
333 G4ThreeVector pos[8];
334
335 if ( partId == 2 ) // west end
336 {
337 if ( numPhi >= 0 && numPhi < 8 * cryInOneSector )
338 m_numPhi = 8 * cryInOneSector - 1 - numPhi;
339 else if ( numPhi >= 8 * cryInOneSector && numPhi < 16 * cryInOneSector )
340 m_numPhi = 16 * cryInOneSector - 1 - numPhi;
341 }
342 else // east end
343 m_numPhi = numPhi;
344
345 if ( numPhi >= 4 * cryInOneSector && numPhi < 5 * cryInOneSector ) // crystal in 4th sector
346 {
347 leftFlag = 1;
348 m_numPhi = 8 * cryInOneSector - 1 - numPhi;
349 }
350 else if ( numPhi >= 11 * cryInOneSector && numPhi < 12 * cryInOneSector ) // crystal in 11th
351 // sector
352 {
353 leftFlag = 1;
354 downFlag = 1;
355 m_numPhi = numPhi - 8 * cryInOneSector;
356 }
357 if ( numPhi >= 12 * cryInOneSector && numPhi < 13 * cryInOneSector ) // crystal in 12th
358 // sector
359 {
360 downFlag = 1;
361 m_numPhi = 16 * cryInOneSector - 1 - numPhi;
362 }
363
364 // translate numTheta and numPhi to sector and cryNb
365 G4int cryNbOffset = 0;
366 for ( G4int i = 0; i < numTheta; i++ ) cryNbOffset += cryNumInOneLayer[i] / 16;
367 if ( cryInOneSector ) sector = m_numPhi / cryInOneSector;
368 cryNb = m_numPhi - cryInOneSector * sector + cryNbOffset;
369 sector += 3;
370 if ( sector > 15 && sector <= 18 ) sector -= 16;
371
372 // sector beside the gap
373 if ( sector == 6 )
374 for ( G4int i = 0; i < 8; i++ ) pos[i] = fPnt1[cryNb][i];
375 else
376 for ( G4int i = 0; i < 8; i++ )
377 {
378 pos[i] = fPnt[cryNb][i];
379 pos[i].rotateZ( -67.5 * deg + sector * 22.5 * deg );
380 }
381
382 // crystal dimension
383 ModifyForCasing( pos, cryNb );
384 G4double A = ( cryPoint[0] - cryPoint[3] ).r();
385 G4double a = ( cryPoint[4] - cryPoint[7] ).r();
386 G4double B = ( cryPoint[1] - cryPoint[2] ).r();
387 G4double b = ( cryPoint[5] - cryPoint[6] ).r();
388 G4double C = ( cryPoint[0] - cryPoint[1] ).r();
389 G4double c = ( cryPoint[4] - cryPoint[5] ).r();
390 G4double D = ( cryPoint[2] - cryPoint[3] ).r();
391 G4double d = ( cryPoint[6] - cryPoint[7] ).r();
392
393 // reflect the axis according to the flag
394 for ( G4int i = 0; i < 8; i++ )
395 {
396 pos[i].setZ( pos[i].getZ() + WorldZPosition ); // give the absolute z coordinate
397 if ( leftFlag == 1 ) pos[i].setX( -pos[i].getX() );
398 if ( downFlag == 1 ) pos[i].setY( -pos[i].getY() );
399 if ( partId == 2 )
400 {
401 pos[i].setX( -pos[i].getX() );
402 pos[i].setZ( -pos[i].getZ() );
403 }
404 }
405
406 // compute the position
407 for ( G4int j = 4; j < 8; j++ ) position += pos[j];
408 position /= 4;
409
410 // compute pentagonal crystal
411 for ( G4int i = 0; i < 5; i++ )
412 {
413 if ( cryNb == pentaInOneSector[i] )
414 {
415 pentaFlag = 1;
416 G4ThreeVector penta[8];
417
418 // sector beside the gap
419 if ( sector == 6 )
420 for ( G4int j = 0; j < 8; j++ ) penta[j] = fPnt1[30 + i][j];
421 else
422 for ( G4int j = 0; j < 8; j++ )
423 {
424 penta[j] = fPnt[30 + i][j];
425 penta[j].rotateZ( -67.5 * deg + sector * 22.5 * deg );
426 }
427
428 // crystal dimension
429 ModifyForCasing( penta, 30 + i );
430 A1 = ( cryPoint[1] - cryPoint[2] ).r();
431 a1 = ( cryPoint[5] - cryPoint[6] ).r();
432 B1 = ( cryPoint[1] - cryPoint[0] ).r();
433 b1 = ( cryPoint[5] - cryPoint[4] ).r();
434 C1 = ( cryPoint[0] - cryPoint[3] ).r() + A;
435 c1 = ( cryPoint[4] - cryPoint[7] ).r() + a;
436 D1 = D;
437 d1 = d;
438 E1 = B;
439 e1 = b;
440
441 // reflect the axis according to the flag
442 for ( G4int j = 0; j < 8; j++ )
443 {
444 penta[j].setZ( penta[j].getZ() + WorldZPosition ); // give the absolute z coordinate
445 if ( leftFlag == 1 ) penta[j].setX( -penta[j].getX() );
446 if ( downFlag == 1 ) penta[j].setY( -penta[j].getY() );
447 if ( partId == 2 )
448 {
449 penta[j].setX( -penta[j].getX() );
450 penta[j].setZ( -penta[j].getZ() );
451 }
452 }
453
454 // compute the position
455 for ( G4int j = 4; j < 8; j++ )
456 {
457 if ( j != 0 && j != 4 ) position += pos[j];
458 if ( j == 0 || j == 1 || j == 4 || j == 5 ) position += penta[j];
459 }
460 position /= 10;
461
462 flag = leftFlag + downFlag;
463 if ( flag == 1 )
464 {
465 G4double temp1 = B1;
466 B1 = D1;
467 D1 = temp1;
468 temp1 = A1;
469 A1 = E1;
470 E1 = temp1;
471 temp1 = b1;
472 b1 = d1;
473 d1 = temp1;
474 temp1 = a1;
475 a1 = e1;
476 e1 = temp1;
477 }
478
479 break;
480 }
481 }
482
483 flag = leftFlag + downFlag + partId / 2;
484 if ( flag == 1 || flag == 3 )
485 {
486 G4double temp = C;
487 C = D;
488 D = temp;
489 temp = c;
490 c = d;
491 d = temp;
492 }
493
494 /*G4cout<<"##########################################################################"<<G4endl;
495 G4cout<<"###sector="<<sector<<",cryNb="<<cryNb<<",left flag="<<leftFlag<<",down
496 flag="<<downFlag<<G4endl; G4cout<<"partId="<<partId<<"\t"<<"theta
497 number="<<numTheta<<"\t"<<"phi number="<<numPhi<<G4endl; G4cout<<"crystal point:"<<G4endl;
498 for(G4int i=0;i<8;i++)
499 G4cout<<pos[i]<<G4endl;
500 G4cout<<"crystal position:"<<"\t"<<position<<"\t"<<"phi="<<position.phi()/deg<<G4endl;
501 G4cout<<"crystal dimension:"<<G4endl;
502 if(pentaFlag==1)
503 G4cout<<"A="<<A1<<",a="<<a1<<",B="<<B1<<",b="<<b1<<",C="<<C1<<",c="<<c1<<",D="<<D1<<",d="<<d1<<",E="<<E1<<",e="<<e1<<G4endl;
504 else
505 G4cout<<"A="<<A<<",a="<<a<<",B="<<B<<",b="<<b<<",C="<<C<<",c="<<c<<",D="<<D<<",d="<<d<<G4endl;
506 G4cout<<"##########################################################################"<<G4endl;
507*/
508 return position;
509}
e1
Double_t e1
Definition DataBase/tau_mode.c:7
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
v
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
Definition KarLud.h:35
C
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations C
Definition RRes.h:29
ExtBesEmcEndGeometry::ModifyForCasing
void ModifyForCasing(G4ThreeVector pos[8], G4int CryNb)
Definition ExtBesEmcEndGeometry.cxx:272
ExtBesEmcEndGeometry::ComputeDimAndPos
G4ThreeVector ComputeDimAndPos(const G4int, const G4int, const G4int)
Definition ExtBesEmcEndGeometry.cxx:319
ExtBesEmcEndGeometry::Exchange
void Exchange(G4int cry1, G4int cry2)
Definition ExtBesEmcEndGeometry.cxx:241
ExtBesEmcParameter::GetPentaInOneSector
G4int GetPentaInOneSector(G4int nb)
Definition ExtBesEmcParameter.h:55
ExtBesEmcParameter::GetCryInOneLayer
G4int GetCryInOneLayer(G4int nb)
Definition ExtBesEmcParameter.h:54