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