Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-08-13 08:31:08

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file LXeDetectorConstruction.cc
0027 /// \brief Implementation of the LXeDetectorConstruction class
0028 
0029 #include "LXeDetectorConstruction.hh"
0030 
0031 #include "LXeDetectorMessenger.hh"
0032 #include "LXeMainVolume.hh"
0033 #include "LXePMTSD.hh"
0034 #include "LXeScintSD.hh"
0035 #include "LXeWLSSlab.hh"
0036 
0037 #include "G4Box.hh"
0038 #include "G4GeometryManager.hh"
0039 #include "G4LogicalBorderSurface.hh"
0040 #include "G4LogicalSkinSurface.hh"
0041 #include "G4LogicalVolume.hh"
0042 #include "G4LogicalVolumeStore.hh"
0043 #include "G4Material.hh"
0044 #include "G4MaterialTable.hh"
0045 #include "G4OpticalSurface.hh"
0046 #include "G4PVPlacement.hh"
0047 #include "G4PhysicalConstants.hh"
0048 #include "G4PhysicalVolumeStore.hh"
0049 #include "G4RunManager.hh"
0050 #include "G4SDManager.hh"
0051 #include "G4SolidStore.hh"
0052 #include "G4Sphere.hh"
0053 #include "G4SystemOfUnits.hh"
0054 #include "G4ThreeVector.hh"
0055 #include "G4Tubs.hh"
0056 #include "G4UImanager.hh"
0057 #include "G4VisAttributes.hh"
0058 #include "globals.hh"
0059 
0060 G4bool LXeDetectorConstruction::fSphereOn = true;
0061 
0062 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0063 
0064 LXeDetectorConstruction::LXeDetectorConstruction()
0065 {
0066   SetDefaults();
0067   DefineMaterials();
0068   fDetectorMessenger = new LXeDetectorMessenger(this);
0069 }
0070 
0071 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0072 
0073 LXeDetectorConstruction::~LXeDetectorConstruction()
0074 {
0075   delete fMainVolume;
0076   delete fLXe_mt;
0077   delete fDetectorMessenger;
0078   delete fMPTPStyrene;
0079 }
0080 
0081 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0082 
0083 void LXeDetectorConstruction::DefineMaterials()
0084 {
0085   G4double a;  // atomic mass
0086   G4double z;  // atomic number
0087   G4double density;
0088 
0089   G4int polyPMMA = 1;
0090   G4int nC_PMMA = 3 + 2 * polyPMMA;
0091   G4int nH_PMMA = 6 + 2 * polyPMMA;
0092 
0093   G4int polyeth = 1;
0094   G4int nC_eth = 2 * polyeth;
0095   G4int nH_eth = 4 * polyeth;
0096 
0097   //***Elements
0098   fH = new G4Element("H", "H", z = 1., a = 1.01 * g / mole);
0099   fC = new G4Element("C", "C", z = 6., a = 12.01 * g / mole);
0100   fN = new G4Element("N", "N", z = 7., a = 14.01 * g / mole);
0101   fO = new G4Element("O", "O", z = 8., a = 16.00 * g / mole);
0102 
0103   //***Materials
0104   // Liquid Xenon
0105   fLXe = new G4Material("LXe", z = 54., a = 131.29 * g / mole, density = 3.020 * g / cm3);
0106   // Aluminum
0107   fAl = new G4Material("Al", z = 13., a = 26.98 * g / mole, density = 2.7 * g / cm3);
0108   // Vacuum
0109   fVacuum = new G4Material("Vacuum", z = 1., a = 1.01 * g / mole, density = universe_mean_density,
0110                            kStateGas, 0.1 * kelvin, 1.e-19 * pascal);
0111   // Air
0112   fAir = new G4Material("Air", density = 1.29 * mg / cm3, 2);
0113   fAir->AddElement(fN, 70 * perCent);
0114   fAir->AddElement(fO, 30 * perCent);
0115   // Glass
0116   fGlass = new G4Material("Glass", density = 1.032 * g / cm3, 2);
0117   fGlass->AddElement(fC, 91.533 * perCent);
0118   fGlass->AddElement(fH, 8.467 * perCent);
0119   // Polystyrene
0120   fPstyrene = new G4Material("Polystyrene", density = 1.03 * g / cm3, 2);
0121   fPstyrene->AddElement(fC, 8);
0122   fPstyrene->AddElement(fH, 8);
0123   // Fiber(PMMA)
0124   fPMMA = new G4Material("PMMA", density = 1190. * kg / m3, 3);
0125   fPMMA->AddElement(fH, nH_PMMA);
0126   fPMMA->AddElement(fC, nC_PMMA);
0127   fPMMA->AddElement(fO, 2);
0128   // Cladding(polyethylene)
0129   fPethylene1 = new G4Material("Pethylene1", density = 1200. * kg / m3, 2);
0130   fPethylene1->AddElement(fH, nH_eth);
0131   fPethylene1->AddElement(fC, nC_eth);
0132   // Double cladding(flourinated polyethylene)
0133   fPethylene2 = new G4Material("Pethylene2", density = 1400. * kg / m3, 2);
0134   fPethylene2->AddElement(fH, nH_eth);
0135   fPethylene2->AddElement(fC, nC_eth);
0136 
0137   //***Material properties tables
0138 
0139   std::vector<G4double> lxe_Energy = {7.0 * eV, 7.07 * eV, 7.14 * eV};
0140 
0141   std::vector<G4double> lxe_SCINT = {0.1, 1.0, 0.1};
0142   std::vector<G4double> lxe_RIND = {1.59, 1.57, 1.54};
0143   std::vector<G4double> lxe_ABSL = {35. * cm, 35. * cm, 35. * cm};
0144   fLXe_mt = new G4MaterialPropertiesTable();
0145   fLXe_mt->AddProperty("SCINTILLATIONCOMPONENT1", lxe_Energy, lxe_SCINT);
0146   fLXe_mt->AddProperty("SCINTILLATIONCOMPONENT2", lxe_Energy, lxe_SCINT);
0147   fLXe_mt->AddProperty("RINDEX", lxe_Energy, lxe_RIND);
0148   fLXe_mt->AddProperty("ABSLENGTH", lxe_Energy, lxe_ABSL);
0149   fLXe_mt->AddConstProperty("SCINTILLATIONYIELD", 12000. / MeV);
0150   fLXe_mt->AddConstProperty("RESOLUTIONSCALE", 1.0);
0151   fLXe_mt->AddConstProperty("SCINTILLATIONTIMECONSTANT1", 20. * ns);
0152   fLXe_mt->AddConstProperty("SCINTILLATIONTIMECONSTANT2", 45. * ns);
0153   fLXe_mt->AddConstProperty("SCINTILLATIONYIELD1", 1.0);
0154   fLXe_mt->AddConstProperty("SCINTILLATIONYIELD2", 0.0);
0155   fLXe->SetMaterialPropertiesTable(fLXe_mt);
0156 
0157   // Set the Birks Constant for the LXe scintillator
0158   fLXe->GetIonisation()->SetBirksConstant(0.126 * mm / MeV);
0159 
0160   std::vector<G4double> glass_AbsLength = {420. * cm, 420. * cm, 420. * cm};
0161   auto glass_mt = new G4MaterialPropertiesTable();
0162   glass_mt->AddProperty("ABSLENGTH", lxe_Energy, glass_AbsLength);
0163   glass_mt->AddProperty("RINDEX", "Fused Silica");
0164   fGlass->SetMaterialPropertiesTable(glass_mt);
0165 
0166   auto vacuum_mt = new G4MaterialPropertiesTable();
0167   vacuum_mt->AddProperty("RINDEX", "Air");
0168   fVacuum->SetMaterialPropertiesTable(vacuum_mt);
0169   fAir->SetMaterialPropertiesTable(vacuum_mt);  // Give air the same rindex
0170 
0171   std::vector<G4double> wls_Energy = {2.00 * eV, 2.87 * eV, 2.90 * eV, 3.47 * eV};
0172 
0173   std::vector<G4double> rIndexPstyrene = {1.5, 1.5, 1.5, 1.5};
0174   std::vector<G4double> absorption1 = {2. * cm, 2. * cm, 2. * cm, 2. * cm};
0175   std::vector<G4double> scintilFast = {0.0, 0.0, 1.0, 1.0};
0176   fMPTPStyrene = new G4MaterialPropertiesTable();
0177   fMPTPStyrene->AddProperty("RINDEX", wls_Energy, rIndexPstyrene);
0178   fMPTPStyrene->AddProperty("ABSLENGTH", wls_Energy, absorption1);
0179   fMPTPStyrene->AddProperty("SCINTILLATIONCOMPONENT1", wls_Energy, scintilFast);
0180   fMPTPStyrene->AddConstProperty("SCINTILLATIONYIELD", 10. / keV);
0181   fMPTPStyrene->AddConstProperty("RESOLUTIONSCALE", 1.0);
0182   fMPTPStyrene->AddConstProperty("SCINTILLATIONTIMECONSTANT1", 10. * ns);
0183   fPstyrene->SetMaterialPropertiesTable(fMPTPStyrene);
0184 
0185   // Set the Birks Constant for the Polystyrene scintillator
0186   fPstyrene->GetIonisation()->SetBirksConstant(0.126 * mm / MeV);
0187 
0188   std::vector<G4double> AbsFiber = {9.0 * m, 9.0 * m, 0.1 * mm, 0.1 * mm};
0189   std::vector<G4double> EmissionFib = {1.0, 1.0, 0.0, 0.0};
0190   auto fiberProperty = new G4MaterialPropertiesTable();
0191   fiberProperty->AddProperty("RINDEX", "PMMA");
0192   fiberProperty->AddProperty("WLSABSLENGTH", wls_Energy, AbsFiber);
0193   fiberProperty->AddProperty("WLSCOMPONENT", wls_Energy, EmissionFib);
0194   fiberProperty->AddConstProperty("WLSTIMECONSTANT", 0.5 * ns);
0195   fPMMA->SetMaterialPropertiesTable(fiberProperty);
0196 
0197   std::vector<G4double> RefractiveIndexClad1 = {1.49, 1.49, 1.49, 1.49};
0198   auto clad1Property = new G4MaterialPropertiesTable();
0199   clad1Property->AddProperty("RINDEX", wls_Energy, RefractiveIndexClad1);
0200   clad1Property->AddProperty("ABSLENGTH", wls_Energy, AbsFiber);
0201   fPethylene1->SetMaterialPropertiesTable(clad1Property);
0202 
0203   std::vector<G4double> RefractiveIndexClad2 = {1.42, 1.42, 1.42, 1.42};
0204   auto clad2Property = new G4MaterialPropertiesTable();
0205   clad2Property->AddProperty("RINDEX", wls_Energy, RefractiveIndexClad2);
0206   clad2Property->AddProperty("ABSLENGTH", wls_Energy, AbsFiber);
0207   fPethylene2->SetMaterialPropertiesTable(clad2Property);
0208 }
0209 
0210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0211 
0212 G4VPhysicalVolume* LXeDetectorConstruction::Construct()
0213 {
0214   // The experimental hall walls are all 1m away from housing walls
0215   G4double expHall_x = fScint_x + fD_mtl + 1. * m;
0216   G4double expHall_y = fScint_y + fD_mtl + 1. * m;
0217   G4double expHall_z = fScint_z + fD_mtl + 1. * m;
0218 
0219   // Create experimental hall
0220   fExperimentalHall_box = new G4Box("expHall_box", expHall_x, expHall_y, expHall_z);
0221   fExperimentalHall_log = new G4LogicalVolume(fExperimentalHall_box, fVacuum, "expHall_log");
0222   fExperimentalHall_phys = new G4PVPlacement(nullptr, G4ThreeVector(), fExperimentalHall_log,
0223                                              "expHall", nullptr, false, 0);
0224 
0225   fExperimentalHall_log->SetVisAttributes(G4VisAttributes::GetInvisible());
0226 
0227   // Place the main volume
0228   if (fMainVolumeOn) {
0229     fMainVolume =
0230       new LXeMainVolume(nullptr, G4ThreeVector(), fExperimentalHall_log, false, 0, this);
0231   }
0232 
0233   // Place the WLS slab
0234   if (fWLSslab) {
0235     G4VPhysicalVolume* slab =
0236       new LXeWLSSlab(nullptr, G4ThreeVector(0., 0., -fScint_z / 2. - fSlab_z - 1. * cm),
0237                      fExperimentalHall_log, false, 0, this);
0238 
0239     // Surface properties for the WLS slab
0240     auto scintWrap = new G4OpticalSurface("ScintWrap");
0241 
0242     new G4LogicalBorderSurface("ScintWrap", slab, fExperimentalHall_phys, scintWrap);
0243 
0244     scintWrap->SetType(dielectric_metal);
0245     scintWrap->SetFinish(polished);
0246     scintWrap->SetModel(glisur);
0247 
0248     std::vector<G4double> pp = {2.0 * eV, 3.5 * eV};
0249     std::vector<G4double> reflectivity = {1.0, 1.0};
0250     std::vector<G4double> efficiency = {0.0, 0.0};
0251 
0252     auto scintWrapProperty = new G4MaterialPropertiesTable();
0253 
0254     scintWrapProperty->AddProperty("REFLECTIVITY", pp, reflectivity);
0255     scintWrapProperty->AddProperty("EFFICIENCY", pp, efficiency);
0256     scintWrap->SetMaterialPropertiesTable(scintWrapProperty);
0257   }
0258 
0259   return fExperimentalHall_phys;
0260 }
0261 
0262 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0263 
0264 void LXeDetectorConstruction::ConstructSDandField()
0265 {
0266   if (!fMainVolume) return;
0267 
0268   // PMT SD
0269 
0270   LXePMTSD* pmt = fPmt_SD.Get();
0271   if (!pmt) {
0272     // Created here so it exists as pmts are being placed
0273     G4cout << "Construction /LXeDet/pmtSD" << G4endl;
0274     auto pmt_SD = new LXePMTSD("/LXeDet/pmtSD");
0275     fPmt_SD.Put(pmt_SD);
0276 
0277     pmt_SD->InitPMTs();
0278     pmt_SD->SetPmtPositions(fMainVolume->GetPmtPositions());
0279   }
0280   else {
0281     pmt->InitPMTs();
0282     pmt->SetPmtPositions(fMainVolume->GetPmtPositions());
0283   }
0284   G4SDManager::GetSDMpointer()->AddNewDetector(fPmt_SD.Get());
0285   // sensitive detector is not actually on the photocathode.
0286   // processHits gets done manually by the stepping action.
0287   // It is used to detect when photons hit and get absorbed & detected at the
0288   // boundary to the photocathode (which doesn't get done by attaching it to a
0289   // logical volume.
0290   // It does however need to be attached to something or else it doesn't get
0291   // reset at the begining of events
0292 
0293   SetSensitiveDetector(fMainVolume->GetLogPhotoCath(), fPmt_SD.Get());
0294 
0295   // Scint SD
0296 
0297   if (!fScint_SD.Get()) {
0298     G4cout << "Construction /LXeDet/scintSD" << G4endl;
0299     auto scint_SD = new LXeScintSD("/LXeDet/scintSD");
0300     fScint_SD.Put(scint_SD);
0301   }
0302   G4SDManager::GetSDMpointer()->AddNewDetector(fScint_SD.Get());
0303   SetSensitiveDetector(fMainVolume->GetLogScint(), fScint_SD.Get());
0304 }
0305 
0306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0307 
0308 void LXeDetectorConstruction::SetDimensions(G4ThreeVector dims)
0309 {
0310   fScint_x = dims[0];
0311   fScint_y = dims[1];
0312   fScint_z = dims[2];
0313   G4RunManager::GetRunManager()->ReinitializeGeometry();
0314 }
0315 
0316 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0317 
0318 void LXeDetectorConstruction::SetHousingThickness(G4double d_mtl)
0319 {
0320   fD_mtl = d_mtl;
0321   G4RunManager::GetRunManager()->ReinitializeGeometry();
0322 }
0323 
0324 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0325 
0326 void LXeDetectorConstruction::SetNX(G4int nx)
0327 {
0328   fNx = nx;
0329   G4RunManager::GetRunManager()->ReinitializeGeometry();
0330 }
0331 
0332 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0333 
0334 void LXeDetectorConstruction::SetNY(G4int ny)
0335 {
0336   fNy = ny;
0337   G4RunManager::GetRunManager()->ReinitializeGeometry();
0338 }
0339 
0340 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0341 
0342 void LXeDetectorConstruction::SetNZ(G4int nz)
0343 {
0344   fNz = nz;
0345   G4RunManager::GetRunManager()->ReinitializeGeometry();
0346 }
0347 
0348 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0349 
0350 void LXeDetectorConstruction::SetPMTRadius(G4double outerRadius_pmt)
0351 {
0352   fOuterRadius_pmt = outerRadius_pmt;
0353   G4RunManager::GetRunManager()->ReinitializeGeometry();
0354 }
0355 
0356 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0357 
0358 void LXeDetectorConstruction::SetDefaults()
0359 {
0360   // Resets to default values
0361   fD_mtl = 0.0635 * cm;
0362 
0363   fScint_x = 17.8 * cm;
0364   fScint_y = 17.8 * cm;
0365   fScint_z = 22.6 * cm;
0366 
0367   fNx = 2;
0368   fNy = 2;
0369   fNz = 3;
0370 
0371   fOuterRadius_pmt = 2.3 * cm;
0372 
0373   fSphereOn = true;
0374   fRefl = 1.0;
0375 
0376   fNfibers = 15;
0377   fWLSslab = false;
0378   fMainVolumeOn = true;
0379   fMainVolume = nullptr;
0380   fSlab_z = 2.5 * mm;
0381 
0382   G4UImanager::GetUIpointer()->ApplyCommand("/LXe/detector/scintYieldFactor 1.");
0383 
0384   if (fLXe_mt) fLXe_mt->AddConstProperty("SCINTILLATIONYIELD", 12000. / MeV);
0385   if (fMPTPStyrene) fMPTPStyrene->AddConstProperty("SCINTILLATIONYIELD", 10. / keV);
0386 }
0387 
0388 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0389 
0390 void LXeDetectorConstruction::SetSphereOn(G4bool b)
0391 {
0392   fSphereOn = b;
0393   G4RunManager::GetRunManager()->ReinitializeGeometry();
0394 }
0395 
0396 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0397 
0398 void LXeDetectorConstruction::SetHousingReflectivity(G4double r)
0399 {
0400   fRefl = r;
0401   G4RunManager::GetRunManager()->ReinitializeGeometry();
0402 }
0403 
0404 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0405 
0406 void LXeDetectorConstruction::SetWLSSlabOn(G4bool b)
0407 {
0408   fWLSslab = b;
0409   G4RunManager::GetRunManager()->ReinitializeGeometry();
0410 }
0411 
0412 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0413 
0414 void LXeDetectorConstruction::SetMainVolumeOn(G4bool b)
0415 {
0416   fMainVolumeOn = b;
0417   G4RunManager::GetRunManager()->ReinitializeGeometry();
0418 }
0419 
0420 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0421 
0422 void LXeDetectorConstruction::SetNFibers(G4int n)
0423 {
0424   fNfibers = n;
0425   G4RunManager::GetRunManager()->ReinitializeGeometry();
0426 }
0427 
0428 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0429 
0430 void LXeDetectorConstruction::SetMainScintYield(G4double y)
0431 {
0432   fLXe_mt->AddConstProperty("SCINTILLATIONYIELD", y / MeV);
0433 }
0434 
0435 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0436 
0437 void LXeDetectorConstruction::SetWLSScintYield(G4double y)
0438 {
0439   fMPTPStyrene->AddConstProperty("SCINTILLATIONYIELD", y / MeV);
0440 }
0441 
0442 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0443 
0444 void LXeDetectorConstruction::SetSaveThreshold(G4int save)
0445 {
0446   // Sets the save threshold for the random number seed. If the number of
0447   // photons generated in an event is lower than this, then save the seed for
0448   // this event in a file called run###evt###.rndm
0449 
0450   fSaveThreshold = save;
0451   G4RunManager::GetRunManager()->SetRandomNumberStore(true);
0452 }
0453 
0454 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......