OpNovice/src/OpNoviceDetectorConstruction.cc

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 OpNoviceDetectorConstruction.cc
0027 /// \brief Implementation of the OpNoviceDetectorConstruction class
0028
0029 #include "OpNoviceDetectorConstruction.hh"
0030
0031 #include "OpNoviceDetectorMessenger.hh"
0032
0033 #include "G4Box.hh"
0034 #include "G4Element.hh"
0035 #include "G4GDMLParser.hh"
0036 #include "G4LogicalBorderSurface.hh"
0037 #include "G4LogicalSkinSurface.hh"
0038 #include "G4LogicalVolume.hh"
0039 #include "G4Material.hh"
0040 #include "G4OpticalSurface.hh"
0041 #include "G4PVPlacement.hh"
0042 #include "G4SystemOfUnits.hh"
0043 #include "G4ThreeVector.hh"
0044
0045 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0046 OpNoviceDetectorConstruction::OpNoviceDetectorConstruction() : G4VUserDetectorConstruction()
0047 {
0048   // create a messenger for this class
0049   fDetectorMessenger = new OpNoviceDetectorMessenger(this);
0050 }
0051
0052 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0053 OpNoviceDetectorConstruction::~OpNoviceDetectorConstruction()
0054 {
0055   delete fDetectorMessenger;
0056 }
0057
0058 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0059 G4VPhysicalVolume* OpNoviceDetectorConstruction::Construct()
0060 {
0061   G4bool checkOverlaps = true;
0062   // ------------- Materials -------------
0063   G4double a, z, density;
0064   G4int nelements;
0065
0066   // Air
0067   auto N = new G4Element("Nitrogen", "N", z = 7, a = 14.01 * g / mole);
0068   auto O = new G4Element("Oxygen", "O", z = 8, a = 16.00 * g / mole);
0069   auto air = new G4Material("Air", density = 1.29 * mg / cm3, nelements = 2);
0070   air->AddElement(N, 70. * perCent);
0071   air->AddElement(O, 30. * perCent);
0072   //
0073   // Water
0074   auto H = new G4Element("Hydrogen", "H", z = 1, a = 1.01 * g / mole);
0075   auto water = new G4Material("Water", density = 1.0 * g / cm3, nelements = 2);
0076   water->AddElement(H, 2);
0077   water->AddElement(O, 1);
0078
0079   // ------------ Generate & Add Material Properties Table ------------
0080   //
0081   std::vector<G4double> photonEnergy = {
0082     2.034 * eV, 2.068 * eV, 2.103 * eV, 2.139 * eV, 2.177 * eV, 2.216 * eV, 2.256 * eV, 2.298 * eV,
0083     2.341 * eV, 2.386 * eV, 2.433 * eV, 2.481 * eV, 2.532 * eV, 2.585 * eV, 2.640 * eV, 2.697 * eV,
0084     2.757 * eV, 2.820 * eV, 2.885 * eV, 2.954 * eV, 3.026 * eV, 3.102 * eV, 3.181 * eV, 3.265 * eV,
0085     3.353 * eV, 3.446 * eV, 3.545 * eV, 3.649 * eV, 3.760 * eV, 3.877 * eV, 4.002 * eV, 4.136 * eV};
0086
0087   // Water
0088   std::vector<G4double> refractiveIndex1 = {
0089     1.3435, 1.344, 1.3445, 1.345,  1.3455, 1.346,  1.3465, 1.347,  1.3475, 1.348,  1.3485,
0090     1.3492, 1.35,  1.3505, 1.351,  1.3518, 1.3522, 1.3530, 1.3535, 1.354,  1.3545, 1.355,
0091     1.3555, 1.356, 1.3568, 1.3572, 1.358,  1.3585, 1.359,  1.3595, 1.36,   1.3608};
0092   std::vector<G4double> absorption = {
0093     3.448 * m,  4.082 * m,  6.329 * m,  9.174 * m,  12.346 * m, 13.889 * m, 15.152 * m, 17.241 * m,
0094     18.868 * m, 20.000 * m, 26.316 * m, 35.714 * m, 45.455 * m, 47.619 * m, 52.632 * m, 52.632 * m,
0095     55.556 * m, 52.632 * m, 52.632 * m, 47.619 * m, 45.455 * m, 41.667 * m, 37.037 * m, 33.333 * m,
0096     30.000 * m, 28.500 * m, 27.000 * m, 24.500 * m, 22.000 * m, 19.500 * m, 17.500 * m, 14.500 * m};
0097
0098   // Material properties can be added as arrays. However, in this case it is
0099   // up to the user to make sure both arrays have the same number of elements.
0100   G4double scintilFastArray[]{1.0, 1.0};
0101   G4double energyArray[]{2.034 * eV, 4.136 * eV};
0102   G4int lenArray = 2;
0103
0104   std::vector<G4double> scintilSlow = {
0105     0.01, 1.00, 2.00, 3.00, 4.00, 5.00, 6.00, 7.00, 8.00, 9.00, 8.00, 7.00, 6.00, 4.00, 3.00, 2.00,
0106     1.00, 0.01, 1.00, 2.00, 3.00, 4.00, 5.00, 6.00, 7.00, 8.00, 9.00, 8.00, 7.00, 6.00, 5.00, 4.00};
0107
0108   auto myMPT1 = new G4MaterialPropertiesTable();
0109
0110   // Values can be added to the material property table individually.
0111   // With this method, spline interpolation cannot be set. Arguments
0112   // createNewKey and spline both take their default values of false.
0113   // Need to specify the number of entries (1) in the arrays, as an argument
0114   // to AddProperty.
0115   G4int numEntries = 1;
0116   myMPT1->AddProperty("RINDEX", &photonEnergy[0], &refractiveIndex1[0], numEntries);
0117
0118   for (size_t i = 1; i < photonEnergy.size(); ++i) {
0119     myMPT1->AddEntry("RINDEX", photonEnergy[i], refractiveIndex1[i]);
0120   }
0121
0122   // Check that group velocity is calculated from RINDEX
0123   if (myMPT1->GetProperty("RINDEX")->GetVectorLength()
0124       != myMPT1->GetProperty("GROUPVEL")->GetVectorLength())
0125   {
0126     G4ExceptionDescription ed;
0127     ed << "Error calculating group velocities. Incorrect number of entries "
0128           "in group velocity material property vector.";
0129     G4Exception("OpNovice::OpNoviceDetectorConstruction", "OpNovice001", FatalException, ed);
0130   }
0131
0132   // Adding a property from two std::vectors. Argument createNewKey is false
0133   // and spline is true.
0134   myMPT1->AddProperty("ABSLENGTH", photonEnergy, absorption, false, true);
0135
0136   // Adding a property using a C-style array.
0137   // Spline interpolation isn't used for scintillation.
0138   // Arguments spline and createNewKey both take default value false.
0139   myMPT1->AddProperty("SCINTILLATIONCOMPONENT1", energyArray, scintilFastArray, lenArray);
0140
0141   myMPT1->AddProperty("SCINTILLATIONCOMPONENT2", photonEnergy, scintilSlow, false, true);
0142   myMPT1->AddConstProperty("SCINTILLATIONYIELD", 50. / MeV);
0143   myMPT1->AddConstProperty("RESOLUTIONSCALE", 1.0);
0144   myMPT1->AddConstProperty("SCINTILLATIONTIMECONSTANT1", 1. * ns);
0145   myMPT1->AddConstProperty("SCINTILLATIONTIMECONSTANT2", 10. * ns);
0146   myMPT1->AddConstProperty("SCINTILLATIONYIELD1", 0.8);
0147   myMPT1->AddConstProperty("SCINTILLATIONYIELD2", 0.2);
0148   std::vector<G4double> energy_water = {
0149     1.56962 * eV, 1.58974 * eV, 1.61039 * eV, 1.63157 * eV, 1.65333 * eV, 1.67567 * eV,
0150     1.69863 * eV, 1.72222 * eV, 1.74647 * eV, 1.77142 * eV, 1.7971 * eV,  1.82352 * eV,
0151     1.85074 * eV, 1.87878 * eV, 1.90769 * eV, 1.93749 * eV, 1.96825 * eV, 1.99999 * eV,
0152     2.03278 * eV, 2.06666 * eV, 2.10169 * eV, 2.13793 * eV, 2.17543 * eV, 2.21428 * eV,
0153     2.25454 * eV, 2.29629 * eV, 2.33962 * eV, 2.38461 * eV, 2.43137 * eV, 2.47999 * eV,
0154     2.53061 * eV, 2.58333 * eV, 2.63829 * eV, 2.69565 * eV, 2.75555 * eV, 2.81817 * eV,
0155     2.88371 * eV, 2.95237 * eV, 3.02438 * eV, 3.09999 * eV, 3.17948 * eV, 3.26315 * eV,
0156     3.35134 * eV, 3.44444 * eV, 3.54285 * eV, 3.64705 * eV, 3.75757 * eV, 3.87499 * eV,
0157     3.99999 * eV, 4.13332 * eV, 4.27585 * eV, 4.42856 * eV, 4.59258 * eV, 4.76922 * eV,
0158     4.95999 * eV, 5.16665 * eV, 5.39129 * eV, 5.63635 * eV, 5.90475 * eV, 6.19998 * eV};
0159
0160   // Rayleigh scattering length is calculated by G4OpRayleigh
0161
0162   // Mie: assume 100 times larger than the rayleigh scattering
0163   std::vector<G4double> mie_water = {
0164     167024.4 * m, 158726.7 * m, 150742 * m,   143062.5 * m, 135680.2 * m, 128587.4 * m,
0165     121776.3 * m, 115239.5 * m, 108969.5 * m, 102958.8 * m, 97200.35 * m, 91686.86 * m,
0166     86411.33 * m, 81366.79 * m, 76546.42 * m, 71943.46 * m, 67551.29 * m, 63363.36 * m,
0167     59373.25 * m, 55574.61 * m, 51961.24 * m, 48527.00 * m, 45265.87 * m, 42171.94 * m,
0168     39239.39 * m, 36462.50 * m, 33835.68 * m, 31353.41 * m, 29010.30 * m, 26801.03 * m,
0169     24720.42 * m, 22763.36 * m, 20924.88 * m, 19200.07 * m, 17584.16 * m, 16072.45 * m,
0170     14660.38 * m, 13343.46 * m, 12117.33 * m, 10977.70 * m, 9920.416 * m, 8941.407 * m,
0171     8036.711 * m, 7202.470 * m, 6434.927 * m, 5730.429 * m, 5085.425 * m, 4496.467 * m,
0172     3960.210 * m, 3473.413 * m, 3032.937 * m, 2635.746 * m, 2278.907 * m, 1959.588 * m,
0173     1675.064 * m, 1422.710 * m, 1200.004 * m, 1004.528 * m, 833.9666 * m, 686.1063 * m};
0174
0175   // Mie: gforward, gbackward, forward backward ratio
0176   G4double mie_water_const[3] = {0.99, 0.99, 0.8};
0177
0178   myMPT1->AddProperty("MIEHG", energy_water, mie_water, false, true);
0179   myMPT1->AddConstProperty("MIEHG_FORWARD", mie_water_const[0]);
0180   myMPT1->AddConstProperty("MIEHG_BACKWARD", mie_water_const[1]);
0181   myMPT1->AddConstProperty("MIEHG_FORWARD_RATIO", mie_water_const[2]);
0182
0183   G4cout << "Water G4MaterialPropertiesTable:" << G4endl;
0184   myMPT1->DumpTable();
0185
0186   water->SetMaterialPropertiesTable(myMPT1);
0187
0188   // Set the Birks Constant for the Water scintillator
0189   water->GetIonisation()->SetBirksConstant(0.126 * mm / MeV);
0190
0191   // Air
0192   std::vector<G4double> refractiveIndex2 = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
0193                                             1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
0194                                             1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
0195
0196   auto myMPT2 = new G4MaterialPropertiesTable();
0197   myMPT2->AddProperty("RINDEX", photonEnergy, refractiveIndex2);
0198
0199   G4cout << "Air G4MaterialPropertiesTable:" << G4endl;
0200   myMPT2->DumpTable();
0201
0202   air->SetMaterialPropertiesTable(myMPT2);
0203
0204   // ------------- Volumes --------------
0205   //
0206   // The world
0207   auto world_box = new G4Box("World", fWorld_x, fWorld_y, fWorld_z);
0208   auto world_log = new G4LogicalVolume(world_box, air, "World");
0209   G4VPhysicalVolume* world_phys = new G4PVPlacement(nullptr, G4ThreeVector(), world_log, "world",
0210                                                     nullptr, false, 0, checkOverlaps);
0211
0212   // The experimental Hall
0213   auto expHall_box = new G4Box("expHall", fExpHall_x, fExpHall_y, fExpHall_z);
0214   auto expHall_log = new G4LogicalVolume(expHall_box, air, "expHall");
0215   G4VPhysicalVolume* expHall_phys =
0216     new G4PVPlacement(nullptr, G4ThreeVector(), expHall_log, "expHall", world_log, false, 0);
0217
0218   // The Water Tank
0219   auto waterTank_box = new G4Box("Tank", fTank_x, fTank_y, fTank_z);
0220   auto waterTank_log = new G4LogicalVolume(waterTank_box, water, "Tank");
0221   G4VPhysicalVolume* waterTank_phys =
0222     new G4PVPlacement(nullptr, G4ThreeVector(), waterTank_log, "Tank", expHall_log, false, 0);
0223
0224   // The Air Bubble
0225   auto bubbleAir_box = new G4Box("Bubble", fBubble_x, fBubble_y, fBubble_z);
0226   auto bubbleAir_log = new G4LogicalVolume(bubbleAir_box, air, "Bubble");
0227   new G4PVPlacement(nullptr, G4ThreeVector(0, 2.5 * m, 0), bubbleAir_log, "Bubble", waterTank_log,
0228                     false, 0);
0229
0230   // ------------- Surfaces --------------
0231
0232   // Water Tank
0233   auto opWaterSurface = new G4OpticalSurface("WaterSurface");
0234   opWaterSurface->SetType(dielectric_LUTDAVIS);
0235   opWaterSurface->SetFinish(Rough_LUT);
0236   opWaterSurface->SetModel(DAVIS);
0237
0238   auto waterSurface =
0239     new G4LogicalBorderSurface("WaterSurface", waterTank_phys, expHall_phys, opWaterSurface);
0240
0241   auto opticalSurface = dynamic_cast<G4OpticalSurface*>(
0242     waterSurface->GetSurface(waterTank_phys, expHall_phys)->GetSurfaceProperty());
0243   if (opticalSurface) opticalSurface->DumpInfo();
0244
0245   // Air Bubble
0246   auto opAirSurface = new G4OpticalSurface("AirSurface");
0247   opAirSurface->SetType(dielectric_dielectric);
0248   opAirSurface->SetFinish(polished);
0249   opAirSurface->SetModel(glisur);
0250
0251   auto airSurface = new G4LogicalSkinSurface("AirSurface", bubbleAir_log, opAirSurface);
0252
0253   opticalSurface =
0254     dynamic_cast<G4OpticalSurface*>(airSurface->GetSurface(bubbleAir_log)->GetSurfaceProperty());
0255   if (opticalSurface) opticalSurface->DumpInfo();
0256
0257   // Generate & Add Material Properties Table attached to the optical surfaces
0258   //
0259   std::vector<G4double> ephoton = {2.034 * eV, 4.136 * eV};
0260
0261   // OpticalAirSurface
0262   std::vector<G4double> reflectivity = {0.3, 0.5};
0263   std::vector<G4double> efficiency = {0.8, 1.0};
0264
0265   auto myST2 = new G4MaterialPropertiesTable();
0266
0267   myST2->AddProperty("REFLECTIVITY", ephoton, reflectivity);
0268   myST2->AddProperty("EFFICIENCY", ephoton, efficiency);
0269   if (fVerbose) {
0270     G4cout << "Air Surface G4MaterialPropertiesTable:" << G4endl;
0271     myST2->DumpTable();
0272   }
0273   opAirSurface->SetMaterialPropertiesTable(myST2);
0274
0275   if (fDumpGdml) {
0276     auto parser = new G4GDMLParser();
0277     parser->Write(fDumpGdmlFileName, world_phys);
0278   }
0279
0280   ////////////////////////////////////////////////////////////////////////////
0281   // test user-defined properties
0282   G4String ed;
0283   if (myMPT1->GetProperty("USERDEFINED") != nullptr) {
0284     ed = "USERDEFINED != nullptr";
0285     PrintError(ed);
0286   }
0287   myMPT1->AddProperty("USERDEFINED", energy_water, mie_water, true, true);
0288   if (myMPT1->GetProperty("USERDEFINED") == nullptr) {
0289     ed = "USERDEFINED == nullptr";
0290     PrintError(ed);
0291   }
0292   [[maybe_unused]] G4int index_userdefined = -1;
0293   if (myMPT1->GetProperty("USERDEFINED") != nullptr) {
0294     index_userdefined = myMPT1->GetPropertyIndex("USERDEFINED");
0295   }
0296   if (!(index_userdefined >= 0
0297         && index_userdefined < (G4int)myMPT1->GetMaterialPropertyNames().size()))
0298   {
0299     ed = "USERDEFINED index out of range";
0300     PrintError(ed);
0301   }
0302   myMPT1->RemoveProperty("USERDEFINED");
0303   if (myMPT1->GetProperty("USERDEFINED") != nullptr) {
0304     ed = "USERDEFINED != nullptr at end";
0305     PrintError(ed);
0306   }
0307
0308   if (myMPT1->ConstPropertyExists("USERDEFINEDCONST") == true) {
0309     ed = "ConstProperty USERDEFINEDCONST already exists.";
0310     PrintError(ed);
0311   }
0312   myMPT1->AddConstProperty("USERDEFINEDCONST", 3.14, true);
0313   if (myMPT1->ConstPropertyExists("USERDEFINEDCONST") == false) {
0314     ed = "ConstProperty USERDEFINEDCONST doesn't exist.";
0315     PrintError(ed);
0316   }
0317   [[maybe_unused]] G4int index_pi = -1;
0318   if (myMPT1->ConstPropertyExists("USERDEFINEDCONST") == true) {
0319     index_pi = myMPT1->GetConstPropertyIndex("USERDEFINEDCONST");
0320   }
0321   if (!(index_pi >= 0 && index_pi < (G4int)myMPT1->GetMaterialConstPropertyNames().size())) {
0322     ed = "ConstProperty USERDEFINEDCONST index out of range.";
0323     PrintError(ed);
0324   }
0325   myMPT1->RemoveConstProperty("USERDEFINEDCONST");
0326   if (myMPT1->ConstPropertyExists("USERDEFINEDCONST") == true) {
0327     ed = "ConstProperty USERDEFINEDCONST still exists.";
0328     PrintError(ed);
0329   }
0330   // done testing user-defined properties
0331   ////////////////////////////////////////////////////////////////////////////
0332
0333   return world_phys;
0334 }
0335
0336 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0337 void OpNoviceDetectorConstruction::SetDumpGdml(G4bool val)
0338 {
0339   fDumpGdml = val;
0340 }
0341
0342 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0343 G4bool OpNoviceDetectorConstruction::IsDumpGdml() const
0344 {
0345   return fDumpGdml;
0346 }
0347
0348 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0349 void OpNoviceDetectorConstruction::SetVerbose(G4bool val)
0350 {
0351   fVerbose = val;
0352 }
0353
0354 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0355 G4bool OpNoviceDetectorConstruction::IsVerbose() const
0356 {
0357   return fVerbose;
0358 }
0359
0360 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0361 void OpNoviceDetectorConstruction::SetDumpGdmlFile(G4String filename)
0362 {
0363   fDumpGdmlFileName = filename;
0364 }
0365
0366 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0367 G4String OpNoviceDetectorConstruction::GetDumpGdmlFile() const
0368 {
0369   return fDumpGdmlFileName;
0370 }
0371
0372 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0373 void OpNoviceDetectorConstruction::PrintError(G4String ed)
0374 {
0375   G4Exception("OpNoviceDetectorConstruction:MaterialProperty test", "op001", FatalException, ed);
0376 }