Hello @guatelli I tried to follow your advice…
I retrieved the positrons
// Retrieve the secondary particles
G4SteppingManager* steppingManager = fpSteppingManager;
fSecondary = steppingManager -> GetfSecondary();
G4String secondaryParticleName;
G4String process;
G4int fSectrackID;
G4int fSsecparID;
G4String volumeName;
for(size_t lp1=0;lp1<(*fSecondary).size(); lp1++)
{
secondaryParticleName = (*fSecondary)[lp1]->GetDefinition() -> GetParticleName(); // name of the secondary
process = (*fSecondary)[lp1]-> GetCreatorProcess()-> GetProcessName(); // process creating it
fSectrackID = (*fSecondary)[lp1]->GetTrackID(); //track id
fSsecparID = (*fSecondary)[lp1]->GetParentID(); //parent id
volumeName = (*fSecondary)[lp1] -> GetVolume() -> GetName(); // volume where secondary was generated
}
Here I selected the positrons outgoing the regenerator porduced by gamma->e+e- and having the same trackID and ParentID of the fSecondary
if(NextVol && ThisVol->GetName()=="Regenerator" && NextVol->GetName()=="World" && step->GetTrack()->GetParentID() == fSsecparID && step->GetTrack()->GetTrackID() == fSectrackID && step->GetTrack()->GetDynamicParticle()->GetPDGcode()==-11 && process == "GammaConversion" ) {
But I don’t get any positrons outgoing the regenerator…
maybe is the problem
process == "GammaConversion"
I couldn’t find a list to read how G4 calls this process…I wrote the string “GammaConversion” because here Low Energy Electromagnetic Physics - Livermore | geant4.web.cern.ch I found the class G4GammaConversion.
Maybe should the string be different?
Thank you
I think you have to check the generator process in fSecondary. not at the level of the Step
To make sure the name of the process is correct, simply print them in the screen. again, there might be more efficient ways…
faca87
July 7, 2021, 6:05am
10
Hello @guatelli
How can I do that?
I did it…
using the selection
if(NextVol && ThisVol->GetName()=="Regenerator" && NextVol->GetName()=="World" && step->GetTrack()->GetParentID() == fSsecparID && step->GetTrack()->GetTrackID() == fSectrackID && step->GetTrack()->GetDynamicParticle()->GetPDGcode()==-11 ) {
it doesn’t get any positrons (this is the selection where fSectrackID and fSecparID are from the fsecondary as you wrote for method 1)…
using the selection without fsecondary
if(NextVol && ThisVol->GetName()=="Regenerator" && NextVol->GetName()=="World" && step->GetTrack()->GetParentID()!=00 && step->GetTrack()->GetDynamicParticle()->GetPDGcode()==-11 ) {
I get only eBrem and eIon processess…I simulated 10000 primary positrons…it’s strange that I don’t get gamma->e+e- processess…
weller
July 7, 2021, 8:21am
11
do you iterate over all secondaries but store only 1 fSectrackID – the last one? Or does your filtering as posted (the long if-clause) happen within the same scope?
I suggest to take a step back and check what is going on in the simulation.
faca87
July 7, 2021, 10:49am
13
Hello @weller
What do you mean? I didn’t understand good…
Hello @guatelli
What should I do? Anyway, in the simulation I do the scoring mesh of the targets, I get particles (primary positrons, secondary muons, electrons, positrons) outgoing the targets, secondary photons outgoing the world and secondary photons entering the regenerator…
weller
July 7, 2021, 11:47am
14
I have edited the post and marked possible positions. My question was: is the if-statement
if(NextVol && ThisVol->GetName()=="Regenerator" && NextVol->GetName()=="World" && step->GetTrack()->GetParentID() == fSsecparID && step->GetTrack()->GetTrackID() == fSectrackID && step->GetTrack()->GetDynamicParticle()->GetPDGcode()==-11 )
(or a different version of it, does not matter) in the same bracket { } or somewhere else?
faca87
July 7, 2021, 3:58pm
15
Ok thank you @weller …now I understood…
The if is out of the
for(size_t lp1=0;lp1<(*fSecondary).size(); lp1++)
The code is
[...]
// Retrieve the secondary particles
G4SteppingManager* steppingManager = fpSteppingManager;
fSecondary = steppingManager -> GetfSecondary();
G4String secondaryParticleName;
G4String process;
G4int fSectrackID;
G4int fSsecparID;
G4String volumeName;
for(size_t lp1=0;lp1<(*fSecondary).size(); lp1++)
{
secondaryParticleName = (*fSecondary)[lp1]->GetDefinition() -> GetParticleName(); // name of the secondary
process = (*fSecondary)[lp1]-> GetCreatorProcess()-> GetProcessName(); // process creating it
fSectrackID = (*fSecondary)[lp1]->GetTrackID(); //track id
fSsecparID = (*fSecondary)[lp1]->GetParentID(); //parent id
volumeName = (*fSecondary)[lp1] -> GetVolume() -> GetName(); // volume where secondary was generated
}
[...]
if(NextVol && ThisVol->GetName()=="Regenerator" && NextVol->GetName()=="World" && step->GetTrack()->GetParentID() == fSsecparID && step->GetTrack()->GetTrackID() == fSectrackID && step->GetTrack()->GetDynamicParticle()->GetPDGcode()==-11 ) {
G4cout << "Process " << process << G4endl;
int pdg = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
analysisManager->FillNtupleDColumn(4,0, pdg);
double kinEnergy = step->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
analysisManager->FillNtupleDColumn(4,1, kinEnergy);
[..]
You can see the stepping action here
B1SteppingAction.cc (24.7 KB)
weller
July 7, 2021, 4:22pm
16
ok, then besides reasons based on physics (I have not thought about whether (*fSecondary).size()
must be 0 or max. 1), this might be limiting your number of “detected” events:
since fSectrackID and fSsecparID are set in the for-loop, their value gets overwritten by the next iteration of the loop. after the loop has finished, only the last values are used for comparison in your if-statement.
I would debug this with some output commands:
G4cerr << "secTrackID: " << (*fSecondary)[lp1]->GetTrackID() << G4endl;
and so on, to see if there are multiple secTrackIDs and so on, for each event. If things don’t work the way you expect, let the program talk
btw: https://geant4.web.cern.ch/node/283:
1.4. Class data members start with a prefix "f" followed with an upper case
letter. This convention makes easier to understand the code.
1.5. Local variables and functions argument names start with a lower case
letter except for the names starting with known acronyms in capital
case letters.
it is tempting to adapt the leading f in variable names
ok, so the simulation is correct.
faca87
July 8, 2021, 5:00am
18
Hello @weller
Here the output file.
I simulated only 100 events…I know they are low, but you see…the file is 74MB…simulating more events I get a bigger file and it’s a problem to upload it …
I hope they are enough!
Notice that in the file you will read, for example,
G4WT0 > Event n°: 1
G4WT0 > secondaryParticleName: gamma
G4WT0 > process: eBrem
G4WT0 > fSectrackID: 0
G4WT0 > fSsecparID: 792
G4WT0 > volumeName: Regenerator
The number of event (1 in this example) is the value of lp1
variable…
i.e. I got these prints in this way:
for(size_t lp1=0;lp1<(*fSecondary).size(); lp1++)
{
secondaryParticleName = (*fSecondary)[lp1]->GetDefinition() -> GetParticleName(); // name of the secondary
process = (*fSecondary)[lp1]-> GetCreatorProcess()-> GetProcessName(); // process creating it
fSectrackID = (*fSecondary)[lp1]->GetTrackID(); //track id
fSsecparID = (*fSecondary)[lp1]->GetParentID(); //parent id
volumeName = (*fSecondary)[lp1] -> GetVolume() -> GetName(); // volume where secondary was generated
G4cerr << G4endl;
G4cerr << "Event n°: " << lp1 << G4endl;
G4cerr << "secondaryParticleName: " << (*fSecondary)[lp1]->GetDefinition() -> GetParticleName() << G4endl;
G4cerr << "process: " << (*fSecondary)[lp1]-> GetCreatorProcess()-> GetProcessName() << G4endl;
G4cerr << "fSectrackID: " << (*fSecondary)[lp1]->GetTrackID() << G4endl;
G4cerr << "fSsecparID: " << (*fSecondary)[lp1]->GetParentID() << G4endl;
G4cerr << "volumeName: " << (*fSecondary)[lp1] -> GetVolume() -> GetName() << G4endl;
G4cerr << G4endl;
}
Yes…I’ve to add the f to local variables…
Hello @guatelli
It looks like so…but I’ve this problem with the regenerator
weller
July 8, 2021, 6:38am
19
I have looked into your output file, and indeed there are often times multiple secondaries. So my statement from before I believe is correct: only the last secondary defines the values for the variables you check in the if-statement.
would it make sense to simply move the closing ´´}´´ in the B1SteppingAction.cc
that you uploaded earlier from line 123 to line 426?
And what purpose has line 427? maybe this line should be moved to line 91?
faca87
July 8, 2021, 7:12am
20
Hello @weller
I did it…but I still don’t get any event (the hystogram is empty)…and I also don’t get photons entering in the regenerator (instead when the if was out of the for cycle I got photons entering in the regenerator)
B1SteppingAction.cc (25.2 KB)
I tried to move this line…but if I move this line, I don’t get the released energy in the targets…
maybe it checks the first time, it isn’t in the volume1 or volume 2 (maybe it is in the world) and it closes without checking other events…?
B1SteppingAction.cc (25.2 KB)
weller
July 8, 2021, 8:46am
21
in the uploaded nohup.out, fSectrackID is always 0
but then you compare step->GetTrack()->GetTrackID() == fSectrackID
, probably GetTrackID is never 0?
faca87
July 8, 2021, 9:37am
22
Hello @weller , I tried to do these prints:
Printing all e+ produced by gamma->e+e- in the regenerator
for(size_t lp1=0;lp1<(*fSecondary).size(); lp1++){
secondaryParticleName = (*fSecondary)[lp1]->GetDefinition() -> GetParticleName(); // name of the secondary
process = (*fSecondary)[lp1]-> GetCreatorProcess()-> GetProcessName(); // process creating it
fSectrackID = (*fSecondary)[lp1]->GetTrackID(); //track id
fSsecparID = (*fSecondary)[lp1]->GetParentID(); //parent id
volumeName = (*fSecondary)[lp1] -> GetVolume() -> GetName(); // volume where secondary was generated
if (process=="conv" && secondaryParticleName =="e+"){
G4cerr << G4endl;
G4cerr << "Event n°: " << lp1 << G4endl;
G4cerr << "secondaryParticleName: " << (*fSecondary)[lp1]->GetDefinition() -> GetParticleName() << G4endl;
G4cerr << "process: " << (*fSecondary)[lp1]-> GetCreatorProcess()-> GetProcessName() << G4endl;
//G4cerr << "fSectrackID: " << (*fSecondary)[lp1]->GetTrackID() << G4endl;
//G4cerr << "fSsecparID: " << (*fSecondary)[lp1]->GetParentID() << G4endl;
G4cerr << "volumeName: " << (*fSecondary)[lp1] -> GetVolume() -> GetName() << G4endl;
G4cerr << G4endl;
}
and I get several particles…
Printing all e+ produced by gamma->e+e- in the regenerator and outgoing it
for(size_t lp1=0;lp1<(*fSecondary).size(); lp1++){
secondaryParticleName = (*fSecondary)[lp1]->GetDefinition() -> GetParticleName(); // name of the secondary
process = (*fSecondary)[lp1]-> GetCreatorProcess()-> GetProcessName(); // process creating it
fSectrackID = (*fSecondary)[lp1]->GetTrackID(); //track id
fSsecparID = (*fSecondary)[lp1]->GetParentID(); //parent id
volumeName = (*fSecondary)[lp1] -> GetVolume() -> GetName(); // volume where secondary was generated
if (process=="conv" && secondaryParticleName =="e+"){
if(NextVol && ThisVol->GetName()=="Regenerator" && NextVol->GetName()=="World"){
G4cerr << G4endl;
G4cerr << "Event n°: " << lp1 << G4endl;
G4cerr << "secondaryParticleName: " << (*fSecondary)[lp1]->GetDefinition() -> GetParticleName() << G4endl;
G4cerr << "process: " << (*fSecondary)[lp1]-> GetCreatorProcess()-> GetProcessName() << G4endl;
//G4cerr << "fSectrackID: " << (*fSecondary)[lp1]->GetTrackID() << G4endl;
//G4cerr << "fSsecparID: " << (*fSecondary)[lp1]->GetParentID() << G4endl;
G4cerr << "volumeName: " << (*fSecondary)[lp1] -> GetVolume() -> GetName() << G4endl;
G4cerr << G4endl;
}
}
but I don’t get any particle…
it looks like that all positrons outgoing the regenerator are produced by other process (i.e. eBrem)…
but it’s strange…because in the project there is the idea to regenerate e+ from gamma conversion…
weller
July 8, 2021, 9:47am
23
in your nohup.out, the very last entry:
G4WT1 > Event n°: 1
G4WT1 > secondaryParticleName: e+
G4WT1 > process: conv
G4WT1 > fSectrackID: 0
G4WT1 > fSsecparID: 588
G4WT1 > volumeName: Regenerator
is that what you are looking for?
faca87
July 8, 2021, 9:52am
24
Yesh I’m looking for the
process =="conv" && volumeName=="Regenerator"
but
Simulating 10000 events (instead, in the file that I uploaded I simulated only 100 events)
I didn’t set selection on fSecparID and fSectrackID…
weller
July 9, 2021, 5:05pm
25
Correct me if I’m wrong, but isn’t eBrem (Bremsstrahlung) a process where a gamma gets created (by decelleration of a charged particle)?
So maybe only the last process a positron/electron went through before leaving the volume under consideration is eBrem
, but the positron could have originally been created via conv
?
Probably it is highly unlikely for a conv
process to happen very close to the surface (ratio of volume close to and far from the surface?), that’s why you did not encounter at least one of them in your limited number of events…
faca87
July 9, 2021, 5:29pm
26
Hello @weller
Yes, when an Electron/positron moves near atomic nuclei, it can emits photons.
Obviously, it’s possible that a photon produces a positron then it emits photon so I get it as an “eBream”…but how can I know if it was produce by a photon or, for example, by positron of primary beam??
Currently both the beam and secondary particles producded in the targets can reich the regeneraror…
So I can’t know if the parent of the eBream is a photon or other particle…I tried to add a magnetic field in order to deflect charged particles (i.e. to have only photons reaching the regenerator), but even if I don’t get errors, the simulation doesn’t run Add a magnetic field - #4 by faca87