diff --git a/app/PrepareGiBUU.cxx b/app/PrepareGiBUU.cxx index 2fd5645..3e2428b 100644 --- a/app/PrepareGiBUU.cxx +++ b/app/PrepareGiBUU.cxx @@ -1,284 +1,286 @@ #include "FitLogger.h" #include "PlotUtils.h" #include "StatUtils.h" #include "TFile.h" #include "TH1D.h" #include "TTree.h" #include #include std::string fInputFiles = ""; std::string fOutputFile = ""; std::string fFluxFile = ""; void PrintOptions(); void ParseOptions(int argc, char *argv[]); void CreateRateHistogram(std::string inputList, std::string flux, std::string output); TH1D* MakeFluxHistFromDatFile(std::string inputDatFile); int main(int argc, char *argv[]) { SETVERBOSITY(FitPar::Config().GetParI("VERBOSITY")); SETERRVERBOSITY(FitPar::Config().GetParI("ERROR")); ParseOptions(argc, argv); NUIS_LOG(FIT, "Running PrepareGiBUU"); CreateRateHistogram(fInputFiles, fFluxFile, fOutputFile); }; TH1D* MakeFluxHistFromDatFile(std::string inputDatFile){ std::vector bin_cent; std::vector bin_vals; // Loop over lines in the file std::string line; std::ifstream in_file(inputDatFile.c_str(), std::ifstream::in); while (std::getline(in_file >> std::ws, line, '\n')) { // Skip any lines that start with a # if (line.at(0) == '#') continue; std::vector entries = GeneralUtils::ParseToDbl(line, "\t"); // Append the bin centers and bin values bin_cent.push_back(entries[0]); bin_vals.push_back(entries[1]); } if (bin_cent.size() < 2 || bin_vals.size() < 2){ NUIS_ABORT(inputDatFile << " is not a well formed file!"); } // Check that the histogram has equal binning (currently required) double bin_stride = bin_cent[1] - bin_cent[0]; // Check stride uniformity for (uint s = 0; s < bin_cent.size()-1; ++s){ double this_stride = bin_cent[s+1] - bin_cent[s]; if (fabs(this_stride-bin_stride) > 1E-8){ NUIS_ABORT(inputDatFile << " does not have regular binning!"); } } // Get bin edges std::vector bin_edges; for (uint s = 0; s < bin_cent.size(); ++s){ bin_edges.push_back(bin_cent[s] - bin_stride/2.); } bin_edges.push_back(bin_cent[bin_cent.size()-1]+bin_stride/2.); // Finally, make the histogram TH1D* flux = new TH1D("flux", "flux; E_{#nu} (GeV)", bin_cent.size(), &bin_edges[0]); for (uint x = 0; x < bin_cent.size(); ++x){ flux -> SetBinContent(x+1, bin_vals[x]); } return flux; } void CreateRateHistogram(std::string inputList, std::string flux, std::string output) { TChain *tn = new TChain("RootTuple"); double E, xsec; tn->SetBranchAddress("lepIn_E", &E); tn->SetBranchAddress("weight", &xsec); std::vector inputs = GeneralUtils::ParseToStr(inputList, ","); for (std::vector::iterator it = inputs.begin(); it != inputs.end(); ++it) { NUIS_LOG(FIT, "Adding " << *it << " to the output"); tn->AddFile((*it).c_str()); } if (inputs.size() > 1 && output.empty()) { NUIS_ABORT("You must provide a new output file name if you want to have " "more than 1 input file!"); } int nevts = tn->GetEntries(); if (!nevts) { NUIS_ABORT("Either the input file is not from GiBUU, or it's empty..."); } // Get flux hist std::vector fluxvect = GeneralUtils::ParseToStr(flux, ","); TH1D *fluxHist = NULL; if (fluxvect.size() > 1) { TFile *fluxfile = new TFile(fluxvect[0].c_str(), "READ"); fluxHist = (TH1D *)fluxfile->Get(fluxvect[1].c_str()); fluxHist->SetDirectory(0); } else if (flux.find(".dat") != std::string::npos){ fluxHist = MakeFluxHistFromDatFile(flux); } else { NUIS_ABORT("NO FLUX SPECIFIED"); } // Make Event Hist TH1D *xsecHist = (TH1D *)fluxHist->Clone(); xsecHist->Reset(); // Make a total cross section hist for shits and giggles TH1D *entryHist = (TH1D *)xsecHist->Clone(); for (int i = 0; i < nevts; ++i) { tn->GetEntry(i); xsecHist->Fill(E, xsec); entryHist->Fill(E); if (i % (nevts / 10) == 0) { NUIS_LOG(FIT, "Processed " << i << "/" << nevts << " GiBUU events." << "(Enu = " << E << ", xsec = " << xsec << ") "); } } NUIS_LOG(FIT, "Processed all events"); - - xsecHist->Scale(1, "width"); - // xsecHist->Divide(entryHist); + + // Somewhat annoyingly, have to include the flux width! + double flux_range = (xsecHist->GetXaxis()->GetBinUpEdge(xsecHist->GetNbinsX()+1) - \ + xsecHist->GetXaxis()->GetBinLowEdge(1)); + xsecHist->Scale(flux_range, "width"); // This will be the evtrt histogram TH1D *evtHist = (TH1D*)xsecHist->Clone(); evtHist->Multiply(fluxHist); // Check whether the overflow is empty. If not, advise that either the wrong // flux histogram or units were used... // If the events were generated with a limited range of the flux histogram, // this may be benign if (evtHist->Integral(0, -1) != evtHist->Integral() || evtHist->Integral(0, -1) == 0) { NUIS_ERR(WRN, "The input file(" << evtHist->Integral(0, -1) << ") and flux histogram provided do not match... "); NUIS_ERR(WRN, "Are the units correct? Did you provide the correct flux file?"); NUIS_ERR(WRN, "Use output with caution..."); } // Pick where the output should go TFile *outFile = NULL; if (!output.empty()) { NUIS_LOG(FIT, "Saving histograms in " << output); outFile = new TFile(output.c_str(), "RECREATE"); } else { NUIS_LOG(FIT, "Saving histograms in " << inputs[0]); outFile = new TFile(inputs[0].c_str(), "UPDATE"); } outFile->cd(); std::string xsec_name = "xsec_PrepareGiBUU"; std::string flux_name = "flux_PrepareGiBUU"; std::string rate_name = "evtrt_PrepareGiBUU"; if (output.empty()) { // Check whether we should overwrite existing histograms std::string input_xsec = PlotUtils::GetObjectWithName(outFile, "xsec"); std::string input_flux = PlotUtils::GetObjectWithName(outFile, "flux"); std::string input_rate = PlotUtils::GetObjectWithName(outFile, "evtrt"); if (!input_xsec.empty()) { NUIS_LOG(FIT, "Updating histogram: " << input_xsec); xsec_name = input_xsec; } if (!input_flux.empty()) { NUIS_LOG(FIT, "Updating histogram: " << input_flux); flux_name = input_flux; } if (!input_rate.empty()) { NUIS_LOG(FIT, "Updating histogram: " << input_rate); rate_name = input_rate; } } else { NUIS_LOG(FIT, "Cloning event tree into output file."); StopTalking(); TTree *newtree = (TTree *)tn->CloneTree(-1, "fast"); StartTalking(); newtree->Write(); } xsecHist->Write(xsec_name.c_str(), TObject::kOverwrite); fluxHist->Write(flux_name.c_str(), TObject::kOverwrite); evtHist->Write(rate_name.c_str(), TObject::kOverwrite); outFile->Close(); return; } void PrintOptions() { std::cout << "PrepareGiBUU NUISANCE app. " << std::endl << "Produces or recalculates evtrt and flux histograms necessary " "for NUISANCE normalization." << std::endl; std::cout << "PrepareGiBUU: " << std::endl; std::cout << " [-h,-help,--h,--help]" << std::endl; std::cout << " -i inputfile1.root,inputfile2.root,inputfile3.root,..." << std::endl; std::cout << " Takes any number of files, but assumes all are " "produced with a single flux" << std::endl; std::cout << " -f flux_root_file.root,flux_hist_name" << std::endl; std::cout << " Path to ROOT or .dat file containing the flux used " "when generating the GiBUU files" << std::endl; std::cout << " [-o outputfile.root] " << std::endl; std::cout << " If an output file is not given, the input file will be used" << std::endl; std::cout << " If more than one input file is given, an output file " "must be given" << std::endl; } void ParseOptions(int argc, char *argv[]) { bool flagopt = false; // If No Arguments print commands for (int i = 1; i < argc; ++i) { if (!std::strcmp(argv[i], "-h")) { flagopt = true; break; } if (i + 1 != argc) { // Cardfile if (!std::strcmp(argv[i], "-h")) { flagopt = true; break; } else if (!std::strcmp(argv[i], "-i")) { fInputFiles = argv[i + 1]; ++i; } else if (!std::strcmp(argv[i], "-o")) { fOutputFile = argv[i + 1]; ++i; } else if (!std::strcmp(argv[i], "-f")) { fFluxFile = argv[i + 1]; ++i; } else { NUIS_ERR(FTL, "ERROR: unknown command line option given! - '" << argv[i] << " " << argv[i + 1] << "'"); PrintOptions(); break; } } } if (fInputFiles == "" && !flagopt) { NUIS_ERR(FTL, "No input file(s) specified!"); flagopt = true; } if (fFluxFile == "" && (!flagopt)) { NUIS_ERR(FTL, "No flux input specified!"); flagopt = true; } if (argc < 1 || flagopt) { PrintOptions(); exit(-1); } return; } diff --git a/data/T2K/CCCOH/C12_Enu_1bin.txt b/data/T2K/CCCOH/C12_Enu_1bin.txt new file mode 100644 index 0000000..5982bdb --- /dev/null +++ b/data/T2K/CCCOH/C12_Enu_1bin.txt @@ -0,0 +1,2 @@ +0.5 2.9E-40 1.304E-40 +1.5 0 0 diff --git a/src/FCN/SampleList.cxx b/src/FCN/SampleList.cxx index 6b28a5b..d8d1eb6 100644 --- a/src/FCN/SampleList.cxx +++ b/src/FCN/SampleList.cxx @@ -1,1576 +1,1582 @@ #include "SampleList.h" #ifndef __NO_ANL__ #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_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" #include "ANL_CC1ppip_Evt_1DWNpi_nu.h" #include "ANL_CC1ppip_Evt_1DWNmu_nu.h" #include "ANL_CC1ppip_Evt_1DWmupi_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" #include "ANL_CC1npip_Evt_1DWNpi_nu.h" #include "ANL_CC1npip_Evt_1DWNmu_nu.h" #include "ANL_CC1npip_Evt_1DWmupi_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" #include "ANL_CC1pi0_Evt_1DWNpi_nu.h" #include "ANL_CC1pi0_Evt_1DWNmu_nu.h" #include "ANL_CC1pi0_Evt_1DWmupi_nu.h" // ANL NC1npip (mm, exotic!) #include "ANL_NC1npip_Evt_1Dppi_nu.h" // ANL NC1ppim (mm, exotic!) #include "ANL_NC1ppim_Evt_1DcosmuStar_nu.h" #include "ANL_NC1ppim_XSec_1DEnu_nu.h" // ANL CC2pi 1pim1pip (mm, even more exotic!) #include "ANL_CC2pi_1pim1pip_Evt_1Dpmu_nu.h" #include "ANL_CC2pi_1pim1pip_Evt_1Dppim_nu.h" #include "ANL_CC2pi_1pim1pip_Evt_1Dppip_nu.h" #include "ANL_CC2pi_1pim1pip_Evt_1Dpprot_nu.h" #include "ANL_CC2pi_1pim1pip_XSec_1DEnu_nu.h" // ANL CC2pi 1pip1pip (mm, even more exotic!) #include "ANL_CC2pi_1pip1pip_Evt_1Dpmu_nu.h" #include "ANL_CC2pi_1pip1pip_Evt_1Dpneut_nu.h" #include "ANL_CC2pi_1pip1pip_Evt_1DppipHigh_nu.h" #include "ANL_CC2pi_1pip1pip_Evt_1DppipLow_nu.h" #include "ANL_CC2pi_1pip1pip_XSec_1DEnu_nu.h" // ANL CC2pi 1pip1pi0 (mm, even more exotic!) #include "ANL_CC2pi_1pip1pi0_Evt_1Dpmu_nu.h" #include "ANL_CC2pi_1pip1pi0_Evt_1Dppi0_nu.h" #include "ANL_CC2pi_1pip1pi0_Evt_1Dppip_nu.h" #include "ANL_CC2pi_1pip1pi0_Evt_1Dpprot_nu.h" #include "ANL_CC2pi_1pip1pi0_XSec_1DEnu_nu.h" #endif #ifndef __NO_ArgoNeuT__ // ArgoNeuT CC1Pi #include "ArgoNeuT_CC1Pi_XSec_1Dpmu_antinu.h" #include "ArgoNeuT_CC1Pi_XSec_1Dpmu_nu.h" #include "ArgoNeuT_CC1Pi_XSec_1Dthetamu_antinu.h" #include "ArgoNeuT_CC1Pi_XSec_1Dthetamu_nu.h" #include "ArgoNeuT_CC1Pi_XSec_1Dthetamupi_antinu.h" #include "ArgoNeuT_CC1Pi_XSec_1Dthetamupi_nu.h" #include "ArgoNeuT_CC1Pi_XSec_1Dthetapi_antinu.h" #include "ArgoNeuT_CC1Pi_XSec_1Dthetapi_nu.h" // ArgoNeuT CC-inclusive #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" #endif #ifndef __NO_BNL__ // 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_1DcosthAdler_nu.h" #include "BNL_CC1ppip_Evt_1Dphi_nu.h" #include "BNL_CC1ppip_XSec_1DEnu_nu.h" #include "BNL_CC1ppip_Evt_1DWNpi_nu.h" #include "BNL_CC1ppip_Evt_1DWNmu_nu.h" #include "BNL_CC1ppip_Evt_1DWmupi_nu.h" // BNL CC1npip #include "BNL_CC1npip_Evt_1DQ2_nu.h" #include "BNL_CC1npip_XSec_1DEnu_nu.h" #include "BNL_CC1npip_Evt_1DWNpi_nu.h" #include "BNL_CC1npip_Evt_1DWNmu_nu.h" #include "BNL_CC1npip_Evt_1DWmupi_nu.h" // BNL CC1pi0 #include "BNL_CC1pi0_Evt_1DQ2_nu.h" #include "BNL_CC1pi0_XSec_1DEnu_nu.h" #include "BNL_CC1pi0_Evt_1DWNpi_nu.h" #include "BNL_CC1pi0_Evt_1DWNmu_nu.h" #include "BNL_CC1pi0_Evt_1DWmupi_nu.h" // BNL multipi #include "BNL_CC2pi_1pim1pip_XSec_1DEnu_nu.cxx" #include "BNL_CC3pi_1pim2pip_XSec_1DEnu_nu.cxx" #include "BNL_CC4pi_2pim2pip_XSec_1DEnu_nu.cxx" #include "BNL_CC2pi_1pim1pip_Evt_1DWpippim_nu.cxx" #include "BNL_CC2pi_1pim1pip_Evt_1DWpippr_nu.cxx" #endif #ifndef __NO_FNAL__ // 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" #endif #ifndef __NO_BEBC__ // 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" #endif #ifndef __NO_GGM__ // GGM CC1ppip #include "GGM_CC1ppip_Evt_1DQ2_nu.h" #include "GGM_CC1ppip_XSec_1DEnu_nu.h" #endif #ifndef __NO_MiniBooNE__ // MiniBooNE CCQE #include "MiniBooNE_CCQE_XSec_1DEnu_nu.h" #include "MiniBooNE_CCQE_XSec_1DQ2_antinu.h" #include "MiniBooNE_CCQE_XSec_1DQ2_nu.h" #include "MiniBooNE_CCQE_XSec_2DTcos_antinu.h" #include "MiniBooNE_CCQE_XSec_2DTcos_nu.h" // MiniBooNE CC1pi+ 1D #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" // MiniBooNE CC1pi+ 2D #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_NC1pi0_XSec_1Dcospi0_antinu.h" #include "MiniBooNE_NC1pi0_XSec_1Dcospi0_nu.h" #include "MiniBooNE_NC1pi0_XSec_1Dppi0_antinu.h" #include "MiniBooNE_NC1pi0_XSec_1Dppi0_nu.h" // MiniBooNE NCEL #include "MiniBooNE_NCEL_XSec_Treco_nu.h" #endif #ifndef __NO_MicroBooNE__ #include "MicroBooNE_CCInc_XSec_2DPcos_nu.h" #include "MicroBooNE_CC1MuNp_XSec_1D_nu.h" #endif #ifndef __NO_MINERvA__ // 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" // 2018 MINERvA CC0pi STV #include "MINERvA_CC0pinp_STV_XSec_1D_nu.h" // 2018 MINERvA CC0pi 2D #include "MINERvA_CC0pi_XSec_1D_2018_nu.h" #include "MINERvA_CC0pi_XSec_2D_nu.h" // #include "MINERvA_CC0pi_XSec_3DptpzTp_nu.h" // 2018 MINERvA CC0pi 2D antinu #include "MINERvA_CC0pi_XSec_2D_antinu.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" // 2017 data update #include "MINERvA_CC1pip_XSec_1D_2017Update.h" // MINERvA CC1pi- #include "MINERvA_CC1pim_XSec_1DEnu_antinu.h" #include "MINERvA_CC1pim_XSec_1DQ2_antinu.h" #include "MINERvA_CC1pim_XSec_1DTpi_antinu.h" #include "MINERvA_CC1pim_XSec_1Dpmu_antinu.h" #include "MINERvA_CC1pim_XSec_1Dth_antinu.h" #include "MINERvA_CC1pim_XSec_1Dthmu_antinu.h" // MINERvA CCNpi+ #include "MINERvA_CCNpip_XSec_1DEnu_nu.h" #include "MINERvA_CCNpip_XSec_1DQ2_nu.h" #include "MINERvA_CCNpip_XSec_1DTpi_nu.h" #include "MINERvA_CCNpip_XSec_1Dpmu_nu.h" #include "MINERvA_CCNpip_XSec_1Dth_nu.h" #include "MINERvA_CCNpip_XSec_1Dthmu_nu.h" // MINERvA CC1pi0 #include "MINERvA_CC1pi0_XSec_1DEnu_antinu.h" #include "MINERvA_CC1pi0_XSec_1DQ2_antinu.h" #include "MINERvA_CC1pi0_XSec_1DTpi0_antinu.h" #include "MINERvA_CC1pi0_XSec_1Dpmu_antinu.h" #include "MINERvA_CC1pi0_XSec_1Dppi0_antinu.h" #include "MINERvA_CC1pi0_XSec_1Dth_antinu.h" #include "MINERvA_CC1pi0_XSec_1Dthmu_antinu.h" // MINERvA CC1pi0 neutrino #include "MINERvA_CC1pi0_XSec_1D_nu.h" // MINERvA CCINC #include "MINERvA_CCinc_XSec_1DEnu_ratio.h" #include "MINERvA_CCinc_XSec_1Dx_ratio.h" #include "MINERvA_CCinc_XSec_2DEavq3_nu.h" // MINERvA CCDIS #include "MINERvA_CCDIS_XSec_1DEnu_ratio.h" #include "MINERvA_CCDIS_XSec_1Dx_ratio.h" // MINERvA CCCOH pion #include "MINERvA_CCCOHPI_XSec_1DEnu_antinu.h" #include "MINERvA_CCCOHPI_XSec_1DEpi_antinu.h" #include "MINERvA_CCCOHPI_XSec_1DQ2_antinu.h" #include "MINERvA_CCCOHPI_XSec_1DEpi_nu.h" #include "MINERvA_CCCOHPI_XSec_1DQ2_nu.h" #include "MINERvA_CCCOHPI_XSec_1Dth_nu.h" #include "MINERvA_CCCOHPI_XSec_joint.h" #include "MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu.h" #include "MINERvA_CC0pi_XSec_1DQ2_Tgt_nu.h" #endif #ifndef __NO_T2K__ // T2K CC0pi 2016 #include "T2K_CC0pi_XSec_2DPcos_nu_I.h" #include "T2K_CC0pi_XSec_2DPcos_nu_II.h" // T2K CC0pi 2020 arXiv:1908.10249 #include "T2K_CC0pi_XSec_H2O_2DPcos_anu.h" // T2K CC0pi 2020 arXiv:2004.05434 #include "T2K_NuMu_CC0pi_OC_XSec_2DPcos.h" #include "T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint.h" // T2K CC0pi 2020 arXiv:2002.09323 #include "T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos.h" #include "T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint.h" // T2K CC-inclusive with full acceptance 2018 #include "T2K_CCinc_XSec_2DPcos_nu_nonuniform.h" // T2K nue CC-inclusive 2019 #include "T2K_nueCCinc_XSec_1Dpe.h" #include "T2K_nueCCinc_XSec_1Dthe.h" #include "T2K_nueCCinc_XSec_1Dpe_joint.h" #include "T2K_nueCCinc_XSec_1Dthe_joint.h" #include "T2K_nueCCinc_XSec_joint.h" // T2K STV CC0pi 2018 #include "T2K_CC0piWithProtons_XSec_2018_multidif_0p_1p_Np.h" #include "T2K_CC0pinp_STV_XSec_1Ddat_nu.h" #include "T2K_CC0pinp_STV_XSec_1Ddphit_nu.h" #include "T2K_CC0pinp_STV_XSec_1Ddpt_nu.h" #include "T2K_CC0pinp_ifk_XSec_3Dinfa_nu.h" #include "T2K_CC0pinp_ifk_XSec_3Dinfip_nu.h" #include "T2K_CC0pinp_ifk_XSec_3Dinfp_nu.h" // T2K CC1pi+ on CH #include "T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.h" #include "T2K_CC1pip_CH_XSec_1DCosThAdler_nu.h" #include "T2K_CC1pip_CH_XSec_1DQ2_nu.h" #include "T2K_CC1pip_CH_XSec_1Dppi_nu.h" #include "T2K_CC1pip_CH_XSec_1Dthmupi_nu.h" #include "T2K_CC1pip_CH_XSec_1Dthpi_nu.h" #include "T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.h" +// T2K CCCOH (single bin) +#include "T2K_CCCOH_C12_XSec_1DEnu_nu.h" + // T2K CC1pi+ on H2O #include "T2K_CC1pip_H2O_XSec_1DEnuDelta_nu.h" #include "T2K_CC1pip_H2O_XSec_1DEnuMB_nu.h" #include "T2K_CC1pip_H2O_XSec_1Dcosmu_nu.h" #include "T2K_CC1pip_H2O_XSec_1Dcosmupi_nu.h" #include "T2K_CC1pip_H2O_XSec_1Dcospi_nu.h" #include "T2K_CC1pip_H2O_XSec_1Dpmu_nu.h" #include "T2K_CC1pip_H2O_XSec_1Dppi_nu.h" // add header here #endif #ifndef __NO_SciBooNE__ // SciBooNE COH studies #include "SciBooNE_CCCOH_1TRK_1DQ2_nu.h" #include "SciBooNE_CCCOH_1TRK_1Dpmu_nu.h" #include "SciBooNE_CCCOH_1TRK_1Dthetamu_nu.h" #include "SciBooNE_CCCOH_MuPiNoVA_1DQ2_nu.h" #include "SciBooNE_CCCOH_MuPiNoVA_1Dpmu_nu.h" #include "SciBooNE_CCCOH_MuPiNoVA_1Dthetamu_nu.h" #include "SciBooNE_CCCOH_MuPiNoVA_1Dthetapi_nu.h" #include "SciBooNE_CCCOH_MuPiNoVA_1Dthetapr_nu.h" #include "SciBooNE_CCCOH_MuPiVA_1DQ2_nu.h" #include "SciBooNE_CCCOH_MuPiVA_1Dpmu_nu.h" #include "SciBooNE_CCCOH_MuPiVA_1Dthetamu_nu.h" #include "SciBooNE_CCCOH_MuPr_1DQ2_nu.h" #include "SciBooNE_CCCOH_MuPr_1Dpmu_nu.h" #include "SciBooNE_CCCOH_MuPr_1Dthetamu_nu.h" #include "SciBooNE_CCCOH_STOPFINAL_1DQ2_nu.h" #include "SciBooNE_CCCOH_STOP_NTrks_nu.h" #include "SciBooNE_CCInc_XSec_1DEnu_nu.h" #endif #ifndef __NO_K2K__ // K2K NC1pi0 #include "K2K_NC1pi0_Evt_1Dppi0_nu.h" #endif // MC Studies #include "ExpMultDist_CCQE_XSec_1DVar_FakeStudy.h" #include "ExpMultDist_CCQE_XSec_2DVar_FakeStudy.h" #include "MCStudy_CCQEHistograms.h" #include "GenericFlux_Tester.h" #include "GenericFlux_Vectors.h" #include "ElectronFlux_FlatTree.h" #include "ElectronScattering_DurhamData.h" #include "MCStudy_KaonPreSelection.h" #include "MCStudy_MuonValidation.h" #include "OfficialNIWGPlots.h" #include "T2K2017_FakeData.h" #include "SigmaEnuHists.h" #include "Simple_Osc.h" #include "Smear_SVDUnfold_Propagation_Osc.h" #include "FitWeight.h" #include "NuisConfig.h" #include "NuisKey.h" #ifdef __USE_DYNSAMPLES__ #include "TRegexp.h" #include // linux #include DynamicSampleFactory::DynamicSampleFactory() : NSamples(0), NManifests(0) { LoadPlugins(); NUIS_LOG(FIT, "Loaded " << NSamples << " from " << NManifests << " shared object libraries."); } DynamicSampleFactory *DynamicSampleFactory::glblDSF = NULL; DynamicSampleFactory::PluginManifest::~PluginManifest() { for (size_t i_it = 0; i_it < Instances.size(); ++i_it) { (*(DSF_DestroySample))(Instances[i_it]); } } std::string EnsureTrailingSlash(std::string const &inp) { if (!inp.length()) { return "/"; } if (inp[inp.length() - 1] == '/') { return inp; } return inp + "/"; } void DynamicSampleFactory::LoadPlugins() { std::vector SearchDirectories; if (Config::HasPar("dynamic_sample.path")) { SearchDirectories = GeneralUtils::ParseToStr(Config::GetParS("dynamic_sample.path"), ":"); } char const *envPath = getenv("NUISANCE_DS_PATH"); if (envPath) { std::vector envPaths = GeneralUtils::ParseToStr(envPath, ":"); for (size_t ep_it = 0; ep_it < envPaths.size(); ++ep_it) { SearchDirectories.push_back(envPaths[ep_it]); } } if (!SearchDirectories.size()) { char const *pwdPath = getenv("PWD"); if (pwdPath) { SearchDirectories.push_back(pwdPath); } } for (size_t sp_it = 0; sp_it < SearchDirectories.size(); ++sp_it) { std::string dirpath = EnsureTrailingSlash(SearchDirectories[sp_it]); NUIS_LOG(FIT, "Searching for dynamic sample manifests in: " << dirpath); Ssiz_t len = 0; DIR *dir; struct dirent *ent; dir = opendir(dirpath.c_str()); if (dir != NULL) { TRegexp matchExp("*.so", true); while ((ent = readdir(dir)) != NULL) { if (matchExp.Index(TString(ent->d_name), &len) != Ssiz_t(-1)) { NUIS_LOG(FIT, "\tFound shared object: " << ent->d_name << " checking for relevant methods..."); void *dlobj = dlopen((dirpath + ent->d_name).c_str(), RTLD_NOW | RTLD_GLOBAL); char const *dlerr_cstr = dlerror(); std::string dlerr; if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { NUIS_ERR(WRN, "\tDL Load Error: " << dlerr); continue; } PluginManifest plgManif; plgManif.dllib = dlobj; plgManif.soloc = (dirpath + ent->d_name); plgManif.DSF_NSamples = reinterpret_cast(dlsym(dlobj, "DSF_NSamples")); dlerr = ""; dlerr_cstr = dlerror(); if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { NUIS_ERR(WRN, "\tFailed to load symbol \"DSF_NSamples\" from " << (dirpath + ent->d_name) << ": " << dlerr); dlclose(dlobj); continue; } plgManif.DSF_GetSampleName = reinterpret_cast( dlsym(dlobj, "DSF_GetSampleName")); dlerr = ""; dlerr_cstr = dlerror(); if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { NUIS_ERR(WRN, "\tFailed to load symbol \"DSF_GetSampleName\" from " << (dirpath + ent->d_name) << ": " << dlerr); dlclose(dlobj); continue; } plgManif.DSF_GetSample = reinterpret_cast( dlsym(dlobj, "DSF_GetSample")); dlerr = ""; dlerr_cstr = dlerror(); if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { NUIS_ERR(WRN, "\tFailed to load symbol \"DSF_GetSample\" from " << (dirpath + ent->d_name) << ": " << dlerr); dlclose(dlobj); continue; } plgManif.DSF_DestroySample = reinterpret_cast( dlsym(dlobj, "DSF_DestroySample")); dlerr = ""; dlerr_cstr = dlerror(); if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { NUIS_ERR(WRN, "Failed to load symbol \"DSF_DestroySample\" from " << (dirpath + ent->d_name) << ": " << dlerr); dlclose(dlobj); continue; } plgManif.NSamples = (*(plgManif.DSF_NSamples))(); NUIS_LOG(FIT, "\tSuccessfully loaded dynamic sample manifest: " << plgManif.soloc << ". Contains " << plgManif.NSamples << " samples."); for (size_t smp_it = 0; smp_it < plgManif.NSamples; ++smp_it) { char const *smp_name = (*(plgManif.DSF_GetSampleName))(smp_it); if (!smp_name) { NUIS_ABORT("Could not load sample " << smp_it << " / " << plgManif.NSamples << " from " << plgManif.soloc); } if (Samples.count(smp_name)) { NUIS_ERR(WRN, "Already loaded a sample named: \"" << smp_name << "\". cannot load duplciates. This " "sample will be skipped."); continue; } plgManif.SamplesProvided.push_back(smp_name); Samples[smp_name] = std::make_pair(plgManif.soloc, smp_it); NUIS_LOG(FIT, "\t\t" << smp_name); } if (plgManif.SamplesProvided.size()) { Manifests[plgManif.soloc] = plgManif; NSamples += plgManif.SamplesProvided.size(); NManifests++; } else { dlclose(dlobj); } } } closedir(dir); } else { NUIS_ERR(WRN, "Tried to open non-existant directory."); } } } DynamicSampleFactory &DynamicSampleFactory::Get() { if (!glblDSF) { glblDSF = new DynamicSampleFactory(); } return *glblDSF; } void DynamicSampleFactory::Print() { std::map > ManifestSamples; for (std::map >::iterator smp_it = Samples.begin(); smp_it != Samples.end(); ++smp_it) { if (!ManifestSamples.count(smp_it->second.first)) { ManifestSamples[smp_it->second.first] = std::vector(); } ManifestSamples[smp_it->second.first].push_back(smp_it->first); } NUIS_LOG(FIT, "Dynamic sample manifest: "); for (std::map >::iterator m_it = ManifestSamples.begin(); m_it != ManifestSamples.end(); ++m_it) { NUIS_LOG(FIT, "\tLibrary " << m_it->first << " contains: "); for (size_t s_it = 0; s_it < m_it->second.size(); ++s_it) { NUIS_LOG(FIT, "\t\t" << m_it->second[s_it]); } } } bool DynamicSampleFactory::HasSample(std::string const &name) { return Samples.count(name); } bool DynamicSampleFactory::HasSample(nuiskey &samplekey) { return HasSample(samplekey.GetS("name")); } MeasurementBase *DynamicSampleFactory::CreateSample(nuiskey &samplekey) { if (!HasSample(samplekey)) { NUIS_ERR(WRN, "Asked to load unknown sample: \"" << samplekey.GetS("name") << "\"."); return NULL; } std::pair sample = Samples[samplekey.GetS("name")]; NUIS_LOG(SAM, "\tLoading sample " << sample.second << " from " << sample.first); return (*(Manifests[sample.first].DSF_GetSample))(sample.second, &samplekey); } DynamicSampleFactory::~DynamicSampleFactory() { Manifests.clear(); } #endif //! 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. MeasurementBase *CreateSample(std::string name, std::string file, std::string type, std::string fkdt, FitWeight *rw) { nuiskey samplekey = Config::CreateKey("sample"); samplekey.Set("name", name); samplekey.Set("input", file); samplekey.Set("type", type); return CreateSample(samplekey); } MeasurementBase *CreateSample(nuiskey samplekey) { #ifdef __USE_DYNSAMPLES__ if (DynamicSampleFactory::Get().HasSample(samplekey)) { NUIS_LOG(SAM, "Instantiating dynamic sample..."); MeasurementBase *ds = DynamicSampleFactory::Get().CreateSample(samplekey); if (ds) { NUIS_LOG(SAM, "Done."); return ds; } NUIS_ABORT("Failed to instantiate dynamic sample."); } #endif FitWeight *rw = FitBase::GetRW(); std::string name = samplekey.GetS("name"); std::string file = samplekey.GetS("input"); std::string type = samplekey.GetS("type"); std::string fkdt = ""; /* ANL CCQE Samples */ #ifndef __NO_ANL__ if (!name.compare("ANL_CCQE_XSec_1DEnu_nu") || !name.compare("ANL_CCQE_XSec_1DEnu_nu_PRD26") || !name.compare("ANL_CCQE_XSec_1DEnu_nu_PRL31") || !name.compare("ANL_CCQE_XSec_1DEnu_nu_PRD16")) { return (new ANL_CCQE_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CCQE_Evt_1DQ2_nu") || !name.compare("ANL_CCQE_Evt_1DQ2_nu_PRL31") || !name.compare("ANL_CCQE_Evt_1DQ2_nu_PRD26") || !name.compare("ANL_CCQE_Evt_1DQ2_nu_PRD16")) { return (new ANL_CCQE_Evt_1DQ2_nu(samplekey)); /* ANL CC1ppip samples */ } else if (!name.compare("ANL_CC1ppip_XSec_1DEnu_nu") || !name.compare("ANL_CC1ppip_XSec_1DEnu_nu_W14Cut") || !name.compare("ANL_CC1ppip_XSec_1DEnu_nu_Uncorr") || !name.compare("ANL_CC1ppip_XSec_1DEnu_nu_W14Cut_Uncorr") || !name.compare("ANL_CC1ppip_XSec_1DEnu_nu_W16Cut_Uncorr")) { return (new ANL_CC1ppip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_XSec_1DQ2_nu")) { return (new ANL_CC1ppip_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1DQ2_nu") || !name.compare("ANL_CC1ppip_Evt_1DQ2_nu_W14Cut")) { return (new ANL_CC1ppip_Evt_1DQ2_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1Dppi_nu")) { return (new ANL_CC1ppip_Evt_1Dppi_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1Dthpr_nu")) { return (new ANL_CC1ppip_Evt_1Dthpr_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1DcosmuStar_nu")) { return (new ANL_CC1ppip_Evt_1DcosmuStar_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1DcosthAdler_nu")) { return (new ANL_CC1ppip_Evt_1DcosthAdler_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1Dphi_nu")) { return (new ANL_CC1ppip_Evt_1Dphi_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1DWNpi_nu")) { return (new ANL_CC1ppip_Evt_1DWNpi_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1DWNmu_nu")) { return (new ANL_CC1ppip_Evt_1DWNmu_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1DWmupi_nu")) { return (new ANL_CC1ppip_Evt_1DWmupi_nu(samplekey)); /* ANL CC1npip sample */ } else if (!name.compare("ANL_CC1npip_XSec_1DEnu_nu") || !name.compare("ANL_CC1npip_XSec_1DEnu_nu_W14Cut") || !name.compare("ANL_CC1npip_XSec_1DEnu_nu_Uncorr") || !name.compare("ANL_CC1npip_XSec_1DEnu_nu_W14Cut_Uncorr") || !name.compare("ANL_CC1npip_XSec_1DEnu_nu_W16Cut_Uncorr")) { return (new ANL_CC1npip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC1npip_Evt_1DQ2_nu") || !name.compare("ANL_CC1npip_Evt_1DQ2_nu_W14Cut")) { return (new ANL_CC1npip_Evt_1DQ2_nu(samplekey)); } else if (!name.compare("ANL_CC1npip_Evt_1Dppi_nu")) { return (new ANL_CC1npip_Evt_1Dppi_nu(samplekey)); } else if (!name.compare("ANL_CC1npip_Evt_1DcosmuStar_nu")) { return (new ANL_CC1npip_Evt_1DcosmuStar_nu(samplekey)); } else if (!name.compare("ANL_CC1npip_Evt_1DWNpi_nu")) { return (new ANL_CC1npip_Evt_1DWNpi_nu(samplekey)); } else if (!name.compare("ANL_CC1npip_Evt_1DWNmu_nu")) { return (new ANL_CC1npip_Evt_1DWNmu_nu(samplekey)); } else if (!name.compare("ANL_CC1npip_Evt_1DWmupi_nu")) { return (new ANL_CC1npip_Evt_1DWmupi_nu(samplekey)); /* ANL CC1pi0 sample */ } else if (!name.compare("ANL_CC1pi0_XSec_1DEnu_nu") || !name.compare("ANL_CC1pi0_XSec_1DEnu_nu_W14Cut") || !name.compare("ANL_CC1pi0_XSec_1DEnu_nu_Uncorr") || !name.compare("ANL_CC1pi0_XSec_1DEnu_nu_W14Cut_Uncorr") || !name.compare("ANL_CC1pi0_XSec_1DEnu_nu_W16Cut_Uncorr")) { return (new ANL_CC1pi0_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC1pi0_Evt_1DQ2_nu") || !name.compare("ANL_CC1pi0_Evt_1DQ2_nu_W14Cut")) { return (new ANL_CC1pi0_Evt_1DQ2_nu(samplekey)); } else if (!name.compare("ANL_CC1pi0_Evt_1DcosmuStar_nu")) { return (new ANL_CC1pi0_Evt_1DcosmuStar_nu(samplekey)); } else if (!name.compare("ANL_CC1pi0_Evt_1DWNpi_nu")) { return (new ANL_CC1pi0_Evt_1DWNpi_nu(samplekey)); } else if (!name.compare("ANL_CC1pi0_Evt_1DWNmu_nu")) { return (new ANL_CC1pi0_Evt_1DWNmu_nu(samplekey)); } else if (!name.compare("ANL_CC1pi0_Evt_1DWmupi_nu")) { return (new ANL_CC1pi0_Evt_1DWmupi_nu(samplekey)); /* ANL NC1npip sample */ } else if (!name.compare("ANL_NC1npip_Evt_1Dppi_nu")) { return (new ANL_NC1npip_Evt_1Dppi_nu(samplekey)); /* ANL NC1ppim sample */ } else if (!name.compare("ANL_NC1ppim_XSec_1DEnu_nu")) { return (new ANL_NC1ppim_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_NC1ppim_Evt_1DcosmuStar_nu")) { return (new ANL_NC1ppim_Evt_1DcosmuStar_nu(samplekey)); /* ANL CC2pi sample */ } else if (!name.compare("ANL_CC2pi_1pim1pip_XSec_1DEnu_nu")) { return (new ANL_CC2pi_1pim1pip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dpmu_nu")) { return (new ANL_CC2pi_1pim1pip_Evt_1Dpmu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dppip_nu")) { return (new ANL_CC2pi_1pim1pip_Evt_1Dppip_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dppim_nu")) { return (new ANL_CC2pi_1pim1pip_Evt_1Dppim_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dpprot_nu")) { return (new ANL_CC2pi_1pim1pip_Evt_1Dpprot_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_XSec_1DEnu_nu")) { return (new ANL_CC2pi_1pip1pip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1Dpmu_nu")) { return (new ANL_CC2pi_1pip1pip_Evt_1Dpmu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1Dpneut_nu")) { return (new ANL_CC2pi_1pip1pip_Evt_1Dpneut_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1DppipHigh_nu")) { return (new ANL_CC2pi_1pip1pip_Evt_1DppipHigh_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1DppipLow_nu")) { return (new ANL_CC2pi_1pip1pip_Evt_1DppipLow_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_XSec_1DEnu_nu")) { return (new ANL_CC2pi_1pip1pi0_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dpmu_nu")) { return (new ANL_CC2pi_1pip1pi0_Evt_1Dpmu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dppip_nu")) { return (new ANL_CC2pi_1pip1pi0_Evt_1Dppip_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dppi0_nu")) { return (new ANL_CC2pi_1pip1pi0_Evt_1Dppi0_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dpprot_nu")) { return (new ANL_CC2pi_1pip1pi0_Evt_1Dpprot_nu(samplekey)); /* ArgoNeut Samples */ } else #endif #ifndef __NO_ArgoNeuT__ if (!name.compare("ArgoNeuT_CCInc_XSec_1Dpmu_antinu")) { return (new ArgoNeuT_CCInc_XSec_1Dpmu_antinu(samplekey)); } else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dpmu_nu")) { return (new ArgoNeuT_CCInc_XSec_1Dpmu_nu(samplekey)); } else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dthetamu_antinu")) { return (new ArgoNeuT_CCInc_XSec_1Dthetamu_antinu(samplekey)); } else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dthetamu_nu")) { return (new ArgoNeuT_CCInc_XSec_1Dthetamu_nu(samplekey)); } else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dpmu_nu")) { return (new ArgoNeuT_CC1Pi_XSec_1Dpmu_nu(samplekey)); } else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetamu_nu")) { return (new ArgoNeuT_CC1Pi_XSec_1Dthetamu_nu(samplekey)); } else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetapi_nu")) { return (new ArgoNeuT_CC1Pi_XSec_1Dthetapi_nu(samplekey)); } else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetamupi_nu")) { return (new ArgoNeuT_CC1Pi_XSec_1Dthetamupi_nu(samplekey)); } else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dpmu_antinu")) { return (new ArgoNeuT_CC1Pi_XSec_1Dpmu_antinu(samplekey)); } else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetamu_antinu")) { return (new ArgoNeuT_CC1Pi_XSec_1Dthetamu_antinu(samplekey)); } else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetapi_antinu")) { return (new ArgoNeuT_CC1Pi_XSec_1Dthetapi_antinu(samplekey)); } else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetamupi_antinu")) { return (new ArgoNeuT_CC1Pi_XSec_1Dthetamupi_antinu(samplekey)); /* BNL Samples */ } else #endif #ifndef __NO_BNL__ if (!name.compare("BNL_CCQE_XSec_1DEnu_nu")) { return (new BNL_CCQE_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BNL_CCQE_Evt_1DQ2_nu")) { return (new BNL_CCQE_Evt_1DQ2_nu(samplekey)); /* BNL CC1ppip samples */ } else if (!name.compare("BNL_CC1ppip_XSec_1DEnu_nu") || !name.compare("BNL_CC1ppip_XSec_1DEnu_nu_Uncorr") || !name.compare("BNL_CC1ppip_XSec_1DEnu_nu_W14Cut") || !name.compare("BNL_CC1ppip_XSec_1DEnu_nu_W14Cut_Uncorr")) { return (new BNL_CC1ppip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BNL_CC1ppip_Evt_1DQ2_nu") || !name.compare("BNL_CC1ppip_Evt_1DQ2_nu_W14Cut")) { return (new BNL_CC1ppip_Evt_1DQ2_nu(samplekey)); } else if (!name.compare("BNL_CC1ppip_Evt_1DcosthAdler_nu")) { return (new BNL_CC1ppip_Evt_1DcosthAdler_nu(samplekey)); } else if (!name.compare("BNL_CC1ppip_Evt_1Dphi_nu")) { return (new BNL_CC1ppip_Evt_1Dphi_nu(samplekey)); } else if (!name.compare("BNL_CC1ppip_Evt_1DWNpi_nu")) { return (new BNL_CC1ppip_Evt_1DWNpi_nu(samplekey)); } else if (!name.compare("BNL_CC1ppip_Evt_1DWNmu_nu")) { return (new BNL_CC1ppip_Evt_1DWNmu_nu(samplekey)); } else if (!name.compare("BNL_CC1ppip_Evt_1DWmupi_nu")) { return (new BNL_CC1ppip_Evt_1DWmupi_nu(samplekey)); /* BNL CC1npip samples */ } else if (!name.compare("BNL_CC1npip_XSec_1DEnu_nu") || !name.compare("BNL_CC1npip_XSec_1DEnu_nu_Uncorr")) { return (new BNL_CC1npip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BNL_CC1npip_Evt_1DQ2_nu")) { return (new BNL_CC1npip_Evt_1DQ2_nu(samplekey)); } else if (!name.compare("BNL_CC1npip_Evt_1DWNpi_nu")) { return (new BNL_CC1npip_Evt_1DWNpi_nu(samplekey)); } else if (!name.compare("BNL_CC1npip_Evt_1DWNmu_nu")) { return (new BNL_CC1npip_Evt_1DWNmu_nu(samplekey)); } else if (!name.compare("BNL_CC1npip_Evt_1DWmupi_nu")) { return (new BNL_CC1npip_Evt_1DWmupi_nu(samplekey)); /* BNL CC1pi0 samples */ } else if (!name.compare("BNL_CC1pi0_XSec_1DEnu_nu")) { return (new BNL_CC1pi0_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BNL_CC1pi0_Evt_1DQ2_nu")) { return (new BNL_CC1pi0_Evt_1DQ2_nu(samplekey)); } else if (!name.compare("BNL_CC1pi0_Evt_1DWNpi_nu")) { return (new BNL_CC1pi0_Evt_1DWNpi_nu(samplekey)); } else if (!name.compare("BNL_CC1pi0_Evt_1DWNmu_nu")) { return (new BNL_CC1pi0_Evt_1DWNmu_nu(samplekey)); } else if (!name.compare("BNL_CC1pi0_Evt_1DWmupi_nu")) { return (new BNL_CC1pi0_Evt_1DWmupi_nu(samplekey)); /* BNL multi-pi */ } else if (!name.compare("BNL_CC2pi_1pim1pip_XSec_1DEnu_nu")) { return (new BNL_CC2pi_1pim1pip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BNL_CC3pi_1pim2pip_XSec_1DEnu_nu")) { return (new BNL_CC3pi_1pim2pip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BNL_CC4pi_2pim2pip_XSec_1DEnu_nu")) { return (new BNL_CC4pi_2pim2pip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BNL_CC2pi_1pim1pip_Evt_1DWpippim_nu")) { return (new BNL_CC2pi_1pim1pip_Evt_1DWpippim_nu(samplekey)); } else if (!name.compare("BNL_CC2pi_1pim1pip_Evt_1DWpippr_nu")) { return (new BNL_CC2pi_1pim1pip_Evt_1DWpippr_nu(samplekey)); /* FNAL Samples */ } else #endif #ifndef __NO_FNAL__ if (!name.compare("FNAL_CCQE_Evt_1DQ2_nu")) { return (new FNAL_CCQE_Evt_1DQ2_nu(samplekey)); /* FNAL CC1ppip */ } else if (!name.compare("FNAL_CC1ppip_XSec_1DEnu_nu")) { return (new FNAL_CC1ppip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("FNAL_CC1ppip_XSec_1DQ2_nu")) { return (new FNAL_CC1ppip_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("FNAL_CC1ppip_Evt_1DQ2_nu")) { return (new FNAL_CC1ppip_Evt_1DQ2_nu(samplekey)); /* FNAL CC1ppim */ } else if (!name.compare("FNAL_CC1ppim_XSec_1DEnu_antinu")) { return (new FNAL_CC1ppim_XSec_1DEnu_antinu(samplekey)); /* BEBC Samples */ } else #endif #ifndef __NO_BEBC__ if (!name.compare("BEBC_CCQE_XSec_1DQ2_nu")) { return (new BEBC_CCQE_XSec_1DQ2_nu(samplekey)); /* BEBC CC1ppip samples */ } else if (!name.compare("BEBC_CC1ppip_XSec_1DEnu_nu")) { return (new BEBC_CC1ppip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BEBC_CC1ppip_XSec_1DQ2_nu")) { return (new BEBC_CC1ppip_XSec_1DQ2_nu(samplekey)); /* BEBC CC1npip samples */ } else if (!name.compare("BEBC_CC1npip_XSec_1DEnu_nu")) { return (new BEBC_CC1npip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BEBC_CC1npip_XSec_1DQ2_nu")) { return (new BEBC_CC1npip_XSec_1DQ2_nu(samplekey)); /* BEBC CC1pi0 samples */ } else if (!name.compare("BEBC_CC1pi0_XSec_1DEnu_nu")) { return (new BEBC_CC1pi0_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BEBC_CC1pi0_XSec_1DQ2_nu")) { return (new BEBC_CC1pi0_XSec_1DQ2_nu(samplekey)); /* BEBC CC1npim samples */ } else if (!name.compare("BEBC_CC1npim_XSec_1DEnu_antinu")) { return (new BEBC_CC1npim_XSec_1DEnu_antinu(samplekey)); } else if (!name.compare("BEBC_CC1npim_XSec_1DQ2_antinu")) { return (new BEBC_CC1npim_XSec_1DQ2_antinu(samplekey)); /* BEBC CC1ppim samples */ } else if (!name.compare("BEBC_CC1ppim_XSec_1DEnu_antinu")) { return (new BEBC_CC1ppim_XSec_1DEnu_antinu(samplekey)); } else if (!name.compare("BEBC_CC1ppim_XSec_1DQ2_antinu")) { return (new BEBC_CC1ppim_XSec_1DQ2_antinu(samplekey)); /* GGM CC1ppip samples */ } else #endif #ifndef __NO_GGM__ if (!name.compare("GGM_CC1ppip_XSec_1DEnu_nu")) { return (new GGM_CC1ppip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("GGM_CC1ppip_Evt_1DQ2_nu")) { return (new GGM_CC1ppip_Evt_1DQ2_nu(samplekey)); /* MiniBooNE Samples */ /* CCQE */ } else #endif #ifndef __NO_MiniBooNE__ if (!name.compare("MiniBooNE_CCQE_XSec_1DQ2_nu") || !name.compare("MiniBooNE_CCQELike_XSec_1DQ2_nu")) { return (new MiniBooNE_CCQE_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("MiniBooNE_CCQE_XSec_1DEnu_nu") || !name.compare("MiniBooNE_CCQELike_XSec_1DEnu_nu")) { return (new MiniBooNE_CCQE_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CCQE_XSec_1DQ2_antinu") || !name.compare("MiniBooNE_CCQELike_XSec_1DQ2_antinu") || !name.compare("MiniBooNE_CCQE_CTarg_XSec_1DQ2_antinu")) { return (new MiniBooNE_CCQE_XSec_1DQ2_antinu(samplekey)); } else if (!name.compare("MiniBooNE_CCQE_XSec_2DTcos_nu") || !name.compare("MiniBooNE_CCQELike_XSec_2DTcos_nu")) { return (new MiniBooNE_CCQE_XSec_2DTcos_nu(samplekey)); } else if (!name.compare("MiniBooNE_CCQE_XSec_2DTcos_antinu") || !name.compare("MiniBooNE_CCQELike_XSec_2DTcos_antinu")) { return (new MiniBooNE_CCQE_XSec_2DTcos_antinu(samplekey)); /* MiniBooNE CC1pi+ */ // 1D } else if (!name.compare("MiniBooNE_CC1pip_XSec_1DEnu_nu")) { return (new MiniBooNE_CC1pip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_1DQ2_nu")) { return (new MiniBooNE_CC1pip_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_1DTpi_nu")) { return (new MiniBooNE_CC1pip_XSec_1DTpi_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_1DTu_nu")) { return (new MiniBooNE_CC1pip_XSec_1DTu_nu(samplekey)); // 2D } else if (!name.compare("MiniBooNE_CC1pip_XSec_2DQ2Enu_nu")) { return (new MiniBooNE_CC1pip_XSec_2DQ2Enu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTpiCospi_nu")) { return (new MiniBooNE_CC1pip_XSec_2DTpiCospi_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTpiEnu_nu")) { return (new MiniBooNE_CC1pip_XSec_2DTpiEnu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTuCosmu_nu")) { return (new MiniBooNE_CC1pip_XSec_2DTuCosmu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTuEnu_nu")) { return (new MiniBooNE_CC1pip_XSec_2DTuEnu_nu(samplekey)); /* MiniBooNE CC1pi0 */ } else if (!name.compare("MiniBooNE_CC1pi0_XSec_1DEnu_nu")) { return (new MiniBooNE_CC1pi0_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pi0_XSec_1DQ2_nu")) { return (new MiniBooNE_CC1pi0_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pi0_XSec_1DTu_nu")) { return (new MiniBooNE_CC1pi0_XSec_1DTu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pi0_XSec_1Dcosmu_nu")) { return (new MiniBooNE_CC1pi0_XSec_1Dcosmu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pi0_XSec_1Dcospi0_nu")) { return (new MiniBooNE_CC1pi0_XSec_1Dcospi0_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pi0_XSec_1Dppi0_nu")) { return (new MiniBooNE_CC1pi0_XSec_1Dppi0_nu(samplekey)); } else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_antinu") || !name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_rhc")) { return (new MiniBooNE_NC1pi0_XSec_1Dcospi0_antinu(samplekey)); } else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_nu") || !name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_fhc")) { return (new MiniBooNE_NC1pi0_XSec_1Dcospi0_nu(samplekey)); } else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_antinu") || !name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_rhc")) { return (new MiniBooNE_NC1pi0_XSec_1Dppi0_antinu(samplekey)); } else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_nu") || !name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_fhc")) { return (new MiniBooNE_NC1pi0_XSec_1Dppi0_nu(samplekey)); /* MiniBooNE NCEL */ } else if (!name.compare("MiniBooNE_NCEL_XSec_Treco_nu")) { return (new MiniBooNE_NCEL_XSec_Treco_nu(samplekey)); } else #endif #ifndef __NO_MicroBooNE__ /* MicroBooNE Samples */ if (!name.compare("MicroBooNE_CCInc_XSec_2DPcos_nu")) { return (new MicroBooNE_CCInc_XSec_2DPcos_nu(samplekey)); } else if (!name.compare("MicroBooNE_CC1MuNp_XSec_1DPmu_nu") || !name.compare("MicroBooNE_CC1MuNp_XSec_1Dcosmu_nu") || !name.compare("MicroBooNE_CC1MuNp_XSec_1DPp_nu") || !name.compare("MicroBooNE_CC1MuNp_XSec_1Dcosp_nu") || !name.compare("MicroBooNE_CC1MuNp_XSec_1Dthetamup_nu")) { return (new MicroBooNE_CC1MuNp_XSec_1D_nu(samplekey)); } else #endif #ifndef __NO_MINERvA__ /* MINERvA Samples */ if (!name.compare("MINERvA_CCQE_XSec_1DQ2_nu") || !name.compare("MINERvA_CCQE_XSec_1DQ2_nu_20deg") || !name.compare("MINERvA_CCQE_XSec_1DQ2_nu_oldflux") || !name.compare("MINERvA_CCQE_XSec_1DQ2_nu_20deg_oldflux")) { return (new MINERvA_CCQE_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("MINERvA_CCQE_XSec_1DQ2_antinu") || !name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_20deg") || !name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_oldflux") || !name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_20deg_oldflux")) { return (new MINERvA_CCQE_XSec_1DQ2_antinu(samplekey)); } else if (!name.compare("MINERvA_CCQE_XSec_1DQ2_joint_oldflux") || !name.compare("MINERvA_CCQE_XSec_1DQ2_joint_20deg_oldflux") || !name.compare("MINERvA_CCQE_XSec_1DQ2_joint") || !name.compare("MINERvA_CCQE_XSec_1DQ2_joint_20deg")) { return (new MINERvA_CCQE_XSec_1DQ2_joint(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DEe_nue")) { return (new MINERvA_CC0pi_XSec_1DEe_nue(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_nue")) { return (new MINERvA_CC0pi_XSec_1DQ2_nue(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DThetae_nue")) { return (new MINERvA_CC0pi_XSec_1DThetae_nue(samplekey)); } else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dpmu_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Dthmu_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Dpprot_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Dthprot_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Dpnreco_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Ddalphat_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Ddpt_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Ddphit_nu")) { return (new MINERvA_CC0pinp_STV_XSec_1D_nu(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_nu_proton")) { return (new MINERvA_CC0pi_XSec_1DQ2_nu_proton(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtC_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtCH_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtFe_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtPb_nu")) { return (new MINERvA_CC0pi_XSec_1DQ2_Tgt_nu(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtRatioC_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtRatioFe_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtRatioPb_nu")) { return (new MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu(samplekey)); // Dan Ruterbories measurements of late 2018 } else if (!name.compare("MINERvA_CC0pi_XSec_2Dptpz_nu")) { return (new MINERvA_CC0pi_XSec_2D_nu(samplekey)); // } else if (!name.compare("MINERvA_CC0pi_XSec_3DptpzTp_nu")) { // return (new MINERvA_CC0pi_XSec_3DptpzTp_nu(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1Dpt_nu") || !name.compare("MINERvA_CC0pi_XSec_1Dpz_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2QE_nu") || !name.compare("MINERvA_CC0pi_XSec_1DEnuQE_nu")) { return (new MINERvA_CC0pi_XSec_1D_2018_nu(samplekey)); // C. Patrick's early 2018 measurements } else if (!name.compare("MINERvA_CC0pi_XSec_2Dptpz_antinu") || !name.compare("MINERvA_CC0pi_XSec_2DQ2QEEnuQE_antinu") || !name.compare("MINERvA_CC0pi_XSec_2DQ2QEEnuTrue_antinu")) { return (new MINERvA_CC0pi_XSec_2D_antinu(samplekey)); /* CC1pi+ */ // DONE } else if (!name.compare("MINERvA_CC1pip_XSec_1DTpi_nu") || !name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_20deg") || !name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_fluxcorr") || !name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_20deg_fluxcorr")) { return (new MINERvA_CC1pip_XSec_1DTpi_nu(samplekey)); // DONE } else if (!name.compare("MINERvA_CC1pip_XSec_1Dth_nu") || !name.compare("MINERvA_CC1pip_XSec_1Dth_nu_20deg") || !name.compare("MINERvA_CC1pip_XSec_1Dth_nu_fluxcorr") || !name.compare("MINERvA_CC1pip_XSec_1Dth_nu_20deg_fluxcorr")) { return (new MINERvA_CC1pip_XSec_1Dth_nu(samplekey)); } else if (!name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_2017") || !name.compare("MINERvA_CC1pip_XSec_1Dth_nu_2017") || !name.compare("MINERvA_CC1pip_XSec_1Dpmu_nu_2017") || !name.compare("MINERvA_CC1pip_XSec_1Dthmu_nu_2017") || !name.compare("MINERvA_CC1pip_XSec_1DQ2_nu_2017") || !name.compare("MINERvA_CC1pip_XSec_1DEnu_nu_2017")) { return (new MINERvA_CC1pip_XSec_1D_2017Update(samplekey)); /* CC1pi- */ } else if (!name.compare("MINERvA_CC1pim_XSec_1DEnu_antinu")) { return (new MINERvA_CC1pim_XSec_1DEnu_antinu(samplekey)); } else if (!name.compare("MINERvA_CC1pim_XSec_1DQ2_antinu")) { return (new MINERvA_CC1pim_XSec_1DQ2_antinu(samplekey)); } else if (!name.compare("MINERvA_CC1pim_XSec_1DTpi_antinu")) { return (new MINERvA_CC1pim_XSec_1DTpi_antinu(samplekey)); } else if (!name.compare("MINERvA_CC1pim_XSec_1Dpmu_antinu")) { return (new MINERvA_CC1pim_XSec_1Dpmu_antinu(samplekey)); } else if (!name.compare("MINERvA_CC1pim_XSec_1Dth_antinu")) { return (new MINERvA_CC1pim_XSec_1Dth_antinu(samplekey)); } else if (!name.compare("MINERvA_CC1pim_XSec_1Dthmu_antinu")) { return (new MINERvA_CC1pim_XSec_1Dthmu_antinu(samplekey)); /* CCNpi+ */ } else if (!name.compare("MINERvA_CCNpip_XSec_1Dth_nu") || !name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015") || !name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2016") || !name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015_20deg") || !name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015_fluxcorr") || !name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015_20deg_fluxcorr")) { return (new MINERvA_CCNpip_XSec_1Dth_nu(samplekey)); } else if (!name.compare("MINERvA_CCNpip_XSec_1DTpi_nu") || !name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2015") || !name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2016") || !name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2015_20deg") || !name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2015_fluxcorr") || !name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2015_20deg_fluxcorr")) { return (new MINERvA_CCNpip_XSec_1DTpi_nu(samplekey)); } else if (!name.compare("MINERvA_CCNpip_XSec_1Dthmu_nu")) { return (new MINERvA_CCNpip_XSec_1Dthmu_nu(samplekey)); } else if (!name.compare("MINERvA_CCNpip_XSec_1Dpmu_nu")) { return (new MINERvA_CCNpip_XSec_1Dpmu_nu(samplekey)); } else if (!name.compare("MINERvA_CCNpip_XSec_1DQ2_nu")) { return (new MINERvA_CCNpip_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("MINERvA_CCNpip_XSec_1DEnu_nu")) { return (new MINERvA_CCNpip_XSec_1DEnu_nu(samplekey)); /* MINERvA CC1pi0 anti-nu */ // Done } else if (!name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu") || !name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2015") || !name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2016") || !name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_fluxcorr") || !name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2015_fluxcorr") || !name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2016_fluxcorr")) { return (new MINERvA_CC1pi0_XSec_1Dth_antinu(samplekey)); } else if (!name.compare("MINERvA_CC1pi0_XSec_1Dppi0_antinu") || !name.compare("MINERvA_CC1pi0_XSec_1Dppi0_antinu_fluxcorr")) { return (new MINERvA_CC1pi0_XSec_1Dppi0_antinu(samplekey)); } else if (!name.compare("MINERvA_CC1pi0_XSec_1DTpi0_antinu")) { return (new MINERvA_CC1pi0_XSec_1DTpi0_antinu(samplekey)); // Done } else if (!name.compare("MINERvA_CC1pi0_XSec_1DQ2_antinu")) { return (new MINERvA_CC1pi0_XSec_1DQ2_antinu(samplekey)); // Done } else if (!name.compare("MINERvA_CC1pi0_XSec_1Dthmu_antinu")) { return (new MINERvA_CC1pi0_XSec_1Dthmu_antinu(samplekey)); // Done } else if (!name.compare("MINERvA_CC1pi0_XSec_1Dpmu_antinu")) { return (new MINERvA_CC1pi0_XSec_1Dpmu_antinu(samplekey)); // Done } else if (!name.compare("MINERvA_CC1pi0_XSec_1DEnu_antinu")) { return (new MINERvA_CC1pi0_XSec_1DEnu_antinu(samplekey)); // MINERvA CC1pi0 nu } else if (!name.compare("MINERvA_CC1pi0_XSec_1DTpi_nu") || !name.compare("MINERvA_CC1pi0_XSec_1Dth_nu") || !name.compare("MINERvA_CC1pi0_XSec_1Dpmu_nu") || !name.compare("MINERvA_CC1pi0_XSec_1Dthmu_nu") || !name.compare("MINERvA_CC1pi0_XSec_1DQ2_nu") || !name.compare("MINERvA_CC1pi0_XSec_1DEnu_nu") || !name.compare("MINERvA_CC1pi0_XSec_1DWexp_nu") || !name.compare("MINERvA_CC1pi0_XSec_1DPPi0Mass_nu") || !name.compare("MINERvA_CC1pi0_XSec_1DPPi0MassDelta_nu") || !name.compare("MINERvA_CC1pi0_XSec_1DCosAdler_nu") || !name.compare("MINERvA_CC1pi0_XSec_1DPhiAdler_nu")) { return (new MINERvA_CC1pi0_XSec_1D_nu(samplekey)); /* CCINC */ } else if (!name.compare("MINERvA_CCinc_XSec_2DEavq3_nu")) { return (new MINERvA_CCinc_XSec_2DEavq3_nu(samplekey)); } 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")) { return (new MINERvA_CCinc_XSec_1Dx_ratio(samplekey)); } 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")) { return (new MINERvA_CCinc_XSec_1DEnu_ratio(samplekey)); /* CCDIS */ } else if (!name.compare("MINERvA_CCDIS_XSec_1Dx_ratio_C12_CH") || !name.compare("MINERvA_CCDIS_XSec_1Dx_ratio_Fe56_CH") || !name.compare("MINERvA_CCDIS_XSec_1Dx_ratio_Pb208_CH")) { return (new MINERvA_CCDIS_XSec_1Dx_ratio(samplekey)); } else if (!name.compare("MINERvA_CCDIS_XSec_1DEnu_ratio_C12_CH") || !name.compare("MINERvA_CCDIS_XSec_1DEnu_ratio_Fe56_CH") || !name.compare("MINERvA_CCDIS_XSec_1DEnu_ratio_Pb208_CH")) { return (new MINERvA_CCDIS_XSec_1DEnu_ratio(samplekey)); /* CC-COH */ } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEnu_nu")) { return (new MINERvA_CCCOHPI_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEpi_nu")) { return (new MINERvA_CCCOHPI_XSec_1DEpi_nu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1Dth_nu")) { return (new MINERvA_CCCOHPI_XSec_1Dth_nu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DQ2_nu")) { return (new MINERvA_CCCOHPI_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEnu_antinu")) { return (new MINERvA_CCCOHPI_XSec_1DEnu_antinu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEpi_antinu")) { return (new MINERvA_CCCOHPI_XSec_1DEpi_antinu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1Dth_antinu")) { return (new MINERvA_CCCOHPI_XSec_1Dth_antinu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DQ2_antinu")) { return (new MINERvA_CCCOHPI_XSec_1DQ2_antinu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEnu_joint")) { return (new MINERvA_CCCOHPI_XSec_joint(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEpi_joint")) { return (new MINERvA_CCCOHPI_XSec_joint(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1Dth_joint")) { return (new MINERvA_CCCOHPI_XSec_joint(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DQ2_joint")) { return (new MINERvA_CCCOHPI_XSec_joint(samplekey)); /* T2K Samples */ } else #endif #ifndef __NO_T2K__ if (!name.compare("T2K_CC0pi_XSec_2DPcos_nu_I")) { return (new T2K_CC0pi_XSec_2DPcos_nu_I(samplekey)); } else if (!name.compare("T2K_CC0pi_XSec_2DPcos_nu_II")) { return (new T2K_CC0pi_XSec_2DPcos_nu_II(samplekey)); } else if (!name.compare("T2K_CCinc_XSec_2DPcos_nu_nonuniform")) { return (new T2K_CCinc_XSec_2DPcos_nu_nonuniform(samplekey)); } else if (!name.compare("T2K_CC0pi_XSec_H2O_2DPcos_anu")) { return (new T2K_CC0pi_XSec_H2O_2DPcos_anu(samplekey)); } else if (!name.compare("T2K_NuMu_CC0pi_O_XSec_2DPcos") || !name.compare("T2K_NuMu_CC0pi_C_XSec_2DPcos")) { return (new T2K_NuMu_CC0pi_OC_XSec_2DPcos(samplekey)); } else if (!name.compare("T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint")) { return (new T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint(samplekey)); } else if (!name.compare("T2K_NuMu_CC0pi_CH_XSec_2DPcos") || !name.compare("T2K_AntiNuMu_CC0pi_CH_XSec_2DPcos")) { return (new T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos(samplekey)); } else if (!name.compare("T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint")) { return (new T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint(samplekey)); } else if (!name.compare("T2K_nueCCinc_XSec_1Dpe_FHC") || !name.compare("T2K_nueCCinc_XSec_1Dpe_RHC") || !name.compare("T2K_nuebarCCinc_XSec_1Dpe_RHC")) { return (new T2K_nueCCinc_XSec_1Dpe(samplekey)); } else if (!name.compare("T2K_nueCCinc_XSec_1Dthe_FHC") || !name.compare("T2K_nueCCinc_XSec_1Dthe_RHC") || !name.compare("T2K_nuebarCCinc_XSec_1Dthe_RHC")) { return (new T2K_nueCCinc_XSec_1Dthe(samplekey)); } else if (!name.compare("T2K_nueCCinc_XSec_1Dpe_joint")) { return (new T2K_nueCCinc_XSec_1Dpe_joint(samplekey)); } else if (!name.compare("T2K_nueCCinc_XSec_1Dthe_joint")) { return (new T2K_nueCCinc_XSec_1Dthe_joint(samplekey)); } else if (!name.compare("T2K_nueCCinc_XSec_joint")) { return (new T2K_nueCCinc_XSec_joint(samplekey)); /* T2K CC1pi+ CH samples */ // Comment these out for now because we don't have the proper data } else if (!name.compare("T2K_CC1pip_CH_XSec_2Dpmucosmu_nu")) { return (new T2K_CC1pip_CH_XSec_2Dpmucosmu_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_CH_XSec_1Dppi_nu")) { return (new T2K_CC1pip_CH_XSec_1Dppi_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthpi_nu")) { return (new T2K_CC1pip_CH_XSec_1Dthpi_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthmupi_nu")) { return (new T2K_CC1pip_CH_XSec_1Dthmupi_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_CH_XSec_1DQ2_nu")) { return (new T2K_CC1pip_CH_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_CH_XSec_1DAdlerPhi_nu")) { return (new T2K_CC1pip_CH_XSec_1DAdlerPhi_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_CH_XSec_1DCosThAdler_nu")) { return (new T2K_CC1pip_CH_XSec_1DCosThAdler_nu(samplekey)); + } else if (!name.compare("T2K_CCCOH_C12_XSec_1DEnu_nu")) { + return (new T2K_CCCOH_C12_XSec_1DEnu_nu(samplekey)); + /* T2K CC1pi+ H2O samples */ } else if (!name.compare("T2K_CC1pip_H2O_XSec_1DEnuDelta_nu")) { return (new T2K_CC1pip_H2O_XSec_1DEnuDelta_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_H2O_XSec_1DEnuMB_nu")) { return (new T2K_CC1pip_H2O_XSec_1DEnuMB_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dcosmu_nu")) { return (new T2K_CC1pip_H2O_XSec_1Dcosmu_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dcosmupi_nu")) { return (new T2K_CC1pip_H2O_XSec_1Dcosmupi_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dcospi_nu")) { return (new T2K_CC1pip_H2O_XSec_1Dcospi_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dpmu_nu")) { return (new T2K_CC1pip_H2O_XSec_1Dpmu_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dppi_nu")) { return (new T2K_CC1pip_H2O_XSec_1Dppi_nu(samplekey)); /* T2K CC0pi + np CH samples */ } else if (!name.compare("T2K_CC0pinp_STV_XSec_1Ddpt_nu")) { return (new T2K_CC0pinp_STV_XSec_1Ddpt_nu(samplekey)); } else if (!name.compare("T2K_CC0pinp_STV_XSec_1Ddphit_nu")) { return (new T2K_CC0pinp_STV_XSec_1Ddphit_nu(samplekey)); } else if (!name.compare("T2K_CC0pinp_STV_XSec_1Ddat_nu")) { return (new T2K_CC0pinp_STV_XSec_1Ddat_nu(samplekey)); } else if (!name.compare("T2K_CC0piWithProtons_XSec_2018_multidif_0p_1p_Np") || !name.compare("T2K_CC0piWithProtons_XSec_2018_multidif_0p_1p") || !name.compare("T2K_CC0piWithProtons_XSec_2018_multidif_0p") || !name.compare("T2K_CC0piWithProtons_XSec_2018_multidif_1p")) { return (new T2K_CC0piWithProtons_XSec_2018_multidif_0p_1p_Np(samplekey)); } else if (!name.compare("T2K_CC0pinp_ifk_XSec_3Dinfp_nu")) { return (new T2K_CC0pinp_ifk_XSec_3Dinfp_nu(samplekey)); } else if (!name.compare("T2K_CC0pinp_ifk_XSec_3Dinfa_nu")) { return (new T2K_CC0pinp_ifk_XSec_3Dinfa_nu(samplekey)); } else if (!name.compare("T2K_CC0pinp_ifk_XSec_3Dinfip_nu")) { return (new T2K_CC0pinp_ifk_XSec_3Dinfip_nu(samplekey)); // SciBooNE COH studies } else #endif #ifndef __NO_SciBooNE__ if (!name.compare("SciBooNE_CCCOH_STOP_NTrks_nu")) { return (new SciBooNE_CCCOH_STOP_NTrks_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_1TRK_1DQ2_nu")) { return (new SciBooNE_CCCOH_1TRK_1DQ2_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_1TRK_1Dpmu_nu")) { return (new SciBooNE_CCCOH_1TRK_1Dpmu_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_1TRK_1Dthetamu_nu")) { return (new SciBooNE_CCCOH_1TRK_1Dthetamu_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPr_1DQ2_nu")) { return (new SciBooNE_CCCOH_MuPr_1DQ2_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPr_1Dpmu_nu")) { return (new SciBooNE_CCCOH_MuPr_1Dpmu_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPr_1Dthetamu_nu")) { return (new SciBooNE_CCCOH_MuPr_1Dthetamu_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiVA_1DQ2_nu")) { return (new SciBooNE_CCCOH_MuPiVA_1DQ2_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiVA_1Dpmu_nu")) { return (new SciBooNE_CCCOH_MuPiVA_1Dpmu_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiVA_1Dthetamu_nu")) { return (new SciBooNE_CCCOH_MuPiVA_1Dthetamu_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1DQ2_nu")) { return (new SciBooNE_CCCOH_MuPiNoVA_1DQ2_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1Dthetapr_nu")) { return (new SciBooNE_CCCOH_MuPiNoVA_1Dthetapr_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1Dthetapi_nu")) { return (new SciBooNE_CCCOH_MuPiNoVA_1Dthetapi_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1Dthetamu_nu")) { return (new SciBooNE_CCCOH_MuPiNoVA_1Dthetamu_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1Dpmu_nu")) { return (new SciBooNE_CCCOH_MuPiNoVA_1Dpmu_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_STOPFINAL_1DQ2_nu")) { return (new SciBooNE_CCCOH_STOPFINAL_1DQ2_nu(samplekey)); } else if (!name.compare("SciBooNE_CCInc_XSec_1DEnu_nu") || !name.compare("SciBooNE_CCInc_XSec_1DEnu_nu_NEUT") || !name.compare("SciBooNE_CCInc_XSec_1DEnu_nu_NUANCE")) { return (new SciBooNE_CCInc_XSec_1DEnu_nu(samplekey)); /* K2K Samples */ /* NC1pi0 */ } else #endif #ifndef __NO_K2K__ if (!name.compare("K2K_NC1pi0_Evt_1Dppi0_nu")) { return (new K2K_NC1pi0_Evt_1Dppi0_nu(samplekey)); /* Fake Studies */ } else #endif if (name.find("ExpMultDist_CCQE_XSec_1D") != std::string::npos && name.find("_FakeStudy") != std::string::npos) { return ( 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) { return ( new ExpMultDist_CCQE_XSec_2DVar_FakeStudy(name, file, rw, type, fkdt)); } else if (name.find("GenericFlux") != std::string::npos) { return (new GenericFlux_Tester(name, file, rw, type, fkdt)); } else if (name.find("GenericVectors") != std::string::npos) { return (new GenericFlux_Vectors(name, file, rw, type, fkdt)); } else if (!name.compare("T2K2017_FakeData")) { return (new T2K2017_FakeData(samplekey)); } else if (!name.compare("MCStudy_CCQE")) { return (new MCStudy_CCQEHistograms(name, file, rw, type, fkdt)); } else if (!name.compare("ElectronFlux_FlatTree")) { return (new ElectronFlux_FlatTree(name, file, rw, type, fkdt)); } else if (name.find("ElectronData_") != std::string::npos) { return new ElectronScattering_DurhamData(samplekey); } else if (name.find("MuonValidation_") != std::string::npos) { return (new MCStudy_MuonValidation(name, file, rw, type, fkdt)); } else if (!name.compare("NIWGOfficialPlots")) { return (new OfficialNIWGPlots(samplekey)); } else if ((name.find("SigmaEnuHists") != std::string::npos) || (name.find("SigmaEnuPerEHists") != std::string::npos)) { return (new SigmaEnuHists(samplekey)); } else if (!name.compare("Simple_Osc")) { return (new Simple_Osc(samplekey)); } else if (!name.compare("Smear_SVDUnfold_Propagation_Osc")) { return (new Smear_SVDUnfold_Propagation_Osc(samplekey)); } else { NUIS_ABORT("Error: No such sample: " << name << std::endl); } // Return NULL if no sample loaded. return NULL; } // namespace SampleUtils } // namespace SampleUtils diff --git a/src/FitBase/Measurement1D.cxx b/src/FitBase/Measurement1D.cxx index 0adce17..256dd07 100644 --- a/src/FitBase/Measurement1D.cxx +++ b/src/FitBase/Measurement1D.cxx @@ -1,2019 +1,2063 @@ // Copyright 2016-2021 L. Pickering, P. Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This ile is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #include "Measurement1D.h" //******************************************************************** Measurement1D::Measurement1D(void) { //******************************************************************** // XSec Scalings fScaleFactor = -1.0; fCurrentNorm = 1.0; // Histograms fDataHist = NULL; fDataTrue = NULL; fMCHist = NULL; fMCFine = NULL; fMCWeighted = NULL; fMaskHist = NULL; // Covar covar = NULL; fFullCovar = NULL; fShapeCovar = NULL; fCovar = NULL; fInvert = NULL; fDecomp = NULL; fResidualHist = NULL; fChi2LessBinHist = NULL; // Fake Data fFakeDataInput = ""; fFakeDataFile = NULL; // Options fDefaultTypes = "FIX/FULL/CHI2"; fAllowedTypes = - "FIX,FREE,SHAPE/FULL,DIAG/CHI2/NORM/ENUCORR/Q2CORR/ENU1D/MASK/NOWIDTH"; + "FIX,FREE,SHAPE/FULL,DIAG/CHI2/NORM/ENUCORR/Q2CORR/ENU1D/MASK/NOWIDTH/SINGLEBIN"; fIsFix = false; fIsShape = false; fIsFree = false; fIsDiag = false; fIsFull = false; + fIsSingleBin = false; fAddNormPen = false; fIsMask = false; fIsChi2SVD = false; fIsRawEvents = false; fIsNoWidth = false; fIsDifXSec = false; fIsEnu1D = false; fIsWriting = false; fSaveFine = true; // Inputs fInput = NULL; fRW = NULL; // Extra Histograms fMCHist_Modes = NULL; + fMCFine_Modes = NULL; } //******************************************************************** Measurement1D::~Measurement1D(void) { //******************************************************************** if (fDataHist) delete fDataHist; if (fDataTrue) delete fDataTrue; if (fMCHist) delete fMCHist; if (fMCFine) delete fMCFine; if (fMCWeighted) delete fMCWeighted; if (fMaskHist) delete fMaskHist; if (covar) delete covar; if (fFullCovar) delete fFullCovar; if (fShapeCovar) delete fShapeCovar; if (fCovar) delete fCovar; if (fInvert) delete fInvert; if (fDecomp) delete fDecomp; delete fResidualHist; delete fChi2LessBinHist; } //******************************************************************** void Measurement1D::FinaliseSampleSettings() { //******************************************************************** MeasurementBase::FinaliseSampleSettings(); // Setup naming + renaming fName = fSettings.GetName(); fSettings.SetS("originalname", fName); if (fSettings.Has("rename")) { fName = fSettings.GetS("rename"); fSettings.SetS("name", fName); } // Setup all other options NUIS_LOG(SAM, "Finalising Sample Settings: " << fName); if ((fSettings.GetS("originalname").find("Evt") != std::string::npos)) { fIsRawEvents = true; NUIS_LOG(SAM, "Found event rate measurement but using poisson likelihoods."); } if (fSettings.GetS("originalname").find("XSec_1DEnu") != std::string::npos) { fIsEnu1D = true; NUIS_LOG(SAM, "Found XSec Enu measurement, applying flux integrated scaling, " << "not flux averaged!"); } if (fIsEnu1D && fIsRawEvents) { NUIS_ERR(FTL, "Found 1D Enu XSec distribution AND fIsRawEvents, is this " "really correct?!"); NUIS_ERR(FTL, "Check experiment constructor for " << fName << " and correct this!"); NUIS_ERR(FTL, "I live in " << __FILE__ << ":" << __LINE__); throw; } if (!fRW) fRW = FitBase::GetRW(); if (!fInput and !fIsJoint) SetupInputs(fSettings.GetS("input")); // Setup options SetFitOptions(fDefaultTypes); // defaults SetFitOptions(fSettings.GetS("type")); // user specified EnuMin = GeneralUtils::StrToDbl(fSettings.GetS("enu_min")); EnuMax = GeneralUtils::StrToDbl(fSettings.GetS("enu_max")); } //******************************************************************** void Measurement1D::CreateDataHistogram(int dimx, double *binx) { //******************************************************************** if (fDataHist) delete fDataHist; fDataHist = new TH1D((fSettings.GetName() + "_data").c_str(), (fSettings.GetFullTitles()).c_str(), dimx, binx); } //******************************************************************** void Measurement1D::SetDataFromTextFile(std::string datafile) { //******************************************************************** NUIS_LOG(SAM, "Reading data from text file: " << datafile); fDataHist = PlotUtils::GetTH1DFromFile( datafile, fSettings.GetName() + "_data", fSettings.GetFullTitles()); } //******************************************************************** void Measurement1D::SetDataFromRootFile(std::string datafile, std::string histname) { //******************************************************************** NUIS_LOG(SAM, "Reading data from root file: " << datafile << ";" << histname); fDataHist = PlotUtils::GetTH1DFromRootFile(datafile, histname); fDataHist->SetNameTitle((fSettings.GetName() + "_data").c_str(), (fSettings.GetFullTitles()).c_str()); return; }; //******************************************************************** void Measurement1D::SetEmptyData() { //******************************************************************** fDataHist = new TH1D("EMPTY_DATA", "EMPTY_DATA", 1, 0.0, 1.0); } //******************************************************************** void Measurement1D::SetPoissonErrors() { //******************************************************************** if (!fDataHist) { NUIS_ERR(FTL, "Need a data hist to setup possion errors! "); NUIS_ERR(FTL, "Setup Data First!"); throw; } for (int i = 0; i < fDataHist->GetNbinsX() + 1; i++) { fDataHist->SetBinError(i + 1, sqrt(fDataHist->GetBinContent(i + 1))); } } //******************************************************************** void Measurement1D::SetCovarFromDiagonal(TH1D *data) { //******************************************************************** if (!data and fDataHist) { data = fDataHist; } if (data) { NUIS_LOG(SAM, "Setting diagonal covariance for: " << data->GetName()); fFullCovar = StatUtils::MakeDiagonalCovarMatrix(data); covar = StatUtils::GetInvert(fFullCovar, true); fDecomp = StatUtils::GetDecomp(fFullCovar); } else { NUIS_ABORT("No data input provided to set diagonal covar from!"); } } //******************************************************************** void Measurement1D::SetCovarFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) { dim = fDataHist->GetNbinsX(); } NUIS_LOG(SAM, "Reading covariance from text file: " << covfile); fFullCovar = StatUtils::GetCovarFromTextFile(covfile, dim); covar = StatUtils::GetInvert(fFullCovar, true); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement1D::SetCovarFromMultipleTextFiles(std::string covfiles, int dim) { //******************************************************************** if (dim == -1) { dim = fDataHist->GetNbinsX(); } std::vector covList = GeneralUtils::ParseToStr(covfiles, ";"); fFullCovar = new TMatrixDSym(dim); for (uint i = 0; i < covList.size(); ++i) { NUIS_LOG(SAM, "Reading covariance from text file: " << covList[i]); TMatrixDSym *temp_cov = StatUtils::GetCovarFromTextFile(covList[i], dim); (*fFullCovar) += (*temp_cov); delete temp_cov; } covar = StatUtils::GetInvert(fFullCovar, true); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement1D::SetCovarFromRootFile(std::string covfile, std::string histname) { //******************************************************************** NUIS_LOG(SAM, "Reading covariance from root file: " << covfile << ";" << histname); fFullCovar = StatUtils::GetCovarFromRootFile(covfile, histname); covar = StatUtils::GetInvert(fFullCovar, true); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement1D::SetCovarInvertFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) { dim = fDataHist->GetNbinsX(); } NUIS_LOG(SAM, "Reading inverted covariance from text file: " << covfile); covar = StatUtils::GetCovarFromTextFile(covfile, dim); fFullCovar = StatUtils::GetInvert(covar, true); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement1D::SetCovarInvertFromRootFile(std::string covfile, std::string histname) { //******************************************************************** NUIS_LOG(SAM, "Reading inverted covariance from text file: " << covfile << ";" << histname); covar = StatUtils::GetCovarFromRootFile(covfile, histname); fFullCovar = StatUtils::GetInvert(covar, true); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement1D::SetCorrelationFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) dim = fDataHist->GetNbinsX(); NUIS_LOG(SAM, "Reading data correlations from text file: " << covfile << ";" << dim); TMatrixDSym *correlation = StatUtils::GetCovarFromTextFile(covfile, dim); if (!fDataHist) { NUIS_ABORT("Trying to set correlations from text file but there is no " "data to build it from. \n" << "In constructor make sure data is set before " "SetCorrelationFromTextFile is called. \n"); } // Fill covar from data errors and correlations fFullCovar = new TMatrixDSym(dim); for (int i = 0; i < fDataHist->GetNbinsX(); i++) { for (int j = 0; j < fDataHist->GetNbinsX(); j++) { (*fFullCovar)(i, j) = (*correlation)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1.E76; } } // Fill other covars. covar = StatUtils::GetInvert(fFullCovar, true); fDecomp = StatUtils::GetDecomp(fFullCovar); delete correlation; } //******************************************************************** void Measurement1D::SetCorrelationFromMultipleTextFiles(std::string corrfiles, int dim) { //******************************************************************** if (dim == -1) { dim = fDataHist->GetNbinsX(); } std::vector corrList = GeneralUtils::ParseToStr(corrfiles, ";"); fFullCovar = new TMatrixDSym(dim); for (uint i = 0; i < corrList.size(); ++i) { NUIS_LOG(SAM, "Reading covariance from text file: " << corrList[i]); TMatrixDSym *temp_cov = StatUtils::GetCovarFromTextFile(corrList[i], dim); for (int i = 0; i < fDataHist->GetNbinsX(); i++) { for (int j = 0; j < fDataHist->GetNbinsX(); j++) { (*temp_cov)(i, j) = (*temp_cov)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1.E76; } } (*fFullCovar) += (*temp_cov); delete temp_cov; } covar = StatUtils::GetInvert(fFullCovar, true); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement1D::SetCorrelationFromRootFile(std::string covfile, std::string histname) { //******************************************************************** NUIS_LOG(SAM, "Reading data correlations from text file: " << covfile << ";" << histname); TMatrixDSym *correlation = StatUtils::GetCovarFromRootFile(covfile, histname); if (!fDataHist) { NUIS_ABORT("Trying to set correlations from text file but there is no " "data to build it from. \n" << "In constructor make sure data is set before " "SetCorrelationFromTextFile is called. \n"); } // Fill covar from data errors and correlations fFullCovar = new TMatrixDSym(fDataHist->GetNbinsX()); for (int i = 0; i < fDataHist->GetNbinsX(); i++) { for (int j = 0; j < fDataHist->GetNbinsX(); j++) { (*fFullCovar)(i, j) = (*correlation)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1.E76; } } // Fill other covars. covar = StatUtils::GetInvert(fFullCovar, true); fDecomp = StatUtils::GetDecomp(fFullCovar); delete correlation; } //******************************************************************** void Measurement1D::SetCholDecompFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) { dim = fDataHist->GetNbinsX(); } NUIS_LOG(SAM, "Reading cholesky from text file: " << covfile); TMatrixD *temp = StatUtils::GetMatrixFromTextFile(covfile, dim, dim); TMatrixD *trans = (TMatrixD *)temp->Clone(); trans->T(); (*trans) *= (*temp); fFullCovar = new TMatrixDSym(dim, trans->GetMatrixArray(), ""); covar = StatUtils::GetInvert(fFullCovar, true); fDecomp = StatUtils::GetDecomp(fFullCovar); delete temp; delete trans; } //******************************************************************** void Measurement1D::SetCholDecompFromRootFile(std::string covfile, std::string histname) { //******************************************************************** NUIS_LOG(SAM, "Reading cholesky decomp from root file: " << covfile << ";" << histname); TMatrixD *temp = StatUtils::GetMatrixFromRootFile(covfile, histname); TMatrixD *trans = (TMatrixD *)temp->Clone(); trans->T(); (*trans) *= (*temp); fFullCovar = new TMatrixDSym(temp->GetNrows(), trans->GetMatrixArray(), ""); covar = StatUtils::GetInvert(fFullCovar, true); fDecomp = StatUtils::GetDecomp(fFullCovar); delete temp; delete trans; } void Measurement1D::SetShapeCovar() { // Return if this is missing any pre-requisites if (!fFullCovar) return; if (!fDataHist) return; // Also return if it's bloody stupid under the circumstances if (fIsDiag) return; fShapeCovar = StatUtils::ExtractShapeOnlyCovar(fFullCovar, fDataHist); return; } //******************************************************************** void Measurement1D::ScaleData(double scale) { //******************************************************************** fDataHist->Scale(scale); } //******************************************************************** void Measurement1D::ScaleDataErrors(double scale) { //******************************************************************** for (int i = 0; i < fDataHist->GetNbinsX(); i++) { fDataHist->SetBinError(i + 1, fDataHist->GetBinError(i + 1) * scale); } } //******************************************************************** void Measurement1D::ScaleCovar(double scale) { //******************************************************************** (*fFullCovar) *= scale; (*covar) *= 1.0 / scale; (*fDecomp) *= sqrt(scale); } //******************************************************************** void Measurement1D::SetBinMask(std::string maskfile) { //******************************************************************** if (!fIsMask) return; NUIS_LOG(SAM, "Reading bin mask from file: " << maskfile); // Create a mask histogram with dim of data int nbins = fDataHist->GetNbinsX(); fMaskHist = new TH1I((fSettings.GetName() + "_BINMASK").c_str(), (fSettings.GetName() + "_BINMASK; Bin; Mask?").c_str(), nbins, 0, nbins); std::string line; std::ifstream mask(maskfile.c_str(), std::ifstream::in); if (!mask.is_open()) { NUIS_ABORT("Cannot find mask file."); } while (std::getline(mask >> std::ws, line, '\n')) { std::vector entries = GeneralUtils::ParseToInt(line, " "); // Skip lines with poorly formatted lines if (entries.size() < 2) { NUIS_LOG(WRN, "Measurement1D::SetBinMask(), couldn't parse line: " << line); continue; } // The first index should be the bin number, the second should be the mask // value. int val = 0; if (entries[1] > 0) val = 1; fMaskHist->SetBinContent(entries[0], val); } // Apply masking by setting masked data bins to zero PlotUtils::MaskBins(fDataHist, fMaskHist); return; } //******************************************************************** void Measurement1D::FinaliseMeasurement() { //******************************************************************** NUIS_LOG(SAM, "Finalising Measurement: " << fName); if (fSettings.GetB("onlymc")) { if (fDataHist) delete fDataHist; fDataHist = new TH1D("empty_data", "empty_data", 1, 0.0, 1.0); } // Make sure data is setup if (!fDataHist) { NUIS_ABORT("No data has been setup inside " << fName << " constructor!"); } // Make sure covariances are setup if (!fFullCovar) { fIsDiag = true; SetCovarFromDiagonal(fDataHist); } else if (fIsDiag) { // Have covariance but also set Diag NUIS_LOG(SAM, "Have full covariance for sample " << GetName() << " but only using diagonal elements for likelihood"); int nbins = fFullCovar->GetNcols(); for (int i = 0; i < nbins; ++i) { for (int j = 0; j < nbins; ++j) { if (i != j) { (*fFullCovar)[i][j] = 0; } } } delete covar; covar = NULL; delete fDecomp; fDecomp = NULL; } if (!covar) { covar = StatUtils::GetInvert(fFullCovar, true); } if (!fDecomp) { fDecomp = StatUtils::GetDecomp(fFullCovar); } // Push the diagonals of fFullCovar onto the data histogram // Comment this out until the covariance/data scaling is consistent! StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar, 1E-38); // If shape only, set covar and fDecomp using the shape-only matrix (if set) if (fIsShape && fShapeCovar && FitPar::Config().GetParB("UseShapeCovar")) { if (covar) delete covar; covar = StatUtils::GetInvert(fShapeCovar, true); if (fDecomp) delete fDecomp; fDecomp = StatUtils::GetDecomp(fFullCovar); fUseShapeNormDecomp = FitPar::Config().GetParB("UseShapeNormDecomp"); if (fUseShapeNormDecomp) { fNormError = 0; // From https://arxiv.org/pdf/2003.00088.pdf for (int i = 0; i < fFullCovar->GetNcols(); ++i) { for (int j = 0; j < fFullCovar->GetNcols(); ++j) { fNormError += (*fFullCovar)[i][j]; } } NUIS_LOG(SAM, "Sample: " << fName << ", using shape/norm decomp with norm error: " << fNormError); } } // Setup fMCHist from data fMCHist = (TH1D *)fDataHist->Clone(); fMCHist->SetNameTitle((fSettings.GetName() + "_MC").c_str(), (fSettings.GetFullTitles()).c_str()); fMCHist->Reset(); - // Setup fMCFine - fMCFine = new TH1D("mcfine", "mcfine", fDataHist->GetNbinsX() * 8, - fMCHist->GetBinLowEdge(1), - fMCHist->GetBinLowEdge(fDataHist->GetNbinsX() + 1)); + // Setup fMCFine, provide option to create it before in the sample constructor first + if (!fMCFine){ + fMCFine = new TH1D("mcfine", "mcfine", fDataHist->GetNbinsX() * 8, + fMCHist->GetBinLowEdge(1), + fMCHist->GetBinLowEdge(fDataHist->GetNbinsX() + 1)); + } fMCFine->SetNameTitle((fSettings.GetName() + "_MC_FINE").c_str(), - (fSettings.GetFullTitles()).c_str()); + (fSettings.GetFullTitles()).c_str()); fMCFine->Reset(); // Setup MC Stat fMCStat = (TH1D *)fMCHist->Clone(); fMCStat->Reset(); // Search drawopts for possible types to include by default std::string drawopts = FitPar::Config().GetParS("drawopts"); if (drawopts.find("MODES") != std::string::npos) { fMCHist_Modes = new TrueModeStack((fSettings.GetName() + "_MODES").c_str(), ("True Channels"), fMCHist); fMCHist_Modes ->SetTitleX(fDataHist->GetXaxis()->GetTitle()); fMCHist_Modes ->SetTitleY(fDataHist->GetYaxis()->GetTitle()); SetAutoProcessTH1(fMCHist_Modes, kCMD_Reset, kCMD_Norm, kCMD_Write); + + fMCFine_Modes = new TrueModeStack((fSettings.GetName() + "_MODES").c_str(), + ("True Channels"), fMCFine); + fMCFine_Modes ->SetTitleX(fDataHist->GetXaxis()->GetTitle()); + fMCFine_Modes ->SetTitleY(fDataHist->GetYaxis()->GetTitle()); + + SetAutoProcessTH1(fMCFine_Modes, kCMD_Reset, kCMD_Norm, kCMD_Write); } if (fSettings.Has("maskfile") && fSettings.Has("maskhist")) { fMaskHist = dynamic_cast(PlotUtils::GetTH1FromRootFile( fSettings.GetS("maskfile"), fSettings.GetS("maskhist"))); fIsMask = bool(fMaskHist); NUIS_LOG(SAM, "Loaded mask histogram: " << fSettings.GetS("maskhist") << " from " << fSettings.GetS("maskfile")); } else if (fIsMask) { // Setup bin masks using sample name std::string curname = fName; std::string origname = fSettings.GetS("originalname"); // Check rename.mask std::string maskloc = FitPar::Config().GetParDIR(curname + ".mask"); // Check origname.mask if (maskloc.empty()) maskloc = FitPar::Config().GetParDIR(origname + ".mask"); // Check database if (maskloc.empty()) { maskloc = FitPar::GetDataBase() + "/masks/" + origname + ".mask"; } // Setup Bin Mask SetBinMask(maskloc); } if (fScaleFactor < 0) { NUIS_ERR(FTL, "I found a negative fScaleFactor in " << __FILE__ << ":" << __LINE__); NUIS_ERR(FTL, "fScaleFactor = " << fScaleFactor); NUIS_ERR(FTL, "EXITING"); throw; } if (fAddNormPen) { if (!fUseShapeNormDecomp) { fNormError = fSettings.GetNormError(); } if (fNormError <= 0.0) { NUIS_ERR(FTL, "Norm error for class " << fName << " is 0.0!"); NUIS_ERR(FTL, "If you want to use it please add fNormError=VAL"); throw; } } // Create and fill Weighted Histogram if (!fMCWeighted) { fMCWeighted = (TH1D *)fMCHist->Clone(); fMCWeighted->SetNameTitle((fName + "_MCWGHTS").c_str(), (fName + "_MCWGHTS" + fPlotTitles).c_str()); fMCWeighted->GetYaxis()->SetTitle("Weighted Events"); } } //******************************************************************** void Measurement1D::SetFitOptions(std::string opt) { //******************************************************************** // Do nothing if default given if (opt == "DEFAULT") return; // CHECK Conflicting Fit Options std::vector fit_option_allow = GeneralUtils::ParseToStr(fAllowedTypes, "/"); for (UInt_t i = 0; i < fit_option_allow.size(); i++) { std::vector fit_option_section = GeneralUtils::ParseToStr(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) { NUIS_ABORT( "ERROR: Conflicting fit options provided: " << opt << std::endl << "Conflicting group = " << fit_option_section.at(i) << std::endl << "You should only supply one of these options in card file."); } } } // Check all options are allowed std::vector fit_options_input = GeneralUtils::ParseToStr(opt, "/"); for (UInt_t i = 0; i < fit_options_input.size(); i++) { if (fAllowedTypes.find(fit_options_input.at(i)) == std::string::npos) { NUIS_ERR(WRN, "ERROR: Fit Option '" << fit_options_input.at(i) << "' Provided is not allowed for this measurement."); NUIS_ERR(WRN, "Fit Options should be provided as a '/' separated list " "(e.g. FREE/DIAG/NORM)"); NUIS_ABORT("Available options for " << fName << " are '" << fAllowedTypes << "'"); } } // Set TYPE fFitType = opt; // FIX,SHAPE,FREE if (opt.find("FIX") != std::string::npos) { fIsFree = fIsShape = false; fIsFix = true; } else if (opt.find("SHAPE") != std::string::npos) { fIsFree = fIsFix = false; fIsShape = true; } else if (opt.find("FREE") != std::string::npos) { fIsFix = fIsShape = false; fIsFree = true; } // DIAG,FULL (or default to full) if (opt.find("DIAG") != std::string::npos) { fIsDiag = true; fIsFull = false; } else if (opt.find("FULL") != std::string::npos) { fIsDiag = false; fIsFull = true; } // CHI2/LL (OTHERS?) if (opt.find("LOG") != std::string::npos) { fIsChi2 = false; NUIS_ERR(FTL, "No other LIKELIHOODS properly supported!"); NUIS_ERR(FTL, "Try to use a chi2!"); throw; } else { fIsChi2 = true; } // EXTRAS if (opt.find("RAW") != std::string::npos) fIsRawEvents = true; if (opt.find("NOWIDTH") != std::string::npos) fIsNoWidth = true; if (opt.find("DIF") != std::string::npos) fIsDifXSec = true; if (opt.find("ENU1D") != std::string::npos) fIsEnu1D = true; if (opt.find("NORM") != std::string::npos) fAddNormPen = true; if (opt.find("MASK") != std::string::npos) fIsMask = true; + if (opt.find("SINGLEBIN") != std::string::npos){ + fIsSingleBin = true; + fIsNoWidth = false; + fIsEnu1D = false; + } return; }; //******************************************************************** void Measurement1D::SetSmearingMatrix(std::string smearfile, int truedim, int recodim) { //******************************************************************** // 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(), std::ifstream::in); // Note that the smearing matrix may be rectangular. fSmearMatrix = new TMatrixD(truedim, recodim); if (smear.is_open()) { NUIS_LOG(SAM, "Reading smearing matrix from file: " << smearfile); } else { NUIS_ABORT("Smearing matrix provided is incorrect: " << smearfile); } while (std::getline(smear >> std::ws, line, '\n')) { int column = 0; std::vector entries = GeneralUtils::ParseToDbl(line, " "); for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { (*fSmearMatrix)(row, column) = (*iter) / 100.; // Convert to fraction from // percentage (this may not be // general enough) column++; } row++; } return; } //******************************************************************** void Measurement1D::ApplySmearingMatrix() { //******************************************************************** if (!fSmearMatrix) { NUIS_ERR(WRN, fName << ": attempted to apply smearing matrix, but none was set"); return; } TH1D *unsmeared = (TH1D *)fMCHist->Clone(); TH1D *smeared = (TH1D *)fMCHist->Clone(); smeared->Reset(); // Loop over reconstructed bins // true = row; reco = column for (int rbin = 0; rbin < fSmearMatrix->GetNcols(); ++rbin) { // Sum up the constributions from all true bins double rBinVal = 0; // Loop over true bins for (int tbin = 0; tbin < fSmearMatrix->GetNrows(); ++tbin) { rBinVal += (*fSmearMatrix)(tbin, rbin) * unsmeared->GetBinContent(tbin + 1); } smeared->SetBinContent(rbin + 1, rBinVal); } fMCHist = (TH1D *)smeared->Clone(); return; } /* Reconfigure LOOP */ //******************************************************************** void Measurement1D::ResetAll() { //******************************************************************** fMCHist->Reset(); fMCFine->Reset(); fMCStat->Reset(); return; }; //******************************************************************** void Measurement1D::FillHistograms() { //******************************************************************** if (Signal) { NUIS_LOG(DEB, "Fill MCHist: " << fXVar << ", " << Weight); - fMCHist->Fill(fXVar, Weight); + // If it's single bin, whatever the limits on the plot are don't apply + if (fIsSingleBin){ + fMCHist->Fill(fMCHist->GetBinCenter(1), Weight); + fMCStat->Fill(fMCStat->GetBinCenter(1), 1.0); + if (fMCHist_Modes) + fMCHist_Modes->Fill(Mode, fMCHist->GetBinCenter(1), Weight); + } else { + fMCHist->Fill(fXVar, Weight); + fMCStat->Fill(fXVar, 1.0); + if (fMCHist_Modes) + fMCHist_Modes->Fill(Mode, fXVar, Weight); + } + fMCFine->Fill(fXVar, Weight); - fMCStat->Fill(fXVar, 1.0); - - if (fMCHist_Modes) - fMCHist_Modes->Fill(Mode, fXVar, Weight); + if (fMCFine_Modes) + fMCFine_Modes->Fill(Mode, fXVar, Weight); } return; }; //******************************************************************** void Measurement1D::ScaleEvents() { //******************************************************************** // Fill MCWeighted; // for (int i = 0; i < fMCHist->GetNbinsX(); i++) { // fMCWeighted->SetBinContent(i + 1, fMCHist->GetBinContent(i + 1)); // fMCWeighted->SetBinError(i + 1, fMCHist->GetBinError(i + 1)); // } // Setup Stat ratios for MC and MC Fine double *statratio = new double[fMCHist->GetNbinsX()]; for (int i = 0; i < fMCHist->GetNbinsX(); i++) { if (fMCHist->GetBinContent(i + 1) != 0) { statratio[i] = fMCHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1); } else { statratio[i] = 0.0; } } double *statratiofine = new double[fMCFine->GetNbinsX()]; for (int i = 0; i < fMCFine->GetNbinsX(); i++) { if (fMCFine->GetBinContent(i + 1) != 0) { statratiofine[i] = fMCFine->GetBinError(i + 1) / fMCFine->GetBinContent(i + 1); } else { statratiofine[i] = 0.0; } } // Scaling for raw event rates if (fIsRawEvents) { double datamcratio = fDataHist->Integral() / fMCHist->Integral(); fMCHist->Scale(datamcratio); fMCFine->Scale(datamcratio); if (fMCHist_Modes) fMCHist_Modes->Scale(datamcratio); + if (fMCFine_Modes) + fMCFine_Modes->Scale(datamcratio); + // Scaling for XSec as function of Enu } else if (fIsEnu1D) { PlotUtils::FluxUnfoldedScaling(fMCHist, GetFluxHistogram(), GetEventHistogram(), fScaleFactor, fNEvents); PlotUtils::FluxUnfoldedScaling(fMCFine, GetFluxHistogram(), GetEventHistogram(), fScaleFactor, fNEvents); if (fMCHist_Modes) { // Loop over the modes fMCHist_Modes->FluxUnfold(GetFluxHistogram(), GetEventHistogram(), fScaleFactor, fNEvents); - // PlotUtils::FluxUnfoldedScaling(fMCHist_Modes, GetFluxHistogram(), - // GetEventHistogram(), fScaleFactor, - // fNEvents); + fMCFine_Modes->FluxUnfold(GetFluxHistogram(), GetEventHistogram(), + fScaleFactor, fNEvents); } } else if (fIsNoWidth) { fMCHist->Scale(fScaleFactor); fMCFine->Scale(fScaleFactor); if (fMCHist_Modes) fMCHist_Modes->Scale(fScaleFactor); + if (fMCFine_Modes) + fMCFine_Modes->Scale(fScaleFactor); + } else if (fIsSingleBin) { + fMCHist->Scale(fScaleFactor); + fMCFine->Scale(fScaleFactor, "width"); + if (fMCHist_Modes) + fMCHist_Modes->Scale(fScaleFactor); + if (fMCFine_Modes) + fMCFine_Modes->Scale(fScaleFactor, "width"); + // Any other differential scaling } else { fMCHist->Scale(fScaleFactor, "width"); fMCFine->Scale(fScaleFactor, "width"); if (fMCHist_Modes) fMCHist_Modes->Scale(fScaleFactor, "width"); + if (fMCFine_Modes) + fMCFine_Modes->Scale(fScaleFactor, "width"); } // Proper error scaling - ROOT Freaks out with xsec weights sometimes for (int i = 0; i < fMCStat->GetNbinsX(); i++) { fMCHist->SetBinError(i + 1, fMCHist->GetBinContent(i + 1) * statratio[i]); } for (int i = 0; i < fMCFine->GetNbinsX(); i++) { fMCFine->SetBinError(i + 1, fMCFine->GetBinContent(i + 1) * statratiofine[i]); } // Clean up delete[] statratio; delete[] statratiofine; return; }; //******************************************************************** void Measurement1D::ApplyNormScale(double norm) { //******************************************************************** fCurrentNorm = norm; fMCHist->Scale(1.0 / norm); fMCFine->Scale(1.0 / norm); return; }; /* Statistic Functions - Outsources to StatUtils */ //******************************************************************** int Measurement1D::GetNDOF() { //******************************************************************** int ndof = fDataHist->GetNbinsX(); if (fMaskHist and fIsMask) ndof -= fMaskHist->Integral(); return ndof; } //******************************************************************** double Measurement1D::GetLikelihood() { //******************************************************************** // If this is for a ratio, there is no data histogram to compare to! if (fNoData || !fDataHist) return 0.; // Apply Masking to MC if Required. if (fIsMask and fMaskHist) { PlotUtils::MaskBins(fMCHist, fMaskHist); } // Sort Shape Scaling double scaleF = 0.0; // TODO Include !fIsRawEvents if (fIsShape) { // Don't renorm based on width if we are using ShapeNormDecomp if (fUseShapeNormDecomp) { if (fMCHist->Integral(1, fMCHist->GetNbinsX())) { scaleF = fDataHist->Integral(1, fDataHist->GetNbinsX()) / fMCHist->Integral(1, fMCHist->GetNbinsX()); fMCHist->Scale(scaleF); fMCFine->Scale(scaleF); } } else { if (fMCHist->Integral(1, fMCHist->GetNbinsX(), "width")) { scaleF = fDataHist->Integral(1, fDataHist->GetNbinsX(), "width") / fMCHist->Integral(1, fMCHist->GetNbinsX(), "width"); fMCHist->Scale(scaleF); fMCFine->Scale(scaleF); } } } // Likelihood Calculation double stat = 0.; if (fIsChi2) { if (fIsRawEvents) { stat = StatUtils::GetChi2FromEventRate(fDataHist, fMCHist, fMaskHist); } else if (fIsDiag) { stat = StatUtils::GetChi2FromDiag(fDataHist, fMCHist, fMaskHist); } else if (!fIsDiag and !fIsRawEvents) { stat = StatUtils::GetChi2FromCov(fDataHist, fMCHist, covar, fMaskHist, 1, 1E76, fIsWriting ? fResidualHist : NULL); if (fChi2LessBinHist && fIsWriting) { for (int xi = 0; xi < fDataHist->GetNbinsX(); ++xi) { TH1I *binmask = fMaskHist ? static_cast(fMaskHist->Clone("mask")) : new TH1I("mask", "", fDataHist->GetNbinsX(), 0, fDataHist->GetNbinsX()); binmask->SetDirectory(NULL); binmask->SetBinContent(xi + 1, 1); fChi2LessBinHist->SetBinContent( xi + 1, StatUtils::GetChi2FromCov(fDataHist, fMCHist, covar, binmask)); delete binmask; } } } } // Sort Penalty Terms if (fAddNormPen) { if (fUseShapeNormDecomp) { // if shape norm, then add the norm penalty from // https://arxiv.org/pdf/2003.00088.pdf TH1 *masked_data = StatUtils::ApplyHistogramMasking(fDataHist, fMaskHist); TH1 *masked_mc = StatUtils::ApplyHistogramMasking(fMCHist, fMaskHist); masked_mc->Scale(scaleF); NUIS_LOG(REC, "Shape Norm Decomp mcinteg: " << masked_mc->Integral() * 1E38 << ", datainteg: " << masked_data->Integral() * 1E38 << ", normerror: " << fNormError); double normpen = std::pow((masked_data->Integral() - masked_mc->Integral()) * 1E38, 2) / fNormError; masked_data->SetDirectory(NULL); delete masked_data; masked_mc->SetDirectory(NULL); delete masked_mc; NUIS_LOG(SAM, "Using Shape/Norm decomposition: Norm penalty " << normpen << " on shape penalty of " << stat); stat += normpen; } else { double penalty = (1. - fCurrentNorm) * (1. - fCurrentNorm) / (fNormError * fNormError); stat += penalty; } } // Return to normal scaling if (fIsShape) { // and !FitPar::Config().GetParB("saveshapescaling")) { fMCHist->Scale(1. / scaleF); fMCFine->Scale(1. / scaleF); } fLikelihood = stat; return stat; } /* Fake Data Functions */ //******************************************************************** void Measurement1D::SetFakeDataValues(std::string fakeOption) { //******************************************************************** // Setup original/datatrue TH1D *tempdata = (TH1D *)fDataHist->Clone(); if (!fIsFakeData) { fIsFakeData = true; // Make a copy of the original data histogram. if (!fDataOrig) fDataOrig = (TH1D *)fDataHist->Clone((fName + "_data_original").c_str()); } else { ResetFakeData(); } // Setup Inputs fFakeDataInput = fakeOption; NUIS_LOG(SAM, "Setting fake data from : " << fFakeDataInput); // From MC if (fFakeDataInput.compare("MC") == 0) { fDataHist = (TH1D *)fMCHist->Clone((fName + "_MC").c_str()); // Fake File } else { if (!fFakeDataFile) fFakeDataFile = new TFile(fFakeDataInput.c_str(), "READ"); fDataHist = (TH1D *)fFakeDataFile->Get((fName + "_MC").c_str()); } // Setup Data Hist fDataHist->SetNameTitle((fName + "_FAKE").c_str(), (fName + fPlotTitles).c_str()); // Replace Data True if (fDataTrue) delete fDataTrue; fDataTrue = (TH1D *)fDataHist->Clone(); fDataTrue->SetNameTitle((fName + "_FAKE_TRUE").c_str(), (fName + fPlotTitles).c_str()); // Make a new covariance for fake data hist. int nbins = fDataHist->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 = fDataHist->GetBinContent(i + 1) / tempdata->GetBinContent(i + 1); alpha_j = fDataHist->GetBinContent(j + 1) / tempdata->GetBinContent(j + 1); (*fFullCovar)(i, j) = alpha_i * alpha_j * (*fFullCovar)(i, j); } } // Setup Covariances if (covar) delete covar; covar = StatUtils::GetInvert(fFullCovar, true); if (fDecomp) delete fDecomp; fDecomp = StatUtils::GetDecomp(fFullCovar); delete tempdata; return; }; //******************************************************************** void Measurement1D::ResetFakeData() { //******************************************************************** if (fIsFakeData) { if (fDataHist) delete fDataHist; fDataHist = (TH1D *)fDataTrue->Clone((fSettings.GetName() + "_FKDAT").c_str()); } } //******************************************************************** void Measurement1D::ResetData() { //******************************************************************** if (fIsFakeData) { if (fDataHist) delete fDataHist; fDataHist = (TH1D *)fDataOrig->Clone((fSettings.GetName() + "_data").c_str()); } fIsFakeData = false; } //******************************************************************** void Measurement1D::ThrowCovariance() { //******************************************************************** // Take a fDecomposition and use it to throw the current dataset. // Requires fDataTrue also be set incase used repeatedly. if (!fDataTrue) fDataTrue = (TH1D *)fDataHist->Clone(); if (fDataHist) delete fDataHist; fDataHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar); return; }; //******************************************************************** void Measurement1D::ThrowDataToy() { //******************************************************************** if (!fDataTrue) fDataTrue = (TH1D *)fDataHist->Clone(); if (fMCHist) delete fMCHist; fMCHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar); } /* Access Functions */ //******************************************************************** TH1D *Measurement1D::GetMCHistogram() { //******************************************************************** if (!fMCHist) return fMCHist; std::ostringstream chi2; chi2 << std::setprecision(5) << this->GetLikelihood(); int linecolor = kRed; int linestyle = 1; int linewidth = 1; int fillcolor = 0; int fillstyle = 1001; // if (fSettings.Has("linecolor")) linecolor = fSettings.GetI("linecolor"); // if (fSettings.Has("linestyle")) linestyle = fSettings.GetI("linestyle"); // if (fSettings.Has("linewidth")) linewidth = fSettings.GetI("linewidth"); // if (fSettings.Has("fillcolor")) fillcolor = fSettings.GetI("fillcolor"); // if (fSettings.Has("fillstyle")) fillstyle = fSettings.GetI("fillstyle"); fMCHist->SetTitle(chi2.str().c_str()); fMCHist->SetLineColor(linecolor); fMCHist->SetLineStyle(linestyle); fMCHist->SetLineWidth(linewidth); fMCHist->SetFillColor(fillcolor); fMCHist->SetFillStyle(fillstyle); return fMCHist; }; //******************************************************************** TH1D *Measurement1D::GetDataHistogram() { //******************************************************************** if (!fDataHist) return fDataHist; int datacolor = kBlack; int datastyle = 1; int datawidth = 1; // if (fSettings.Has("datacolor")) datacolor = fSettings.GetI("datacolor"); // if (fSettings.Has("datastyle")) datastyle = fSettings.GetI("datastyle"); // if (fSettings.Has("datawidth")) datawidth = fSettings.GetI("datawidth"); fDataHist->SetLineColor(datacolor); fDataHist->SetLineWidth(datawidth); fDataHist->SetMarkerStyle(datastyle); return fDataHist; }; /* Write Functions */ // Save all the histograms at once //******************************************************************** void Measurement1D::Write(std::string drawOpt) { //******************************************************************** // Get Draw Options drawOpt = FitPar::Config().GetParS("drawopts"); // Write Settigns if (drawOpt.find("SETTINGS") != std::string::npos) { fSettings.Set("#chi^{2}", fLikelihood); fSettings.Set("NDOF", this->GetNDOF()); fSettings.Set("#chi^{2}/NDOF", fLikelihood / this->GetNDOF()); fSettings.Write(); } // Write Data/MC if (drawOpt.find("DATA") != std::string::npos) GetDataList().at(0)->Write(); if (drawOpt.find("MC") != std::string::npos) { GetMCList().at(0)->Write(); if ((fEvtRateScaleFactor != 0xdeadbeef) && GetMCList().at(0)) { TH1D *PredictedEvtRate = static_cast(GetMCList().at(0)->Clone()); PredictedEvtRate->Scale(fEvtRateScaleFactor); PredictedEvtRate->GetYaxis()->SetTitle("Predicted event rate"); PredictedEvtRate->Write(); } } // Write Fine Histogram if (fSaveFine && drawOpt.find("FINE") != std::string::npos) GetFineList().at(0)->Write(); // Write Weighted Histogram if (drawOpt.find("WEIGHTS") != std::string::npos && fMCWeighted) fMCWeighted->Write(); // Save Flux/Evt if no event manager if (!FitPar::Config().GetParB("EventManager")) { if (drawOpt.find("FLUX") != std::string::npos && GetFluxHistogram()) GetFluxHistogram()->Write(); if (drawOpt.find("EVT") != std::string::npos && GetEventHistogram()) GetEventHistogram()->Write(); if (drawOpt.find("XSEC") != std::string::npos && GetEventHistogram()) GetXSecHistogram()->Write(); } // Write Mask if (fIsMask && (drawOpt.find("MASK") != std::string::npos)) { fMaskHist->Write(); } // Write Covariances if (drawOpt.find("COV") != std::string::npos && fFullCovar) { PlotUtils::GetFullCovarPlot(fFullCovar, fSettings.GetName())->Write(); } if (drawOpt.find("INVCOV") != std::string::npos && covar) { PlotUtils::GetInvCovarPlot(covar, fSettings.GetName())->Write(); } if (drawOpt.find("DECOMP") != std::string::npos && fDecomp) { PlotUtils::GetDecompCovarPlot(fDecomp, fSettings.GetName())->Write(); } // // Likelihood residual plots // if (drawOpt.find("RESIDUAL") != std::string::npos) { // WriteResidualPlots(); // } // Ratio and Shape Plots if (drawOpt.find("RATIO") != std::string::npos) { WriteRatioPlot(); } if (drawOpt.find("SHAPE") != std::string::npos) { WriteShapePlot(); if (drawOpt.find("RATIO") != std::string::npos) WriteShapeRatioPlot(); } // // RATIO // if (drawOpt.find("CANVMC") != std::string::npos) { // TCanvas* c1 = WriteMCCanvas(fDataHist, fMCHist); // c1->Write(); // delete c1; // } // // PDG // if (drawOpt.find("CANVPDG") != std::string::npos && fMCHist_Modes) { // TCanvas* c2 = WritePDGCanvas(fDataHist, fMCHist, fMCHist_Modes); // c2->Write(); // delete c2; // } if (fIsChi2 && !fIsDiag) { fResidualHist = (TH1D *)fMCHist->Clone((fName + "_RESIDUAL").c_str()); fResidualHist->GetYaxis()->SetTitle("#Delta#chi^{2}"); fResidualHist->Reset(); fChi2LessBinHist = (TH1D *)fMCHist->Clone((fName + "_Chi2NMinusOne").c_str()); fChi2LessBinHist->GetYaxis()->SetTitle("Total #chi^{2} without bin_{i}"); fChi2LessBinHist->Reset(); fIsWriting = true; (void)GetLikelihood(); fIsWriting = false; fResidualHist->Write((fName + "_RESIDUAL").c_str()); fChi2LessBinHist->Write((fName + "_Chi2NMinusOne").c_str()); } // Write Extra Histograms AutoWriteExtraTH1(); WriteExtraHistograms(); // Returning NUIS_LOG(SAM, "Written Histograms: " << fName); return; } //******************************************************************** void Measurement1D::WriteRatioPlot() { //******************************************************************** // Setup mc data ratios TH1D *dataRatio = (TH1D *)fDataHist->Clone((fName + "_data_RATIO").c_str()); TH1D *mcRatio = (TH1D *)fMCHist->Clone((fName + "_MC_RATIO").c_str()); // Extra MC Data Ratios for (int i = 0; i < mcRatio->GetNbinsX(); i++) { dataRatio->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) / fMCHist->GetBinContent(i + 1)); dataRatio->SetBinError(i + 1, fDataHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1)); mcRatio->SetBinContent(i + 1, fMCHist->GetBinContent(i + 1) / fMCHist->GetBinContent(i + 1)); mcRatio->SetBinError(i + 1, fMCHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1)); } // Write ratios mcRatio->Write(); dataRatio->Write(); delete mcRatio; delete dataRatio; } //******************************************************************** void Measurement1D::WriteShapePlot() { //******************************************************************** TH1D *mcShape = (TH1D *)fMCHist->Clone((fName + "_MC_SHAPE").c_str()); TH1D *dataShape = (TH1D *)fDataHist->Clone((fName + "_data_SHAPE").c_str()); // Set the shape covariance to calculate the chi2 if (!fShapeCovar) SetShapeCovar(); // Don't check error if (fShapeCovar) StatUtils::SetDataErrorFromCov(dataShape, fShapeCovar, 1E-38, false); double shapeScale = 1.0; if (fIsRawEvents) { shapeScale = fDataHist->Integral() / fMCHist->Integral(); } else { shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width"); } mcShape->Scale(shapeScale); std::stringstream ss; ss << shapeScale; mcShape->SetTitle(ss.str().c_str()); mcShape->SetLineWidth(3); mcShape->SetLineStyle(7); mcShape->Write(); dataShape->Write(); delete mcShape; } //******************************************************************** void Measurement1D::WriteShapeRatioPlot() { //******************************************************************** // Get a mcshape histogram TH1D *mcShape = (TH1D *)fMCHist->Clone((fName + "_MC_SHAPE").c_str()); double shapeScale = 1.0; if (fIsRawEvents) { shapeScale = fDataHist->Integral() / fMCHist->Integral(); } else { shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width"); } mcShape->Scale(shapeScale); // Create shape ratio histograms TH1D *mcShapeRatio = (TH1D *)mcShape->Clone((fName + "_MC_SHAPE_RATIO").c_str()); TH1D *dataShapeRatio = (TH1D *)fDataHist->Clone((fName + "_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); mcShapeRatio->Write(); dataShapeRatio->Write(); delete mcShapeRatio; delete dataShapeRatio; } //// CRAP TO BE REMOVED //******************************************************************** void Measurement1D::SetupMeasurement(std::string inputfile, std::string type, FitWeight *rw, std::string fkdt) { //******************************************************************** nuiskey samplekey = Config::CreateKey("sample"); samplekey.Set("name", fName); samplekey.Set("type", type); samplekey.Set("input", inputfile); fSettings = LoadSampleSettings(samplekey); // 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 fIsRawEvents = false; if ((fName.find("Evt") != std::string::npos) && fIsRawEvents == false) { fIsRawEvents = true; NUIS_LOG(SAM, "Found event rate measurement but fIsRawEvents == false!"); NUIS_LOG(SAM, "Overriding this and setting fIsRawEvents == true!"); } fIsEnu1D = false; if (fName.find("XSec_1DEnu") != std::string::npos) { fIsEnu1D = true; NUIS_LOG(SAM, "Found XSec Enu measurement, applying flux integrated scaling, " "not flux averaged!"); } if (fIsEnu1D && fIsRawEvents) { NUIS_ERR(FTL, "Found 1D Enu XSec distribution AND fIsRawEvents, is this " "really correct?!"); NUIS_ERR(FTL, "Check experiment constructor for " << fName << " and correct this!"); NUIS_ERR(FTL, "I live in " << __FILE__ << ":" << __LINE__); throw; } fRW = rw; if (!fInput and !fIsJoint) SetupInputs(inputfile); // Set Default Options SetFitOptions(fDefaultTypes); // Set Passed Options SetFitOptions(type); // Still adding support for flat flux inputs // // Set Enu Flux Scaling // if (isFlatFluxFolding) this->Input()->ApplyFluxFolding( // this->defaultFluxHist ); // FinaliseMeasurement(); } //******************************************************************** void Measurement1D::SetupDefaultHist() { //******************************************************************** // Setup fMCHist fMCHist = (TH1D *)fDataHist->Clone(); fMCHist->SetNameTitle((fName + "_MC").c_str(), (fName + "_MC" + fPlotTitles).c_str()); // Setup fMCFine Int_t nBins = fMCHist->GetNbinsX(); fMCFine = new TH1D( (fName + "_MC_FINE").c_str(), (fName + "_MC_FINE" + fPlotTitles).c_str(), nBins * 6, fMCHist->GetBinLowEdge(1), fMCHist->GetBinLowEdge(nBins + 1)); fMCStat = (TH1D *)fMCHist->Clone(); fMCStat->Reset(); fMCHist->Reset(); fMCFine->Reset(); // Setup the NEUT Mode Array PlotUtils::CreateNeutModeArray((TH1D *)fMCHist, (TH1 **)fMCHist_PDG); PlotUtils::ResetNeutModeArray((TH1 **)fMCHist_PDG); // Setup bin masks using sample name if (fIsMask) { std::string maskloc = FitPar::Config().GetParDIR(fName + ".mask"); if (maskloc.empty()) { maskloc = FitPar::GetDataBase() + "/masks/" + fName + ".mask"; } SetBinMask(maskloc); } fMCHist_Modes = new TrueModeStack((fName + "_MODES").c_str(), ("True Channels"), fMCHist); SetAutoProcessTH1(fMCHist_Modes, kCMD_Reset, kCMD_Norm, kCMD_Write); + fMCFine_Modes = + new TrueModeStack((fName + "_MODES").c_str(), ("True Channels"), fMCHist); + SetAutoProcessTH1(fMCFine_Modes, kCMD_Reset, kCMD_Norm, kCMD_Write); + return; } //******************************************************************** void Measurement1D::SetDataValues(std::string dataFile) { //******************************************************************** // Override this function if the input file isn't in a suitable format NUIS_LOG(SAM, "Reading data from: " << dataFile.c_str()); fDataHist = PlotUtils::GetTH1DFromFile(dataFile, (fName + "_data"), fPlotTitles); fDataTrue = (TH1D *)fDataHist->Clone(); // Number of data points is number of bins fNDataPointsX = fDataHist->GetXaxis()->GetNbins(); return; }; //******************************************************************** void Measurement1D::SetDataFromDatabase(std::string inhistfile, std::string histname) { //******************************************************************** NUIS_LOG(SAM, "Filling histogram from " << inhistfile << "->" << histname); fDataHist = PlotUtils::GetTH1DFromRootFile( (GeneralUtils::GetTopLevelDir() + "/data/" + inhistfile), histname); fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_data").c_str()); return; }; //******************************************************************** void Measurement1D::SetDataFromFile(std::string inhistfile, std::string histname) { //******************************************************************** NUIS_LOG(SAM, "Filling histogram from " << inhistfile << "->" << histname); fDataHist = PlotUtils::GetTH1DFromRootFile((inhistfile), histname); fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_data").c_str()); 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 *fFullCovarPlot = new TH2D(); std::string covName = ""; std::string covOption = FitPar::Config().GetParS("thrown_covariance"); if (fIsShape || fIsFree) covName = "shp_"; if (fIsDiag) covName += "diag"; else covName += "full"; covarPlot = (TH2D *)tempFile->Get((covName + "cov").c_str()); if (!covOption.compare("SUB")) fFullCovarPlot = (TH2D *)tempFile->Get((covName + "cov").c_str()); else if (!covOption.compare("FULL")) fFullCovarPlot = (TH2D *)tempFile->Get("fullcov"); else { NUIS_ERR(WRN, "Incorrect thrown_covariance option in parameters."); } int dim = int(fDataHist->GetNbinsX()); //-this->masked->Integral()); int covdim = int(fDataHist->GetNbinsX()); this->covar = new TMatrixDSym(dim); fFullCovar = new TMatrixDSym(dim); fDecomp = 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); (*fFullCovar)(row, column) = fFullCovarPlot->GetBinContent(i + 1, j + 1); column++; } column = 0; row++; } // Set bin errors on data if (!fIsDiag) { StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar); } // Get Deteriminant and inverse matrix // fCovDet = this->covar->Determinant(); TDecompSVD LU = TDecompSVD(*this->covar); this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), ""); return; }; //******************************************************************** // Sets the covariance matrix from a provided file in a text format // scale is a multiplicative pre-factor to apply in the case where the // covariance is given in some unit (e.g. 1E-38) void Measurement1D::SetCovarMatrixFromText(std::string covarFile, int dim, double scale) { //******************************************************************** // Make a counter to track the line number int row = 0; std::string line; std::ifstream covarread(covarFile.c_str(), std::ifstream::in); this->covar = new TMatrixDSym(dim); fFullCovar = new TMatrixDSym(dim); if (covarread.is_open()) { NUIS_LOG(SAM, "Reading covariance matrix from file: " << covarFile); } else { NUIS_ABORT("Covariance matrix provided is incorrect: " << covarFile); } // Loop over the lines in the file while (std::getline(covarread >> std::ws, line, '\n')) { int column = 0; // Loop over entries and insert them into matrix std::vector entries = GeneralUtils::ParseToDbl(line, " "); if (entries.size() <= 1) { NUIS_ERR(WRN, "SetCovarMatrixFromText -> Covariance matrix only has <= 1 " "entries on this line: " << row); } for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { (*covar)(row, column) = *iter; (*fFullCovar)(row, column) = *iter; column++; } row++; } covarread.close(); // Scale the actualy covariance matrix by some multiplicative factor (*fFullCovar) *= scale; // Robust matrix inversion method TDecompSVD LU = TDecompSVD(*this->covar); // THIS IS ACTUALLY THE INVERSE COVARIANCE MATRIXA AAAAARGH delete this->covar; this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), ""); // Now need to multiply by the scaling factor // If the covariance (*this->covar) *= 1. / (scale); return; }; //******************************************************************** void Measurement1D::SetCovarMatrixFromCorrText(std::string corrFile, int dim) { //******************************************************************** // Make a counter to track the line number int row = 0; std::string line; std::ifstream corr(corrFile.c_str(), std::ifstream::in); this->covar = new TMatrixDSym(dim); this->fFullCovar = new TMatrixDSym(dim); if (corr.is_open()) { NUIS_LOG(SAM, "Reading and converting correlation matrix from file: " << corrFile); } else { NUIS_ABORT("Correlation matrix provided is incorrect: " << corrFile); } while (std::getline(corr >> std::ws, line, '\n')) { int column = 0; // Loop over entries and insert them into matrix // Multiply by the errors to get the covariance, rather than the correlation // matrix std::vector entries = GeneralUtils::ParseToDbl(line, " "); for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { double val = (*iter) * this->fDataHist->GetBinError(row + 1) * 1E38 * this->fDataHist->GetBinError(column + 1) * 1E38; if (val == 0) { NUIS_ABORT("Found a zero value in the covariance matrix, assuming " "this is an error!"); } (*this->covar)(row, column) = val; (*this->fFullCovar)(row, column) = val; column++; } row++; } // Robust matrix inversion method TDecompSVD LU = TDecompSVD(*this->covar); delete this->covar; this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), ""); return; }; //******************************************************************** // FullUnits refers to if we have "real" unscaled units in the covariance // matrix, e.g. 1E-76. If this is the case we need to scale it so that the chi2 // contribution is correct NUISANCE internally assumes the covariance matrix has // units of 1E76 void Measurement1D::SetCovarFromDataFile(std::string covarFile, std::string covName, bool FullUnits) { //******************************************************************** NUIS_LOG(SAM, "Getting covariance from " << covarFile << "->" << covName); TFile *tempFile = new TFile(covarFile.c_str(), "READ"); TH2D *covPlot = (TH2D *)tempFile->Get(covName.c_str()); covPlot->SetDirectory(0); // Scale the covariance matrix if it comes in normal units if (FullUnits) { covPlot->Scale(1.E76); } int dim = covPlot->GetNbinsX(); fFullCovar = new TMatrixDSym(dim); for (int i = 0; i < dim; i++) { for (int j = 0; j < dim; j++) { (*fFullCovar)(i, j) = covPlot->GetBinContent(i + 1, j + 1); } } this->covar = (TMatrixDSym *)fFullCovar->Clone(); fDecomp = (TMatrixDSym *)fFullCovar->Clone(); TDecompSVD LU = TDecompSVD(*this->covar); this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), ""); TDecompChol LUChol = TDecompChol(*fDecomp); LUChol.Decompose(); fDecomp = new TMatrixDSym(dim, LU.GetU().GetMatrixArray(), ""); return; }; // //******************************************************************** // void Measurement1D::SetBinMask(std::string maskFile) { // //******************************************************************** // // Create a mask histogram. // int nbins = fDataHist->GetNbinsX(); // fMaskHist = // new TH1I((fName + "_fMaskHist").c_str(), // (fName + "_fMaskHist; Bin; Mask?").c_str(), nbins, 0, nbins); // std::string line; // std::ifstream mask(maskFile.c_str(), std::ifstream::in); // if (mask.is_open()) // LOG(SAM) << "Reading bin mask from file: " << maskFile << std::endl; // else // LOG(FTL) << " Cannot find mask file." << std::endl; // while (std::getline(mask >> std::ws, line, '\n')) { // std::vector entries = GeneralUtils::ParseToInt(line, " "); // // Skip lines with poorly formatted lines // if (entries.size() < 2) { // LOG(WRN) << "Measurement1D::SetBinMask(), couldn't parse line: " << // line // << std::endl; // continue; // } // // The first index should be the bin number, the second should be the // mask // // value. // fMaskHist->SetBinContent(entries[0], entries[1]); // } // // Set masked data bins to zero // PlotUtils::MaskBins(fDataHist, fMaskHist); // return; // } // //******************************************************************** // void Measurement1D::GetBinContents(std::vector& cont, // std::vector& err) { // //******************************************************************** // // Return a vector of the main bin contents // for (int i = 0; i < fMCHist->GetNbinsX(); i++) { // cont.push_back(fMCHist->GetBinContent(i + 1)); // err.push_back(fMCHist->GetBinError(i + 1)); // } // 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"); // fFluxHist = // new TH1D((fName + "_flux").c_str(), (fName + "; E_{#nu} (GeV)").c_str(), // f.GetN() - 1, minE, maxE); // Double_t* yVal = f.GetY(); // for (int i = 0; i < fFluxHist->GetNbinsX(); ++i) // fFluxHist->SetBinContent(i + 1, yVal[i] * fluxNorm); // }; // //******************************************************************** // double Measurement1D::TotalIntegratedFlux(std::string intOpt, double low, // double high) { // //******************************************************************** // if (fInput->GetType() == kGiBUU) { // return 1.0; // } // // The default case of low = -9999.9 and high = -9999.9 // if (low == -9999.9) low = this->EnuMin; // if (high == -9999.9) high = this->EnuMax; // int minBin = fFluxHist->GetXaxis()->FindBin(low); // int maxBin = fFluxHist->GetXaxis()->FindBin(high); // // Get integral over custom range // double integral = fFluxHist->Integral(minBin, maxBin + 1, intOpt.c_str()); // return integral; // }; diff --git a/src/FitBase/Measurement1D.h b/src/FitBase/Measurement1D.h index ef53680..cbb7ee2 100644 --- a/src/FitBase/Measurement1D.h +++ b/src/FitBase/Measurement1D.h @@ -1,659 +1,660 @@ // Copyright 2016-2021 L. Pickering, P towell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #ifndef MEASUREMENT_1D_H_SEEN #define MEASUREMENT_1D_H_SEEN /*! * \addtogroup FitBase * @{ */ #include #include #include #include #include #include #include #include // ROOT includes #include #include #include #include #include #include #include #include #include #include #include #include // External data fit includes #include "FitEvent.h" #include "FitUtils.h" #include "MeasurementBase.h" #include "PlotUtils.h" #include "StatUtils.h" #include "SignalDef.h" #include "MeasurementVariableBox.h" #include "MeasurementVariableBox1D.h" namespace NUISANCE { namespace FitBase { } } //******************************************************************** /// 1D Measurement base class. Histogram handling is done in this base layer. class Measurement1D : public MeasurementBase { //******************************************************************** public: /* Constructor/Deconstuctor */ Measurement1D(void); virtual ~Measurement1D(void); /* Setup Functions */ /// \brief Setup all configs once initialised /// /// Should be called after all configs have been setup inside fSettings container. /// Handles the processing of inputs and setting up of types. /// Replaces the old 'SetupMeasurement' function. void FinaliseSampleSettings(); /// \brief Creates the 1D data distribution given the binning provided. virtual void CreateDataHistogram(int dimx, double* binx); /// \brief Read 1D data inputs from a text file. /// /// Inputfile should have the format: \n /// low_binedge_1 bin_content_1 bin_error_1 \n /// low_binedge_2 bin_content_2 bin_error_2 \n /// .... .... .... \n /// high_bin_edge_N 0.0 0.0 virtual void SetDataFromTextFile(std::string datafile); /// \brief Read 1D data inputs from a root file. /// /// inhistfile specifies the path to the root file /// histname specifies the name of the histogram. /// /// If no histogram name is given the inhistfile value /// is automatically parsed with ';' so that: \n /// 'myhistfile.root;myhistname' \n /// will also work. virtual void SetDataFromRootFile(std::string inhistfile, std::string histname = ""); /// \brief Setup a default empty data histogram /// /// Only used for flattree creators. virtual void SetEmptyData(); /// \brief Set data bin errors to sqrt(entries) /// /// \warning REQUIRES DATA HISTOGRAM TO BE SET FIRST /// /// Sets the data errors as the sqrt of the bin contents /// Should be use for counting experiments virtual void SetPoissonErrors(); /// \brief Make diagonal covariance from data /// /// \warning If no histogram passed, data must be setup first! /// Setup the covariance inputs by taking the data histogram /// errors and setting up a diagonal covariance matrix. /// /// If no data is supplied, fDataHist is used if already set. virtual void SetCovarFromDiagonal(TH1D* data = NULL); /// \brief Read the data covariance from a text file. /// /// Inputfile should have the format: \n /// covariance_11 covariance_12 covariance_13 ... \n /// covariance_21 covariance_22 covariance_23 ... \n /// ... ... ... ... \n /// /// If no dimensions are given, it is assumed from the number /// entries in the first line of covfile. virtual void SetCovarFromTextFile(std::string covfile, int dim = -1); virtual void SetCovarFromMultipleTextFiles(std::string covfiles, int dim = -1); /// \brief Read the data covariance from a ROOT file. /// /// - covfile specifies the full path to the file /// - histname specifies the name of the covariance object. Both TMatrixDSym and TH2D are supported. /// /// If no histogram name is given the inhistfile value /// is automatically parsed with ; so that: \n /// mycovfile.root;myhistname \n /// will also work. virtual void SetCovarFromRootFile(std::string covfile, std::string histname=""); /// \brief Read the inverted data covariance from a text file. /// /// Inputfile should have similar format to that shown /// in SetCovarFromTextFile. /// /// If no dimensions are given, it is assumed from the number /// entries in the first line of covfile. virtual void SetCovarInvertFromTextFile(std::string covfile, int dim = -1); /// \brief Read the inverted data covariance from a ROOT file. /// /// Inputfile should have similar format to that shown /// in SetCovarFromRootFile. /// /// If no histogram name is given the inhistfile value /// is automatically parsed with ; so that: \n /// mycovfile.root;myhistname \n /// will also work. virtual void SetCovarInvertFromRootFile(std::string covfile, std::string histname=""); /// \brief Read the data correlations from a text file. /// /// \warning REQUIRES DATA HISTOGRAM TO BE SET FIRST /// /// Inputfile should have similar format to that shown /// in SetCovarFromTextFile. /// /// If no dimensions are given, it is assumed from the number /// entries in the first line of covfile. virtual void SetCorrelationFromTextFile(std::string covfile, int dim = -1); /// \brief Read the data correlations from multiple text files. /// /// \warning REQUIRES DATA HISTOGRAM TO BE SET FIRST /// /// Inputfile should have similar format to that shown /// in SetCovarFromTextFile. /// /// If no dimensions are given, it is assumed from the number /// entries in the first line of the first corrfile. virtual void SetCorrelationFromMultipleTextFiles(std::string corrfiles, int dim = -1); /// \brief Read the data correlations from a ROOT file. /// /// \warning REQUIRES DATA TO BE SET FIRST /// /// Inputfile should have similar format to that shown /// in SetCovarFromRootFile. /// /// If no histogram name is given the inhistfile value /// is automatically parsed with ; so that: \n /// mycovfile.root;myhistname \n /// will also work. virtual void SetCorrelationFromRootFile(std::string covfile, std::string histname=""); /// \brief Read the cholesky decomposed covariance from a text file and turn it into a covariance /// /// Inputfile should have similar format to that shown /// in SetCovarFromTextFile. /// /// If no dimensions are given, it is assumed from the number /// entries in the first line of covfile. virtual void SetCholDecompFromTextFile(std::string covfile, int dim = -1); /// \brief Read the cholesky decomposed covariance from a ROOT file and turn it into a covariance /// /// Inputfile should have similar format to that shown /// in SetCovarFromRootFile. /// /// If no histogram name is given the inhistfile value /// is automatically parsed with ; so that: \n /// mycovfile.root;myhistname \n /// will also work. virtual void SetCholDecompFromRootFile(std::string covfile, std::string histname=""); /// \brief Try to extract a shape-only matrix from the existing covariance virtual void SetShapeCovar(); /// \brief Scale the data by some scale factor virtual void ScaleData(double scale); /// \brief Scale the data error bars by some scale factor virtual void ScaleDataErrors(double scale); /// \brief Scale the covariaince and its invert/decomp by some scale factor. virtual void ScaleCovar(double scale); /// \brief Setup a bin masking histogram and apply masking to data /// /// \warning REQUIRES DATA HISTOGRAM TO BE SET FIRST /// /// Reads in a list of bins in a text file to be masked. Format is: \n /// bin_index_1 1 \n /// bin_index_2 1 \n /// bin_index_3 1 \n /// /// If 0 is given then a bin entry will NOT be masked. So for example: \n\n /// 1 1 \n /// 2 1 \n /// 3 0 \n /// 4 1 \n\n /// Will mask only the 1st, 2nd, and 4th bins. /// /// Masking can be turned on by specifiying the MASK option when creating a sample. /// When this is passed NUISANCE will look in the following locations for the mask file: /// - FitPar::Config().GetParS(fName + ".mask") /// - "data/masks/" + fName + ".mask"; virtual void SetBinMask(std::string maskfile); /// \brief Set the current fit options from a string. /// /// This is called twice for each sample, once to set the default /// and once to set the current setting (if anything other than default given) /// /// For this to work properly it requires the default and allowed types to be /// set correctly. These should be specified as a string listing options. /// /// To split up options so that NUISANCE can automatically detect ones that /// are conflicting. Any options separated with the '/' symbol are non conflicting /// and can be given together, whereas any separated with the ',' symbol cannot /// be specified by the end user at the same time. /// /// Default Type Examples: /// - DIAG/FIX = Default option will be a diagonal covariance, with FIXED norm. /// - MASK/SHAPE = Default option will be a masked hist, with SHAPE always on. /// /// Allowed Type examples: /// - 'FULL/DIAG/NORM/MASK' = Any of these options can be specified. /// - 'FULL,FREE,SHAPE/MASK/NORM' = User can give either FULL, FREE, or SHAPE as on option. /// MASK and NORM can also be included as options. virtual void SetFitOptions(std::string opt); /// \brief Final constructor setup /// \warning Should be called right at the end of the constructor. /// /// Contains a series of checks to ensure the data and inputs have been setup. /// Also creates the MC histograms needed for fitting. void FinaliseMeasurement(); /* Smearing */ /// \brief Read in smearing matrix from file /// /// Set the smearing matrix from a text file given the size of the matrix virtual void SetSmearingMatrix(std::string smearfile, int truedim, int recodim); /// \brief Apply smearing to MC true to get MC reco /// /// Apply smearing matrix to fMCHist using fSmearingMatrix virtual void ApplySmearingMatrix(void); /* Reconfigure Functions */ /// \brief Create a Measurement1D box /// /// Creates a new 1D variable box containing just fXVar. /// /// This box is the bare minimum required by the JointFCN when /// running fast reconfigures during a routine. /// If for some reason a sample needs extra variables to be saved then /// it should override this function creating its own MeasurementVariableBox /// that contains the extra variables. virtual MeasurementVariableBox* CreateBox() {return new MeasurementVariableBox1D();}; /// \brief Reset all MC histograms /// /// Resets all standard histograms and those registered to auto /// process to zero. /// /// If extra histograms are not included in auto processing, then they must be reset /// by overriding this function and doing it manually if required. virtual void ResetAll(void); /// \brief Fill MC Histograms from XVar /// /// Fill standard histograms using fXVar, Weight read from the variable box. /// /// WARNING : Any extra MC histograms need to be filled by overriding this function, /// even if they have been set to auto process. virtual void FillHistograms(void); // \brief Convert event rates to final histogram /// /// Apply standard scaling procedure to standard mc histograms to convert from /// raw events to xsec prediction. /// /// If any distributions have been set to auto process /// that is done during this function call, and a differential xsec is assumed. /// If that is not the case this function must be overriden. virtual void ScaleEvents(void); /// \brief Scale MC by a factor=1/norm /// /// Apply a simple normalisation scaling if the option FREE or a norm_parameter /// has been specified in the NUISANCE routine. virtual void ApplyNormScale(double norm); /* Statistical Functions */ /// \brief Get Number of degrees of freedom /// /// Returns the number bins inside the data histogram accounting for /// any bin masking applied. virtual int GetNDOF(void); /// \brief Return Data/MC Likelihood at current state /// /// Returns the likelihood of the data given the current MC prediction. /// Diferent likelihoods definitions are used depending on the FitOptions. virtual double GetLikelihood(void); /* Fake Data */ /// \brief Set the fake data values from either a file, or MC /// /// - Setting from a file "path": \n /// When reading from a file the full path must be given to a standard /// nuisance output. The standard MC histogram should have a name that matches /// this sample for it to be read in. /// \n\n /// - Setting from "MC": \n /// If the MC option is given the current MC prediction is used as fake data. virtual void SetFakeDataValues(std::string fakeOption); /// \brief Reset fake data back to starting fake data /// /// Reset the fake data back to original fake data (Reset back to before /// ThrowCovariance was first called) virtual void ResetFakeData(void); /// \brief Reset fake data back to original data /// /// Reset the data histogram back to the true original dataset for this sample /// before any fake data was defined. virtual void ResetData(void); /// \brief Generate fake data by throwing the covariance. /// /// Can be used on fake MC data or just the original dataset. /// Call ResetFakeData or ResetData to return to values before the throw. virtual void ThrowCovariance(void); /// \brief Throw the data by its assigned errors and assign this to MC /// /// Used when creating data toys by assign the MC to this thrown data /// so that the likelihood is calculated between data and thrown data virtual void ThrowDataToy(void); /* Access Functions */ /// \brief Returns nicely formatted MC Histogram /// /// Format options can also be given in the samplesettings: /// - linecolor /// - linestyle /// - linewidth /// - fillcolor /// - fillstyle /// /// So to have a sample line colored differently in the xml cardfile put: \n /// virtual TH1D* GetMCHistogram(void); /// \brief Returns nicely formatted data Histogram /// /// Format options can also be given in the samplesettings: /// - datacolor /// - datastyle /// - datawidth /// /// So to have a sample data colored differently in the xml cardfile put: \n /// virtual TH1D* GetDataHistogram(void); /// \brief Returns a list of all MC histograms. /// /// Override this if you have extra histograms that need to be /// accessed outside of the Measurement1D class. inline virtual std::vector GetMCList(void) { return std::vector(1, GetMCHistogram()); } /// \brief Returns a list of all Data histograms. /// /// Override this if you have extra histograms that need to be /// accessed outside of the Measurement1D class. inline virtual std::vector GetDataList(void) { return std::vector(1, GetDataHistogram()); } /// \brief Returns a list of all Mask histograms. /// /// Override this if you have extra histograms that need to be /// accessed outside of the Measurement1D class. inline virtual std::vector GetMaskList(void) { return std::vector(1, fMaskHist); }; /// \brief Returns a list of all Fine histograms. /// /// Override this if you have extra histograms that need to be /// accessed outside of the Measurement1D class. inline virtual std::vector GetFineList(void) { return std::vector(1, fMCFine); }; /* Write Functions */ /// \brief Save the current state to the current TFile directory \n /// /// Data/MC are both saved by default. /// A range of other histograms can be saved by setting the /// config option 'drawopts'. /// /// Possible options: \n /// - FINE = Write Fine Histogram \n /// - WEIGHTS = Write Weighted MC Histogram (before scaling) \n /// - FLUX = Write Flux Histogram from MC Input \n /// - EVT = Write Event Histogram from MC Input \n /// - XSEC = Write XSec Histogram from MC Input \n /// - MASK = Write Mask Histogram \n /// - COV = Write Covariance Histogram \n /// - INVCOV = Write Inverted Covariance Histogram \n /// - DECMOP = Write Decomp. Covariance Histogram \n /// - RESIDUAL= Write Resudial Histograms \n /// - RATIO = Write Data/MC Ratio Histograms \n /// - SHAPE = Write MC Shape Histograms norm. to Data \n /// - CANVMC = Build MC Canvas Showing Data, MC, Shape \n /// - MODES = Write PDG Stack \n /// - CANVPDG = Build MC Canvas Showing Data, PDGStack \n /// /// So to save a range of these in parameters/config.xml set: \n /// virtual void Write(std::string drawOpt); virtual void WriteRatioPlot(); virtual void WriteShapePlot(); virtual void WriteShapeRatioPlot(); //* // OLD DEFUNCTIONS // /// OLD FUNCTION virtual void SetupMeasurement(std::string input, std::string type, FitWeight* rw, std::string fkdt); /// OLD FUNCTION virtual void SetupDefaultHist(void); /// OLD FUNCTION virtual void SetDataValues(std::string dataFile); /// OLD FUNCTION virtual void SetDataFromFile(std::string inhistfile, std::string histname); /// OLD FUNCTION virtual void SetDataFromDatabase(std::string inhistfile, std::string histname); /// OLD FUNCTION virtual void SetCovarMatrix(std::string covarFile); /// OLD FUNCTION virtual void SetCovarMatrixFromText(std::string covarFile, int dim, double scale = 1.0); /// OLD FUNCTION virtual void SetCovarMatrixFromCorrText(std::string covarFile, int dim); /// OLD FUNCTION virtual void SetCovarFromDataFile(std::string covarFile, std::string covName, bool FullUnits = false); /// OLD FUNCTION // virtual THStack GetModeStack(void); protected: // Data TH1D* fDataHist; ///< default data histogram TH1D* fDataOrig; ///< histogram to store original data before throws. TH1D* fDataTrue; ///< histogram to store true dataset std::string fPlotTitles; ///< Plot title x and y for the histograms // MC TH1D* fMCHist; ///< default MC Histogram used in the chi2 fits TH1D* fMCFine; ///< finely binned MC histogram TH1D* fMCStat; ///< histogram with unweighted events to properly calculate TH1D* fMCWeighted; ///< Weighted histogram before xsec scaling TH1I* fMaskHist; ///< Mask histogram for neglecting specific bins TMatrixD* fSmearMatrix; ///< Smearing matrix (note, this is not symmetric) TH1D *fResidualHist; TH1D *fChi2LessBinHist; TrueModeStack* fMCHist_Modes; ///< Optional True Mode Stack - + TrueModeStack* fMCFine_Modes; ///< Optional True Mode Stack // Statistical TMatrixDSym* covar; ///< Inverted Covariance TMatrixDSym* fFullCovar; ///< Full Covariance TMatrixDSym* fDecomp; ///< Decomposed Covariance TMatrixDSym* fCorrel; ///< Correlation Matrix TMatrixDSym* fShapeCovar; ///< Shape-only covariance TMatrixDSym* fShapeDecomp; ///< Decomposed shape-only covariance TMatrixDSym* fShapeInvert; ///< Inverted shape-only covariance TMatrixDSym* fCovar; ///< New FullCovar TMatrixDSym* fInvert; ///< New covar double fNormError; ///< Sample norm error double fLikelihood; ///< Likelihood value // Fake Data bool fIsFakeData; ///< Flag: is the current data fake from MC std::string fFakeDataInput; ///< Input fake data file path TFile* fFakeDataFile; ///< Input fake data file // Fit specific flags std::string fFitType; ///< Current fit type std::string fAllowedTypes; ///< Fit Types Possible std::string fDefaultTypes; ///< Starting Default Fit Types bool fIsShape; ///< Flag : Perform Shape-only fit bool fUseShapeNormDecomp; bool fIsFree; ///< Flag : Perform normalisation free fit bool fIsDiag; ///< Flag : only include uncorrelated diagonal errors bool fIsMask; ///< Flag : Apply bin masking bool fIsRawEvents; ///< Flag : Are bin contents just event rates bool fIsEnu1D; ///< Flag : Perform Flux Unfolded Scaling bool fIsChi2SVD; ///< Flag : Use alternative Chi2 SVD Method (Do not use) bool fAddNormPen; ///< Flag : Add a normalisation penalty term to the chi2. bool fIsFix; ///< Flag : keeping norm fixed bool fIsFull; ///< Flag : using full covariaince bool fIsDifXSec; ///< Flag : creating a dif xsec bool fIsChi2; ///< Flag : using Chi2 over LL methods bool fIsSmeared; ///< Flag : Apply smearing? + bool fIsSingleBin; ///< Flag : Is the data and MC single bin? bool fIsWriting; bool fSaveFine; /// OLD STUFF TO REMOVE TH1D* fMCHist_PDG[61]; ///< REMOVE OLD MC PDG Plot // Arrays for data entries Double_t* fXBins; ///< REMOVE xBin edges Double_t* fDataValues; ///< REMOVE data bin contents Double_t* fDataErrors; ///< REMOVE data bin errors Int_t fNDataPointsX; ///< REMOVE number of data points }; /*! @} */ #endif diff --git a/src/InputHandler/GiBUUNativeInputHandler.cxx b/src/InputHandler/GiBUUNativeInputHandler.cxx index f5c7335..3e9c832 100644 --- a/src/InputHandler/GiBUUNativeInputHandler.cxx +++ b/src/InputHandler/GiBUUNativeInputHandler.cxx @@ -1,511 +1,512 @@ #ifdef __GiBUU_ENABLED__ #include "GiBUUNativeInputHandler.h" #include "InputUtils.h" #include "PhysConst.h" bool GiBUUEventReader::SetBranchAddresses(TChain* tn){ tn->SetBranchAddress("weight", &weight); tn->SetBranchAddress("evType", &mode); tn->SetBranchAddress("process_ID", &process_ID); tn->SetBranchAddress("flavor_ID", &flavor_ID); tn->SetBranchAddress("numEnsembles", &num_ensembles); tn->SetBranchAddress("numRuns", &num_runs); tn->SetBranchAddress("nucleus_A", &nucleus_A); tn->SetBranchAddress("nucleus_Z", &nucleus_Z); tn->SetBranchAddress("barcode", &pdg, &b_pdg); tn->SetBranchAddress("Px", &Px, &b_Px); tn->SetBranchAddress("Py", &Py, &b_Py); tn->SetBranchAddress("Pz", &Pz, &b_Pz); tn->SetBranchAddress("E", &E, &b_E); tn->SetBranchAddress("x", &x, &b_x); tn->SetBranchAddress("y", &y, &b_y); tn->SetBranchAddress("z", &z, &b_z); tn->SetBranchAddress("lepIn_Px", &lepIn_Px); tn->SetBranchAddress("lepIn_Py", &lepIn_Py); tn->SetBranchAddress("lepIn_Pz", &lepIn_Pz); tn->SetBranchAddress("lepIn_E", &lepIn_E); tn->SetBranchAddress("lepOut_Px", &lepOut_Px); tn->SetBranchAddress("lepOut_Py", &lepOut_Py); tn->SetBranchAddress("lepOut_Pz", &lepOut_Pz); tn->SetBranchAddress("lepOut_E", &lepOut_E); tn->SetBranchAddress("nuc_Px", &nuc_Px); tn->SetBranchAddress("nuc_Py", &nuc_Py); tn->SetBranchAddress("nuc_Pz", &nuc_Pz); tn->SetBranchAddress("nuc_E", &nuc_E); tn->SetBranchAddress("nuc_charge", &nuc_chrg); return 0; } int ConvertModeGiBUUtoNEUT(int gibuu_mode, int process_ID, int struck_nucleon_pdg, int first_part_pdg){ bool IsCC = (abs(process_ID) == 2) ? true : false; bool IsNu = (process_ID > 0) ? true : false; switch(gibuu_mode){ // QE/elastic case 1: if (IsCC) { return (IsNu ? 1 : -1); } else { return (IsNu ? 1 : -1) * ((struck_nucleon_pdg == 2212) ? 51 : 52); } // Different resonances case 2: case 3: case 4: case 5: case 6: case 7: case 8: case 9: case 10: case 11: case 12: case 13: case 14: case 15: case 16: case 17: case 18: case 19: case 20: case 21: case 22: case 23: case 24: case 25: case 26: case 27: case 28: case 29: case 30: case 31: { if (IsCC){ if (IsNu){ if (struck_nucleon_pdg == 2212) return 11; if (first_part_pdg == 111) return 12; return 13; } else { if (struck_nucleon_pdg == 2112) return -11; if (first_part_pdg == 111) return -12; return -13; } } else { if (struck_nucleon_pdg == 2212) { if (first_part_pdg == 111) return 32 * (IsNu ? 1 : -1); return 34 * (IsNu ? 1 : -1); } else { if (first_part_pdg == 111) return 31 * (IsNu ? 1 : -1); return 33 * (IsNu ? 1 : -1); } } } case 32: case 33: { // 1Pi Bkg return (IsNu ? 1 : -1) * (10 + 20 * (!IsCC)); } case 34: { // DIS return (IsNu ? 1 : -1) * (26 + 20 * (!IsCC)); } case 35: case 36: { // MEC/2p-2h return (IsNu ? 1 : -1) * (2 + 40 * (!IsCC)); } case 37: { // MultiPi return (IsNu ? 1 : -1) * (21 + 20 * (!IsCC)); } default: NUIS_ERR(WRN, "Unable to map GiBUU code " << gibuu_mode << " to a NEUT mode"); } return 0; } int GetGiBUUNuPDG(int flavor_ID, int process_ID){ int pdg = 0; switch(flavor_ID){ case 1: pdg = 12; break; case 2: pdg = 14; break; case 3: pdg = 16; break; } if (process_ID < 0) pdg*=-1; return pdg; } int GetGiBUULepPDG(int flavor_ID, int process_ID){ int nuPDG = GetGiBUUNuPDG(flavor_ID, process_ID); // Check for EM if (abs(process_ID) == 1){ NUIS_ABORT("GiBUU file includes electron scattering events, not currently supported!"); } else if (abs(process_ID) == 3){ return nuPDG; } else if (abs(process_ID) == 2){ return (nuPDG > 0) ? nuPDG - 1 : nuPDG + 1; } NUIS_ABORT("Unknown GiBUU process_ID " << process_ID); } // Check whether the particle is on-shell or not -int CheckGiBUUParticleStatus(double E, int pdg, double dist){ +int CheckGiBUUParticleStatus(double E, int pdg, double dist, int targetA){ double onshell_mass = PhysConst::GetMass(pdg); // Hard code some other values used in GiBUU... if (pdg == 2212 || pdg == 2112) onshell_mass = 0.938; if (abs(pdg) == 211 || pdg == 111) onshell_mass = 0.138; // If the particle is still within the nucleus, it's not final state // Not that this depends on the nucleus, and too few time steps could cause an issue // 6 fm is Ulrich's guess for Argon (and will be fine for smaller nuclei) - if (dist < 6) { + // Note that hydrogen is an exception because of how the timesteps work + if (dist < 6 && targetA > 1) { return kFSIState; } // Check for unknown particles and default to on shell if (onshell_mass == -1) return kFinalState; // Ulrich's second criteria if (E < onshell_mass) return kFSIState; return kFinalState; } GiBUUNativeInputHandler::~GiBUUNativeInputHandler() { jointtype .clear(); jointrequested.clear(); } GiBUUNativeInputHandler::GiBUUNativeInputHandler(std::string const &handle, std::string const &rawinputs) { NUIS_LOG(SAM, "Creating GiBUUNativeInputHandler : " << handle); // Run a joint input handling fName = handle; fEventType = kGiBUU; fGiBUUTree = new TChain("RootTuple"); // Loop over all inputs and grab flux, eventhist, and nevents std::vector inputs = InputUtils::ParseInputFileList(rawinputs); for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) { // Open File for histogram access NUIS_LOG(SAM, "Opening event file " << inputs[inp_it]); TFile *inp_file = new TFile(inputs[inp_it].c_str(), "READ"); if ((!inp_file) || (!inp_file->IsOpen())) { NUIS_ABORT( "GiBUU file !IsOpen() at : '" << inputs[inp_it] << "'" << std::endl << "Check that your file paths are correct and the file exists!"); } // Get Flux/Event hist TH1D *fluxhist = (TH1D*)inp_file->Get(PlotUtils::GetObjectWithName(inp_file, "flux").c_str()); TH1D *eventhist = (TH1D*)inp_file->Get(PlotUtils::GetObjectWithName(inp_file, "evtrt").c_str()); if (!fluxhist || !eventhist) { NUIS_ERR(FTL, "Input File Contents: " << inputs[inp_it]); inp_file->ls(); NUIS_ABORT("GiBUU file doesn't contain flux/xsec info. You may have to " "use PrepareGiBUU!"); } // Get N Events TChain *tempChain = new TChain("RootTuple"); tempChain->Add(inputs[inp_it].c_str()); if (!tempChain){ NUIS_ERR(FTL, "RootTuple tree missing from GiBUU file " << inputs[inp_it]); NUIS_ABORT( "Check your inputs, they may need to be completely regenerated!"); } int nevents = tempChain->GetEntries(); if (nevents <= 0) { NUIS_ABORT("Trying to a TTree with " << nevents << " to TChain from : " << inputs[inp_it]); } // Get specific information about the config from the file GiBUUEventReader *tempReader = new GiBUUEventReader(); tempReader->SetBranchAddresses(tempChain); tempChain ->GetEntry(0); // Register input to form flux/event rate hists RegisterJointInput(inputs[inp_it], tempReader->process_ID, tempReader->flavor_ID, tempReader->nucleus_A, nevents, tempReader->nucleus_A*tempReader->num_runs*tempReader->num_ensembles, fluxhist, eventhist); // Add To TChain fGiBUUTree->AddFile(inputs[inp_it].c_str()); delete tempChain; delete tempReader; } // Registor all our file inputs SetupJointInputs(); // Create Fit Event fNUISANCEEvent = new FitEvent(); fGiReader = new GiBUUEventReader(); fGiReader->SetBranchAddresses(fGiBUUTree); fNUISANCEEvent->HardReset(); }; FitEvent *GiBUUNativeInputHandler::GetNuisanceEvent(const UInt_t ent, const bool lightweight) { UInt_t entry = ent + fSkip; // Check out of bounds if (entry >= (UInt_t)fNEvents) return NULL; fGiBUUTree->GetEntry(entry); fNUISANCEEvent->ResetEvent(); fNUISANCEEvent->fEventNo = entry; // Run NUISANCE Vector Filler if (!lightweight) { CalcNUISANCEKinematics(); } #ifdef __PROB3PP_ENABLED__ else { fNUISANCEEvent->probe_E = fGiReader->lepIn_E*1E3; fNUISANCEEvent->probe_pdg = GetGiBUUNuPDG(fGiReader->flavor_ID, fGiReader->process_ID); } #endif fNUISANCEEvent->InputWeight *= GetInputWeight(entry); return fNUISANCEEvent; } void GiBUUNativeInputHandler::CalcNUISANCEKinematics() { // Reset all variables // fNUISANCEEvent->ResetEvent(); FitEvent *evt = fNUISANCEEvent; evt->Mode = ConvertModeGiBUUtoNEUT(fGiReader->mode, fGiReader->process_ID, (fGiReader->nuc_chrg) ? 2212 : 2112, (*fGiReader->pdg)[0]); // evt->fEventNo = 0.0; evt->fTotCrs = 0; evt->fTargetA = fGiReader->nucleus_A; evt->fTargetZ = fGiReader->nucleus_Z; evt->fTargetH = 0; evt->fBound = 0.0; // Extra GiBUU Input Weight evt->InputWeight = fGiReader->weight; // Check number of particles in the stack. Note the additional 3 particles! uint npart = fGiReader->pdg->size(); int kmax = evt->kMaxParticles; if ((UInt_t)npart+3 > (UInt_t)kmax) { NUIS_ERR(WRN, "GiBUU has too many particles (" << npart << ") Expanding Stack."); fNUISANCEEvent->ExpandParticleStack(npart+2); } // Stuff the leptons in evt->fParticleState[0] = kInitialState; evt->fParticleMom[0][0] = fGiReader->lepIn_Px * 1E3; evt->fParticleMom[0][1] = fGiReader->lepIn_Py * 1E3; evt->fParticleMom[0][2] = fGiReader->lepIn_Pz * 1E3; evt->fParticleMom[0][3] = fGiReader->lepIn_E * 1E3; evt->fParticlePDG[0] = GetGiBUUNuPDG(fGiReader->flavor_ID, fGiReader->process_ID); evt->fParticleState[1] = kFinalState; evt->fParticleMom[1][0] = fGiReader->lepOut_Px * 1E3; evt->fParticleMom[1][1] = fGiReader->lepOut_Py * 1E3; evt->fParticleMom[1][2] = fGiReader->lepOut_Pz * 1E3; evt->fParticleMom[1][3] = fGiReader->lepOut_E * 1E3; evt->fParticlePDG[1] = GetGiBUULepPDG(fGiReader->flavor_ID, fGiReader->process_ID); // Also the struck nucleon (note only one!) evt->fParticleState[2] = kInitialState; evt->fParticleMom[2][0] = fGiReader->nuc_Px * 1E3; evt->fParticleMom[2][1] = fGiReader->nuc_Py * 1E3; evt->fParticleMom[2][2] = fGiReader->nuc_Pz * 1E3; evt->fParticleMom[2][3] = fGiReader->nuc_E * 1E3; evt->fParticlePDG[2] = (fGiReader->nuc_chrg) ? 2212 : 2112; // Create Stack evt->fNParticles = 3; for (uint i = 0; i < npart; i++) { int curpart = evt->fNParticles; double dist = sqrt((*fGiReader->x)[i]*(*fGiReader->x)[i] + (*fGiReader->y)[i]*(*fGiReader->y)[i] + (*fGiReader->z)[i]*(*fGiReader->z)[i]); // Set State - evt->fParticleState[curpart] = CheckGiBUUParticleStatus((*fGiReader->E)[i], (*fGiReader->pdg)[i], dist); + evt->fParticleState[curpart] = CheckGiBUUParticleStatus((*fGiReader->E)[i], (*fGiReader->pdg)[i], dist, evt->fTargetA); // Mom evt->fParticleMom[curpart][0] = (*fGiReader->Px)[i] * 1E3; evt->fParticleMom[curpart][1] = (*fGiReader->Py)[i] * 1E3; evt->fParticleMom[curpart][2] = (*fGiReader->Pz)[i] * 1E3; evt->fParticleMom[curpart][3] = (*fGiReader->E)[i] * 1E3; // PDG evt->fParticlePDG[curpart] = (*fGiReader->pdg)[i]; // Add to total particles evt->fNParticles++; } // Run Initial, FSI, Final, Other ordering. fNUISANCEEvent->OrderStack(); FitParticle *ISAnyLepton = fNUISANCEEvent->GetHMISAnyLeptons(); if (ISAnyLepton) { fNUISANCEEvent->probe_E = ISAnyLepton->E(); fNUISANCEEvent->probe_pdg = ISAnyLepton->PDG(); } return; } void GiBUUNativeInputHandler::Print() {} void GiBUUNativeInputHandler::RegisterJointInput(std::string input, int process_ID, int flavor_ID, int nnucleons, int n, int nrequested, TH1D *f, TH1D *e){ if (jointfluxinputs.size() == 0) { jointindexswitch = 0; fNEvents = 0; } // Push into individual input vectors jointfluxinputs.push_back((TH1D *)f->Clone()); jointeventinputs.push_back((TH1D *)e->Clone()); // This has to correspond to the real number of events jointindexlow .push_back(fNEvents); jointindexhigh.push_back(fNEvents + n); fNEvents += n; // This should be nnucleons*num_ensembles*num_runs (the number of requested events) jointrequested .push_back(nrequested); // To keep track of how many of each type there are jointtype .push_back(nnucleons*1000+process_ID*10 + flavor_ID); } void GiBUUNativeInputHandler::SetupJointInputs() { // Always apply the scaling jointinput = true; if (jointeventinputs.size() > 1) { jointindexswitch = 0; } fMaxEvents = FitPar::Config().GetParI("MAXEVENTS"); if (fMaxEvents != -1 and jointeventinputs.size() > 1) { NUIS_ABORT("Can only handle joint inputs when config MAXEVENTS = -1!"); } // Need to know what the total number of nucleons is for correct normalization... int total_unique_nucl = 0; std::vector unique_nucl; std::vector nsame_vect; for (size_t i = 0; i < jointeventinputs.size(); i++) { // Assume all of the fluxes are the same if (!fFluxHist) fFluxHist = (TH1D*)jointfluxinputs[i]->Clone(); // Get the list of unique nuclei nuclei int this_nucl = jointtype[i]/1000; if (std::find(unique_nucl.begin(), unique_nucl.end(), this_nucl) == unique_nucl.end()){ unique_nucl.push_back(this_nucl); total_unique_nucl += this_nucl; } // Get the total event rate and number of requested events TH1D *this_evt_total = NULL; int nsame = 0; // How many files were set up in the same way? for (size_t j = 0; j < jointeventinputs.size(); j++) { if (jointtype[j] != jointtype[i]) continue; nsame++; if (!this_evt_total) this_evt_total = (TH1D*)jointeventinputs[j]->Clone(); else this_evt_total->Add(jointeventinputs[j]); } nsame_vect.push_back(nsame); // Need to average events with the same run config, and add others this_evt_total->Scale(jointtype[i]/1000/double(nsame)); if (!fEventHist) fEventHist = (TH1D*)this_evt_total->Clone(); else fEventHist ->Add(this_evt_total); delete this_evt_total; } // Get the correct / nucleon event rate fEventHist->Scale(1/double(total_unique_nucl)); // Loop over the samples to get the event scaling correct for (size_t i = 0; i < jointeventinputs.size(); i++) { // This simply reverses the usual scaling to get to a flux averaged XSEC used everywhere else in NUISANCE double scale = fFluxHist->Integral("width")*fNEvents/fEventHist->Integral("width"); // Need to correctly weight by the number of nucleons to arrive at the /nucleon XSEC for a compound target scale *= jointtype[i]/1000/double(total_unique_nucl)/nsame_vect[i]; jointindexscale.push_back(scale); } fEventHist->SetNameTitle((fName + "_EVT").c_str(), (fName + "_EVT").c_str()); fFluxHist->SetNameTitle((fName + "_FLUX").c_str(), (fName + "_FLUX").c_str()); // Setup Max Events if (fMaxEvents > 1 && fMaxEvents < fNEvents) { if (LOG_LEVEL(SAM)) { std::cout << "\t\t|-> Read Max Entries : " << fMaxEvents << std::endl; } fNEvents = fMaxEvents; } // Print out Status if (LOG_LEVEL(SAM)) { std::cout << "\t\t|-> Total Entries : " << fNEvents << std::endl << "\t\t|-> Event Integral : " << fEventHist->Integral("width") * 1.E-38 << " events/nucleon" << std::endl << "\t\t|-> Flux Integral : " << fFluxHist->Integral("width") << " /cm2" << std::endl << "\t\t|-> Event/Flux : " << fEventHist->Integral("width") * 1.E-38 / fFluxHist->Integral("width") << " cm2/nucleon" << std::endl; } } #endif diff --git a/src/InputHandler/GiBUUNativeInputHandler.h b/src/InputHandler/GiBUUNativeInputHandler.h index ad384eb..d6cd594 100644 --- a/src/InputHandler/GiBUUNativeInputHandler.h +++ b/src/InputHandler/GiBUUNativeInputHandler.h @@ -1,110 +1,110 @@ #ifndef GIBUUNATIVEINPUTHANDLER_H #define GIBUUNATIVEINPUTHANDLER_H #ifdef __GiBUU_ENABLED__ /*! * \addtogroup InputHandler * @{ */ #include "InputHandler.h" #include "PlotUtils.h" struct GiBUUEventReader { public: // Event weight double weight; // The full particle stack TBranch *b_pdg = 0; TBranch *b_Px = 0; TBranch *b_Py = 0; TBranch *b_Pz = 0; TBranch *b_E = 0; // These are only if you build with the right options TBranch *b_x = 0; TBranch *b_y = 0; TBranch *b_z = 0; std::vector *pdg = 0; std::vector *Px = 0; std::vector *Py = 0; std::vector *Pz = 0; std::vector *E = 0; std::vector *x = 0; std::vector *y = 0; std::vector *z = 0; int mode; int process_ID; int flavor_ID; int num_ensembles; int num_runs; int nucleus_A; int nucleus_Z; double lepIn_E; double lepIn_Px; double lepIn_Py; double lepIn_Pz; double lepOut_E; double lepOut_Px; double lepOut_Py; double lepOut_Pz; double nuc_E; double nuc_Px; double nuc_Py; double nuc_Pz; int nuc_chrg; bool SetBranchAddresses(TChain* tn); }; int GetGiBUUNuPDG(int flavor_ID, int process_ID); int GetGiBUULepPDG(int flavor_ID, int process_ID); -int CheckGiBUUParticleStatus(double E, int pdg, double dist); +int CheckGiBUUParticleStatus(double E, int pdg, double dist, int targetA); int ConvertModeGiBUUtoNEUT(int gibuu_mode, int process_ID, int struck_nucleon_pdg, int first_part_pdg); /// GiBUU Handler to read in GiBUU's ROOT format class GiBUUNativeInputHandler : public InputHandlerBase { public: /// Standard constructor given name and inputs GiBUUNativeInputHandler(std::string const& handle, std::string const& rawinputs); /// Actually need to clean up some GiBUU-specific things ~GiBUUNativeInputHandler(); /// Returns NUISANCE Format Event from the GiBUU stack FitEvent* GetNuisanceEvent(const UInt_t entry, const bool lightweight = false); /// Override to handle the fact that GiBUU already averages over nevents. void SetupJointInputs(); /// Special info is required for GiBUU... void RegisterJointInput(std::string input, int process_ID, int flavor_ID, int nnucleons, int n, int nrequested, TH1D *f, TH1D *e); /// Fill NUISANCE Particle Stack void CalcNUISANCEKinematics(); /// Print event information void Print(); private: GiBUUEventReader* fGiReader; TChain* fGiBUUTree; ///< GiBUU Event Tree // For GiBUU specific normalization std::vector jointrequested; std::vector jointtype; }; #endif #endif diff --git a/src/InputHandler/NEUTInputHandler.cxx b/src/InputHandler/NEUTInputHandler.cxx index 58cb9b6..337b153 100644 --- a/src/InputHandler/NEUTInputHandler.cxx +++ b/src/InputHandler/NEUTInputHandler.cxx @@ -1,556 +1,557 @@ #if defined(__NEUT_ENABLED__) || defined(NEUT_EVENT_ENABLED) #include "NEUTInputHandler.h" #include "InputUtils.h" #include "PlotUtils.h" #include "TTreePerfStats.h" #include "fsihistC.h" #include "necardC.h" #include "nefillverC.h" #include "neutcrsC.h" #include "neutfsipart.h" #include "neutfsivert.h" #include "neutmodelC.h" #include "neutparamsC.h" #include "neutpart.h" #include "neutrootTreeSingleton.h" #include "neutvect.h" #include "neworkC.h" #include "vcworkC.h" #include "posinnucC.h" #ifdef __NEUT_NUCFSI_ENABLED__ #include "neutnucfsistep.h" #include "neutnucfsivert.h" #include "nucleonfsihistC.h" #endif NEUTGeneratorInfo::~NEUTGeneratorInfo() { DeallocateParticleStack(); } void NEUTGeneratorInfo::AddBranchesToTree(TTree *tn) { tn->Branch("NEUTParticleN", fNEUTParticleN, "NEUTParticleN/I"); tn->Branch("NEUTParticleStatusCode", fNEUTParticleStatusCode, "NEUTParticleStatusCode[NEUTParticleN]/I"); tn->Branch("NEUTParticleAliveCode", fNEUTParticleAliveCode, "NEUTParticleAliveCode[NEUTParticleN]/I"); } void NEUTGeneratorInfo::SetBranchesFromTree(TTree *tn) { tn->SetBranchAddress("NEUTParticleN", &fNEUTParticleN); tn->SetBranchAddress("NEUTParticleStatusCode", &fNEUTParticleStatusCode); tn->SetBranchAddress("NEUTParticleAliveCode", &fNEUTParticleAliveCode); } void NEUTGeneratorInfo::AllocateParticleStack(int stacksize) { fNEUTParticleN = 0; fNEUTParticleStatusCode = new int[stacksize]; fNEUTParticleStatusCode = new int[stacksize]; } void NEUTGeneratorInfo::DeallocateParticleStack() { delete fNEUTParticleStatusCode; delete fNEUTParticleAliveCode; } void NEUTGeneratorInfo::FillGeneratorInfo(NeutVect *nevent) { Reset(); for (int i = 0; i < nevent->Npart(); i++) { fNEUTParticleStatusCode[i] = nevent->PartInfo(i)->fStatus; fNEUTParticleAliveCode[i] = nevent->PartInfo(i)->fIsAlive; fNEUTParticleN++; } } void NEUTGeneratorInfo::Reset() { for (int i = 0; i < fNEUTParticleN; i++) { fNEUTParticleStatusCode[i] = -1; fNEUTParticleAliveCode[i] = 9; } fNEUTParticleN = 0; } NEUTInputHandler::NEUTInputHandler(std::string const &handle, std::string const &rawinputs) { NUIS_LOG(SAM, "Creating NEUTInputHandler : " << handle); // Run a joint input handling fName = handle; // Setup the TChain fNEUTTree = new TChain("neuttree"); fSaveExtra = FitPar::Config().GetParB("SaveExtraNEUT"); fCacheSize = FitPar::Config().GetParI("CacheSize"); fMaxEvents = FitPar::Config().GetParI("MAXEVENTS"); // Loop over all inputs and grab flux, eventhist, and nevents std::vector inputs = InputUtils::ParseInputFileList(rawinputs); for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) { // Open File for histogram access TFile *inp_file = new TFile(inputs[inp_it].c_str(), "READ"); if (!inp_file or inp_file->IsZombie()) { NUIS_ABORT( "NEUT File IsZombie() at : '" << inputs[inp_it] << "'" << std::endl << "Check that your file paths are correct and the file exists!" << std::endl << "$ ls -lh " << inputs[inp_it]); } // Get Flux/Event hist TH1D *fluxhist = (TH1D *)inp_file->Get( (PlotUtils::GetObjectWithName(inp_file, "flux")).c_str()); TH1D *eventhist = (TH1D *)inp_file->Get( (PlotUtils::GetObjectWithName(inp_file, "evt")).c_str()); if (!fluxhist or !eventhist) { NUIS_ERR(FTL, "Input File Contents: " << inputs[inp_it]); inp_file->ls(); NUIS_ABORT("NEUT FILE doesn't contain flux/xsec info. You may have to " "regenerate your MC!"); } // Get N Events TTree *neuttree = (TTree *)inp_file->Get("neuttree"); if (!neuttree) { NUIS_ERR(FTL, "neuttree not located in NEUT file: " << inputs[inp_it]); NUIS_ABORT( "Check your inputs, they may need to be completely regenerated!"); throw; } int nevents = neuttree->GetEntries(); if (nevents <= 0) { NUIS_ABORT("Trying to a TTree with " << nevents << " to TChain from : " << inputs[inp_it]); } // Register input to form flux/event rate hists RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist); // Add To TChain fNEUTTree->AddFile(inputs[inp_it].c_str()); } // Registor all our file inputs SetupJointInputs(); // Assign to tree fEventType = kNEUT; fNeutVect = NULL; fNEUTTree->SetBranchAddress("vectorbranch", &fNeutVect); #if defined(ROOT6) && defined(__NEUT_VERSION__) && (__NEUT_VERSION__ >= 541) fNEUTTree->SetAutoDelete(true); #endif fNEUTTree->GetEntry(0); // Create Fit Event fNUISANCEEvent = new FitEvent(); fNUISANCEEvent->SetNeutVect(fNeutVect); if (fSaveExtra) { fNeutInfo = new NEUTGeneratorInfo(); fNUISANCEEvent->AddGeneratorInfo(fNeutInfo); } fNUISANCEEvent->HardReset(); }; NEUTInputHandler::~NEUTInputHandler(){ // if (fNEUTTree) delete fNEUTTree; // if (fNeutVect) delete fNeutVect; // if (fNeutInfo) delete fNeutInfo; }; void NEUTInputHandler::CreateCache() { if (fCacheSize > 0) { // fNEUTTree->SetCacheEntryRange(0, fNEvents); fNEUTTree->AddBranchToCache("vectorbranch", 1); fNEUTTree->SetCacheSize(fCacheSize); } } void NEUTInputHandler::RemoveCache() { // fNEUTTree->SetCacheEntryRange(0, fNEvents); fNEUTTree->AddBranchToCache("vectorbranch", 0); fNEUTTree->SetCacheSize(0); } FitEvent *NEUTInputHandler::GetNuisanceEvent(const UInt_t ent, const bool lightweight) { UInt_t entry = ent + fSkip; // Catch too large entries if (entry >= (UInt_t)fNEvents) return NULL; // Read Entry from TTree to fill NEUT Vect in BaseFitEvt; fNEUTTree->GetEntry(entry); // Run NUISANCE Vector Filler if (!lightweight) { CalcNUISANCEKinematics(); } #ifdef __PROB3PP_ENABLED__ else { UInt_t npart = fNeutVect->Npart(); for (size_t i = 0; i < npart; i++) { NeutPart *part = fNUISANCEEvent->fNeutVect->PartInfo(i); if ((part->fIsAlive == false) && (part->fStatus == -1) && std::count(PhysConst::pdg_neutrinos, PhysConst::pdg_neutrinos + 4, part->fPID)) { fNUISANCEEvent->probe_E = part->fP.T(); fNUISANCEEvent->probe_pdg = part->fPID; break; } else { continue; } } } #endif // Setup Input scaling for joint inputs fNUISANCEEvent->InputWeight = GetInputWeight(entry); // Return event pointer return fNUISANCEEvent; } // From NEUT neutclass/neutpart.h // Bool_t fIsAlive; // Particle should be tracked or not // ( in the detector simulator ) // // Int_t fStatus; // Status flag of this particle // -2: Non existing particle // -1: Initial state particle // 0: Normal // 1: Decayed to the other particle // 2: Escaped from the detector // 3: Absorped // 4: Charge exchanged // 5: Pauli blocked // 6: N/A // 7: Produced child particles // 8: Inelastically scattered // int NEUTInputHandler::GetNeutParticleStatus(NeutPart *part) { // State int state = kUndefinedState; // Remove Pauli blocked events, probably just single pion events if (part->fStatus == 5) { state = kFSIState; // fStatus == -1 means initial state } else if (part->fIsAlive == false && part->fStatus == -1) { state = kInitialState; // NEUT has a bit of a strange convention for fIsAlive and fStatus // combinations // for NC and neutrino particle isAlive true/false and status 2 means // final state particle // for other particles in NC status 2 means it's an FSI particle // for CC it means it was an FSI particle } else if (part->fStatus == 2) { // NC case is a little strange... The outgoing neutrino might be alive or // not alive. Remaining particles with status 2 are FSI particles that // reinteracted if (abs(fNeutVect->Mode) > 30 && (abs(part->fPID) == 16 || abs(part->fPID) == 14 || abs(part->fPID) == 12)) { state = kFinalState; // The usual CC case } else if (part->fIsAlive == true) { state = kFSIState; } } else if (part->fIsAlive == true && part->fStatus == 2 && (abs(part->fPID) == 16 || abs(part->fPID) == 14 || abs(part->fPID) == 12)) { state = kFinalState; } else if (part->fIsAlive == true && part->fStatus == 0) { state = kFinalState; } else if (!part->fIsAlive && (part->fStatus == 1 || part->fStatus == 3 || part->fStatus == 4 || part->fStatus == 7 || part->fStatus == 8)) { state = kFSIState; // There's one hyper weird case where fStatus = -3. This apparently // corresponds to a nucleon being ejected via pion FSI when there is "data // available" } else if (!part->fIsAlive && (part->fStatus == -3)) { state = kUndefinedState; // NC neutrino outgoing } else if (!part->fIsAlive && part->fStatus == 0 && (abs(part->fPID) == 16 || abs(part->fPID) == 14 || abs(part->fPID) == 12)) { state = kFinalState; // Warn if we still find alive particles without classifying them } else if (part->fIsAlive == true) { NUIS_ABORT("Undefined NEUT state " << " Alive: " << part->fIsAlive << " Status: " << part->fStatus << " PDG: " << part->fPID << " Mode: " << fNeutVect->Mode); } else if (abs(fNeutVect->Mode) == 35){ NUIS_ERR(WRN, "Marking nonsensical CC difractive event as undefined " << " Alive: " << part->fIsAlive << " Status: " << part->fStatus << " PDG: " << part->fPID << " Mode: " << fNeutVect->Mode); state = kUndefinedState; // Warn if we find dead particles that we haven't classified } else { NUIS_ABORT("Undefined NEUT state " << " Alive: " << part->fIsAlive << " Status: " << part->fStatus << " PDG: " << part->fPID << " Mode: " << fNeutVect->Mode); } return state; } void NEUTInputHandler::CalcNUISANCEKinematics() { // Reset all variables fNUISANCEEvent->ResetEvent(); // Fill Globals fNUISANCEEvent->Mode = fNeutVect->Mode; fNUISANCEEvent->fEventNo = fNeutVect->EventNo; fNUISANCEEvent->fTargetA = fNeutVect->TargetA; fNUISANCEEvent->fTargetZ = fNeutVect->TargetZ; fNUISANCEEvent->fTargetH = fNeutVect->TargetH; fNUISANCEEvent->fBound = bool(fNeutVect->Ibound); - if (fNUISANCEEvent->fBound) { + if (fNUISANCEEvent->fBound || + (!fNUISANCEEvent->fBound && abs(fNUISANCEEvent->Mode) == 16)) { // Make special exception for coherent events (mode 16) fNUISANCEEvent->fTargetPDG = TargetUtils::GetTargetPDGFromZA( fNUISANCEEvent->fTargetZ, fNUISANCEEvent->fTargetA); } else { fNUISANCEEvent->fTargetPDG = 1000010010; } // Check Particle Stack UInt_t npart = fNeutVect->Npart(); UInt_t kmax = fNUISANCEEvent->kMaxParticles; if (npart > kmax) { NUIS_ERR(WRN, "NEUT has too many particles. Expanding stack."); fNUISANCEEvent->ExpandParticleStack(npart); } UInt_t nprimary = fNeutVect->Nprimary(); // Fill Particle Stack for (size_t i = 0; i < npart; i++) { // Get Current Count int curpart = fNUISANCEEvent->fNParticles; // Get NEUT Particle NeutPart *part = fNeutVect->PartInfo(i); // State int state = GetNeutParticleStatus(part); // Remove Undefined if (kRemoveUndefParticles && state == kUndefinedState) continue; // Remove FSI if (kRemoveFSIParticles && state == kFSIState) continue; // Remove Nuclear if (kRemoveNuclearParticles && (state == kNuclearInitial || state == kNuclearRemnant)) continue; // State fNUISANCEEvent->fParticleState[curpart] = state; // Is the paricle associated with the primary vertex? bool primary = false; // NEUT events are just popped onto the stack as primary, then continues to // be non-primary if (i < nprimary) primary = true; fNUISANCEEvent->fPrimaryVertex[curpart] = primary; // Mom fNUISANCEEvent->fParticleMom[curpart][0] = part->fP.X(); fNUISANCEEvent->fParticleMom[curpart][1] = part->fP.Y(); fNUISANCEEvent->fParticleMom[curpart][2] = part->fP.Z(); fNUISANCEEvent->fParticleMom[curpart][3] = part->fP.T(); // PDG fNUISANCEEvent->fParticlePDG[curpart] = part->fPID; // Add up particle count fNUISANCEEvent->fNParticles++; } // Save Extra Generator Info if (fSaveExtra) { fNeutInfo->FillGeneratorInfo(fNeutVect); } // Run Initial, FSI, Final, Other ordering. fNUISANCEEvent->OrderStack(); FitParticle *ISAnyLepton = fNUISANCEEvent->GetHMISAnyLeptons(); if (ISAnyLepton) { fNUISANCEEvent->probe_E = ISAnyLepton->E(); fNUISANCEEvent->probe_pdg = ISAnyLepton->PDG(); } return; } #endif #ifdef NEED_FILL_NEUT_COMMONS void NEUTUtils::FillNeutCommons(NeutVect *nvect) { // WARNING: This has only been implemented for a neuttree and not GENIE // This should be kept in sync with T2KNIWGUtils::GetNIWGEvent(TTree) // NEUT version info. Can't get it to compile properly with this yet // neutversion_.corev = nvect->COREVer; // neutversion_.nucev = nvect->NUCEVer; // neutversion_.nuccv = nvect->NUCCVer; // Documentation: See nework.h nework_.modene = nvect->Mode; nework_.numne = nvect->Npart(); #ifdef NEUT_COMMON_QEAV nemdls_.mdlqeaf = nvect->QEAVForm; #else nemdls_.mdlqeaf = nvect->QEVForm; #endif nemdls_.mdlqe = nvect->QEModel; nemdls_.mdlspi = nvect->SPIModel; nemdls_.mdldis = nvect->DISModel; nemdls_.mdlcoh = nvect->COHModel; neutcoh_.necohepi = nvect->COHModel; nemdls_.xmaqe = nvect->QEMA; nemdls_.xmvqe = nvect->QEMV; nemdls_.kapp = nvect->KAPPA; // nemdls_.sccfv = SCCFVdef; // nemdls_.sccfa = SCCFAdef; // nemdls_.fpqe = FPQEdef; nemdls_.xmaspi = nvect->SPIMA; nemdls_.xmvspi = nvect->SPIMV; nemdls_.xmares = nvect->RESMA; nemdls_.xmvres = nvect->RESMV; neut1pi_.xmanffres = nvect->SPIMA; neut1pi_.xmvnffres = nvect->SPIMV; neut1pi_.xmarsres = nvect->RESMA; neut1pi_.xmvrsres = nvect->RESMV; neut1pi_.neiff = nvect->SPIForm; neut1pi_.nenrtype = nvect->SPINRType; neut1pi_.rneca5i = nvect->SPICA5I; neut1pi_.rnebgscl = nvect->SPIBGScale; nemdls_.xmacoh = nvect->COHMA; nemdls_.rad0nu = nvect->COHR0; // nemdls_.fa1coh = nvect->COHA1err; // nemdls_.fb1coh = nvect->COHb1err; // neutdis_.nepdf = NEPDFdef; // neutdis_.nebodek = NEBODEKdef; neutcard_.nefrmflg = nvect->FrmFlg; neutcard_.nepauflg = nvect->PauFlg; neutcard_.nenefo16 = nvect->NefO16; neutcard_.nemodflg = nvect->ModFlg; // neutcard_.nenefmodl = 1; // neutcard_.nenefmodh = 1; // neutcard_.nenefkinh = 1; // neutpiabs_.neabspiemit = 1; nenupr_.iformlen = nvect->FormLen; neutpiless_.ipilessdcy = nvect->IPilessDcy; neutpiless_.rpilessdcy = nvect->RPilessDcy; neutpiless_.ipilessdcy = nvect->IPilessDcy; neutpiless_.rpilessdcy = nvect->RPilessDcy; neffpr_.fefqe = nvect->NuceffFactorPIQE; neffpr_.fefqeh = nvect->NuceffFactorPIQEH; neffpr_.fefinel = nvect->NuceffFactorPIInel; neffpr_.fefabs = nvect->NuceffFactorPIAbs; neffpr_.fefcx = nvect->NuceffFactorPICX; neffpr_.fefcxh = nvect->NuceffFactorPICXH; neffpr_.fefcoh = nvect->NuceffFactorPICoh; neffpr_.fefqehf = nvect->NuceffFactorPIQEHKin; neffpr_.fefcxhf = nvect->NuceffFactorPICXKin; neffpr_.fefcohf = nvect->NuceffFactorPIQELKin; for (int i = 0; i < nework_.numne; i++) { nework_.ipne[i] = nvect->PartInfo(i)->fPID; nework_.pne[i][0] = (float)nvect->PartInfo(i)->fP.X() / 1000; // VC(NE)WORK in M(G)eV nework_.pne[i][1] = (float)nvect->PartInfo(i)->fP.Y() / 1000; // VC(NE)WORK in M(G)eV nework_.pne[i][2] = (float)nvect->PartInfo(i)->fP.Z() / 1000; // VC(NE)WORK in M(G)eV } // fsihist.h // neutroot fills a dummy object for events with no FSI to prevent memory leak // when // reading the TTree, so check for it here if ((int)nvect->NfsiVert() == 1) { // An event with FSI must have at least two vertices // if (nvect->NfsiPart()!=1 || nvect->Fsiprob!=-1) // ERR(WRN) << "T2KNeutUtils::fill_neut_commons(TTree) NfsiPart!=1 or // Fsiprob!=-1 when NfsiVert==1" << std::endl; fsihist_.nvert = 0; fsihist_.nvcvert = 0; fsihist_.fsiprob = 1; } else { // Real FSI event fsihist_.nvert = (int)nvect->NfsiVert(); for (int ivert = 0; ivert < fsihist_.nvert; ivert++) { fsihist_.iflgvert[ivert] = nvect->FsiVertInfo(ivert)->fVertID; fsihist_.posvert[ivert][0] = (float)nvect->FsiVertInfo(ivert)->fPos.X(); fsihist_.posvert[ivert][1] = (float)nvect->FsiVertInfo(ivert)->fPos.Y(); fsihist_.posvert[ivert][2] = (float)nvect->FsiVertInfo(ivert)->fPos.Z(); } fsihist_.nvcvert = nvect->NfsiPart(); for (int ip = 0; ip < fsihist_.nvcvert; ip++) { fsihist_.abspvert[ip] = (float)nvect->FsiPartInfo(ip)->fMomLab; fsihist_.abstpvert[ip] = (float)nvect->FsiPartInfo(ip)->fMomNuc; fsihist_.ipvert[ip] = nvect->FsiPartInfo(ip)->fPID; fsihist_.iverti[ip] = nvect->FsiPartInfo(ip)->fVertStart; fsihist_.ivertf[ip] = nvect->FsiPartInfo(ip)->fVertEnd; fsihist_.dirvert[ip][0] = (float)nvect->FsiPartInfo(ip)->fDir.X(); fsihist_.dirvert[ip][1] = (float)nvect->FsiPartInfo(ip)->fDir.Y(); fsihist_.dirvert[ip][2] = (float)nvect->FsiPartInfo(ip)->fDir.Z(); } fsihist_.fsiprob = nvect->Fsiprob; } neutcrscom_.crsx = nvect->Crsx; neutcrscom_.crsy = nvect->Crsy; neutcrscom_.crsz = nvect->Crsz; neutcrscom_.crsphi = nvect->Crsphi; neutcrscom_.crsq2 = nvect->Crsq2; neuttarget_.numbndn = nvect->TargetA - nvect->TargetZ; neuttarget_.numbndp = nvect->TargetZ; neuttarget_.numfrep = nvect->TargetH; neuttarget_.numatom = nvect->TargetA; posinnuc_.ibound = nvect->Ibound; // put empty nucleon FSI history (since it is not saved in the NeutVect // format) // Comment out as NEUT does not have the necessary proton FSI information yet // nucleonfsihist_.nfnvert = 0; // nucleonfsihist_.nfnstep = 0; } #endif diff --git a/src/Reweight/GENIEWeightEngine.cxx b/src/Reweight/GENIEWeightEngine.cxx index c248077..30aa79c 100644 --- a/src/Reweight/GENIEWeightEngine.cxx +++ b/src/Reweight/GENIEWeightEngine.cxx @@ -1,495 +1,502 @@ #include "GENIEWeightEngine.h" #ifdef __GENIE_ENABLED__ #ifdef GENIE_PRE_R3 #include "Algorithm/AlgConfigPool.h" #include "EVGCore/EventRecord.h" #include "GHEP/GHepRecord.h" #include "Ntuple/NtpMCEventRecord.h" #ifndef __NO_REWEIGHT__ #include "ReWeight/GReWeightAGKY.h" #include "ReWeight/GReWeightDISNuclMod.h" #include "ReWeight/GReWeightFGM.h" #include "ReWeight/GReWeightFZone.h" #include "ReWeight/GReWeightINuke.h" #include "ReWeight/GReWeightNonResonanceBkg.h" #include "ReWeight/GReWeightNuXSecCCQE.h" #include "ReWeight/GReWeightNuXSecCCQEvec.h" #include "ReWeight/GReWeightNuXSecCCRES.h" #include "ReWeight/GReWeightNuXSecCOH.h" #include "ReWeight/GReWeightNuXSecDIS.h" #include "ReWeight/GReWeightNuXSecNC.h" #include "ReWeight/GReWeightNuXSecNCEL.h" #include "ReWeight/GReWeightNuXSecNCRES.h" #include "ReWeight/GReWeightResonanceDecay.h" #include "ReWeight/GSystUncertainty.h" #ifdef __GENIE_EMP_MECRW_ENABLED #include "ReWeight/GReWeightXSecEmpiricalMEC.h" #endif #endif #if __GENIE_VERSION__ >= 212 #include "ReWeight/GReWeightNuXSecCCQEaxial.h" #endif #else // GENIE v3 #include "Framework/Algorithm/AlgConfigPool.h" #include "Framework/EventGen/EventRecord.h" #include "Framework/GHEP/GHepParticle.h" #include "Framework/GHEP/GHepRecord.h" #include "Framework/Ntuple/NtpMCEventRecord.h" #include "Framework/Utils/AppInit.h" #include "Framework/Utils/RunOpt.h" using namespace genie; #ifndef __NO_REWEIGHT__ #include "RwCalculators/GReWeightAGKY.h" #include "RwCalculators/GReWeightDISNuclMod.h" #include "RwCalculators/GReWeightFGM.h" #include "RwCalculators/GReWeightFZone.h" #include "RwCalculators/GReWeightINuke.h" #include "RwCalculators/GReWeightNonResonanceBkg.h" #include "RwCalculators/GReWeightNuXSecCCQE.h" #include "RwCalculators/GReWeightNuXSecCCQEaxial.h" #include "RwCalculators/GReWeightNuXSecCCQEvec.h" #include "RwCalculators/GReWeightNuXSecCCRES.h" #include "RwCalculators/GReWeightNuXSecCOH.h" #include "RwCalculators/GReWeightNuXSecDIS.h" #include "RwCalculators/GReWeightNuXSecNC.h" #include "RwCalculators/GReWeightNuXSecNCEL.h" #include "RwCalculators/GReWeightNuXSecNCRES.h" #ifdef USE_GENIE_XSECMEC #include "RwCalculators/GReWeightXSecMEC.h" #endif #include "RwCalculators/GReWeightResonanceDecay.h" #include "RwCalculators/GReWeightXSecEmpiricalMEC.h" #include "RwFramework/GSystUncertainty.h" using namespace genie::rew; #endif #endif #endif #include "FitLogger.h" GENIEWeightEngine::GENIEWeightEngine(std::string name) { #ifdef __GENIE_ENABLED__ #ifndef __NO_REWEIGHT__ // Setup the NEUT Reweight engien fCalcName = name; NUIS_LOG(DEB, "Setting up GENIE RW : " << fCalcName); // Create RW Engine suppressing cout StopTalking(); fGenieRW = new genie::rew::GReWeight(); // Get List of Vetos (Just for debugging) std::string rw_engine_list = FitPar::Config().GetParS("FitWeight_fGenieRW_veto"); bool xsec_ncel = rw_engine_list.find("xsec_ncel") == std::string::npos; bool xsec_ccqe = rw_engine_list.find("xsec_ccqe") == std::string::npos; bool xsec_coh = rw_engine_list.find("xsec_coh") == std::string::npos; bool xsec_nnres = rw_engine_list.find("xsec_nonresbkg") == std::string::npos; bool xsec_nudis = rw_engine_list.find("nuclear_dis") == std::string::npos; bool xsec_resdec = rw_engine_list.find("hadro_res_decay") == std::string::npos; bool xsec_fzone = rw_engine_list.find("hadro_intranuke") == std::string::npos; bool xsec_intra = rw_engine_list.find("hadro_fzone") == std::string::npos; bool xsec_agky = rw_engine_list.find("hadro_agky") == std::string::npos; bool xsec_qevec = rw_engine_list.find("xsec_ccqe_vec") == std::string::npos; bool xsec_dis = rw_engine_list.find("xsec_dis") == std::string::npos; bool xsec_nc = rw_engine_list.find("xsec_nc") == std::string::npos; bool xsec_ccres = rw_engine_list.find("xsec_ccres") == std::string::npos; bool xsec_ncres = rw_engine_list.find("xsec_ncres") == std::string::npos; bool xsec_nucqe = rw_engine_list.find("nuclear_qe") == std::string::npos; bool xsec_qeaxial = rw_engine_list.find("xsec_ccqe_axial") == std::string::npos; #ifdef __GENIE_EMP_MECRW_ENABLED bool xsec_empMEC = rw_engine_list.find("xsec_empMEC") == std::string::npos; #endif #ifdef USE_GENIE_XSECMEC bool xsec_MEC = rw_engine_list.find("xsec_MEC") == std::string::npos; #endif #ifndef GENIE_PRE_R3 genie::RunOpt *grunopt = genie::RunOpt::Instance(); grunopt->EnableBareXSecPreCalc(true); grunopt->SetEventGeneratorList(Config::GetParS("GENIEEventGeneratorList")); + if (!Config::HasPar("GENIETune")) { + StartTalking(); NUIS_ABORT( "GENIE tune was not specified, this is required when reweighting GENIE " "V3+ events. Add a config parameter like: to your nuisance card."); } grunopt->SetTuneName(Config::GetParS("GENIETune")); grunopt->BuildTune(); std::string genv = std::string(getenv("GENIE")) + "/config/Messenger_laconic.xml"; genie::utils::app_init::MesgThresholds(genv); #endif // Now actually add the RW Calcs if (xsec_ncel) fGenieRW->AdoptWghtCalc("xsec_ncel", new genie::rew::GReWeightNuXSecNCEL); if (xsec_ccqe) { fGenieRW->AdoptWghtCalc("xsec_ccqe", new genie::rew::GReWeightNuXSecCCQE); // (dynamic_cast (fGenieRW->WghtCalc("xsec_ccqe"))) // ->SetXSecModel( FitPar::Config().GetParS("GENIEXSecModelCCQE") ); } #ifdef __GENIE_EMP_MECRW_ENABLED if (xsec_empMEC) { fGenieRW->AdoptWghtCalc("xsec_empMEC", new genie::rew::GReWeightXSecEmpiricalMEC); } #endif #ifdef USE_GENIE_XSECMEC if (xsec_MEC) { fGenieRW->AdoptWghtCalc("xsec_MEC", new genie::rew::GReWeightXSecMEC); } #endif if (xsec_coh) { fGenieRW->AdoptWghtCalc("xsec_coh", new genie::rew::GReWeightNuXSecCOH()); // (dynamic_cast (fGenieRW->WghtCalc("xsec_coh"))) // ->SetXSecModel( FitPar::Config().GetParS("GENIEXSecModelCOH") ); } if (xsec_nnres) fGenieRW->AdoptWghtCalc("xsec_nonresbkg", new genie::rew::GReWeightNonResonanceBkg); if (xsec_nudis) fGenieRW->AdoptWghtCalc("nuclear_dis", new genie::rew::GReWeightDISNuclMod); if (xsec_resdec) fGenieRW->AdoptWghtCalc("hadro_res_decay", new genie::rew::GReWeightResonanceDecay); if (xsec_fzone) fGenieRW->AdoptWghtCalc("hadro_fzone", new genie::rew::GReWeightFZone); if (xsec_intra) fGenieRW->AdoptWghtCalc("hadro_intranuke", new genie::rew::GReWeightINuke); if (xsec_agky) fGenieRW->AdoptWghtCalc("hadro_agky", new genie::rew::GReWeightAGKY); if (xsec_qevec) fGenieRW->AdoptWghtCalc("xsec_ccqe_vec", new genie::rew::GReWeightNuXSecCCQEvec); #if __GENIE_VERSION__ >= 212 if (xsec_qeaxial) fGenieRW->AdoptWghtCalc("xsec_ccqe_axial", new genie::rew::GReWeightNuXSecCCQEaxial); #endif if (xsec_dis) fGenieRW->AdoptWghtCalc("xsec_dis", new genie::rew::GReWeightNuXSecDIS); if (xsec_nc) fGenieRW->AdoptWghtCalc("xsec_nc", new genie::rew::GReWeightNuXSecNC); if (xsec_ccres) { #if __GENIE_VERSION__ < 213 fGenieRW->AdoptWghtCalc("xsec_ccres", new genie::rew::GReWeightNuXSecCCRES); #else fGenieRW->AdoptWghtCalc( "xsec_ccres", new genie::rew::GReWeightNuXSecCCRES( FitPar::Config().GetParS("GENIEXSecModelCCRES"), "Default")); #endif } if (xsec_ncres) fGenieRW->AdoptWghtCalc("xsec_ncres", new genie::rew::GReWeightNuXSecNCRES); if (xsec_nucqe) fGenieRW->AdoptWghtCalc("nuclear_qe", new genie::rew::GReWeightFGM); // Set the CCQE reweighting style GReWeightNuXSecCCQE *rwccqe = dynamic_cast(fGenieRW->WghtCalc("xsec_ccqe")); // For MaCCQE reweighting std::string ccqetype = FitPar::Config().GetParS("GENIEWeightEngine_CCQEMode"); if (ccqetype == "kModeMa") { NUIS_LOG(DEB, "Setting GENIE ReWeight CCQE to kModeMa"); rwccqe->SetMode(GReWeightNuXSecCCQE::kModeMa); } else if (ccqetype == "kModeNormAndMaShape") { NUIS_LOG(DEB, "Setting GENIE ReWeight CCQE to kModeNormAndMaShape"); rwccqe->SetMode(GReWeightNuXSecCCQE::kModeNormAndMaShape); // For z-expansion reweighting, only available after 2.10 #if __GENIE_VERSION__ >= 210 } else if (ccqetype == "kModeZExp") { NUIS_LOG(DEB, "Setting GENIE ReWeight CCQE to kModeZExp"); rwccqe->SetMode(GReWeightNuXSecCCQE::kModeZExp); #endif } else { NUIS_ERR(FTL, "Did not find specified GENIE ReWeight CCQE mode"); NUIS_ABORT("You provided: " << ccqetype << " in parameters/config.xml"); } #if (__GENIE_VERSION__ >= 212) and \ (__GENIE_VERSION__ <= \ 300) // This doesn't currently work as is for GENIE v3, but the reweighting // in v3 supposedly does similar checks anyway. // Check the UserPhysicsOptions too! AlgConfigPool *Pool = genie::AlgConfigPool::Instance(); Registry *full = Pool->GlobalParameterList(); std::string name_ax = full->GetAlg("AxialFormFactorModel").name; std::string config_ax = full->GetAlg("AxialFormFactorModel").config; if (name_ax == "genie::DipoleAxialFormFactorModel" && ccqetype == "kModeZExp") { + StartTalking(); NUIS_ERR( FTL, "Trying to run Z Expansion reweighting with Llewelyn-Smith model."); NUIS_ERR(FTL, "Please check your " << std::getenv("GENIE") << "/config/UserPhysicsOptions.xml to match generated"); NUIS_ERR(FTL, "You're telling me " << name_ax << "/" << config_ax); NUIS_ABORT("Also check your " << std::getenv("NUISANCE") << "/parameters/config.xml GENIEWeightEngine_CCQEMode: " << ccqetype); } if (name_ax == "genie::ZExpAxialFormFactorModel" && ccqetype != "kModeZExp") { + StartTalking(); NUIS_ERR( FTL, "Trying to run Llewelyn-Smith reweighting with Z Expansion model."); NUIS_ERR(FTL, "Please change your " << std::getenv("GENIE") << "/config/UserPhysicsOptions.xml to match generated"); NUIS_ERR(FTL, "You're telling me " << name_ax << "/" << config_ax); NUIS_ABORT("Also check your " << std::getenv("NUISANCE") << "/parameters/config.xml GENIEWeightEngine_CCQEMode: " << ccqetype); } std::string name_qelcc = full->GetAlg("XSecModel@genie::EventGenerator/QEL-CC").name; std::string config_qelcc = full->GetAlg("XSecModel@genie::EventGenerator/QEL-CC").config; if (config_qelcc == "Default" && ccqetype == "kModeZExp") { + StartTalking(); NUIS_ERR( FTL, "Trying to run Z Expansion reweighting with Llewelyn-Smith model."); NUIS_ERR(FTL, "Please change your " << std::getenv("GENIE") << "/config/UserPhysicsOptions.xml to match generated"); NUIS_ERR(FTL, "You're telling me " << name_qelcc << "/" << config_qelcc); NUIS_ABORT("Also check your " << std::getenv("NUISANCE") << "/parameters/config.xml GENIEWeightEngine_CCQEMode: " << ccqetype); } if (config_qelcc == "ZExp" && ccqetype != "kModeZExp") { + StartTalking(); NUIS_ERR( FTL, "Trying to run Llewelyn-Smith reweighting with Z Expansion model."); NUIS_ERR(FTL, "Please change your " << std::getenv("GENIE") << "/config/UserPhysicsOptions.xml to match generated"); NUIS_ERR(FTL, "You're telling me " << name_qelcc << "/" << config_qelcc); NUIS_ABORT("Also check your " << std::getenv("NUISANCE") << "/parameters/config.xml GENIEWeightEngine_CCQEMode: " << ccqetype); } #endif if (xsec_ccres) { // Default to include shape and normalization changes for CCRES (can be // changed downstream if desired) GReWeightNuXSecCCRES *rwccres = dynamic_cast(fGenieRW->WghtCalc("xsec_ccres")); std::string marestype = FitPar::Config().GetParS("GENIEWeightEngine_CCRESMode"); if (!marestype.compare("kModeNormAndMaMvShape")) { rwccres->SetMode(GReWeightNuXSecCCRES::kModeNormAndMaMvShape); } else if (!marestype.compare("kModeMaMv")) { rwccres->SetMode(GReWeightNuXSecCCRES::kModeMaMv); } else { + StartTalking(); NUIS_ABORT("Unkown MARES Mode in GENIE Weight Engine : " << marestype); } } if (xsec_ncres) { // Default to include shape and normalization changes for NCRES (can be // changed downstream if desired) GReWeightNuXSecNCRES *rwncres = dynamic_cast(fGenieRW->WghtCalc("xsec_ncres")); rwncres->SetMode(GReWeightNuXSecNCRES::kModeMaMv); } if (xsec_dis) { // Default to include shape and normalization changes for DIS (can be // changed downstream if desired) GReWeightNuXSecDIS *rwdis = dynamic_cast(fGenieRW->WghtCalc("xsec_dis")); rwdis->SetMode(GReWeightNuXSecDIS::kModeABCV12u); // Set Abs Twk Config fIsAbsTwk = (FitPar::Config().GetParB("setabstwk")); } // allow cout again StartTalking(); #else NUIS_ERR(FTL, "GENIE ReWeight is __NOT ENABLED__ in GENIE and you're trying to " "run NUISANCE with it enabled"); NUIS_ERR(FTL, "Check your genie-config --libs for reweighting"); NUIS_ERR(FTL, "You might also have run with -DUSE_REWEIGHT=0 (false) in your cmake build?"); NUIS_ERR(FTL, "If not present you need to recompile GENIE"); NUIS_ABORT("If present you need to contact NUISANCE authors"); #endif #endif }; void GENIEWeightEngine::IncludeDial(std::string name, double startval) { #ifdef __GENIE_ENABLED__ #ifndef __NO_REWEIGHT__ // Get First enum int nuisenum = Reweight::ConvDial(name, kGENIE); // Check ZExp sillyness in GENIE // If ZExpansion parameters are used we need to set a different mode in GENIE // ReWeight... GENIE doesn't have a setter either... #if (__GENIE_VERSION__ >= 212) and (__GENIE_VERSION__ <= 300) std::string ccqetype = FitPar::Config().GetParS("GENIEWeightEngine_CCQEMode"); if (ccqetype != "kModeZExp" && (name == "ZExpA1CCQE" || name == "ZExpA2CCQE" || name == "ZExpA3CCQE" || name == "ZExpA4CCQE")) { NUIS_ERR(FTL, "Found a Z-expansion parameter in GENIE although the GENIE " "ReWeighting engine is set to use Llewelyn-Smith and MaQE!"); NUIS_ABORT("Change your GENIE UserPhysicsOptions.xml in " << std::getenv("GENIE") << "/config/UserPhysicsOptions.xml to match requirements"); } if ((ccqetype != "kModeMa" && ccqetype != "kModeMaNormAndMaShape") && (name == "MaCCQE")) { NUIS_ERR(FTL, "Found MaCCQE parameter in GENIE although the GENIE " "ReWeighting engine is set to not use this!"); NUIS_ABORT("Change your GENIE UserPhysicsOptions.xml in " << std::getenv("GENIE") << "/config/UserPhysicsOptions.xml to match requirements"); } #endif // Setup Maps fEnumIndex[nuisenum]; // = std::vector(0); fNameIndex[name]; // = std::vector(0); // Split by commas std::vector allnames = GeneralUtils::ParseToStr(name, ","); for (uint i = 0; i < allnames.size(); i++) { std::string singlename = allnames[i]; // Get RW genie::rew::GSyst_t rwsyst = GSyst::FromString(singlename); // Fill Maps int index = fValues.size(); fValues.push_back(0.0); fGENIESysts.push_back(rwsyst); // Initialize dial NUIS_LOG(DEB, "Registering " << singlename << " from " << name); fGenieRW->Systematics().Init(fGENIESysts[index]); // If Absolute if (fIsAbsTwk) { GSystUncertainty::Instance()->SetUncertainty(rwsyst, 1.0, 1.0); } // Setup index fEnumIndex[nuisenum].push_back(index); fNameIndex[name].push_back(index); } // Set Value if given if (startval != _UNDEF_DIAL_VALUE_) { SetDialValue(nuisenum, startval); } #endif #endif } void GENIEWeightEngine::SetDialValue(int nuisenum, double val) { #ifdef __GENIE_ENABLED__ #ifndef __NO_REWEIGHT__ std::vector indices = fEnumIndex[nuisenum]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; fGenieRW->Systematics().Set(fGENIESysts[indices[i]], val); } #endif #endif } void GENIEWeightEngine::SetDialValue(std::string name, double val) { #ifdef __GENIE_ENABLED__ #ifndef __NO_REWEIGHT__ std::vector indices = fNameIndex[name]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; fGenieRW->Systematics().Set(fGENIESysts[indices[i]], val); } #endif #endif } void GENIEWeightEngine::Reconfigure(bool silent) { #ifdef __GENIE_ENABLED__ #ifndef __NO_REWEIGHT__ // Hush now... if (silent) StopTalking(); // Reconf fGenieRW->Reconfigure(); fGenieRW->Print(); // Shout again if (silent) StartTalking(); #endif #endif } double GENIEWeightEngine::CalcWeight(BaseFitEvt *evt) { double rw_weight = 1.0; #ifdef __GENIE_ENABLED__ #ifndef __NO_REWEIGHT__ // Make nom weight if (!evt) { NUIS_ABORT("evt not found : " << evt); } // Skip Non GENIE if (evt->fType != kGENIE) return 1.0; if (!(evt->genie_event)) { NUIS_ABORT("evt->genie_event not found!" << evt->genie_event); } if (!(evt->genie_event->event)) { NUIS_ABORT("evt->genie_event->event GHepRecord not found!" << (evt->genie_event->event)); } if (!fGenieRW) { NUIS_ABORT("GENIE RW Not Found!" << fGenieRW); } rw_weight = fGenieRW->CalcWeight(*(evt->genie_event->event)); // std::cout << "Returning GENIE Weight for electron scattering = " << // rw_weight << std::endl; // if (rw_weight != 1.0 )std::cout << "mode=" << evt->Mode << " rw_weight = " // << rw_weight << std::endl; #endif #endif // Return rw_weight return rw_weight; } diff --git a/src/Routines/SystematicRoutines.cxx b/src/Routines/SystematicRoutines.cxx index a89fe2f..cf53d85 100755 --- a/src/Routines/SystematicRoutines.cxx +++ b/src/Routines/SystematicRoutines.cxx @@ -1,1492 +1,1493 @@ // Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #include "SystematicRoutines.h" void SystematicRoutines::Init() { fInputFile = ""; fInputRootFile = NULL; fOutputFile = ""; fOutputRootFile = NULL; fCovar = fCovarFree = NULL; fCorrel = fCorrelFree = NULL; fDecomp = fDecompFree = NULL; fStrategy = "ErrorBands"; fRoutines.clear(); fRoutines.push_back("ErrorBands"); fCardFile = ""; fFakeDataInput = ""; fSampleFCN = NULL; fAllowedRoutines = ("ErrorBands,PlotLimits"); }; SystematicRoutines::~SystematicRoutines(){}; SystematicRoutines::SystematicRoutines(int argc, char *argv[]) { // Initialise Defaults Init(); nuisconfig configuration = Config::Get(); // Default containers std::string cardfile = ""; std::string maxevents = "-1"; int errorcount = 0; int verbocount = 0; std::vector xmlcmds; std::vector configargs; fNThrows = 250; fStartThrows = 0; fThrowString = ""; // Make easier to handle arguments. std::vector args = GeneralUtils::LoadCharToVectStr(argc, argv); ParserUtils::ParseArgument(args, "-c", fCardFile, true); ParserUtils::ParseArgument(args, "-o", fOutputFile, false, false); ParserUtils::ParseArgument(args, "-n", maxevents, false, false); ParserUtils::ParseArgument(args, "-f", fStrategy, false, false); ParserUtils::ParseArgument(args, "-d", fFakeDataInput, false, false); ParserUtils::ParseArgument(args, "-s", fStartThrows, false, false); ParserUtils::ParseArgument(args, "-t", fNThrows, false, false); ParserUtils::ParseArgument(args, "-p", fThrowString, false, false); ParserUtils::ParseArgument(args, "-i", xmlcmds); ParserUtils::ParseArgument(args, "-q", configargs); ParserUtils::ParseCounter(args, "e", errorcount); ParserUtils::ParseCounter(args, "v", verbocount); ParserUtils::CheckBadArguments(args); // Add extra defaults if none given if (fCardFile.empty() and xmlcmds.empty()) { NUIS_ABORT("No input supplied!"); } if (fOutputFile.empty() and !fCardFile.empty()) { fOutputFile = fCardFile + ".root"; NUIS_ERR(WRN, "No output supplied so saving it to: " << fOutputFile); } else if (fOutputFile.empty()) { NUIS_ABORT("No output file or cardfile supplied!"); } // Configuration Setup ============================= // Check no comp key is available if (Config::Get().GetNodes("nuiscomp").empty()) { fCompKey = Config::Get().CreateNode("nuiscomp"); } else { fCompKey = Config::Get().GetNodes("nuiscomp")[0]; } if (!fCardFile.empty()) fCompKey.Set("cardfile", fCardFile); if (!fOutputFile.empty()) fCompKey.Set("outputfile", fOutputFile); if (!fStrategy.empty()) fCompKey.Set("strategy", fStrategy); // Load XML Cardfile configuration.LoadSettings(fCompKey.GetS("cardfile"), ""); // Add Config Args for (size_t i = 0; i < configargs.size(); i++) { configuration.OverrideConfig(configargs[i]); } if (maxevents.compare("-1")) { configuration.OverrideConfig("MAXEVENTS=" + maxevents); } // Finish configuration XML configuration.FinaliseSettings(fCompKey.GetS("outputfile") + ".xml"); // Add Error Verbo Lines verbocount += Config::GetParI("VERBOSITY"); errorcount += Config::GetParI("ERROR"); std::cout << "[ NUISANCE ]: Setting VERBOSITY=" << verbocount << std::endl; std::cout << "[ NUISANCE ]: Setting ERROR=" << errorcount << std::endl; SETVERBOSITY(verbocount); // Proper Setup if (fStrategy.find("ErrorBands") != std::string::npos || fStrategy.find("MergeErrors") != std::string::npos) { fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE"); } // fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), // "RECREATE"); SetupSystematicsFromXML(); SetupCovariance(); SetupRWEngine(); SetupFCN(); GetCovarFromFCN(); // Run(); return; }; void SystematicRoutines::SetupSystematicsFromXML() { NUIS_LOG(FIT, "Setting up nuismin"); // Setup Parameters ------------------------------------------ std::vector parkeys = Config::QueryKeys("parameter"); if (!parkeys.empty()) { NUIS_LOG(FIT, "Number of parameters : " << parkeys.size()); } for (size_t i = 0; i < parkeys.size(); i++) { nuiskey key = parkeys.at(i); // Check for type,name,nom if (!key.Has("type")) { NUIS_ABORT("No type given for parameter " << i); } else if (!key.Has("name")) { NUIS_ABORT("No name given for parameter " << i); } else if (!key.Has("nominal")) { NUIS_ABORT("No nominal given for parameter " << i); } // Get Inputs std::string partype = key.GetS("type"); std::string parname = key.GetS("name"); double parnom = key.GetD("nominal"); double parlow = parnom - 1; double parhigh = parnom + 1; double parstep = 1; // Override if state not given if (!key.Has("state")) { key.SetS("state", "FIX"); } std::string parstate = key.GetS("state"); // Extra limits if (key.Has("low")) { parlow = key.GetD("low"); parhigh = key.GetD("high"); parstep = key.GetD("step"); NUIS_LOG(FIT, "Read " << partype << " : " << parname << " = " << parnom << " : " << parlow << " < p < " << parhigh << " : " << parstate); } else { NUIS_LOG(FIT, "Read " << partype << " : " << parname << " = " << parnom << " : " << parstate); } // Run Parameter Conversion if needed if (parstate.find("ABS") != std::string::npos) { parnom = FitBase::RWAbsToSigma(partype, parname, parnom); parlow = FitBase::RWAbsToSigma(partype, parname, parlow); parhigh = FitBase::RWAbsToSigma(partype, parname, parhigh); parstep = FitBase::RWAbsToSigma(partype, parname, parstep); } else if (parstate.find("FRAC") != std::string::npos) { parnom = FitBase::RWFracToSigma(partype, parname, parnom); parlow = FitBase::RWFracToSigma(partype, parname, parlow); parhigh = FitBase::RWFracToSigma(partype, parname, parhigh); parstep = FitBase::RWFracToSigma(partype, parname, parstep); } // Push into vectors fParams.push_back(parname); fTypeVals[parname] = FitBase::ConvDialType(partype); fStartVals[parname] = parnom; fCurVals[parname] = parnom; fErrorVals[parname] = 0.0; fStateVals[parname] = parstate; bool fixstate = parstate.find("FIX") != std::string::npos; fFixVals[parname] = fixstate; fStartFixVals[parname] = fFixVals[parname]; fMinVals[parname] = parlow; fMaxVals[parname] = parhigh; fStepVals[parname] = parstep; } // Setup Samples ---------------------------------------------- std::vector samplekeys = Config::QueryKeys("sample"); if (!samplekeys.empty()) { NUIS_LOG(FIT, "Number of samples : " << samplekeys.size()); } for (size_t i = 0; i < samplekeys.size(); i++) { nuiskey key = samplekeys.at(i); // Get Sample Options std::string samplename = key.GetS("name"); std::string samplefile = key.GetS("input"); std::string sampletype = key.Has("type") ? key.GetS("type") : "DEFAULT"; double samplenorm = key.Has("norm") ? key.GetD("norm") : 1.0; // Print out NUIS_LOG(FIT, "Read sample info " << i << " : " << samplename << std::endl << "\t\t input -> " << samplefile << std::endl << "\t\t state -> " << sampletype << std::endl << "\t\t norm -> " << samplenorm); // If FREE add to parameters otherwise continue if (sampletype.find("FREE") == std::string::npos) { continue; } // Form norm dial from samplename + sampletype + "_norm"; std::string normname = samplename + "_norm"; // Check normname not already present if (fTypeVals.find(normname) != fTypeVals.end()) { continue; } // Add new norm dial to list if its passed above checks fParams.push_back(normname); fTypeVals[normname] = kNORM; fStateVals[normname] = sampletype; fCurVals[normname] = samplenorm; fErrorVals[normname] = 0.0; fMinVals[normname] = 0.1; fMaxVals[normname] = 10.0; fStepVals[normname] = 0.5; bool state = sampletype.find("FREE") == std::string::npos; fFixVals[normname] = state; fStartFixVals[normname] = state; } // Setup Fake Parameters ----------------------------- std::vector fakekeys = Config::QueryKeys("fakeparameter"); if (!fakekeys.empty()) { NUIS_LOG(FIT, "Number of fake parameters : " << fakekeys.size()); } for (size_t i = 0; i < fakekeys.size(); i++) { nuiskey key = fakekeys.at(i); // Check for type,name,nom if (!key.Has("name")) { NUIS_ABORT("No name given for fakeparameter " << i); } else if (!key.Has("nom")) { NUIS_ABORT("No nominal given for fakeparameter " << i); } // Get Inputs std::string parname = key.GetS("name"); double parnom = key.GetD("nom"); // Push into vectors fFakeVals[parname] = parnom; } } /* Setup Functions */ //************************************* void SystematicRoutines::SetupRWEngine() { //************************************* for (UInt_t i = 0; i < fParams.size(); i++) { std::string name = fParams[i]; FitBase::GetRW()->IncludeDial(name, fTypeVals.at(name)); } UpdateRWEngine(fStartVals); return; } //************************************* void SystematicRoutines::SetupFCN() { //************************************* NUIS_LOG(FIT, "Making the jointFCN"); if (fSampleFCN) delete fSampleFCN; fSampleFCN = new JointFCN(fOutputRootFile); + fSampleFCN->SetNParams(int(fParams.size())); SetFakeData(); return; } //************************************* void SystematicRoutines::SetFakeData() { //************************************* if (fFakeDataInput.empty()) return; if (fFakeDataInput.compare("MC") == 0) { NUIS_LOG(FIT, "Setting fake data from MC starting prediction."); UpdateRWEngine(fFakeVals); FitBase::GetRW()->Reconfigure(); fSampleFCN->ReconfigureAllEvents(); fSampleFCN->SetFakeData("MC"); UpdateRWEngine(fCurVals); NUIS_LOG(FIT, "Set all data to fake MC predictions."); } else { fSampleFCN->SetFakeData(fFakeDataInput); } return; } //***************************************** // Setup the parameter covariances from the FCN void SystematicRoutines::GetCovarFromFCN() { //***************************************** NUIS_LOG(FIT, "Loading ParamPull objects from FCN to build covariance..."); // Make helperstring std::ostringstream helperstr; // Keep track of what is being thrown std::map dialthrowhandle; // Get Covariance Objects from FCN std::list inputpulls = fSampleFCN->GetPullList(); for (PullListConstIter iter = inputpulls.begin(); iter != inputpulls.end(); iter++) { ParamPull *pull = (*iter); if (pull->GetType().find("THROW") != std::string::npos) { fInputThrows.push_back(pull); fInputCovar.push_back(pull->GetFullCovarMatrix()); fInputDials.push_back(pull->GetDataHist()); NUIS_LOG(FIT, "Read ParamPull: " << pull->GetName() << " " << pull->GetType()); } TH1D dialhist = pull->GetDataHist(); TH1D minhist = pull->GetMinHist(); TH1D maxhist = pull->GetMaxHist(); TH1I typehist = pull->GetDialTypes(); for (int i = 0; i < dialhist.GetNbinsX(); i++) { std::string name = std::string(dialhist.GetXaxis()->GetBinLabel(i + 1)); dialthrowhandle[name] = pull->GetName(); // Add to Containers // Set the starting values etc to the postfit fParams.push_back(name); fCurVals[name] = dialhist.GetBinContent(i + 1); // Set the starting values to be nominal MC fStartVals[name] = 0.0; fMinVals[name] = minhist.GetBinContent(i + 1); fMaxVals[name] = maxhist.GetBinContent(i + 1); fStepVals[name] = 1.0; fFixVals[name] = false; fStartFixVals[name] = false; fTypeVals[name] = typehist.GetBinContent(i + 1); fStateVals[name] = "FREE," + pull->GetType(); // If we find the string if (fCurVals.find(name) == fCurVals.end()) { // Maker Helper helperstr << std::string(16, ' ') << FitBase::ConvDialType(fTypeVals[name]) << " " << name << " " << fMinVals[name] << " " << fMaxVals[name] << " " << fStepVals[name] << " " << fStateVals[name] << std::endl; } } } // Check if no throws given if (fInputThrows.empty()) { NUIS_ERR(WRN, "No covariances given to nuissyst"); NUIS_ERR(WRN, "Pushing back an uncorrelated gaussian throw error for each " "free parameter using step size"); for (UInt_t i = 0; i < fParams.size(); i++) { std::string syst = fParams[i]; if (fFixVals[syst]) continue; // Make Terms std::string name = syst + "_pull"; std::ostringstream pullterm; pullterm << "DIAL:" << syst << ";" << fStartVals[syst] << ";" << fStepVals[syst]; std::string type = "GAUSTHROW/NEUT"; // Push Back Pulls ParamPull *pull = new ParamPull(name, pullterm.str(), type); fInputThrows.push_back(pull); fInputCovar.push_back(pull->GetFullCovarMatrix()); fInputDials.push_back(pull->GetDataHist()); // Print Whats added NUIS_LOG(FIT, "Added ParamPull : " << name << " " << pullterm.str() << " " << type); // Add helper string for future fits helperstr << std::string(16, ' ') << "covar " << name << " " << pullterm.str() << " " << type << std::endl; // Keep Track of Throws dialthrowhandle[syst] = pull->GetName(); } } // Print Helper String if (!helperstr.str().empty()) { NUIS_LOG(FIT, "To remove these statements in future studies, add the lines " "below to your card:"); // Can't use the logger properly because this can be multi-line. Use cout // and added spaces to look better! std::cout << helperstr.str(); sleep(2); } // Print Throw State for (UInt_t i = 0; i < fParams.size(); i++) { std::string syst = fParams[i]; if (dialthrowhandle.find(syst) != dialthrowhandle.end()) { NUIS_LOG(FIT, "Dial " << i << ". " << setw(20) << syst << " = THROWing with " << dialthrowhandle[syst]); } else { NUIS_LOG(FIT, "Dial " << i << ". " << setw(20) << syst << " = is FIXED"); } } } /* Fitting Functions */ //************************************* void SystematicRoutines::UpdateRWEngine( std::map &updateVals) { //************************************* for (UInt_t i = 0; i < fParams.size(); i++) { std::string name = fParams[i]; if (updateVals.find(name) == updateVals.end()) continue; FitBase::GetRW()->SetDialValue(name, updateVals.at(name)); } FitBase::GetRW()->Reconfigure(); } //************************************* void SystematicRoutines::PrintState() { //************************************* NUIS_LOG(FIT, "------------"); // Count max size int maxcount = 0; for (UInt_t i = 0; i < fParams.size(); i++) { maxcount = max(int(fParams[i].size()), maxcount); } // Header NUIS_LOG(FIT, " # " << left << setw(maxcount) << "Parameter " << " = " << setw(10) << "Value" << " +- " << setw(10) << "Error" << " " << setw(8) << "(Units)" << " " << setw(10) << "Conv. Val" << " +- " << setw(10) << "Conv. Err" << " " << setw(8) << "(Units)"); // Parameters for (UInt_t i = 0; i < fParams.size(); i++) { std::string syst = fParams.at(i); std::string typestr = FitBase::ConvDialType(fTypeVals[syst]); std::string curunits = "(sig.)"; double curval = fCurVals[syst]; double curerr = fErrorVals[syst]; if (fStateVals[syst].find("ABS") != std::string::npos) { curval = FitBase::RWSigmaToAbs(typestr, syst, curval); curerr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) - FitBase::RWSigmaToAbs(typestr, syst, 0.0)); curunits = "(Abs.)"; } else if (fStateVals[syst].find("FRAC") != std::string::npos) { curval = FitBase::RWSigmaToFrac(typestr, syst, curval); curerr = (FitBase::RWSigmaToFrac(typestr, syst, curerr) - FitBase::RWSigmaToFrac(typestr, syst, 0.0)); curunits = "(Frac)"; } std::string convunits = "(" + FitBase::GetRWUnits(typestr, syst) + ")"; double convval = FitBase::RWSigmaToAbs(typestr, syst, curval); double converr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) - FitBase::RWSigmaToAbs(typestr, syst, 0.0)); std::ostringstream curparstring; curparstring << " " << setw(3) << left << i << ". " << setw(maxcount) << syst << " = " << setw(10) << curval << " +- " << setw(10) << curerr << " " << setw(8) << curunits << " " << setw(10) << convval << " +- " << setw(10) << converr << " " << setw(8) << convunits; NUIS_LOG(FIT, curparstring.str()); } NUIS_LOG(FIT, "------------"); double like = fSampleFCN->GetLikelihood(); NUIS_LOG(FIT, std::left << std::setw(46) << "Likelihood for JointFCN: " << like); NUIS_LOG(FIT, "------------"); } /* Write Functions */ //************************************* void SystematicRoutines::SaveResults() { //************************************* if (!fOutputRootFile) fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE"); fOutputRootFile->cd(); SaveCurrentState(); } //************************************* void SystematicRoutines::SaveCurrentState(std::string subdir) { //************************************* NUIS_LOG(FIT, "Saving current full FCN predictions"); // Setup DIRS TDirectory *curdir = gDirectory; if (!subdir.empty()) { TDirectory *newdir = (TDirectory *)gDirectory->mkdir(subdir.c_str()); newdir->cd(); } FitBase::GetRW()->Reconfigure(); fSampleFCN->ReconfigureAllEvents(); fSampleFCN->Write(); // Change back to current DIR curdir->cd(); return; } //************************************* void SystematicRoutines::SaveNominal() { //************************************* if (!fOutputRootFile) fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE"); fOutputRootFile->cd(); NUIS_LOG(FIT, "Saving Nominal Predictions (be cautious with this)"); FitBase::GetRW()->Reconfigure(); SaveCurrentState("nominal"); }; //************************************* void SystematicRoutines::SavePrefit() { //************************************* if (!fOutputRootFile) fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE"); fOutputRootFile->cd(); NUIS_LOG(FIT, "Saving Prefit Predictions"); UpdateRWEngine(fStartVals); SaveCurrentState("prefit"); UpdateRWEngine(fCurVals); }; /* MISC Functions */ //************************************* int SystematicRoutines::GetStatus() { //************************************* return 0; } //************************************* void SystematicRoutines::SetupCovariance() { //************************************* // Remove covares if they exist if (fCovar) delete fCovar; if (fCovarFree) delete fCovarFree; if (fCorrel) delete fCorrel; if (fCorrelFree) delete fCorrelFree; if (fDecomp) delete fDecomp; if (fDecompFree) delete fDecompFree; int NFREE = 0; int NDIM = 0; // Get NFREE from min or from vals (for cases when doing throws) NDIM = fParams.size(); for (UInt_t i = 0; i < fParams.size(); i++) { if (!fFixVals[fParams[i]]) NFREE++; } if (NDIM == 0) return; fCovar = new TH2D("covariance", "covariance", NDIM, 0, NDIM, NDIM, 0, NDIM); if (NFREE > 0) { fCovarFree = new TH2D("covariance_free", "covariance_free", NFREE, 0, NFREE, NFREE, 0, NFREE); } // Set Bin Labels int countall = 0; int countfree = 0; for (UInt_t i = 0; i < fParams.size(); i++) { fCovar->GetXaxis()->SetBinLabel(countall + 1, fParams[i].c_str()); fCovar->GetYaxis()->SetBinLabel(countall + 1, fParams[i].c_str()); countall++; if (!fFixVals[fParams[i]] and NFREE > 0) { fCovarFree->GetXaxis()->SetBinLabel(countfree + 1, fParams[i].c_str()); fCovarFree->GetYaxis()->SetBinLabel(countfree + 1, fParams[i].c_str()); countfree++; } } fCorrel = PlotUtils::GetCorrelationPlot(fCovar, "correlation"); fDecomp = PlotUtils::GetDecompPlot(fCovar, "decomposition"); if (NFREE > 0) fCorrelFree = PlotUtils::GetCorrelationPlot(fCovarFree, "correlation_free"); if (NFREE > 0) fDecompFree = PlotUtils::GetDecompPlot(fCovarFree, "decomposition_free"); return; }; //************************************* void SystematicRoutines::ThrowCovariance(bool uniformly) { //************************************* // Set fThrownVals to all values in currentVals for (UInt_t i = 0; i < fParams.size(); i++) { std::string name = fParams.at(i); fThrownVals[name] = fCurVals[name]; } for (PullListConstIter iter = fInputThrows.begin(); iter != fInputThrows.end(); iter++) { ParamPull *pull = *iter; pull->ThrowCovariance(); TH1D dialhist = pull->GetDataHist(); for (int i = 0; i < dialhist.GetNbinsX(); i++) { std::string name = std::string(dialhist.GetXaxis()->GetBinLabel(i + 1)); if (fCurVals.find(name) != fCurVals.end()) { fThrownVals[name] = dialhist.GetBinContent(i + 1); } } // Reset throw in case pulls are calculated. pull->ResetToy(); } }; //************************************* void SystematicRoutines::PlotLimits() { //************************************* std::cout << "Plotting Limits" << std::endl; if (!fOutputRootFile) fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE"); TDirectory *limfolder = (TDirectory *)fOutputRootFile->mkdir("Limits"); limfolder->cd(); // Set all parameters at their starting values for (UInt_t i = 0; i < fParams.size(); i++) { fCurVals[fParams[i]] = fStartVals[fParams[i]]; } TDirectory *nomfolder = (TDirectory *)limfolder->mkdir("nominal"); nomfolder->cd(); UpdateRWEngine(fCurVals); fSampleFCN->ReconfigureAllEvents(); fSampleFCN->Write(); limfolder->cd(); std::vector allfolders; // Loop through each parameter for (UInt_t i = 0; i < fParams.size(); i++) { std::string syst = fParams[i]; std::cout << "Starting Param " << syst << std::endl; if (fFixVals[syst]) continue; // Loop Downwards while (fCurVals[syst] > fMinVals[syst]) { fCurVals[syst] = fCurVals[syst] - fStepVals[syst]; // Check Limit if (fCurVals[syst] < fMinVals[syst]) fCurVals[syst] = fMinVals[syst]; // Check folder exists std::string curvalstring = std::string(Form((syst + "_%f").c_str(), fCurVals[syst])); if (std::find(allfolders.begin(), allfolders.end(), curvalstring) != allfolders.end()) break; // Make new folder for variation TDirectory *minfolder = (TDirectory *)limfolder->mkdir( Form((syst + "_%f").c_str(), fCurVals[syst])); minfolder->cd(); allfolders.push_back(curvalstring); // Update Iterations double *vals = FitUtils::GetArrayFromMap(fParams, fCurVals); fSampleFCN->DoEval(vals); delete vals; // Save to folder fSampleFCN->Write(); } // Reset before next loop fCurVals[syst] = fStartVals[syst]; // Loop Upwards now while (fCurVals[syst] < fMaxVals[syst]) { fCurVals[syst] = fCurVals[syst] + fStepVals[syst]; // Check Limit if (fCurVals[syst] > fMaxVals[syst]) fCurVals[syst] = fMaxVals[syst]; // Check folder exists std::string curvalstring = std::string(Form((syst + "_%f").c_str(), fCurVals[syst])); if (std::find(allfolders.begin(), allfolders.end(), curvalstring) != allfolders.end()) break; // Make new folder TDirectory *maxfolder = (TDirectory *)limfolder->mkdir( Form((syst + "_%f").c_str(), fCurVals[syst])); maxfolder->cd(); allfolders.push_back(curvalstring); // Update Iterations double *vals = FitUtils::GetArrayFromMap(fParams, fCurVals); fSampleFCN->DoEval(vals); delete vals; // Save to file fSampleFCN->Write(); } // Reset before leaving fCurVals[syst] = fStartVals[syst]; UpdateRWEngine(fCurVals); } return; } //************************************* void SystematicRoutines::Run() { //************************************* fRoutines = GeneralUtils::ParseToStr(fStrategy, ","); for (UInt_t i = 0; i < fRoutines.size(); i++) { std::string routine = fRoutines.at(i); int fitstate = kFitUnfinished; NUIS_LOG(FIT, "Running Routine: " << routine); if (routine.compare("PlotLimits") == 0) PlotLimits(); else if (routine.compare("ErrorBands") == 0) GenerateErrorBands(); else if (routine.compare("ThrowErrors") == 0) GenerateThrows(); else if (routine.compare("MergeErrors") == 0) MergeThrows(); else if (routine.compare("EigenErrors") == 0) EigenErrors(); else { NUIS_ABORT("Unknown ROUTINE : " << routine); } // If ending early break here if (fitstate == kFitFinished || fitstate == kNoChange) { NUIS_LOG(FIT, "Ending fit routines loop."); break; } } } void SystematicRoutines::GenerateErrorBands() { GenerateThrows(); MergeThrows(); } //************************************* void SystematicRoutines::GenerateThrows() { //************************************* TFile *tempfile = new TFile((fOutputFile + ".throws.root").c_str(), "RECREATE"); tempfile->cd(); // For generating throws we check with the config int nthrows = Config::GetParI("error_throws"); int startthrows = fStartThrows; int endthrows = startthrows + nthrows; if (nthrows < 0) nthrows = endthrows; if (startthrows < 0) startthrows = 0; if (endthrows < 0) endthrows = startthrows + nthrows; // Setting Seed // Matteo Mazzanti's Fix struct timeval mytime; gettimeofday(&mytime, NULL); Double_t seed = time(NULL) + int(getpid()) + (mytime.tv_sec * 1000.) + (mytime.tv_usec / 1000.); gRandom->SetSeed(seed); // int seed = (gRandom->Uniform(0.0,1.0)*100000 + 100000000*(startthrows + // endthrows) + time(NULL) + int(getpid()) ); gRandom->SetSeed(seed); NUIS_LOG(FIT, "Using Seed : " << seed); NUIS_LOG(FIT, "nthrows = " << nthrows); NUIS_LOG(FIT, "startthrows = " << startthrows); NUIS_LOG(FIT, "endthrows = " << endthrows); UpdateRWEngine(fStartVals); fSampleFCN->ReconfigureAllEvents(); // Make the nominal if (startthrows == 0) { NUIS_LOG(FIT, "Making nominal "); TDirectory *nominal = (TDirectory *)tempfile->mkdir("nominal"); nominal->cd(); fSampleFCN->Write(); // Make the postfit reading from the pull NUIS_LOG(FIT, "Making postfit "); TDirectory *postfit = (TDirectory *)tempfile->mkdir("postfit"); postfit->cd(); UpdateRWEngine(fCurVals); fSampleFCN->ReconfigureSignal(); fSampleFCN->Write(); } fSampleFCN->CreateIterationTree("error_iterations", FitBase::GetRW()); // Would anybody actually want to do uniform throws of any parameter?? bool uniformly = FitPar::Config().GetParB("error_uniform"); // Run Throws and save for (Int_t i = 0; i < endthrows + 1; i++) { // Generate Random Parameter Throw ThrowCovariance(uniformly); if (i < startthrows) continue; if (i == 0) continue; NUIS_LOG(FIT, "Throw " << i << "/" << endthrows << " ================================"); TDirectory *throwfolder = (TDirectory *)tempfile->mkdir(Form("throw_%i", i)); throwfolder->cd(); // Run Eval double *vals = FitUtils::GetArrayFromMap(fParams, fThrownVals); fSampleFCN->DoEval(vals); delete vals; // Save the FCN fSampleFCN->Write(); } tempfile->cd(); fSampleFCN->WriteIterationTree(); tempfile->Close(); } // Merge throws together into one summary void SystematicRoutines::MergeThrows() { fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE"); fOutputRootFile->cd(); // Make a container folder TDirectory *errorDIR = (TDirectory *)fOutputRootFile->mkdir("error_bands"); errorDIR->cd(); TDirectory *outnominal = (TDirectory *)fOutputRootFile->mkdir("nominal_throw"); outnominal->cd(); // Split Input Files if (!fThrowString.empty()) fThrowList = GeneralUtils::ParseToStr(fThrowString, ","); // Add default if no throwlist given if (fThrowList.size() < 1) fThrowList.push_back(fOutputFile + ".throws.root"); /// Save location of file containing nominal std::string nominalfile; bool nominalfound; // Loop over files and check they exist. for (uint i = 0; i < fThrowList.size(); i++) { std::string file = fThrowList[i]; bool found = false; // normal std::string newfile = file; TFile *throwfile = new TFile(file.c_str(), "READ"); if (throwfile and !throwfile->IsZombie()) { found = true; } // normal.throws.root if (!found) { newfile = file + ".throws.root"; throwfile = new TFile((file + ".throws.root").c_str(), "READ"); if (throwfile and !throwfile->IsZombie()) { found = true; } } // If its found save to throwlist, else save empty. // Also search for nominal if (found) { fThrowList[i] = newfile; NUIS_LOG(FIT, "Throws File :" << newfile); // Find input which contains nominal if (throwfile->Get("nominal")) { nominalfound = true; nominalfile = newfile; } throwfile->Close(); } else { fThrowList[i] = ""; } delete throwfile; } // Make sure we have a nominal file if (!nominalfound or nominalfile.empty()) { NUIS_ABORT("No nominal found when merging! Exiting!"); } // Get the nominal throws file TFile *tempfile = new TFile((nominalfile).c_str(), "READ"); tempfile->cd(); TDirectory *nominal = (TDirectory *)tempfile->Get("nominal"); bool uniformly = FitPar::Config().GetParB("error_uniform"); // Check percentage of bad files is okay. int badfilecount = 0; for (uint i = 0; i < fThrowList.size(); i++) { if (!fThrowList[i].empty()) { NUIS_LOG(FIT, "Loading Throws From File " << i << " : " << fThrowList[i]); } else { badfilecount++; } } // Check we have at least one good file if ((uint)badfilecount == fThrowList.size()) { NUIS_ABORT("Found no good throw files for MergeThrows"); throw; } else if (badfilecount > (fThrowList.size() * 0.25)) { NUIS_ERR( WRN, "Over 25% of your throw files are dodgy. Please check this is okay!"); NUIS_ERR(WRN, "Will continue for the time being..."); sleep(5); } // Now go through the keys in the temporary file and look for TH1D, and TH2D // plots TIter next(nominal->GetListOfKeys()); TKey *key; while ((key = (TKey *)next())) { TClass *cl = gROOT->GetClass(key->GetClassName()); if (!cl->InheritsFrom("TH1D") and !cl->InheritsFrom("TH2D")) continue; TH1 *baseplot = (TH1D *)key->ReadObj(); std::string plotname = std::string(baseplot->GetName()); NUIS_LOG(FIT, "Creating error bands for " << plotname); if (LOG_LEVEL(FIT)) { if (!uniformly) { NUIS_LOG(FIT, " : Using COVARIANCE Throws! "); } else { NUIS_LOG(FIT, " : Using UNIFORM THROWS!!! "); } } int nbins = 0; if (cl->InheritsFrom("TH1D")) nbins = ((TH1D *)baseplot)->GetNbinsX(); else nbins = ((TH1D *)baseplot)->GetNbinsX() * ((TH1D *)baseplot)->GetNbinsY(); // Setup TProfile with RMS option TProfile *tprof = new TProfile((plotname + "_prof").c_str(), (plotname + "_prof").c_str(), nbins, 0, nbins, "S"); // Setup The TTREE double *bincontents; bincontents = new double[nbins]; double *binlowest; binlowest = new double[nbins]; double *binhighest; binhighest = new double[nbins]; errorDIR->cd(); TTree *bintree = new TTree((plotname + "_tree").c_str(), (plotname + "_tree").c_str()); for (Int_t i = 0; i < nbins; i++) { bincontents[i] = 0.0; binhighest[i] = 0.0; binlowest[i] = 0.0; bintree->Branch(Form("content_%i", i), &bincontents[i], Form("content_%i/D", i)); } // Make new throw plot TH1 *newplot; // Run Throw Merging. for (UInt_t i = 0; i < fThrowList.size(); i++) { TFile *throwfile = new TFile(fThrowList[i].c_str(), "READ"); // Loop over all throws in a folder TIter nextthrow(throwfile->GetListOfKeys()); TKey *throwkey; while ((throwkey = (TKey *)nextthrow())) { // Skip non throw folders if (std::string(throwkey->GetName()).find("throw_") == std::string::npos) continue; // Get Throw DIR TDirectory *throwdir = (TDirectory *)throwkey->ReadObj(); // Get Plot From Throw newplot = (TH1 *)throwdir->Get(plotname.c_str()); if (!newplot) continue; // Loop Over Plot for (Int_t j = 0; j < nbins; j++) { tprof->Fill(j + 0.5, newplot->GetBinContent(j + 1)); bincontents[j] = newplot->GetBinContent(j + 1); if (bincontents[j] < binlowest[j] or i == 0) binlowest[j] = bincontents[j]; if (bincontents[j] > binhighest[j] or i == 0) binhighest[j] = bincontents[j]; } errorDIR->cd(); bintree->Fill(); } throwfile->Close(); delete throwfile; } errorDIR->cd(); if (uniformly) { NUIS_LOG(FIT, "Uniformly Calculating Plot Errors!"); } TH1 *statplot = (TH1 *)baseplot->Clone(); for (Int_t j = 0; j < nbins; j++) { if (!uniformly) { // if ((baseplot->GetBinError(j+1)/baseplot->GetBinContent(j+1)) //< 1.0) { // baseplot->SetBinError(j+1,sqrt(pow(tprof->GetBinError(j+1),2) //+ pow(baseplot->GetBinError(j+1),2))); } else { baseplot->SetBinContent(j + 1, tprof->GetBinContent(j + 1)); baseplot->SetBinError(j + 1, tprof->GetBinError(j + 1)); // } } else { baseplot->SetBinContent(j + 1, 0.0); //(binlowest[j] + binhighest[j]) / 2.0); baseplot->SetBinError(j + 1, 0.0); //(binhighest[j] - binlowest[j])/2.0); } } baseplot->SetTitle("Profiled throws"); errorDIR->cd(); baseplot->Write(); tprof->Write(); bintree->Write(); outnominal->cd(); for (int i = 0; i < nbins; i++) { baseplot->SetBinError(i + 1, sqrt(pow(statplot->GetBinError(i + 1), 2) + pow(baseplot->GetBinError(i + 1), 2))); } baseplot->Write(); delete statplot; delete baseplot; delete tprof; delete bintree; delete[] bincontents; } fOutputRootFile->Write(); fOutputRootFile->Close(); }; void SystematicRoutines::EigenErrors() { fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE"); fOutputRootFile->cd(); // Make Covariance TMatrixDSym *fullcovar = new TMatrixDSym(fParams.size()); // Extract covariance from all loaded ParamPulls for (PullListConstIter iter = fInputThrows.begin(); iter != fInputThrows.end(); iter++) { ParamPull *pull = *iter; // Check pull is actualyl Gaussian std::string pulltype = pull->GetType(); if (pulltype.find("GAUSTHROW") == std::string::npos) { NUIS_ABORT("Can only calculate EigenErrors for Gaussian pulls!"); } // Get data and covariances TH1D dialhist = pull->GetDataHist(); TH2D covhist = pull->GetFullCovar(); // Loop over all dials and compare names for (size_t pari = 0; pari < fParams.size(); pari++) { for (size_t parj = 0; parj < fParams.size(); parj++) { std::string name_pari = fParams[pari]; std::string name_parj = fParams[parj]; // Compare names to those in the pull for (int pulli = 0; pulli < dialhist.GetNbinsX(); pulli++) { for (int pullj = 0; pullj < dialhist.GetNbinsX(); pullj++) { std::string name_pulli = dialhist.GetXaxis()->GetBinLabel(pulli + 1); std::string name_pullj = dialhist.GetXaxis()->GetBinLabel(pullj + 1); if (name_pulli == name_pari && name_pullj == name_parj) { (*fullcovar)[pari][parj] = covhist.GetBinContent(pulli + 1, pullj + 1); fCurVals[name_pari] = dialhist.GetBinContent(pulli + 1); fCurVals[name_parj] = dialhist.GetBinContent(pullj + 1); } } } } } } /* TFile* test = new TFile("testingcovar.root","RECREATE"); test->cd(); TH2D* joinedcov = new TH2D("COVAR","COVAR", fullcovar->GetNrows(), 0.0, float(fullcovar->GetNrows()), fullcovar->GetNrows(), 0.0, float(fullcovar->GetNrows())); for (int i = 0; i < fullcovar->GetNrows(); i++){ for (int j = 0; j < fullcovar->GetNcols(); j++){ joinedcov->SetBinContent(i+1, j+1, (*fullcovar)[i][j]); } } joinedcov->Write("COVAR"); test->Close(); */ // Calculator all EigenVectors and EigenValues TMatrixDSymEigen *eigen = new TMatrixDSymEigen(*fullcovar); const TVectorD eigenVals = eigen->GetEigenValues(); const TMatrixD eigenVect = eigen->GetEigenVectors(); eigenVals.Print(); eigenVect.Print(); TDirectory *outnominal = (TDirectory *)fOutputRootFile->mkdir("nominal"); outnominal->cd(); // double *valst = FitUtils::GetArrayFromMap(fParams, fCurVals); // double chi2 = fSampleFCN->DoEval( valst ); // delete valst; fSampleFCN->Write(); // Loop over all throws TDirectory *throwsdir = (TDirectory *)fOutputRootFile->mkdir("throws"); throwsdir->cd(); int count = 0; // Produce all error throws. for (int i = 0; i < eigenVect.GetNrows(); i++) { TDirectory *throwfolder = (TDirectory *)throwsdir->mkdir(Form("throw_%i", count)); throwfolder->cd(); // Get New Parameter Vector NUIS_LOG(FIT, "Parameter Set " << count); for (int j = 0; j < eigenVect.GetNrows(); j++) { std::string param = fParams[j]; NUIS_LOG(FIT, " " << j << ". " << param << " : " << fCurVals[param] + sqrt(eigenVals[i]) * eigenVect[j][i]); fThrownVals[param] = fCurVals[param] + sqrt(eigenVals[i]) * eigenVect[j][i]; } // Run Eval double *vals = FitUtils::GetArrayFromMap(fParams, fThrownVals); double chi2 = fSampleFCN->DoEval(vals); NUIS_LOG(DEB, "Chi2 = " << chi2); delete vals; count++; fSampleFCN->Write(); throwfolder = (TDirectory *)throwsdir->mkdir(Form("throw_%i", count)); throwfolder->cd(); // Get New Parameter Vector NUIS_LOG(FIT, "Parameter Set " << count); for (int j = 0; j < eigenVect.GetNrows(); j++) { std::string param = fParams[j]; NUIS_LOG(FIT, " " << j << ". " << param << " : " << fCurVals[param] - sqrt(eigenVals[i]) * eigenVect[j][i]); fThrownVals[param] = fCurVals[param] - sqrt(eigenVals[i]) * eigenVect[j][i]; } // Run Eval double *vals2 = FitUtils::GetArrayFromMap(fParams, fThrownVals); chi2 = fSampleFCN->DoEval(vals2); NUIS_LOG(DEB, "Chi2 = " << chi2); delete vals2; count++; // Save the FCN fSampleFCN->Write(); } fOutputRootFile->Close(); fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "UPDATE"); fOutputRootFile->cd(); throwsdir = (TDirectory *)fOutputRootFile->Get("throws"); outnominal = (TDirectory *)fOutputRootFile->Get("nominal"); // Loop through Error DIR TDirectory *outerr = (TDirectory *)fOutputRootFile->mkdir("errors"); outerr->cd(); TIter next(outnominal->GetListOfKeys()); TKey *key; while ((key = (TKey *)next())) { TClass *cl = gROOT->GetClass(key->GetClassName()); if (!cl->InheritsFrom("TH1D") and !cl->InheritsFrom("TH2D")) continue; NUIS_LOG(FIT, "Creating error bands for " << key->GetName()); std::string plotname = std::string(key->GetName()); if (plotname.find("_EVT") != std::string::npos) continue; if (plotname.find("_FLUX") != std::string::npos) continue; if (plotname.find("_FLX") != std::string::npos) continue; TH1 *baseplot = (TH1D *)key->ReadObj()->Clone(Form("%s_ORIGINAL", key->GetName())); TH1 *errorplot_upper = (TH1D *)baseplot->Clone(Form("%s_ERROR_UPPER", key->GetName())); TH1 *errorplot_lower = (TH1D *)baseplot->Clone(Form("%s_ERROR_LOWER", key->GetName())); TH1 *meanplot = (TH1D *)baseplot->Clone(Form("%s_SET_MEAN", key->GetName())); TH1 *systplot = (TH1D *)baseplot->Clone(Form("%s_SYST", key->GetName())); TH1 *statplot = (TH1D *)baseplot->Clone(Form("%s_STAT", key->GetName())); TH1 *totlplot = (TH1D *)baseplot->Clone(Form("%s_TOTAL", key->GetName())); int nbins = 0; if (cl->InheritsFrom("TH1D")) nbins = ((TH1D *)baseplot)->GetNbinsX(); else nbins = ((TH1D *)baseplot)->GetNbinsX() * ((TH1D *)baseplot)->GetNbinsY(); meanplot->Reset(); errorplot_upper->Reset(); errorplot_lower->Reset(); for (int j = 0; j < nbins; j++) { errorplot_upper->SetBinError(j + 1, 0.0); errorplot_lower->SetBinError(j + 1, 0.0); } // Loop over throws and calculate mean and error for +- throws int addcount = 0; // Add baseplot first to slightly bias to central value meanplot->Add(baseplot); addcount++; for (int i = 0; i < count; i++) { TH1 *newplot = (TH1D *)throwsdir->Get(Form("throw_%i/%s", i, plotname.c_str())); if (!newplot) { NUIS_ERR(WRN, "Cannot find new plot : " << Form("throw_%i/%s", i, plotname.c_str())); NUIS_ERR(WRN, "This plot will not have the correct errors!"); continue; } newplot->SetDirectory(0); nbins = newplot->GetNbinsX(); for (int j = 0; j < nbins; j++) { if (i % 2 == 0) { errorplot_upper->SetBinContent( j + 1, errorplot_upper->GetBinContent(j + 1) + pow(baseplot->GetBinContent(j + 1) - newplot->GetBinContent(j + 1), 2)); } else { errorplot_lower->SetBinContent( j + 1, errorplot_lower->GetBinContent(j + 1) + pow(baseplot->GetBinContent(j + 1) - newplot->GetBinContent(j + 1), 2)); } meanplot->SetBinContent(j + 1, meanplot->GetBinContent(j + 1) + baseplot->GetBinContent(j + 1)); } delete newplot; addcount++; } // Get mean Average for (int j = 0; j < nbins; j++) { meanplot->SetBinContent(j + 1, meanplot->GetBinContent(j + 1) / double(addcount)); } for (int j = 0; j < nbins; j++) { errorplot_upper->SetBinContent( j + 1, sqrt(errorplot_upper->GetBinContent(j + 1))); errorplot_lower->SetBinContent( j + 1, sqrt(errorplot_lower->GetBinContent(j + 1))); statplot->SetBinError(j + 1, baseplot->GetBinError(j + 1)); systplot->SetBinError(j + 1, (errorplot_upper->GetBinContent(j + 1) + errorplot_lower->GetBinContent(j + 1)) / 2.0); totlplot->SetBinError(j + 1, sqrt(pow(statplot->GetBinError(j + 1), 2) + pow(systplot->GetBinError(j + 1), 2))); meanplot->SetBinError(j + 1, sqrt(pow(statplot->GetBinError(j + 1), 2) + pow(systplot->GetBinError(j + 1), 2))); } outerr->cd(); errorplot_upper->Write(); errorplot_lower->Write(); baseplot->Write(); meanplot->Write(); statplot->Write(); systplot->Write(); totlplot->Write(); delete errorplot_upper; delete errorplot_lower; delete baseplot; delete meanplot; delete statplot; delete systplot; delete totlplot; } } diff --git a/src/T2K/CMakeLists.txt b/src/T2K/CMakeLists.txt index e41ad31..e223e00 100644 --- a/src/T2K/CMakeLists.txt +++ b/src/T2K/CMakeLists.txt @@ -1,122 +1,126 @@ # Copyright 2016-2021 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 . ################################################################################ set(IMPLFILES T2K_CC0pi_XSec_H2O_2DPcos_anu.cxx T2K_NuMu_CC0pi_OC_XSec_2DPcos.cxx T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint.cxx T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos.cxx T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint.cxx T2K_CC0pi_XSec_2DPcos_nu_I.cxx T2K_CC0pi_XSec_2DPcos_nu_II.cxx T2K_CCinc_XSec_2DPcos_nu_nonuniform.cxx T2K_nueCCinc_XSec_1Dpe.cxx T2K_nueCCinc_XSec_1Dthe.cxx T2K_nueCCinc_XSec_1Dpe_joint.cxx T2K_nueCCinc_XSec_1Dthe_joint.cxx T2K_nueCCinc_XSec_joint.cxx T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.cxx T2K_CC1pip_CH_XSec_1Dppi_nu.cxx T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx T2K_CC1pip_CH_XSec_1DQ2_nu.cxx T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.cxx T2K_CC1pip_CH_XSec_1DCosThAdler_nu.cxx +T2K_CCCOH_C12_XSec_1DEnu_nu.cxx + T2K_CC1pip_H2O_XSec_1DEnuDelta_nu.cxx T2K_CC1pip_H2O_XSec_1DEnuMB_nu.cxx T2K_CC1pip_H2O_XSec_1Dcosmu_nu.cxx T2K_CC1pip_H2O_XSec_1Dcosmupi_nu.cxx T2K_CC1pip_H2O_XSec_1Dcospi_nu.cxx T2K_CC1pip_H2O_XSec_1Dpmu_nu.cxx T2K_CC1pip_H2O_XSec_1Dppi_nu.cxx T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx T2K_CC0piWithProtons_XSec_2018_multidif_0p_1p_Np.cxx T2K_CC0pinp_ifk_XSec_3Dinfp_nu.cxx T2K_CC0pinp_ifk_XSec_3Dinfa_nu.cxx T2K_CC0pinp_ifk_XSec_3Dinfip_nu.cxx T2K_SignalDef.cxx ) set(HEADERFILES T2K_CC0pi_XSec_H2O_2DPcos_anu.h T2K_NuMu_CC0pi_OC_XSec_2DPcos.h T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint.h T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos.h T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint.h T2K_CC0pi_XSec_2DPcos_nu_I.h T2K_CC0pi_XSec_2DPcos_nu_II.h T2K_CCinc_XSec_2DPcos_nu_nonuniform.h T2K_nueCCinc_XSec_1Dpe.h T2K_nueCCinc_XSec_1Dthe.h T2K_nueCCinc_XSec_1Dpe_joint.h T2K_nueCCinc_XSec_1Dthe_joint.h T2K_nueCCinc_XSec_joint.h T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.h T2K_CC1pip_CH_XSec_1Dppi_nu.h T2K_CC1pip_CH_XSec_1Dthpi_nu.h T2K_CC1pip_CH_XSec_1Dthmupi_nu.h T2K_CC1pip_CH_XSec_1DQ2_nu.h T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.h T2K_CC1pip_CH_XSec_1DCosThAdler_nu.h +T2K_CCCOH_C12_XSec_1DEnu_nu.h + T2K_CC1pip_H2O_XSec_1DEnuDelta_nu.h T2K_CC1pip_H2O_XSec_1DEnuMB_nu.h T2K_CC1pip_H2O_XSec_1Dcosmu_nu.h T2K_CC1pip_H2O_XSec_1Dcosmupi_nu.h T2K_CC1pip_H2O_XSec_1Dcospi_nu.h T2K_CC1pip_H2O_XSec_1Dpmu_nu.h T2K_CC1pip_H2O_XSec_1Dppi_nu.h T2K_CC0pinp_STV_XSec_1Ddpt_nu.h T2K_CC0pinp_STV_XSec_1Ddphit_nu.h T2K_CC0pinp_STV_XSec_1Ddat_nu.h T2K_CC0piWithProtons_XSec_2018_multidif_0p_1p_Np.h T2K_CC0pinp_ifk_XSec_3Dinfp_nu.h T2K_CC0pinp_ifk_XSec_3Dinfa_nu.h T2K_CC0pinp_ifk_XSec_3Dinfip_nu.h T2K_SignalDef.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(${MINIMUM_INCLUDE_DIRECTORIES}) set_target_properties(${LIBNAME} PROPERTIES VERSION "${NUISANCE_VERSION_MAJOR}.${NUISANCE_VERSION_MINOR}.${NUISANCE_VERSION_REVISION}") #set_target_properties(${LIBNAME} PROPERTIES LINK_FLAGS ${ROOT_LD_FLAGS}) if(DEFINED PROJECTWIDE_EXTRA_DEPENDENCIES) add_dependencies(${LIBNAME} ${PROJECTWIDE_EXTRA_DEPENDENCIES}) endif() install(TARGETS ${LIBNAME} DESTINATION lib) #Can uncomment this to install the headers... but is it really neccessary? install(FILES ${HEADERFILES} DESTINATION include/T2K) set(MODULETargets ${MODULETargets} ${LIBNAME} PARENT_SCOPE) diff --git a/src/T2K/T2K_CCCOH_C12_XSec_1DEnu_nu.cxx b/src/T2K/T2K_CCCOH_C12_XSec_1DEnu_nu.cxx new file mode 100644 index 0000000..8b78ac0 --- /dev/null +++ b/src/T2K/T2K_CCCOH_C12_XSec_1DEnu_nu.cxx @@ -0,0 +1,67 @@ +#include "T2K_CCCOH_C12_XSec_1DEnu_nu.h" + +// The constructor +T2K_CCCOH_C12_XSec_1DEnu_nu::T2K_CCCOH_C12_XSec_1DEnu_nu(nuiskey samplekey) { + + // Sample overview --------------------------------------------------- + std::string descrip = "T2K_CCCOH_C12_XSec_nu sample. \n" + "Target: C12 \n" + "Flux: T2K FHC numu \n" + "Signal: CCCOHpi+, p_mu,pi > 0.18 GeV, theta_mu,pi < 70 deg\n" + ", p_pi < 1.6 GeV\n" + "https://arxiv.org/abs/1604.04406"; + + // Setup common settings + fSettings = LoadSampleSettings(samplekey); + fSettings.SetTitle("T2K_CCCOH_C12_XSec_1DEnu_nu"); + fSettings.SetDescription(descrip); + fSettings.SetXTitle("E_{#nu} (GeV)"); + fSettings.SetYTitle("d#sigma/dE_{#nu} (cm^{2}/^{12}C)"); + fSettings.SetEnuRange(0.0, 100.0); + fSettings.DefineAllowedTargets("C"); + fSettings.DefineAllowedSpecies("numu"); + FinaliseSampleSettings(); + + // Override the default + fMCFine = new TH1D("mcfine", "mcfine", 50, 0, 5); + + // Scaling Setup --------------------------------------------------- + // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon + fScaleFactor = 12*(GetEventHistogram()->Integral("width") * 1E-38) / + double(fNEvents) / TotalIntegratedFlux("width"); + + // Plot Setup ------------------------------------------------------- + SetDataFromTextFile(GeneralUtils::GetTopLevelDir()+"/data/T2K/CCCOH/C12_Enu_1bin.txt"); + + // Slight hack to get single bin comparisons recognised + this->SetFitOptions("SINGLEBIN"); + + FinaliseMeasurement(); +}; + +void T2K_CCCOH_C12_XSec_1DEnu_nu::FillEventVariables(FitEvent *event) { + TLorentzVector Pnu = event->GetNeutrinoIn()->fP; + fXVar = Pnu.E()/1000.; + return; +}; + + +bool T2K_CCCOH_C12_XSec_1DEnu_nu::isSignal(FitEvent *event) { + + if (!SignalDef::isCCCOH(event, 14, 211)) + return false; + + TLorentzVector Pnu = event->GetHMISParticle(14)->fP; + TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; + TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; + + double p_mu = FitUtils::p(Pmu); + double p_pi = FitUtils::p(Ppip); + double th_mu = FitUtils::th(Pnu, Pmu)* 180 / TMath::Pi(); + double th_pi = FitUtils::th(Pnu, Ppip)* 180 / TMath::Pi(); + + if (p_mu < 0.18 || p_pi < 0.18 || p_pi > 1.6 || th_mu > 70 || th_pi > 70) { + return false; + } + return true; +} diff --git a/src/T2K/T2K_CCCOH_C12_XSec_1DEnu_nu.h b/src/T2K/T2K_CCCOH_C12_XSec_1DEnu_nu.h new file mode 100644 index 0000000..8d4d941 --- /dev/null +++ b/src/T2K/T2K_CCCOH_C12_XSec_1DEnu_nu.h @@ -0,0 +1,16 @@ +#ifndef T2K_CCCOH_C12_XSEC_1DENU_NU_H_SEEN +#define T2K_CCCOH_C12_XSEC_1DENU_NU_H_SEEN + +#include "Measurement1D.h" + +class T2K_CCCOH_C12_XSec_1DEnu_nu : public Measurement1D { +public: + T2K_CCCOH_C12_XSec_1DEnu_nu(nuiskey samplekey); + virtual ~T2K_CCCOH_C12_XSec_1DEnu_nu() {}; + + void FillEventVariables(FitEvent *event); + bool isSignal(FitEvent *event); + +}; + +#endif diff --git a/src/Utils/SignalDef.cxx b/src/Utils/SignalDef.cxx index 03456ad..13f4204 100644 --- a/src/Utils/SignalDef.cxx +++ b/src/Utils/SignalDef.cxx @@ -1,289 +1,291 @@ // Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #include "FitUtils.h" #include "SignalDef.h" bool SignalDef::isCCINC(FitEvent *event, int nuPDG, double EnuMin, double EnuMax) { // Check for the desired PDG code if (!event->HasISParticle(nuPDG)) return false; // Check that it's within the allowed range if set if (EnuMin != EnuMax) { if (!SignalDef::IsEnuInRange(event, EnuMin*1000, EnuMax*1000)) { return false; } } // Check that the charged lepton we expect has been produced if (!event->HasFSParticle(nuPDG > 0 ? nuPDG-1 : nuPDG+1)) return false; return true; } bool SignalDef::isNCINC(FitEvent *event, int nuPDG, double EnuMin, double EnuMax) { // Check for the desired PDG code before and after the interaction if (!event->HasISParticle(nuPDG) || !event->HasFSParticle(nuPDG)) return false; // Check that it's within the allowed range if set if (EnuMin != EnuMax) if (!SignalDef::IsEnuInRange(event, EnuMin*1000, EnuMax*1000)) return false; return true; } bool SignalDef::isCC0pi(FitEvent *event, int nuPDG, double EnuMin, double EnuMax){ // Check it's CCINC if (!SignalDef::isCCINC(event, nuPDG, EnuMin, EnuMax)) return false; // Veto event with mesons if (event->NumFSMesons() != 0) return false; // Veto events which don't have exactly 1 outgoing charged lepton if (event->NumFSLeptons() != 1) return false; return true; } bool SignalDef::isNC0pi(FitEvent *event, int nuPDG, double EnuMin, double EnuMax){ // Check it's NCINC if (!SignalDef::isNCINC(event, nuPDG, EnuMin, EnuMax)) return false; // Veto event with mesons if (event->NumFSMesons() != 0) return false; // Veto events with a charged lepton if (event->NumFSLeptons() != 0) return false; return true; } bool SignalDef::isCCQE(FitEvent *event, int nuPDG, double EnuMin, double EnuMax){ // Check if it's CCINC if (!SignalDef::isCCINC(event, nuPDG, EnuMin, EnuMax)) return false; // Require modes 1 (CCQE) if (abs(event->Mode) != 1) return false; return true; } bool SignalDef::isNCEL(FitEvent *event, int nuPDG, double EnuMin, double EnuMax){ // Check if it's NCINC if (!SignalDef::isNCINC(event, nuPDG, EnuMin, EnuMax)) return false; // Require modes 51/52 (NCEL) if (abs(event->Mode) != 51 && abs(event->Mode) != 52) return false; return true; } bool SignalDef::isCCQELike(FitEvent *event, int nuPDG, double EnuMin, double EnuMax){ // Check if it's CCINC if (!SignalDef::isCCINC(event, nuPDG, EnuMin, EnuMax)) return false; // Require modes 1/2 (CCQE and MEC) if (abs(event->Mode) != 1 && abs(event->Mode) != 2) return false; return true; } // Require one meson, one charged lepton. types specified in the arguments bool SignalDef::isCC1pi(FitEvent *event, int nuPDG, int piPDG, double EnuMin, double EnuMax){ // First, make sure it's CCINC if (!SignalDef::isCCINC(event, nuPDG, EnuMin, EnuMax)) return false; int nMesons = event->NumFSMesons(); int nLeptons = event->NumFSLeptons(); int nPion = event->NumFSParticle(piPDG); // Check that the desired pion exists and is the only meson if (nPion != 1 || nMesons != 1) return false; // Check that there is only one final state lepton if (nLeptons != 1) return false; return true; } // Require one meson, one neutrino. Types specified as arguments bool SignalDef::isNC1pi(FitEvent *event, int nuPDG, int piPDG, double EnuMin, double EnuMax){ // First, make sure it's NCINC if (!SignalDef::isNCINC(event, nuPDG, EnuMin, EnuMax)) return false; int nMesons = event->NumFSMesons(); int nLeptons = event->NumFSLeptons(); int nPion = event->NumFSParticle(piPDG); // Check that the desired pion exists and is the only meson if (nPion != 1 || nMesons != 1) return false; // Check that there are no charged leptons if (nLeptons != 0) return false; return true; } // A slightly ugly function to replace the BC 2pi channels. // All particles which are allowed in the final state are specified bool SignalDef::isCCWithFS(FitEvent *event, int nuPDG, std::vector pdgs, double EnuMin, double EnuMax){ // Check it's CCINC if (!SignalDef::isCCINC(event, nuPDG, EnuMin, EnuMax)) return false; // Remove events where the number of final state particles // do not match the number specified in the signal definition if ((int)pdgs.size() != event->NumFSParticle()) return false; // For every particle in the list, check the number in the FS for (std::vector::iterator it = pdgs.begin(); it != pdgs.end(); ++it){ // Check how many times this pdg is in the vector int nEntries = std::count (pdgs.begin(), pdgs.end(), *it); if (event->NumFSParticle(*it) != nEntries) return false; } return true; } // Require one meson, one charged lepton, AND specify the only other final state particle // This is only suitable for bubble chambers. Types specified in the arguments bool SignalDef::isCC1pi3Prong(FitEvent *event, int nuPDG, int piPDG, int thirdPDG, double EnuMin, double EnuMax){ // First, make sure it's CC1pi if (!SignalDef::isCC1pi(event, nuPDG, piPDG, EnuMin, EnuMax)) return false; // Check we have the third prong if (event->NumFSParticle(thirdPDG) == 0) return false; // Check that there are only three FS particles if (event->NumFSParticle() != 3) return false; return true; } // Require one meson, one neutrino, AND specify the only other final state particle // This is only suitable for bubble chambers. Types specified in the arguments bool SignalDef::isNC1pi3Prong(FitEvent *event, int nuPDG, int piPDG, int thirdPDG, double EnuMin, double EnuMax){ // First, make sure it's NC1pi if (!SignalDef::isNC1pi(event, nuPDG, piPDG, EnuMin, EnuMax)) return false; // Check we have the third prong if (event->NumFSParticle(thirdPDG) == 0) return false; // Check that there are only three FS particles if (event->NumFSParticle() != 3) return false; return true; } bool SignalDef::isCCCOH(FitEvent *event, int nuPDG, int piPDG, double EnuMin, double EnuMax){ // Check this is a CCINC event if (!SignalDef::isCCINC(event, nuPDG, EnuMin, EnuMax)) return false; int nLepton = event->NumFSParticle(nuPDG > 0 ? nuPDG-1 : nuPDG+1); int nPion = event->NumFSParticle(piPDG); int nFS = event->NumFSParticle(); if (nLepton != 1 || nPion != 1) return false; - if (nFS != 2) return false; + // if (nFS != 2) return false; + // GENIE v3 includes the nucleus in the final state, so modify the definition for now... + if (abs(event->Mode) != 16) return false; return true; } bool SignalDef::isNCCOH(FitEvent *event, int nuPDG, int piPDG, double EnuMin, double EnuMax){ // Check this is an NCINC event if (!SignalDef::isNCINC(event, nuPDG, EnuMin, EnuMax)) return false; int nLepton = event->NumFSParticle(nuPDG); int nPion = event->NumFSParticle(piPDG); int nFS = event->NumFSParticle(); if (nLepton != 1 || nPion != 1) return false; if (nFS != 2) return false; return true; } bool SignalDef::HasProtonKEAboveThreshold(FitEvent* event, double threshold){ for (uint i = 0; i < event->Npart(); i++){ FitParticle* p = event->PartInfo(i); if (!p->IsFinalState()) continue; if (p->fPID != 2212) continue; if (FitUtils::T(p->fP) > threshold / 1000.0) return true; } return false; } bool SignalDef::HasProtonMomAboveThreshold(FitEvent* event, double threshold){ for (uint i = 0; i < event->Npart(); i++){ FitParticle* p = event->PartInfo(i); if (!p->IsFinalState()) continue; if (p->fPID != 2212) continue; if (p->fP.Vect().Mag() > threshold) return true; } return false; } // Calculate the angle between the neutrino and an outgoing particle, apply a cut bool SignalDef::IsRestrictedAngle(FitEvent* event, int nuPDG, int otherPDG, double angle){ // If the particles don't exist, return false if (!event->HasISParticle(nuPDG) || !event->HasFSParticle(otherPDG)) return false; // Get Mom TVector3 pnu = event->GetHMISParticle(nuPDG)->fP.Vect(); TVector3 p2 = event->GetHMFSParticle(otherPDG)->fP.Vect(); double theta = pnu.Angle(p2) * 180. / TMath::Pi(); return (theta < angle); } bool SignalDef::IsEnuInRange(FitEvent* event, double emin, double emax){ return (event->GetNeutrinoIn()->fP.E() > emin && event->GetNeutrinoIn()->fP.E() < emax); }