diff --git a/app/PrepareGENIE.cxx b/app/PrepareGENIE.cxx index 939f87b..442326b 100644 --- a/app/PrepareGENIE.cxx +++ b/app/PrepareGENIE.cxx @@ -1,554 +1,559 @@ #include #include #include "FitLogger.h" #include "PlotUtils.h" #include "TFile.h" #include "TH1D.h" #include "TTree.h" #ifdef __GENIE_ENABLED__ #include "Conventions/Units.h" #include "GHEP/GHepParticle.h" #include "PDG/PDGUtils.h" #endif std::string gInputFiles = ""; std::string gOutputFile = ""; std::string gFluxFile = ""; std::string gTarget = ""; double MonoEnergy; bool IsMonoE = false; void PrintOptions(); void ParseOptions(int argc, char* argv[]); void RunGENIEPrepareMono(std::string input, std::string target, std::string output); void RunGENIEPrepare(std::string input, std::string flux, std::string target, std::string output); int main(int argc, char* argv[]) { ParseOptions(argc, argv); if (IsMonoE) { RunGENIEPrepareMono(gInputFiles, gTarget, gOutputFile); } else { RunGENIEPrepare(gInputFiles, gFluxFile, gTarget, gOutputFile); } } void RunGENIEPrepareMono(std::string input, std::string target, std::string output) { // Setup TTree TChain* tn = new TChain("gtree"); tn->AddFile(input.c_str()); int nevt = tn->GetEntries(); NtpMCEventRecord* genientpl = NULL; tn->SetBranchAddress("gmcrec", &genientpl); - TH1D *fluxhist = new TH1D("flux", "flux", 1000, 0, 10); + TH1D* fluxhist = new TH1D("flux", "flux", 1000, 0, 10); fluxhist->Fill(MonoEnergy); fluxhist->Scale(1, "width"); // Make Event Hist TH1D* eventhist = (TH1D*)fluxhist->Clone(); eventhist->Reset(); TH1D* xsechist = (TH1D*)eventhist->Clone(); // Create maps std::map modexsec; std::map modecount; std::vector genieids; std::vector targetids; std::vector interids; // Loop over all events for (int i = 0; i < nevt; i++) { tn->GetEntry(i); StopTalking(); EventRecord& event = *(genientpl->event); GHepParticle* neu = event.Probe(); StartTalking(); // Get XSec From Spline GHepRecord genie_record = static_cast(event); double xsec = (genie_record.XSec() / (1E-38 * genie::units::cm2)); // Parse Interaction String std::string mode = genie_record.Summary()->AsString(); std::vector modevec = GeneralUtils::ParseToStr(mode, ";"); std::string targ = (modevec[0] + ";" + modevec[1]); std::string inter = mode; // Fill lists of Unique IDS if (std::find(targetids.begin(), targetids.end(), targ) == targetids.end()) { targetids.push_back(targ); } if (std::find(interids.begin(), interids.end(), inter) == interids.end()) { interids.push_back(inter); } // Create entries Mode Maps if (modexsec.find(mode) == modexsec.end()) { genieids.push_back(mode); modexsec[mode] = (TH1D*)xsechist->Clone(); modecount[mode] = (TH1D*)xsechist->Clone(); } // Fill XSec Histograms modexsec[mode]->Fill(neu->E(), xsec); modecount[mode]->Fill(neu->E()); // Fill total event hist eventhist->Fill(neu->E()); // Clear Event genientpl->Clear(); if (i % (nevt / 20) == 0) { LOG(FIT) << "Processed " << i << "/" << nevt << " GENIE events." << std::endl; } } LOG(FIT) << "Processed all events" << std::endl; TFile* outputfile = new TFile(input.c_str(), "UPDATE"); outputfile->cd(); LOG(FIT) << "Getting splines " << std::endl; // Save each of the reconstructed splines to file std::map modeavg; TDirectory* inddir = (TDirectory*)outputfile->Get("IndividualGENIESplines"); if (!inddir) inddir = (TDirectory*)outputfile->mkdir("IndividualGENIESplines"); inddir->cd(); // Loop over GENIE ID's and get MEC count int MECcount = 0; bool MECcorrect = FitPar::Config().GetParB("CorrectGENIEMECNorm"); for (UInt_t i = 0; i < genieids.size(); i++) { if (genieids[i].find("MEC") != std::string::npos) { MECcount++; } } LOG(FIT) << "Found " << MECcount << " repeated MEC instances." << std::endl; for (UInt_t i = 0; i < genieids.size(); i++) { std::string mode = genieids[i]; modexsec[mode]->Write((mode + "_summed_xsec").c_str(), TObject::kOverwrite); modecount[mode]->Write((mode + "_summed_evt").c_str(), TObject::kOverwrite); // Form extra avg xsec map -> Reconstructed spline modeavg[mode] = (TH1D*)modexsec[mode]->Clone(); modeavg[mode]->Divide(modecount[mode]); if (MECcorrect && (mode.find("MEC") != std::string::npos)) { modeavg[mode]->Scale(1.0 / double(MECcount)); } modeavg[mode]->Write((mode + "_rec_spline").c_str(), TObject::kOverwrite); } TDirectory* targdir = (TDirectory*)outputfile->Get("TargetGENIESplines"); if (!targdir) targdir = (TDirectory*)outputfile->mkdir("TargetGENIESplines"); targdir->cd(); LOG(FIT) << "Getting Target Splines" << std::endl; // For each target save a total spline std::map targetsplines; for (uint i = 0; i < targetids.size(); i++) { LOG(FIT) << "Getting target " << i << std::endl; std::string targ = targetids[i]; targetsplines[targ] = (TH1D*)xsechist->Clone(); LOG(FIT) << "Created target spline for " << targ << std::endl; for (uint j = 0; j < genieids.size(); j++) { std::string mode = genieids[j]; if (mode.find(targ) != std::string::npos) { LOG(FIT) << "Mode " << mode << " contains " << targ << " target!" << std::endl; targetsplines[targ]->Add(modeavg[mode]); LOG(FIT) << "Finished with Mode " << mode << " " << modeavg[mode]->Integral() << std::endl; } } LOG(FIT) << "Saving target spline:" << targ << std::endl; targetsplines[targ]->Write(("Total" + targ).c_str(), TObject::kOverwrite); } LOG(FIT) << "Getting total splines" << std::endl; // Now we have each of the targets we need to create a total cross-section. int totalnucl = 0; std::vector targprs = GeneralUtils::ParseToStr(target, ","); TH1D* totalxsec = (TH1D*)xsechist->Clone(); for (uint i = 0; i < targprs.size(); i++) { std::string targpdg = targprs[i]; for (std::map::iterator iter = targetsplines.begin(); iter != targetsplines.end(); iter++) { std::string targstr = iter->first; TH1D* xsec = iter->second; if (targstr.find(targpdg) != std::string::npos) { LOG(FIT) << "Adding target spline " << targstr << " Integral = " << xsec->Integral("width") << std::endl; totalxsec->Add(xsec); int nucl = atoi(targpdg.c_str()); totalnucl += int((nucl % 10000) / 10); } } } outputfile->cd(); totalxsec->Write("nuisance_Xsec", TObject::kOverwrite); eventhist = (TH1D*)totalxsec->Clone(); eventhist->Multiply(fluxhist); eventhist->Write("nuisance_events", TObject::kOverwrite); fluxhist->Write("nuisance_flux", TObject::kOverwrite); LOG(FIT) << "Inclusive XSec Per Nucleon = " << eventhist->Integral("width") * 1E-38 / fluxhist->Integral("width") << std::endl; std::cout << "XSec Hist Integral = " << xsechist->Integral("width") << std::endl; return; } void RunGENIEPrepare(std::string input, std::string flux, std::string target, std::string output) { LOG(FIT) << "Running GENIE Prepare" << std::endl; // 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"); if (!fluxfile->IsZombie()) { - fluxhist = (TH1D*)fluxfile->Get(fluxvect[1].c_str()); + fluxhist = dynamic_cast(fluxfile->Get(fluxvect[1].c_str())); + if (!fluxhist) { + ERR(FTL) << "Couldn't find histogram named: \"" << fluxvect[1] + << "\" in file: \"" << fluxvect[0] << std::endl; + throw; + } fluxhist->SetDirectory(0); } else { // Function with EnuRange if (fluxvect.size() == 3) { ERR(FTL) << "FUNCTION WITH ENU RANGE NOT SUPPORTED SORRY!" << std::endl; throw; } } } else if (fluxvect.size() == 1) { MonoEnergy = GeneralUtils::StrToDbl(fluxvect[0]); RunGENIEPrepareMono(input, target, output); return; } else { LOG(FTL) << "NO FLUX SPECIFIED" << std::endl; throw; } // Setup TTree TChain* tn = new TChain("gtree"); tn->AddFile(input.c_str()); int nevt = tn->GetEntries(); NtpMCEventRecord* genientpl = NULL; tn->SetBranchAddress("gmcrec", &genientpl); // Make Event Hist TH1D* eventhist = (TH1D*)fluxhist->Clone(); eventhist->Reset(); TH1D* xsechist = (TH1D*)eventhist->Clone(); // Create maps std::map modexsec; std::map modecount; std::vector genieids; std::vector targetids; std::vector interids; // Loop over all events for (int i = 0; i < nevt; i++) { tn->GetEntry(i); StopTalking(); EventRecord& event = *(genientpl->event); GHepParticle* neu = event.Probe(); StartTalking(); // Get XSec From Spline GHepRecord genie_record = static_cast(event); double xsec = (genie_record.XSec() / (1E-38 * genie::units::cm2)); // Parse Interaction String std::string mode = genie_record.Summary()->AsString(); std::vector modevec = GeneralUtils::ParseToStr(mode, ";"); std::string targ = (modevec[0] + ";" + modevec[1]); std::string inter = mode; // Fill lists of Unique IDS if (std::find(targetids.begin(), targetids.end(), targ) == targetids.end()) { targetids.push_back(targ); } if (std::find(interids.begin(), interids.end(), inter) == interids.end()) { interids.push_back(inter); } // Create entries Mode Maps if (modexsec.find(mode) == modexsec.end()) { genieids.push_back(mode); modexsec[mode] = (TH1D*)xsechist->Clone(); modecount[mode] = (TH1D*)xsechist->Clone(); } // Fill XSec Histograms modexsec[mode]->Fill(neu->E(), xsec); modecount[mode]->Fill(neu->E()); // Fill total event hist eventhist->Fill(neu->E()); // Clear Event genientpl->Clear(); if (i % (nevt / 20) == 0) { LOG(FIT) << "Processed " << i << "/" << nevt << " GENIE events." << std::endl; } } LOG(FIT) << "Processed all events" << std::endl; // Once event loop is done we can start saving stuff into the file TFile* outputfile = new TFile(input.c_str(), "UPDATE"); outputfile->cd(); LOG(FIT) << "Getting splines " << std::endl; // Save each of the reconstructed splines to file std::map modeavg; TDirectory* inddir = (TDirectory*)outputfile->Get("IndividualGENIESplines"); if (!inddir) inddir = (TDirectory*)outputfile->mkdir("IndividualGENIESplines"); inddir->cd(); // Loop over GENIE ID's and get MEC count int MECcount = 0; bool MECcorrect = FitPar::Config().GetParB("CorrectGENIEMECNorm"); for (UInt_t i = 0; i < genieids.size(); i++) { if (genieids[i].find("MEC") != std::string::npos) { MECcount++; } } LOG(FIT) << "Found " << MECcount << " repeated MEC instances." << std::endl; for (UInt_t i = 0; i < genieids.size(); i++) { std::string mode = genieids[i]; modexsec[mode]->Write((mode + "_summed_xsec").c_str(), TObject::kOverwrite); modecount[mode]->Write((mode + "_summed_evt").c_str(), TObject::kOverwrite); // Form extra avg xsec map -> Reconstructed spline modeavg[mode] = (TH1D*)modexsec[mode]->Clone(); modeavg[mode]->Divide(modecount[mode]); if (MECcorrect && (mode.find("MEC") != std::string::npos)) { modeavg[mode]->Scale(1.0 / double(MECcount)); } modeavg[mode]->Write((mode + "_rec_spline").c_str(), TObject::kOverwrite); } TDirectory* targdir = (TDirectory*)outputfile->Get("TargetGENIESplines"); if (!targdir) targdir = (TDirectory*)outputfile->mkdir("TargetGENIESplines"); targdir->cd(); LOG(FIT) << "Getting Target Splines" << std::endl; // For each target save a total spline std::map targetsplines; for (uint i = 0; i < targetids.size(); i++) { LOG(FIT) << "Getting target " << i << std::endl; std::string targ = targetids[i]; targetsplines[targ] = (TH1D*)xsechist->Clone(); LOG(FIT) << "Created target spline for " << targ << std::endl; for (uint j = 0; j < genieids.size(); j++) { std::string mode = genieids[j]; // Look at all matching modes/targets if (mode.find(targ) != std::string::npos) { LOG(FIT) << "Mode " << mode << " contains " << targ << " target!" << std::endl; // modeavg[mode]->Write( (mode + "_cont_" + targ).c_str() , // TObject::kOverwrite); targetsplines[targ]->Add(modeavg[mode]); LOG(FIT) << "Finished with Mode " << mode << " " << modeavg[mode]->Integral() << std::endl; } } LOG(FIT) << "Saving target spline:" << targ << std::endl; targetsplines[targ]->Write(("Total" + targ).c_str(), TObject::kOverwrite); } LOG(FIT) << "Getting total splines" << std::endl; // Now we have each of the targets we need to create a total cross-section. int totalnucl = 0; std::vector targprs = GeneralUtils::ParseToStr(target, ","); TH1D* totalxsec = (TH1D*)xsechist->Clone(); for (uint i = 0; i < targprs.size(); i++) { std::string targpdg = targprs[i]; for (std::map::iterator iter = targetsplines.begin(); iter != targetsplines.end(); iter++) { std::string targstr = iter->first; TH1D* xsec = iter->second; if (targstr.find(targpdg) != std::string::npos) { LOG(FIT) << "Adding target spline " << targstr << " Integral = " << xsec->Integral("width") << std::endl; totalxsec->Add(xsec); int nucl = atoi(targpdg.c_str()); totalnucl += int((nucl % 10000) / 10); } } } LOG(FIT) << "Total XSec Integral = " << totalxsec->Integral("width") << std::endl; outputfile->cd(); totalxsec->Write("nuisance_xsec", TObject::kOverwrite); eventhist = (TH1D*)fluxhist->Clone(); eventhist->Multiply(totalxsec); LOG(FIT) << "Dividing by Total Nucl = " << totalnucl << std::endl; eventhist->Scale(1.0 / double(totalnucl)); eventhist->Write("nuisance_events", TObject::kOverwrite); fluxhist->Write("nuisance_flux", TObject::kOverwrite); LOG(FIT) << "Inclusive XSec Per Nucleon = " << eventhist->Integral("width") * 1E-38 / fluxhist->Integral("width") << std::endl; std::cout << "XSec Hist Integral = " << xsechist->Integral("width") << std::endl; return; }; void PrintOptions() { std::cout << "PrepareGENIEEvents NUISANCE app. " << std::endl << "Takes GHep Outputs and prepares events for NUISANCE." << std::endl << std::endl << "PrepareGENIEEvents [-h,-help,--h,--help] [-i " "inputfile1.root,inputfile2.root,inputfile3.root,...] " << "[-f flux_root_file.root,flux_hist_name] [-t " "target1[frac1],target2[frac2],...]" << std::endl << std::endl; std::cout << "Prepare Mode [Default] : Takes a single GHep file, " "reconstructs the original GENIE splines, " << " and creates a duplicate file that also contains the flux, " "event rate, and xsec predictions that NUISANCE needs. " << std::endl; std::cout << "Following options are required for Prepare Mode:" << std::endl; std::cout << " [ -i inputfile.root ] : Reads in a single GHep input file " "that needs the xsec calculation ran on it. " << std::endl; std::cout << " [ -f flux_file.root,hist_name ] : Path to root file " "containing the flux histogram the GHep records were generated " "with." << " A simple method is to point this to the flux histogram genie " "generatrs '-f /path/to/events/input-flux.root,spectrum'. " << std::endl; std::cout << " [ -t target ] : Target that GHepRecords were generated with. " "Comma seperated list. E.g. for CH2 " "target=1000060120,1000010010,1000010010" << std::endl; std::cout << " [ -o outputfile.root ] : File to write prepared input file to." << std::endl; std::cout << " [ -m Mono_E_nu_GeV ] : Run in mono-energetic mode." << 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")) { gInputFiles = argv[i + 1]; ++i; } else if (!std::strcmp(argv[i], "-o")) { gOutputFile = argv[i + 1]; ++i; } else if (!std::strcmp(argv[i], "-f")) { gFluxFile = argv[i + 1]; ++i; } else if (!std::strcmp(argv[i], "-t")) { gTarget = argv[i + 1]; ++i; } else if (!std::strcmp(argv[i], "-m")) { MonoEnergy = GeneralUtils::StrToDbl(argv[i + 1]); IsMonoE = true; ++i; } else { ERR(FTL) << "ERROR: unknown command line option given! - '" << argv[i] << " " << argv[i + 1] << "'" << std::endl; PrintOptions(); break; } } } if (gInputFiles == "" && !flagopt) { ERR(FTL) << "No input file(s) specified!" << std::endl; flagopt = true; } if (gFluxFile == "" && !flagopt) { ERR(FTL) << "No flux input specified for Prepare Mode" << std::endl; flagopt = true; } if (gTarget == "" && !flagopt) { ERR(FTL) << "No target specified for Prepare Mode" << std::endl; flagopt = true; } if (argc < 1 || flagopt) { PrintOptions(); exit(-1); } return; } diff --git a/cmake/c++CompilerSetup.cmake b/cmake/c++CompilerSetup.cmake index ebdd76b..89a51eb 100644 --- a/cmake/c++CompilerSetup.cmake +++ b/cmake/c++CompilerSetup.cmake @@ -1,103 +1,108 @@ # Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret ################################################################################ # This file is part of NUISANCE. # # NUISANCE is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # NUISANCE is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with NUISANCE. If not, see . ################################################################################ if(USE_OMP) LIST(APPEND EXTRA_CXX_FLAGS -fopenmp) endif() if(USE_DYNSAMPLES) LIST(APPEND EXTRA_LIBS dl) LIST(APPEND EXTRA_CXX_FLAGS -D__USE_DYNSAMPLES__) endif() set(CXX_WARNINGS -Wall ) cmessage(DEBUG "EXTRA_CXX_FLAGS: ${EXTRA_CXX_FLAGS}") string(REPLACE ";" " " STR_EXTRA_CXX_FLAGS "${EXTRA_CXX_FLAGS}") set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${STR_EXTRA_CXX_FLAGS} ${CXX_WARNINGS}") set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -O0") + +if(USE_DYNSAMPLES) + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fPIC") +endif() + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -fPIC -O3") if(CMAKE_BUILD_TYPE MATCHES DEBUG) set(CURRENT_CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS_DEBUG}) elseif(CMAKE_BUILD_TYPE MATCHES RELEASE) set(CURRENT_CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS_RELEASE}) else() cmessage(FATAL_ERROR "[ERROR]: Unknown CMAKE_BUILD_TYPE (\"${CMAKE_BUILD_TYPE}\"): Should be \"DEBUG\" or \"RELEASE\".") endif() SET(STR_EXTRA_LINK_DIRS) if(NOT EXTRA_LINK_DIRS STREQUAL "") string(REPLACE ";" " -L" STR_EXTRA_LINK_DIRS "-L${EXTRA_LINK_DIRS}") endif() SET(STR_EXTRA_LIBS) if(NOT EXTRA_LIBS STREQUAL "") string(REPLACE ";" " -l" STR_EXTRA_LIBS "-l${EXTRA_LIBS}") endif() SET(STR_EXTRA_SHAREDOBJS) if(NOT EXTRA_SHAREDOBJS STREQUAL "") string(REPLACE ";" " " STR_EXTRA_SHAREDOBJS "${EXTRA_SHAREDOBJS}") endif() SET(STR_EXTRA_LINK_FLAGS) if(NOT EXTRA_LINK_FLAGS STREQUAL "") string(REPLACE ";" " " STR_EXTRA_LINK_FLAGS "${EXTRA_LINK_FLAGS}") endif() cmessage(DEBUG "EXTRA_LINK_DIRS: ${STR_EXTRA_LINK_DIRS}") cmessage(DEBUG "EXTRA_LIBS: ${STR_EXTRA_LIBS}") cmessage(DEBUG "EXTRA_SHAREDOBJS: ${STR_EXTRA_SHAREDOBJS}") cmessage(DEBUG "EXTRA_LINK_FLAGS: ${STR_EXTRA_LINK_FLAGS}") if(NOT STR_EXTRA_LINK_DIRS STREQUAL "" AND NOT STR_EXTRA_LIBS STREQUAL "") SET(CMAKE_DEPENDLIB_FLAGS "${STR_EXTRA_LINK_DIRS} ${STR_EXTRA_LIBS}") endif() if(NOT EXTRA_SHAREDOBJS STREQUAL "") if(NOT STR_EXTRA_LINK_FLAGS STREQUAL "") SET(STR_EXTRA_LINK_FLAGS "${STR_EXTRA_SHAREDOBJS} ${STR_EXTRA_LINK_FLAGS}") else() SET(STR_EXTRA_LINK_FLAGS "${STR_EXTRA_SHAREDOBJS}") endif() endif() if(NOT EXTRA_LINK_FLAGS STREQUAL "") if(NOT CMAKE_LINK_FLAGS STREQUAL "") SET(CMAKE_LINK_FLAGS "${CMAKE_LINK_FLAGS} ${STR_EXTRA_LINK_FLAGS}") else() SET(CMAKE_LINK_FLAGS "${STR_EXTRA_LINK_FLAGS}") endif() endif() if(USE_OMP) cmessage(FATAL_ERROR "No OMP features currently enabled so this is a FATAL_ERROR to let you know that you don't gain anything with this declaration.") endif() if (VERBOSE) cmessage (STATUS "C++ Compiler : ${CXX_COMPILER_NAME}") cmessage (STATUS " flags : ${CMAKE_CXX_FLAGS}") cmessage (STATUS " Release flags : ${CMAKE_CXX_FLAGS_RELEASE}") cmessage (STATUS " Debug flags : ${CMAKE_CXX_FLAGS_DEBUG}") cmessage (STATUS " Link Flags : ${CMAKE_LINK_FLAGS}") cmessage (STATUS " Lib Flags : ${CMAKE_DEPENDLIB_FLAGS}") endif() diff --git a/src/FitBase/Measurement1D.cxx b/src/FitBase/Measurement1D.cxx index 0028d0f..cb2cdd3 100644 --- a/src/FitBase/Measurement1D.cxx +++ b/src/FitBase/Measurement1D.cxx @@ -1,1903 +1,1903 @@ // Copyright 2016 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; // Fake Data fFakeDataInput = ""; fFakeDataFile = NULL; // Options fDefaultTypes = "FIX/FULL/CHI2"; fAllowedTypes = "FIX,FREE,SHAPE/FULL,DIAG/CHI2/NORM/ENUCORR/Q2CORR/ENU1D/MASK/NOWIDTH"; fIsFix = false; fIsShape = false; fIsFree = false; fIsDiag = false; fIsFull = false; fAddNormPen = false; fIsMask = false; fIsChi2SVD = false; fIsRawEvents = false; fIsNoWidth = false; fIsDifXSec = false; fIsEnu1D = false; // Inputs fInput = NULL; fRW = NULL; // Extra Histograms fMCHist_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; } //******************************************************************** 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 LOG(SAM) << "Finalising Sample Settings: " << fName << std::endl; if ((fSettings.GetS("originalname").find("Evt") != std::string::npos)) { fIsRawEvents = true; LOG(SAM) << "Found event rate measurement but using poisson likelihoods." << std::endl; } if (fSettings.GetS("originalname").find("XSec_1DEnu") != std::string::npos) { fIsEnu1D = true; LOG(SAM) << "::" << fName << "::" << std::endl; LOG(SAM) << "Found XSec Enu measurement, applying flux integrated scaling, " << "not flux averaged!" << std::endl; } if (fIsEnu1D && fIsRawEvents) { LOG(SAM) << "Found 1D Enu XSec distribution AND fIsRawEvents, is this " "really correct?!" << std::endl; LOG(SAM) << "Check experiment constructor for " << fName << " and correct this!" << std::endl; LOG(SAM) << "I live in " << __FILE__ << ":" << __LINE__ << std::endl; exit(-1); } 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")); if (fAddNormPen) { if (fNormError <= 0.0) { ERR(WRN) << "Norm error for class " << fName << " is 0.0!" << std::endl; ERR(WRN) << "If you want to use it please add fNormError=VAL" << std::endl; throw; } } } //******************************************************************** 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) { //******************************************************************** LOG(SAM) << "Reading data from text file: " << datafile << std::endl; fDataHist = PlotUtils::GetTH1DFromFile(datafile, fSettings.GetName() + "_data", fSettings.GetFullTitles()); } //******************************************************************** void Measurement1D::SetDataFromRootFile(std::string datafile, std::string histname) { //******************************************************************** LOG(SAM) << "Reading data from root file: " << datafile << ";" << histname << std::endl; 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) { ERR(FTL) << "Need a data hist to setup possion errors! " << std::endl; ERR(FTL) << "Setup Data First!" << std::endl; 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) { LOG(SAM) << "Setting diagonal covariance for: " << data->GetName() << std::endl; fFullCovar = StatUtils::MakeDiagonalCovarMatrix(data); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); } else { ERR(FTL) << "No data input provided to set diagonal covar from!" << std::endl; } // if (!fIsDiag) { // ERR(FTL) << "SetCovarMatrixFromDiag called for measurement " // << "that is not set as diagonal." << std::endl; // throw; // } } //******************************************************************** void Measurement1D::SetCovarFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) { dim = fDataHist->GetNbinsX(); } LOG(SAM) << "Reading covariance from text file: " << covfile << std::endl; fFullCovar = StatUtils::GetCovarFromTextFile(covfile, dim); covar = StatUtils::GetInvert(fFullCovar); 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){ LOG(SAM) << "Reading covariance from text file: " << covList[i] << std::endl; TMatrixDSym* temp_cov = StatUtils::GetCovarFromTextFile(covList[i], dim); (*fFullCovar) += (*temp_cov); delete temp_cov; } covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement1D::SetCovarFromRootFile(std::string covfile, std::string histname) { //******************************************************************** LOG(SAM) << "Reading covariance from text file: " << covfile << ";" << histname << std::endl; fFullCovar = StatUtils::GetCovarFromRootFile(covfile, histname); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement1D::SetCovarInvertFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) { dim = fDataHist->GetNbinsX(); } LOG(SAM) << "Reading inverted covariance from text file: " << covfile << std::endl; covar = StatUtils::GetCovarFromTextFile(covfile, dim); fFullCovar = StatUtils::GetInvert(covar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement1D::SetCovarInvertFromRootFile(std::string covfile, std::string histname) { //******************************************************************** LOG(SAM) << "Reading inverted covariance from text file: " << covfile << ";" << histname << std::endl; covar = StatUtils::GetCovarFromRootFile(covfile, histname); fFullCovar = StatUtils::GetInvert(covar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement1D::SetCorrelationFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) dim = fDataHist->GetNbinsX(); LOG(SAM) << "Reading data correlations from text file: " << covfile << ";" << dim << std::endl; TMatrixDSym* correlation = StatUtils::GetCovarFromTextFile(covfile, dim); if (!fDataHist) { ERR(FTL) << "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" << std::endl; throw; } // 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); 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){ LOG(SAM) << "Reading covariance from text file: " << corrList[i] << std::endl; 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); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement1D::SetCorrelationFromRootFile(std::string covfile, std::string histname) { //******************************************************************** LOG(SAM) << "Reading data correlations from text file: " << covfile << ";" << histname << std::endl; TMatrixDSym* correlation = StatUtils::GetCovarFromRootFile(covfile, histname); if (!fDataHist) { ERR(FTL) << "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" << std::endl; throw; } // 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); fDecomp = StatUtils::GetDecomp(fFullCovar); delete correlation; } //******************************************************************** void Measurement1D::SetCholDecompFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) { dim = fDataHist->GetNbinsX(); } LOG(SAM) << "Reading cholesky from text file: " << covfile << std::endl; 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); fDecomp = StatUtils::GetDecomp(fFullCovar); delete temp; delete trans; } //******************************************************************** void Measurement1D::SetCholDecompFromRootFile(std::string covfile, std::string histname) { //******************************************************************** LOG(SAM) << "Reading cholesky decomp from root file: " << covfile << ";" << histname << std::endl; 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); 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; LOG(SAM) << "Reading bin mask from file: " << maskfile << std::endl; // 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(), ifstream::in); if (!mask.is_open()) { LOG(FTL) << " Cannot find mask file." << std::endl; throw; } 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. 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() { //******************************************************************** LOG(SAM) << "Finalising Measurement: " << fName << std::endl; 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) { ERR(FTL) << "No data has been setup inside " << fName << " constructor!" << std::endl; throw; } // Make sure covariances are setup if (!fFullCovar) { fIsDiag = true; SetCovarFromDiagonal(fDataHist); } if (!covar) { covar = StatUtils::GetInvert(fFullCovar); } 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 and FitPar::Config().GetParB("UseShapeCovar")){ + if (fIsShape && fShapeCovar and FitPar::Config().GetParB("UseShapeCovar")){ if (covar) delete covar; covar = StatUtils::GetInvert(fShapeCovar); if (fDecomp) delete fDecomp; fDecomp = StatUtils::GetDecomp(fFullCovar); } // 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)); fMCFine->SetNameTitle((fSettings.GetName() + "_MC_FINE").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); SetAutoProcessTH1(fMCHist_Modes, kCMD_Reset, kCMD_Norm, kCMD_Write); } // Setup bin masks using sample name if (fIsMask) { 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) { ERR(FTL) << "I found a negative fScaleFactor in " << __FILE__ << ":" << __LINE__ << std::endl; ERR(FTL) << "fScaleFactor = " << fScaleFactor << std::endl; ERR(FTL) << "EXITING" << std::endl; 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) { ERR(FTL) << "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." << std::endl; throw; } } } // 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) { ERR(FTL) << "ERROR: Fit Option '" << fit_options_input.at(i) << "' Provided is not allowed for this measurement." << std::endl; ERR(FTL) << "Fit Options should be provided as a '/' seperated list " "(e.g. FREE/DIAG/NORM)" << std::endl; ERR(FTL) << "Available options for " << fName << " are '" << fAllowedTypes << "'" << std::endl; throw; } } // 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; ERR(FTL) << "No other LIKELIHOODS properly supported!" << std::endl; ERR(FTL) << "Try to use a chi2!" << std::endl; 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; 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(), ifstream::in); // Note that the smearing matrix may be rectangular. fSmearMatrix = new TMatrixD(truedim, recodim); if (smear.is_open()) LOG(SAM) << "Reading smearing matrix from file: " << smearfile << std::endl; else ERR(FTL) << "Smearing matrix provided is incorrect: " << smearfile << std::endl; while (std::getline(smear >> 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) { ERR(WRN) << fName << ": attempted to apply smearing matrix, but none was set" << std::endl; 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) { fMCHist->Fill(fXVar, Weight); fMCFine->Fill(fXVar, Weight); fMCStat->Fill(fXVar, 1.0); if (fMCHist_Modes) fMCHist_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); // 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) { // PlotUtils::FluxUnfoldedScaling(fMCHist_Modes, GetFluxHistogram(), // GetEventHistogram(), fScaleFactor, // fNEvents); // } } else if (fIsNoWidth) { fMCHist->Scale(fScaleFactor); fMCFine->Scale(fScaleFactor); if (fMCHist_Modes) fMCHist_Modes->Scale(fScaleFactor); // Any other differential scaling } else { fMCHist->Scale(fScaleFactor, "width"); fMCFine->Scale(fScaleFactor, "width"); if (fMCHist_Modes) fMCHist_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) { 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); } } // Sort Penalty Terms if (fAddNormPen) { 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; LOG(SAM) << "Setting fake data from : " << fFakeDataInput << std::endl; // 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); if (fDecomp) delete fDecomp; fDecomp = StatUtils::GetInvert(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 GetDataList().at(0)->Write(); GetMCList().at(0)->Write(); - if(fEvtRateScaleFactor != 0xdeadbeef){ + 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 (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()); } if (drawOpt.find("INVCOV") != std::string::npos && covar) { PlotUtils::GetInvCovarPlot(covar, fSettings.GetName()); } if (drawOpt.find("DECOMP") != std::string::npos && fDecomp) { PlotUtils::GetDecompCovarPlot(fDecomp, fSettings.GetName()); } // // 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; // } // Write Extra Histograms AutoWriteExtraTH1(); WriteExtraHistograms(); // Returning LOG(SAM) << "Written Histograms: " << fName << std::endl; 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()); if (fShapeCovar) StatUtils::SetDataErrorFromCov(dataShape, fShapeCovar, 1E-38); 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; LOG(SAM) << "Found event rate measurement but fIsRawEvents == false!" << std::endl; LOG(SAM) << "Overriding this and setting fIsRawEvents == true!" << std::endl; } fIsEnu1D = false; if (fName.find("XSec_1DEnu") != std::string::npos) { fIsEnu1D = true; LOG(SAM) << "::" << fName << "::" << std::endl; LOG(SAM) << "Found XSec Enu measurement, applying flux integrated scaling, " "not flux averaged!" << std::endl; } if (fIsEnu1D && fIsRawEvents) { LOG(SAM) << "Found 1D Enu XSec distribution AND fIsRawEvents, is this " "really correct?!" << std::endl; LOG(SAM) << "Check experiment constructor for " << fName << " and correct this!" << std::endl; LOG(SAM) << "I live in " << __FILE__ << ":" << __LINE__ << std::endl; exit(-1); } 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); return; } //******************************************************************** void Measurement1D::SetDataValues(std::string dataFile) { //******************************************************************** // Override this function if the input file isn't in a suitable format LOG(SAM) << "Reading data from: " << dataFile.c_str() << std::endl; 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) { //******************************************************************** LOG(SAM) << "Filling histogram from " << inhistfile << "->" << histname << std::endl; 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) { //******************************************************************** LOG(SAM) << "Filling histogram from " << inhistfile << "->" << histname << std::endl; 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* decmpPlot = new TH2D(); TH2D* covarInvPlot = 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()); covarInvPlot = (TH2D*)tempFile->Get((covName + "covinv").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 ERR(WRN) << "Incorrect thrown_covariance option in parameters." << std::endl; 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(), ifstream::in); this->covar = new TMatrixDSym(dim); fFullCovar = new TMatrixDSym(dim); if (covarread.is_open()) LOG(SAM) << "Reading covariance matrix from file: " << covarFile << std::endl; else ERR(FTL) << "Covariance matrix provided is incorrect: " << covarFile << std::endl; // 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) { ERR(WRN) << "SetCovarMatrixFromText -> Covariance matrix only has <= 1 " "entries on this line: " << row << std::endl; } 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(), ifstream::in); this->covar = new TMatrixDSym(dim); this->fFullCovar = new TMatrixDSym(dim); if (corr.is_open()) LOG(SAM) << "Reading and converting correlation matrix from file: " << corrFile << std::endl; else { ERR(FTL) << "Correlation matrix provided is incorrect: " << corrFile << std::endl; exit(-1); } 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) { ERR(FTL) << "Found a zero value in the covariance matrix, assuming " "this is an error!" << std::endl; exit(-1); } (*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) { //******************************************************************** LOG(SAM) << "Getting covariance from " << covarFile << "->" << covName << std::endl; TFile* tempFile = new TFile(covarFile.c_str(), "READ"); TH2D* covPlot = (TH2D*)tempFile->Get(covName.c_str()); covPlot->SetDirectory(0); // 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(), 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/MCStudies/Smearceptance_Tester.cxx b/src/MCStudies/Smearceptance_Tester.cxx index db06681..072f2bd 100644 --- a/src/MCStudies/Smearceptance_Tester.cxx +++ b/src/MCStudies/Smearceptance_Tester.cxx @@ -1,790 +1,884 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #include "Smearceptance_Tester.h" #include "SmearceptanceUtils.h" #include "Smearcepterton.h" //#define DEBUG_SMEARTESTER //******************************************************************** /// @brief Class to perform smearceptance MC Studies on a custom measurement Smearceptance_Tester::Smearceptance_Tester(nuiskey samplekey) { //******************************************************************** samplekey.Print(); // Sample overview --------------------------------------------------- std::string descrip = "Simple measurement class for producing an event summary tree of smeared " "events.\n"; if (Config::HasPar("NPOT")) { samplekey.SetS("NPOT", Config::GetParS("NPOT")); } if (Config::HasPar("FluxIntegralOverride")) { samplekey.SetS("FluxIntegralOverride", Config::GetParS("FluxIntegralOverride")); } if (Config::HasPar("TargetVolume")) { samplekey.SetS("TargetVolume", Config::GetParS("TargetVolume")); } if (Config::HasPar("TargetMaterialDensity")) { samplekey.SetS("TargetMaterialDensity", Config::GetParS("TargetMaterialDensity")); } OutputSummaryTree = true; if (Config::HasPar("smear.OutputSummaryTree")) { OutputSummaryTree = Config::GetParI("smear.OutputSummaryTree"); } // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetTitle("Smearceptance Studies"); fSettings.SetDescription(descrip); fSettings.SetXTitle("XXX"); fSettings.SetYTitle("Number of events"); fSettings.SetEnuRange(0.0, 1E5); fSettings.SetAllowedTypes("EVT/SHAPE/DIAG", "EVT/SHAPE/DIAG"); fSettings.DefineAllowedTargets("*"); fSettings.DefineAllowedSpecies("*"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) / TotalIntegratedFlux(); // Measurement Details std::vector splitName = GeneralUtils::ParseToStr(fName, "_"); size_t firstUS = fName.find_first_of("_"); std::string smearceptorName = samplekey.GetS("smearceptor"); QLOG(SAM, "Using smearceptor: " << smearceptorName << " (parsed from: " << fName << ")."); fDataHist = new TH1D(("empty_data"), ("empty-data"), 1, 0, 1); SetupDefaultHist(); fFullCovar = StatUtils::MakeDiagonalCovarMatrix(fDataHist); covar = StatUtils::GetInvert(fFullCovar); eventVariables = NULL; QLOG(SAM, "Smearceptance Flux Scaling Factor = " << fScaleFactor); if (fScaleFactor <= 0.0) { ERROR(WRN, "SCALE FACTOR TOO LOW "); sleep(20); } // Setup our TTrees AddEventVariablesToTree(); smearceptor = &Smearcepterton::Get().GetSmearcepter(smearceptorName); Int_t RecNBins = 20, TrueNBins = 20; double RecBinL = 0xdeadbeef, TrueBinL = 0, RecBinH = 10, TrueBinH = 10; if (Config::HasPar("smear.reconstructed.binning")) { std::vector args = GeneralUtils::ParseToStr( Config::GetParS("smear.reconstructed.binning"), ","); RecNBins = GeneralUtils::StrToInt(args[0]); RecBinL = GeneralUtils::StrToDbl(args[1]); RecBinH = GeneralUtils::StrToDbl(args[2]); TrueNBins = RecNBins; TrueBinL = RecBinL; TrueBinH = RecBinH; } if (Config::HasPar("smear.true.binning")) { std::vector args = GeneralUtils::ParseToStr(Config::GetParS("smear.true.binning"), ","); TrueNBins = GeneralUtils::StrToInt(args[0]); TrueBinL = GeneralUtils::StrToDbl(args[1]); TrueBinH = GeneralUtils::StrToDbl(args[2]); } SVDTruncation = 0; if (Config::HasPar("smear.true.binning")) { SVDTruncation = Config::GetParI("smear.SVD.truncation"); QLOG(SAM, "Applying SVD truncation of: " << SVDTruncation) } ETrueDistrib = NULL; ETrueDistrib_noweight = NULL; ERecDistrib = NULL; RecoSmear = NULL; if (RecBinL != 0xdeadbeef) { QLOG(SAM, "Using binning True: " << TrueNBins << ", [" << TrueBinL << " -- " << TrueBinH << "], Rec: " << RecNBins << ", [" << RecBinL << " -- " << RecBinH << "]"); ETrueDistrib = new TH1D("ELep_rate", ";True E_{#nu};Count", TrueNBins, TrueBinL, TrueBinH); ETrueDistrib_noweight = new TH1D("ELep_rate_noweight", ";True E_{#nu};Count", TrueNBins, TrueBinL, TrueBinH); ERecDistrib = new TH1D("ELepRec_rate", ";Rec E_{#nu};Count", RecNBins, RecBinL, RecBinH); ETrueDistrib->Sumw2(); ERecDistrib->Sumw2(); RecoSmear = new TH2D("ELepHadVis_Recon", ";True E_{#nu};Recon. E_{#nu}", RecNBins, RecBinL, RecBinH, TrueNBins, TrueBinL, TrueBinH); RecoSmear->Sumw2(); } // Final setup --------------------------------------------------- FinaliseMeasurement(); } void Smearceptance_Tester::AddEventVariablesToTree() { if (OutputSummaryTree) { // Setup the TTree to save everything if (!eventVariables) { Config::Get().out->cd(); eventVariables = new TTree((fName + "_VARS").c_str(), (fName + "_VARS").c_str()); } LOG(SAM) << "Adding Event Variables" << std::endl; eventVariables->Branch("Omega_true", &Omega_true, "Omega_true/F"); eventVariables->Branch("Q2_true", &Q2_true, "Q2_true/F"); eventVariables->Branch("Mode_true", &Mode_true, "Mode_true/I"); eventVariables->Branch("EISLep_true", &EISLep_true, "EISLep_true/F"); - eventVariables->Branch("HMFS_mu_true", &HMFS_mu_true); + eventVariables->Branch("HMFS_clep_true", &HMFS_clep_true); eventVariables->Branch("HMFS_pip_true", &HMFS_pip_true); eventVariables->Branch("HMFS_pim_true", &HMFS_pim_true); eventVariables->Branch("HMFS_cpi_true", &HMFS_cpi_true); + eventVariables->Branch("HMFS_pi0_true", &HMFS_pi0_true); + eventVariables->Branch("HMFS_cK_true", &HMFS_cK_true); + eventVariables->Branch("HMFS_K0_true", &HMFS_K0_true); eventVariables->Branch("HMFS_p_true", &HMFS_p_true); eventVariables->Branch("KEFSHad_cpip_true", &KEFSHad_cpip_true, "KEFSHad_cpip_true/F"); eventVariables->Branch("KEFSHad_cpim_true", &KEFSHad_cpim_true, "KEFSHad_cpim_true/F"); eventVariables->Branch("KEFSHad_cpi_true", &KEFSHad_cpi_true, "KEFSHad_cpi_true/F"); eventVariables->Branch("TEFSHad_pi0_true", &TEFSHad_pi0_true, "TEFSHad_pi0_true/F"); + eventVariables->Branch("KEFSHad_cK_true", &KEFSHad_cK_true, + "KEFSHad_cK_true/F"); + eventVariables->Branch("KEFSHad_K0_true", &KEFSHad_K0_true, + "KEFSHad_K0_true/F"); eventVariables->Branch("KEFSHad_p_true", &KEFSHad_p_true, "KEFSHad_p_true/F"); eventVariables->Branch("KEFSHad_n_true", &KEFSHad_n_true, "KEFSHad_n_true/F"); eventVariables->Branch("EFSHad_true", &EFSHad_true, "EFSHad_true/F"); eventVariables->Branch("EFSChargedEMHad_true", &EFSChargedEMHad_true, "EFSChargedEMHad_true/F"); eventVariables->Branch("EFSLep_true", &EFSLep_true, "EFSLep_true/F"); eventVariables->Branch("EFSgamma_true", &EFSgamma_true, "EFSgamma_true/F"); eventVariables->Branch("PDGISLep_true", &PDGISLep_true, "PDGISLep_true/I"); eventVariables->Branch("PDGFSLep_true", &PDGFSLep_true, "PDGFSLep_true/I"); eventVariables->Branch("Nprotons_true", &Nprotons_true, "Nprotons_true/I"); eventVariables->Branch("Nneutrons_true", &Nneutrons_true, "Nneutrons_true/I"); eventVariables->Branch("Ncpiplus_true", &Ncpiplus_true, "Ncpiplus_true/I"); eventVariables->Branch("Ncpiminus_true", &Ncpiminus_true, "Ncpiminus_true/I"); eventVariables->Branch("Ncpi_true", &Ncpi_true, "Ncpi_true/I"); eventVariables->Branch("Npi0_true", &Npi0_true, "Npi0_true/I"); + eventVariables->Branch("NcK_true", &NcK_true, "NcK_true/I"); + eventVariables->Branch("NK0_true", &NK0_true, "NK0_true/I"); - eventVariables->Branch("HMFS_mu_rec", &HMFS_mu_rec); + eventVariables->Branch("HMFS_clep_rec", &HMFS_clep_rec); eventVariables->Branch("HMFS_pip_rec", &HMFS_pip_rec); eventVariables->Branch("HMFS_pim_rec", &HMFS_pim_rec); eventVariables->Branch("HMFS_cpi_rec", &HMFS_cpi_rec); + eventVariables->Branch("HMFS_pi0_rec", &HMFS_pi0_rec); + eventVariables->Branch("HMFS_cK_rec", &HMFS_cK_rec); + eventVariables->Branch("HMFS_K0_rec", &HMFS_K0_rec); eventVariables->Branch("HMFS_p_rec", &HMFS_p_rec); eventVariables->Branch("KEFSHad_cpip_rec", &KEFSHad_cpip_rec, "KEFSHad_cpip_rec/F"); eventVariables->Branch("KEFSHad_cpim_rec", &KEFSHad_cpim_rec, "KEFSHad_cpim_rec/F"); eventVariables->Branch("KEFSHad_cpi_rec", &KEFSHad_cpi_rec, "KEFSHad_cpi_rec/F"); eventVariables->Branch("TEFSHad_pi0_rec", &TEFSHad_pi0_rec, "TEFSHad_pi0_rec/F"); + eventVariables->Branch("KEFSHad_cK_rec", &KEFSHad_cK_rec, + "KEFSHad_cK_rec/F"); + eventVariables->Branch("KEFSHad_K0_rec", &KEFSHad_K0_rec, + "KEFSHad_K0_rec/F"); eventVariables->Branch("KEFSHad_p_rec", &KEFSHad_p_rec, "KEFSHad_p_rec/F"); eventVariables->Branch("KEFSHad_n_rec", &KEFSHad_n_rec, "KEFSHad_n_rec/F"); eventVariables->Branch("EFSHad_rec", &EFSHad_rec, "EFSHad_rec/F"); eventVariables->Branch("EFSLep_rec", &EFSLep_rec, "EFSLep_rec/F"); eventVariables->Branch("EFSVis_cpip", &EFSVis_cpip, "EFSVis_cpip/F"); eventVariables->Branch("EFSVis_cpim", &EFSVis_cpim, "EFSVis_cpim/F"); eventVariables->Branch("EFSVis_cpi", &EFSVis_cpi, "EFSVis_cpi/F"); eventVariables->Branch("EFSVis_pi0", &EFSVis_pi0, "EFSVis_pi0/F"); + eventVariables->Branch("EFSVis_cK", &EFSVis_cK, "EFSVis_cK/F"); + eventVariables->Branch("EFSVis_K0", &EFSVis_K0, "EFSVis_K0/F"); eventVariables->Branch("EFSVis_p", &EFSVis_p, "EFSVis_p/F"); eventVariables->Branch("EFSVis_n", &EFSVis_n, "EFSVis_n/F"); eventVariables->Branch("EFSVis_gamma", &EFSVis_gamma, "EFSVis_gamma/F"); eventVariables->Branch("EFSVis_other", &EFSVis_other, "EFSVis_other/F"); eventVariables->Branch("EFSVis", &EFSVis, "EFSVis/F"); eventVariables->Branch("FSCLep_seen", &FSCLep_seen, "FSCLep_seen/I"); eventVariables->Branch("Nprotons_seen", &Nprotons_seen, "Nprotons_seen/I"); eventVariables->Branch("Nneutrons_seen", &Nneutrons_seen, "Nneutrons_seen/I"); eventVariables->Branch("Ncpip_seen", &Ncpip_seen, "Ncpip_seen/I"); eventVariables->Branch("Ncpim_seen", &Ncpim_seen, "Ncpim_seen/I"); eventVariables->Branch("Ncpi_seen", &Ncpi_seen, "Ncpi_seen/I"); eventVariables->Branch("Npi0_seen", &Npi0_seen, "Npi0_seen/I"); + eventVariables->Branch("NcK_seen", &NcK_seen, "NcK_seen/I"); + eventVariables->Branch("NK0_seen", &NK0_seen, "NK0_seen/I"); eventVariables->Branch("Nothers_seen", &Nothers_seen, "Nothers_seen/I"); eventVariables->Branch("EISLep_QE_rec", &EISLep_QE_rec, "EISLep_QE_rec/F"); eventVariables->Branch("EISLep_LepHad_rec", &EISLep_LepHad_rec, "EISLep_LepHad_rec/F"); eventVariables->Branch("EISLep_LepHadVis_rec", &EISLep_LepHadVis_rec, "EISLep_LepHadVis_rec/F"); eventVariables->Branch("Nprotons_contributed", &Nprotons_contributed, "Nprotons_contributed/I"); eventVariables->Branch("Nneutrons_contributed", &Nneutrons_contributed, "Nneutrons_contributed/I"); eventVariables->Branch("Ncpip_contributed", &Ncpip_contributed, "Ncpip_contributed/I"); eventVariables->Branch("Ncpim_contributed", &Ncpim_contributed, "Ncpim_contributed/I"); eventVariables->Branch("Ncpi_contributed", &Ncpi_contributed, "Ncpi_contributed/I"); eventVariables->Branch("Npi0_contributed", &Npi0_contributed, "Npi0_contributed/I"); + eventVariables->Branch("NcK_contributed", &NcK_contributed, + "NcK_contributed/I"); + eventVariables->Branch("NK0_contributed", &NK0_contributed, + "NK0_contributed/I"); eventVariables->Branch("Ngamma_contributed", &Ngamma_contributed, "Ngamma_contributed/I"); eventVariables->Branch("Nothers_contibuted", &Nothers_contibuted, "Nothers_contibuted/I"); eventVariables->Branch("Weight", &Weight, "Weight/F"); eventVariables->Branch("RWWeight", &RWWeight, "RWWeight/F"); eventVariables->Branch("InputWeight", &InputWeight, "InputWeight/F"); eventVariables->Branch("FluxWeight", &FluxWeight, "FluxWeight/F"); eventVariables->Branch("EffWeight", &EffWeight, "EffWeight/F"); xsecScaling = fScaleFactor; eventVariables->Branch("xsecScaling", &xsecScaling, "xsecScaling/F"); eventVariables->Branch("flagCCINC_true", &flagCCINC_true, "flagCCINC_true/O"); + eventVariables->Branch("flagCC0K_true", &flagCC0K_true, + "flagCC0K_true/O"); eventVariables->Branch("flagCC0Pi_true", &flagCC0Pi_true, "flagCC0Pi_true/O"); eventVariables->Branch("flagCC1Pi_true", &flagCC1Pi_true, "flagCC1Pi_true/O"); eventVariables->Branch("flagCCINC_rec", &flagCCINC_rec, "flagCCINC_rec/O"); + eventVariables->Branch("flagCC0K_rec", &flagCC0K_rec, "flagCC0K_rec/O"); eventVariables->Branch("flagCC0Pi_rec", &flagCC0Pi_rec, "flagCC0Pi_rec/O"); eventVariables->Branch("flagCC1Pi_rec", &flagCC1Pi_rec, "flagCC1Pi_rec/O"); } PredEvtRateWeight = 1; if (fEvtRateScaleFactor != 0xdeadbeef) { if (OutputSummaryTree) { eventVariables->Branch("PredEvtRateWeight", &PredEvtRateWeight, "PredEvtRateWeight/F"); } PredEvtRateWeight = fScaleFactor * fEvtRateScaleFactor; } } template int CountNPdgsSeen(RecoInfo ri, int const (&pdgs)[N]) { int sum = 0; for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) { sum += std::count(ri.RecObjClass.begin(), ri.RecObjClass.end(), pdgs[pdg_it]); } return sum; } template int CountNNotPdgsSeen(RecoInfo ri, int const (&pdgs)[N]) { int sum = 0; for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) { sum += (std::count(ri.RecObjClass.begin(), ri.RecObjClass.end(), pdgs[pdg_it]) ? 0 : 1); } return sum; } template int CountNPdgsContributed(RecoInfo ri, int const (&pdgs)[N]) { int sum = 0; for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) { sum += std::count(ri.TrueContribPDGs.begin(), ri.TrueContribPDGs.end(), pdgs[pdg_it]); } return sum; } template int CountNNotPdgsContributed(RecoInfo ri, int const (&pdgs)[N]) { int sum = 0; for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) { sum += (std::count(ri.TrueContribPDGs.begin(), ri.TrueContribPDGs.end(), pdgs[pdg_it]) ? 0 : 1); } return sum; } TLorentzVector GetHMFSRecParticles(RecoInfo ri, int pdg) { TLorentzVector mom(0, 0, 0, 0); for (size_t p_it = 0; p_it < ri.RecObjMom.size(); ++p_it) { if ((ri.RecObjClass[p_it] == pdg) && (mom.Mag() < ri.RecObjMom[p_it].Mag())) { mom.SetXYZM(ri.RecObjMom[p_it].X(), ri.RecObjMom[p_it].Y(), ri.RecObjMom[p_it].Z(), PhysConst::GetMass(ri.RecObjClass[p_it]) * 1.0E3); } } return mom; } template double SumKE_RecoInfo(RecoInfo ri, int const (&pdgs)[N], double mass) { double sum = 0; for (size_t p_it = 0; p_it < ri.RecObjMom.size(); ++p_it) { if (!std::count(pdgs, pdgs + N, ri.RecObjClass[p_it])) { // If we don't care about this // particle type. continue; } sum += sqrt(ri.RecObjMom[p_it].Mag2() + mass * mass) - mass; } return sum; } template double SumTE_RecoInfo(RecoInfo ri, int const (&pdgs)[N], double mass) { double sum = 0; for (size_t p_it = 0; p_it < ri.RecObjMom.size(); ++p_it) { if (!std::count(pdgs, pdgs + N, ri.RecObjClass[p_it])) { // If we don't care about this // particle type. continue; } sum += sqrt(ri.RecObjMom[p_it].Mag2() + mass * mass); } return sum; } template double SumVisE_RecoInfo(RecoInfo ri, int const (&pdgs)[N]) { double sum = 0; for (size_t p_it = 0; p_it < ri.RecVisibleEnergy.size(); ++p_it) { if (!std::count(pdgs, pdgs + N, ri.TrueContribPDGs[p_it])) { // If we don't care about this // particle type. continue; } sum += ri.RecVisibleEnergy[p_it]; } return sum; } template double SumVisE_RecoInfo_NotPdgs(RecoInfo ri, int const (&pdgs)[N]) { double sum = 0; for (size_t p_it = 0; p_it < ri.RecVisibleEnergy.size(); ++p_it) { if (std::count(pdgs, pdgs + N, ri.TrueContribPDGs[p_it])) { // If we know about this // particle type. continue; } sum += ri.RecVisibleEnergy[p_it]; } return sum; } //******************************************************************** void Smearceptance_Tester::FillEventVariables(FitEvent *event) { //******************************************************************** static int const cpipPDG[] = {211}; static int const cpimPDG[] = {-211}; static int const pi0PDG[] = {111}; + + static int const cKPDG[] = {321, -321}; + static int const K0PDG[] = {311, 310, 130}; + static int const ProtonPDG[] = {2212}; static int const NeutronPDG[] = {2112}; static int const GammaPDG[] = {22}; - static int const CLeptonPDGs[] = {11, 13, 15}; - static int const ExplicitPDGs[] = {211, -211, 111, 2212, 2112, 22, - 11, 13, 15, 12, 14, 16}; + static int const CLeptonPDGs[] = {11, 13, 15, -11, -13, -15}; + static int const ExplicitPDGs[] = {211, -211, 111, 321, -321, 311, 310, 130, + 2212, 2112, 22, 11, 13, 15, 12, 14, + 16, -11, -13, -15, -12, -14, -16}; RecoInfo *ri = smearceptor->Smearcept(event); - HMFS_mu_true = TLorentzVector(0, 0, 0, 0); - HMFS_mu_rec = TLorentzVector(0, 0, 0, 0); - FitParticle *fsMu = event->GetHMFSMuon(); - if (fsMu) { - HMFS_mu_true = fsMu->P4(); - HMFS_mu_rec = GetHMFSRecParticles(*ri, 13); + //** START Pions + + HMFS_clep_true = TLorentzVector(0, 0, 0, 0); + HMFS_clep_rec = TLorentzVector(0, 0, 0, 0); + FitParticle *fsCLep = event->GetHMFSParticle(CLeptonPDGs); + if (fsCLep) { + HMFS_clep_true = fsCLep->P4(); + HMFS_clep_rec = GetHMFSRecParticles(*ri, fsCLep->PDG()); } + //** END Charged leptons + + //** START Pions + HMFS_pip_true = TLorentzVector(0, 0, 0, 0); HMFS_pip_rec = TLorentzVector(0, 0, 0, 0); FitParticle *fsPip = event->GetHMFSPiPlus(); if (fsPip) { HMFS_pip_true = fsPip->P4(); - HMFS_pip_rec = GetHMFSRecParticles(*ri, 211); + HMFS_pip_rec = GetHMFSRecParticles(*ri, fsPip->PDG()); } HMFS_pim_true = TLorentzVector(0, 0, 0, 0); HMFS_pim_rec = TLorentzVector(0, 0, 0, 0); FitParticle *fsPim = event->GetHMFSPiMinus(); if (fsPim) { HMFS_pim_true = fsPim->P4(); - HMFS_pim_rec = GetHMFSRecParticles(*ri, -211); + HMFS_pim_rec = GetHMFSRecParticles(*ri, fsPim->PDG()); } HMFS_cpi_true = TLorentzVector(0, 0, 0, 0); HMFS_cpi_rec = TLorentzVector(0, 0, 0, 0); if (fsPip || fsPim) { if (!fsPip) { HMFS_cpi_true = HMFS_pim_true; HMFS_cpi_rec = HMFS_pim_rec; } else if (!fsPim) { HMFS_cpi_true = HMFS_pip_true; HMFS_cpi_rec = HMFS_pip_rec; } else { HMFS_cpi_true = (fsPip->p2() > fsPim->p2()) ? HMFS_pip_true : HMFS_pim_true; HMFS_cpi_rec = (fsPip->p2() > fsPim->p2()) ? HMFS_pip_rec : HMFS_pim_rec; } } + HMFS_pi0_true = TLorentzVector(0, 0, 0, 0); + HMFS_pi0_rec = TLorentzVector(0, 0, 0, 0); + FitParticle *fsPi0 = event->GetHMFSPiZero(); + if (fsPi0) { + HMFS_pi0_true = fsPi0->P4(); + HMFS_pi0_rec = GetHMFSRecParticles(*ri, fsPi0->PDG()); + } + + //** END Pions + + //** START Kaons + + HMFS_cK_true = TLorentzVector(0, 0, 0, 0); + HMFS_cK_rec = TLorentzVector(0, 0, 0, 0); + FitParticle *fscK = event->GetHMFSParticle(cKPDG); + if (fscK) { + HMFS_cK_true = fscK->P4(); + HMFS_cK_rec = GetHMFSRecParticles(*ri, fscK->PDG()); + } + + HMFS_K0_true = TLorentzVector(0, 0, 0, 0); + HMFS_K0_rec = TLorentzVector(0, 0, 0, 0); + FitParticle *fsK0 = event->GetHMFSParticle(K0PDG); + if (fsK0) { + HMFS_K0_true = fsK0->P4(); + HMFS_K0_rec = GetHMFSRecParticles(*ri, fsK0->PDG()); + } + + //** END Kaons + + //** START Nucleons + HMFS_p_true = TLorentzVector(0, 0, 0, 0); HMFS_p_rec = TLorentzVector(0, 0, 0, 0); FitParticle *fsP = event->GetHMFSProton(); if (fsP) { HMFS_p_true = fsP->P4(); - HMFS_p_rec = GetHMFSRecParticles(*ri, 2212); + HMFS_p_rec = GetHMFSRecParticles(*ri, fsP->PDG()); } TLorentzVector FourMomentumTransfer = (event->GetHMISAnyLeptons()->P4() - event->GetHMFSAnyLeptons()->P4()); Omega_true = FourMomentumTransfer.E(); Q2_true = -1 * FourMomentumTransfer.Mag2(); Mode_true = event->Mode; EISLep_true = event->GetHMISAnyLeptons()->E(); KEFSHad_cpip_true = FitUtils::SumTE_PartVect(event->GetAllFSPiPlus()); KEFSHad_cpim_true = FitUtils::SumTE_PartVect(event->GetAllFSPiMinus()); KEFSHad_cpi_true = KEFSHad_cpip_true + KEFSHad_cpim_true; TEFSHad_pi0_true = FitUtils::SumTE_PartVect(event->GetAllFSPiZero()); + KEFSHad_cK_true = FitUtils::SumTE_PartVect(event->GetAllFSParticle(cKPDG)); + KEFSHad_K0_true = FitUtils::SumTE_PartVect(event->GetAllFSParticle(K0PDG)); KEFSHad_p_true = FitUtils::SumKE_PartVect(event->GetAllFSProton()); KEFSHad_n_true = FitUtils::SumKE_PartVect(event->GetAllFSNeutron()); EFSHad_true = KEFSHad_cpi_true + TEFSHad_pi0_true + KEFSHad_p_true + KEFSHad_n_true; - EFSChargedEMHad_true = KEFSHad_cpi_true + TEFSHad_pi0_true + KEFSHad_p_true; + EFSChargedEMHad_true = KEFSHad_cpi_true + TEFSHad_pi0_true + KEFSHad_p_true + + KEFSHad_cK_true + KEFSHad_K0_true; EFSLep_true = event->GetHMFSAnyLeptons()->E(); EFSgamma_true = FitUtils::SumTE_PartVect(event->GetAllFSPhoton()); PDGISLep_true = event->GetHMISAnyLeptons()->PDG(); PDGFSLep_true = event->GetHMFSAnyLeptons()->PDG(); Nprotons_true = event->GetAllFSProton().size(); Nneutrons_true = event->GetAllFSNeutron().size(); Ncpiplus_true = event->GetAllFSPiPlus().size(); Ncpiminus_true = event->GetAllFSPiMinus().size(); Ncpi_true = Ncpiplus_true + Ncpiminus_true; Npi0_true = event->GetAllFSPiZero().size(); + NcK_true = event->GetAllFSParticle(cKPDG).size(); + NK0_true = event->GetAllFSParticle(K0PDG).size(); KEFSHad_cpip_rec = SumKE_RecoInfo(*ri, cpipPDG, PhysConst::mass_cpi * PhysConst::mass_MeV); KEFSHad_cpim_rec = SumKE_RecoInfo(*ri, cpimPDG, PhysConst::mass_cpi * PhysConst::mass_MeV); KEFSHad_cpi_rec = KEFSHad_cpip_rec + KEFSHad_cpim_rec; TEFSHad_pi0_rec = SumTE_RecoInfo(*ri, pi0PDG, PhysConst::mass_pi0 * PhysConst::mass_MeV); + + KEFSHad_cK_rec = + SumKE_RecoInfo(*ri, cKPDG, PhysConst::mass_cK * PhysConst::mass_MeV); + KEFSHad_K0_rec = + SumKE_RecoInfo(*ri, K0PDG, PhysConst::mass_K0 * PhysConst::mass_MeV); + KEFSHad_p_rec = SumKE_RecoInfo(*ri, ProtonPDG, PhysConst::mass_proton * PhysConst::mass_MeV); KEFSHad_n_rec = SumKE_RecoInfo(*ri, NeutronPDG, PhysConst::mass_neutron * PhysConst::mass_MeV); EFSHad_rec = KEFSHad_cpi_rec + TEFSHad_pi0_rec + KEFSHad_p_rec + KEFSHad_n_rec; TLorentzVector FSLepMom_rec(0, 0, 0, 0); if (event->GetHMFSAnyLeptons()) { FSLepMom_rec = GetHMFSRecParticles(*ri, event->GetHMFSAnyLeptons()->PDG()); EFSLep_rec = FSLepMom_rec.E(); } else { EFSLep_rec = 0; } EFSVis_cpip = SumVisE_RecoInfo(*ri, cpipPDG); EFSVis_cpim = SumVisE_RecoInfo(*ri, cpimPDG); EFSVis_cpi = EFSVis_cpip + EFSVis_cpim; EFSVis_pi0 = SumVisE_RecoInfo(*ri, pi0PDG); + EFSVis_cK = SumVisE_RecoInfo(*ri, cKPDG); + EFSVis_K0 = SumVisE_RecoInfo(*ri, K0PDG); EFSVis_p = SumVisE_RecoInfo(*ri, ProtonPDG); EFSVis_n = SumVisE_RecoInfo(*ri, NeutronPDG); EFSVis_gamma = SumVisE_RecoInfo(*ri, GammaPDG); EFSVis_other = SumVisE_RecoInfo_NotPdgs(*ri, ExplicitPDGs); - EFSVis = EFSVis_cpi + EFSVis_pi0 + EFSVis_p + EFSVis_n + EFSVis_gamma; + EFSVis = EFSVis_cpi + EFSVis_pi0 + EFSVis_p + EFSVis_n + EFSVis_gamma + + EFSVis_cK + EFSVis_K0; FSCLep_seen = CountNPdgsSeen(*ri, CLeptonPDGs); Nprotons_seen = CountNPdgsSeen(*ri, ProtonPDG); Nneutrons_seen = CountNPdgsSeen(*ri, NeutronPDG); Ncpip_seen = CountNPdgsSeen(*ri, cpipPDG); Ncpim_seen = CountNPdgsSeen(*ri, cpimPDG); Ncpi_seen = Ncpip_seen + Ncpim_seen; Npi0_seen = CountNPdgsSeen(*ri, pi0PDG); + NcK_seen = CountNPdgsSeen(*ri, cKPDG); + NK0_seen = CountNPdgsSeen(*ri, K0PDG); Nothers_seen = CountNNotPdgsSeen(*ri, ExplicitPDGs); if (FSCLep_seen && (FSLepMom_rec.Mag() > 1E-8)) { EISLep_QE_rec = FitUtils::EnuQErec(FSLepMom_rec.Mag() / 1000.0, FSLepMom_rec.CosTheta(), 34, PDGFSLep_true > 0) * 1000.0; } else { EISLep_QE_rec = 0; } EISLep_LepHad_rec = EFSLep_rec + EFSHad_rec; EISLep_LepHadVis_rec = EFSLep_rec + EFSHad_rec + EFSVis; Nprotons_contributed = CountNPdgsContributed(*ri, ProtonPDG); Nneutrons_contributed = CountNPdgsContributed(*ri, NeutronPDG); Ncpip_contributed = CountNPdgsContributed(*ri, cpipPDG); Ncpim_contributed = CountNPdgsContributed(*ri, cpimPDG); Ncpi_contributed = Ncpip_contributed + Ncpim_contributed; Npi0_contributed = CountNPdgsContributed(*ri, pi0PDG); + NcK_contributed = CountNPdgsContributed(*ri, cKPDG); + NK0_contributed = CountNPdgsContributed(*ri, K0PDG); Ngamma_contributed = CountNPdgsContributed(*ri, GammaPDG); Nothers_contibuted = CountNNotPdgsContributed(*ri, ExplicitPDGs); Weight = event->RWWeight * event->InputWeight; RWWeight = event->RWWeight; InputWeight = event->InputWeight; FluxWeight = GetFluxHistogram()->GetBinContent( GetFluxHistogram()->FindBin(EISLep_true)) / GetFluxHistogram()->Integral(); EffWeight = ri->Weight; flagCCINC_true = PDGFSLep_true & 1; - flagCC0Pi_true = (Ncpi_true + Npi0_true) == 0; - flagCC1Pi_true = (Ncpi_true + Npi0_true) == 1; - - flagCCINC_rec = FSCLep_seen && PDGFSLep_true & 1; - flagCC0Pi_rec = ((Ncpi_seen + Npi0_seen) == 0) && flagCCINC_rec; - flagCC1Pi_rec = ((Ncpi_seen + Npi0_seen) == 1) && flagCCINC_rec; + flagCC0K_true = (NcK_true + NK0_true) == 0; + flagCC0Pi_true = + flagCCINC_true && flagCC0K_true && ((Ncpi_true + Npi0_true) == 0); + flagCC1Pi_true = + flagCCINC_true && flagCC0K_true && ((Ncpi_true + Npi0_true) == 1); + + flagCCINC_rec = FSCLep_seen && flagCCINC_true; + flagCC0K_rec = (NcK_seen + NK0_seen) == 0; + flagCC0Pi_rec = + flagCCINC_rec && flagCC0K_rec && ((Ncpi_seen + Npi0_seen) == 0); + flagCC1Pi_rec = + flagCCINC_rec && flagCC0K_rec && ((Ncpi_seen + Npi0_seen) == 1); if (OutputSummaryTree) { // Fill the eventVariables Tree eventVariables->Fill(); } if (RecoSmear) { RecoSmear->Fill(EISLep_true / 1000.0, flagCCINC_rec ? EISLep_LepHadVis_rec / 1000.0 : -1, Weight); ETrueDistrib_noweight->Fill(EISLep_true / 1000.0, flagCCINC_true ? Weight : 0); ETrueDistrib->Fill(EISLep_true / 1000.0, flagCCINC_true ? Weight * PredEvtRateWeight : 0); ERecDistrib->Fill(EISLep_LepHadVis_rec / 1000.0, flagCCINC_rec ? Weight * PredEvtRateWeight : 0); } }; //******************************************************************** void Smearceptance_Tester::Write(std::string drawOpt) { //******************************************************************** if (OutputSummaryTree) { // First save the TTree eventVariables->Write(); } // Save Flux and Event Histograms too GetInput()->GetFluxHistogram()->Write(); GetInput()->GetEventHistogram()->Write(); if (!RecoSmear) { return; } TH2D *SmearMatrix_ev = static_cast(RecoSmear->Clone("ELepHadVis_Smear_ev")); for (Int_t trueAxis_it = 1; trueAxis_it < RecoSmear->GetXaxis()->GetNbins() + 1; ++trueAxis_it) { double NEISLep = ETrueDistrib_noweight->GetBinContent(trueAxis_it); for (Int_t recoAxis_it = 1; recoAxis_it < RecoSmear->GetYaxis()->GetNbins() + 1; ++recoAxis_it) { if (NEISLep > std::numeric_limits::epsilon()) { SmearMatrix_ev->SetBinContent( trueAxis_it, recoAxis_it, SmearMatrix_ev->GetBinContent(trueAxis_it, recoAxis_it) / NEISLep); } } } ETrueDistrib_noweight->Write(); ETrueDistrib->Write(); ERecDistrib->Write(); RecoSmear->Write(); SmearMatrix_ev->Write(); TH2D *ResponseMatrix_ev = SmearceptanceUtils::SVDGetInverse(SmearMatrix_ev, SVDTruncation); ResponseMatrix_ev->SetName("ResponseMatrix_ev"); ResponseMatrix_ev->Write(); #ifdef DEBUG_SMEARTESTER TMatrixD SmearMatrix_ev_md = SmearceptanceUtils::GetMatrix(SmearMatrix_ev); TH1D *SmearedEvt = static_cast(ERecDistrib->Clone()); SmearedEvt->SetNameTitle("SmearedEvt", ";Rec E_{#nu}; count"); SmearceptanceUtils::PushTH1ThroughMatrixWithErrors( ETrueDistrib, SmearedEvt, SmearMatrix_ev_md, 5000, false); SmearedEvt->Write(); SmearedEvt->Scale(1, "width"); SmearedEvt->SetName("SmearedEvt_bw"); SmearedEvt->Write(); #endif #ifdef __PROB3PP_ENABLED__ FitWeight *fw = FitBase::GetRW(); if (fw->HasRWEngine(kOSCILLATION)) { OscWeightEngine *oscWE = dynamic_cast(fw->GetRWEngine(kOSCILLATION)); TGraph POsc; POsc.Set(1E4 - 1); double min = ETrueDistrib->GetXaxis()->GetBinLowEdge(1); double step = (ETrueDistrib->GetXaxis()->GetBinUpEdge( ETrueDistrib->GetXaxis()->GetNbins()) - ETrueDistrib->GetXaxis()->GetBinLowEdge(1)) / double(1E4); for (size_t i = 1; i < 1E4; ++i) { double enu = min + i * step; double ow = oscWE->CalcWeight(enu, 14); if (ow != ow) { std::cout << "Bad osc weight for ENu: " << enu << std::endl; } POsc.SetPoint(i - 1, enu, ow); } POsc.Write("POsc", TObject::kOverwrite); } #endif TMatrixD ResponseMatrix_evt_md = SmearceptanceUtils::GetMatrix(ResponseMatrix_ev); TH1D *Unfolded_enu_obs = static_cast(ETrueDistrib->Clone()); Unfolded_enu_obs->SetNameTitle("UnfoldedENu_evt", ";True E_{#nu};count"); SmearceptanceUtils::PushTH1ThroughMatrixWithErrors( ERecDistrib, Unfolded_enu_obs, ResponseMatrix_evt_md, 5000, false); Unfolded_enu_obs->Write(); Unfolded_enu_obs->Scale(1, "width"); Unfolded_enu_obs->SetName("UnfoldedENu_evt_bw"); Unfolded_enu_obs->Write(); ETrueDistrib->Scale(1, "width"); ETrueDistrib->SetName("ELep_rate_bw"); ETrueDistrib->Write(); ERecDistrib->Scale(1, "width"); ERecDistrib->SetName("ELepRec_rate_bw"); ERecDistrib->Write(); } // ------------------------------------------------------------------- // Purely MC Plot // Following functions are just overrides to handle this // ------------------------------------------------------------------- //******************************************************************** /// Everything is classed as signal... bool Smearceptance_Tester::isSignal(FitEvent *event) { //******************************************************************** (void)event; return true; }; //******************************************************************** void Smearceptance_Tester::ScaleEvents() { //******************************************************************** // Saving everything to a TTree so no scaling required return; } //******************************************************************** void Smearceptance_Tester::ApplyNormScale(float norm) { //******************************************************************** // Saving everything to a TTree so no scaling required fCurrentNorm = norm; return; } //******************************************************************** void Smearceptance_Tester::FillHistograms() { //******************************************************************** // No Histograms need filling........ return; } //******************************************************************** void Smearceptance_Tester::ResetAll() { //******************************************************************** if (OutputSummaryTree) { eventVariables->Reset(); } return; } //******************************************************************** float Smearceptance_Tester::GetChi2() { //******************************************************************** // No Likelihood to test, purely MC return 0.0; } diff --git a/src/MCStudies/Smearceptance_Tester.h b/src/MCStudies/Smearceptance_Tester.h index dc471db..aa97221 100644 --- a/src/MCStudies/Smearceptance_Tester.h +++ b/src/MCStudies/Smearceptance_Tester.h @@ -1,181 +1,201 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #ifndef Smearceptance_Tester_H_SEEN #define Smearceptance_Tester_H_SEEN #include "Measurement1D.h" #include "ISmearcepter.h" #ifdef __PROB3PP_ENABLED__ #include "OscWeightEngine.h" #endif //******************************************************************** class Smearceptance_Tester : public Measurement1D { //******************************************************************** public: Smearceptance_Tester(nuiskey samplekey); virtual ~Smearceptance_Tester(){}; //! Grab info from event void FillEventVariables(FitEvent *event); //! Fill Custom Histograms void FillHistograms(); //! ResetAll void ResetAll(); //! Scale void ScaleEvents(); //! Norm void ApplyNormScale(float norm); //! Define this samples signal bool isSignal(FitEvent *nvect); //! Write Files void Write(std::string drawOpt); //! Get Chi2 float GetChi2(); void AddEventVariablesToTree(); private: ISmearcepter *smearceptor; TTree *eventVariables; float Omega_true; float Q2_true; int Mode_true; float EISLep_true; - TLorentzVector HMFS_mu_true; + TLorentzVector HMFS_clep_true; TLorentzVector HMFS_pip_true; TLorentzVector HMFS_pim_true; TLorentzVector HMFS_cpi_true; + TLorentzVector HMFS_pi0_true; + TLorentzVector HMFS_cK_true; + TLorentzVector HMFS_K0_true; TLorentzVector HMFS_p_true; float KEFSHad_cpip_true; float KEFSHad_cpim_true; float KEFSHad_cpi_true; float TEFSHad_pi0_true; + float KEFSHad_cK_true; + float KEFSHad_K0_true; float KEFSHad_p_true; float KEFSHad_n_true; float EFSHad_true; float EFSChargedEMHad_true; float EFSLep_true; float EFSgamma_true; int PDGISLep_true; int PDGFSLep_true; int Nprotons_true; int Nneutrons_true; int Ncpiplus_true; int Ncpiminus_true; int Ncpi_true; int Npi0_true; + int NcK_true; + int NK0_true; - TLorentzVector HMFS_mu_rec; + TLorentzVector HMFS_clep_rec; TLorentzVector HMFS_pip_rec; TLorentzVector HMFS_pim_rec; TLorentzVector HMFS_cpi_rec; + TLorentzVector HMFS_pi0_rec; + TLorentzVector HMFS_cK_rec; + TLorentzVector HMFS_K0_rec; TLorentzVector HMFS_p_rec; float KEFSHad_cpip_rec; float KEFSHad_cpim_rec; float KEFSHad_cpi_rec; float TEFSHad_pi0_rec; + float KEFSHad_cK_rec; + float KEFSHad_K0_rec; float KEFSHad_p_rec; float KEFSHad_n_rec; float EFSHad_rec; float EFSLep_rec; float EFSVis_cpip; float EFSVis_cpim; float EFSVis_cpi; float EFSVis_pi0; + float EFSVis_cK; + float EFSVis_K0; float EFSVis_p; float EFSVis_n; float EFSVis_gamma; float EFSVis_other; float EFSVis; int FSCLep_seen; int Nprotons_seen; int Nneutrons_seen; int Ncpip_seen; int Ncpim_seen; int Ncpi_seen; int Npi0_seen; + int NcK_seen; + int NK0_seen; int Nothers_seen; float EISLep_QE_rec; float EISLep_LepHad_rec; float EISLep_LepHadVis_rec; int Nprotons_contributed; int Nneutrons_contributed; int Ncpip_contributed; int Ncpim_contributed; int Ncpi_contributed; int Npi0_contributed; + int NcK_contributed; + int NK0_contributed; int Ngamma_contributed; int Nothers_contibuted; float Weight; float RWWeight; float InputWeight; float FluxWeight; float EffWeight; float PredEvtRateWeight; float xsecScaling; bool flagCCINC_true; + bool flagCC0K_true; bool flagCC0Pi_true; bool flagCC1Pi_true; bool flagCCINC_rec; + bool flagCC0K_rec; bool flagCC0Pi_rec; bool flagCC1Pi_rec; bool OutputSummaryTree; int SVDTruncation; TH2D *RecoSmear; TH1D *ETrueDistrib; TH1D *ETrueDistrib_noweight; TH1D *ERecDistrib; }; #endif diff --git a/src/Utils/PhysConst.h b/src/Utils/PhysConst.h index ecd9772..a88a447 100644 --- a/src/Utils/PhysConst.h +++ b/src/Utils/PhysConst.h @@ -1,105 +1,112 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #ifndef PHYSCONST_H_SEEN #define PHYSCONST_H_SEEN #include #include #include #include #include #include #include #include #include #include #include #include #include "FitLogger.h" /// namespace to contain all physical constants used by NUISANCE namespace PhysConst { -const double mass_proton = 0.93827203; // Proton mass in GeV -const double mass_neutron = 0.93956536; // Neutron mass in GeV -const double mass_nucleon = (mass_proton + mass_neutron)/2.; -const double mass_proton_kg = 1.6727E-27; // Proton mass in kg -const double mass_neutron_kg = 1.6750E-27; // Neutron mass in kg -const double mass_nucleon_kg = (mass_proton_kg + mass_neutron_kg)/2.; +const double mass_proton = 0.93827203; // Proton mass in GeV +const double mass_neutron = 0.93956536; // Neutron mass in GeV +const double mass_nucleon = (mass_proton + mass_neutron) / 2.; +const double mass_proton_kg = 1.6727E-27; // Proton mass in kg +const double mass_neutron_kg = 1.6750E-27; // Neutron mass in kg +const double mass_nucleon_kg = (mass_proton_kg + mass_neutron_kg) / 2.; const double mass_delta = 1.232; // Delta mass in GeV const double mass_muon = 0.10565837; // Muon mass in GeV const double mass_electron = 0.000510998928; // Electron mass in GeV const double mass_cpi = 0.13957; // charged pion mass in GeV const double mass_pi0 = 0.13498; // neutral pion mass in GeV +const double mass_cK = 0.493677; // charged kaon mass in GeV +const double mass_K0 = 0.497611; // neutral kaon mass in GeV inline double GetMass(int pdg) { - switch (pdg) { + switch (abs(pdg)) { case 11: return mass_electron; case 13: return mass_muon; case 111: return mass_pi0; case 211: - case -211: return mass_cpi; + case 321: + return mass_cK; + case 311: + case 310: + case 130: + return mass_K0; case 2112: return mass_neutron; case 2212: return mass_proton; default: { ERROR(WRN, "Attempted to get mass for PDG: " << pdg << ", but it is not catered for. Please add it to " "src/Utils/PhysConst.h"); return -1; } } } const double mass_MeV = 1000; // MeV/GeV const int pdg_neutrinos[] = {12, -12, 14, -14 /*, 16, -16*/}; const int pdg_muons[] = {13, -13}; const int pdg_leptons[] = {11, -11, 13, -13, 15, -15}; const int pdg_all_leptons[] = {11, -11, 13, -13, 15, -15, 12, -12, 14, -14, 16, -16}; const int pdg_pions[] = {211, -211, 111}; const int pdg_charged_pions[] = {211, -211}; const int pdg_strangemesons[] = { 130, 310, 311, 321, 9000311, 9000321, 10311, 10321, 100311, 100321, 9010311, 9010321, 9020311, 9020321, 313, 323, 10313, 10323, 20313, 20323, 100313, 100323, 9000313, 9000323, 30313, 30323, 315, 325, 9000315, 9000325, 10315, 10325, 20315, 20325, 9010315, 9010325, 9020315, 9020325, 317, 327, 9010317, 9010327}; // Just *-1 to cover possibility const int pdg_kplus = 321; const int pdg_antistrangemesons[] = { -130, -310, -311, -321, -9000311, -9000321, -10311, -10321, -100311, -100321, -9010311, -9010321, -9020311, -9020321, -313, -323, -10313, -10323, -20313, -20323, -100313, -100323, -9000313, -9000323, -30313, -30323, -315, -325, -9000315, -9000325, -10315, -10325, -20315, -20325, -9010315, -9010325, -9020315, -9020325, -317, -327, -9010317, -9010327}; } /*! @} */ #endif