source: trunk/examples/extended/biasing/B01/src/B01DetectorConstruction.cc

Last change on this file was 1337, checked in by garnier, 14 years ago

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