2 circular targets

[faca87]

Hi @anna thank you.

1.The root file wasn’t created in the C:/ directory, I moved it in that directory to have a shorter path to write in ROOT.
2. In the B1RunAction I modified the
analysisManager->CreateNtupleDColumn("Output Energy");
to
analysisManager->CreateNtupleDColumn("OutputEnergy");
B1RunAction.cc (7.4 KB)

3.I compiled and runned again the file but it looks like I don’t get error as you can see in the log

log.txt (96.1 KB)

but I stil get a 156b root file and I can’t open it by Root (I read the same error that I showed you in my previous message)

[anna]

Is the file of the same size in the directory that you run your analysis? (to eliminate problems with file movement)
Did you change anything else in the B1 example other than the detector and the particle source (primary generator)?

[faca87]

Hi @anna thank you

1. Yes it’s the same size …156b
2. I just modified what I wrote in the B1eventaction and in B1runaction.

Anyway, here the full B1 file. You can also find the root file in the
B1\B1-build\Release directory

[anna]

Thanks.
What is missing from your log.txt is information about file B1.root being written and closed (identical to the one about it being opened).
Also, you have some printouts in B1RunAction::EndOfRunAction (lines 173-182) that are not printed. It means that run is not correctly finished.

When I run your example, what I get is a "exampleB1" received signal SIGSEGV, Segmentation fault. error at the end. It is weird you do not see it on windows.

To see what happens I’ve run gdb (a debugging tool):

gdb --args ./example1 ../run1.mac
(gdb) r
(gdb) bt 2

And what it tells is where exactly is the problem:

0x00007ffff14237cf in G4LogicalVolume::GetMass (this=0x0, forced=false, propagate=true, parMaterial=0x0)
    at /builds/geant4-dev/source/geometry/management/src/G4LogicalVolume.cc:558
558	  if ( (G4MT_mass) && (!forced) )  { return G4MT_mass; }
(gdb) bt 2
#0  0x00007ffff14237cf in G4LogicalVolume::GetMass (this=0x0, forced=false, propagate=true, parMaterial=0x0)
    at /builds/geant4-dev/source/geometry/management/src/G4LogicalVolume.cc:558
#1  0x0000555555560a0f in B1RunAction::EndOfRunAction (this=0x7fffe4088cc0, run=0x7fffe4371280) at /tmp/B1/B1RunAction.cc:140

Before you close and save your file you try to get mass from fScoringVolume in your detector. As you modified this part, it is currently a nullptr, it is not being set. So there is one thing about fixing the file (run action), but then another about fixing your detector. I’d propose you do it in steps:

1. Comment out the part in the run action you do not use (you only need what you’ve introduced). Now the file should be created, but as you’ll see it will be filled with zeros.

2. Now fix the energy that is saved. For this you will need to re-introduce fScoringVolume in your Detector construction (e.g. point to logicEnv), and also uncomment the SteppingAction. It is here that for each step volume is checked and if the step happens to be within fScoringVolume - the energy is added to fEdep in your event action. Now you should see in the file non-zero energy.

3. You have 2 targets so you probably need 2 ntuple columns, fEdep2 and AddEnergy2(double) method in the event action, another check in the stepping action, and fScoringVolume2 = logicEnv2 in your detector, etc. Or you can skip some of it and accumulate energy in the same ntuple column, depends what you need.

[faca87]

Hi @anna

1. I started following the first step.

a. I commented the B1runaction file in this way

B1RunAction.cc (7.4 KB)

b. I compiled the file and I didn’t get errors, but I’m not sure do not get error during the running because at the end of the run I read

    Graphics systems deleted.
    Visualization Manager deleting...
    G4 kernel has come to Quit state.
    UserDetectorConstruction deleted.
    UserPhysicsList deleted.
    UserActionInitialization deleted.
    UserRunAction deleted.
    UserPrimaryGenerator deleted.
    RunManager is deleting RunManagerKernel.
    EventManager deleted.
    Units table cleared.
    TransportationManager deleted.
    Total navigation history collections cleaned: 8
 Deleting memory pools 
    Pool ID 'class G4NavigationLevelRep', size : 0.00961 MB
    Pool ID 'class G4ReferenceCountedHandle<void>', size : 0.000961 MB
    Pool ID 'class G4Event', size : 0.000961 MB
    Pool ID 'class G4PrimaryVertex', size : 0.000961 MB
    Pool ID 'class G4PrimaryParticle', size : 0.000961 MB
    Pool ID 'class G4DynamicParticle', size : 0.00288 MB
    Pool ID 'class G4Track', size : 0.00577 MB
    Pool ID 'class G4TouchableHistory', size : 0.000961 MB
    Pool ID 'class G4CountedObject<void>', size : 0.000961 MB
    Number of memory pools allocated: 9; of which, static: 0
    Dynamic pools deleted: 9 / Total memory freed: 0.024 MB
  
    G4Allocator objects are deleted.
    UImanager deleted.
    StateManager deleted.
    RunManagerKernel is deleted. Good bye :)

This is the full log
log.txt (97.5 KB)

c. I got a 31.1kB Root File and I tried to open it

and this is the plot I got

then I think it worked because I got a histogram filled with zero energies (as you said)…

2. Second step
   a. In the Detectorconstrucion.cc I uncommented the lines

   // Shape 2
   //
   G4Material* shape2_mat = nist->FindOrBuildMaterial(“G4_BONE_COMPACT_ICRU”);
   G4ThreeVector pos2 = G4ThreeVector(0, -1cm, 7cm);

   // Trapezoid shape
   G4double shape2_dxa = 12cm, shape2_dxb = 12cm;
   G4double shape2_dya = 10cm, shape2_dyb = 16cm;
   G4double shape2_dz = 6cm;
   G4Trd solidShape2 =
   new G4Trd(“Shape2”, //its name
   0.5shape2_dxa, 0.5shape2_dxb,
   0.5shape2_dya, 0.5shape2_dyb, 0.5*shape2_dz); //its size

   G4LogicalVolume* logicShape2 =
   new G4LogicalVolume(solidShape2, //its solid
   shape2_mat, //its material
   “Shape2”); //its name

   new G4PVPlacement(0, //no rotation
   pos2, //at position
   logicShape2, //its logical volume
   “Shape2”, //its name
   logicEnv, //its mother volume
   false, //no boolean operation
   0, //copy number
   checkOverlaps); //overlaps checking

   // Set Shape2 as scoring volume
   //

   fScoringVolume = logicShape2;

so that this is my new Detectorconstruction file

B1DetectorConstruction.cc (9.7 KB)

I uncommented the B1SteppingAction file, then this is the file

B1SteppingAction.cc (3.2 KB)

Then I compiled and runned again and this is the new full log, but also this time, I don’t know if the message

Graphics systems deleted.
Visualization Manager deleting...
G4 kernel has come to Quit state.
UserDetectorConstruction deleted.
UserPhysicsList deleted.
UserActionInitialization deleted.
UserRunAction deleted.
UserPrimaryGenerator deleted.
RunManager is deleting RunManagerKernel.
EventManager deleted.
Units table cleared.
TransportationManager deleted.
Total navigation history collections cleaned: 9
 Deleting memory pools
Pool ID 'class G4NavigationLevelRep', size : 0.0115 MB
Pool ID 'class G4ReferenceCountedHandle<void>', size : 0.000961 MB
Pool ID 'class G4Event', size : 0.000961 MB
Pool ID 'class G4PrimaryVertex', size : 0.000961 MB
Pool ID 'class G4PrimaryParticle', size : 0.000961 MB
Pool ID 'class G4DynamicParticle', size : 0.00288 MB
Pool ID 'class G4Track', size : 0.00481 MB
Pool ID 'class G4TouchableHistory', size : 0.000961 MB
Pool ID 'class G4CountedObject<void>', size : 0.000961 MB
Number of memory pools allocated: 9; of which, static: 0
Dynamic pools deleted: 9 / Total memory freed: 0.025 MB
G4Allocator objects are deleted.
UImanager deleted.
StateManager deleted.
RunManagerKernel is deleted. Good bye :) 

is an error! Anyway, I got a new B1.root file I tried to open it

but I got again a zero energy histogram

Why do I still get zero energy?

PS. It’s a good idea to get the enery the energy at the end of the first target and later after the second one, butI’m not sure to understand what I should do.

This is the B1 file You can also find the root file in the
B1\B1-build\Release directory

[anna]

You do not need to uncomment the full DetectorConstruction (step 2). You only need to set fScoringVolume to one of your G4LogicalVolumes. What you did creates the detector volumes from the example (that most likely are not on the path of the particle hence 0s in the histogram).

[faca87]

Hi @anna thank you…I commented again the two shapes and I just wrote:

fScoringVolume = logicEnv2;

B1DetectorConstruction.cc (9.8 KB)

I compiled and runned the file
log.txt (97.3 KB)

but I still get zero energy…

[anna]

Did you uncomment B1SteppingAction also in B1ActionInitialization? It is where it is attached to the list of actions to be invoked in your simulation.

[faca87]

Hi @anna thank you!
I uncommented the line

  SetUserAction(new B1SteppingAction(eventAction));

in B1ActionInitialization

B1ActionInitialization.cc (2.9 KB)

and now I get no-zero energy!

Given that I setted

fScoringVolume = logicEnv2;

in the B1DetectorConstruction, I guess this is the energy after the second target…true?

What shoud I do to get both the released energies after the first and the second targets?

[anna]

You need to create analogous fScoringVolume2 with pointer to logicEnv in detector construction. It means you also need to add it to Stepping Action and check if particle is within volume 1 or volume 2. And as you want to store it separately - invoke in both cases two different methods from event action, the first one adding to fEdep, the second adding to a (to be introduced) fEdep2. Then, you need another ntuple to be created in a run action and to fill it from the end of the event action.

[faca87]

Hi @anna thank you. Following your message

1. In B1detectoconstruction.cc

a. I defined fScoringVolume2(0)

b. I setted

 fScoringVolume = logicEnv;
  fScoringVolume2 = logicEnv2;

B1DetectorConstruction.cc (10.0 KB)

2. In B1SteppingAction.cc
   a. I declared

    fEventAction2(eventAction),
      fScoringVolume2(0)
   

but I’m not sure if I’ve to set fEventAction2(eventAction2),

b. In the B1SteppingAction::UserSteppingAction(const G4Step* step) function I added

 if (!fScoringVolume2) { 
const B1DetectorConstruction* detectorConstruction2
  = static_cast<const B1DetectorConstruction*>
    (G4RunManager::GetRunManager()->GetUserDetectorConstruction());
fScoringVolume2 = detectorConstruction2->GetScoringVolume();   
  }

  // get volume of the current step
  G4LogicalVolume* volume2 
= step->GetPreStepPoint()->GetTouchableHandle()
  ->GetVolume()->GetLogicalVolume();
  
  // check if we are in scoring volume
  if (volume2 != fScoringVolume2) return;

  // collect energy deposited in this step
  G4double edepStep2 = step->GetTotalEnergyDeposit();
  fEventAction2->AddEdep(edepStep2);  

but I’m not sure if I’ve to define step2

B1SteppingAction.cc (4.2 KB)

2. In B1EventAction.cc

a. I defined fEdep2Tar(0.) and fRunAction2(runAction),

but I’m not sure if I had to write fRunAction2(runAction2),

b. In B1EventAction::BeginOfEventAction(const G4Event*) function I setted
fEdep2Tar = 0.;

c. in B1EventAction::EndOfEventAction(const G4Event*) function I setted

fRunAction2->AddEdep(fEdep2Tar);

and

analysisManager->FillNtupleDColumn(1, fEdep2Tar);

B1EventAction.cc (3.9 KB)

3. In B1RunAction
   a. I declared fEdep2Tar(0.)
   b. In B1RunAction::B1RunAction() I added the line
   accumulableManager->RegisterAccumulable(fEdep2Tar);
   c. In B1RunAction::B1RunAction() I added the line
   analysisManager->CreateNtupleDColumn("OutputEnergy2");

d. I modified the B1RunAction::AddEdep(G4double edep) function to B1RunAction::AddEdep(G4double edep, G4double edep2) and I added in it the line

fEdep2Tar += edep2;

B1RunAction.cc (8.0 KB)

Lastly I tried to compile it but

1. I get so much errors, here the log

log.txt (3.8 KB)

2. As I said previous, I’m not sure if I had to define step2 and frunction2.

Here the full file

[anna]

1. Make sure you’ve added fScoringVolume2 to B1DetectorConstruction.hh
   and that you have a method to retrieve it.

2. a) There should be only one event action.
   b) But - it should have 2 methods to add energy.
   fEventAction2->AddEdep(edepStep2); -> fEventAction->AddEdep2(edepStep);

Also the detector should have 2 methods - to retrieve the first and the second scoring volume (method GetScoringVolume2()):

fScoringVolume2 = detectorConstruction2->GetScoringVolume2();   

Please notice that you do not need to re-define energy deposit or the volume as they are the same as the original code, you just give new names (edepStep2 and volume2) and be careful with return statements to avoid the code to exit prematurely (if (volume != fScoringVolume2 && volume != fScoringVolume) return;)

2. (second one) there should be only one run action,
   c) you do not use the energy in the run action, just skip fRunAction2->AddEdep(fEdep2Tar);
   Remember to add method AddEdep2(...) so it can be used in the stepping action.

3. a)b)d) - not needed, you do not use that energy here, it is sufficient to have it in the event action (in the example it is used also in the run action, but it is not needed to save it into the file).

And an overall comment - remember to declare class members and methods in the header files as well (identificatore non dichiarato)

[faca87]

Hi @anna thank you

Hi @anna thank you. So

1. In B1detectorconstruction.hh file , in B1DetectorConstruction : public G4VUserDetectorConstruction function I added the line

G4LogicalVolume* fScoringVolume2;

B1DetectorConstruction.hh (2.5 KB)

I retrieve it, because in B1detectorconstruction.cc there is the line

#include "B1DetectorConstruction.hh"

2. In B1SteppingAction.cc file

a. I deleted the line fEventAction2(eventAction),
b. I modified the line fEventAction2->AddEdep(edepStep2); to
fEventAction2->AddEdep2(edepStep);
c. I modified the last lines writing

  if (!fScoringVolume2) { 
    const B1DetectorConstruction* detectorConstruction2
      = static_cast<const B1DetectorConstruction*>
        (G4RunManager::GetRunManager()->GetUserDetectorConstruction());
    fScoringVolume2 = detectorConstruction2->GetScoringVolume2();   
  }

  // get volume of the current step
  G4LogicalVolume* volume2 
    = step->GetPreStepPoint()->GetTouchableHandle()
      ->GetVolume()->GetLogicalVolume();
      
  // check if we are in scoring volume

if (volume != fScoringVolume2 && volume != fScoringVolume) return;

  // collect energy deposited in this step
  G4double edepStep2 = step->GetTotalEnergyDeposit();
  fEventAction->AddEdep2(edepStep);  

but I’m not sure to understand what you mean in this line

if (volume != fScoringVolume2 && volume != fScoringVolume) return;

because in this way I think we just check about volume, but not about volume2 …or not?

B1SteppingAction.cc (4.4 KB)

3. In B1eventaction.cc file
   a. I deleted the lines

   fRunAction2(runAction),
   fEdep2Tar(0.)

b. I modified the line fRunAction2->AddEdep(fEdep2Tar); to fRunAction->AddEdep2(fEdep);

B1EventAction.cc (3.7 KB)

3. In B1Runaction.cc file
   a. I deleted the line
   fEdep2Tar(0.)
   b. I deleted the line accumulableManager->RegisterAccumulable(fEdep2Tar);
   c. I modified from B1RunAction::AddEdep(G4double edep, G4double edep2) to B1RunAction::AddEdep(G4double edep)
   d. I deleted the line fEdep2Tar += edep2;
   e. I added the lines
   void B1RunAction::AddEdep2(G4double edep)
   {
   fEdep += edep;

   }

but I’m not sure if I had to write fEdep2Tar += edep; instead of fEdep += edep;

B1RunAction.cc (7.7 KB)

4. In B1eventAction.hh file
   a. I added the line
   void AddEdep2(G4double edep) { fEdep2Tar += edep; }
   b. I added the line
   G4double fEdep2tAR;
   B1EventAction.hh (2.7 KB)

5. In B1SteppingAction.hh
   a. I added the line G4LogicalVolume* fScoringVolume2;

B1SteppingAction.hh (2.5 KB)

6. In B1runaction.hh file
   a. I added the line
   void AddEdep2 (G4double edep);

Notice that I also added the line G4Accumulable<G4double> fEdep2Tar; but adding that line I got errors, so I commented it, but I’m not sure if I need it

B1RunAction.hh (2.8 KB)

7. In B1DetectorConstruction.hh file
   a. I added the line
   G4LogicalVolume* GetScoringVolume2() const { return fScoringVolume2; }
   B1DetectorConstruction.hh (2.8 KB)

Lastly, I tried to run and compile. It looks like I didnt’ get errors during running and compiling. This is the full log
log.txt (97.7 KB)

I got a 31,3kB Root File, so I tried to open it and it looks like I have OutputEnergy and OtputEnergy2 in the B1 TTree

This is the full Root log

and here the 2 plots

OutputEnergy

and OutputEnergy2

you see that the 2 plots are identical, then I guess that I’m plotting the same energy, not the first one after the first target and the second one after the second target.

Here the full B1 file. You can also find the root file in B1\B1-build\Release directory
I don’t know what I did wrong…

[anna]

You need to rethink what do you put inside the stepping action, and the order in which you do it. At the moment, you:

1. retrieve the first target from the detector construction,
2. retrieve the current volume in which the step is made,
3. if current volume is not the first target - you exit
4. you add twice (to both histograms) energy deposit in the event action

Number 3) is is why you never see the second target (that’s what I meant by the change to the return statement).

What you need is 1) and 2) and also here to retrieve the second target from the detector construction.
And only then you do comparisons and actions:
a) if current volume == the first target then eventAction->AddEdep(…)
b) else if current volume == the second target then eventAction->AddEdep2(…)
c) else return (if step is outside of both volumes) - optionally, it’s anyway the end

[faca87]

I @anna thank you but I solved the problem yesterday!
This is the new B1steppingaction
B1SteppingAction.cc (4.6 KB)
So I get two different energies…
Now I must add a 2 new culomns

1. The particle id
2. the output angles
   Is it possible?

[anna]

Yes, have a look what can you access from G4Step. E.g. G4Step::GetTrack() will give you access to G4Track (momentum, direction, …), and further from this track you can access G4Track::GetParticleDefinition() to get information on the particle (like PDG code, etc.)

[faca87]

Hi @anna thank you. Are there any examples to have a look to? I would like go forward step by step (then first adding the particle id information)

[anna]

By particle id do you mean e.g. the PDG code?
I’m sure it is in one of the examples, but I do not recall it. All the intrinsic properties of particle are accessible from G4DynamicParticle (so one can see what methods it has). For instance, you can do:

int pdg = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();

and later store pdg in your ntuple, or on a histogram.

[faca87]

Hi @anna thank you. By particle ID I mean the t number related to the particle (i.e. number=particle)! I must store it after the first and the second target.

Currently the simulation is able to describe the released energy in the two targets. But I nee these information too

1. ParticleID: Which particle relased that energy? A primary positron? a secondary muon? A secondary electron?
2. Particle Kinetic energy: How much is the kinetic energy of the particle that released the energy?
3. Particle Direction: What is the direction of particle that released the energy?

** PARTICLE ID**
Following your message I did:

1. In B1SteppingAction.cc
   a. I added
   #include "g4analysis.hh"

b.I added
auto analysisManager = G4AnalysisManager::Instance();

c. added
if (volume == fScoringVolume) {
int pdg = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
analysisManager->FillNtupleDColumn(2, pdg);

	}
	if (volume == fScoringVolume2) {
		int pdg2 = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
		analysisManager->FillNtupleDColumn(3, pdg2);
	}

B1SteppingAction.cc (5.3 KB)

2. In B1RunAction.cc
   a. I wrote
   analysisManager->CreateNtupleDColumn(“ParticleID”);
   analysisManager->CreateNtupleDColumn(“ParticleID2”);

The I compiled and runned without errors

logB1.txt (98.1 KB)

and I got a B1.root file

So I opened the root file

and I get two graphs

the first one for the target 1

and the second one for the target 2

but the 2 plots are identical…thereare two partice having ID=-11 ( I think they are positrons) and one particle having ID=11 (I think it’s an electron).
Then…I’m not sure if the graph are identical because I really have the same particles after each target or because I made mistakes.…

PARTICLE KINETIC ENERGY

1. In B1SteppingAction.cc
   a. I added
   #include "G4Track.hh"
   b. I added
   if (volume == fScoringVolume) {
   double kinEnergy = step->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
   analysisManager->FillNtupleDColumn(2, kinEnergy);

    }
    if (volume == fScoringVolume2) {
    	 double kinEnergy2 = step->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
    	 analysisManager->FillNtupleDColumn(3, kinEnergy2);
    }
   

B1PrimaryGeneratorAction.cc (4.3 KB)

2. In B1RunAction.cc
   a. I added
   analysisManager->CreateNtupleDColumn(“kinEnergy”);[B1RunAction.cc|attachment]
   analysisManager->CreateNtupleDColumn(“kinEnergy2”);

(upload://uGBGEMCRdDK7eaX13LxInnaZYPV.cc) (7.9 KB)

Then I compiled and runned without errors

logB1-2.txt (97.9 KB)

I got a B1.root file, so I opened it and here the 2 graphs

again the twro graphs are identical, then I’m not sure if really secondary particles have the same kinetic energy or I’m plotting the just particles after one of the two targets…

Moreover I notice that:
1. I simulated 45GeV positrons
2. By the particleID plots I see there are 4 positrons (particle id =-11) and one electron both after first and second target.
**3. By the kinetic energy plots I see there are 4 particles having 45GeV kinetic and one having a very low enegy. **
Given that I simulate 45GeV positrons hitting the targets, it looks like that 4 positrons didn’t interact in the targets and just one interacted…but, how is it possible that I get 5 particles? If just one positron interacted, I should have a couple electron/positron and the 4 primary positrons…not the 4 primary positrons and just one electron!

PARTICLE DIRECTIONS
I used the MomentumDirection functions…so

1. In B1SteppingAction:
   a. I added
   #include “G4ParticleMomentum.hh”

b.I added

if (volume == fScoringVolume) {
double MomDirx = step->GetTrack()->GetMomentumDirection().x();
analysisManager->FillNtupleDColumn(4, MomDirx);
double MomDiry = step->GetTrack()->GetMomentumDirection().y();
analysisManager->FillNtupleDColumn(5, MomDiry);
double MomDirz = step->GetTrack()->GetMomentumDirection().z();
analysisManager->FillNtupleDColumn(6, MomDirz);

	}
	if (volume == fScoringVolume2) {
		 double MomDirx2 = step->GetTrack()->GetMomentumDirection().x();
		 analysisManager->FillNtupleDColumn(7, MomDirx2);
		 double MomDiry2 = step->GetTrack()->GetMomentumDirection().y();
		 analysisManager->FillNtupleDColumn(8, MomDiry2);
		 double MomDirz2 = step->GetTrack()->GetMomentumDirection().z();
		 analysisManager->FillNtupleDColumn(9, MomDirz2);

	}

B1SteppingAction.cc (6.5 KB)

2. In B1RunAction
   a.I added
   analysisManager->CreateNtupleDColumn(“MomentumDirection.x”);
   analysisManager->CreateNtupleDColumn(“MomentumDirection.y”);
   analysisManager->CreateNtupleDColumn(“MomentumDirection.z”);
   analysisManager->CreateNtupleDColumn(“MomentumDirection.x2”);
   analysisManager->CreateNtupleDColumn(“MomentumDirection.y2”);
   analysisManager->CreateNtupleDColumn(“MomentumDirection.z2”);

B1RunAction.cc (8.3 KB)

Then I compiled and runned without errors, I opened the Root file, so I get the three momentum directions after the first target

And after the second one

In this case it looks like that the momentum directions after the two targets are different…

Do you think that I did mistakes or not?

[anna]

The code looks reasonable.
Remember that you store only information about what happens inside either the first or the second target. All the electrons that are created (and deposit energy) in between or outside the targets - you do not see. If you want to check what happens to those particles add /tracking/verbose 1 to your .mac macro and analyse the output. Initial particle will have parentID=0 and trackID=1. Next, all the particles will have an incremented trackID and parentID pointing to the parent.