00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "G4ThreeVector.hh"
00011
00012 #include "BesEmcGeometry.hh"
00013 #include "BesEmcParameter.hh"
00014 #include "G4ios.hh"
00015
00016 BesEmcGeometry::BesEmcGeometry()
00017 {
00018 verboseLevel = 0;
00019 }
00020
00021 BesEmcGeometry::~BesEmcGeometry()
00022 {;}
00023
00024 void BesEmcGeometry::ReadEMCParameters()
00025 {
00026
00027 BesEmcParameter& emcPara=BesEmcParameter::GetInstance();
00028
00029 TaperRingDz = emcPara.GetTaperRingDz();
00030 TaperRingThickness1 = emcPara.GetTaperRingThickness1();
00031 TaperRingThickness2 = emcPara.GetTaperRingThickness2();
00032 TaperRingThickness3 = emcPara.GetTaperRingThickness3();
00033 TaperRingTheta = emcPara.GetTaperRingTheta()*deg;
00034 TaperRingInnerLength = emcPara.GetTaperRingInnerLength();
00035 TaperRingOuterLength = emcPara.GetTaperRingOuterLength();
00036
00037 BSCRmin = emcPara.GetBSCRmin();
00038 BSCDz = emcPara.GetBSCDz()-TaperRingThickness3;
00039
00040 BSCRmin1 = emcPara.GetBSCRmin1();
00041 BSCRmax1 = emcPara.GetBSCRmax1();
00042 BSCRmin2 = emcPara.GetBSCRmin2();
00043
00044 BSCDz1 = emcPara.GetBSCDz1()-TaperRingThickness3;
00045
00046 BSCAngleRotat = emcPara.GetBSCAngleRotat()*deg;
00047 BSCNbPhi = emcPara.GetBSCNbPhi();
00048 BSCNbTheta = emcPara.GetBSCNbTheta();
00049
00050 BSCCryLength = emcPara.GetCrystalLength();
00051 BSCCryLength1 = emcPara.GetCrystalLength1();
00052 BSCYFront0 = emcPara.GetBSCYFront0();
00053 BSCYFront = emcPara.GetBSCYFront();
00054 BSCYFront1 = emcPara.GetBSCYFront1();
00055 BSCPosition0 = emcPara.GetBSCPosition0();
00056 BSCPosition1 = emcPara.GetBSCPosition1();
00057
00058 fTyvekThickness = emcPara.GetTyvekThickness();
00059 fAlThickness = emcPara.GetAlThickness();
00060 fMylarThickness = emcPara.GetMylarThickness();
00061
00062 rearBoxLength = emcPara.GetRearBoxLength();
00063 rearBoxDz = emcPara.GetRearBoxDz();
00064
00065 HangingPlateDz = emcPara.GetHangingPlateDz();
00066 OCGirderAngle = emcPara.GetOCGirderAngle()*deg;
00067
00068 rearCasingThickness = emcPara.GetRearCasingThickness();
00069
00070 orgGlassLengthX = emcPara.GetOrgGlassLengthX();
00071 orgGlassLengthY = emcPara.GetOrgGlassLengthY();
00072 orgGlassLengthZ = emcPara.GetOrgGlassLengthZ();
00073
00074 PDLengthX = emcPara.GetPDLengthX();
00075 PDLengthY = emcPara.GetPDLengthY();
00076 PDLengthZ = emcPara.GetPDLengthZ();
00077
00078 AlPlateDz = emcPara.GetAlPlateDz();
00079 PABoxDz = emcPara.GetPABoxDz();
00080 PABoxThickness = emcPara.GetPABoxThickness();
00081
00082 cableDr = emcPara.GetCableDr();
00083 waterPipeDr = emcPara.GetWaterPipeDr();
00084 waterPipeThickness= emcPara.GetWaterPipeThickness();
00085
00086 SPBarThickness = emcPara.GetSPBarThickness();
00087 SPBarThickness1 = emcPara.GetSPBarThickness1();
00088 SPBarwidth = emcPara.GetSPBarwidth();
00089
00090 EndRingDz = emcPara.GetEndRingDz();
00091 EndRingDr = emcPara.GetEndRingDr();
00092 EndRingRmin = emcPara.GetEndRingRmin();
00093
00094 for(G4int i=0;i<BSCNbTheta*6;i++)
00095 {
00096 zHalfLength[i] =0;
00097 thetaAxis[i] =0;
00098 phiAxis [i] =0;
00099 yHalfLength1[i] =0;
00100 xHalfLength1[i] =0;
00101 xHalfLength2[i] =0;
00102 tanAlpha1[i] =0;
00103 yHalfLength2[i] =0;
00104 xHalfLength3[i] =0;
00105 xHalfLength4[i] =0;
00106 tanAlpha2[i] =0;
00107 thetaPosition[i]=0;
00108 xPosition[i] =0;
00109 yPosition[i] =0;
00110 zPosition[i] =0;
00111 }
00112 }
00113
00114 void BesEmcGeometry::PrintEMCParameters()
00115 {
00116
00117 ;
00118 }
00119
00120 void BesEmcGeometry::ComputeEMCParameters()
00121 {
00122 ReadEMCParameters();
00123
00124
00125 G4double theta0;
00126 G4int i;
00127
00128
00129
00130
00131
00132 const G4double delta=1.*mm;
00133 TaperRingRmin1 = BSCRmax1-TaperRingThickness1;
00134 TaperRingRmin2 = TaperRingRmin1+TaperRingDz*tan(TaperRingTheta);
00135 TaperRingDr = TaperRingThickness2/cos(TaperRingTheta);
00136 TaperRingOuterLength1 = TaperRingOuterLength+TaperRingThickness3*tan(TaperRingTheta);
00137
00138 BSCRmax2 = TaperRingRmin2+TaperRingDr-TaperRingThickness3*tan(TaperRingTheta);
00139 BSCRmax = BSCRmin+33.8*cm;
00140 BSCPhiDphi= 360.*deg/BSCNbPhi;
00141 BSCPhiSphi= -BSCPhiDphi/2.;
00142 BSCPhiDz = BSCDz;
00143 BSCPhiRmax= BSCRmax-0.5*cm;
00144 BSCPhiRmin= BSCRmin/sin(90.*deg+BSCPhiDphi/2.)
00145 *sin(90.*deg+BSCPhiDphi/2.-BSCAngleRotat);
00146 SPBarDphi=2*atan(SPBarwidth/(2*BSCRmax));
00147
00148
00149 zHalfLength[0] = BSCCryLength/2.;
00150 yHalfLength1[0]= BSCYFront0/2.;
00151 xHalfLength1[0]= BSCPhiRmin*tan(BSCPhiDphi)/2.;
00152 xHalfLength2[0]= xHalfLength1[0];
00153
00154
00155 G4double rminprojected=BSCRmin*cos(BSCAngleRotat);
00156 rminprojected=BSCRmin;
00157
00158 yHalfLength2[0]= yHalfLength1[0]
00159 +(BSCYFront0-BSCPosition0)/rminprojected*BSCCryLength/2.;
00160 xHalfLength3[0]= xHalfLength1[0]*(BSCPhiRmin+BSCCryLength)/BSCPhiRmin;
00161 xHalfLength4[0]= xHalfLength2[0]*(BSCPhiRmin+BSCCryLength)/BSCPhiRmin;
00162
00163 thetaPosition[0]= 90.*deg;
00164 xPosition[0] = BSCPhiRmin+BSCCryLength/2.;
00165 yPosition[0] =
00166 -(xHalfLength1[0]+xHalfLength3[0]+xHalfLength2[0]+xHalfLength4[0])/4.;
00167 zPosition[0] = (yHalfLength1[0]+yHalfLength2[0])/2.;
00168
00169 theta0= 90.*deg-atan((BSCYFront0-BSCPosition0)/rminprojected);
00170 if(verboseLevel>1)
00171 {
00172 G4cout << "------------------------------------>1"<< G4endl
00173 << "The direction of crystal is:"
00174 << theta0/deg <<"(deg)."<< G4endl;
00175 }
00176
00177 rearBoxPosX[0] = xPosition[0]+BSCCryLength/2.+rearBoxDz/2.;
00178 rearBoxPosY[0] = -(xHalfLength3[0]+xHalfLength4[0])/2.;
00179 rearBoxPosZ[0] = yHalfLength2[0];
00180
00181 OCGirderRmin1[0] = rearBoxPosX[0]+rearBoxDz/2.+delta;
00182 OCGirderRmin2[0] = rearBoxPosX[0]+rearBoxDz/2.+delta;
00183 OCGirderDz[0] = rearBoxLength;
00184 OCGirderPosZ[0] = rearBoxLength/2;
00185
00186 cableLength[0] = BSCDz-rearBoxPosZ[0];
00187 cablePosX[0] = BSCPhiRmax-cableDr-2.*mm;
00188 cablePosY[0] = rearBoxPosY[0]-3*cableDr;
00189 cablePosZ[0] = BSCDz-cableLength[0]/2;
00190
00191
00192 zHalfLength[1]= BSCCryLength/2.;
00193 yHalfLength1[1]= BSCYFront0/2.;
00194 G4double deltaR= BSCYFront0*cos(theta0);
00195 xHalfLength1[1]= xHalfLength1[0];
00196 xHalfLength2[1]= xHalfLength1[1]*(BSCPhiRmin+deltaR)/BSCPhiRmin;
00197 yHalfLength2[1]= yHalfLength1[1]
00198 +tan(theta0-atan(rminprojected/(BSCYFront0*(1.+1./sin(theta0))-
00199 BSCPosition1)))*BSCCryLength/2.;
00200 deltaR=deltaR+BSCCryLength*sin(theta0);
00201 xHalfLength4[1]= xHalfLength1[1]*(BSCPhiRmin+deltaR)/BSCPhiRmin;
00202 deltaR=deltaR-yHalfLength2[1]*2.*cos(theta0);
00203 xHalfLength3[1]= xHalfLength1[1]*(BSCPhiRmin+deltaR)/BSCPhiRmin;
00204
00205 thetaPosition[1]= theta0;
00206 xPosition[1] = BSCPhiRmin+BSCCryLength/2.*sin(theta0)
00207 + (3.*yHalfLength1[1]-yHalfLength2[1])/2.*cos(theta0);
00208 yPosition[1] =
00209 -(xHalfLength1[1]+xHalfLength3[1]+xHalfLength2[1]+xHalfLength4[1])/4.;
00210 zPosition[1] = yHalfLength1[0]*2.
00211 + (yHalfLength1[1]*2./tan(theta0)+BSCCryLength/2.)*cos(theta0)
00212 + (yHalfLength1[1]+yHalfLength2[1])/2.*sin(theta0);
00213
00214 rearBoxPosX[1] = xPosition[1]+(BSCCryLength/2.+rearBoxDz/2.)*sin(theta0);
00215 rearBoxPosY[1] = -(xHalfLength3[1]+xHalfLength4[1])/2.;
00216 rearBoxPosZ[1] = zPosition[1]+(BSCCryLength/2.+rearBoxDz/2.)*cos(theta0);
00217
00218 OCGirderRmin1[1] = rearBoxPosX[1]+rearBoxDz*sin(theta0)/2.+rearBoxLength*cos(theta0)/2+delta;
00219 OCGirderRmin2[1] = rearBoxPosX[1]+rearBoxDz*sin(theta0)/2.-rearBoxLength*cos(theta0)/2+delta;
00220 OCGirderDz[1] = rearBoxLength*sin(theta0);
00221 OCGirderPosZ[1] = rearBoxPosZ[1]+rearBoxDz*cos(theta0)/2.;
00222
00223 cableLength[1] = BSCDz-rearBoxPosZ[1];
00224 cablePosX[1] = cablePosX[0];
00225 cablePosY[1] = cablePosY[0]+2*cableDr;
00226 cablePosZ[1] = BSCDz-cableLength[1]/2;
00227
00228 theta0= theta0-atan((yHalfLength2[1]-yHalfLength1[1])*2./BSCCryLength);
00229
00230 for(i=2;i<BSCNbTheta;i++)
00231 {
00232 if(verboseLevel>1)
00233 {
00234 G4cout << "------------------------------------>"<<i<< G4endl
00235 << "The direction of crystal is:"
00236 << theta0/deg <<"(deg)." << G4endl;
00237 }
00238 G4double CryLength;
00239 if(i==BSCNbTheta-1)
00240 {
00241 CryLength=BSCCryLength1;
00242 yHalfLength1[i]=BSCYFront1/2.;
00243 } else {
00244
00245 CryLength=BSCCryLength;
00246 yHalfLength1[i]=BSCYFront/2.;
00247 }
00248 zHalfLength[i]=CryLength/2;
00249
00250 deltaR=yHalfLength1[i]*2.*cos(theta0);
00251 xHalfLength1[i]=xHalfLength1[0];
00252 xHalfLength2[i]=xHalfLength1[i]/BSCPhiRmin*(BSCPhiRmin+deltaR);
00253 yHalfLength2[i]=yHalfLength1[i]
00254 *(1.+CryLength/(rminprojected/sin(theta0)
00255 +yHalfLength1[i]*2./tan(theta0)));
00256 deltaR=deltaR+CryLength*sin(theta0);
00257 xHalfLength4[i]=xHalfLength1[i]/BSCPhiRmin*(BSCPhiRmin+deltaR);
00258 deltaR=deltaR-yHalfLength2[i]*2.*cos(theta0);
00259 xHalfLength3[i]=xHalfLength1[i]/BSCPhiRmin*(BSCPhiRmin+deltaR);
00260
00261 thetaPosition[i]=theta0;
00262 xPosition[i]=BSCPhiRmin+zHalfLength[i]*sin(theta0)
00263 + (3.*yHalfLength1[i]-yHalfLength2[i])/2.*cos(theta0);
00264 yPosition[i]=
00265 -(xHalfLength1[i]+xHalfLength3[i]+xHalfLength2[i]+xHalfLength4[i])/4.;
00266 zPosition[i]=BSCPosition1+rminprojected/tan(theta0)
00267 +(2.*yHalfLength1[i]/tan(theta0)+zHalfLength[i])*cos(theta0)
00268 +(yHalfLength1[i]+yHalfLength2[i])/2.*sin(theta0);
00269
00270 rearBoxPosX[i] = xPosition[i]+(CryLength/2.+rearBoxDz/2.)*sin(theta0);
00271 rearBoxPosY[i] = -(xHalfLength3[i]+xHalfLength4[i])/2.;
00272 rearBoxPosZ[i] = zPosition[i]+(CryLength/2.+rearBoxDz/2.)*cos(theta0);
00273
00274 OCGirderRmin1[i] = rearBoxPosX[i]+rearBoxDz*sin(theta0)/2.+rearBoxLength*cos(theta0)/2+delta;
00275 OCGirderRmin2[i] = rearBoxPosX[i]+rearBoxDz*sin(theta0)/2.-rearBoxLength*cos(theta0)/2+delta;
00276 OCGirderDz[i] = rearBoxLength*sin(theta0);
00277 OCGirderPosZ[i] = rearBoxPosZ[i]+rearBoxDz*cos(theta0)/2.;
00278
00279 G4int xCable,yCable;
00280 xCable = i/4;
00281 yCable = i-4*(G4int)(i/4);
00282 cableLength[i] = BSCDz-(rearBoxPosZ[i]+rearBoxDz/2.*cos(theta0));
00283 cablePosX[i] = cablePosX[0]-xCable*2*cableDr;
00284 cablePosY[i] = cablePosY[0]+yCable*2*cableDr;
00285 cablePosZ[i] = BSCDz-cableLength[i]/2;
00286
00287 theta0=theta0-atan(2.*yHalfLength1[i]/(rminprojected/sin(theta0)
00288 +2.*yHalfLength1[i]/tan(theta0)));
00289
00290 }
00291 thetaPosition[BSCNbTheta]=theta0;
00292 if(verboseLevel>1)
00293 {
00294 G4cout << "------------------------------------>"<<i<< G4endl
00295 << "The direction of crystal is:"
00296 << theta0/deg <<"(deg)." << G4endl;
00297 }
00298
00299 G4double oop;
00300 G4double dx=0.001*mm;
00301 for(i=0;i<BSCNbTheta;i++)
00302 {
00303 xHalfLength1[i]-=dx;
00304 xHalfLength2[i]-=dx;
00305 xHalfLength3[i]-=dx;
00306 xHalfLength4[i]-=dx;
00307 yHalfLength1[i]-=dx;
00308 yHalfLength2[i]-=dx;
00309
00310 G4double CryLength;
00311 if(i==BSCNbTheta-1)
00312 {
00313 CryLength=BSCCryLength1;
00314 } else {
00315 CryLength=BSCCryLength;
00316 }
00317 oop=sqrt((yHalfLength2[i]-yHalfLength1[i])*(yHalfLength2[i]
00318 -yHalfLength1[i])
00319 +(xHalfLength3[i]+xHalfLength4[i]-xHalfLength1[i]
00320 -xHalfLength2[i])*(xHalfLength3[i]+xHalfLength4[i]
00321 -xHalfLength1[i]-xHalfLength2[i])/4);
00322 thetaAxis[i]=atan(oop/CryLength);
00323 phiAxis[i] =180.*deg+atan((yHalfLength2[i]-yHalfLength1[i])
00324 /(xHalfLength3[i]+xHalfLength4[i]
00325 -xHalfLength1[i]-xHalfLength2[i])*2.);
00326 tanAlpha2[i]=-(xHalfLength4[i]-xHalfLength3[i])/yHalfLength2[i]/2.;
00327 tanAlpha1[i]=-(xHalfLength2[i]-xHalfLength1[i])/yHalfLength1[i]/2.;
00328
00329
00330
00331
00332 }
00333
00334 }
00335
00336 void BesEmcGeometry::ModifyForCasing()
00337 {
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356 G4double totalThickness=fTyvekThickness+fAlThickness+fMylarThickness;
00357 G4double delta=0.,angle1=0.*deg,angle2=0.*deg;
00358 G4double oop;
00359
00360 G4double rminprojected=BSCRmin*cos(BSCAngleRotat);
00361 rminprojected=BSCRmin;
00362
00363
00364 G4int i;
00365 for(i=0;i<BSCNbTheta;i++)
00366 {
00367 G4double CryLength;
00368 if(i==BSCNbTheta-1)
00369 {
00370 CryLength=BSCCryLength1;
00371 } else {
00372 CryLength=BSCCryLength;
00373 }
00374 zHalfLength[BSCNbTheta+i]=totalThickness/2.;
00375 yHalfLength1[BSCNbTheta+i]=yHalfLength1[i];
00376 yHalfLength2[BSCNbTheta+i]=yHalfLength1[i]
00377 +(yHalfLength2[i]-yHalfLength1[i])*totalThickness/CryLength;
00378 xHalfLength1[BSCNbTheta+i]=xHalfLength1[i];
00379 xHalfLength2[BSCNbTheta+i]=xHalfLength2[i];
00380 xHalfLength1[BSCNbTheta*2+i]=xHalfLength3[BSCNbTheta+i]=
00381 xHalfLength1[i]*(CryLength-totalThickness)/CryLength
00382 +xHalfLength3[i]*totalThickness/CryLength;
00383 xHalfLength2[BSCNbTheta*4+i]=xHalfLength4[BSCNbTheta+i]=
00384 xHalfLength2[i]*(CryLength-totalThickness)/CryLength
00385 +xHalfLength4[i]*totalThickness/CryLength;
00386
00387 zHalfLength[BSCNbTheta*5+i]=zHalfLength[BSCNbTheta*4+i]=
00388 zHalfLength[BSCNbTheta*3+i]=zHalfLength[BSCNbTheta*2+i]=
00389 zHalfLength[i]-totalThickness/2.;
00390
00391 yHalfLength2[BSCNbTheta*2+i]=yHalfLength1[BSCNbTheta*2+i]=
00392 totalThickness/cos(thetaPosition[i+1]-thetaPosition[i])/2.;
00393 xHalfLength3[BSCNbTheta*2+i]=xHalfLength3[i];
00394 xHalfLength4[BSCNbTheta*2+i]=xHalfLength3[i]
00395 +(xHalfLength4[i]-xHalfLength3[i])*yHalfLength2[BSCNbTheta*2+i]
00396 /yHalfLength2[i];
00397 xHalfLength2[BSCNbTheta*2+i]=xHalfLength3[BSCNbTheta+i]
00398 +(xHalfLength4[BSCNbTheta+i]-xHalfLength3[BSCNbTheta+i])
00399 *yHalfLength1[BSCNbTheta*2+i]/yHalfLength2[BSCNbTheta*1+i];
00400
00401 yHalfLength2[BSCNbTheta*4+i]=yHalfLength1[BSCNbTheta*4+i]=
00402 totalThickness/2.;
00403 xHalfLength4[BSCNbTheta*4+i]=xHalfLength4[i];
00404 xHalfLength3[BSCNbTheta*4+i]=xHalfLength4[i]
00405 -(xHalfLength4[i]-xHalfLength3[i])*yHalfLength2[BSCNbTheta*4+i]
00406 /yHalfLength2[i];
00407 xHalfLength1[BSCNbTheta*4+i]=xHalfLength4[BSCNbTheta+i]
00408 -(xHalfLength4[BSCNbTheta+i]-xHalfLength3[BSCNbTheta+i])
00409 *yHalfLength1[BSCNbTheta*4+i]/yHalfLength2[BSCNbTheta*1+i];
00410
00411 delta=totalThickness/2.+yHalfLength1[BSCNbTheta*2+i];
00412 angle1=atan(yHalfLength2[i]/(xHalfLength4[i]-xHalfLength3[i]));
00413 angle2=atan(2.*(xHalfLength4[i]-xHalfLength2[i])*sin(angle1)
00414 /CryLength);
00415
00416 yHalfLength1[BSCNbTheta*5+i]=yHalfLength1[BSCNbTheta*3+i]=
00417 yHalfLength1[i]-delta;
00418 yHalfLength2[BSCNbTheta*5+i]=yHalfLength2[BSCNbTheta*3+i]=
00419 yHalfLength2[i]-delta;
00420 xHalfLength4[BSCNbTheta*3+i]=xHalfLength3[BSCNbTheta*3+i]=
00421 xHalfLength2[BSCNbTheta*3+i]=xHalfLength1[BSCNbTheta*3+i]=
00422 totalThickness/cos(angle2)/sin(angle1)/2.;
00423 xHalfLength4[BSCNbTheta*5+i]=xHalfLength3[BSCNbTheta*5+i]=
00424 xHalfLength2[BSCNbTheta*5+i]=xHalfLength1[BSCNbTheta*5+i]=
00425 totalThickness/2.;
00426
00427 zHalfLength[i]=zHalfLength[i]-totalThickness/2.;
00428 yHalfLength1[i]=yHalfLength1[i]-delta;
00429 yHalfLength2[i]=yHalfLength2[i]-delta;
00430 delta=totalThickness*(1.+1./cos(angle2)/sin(angle1))/2.;
00431 xHalfLength1[i]=xHalfLength1[i]-delta;
00432 xHalfLength2[i]=xHalfLength2[i]-delta;
00433 xHalfLength3[i]=xHalfLength3[i]-delta;
00434 xHalfLength4[i]=xHalfLength4[i]-delta;
00435
00436 oop=sqrt((yHalfLength2[i]-yHalfLength1[i])*(yHalfLength2[i]
00437 -yHalfLength1[i])
00438 +(xHalfLength3[i]+xHalfLength4[i]-xHalfLength1[i]
00439 -xHalfLength2[i])*(xHalfLength3[i]+xHalfLength4[i]
00440 -xHalfLength1[i]-xHalfLength2[i])/4);
00441 thetaAxis[i]=atan(oop/CryLength);
00442
00443 phiAxis[i] =180.*deg+atan((yHalfLength2[i]-yHalfLength1[i])
00444 /(xHalfLength3[i]+xHalfLength4[i]
00445 -xHalfLength1[i]-xHalfLength2[i])*2.);
00446
00447 oop=sqrt((yHalfLength2[BSCNbTheta+i]-yHalfLength1[BSCNbTheta+i])
00448 *(yHalfLength2[BSCNbTheta+i]-yHalfLength1[BSCNbTheta+i])
00449 +(xHalfLength3[BSCNbTheta+i]+xHalfLength4[BSCNbTheta+i]
00450 -xHalfLength1[BSCNbTheta+i]-xHalfLength2[BSCNbTheta+i])
00451 *(xHalfLength3[BSCNbTheta+i]+xHalfLength4[BSCNbTheta+i]
00452 -xHalfLength1[BSCNbTheta+i]-xHalfLength2[BSCNbTheta+i])/4);
00453 thetaAxis[BSCNbTheta+i]=atan(oop/totalThickness);
00454 phiAxis [BSCNbTheta+i]=
00455 -atan((yHalfLength2[BSCNbTheta+i]-yHalfLength1[BSCNbTheta+i])
00456 /(xHalfLength3[BSCNbTheta+i]+xHalfLength4[BSCNbTheta+i]
00457 -xHalfLength1[BSCNbTheta+i]-xHalfLength2[BSCNbTheta+i])*2.);
00458
00459 oop=sqrt((yHalfLength2[i]-yHalfLength1[i])*(yHalfLength2[i]
00460 -yHalfLength1[i])*4
00461 +(xHalfLength3[BSCNbTheta*2+i]+xHalfLength4[BSCNbTheta*2+i]
00462 -xHalfLength1[BSCNbTheta*2+i]-xHalfLength2[BSCNbTheta*2+i])
00463 *(xHalfLength3[BSCNbTheta*2+i]+xHalfLength4[BSCNbTheta*2+i]
00464 -xHalfLength1[BSCNbTheta*2+i]-xHalfLength2[BSCNbTheta*2+i])
00465 /4);
00466 thetaAxis[BSCNbTheta*2+i]=atan(oop/(zHalfLength[BSCNbTheta*2+i]*2));
00467 phiAxis [BSCNbTheta*2+i]=
00468 -atan((yHalfLength2[i]-yHalfLength1[i])
00469 /(xHalfLength3[BSCNbTheta*2+i]+xHalfLength4[BSCNbTheta*2+i]
00470 -xHalfLength1[BSCNbTheta*2+i]-xHalfLength2[BSCNbTheta*2+i])*4);
00471
00472 oop=sqrt((yHalfLength2[i]-yHalfLength1[i])*(yHalfLength2[i]
00473 -yHalfLength1[i])*4
00474 +(xHalfLength4[i]-xHalfLength2[i])
00475 *(xHalfLength4[i]-xHalfLength2[i])*4);
00476 thetaAxis[BSCNbTheta*3+i]=atan(oop/(zHalfLength[BSCNbTheta*3+i]*2));
00477 phiAxis [BSCNbTheta*3+i]=-atan((yHalfLength2[i]-yHalfLength1[i])
00478 /(xHalfLength4[i]-xHalfLength2[i]));
00479
00480 thetaAxis[BSCNbTheta*4+i]=
00481 atan((xHalfLength4[BSCNbTheta*4+i]+xHalfLength3[BSCNbTheta*4+i]
00482 -xHalfLength2[BSCNbTheta*4+i]-xHalfLength1[BSCNbTheta*4+i])/2.
00483 /(zHalfLength[BSCNbTheta*4+i]*2));
00484 phiAxis [BSCNbTheta*4+i]=0;
00485
00486 thetaAxis[BSCNbTheta*5+i]=atan((xHalfLength3[BSCNbTheta*5+i]
00487 -xHalfLength1[BSCNbTheta*5+i])
00488 /(zHalfLength[BSCNbTheta*5+i]*2));
00489 phiAxis [BSCNbTheta*5+i]=-90.*deg;
00490
00491 tanAlpha2[BSCNbTheta+i]=tanAlpha1[BSCNbTheta+i]=tanAlpha1[i]=
00492 -(xHalfLength2[i]-xHalfLength1[i])/yHalfLength1[i]/2.;
00493 tanAlpha1[BSCNbTheta*2+i]=(xHalfLength2[BSCNbTheta*2+i]
00494 -xHalfLength1[BSCNbTheta*2+i])
00495 /yHalfLength1[BSCNbTheta*2+i]/2.;
00496 tanAlpha1[BSCNbTheta*3+i]=tanAlpha1[i]*2.;
00497 tanAlpha1[BSCNbTheta*4+i]=(xHalfLength2[BSCNbTheta*4+i]
00498 -xHalfLength1[BSCNbTheta*4+i])
00499 /yHalfLength1[BSCNbTheta*4+i]/2.;
00500 tanAlpha1[BSCNbTheta*5+i]=(xHalfLength2[BSCNbTheta*5+i]
00501 -xHalfLength1[BSCNbTheta*5+i])
00502 /yHalfLength1[BSCNbTheta*5+i]/2.;
00503
00504 tanAlpha2[i]=-(xHalfLength4[i]-xHalfLength3[i])/yHalfLength2[i]/2.;
00505
00506 tanAlpha2[BSCNbTheta*2+i]=(xHalfLength4[BSCNbTheta*2+i]
00507 -xHalfLength3[BSCNbTheta*2+i])
00508 /yHalfLength2[BSCNbTheta*2+i]/2.;
00509 tanAlpha2[BSCNbTheta*3+i]=tanAlpha2[i]*2.;
00510 tanAlpha2[BSCNbTheta*4+i]=(xHalfLength4[BSCNbTheta*4+i]
00511 -xHalfLength3[BSCNbTheta*4+i])
00512 /yHalfLength2[BSCNbTheta*4+i]/2.;
00513 tanAlpha2[BSCNbTheta*5+i]=(xHalfLength4[BSCNbTheta*5+i]
00514 -xHalfLength3[BSCNbTheta*5+i])
00515 /yHalfLength2[BSCNbTheta*5+i]/2.;
00516
00517 zPosition[BSCNbTheta*5+i]=zPosition[BSCNbTheta*3+i]=zPosition[i]=
00518 zPosition[i]+totalThickness/2.*cos(thetaPosition[i])
00519 -yHalfLength1[BSCNbTheta*2+i]*sin(thetaPosition[i]);
00520 zPosition[i]=totalThickness/2.;
00521 xPosition[BSCNbTheta*5+i]=xPosition[BSCNbTheta*3+i]=xPosition[i]=
00522 xPosition[i]+totalThickness/2.*sin(thetaPosition[i])
00523 +totalThickness*(1./cos(thetaPosition[i+1]-thetaPosition[i])-1)/2.
00524 *cos(thetaPosition[i]);
00525 xPosition[i]=totalThickness*(1.-1./cos(angle2)/sin(angle1))/2.;
00526
00527 yPosition[i]=yPosition[i]
00528 +totalThickness*(1.-1./cos(angle2)/sin(angle1))/2.;
00529 yPosition[i]=yHalfLength1[BSCNbTheta*2+i]-totalThickness/2.;
00530 yPosition[BSCNbTheta*3+i]=yPosition[i]*2.+xHalfLength1[BSCNbTheta*3+i];
00531 yPosition[BSCNbTheta*5+i]=xHalfLength1[BSCNbTheta*5+i];
00532
00533 xPosition[BSCNbTheta+i]=BSCPhiRmin
00534 +zHalfLength[BSCNbTheta+i]*sin(thetaPosition[i])
00535 +(3.*yHalfLength1[BSCNbTheta+i]-yHalfLength2[BSCNbTheta+i])/2.
00536 *cos(thetaPosition[i]);
00537 yPosition[BSCNbTheta+i]=(xHalfLength1[BSCNbTheta+i]
00538 +xHalfLength3[BSCNbTheta+i]
00539 +xHalfLength2[BSCNbTheta+i]
00540 +xHalfLength4[BSCNbTheta+i])/4.;
00541 zPosition[BSCNbTheta+i]=BSCPosition1+rminprojected/tan(thetaPosition[i])
00542 +(2.*yHalfLength1[BSCNbTheta+i]/tan(thetaPosition[i])
00543 +zHalfLength[BSCNbTheta+i])*cos(thetaPosition[i])
00544 +(yHalfLength1[BSCNbTheta+i]+yHalfLength2[BSCNbTheta+i])/2.
00545 *sin(thetaPosition[i]);
00546
00547 xPosition[BSCNbTheta*2+i]=xPosition[i]
00548 +((yHalfLength1[i]+yHalfLength2[i])/2.+yHalfLength1[BSCNbTheta*2+i])
00549 *cos(thetaPosition[i]);
00550 zPosition[BSCNbTheta*2+i]=zPosition[i]
00551 -((yHalfLength1[i]+yHalfLength2[i])/2.+yHalfLength1[BSCNbTheta*2+i])
00552 *sin(thetaPosition[i]);
00553 yPosition[BSCNbTheta*2+i]=(xHalfLength1[BSCNbTheta*2+i]
00554 +xHalfLength3[BSCNbTheta*2+i]
00555 +xHalfLength2[BSCNbTheta*2+i]
00556 +xHalfLength4[BSCNbTheta*2+i])/4.;
00557
00558 xPosition[BSCNbTheta*4+i]=xPosition[i]
00559 -((yHalfLength1[i]+yHalfLength2[i])/2.+yHalfLength1[BSCNbTheta*4+i])
00560 *cos(thetaPosition[i]);
00561 zPosition[BSCNbTheta*4+i]=zPosition[i]
00562 -((yHalfLength1[i]+yHalfLength2[i])/2.+yHalfLength1[BSCNbTheta*4+i])
00563 *sin(thetaPosition[i]);
00564 yPosition[BSCNbTheta*4+i]=(xHalfLength1[BSCNbTheta*4+i]
00565 +xHalfLength3[BSCNbTheta*4+i]
00566 +xHalfLength2[BSCNbTheta*4+i]
00567 +xHalfLength4[BSCNbTheta*4+i])/4.;
00568
00569 }
00570
00571 if(verboseLevel>1)
00572 for(i=0;i<BSCNbTheta*6;i++)
00573 {
00574 G4cout << "The sizes of the " << i+1 << " crystal are:" << G4endl
00575 << "zHalfLength =" << zHalfLength[i]/cm << "(cm)," << G4endl
00576 << "thetaAxis =" << thetaAxis[i]/deg << "(deg),"<< G4endl
00577 << "phiAxis =" << phiAxis[i]/deg << "(deg),"<< G4endl
00578 << "yHalfLength1=" << yHalfLength1[i]/cm << "(cm)," << G4endl
00579 << "xHalfLength1=" << xHalfLength1[i]/cm << "(cm)," << G4endl
00580 << "xHalfLength2=" << xHalfLength2[i]/cm << "(cm)," << G4endl
00581 << "tanAlpha1 =" << tanAlpha1[i] << G4endl
00582 << "yHalfLength2=" << yHalfLength2[i]/cm << "(cm)," << G4endl
00583 << "xHalfLength3=" << xHalfLength3[i]/cm << "(cm)," << G4endl
00584 << "xHalfLength4=" << xHalfLength4[i]/cm << "(cm)," << G4endl
00585 << "tanAlpha2 =" << tanAlpha2[i] << "." << G4endl
00586 << "The position of the " << i+1 << " crystal is:" << G4endl
00587 << "(" << xPosition[i]/cm << ","
00588 << yPosition[i]/cm << ","
00589 << zPosition[i]/cm << ")cm" << G4endl;
00590 }
00591
00592 }
00593
00594 void BesEmcGeometry::SetCasingThickness(G4ThreeVector val)
00595 {
00596
00597 fTyvekThickness = val('X');
00598 fAlThickness = val('Y');
00599 fMylarThickness = val('Z');
00600 }
00601
00602 G4double BesEmcGeometry::GetXPosition(G4int NbCrystal)
00603 {
00604 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
00605 {
00606 return xPosition[NbCrystal];
00607 }
00608 else
00609 {
00610 return 0.;
00611 }
00612 }
00613
00614 G4double BesEmcGeometry::GetYPosition(G4int NbCrystal)
00615 {
00616 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
00617 {
00618 return yPosition[NbCrystal];
00619 }
00620 else
00621 {
00622 return 0.;
00623 }
00624 }
00625
00626 G4double BesEmcGeometry::GetZPosition(G4int NbCrystal)
00627 {
00628 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
00629 {
00630 return zPosition[NbCrystal];
00631 }
00632 else
00633 {
00634 return 0.;
00635 }
00636 }
00637
00638 G4double BesEmcGeometry::GetThetaPosition(G4int NbCrystal)
00639 {
00640 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
00641 {
00642 return thetaPosition[NbCrystal];
00643 }
00644 else
00645 {
00646 return 0.;
00647 }
00648 }
00649
00650 G4double BesEmcGeometry::GetZHalfLength(G4int NbCrystal)
00651 {
00652 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
00653 {
00654 return zHalfLength[NbCrystal];
00655 }
00656 else
00657 {
00658 return 0.;
00659 }
00660 }
00661
00662 G4double BesEmcGeometry::GetThetaAxis(G4int NbCrystal)
00663 {
00664 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
00665 {
00666 return thetaAxis[NbCrystal];
00667 }
00668 else
00669 {
00670 return 0.;
00671 }
00672 }
00673
00674 G4double BesEmcGeometry::GetPhiAxis(G4int NbCrystal)
00675 {
00676 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
00677 {
00678 return phiAxis[NbCrystal];
00679 }
00680 else
00681 {
00682 return 0.;
00683 }
00684 }
00685
00686 G4double BesEmcGeometry::GetYHalfLength1(G4int NbCrystal)
00687 {
00688 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
00689 {
00690 return yHalfLength1[NbCrystal];
00691 }
00692 else
00693 {
00694 return 0.;
00695 }
00696 }
00697
00698 G4double BesEmcGeometry::GetXHalfLength1(G4int NbCrystal)
00699 {
00700 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
00701 {
00702 return xHalfLength1[NbCrystal];
00703 }
00704 else
00705 {
00706 return 0.;
00707 }
00708 }
00709
00710 G4double BesEmcGeometry::GetXHalfLength2(G4int NbCrystal)
00711 {
00712 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
00713 {
00714 return xHalfLength2[NbCrystal];
00715 }
00716 else
00717 {
00718 return 0.;
00719 }
00720 }
00721
00722 G4double BesEmcGeometry::GetTanAlpha1(G4int NbCrystal)
00723 {
00724 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
00725 {
00726 return tanAlpha1[NbCrystal];
00727 }
00728 else
00729 {
00730 return 0.;
00731 }
00732 }
00733
00734 G4double BesEmcGeometry::GetYHalfLength2(G4int NbCrystal)
00735 {
00736 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
00737 {
00738 return yHalfLength2[NbCrystal];
00739 }
00740 else
00741 {
00742 return 0.;
00743 }
00744 }
00745
00746 G4double BesEmcGeometry::GetXHalfLength3(G4int NbCrystal)
00747 {
00748 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
00749 {
00750 return xHalfLength3[NbCrystal];
00751 }
00752 else
00753 {
00754 return 0.;
00755 }
00756 }
00757
00758 G4double BesEmcGeometry::GetXHalfLength4(G4int NbCrystal)
00759 {
00760 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
00761 {
00762 return xHalfLength4[NbCrystal];
00763 }
00764 else
00765 {
00766 return 0.;
00767 }
00768 }
00769
00770 G4double BesEmcGeometry::GetTanAlpha2(G4int NbCrystal)
00771 {
00772 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
00773 {
00774 return tanAlpha2[NbCrystal];
00775 }
00776 else
00777 {
00778 return 0.;
00779 }
00780 }
00781
00782 G4VPhysicalVolume* BesEmcGeometry::GetPhysiBSCCrystal(G4int NbCrystal)
00783 {
00784 if(NbCrystal>=0&&NbCrystal<BSCNbTheta*6)
00785 {
00786 return physiBSCCrystal[NbCrystal];
00787 }
00788 else
00789 {
00790 return NULL;
00791 }
00792 }
