diff --git a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx index a403098..a4818f0 100644 --- a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx +++ b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx @@ -1,373 +1,378 @@ // 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_SignalDef.h" #include "T2K_CC0pi_XSec_2DPcos_nu_I.h" static size_t nangbins = 9; static double angular_binning_costheta[] = {-1, 0, 0.6, 0.7, 0.8, 0.85, 0.9, 0.94, 0.98, 1}; //******************************************************************** T2K_CC0pi_XSec_2DPcos_nu_I::T2K_CC0pi_XSec_2DPcos_nu_I(nuiskey samplekey) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "T2K_CC0pi_XSec_2DPcos_nu_I sample. \n" "Target: CH \n" "Flux: T2K 2.5 degree off-axis (ND280) \n" "Signal: CC0pi\n" "https://journals.aps.org/prd/abstract/10.1103/" "PhysRevD.93.112012 Analysis I"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); - fSettings.SetXTitle("P_{#mu} (GeV)"); - fSettings.SetYTitle("cos#theta_{#mu}"); - fSettings.SetZTitle("d^{2}#sigma/dP_{#mu}dcos#theta_{#mu} (cm^{2}/GeV)"); + fSettings.SetXTitle("cos#theta_{#mu}-p_{#mu} (GeV)"); + fSettings.SetYTitle("d^{2}#sigma/dp_{#mu}dcos#theta_{#mu} (cm^{2}/GeV)"); fSettings.SetAllowedTypes("FULL,DIAG/FREE,SHAPE,FIX/SYSTCOV/STATCOV", "FIX/FULL"); fSettings.DefineAllowedTargets("C,H"); fSettings.SetEnuRangeFromFlux(fFluxHist); // CCQELike plot information fSettings.SetTitle("T2K_CC0pi_XSec_2DPcos_nu_I"); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); fMaskMomOverflow = false; if (samplekey.Has("mask_mom_overflow")) { fMaskMomOverflow = samplekey.GetB("mask_mom_overflow"); } // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = ((GetEventHistogram()->Integral("width") / (fNEvents + 0.)) * 1E-38 / (TotalIntegratedFlux())); assert(std::isnormal(fScaleFactor)); // Setup Histograms SetHistograms(); // Final setup --------------------------------------------------- FinaliseMeasurement(); + fSaveFine = false; }; bool T2K_CC0pi_XSec_2DPcos_nu_I::isSignal(FitEvent *event) { return SignalDef::isT2K_CC0pi(event, EnuMin, EnuMax, SignalDef::kAnalysis_I); }; void T2K_CC0pi_XSec_2DPcos_nu_I::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; fMode = event->Mode; return; }; void T2K_CC0pi_XSec_2DPcos_nu_I::FillHistograms() { Measurement1D::FillHistograms(); if (Signal) { fMCHist_Fine2D->Fill(fXVar, fYVar, Weight); FillMCSlice(fXVar, fYVar, Weight); } } // Modification is needed after the full reconfigure to move bins around // Otherwise this would need to be replaced by a TH2Poly which is too awkward. void T2K_CC0pi_XSec_2DPcos_nu_I::ConvertEventRates() { // Do standard conversion. Measurement1D::ConvertEventRates(); // Scale MC slices by their bin area for (size_t i = 0; i < nangbins; ++i) { fMCHist_Slices[i]->Scale(fScaleFactor / (angular_binning_costheta[i + 1] - angular_binning_costheta[i]), "width"); } // Now Convert into 1D list fMCHist->Reset(); int bincount = 0; for (size_t i = 0; i < nangbins; i++) { for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) { fMCHist->SetBinContent(bincount + 1, fMCHist_Slices[i]->GetBinContent(j + 1)); fMCHist->SetBinError(bincount + 1, fMCHist_Slices[i]->GetBinError(j + 1)); bincount++; } } if (fMCHist_Modes) { fMCHist_Modes->Reset(); for (std::map >::iterator it = fMCModeHists_Slices.begin(); it != fMCModeHists_Slices.end(); ++it) { bincount = 0; for (size_t i = 0; i < nangbins; i++) { it->second[i]->Scale(fScaleFactor / (angular_binning_costheta[i + 1] - angular_binning_costheta[i]), "width"); for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) { fMCHist_Modes->SetBinContent(it->first, bincount + 1, 0, 0, it->second[i]->GetBinContent(j + 1)); fMCHist_Modes->SetBinError(it->first, bincount + 1, 0, 0, it->second[i]->GetBinError(j + 1)); bincount++; } } } } return; } void T2K_CC0pi_XSec_2DPcos_nu_I::FillMCSlice(double x, double y, double w) { if (fMCHist_Modes && !fMCModeHists_Slices.count(fMode)) { for (size_t i = 0; i < nangbins; ++i) { std::stringstream ss(""); ss << "T2K_CC0pi_XSec_2DPcos_nu_I_MODE_ " << fMode << "_slice" << i << std::endl; fMCModeHists_Slices[fMode].push_back( static_cast(fMCHist_Slices[i]->Clone(ss.str().c_str()))); fMCModeHists_Slices[fMode].back()->Reset(); } } for (size_t i = 0; i < nangbins; ++i) { if ((y >= angular_binning_costheta[i]) && (y < angular_binning_costheta[i + 1])) { fMCHist_Slices[i]->Fill(x, w); if (fMCHist_Modes) { fMCModeHists_Slices[fMode][i]->Fill(x, w); } } } } void T2K_CC0pi_XSec_2DPcos_nu_I::SetHistograms() { // Read in 1D Data Histograms TFile input( (FitPar::GetDataBase() + "/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root") - .c_str(), - "READ"); + .c_str(), "READ"); fMCHist_Fine2D = new TH2D("T2K_CC0pi_XSec_2DPcos_nu_I_Fine2D", - "T2K_CC0pi_XSec_2DPcos_nu_I_Fine2D", 400, 0.0, 30.0, - 100, -1.0, 1.0); + "T2K_CC0pi_XSec_2DPcos_nu_I_Fine2D", + 400, 0.0, 30.0, 100, -1.0, 1.0); + fMCHist_Fine2D->GetXaxis()->SetTitle("p_{#mu} (GeV)"); + fMCHist_Fine2D->GetYaxis()->SetTitle("cos#theta_{#mu}"); + fMCHist_Fine2D->GetZaxis()->SetTitle(fSettings.GetYTitle().c_str()); + fMCHist_Fine2D->SetDirectory(NULL); SetAutoProcessTH1(fMCHist_Fine2D); TH2D *tempcov = (TH2D *)input.Get("analysis1_totcov"); fFullCovar = new TMatrixDSym(tempcov->GetNbinsX()); for (int i = 0; i < tempcov->GetNbinsX(); i++) { for (int j = 0; j < tempcov->GetNbinsX(); j++) { (*fFullCovar)(i, j) = tempcov->GetBinContent(i + 1, j + 1); } } // Read in 2D Data Slices and Make MC Slices int bincount = 0; for (size_t i = 0; i < nangbins; i++) { fDataHist_Slices.push_back( (TH1D *)input.Get(Form("dataslice_%lu", i))->Clone()); fDataHist_Slices[i]->SetNameTitle( Form("T2K_CC0pi_XSec_2DPcos_nu_I_data_Slice%lu", i), (Form("T2K_CC0pi_XSec_2DPcos_nu_I_data_Slice%lu, cos(#theta) [%f,%f] ", i, angular_binning_costheta[i], angular_binning_costheta[i + 1]))); fDataHist_Slices.back()->SetDirectory(NULL); + fDataHist_Slices.back()->GetXaxis()->SetTitle("p_{#mu} (GeV)"); + fDataHist_Slices.back()->GetYaxis()->SetTitle(fSettings.GetYTitle().c_str()); // Loop over nbins and set errors from covar for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) { fDataHist_Slices[i]->SetBinError( j + 1, sqrt((*fFullCovar)(bincount, bincount)) * 1E-38); bincount++; } } assert(bincount == tempcov->GetNbinsX()); if (fMaskMomOverflow) { MaskMomOverflow(); bincount = fFullCovar->GetNcols(); } std::vector > data_slice_bcbes; for (size_t i = 0; i < nangbins; i++) { for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) { data_slice_bcbes.push_back( std::make_pair(fDataHist_Slices[i]->GetBinContent(j + 1), fDataHist_Slices[i]->GetBinError(j + 1))); } } for (size_t i = 0; i < nangbins; i++) { fMCHist_Slices.push_back((TH1D *)fDataHist_Slices[i]->Clone()); fMCHist_Slices.back()->SetDirectory(NULL); fMCHist_Slices.back()->Reset(); fMCHist_Slices.back()->SetLineColor(kRed); fMCHist_Slices[i]->SetNameTitle( Form("T2K_CC0pi_XSec_2DPcos_nu_I_MC_Slice%lu", i), (Form("T2K_CC0pi_XSec_2DPcos_nu_I_MC_Slice%lu, cos(#theta) [%f,%f] ", i, angular_binning_costheta[i], angular_binning_costheta[i + 1]))); } covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); - fDataHist = - new TH1D("T2K_CC0pi_XSec_2DPcos_nu_I_DATA_1D", - "T2K_CC0pi_XSec_2DPcos_nu_I_DATA_1D", bincount, 0, bincount); + fDataHist = new TH1D("T2K_CC0pi_XSec_2DPcos_nu_I_data_1D", + ("T2K_CC0pi_XSec_2DPcos_nu_I_data_1D"+fSettings.PlotTitles()).c_str(), + bincount, 0, bincount); fDataHist->SetDirectory(NULL); for (size_t i = 0; i < data_slice_bcbes.size(); ++i) { fDataHist->SetBinContent(i + 1, data_slice_bcbes[i].first); fDataHist->SetBinError(i + 1, data_slice_bcbes[i].second); } SetShapeCovar(); fMCHist = (TH1D *)fDataHist->Clone("T2K_CC0pi_XSec_2DPcos_nu_I_MC_1D"); fMCHist->SetDirectory(NULL); return; } void T2K_CC0pi_XSec_2DPcos_nu_I::MaskMomOverflow() { std::vector bins_to_cut; size_t nallbins = 0; for (size_t i = 0; i < nangbins; i++) { std::vector slice_bin_edges; slice_bin_edges.push_back( fDataHist_Slices[i]->GetXaxis()->GetBinLowEdge(1)); for (int j = 0; j < (fDataHist_Slices[i]->GetNbinsX() - 1); j++) { slice_bin_edges.push_back( fDataHist_Slices[i]->GetXaxis()->GetBinUpEdge(j + 1)); nallbins++; } bins_to_cut.push_back(nallbins++); TH1D *tmp = new TH1D(fDataHist_Slices[i]->GetName(), fDataHist_Slices[i]->GetTitle(), slice_bin_edges.size() - 1, slice_bin_edges.data()); tmp->SetDirectory(NULL); for (int j = 0; j < (fDataHist_Slices[i]->GetNbinsX() - 1); j++) { tmp->SetBinContent(j + 1, fDataHist_Slices[i]->GetBinContent(j + 1)); tmp->SetBinError(j + 1, fDataHist_Slices[i]->GetBinError(j + 1)); } delete fDataHist_Slices[i]; fDataHist_Slices[i] = tmp; } TMatrixDSym *tmpcovar = new TMatrixDSym(nallbins - bins_to_cut.size()); int icut = 0; for (int ifull = 0; ifull < fFullCovar->GetNcols(); ifull++) { if (std::find(bins_to_cut.begin(), bins_to_cut.end(), ifull) != bins_to_cut.end()) { continue; } int jcut = 0; for (int jfull = 0; jfull < fFullCovar->GetNcols(); jfull++) { if (std::find(bins_to_cut.begin(), bins_to_cut.end(), jfull) != bins_to_cut.end()) { continue; } (*tmpcovar)[icut][jcut] = (*fFullCovar)[ifull][jfull]; jcut++; } icut++; } delete fFullCovar; fFullCovar = tmpcovar; } void T2K_CC0pi_XSec_2DPcos_nu_I::Write(std::string drawopt) { this->Measurement1D::Write(drawopt); for (size_t i = 0; i < nangbins; i++) { fMCHist_Slices[i]->Write(); fDataHist_Slices[i]->Write(); } if (fResidualHist) { std::vector MCResidual_Slices; size_t tb_it = 0; for (size_t i = 0; i < fMCHist_Slices.size(); ++i) { std::string name = Form("T2K_CC0pi_XSec_2DPcos_nu_I_RESIDUAL_Slice%lu", i); MCResidual_Slices.push_back( static_cast(fMCHist_Slices[i]->Clone(name.c_str()))); MCResidual_Slices.back()->Reset(); for (int j = 0; j < fMCHist_Slices[i]->GetXaxis()->GetNbins(); ++j) { double bc = fResidualHist->GetBinContent(tb_it + 1); MCResidual_Slices.back()->SetBinContent(j + 1, bc); tb_it++; } MCResidual_Slices.back()->Write(); } } if (fChi2LessBinHist) { std::vector MCChi2LessBin_Slices; size_t tb_it = 0; for (size_t i = 0; i < fMCHist_Slices.size(); ++i) { std::string name = Form("T2K_CC0pi_XSec_2DPcos_nu_I_Chi2NMinusOne_Slice%lu", i); MCChi2LessBin_Slices.push_back( static_cast(fMCHist_Slices[i]->Clone(name.c_str()))); MCChi2LessBin_Slices.back()->Reset(); for (int j = 0; j < fMCHist_Slices[i]->GetXaxis()->GetNbins(); ++j) { double bc = fChi2LessBinHist->GetBinContent(tb_it + 1); MCChi2LessBin_Slices.back()->SetBinContent(j + 1, bc); tb_it++; } MCChi2LessBin_Slices.back()->Write(); } } } diff --git a/src/T2K/T2K_CC0pi_XSec_H2O_2DPcos_anu.cxx b/src/T2K/T2K_CC0pi_XSec_H2O_2DPcos_anu.cxx index bb7dfd5..5a3ba37 100755 --- a/src/T2K/T2K_CC0pi_XSec_H2O_2DPcos_anu.cxx +++ b/src/T2K/T2K_CC0pi_XSec_H2O_2DPcos_anu.cxx @@ -1,190 +1,199 @@ // 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_SignalDef.h" #include "T2K_CC0pi_XSec_H2O_2DPcos_anu.h" static size_t nmombins = 7; static double mom_binning[] = { 400., 530., 670., 800., 1000., 1380., 2010., 3410. }; static int ncosbins[] = { 2, 3, 3, 3, 3, 3, 2 }; static double costheta_binning[][7] = { { 0.84, 0.94, 1. }, { 0.85, 0.92, 0.96, 1. }, { 0.88, 0.93, 0.97, 1. }, { 0.90, 0.94, 0.97, 1. }, { 0.91, 0.95, 0.97, 1. }, { 0.92, 0.96, 0.98, 1. }, { 0.95, 0.98, 1. } }; //******************************************************************** T2K_CC0pi_XSec_H2O_2DPcos_anu::T2K_CC0pi_XSec_H2O_2DPcos_anu(nuiskey samplekey) { //******************************************************************** + // Don't use the FINE histogram + fSaveFine = false; + // Sample overview --------------------------------------------------- std::string descrip = "T2K_CC0pi_XSec_H2O_2DPcos_anu. \n" "Target: H2O \n" "Flux: T2K 2.5 degree off-axis (ND280-P0D) \n" "Signal: CC0pi\n" "arXiv:1908.10249"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); - fSettings.SetXTitle("cos#theta_{#mu}-p_{#mu}"); + fSettings.SetXTitle("cos#theta_{#mu}-p_{#mu} (GeV)"); fSettings.SetYTitle("d^{2}#sigma/dp_{#mu}dcos#theta_{#mu} (cm^{2}/GeV)"); fSettings.SetAllowedTypes("FULL,DIAG/FREE,SHAPE,FIX/SYSTCOV/STATCOV","FIX/FULL"); fSettings.SetEnuRange(0.0, 10.0); fSettings.DefineAllowedTargets("H,O"); // CCQELike plot information fSettings.SetTitle("T2K_CC0pi_XSec_H2O_2DPcos_anu"); fSettings.DefineAllowedSpecies("numub"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/water molecule fScaleFactor = ((GetEventHistogram()->Integral("weight")/(fNEvents+0.)) * 18E-38 / (TotalIntegratedFlux(""))); // Setup Histograms SetHistograms(); // Final setup FinaliseMeasurement(); }; bool T2K_CC0pi_XSec_H2O_2DPcos_anu::isSignal(FitEvent *event){ return SignalDef::isT2K_CC0piAnuP0D(event, EnuMin, EnuMax); }; void T2K_CC0pi_XSec_H2O_2DPcos_anu::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(); double CosThetaMu = cos(Pnu.Vect().Angle(Pmu.Vect())); fXVar = pmu; fYVar = CosThetaMu; return; }; void T2K_CC0pi_XSec_H2O_2DPcos_anu::FillHistograms(){ Measurement1D::FillHistograms(); if (Signal){ FillMCSlice( fXVar, fYVar, Weight ); } } // Modification is needed after the full reconfigure to move bins around // Otherwise this would need to be replaced by a TH2Poly which is too awkward. void T2K_CC0pi_XSec_H2O_2DPcos_anu::ConvertEventRates(){ for (size_t i = 0; i < nmombins; i++){ fMCHist_Slices[i]->GetSumw2(); } // Do standard conversion. Measurement1D::ConvertEventRates(); for (size_t i = 0; i < nmombins; ++i) { fMCHist_Slices[i]->Scale(1. / (mom_binning[i + 1] - mom_binning[i])); } // Now Convert into 1D list fMCHist->Reset(); int bincount = 0; for (size_t i = 0; i < nmombins; i++){ for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++){ fMCHist->SetBinContent(bincount+1, fMCHist_Slices[i]->GetBinContent(j+1)); bincount++; } } return; } void T2K_CC0pi_XSec_H2O_2DPcos_anu::FillMCSlice(double x, double y, double w){ for (size_t i = 0; i < nmombins; ++i) { if ((x > mom_binning[i]) && (x <= mom_binning[i + 1])) { fMCHist_Slices[i]->Fill(y, w); } } } void T2K_CC0pi_XSec_H2O_2DPcos_anu::SetHistograms(){ // Read in 1D Data Histograms fInputFile = new TFile( (FitPar::GetDataBase() + "/T2K/CC0pi/AntiNuMuH2O/AntiNuMuOnH2O_unreg.root").c_str(),"READ"); + fPlotTitles = fSettings.GetFullTitles(); // Read in 1D Data fDataHist = (TH1D*) fInputFile->Get("xsecDataRelease"); fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_data" + fPlotTitles).c_str()); int Nbins = fDataHist->GetNbinsX(); // Read relative covariance matrix TH2D* tempcov = (TH2D*) fInputFile->Get("covDataRelease"); // Make absolute covariance matrix fFullCovar = new TMatrixDSym(Nbins); for (int i = 0; i < Nbins; i++){ for (int j = 0; j < Nbins; j++){ (*fFullCovar)(i,j) = tempcov->GetBinContent(i+1, j+1)*fDataHist->GetBinContent(i+1)*fDataHist->GetBinContent(j+1)*1E76; } } covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); // Read in 2D Data Slices and Make MC Slices int bincount = 0; for (uint i = 0; i < nmombins; ++i) { - fDataHist_Slices.push_back(new TH1D(Form("T2K_CC0pi_XSec_H2O_2DPcos_anu_data_Slice%u",i),Form("T2K_CC0pi_XSec_H2O_2DPcos_anu_data_Slice%u",i),ncosbins[i],costheta_binning[i])); + fDataHist_Slices.push_back(new TH1D(Form("T2K_CC0pi_XSec_H2O_2DPcos_anu_data_Slice%u",i), + Form("T2K_CC0pi_XSec_H2O_2DPcos_anu_data_Slice%u",i), + ncosbins[i],costheta_binning[i])); + fDataHist_Slices[i]->GetXaxis()->SetTitle("cos#theta_{#mu}"); + fDataHist_Slices[i]->GetYaxis()->SetTitle(fSettings.GetYTitle().c_str()); + for (int j = 0; j < ncosbins[i]; ++j) { fDataHist->SetBinError(bincount+1,sqrt((*fFullCovar)(bincount,bincount))*1E-38); fDataHist_Slices[i]->SetBinContent(j+1, fDataHist->GetBinContent(bincount+1)); fDataHist_Slices[i]->SetBinError(j+1, sqrt((*fFullCovar)(bincount,bincount))*1E-38); bincount++; } fMCHist_Slices.push_back((TH1D*) fDataHist_Slices[i]->Clone()); fMCHist_Slices[i]->SetNameTitle(Form("T2K_CC0pi_XSec_H2O_2DPcos_anu_MC_Slice%u",i), - (Form("T2K_CC0pi_XSec_H2O_2DPcos_anu_MC_Slice%u",i))); + (Form("T2K_CC0pi_XSec_H2O_2DPcos_anu_MC_Slice%u",i))); SetAutoProcessTH1(fDataHist_Slices[i],kCMD_Write); SetAutoProcessTH1(fMCHist_Slices[i]); } return; }; diff --git a/src/T2K/T2K_CCinc_XSec_2DPcos_nu_nonuniform.cxx b/src/T2K/T2K_CCinc_XSec_2DPcos_nu_nonuniform.cxx index 06fbefb..96f4f52 100644 --- a/src/T2K/T2K_CCinc_XSec_2DPcos_nu_nonuniform.cxx +++ b/src/T2K/T2K_CCinc_XSec_2DPcos_nu_nonuniform.cxx @@ -1,257 +1,254 @@ #include "T2K_CCinc_XSec_2DPcos_nu_nonuniform.h" #include "T2K_SignalDef.h" // *********************************** // Implemented by Alfonso Garcia, Barcelona (now NIKHEF) // Clarence Wret, Rochester // (Alfonso was the T2K analyser) // *********************************** //******************************************************************** T2K_CCinc_XSec_2DPcos_nu_nonuniform::T2K_CCinc_XSec_2DPcos_nu_nonuniform( nuiskey samplekey) { //******************************************************************** + fSaveFine = false; fAllowedTypes += "/GENIE/NEUT"; // Sample overview --------------------------------------------------- std::string descrip = "T2K_CCinc_XSec_2DPcos_nu_nonuniform sample. \n" "Target: CH \n" "Flux: T2K FHC numu \n" "Signal: CC-inclusive \n" "https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.012004"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); - fSettings.SetXTitle(""); - // fSettings.SetXTitle("p_{#mu} (GeV)"); - // fSettings.SetYTitle("cos#theta_{#mu}"); + fSettings.SetXTitle("p_{#mu}-cos#theta_{#mu}"); fSettings.SetYTitle("#frac{d^{2}#sigma}{dp_{#mu}dcos#theta_{#mu}} " - "[#frac{cm^{2}}{nucleon/GeV/c}]"); + "(cm^{2}/nucleon/GeV)"); fSettings.SetEnuRange(0.0, 30.0); fSettings.DefineAllowedTargets("C,H"); // CCQELike plot information fSettings.SetTitle("T2K CC-inclusive p_{#mu} cos#theta_{#mu}"); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = GetEventHistogram()->Integral("width") * 1E-38 / fNEvents / TotalIntegratedFlux(); // Default to using the NEUT unfolded data UnfoldWithGENIE = false; // Check option if (fSettings.Found("type", "GENIE")) UnfoldWithGENIE = true; // Tell user what's happening if (UnfoldWithGENIE) { NUIS_LOG(SAM, fName << " is using GENIE unfolded data. Want NEUT? Specify " "type=\"NEUT\" in your config file"); } else { NUIS_LOG(SAM, fName << " is using NEUT unfolded data. Want GENIE? Specify " "type=\"GENIE\" in your config file"); } // Setup Histograms SetHistograms(); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; // Signal is simply a CC inclusive without any angular/momentum cuts bool T2K_CCinc_XSec_2DPcos_nu_nonuniform::isSignal(FitEvent *event) { return SignalDef::isCCINC(event, 14, EnuMin, EnuMax); }; void T2K_CCinc_XSec_2DPcos_nu_nonuniform::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; }; // Fill up the MCSlice void T2K_CCinc_XSec_2DPcos_nu_nonuniform::FillHistograms() { if (Signal) FillMCSlice(fXVar, fYVar, Weight); } // Modification is needed after the full reconfigure to move bins around void T2K_CCinc_XSec_2DPcos_nu_nonuniform::ConvertEventRates() { // Do standard conversion. Measurement1D::ConvertEventRates(); // First scale MC slices also by their width in Y fMCHist_Slices[0]->Scale(1.0 / 0.25); fMCHist_Slices[1]->Scale(1.0 / 0.50); fMCHist_Slices[2]->Scale(1.0 / 0.20); fMCHist_Slices[3]->Scale(1.0 / 0.15); fMCHist_Slices[4]->Scale(1.0 / 0.11); fMCHist_Slices[5]->Scale(1.0 / 0.09); fMCHist_Slices[6]->Scale(1.0 / 0.07); fMCHist_Slices[7]->Scale(1.0 / 0.05); fMCHist_Slices[8]->Scale(1.0 / 0.04); fMCHist_Slices[9]->Scale(1.0 / 0.025); fMCHist_Slices[10]->Scale(1.0 / 0.015); // Now Convert into 1D list fMCHist->Reset(); int bincount = 0; for (int i = 0; i < nSlices; i++) { for (int j = 0; j < fMCHist_Slices[i]->GetNbinsX(); j++) { fMCHist->SetBinContent(bincount + 1, fMCHist_Slices[i]->GetBinContent(j + 1)); fMCHist->SetBinError(bincount + 1, fMCHist_Slices[i]->GetBinError(j + 1)); bincount++; } } }; void T2K_CCinc_XSec_2DPcos_nu_nonuniform::FillMCSlice(double x, double y, double w) { if (y >= -1.0 && y < -0.25) fMCHist_Slices[0]->Fill(x, w); else if (y >= -0.25 && y < 0.25) fMCHist_Slices[1]->Fill(x, w); else if (y >= 0.25 && y < 0.45) fMCHist_Slices[2]->Fill(x, w); else if (y >= 0.45 && y < 0.6) fMCHist_Slices[3]->Fill(x, w); else if (y >= 0.6 && y < 0.71) fMCHist_Slices[4]->Fill(x, w); else if (y >= 0.71 && y < 0.80) fMCHist_Slices[5]->Fill(x, w); else if (y >= 0.80 && y < 0.87) fMCHist_Slices[6]->Fill(x, w); else if (y >= 0.87 && y < 0.92) fMCHist_Slices[7]->Fill(x, w); else if (y >= 0.92 && y < 0.96) fMCHist_Slices[8]->Fill(x, w); else if (y >= 0.96 && y <= 0.985) fMCHist_Slices[9]->Fill(x, w); else if (y >= 0.985 && y <= 1.0) fMCHist_Slices[10]->Fill(x, w); }; void T2K_CCinc_XSec_2DPcos_nu_nonuniform::SetHistograms() { // Read in 1D Data Histograms - TFile *fInputFile = - new TFile((FitPar::GetDataBase() + + TFile *fInputFile = new TFile((FitPar::GetDataBase() + "T2K/CCinc/nd280data-numu-cc-inc-xs-on-c-2018/histograms.root") - .c_str(), - "OPEN"); + .c_str(), "OPEN"); // Number of theta slices in the release nSlices = 11; // Data release includes unfolding with NEUT or GENIE as prior // Choose whichever the user specifies std::string basename; if (UnfoldWithGENIE) basename = "hist_xsec_data_prior_neut_cthbin"; else basename = "hist_xsec_data_prior_neut_cthbin"; // Read in 2D Data Slices and Make MC Slices // Count the number of bins we have in total so we can match covariance matrix int bincount = 0; for (int i = 0; i < nSlices; i++) { // Get Data Histogram - // fDataHist_Slices.push_back((TH1D*)fInputFile->Get(Form("dataslice_%i",i))->Clone()); fDataHist_Slices.push_back( (TH1D *)fInputFile->Get(Form("%s%i", basename.c_str(), i)) ->Clone( - Form("T2K_CCinc_XSec_2DPcos_nu_nonuniform_slice%i_data", i))); + Form("T2K_CCinc_XSec_2DPcos_nu_nonuniform_Slice%i_data", i))); fDataHist_Slices[i]->SetDirectory(0); fDataHist_Slices[i]->Scale(1E-39); - fDataHist_Slices[i]->GetYaxis()->SetTitle(fSettings.GetS("ytitle").c_str()); + fDataHist_Slices[i]->GetXaxis()->SetTitle("p_{#mu} (GeV)"); + fDataHist_Slices[i]->GetYaxis()->SetTitle(fSettings.GetYTitle().c_str()); // Count up the bins bincount += fDataHist_Slices.back()->GetNbinsX(); // Make MC Clones fMCHist_Slices.push_back((TH1D *)fDataHist_Slices[i]->Clone( Form("T2K_CCinc_XSec_2DPcos_nu_nonuniform_Slice%i_MC", i))); fMCHist_Slices[i]->Reset(); fMCHist_Slices[i]->SetDirectory(0); fMCHist_Slices[i]->SetLineColor(kRed); - fMCHist_Slices[i]->GetYaxis()->SetTitle(fSettings.GetS("ytitle").c_str()); + fMCHist_Slices[i]->GetYaxis()->SetTitle(fSettings.GetYTitle().c_str()); SetAutoProcessTH1(fDataHist_Slices[i], kCMD_Write); SetAutoProcessTH1(fMCHist_Slices[i]); } fDataHist = new TH1D((fSettings.GetName() + "_data").c_str(), (fSettings.GetFullTitles()).c_str(), bincount, 0, bincount); fDataHist->SetDirectory(0); int counter = 0; for (int i = 0; i < nSlices; ++i) { // Set a nice title std::string costitle = fDataHist_Slices[i]->GetTitle(); costitle = costitle.substr(costitle.find("-> ") + 3, costitle.size()); std::string found = costitle.substr(0, costitle.find(" < ")); std::string comp = costitle.substr(costitle.find(found) + found.size() + 3, costitle.size()); comp = comp.substr(comp.find(" < ") + 3, comp.size()); costitle = "cos#theta_{#mu}=" + found + "-" + comp; for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) { fDataHist->SetBinContent(counter + 1, fDataHist_Slices[i]->GetBinContent(j + 1)); fDataHist->SetBinError(counter + 1, fDataHist_Slices[i]->GetBinError(j + 1)); // Set a nice axis if (j == 0) fDataHist->GetXaxis()->SetBinLabel( counter + 1, Form("%s, p_{#mu}=%.1f-%.1f", costitle.c_str(), fDataHist_Slices[i]->GetBinLowEdge(j + 1), fDataHist_Slices[i]->GetBinLowEdge(j + 2))); else fDataHist->GetXaxis()->SetBinLabel( counter + 1, Form("p_{#mu}=%.1f-%.1f", fDataHist_Slices[i]->GetBinLowEdge(j + 1), fDataHist_Slices[i]->GetBinLowEdge(j + 2))); counter++; } } // The correlation matrix // Major in angular bins, minor in momentum bins: runs theta1, pmu1, pmu2, // theta2, pmu1, pmu2, theta3, pmu1, pmu2 etc The correlation matrix TH2D *tempcov = NULL; if (UnfoldWithGENIE) tempcov = (TH2D *)fInputFile->Get("covariance_matrix_genie"); else tempcov = (TH2D *)fInputFile->Get("covariance_matrix_neut"); fFullCovar = new TMatrixDSym(bincount); for (int i = 0; i < fDataHist->GetNbinsX(); i++) { for (int j = 0; j < fDataHist->GetNbinsX(); j++) { (*fFullCovar)(i, j) = tempcov->GetBinContent(i + 1, j + 1); } } covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); (*fFullCovar) *= 1E38 * 1E38; (*covar) *= 1E-38 * 1E-38; fInputFile->Close(); };