// // ******************************************************************** // * License and Disclaimer * // * * // * The Geant4 software is copyright of the Copyright Holders of * // * the Geant4 Collaboration. It is provided under the terms and * // * conditions of the Geant4 Software License, included in the file * // * LICENSE and available at http://cern.ch/geant4/license . These * // * include a list of copyright holders. * // * * // * Neither the authors of this software system, nor their employing * // * institutes,nor the agencies providing financial support for this * // * work make any representation or warranty, express or implied, * // * regarding this software system or assume any liability for its * // * use. Please see the license in the file LICENSE and URL above * // * for the full disclaimer and the limitation of liability. * // * * // * This code implementation is the result of the scientific and * // * technical work of the GEANT4 collaboration. * // * By using, copying, modifying or distributing the software (or * // * any work based on the software) you agree to acknowledge its * // * use in resulting scientific publications, and indicate your * // * acceptance of all terms of the Geant4 Software license. * // ******************************************************************** // #include "DetectorConstruction.hh" #include "G4NistManager.hh" #include "G4Box.hh" #include "G4Tubs.hh" #include "G4LogicalVolume.hh" #include "G4PVPlacement.hh" #include "G4RotationMatrix.hh" #include "G4Transform3D.hh" #include "G4MultiFunctionalDetector.hh" #include "G4VPrimitiveScorer.hh" #include "G4VisAttributes.hh" #include "G4PhysicalConstants.hh" #include "G4SystemOfUnits.hh" #include "G4GeometryManager.hh" #include "G4PhysicalVolumeStore.hh" #include "G4LogicalVolumeStore.hh" #include "G4SolidStore.hh" #include "G4RunManager.hh" #include "G4GDMLParser.hh" #include "G4UniformMagField.hh" #include "G4MonopoleFieldSetup.hh" #include "G4GlobalMagFieldMessenger.hh" #include "G4AutoDelete.hh" #include "G4UserLimits.hh" #include "G4Material.hh" #include "G4Element.hh" #include "G4MaterialTable.hh" #include "G4NistManager.hh" #include "G4VisAttributes.hh" #include "G4Region.hh" #include "G4ProductionCuts.hh" #include "G4LogicalVolume.hh" #include "G4LogicalVolumeStore.hh" #include "G4RegionStore.hh" #include "atmosUnits.hh" DetectorConstruction::DetectorConstruction() : G4VUserDetectorConstruction(), /*fCheckOverlaps(false),*/ fMonFieldSetup() { //static const G4double re=6371.2*CLHEP::km; DefineMaterials(); } DetectorConstruction::~DetectorConstruction() { } void DetectorConstruction::DefineMaterials() { G4NistManager* man = G4NistManager::Instance(); G4bool isotopes = false; G4Element* O = man->FindOrBuildElement("O" , isotopes); G4Element* Si = man->FindOrBuildElement("Si", isotopes); G4Element* Lu = man->FindOrBuildElement("Lu", isotopes); G4Material* LSO = new G4Material("Lu2SiO5", 7.4*g/cm3, 3); LSO->AddElement(Lu, 2); LSO->AddElement(Si, 1); LSO->AddElement(O , 5); } G4VPhysicalVolume* DetectorConstruction::Construct() { G4GeometryManager::GetInstance()->OpenGeometry(); G4PhysicalVolumeStore::GetInstance()->Clean(); G4LogicalVolumeStore::GetInstance()->Clean(); G4SolidStore::GetInstance()->Clean(); G4double z, a,density, fractionmass; G4String name, symbol; G4int nel, ncomponents, natoms; a = 1.*g/mole; density = 1.e-10*g/cm3; G4Element* elN = new G4Element( "Nitrogen", "N", 7. , 14.00674*g/mole ); G4Element* elO = new G4Element( "Oxygen", "O", 8. , 15.9994*g/mole ); G4Element* elC = new G4Element(name = "Carbon", symbol = "C", z = 6, a = 12.01*g/mole); G4Element* elH = new G4Element(name = "Hydrogen", symbol = "H", z = 1, a = 1.01*g/mole); G4Element* elHe = new G4Element(name = "Helium", symbol="He", z=2, a =4.0026*g/mole); G4Element* elAr = new G4Element(name = "Argon", symbol = "Ar", z = 18, a = 39.948*g/mole); G4bool checkOverlaps = true; G4Material* Vacuum =new G4Material(name="Vacuum",density, nel=2,kStateSolid, 293.0*kelvin, 1.e-6*atmosphere ); Vacuum->AddElement(elN, .7); Vacuum->AddElement(elO, .3); auto nistManager = G4NistManager::Instance(); nistManager->FindOrBuildMaterial("G4_AIR"); G4Material* Air = G4Material::GetMaterial("G4_AIR"); nistManager->FindOrBuildMaterial("G4_Al"); G4Material* G4Al = G4Material::GetMaterial("G4_Al"); G4Element* elCa = new G4Element(name = "Calcium", symbol = "Ca", z = 20, a = 40.078*g/mole); G4Element* elMg = new G4Element(name = "Magnesium", symbol = "Mg", z = 12, a = 24.305*g/mole); G4Material* CaCO3 = new G4Material(name="CaCO3", density = 2.71*g/cm3, ncomponents = 3); CaCO3 -> AddElement(elO, natoms = 3); CaCO3 -> AddElement(elC, natoms = 1); CaCO3 -> AddElement(elCa, natoms = 1); G4Material* MgCO3 = new G4Material(name="MgCO3", density = 2.96*g/cm3, ncomponents = 3); MgCO3 -> AddElement(elO, natoms = 3); MgCO3 -> AddElement(elC, natoms = 1); MgCO3 -> AddElement(elMg, natoms = 1); G4Material* limestone = new G4Material("Limestone", density = 2.711*g/cm3, ncomponents=2); limestone->AddMaterial(CaCO3, fractionmass=90*perCent); limestone->AddMaterial(MgCO3, fractionmass=10*perCent); //// ICE density = 0.917*g/cm3; G4double temp = 273.15 *kelvin; G4double pressure = 1.0 * atmosphere; G4double meanExcitationEnergy = 75.0 * eV; G4Material* ice = new G4Material("ice", density, 2); ice->AddElement(elH, 2); ice->AddElement(elO, 1); G4MaterialPropertiesTable* iceProp = new G4MaterialPropertiesTable(); iceProp->AddConstProperty("TEMPERATURE", temp); iceProp->AddConstProperty("PRESSURE", pressure); ice->SetMaterialPropertiesTable(iceProp); ice->GetIonisation()->SetMeanExcitationEnergy(meanExcitationEnergy); G4VisAttributes * VisAttEarth = new G4VisAttributes(G4Colour(0.,0.,1.)); VisAttEarth->SetVisibility(true); G4VisAttributes * VisAttWorld = new G4VisAttributes(); VisAttWorld->SetVisibility(false); //// World // G4double MagSizeX=200*re; G4double MagSizeY=200*re; G4double MagSizeZ=200*re; solidWorld= new G4Box("World",MagSizeX/2.,MagSizeY/2.,MagSizeZ/2.); logicWorld= new G4LogicalVolume( solidWorld, Vacuum, "World", 0, 0, 0); maxStepSize=1.*km; theWorldUserLimits=new G4UserLimits(maxStepSize); theWorldUserLimits->SetUserMaxTrackLength(200.*re); logicWorld->SetUserLimits(theWorldUserLimits); logicWorld->SetVisAttributes(VisAttWorld); physiWorld= new G4PVPlacement(0, // no rotation G4ThreeVector(), // at (0,0,0) logicWorld, // its logical volume "World", 0, // its mother volume false, // no boolean operations 0, checkOverlaps); // no field specific to volume G4VisAttributes * VisAttSpace = new G4VisAttributes(); VisAttSpace->SetVisibility(true); G4double radSpace = 20*re; G4Orb* Space = new G4Orb("Space", radSpace); G4LogicalVolume* logSpace = new G4LogicalVolume(Space, Vacuum, "Space"); G4UserLimits* theSpaceUserLimits=new G4UserLimits(1*km);//100*km); logSpace->SetUserLimits(theSpaceUserLimits); theSpaceUserLimits->SetUserMaxTrackLength(200*re); logSpace->SetVisAttributes(VisAttSpace); G4VPhysicalVolume* fPhysSpace = new G4PVPlacement(0, G4ThreeVector(0, 0, 0), logSpace, "Space", logicWorld, false, 0, checkOverlaps); for (G4int i=0; i<18; i++) { upperAirMat[i] = new G4Material(name=upperVolName[i], upperDensityAtmos[i]*g/cm3, ncomponents=5, kStateGas, upperTempAtmos[i]*kelvin); upperAirMat[i]->AddElement(elN, upper_N2_percent[i]*perCent); upperAirMat[i]->AddElement(elO, upper_O2_percent[i]*perCent); upperAirMat[i]->AddElement(elHe, upper_He_percent[i]*perCent); upperAirMat[i]->AddElement(elAr, upper_Ar_percent[i]*perCent); upperAirMat[i]->AddElement(elH, upper_H_percent[i]*perCent); } G4int l = 17; upperAtmos[l] = new G4Orb("upperAtmos", re+(1000.000000)*km); upperLogAtmos[l] = new G4LogicalVolume(upperAtmos[l], upperAirMat[l], "upperAtmos"); upperAtmosUserLimits[l] = new G4UserLimits(100*m); upperLogAtmos[l]->SetUserLimits(upperAtmosUserLimits[l]); upperAtmosUserLimits[l]->SetUserMaxTrackLength(200*re); upperPhysAtmos[l] = new G4PVPlacement(0, G4ThreeVector(0, 0, 0), upperLogAtmos[l], "upperAtmos", logSpace, false, 0, checkOverlaps); for (G4int i = 16; i>=0; i--) { upperAtmos[i] = new G4Orb("upperAtmos", re+(i*50.000000+151.000000)*km); upperLogAtmos[i] = new G4LogicalVolume(upperAtmos[i], upperAirMat[i], "upperAtmos"); upperAtmosUserLimits[i] = new G4UserLimits(100*m); upperLogAtmos[i]->SetUserLimits(upperAtmosUserLimits[i]); upperAtmosUserLimits[i]->SetUserMaxTrackLength(200*re); upperPhysAtmos[i] = new G4PVPlacement(0, G4ThreeVector(0, 0, 0), upperLogAtmos[i], "upperAtmos", upperLogAtmos[i+1], false, 0, checkOverlaps); } for (G4int i=0; i<150; i++) { airMat[i] = new G4Material(name=volName[i], densityAtmos[i]*g/cm3, ncomponents=5, kStateGas, tempAtmos[i]*kelvin); airMat[i]->AddElement(elN, N2_percent[i]*perCent); airMat[i]->AddElement(elO, O2_percent[i]*perCent); airMat[i]->AddElement(elHe, He_percent[i]*perCent); airMat[i]->AddElement(elAr, Ar_percent[i]*perCent); airMat[i]->AddElement(elH, H_percent[i]*perCent); } G4int k = 149; atmos[k] = new G4Orb("atmos", re+(k+1)*km); logAtmos[k] = new G4LogicalVolume(atmos[k], airMat[k], "atmos"); atmosUserLimits[k] = new G4UserLimits(100*m); //logAtmos[k]->SetUserLimits(atmosUserLimits[k]); //atmosUserLimits[k]->SetUserMaxTrackLength(200*re); physAtmos[k] = new G4PVPlacement(0, G4ThreeVector(0, 0, 0), logAtmos[k], "atmos", upperLogAtmos[0], false, 0, checkOverlaps); for (G4int i = 148; i>=0; i--) { atmos[i] = new G4Orb("atmos", re+(i+1)*km); logAtmos[i] = new G4LogicalVolume(atmos[i], airMat[i], "atmos"); atmosUserLimits[i] = new G4UserLimits(100*m); logAtmos[i]->SetUserLimits(atmosUserLimits[i]); //atmosUserLimits[i]->SetUserMaxTrackLength(200*re); physAtmos[i] = new G4PVPlacement(0, G4ThreeVector(0, 0, 0), logAtmos[i], "atmos", logAtmos[i+1], false, 0, checkOverlaps); } G4VisAttributes * VisAttRock = new G4VisAttributes(); VisAttRock->SetVisibility(true); G4Orb* rock = new G4Orb("Rock", re); G4LogicalVolume* logRock = new G4LogicalVolume(rock, limestone, "Rock"); G4UserLimits* theRockUserLimits=new G4UserLimits(100*m);//100*km); logRock->SetUserLimits(theRockUserLimits); logRock->SetUserLimits(theRockUserLimits); logRock->SetVisAttributes(VisAttRock); G4VPhysicalVolume* fPhysRock = new G4PVPlacement(0, G4ThreeVector(0, 0, 0), logRock, "Rock", logAtmos[0], false, 0, checkOverlaps); if (!RemoveEarth){ solidEarth=new G4Orb("Earth",re-(0.2)*km/*+TopOfAtmosphere*/); logicEarth=new G4LogicalVolume(solidEarth,Vacuum,"Earth",0,0,0); theEarthUserLimits=new G4UserLimits(1.*km); logicEarth->SetUserLimits(theEarthUserLimits); theEarthUserLimits->SetUserMaxTrackLength(200.*re); logicEarth->SetVisAttributes(VisAttEarth); physiEarth=new G4PVPlacement(0, // no rotation G4ThreeVector(0.,0.,0.), // at (0,0,0) logicEarth, // its logical volume "Earth", logRock, //logAtmos[0], //iceLV // its mother volume false, // no boolean operations 0, checkOverlaps); } return physiWorld; } void DetectorConstruction::ConstructSDandField() { if (! fMonFieldSetup.Get()) { G4MonopoleFieldSetup* fieldSetup = new G4MonopoleFieldSetup(); G4AutoDelete::Register(fieldSetup); fMonFieldSetup.Put(fieldSetup); } fMonFieldSetup.Get()->ConstructMagField(); } void DetectorConstruction::SetMaxStepSize(G4double step) { fMaxStepSize = step; /*G4int daughters = lWorldLogVol->GetNoDaughters(); for (G4int i=0; iGetDaughter(i); G4cout << daughter[i]->GetName() << G4endl; if(daughter[i]->GetName().at(0)=='B') { daughter[i]->GetLogicalVolume()->SetUserLimits(new G4UserLimits(fMaxStepSize)); } }*/ }