Get interaction of a secondary particle

Hello @guatelli I tried to follow your advice…

  1. 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…

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…

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.

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…

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?

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)

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 :stuck_out_tongue:

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 :slight_smile:

ok, so the simulation is correct.

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

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?

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)

in the uploaded nohup.out, fSectrackID is always 0
but then you compare step->GetTrack()->GetTrackID() == fSectrackID, probably GetTrackID is never 0?

Hello @weller, I tried to do these prints:

  1. 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…

    1. 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…

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?

Yesh I’m looking for the
process =="conv" && volumeName=="Regenerator"

but

  1. Simulating 10000 events (instead, in the file that I uploaded I simulated only 100 events)
  2. I didn’t set selection on fSecparID and fSectrackID…

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…

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