Dear Philip,
Thank you for your immediate response!
>With the total cross section given by the spline, along with your flux, you can calculate the expected total number of events as >flux*(cross section), and use it to scale your event histogram.
I see, I think now I start to understand how to treat my problem.
EventRecord holds cross section that the "selected" event. It also hold an interaction probability that the specific event make happen.
So, when I pick up one event and ask "what is the probability that this event happens?", it's simply obtained with
prob = evt->Probability().
The sum of prob of all event will be how many neutrino actually interacted, when 2000 primary neutrino is injected
(since I generated 2000 events). In other words, this is the event weight.
Am I correct?
(It looks very simple, and I'm a bit ashamed that I didn't understand it earlier)
>This gets a little tricky when your detector contains multiple elements, because you have to get multiple splines, and weight >them correctly, so I wonder if any of the GENIE developers know an easier way (or some existing tools that do this)
Here I think you are trying to explain another way to weight event, which I'm using for our group's neutrino MC too.
weight = 1 - exp (-N * total cross section * path length) (= N * total cross section * path length @ low E)
In this case ALL injected neutrino will make interaction.
Yes that would be nice if we have any good tools to do it easier! (I already wrote the code in my test program, though).
I hope now I understand how to weight genie event correctly, but weighting is always tricky and dangerous if users don't know how the events are generated.
So it would be nice if experts can add one short section "how to weight genie event" in genie user's manual in near future :).
Thank you again,
Kotoyo
On Feb 20, 2013, at 9:49 AM, Philip Rodrigues <[log in to unmask]> wrote:
Dear Kotoyo,
It sounds like you want to make cross section plots (or something that you can easily calculate given the cross section, like events in a given exposure in your detector). For that you'll need the cross section spline file. I found it easier to work with a ROOT version of the spline file, which you can create from the XML spline file using gspl2root. With the total cross section given by the spline, along with your flux, you can calculate the expected total number of events as flux*(cross section), and use it to scale your event histogram.
This gets a little tricky when your detector contains multiple elements, because you have to get multiple splines, and weight them correctly, so I wonder if any of the GENIE developers know an easier way (or some existing tools that do this)
Hope that helps,
Phil
On 02/19/2013 07:17 PM, Kotoyo Hoshina wrote:
Hello,
I'm a new beginner of genie neutrino generator, and have a question about how to weight genie events.
The "weight" I mean is NOT same thing as the field I can get from EventRecord with following function
evt->Weight();
which for now I set always 1.
I generated 2000 genie event with gevgen command, and converted result .root file to gst format root file.
When I draw 1D histogram of El parameter (energy of final-state lepton), the vertical axis will be number of generated events per energy bin. This is just unweighted MC histogram and can't be used for physics analysis.
Now I want to convert this number to physics plot, number of interactions actually happened per 2000 primary neutrinos. Due to the small cross section, the vertical axis should become very small (~1e-36).
According to Genie user manual, genie is using internal global probability scale which is chosen from highest input energy and longest path length. At first, I thought this number will be the constant event weight to get actual number of events :
event_weight = gmcj_driver->GlobProbScale()
...
gst->Draw("El", event_weight);
but after I read GMCJDriver carefully I noticed it won't work, because the event loop inside GMCJDriver actually loops much more than 2000. There is no way to know how many neutrinos are actually injected, thus no way to get normalization factor.
It is always simple to calculate event weight of a specific event : it's just the interaction probability (and easy to obtain from EventRecord), but weighting entire MC sets strongly depends on how these MC sets are generated...
and I couldn't find any document about weighting procedures, instead, I found re-weighting staff (which currently I don't need.)
Is there any document or example code?
thank you,
kotoyo
|