| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26
27 // (adapted from B1DetectorConstruction)
28 // Author: A.Knaian (ara@nklabs.com), N.MacFad
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 po
68
69 DetectorConstruction::DetectorConstruction()
70 : G4VUserDetectorConstruction(),
71 fScoringVolume(nullptr)
72 {
73 fMessenger = new DetectorConstructionMessenger
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::Const
86 {
87 //
88 // Check cloud build settings
89 //
90 if (fFastAerosolCloud + fParameterisedCloud +
91 {
92 std::ostringstream message;
93 message << "Must select at most one build ty
94 << " fFastAerosolCloud = " << fFastAer
95 << " fParameterisedCloud = " << fParam
96 << " fSmoothCloud = " << fSmoothCloud
97 G4Exception("DetectorConstruction::Constr
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 volume
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*cl
123 }
124 else if (fCloudShapeStr == "ellipsoid")
125 {
126 G4cout << "Cloud shape = ellipsoid" << G4en
127 fCloudShape = new G4Ellipsoid("cloudShape",
128 }
129 else if (fCloudShapeStr == "cylinder")
130 {
131 G4cout << "Cloud shape = cylinder" << G4endl;
132 fCloudShape = new G4Tubs("cloudShape", 0.0,
133 }
134 else if (fCloudShapeStr == "pipe")
135 {
136 G4cout << "Cloud shape = pipe" << G4endl;
137 fCloudShape = new G4Tubs("cloudShape", 0.25*
138 }
139 else
140 {
141 std::ostringstream message;
142 message << "Invalid cloud shape = " << fClou
143 G4Exception("DetectorConstruction::Construct
144 FatalException, message);
145 }
146
147 //
148 // Droplet Shape
149 //
150
151 // The difference in radii of the maximal sphe
152 G4double sphericalUncertainty = 0.0;
153 if (fDropletShapeStr == "sphere")
154 {
155 G4cout << "Droplet shape = sphere" << G4endl
156 fDropletShape = new G4Orb("dropletSV", fDrop
157 sphericalUncertainty = 0.0;
158 }
159 else if (fDropletShapeStr == "halfSphere")
160 {
161 G4cout << "Droplet shape = halfSphere" << G4
162 fDropletShape = new G4Sphere("dropletSV", 0.
163 sphericalUncertainty = fDropletR;
164 }
165 else if (fDropletShapeStr == "cylinder")
166 {
167 G4cout << "Droplet shape = cylinder" << G4en
168 fDropletShape = new G4Tubs("dropletSV", 0, f
169 sphericalUncertainty = fDropletR*(1-1/std::s
170 }
171 else if (fDropletShapeStr == "box")
172 {
173 G4cout << "Droplet shape = box" << G4endl;
174 fDropletShape = new G4Box("dropletSV", fDrop
175 sphericalUncertainty = fDropletR*(1-1/std::s
176 }
177 else
178 {
179 std::ostringstream message;
180 message << "Invalid droplet shape = " << fC
181 G4Exception("DetectorConstruction::Constru
182 FatalException, message);
183 }
184
185 //
186 // Materials
187 //
188 // Compute the density of air at 14 km using t
189 // see, e.g., https://en.wikipedia.org/wiki/De
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*L
199 G4double air_density = p*M/(R*T);
200
201 // make materials and set densities
202 G4Material* air_mat = nist->BuildMaterialWithN
203 G4Material* water_mat = nist->FindOrBuildMater
204
205 G4double water_density = water_mat->GetDensity
206 G4double ice_density = 0.9168*g/cm3;
207
208 G4Material* ice_mat = new G4Material("Water ic
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*(fClo
218
219 G4double droplet_volume = fDropletShape->GetCu
220 G4double droplet_total_volume = droplet_count*
221
222 G4double droplet_total_mass = droplet_total_v
223
224 //
225 // Cloud macroscopic quantities
226 //
227 G4double cloud_volume = fCloudShape->GetCubicV
228 G4double cloud_air_volume = cloud_volume - dro
229 G4double cloud_air_mass = air_density*cloud_ai
230
231 //
232 // Step limit
233 //
234 fStepLimits = new G4UserLimits(fStepLim);
235
236 //
237 // Build world
238 //
239 G4Box* solidWorld =new G4Box("World",//its nam
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 soli
246 air_mat, //its materia
247 "World");//its name
248
249 logicWorld->SetUserLimits(fStepLimits);
250
251 G4VPhysicalVolume* physWorld = new G4PVPlaceme
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 geome
268 //
269 // *******************************************
270
271 if (fFastAerosolCloud)
272 {
273 G4cout << "\nFastAerosol geometry with n=" <<
274
275 fCloud = new FastAerosol("cloud",
276 fDropletR, //bounding radius of droplets
277 fMinSpacing,//minimum spaci
278 fDropletNumDe
279 sphericalUncertainty); //un
280 fCloud->SetDropletsPerVoxel(4);
281
282 FastAerosolSolid* solidCloud = new FastAerosol
283 fCloud, //its shape
284 fDropletShape); //its droplets
285
286 solidCloud->SetStepLim(fStepLim);
287 //FastAerosol can use step limit to speed calc
288
289 logicCloud = new G4LogicalVolume(solidCloud,//
290 droplet_mat,//its m
291 "cloudLV");//its name
292 logicCloud->SetUserLimits(fStepLimits);
293 logicCloud->SetVisAttributes(G4VisAttributes(G
294
295 new G4PVPlacement(nullptr, G4ThreeVector(), lo
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 voxe
304 if (fPrePopulate)
305 {
306 // populate (proving it to the user by printi
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
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,
327 replace( rStr.begin(), rStr.end(), '.', 'p');
328 if (rStr.back() == 'p') { rStr.pop_back(); }
329
330 // want to represent the number density as 1E
331 G4int order10 = (G4int) -round(10*std::log10(f
332 G4int leading = order10 / 10; // first number
333 G4int trailing = order10 % 10; // second numb
334 G4String nStr = "1E-" + std::to_string(leading
335
336 // save population time
337 std::ofstream file;
338 file.open("popTime_r" + rStr + "mm_n" + nStr +
339 file << ((float)t)/CLOCKS_PER_SEC;
340 file.close();
341
342 // save distribution
343 G4String fName = "distribution_r" + rStr + "mm
344 fCloud->SaveToFile(fName);
345 }
346 }
347 // *******************************************
348 //
349 // (For comparision/benchmarking) Build the cl
350 //
351 // *******************************************
352 // the droplet positions for this cloud are th
353 // this is to make comparable simulations betw
354 // this requires that we first simulate FastAe
355
356 else if (fParameterisedCloud)
357 {
358 G4cout << "\nParameterised geometry with n="
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') +
366 replace( rStr.begin(), rStr.end(), '.', 'p
367 if (rStr.back() == 'p') { rStr.pop_back();
368
369 // want to represent the number density as 1E
370 G4int order10 = (G4int) -round(10*std::log10(
371 G4int leading = order10 / 10; // first numbe
372 G4int trailing = order10 % 10; // second numb
373 G4String nStr = "1E-" + std::to_string(leadin
374
375 fName = "distribution_r" + rStr + "mm_n" + nS
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("cloudBoundi
393 0.5*cloud_sizeXY,//half
394 0.5*cloud_sizeXY, //half y-span
395 0.5*cloud_sizeZ); //half z-span
396
397 logicCloud = new G4LogicalVolume(cloudBoundin
398 air_mat,//its material
399 "cloudLV");//its name
400
401 logicCloud->SetSmartless(fSmartless);
402 logicCloud->SetUserLimits(fStepLimits);
403 logicCloud->SetVisAttributes(G4VisAttributes(
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
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
425 positions.size(),//number of dropl
426 cloudParam);//the parametrisation
427 }
428 // *******************************************
429 //
430 // (For comparision/benchmarking) Simulate the
431 //
432 // *******************************************
433 else if (fSmoothCloud)
434 {
435 G4cout << "\nSmooth geometry based on a clou
436 // build cloud by smearing the droplets unifor
437 G4Material* cloud_mat = new G4Material("Cloud"
438
439 cloud_mat->AddMaterial(droplet_mat, droplet_to
440
441 cloud_mat->AddMaterial(air_mat, cloud_air_mass
442
443 logicCloud = new G4LogicalVolume(fCloudShape,
444 cloud_mat, //its material
445 "cloudLV"); //its name
446
447 logicCloud->SetUserLimits(fStepLimits);
448 logicCloud->SetVisAttributes(G4VisAttributes(
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->FindOrBuildM
470 G4ThreeVector detector_pos = G4ThreeVector(0,
471
472 G4Box* soldDetector = new G4Box("detectorSV"
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 so
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 mothervolu
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