1#include "G4ThreeVector.hh"
3#include "BesEmcEndGeometry.hh"
4#include "BesEmcParameter.hh"
34 SectorRmin1 = 489.27*mm;
35 SectorRmax1 = 880.*mm;
36 SectorRmin2 = 621.57*mm;
37 SectorRmax2 = 1113.53*mm;
38 SectorDz = WorldDz-44.*mm;
39 SectorZPosition = -18.*mm;
41 for(G4int i=0;i<6;i++)
43 for(G4int i=0;i<5;i++)
49 totalThickness=fTyvekThickness+fAlThickness+fMylarThickness;
51 G4String ParaPath = getenv(
"EMCSIMROOT");
53 G4Exception(
"BOOST environment not set!");
55 ParaPath +=
"/dat/EmcEndGeometry.dat";
60 for(G4int i=0;i<30;i++)
61 for (G4int j=0;j<24;j++)
63 for(G4int i=0;i<5;i++)
64 for (G4int j=0;j<6;j++)
96 for(G4int i=0;i<30;i++)
98 for (G4int j=0;j<8;j++)
101 G4ThreeVector* pPnt =
new G4ThreeVector(param[i][j],param[i][j+8],param[i][j+16]-WorldZPosition-SectorZPosition);
105 if(i==5||i==6||i==14||i==15||i==16)
108 for(G4int j=0;j<8;j++)
109 fPnt[30+pentaNb][j] = fPnt[i][j];
112 G4double y0,y1,y4,y5;
113 G4ThreeVector v0,v1,v4,v5;
114 y0 = penta[pentaNb][0];
115 y1 = penta[pentaNb][1];
116 y4 = penta[pentaNb][2];
117 y5 = penta[pentaNb][3];
118 v0 = (y0-fPnt[i][3].getY())*(fPnt[i][0]-fPnt[i][3])/(fPnt[i][0].getY()-fPnt[i][3].getY())
120 v1 = (y1-fPnt[i][2].getY())*(fPnt[i][1]-fPnt[i][2])/(fPnt[i][1].getY()-fPnt[i][2].getY())
122 v4 = (y4-fPnt[i][7].getY())*(fPnt[i][4]-fPnt[i][7])/(fPnt[i][4].getY()-fPnt[i][7].getY())
124 v5 = (y5-fPnt[i][6].getY())*(fPnt[i][5]-fPnt[i][6])/(fPnt[i][5].getY()-fPnt[i][6].getY())
128 x1 = penta[pentaNb][4];
129 x5 = penta[pentaNb][5];
130 v1 = (x1-v0.getX())*(v1-v0)/(v1.getX()-v0.getX())+v0;
131 v5 = (x5-v4.getX())*(v5-v4)/(v5.getX()-v4.getX())+v4;
138 fPnt[30+pentaNb][2] = v1;
139 fPnt[30+pentaNb][3] = v0;
140 fPnt[30+pentaNb][6] = v5;
141 fPnt[30+pentaNb][7] = v4;
148 G4ThreeVector temp[35][8];
149 for(G4int i=0;i<35;i++)
151 for (G4int j=0;j<8;j++)
153 temp[i][j]=fPnt[i][j];
154 fPnt[i][j].rotateZ(157.5*deg);
155 fPnt[i][j].setX(-fPnt[i][j].getX());
159 for (G4int j=0;j<8;j++)
163 G4ThreeVector
v=fPnt[i][j];
164 fPnt[i][j]=fPnt[i][3-j];
169 G4ThreeVector
v=fPnt[i][j];
170 fPnt[i][j]=fPnt[i][11-j];
194 for(G4int i=0;i<35;i++)
196 for (G4int j=0;j<8;j++)
198 G4ThreeVector pPnt1 = temp[i][j];
199 fPnt1[i][j] = pPnt1.rotateZ(67.5*deg);
201 if((i==3)||(i==7)||(i==12)||(i==17)||(i==23)||(i==29))
203 fPnt1[i][0].setX(10);
204 fPnt1[i][1].setX(10);
205 fPnt1[i][4].setX(10);
206 fPnt1[i][5].setX(10);
208 G4double y0 = fPnt1[i][0].getY()+10*(fPnt1[i][3].getY()-fPnt1[i][0].getY())/fPnt1[i][3].getX();
209 G4double y1 = fPnt1[i][1].getY()+10*(fPnt1[i][2].getY()-fPnt1[i][1].getY())/fPnt1[i][2].getX();
210 G4double y4 = fPnt1[i][4].getY()+10*(fPnt1[i][7].getY()-fPnt1[i][4].getY())/fPnt1[i][7].getX();
211 G4double y5 = fPnt1[i][5].getY()+10*(fPnt1[i][6].getY()-fPnt1[i][5].getY())/fPnt1[i][6].getX();
213 G4double z0 = fPnt1[i][0].getZ()+10*(fPnt1[i][3].getZ()-fPnt1[i][0].getZ())/fPnt1[i][3].getX();
214 G4double z1 = fPnt1[i][1].getZ()+10*(fPnt1[i][2].getZ()-fPnt1[i][1].getZ())/fPnt1[i][2].getX();
215 G4double z4 = fPnt1[i][4].getZ()+10*(fPnt1[i][7].getZ()-fPnt1[i][4].getZ())/fPnt1[i][7].getX();
216 G4double z5 = fPnt1[i][5].getZ()+10*(fPnt1[i][6].getZ()-fPnt1[i][5].getZ())/fPnt1[i][6].getX();
218 fPnt1[i][0].setY(y0);
219 fPnt1[i][1].setY(y1);
220 fPnt1[i][4].setY(y4);
221 fPnt1[i][5].setY(y5);
223 fPnt1[i][0].setZ(z0);
224 fPnt1[i][1].setZ(z1);
225 fPnt1[i][4].setZ(z4);
226 fPnt1[i][5].setZ(z5);
234 for(G4int i=0;i<8;i++)
237 fPnt[cry1][i] = fPnt[cry2][i];
245 for(G4int i=0;i<8;i++)
248 fPnt1[cry1][i] = fPnt1[cry2][i];
256 for(G4int i=0;i<35;i++)
258 for(G4int j=0;j<8;j++)
260 fPnt1[i][j].setX(-fPnt1[i][j].getX());
264 for (G4int j=0;j<8;j++)
268 G4ThreeVector
v=fPnt1[i][j];
269 fPnt1[i][j]=fPnt1[i][3-j];
274 G4ThreeVector
v=fPnt1[i][j];
275 fPnt1[i][j]=fPnt1[i][11-j];
302 G4ThreeVector center1(0,0,0);
303 G4ThreeVector center2(0,0,0);
304 for(G4int i=0;i<8;i++)
307 for(G4int i=0;i<8;i++)
309 if(i<4) center1+=pos[i];
310 else center2+=pos[i];
315 for(G4int i=0;i<8;i++)
317 if(i<4) zoomPoint[i]=factor*pos[i]+(1-factor)*center1;
318 else zoomPoint[i]=factor*pos[i]+(1-factor)*center2;
325 G4ThreeVector center1=0;
326 G4ThreeVector center2=0;
328 const G4double
dt=1e-5;
330 if(CryNb==5||CryNb==6||CryNb==14||CryNb==15||CryNb==16)
332 center1 = (pos[0]+pos[1])*(1-
dt)/2+(pos[2]+pos[3])*
dt/2;
333 center2 = (pos[4]+pos[5])*(1-
dt)/2+(pos[6]+pos[7])*
dt/2;
335 else if(CryNb>=30&&CryNb<35)
337 center1 = (pos[0]+pos[1])*
dt/2+(pos[2]+pos[3])*(1-
dt)/2;
338 center2 = (pos[4]+pos[5])*
dt/2+(pos[6]+pos[7])*(1-
dt)/2;
342 center1 = (pos[0]+pos[1]+pos[2]+pos[3])/4;
343 center2 = (pos[4]+pos[5]+pos[6]+pos[7])/4;
346 G4double r1=(pos[1]-center1).r();
347 G4double r2=(pos[2]-center1).r();
348 G4double r12=(pos[1]-pos[2]).r();
349 G4double theta=acos((r2*r2+r12*r12-r1*r1)/(2*r2*r12));
350 G4double h=r2*
sin(theta);
351 G4double t1=totalThickness/h;
353 r1=(pos[5]-center2).r();
354 r2=(pos[6]-center2).r();
355 r12=(pos[5]-pos[6]).r();
356 theta=acos((r2*r2+r12*r12-r1*r1)/(2*r2*r12));
358 G4double t2=totalThickness/h;
360 for(G4int i=0;i<8;i++)
364 cryPoint[i] = (1-t1)*pos[i]+t1*center1;
368 G4ThreeVector temp = (1-t2)*pos[i]+t2*center2;
369 cryPoint[i] = (1-totalThickness/CrystalLength)*temp+(totalThickness/CrystalLength)*pos[i-4];
376(
const G4int partId,
const G4int numTheta,
const G4int numPhi)
379 G4int sector=-1, cryNb=-1;
384 G4double A1=0,a1=0,B1=0,b1=0,C1=0,c1=0,D1=0,d1=0,E1=0,
e1=0;
386 G4ThreeVector position=0;
387 G4int cryInOneSector = cryNumInOneLayer[numTheta]/16;
388 G4ThreeVector pos[8];
392 if(numPhi>=0&&numPhi<8*cryInOneSector)
393 m_numPhi=8*cryInOneSector-1-numPhi;
394 else if(numPhi>=8*cryInOneSector&&numPhi<16*cryInOneSector)
395 m_numPhi=16*cryInOneSector-1-numPhi;
400 if(numPhi>=4*cryInOneSector&&numPhi<5*cryInOneSector)
403 m_numPhi=8*cryInOneSector-1-numPhi;
405 else if(numPhi>=11*cryInOneSector&&numPhi<12*cryInOneSector)
409 m_numPhi=numPhi-8*cryInOneSector;
411 if(numPhi>=12*cryInOneSector&&numPhi<13*cryInOneSector)
414 m_numPhi=16*cryInOneSector-1-numPhi;
418 G4int cryNbOffset = 0;
419 for(G4int i=0;i<numTheta;i++)
420 cryNbOffset += cryNumInOneLayer[i]/16;
422 sector = m_numPhi/cryInOneSector;
423 cryNb = m_numPhi-cryInOneSector*sector+cryNbOffset;
425 if(sector>15&§or<=18)
430 for(G4int i=0;i<8;i++)
431 pos[i]=fPnt1[cryNb][i];
433 for(G4int i=0;i<8;i++)
435 pos[i]=fPnt[cryNb][i];
436 pos[i].rotateZ(-67.5*deg+sector*22.5*deg);
442 G4double A = (cryPoint[0]-cryPoint[3]).r();
443 G4double a = (cryPoint[4]-cryPoint[7]).r();
444 G4double B = (cryPoint[1]-cryPoint[2]).r();
445 G4double
b = (cryPoint[5]-cryPoint[6]).r();
446 G4double
C = (cryPoint[0]-cryPoint[1]).r();
447 G4double c = (cryPoint[4]-cryPoint[5]).r();
448 G4double D = (cryPoint[2]-cryPoint[3]).r();
449 G4double d = (cryPoint[6]-cryPoint[7]).r();
452 for(G4int i=0;i<8;i++)
454 pos[i].setZ(pos[i].getZ()+WorldZPosition+SectorZPosition);
456 pos[i].setX(-pos[i].getX());
458 pos[i].setY(-pos[i].getY());
461 pos[i].setX(-pos[i].getX());
462 pos[i].setZ(-pos[i].getZ());
467 for(G4int j=4;j<8;j++)
472 for(G4int i=0;i<5;i++)
474 if(cryNb==pentaInOneSector[i])
477 G4ThreeVector penta[8];
481 for(G4int j=0;j<8;j++)
482 penta[j]=fPnt1[30+i][j];
484 for(G4int j=0;j<8;j++)
486 penta[j]=fPnt[30+i][j];
487 penta[j].rotateZ(-67.5*deg+sector*22.5*deg);
492 A1 = (cryPoint[1]-cryPoint[2]).r();
493 a1 = (cryPoint[5]-cryPoint[6]).r();
494 B1 = (cryPoint[1]-cryPoint[0]).r();
495 b1 = (cryPoint[5]-cryPoint[4]).r();
496 C1 = (cryPoint[0]-cryPoint[3]).r()+A;
497 c1 = (cryPoint[4]-cryPoint[7]).r()+a;
504 for(G4int j=0;j<8;j++)
506 penta[j].setZ(penta[j].getZ()+WorldZPosition+SectorZPosition);
508 penta[j].setX(-penta[j].getX());
510 penta[j].setY(-penta[j].getY());
513 penta[j].setX(-penta[j].getX());
514 penta[j].setZ(-penta[j].getZ());
519 for(G4int j=4;j<8;j++)
523 if(j==0||j==1||j==4||j==5)
524 position += penta[j];
528 flag = leftFlag+downFlag;
531 G4double temp1 = B1; B1 = D1; D1 = temp1;
532 temp1 = A1; A1 = E1; E1 = temp1;
533 temp1 = b1; b1 = d1; d1 = temp1;
534 temp1 = a1; a1 =
e1;
e1 = temp1;
541 flag = leftFlag+downFlag+partId/2;
544 G4double temp =
C;
C = D; D = temp;
545 temp = c; c = d; d = temp;
double sin(const BesAngle a)
**********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
***************************************************************************************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
void Exchange(G4int cry1, G4int cry2)
void ExchangeSector7(G4int cry1, G4int cry2)
void ModifyForCasing(G4ThreeVector pos[8], G4int CryNb)
void Zoom(const G4ThreeVector pos[8], const G4double factor)
G4ThreeVector ComputeDimAndPos(const G4int, const G4int, const G4int)
G4double GetWorldZPosition()
G4double GetCrystalLength()
G4int GetCryInOneLayer(G4int nb)
static BesEmcParameter & GetInstance()
G4double GetTyvekThickness()
G4double GetMylarThickness()
G4int GetPentaInOneSector(G4int nb)
G4double GetAlThickness()