File indexing completed on 2025-04-04 08:04:42
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042 #include "DetectorConstruction.hh"
0043
0044 #include "DetectorMessenger.hh"
0045 #include "TargetSD.hh"
0046 #include "TestParameters.hh"
0047
0048 #include "G4Colour.hh"
0049 #include "G4GeometryManager.hh"
0050 #include "G4LogicalVolume.hh"
0051 #include "G4LogicalVolumeStore.hh"
0052 #include "G4Material.hh"
0053 #include "G4NistManager.hh"
0054 #include "G4PVPlacement.hh"
0055 #include "G4PhysicalConstants.hh"
0056 #include "G4PhysicalVolumeStore.hh"
0057 #include "G4PionPlus.hh"
0058 #include "G4ProductionCuts.hh"
0059 #include "G4Region.hh"
0060 #include "G4RegionStore.hh"
0061 #include "G4RunManager.hh"
0062 #include "G4SDManager.hh"
0063 #include "G4SolidStore.hh"
0064 #include "G4SystemOfUnits.hh"
0065 #include "G4Tubs.hh"
0066 #include "G4UnitsTable.hh"
0067 #include "G4VisAttributes.hh"
0068 #include "G4ios.hh"
0069
0070 #include <vector>
0071
0072
0073
0074 DetectorConstruction::DetectorConstruction()
0075 {
0076 fGasThickness = 23.0 * mm;
0077 fGasRadius = 10. * cm;
0078 fMaxStep = DBL_MAX;
0079
0080 fWindowThick = 51.0 * micrometer;
0081
0082 DefineMaterials();
0083
0084 fDetectorMessenger = new DetectorMessenger(this);
0085
0086 G4double cut = 0.7 * mm;
0087 fGasDetectorCuts = new G4ProductionCuts();
0088 fGasDetectorCuts->SetProductionCut(cut, "gamma");
0089 fGasDetectorCuts->SetProductionCut(cut, "e-");
0090 fGasDetectorCuts->SetProductionCut(cut, "e+");
0091 fGasDetectorCuts->SetProductionCut(cut, "proton");
0092 }
0093
0094
0095
0096 DetectorConstruction::~DetectorConstruction()
0097 {
0098 delete fDetectorMessenger;
0099 }
0100
0101
0102
0103 void DetectorConstruction::DefineMaterials()
0104 {
0105
0106 G4String name, symbol;
0107 G4double density;
0108 G4int nel;
0109 G4int ncomponents;
0110 G4double fractionmass;
0111
0112 G4NistManager* manager = G4NistManager::Instance();
0113
0114
0115
0116 G4Element* elH = manager->FindOrBuildElement(1);
0117 G4Element* elC = manager->FindOrBuildElement(6);
0118 G4Element* elO = manager->FindOrBuildElement(8);
0119 G4Element* elF = manager->FindOrBuildElement(9);
0120 G4Element* elNe = manager->FindOrBuildElement(10);
0121 G4Element* elXe = manager->FindOrBuildElement(54);
0122
0123
0124
0125 G4Material* Argon = manager->FindOrBuildMaterial("G4_Ar");
0126 G4Material* Kr = manager->FindOrBuildMaterial("G4_Kr");
0127 G4Material* Xe = manager->FindOrBuildMaterial("G4_Xe");
0128
0129
0130
0131 G4Material* CarbonDioxide = manager->FindOrBuildMaterial("G4_CARBON_DIOXIDE");
0132 G4Material* Mylar = manager->FindOrBuildMaterial("G4_MYLAR");
0133 G4Material* Methane = manager->FindOrBuildMaterial("G4_METHANE");
0134 G4Material* Propane = manager->FindOrBuildMaterial("G4_PROPANE");
0135
0136
0137 manager->ConstructNewGasMaterial("Propane10", "G4_PROPANE", NTP_Temperature, 10. * atmosphere);
0138
0139 G4Material* empty = manager->FindOrBuildMaterial("G4_Galactic");
0140
0141
0142 density = 3.491 * mg / cm3;
0143 G4Material* Kr7CH4 = new G4Material(name = "Kr7CH4", density, ncomponents = 2);
0144 Kr7CH4->AddMaterial(Kr, fractionmass = 0.986);
0145 Kr7CH4->AddMaterial(Methane, fractionmass = 0.014);
0146
0147 G4double TRT_Xe_density = 5.485 * mg / cm3;
0148 G4Material* TRT_Xe = new G4Material(name = "TRT_Xe", TRT_Xe_density, nel = 1, kStateGas,
0149 293.15 * kelvin, 1. * atmosphere);
0150 TRT_Xe->AddElement(elXe, 1);
0151
0152 G4double TRT_CO2_density = 1.842 * mg / cm3;
0153 G4Material* TRT_CO2 = new G4Material(name = "TRT_CO2", TRT_CO2_density, nel = 2, kStateGas,
0154 293.15 * kelvin, 1. * atmosphere);
0155 TRT_CO2->AddElement(elC, 1);
0156 TRT_CO2->AddElement(elO, 2);
0157
0158
0159 std::vector<G4String> trtatom = {"C", "O"};
0160 std::vector<G4int> trtnum = {1, 2};
0161 manager->ConstructNewMaterial("TRT_CO2p", trtatom, trtnum, TRT_CO2_density, true, kStateGas,
0162 NTP_Temperature, atmosphere);
0163
0164 G4double TRT_CF4_density = 3.9 * mg / cm3;
0165 G4Material* TRT_CF4 = new G4Material(name = "TRT_CF4", TRT_CF4_density, nel = 2, kStateGas,
0166 293.15 * kelvin, 1. * atmosphere);
0167 TRT_CF4->AddElement(elC, 1);
0168 TRT_CF4->AddElement(elF, 4);
0169
0170
0171 G4double XeCO2CF4_density = 4.76 * mg / cm3;
0172 G4Material* XeCO2CF4 = new G4Material(name = "XeCO2CF4", XeCO2CF4_density, ncomponents = 3,
0173 kStateGas, 293.15 * kelvin, 1. * atmosphere);
0174 XeCO2CF4->AddMaterial(TRT_Xe, 0.807);
0175 XeCO2CF4->AddMaterial(TRT_CO2, 0.039);
0176 XeCO2CF4->AddMaterial(TRT_CF4, 0.154);
0177
0178
0179 density = 3.758 * mg / cm3;
0180 G4Material* C3H8 =
0181 new G4Material(name = "C3H8", density, nel = 2, kStateGas, 293.15 * kelvin, 2. * atmosphere);
0182 C3H8->AddElement(elC, 3);
0183 C3H8->AddElement(elH, 8);
0184
0185
0186 std::vector<G4String> elmname = {"C", "H"};
0187 std::vector<G4int> atomnum = {3, 8};
0188 manager->ConstructNewIdealGasMaterial("C3H8p", elmname, atomnum, true, 293.15 * kelvin,
0189 2. * atmosphere);
0190
0191
0192 density = 4.9196 * mg / cm3;
0193 G4Material* XeCH4C3H8 = new G4Material(name = "XeCH4C3H8", density, ncomponents = 3, kStateGas,
0194 NTP_Temperature, 1. * atmosphere);
0195 XeCH4C3H8->AddMaterial(Xe, fractionmass = 0.971);
0196 XeCH4C3H8->AddMaterial(Methane, fractionmass = 0.010);
0197 XeCH4C3H8->AddMaterial(Propane, fractionmass = 0.019);
0198
0199
0200 density = 1.709 * mg / cm3;
0201 G4Material* Ar7CH4 = new G4Material(name = "Ar7CH4", density, ncomponents = 2, kStateGas,
0202 STP_Temperature, STP_Pressure);
0203 Ar7CH4->AddMaterial(Argon, fractionmass = 0.971);
0204 Ar7CH4->AddMaterial(Methane, fractionmass = 0.029);
0205
0206
0207 density = 1.8223 * mg / cm3;
0208 G4Material* Ar_80CO2_20 = new G4Material(name = "ArCO2", density, ncomponents = 2, kStateGas,
0209 STP_Temperature, STP_Pressure);
0210 Ar_80CO2_20->AddMaterial(Argon, fractionmass = 0.783);
0211 Ar_80CO2_20->AddMaterial(CarbonDioxide, fractionmass = 0.217);
0212
0213
0214 density = 5.0818 * mg / cm3;
0215 G4Material* Xe20CO2 = new G4Material(name = "Xe20CO2", density, ncomponents = 2, kStateGas,
0216 STP_Temperature, STP_Pressure);
0217 Xe20CO2->AddMaterial(Xe, fractionmass = 0.922);
0218 Xe20CO2->AddMaterial(CarbonDioxide, fractionmass = 0.078);
0219
0220
0221 density = 3.601 * mg / cm3;
0222 G4Material* Kr20CO2 = new G4Material(name = "Kr20CO2", density, ncomponents = 2, kStateGas,
0223 STP_Temperature, STP_Pressure);
0224 Kr20CO2->AddMaterial(Kr, fractionmass = 0.89);
0225 Kr20CO2->AddMaterial(CarbonDioxide, fractionmass = 0.11);
0226
0227
0228 density = 0.939 * mg / cm3;
0229 G4Material* NeCO2 = new G4Material(name = "TPC_Ne-CO2-2", density, ncomponents = 3, kStateGas,
0230 NTP_Temperature, 1. * atmosphere);
0231 NeCO2->AddElement(elNe, fractionmass = 0.8039);
0232 NeCO2->AddElement(elO, fractionmass = 0.1426);
0233 NeCO2->AddElement(elC, fractionmass = 0.0535);
0234
0235
0236 std::vector<G4String> neatom = {"Ne", "O", "C"};
0237 std::vector<G4double> nefr = {0.8039, 0.1426, 0.0536};
0238 manager->ConstructNewMaterial("TPC_Ne-CO2-2p", neatom, nefr, density, true, kStateGas,
0239 NTP_Temperature, atmosphere);
0240
0241
0242 density = 4.9389 * mg / cm3;
0243 G4Material* Xe15CO2 = new G4Material(name = "Xe15CO2", density, ncomponents = 2, kStateGas,
0244 NTP_Temperature, atmosphere);
0245 Xe15CO2->AddMaterial(Xe, fractionmass = 0.944);
0246 Xe15CO2->AddMaterial(CarbonDioxide, fractionmass = 0.056);
0247
0248 fGasMat = XeCH4C3H8;
0249 fWindowMat = Mylar;
0250 fWorldMaterial = empty;
0251
0252 G4cout << *(G4Material::GetMaterialTable()) << G4endl;
0253 }
0254
0255
0256
0257 G4VPhysicalVolume* DetectorConstruction::Construct()
0258 {
0259 if (nullptr != fPhysWorld) {
0260 return fPhysWorld;
0261 }
0262
0263 G4double contThick = fWindowThick * 2 + fGasThickness;
0264 G4double contR = fWindowThick * 2 + fGasRadius;
0265
0266 G4double worldSizeZ = contThick * 1.2;
0267 G4double worldSizeR = contR * 1.2;
0268
0269 TestParameters::GetPointer()->SetPositionZ(-0.55 * contThick);
0270
0271
0272 G4cout << "\n The WORLD is made of " << worldSizeZ / mm << "mm of "
0273 << fWorldMaterial->GetName();
0274 G4cout << ", the transverse size (R) of the world is " << worldSizeR / mm << " mm. " << G4endl;
0275 G4cout << " The CONTAINER is made of " << fWindowThick / mm << "mm of " << fWindowMat->GetName()
0276 << G4endl;
0277 G4cout << " The TARGET is made of " << fGasThickness / mm << "mm of " << fGasMat->GetName();
0278 G4cout << ", the transverse size (R) is " << fGasRadius / mm << " mm. " << G4endl;
0279 G4cout << G4endl;
0280
0281
0282 fSolidWorld = new G4Tubs("World", 0., worldSizeR, worldSizeZ / 2., 0., CLHEP::twopi);
0283
0284 fLogicWorld = new G4LogicalVolume(fSolidWorld, fWorldMaterial, "World");
0285
0286 fPhysWorld = new G4PVPlacement(0, G4ThreeVector(0., 0., 0.), "World", fLogicWorld, 0, false, 0);
0287
0288
0289 fSolidContainer = new G4Tubs("Absorber", 0., contR, contThick / 2., 0., CLHEP::twopi);
0290
0291 fLogicContainer = new G4LogicalVolume(fSolidContainer, fWindowMat, "Window");
0292
0293 G4PVPlacement* PhysWind = new G4PVPlacement(0, G4ThreeVector(0., 0., 0.), "Window",
0294 fLogicContainer, fPhysWorld, false, 0);
0295
0296
0297 fSolidDetector = new G4Tubs("Gas", 0., fGasRadius, fGasThickness / 2., 0., CLHEP::twopi);
0298
0299 fLogicDetector = new G4LogicalVolume(fSolidDetector, fGasMat, "Gas");
0300
0301 new G4PVPlacement(0, G4ThreeVector(0., 0., 0.), "Gas", fLogicDetector, PhysWind, false, 0);
0302
0303
0304 fRegGasDet = new G4Region("GasDetector");
0305 fRegGasDet->SetProductionCuts(fGasDetectorCuts);
0306 fRegGasDet->AddRootLogicalVolume(fLogicDetector);
0307
0308
0309 fLogicWorld->SetVisAttributes(G4VisAttributes::GetInvisible());
0310 G4VisAttributes* color1 = new G4VisAttributes(G4Colour(0.3, 0.3, 0.3));
0311 fLogicContainer->SetVisAttributes(color1);
0312 G4VisAttributes* color2 = new G4VisAttributes(G4Colour(0.0, 0.3, 0.7));
0313 fLogicDetector->SetVisAttributes(color2);
0314
0315 if (0.0 == fGasMat->GetIonisation()->GetMeanEnergyPerIonPair()) {
0316 SetPairEnergy(20 * eV);
0317 }
0318 return fPhysWorld;
0319 }
0320
0321
0322
0323 void DetectorConstruction::ConstructSDandField()
0324 {
0325 auto sd = new TargetSD("GasSD");
0326 G4SDManager::GetSDMpointer()->AddNewDetector(sd);
0327 SetSensitiveDetector(fLogicDetector, sd);
0328 }
0329
0330
0331
0332 void DetectorConstruction::SetGasMaterial(const G4String& name)
0333 {
0334
0335 G4Material* mat = G4Material::GetMaterial(name, false);
0336
0337
0338 if (!mat) {
0339 mat = G4NistManager::Instance()->FindOrBuildMaterial(name);
0340 }
0341
0342 if (mat && mat != fGasMat) {
0343 G4cout << "### New target material: " << mat->GetName() << G4endl;
0344 fGasMat = mat;
0345 if (fLogicDetector) {
0346 fLogicDetector->SetMaterial(mat);
0347 G4RunManager::GetRunManager()->PhysicsHasBeenModified();
0348 }
0349 }
0350 }
0351
0352
0353
0354 void DetectorConstruction::SetContainerMaterial(const G4String& name)
0355 {
0356
0357 G4Material* mat = G4Material::GetMaterial(name, false);
0358
0359
0360 if (!mat) {
0361 mat = G4NistManager::Instance()->FindOrBuildMaterial(name);
0362 }
0363
0364 if (mat && mat != fWindowMat) {
0365 G4cout << "### New material for container: " << mat->GetName() << G4endl;
0366 fWindowMat = mat;
0367 if (fLogicContainer) {
0368 fLogicContainer->SetMaterial(mat);
0369 G4RunManager::GetRunManager()->PhysicsHasBeenModified();
0370 }
0371 }
0372 }
0373
0374
0375
0376 void DetectorConstruction::SetWorldMaterial(const G4String& name)
0377 {
0378
0379 G4Material* mat = G4Material::GetMaterial(name, false);
0380
0381
0382 if (!mat) {
0383 mat = G4NistManager::Instance()->FindOrBuildMaterial(name);
0384 }
0385
0386 if (mat && mat != fWorldMaterial) {
0387 G4cout << "### New World material: " << mat->GetName() << G4endl;
0388 fWorldMaterial = mat;
0389 if (fLogicWorld) {
0390 fLogicWorld->SetMaterial(mat);
0391 G4RunManager::GetRunManager()->PhysicsHasBeenModified();
0392 }
0393 }
0394 }
0395
0396
0397
0398 void DetectorConstruction::SetGasThickness(G4double val)
0399 {
0400 if (val <= 0.0) {
0401 return;
0402 }
0403 fGasThickness = val;
0404 if (nullptr != fPhysWorld) {
0405 ChangeGeometry();
0406 }
0407 }
0408
0409
0410
0411 void DetectorConstruction::SetGasRadius(G4double val)
0412 {
0413 if (val <= 0.0) {
0414 return;
0415 }
0416 fGasRadius = val;
0417 if (nullptr != fPhysWorld) {
0418 ChangeGeometry();
0419 }
0420 }
0421
0422
0423
0424 void DetectorConstruction::SetContainerThickness(G4double val)
0425 {
0426 if (val <= 0.0) {
0427 return;
0428 }
0429 fWindowThick = val;
0430 if (nullptr != fPhysWorld) {
0431 ChangeGeometry();
0432 }
0433 }
0434
0435
0436
0437 void DetectorConstruction::SetPairEnergy(G4double val)
0438 {
0439 if (val > 0.0) {
0440 fGasMat->GetIonisation()->SetMeanEnergyPerIonPair(val);
0441 }
0442 }
0443
0444
0445
0446 void DetectorConstruction::ChangeGeometry()
0447 {
0448 G4double contThick = fWindowThick * 2 + fGasThickness;
0449 G4double contR = fWindowThick * 2 + fGasRadius;
0450
0451 G4double worldSizeZ = contThick * 1.2;
0452 G4double worldSizeR = contR * 1.2;
0453
0454 TestParameters::GetPointer()->SetPositionZ(-0.55 * contThick);
0455
0456 fSolidWorld->SetOuterRadius(worldSizeR);
0457 fSolidWorld->SetZHalfLength(worldSizeZ * 0.5);
0458
0459 fSolidContainer->SetOuterRadius(contR);
0460 fSolidContainer->SetZHalfLength(contThick * 0.5);
0461
0462 fSolidDetector->SetOuterRadius(fGasRadius);
0463 fSolidDetector->SetZHalfLength(fGasThickness * 0.5);
0464 }
0465
0466
