Geant4_10
B01DetectorConstruction.cc
Go to the documentation of this file.
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 //
28 //
29 //
30 // $Id: B01DetectorConstruction.cc 77475 2013-11-25 09:38:51Z gcosmo $
31 //
32 
33 #include "G4Types.hh"
34 #include <sstream>
35 #include <set>
36 #include "globals.hh"
37 
39 
40 #include "G4Material.hh"
41 #include "G4Box.hh"
42 #include "G4Tubs.hh"
43 #include "G4LogicalVolume.hh"
44 #include "G4ThreeVector.hh"
45 #include "G4PVPlacement.hh"
46 #include "G4VisAttributes.hh"
47 #include "G4Colour.hh"
48 #include "G4PhysicalConstants.hh"
49 #include "G4SystemOfUnits.hh"
50 
51 // For Primitive Scorers
52 #include "G4SDManager.hh"
54 #include "G4SDParticleFilter.hh"
55 #include "G4PSNofCollision.hh"
56 #include "G4PSPopulation.hh"
57 #include "G4PSTrackCounter.hh"
58 #include "G4PSTrackLength.hh"
59 
60 // for importance biasing
61 #include "G4IStore.hh"
62 
63 // for weight window technique
64 #include "G4WeightWindowStore.hh"
65 
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
67 
70  fLogicalVolumeVector(),fPhysicalVolumeVector()
71 {;}
72 
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
74 
76 {
77  fLogicalVolumeVector.clear();
78  fPhysicalVolumeVector.clear();
79 }
80 
81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82 
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 
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 =
131  new G4Material(name="Galactic", z=1., A=1.01*g/mole, density,
132  kStateGas,temperature,pressure);
133 
134  density = 2.03*g/cm3;
135  G4Material* Concrete = new G4Material("Concrete", density, 10);
136  Concrete->AddElement(elH , fractionmass= 0.01);
137  Concrete->AddElement(elO , fractionmass= 0.529);
138  Concrete->AddElement(elNa , fractionmass= 0.016);
139  Concrete->AddElement(elHg , fractionmass= 0.002);
140  Concrete->AddElement(elAl , fractionmass= 0.034);
141  Concrete->AddElement(elSi , fractionmass= 0.337);
142  Concrete->AddElement(elK , fractionmass= 0.013);
143  Concrete->AddElement(elCa , fractionmass= 0.044);
144  Concrete->AddElement(elFe , fractionmass= 0.014);
145  Concrete->AddElement(elC , fractionmass= 0.001);
146 
148  // world cylinder volume
150 
151  // world solid
152 
153  G4double innerRadiusCylinder = 0*cm;
154  G4double outerRadiusCylinder = 100*cm;
155  G4double heightCylinder = 100*cm;
156  G4double startAngleCylinder = 0*deg;
157  G4double spanningAngleCylinder = 360*deg;
158 
159  G4Tubs *worldCylinder = new G4Tubs("worldCylinder",
160  innerRadiusCylinder,
161  outerRadiusCylinder,
162  heightCylinder,
163  startAngleCylinder,
164  spanningAngleCylinder);
165 
166  // logical world
167 
168  G4LogicalVolume *worldCylinder_log =
169  new G4LogicalVolume(worldCylinder, Galactic, "worldCylinder_log");
170  fLogicalVolumeVector.push_back(worldCylinder_log);
171 
172  name = "shieldWorld";
173  fWorldVolume = new
174  G4PVPlacement(0, G4ThreeVector(0,0,0), worldCylinder_log,
175  name, 0, false, 0);
176 
177 
178  fPhysicalVolumeVector.push_back(fWorldVolume);
179 
180  // creating 18 slabs of 10 cm thick concrete
181 
182  G4double innerRadiusShield = 0*cm;
183  G4double outerRadiusShield = 100*cm;
184  G4double heightShield = 5*cm;
185  G4double startAngleShield = 0*deg;
186  G4double spanningAngleShield = 360*deg;
187 
188  G4Tubs *aShield = new G4Tubs("aShield",
189  innerRadiusShield,
190  outerRadiusShield,
191  heightShield,
192  startAngleShield,
193  spanningAngleShield);
194 
195  // logical shield
196 
197  G4LogicalVolume *aShield_log =
198  new G4LogicalVolume(aShield, Concrete, "aShield_log");
199  fLogicalVolumeVector.push_back(aShield_log);
200 
201  G4VisAttributes* pShieldVis = new G4VisAttributes(G4Colour(0.0,0.0,1.0));
202  pShieldVis->SetForceSolid(true);
203  aShield_log->SetVisAttributes(pShieldVis);
204 
205  // physical shields
206 
207  G4int i;
208  G4double startz = -85*cm;
209  for (i=1; i<=18; i++)
210  {
211  name = GetCellName(i);
212  pos_x = 0*cm;
213  pos_y = 0*cm;
214  pos_z = startz + (i-1) * (2*heightShield);
215  G4VPhysicalVolume *pvol =
216  new G4PVPlacement(0,
217  G4ThreeVector(pos_x, pos_y, pos_z),
218  aShield_log,
219  name,
220  worldCylinder_log,
221  false,
222  i);
223  fPhysicalVolumeVector.push_back(pvol);
224  }
225 
226  // filling the rest of the world volume behind the concrete with
227  // another slab which should get the same importance value
228  // or lower weight bound as the last slab
229  //
230  innerRadiusShield = 0*cm;
231  outerRadiusShield = 100*cm;
232  heightShield = 5*cm;
233  startAngleShield = 0*deg;
234  spanningAngleShield = 360*deg;
235 
236  G4Tubs *aRest = new G4Tubs("Rest",
237  innerRadiusShield,
238  outerRadiusShield,
239  heightShield,
240  startAngleShield,
241  spanningAngleShield);
242 
243  G4LogicalVolume *aRest_log =
244  new G4LogicalVolume(aRest, Galactic, "aRest_log");
245  fLogicalVolumeVector.push_back(aRest_log);
246  name = "rest";
247 
248  pos_x = 0*cm;
249  pos_y = 0*cm;
250  pos_z = 95*cm;
251  G4VPhysicalVolume *pvol_rest =
252  new G4PVPlacement(0,
253  G4ThreeVector(pos_x, pos_y, pos_z),
254  aRest_log,
255  name,
256  worldCylinder_log,
257  false,
258  19); // i=19
259 
260  fPhysicalVolumeVector.push_back(pvol_rest);
261 
262  SetSensitive();
263  return fWorldVolume;
264 }
265 
266 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
267 
269 {
270  G4cout << " B01DetectorConstruction:: Creating Importance Store " << G4endl;
271  if (!fPhysicalVolumeVector.size())
272  {
273  G4Exception("B01DetectorConstruction::CreateImportanceStore"
274  ,"exampleB01_0001",RunMustBeAborted
275  ,"no physical volumes created yet!");
276  }
277 
278  fWorldVolume = fPhysicalVolumeVector[0];
279 
280  // creating and filling the importance store
281 
282  G4IStore *istore = G4IStore::GetInstance();
283 
284  G4int n = 0;
285  G4double imp =1;
286  istore->AddImportanceGeometryCell(1, *fWorldVolume);
287  for (std::vector<G4VPhysicalVolume *>::iterator
288  it = fPhysicalVolumeVector.begin();
289  it != fPhysicalVolumeVector.end() - 1; it++)
290  {
291  if (*it != fWorldVolume)
292  {
293  imp = std::pow(2., n++);
294  G4cout << "Going to assign importance: " << imp << ", to volume: "
295  << (*it)->GetName() << G4endl;
296  istore->AddImportanceGeometryCell(imp, *(*it),n);
297  }
298  }
299 
300  // the remaining part pf the geometry (rest) gets the same
301  // importance as the last conrete cell
302  //
303  istore->AddImportanceGeometryCell(imp,
304  *(fPhysicalVolumeVector[fPhysicalVolumeVector.size()-1]),++n);
305 
306  return istore;
307 }
308 
309 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
310 
312 {
313  if (!fPhysicalVolumeVector.size())
314  {
315  G4Exception("B01DetectorConstruction::CreateWeightWindowStore"
316  ,"exampleB01_0002",RunMustBeAborted
317  ,"no physical volumes created yet!");
318  }
319 
320  fWorldVolume = fPhysicalVolumeVector[0];
321 
322  // creating and filling the weight window store
323 
325 
326  // create one energy region covering the energies of the problem
327  //
328  std::set<G4double, std::less<G4double> > enBounds;
329  enBounds.insert(1 * GeV);
330  wwstore->SetGeneralUpperEnergyBounds(enBounds);
331 
332  G4int n = 0;
333  G4double lowerWeight =1;
334  std::vector<G4double> lowerWeights;
335 
336  lowerWeights.push_back(1);
337  G4GeometryCell gWorldCell(*fWorldVolume,0);
338  wwstore->AddLowerWeights(gWorldCell, lowerWeights);
339 
340  for (std::vector<G4VPhysicalVolume *>::iterator
341  it = fPhysicalVolumeVector.begin();
342  it != fPhysicalVolumeVector.end() - 1; it++)
343  {
344  if (*it != fWorldVolume)
345  {
346  lowerWeight = 1./std::pow(2., n++);
347  G4cout << "Going to assign lower weight: " << lowerWeight
348  << ", to volume: "
349  << (*it)->GetName() << G4endl;
350  G4GeometryCell gCell(*(*it),n);
351  lowerWeights.clear();
352  lowerWeights.push_back(lowerWeight);
353  wwstore->AddLowerWeights(gCell, lowerWeights);
354  }
355  }
356 
357  // the remaining part pf the geometry (rest) gets the same
358  // lower weight bound as the last conrete cell
359  //
361  gRestCell(*(fPhysicalVolumeVector[fPhysicalVolumeVector.size()-1]), ++n);
362  wwstore->AddLowerWeights(gRestCell, lowerWeights);
363 
364  return wwstore;
365 }
366 
367 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
368 
370 {
371  std::ostringstream os;
372  os << "cell_";
373  if (i<10)
374  {
375  os << "0";
376  }
377  os << i ;
378  G4String name = os.str();
379  return name;
380 }
381 
383  return fWorldVolume;
384 }
385 
386 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
387 
389 
390  // -------------------------------------------------
391  // The collection names of defined Primitives are
392  // 0 ConcreteSD/Collisions
393  // 1 ConcreteSD/CollWeight
394  // 2 ConcreteSD/Population
395  // 3 ConcreteSD/TrackEnter
396  // 4 ConcreteSD/SL
397  // 5 ConcreteSD/SLW
398  // 6 ConcreteSD/SLWE
399  // 7 ConcreteSD/SLW_V
400  // 8 ConcreteSD/SLWE_V
401  // -------------------------------------------------
402 
403  // moved to ConstructSDandField() for MT compliance
404 
405 }
406 
407 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
408 
410 {
411 
412  // Sensitive Detector Manager.
414  // Sensitive Detector Name
415  G4String concreteSDname = "ConcreteSD";
416 
417  //------------------------
418  // MultiFunctionalDetector
419  //------------------------
420  //
421  // Define MultiFunctionalDetector with name.
422  G4MultiFunctionalDetector* MFDet =
423  new G4MultiFunctionalDetector(concreteSDname);
424  SDman->AddNewDetector( MFDet ); // Register SD to SDManager
425 
426 
427  G4String fltName,particleName;
428  G4SDParticleFilter* neutronFilter =
429  new G4SDParticleFilter(fltName="neutronFilter", particleName="neutron");
430 
431  MFDet->SetFilter(neutronFilter);
432 
433 
434  for (std::vector<G4LogicalVolume *>::iterator it =
435  fLogicalVolumeVector.begin();
436  it != fLogicalVolumeVector.end(); it++){
437  // (*it)->SetSensitiveDetector(MFDet);
438  SetSensitiveDetector((*it)->GetName(), MFDet);
439  }
440 
441  G4String psName;
442  G4PSNofCollision* scorer0 = new G4PSNofCollision(psName="Collisions");
443  MFDet->RegisterPrimitive(scorer0);
444 
445 
446  G4PSNofCollision* scorer1 = new G4PSNofCollision(psName="CollWeight");
447  scorer1->Weighted(true);
448  MFDet->RegisterPrimitive(scorer1);
449 
450 
451  G4PSPopulation* scorer2 = new G4PSPopulation(psName="Population");
452  MFDet->RegisterPrimitive(scorer2);
453 
454  G4PSTrackCounter* scorer3 = new G4PSTrackCounter(psName="TrackEnter"
455  ,fCurrent_In);
456  MFDet->RegisterPrimitive(scorer3);
457 
458  G4PSTrackLength* scorer4 = new G4PSTrackLength(psName="SL");
459  MFDet->RegisterPrimitive(scorer4);
460 
461  G4PSTrackLength* scorer5 = new G4PSTrackLength(psName="SLW");
462  scorer5->Weighted(true);
463  MFDet->RegisterPrimitive(scorer5);
464 
465  G4PSTrackLength* scorer6 = new G4PSTrackLength(psName="SLWE");
466  scorer6->Weighted(true);
467  scorer6->MultiplyKineticEnergy(true);
468  MFDet->RegisterPrimitive(scorer6);
469 
470  G4PSTrackLength* scorer7 = new G4PSTrackLength(psName="SLW_V");
471  scorer7->Weighted(true);
472  scorer7->DivideByVelocity(true);
473  MFDet->RegisterPrimitive(scorer7);
474 
475  G4PSTrackLength* scorer8 = new G4PSTrackLength(psName="SLWE_V");
476  scorer8->Weighted(true);
477  scorer8->MultiplyKineticEnergy(true);
478  scorer8->DivideByVelocity(true);
479  MFDet->RegisterPrimitive(scorer8);
480 
481 }
482 
483 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4bool RegisterPrimitive(G4VPrimitiveScorer *)
G4VPhysicalVolume * GetWorldVolume()
G4String symbol
Definition: TRTMaterials.hh:40
CLHEP::Hep3Vector G4ThreeVector
int universe_mean_density
Definition: hepunit.py:307
void AddLowerWeights(const G4GeometryCell &gCell, const std::vector< G4double > &lowerWeights)
virtual G4VPhysicalVolume * Construct()
Definition: G4Tubs.hh:84
G4Element * elC
Definition: TRTMaterials.hh:48
const XML_Char * name
Definition: expat.h:151
void SetForceSolid(G4bool)
Definition of the B01DetectorConstruction class.
int G4int
Definition: G4Types.hh:78
void Weighted(G4bool flg=true)
G4Element * elH
Definition: TRTMaterials.hh:50
G4double density
Definition: TRTMaterials.hh:39
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
Char_t n[5]
void Weighted(G4bool flg=true)
G4GLOB_DLL std::ostream G4cout
G4VWeightWindowStore * CreateWeightWindowStore()
G4Element * elO
Definition: TRTMaterials.hh:46
Float_t Z
Definition: plot.C:39
#define pascal
void AddImportanceGeometryCell(G4double importance, const G4GeometryCell &gCell)
Definition: G4IStore.cc:104
void SetFilter(G4VSDFilter *value)
static G4WeightWindowStore * GetInstance()
void SetSensitiveDetector(const G4String &logVolName, G4VSensitiveDetector *aSD, G4bool multi=false)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
void AddNewDetector(G4VSensitiveDetector *aSD)
Definition: G4SDManager.cc:67
static G4IStore * GetInstance()
Definition: G4IStore.cc:211
static G4SDManager * GetSDMpointer()
Definition: G4SDManager.cc:40
void MultiplyKineticEnergy(G4bool flg=true)
void SetGeneralUpperEnergyBounds(const std::set< G4double, std::less< G4double > > &enBounds)
tuple z
Definition: test.py:28
#define G4endl
Definition: G4ios.hh:61
void AddElement(G4Element *element, G4int nAtoms)
Definition: G4Material.cc:345
double G4double
Definition: G4Types.hh:76
void DivideByVelocity(G4bool flg=true)
void SetVisAttributes(const G4VisAttributes *pVA)