Hi all,
I am building a scintillator of a custom made material(structure) and my setup is like this:
- A sheet absorbs neutrons
- Sometimes the Boron in the sheet decays to Li and alpha, when this happens in the next event a Li and alpha are created in a microscopic moduleSphere outside of the scintillator sheet.
- The moduleSphere consists of randomly positioned spheres made out of phosphor material.
- These phosphor spheres are sensitive detectors giving me the absorbed energy from which I can deduce the amount of photons emitted.
I am not using the built-in G4Scintillator because in my situation this will give me a more realistic photon yield and allows me to change and optimize the phosphor particle size. However, I want to repopulate the module with new randomly positioned spheres every N events such that it simulates different places in the sheet.
PROBLEM: I’ve tried this with building a PopulateModule() method in my EventAction and call it every N loops but now there seem to be issues rendering the program unable to do a beamOn N_beam where N_beam > 2N (maybe memory leak or something?), and the program fails due to segmentation fault. However if I execute beamOn N and then again beamOn N, it works.
Below is my PopulateModule method. It uses given mass ratios to calculate the volume ratio of the scintillator/totalVolume. Then it tries to fill the module with as many scintOrbs (phosphorSpheres) as possible to reach this ratio, if not Solvent is added to get to the ratio. Then, the materials of the sheet and the modulesphere are changed according to these ratios.
void B1EventAction::PopulateModuleSphere() {
if (!G4PhysicalVolumeStore::GetInstance()->GetVolume("ModuleSphere")){
G4cout << "No ModuleSphere found" << G4endl;
return;
}
G4VPhysicalVolume* lModuleSphere = G4PhysicalVolumeStore::GetInstance()->GetVolume("ModuleSphere");
lModuleSphere->GetLogicalVolume()->ClearDaughters();
// set Solvent Ratio
G4double SolventN = B1DetectorConstruction::fSolvent_target;
G4double AbsorberN = B1DetectorConstruction::fAbsorberN;
G4double ScintillatorN = B1DetectorConstruction::fScintillatorN;
G4double BinderN = B1DetectorConstruction::fBinderN;
// Get modulesphere properties
G4Orb* solidModuleSphere = dynamic_cast<G4Orb*>(lModuleSphere->GetLogicalVolume()->GetSolid());
G4double ModuleSphere_radius = solidModuleSphere->GetRadius();
// Set grain sizes --> have to be set as input parameters
G4double R_nominal_grainsize_scint = 6*mm/2;
G4double R_nominal_grainsize_abs = 3*mm/2;
// Making absorber orb at origin where particles originate from
auto absOrbSolid = new G4Orb("absOrb", R_nominal_grainsize_abs);
auto absOrbLogic = new G4LogicalVolume(absOrbSolid, B1DetectorConstruction::fAbsorber_mat, "absOrb");
new G4PVPlacement(nullptr,G4ThreeVector (0.,0.,0.),absOrbLogic, "absOrb", lModuleSphere->GetLogicalVolume(), false, 0, false );
// Calculate Volume ratio's
G4double absVolN = AbsorberN/GetDensityInGramsPerCubicMillimeter(B1DetectorConstruction::fAbsorber_mat);
G4double scintVolN = ScintillatorN/GetDensityInGramsPerCubicMillimeter(B1DetectorConstruction::fScintillator_mat);
G4double bindVolN = BinderN/GetDensityInGramsPerCubicMillimeter(B1DetectorConstruction::fBinder_mat);
G4double solvVolN = SolventN/GetDensityInGramsPerCubicMillimeter(B1DetectorConstruction::fSolvent_mat);
G4double totVolN = absVolN+scintVolN+bindVolN+solvVolN;
G4double ScintTargetVolRatio = scintVolN/totVolN;
// Generate positions and radii of spheres
G4double rvar = 0.2; //Radius variation of 20 percent
std::vector<G4ThreeVector> pos; // set position vector
std::vector<G4double> radii; // set radius vector
for (G4int i = 0; i<100; ++i){
// custom function which generates positions and radii of nonoverlapping spheres (tries 100k times)
addSphere(R_nominal_grainsize_abs, ModuleSphere_radius, R_nominal_grainsize_scint, rvar, pos, radii);
}
// Orb placement and cut when outside of modulesphere
G4double scintVol;
G4double volRatio;
G4int npos = pos.size();
G4String message;
G4VSensitiveDetector* SDorb = G4SDManager::GetSDMpointer()->FindSensitiveDetector("ScintOrbSD");
for (G4int k=0; k<npos; ++k){
auto scintOrbSolid = new G4Orb("ScintOrb", radii[k]);
auto intersect = new G4IntersectionSolid("ScintOrb",lModuleSphere->GetLogicalVolume()->GetSolid(), scintOrbSolid,nullptr, pos[k]);
auto scintOrbLogic = new G4LogicalVolume(intersect, B1DetectorConstruction::fScintillator_mat, "ScintOrb");
new G4PVPlacement(nullptr,G4ThreeVector(0,0,0),scintOrbLogic, "ScintOrb", lModuleSphere->GetLogicalVolume(), false, k, false );
// Calculating volume
scintVol = scintVol + intersect->GetCubicVolume();
volRatio = scintVol/lModuleSphere->GetLogicalVolume()->GetSolid()->GetCubicVolume();
if ( volRatio >= ScintTargetVolRatio){
message = "[ModuleSphere] >> Volume ratio of " + std::to_string(ScintTargetVolRatio) + " reached! Current is: " + std::to_string(volRatio);
break;
}
G4LogicalVolumeStore::GetInstance()->GetVolume("ScintOrb")->SetSensitiveDetector(SDorb);
}
if (volRatio < ScintTargetVolRatio){
message = "[ModuleSphere] >> Could not reach desired volume ratio (scint/(scint+abs+bind)) of " + std::to_string(ScintTargetVolRatio) +
", current is " + std::to_string(volRatio) + "\nAdded solvent to reach volumetric ratio, CHECK IF REALISTIC";
// Define New matrix material and add to modulesphere and scintillator sheet
solvVolN = scintVolN/volRatio - (absVolN + scintVolN + bindVolN); // calculated solvent volume to remain in volume
G4Material* newmatrixmat = MakeMatrixMaterial(B1DetectorConstruction::fAbsorberN, B1DetectorConstruction::fBinderN, solvVolN*GetDensityInGramsPerCubicMillimeter(B1DetectorConstruction::fSolvent_mat));
G4LogicalVolumeStore::GetInstance()->GetVolume("Shape3")->SetMaterial(newmatrixmat); // Scintillator sheet
G4LogicalVolumeStore::GetInstance()->GetVolume("ModuleSphere")->SetMaterial(newmatrixmat);
}
// Notify run manager that geometry has changed
G4RunManager::GetRunManager()->GeometryHasBeenModified();
G4cerr << message << G4endl;}
What am I doing wrong? Is there a better way to do this?