SoilConstruction.cc
Go to the documentation of this file.
1 #include "SoilConstruction.h"
2 
3 #include <fwk/CentralConfig.h>
4 #include <utl/ErrorLogger.h>
5 #include <utl/Reader.h>
6 #include <utl/MathConstants.h>
7 
8 //G4 Classes
9 #include <G4RunManager.hh>
10 #include <G4RegionStore.hh>
11 
12 #include <G4Box.hh>
13 #include <G4Tubs.hh>
14 #include <G4LogicalVolume.hh>
15 #include <G4ThreeVector.hh>
16 #include <G4PVPlacement.hh>
17 #include <G4Transform3D.hh>
18 #include <G4RotationMatrix.hh>
19 #include <G4VisAttributes.hh>
20 #include <G4Colour.hh>
21 #include <G4Material.hh>
22 #include <G4Element.hh>
23 #include <G4MaterialPropertiesTable.hh>
24 
25 //#warning "MD_THETASIM_MAX and MD_SIM_AREA are brutally hard coded!!!"
26 #define MD_THETASIM_MAX 80*utl::degree
27 #define MD_SIM_AREA 30*utl::m2
28 #define MD_LARGE_SIDE 12*utl::m
29 
30 using namespace fwk;
31 using namespace utl;
32 using namespace GroundPropagatorAG;
33 
34 SoilConstruction::SoilConstruction ( G4double grammage ) :
35 fSoilGrammage( grammage )
36 {
37  fLogicalWorld = 0;
38  fPhysicalWorld = 0;
39  fLogicalSoil = 0;
40  fPhysicalSoil = 0;
41 
42  INFO("\n++++\nGetting parameters from XML file...");
43  CentralConfig* theConfig = CentralConfig::GetInstance();
44  Branch topB = theConfig->GetTopBranch("MdGroundPropagator");
45 
46  //Load X-Y world dimensions
47  std::vector<G4double> vWorldSize;
48  Branch worldB = topB.GetChild("world");
49  worldB.GetChild("size").GetData(vWorldSize);
50 
51  std::vector<G4double> vSoilSampleA;
52  std::vector<G4double> vSoilSampleB;
53  std::vector<G4double> vSoilSampleC;
54 
55  G4double SampleA=0;
56  G4double SampleB=0;
57  G4double SampleC=0;
58 
59  Branch soilB = topB.GetChild("soil");
60 
61  soilB.GetChild("densityA").GetData(vSoilSampleA);
62  for( size_t i = 0; i < vSoilSampleA.size(); ++i )
63  SampleA += vSoilSampleA[i];
64  SampleA /= vSoilSampleA.size();
65 
66  soilB.GetChild("densityB").GetData(vSoilSampleB);
67  for( size_t i = 0; i < vSoilSampleB.size(); ++i )
68  SampleB += vSoilSampleB[i];
69  SampleB /= vSoilSampleB.size();
70 
71  soilB.GetChild("densityC").GetData(vSoilSampleC);
72  for( size_t i = 0; i < vSoilSampleC.size(); ++i )
73  SampleC += vSoilSampleC[i];
74  SampleC /= vSoilSampleC.size();
75 
76  G4cerr << "fSoilGrammage " << fSoilGrammage/(utl::g/utl::cm2) << " g/cm2"<< G4endl;
77  G4cerr << "SampleA " << SampleA/(utl::g/utl::cm3) << " g/cm3" << G4endl;
78  G4cerr << "SampleB " << SampleB/(utl::g/utl::cm3) << " g/cm3" << G4endl;
79  G4cerr << "SampleC " << SampleC/(utl::g/utl::cm3) << " g/cm3" << G4endl;
80 
81  fSoilDensity = (SampleA + SampleB + SampleC)/3;
82 
83  G4cerr << "fSoilDensity " << fSoilDensity/(utl::g/utl::cm3) << " g/cm3"<< G4endl;
84  G4cerr << "fSoilGrammage/fSoilDensity " << (fSoilGrammage/fSoilDensity)/utl::cm << " cm " << G4endl;
85  G4cerr << "MD_THETASIM_MAX " << MD_THETASIM_MAX/utl::degree << G4endl;
86  G4cerr << "MD_SIM_AREA " << MD_SIM_AREA/utl::m2 << " (m2)" << G4endl;
87 
88  //Nominal area of muon counters (defined in MDetector.h)
89  const double Radius = MD_LARGE_SIDE;//15*utl::m;//sqrt( MD_SIM_AREA / kPi );
90  G4cerr << "Radius " << Radius/utl::m << " (m)" << G4endl;
91  fWorld_Side_Z = fSoilGrammage/fSoilDensity;
92  //X-Y dimensions are defined to allow up to MD_THETASIM_MAX inclided particles
93  fWorld_Side_Y = Radius + fWorld_Side_Z * tan( MD_THETASIM_MAX );
94  fWorld_Side_X = Radius + fWorld_Side_Z * tan( MD_THETASIM_MAX );
95 
96  fSoil_Side_X = fWorld_Side_X;
97  fSoil_Side_Y = fWorld_Side_Y;
98  fSoil_Side_Z = fWorld_Side_Z/2;
99  //Convert from AugerUnits to CLHEP units
100  fSoilGrammage *= ((CLHEP::g/CLHEP::cm2)/(utl::g/utl::cm2));
101  fSoilDensity *= ((CLHEP::g/CLHEP::cm3)/(utl::g/utl::cm3));
102 
103  fWorld_Side_X *= (CLHEP::m / utl::m);
104  fWorld_Side_Y *= (CLHEP::m / utl::m);
105  fWorld_Side_Z *= (CLHEP::m / utl::m);
106 
107  fSoil_Side_X *= (CLHEP::m / utl::m);
108  fSoil_Side_Y *= (CLHEP::m / utl::m);
109  fSoil_Side_Z *= (CLHEP::m / utl::m);
110 
111  G4cerr << "fSoilGrammage " << fSoilGrammage/(CLHEP::g/CLHEP::cm2) << " g/cm2"<< G4endl;
112  G4cerr << "fSoilDensity " << fSoilDensity/(CLHEP::g/CLHEP::cm3) << " g/cm3"<< G4endl;
113  G4cerr << "Depth " << (fSoilGrammage/fSoilDensity)/CLHEP::cm << " cm " << G4endl;
114 
115  G4cerr << "X World size " << fWorld_Side_X/CLHEP::m << " (m)" << G4endl;
116  G4cerr << "Y World size " << fWorld_Side_Y/CLHEP::m << " (m)" << G4endl;
117  G4cerr << "Depth " << fWorld_Side_Z/CLHEP::m << " (m)" << G4endl;
118 
119  //Creates materials after unit conversion
120  CreateMaterials();
121 
122 }
123 
124 SoilConstruction::~SoilConstruction ( ){}
125 
126 G4VPhysicalVolume *
127 SoilConstruction::Construct ( )
128 {
129  if(!fPhysicalWorld){
130  G4ThreeVector placement;
131  INFO("\n++++\nCreating World");
132  //World Box
133 
134  G4double pRMax = fWorld_Side_X;
135 
136  G4Tubs *WorldBox = new G4Tubs("WorldBox", 0, pRMax, fWorld_Side_Z, 0*CLHEP::deg, 360*CLHEP::deg);
137 
138  //G4Box *WorldBox = new G4Box("WorldBox", fWorld_Side_X, fWorld_Side_Y, fWorld_Side_Z);
139  fLogicalWorld = new G4LogicalVolume ( WorldBox, fWorldMaterial, "LogWorld", 0, 0, 0 );
140  fPhysicalWorld = new G4PVPlacement ( 0, G4ThreeVector ( ), fLogicalWorld, "WorldPhys", 0, false, 0 );
141  G4Colour blue ( 0, 0, -1 );
142  G4VisAttributes *fLogicalWorldVisAtt = new G4VisAttributes ( blue );
143  fLogicalWorld->SetVisAttributes ( fLogicalWorldVisAtt );
144  //fLogicalWorld->SetVisAttributes ( G4VisAttributes::Invisible );
145 
146  INFO("\n++++\nCreating Soil");
147  //Soil Box
148  //G4Box *SoilBox = new G4Box ( "SoilBox", fSoil_Side_X, fSoil_Side_Y, fSoil_Side_Z );
149  G4Tubs *SoilBox = new G4Tubs ( "SoilBox", 0, pRMax, fSoil_Side_Z, 0*CLHEP::deg, 360*CLHEP::deg );
150  fLogicalSoil = new G4LogicalVolume ( SoilBox, fSoilMaterial, "SoilLog" );
151  G4double soil_pos_x = 0.0 * utl::m;
152  G4double soil_pos_y = 0.0 * utl::m;
153  G4double soil_pos_z = -fSoil_Side_Z;
154  placement = G4ThreeVector ( soil_pos_x, soil_pos_y, soil_pos_z );
155  fPhysicalSoil = new G4PVPlacement ( 0, placement, fLogicalSoil, "SoilPhys", fLogicalWorld, false, 0 );
156  G4Colour brown ( 0.7, 0.4, 0.1 );
157  G4VisAttributes *fLogicalSoilVisAtt = new G4VisAttributes ( brown );
158  fLogicalSoil->SetVisAttributes ( fLogicalSoilVisAtt );
159  }
160  return fPhysicalWorld;
161 }
162 
163 void
164 SoilConstruction::CreateMaterials ( )
165 {
166  INFO("\n++++\nCreating materials...");
167  G4double a, z;
168  G4int nel, natoms;
169 
170  H = new G4Element ( "Hydrogen" , "H", z = 1., a = 1.01 * CLHEP::g / CLHEP::mole );
171  C = new G4Element ( "Carbon" , "C", z = 6., a = 12.01 * CLHEP::g / CLHEP::mole );
172  N = new G4Element ( "Nitrogen" , "N", z = 7., a = 14.01 * CLHEP::g / CLHEP::mole );
173  O = new G4Element ( "Oxygen" , "O", z = 8., a = 15.99 * CLHEP::g / CLHEP::mole );
174  Na = new G4Element ( "Sodium" , "Na", z = 11, a = 22.98 * CLHEP::g / CLHEP::mole );
175  Mg = new G4Element ( "Magnesium", "Mg", z = 12, a = 24.30 * CLHEP::g / CLHEP::mole );
176  Al = new G4Element ( "Aluminum" , "Al", z = 13, a = 26.98 * CLHEP::g / CLHEP::mole );
177  Si = new G4Element ( "Silicon" , "Si", z = 14, a = 28.08 * CLHEP::g / CLHEP::mole );
178  Cl = new G4Element ( "Chlorine" , "C", z = 17, a = 35.45 * CLHEP::g / CLHEP::mole );
179  K = new G4Element ( "Potassium", "K", z = 19, a = 39.09 * CLHEP::g / CLHEP::mole );
180  Ti = new G4Element ( "Titanium" , "Ti", z = 22, a = 47.86 * CLHEP::g / CLHEP::mole );
181  Ca = new G4Element ( "Calcium" , "Ca", z = 20, a = 40.07 * CLHEP::g / CLHEP::mole );
182  Fe = new G4Element ( "Iron" , "Fe", z = 26, a = 55.85 * CLHEP::g / CLHEP::mole );
183 
184  SiO2 = new G4Material ( "SiO2", fSoilDensity, nel = 2 );
185  SiO2->AddElement ( Si, natoms = 1 );
186  SiO2->AddElement ( O , natoms = 2 );
187 
188  Al2O3 = new G4Material ( "Al2O3", fSoilDensity, nel = 2 );
189  Al2O3->AddElement ( Al, natoms = 2 );
190  Al2O3->AddElement ( O , natoms = 3 );
191 
192  Fe2O3 = new G4Material ( "Fe2O3", fSoilDensity, nel = 2 );
193  Fe2O3->AddElement ( Fe, natoms = 2 );
194  Fe2O3->AddElement ( O , natoms = 3 );
195 
196  TiO2 = new G4Material ( "TiO2", fSoilDensity, nel = 2 );
197  TiO2->AddElement ( Ti, natoms = 1 );
198  TiO2->AddElement ( O , natoms = 2 );
199 
200  CaO = new G4Material ( "CaO", fSoilDensity, nel = 2 );
201  CaO->AddElement ( Ca, natoms = 1 );
202  CaO->AddElement ( O , natoms = 1 );
203 
204  MgO = new G4Material ( "MgO", fSoilDensity, nel = 2 );
205  MgO->AddElement ( Mg, natoms = 1 );
206  MgO->AddElement ( O , natoms = 1 );
207 
208  K2O = new G4Material ( "K2O", fSoilDensity, nel = 2 );
209  K2O->AddElement ( K, natoms = 2 );
210  K2O->AddElement ( O , natoms = 1 );
211 
212  Na2O = new G4Material ( "Na2O", fSoilDensity, nel = 2 );
213  Na2O->AddElement ( Na, natoms = 2 );
214  Na2O->AddElement ( O , natoms = 1 );
215 
216  //Malargue Soil
217  fSoilMaterial = new G4Material ( "MalargueSoil", fSoilDensity , nel = 8 );
218  fSoilMaterial->AddMaterial ( SiO2, 66.0 * CLHEP::perCent );
219  fSoilMaterial->AddMaterial ( Al2O3, 13.0 * CLHEP::perCent );
220  fSoilMaterial->AddMaterial ( Fe2O3, 6.00 * CLHEP::perCent );
221  fSoilMaterial->AddMaterial ( TiO2, 0.70 * CLHEP::perCent );
222  fSoilMaterial->AddMaterial ( CaO, 6.00 * CLHEP::perCent );
223  fSoilMaterial->AddMaterial ( MgO, 1.00 * CLHEP::perCent );
224  fSoilMaterial->AddMaterial ( K2O, 2.60 * CLHEP::perCent );
225  fSoilMaterial->AddMaterial ( Na2O, 4.70 * CLHEP::perCent );
226  //World material (Air)
227  fWorldMaterial = new G4Material ( "Air_MdGround", 1.29 * CLHEP::mg / CLHEP::cm3, nel = 2 );
228  fWorldMaterial->AddElement ( N, 70. * CLHEP::perCent );
229  fWorldMaterial->AddElement ( O, 30. * CLHEP::perCent );
230 
231  return;
232 }
233 
234 G4double
235 SoilConstruction::GetDepth ( )
236 {
237  return fWorld_Side_Z;
238 }
utl::perCent
constexpr double perCent
Definition: AugerUnits.h:282
utl::cm3
constexpr double cm3
Definition: AugerUnits.h:119
MD_SIM_AREA
#define MD_SIM_AREA
Definition: SoilConstruction.cc:27
INFO
#define INFO(message)
Macro for logging informational messages.
Definition: ErrorLogger.h:161
utl::m2
constexpr double m2
Definition: AugerUnits.h:122
utl::Branch::GetChild
Branch GetChild(const std::string &childName) const
Get child of this Branch by child name.
Definition: Branch.cc:211
utl::deg
constexpr double deg
Definition: AugerUnits.h:140
utl::Branch
Class representing a document branch.
Definition: Branch.h:107
utl::mole
constexpr double mole
Definition: AugerUnits.h:262
MD_THETASIM_MAX
#define MD_THETASIM_MAX
Definition: SoilConstruction.cc:26
utl::degree
constexpr double degree
Definition: Tools/CIC/utl/Math.h:18
utl::g
constexpr double g
Definition: AugerUnits.h:200
MD_LARGE_SIDE
#define MD_LARGE_SIDE
Definition: SoilConstruction.cc:28
utl::Branch::GetData
void GetData(bool &b) const
Overloads of the GetData member template function.
Definition: Branch.cc:644
utl::mg
constexpr double mg
Definition: AugerUnits.h:201
utl::cm
constexpr double cm
Definition: AugerUnits.h:117
fwk::CentralConfig
Main configuration utility.
Definition: CentralConfig.h:51
GroundPropagatorAG::SoilConstruction::fLogicalWorld
G4LogicalVolume * fLogicalWorld
Volumes.
Definition: SoilConstruction.h:28
utl::m
constexpr double m
Definition: AugerUnits.h:121
fwk::CentralConfig::GetTopBranch
utl::Branch GetTopBranch(const std::string &id)
Get top branch for moduleConfigLink with given id (XML files)
Definition: CentralConfig.cc:229
TelescopeSimulatorLX::blue
static const G4Colour blue(0.0, 0.0, 1.0)
utl::cm2
constexpr double cm2
Definition: AugerUnits.h:118
, generated on Tue Sep 26 2023.