diff --git a/src/T2K/T2K_NuMu_CC0pi_OC_XSec_2DPcos.cxx b/src/T2K/T2K_NuMu_CC0pi_OC_XSec_2DPcos.cxx index 75e8746..f3c1c68 100644 --- a/src/T2K/T2K_NuMu_CC0pi_OC_XSec_2DPcos.cxx +++ b/src/T2K/T2K_NuMu_CC0pi_OC_XSec_2DPcos.cxx @@ -1,274 +1,215 @@ // Copyright 2016-2021 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 "T2K_NuMu_CC0pi_OC_XSec_2DPcos.h" static size_t nangbins = 6; static double angular_binning_costheta[] = {-1, 0., 0.6, 0.75, 0.86, 0.93, 1.}; //******************************************************************** T2K_NuMu_CC0pi_OC_XSec_2DPcos::T2K_NuMu_CC0pi_OC_XSec_2DPcos(nuiskey samplekey) { //******************************************************************** // Samples overview --------------------------------------------------- fSettings = LoadSampleSettings(samplekey); - std::string name = fSettings.GetS("name"); + std::string name = fSettings.GetName(); std::string descrip = ""; // This has to deal with NuMu on O and on C if (!name.compare("T2K_NuMu_CC0pi_O_XSec_2DPcos")){ descrip = name +". \n" "Target: Oxygen \n" "Flux: T2K 2.5 degree off-axis (ND280) \n" "Signal: CC0pi\n" "DOI:10.1103/PhysRevD.101.112004"; fSettings.SetTitle("T2K_NuMu_CC0pi_O_XSec_2DPcos"); fSettings.DefineAllowedTargets("O"); Target = "O"; } else if (!name.compare("T2K_NuMu_CC0pi_C_XSec_2DPcos")){ descrip = name +". \n" "Target: Carbon \n" "Flux: T2K 2.5 degree off-axis (ND280) \n" "Signal: CC0pi\n" "DOI:10.1103/PhysRevD.101.112004"; fSettings.SetTitle("T2K_NuMu_CC0pi_C_XSec_2DPcos"); fSettings.DefineAllowedTargets("C"); Target = "C"; } // Setup common settings fSettings.SetDescription(descrip); fSettings.SetXTitle("p_{#mu}-cos#theta_{#mu}"); - fSettings.SetYTitle("d^{2}#sigma/dP_{#mu}dcos#theta_{#mu} (cm^{2}/GeV)"); + fSettings.SetYTitle("d^{2}#sigma/dp_{#mu}dcos#theta_{#mu} (cm^{2}/GeV)"); fSettings.SetAllowedTypes("DIAG,FULL/FREE,SHAPE,FIX/SYSTCOV/STATCOV","FIX"); fSettings.SetEnuRangeFromFlux(fFluxHist); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = ((GetEventHistogram()->Integral("width")/(fNEvents+0.)) * 1E-38 / (TotalIntegratedFlux())); // Setup Histograms SetHistograms(); // Final setup --------------------------------------------------- FinaliseMeasurement(); - + fSaveFine = false; }; bool T2K_NuMu_CC0pi_OC_XSec_2DPcos::isSignal(FitEvent *event){ return SignalDef::isCC0pi(event, 14, EnuMin, EnuMax); }; void T2K_NuMu_CC0pi_OC_XSec_2DPcos::FillEventVariables(FitEvent* event){ if (event->NumFSParticle(13) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; double pmu = Pmu.Vect().Mag()/1000.; double CosThetaMu = cos(Pnu.Vect().Angle(Pmu.Vect())); fXVar = pmu; fYVar = CosThetaMu; return; }; void T2K_NuMu_CC0pi_OC_XSec_2DPcos::FillHistograms(){ Measurement1D::FillHistograms(); if (Signal){ FillMCSlice( fXVar, fYVar, Weight ); } } void T2K_NuMu_CC0pi_OC_XSec_2DPcos::ConvertEventRates(){ for (size_t i = 0; i < nangbins; i++){ - if(Target=="O") fMCHistNuMuO_Slices[i]->GetSumw2(); - else if(Target=="C") fMCHistNuMuC_Slices[i]->GetSumw2(); + fMCHist_Slices[i]->GetSumw2(); } // Do standard conversion. Measurement1D::ConvertEventRates(); // Scale MC slices by their bin width for (size_t i = 0; i < nangbins; ++i) { - if(Target=="O") fMCHistNuMuO_Slices[i]->Scale(1. / (angular_binning_costheta[i + 1] - angular_binning_costheta[i])); - else if(Target=="C") fMCHistNuMuC_Slices[i]->Scale(1. / (angular_binning_costheta[i + 1] - angular_binning_costheta[i])); + fMCHist_Slices[i]->Scale(1. / (angular_binning_costheta[i + 1] - angular_binning_costheta[i])); } // Now Convert into 1D histogram fMCHist->Reset(); int bincount = 0; for (size_t i = 0; i < nangbins; i++){ - if(Target=="O"){ - for (int j = 0; j < fMCHistNuMuO_Slices[i]->GetNbinsX(); j++){ - fMCHist->SetBinContent(bincount+1, fMCHistNuMuO_Slices[i]->GetBinContent(j+1)); - bincount++; - } - } - else if(Target=="C"){ - for (int j = 0; j < fMCHistNuMuC_Slices[i]->GetNbinsX(); j++){ - fMCHist->SetBinContent(bincount+1, fMCHistNuMuC_Slices[i]->GetBinContent(j+1)); - bincount++; - } + for (int j = 0; j < fMCHist_Slices[i]->GetNbinsX(); j++){ + fMCHist->SetBinContent(bincount+1, fMCHist_Slices[i]->GetBinContent(j+1)); + bincount++; } } return; } void T2K_NuMu_CC0pi_OC_XSec_2DPcos::FillMCSlice(double x, double y, double w){ for (size_t i = 0; i < nangbins; ++i) { if ((y > angular_binning_costheta[i]) && (y <= angular_binning_costheta[i + 1])) { - if(Target=="O") fMCHistNuMuO_Slices[i]->Fill(x, w); - else if(Target=="C") fMCHistNuMuC_Slices[i]->Fill(x, w); + fMCHist_Slices[i]->Fill(x, w); } } } void T2K_NuMu_CC0pi_OC_XSec_2DPcos::SetHistograms(){ // Read covariance matrix fInputFileCov = new TFile( (FitPar::GetDataBase() + "/T2K/CC0pi/JointO-C/covmatrix_noreg.root").c_str(),"READ"); - TH1D* hLinearResult; - int Nbins; - - if(Target=="O"){ - fInputFile = new TFile( (FitPar::GetDataBase() + "/T2K/CC0pi/JointO-C/linear_unreg_results_O_nuisance.root").c_str(),"READ"); - // Read 1D data histogram - hLinearResult = (TH1D*) fInputFile->Get("LinResult"); - - Nbins = hLinearResult->GetNbinsX(); - - // Make covariance matrix - fFullCovar = new TMatrixDSym(Nbins); - TMatrixDSym* tempcov = (TMatrixDSym*) fInputFileCov->Get("covmatrixObin"); - - for(int ibin=0; ibinSetBinContent(bin+1, hLinearResult->GetBinContent(bin+1)); - } + fInputFile = new TFile(input_file_name.c_str(), "READ"); + TH1D* hLinearResult = (TH1D*) fInputFile->Get("LinResult"); + int Nbins = hLinearResult->GetNbinsX(); - fDataHist->Reset(); - int bincount = 0; - for (size_t i = 0; i < nangbins; i++){ - // Make slices for data - fDataHistNuMuO_Slices.push_back((TH1D*) fInputFile->Get(Form("dataslice_%zu",i))->Clone()); - fDataHistNuMuO_Slices[i]->SetNameTitle(Form("T2K_NuMu_CC0pi_O_2DPcos_data_Slice%zu",i), - (Form("T2K_NuMu_CC0pi_O_2DPcos_data_Slice%zu",i))); - fDataHistNuMuO_Slices[i]->Scale(1E-39); - // Loop over nbins and set errors from covariance - for (int j = 0; j < fDataHistNuMuO_Slices[i]->GetNbinsX(); j++){ - fDataHistNuMuO_Slices[i]->SetBinError(j+1, sqrt((*fFullCovar)(bincount,bincount))*1E-38); - - fDataHist->SetBinContent(bincount+1, fDataHistNuMuO_Slices[i]->GetBinContent(j+1)); - fDataHist->SetBinError(bincount+1, fDataHistNuMuO_Slices[i]->GetBinError(j+1)); - bincount++; - } - - // Save MC slices - fMCHistNuMuO_Slices.push_back((TH1D*) fDataHistNuMuO_Slices[i]->Clone()); - fMCHistNuMuO_Slices[i]->SetNameTitle(Form("T2K_NuMu_CC0pi_O_2DPcos_MC_Slice%zu",i), (Form("T2K_NuMu_CC0pi_O_2DPcos_MC_Slice%zu",i))); - - SetAutoProcessTH1(fDataHistNuMuO_Slices[i],kCMD_Write); - SetAutoProcessTH1(fMCHistNuMuO_Slices[i]); - } - } - else if(Target=="C"){ - - fInputFile = new TFile( (FitPar::GetDataBase() + "/T2K/CC0pi/JointO-C/linear_unreg_results_C_nuisance.root").c_str(),"READ"); - - // Read 1D data histogram - hLinearResult = (TH1D*) fInputFile->Get("LinResult"); - - Nbins = hLinearResult->GetNbinsX(); - - // Make the covariance matrix - fFullCovar = new TMatrixDSym(Nbins); - TMatrixDSym* tempcov = (TMatrixDSym*) fInputFileCov->Get("covmatrixCbin"); - - for(int ibin=0; ibinGet(covar_name.c_str()); - // Now Convert into 1D histrogram - fDataHist = new TH1D("LinearResultCarbon","LinearResultCarbon",Nbins,0,Nbins); - for (int bin = 0; bin < Nbins; bin++){ - fDataHist->SetBinContent(bin+1, hLinearResult->GetBinContent(bin+1)); + for(int ibin=0; ibinSetBinContent(bin+1, hLinearResult->GetBinContent(bin+1)); + } - fDataHist->Reset(); - - int bincount=0; - for (size_t i = 0; i < nangbins; i++){ - // Make slices for data - fDataHistNuMuC_Slices.push_back((TH1D*) fInputFile->Get(Form("dataslice_%zu",i))->Clone()); - fDataHistNuMuC_Slices[i]->SetNameTitle(Form("T2K_NuMu_CC0pi_C_2DPcos_data_Slice%zu",i), - (Form("T2K_NuMu_CC0pi_C_2DPcos_data_Slice%zu",i))); - - fDataHistNuMuC_Slices[i]->Scale(1E-39); - //Loop over nbins and set errors from covariance - for (int j = 0; j < fDataHistNuMuC_Slices[i]->GetNbinsX(); j++){ - fDataHistNuMuC_Slices[i]->SetBinError(j+1, sqrt((*fFullCovar)(bincount,bincount))*1E-38); - - fDataHist->SetBinContent(bincount+1, fDataHistNuMuC_Slices[i]->GetBinContent(j+1)); - fDataHist->SetBinError(bincount+1, fDataHistNuMuC_Slices[i]->GetBinError(j+1)); - bincount++; - } - - //Save MC slices - fMCHistNuMuC_Slices.push_back((TH1D*) fDataHistNuMuC_Slices[i]->Clone()); - fMCHistNuMuC_Slices[i]->SetNameTitle(Form("T2K_NuMu_CC0pi_C_2DPcos_MC_Slice%zu",i), (Form("T2K_NuMu_CC0pi_C_2DPcos_MC_Slice%zu",i))); - - SetAutoProcessTH1(fDataHistNuMuC_Slices[i],kCMD_Write); - SetAutoProcessTH1(fMCHistNuMuC_Slices[i]); + fDataHist->Reset(); + int bincount = 0; + for (size_t i = 0; i < nangbins; i++){ + // Make slices for data + fDataHist_Slices.push_back((TH1D*) fInputFile->Get(Form("dataslice_%zu",i))->Clone()); + fDataHist_Slices[i]->SetNameTitle(Form("%s_data_Slice%zu",name.c_str(), i), + (Form("%s_data_Slice%zu%s",name.c_str(),i,titles.c_str()))); + fDataHist_Slices[i]->Scale(1E-39); + // Loop over nbins and set errors from covariance + for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++){ + fDataHist_Slices[i]->SetBinError(j+1, sqrt((*fFullCovar)(bincount,bincount))*1E-38); + + fDataHist->SetBinContent(bincount+1, fDataHist_Slices[i]->GetBinContent(j+1)); + fDataHist->SetBinError(bincount+1, fDataHist_Slices[i]->GetBinError(j+1)); + bincount++; } - } - return; + // Save MC slices + fMCHist_Slices.push_back((TH1D*) fDataHist_Slices[i]->Clone()); + fMCHist_Slices[i]->SetNameTitle(Form("%s_MC_Slice%zu",name.c_str(),i), + Form("%s_MC_Slice%zu%s",name.c_str(),i,titles.c_str())); + SetAutoProcessTH1(fDataHist_Slices[i],kCMD_Write); + SetAutoProcessTH1(fMCHist_Slices[i]); + } + return; }; diff --git a/src/T2K/T2K_NuMu_CC0pi_OC_XSec_2DPcos.h b/src/T2K/T2K_NuMu_CC0pi_OC_XSec_2DPcos.h index 4fb080a..fe630d2 100644 --- a/src/T2K/T2K_NuMu_CC0pi_OC_XSec_2DPcos.h +++ b/src/T2K/T2K_NuMu_CC0pi_OC_XSec_2DPcos.h @@ -1,68 +1,65 @@ // Copyright 2016-2021 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 T2K_NUMU_CC0PI_OC_XSEC_2DPCOS_H_SEEN #define T2K_NUMU_CC0PI_OC_XSEC_2DPCOS_H_SEEN #include "Measurement1D.h" #include "TH2Poly.h" class T2K_NuMu_CC0pi_OC_XSec_2DPcos : public Measurement1D { public: /// Basic Constructor. /// /brief Parses two different measurements. /// T2K_NuMu_CC0pi_OC_XSec_2DPcos(nuiskey samplekey); /// Virtual Destructor ~T2K_NuMu_CC0pi_OC_XSec_2DPcos() {}; /// Signal definition bool isSignal(FitEvent *nvect); /// Read histograms void SetHistograms(); /// Bin Tmu CosThetaMu void FillEventVariables(FitEvent* customEvent); // Fill Histograms void FillHistograms(); /// Event scaling void ConvertEventRates(); private: TFile* fInputFile; TFile* fInputFileCov; - std::vector fMCHistNuMuO_Slices; - std::vector fDataHistNuMuO_Slices; - std::vector fMCHistNuMuC_Slices; - std::vector fDataHistNuMuC_Slices; + std::vector fMCHist_Slices; + std::vector fDataHist_Slices; double pmu, CosThetaMu; std::string Target; void FillMCSlice(double x, double y, double w); - }; #endif diff --git a/src/T2K/T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint.cxx b/src/T2K/T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint.cxx index 1ad9e6a..1209f32 100644 --- a/src/T2K/T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint.cxx +++ b/src/T2K/T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint.cxx @@ -1,161 +1,162 @@ #include "T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint.h" //******************************************************************** T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint::T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint(nuiskey samplekey){ //******************************************************************** fSettings = LoadSampleSettings(samplekey); std::string descrip = "T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint. \n" "Target: O-C \n" "Flux: T2K 2.5 degree off-axis (ND280) \n" "Signal: CC0pi\n" "DOI:10.1103/PhysRevD.101.112004"; fSettings.SetTitle("T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint"); fSettings.DefineAllowedSpecies("numu"); fSettings.SetDescription(descrip); fSettings.SetXTitle("p_{#mu}-cos#theta_{#mu}"); - fSettings.SetYTitle("d^{2}#sigma/dP_{#mu}dcos#theta_{#mu} (cm^{2}/GeV)"); + fSettings.SetYTitle("d^{2}#sigma/dp_{#mu}dcos#theta_{#mu} (cm^{2}/GeV)"); fSettings.SetAllowedTypes("DIAG,FULL/FREE,SHAPE,FIX/SYSTCOV/STATCOV","FIX"); fSettings.SetEnuRange(0.0, 30.0); fSettings.DefineAllowedTargets("O,C"); FinaliseSampleSettings(); if (fSubInFiles.size() != 2) { NUIS_ABORT("T2K NuMuCC0pi O-C joint requires input files in format: NuMuCC0pi on O and NuMuCC0pi on C"); } std::string inFileNuMuO = fSubInFiles.at(0); std::string inFileNuMuC = fSubInFiles.at(1); // Create some config keys nuiskey NuMuCC0piOKey = Config::CreateKey("sample"); NuMuCC0piOKey.SetS("input", inFileNuMuO); NuMuCC0piOKey.SetS("type", fSettings.GetS("type")); NuMuCC0piOKey.SetS("name", "T2K_NuMu_CC0pi_O_XSec_2DPcos"); NuMuCC0piO = new T2K_NuMu_CC0pi_OC_XSec_2DPcos(NuMuCC0piOKey); nuiskey NuMuCC0piCKey = Config::CreateKey("sample"); NuMuCC0piCKey.SetS("input", inFileNuMuC); NuMuCC0piCKey.SetS("type", fSettings.GetS("type")); NuMuCC0piCKey.SetS("name", "T2K_NuMu_CC0pi_C_XSec_2DPcos"); NuMuCC0piC = new T2K_NuMu_CC0pi_OC_XSec_2DPcos(NuMuCC0piCKey); // Sort out the data hist this->CombineDataHists(); // Set the covariance SetCovariance(); // Add to chain for processing fSubChain.clear(); fSubChain.push_back(NuMuCC0piO); fSubChain.push_back(NuMuCC0piC); // This saves information from the sub-measurements fSaveSubMeas = true; + fSaveFine = false; FinaliseMeasurement(); }; //******************************************************************** void T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint::SetCovariance(){ //******************************************************************** // Read covariance matrix fInputFileCov = new TFile( (FitPar::GetDataBase() + "/T2K/CC0pi/JointO-C/covmatrix_noreg.root").c_str(),"READ"); TMatrixDSym* tempcov = (TMatrixDSym*) fInputFileCov->Get("covmatrixOeCbin"); fFullCovar = new TMatrixDSym(tempcov->GetNrows()); for(int ibin=0; ibinGetNrows(); ibin++) { for(int jbin=0; jbinGetNrows(); jbin++) { // The factor 1E-2 needed since the covariance matrix in the // data release is divided by 1E-78 (*fFullCovar)(ibin,jbin) = (*tempcov)(ibin,jbin)*1E-2; } } covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); return; } //******************************************************************** void T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint::CombineDataHists(){ //******************************************************************** TH1D *hNuMuOData = (TH1D*)NuMuCC0piO->GetDataHistogram(); TH1D *hNuMuCData = (TH1D*)NuMuCC0piC->GetDataHistogram(); int nbins = hNuMuOData->GetNbinsX() + hNuMuCData->GetNbinsX(); fDataHist = new TH1D((fSettings.GetName() + "_data").c_str(), (fSettings.GetFullTitles()).c_str(), nbins, 0, nbins); fDataHist->SetDirectory(0); int count = 0; for (int x=0; xGetNbinsX(); ++x){ fDataHist->SetBinContent(count+1, hNuMuOData->GetBinContent(x+1)); fDataHist->SetBinError(count+1, hNuMuOData->GetBinError(x+1)); fDataHist->GetXaxis()->SetBinLabel(count+1, Form("NuMu CC0pi %.1f-%.1f", hNuMuOData->GetXaxis()->GetBinLowEdge(x+1), hNuMuOData->GetXaxis()->GetBinUpEdge(x+1))); count++; } for (int x=0; xGetNbinsX(); ++x){ fDataHist->SetBinContent(count+1, hNuMuCData->GetBinContent(x+1)); fDataHist->SetBinError(count+1, hNuMuCData->GetBinError(x+1)); fDataHist->GetXaxis()->SetBinLabel(count+1, Form("AntiNuMu CC0pi %.1f-%.1f", hNuMuCData->GetXaxis()->GetBinLowEdge(x+1), hNuMuCData->GetXaxis()->GetBinUpEdge(x+1))); count++; } } //******************************************************************** void T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint::SetHistograms() { //******************************************************************** NuMuCC0piO->SetHistograms(); NuMuCC0piC->SetHistograms(); return; } //******************************************************************** void T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint::FillHistograms() { //******************************************************************** NuMuCC0piO->FillHistograms(); NuMuCC0piC->FillHistograms(); return; } //******************************************************************** void T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint::ConvertEventRates() { //******************************************************************** NuMuCC0piO->ConvertEventRates(); NuMuCC0piC->ConvertEventRates(); TH1D* hNuMuCC0piO = (TH1D*)NuMuCC0piO->GetMCHistogram(); TH1D* hNuMuCC0piC = (TH1D*)NuMuCC0piC->GetMCHistogram(); int count = 0; for (int i = 0; i < hNuMuCC0piO->GetNbinsX(); ++i) { fMCHist->SetBinContent(count + 1, hNuMuCC0piO->GetBinContent(i + 1)); fMCHist->SetBinError(count + 1, hNuMuCC0piO->GetBinError(i + 1)); count++; } for (int i = 0; i < hNuMuCC0piC->GetNbinsX(); ++i) { fMCHist->SetBinContent(count + 1, hNuMuCC0piC->GetBinContent(i + 1)); fMCHist->SetBinError(count + 1, hNuMuCC0piC->GetBinError(i + 1)); count++; } return; }