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