B01/src/B01DetectorConstruction.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 biasing/B01/src/B01DetectorConstruction.cc
0027 /// \brief Implementation of the B01DetectorConstruction class
0028 //
0029 //
0030 //
0031
0032 #include "B01DetectorConstruction.hh"
0033
0034 #include "G4Box.hh"
0035 #include "G4Colour.hh"
0036 #include "G4LogicalVolume.hh"
0037 #include "G4Material.hh"
0038 #include "G4PVPlacement.hh"
0039 #include "G4PhysicalConstants.hh"
0040 #include "G4SystemOfUnits.hh"
0041 #include "G4ThreeVector.hh"
0042 #include "G4Tubs.hh"
0043 #include "G4Types.hh"
0044 #include "G4VisAttributes.hh"
0045 #include "globals.hh"
0046
0047 #include <set>
0048 #include <sstream>
0049
0050 // For Primitive Scorers
0051 #include "G4MultiFunctionalDetector.hh"
0052 #include "G4PSNofCollision.hh"
0053 #include "G4PSPopulation.hh"
0054 #include "G4PSTrackCounter.hh"
0055 #include "G4PSTrackLength.hh"
0056 #include "G4SDManager.hh"
0057 #include "G4SDParticleFilter.hh"
0058
0059 // for importance biasing
0060 #include "G4IStore.hh"
0061
0062 // for weight window technique
0063 #include "G4WeightWindowStore.hh"
0064
0065 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0066
0067 B01DetectorConstruction::B01DetectorConstruction()
0068   : G4VUserDetectorConstruction(), fLogicalVolumeVector(), fPhysicalVolumeVector()
0069 {
0070   ;
0071 }
0072
0073 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0074
0075 B01DetectorConstruction::~B01DetectorConstruction()
0076 {
0077   fLogicalVolumeVector.clear();
0078   fPhysicalVolumeVector.clear();
0079 }
0080
0081 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0082
0083 G4VPhysicalVolume* B01DetectorConstruction::Construct()
0084 {
0085   G4double pos_x;
0086   G4double pos_y;
0087   G4double pos_z;
0088
0089   G4double density, pressure, temperature;
0090   G4double A;
0091   G4int Z;
0092
0093   G4String name, symbol;
0094   G4double z;
0095   G4double fractionmass;
0096
0097   A = 1.01 * g / mole;
0098   G4Element* elH = new G4Element(name = "Hydrogen", symbol = "H", Z = 1, A);
0099
0100   A = 12.01 * g / mole;
0101   G4Element* elC = new G4Element(name = "Carbon", symbol = "C", Z = 6, A);
0102
0103   A = 16.00 * g / mole;
0104   G4Element* elO = new G4Element(name = "Oxygen", symbol = "O", Z = 8, A);
0105
0106   A = 22.99 * g / mole;
0107   G4Element* elNa = new G4Element(name = "Natrium", symbol = "Na", Z = 11, A);
0108
0109   A = 200.59 * g / mole;
0110   G4Element* elHg = new G4Element(name = "Hg", symbol = "Hg", Z = 80, A);
0111
0112   A = 26.98 * g / mole;
0113   G4Element* elAl = new G4Element(name = "Aluminium", symbol = "Al", Z = 13, A);
0114
0115   A = 28.09 * g / mole;
0116   G4Element* elSi = new G4Element(name = "Silicon", symbol = "Si", Z = 14, A);
0117
0118   A = 39.1 * g / mole;
0119   G4Element* elK = new G4Element(name = "K", symbol = "K", Z = 19, A);
0120
0121   A = 69.72 * g / mole;
0122   G4Element* elCa = new G4Element(name = "Calzium", symbol = "Ca", Z = 31, A);
0123
0124   A = 55.85 * g / mole;
0125   G4Element* elFe = new G4Element(name = "Iron", symbol = "Fe", Z = 26, A);
0126
0127   density = universe_mean_density;  // from PhysicalConstants.h
0128   pressure = 3.e-18 * pascal;
0129   temperature = 2.73 * kelvin;
0130   G4Material* Galactic = new G4Material(name = "Galactic", z = 1., A = 1.01 * g / mole, density,
0131                                         kStateGas, temperature, pressure);
0132
0133   density = 2.03 * g / cm3;
0134   G4Material* Concrete = new G4Material("Concrete", density, 10);
0135   Concrete->AddElement(elH, fractionmass = 0.01);
0136   Concrete->AddElement(elO, fractionmass = 0.529);
0137   Concrete->AddElement(elNa, fractionmass = 0.016);
0138   Concrete->AddElement(elHg, fractionmass = 0.002);
0139   Concrete->AddElement(elAl, fractionmass = 0.034);
0140   Concrete->AddElement(elSi, fractionmass = 0.337);
0141   Concrete->AddElement(elK, fractionmass = 0.013);
0142   Concrete->AddElement(elCa, fractionmass = 0.044);
0143   Concrete->AddElement(elFe, fractionmass = 0.014);
0144   Concrete->AddElement(elC, fractionmass = 0.001);
0145
0146   /////////////////////////////
0147   // world cylinder volume
0148   ////////////////////////////
0149
0150   // world solid
0151
0152   G4double innerRadiusCylinder = 0 * cm;
0153   G4double outerRadiusCylinder = 100 * cm;
0154   G4double heightCylinder = 100 * cm;
0155   G4double startAngleCylinder = 0 * deg;
0156   G4double spanningAngleCylinder = 360 * deg;
0157
0158   G4Tubs* worldCylinder = new G4Tubs("worldCylinder", innerRadiusCylinder, outerRadiusCylinder,
0159                                      heightCylinder, startAngleCylinder, spanningAngleCylinder);
0160
0161   // logical world
0162
0163   G4LogicalVolume* worldCylinder_log =
0164     new G4LogicalVolume(worldCylinder, Galactic, "worldCylinder_log");
0165   fLogicalVolumeVector.push_back(worldCylinder_log);
0166
0167   name = "shieldWorld";
0168   fWorldVolume = new G4PVPlacement(0, G4ThreeVector(0, 0, 0), worldCylinder_log, name, 0, false, 0);
0169
0170   fPhysicalVolumeVector.push_back(fWorldVolume);
0171
0172   // creating 18 slabs of 10 cm thick concrete
0173
0174   G4double innerRadiusShield = 0 * cm;
0175   G4double outerRadiusShield = 100 * cm;
0176   G4double heightShield = 5 * cm;
0177   G4double startAngleShield = 0 * deg;
0178   G4double spanningAngleShield = 360 * deg;
0179
0180   G4Tubs* aShield = new G4Tubs("aShield", innerRadiusShield, outerRadiusShield, heightShield,
0181                                startAngleShield, spanningAngleShield);
0182
0183   // logical shield
0184
0185   G4LogicalVolume* aShield_log = new G4LogicalVolume(aShield, Concrete, "aShield_log");
0186   fLogicalVolumeVector.push_back(aShield_log);
0187
0188   G4VisAttributes* pShieldVis = new G4VisAttributes(G4Colour(0.0, 0.0, 1.0));
0189   pShieldVis->SetForceSolid(true);
0190   aShield_log->SetVisAttributes(pShieldVis);
0191
0192   // physical shields
0193
0194   G4int i;
0195   G4double startz = -85 * cm;
0196   for (i = 1; i <= 18; i++) {
0197     name = GetCellName(i);
0198     pos_x = 0 * cm;
0199     pos_y = 0 * cm;
0200     pos_z = startz + (i - 1) * (2 * heightShield);
0201     G4VPhysicalVolume* pvol = new G4PVPlacement(0, G4ThreeVector(pos_x, pos_y, pos_z), aShield_log,
0202                                                 name, worldCylinder_log, false, i);
0203     fPhysicalVolumeVector.push_back(pvol);
0204   }
0205
0206   // filling the rest of the world volume behind the concrete with
0207   // another slab which should get the same importance value
0208   // or lower weight bound as the last slab
0209   //
0210   innerRadiusShield = 0 * cm;
0211   outerRadiusShield = 100 * cm;
0212   heightShield = 5 * cm;
0213   startAngleShield = 0 * deg;
0214   spanningAngleShield = 360 * deg;
0215
0216   G4Tubs* aRest = new G4Tubs("Rest", innerRadiusShield, outerRadiusShield, heightShield,
0217                              startAngleShield, spanningAngleShield);
0218
0219   G4LogicalVolume* aRest_log = new G4LogicalVolume(aRest, Galactic, "aRest_log");
0220   fLogicalVolumeVector.push_back(aRest_log);
0221   name = "rest";
0222
0223   pos_x = 0 * cm;
0224   pos_y = 0 * cm;
0225   pos_z = 95 * cm;
0226   G4VPhysicalVolume* pvol_rest = new G4PVPlacement(0, G4ThreeVector(pos_x, pos_y, pos_z), aRest_log,
0227                                                    name, worldCylinder_log, false,
0228                                                    19);  // i=19
0229
0230   fPhysicalVolumeVector.push_back(pvol_rest);
0231
0232   SetSensitive();
0233   return fWorldVolume;
0234 }
0235
0236 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0237
0238 G4VIStore* B01DetectorConstruction::CreateImportanceStore()
0239 {
0240   G4cout << " B01DetectorConstruction:: Creating Importance Store " << G4endl;
0241   if (!fPhysicalVolumeVector.size()) {
0242     G4Exception("B01DetectorConstruction::CreateImportanceStore", "exampleB01_0001",
0243                 RunMustBeAborted, "no physical volumes created yet!");
0244   }
0245
0246   fWorldVolume = fPhysicalVolumeVector[0];
0247
0248   // creating and filling the importance store
0249
0250   G4IStore* istore = G4IStore::GetInstance();
0251
0252   G4int n = 0;
0253   G4double imp = 1;
0254   istore->AddImportanceGeometryCell(1, *fWorldVolume);
0255   for (std::vector<G4VPhysicalVolume*>::iterator it = fPhysicalVolumeVector.begin();
0256        it != fPhysicalVolumeVector.end() - 1; it++)
0257   {
0258     if (*it != fWorldVolume) {
0259       imp = std::pow(2., n++);
0260       G4cout << "Going to assign importance: " << imp << ", to volume: " << (*it)->GetName()
0261              << G4endl;
0262       istore->AddImportanceGeometryCell(imp, *(*it), n);
0263     }
0264   }
0265
0266   // the remaining part pf the geometry (rest) gets the same
0267   // importance as the last conrete cell
0268   //
0269   istore->AddImportanceGeometryCell(imp, *(fPhysicalVolumeVector[fPhysicalVolumeVector.size() - 1]),
0270                                     ++n);
0271
0272   return istore;
0273 }
0274
0275 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0276
0277 G4VWeightWindowStore* B01DetectorConstruction::CreateWeightWindowStore()
0278 {
0279   if (!fPhysicalVolumeVector.size()) {
0280     G4Exception("B01DetectorConstruction::CreateWeightWindowStore", "exampleB01_0002",
0281                 RunMustBeAborted, "no physical volumes created yet!");
0282   }
0283
0284   fWorldVolume = fPhysicalVolumeVector[0];
0285
0286   // creating and filling the weight window store
0287
0288   G4WeightWindowStore* wwstore = G4WeightWindowStore::GetInstance();
0289
0290   // create one energy region covering the energies of the problem
0291   //
0292   std::set<G4double, std::less<G4double>> enBounds;
0293   enBounds.insert(1 * GeV);
0294   wwstore->SetGeneralUpperEnergyBounds(enBounds);
0295
0296   G4int n = 0;
0297   G4double lowerWeight = 1;
0298   std::vector<G4double> lowerWeights;
0299
0300   lowerWeights.push_back(1);
0301   G4GeometryCell gWorldCell(*fWorldVolume, 0);
0302   wwstore->AddLowerWeights(gWorldCell, lowerWeights);
0303
0304   for (std::vector<G4VPhysicalVolume*>::iterator it = fPhysicalVolumeVector.begin();
0305        it != fPhysicalVolumeVector.end() - 1; it++)
0306   {
0307     if (*it != fWorldVolume) {
0308       lowerWeight = 1. / std::pow(2., n++);
0309       G4cout << "Going to assign lower weight: " << lowerWeight
0310              << ", to volume: " << (*it)->GetName() << G4endl;
0311       G4GeometryCell gCell(*(*it), n);
0312       lowerWeights.clear();
0313       lowerWeights.push_back(lowerWeight);
0314       wwstore->AddLowerWeights(gCell, lowerWeights);
0315     }
0316   }
0317
0318   // the remaining part pf the geometry (rest) gets the same
0319   // lower weight bound  as the last conrete cell
0320   //
0321   G4GeometryCell gRestCell(*(fPhysicalVolumeVector[fPhysicalVolumeVector.size() - 1]), ++n);
0322   wwstore->AddLowerWeights(gRestCell, lowerWeights);
0323
0324   return wwstore;
0325 }
0326
0327 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0328
0329 G4String B01DetectorConstruction::GetCellName(G4int i)
0330 {
0331   std::ostringstream os;
0332   os << "cell_";
0333   if (i < 10) {
0334     os << "0";
0335   }
0336   os << i;
0337   G4String name = os.str();
0338   return name;
0339 }
0340
0341 G4VPhysicalVolume* B01DetectorConstruction::GetWorldVolume()
0342 {
0343   return fWorldVolume;
0344 }
0345
0346 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0347
0348 void B01DetectorConstruction::SetSensitive()
0349 {
0350   //  -------------------------------------------------
0351   //   The collection names of defined Primitives are
0352   //   0       ConcreteSD/Collisions
0353   //   1       ConcreteSD/CollWeight
0354   //   2       ConcreteSD/Population
0355   //   3       ConcreteSD/TrackEnter
0356   //   4       ConcreteSD/SL
0357   //   5       ConcreteSD/SLW
0358   //   6       ConcreteSD/SLWE
0359   //   7       ConcreteSD/SLW_V
0360   //   8       ConcreteSD/SLWE_V
0361   //  -------------------------------------------------
0362
0363   // moved to ConstructSDandField() for MT compliance
0364 }
0365
0366 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0367
0368 void B01DetectorConstruction::ConstructSDandField()
0369 {
0370   //  Sensitive Detector Manager.
0371   G4SDManager* SDman = G4SDManager::GetSDMpointer();
0372   // Sensitive Detector Name
0373   G4String concreteSDname = "ConcreteSD";
0374
0375   //------------------------
0376   // MultiFunctionalDetector
0377   //------------------------
0378   //
0379   // Define MultiFunctionalDetector with name.
0380   G4MultiFunctionalDetector* MFDet = new G4MultiFunctionalDetector(concreteSDname);
0381   SDman->AddNewDetector(MFDet);  // Register SD to SDManager
0382
0383   G4String fltName, particleName;
0384   G4SDParticleFilter* neutronFilter =
0385     new G4SDParticleFilter(fltName = "neutronFilter", particleName = "neutron");
0386
0387   MFDet->SetFilter(neutronFilter);
0388
0389   for (std::vector<G4LogicalVolume*>::iterator it = fLogicalVolumeVector.begin();
0390        it != fLogicalVolumeVector.end(); it++)
0391   {
0392     //      (*it)->SetSensitiveDetector(MFDet);
0393     SetSensitiveDetector((*it)->GetName(), MFDet);
0394   }
0395
0396   G4String psName;
0397   G4PSNofCollision* scorer0 = new G4PSNofCollision(psName = "Collisions");
0398   MFDet->RegisterPrimitive(scorer0);
0399
0400   G4PSNofCollision* scorer1 = new G4PSNofCollision(psName = "CollWeight");
0401   scorer1->Weighted(true);
0402   MFDet->RegisterPrimitive(scorer1);
0403
0404   G4PSPopulation* scorer2 = new G4PSPopulation(psName = "Population");
0405   MFDet->RegisterPrimitive(scorer2);
0406
0407   G4PSTrackCounter* scorer3 = new G4PSTrackCounter(psName = "TrackEnter", fCurrent_In);
0408   MFDet->RegisterPrimitive(scorer3);
0409
0410   G4PSTrackLength* scorer4 = new G4PSTrackLength(psName = "SL");
0411   MFDet->RegisterPrimitive(scorer4);
0412
0413   G4PSTrackLength* scorer5 = new G4PSTrackLength(psName = "SLW");
0414   scorer5->Weighted(true);
0415   MFDet->RegisterPrimitive(scorer5);
0416
0417   G4PSTrackLength* scorer6 = new G4PSTrackLength(psName = "SLWE");
0418   scorer6->Weighted(true);
0419   scorer6->MultiplyKineticEnergy(true);
0420   MFDet->RegisterPrimitive(scorer6);
0421
0422   G4PSTrackLength* scorer7 = new G4PSTrackLength(psName = "SLW_V");
0423   scorer7->Weighted(true);
0424   scorer7->DivideByVelocity(true);
0425   MFDet->RegisterPrimitive(scorer7);
0426
0427   G4PSTrackLength* scorer8 = new G4PSTrackLength(psName = "SLWE_V");
0428   scorer8->Weighted(true);
0429   scorer8->MultiplyKineticEnergy(true);
0430   scorer8->DivideByVelocity(true);
0431   MFDet->RegisterPrimitive(scorer8);
0432 }
0433
0434 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......