Skip to content

Commit

Permalink
Merge pull request #47 from DUNE/feature/dbarrow257/OscProb_DUNE_OldCore
Browse files Browse the repository at this point in the history 
Feature/dbarrow257/osc prob dune old core
  • Loading branch information
@dbarrow257
dbarrow257 authored Feb 12, 2025
2 parents ae83aad + aea6dfb commit 6a0be3f
Show file tree
Hide file tree
Showing 4 changed files with 269 additions and 1 deletion.
128 changes: 128 additions & 0 deletions configs/Variations_Atmospherics.yaml
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,128 @@
General:
OutputFile: "parameter_study_sterile/test_2d_var/test_var_sterile_all_param_thetav1_dm41_1e-3_nue_numu_selec_recoE_recoCos.root"

DUNESamples: ["configs/Samples/AtmSample_nueselec.yaml", "configs/Samples/AtmSample_numuselec.yaml"]

#Nu-FIT
#OscillationParameters: [0.310, 0.582, 0.224, 7.39E-5, 2.5254E-3, -2.498, 25]

#Nu-FIT 5.2 w. SK atm
#OscillationParameters: [0.303, 0.451, 0.02225, 7.41E-5, 2.45E-3, -2.233]
OscillationParameters: [0.303, 0.451, 0.02225, 7.41E-5, 2.45E-3, -2.233, 0., 0., 0., 0.001, 0., 0.]
#OscillationParameters: [0.303, 0.451, 0.02225, 7.41E-5, 2.45E-3, -2.233, -0.5, 0.2, 0.4, 0.3, 0.5, 0.25, 0., 0., 0., 0.2, 0.4, 0.4]

#T2K-like best-fit
#OscillationParameters: [0.307, 0.528, 0.0218, 7.53e-5, 2.509e-3, -1.601, 25]

Systematics:
XsecCovFile: ["configs/CovObjs/xsec_covariance_DUNE_systs_2022a_FD_v3.yaml"]
XsecCovName: "xsec_cov"
XsecStepScale: 0.1
XsecAtGen: false
OscCovFile: ["configs/CovObjs/OscProb_test.yaml"]
OscCovName: "osc_cov"

Fitter:
FitTestLikelihood: false

MCMC:
NSteps: 2000
AutoSave: 10000

Output:
FileName: "TestEventRates.root"
OUTPUTNAME: "TestLLH.root"

ProcessMCMC: No
Seed: 0
Debug: No

"Variations": [
{
"Name": "sin2th_14",
#"OscParDefault": [0.303, 0.451, 0.02225, 7.41E-5, 2.45E-3, -2.233, 0., 0., 0., 0.001, 0., 0.],
"VarValues": [0., 0.001, 0.005, 0.01, 0.05, 0.1, 0.5],
},

{
"Name": "sin2th_24",
#"OscParDefault": [0.303, 0.451, 0.02225, 7.41E-5, 2.45E-3, -2.233, 0., 0., 0., 0.001, 0., 0.],
"VarValues": [0., 0.001, 0.005, 0.01, 0.05, 0.1, 0.5],
},

{
"Name": "sin2th_34",
#"OscParDefault": [0.303, 0.451, 0.02225, 7.41E-5, 2.45E-3, -2.233, 0., 0., 0., 0.001, 0., 0.],
"VarValues": [0., 0.001, 0.005, 0.01, 0.05, 0.1, 0.5],
},

{
"Name": "delm2_14",
"OscParDefault": [0.303, 0.451, 0.02225, 7.41E-5, 2.45E-3, -2.233, 0.03, 0.06, 0.05, 0.001, 0., 0.],
"VarValues": [0., 0.00001, 0.0001, 0.001, 0.01, 0.1, 1.],
},

{
"Name": "delta_14",
"OscParDefault": [0.303, 0.451, 0.02225, 7.41E-5, 2.45E-3, -2.233, 0.03, 0.06, 0.05, 0.001, 0., 0.],
"VarValues": [-2.35619449019, -1.57079632679, -0.78539816339, 0., 0.78539816339, 1.57079632679, 2.35619449019, 3.14159265],
},

{
"Name": "delta_24",
"OscParDefault": [0.303, 0.451, 0.02225, 7.41E-5, 2.45E-3, -2.233, 0.03, 0.06, 0.05, 0.001, 0., 0.],
"VarValues": [-2.35619449019, -1.57079632679, -0.78539816339, 0., 0.78539816339, 1.57079632679, 2.35619449019, 3.14159265],
}

]

"Projections": [
{
"Name": "RecoNuEnergy",
"VarString": "RecoNeutrinoEnergy",
"VarBins": [0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.2,1.4],
"KinematicCuts": [
{
"Range": [0.1,0.8],
"Name": "TrueNuEnergy",
"VarString": "TrueNeutrinoEnergy"
},
{
"Range": [0.4,1.0],
"Name": "TrueCosZ",
"VarString": "TrueCosineZ"
}
],
"CategoryCuts": [
{
"Name": "OscillationChannel_Single",
"VarString": "OscChannel",
"Breakdown": [[0.0], [1.0], [2.0], [3.0], [4.0], [5.0], [6.0], [7.0], [8.0], [9.0], [10.0], [11.0]],
"Names": ["nue_x_nue","nue_x_numu","nue_x_nutau","numu_x_nue","numu_x_numu","numu_x_nutau","nuebar_x_nuebar","nuebar_x_numubar","nuebar_x_nutaubar","numubar_x_nuebar","numubar_x_numubar","numubar_x_nutaubar"],
},
{
"Name": "OscillationChannel_Group",
"VarString": "OscChannel",
"Breakdown": [[0.0, 1.0, 2.0, 3.0, 4.0, 5.0], [6.0, 7.0, 8.0, 9.0, 10.0, 11.0]],
"Names": ["Nu","Nubar"],
},
{
"Name": "Mode_Single",
"VarString": "Mode",
"Breakdown": [[0.0] , [1.0], [2.0], [3.0], [4.0], [5.0], [6.0], [7.0], [8.0], [9.0], [10.0], [11.0], [12.0], [13.0], [14.0], [15.0], [16.0], [17.0], [18.0], [19.0], [20.0], [21.0], [22.0], [23.0], [24.0], [25.0], [26.0]],
},
{
"Name": "Mode_Group",
"VarString": "Mode",
"Breakdown": [ [0.0], [2.0], [3.0], [9.0], [15.0], [16.0], [1.0, 4.0, 5.0, 6.0, 7.0, 8.0, 10.0, 11.0, 12.0, 13.0], [14.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0, 26.0]],
"Names": ["CCQE", "CCDIS", "CCRES", "CCMEC", "NCDIS", "NCRES", "CCOth", "NCOth"],
}
],
},

{
"Name": "TrueNuEnergy",
"VarString": "TrueNeutrinoEnergy",
"VarBins": [0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.0,3.5,4.0,4.5,5.0,6.0,7.0,8.0,9.0,10.0],
}
]
1 change: 1 addition & 0 deletions src/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@ foreach(app
Projections
SigmaVariation
LikelihoodScan
Variations
)

add_executable(${app} ${app}.cpp)
Expand Down
2 changes: 1 addition & 1 deletion src/LikelihoodScan.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ int main(int argc, char * argv[]) {

MACH3LOG_INFO("Event rate for {} : {:<5.2f}", Sample->GetName(), Sample->get1DHist()->Integral());
if (Sample->GetNDim() == 1)
Sample->addData((TH2D*)DUNEHists.back());
Sample->addData((TH1D*)DUNEHists.back());
else if (Sample->GetNDim() == 2)
Sample->addData((TH2D*)DUNEHists.back());
}
Expand Down
139 changes: 139 additions & 0 deletions src/Variations.cpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,139 @@
#include <iostream>
#include <chrono>
#include <iomanip>
#include <vector>

#include <TH1D.h>
#include <TH2D.h>
#include <THStack.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TRint.h>
#include <TLegend.h>
#include <TColor.h>
#include <TMath.h>

#include "mcmc/mcmc.h"
#include "samplePDFDUNE/MaCh3DUNEFactory.h"
#include "samplePDFDUNE/StructsDUNE.h"

//CS YAML "Variations" node expects to be given each parameter to vary with :
// - "Name": name of the parameter as written in the CovObjs YAML file
// - "OscParDefault" array: values that oscillation parameters should be given as a default for this specific parameter variation
// ex : put sterile mixing angles to non zero values to test variations of sterile cp phases
// If "OscParDefault" not given, will use the "OscillationParameters" values specified in "General" node
// - "VarValues" array: values we want the parameter to take

//TODO: Consider merging with SigmaVariations app at some point

int main(int argc, char * argv[]) {
if(argc == 1){
MACH3LOG_ERROR("Usage: bin/EventRatesDUNEBeam config.cfg");
return 1;
}
manager* FitManager = new manager(argv[1]);


//###############################################################################################################################
//Create samplePDFFD objects

covarianceXsec* xsec = nullptr;
covarianceOsc* osc = nullptr;

std::vector<samplePDFFDBase*> DUNEPdfs;
MakeMaCh3DuneInstance(FitManager, DUNEPdfs, xsec, osc);

std::vector<double> oscpars = FitManager->raw()["General"]["OscillationParameters"].as<std::vector<double>>();

//###############################################################################################################################
//Perform reweight and print total integral

MACH3LOG_INFO("=======================================================");
for(samplePDFFDBase* Sample: DUNEPdfs){
Sample->reweight();
MACH3LOG_INFO("Event rate for {} : {:<5.2f}", Sample->GetName(), Sample->get1DHist()->Integral());
}

//###############################################################################################################################
//DB Can't use the core sigma variations as it's entirely set up around the concept of multiple selections per samplePDF object
// Thats not the case in the FD code, which has one selection per samplePDF object
// Consequently have to write out own code

std::vector<covarianceBase*> CovObjs;
CovObjs.emplace_back(xsec);
CovObjs.emplace_back(osc);

MACH3LOG_INFO("=======================================================");

std::string OutputFileName = FitManager->raw()["General"]["OutputFile"].as<std::string>();
TFile* File = TFile::Open(OutputFileName.c_str(),"RECREATE");


for (covarianceBase* CovObj: CovObjs) {
MACH3LOG_INFO("Starting Variations for covarianceBase Object: {}",CovObj->getName());

int nPars = CovObj->getNpars();

for (int iPar=0;iPar<nPars;iPar++) {
std::string ParName = CovObj->GetParName(iPar);

for (auto const &param : FitManager->raw()["Variations"]) {

std::string VarName = (param["Name"].as<std::string>());

if(ParName == VarName) {

MACH3LOG_INFO("\tParameter : {:<30}",ParName);

if(!param["OscParDefault"]){ //if specific default values not specified for the parameter then use global default ones
CovObj->setParameters(oscpars);
}
else {
CovObj->setParameters((param["OscParDefault"].as<std::vector<double>>()));
}

std::vector<double> valVariations = (param["VarValues"].as<std::vector<double>>());

File->cd();
File->mkdir(ParName.c_str());
File->cd(ParName.c_str());

for (size_t iSigVar=0;iSigVar<valVariations.size();iSigVar++) {
double VarVal = valVariations[iSigVar];

MACH3LOG_INFO("\t\tParameter Value : {:<10.7f}",VarVal);
CovObj->setParProp(iPar,VarVal);

for (size_t iSample=0;iSample<DUNEPdfs.size();iSample++) {
std::string SampleName = DUNEPdfs[iSample]->GetName();

File->cd(ParName.c_str());
if (iSigVar == 0) {
File->mkdir((ParName+"/"+SampleName).c_str());
}
File->cd((ParName+"/"+SampleName).c_str());

DUNEPdfs[iSample]->reweight();
TH1* Hist;
if(DUNEPdfs[iSample]->GetNDim() == 1)
Hist = DUNEPdfs[iSample]->get1DHist();
else if(DUNEPdfs[iSample]->GetNDim() == 2)
Hist = DUNEPdfs[iSample]->get2DHist();

MACH3LOG_INFO("\t\t\tSample : {:<30} - Integral : {:<10}",SampleName,Hist->Integral());

Hist->Write(Form("Variation_%.2e",VarVal));
}
}

}
}
}

MACH3LOG_INFO("=======================================================");
}


File->Close();

}

0 comments on commit 6a0be3f

Please sign in to comment.