I am new to geant4 and am trying to create a simulation of an integrating sphere with liquid (currently water) and a detector inside it.
The structure of the geometry is almost fully completed. My idea is to fire a certain number of particles (fParticleGun already implemented) from the edge of one of the two conduits and count the number of bounces of the optical photon on the surface of the sphere, the number of photons reaching the yellow detector and the time it takes for the photon to reach that detector.
As far as the two surfaces (sphere and detector) are concerned, I think I have already made them sensitive:
G4MultiFunctionalDetector* SferaIntegratrice_SD = new G4MultiFunctionalDetector("SferaIntegratrice_SD");
G4SDManager::GetSDMpointer()->AddNewDetector(SferaIntegratrice_SD);
//G4VPrimitiveScorer* primitiv1 = new G4PSDoseDeposit("dose");
//SferaIntegratrice_SD->RegisterPrimitive(primitiv1);
SetSensitiveDetector("SferaIntegratriceLV",SferaIntegratrice_SD);
G4MultiFunctionalDetector* SiPM_SD = new G4MultiFunctionalDetector("SiPM_SD");
G4SDManager::GetSDMpointer()->AddNewDetector(SiPM_SD);
//G4VPrimitiveScorer* primitiv2 = new G4PSDoseDeposit("dose");
//SiPM_SD->RegisterPrimitive(primitiv2);
SetSensitiveDetector("SiPMLV",SiPM_SD);
Correct me if I am wrong. To count the number of bounces, for simplicityâs sake letâs consider firing a single optical photon, I think we need to use the class G4PSSphereSurfaceCurrent.hh (reading from the developerâs manual). I couldnât find any examples on it unfortunatelyâŚ
So I ask, how can I count the number of bounces, the number of photons that reach the yellow detector and then disappear, and the time it takes for the photon to reach that detector once it has been fired?
alternatively, you can attach the scorers via commands/macro. there you can also dump results after each run. think you will find suitable examples for that
Edit: sorry, missed that it is a sphere⌠not so sure about that then
but as mentioned, the SphereSurfaceCurrent is not in the list of available primitive scorers for command based scoring, as you will see when you study that page
an idea could be to define an additional hollow sphere, and attach the volumeFlux scorer to it?
I was able to resolve part of my problem in this way:
void SteppingAction::UserSteppingAction(const G4Step* step)
{
// Collect energy and track length step by step
// get volume of the current step
auto volume_pre = step->GetPreStepPoint()->GetPhysicalVolume()->GetName();
auto volume_post = step->GetPostStepPoint()->GetPhysicalVolume()->GetName();
G4double stepLength = 0.;
G4double stepTime = 0.;
if ( step->GetTrack()->GetTrackID() == 1. ) {
stepLength = step->GetStepLength()/CLHEP::cm;
stepTime = step->GetTrack()->GetGlobalTime()/CLHEP::ns;
if ( stepLength == 0 ) {
++n_bounce;
}
if (step->GetTrack()->GetTrackStatus() == fStopAndKill) {
if (volume_post == "SiPMPhys"){
G4cout << "Particle is died after " << stepTime << " [ns], the # of bouce before to exit was " << n_bounce << "\t" << G4endl;
}
}
}
Now i want to save in a .root file the info about stepTime and n_bounce in order to do two different histo.
But i donât understand because if I shot, hard coded, a number of photon > 1, the .root file is not generated:
/* //with this part of code, the simulation doesn't create the .root file
G4int rate = 220;
G4int time = 60;
G4double energy = 45*10^(-6);
G4int n_particle = rate*time*(energy/Energia_fotone);
*/
fParticleGun = new G4ParticleGun(1);
But works with the command: /run/beamOn 10000000. Why this happens?
@samuelemillesoli In C++, ^ is the bitwise XOR operator. So in the energyA expression, the OP is computing `(45*10) ^ (-6) = 0x000001C2 xor 0xFFFFFFFA = 0xFFFFFE38 = -456.
C++ doesnât have an exponential operator. You could (foolishly ) use std::pow(10.,-6) for this, but 45e-6 is the correct notation.
what is the actual value of n_particle during execution (G4cerr << n_particle << G4endl;)? does it work if you put in fParticleGun = new G4ParticleGun(5000);
?
btw: the particle gun that way will produce 5000 identical primaries, at least same origin and same direction (not sure about energy), so if you expect a distribution in any of these properties there you need to be cautious.
also i think it is favorable to have multiple events with few primaries instead of few events with massively many primaries, as each event is independent from the others and can be computed in parallel. /run/beamOn 5000 on the other hand creates 5k events with the number of primaries you specify in the respective action
regarding the gui buttons: you find instruction to append your own button in the documentation, and also in the examples. not sure if it is as straight forward to adjust the default buttons, but please report back if you find out
Running the simulation I realised that the reflectivity parameter seems to have no weight, in fact if I enter 0.01 or 0.99, the number of bounces for the same number of photons sent is the same:
G4OpticalSurface* OpSurface = new G4OpticalSurface("Superfice_Sfera_Integratrice");
G4LogicalBorderSurface* Surface = new
G4LogicalBorderSurface("Superfice_Sfera_Integratrice", physSferaIntegratrice, physCilindro, OpSurface);
OpSurface->SetType(dielectric_dielectric);
OpSurface->SetModel(unified);
OpSurface->SetFinish(groundfrontpainted);
//OpSurface->SetSigmaAlpha(0.1);
#include "Energy.hh"
G4MaterialPropertiesTable* TEFLON_MPT = new G4MaterialPropertiesTable();
std::vector<G4double> energy = {Energia_fotone*eV};
std::vector<G4double> rindex_teflon = {1.30};
std::vector<G4double> transmission = {0.9};
std::vector<G4double> reflectivity = {0.1};
//
//properties of teflon that depend on energy
//
TEFLON_MPT->AddProperty("RINDEX", energy, rindex_teflon);
TEFLON_MPT->AddProperty("TRANSMITTANCE", energy, transmission);
TEFLON_MPT->AddProperty("REFLECTIVITY", energy, reflectivity);
teflon_mat->SetMaterialPropertiesTable(TEFLON_MPT);
OpSurface->SetMaterialPropertiesTable(TEFLON_MPT);
G4MaterialPropertiesTable* WATER_MPT = new G4MaterialPropertiesTable();
//std::vector<G4double> energy1 = {Energia_fotone*eV};
std::vector<G4double> rindex_water = {1.34};
std::vector<G4double> absorption1 = {344.8*cm};
//std::vector<G4double> reflectivity1 = {0.};
//
//properties of water that depend on energy
//
WATER_MPT->AddProperty("RINDEX", energy, rindex_water);
WATER_MPT->AddProperty("ABSLENGTH", energy, absorption1);
//WATER_MPT->AddProperty("REFLECTIVITY", energy1, reflectivity1);
water_mat->SetMaterialPropertiesTable(WATER_MPT);
OpSurface->SetMaterialPropertiesTable(WATER_MPT);
Yet I followed the developerâs guide, following the same pattern in the code described in pag. 257 (listing 5.9).