diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.cxx index e7fe2e8..418b4f3 100644 --- a/src/T2K/T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.cxx +++ b/src/T2K/T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.cxx @@ -1,111 +1,111 @@ #include #include "T2K_SignalDef.h" #include "T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.h" // The constructor T2K_CC1pip_CH_XSec_1DAdlerPhi_nu::T2K_CC1pip_CH_XSec_1DAdlerPhi_nu( nuiskey samplekey) { // Sample overview --------------------------------------------------- std::string descrip = "T2K_CC1pip_CH_XSec_nu sample. \n" "Target: CH \n" "Flux: T2K FHC numu \n" "Signal: CC1pi+, p_mu > 200 MeV, p_pi > 200 MeV\n" ", costheta_mu > 0.2, costheta_pi > 0.2\n" "https://arxiv.org/abs/1909.03936"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetTitle("T2K_CC1pip_CH_XSec_1DAdlerPhi_nu"); fSettings.SetDescription(descrip); - fSettings.SetXTitle("#phi_{Adler} (radians)"); - fSettings.SetYTitle("d#sigma/d#phi_{Adler} (cm^{2}/rad/nucleon)"); + fSettings.SetXTitle("#phi_{Adler} (rad.)"); + fSettings.SetYTitle("d#sigma/d#phi_{Adler} (cm^{2}/rad./nucleon)"); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG"); fSettings.SetEnuRange(0.0, 100.0); fSettings.DefineAllowedTargets("C,H"); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) / double(fNEvents) / TotalIntegratedFlux("width"); // Plot Setup ------------------------------------------------------- SetDataFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/phi_adler.root", "Phi_Adler"); SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/phi_adler.root", "Phi_AdlerCov"); SetShapeCovar(); fDataHist->Scale(1E-38); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; void T2K_CC1pip_CH_XSec_1DAdlerPhi_nu::FillEventVariables(FitEvent *event) { if (event->NumFSParticle(13) == 0 || event->NumFSParticle(211) == 0) return; // Reconstruct the neutrino TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; double rEnu = FitUtils::EnuCC1piprec_T2K_eMB(Pnu, Pmu, Ppip) * 1000.; // Now make the reconstructed neutrino // Has the same direction as the predicted neutrino TLorentzVector PnuReco(Pnu.Vect().Unit().X() * rEnu, Pnu.Vect().Unit().Y() * rEnu, Pnu.Vect().Unit().Z() * rEnu, rEnu); // Reconstruct the initial state assuming still nucleon TLorentzVector Pinit(0, 0, 0, PhysConst::mass_proton * 1000.); // Pretty much a copy of FitUtils::PhiAdler but using the reconstructed // neutrino instead of true neutrino{ // Get the "resonance" lorentz vector (pion proton system) reconstructed from // the variables TLorentzVector Pres = PnuReco + Pinit - Pmu; // Boost the particles into the resonance rest frame so we can define the // x,y,z axis PnuReco.Boost(-Pres.BoostVector()); Pmu.Boost(-Pres.BoostVector()); Ppip.Boost(-Pres.BoostVector()); // The z-axis is defined as the axis of three-momentum transfer, \vec{k} // Also unit normalise the axis TVector3 zAxis = (PnuReco.Vect() - Pmu.Vect()) * (1.0 / ((PnuReco.Vect() - Pmu.Vect()).Mag())); // The y-axis is then defined perpendicular to z and muon momentum in the // resonance frame TVector3 yAxis = PnuReco.Vect().Cross(Pmu.Vect()); yAxis *= 1.0 / double(yAxis.Mag()); // And the x-axis is then simply perpendicular to z and x TVector3 xAxis = yAxis.Cross(zAxis); xAxis *= 1.0 / double(xAxis.Mag()); double x = Ppip.Vect().Dot(xAxis); double y = Ppip.Vect().Dot(yAxis); double newphi = atan2(y, x); fXVar = newphi; return; }; //******************************************************************** bool T2K_CC1pip_CH_XSec_1DAdlerPhi_nu::isSignal(FitEvent *event) { //******************************************************************** return SignalDef::isCC1pip_T2K_arxiv1909_03936( event, EnuMin, EnuMax, SignalDef::kMuonHighEff | SignalDef::kPionVFwd | SignalDef::kPionHighMom); } diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1DCosThAdler_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1DCosThAdler_nu.cxx index 7f79aba..376b9c3 100644 --- a/src/T2K/T2K_CC1pip_CH_XSec_1DCosThAdler_nu.cxx +++ b/src/T2K/T2K_CC1pip_CH_XSec_1DCosThAdler_nu.cxx @@ -1,102 +1,102 @@ #include #include "T2K_SignalDef.h" #include "T2K_CC1pip_CH_XSec_1DCosThAdler_nu.h" // The constructor T2K_CC1pip_CH_XSec_1DCosThAdler_nu::T2K_CC1pip_CH_XSec_1DCosThAdler_nu( nuiskey samplekey) { // Sample overview --------------------------------------------------- std::string descrip = "T2K_CC1pip_CH_XSec_nu sample. \n" "Target: CH \n" "Flux: T2K FHC numu \n" "Signal: CC1pi+, p_mu > 200 MeV, p_pi > 200 MeV\n" ", costheta_mu > 0.2, costheta_pi > 0.2\n" "https://arxiv.org/abs/1909.03936"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetTitle("T2K_CC1pip_CH_XSec_1DCosThAdler_nu"); fSettings.SetDescription(descrip); - fSettings.SetXTitle("cos#theta_{Adler} (radians)"); + fSettings.SetXTitle("cos#theta_{Adler} (rad.)"); fSettings.SetYTitle("d#sigma/dcos#theta_{Adler} (cm^{2}/1/nucleon)"); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG"); fSettings.SetEnuRange(0.0, 100.0); fSettings.DefineAllowedTargets("C,H"); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) / double(fNEvents) / TotalIntegratedFlux("width"); // Plot Setup ------------------------------------------------------- SetDataFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/theta_adler.root", "Theta_Adler"); SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/theta_adler.root", "Theta_AdlerCov"); SetShapeCovar(); fDataHist->Scale(1E-38); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; void T2K_CC1pip_CH_XSec_1DCosThAdler_nu::FillEventVariables(FitEvent *event) { if (event->NumFSParticle(13) == 0 || event->NumFSParticle(211) == 0) return; // Essentially the same as FitUtils::CosThAdler but uses the reconsturcted // neutrino instead of the true neutrino Reconstruct the neutrino TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; double rEnu = FitUtils::EnuCC1piprec_T2K_eMB(Pnu, Pmu, Ppip) * 1000.; // Now make the reconstructed neutrino // Has the same direction as the predicted neutrino TLorentzVector PnuReco(Pnu.Vect().Unit().X() * rEnu, Pnu.Vect().Unit().Y() * rEnu, Pnu.Vect().Unit().Z() * rEnu, rEnu); // Reconstruct the initial state assuming still nucleon TLorentzVector Pinit(0, 0, 0, PhysConst::mass_proton * 1000.); // Pretty much a copy of FitUtils::PhiAdler but using the reconstructed // neutrino instead of true neutrino{ // Get the "resonance" lorentz vector (pion proton system) reconstructed from // the variables TLorentzVector Pres = PnuReco + Pinit - Pmu; // Boost the particles into the resonance rest frame so we can define the // x,y,z axis PnuReco.Boost(-Pres.BoostVector()); Pmu.Boost(-Pres.BoostVector()); Ppip.Boost(-Pres.BoostVector()); // The z-axis is defined as the axis of three-momentum transfer, \vec{k} // Also unit normalise the axis TVector3 zAxis = (PnuReco.Vect() - Pmu.Vect()) * (1.0 / ((PnuReco.Vect() - Pmu.Vect()).Mag())); // Then the angle between the pion in the "resonance" rest-frame and the // z-axis is the theta Adler angle double costhAdler = cos(Ppip.Vect().Angle(zAxis)); fXVar = costhAdler; return; }; //******************************************************************** bool T2K_CC1pip_CH_XSec_1DCosThAdler_nu::isSignal(FitEvent *event) { //******************************************************************** return SignalDef::isCC1pip_T2K_arxiv1909_03936( event, EnuMin, EnuMax, SignalDef::kMuonHighEff | SignalDef::kPionVFwd | SignalDef::kPionHighMom); } diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1DQ2_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1DQ2_nu.cxx index 0ef525e..879a4b4 100644 --- a/src/T2K/T2K_CC1pip_CH_XSec_1DQ2_nu.cxx +++ b/src/T2K/T2K_CC1pip_CH_XSec_1DQ2_nu.cxx @@ -1,102 +1,102 @@ #include #include "T2K_SignalDef.h" #include "T2K_CC1pip_CH_XSec_1DQ2_nu.h" // The constructor T2K_CC1pip_CH_XSec_1DQ2_nu::T2K_CC1pip_CH_XSec_1DQ2_nu(nuiskey samplekey) { // Sample overview --------------------------------------------------- std::string descrip = "T2K_CC1pip_CH_XSec_nu sample. \n" "Target: CH \n" "Flux: T2K FHC numu \n" "Signal: CC1pi+, p_mu > 200 MeV, p_pi > 200 MeV\n" ", costheta_mu > 0.2, costheta_pi > 0.2\n" "https://arxiv.org/abs/1909.03936"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetTitle("T2K_CC1pip_CH_XSec_1DQ2_nu"); fSettings.SetDescription(descrip); - fSettings.SetXTitle("Q^{2} (GeV/c)^{2}"); + fSettings.SetXTitle("Q^{2} (GeV^{2})"); fSettings.SetYTitle("d#sigma/dQ^{2} (cm^{2}/GeV^{2}/nucleon)"); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG"); fSettings.SetEnuRange(0.0, 100.0); fSettings.DefineAllowedTargets("C,H"); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) / double(fNEvents) / TotalIntegratedFlux("width"); // Plot Setup ------------------------------------------------------- SetDataFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Q2.root", "Q2"); SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Q2.root", "Q2Cov"); SetShapeCovar(); fDataHist->Scale(1E-38); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; void T2K_CC1pip_CH_XSec_1DQ2_nu::FillEventVariables(FitEvent *event) { if (event->NumFSParticle(13) == 0 || event->NumFSParticle(211) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; double q2 = -999; // switch(fT2KSampleType) { // First int refers to how we reconstruct Enu // 0 uses true neutrino energy (not used here but common for other analyses // when they unfold to true Enu from reconstructed Enu) 1 uses "extended // MiniBooNE" method 2 uses "MiniBooNE reconstructed" method 3 uses Delta // resonance mass for reconstruction // // The last bool refers to if pion directional information was used // // Use MiniBooNE reconstructed Enu; uses Michel tag so no pion direction // information // case kMB: // q2 = FitUtils::Q2CC1piprec(Pnu, Pmu, Ppip, 2, false); // break; // Use Extended MiniBooNE reconstructed Enu // Needs pion information to reconstruct so bool is true (did not include // Michel e tag) // case keMB: q2 = FitUtils::Q2CC1piprec(Pnu, Pmu, Ppip, 1, true); // break; // Use Delta resonance reconstructed Enu // Uses Michel electron so don't have pion directional information // case kDelta: // q2 = FitUtils::Q2CC1piprec(Pnu, Pmu, Ppip, 3, false); // break; //} fXVar = q2; return; }; //******************************************************************** bool T2K_CC1pip_CH_XSec_1DQ2_nu::isSignal(FitEvent *event) { //******************************************************************** return SignalDef::isCC1pip_T2K_arxiv1909_03936( event, EnuMin, EnuMax, SignalDef::kMuonHighEff | SignalDef::kPionVFwd | SignalDef::kPionHighMom); } diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dppi_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1Dppi_nu.cxx index 041ea16..718a11f 100644 --- a/src/T2K/T2K_CC1pip_CH_XSec_1Dppi_nu.cxx +++ b/src/T2K/T2K_CC1pip_CH_XSec_1Dppi_nu.cxx @@ -1,71 +1,71 @@ #include #include "T2K_CC1pip_CH_XSec_1Dppi_nu.h" #include "T2K_SignalDef.h" //******************************************************************** T2K_CC1pip_CH_XSec_1Dppi_nu::T2K_CC1pip_CH_XSec_1Dppi_nu(nuiskey samplekey) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "T2K_CC1pip_CH_XSec_nu sample. \n" "Target: CH \n" "Flux: T2K FHC numu \n" "Signal: CC1pi+, p_mu > 200 MeV, p_pi > 200 MeV\n" ", costheta_mu > 0.2, costheta_pi > 0.2\n" "https://arxiv.org/abs/1909.03936"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetTitle("T2K_CC1pip_CH_XSec_1Dppi_nu"); fSettings.SetDescription(descrip); - fSettings.SetXTitle("p_{#pi} (GeV/c)"); - fSettings.SetYTitle("d#sigma/dp_{#pi} (cm^{2}/(GeV/c)/nucleon)"); + fSettings.SetXTitle("p_{#pi} (GeV)"); + fSettings.SetYTitle("d#sigma/dp_{#pi} (cm^{2}/GeV/nucleon)"); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG"); fSettings.SetEnuRange(0.0, 100.0); fSettings.DefineAllowedTargets("C,H"); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) / double(fNEvents) / TotalIntegratedFlux("width"); // Plot Setup ------------------------------------------------------- SetDataFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/MomentumPion.root", "Momentum_pion"); SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/MomentumPion.root", "Momentum_pionCov"); SetShapeCovar(); fDataHist->Scale(1E-38); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; void T2K_CC1pip_CH_XSec_1Dppi_nu::FillEventVariables(FitEvent *event) { if (event->NumFSParticle(13) == 0 || event->NumFSParticle(211) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; double ppip = FitUtils::p(Ppip); fXVar = ppip; return; }; //******************************************************************** bool T2K_CC1pip_CH_XSec_1Dppi_nu::isSignal(FitEvent *event) { //******************************************************************** return SignalDef::isCC1pip_T2K_arxiv1909_03936( event, EnuMin, EnuMax, SignalDef::kMuonHighEff | SignalDef::kPionVFwd); } diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx index fd4826a..db2a4c5 100644 --- a/src/T2K/T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx +++ b/src/T2K/T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx @@ -1,72 +1,72 @@ #include #include "T2K_CC1pip_CH_XSec_1Dthmupi_nu.h" #include "T2K_SignalDef.h" // The constructor T2K_CC1pip_CH_XSec_1Dthmupi_nu::T2K_CC1pip_CH_XSec_1Dthmupi_nu( nuiskey samplekey) { // Sample overview --------------------------------------------------- std::string descrip = "T2K_CC1pip_CH_XSec_nu sample. \n" "Target: CH \n" "Flux: T2K FHC numu \n" "Signal: CC1pi+, p_mu > 200 MeV, p_pi > 200 MeV\n" ", costheta_mu > 0.2, costheta_pi > 0.2\n" "https://arxiv.org/abs/1909.03936"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetTitle("T2K_CC1pip_CH_XSec_1Dthmupi_nu"); fSettings.SetDescription(descrip); - fSettings.SetXTitle("#theta_{#mu,#pi} (radians)"); - fSettings.SetYTitle("d#sigma/d#theta_{#mu,#pi} (cm^{2}/radians/nucleon)"); + fSettings.SetXTitle("#theta_{#mu,#pi} (rad.)"); + fSettings.SetYTitle("d#sigma/d#theta_{#mu,#pi} (cm^{2}/rad./nucleon)"); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG"); fSettings.SetEnuRange(0.0, 100.0); fSettings.DefineAllowedTargets("C,H"); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) / double(fNEvents) / TotalIntegratedFlux("width"); // Plot Setup ------------------------------------------------------- SetDataFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Thetapimu.root", "Theta(pi,mu)(rads)"); SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Thetapimu.root", "Theta(pi,mu)(rads)Cov"); SetShapeCovar(); fDataHist->Scale(1E-38); FinaliseMeasurement(); }; void T2K_CC1pip_CH_XSec_1Dthmupi_nu::FillEventVariables(FitEvent *event) { if (event->NumFSParticle(13) == 0 || event->NumFSParticle(211) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; double thmupi = FitUtils::th(Pmu, Ppip); fXVar = thmupi; // std::cout << thmupi << std::endl; return; }; //******************************************************************** bool T2K_CC1pip_CH_XSec_1Dthmupi_nu::isSignal(FitEvent *event) { //******************************************************************** return SignalDef::isCC1pip_T2K_arxiv1909_03936( event, EnuMin, EnuMax, SignalDef::kMuonHighEff | SignalDef::kPionVFwd | SignalDef::kPionHighMom); } diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx index abc4ea1..9244b8c 100644 --- a/src/T2K/T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx +++ b/src/T2K/T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx @@ -1,67 +1,67 @@ #include #include "T2K_CC1pip_CH_XSec_1Dthpi_nu.h" #include "T2K_SignalDef.h" // The constructor T2K_CC1pip_CH_XSec_1Dthpi_nu::T2K_CC1pip_CH_XSec_1Dthpi_nu(nuiskey samplekey) { // Sample overview --------------------------------------------------- std::string descrip = "T2K_CC1pip_CH_XSec_nu sample. \n" "Target: CH \n" "Flux: T2K FHC numu \n" "Signal: CC1pi+, p_mu > 200 MeV, p_pi > 200 MeV\n" ", costheta_mu > 0.2, costheta_pi > 0\n" "https://arxiv.org/abs/1909.03936"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetTitle("T2K_CC1pip_CH_XSec_1Dthpi_nu"); fSettings.SetDescription(descrip); - fSettings.SetXTitle("#theta_{#pi} (radians)"); - fSettings.SetYTitle("d#sigma/d#theta_{#pi} (cm^{2}/radians/nucleon)"); + fSettings.SetXTitle("#theta_{#pi} (rad.)"); + fSettings.SetYTitle("d#sigma/d#theta_{#pi} (cm^{2}/rad./nucleon)"); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG"); fSettings.SetEnuRange(0.0, 100.0); fSettings.DefineAllowedTargets("C,H"); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) / double(fNEvents) / TotalIntegratedFlux("width"); // Plot Setup ------------------------------------------------------- SetDataFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Thetapion.root", "Theta_pion"); SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Thetapion.root", "Theta_pionCov"); SetShapeCovar(); fDataHist->Scale(1E-38); FinaliseMeasurement(); }; void T2K_CC1pip_CH_XSec_1Dthpi_nu::FillEventVariables(FitEvent *event) { if (event->NumFSParticle(211) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; double thpi = FitUtils::th(Pnu, Ppip); fXVar = thpi; }; //******************************************************************** bool T2K_CC1pip_CH_XSec_1Dthpi_nu::isSignal(FitEvent *event) { //******************************************************************** return SignalDef::isCC1pip_T2K_arxiv1909_03936( event, EnuMin, EnuMax, SignalDef::kMuonHighEff | SignalDef::kPionFwd | SignalDef::kPionHighMom); } diff --git a/src/T2K/T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.cxx index cec28d0..1d41d6d 100644 --- a/src/T2K/T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.cxx +++ b/src/T2K/T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.cxx @@ -1,211 +1,208 @@ #include "T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.h" #include "T2K_SignalDef.h" //******************************************************************** T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::T2K_CC1pip_CH_XSec_2Dpmucosmu_nu( nuiskey samplekey) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "T2K_CC1pip_CH_XSec_nu sample. \n" "Target: CH \n" "Flux: T2K FHC numu \n" "Signal: CC1pi+, costheta_mu > 0" "https://arxiv.org/abs/1909.03936"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); - fSettings.SetXTitle(" "); - fSettings.SetYTitle( - "d^{2}#sigma/dp_{#mu}dcos#theta_{#mu} (cm^{2}/(GeV/c)/nucleon)"); + fSettings.SetXTitle("p_{#mu}-cos#theta_{#mu}"); + fSettings.SetYTitle("d^{2}#sigma/dp_{#mu}dcos#theta_{#mu} (cm^{2}/GeV/nucleon)"); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG"); fSettings.SetEnuRange(0.0, 100.0); fSettings.DefineAllowedTargets("C,H"); // CCQELike plot information fSettings.SetTitle("T2K_CC1pip_CH_XSec_2Dpmucosmu_nu"); fSettings.SetDataInput(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/PmuThetamu.root"); fSettings.SetCovarInput(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/PmuThetamu.root"); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) / double(fNEvents) / TotalIntegratedFlux("width"); // Plot Setup ------------------------------------------------------- SetHistograms(); covar = StatUtils::GetInvert(fFullCovar, true); fDecomp = StatUtils::GetDecomp(fFullCovar); SetShapeCovar(); // Final setup --------------------------------------------------- FinaliseMeasurement(); + fSaveFine = false; }; void T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::SetHistograms() { TFile *data = new TFile(fSettings.GetDataInput().c_str(), "open"); std::string dataname = "p_mu_theta_mu"; // Number of slices we have const int nslices = 4; int nbins = 0; for (int i = 0; i < nslices; ++i) { TH1D *slice = (TH1D *)data->Get(Form("%s_%i", dataname.c_str(), i)); - slice = (TH1D *)slice->Clone((fName + Form("_data_slice%i", i)).c_str()); + slice = (TH1D *)slice->Clone((fName + Form("_data_Slice%i", i)).c_str()); slice->Scale(1E-38); - slice->GetXaxis()->SetTitle(fSettings.GetS("xtitle").c_str()); - slice->GetYaxis()->SetTitle(fSettings.GetS("ytitle").c_str()); + slice->GetXaxis()->SetTitle("p_{#mu}"); + slice->GetYaxis()->SetTitle(fSettings.GetYTitle().c_str()); fDataHist_Slices.push_back(slice); fMCHist_Slices.push_back( - (TH1D *)slice->Clone((fName + Form("_mc_slice%i", i)).c_str())); + (TH1D *)slice->Clone((fName + Form("_MC_Slice%i", i)).c_str())); SetAutoProcessTH1(fDataHist_Slices[i], kCMD_Write); SetAutoProcessTH1(fMCHist_Slices[i]); fMCHist_Slices[i]->Reset(); fMCHist_Slices[i]->SetLineColor(kRed); //nbins += slice->GetXaxis()->GetNbins(); // Skip the highest momentum bin because it's rubbish nbins += slice->GetXaxis()->GetNbins() - 1; } fDataHist = new TH1D(dataname.c_str(), dataname.c_str(), nbins, 0, nbins); fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_data").c_str()); int bincount = 1; for (int i = 0; i < nslices; ++i) { // Skip the highest momentum bin because it's rubbish for (int j = 0; j < fDataHist_Slices[i]->GetXaxis()->GetNbins() - 1; ++j) { //for (int j = 0; j < fDataHist_Slices[i]->GetXaxis()->GetNbins(); ++j) { fDataHist->SetBinContent(bincount, fDataHist_Slices[i]->GetBinContent(j + 1)); fDataHist->SetBinError(bincount, fDataHist_Slices[i]->GetBinError(j + 1)); TString title; if (j == 0) { title = "cos#theta_{#mu}="; if (i == 0) { title += "0-0.8, "; } else if (i == 1) { title += "0.8-0.85, "; } else if (i == 2) { title += "0.85-0.90, "; } else if (i == 3) { title += "0.90-1.00, "; } } title += Form("p_{#mu}=%.2f-%.2f", fDataHist_Slices[i]->GetBinLowEdge(j + 1), fDataHist_Slices[i]->GetBinLowEdge(j + 2)); fDataHist->GetXaxis()->SetBinLabel(bincount, title); bincount++; } } fDataHist->GetXaxis()->SetTitle(fSettings.GetS("xtitle").c_str()); fDataHist->GetYaxis()->SetTitle(fSettings.GetS("ytitle").c_str()); // Get the covariance TMatrixDSym *temp = StatUtils::GetCovarFromRootFile(fSettings.GetCovarInput(), "Covariance_pmu_thetamu"); int ncovbins = temp->GetNrows(); // Skip the highest momentum bin because it's rubbish fFullCovar = new TMatrixDSym(ncovbins - 4); //fFullCovar = new TMatrixDSym(ncovbins); if (fFullCovar->GetNrows() != fDataHist->GetXaxis()->GetNbins()*fDataHist->GetYaxis()->GetNbins()) { NUIS_ERR(FTL, "Number of bins in covariance matrix does not match data"); } // Number of costhetamu slices is nslices // Number of pmu slices is int count1 = 0; // Skip the highest momentum bin because it's rubbish for (int i = 0; i < ncovbins - 4; ++i) { //for (int i = 0; i < ncovbins; ++i) { int count2 = 0; // Skip the highest momentum bin because it's rubbish for (int j = 0; j < ncovbins - 4; ++j) { //for (int j = 0; j < ncovbins; ++j) { // 1E79 matched to diagonal error (*fFullCovar)(count1, count2) = (*temp)(i, j); count2++; } count1++; } - // Now reorganise the rows - delete temp; }; void T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::FillEventVariables(FitEvent *event) { if (event->NumFSParticle(13) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; double pmu = FitUtils::p(Pmu); double costhmu = cos(FitUtils::th(Pnu, Pmu)); fXVar = pmu; fYVar = costhmu; }; void T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::FillHistograms() { Measurement1D::FillHistograms(); if (Signal) { FillMCSlice(fXVar, fYVar, Weight); } }; void T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::ConvertEventRates() { const int nslices = 4; for (int i = 0; i < nslices; i++) { fMCHist_Slices[i]->GetSumw2(); } // Do standard conversion. Measurement1D::ConvertEventRates(); // First scale MC slices also by their width in Y and Z fMCHist_Slices[0]->Scale(1.0 / 0.80); fMCHist_Slices[1]->Scale(1.0 / 0.05); fMCHist_Slices[2]->Scale(1.0 / 0.05); fMCHist_Slices[3]->Scale(1.0 / 0.10); // Now Convert into 1D list fMCHist->Reset(); int bincount = 1; for (int i = 0; i < nslices; i++) { for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX() - 1; j++) { - //for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) { fMCHist->SetBinContent(bincount, fMCHist_Slices[i]->GetBinContent(j + 1)); bincount++; } } } void T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::FillMCSlice(double pmu, double cosmu, double weight) { // Hard code the bin edges in here if (cosmu < 0.8) { fMCHist_Slices[0]->Fill(pmu, weight); } else if (cosmu > 0.8 && cosmu < 0.85) { fMCHist_Slices[1]->Fill(pmu, weight); } else if (cosmu > 0.85 && cosmu < 0.90) { fMCHist_Slices[2]->Fill(pmu, weight); } else if (cosmu > 0.90 && cosmu < 1.00) { fMCHist_Slices[3]->Fill(pmu, weight); } }; //******************************************************************** bool T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::isSignal(FitEvent *event) { //******************************************************************** return SignalDef::isCC1pip_T2K_arxiv1909_03936(event, EnuMin, EnuMax, SignalDef::kMuonFwd); }; diff --git a/src/T2K/T2K_CC1pip_H2O_XSec_1Dppi_nu.cxx b/src/T2K/T2K_CC1pip_H2O_XSec_1Dppi_nu.cxx index 52402e1..1ecb208 100644 --- a/src/T2K/T2K_CC1pip_H2O_XSec_1Dppi_nu.cxx +++ b/src/T2K/T2K_CC1pip_H2O_XSec_1Dppi_nu.cxx @@ -1,79 +1,79 @@ #include "T2K_CC1pip_H2O_XSec_1Dppi_nu.h" // The momentum of the (positive) pion //******************************************************************** T2K_CC1pip_H2O_XSec_1Dppi_nu::T2K_CC1pip_H2O_XSec_1Dppi_nu(nuiskey samplekey) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "T2K_CC1pip_H2O_XSec_nu sample. \n" "Target: H20 \n" "Flux: T2K FHC numu \n" "Signal: CC1pi+, p_mu > 200 MeV, p_pi > 200 MeV\n" ", costheta_mu > 0.3, costheta_pi > 0.3\n" "https://journals.aps.org/prd/pdf/10.1103/PhysRevD.95.012010"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetTitle("T2K_CC1pip_H2O_XSec_1Dppi_nu"); fSettings.SetDescription(descrip); - fSettings.SetXTitle("p_{#pi^{+}} (GeV/c)"); - fSettings.SetYTitle("d#sigma/dp_{#pi^{+}} (cm^{2}/(GeV/c)/nucleon)"); + fSettings.SetXTitle("p_{#pi^{+}} (GeV)"); + fSettings.SetYTitle("d#sigma/dp_{#pi^{+}} (cm^{2}/GeV/nucleon)"); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL"); fSettings.SetEnuRange(0.0, 100.0); fSettings.DefineAllowedTargets("C,H"); fSettings.DefineAllowedSpecies("numu"); fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/H2O/" "nd280data-numu-cc1pi-xs-on-h2o-2015.root;PosPionMom/hResultTot"); fSettings.SetCovarInput( GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/H2O/" "nd280data-numu-cc1pi-xs-on-h2o-2015.root;PosPionMom/TotalCovariance"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) / double(fNEvents) / TotalIntegratedFlux("width"); // Plot Setup ------------------------------------------------------- SetDataFromRootFile(fSettings.GetDataInput()); SetCovarFromRootFile(fSettings.GetCovarInput()); ScaleCovar(1E76); SetShapeCovar(); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; //******************************************************************** // Find the momentum of the (positively charged) pion void T2K_CC1pip_H2O_XSec_1Dppi_nu::FillEventVariables(FitEvent *event) { //******************************************************************** // Need to make sure there's a muon if (event->NumFSParticle(211) == 0) return; // Get the pion TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; double p_pi = FitUtils::p(Ppip); fXVar = p_pi; return; }; //******************************************************************** // Beware: The H2O analysis has different signal definition to the CH analysis! bool T2K_CC1pip_H2O_XSec_1Dppi_nu::isSignal(FitEvent *event) { //******************************************************************** return SignalDef::isCC1pip_T2K_PRD97_012001(event, EnuMin, EnuMax); }