Good afternoon,
I have been trying lately to simulate a proton beam hitting a series of modules containing silicon sensors and scintillators and I’d like to detect all the visible light (optical photons) produced in my World. While I have no trouble detecting gamma rays (photons of energies ~keV), electrons, muons and other particles, I detect 0 optical photons after running my program.
In the executable I have activated the processes:
/gps/particle proton
/gps/direction 0 0 -1
/gps/energy 400 GeV
/process/activate Cerenkov
/process/activate Scintillation
In DetectorConstruction.cc I have constructed the scintillator:
SetMaterialScint(“Scintillator”);
…
G4double wls_Energy[] = {2.00eV,2.87eV,2.90eV,3.47eV};
const G4int wlsnum = sizeof(wls_Energy)/sizeof(G4double);
G4Material* Scintillator =
new G4Material(“Scintillator”, density= 1.032*g/cm3, ncomponents=2);
Scintillator->AddElement(C, natoms=8);
Scintillator->AddElement(H, natoms=8);
G4double rIndexPstyrene[]={ 1.5916, 1.5916, 1.5916, 1.5916};
assert(sizeof(rIndexPstyrene) == sizeof(wls_Energy));
G4double absorption1[]={2.*cm, 2.*cm, 2.*cm, 2.*cm};
assert(sizeof(absorption1) == sizeof(wls_Energy));
G4double scintilFast[]={0.00, 0.00, 1.00, 1.00};
assert(sizeof(scintilFast) == sizeof(wls_Energy));
fMPTPStyrene = new G4MaterialPropertiesTable();
fMPTPStyrene->AddProperty(“RINDEX”,wls_Energy,rIndexPstyrene,wlsnum);
fMPTPStyrene->AddProperty(“ABSLENGTH”,wls_Energy,absorption1,wlsnum);
fMPTPStyrene->AddProperty(“FASTCOMPONENT”,wls_Energy, scintilFast,wlsnum);
fMPTPStyrene->AddConstProperty(“SCINTILLATIONYIELD”,10./keV);
fMPTPStyrene->AddConstProperty(“RESOLUTIONSCALE”,1.0);
fMPTPStyrene->AddConstProperty(“FASTTIMECONSTANT”, 10.*ns);
fMPTPStyrene->AddConstProperty(“YIELDRATIO”, 1.0);
Scintillator->SetMaterialPropertiesTable(fMPTPStyrene);
// Set the Birks Constant for the Polystyrene scintillator
Scintillator->GetIonisation()->SetBirksConstant(0.126*mm/MeV);
…
//Scintillator B 1
G4Box*
sBox1SB = new G4Box(“ScintillatorB1”, //its name
fScintBSizeX/2,fScintBSizeY/2,fScintBSizeZ/2); //its dimensions
fLBox1SB = new G4LogicalVolume(sBox1SB, //its shape
fMaterialScint, //its material
fMaterialScint->GetName()); //its name
G4double pos_x1SB = 0mm;
G4double pos_y1SB = -(13.50mm - 11.00*mm)/2;
G4double pos_z1SB = pos_z1 + fModSizeZ/2 - fScintBSizeZ/2;
fPBox1SB = new G4PVPlacement(0, //rotation
G4ThreeVector(pos_x1SB, pos_y1SB, pos_z1SB),
“ScintillatorB1”, //its name
fLBox1SB, //its logical volume
fPBoxW, //its mother volume
false, //no boolean operation
0); //copy number
[All other copies are defined similarly.]
…
//Surface of the scintillators
G4OpticalSurface* scintWrapB1 = new G4OpticalSurface(“ScintWrapB1”);
new G4LogicalBorderSurface(“ScintWrapB1”, fPBoxW,
fPBox1SB,
scintWrapB1);
scintWrapB1->SetType(dielectric_metal);
scintWrapB1->SetFinish(polished);
scintWrapB1->SetModel(glisur);
scintWrapB1->SetMaterialPropertiesTable(scintWrapProperty);
In PhysicsList.cc I have registered Optical Physics:
// Optical Physics
G4OpticalPhysics* opticalPhysics = new G4OpticalPhysics();
RegisterPhysics( opticalPhysics );
opticalPhysics->SetWLSTimeProfile(“delta”);
opticalPhysics->SetScintillationYieldFactor(1.0);
opticalPhysics->SetScintillationExcitationRatio(0.0);
opticalPhysics->SetMaxNumPhotonsPerStep(100);
opticalPhysics->SetMaxBetaChangePerStep(10.0);
opticalPhysics->SetTrackSecondariesFirst(kCerenkov,true);
opticalPhysics->SetTrackSecondariesFirst(kScintillation,true);
And added the Cerenkov and Scintillation processes separately:
//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…
void PhysicsList::AddCerenkovProcess()
{
G4Cerenkov* theCerenkovProcess = new G4Cerenkov(“Cerenkov”);
theCerenkovProcess->SetTrackSecondariesFirst(true);
G4int MaxNumPhotons = 300;
theCerenkovProcess->SetMaxNumPhotonsPerStep(MaxNumPhotons);
theCerenkovProcess->SetMaxBetaChangePerStep(10.0);
auto particleIterator=GetParticleIterator();
particleIterator->reset();
while ((particleIterator)()){
G4ParticleDefinition particle = particleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager();
if (theCerenkovProcess->IsApplicable(*particle))
{
pmanager->AddProcess(theCerenkovProcess);
}
}
}
//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…
void PhysicsList::AddScintillation()
{
G4Scintillation* fScintillationProcess = new G4Scintillation(“Scintillation”);
fScintillationProcess->SetScintillationYieldFactor(1.);
fScintillationProcess->SetTrackSecondariesFirst(true);
// Use Birks Correction in the Scintillation process
if(G4Threading::IsMasterThread())
{
G4EmSaturation* emSaturation =
G4LossTableManager::Instance()->EmSaturation();
fScintillationProcess->AddSaturation(emSaturation);
}
auto particleIterator=GetParticleIterator();
particleIterator->reset();
while ((particleIterator)()){
G4ParticleDefinition particle = particleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager();
if (fScintillationProcess->IsApplicable(*particle))
{
pmanager->AddProcess(fScintillationProcess);
}
}
}
In StackingAction.cc I do:
if(particleDefinition == G4OpticalPhoton::Definition()) {
fRunaction->AddNbPhotonsTotal();
fEventaction->AddNbPhotonsPerEventTotal();
G4cout
<< "\n Optical photon created with energy: "
<< G4BestUnit(energy, “Energy”)
<< G4endl;
}
If I replace G4OpticalPhoton with G4Gamma or G4Electron or any other type of particle, it just works fine, but with optical photons always 0 particles are printed.
Do you have any suggestions what could I be missing?
Thank you very much in advance!
