diff --git a/src/MINERvA/MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu.cxx b/src/MINERvA/MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu.cxx index cb828e2..223d162 100644 --- a/src/MINERvA/MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu.cxx +++ b/src/MINERvA/MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu.cxx @@ -1,149 +1,140 @@ // Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #include "MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu.h" #include "MINERvA_SignalDef.h" //******************************************************************** MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu::MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu( nuiskey samplekey) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu sample. \n" "Target: Target;CH (2 INPUTS)\n" "Flux: MINERvA Forward Horn Current Numu \n" "Signal: Any event with 1 muon, 1 proton p>450, no pions"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); fSettings.SetXTitle("Q^{2}_{QE} (GeV^{2})"); fSettings.SetYTitle(" d#sigma/dQ^{2}_{QE} (cm^{2}/GeV^{2}/nucleon)"); fSettings.SetAllowedTypes("FIX/DIAG,FULL/MASK", "FIX/FULL"); fSettings.SetEnuRange(0.0, 100.0); // CCQELike plot information fSettings.SetTitle("MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu"); fIsRatio = true; nBins = 5; target = ""; if (fSettings.Found("name", "TgtRatioFe")) target = "Fe"; else if (fSettings.Found("name", "TgtRatioPb")) target = "Pb"; else if (fSettings.Found("name", "TgtRatioC")) target = "C"; else { NUIS_ABORT("target " << target << " was not found!"); } std::string basedir = FitPar::GetDataBase() + "/MINERvA/CC0pi/"; fSettings.SetDataInput(basedir + "Q2_TgtRatio_" + target + "_data.txt"); fSettings.SetCovarInput(basedir + "Q2_TgtRatio_" + target + "_covar.txt"); FinaliseSampleSettings(); // Get parsed input files if (fSubInFiles.size() != 2) { NUIS_ABORT("MINERvA CC0pi ratio requires input files in format: " "NUMERATOR;DENOMINATOR"); } std::string inFileNUM = fSubInFiles.at(0); std::string inFileDEN = fSubInFiles.at(1); // Scaling Setup --------------------------------------------------- // Ratio of sub classes so non needed // Plot Setup ------------------------------------------------------- SetDataFromTextFile(fSettings.GetDataInput()); SetCorrelationFromTextFile(fSettings.GetCovarInput()); // Setup Experiments ------------------------------------------------------- std::string type = samplekey.GetS("type"); nuiskey samplekey_num = nuiskey("sample"); samplekey_num.Set("name", "MINERvA_CC0pi_XSec_1DQ2_Tgt" + target + "_nu"); samplekey_num.Set("input", inFileNUM); samplekey_num.Set("type", type); nuiskey samplekey_den = nuiskey("sample"); samplekey_den.Set("name", "MINERvA_CC0pi_XSec_1DQ2_TgtCH_nu"); samplekey_den.Set("input", inFileDEN); samplekey_den.Set("type", type); NUM = new MINERvA_CC0pi_XSec_1DQ2_Tgt_nu(samplekey_num); DEN = new MINERvA_CC0pi_XSec_1DQ2_Tgt_nu(samplekey_den); NUM->SetNoData(); DEN->SetNoData(); // Add to chain for processing this->fSubChain.clear(); this->fSubChain.push_back(NUM); this->fSubChain.push_back(DEN); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; //******************************************************************** void MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu::MakePlots() { //******************************************************************** UInt_t sample = 0; for (std::vector::const_iterator expIter = this->fSubChain.begin(); expIter != this->fSubChain.end(); expIter++) { MeasurementBase *exp = static_cast(*expIter); if (sample == 0) this->NUM = static_cast(exp); else if (sample == 1) this->DEN = static_cast(exp); else break; sample++; } // Now make the ratio histogram - TH1D *NUM_MC = (TH1D *)this->NUM->GetMCList().at(0)->Clone(); - TH1D *DEN_MC = (TH1D *)this->DEN->GetMCList().at(0)->Clone(); - - for (int i = 0; i < nBins; ++i) { - double binVal = 0; - double binErr = 0; - - if (DEN_MC->GetBinContent(i + 1) && NUM_MC->GetBinContent(i + 1)) { - binVal = NUM_MC->GetBinContent(i + 1) / DEN_MC->GetBinContent(i + 1); - double fractErrNUM = - NUM_MC->GetBinError(i + 1) / NUM_MC->GetBinContent(i + 1); - double fractErrDEN = - DEN_MC->GetBinError(i + 1) / DEN_MC->GetBinContent(i + 1); - binErr = - binVal * sqrt(fractErrNUM * fractErrNUM + fractErrDEN * fractErrDEN); - } - - this->fMCHist->SetBinContent(i + 1, binVal); - this->fMCHist->SetBinError(i + 1, binErr); - } + TH1D *NUM_MC = (TH1D *)this->NUM->GetMCList().at(0)->Clone(); + TH1D *DEN_MC = (TH1D *)this->DEN->GetMCList().at(0)->Clone(); + this->fMCHist = PlotUtils::GetRatioPlot(NUM_MC, DEN_MC, this->fMCHist); + + // Also make fine histograms for good measure + TH1D *NUM_FINE = (TH1D*)this->NUM->GetFineList().at(0)->Clone(); + TH1D *DEN_FINE = (TH1D*)this->DEN->GetFineList().at(0)->Clone(); + this->fMCFine = PlotUtils::GetRatioPlot(NUM_FINE, DEN_FINE); + this->fMCFine ->SetNameTitle(Form("%s_MC_FINE", fSettings.GetName().c_str()), + Form("%s_MC_FINE%s", fSettings.GetName().c_str(), + fSettings.GetFullTitles().c_str())); return; } diff --git a/src/MINERvA/MINERvA_CCDIS_XSec_1DEnu_nu.cxx b/src/MINERvA/MINERvA_CCDIS_XSec_1DEnu_nu.cxx index e6336f6..ff58ea8 100644 --- a/src/MINERvA/MINERvA_CCDIS_XSec_1DEnu_nu.cxx +++ b/src/MINERvA/MINERvA_CCDIS_XSec_1DEnu_nu.cxx @@ -1,117 +1,111 @@ // Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #include "MINERvA_CCDIS_XSec_1DEnu_nu.h" #include "MINERvA_SignalDef.h" //******************************************************************** MINERvA_CCDIS_XSec_1DEnu_nu::MINERvA_CCDIS_XSec_1DEnu_nu(std::string name, std::string inputfile, std::string type) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "MINERvA_CCDIS_XSec_1DEnu_nu sample."; // Setup common settings fSettings = LoadSampleSettings(name, inputfile, type); fSettings.SetDescription(descrip); fSettings.SetXTitle("E_{#nu} (GeV)"); fSettings.SetYTitle("d#sigma/dE_{#nu} (cm^{2}/GeV/nucleon)"); fSettings.SetAllowedTypes("FIX/DIAG/MASK", "FIX/DIAG"); fSettings.SetEnuRange(5.0, 50.0); fSettings.DefineAllowedTargets("Pb,Fe,C,H"); fSettings.DefineAllowedSpecies("numu"); fSettings.SetTitle("MINERvA_CCDIS_XSec_1DEnu_nu"); target = ""; if (name.find("C12") != std::string::npos) target = "C12"; else if (name.find("Fe56") != std::string::npos) target = "Fe56"; else if (name.find("Pb208") != std::string::npos) target = "Pb208"; if (name.find("DEN") != std::string::npos) target = "CH"; if (target == "") { NUIS_ERR(WRN, "target " << target << " was not found!"); } FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) / this->TotalIntegratedFlux(); // Plot Setup ------------------------------------------------------- double binsx[8] = {5.00, 10.00, 15.00, 20.00, 25.00, 30.00, 40.00, 50.00}; CreateDataHistogram(7, binsx); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; //******************************************************************** void MINERvA_CCDIS_XSec_1DEnu_nu::FillEventVariables(FitEvent *event) { //******************************************************************** if (event->NumFSParticle(13) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; // Need to calculate Q2 and W using the MINERvA method double Enu = Pnu.E() / 1000.; double Emu = Pmu.E() / 1000.; double theta = FitUtils::th(Pnu, Pmu); Q2 = 4 * Enu * Emu * sin(theta / 2.) * sin(theta / 2.); W = sqrt(PhysConst::mass_nucleon * PhysConst::mass_nucleon + 2 * PhysConst::mass_nucleon * (Enu - Emu) - Q2); fXVar = Enu; return; } //******************************************************************** bool MINERvA_CCDIS_XSec_1DEnu_nu::isSignal(FitEvent *event) { //******************************************************************* if (!SignalDef::isCCINC(event, 14, this->EnuMin, this->EnuMax)) return false; // Restrict the phase space to theta < 17 degrees if (!SignalDef::IsRestrictedAngle(event, 14, 13, 17)) return false; if (Q2 < 1.0) return false; if (W < 2.0) return false; return true; }; - -//******************************************************************** -void MINERvA_CCDIS_XSec_1DEnu_nu::ScaleEvents() { - //******************************************************************** - Measurement1D::ScaleEvents(); -} diff --git a/src/MINERvA/MINERvA_CCDIS_XSec_1DEnu_nu.h b/src/MINERvA/MINERvA_CCDIS_XSec_1DEnu_nu.h index b4bb462..daf72b5 100644 --- a/src/MINERvA/MINERvA_CCDIS_XSec_1DEnu_nu.h +++ b/src/MINERvA/MINERvA_CCDIS_XSec_1DEnu_nu.h @@ -1,49 +1,47 @@ // Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #ifndef MINERVA_CCDIS_1DEnu_nu_H_SEEN #define MINERVA_CCDIS_1DEnu_nu_H_SEEN #include "Measurement1D.h" //******************************************************************** class MINERvA_CCDIS_XSec_1DEnu_nu : public Measurement1D { //******************************************************************** public: MINERvA_CCDIS_XSec_1DEnu_nu(std::string name, std::string inputfile, std::string type); virtual ~MINERvA_CCDIS_XSec_1DEnu_nu() {}; // Functions for handling each neut event void FillEventVariables(FitEvent *event); bool isSignal(FitEvent *event); - void ScaleEvents(); - double GetChi2(){return 0.0;}; private: double Q2, W; std::string target; }; #endif diff --git a/src/MINERvA/MINERvA_CCDIS_XSec_1DEnu_ratio.cxx b/src/MINERvA/MINERvA_CCDIS_XSec_1DEnu_ratio.cxx index 135616f..273b04e 100644 --- a/src/MINERvA/MINERvA_CCDIS_XSec_1DEnu_ratio.cxx +++ b/src/MINERvA/MINERvA_CCDIS_XSec_1DEnu_ratio.cxx @@ -1,149 +1,139 @@ // Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #include "MINERvA_CCDIS_XSec_1DEnu_ratio.h" #include "MINERvA_SignalDef.h" //******************************************************************** MINERvA_CCDIS_XSec_1DEnu_ratio::MINERvA_CCDIS_XSec_1DEnu_ratio( nuiskey samplekey) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "MINERvA_CCDIS_XSec_1DEnu_ratio sample."; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); fSettings.SetXTitle("E_{#nu} (GeV)"); fSettings.SetYTitle(" d#sigma/dE_{#nu} (cm^{2}/GeV/nucleon)"); fSettings.SetAllowedTypes("FIX/DIAG,FULL/MASK", "FIX/FULL"); fSettings.SetEnuRange(5.0, 50.0); // CCQELike plot information fSettings.SetTitle("MINERvA_CCDIS_XSec_1DEnu_ratio"); fIsRatio = true; nBins = 7; target = ""; if (fSettings.Found("name", "C12")) target = "C12"; else if (fSettings.Found("name", "Fe56")) target = "Fe56"; else if (fSettings.Found("name", "Pb208")) target = "Pb208"; else { NUIS_ABORT("target " << target << " was not found!"); } fSettings.SetDataInput(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCDIS/CCDIS_" + target + "_CH_ratio_Enu_data.csv"); fSettings.SetCovarInput( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCDIS/CCDIS_" + target + "_CH_ratio_Enu_stat.csv;" + GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCDIS/CCDIS_" + target + "_CH_ratio_Enu_syst.csv"); FinaliseSampleSettings(); // Get parsed input files if (fSubInFiles.size() != 2) { NUIS_ABORT("MINERvA CCDIS ratio requires input files in format: " "NUMERATOR;DENOMINATOR"); } std::string inFileNUM = fSubInFiles.at(0); std::string inFileDEN = fSubInFiles.at(1); // Scaling Setup --------------------------------------------------- // Ratio of sub classes so non needed // Plot Setup ------------------------------------------------------- SetDataFromTextFile(fSettings.GetDataInput()); SetCovarFromMultipleTextFiles(fSettings.GetCovarInput()); // Need to overlay the sqrt covariance diagonals onto the data histogram // False stops it complaining that the data hist errors aren't set yet StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar, 1, false); // Need to scale the covariance by 1E-76... this cancels with the factor of // 1E76 introduced in StatUtils::GetChi2FromCov Who says two wrongs don't make // a right ScaleCovar(1E76); // Setup Experiments ------------------------------------------------------- std::string type = samplekey.GetS("type"); NUM = new MINERvA_CCDIS_XSec_1DEnu_nu( "MINERvA_CCDIS_XSec_1DEnu_" + target + "_CH_NUM", inFileNUM, type); DEN = new MINERvA_CCDIS_XSec_1DEnu_nu( "MINERvA_CCDIS_XSec_1DEnu_" + target + "_CH_DEN", inFileDEN, type); NUM->SetNoData(); DEN->SetNoData(); // Add to chain for processing this->fSubChain.clear(); this->fSubChain.push_back(NUM); this->fSubChain.push_back(DEN); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; //******************************************************************** void MINERvA_CCDIS_XSec_1DEnu_ratio::MakePlots() { //******************************************************************** UInt_t sample = 0; for (std::vector::const_iterator expIter = this->fSubChain.begin(); expIter != this->fSubChain.end(); expIter++) { MeasurementBase *exp = static_cast(*expIter); if (sample == 0) this->NUM = static_cast(exp); else if (sample == 1) this->DEN = static_cast(exp); else break; sample++; } // Now make the ratio histogram - TH1D *NUM_MC = (TH1D *)this->NUM->GetMCList().at(0)->Clone(); - TH1D *DEN_MC = (TH1D *)this->DEN->GetMCList().at(0)->Clone(); - - for (int i = 0; i < nBins; ++i) { - double binVal = 0; - double binErr = 0; - - if (DEN_MC->GetBinContent(i + 1) && NUM_MC->GetBinContent(i + 1)) { - binVal = NUM_MC->GetBinContent(i + 1) / DEN_MC->GetBinContent(i + 1); - double fractErrNUM = - NUM_MC->GetBinError(i + 1) / NUM_MC->GetBinContent(i + 1); - double fractErrDEN = - DEN_MC->GetBinError(i + 1) / DEN_MC->GetBinContent(i + 1); - binErr = - binVal * sqrt(fractErrNUM * fractErrNUM + fractErrDEN * fractErrDEN); - } - - this->fMCHist->SetBinContent(i + 1, binVal); - this->fMCHist->SetBinError(i + 1, binErr); - } + TH1D *NUM_MC = (TH1D *)this->NUM->GetMCList().at(0)->Clone(); + TH1D *DEN_MC = (TH1D *)this->DEN->GetMCList().at(0)->Clone(); + this->fMCHist = PlotUtils::GetRatioPlot(NUM_MC, DEN_MC, this->fMCHist); + + TH1D *NUM_FINE = (TH1D *)this->NUM->GetFineList().at(0)->Clone(); + TH1D *DEN_FINE = (TH1D *)this->DEN->GetFineList().at(0)->Clone(); + this->fMCFine = PlotUtils::GetRatioPlot(NUM_FINE, DEN_FINE); + this->fMCFine ->SetNameTitle(Form("%s_MC_FINE", fSettings.GetName().c_str()), + Form("%s_MC_FINE%s", fSettings.GetName().c_str(), + fSettings.GetFullTitles().c_str())); return; } diff --git a/src/MINERvA/MINERvA_CCDIS_XSec_1Dx_nu.cxx b/src/MINERvA/MINERvA_CCDIS_XSec_1Dx_nu.cxx index 1128bd1..07040b7 100644 --- a/src/MINERvA/MINERvA_CCDIS_XSec_1Dx_nu.cxx +++ b/src/MINERvA/MINERvA_CCDIS_XSec_1Dx_nu.cxx @@ -1,118 +1,111 @@ // Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #include "MINERvA_CCDIS_XSec_1Dx_nu.h" #include "MINERvA_SignalDef.h" //******************************************************************** MINERvA_CCDIS_XSec_1Dx_nu::MINERvA_CCDIS_XSec_1Dx_nu(std::string name, std::string inputfile, std::string type) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "MINERvA_CCDIS_XSec_1Dx_nu sample."; // Setup common settings fSettings = LoadSampleSettings(name, inputfile, type); fSettings.SetDescription(descrip); fSettings.SetXTitle("Reconstructed Bjorken x"); fSettings.SetYTitle("d#sigma/dx (cm^{2}/nucleon)"); fSettings.SetAllowedTypes("FIX/DIAG/MASK", "FIX/DIAG"); fSettings.SetEnuRange(5.0, 50.0); fSettings.DefineAllowedTargets("Pb,Fe,C,H"); fSettings.DefineAllowedSpecies("numu"); fSettings.SetTitle("MINERvA_CCDIS_XSec_1Dx_nu"); target = ""; if (name.find("C12") != std::string::npos) target = "C12"; else if (name.find("Fe56") != std::string::npos) target = "Fe56"; else if (name.find("Pb208") != std::string::npos) target = "Pb208"; if (name.find("DEN") != std::string::npos) target = "CH"; if (target == "") { NUIS_ABORT("target " << target << " was not found!"); } FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) / this->TotalIntegratedFlux(); // Plot Setup ------------------------------------------------------- double binsx[6] = {0.00, 0.10, 0.20, 0.30, 0.40, 0.75}; CreateDataHistogram(5, binsx); FinaliseMeasurement(); }; //******************************************************************** void MINERvA_CCDIS_XSec_1Dx_nu::FillEventVariables(FitEvent *event) { //******************************************************************** if (event->NumFSParticle(13) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; // Need to calculate Q2 and W using the MINERvA method double Enu = Pnu.E() / 1000.; double Emu = Pmu.E() / 1000.; double theta = FitUtils::th(Pnu, Pmu); Q2 = 4 * Enu * Emu * sin(theta / 2.) * sin(theta / 2.); W = sqrt(PhysConst::mass_nucleon * PhysConst::mass_nucleon + 2 * PhysConst::mass_nucleon * (Enu - Emu) - Q2); fXVar = Q2 / (2 * PhysConst::mass_nucleon * (Enu - Emu)); return; } //******************************************************************** bool MINERvA_CCDIS_XSec_1Dx_nu::isSignal(FitEvent *event) { //******************************************************************* // Only look at numu events if (!SignalDef::isCCINC(event, 14, EnuMin, EnuMax)) return false; // Restrict the phase space to theta < 17 degrees if (!SignalDef::IsRestrictedAngle(event, 14, 13, 17)) return false; if (Q2 < 1.0) return false; if (W < 2.0) return false; return true; }; - -//******************************************************************** -void MINERvA_CCDIS_XSec_1Dx_nu::ScaleEvents() { - //******************************************************************** - - Measurement1D::ScaleEvents(); -} diff --git a/src/MINERvA/MINERvA_CCDIS_XSec_1Dx_nu.h b/src/MINERvA/MINERvA_CCDIS_XSec_1Dx_nu.h index a5faa56..a984703 100644 --- a/src/MINERvA/MINERvA_CCDIS_XSec_1Dx_nu.h +++ b/src/MINERvA/MINERvA_CCDIS_XSec_1Dx_nu.h @@ -1,48 +1,46 @@ // Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #ifndef MINERVA_CCDIS_1Dx_nu_H_SEEN #define MINERVA_CCDIS_1Dx_nu_H_SEEN #include "Measurement1D.h" //******************************************************************** class MINERvA_CCDIS_XSec_1Dx_nu : public Measurement1D { //******************************************************************** public: MINERvA_CCDIS_XSec_1Dx_nu(std::string name, std::string inputfile, std::string type); virtual ~MINERvA_CCDIS_XSec_1Dx_nu() {}; // Functions for handling each neut event void FillEventVariables(FitEvent *event); bool isSignal(FitEvent *event); - void ScaleEvents(); - double GetChi2(){return 0.0;}; private: double Q2, W; std::string target; }; #endif diff --git a/src/MINERvA/MINERvA_CCDIS_XSec_1Dx_ratio.cxx b/src/MINERvA/MINERvA_CCDIS_XSec_1Dx_ratio.cxx index 7a3ff90..cd9b137 100644 --- a/src/MINERvA/MINERvA_CCDIS_XSec_1Dx_ratio.cxx +++ b/src/MINERvA/MINERvA_CCDIS_XSec_1Dx_ratio.cxx @@ -1,147 +1,137 @@ // Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #include "MINERvA_CCDIS_XSec_1Dx_ratio.h" #include "MINERvA_SignalDef.h" //******************************************************************** MINERvA_CCDIS_XSec_1Dx_ratio::MINERvA_CCDIS_XSec_1Dx_ratio(nuiskey samplekey) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "MINERvA_CCDIS_XSec_1Dx_ratio sample."; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); fSettings.SetXTitle("Reconstructed Bjorken x"); fSettings.SetYTitle("d#sigma/dx (cm^{2}/nucleon)"); fSettings.SetAllowedTypes("FIX/DIAG,FULL/MASK", "FIX/FULL"); fSettings.SetEnuRange(5.0, 50.0); // CCQELike plot information fSettings.SetTitle("MINERvA_CCDIS_XSec_1Dx_ratio"); fIsRatio = true; nBins = 5; target = ""; if (fSettings.Found("name", "C12")) target = "C12"; else if (fSettings.Found("name", "Fe56")) target = "Fe56"; else if (fSettings.Found("name", "Pb208")) target = "Pb208"; else { NUIS_ABORT("target " << target << " was not found!"); } fSettings.SetDataInput(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCDIS/CCDIS_" + target + "_CH_ratio_x_data.csv"); fSettings.SetCovarInput( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCDIS/CCDIS_" + target + "_CH_ratio_x_stat.csv;" + GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCDIS/CCDIS_" + target + "_CH_ratio_x_syst.csv"); FinaliseSampleSettings(); // Get parsed input files if (fSubInFiles.size() != 2) { NUIS_ABORT("MINERvA CCDIS ratio requires input files in format: " "NUMERATOR;DENOMINATOR"); } std::string inFileNUM = fSubInFiles.at(0); std::string inFileDEN = fSubInFiles.at(1); // Scaling Setup --------------------------------------------------- // Ratio of sub classes so non needed // Plot Setup ------------------------------------------------------- SetDataFromTextFile(fSettings.GetDataInput()); SetCovarFromMultipleTextFiles(fSettings.GetCovarInput()); // Need to overlay the sqrt covariance diagonals (*1E-38) onto the data hist // False stops it complaining that the data hist errors aren't set yet StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar, 1, false); // Need to scale the covariance by 1E-76... this cancels with the factor of // 1E76 introduced in StatUtils::GetChi2FromCov Who says two wrongs don't make // a right ScaleCovar(1E76); // Setup Experiments ------------------------------------------------------- std::string type = samplekey.GetS("type"); NUM = new MINERvA_CCDIS_XSec_1Dx_nu( "MINERvA_CCDIS_XSec_1Dx_" + target + "_CH_NUM", inFileNUM, type); DEN = new MINERvA_CCDIS_XSec_1Dx_nu( "MINERvA_CCDIS_XSec_1Dx_" + target + "_CH_DEN", inFileDEN, type); NUM->SetNoData(); DEN->SetNoData(); // Add to chain for processing this->fSubChain.clear(); this->fSubChain.push_back(NUM); this->fSubChain.push_back(DEN); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; //******************************************************************** void MINERvA_CCDIS_XSec_1Dx_ratio::MakePlots() { //******************************************************************** UInt_t sample = 0; for (std::vector::const_iterator expIter = this->fSubChain.begin(); expIter != this->fSubChain.end(); expIter++) { MeasurementBase *exp = static_cast(*expIter); if (sample == 0) this->NUM = static_cast(exp); else if (sample == 1) this->DEN = static_cast(exp); else break; sample++; } // Now make the ratio histogram TH1D *NUM_MC = (TH1D *)this->NUM->GetMCList().at(0)->Clone(); TH1D *DEN_MC = (TH1D *)this->DEN->GetMCList().at(0)->Clone(); - - for (int i = 0; i < nBins; ++i) { - double binVal = 0; - double binErr = 0; - - if (DEN_MC->GetBinContent(i + 1) && NUM_MC->GetBinContent(i + 1)) { - binVal = NUM_MC->GetBinContent(i + 1) / DEN_MC->GetBinContent(i + 1); - double fractErrNUM = - NUM_MC->GetBinError(i + 1) / NUM_MC->GetBinContent(i + 1); - double fractErrDEN = - DEN_MC->GetBinError(i + 1) / DEN_MC->GetBinContent(i + 1); - binErr = - binVal * sqrt(fractErrNUM * fractErrNUM + fractErrDEN * fractErrDEN); - } - - this->fMCHist->SetBinContent(i + 1, binVal); - this->fMCHist->SetBinError(i + 1, binErr); - } - + this->fMCHist = PlotUtils::GetRatioPlot(NUM_MC, DEN_MC, this->fMCHist); + + // Also the fine hist + TH1D *NUM_FINE = (TH1D *)this->NUM->GetFineList().at(0)->Clone(); + TH1D *DEN_FINE = (TH1D *)this->DEN->GetFineList().at(0)->Clone(); + this->fMCFine = PlotUtils::GetRatioPlot(NUM_FINE, DEN_FINE); + this->fMCFine ->SetNameTitle(Form("%s_MC_FINE", fSettings.GetName().c_str()), + Form("%s_MC_FINE%s", fSettings.GetName().c_str(), + fSettings.GetFullTitles().c_str())); return; } diff --git a/src/MINERvA/MINERvA_CCinc_XSec_1DEnu_nu.cxx b/src/MINERvA/MINERvA_CCinc_XSec_1DEnu_nu.cxx index 118f09b..be5ff81 100644 --- a/src/MINERvA/MINERvA_CCinc_XSec_1DEnu_nu.cxx +++ b/src/MINERvA/MINERvA_CCinc_XSec_1DEnu_nu.cxx @@ -1,125 +1,106 @@ // Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #include "MINERvA_SignalDef.h" #include "MINERvA_CCinc_XSec_1DEnu_nu.h" //******************************************************************** MINERvA_CCinc_XSec_1DEnu_nu::MINERvA_CCinc_XSec_1DEnu_nu(std::string name, std::string inputfile, std::string type) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "MINERvA_CCinc_XSec_1DEnu_nu sample. \n" "Target: CH \n" "Flux: MiniBooNE Forward Horn Current nue + nuebar \n" "Signal: Any event with 1 muon, any nucleons, and no " "other FS particles \n"; // Setup common settings fSettings = LoadSampleSettings(name, inputfile, type); fSettings.SetDescription(descrip); fSettings.SetXTitle("E_{#nu} (GeV)"); fSettings.SetYTitle("d#sigma/dE_{#nu} (cm^{2}/GeV/nucleon)"); fSettings.SetAllowedTypes("FIX/DIAG/MASK", "FIX/DIAG"); fSettings.SetEnuRange(2.0, 20.0); fSettings.DefineAllowedTargets("C,H"); fSettings.DefineAllowedSpecies("numu"); fSettings.SetTitle("MINERvA_CCinc_XSec_1DEnu_nu"); target = ""; if (name.find("C12") != std::string::npos) target = "C12"; else if (name.find("Fe56") != std::string::npos) target = "Fe56"; else if (name.find("Pb208") != std::string::npos) target = "Pb208"; if (name.find("DEN") != std::string::npos) target = "CH"; if (target == "") { NUIS_ABORT("target " << target << " was not found!"); } // fSettings.SetSmearingInput( FitPar::GetDataBase() + // "/MINERvA/CCinc/CCinc_"+target+"_x_smear.csv" ); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) / this->TotalIntegratedFlux(); // Plot Setup ------------------------------------------------------- double binsx[9] = {2, 3, 4, 5, 6, 8, 10, 15, 20}; CreateDataHistogram(8, binsx); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; //******************************************************************** void MINERvA_CCinc_XSec_1DEnu_nu::FillEventVariables(FitEvent *event) { //******************************************************************** if (event->NumFSParticle(13) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; fXVar = Pnu.E() / 1000.; ThetaMu = Pnu.Vect().Angle(Pmu.Vect()); return; } //******************************************************************** bool MINERvA_CCinc_XSec_1DEnu_nu::isSignal(FitEvent *event) { //******************************************************************* if (!SignalDef::isCCINC(event, 14, this->EnuMin, this->EnuMax)) return false; // Restrict the phase space to theta < 17 degrees if (ThetaMu > 0.296706) return false; return true; }; - -//******************************************************************** -void MINERvA_CCinc_XSec_1DEnu_nu::ScaleEvents() { - //******************************************************************** - - // Get rid of this because it causes odd behaviour - Measurement1D::ScaleEvents(); - - // this->fMCHist->Scale(this->fScaleFactor, "width"); - - // // Proper error scaling - ROOT Freaks out with xsec weights sometimes - // for(int i=0; ifMCStat->GetNbinsX();i++) { - - // if (this->fMCStat->GetBinContent(i+1) != 0) - // this->fMCHist->SetBinError(i+1, this->fMCHist->GetBinContent(i+1) * - // this->fMCStat->GetBinError(i+1) / this->fMCStat->GetBinContent(i+1) ); - // else this->fMCHist->SetBinError(i+1, this->fMCHist->Integral()); - // } -} diff --git a/src/MINERvA/MINERvA_CCinc_XSec_1DEnu_nu.h b/src/MINERvA/MINERvA_CCinc_XSec_1DEnu_nu.h index af59e6a..e7cf09e 100644 --- a/src/MINERvA/MINERvA_CCinc_XSec_1DEnu_nu.h +++ b/src/MINERvA/MINERvA_CCinc_XSec_1DEnu_nu.h @@ -1,49 +1,47 @@ // Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #ifndef MINERVA_CCinc_1DEnu_nu_H_SEEN #define MINERVA_CCinc_1DEnu_nu_H_SEEN #include "Measurement1D.h" //******************************************************************** class MINERvA_CCinc_XSec_1DEnu_nu : public Measurement1D { //******************************************************************** public: MINERvA_CCinc_XSec_1DEnu_nu(std::string name, std::string inputfile, std::string type); virtual ~MINERvA_CCinc_XSec_1DEnu_nu() {}; // Functions for handling each neut event void FillEventVariables(FitEvent *event); bool isSignal(FitEvent *event); - void ScaleEvents(); - double GetChi2(){return 0.0;}; private: double ThetaMu; std::string target; }; #endif diff --git a/src/MINERvA/MINERvA_CCinc_XSec_1DEnu_ratio.cxx b/src/MINERvA/MINERvA_CCinc_XSec_1DEnu_ratio.cxx index a68660d..00671ba 100644 --- a/src/MINERvA/MINERvA_CCinc_XSec_1DEnu_ratio.cxx +++ b/src/MINERvA/MINERvA_CCinc_XSec_1DEnu_ratio.cxx @@ -1,145 +1,136 @@ // Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #include "MINERvA_CCinc_XSec_1DEnu_ratio.h" #include "MINERvA_SignalDef.h" //******************************************************************** MINERvA_CCinc_XSec_1DEnu_ratio::MINERvA_CCinc_XSec_1DEnu_ratio( nuiskey samplekey) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "MINERvA_CCinc_XSec_1DEnu_ratio sample. \n" "Target: CH \n" "Flux: MINERvA Forward Horn Current nue + nuebar \n" "Signal: Any event with 1 electron, any nucleons, and " "no other FS particles \n"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); fSettings.SetXTitle("E_{#nu} (GeV)"); fSettings.SetYTitle(" d#sigma/dE_{#nu} (cm^{2}/GeV/nucleon)"); fSettings.SetAllowedTypes("FIX/DIAG,FULL/MASK", "FIX/FULL"); fSettings.SetEnuRange(0.0, 20.0); // CCQELike plot information fSettings.SetTitle("MINERvA_CCinc_XSec_1DEnu_ratio"); fIsRatio = true; nBins = 8; target = ""; if (fSettings.Found("name", "C12")) target = "C12"; else if (fSettings.Found("name", "Fe56")) target = "Fe56"; else if (fSettings.Found("name", "Pb208")) target = "Pb208"; else { NUIS_ABORT("target " << target << " was not found!"); } std::string basedir = FitPar::GetDataBase() + "/MINERvA/CCinc/"; fSettings.SetDataInput(basedir + "CCinc_" + target + "_CH_ratio_Enu_data.csv"); fSettings.SetCovarInput(basedir + "CCinc_" + target + "_CH_ratio_Enu_covar.csv"); FinaliseSampleSettings(); // Get parsed input files if (fSubInFiles.size() != 2) { NUIS_ABORT("MINERvA CCinc ratio requires input files in format: " "NUMERATOR;DENOMINATOR"); } std::string inFileNUM = fSubInFiles.at(0); std::string inFileDEN = fSubInFiles.at(1); // Scaling Setup --------------------------------------------------- // Ratio of sub classes so non needed // Plot Setup ------------------------------------------------------- SetDataFromTextFile(fSettings.GetDataInput()); // This function forces in a factor of 1E76 to the covariance. // This cancels with the factor of 1E38 which is added to the data in the chi2 // calculation... Who said two wrongs don't make a right? SetCorrelationFromTextFile(fSettings.GetCovarInput()); // Setup Experiments ------------------------------------------------------- std::string type = samplekey.GetS("type"); NUM = new MINERvA_CCinc_XSec_1DEnu_nu( "MINERvA_CCinc_XSec_1DEnu_" + target + "_CH_NUM", inFileNUM, type); DEN = new MINERvA_CCinc_XSec_1DEnu_nu( "MINERvA_CCinc_XSec_1DEnu_" + target + "_CH_DEN", inFileDEN, type); NUM->SetNoData(); DEN->SetNoData(); // Add to chain for processing this->fSubChain.clear(); this->fSubChain.push_back(NUM); this->fSubChain.push_back(DEN); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; //******************************************************************** void MINERvA_CCinc_XSec_1DEnu_ratio::MakePlots() { //******************************************************************** UInt_t sample = 0; for (std::vector::const_iterator expIter = this->fSubChain.begin(); expIter != this->fSubChain.end(); expIter++) { MeasurementBase *exp = static_cast(*expIter); if (sample == 0) this->NUM = static_cast(exp); else if (sample == 1) this->DEN = static_cast(exp); else break; sample++; } // Now make the ratio histogram TH1D *NUM_MC = (TH1D *)this->NUM->GetMCList().at(0)->Clone(); TH1D *DEN_MC = (TH1D *)this->DEN->GetMCList().at(0)->Clone(); - - for (int i = 0; i < nBins; ++i) { - double binVal = 0; - double binErr = 0; - - if (DEN_MC->GetBinContent(i + 1) && NUM_MC->GetBinContent(i + 1)) { - binVal = NUM_MC->GetBinContent(i + 1) / DEN_MC->GetBinContent(i + 1); - double fractErrNUM = - NUM_MC->GetBinError(i + 1) / NUM_MC->GetBinContent(i + 1); - double fractErrDEN = - DEN_MC->GetBinError(i + 1) / DEN_MC->GetBinContent(i + 1); - binErr = - binVal * sqrt(fractErrNUM * fractErrNUM + fractErrDEN * fractErrDEN); - } - - this->fMCHist->SetBinContent(i + 1, binVal); - this->fMCHist->SetBinError(i + 1, binErr); - } + this->fMCHist = PlotUtils::GetRatioPlot(NUM_MC, DEN_MC, this->fMCHist); + + // Fine ratio for good measure + TH1D *NUM_FINE = (TH1D *)this->NUM->GetFineList().at(0)->Clone(); + TH1D *DEN_FINE = (TH1D *)this->DEN->GetFineList().at(0)->Clone(); + this->fMCFine = PlotUtils::GetRatioPlot(NUM_FINE, DEN_FINE); + this->fMCFine ->SetNameTitle(Form("%s_MC_FINE", fSettings.GetName().c_str()), + Form("%s_MC_FINE%s", fSettings.GetName().c_str(), + fSettings.GetFullTitles().c_str())); return; }