3#include "G4ThreeVector.hh"
4#include "G4PVPlacement.hh"
5#include "G4LogicalVolume.hh"
6#include "G4VPhysicalVolume.hh"
9#include "G4ReflectionFactory.hh"
19#include "G4IrregBox.hh"
21#include "G4Transform3D.hh"
25#include "G4UnionSolid.hh"
26#include "G4SubtractionSolid.hh"
27#include "G4Polyhedra.hh"
30#include "G4Material.hh"
32#include "G4PVParameterised.hh"
33#include "G4PVReplica.hh"
35#include "G4UniformMagField.hh"
36#include "G4FieldManager.hh"
37#include "G4TransportationManager.hh"
38#include "G4SDManager.hh"
39#include "G4RunManager.hh"
40#include "G4VisAttributes.hh"
52 solidEMC(0),logicEMC(0),physiEMC(0),
53 solidBSCPhi(0),logicBSCPhi(0),physiBSCPhi(0),
54 solidBSCTheta(0),logicBSCTheta(0),physiBSCTheta(0),
55 solidBSCCrystal(0),logicBSCCrystal(0),physiBSCCrystal(0),
56 magField(0),detectorMessenger(0),
57 besEMCSD(0),crystalParam(0)
59 if(!fBesEmcConstruction) fBesEmcConstruction=
this;
70 if(crystalParam)
delete crystalParam;
71 if(besEMCGeometry)
delete besEMCGeometry;
72 if(emcEnd)
delete emcEnd;
99 if(logicEMC) physiEMC =
new G4PVPlacement(0,G4ThreeVector(0.0 ,0.0 ,0.0),logicEMC,
"physicalEMC",logicBes,
false, 0);
851void ExtBesEmcConstruction::DefineMaterials()
853 G4String name, symbol;
854 G4double a, z, density;
858 G4int ncomponents, natoms;
859 G4double fractionmass;
868 G4Element*
H=G4Element::GetElement(
"Hydrogen");
872 H =
new G4Element(name=
"Hydrogen",symbol=
"H" , z= 1., a);
874 G4Element*
C=G4Element::GetElement(
"Carbon");
878 C =
new G4Element(name=
"Carbon" ,symbol=
"C" , z= 6., a);
880 G4Element* O=G4Element::GetElement(
"Oxygen");
884 O =
new G4Element(name=
"Oxygen" ,symbol=
"O" , z= 8., a);
887 density = 0.344*g/cm3;
888 G4Material* Tyvek =
new G4Material(name=
"Polyethylene", density, ncomponents=2);
889 Tyvek->AddElement(
C, natoms=1);
890 Tyvek->AddElement(
H, natoms=2);
892 density = 1.39*g/cm3;
893 G4Material* Mylar =
new G4Material(name=
"PolyethyleneTerephthlate", density, ncomponents=3);
894 Mylar->AddElement(
C, natoms=5);
895 Mylar->AddElement(
H, natoms=4);
896 Mylar->AddElement(O, natoms=2);
898 density = 1.18*g/cm3;
899 organicGlass =
new G4Material(name=
"OrganicGlass", density, ncomponents=3);
900 organicGlass->AddElement(
C, natoms=5);
901 organicGlass->AddElement(
H, natoms=7);
902 organicGlass->AddElement(O, natoms=2);
904 G4Material *Fe =
new G4Material(name=
"Iron", z=26., a=55.85*g/mole, density=7.87*g/cm3);
905 G4Material *Cr =
new G4Material(name=
"Chromium", z=24., a=52.00*g/mole, density=8.72*g/cm3);
906 G4Material *Ni =
new G4Material(name=
"Nickel", z=28., a=58.69*g/mole, density=8.72*g/cm3);
908 stainlessSteel =
new G4Material(name=
"0Cr18Ni9", density=7.85*g/cm3, ncomponents=3);
909 stainlessSteel->AddMaterial(Fe, fractionmass=73.*perCent);
910 stainlessSteel->AddMaterial(Cr, fractionmass=18.*perCent);
911 stainlessSteel->AddMaterial(Ni, fractionmass=9.*perCent);
913 G4Material *H2O = G4Material::GetMaterial(
"Water");
914 G4Material *Cu = G4Material::GetMaterial(
"Copper");
915 G4double dWater = 1.*g/cm3;
916 G4double dCopper = 8.96*g/cm3;
917 G4double aWater = ((*besEMCGeometry).waterPipeDr-(*besEMCGeometry).waterPipeThickness)
918 *((*besEMCGeometry).waterPipeDr-(*besEMCGeometry).waterPipeThickness);
919 G4double aCopper = (*besEMCGeometry).waterPipeDr*(*besEMCGeometry).waterPipeDr-aWater;
920 density = (dWater*aWater+dCopper*aCopper)/(aWater+aCopper);
922 waterPipe =
new G4Material(name=
"WaterPipe", density, ncomponents=2);
923 fractionmass = dWater*aWater/(dWater*aWater+dCopper*aCopper);
924 waterPipe->AddMaterial(H2O, fractionmass);
925 fractionmass = dCopper*aCopper/(dWater*aWater+dCopper*aCopper);
926 waterPipe->AddMaterial(Cu, fractionmass);
928 cable =
new G4Material(name=
"Cable", density=4.*g/cm3, ncomponents=1);
929 cable->AddMaterial(Cu,1);
937 G4Material*
Al=G4Material::GetMaterial(
"Aluminium");
940 Al =
new G4Material(name=
"Aluminium", z=13., a=26.98*g/mole, density=2.700*g/cm3);
943 G4Material *Si=G4Material::GetMaterial(
"Silicon");
946 Si =
new G4Material(name=
"Silicon", z=14., a=28.0855*g/mole, density=2.33*g/cm3);
950 G4double totalThickness=(*besEMCGeometry).fTyvekThickness
951 +(*besEMCGeometry).fAlThickness+(*besEMCGeometry).fMylarThickness;
952 density = (Tyvek->GetDensity()*(*besEMCGeometry).fTyvekThickness+
953 Al->GetDensity()*(*besEMCGeometry).fAlThickness+
954 Mylar->GetDensity()*(*besEMCGeometry).fMylarThickness)
956 G4Material* Casing =
new G4Material(name=
"Casing", density, ncomponents=3);
959 fractionmass=Tyvek->GetDensity()/density
960 *(*besEMCGeometry).fTyvekThickness
964 fractionmass=
Al->GetDensity()/density
965 *(*besEMCGeometry).fAlThickness
969 fractionmass=Mylar->GetDensity()/density
970 *(*besEMCGeometry).fMylarThickness
972 fCasingMaterial = Casing;
973 rearCasingMaterial = Tyvek;
976 fCrystalMaterial = G4Material::GetMaterial(
"Cesiumiodide");
982 G4Material* fCrystalMaterial = G4Material::GetMaterial(
"Cesiumiodide");
990 solidEnd =
new G4Cons(
"EndWorld",(*emcEnd).WorldRmin1,(*emcEnd).WorldRmax1,(*emcEnd).WorldRmin2,(*emcEnd).WorldRmax2,
991 (*emcEnd).WorldDz/2,0.*deg,360.*deg);
992 logicEnd =
new G4LogicalVolume(solidEnd, G4Material::GetMaterial(
"Air"),
"EndWorld", 0, 0, 0);
993 physiEnd =
new G4PVPlacement(0,
994 G4ThreeVector(0,0,(*emcEnd).WorldZPosition),
1005 G4RotationMatrix *rotateEnd =
new G4RotationMatrix();
1006 rotateEnd->rotateY(180.*deg);
1007 physiEnd =
new G4PVPlacement(rotateEnd,
1008 G4ThreeVector(0,0,-(*emcEnd).WorldZPosition),
1040 solidEndPhi =
new G4Cons(
"EndPhi",(*emcEnd).WorldRmin1,(*emcEnd).WorldRmax1,(*emcEnd).WorldRmin2,(*emcEnd).WorldRmax2,
1041 (*emcEnd).WorldDz/2,0*deg,22.5*deg);
1042 logicEndPhi =
new G4LogicalVolume(solidEndPhi, G4Material::GetMaterial(
"Air"),
"EndPhi", 0, 0, 0);
1043 for(G4int i=0;i<14;i++)
1047 G4RotationMatrix *rotatePhi =
new G4RotationMatrix();
1048 rotatePhi->rotateZ(-i*22.5*deg+67.5*deg);
1049 physiEndPhi =
new G4PVPlacement(rotatePhi,0,logicEndPhi,
"EndPhi",logicEnd,
false,i);
1055 for(G4int i=0;i<35;i++)
1059 solidEndCasing =
new G4IrregBox(
"EndCasing",(*emcEnd).fPnt[i]);
1060 logicEndCasing =
new G4LogicalVolume(solidEndCasing,fCasingMaterial,
"EndCasing");
1061 physiEndCasing =
new G4PVPlacement(0,0,logicEndCasing,
"EndCasing",logicEndPhi,
false,copyNb);
1064 solidEndCrystal =
new G4IrregBox(
"EndCrystal",(*emcEnd).cryPoint);
1065 logicEndCrystal =
new G4LogicalVolume(solidEndCrystal,fCrystalMaterial,
"EndCrystal");
1066 physiEndCrystal =
new G4PVPlacement(0,0,logicEndCrystal,
"EndCrystal",logicEndCasing,
false,copyNb);
1076 solidEndPhi =
new G4Cons(
"EndPhi",(*emcEnd).WorldRmin1,(*emcEnd).WorldRmax1,(*emcEnd).WorldRmin2,(*emcEnd).WorldRmax2,
1077 (*emcEnd).WorldDz/2,67.5*deg,22.5*deg);
1078 logicEndPhi =
new G4LogicalVolume(solidEndPhi, G4Material::GetMaterial(
"Air"),
"EndPhi", 0, 0, 0);
1079 for(G4int i=0;i<2;i++)
1081 G4RotationMatrix *rotatePhi =
new G4RotationMatrix();
1082 rotatePhi->rotateZ(-i*180.*deg);
1083 physiEndPhi =
new G4PVPlacement(rotatePhi,0,logicEndPhi,
"EndPhi",logicEnd,
false,i*8+6);
1088 for(G4int i=0;i<35;i++)
1091 solidEndCasing =
new G4IrregBox(
"EndCasing",(*emcEnd).fPnt1[i]);
1092 logicEndCasing =
new G4LogicalVolume(solidEndCasing,fCasingMaterial,
"EndCasing");
1093 physiEndCasing =
new G4PVPlacement(0,0,logicEndCasing,
"EndCasing",logicEndPhi,
false,copyNb);
1096 solidEndCrystal =
new G4IrregBox(
"EndCrystal",(*emcEnd).cryPoint);
1097 logicEndCrystal =
new G4LogicalVolume(solidEndCrystal,fCrystalMaterial,
"EndCrystal");
1098 physiEndCrystal =
new G4PVPlacement(0,0,logicEndCrystal,
"EndCrystal",logicEndCasing,
false,copyNb);
1105 (*emcEnd).ReflectX();
1108 for(G4int i=0;i<35;i++)
1109 for (G4int j=0;j<8;j++)
1110 (*emcEnd).fPnt1[i][j].rotateZ(-90.*deg);
1112 solidEndPhi =
new G4Cons(
"EndPhi",(*emcEnd).WorldRmin1,(*emcEnd).WorldRmax1,(*emcEnd).WorldRmin2,(*emcEnd).WorldRmax2,
1113 (*emcEnd).WorldDz/2,0*deg,22.5*deg);
1114 logicEndPhi =
new G4LogicalVolume(solidEndPhi, G4Material::GetMaterial(
"Air"),
"EndPhi", 0, 0, 0);
1115 for(G4int i=0;i<2;i++)
1117 G4RotationMatrix *rotatePhi =
new G4RotationMatrix();
1118 rotatePhi->rotateZ(-i*180.*deg-90.*deg);
1119 physiEndPhi =
new G4PVPlacement(rotatePhi,0,logicEndPhi,
"EndPhi",logicEnd,
false,i*8+7);
1124 for(G4int i=0;i<35;i++)
1127 solidEndCasing =
new G4IrregBox(
"EndCasing",(*emcEnd).fPnt1[i]);
1128 logicEndCasing =
new G4LogicalVolume(solidEndCasing,fCrystalMaterial,
"EndCasing");
1129 physiEndCasing =
new G4PVPlacement(0,0,logicEndCasing,
"EndCasing",logicEndPhi,
false,copyNb);
1132 solidEndCrystal =
new G4IrregBox(
"EndCrystal",(*emcEnd).cryPoint);
1133 logicEndCrystal =
new G4LogicalVolume(solidEndCrystal,fCrystalMaterial,
"EndCrystal");
1134 physiEndCrystal =
new G4PVPlacement(0,0,logicEndCrystal,
"EndCrystal",logicEndCasing,
false,copyNb);
1148 solidSupportBar =
new G4Tubs(
"SupportBar",
1149 (*besEMCGeometry).BSCRmax+(*besEMCGeometry).SPBarThickness1,
1150 (*besEMCGeometry).BSCRmax+(*besEMCGeometry).SPBarThickness+(*besEMCGeometry).SPBarThickness1,
1151 (*besEMCGeometry).BSCDz+(*besEMCGeometry).TaperRingThickness3+(*besEMCGeometry).EndRingDz,
1155 logicSupportBar =
new G4LogicalVolume(solidSupportBar,stainlessSteel,
"SupportBar");
1157 physiSupportBar =
new G4PVPlacement(0,0,logicSupportBar,
"SupportBar",logicEMC,
false,0);
1159 solidSupportBar1 =
new G4Tubs(
"SupportBar1",
1160 (*besEMCGeometry).BSCRmax,
1161 (*besEMCGeometry).BSCRmax+(*besEMCGeometry).SPBarThickness1,
1162 (*besEMCGeometry).BSCDz+(*besEMCGeometry).TaperRingThickness3,
1163 (*besEMCGeometry).BSCPhiDphi-(*besEMCGeometry).SPBarDphi/2,
1164 (*besEMCGeometry).SPBarDphi);
1166 logicSupportBar1 =
new G4LogicalVolume(solidSupportBar1,stainlessSteel,
"SupportBar1");
1168 for(G4int i=0;i<(*besEMCGeometry).BSCNbPhi/2;i++)
1170 G4RotationMatrix *rotateSPBar =
new G4RotationMatrix();
1171 rotateSPBar->rotateZ((*besEMCGeometry).BSCPhiDphi-i*2*(*besEMCGeometry).BSCPhiDphi);
1172 physiSupportBar1 =
new G4PVPlacement(rotateSPBar,0,logicSupportBar1,
"SupportBar1",logicEMC,
false,0);
1176 solidEndRing =
new G4Tubs(
"EndRing",
1177 (*besEMCGeometry).EndRingRmin,
1178 (*besEMCGeometry).EndRingRmin+(*besEMCGeometry).EndRingDr/2,
1179 (*besEMCGeometry).EndRingDz/2,
1183 solidGear =
new G4Tubs(
"Gear",
1184 (*besEMCGeometry).EndRingRmin+(*besEMCGeometry).EndRingDr/2,
1185 (*besEMCGeometry).EndRingRmin+(*besEMCGeometry).EndRingDr,
1186 (*besEMCGeometry).EndRingDz/2,
1188 (*besEMCGeometry).BSCPhiDphi);
1191 solidTaperRing1 =
new G4Tubs(
"TaperRing1",
1192 (*besEMCGeometry).TaperRingRmin1,
1193 (*besEMCGeometry).TaperRingRmin1+(*besEMCGeometry).TaperRingThickness1,
1194 (*besEMCGeometry).TaperRingInnerLength/2,
1198 solidTaperRing2 =
new G4Cons(
"TaperRing2",
1199 (*besEMCGeometry).TaperRingRmin1,
1200 (*besEMCGeometry).TaperRingRmin1+(*besEMCGeometry).TaperRingDr,
1201 (*besEMCGeometry).TaperRingRmin2,
1202 (*besEMCGeometry).TaperRingRmin2+(*besEMCGeometry).TaperRingDr,
1203 (*besEMCGeometry).TaperRingDz/2,
1207 solidTaperRing3 =
new G4Cons(
"TaperRing3",
1208 (*besEMCGeometry).BSCRmax2,
1209 (*besEMCGeometry).BSCRmax2+(*besEMCGeometry).TaperRingOuterLength1,
1210 (*besEMCGeometry).TaperRingRmin2+(*besEMCGeometry).TaperRingDr,
1211 (*besEMCGeometry).TaperRingRmin2+(*besEMCGeometry).TaperRingDr+(*besEMCGeometry).TaperRingOuterLength,
1212 (*besEMCGeometry).TaperRingThickness3/2,
1216 logicEndRing =
new G4LogicalVolume(solidEndRing,stainlessSteel,
"EmcEndRing");
1217 logicGear =
new G4LogicalVolume(solidGear,stainlessSteel,
"Gear");
1218 logicTaperRing1 =
new G4LogicalVolume(solidTaperRing1,stainlessSteel,
"TaperRing1");
1219 logicTaperRing2 =
new G4LogicalVolume(solidTaperRing2,stainlessSteel,
"TaperRing2");
1220 logicTaperRing3 =
new G4LogicalVolume(solidTaperRing3,stainlessSteel,
"TaperRing3");
1222 for(G4int i=0;i<2;i++)
1224 G4RotationMatrix *rotateSPRing =
new G4RotationMatrix();
1225 G4double zEndRing,z1,z2,z3;
1228 zEndRing = (*besEMCGeometry).BSCDz+(*besEMCGeometry).TaperRingThickness3+(*besEMCGeometry).EndRingDz/2;
1229 z1 = (*besEMCGeometry).BSCDz+(*besEMCGeometry).TaperRingThickness3
1230 -(*besEMCGeometry).TaperRingDz-(*besEMCGeometry).TaperRingInnerLength/2;
1231 z2 = (*besEMCGeometry).BSCDz+(*besEMCGeometry).TaperRingThickness3-(*besEMCGeometry).TaperRingDz/2;
1232 z3 = (*besEMCGeometry).BSCDz+(*besEMCGeometry).TaperRingThickness3/2;
1236 rotateSPRing->rotateY(180.*deg);
1237 zEndRing = -((*besEMCGeometry).BSCDz+(*besEMCGeometry).TaperRingThickness3+(*besEMCGeometry).EndRingDz/2);
1238 z1 = -((*besEMCGeometry).BSCDz+(*besEMCGeometry).TaperRingThickness3
1239 -(*besEMCGeometry).TaperRingDz-(*besEMCGeometry).TaperRingInnerLength/2);
1240 z2 = -((*besEMCGeometry).BSCDz+(*besEMCGeometry).TaperRingThickness3-(*besEMCGeometry).TaperRingDz/2);
1241 z3 = -((*besEMCGeometry).BSCDz+(*besEMCGeometry).TaperRingThickness3/2);
1244 physiEndRing =
new G4PVPlacement(rotateSPRing,G4ThreeVector(0,0,zEndRing),
1245 logicEndRing,
"EndRing",logicEMC,
false,0);
1247 for(G4int j=0;j<(*besEMCGeometry).BSCNbPhi/2;j++)
1249 G4RotationMatrix *rotateGear =
new G4RotationMatrix();
1250 rotateGear->rotateZ((*besEMCGeometry).BSCPhiDphi/2-j*2*(*besEMCGeometry).BSCPhiDphi);
1251 physiGear =
new G4PVPlacement(rotateGear,G4ThreeVector(0,0,zEndRing),
1252 logicGear,
"Gear",logicEMC,
false,0);
1255 physiTaperRing1 =
new G4PVPlacement(rotateSPRing,G4ThreeVector(0,0,z1),
1256 logicTaperRing1,
"TaperRing1",logicEMC,
false,0);
1258 physiTaperRing2 =
new G4PVPlacement(rotateSPRing,G4ThreeVector(0,0,z2),
1259 logicTaperRing2,
"TaperRing2",logicEMC,
false,0);
1261 physiTaperRing3 =
new G4PVPlacement(rotateSPRing,G4ThreeVector(0,0,z3),
1262 logicTaperRing3,
"TaperRing3",logicEMC,
false,0);
1304 G4cout <<
"-------------------------------------------------------"<< G4endl
1305 <<
"---> There are "
1306 << phiNbCrystals <<
"(max=" << (*besEMCGeometry).BSCNbPhi
1307 <<
") crystals along phi direction and "
1308 << thetaNbCrystals <<
"(max=" << (*besEMCGeometry).BSCNbTheta
1309 <<
") crystals along theta direction."<< G4endl
1310 <<
"The crystals have sizes of "
1311 << (*besEMCGeometry).BSCCryLength/cm <<
"cm(L) and "
1312 << (*besEMCGeometry).BSCYFront/cm <<
"cm(Y) with "
1313 << fCrystalMaterial->GetName() <<
"."<< G4endl
1314 <<
"The casing is layer of "
1315 << (*besEMCGeometry).fTyvekThickness/mm <<
"mm tyvek,"
1316 << (*besEMCGeometry).fAlThickness/mm <<
"mm aluminum and"
1317 << (*besEMCGeometry).fMylarThickness/mm <<
"mm mylar."<< G4endl
1318 <<
"-------------------------------------------------------"<< G4endl;
1319 G4cout << G4Material::GetMaterial(
"PolyethyleneTerephthlate") << G4endl
1320 << G4Material::GetMaterial(
"Casing") << G4endl
1321 << G4Material::GetMaterial(
"Polyethylene") << G4endl
1322 <<
"-------------------------------------------------------"<< G4endl;
1330 G4Material* pttoMaterial = G4Material::GetMaterial(materialChoice);
1332 {fCrystalMaterial = pttoMaterial;
1333 logicBSCCrystal->SetMaterial(pttoMaterial);
1342 G4Material* pttoMaterial = G4Material::GetMaterial(materialChoice);
1344 {fCasingMaterial = pttoMaterial;
1345 logicBSCTheta->SetMaterial(pttoMaterial);
1355 (*besEMCGeometry).fTyvekThickness = val(
'X');
1356 (*besEMCGeometry).fAlThickness = val(
'Y');
1357 (*besEMCGeometry).fMylarThickness = val(
'Z');
1364 (*besEMCGeometry).BSCRmin = val;
1371 (*besEMCGeometry).BSCNbPhi = val;
1378 (*besEMCGeometry).BSCNbTheta = val;
1391 (*besEMCGeometry).BSCCryLength = val;
1398 (*besEMCGeometry).BSCYFront0 = val;
1406 (*besEMCGeometry).BSCYFront = val;
1413 (*besEMCGeometry).BSCPosition0 = val;
1420 (*besEMCGeometry).BSCPosition1 = val;
1428 G4FieldManager* fieldMgr
1429 = G4TransportationManager::GetTransportationManager()->GetFieldManager();
1431 if(magField)
delete magField;
1434 { magField =
new G4UniformMagField(G4ThreeVector(0.,0.,fieldValue));
1435 fieldMgr->SetDetectorField(magField);
1436 fieldMgr->CreateChordFinder(magField);
1437 fmagField=fieldValue;
1440 fieldMgr->SetDetectorField(magField);
***************************************************************************************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
static EmcG4Geo * Instance()
Get a pointer to the single instance of EmcG4Geo.
void PrintEMCParameters()
void SetBSCRmin(G4double)
void SetBSCPosition0(G4double)
void ConstructEndGeometry(G4LogicalVolume *)
void SetStartIDTheta(G4int)
void SetBSCPosition1(G4double)
void SetCasingThickness(G4ThreeVector)
G4int ComputeEndCopyNb(G4int)
void SetMagField(G4double)
void SetBSCCrystalLength(G4double)
void SetBSCYFront0(G4double)
void SetCrystalMaterial(G4String)
void SetCasingMaterial(G4String)
void Construct(G4LogicalVolume *)
void SetBSCNbTheta(G4int)
void SetBSCYFront(G4double)
void ConstructSPFrame(G4LogicalVolume *, ExtBesEmcGeometry *)
void ModifyForCasing(G4ThreeVector pos[8], G4int CryNb)
G4LogicalVolume * GetTopVolume()
Get the top(world) volume;.