diff --git a/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_enuqe_qelike.root b/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_enuqe_qelike.root
new file mode 100644
index 0000000..07d1e73
Binary files /dev/null and b/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_enuqe_qelike.root differ
diff --git a/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_enuqe_qelike.txt b/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_enuqe_qelike.txt
new file mode 100644
index 0000000..344bb51
--- /dev/null
+++ b/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_enuqe_qelike.txt
@@ -0,0 +1,144 @@
+1,1,2,2,7.879204e-79
+1,2,2,3,4.045300e-79
+1,3,2,4,3.388994e-79
+1,4,2,5,1.562272e-79
+1,5,2,6,-8.980128e-80
+1,6,2,7,-1.964155e-79
+1,7,2,8,-1.926981e-79
+1,8,2,9,-5.912014e-80
+1,9,2,10,8.045190e-80
+1,10,2,11,6.959615e-80
+1,11,2,12,1.757468e-79
+1,12,2,13,-1.849971e-79
+2,1,3,2,4.045300e-79
+2,2,3,3,3.115976e-79
+2,3,3,4,2.789322e-79
+2,4,3,5,1.951562e-79
+2,5,3,6,7.641406e-80
+2,6,3,7,1.974635e-80
+2,7,3,8,2.377597e-80
+2,8,3,9,1.039844e-79
+2,9,3,10,2.049484e-79
+2,10,3,11,1.974055e-79
+2,11,3,12,2.421120e-79
+2,12,3,13,1.057254e-79
+3,1,4,2,3.388994e-79
+3,2,4,3,2.789322e-79
+3,3,4,4,2.855821e-79
+3,4,4,5,2.140941e-79
+3,5,4,6,1.196706e-79
+3,6,4,7,7.289655e-80
+3,7,4,8,7.453811e-80
+3,8,4,9,1.392040e-79
+3,9,4,10,2.285917e-79
+3,10,4,11,2.232450e-79
+3,11,4,12,2.574354e-79
+3,12,4,13,1.630178e-79
+4,1,5,2,1.562272e-79
+4,2,5,3,1.951562e-79
+4,3,5,4,2.140941e-79
+4,4,5,5,2.174830e-79
+4,5,5,6,1.959540e-79
+4,6,5,7,1.912191e-79
+4,7,5,8,2.013867e-79
+4,8,5,9,2.120246e-79
+4,9,5,10,2.412341e-79
+4,10,5,11,2.319603e-79
+4,11,5,12,2.387827e-79
+4,12,5,13,2.413977e-79
+5,1,6,2,-8.980128e-80
+5,2,6,3,7.641406e-80
+5,3,6,4,1.196706e-79
+5,4,6,5,1.959540e-79
+5,5,6,6,3.025884e-79
+5,6,6,7,3.670362e-79
+5,7,6,8,3.935535e-79
+5,8,6,9,3.226150e-79
+5,9,6,10,2.591500e-79
+5,10,6,11,2.425632e-79
+5,11,6,12,2.088502e-79
+5,12,6,13,3.466435e-79
+6,1,7,2,-1.964155e-79
+6,2,7,3,1.974635e-80
+6,3,7,4,7.289655e-80
+6,4,7,5,1.912191e-79
+6,5,7,6,3.670362e-79
+6,6,7,7,5.109811e-79
+6,7,7,8,5.627839e-79
+6,8,7,9,4.211592e-79
+6,9,7,10,2.840694e-79
+6,10,7,11,2.571972e-79
+6,11,7,12,2.030374e-79
+6,12,7,13,4.078384e-79
+7,1,8,2,-1.926981e-79
+7,2,8,3,2.377597e-80
+7,3,8,4,7.453811e-80
+7,4,8,5,2.013867e-79
+7,5,8,6,3.935535e-79
+7,6,8,7,5.627839e-79
+7,7,8,8,6.641857e-79
+7,8,8,9,5.019325e-79
+7,9,8,10,3.360928e-79
+7,10,8,11,3.033802e-79
+7,11,8,12,2.463045e-79
+7,12,8,13,4.633835e-79
+8,1,9,2,-5.912014e-80
+8,2,9,3,1.039844e-79
+8,3,9,4,1.392040e-79
+8,4,9,5,2.120246e-79
+8,5,9,6,3.226150e-79
+8,6,9,7,4.211592e-79
+8,7,9,8,5.019325e-79
+8,8,9,9,4.776774e-79
+8,9,9,10,4.025047e-79
+8,10,9,11,3.748314e-79
+8,11,9,12,3.420851e-79
+8,12,9,13,4.995705e-79
+9,1,10,2,8.045190e-80
+9,2,10,3,2.049484e-79
+9,3,10,4,2.285917e-79
+9,4,10,5,2.412341e-79
+9,5,10,6,2.591500e-79
+9,6,10,7,2.840694e-79
+9,7,10,8,3.360928e-79
+9,8,10,9,4.025047e-79
+9,9,10,10,5.369403e-79
+9,10,10,11,5.104162e-79
+9,11,10,12,4.963167e-79
+9,12,10,13,6.085889e-79
+10,1,11,2,6.959615e-80
+10,2,11,3,1.974055e-79
+10,3,11,4,2.232450e-79
+10,4,11,5,2.319603e-79
+10,5,11,6,2.425632e-79
+10,6,11,7,2.571972e-79
+10,7,11,8,3.033802e-79
+10,8,11,9,3.748314e-79
+10,9,11,10,5.104162e-79
+10,10,11,11,5.719368e-79
+10,11,11,12,5.198401e-79
+10,12,11,13,6.381431e-79
+11,1,12,2,1.757468e-79
+11,2,12,3,2.421120e-79
+11,3,12,4,2.574354e-79
+11,4,12,5,2.387827e-79
+11,5,12,6,2.088502e-79
+11,6,12,7,2.030374e-79
+11,7,12,8,2.463045e-79
+11,8,12,9,3.420851e-79
+11,9,12,10,4.963167e-79
+11,10,12,11,5.198401e-79
+11,11,12,12,5.812268e-79
+11,12,12,13,6.182086e-79
+12,1,13,2,-1.849971e-79
+12,2,13,3,1.057254e-79
+12,3,13,4,1.630178e-79
+12,4,13,5,2.413977e-79
+12,5,13,6,3.466435e-79
+12,6,13,7,4.078384e-79
+12,7,13,8,4.633835e-79
+12,8,13,9,4.995705e-79
+12,9,13,10,6.085889e-79
+12,10,13,11,6.381431e-79
+12,11,13,12,6.182086e-79
+12,12,13,13,1.055083e-78
diff --git a/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_ptmu_qelike.root b/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_ptmu_qelike.root
new file mode 100644
index 0000000..f4feb6d
Binary files /dev/null and b/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_ptmu_qelike.root differ
diff --git a/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_ptmu_qelike.txt b/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_ptmu_qelike.txt
new file mode 100644
index 0000000..13fd4fb
--- /dev/null
+++ b/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_ptmu_qelike.txt
@@ -0,0 +1,169 @@
+1,1,2,2,4.090292e-81
+1,2,2,3,2.562269e-81
+1,3,2,4,6.989749e-81
+1,4,2,5,1.183908e-80
+1,5,2,6,1.550956e-80
+1,6,2,7,1.685813e-80
+1,7,2,8,1.670301e-80
+1,8,2,9,1.213871e-80
+1,9,2,10,5.493265e-81
+1,10,2,11,1.819849e-81
+1,11,2,12,2.646811e-82
+1,12,2,13,-4.063299e-83
+1,13,2,14,-8.355687e-84
+2,1,3,2,2.562269e-81
+2,2,3,3,2.273044e-80
+2,3,3,4,3.004625e-80
+2,4,3,5,5.078005e-80
+2,5,3,6,6.619245e-80
+2,6,3,7,7.140398e-80
+2,7,3,8,7.007077e-80
+2,8,3,9,5.070531e-80
+2,9,3,10,2.215350e-80
+2,10,3,11,7.260039e-81
+2,11,3,12,8.589517e-82
+2,12,3,13,-2.546520e-82
+2,13,3,14,-4.703057e-83
+3,1,4,2,6.989749e-81
+3,2,4,3,3.004625e-80
+3,3,4,4,1.024760e-79
+3,4,4,5,1.418386e-79
+3,5,4,6,1.787637e-79
+3,6,4,7,1.930432e-79
+3,7,4,8,1.897062e-79
+3,8,4,9,1.369680e-79
+3,9,4,10,5.949686e-80
+3,10,4,11,1.778499e-80
+3,11,4,12,1.686449e-81
+3,12,4,13,-9.768174e-82
+3,13,4,14,-1.352166e-82
+4,1,5,2,1.183908e-80
+4,2,5,3,5.078005e-80
+4,3,5,4,1.418386e-79
+4,4,5,5,2.700500e-79
+4,5,5,6,3.013533e-79
+4,6,5,7,3.272854e-79
+4,7,5,8,3.239328e-79
+4,8,5,9,2.363072e-79
+4,9,5,10,1.057790e-79
+4,10,5,11,3.451906e-80
+4,11,5,12,4.628485e-81
+4,12,5,13,-1.127942e-81
+4,13,5,14,-1.906346e-82
+5,1,6,2,1.550956e-80
+5,2,6,3,6.619245e-80
+5,3,6,4,1.787637e-79
+5,4,6,5,3.013533e-79
+5,5,6,6,4.451882e-79
+5,6,6,7,4.400985e-79
+5,7,6,8,4.357062e-79
+5,8,6,9,3.218683e-79
+5,9,6,10,1.475811e-79
+5,10,6,11,5.300167e-80
+5,11,6,12,9.335704e-81
+5,12,6,13,-3.024653e-82
+5,13,6,14,-1.327143e-82
+6,1,7,2,1.685813e-80
+6,2,7,3,7.140398e-80
+6,3,7,4,1.930432e-79
+6,4,7,5,3.272854e-79
+6,5,7,6,4.400985e-79
+6,6,7,7,5.290281e-79
+6,7,7,8,4.977330e-79
+6,8,7,9,3.710658e-79
+6,9,7,10,1.857477e-79
+6,10,7,11,7.816982e-80
+6,11,7,12,1.864065e-80
+6,12,7,13,2.276056e-81
+6,13,7,14,8.247696e-83
+7,1,8,2,1.670301e-80
+7,2,8,3,7.007077e-80
+7,3,8,4,1.897062e-79
+7,4,8,5,3.239328e-79
+7,5,8,6,4.357062e-79
+7,6,8,7,4.977330e-79
+7,7,8,8,5.328125e-79
+7,8,8,9,3.885338e-79
+7,9,8,10,2.056601e-79
+7,10,8,11,9.174284e-80
+7,11,8,12,2.471262e-80
+7,12,8,13,4.092175e-81
+7,13,8,14,2.268682e-82
+8,1,9,2,1.213871e-80
+8,2,9,3,5.070531e-80
+8,3,9,4,1.369680e-79
+8,4,9,5,2.363072e-79
+8,5,9,6,3.218683e-79
+8,6,9,7,3.710658e-79
+8,7,9,8,3.885338e-79
+8,8,9,9,3.378429e-79
+8,9,9,10,1.843751e-79
+8,10,9,11,8.456050e-80
+8,11,9,12,2.387651e-80
+8,12,9,13,4.404858e-81
+8,13,9,14,3.023314e-82
+9,1,10,2,5.493265e-81
+9,2,10,3,2.215350e-80
+9,3,10,4,5.949686e-80
+9,4,10,5,1.057790e-79
+9,5,10,6,1.475811e-79
+9,6,10,7,1.857477e-79
+9,7,10,8,2.056601e-79
+9,8,10,9,1.843751e-79
+9,9,10,10,1.397019e-79
+9,10,10,11,7.572276e-80
+9,11,10,12,2.597400e-80
+9,12,10,13,6.617491e-81
+9,13,10,14,5.208417e-82
+10,1,11,2,1.819849e-81
+10,2,11,3,7.260039e-81
+10,3,11,4,1.778499e-80
+10,4,11,5,3.451906e-80
+10,5,11,6,5.300167e-80
+10,6,11,7,7.816982e-80
+10,7,11,8,9.174284e-80
+10,8,11,9,8.456050e-80
+10,9,11,10,7.572276e-80
+10,10,11,11,5.714153e-80
+10,11,11,12,2.084385e-80
+10,12,11,13,5.833001e-81
+10,13,11,14,4.831035e-82
+11,1,12,2,2.646811e-82
+11,2,12,3,8.589517e-82
+11,3,12,4,1.686449e-81
+11,4,12,5,4.628485e-81
+11,5,12,6,9.335704e-81
+11,6,12,7,1.864065e-80
+11,7,12,8,2.471262e-80
+11,8,12,9,2.387651e-80
+11,9,12,10,2.597400e-80
+11,10,12,11,2.084385e-80
+11,11,12,12,9.553493e-81
+11,12,12,13,2.776667e-81
+11,13,12,14,2.261623e-82
+12,1,13,2,-4.063299e-83
+12,2,13,3,-2.546520e-82
+12,3,13,4,-9.768174e-82
+12,4,13,5,-1.127942e-81
+12,5,13,6,-3.024653e-82
+12,6,13,7,2.276056e-81
+12,7,13,8,4.092175e-81
+12,8,13,9,4.404858e-81
+12,9,13,10,6.617491e-81
+12,10,13,11,5.833001e-81
+12,11,13,12,2.776667e-81
+12,12,13,13,1.041425e-81
+12,13,13,14,8.252997e-83
+13,1,14,2,-8.355687e-84
+13,2,14,3,-4.703057e-83
+13,3,14,4,-1.352166e-82
+13,4,14,5,-1.906346e-82
+13,5,14,6,-1.327143e-82
+13,6,14,7,8.247696e-83
+13,7,14,8,2.268682e-82
+13,8,14,9,3.023314e-82
+13,9,14,10,5.208417e-82
+13,10,14,11,4.831035e-82
+13,11,14,12,2.261623e-82
+13,12,14,13,8.252997e-83
+13,13,14,14,1.055466e-83
diff --git a/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_pzmu_qelike.root b/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_pzmu_qelike.root
new file mode 100644
index 0000000..a5790c2
Binary files /dev/null and b/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_pzmu_qelike.root differ
diff --git a/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_pzmu_qelike.txt b/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_pzmu_qelike.txt
new file mode 100644
index 0000000..8d7ce8d
--- /dev/null
+++ b/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_pzmu_qelike.txt
@@ -0,0 +1,144 @@
+1,1,2,2,1.828348e-80
+1,2,2,3,1.738194e-80
+1,3,2,4,1.453120e-80
+1,4,2,5,3.968928e-81
+1,5,2,6,-5.841342e-82
+1,6,2,7,-6.917475e-82
+1,7,2,8,1.497378e-82
+1,8,2,9,6.892747e-82
+1,9,2,10,6.066394e-82
+1,10,2,11,4.239844e-82
+1,11,2,12,3.436567e-82
+1,12,2,13,4.418961e-84
+2,1,3,2,1.738194e-80
+2,2,3,3,2.115946e-80
+2,3,3,4,1.944783e-80
+2,4,3,5,1.050013e-80
+2,5,3,6,4.240760e-81
+2,6,3,7,2.058431e-81
+2,7,3,8,1.659607e-81
+2,8,3,9,1.575465e-81
+2,9,3,10,1.115327e-81
+2,10,3,11,7.349264e-82
+2,11,3,12,4.990915e-82
+2,12,3,13,1.233453e-82
+3,1,4,2,1.453120e-80
+3,2,4,3,1.944783e-80
+3,3,4,4,2.206542e-80
+3,4,4,5,1.534332e-80
+3,5,4,6,8.395976e-81
+3,6,4,7,4.508937e-81
+3,7,4,8,2.905860e-81
+3,8,4,9,2.178415e-81
+3,9,4,10,1.409052e-81
+3,10,4,11,9.157068e-82
+3,11,4,12,5.686401e-82
+3,12,4,13,2.140252e-82
+4,1,5,2,3.968928e-81
+4,2,5,3,1.050013e-80
+4,3,5,4,1.534332e-80
+4,4,5,5,1.750631e-80
+4,5,5,6,1.190941e-80
+4,6,5,7,6.688076e-81
+4,7,5,8,3.799810e-81
+4,8,5,9,2.344579e-81
+4,9,5,10,1.369806e-81
+4,10,5,11,8.637200e-82
+4,11,5,12,4.604361e-82
+4,12,5,13,2.870174e-82
+5,1,6,2,-5.841342e-82
+5,2,6,3,4.240760e-81
+5,3,6,4,8.395976e-81
+5,4,6,5,1.190941e-80
+5,5,6,6,9.281508e-81
+5,6,6,7,5.211823e-81
+5,7,6,8,2.827167e-81
+5,8,6,9,1.626820e-81
+5,9,6,10,9.106072e-82
+5,10,6,11,5.670440e-82
+5,11,6,12,2.768028e-82
+5,12,6,13,2.149749e-82
+6,1,7,2,-6.917475e-82
+6,2,7,3,2.058431e-81
+6,3,7,4,4.508937e-81
+6,4,7,5,6.688076e-81
+6,5,7,6,5.211823e-81
+6,6,7,7,3.160368e-81
+6,7,7,8,1.703479e-81
+6,8,7,9,9.524691e-82
+6,9,7,10,5.367817e-82
+6,10,7,11,3.340519e-82
+6,11,7,12,1.606655e-82
+6,12,7,13,1.274759e-82
+7,1,8,2,1.497378e-82
+7,2,8,3,1.659607e-81
+7,3,8,4,2.905860e-81
+7,4,8,5,3.799810e-81
+7,5,8,6,2.827167e-81
+7,6,8,7,1.703479e-81
+7,7,8,8,1.039154e-81
+7,8,8,9,5.820767e-82
+7,9,8,10,3.324698e-82
+7,10,8,11,2.093137e-82
+7,11,8,12,1.056706e-82
+7,12,8,13,7.373166e-83
+8,1,9,2,6.892747e-82
+8,2,9,3,1.575465e-81
+8,3,9,4,2.178415e-81
+8,4,9,5,2.344579e-81
+8,5,9,6,1.626820e-81
+8,6,9,7,9.524691e-82
+8,7,9,8,5.820767e-82
+8,8,9,9,4.039794e-82
+8,9,9,10,2.269601e-82
+8,10,9,11,1.427983e-82
+8,11,9,12,7.687687e-83
+8,12,9,13,4.481408e-83
+9,1,10,2,6.066394e-82
+9,2,10,3,1.115327e-81
+9,3,10,4,1.409052e-81
+9,4,10,5,1.369806e-81
+9,5,10,6,9.106072e-82
+9,6,10,7,5.367817e-82
+9,7,10,8,3.324698e-82
+9,8,10,9,2.269601e-82
+9,9,10,10,1.572030e-82
+9,10,10,11,9.337739e-83
+9,11,10,12,5.086832e-83
+9,12,10,13,2.729637e-83
+10,1,11,2,4.239844e-82
+10,2,11,3,7.349264e-82
+10,3,11,4,9.157068e-82
+10,4,11,5,8.637200e-82
+10,5,11,6,5.670440e-82
+10,6,11,7,3.340519e-82
+10,7,11,8,2.093137e-82
+10,8,11,9,1.427983e-82
+10,9,11,10,9.337739e-83
+10,10,11,11,6.812706e-83
+10,11,11,12,3.309201e-83
+10,12,11,13,1.709428e-83
+11,1,12,2,3.436567e-82
+11,2,12,3,4.990915e-82
+11,3,12,4,5.686401e-82
+11,4,12,5,4.604361e-82
+11,5,12,6,2.768028e-82
+11,6,12,7,1.606655e-82
+11,7,12,8,1.056706e-82
+11,8,12,9,7.687687e-83
+11,9,12,10,5.086832e-83
+11,10,12,11,3.309201e-83
+11,11,12,12,2.100304e-83
+11,12,12,13,8.605735e-84
+12,1,13,2,4.418961e-84
+12,2,13,3,1.233453e-82
+12,3,13,4,2.140252e-82
+12,4,13,5,2.870174e-82
+12,5,13,6,2.149749e-82
+12,6,13,7,1.274759e-82
+12,7,13,8,7.373166e-83
+12,8,13,9,4.481408e-83
+12,9,13,10,2.729637e-83
+12,10,13,11,1.709428e-83
+12,11,13,12,8.605735e-84
+12,12,13,13,6.874882e-84
diff --git a/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_q2qe_qelike.root b/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_q2qe_qelike.root
new file mode 100644
index 0000000..8b62cc4
Binary files /dev/null and b/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_q2qe_qelike.root differ
diff --git a/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_q2qe_qelike.txt b/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_q2qe_qelike.txt
new file mode 100644
index 0000000..6dd7b66
--- /dev/null
+++ b/data/MINERvA/CC0pi_1D/FixedBinWidthPub/cov_q2qe_qelike.txt
@@ -0,0 +1,256 @@
+1,1,2,2,3.364235e-79
+1,2,2,3,1.578055e-79
+1,3,2,4,2.189488e-79
+1,4,2,5,2.856008e-79
+1,5,2,6,3.109782e-79
+1,6,2,7,3.337411e-79
+1,7,2,8,3.253582e-79
+1,8,2,9,3.086935e-79
+1,9,2,10,2.372775e-79
+1,10,2,11,1.623871e-79
+1,11,2,12,9.632743e-80
+1,12,2,13,4.834720e-80
+1,13,2,14,2.247655e-80
+1,14,2,15,1.041239e-80
+1,15,2,16,2.112753e-81
+1,16,2,17,-6.518251e-83
+2,1,3,2,1.578055e-79
+2,2,3,3,3.179707e-79
+2,3,3,4,3.149551e-79
+2,4,3,5,3.268118e-79
+2,5,3,6,3.629008e-79
+2,6,3,7,3.965343e-79
+2,7,3,8,3.806013e-79
+2,8,3,9,3.525447e-79
+2,9,3,10,2.630355e-79
+2,10,3,11,1.763856e-79
+2,11,3,12,9.937935e-80
+2,12,3,13,4.797530e-80
+2,13,3,14,1.980317e-80
+2,14,3,15,8.130958e-81
+2,15,3,16,7.276350e-82
+2,16,3,17,-2.832949e-82
+3,1,4,2,2.189488e-79
+3,2,4,3,3.149551e-79
+3,3,4,4,5.374580e-79
+3,4,4,5,5.251322e-79
+3,5,4,6,5.415394e-79
+3,6,4,7,5.915613e-79
+3,7,4,8,5.683775e-79
+3,8,4,9,5.270431e-79
+3,9,4,10,3.895419e-79
+3,10,4,11,2.593050e-79
+3,11,4,12,1.457446e-79
+3,12,4,13,6.859378e-80
+3,13,4,14,2.666378e-80
+3,14,4,15,1.015332e-80
+3,15,4,16,3.524381e-82
+3,16,4,17,-5.357672e-82
+4,1,5,2,2.856008e-79
+4,2,5,3,3.268118e-79
+4,3,5,4,5.251322e-79
+4,4,5,5,6.987978e-79
+4,5,5,6,6.884911e-79
+4,6,5,7,7.014411e-79
+4,7,5,8,6.815112e-79
+4,8,5,9,6.356982e-79
+4,9,5,10,4.756715e-79
+4,10,5,11,3.199762e-79
+4,11,5,12,1.821589e-79
+4,12,5,13,8.673351e-80
+4,13,5,14,3.583119e-80
+4,14,5,15,1.449002e-80
+4,15,5,16,1.405239e-81
+4,16,5,17,-5.071381e-82
+5,1,6,2,3.109782e-79
+5,2,6,3,3.629008e-79
+5,3,6,4,5.415394e-79
+5,4,6,5,6.884911e-79
+5,5,6,6,8.141007e-79
+5,6,6,7,7.738982e-79
+5,7,6,8,7.341207e-79
+5,8,6,9,6.921160e-79
+5,9,6,10,5.145867e-79
+5,10,6,11,3.427129e-79
+5,11,6,12,1.982453e-79
+5,12,6,13,9.364193e-80
+5,13,6,14,3.630997e-80
+5,14,6,15,1.387735e-80
+5,15,6,16,8.677617e-82
+5,16,6,17,-6.424076e-82
+6,1,7,2,3.337411e-79
+6,2,7,3,3.965343e-79
+6,3,7,4,5.915613e-79
+6,4,7,5,7.014411e-79
+6,5,7,6,7.738982e-79
+6,6,7,7,8.554283e-79
+6,7,7,8,7.879186e-79
+6,8,7,9,7.433642e-79
+6,9,7,10,5.656739e-79
+6,10,7,11,3.837563e-79
+6,11,7,12,2.236046e-79
+6,12,7,13,1.084635e-79
+6,13,7,14,4.679260e-80
+6,14,7,15,2.016688e-80
+6,15,7,16,2.959829e-81
+6,16,7,17,-4.369278e-82
+7,1,8,2,3.253582e-79
+7,2,8,3,3.806013e-79
+7,3,8,4,5.683775e-79
+7,4,8,5,6.815112e-79
+7,5,8,6,7.341207e-79
+7,6,8,7,7.879186e-79
+7,7,8,8,8.088892e-79
+7,8,8,9,7.367133e-79
+7,9,8,10,5.633745e-79
+7,10,8,11,3.920199e-79
+7,11,8,12,2.340422e-79
+7,12,8,13,1.177491e-79
+7,13,8,14,5.536190e-80
+7,14,8,15,2.594136e-80
+7,15,8,16,5.464148e-81
+7,16,8,17,-8.867431e-83
+8,1,9,2,3.086935e-79
+8,2,9,3,3.525447e-79
+8,3,9,4,5.270431e-79
+8,4,9,5,6.356982e-79
+8,5,9,6,6.921160e-79
+8,6,9,7,7.433642e-79
+8,7,9,8,7.367133e-79
+8,8,9,9,7.379343e-79
+8,9,9,10,5.547608e-79
+8,10,9,11,3.855583e-79
+8,11,9,12,2.370350e-79
+8,12,9,13,1.208208e-79
+8,13,9,14,5.987985e-80
+8,14,9,15,2.927504e-80
+8,15,9,16,7.101360e-81
+8,16,9,17,1.339414e-82
+9,1,10,2,2.372775e-79
+9,2,10,3,2.630355e-79
+9,3,10,4,3.895419e-79
+9,4,10,5,4.756715e-79
+9,5,10,6,5.145867e-79
+9,6,10,7,5.656739e-79
+9,7,10,8,5.633745e-79
+9,8,10,9,5.547608e-79
+9,9,10,10,4.682480e-79
+9,10,10,11,3.277024e-79
+9,11,10,12,2.069763e-79
+9,12,10,13,1.123189e-79
+9,13,10,14,6.156802e-80
+9,14,10,15,3.261566e-80
+9,15,10,16,9.628625e-81
+9,16,10,17,6.354506e-82
+10,1,11,2,1.623871e-79
+10,2,11,3,1.763856e-79
+10,3,11,4,2.593050e-79
+10,4,11,5,3.199762e-79
+10,5,11,6,3.427129e-79
+10,6,11,7,3.837563e-79
+10,7,11,8,3.920199e-79
+10,8,11,9,3.855583e-79
+10,9,11,10,3.277024e-79
+10,10,11,11,2.621492e-79
+10,11,11,12,1.667422e-79
+10,12,11,13,9.187279e-80
+10,13,11,14,5.292266e-80
+10,14,11,15,2.848306e-80
+10,15,11,16,8.861856e-81
+10,16,11,17,6.954773e-82
+11,1,12,2,9.632743e-80
+11,2,12,3,9.937935e-80
+11,3,12,4,1.457446e-79
+11,4,12,5,1.821589e-79
+11,5,12,6,1.982453e-79
+11,6,12,7,2.236046e-79
+11,7,12,8,2.340422e-79
+11,8,12,9,2.370350e-79
+11,9,12,10,2.069763e-79
+11,10,12,11,1.667422e-79
+11,11,12,12,1.193502e-79
+11,12,12,13,6.598504e-80
+11,13,12,14,3.853129e-80
+11,14,12,15,2.144629e-80
+11,15,12,16,7.141230e-81
+11,16,12,17,6.409432e-82
+12,1,13,2,4.834720e-80
+12,2,13,3,4.797530e-80
+12,3,13,4,6.859378e-80
+12,4,13,5,8.673351e-80
+12,5,13,6,9.364193e-80
+12,6,13,7,1.084635e-79
+12,7,13,8,1.177491e-79
+12,8,13,9,1.208208e-79
+12,9,13,10,1.123189e-79
+12,10,13,11,9.187279e-80
+12,11,13,12,6.598504e-80
+12,12,13,13,4.214454e-80
+12,13,13,14,2.596900e-80
+12,14,13,15,1.461251e-80
+12,15,13,16,5.186746e-81
+12,16,13,17,5.304958e-82
+13,1,14,2,2.247655e-80
+13,2,14,3,1.980317e-80
+13,3,14,4,2.666378e-80
+13,4,14,5,3.583119e-80
+13,5,14,6,3.630997e-80
+13,6,14,7,4.679260e-80
+13,7,14,8,5.536190e-80
+13,8,14,9,5.987985e-80
+13,9,14,10,6.156802e-80
+13,10,14,11,5.292266e-80
+13,11,14,12,3.853129e-80
+13,12,14,13,2.596900e-80
+13,13,14,14,2.003831e-80
+13,14,14,15,1.170769e-80
+13,15,14,16,4.245816e-81
+13,16,14,17,4.787901e-82
+14,1,15,2,1.041239e-80
+14,2,15,3,8.130958e-81
+14,3,15,4,1.015332e-80
+14,4,15,5,1.449002e-80
+14,5,15,6,1.387735e-80
+14,6,15,7,2.016688e-80
+14,7,15,8,2.594136e-80
+14,8,15,9,2.927504e-80
+14,9,15,10,3.261566e-80
+14,10,15,11,2.848306e-80
+14,11,15,12,2.144629e-80
+14,12,15,13,1.461251e-80
+14,13,15,14,1.170769e-80
+14,14,15,15,7.702367e-81
+14,15,15,16,2.862324e-81
+14,16,15,17,3.233580e-82
+15,1,16,2,2.112753e-81
+15,2,16,3,7.276350e-82
+15,3,16,4,3.524381e-82
+15,4,16,5,1.405239e-81
+15,5,16,6,8.677617e-82
+15,6,16,7,2.959829e-81
+15,7,16,8,5.464148e-81
+15,8,16,9,7.101360e-81
+15,9,16,10,9.628625e-81
+15,10,16,11,8.861856e-81
+15,11,16,12,7.141230e-81
+15,12,16,13,5.186746e-81
+15,13,16,14,4.245816e-81
+15,14,16,15,2.862324e-81
+15,15,16,16,1.247998e-81
+15,16,16,17,1.468619e-82
+16,1,17,2,-6.518251e-83
+16,2,17,3,-2.832949e-82
+16,3,17,4,-5.357672e-82
+16,4,17,5,-5.071381e-82
+16,5,17,6,-6.424076e-82
+16,6,17,7,-4.369278e-82
+16,7,17,8,-8.867431e-83
+16,8,17,9,1.339414e-82
+16,9,17,10,6.354506e-82
+16,10,17,11,6.954773e-82
+16,11,17,12,6.409432e-82
+16,12,17,13,5.304958e-82
+16,13,17,14,4.787901e-82
+16,14,17,15,3.233580e-82
+16,15,17,16,1.468619e-82
+16,16,17,17,2.297338e-83
diff --git a/data/MINERvA/CC0pi_2D/cov_ptpl_2D_qelike.root b/data/MINERvA/CC0pi_2D/cov_ptpl_2D_qelike.root
new file mode 100644
index 0000000..d8788ef
Binary files /dev/null and b/data/MINERvA/CC0pi_2D/cov_ptpl_2D_qelike.root differ
diff --git a/parameters/config.xml b/parameters/config.xml
index 49dc841..b4a01c0 100644
--- a/parameters/config.xml
+++ b/parameters/config.xml
@@ -1,212 +1,215 @@
+
+
+
diff --git a/src/InputHandler/GENIEInputHandler.cxx b/src/InputHandler/GENIEInputHandler.cxx
index a81cd28..a5eb414 100644
--- a/src/InputHandler/GENIEInputHandler.cxx
+++ b/src/InputHandler/GENIEInputHandler.cxx
@@ -1,504 +1,503 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifdef __GENIE_ENABLED__
#include "GENIEInputHandler.h"
#include "InputUtils.h"
GENIEGeneratorInfo::~GENIEGeneratorInfo() { DeallocateParticleStack(); }
void GENIEGeneratorInfo::AddBranchesToTree(TTree* tn) {
tn->Branch("GenieParticlePDGs", &fGenieParticlePDGs, "GenieParticlePDGs/I");
}
void GENIEGeneratorInfo::SetBranchesFromTree(TTree* tn) {
tn->SetBranchAddress("GenieParticlePDGs", &fGenieParticlePDGs);
}
void GENIEGeneratorInfo::AllocateParticleStack(int stacksize) {
fGenieParticlePDGs = new int[stacksize];
}
void GENIEGeneratorInfo::DeallocateParticleStack() {
delete fGenieParticlePDGs;
}
void GENIEGeneratorInfo::FillGeneratorInfo(NtpMCEventRecord* ntpl) {
Reset();
// Check for GENIE Event
if (!ntpl) return;
if (!ntpl->event) return;
// Cast Event Record
GHepRecord* ghep = static_cast(ntpl->event);
if (!ghep) return;
// Fill Particle Stack
GHepParticle* p = 0;
TObjArrayIter iter(ghep);
// Loop over all particles
int i = 0;
while ((p = (dynamic_cast((iter).Next())))) {
if (!p) continue;
// Get PDG
fGenieParticlePDGs[i] = p->Pdg();
i++;
}
}
void GENIEGeneratorInfo::Reset() {
for (int i = 0; i < kMaxParticles; i++) {
fGenieParticlePDGs[i] = 0;
}
}
GENIEInputHandler::GENIEInputHandler(std::string const& handle,
std::string const& rawinputs) {
LOG(SAM) << "Creating GENIEInputHandler : " << handle << std::endl;
// Run a joint input handling
fName = handle;
// Setup the TChain
fGENIETree = new TChain("gtree");
fSaveExtra = FitPar::Config().GetParB("SaveExtraGenie");
fCacheSize = FitPar::Config().GetParI("CacheSize");
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
-
// Are we running with NOvA weights
- bool nova_wgt = true;
+ fNOvAWeights = FitPar::Config().GetParB("NOvA_Weights");
MAQEw = 1.0;
NonResw = 1.0;
RPAQEw = 1.0;
RPARESw = 1.0;
MECw = 1.0;
DISw = 1.0;
NOVAw = 1.0;
// Loop over all inputs and grab flux, eventhist, and nevents
std::vector inputs = InputUtils::ParseInputFileList(rawinputs);
for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) {
// Open File for histogram access
TFile* inp_file = new TFile(
InputUtils::ExpandInputDirectories(inputs[inp_it]).c_str(), "READ");
if (!inp_file or inp_file->IsZombie()) {
THROW("GENIE File IsZombie() at : '"
<< inputs[inp_it] << "'" << std::endl
<< "Check that your file paths are correct and the file exists!"
<< std::endl
<< "$ ls -lh " << inputs[inp_it]);
}
// Get Flux/Event hist
TH1D* fluxhist = (TH1D*)inp_file->Get("nuisance_flux");
TH1D* eventhist = (TH1D*)inp_file->Get("nuisance_events");
if (!fluxhist or !eventhist) {
ERROR(FTL, "Input File Contents: " << inputs[inp_it]);
inp_file->ls();
THROW("GENIE FILE doesn't contain flux/xsec info."
<< std::endl
<< "Try running the app PrepareGENIE first on :" << inputs[inp_it]
<< std::endl
<< "$ PrepareGENIE -h");
}
// Get N Events
TTree* genietree = (TTree*)inp_file->Get("gtree");
if (!genietree) {
ERROR(FTL, "gtree not located in GENIE file: " << inputs[inp_it]);
THROW("Check your inputs, they may need to be completely regenerated!");
throw;
}
int nevents = genietree->GetEntries();
if (nevents <= 0) {
THROW("Trying to a TTree with "
<< nevents << " to TChain from : " << inputs[inp_it]);
}
// Check for precomputed weights
TTree *weighttree = (TTree*)inp_file->Get("nova_wgts");
- if (!weighttree) {
- LOG(FIT) << "Did not find nova_wgts tree in file " << inputs[inp_it] << std::endl;
- nova_wgt = false;
- } else {
- LOG(FIT) << "Found nova_wgts tree in file " << inputs[inp_it] << std::endl;
- nova_wgt = true;
+ if (fNOvAWeights) {
+ if (!weighttree) {
+ THROW("Did not find nova_wgts tree in file " << inputs[inp_it] << " but you specified it" << std::endl);
+ } else {
+ LOG(FIT) << "Found nova_wgts tree in file " << inputs[inp_it] << std::endl;
+ }
}
// Register input to form flux/event rate hists
RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist);
// Add To TChain
fGENIETree->AddFile(inputs[inp_it].c_str());
if (weighttree != NULL) fGENIETree->AddFriend(weighttree);
}
// Registor all our file inputs
SetupJointInputs();
// Assign to tree
fEventType = kGENIE;
fGenieNtpl = NULL;
fGENIETree->SetBranchAddress("gmcrec", &fGenieNtpl);
// Set up the custom weights
- if (nova_wgt) {
+ if (fNOvAWeights) {
fGENIETree->SetBranchAddress("MAQEwgt", &MAQEw);
fGENIETree->SetBranchAddress("nonResNormWgt", &NonResw);
fGENIETree->SetBranchAddress("RPAQEWgt", &RPAQEw);
fGENIETree->SetBranchAddress("RPARESWgt", &RPARESw);
fGENIETree->SetBranchAddress("MECWgt", &MECw);
fGENIETree->SetBranchAddress("DISWgt", &DISw);
fGENIETree->SetBranchAddress("nova2018CVWgt", &NOVAw);
}
// Libraries should be seen but not heard...
StopTalking();
fGENIETree->GetEntry(0);
StartTalking();
// Create Fit Event
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->SetGenieEvent(fGenieNtpl);
if (fSaveExtra) {
fGenieInfo = new GENIEGeneratorInfo();
fNUISANCEEvent->AddGeneratorInfo(fGenieInfo);
}
fNUISANCEEvent->HardReset();
};
GENIEInputHandler::~GENIEInputHandler() {
// if (fGenieGHep) delete fGenieGHep;
// if (fGenieNtpl) delete fGenieNtpl;
// if (fGENIETree) delete fGENIETree;
// if (fGenieInfo) delete fGenieInfo;
}
void GENIEInputHandler::CreateCache() {
if (fCacheSize > 0) {
// fGENIETree->SetCacheEntryRange(0, fNEvents);
fGENIETree->AddBranchToCache("*", 1);
fGENIETree->SetCacheSize(fCacheSize);
}
}
void GENIEInputHandler::RemoveCache() {
// fGENIETree->SetCacheEntryRange(0, fNEvents);
fGENIETree->AddBranchToCache("*", 0);
fGENIETree->SetCacheSize(0);
}
FitEvent* GENIEInputHandler::GetNuisanceEvent(const UInt_t entry, const bool lightweight) {
if (entry >= (UInt_t)fNEvents) return NULL;
// Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
fGENIETree->GetEntry(entry);
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
#ifdef __PROB3PP_ENABLED__
else {
// Check for GENIE Event
if (!fGenieNtpl) return NULL;
if (!fGenieNtpl->event) return NULL;
// Cast Event Record
fGenieGHep = static_cast(fGenieNtpl->event);
if (!fGenieGHep) return NULL;
TObjArrayIter iter(fGenieGHep);
genie::GHepParticle* p;
while ((p = (dynamic_cast((iter).Next())))) {
if (!p) {
continue;
}
// Get Status
int state = GetGENIEParticleStatus(p, fNUISANCEEvent->Mode);
if (state != genie::kIStInitialState) {
continue;
}
fNUISANCEEvent->probe_E = p->E() * 1.E3;
fNUISANCEEvent->probe_pdg = p->Pdg();
break;
}
}
#endif
// Setup Input scaling for joint inputs
fNUISANCEEvent->InputWeight = GetInputWeight(entry);
return fNUISANCEEvent;
}
int GENIEInputHandler::GetGENIEParticleStatus(genie::GHepParticle* p, int mode) {
/*
kIStUndefined = -1,
kIStInitialState = 0, / generator-level initial state /
kIStStableFinalState = 1, / generator-level final state:
particles to be tracked by detector-level MC /
kIStIntermediateState = 2,
kIStDecayedState = 3,
kIStCorrelatedNucleon = 10,
kIStNucleonTarget = 11,
kIStDISPreFragmHadronicState = 12,
kIStPreDecayResonantState = 13,
kIStHadronInTheNucleus = 14, / hadrons inside the nucleus: marked
for hadron transport modules to act on /
kIStFinalStateNuclearRemnant = 15, / low energy nuclear fragments
entering the record collectively as a 'hadronic blob' pseudo-particle /
kIStNucleonClusterTarget = 16, // for composite nucleons before
phase space decay
*/
int state = kUndefinedState;
switch (p->Status()) {
case genie::kIStNucleonTarget:
case genie::kIStInitialState:
case genie::kIStCorrelatedNucleon:
case genie::kIStNucleonClusterTarget:
state = kInitialState;
break;
case genie::kIStStableFinalState:
state = kFinalState;
break;
case genie::kIStHadronInTheNucleus:
if (abs(mode) == 2)
state = kInitialState;
else
state = kFSIState;
break;
case genie::kIStPreDecayResonantState:
case genie::kIStDISPreFragmHadronicState:
case genie::kIStIntermediateState:
state = kFSIState;
break;
case genie::kIStFinalStateNuclearRemnant:
case genie::kIStUndefined:
case genie::kIStDecayedState:
default:
break;
}
// Flag to remove nuclear part in genie
if (p->Pdg() > 1000000) {
if (state == kInitialState)
state = kNuclearInitial;
else if (state == kFinalState)
state = kNuclearRemnant;
}
return state;
}
#endif
#ifdef __GENIE_ENABLED__
int GENIEInputHandler::ConvertGENIEReactionCode(GHepRecord* gheprec) {
// Electron Scattering
if (gheprec->Summary()->ProcInfo().IsEM()) {
if (gheprec->Summary()->InitState().ProbePdg() == 11) {
if (gheprec->Summary()->ProcInfo().IsQuasiElastic())
return 1;
else if (gheprec->Summary()->ProcInfo().IsMEC())
return 2;
else if (gheprec->Summary()->ProcInfo().IsResonant())
return 13;
else if (gheprec->Summary()->ProcInfo().IsDeepInelastic())
return 26;
else {
ERROR(WRN,
"Unknown GENIE Electron Scattering Mode!"
<< std::endl
<< "ScatteringTypeId = "
<< gheprec->Summary()->ProcInfo().ScatteringTypeId() << " "
<< "InteractionTypeId = "
<< gheprec->Summary()->ProcInfo().InteractionTypeId()
<< std::endl
<< genie::ScatteringType::AsString(
gheprec->Summary()->ProcInfo().ScatteringTypeId())
<< " "
<< genie::InteractionType::AsString(
gheprec->Summary()->ProcInfo().InteractionTypeId())
<< " " << gheprec->Summary()->ProcInfo().IsMEC());
return 0;
}
}
// Weak CC
} else if (gheprec->Summary()->ProcInfo().IsWeakCC()) {
// CC MEC
if (gheprec->Summary()->ProcInfo().IsMEC()) {
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return 2;
else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return -2;
// CC OTHER
} else {
return utils::ghep::NeutReactionCode(gheprec);
}
// Weak NC
} else if (gheprec->Summary()->ProcInfo().IsWeakNC()) {
// NC MEC
if (gheprec->Summary()->ProcInfo().IsMEC()) {
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return 32;
else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return -32;
// NC OTHER
} else {
return utils::ghep::NeutReactionCode(gheprec);
}
}
return 0;
}
void GENIEInputHandler::CalcNUISANCEKinematics() {
// Reset all variables
fNUISANCEEvent->ResetEvent();
// Check for GENIE Event
if (!fGenieNtpl) return;
if (!fGenieNtpl->event) return;
// Cast Event Record
fGenieGHep = static_cast(fGenieNtpl->event);
if (!fGenieGHep) return;
// Convert GENIE Reaction Code
fNUISANCEEvent->Mode = ConvertGENIEReactionCode(fGenieGHep);
// Set Event Info
fNUISANCEEvent->fEventNo = 0.0;
fNUISANCEEvent->fTotCrs = fGenieGHep->XSec();
fNUISANCEEvent->fTargetA = 0.0;
fNUISANCEEvent->fTargetZ = 0.0;
fNUISANCEEvent->fTargetH = 0;
fNUISANCEEvent->fBound = 0.0;
fNUISANCEEvent->InputWeight = 1.0; //(1E+38 / genie::units::cm2) * fGenieGHep->XSec();
// And the custom weights
- fNUISANCEEvent->CustomWeight = NOVAw;
- fNUISANCEEvent->CustomWeightArray[0] = MAQEw;
- fNUISANCEEvent->CustomWeightArray[1] = NonResw;
- fNUISANCEEvent->CustomWeightArray[2] = RPAQEw;
- fNUISANCEEvent->CustomWeightArray[3] = RPARESw;
- fNUISANCEEvent->CustomWeightArray[4] = MECw;
- fNUISANCEEvent->CustomWeightArray[5] = NOVAw;
-
- /*
- fNUISANCEEvent->CustomWeight = 1.0;
- fNUISANCEEvent->CustomWeightArray[0] = 1.0;
- fNUISANCEEvent->CustomWeightArray[1] = 1.0;
- fNUISANCEEvent->CustomWeightArray[2] = 1.0;
- fNUISANCEEvent->CustomWeightArray[3] = 1.0;
- fNUISANCEEvent->CustomWeightArray[4] = 1.0;
- fNUISANCEEvent->CustomWeightArray[5] = 1.0;
- */
+ if (fNOvAWeights) {
+ fNUISANCEEvent->CustomWeight = NOVAw;
+ fNUISANCEEvent->CustomWeightArray[0] = MAQEw;
+ fNUISANCEEvent->CustomWeightArray[1] = NonResw;
+ fNUISANCEEvent->CustomWeightArray[2] = RPAQEw;
+ fNUISANCEEvent->CustomWeightArray[3] = RPARESw;
+ fNUISANCEEvent->CustomWeightArray[4] = MECw;
+ fNUISANCEEvent->CustomWeightArray[5] = NOVAw;
+ } else {
+ fNUISANCEEvent->CustomWeight = 1.0;
+ fNUISANCEEvent->CustomWeightArray[0] = 1.0;
+ fNUISANCEEvent->CustomWeightArray[1] = 1.0;
+ fNUISANCEEvent->CustomWeightArray[2] = 1.0;
+ fNUISANCEEvent->CustomWeightArray[3] = 1.0;
+ fNUISANCEEvent->CustomWeightArray[4] = 1.0;
+ fNUISANCEEvent->CustomWeightArray[5] = 1.0;
+ }
// Get N Particle Stack
unsigned int npart = fGenieGHep->GetEntries();
unsigned int kmax = fNUISANCEEvent->kMaxParticles;
if (npart > kmax) {
ERR(WRN) << "GENIE has too many particles, expanding stack." << std::endl;
fNUISANCEEvent->ExpandParticleStack(npart);
}
// Fill Particle Stack
GHepParticle* p = 0;
TObjArrayIter iter(fGenieGHep);
fNUISANCEEvent->fNParticles = 0;
// Loop over all particles
while ((p = (dynamic_cast((iter).Next())))) {
if (!p) continue;
// Get Status
int state = GetGENIEParticleStatus(p, fNUISANCEEvent->Mode);
// Remove Undefined
if (kRemoveUndefParticles && state == kUndefinedState) continue;
// Remove FSI
if (kRemoveFSIParticles && state == kFSIState) continue;
if (kRemoveNuclearParticles &&
(state == kNuclearInitial || state == kNuclearRemnant))
continue;
// Fill Vectors
int curpart = fNUISANCEEvent->fNParticles;
fNUISANCEEvent->fParticleState[curpart] = state;
// Mom
fNUISANCEEvent->fParticleMom[curpart][0] = p->Px() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][1] = p->Py() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][2] = p->Pz() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][3] = p->E() * 1.E3;
// PDG
fNUISANCEEvent->fParticlePDG[curpart] = p->Pdg();
// Add to N particle count
fNUISANCEEvent->fNParticles++;
// Extra Check incase GENIE fails.
if ((UInt_t)fNUISANCEEvent->fNParticles == kmax) {
ERR(WRN) << "Number of GENIE Particles exceeds maximum!" << std::endl;
ERR(WRN) << "Extend kMax, or run without including FSI particles!"
<< std::endl;
break;
}
}
// Fill Extra Stack
if (fSaveExtra) fGenieInfo->FillGeneratorInfo(fGenieNtpl);
// Run Initial, FSI, Final, Other ordering.
fNUISANCEEvent->OrderStack();
FitParticle* ISNeutralLepton =
fNUISANCEEvent->GetHMISParticle(PhysConst::pdg_neutrinos);
if (ISNeutralLepton) {
fNUISANCEEvent->probe_E = ISNeutralLepton->E();
fNUISANCEEvent->probe_pdg = ISNeutralLepton->PDG();
}
return;
}
void GENIEInputHandler::Print() {}
#endif
diff --git a/src/InputHandler/GENIEInputHandler.h b/src/InputHandler/GENIEInputHandler.h
index 818fda8..d2a1d18 100644
--- a/src/InputHandler/GENIEInputHandler.h
+++ b/src/InputHandler/GENIEInputHandler.h
@@ -1,117 +1,119 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifndef GENIEINPUTHANDLER_H
#define GENIEINPUTHANDLER_H
/*!
* \addtogroup InputHandler
* @{
*/
#ifdef __GENIE_ENABLED__
#include "InputHandler.h"
#include "InputUtils.h"
#include "PlotUtils.h"
#include "GHEP/GHepParticle.h"
#include "PDG/PDGUtils.h"
#include "GHEP/GHepUtils.h"
#include "Conventions/Units.h"
#include "EVGCore/EventRecord.h"
#include "GHEP/GHepRecord.h"
#include "Ntuple/NtpMCEventRecord.h"
using namespace genie;
/// GENIE Generator Container to save extra particle status codes.
class GENIEGeneratorInfo : public GeneratorInfoBase {
public:
GENIEGeneratorInfo() {};
virtual ~GENIEGeneratorInfo();
/// Assigns information to branches
void AddBranchesToTree(TTree* tn);
/// Setup reading information from branches
void SetBranchesFromTree(TTree* tn);
/// Allocate any dynamic arrays for a new particle stack size
void AllocateParticleStack(int stacksize);
/// Clear any dynamic arrays
void DeallocateParticleStack();
/// Read extra genie information from the event
void FillGeneratorInfo(NtpMCEventRecord* ntpl);
/// Reset extra information to default/empty values
void Reset();
int kMaxParticles; ///< Number of particles in stack
int* fGenieParticlePDGs; ///< GENIE Particle PDGs (example)
};
/// Main GENIE InputHandler
class GENIEInputHandler : public InputHandlerBase {
public:
/// Standard constructor given a name and input files
GENIEInputHandler(std::string const& handle, std::string const& rawinputs);
virtual ~GENIEInputHandler();
/// Create a TTree Cache to speed up file read
void CreateCache();
/// Remove TTree Cache to save memory
void RemoveCache();
/// Returns a NUISANCE format event from the GENIE TTree. If !lightweight
/// then CalcNUISANCEKinematics() is called to convert the GENIE event into
/// a standard NUISANCE format.
FitEvent* GetNuisanceEvent(const UInt_t entry, const bool lightweight = false);
/// Converts GENIE event into standard NUISANCE FitEvent by looping over all
/// particles in the event and adding them to stack in fNUISANCEEvent.
void CalcNUISANCEKinematics();
/// Placeholder for GENIE related event printing.
void Print();
/// Converts GENIE particle status codes into NUISANCE status codes.
int GetGENIEParticleStatus(genie::GHepParticle* part, int mode = 0);
/// Converts GENIE event reaction codes into NUISANCE reaction codes.
int ConvertGENIEReactionCode(GHepRecord* gheprec);
GHepRecord* fGenieGHep; ///< Pointer to actual event record
NtpMCEventRecord* fGenieNtpl; ///< Ntpl Wrapper Class
TChain* fGENIETree; ///< Main GENIE Event TTree
bool fSaveExtra; ///< Flag to save Extra GENIE info into Nuisance Event
GENIEGeneratorInfo* fGenieInfo; ///< Extra GENIE Generator Info Writer
- // Extra weights from Jeremgy
+ bool fNOvAWeights; ///< Flag to save nova weights or not
+
+ // Extra weights from Jeremy for NOvA weights
double MAQEw;
double NonResw;
double RPAQEw;
double RPARESw;
double MECw;
double DISw;
double NOVAw;
};
/*! @} */
#endif
#endif
diff --git a/src/MCStudies/GenericFlux_Vectors.cxx b/src/MCStudies/GenericFlux_Vectors.cxx
index 7b65fd4..045f0c0 100644
--- a/src/MCStudies/GenericFlux_Vectors.cxx
+++ b/src/MCStudies/GenericFlux_Vectors.cxx
@@ -1,337 +1,339 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "GenericFlux_Vectors.h"
GenericFlux_Vectors::GenericFlux_Vectors(std::string name,
std::string inputfile, FitWeight *rw,
std::string type,
std::string fakeDataFile) {
// Measurement Details
fName = name;
eventVariables = NULL;
// Define our energy range for flux calcs
EnuMin = 0.;
EnuMax = 1E10; // Arbritrarily high energy limit
if (Config::HasPar("EnuMin")) {
EnuMin = Config::GetParD("EnuMin");
}
if (Config::HasPar("EnuMax")) {
EnuMax = Config::GetParD("EnuMax");
}
// Set default fitter flags
fIsDiag = true;
fIsShape = false;
fIsRawEvents = false;
// This function will sort out the input files automatically and parse all the
// inputs,flags,etc.
// There may be complex cases where you have to do this by hand, but usually
// this will do.
Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile);
eventVariables = NULL;
// Setup fDataHist as a placeholder
this->fDataHist = new TH1D(("empty_data"), ("empty-data"), 1, 0, 1);
this->SetupDefaultHist();
fFullCovar = StatUtils::MakeDiagonalCovarMatrix(fDataHist);
covar = StatUtils::GetInvert(fFullCovar);
// 1. The generator is organised in SetupMeasurement so it gives the
// cross-section in "per nucleon" units.
// So some extra scaling for a specific measurement may be required. For
// Example to get a "per neutron" measurement on carbon
// which we do here, we have to multiple by the number of nucleons 12 and
// divide by the number of neutrons 6.
// N.B. MeasurementBase::PredictedEventRate includes the 1E-38 factor that is
// often included here in other classes that directly integrate the event
// histogram. This method is used here as it now respects EnuMin and EnuMax
// correctly.
this->fScaleFactor =
(this->PredictedEventRate("width", 0, EnuMax) / double(fNEvents)) /
this->TotalIntegratedFlux();
LOG(SAM) << " Generic Flux Scaling Factor = " << fScaleFactor
<< " [= " << (GetEventHistogram()->Integral("width") * 1E-38) << "/("
<< (fNEvents + 0.) << "*" << this->TotalIntegratedFlux() << ")]"
<< std::endl;
if (fScaleFactor <= 0.0) {
ERR(WRN) << "SCALE FACTOR TOO LOW " << std::endl;
throw;
}
// Setup our TTrees
this->AddEventVariablesToTree();
this->AddSignalFlagsToTree();
}
void GenericFlux_Vectors::AddEventVariablesToTree() {
// Setup the TTree to save everything
if (!eventVariables) {
Config::Get().out->cd();
eventVariables = new TTree((this->fName + "_VARS").c_str(),
(this->fName + "_VARS").c_str());
}
LOG(SAM) << "Adding Event Variables" << std::endl;
eventVariables->Branch("Mode", &Mode, "Mode/I");
eventVariables->Branch("cc", &cc, "cc/B");
eventVariables->Branch("PDGnu", &PDGnu, "PDGnu/I");
eventVariables->Branch("Enu_true", &Enu_true, "Enu_true/F");
eventVariables->Branch("tgt", &tgt, "tgt/I");
eventVariables->Branch("PDGLep", &PDGLep, "PDGLep/I");
eventVariables->Branch("ELep", &ELep, "ELep/F");
eventVariables->Branch("CosLep", &CosLep, "CosLep/F");
// Basic interaction kinematics
eventVariables->Branch("Q2", &Q2, "Q2/F");
eventVariables->Branch("q0", &q0, "q0/F");
eventVariables->Branch("q3", &q3, "q3/F");
eventVariables->Branch("Enu_QE", &Enu_QE, "Enu_QE/F");
eventVariables->Branch("Q2_QE", &Q2_QE, "Q2_QE/F");
eventVariables->Branch("W_nuc_rest", &W_nuc_rest, "W_nuc_rest/F");
eventVariables->Branch("W", &W, "W/F");
eventVariables->Branch("x", &x, "x/F");
eventVariables->Branch("y", &y, "y/F");
eventVariables->Branch("Eav", &Eav, "Eav/F");
+ eventVariables->Branch("EavAlt", &EavAlt, "EavAlt/F");
// Save outgoing particle vectors
eventVariables->Branch("nfsp", &nfsp, "nfsp/I");
eventVariables->Branch("px", px, "px[nfsp]/F");
eventVariables->Branch("py", py, "py[nfsp]/F");
eventVariables->Branch("pz", pz, "pz[nfsp]/F");
eventVariables->Branch("E", E, "E[nfsp]/F");
eventVariables->Branch("pdg", pdg, "pdg[nfsp]/I");
// Event Scaling Information
eventVariables->Branch("Weight", &Weight, "Weight/F");
eventVariables->Branch("InputWeight", &InputWeight, "InputWeight/F");
eventVariables->Branch("RWWeight", &RWWeight, "RWWeight/F");
// Should be a double because may be 1E-39 and less
eventVariables->Branch("fScaleFactor", &fScaleFactor, "fScaleFactor/D");
// The customs
eventVariables->Branch("CustomWeight", &CustomWeight, "CustomWeight/F");
eventVariables->Branch("CustomWeightArray", CustomWeightArray, "CustomWeightArray[6]/F");
return;
}
void GenericFlux_Vectors::FillEventVariables(FitEvent *event) {
ResetVariables();
// Fill Signal Variables
FillSignalFlags(event);
LOG(DEB) << "Filling signal" << std::endl;
// Now fill the information
Mode = event->Mode;
cc = (abs(event->Mode) < 30);
// Get the incoming neutrino and outgoing lepton
FitParticle *nu = event->GetNeutrinoIn();
FitParticle *lep = event->GetHMFSAnyLepton();
PDGnu = nu->fPID;
Enu_true = nu->fP.E() / 1E3;
tgt = event->fTargetPDG;
if (lep != NULL) {
PDGLep = lep->fPID;
ELep = lep->fP.E() / 1E3;
CosLep = cos(nu->fP.Vect().Angle(lep->fP.Vect()));
// Basic interaction kinematics
Q2 = -1 * (nu->fP - lep->fP).Mag2() / 1E6;
q0 = (nu->fP - lep->fP).E() / 1E3;
q3 = (nu->fP - lep->fP).Vect().Mag() / 1E3;
// These assume C12 binding from MINERvA... not ideal
Enu_QE = FitUtils::EnuQErec(lep->fP, CosLep, 34., true);
Q2_QE = FitUtils::Q2QErec(lep->fP, CosLep, 34., true);
Eav = FitUtils::GetErecoil_MINERvA_LowRecoil(event)/1.E3;
+ EavAlt = FitUtils::Eavailable(event)/1.E3;
// Get W_true with assumption of initial state nucleon at rest
float m_n = (float)PhysConst::mass_proton;
// Q2 assuming nucleon at rest
W_nuc_rest = sqrt(-Q2 + 2 * m_n * q0 + m_n * m_n);
// True Q2
W = sqrt(-Q2 + 2 * m_n * q0 + m_n * m_n);
x = Q2 / (2 * m_n * q0);
y = 1 - ELep / Enu_true;
}
// Loop over the particles and store all the final state particles in a vector
for (UInt_t i = 0; i < event->Npart(); ++i) {
bool part_alive = event->PartInfo(i)->fIsAlive &&
event->PartInfo(i)->Status() == kFinalState;
if (!part_alive) continue;
partList.push_back(event->PartInfo(i));
}
// Save outgoing particle vectors
nfsp = (int)partList.size();
for (int i = 0; i < nfsp; ++i) {
px[i] = partList[i]->fP.X() / 1E3;
py[i] = partList[i]->fP.Y() / 1E3;
pz[i] = partList[i]->fP.Z() / 1E3;
E[i] = partList[i]->fP.E() / 1E3;
pdg[i] = partList[i]->fPID;
}
// Fill event weights
Weight = event->RWWeight * event->InputWeight;
RWWeight = event->RWWeight;
InputWeight = event->InputWeight;
// And the Customs
CustomWeight = event->CustomWeight;
for (int i = 0; i < 6; ++i) {
CustomWeightArray[i] = event->CustomWeightArray[i];
}
// Fill the eventVariables Tree
eventVariables->Fill();
return;
};
//********************************************************************
void GenericFlux_Vectors::ResetVariables() {
//********************************************************************
cc = false;
// Reset all Function used to extract any variables of interest to the event
Mode = PDGnu = tgt = PDGLep = 0;
- Enu_true = ELep = CosLep = Q2 = q0 = q3 = Enu_QE = Q2_QE = W_nuc_rest = W = x = y = Eav = -999.9;
+ Enu_true = ELep = CosLep = Q2 = q0 = q3 = Enu_QE = Q2_QE = W_nuc_rest = W = x = y = Eav = EavAlt = -999.9;
nfsp = 0;
for (int i = 0; i < kMAX; ++i){
px[i] = py[i] = pz[i] = E[i] = -999;
pdg[i] = 0;
}
Weight = InputWeight = RWWeight = 0.0;
CustomWeight = 0.0;
for (int i = 0; i < 6; ++i) CustomWeightArray[i] = 0.0;
partList.clear();
flagCCINC = flagNCINC = flagCCQE = flagCC0pi = flagCCQELike = flagNCEL = flagNC0pi = flagCCcoh = flagNCcoh = flagCC1pip = flagNC1pip = flagCC1pim = flagNC1pim = flagCC1pi0 = flagNC1pi0 = false;
}
//********************************************************************
void GenericFlux_Vectors::FillSignalFlags(FitEvent *event) {
//********************************************************************
// Some example flags are given from SignalDef.
// See src/Utils/SignalDef.cxx for more.
int nuPDG = event->PartInfo(0)->fPID;
// Generic signal flags
flagCCINC = SignalDef::isCCINC(event, nuPDG);
flagNCINC = SignalDef::isNCINC(event, nuPDG);
flagCCQE = SignalDef::isCCQE(event, nuPDG);
flagCCQELike = SignalDef::isCCQELike(event, nuPDG);
flagCC0pi = SignalDef::isCC0pi(event, nuPDG);
flagNCEL = SignalDef::isNCEL(event, nuPDG);
flagNC0pi = SignalDef::isNC0pi(event, nuPDG);
flagCCcoh = SignalDef::isCCCOH(event, nuPDG, 211);
flagNCcoh = SignalDef::isNCCOH(event, nuPDG, 111);
flagCC1pip = SignalDef::isCC1pi(event, nuPDG, 211);
flagNC1pip = SignalDef::isNC1pi(event, nuPDG, 211);
flagCC1pim = SignalDef::isCC1pi(event, nuPDG, -211);
flagNC1pim = SignalDef::isNC1pi(event, nuPDG, -211);
flagCC1pi0 = SignalDef::isCC1pi(event, nuPDG, 111);
flagNC1pi0 = SignalDef::isNC1pi(event, nuPDG, 111);
}
void GenericFlux_Vectors::AddSignalFlagsToTree() {
if (!eventVariables) {
Config::Get().out->cd();
eventVariables = new TTree((this->fName + "_VARS").c_str(),
(this->fName + "_VARS").c_str());
}
LOG(SAM) << "Adding signal flags" << std::endl;
// Signal Definitions from SignalDef.cxx
eventVariables->Branch("flagCCINC", &flagCCINC, "flagCCINC/O");
eventVariables->Branch("flagNCINC", &flagNCINC, "flagNCINC/O");
eventVariables->Branch("flagCCQE", &flagCCQE, "flagCCQE/O");
eventVariables->Branch("flagCC0pi", &flagCC0pi, "flagCC0pi/O");
eventVariables->Branch("flagCCQELike", &flagCCQELike, "flagCCQELike/O");
eventVariables->Branch("flagNCEL", &flagNCEL, "flagNCEL/O");
eventVariables->Branch("flagNC0pi", &flagNC0pi, "flagNC0pi/O");
eventVariables->Branch("flagCCcoh", &flagCCcoh, "flagCCcoh/O");
eventVariables->Branch("flagNCcoh", &flagNCcoh, "flagNCcoh/O");
eventVariables->Branch("flagCC1pip", &flagCC1pip, "flagCC1pip/O");
eventVariables->Branch("flagNC1pip", &flagNC1pip, "flagNC1pip/O");
eventVariables->Branch("flagCC1pim", &flagCC1pim, "flagCC1pim/O");
eventVariables->Branch("flagNC1pim", &flagNC1pim, "flagNC1pim/O");
eventVariables->Branch("flagCC1pi0", &flagCC1pi0, "flagCC1pi0/O");
eventVariables->Branch("flagNC1pi0", &flagNC1pi0, "flagNC1pi0/O");
};
void GenericFlux_Vectors::Write(std::string drawOpt) {
// First save the TTree
eventVariables->Write();
// Save Flux and Event Histograms too
GetInput()->GetFluxHistogram()->Write();
GetInput()->GetEventHistogram()->Write();
return;
}
// Override functions which aren't really necessary
bool GenericFlux_Vectors::isSignal(FitEvent *event) {
(void)event;
return true;
};
void GenericFlux_Vectors::ScaleEvents() { return; }
void GenericFlux_Vectors::ApplyNormScale(float norm) {
this->fCurrentNorm = norm;
return;
}
void GenericFlux_Vectors::FillHistograms() { return; }
void GenericFlux_Vectors::ResetAll() {
eventVariables->Reset();
return;
}
float GenericFlux_Vectors::GetChi2() { return 0.0; }
diff --git a/src/MCStudies/GenericFlux_Vectors.h b/src/MCStudies/GenericFlux_Vectors.h
index 7b76a57..921b247 100644
--- a/src/MCStudies/GenericFlux_Vectors.h
+++ b/src/MCStudies/GenericFlux_Vectors.h
@@ -1,126 +1,127 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifndef GenericFlux_Vectors_H_SEEN
#define GenericFlux_Vectors_H_SEEN
#include "Measurement1D.h"
class GenericFlux_Vectors : public Measurement1D {
public:
GenericFlux_Vectors(std::string name, std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile);
virtual ~GenericFlux_Vectors() {};
//! Grab info from event
void FillEventVariables(FitEvent *event);
//! Fill signal flags
void FillSignalFlags(FitEvent *event);
void ResetVariables();
//! Fill Custom Histograms
void FillHistograms();
//! ResetAll
void ResetAll();
//! Scale
void ScaleEvents();
//! Norm
void ApplyNormScale(float norm);
//! Define this samples signal
bool isSignal(FitEvent *nvect);
//! Write Files
void Write(std::string drawOpt);
//! Get Chi2
float GetChi2();
void AddEventVariablesToTree();
void AddSignalFlagsToTree();
private:
TTree* eventVariables;
std::vector partList;
int Mode;
bool cc;
int PDGnu;
int tgt;
int PDGLep;
float ELep;
float CosLep;
// Basic interaction kinematics
float Q2;
float q0;
float q3;
float Enu_QE;
float Enu_true;
float Q2_QE;
float W_nuc_rest;
float W;
float x;
float y;
float Eav;
+ float EavAlt;
// Save outgoing particle vectors
int nfsp;
static const int kMAX = 200;
float px[kMAX];
float py[kMAX];
float pz[kMAX];
float E[kMAX];
int pdg[kMAX];
// Basic event info
float Weight;
float InputWeight;
float RWWeight;
double fScaleFactor;
// Custom weights
float CustomWeight;
float CustomWeightArray[6];
// Generic signal flags
bool flagCCINC;
bool flagNCINC;
bool flagCCQE;
bool flagCC0pi;
bool flagCCQELike;
bool flagNCEL;
bool flagNC0pi;
bool flagCCcoh;
bool flagNCcoh;
bool flagCC1pip;
bool flagNC1pip;
bool flagCC1pim;
bool flagNC1pim;
bool flagCC1pi0;
bool flagNC1pi0;
};
#endif
diff --git a/src/Reweight/weightRPA.h b/src/Reweight/weightRPA.h
index c44acdd..ca8e635 100644
--- a/src/Reweight/weightRPA.h
+++ b/src/Reweight/weightRPA.h
@@ -1,409 +1,419 @@
#ifndef weightRPA_h
#define weightRPA_h
#include //ifstream
#include //cout
#include
#include
#include
#include
#include "math.h"
#include "assert.h"
//For Compatibility with ROOT compiler
//uncomment the following:
//#define ROOT
/*!
* Code example to extract the RPA effect central weight
* and its uncertainties from the prepared files.
* Heidi Schellman (Oregon State) and Rik Gran (Minnesota Duluth)
* for use in MINERvA experiment analysis
* must compile with the ROOT libraries
* g++ `root-config --glibs --cflags` -O3 weightRPAtest.cxx -o weightRPA
* The underlying model is from the IFIC Valencia group
* see (public) minerva docdb:12688 for the full physics discussion
*/
// NOTE UNITS ARE GEV in the calculation
// make sure you convert MeV to GeV before calling these functions
// Class for getting RPA paramaters inputs from a given file
// Includes methods that return all five RPA weights at once
// (is the most cpu efficient way to get them)
// Or return just the central value
// (skipping the uncertainties code completely if only cv wanted)
// Or return each CV and uncertainty one at a time
// (is the least cpu efficient, repeats some calculations 5 times)
class weightRPA {
public:
//Constructor: Read in params from a filename
weightRPA(const TString f) { read(f); } //Read in params from file
TString filename;
TFile* fRPAratio;
TH2D *hRPArelratio;
TH2D *hRPAnonrelratio;
TH1D *hQ2relratio;
TH1D *hQ2nonrelratio;
TArrayD *TADrpapolyrel;
Double_t *rpapolyrel ;
TArrayD *TADrpapolynonrel;
Double_t *rpapolynonrel;
static const int CENTRAL=0;
static const int LOWQ2 = 1;
static const int HIGHQ2 = 2;
// true to take from histogram, false use parameterization
static const bool Q2histORparam = true;
// MINERvA holds kinematics in MeV, but all these functions require GeV
// So make sure you pass them in GeV.
inline double getWeight(const double mc_q0, const double mc_q3, double * weights); //in GeV
inline double getWeight(const double mc_q0, const double mc_q3); //in GeV
inline double getWeight(const double mc_q0, const double mc_q3, int type, int sign); //in GeV
inline double getWeight(const double Q2); //in GeV^2
inline double getWeightLowQ2(const double mc_q0, const double mc_q3, const int sign);
inline double getWeightHighQ2(const double mc_q0, const double mc_q3, const int sign);
inline double getWeightQ2(const double mc_Q2, const bool relORnonrel=true);
//Initializer
inline void read(const TString f);
// q0 and q3 in GeV, type = 1 for low Q2, 2 for high Q2, 0 for central
//double getWeightInternal(const double mc_q0, const double mc_q3,int type, int sign);
private:
inline double getWeightInternal(const double mc_Q2);
inline double getWeightInternal(const double mc_q0, const double mc_q3, const int type, const int sign);
inline double getWeightInternal(const double mc_q0, const double mc_q3, double *weights=0);
inline double getWeightQ2parameterization(const double mc_Q2, const bool relORnonrel);
inline double getWeightQ2fromhistogram(const double mc_Q2, const bool relORnonrel);
};
void weightRPA::read(const TString f)
//Read in the params doubles from a file
//argument: valid filename
{
fRPAratio = TFile::Open(f,"READONLY");
if (fRPAratio){
hRPArelratio = (TH2D*)fRPAratio->Get("hrelratio");
hRPAnonrelratio = (TH2D*)fRPAratio->Get("hnonrelratio");
hQ2relratio = (TH1D*)fRPAratio->Get("hQ2relratio");
hQ2nonrelratio = (TH1D*)fRPAratio->Get("hQ2nonrelratio");
TADrpapolyrel = (TArrayD*)fRPAratio->Get("rpapolyrel");
rpapolyrel = TADrpapolyrel->GetArray();
TADrpapolynonrel = (TArrayD*)fRPAratio->Get("rpapolynonrel");
rpapolynonrel = TADrpapolynonrel->GetArray();
hRPArelratio->Print();
std::cout << "have read in ratios from file " << f <= gevlimit) q0bin = rpamevlimit - 1;
if(mc_q3 >= gevlimit) q3bin = rpamevlimit - 1;
// Nieves does not calculate anything below binding energy.
// I don't know what GENIE does, but lets be soft about this.
// Two things lurking here at once.
// One, we are taking out a 10 MeV offset between GENIE and Valencia.
// Second, we are protecting against being asked for a weight that is too low in q0.
// It actually shouldn't happen for real GENIE events,
// but this protection does something that doesn't suck, just in case.
// you would see the artifact in a plot for sure, but better than writing 1.0.
+
+ // CWret after talking to Rik in Nov 2018 at MINERvA CM
+ // 2.12.8 sets this to 27 because change of Q definition: need to offset GENIE and Nieves Eb even more
+#if __GENIE_VERSION__ >= 210
+ Int_t q0offsetValenciaGENIE = 27;
+#else
Int_t q0offsetValenciaGENIE = 10;
+#endif
+ std::cout << "q0offsetValenciaGENIE: " << q0offsetValenciaGENIE << std::endl;
+
+
if(mc_q0 < 0.018) q0bin = 18+q0offsetValenciaGENIE;
Double_t thisrwtemp = hRPArelratio->GetBinContent(q3bin,q0bin-q0offsetValenciaGENIE);
// now trap bogus entries. Not sure why they happen, but set to 1.0 not 0.0
if(thisrwtemp <= 0.001)thisrwtemp = 1.0;
// events in genie but not in valencia should get a weight
// related to a similar q0 from the bulk distribution.
if(mc_q0 < 0.15 && thisrwtemp > 0.9){
thisrwtemp = hRPArelratio->GetBinContent(q3bin+150, q0bin-q0offsetValenciaGENIE);
}
//Double_t *mypoly;
//mypoly = rpapoly;
if(Q2gev >= 9.0){
thisrwtemp = 1.0;
}
else if(Q2gev > 3.0) {
// hiding option, use the Q2 parameterization everywhere
// } else if(Q2gev > 3.0 || rwRPAQ2) {
// turn rwRPAQ2 all the way on to override the 2D histogram
// illustrates the old-style Q2 suppression many folks still use.
thisrwtemp = getWeightQ2(Q2gev,true);
// double powerQ2 = 1.0;
//thisrwtemp = 0.0;
//for(int ii=0; ii<10; ii++){
// thisrwtemp += rpapoly[ii]*powerQ2;
// powerQ2 *= Q2gev;
//}
//std::cout << "test temp " << thisrwtemp << " " << rpamypoly[2] << std::endl;
}
if(!(thisrwtemp >= 0.001 && thisrwtemp <= 2.0))thisrwtemp = 1.0;
// hiding option, turn off the enhancement.
//if(rwoffSRC && thisrwtemp > 1.0)thisrwtemp = 1.0;
if (0 == weights) return thisrwtemp;
// if this was called without passing an array,
// the user didn't want us to calculate the +/- 1-sigma bounds
// so the above line returned before even trying.
weights[0] = thisrwtemp;
//if (type == 0) return thisrwtemp;
//if (type == 1) {
// Construct the error bands on the low Q2 suppression.
// Double_t thisrwSupP1 = 1.0;
// Double_t thisrwSupM1 = 1.0;
if( thisrwtemp < 1.0){
// make the suppression stronger or weaker to muon capture uncertainty
// rwRPAonesig is either +1 or -1, which is 0.25 (25%).
// possible to be re-written to produce 2 and 3 sigma.
weights[1] = thisrwtemp + 1.0 * (0.25)*(1.0 - thisrwtemp);
weights[2] = thisrwtemp - 1.0 * (0.25)*(1.0 - thisrwtemp);
}
else{
weights[1] = thisrwtemp;
weights[2] = thisrwtemp;
}
//std::cout << "check " << thisrwtemp << " " << weights[1] << " " << weights[2] << std::endl;
// Construct the rest of the error bands on the low Q2 suppression.
// this involves getting the weight from the non-relativistic ratio
//if (type == 2){
Double_t thisrwEnhP1 = 1.0;
Double_t thisrwEnhM1 = 1.0;
// make enhancement stronger or weaker to Federico Sanchez uncertainty
// this does NOT mean two sigma, its overloading the option.
Double_t thisrwextreme = hRPAnonrelratio->GetBinContent(q3bin,q0bin-q0offsetValenciaGENIE);
// now trap bogus entries. Not sure why they happen, but set to 1.0 not 0.0
if(thisrwextreme <= 0.001)thisrwextreme = 1.0;
if(mc_q0 < 0.15 && thisrwextreme > 0.9){
thisrwextreme = hRPAnonrelratio->GetBinContent(q3bin+150, q0bin-q0offsetValenciaGENIE);
}
//std::cout << "ext " << thisrwextreme << " " << thisrwtemp << std::endl;
// get the same for the Q2 dependent thing,
// but from the nonrelativistic polynomial
if(Q2gev >= 9.0){
thisrwextreme = 1.0;
}
else if(Q2gev > 3.0 ) {
thisrwextreme = getWeightQ2(Q2gev,false);
//double powerQ2 = 1.0;
//thisrwextreme = 0.0;
//for(int ii=0; ii<10; ii++){
// thisrwextreme += rpapolynonrel[ii]*powerQ2;
// powerQ2 *= Q2gev;
//}
//std::cout << "test extreme " << thisrwextreme << " " << mypolynonrel[2] << std::endl;
}
if(!(thisrwextreme >= 0.001 && thisrwextreme <= 2.0))thisrwextreme = 1.0;
//std::cout << "test extreme " << Q2gev << " " << thisrwextreme << " " << thisrwtemp << std::endl;
Double_t RelToNonRel = 0.6;
// based on some distance between central value and extreme
thisrwEnhP1 = thisrwtemp + RelToNonRel * (thisrwextreme-thisrwtemp);
Double_t thisrwEnhP1max = thisrwextreme;
if(Q2gev < 0.9)thisrwEnhP1 += 1.5*(0.9 - Q2gev)*(thisrwEnhP1max - thisrwEnhP1);
// sanity check, don't let the upper error bound go above the nonrel limit.
if(thisrwEnhP1 > thisrwEnhP1max)thisrwEnhP1 = thisrwEnhP1max;
// don't let positive error bound be closer than 3% above the central value
// will happen at very high Q2 and very close to Q2 = 0
if(thisrwEnhP1 < thisrwtemp + 0.03)thisrwEnhP1 = thisrwtemp + 0.03;
thisrwEnhM1 = thisrwtemp - RelToNonRel * (thisrwextreme-thisrwtemp);
// don't let negative error bound be closer than 3% below the central value
if(thisrwEnhM1 > thisrwtemp - 0.03)thisrwEnhM1 = thisrwtemp - 0.03;
// even still, don't let the lower error bound go below 1.0 at high-ish Q2
if(Q2gev > 1.0 && thisrwEnhM1 < 1.0)thisrwEnhM1 = 1.0;
// whew. so now return the main weight
// and return the array of all five weights in some array
// thisrwtemp, thisrwSupP1, thisrwSupM1, thisrwEnhP1, thisrwEnhM1
//if (sign == 1) return thisrwEnhP1;
//if (sign == -1) return thisrwEnhM1;
weights[3] = thisrwEnhP1;
weights[4] = thisrwEnhM1;
// still return the central value
return thisrwtemp;
}
double weightRPA::getWeight(const double mc_q0, const double mc_q3){
return getWeightInternal(mc_q0, mc_q3);
}
double weightRPA::getWeight(const double mc_q0, const double mc_q3, double *weights){
return getWeightInternal(mc_q0, mc_q3, weights);
}
double weightRPA::getWeight(const double mc_q0, const double mc_q3, int type, int sign){
return getWeightInternal(mc_q0, mc_q3, type, sign);
}
double weightRPA::getWeight(const double mc_Q2){
return getWeightQ2(mc_Q2);
}
double weightRPA::getWeightInternal(const double mc_q0, const double mc_q3, int type, int sign){
double weights[5] = {1., 1., 1., 1., 1.};
double cv = getWeightInternal(mc_q0, mc_q3, weights);
if(type==0)return cv;
else if(type==weightRPA::LOWQ2 && sign == 1)return weights[1];
else if(type==weightRPA::LOWQ2 && sign == -1)return weights[2];
else if(type==weightRPA::HIGHQ2 && sign == 1)return weights[3];
else if(type==weightRPA::HIGHQ2 && sign == -1)return weights[4];
//else {
// // should never happen? Bork ?
// return cv; //?
//}
return cv;
}
double weightRPA::getWeightQ2(const double mc_Q2, const bool relORnonrel){
if(mc_Q2 < 0.0)return 1.0; // this is Q2 actually, not sure hw
if(mc_Q2 > 9.0)return 1.0;
// this function needs to know two options.
// does user want rel (cv) or nonrel
// does user want to use the histogram or parameterization
if(Q2histORparam)return getWeightQ2fromhistogram(mc_Q2, relORnonrel);
else return getWeightQ2parameterization(mc_Q2, relORnonrel);
}
double weightRPA::getWeightQ2parameterization(const double mc_Q2, const bool relORnonrel){
if(mc_Q2 < 0.0)return 1.0;
if(mc_Q2 > 9.0)return 1.0;
// this one returns just the polynomial Q2 version
// for special tests. Poor answer for baseline MINERvA QE events.
// No uncertainty assigned to this usecase at this time.
//double gevmev = 0.001; // minerva sends in MeV.
double Q2gev = mc_Q2;
double powerQ2 = 1.0;
double thisrwtemp = 0.0;
thisrwtemp = 0.0;
for(int ii=0; ii<10; ii++){
if(relORnonrel)thisrwtemp += rpapolyrel[ii]*powerQ2;
else thisrwtemp += rpapolynonrel[ii]*powerQ2;
powerQ2 *= Q2gev;
}
return thisrwtemp;
}
double weightRPA::getWeightQ2fromhistogram(const double mc_Q2, const bool relORnonrel){
if(mc_Q2 < 0.0)return 1.0;
if(mc_Q2 > 9.0) return 1.0;
if(relORnonrel)return hQ2relratio->GetBinContent( hQ2relratio->FindBin(mc_Q2) );
else return hQ2nonrelratio->GetBinContent( hQ2nonrelratio->FindBin(mc_Q2) );
// interpolation might be overkill for such a finely binned histogram, 0.01%
// but the extra cpu cycles maybe small.
// save it here for some future use.
//if(relORnonrel)return hQ2relratio->Interpolate(mc_Q2);
//else return hQ2nonrelratio->Interpolate(mc_Q2);
}
double weightRPA::getWeightLowQ2(const double mc_q0, const double mc_q3, int sign){
return getWeightInternal(mc_q0,mc_q3,weightRPA::LOWQ2,sign);
}
double weightRPA::getWeightHighQ2(const double mc_q0, const double mc_q3, int sign){
return getWeightInternal(mc_q0,mc_q3,weightRPA::HIGHQ2,sign);
}
#endif
diff --git a/src/Utils/FitUtils.cxx b/src/Utils/FitUtils.cxx
index f157e48..9e38b42 100644
--- a/src/Utils/FitUtils.cxx
+++ b/src/Utils/FitUtils.cxx
@@ -1,1056 +1,1086 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "FitUtils.h"
/*
MISC Functions
*/
//********************************************************************
double *FitUtils::GetArrayFromMap(std::vector invals,
std::map inmap) {
//********************************************************************
double *outarr = new double[invals.size()];
int count = 0;
for (size_t i = 0; i < invals.size(); i++) {
outarr[count++] = inmap[invals.at(i)];
}
return outarr;
}
/*
MISC Event
*/
//********************************************************************
// Returns the kinetic energy of a particle in GeV
double FitUtils::T(TLorentzVector part) {
//********************************************************************
double E_part = part.E() / 1000.;
double p_part = part.Vect().Mag() / 1000.;
double m_part = sqrt(E_part * E_part - p_part * p_part);
double KE_part = E_part - m_part;
return KE_part;
};
//********************************************************************
// Returns the momentum of a particle in GeV
double FitUtils::p(TLorentzVector part) {
//********************************************************************
double p_part = part.Vect().Mag() / 1000.;
return p_part;
};
double FitUtils::p(FitParticle *part) { return FitUtils::p(part->fP); };
//********************************************************************
// Returns the angle between two particles in radians
double FitUtils::th(TLorentzVector part1, TLorentzVector part2) {
//********************************************************************
double th = part1.Vect().Angle(part2.Vect());
return th;
};
double FitUtils::th(FitParticle *part1, FitParticle *part2) {
return FitUtils::th(part1->fP, part2->fP);
};
// T2K CC1pi+ helper functions
//
//********************************************************************
// Returns the angle between q3 and the pion defined in Raquel's CC1pi+ on CH
// paper
// Uses "MiniBooNE formula" for Enu, here called EnuCC1pip_T2K_MB
//********************************************************************
double FitUtils::thq3pi_CC1pip_T2K(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi) {
// Want this in GeV
TVector3 p_mu = pmu.Vect() * (1. / 1000.);
// Get the reconstructed Enu
// We are not using Michel e sample, so we have pion kinematic information
double Enu = EnuCC1piprec(pnu, pmu, ppi, true);
// Get neutrino unit direction, multiply by reconstructed Enu
TVector3 p_nu = pnu.Vect() * (1. / (pnu.Vect().Mag())) * Enu;
TVector3 p_pi = ppi.Vect() * (1. / 1000.);
// This is now in GeV
TVector3 q3 = (p_nu - p_mu);
// Want this in GeV
double th_q3_pi = q3.Angle(p_pi);
return th_q3_pi;
}
//********************************************************************
// Returns the q3 defined in Raquel's CC1pi+ on CH paper
// Uses "MiniBooNE formula" for Enu
//********************************************************************
double FitUtils::q3_CC1pip_T2K(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi) {
// Can't use the true Enu here; need to reconstruct it
// We do have Michel e- here so reconstruct Enu by "MiniBooNE formula" without
// pion kinematics
// The bool false refers to that we don't have pion kinematics
// Last bool refers to if we have pion kinematic information or not
double Enu = EnuCC1piprec(pnu, pmu, ppi, false);
TVector3 p_mu = pmu.Vect() * (1. / 1000.);
TVector3 p_nu = pnu.Vect() * (1. / pnu.Vect().Mag()) * Enu;
double q3 = (p_nu - p_mu).Mag();
return q3;
}
//********************************************************************
// Returns the W reconstruction from Raquel CC1pi+ CH thesis
// Uses the MiniBooNE formula Enu
//********************************************************************
double FitUtils::WrecCC1pip_T2K_MB(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi) {
double E_mu = pmu.E() / 1000.;
double p_mu = pmu.Vect().Mag() / 1000.;
double E_nu = EnuCC1piprec(pnu, pmu, ppi, false);
double a1 = (E_nu + PhysConst::mass_neutron) - E_mu;
double a2 = E_nu - p_mu;
double wrec = sqrt(a1 * a1 - a2 * a2);
return wrec;
}
//********************************************************
double FitUtils::ProtonQ2QErec(double pE, double binding) {
//********************************************************
const double V = binding / 1000.; // binding potential
const double mn = PhysConst::mass_neutron; // neutron mass
const double mp = PhysConst::mass_proton; // proton mass
const double mn_eff = mn - V; // effective proton mass
const double pki = (pE / 1000.0) - mp;
double q2qe = mn_eff * mn_eff - mp * mp + 2 * mn_eff * (pki + mp - mn_eff);
return q2qe;
};
//********************************************************************
double FitUtils::EnuQErec(TLorentzVector pmu, double costh, double binding,
bool neutrino) {
//********************************************************************
// Convert all values to GeV
const double V = binding / 1000.; // binding potential
const double mn = PhysConst::mass_neutron; // neutron mass
const double mp = PhysConst::mass_proton; // proton mass
double mN_eff = mn - V;
double mN_oth = mp;
if (!neutrino) {
mN_eff = mp - V;
mN_oth = mn;
}
double el = pmu.E() / 1000.;
double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton
double ml = sqrt(el * el - pl * pl); // lepton mass
double rEnu =
(2 * mN_eff * el - ml * ml + mN_oth * mN_oth - mN_eff * mN_eff) /
(2 * (mN_eff - el + pl * costh));
return rEnu;
};
double FitUtils::Q2QErec(TLorentzVector pmu, double costh, double binding, bool neutrino) {
double el = pmu.E() / 1000.;
double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton
double ml = sqrt(el * el - pl * pl); // lepton mass
double rEnu = EnuQErec(pmu, costh, binding, neutrino);
double q2 = -ml * ml + 2. * rEnu * (el - pl * costh);
return q2;
};
double FitUtils::Q2QErec(TLorentzVector Pmu, TLorentzVector Pnu, double binding, bool neutrino) {
double q2qe = Q2QErec(Pmu, cos(Pnu.Vect().Angle(Pmu.Vect())), binding, neutrino);
return q2qe;
}
double FitUtils::EnuQErec(double pl, double costh, double binding,
bool neutrino) {
if (pl < 0) return 0.; // Make sure nobody is silly
double mN_eff = PhysConst::mass_neutron - binding / 1000.;
double mN_oth = PhysConst::mass_proton;
if (!neutrino) {
mN_eff = PhysConst::mass_proton - binding / 1000.;
mN_oth = PhysConst::mass_neutron;
}
double ml = PhysConst::mass_muon;
double el = sqrt(pl * pl + ml * ml);
double rEnu =
(2 * mN_eff * el - ml * ml + mN_oth * mN_oth - mN_eff * mN_eff) /
(2 * (mN_eff - el + pl * costh));
return rEnu;
};
double FitUtils::Q2QErec(double pl, double costh, double binding,
bool neutrino) {
if (pl < 0) return 0.; // Make sure nobody is silly
double ml = PhysConst::mass_muon;
double el = sqrt(pl * pl + ml * ml);
double rEnu = EnuQErec(pl, costh, binding, neutrino);
double q2 = -ml * ml + 2. * rEnu * (el - pl * costh);
return q2;
};
//********************************************************************
// Reconstructs Enu for CC1pi0
// Very similar for CC1pi+ reconstruction
double FitUtils::EnuCC1pi0rec(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi0) {
//********************************************************************
double E_mu = pmu.E() / 1000;
double p_mu = pmu.Vect().Mag() / 1000;
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
double E_pi0 = ppi0.E() / 1000;
double p_pi0 = ppi0.Vect().Mag() / 1000;
double m_pi0 = sqrt(E_pi0 * E_pi0 - p_pi0 * p_pi0);
double th_nu_pi0 = pnu.Vect().Angle(ppi0.Vect());
const double m_n = PhysConst::mass_neutron; // neutron mass
const double m_p = PhysConst::mass_proton; // proton mass
double th_pi0_mu = ppi0.Vect().Angle(pmu.Vect());
double rEnu = (m_mu * m_mu + m_pi0 * m_pi0 + m_n * m_n - m_p * m_p -
2 * m_n * (E_pi0 + E_mu) + 2 * E_pi0 * E_mu -
2 * p_pi0 * p_mu * cos(th_pi0_mu)) /
(2 * (E_pi0 + E_mu - p_pi0 * cos(th_nu_pi0) -
p_mu * cos(th_nu_mu) - m_n));
return rEnu;
};
//********************************************************************
// Reconstruct Q2 for CC1pi0
// Beware: uses true Enu, not reconstructed Enu
double FitUtils::Q2CC1pi0rec(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi0) {
//********************************************************************
double E_mu = pmu.E() / 1000.; // energy of lepton in GeV
double p_mu = pmu.Vect().Mag() / 1000.; // momentum of lepton
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); // lepton mass
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
// double rEnu = EnuCC1pi0rec(pnu, pmu, ppi0); //reconstructed neutrino energy
// Use true neutrino energy
double rEnu = pnu.E() / 1000.;
double q2 = -m_mu * m_mu + 2. * rEnu * (E_mu - p_mu * cos(th_nu_mu));
return q2;
};
//********************************************************************
// Reconstruct Enu for CC1pi+
// pionInfo reflects if we're using pion kinematics or not
// In T2K CC1pi+ CH the Michel tag is used for pion in which pion kinematic info
// is lost and Enu is reconstructed without pion kinematics
double FitUtils::EnuCC1piprec(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi, bool pionInfo) {
//********************************************************************
double E_mu = pmu.E() / 1000.;
double p_mu = pmu.Vect().Mag() / 1000.;
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
double E_pi = ppi.E() / 1000.;
double p_pi = ppi.Vect().Mag() / 1000.;
double m_pi = sqrt(E_pi * E_pi - p_pi * p_pi);
const double m_n = PhysConst::mass_neutron; // neutron/proton mass
// should really take proton mass for proton interaction, neutron for neutron
// interaction. However, difference is pretty much negligable here!
// need this angle too
double th_nu_pi = pnu.Vect().Angle(ppi.Vect());
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
double th_pi_mu = ppi.Vect().Angle(pmu.Vect());
double rEnu = -999;
// If experiment doesn't have pion kinematic information (T2K CC1pi+ CH
// example of this; Michel e sample has no directional information on pion)
// We'll need to modify the reconstruction variables
if (!pionInfo) {
// Assumes that pion angle contribution contributes net zero
// Also assumes the momentum of reconstructed pions when using Michel
// electrons is 130 MeV
// So need to redefine E_pi and p_pi
// It's a little unclear what happens to the pmu . ppi term (4-vector); do
// we include the angular contribution?
// This below doesn't
th_nu_pi = M_PI / 2.;
p_pi = 0.130;
E_pi = sqrt(p_pi * p_pi + m_pi * m_pi);
// rEnu = (m_mu*m_mu + m_pi*m_pi - 2*m_n*(E_mu + E_pi) + 2*E_mu*E_pi -
// 2*p_mu*p_pi) / (2*(E_mu + E_pi - p_mu*cos(th_nu_mu) - m_n));
}
rEnu =
(m_mu * m_mu + m_pi * m_pi - 2 * m_n * (E_pi + E_mu) + 2 * E_pi * E_mu -
2 * p_pi * p_mu * cos(th_pi_mu)) /
(2 * (E_pi + E_mu - p_pi * cos(th_nu_pi) - p_mu * cos(th_nu_mu) - m_n));
return rEnu;
};
//********************************************************************
// Reconstruct neutrino energy from outgoing particles; will differ from the
// actual neutrino energy. Here we use assumption of a Delta resonance
double FitUtils::EnuCC1piprecDelta(TLorentzVector pnu, TLorentzVector pmu) {
//********************************************************************
const double m_Delta =
PhysConst::mass_delta; // PDG value for Delta mass in GeV
const double m_n = PhysConst::mass_neutron; // neutron/proton mass
// should really take proton mass for proton interaction, neutron for neutron
// interaction. However, difference is pretty much negligable here!
double E_mu = pmu.E() / 1000;
double p_mu = pmu.Vect().Mag() / 1000;
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
double rEnu = (m_Delta * m_Delta - m_n * m_n - m_mu * m_mu + 2 * m_n * E_mu) /
(2 * (m_n - E_mu + p_mu * cos(th_nu_mu)));
return rEnu;
};
// MOVE TO T2K UTILS!
//********************************************************************
// Reconstruct Enu using "extended MiniBooNE" as defined in Raquel's T2K TN
//
// Supposedly includes pion direction and binding energy of target nucleon
// I'm not convinced (yet), maybe
double FitUtils::EnuCC1piprec_T2K_eMB(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi) {
//********************************************************************
// Unit vector for neutrino momentum
TVector3 p_nu_vect_unit = pnu.Vect() * (1. / pnu.E());
double E_mu = pmu.E() / 1000.;
TVector3 p_mu_vect = pmu.Vect() * (1. / 1000.);
double E_pi = ppi.E() / 1000.;
TVector3 p_pi_vect = ppi.Vect() * (1. / 1000.);
double E_bind =
27. / 1000.; // This should be roughly correct for CH; but not clear!
double m_p = PhysConst::mass_proton;
// Makes life a little easier, gonna square this one
double a1 = m_p - E_bind - E_mu - E_pi;
// Just gets the magnitude square of the muon and pion momentum vectors
double a2 = (p_mu_vect + p_pi_vect).Mag2();
// Gets the somewhat complicated scalar product between neutrino and
// (p_mu+p_pi), scaled to Enu
double a3 = p_nu_vect_unit * (p_mu_vect + p_pi_vect);
double rEnu =
(m_p * m_p - a1 * a1 + a2) / (2 * (m_p - E_bind - E_mu - E_pi + a3));
return rEnu;
}
//********************************************************************
// Reconstructed Q2 for CC1pi+
//
// enuType describes how the neutrino energy is reconstructed
// 0 uses true neutrino energy; case for MINERvA and MiniBooNE
// 1 uses "extended MiniBooNE" reconstructed (T2K CH)
// 2 uses "MiniBooNE" reconstructed (EnuCC1piprec with pionInfo = true) (T2K CH)
// "MiniBooNE" reconstructed (EnuCC1piprec with pionInfo = false, the
// case for T2K when using Michel tag) (T2K CH)
// 3 uses Delta for reconstruction (T2K CH)
double FitUtils::Q2CC1piprec(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi, int enuType, bool pionInfo) {
//********************************************************************
double E_mu = pmu.E() / 1000.; // energy of lepton in GeV
double p_mu = pmu.Vect().Mag() / 1000.; // momentum of lepton
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); // lepton mass
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
// Different ways of reconstructing the neutrino energy; default is just to
// use the truth (case 0)
double rEnu = -999;
switch (enuType) {
case 0: // True neutrino energy, default; this is the case for when Q2
// defined as Q2 true (MiniBooNE, MINERvA)
rEnu = pnu.E() / 1000.;
break;
case 1: // Extended MiniBooNE reconstructed, as defined by Raquel's CC1pi+
// CH T2K analysis
// Definitely uses pion info :)
rEnu = EnuCC1piprec_T2K_eMB(pnu, pmu, ppi);
break;
case 2: // MiniBooNE reconstructed, as defined by MiniBooNE and Raquel's
// CC1pi+ CH T2K analysis and Linda's CC1pi+ H2O
// Can have this with and without pion info, depending on if Michel tag
// used (Raquel)
rEnu = EnuCC1piprec(pnu, pmu, ppi, pionInfo);
break;
case 3:
rEnu = EnuCC1piprecDelta(pnu, pmu);
break;
} // No need for default here since default value of enuType = 0, defined in
// header
double q2 = -m_mu * m_mu + 2. * rEnu * (E_mu - p_mu * cos(th_nu_mu));
return q2;
};
//********************************************************************
// Returns the reconstructed W from a nucleon and an outgoing pion
//
// Could do this a lot more clever (pp + ppi).Mag() would do the job, but this
// would be less instructive
//********************************************************************
double FitUtils::MpPi(TLorentzVector pp, TLorentzVector ppi) {
double E_p = pp.E();
double p_p = pp.Vect().Mag();
double m_p = sqrt(E_p * E_p - p_p * p_p);
double E_pi = ppi.E();
double p_pi = ppi.Vect().Mag();
double m_pi = sqrt(E_pi * E_pi - p_pi * p_pi);
double th_p_pi = pp.Vect().Angle(ppi.Vect());
// fairly easy thing to derive since bubble chambers measure the proton!
double invMass = sqrt(m_p * m_p + m_pi * m_pi + 2 * E_p * E_pi -
2 * p_pi * p_p * cos(th_p_pi));
return invMass;
};
//********************************************************
// Reconstruct the hadronic mass using neutrino and muon
// Could technically do E_nu = EnuCC1pipRec(pnu,pmu,ppi) too, but this wwill be
// reconstructed Enu; so gives reconstructed Wrec which most of the time isn't
// used!
//
// Only MINERvA uses this so far; and the Enu is Enu_true
// If we want W_true need to take initial state nucleon motion into account
// Return value is in MeV!!!
double FitUtils::Wrec(TLorentzVector pnu, TLorentzVector pmu) {
//********************************************************
double E_mu = pmu.E();
double p_mu = pmu.Vect().Mag();
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
// The factor of 1000 is necessary for downstream functions
const double m_p = PhysConst::mass_proton * 1000;
// MINERvA cut on W_exp which is tuned to W_true; so use true Enu from
// generators
double E_nu = pnu.E();
double w_rec = sqrt(m_p * m_p + m_mu * m_mu - 2 * m_p * E_mu +
2 * E_nu * (m_p - E_mu + p_mu * cos(th_nu_mu)));
return w_rec;
};
//********************************************************
// Reconstruct the true hadronic mass using the initial state and muon
// Could technically do E_nu = EnuCC1pipRec(pnu,pmu,ppi) too, but this wwill be
// reconstructed Enu; so gives reconstructed Wrec which most of the time isn't
// used!
//
// No one seems to use this because it's fairly MC dependent!
// Return value is in MeV!!!
double FitUtils::Wtrue(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector pnuc) {
//********************************************************
// Could simply do the TLorentzVector operators here but this is more
// instructive?
// ... and prone to errors ...
double E_mu = pmu.E();
double p_mu = pmu.Vect().Mag();
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
double E_nuc = pnuc.E();
double p_nuc = pnuc.Vect().Mag();
double m_nuc = sqrt(E_nuc * E_nuc - p_nuc * p_nuc);
double th_nuc_mu = pmu.Vect().Angle(pnuc.Vect());
double th_nu_nuc = pnu.Vect().Angle(pnuc.Vect());
double E_nu = pnu.E();
double w_rec = sqrt(m_nuc * m_nuc + m_mu * m_mu - 2 * E_nu * E_mu +
2 * E_nu * p_mu * cos(th_nu_mu) - 2 * E_nuc * E_mu +
2 * p_nuc * p_mu * cos(th_nuc_mu) + 2 * E_nu * E_nuc -
2 * E_nu * p_nuc * cos(th_nu_nuc));
return w_rec;
};
double FitUtils::SumKE_PartVect(std::vector const fps) {
double sum = 0.0;
for (size_t p_it = 0; p_it < fps.size(); ++p_it) {
sum += fps[p_it]->KE();
}
return sum;
}
double FitUtils::SumTE_PartVect(std::vector const fps) {
double sum = 0.0;
for (size_t p_it = 0; p_it < fps.size(); ++p_it) {
sum += fps[p_it]->E();
}
return sum;
}
/*
E Recoil
*/
double FitUtils::GetErecoil_TRUE(FitEvent *event) {
// Get total energy of hadronic system.
double Erecoil = 0.0;
for (unsigned int i = 2; i < event->Npart(); i++) {
// Only final state
if (!event->PartInfo(i)->fIsAlive) continue;
if (event->PartInfo(i)->fNEUTStatusCode != 0) continue;
// Skip Lepton
if (abs(event->PartInfo(i)->fPID) == abs(event->PartInfo(0)->fPID) - 1)
continue;
// Add Up KE of protons and TE of everything else
if (event->PartInfo(i)->fPID == 2212 || event->PartInfo(i)->fPID == 2112) {
Erecoil +=
fabs(event->PartInfo(i)->fP.E()) - fabs(event->PartInfo(i)->fP.Mag());
} else {
Erecoil += event->PartInfo(i)->fP.E();
}
}
return Erecoil;
}
double FitUtils::GetErecoil_CHARGED(FitEvent *event) {
// Get total energy of hadronic system.
double Erecoil = 0.0;
for (unsigned int i = 2; i < event->Npart(); i++) {
// Only final state
if (!event->PartInfo(i)->fIsAlive) continue;
if (event->PartInfo(i)->fNEUTStatusCode != 0) continue;
// Skip Lepton
if (abs(event->PartInfo(i)->fPID) == abs(event->PartInfo(0)->fPID) - 1)
continue;
// Skip Neutral particles
if (event->PartInfo(i)->fPID == 2112 || event->PartInfo(i)->fPID == 111 ||
event->PartInfo(i)->fPID == 22)
continue;
// Add Up KE of protons and TE of everything else
if (event->PartInfo(i)->fPID == 2212) {
Erecoil +=
fabs(event->PartInfo(i)->fP.E()) - fabs(event->PartInfo(i)->fP.Mag());
} else {
Erecoil += event->PartInfo(i)->fP.E();
}
}
return Erecoil;
}
// MOVE TO MINERVA Utils!
double FitUtils::GetErecoil_MINERvA_LowRecoil(FitEvent *event) {
// Get total energy of hadronic system.
double Erecoil = 0.0;
for (unsigned int i = 2; i < event->Npart(); i++) {
// Only final state
if (!event->PartInfo(i)->fIsAlive) continue;
if (event->PartInfo(i)->fNEUTStatusCode != 0) continue;
// Skip Lepton
if (abs(event->PartInfo(i)->fPID) == 13) continue;
// Skip Neutrons particles
if (event->PartInfo(i)->fPID == 2112) continue;
int PID = event->PartInfo(i)->fPID;
// KE of Protons and charged pions
if (PID == 2212 or PID == 211 or PID == -211) {
// Erecoil += FitUtils::T(event->PartInfo(i)->fP);
Erecoil +=
fabs(event->PartInfo(i)->fP.E()) - fabs(event->PartInfo(i)->fP.Mag());
// Total Energy of non-neutrons
// } else if (PID != 2112 and PID < 999 and PID != 22 and abs(PID) !=
// 14) {
} else if (PID == 111 || PID == 11 || PID == -11 || PID == 22) {
Erecoil += (event->PartInfo(i)->fP.E());
}
}
return Erecoil;
}
+// MOVE TO MINERVA Utils!
+// The alternative Eavailble definition takes true q0 and subtracts the kinetic energy of neutrons and pion masses
+// returns in MeV
+double FitUtils::Eavailable(FitEvent *event) {
+ double Eav = 0.0;
+
+ // Now take q0 and subtract Eav
+ double q0 = event->GetNeutrinoIn()->fP.E();
+ if (event->GetHMFSMuon()) q0 -= event->GetHMFSMuon()->fP.E();
+ else if (!event->GetHMFSNuMuon()) q0 -= event->GetHMFSNuMuon()->fP.E();
+
+ for (unsigned int i = 2; i < event->Npart(); i++) {
+ // Only final state
+ if (!event->PartInfo(i)->fIsAlive) continue;
+ if (event->PartInfo(i)->fNEUTStatusCode != 0) continue;
+ int PID = event->PartInfo(i)->fPID;
+
+ // Neutrons
+ if (PID == 2112) {
+ // Adding kinetic energy of neutron
+ Eav += FitUtils::T(event->PartInfo(i)->fP)*1000.;
+ // All pion masses
+ } else if (abs(PID) == 211 || PID == 111) {
+ Eav += event->PartInfo(i)->fP.M();
+ }
+ }
+
+ return q0-Eav;
+}
+
TVector3 GetVectorInTPlane(const TVector3 &inp, const TVector3 &planarNormal) {
TVector3 pnUnit = planarNormal.Unit();
double inpProjectPN = inp.Dot(pnUnit);
TVector3 InPlanarInput = inp - (inpProjectPN * pnUnit);
return InPlanarInput;
}
TVector3 GetUnitVectorInTPlane(const TVector3 &inp,
const TVector3 &planarNormal) {
return GetVectorInTPlane(inp, planarNormal).Unit();
}
Double_t GetDeltaPhiT(TVector3 const &V_lepton, TVector3 const &V_other,
TVector3 const &Normal, bool PiMinus = false) {
TVector3 V_lepton_T = GetUnitVectorInTPlane(V_lepton, Normal);
TVector3 V_other_T = GetUnitVectorInTPlane(V_other, Normal);
return PiMinus ? acos(V_lepton_T.Dot(V_other_T))
: (M_PI - acos(V_lepton_T.Dot(V_other_T)));
}
TVector3 GetDeltaPT(TVector3 const &V_lepton, TVector3 const &V_other,
TVector3 const &Normal) {
TVector3 V_lepton_T = GetVectorInTPlane(V_lepton, Normal);
TVector3 V_other_T = GetVectorInTPlane(V_other, Normal);
return V_lepton_T + V_other_T;
}
Double_t GetDeltaAlphaT(TVector3 const &V_lepton, TVector3 const &V_other,
TVector3 const &Normal, bool PiMinus = false) {
TVector3 DeltaPT = GetDeltaPT(V_lepton, V_other, Normal);
return GetDeltaPhiT(V_lepton, DeltaPT, Normal, PiMinus);
}
double FitUtils::Get_STV_dpt(FitEvent *event, int ISPDG, bool Is0pi) {
// Check that the neutrino exists
if (event->NumISParticle(ISPDG) == 0) {
return -9999;
}
// Return 0 if the proton or muon are missing
if (event->NumFSParticle(2212) == 0 ||
event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) {
return -9999;
}
// Now get the TVector3s for each particle
TVector3 const &NuP = event->GetHMISParticle(14)->fP.Vect();
TVector3 const &LeptonP =
event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect();
// Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
int HMFSProton = 0;
double HighestMomentum = 0.0;
// Get the stack of protons
std::vector Protons = event->GetAllFSProton();
for (size_t i = 0; i < Protons.size(); ++i) {
if (Protons[i]->p() > 450 &&
Protons[i]->p() < 1200 &&
Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 &&
Protons[i]->p() > HighestMomentum) {
HighestMomentum = Protons[i]->p();
HMFSProton = i;
}
}
// Now get the proton
TLorentzVector Pprot = Protons[HMFSProton]->fP;
// Get highest momentum proton in allowed proton range
TVector3 HadronP = Pprot.Vect();
// If we don't have a CC0pi signal definition we also add in pion momentum
if (!Is0pi) {
if (event->NumFSParticle(PhysConst::pdg_pions) == 0) {
return -9999;
}
// Count up pion momentum
TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP;
HadronP += pp.Vect();
}
return GetDeltaPT(LeptonP, HadronP, NuP).Mag();
}
double FitUtils::Get_STV_dphit(FitEvent *event, int ISPDG, bool Is0pi) {
// Check that the neutrino exists
if (event->NumISParticle(ISPDG) == 0) {
return -9999;
}
// Return 0 if the proton or muon are missing
if (event->NumFSParticle(2212) == 0 ||
event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) {
return -9999;
}
// Now get the TVector3s for each particle
TVector3 const &NuP = event->GetHMISParticle(ISPDG)->fP.Vect();
TVector3 const &LeptonP =
event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect();
// Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
int HMFSProton = 0;
double HighestMomentum = 0.0;
// Get the stack of protons
std::vector Protons = event->GetAllFSProton();
for (size_t i = 0; i < Protons.size(); ++i) {
if (Protons[i]->p() > 450 &&
Protons[i]->p() < 1200 &&
Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 &&
Protons[i]->p() > HighestMomentum) {
HighestMomentum = Protons[i]->p();
HMFSProton = i;
}
}
// Now get the proton
TLorentzVector Pprot = Protons[HMFSProton]->fP;
// Get highest momentum proton in allowed proton range
TVector3 HadronP = Pprot.Vect();
if (!Is0pi) {
if (event->NumFSParticle(PhysConst::pdg_pions) == 0) {
return -9999;
}
TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP;
HadronP += pp.Vect();
}
return GetDeltaPhiT(LeptonP, HadronP, NuP);
}
double FitUtils::Get_STV_dalphat(FitEvent *event, int ISPDG, bool Is0pi) {
// Check that the neutrino exists
if (event->NumISParticle(ISPDG) == 0) {
return -9999;
}
// Return 0 if the proton or muon are missing
if (event->NumFSParticle(2212) == 0 ||
event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) {
return -9999;
}
// Now get the TVector3s for each particle
TVector3 const &NuP = event->GetHMISParticle(ISPDG)->fP.Vect();
TVector3 const &LeptonP =
event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect();
// Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
int HMFSProton = 0;
double HighestMomentum = 0.0;
// Get the stack of protons
std::vector Protons = event->GetAllFSProton();
for (size_t i = 0; i < Protons.size(); ++i) {
if (Protons[i]->p() > 450 &&
Protons[i]->p() < 1200 &&
Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 &&
Protons[i]->p() > HighestMomentum) {
HighestMomentum = Protons[i]->p();
HMFSProton = i;
}
}
// Now get the proton
TLorentzVector Pprot = Protons[HMFSProton]->fP;
// Get highest momentum proton in allowed proton range
TVector3 HadronP = Pprot.Vect();
if (!Is0pi) {
if (event->NumFSParticle(PhysConst::pdg_pions) == 0) {
return -9999;
}
TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP;
HadronP += pp.Vect();
}
return GetDeltaAlphaT(LeptonP, HadronP, NuP);
}
// As defined in PhysRevC.95.065501
// Using prescription from arXiv 1805.05486
// Returns in GeV
double FitUtils::Get_pn_reco_C(FitEvent *event, int ISPDG, bool Is0pi) {
const double mn = PhysConst::mass_neutron; // neutron mass
const double mp = PhysConst::mass_proton; // proton mass
// Check that the neutrino exists
if (event->NumISParticle(ISPDG) == 0) {
return -9999;
}
// Return 0 if the proton or muon are missing
if (event->NumFSParticle(2212) == 0 ||
event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) {
return -9999;
}
// Now get the TVector3s for each particle
TVector3 const &NuP = event->GetHMISParticle(14)->fP.Vect();
TVector3 const &LeptonP =
event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect();
// Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
int HMFSProton = 0;
double HighestMomentum = 0.0;
// Get the stack of protons
std::vector Protons = event->GetAllFSProton();
for (size_t i = 0; i < Protons.size(); ++i) {
// Update the highest momentum particle
if (Protons[i]->p() > 450 &&
Protons[i]->p() < 1200 &&
Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 &&
Protons[i]->p() > HighestMomentum) {
HighestMomentum = Protons[i]->p();
HMFSProton = i;
}
}
// Now get the proton
TLorentzVector Pprot = Protons[HMFSProton]->fP;
// Get highest momentum proton in allowed proton range
TVector3 HadronP = Pprot.Vect();
//TVector3 HadronP = event->GetHMFSParticle(2212)->fP.Vect();
double const el = event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->E()/1000.;
double const eh = Pprot.E()/1000.;
if (!Is0pi) {
if (event->NumFSParticle(PhysConst::pdg_pions) == 0) {
return -9999;
}
TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP;
HadronP += pp.Vect();
}
TVector3 dpt = GetDeltaPT(LeptonP, HadronP, NuP);
double dptMag = dpt.Mag()/1000.;
double ma = 6*mn + 6*mp - 0.09216; // target mass (E is from PhysRevC.95.065501)
double map = ma - mn + 0.02713; // remnant mass
double pmul = LeptonP.Dot(NuP.Unit())/1000.;
double phl = HadronP.Dot(NuP.Unit())/1000.;
//double pmul = GetVectorInTPlane(LeptonP, dpt).Mag()/1000.;
//double phl = GetVectorInTPlane(HadronP, dpt).Mag()/1000.;
double R = ma + pmul + phl - el - eh;
double dpl = 0.5*R - (map*map + dptMag*dptMag)/(2*R);
//double dpl = ((R*R)-(dptMag*dptMag)-(map*map))/(2*R); // as in in PhysRevC.95.065501 - gives same result
double pn_reco = sqrt((dptMag*dptMag) + (dpl*dpl));
//std::cout << "Diagnostics: " << std::endl;
//std::cout << "mn: " << mn << std::endl;
//std::cout << "ma: " << ma << std::endl;
//std::cout << "map: " << map << std::endl;
//std::cout << "pmu: " << LeptonP.Mag()/1000. << std::endl;
//std::cout << "ph: " << HadronP.Mag()/1000. << std::endl;
//std::cout << "pmul: " << pmul << std::endl;
//std::cout << "phl: " << phl << std::endl;
//std::cout << "el: " << el << std::endl;
//std::cout << "eh: " << eh << std::endl;
//std::cout << "R: " << R << std::endl;
//std::cout << "dptMag: " << dptMag << std::endl;
//std::cout << "dpl: " << dpl << std::endl;
//std::cout << "pn_reco: " << pn_reco << std::endl;
return pn_reco;
}
// Get Cos theta with Adler angles
double FitUtils::CosThAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi, TLorentzVector Pprot) {
// Get the "resonance" lorentz vector (pion proton system)
TLorentzVector Pres = Pprot + Ppi;
// Boost the particles into the resonance rest frame so we can define the x,y,z axis
Pnu.Boost(Pres.BoostVector());
Pmu.Boost(-Pres.BoostVector());
Ppi.Boost(-Pres.BoostVector());
// The z-axis is defined as the axis of three-momentum transfer, \vec{k}
// Also unit normalise the axis
TVector3 zAxis = (Pnu.Vect()-Pmu.Vect())*(1.0/((Pnu.Vect()-Pmu.Vect()).Mag()));
// Then the angle between the pion in the "resonance" rest-frame and the z-axis is the theta Adler angle
double costhAdler = cos(Ppi.Vect().Angle(zAxis));
return costhAdler;
}
// Get phi with Adler angles, a bit more complicated...
double FitUtils::PhiAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi, TLorentzVector Pprot) {
// Get the "resonance" lorentz vector (pion proton system)
TLorentzVector Pres = Pprot + Ppi;
// Boost the particles into the resonance rest frame so we can define the x,y,z axis
Pnu.Boost(Pres.BoostVector());
Pmu.Boost(-Pres.BoostVector());
Ppi.Boost(-Pres.BoostVector());
// The z-axis is defined as the axis of three-momentum transfer, \vec{k}
// Also unit normalise the axis
TVector3 zAxis = (Pnu.Vect()-Pmu.Vect())*(1.0/((Pnu.Vect()-Pmu.Vect()).Mag()));
// The y-axis is then defined perpendicular to z and muon momentum in the resonance frame
TVector3 yAxis = Pnu.Vect().Cross(Pmu.Vect());
yAxis *= 1.0/double(yAxis.Mag());
// And the x-axis is then simply perpendicular to z and x
TVector3 xAxis = yAxis.Cross(zAxis);
xAxis *= 1.0/double(xAxis.Mag());
double x = Ppi.Vect().Dot(xAxis);
double y = Ppi.Vect().Dot(yAxis);
//double z = Ppi.Vect().Dot(zAxis);
double newphi = atan2(y, x)*(180./M_PI);
// Convert negative angles to positive
if (newphi < 0.0) newphi += 360.0;
// Old silly method before atan2
/*
// Then finally construct phi as the angle between pion projection and x axis
// Get the project of the pion momentum on to the zaxis
TVector3 PiVectZ = zAxis*Ppi.Vect().Dot(zAxis);
// The subtract the projection off the pion vector to get to get the plane
TVector3 PiPlane = Ppi.Vect() - PiVectZ;
double phi = -999.99;
if (PiPlane.Y() > 0) {
phi = (180./M_PI)*PiPlane.Angle(xAxis);
} else if (PiPlane.Y() < 0) {
phi = (180./M_PI)*(2*M_PI-PiPlane.Angle(xAxis));
} else if (PiPlane.Y() == 0) {
TRandom3 rand;
double randNo = rand.Rndm();
if (randNo > 0.5) {
phi = (180./M_PI)*PiPlane.Angle(xAxis);
} else {
phi = (180./M_PI)*(2*M_PI-PiPlane.Angle(xAxis));
}
}
*/
return newphi;
}
//********************************************************************
double FitUtils::ppInfK(TLorentzVector pmu, double costh, double binding,
bool neutrino) {
//********************************************************************
// Convert all values to GeV
//const double V = binding / 1000.; // binding potential
//const double mn = PhysConst::mass_neutron; // neutron mass
const double mp = PhysConst::mass_proton; // proton mass
double el = pmu.E() / 1000.;
//double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton
double enu = EnuQErec(pmu, costh, binding, neutrino);
double ep_inf = enu - el + mp;
double pp_inf = sqrt(ep_inf * ep_inf - mp * mp);
return pp_inf;
};
//********************************************************************
TVector3 FitUtils::tppInfK(TLorentzVector pmu, double costh, double binding,
bool neutrino) {
//********************************************************************
// Convert all values to GeV
//const double V = binding / 1000.; // binding potential
//const double mn = PhysConst::mass_neutron; // neutron mass
//const double mp = PhysConst::mass_proton; // proton mass
double pl_x = pmu.X() / 1000.;
double pl_y = pmu.Y() / 1000.;
double pl_z= pmu.Z() / 1000.;
double enu = EnuQErec(pmu, costh, binding, neutrino);
TVector3 tpp_inf(-pl_x, -pl_y, -pl_z+enu);
return tpp_inf;
};
//********************************************************************
double FitUtils::cthpInfK(TLorentzVector pmu, double costh, double binding,
bool neutrino) {
//********************************************************************
// Convert all values to GeV
//const double V = binding / 1000.; // binding potential
//const double mn = PhysConst::mass_neutron; // neutron mass
const double mp = PhysConst::mass_proton; // proton mass
double el = pmu.E() / 1000.;
double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton
double enu = EnuQErec(pmu, costh, binding, neutrino);
double ep_inf = enu - el + mp;
double pp_inf = sqrt(ep_inf * ep_inf - mp * mp);
double cth_inf = (enu*enu + pp_inf*pp_inf - pl*pl)/(2*enu*pp_inf);
return cth_inf;
};
diff --git a/src/Utils/FitUtils.h b/src/Utils/FitUtils.h
index 11e2d5f..9510e20 100644
--- a/src/Utils/FitUtils.h
+++ b/src/Utils/FitUtils.h
@@ -1,185 +1,186 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifndef FITUTILS_H_SEEN
#define FITUTILS_H_SEEN
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include "FitEvent.h"
#include "TGraph.h"
#include "TH2Poly.h"
#include "FitEvent.h"
#include "FitLogger.h"
/*!
* \addtogroup Utils
* @{
*/
/// Functions needed by individual samples for calculating kinematic quantities.
namespace FitUtils {
/// Return a vector of all values saved in map
double *GetArrayFromMap(std::vector invals,
std::map inmap);
/// Returns kinetic energy of particle
double T(TLorentzVector part);
/// Returns momentum of particle
double p(TLorentzVector part);
double p(FitParticle* part);
/// Returns angle between particles (_NOT_ cosine!)
double th(TLorentzVector part, TLorentzVector part2);
double th(FitParticle* part1, FitParticle* part2);
/// Hadronic mass reconstruction
double Wrec(TLorentzVector pnu, TLorentzVector pmu);
/// Hadronic mass true from initial state particles and muon; useful if the full
/// FSI vectors aren't not saved and we for some reasons need W_true
double Wtrue(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector pnuc);
double SumKE_PartVect(std::vector const fps);
double SumTE_PartVect(std::vector const fps);
/// Return E Hadronic for all FS Particles in Hadronic System
double GetErecoil_TRUE(FitEvent *event);
/// Return E Hadronic for all Charged FS Particles in Hadronic System
double GetErecoil_CHARGED(FitEvent *event);
+double Eavailable(FitEvent *event);
/*
CCQE MiniBooNE/MINERvA
*/
/// Function to calculate the reconstructed Q^{2}_{QE}
double Q2QErec(TLorentzVector pmu, double costh, double binding, bool neutrino = true);
/// Function returns the reconstructed E_{nu} values
double EnuQErec(TLorentzVector pmu, double costh, double binding, bool neutrino = true);
//! Function to calculate the reconstructed Q^{2}_{QE}
double Q2QErec(double pl, double costh, double binding, bool neutrino = true);
//! Function to calculate the reconstructed Q^{2}_{QE}
double Q2QErec(TLorentzVector Pmu, TLorentzVector Pnu, double binding, bool neutrino = true);
//! Function returns the reconstructed E_{nu} values
double EnuQErec(double pl, double costh, double binding,
bool neutrino = true);
/*
CCQE1p MINERvA
*/
/// Reconstruct Q2QE given just the maximum energy proton.
double ProtonQ2QErec(double pE, double binding);
/*
E Recoil MINERvA
*/
double GetErecoil_MINERvA_LowRecoil(FitEvent *event);
/*
CC1pi0 MiniBooNE
*/
/// Reconstruct Enu from CCpi0 vectors and binding energy
double EnuCC1pi0rec(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi0 = TLorentzVector(0, 0, 0, 0));
/// Reconstruct Q2 from CCpi0 vectors and binding energy
double Q2CC1pi0rec(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi0 = TLorentzVector(0, 0, 0, 0));
/*
CC1pi+ MiniBooNE
*/
/// returns reconstructed Enu a la MiniBooNE CCpi+
/// returns reconstructed Enu a la MiniBooNE CCpi+
// Also for when not having pion info (so when we have a Michel tag in T2K)
double EnuCC1piprec(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppip,
bool pionInfo = true);
/// returns reconstructed Enu assumming resonance interaction where intermediate
/// resonance was a Delta
double EnuCC1piprecDelta(TLorentzVector pnu, TLorentzVector pmu);
/// returns reconstructed in a variety of flavours
double Q2CC1piprec(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppip,
int enuType = 0, bool pionInfo = true);
/*
T2K CC1pi+ on CH
*/
double thq3pi_CC1pip_T2K(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi);
double q3_CC1pip_T2K(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi);
double WrecCC1pip_T2K_MB(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppip);
double EnuCC1piprec_T2K_eMB(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi);
/*
nucleon single pion
*/
double MpPi(TLorentzVector pp, TLorentzVector ppi);
/// Gets delta p T as defined in Phys.Rev. C94 (2016) no.1, 015503
double Get_STV_dpt(FitEvent *event, int ISPDG, bool Is0pi);
/// Gets delta phi T as defined in Phys.Rev. C94 (2016) no.1, 015503
double Get_STV_dphit(FitEvent *event, int ISPDG, bool Is0pi);
/// Gets delta alpha T as defined in Phys.Rev. C94 (2016) no.1, 015503
double Get_STV_dalphat(FitEvent *event, int ISPDG, bool Is0pi);
// As defined in PhysRevC.95.065501
// Using prescription from arXiv 1805.05486
double Get_pn_reco_C(FitEvent *event, int ISPDG, bool Is0pi);
//For T2K inferred kinematics analyis - variables defined as on page 7 of T2K TN287v11 (and now arXiv 1802.05078)
double ppInfK(TLorentzVector pmu, double costh, double binding, bool neutrino);
TVector3 tppInfK(TLorentzVector pmu, double costh, double binding, bool neutrino);
double cthpInfK(TLorentzVector pmu, double costh, double binding, bool neutrino);
double CosThAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi, TLorentzVector Pprot);
double PhiAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi, TLorentzVector Pprot);
}
/*! @} */
#endif