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;
}