Potential issue with GENIE on atmospheric events.

[candreop]

Reported by Xianguo Lu Xianguo.Lu@warwick.ac.uk, and his student Qiyu Yan yanqiyu17@mails.ucas.edu.cn

They report that gevgen and gevgen_atmo generates different output for same flux input.

Gevgen and gevgen_atmo.pdf

[nusense]

I haven't checked but by chance is the histogram one with variable width bins? In which case I think they may need to use the WIDTH indicator depending on how the bin content are normalized.

 A 1-D ROOT histogram (TH1D):
                 The general syntax is `-f /full/path/file.root,object_name<,[NO]WIDTH>'
                 by default variable bin histograms are multiplied by bin width
                 if ,NOWIDTH then they are not
                 if ,WIDTH then they are alway multiplied by bin width

I'm not setup for testing this at the moment, but if you @candreop are, or they would like to try it someone should.

[candreop]

Not setup for testing at the moment either, but will forward the above to Xianguo et al. I should be able to test next week, if this persists. Thanks @nusense

[nusense]

Yeah, I think the histogram file is variable binned, and filled with E^{-2.7}; but it need to be filled with E^{-2.7)*{bin-width} for ROOT's sampling to correctly work, so needs the WIDTH flag to adjust for this. Ala, gevgen -p 14 ...

-f /path/to/honda.root,fluxnumu,WIDTH

or

-f /path/to/flux1d.root,nu_mu_px,WIDTH

[karuboniru]

Hi, I am Qiyu Yan. Thank you for your following up!

I tried to add ,WIDTH to every gevgen command's flux input, and result seems to be exact as before. (I suppose that by not specifying the option, default for a variable bin hist will result as if with ,WIDTH?)

by default variable bin histograms are multiplied by bin width.

Compare plots with .WIDTH in commands

[nusense]

It's also not clear whether this was done with R-3_00_06 or R-3_02_00.

I ~think~ am sure the ,WIDTH stuff is only available in R-3_02 and newer.

[karuboniru]

I am on R-3_02_00 tag.

[nusense]

Oh, drat. I was kind of hoping that you were using R-3_00_06 and this was fixed in R-3_02_00.

Here's another "easy" suggestion to help identify the problem: what if you used ,NOWIDTH in gevgen. Does that bring them in alignment with gevgen_atmo? (If so, then it probably means that gevgen_atmo is doing things wrong and not accounting for the variable bin widths)

[karuboniru]

with NOWIDTH.pdf

With NOWIDTH option the flux changed a lot thus changed the results a lot (as the difference between gevgen_atmo and gevgen is so small), so no luck.

[nusense]

Double drat. Okay, this is going to take some more investigation.

[nusense]

@karuboniru Could you supply the honda.dummy.txt, flux_14.root, honda.root and flux1d.root files so one doesn't have to go through the steps of patching for the extraction. I'd like to try running with an unmodified existing code base.

I've got my own XML and ROOT spline files (those shouldn't matter, hopefully)

[karuboniru]

files.tar.gz Yes, they are in the tarball above.

[nusense]

Okay, I ran case (1) and (2) for myself. What I noticed is that we're making an energy restriction [0.0,20.0] where the 20 does not correspond directly to a flux bin boundary. If I plot the nu energy distribution for the atmospherics it cuts off exactly at 20 GeV, but the gevgen distribution seems to be cut closer to 19 GeV (which isn't a boundary either since the bins are ... 17.783, 19.953, 22.387 ...). But it does mean that two have a slightly different energy distribution.

[nusense]

[nusense]

Hmm, the bin boundaries on honda.root's fluxnumu don't seem to correspond with the ones in honda.dummy.txt with bin boundaries ... 16.78823, 18.83685, 21.13469, ... that middle one seemingly corresponds with the cut-off plot above.

The bin skew appears to be due to

       energy.push_back(energy_temp / std::pow(10., 1 / 40.));

in honda2root.cxx, where std::pow(10., 1 / 40.) -> (double) 1.0592537

I'm not sure I understand how to get from three .gtrac.root files to the 3 curves on the plot. I though the kaon program was going to do that, but it apparently outputs a single TH1D ... which I also don't understand what it represents (I was expecting N objects each corresponding to an input .gtrac.root file).

[karuboniru]

If I plot the nu energy distribution for the atmospherics it cuts off exactly at 20 GeV, but the gevgen distribution seems to be cut closer to 19 GeV

Yeah, I never checked that, which means the boundary handling is somehow different between 2 commands? This can cause the difference I see. I will rerun the whole simulation/analysis and see what if I only select events under some certain threshold.

the bin boundaries on honda.root's fluxnumu don't seem to correspond with the ones in honda.dummy.txt with bin boundaries

I tried to do the shift to make data points given in .txt lays in the middle of the bin (instead of on edge, when looking from log scale), I also had checked that the shift don't change the result so much.

Also, the step of extracting internal histogram from gevgen_atmo step is to rule out any problem that may be caused from binning.

it apparently outputs a single TH1D

The plot is done by drawing all 3 single TH1Ds in same plot, I could provide code doing the plotting.

which I also don't understand what it represents

Store per particle $K^+$ kinetic energy in a histogram and normalize it as if normalizing cross section plot.

[nusense]

Okay, so GAtmoFlux.cxx imposes hard cuts based on the range set by the user, where sometime getting flux energies out of this range and then rejecting them (not passing them on for event generation) before going to generate the next flux ray. Thus it will always exactly represent the chosen -E emin,emax range.

On the other hand, gevgen when given a root file with a histogram does the following:

    // Copy in the flux histogram from the root file and remove bins outside the emin,emax range
    spectrum = (TH1D*)hst->Clone();
    spectrum->SetNameTitle("spectrum","neutrino_flux");
    spectrum->SetDirectory(0);
    for(int ibin = 1; ibin <= hst->GetNbinsX(); ibin++) {
      if(hst->GetBinLowEdge(ibin) + hst->GetBinWidth(ibin) > emax ||
         hst->GetBinLowEdge(ibin) < emin) {
        spectrum->SetBinContent(ibin, 0);
      }
    }

So if you don't pick a -e emin,emax corresponding to bin boundaries you can have the bin containing the emax zero'ed out -- effectively cutting off the spectrum a the low edge bin that contain emax. Thus it is best to pick a bin boundary for the -e flag.

[karuboniru]

Seems that the high energy part don't explain this. I tried to place a $E_{\nu} < 10~\mathrm{GeV}$ cut on my analysis code. And the difference we see still exist. So problem doesn't come from the edge handling.

Enu<10GeV cut.pdf

[nusense]

But did you change the binning back to match? Rebinning everything down by 95% but keeping the contents the same (as is done to generate honda.root 's fluxnumu ) is going to soften the spectrum.

I'm just having a very, very hard time coming up with a hypothesis for how using the atmospheric flux driver (vs. a supposedly identical histogram) can change the underlying physics. Fundamentally they both get configured identically, so it's got to be either the cuts or the actual spectrum used. I guess one other possibility is smoothing of the steeply falling spectrum. I'd have to look to see if atmospheric flux driver does something special to account for that.

[karuboniru]

Add_aligned_bin.pdf

I tried to remove / std::pow(10., 1 / 40.) part when generating flux histogram. The result did not change.

[nusense]

Okay, I'm still not at the level of calculating / graphing sigma/dp_K but I did exercise the flux drivers (with aligned bins, and cutoffs at bin boundaries) and to a very significant degree the two give exactly the same spectrum.

I really am at a loss of what else there can be. The underlying physics should all be the same (configured via a tune). Events is flux * cross-section.

[karuboniru]

I really am at a loss of what else there can be. The underlying physics should all be the same (configured via a tune). Events is flux * cross-section.

This is the exact questions I am having, it is hard to imagine whatever could go different except flux driver, but I can't find anything wrong with flux driver.

I made a guide containing whatever tool and steps to generate exact plots I am having.

[nusense]

Unfortunately I can't get kaon_plot 's genie_kaon.cxx to build with ROOT v6_22_08 and C++17 [sigh]. Among other things it fails on:

...
/genie/app/rhatcher/genie-kaon-spectrum-simple-master/kaon_plot/genie_kaon.cxx:109:37:   required from here
/cvmfs/larsoft.opensciencegrid.org/products/gcc/v9_3_0/Linux64bit+3.10-2.17/include/c++/9.3.0/ext/new_allocator.h:145:20: 
error: no matching function for call to ‘ROOT::RDF::RResultPtr::RResultPtr(ROOT::RDF::RResultPtr)’
  145 |  noexcept(noexcept(::new((void *)__p)
      |                    ^~~~~~~~~~~~~~~~~~
  146 |        _Up(std::forward<_Args>(__args)...)))
      |        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
...
/cvmfs/larsoft.opensciencegrid.org/products/root/v6_22_08d/Linux64bit+3.10-2.17-e20-p392-debug/include/ROOT/RResultPtr.hxx:135:4: note: candidate: ‘ROOT::RDF::RResultPtr::RResultPtr(std::shared_ptr<_Tp>, ROOT::Detail::RDF::RLoopManager*, std::shared_ptr) [with T = TH1]’
  135 |    RResultPtr(std::shared_ptr objPtr, RDFDetail::RLoopManager *lm,
      |    ^~~~~~~~~~
/cvmfs/larsoft.opensciencegrid.org/products/root/v6_22_08d/Linux64bit+3.10-2.17-e20-p392-debug/include/ROOT/
RResultPtr.hxx:135:4: note:   candidate expects 3 arguments, 1 provided

Also how many jobs of 200,000 events goes into your samples for atmo / gevgen each to get the statistics you're showing?

[karuboniru]

Unfortunately I can't get kaon_plot 's genie_kaon.cxx to build with ROOT v6_22_08 and C++17 [sigh]. Among other things it fails on:

I have done some modification and now it should work on your exact GCC and ROOT version.

Also how many jobs of 200,000 events goes into your samples for atmo / gevgen each to get the statistics you're showing?

I am running 32 * 200,000 for each case.

[nusense]

Well I had to run the jobs on the grid, but I got 32 jobs of gevgen_atmo and gevgen (with honda.root ) both generated Enu 0,22.387 (ie bin boundary). which shows the same atmo enhancement below 0.3 GeV.

[karuboniru]

Yes, seems we are at same place, I am also very confused when guessing where the problem could be apart from different flux driver used.

Inspecting GENIE code is hard for me since I don't have any knowledge about the whole picture. So I think I need some help from you to find where where the real different is.

[nusense]

Trying to record here what differences there are between the two cases.

Atmos events are isotropic in direction; gevgen events are all along z axis (have verified that atmospheric come from different directions, but not that it is isotropic; verified that gevgen are p_z=E). This might matter if the "analysis" code were direction dependent. At first glance, it doesn't appear to be but I state this difference between the samples for the record.

gevgen calls mcj_driver->SetUnphysEventMask(*RunOpt::Instance()->UnphysEventMask()); (though no mask is input to gevgen) but the gAtmoEvGen doesn't seem to. But I think the default should be to reject unphysical events in both cases, but this initialization might be missing. A gevdump of some of the atmospheric files didn't turn up any flagged as unphysical.

Atmospheric flux driver was run with Rl``Rt forced to 1000 in gAtmoEvGen but left at 0,0 for the flux exercising shown above. atmospheric_flux_driver->SetRadii(gOptRL, gOptRT). I think this would only matter if we were running with a real geometry with "extent" rather than a point geometry of just C12.

The code for initialization of the tune is essentially identical for the two. The initialization of the point geometry driver is essentially identical for the two.

[nusense]

Also, the GHAKKMAtmoFlux driver expects/assumes 20 cosTheta bins (over -1,1), 12 phi bins and 101 bins in energy. The honda.dummy.txt file doesn't have but one cosTheta / one phi bin's worth of data (101 values). So when thrown it only generates in a range cosTheta [-1,-0.9] and phi [-pi,-pi+2pi/12]. So, the atmos flux is almost the opposite direction as the histogram flux (which has pz = +E).

[nusense]

[nusense]

Okay, so I made a modified gevgen_atmo that did two things: First, explicitly set the unphysical mask from RunOpt Secondly, had the option (which I used) to force the flux along the pz = +E axis. One (or both?) of those made the atmospheric modified curve look like gevgen This is the plot I got:

[karuboniru]

So the root cause gets narrowed down to only 2 reason? One is unphysical mask and another is angular distribution in one bin?

Is is hard for me to imagine what the direction will change the $K^+$ energy distribution. So maybe looking into unphysical mask is a good start?

[nusense]

The next trial was to run with the unphysical mask enabled and back to the flux from (roughly) -pz just because that was easier with the modified atmo executable. If it was the lack of (re) setting the unphysical mask in the baseline gevgen_atmo then the curve should stay at the gevgen level when run with the original -pz flux. Alas, it didn't:

So, that means that either the "analysis" is somehow sensitive to the flux direction, or GENIE has an internal inconsistency with respect to flux isotropy. The next thing I can try is modifying gevgen to allow setting a flux direction other than along+z (the driver supports it, there's just no command line interface for setting it ... but I think gevgen_fnal supports it so I might either try that or carry over the code from that to plain ol' gevgen

[nusense]

Next I modified gevgen to accept flags so one could specify the GCylinderTH1Flux direction, radius, start spot. Then I ran with direction (0,0,D)where D=+1 for "pluspz" and D=-1 for "minuspz" keeping radius=-1, spot (0,0,0). I would have expected the "pluspz" to exactly reproduce the original gevgen running as it should have been the same conditions, but it didn't (in fact I got two run failures due to the assert(s>0) failure; I also got two runs failed for "minuspz" running for the same assert). Nevertheless, the pattern matched: higher rate at low T_K for "minuspz" running relative to "pluspz" running.

I looked at the "analysis" code. I don't speak ROOT's RDataFrame language, but best I can work out I don't see anything there that should depend on the flux's directionality. I haven't worked out in my mind the normalization (something involving the rate and the underlying spline cross-sections) but again I don't see anything that should depend on direction. Which leads to the uncomfortable conclusion that there's something in GENIE itself that leads to a rate of Kaon production (as a function of Kaon KE) that depends on the flux direction. This is super subtle, and pulling it out depends on lots of statistics (on order 6.4 million events, with significant wiggle in the curve), so I'm wondering if there's something more common / general that can be identified that also has a behavior that depends on flux direction.

[nusense]

Just to put to rest any concern about the ghep.root to rootracker format conversion and use, or the "analysis" cross-section normalization: I've plotted the bare KE of kaons in the original "minuspz" / "pluspz" .ghep.root files, the only cut being that the event be "Weak[CC]", normalized solely to the number of events.

[nusense]

The E^{-2.7} spectrum is heavily weighted towards lower energies. I wondered if this behavior held for higher energies. Here's the same flux restricted from [17.783,22.387] GeV.

[karuboniru]

I don't speak ROOT's RDataFrame language, but best I can work out I don't see anything there that should depend on the flux's directionality.

There should not be any dependencies on direction in analysis code. Although complexed at first look, it just reading rootracker file, fill 4-momenta to TLorentzVectors and calculate kinetic energy by TLorentzVector's E()-M() method.

[karuboniru]

Any updates?

[karuboniru]

More investigation results: the problem comes from leading kaon production at nucleon level (no fermi motion or FSI involved) kaon-update.pdf

[karuboniru]

GENIE direction dependence in DIS channel.pdf

Some new observation: The direction dependence exists on single nucleon level in DIS channel, and also affect $\pi^+$ production. For details please refer to the PDF file attached.

[candreop]

That's very interesting @karuboniru. It seems to suggest that there is something in the DIS generator specifically (you have narrowed it to very few possible event generation modules) that has a built-in assumption about a direction (+z?) and doesn't actually read the direction of the incoming neutrino (or something else of that similar flavour). Where that may happen still evades us. Can I suggest a couple more checks to try and diagnose this:

• Can we take gevgen_atmo and the atmospheric flux drivers out of the test? For example, if we change the simple GMonoEnergeticFlux driver used in gevgen, so that the flux is along +x or +y instead of +z (change line 126 in GMonoEnergeticFlux.cxx and remember to recompile) what do we get for dsigma/dT for pions and Kaons? If it changes from the previous gevgen result then (without doubt) it is a problem a hidden +z assumption in the DIS event generation modules
• If we use the same random number seed for a run with a simple flux in +x direction and for a run with a simple flux in +z direction, then we may be able to relate the events in the 2 runs? For example, if there is a missing or extra boost or rotation somewhere, we may be able to understand what that boost or rotation is.

[karuboniru]

6-dirs.pdf Tried to run on 6 different directions[^1], it seems that it is due to some hidden +z (or -z) assumption.

And I looked at

• $T_{\mu^-}$ distribution here, where 2 directions along z axis agreed,
• $T_{p}$ distribution here, where the problem exist.

This suggests that the problem have something to do with hadronic system and not with leptonic kinematics. This could be somewhere around PYTHIA where hadronization happen?

About the correlation check, I may need sometime to find out a proper way to do this.

---

---

Robert W Hatcher's code

``` diff --git a/src/Apps/gEvGen.cxx b/src/Apps/gEvGen.cxx index 67c6338..5e99296 100644 --- a/src/Apps/gEvGen.cxx +++ b/src/Apps/gEvGen.cxx @@ -77,6 +77,10 @@ Syntax: by default variable bin histograms are multiplied by bin width if ,NOWIDTH then they are not if ,WIDTH then they are alway multiplied by bin width + -F + Specifies direction,size,start point of histogram flux + dircosx,dircosy,dircosz,Radius,spotx,spoty,spotz + -o Specifies the name of the output file events will be saved in. -w @@ -217,6 +221,7 @@ int gOptNuPdgCode; // neutrino PDG code map gOptTgtMix; // target mix (each with its relative weight) Long_t gOptRunNu; // run number string gOptFlux; // +string gOptFluxFactors; // bool gOptWeighted; // bool gOptUsingFluxOrTgtMix = false; long int gOptRanSeed; // random number seed @@ -572,11 +577,30 @@ GFluxI * TH1FluxDriver(void) f.Close(); TVector3 bdir (0,0,1); + double radius = -1; TVector3 bspot(0,0,0); + if ( gOptFluxFactors != "" ) { + vector fv = utils::str::Split(gOptFluxFactors,","); + if ( fv.size() >= 3 ) { + bdir.SetX(atoi(fv[0].c_str())); + bdir.SetY(atoi(fv[1].c_str())); + bdir.SetZ(atoi(fv[2].c_str())); + if ( fv.size() >= 4 ) { + radius = atoi(fv[3].c_str()); + if ( fv.size() >= 7 ) { + bspot.SetX(atoi(fv[4].c_str())); + bspot.SetY(atoi(fv[5].c_str())); + bspot.SetZ(atoi(fv[6].c_str())); + } + } + } + } + flux->SetNuDirection (bdir); + flux->SetTransverseRadius (radius); flux->SetBeamSpot (bspot); - flux->SetTransverseRadius (-1); + flux->AddEnergySpectrum (gOptNuPdgCode, spectrum); GFluxI * flux_driver = dynamic_cast(flux); @@ -644,6 +668,10 @@ void GetCommandLineArgs(int argc, char ** argv) gOptFlux = parser.ArgAsString('f'); using_flux = true; } + if (parser.OptionExists('F') ) { + LOG("gevgen", pINFO) << "Reading flux factors"; + gOptFluxFactors = parser.ArgAsString('F'); + } if(parser.OptionExists('s')) { LOG("gevgen", pWARN) @@ -779,7 +807,7 @@ void GetCommandLineArgs(int argc, char ** argv) << "No input cross-section spline file"; } LOG("gevgen", pNOTICE) - << "Flux: " << gOptFlux; + << "Flux: " << gOptFlux << " factors " << gOptFluxFactors; LOG("gevgen", pNOTICE) << "Generate weighted events? " << gOptWeighted; if(gOptNuEnergyRange>0) { ```

[^1]: This is done with the patch by Robert W Hatcher sent to me via e-mail

[candreop]

Thanks @karuboniru. A potential issue I see in the PYTHIA6 and PYTHIA8 interface code (but not in the empirical low mass AGKY or in resonances) is summarized below.

Take the PYTHIA6 code, for example: https://github.com/GENIE-MC/Generator/blob/master/src/Physics/Hadronization/Pythia6Hadro2019.cxx

The hadronic system is generated at a special reference frame (CM frame with z aligned with the direction of the total hadronic momentum). Then, the system is boosted back to a rotated Lab frame (line 159) and then rotated to the proper Lab frame (line 160). The rotation is determined by the unitvq TVector3 in line 116 (TVector3 unitvq = p4Had.Vect().Unit()). That hadronic system 4-momentum (p4Had) is obtained in line 113 (TLorentzVector p4Had = kinematics.HadSystP4())

A problem is that the Kinematics object does not have the full event information to determine the hadronic system 4-momentum. For example, it does not have the full neutrino or hit nucleon 4-momenta. It is just a little bag we use to carry information around, so that we don't have to extract it from the full event record again and again. If we want to extract info from kinematics.HadSystP4(), then we must set it first. However, running a grep command, I cannot find an instance where it is set (in the code run before PYTHIA6).

I am wondering if the following solves the problem. Edit src/Physics/Hadronization/Pythia6Hadro2019.cxx, comment out line 113 (TLorentzVector p4Had = kinematics.HadSystP4();) and replace it with the following: TLorentzVector p4v ( event->Probe()->P4() ); TLorentzVector p4N ( event->HitNucleon()->P4() ); TLorentzVector p4l ( * event->FinalStatePrimaryLepton()->P4() ); TLorentzVector p4Had = p4v + p4N - p4l;

Does that improve the comparisons?

[karuboniru]

The result with your patch, seems the result did not change.

---

I am also looking at the correlation from samples before (without your modification, using +z, -z, +x and -x sample).

I found a set of events below (just by looking at genie-mcjob-0.status generated during simulation), if you compare them by ignoring difference in quadrant, you can see that they are exact event until final state hadron appears (namely Idx=4 nucleon here). This agreed with the initial guess that problem is around hadronization and have nothing to do with lepton kinematics. Should I look at more such comparison and see if there is some underlying pattern?

And if compare +x and -x, we can see the kinetic energy of final state particles are not still the same. But we can't see statistically difference in https://github.com/GENIE-MC/Generator/issues/226#issuecomment-1454746259 could it be

• Different direction in x axis do involve some buggy boost, but it don't change statistical distribution for final state particles? OR
• PYTHIA have some internal random engine that is not under control of GENIE (unlikely, but I am unfamiliar with PYTHIA coding, just an uneducated guess)

---

+z

|------------------------------------------------------------------------------------------------------------------|
|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 |   2 |   2 |   0.000 |   0.000 |   6.733 |   6.733 |   0.000 | 
|   1 |        proton |   0 |       2212 |  -1 |  -1 |   3 |   3 |   0.000 |   0.000 |   0.000 |   0.938 |   0.938 | 
|   2 |           mu- |   1 |         13 |   0 |  -1 |  -1 |  -1 |   1.337 |  -0.033 |   2.886 |   3.183 |   0.106 | P = (-0.420,0.010,-0.907)
|   3 |      HadrSyst |  12 | 2000000001 |   1 |  -1 |   4 |   6 |  -1.337 |   0.033 |   3.847 |   4.488 | **0.000 | M = 1.887 
|   4 |       neutron |   1 |       2112 |   3 |  -1 |  -1 |  -1 |  -1.020 |  -0.340 |   2.182 |   2.607 |   0.940 | 
|   5 |           pi+ |   1 |        211 |   3 |  -1 |  -1 |  -1 |  -0.427 |   0.350 |   0.715 |   0.915 |   0.140 | 
|   6 |           pi+ |   1 |        211 |   3 |  -1 |  -1 |  -1 |   0.110 |   0.023 |   0.950 |   0.967 |   0.140 | 
|------------------------------------------------------------------------------------------------------------------|
|       Fin-Init:                                                |   0.000 |   0.000 |  -0.000 |  -0.000 |         | 
|------------------------------------------------------------------------------------------------------------------|
|       Vertex:          nu_mu @ (x =     0.00000 m, y =     0.00000 m, z =     0.00000 m, t =    0.000000e+00 s)  |
|------------------------------------------------------------------------------------------------------------------|
| Err flag [bits:15->0] : 0000000000000000    |  1st set:                                                     none | 
| Err mask [bits:15->0] : 1111111111111111    |  Is unphysical:    NO |   Accepted:   YES                          |
|------------------------------------------------------------------------------------------------------------------|
| sig(Ev) =       1.96305e-38 cm^2  | d2sig(x,y;E)/dxdy =       1.54265e-38 cm^2       | Weight =          1.00000 |
|------------------------------------------------------------------------------------------------------------------|

-z

|------------------------------------------------------------------------------------------------------------------|
|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 |   2 |   2 |   0.000 |   0.000 |  -6.733 |   6.733 |   0.000 | 
|   1 |        proton |   0 |       2212 |  -1 |  -1 |   3 |   3 |   0.000 |   0.000 |   0.000 |   0.938 |   0.938 | 
|   2 |           mu- |   1 |         13 |   0 |  -1 |  -1 |  -1 |  -1.337 |  -0.033 |  -2.886 |   3.183 |   0.106 | P = (0.420,0.010,0.907)
|   3 |      HadrSyst |  12 | 2000000001 |   1 |  -1 |   4 |   6 |   1.337 |   0.033 |  -3.847 |   4.488 | **0.000 | M = 1.887 
|   4 |       neutron |   1 |       2112 |   3 |  -1 |  -1 |  -1 |   0.449 |  -0.342 |  -2.043 |   2.319 |   0.940 | 
|   5 |           pi+ |   1 |        211 |   3 |  -1 |  -1 |  -1 |   0.154 |   0.351 |  -0.956 |   1.040 |   0.140 | 
|   6 |           pi+ |   1 |        211 |   3 |  -1 |  -1 |  -1 |   0.735 |   0.024 |  -0.847 |   1.130 |   0.140 | 
|------------------------------------------------------------------------------------------------------------------|
|       Fin-Init:                                                |  -0.000 |   0.000 |   0.000 |  -0.000 |         | 
|------------------------------------------------------------------------------------------------------------------|
|       Vertex:          nu_mu @ (x =     0.00000 m, y =     0.00000 m, z =     0.00000 m, t =    0.000000e+00 s)  |
|------------------------------------------------------------------------------------------------------------------|
| Err flag [bits:15->0] : 0000000000000000    |  1st set:                                                     none | 
| Err mask [bits:15->0] : 1111111111111111    |  Is unphysical:    NO |   Accepted:   YES                          |
|------------------------------------------------------------------------------------------------------------------|
| sig(Ev) =       1.96305e-38 cm^2  | d2sig(x,y;E)/dxdy =       1.54265e-38 cm^2       | Weight =          1.00000 |
|------------------------------------------------------------------------------------------------------------------|

+x

|------------------------------------------------------------------------------------------------------------------|
|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 |   2 |   2 |   6.733 |   0.000 |   0.000 |   6.733 |   0.000 | 
|   1 |        proton |   0 |       2212 |  -1 |  -1 |   3 |   3 |   0.000 |   0.000 |   0.000 |   0.938 |   0.938 | 
|   2 |           mu- |   1 |         13 |   0 |  -1 |  -1 |  -1 |   2.886 |  -0.033 |  -1.337 |   3.183 |   0.106 | P = (-0.907,0.010,0.420)
|   3 |      HadrSyst |  12 | 2000000001 |   1 |  -1 |   4 |   6 |   3.847 |   0.033 |   1.337 |   4.488 | **0.000 | M = 1.887 
|   4 |       neutron |   1 |       2112 |   3 |  -1 |  -1 |  -1 |   1.535 |  -0.344 |   0.605 |   1.930 |   0.940 | 
|   5 |           pi+ |   1 |        211 |   3 |  -1 |  -1 |  -1 |   0.466 |   0.348 |   0.131 |   0.613 |   0.140 | 
|   6 |           pi+ |   1 |        211 |   3 |  -1 |  -1 |  -1 |   1.845 |   0.028 |   0.601 |   1.946 |   0.140 | 
|------------------------------------------------------------------------------------------------------------------|
|       Fin-Init:                                                |  -0.000 |   0.000 |  -0.000 |  -0.000 |         | 
|------------------------------------------------------------------------------------------------------------------|
|       Vertex:          nu_mu @ (x =     0.00000 m, y =     0.00000 m, z =     0.00000 m, t =    0.000000e+00 s)  |
|------------------------------------------------------------------------------------------------------------------|
| Err flag [bits:15->0] : 0000000000000000    |  1st set:                                                     none | 
| Err mask [bits:15->0] : 1111111111111111    |  Is unphysical:    NO |   Accepted:   YES                          |
|------------------------------------------------------------------------------------------------------------------|
| sig(Ev) =       1.96305e-38 cm^2  | d2sig(x,y;E)/dxdy =       1.54265e-38 cm^2       | Weight =          1.00000 |
|------------------------------------------------------------------------------------------------------------------|

-x

|------------------------------------------------------------------------------------------------------------------|
|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 |   2 |   2 |  -6.733 |   0.000 |   0.000 |   6.733 |   0.000 | 
|   1 |        proton |   0 |       2212 |  -1 |  -1 |   3 |   3 |   0.000 |   0.000 |   0.000 |   0.938 |   0.938 | 
|   2 |           mu- |   1 |         13 |   0 |  -1 |  -1 |  -1 |  -2.886 |   0.033 |  -1.337 |   3.183 |   0.106 | P = (0.907,-0.010,0.420)
|   3 |      HadrSyst |  12 | 2000000001 |   1 |  -1 |   4 |   6 |  -3.847 |  -0.033 |   1.337 |   4.488 | **0.000 | M = 1.887 
|   4 |       neutron |   1 |       2112 |   3 |  -1 |  -1 |  -1 |  -2.682 |  -0.380 |   0.825 |   2.983 |   0.940 | 
|   5 |           pi+ |   1 |        211 |   3 |  -1 |  -1 |  -1 |  -1.128 |   0.335 |   0.257 |   1.212 |   0.140 | 
|   6 |           pi+ |   1 |        211 |   3 |  -1 |  -1 |  -1 |  -0.037 |   0.012 |   0.255 |   0.293 |   0.140 | 
|------------------------------------------------------------------------------------------------------------------|
|       Fin-Init:                                                |   0.000 |   0.000 |  -0.000 |  -0.000 |         | 
|------------------------------------------------------------------------------------------------------------------|
|       Vertex:          nu_mu @ (x =     0.00000 m, y =     0.00000 m, z =     0.00000 m, t =    0.000000e+00 s)  |
|------------------------------------------------------------------------------------------------------------------|
| Err flag [bits:15->0] : 0000000000000000    |  1st set:                                                     none | 
| Err mask [bits:15->0] : 1111111111111111    |  Is unphysical:    NO |   Accepted:   YES                          |
|------------------------------------------------------------------------------------------------------------------|
| sig(Ev) =       1.96305e-38 cm^2  | d2sig(x,y;E)/dxdy =       1.54265e-38 cm^2       | Weight =          1.00000 |
|------------------------------------------------------------------------------------------------------------------|

[karuboniru]

Hi, I think I am close to the answer: I tried some event by event check and found that some events are affected by the bug and others not. Then I checked that problem events (I only checked a few) are from AGKYLowW2019 (sorry for confusion by guesses before, the underlying problem is not pythia calling)

Then I checked the boost in AGKYLowW2019, seems it is doing something wrong. From my understanding, the before-boost hadron system from hadronization should first boost along z axis then rotate by the direction of hadronic system momentum. This is because hadronization code itself calculate things in center-of-mass frame of hadronic system and assume z axis as momentum transfer direction. The first boost boost to rotated LAB frame which remains z axis as momentum transfer direction. Then the rotation gives final momentum in LAB frame.

But in AGKYLowW2019 there is only one boost directly from hadronic system frame to lab frame, which break some angular distribution against momentum transfer direction hence caused the difference we see.

So I modified code and checked some event by event result, energy matched exactly but I am seeing problem like angle against initial neutron changed in event-by-event comparison (but angle against HadrSyst, or momentum transfer $q$ remains same in event-by-event basis in 2 neutrino direction, I believe this is expected and this won't affect final distribution of momentum and angle? Please correct me if I am wrong).

I haven’t looked into distributions yet since MC needs some time, but I will update as soon as plots are ready.

From bc99f3605d4df9e5111d41e864544ff29fecb23b Mon Sep 17 00:00:00 2001
From: Qiyu Yan <yanqiyu17@mails.ucas.ac.cn>
Date: Thu, 9 Mar 2023 00:41:12 +0800
Subject: [PATCH] fix a dir bug

---
 src/Physics/Hadronization/AGKYLowW2019.cxx | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/Physics/Hadronization/AGKYLowW2019.cxx b/src/Physics/Hadronization/AGKYLowW2019.cxx
index 054208f71..3c36816ba 100644
--- a/src/Physics/Hadronization/AGKYLowW2019.cxx
+++ b/src/Physics/Hadronization/AGKYLowW2019.cxx
@@ -121,7 +121,8 @@ void AGKYLowW2019::ProcessEventRecord(GHepRecord * event) const {
   // retrieve the hadronic blob lorentz boost
   // Because Hadronize() returned particles not in the LAB reference frame
   const TLorentzVector * had_syst = event -> Particle(mom) -> P4() ;
-  TVector3 boost = had_syst -> BoostVector() ;
+  TVector3 beta = TVector3(0,0,had_syst->P()/had_syst->E());
+  TVector3 unitvq = had_syst->Vect().Unit();

   GHepParticle * neutrino  = event->Probe();
   const TLorentzVector & vtx = *(neutrino->X4());
@@ -133,8 +134,8 @@ void AGKYLowW2019::ProcessEventRecord(GHepRecord * event) const {
     int pdgc = particle -> Pdg() ;

     //  bring the particle in the LAB reference frame
-    particle -> P4() -> Boost( boost ) ;
-
+    particle -> P4() -> Boost (beta);
+    particle -> P4() -> RotateUz(unitvq);
     // set the proper status according to a number of things:
     // interaction on a nucleaus or nucleon, particle type
     GHepStatus_t ist = ( particle -> Status() ==1 ) ? istfin : kIStDISPreFragmHadronicState;
-- 
2.39.1

[luxianguo]

Reported by Xianguo Lu Xianguo.Lu@warwick.ac.uk, and his student Qiyu Yan yanqiyu17@mails.ucas.edu.cn

They report that gevgen and gevgen_atmo generates different output for same flux input.

Gevgen and gevgen_atmo.pdf

TL;DR: this is an Executive Summary including the discussions and solution: 20230316_GENIEAGKYLowWDISBugFix_QYan4.pdf