rjones30
commented
3 years ago Mark shared with me the location of a subset of his genBH sample, which I found at this location: /work/halld/home/dalton/data/GDH/mcwrapper/hddm/dana_rest_genBH_sampweight_1e7041563[000-499].hddm I found 20k events in each hddm file x 500 files = 10M events. The following code was used to generate the same 1D plots as Mark shows above. My results are very different from Mark's, as shown in the 6 plots below that mirror his. Here are my primary observations.
- x is always 1, as it must be for elastic scattering.
- because x=1 we must have Q2=2 m_proton nu, so Q2 and nu are identical apart from the units
- the sign of nu and Q2 is always the same, and always positive, unlike what Mark saw with his code
- there is an anomaly with these special events piling up at Q2=1e-12 GeV^2, but this seems to be unrelated to what Mark showed.
- there is no cutoff in Q2, as Simon said, the Q2 extends well above 1 GeV^2 and extends up to the kinematic limit for this beam energy.
#!/usr/bin/python3
#
# md_study.py - utilities to study the genBH generator bugs reported by
# Mark Dalton in January, 2021.
#
# author: richard.t.jones at uconn.edu
# version: february 3, 2021
import ROOT
import glob
import numpy as np
import hddm_r
mElectron = 0.51099895e-3 # GeV/c^2
mProton = 0.93827208816 # GeV/c^2
hQ2 = ROOT.TH1D("hQ2", "", 100, -15, 2)
hnu = ROOT.TH1D("hnu", "", 100, -15, 2)
hx = ROOT.TH1D("hx", "", 100, -15, 2)
wQ2 = ROOT.TH1D("wQ2", "", 100, -15, 2)
wnu = ROOT.TH1D("wnu", "", 100, -15, 2)
wx = ROOT.TH1D("wx", "", 100, -15, 2)
log10 = np.log(10)
for infile in glob.glob("MD-1-2021/*"):
for rec in hddm_r.istream(infile):
rea = rec.getReactions()
prods = rec.getProducts()
if len(prods) != 3:
continue
nppho = np.array([rea[0].Ebeam, 0, 0, rea[0].Ebeam])
nptar = np.array([mProton, 0, 0, 0])
pele = prods[0].getMomentum()
npele = np.array([pele.E, pele.px, pele.py, pele.pz])
ppos = prods[1].getMomentum()
nppos = np.array([ppos.E, ppos.px, ppos.py, ppos.pz])
prec = prods[2].getMomentum()
nprec = np.array([prec.E, prec.px, prec.py, prec.pz])
Q2 = ((nptar[1] - nprec[1])**2 + (nptar[2] - nprec[2])**2 +
(nptar[3] - nprec[3])**2 - (nptar[0] - nprec[0])**2)
prec2 = nprec[1]**2 + nprec[2]**2 + nprec[3]**2
if prec2 / mProton < 1e-3:
nu = (nprec[1]**2 + nprec[2]**2 + nprec[3]**2) / (2 * mProton)
else:
nu = nprec[0] - nptar[0]
x = Q2 / (2 * mProton * nu)
hQ2.Fill(np.log(Q2)/log10)
hnu.Fill(np.log(nu)/log10)
hx.Fill(np.log(x)/log10)
wgt = rea[0].weight
wQ2.Fill(np.log(Q2)/log10, wgt)
wnu.Fill(np.log(nu)/log10, wgt)
wx.Fill(np.log(x)/log10, wgt)
From Mark Dalton, Jan 20, 2021, 10:47 AM
I am trying to look at the kinematics of the events from the genBH generator to determine, amongst other things, which form factor parametrization is most appropriate. Here I look at the values of Q2, ν and x and the correlations.
I started with 50M events that I had already simulated some months ago using genBH. Q2 and ν are both very small and I get different values if I use the target proton and recoil proton or if I use the beam photon and the 2 leptons. Both ν and Q2 have both positive and negative values. Might these issues be due to round off in the event storage?
Most of the events came with negative ν, so I flipped the sign of ν. I don't know why I needed to flip the sign perhaps a mistake on my part—but I did double check the algebra. My default algebra was done with the root class TLorentzVector which makes it easy but I cannot check for roundoff issues internally. Ι then eliminated any events where Q2 or ν were negative and plotted log of them in an effort to get an idea of our distribution. The first plot shows 2 correlations, Q2 vs ν and Q2 vs x, on top are the number of generated events and on the bottom are the events after weighting. The second plot is the distributions of Q2, ν and x unweighted (top) and weighted (bottom.)
As you can see, x has values greater than 1, which is not physical for a proton. The weighted histogram shows that the events are clustered "around” x=1, which would make sense if they are elastic and we have some roundoff or resolution issues somewhere. I would therefore say that most of the events come with Q2 ~ 4x10-7 and ν ~ 2x10-7. The abrupt cutoff at values a bit higher than this is something that I cannot explain either.