Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/hadronic/Hadr05/src/DetectorConstruction.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 DetectorConstruction.cc
 27 /// \brief Implementation of the DetectorConstruction class
 28 //
 29 //
 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 32 
 33 #include "DetectorConstruction.hh"
 34 
 35 #include "DetectorMessenger.hh"
 36 
 37 #include "G4Box.hh"
 38 #include "G4GeometryManager.hh"
 39 #include "G4LogicalVolume.hh"
 40 #include "G4LogicalVolumeStore.hh"
 41 #include "G4Material.hh"
 42 #include "G4NistManager.hh"
 43 #include "G4PVPlacement.hh"
 44 #include "G4PVReplica.hh"
 45 #include "G4PhysicalConstants.hh"
 46 #include "G4PhysicalVolumeStore.hh"
 47 #include "G4RunManager.hh"
 48 #include "G4SolidStore.hh"
 49 #include "G4SystemOfUnits.hh"
 50 #include "G4UnitsTable.hh"
 51 
 52 #include <iomanip>
 53 
 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 55 
 56 DetectorConstruction::DetectorConstruction()
 57 {
 58   for (G4int i = 0; i < kMaxAbsor; ++i) {
 59     fAbsorMaterial[i] = nullptr;
 60     fAbsorThickness[i] = 0.0;
 61     fSolidAbsor[i] = nullptr;
 62     fLogicAbsor[i] = nullptr;
 63     fPhysiAbsor[i] = nullptr;
 64   }
 65 
 66   // default parameter values of the calorimeter
 67   fNbOfAbsor = 2;
 68   fAbsorThickness[1] = 36 * mm;
 69   fAbsorThickness[2] = 4 * mm;
 70   fNbOfLayers = 50;
 71   fCalorSizeYZ = 1.5 * m;
 72   ComputeCalorParameters();
 73 
 74   // materials
 75   DefineMaterials();
 76   SetWorldMaterial("Galactic");
 77   SetAbsorMaterial(1, "Iron");
 78   SetAbsorMaterial(2, "Scintillator");
 79 
 80   // create commands for interactive definition of the calorimeter
 81   fDetectorMessenger = new DetectorMessenger(this);
 82 }
 83 
 84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 85 
 86 DetectorConstruction::~DetectorConstruction()
 87 {
 88   delete fDetectorMessenger;
 89 }
 90 
 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 92 
 93 void DetectorConstruction::DefineMaterials()
 94 {
 95   // This function illustrates the possible ways to define materials using
 96   // G4 database on G4Elements
 97   G4NistManager* manager = G4NistManager::Instance();
 98   manager->SetVerbose(0);
 99   //
100   // define Elements
101   //
102   G4Element* H = manager->FindOrBuildElement(1);
103   G4Element* C = manager->FindOrBuildElement(6);
104   G4Element* O = manager->FindOrBuildElement(8);
105   //
106   // define an Element from isotopes, by relative abundance
107   //
108   G4int iz, n;  // iz=number of protons  in an isotope;
109                 //  n=number of nucleons in an isotope;
110   G4int ncomponents;
111   G4double z, a;
112   G4double abundance;
113 
114   G4Isotope* U5 = new G4Isotope("U235", iz = 92, n = 235, a = 235.01 * g / mole);
115   G4Isotope* U8 = new G4Isotope("U238", iz = 92, n = 238, a = 238.03 * g / mole);
116 
117   G4Element* U = new G4Element("enriched Uranium", "U", ncomponents = 2);
118   U->AddIsotope(U5, abundance = 90. * perCent);
119   U->AddIsotope(U8, abundance = 10. * perCent);
120 
121   //
122   // define simple materials
123   //
124   G4double density;
125 
126   new G4Material("liquidH2", z = 1., a = 1.008 * g / mole, density = 70.8 * mg / cm3);
127   new G4Material("Aluminium", z = 13., a = 26.98 * g / mole, density = 2.700 * g / cm3);
128   new G4Material("liquidArgon", z = 18, a = 39.948 * g / mole, density = 1.396 * g / cm3);
129   new G4Material("Titanium", z = 22., a = 47.867 * g / mole, density = 4.54 * g / cm3);
130   new G4Material("Iron", z = 26., a = 55.85 * g / mole, density = 7.870 * g / cm3);
131   new G4Material("Copper", z = 29., a = 63.55 * g / mole, density = 8.960 * g / cm3);
132   new G4Material("Tungsten", z = 74., a = 183.85 * g / mole, density = 19.30 * g / cm3);
133   new G4Material("Gold", z = 79., a = 196.97 * g / mole, density = 19.32 * g / cm3);
134   new G4Material("Lead", z = 82., a = 207.20 * g / mole, density = 11.35 * g / cm3);
135   new G4Material("Uranium", z = 92., a = 238.03 * g / mole, density = 18.95 * g / cm3);
136 
137   //
138   // define a material from elements.   case 1: chemical molecule
139   //
140   G4int natoms;
141 
142   G4Material* H2O = new G4Material("Water", density = 1.000 * g / cm3, ncomponents = 2);
143   H2O->AddElement(H, natoms = 2);
144   H2O->AddElement(O, natoms = 1);
145   H2O->GetIonisation()->SetMeanExcitationEnergy(78.0 * eV);
146   H2O->SetChemicalFormula("H_2O");
147 
148   G4Material* CH = new G4Material("Polystyrene", density = 1.032 * g / cm3, ncomponents = 2);
149   CH->AddElement(C, natoms = 1);
150   CH->AddElement(H, natoms = 1);
151 
152   G4Material* Sci = new G4Material("Scintillator", density = 1.032 * g / cm3, ncomponents = 2);
153   Sci->AddElement(C, natoms = 9);
154   Sci->AddElement(H, natoms = 10);
155 
156   Sci->GetIonisation()->SetBirksConstant(0.126 * mm / MeV);
157 
158   //
159   // examples of gas in non STP conditions
160   //
161   G4double temperature, pressure;
162 
163   G4Material* CO2 =
164     new G4Material("CarbonicGas", density = 27. * mg / cm3, ncomponents = 2, kStateGas,
165                    temperature = 325. * kelvin, pressure = 50. * atmosphere);
166   CO2->AddElement(C, natoms = 1);
167   CO2->AddElement(O, natoms = 2);
168 
169   new G4Material("ArgonGas", z = 18, a = 39.948 * g / mole, density = 1.782 * mg / cm3, kStateGas,
170                  273.15 * kelvin, 1 * atmosphere);
171   //
172   // example of vacuum
173   //
174   density = universe_mean_density;  // from PhysicalConstants.h
175   pressure = 3.e-18 * pascal;
176   temperature = 2.73 * kelvin;
177   new G4Material("Galactic", z = 1., a = 1.008 * g / mole, density, kStateGas, temperature,
178                  pressure);
179 
180   //  G4cout << *(G4Material::GetMaterialTable()) << G4endl;
181 }
182 
183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
184 
185 G4Material* DetectorConstruction::MaterialWithSingleIsotope(G4String name, G4String symbol,
186                                                             G4double density, G4int Z, G4int A)
187 {
188   // define a material from an isotope
189   //
190   G4int ncomponents;
191   G4double abundance, massfraction;
192 
193   G4Isotope* isotope = new G4Isotope(symbol, Z, A);
194 
195   G4Element* element = new G4Element(name, symbol, ncomponents = 1);
196   element->AddIsotope(isotope, abundance = 100. * perCent);
197 
198   G4Material* material = new G4Material(name, density, ncomponents = 1);
199   material->AddElement(element, massfraction = 100. * perCent);
200 
201   return material;
202 }
203 
204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
205 
206 void DetectorConstruction::ComputeCalorParameters()
207 {
208   // Compute derived parameters of the calorimeter
209   fLayerThickness = 0.;
210   for (G4int iAbs = 1; iAbs <= fNbOfAbsor; iAbs++) {
211     fLayerThickness += fAbsorThickness[iAbs];
212   }
213   fCalorThickness = fNbOfLayers * fLayerThickness;
214   fWorldSizeX = 1.2 * fCalorThickness;
215   fWorldSizeYZ = 1.2 * fCalorSizeYZ;
216 }
217 
218 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
219 
220 G4VPhysicalVolume* DetectorConstruction::Construct()
221 {
222   if (fPhysiWorld) {
223     return fPhysiWorld;
224   }
225   // complete the Calor parameters definition
226   ComputeCalorParameters();
227 
228   //
229   // World
230   //
231   fSolidWorld = new G4Box("World",  // its name
232                           fWorldSizeX / 2, fWorldSizeYZ / 2, fWorldSizeYZ / 2);  // its size
233 
234   fLogicWorld = new G4LogicalVolume(fSolidWorld,  // its solid
235                                     fWorldMaterial,  // its material
236                                     "World");  // its name
237 
238   fPhysiWorld = new G4PVPlacement(0,  // no rotation
239                                   G4ThreeVector(),  // at (0,0,0)
240                                   fLogicWorld,  // its fLogical volume
241                                   "World",  // its name
242                                   0,  // its mother  volume
243                                   false,  // no boolean operation
244                                   0);  // copy number
245   //
246   // Calorimeter
247   //
248 
249   fSolidCalor = new G4Box("Calorimeter", fCalorThickness / 2, fCalorSizeYZ / 2, fCalorSizeYZ / 2);
250 
251   fLogicCalor = new G4LogicalVolume(fSolidCalor, fWorldMaterial, "Calorimeter");
252 
253   fPhysiCalor = new G4PVPlacement(0,  // no rotation
254                                   G4ThreeVector(),  // at (0,0,0)
255                                   fLogicCalor,  // its fLogical volume
256                                   "Calorimeter",  // its name
257                                   fLogicWorld,  // its mother  volume
258                                   false,  // no boolean operation
259                                   0);  // copy number
260 
261   //
262   // Layers
263   //
264 
265   fSolidLayer = new G4Box("Layer", fLayerThickness / 2, fCalorSizeYZ / 2, fCalorSizeYZ / 2);
266 
267   fLogicLayer = new G4LogicalVolume(fSolidLayer, fWorldMaterial, "Layer");
268   if (fNbOfLayers > 1) {
269     fPhysiLayer =
270       new G4PVReplica("Layer", fLogicLayer, fLogicCalor, kXAxis, fNbOfLayers, fLayerThickness);
271   }
272   else {
273     fPhysiLayer =
274       new G4PVPlacement(0, G4ThreeVector(), fLogicLayer, "Layer", fLogicCalor, false, 0);
275   }
276   //
277   // Absorbers
278   //
279 
280   G4double xfront = -0.5 * fLayerThickness;
281   for (G4int k = 1; k <= fNbOfAbsor; ++k) {
282     fSolidAbsor[k] = new G4Box("Absorber",  // its name
283                                fAbsorThickness[k] / 2, fCalorSizeYZ / 2, fCalorSizeYZ / 2);
284 
285     fLogicAbsor[k] = new G4LogicalVolume(fSolidAbsor[k],  // its solid
286                                          fAbsorMaterial[k],  // its material
287                                          fAbsorMaterial[k]->GetName());
288 
289     G4double xcenter = xfront + 0.5 * fAbsorThickness[k];
290     xfront += fAbsorThickness[k];
291     fPhysiAbsor[k] = new G4PVPlacement(0, G4ThreeVector(xcenter, 0., 0.), fLogicAbsor[k],
292                                        fAbsorMaterial[k]->GetName(), fLogicLayer, false,
293                                        k);  // copy number
294   }
295 
296   PrintCalorParameters();
297 
298   // always return the fPhysical World
299   //
300   return fPhysiWorld;
301 }
302 
303 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
304 
305 void DetectorConstruction::PrintCalorParameters()
306 {
307   G4int prec = 4, wid = prec + 2;
308   G4int dfprec = G4cout.precision(prec);
309 
310   G4double totLength(0.), totRadl(0.), totNuclear(0.);
311 
312   G4cout << "\n-------------------------------------------------------------"
313          << "\n ---> The calorimeter is " << fNbOfLayers << " layers of:";
314   for (G4int i = 1; i <= fNbOfAbsor; ++i) {
315     G4Material* material = fAbsorMaterial[i];
316     G4double radl = material->GetRadlen();
317     G4double nuclearl = material->GetNuclearInterLength();
318     G4double sumThickness = fNbOfLayers * fAbsorThickness[i];
319     G4double nbRadl = sumThickness / radl;
320     G4double nbNuclearl = sumThickness / nuclearl;
321     totLength += sumThickness;
322     totRadl += nbRadl;
323     totNuclear += nbNuclearl;
324     G4cout << "\n   " << std::setw(12) << fAbsorMaterial[i]->GetName() << ": " << std::setw(wid)
325            << G4BestUnit(fAbsorThickness[i], "Length") << "  --->  sum = " << std::setw(wid)
326            << G4BestUnit(sumThickness, "Length") << " = " << std::setw(wid) << nbRadl << " Radl "
327            << " = " << std::setw(wid) << nbNuclearl << " NuclearInteractionLength ";
328   }
329   G4cout << "\n\n                       total thickness = " << std::setw(wid)
330          << G4BestUnit(totLength, "Length") << " = " << std::setw(wid) << totRadl << " Radl "
331          << " = " << std::setw(wid) << totNuclear << " NuclearInteractionLength " << G4endl;
332 
333   G4cout << "                     transverse sizeYZ = " << std::setw(wid)
334          << G4BestUnit(fCalorSizeYZ, "Length") << G4endl;
335   G4cout << "-------------------------------------------------------------\n";
336 
337   G4cout << "\n" << fWorldMaterial << G4endl;
338   for (G4int j = 1; j <= fNbOfAbsor; ++j) {
339     G4cout << "\n" << fAbsorMaterial[j] << G4endl;
340   }
341   G4cout << "\n-------------------------------------------------------------\n";
342 
343   // restore default format
344   G4cout.precision(dfprec);
345 }
346 
347 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
348 
349 void DetectorConstruction::SetWorldMaterial(const G4String& material)
350 {
351   // search the material by its name
352   G4Material* pttoMaterial = G4NistManager::Instance()->FindOrBuildMaterial(material);
353   if (pttoMaterial) {
354     fWorldMaterial = pttoMaterial;
355     if (fLogicWorld) {
356       fLogicWorld->SetMaterial(fWorldMaterial);
357       fLogicLayer->SetMaterial(fWorldMaterial);
358       G4RunManager::GetRunManager()->PhysicsHasBeenModified();
359     }
360   }
361 }
362 
363 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
364 
365 void DetectorConstruction::SetNbOfLayers(G4int ival)
366 {
367   // set the number of Layers
368   //
369   if (ival < 1) {
370     G4cout << "\n --->warning from SetfNbOfLayers: " << ival
371            << " must be at least 1. Command refused" << G4endl;
372     return;
373   }
374   fNbOfLayers = ival;
375 }
376 
377 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
378 
379 void DetectorConstruction::SetNbOfAbsor(G4int ival)
380 {
381   // set the number of Absorbers
382   //
383   if (ival < 1 || ival > (kMaxAbsor - 1)) {
384     G4cout << "\n ---> warning from SetfNbOfAbsor: " << ival << " must be at least 1 and and most "
385            << kMaxAbsor - 1 << ". Command refused" << G4endl;
386     return;
387   }
388   fNbOfAbsor = ival;
389 }
390 
391 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
392 
393 void DetectorConstruction::SetAbsorMaterial(G4int ival, const G4String& material)
394 {
395   // search the material by its name
396   //
397   if (ival > fNbOfAbsor || ival <= 0) {
398     G4cout << "\n --->warning from SetAbsorMaterial: absor number " << ival
399            << " out of range. Command refused" << G4endl;
400     return;
401   }
402 
403   G4Material* pttoMaterial = G4NistManager::Instance()->FindOrBuildMaterial(material);
404   if (pttoMaterial) {
405     fAbsorMaterial[ival] = pttoMaterial;
406     if (fLogicAbsor[ival]) {
407       fLogicAbsor[ival]->SetMaterial(pttoMaterial);
408       G4RunManager::GetRunManager()->PhysicsHasBeenModified();
409     }
410   }
411 }
412 
413 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
414 
415 void DetectorConstruction::SetAbsorThickness(G4int ival, G4double val)
416 {
417   // change Absorber thickness
418   //
419   if (ival > fNbOfAbsor || ival <= 0) {
420     G4cout << "\n --->warning from SetAbsorThickness: absor number " << ival
421            << " out of range. Command refused" << G4endl;
422     return;
423   }
424   if (val <= DBL_MIN) {
425     G4cout << "\n --->warning from SetAbsorThickness: thickness " << val
426            << " out of range. Command refused" << G4endl;
427     return;
428   }
429   fAbsorThickness[ival] = val;
430 }
431 
432 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
433 
434 void DetectorConstruction::SetCalorSizeYZ(G4double val)
435 {
436   // change the transverse size
437   //
438   if (val <= DBL_MIN) {
439     G4cout << "\n --->warning from SetfCalorSizeYZ: thickness " << val
440            << " out of range. Command refused" << G4endl;
441     return;
442   }
443   fCalorSizeYZ = val;
444 }
445 
446 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
447 
448 #include "G4AutoDelete.hh"
449 #include "G4GlobalMagFieldMessenger.hh"
450 
451 void DetectorConstruction::ConstructSDandField()
452 {
453   if (fFieldMessenger.Get() == nullptr) {
454     // Create global magnetic field messenger.
455     // Uniform magnetic field is then created automatically if
456     // the field value is not zero.
457     G4ThreeVector fieldValue = G4ThreeVector();
458     G4GlobalMagFieldMessenger* msg = new G4GlobalMagFieldMessenger(fieldValue);
459     // msg->SetVerboseLevel(1);
460     G4AutoDelete::Register(msg);
461     fFieldMessenger.Put(msg);
462   }
463 }
464 
465 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
466