And This is all my PhysicsList:
#include "NPPhysics.h"
#include "DetectorConstruction.h"
#include "G4RunManager.hh"
#include "G4SystemOfUnits.hh"
#include "G4UnitsTable.hh"
#include "G4PhysicsConstructorRegistry.hh"
#include "G4PhysicsListHelper.hh"
#include "G4VModularPhysicsList.hh"
#include "G4Gamma.hh"
#include "G4Positron.hh"
#include "G4GenericIon.hh"
#include "G4DNAElectronSolvation.hh"
#include "G4DNAElastic.hh"
#include "G4DNAELSEPAElasticModel.hh"
#include "G4DNAExcitation.hh"
#include "G4DNAEmfietzoglouExcitationModel.hh"
#include "G4DNAIonisation.hh"
#include "G4DNAEmfietzoglouIonisationModel.hh"
#include "G4DNAAttachment.hh"
#include "G4DNAMeltonAttachmentModel.hh"
#include "G4DNAVibExcitation.hh"
#include "G4DNASancheExcitationModel.hh"
#include "G4DNAChargeDecrease.hh"
#include "G4WentzelVIModel.hh"
#include "G4CoulombScattering.hh"
#include "G4eCoulombScatteringModel.hh"
#include "G4BraggIonModel.hh"
#include "G4BraggModel.hh"
#include "G4BetheBlochModel.hh"
#include "G4DNAMillerGreenExcitationModel.hh"
#include "G4DNABornExcitationModel.hh"
// Hydrogen
#include "G4DNAChargeIncrease.hh"
// gamma
#include "G4PhotoElectricEffect.hh"
#include "G4LivermorePhotoElectricModel.hh"
#include "G4ComptonScattering.hh"
#include "G4LivermoreComptonModel.hh"
#include "G4GammaConversion.hh"
#include "G4LivermoreGammaConversionModel.hh"
#include "G4RayleighScattering.hh"
#include "G4LivermoreRayleighModel.hh"
// e+
#include "G4eMultipleScattering.hh"
#include "G4eIonisation.hh"
#include "G4eBremsstrahlung.hh"
#include "G4eplusAnnihilation.hh"
#include "G4LossTableManager.hh"
#include "G4UAtomicDeexcitation.hh"
#include "G4UrbanMscModel.hh"
#include "G4LivermoreIonisationModel.hh"
#include "G4MollerBhabhaModel.hh"
#include "G4LivermoreComptonModel.hh"
#include "G4IonFluctuations.hh"
#include "G4LivermoreRayleighModel.hh"
#include "G4LivermoreBremsstrahlungModel.hh"
#include "G4UniversalFluctuation.hh"
#include "G4EmConfigurator.hh"
#include "G4EmParameters.hh"
#include "G4EmDNAChemistry_option3.hh"
#include "G4ProcessTable.hh"
#include "G4PhysicsConstructorFactory.hh"
#include "G4VPhysicsConstructor.hh"
#include "G4BuilderType.hh"
#include "G4StepLimiter.hh"
G4_DECLARE_PHYSCONSTR_FACTORY(NPPhysics);
NPPhysics::NPPhysics(G4int ver) : G4VPhysicsConstructor("NPPhysics"), verbose(ver) {
G4EmParameters* param = G4EmParameters::Instance();
param->SetDefaults();
param->SetFluo(true);
param->SetAuger(true);
param->SetAugerCascade(true);
param->SetDeexcitationIgnoreCut(true);
param->ActivateDNA();
param->SetVerbose(verbose);
//enum sub-type of electro-magetic in G4BUilderType.hh([0,bUnknown],[1, bTransportation], [2, bElectromagnetic], [3, bEmExtra],[4, bDecay], [5, bHadronElastic], [6, bHadronInelastic], [7,bStopping], [8, bIons])
SetPhysicsType(bElectromagnetic);
}
void NPPhysics::ConstructParticle() {
G4Gamma::Gamma();
G4Electron::Electron();
G4Positron::Positron();
G4Proton::Proton();
G4GenericIon::GenericIonDefinition();
G4DNAGenericIonsManager * genericIonsManager;
genericIonsManager=G4DNAGenericIonsManager::Instance();
genericIonsManager->GetIon("alpha++");
genericIonsManager->GetIon("alpha+");
genericIonsManager->GetIon("helium");
genericIonsManager->GetIon("hydrogen");
}
void NPPhysics::ConstructProcess() {
G4PhysicsListHelper *ph = G4PhysicsListHelper::GetPhysicsListHelper();
auto theParticleIterator = GetParticleIterator();
theParticleIterator->reset();
while( (*theParticleIterator)() ) {
G4ParticleDefinition* particle = theParticleIterator->value();
G4String particleName = particle->GetParticleName();
G4ProcessManager* pm = particle->GetProcessManager();
if (particleName == "e-") {
G4DNAElectronSolvation *solvation = new G4DNAElectronSolvation("e-_G4DNAElectronSolvation");
// use macro file control as: /process/dna/e-SolvationSubType Meesungnoen2002
//candidate: Ritchie1994, Terrisol1990,Kreipl2009,Meesungnoen2002_amorphous,Meesungnoen2002
auto therm = G4DNASolvationModelFactory::GetMacroDefinedModel();
therm->SetHighEnergyLimit(10.*eV);// the limit of Champion Elastic model:7.4 eV , ELSEPA model:10eV,
solvation->SetEmModel(therm);
pm->AddDiscreteProcess(solvation);
G4DNAElastic* theDNAElasticProcess = new G4DNAElastic("e-_G4DNAElastic");
theDNAElasticProcess->SetEmModel(new G4DNAELSEPAElasticModel());
pm->AddDiscreteProcess(theDNAElasticProcess);
G4DNAExcitation* theDNAExcitationProcess = new G4DNAExcitation("e-_G4DNAExcitation");
theDNAExcitationProcess->SetEmModel(new G4DNAEmfietzoglouExcitationModel());
pm->AddDiscreteProcess(theDNAExcitationProcess);
G4DNAIonisation* theDNAIonisationProcess = new G4DNAIonisation("e-_G4DNAIonisation");
theDNAIonisationProcess->SetEmModel(new G4DNAEmfietzoglouIonisationModel());
pm->AddDiscreteProcess(theDNAIonisationProcess);
G4DNAVibExcitation* theDNAVibExcProcess = new G4DNAVibExcitation("e-_G4DNAVibExcitation");
pm->AddDiscreteProcess(theDNAVibExcProcess);
G4eBremsstrahlung* theDNABremProcess = new G4eBremsstrahlung("e-_G4DNABremsstrahlung");
pm->AddDiscreteProcess(theDNABremProcess);
}
else if ( particleName == "proton" ) {
ph->RegisterProcess(new G4DNAElastic("proton_G4DNAElastic"), particle);
ph->RegisterProcess(new G4DNAExcitation("proton_G4DNAExcitation"), particle);
ph->RegisterProcess(new G4DNAIonisation("proton_G4DNAIonisation"), particle);
ph->RegisterProcess(new G4DNAChargeDecrease("proton_G4DNAChargeDecrease"), particle);
}
else if ( particleName == "hydrogen" ) {
ph->RegisterProcess(new G4DNAElastic("hydrogen_G4DNAElastic"), particle);
ph->RegisterProcess(new G4DNAExcitation("hydrogen_G4DNAExcitation"), particle);
ph->RegisterProcess(new G4DNAIonisation("hydrogen_G4DNAIonisation"), particle);
ph->RegisterProcess(new G4DNAChargeIncrease("hydrogen_G4DNAChargeIncrease"), particle);
}
else if ( particleName == "gamma") {
G4PhotoElectricEffect* thePhotoElectricEffect = new G4PhotoElectricEffect();
thePhotoElectricEffect->SetEmModel(new G4LivermorePhotoElectricModel());
pm->AddDiscreteProcess(thePhotoElectricEffect);
G4ComptonScattering* theComptonScattering = new G4ComptonScattering();
theComptonScattering->SetEmModel(new G4LivermoreComptonModel());
pm->AddDiscreteProcess(theComptonScattering);
G4GammaConversion* theGammaConversion = new G4GammaConversion();
theGammaConversion->SetEmModel(new G4LivermoreGammaConversionModel());
pm->AddDiscreteProcess(theGammaConversion);
G4RayleighScattering* theRayleigh = new G4RayleighScattering();
pm->AddDiscreteProcess(theRayleigh);
}
else if ( particleName == "alpha" ) {
ph->RegisterProcess(new G4DNAElastic("alpha_G4DNAElastic"), particle);
ph->RegisterProcess(new G4DNAExcitation("alpha_G4DNAExcitation"), particle);
ph->RegisterProcess(new G4DNAIonisation("alpha_G4DNAIonisation"), particle);
ph->RegisterProcess(new G4DNAChargeDecrease("alpha_G4DNAChargeDecrease"), particle);
}
else if ( particleName == "alpha+" ) {
ph->RegisterProcess(new G4DNAElastic("alpha+_G4DNAElastic"), particle);
ph->RegisterProcess(new G4DNAExcitation("alpha+_G4DNAExcitation"), particle);
ph->RegisterProcess(new G4DNAIonisation("alpha+_G4DNAIonisation"), particle);
ph->RegisterProcess(new G4DNAChargeDecrease("alpha+_G4DNAChargeDecrease"), particle);
ph->RegisterProcess(new G4DNAChargeIncrease("alpha+_G4DNAChargeIncrease"), particle);
}
else if ( particleName == "helium" ) {
ph->RegisterProcess(new G4DNAElastic("helium_G4DNAElastic"), particle);
ph->RegisterProcess(new G4DNAExcitation("helium_G4DNAExcitation"), particle);
ph->RegisterProcess(new G4DNAIonisation("helium_G4DNAIonisation"), particle);
ph->RegisterProcess(new G4DNAChargeIncrease("helium_G4DNAChargeIncrease"), particle);
}
else if (particleName == "GenericIon") {
ph->RegisterProcess(new G4DNAIonisation("GenericIon_G4DNAIonisation"), particle);
}
else if (particleName == "e+") {
G4eMultipleScattering* msc = new G4eMultipleScattering();
msc->SetStepLimitType(fUseDistanceToBoundary);
G4eIonisation* eIoni = new G4eIonisation();
eIoni->SetStepFunction(0.2, 100*um);
ph->RegisterProcess(msc, particle);
ph->RegisterProcess(eIoni, particle);
ph->RegisterProcess(new G4eBremsstrahlung(), particle);
ph->RegisterProcess(new G4eplusAnnihilation(), particle);
}
}
G4EmConfigurator* em_config = G4LossTableManager::Instance()->EmConfigurator();
G4VEmModel* mod;
G4double standEnergyLimit = 9.9*MeV;
G4double massFactor = 1.0079/4.0026;
mod = new G4UrbanMscModel();
//mod = new G4GoudsmitSaundersonMscModel();
em_config->SetExtraEmModel("e-","msc", mod,"NP", 0*eV,100*MeV, new G4UniversalFluctuation());
mod = new G4LivermoreIonisationModel();
em_config->SetExtraEmModel("e-","e-_G4LivermoreIoni", mod,"NP", 0*eV,1.0*MeV, new G4UniversalFluctuation());
mod = new G4LivermoreBremsstrahlungModel();
em_config->SetExtraEmModel("e-","e-_G4LivermoreBrem", mod,"NP", 0*eV,1*GeV , new G4UniversalFluctuation());
mod = new G4UrbanMscModel();
em_config->SetExtraEmModel("gamma", "msc", mod, "NP", 0*eV, 100*MeV, new G4UniversalFluctuation());
mod = new G4LivermorePhotoElectricModel();
em_config->SetExtraEmModel("gamma", "phot", mod, "NP", 0*eV, 100*MeV, new G4UniversalFluctuation());
mod = new G4LivermoreComptonModel();
em_config->SetExtraEmModel("gamma", "compt", mod, "NP", 0*eV, 100*MeV, new G4UniversalFluctuation());
mod = new G4LivermoreGammaConversionModel();
em_config->SetExtraEmModel("gamma", "conv", mod, "NP", 0*eV, 100*MeV, new G4UniversalFluctuation());
mod = new G4WentzelVIModel();
mod->SetActivationLowEnergyLimit(1.*MeV);
em_config->SetExtraEmModel("proton","msc",mod,"NP", 0, 100*TeV);
mod = new G4eCoulombScatteringModel();
mod->SetActivationLowEnergyLimit(1.*MeV);
em_config->SetExtraEmModel("proton","CoulombScat",mod,"NP", 0, 100*TeV);
mod = new G4BraggModel();
mod->SetActivationLowEnergyLimit(standEnergyLimit);
em_config->SetExtraEmModel("proton","hIoni", mod,"NP",0.0,2*MeV, new G4IonFluctuations());
mod = new G4BetheBlochModel();
mod->SetActivationLowEnergyLimit(standEnergyLimit);
em_config->SetExtraEmModel("proton","hIoni",mod,"NP",2*MeV,100*TeV,new G4UniversalFluctuation());
mod = new G4UrbanMscModel();
mod->SetActivationLowEnergyLimit(1.*MeV);
em_config->SetExtraEmModel("alpha","msc",mod,"NP", 0, 100*TeV);
mod = new G4BraggIonModel();
mod->SetActivationLowEnergyLimit(standEnergyLimit);
em_config->SetExtraEmModel("alpha","ionIoni", mod,"NP",0.0,2*MeV/massFactor, new G4IonFluctuations());
mod = new G4BetheBlochModel();
mod->SetActivationLowEnergyLimit(standEnergyLimit);
em_config->SetExtraEmModel("alpha","ionIoni", mod,"NP",2*MeV/massFactor,100*TeV, new G4UniversalFluctuation());
mod = new G4UrbanMscModel();
mod->SetActivationLowEnergyLimit(1.*MeV);
em_config->SetExtraEmModel("alpha+","msc",mod,"NP", 0, 100*TeV);
mod = new G4BraggIonModel();
mod->SetActivationLowEnergyLimit(standEnergyLimit);
em_config->SetExtraEmModel("alpha+","ionIoni", mod,"NP",0.0,2*MeV/massFactor, new G4IonFluctuations());
mod = new G4BetheBlochModel();
mod->SetActivationLowEnergyLimit(standEnergyLimit);
em_config->SetExtraEmModel("alpha+","ionIoni", mod,"NP",2*MeV/massFactor,100*TeV, new G4UniversalFluctuation());
// mod = new G4DNAELSEPAElasticModel();
// em_config->SetExtraEmModel("e-","e-_G4DNAElastic", mod,"NP",10.0*eV,1.*MeV);
// mod = new G4DNAEmfietzoglouIonisationModel();
// em_config->SetExtraEmModel("e-","e-_G4DNAIonisation", mod,"NP",11.*eV,1.*MeV);
// mod = new G4DNAEmfietzoglouExcitationModel();
// em_config->SetExtraEmModel("e-","e-_G4DNAExcitation", mod,"NP",9.*eV,1.*MeV);
// mod = new G4DNAMeltonAttachmentModel();
// em_config->SetExtraEmModel("e-","e-_G4DNAAttachment", mod,"NP",4.*eV,13.*eV);
// mod = new G4DNASancheExcitationModel();
// em_config->SetExtraEmModel("e-","e-_G4DNAVibExcitation", mod,"NP",2.*eV,100.*eV);
// mod = new G4DNAIonElasticModel();
// em_config->SetExtraEmModel("proton","proton_G4DNAElastic", mod,"NP",100*eV,1.*MeV);
// mod = new G4DNARuddIonisationModel();
// em_config->SetExtraEmModel("proton","proton_G4DNAIonisation", mod,"NP",100*eV,0.5*MeV);
// mod = new G4DNABornIonisationModel();
// em_config->SetExtraEmModel("proton","proton_G4DNAIonisation",mod,"NP",0.5*MeV,10*MeV);
// mod = new G4DNAMillerGreenExcitationModel();
// em_config->SetExtraEmModel("proton","proton_G4DNAExcitation", mod,"NP",10*eV,0.5*MeV);
// mod = new G4DNABornExcitationModel();
// em_config->SetExtraEmModel("proton","proton_G4DNAExcitation", mod,"NP",0.5*MeV,10*MeV);
// mod = new G4DNADingfelderChargeDecreaseModel();
// em_config->SetExtraEmModel("proton","proton_G4DNAChargeDecrease", mod,"NP",100*eV,10*MeV);
// mod = new G4DNAIonElasticModel();
// em_config->SetExtraEmModel("alpha","alpha_G4DNAElastic", mod,"NP",100*eV,1.*MeV);
// mod = new G4DNARuddIonisationModel();
// em_config->SetExtraEmModel("alpha","alpha_G4DNAIonisation", mod,"NP",1*keV,10*MeV);
// mod = new G4DNAMillerGreenExcitationModel();
// em_config->SetExtraEmModel("alpha","alpha_G4DNAExcitation", mod,"NP",1*keV,10*MeV);
// mod = new G4DNADingfelderChargeDecreaseModel();
// em_config->SetExtraEmModel("alpha","alpha_G4DNAChargeDecrease", mod,"NP",1*keV,10*MeV);
// mod = new G4DNAIonElasticModel();
// em_config->SetExtraEmModel("hydrogen","hydrogen_G4DNAElastic", mod,"NP",100.*eV,1.*MeV);
// mod = new G4DNARuddIonisationModel();
// em_config->SetExtraEmModel("hydrogen","hydrogen_G4DNAIonisation", mod,"NP",100*eV,10*MeV);
// mod = new G4DNAMillerGreenExcitationModel();
// em_config->SetExtraEmModel("hydrogen","hydrogen_G4DNAExcitation", mod,"NP",10*eV,0.5*MeV);
// mod = new G4DNADingfelderChargeIncreaseModel();
// em_config->SetExtraEmModel("hydrogen","hydrogen_G4DNAChargeIncrease", mod,"NP",100*eV,10*MeV);
// mod = new G4DNAIonElasticModel();
// em_config->SetExtraEmModel("helium","helium_G4DNAElastic", mod,"NP",100*eV,1.*MeV);
// mod = new G4DNARuddIonisationModel();
// em_config->SetExtraEmModel("helium","helium_G4DNAIonisation", mod,"NP",1*keV,10*MeV);
// mod = new G4DNAMillerGreenExcitationModel();
// em_config->SetExtraEmModel("helium","helium_G4DNAExcitation", mod,"NP",1*keV,10*MeV);
// mod = new G4DNADingfelderChargeIncreaseModel();
// em_config->SetExtraEmModel("helium","helium_G4DNAChargeIncrease", mod,"NP",1*keV,10*MeV);
// mod = new G4DNAIonElasticModel();
// em_config->SetExtraEmModel("alpha+","alpha+_G4DNAElastic", mod,"NP",100*eV,1.*MeV);
// mod = new G4DNARuddIonisationModel();
// em_config->SetExtraEmModel("alpha+","alpha+_G4DNAIonisation", mod,"NP",1*keV,10*MeV);
// mod = new G4DNAMillerGreenExcitationModel();
// em_config->SetExtraEmModel("alpha+","alpha+_G4DNAExcitation", mod,"NP",1*keV,10*MeV);
// mod = new G4DNADingfelderChargeIncreaseModel();
// em_config->SetExtraEmModel("alpha+","alpha+_G4DNAChargeIncrease", mod,"NP",1*keV,10*MeV);
// mod = new G4DNADingfelderChargeDecreaseModel();
// em_config->SetExtraEmModel("alpha+","alpha+_G4DNAChargeDecrease", mod,"NP",1*keV,10*MeV);
G4VAtomDeexcitation* de = new G4UAtomicDeexcitation();
G4LossTableManager::Instance()->SetAtomDeexcitation(de);
de->SetFluo(true);
de->SetPIXE(true);
de->SetAuger(true);
de->SetAugerCascade(true);
}
The DetectorConstruction.cpp:
G4VPhysicalVolume* DetectorConstruction::Construct() {
G4NistManager *nist = G4NistManager::Instance();
// cell: Water, density = 1.0 * g /cm3
G4Material *H2O = nist->FindOrBuildMaterial("G4_WATER");
world_mat_ = H2O;
cell_material_ = H2O;
// nucleus material
// G4Material *water= nist->BuildMaterialWithNewDensity("G4_WATER_MODIFIED", "G4_WATER", 1.407 * g/cm3);
nucleus_material_ = H2O;
// NP materials
// W_{18}O_{49}
G4Element *W = nist->FindOrBuildElement("W");
G4Element *O = nist->FindOrBuildElement("O");
G4Material *W18O49 = new G4Material("W18O49", 7.16 * g/cm3, 2);
W18O49->AddElement(W, 18);
W18O49->AddElement(O, 49);
// Cu@W18O49
G4Material *Cu = nist->FindOrBuildMaterial("G4_Cu");
auto *Cu_W18O49 = new G4Material("CuW18O49", 7.75 * g/cm3, 2);
Cu_W18O49->AddMaterial(W18O49, 0.9);
Cu_W18O49->AddMaterial(Cu, 0.1);
// Ir@W18O49
G4Material *Ir = nist->FindOrBuildMaterial("G4_Ir");
G4Material *Ir_W18O49 = new G4Material("IrW18O49", 8.01 * g/cm3, 3);
Ir_W18O49->AddElement(O, 0.3275);
Ir_W18O49->AddElement(W, 0.518);
Ir_W18O49->AddMaterial(Ir, 0.1545);
NP_material_ = nist->FindOrBuildMaterial("G4_Au");
/************************************************
* Material : H_{2}O
* Box: 60 * 60 * 60 um3
* default: Invisible
***************************************************/
const G4double world_half_x = 30 * um, world_half_y = 30 * um, world_half_z = 30 * um;
G4Box *world_solid = new G4Box("World_Solid", world_half_x, world_half_y, world_half_z);
G4LogicalVolume *world_logic = new G4LogicalVolume(world_solid, H2O, "World_Logic");
G4VPhysicalVolume *world_phys = new G4PVPlacement(nullptr, G4ThreeVector(), world_logic, "World_Phys", nullptr, false, 0, false);
world_logic->SetVisAttributes(G4VisAttributes::GetInvisible());
/*****************************************************
* Spherical Cell
* radius: 15 um, material : H_{2}O
* Color: #7A67EE
******************************************************/
G4Orb *cell_solid = new G4Orb("Cell_Solid", cell_radius);
G4LogicalVolume *cell_logic = new G4LogicalVolume(cell_solid, H2O, "Cell_Logic");
new G4PVPlacement(nullptr, G4ThreeVector(), cell_logic, "Cell_Phys", world_logic, false, 0, true);
// G4VisAttributes *cell_vis = new G4VisAttributes(G4Color(0.56470588, 0.79215686, 0.97647059, 0.3));
// cell_vis->SetForceSolid(true);
// cell_logic->SetVisAttributes(cell_vis);
G4double max_step;
/***********************************************************
* Spherical nucleus
* radius: 5.5 * um, material: Water, density= 1.407 * g/cm3
* Color:#D2691E
***********************************************************/
G4Orb *nucleus_solid = new G4Orb("Nucleus_Solid", nucleus_radius);
G4LogicalVolume *nucleus_logic = new G4LogicalVolume(nucleus_solid, nucleus_material_, "Nucleus_Logic");
new G4PVPlacement(nullptr, G4ThreeVector(), nucleus_logic, "Nucleus_Phys", cell_logic, false, 0, true);
fScoringVolume = nucleus_logic;
G4VisAttributes *nucleus_vis = new G4VisAttributes(G4Color(1, 0.733333333, 0.81568627, 0.4));
nucleus_vis->SetForceSolid(true);
nucleus_logic->SetVisAttributes(nucleus_vis);
/*********************************************************
* NP
* radis: 10-15 nm, length = 240 * nm;
* Color: #7FFFD4
**********************************************************/
const G4double NP_length = 240 * nm;
NP_number = 200;
G4double minDiameter = 20 * nm;
G4double maxDiameter = 30 * nm;
std::vector<G4VPhysicalVolume*> NP_Phys;
G4LogicalVolume *NP_logic[NP_number];
max_step = 15 * nm / 2.;
for(size_t i = 0; i < NP_number; ++i) {
G4double NP_diameter = GenerateRandomDiameter(minDiameter, maxDiameter);
G4double NP_radius = NP_diameter / 2.0;
G4Tubs *NP_solid = new G4Tubs("NP_Solid", 0., NP_radius, NP_length / 2., 0., 360 * deg);
G4ThreeVector position;
G4RotationMatrix *rot = new G4RotationMatrix();
CylinderInfo cylinder;
std::vector<CylinderInfo> placed_cylinders;
G4bool isOverlap;
do {
cylinder.radius = NP_radius;
cylinder.length = NP_length;
// make NP all in the cell
position = GenerateRandomPosition(nucleus_radius + NP_length, cell_radius - NP_length);
rot = GenerateRandomEulerAngles();
cylinder.position = position;
cylinder.rotation = (*rot);
isOverlap = false;
if (!isOverlap && octree->CheckOverlap(cylinder)) {
isOverlap = true;
}
} while(isOverlap);
placed_cylinders.emplace_back(cylinder);
octree->Insert(cylinder);
G4String NP_name = "NP_" + std::to_string(i);
NP_logic[i] = new G4LogicalVolume(NP_solid, NP_material_, "NP_Logic" + std::to_string(i));
G4VisAttributes *NP_vis = new G4VisAttributes(G4Color(0.847058882, 0.2627451, 0.08235294, 1));
NP_vis->SetForceSolid(true);
NP_logic[i]->SetVisAttributes(NP_vis);
NP_Phys[i] = new G4PVPlacement(rot, position, NP_logic[i], NP_name, cell_logic, false, i);
if(check_overlap) {
if(NP_Phys[i]->CheckOverlaps(1000, 0, false)) {
G4cout << "Have Overlap" << G4endl;
}
}
}
auto step_limit = new G4UserLimits();
step_limit->SetMaxAllowedStep(max_step);
// /***************************************
// * Set Region for Physics List
// *****************************************/
fRegion = new G4Region("NP");
G4ProductionCuts *cuts = new G4ProductionCuts();
G4double defCut = 0.1 * nm;
cuts->SetProductionCut(defCut, "gamma");
cuts->SetProductionCut(defCut, "e-");
cuts->SetProductionCut(defCut, "e+");
cuts->SetProductionCut(defCut, "proton");
fRegion->SetProductionCuts(cuts);
for(size_t i = 0; i < NP_number; ++i) {
fRegion->AddRootLogicalVolume(NP_logic[i]);
NP_logic[i]->SetUserLimits(new G4UserLimits(max_step, DBL_MAX, DBL_MAX, tracking_cut));
}
world_logic->SetUserLimits(step_limit);
cell_logic->SetUserLimits(step_limit);
nucleus_logic->SetUserLimits(step_limit);
return world_phys;
}
/***************************************************************
*@function: Creater an NP diameter
************************************************************/
G4double DetectorConstruction::GenerateRandomDiameter(G4double minDiameter, G4double maxDiameter) {
return minDiameter + (maxDiameter - minDiameter) * G4UniformRand();
}
/***************************************************************************************
* @function: Create the random position
***********************************************************************************************/
G4ThreeVector DetectorConstruction::GenerateRandomPosition(G4double innerRadius, G4double outerRadius) {
G4double r = G4UniformRand() * (outerRadius - innerRadius) + innerRadius;
G4double cosTheta = G4UniformRand() * 2 - 1;
G4double sinTheta = std::sqrt(1 - cosTheta * cosTheta);
G4double phi = G4UniformRand() * 2 * M_PI;
G4double x = r * sinTheta * std::cos(phi);
G4double y = r * sinTheta * std::sin(phi);
G4double z = r * cosTheta;
return G4ThreeVector(x, y, z);
}
/***************************************************************************************
* @function: Create the random rotation
****************************************************************************************/
G4RotationMatrix* DetectorConstruction::GenerateRandomEulerAngles() {
G4double cosTheta = G4UniformRand() * 2 - 1;
G4double sinTheta = std::sqrt(1 - cosTheta * cosTheta);
G4double phi = G4UniformRand() * 2 * M_PI;
G4double cosPhi = std::cos(phi);
G4double sinPhi = std::sin(phi);
G4RotationMatrix *rotation = new G4RotationMatrix();
rotation->rotateX(sinTheta * cosPhi);
rotation->rotateY(sinTheta * sinPhi);
rotation->rotateZ(cosTheta);
return rotation;
}