diff --git a/data/T2K/CC1pip/CH/MomentumPion.rootout.root b/data/T2K/CC1pip/CH/MomentumPion.rootout.root deleted file mode 100644 index 5b2c080..0000000 Binary files a/data/T2K/CC1pip/CH/MomentumPion.rootout.root and /dev/null differ diff --git a/data/T2K/CC1pip/CH/PmuThetamu.root b/data/T2K/CC1pip/CH/PmuThetamu.root index 041dd81..89585a2 100644 Binary files a/data/T2K/CC1pip/CH/PmuThetamu.root and b/data/T2K/CC1pip/CH/PmuThetamu.root differ diff --git a/data/T2K/CC1pip/CH/PmuThetamu.rootout.root b/data/T2K/CC1pip/CH/PmuThetamu.rootout.root deleted file mode 100644 index f60fbd2..0000000 Binary files a/data/T2K/CC1pip/CH/PmuThetamu.rootout.root and /dev/null differ diff --git a/data/T2K/CC1pip/CH/Q2.rootout.root b/data/T2K/CC1pip/CH/Q2.rootout.root deleted file mode 100644 index f23e9b1..0000000 Binary files a/data/T2K/CC1pip/CH/Q2.rootout.root and /dev/null differ diff --git a/data/T2K/CC1pip/CH/Thetapimu.rootout.root b/data/T2K/CC1pip/CH/Thetapimu.rootout.root deleted file mode 100644 index 154bf5a..0000000 Binary files a/data/T2K/CC1pip/CH/Thetapimu.rootout.root and /dev/null differ diff --git a/data/T2K/CC1pip/CH/Thetapion.rootout.root b/data/T2K/CC1pip/CH/Thetapion.rootout.root deleted file mode 100644 index 53e2b5a..0000000 Binary files a/data/T2K/CC1pip/CH/Thetapion.rootout.root and /dev/null differ diff --git a/data/T2K/CC1pip/CH/phi_adler.rootout.root b/data/T2K/CC1pip/CH/phi_adler.rootout.root deleted file mode 100644 index 2ca4a16..0000000 Binary files a/data/T2K/CC1pip/CH/phi_adler.rootout.root and /dev/null differ diff --git a/data/T2K/CC1pip/CH/theta_adler.rootout.root b/data/T2K/CC1pip/CH/theta_adler.rootout.root deleted file mode 100644 index 544ad36..0000000 Binary files a/data/T2K/CC1pip/CH/theta_adler.rootout.root and /dev/null differ diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.cxx index 19e493d..e7fe2e8 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.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.rootout.root", + "/data/T2K/CC1pip/CH/phi_adler.root", "Phi_Adler"); SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() + - "/data/T2K/CC1pip/CH/phi_adler.rootout.root", + "/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 8b1ccac..7f79aba 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.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.rootout.root", + "/data/T2K/CC1pip/CH/theta_adler.root", "Theta_Adler"); SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() + - "/data/T2K/CC1pip/CH/theta_adler.rootout.root", + "/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 5da96bc..0ef525e 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.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.rootout.root", + "/data/T2K/CC1pip/CH/Q2.root", "Q2"); SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() + - "/data/T2K/CC1pip/CH/Q2.rootout.root", + "/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 9c7c211..041ea16 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.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.rootout.root", + "/data/T2K/CC1pip/CH/MomentumPion.root", "Momentum_pion"); SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() + - "/data/T2K/CC1pip/CH/MomentumPion.rootout.root", + "/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 ed4002c..fd4826a 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.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.rootout.root", + "/data/T2K/CC1pip/CH/Thetapimu.root", "Theta(pi,mu)(rads)"); SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() + - "/data/T2K/CC1pip/CH/Thetapimu.rootout.root", + "/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 9028886..abc4ea1 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.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.rootout.root", + "/data/T2K/CC1pip/CH/Thetapion.root", "Theta_pion"); SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() + - "/data/T2K/CC1pip/CH/Thetapion.rootout.root", + "/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 161a179..0211301 100644 --- a/src/T2K/T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.cxx +++ b/src/T2K/T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.cxx @@ -1,215 +1,211 @@ #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.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 ------------------------------------------------------- - // SetDataValues( fSettings.GetDataInput() ); - // SetCovarMatrix( fSettings.GetCovarInput() ); SetHistograms(); - // fFullCovar = StatUtils::GetCovarFromRootFile(fSettings.GetCovarInput(), - //"Covariance_pmu_thetamu"); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); SetShapeCovar(); - /* - for (int i = 0; i < covar->GetNrows(); ++i) { - for (int j = 0; j < covar->GetNrows(); ++j) { - if (i == j) std::cout << i << " " << j << " = " << 1/sqrt((*covar)(i,j)) - << std::endl; - } - } - throw; - */ // Final setup --------------------------------------------------- FinaliseMeasurement(); }; void T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::SetHistograms() { TFile *data = new TFile(fSettings.GetDataInput().c_str(), "open"); - // std::string dataname = fSettings.Get 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->Scale(1E-38); slice->GetXaxis()->SetTitle(fSettings.GetS("xtitle").c_str()); slice->GetYaxis()->SetTitle(fSettings.GetS("ytitle").c_str()); fDataHist_Slices.push_back(slice); fMCHist_Slices.push_back( (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(); - nbins += slice->GetXaxis()->GetNbins() - 1; + nbins += slice->GetXaxis()->GetNbins(); + //nbins += slice->GetXaxis()->GetNbins() - 1; + std::cout << "slice " << i << " nbins: " << slice->GetXaxis()->GetNbins() << std::endl; + for (int j = 0; j < slice->GetXaxis()->GetNbins()+1; ++j) { + std::cout << slice->GetXaxis()->GetBinLowEdge(j+1) << "-"; + } + std::cout << std::endl; } 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) { - for (int j = 0; j < fDataHist_Slices[i]->GetXaxis()->GetNbins() - 1; ++j) { + //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(); - fFullCovar = new TMatrixDSym(ncovbins - 4); + //fFullCovar = new TMatrixDSym(ncovbins - 4); + fFullCovar = new TMatrixDSym(ncovbins); if (ncovbins != fDataHist->GetXaxis()->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; - for (int i = 0; i < ncovbins - 4; ++i) { + //for (int i = 0; i < ncovbins - 4; ++i) { + for (int i = 0; i < ncovbins; ++i) { int count2 = 0; - for (int j = 0; j < ncovbins - 4; ++j) { + //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() - 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); };