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