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