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