ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DetectorConstruction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DetectorConstruction.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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 
33 #include "DetectorConstruction.hh"
34 #include "DetectorMessenger.hh"
35 
36 #include "G4NistManager.hh"
37 #include "G4Material.hh"
38 #include "G4Box.hh"
39 #include "G4LogicalVolume.hh"
40 #include "G4PVPlacement.hh"
41 #include "G4PVReplica.hh"
42 
43 #include "G4GeometryManager.hh"
44 #include "G4PhysicalVolumeStore.hh"
45 #include "G4LogicalVolumeStore.hh"
46 #include "G4SolidStore.hh"
47 
48 #include "G4RunManager.hh"
49 #include "G4SystemOfUnits.hh"
50 #include "G4UnitsTable.hh"
51 #include "G4PhysicalConstants.hh"
52 
53 #include <iomanip>
54 
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 
57 DetectorConstruction::DetectorConstruction()
58  :G4VUserDetectorConstruction(),
59  fWorldMaterial(nullptr),fSolidWorld(nullptr),fLogicWorld(nullptr),
60  fPhysiWorld(nullptr),fSolidCalor(nullptr),fLogicCalor(nullptr),
61  fPhysiCalor(nullptr),fSolidLayer(nullptr),fLogicLayer(nullptr),
62  fPhysiLayer(nullptr)
63 {
64  for(G4int i=0; i<kMaxAbsor; ++i) {
65  fAbsorMaterial[i] = nullptr;
66  fAbsorThickness[i] = 0.0;
67  fSolidAbsor[i] = nullptr;
68  fLogicAbsor[i] = nullptr;
69  fPhysiAbsor[i] = nullptr;
70  }
71 
72  // default parameter values of the calorimeter
73  fNbOfAbsor = 2;
74  fAbsorThickness[1] = 2.3*mm;
75  fAbsorThickness[2] = 5.7*mm;
76  fNbOfLayers = 50;
77  fCalorSizeYZ = 40.*cm;
78  ComputeCalorParameters();
79 
80  // materials
81  DefineMaterials();
82  SetWorldMaterial("Galactic");
83  SetAbsorMaterial(1,"G4_Pb");
84  SetAbsorMaterial(2,"G4_lAr");
85 
86  // create commands for interactive definition of the calorimeter
87  fDetectorMessenger = new DetectorMessenger(this);
88 }
89 
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
91 
92 DetectorConstruction::~DetectorConstruction()
93 {
94  delete fDetectorMessenger;
95 }
96 
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
98 
99 void DetectorConstruction::DefineMaterials()
100 {
101  // This function illustrates the possible ways to define materials using
102  // G4 database on G4Elements
103  G4NistManager* manager = G4NistManager::Instance();
104  manager->SetVerbose(0);
105  //
106  // define Elements
107  //
108  G4double z,a;
109 
110  G4Element* H = manager->FindOrBuildElement(1);
111  G4Element* C = manager->FindOrBuildElement(6);
112  G4Element* N = manager->FindOrBuildElement(7);
113  G4Element* O = manager->FindOrBuildElement(8);
114  G4Element* Si = manager->FindOrBuildElement(14);
115  G4Element* Ge = manager->FindOrBuildElement(32);
116  G4Element* Sb = manager->FindOrBuildElement(51);
117  G4Element* I = manager->FindOrBuildElement(53);
118  G4Element* Cs = manager->FindOrBuildElement(55);
119  G4Element* Pb = manager->FindOrBuildElement(82);
120  G4Element* Bi = manager->FindOrBuildElement(83);
121 
122  //
123  // define an Element from isotopes, by relative abundance
124  //
125  G4int iz, n; //iz=number of protons in an isotope;
126  // n=number of nucleons in an isotope;
127  G4int ncomponents;
128  G4double abundance;
129 
130  G4Isotope* U5 = new G4Isotope("U235", iz=92, n=235, a=235.01*g/mole);
131  G4Isotope* U8 = new G4Isotope("U238", iz=92, n=238, a=238.03*g/mole);
132 
133  G4Element* U = new G4Element("enriched Uranium", "U", ncomponents=2);
134  U->AddIsotope(U5, abundance= 90.*perCent);
135  U->AddIsotope(U8, abundance= 10.*perCent);
136 
137  //
138  // define simple materials
139  //
140  G4double density;
141 
142  new G4Material("liquidH2", z=1., a= 1.008*g/mole, density= 70.8*mg/cm3);
143  new G4Material("Aluminium", z=13., a= 26.98*g/mole, density= 2.700*g/cm3);
144  new G4Material("Titanium", z=22., a= 47.867*g/mole, density= 4.54*g/cm3);
145  new G4Material("Iron", z=26., a= 55.85*g/mole, density= 7.870*g/cm3);
146  new G4Material("Copper", z=29., a= 63.55*g/mole, density= 8.960*g/cm3);
147  new G4Material("Tungsten", z=74., a= 183.85*g/mole, density= 19.30*g/cm3);
148  new G4Material("Gold", z=79., a= 196.97*g/mole, density= 19.32*g/cm3);
149  new G4Material("Uranium", z=92., a= 238.03*g/mole, density= 18.95*g/cm3);
150 
151  //
152  // define a material from elements. case 1: chemical molecule
153  //
154  G4int natoms;
155 
156  G4Material* H2O =
157  new G4Material("Water", density= 1.000*g/cm3, ncomponents=2);
158  H2O->AddElement(H, natoms=2);
159  H2O->AddElement(O, natoms=1);
160  H2O->GetIonisation()->SetMeanExcitationEnergy(78.0*eV);
161  H2O->SetChemicalFormula("H_2O");
162 
163  G4Material* CH =
164  new G4Material("Polystyrene", density= 1.032*g/cm3, ncomponents=2);
165  CH->AddElement(C, natoms=1);
166  CH->AddElement(H, natoms=1);
167 
168  G4Material* Sci =
169  new G4Material("Scintillator", density= 1.032*g/cm3, ncomponents=2);
170  Sci->AddElement(C, natoms=9);
171  Sci->AddElement(H, natoms=10);
172 
173  Sci->GetIonisation()->SetBirksConstant(0.126*mm/MeV);
174 
175  G4Material* Lct =
176  new G4Material("Lucite", density= 1.185*g/cm3, ncomponents=3);
177  Lct->AddElement(C, 59.97*perCent);
178  Lct->AddElement(H, 8.07*perCent);
179  Lct->AddElement(O, 31.96*perCent);
180 
181  G4Material* Sili =
182  new G4Material("Silicon", density= 2.330*g/cm3, ncomponents=1);
183  Sili->AddElement(Si, natoms=1);
184 
185  G4Material* SiO2 =
186  new G4Material("quartz", density= 2.200*g/cm3, ncomponents=2);
187  SiO2->AddElement(Si, natoms=1);
188  SiO2->AddElement(O , natoms=2);
189 
190  G4Material* G10 =
191  new G4Material("NemaG10", density= 1.700*g/cm3, ncomponents=4);
192  G10->AddElement(Si, natoms=1);
193  G10->AddElement(O , natoms=2);
194  G10->AddElement(C , natoms=3);
195  G10->AddElement(H , natoms=3);
196 
197  G4Material* CsI =
198  new G4Material("CsI", density= 4.534*g/cm3, ncomponents=2);
199  CsI->AddElement(Cs, natoms=1);
200  CsI->AddElement(I , natoms=1);
201  CsI->GetIonisation()->SetMeanExcitationEnergy(553.1*eV);
202 
203  G4Material* BGO =
204  new G4Material("BGO", density= 7.10*g/cm3, ncomponents=3);
205  BGO->AddElement(O , natoms=12);
206  BGO->AddElement(Ge, natoms= 3);
207  BGO->AddElement(Bi, natoms= 4);
208 
209  //SiNx
210  density= 3.1 *g/cm3;
211  G4Material* SiNx= new G4Material("SiNx", density, ncomponents=3);
212  SiNx-> AddElement(Si, 300);
213  SiNx-> AddElement(N, 310);
214  SiNx-> AddElement(H, 6);
215 
216  //
217  // define gaseous materials using G4 NIST database
218  //
219  G4double fractionmass;
220 
221  G4Material* Air = manager->FindOrBuildMaterial("G4_AIR");
222  manager->ConstructNewGasMaterial("Air20","G4_AIR",293.*kelvin,1.*atmosphere);
223 
224  G4Material* lAr = manager->FindOrBuildMaterial("G4_lAr");
225  G4Material* lArEm3 = new G4Material("liquidArgon", density= 1.390*g/cm3,
226  ncomponents=1);
227  lArEm3->AddMaterial(lAr, fractionmass=1.0);
228 
229  //
230  // define a material from elements and others materials (mixture of mixtures)
231  //
232 
233  G4Material* Lead = new G4Material("Lead",density=11.35*g/cm3,ncomponents=1);
234  Lead->AddElement(Pb, fractionmass=1.0);
235 
236  G4Material* LeadSb = new G4Material("LeadSb", density=11.35*g/cm3,
237  ncomponents=2);
238  LeadSb->AddElement(Sb, fractionmass=4.*perCent);
239  LeadSb->AddElement(Pb, fractionmass=96.*perCent);
240 
241  G4Material* Aerog = new G4Material("Aerogel", density= 0.200*g/cm3,
242  ncomponents=3);
243  Aerog->AddMaterial(SiO2, fractionmass=62.5*perCent);
244  Aerog->AddMaterial(H2O , fractionmass=37.4*perCent);
245  Aerog->AddElement (C , fractionmass= 0.1*perCent);
246 
247  //
248  // examples of gas in non STP conditions
249  //
250  G4double temperature, pressure;
251 
252  G4Material* CO2 =
253  new G4Material("CarbonicGas", density= 27.*mg/cm3, ncomponents=2,
254  kStateGas, temperature= 325.*kelvin, pressure= 50.*atmosphere);
255  CO2->AddElement(C, natoms=1);
256  CO2->AddElement(O, natoms=2);
257 
258  G4Material* steam =
259  new G4Material("WaterSteam", density= 1.0*mg/cm3, ncomponents=1,
260  kStateGas, temperature= 273*kelvin, pressure= 1*atmosphere);
261  steam->AddMaterial(H2O, fractionmass=1.);
262 
263  new G4Material("ArgonGas", z=18, a=39.948*g/mole, density= 1.782*mg/cm3,
264  kStateGas, 273.15*kelvin, 1*atmosphere);
265  //
266  // examples of vacuum
267  //
268 
269  density = universe_mean_density; //from PhysicalConstants.h
270  pressure = 3.e-18*pascal;
271  temperature = 2.73*kelvin;
272  new G4Material("Galactic", z=1., a=1.008*g/mole, density,
273  kStateGas,temperature,pressure);
274 
275  density = 1.e-5*g/cm3;
276  pressure = 2.e-2*bar;
277  temperature = STP_Temperature; //from PhysicalConstants.h
278  G4Material* beam =
279  new G4Material("Beam", density, ncomponents=1,
280  kStateGas,temperature,pressure);
281  beam->AddMaterial(Air, fractionmass=1.);
282 
283  // G4cout << *(G4Material::GetMaterialTable()) << G4endl;
284 }
285 
286 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
287 
288 void DetectorConstruction::ComputeCalorParameters()
289 {
290  // Compute derived parameters of the calorimeter
291  fLayerThickness = 0.;
292  for (G4int iAbs=1; iAbs<=fNbOfAbsor; iAbs++) {
293  fLayerThickness += fAbsorThickness[iAbs];
294  }
295  fCalorThickness = fNbOfLayers*fLayerThickness;
296  fWorldSizeX = 1.2*fCalorThickness;
297  fWorldSizeYZ = 1.2*fCalorSizeYZ;
298 }
299 
300 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
301 
302 G4VPhysicalVolume* DetectorConstruction::Construct()
303 {
304  if(fPhysiWorld) { return fPhysiWorld; }
305  // complete the Calor parameters definition
306  ComputeCalorParameters();
307 
308  //
309  // World
310  //
311  fSolidWorld = new G4Box("World", //its name
312  fWorldSizeX/2,fWorldSizeYZ/2,fWorldSizeYZ/2); //its size
313 
314  fLogicWorld = new G4LogicalVolume(fSolidWorld, //its solid
315  fWorldMaterial, //its material
316  "World"); //its name
317 
318  fPhysiWorld = new G4PVPlacement(0, //no rotation
319  G4ThreeVector(), //at (0,0,0)
320  fLogicWorld, //its fLogical volume
321  "World", //its name
322  0, //its mother volume
323  false, //no boolean operation
324  0); //copy number
325  //
326  // Calorimeter
327  //
328 
329  fSolidCalor = new G4Box("Calorimeter",
330  fCalorThickness/2,fCalorSizeYZ/2,fCalorSizeYZ/2);
331 
332  fLogicCalor = new G4LogicalVolume(fSolidCalor,
333  fWorldMaterial,
334  "Calorimeter");
335 
336  fPhysiCalor = new G4PVPlacement(0, //no rotation
337  G4ThreeVector(), //at (0,0,0)
338  fLogicCalor, //its fLogical volume
339  "Calorimeter", //its name
340  fLogicWorld, //its mother volume
341  false, //no boolean operation
342  0); //copy number
343 
344  //
345  // Layers
346  //
347 
348  fSolidLayer = new G4Box("Layer",
349  fLayerThickness/2,fCalorSizeYZ/2,fCalorSizeYZ/2);
350 
351  fLogicLayer = new G4LogicalVolume(fSolidLayer,
352  fWorldMaterial,
353  "Layer");
354  if (fNbOfLayers > 1) {
355  fPhysiLayer = new G4PVReplica("Layer",
356  fLogicLayer,
357  fLogicCalor,
358  kXAxis,
359  fNbOfLayers,
360  fLayerThickness);
361  } else {
362  fPhysiLayer = new G4PVPlacement(0,
363  G4ThreeVector(),
364  fLogicLayer,
365  "Layer",
366  fLogicCalor,
367  false,
368  0);
369  }
370  //
371  // Absorbers
372  //
373 
374  G4double xfront = -0.5*fLayerThickness;
375  for (G4int k=1; k<=fNbOfAbsor; ++k) {
376  fSolidAbsor[k] = new G4Box("Absorber", //its name
377  fAbsorThickness[k]/2,fCalorSizeYZ/2,fCalorSizeYZ/2);
378 
379  fLogicAbsor[k] = new G4LogicalVolume(fSolidAbsor[k], //its solid
380  fAbsorMaterial[k], //its material
381  fAbsorMaterial[k]->GetName());
382 
383  G4double xcenter = xfront+0.5*fAbsorThickness[k];
384  xfront += fAbsorThickness[k];
385  fPhysiAbsor[k] = new G4PVPlacement(0,
386  G4ThreeVector(xcenter,0.,0.),
387  fLogicAbsor[k],
388  fAbsorMaterial[k]->GetName(),
389  fLogicLayer,
390  false,
391  k); //copy number
392 
393  }
394 
395  PrintCalorParameters();
396 
397  //always return the fPhysical World
398  //
399  return fPhysiWorld;
400 }
401 
402 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
403 
404 void DetectorConstruction::PrintCalorParameters()
405 {
406  G4cout << "\n-------------------------------------------------------------"
407  << "\n ---> The calorimeter is " << fNbOfLayers << " layers of:";
408  for (G4int i=1; i<=fNbOfAbsor; ++i) {
409  G4cout << "\n \t" << std::setw(12) << fAbsorMaterial[i]->GetName() <<": "
410  << std::setw(6) << G4BestUnit(fAbsorThickness[i],"Length");
411  }
412  G4cout << "\n-------------------------------------------------------------\n";
413 
414  G4cout << "\n" << fWorldMaterial << G4endl;
415  for (G4int j=1; j<=fNbOfAbsor; ++j) {
416  G4cout << "\n" << fAbsorMaterial[j] << G4endl;
417  }
418  G4cout << "\n-------------------------------------------------------------\n";
419 }
420 
421 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
422 
423 void DetectorConstruction::SetWorldMaterial(const G4String& material)
424 {
425  // search the material by its name
426  G4Material* pttoMaterial =
427  G4NistManager::Instance()->FindOrBuildMaterial(material);
428  if(pttoMaterial) {
429  fWorldMaterial = pttoMaterial;
430  if(fLogicWorld) {
431  fLogicWorld->SetMaterial(fWorldMaterial);
432  fLogicLayer->SetMaterial(fWorldMaterial);
433  G4RunManager::GetRunManager()->PhysicsHasBeenModified();
434  }
435  }
436 }
437 
438 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
439 
440 void DetectorConstruction::SetNbOfLayers(G4int ival)
441 {
442  // set the number of Layers
443  //
444  if (ival < 1)
445  { G4cout << "\n --->warning from SetfNbOfLayers: "
446  << ival << " must be at least 1. Command refused" << G4endl;
447  return;
448  }
449  fNbOfLayers = ival;
450 }
451 
452 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
453 
454 void DetectorConstruction::SetNbOfAbsor(G4int ival)
455 {
456  // set the number of Absorbers
457  //
458  if (ival < 1 || ival > (kMaxAbsor-1))
459  { G4cout << "\n ---> warning from SetfNbOfAbsor: "
460  << ival << " must be at least 1 and and most " << kMaxAbsor-1
461  << ". Command refused" << G4endl;
462  return;
463  }
464  fNbOfAbsor = ival;
465 }
466 
467 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
468 
469 void DetectorConstruction::SetAbsorMaterial(G4int ival,
470  const G4String& material)
471 {
472  // search the material by its name
473  //
474  if (ival > fNbOfAbsor || ival <= 0)
475  { G4cout << "\n --->warning from SetAbsorMaterial: absor number "
476  << ival << " out of range. Command refused" << G4endl;
477  return;
478  }
479 
480  G4Material* pttoMaterial =
481  G4NistManager::Instance()->FindOrBuildMaterial(material);
482  if (pttoMaterial) {
483  fAbsorMaterial[ival] = pttoMaterial;
484  if(fLogicAbsor[ival]) {
485  fLogicAbsor[ival]->SetMaterial(pttoMaterial);
486  G4RunManager::GetRunManager()->PhysicsHasBeenModified();
487  }
488  }
489 }
490 
491 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
492 
493 void DetectorConstruction::SetAbsorThickness(G4int ival, G4double val)
494 {
495  // change Absorber thickness
496  //
497  if (ival > fNbOfAbsor || ival <= 0)
498  { G4cout << "\n --->warning from SetAbsorThickness: absor number "
499  << ival << " out of range. Command refused" << G4endl;
500  return;
501  }
502  if (val <= DBL_MIN)
503  { G4cout << "\n --->warning from SetAbsorThickness: thickness "
504  << val << " out of range. Command refused" << G4endl;
505  return;
506  }
507  fAbsorThickness[ival] = val;
508 }
509 
510 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
511 
512 void DetectorConstruction::SetCalorSizeYZ(G4double val)
513 {
514  // change the transverse size
515  //
516  if (val <= DBL_MIN)
517  { G4cout << "\n --->warning from SetfCalorSizeYZ: thickness "
518  << val << " out of range. Command refused" << G4endl;
519  return;
520  }
521  fCalorSizeYZ = val;
522 }
523 
524 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
525 
526 #include "G4GlobalMagFieldMessenger.hh"
527 #include "G4AutoDelete.hh"
528 
529 void DetectorConstruction::ConstructSDandField()
530 {
531  if ( fFieldMessenger.Get() == nullptr ) {
532  // Create global magnetic field messenger.
533  // Uniform magnetic field is then created automatically if
534  // the field value is not zero.
535  G4ThreeVector fieldValue = G4ThreeVector();
536  G4GlobalMagFieldMessenger* msg =
537  new G4GlobalMagFieldMessenger(fieldValue);
538  //msg->SetVerboseLevel(1);
539  G4AutoDelete::Register(msg);
540  fFieldMessenger.Put( msg );
541  }
542 }
543 
544 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......