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