diff --git a/src/MINERvA/MINERvA_CCinc_XSec_1DEnu_nu.cxx b/src/MINERvA/MINERvA_CCinc_XSec_1DEnu_nu.cxx index a13dc29..18e7e80 100644 --- a/src/MINERvA/MINERvA_CCinc_XSec_1DEnu_nu.cxx +++ b/src/MINERvA/MINERvA_CCinc_XSec_1DEnu_nu.cxx @@ -1,117 +1,115 @@ // 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_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 == "") ERR(WRN) << "target " << target << " was not found!" << std::endl; // 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; - Enu = Pnu.E()/1000.; + fXVar = Pnu.E()/1000.; ThetaMu = Pnu.Vect().Angle(Pmu.Vect()); - - fXVar = Enu; 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(); + Measurement1D::ScaleEvents(); - this->fMCHist->Scale(this->fScaleFactor, "width"); + // this->fMCHist->Scale(this->fScaleFactor, "width"); - // Proper error scaling - ROOT Freaks out with xsec weights sometimes - for(int i=0; ifMCStat->GetNbinsX();i++) { + // // 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()); - } + // 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 56a993c..47c0829 100644 --- a/src/MINERvA/MINERvA_CCinc_XSec_1DEnu_nu.h +++ b/src/MINERvA/MINERvA_CCinc_XSec_1DEnu_nu.h @@ -1,49 +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_1DEnu_nu_H_SEEN -#define MINERVA_1DEnu_nu_H_SEEN +#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 Enu, ThetaMu; + 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 a274bd7..4acf8e5 100644 --- a/src/MINERvA/MINERvA_CCinc_XSec_1DEnu_ratio.cxx +++ b/src/MINERvA/MINERvA_CCinc_XSec_1DEnu_ratio.cxx @@ -1,167 +1,129 @@ // 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_CCinc_XSec_1DEnu_ratio.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 { ERR(FTL) << "target " << target << " was not found!" << std::endl; exit(-1); } 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) ERR(FTL) << "MINERvA CCinc 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() ); - SetCovarFromTextFile(fSettings.GetCovarInput()); - + // 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); } return; } - -//******************************************************************** -void MINERvA_CCinc_XSec_1DEnu_ratio::SetCovarMatrixFromText(std::string covarFile, int dim) { -//******************************************************************** - - // WARNING this reads in the data CORRELATIONS - // Make a counter to track the line number - int row = 0; - - std::string line; - std::ifstream covar(covarFile.c_str(), ifstream::in); - - this->fFullCovar = new TMatrixDSym(dim); - if (covar.is_open()) LOG(SAM) << "Reading covariance matrix from file: " << covarFile << std::endl; - else ERR(FTL) << "Covariance matrix provided is incorrect: " << covarFile << std::endl; - - while (std::getline(covar >> std::ws, line, '\n')) { - int column = 0; - - // Loop over entries and insert them into matrix - // Multiply by the errors to get the covariance, rather than the correlation matrix - std::vector entries = GeneralUtils::ParseToDbl(line, " "); - for (std::vector::iterator iter = entries.begin(); - iter != entries.end(); iter++) { - - double val = (*iter) * this->fDataHist->GetBinError(row + 1) * this->fDataHist->GetBinError(column + 1); - - (*this->fFullCovar)(row, column) = val; - column++; - } - row++; - } - - // Robust matrix inversion method - this->covar = StatUtils::GetInvert(fFullCovar); - - return; -}; - - diff --git a/src/MINERvA/MINERvA_CCinc_XSec_1DEnu_ratio.h b/src/MINERvA/MINERvA_CCinc_XSec_1DEnu_ratio.h index 357c14a..c7de8b2 100644 --- a/src/MINERvA/MINERvA_CCinc_XSec_1DEnu_ratio.h +++ b/src/MINERvA/MINERvA_CCinc_XSec_1DEnu_ratio.h @@ -1,54 +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_CCinc_XSec_1DEnu_ratio_H_SEEN #define MINERVA_CCinc_XSec_1DEnu_ratio_H_SEEN // Fit Includes #include "MeasurementBase.h" #include "JointMeas1D.h" #include "MINERvA_CCinc_XSec_1DEnu_nu.h" class MINERvA_CCinc_XSec_1DEnu_ratio : public JointMeas1D { public: MINERvA_CCinc_XSec_1DEnu_ratio(nuiskey samplekey); virtual ~MINERvA_CCinc_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;};; - void SetCovarMatrixFromText(std::string covarFile, int dim); - 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_CCinc_XSec_1DEnu_nu * NUM; MINERvA_CCinc_XSec_1DEnu_nu * DEN; Int_t nBins; std::string target; }; #endif diff --git a/src/MINERvA/MINERvA_CCinc_XSec_1Dx_nu.cxx b/src/MINERvA/MINERvA_CCinc_XSec_1Dx_nu.cxx index 6151e1f..fcc7820 100644 --- a/src/MINERvA/MINERvA_CCinc_XSec_1Dx_nu.cxx +++ b/src/MINERvA/MINERvA_CCinc_XSec_1Dx_nu.cxx @@ -1,128 +1,114 @@ // 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_CCinc_XSec_1Dx_nu.h" //******************************************************************** MINERvA_CCinc_XSec_1Dx_nu::MINERvA_CCinc_XSec_1Dx_nu(std::string name, std::string inputfile, std::string type){ //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "MINERvA_CCinc_XSec_1Dx_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("Reconstructed Bjorken x"); fSettings.SetYTitle("d#sigma/dx (cm^{2}/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_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[7] = {0, 0.1, 0.3, 0.7, 0.9, 1.1, 1.5}; CreateDataHistogram(6, binsx); std::string basedir = FitPar::GetDataBase()+"/MINERvA/CCinc/"; std::string smearfilename = "CCinc_"+target+"_x_smear.csv"; SetSmearingMatrix(basedir + smearfilename, 6, 7); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; //******************************************************************** void MINERvA_CCinc_XSec_1Dx_nu::FillEventVariables(FitEvent *event){ //******************************************************************** if (event->NumFSParticle(13) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; ThetaMu = Pnu.Vect().Angle(Pmu.Vect()); TLorentzVector q = Pnu - Pmu; double q0 = q.E()/1000.0; double Emu = (Pmu.E())/1000.0; Enu_rec = Emu + q0; double Q2 = 4*Enu_rec*Emu*sin(ThetaMu/2)*sin(ThetaMu/2); bjork_x = Q2/2./q0/((PhysConst::mass_proton+PhysConst::mass_neutron)/2.); // Average nucleon masses - // if (fIsNumerator) - fXVar = bjork_x; return; } //******************************************************************** bool MINERvA_CCinc_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 (ThetaMu > 0.296706) return false; return true; }; //******************************************************************** void MINERvA_CCinc_XSec_1Dx_nu::ScaleEvents(){ //******************************************************************** - this->fDataHist = (TH1D*)this->GetMCList().at(0)->Clone(); - this->fDataHist->SetNameTitle((this->fName+"_unsmear").c_str(), (this->fName+"_unsmear"+this->fPlotTitles).c_str()); this->ApplySmearingMatrix(); - - 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()); - } - + Measurement1D::ScaleEvents(); } diff --git a/src/MINERvA/MINERvA_CCinc_XSec_1Dx_nu.h b/src/MINERvA/MINERvA_CCinc_XSec_1Dx_nu.h index 7e3d49f..95f88f4 100644 --- a/src/MINERvA/MINERvA_CCinc_XSec_1Dx_nu.h +++ b/src/MINERvA/MINERvA_CCinc_XSec_1Dx_nu.h @@ -1,48 +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_1Dx_nu_H_SEEN -#define MINERVA_1Dx_nu_H_SEEN +#ifndef MINERVA_CCinc_1Dx_nu_H_SEEN +#define MINERVA_CCinc_1Dx_nu_H_SEEN #include "Measurement1D.h" //******************************************************************** class MINERvA_CCinc_XSec_1Dx_nu : public Measurement1D { //******************************************************************** public: MINERvA_CCinc_XSec_1Dx_nu(std::string name, std::string inputfile, std::string type); virtual ~MINERvA_CCinc_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 Enu, Enu_rec, ThetaMu, bjork_x; + double Enu_rec, ThetaMu, bjork_x; std::string target; }; #endif diff --git a/src/MINERvA/MINERvA_CCinc_XSec_1Dx_ratio.cxx b/src/MINERvA/MINERvA_CCinc_XSec_1Dx_ratio.cxx index 6a803c5..3d22215 100644 --- a/src/MINERvA/MINERvA_CCinc_XSec_1Dx_ratio.cxx +++ b/src/MINERvA/MINERvA_CCinc_XSec_1Dx_ratio.cxx @@ -1,198 +1,128 @@ // 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_CCinc_XSec_1Dx_ratio.h" //******************************************************************** MINERvA_CCinc_XSec_1Dx_ratio::MINERvA_CCinc_XSec_1Dx_ratio(nuiskey samplekey) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "MINERvA_CCinc_XSec_1Dx_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("Reconstructed Bjorken x"); fSettings.SetYTitle("d#sigma/dx (cm^{2}/nucleon)"); fSettings.SetAllowedTypes("FIX/DIAG,FULL/MASK", "FIX/FULL"); fSettings.SetEnuRange(0.0, 20.0); // CCQELike plot information fSettings.SetTitle("MINERvA_CCinc_XSec_1Dx_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 { ERR(FTL) << "target " << target << " was not found!" << std::endl; exit(-1); } std::string basedir = FitPar::GetDataBase() + "/MINERvA/CCinc/"; fSettings.SetDataInput( basedir + "CCinc_" + target + "_CH_ratio_x_data.csv" ); fSettings.SetCovarInput( basedir + "CCinc_" + target + "_CH_ratio_x_covar.csv" ); FinaliseSampleSettings(); // Get parsed input files if (fSubInFiles.size() != 2) ERR(FTL) << "MINERvA CCinc 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() ); - SetCovarFromTextFile(fSettings.GetCovarInput()); + // 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_1Dx_nu("MINERvA_CCinc_XSec_1Dx_" + target + "_CH_NUM", inFileNUM, type); DEN = new MINERvA_CCinc_XSec_1Dx_nu("MINERvA_CCinc_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_CCinc_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(); - std::cout << "MakingPlots CCINC X = " << NUM_MC->Integral() << " " << DEN_MC->Integral() << std::endl; - sleep(10); - 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; } - - -//******************************************************************** -void MINERvA_CCinc_XSec_1Dx_ratio::SetCovarMatrixFromText(std::string covarFile, int dim) { -//******************************************************************** - - // WARNING this reads in the data CORRELATIONS - // Make a counter to track the line number - int row = 0; - - std::string line; - std::ifstream covar(covarFile.c_str(), ifstream::in); - - this->covar = new TMatrixDSym(dim); - this->fFullCovar = new TMatrixDSym(dim); - if (covar.is_open()) LOG(SAM) << "Reading covariance matrix from file: " << covarFile << std::endl; - else ERR(FTL) << "Covariance matrix provided is incorrect: " << covarFile << std::endl; - - while (std::getline(covar >> std::ws, line, '\n')) { - int column = 0; - - // Loop over entries and insert them into matrix - // Multiply by the errors to get the covariance, rather than the correlation matrix - std::vector entries = GeneralUtils::ParseToDbl(line, " "); - for (std::vector::iterator iter = entries.begin(); - iter != entries.end(); iter++) { - - double val = (*iter) * this->fDataHist->GetBinError(row + 1) * this->fDataHist->GetBinError(column + 1); - - (*this->covar)(row, column) = val; - (*this->fFullCovar)(row, column) = val; - column++; - } - row++; - } - - // Robust matrix inversion method - TDecompSVD LU = TDecompSVD(*this->covar); - this->covar = new TMatrixDSym(dim, LU .Invert().GetMatrixArray(), ""); - - return; -}; - - - -//******************************************************************** -void MINERvA_CCinc_XSec_1Dx_ratio::Write(std::string drawOpt) { -//******************************************************************** - - LOG(SAM) << "Writing Normal Plots in MINERvA_CCinc_XSec_1Dx_ratio::Write()" << std::endl; - JointMeas1D::Write(drawOpt); - return; - - //this->GetDataList().at(0)->Write(); - // this->GetMCList() .at(0)->Write(); - - if (this->fFullCovar) { - TH2D cov = TH2D((*this->fFullCovar)); - cov.SetNameTitle((this->fName + "_cov").c_str(), (this->fName + "_cov;Bins; Bins;").c_str()); - cov.Write(); - } - - if (this->covar) { - TH2D covinv = TH2D((*this->covar)); - covinv.SetNameTitle((this->fName + "_covinv").c_str(), (this->fName + "_covinv;Bins; Bins;").c_str()); - covinv.Write(); - } - - return; -} diff --git a/src/MINERvA/MINERvA_CCinc_XSec_1Dx_ratio.h b/src/MINERvA/MINERvA_CCinc_XSec_1Dx_ratio.h index fa4e101..d2bae7f 100644 --- a/src/MINERvA/MINERvA_CCinc_XSec_1Dx_ratio.h +++ b/src/MINERvA/MINERvA_CCinc_XSec_1Dx_ratio.h @@ -1,57 +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_CCinc_XSec_1Dx_ratio_H_SEEN #define MINERVA_CCinc_XSec_1Dx_ratio_H_SEEN // Fit Includes #include "MeasurementBase.h" #include "JointMeas1D.h" #include "MINERvA_CCinc_XSec_1Dx_nu.h" class MINERvA_CCinc_XSec_1Dx_ratio : public JointMeas1D { public: MINERvA_CCinc_XSec_1Dx_ratio(nuiskey samplekey); virtual ~MINERvA_CCinc_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;};; - void SetCovarMatrixFromText(std::string covarFile, int dim); - - // Custom write function because the ratio has much more limited information available than normal. - void Write(std::string drawOpt); - 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_CCinc_XSec_1Dx_nu * NUM; MINERvA_CCinc_XSec_1Dx_nu * DEN; Int_t nBins; std::string target; }; #endif