diff --git a/src/MINERvA/MINERvA_CCDIS_XSec_1DEnu_nu.cxx b/src/MINERvA/MINERvA_CCDIS_XSec_1DEnu_nu.cxx new file mode 100644 index 0000000..cc4a31c --- /dev/null +++ b/src/MINERvA/MINERvA_CCDIS_XSec_1DEnu_nu.cxx @@ -0,0 +1,119 @@ +// 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 "MINERvA_SignalDef.h" +#include "MINERvA_CCDIS_XSec_1DEnu_nu.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 == "") ERR(WRN) << "target " << target << " was not found!" << std::endl; + + 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(){ +//******************************************************************** + + // 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_CCDIS_XSec_1DEnu_nu.h b/src/MINERvA/MINERvA_CCDIS_XSec_1DEnu_nu.h new file mode 100644 index 0000000..bca64dd --- /dev/null +++ b/src/MINERvA/MINERvA_CCDIS_XSec_1DEnu_nu.h @@ -0,0 +1,49 @@ +// 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 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 new file mode 100644 index 0000000..fc31aa8 --- /dev/null +++ b/src/MINERvA/MINERvA_CCDIS_XSec_1DEnu_ratio.cxx @@ -0,0 +1,122 @@ +// 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 "MINERvA_SignalDef.h" +#include "MINERvA_CCDIS_XSec_1DEnu_ratio.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 { + ERR(FTL) << "target " << target << " was not found!" << std::endl; + exit(-1); + } + + 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) ERR(FTL) << "MINERvA CCDIS ratio requires input files in format: NUMERATOR;DENOMINATOR" << std::endl; + 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()); + + // 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); + } + + return; +} diff --git a/src/MINERvA/MINERvA_CCDIS_XSec_1DEnu_ratio.h b/src/MINERvA/MINERvA_CCDIS_XSec_1DEnu_ratio.h new file mode 100644 index 0000000..5fe5ee8 --- /dev/null +++ b/src/MINERvA/MINERvA_CCDIS_XSec_1DEnu_ratio.h @@ -0,0 +1,52 @@ +// 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 MINERVA_CCDIS_XSec_1DEnu_ratio_H_SEEN +#define MINERVA_CCDIS_XSec_1DEnu_ratio_H_SEEN + +// Fit Includes +#include "MeasurementBase.h" +#include "JointMeas1D.h" +#include "MINERvA_CCDIS_XSec_1DEnu_nu.h" + +class MINERvA_CCDIS_XSec_1DEnu_ratio : public JointMeas1D { +public: + + MINERvA_CCDIS_XSec_1DEnu_ratio(nuiskey samplekey); + virtual ~MINERvA_CCDIS_XSec_1DEnu_ratio() {}; + + void MakePlots(); + + // This is a dummy function as it is not required for the ratio (and does bad bad things) + void ScaleEvents(){return;};; + + private: + + // This is a dummy, the signal is defined separately for each sample! + bool isSignal(){return false;}; + + // Need to have the distributions for the numerator and denominator stored separately + MINERvA_CCDIS_XSec_1DEnu_nu * NUM; + MINERvA_CCDIS_XSec_1DEnu_nu * DEN; + Int_t nBins; + std::string target; + +}; + +#endif diff --git a/src/MINERvA/MINERvA_CCDIS_XSec_1Dx_nu.cxx b/src/MINERvA/MINERvA_CCDIS_XSec_1Dx_nu.cxx new file mode 100644 index 0000000..f1cbd25 --- /dev/null +++ b/src/MINERvA/MINERvA_CCDIS_XSec_1Dx_nu.cxx @@ -0,0 +1,120 @@ +// 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 "MINERvA_SignalDef.h" +#include "MINERvA_CCDIS_XSec_1Dx_nu.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 == "") ERR(WRN) << "target " << target << " was not found!" << std::endl; + + 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(); + this->fDataHist = (TH1D*)this->GetMCList().at(0)->Clone(); + + 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_CCDIS_XSec_1Dx_nu.h b/src/MINERvA/MINERvA_CCDIS_XSec_1Dx_nu.h new file mode 100644 index 0000000..c774cbb --- /dev/null +++ b/src/MINERvA/MINERvA_CCDIS_XSec_1Dx_nu.h @@ -0,0 +1,48 @@ +// 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 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 new file mode 100644 index 0000000..74acc1a --- /dev/null +++ b/src/MINERvA/MINERvA_CCDIS_XSec_1Dx_ratio.cxx @@ -0,0 +1,121 @@ +// 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 "MINERvA_SignalDef.h" +#include "MINERvA_CCDIS_XSec_1Dx_ratio.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 { + ERR(FTL) << "target " << target << " was not found!" << std::endl; + exit(-1); + } + + 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) ERR(FTL) << "MINERvA CCDIS ratio requires input files in format: NUMERATOR;DENOMINATOR" << std::endl; + 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()); + + // 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); + } + + return; +} diff --git a/src/MINERvA/MINERvA_CCDIS_XSec_1Dx_ratio.h b/src/MINERvA/MINERvA_CCDIS_XSec_1Dx_ratio.h new file mode 100644 index 0000000..5686c8a --- /dev/null +++ b/src/MINERvA/MINERvA_CCDIS_XSec_1Dx_ratio.h @@ -0,0 +1,52 @@ +// 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 MINERVA_CCDIS_XSec_1Dx_ratio_H_SEEN +#define MINERVA_CCDIS_XSec_1Dx_ratio_H_SEEN + +// Fit Includes +#include "MeasurementBase.h" +#include "JointMeas1D.h" +#include "MINERvA_CCDIS_XSec_1Dx_nu.h" + +class MINERvA_CCDIS_XSec_1Dx_ratio : public JointMeas1D { +public: + + MINERvA_CCDIS_XSec_1Dx_ratio(nuiskey samplekey); + virtual ~MINERvA_CCDIS_XSec_1Dx_ratio() {}; + + void MakePlots(); + + // This is a dummy function as it is not required for the ratio (and does bad bad things) + void ScaleEvents(){return;}; + + private: + + // This is a dummy, the signal is defined separately for each sample! + bool isSignal(){return false;}; + + // Need to have the distributions for the numerator and denominator stored separately + MINERvA_CCDIS_XSec_1Dx_nu * NUM; + MINERvA_CCDIS_XSec_1Dx_nu * DEN; + Int_t nBins; + std::string target; + +}; + +#endif