weller
January 21, 2022, 7:24am
9
faca87:
G4String proc_name = step->GetPreStepPoint()->GetProcessDefinedStep()->GetProcessName();
G4int n_scint = 0;
if(proc_name == "Scintillation"){
G4Scintillation* scint = new G4Scintillation("Scintillation");
n_scint = scint->GetNumPhotons();
analysisManager->FillNtupleDColumn(6,17, n_scint);
}
this is all in the bracket of the if-case
if(NextVol && ThisVol->GetName()=="World" && NextVol->GetName()=="PlasticScintillator") {
when the track crosses a geometric boundary, I don’t think the GetProcessDefinedStep will ever be scintillation, it will always be the one for crossing boundaries.
faca87
January 21, 2022, 7:26am
10
Thank you @weller , then what should I do to compare results of the simulation (i.e. 1.210^7 particles entering the scintillator) with the measurement by the ArduSiPM (i.e. 4.4 10^5 particles)?
What should I do to get scintillation?
As I wrote, I tried to look only in the scintillator i.e.
if(NextVol && ThisVol->GetName()=="PlastiScintillator" && NextVol->GetName()=="PlasticScintillator") {
G4cout << "Process " << proc_name << G4endl;
G4int n_scint = 0;
if(proc_name == "Scintillation"){
G4Scintillation* scint = new G4Scintillation("Scintillation");
n_scint = scint->GetNumPhotons();
analysisManager->FillNtupleDColumn(6,17, n_scint);
}
}
but after some runs it doesn’t make anything…i.e. there aren’t errors messages, but it looks like that the simulation doens’t simulate all the primary particles
Here the example…it simulates 1.2*10^6 events, then stops without errors
log.txt (3.0 KB)
weller
January 21, 2022, 8:17am
11
is this a typo? missing a c there…
faca87:
but after some runs it doesn’t make anything…i.e. there aren’t errors messages, but it looks like that the simulation doens’t simulate all the primary particles
Here the example…it simulates 1.2*10^6 events, then stops without errors
log.txt (3.0 KB)
the log looks ok to me, note that the different workers are processing the events independently. so after the 1.2e6 event start the remaining work load is simply finished.
fyi: Just checked because I was unsure: the process for crossing geometric boundaries is indeed CoupledTransportation
.
faca87
January 21, 2022, 8:35am
12
Yes…I checked the stepping action, there is the “c”
but for other run, I read alll the current runs and the log ends with the root file closing!
Yest it is the crossing geometry CoupledTransportation
Now I tried to write
auto proc_man = step->GetTrack()->GetParticleDefinition()->GetProcessManager();
G4String proc_name = step->GetPreStepPoint()->GetProcessDefinedStep()->GetProcessName();
// G4cout << "Process " << proc_name << G4endl;
G4int n_scint = 0;
if(proc_name == "Scintillation"){
G4Scintillation* scint = new G4Scintillation("Scintillation");
n_scint = scint->GetNumPhotons();
analysisManager->FillNtupleDColumn(8,0, n_scint);
analysisManager->AddNtupleRow(8);
}
out of the “if”… and to save the scintillation photons in a new root ttree…I simulated 3*10^6 electrons, but I don’t get scintillation photons…
I really don’t understand…
weller
January 21, 2022, 8:55am
13
have you compared your code to the OpNovice2 example that was linked here before? mabe you could try to adapt your code to how it was implemented there?
the main difference I see at first glance is that they iterate over a vector
G4ProcessVector* proc_vec = proc_man->GetPostStepProcessVector(typeDoIt);
instead of using
step->GetPreStepPoint()->GetProcessDefinedStep()->GetProcessName();
don’t know if that could solve your issue?
faca87
January 21, 2022, 8:59am
14
yes @weller ! I just did it! I wrote in the stepping action
auto proc_man = step->GetTrack()->GetDynamicParticle()->GetParticleDefinition()->GetProcessManager();
G4ProcessVector* proc_vec = proc_man->GetPostStepProcessVector(typeDoIt);
G4int n_proc = proc_vec->entries();
G4int n_scint = 0;
for(G4int i = 0; i < n_proc; ++i) {
G4String proc_name = (*proc_vec)[i]->GetProcessName();
if(proc_name.compare("Scintillation") == 0) {
auto scint = (G4Scintillation*) (*proc_vec)[i];
n_scint = scint->GetNumPhotons();
}
G4cout << "In this step, " << n_scint << " scintillation photons were produced." << G4endl;
}
Then it would write something about n_scint…i.e. I didn’t write if n_scint>0,…so even if n_scint==0 it should write
In this step, 0 scintillation photons were produced.
but look the log…it didn’t write anything …
log.txt (5.0 KB)
B1SteppingAction.cc (46.8 KB)
weller
January 21, 2022, 9:14am
15
faca87:
In my DetecorConstruction
G4Material* scintillator = nist->FindOrBuildMaterial("G4_PLASTIC_SC_VINYLTOLUENE");
G4Box* PlastScint;
G4Material* ScintMat = scintillator;
PlastScint = new G4Box("PlastScint", 0.5*fScint_sizeX, 0.5*fScint_sizeY, 0.5*fScint_sizeZ);
logicScint = new G4LogicalVolume(PlastScint, ScintMat, "PlasticScintillator");
https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/TrackingAndPhysics/physicsProcess.html#pre-defined-properties
“Starting in version 11.0, some optical material properties are defined in Geant4.”
→ scintillation seems to not be defined in most / all the predefined materials. I think you have to add this yoursel
example: Physics Processes — Book For Application Developers 11.2 documentation
faca87
January 21, 2022, 10:20am
17
Hello @wellr ,
firstly I fixed an issue…it didn’t write the sentence, because I wrote the G4couet in other if …sorry! now I fixed it
B1SteppingAction.cc (46.9 KB)
I added the scintillation property in the detector construction
//Scintillation process
G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
MPT->AddConstProperty("SCINTILLATIONYIELD", 11800./MeV);
scintillator->SetMaterialPropertiesTable(MPT);
B1DetectorConstruction.cc (29.7 KB)
I simulated 3*10^6 primary electrons, but it didn’t get scintillator photons
even if 39 particles entered in the scintillator
PS
This line
G4int n_proc = proc_vec->entries();
produces this warinng
C:\B1\src\B1SteppingAction.cc(784,41): warning C4267: 'inizializzazione': conversione da 'size_t' a 'G4int'. Possibile
perdita di dati. [C:\B1\B1-build\exampleB1.vcxproj]
PPS. to have more particles, I set the scintillator material to the target…then 3*10^6 primary electrons (800MeV) cross a scintillator material…but I still don’t get scintillator photons…
weller
January 21, 2022, 10:34am
18
what physics list do you use in your code?
does it somehow resemble this from the examples:
G4VModularPhysicsList* physicsList = new FTFP_BERT;
physicsList->ReplacePhysics(new G4EmStandardPhysics_option4());
G4OpticalPhysics* opticalPhysics = new G4OpticalPhysics();
physicsList->RegisterPhysics(opticalPhysics);
runManager->SetUserInitialization(physicsList);
the warning about size_t type not being G4int can be ignored.
faca87
January 21, 2022, 10:57am
19
This is the physicslist:
PhysicsList.cc (8.4 KB)
weller:
does it somehow resemble this from the examples:
G4VModularPhysicsList* physicsList = new FTFP_BERT;
physicsList->ReplacePhysics(new G4EmStandardPhysics_option4());
G4OpticalPhysics* opticalPhysics = new G4OpticalPhysics();
physicsList->RegisterPhysics(opticalPhysics);
runManager->SetUserInitialization(physicsList);
Can you upload the file to see where adding the code? in the example
examples\extended\optical\OpNovice2\src
I don’t see the physiclist file…
weller
January 21, 2022, 11:04am
20
it is the few lines that I quoted for you. You find it here:
https://geant4.kek.jp/lxr/source/examples/extended/optical/OpNovice2/OpNovice2.cc#L59
no trace of G4OpticalPhysics in this file. so that could explain why there is no optical photons
faca87
January 21, 2022, 11:29am
21
Thank you…then I modified in the main file
from
// Physics list
G4VModularPhysicsList* physicsList = new QBBC;
physicsList->SetVerboseLevel(1);
runManager->SetUserInitialization(physicsList);
to
G4VModularPhysicsList* physicsList = new FTFP_BERT;
physicsList->ReplacePhysics(new G4EmStandardPhysics_option4());
G4OpticalPhysics* opticalPhysics = new G4OpticalPhysics();
physicsList->RegisterPhysics(opticalPhysics);
physicsList->SetVerboseLevel(1);
runManager->SetUserInitialization(physicsList);
exampleB1.cc (20.4 KB)
but I still don’t get scintillation photon, even trying to set the target with scintillator material…
weller
January 21, 2022, 12:12pm
22
that is indeed weird. I just implemented the steps as described, and it seems to work.
just for completeness, the steps:
Register G4OpticalPhysics to the physics list by adding this line:
physicsList->RegisterPhysics(new G4OpticalPhysics());
Add Scintillation properties to the material in question:
auto* silicium = nistManager->FindOrBuildMaterial("G4_Si");
G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
MPT->AddConstProperty("SCINTILLATIONYIELD", 11800./MeV);
silicium->SetMaterialPropertiesTable(MPT);
Print/Filter out Scintillation events in the stepping action:
if(aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName() == "Scintillation")
G4cerr << "Process: " << aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName() << G4endl;
Rebuild code
Check output:
...
G4WT1 > Process: Scintillation
...
faca87
January 21, 2022, 12:55pm
23
I added it in my main file
// Physics list
G4VModularPhysicsList* physicsList = new FTFP_BERT;
physicsList->ReplacePhysics(new G4EmStandardPhysics_option4());
G4OpticalPhysics* opticalPhysics = new G4OpticalPhysics();
physicsList->RegisterPhysics(new G4OpticalPhysics());
physicsList->RegisterPhysics(opticalPhysics);
physicsList->SetVerboseLevel(1);
runManager->SetUserInitialization(physicsList);
exampleB1.cc (20.5 KB)
it should be my
G4Material* scintillator = nist->FindOrBuildMaterial("G4_PLASTIC_SC_VINYLTOLUENE");
G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
MPT->AddConstProperty("SCINTILLATIONYIELD", 11800./MeV);
scintillator->SetMaterialPropertiesTable(MPT);
in the detector construction
B1DetectorConstruction.cc (29.7 KB)
I added those lines in the stepping action
B1SteppingAction.cc (47.1 KB)
When I run, I get many Process: Scintillation
printed by the line
G4cerr << "Process: " << aStep->GetPostStepPoint()->GetProcessDefinedStep()-
but I don’t get any scintillation photons (that should be no phyisically possible…if there is scintillation light…there are photons)…
if it can help, this is the full simulation code WeTransfer - Send Large Files & Share Photos Online - Up to 2GB Free
weller
January 21, 2022, 1:34pm
24
ok, think I got it:
https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/TrackingAndPhysics/physicsProcess.html?#example-physproc-5
the material presumably needs some more properties being defined.
auto* silicium = nistManager->FindOrBuildMaterial("G4_Si");
std::vector<G4double> energy = {2.034*eV, 3.*eV, 4.136*eV};
std::vector<G4double> rindex = {1.3435, 1.351, 1.3608};
std::vector<G4double> absorption = {344.8*cm, 850.*cm, 1450.0*cm};
G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
// property independent of energy
MPT->AddConstProperty("SCINTILLATIONYIELD", 11800./MeV);
// properties that depend on energy
MPT->AddProperty("RINDEX", energy, rindex);
MPT->AddProperty("ABSLENGTH", energy, absorption);
silicium->SetMaterialPropertiesTable(MPT);
in the stepping action:
static G4ParticleDefinition* opticalphoton =
G4OpticalPhoton::OpticalPhotonDefinition();
const G4ParticleDefinition* particleDef =
aStep->GetTrack()->GetDynamicParticle()->GetParticleDefinition();
if(particleDef == opticalphoton)
G4cerr << "I am a photon :)" << G4endl;
→
...
G4WT0 > Process: Scintillation
G4WT0 > I am a photon :)
G4WT0 > I am a photon :)
G4WT0 > I am a photon :)
...
edit: in a way this makes sense… without a defined refractive index, what should the wavelength be? defaulting to zero would yield division by zero.
faca87
January 21, 2022, 2:48pm
25
Hello @weller , thank you
I tried to add the code
DetectorConstruction
//Scintillation process
std::vector<G4double> energy = {2.034*eV, 3.*eV, 4.136*eV};
std::vector<G4double> rindex = {1.3435, 1.351, 1.3608};
std::vector<G4double> absorption = {344.8*cm, 850.*cm, 1450.0*cm};
G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
MPT->AddProperty("RINDEX", energy, rindex);
MPT->AddProperty("ABSLENGTH", energy, absorption);
MPT->AddConstProperty("SCINTILLATIONYIELD", 11800./MeV);
scintillator->SetMaterialPropertiesTable(MPT);
and steppingaction
> static G4ParticleDefinition* opticalphoton = G4OpticalPhoton::OpticalPhotonDefinition();
> const G4ParticleDefinition* particleDef = track->GetDynamicParticle()->GetParticleDefinition();
> auto proc_man = track->GetDynamicParticle()->GetParticleDefinition()->GetProcessManager();
> G4ProcessVector* proc_vec = proc_man->GetPostStepProcessVector(typeDoIt);
> G4int n_proc = proc_vec->entries();
> G4int n_scint = 0;
> for(G4int i = 0; i < n_proc; ++i) {
> G4String proc_name = (*proc_vec)[i]->GetProcessName();
> if(proc_name.compare("Scintillation") == 0){
> auto scint = (G4Scintillation*) (*proc_vec)[i];
> n_scint = scint->GetNumPhotons();
> analysisManager->FillNtupleDColumn(8,0, n_scint);
> analysisManager->AddNtupleRow(8);
> }
> }
> if( n_scint > 0){
> G4cout << "In this step, "<< n_scint << " scintillation photons were produced." << G4endl;
> }
but I still don’t get scintillation photons
B1SteppingAction.cc (47.2 KB)
B1DetectorConstruction.cc (30.0 KB)
weller
January 21, 2022, 6:15pm
26
what if you switch back to the variant without iterating over the proc_vec?
and what is values is n_proc when there is scintillation?
faca87
January 22, 2022, 7:28am
27
in which way?
and what is values is n_proc when there is scintillation?
[/quote]
I tried to print
G4cout << "n_proc==" <<n_proc<< G4endl;
both with and without the filter
if(step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName() == "Scintillation"){
and it always writes
weller
January 22, 2022, 7:36am
28
faca87:
in which way?
Yes…I checked the stepping action, there is the “c”
but for other run, I read alll the current runs and the log ends with the root file closing!
Yest it is the crossing geometry https://geant4-forum.web.cern.ch/t/coupledtransportation/2363
Now I tried to write
auto proc_man = step->GetTrack()->GetParticleDefinition()->GetProcessManager();
G4String proc_name = step->GetPreStepPoint()->GetProcessDefinedStep()->GetProcessName();
// G4cout << "Process " << proc_name << G4endl;
…
in a post before (Scintillation light - #14 by faca87 ) you mentioned the issue that it did not print anything, as if n_proc was zero. that is weird now that it is always 10…
do you strictly copy paste the code from your program to your posts?
faca87
January 22, 2022, 8:14am
29
Hello @weller , it didn’t print anything because for error, I wrote the printing line in other “if”…
I fixed that issue and this is my stepping action but I still don’t get scintillation photons
B1SteppingAction.cc (48.0 KB)