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