fastAerosol/src/FADetectorConstruction.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026
0027 // (adapted from B1DetectorConstruction)
0028 // Author: A.Knaian (ara@nklabs.com), N.MacFadden (natemacfadden@gmail.com)
0029
0030 #include "FADetectorConstruction.hh"
0031 #include "FADetectorConstructionMessenger.hh"
0032
0033 #include "G4RunManager.hh"
0034 #include "G4NistManager.hh"
0035
0036 #include "G4LogicalVolume.hh"
0037 #include "G4PVPlacement.hh"
0038 #include "G4SystemOfUnits.hh"
0039
0040 // shapes
0041 #include "G4Box.hh"
0042 #include "G4Cons.hh"
0043 #include "G4Orb.hh"
0044 #include "G4Sphere.hh"
0045 #include "G4Trd.hh"
0046 #include "G4Tubs.hh"
0047 #include "G4Ellipsoid.hh"
0048
0049 // to build FastAerosol cloud
0050 #include "FastAerosolSolid.hh"
0051
0052 // to build parameterised cloud
0053 #include "FACloudParameterisation.hh"
0054 #include "G4PVParameterised.hh"
0055 #include <fstream>
0056
0057 // step limits
0058 #include "G4UserLimits.hh"
0059
0060 // visualization
0061 #include "G4VisAttributes.hh"
0062 #include "G4Colour.hh"
0063
0064 // to save distribution
0065 #include <sys/stat.h>
0066 #include <ctime>
0067 // for measuring FastAerosol droplet center population time
0068
0069 DetectorConstruction::DetectorConstruction()
0070 : G4VUserDetectorConstruction(),
0071 fScoringVolume(nullptr)
0072 {
0073 fMessenger = new DetectorConstructionMessenger(this);
0074 }
0075
0076 DetectorConstruction::~DetectorConstruction()
0077 {
0078  delete fMessenger;
0079  delete fStepLimits;
0080  delete fCloudShape;
0081  delete fDropletShape;
0082  delete fCloud;
0083 }
0084
0085 G4VPhysicalVolume* DetectorConstruction::Construct()
0086 {
0087  //
0088  // Check cloud build settings
0089  //
0090  if (fFastAerosolCloud + fParameterisedCloud + fSmoothCloud > 1)
0091  {
0092  std::ostringstream message;
0093   message << "Must select at most one build type! Selections:" << G4endl
0094       << "     fFastAerosolCloud = " << fFastAerosolCloud << G4endl
0095       << "     fParameterisedCloud = " << fParameterisedCloud << G4endl
0096       << "     fSmoothCloud = " << fSmoothCloud << G4endl;
0097        G4Exception("DetectorConstruction::Construct()", "GeomSolids0002",
0098                     FatalException, message);
0099     }
0100  //
0101  // Get nist material manager
0102  //
0103  G4NistManager* nist = G4NistManager::Instance();
0104  //
0105  // Option to switch on/off checking of volumes overlaps
0106  //
0107  G4bool checkOverlaps = false;
0108  //
0109  // Large scale geometry dimensions
0110  //
0111  G4double cloud_sizeXY = 0.5*m;
0112  G4double cloud_sizeZ = 5.0*m;
0113
0114  G4double world_sizeXY = 1.1*(cloud_sizeXY);
0115  G4double world_sizeZ= 1.1*(cloud_sizeZ);
0116  //
0117  // Cloud shape
0118  //
0119  if (fCloudShapeStr == "box")
0120  {
0121  G4cout << "Cloud shape = box" << G4endl;
0122   fCloudShape = new G4Box("cloudShape", 0.5*cloud_sizeXY, 0.5*cloud_sizeXY, 0.5*cloud_sizeZ);
0123  }
0124  else if (fCloudShapeStr == "ellipsoid")
0125  {
0126    G4cout << "Cloud shape = ellipsoid" << G4endl;
0127    fCloudShape = new G4Ellipsoid("cloudShape", 0.5*cloud_sizeXY, 0.5*cloud_sizeXY, 0.5*cloud_sizeZ, 0, 0);
0128  }
0129  else if (fCloudShapeStr == "cylinder")
0130  {
0131  G4cout << "Cloud shape = cylinder" << G4endl;
0132   fCloudShape = new G4Tubs("cloudShape", 0.0, 0.5*cloud_sizeXY, 0.5*cloud_sizeZ, 0, 360*deg);
0133     }
0134  else if (fCloudShapeStr == "pipe")
0135  {
0136   G4cout << "Cloud shape = pipe" << G4endl;
0137   fCloudShape = new G4Tubs("cloudShape", 0.25*cloud_sizeXY, 0.5*cloud_sizeXY, 0.5*cloud_sizeZ, 0, 360.*deg);
0138  }
0139  else
0140  {
0141   std::ostringstream message;
0142   message << "Invalid cloud shape = " << fCloudShapeStr << "!";
0143   G4Exception("DetectorConstruction::Construct()", "GeomSolids0002",
0144   FatalException, message);
0145     }
0146
0147  //
0148  // Droplet Shape
0149  //
0150
0151 // 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
0152  G4double sphericalUncertainty = 0.0;
0153  if (fDropletShapeStr == "sphere")
0154  {
0155   G4cout << "Droplet shape = sphere" << G4endl;
0156   fDropletShape = new G4Orb("dropletSV", fDropletR);
0157   sphericalUncertainty = 0.0;
0158  }
0159  else if (fDropletShapeStr == "halfSphere")
0160  {
0161   G4cout << "Droplet shape = halfSphere" << G4endl;
0162   fDropletShape = new G4Sphere("dropletSV", 0.0, fDropletR,                                 0.0, 180.*deg, 0.0, 180.*deg);
0163   sphericalUncertainty = fDropletR;
0164  }
0165  else if (fDropletShapeStr == "cylinder")
0166  {
0167   G4cout << "Droplet shape = cylinder" << G4endl;
0168   fDropletShape = new G4Tubs("dropletSV", 0, fDropletR/std::sqrt(3), fDropletR/std::sqrt(3), 0, 360.*deg);
0169   sphericalUncertainty = fDropletR*(1-1/std::sqrt(3));
0170  }
0171  else if (fDropletShapeStr == "box")
0172  {
0173   G4cout << "Droplet shape = box" << G4endl;
0174   fDropletShape = new G4Box("dropletSV", fDropletR/std::sqrt(3), fDropletR/std::sqrt(3), fDropletR/std::sqrt(3));
0175   sphericalUncertainty = fDropletR*(1-1/std::sqrt(3));
0176  }
0177  else
0178  {
0179   std::ostringstream message;
0180    message << "Invalid droplet shape = " << fCloudShapeStr << "!";
0181         G4Exception("DetectorConstruction::Construct()", "GeomSolids0002",
0182                     FatalException, message);
0183     }
0184
0185 //
0186 // Materials
0187 //
0188 // Compute the density of air at 14 km using the Barometric formula
0189 // see, e.g., https://en.wikipedia.org/wiki/Density_of_air
0190 G4double h = 14.0*km;
0191 G4double p0 = 101325*hep_pascal;
0192 G4double T0 = 288.15*kelvin;
0193 G4double grav = 9.80665*m/(s*s);
0194 G4double La = 0.0065*kelvin/m;
0195 G4double R = 8.31447*joule/(mole*kelvin);
0196 G4double M = 0.0289644*kg/mole;
0197 G4double T = T0 - La*h;
0198 G4double p = p0*std::pow(1-La*h/T0,grav*M/(R*La));
0199 G4double air_density = p*M/(R*T);
0200
0201 // make materials and set densities
0202 G4Material* air_mat = nist->BuildMaterialWithNewDensity("Atmosphere","G4_AIR",air_density);
0203 G4Material* water_mat = nist->FindOrBuildMaterial("G4_WATER");
0204
0205 G4double water_density = water_mat->GetDensity();
0206 G4double ice_density = 0.9168*g/cm3;
0207
0208 G4Material* ice_mat = new G4Material("Water ice ", ice_density, 1, kStateSolid, T, p);
0209 ice_mat->AddMaterial(water_mat, 1.);
0210
0211 //
0212 // Droplets
0213 //
0214 G4double droplet_density = water_density;
0215 G4Material* droplet_mat = water_mat;
0216
0217 G4double droplet_count = fDropletNumDens*(fCloudShape->GetCubicVolume());
0218
0219 G4double droplet_volume = fDropletShape->GetCubicVolume();
0220 G4double droplet_total_volume = droplet_count*droplet_volume;
0221
0222 G4double droplet_total_mass =  droplet_total_volume*droplet_density;
0223
0224 //
0225 // Cloud macroscopic quantities
0226 //
0227 G4double cloud_volume = fCloudShape->GetCubicVolume();
0228 G4double cloud_air_volume = cloud_volume - droplet_total_volume;
0229 G4double cloud_air_mass = air_density*cloud_air_volume;
0230
0231 //
0232 // Step limit
0233 //
0234 fStepLimits = new G4UserLimits(fStepLim);
0235
0236 //
0237 // Build world
0238 //
0239 G4Box* solidWorld =new G4Box("World",//its name
0240                   0.5*world_sizeXY,//half x-span
0241                   0.5*world_sizeXY,//half y-span
0242                   0.5*world_sizeZ);//half z-span
0243
0244 G4LogicalVolume* logicWorld =
0245         new G4LogicalVolume(solidWorld, //its solid
0246                             air_mat, //its material
0247                             "World");//its name
0248
0249 logicWorld->SetUserLimits(fStepLimits);
0250
0251 G4VPhysicalVolume* physWorld = new G4PVPlacement(nullptr,//no rotation
0252                          G4ThreeVector(),//at (0,0,0)
0253  logicWorld,//its logical volume
0254  "World",//its name
0255 nullptr,//its mothervolume
0256 false,//no boolean operation
0257 0, //copy number
0258 checkOverlaps); //overlaps checking
0259
0260 //
0261 // Build cloud
0262 //
0263 G4LogicalVolume* logicCloud;
0264
0265 // **********************************************************
0266 //
0267 // Build the cloud using the FastAerosol geometry class
0268 //
0269 // ***********************************************************
0270
0271 if (fFastAerosolCloud)
0272 {
0273 G4cout << "\nFastAerosol geometry with n=" << fDropletNumDens*mm3 << "/mm3, r=" << fDropletR/mm << "mm spheres.\n" << G4endl;
0274
0275 fCloud = new FastAerosol("cloud",                                  fCloudShape,//cloud shape
0276 fDropletR, //bounding radius of droplets
0277                                    fMinSpacing,//minimum spacing between droplets
0278                                                            fDropletNumDens, //approximate number of droplets in cloud
0279                                    sphericalUncertainty); //uncertainty in distance to droplet surface from outside using just droplet's origin as info
0280 fCloud->SetDropletsPerVoxel(4);
0281
0282 FastAerosolSolid* solidCloud = new FastAerosolSolid("cloudSV",//its name
0283 fCloud, //its shape
0284 fDropletShape); //its droplets
0285
0286 solidCloud->SetStepLim(fStepLim);
0287 //FastAerosol can use step limit to speed calculations
0288
0289 logicCloud = new G4LogicalVolume(solidCloud,//its solid
0290                              droplet_mat,//its material
0291                 "cloudLV");//its name
0292 logicCloud->SetUserLimits(fStepLimits);
0293 logicCloud->SetVisAttributes(G4VisAttributes(G4Colour(0.0,0.0,1.0,0.4)));
0294
0295 new G4PVPlacement(nullptr, G4ThreeVector(), logicCloud,                               "cloudPV", //its name
0296 logicWorld,//its mother volume
0297 false, //no boolean operation
0298 0,//copy number
0299 checkOverlaps);//overlaps checking
0300
0301 fCloud->SetSeed(fCloudSeed);
0302
0303 // fPrePopulate = whether to populate all voxels at the beginning or on the fly
0304 if (fPrePopulate)
0305 {
0306  // populate (proving it to the user by printing population reports)
0307  clock_t t;
0308  t = clock();
0309
0310  G4cout << "\nBefore populating" << G4endl;
0311  G4cout <<   "=================" << G4endl;
0312  fCloud->PrintPopulationReport();
0313  G4cout << "\nPopulating..." << G4endl;
0314  fCloud->PopulateAllGrids();
0315  G4cout << "\nAfter populating" << G4endl;
0316  G4cout <<   "================" << G4endl;
0317  fCloud->PrintPopulationReport();
0318  G4cout << G4endl;
0319
0320  t = clock() - t;
0321
0322  G4cout << "\nThis took " << ((float)t)/CLOCKS_PER_SEC << "s\n" << G4endl;
0323
0324  // make filename variables to save data
0325  G4String rStr = std::to_string(fDropletR/mm);
0326  rStr.erase ( rStr.find_last_not_of('0') + 1, std::string::npos );  // drop trailing 0
0327  replace( rStr.begin(), rStr.end(), '.', 'p');
0328  if (rStr.back() == 'p') { rStr.pop_back(); }   // don't write "3p" for 3.0, just write "3"
0329
0330  // want to represent the number density as 1E-ApB for some A, B
0331 G4int order10 = (G4int) -round(10*std::log10(fDropletNumDens*mm3)); // gives 10x the exponent rounded to the int (10x so we get two decimals)
0332 G4int leading = order10 / 10;   // first number
0333 G4int trailing = order10 % 10;  // second number
0334 G4String nStr = "1E-" + std::to_string(leading) + "p" + std::to_string(trailing);
0335
0336 // save population time
0337 std::ofstream file;
0338 file.open("popTime_r" + rStr + "mm_n" + nStr + "mm-3.csv");
0339 file << ((float)t)/CLOCKS_PER_SEC;
0340 file.close();
0341
0342 // save distribution
0343 G4String fName = "distribution_r" + rStr + "mm_n" + nStr + "mm-3.csv";
0344 fCloud->SaveToFile(fName);
0345  }
0346 }
0347 // **********************************************************
0348 //
0349 // (For comparision/benchmarking) Build the cloud using G4VParameterized (does not use FastAerosol)
0350 //
0351 // ***********************************************************
0352 // the droplet positions for this cloud are those saved in the "distribution" folder of our data
0353 // this is to make comparable simulations between FastAerosol and parameterised clouds
0354 // this requires that we first simulate FastAerosol (pre-populated) to generate the positions
0355
0356 else if (fParameterisedCloud)
0357  {
0358  G4cout << "\nParameterised geometry with n=" << fDropletNumDens*mm3 << "/mm3 and r=" << fDropletR/mm << "mm spheres.\n" << G4endl;
0359         std::vector<G4ThreeVector> positions;
0360         G4double x,y,z;
0361
0362  // load distribution file
0363  G4String fName;
0364  G4String rStr = std::to_string(fDropletR/mm);
0365         rStr.erase ( rStr.find_last_not_of('0') + 1, std::string::npos );   // drop trailing 0
0366         replace( rStr.begin(), rStr.end(), '.', 'p');
0367         if (rStr.back() == 'p') { rStr.pop_back(); }    // don't write "3p" for 3.0, just write "3"
0368
0369  // want to represent the number density as 1E-ApB for some A, B
0370  G4int order10 = (G4int) -round(10*std::log10(fDropletNumDens*mm3));    // gives 10x the exponent rounded to the int (10x so we get two decimals)
0371  G4int leading = order10 / 10;  // first number
0372  G4int trailing = order10 % 10; // second number
0373  G4String nStr = "1E-" + std::to_string(leading) + "p" + std::to_string(trailing);
0374
0375  fName = "distribution_r" + rStr + "mm_n" + nStr + "mm-3.csv";
0376  std::ifstream infile(fName);
0377  std::string line;
0378
0379  while (getline(infile,line)) {
0380  std::istringstream stream(line);
0381  std::string field;
0382  getline(stream,field,','); x = stod(field)*mm;
0383  getline(stream,field,','); y = stod(field)*mm;
0384  getline(stream,field,','); z = stod(field)*mm;
0385
0386  positions.push_back(G4ThreeVector(x,y,z));
0387         }
0388
0389  G4VPVParameterisation* cloudParam =
0390             new CloudParameterisation(positions);
0391
0392  G4Box* cloudBounding = new G4Box("cloudBounding",//its name
0393                           0.5*cloud_sizeXY,//half x-span
0394                   0.5*cloud_sizeXY, //half y-span
0395                   0.5*cloud_sizeZ); //half z-span
0396
0397  logicCloud = new G4LogicalVolume(cloudBounding,    //its solid
0398                     air_mat,//its material
0399                 "cloudLV");//its name
0400
0401  logicCloud->SetSmartless(fSmartless);
0402  logicCloud->SetUserLimits(fStepLimits);
0403  logicCloud->SetVisAttributes(G4VisAttributes(false));
0404
0405  new G4PVPlacement(nullptr,//no rotation
0406            G4ThreeVector(),//at position
0407            logicCloud,//its logical volume
0408            "cloudPV",//its name
0409            logicWorld,  //its mothervolume
0410            false, //no boolean operation
0411             0,//copy number
0412            checkOverlaps);//overlaps checking
0413
0414  G4LogicalVolume* logicDroplet =
0415             new G4LogicalVolume(fDropletShape,//its solid
0416             droplet_mat,//its material
0417             "dropletLV");//its name
0418
0419  logicDroplet->SetUserLimits(fStepLimits);
0420
0421  new G4PVParameterised("droplets",//its name
0422             logicDroplet,//droplet logical volume
0423             logicCloud,//mother logical volume
0424                 kUndefined,//droplets placed along this axis
0425                 positions.size(),//number of droplets
0426                  cloudParam);//the parametrisation
0427  }
0428 // **********************************************************
0429 //
0430 // (For comparision/benchmarking) Simulate the cloud by smearing droplets out into a single solid (does not use FastAerosol)
0431 //
0432 // ***********************************************************
0433 else if (fSmoothCloud)
0434  {
0435   G4cout << "\nSmooth geometry based on a cloud of n=" << fDropletNumDens*mm3 << "/mm3 and r=" << fDropletR/mm << "mm spheres.\n" << G4endl;
0436 // build cloud by smearing the droplets uniformly across the cloud volume, for comparison/benchmarking purposes (does not use FastAerosol)
0437 G4Material* cloud_mat = new G4Material("Cloud", (droplet_total_mass+cloud_air_mass)/cloud_volume, 2);
0438
0439 cloud_mat->AddMaterial(droplet_mat, droplet_total_mass/(cloud_air_mass+droplet_total_mass));
0440
0441 cloud_mat->AddMaterial(air_mat, cloud_air_mass/(cloud_air_mass+droplet_total_mass));
0442
0443  logicCloud = new G4LogicalVolume(fCloudShape,  //its solid
0444                   cloud_mat,    //its material
0445                   "cloudLV");   //its name
0446
0447  logicCloud->SetUserLimits(fStepLimits);
0448  logicCloud->SetVisAttributes(G4VisAttributes(G4Colour(0.0,0.0,1.0,0.4)));
0449
0450  new G4PVPlacement(nullptr, //no rotation
0451            G4ThreeVector(), //at position
0452            logicCloud,//its logical volume
0453            "cloudPV", //its name
0454             logicWorld, //its mothervolume
0455                 false,  //no boolean operation
0456             0,//copy number
0457             checkOverlaps); //overlaps checking
0458  }
0459  else
0460  {
0461  G4cout << "\nNo cloud.\n" << G4endl;
0462  }
0463
0464  //
0465  // Build detector
0466  //
0467  G4double detector_sizeXY = cloud_sizeXY;
0468  G4double detector_sizeZ = 0.05*m;
0469  G4Material* detector_mat = nist->FindOrBuildMaterial("G4_Al");
0470  G4ThreeVector detector_pos = G4ThreeVector(0, 0, 0.5*1.05*cloud_sizeZ);
0471
0472  G4Box* soldDetector =  new G4Box("detectorSV", //its name
0473                    0.5*detector_sizeXY, //half x-span
0474                    0.5*detector_sizeXY, //half y-span
0475                    0.5*detector_sizeZ); //half z-span
0476
0477  G4LogicalVolume* logicDetector =
0478         new G4LogicalVolume(soldDetector,   //its solid
0479                     detector_mat,   //its material
0480                     "detectorLV");  //its name
0481
0482  logicDetector->SetUserLimits(fStepLimits);
0483
0484  new G4PVPlacement(nullptr,//no rotation
0485            detector_pos,        //at position
0486            logicDetector,       //its logical volume
0487            "detectorPV",        //its name
0488                logicWorld,          //its mothervolume
0489            false,    //no boolean operation
0490             0,  //copy number
0491            checkOverlaps);      //overlaps checking
0492 //
0493 // Scoring Volume
0494 //
0495 fScoringVolume = logicDetector;
0496
0497 return physWorld;
0498 }
0499