| Geant4 Cross Reference |
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 // (adapted from B1DetectorConstruction)
28 // Author: A.Knaian (ara@nklabs.com), N.MacFadden (natemacfadden@gmail.com)
29
30 #include "FADetectorConstruction.hh"
31 #include "FADetectorConstructionMessenger.hh"
32
33 #include "G4RunManager.hh"
34 #include "G4NistManager.hh"
35
36 #include "G4LogicalVolume.hh"
37 #include "G4PVPlacement.hh"
38 #include "G4SystemOfUnits.hh"
39
40 // shapes
41 #include "G4Box.hh"
42 #include "G4Cons.hh"
43 #include "G4Orb.hh"
44 #include "G4Sphere.hh"
45 #include "G4Trd.hh"
46 #include "G4Tubs.hh"
47 #include "G4Ellipsoid.hh"
48
49 // to build FastAerosol cloud
50 #include "FastAerosolSolid.hh"
51
52 // to build parameterised cloud
53 #include "FACloudParameterisation.hh"
54 #include "G4PVParameterised.hh"
55 #include <fstream>
56
57 // step limits
58 #include "G4UserLimits.hh"
59
60 // visualization
61 #include "G4VisAttributes.hh"
62 #include "G4Colour.hh"
63
64 // to save distribution
65 #include <sys/stat.h>
66 #include <ctime>
67 // for measuring FastAerosol droplet center population time
68
69 DetectorConstruction::DetectorConstruction()
70 : G4VUserDetectorConstruction(),
71 fScoringVolume(nullptr)
72 {
73 fMessenger = new DetectorConstructionMessenger(this);
74 }
75
76 DetectorConstruction::~DetectorConstruction()
77 {
78 delete fMessenger;
79 delete fStepLimits;
80 delete fCloudShape;
81 delete fDropletShape;
82 delete fCloud;
83 }
84
85 G4VPhysicalVolume* DetectorConstruction::Construct()
86 {
87 //
88 // Check cloud build settings
89 //
90 if (fFastAerosolCloud + fParameterisedCloud + fSmoothCloud > 1)
91 {
92 std::ostringstream message;
93 message << "Must select at most one build type! Selections:" << G4endl
94 << " fFastAerosolCloud = " << fFastAerosolCloud << G4endl
95 << " fParameterisedCloud = " << fParameterisedCloud << G4endl
96 << " fSmoothCloud = " << fSmoothCloud << G4endl;
97 G4Exception("DetectorConstruction::Construct()", "GeomSolids0002",
98 FatalException, message);
99 }
100 //
101 // Get nist material manager
102 //
103 G4NistManager* nist = G4NistManager::Instance();
104 //
105 // Option to switch on/off checking of volumes overlaps
106 //
107 G4bool checkOverlaps = false;
108 //
109 // Large scale geometry dimensions
110 //
111 G4double cloud_sizeXY = 0.5*m;
112 G4double cloud_sizeZ = 5.0*m;
113
114 G4double world_sizeXY = 1.1*(cloud_sizeXY);
115 G4double world_sizeZ= 1.1*(cloud_sizeZ);
116 //
117 // Cloud shape
118 //
119 if (fCloudShapeStr == "box")
120 {
121 G4cout << "Cloud shape = box" << G4endl;
122 fCloudShape = new G4Box("cloudShape", 0.5*cloud_sizeXY, 0.5*cloud_sizeXY, 0.5*cloud_sizeZ);
123 }
124 else if (fCloudShapeStr == "ellipsoid")
125 {
126 G4cout << "Cloud shape = ellipsoid" << G4endl;
127 fCloudShape = new G4Ellipsoid("cloudShape", 0.5*cloud_sizeXY, 0.5*cloud_sizeXY, 0.5*cloud_sizeZ, 0, 0);
128 }
129 else if (fCloudShapeStr == "cylinder")
130 {
131 G4cout << "Cloud shape = cylinder" << G4endl;
132 fCloudShape = new G4Tubs("cloudShape", 0.0, 0.5*cloud_sizeXY, 0.5*cloud_sizeZ, 0, 360*deg);
133 }
134 else if (fCloudShapeStr == "pipe")
135 {
136 G4cout << "Cloud shape = pipe" << G4endl;
137 fCloudShape = new G4Tubs("cloudShape", 0.25*cloud_sizeXY, 0.5*cloud_sizeXY, 0.5*cloud_sizeZ, 0, 360.*deg);
138 }
139 else
140 {
141 std::ostringstream message;
142 message << "Invalid cloud shape = " << fCloudShapeStr << "!";
143 G4Exception("DetectorConstruction::Construct()", "GeomSolids0002",
144 FatalException, message);
145 }
146
147 //
148 // Droplet Shape
149 //
150
151 // The difference in radii of the maximal sphere (centered at the origin) contained in the droplet and the minimal sphere (centered at the origin) containing the droplet
152 G4double sphericalUncertainty = 0.0;
153 if (fDropletShapeStr == "sphere")
154 {
155 G4cout << "Droplet shape = sphere" << G4endl;
156 fDropletShape = new G4Orb("dropletSV", fDropletR);
157 sphericalUncertainty = 0.0;
158 }
159 else if (fDropletShapeStr == "halfSphere")
160 {
161 G4cout << "Droplet shape = halfSphere" << G4endl;
162 fDropletShape = new G4Sphere("dropletSV", 0.0, fDropletR, 0.0, 180.*deg, 0.0, 180.*deg);
163 sphericalUncertainty = fDropletR;
164 }
165 else if (fDropletShapeStr == "cylinder")
166 {
167 G4cout << "Droplet shape = cylinder" << G4endl;
168 fDropletShape = new G4Tubs("dropletSV", 0, fDropletR/std::sqrt(3), fDropletR/std::sqrt(3), 0, 360.*deg);
169 sphericalUncertainty = fDropletR*(1-1/std::sqrt(3));
170 }
171 else if (fDropletShapeStr == "box")
172 {
173 G4cout << "Droplet shape = box" << G4endl;
174 fDropletShape = new G4Box("dropletSV", fDropletR/std::sqrt(3), fDropletR/std::sqrt(3), fDropletR/std::sqrt(3));
175 sphericalUncertainty = fDropletR*(1-1/std::sqrt(3));
176 }
177 else
178 {
179 std::ostringstream message;
180 message << "Invalid droplet shape = " << fCloudShapeStr << "!";
181 G4Exception("DetectorConstruction::Construct()", "GeomSolids0002",
182 FatalException, message);
183 }
184
185 //
186 // Materials
187 //
188 // Compute the density of air at 14 km using the Barometric formula
189 // see, e.g., https://en.wikipedia.org/wiki/Density_of_air
190 G4double h = 14.0*km;
191 G4double p0 = 101325*hep_pascal;
192 G4double T0 = 288.15*kelvin;
193 G4double grav = 9.80665*m/(s*s);
194 G4double La = 0.0065*kelvin/m;
195 G4double R = 8.31447*joule/(mole*kelvin);
196 G4double M = 0.0289644*kg/mole;
197 G4double T = T0 - La*h;
198 G4double p = p0*std::pow(1-La*h/T0,grav*M/(R*La));
199 G4double air_density = p*M/(R*T);
200
201 // make materials and set densities
202 G4Material* air_mat = nist->BuildMaterialWithNewDensity("Atmosphere","G4_AIR",air_density);
203 G4Material* water_mat = nist->FindOrBuildMaterial("G4_WATER");
204
205 G4double water_density = water_mat->GetDensity();
206 G4double ice_density = 0.9168*g/cm3;
207
208 G4Material* ice_mat = new G4Material("Water ice ", ice_density, 1, kStateSolid, T, p);
209 ice_mat->AddMaterial(water_mat, 1.);
210
211 //
212 // Droplets
213 //
214 G4double droplet_density = water_density;
215 G4Material* droplet_mat = water_mat;
216
217 G4double droplet_count = fDropletNumDens*(fCloudShape->GetCubicVolume());
218
219 G4double droplet_volume = fDropletShape->GetCubicVolume();
220 G4double droplet_total_volume = droplet_count*droplet_volume;
221
222 G4double droplet_total_mass = droplet_total_volume*droplet_density;
223
224 //
225 // Cloud macroscopic quantities
226 //
227 G4double cloud_volume = fCloudShape->GetCubicVolume();
228 G4double cloud_air_volume = cloud_volume - droplet_total_volume;
229 G4double cloud_air_mass = air_density*cloud_air_volume;
230
231 //
232 // Step limit
233 //
234 fStepLimits = new G4UserLimits(fStepLim);
235
236 //
237 // Build world
238 //
239 G4Box* solidWorld =new G4Box("World",//its name
240 0.5*world_sizeXY,//half x-span
241 0.5*world_sizeXY,//half y-span
242 0.5*world_sizeZ);//half z-span
243
244 G4LogicalVolume* logicWorld =
245 new G4LogicalVolume(solidWorld, //its solid
246 air_mat, //its material
247 "World");//its name
248
249 logicWorld->SetUserLimits(fStepLimits);
250
251 G4VPhysicalVolume* physWorld = new G4PVPlacement(nullptr,//no rotation
252 G4ThreeVector(),//at (0,0,0)
253 logicWorld,//its logical volume
254 "World",//its name
255 nullptr,//its mothervolume
256 false,//no boolean operation
257 0, //copy number
258 checkOverlaps); //overlaps checking
259
260 //
261 // Build cloud
262 //
263 G4LogicalVolume* logicCloud;
264
265 // **********************************************************
266 //
267 // Build the cloud using the FastAerosol geometry class
268 //
269 // ***********************************************************
270
271 if (fFastAerosolCloud)
272 {
273 G4cout << "\nFastAerosol geometry with n=" << fDropletNumDens*mm3 << "/mm3, r=" << fDropletR/mm << "mm spheres.\n" << G4endl;
274
275 fCloud = new FastAerosol("cloud", fCloudShape,//cloud shape
276 fDropletR, //bounding radius of droplets
277 fMinSpacing,//minimum spacing between droplets
278 fDropletNumDens, //approximate number of droplets in cloud
279 sphericalUncertainty); //uncertainty in distance to droplet surface from outside using just droplet's origin as info
280 fCloud->SetDropletsPerVoxel(4);
281
282 FastAerosolSolid* solidCloud = new FastAerosolSolid("cloudSV",//its name
283 fCloud, //its shape
284 fDropletShape); //its droplets
285
286 solidCloud->SetStepLim(fStepLim);
287 //FastAerosol can use step limit to speed calculations
288
289 logicCloud = new G4LogicalVolume(solidCloud,//its solid
290 droplet_mat,//its material
291 "cloudLV");//its name
292 logicCloud->SetUserLimits(fStepLimits);
293 logicCloud->SetVisAttributes(G4VisAttributes(G4Colour(0.0,0.0,1.0,0.4)));
294
295 new G4PVPlacement(nullptr, G4ThreeVector(), logicCloud, "cloudPV", //its name
296 logicWorld,//its mother volume
297 false, //no boolean operation
298 0,//copy number
299 checkOverlaps);//overlaps checking
300
301 fCloud->SetSeed(fCloudSeed);
302
303 // fPrePopulate = whether to populate all voxels at the beginning or on the fly
304 if (fPrePopulate)
305 {
306 // populate (proving it to the user by printing population reports)
307 clock_t t;
308 t = clock();
309
310 G4cout << "\nBefore populating" << G4endl;
311 G4cout << "=================" << G4endl;
312 fCloud->PrintPopulationReport();
313 G4cout << "\nPopulating..." << G4endl;
314 fCloud->PopulateAllGrids();
315 G4cout << "\nAfter populating" << G4endl;
316 G4cout << "================" << G4endl;
317 fCloud->PrintPopulationReport();
318 G4cout << G4endl;
319
320 t = clock() - t;
321
322 G4cout << "\nThis took " << ((float)t)/CLOCKS_PER_SEC << "s\n" << G4endl;
323
324 // make filename variables to save data
325 G4String rStr = std::to_string(fDropletR/mm);
326 rStr.erase ( rStr.find_last_not_of('0') + 1, std::string::npos ); // drop trailing 0
327 replace( rStr.begin(), rStr.end(), '.', 'p');
328 if (rStr.back() == 'p') { rStr.pop_back(); } // don't write "3p" for 3.0, just write "3"
329
330 // want to represent the number density as 1E-ApB for some A, B
331 G4int order10 = (G4int) -round(10*std::log10(fDropletNumDens*mm3)); // gives 10x the exponent rounded to the int (10x so we get two decimals)
332 G4int leading = order10 / 10; // first number
333 G4int trailing = order10 % 10; // second number
334 G4String nStr = "1E-" + std::to_string(leading) + "p" + std::to_string(trailing);
335
336 // save population time
337 std::ofstream file;
338 file.open("popTime_r" + rStr + "mm_n" + nStr + "mm-3.csv");
339 file << ((float)t)/CLOCKS_PER_SEC;
340 file.close();
341
342 // save distribution
343 G4String fName = "distribution_r" + rStr + "mm_n" + nStr + "mm-3.csv";
344 fCloud->SaveToFile(fName);
345 }
346 }
347 // **********************************************************
348 //
349 // (For comparision/benchmarking) Build the cloud using G4VParameterized (does not use FastAerosol)
350 //
351 // ***********************************************************
352 // the droplet positions for this cloud are those saved in the "distribution" folder of our data
353 // this is to make comparable simulations between FastAerosol and parameterised clouds
354 // this requires that we first simulate FastAerosol (pre-populated) to generate the positions
355
356 else if (fParameterisedCloud)
357 {
358 G4cout << "\nParameterised geometry with n=" << fDropletNumDens*mm3 << "/mm3 and r=" << fDropletR/mm << "mm spheres.\n" << G4endl;
359 std::vector<G4ThreeVector> positions;
360 G4double x,y,z;
361
362 // load distribution file
363 G4String fName;
364 G4String rStr = std::to_string(fDropletR/mm);
365 rStr.erase ( rStr.find_last_not_of('0') + 1, std::string::npos ); // drop trailing 0
366 replace( rStr.begin(), rStr.end(), '.', 'p');
367 if (rStr.back() == 'p') { rStr.pop_back(); } // don't write "3p" for 3.0, just write "3"
368
369 // want to represent the number density as 1E-ApB for some A, B
370 G4int order10 = (G4int) -round(10*std::log10(fDropletNumDens*mm3)); // gives 10x the exponent rounded to the int (10x so we get two decimals)
371 G4int leading = order10 / 10; // first number
372 G4int trailing = order10 % 10; // second number
373 G4String nStr = "1E-" + std::to_string(leading) + "p" + std::to_string(trailing);
374
375 fName = "distribution_r" + rStr + "mm_n" + nStr + "mm-3.csv";
376 std::ifstream infile(fName);
377 std::string line;
378
379 while (getline(infile,line)) {
380 std::istringstream stream(line);
381 std::string field;
382 getline(stream,field,','); x = stod(field)*mm;
383 getline(stream,field,','); y = stod(field)*mm;
384 getline(stream,field,','); z = stod(field)*mm;
385
386 positions.push_back(G4ThreeVector(x,y,z));
387 }
388
389 G4VPVParameterisation* cloudParam =
390 new CloudParameterisation(positions);
391
392 G4Box* cloudBounding = new G4Box("cloudBounding",//its name
393 0.5*cloud_sizeXY,//half x-span
394 0.5*cloud_sizeXY, //half y-span
395 0.5*cloud_sizeZ); //half z-span
396
397 logicCloud = new G4LogicalVolume(cloudBounding, //its solid
398 air_mat,//its material
399 "cloudLV");//its name
400
401 logicCloud->SetSmartless(fSmartless);
402 logicCloud->SetUserLimits(fStepLimits);
403 logicCloud->SetVisAttributes(G4VisAttributes(false));
404
405 new G4PVPlacement(nullptr,//no rotation
406 G4ThreeVector(),//at position
407 logicCloud,//its logical volume
408 "cloudPV",//its name
409 logicWorld, //its mothervolume
410 false, //no boolean operation
411 0,//copy number
412 checkOverlaps);//overlaps checking
413
414 G4LogicalVolume* logicDroplet =
415 new G4LogicalVolume(fDropletShape,//its solid
416 droplet_mat,//its material
417 "dropletLV");//its name
418
419 logicDroplet->SetUserLimits(fStepLimits);
420
421 new G4PVParameterised("droplets",//its name
422 logicDroplet,//droplet logical volume
423 logicCloud,//mother logical volume
424 kUndefined,//droplets placed along this axis
425 positions.size(),//number of droplets
426 cloudParam);//the parametrisation
427 }
428 // **********************************************************
429 //
430 // (For comparision/benchmarking) Simulate the cloud by smearing droplets out into a single solid (does not use FastAerosol)
431 //
432 // ***********************************************************
433 else if (fSmoothCloud)
434 {
435 G4cout << "\nSmooth geometry based on a cloud of n=" << fDropletNumDens*mm3 << "/mm3 and r=" << fDropletR/mm << "mm spheres.\n" << G4endl;
436 // build cloud by smearing the droplets uniformly across the cloud volume, for comparison/benchmarking purposes (does not use FastAerosol)
437 G4Material* cloud_mat = new G4Material("Cloud", (droplet_total_mass+cloud_air_mass)/cloud_volume, 2);
438
439 cloud_mat->AddMaterial(droplet_mat, droplet_total_mass/(cloud_air_mass+droplet_total_mass));
440
441 cloud_mat->AddMaterial(air_mat, cloud_air_mass/(cloud_air_mass+droplet_total_mass));
442
443 logicCloud = new G4LogicalVolume(fCloudShape, //its solid
444 cloud_mat, //its material
445 "cloudLV"); //its name
446
447 logicCloud->SetUserLimits(fStepLimits);
448 logicCloud->SetVisAttributes(G4VisAttributes(G4Colour(0.0,0.0,1.0,0.4)));
449
450 new G4PVPlacement(nullptr, //no rotation
451 G4ThreeVector(), //at position
452 logicCloud,//its logical volume
453 "cloudPV", //its name
454 logicWorld, //its mothervolume
455 false, //no boolean operation
456 0,//copy number
457 checkOverlaps); //overlaps checking
458 }
459 else
460 {
461 G4cout << "\nNo cloud.\n" << G4endl;
462 }
463
464 //
465 // Build detector
466 //
467 G4double detector_sizeXY = cloud_sizeXY;
468 G4double detector_sizeZ = 0.05*m;
469 G4Material* detector_mat = nist->FindOrBuildMaterial("G4_Al");
470 G4ThreeVector detector_pos = G4ThreeVector(0, 0, 0.5*1.05*cloud_sizeZ);
471
472 G4Box* soldDetector = new G4Box("detectorSV", //its name
473 0.5*detector_sizeXY, //half x-span
474 0.5*detector_sizeXY, //half y-span
475 0.5*detector_sizeZ); //half z-span
476
477 G4LogicalVolume* logicDetector =
478 new G4LogicalVolume(soldDetector, //its solid
479 detector_mat, //its material
480 "detectorLV"); //its name
481
482 logicDetector->SetUserLimits(fStepLimits);
483
484 new G4PVPlacement(nullptr,//no rotation
485 detector_pos, //at position
486 logicDetector, //its logical volume
487 "detectorPV", //its name
488 logicWorld, //its mothervolume
489 false, //no boolean operation
490 0, //copy number
491 checkOverlaps); //overlaps checking
492 //
493 // Scoring Volume
494 //
495 fScoringVolume = logicDetector;
496
497 return physWorld;
498 }
499
500