Unexpected GRIFFIN gamma-ray efficiency when using SCEPTAR

[mwinokan]

When simulating GRIFFIN efficiency with SCEPTAR paddles present unexpected and unphysical efficiencies are observed above about 250 keV. Adding SCEPTAR increases amount of high energy gamma-rays deposited in the GRIFFIN crystals.

Below are the efficiency curves that I get:

In the following energy deposition histograms the last bin has more counts with SCEPTAR than without. (N.B. y-axis scale is different)

James Smallcombe made some suggestions regarding this issue:

• When commenting out the four MakeImprint calls in DetectionSystemSceptar.cc and rebuilding the behaviour does not go away, although the volumes disappear from the scene tree and events with depNumber > 16 are added to the tree.

• It is possible that another file is changing material or physics properties based on the SCEPTAR flag being set.

Help is appreciated, I'm happy to provide more information and histogram examples if needed.

Thanks, Max

[AndrewMac1]

Hey Max,

A couple questions, first what version of Geant4 are you using? Just so I can try to reproduce this. How are you getting your energy histograms? Your binning is quite large. Why is this a projection? For the efficiency you should just be looking at a singles spectra so I am a little confused what you are projecting. For completeness when you add SCEPTAR do you also have the vacuum chamber or no?

Cheers, Andrew

[mwinokan]

Hi Andrew,

Thanks for getting back to me. I'm running Geant4.10.02.

Essentially the two main histograms I've been making are:

tree1->Draw("depEnergy:detNumber>>h0()"); tree1->Draw("depEnergy>>h1(200)","systemID==1000 && detNumber <= 16");

Which produce something like this:

I then take the last bin of the latter histogram to get the top 0.5% of energies. I know this is a very loose definition of the photo-peak but I think that is a different problem all-together, there are still way too many counts at high energies.

And yes I do use the 8pi vacuum chamber with SCEPTAR as well.

Thanks again,

Max

[AndrewMac1]

Hey Max,

You should throw some binning into the first histogram and you don't need to make a histogram of channel vs energy, just do something like: tree1->Draw("depEnergy>>h0(5000,0.,5000.)","systemID == 1000");. If you don't bin the histogram you can end up having Compton in thick bins effecting the overall efficiency. You can also use NTuple2EventTree to sort the simulated data. This gives you an analysis tree just like experimental data and it can be analyzed in the same way. I have a summer student here who is going to run some similar simulations but only at 2 MeV to see if we see the same effects. We will get back to you once they are done.

Cheers, Andrew

[mwinokan]

Thank you very much, I will give this a try.

Max

[mwinokan]

I'm afraid that using t_withSCEP->Draw("depEnergy>>h_withSCEP(1100,0.,1100.)","systemID==1000"); didn't solve my problem. See the histograms below. The photopeak is still more populated with SCEPTAR.

The efficiency, ratio of gammas in photopeak to the 1 million beam gammas, for 1 MeV:

no SCEPTAR: 11.5736%
with SCEPTAR: 16.2839%

I will take a look at using NTuple2EventTree now.

[r3dunlop]

What do the top vs the bottom plots represent?

[mwinokan]

My apologies for the lack of clarity.

Both plots should show the energy spectrum of gamma rays deposited in the Ge crystals for an isotropic distribution of 1 MeV gamma rays originating at the centre of GRIFFIN (which I select for by cutting for systemID==1000 as is assigned to steps inside the crystals in SteppingAction.cc). The bottom plots are a zoomed in section to show the low energy (below photopeak) behaviour.

Hope this helps.

[r3dunlop]

The first thing I notice is that even though the efficiency is higher, there is an overall fewer number of counts in the plot with SCEPTAR. Also, the efficiencies look somewhat like the difference between addback with sceptar in and non-addback with sceptar out. Any chance the SCEPTAR data is using clover addback while the non sceptar version isn't?

[mwinokan]

I had attributed the difference in overall counts to gamma rays absorbed in the scintillators but it seems that these do not account for the over 100k difference in counts, so maybe you are onto something there:

(the above is produced by cutting for systemID==5000)

I will have to look into whether addback is used in these, have not used this functionality yet. The only difference between my Geant4 macros is the line /DetSys/det/addSceptar 20 which I comment out for the data without SCEPTAR.

How would one usually enable the addback functionality?

[AndrewMac1]

Add back functionality is not really accessed until after the raw simulation is sorted. If you are drawing from the raw simulation it should have no idea what addback is, as it is constructed from the event information after running the simulation. Can you confirm this by sorting your simulations and look at a singles vs addback spectra please? Once the simulations of the student finishes tomorrow we will be able to compare our results to this.

[mwinokan]

Thanks for the information, Andrew. I have just been using the raw simulation data so far. I will see about getting GRSISort and NTuple2EventTree set up on my PC, hopefully can have some basic data to compare tomorrow.

[AndrewMac1]

I will double check this today to make sure I am not giving you the wrong information. Sorting the data should help I think.

[AntonNI8]

Hey everyone. I ran a few quick simulations, each with 10^7 events, at 100 KeV and 2000 KeV in order to compare to the first graph. I've provided an overall plot, as well as two others so that you can see the data at both energies more closely. I ran simulations with both the vacuum chamber and SCEPTAR (blue), with only the vacuum (yellow), with only SCEPTAR (orange), and with neither (green). My results at 2000.0 KeV are particularly different from yours, with no appearance of the aforementioned unexpected efficiency.

[AndrewMac1]

Just to add to this, this was sorted with NTuple2EventTree and not drawn from the raw simulation. @mwinokan can you try to sort your simulations and let us know if you everything is working alright?

[mwinokan]

Thank you Anton, Andrew. I'm glad the discrepancies did not appear in your results. I'm working on sorting my data now. Will keep you updated on my progress.

[mwinokan]

@AntonNaimIbrahim Could you please run me through your procedure to sort the data? A settings file for NTuple2EventTree would be particularly helpful.

Thanks,

Max

[AndrewMac1]

There is a Settings.dat file that comes when you download the code in the NTuple2EventTree directory, it should have all the nominal settings I believe. For his analysis, to make it easier, he turned off the energy resolution. You can do in the settings file such that the energy isn't smeared and the peak lies in a single bin. Then just integrate that single bin for the number of counts. Is this what you were looking for? Or more of a set of commands?

[mwinokan]

Thanks Andrew. I had been using the provided Settings.dat file and had come across the error:

terminate called after throwing an instance of 'std::out_of_range'
  what():  vector::_M_range_check
Aborted (core dumped)

after running ./NTuple2EventTree -sf Settings.dat -if withSCEPTAR.root

I wanted to see if this issue was possibly arising from using an incorrect settings file, but if Anton is only removing energy resolution definitions then I would guess that this is not the case.

[r3dunlop]

Are there any other messages before this exception is thrown? That will help someone pin it down.

[mwinokan]

Here's everything I see:

[mwinokan@udon NTuple2EventTree]$ ./NTuple2EventTree -sf Settings.dat -if withSCEPTAR.root 
will read from 1 files
terminate called after throwing an instance of 'std::out_of_range'
  what():  vector::_M_range_check
Aborted (core dumped)
[mwinokan@udon NTuple2EventTree]$ 

[r3dunlop]

Can you put cout lines in Converter and recompile? If we know around what line it’s happening on, we can figure out which vector is throwing.

On May 17, 2018, at 5:47 PM, mwinokan notifications@github.com<mailto:notifications@github.com> wrote:

Here's everything I see:

[mwinokan@udon NTuple2EventTree]$ ./NTuple2EventTree -sf Settings.dat -if withSCEPTAR.root will read from 1 files terminate called after throwing an instance of 'std::out_of_range' what(): vector::_M_range_check Aborted (core dumped) [mwinokan@udon NTuple2EventTree]$

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/GRIFFINCollaboration/detectorSimulations_v10/issues/55#issuecomment-390023790, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AFLTwAFhIuAJ83x-7bDHwAKqh7jFkglcks5tzfAEgaJpZM4UBy35.

[r3dunlop]

Actually you might want to start with NTuple2EventTree in general. Either way, finding the line this issue is happening on will make fixing it much easier. Don’t forget to but an endl after the out statements.

[mwinokan]

Thanks Ryan. Editing the files now.

[mwinokan]

It appears that the call to fRandom.Gaus() in Converter::Run() seems to be causing the error to be thrown. The line in question is:

smearedEnergy = fRandom.Gaus(fDepEnergy,fSettings->Resolution(fSystemID,fDetNumber,fCryNumber,fDepEnergy));

Around line 215 of Converter.cc

[r3dunlop]

Any chance you have some system ID, detector number or crystal number that doesn’t exist?

That line calls fSettings->Resolution which does (in Settings.hh)

 double Resolution(int systemID, int detectorID, int crystalID, double en) {
        if(fResolution.find(systemID) != fResolution.end()) {
            return fResolution[systemID].at(detectorID).at(crystalID).Eval(en);
        }
        return 0.;
    }

My guess is for the system ID being passed there is no correct detectorID or crystal ID, and the at call is throwing. How long does it run for before throwing the exception? Can you change the code to

 double Resolution(int systemID, int detectorID, int crystalID, double en) {
        if(fResolution.find(systemID) != fResolution.end()) {
            try{
                  return fResolution[systemID].at(detectorID).at(crystalID).Eval(en);
            }
            catch (const std::out_of_range& e) {
                 std::cout << "Out of Range error throwing " << systemID << " " << detectorID << " " << crystalID << std::endl;
            }
        }
        return 0.;
    }

You may have to include <'stdexcept'> (It won't let me type this without breaking up the brackets...?)

It looks to me that there is no crystal vector for each SCEPTAR detector and that's what's causing it to fail. I don't really know this code, so this is a guess. Try the above, and we'll take it from there.

[mwinokan]

Hi Ryan, thanks for your help. Below are some of the things I tried this morning:

I did modify SteppingAction.cc to change system-, detector-, and crystalID's for generic targets that I was using to do some tests for new detector designs, have commented these lines out and reproduced a new root file, which still has the same problem.

Running with your error trap I see now that the out of range ID's are:

systemID = 1000
detectorID = 8
crystalID = 7

Seems to me that the crystalID of 7 is the problem. Looking at which crystalID's are present in the ntuple it seems that somewhere 4 is added to the crystalID (5,6,7,8 are present instead of 0,1,2,3).

Running the macro with a clean build of detectorSimulations, yields crystalID's of 5,6,7,8 and -1.

For completeness, a root tree generated from a SCEPTAR-less macro generates no errors.

I tried adding some lines to set fCry = 0 before SCEPTAR hit-trackers are added in SteppingAction.cc. This didn't fix the problem.

I'm going to continue trying to identify where in detectorSimulations fCry is being incorrectly set.

Some other info:

Time statistics of running time ./NTuple2EventTree -sf Settings.dat -if withSCEPTAR_x.root:

real    0m1.831s
user    0m0.425s
sys 0m0.317s

Full output:

[mwinokan@udon NTuple2EventTree]$ ./NTuple2EventTree -sf Settings.dat -if withSCEPTAR_x.root 
### Converter::Converter()
will read from 1 files
### Converter::Converter() FIN
### Converter::Run()
### Converter::Run() CHECK1
### Converter::Run() CHECK2
### Converter::FillDetectors()
### Converter::FillDetectors() FIN
### Converter::Run() CHECK3
### Converter::Run() CHECK3.1
### Converter::Run() CHECK3.3
Out of Range error throwing 1000 8 7
### Converter::Run() CHECK3.4
### Converter::Run() CHECK4
### Converter::AboveThreshold()
terminate called after throwing an instance of 'std::out_of_range'
  what():  vector::_M_range_check
Aborted (core dumped)
[mwinokan@udon NTuple2EventTree]$ 

[r3dunlop]

So all GRIFFIN crystals are either 5,6,7,8 or -1 ? There are no proper crystals? Interesting.

[mwinokan]

Yes that's correct.

[r3dunlop]

Interesting, so adding SCEPTAR changes the GRIFFIN crystal numbers in the Raw Simulation data. @VinzenzBildstein do you have any thoughts why this might be happening? It might take me a while to dig through and find that error.

[mwinokan]

Curiously, the process types seem to be incorrect as well: no scintillation is present and lots of process-types of 0 and 1. (x is crystalID, y is processType)

[AndrewMac1]

Hmmm...just as a side note to this, I just looked at a simulation I ran around 1-2 months ago with SCEPTAR and when I made a hit pattern it looks ok. All the GRIFFIN channels seem to be properly assigned. This data wasn't drawn from directly from the tree, but it would be weird if things were improperly assigned running the simulations, then properly assigned in sorting.

[AndrewMac1]

Drawing cryNumber from the tree I get a histogram going from 0 to 4 and detNumber goes from 0 to 16. @mwinokan are you axes backwards from what you said in the previous plot? (y-axis is cryNumber x-axis is processType)

[r3dunlop]

I think the throw error shows there are crystal numbers that are wrong, no?

[mwinokan]

Axes are definitely correct:

[AndrewMac1]

Yeah it looks like that, just wondering why his is different than mine, more a note for Vinzenz for a timeline if something has changed.

[AndrewMac1]

Sorry I didn't mean correct in terms of variables I just meant they were backwards.

[AndrewMac1]

Either way something is wrong here.

[mwinokan]

If somebody could run the following macro on their Griffinv10 build and take a look at the detNumber and processType branches that would be great. Just to do a sanity check and see if maybe there is just something wrong with how I am running things.

/run/initialize
/control/verbose 2
/run/verbose 2
/event/verbose 0
/tracking/verbose 0

/DetSys/world/material Vacuum

/DetSys/det/addSceptar 20
/DetSys/app/add8piVacuumChamber
# /DetSys/det/addPaces 1

/DetSys/det/SetCustomShieldsPresent 0
/DetSys/det/SetCustomRadialDistance 11 cm
/DetSys/det/SetCustomExtensionSuppressorLocation 0
/DetSys/det/includeGriffinHevimet 0

/DetSys/det/SetCustomDeadLayer 1 1 0
/DetSys/det/addGriffinCustomDetector 1
/DetSys/det/SetCustomDeadLayer 2 2 0
/DetSys/det/addGriffinCustomDetector 2
/DetSys/det/SetCustomDeadLayer 3 3 0
/DetSys/det/addGriffinCustomDetector 3
/DetSys/det/SetCustomDeadLayer 4 4 0
/DetSys/det/addGriffinCustomDetector 4
/DetSys/det/SetCustomDeadLayer 5 5 0
/DetSys/det/addGriffinCustomDetector 5
/DetSys/det/SetCustomDeadLayer 6 6 0
/DetSys/det/addGriffinCustomDetector 6
/DetSys/det/SetCustomDeadLayer 7 7 0
/DetSys/det/addGriffinCustomDetector 7
/DetSys/det/SetCustomDeadLayer 8 8 0
/DetSys/det/addGriffinCustomDetector 8
/DetSys/det/SetCustomDeadLayer 9 9 0
/DetSys/det/addGriffinCustomDetector 9
/DetSys/det/SetCustomDeadLayer 10 10 0
/DetSys/det/addGriffinCustomDetector 10
/DetSys/det/SetCustomDeadLayer 11 11 0
/DetSys/det/addGriffinCustomDetector 11
/DetSys/det/SetCustomDeadLayer 12 12 0
/DetSys/det/addGriffinCustomDetector 12
/DetSys/det/SetCustomDeadLayer 13 13 0
/DetSys/det/addGriffinCustomDetector 13
/DetSys/det/SetCustomDeadLayer 14 14 0
/DetSys/det/addGriffinCustomDetector 14
/DetSys/det/SetCustomDeadLayer 15 15 0
/DetSys/det/addGriffinCustomDetector 15
/DetSys/det/SetCustomDeadLayer 16 16 0
/DetSys/det/addGriffinCustomDetector 16
# /DetSys/det/update

# comment out the following /vis/ commands for a batch process
# /vis/open OGL
# /vis/viewer/set/viewpointThetaPhi 70 20 deg
# /vis/viewer/zoomTo 0.3
# /vis/viewer/zoomTo 0.5
# /vis/viewer/set/style surface
# /vis/drawVolume
# /vis/scene/endOfEventAction accumulate 100
# /vis/scene/add/trajectories smooth

/DetSys/gun/particle gamma
/DetSys/gun/efficiencyEnergy 1000.00 keV
/run/beamOn 1000000

[AndrewMac1]

Hey @mwinokan I got the same as you did, can you try to run with the following and see if you get the same thing? This one worked for me where your run macro did not.

###################### PHYSICS LIST OPTIONS ############################################# /DetSys/phys/SelectPhysics emlivermore /run/initialize /process/em/fluo true /process/em/auger true /process/em/pixe true

/DetSys/phys/SelectPhysics QGSP_BERT_HP

/run/initialize

###################### GRIFFIN DETECTOR PROPERTIES ######################################

SetCustomShieldsPresent 1 (include suppressors)

SetCustomShieldsPresent 0 (do NOT include suppressors)

SetCustomShieldsPresent -1 (only include side and back suppressors, ie. no extension)

SetCustomRadialDistance 11 cm (leave this at 11 cm, even in back mode)

SetCustomExtensionSuppressorLocation 0 (forward mode)

SetCustomExtensionSuppressorLocation 1 (back mode)

includeGriffinHevimet 0 (no)

includeGriffinHevimet 1 (yes)

######################################################################################### /DetSys/det/SetCustomShieldsPresent 0 /DetSys/det/SetCustomRadialDistance 11 cm /DetSys/det/SetCustomExtensionSuppressorLocation 0 /DetSys/det/includeGriffinHevimet 0

/DetSys/det/SetCustomDeadLayer 1 1 0 /DetSys/det/addGriffinCustomDetector 1 /DetSys/det/SetCustomDeadLayer 2 2 0 /DetSys/det/addGriffinCustomDetector 2 /DetSys/det/SetCustomDeadLayer 3 3 0 /DetSys/det/addGriffinCustomDetector 3 /DetSys/det/SetCustomDeadLayer 4 4 0 /DetSys/det/addGriffinCustomDetector 4 /DetSys/det/SetCustomDeadLayer 5 5 0 /DetSys/det/addGriffinCustomDetector 5 /DetSys/det/SetCustomDeadLayer 6 6 0 /DetSys/det/addGriffinCustomDetector 6 /DetSys/det/SetCustomDeadLayer 7 7 0 /DetSys/det/addGriffinCustomDetector 7 /DetSys/det/SetCustomDeadLayer 8 8 0 /DetSys/det/addGriffinCustomDetector 8 /DetSys/det/SetCustomDeadLayer 9 9 0 /DetSys/det/addGriffinCustomDetector 9 /DetSys/det/SetCustomDeadLayer 10 10 0 /DetSys/det/addGriffinCustomDetector 10 /DetSys/det/SetCustomDeadLayer 11 11 0 /DetSys/det/addGriffinCustomDetector 11 /DetSys/det/SetCustomDeadLayer 12 12 0 /DetSys/det/addGriffinCustomDetector 12 /DetSys/det/SetCustomDeadLayer 13 13 0 /DetSys/det/addGriffinCustomDetector 13 /DetSys/det/SetCustomDeadLayer 14 14 0 /DetSys/det/addGriffinCustomDetector 14 /DetSys/det/SetCustomDeadLayer 15 15 0 /DetSys/det/addGriffinCustomDetector 15 /DetSys/det/SetCustomDeadLayer 16 16 0 /DetSys/det/addGriffinCustomDetector 16

/DetSys/det/addDescant 70

###################### LaBr3 DETECTOR PROPERTIES ########################################

addLanthanumBromide #_of_dets radial_pos_in_cm null

addAncillaryBGO #_of_dets radial_pos_in_cm include_hevimet

/DetSys/det/addLanthanumBromide 8 16.5 0

/DetSys/det/addAncillaryBGO 8 16.5 0

###################### SCEPTAR and Others ############################################### /DetSys/det/addSceptar 20 /DetSys/app/add8piVacuumChamber

/DetSys/app/add8piVacuumChamberAuxMatShell 20 mm

###################### Griffin Structure ###############################################

/DetSys/app/addGriffinStructure 0 (include both up- and down-stream halves)

/DetSys/app/addGriffinStructure 1 (include both upstream half)

/DetSys/app/addGriffinStructure 2 (include both downstream half)

/DetSys/app/addGriffinStructure 1

###################### VERBOSE LEVELS ################################################### /control/verbose 1 /run/verbose 0 /event/verbose 0 /tracking/verbose 0

/DetSys/gun/particle gamma /DetSys/gun/efficiencyEnergy 1000 keV

###################### BEAM ON ##########################################################

/run/beamOn 100000

[skwzrd]

Hello,

I have been working on a project using detectorSimulations_v10 myself, and I noticed specifying "Vacuum" for a world material does not work. It crashes my executable. I have to use:

/DetSys/world/material G4_Galactic

Perhaps you are not actually specifying a vacuum in your simulations?

Michael

[mwinokan]

Thanks @AndrewMac1, just ran this and upon initial inspection in root the data looks good. (normal crystal numbers). I will now see about the efficiency of the raw data and see if there is still a problem. Then will proceed to sorting.

I suppose including the extra physics processes at the start of your macro might be making the difference, if not the ordering of volume creation.

@Michael-Hladun, thank you for the heads up, I will give this a try as well.

[skwzrd]

@mwinokan I also declared fMatWorldName = "G4_Galactic"; in DetectorConstruction.cc. Previously it was G4_AIR.

Edit: I tried specifying /DetSys/world/material Vacuum in my .mac files and it seems to work now. My executable does not crash anymore.

[jsmallcombe]

That is precisely what the /DetSys/world/material does. /DetSys/world/material Vacuum definitely works for me, I only get charged particle detection after calling it (as the default is indeed Air).

[mwinokan]

So I'm pleased to say that from some basic manipulation of the raw data at 1 MeV it looks like the problem has been fixed. Here are my efficiencies:

w/ 8pi chamber: 9.858
w/ SCEPTAR + 8pi chamber: 9.748 %
w/ SCEPTAR + 8pi chamber (em processes and emlivermore commented out): 11.484 % 

I am going to run the simulation with higher beam numbers and for more energy points over the weekend, and will keep the issue open until I've checked those results.

Thank you everyone for your help on this!

[VinzenzBildstein]

Apparently this issue happened when setting up SCEPTAR before GRIFFIN. This is due to the fact that the detector and crystal numbers for GRIFFIN were calculated from the assembly volume index (that was extracted from the volume name), which changes when other detectors are implemented first. There is no reason to do that, the name also contains the detector and crystal number, so we can just use that. This will be fixed in response to issue #56.