Bug #18025
Need of update translation of GENIE interaction type
Description
While looking at the interaction type of MC neutrino event (generated with icaruscode fcl/gen/genie/simulation_genie_icarus_simple_workshop.fcl), around the 20% of events have an InteractionType()==1000, that is defined as an offset and do not correspond to any physical interaction.
Same error happens also in sbndcode.
History
#1 Updated by Gianluca Petrillo over 1 year ago
The module is GENIEGen, the point where this happens is nutools:source:nutools/EventGeneratorBase/GENIE/GENIEHelper.cxx line 1593:
int itype = simb::kNuanceOffset + genie::utils::ghep::NuanceReactionCode(record);There may be two cases:
- GENIE helper is simulating some interactions for which it does not map the code to Nuance; in this case, GENIE shoudl update their mapping (
genie::utils::ghep::NuanceReactionCode()) - GENIE is simulating interactions that have no Nuance code; in that case, nutools helper should create specific codes for these GENIE processes
LArSoft code does not get access to the code from GENIE, so we can't inform you on which interaction codes are missing. The FHiCL configuration provided by the reporter should help reproducing the issue, though.
#2 Updated by Robert Hatcher over 1 year ago
Some of GENIE's underlying information is also recorded in a unsullied version in a nusimdata/SimulationBase/GTruth that is associated with every MCTruth object. It would be instructive to see if there is a systematic combination of GTruth.fGint & GTruth.fGscatter that is repeated when this Nuance mapping fails. This would help identify what's going on.
Also, if one could enable GENIEHelper's fcl parameter GHepPrintLevel=13, ala:
https://cdcvs.fnal.gov/redmine/projects/nutools/wiki/GENIEHelper#Printing-GENIE-Event-Record
then we can go back an correlate full GENIE records with MCTruth entries with this nonsense value for simb::MCNeutrino::InteractionType()
For the record, the genie::utils::ghep::NuanceReactionCode(record) code lives here:
https://genie.hepforge.org/trac/browser/generator/trunk/src/GHEP/GHepUtils.cxx#L271
Also a cursory look seems to me that it doesn't properly classify MEC events (but those shouldn't be 20% of a normal sample, I think), so this might be be a "GENIE" issue.
#3 Updated by Andrea Falcone over 1 year ago
I looked for the variable simb::GTruth::fGscatter.
I have found that:
if fGscatter = 10, then InteractionType()=1000,
if InteractionType()=1000 , then fGscatter = 10 or fGscatter = 4.
The event I'm using are in
/pnfs/icarus/scratch/users/cfarnese/workshop/v06_54_00/gen/simulation_genie_icarus_simple_workshop/out/
#4 Updated by Robert Hatcher over 1 year ago
Hmm, that's interesting
Here's GENIE's enumeration of fGscatter
https://genie.hepforge.org/trac/browser/generator/trunk/src/Interaction/ScatteringType.h
4=kScResonant
10=kScMEC,
So, as I suspected GENIE's MEC has no translation to a NUANCE code. This would need a "slot" in MCNeutrino's enumeration to put the value in, and a fix to GENIE's code (and/or a hack to NuTools before GENIE gets a release cut with this fix).
I'm getting a copy of NUANCE and it's scatter type enumeration (soon hopefully) and will look to see where I can start "extending" it to slot in MEC.
The case of fGscatter=4 is a bit more surprising. GENIE's classification of "resonant" scatters are broken down more finely in the NUANCE scheme. GENIE does this translation by looking at the whole event and seeing some topology (#'s of pions, etc). Looks like there might be some case there that isn't handled completely. But to identify the case we would need more full dumps of the event records themselves to figure it out.
This could be done outside of "art"/"larsoft" with basic GENIE functionality. But I'm seriously swamped so if there is someone here at Fermilab who could temporarily be my minion (doing work at my direction) it would happen infinitely faster. They need to be able to setup larsoft (which get's them GENIE), run some gevgen_fnal files (again not hard), and write a bit of C++ as a macro (again not hard).
Volunteers?
#5 Updated by Robert Hatcher over 1 year ago
Also, could I get clarification from the originator which version of GENIE was used, and which set of splines (which genie_xsec UPS product was setup) & the setting of the EventGeneratorList fcl parameter?
For reference:
https://cdcvs.fnal.gov/redmine/projects/nutools/wiki/GENIE_Configuration_Files
#6 Updated by Gianluca Petrillo over 1 year ago
You get the complete set of dependencies with
source /cvmfs/fermilab.opensciencegrid.org/products/larsoft/setup source /cvmfs/icarus.opensciencegrid.org/products/icarus/setup ups depend icaruscode v06_54_00 -q e14:prof -f 'Linux64bit+2.6-2.12'which tells me
genie v2_12_8 and genie_xsec v2_12_6. The EventGeneratorList is "Default+CCMEC+NCMEC": the complete GENIEGen module configuration is generator: {
BeamCenter: [
-1400,
-350,
0
]
BeamDirection: [
0,
0,
1
]
BeamName: "booster"
BeamRadius: 3
DebugFlags: 0
DefinedVtxHistRange: false
DetectorLocation: "MINOS-NearDet"
Environment: [
"GPRODMODE",
"YES"
]
EventGeneratorList: "Default+CCMEC+NCMEC"
EventsPerSpill: 1
FiducialCut: "none"
FluxFiles: [
"converted_beammc_icarus_*.root"
]
FluxSearchPaths: "/pnfs/icarus/persistent/bnb_gsimple/fluxes_Oct2017/"
FluxType: "simple_flux"
GHepPrintLevel: 13
GenFlavors: [
14
]
GlobalTimeOffset: 10000
MixerBaseline: 0
MixerConfig: "none"
MonoEnergy: 2
POTPerSpill: 5e12
PassEmptySpills: false
ProductionMode: "yes"
RandomTimeOffset: 10000
SurroundingMass: 0
TopVolume: "volTPCActive"
VtxPosHistRange: [
0,
0,
0,
0,
0,
0
]
XSecTable: "gxspl-FNALsmall.xml"
module_type: "GENIEGen"
}#7 Updated by Lynn Garren over 1 year ago
- Status changed from New to Assigned
- Assignee set to Robert Hatcher
Gianluca will work with Robert.
#8 Updated by Gianluca Petrillo over 1 year ago
- File GENIEmissingProcessMapping.log.bz2 GENIEmissingProcessMapping.log.bz2 added
- Assignee deleted (
Robert Hatcher)
I have instrumented our generator module to dump all the information that LArSoft receives about the generated event, when its interaction type is kNuanceOffset (1000).
Is this information going to be enough? (needless to say, it would make it easier for us to provide the needed information).
I attach the art output from 1000 events generated with a MicroBooNE BNB configuration (GENIEmissingProcessMapping.log.bz2).
For each event, the format is:
%MSG-w GENIEmissingProcessMapping: GENIEGen:generator@ run: 1 subRun: 0 event: 1
Found an interaction that is not represented by the interaction type code in GENIE:
MCTruth record:
16 particles from neutrinos from beam
neutrino information: from charged weak current, <nuance offset>, mode: 1 (resonant), target: 1000180400 (Ar40), hit nucleon: 2112 (neutron), hit quark: 0
x=0.136345 y=0.780693 w=1.53629 Q^2=0.233431 GeV^2; theta=0.864647 rad pT=0.00456249 GeV/c
neutrino: ID=0 PDG ID 14 mass=0 GeV/c2 status=0 (initial state) weight=0 rescattered at vertex (0.102213, 1.70348, 0.426736; 0)
created at (282.29, 100.544, 331.583; 4492.84) cm by primary from ID=-1 with momentum (0.00359235, 0.0028127, 1.16783; 1.16784) GeV/c
still alive!
comes with a trajectory 0 cm long in 1 points
outgoing lepton: ID=4 PDG ID 13 mass=0.10566 GeV/c2 status=1 (stable final state) weight=0 rescattered at vertex (0.102213, 1.70348, 0.426736; 0)
created at (282.29, 100.544, 331.583; 4492.84) cm by the gods with momentum (0.165356, -0.0653761, 0.151043; 0.256114) GeV/c
still alive!
comes with a trajectory 0 cm long in 1 points
[#0] ID=0 PDG ID 14 mass=0 GeV/c2 status=0 (initial state) weight=0 rescattered at vertex (0.102213, 1.70348, 0.426736; 0)
created at (282.29, 100.544, 331.583; 4492.84) cm by primary from ID=-1 with momentum (0.00359235, 0.0028127, 1.16783; 1.16784) GeV/c
still alive!
comes with a trajectory 0 cm long in 1 points:
(282.29, 100.544, 331.583; 4492.84)
[... all other particles in the record follow]
GENIE truth record:
interaction code: 2, neutrino scattering code: 4 at (2.8229, 1.00544, 3.31583; 4.44156e-08)
probe: nu_mu with cp=(0.124363, -0.216951, 1.21553; 1.24099) hit nucleon with cp=(-0.092137, 0.167659, -0.0363944; 0.929174) GeV (not a sea quark) in target: Ar40 (Z: 18, A: 40)
event interaction weight (genie internal): 1, interaction probability: 2.89285e-13, cross section: 3.22351e-11, differential cross section: 3.9923e-10
particles after reaction, before FSI: 1 pi+, 0 pi-, 1 pi0, 0 p/pbar, 1 n/nbar
total 3 particles after reaction before FSI: 1/2 charged/neutral, 2 pions, 1 nucleons
process without charmed hadron, resonance: #1
internal (on shell) genie kinematics: Q^2: 0.233431 GeV^2 q^2: -0.233431 GeV^2, w: 1.51909 GeV^2, t: -99999 GeV^2, x: 0.136061, y: 0.760823
FShadSyst: (0, 0, 0; 0)
%MSG(I have no idea about what
FShadSyst might be)#9 Updated by Robert Hatcher over 1 year ago
- File generate_icarus_events.sh generate_icarus_events.sh added
- File find_bad_nuance.C find_bad_nuance.C added
Okay, I've done most of what I was hoping a minion would do. But now I need someone to "think" about the results to identify the condition(s) that are failing in the genie::utils::ghep::NanceReactionCode() function's logic.
I've attached to scripts to document and help with the work.
The generate_icarus_events.sh script is how I generated some events. If whomever can't access my /dune/app/users/rhatcher/fix_nuance_codes area and get a copy of icarus.4242.ghep.root, or if they'd like to make more statistics on their own, this is how they could create their own file to work on.
Once a file of events in GENIE format is available then the second script can be invoked via:
# assumes genie has been setup (ie. dunetpc, icaruscode, larsoft, whatever ...) genie -b -q find_bad_nuance.C\(\"icarus.4242.ghep.root\"\)
This script calls the GENIE code for translating event records into Nuance code and if that returns zero, then it tests whether the event is MEC (a known hole in this scheme). If it isn't MEC then it runs a copy of that routine and prints out the GENIE event record so one can look to see how it's failing; i.e. where's the hole in the NuanceReactionCode() logic.
My run of this yielded:
... [lots of printout ...] Final results: 10000 checked, 1822 MEC, 164 Resonant, 0 other
So, there are 164 Resonant interactions that don't get classified. But I don't have the time to look for the pattern(s) so here I plead for someone else to take this on. The modified copy of the code is in the script and should be fairly easy to read (also to allow you to test whether your "fix" catches everything). if you don't understand how GENIE's event records are formatted let me know and I can help with this.
Once we have a pattern, we can insert code in to nutools to "fix up" this classification until such time that the fix moves upstream to GENIE itself.
Here's an example of a failure:
*** Event #: 9857 |------------------------------------------------------------------------------------------------------------------| |GENIE GHEP Event Record [print level: 3] | |------------------------------------------------------------------------------------------------------------------| | Idx | Name | Ist | PDG | Mother | Daughter | Px | Py | Pz | E | m | |------------------------------------------------------------------------------------------------------------------| | 0 | nu_mu | 0 | 14 | -1 | -1 | 4 | 4 | 0.002 | 0.002 | 1.391 | 1.391 | 0.000 | | 1 | Ar40 | 0 | 1000180400 | -1 | -1 | 2 | 3 | 0.000 | 0.000 | 0.000 | 37.216 | 37.216 | | 2 | neutron | 11 | 2112 | 1 | -1 | 5 | 5 | -0.391 | -0.060 | -0.179 | 0.927 | **0.940 | M = 0.819 | 3 | Ar39 | 2 | 1000180390 | 1 | -1 | 14 | 14 | 0.391 | 0.060 | 0.179 | 36.288 | 36.286 | | 4 | mu- | 1 | 13 | 0 | -1 | -1 | -1 | 0.242 | -0.343 | 0.234 | 0.493 | 0.106 | P = (-0.503,0.714,-0.487) | 5 | N+(1440) | 3 | 12212 | 2 | -1 | 6 | 7 | -0.630 | 0.285 | 0.978 | 1.826 | **1.440 | M = 1.378 | 6 | Delta+ | 3 | 2214 | 5 | -1 | 8 | 9 | -0.599 | 0.294 | 0.863 | 1.645 | 1.232 | | 7 | pi0 | 14 | 111 | 5 | -1 | 10 | 10 | -0.032 | -0.009 | 0.115 | 0.180 | 0.135 | FSI = 3 | 8 | neutron | 14 | 2112 | 6 | -1 | 11 | 12 | -0.456 | -0.002 | 0.647 | 1.229 | 0.940 | FSI = 4 | 9 | pi+ | 14 | 211 | 6 | -1 | 13 | 13 | -0.142 | 0.295 | 0.216 | 0.417 | 0.140 | FSI = 3 | 10 | pi0 | 1 | 111 | 7 | -1 | -1 | -1 | 0.094 | 0.010 | 0.076 | 0.181 | 0.135 | | 11 | neutron | 1 | 2112 | 8 | -1 | -1 | -1 | -0.499 | 0.112 | 0.538 | 1.197 | 0.940 | | 12 | proton | 1 | 2212 | 8 | -1 | -1 | -1 | 0.125 | -0.045 | 0.082 | 0.951 | 0.938 | | 13 | pi+ | 1 | 211 | 9 | -1 | -1 | -1 | 0.057 | 0.185 | 0.343 | 0.418 | 0.140 | | 14 | HadrBlob | 15 | 2000000002 | 3 | -1 | -1 | -1 | -0.017 | 0.084 | 0.118 | 35.367 | **0.000 | M = 35.366 |------------------------------------------------------------------------------------------------------------------| | Fin-Init: | 0.000 | -0.000 | -0.000 | -0.000 | | |------------------------------------------------------------------------------------------------------------------| | Vertex: nu_mu @ (x = 1.29272 m, y = 1.43351 m, z = 15.72722 m, t = 1.191694e-07 s) | |------------------------------------------------------------------------------------------------------------------| | Err flag [bits:15->0] : 0000000000000000 | 1st set: none | | Err mask [bits:15->0] : 1111111111111111 | Is unphysical: NO | Accepted: YES | |------------------------------------------------------------------------------------------------------------------| | sig(Ev) = 7.17431e-39 cm^2 | d2sig(W,Q2;E)/dWdQ2 = 3.35100e-38 cm^2/GeV^3 | Weight = 1.00000 | |------------------------------------------------------------------------------------------------------------------| -------------------------------------------------------------------------------------------------------------- GENIE Interaction Summary -------------------------------------------------------------------------------------------------------------- [-] [Init-State] |--> probe : PDG-code = 14 (nu_mu) |--> nucl. target : Z = 18, A = 40, PDG-Code = 1000180400 (Ar40) |--> hit nucleon : PDC-Code = 2112 (neutron) |--> hit quark : no set |--> probe 4P : (E = 1.391272, Px = 0.002316, Py = 0.002085, Pz = 1.391269) |--> target 4P : (E = 37.215526, Px = 0.000000, Py = 0.000000, Pz = 0.000000) |--> nucleon 4P : (E = 0.927105, Px = -0.390514, Py = -0.060287, Pz = -0.178662) [-] [Process-Info] |--> Interaction : Weak[CC] |--> Scattering : RES [-] [Kinematics] |--> *Selected* Bjorken x = 0.365773 |--> *Selected* Inelasticity y = 0.628192 |--> *Selected* Momentum transfer Q2 (>0) = 0.707456 |--> *Selected* Hadronic invariant mass W = 1.377735 [-] [Exclusive Process Info] |--> charm prod. : false |--> strange prod. : false |--> f/s nucleons : N(p) = 0 N(n) = 0 |--> f/s pions : N(pi^0) = 0 N(pi^+) = 0 N(pi^-) = 0 |--> resonance : P11(1440) -------------------------------------------------------------------------------------------------------------- FOR CONVENIENCE: IsResonant 1 np 0 nn 1 npip 1 npim 0 npi0 1 IsWeakNC() 0 IsWeakCC() 1 IsNuP() 0 IsNuN() 1 IsNuBarP() 0 IsNuBarN() 0
#10 Updated by Robert Hatcher over 1 year ago
- File nuance-channels.pdf nuance-channels.pdf added
Also attaching a PDF of the Nuance channel encoding as I got it from Sam Zeller
#11 Updated by Lynn Garren over 1 year ago
Andrea, would you take a look and see if you can identify a pattern in these events?
#12 Updated by Katherine Lato about 1 year ago
- Status changed from Assigned to Work in progress
- Assignee set to Andrea Falcone
#13 Updated by Andrea Falcone about 1 year ago
After a first look at the resonant interactions that don't get classified by Robert routine:
- they seem to be "normal" interaction without strange stuff.
- the resonance particle (delta, eta, lambda) decay in pi or gamma -> can the "strange" number of pions in final stat led to the non classification? Strange
I'll look more in details