B02/src/B02ImportanceDetectorConstruction.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/B02/src/B02ImportanceDetectorConstruction.cc
0027 /// \brief Implementation of the B02ImportanceDetectorConstruction class
0028 //
0029 //
0030 //
0031
0032 #include "B02ImportanceDetectorConstruction.hh"
0033
0034 #include "G4LogicalVolume.hh"
0035 #include "G4Material.hh"
0036 #include "G4PVPlacement.hh"
0037 #include "G4PhysicalConstants.hh"
0038 #include "G4SystemOfUnits.hh"
0039 #include "G4ThreeVector.hh"
0040 #include "G4Tubs.hh"
0041 #include "globals.hh"
0042
0043 #include <sstream>
0044
0045 // For Primitive Scorers
0046 #include "G4MultiFunctionalDetector.hh"
0047 #include "G4PSNofCollision.hh"
0048 #include "G4PSPopulation.hh"
0049 #include "G4PSTrackCounter.hh"
0050 #include "G4PSTrackLength.hh"
0051 #include "G4SDManager.hh"
0052 #include "G4SDParticleFilter.hh"
0053
0054 // for importance biasing
0055 #include "G4IStore.hh"
0056
0057 // for weight window technique
0058 #include "G4WeightWindowStore.hh"
0059
0060 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0061
0062 B02ImportanceDetectorConstruction::B02ImportanceDetectorConstruction(G4String worldName)
0063   : G4VUserParallelWorld(worldName), fLogicalVolumeVector()
0064 {
0065   //  Construct();
0066 }
0067
0068 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0069
0070 B02ImportanceDetectorConstruction::~B02ImportanceDetectorConstruction()
0071 {
0072   fLogicalVolumeVector.clear();
0073 }
0074
0075 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0076
0077 void B02ImportanceDetectorConstruction::Construct()
0078 {
0079   G4cout << " constructing parallel world " << G4endl;
0080
0081   G4Material* dummyMat = 0;
0082
0083   // GetWorld methods create a clone of the mass world to the parallel world (!)
0084   //  via the transportation manager
0085   fGhostWorld = GetWorld();
0086   G4cout << " B02ImportanceDetectorConstruction:: ghostWorldName = " << fGhostWorld->GetName()
0087          << G4endl;
0088   G4LogicalVolume* worldLogical = fGhostWorld->GetLogicalVolume();
0089   fLogicalVolumeVector.push_back(worldLogical);
0090
0091   //  fPVolumeStore.AddPVolume(G4GeometryCell(*pWorldVolume, 0));
0092   fPVolumeStore.AddPVolume(G4GeometryCell(*fGhostWorld, 0));
0093
0094   // creating 18 slobs of 10 cm thicknes
0095
0096   G4double innerRadiusShield = 0 * cm;
0097   G4double outerRadiusShield = 100 * cm;
0098   G4double heightShield = 5 * cm;
0099   G4double startAngleShield = 0 * deg;
0100   G4double spanningAngleShield = 360 * deg;
0101
0102   G4Tubs* aShield = new G4Tubs("aShield", innerRadiusShield, outerRadiusShield, heightShield,
0103                                startAngleShield, spanningAngleShield);
0104
0105   // logical parallel cells
0106
0107   G4LogicalVolume* aShield_log_imp = new G4LogicalVolume(aShield, dummyMat, "aShield_log_imp");
0108   fLogicalVolumeVector.push_back(aShield_log_imp);
0109
0110   // physical parallel cells
0111   G4String name = "none";
0112   G4int i = 1;
0113   G4double startz = -85 * cm;
0114   //  for (i=1; i<=18; ++i) {
0115   for (i = 1; i <= 18; i++) {
0116     name = GetCellName(i);
0117
0118     G4double pos_x = 0 * cm;
0119     G4double pos_y = 0 * cm;
0120     G4double pos_z = startz + (i - 1) * (2 * heightShield);
0121     G4VPhysicalVolume* pvol = new G4PVPlacement(0, G4ThreeVector(pos_x, pos_y, pos_z),
0122                                                 aShield_log_imp, name, worldLogical, false, i);
0123     //                        0);
0124     G4GeometryCell cell(*pvol, i);
0125     //    G4GeometryCell cell(*pvol, 0);
0126     fPVolumeStore.AddPVolume(cell);
0127   }
0128
0129   // filling the rest of the world volumr behind the concrete with
0130   // another slob which should get the same importance value as the
0131   // last slob
0132   innerRadiusShield = 0 * cm;
0133   //  outerRadiusShield = 110*cm; exceeds world volume!!!!
0134   outerRadiusShield = 100 * cm;
0135   //  heightShield       = 10*cm;
0136   heightShield = 5 * cm;
0137   startAngleShield = 0 * deg;
0138   spanningAngleShield = 360 * deg;
0139
0140   G4Tubs* aRest = new G4Tubs("Rest", innerRadiusShield, outerRadiusShield, heightShield,
0141                              startAngleShield, spanningAngleShield);
0142
0143   G4LogicalVolume* aRest_log = new G4LogicalVolume(aRest, dummyMat, "aRest_log");
0144
0145   fLogicalVolumeVector.push_back(aRest_log);
0146
0147   name = GetCellName(19);
0148
0149   G4double pos_x = 0 * cm;
0150   G4double pos_y = 0 * cm;
0151   //  G4double pos_z = 100*cm;
0152   G4double pos_z = 95 * cm;
0153   G4VPhysicalVolume* pvol = new G4PVPlacement(0, G4ThreeVector(pos_x, pos_y, pos_z), aRest_log,
0154                                               name, worldLogical, false, 19);
0155   //                      0);
0156   G4GeometryCell cell(*pvol, 19);
0157   //  G4GeometryCell cell(*pvol, 0);
0158   fPVolumeStore.AddPVolume(cell);
0159
0160   SetSensitive();
0161 }
0162
0163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0164
0165 const G4VPhysicalVolume&
0166 B02ImportanceDetectorConstruction::GetPhysicalVolumeByName(const G4String& name) const
0167 {
0168   return *fPVolumeStore.GetPVolume(name);
0169 }
0170
0171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0172
0173 G4String B02ImportanceDetectorConstruction::ListPhysNamesAsG4String()
0174 {
0175   G4String names(fPVolumeStore.GetPNames());
0176   return names;
0177 }
0178
0179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0180
0181 G4String B02ImportanceDetectorConstruction::GetCellName(G4int i)
0182 {
0183   std::ostringstream os;
0184   os << "cell_";
0185   if (i < 10) {
0186     os << "0";
0187   }
0188   os << i;
0189   G4String name = os.str();
0190   return name;
0191 }
0192
0193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0194
0195 G4GeometryCell B02ImportanceDetectorConstruction::GetGeometryCell(G4int i)
0196 {
0197   G4String name(GetCellName(i));
0198   const G4VPhysicalVolume* p = 0;
0199   p = fPVolumeStore.GetPVolume(name);
0200   if (p) {
0201     return G4GeometryCell(*p, 0);
0202   }
0203   else {
0204     G4cout << "B02ImportanceDetectorConstruction::GetGeometryCell: " << G4endl
0205            << " couldn't get G4GeometryCell" << G4endl;
0206     return G4GeometryCell(*fGhostWorld, -2);
0207   }
0208 }
0209
0210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0211
0212 G4VPhysicalVolume& B02ImportanceDetectorConstruction::GetWorldVolumeAddress() const
0213 {
0214   return *fGhostWorld;
0215 }
0216
0217 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0218
0219 G4VPhysicalVolume* B02ImportanceDetectorConstruction::GetWorldVolume()
0220 {
0221   return fGhostWorld;
0222 }
0223
0224 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0225
0226 void B02ImportanceDetectorConstruction::SetSensitive()
0227 {
0228   //  -------------------------------------------------
0229   //   The collection names of defined Primitives are
0230   //   0       ConcreteSD/Collisions
0231   //   1       ConcreteSD/CollWeight
0232   //   2       ConcreteSD/Population
0233   //   3       ConcreteSD/TrackEnter
0234   //   4       ConcreteSD/SL
0235   //   5       ConcreteSD/SLW
0236   //   6       ConcreteSD/SLWE
0237   //   7       ConcreteSD/SLW_V
0238   //   8       ConcreteSD/SLWE_V
0239   //  -------------------------------------------------
0240
0241   // moved to ConstructSD() for MT compliance
0242 }
0243
0244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0245 void B02ImportanceDetectorConstruction::ConstructSD()
0246 {
0247   G4SDManager* SDman = G4SDManager::GetSDMpointer();
0248   //
0249   // Sensitive Detector Name
0250   G4String concreteSDname = "ConcreteSD";
0251
0252   //------------------------
0253   // MultiFunctionalDetector
0254   //------------------------
0255   //
0256   // Define MultiFunctionalDetector with name.
0257   G4MultiFunctionalDetector* MFDet = new G4MultiFunctionalDetector(concreteSDname);
0258   SDman->AddNewDetector(MFDet);  // Register SD to SDManager
0259
0260   G4String fltName, particleName;
0261   G4SDParticleFilter* neutronFilter =
0262     new G4SDParticleFilter(fltName = "neutronFilter", particleName = "neutron");
0263
0264   MFDet->SetFilter(neutronFilter);
0265
0266   for (std::vector<G4LogicalVolume*>::iterator it = fLogicalVolumeVector.begin();
0267        it != fLogicalVolumeVector.end(); it++)
0268   {
0269     //      (*it)->SetSensitiveDetector(MFDet);
0270     SetSensitiveDetector((*it)->GetName(), MFDet);
0271   }
0272
0273   G4String psName;
0274   G4PSNofCollision* scorer0 = new G4PSNofCollision(psName = "Collisions");
0275   MFDet->RegisterPrimitive(scorer0);
0276
0277   G4PSNofCollision* scorer1 = new G4PSNofCollision(psName = "CollWeight");
0278   scorer1->Weighted(true);
0279   MFDet->RegisterPrimitive(scorer1);
0280
0281   G4PSPopulation* scorer2 = new G4PSPopulation(psName = "Population");
0282   MFDet->RegisterPrimitive(scorer2);
0283
0284   G4PSTrackCounter* scorer3 = new G4PSTrackCounter(psName = "TrackEnter", fCurrent_In);
0285   MFDet->RegisterPrimitive(scorer3);
0286
0287   G4PSTrackLength* scorer4 = new G4PSTrackLength(psName = "SL");
0288   MFDet->RegisterPrimitive(scorer4);
0289
0290   G4PSTrackLength* scorer5 = new G4PSTrackLength(psName = "SLW");
0291   scorer5->Weighted(true);
0292   MFDet->RegisterPrimitive(scorer5);
0293
0294   G4PSTrackLength* scorer6 = new G4PSTrackLength(psName = "SLWE");
0295   scorer6->Weighted(true);
0296   scorer6->MultiplyKineticEnergy(true);
0297   MFDet->RegisterPrimitive(scorer6);
0298
0299   G4PSTrackLength* scorer7 = new G4PSTrackLength(psName = "SLW_V");
0300   scorer7->Weighted(true);
0301   scorer7->DivideByVelocity(true);
0302   MFDet->RegisterPrimitive(scorer7);
0303
0304   G4PSTrackLength* scorer8 = new G4PSTrackLength(psName = "SLWE_V");
0305   scorer8->Weighted(true);
0306   scorer8->MultiplyKineticEnergy(true);
0307   scorer8->DivideByVelocity(true);
0308   MFDet->RegisterPrimitive(scorer8);
0309 }
0310
0311 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0312 G4VIStore* B02ImportanceDetectorConstruction::CreateImportanceStore()
0313 {
0314   G4cout << " B02ImportanceDetectorConstruction:: Creating Importance Store " << G4endl;
0315   if (!fPVolumeStore.Size()) {
0316     G4Exception("B02ImportanceDetectorConstruction::CreateImportanceStore", "exampleB02_0001",
0317                 RunMustBeAborted, "no physical volumes created yet!");
0318   }
0319
0320   // creating and filling the importance store
0321
0322   //  G4IStore *istore = new G4IStore(*fWorldVolume);
0323
0324   G4IStore* istore = G4IStore::GetInstance(GetName());
0325
0326   G4GeometryCell gWorldVolumeCell(GetWorldVolumeAddress(), 0);
0327
0328   G4double imp = 1;
0329
0330   istore->AddImportanceGeometryCell(1, gWorldVolumeCell);
0331
0332   // set importance values and create scorers
0333   G4int cell(1);
0334   for (cell = 1; cell <= 18; cell++) {
0335     G4GeometryCell gCell = GetGeometryCell(cell);
0336     G4cout << " adding cell: " << cell << " replica: " << gCell.GetReplicaNumber()
0337            << " name: " << gCell.GetPhysicalVolume().GetName() << G4endl;
0338     imp = std::pow(2.0, cell - 1);
0339
0340     G4cout << "Going to assign importance: " << imp
0341            << ", to volume: " << gCell.GetPhysicalVolume().GetName() << G4endl;
0342     // x    aIstore.AddImportanceGeometryCell(imp, gCell);
0343     istore->AddImportanceGeometryCell(imp, gCell.GetPhysicalVolume(), cell);
0344   }
0345
0346   // creating the geometry cell and add both to the store
0347   //  G4GeometryCell gCell = GetGeometryCell(18);
0348
0349   // create importance geometry cell pair for the "rest"cell
0350   // with the same importance as the last concrete cell
0351   G4GeometryCell gCell = GetGeometryCell(19);
0352   //  G4double imp = std::pow(2.0,18);
0353   imp = std::pow(2.0, 17);
0354   istore->AddImportanceGeometryCell(imp, gCell.GetPhysicalVolume(), 19);
0355
0356   return istore;
0357 }
0358
0359 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0360
0361 G4VWeightWindowStore* B02ImportanceDetectorConstruction::CreateWeightWindowStore()
0362 {
0363   G4cout << " B02ImportanceDetectorConstruction:: Creating Importance Store " << G4endl;
0364   if (!fPVolumeStore.Size()) {
0365     G4Exception("B02ImportanceDetectorConstruction::CreateWeightWindowStore", "exampleB02_0002",
0366                 RunMustBeAborted, "no physical volumes created yet!");
0367   }
0368
0369   // creating and filling the importance store
0370
0371   //  G4IStore *istore = new G4IStore(*fWorldVolume);
0372
0373   G4WeightWindowStore* wwstore = G4WeightWindowStore::GetInstance(GetName());
0374
0375   // create one energy region covering the energies of the problem
0376   //
0377   std::set<G4double, std::less<G4double>> enBounds;
0378   enBounds.insert(1 * GeV);
0379   wwstore->SetGeneralUpperEnergyBounds(enBounds);
0380
0381   G4int n = 0;
0382   G4double lowerWeight = 1;
0383   std::vector<G4double> lowerWeights;
0384
0385   lowerWeights.push_back(1);
0386   G4GeometryCell gWorldCell(GetWorldVolumeAddress(), 0);
0387   wwstore->AddLowerWeights(gWorldCell, lowerWeights);
0388
0389   G4int cell(1);
0390   for (cell = 1; cell <= 18; cell++) {
0391     G4GeometryCell gCell = GetGeometryCell(cell);
0392     G4cout << " adding cell: " << cell << " replica: " << gCell.GetReplicaNumber()
0393            << " name: " << gCell.GetPhysicalVolume().GetName() << G4endl;
0394
0395     lowerWeight = 1. / std::pow(2., n++);
0396     G4cout << "Going to assign lower weight: " << lowerWeight
0397            << ", to volume: " << gCell.GetPhysicalVolume().GetName() << G4endl;
0398     lowerWeights.clear();
0399     lowerWeights.push_back(lowerWeight);
0400     wwstore->AddLowerWeights(gCell, lowerWeights);
0401   }
0402
0403   // the remaining part pf the geometry (rest) gets the same
0404   // lower weight bound  as the last conrete cell
0405   //
0406
0407   // create importance geometry cell pair for the "rest"cell
0408   // with the same importance as the last concrete cell
0409   G4GeometryCell gCell = GetGeometryCell(19);
0410   wwstore->AddLowerWeights(gCell, lowerWeights);
0411
0412   return wwstore;
0413 }