diff --git a/src/MINERvA/MINERvA_CC1pi0_XSec_1D_nu.cxx b/src/MINERvA/MINERvA_CC1pi0_XSec_1D_nu.cxx index d3a797b..7bd1946 100644 --- a/src/MINERvA/MINERvA_CC1pi0_XSec_1D_nu.cxx +++ b/src/MINERvA/MINERvA_CC1pi0_XSec_1D_nu.cxx @@ -1,320 +1,320 @@ // 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_CC1pi0_XSec_1D_nu.h" // Implementation of 2017 MINERvA numu CC1pi0 // arxiv:1708.03723v1 hep-ex // c.wret14@imperial.ac.uk //******************************************************************** void MINERvA_CC1pi0_XSec_1D_nu::SetupDataSettings(){ //******************************************************************** // Set Distribution // See header file for enum and some descriptions std::string name = fSettings.GetS("name"); if (!name.compare("MINERvA_CC1pi0_XSec_1DTpi_nu")) fDist = kTpi; else if (!name.compare("MINERvA_CC1pi0_XSec_1Dth_nu")) fDist= kth; else if (!name.compare("MINERvA_CC1pi0_XSec_1Dpmu_nu")) fDist= kpmu; else if (!name.compare("MINERvA_CC1pi0_XSec_1Dthmu_nu")) fDist= kthmu; else if (!name.compare("MINERvA_CC1pi0_XSec_1DQ2_nu")) fDist= kQ2; else if (!name.compare("MINERvA_CC1pi0_XSec_1DEnu_nu")) fDist= kEnu; else if (!name.compare("MINERvA_CC1pi0_XSec_1DWexp_nu")) fDist= kWexp; else if (!name.compare("MINERvA_CC1pi0_XSec_1DPPi0Mass_nu")) fDist= kPPi0Mass; else if (!name.compare("MINERvA_CC1pi0_XSec_1DPPi0MassDelta_nu")) fDist= kPPi0MassDelta; else if (!name.compare("MINERvA_CC1pi0_XSec_1DCosAdler_nu")) fDist= kCosAdler; else if (!name.compare("MINERvA_CC1pi0_XSec_1DPhiAdler_nu")) fDist= kPhiAdler; // Define what files to use from the dist std::string datafile = ""; std::string corrfile = ""; std::string titles = ""; std::string distdescript = ""; // Set the default to essentially not be a cut on proton kinetic energy // The Adler angles and reconstructed p,pi0 invariant mass have cuts on these ProtonCut = 100; // W exp is 1.8 GeV or lower (dealt with below) WexpCut = 1.8; // Load up the data switch (fDist) { case (kTpi): datafile = "data/XSec_Table_pi0_KE_xsec.csv"; corrfile = "corr/Correlation_Table_pi0_KE_xsec.csv"; titles = "CC1#pi^{0};T_{#pi} (GeV);d#sigma/dT_{#pi} (cm^{2}/nucleon/GeV)"; break; case (kth): datafile = "data/XSec_Table_pi0_theta_xsec.csv"; corrfile = "corr/Correlation_Table_pi0_theta_xsec.csv"; titles = "CC1#pi^{0};#theta_{#pi} (degrees); d#sigma/d#theta_{#pi} (cm^{2}/nucleon/degree)"; break; case (kpmu): datafile = "data/XSec_Table_muon_P_xsec.csv"; corrfile = "corr/Correlation_Table_muon_P_xsec.csv"; titles = "CC1#pi^{0};p_{#mu} (GeV);d#sigma/dp_{#mu} (cm^{2}/nucleon/GeV)"; break; case (kthmu): datafile = "data/XSec_Table_muon_theta_xsec.csv"; corrfile = "corr/Correlation_Table_muon_theta_xsec.csv"; titles = "CC1#pi^{0};#theta_{#mu} (degrees);d#sigma/d#theta_{#mu} (cm^{2}/nucleon/degree)"; break; case (kQ2): datafile = "data/XSec_Table_QSq_xsec.csv"; corrfile = "corr/Correlation_Table_QSq_xsec.csv"; titles = "CC1#pi^{0};Q^{2} (GeV^{2});d#sigma/dQ^{2} (cm^{2}/nucleon/GeV^{2})"; break; case (kEnu): datafile = "data/XSec_Table_Enu_xsec.csv"; corrfile = "corr/Correlation_Table_Enu_xsec.csv"; titles = "CC1#pi^{0};E_{#nu} (GeV);#sigma(E_#nu) (cm^{2}/nucleon)"; break; case (kWexp): datafile = "data/XSec_Table_W_xsec.csv"; corrfile = "corr/Correlation_Table_W_xsec.csv"; titles = "CC1#pi^{0};W_{exp} (GeV);d#sigma/dW_{exp} (cm^{2}/nucleon/GeV)"; break; case (kPPi0Mass): datafile = "data/XSec_Table_deltaInvMass_xsec.csv"; corrfile = "corr/Correlation_Table_deltaInvMass_xsec.csv"; titles = "CC1#pi^{0}; M_{p#pi^{0}} (GeV); d#sigma/dM_{p#pi^{0}} (cm^{2}/nucleon/GeV)"; break; case (kPPi0MassDelta): datafile = "data/XSec_Table_deltaInvMass_xsec_DeltaRich.csv"; corrfile = "corr/Correlation_Table_deltaInvMass_xsec_DeltaRich.csv"; titles = "CC1#pi^{0}; M_{p#pi^{0}} W_{exp} < 1.4 (GeV); d#sigma/dM_{p#pi^{0}} (cm^{2}/nucleon/GeV)"; break; case (kCosAdler): datafile = "data/XSec_Table_Delta_pi_theta_xsec.csv"; corrfile = "corr/Correlation_Table_Delta_pi_theta_xsec.csv"; titles = "CC1#pi^{0}; cos#theta_{Adler}; d#sigma/dcos#theta_{Adler} (cm^{2}/nucleon/0.1)"; break; case (kPhiAdler): datafile = "data/XSec_Table_Delta_pi_phi_xsec.csv"; corrfile = "corr/Correlation_Table_Delta_pi_phi_xsec.csv"; titles = "CC1#pi^{0}; #phi_{Adler} (degrees); d#sigma/d#phi_{Adler} (cm^{2}/nucleon/degree)"; break; default: THROW("Unknown Analysis Distribution : " << fDist); } // Set the Wexp and proton kinetic energy cuts depending on sample // for Ppi0Mass distributions and Adler angles we require a proton of at least 100 MeV kinetic energy if (fDist >= kPPi0Mass) { ProtonCut = 0.1; // 0.1 GeV proton kinetic energy cut // Some distributions have a Wexp cut at 1.4 GeV which attempts to isolate delta production if (fDist >= kPPi0MassDelta) { WexpCut = 1.4; } } // Only have xsec covariance (not shape only) // Now setup each data distribution and description. std::string descrip = distdescript + \ "Target: CH \n" \ "Flux: MINERvA Forward Horn Current numu ONLY \n" \ - "Signal: Any event with 1 muon, and 1pi0 in FS, no mesons, any nucleon(s). W < 1.8" \ + "Signal: Any event with 1 muon with #theta_{#mu,#nu}<25#degree, and 1pi0 in FS, no mesons, any nucleon(s). W < 1.8" \ "Alt Signal: Add in requirement of 1 proton with 100 MeV and sometimes W < 1.4"; fSettings.SetDescription(descrip); // Specify the data fSettings.SetDataInput( GeneralUtils::GetTopLevelDir()+"/data/MINERvA/CC1pi0/2017_nu/" + datafile); // And the correlations fSettings.SetCovarInput(GeneralUtils::GetTopLevelDir()+"/data/MINERvA/CC1pi0/2017_nu/" + corrfile); fSettings.SetTitle( GeneralUtils::ParseToStr(titles,";")[0] ); fSettings.SetXTitle( GeneralUtils::ParseToStr(titles,";")[1] ); fSettings.SetYTitle( GeneralUtils::ParseToStr(titles,";")[2] ); return; } //******************************************************************** MINERvA_CC1pi0_XSec_1D_nu::MINERvA_CC1pi0_XSec_1D_nu(nuiskey samplekey) { //******************************************************************** // Define Sample Settings common to all data distributions fSettings = LoadSampleSettings(samplekey); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL"); // From 1.5 to 20 GeV Enu fSettings.SetEnuRange(1.5, 20.0); fSettings.DefineAllowedTargets("C,H"); fSettings.DefineAllowedSpecies("numu"); SetupDataSettings(); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // If Enu setup scale factor for Enu Unfolded, otherwise use differential if (fDist == kEnu) fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / double(fNEvents); else fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / double(fNEvents) / TotalIntegratedFlux("width"); // Plot Setup ------------------------------------------------------- SetDataFromTextFile( fSettings.GetDataInput() ); // The errors come as a percent of the cross-section for this measurement so need rescaling // Have set the fDataHist above so just loop over and set new errors for (int i = 0; i < fDataHist->GetNbinsX(); ++i) { fDataHist->SetBinError(i+1, (fDataHist->GetBinError(i+1)/100.0)*fDataHist->GetBinContent(i+1)); } // And finally all the numbers are in units of 1E-40 so scale the histograms as such fDataHist->Scale(1.0E-40); // This measurement gives us a correlation matrix, so should set it up as such SetCorrelationFromTextFile( fSettings.GetCovarInput() ); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; //******************************************************************** void MINERvA_CC1pi0_XSec_1D_nu::FillEventVariables(FitEvent *event) { //******************************************************************** fXVar = -999.9; // Need a neutral pion and a muon if (event->NumFSParticle(111) == 0 || event->NumFSParticle(13) == 0) { return; } // Get the TLorentzVectors from the event TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Ppi0 = event->GetHMFSParticle(111)->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; // Pion kinetic energy double Tpi = (Ppi0.E() - Ppi0.Mag())/1.E3; // Pion-neutrino angle double th = (180./M_PI)*FitUtils::th(Pnu, Ppi0); // Muon momentum double pmu = Pmu.Vect().Mag()/1.E3; // Muon-neutrino angle double thmu = (180.0/M_PI)*FitUtils::th(Pnu, Pmu); // True Q2 double Q2 = fabs((Pmu - Pnu).Mag2()) / 1.E6; // True Enu double Enu = Pnu.E() / 1.E3; // Wexp (derived from "kinematic quantities" but uses EnuTrue) double Wexp = FitUtils::Wrec(Pnu, Pmu)/1.E3; // Wexp cut of 1.8 GeV in signal definition // N.B. the Adler angles and PPi0 mass requires this to be 1400 if (Wexp > WexpCut) return; + // There's a theta mu cut of 25 degrees + if (thmu > 25) return; + // Some distributions require the final state proton: check that it exists if (fDist >= kPPi0Mass && event->NumFSParticle(2212) == 0) return; // Fill the variables depending on the enums switch (fDist) { case kTpi: fXVar = Tpi; break; case kth: fXVar = th; break; case kpmu: // Pmu has a theta_mu < 25degree cut - if (thmu > 25) return; - else fXVar = pmu; + fXVar = pmu; break; case kthmu: // thmu has a theta_mu < 25degree cut - if (thmu > 25) return; - else fXVar = pmu; fXVar = thmu; break; case kQ2: fXVar = Q2; break; case kEnu: fXVar = Enu; break; case kWexp: fXVar = Wexp; break; // p, pi0 invariant mass with Wexp < 1.8 or Wexp < 1.4: already checked these above case kPPi0Mass: case kPPi0MassDelta: { // Get the proton TLorentzVector Pprot = event->GetHMFSParticle(2212)->fP; double Ppi0Mass = (Ppi0+Pprot).Mag()/1.E3; fXVar = Ppi0Mass; break; } // Cos theta Adler angle case kCosAdler: { TLorentzVector Pprot = event->GetHMFSParticle(2212)->fP; double CosThAdler = FitUtils::CosThAdler(Pnu, Pmu, Ppi0, Pprot); fXVar = CosThAdler; break; } // Phi Adler angle case kPhiAdler: { TLorentzVector Pprot = event->GetHMFSParticle(2212)->fP; double PhiAdler = FitUtils::PhiAdler(Pnu, Pmu, Ppi0, Pprot); fXVar = PhiAdler; break; } default: THROW("DIST NOT FOUND : " << fDist); break; } }; //******************************************************************** bool MINERvA_CC1pi0_XSec_1D_nu::isSignal(FitEvent *event) { //******************************************************************** // Some of the distributions require a proton with at least 100 MeV KE if (fDist >= kPPi0Mass) { // First of needs a proton in the final state if (event->NumFSParticle(2212) == 0) return false; // Needs to pass CC1pi0 signal definition bool pass_cc1pi0 = SignalDef::isCC1pi0_MINERvA_nu(event, EnuMin, EnuMax); if (!pass_cc1pi0) return false; // And the proton needs at least 100 MeV kinetic energy TLorentzVector Pprot = event->GetHMFSParticle(2212)->fP; double ke = (Pprot.E() - Pprot.Mag())/1.E3; if (pass_cc1pi0 && ke > ProtonCut) { return true; } else { return false; } // The other distributions ahve a more generic "1mu, 1pi0, no mesons, any nucleon(s)" // The W cut is instead made in FillEventVariables } else { return SignalDef::isCC1pi0_MINERvA_nu(event, EnuMin, EnuMax); } }