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