Performing cuts on restG4 simulated data

Newbie
1

Hi everyone,

I have a naïve question. I looked for an answer but couldn’t find any complete one.

I am generating particles outside of my detector, neutrons for example.
I would like to access the primary interaction point of these neutrons, and select those who first interacted in a given physical volume of my geometry.

I looked in the documentation for something around the TRestRun metadata object, but I’m not sure what I should look for.
If you have any example, advice, or even place where I could find some answers, I would be really glad.

Thanks,
Cloe

2

Hi Cloe,

I believe what you may need is related with the following methods implemented inside TRestGeant4Event.

    TVector3 GetMeanPositionInVolume(Int_t volID) const;
    TVector3 GetFirstPositionInVolume(Int_t volID) const; 
    TVector3 GetLastPositionInVolume(Int_t volID) const;
    TVector3 GetPositionDeviationInVolume(Int_t volID) const;

You might access your file and retrieve the event pointer.

restRoot YourRestG4File.root
run0->GetEntry(N);
TVector3 pos = ev0->GetFirstPositionInVolume(0);

It will return a TVector3 position where you can access the individual components as pos.X(), pos.Y() and pos.Z().

The relation between the volume Id and the name on the geometry is stored inside TRestGeant4Metadata.

I added a printout at TRestGeant4Metadata::PrintMetadata

====

The methods at TRestGeant4Event will allow you to test or to implement a macro to get some results. However, it would be interesting to implement an observable inside perhaps a dedicated analysis process, or inside TRestGeant4AnalysisProcess itself. So that this observable gets into the main event data processing and can be added to the analysis tree on demand.

There are already some observables such as volumeMeanPos inside TRestGeant4AnalysisProcess.

One could think of a new observable volumeFirstPosition.

https://sultan.unizar.es/rest/classTRestGeant4AnalysisProcess.html

1 Like
3

Hi @jgalan,

Thank you for your quick answer.

The methods at TRestGeant4Event will allow you to test or to implement a macro to get some results.

I think I will use them in a first time.
But as you say, in the end, it would probably be a better idea to implement some observables in the event processing pipeline, to have them in the RestG4 Root output file.

I will try and probably come back here to ask some more questions.

Thank you again,
Cloe

1 Like
4

It seems @lobis is working already on something related.

I think it would be better to differentiate different analysis processes.

We could have

  • TRestGeant4AnalysisProcess for basic observables.
  • TRestGeant4VolumeAnalysisProcess for observables related with geometry volumes.
  • TRestGeant4ParticleAnalysisProcess for observables related with particle names.

Perhaps there is a better structure/naming. But it is good that processes do not get huge. Also for computational reasons.

5

For the moment I’ve done something like that, using TRestGeant4Event:

TRestRun run("file.root");
  
TRestGeant4Event* event = run.GetInputEvent<TRestGeant4Event>();

for (size_t i = 0; i < run.GetEntries(); i++) {
  run.GetEntry(i);
  cout << i << endl;
  cout << "First position in sensitive volume (" << event->GetFirstPositionInVolume(0).x() << ", " << event->GetFirstPositionInVolume(0).y() << ", " << event->GetFirstPositionInVolume(0).z() << ") mm" << endl;
  cout << "First position in active volume (" << event->GetFirstPositionInVolume(1).x() << ", " << event->GetFirstPositionInVolume(1).y() << ", " << event->GetFirstPositionInVolume(1).z() << ") mm" << endl;
    
 }

But as output I only have nan for the active volume (ID 1), while I have coherent positions for the sensitive volume (ID 0).

If I do

TRestGeant4Metadata* geant4Metadata = (TRestGeant4Metadata*)run.GetMetadataClass("TRestGeant4Metadata");
  
for (size_t i = 0; i < geant4Metadata->GetNumberOfActiveVolumes(); i++) {
  cout << geant4Metadata->GetActiveVolumeID(geant4Metadata->GetActiveVolumes()[i]) << endl;
}

I have the output IDs 0 and 1 as expected.

I checked with the visu REST_Geant4_ViewEvent and I have events depositing energy in my volume 1.

6

What is the output of event->PrintEvent()?

7

You will find attached an example of a PrintEvent for an event depositing energy in my active volume with ID 1 (name “ID_lsAMT”)
test.txt (248.5 KB)

8

I think there are some issues identifying the volumes.

I am trying with the example 04.MuonScan where we actually define 2 active volumes in the geometry.

TRestGeant4Metadata produces the following output:

              ||                                ++++++++++ Detector +++++++++++                                 ||              
              ||                                 Energy range (keV): (0, 1e+20)                                 ||              
              ||                                 Number of sensitive volumes: 1                                 ||              
              ||                                  Sensitive volume: det_dw_01                                   ||              
              ||                                  Number of active volumes: 2                                   ||              
              ||                        Name: det_dw_01, ID: 0, maxStep: 1mm , chance: 1                        ||              
              ||                        Name: det_up_01, ID: 1, maxStep: 1mm , chance: 1                        ||              
              ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

I print the track #7 here. The volume name is det_dw_01, thus the id should be 0.

 * TrackID: 7 - Particle: e- - ParentID: 1 - Parent particle: mu- - Created by 'muIoni' with initial KE of 17.12 keV
   Initial position (-224.115, -30.444, -321.290) mm at time 2.91 ns - Time length of 0.16 ns and spatial length of 9.98 mm
      - Hit 0 - Energy: 0 - Process: Init - Volume: det_dw_01 - Position: (-224.115, -30.444, -321.290) mm - Time: 2.91 ns - KE: 17.12 keV
      - Hit 1 - Energy: 228.69 eV - Process: msc - Volume: det_dw_01 - Position: (-224.063, -31.162, -321.946) mm - Time: 2.93 ns - KE: 16.89 keV
      - Hit 2 - Energy: 492.31 eV - Process: msc - Volume: det_dw_01 - Position: (-223.678, -31.773, -322.564) mm - Time: 2.94 ns - KE: 16.40 keV
      - Hit 3 - Energy: 1.43 keV - Process: msc - Volume: det_dw_01 - Position: (-223.406, -32.124, -323.359) mm - Time: 2.95 ns - KE: 14.98 keV
      - Hit 4 - Energy: 738.05 eV - Process: msc - Volume: det_dw_01 - Position: (-223.055, -32.513, -323.911) mm - Time: 2.96 ns - KE: 14.24 keV
      - Hit 5 - Energy: 792.08 eV - Process: msc - Volume: det_dw_01 - Position: (-222.858, -32.947, -324.430) mm - Time: 2.97 ns - KE: 13.45 keV
      - Hit 6 - Energy: 91.86 eV - Process: msc - Volume: det_dw_01 - Position: (-222.601, -33.057, -324.606) mm - Time: 2.98 ns - KE: 13.35 keV
      - Hit 7 - Energy: 581.28 eV - Process: msc - Volume: det_dw_01 - Position: (-222.426, -32.868, -324.899) mm - Time: 2.98 ns - KE: 12.77 keV

However, if I add a bit of debugging output in the method TRestGeant4Track::GetEnergyInVolume() that is used internally in the method that we are testing, I get the following.

TRestGeant4Track::GetEnergyInVolume(). Track : 7
Hit: 0 volId: 1
Hit: 1 volId: 1
Hit: 2 volId: 1
Hit: 3 volId: 1
Hit: 4 volId: 1
Hit: 5 volId: 1
Hit: 6 volId: 1
Hit: 7 volId: 1
Hit: 8 volId: 1

Thus, there is no hit with id 0, as it should be. I hope @lobis could have a quick solution for this.

I opened an issue Problem retrieving volume id at `GetEnergyInVolume` · Issue #99 · rest-for-physics/geant4lib · GitHub

1 Like
9

I will take a look tomorrow

1 Like
10

I may have found a hint on what is not working.

I looks like the IDs of your physical volumes does depend on the order you first defined them in the setup.gdml file.

I’m in a hurry and I don’t have time to elaborate what I found, but I will edit later my comment. Just in case it can help.

Cloe

11

Right! It seems there is some confusion between active volume Ids and geometry volumes.

Some discussion run out in parallel at