Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/biasing/B01/src/B01DetectorConstruction.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 /// \file biasing/B01/src/B01DetectorConstruction.cc
 27 /// \brief Implementation of the B01DetectorConstruction class
 28 //
 29 //
 30 //
 31 
 32 #include "B01DetectorConstruction.hh"
 33 
 34 #include "G4Box.hh"
 35 #include "G4Colour.hh"
 36 #include "G4LogicalVolume.hh"
 37 #include "G4Material.hh"
 38 #include "G4PVPlacement.hh"
 39 #include "G4PhysicalConstants.hh"
 40 #include "G4SystemOfUnits.hh"
 41 #include "G4ThreeVector.hh"
 42 #include "G4Tubs.hh"
 43 #include "G4Types.hh"
 44 #include "G4VisAttributes.hh"
 45 #include "globals.hh"
 46 
 47 #include <set>
 48 #include <sstream>
 49 
 50 // For Primitive Scorers
 51 #include "G4MultiFunctionalDetector.hh"
 52 #include "G4PSNofCollision.hh"
 53 #include "G4PSPopulation.hh"
 54 #include "G4PSTrackCounter.hh"
 55 #include "G4PSTrackLength.hh"
 56 #include "G4SDManager.hh"
 57 #include "G4SDParticleFilter.hh"
 58 
 59 // for importance biasing
 60 #include "G4IStore.hh"
 61 
 62 // for weight window technique
 63 #include "G4WeightWindowStore.hh"
 64 
 65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 66 
 67 B01DetectorConstruction::B01DetectorConstruction()
 68   : G4VUserDetectorConstruction(), fLogicalVolumeVector(), fPhysicalVolumeVector()
 69 {
 70   ;
 71 }
 72 
 73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 74 
 75 B01DetectorConstruction::~B01DetectorConstruction()
 76 {
 77   fLogicalVolumeVector.clear();
 78   fPhysicalVolumeVector.clear();
 79 }
 80 
 81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 82 
 83 G4VPhysicalVolume* B01DetectorConstruction::Construct()
 84 {
 85   G4double pos_x;
 86   G4double pos_y;
 87   G4double pos_z;
 88 
 89   G4double density, pressure, temperature;
 90   G4double A;
 91   G4int Z;
 92 
 93   G4String name, symbol;
 94   G4double z;
 95   G4double fractionmass;
 96 
 97   A = 1.01 * g / mole;
 98   G4Element* elH = new G4Element(name = "Hydrogen", symbol = "H", Z = 1, A);
 99 
100   A = 12.01 * g / mole;
101   G4Element* elC = new G4Element(name = "Carbon", symbol = "C", Z = 6, A);
102 
103   A = 16.00 * g / mole;
104   G4Element* elO = new G4Element(name = "Oxygen", symbol = "O", Z = 8, A);
105 
106   A = 22.99 * g / mole;
107   G4Element* elNa = new G4Element(name = "Natrium", symbol = "Na", Z = 11, A);
108 
109   A = 200.59 * g / mole;
110   G4Element* elHg = new G4Element(name = "Hg", symbol = "Hg", Z = 80, A);
111 
112   A = 26.98 * g / mole;
113   G4Element* elAl = new G4Element(name = "Aluminium", symbol = "Al", Z = 13, A);
114 
115   A = 28.09 * g / mole;
116   G4Element* elSi = new G4Element(name = "Silicon", symbol = "Si", Z = 14, A);
117 
118   A = 39.1 * g / mole;
119   G4Element* elK = new G4Element(name = "K", symbol = "K", Z = 19, A);
120 
121   A = 69.72 * g / mole;
122   G4Element* elCa = new G4Element(name = "Calzium", symbol = "Ca", Z = 31, A);
123 
124   A = 55.85 * g / mole;
125   G4Element* elFe = new G4Element(name = "Iron", symbol = "Fe", Z = 26, A);
126 
127   density = universe_mean_density;  // from PhysicalConstants.h
128   pressure = 3.e-18 * pascal;
129   temperature = 2.73 * kelvin;
130   G4Material* Galactic = new G4Material(name = "Galactic", z = 1., A = 1.01 * g / mole, density,
131                                         kStateGas, temperature, pressure);
132 
133   density = 2.03 * g / cm3;
134   G4Material* Concrete = new G4Material("Concrete", density, 10);
135   Concrete->AddElement(elH, fractionmass = 0.01);
136   Concrete->AddElement(elO, fractionmass = 0.529);
137   Concrete->AddElement(elNa, fractionmass = 0.016);
138   Concrete->AddElement(elHg, fractionmass = 0.002);
139   Concrete->AddElement(elAl, fractionmass = 0.034);
140   Concrete->AddElement(elSi, fractionmass = 0.337);
141   Concrete->AddElement(elK, fractionmass = 0.013);
142   Concrete->AddElement(elCa, fractionmass = 0.044);
143   Concrete->AddElement(elFe, fractionmass = 0.014);
144   Concrete->AddElement(elC, fractionmass = 0.001);
145 
146   /////////////////////////////
147   // world cylinder volume
148   ////////////////////////////
149 
150   // world solid
151 
152   G4double innerRadiusCylinder = 0 * cm;
153   G4double outerRadiusCylinder = 100 * cm;
154   G4double heightCylinder = 100 * cm;
155   G4double startAngleCylinder = 0 * deg;
156   G4double spanningAngleCylinder = 360 * deg;
157 
158   G4Tubs* worldCylinder = new G4Tubs("worldCylinder", innerRadiusCylinder, outerRadiusCylinder,
159                                      heightCylinder, startAngleCylinder, spanningAngleCylinder);
160 
161   // logical world
162 
163   G4LogicalVolume* worldCylinder_log =
164     new G4LogicalVolume(worldCylinder, Galactic, "worldCylinder_log");
165   fLogicalVolumeVector.push_back(worldCylinder_log);
166 
167   name = "shieldWorld";
168   fWorldVolume = new G4PVPlacement(0, G4ThreeVector(0, 0, 0), worldCylinder_log, name, 0, false, 0);
169 
170   fPhysicalVolumeVector.push_back(fWorldVolume);
171 
172   // creating 18 slabs of 10 cm thick concrete
173 
174   G4double innerRadiusShield = 0 * cm;
175   G4double outerRadiusShield = 100 * cm;
176   G4double heightShield = 5 * cm;
177   G4double startAngleShield = 0 * deg;
178   G4double spanningAngleShield = 360 * deg;
179 
180   G4Tubs* aShield = new G4Tubs("aShield", innerRadiusShield, outerRadiusShield, heightShield,
181                                startAngleShield, spanningAngleShield);
182 
183   // logical shield
184 
185   G4LogicalVolume* aShield_log = new G4LogicalVolume(aShield, Concrete, "aShield_log");
186   fLogicalVolumeVector.push_back(aShield_log);
187 
188   G4VisAttributes* pShieldVis = new G4VisAttributes(G4Colour(0.0, 0.0, 1.0));
189   pShieldVis->SetForceSolid(true);
190   aShield_log->SetVisAttributes(pShieldVis);
191 
192   // physical shields
193 
194   G4int i;
195   G4double startz = -85 * cm;
196   for (i = 1; i <= 18; i++) {
197     name = GetCellName(i);
198     pos_x = 0 * cm;
199     pos_y = 0 * cm;
200     pos_z = startz + (i - 1) * (2 * heightShield);
201     G4VPhysicalVolume* pvol = new G4PVPlacement(0, G4ThreeVector(pos_x, pos_y, pos_z), aShield_log,
202                                                 name, worldCylinder_log, false, i);
203     fPhysicalVolumeVector.push_back(pvol);
204   }
205 
206   // filling the rest of the world volume behind the concrete with
207   // another slab which should get the same importance value
208   // or lower weight bound as the last slab
209   //
210   innerRadiusShield = 0 * cm;
211   outerRadiusShield = 100 * cm;
212   heightShield = 5 * cm;
213   startAngleShield = 0 * deg;
214   spanningAngleShield = 360 * deg;
215 
216   G4Tubs* aRest = new G4Tubs("Rest", innerRadiusShield, outerRadiusShield, heightShield,
217                              startAngleShield, spanningAngleShield);
218 
219   G4LogicalVolume* aRest_log = new G4LogicalVolume(aRest, Galactic, "aRest_log");
220   fLogicalVolumeVector.push_back(aRest_log);
221   name = "rest";
222 
223   pos_x = 0 * cm;
224   pos_y = 0 * cm;
225   pos_z = 95 * cm;
226   G4VPhysicalVolume* pvol_rest = new G4PVPlacement(0, G4ThreeVector(pos_x, pos_y, pos_z), aRest_log,
227                                                    name, worldCylinder_log, false,
228                                                    19);  // i=19
229 
230   fPhysicalVolumeVector.push_back(pvol_rest);
231 
232   SetSensitive();
233   return fWorldVolume;
234 }
235 
236 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
237 
238 G4VIStore* B01DetectorConstruction::CreateImportanceStore()
239 {
240   G4cout << " B01DetectorConstruction:: Creating Importance Store " << G4endl;
241   if (!fPhysicalVolumeVector.size()) {
242     G4Exception("B01DetectorConstruction::CreateImportanceStore", "exampleB01_0001",
243                 RunMustBeAborted, "no physical volumes created yet!");
244   }
245 
246   fWorldVolume = fPhysicalVolumeVector[0];
247 
248   // creating and filling the importance store
249 
250   G4IStore* istore = G4IStore::GetInstance();
251 
252   G4int n = 0;
253   G4double imp = 1;
254   istore->AddImportanceGeometryCell(1, *fWorldVolume);
255   for (std::vector<G4VPhysicalVolume*>::iterator it = fPhysicalVolumeVector.begin();
256        it != fPhysicalVolumeVector.end() - 1; it++)
257   {
258     if (*it != fWorldVolume) {
259       imp = std::pow(2., n++);
260       G4cout << "Going to assign importance: " << imp << ", to volume: " << (*it)->GetName()
261              << G4endl;
262       istore->AddImportanceGeometryCell(imp, *(*it), n);
263     }
264   }
265 
266   // the remaining part pf the geometry (rest) gets the same
267   // importance as the last conrete cell
268   //
269   istore->AddImportanceGeometryCell(imp, *(fPhysicalVolumeVector[fPhysicalVolumeVector.size() - 1]),
270                                     ++n);
271 
272   return istore;
273 }
274 
275 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
276 
277 G4VWeightWindowStore* B01DetectorConstruction::CreateWeightWindowStore()
278 {
279   if (!fPhysicalVolumeVector.size()) {
280     G4Exception("B01DetectorConstruction::CreateWeightWindowStore", "exampleB01_0002",
281                 RunMustBeAborted, "no physical volumes created yet!");
282   }
283 
284   fWorldVolume = fPhysicalVolumeVector[0];
285 
286   // creating and filling the weight window store
287 
288   G4WeightWindowStore* wwstore = G4WeightWindowStore::GetInstance();
289 
290   // create one energy region covering the energies of the problem
291   //
292   std::set<G4double, std::less<G4double>> enBounds;
293   enBounds.insert(1 * GeV);
294   wwstore->SetGeneralUpperEnergyBounds(enBounds);
295 
296   G4int n = 0;
297   G4double lowerWeight = 1;
298   std::vector<G4double> lowerWeights;
299 
300   lowerWeights.push_back(1);
301   G4GeometryCell gWorldCell(*fWorldVolume, 0);
302   wwstore->AddLowerWeights(gWorldCell, lowerWeights);
303 
304   for (std::vector<G4VPhysicalVolume*>::iterator it = fPhysicalVolumeVector.begin();
305        it != fPhysicalVolumeVector.end() - 1; it++)
306   {
307     if (*it != fWorldVolume) {
308       lowerWeight = 1. / std::pow(2., n++);
309       G4cout << "Going to assign lower weight: " << lowerWeight
310              << ", to volume: " << (*it)->GetName() << G4endl;
311       G4GeometryCell gCell(*(*it), n);
312       lowerWeights.clear();
313       lowerWeights.push_back(lowerWeight);
314       wwstore->AddLowerWeights(gCell, lowerWeights);
315     }
316   }
317 
318   // the remaining part pf the geometry (rest) gets the same
319   // lower weight bound  as the last conrete cell
320   //
321   G4GeometryCell gRestCell(*(fPhysicalVolumeVector[fPhysicalVolumeVector.size() - 1]), ++n);
322   wwstore->AddLowerWeights(gRestCell, lowerWeights);
323 
324   return wwstore;
325 }
326 
327 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
328 
329 G4String B01DetectorConstruction::GetCellName(G4int i)
330 {
331   std::ostringstream os;
332   os << "cell_";
333   if (i < 10) {
334     os << "0";
335   }
336   os << i;
337   G4String name = os.str();
338   return name;
339 }
340 
341 G4VPhysicalVolume* B01DetectorConstruction::GetWorldVolume()
342 {
343   return fWorldVolume;
344 }
345 
346 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
347 
348 void B01DetectorConstruction::SetSensitive()
349 {
350   //  -------------------------------------------------
351   //   The collection names of defined Primitives are
352   //   0       ConcreteSD/Collisions
353   //   1       ConcreteSD/CollWeight
354   //   2       ConcreteSD/Population
355   //   3       ConcreteSD/TrackEnter
356   //   4       ConcreteSD/SL
357   //   5       ConcreteSD/SLW
358   //   6       ConcreteSD/SLWE
359   //   7       ConcreteSD/SLW_V
360   //   8       ConcreteSD/SLWE_V
361   //  -------------------------------------------------
362 
363   // moved to ConstructSDandField() for MT compliance
364 }
365 
366 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
367 
368 void B01DetectorConstruction::ConstructSDandField()
369 {
370   //  Sensitive Detector Manager.
371   G4SDManager* SDman = G4SDManager::GetSDMpointer();
372   // Sensitive Detector Name
373   G4String concreteSDname = "ConcreteSD";
374 
375   //------------------------
376   // MultiFunctionalDetector
377   //------------------------
378   //
379   // Define MultiFunctionalDetector with name.
380   G4MultiFunctionalDetector* MFDet = new G4MultiFunctionalDetector(concreteSDname);
381   SDman->AddNewDetector(MFDet);  // Register SD to SDManager
382 
383   G4String fltName, particleName;
384   G4SDParticleFilter* neutronFilter =
385     new G4SDParticleFilter(fltName = "neutronFilter", particleName = "neutron");
386 
387   MFDet->SetFilter(neutronFilter);
388 
389   for (std::vector<G4LogicalVolume*>::iterator it = fLogicalVolumeVector.begin();
390        it != fLogicalVolumeVector.end(); it++)
391   {
392     //      (*it)->SetSensitiveDetector(MFDet);
393     SetSensitiveDetector((*it)->GetName(), MFDet);
394   }
395 
396   G4String psName;
397   G4PSNofCollision* scorer0 = new G4PSNofCollision(psName = "Collisions");
398   MFDet->RegisterPrimitive(scorer0);
399 
400   G4PSNofCollision* scorer1 = new G4PSNofCollision(psName = "CollWeight");
401   scorer1->Weighted(true);
402   MFDet->RegisterPrimitive(scorer1);
403 
404   G4PSPopulation* scorer2 = new G4PSPopulation(psName = "Population");
405   MFDet->RegisterPrimitive(scorer2);
406 
407   G4PSTrackCounter* scorer3 = new G4PSTrackCounter(psName = "TrackEnter", fCurrent_In);
408   MFDet->RegisterPrimitive(scorer3);
409 
410   G4PSTrackLength* scorer4 = new G4PSTrackLength(psName = "SL");
411   MFDet->RegisterPrimitive(scorer4);
412 
413   G4PSTrackLength* scorer5 = new G4PSTrackLength(psName = "SLW");
414   scorer5->Weighted(true);
415   MFDet->RegisterPrimitive(scorer5);
416 
417   G4PSTrackLength* scorer6 = new G4PSTrackLength(psName = "SLWE");
418   scorer6->Weighted(true);
419   scorer6->MultiplyKineticEnergy(true);
420   MFDet->RegisterPrimitive(scorer6);
421 
422   G4PSTrackLength* scorer7 = new G4PSTrackLength(psName = "SLW_V");
423   scorer7->Weighted(true);
424   scorer7->DivideByVelocity(true);
425   MFDet->RegisterPrimitive(scorer7);
426 
427   G4PSTrackLength* scorer8 = new G4PSTrackLength(psName = "SLWE_V");
428   scorer8->Weighted(true);
429   scorer8->MultiplyKineticEnergy(true);
430   scorer8->DivideByVelocity(true);
431   MFDet->RegisterPrimitive(scorer8);
432 }
433 
434 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
435