Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add -T option to gevgen_atmo to specify livetime instead of number of events #92

Merged
merged 1 commit into from Oct 9, 2020

Conversation

tlatorre-uchicago
Copy link
Contributor

@tlatorre-uchicago tlatorre-uchicago commented Feb 19, 2020

This commit adds the -T command line option to gevgen_atmo to specify a
livetime in seconds instead of a total number of events using the -n option.

To implement this gevgen_atmo first integrates the 3D flux histogram and
multiplies it by the livetime to get the total number of neutrinos that should
be simulated per square meter. This number is then multiplied by the area of
the flux surface to determine the total number of neutrinos that should be
simulated.

Instead of looping over a fixed number of events, gevgen_atmo then continues
throwing neutrinos until the number of neutrinos thrown divided by the global
probability scale is greater than the total number of neutrinos that should be
simulated.

Some potential issues that should be considered before merging:

  1. The method I'm using definitely is dimensionally correct and I verified that it behaves as expected using a very simple Monte Carlo I wrote in python, however the actual flux driver and GMCJDriver are a bit more complicated and so I'm not 100% sure that there is something I'm missing and whether it is OK to simply scale the number of flux neutrinos thrown by the global probability scale to get the effective total number of neutrinos simulated.

  2. Technically the code in the flux driver allows one to do a scheme whereby the neutrinos are thrown not according to a global probability rescaling but by rescaling each bin and then assigning weights. The whole method I use in this pull request will not work if this method is activated. Currently in gevgen_atmo it's not activated and I added a comment explaining that if it is ever activated then the exposure time method I use will not work.

  3. This pull request has not gone through extensive testing. I did some very basic testing with a similar commit applied to version 2.12.8 which is what I'm currently using, and I just build tested it against master. But it probably needs much more testing before being applied.

@candreop candreop requested a review from nusense Feb 26, 2020
@candreop
Copy link
Member

candreop commented Feb 26, 2020

Without thinking too deeply about it, it seems correct. It is broadly equivalent with what we do for accelerator neutrinos. There, each flux ray thrown incremements a POT counter. The rate of increase of the POT counter per flux ray is determined by the flux simulation. Here, each flux ray thrown increments a time counter. The rate of increase is determined by the input flux histograms. The input flux is given in terms of number of flux neutrinos per unit energy per unit solid angle and per unit time. If you integrate over \theta, \phi and E, that histogram tells you, for each neutrino species, how many neutrinos we have per unit time. As for acceleator neutrinos, the probability scale affects that exposure unit (at each turn we don't simulate single POTs or single seconds, but rather larger chunks that are more likely to yield an interaction).

@candreop
Copy link
Member

candreop commented Feb 26, 2020

I agree this will not work with the option to produce weighted events, but it seems to be that it can be made to work with that option too. Another thing that needs attention is that gevgen_app allows a user to specify energy limits (in which case, I guess, the input flux simulation will need to be integrated within that specified range in order to yield the appropriate time increment per flux ray)

@tlatorre-uchicago
Copy link
Contributor Author

tlatorre-uchicago commented Mar 4, 2020

Thanks for the review! I will update this to take into account the energy bounds, but I'm still not sure how it can be made to work with weighted events. Also, to be clear, currently the gevgen_atmo script does not allow one to use weighted events, although this might have been something planned for the future. Or am I misunderstanding what you were saying?

@tlatorre-uchicago
Copy link
Contributor Author

tlatorre-uchicago commented Apr 1, 2020

I've updated this to hopefully take into account the energy range. I also added two very simple tests to make sure that the integration routines match what we expect. You can run them by doing:

$ cd src/contrib/test
$ make gtestGAtmoFlux
$ cd ../../
$ ./bin/gtestGAtmoFlux

which will output:

$ ./bin/gtestGAtmoFlux
...
[ok] testGetTotalFlux
[ok] testGetTotalFluxInEnergyRange

I'd like to add a more comprehensive test to test the whole scheme, but haven't come up with a clean way to do that yet.

@tlatorre-uchicago tlatorre-uchicago force-pushed the exposure-squashed branch 2 times, most recently from 0fd7c45 to 5bb353e Compare Apr 14, 2020
… events

This commit adds the -T command line option to gevgen_atmo to specify a
livetime in seconds instead of a total number of events using the -n option.

To implement this gevgen_atmo first integrates the 3D flux histogram and
multiplies it by the livetime to get the total number of neutrinos that should
be simulated per square meter. This number is then multiplied by the area of
the flux surface to determine the total number of neutrinos that should be
simulated.

Instead of looping over a fixed number of events, gevgen_atmo then continues
throwing neutrinos until the number of neutrinos thrown divided by the global
probability scale is greater than the total number of neutrinos that should be
simulated.
// event loop
for(int iev = 0; iev < gOptNev; iev++) {
for(int iev = 0; gOptNev > 0 ? iev < gOptNev : 1; iev++) {
Copy link
Member

@nusense nusense Apr 21, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I suggest some ()'s here to help with readability:

 for (int iev = 0; ( (gOptNev > 0) ? (iev < gOptNev) : true) ; iev++)

Copy link
Contributor Author

@tlatorre-uchicago tlatorre-uchicago Apr 22, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What do you think about just turning it into a while loop to make things easier to understand? I originally thought about doing that but didn't want to change too much, but now that you mention it I think it would be much easier to understand that way.

atmo_flux_driver = GetFlux();

// Cast to GFluxI, the generic flux driver interface
GFluxI *flux_driver = dynamic_cast<GFluxI *>(atmo_flux_driver);
Copy link
Member

@nusense nusense Apr 21, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

pedantically I'm not sure this is ever necessary to down(?) cast to a base class unless the derived class hides some base class method. You should be able to just call those base class methods on the derived class or pass derived object pointer to methods that expect the base class. But that was an issue with the original code, so not a required change here.

Copy link
Contributor Author

@tlatorre-uchicago tlatorre-uchicago Apr 22, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good point, I will update this and push a new set of changes.

/* Note: For the method of calculating the total number of events using a
* livetime we *need* to force a single probability scale. So if you change
* the next line you need to make sure that the user didn't specify the -T
* option. */
mcj_driver->ForceSingleProbScale();
Copy link
Member

@nusense nusense Apr 21, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I presume it's thus safe to call ForceSingleProbScale() in all cases, but unsafe to not call it if -T was specified. Am I reading that right?

Copy link
Contributor Author

@tlatorre-uchicago tlatorre-uchicago Apr 22, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes that's right. Do you have any suggestions on how I can update that comment or if it belongs somewhere else? I added that just as kind of a future warning to anyone editing it that they shouldn't change that without adding a guard to check against the -T option.

@candreop candreop added this to the 3.2.0 release milestone Jul 6, 2020
@mroda88 mroda88 merged commit f4d15d7 into GENIE-MC:master Oct 9, 2020
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

4 participants