Charm hadronization: GENIE assumes that hadr. system momentum is parallel to LAB z axis, which causes charm events to be rejected when it's not the case
marialiubarska opened this issue · 2 comments
Hi! I initially found this issue in GENIE 2.12.8, but I think it is still there in the master:
in 2.12.8:
Generator/src/Fragmentation/CharmHadronization.cxx
Lines 168 to 188 in d154450
in master:
Generator/src/Physics/Hadronization/AGCharm2019.cxx
Lines 251 to 272 in dc1ef00
The issue causes charm events to be often rejected as unphysical when incoming neutrino direction is not oriented along LAB's z axis and leads to the incorrect charm/not charm ratio in the simulation.
Dear Maria.
I am not sure whether your comment is correct. I seem to think that AGCharm2019::Hadronize() initially generates the hadronic system in a special frame: the hadronic CM frame, where z is in the direction of the momentum transfer. I do not think that the special status of the z axis remains in the final event, and we perform all the boosts and rotations to bring the particle 4-momenta back in the Lab frame. Notice the `particle -> P4() -> RotateUz( unitvq );' line (line 133) in AGCharm2019.cxx. Now, my memory of that code is a bit sketchy and it may well be that I missed something when I briefly looked at it again today. Please, if you can, pass some more information on what you think is wrong and we can look at it again. Do you have examples of events that are rejected as unphysical? How does that bias the charm fraction / what charm fractions do you obtain?
Kind regards
Costas
Dear Costas,
Thank you for the quick response!
I think you are right about AGCharm2019 - I only tested this in GENIE 2_12_8, which was used for atmospheric neutrino simulation after I noticed lack of charm events in Monte-Carlo. I then saw that the Hadronize() function is the same as in src/Fragmentation/CharmHadronization.cxx in 2_12_8 and also did not notice anything similar to this in the release notes, so I assumed the issue is still there in master. Now I see that the issue was probably fixed by this commit in 2019. Sorry!
However, the issue is still there in release 2_12_8 (I understood that it was fixed in R-3_00_06_NT), which causes almost all charm events for neutrinos with momentum not along LAB z axis to be rejected as unphysical with this kind of error:
DEBUG (genie-icetray): GENIE said: 1669332045 INFO HadronicVtx : [n] <HadronicSystemGenerator.cxx::Hadronic4pLAB (151)> :
v [LAB]: (E = 175.746, px = -128.395, py = 80.3133, pz = 89.1684) / M^2 = -7.27596e-12 / P = 175.746
N [LAB]: (E = 0.925253, px = -0.0325463, py = -0.137003, pz = 0.0713208) / M = 0.911689 / P = 0.157847
l [LAB]: (E = 39.121, px = -26.8493, py = 17.9183, pz = 22.1019) / M = 0.10566 / P = 39.1208
(I3LoggingAppender.cxx:68 in virtual void log4cpp::I3LoggingAppender::_append(const log4cpp::LoggingEvent&))
DEBUG (genie-icetray): GENIE said: 1669332045 INFO HadronicVtx : [n] <HadronicSystemGenerator.cxx::Hadronic4pLAB (164)> :
HadrSyst [LAB]: (E = 137.55, px = -101.578, py = 62.258, pz = 67.1378) / M = 14.7736 / P = 136.754
(I3LoggingAppender.cxx:68 in virtual void log4cpp::I3LoggingAppender::_append(const log4cpp::LoggingEvent&))
TRACE (genie-icetray): GENIE said: 1669332045 DEBUG Interaction : [n] <Kinematics.cxx::SetKV (347)> : Setting *Running* Hadronic invariant mass W to 14.7736
(I3LoggingAppender.cxx:76 in virtual void log4cpp::I3LoggingAppender::_append(const log4cpp::LoggingEvent&))
DEBUG (genie-icetray): GENIE said: 1669332045 NOTICE CharmHad : [n] <CharmHadronization.cxx::Hadronize (92)> : ** Running CHARM hadronizer
(I3LoggingAppender.cxx:68 in virtual void log4cpp::I3LoggingAppender::_append(const log4cpp::LoggingEvent&))
DEBUG (genie-icetray): GENIE said: 1669332045 NOTICE CharmHad : [n] <CharmHadronization.cxx::Hadronize (114)> : Ehad (LAB) = 137.55, W = 14.7736
(I3LoggingAppender.cxx:68 in virtual void log4cpp::I3LoggingAppender::_append(const log4cpp::LoggingEvent&))
DEBUG (genie-icetray): GENIE said: 1669332045 NOTICE CharmHad : [n] <CharmHadronization.cxx::Hadronize (141)> : Trying charm hadron = -421(m = 1.8645)
(I3LoggingAppender.cxx:68 in virtual void log4cpp::I3LoggingAppender::_append(const log4cpp::LoggingEvent&))
DEBUG (genie-icetray): GENIE said: 1669332045 INFO CharmHad : [n] <CharmHadronization.cxx::Hadronize (154)> : Trying charm hadron z = 0.711458, E = 97.861
(I3LoggingAppender.cxx:68 in virtual void log4cpp::I3LoggingAppender::_append(const log4cpp::LoggingEvent&))
DEBUG (genie-icetray): GENIE said: 1669332045 INFO CharmHad : [n] <CharmHadronization.cxx::Hadronize (163)> : Trying charm hadron pT^2 (tranv to pHad) = 0.0525861
(I3LoggingAppender.cxx:68 in virtual void log4cpp::I3LoggingAppender::_append(const log4cpp::LoggingEvent&))
TRACE (genie-icetray): GENIE said: 1669332045 DEBUG CharmHad : [n] <CharmHadronization.cxx::Hadronize (179)> : Charm hadron p4 (@LAB') = (E = 97.861, px = 0.0336402, py = 0.226836, pz = 97.8429) / M = 1.8645 / P = 97.8432
(I3LoggingAppender.cxx:76 in virtual void log4cpp::I3LoggingAppender::_append(const log4cpp::LoggingEvent&))
TRACE (genie-icetray): GENIE said: 1669332045 DEBUG CharmHad : [n] <CharmHadronization.cxx::Hadronize (184)> : Charm hadron p4 (@HCM') = (E = 465.771, px = 375.898, py = -230.142, pz = -150.583) / M = 1.8645 / P = 465.768
(I3LoggingAppender.cxx:76 in virtual void log4cpp::I3LoggingAppender::_append(const log4cpp::LoggingEvent&))
DEBUG (genie-icetray): GENIE said: 1669332045 INFO CharmHad : [n] <CharmHadronization.cxx::Hadronize (190)> : Invariant mass of remnant hadronic system= -116.364
(I3LoggingAppender.cxx:68 in virtual void log4cpp::I3LoggingAppender::_append(const log4cpp::LoggingEvent&))
DEBUG (genie-icetray): GENIE said: 1669332045 INFO CharmHad : [n] <CharmHadronization.cxx::Hadronize (193)> : Too small hadronic remnant mass!
...
INFO (genie-icetray): GENIE said: 1669313326 WARN CharmHad : [n] <CharmHadronization.cxx::Hadronize (302)> : Couldn't generate charm hadron for:
--------------------------------------------------------------------------------------------------------------
GENIE Interaction Summary
--------------------------------------------------------------------------------------------------------------
[-] [Init-State]
|--> probe : PDG-code = -14 (nu_mu_bar)
|--> nucl. target : Z = 8, A = 16, PDG-Code = 1000080160 (O16)
|--> hit nucleon : PDC-Code = 2112 (neutron)
|--> hit quark : PDC-Code = -3 (s_bar) [sea]
|--> probe 4P : (E = 251.092, Px = 43.0619, Py = 6.22431, Pz = -247.294)
|--> target 4P : (E = 14.8951, Px = 0, Py = 0, Pz = 0)
|--> nucleon 4P : (E = 0.918839, Px = 0.0328008, Py = -0.27124, Pz = 0.258508)
[-] [Process-Info]
|--> Interaction : Weak[CC]
|--> Scattering : DIS
[-] [Kinematics]
|--> *Running* Hadronic invariant mass W = 21.1162
|--> *Selected* Bjorken x = 0.0376547
|--> *Selected* Inelasticity y = 0.784309
|--> *Selected* Momentum transfer Q2 (>0) = 17.4195
|--> *Selected* Hadronic invariant mass W = 21.1162
[-] [Exclusive Process Info]
|--> charm prod. : true [inclusive] |--> 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 : [not set]
--------------------------------------------------------------------------------------------------------------
Unfortunately we are stuck with release 2_12_8 because the bulk of Monte-Carlo was produced few years ago and takes a long time to process. I'm looking for a way to produce an additional sample which can be used to correct existing MC. I don't know what is your policy on this, but is there a chance this can be patched in previous releases? Or would you suggest to go for a local fix?
Best regards,
Maria