Get interaction of a secondary particle

Hello, in my simulation I simulate a 45 GeV positron beam hitting 2 targets of several shapes. Now I would after the multi-target, a block of tungsten and I would study the positron regeneration from secondary photons. I mean

  1. The positron beam hits the multi-target and it produce several particles (electrons, positrons, photons and muons)
  2. After the target there is the tungsten regenerator, so I want to get the positron outgoing the tungsten produced by photons in the tungsten)…
    I can’t simply put the tungsten after the multi-target, otherwise, the tungsten will be hitted also by other particles (Especially by the positron beam)…

Then I’ve to select just the processes produce by photons…how should I do this selection?

I tried in the stepping action

if(NextVol && ThisVol->GetName()=="World" && NextVol->GetName()=="Tungsten" && step->GetTrack()->GetParentID()!=00 && step->GetTrack()->GetDynamicParticle()->GetPDGcode()==22) { 
if(NextVol && ThisVol->GetName()=="Tungsten" && NextVol->GetName()=="World" && step->GetTrack()->GetParentID()!=00 && step->GetTrack()->GetDynamicParticle()->GetPDGcode()==-11) {

but I don’t get any particle

Hallo,

I suggest to analyze the secondary particles, using fSecondary = steppingManager → GetfSecondary();

Please look at the Stepping Action of the radioprotection example.

Hello @guatelli thank you… I followed your example so I wrote:

// Retrieve the secondary particles
  	G4SteppingManager*  steppingManager = fpSteppingManager;
    fSecondary = steppingManager -> GetfSecondary();
    G4String secondaryParticleName;
     for(size_t lp1=0;lp1<(*fSecondary).size(); lp1++)
     { 
     secondaryParticleName =  (*fSecondary)[lp1]->GetDefinition() -> GetParticleName();  // name of the secondary
 }
if (fRegeneratedPositronFlag){
	if(secondaryParticleName == "gamma" ) {
 
	
		if(NextVol && ThisVol->GetName()=="World" && NextVol->GetName()=="Regenerator" && step->GetTrack()->GetParentID()!=00 && step->GetTrack()->GetDynamicParticle()->GetPDGcode()==22 ) { 
			 		    int pdg = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
			 			analysisManager->FillNtupleDColumn(3,0, pdg);
			 			double kinEnergy = step->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
			 			analysisManager->FillNtupleDColumn(3,1, kinEnergy);
			 			double MomDirx = step->GetTrack()->GetMomentumDirection().x();
			 			analysisManager->FillNtupleDColumn(3,2, MomDirx);
			 			double MomDiry = step->GetTrack()->GetMomentumDirection().y();
			 			analysisManager->FillNtupleDColumn(3,3, MomDiry);
			 			double MomDirz = step->GetTrack()->GetMomentumDirection().z();
			 			analysisManager->FillNtupleDColumn(3,4, MomDirz);
			 			double Vertx = step->GetTrack()->GetVertexPosition().x();
			 			analysisManager->FillNtupleDColumn(3,5, Vertx);
						double Verty = step->GetTrack()->GetVertexPosition().y();
						analysisManager->FillNtupleDColumn(3,6, Verty);
						double Vertz = step->GetTrack()->GetVertexPosition().z();
						analysisManager->FillNtupleDColumn(3,7, Vertz);
						double Mom = step->GetTrack()->GetDynamicParticle()->GetTotalMomentum();
						analysisManager->FillNtupleDColumn(3,8, Mom);
						double Dirxsp = step->GetTrack()->GetPosition().x();
			 			analysisManager->FillNtupleDColumn(3,9, Dirxsp); 
			 			double Dirysp = step->GetTrack()->GetPosition().y();
			 			analysisManager->FillNtupleDColumn(3,10, Dirysp); 
			 			double Dirzsp = step->GetTrack()->GetPosition().z();
			 			analysisManager->FillNtupleDColumn(3,11, Dirzsp); 
			 			double Momx = step->GetTrack()->GetDynamicParticle()->GetMomentum().x();
						analysisManager->FillNtupleDColumn(3,12, Momx);
						double Momy = step->GetTrack()->GetDynamicParticle()->GetMomentum().y();
						analysisManager->FillNtupleDColumn(3,13, Momy);
						double Momz = step->GetTrack()->GetDynamicParticle()->GetMomentum().z();
						analysisManager->FillNtupleDColumn(3,14, Momz);
						double xprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().x())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
						analysisManager->FillNtupleDColumn(3,15, xprime);
						double yprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().y())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
						analysisManager->FillNtupleDColumn(3,16, yprime);
						double ParID=step->GetTrack()->GetParentID();
						analysisManager->FillNtupleDColumn(3,17, ParID);
			   			analysisManager->AddNtupleRow(3);  	
		}
	
	if(NextVol && ThisVol->GetName()=="Regenerator" && NextVol->GetName()=="World" && step->GetTrack()->GetParentID()!=00 && step->GetTrack()->GetDynamicParticle()->GetPDGcode()==-11 ) { 
			 		    int pdg = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
			 			analysisManager->FillNtupleDColumn(4,0, pdg);
			 			double kinEnergy = step->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
			 			analysisManager->FillNtupleDColumn(4,1, kinEnergy);
			 			double MomDirx = step->GetTrack()->GetMomentumDirection().x();
			 			analysisManager->FillNtupleDColumn(4,2, MomDirx);
			 			double MomDiry = step->GetTrack()->GetMomentumDirection().y();
			 			analysisManager->FillNtupleDColumn(4,3, MomDiry);
			 			double MomDirz = step->GetTrack()->GetMomentumDirection().z();
			 			analysisManager->FillNtupleDColumn(4,4, MomDirz);
			 			double Vertx = step->GetTrack()->GetVertexPosition().x();
			 			analysisManager->FillNtupleDColumn(4,5, Vertx);
						double Verty = step->GetTrack()->GetVertexPosition().y();
						analysisManager->FillNtupleDColumn(4,6, Verty);
						double Vertz = step->GetTrack()->GetVertexPosition().z();
						analysisManager->FillNtupleDColumn(4,7, Vertz);
						double Mom = step->GetTrack()->GetDynamicParticle()->GetTotalMomentum();
						analysisManager->FillNtupleDColumn(4,8, Mom);
						double Dirxsp = step->GetTrack()->GetPosition().x();
			 			analysisManager->FillNtupleDColumn(4,9, Dirxsp); 
			 			double Dirysp = step->GetTrack()->GetPosition().y();
			 			analysisManager->FillNtupleDColumn(4,10, Dirysp); 
			 			double Dirzsp = step->GetTrack()->GetPosition().z();
			 			analysisManager->FillNtupleDColumn(4,11, Dirzsp); 
			 			double Momx = step->GetTrack()->GetDynamicParticle()->GetMomentum().x();
						analysisManager->FillNtupleDColumn(4,12, Momx);
						double Momy = step->GetTrack()->GetDynamicParticle()->GetMomentum().y();
						analysisManager->FillNtupleDColumn(4,13, Momy);
						double Momz = step->GetTrack()->GetDynamicParticle()->GetMomentum().z();
						analysisManager->FillNtupleDColumn(4,14, Momz);
						double xprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().x())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
						analysisManager->FillNtupleDColumn(4,15, xprime);
						double yprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().y())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
						analysisManager->FillNtupleDColumn(4,16, yprime);
						double ParID=step->GetTrack()->GetParentID();
						analysisManager->FillNtupleDColumn(4,17, ParID);
			   			analysisManager->AddNtupleRow(4);  	
		}
	}
}

It runs but in the ntupla 4 I get the positrons outgoing the regenerator, instead the ntupla 3 is empty (i.e. I don't get photons entering the regeneator)...so I don't think that I'm getting positrons produce by  by gamma ->e+ e?

Here the stepping action

[B1SteppingAction.cc|attachment](upload://s0qPmLZrBtoyQI8NSBXDn7cstBp.cc) (24.0 KB)

ok. Just to clarify.
If you look at the Stepping Action in radioprotection, the secondaries are analysed when the parent of those secondaries is not tracked anymore (absorbed or going outside the world volume). You can retrieve particle type, process, position, volume of generation of the secondary, energy, etc.

If instead you would like to have information about particles exiting a specific volume, then you need to check the PreStep and the PostStep point of a Step.

Hello @guatelli

Ok then it isn’t my case.

I’ve a positrons beam hitting 2 targets and I look for particles outgoing the targets in this way (in this example I’m looking for muons):


if(NextVol && ThisVol->GetName()==TargetExt && NextVol->GetName()=="World" && step->GetTrack()->GetParentID()!=00 && (step->GetTrack()->GetDynamicParticle()->GetPDGcode()==13 || step->GetTrack()->GetDynamicParticle()->GetPDGcode()==-13)) {

then I know how to look for particles outgoing a volume…

Now I added a new volume called “Regenerator”…I want to select positrons produced in the regenerator by gamma ->e+e- process where the gamma have been produced in the two targets.

I can’t simply do something like

if(NextVol && ThisVol->GetName()=="Regenerator" && NextVol->GetName()=="World" && step->GetTrack()->GetParentID()!=00 && (step->GetTrack()->GetDynamicParticle()->GetPDGcode()==13 || step->GetTrack()->GetDynamicParticle()->GetPDGcode()==-11)) {

because, in this way I will get all positrons outgoing the regenetor, not only the positrons produce by gamma ->e+e- process .
Then, how should I do to select only the positrons produced by this process?
Thank you

What I would do is the following.
Method 1: Retrieve all the positrons generated from the process of interest using the fSecondary. I think you can also retrieve the trackID and the parent track ID (please check).

Then, Method 2: For that event, retrieve the positrons that exit your volume of interest and retrieve their trackID and parent ID. Finally compare the trackID and parent trackID from method 1 and 2. when the TrackID and parent track ID are the same, you find your positron.

I am not sure this is the most efficient method but we use it pretty successfully in our lab.

Cheers

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)