diff --git a/configs/Variations_Atmospherics.yaml b/configs/Variations_Atmospherics.yaml new file mode 100644 index 0000000..d954768 --- /dev/null +++ b/configs/Variations_Atmospherics.yaml @@ -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], + } +] diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 6590947..ad23225 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -6,6 +6,7 @@ foreach(app Projections SigmaVariation LikelihoodScan + Variations ) add_executable(${app} ${app}.cpp) diff --git a/src/LikelihoodScan.cpp b/src/LikelihoodScan.cpp index d38516b..94ef3d5 100644 --- a/src/LikelihoodScan.cpp +++ b/src/LikelihoodScan.cpp @@ -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()); } diff --git a/src/Variations.cpp b/src/Variations.cpp new file mode 100644 index 0000000..7f2c7be --- /dev/null +++ b/src/Variations.cpp @@ -0,0 +1,139 @@ +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#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 DUNEPdfs; + MakeMaCh3DuneInstance(FitManager, DUNEPdfs, xsec, osc); + + std::vector oscpars = FitManager->raw()["General"]["OscillationParameters"].as>(); + + //############################################################################################################################### + //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 CovObjs; + CovObjs.emplace_back(xsec); + CovObjs.emplace_back(osc); + + MACH3LOG_INFO("======================================================="); + + std::string OutputFileName = FitManager->raw()["General"]["OutputFile"].as(); + 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;iParGetParName(iPar); + + for (auto const ¶m : FitManager->raw()["Variations"]) { + + std::string VarName = (param["Name"].as()); + + 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 valVariations = (param["VarValues"].as>()); + + File->cd(); + File->mkdir(ParName.c_str()); + File->cd(ParName.c_str()); + + for (size_t iSigVar=0;iSigVarsetParProp(iPar,VarVal); + + for (size_t iSample=0;iSampleGetName(); + + 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(); + +}