Page MenuHomeHEPForge
Files F9514778

No OneTemporary
Actions

View Options

diff --git a/data/ArgoNeuT/CCInc_dsig_dmumom_nu.dat b/data/ArgoNeuT/CCInc_dsig_dmumom_nu.dat
index e24fa9c..9d7f7fa 100644
--- a/data/ArgoNeuT/CCInc_dsig_dmumom_nu.dat
+++ b/data/ArgoNeuT/CCInc_dsig_dmumom_nu.dat
@@ -1,22 +1,22 @@
-# Bin (GeV/c), dsig/dp_mu (10^-38 cm^2 / GeV/c), Error
+# Bin (GeV/c), dsig/dp_mu (10^-38 cm^2 /Ar /GeV /c), Error
# http://arxiv.org/pdf/1404.4809v3.pdf
0 20.3 7.3
1.25 27.1 4.7
2.5 18.9 3.7
3.75 22.5 3.9
5 17.3 3.1
6.25 15.1 2.9
7.5 10.8 2.1
8.75 12.0 3.0
10 7.64 1.74
11.2 7.43 1.66
12.5 6.36 1.75
13.8 3.80 1.10
15 3.75 1.06
16.2 2.84 0.89
17.5 2.01 0.88
18.8 3.96 1.19
20 2.83 1.11
21.2 2.65 0.96
22.5 1.63 0.73
23.8 1.40 0.67
diff --git a/data/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.dat b/data/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.dat
new file mode 100644
index 0000000..942baa5
--- /dev/null
+++ b/data/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.dat
@@ -0,0 +1,13 @@
+# Bin (GeV/c), dsig/ddpt (10^-38 cm^2 /GeV /neutron /c), Error
+# PRELIMINARY Fake NuWro data in proposed binning
+# Binning taken from: https://arxiv.org/pdf/1605.00179v1.pdf Fig.3
+0.0 0.746 0.1
+0.08 1.09 0.35
+0.12 0.975 0.35
+0.155 0.754 0.3
+0.2 0.396 0.125
+0.26 0.161 0.1
+0.36 0.114 0.075
+0.51 0.021 0.005
+1.1 0 0
+
diff --git a/event_gen/nuwro/ND280/ND280_CH_FHC.params b/event_gen/nuwro/ND280/ND280_CH_FHC.params
new file mode 100644
index 0000000..22b3124
--- /dev/null
+++ b/event_gen/nuwro/ND280/ND280_CH_FHC.params
@@ -0,0 +1,136 @@
+random_seed = 0
+number_of_events = 50000
+number_of_test_events = 1000000
+save_test_events = 0
+OutputEvtHistogram = 1
+user_events = 0
+user_params =
+beam_type = 1
+beam_energy =1000
+beam_particle = 14
+beam_direction = 0 0 1
+beam_content = 14 92.8% 50 9450 38571.4 82467.5 134286 208016 291786 389464 496786 601825 747946 945238 1.09036e+06 1.13119e+06 1.07349e+06 915000 676786 455000 293571 194048 141464 110714 91904.8 78285.7 68015.9 60918.4 53954.1 48452.4 43035.7 39365.1 35754 33992.3 32251.3 30510.2 28769.1 27028.1 25287 24085.1 22989.6 21894.1 20798.6 19703 18607.5 17512 16416.5 15321 14266.5 13917.1 13567.6 13218.1 12868.7 12519.2 12169.7 11820.3 11470.8 11121.4 10771.9 10507.9 10260.8 10013.6 9766.52 9519.39 9272.26 9025.14 8778.01 8530.88 8308.78 8206.78 8104.78 8002.77 7900.77 7798.77 7696.77 7594.77 7492.77 7390.76 7288.76 7196.35 7105.94 7015.53 6925.12 6834.71 6744.3 6653.89 6563.47 6473.06 6382.65 6245.64 6098.92 5952.2 5805.48 5658.75 5512.03 5365.31 5218.59 5071.87 4925.15 4787.18 4651.05 4514.91 4378.77 4242.63 4106.49 3970.35 3834.22 3698.08 3561.94 3425.8 3289.66 3153.53 3017.39 2881.25 2767.77 2704.64 2641.51 2578.38 2515.25 2452.12 2389 2325.87 2262.74 2199.61 2136.48 2073.35 2010.22 1947.09 1883.97 1820.84 1757.71 1694.58 1631.45 1568.32 1524.02 1489.62 1455.22 1420.82 1386.42 1352.02 1317.63 1283.23 1248.83 1214.43 1180.03 1145.63 1111.24 1076.84 1042.44 1008.04 973.643 939.245 904.847 870.448 847.568 830.749 813.93 797.111 780.292 763.473 746.654 729.835 713.016 696.197 679.377 662.558 645.739 628.92 612.101 595.282 578.463 561.644 544.825 528.006 514.286 502.196 490.106 478.017 465.927 453.837 441.748 429.658 417.568 405.478 393.389 381.299 369.209 357.12 345.03 332.94 320.85 308.761 296.671
+beam_content += -14 5.9% 50 9450 30264.9 34305.2 34028.6 30703.8 26975.1 27044.5 27369.6 27746.9 28126.5 27779.2 27413.1 26882.3 26331.1 24512 22578 20644.1 18710.1 17191.8 15675.9 14695.9 13768.4 12933.3 12098.2 11433.4 10770.2 9830.67 8891.1 8107.64 7336.35 6880.7 6439.43 6213.76 5990.46 5520.28 5050.11 4809.09 4587.31 4365.54 4143.76 3868.78 3592.18 3381.59 3171.36 2938.91 2705.28 2495.05 2284.81 2307.97 2341.74 2375.5 2270.61 2164.19 2057.77 1951.34 1840.22 1728.04 1615.85 1503.66 1491 1494.39 1497.78 1482.94 1388.2 1293.45 1198.71 1138.77 1134.57 1130.36 1126.16 1108.57 1048.01 987.446 926.886 869.666 845.929 822.193 798.457 774.72 749.48 724.151 698.822 673.492 648.163 622.833 597.859 583.501 569.143 554.785 540.427 526.068 511.71 497.352 482.994 468.636 454.31 443.006 431.702 420.398 409.093 397.789 386.485 375.181 363.877 352.572 341.268 329.964 318.66 307.355 296.051 285.453 278.819 272.185 265.551 258.917 252.283 245.648 239.014 232.38 225.746 219.112 212.478 205.844 199.21 192.575 185.941 179.307 172.673 166.039 159.405 152.771 147.274 143.367 139.46 135.554 131.647 127.741 123.834 119.928 116.021 112.114 108.208 104.301 100.395 96.4881 92.5816 88.675 84.7684 80.8618 76.9553 73.0487 70.0752 68.666 67.2567 65.8475 64.4382 63.029 61.6197 60.2105 58.8012 57.392 55.9827 54.5735 53.1642 51.755 50.3457 48.9365 47.5272 46.118 44.7087 43.2995 42.5277 42.0615 41.5953 41.129 40.6628 40.1966 39.7303 39.2641 38.7979 38.3316 37.8654 37.3991 36.9329 36.4667 36.0004 35.5342 35.068 34.6017 34.1355
+beam_content += 12 1.1% 50 9450 682.68 1614.26 2547.81 3539.18 4497.06 4879.09 5261.12 5389.59 5493.08 5447.68 5400.58 5199.56 4992.35 4603.48 4214.62 3951.77 3690.69 3362.89 3035.08 2856.57 2678.66 2405.14 2135.59 1982.91 1833.19 1734.37 1635.54 1523.59 1411.2 1357.52 1304.49 1222.73 1139.13 1093.73 1048.33 1030.99 1015.27 983.975 952.147 928.287 904.697 897.463 890.229 866.417 841.434 820.08 798.725 782.205 766.49 750.775 731.959 712.76 693.561 674.362 651.394 627.623 603.851 580.08 558.781 537.857 516.933 496.01 477.714 460.207 442.699 425.192 409.079 393.749 378.419 363.089 352.524 342.774 333.024 323.273 310.985 297.728 284.472 271.215 259.332 251.523 243.713 235.904 228.094 220.285 212.476 204.683 197.583 190.483 183.383 176.282 169.182 162.082 154.982 147.882 140.782 135.001 131.267 127.533 123.799 120.065 116.331 112.597 108.863 105.129 101.396 97.6616 93.9277 90.1938 86.4599 82.7259 79.3897 77.2714 75.1531 73.0348 70.9164 68.7981 66.6798 64.5615 62.4431 60.3248 58.2065 56.0882 53.9698 51.8515 49.7332 47.6149 45.4965 43.3782 41.2599 39.1416 37.3594 36.4119 35.4644 34.5169 33.5694 32.6219 31.6744 30.7269 29.7794 28.8319 27.8844 26.9369 25.9894 25.042 24.0945 23.147 22.1995 21.252 20.3045 19.357 18.4095 17.8893 17.38 16.8706 16.3612 15.8519 15.3425 14.8332 14.3238 13.8145 13.3051 12.7958 12.2864 11.7771 11.2677 10.7583 10.249 9.73964 9.23028 8.72093 8.34094 8.18991 8.03887 7.88783 7.73679 7.58575 7.43471 7.28368 7.13264 6.9816 6.83056 6.67952 6.52848 6.37745 6.22641 6.07537 5.92433 5.77329 5.62226
+beam_content += -12 0.2% 50 9450 69.169 148.438 227.629 291.149 353.265 382.617 411.969 463.842 510.51 499.004 487.645 486.89 486.136 481.116 475.782 511.634 539.866 479.616 419.964 403.28 387.281 387.244 385.393 358.873 333.381 317.574 302.12 312.118 322.116 311.572 299.771 291.411 283.23 276.741 270.206 260.333 250.926 252.401 252.375 231.908 213.194 211.007 208.89 211.879 214.122 198.952 184.749 183.703 182.657 181.116 175.558 169.999 164.44 159.004 155.67 152.336 149.002 145.669 142.772 139.884 136.996 134.108 130.574 126.825 123.076 119.326 113.526 107.593 101.661 95.7284 94.1626 93.803 93.4434 93.0838 90.1456 87.1387 84.1318 81.1249 78.5051 76.3185 74.1319 71.9453 69.7588 67.5722 65.3856 63.899 62.4974 61.0958 59.6942 58.2926 56.891 55.4894 54.0878 52.6862 51.2846 49.8847 48.4864 47.0881 45.6899 44.2916 42.8933 41.4951 40.0968 38.6985 37.3003 35.902 34.5037 33.1055 31.7072 30.3089 29.547 28.8139 28.0808 27.3477 26.6146 25.8816 25.1485 24.4154 23.6823 22.9492 22.2161 21.483 20.7499 20.0168 19.2837 18.5506 17.8175 17.0844 16.3513 15.6873 15.2732 14.8591 14.445 14.031 13.6169 13.2028 12.7887 12.3747 11.9606 11.5465 11.1324 10.7184 10.3043 9.89021 9.47614 9.06206 8.64799 8.23391 7.81984 7.40576 7.17193 6.98789 6.80385 6.61981 6.43577 6.25173 6.06769 5.88365 5.69961 5.51557 5.33153 5.14749 4.96345 4.77941 4.59537 4.41133 4.22728 4.04324 3.8592 3.68106 3.63326 3.58546 3.53767 3.48987 3.44207 3.39427 3.34647 3.29868 3.25088 3.20308 3.15528 3.10748 3.05968 3.01189 2.96409 2.91629 2.86849 2.82069 2.77289
+beam_folder = flux
+beam_file = 1
+beam_file_first = 1
+beam_file_limit = 0
+beam_weighted = 0
+beam_offset = 0 0 0
+beam_placement = 0
+beam_test_only = 0
+target_type = 1
+nucleus_p = 6
+nucleus_n = 6
+nucleus_E_b = 34
+nucleus_kf = 220
+target_content = 1 0 1x 0 0 0
+target_content += 6 6 1x 27 220 2
+geo_file = target/ND280_v9r7p5.root
+geo_name = ND280Geometry_v9r7p5
+geo_volume =
+geo_o = 0 0 0
+geo_d = 2000 2000 5000
+nucleus_target = 2
+nucleus_model = 1
+nucleus_LFGkfScale = 1
+dyn_qel_cc = 1
+dyn_qel_nc = 0
+dyn_res_cc = 1
+dyn_res_nc = 0
+dyn_dis_cc = 1
+dyn_dis_nc = 0
+dyn_coh_cc = 0
+dyn_coh_nc = 0
+dyn_mec_cc = 1
+dyn_mec_nc = 0
+qel_kinematics = 0
+qel_vector_ff_set = 2
+bba07_AEp1 = 1
+bba07_AEp2 = 0.9927
+bba07_AEp3 = 0.9898
+bba07_AEp4 = 0.9975
+bba07_AEp5 = 0.9812
+bba07_AEp6 = 0.934
+bba07_AEp7 = 1
+bba07_AMp1 = 1
+bba07_AMp2 = 1.0011
+bba07_AMp3 = 0.9992
+bba07_AMp4 = 0.9974
+bba07_AMp5 = 1.001
+bba07_AMp6 = 1.0003
+bba07_AMp7 = 1
+bba07_AEn1 = 1
+bba07_AEn2 = 1.011
+bba07_AEn3 = 1.1392
+bba07_AEn4 = 1.0203
+bba07_AEn5 = 1.1093
+bba07_AEn6 = 1.5429
+bba07_AEn7 = 0.9706
+bba07_AMn1 = 1
+bba07_AMn2 = 0.9958
+bba07_AMn3 = 0.9877
+bba07_AMn4 = 1.0193
+bba07_AMn5 = 1.035
+bba07_AMn6 = 0.9164
+bba07_AMn7 = 0.73
+bba07_AAx1 = 1
+bba07_AAx2 = 0.9958
+bba07_AAx3 = 0.9877
+bba07_AAx4 = 1.0193
+bba07_AAx5 = 1.035
+bba07_AAx6 = 0.9164
+bba07_AAx7 = 0.73
+qel_axial_ff_set = 1
+qel_rpa = 0
+qel_strange = 1
+qel_strangeEM = 0
+delta_s = -0.15
+qel_cc_vector_mass = 1000
+qel_cc_axial_mass = 1200
+qel_nc_axial_mass = 1350
+qel_s_axial_mass = 1200
+qel_axial_2comp_gamma = 0.15
+qel_axial_2comp_alpha = 2
+qel_axial_3comp_theta = 0.15
+qel_axial_3comp_beta = 2
+flux_correction = 0
+sf_method = 0
+sf_pb = 0
+cc_smoothing = 0
+delta_FF_set = 1
+pion_axial_mass = 0.94
+pion_C5A = 1.19
+SPPBkgScale = 0
+delta_angular = 0
+spp_precision = 500
+res_dis_cut = 1600
+bkgrscaling = 0
+coh_mass_correction = 1
+coh_new = 1
+coh_kind = 2
+mec_kind = 3
+mec_ratio_pp = 0.6
+mec_ratio_ppp = 0.8
+mec_central_motion = 0
+mec_back_to_back_smearing = 0
+mec_pb_trials = 25
+MEC_pauli_blocking = 1
+MEC_cm_direction = 0
+kaskada_on = 1
+kaskada_w = 7
+kaskada_redo = 0
+kaskada_writeall = 0
+kaskada_newangle = 0
+formation_zone = fz
+tau = 8
+formation_length = 1
+first_step = 1
+step = 0.2
+xsec = 1
+pauli_blocking = 1
+mixed_order = 1
+rmin = 1
+rmax = 1
diff --git a/src/ArgoNeuT/ArgoNeuT_CCInc_XSec_1Dpmu_antinu.cxx b/src/ArgoNeuT/ArgoNeuT_CCInc_XSec_1Dpmu_antinu.cxx
index 25018ce..eb02d50 100644
--- a/src/ArgoNeuT/ArgoNeuT_CCInc_XSec_1Dpmu_antinu.cxx
+++ b/src/ArgoNeuT/ArgoNeuT_CCInc_XSec_1Dpmu_antinu.cxx
@@ -1,37 +1,37 @@
#include "ArgoNeuT_CCInc_XSec_1Dpmu_antinu.h"
//********************************************************************
ArgoNeuT_CCInc_XSec_1Dpmu_antinu::ArgoNeuT_CCInc_XSec_1Dpmu_antinu(
std::string inputfile, FitWeight *rw, std::string type,
std::string fakeDataFile)
//********************************************************************
{
measurementName = "ArgoNeuT_CCInc_XSec_1Dpmu_antinu";
default_types = "FIX/DIAG/CHI2";
plotTitles = "; p_{#mu} (GeV); d#sigma/dp_{#mu} (cm^{2} Ar^{-1} GeV^{-1})";
EnuMin = 0;
EnuMax = 50;
isDiag = true;
Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile);
SetDataValues(std::string(std::getenv("EXT_FIT")) +
"/data/ArgoNeuT/CCInc_dsig_dmumom_nubar.dat");
SetupDefaultHist();
scaleFactor = eventHist->Integral("width") * double(1E-38) / double(nevents) *
(40.0 /*Data is /Ar */) / TotalIntegratedFlux("width");
};
void ArgoNeuT_CCInc_XSec_1Dpmu_antinu::FillEventVariables(FitEvent *event) {
- X_VAR = FitUtils::GetHMPDG_4Mom(-13, event).Vect().Mag();
+ X_VAR = FitUtils::GetHMPDG_4Mom(-13, event).first.Vect().Mag();
return;
};
//********************************************************************
bool ArgoNeuT_CCInc_XSec_1Dpmu_antinu::isSignal(FitEvent *event)
//********************************************************************
{
return SignalDef::isCCInc_ArgoNeuT(event, true);
}
diff --git a/src/ArgoNeuT/ArgoNeuT_CCInc_XSec_1Dpmu_nu.cxx b/src/ArgoNeuT/ArgoNeuT_CCInc_XSec_1Dpmu_nu.cxx
index e9d035b..fc53735 100644
--- a/src/ArgoNeuT/ArgoNeuT_CCInc_XSec_1Dpmu_nu.cxx
+++ b/src/ArgoNeuT/ArgoNeuT_CCInc_XSec_1Dpmu_nu.cxx
@@ -1,42 +1,60 @@
+// 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 <http://www.gnu.org/licenses/>.
+*******************************************************************************/
#include "ArgoNeuT_CCInc_XSec_1Dpmu_nu.h"
//********************************************************************
ArgoNeuT_CCInc_XSec_1Dpmu_nu::ArgoNeuT_CCInc_XSec_1Dpmu_nu(
std::string inputfile, FitWeight *rw, std::string type,
std::string fakeDataFile)
//********************************************************************
{
measurementName = "ArgoNeuT_CCInc_XSec_1Dpmu_nu";
default_types = "FIX/DIAG/CHI2";
plotTitles = "; p_{#mu} (GeV); d#sigma/dp_{#mu} (cm^{2} Ar^{-1} GeV^{-1})";
EnuMin = 0;
EnuMax = 50;
isDiag = true;
Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile);
SetDataValues(std::string(std::getenv("EXT_FIT")) +
"/data/ArgoNeuT/CCInc_dsig_dmumom_nu.dat");
dataHist->Scale(1E-38);
dataTrue->Scale(1E-38);
SetupDefaultHist();
scaleFactor = eventHist->Integral("width") * double(1E-38) / double(nevents) *
(40.0 /*Data is /Ar */) / TotalIntegratedFlux("width");
- std::cout << "SF: " << eventHist->Integral("width") << " " << nevents << " " << TotalIntegratedFlux("width") << std::endl;
+ std::cout << "SF: " << eventHist->Integral("width") << " " << nevents << " "
+ << TotalIntegratedFlux("width") << std::endl;
};
void ArgoNeuT_CCInc_XSec_1Dpmu_nu::FillEventVariables(FitEvent *event) {
- X_VAR = FitUtils::GetHMPDG_4Mom(13, event).Vect().Mag()/1000.0;
+ X_VAR = FitUtils::GetHMPDG_4Mom(13, event).first.Vect().Mag() / 1000.0;
return;
};
//********************************************************************
bool ArgoNeuT_CCInc_XSec_1Dpmu_nu::isSignal(FitEvent *event)
//********************************************************************
{
return SignalDef::isCCInc_ArgoNeuT(event);
}
-
diff --git a/src/ArgoNeuT/ArgoNeuT_CCInc_XSec_1Dpmu_nu.h b/src/ArgoNeuT/ArgoNeuT_CCInc_XSec_1Dpmu_nu.h
index 2983de6..924001f 100644
--- a/src/ArgoNeuT/ArgoNeuT_CCInc_XSec_1Dpmu_nu.h
+++ b/src/ArgoNeuT/ArgoNeuT_CCInc_XSec_1Dpmu_nu.h
@@ -1,17 +1,36 @@
+// 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 <http://www.gnu.org/licenses/>.
+*******************************************************************************/
+
#ifndef ArgoNeuT_CCInc_XSec_1Dpmu_nu_H_SEEN
#define ArgoNeuT_CCInc_XSec_1Dpmu_nu_H_SEEN
#include "Measurement1D.h"
class ArgoNeuT_CCInc_XSec_1Dpmu_nu : public Measurement1D {
public:
ArgoNeuT_CCInc_XSec_1Dpmu_nu(std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile);
virtual ~ArgoNeuT_CCInc_XSec_1Dpmu_nu() {};
void FillEventVariables(FitEvent *event);
bool isSignal(FitEvent *event);
private:
};
#endif
diff --git a/src/ArgoNeuT/ArgoNeuT_CCInc_XSec_1Dthetamu_antinu.cxx b/src/ArgoNeuT/ArgoNeuT_CCInc_XSec_1Dthetamu_antinu.cxx
index 055f88f..ef64a3e 100644
--- a/src/ArgoNeuT/ArgoNeuT_CCInc_XSec_1Dthetamu_antinu.cxx
+++ b/src/ArgoNeuT/ArgoNeuT_CCInc_XSec_1Dthetamu_antinu.cxx
@@ -1,37 +1,40 @@
#include "ArgoNeuT_CCInc_XSec_1Dthetamu_antinu.h"
//********************************************************************
ArgoNeuT_CCInc_XSec_1Dthetamu_antinu::ArgoNeuT_CCInc_XSec_1Dthetamu_antinu(
std::string inputfile, FitWeight *rw, std::string type,
std::string fakeDataFile)
//********************************************************************
{
measurementName = "ArgoNeuT_CCInc_XSec_1Dthetamu_antinu";
default_types = "FIX/DIAG/CHI2";
- plotTitles = "; theta_{#mu} (degrees); d#sigma/d#theta_{#mu} (cm^{2} Ar^{-1} degrees^{-1})";
+ plotTitles =
+ "; theta_{#mu} (degrees); d#sigma/d#theta_{#mu} (cm^{2} Ar^{-1} "
+ "degrees^{-1})";
EnuMin = 0;
EnuMax = 50;
isDiag = true;
Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile);
SetDataValues(std::string(std::getenv("EXT_FIT")) +
"/data/ArgoNeuT/CCInc_dsig_dthetamu_nubar.dat");
SetupDefaultHist();
scaleFactor = eventHist->Integral("width") * double(1E-38) / double(nevents) *
(40.0 /*Data is /Ar */) / TotalIntegratedFlux("width");
};
void ArgoNeuT_CCInc_XSec_1Dthetamu_antinu::FillEventVariables(FitEvent *event) {
- X_VAR = (FitUtils::GetHMPDG_4Mom(-13, event).Vect().Theta()/TMath::Pi())*180.;
+ X_VAR =
+ (FitUtils::GetHMPDG_4Mom(-13, event).first.Vect().Theta() / TMath::Pi()) *
+ 180.;
return;
};
//********************************************************************
bool ArgoNeuT_CCInc_XSec_1Dthetamu_antinu::isSignal(FitEvent *event)
//********************************************************************
{
return SignalDef::isCCInc_ArgoNeuT(event, true);
}
-
diff --git a/src/ArgoNeuT/ArgoNeuT_CCInc_XSec_1Dthetamu_nu.cxx b/src/ArgoNeuT/ArgoNeuT_CCInc_XSec_1Dthetamu_nu.cxx
index 98fcb88..aab2a92 100644
--- a/src/ArgoNeuT/ArgoNeuT_CCInc_XSec_1Dthetamu_nu.cxx
+++ b/src/ArgoNeuT/ArgoNeuT_CCInc_XSec_1Dthetamu_nu.cxx
@@ -1,36 +1,39 @@
#include "ArgoNeuT_CCInc_XSec_1Dthetamu_nu.h"
//********************************************************************
ArgoNeuT_CCInc_XSec_1Dthetamu_nu::ArgoNeuT_CCInc_XSec_1Dthetamu_nu(
std::string inputfile, FitWeight *rw, std::string type,
std::string fakeDataFile)
//********************************************************************
{
measurementName = "ArgoNeuT_CCInc_XSec_1Dthetamu_nu";
default_types = "FIX/DIAG/CHI2";
- plotTitles = "; theta_{#mu} (degrees); d#sigma/d#theta_{#mu} (cm^{2} Ar^{-1} degrees^{-1})";
+ plotTitles =
+ "; theta_{#mu} (degrees); d#sigma/d#theta_{#mu} (cm^{2} Ar^{-1} "
+ "degrees^{-1})";
EnuMin = 0;
EnuMax = 50;
Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile);
SetDataValues(std::string(std::getenv("EXT_FIT")) +
"/data/ArgoNeuT/CCInc_dsig_dthetamu_nu.dat");
SetupDefaultHist();
scaleFactor = eventHist->Integral("width") * double(1E-38) / double(nevents) *
(40.0 /*Data is /Ar */) / TotalIntegratedFlux("width");
};
void ArgoNeuT_CCInc_XSec_1Dthetamu_nu::FillEventVariables(FitEvent *event) {
- X_VAR = (FitUtils::GetHMPDG_4Mom(13, event).Vect().Theta()/TMath::Pi())*180.;
+ X_VAR =
+ (FitUtils::GetHMPDG_4Mom(13, event).first.Vect().Theta() / TMath::Pi()) *
+ 180.;
return;
};
//********************************************************************
bool ArgoNeuT_CCInc_XSec_1Dthetamu_nu::isSignal(FitEvent *event)
//********************************************************************
{
return SignalDef::isCCInc_ArgoNeuT(event, true);
}
-
diff --git a/src/FCN/SampleList.cxx b/src/FCN/SampleList.cxx
index 43face5..e55dd58 100644
--- a/src/FCN/SampleList.cxx
+++ b/src/FCN/SampleList.cxx
@@ -1,513 +1,519 @@
#include "SampleList.h"
//! Functions to make it easier for samples to be created and handled.
namespace SampleUtils {
//! Create a given sample given its name, file, type, fakdata(fkdt) file and the
//! current rw engine and push it back into the list fChain.
bool LoadSample(std::list<MeasurementBase*>* fChain, std::string name,
std::string file, std::string type, std::string fkdt,
FitWeight* rw) {
/*
ANL CCQE Samples
*/
if (!name.compare("ANL_CCQE_XSec_1DEnu_nu") ||
!name.compare("ANL_CCQE_XSec_1DEnu_nu_PRL31") ||
!name.compare("ANL_CCQE_XSec_1DEnu_nu_PRD16")) {
fChain->push_back(new ANL_CCQE_XSec_1DEnu_nu(name, file, rw, type, fkdt));
} else if (!name.compare("ANL_CCQE_Evt_1DQ2_nu") ||
!name.compare("ANL_CCQE_Evt_1DQ2_nu_PRL31") ||
!name.compare("ANL_CCQE_Evt_1DQ2_nu_PRD16")) {
fChain->push_back(new ANL_CCQE_Evt_1DQ2_nu(name, file, rw, type, fkdt));
/*
ANL CC1ppip samples
*/
} else if (!name.compare("ANL_CC1ppip_XSec_1DEnu_nu")) {
fChain->push_back(new ANL_CC1ppip_XSec_1DEnu_nu(file, rw, type, fkdt));
} else if (!name.compare("ANL_CC1ppip_XSec_1DQ2_nu")) {
fChain->push_back(new ANL_CC1ppip_XSec_1DQ2_nu(file, rw, type, fkdt));
} else if (!name.compare("ANL_CC1ppip_Evt_1DQ2_nu")) {
fChain->push_back(new ANL_CC1ppip_Evt_1DQ2_nu(file, rw, type, fkdt));
} else if (!name.compare("ANL_CC1ppip_Evt_1Dppi_nu")) {
fChain->push_back(new ANL_CC1ppip_Evt_1Dppi_nu(file, rw, type, fkdt));
} else if (!name.compare("ANL_CC1ppip_Evt_1Dthpr_nu")) {
fChain->push_back(new ANL_CC1ppip_Evt_1Dthpr_nu(file, rw, type, fkdt));
} else if (!name.compare("ANL_CC1ppip_Evt_1DcosmuStar_nu")) {
fChain->push_back(new ANL_CC1ppip_Evt_1DcosmuStar_nu(file, rw, type, fkdt));
} else if (!name.compare("ANL_CC1ppip_Evt_1DcosthAdler_nu")) {
fChain->push_back(
new ANL_CC1ppip_Evt_1DcosthAdler_nu(file, rw, type, fkdt));
} else if (!name.compare("ANL_CC1ppip_Evt_1Dphi_nu")) {
fChain->push_back(new ANL_CC1ppip_Evt_1Dphi_nu(file, rw, type, fkdt));
/*
ANL CC1npip sample
*/
} else if (!name.compare("ANL_CC1npip_XSec_1DEnu_nu")) {
fChain->push_back(new ANL_CC1npip_XSec_1DEnu_nu(file, rw, type, fkdt));
} else if (!name.compare("ANL_CC1npip_Evt_1DQ2_nu")) {
fChain->push_back(new ANL_CC1npip_Evt_1DQ2_nu(file, rw, type, fkdt));
} else if (!name.compare("ANL_CC1npip_Evt_1Dppi_nu")) {
fChain->push_back(new ANL_CC1npip_Evt_1Dppi_nu(file, rw, type, fkdt));
} else if (!name.compare("ANL_CC1npip_Evt_1DcosmuStar_nu")) {
fChain->push_back(new ANL_CC1npip_Evt_1DcosmuStar_nu(file, rw, type, fkdt));
/*
ANL CC1pi0 sample
*/
} else if (!name.compare("ANL_CC1pi0_XSec_1DEnu_nu")) {
fChain->push_back(new ANL_CC1pi0_XSec_1DEnu_nu(file, rw, type, fkdt));
} else if (!name.compare("ANL_CC1pi0_Evt_1DQ2_nu")) {
fChain->push_back(new ANL_CC1pi0_Evt_1DQ2_nu(file, rw, type, fkdt));
} else if (!name.compare("ANL_CC1pi0_Evt_1DcosmuStar_nu")) {
fChain->push_back(new ANL_CC1pi0_Evt_1DcosmuStar_nu(file, rw, type, fkdt));
/*
ANL NC1npip sample
*/
} else if (!name.compare("ANL_NC1npip_Evt_1Dppi_nu")) {
fChain->push_back(new ANL_NC1npip_Evt_1Dppi_nu(file, rw, type, fkdt));
} else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dpmu_antinu")) {
fChain->push_back(
new ArgoNeuT_CCInc_XSec_1Dpmu_antinu(file, rw, type, fkdt));
} else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dpmu_nu")) {
fChain->push_back(new ArgoNeuT_CCInc_XSec_1Dpmu_nu(file, rw, type, fkdt));
} else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dthetamu_antinu")) {
fChain->push_back(
new ArgoNeuT_CCInc_XSec_1Dthetamu_antinu(file, rw, type, fkdt));
} else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dthetamu_nu")) {
fChain->push_back(
new ArgoNeuT_CCInc_XSec_1Dthetamu_nu(file, rw, type, fkdt));
/*
BNL Samples
*/
} else if (!name.compare("BNL_CCQE_XSec_1DEnu_nu")) {
fChain->push_back(new BNL_CCQE_XSec_1DEnu_nu(file, rw, type, fkdt));
} else if (!name.compare("BNL_CCQE_Evt_1DQ2_nu")) {
fChain->push_back(new BNL_CCQE_Evt_1DQ2_nu(file, rw, type, fkdt));
/*
BNL CC1ppip samples
*/
} else if (!name.compare("BNL_CC1ppip_XSec_1DEnu_nu")) {
fChain->push_back(new BNL_CC1ppip_XSec_1DEnu_nu(file, rw, type, fkdt));
} else if (!name.compare("BNL_CC1ppip_Evt_1DQ2_nu")) {
fChain->push_back(new BNL_CC1ppip_Evt_1DQ2_nu(file, rw, type, fkdt));
} else if (!name.compare("BNL_CC1ppip_Evt_1DcosthAdler_nu")) {
fChain->push_back(
new BNL_CC1ppip_Evt_1DcosthAdler_nu(file, rw, type, fkdt));
} else if (!name.compare("BNL_CC1ppip_Evt_1Dphi_nu")) {
fChain->push_back(new BNL_CC1ppip_Evt_1Dphi_nu(file, rw, type, fkdt));
/*
BNL CC1npip samples
*/
} else if (!name.compare("BNL_CC1npip_XSec_1DEnu_nu")) {
fChain->push_back(new BNL_CC1npip_XSec_1DEnu_nu(file, rw, type, fkdt));
} else if (!name.compare("BNL_CC1npip_Evt_1DQ2_nu")) {
fChain->push_back(new BNL_CC1npip_Evt_1DQ2_nu(file, rw, type, fkdt));
/*
BNL CC1pi0 samples
*/
} else if (!name.compare("BNL_CC1pi0_XSec_1DEnu_nu")) {
fChain->push_back(new BNL_CC1pi0_XSec_1DEnu_nu(file, rw, type, fkdt));
} else if (!name.compare("BNL_CC1pi0_Evt_1DQ2_nu")) {
fChain->push_back(new BNL_CC1pi0_Evt_1DQ2_nu(file, rw, type, fkdt));
/*
FNAL Samples
*/
} else if (!name.compare("FNAL_CCQE_Evt_1DQ2_nu")) {
fChain->push_back(new FNAL_CCQE_Evt_1DQ2_nu(file, rw, type, fkdt));
/*
FNAL CC1ppip
*/
} else if (!name.compare("FNAL_CC1ppip_XSec_1DEnu_nu")) {
fChain->push_back(new FNAL_CC1ppip_XSec_1DEnu_nu(file, rw, type, fkdt));
} else if (!name.compare("FNAL_CC1ppip_XSec_1DQ2_nu")) {
fChain->push_back(new FNAL_CC1ppip_XSec_1DQ2_nu(file, rw, type, fkdt));
} else if (!name.compare("FNAL_CC1ppip_Evt_1DQ2_nu")) {
fChain->push_back(new FNAL_CC1ppip_Evt_1DQ2_nu(file, rw, type, fkdt));
/*
FNAL CC1ppim
*/
// } else if (!name.compare("FNAL_CC1ppim_XSec_1DEnu_antinu")) {
// fChain->push_back(new FNAL_CC1ppim_XSec_1DEnu_antinu(file, rw, type,
// fkdt));
/*
BEBC Samples
*/
} else if (!name.compare("BEBC_CCQE_XSec_1DQ2_nu")) {
fChain->push_back(new BEBC_CCQE_XSec_1DQ2_nu(name, file, rw, type, fkdt));
/*
BEBC CC1ppip samples
*/
} else if (!name.compare("BEBC_CC1ppip_XSec_1DEnu_nu")) {
fChain->push_back(new BEBC_CC1ppip_XSec_1DEnu_nu(file, rw, type, fkdt));
} else if (!name.compare("BEBC_CC1ppip_XSec_1DQ2_nu")) {
fChain->push_back(new BEBC_CC1ppip_XSec_1DQ2_nu(file, rw, type, fkdt));
/*
BEBC CC1npip samples
*/
} else if (!name.compare("BEBC_CC1npip_XSec_1DEnu_nu")) {
fChain->push_back(new BEBC_CC1npip_XSec_1DEnu_nu(file, rw, type, fkdt));
} else if (!name.compare("BEBC_CC1npip_XSec_1DQ2_nu")) {
fChain->push_back(new BEBC_CC1npip_XSec_1DQ2_nu(file, rw, type, fkdt));
/*
BEBC CC1pi0 samples
*/
} else if (!name.compare("BEBC_CC1pi0_XSec_1DEnu_nu")) {
fChain->push_back(new BEBC_CC1pi0_XSec_1DEnu_nu(file, rw, type, fkdt));
} else if (!name.compare("BEBC_CC1pi0_XSec_1DQ2_nu")) {
fChain->push_back(new BEBC_CC1pi0_XSec_1DQ2_nu(file, rw, type, fkdt));
/*
BEBC CC1npim samples
*/
} else if (!name.compare("BEBC_CC1npim_XSec_1DEnu_antinu")) {
fChain->push_back(new BEBC_CC1npim_XSec_1DEnu_antinu(file, rw, type, fkdt));
} else if (!name.compare("BEBC_CC1npim_XSec_1DQ2_antinu")) {
fChain->push_back(new BEBC_CC1npim_XSec_1DQ2_antinu(file, rw, type, fkdt));
/*
BEBC CC1ppim samples
*/
} else if (!name.compare("BEBC_CC1ppim_XSec_1DEnu_antinu")) {
fChain->push_back(new BEBC_CC1ppim_XSec_1DEnu_antinu(file, rw, type, fkdt));
} else if (!name.compare("BEBC_CC1ppim_XSec_1DQ2_antinu")) {
fChain->push_back(new BEBC_CC1ppim_XSec_1DQ2_antinu(file, rw, type, fkdt));
/*
GGM CC1ppip samples
*/
} else if (!name.compare("GGM_CC1ppip_XSec_1DEnu_nu")) {
fChain->push_back(new GGM_CC1ppip_XSec_1DEnu_nu(file, rw, type, fkdt));
} else if (!name.compare("GGM_CC1ppip_Evt_1DQ2_nu")) {
fChain->push_back(new GGM_CC1ppip_Evt_1DQ2_nu(file, rw, type, fkdt));
/*
MiniBooNE Samples
*/
/*
CCQE
*/
} else if (!name.compare("MiniBooNE_CCQE_XSec_1DQ2_nu") ||
!name.compare("MiniBooNE_CCQELike_XSec_1DQ2_nu")) {
fChain->push_back(
new MiniBooNE_CCQE_XSec_1DQ2_nu(name, file, rw, type, fkdt));
} else if (!name.compare("MiniBooNE_CCQE_XSec_1DQ2_antinu") ||
!name.compare("MiniBooNE_CCQELike_XSec_1DQ2_antinu")) {
fChain->push_back(
new MiniBooNE_CCQE_XSec_1DQ2_antinu(name, file, rw, type, fkdt));
} else if (!name.compare("MiniBooNE_CCQE_XSec_2DTcos_nu") ||
!name.compare("MiniBooNE_CCQELike_XSec_2DTcos_nu")) {
fChain->push_back(
new MiniBooNE_CCQE_XSec_2DTcos_nu(name, file, rw, type, fkdt));
} else if (!name.compare("MiniBooNE_CCQE_XSec_2DTcos_antinu") ||
!name.compare("MiniBooNE_CCQELike_XSec_2DTcos_antinu")) {
fChain->push_back(
new MiniBooNE_CCQE_XSec_2DTcos_antinu(name, file, rw, type, fkdt));
/*
MiniBooNE CC1pi+
*/
// 1D
} else if (!name.compare("MiniBooNE_CC1pip_XSec_1DEnu_nu")) {
fChain->push_back(new MiniBooNE_CC1pip_XSec_1DEnu_nu(file, rw, type, fkdt));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_1DQ2_nu")) {
fChain->push_back(new MiniBooNE_CC1pip_XSec_1DQ2_nu(file, rw, type, fkdt));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_1DTpi_nu")) {
fChain->push_back(new MiniBooNE_CC1pip_XSec_1DTpi_nu(file, rw, type, fkdt));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_1DTu_nu")) {
fChain->push_back(new MiniBooNE_CC1pip_XSec_1DTu_nu(file, rw, type, fkdt));
// 2D
} else if (!name.compare("MiniBooNE_CC1pip_XSec_2DQ2Enu_nu")) {
fChain->push_back(
new MiniBooNE_CC1pip_XSec_2DQ2Enu_nu(file, rw, type, fkdt));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTpiCospi_nu")) {
fChain->push_back(
new MiniBooNE_CC1pip_XSec_2DTpiCospi_nu(file, rw, type, fkdt));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTpiEnu_nu")) {
fChain->push_back(
new MiniBooNE_CC1pip_XSec_2DTpiEnu_nu(file, rw, type, fkdt));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTuCosmu_nu")) {
fChain->push_back(
new MiniBooNE_CC1pip_XSec_2DTuCosmu_nu(file, rw, type, fkdt));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTuEnu_nu")) {
fChain->push_back(
new MiniBooNE_CC1pip_XSec_2DTuEnu_nu(file, rw, type, fkdt));
/*
MiniBooNE CC1pi0
*/
} else if (!name.compare("MiniBooNE_CC1pi0_XSec_1DEnu_nu")) {
fChain->push_back(new MiniBooNE_CC1pi0_XSec_1DEnu_nu(file, rw, type, fkdt));
} else if (!name.compare("MiniBooNE_CC1pi0_XSec_1DQ2_nu")) {
fChain->push_back(new MiniBooNE_CC1pi0_XSec_1DQ2_nu(file, rw, type, fkdt));
} else if (!name.compare("MiniBooNE_CC1pi0_XSec_1DTu_nu")) {
fChain->push_back(new MiniBooNE_CC1pi0_XSec_1DTu_nu(file, rw, type, fkdt));
} else if (!name.compare("MiniBooNE_CC1pi0_XSec_1Dcosmu_nu")) {
fChain->push_back(
new MiniBooNE_CC1pi0_XSec_1Dcosmu_nu(file, rw, type, fkdt));
} else if (!name.compare("MiniBooNE_CC1pi0_XSec_1Dcospi0_nu")) {
fChain->push_back(
new MiniBooNE_CC1pi0_XSec_1Dcospi0_nu(file, rw, type, fkdt));
} else if (!name.compare("MiniBooNE_CC1pi0_XSec_1Dppi0_nu")) {
fChain->push_back(
new MiniBooNE_CC1pi0_XSec_1Dppi0_nu(file, rw, type, fkdt));
/*
MiniBooNE NCEL
*/
} else if (!name.compare("MiniBooNE_NCEL_XSec_Treco_nu")) {
std::cerr
<< "MiniBooNE_NCEL_XSec_Treco_nu not implemented in current interface."
<< std::endl;
throw 5;
// fChain->push_back(new MiniBooNE_NCEL_XSec_Treco_nu(file, rw, type,
// fkdt));
/*
MINERvA Samples
*/
} else if (!name.compare("MINERvA_CCQE_XSec_1DQ2_nu") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_nu_20deg") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_nu_newflux") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_nu_20deg_newflux")) {
fChain->push_back(
new MINERvA_CCQE_XSec_1DQ2_nu(name, file, rw, type, fkdt));
} else if (!name.compare("MINERvA_CCQE_XSec_1DQ2_antinu") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_20deg") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_newflux") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_20deg_newflux")) {
fChain->push_back(
new MINERvA_CCQE_XSec_1DQ2_antinu(name, file, rw, type, fkdt));
} else if (!name.compare("MINERvA_CCQE_XSec_1DQ2_joint") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_joint_20deg") ||
// !name.compare("MINERvA_CCQE_XSec_1DQ2_joint_newflux") ||
// //< Not Current Supported/Stable
!name.compare("MINERvA_CCQE_XSec_1DQ2_joint_20deg_newflux")) {
fChain->push_back(
new MINERvA_CCQE_XSec_1DQ2_joint(name, file, rw, type, fkdt));
} else if (!name.compare("MINERvA_CC0pi_XSec_1DEe_nue")) {
fChain->push_back(new MINERvA_CC0pi_XSec_1DEe_nue(file, rw, type, fkdt));
} else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_nue")) {
fChain->push_back(new MINERvA_CC0pi_XSec_1DQ2_nue(file, rw, type, fkdt));
} else if (!name.compare("MINERvA_CC0pi_XSec_1DThetae_nue")) {
fChain->push_back(
new MINERvA_CC0pi_XSec_1DThetae_nue(file, rw, type, fkdt));
} else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_nu_proton")) {
fChain->push_back(
new MINERvA_CC0pi_XSec_1DQ2_nu_proton(file, rw, type, fkdt));
/*
CC1pi+
*/
} else if (!name.compare("MINERvA_CC1pip_XSec_1DTpi_20deg_nu")) {
fChain->push_back(
new MINERvA_CC1pip_XSec_1DTpi_20deg_nu(file, rw, type, fkdt));
} else if (!name.compare("MINERvA_CC1pip_XSec_1DTpi_nu")) {
fChain->push_back(new MINERvA_CC1pip_XSec_1DTpi_nu(file, rw, type, fkdt));
} else if (!name.compare("MINERvA_CC1pip_XSec_1Dth_20deg_nu")) {
fChain->push_back(
new MINERvA_CC1pip_XSec_1Dth_20deg_nu(file, rw, type, fkdt));
} else if (!name.compare("MINERvA_CC1pip_XSec_1Dth_nu")) {
fChain->push_back(new MINERvA_CC1pip_XSec_1Dth_nu(file, rw, type, fkdt));
/*
CCNpi+
*/
} else if (!name.compare("MINERvA_CCNpip_XSec_1Dth_nu")) {
fChain->push_back(new MINERvA_CCNpip_XSec_1Dth_nu(file, rw, type, fkdt));
} else if (!name.compare("MINERvA_CCNpip_XSec_1Dth_20deg_nu")) {
fChain->push_back(
new MINERvA_CCNpip_XSec_1Dth_20deg_nu(file, rw, type, fkdt));
} else if (!name.compare("MINERvA_CCNpip_XSec_1DTpi_nu")) {
fChain->push_back(new MINERvA_CCNpip_XSec_1DTpi_nu(file, rw, type, fkdt));
} else if (!name.compare("MINERvA_CCNpip_XSec_1DTpi_20deg_nu")) {
fChain->push_back(
new MINERvA_CCNpip_XSec_1DTpi_20deg_nu(file, rw, type, fkdt));
} else if (!name.compare("MINERvA_CCNpip_XSec_1Dthmu_nu")) {
fChain->push_back(new MINERvA_CCNpip_XSec_1Dthmu_nu(file, rw, type, fkdt));
} else if (!name.compare("MINERvA_CCNpip_XSec_1Dpmu_nu")) {
fChain->push_back(new MINERvA_CCNpip_XSec_1Dpmu_nu(file, rw, type, fkdt));
} else if (!name.compare("MINERvA_CCNpip_XSec_1DQ2_nu")) {
fChain->push_back(new MINERvA_CCNpip_XSec_1DQ2_nu(file, rw, type, fkdt));
} else if (!name.compare("MINERvA_CCNpip_XSec_1DEnu_nu")) {
fChain->push_back(new MINERvA_CCNpip_XSec_1DEnu_nu(file, rw, type, fkdt));
/*
CC1pi0
*/
} else if (!name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu")) {
fChain->push_back(
new MINERvA_CC1pi0_XSec_1Dth_antinu(file, rw, type, fkdt));
} else if (!name.compare("MINERvA_CC1pi0_XSec_1Dppi0_antinu")) {
fChain->push_back(
new MINERvA_CC1pi0_XSec_1Dppi0_antinu(file, rw, type, fkdt));
} else if (!name.compare("MINERvA_CC1pi0_XSec_1DQ2_antinu")) {
fChain->push_back(
new MINERvA_CC1pi0_XSec_1DQ2_antinu(file, rw, type, fkdt));
} else if (!name.compare("MINERvA_CC1pi0_XSec_1Dthmu_antinu")) {
fChain->push_back(
new MINERvA_CC1pi0_XSec_1Dthmu_antinu(file, rw, type, fkdt));
} else if (!name.compare("MINERvA_CC1pi0_XSec_1Dpmu_antinu")) {
fChain->push_back(
new MINERvA_CC1pi0_XSec_1Dpmu_antinu(file, rw, type, fkdt));
} else if (!name.compare("MINERvA_CC1pi0_XSec_1DEnu_antinu")) {
fChain->push_back(
new MINERvA_CC1pi0_XSec_1DEnu_antinu(file, rw, type, fkdt));
/*
CCINC
*/
} else if (!name.compare("MINERvA_CCinc_XSec_2DEavq3_nu")) {
fChain->push_back(new MINERvA_CCinc_XSec_2DEavq3_nu(file, rw, type, fkdt));
} else if (!name.compare("MINERvA_CCinc_XSec_1Dx_ratio_C12_CH") ||
!name.compare("MINERvA_CCinc_XSec_1Dx_ratio_Fe56_CH") ||
!name.compare("MINERvA_CCinc_XSec_1Dx_ratio_Pb208_CH")) {
fChain->push_back(
new MINERvA_CCinc_XSec_1Dx_ratio(name, file, rw, type, fkdt));
} else if (!name.compare("MINERvA_CCinc_XSec_1DEnu_ratio_C12_CH") ||
!name.compare("MINERvA_CCinc_XSec_1DEnu_ratio_Fe56_CH") ||
!name.compare("MINERvA_CCinc_XSec_1DEnu_ratio_Pb208_CH")) {
fChain->push_back(
new MINERvA_CCinc_XSec_1DEnu_ratio(name, file, rw, type, fkdt));
/*
T2K Samples
*/
} else if (!name.compare("T2K_CC0pi_XSec_2DPcos_nu") ||
!name.compare("T2K_CC0pi_XSec_2DPcos_nu_I") ||
!name.compare("T2K_CC0pi_XSec_2DPcos_nu_II")) {
fChain->push_back(new T2K_CC0pi_XSec_2DPcos_nu(name, file, rw, type, fkdt));
/*
T2K CC1pi+ CH samples
*/
} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dpmu_nu")) {
fChain->push_back(new T2K_CC1pip_CH_XSec_1Dpmu_nu(file, rw, type, fkdt));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dppi_nu")) {
fChain->push_back(new T2K_CC1pip_CH_XSec_1Dppi_nu(file, rw, type, fkdt));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1DQ2_nu")) {
fChain->push_back(new T2K_CC1pip_CH_XSec_1DQ2_nu(file, rw, type, fkdt));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dq3_nu")) {
fChain->push_back(new T2K_CC1pip_CH_XSec_1Dq3_nu(file, rw, type, fkdt));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthmupi_nu")) {
fChain->push_back(new T2K_CC1pip_CH_XSec_1Dthmupi_nu(file, rw, type, fkdt));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthpi_nu")) {
fChain->push_back(new T2K_CC1pip_CH_XSec_1Dthpi_nu(file, rw, type, fkdt));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthq3pi_nu")) {
fChain->push_back(new T2K_CC1pip_CH_XSec_1Dthq3pi_nu(file, rw, type, fkdt));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1DWrec_nu")) {
fChain->push_back(new T2K_CC1pip_CH_XSec_1DWrec_nu(file, rw, type, fkdt));
/*
+ T2K CC0pi + np CH samples
+ */
+ } else if (!name.compare("T2K_CC0pinp_STV_XSec_1Ddpt_nu")) {
+ fChain->push_back(new T2K_CC0pinp_STV_XSec_1Ddpt_nu(file, rw, type, fkdt));
+
+ /*
K2K Samples
*/
} else if (!name.compare("K2K_CC0pi_XSec_1DCosThetaMu_nu_1trk") ||
!name.compare("K2K_CC0pi_XSec_1DCosThetaMu_nu_2trkQE") ||
!name.compare("K2K_CC0pi_XSec_1DCosThetaMu_nu_2trkNonQE")) {
fChain->push_back(
new K2K_CC0pi_XSec_1DCosThetaMu_nu_subtrks(name, file, rw, type, fkdt));
} else if (!name.compare("K2KI_CC0pi_XSec_1DQ2_nu_1trk") ||
!name.compare("K2KI_CC0pi_XSec_1DQ2_nu_2trkQE") ||
!name.compare("K2KI_CC0pi_XSec_1DQ2_nu_2trkNonQE") ||
!name.compare("K2KIIa_CC0pi_XSec_1DQ2_nu_1trk") ||
!name.compare("K2KIIa_CC0pi_XSec_1DQ2_nu_2trkQE") ||
!name.compare("K2KIIa_CC0pi_XSec_1DQ2_nu_2trkNonQE")) {
fChain->push_back(
new K2K_CC0pi_XSec_1DQ2_nu_subtrks(name, file, rw, type, fkdt));
} else if (!name.compare("K2K_CC0pi_XSec_1DThetaMu_nu_Ntrks")) {
fChain->push_back(
new K2K_CC0pi_XSec_1DThetaMu_nu_Ntrks(file, rw, type, fkdt));
} else if (!name.compare("K2K_CC0pi_XSec_1DPmu_nu_Ntrks")) {
fChain->push_back(new K2K_CC0pi_XSec_1DPmu_nu_Ntrks(file, rw, type, fkdt));
} else if (!name.compare("K2K_CC0pi_XSec_1DDelPhi_nu_Ntrks")) {
fChain->push_back(
new K2K_CC0pi_XSec_1DDelPhi_nu_Ntrks(file, rw, type, fkdt));
/*
NC1pi0
*/
} else if (!name.compare("K2K_NC1pi0_Evt_1Dppi0_nu")) {
fChain->push_back(new K2K_NC1pi0_Evt_1Dppi0_nu(file, rw, type, fkdt));
/*
Stat Samples (Pulls etc)
*/
} else if (name.find("parameter_pulls") != std::string::npos) {
fChain->push_back(new parameter_pulls(name, file, rw, type, fkdt));
/*
Fake Studies
*/
} else if (name.find("ExpMultDist_CCQE_XSec_1D") != std::string::npos &&
name.find("_FakeStudy") != std::string::npos) {
fChain->push_back(
new ExpMultDist_CCQE_XSec_1DVar_FakeStudy(name, file, rw, type, fkdt));
} else if (name.find("ExpMultDist_CCQE_XSec_2D") != std::string::npos &&
name.find("_FakeStudy") != std::string::npos) {
fChain->push_back(
new ExpMultDist_CCQE_XSec_2DVar_FakeStudy(name, file, rw, type, fkdt));
} else if (name.find("GenericFlux_") != std::string::npos) {
fChain->push_back(new GenericFlux_Tester(name, file, rw, type, fkdt));
} else {
std::cerr << "Error: could not find " << name << std::endl;
exit(-1);
return false;
}
// Return if sample was loaded correctly;
return true;
}
}
diff --git a/src/FCN/SampleList.h b/src/FCN/SampleList.h
index 185ef3f..4dadf04 100644
--- a/src/FCN/SampleList.h
+++ b/src/FCN/SampleList.h
@@ -1,195 +1,197 @@
#ifndef _SAMPLE_LIST_H_
#define _SAMPLE_LIST_H_
/*!
* \addtogroup FCN
* @{
*/
#include "ANL_CCQE_Evt_1DQ2_nu.h"
#include "ANL_CCQE_XSec_1DEnu_nu.h"
// ANL CC1ppip
#include "ANL_CC1ppip_Evt_1DQ2_nu.h"
#include "ANL_CC1ppip_Evt_1DcosmuStar_nu.h"
#include "ANL_CC1ppip_Evt_1DcosmuStar_nu.h"
#include "ANL_CC1ppip_Evt_1DcosthAdler_nu.h"
#include "ANL_CC1ppip_Evt_1Dphi_nu.h"
#include "ANL_CC1ppip_Evt_1Dppi_nu.h"
#include "ANL_CC1ppip_Evt_1Dthpr_nu.h"
#include "ANL_CC1ppip_XSec_1DEnu_nu.h"
#include "ANL_CC1ppip_XSec_1DQ2_nu.h"
// ANL CC1npip
#include "ANL_CC1npip_Evt_1DQ2_nu.h"
#include "ANL_CC1npip_Evt_1DcosmuStar_nu.h"
#include "ANL_CC1npip_Evt_1Dppi_nu.h"
#include "ANL_CC1npip_XSec_1DEnu_nu.h"
// ANL CC1pi0
#include "ANL_CC1pi0_Evt_1DQ2_nu.h"
#include "ANL_CC1pi0_Evt_1DcosmuStar_nu.h"
#include "ANL_CC1pi0_XSec_1DEnu_nu.h"
// ANL NC1npip (mm, exotic!)
#include "ANL_NC1npip_Evt_1Dppi_nu.h"
#include "ArgoNeuT_CCInc_XSec_1Dpmu_antinu.h"
#include "ArgoNeuT_CCInc_XSec_1Dpmu_nu.h"
#include "ArgoNeuT_CCInc_XSec_1Dthetamu_antinu.h"
#include "ArgoNeuT_CCInc_XSec_1Dthetamu_nu.h"
// BNL CCQE
#include "BNL_CCQE_Evt_1DQ2_nu.h"
#include "BNL_CCQE_XSec_1DEnu_nu.h"
// BNL CC1ppip
#include "BNL_CC1ppip_Evt_1DQ2_nu.h"
#include "BNL_CC1ppip_Evt_1DQ2_nu.h"
#include "BNL_CC1ppip_Evt_1DcosthAdler_nu.h"
#include "BNL_CC1ppip_Evt_1Dphi_nu.h"
#include "BNL_CC1ppip_XSec_1DEnu_nu.h"
// BNL CC1npip
#include "BNL_CC1npip_Evt_1DQ2_nu.h"
#include "BNL_CC1npip_XSec_1DEnu_nu.h"
// BNL CC1pi0
#include "BNL_CC1pi0_Evt_1DQ2_nu.h"
#include "BNL_CC1pi0_XSec_1DEnu_nu.h"
// FNAL CCQE
#include "FNAL_CCQE_Evt_1DQ2_nu.h"
// FNAL CC1ppip
#include "FNAL_CC1ppip_Evt_1DQ2_nu.h"
#include "FNAL_CC1ppip_XSec_1DEnu_nu.h"
#include "FNAL_CC1ppip_XSec_1DQ2_nu.h"
// FNAL CC1ppim
//#include "FNAL_CC1ppim_XSec_1DEnu_antinu.h"
// BEBC CCQE
#include "BEBC_CCQE_XSec_1DQ2_nu.h"
// BEBC CC1ppip
#include "BEBC_CC1ppip_XSec_1DEnu_nu.h"
#include "BEBC_CC1ppip_XSec_1DQ2_nu.h"
// BEBC CC1npip
#include "BEBC_CC1npip_XSec_1DEnu_nu.h"
#include "BEBC_CC1npip_XSec_1DQ2_nu.h"
// BEBC CC1pi0
#include "BEBC_CC1pi0_XSec_1DEnu_nu.h"
#include "BEBC_CC1pi0_XSec_1DQ2_nu.h"
// BEBC CC1npim
#include "BEBC_CC1npim_XSec_1DEnu_antinu.h"
#include "BEBC_CC1npim_XSec_1DQ2_antinu.h"
// BEBC CC1ppim
#include "BEBC_CC1ppim_XSec_1DEnu_antinu.h"
#include "BEBC_CC1ppim_XSec_1DQ2_antinu.h"
// GGM CC1ppip
#include "GGM_CC1ppip_Evt_1DQ2_nu.h"
#include "GGM_CC1ppip_XSec_1DEnu_nu.h"
// MiniBooNE CCQE
#include "MiniBooNE_CCQE_XSec_1DQ2_nu.h"
#include "MiniBooNE_CCQE_XSec_1DQ2_antinu.h"
#include "MiniBooNE_CCQE_XSec_2DTcos_antinu.h"
#include "MiniBooNE_CCQE_XSec_2DTcos_antinu.h"
#include "MiniBooNE_CCQE_XSec_2DTcos_nu.h"
// MiniBooNE CC1pi+
#include "MiniBooNE_CC1pip_XSec_1DEnu_nu.h"
#include "MiniBooNE_CC1pip_XSec_1DQ2_nu.h"
#include "MiniBooNE_CC1pip_XSec_1DTpi_nu.h"
#include "MiniBooNE_CC1pip_XSec_1DTu_nu.h"
#include "MiniBooNE_CC1pip_XSec_2DQ2Enu_nu.h"
#include "MiniBooNE_CC1pip_XSec_2DTpiCospi_nu.h"
#include "MiniBooNE_CC1pip_XSec_2DTpiEnu_nu.h"
#include "MiniBooNE_CC1pip_XSec_2DTuCosmu_nu.h"
#include "MiniBooNE_CC1pip_XSec_2DTuEnu_nu.h"
// MiniBooNE CC1pi0
#include "MiniBooNE_CC1pi0_XSec_1DEnu_nu.h"
#include "MiniBooNE_CC1pi0_XSec_1DQ2_nu.h"
#include "MiniBooNE_CC1pi0_XSec_1DTu_nu.h"
#include "MiniBooNE_CC1pi0_XSec_1Dcosmu_nu.h"
#include "MiniBooNE_CC1pi0_XSec_1Dcospi0_nu.h"
#include "MiniBooNE_CC1pi0_XSec_1Dppi0_nu.h"
//#include "MiniBooNE_NCpi0_XSec_1Dppi0_nu.h"
// MiniBooNE NCEL
// #include "MiniBooNE_NCEL_XSec_Treco_nu.h"
// MINERvA CCQE
#include "MINERvA_CCQE_XSec_1DQ2_antinu.h"
#include "MINERvA_CCQE_XSec_1DQ2_joint.h"
#include "MINERvA_CCQE_XSec_1DQ2_nu.h"
// MINERvA CC0pi
#include "MINERvA_CC0pi_XSec_1DEe_nue.h"
#include "MINERvA_CC0pi_XSec_1DQ2_nu_proton.h"
#include "MINERvA_CC0pi_XSec_1DQ2_nue.h"
#include "MINERvA_CC0pi_XSec_1DThetae_nue.h"
// MINERvA CC1pi+
#include "MINERvA_CC1pip_XSec_1DTpi_20deg_nu.h"
#include "MINERvA_CC1pip_XSec_1DTpi_nu.h"
#include "MINERvA_CC1pip_XSec_1Dth_20deg_nu.h"
#include "MINERvA_CC1pip_XSec_1Dth_nu.h"
// MINERvA CCNpi+
#include "MINERvA_CCNpip_XSec_1Dth_nu.h"
#include "MINERvA_CCNpip_XSec_1Dth_20deg_nu.h"
#include "MINERvA_CCNpip_XSec_1DTpi_nu.h"
#include "MINERvA_CCNpip_XSec_1DTpi_20deg_nu.h"
#include "MINERvA_CCNpip_XSec_1Dpmu_nu.h"
#include "MINERvA_CCNpip_XSec_1Dthmu_nu.h"
#include "MINERvA_CCNpip_XSec_1DQ2_nu.h"
#include "MINERvA_CCNpip_XSec_1DEnu_nu.h"
// MINERvA CC1pi0
#include "MINERvA_CC1pi0_XSec_1Dth_antinu.h"
#include "MINERvA_CC1pi0_XSec_1Dppi0_antinu.h"
#include "MINERvA_CC1pi0_XSec_1Dthmu_antinu.h"
#include "MINERvA_CC1pi0_XSec_1Dpmu_antinu.h"
#include "MINERvA_CC1pi0_XSec_1DQ2_antinu.h"
#include "MINERvA_CC1pi0_XSec_1DEnu_antinu.h"
// MINERvA CCINC
#include "MINERvA_CCinc_XSec_2DEavq3_nu.h"
#include "MINERvA_CCinc_XSec_1Dx_ratio.h"
#include "MINERvA_CCinc_XSec_1DEnu_ratio.h"
// T2K CC0pi
#include "T2K_CC0pi_XSec_2DPcos_nu.h"
// T2K CC1pi+ on CH
#include "T2K_CC1pip_CH_XSec_1Dpmu_nu.h"
#include "T2K_CC1pip_CH_XSec_1Dppi_nu.h"
#include "T2K_CC1pip_CH_XSec_1Dthpi_nu.h"
#include "T2K_CC1pip_CH_XSec_1Dthmupi_nu.h"
#include "T2K_CC1pip_CH_XSec_1DQ2_nu.h"
#include "T2K_CC1pip_CH_XSec_1Dq3_nu.h"
#include "T2K_CC1pip_CH_XSec_1Dthq3pi_nu.h"
#include "T2K_CC1pip_CH_XSec_1DWrec_nu.h"
+// T2K STV CC0pi
+#include "T2K_CC0pinp_STV_XSec_1Ddpt_nu.h"
// K2K CC0pi
#include "K2K_CC0pi_XSec_1DDelPhi_nu_Ntrks.h"
#include "K2K_CC0pi_XSec_1DPmu_nu_Ntrks.h"
#include "K2K_CC0pi_XSec_1DQ2_nu_subtrks.h"
#include "K2K_CC0pi_XSec_1DCosThetaMu_nu_subtrks.h"
#include "K2K_CC0pi_XSec_1DThetaMu_nu_Ntrks.h"
// K2K NC1pi0
#include "K2K_NC1pi0_Evt_1Dppi0_nu.h"
#include "ExpMultDist_CCQE_XSec_1DVar_FakeStudy.h"
#include "ExpMultDist_CCQE_XSec_2DVar_FakeStudy.h"
#include "GenericFlux_Tester.h"
#include "FitWeight.h"
#include "parameter_pulls.h"
//! Functions to make it easier for samples to be created and handled.
namespace SampleUtils {
//! Create a given sample given its name, file, type, fakdata(fkdt) file and the
//! current rw engine and push it back into the list fChain.
bool LoadSample(std::list<MeasurementBase*>* fChain, std::string name,
std::string file, std::string type, std::string fkdt,
FitWeight* rw);
}
/*! @} */
#endif
diff --git a/src/FitBase/InputHandler.cxx b/src/FitBase/InputHandler.cxx
index d8f575e..a48099d 100644
--- a/src/FitBase/InputHandler.cxx
+++ b/src/FitBase/InputHandler.cxx
@@ -1,1011 +1,1016 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "InputHandler.h"
//***********************************
InputHandler::InputHandler(std::string handle, std::string infile_name) {
//***********************************
LOG(SAM) << "Creating InputHandler for " << handle << "..." << std::endl;
LOG(SAM) << " -> [" << infile_name << "]" << std::endl;
// Read in parameters for handler
this->maxEvents = FitPar::Config().GetParI("MAXEVENTS");
isJointInput = false;
// Setup a custom Event class
this->cust_event = new FitEvent();
this->signal_event = new BaseFitEvt();
this->inFile = infile_name;
this->handleName = handle;
// Parse Infile to allow enviornmental flags
this->inFile = this->ParseInputFile(this->inFile);
LOG(SAM) << " -> Type = " << inType << std::endl;
LOG(SAM) << " -> Input = " << inFile << std::endl;
// Automatically check what sort of event file it is
if (inType.compare("JOINT"))
this->inRootFile = new TFile(this->inFile.c_str(), "READ");
// Check file exists
if (this->inRootFile->IsZombie()){
ERR(FTL) << "Cannot find InputFile!" << endl;
throw;
}
-
+
// Setup the handler for each type
if (!inType.compare("NEUT"))
this->ReadNeutFile();
else if (!inType.compare("NUWRO"))
this->ReadNuWroFile();
else if (!inType.compare("GENIE"))
this->ReadGenieFile();
else if (!inType.compare("GiBUU_nu"))
this->ReadGiBUUFile(false);
else if (!inType.compare("GiBUU_nub"))
this->ReadGiBUUFile(true);
else if (!inType.compare("HIST"))
this->ReadHistogramFile();
else if (!inType.compare("BNSPLN"))
this->ReadBinSplineFile();
else if (!inType.compare("EVSPLN"))
this->ReadEventSplineFile();
else if (!inType.compare("NUANCE"))
this->ReadNuanceFile();
else if (!inType.compare("JOINT"))
this->ReadJointFile();
else {
LOG(FTL) << " -> ERROR: Invalid Event File Type" << std::endl;
inRootFile->ls();
exit(-1);
}
// Setup MaxEvents After setup of ttree
if (maxEvents > 1 and maxEvents < nEvents) {
LOG(SAM) << " -> Reading only " << maxEvents << " events from total."
<< std::endl;
nEvents = maxEvents;
}
this->fluxList.push_back(this->fluxHist);
this->eventList.push_back(this->eventHist);
this->xsecList.push_back(this->xsecHist);
LOG(SAM) << " -> Finished handler initialisation." << std::endl;
return;
};
//********************************************************************
std::string InputHandler::ParseInputFile(std::string inputstring) {
//********************************************************************
// Parse out the input_type
const int nfiletypes = 8;
const std::string filetypes[nfiletypes] = {"NEUT", "NUWRO", "GENIE",
"EVSPLN", "JOINT", "NUANCE",
"GiBUU_nu", "GiBUU_nub"};
for (int i = 0; i < nfiletypes; i++) {
std::string tempTypes = filetypes[i] + ":";
if (inputstring.find(tempTypes) != std::string::npos) {
inType = filetypes[i];
inputstring.replace(inputstring.find(tempTypes), tempTypes.size(), "");
break;
}
}
// If no input type ERROR!
if (inType.empty()){
ERR(FTL) << "No input type supplied for InputHandler!" << endl;
ERR(FTL) << "Problematic Input: " << inputstring << endl;
throw;
}
-
+
// Parse out envir flags
const int nfiledir = 5;
const std::string filedir[nfiledir] = {"NEUT_DIR", "NUWRO_DIR", "GENIE_DIR",
"NUANCE_DIR", "EVSPLN_DIR"};
for (int i = 0; i < nfiledir; i++) {
std::string tempDir = "@" + filedir[i];
if (inputstring.find(tempDir) != std::string::npos) {
std::string event_folder = FitPar::Config().GetParS(filedir[i]);
inputstring.replace(inputstring.find(tempDir), tempDir.size(),
event_folder);
break;
}
}
return inputstring;
}
//********************************************************************
bool InputHandler::CanIGoFast() {
//********************************************************************
if (eventType == 6) {
return true;
}
return false;
}
//********************************************************************
void InputHandler::ReadEventSplineFile() {
//********************************************************************
LOG(SAM) << " -> Setting up SPLINE inputs" << std::endl;
// Event Type 7 SPLINES
this->eventType = 6;
// Get flux histograms NEUT supplies
this->fluxHist = (TH1D*)inRootFile->Get((this->handleName + "_FLUX").c_str());
this->eventHist = (TH1D*)inRootFile->Get((this->handleName + "_EVT").c_str());
this->xsecHist = (TH1D*)inRootFile->Get((this->handleName + "_XSEC").c_str());
// Setup Spline Stuff
this->splhead = (FitSplineHead*)inRootFile->Get(
(this->handleName + "_splineHead").c_str());
tn = new TChain(Form("%s", (this->handleName + "_splineEvents").c_str()), "");
tn->Add(Form("%s/%s", this->inFile.c_str(),
(this->handleName + "_splineEvents").c_str()));
// Assign nvect
nEvents = tn->GetEntries();
cust_event = NULL;
tn->SetBranchAddress("FitEvent", &cust_event);
// Load Dial Coeffs into vector
for (int i = 0; i < nEvents; i++) {
tn->GetEntry(i);
tn->Show(i);
spline_list.push_back(*cust_event->dial_coeff);
}
sleep(5);
// Set MAXEVENTS CALC Here before we load in splines
if (maxEvents > 1 and maxEvents < nEvents) {
LOG(SAM) << " -> Reading only " << maxEvents
<< " events from total spline events." << std::endl;
nEvents = maxEvents;
}
// Load all the splines into signal memory
// for (int i = 0; i < nEvents; i++){
// tn->GetEntry(i);
// BaseFitEvt* base_event = (new BaseFitEvt(cust_event));
// base_event->fType=6;
// signal_events.push_back( base_event );
// }
// Print out what was read in
LOG(SAM) << " -> Successfully Read SPLINE file" << std::endl;
if (LOG_LEVEL(SAM)) this->PrintStartInput();
int cnt = 1;
std::list<FitSpline*>::iterator spl_iter =
this->splhead->SplineObjects.begin();
for (; spl_iter != this->splhead->SplineObjects.end(); spl_iter++) {
FitSpline* spl = (*spl_iter);
LOG(SAM) << " -> Spline " << cnt << ". " << spl->id << " " << spl->form
<< " "
<< "NDIM(" << spl->ndim << ") "
<< "NPAR(" << spl->npar << ") "
<< "PTS(" << spl->points << ") " << std::endl;
cnt++;
}
}
//********************************************************************
FitSplineHead* InputHandler::GetSplineHead() {
//********************************************************************
return this->splhead;
}
//********************************************************************
void InputHandler::ReadJointFile() {
//********************************************************************
LOG(SAM) << " -> Reading list of inputs from file" << std::endl;
isJointInput = true;
// Parse Input File
std::string line;
std::ifstream card(inFile.c_str(), ifstream::in);
std::vector<std::string> input_lines;
while (std::getline(card, line, '\n')) {
std::istringstream stream(line);
// Add normalisation option for second line
input_lines.push_back(line);
// Split to get normalisation
}
card.close();
// Loop over input and get the flux files
// Using a temporary input handler to do this, which is a bit dodge.
int count_low = 0;
int temp_type = -1;
for (UInt_t i = 0; i < input_lines.size(); i++) {
// Create Temporary InputHandlers inside
InputHandler* temp_input = new InputHandler(
std::string(Form("temp_input_%i", i)), input_lines.at(i));
if (temp_type != temp_input->GetType() and i > 0) {
ERR(FTL) << " Can't use joint events with mismatched trees yet!"
<< std::endl;
ERR(FTL) << " Make them all the same type!" << std::endl;
}
temp_type = temp_input->GetType();
TH1D* temp_flux = (TH1D*)temp_input->GetFluxHistogram()->Clone();
TH1D* temp_evts = (TH1D*)temp_input->GetEventHistogram()->Clone();
TH1D* temp_xsec = (TH1D*)temp_input->GetXSecHistogram()->Clone();
int temp_events = temp_input->GetNEvents();
temp_flux->SetName(
(this->handleName + "_" + temp_input->GetInputStateString() + "_FLUX")
.c_str());
temp_evts->SetName(
(this->handleName + "_" + temp_input->GetInputStateString() + "_EVT")
.c_str());
temp_xsec->SetName(
(this->handleName + "_" + temp_input->GetInputStateString() + "_XSEC")
.c_str());
this->fluxList.push_back(temp_flux);
this->eventList.push_back(temp_evts);
this->xsecList.push_back(temp_xsec);
this->joint_index_low.push_back(count_low);
this->joint_index_high.push_back(count_low + temp_events);
this->joint_index_hist.push_back((TH1D*)temp_evts->Clone());
count_low += temp_events;
if (i == 0) {
this->fluxHist = (TH1D*)temp_flux->Clone();
this->eventHist = (TH1D*)temp_evts->Clone();
} else {
this->fluxHist->Add(temp_flux);
this->eventHist->Add(temp_evts);
}
std::cout << "Added Input File " << input_lines.at(i) << std::endl
<< " with " << temp_events << std::endl;
}
// Now have all correctly normalised histograms all we need to do is setup the
// TChains
// Input Assumes all the same type
std::string tree_name = "";
if (temp_type == 0)
tree_name = "neuttree";
else if (temp_type == 1)
tree_name = "treeout";
// Add up the TChains
tn = new TChain(tree_name.c_str());
for (UInt_t i = 0; i < input_lines.size(); i++) {
// PARSE INPUT
std::cout << "Adding new tchain " << input_lines.at(i) << std::endl;
std::string temp_file = this->ParseInputFile(input_lines.at(i));
tn->Add(temp_file.c_str());
}
// Setup Events
nEvents = tn->GetEntries();
if (temp_type == 0) {
#ifdef __NEUT_ENABLED__
eventType = 0;
neut_event = NULL;
tn->SetBranchAddress("vectorbranch", &neut_event);
this->cust_event->SetEventAddress(&neut_event);
#endif
} else if (temp_type == 1) {
#ifdef __NUWRO_ENABLED__
eventType = 1;
nuwro_event = NULL;
tn->SetBranchAddress("e", &nuwro_event);
this->cust_event->SetEventAddress(&nuwro_event);
#endif
}
// Normalise event histogram PDFS for weights
for (UInt_t i = 0; i < input_lines.size(); i++) {
TH1D* temp_hist = (TH1D*)joint_index_hist.at(i)->Clone();
joint_index_weight.push_back(
double(nEvents) / eventHist->Integral("width") *
joint_index_hist.at(i)->Integral("width") /
double(joint_index_high.at(i) - joint_index_low.at(i)));
temp_hist->Scale(double(nEvents) / eventHist->Integral("width"));
temp_hist->Scale(joint_index_hist.at(i)->Integral("width") /
double(joint_index_high.at(i)));
this->joint_index_hist.at(i) = temp_hist;
}
this->eventHist->SetNameTitle((this->handleName + "_EVT").c_str(),
(this->handleName + "_EVT").c_str());
this->fluxHist->SetNameTitle((this->handleName + "_FLUX").c_str(),
(this->handleName + "_FLUX").c_str());
return;
}
//********************************************************************
void InputHandler::ReadNeutFile() {
//********************************************************************
#ifdef __NEUT_ENABLED__
LOG(SAM) << " -> Setting up NEUT inputs" << std::endl;
// Event Type 0 Neut
this->eventType = 0;
// Get flux histograms NEUT supplies
this->fluxHist = (TH1D*)inRootFile->Get(
(PlotUtils::GetObjectWithName(inRootFile, "flux")).c_str());
this->fluxHist->SetNameTitle((this->handleName + "_FLUX").c_str(),
(this->handleName + "; E_{#nu} (GeV)").c_str());
this->eventHist = (TH1D*)inRootFile->Get(
(PlotUtils::GetObjectWithName(inRootFile, "evtrt")).c_str());
this->eventHist->SetNameTitle(
(this->handleName + "_EVT").c_str(),
(this->handleName + "; E_{#nu} (GeV); Event Rate").c_str());
this->xsecHist = (TH1D*)eventHist->Clone();
this->xsecHist->Divide(this->fluxHist);
this->xsecHist->SetNameTitle(
(this->handleName + "_XSEC").c_str(),
(this->handleName + "_XSEC;E_{#nu} (GeV); XSec (1#times10^{-38} cm^{2})")
.c_str());
// Read in the file once only
tn = new TChain("neuttree", "");
tn->Add(Form("%s/neuttree", this->inFile.c_str()));
// Assign nvect
nEvents = tn->GetEntries();
neut_event = NULL;
tn->SetBranchAddress("vectorbranch", &neut_event);
// Make the custom event read in nvect when calling CalcKinematics
this->cust_event->SetEventAddress(&neut_event);
// Print out what was read in
LOG(SAM) << " -> Successfully Read NEUT file" << std::endl;
if (LOG_LEVEL(SAM)) this->PrintStartInput();
#else
ERR(FTL) << "ERROR: Invalid Event File Provided" << std::endl;
ERR(FTL) << "NEUT Input Not Enabled." << std::endl;
ERR(FTL) << "Rebuild with --enable-neut or check FitBuild.h!" << std::endl;
exit(-1);
#endif
return;
}
//********************************************************************
void InputHandler::ReadNuWroFile() {
//********************************************************************
#ifdef __NUWRO_ENABLED__
LOG(SAM) << " -> Setting up Nuwro inputs" << std::endl;
// Event Type 1 == NuWro
this->eventType = 1;
// Setup the TChain for nuwro event tree
tn = new TChain("treeout");
tn->AddFile(this->inFile.c_str());
// Get entries and nuwro_event
nEvents = tn->GetEntries();
nuwro_event = NULL;
tn->SetBranchAddress("e", &nuwro_event);
this->cust_event->SetEventAddress(&nuwro_event);
// Check if we have saved an xsec histogram before
this->fluxHist = (TH1D*)inRootFile->Get(
(PlotUtils::GetObjectWithName(inRootFile, "FluxHist")).c_str());
this->eventHist = (TH1D*)inRootFile->Get(
(PlotUtils::GetObjectWithName(inRootFile, "EvtHist")).c_str());
// Check if we are forcing plot generation (takes time)
bool regenFlux = FitPar::Config().GetParB("input.regen_nuwro_plots");
if (regenFlux)
LOG(SAM)
<< " -> Forcing NuWro XSec/Flux plots to be generated at the start. "
<< std::endl;
// Already generated flux and event histograms
if (fluxHist and eventHist and !regenFlux) {
this->xsecHist = (TH1D*)inRootFile->Get(
(PlotUtils::GetObjectWithName(inRootFile, "xsec")).c_str());
this->fluxHist->SetNameTitle((this->handleName + "_FLUX").c_str(),
(this->handleName + "_FLUX").c_str());
this->eventHist->SetNameTitle((this->handleName + "_EVT").c_str(),
(this->handleName + "_EVT").c_str());
this->xsecHist->SetNameTitle((this->handleName + "_XSEC").c_str(),
(this->handleName + "_XSEC").c_str());
// Need to regenerate if not found
} else {
LOG(SAM)
<< " -> No NuWro XSec or Flux Histograms found, need to regenerate!"
<< std::endl;
// Can grab flux histogram from the pars
tn->GetEntry(0);
int beamtype = nuwro_event->par.beam_type;
if (beamtype == 0) {
std::string fluxstring = nuwro_event->par.beam_energy;
std::vector<double> fluxvals =
PlotUtils::FillVectorDFromString(fluxstring, " ");
int pdg = nuwro_event->par.beam_particle;
double Elow = double(fluxvals[0]) / 1000.0;
double Ehigh = double(fluxvals[1]) / 1000.0;
std::cout << " - Adding new nuwro flux "
<< "pdg: " << pdg << "Elow: " << Elow << "Ehigh: " << Ehigh
<< std::endl;
fluxHist = new TH1D("fluxplot", "fluxplot", fluxvals.size() - 2, Elow, Ehigh);
for (int j = 2; j < fluxvals.size(); j++) {
cout << j << " " << fluxvals[j] << endl;
fluxHist->SetBinContent(j - 1, fluxvals[j]);
}
} else if (beamtype == 1) {
std::string fluxstring = nuwro_event->par.beam_content;
std::vector<std::string> fluxlines =
PlotUtils::FillVectorSFromString(fluxstring, "\n");
for (int i = 0; i < fluxlines.size(); i++) {
std::vector<double> fluxvals =
PlotUtils::FillVectorDFromString(fluxlines[i], " ");
int pdg = int(fluxvals[0]);
double pctg = double(fluxvals[1]) / 100.0;
double Elow = double(fluxvals[2]) / 1000.0;
double Ehigh = double(fluxvals[3]) / 1000.0;
std::cout << " - Adding new nuwro flux "
<< "pdg: " << pdg << "pctg: " << pctg << "Elow: " << Elow
<< "Ehigh: " << Ehigh << std::endl;
TH1D* fluxplot =
new TH1D("fluxplot", "fluxplot", fluxvals.size() - 4, Elow, Ehigh);
for (int j = 4; j < fluxvals.size(); j++) {
fluxplot->SetBinContent(j + 1, fluxvals[j]);
}
if (this->fluxHist)
fluxHist->Add(fluxplot);
else
this->fluxHist = (TH1D*)fluxplot->Clone();
}
}
this->fluxHist->SetNameTitle("nuwro_flux",
"nuwro_flux;E_{#nu} (GeV); Flux");
this->eventHist = (TH1D*)this->fluxHist->Clone();
this->eventHist->Reset();
this->eventHist->SetNameTitle("nuwro_evt", "nuwro_evt");
this->xsecHist = (TH1D*)this->fluxHist->Clone();
this->xsecHist->Reset();
this->xsecHist->SetNameTitle("nuwro_xsec", "nuwro_xsec");
// Start Processing
LOG(SAM) << " -> Processing NuWro Input Flux for " << nEvents
<< " events (This can take a while...) " << std::endl;
double Enu = 0.0;
double TotXSec = 0.0;
double totaleventmode = 0.0;
double totalevents = 0.0;
// --- loop
for (int i = 0; i < nEvents; i++) {
tn->GetEntry(i);
if (i % 100000 == 0) cout << " i " << i << std::endl;
// Get Variables
Enu = nuwro_event->in[0].E() / 1000.0;
TotXSec = nuwro_event->weight;
// Fill a flux and xsec histogram
this->eventHist->Fill(Enu);
this->xsecHist->Fill(Enu, TotXSec);
// Keep Tally
totaleventmode += TotXSec;
totalevents++;
};
LOG(SAM) << " -> Flux Processing Loop Finished." << std::endl;
if (this->eventHist->Integral() == 0.0) {
std::cout << "ERROR NO EVENTS FOUND IN RANGE! " << std::endl;
exit(-1);
}
// Sort out plot scaling
double AvgXSec = (totaleventmode * 1.0E38 / (totalevents + 0.));
LOG(SAM) << " -> Average XSec = " << AvgXSec << std::endl;
this->eventHist->Scale(1.0 / eventHist->Integral()); // Convert to PDF
this->eventHist->Scale(this->fluxHist->Integral() *
AvgXSec); // Convert to Proper Event Rate
this->xsecHist->Add(eventHist); // Get Event Rate Plot
this->xsecHist->Divide(this->fluxHist); // Make XSec Plot
// this->eventHist = (TH1D*)this->fluxHist->Clone();
// this->eventHist->Multiply(this->xsecHist);
// Clear over/underflows incase they mess with integrals later.
this->fluxHist->SetBinContent(0, 0.0);
this->fluxHist->SetBinContent(this->fluxHist->GetNbinsX() + 2, 0.0);
this->eventHist->SetBinContent(0, 0.0);
this->eventHist->SetBinContent(this->eventHist->GetNbinsX() + 2, 0.0);
LOG(SAM)
<< " -> Finished making NuWro event plots. Saving them for next time..."
<< std::endl;
TFile* temp_save_file = new TFile(this->inFile.c_str(), "UPDATE");
temp_save_file->cd();
this->fluxHist->Write("nuwro_flux", TObject::kOverwrite);
this->eventHist->Write("nuwro_evtrt", TObject::kOverwrite);
this->xsecHist->Write("nuwro_xsec", TObject::kOverwrite);
temp_save_file->ls();
temp_save_file->Close();
delete temp_save_file;
this->fluxHist->SetNameTitle((this->handleName + "_FLUX").c_str(),
(this->handleName + "_FLUX").c_str());
this->eventHist->SetNameTitle((this->handleName + "_EVT").c_str(),
(this->handleName + "_EVT").c_str());
this->xsecHist->SetNameTitle((this->handleName + "_XSEC").c_str(),
(this->handleName + "_XSEC").c_str());
}
// Print out what was read in
LOG(SAM) << " -> Successfully Read NUWRO file" << std::endl;
if (LOG_LEVEL(SAM)) this->PrintStartInput();
#else
ERR(FTL) << "ERROR: Invalid Event File Provided" << std::endl;
ERR(FTL) << "NuWro Input Not Enabled." << std::endl;
ERR(FTL) << "Rebuild with --enable-nuwro or check FitBuild.h!" << std::endl;
exit(-1);
#endif
return;
}
//********************************************************************
void InputHandler::ReadGenieFile() {
//********************************************************************
#ifdef __GENIE_ENABLED__
// Event Type 1 NuWro
this->eventType = 5;
// Open Root File
LOG(SAM) << "Reading event file " << this->inFile << std::endl;
// Get flux histograms NEUT supplies
this->fluxHist = (TH1D*)inRootFile->Get(
(PlotUtils::GetObjectWithName(inRootFile, "spectrum")).c_str());
this->fluxHist->SetNameTitle((this->handleName + "_FLUX").c_str(),
(this->handleName + "; E_{#nu} (GeV)").c_str());
this->eventHist = (TH1D*)inRootFile->Get(
(PlotUtils::GetObjectWithName(inRootFile, "spectrum")).c_str());
this->eventHist->SetNameTitle(
(this->handleName + "_EVT").c_str(),
(this->handleName + "; E_{#nu} (GeV); Event Rate").c_str());
this->xsecHist = (TH1D*)inRootFile->Get(
(PlotUtils::GetObjectWithName(inRootFile, "spectrum")).c_str());
this->xsecHist->SetNameTitle(
(this->handleName + "_XSEC").c_str(),
(this->handleName + "; E_{#nu} (GeV); Event Rate").c_str());
double average_xsec = 0.0;
int total_events = 0;
// Setup the TChain for nuwro event tree
tn = new TChain("gtree");
tn->AddFile(this->inFile.c_str());
nEvents = tn->GetEntries();
LOG(SAM) << "Number of GENIE Eevents " << tn->GetEntries() << std::endl;
genie_event = NULL;
mcrec = NULL;
// NtpMCEventRecord * mcrec = 0; tree->SetBranchAddress(gmrec, &mcrec);
tn->SetBranchAddress("gmcrec", &mcrec);
this->eventHist->Reset();
// Make the custom event read in nvect when calling CalcKinematics
this->cust_event->SetEventAddress(&mcrec);
LOG(SAM) << "Processing GENIE flux events." << std::endl;
for (int i = 0; i < nEvents; i++) {
tn->GetEntry(i);
EventRecord& event = *(mcrec->event);
GHepParticle* neu = event.Probe();
GHepRecord genie_record = static_cast<GHepRecord>(event);
double xsec = (genie_record.XSec() / (1E-38 * genie::units::cm2));
average_xsec += xsec;
total_events += 1;
this->eventHist->Fill(neu->E());
this->xsecHist->Fill(neu->E(), xsec);
mcrec->Clear();
}
// Sort xsec hist
this->xsecHist->Divide(this->eventHist);
// Sort eventHist as xsecHist * eventHist
this->eventHist = (TH1D*)this->fluxHist->Clone();
this->eventHist->Multiply(this->xsecHist);
this->eventHist->SetNameTitle(
(this->handleName + "_EVT").c_str(),
(this->handleName + "_EVT;E_{#nu} (GeV); Events (1#times10^{-38})")
.c_str());
this->xsecHist->SetNameTitle(
(this->handleName + "_XSEC").c_str(),
(this->handleName + "_XSEC;E_{#nu} (GeV); XSec (1#times10^{-38} cm^{2})")
.c_str());
#else
ERR(FTL) << "ERROR: Invalid Event File Provided" << std::endl;
ERR(FTL) << "GENIE Input Not Enabled." << std::endl;
ERR(FTL) << "Rebuild with --enable-genie or check FitBuild.h!" << std::endl;
exit(-1);
#endif
return;
}
//********************************************************************
void InputHandler::ReadGiBUUFile(bool IsNuBarDominant) {
//********************************************************************
#ifdef __GiBUU_ENABLED__
this->eventType = kGiBUU;
// Open Root File
LOG(SAM) << "Opening event file " << this->inFile << std::endl;
TFile* rootFile = new TFile(this->inFile.c_str(), "READ");
// Get flux histograms NEUT supplies
TH1D* numuFlux = dynamic_cast<TH1D*>(rootFile->Get("numu_flux"));
TH1D* numubFlux = dynamic_cast<TH1D*>(rootFile->Get("numub_flux"));
if (numuFlux) {
numuFlux = static_cast<TH1D*>(numuFlux->Clone());
numuFlux->SetDirectory(NULL);
numuFlux->SetNameTitle(
(this->handleName + "_numu_FLUX").c_str(),
(this->handleName + "; E_{#nu} (GeV); #Phi_{#nu} (A.U.)").c_str());
fluxList.push_back(numuFlux);
}
if (numubFlux) {
numubFlux = static_cast<TH1D*>(numubFlux->Clone());
numubFlux->SetDirectory(NULL);
numubFlux->SetNameTitle(
(this->handleName + "_numub_FLUX").c_str(),
(this->handleName + "; E_{#nu} (GeV); #Phi_{#bar{#nu}} (A.U.)")
.c_str());
fluxList.push_back(numubFlux);
}
rootFile->Close();
// Set flux hist to the dominant mode
fluxHist = IsNuBarDominant ? numubFlux : numuFlux;
if (!fluxHist) {
ERR(FTL) << "Couldn't find: "
<< (IsNuBarDominant ? "numub_flux" : "numu_flux")
<< " in input file: " << inRootFile->GetName() << std::endl;
exit(1);
}
fluxHist->SetNameTitle(
(this->handleName + "_FLUX").c_str(),
(this->handleName + "; E_{#nu} (GeV);" +
(IsNuBarDominant ? "#Phi_{#bar{#nu}} (A.U.)" : "#Phi_{#nu} (A.U.)"))
.c_str());
tn = new TChain("giRooTracker");
tn->AddFile(this->inFile.c_str());
eventHist =
static_cast<TH1D*>(fluxHist->Clone((this->handleName + "_EVT").c_str()));
eventHist->Reset();
nEvents = tn->GetEntries();
eventHist->SetBinContent(
1, double(nEvents) / eventHist->GetXaxis()->GetBinWidth(1));
GiBUUStdHepReader* giRead = new GiBUUStdHepReader();
giRead->SetBranchAddresses(tn);
cust_event->SetEventAddress(giRead);
+ #else
+ ERR(FTL) << "ERROR: Invalid Event File Provided" << std::endl;
+ ERR(FTL) << "GiBUU Input Not Enabled." << std::endl;
+ ERR(FTL) << "Rebuild with -DUSE_GiBUU=1." << std::endl;
+ exit(-1);
#endif
}
//********************************************************************
void InputHandler::ReadBinSplineFile() {
//********************************************************************
// Bin Splines are saved as one event for each histogram bin.
// So just read in as normal event splines and it'll all get sorted easily.
}
//********************************************************************
void InputHandler::ReadHistogramFile() {
//********************************************************************
// Convert the raw histogram into a series of events with X variables
// So we don't have to pass stuff upsteam
}
//********************************************************************
void InputHandler::ReadNuanceFile() {
//********************************************************************
#ifdef __NUANCE_ENABLED__
// Read in Nuance output ROOT file (converted from hbook)
LOG(SAM) << " Reading NUANCE " << std::endl;
eventType = kNUANCE;
// Read in NUANCE Tree
tn = new TChain("h3");
tn->AddFile(this->inFile.c_str());
// Get entries and nuwro_event
nEvents = tn->GetEntries();
nuance_event = new NuanceEvent();
// SetBranchAddress for Nuance
// tn->SetBranchAddress("cc",&nuance_event->cc);
// tn->SetBranchAddress("bound",&nuance_event->bound);
tn->SetBranchAddress("neutrino", &nuance_event->neutrino);
tn->SetBranchAddress("target", &nuance_event->target);
tn->SetBranchAddress("channel", &nuance_event->channel);
// tn->SetBranchAddress("iniQ", &nuance_event->iniQ);
// tn->SetBranchAddress("finQ", &nuance_event->finQ);
// tn->SetBranchAddress("lepton0", &nuance_event->lepton0);
// tn->SetBranchAddress("polar", &nuance_event->polar);
// tn->SetBranchAddress("qsq", &nuance_event->qsq);
// tn->SetBranchAddress("w", &nuance_event->w);
// tn->SetBranchAddress("x",&nuance_event->x);
// tn->SetBranchAddress("y",&nuance_event->y);
tn->SetBranchAddress("p_neutrino", &nuance_event->p_neutrino);
tn->SetBranchAddress("p_targ", &nuance_event->p_targ);
// tn->SetBranchAddress("vertex", &nuance_event->vertex);
// tn->SetBranchAddress("start",&nuance_event->start);
// tn->SetBranchAddress("depth",&nuance_event->depth);
// tn->SetBranchAddress("flux",&nuance_event->flux);
tn->SetBranchAddress("n_leptons", &nuance_event->n_leptons);
tn->SetBranchAddress("p_ltot", &nuance_event->p_ltot);
tn->SetBranchAddress("lepton", &nuance_event->lepton);
tn->SetBranchAddress("p_lepton", &nuance_event->p_lepton);
tn->SetBranchAddress("n_hadrons", &nuance_event->n_hadrons);
tn->SetBranchAddress("p_htot", &nuance_event->p_htot);
tn->SetBranchAddress("hadron", &nuance_event->hadron);
tn->SetBranchAddress("p_hadron", &nuance_event->p_hadron);
this->cust_event->SetEventAddress(&nuance_event);
this->fluxHist = new TH1D((this->handleName + "_FLUX").c_str(),
(this->handleName + "_FLUX").c_str(), 1, 0.0, 1.0);
this->fluxHist->SetBinContent(1, 1.0);
this->eventHist = new TH1D((this->handleName + "_EVT").c_str(),
(this->handleName + "_EVT").c_str(), 1, 0.0, 1.0);
this->eventHist->SetBinContent(1, nEvents);
#else
ERR(FTL) << "ERROR: Invalid Event File Provided" << std::endl;
ERR(FTL) << "NUANCE Input Not Enabled." << std::endl;
ERR(FTL) << "Rebuild with -DUSE_NUANCE=1!" << std::endl;
exit(-1);
#endif
}
//********************************************************************
void InputHandler::PrintStartInput() {
//********************************************************************
LOG(SAM) << " -> Total events = " << nEvents << std::endl;
LOG(SAM) << " -> Energy Range = " << fluxHist->GetXaxis()->GetXmin() << "-"
<< fluxHist->GetXaxis()->GetXmax() << " GeV" << std::endl;
LOG(SAM) << " -> Integrated Flux Hist = "
<< fluxHist->Integral(0, fluxHist->GetNbinsX(), "width")
<< std::endl;
LOG(SAM) << " -> Integrated Event Hist = "
<< eventHist->Integral(0, eventHist->GetNbinsX(), "width")
<< std::endl;
LOG(SAM) << " -> Integrated Inclusive XSec = "
<< (eventHist->Integral(0, eventHist->GetNbinsX(), "width") /
fluxHist->Integral(0, fluxHist->GetNbinsX(), "width")) *
1E-38
<< std::endl;
if (eventType == kEVTSPLINE) return;
// Get First event info
tn->GetEntry(0);
cust_event->CalcKinematics();
LOG(SAM) << " -> Event 0. Neutrino PDG = " << cust_event->PartInfo(0)->fPID
<< std::endl;
LOG(SAM) << " Target A = " << cust_event->TargetA
<< std::endl;
LOG(SAM) << " Target Z = " << cust_event->TargetZ
<< std::endl;
}
//********************************************************************
std::string InputHandler::GetInputStateString() {
//********************************************************************
-
+
tn->GetEntry(0);
cust_event->CalcKinematics();
std::ostringstream state;
state << "T" << eventType << "_PDG" << cust_event->PartInfo(0)->fPID << "_Z"
<< cust_event->TargetZ << "_A" << cust_event->TargetA;
return state.str();
}
//********************************************************************
void InputHandler::ReadEvent(unsigned int i) {
//********************************************************************
bool using_events =
(eventType == 0 or eventType == 5 or eventType == 1 or
eventType == kEVTSPLINE or eventType == kNUANCE or eventType == kGiBUU);
if (using_events) {
tn->GetEntry(i);
if (eventType != kEVTSPLINE) cust_event->CalcKinematics();
cust_event->Index = i;
cur_entry = i;
cust_event->InputWeight = GetInputWeight(i);
} else {
this->GetTreeEntry(i);
}
}
//********************************************************************
void InputHandler::GetTreeEntry(const Long64_t i) {
//********************************************************************
if (eventType != kEVTSPLINE)
tn->GetEntry(i);
else
(*(cust_event->dial_coeff)) = spline_list.at(i);
cur_entry = i;
cust_event->InputWeight = GetInputWeight(i);
}
//********************************************************************
double InputHandler::GetInputWeight(const int entry) {
//********************************************************************
if (eventType == kGiBUU) {
return cust_event->InputWeight;
}
if (!isJointInput) {
return 1.0;
}
double weight = 1.0;
// Find Histogram
for (UInt_t j = 0; j < joint_index_low.size(); j++) {
if (entry >= joint_index_low.at(j) and entry < joint_index_high.at(j)) {
weight *= joint_index_weight.at(j);
break;
}
}
return weight;
}
//********************************************************************
int InputHandler::GetGenEvents() {
//********************************************************************
if (eventType == 6)
return this->splhead->ngen_events;
else
return this->GetNEvents();
}
//********************************************************************
double InputHandler::TotalIntegratedFlux(double low, double high,
std::string intOpt) {
//********************************************************************
throw;
int minBin = this->fluxHist->GetXaxis()->FindBin(low);
int maxBin = this->fluxHist->GetXaxis()->FindBin(high);
double integral =
this->fluxHist->Integral(minBin, maxBin + 1, intOpt.c_str());
return integral;
};
//********************************************************************
double InputHandler::PredictedEventRate(double low, double high,
std::string intOpt) {
//********************************************************************
int minBin = this->fluxHist->GetXaxis()->FindBin(low);
int maxBin = this->fluxHist->GetXaxis()->FindBin(high);
return this->eventHist->Integral(minBin, maxBin + 1, intOpt.c_str());
}
diff --git a/src/FitBase/Measurement1D.cxx b/src/FitBase/Measurement1D.cxx
index 8ef36f1..3cbc19b 100644
--- a/src/FitBase/Measurement1D.cxx
+++ b/src/FitBase/Measurement1D.cxx
@@ -1,1202 +1,1202 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "Measurement1D.h"
/*
Constructor/destructor Functions
*/
//********************************************************************
Measurement1D::Measurement1D() {
//********************************************************************
this->currentNorm = 1.0;
this->mcHist = NULL;
this->dataHist = NULL;
this->mcFine = NULL;
this->maskHist = NULL;
this->covar = NULL;
this->fullcovar = NULL;
this->decomp = NULL;
this->fakeDataFile = "";
this->fluxHist = NULL;
this->eventHist = NULL;
this->xsecHist = NULL;
default_types = "FIX/FULL/CHI2";
allowed_types = "FIX,FREE,SHAPE/FULL,DIAG/CHI2/NORM/ENUCORR/Q2CORR/ENU1D/MASK";
isFix = false;
isShape = false;
isFree = false;
isDiag = false;
isFull = false;
addNormPenalty = false;
isMask = false;
isChi2SVD = false;
isRawEvents = false;
isDifXSec = false;
isEnu1D = false;
}
//********************************************************************
Measurement1D::~Measurement1D(){
//********************************************************************
}
//********************************************************************
void Measurement1D::Init(){
//********************************************************************
}
/*
Setup Functions
-- All these are used only at the start of the Measurement
*/
//********************************************************************
void Measurement1D::SetupMeasurement(std::string inputfile, std::string type, FitWeight *rw, std::string fkdt){
//********************************************************************
// Reset everything to NULL
Init();
// Check if name contains Evt, indicating that it is a raw number of events measurements
// and should thus be treated as once
isRawEvents = false;
if ((measurementName.find("Evt") != std::string::npos) && isRawEvents == false) {
isRawEvents = true;
LOG(SAM) << "Found event rate measurement but isRawEvents == false!" << std::endl;
LOG(SAM) << "Overriding this and setting isRawEvents == true!" << std::endl;
}
isEnu1D = false;
if (measurementName.find("XSec_1DEnu") != std::string::npos) {
isEnu1D = true;
LOG(SAM) << "::" << measurementName << "::" << std::endl;
LOG(SAM) << "Found XSec Enu measurement, applying flux integrated scaling, not flux averaged!" << std::endl;
}
if (isEnu1D && isRawEvents) {
LOG(SAM) << "Found 1D Enu XSec distribution AND isRawEvents, is this really correct?!" << std::endl;
LOG(SAM) << "Check experiment constructor for " << measurementName << " and correct this!" << std::endl;
LOG(SAM) << "I live in " << __FILE__ << ":" << __LINE__ << std::endl;
exit(-1);
}
rw_engine = rw;
this->SetupInputs(inputfile);
// Set Default Options
SetFitOptions( this->default_types );
// Set Passed Options
SetFitOptions(type);
// Still adding support for flat flux inputs
// // Set Enu Flux Scaling
// if (isFlatFluxFolding) this->Input()->ApplyFluxFolding( this->defaultFluxHist );
}
//********************************************************************
void Measurement1D::SetupDefaultHist(){
//********************************************************************
// Setup mcHist
this->mcHist = (TH1D*) this->dataHist->Clone();
this->mcHist->SetNameTitle( (this->measurementName + "_MC").c_str(), (this->measurementName + "_MC" + this->plotTitles).c_str() );
// Setup mcFine
Int_t nBins = this->mcHist->GetNbinsX();
this->mcFine = new TH1D( (this->measurementName + "_MC_FINE").c_str(), (this->measurementName + "_MC_FINE" + this->plotTitles).c_str(),
nBins*6, this->mcHist->GetBinLowEdge(1), this->mcHist->GetBinLowEdge(nBins+1) );
mcStat = (TH1D*) mcHist->Clone();
mcStat->Reset();
this->mcHist->Reset();
this->mcFine->Reset();
// Setup the NEUT Mode Array
PlotUtils::CreateNeutModeArray((TH1D*)this->mcHist,(TH1**)this->mcHist_PDG);
PlotUtils::ResetNeutModeArray((TH1**)this->mcHist_PDG);
return;
}
//********************************************************************
void Measurement1D::SetFitOptions(std::string opt){
//********************************************************************
// Do nothing if default given
if (opt == "DEFAULT") return;
-
+
// CHECK Conflicting Fit Options
std::vector<std::string> fit_option_allow = PlotUtils::FillVectorSFromString(allowed_types, "/");
for (UInt_t i = 0; i < fit_option_allow.size(); i++){
std::vector<std::string> fit_option_section = PlotUtils::FillVectorSFromString(fit_option_allow.at(i), ",");
bool found_option = false;
for (UInt_t j = 0; j < fit_option_section.size(); j++){
std::string av_opt = fit_option_section.at(j);
if (!found_option and opt.find(av_opt) != std::string::npos) {
found_option = true;
} else if (found_option and opt.find(av_opt) != std::string::npos){
ERR(FTL) << "ERROR: Conflicting fit options provided: "<<opt<<std::endl;
ERR(FTL) << "Conflicting group = "<<fit_option_section.at(i)<<std::endl;
ERR(FTL) << "You should only supply one of these options in card file."<<std::endl;
exit(-1);
}
}
}
// Check all options are allowed
std::vector<std::string> fit_options_input = PlotUtils::FillVectorSFromString(opt,"/");
for (UInt_t i = 0; i < fit_options_input.size(); i++){
if (allowed_types.find(fit_options_input.at(i)) == std::string::npos){
ERR(FTL) <<"ERROR: Fit Option '"<<fit_options_input.at(i)<<"' Provided is not allowed for this measurement."<<std::endl;
ERR(FTL) <<"Fit Options should be provided as a '/' seperated list (e.g. FREE/DIAG/NORM)" << std::endl;
ERR(FTL) <<"Available options for "<<measurementName<<" are '"<< allowed_types <<"'"<<std::endl;
exit(-1);
}
}
// Set TYPE
this->fitType = opt;
// FIX,SHAPE,FREE
if (opt.find("FIX") != std::string::npos){
isFree = isShape = false;
isFix = true;
} else if (opt.find("SHAPE") != std::string::npos){
isFree = isFix = false;
isShape = true;
} else if (opt.find("FREE") != std::string::npos){
isFix = isShape = false;
isFree = true;
}
// DIAG,FULL (or default to full)
if (opt.find("DIAG") != std::string::npos){
isDiag = true;
isFull = false;
} else if (opt.find("FULL") != std::string::npos){
isDiag = false;
isFull = true;
}
// CHI2/LL (OTHERS?)
if (opt.find("LOG") != std::string::npos) isChi2 = false;
else isChi2 = true;
// EXTRAS
if (opt.find("RAW") != std::string::npos) isRawEvents = true;
if (opt.find("DIF") != std::string::npos) isDifXSec = true;
if (opt.find("ENU1D") != std::string::npos) isEnu1D = true;
if (opt.find("NORM") != std::string::npos) addNormPenalty = true;
if (opt.find("MASK") != std::string::npos) isMask = true;
return;
};
//********************************************************************
void Measurement1D::SetDataValues(std::string dataFile) {
//********************************************************************
// Override this function if the input file isn't in a suitable format
LOG(SAM) << "Reading data from: " << dataFile.c_str() << std::endl;
this->dataHist = PlotUtils::GetTH1DFromFile(dataFile, (this->measurementName+"_data"), this->plotTitles);
this->dataTrue = (TH1D*)this->dataHist->Clone();
return;
};
//********************************************************************
void Measurement1D::SetDataFromDatabase(std::string inhistfile, std::string histname){
//********************************************************************
LOG(SAM) << "Filling histogram from "<< inhistfile << "->"<< histname <<std::endl;
this->dataHist = PlotUtils::GetTH1DFromRootFile((std::string(std::getenv("EXT_FIT")) + "/data/" + inhistfile), histname);
this->dataHist->SetNameTitle((this->measurementName+"_data").c_str(), (this->measurementName + "_data").c_str());
return;
};
//********************************************************************
void Measurement1D::SetDataFromFile(std::string inhistfile, std::string histname){
//********************************************************************
LOG(SAM) << "Filling histogram from "<< inhistfile << "->"<< histname <<std::endl;
this->dataHist = PlotUtils::GetTH1DFromRootFile((inhistfile), histname);
return;
};
//********************************************************************
void Measurement1D::SetCovarMatrix(std::string covarFile){
//********************************************************************
// Covariance function, only really used when reading in the MB Covariances.
TFile* tempFile = new TFile(covarFile.c_str(),"READ");
TH2D* covarPlot = new TH2D();
// TH2D* decmpPlot = new TH2D();
TH2D* covarInvPlot = new TH2D();
TH2D* fullcovarPlot = new TH2D();
std::string covName = "";
std::string covOption = FitPar::Config().GetParS("thrown_covariance");
if (this->isShape || this->isFree) covName = "shp_";
if (this->isDiag) covName += "diag";
else covName += "full";
covarPlot = (TH2D*) tempFile->Get((covName + "cov").c_str());
covarInvPlot = (TH2D*) tempFile->Get((covName + "covinv").c_str());
if (!covOption.compare("SUB")) fullcovarPlot = (TH2D*) tempFile->Get((covName + "cov").c_str());
else if (!covOption.compare("FULL")) fullcovarPlot = (TH2D*) tempFile->Get("fullcov");
else ERR(WRN) <<"Incorrect thrown_covariance option in parameters."<<std::endl;
int dim = int(this->dataHist->GetNbinsX());//-this->masked->Integral());
int covdim = int(this->dataHist->GetNbinsX());
this->covar = new TMatrixDSym(dim);
this->fullcovar = new TMatrixDSym(dim);
this->decomp = new TMatrixDSym(dim);
int row,column = 0;
row = 0;
column = 0;
for (Int_t i = 0; i < covdim; i++){
// if (this->masked->GetBinContent(i+1) > 0) continue;
for (Int_t j = 0; j < covdim; j++){
// if (this->masked->GetBinContent(j+1) > 0) continue;
(*this->covar)(row, column) = covarPlot->GetBinContent(i+1,j+1);
(*this->fullcovar)(row, column) = fullcovarPlot->GetBinContent(i+1,j+1);
column++;
}
column = 0;
row++;
}
// Set bin errors on data
if (!this->isDiag){
StatUtils::SetDataErrorFromCov(dataHist, fullcovar);
}
// Get Deteriminant and inverse matrix
covDet = this->covar->Determinant();
TDecompSVD LU = TDecompSVD(*this->covar);
this->covar = new TMatrixDSym(dim, LU .Invert().GetMatrixArray(), "");
return;
};
//********************************************************************
void Measurement1D::SetCovarMatrixFromText(std::string covarFile, int dim){
//********************************************************************
// WARNING this reads in the data CORRELATIONS
// Make a counter to track the line number
int row = 0;
std::string line;
std::ifstream covar(covarFile.c_str(),ifstream::in);
this->covar = new TMatrixDSym(dim);
this->fullcovar = new TMatrixDSym(dim);
if(covar.is_open()) LOG(SAM) << "Reading covariance matrix from file: " << covarFile << std::endl;
else ERR(FTL) <<"Covariance matrix provided is incorrect: "<<covarFile<<std::endl;
// MINERvA CC1pip needs slightly different method
// Only half the covariance matrix is given and I'm too lazy to write the full one so let the code do it instead
if (measurementName.find("MINERvA_CC1pip") == std::string::npos && measurementName.find("MINERvA_CCNpip") == std::string::npos) {
while(std::getline(covar, line, '\n')){
std::istringstream stream(line);
double entry;
int column = 0;
// Loop over entries and insert them into matrix
// Multiply by the errors to get the covariance, rather than the correlation matrix
while(stream >> entry){
double val = entry * this->dataHist->GetBinError(row+1)*1E38*this->dataHist->GetBinError(column+1)*1E38;
(*this->covar)(row, column) = val;
(*this->fullcovar)(row, column) = val;
column++;
}
row++;
}
// Robust matrix inversion method
TDecompSVD LU = TDecompSVD(*this->covar);
this->covar = new TMatrixDSym(dim, LU .Invert().GetMatrixArray(), "");
} else { // Here's the MINERvA CC1pip method; very similar
while(std::getline(covar, line, '\n')){
std::istringstream stream(line);
double entry;
int column = 0;
while (column < dim) {
if (column < row) {
(*this->covar)(row,column) = (*this->covar)(column,row);
column++;
} else {
while(stream >> entry){
double val = entry*(this->dataHist->GetBinError(row+1)*1E38*this->dataHist->GetBinError(column+1)*1E38); // need in these units to do Cholesky
(*this->covar)(row, column) = val;
(*this->fullcovar)(row, column) = val;
column++;
}
}
}
row++;
}
// Robust matrix inversion method
TDecompChol a = TDecompChol(*this->covar);
this->covar = new TMatrixDSym(dim, a.Invert().GetMatrixArray(), "");
} // end special treatment for MINERvA CC1pi+
return;
};
//********************************************************************
void Measurement1D::SetSmearingMatrix(std::string smearFile, int true_dim, int reco_dim){
//********************************************************************
// The smearing matrix describes the migration from true bins (rows) to reco bins (columns)
// Counter over the true bins!
int row = 0;
std::string line;
std::ifstream smear(smearFile.c_str(),ifstream::in);
// Note that the smearing matrix may be rectangular.
this->smear = new TMatrixD(true_dim, reco_dim);
if(smear.is_open()) LOG(SAM) << "Reading smearing matrix from file: " << smearFile << std::endl;
else ERR(FTL) <<"Smearing matrix provided is incorrect: "<< smearFile <<std::endl;
while(std::getline(smear, line, '\n')){
std::istringstream stream(line);
double entry;
int column = 0;
while(stream >> entry){
double val = entry;
(*this->smear)(row, column) = val/100.; // Convert to fraction from percentage (this may not be general enough)
column++;
}
row++;
}
return;
}
//********************************************************************
void Measurement1D::SetCovarFromDataFile(std::string covarFile, std::string covName){
//********************************************************************
LOG(SAM) << "Getting covariance from "<<covarFile<<"->"<<covName<<std::endl;
TFile* tempFile = new TFile(covarFile.c_str(),"READ");
TH2D* covPlot = (TH2D*) tempFile->Get(covName.c_str());
covPlot->SetDirectory(0);
int dim = covPlot->GetNbinsX();
this->fullcovar = new TMatrixDSym(dim);
for (int i = 0; i < dim; i++){
for (int j = 0; j < dim; j++){
(*this->fullcovar)(i,j) = covPlot->GetBinContent(i+1,j+1);
}
}
this->covar = (TMatrixDSym*) this->fullcovar->Clone();
this->decomp =(TMatrixDSym*) this->fullcovar->Clone();
TDecompSVD LU = TDecompSVD(*this->covar);
this->covar = new TMatrixDSym(dim, LU .Invert().GetMatrixArray(), "");
TDecompChol LUChol = TDecompChol(*this->decomp);
LUChol.Decompose();
this->decomp = new TMatrixDSym(dim, LU .GetU().GetMatrixArray(), "");
return;
};
//********************************************************************
void Measurement1D::SetBinMask(std::string maskFile){
//********************************************************************
// Create a mask histogram.
int nbins = this->dataHist->GetNbinsX();
this->maskHist = new TH1I((this->measurementName+"_maskHist").c_str(),(this->measurementName+"_maskHist; Bin; Mask?").c_str(),nbins,0,nbins);
std::string line;
std::ifstream mask(maskFile.c_str(),ifstream::in);
if (mask.is_open()) LOG(SAM) <<"Reading bin mask from file: "<<maskFile <<std::endl;
else std::cerr <<" Cannot find mask file."<<std::endl;
while(std::getline(mask, line, '\n')){
std::istringstream stream(line);
int column = 0;
double entry;
int bin;
while (stream >> entry){
if (column == 0) bin = int(entry);
if (column > 1) break;
column++;
}
this->maskHist->SetBinContent(bin,entry);
}
// Set masked data bins to zero
PlotUtils::MaskBins(this->dataHist, this->maskHist);
return;
}
/*
XSec Functions
*/
//********************************************************************
void Measurement1D::SetFluxHistogram(std::string fluxFile, int minE, int maxE, double fluxNorm){
//********************************************************************
// Note this expects the flux bins to be given in terms of MeV
LOG(SAM) << "Reading flux from file: " << fluxFile << std::endl;
TGraph f(fluxFile.c_str(),"%lg %lg");
this->fluxHist = new TH1D((this->measurementName+"_flux").c_str(), (this->measurementName+"; E_{#nu} (GeV)").c_str(), f.GetN()-1, minE, maxE);
Double_t *yVal = f.GetY();
for (int i = 0; i<fluxHist->GetNbinsX(); ++i)
this->fluxHist->SetBinContent(i+1, yVal[i]*fluxNorm);
};
//********************************************************************
double Measurement1D::TotalIntegratedFlux(std::string intOpt, double low, double high){
//********************************************************************
if(GetInput()->GetType() == kGiBUU){
return 1.0;
}
// Set Energy Limits
if (low == -9999.9) low = this->EnuMin;
if (high == -9999.9) high = this->EnuMax;
int minBin = this->fluxHist->GetXaxis()->FindBin(low);
int maxBin = this->fluxHist->GetXaxis()->FindBin(high);
// Get integral over custom range
double integral = this->fluxHist->Integral(minBin, maxBin+1, intOpt.c_str());
return integral;
};
/*
Reconfigure LOOP
*/
//********************************************************************
void Measurement1D::ResetAll(){
//********************************************************************
// Simple function to reset the mc Histograms incase that is all that is needed.
// Clear histograms
this->mcHist->Reset();
this->mcFine->Reset();
this->mcStat->Reset();
PlotUtils::ResetNeutModeArray((TH1**)this->mcHist_PDG);
return;
};
//********************************************************************
void Measurement1D::FillHistograms(){
//********************************************************************
if (Signal){
this->mcHist->Fill(X_VAR, Weight);
this->mcFine->Fill(X_VAR, Weight);
this->mcStat->Fill(X_VAR, 1.0);
PlotUtils::FillNeutModeArray(mcHist_PDG, Mode, X_VAR, Weight);
}
return;
};
//********************************************************************
void Measurement1D::ScaleEvents(){
//********************************************************************
// Simple function to scale to xsec result if this is all that is needed.
// Scale bin errors correctly
TH1D* tempFine = (TH1D*) mcFine->Clone();
// Should apply different scaling for:
// 1D Enu distributions -- need bin by bin flux unfolding (bin by bin flux integration)
// 1D count distributions -- need shape scaling to data
// anything else -- flux averages
// Scaling for raw event rates
if (isRawEvents) {
PlotUtils::ScaleNeutModeArray((TH1**)this->mcHist_PDG, (dataHist->Integral()/mcHist->Integral()), "width");
this->mcHist->Scale(dataHist->Integral()/mcHist->Integral());
this->mcFine->Scale(dataHist->Integral()/mcFine->Integral());
// Scaling for XSec as function of Enu
} else if (isEnu1D) {
PlotUtils::FluxUnfoldedScaling(mcHist, fluxHist);
PlotUtils::FluxUnfoldedScaling(mcFine, fluxHist);
mcHist->Scale(scaleFactor);
mcFine->Scale(scaleFactor);
// Any other differential scaling
} else {
this->mcHist->Scale(this->scaleFactor, "width");
this->mcFine->Scale(this->scaleFactor, "width");
PlotUtils::ScaleNeutModeArray((TH1**)this->mcHist_PDG, this->scaleFactor, "width");
}
// Proper error scaling - ROOT Freaks out with xsec weights sometimes
// Scale the MC histogram
for(int i=0; i<mcStat->GetNbinsX();i++) {
if (mcStat->GetBinContent(i+1) != 0) {
this->mcHist->SetBinError(i+1, this->mcHist->GetBinContent(i+1) * mcStat->GetBinError(i+1) / mcStat->GetBinContent(i+1) );
} else {
this->mcHist->SetBinError(i+1, this->mcHist->Integral());
}
}
// Scale the fine MC histogram
for(int i=0; i<tempFine->GetNbinsX();i++) {
if (tempFine->GetBinContent(i+1) != 0) {
this->mcFine->SetBinError(i+1, this->mcFine->GetBinContent(i+1) * tempFine->GetBinError(i+1) / tempFine->GetBinContent(i+1) );
} else {
this->mcFine->SetBinError(i+1, this->mcFine->Integral());
}
}
return;
};
//********************************************************************
void Measurement1D::ApplyNormScale(double norm){
//********************************************************************
this->currentNorm = norm;
this->mcHist->Scale(1.0/norm);
this->mcFine->Scale(1.0/norm);
PlotUtils::ScaleNeutModeArray((TH1**)this->mcHist_PDG, 1.0/norm);
return;
};
//********************************************************************
void Measurement1D::ApplySmearingMatrix(){
//********************************************************************
if (!this->smear){
ERR(WRN) << this->measurementName <<": attempted to apply smearing matrix, but none was set"<<std::endl;
return;
}
TH1D* unsmeared = (TH1D*)mcHist->Clone();
TH1D* smeared = (TH1D*)mcHist->Clone();
smeared->Reset();
// Loop over reconstructed bins
// true = row; reco = column
for (int rbin=0; rbin < this->smear->GetNcols(); ++rbin){
// Sum up the constributions from all true bins
double rBinVal = 0;
// Loop over true bins
for (int tbin=0; tbin < this->smear->GetNrows(); ++tbin){
rBinVal += (*this->smear)(tbin,rbin)*unsmeared->GetBinContent(tbin+1);
}
smeared->SetBinContent(rbin+1, rBinVal);
}
this->mcHist = (TH1D*)smeared->Clone();
return;
}
/*
Statistic Functions - Outsources to StatUtils
*/
//********************************************************************
int Measurement1D::GetNDOF(){
//********************************************************************
return this->dataHist->GetNbinsX(); // - this->maskHist->Integral();
}
//********************************************************************
double Measurement1D::GetLikelihood(){
//********************************************************************
double stat = 0.0;
// Fix weird masking bug
if (!isMask){
if (maskHist){
maskHist = NULL;
}
}
// Sort Initial Scaling
double scaleF = this->dataHist->Integral(1,this->dataHist->GetNbinsX(),"width")/this->mcHist->Integral(1, this->mcHist->GetNbinsX(), "width");
if (isShape){
this->mcHist->Scale(scaleF);
this->mcFine->Scale(scaleF);
}
// Get Chi2
if (isChi2){
if (!isDiag){
if (!isChi2SVD) {
stat = StatUtils::GetChi2FromCov(dataHist, mcHist, covar, maskHist);
} else {
stat = StatUtils::GetChi2FromSVD(dataHist,mcHist, fullcovar, maskHist);
}
} else {
if (isRawEvents) {
stat = StatUtils::GetChi2FromEventRate(dataHist, mcHist, maskHist);
} else {
stat = StatUtils::GetChi2FromDiag(dataHist, mcHist, maskHist);
}
}
} else {
if (!this->isDiag){
if (!isChi2SVD) stat = StatUtils::GetLikelihoodFromCov(dataHist, mcHist, covar, maskHist);
else stat = StatUtils::GetLikelihoodFromSVD(dataHist,mcHist, fullcovar, maskHist);
} else {
if (this->isRawEvents) stat = StatUtils::GetLikelihoodFromEventRate(dataHist, mcHist, maskHist);
else stat = StatUtils::GetLikelihoodFromDiag(dataHist, mcHist, maskHist);
}
}
// Sort Penalty Terms
if (this->addNormPenalty){
double penalty = (1. - this->currentNorm)*(1. - this->currentNorm)/(this->normError*this->normError);
stat += penalty;
}
LOG(REC) << this->measurementName<<": Sample Chi^2 = " << stat <<std::endl;
// Return to normal scaling
if (this->isShape){
this->mcHist->Scale(1./scaleF);
this->mcFine->Scale(1./scaleF);
}
return stat;
}
//********************************************************************
void Measurement1D::SetFakeDataValues(std::string fakeOption) {
//********************************************************************
// Reset things
if (usingfakedata){
this->ResetFakeData();
} else {
usingfakedata = true;
}
// Make a copy of the original data histogram.
if (!(this->dataOrig)) this->dataOrig = (TH1D*)this->dataHist->Clone((this->measurementName+"_data_original").c_str());
TH1D *tempData = (TH1D*)this->dataHist->Clone();
TFile *fake = new TFile();
if (fakeOption.compare("MC")==0){
LOG(SAM) << "Setting fake data from MC "<<std::endl;
this->dataHist = (TH1D*)this->mcHist->Clone((this->measurementName+"_MC").c_str());
if (this->mcHist->Integral() == 0.0) ERR(WRN) << this->measurementName <<": Invalid histogram"<<std::endl;
}
else {
fake = new TFile(fakeOption.c_str());
this->dataHist = (TH1D*)fake->Get((this->measurementName+"_MC").c_str());
}
this->dataHist ->SetNameTitle((this->measurementName+"_FAKE").c_str(), (this->measurementName+this->plotTitles).c_str());
this->dataTrue = (TH1D*)this->dataHist->Clone();
this->dataTrue ->SetNameTitle((this->measurementName+"_FAKE_TRUE").c_str(), (this->measurementName+this->plotTitles).c_str());
int nbins = this->dataHist->GetNbinsX();
double alpha_i = 0.0;
double alpha_j = 0.0;
for (int i = 0; i < nbins; i++){
for (int j = 0; j < nbins; j++){
alpha_i = this->dataHist->GetBinContent(i+1)/tempData->GetBinContent(i+1);
alpha_j = this->dataHist->GetBinContent(j+1)/tempData->GetBinContent(j+1);
(*this->covar)(i,j) = (1.0/(alpha_i*alpha_j))*(*this->covar)(i,j);
(*this->fullcovar)(i,j) = alpha_i*alpha_j*(*this->fullcovar)(i,j);
}
}
(this->covar) = (TMatrixDSym*) this->fullcovar->Clone();
TDecompSVD LU = TDecompSVD(*this->covar);
this->covar = new TMatrixDSym(nbins, LU .Invert().GetMatrixArray(), "");
return;
};
//********************************************************************
void Measurement1D::ResetFakeData(){
//********************************************************************
if (usingfakedata)
if (this->dataHist) delete dataHist;
this->dataHist = (TH1D*) this->dataTrue->Clone((this->measurementName+"_FKDAT").c_str());
return;
}
//********************************************************************
void Measurement1D::ResetData(){
//********************************************************************
if (usingfakedata)
if (this->dataHist) delete dataHist;
this->dataHist = (TH1D*) this->dataTrue->Clone((this->measurementName+"_Data").c_str());
usingfakedata = false;
}
//********************************************************************
void Measurement1D::ThrowCovariance(){
//********************************************************************
// Take a decomposition and use it to throw the current dataset.
// Requires dataTrue also be set incase used repeatedly.
delete dataHist;
this->dataHist = StatUtils::ThrowHistogram(this->dataTrue, this->fullcovar);
return;
};
/*
Access Functions
*/
//********************************************************************
std::vector<TH1*> Measurement1D::GetMCList(){
//********************************************************************
// If this isn't a NULL pointer, make the plot pretty!
if (!this->mcHist) return std::vector<TH1*> (1, this->mcHist);
std::ostringstream chi2;
chi2 << std::setprecision(5) << this->GetLikelihood();
int plotcolor = kRed;
if (FitPar::Config().GetParI("linecolour") > 0){
plotcolor = FitPar::Config().GetParI("linecolour");
}
int plotstyle = 1;
if (FitPar::Config().GetParI("linestyle") > 0){
plotstyle = FitPar::Config().GetParI("linestyle");
}
int plotfillstyle=0;
if (FitPar::Config().GetParI("fillstyle") > 0){
plotfillstyle = FitPar::Config().GetParI("fillstyle");
}
std::cout << measurementName << " chi2 = " << GetLikelihood() << std::endl;
this->mcHist->SetTitle(chi2.str().c_str());
this->mcHist->SetLineWidth(3);
this->mcHist->SetLineColor(plotcolor);
this->mcHist->SetFillColor(plotcolor);
this->mcHist->SetLineStyle(plotstyle);
this->mcHist->SetFillStyle(plotfillstyle);
return std::vector<TH1*> (1, this->mcHist);
};
//********************************************************************
std::vector<TH1*> Measurement1D::GetDataList(){
//********************************************************************
// If this isn't a NULL pointer, make the plot pretty!
if (!this->dataHist) return std::vector<TH1*> (1, this->dataHist);
this->dataHist->SetLineWidth(2);
this->dataHist->SetMarkerStyle(8);
this->dataHist->SetLineColor(kBlack);
return std::vector<TH1*> (1, this->dataHist);
};
//********************************************************************
void Measurement1D::GetBinContents(std::vector<double>& cont, std::vector<double>& err){
//********************************************************************
// Return a vector of the main bin contents
for (int i = 0; i < this->mcHist->GetNbinsX(); i++){
cont.push_back(this->mcHist->GetBinContent(i+1));
err.push_back(this->mcHist->GetBinError(i+1));
}
return;
};
//********************************************************************
std::vector<double> Measurement1D::GetXSec(std::string option){
//********************************************************************
std::vector<double> vals;
vals.push_back(0.0);
vals.push_back(0.0);
bool getMC = !option.compare("MC");
bool getDT = !option.compare("DATA");
for (int i = 0; i < this->mcHist->GetNbinsX(); i++){
if (this->dataHist->GetBinContent(i+1) == 0.0 and this->dataHist->GetBinError(i+1) == 0.0) continue;
if (getMC){
vals[0] += this->mcHist->GetBinContent(i+1) * this->mcHist->GetXaxis()->GetBinWidth(i+1);
vals[1] += this->mcHist->GetBinError(i+1) * this->mcHist->GetBinError(i+1) * this->mcHist->GetXaxis()->GetBinWidth(i+1) * this->mcHist->GetXaxis()->GetBinWidth(i+1);
} else if (getDT){
vals[0] += this->dataHist->GetBinContent(i+1) * this->dataHist->GetXaxis()->GetBinWidth(i+1);
vals[1] += this->dataHist->GetBinError(i+1) * this->dataHist->GetBinError(i+1) * this->dataHist->GetXaxis()->GetBinWidth(i+1) * this->dataHist->GetXaxis()->GetBinWidth(i+1);
}
}
// If not diag Get the total error from the covariance
if (!this->isDiag and !this->isRawEvents and getDT and fullcovar){
vals[1] = 0.0;
for (int i = 0; i < this->dataHist->GetNbinsX(); i++){
for(int j = 0; j < this->dataHist->GetNbinsX(); j++){
vals[1] += (*fullcovar)(i,j);
}
}
vals[1] = sqrt(vals[1]) * 1E-38;
}
return vals;
}
/*
Write Functions
*/
// Save all the histograms at once
//********************************************************************
void Measurement1D::Write(std::string drawOpt){
//********************************************************************
// If null pointer return
if (!this->mcHist and !this->dataHist){
LOG(SAM) << this->measurementName <<"Incomplete histogram set!"<<std::endl;
return;
}
// Get Draw Options
drawOpt = FitPar::Config().GetParS("drawopts");
bool drawData = (drawOpt.find("DATA") != std::string::npos);
bool drawNormal = (drawOpt.find("MC") != std::string::npos);
bool drawEvents = (drawOpt.find("EVT") != std::string::npos);
bool drawFine = (drawOpt.find("FINE") != std::string::npos);
bool drawRatio = (drawOpt.find("RATIO") != std::string::npos);
bool drawModes = (drawOpt.find("MODES") != std::string::npos);
bool drawShape = (drawOpt.find("SHAPE") != std::string::npos);
bool residual = (drawOpt.find("RESIDUAL") != std::string::npos);
// bool drawMatrix = (drawOpt.find("MATRIX") != std::string::npos);
// bool drawXSec = (drawOpt.find("XSEC") != std::string::npos);
bool drawFlux = (drawOpt.find("FLUX") != std::string::npos);
bool drawMask = (drawOpt.find("MASK") != std::string::npos);
bool drawCov = (drawOpt.find("COV") != std::string::npos);
bool drawInvCov = (drawOpt.find("INVCOV") != std::string::npos);
bool drawDecomp = (drawOpt.find("DECOMP") != std::string::npos);
bool drawCanvPDG = (drawOpt.find("CANVPDG") != std::string::npos);
bool drawCanvMC = (drawOpt.find("CANVMC") != std::string::npos);
-
+
LOG(SAM)<<"Writing Normal Plots" <<std::endl;
// Save standard plots
if (drawData) this->GetDataList().at(0)->Write();
if (drawNormal) this->GetMCList() .at(0)->Write();
// Save only mc and data if splines
if(this->eventType == 4 or this->eventType==3){ return; }
// Draw Extra plots
LOG(SAM)<<"Writing Fine List"<<std::endl;
if (drawFine) this->GetFineList().at(0)->Write();
LOG(SAM)<<"Writing events"<<std::endl;
if (drawFlux) this->fluxHist->Write();
LOG(SAM)<<"Writing true events"<<std::endl;
// if (drawXSec) this->xsecHist->Write();
if (drawEvents) this->eventHist->Write();
if (isMask and drawMask) this->maskHist->Write( (this->measurementName + "_MSK").c_str() ); //< save mask
// Save neut stack
if (drawModes){
LOG(SAM) << "Writing MC Hist PDG"<<std::endl;
THStack combo_mcHist_PDG = PlotUtils::GetNeutModeStack((this->measurementName + "_MC_PDG").c_str(), (TH1**)this->mcHist_PDG, 0);
combo_mcHist_PDG.Write();
}
// Save Matrix plots
if (!isRawEvents and !isDiag){
if (drawCov and fullcovar){
TH2D cov = TH2D((*this->fullcovar));
cov.SetNameTitle((this->measurementName+"_cov").c_str(),(this->measurementName+"_cov;Bins; Bins;").c_str());
cov.Write();
}
if (drawInvCov and covar){
TH2D covinv = TH2D((*this->covar));
covinv.SetNameTitle((this->measurementName+"_covinv").c_str(),(this->measurementName+"_cov;Bins; Bins;").c_str());
covinv.Write();
}
if (drawDecomp and decomp){
TH2D covdec = TH2D((*this->decomp));
covdec.SetNameTitle((this->measurementName+"_covdec").c_str(),(this->measurementName+"_cov;Bins; Bins;").c_str());
covdec.Write();
}
}
// Save ratio plots if required
if (drawRatio){
// Needed for error bars
for(int i = 0; i < this->mcHist->GetNbinsX()*this->mcHist->GetNbinsY(); i++)
this->mcHist->SetBinError(i+1,0.0);
this->dataHist->GetSumw2();
this->mcHist->GetSumw2();
// Create Ratio Histograms
TH1D* dataRatio = (TH1D*) this->dataHist->Clone((this->measurementName + "_data_RATIO").c_str());
TH1D* mcRatio = (TH1D*) this->mcHist->Clone((this->measurementName + "_MC_RATIO").c_str());
mcRatio->Divide(this->mcHist);
dataRatio->Divide(this->mcHist);
// Cancel bin errors on MC
for(int i = 0; i < mcRatio->GetNbinsX(); i++)
mcRatio->SetBinError(i+1,this->mcHist->GetBinError(i+1) / this->mcHist->GetBinContent(i+1));
mcRatio->SetMinimum(0);
mcRatio->SetMaximum(2);
dataRatio->SetMinimum(0);
dataRatio->SetMaximum(2);
mcRatio->Write();
dataRatio->Write();
delete mcRatio;
delete dataRatio;
}
// Save Shape Plots if required
if (drawShape){
// Create Shape Histogram
TH1D* mcShape = (TH1D*) this->mcHist->Clone((this->measurementName + "_MC_SHAPE").c_str());
double shapeScale = dataHist->Integral("width")/mcHist->Integral("width");
mcShape->Scale(shapeScale);
std::stringstream ss;
ss << shapeScale;
mcShape->SetTitle(ss.str().c_str());
mcShape->SetLineWidth(3);
mcShape->SetLineStyle(7); //dashes
mcShape->Write();
// Save shape ratios
if (drawRatio) {
// Needed for error bars
mcShape->GetSumw2();
// Create shape ratio histograms
TH1D* mcShapeRatio = (TH1D*)mcShape->Clone((this->measurementName + "_MC_SHAPE_RATIO").c_str());
TH1D* dataShapeRatio = (TH1D*)dataHist->Clone((this->measurementName + "_data_SHAPE_RATIO").c_str());
// Divide the histograms
mcShapeRatio ->Divide(mcShape);
dataShapeRatio ->Divide(mcShape);
// Colour the shape ratio plots
mcShapeRatio ->SetLineWidth(3);
mcShapeRatio ->SetLineStyle(7); // dashes
mcShapeRatio ->Write();
dataShapeRatio->Write();
delete mcShapeRatio;
delete dataShapeRatio;
}
delete mcShape;
}
// Save residual calculations of what contributed to the chi2 values.
if (residual){
}
// Make a pretty PDG Canvas
if (drawCanvPDG or true){
TCanvas* c1 = new TCanvas((this->measurementName + "_PDG_CANV").c_str(),
(this->measurementName + "_PDG_CANV").c_str(),
800,600);
dataHist->Draw("E1");
mcHist->Draw("HIST SAME");
-
+
THStack combo_mcHist_PDG = PlotUtils::GetNeutModeStack((this->measurementName + "_MC_PDG").c_str(),
(TH1**)this->mcHist_PDG, 0);
combo_mcHist_PDG.Draw("HIST SAME");
TLegend leg = PlotUtils::GenerateStackLegend(combo_mcHist_PDG, 0.6,0.6,0.9,0.9);
dataHist->Draw("E1 SAME");
-
+
//leg.Draw("SAME");
c1->Write();
}
if (drawCanvMC or true){
TCanvas* c1 = new TCanvas((this->measurementName + "_MC_CANV").c_str(),
(this->measurementName + "_MC_CANV").c_str(),
800,600);
c1->cd();
dataHist->Draw("E1");
mcHist->Draw("SAME HIST C");
-
+
TH1D* mcShape = (TH1D*) this->mcHist->Clone((this->measurementName + "_MC_SHAPE").c_str());
double shapeScale = dataHist->Integral("width")/mcHist->Integral("width");
mcShape->Scale(shapeScale);
mcShape->SetLineStyle(7);
mcShape->Draw("SAME HIST C");
TLegend* leg = new TLegend(0.6,0.6,0.9,0.9);
leg->AddEntry(dataHist, (this->measurementName + " Data").c_str(), "ep");
leg->AddEntry(mcHist, (this->measurementName + " MC").c_str(), "l");
leg->AddEntry(mcShape, (this->measurementName + " Shape").c_str(), "l");
}
-
+
// Returning
LOG(SAM) << this->measurementName << "Written Histograms: "<<this->measurementName<<std::endl;
return;
};
THStack Measurement1D::GetModeStack(){
THStack combo_hist = PlotUtils::GetNeutModeStack((this->measurementName + "_MC_PDG").c_str(), (TH1**)this->mcHist_PDG, 0);
return combo_hist;
}
diff --git a/src/T2K/CMakeLists.txt b/src/T2K/CMakeLists.txt
index 6b5b00a..1cb5a86 100644
--- a/src/T2K/CMakeLists.txt
+++ b/src/T2K/CMakeLists.txt
@@ -1,64 +1,66 @@
# 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 <http://www.gnu.org/licenses/>.
################################################################################
set(IMPLFILES
T2K_CC0pi_XSec_2DPcos_nu.cxx
T2K_CC1pip_CH_XSec_1DQ2_nu.cxx
T2K_CC1pip_CH_XSec_1DWrec_nu.cxx
T2K_CC1pip_CH_XSec_1Dpmu_nu.cxx
T2K_CC1pip_CH_XSec_1Dppi_nu.cxx
T2K_CC1pip_CH_XSec_1Dq3_nu.cxx
T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx
T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx
T2K_CC1pip_CH_XSec_1Dthq3pi_nu.cxx
+T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx
)
set(HEADERFILES
T2K_CC0pi_XSec_2DPcos_nu.h
T2K_CC1pip_CH_XSec_1DQ2_nu.h
T2K_CC1pip_CH_XSec_1DWrec_nu.h
T2K_CC1pip_CH_XSec_1Dpmu_nu.h
T2K_CC1pip_CH_XSec_1Dppi_nu.h
T2K_CC1pip_CH_XSec_1Dq3_nu.h
T2K_CC1pip_CH_XSec_1Dthmupi_nu.h
T2K_CC1pip_CH_XSec_1Dthpi_nu.h
T2K_CC1pip_CH_XSec_1Dthq3pi_nu.h
+T2K_CC0pinp_STV_XSec_1Ddpt_nu.h
)
set(LIBNAME expT2K)
if(CMAKE_BUILD_TYPE MATCHES DEBUG)
add_library(${LIBNAME} STATIC ${IMPLFILES})
else(CMAKE_BUILD_TYPE MATCHES RELEASE)
add_library(${LIBNAME} SHARED ${IMPLFILES})
endif()
include_directories(${RWENGINE_INCLUDE_DIRECTORIES})
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
include_directories(${CMAKE_SOURCE_DIR}/src/Utils)
set_target_properties(${LIBNAME} PROPERTIES VERSION
"${ExtFit_VERSION_MAJOR}.${ExtFit_VERSION_MINOR}.${ExtFit_VERSION_REVISION}")
set_target_properties(${LIBNAME} PROPERTIES LINK_FLAGS ${ROOT_LD_FLAGS})
install(TARGETS ${LIBNAME} DESTINATION lib)
#Can uncomment this to install the headers... but is it really neccessary?
#install(FILES ${HEADERFILES} DESTINATION include)
set(MODULETargets ${MODULETargets} ${LIBNAME} PARENT_SCOPE)
diff --git a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx
new file mode 100644
index 0000000..0d90ccb
--- /dev/null
+++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx
@@ -0,0 +1,61 @@
+// 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 <http://www.gnu.org/licenses/>.
+*******************************************************************************/
+
+#include "T2K_CC0pinp_STV_XSec_1Ddpt_nu.h"
+
+//********************************************************************
+T2K_CC0pinp_STV_XSec_1Ddpt_nu::T2K_CC0pinp_STV_XSec_1Ddpt_nu(
+ std::string inputfile, FitWeight *rw, std::string type,
+ std::string fakeDataFile)
+//********************************************************************
+{
+ measurementName = "T2K_CC0pinp_STV_XSec_1Ddpt_nu";
+ default_types = "FIX/DIAG/CHI2";
+ plotTitles =
+ "; #delta#it{p}_{T} (GeV c^{-1}); #frac{d#sigma}{d#delta#it{p}_{T}} "
+ "(cm^{2} neutron^{-1} GeV^{-1} c)";
+ EnuMin = 0;
+ EnuMax = 50;
+ isDiag = true;
+ Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile);
+
+ SetDataValues(std::string(std::getenv("EXT_FIT")) +
+ "/data/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.dat");
+
+ dataHist->Scale(1E-38);
+ dataTrue->Scale(1E-38);
+
+ SetupDefaultHist();
+
+ scaleFactor = eventHist->Integral("width") * double(1E-38) *
+ (13.0 / 6.0 /*Data is /neutron */) /
+ (double(nevents) * TotalIntegratedFlux("width"));
+};
+
+void T2K_CC0pinp_STV_XSec_1Ddpt_nu::FillEventVariables(FitEvent *event) {
+ X_VAR = FitUtils::Get_STV_dpt(event, true) / 1000.0;
+ return;
+};
+
+//********************************************************************
+bool T2K_CC0pinp_STV_XSec_1Ddpt_nu::isSignal(FitEvent *event)
+//********************************************************************
+{
+ return SignalDef::isT2K_CC0pi_STV(event, EnuMin, EnuMax);
+}
diff --git a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.h b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.h
new file mode 100644
index 0000000..5acc0ae
--- /dev/null
+++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.h
@@ -0,0 +1,36 @@
+// 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 <http://www.gnu.org/licenses/>.
+*******************************************************************************/
+
+#ifndef T2K_CC0pinp_STV_XSec_1Ddpt_nu_H_SEEN
+#define T2K_CC0pinp_STV_XSec_1Ddpt_nu_H_SEEN
+
+#include "Measurement1D.h"
+
+class T2K_CC0pinp_STV_XSec_1Ddpt_nu : public Measurement1D {
+public:
+ T2K_CC0pinp_STV_XSec_1Ddpt_nu(std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile);
+ virtual ~T2K_CC0pinp_STV_XSec_1Ddpt_nu() {};
+
+ void FillEventVariables(FitEvent *event);
+ bool isSignal(FitEvent *event);
+
+ private:
+};
+
+#endif
diff --git a/src/Utils/FitUtils.cxx b/src/Utils/FitUtils.cxx
index e7db207..8214e57 100644
--- a/src/Utils/FitUtils.cxx
+++ b/src/Utils/FitUtils.cxx
@@ -1,634 +1,774 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NuFiX.
*
* NuFiX 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.
*
* NuFiX 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 NuFiX. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "FitUtils.h"
/*
MISC Functions
*/
//********************************************************
-TH2D* FitUtils::CalculateQ3Cut(std::string inFile, TH2D* data, double qCut,
+TH2D *FitUtils::CalculateQ3Cut(std::string inFile, TH2D *data, double qCut,
int nuPDG, double eMin, double eMax) {
//********************************************************
// This function calculates the bins in the 2D distribution which should be
// excluded because > 50% of the events in that bin have q^2 < qCut
- TH2D* qCutHist = (TH2D*)data->Clone("qCutHist");
+ TH2D *qCutHist = (TH2D *)data->Clone("qCutHist");
qCutHist->Reset();
// DEPRECATED CODE
// - Needs to be updated to fit in with new InputHandler format
(void)inFile;
(void)qCut;
(void)nuPDG;
(void)eMin;
(void)eMax;
// This should instead use VEXperimentBase to allow multiple inputs.
// TH2D *qEvts = (TH2D*)data->Clone("qEvts");
// qEvts ->Reset();
// TH2D *qAll = (TH2D*)data->Clone("qAll");
// qAll ->Reset();
// int muPDG = 0;
// if (nuPDG == 14) muPDG = 13;
// else if (nuPDG == -14) muPDG = -13;
// // Now iterate through the events and plot the qCut
// TChain *tn = new TChain("neuttree", "");
// tn->Add(Form("%s/neuttree", inFile.c_str()));
// Int_t nevents = tn->GetEntries();
// NeutVect *nvect = NULL;
// tn->SetBranchAddress("vectorbranch", &nvect);
// for (int i = 0; i < nevents; ++i){
// tn->GetEntry(i);
// if (!isSignal(nvect, nuPDG, eMin, eMax, false, 0)) continue;
// // Now calculate the 3-momentum transfer
// // Loop over the particle stack
// for (int j = 2; j < nvect->Npart(); ++j){
// // Look for the outgoing muon
// if ((nvect->PartInfo(j))->fPID != muPDG) continue;
// double q3 = sqrt(((nvect->PartInfo(j))->fP -
// (nvect->PartInfo(0))->fP).Vect().Mag2());
// // Now find the kinematic values and fill the histogram
// double Ekmu = (nvect->PartInfo(j))->fP.E()/1000 - 0.105658367;
// double costheta =
// cos(((nvect->PartInfo(0))->fP.Vect().Angle((nvect->PartInfo(j))->fP.Vect())));
// qAll->Fill(Ekmu, costheta);
// if (q3 < qCut) qEvts ->Fill(Ekmu, costheta);
// }
// }
// for (int xBin = 1; xBin < qCutHist->GetNbinsX()+1; ++xBin){
// for (int yBin = 1; yBin < qCutHist->GetNbinsY()+1; ++yBin){
// double fract = 0;
// if (qAll->GetBinContent(xBin, yBin) != 0) fract =
// qEvts->GetBinContent(xBin, yBin)/(qAll->GetBinContent(xBin, yBin) +
// 0.);
// qCutHist->SetBinContent(xBin, yBin, fract);
// }
// }
return qCutHist;
};
/*
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;
};
//********************************************************************
// Returns the angle between two particles in radians
double FitUtils::th(TLorentzVector part1, TLorentzVector part2) {
//********************************************************************
double th = part1.Vect().Angle(part2.Vect());
return th;
};
// 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);
const double m_n = 0.93956536; // neutron/proton mass
double a1 = (E_nu + m_n) - 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 = 0.93956536; // neutron mass
const double mp = 0.93827203; // 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) {
//********************************************************************
double momshift = 0.0;
double temp = FitPar::Config().GetParD("muon_momentum_shift");
if (temp != -999.9 and temp != 0.0) {
if (FitPar::Config().GetParI("muon_momentum_throw") == 0)
momshift = temp;
else if (FitPar::Config().GetParI("muon_momentum_throw") == 1) {
momshift = gRandom->Gaus(0.0, 1.0) * temp;
}
}
// std::cout<<"Current Momentum Shift = "<<momshift<<std::endl;
// Convert all values to GeV
const double V = binding / 1000.; // binding potential
const double mn = 0.93956536; // neutron mass
const double mp = 0.93827203; // 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
pl += momshift;
double rEnu =
(2 * mN_eff * el - ml * ml + mN_oth * mN_oth - mN_eff * mN_eff) /
(2 * (mN_eff - el + pl * costh));
// std::cout<<"Enu = "<<rEnu<<std::endl;
return rEnu;
};
double FitUtils::Q2QErec(TLorentzVector pmu, double costh, double binding,
bool neutrino) {
double momshift = 0.0;
double temp = FitPar::Config().GetParD("muon_momentum_shift");
if (temp != -999.9 and temp != 0.0) {
if (FitPar::Config().GetParI("muon_momentum_throw") == 0)
momshift = temp;
else if (FitPar::Config().GetParI("muon_momentum_throw") == 1) {
momshift = gRandom->Gaus(0.0, 1.0) * temp;
}
}
// std::cout<<"Current Q2QE Momentum Shift = "<<momshift<<std::endl;
double el = pmu.E() / 1000.;
double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton
double ml = sqrt(el * el - pl * pl); // lepton mass
pl += momshift / 1000.;
double rEnu = EnuQErec(pmu, 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 = 0.93956536; // neutron mass
const double m_p = 0.93827203; // 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 = 0.93956536; // 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 = 1.232; // PDG value for Delta mass in GeV
const double m_n = 0.93956536; // 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;
};
//********************************************************************
// 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.;
// double p_mu = pmu.Vect().Mag()/1000.;
TVector3 p_mu_vect = pmu.Vect() * (1. / 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.;
TVector3 p_pi_vect = ppi.Vect() * (1. / 1000.);
// double m_pi = sqrt(E_pi*E_pi - p_pi*p_pi);
double E_bind =
27. / 1000.; // This should be roughly correct for CH; but not clear!
double m_p = 0.939;
// 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 all outgoing particles
// Requires pion vector for reconstructing the neutrino energy
// 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
// Return value is in MeV!!!
double FitUtils::Wrec(TLorentzVector pnu, TLorentzVector pmu) {
//********************************************************
// Reconstruct the hadronic mass using all outgoing particles
// Requires pion vector for reconstructing the neutrino energy
// Could technically do E_nu = pnu.E() too, but this won't be reconstructed
// Enu; it's true Enu
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());
const double m_p = 938.27203;
// MINERvA cut on W_exp which is tuned to W_true; so use true Enu from
// generators
double E_nu = pnu.E();
// double E_nu = FitUtils::EnuCC1piprec(pnu, pmu, ppi)*1000.;
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;
};
/*
E Recoil
*/
-double FitUtils::GetErecoil_TRUE(FitEvent* event) {
+double FitUtils::GetErecoil_TRUE(FitEvent *event) {
// Get total energy of hadronic system.
double Erecoil;
- for (int i = 0; i < event->Npart(); i++) {
+ for (unsigned int i = 0; i < event->Npart(); i++) {
// Only final state
if (!event->PartInfo(i)->fIsAlive) continue;
// Skip Lepton
if (abs(event->PartInfo(i)->fPID) == abs(event->PartInfo(0)->fPID) + 1)
continue;
// Add up Erecoil
Erecoil += event->PartInfo(i)->fP.E();
}
return Erecoil;
}
-double FitUtils::GetErecoil_CHARGED(FitEvent* event) {
+double FitUtils::GetErecoil_CHARGED(FitEvent *event) {
// Get total energy of hadronic system.
double Erecoil;
- for (int i = 0; i < event->Npart(); i++) {
+ for (unsigned int i = 0; i < event->Npart(); i++) {
// Only final state
if (!event->PartInfo(i)->fIsAlive) 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 Erecoil
Erecoil += event->PartInfo(i)->fP.E();
}
return Erecoil;
}
-double FitUtils::GetErecoil_MINERvA_LowRecoil(FitEvent* event) {
+double FitUtils::GetErecoil_MINERvA_LowRecoil(FitEvent *event) {
// Get total energy of hadronic system.
double Erecoil;
- for (int i = 0; i < event->Npart(); i++) {
+ for (unsigned int i = 0; i < event->Npart(); i++) {
// Only final state
if (!event->PartInfo(i)->fIsAlive) continue;
// Skip Lepton
if (abs(event->PartInfo(i)->fPID) == 14) 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);
// Total Energy of non-neutrons
} else if (PID != 2112 and PID < 999 and PID != 22 and abs(PID) != 14) {
Erecoil += (event->PartInfo(i)->fP.E()) / 1000.0;
}
}
return Erecoil;
}
-TLorentzVector FitUtils::GetHMPDG_4Mom(int pdg, FitEvent *event) {
+std::pair<TLorentzVector, int> FitUtils::GetHMPDG_4Mom(int pdg,
+ FitEvent *event) {
int pdgs[] = {pdg};
- return GetHMPDG_4Mom(pdgs, event).first;
+ return GetHMPDG_4Mom(pdgs, event);
+}
+
+std::pair<TLorentzVector, int> FitUtils::GetHMFSCLepton_4Mom(FitEvent *event) {
+ int pdgs[] = {11, 13, 15, -11, -13, -15};
+ return GetHMPDG_4Mom(pdgs, event);
+}
+std::pair<TLorentzVector, int> FitUtils::GetHMFSMuon_4Mom(FitEvent *event) {
+ int pdgs[] = {13};
+ return GetHMPDG_4Mom(pdgs, event);
+}
+std::pair<TLorentzVector, int> FitUtils::GetHMFSAnitMuon_4Mom(FitEvent *event) {
+ int pdgs[] = {-13};
+ return GetHMPDG_4Mom(pdgs, event);
+}
+std::pair<TLorentzVector, int> FitUtils::GetHMISNLepton_4Mom(FitEvent *event) {
+ int pdgs[] = {12, 14, 16, -12, -14, -16};
+ return GetHMPDG_4Mom(pdgs, event, true);
+}
+std::pair<TLorentzVector, int> FitUtils::GetHMISMuonNeutrino_4Mom(
+ FitEvent *event) {
+ int pdgs[] = {14};
+ return GetHMPDG_4Mom(pdgs, event, true);
+}
+std::pair<TLorentzVector, int> FitUtils::GetHMISAntiMuonNeutrino_4Mom(
+ FitEvent *event) {
+ int pdgs[] = {-14};
+ return GetHMPDG_4Mom(pdgs, event, true);
+}
+std::pair<TLorentzVector, int> FitUtils::GetHMFSProton_4Mom(FitEvent *event) {
+ int pdgs[] = {2212};
+ return GetHMPDG_4Mom(pdgs, event);
+}
+std::pair<TLorentzVector, int> FitUtils::GetHMFSCPion_4Mom(FitEvent *event) {
+ int pdgs[] = {211, -211};
+ return GetHMPDG_4Mom(pdgs, event);
+}
+std::pair<TLorentzVector, int> FitUtils::GetHMFSPion_4Mom(FitEvent *event) {
+ int pdgs[] = {211, -211, 111};
+ return GetHMPDG_4Mom(pdgs, event);
+}
+
+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, bool Is0pi) {
+ std::pair<TLorentzVector, int> pmu = FitUtils::GetHMFSMuon_4Mom(event);
+ std::pair<TLorentzVector, int> pp = FitUtils::GetHMFSProton_4Mom(event);
+ std::pair<TLorentzVector, int> pnu = FitUtils::GetHMISNLepton_4Mom(event);
+ if (!pnu.second) {
+ std::cerr << "[ERROR]: Couldn't find initial state neutrino." << std::endl;
+ throw;
+ }
+ if (!pnu.second || !pp.second) {
+ return 0;
+ }
+ TVector3 const &NuP = pnu.first.Vect();
+ TVector3 const &LeptonP = pmu.first.Vect();
+ TVector3 HadronP = pp.first.Vect();
+ if (!Is0pi) {
+ std::pair<TLorentzVector, int> pp = FitUtils::GetHMFSPion_4Mom(event);
+ HadronP += pp.first.Vect();
+ }
+ return GetDeltaPT(LeptonP, HadronP, NuP).Mag();
+}
+double FitUtils::Get_STV_dphit(FitEvent *event, bool Is0pi) {
+ std::pair<TLorentzVector, int> pmu = FitUtils::GetHMFSMuon_4Mom(event);
+ std::pair<TLorentzVector, int> pp = FitUtils::GetHMFSProton_4Mom(event);
+ std::pair<TLorentzVector, int> pnu = FitUtils::GetHMISNLepton_4Mom(event);
+ if (!pnu.second) {
+ std::cerr << "[ERROR]: Couldn't find initial state neutrino." << std::endl;
+ throw;
+ }
+ if (!pnu.second || !pp.second) {
+ return 0;
+ }
+ TVector3 const &NuP = pnu.first.Vect();
+ TVector3 const &LeptonP = pmu.first.Vect();
+ TVector3 HadronP = pp.first.Vect();
+ if (!Is0pi) {
+ std::pair<TLorentzVector, int> pp = FitUtils::GetHMFSPion_4Mom(event);
+ HadronP += pp.first.Vect();
+ }
+ return GetDeltaPhiT(LeptonP, HadronP, NuP);
+}
+double FitUtils::Get_STV_dalphat(FitEvent *event, bool Is0pi) {
+ std::pair<TLorentzVector, int> pmu = FitUtils::GetHMFSMuon_4Mom(event);
+ std::pair<TLorentzVector, int> pp = FitUtils::GetHMFSProton_4Mom(event);
+ std::pair<TLorentzVector, int> pnu = FitUtils::GetHMISNLepton_4Mom(event);
+ if (!pnu.second) {
+ std::cerr << "[ERROR]: Couldn't find initial state neutrino." << std::endl;
+ throw;
+ }
+ if (!pnu.second || !pp.second) {
+ return 0;
+ }
+ TVector3 const &NuP = pnu.first.Vect();
+ TVector3 const &LeptonP = pmu.first.Vect();
+ TVector3 HadronP = pp.first.Vect();
+ if (!Is0pi) {
+ std::pair<TLorentzVector, int> pp = FitUtils::GetHMFSPion_4Mom(event);
+ HadronP += pp.first.Vect();
+ }
+ return GetDeltaAlphaT(LeptonP, HadronP, NuP);
}
diff --git a/src/Utils/FitUtils.h b/src/Utils/FitUtils.h
index 0b52a49..963ce79 100644
--- a/src/Utils/FitUtils.h
+++ b/src/Utils/FitUtils.h
@@ -1,192 +1,217 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NuFiX.
*
* NuFiX 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.
*
* NuFiX 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 NuFiX. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#ifndef FITUTILS_H_SEEN
#define FITUTILS_H_SEEN
// C Includes
#include <math.h>
#include <stdlib.h>
#include <unistd.h>
#include <ctime>
#include <iostream>
#include <numeric>
// ROOT includes
#include <TChain.h>
#include <TFile.h>
#include <TH1D.h>
#include <TH2D.h>
#include <THStack.h>
#include <TKey.h>
#include <TLegend.h>
#include <TList.h>
#include <TLorentzVector.h>
#include <TObjArray.h>
#include <TROOT.h>
#include <TRandom3.h>
#include <TTree.h>
#include "FitEvent.h"
#include "TGraph.h"
#include "TH2Poly.h"
// Fit includes
#include "FitParameters.h"
/*!
* \addtogroup Utils
* @{
*/
//! Functions needed by individual samples for calculating kinematic quantities.
namespace FitUtils {
/*
MISC
*/
//! Function to produce a histogram indicating the fraction of events which have
//! q3 < qCut
-TH2D* CalculateQ3Cut(std::string inFile, TH2D* data, double qCut, int nuPDG,
+TH2D *CalculateQ3Cut(std::string inFile, TH2D *data, double qCut, int nuPDG,
double eMin, double eMax);
/*
MISC Event
*/
//! Returns kinetic energy of particle
double T(TLorentzVector part);
//! Returns momentum of particle
double p(TLorentzVector part);
//! Returns angle between particles (_NOT_ cosine!)
double th(TLorentzVector part, TLorentzVector part2);
//! Hadronic mass reconstruction
double Wrec(TLorentzVector pnu, TLorentzVector pmu);
/*
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);
/*
CCQE1p MINERvA
*/
//! Reconstruct Q2QE given just the maximum energy proton.
double ProtonQ2QErec(double pE, double binding);
/*
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);
/*
E Recoil Calculations
*/
-double GetErecoil_TRUE(FitEvent* event);
-double GetErecoil_CHARGED(FitEvent* event);
-double GetErecoil_MINERvA_LowRecoil(FitEvent* event);
+double GetErecoil_TRUE(FitEvent *event);
+double GetErecoil_CHARGED(FitEvent *event);
+double GetErecoil_MINERvA_LowRecoil(FitEvent *event);
/// Gets the 4 mom of the highest momentum alive particle with a PDG
/// contained within pdgs.
template <size_t N>
std::pair<TLorentzVector, int> GetHMPDG_4Mom(int const (&pdgs)[N],
- FitEvent *event) {
+ FitEvent *event,
+ bool IsInitState = false) {
TLorentzVector CHM(0, 0, 0, 0);
int HMPDG = 0;
for (int i = 0; i < event->Npart(); i++) {
// Only final state
- if (!event->PartInfo(i)->fIsAlive) continue;
+ if ((IsInitState && event->PartInfo(i)->fIsAlive) ||
+ (!IsInitState && !event->PartInfo(i)->fIsAlive)) {
+ continue;
+ }
bool found = false;
for (size_t it = 0; it < N; ++it) {
if (event->PartInfo(i)->fPID != pdgs[it]) continue;
found = true;
break;
}
if (!found) continue;
if (event->PartInfo(i)->fP.Vect().Mag2() > CHM.Vect().Mag2()) {
CHM = event->PartInfo(i)->fP;
HMPDG = event->PartInfo(i)->fPID;
}
}
return std::make_pair(CHM, HMPDG);
}
/// Gets the 4 mom of the highest momentum alive particle with a PDG
/// contained equal to pdg.
-TLorentzVector GetHMPDG_4Mom(int pdg, FitEvent *event);
-
+std::pair<TLorentzVector, int> GetHMPDG_4Mom(int pdg, FitEvent *event);
+/// Gets the 4 mom of the highest momentum final state charged lepton.
+std::pair<TLorentzVector, int> GetHMFSCLepton_4Mom(FitEvent *event);
+/// Gets the 4 mom of the highest momentum final state muon.
+std::pair<TLorentzVector, int> GetHMFSMuon_4Mom(FitEvent *event);
+/// Gets the 4 mom of the highest momentum final state anti-muon.
+std::pair<TLorentzVector, int> GetHMFSAnitMuon_4Mom(FitEvent *event);
+/// Gets the 4 mom of the highest momentum initial state neutral lepton.
+std::pair<TLorentzVector, int> GetHMISNLepton_4Mom(FitEvent *event);
+/// Gets the 4 mom of the highest momentum initial state muon neutrino.
+std::pair<TLorentzVector, int> GetHMISMuonNeutrino_4Mom(FitEvent *event);
+/// Gets the 4 mom of the highest momentum initial state anti-muon neutrino.
+std::pair<TLorentzVector, int> GetHMISAntiMuonNeutrino_4Mom(FitEvent *event);
+/// Gets the 4 mom of the highest momentum final state proton.
+std::pair<TLorentzVector, int> GetHMFSProton_4Mom(FitEvent *event);
+/// Gets the 4 mom of the highest momentum final state charged pion.
+std::pair<TLorentzVector, int> GetHMFSCPion_4Mom(FitEvent *event);
+/// Gets the 4 mom of the highest momentum final state pion.
+std::pair<TLorentzVector, int> GetHMFSPion_4Mom(FitEvent *event);
+
+double Get_STV_dpt(FitEvent *event, bool Is0pi);
+double Get_STV_dphit(FitEvent *event, bool Is0pi);
+double Get_STV_dalphat(FitEvent *event, bool Is0pi);
}
/*! @} */
#endif
diff --git a/src/Utils/SignalDef.cxx b/src/Utils/SignalDef.cxx
index 15fa794..a6265ab 100644
--- a/src/Utils/SignalDef.cxx
+++ b/src/Utils/SignalDef.cxx
@@ -1,997 +1,1022 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NuFiX.
*
* NuFiX 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.
*
* NuFiX 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 NuFiX. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "FitUtils.h"
#include "SignalDef.h"
bool SignalDef::isCCQE(FitEvent *event, double EnuMin, double EnuMax,
bool isRestricted) {
if (event->Mode != 1 && event->Mode != 2) return false;
if (event->PartInfo(0)->fPID != 14) return false;
if ((event->PartInfo(0)->fP.E() < EnuMin * 1000.) ||
(event->PartInfo(0)->fP.E() > EnuMax * 1000.))
return false;
// NO MESONS and ONE MUON
int lepCnt = 0;
TLorentzVector pnu = (event->PartInfo(0))->fP;
TLorentzVector pmu;
for (unsigned int j = 2; j < event->Npart(); j++) {
if (!((event->PartInfo(j))->fIsAlive) && (event->PartInfo(j))->fStatus != 0)
continue; // move on if NOT ALIVE and NOT NORMAL
int PID = (event->PartInfo(j))->fPID;
if (abs(PID) >= 111 && abs(PID) <= 557)
return false;
else if (abs(PID) == 11 || abs(PID) == 13 || abs(PID) == 15 ||
abs(PID) == 17) {
lepCnt++;
if (isRestricted) pmu = (event->PartInfo(j))->fP;
}
}
if (lepCnt != 1) return false;
if (isRestricted) {
double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI;
if (th_nu_mu >= 20) return false;
}
return true;
};
bool SignalDef::isCCQEBar(FitEvent *event, double EnuMin, double EnuMax,
bool isRestricted) {
// checks mode is correct
if (event->Mode != -1 && event->Mode != -2) return false;
if (event->PartInfo(0)->fPID != -14) return false;
if ((event->PartInfo(0)->fP.E() < EnuMin * 1000.) ||
(event->PartInfo(0)->fP.E() > EnuMax * 1000.))
return false;
// checks that we have NO mesons and ONE muon
int lepCnt = 0;
TLorentzVector pnu = (event->PartInfo(0))->fP;
TLorentzVector pmu;
for (unsigned int j = 2; j < event->Npart(); j++) {
if (!((event->PartInfo(j))->fIsAlive) && (event->PartInfo(j))->fStatus != 0)
continue; // move on if NOT ALIVE and NOT NORMAL
int PID = (event->PartInfo(j))->fPID;
if (abs(PID) >= 111 && abs(PID) <= 557)
return false;
else if (abs(PID) == 11 || abs(PID) == 13 || abs(PID) == 15 ||
abs(PID) == 17) {
lepCnt++;
if (isRestricted) pmu = (event->PartInfo(j))->fP;
}
}
if (lepCnt != 1) return false;
if (isRestricted) {
double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI;
if (th_nu_mu >= 20) return false;
}
return true;
};
bool SignalDef::isCCQELike(FitEvent *event, double EnuMin, double EnuMax) {
if (event->PartInfo(0)->fPID != 14) return false;
if ((event->PartInfo(0)->fP.E() < EnuMin * 1000.) ||
(event->PartInfo(0)->fP.E() > EnuMax * 1000.))
return false;
if (((event->PartInfo(2))->fPID != 13) && ((event->PartInfo(3))->fPID != 13))
return false;
int lepCnt = 0;
for (unsigned int j = 2; j < event->Npart(); j++) {
if (!(event->PartInfo(j))->fIsAlive && (event->PartInfo(j))->fStatus != 0)
continue;
int PID = (event->PartInfo(j))->fPID;
if (abs(PID) >= 110 && abs(PID) <= 557) return false;
// else if (abs(PID) == 1114 || abs(PID) == 2114 || (abs(PID) >= 2214 &&
// abs(PID) <= 5554)) return false;
else if (abs(PID) == 11 || abs(PID) == 13 || abs(PID) == 15 ||
abs(PID) == 17)
lepCnt++;
}
if (lepCnt != 1) return false;
return true;
};
bool SignalDef::isCCQELikeBar(FitEvent *event, double EnuMin, double EnuMax) {
if (event->PartInfo(0)->fPID != -14) return false;
if ((event->PartInfo(0)->fP.E() < EnuMin * 1000.) ||
(event->PartInfo(0)->fP.E() > EnuMax * 1000.))
return false;
for (unsigned int j = 2; j < event->Npart(); j++) {
if (!(event->PartInfo(j))->fIsAlive && (event->PartInfo(j))->fStatus != 0)
continue; // maybe need not 2212 2112?
if ((event->PartInfo(j))->fPID != 22 && // photon OK
(event->PartInfo(j))->fPID != 2212 && // neutron OK
(event->PartInfo(j))->fPID != 2112 && // proton OK
(event->PartInfo(j))->fPID != -13) // muon OK
return false;
}
return true;
};
bool SignalDef::isMiniBooNE_CCQELike(FitEvent *event, double EnuMin,
double EnuMax) {
if (abs(event->PartInfo(0)->fPID) != 14) return false;
if ((event->PartInfo(0)->fP.E() < EnuMin * 1000.) ||
(event->PartInfo(0)->fP.E() > EnuMax * 1000.))
return false;
int lepCnt = 0;
for (unsigned int j = 2; j < event->Npart(); j++) {
if (!(event->PartInfo(j))->fIsAlive && (event->PartInfo(j))->fStatus != 0)
continue; // maybe need not 2212 2112?
int PID = event->PartInfo(j)->fPID;
if (abs(PID) == 11 || abs(PID) == 13 || abs(PID) == 15 || abs(PID) == 17)
lepCnt++;
if ((event->PartInfo(j))->fPID != 22 && // photon OK
(event->PartInfo(j))->fPID != 2212 && // neutron OK
(event->PartInfo(j))->fPID != 2112 && // proton OK
abs((event->PartInfo(j))->fPID) != 13) // muon OK
return false;
}
if (lepCnt != 1) return false;
return true;
};
bool SignalDef::isMiniBooNE_CCQE(FitEvent *event, double EnuMin,
double EnuMax) {
// Only NUMU
if (event->PartInfo(0)->fPID != 14) return false;
// E Within Range
if ((event->PartInfo(0)->fP.E() < EnuMin * 1000.) ||
(event->PartInfo(0)->fP.E() > EnuMax * 1000.))
return false;
// Mode == 1 or 2
if (event->Mode != 2 and event->Mode != 1) return false;
return true;
}
bool SignalDef::isMiniBooNE_CCQEBar(FitEvent *event, double EnuMin,
double EnuMax) {
if (event->PartInfo(0)->fPID != -14) return false;
if ((event->PartInfo(0)->fP.E() < EnuMin * 1000.) ||
(event->PartInfo(0)->fP.E() > EnuMax * 1000.))
return false;
if (event->Mode != -2 and event->Mode != -1) return false;
return true;
}
// *********************************************
// MiniBooNE CC1pi+ signal definition
// Warning: This one is a little scary because there's a W = 1.35 GeV cut for
// signal in the selection
// Although this is unfolded over and is filled up with NUANCE
// So there is actually no W cut applied, but everything above W = 1.35
// GeV is NUANCE...
//
// The signal definition is:
// Exactly one negative muon
// Exactly one positive pion
// No other mesons
// No requirements on photons, nucleons and
// multi-nucleons
// Doesn't mention other leptons
//
// Additionally, it asks for 2 Michel e- from decay of muon and pion
// So there is good purity and we can be fairly sure that the positive pion is a
// positive pion
bool SignalDef::isCC1pip_MiniBooNE(FitEvent *event, double EnuMin,
double EnuMax) {
// *********************************************
// Do some initial checks on the incoming neutrino
// Make sure it's a muon neutrino
if ((event->PartInfo(0))->fPID != 14) return false;
// Make sure the muon neutrino is within the E_nu defined at the experiment
if (((event->PartInfo(0))->fP.E() < EnuMin * 1000.) ||
((event->PartInfo(0))->fP.E() > EnuMax * 1000.))
return false;
// Make sure the outgoing lepton is a muon
if (((event->PartInfo(2))->fPID != 13) && ((event->PartInfo(3))->fPID != 13))
return false;
int pipCnt = 0; // Counts number of pions
int lepCnt = 0; // Counts number of muons
for (unsigned int j = 2; j < event->Npart(); j++) {
if (!((event->PartInfo(j))->fIsAlive) && (event->PartInfo(j))->fStatus != 0)
continue; // Move on if NOT ALIVE and NOT NORMAL
int PID = (event->PartInfo(j))->fPID;
// Reject other mesons than pi+
if ((abs(PID) >= 111 && abs(PID) <= 210) ||
(abs(PID) >= 212 && abs(PID) <= 557) || PID == -211) {
return false;
// Reject other leptons
} else if (abs(PID) == 11 || PID == -13 || abs(PID) == 15 ||
abs(PID) == 17) {
return false;
// Count the number of muons
} else if (PID == 13) {
lepCnt++;
// Count the number of pions
} else if (PID == 211) {
pipCnt++;
}
}
// Make sure there's only one pion
if (pipCnt != 1) {
return false;
}
// Make sure there's only one muon
if (lepCnt != 1) {
return false;
}
// If it's passed all of the above we're good and have passed the selection
return true;
};
// **************************************************
// MiniBooNE CC1pi0 signal definition
//
// The signal definition is:
// Exactly one negative muon
// Exactly one pi0
// No additional mesons
// Any number of nucleons
//
// Does a few clever cuts on the likelihood to reduce CCQE contamination by
// looking at "fuzziness" of the ring; CCQE events are sharp, CC1pi0 are fuzzy
// (because of the pi0->2 gamma collinearity)
bool SignalDef::isCC1pi0_MiniBooNE(FitEvent *event, double EnuMin,
double EnuMax) {
// **************************************************
if ((event->PartInfo(0))->fPID != 14) return false;
if (((event->PartInfo(0))->fP.E() < EnuMin * 1000.) ||
((event->PartInfo(0))->fP.E() > EnuMax * 1000.))
return false;
if (((event->PartInfo(2))->fPID != 13) && ((event->PartInfo(3))->fPID != 13))
return false;
int pi0Cnt = 0;
int lepCnt = 0;
for (unsigned int j = 2; j < event->Npart(); j++) {
if (!((event->PartInfo(j))->fIsAlive) && (event->PartInfo(j))->fStatus != 0)
continue; // move to next particle if NOT ALIVE and NOT NORMAL
int PID = (event->PartInfo(j))->fPID;
// Reject if any other mesons
if (abs(PID) >= 113 && abs(PID) <= 557) {
return false;
// Reject other leptons
} else if (abs(PID) == 11 || PID == -13 || abs(PID) == 15) {
return false;
// Count number of leptons
} else if (PID == 13) {
lepCnt++;
// Count number of pi0
} else if (PID == 111) {
pi0Cnt++;
}
}
if (pi0Cnt != 1) return false;
if (lepCnt != 1) return false;
return true;
};
// **************************************
// MINERvA CC1pi0 in anti-neutrino mode
// Unfortunately there's no information on the neutrino component which is
// subtracted off
//
// 2014 analysis:
// Exactly one positive muon
// Exactly one observed pi0
// No pi+/pi allowed
// No information on what is done with mesons, oops?
// No information on what is done with nucleons, oops?
//
// 2016 analysis:
// Exactly one positive muon
// Exactly one observed pi0
// No other mesons
// No restriction on number of nucleons
//
bool SignalDef::isCC1pi0Bar_MINERvA(FitEvent *event, double EnuMin,
double EnuMax) {
// **************************************
if ((event->PartInfo(0))->fPID != -14) return false;
if (((event->PartInfo(0))->fP.E() < EnuMin * 1000.) ||
((event->PartInfo(0))->fP.E() > EnuMax * 1000.))
return false;
if (((event->PartInfo(2))->fPID != -13) &&
((event->PartInfo(3))->fPID != -13))
return false;
int pi0Cnt = 0;
int lepCnt = 0;
for (unsigned int j = 2; j < event->Npart(); j++) {
if (!((event->PartInfo(j))->fIsAlive) && (event->PartInfo(j))->fStatus != 0)
continue; // Move to next particle if NOT ALIVE and NOT NORMAL
int PID = (event->PartInfo(j))->fPID;
// No other mesons
if (abs(PID) >= 113 && abs(PID) <= 557) {
return false;
} else if (PID == -13) {
lepCnt++;
} else if (PID == 111) {
pi0Cnt++;
}
}
if (pi0Cnt != 1) return false;
if (lepCnt != 1) return false;
return true;
};
bool SignalDef::isNC1pi0_MiniBooNE(FitEvent *event, double EnuMin,
double EnuMax) {
if ((event->PartInfo(0))->fPID != 14) return false;
if (((event->PartInfo(0))->fP.E() < EnuMin * 1000.) ||
((event->PartInfo(0))->fP.E() > EnuMax * 1000.))
return false;
if (((event->PartInfo(2))->fPID != 14) && ((event->PartInfo(3))->fPID != 14))
return false;
int pi0Cnt = 0;
for (unsigned int j = 2; j < event->Npart(); j++) {
if (!((event->PartInfo(j))->fIsAlive) && (event->PartInfo(j))->fStatus != 0)
continue; // move to next particle if NOT ALIVE and NOT NORMAL
int PID = (event->PartInfo(j))->fPID;
if (abs(PID) >= 113 && abs(PID) <= 557) return false;
// else if (abs(PID) == 1114 || abs(PID) == 2114 || (abs(PID) >= 2214 &&
// abs(PID) <= 5554)) return false;
else if (abs(PID) == 11 || abs(PID) == 13 || abs(PID) == 15 ||
abs(PID) == 17)
return false;
else if (PID == 111)
pi0Cnt++;
}
if (pi0Cnt != 1) return false;
return true;
};
bool SignalDef::isNC1pi0Bar_MiniBooNE(FitEvent *event, double EnuMin,
double EnuMax) {
if ((event->PartInfo(0))->fPID != -14) return false;
if (((event->PartInfo(0))->fP.E() < EnuMin * 1000.) ||
((event->PartInfo(0))->fP.E() > EnuMax * 1000.))
return false;
if (((event->PartInfo(2))->fPID != -14) &&
((event->PartInfo(3))->fPID != -14))
return false;
int pi0Cnt = 0;
for (unsigned int j = 2; j < event->Npart(); j++) {
if (!((event->PartInfo(j))->fIsAlive) && (event->PartInfo(j))->fStatus != 0)
continue; // move to next particle if NOT ALIVE and NOT NORMAL
int PID = (event->PartInfo(j))->fPID;
if (abs(PID) >= 113 && abs(PID) <= 557) return false;
// else if (abs(PID) == 1114 || abs(PID) == 2114 || (abs(PID) >= 2214 &&
// abs(PID) <= 5554)) return false;
else if (abs(PID) == 11 || abs(PID) == 13 || abs(PID) == 15 ||
abs(PID) == 17)
return false;
else if (PID == 111)
pi0Cnt++;
}
if (pi0Cnt != 1) return false;
return true;
};
bool SignalDef::isCCcoh_MINERvA(FitEvent *event, double EnuMin, double EnuMax) {
if ((event->PartInfo(0))->fPID != 14) return false;
if (((event->PartInfo(0))->fP.E() < EnuMin * 1000.) ||
((event->PartInfo(0))->fP.E() > EnuMax * 1000.))
return false;
if (((event->PartInfo(2))->fPID != 13) && ((event->PartInfo(3))->fPID != 13))
return false;
int pipCnt = 0; // counts number of pions
int lepCnt = 0;
// double vertexE = 0;
for (unsigned int j = 2; j < event->Npart(); j++) {
if (!((event->PartInfo(j))->fIsAlive) && (event->PartInfo(j))->fStatus != 0)
continue; // move on if NOT ALIVE and NOT NORMAL
int PID = (event->PartInfo(j))->fPID;
if (PID == 13)
lepCnt++;
else if (PID == 211)
pipCnt++;
else
return false; // CCcoh definition is only 1 pi, only 1 lep
}
if (pipCnt != 1) return false;
if (lepCnt != 1) return false;
return true;
};
bool SignalDef::isCCcohBar_MINERvA(FitEvent *event, double EnuMin,
double EnuMax) {
if ((event->PartInfo(0))->fPID != -14) return false;
if (((event->PartInfo(0))->fP.E() < EnuMin * 1000.) ||
((event->PartInfo(0))->fP.E() > EnuMax * 1000.))
return false;
if (((event->PartInfo(2))->fPID != -13) &&
((event->PartInfo(3))->fPID != -13))
return false;
int pipCnt = 0; // counts number of pions
int lepCnt = 0;
// double vertexE = 0;
for (unsigned int j = 2; j < event->Npart(); j++) {
if (!((event->PartInfo(j))->fIsAlive) && (event->PartInfo(j))->fStatus != 0)
continue; // move on if NOT ALIVE and NOT NORMAL
int PID = (event->PartInfo(j))->fPID;
if (PID == -13)
lepCnt++;
else if (PID == -211)
pipCnt++;
else
return false; // CCcoh definition is only 1 pi, only 1 lep
}
if (pipCnt != 1) return false;
if (lepCnt != 1) return false;
return true;
};
// *********************************
// MINERvA CC1pi+/- signal definition (2015 release)
// Note: There is a 2016 release which is different to this (listed below), but
// it is CCNpi+ and has a different W cut
// Note2: The W cut is implemented in the class implementation in MINERvA/
// rather than here so we can draw events that don't pass the W cut (W cut is
// 1.4 GeV)
// Could possibly be changed for slight speed increase since less events
// would be used
//
// MINERvA signal is slightly different to MiniBooNE
//
// Exactly one negative muon
// Exactly one charged pion (both + and -); however, there is a Michel e-
// requirement but UNCLEAR IF UNFOLDED OR NOT (so don't know if should be
// applied)
// Exactly 1 charged pion exits (so + and - charge), however, has Michel
// electron requirement, so look for + only here?
// No restriction on neutral pions or other mesons
// MINERvA has unfolded and not unfolded muon phase space for 2015
//
// Possible problems:
// 1) Should there be a pi+ only cut implemented due to Michel requirement, or
// is pi- events filled from MC?
// 2) There is a T_pi < 350 MeV cut coming from requiring a stopping pion so the
// Michel e is seen, this is also unclear if it's unfolded so any pion is OK
//
// Nice things:
// Much data given: with and without muon angle cuts and with and without shape
// only data + covariance
//
bool SignalDef::isCC1pip_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
bool isRestricted) {
// *********************************
if ((event->PartInfo(0))->fPID != 14) return false;
if (((event->PartInfo(0))->fP.E() < EnuMin * 1000.) ||
((event->PartInfo(0))->fP.E() > EnuMax * 1000.))
return false;
if (((event->PartInfo(2))->fPID != 13) && ((event->PartInfo(3))->fPID != 13))
return false;
int pipCnt = 0; // Counts number of pi+
int lepCnt = 0; // Counts number of leptons
TLorentzVector pnu = (event->PartInfo(0))->fP;
TLorentzVector pmu;
TLorentzVector ppi;
for (unsigned int j = 2; j < event->Npart(); j++) {
if (!((event->PartInfo(j))->fIsAlive) && (event->PartInfo(j))->fStatus != 0)
continue; // Move on if NOT ALIVE and NOT NORMAL
int PID = (event->PartInfo(j))->fPID;
if (PID == 13) {
lepCnt++;
// if restricted we need the muon to find it's angle and if it's visible
// in MINOS
if (isRestricted) {
pmu = (event->PartInfo(j))->fP;
}
// Signal is both pi+ and pi-
// WARNING: PI- CONTAMINATION IS FULLY GENIE BECAUSE THE MICHEL TAG
} else if (abs(PID) == 211) {
ppi = event->PartInfo(j)->fP;
pipCnt++;
}
}
// Only one pion-like track
if (pipCnt != 1) return false;
// only one lepton
if (lepCnt != 1) return false;
// Pion kinetic energy requirement for Michel tag, leave commented for now
// if (FitUtils::T(ppi)*1000. > 350 || FitUtils::T(ppi)*1000. < 35) return
// false; // Need to have a 35 to 350 MeV pion kinetic energy requirement
// MINERvA released another CC1pi+ xsec without muon unfolding!
// here the muon angle is < 20 degrees (seen in MINOS)
if (isRestricted) {
double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI;
if (th_nu_mu >= 20) return false;
}
return true;
};
// *********************************
// MINERvA CCNpi+/- signal definition from 2016 publication
// Different to CC1pi+/- listed above; additional has W < 1.8 GeV
//
// Still asks for a Michel e and still unclear if this is unfolded or not
// Says stuff like "requirement that a Michel e isolates a subsample that is
// more nearly a pi+ prodution", yet the signal definition is both pi+ and pi-?
//
// One negative muon
// At least one charged pion
// 1.5 < Enu < 10
// No restrictions on pi0 or other mesons or baryons
//
// Also writes number of pions (nPions) if studies on this want to be done...
bool SignalDef::isCCNpip_MINERvA(FitEvent *event, int &nPions, double EnuMin,
double EnuMax, bool isRestricted) {
// *********************************
if ((event->PartInfo(0))->fPID != 14) return false;
if (((event->PartInfo(0))->fP.E() < EnuMin * 1000.) ||
((event->PartInfo(0))->fP.E() > EnuMax * 1000.))
return false;
if (((event->PartInfo(2))->fPID != 13) && ((event->PartInfo(3))->fPID != 13))
return false;
int pipCnt = 0; // Counts number of pions
int lepCnt = 0; // Counts number of leptons
TLorentzVector pnu = (event->PartInfo(0))->fP;
TLorentzVector pmu;
for (unsigned int j = 2; j < event->Npart(); j++) {
if (!((event->PartInfo(j))->fIsAlive) && (event->PartInfo(j))->fStatus != 0)
continue; // Move on if NOT ALIVE and NOT NORMAL
int PID = (event->PartInfo(j))->fPID;
if (PID == 13) {
lepCnt++;
if (isRestricted) {
pmu = (event->PartInfo(j))->fP;
}
} else if (abs(PID) == 211) {
pipCnt++; // technically also allows for a pi- to be ID as pi+, but
// there's Michel electron criteria too which eliminates this
}
// Has no restrictions on mesons, neutral pions or baryons
}
// Any number of pions greater than 0
if (pipCnt == 0) return false;
// Only one lepton
if (lepCnt != 1) return false;
// Just write the number of pions so the experiment class can use it
nPions = pipCnt;
// MINERvA released another CC1pi+ xsec without muon unfolding!
// Here the muon angle is < 20 degrees (seen in MINOS)
if (isRestricted) {
double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI;
if (th_nu_mu >= 20) return false;
}
return true;
};
// T2K H2O signal definition
bool SignalDef::isCC1pip_T2K_H2O(FitEvent *event, double EnuMin,
double EnuMax) {
if ((event->PartInfo(0))->fPID != 14) return false;
if (((event->PartInfo(0))->fP.E() < EnuMin * 1000.) ||
((event->PartInfo(0))->fP.E() > EnuMax * 1000.))
return false;
if (((event->PartInfo(2))->fPID != 13) && ((event->PartInfo(3))->fPID != 13))
return false;
int pipCnt = 0; // counts number of pions
int lepCnt = 0;
TLorentzVector Pnu = (event->PartInfo(0))->fP;
TLorentzVector Ppip;
TLorentzVector Pmu;
for (unsigned int j = 2; j < event->Npart(); j++) {
if (!((event->PartInfo(j))->fIsAlive) && (event->PartInfo(j))->fStatus != 0)
continue; // move on if NOT ALIVE and NOT NORMAL
int PID = (event->PartInfo(j))->fPID;
if ((abs(PID) >= 111 && abs(PID) <= 210) ||
(abs(PID) >= 212 && abs(PID) <= 557) || PID == -211)
return false;
else if (abs(PID) == 11 || abs(PID) == 13 || abs(PID) == 15 ||
abs(PID) == 17) {
lepCnt++;
Pmu = (event->PartInfo(j))->fP;
} else if (PID == 211) {
pipCnt++;
Ppip = (event->PartInfo(j))->fP;
}
}
if (pipCnt != 1) return false;
if (lepCnt != 1) return false;
// relatively generic CC1pi+ definition done
// now measurement specific (p_mu > 200 MeV, p_pi > 200 MeV, cos th_mu > 0.3,
// cos th_pi > 0.3 in TRUE AND RECONSTRUCTED!
double p_mu = FitUtils::p(Pmu) * 1000;
double p_pi = FitUtils::p(Ppip) * 1000;
double cos_th_mu = cos(FitUtils::th(Pnu, Pmu));
double cos_th_pi = cos(FitUtils::th(Pnu, Ppip));
if (p_mu <= 200 || p_pi <= 200 || cos_th_mu <= 0.3 || cos_th_pi <= 0.3)
return false;
return true;
};
// ******************************************************
// T2K CC1pi+ CH analysis (Raquel's thesis)
// Has different phase space cuts depending on if using Michel tag or not
//
// Essentially consists of two samples: one sample which has Michel e (which we
// can't get pion direction from); this covers backwards angles quite well.
// Measurements including this sample does not have include pion kinematics cuts
// one sample which has PID in FGD and TPC
// and no Michel e. These are mostly
// forward-going so require a pion
// kinematics cut
//
// Essentially, cuts are:
// 1 negative muon
// 1 positive pion
// Any number of nucleons
// No other particles in the final state
//
bool SignalDef::isCC1pip_T2K_CH(FitEvent *event, double EnuMin, double EnuMax,
bool MichelElectron) {
// ******************************************************
if ((event->PartInfo(0))->fPID != 14) return false;
if (((event->PartInfo(0))->fP.E() < EnuMin * 1000.) ||
((event->PartInfo(0))->fP.E() > EnuMax * 1000.))
return false;
if (((event->PartInfo(2))->fPID != 13) && ((event->PartInfo(3))->fPID != 13))
return false;
int pipCnt = 0; // Counts number of positive pions
int lepCnt = 0; // Counts number of muons
TLorentzVector Pnu = (event->PartInfo(0))->fP;
TLorentzVector Ppip;
TLorentzVector Pmu;
for (unsigned int j = 2; j < event->Npart(); j++) {
if (!(event->PartInfo(j)->fIsAlive) && (event->PartInfo(j))->fStatus != 0)
continue; // Move on if NOT ALIVE and NOT NORMAL
int PID = (event->PartInfo(j))->fPID;
// No additional mesons in the final state
if ((abs(PID) >= 111 && abs(PID) <= 210) ||
(abs(PID) >= 212 && abs(PID) <= 557) || PID == -211)
return false;
// Count leptons
else if (PID == 13) {
lepCnt++;
Pmu = (event->PartInfo(j))->fP;
// Count the pi+
} else if (PID == 211) {
pipCnt++;
Ppip = (event->PartInfo(j))->fP;
}
}
// Make the cuts on the final state particles
if (pipCnt != 1) return false;
if (lepCnt != 1) return false;
// Done with relatively generic CC1pi+ definition
// If this event passes the criteria on particle counting, enforce the T2K
// ND280 phase space constraints
// Will be different if Michel tag sample is included or not
// Essentially, if there's a Michel tag we don't cut on the pion variables
double p_mu = FitUtils::p(Pmu) * 1000;
double p_pi = FitUtils::p(Ppip) * 1000;
double cos_th_mu = cos(FitUtils::th(Pnu, Pmu));
double cos_th_pi = cos(FitUtils::th(Pnu, Ppip));
// If we're using Michel e- requirement we only care about the muon restricted
// phase space and use full pion phase space
if (MichelElectron) {
// Make the cuts on the muon variables
if (p_mu <= 200 || cos_th_mu <= 0.2) {
return false;
} else {
return true;
}
// If we aren't using Michel e- (i.e. we use directional information from
// pion) we need to impose the phase space cuts on the muon AND the pion)
} else {
// Make the cuts on muon and pion variables
if (p_mu <= 200 || p_pi <= 200 || cos_th_mu <= 0.2 || cos_th_pi <= 0.2) {
return false;
} else {
return true;
}
}
// Default to false; should never fire
return false;
};
//********************************************************************
bool SignalDef::isCCQEnumu_MINERvA(FitEvent *event, double EnuMin,
double EnuMax, bool fullphasespace) {
//********************************************************************
// For now, define as the true mode being CCQE or npnh
if (event->Mode != 1 and event->Mode != 2) return false;
// Only look at numu events
if ((event->PartInfo(0))->fPID != 14) return false;
// Get Theta Variables
double ThetaMu = 999.9;
double Enu_rec = -1.0;
for (UInt_t i = 2; i < event->Npart(); i++) {
if (event->PartInfo(i)->fPID == 13) {
ThetaMu =
(event->PartInfo(0)->fP.Vect().Angle(event->PartInfo(i)->fP.Vect()));
Enu_rec =
FitUtils::EnuQErec((event->PartInfo(i))->fP, cos(ThetaMu), 34., true);
break;
}
}
// If Restricted phase space
if (!fullphasespace && ThetaMu > 0.34906585) return false;
// restrict energy range
if (event->Enu() / 1000.0 < EnuMin || event->Enu() / 1000.0 > EnuMax)
return false;
if (Enu_rec < EnuMin || Enu_rec > EnuMax) return false;
return true;
};
//********************************************************************
bool SignalDef::isCCQEnumubar_MINERvA(FitEvent *event, double EnuMin,
double EnuMax, bool fullphasespace) {
//********************************************************************
// For now, define as the true mode being CCQE or npnh
if (event->Mode != -1 and event->Mode != -2) return false;
// Only look at numu events
if ((event->PartInfo(0))->fPID != -14) return false;
// Get Theta Variables
double ThetaMu = 999.9;
double Enu_rec = -1.0;
for (UInt_t i = 2; i < event->Npart(); i++) {
if (event->PartInfo(i)->fPID == -13) {
ThetaMu =
(event->PartInfo(0)->fP.Vect().Angle(event->PartInfo(i)->fP.Vect()));
Enu_rec = FitUtils::EnuQErec((event->PartInfo(i))->fP, cos(ThetaMu), 30.,
false);
break;
}
}
// If Restricted phase space
if (!fullphasespace && ThetaMu > 0.34906585) return false;
// restrict energy range
if (event->Enu() / 1000.0 < EnuMin || event->Enu() / 1000.0 > EnuMax)
return false;
if (Enu_rec < EnuMin || Enu_rec > EnuMax) return false;
return true;
}
//********************************************************************
bool SignalDef::isCCincLowRecoil_MINERvA(FitEvent *event, double EnuMin,
double EnuMax, bool hadroncut) {
//********************************************************************
// Only look at numu events
if ((event->PartInfo(0))->fPID != 14) return false;
// Restrict true energy range
if (event->Enu() / 1000.0 < EnuMin || event->Enu() / 1000.0 > EnuMax)
return false;
// Loop Particles
int nhadrons = 0;
int nmuons = 0;
double ThetaMu = 0.0;
double Emu = 0.0;
for (UInt_t i = 2; i < event->Npart(); i++) {
if (!(event->PartInfo(i))->fIsAlive) continue;
if (event->PartInfo(i)->fStatus != 0) continue;
int PID = event->PartInfo(i)->fPID;
if (PID == 13) {
nmuons++;
ThetaMu =
event->PartInfo(i)->fP.Vect().Angle(event->PartInfo(0)->fP.Vect());
Emu = event->PartInfo(i)->fP.E() / 1000.0;
} else if (PID != 2112 and PID < 999 and PID != 22 and abs(PID) != 14) {
nhadrons++;
}
}
// Need at least one muon
if (nmuons < 1) return false;
// Require Eav > 0.0
if (hadroncut and nhadrons < 1) return false;
// Cut on muon angle greated than 20deg
if (cos(ThetaMu) < 0.93969262078) return false;
// Cut on muon energy < 1.5 GeV
if (Emu < 1.5) return false;
return true;
}
bool SignalDef::isT2K_CC0pi(FitEvent *event, double EnuMin, double EnuMax,
bool forwardgoing) {
// Only Numu
if (!event->PartInfo(0)->fPID == 14) return false;
// Cut on Energy
if (event->Enu() / 1000.0 < EnuMin || event->Enu() / 1000.0 > EnuMax)
return false;
// Particle Checks
bool only_allowed_particles = true;
bool muon_found = false;
double CosThetaMu = -9999.9;
// Loop over all particles
for (UInt_t j = 2; j < event->Npart(); ++j) {
// Get only final state
if (!(event->PartInfo(j))->fIsAlive or (event->PartInfo(j))->fStatus != 0)
continue;
// Get PDG
int particle_pdg = (event->PartInfo(j))->fPID;
// Muon section
if (particle_pdg == 13) {
CosThetaMu = cos(((event->PartInfo(0))
->fP.Vect()
.Angle((event->PartInfo(j))->fP.Vect())));
muon_found = true;
} else if (particle_pdg != 22 && particle_pdg != 2212 &&
particle_pdg != 2112 && particle_pdg != 13) {
only_allowed_particles = false;
}
}
// CC0PI Cut
if (!only_allowed_particles or !muon_found) return false;
// restricted phase space
if (forwardgoing and CosThetaMu < 0.0 and CosThetaMu != -9999.9) return false;
}
+bool SignalDef::isT2K_CC0pi_STV(FitEvent *event, double EnuMin, double EnuMax) {
+ std::pair<TLorentzVector, int> ppi = FitUtils::GetHMFSPion_4Mom(event);
+ if (ppi.second) { // 0 pi
+ return false;
+ }
+
+ std::pair<TLorentzVector, int> pmu = FitUtils::GetHMFSMuon_4Mom(event);
+ std::pair<TLorentzVector, int> pp = FitUtils::GetHMFSProton_4Mom(event);
+
+ if (!pmu.second || !pp.second) { // need mu and p
+ return false;
+ }
+
+ // mu phase space
+ if ((pmu.first.Vect().Mag() < 250) || (pmu.first.Vect().CosTheta() < -0.6)) {
+ return false;
+ }
+ // p phase space
+ if ((pp.first.Vect().Mag() < 250) || (pp.first.Vect().Mag() > 1E3) ||
+ (pp.first.Vect().CosTheta() < 0.4)) {
+ return false;
+ }
+ return true;
+}
+
bool SignalDef::isCCInc_ArgoNeuT(FitEvent *event, bool IsAnti) {
- TLorentzVector pmu = FitUtils::GetHMPDG_4Mom(IsAnti ? -13 : 13, event);
- //std::cout << pmu.Vect().Theta()*18 << std::endl;
+ TLorentzVector pmu = FitUtils::GetHMPDG_4Mom(IsAnti ? -13 : 13, event).first;
+ // std::cout << pmu.Vect().Theta()*18 << std::endl;
return (pmu.Vect().Mag2() > 0) && (pmu.E() < 25E3) &&
- ((pmu.Vect().Theta()*180./TMath::Pi()) < 36);
+ ((pmu.Vect().Theta() * 180. / TMath::Pi()) < 36);
}
diff --git a/src/Utils/SignalDef.h b/src/Utils/SignalDef.h
index e4ff3da..3569a32 100644
--- a/src/Utils/SignalDef.h
+++ b/src/Utils/SignalDef.h
@@ -1,80 +1,81 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NuFiX.
*
* NuFiX 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.
*
* NuFiX 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 NuFiX. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
// contains signal definitions for various interactions
#ifndef SIGNALDEF_H_SEEN
#define SIGNALDEF_H_SEEN
// C/C++ includes
#include <math.h>
// ROOT includes
#include <TLorentzVector.h>
// ExtFit includes
#include "FitEvent.h"
// make it a namespace
namespace SignalDef {
// These are all the interaction modes for neutrinos
bool isCCQE(FitEvent *event, double EnuMin, double EnuMax, bool isRestricted = false);
bool isCCQEBar(FitEvent *event, double EnuMin, double EnuMax, bool isRestricted = false);
bool isCCQELike(FitEvent *event, double EnuMin, double EnuMax);
bool isCCQEBar_res(FitEvent *event, double EnuMin, double EnuMax);
bool isCCQELikeBar(FitEvent *event, double EnuMin, double EnuMax);
// MiniBooNE CC1pi+ differs from MINERvA CC1pi+ differs from T2K CC1pi+!
bool isCC1pip_MiniBooNE(FitEvent *event, double EnuMin, double EnuMax);
// MINERvA has unfolded and not unfolded muon phase space
bool isCC1pip_MINERvA (FitEvent *event, double EnuMin, double EnuMax, bool isRestricted = false);
bool isCCNpip_MINERvA (FitEvent *event, int &nPions, double EnuMin, double EnuMax, bool isRestricted = false);
// T2K not unfolded phase space restrictions
bool isCC1pip_T2K_H2O(FitEvent *event, double EnuMin, double EnuMax);
bool isCC1pip_T2K_CH(FitEvent *event, double EnuMin, double EnuMax, bool Michel);
bool isCC1pi0_MiniBooNE (FitEvent *event, double EnuMin, double EnuMax);
bool isCC1pi0Bar_MINERvA (FitEvent *event, double EnuMin, double EnuMax);
bool isNC1pi0_MiniBooNE (FitEvent *event, double EnuMin, double EnuMax);
bool isNC1pi0Bar_MiniBooNE(FitEvent *event, double EnuMin, double EnuMax);
bool isCCcoh_MINERvA (FitEvent *event, double EnuMin, double EnuMax);
bool isCCcohBar_MINERvA (FitEvent *event, double EnuMin, double EnuMax);
bool isCCQEnumu_MINERvA(FitEvent* event, double EnuMin, double EnuMax, bool fullphasespace=true);
bool isCCQEnumubar_MINERvA(FitEvent* event, double EnuMin, double EnuMax, bool fullphasespace=true);
bool isCCincLowRecoil_MINERvA(FitEvent *event, double EnuMin, double EnuMax, bool hadroncut);
bool isMiniBooNE_CCQELike(FitEvent *event, double EnuMin, double EnuMax);
bool isMiniBooNE_CCQE(FitEvent *event, double EnuMin, double EnuMax);
bool isMiniBooNE_CCQEBar(FitEvent *event, double EnuMin, double EnuMax);
bool isT2K_CC0pi(FitEvent* event, double EnuMin, double EnuMax, bool forwardgoing);
+ bool isT2K_CC0pi_STV(FitEvent* event, double EnuMin, double EnuMax);
bool isCCInc_ArgoNeuT(FitEvent* event, bool IsAnti=false);
}
#endif

File Metadata

Mime Type
text/x-diff
Expires
Mon, Feb 24, 6:43 AM (1 d, 14 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4494607
Default Alt Text
(203 KB)

Event Timeline