feat: introduce plotting of particle-level kinematic parameters and make plotting with KinemDistribution configurable

README.md

@@ -244,7 +244,7 @@ You can plot kinematic distribution of your sample using
 ```bash
 ./bin/KinemDistributionTutorial TutorialConfigs/FitterConfig.yaml
 ```
-Notice same config being used. In other words you can add or disable sample, modify cuts in same way as was discussed.
+Notice same config being used. In other words you can add or disable sample, modify cuts in same way as was discussed. At the bottom of the config, you can specify any individual variables you would like to plot with this executable, along with any selection cuts for each plot.
 Example of plot you can see here:
 
 <img width="350" alt="Kinematic example" src="https://github.com/user-attachments/assets/534bcb17-f26c-4fc2-a77a-5d253b0ed241">

SamplesTutorial/SampleHandlerTutorial.cpp

@@ -8,6 +8,11 @@ SampleHandlerTutorial::SampleHandlerTutorial(const std::string& config_name, Par
   KinematicParameters = &KinematicParametersTutorial;
   ReversedKinematicParameters = &ReversedKinematicParametersTutorial;
 
+  // === JM assign kinematic vector maps ===
+  KinematicVectors = &KinematicVectorsTutorial;
+  ReversedKinematicVectors = &ReversedKinematicVectorsTutorial;
+  // =======================================
+
   isATM = false;
   Initialise();
 }
@@ -103,6 +108,10 @@ void SampleHandlerTutorial::SetupWeightPointers() {
   }
 }
 
+void SampleHandlerTutorial::CleanMemoryBeforeFit() {
+  CleanVector(TutorialPlottingSamples);
+}
+
 // ************************************************
 int SampleHandlerTutorial::SetupExperimentMC() {
 // ************************************************
@@ -115,6 +124,7 @@ int SampleHandlerTutorial::SetupExperimentMC() {
   delete _Chain;
 
   TutorialSamples.resize(nEntries);
+  TutorialPlottingSamples.resize(nEntries);
 
   int TotalEventCounter = 0.;
   for(int iChannel = 0 ; iChannel < static_cast<int>(OscChannels.size()); iChannel++) {
@@ -182,8 +192,25 @@ int SampleHandlerTutorial::SetupExperimentMC() {
     */
 
     _data->GetEntry(0);
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    std::uniform_int_distribution<> unif(0,3);
+    std::normal_distribution<> mu_angle(0,M_PI/8);
+    std::normal_distribution<> pi_angle(0,M_PI/2);
+    std::uniform_real_distribution<> nucl_angle(-M_PI,M_PI);
+    std::exponential_distribution<> pi_energy(1./0.5);
+    std::exponential_distribution<> nucl_energy(1./2);
+
     for (int i = 0; i < nEntries; ++i) { // Loop through tree
       _data->GetEntry(i);
+
+      // === JM resize particle-level vectors ===
+      // JM: We don't have particle-level info in the tutorial sample, so will fake it
+      int nParticles = 5; //fake number of particles in event
+      TutorialPlottingSamples[TotalEventCounter].particle_energy.resize(nParticles);
+      TutorialPlottingSamples[TotalEventCounter].particle_pdg.resize(nParticles);
+      TutorialPlottingSamples[TotalEventCounter].particle_beamangle.resize(nParticles);
+      // ========================================
 
       TutorialSamples[TotalEventCounter].TrueEnu = Enu_true;
       // HH: We don't have Erec in the tutorial sample, so we set it to the true energy
@@ -206,6 +233,51 @@ int SampleHandlerTutorial::SetupExperimentMC() {
       if(isATM) {
         TutorialSamples[TotalEventCounter].TrueCosZenith = trueCZ;
       }
+
+      // === JM loop through particles in event ===
+      for (int iParticle = 0; iParticle < nParticles; ++iParticle) {
+        //JM: No particle-level data in sample, so fake it
+        if (iParticle==0) {
+          TutorialPlottingSamples[TotalEventCounter].particle_pdg[iParticle] = PDGLep;
+          TutorialPlottingSamples[TotalEventCounter].particle_energy[iParticle] = ELep;
+          TutorialPlottingSamples[TotalEventCounter].particle_beamangle[iParticle] = mu_angle(gen);
+        }
+        else {
+          int particle_seed = unif(gen);
+          double angle, energy;
+          int pdg;
+          switch (particle_seed) {
+            case 0: 
+              pdg = 211;
+              energy = pi_energy(gen);
+              angle = pi_angle(gen);
+              while (angle>M_PI || angle<-M_PI) angle = pi_angle(gen);
+              break;
+            case 1:
+              pdg = -211;
+              energy = pi_energy(gen);
+              angle = pi_angle(gen);
+              while (angle>M_PI || angle<-M_PI) angle = pi_angle(gen);
+              break;
+            case 2: 
+              pdg = 2212;
+              energy = nucl_energy(gen);
+              angle = nucl_angle(gen);
+              break;
+            case 3:
+              pdg = 2112;
+              energy = nucl_energy(gen);
+              angle = nucl_angle(gen);
+              break;
+            default:
+              break;
+          }
+          TutorialPlottingSamples[TotalEventCounter].particle_energy[iParticle] = energy;
+          TutorialPlottingSamples[TotalEventCounter].particle_beamangle[iParticle] = angle;
+          TutorialPlottingSamples[TotalEventCounter].particle_pdg[iParticle] = pdg;
+        }
+      }
+      // ==========================================
       TotalEventCounter++;
     }
     _sampleFile->Close();
@@ -230,6 +302,32 @@ double SampleHandlerTutorial::ReturnKinematicParameter(std::string KinematicPara
   return ReturnKinematicParameter(KinPar, iEvent);
 }
 
+// === JM Define ReturnKinematicVector functions === 
+std::vector<double> SampleHandlerTutorial::ReturnKinematicVector(KinematicParticleVecs KinVec, int iEvent) {
+  switch (KinVec) {
+    case kParticleEnergy:
+      return TutorialPlottingSamples[iEvent].particle_energy;
+    case kParticlePDG:
+      return TutorialPlottingSamples[iEvent].particle_pdg;
+    case kParticleBeamAngle:
+      return TutorialPlottingSamples[iEvent].particle_beamangle;
+    default:
+      MACH3LOG_ERROR("Unrecognized Kinematic Vector: {}", static_cast<int>(KinVec));
+      throw MaCh3Exception(__FILE__, __LINE__);
+  }
+}
+
+std::vector<double> SampleHandlerTutorial::ReturnKinematicVector(int KinematicVector, int iEvent) {
+  KinematicParticleVecs KinVec = static_cast<KinematicParticleVecs>(std::round(KinematicVector));
+  return ReturnKinematicVector(KinVec, iEvent);
+}
+
+std::vector<double> SampleHandlerTutorial::ReturnKinematicVector(std::string KinematicVector, int iEvent) {
+  KinematicParticleVecs KinVec = static_cast<KinematicParticleVecs>(ReturnKinematicVectorFromString(KinematicVector));
+  return ReturnKinematicVector(KinVec, iEvent);
+}
+// =================================================
+
 const double* SampleHandlerTutorial::GetPointerToKinematicParameter(KinematicTypes KinPar, int iEvent) {
   switch (KinPar) {
     case kTrueNeutrinoEnergy:
@@ -270,6 +368,8 @@ void SampleHandlerTutorial::SetupFDMC() {
 }
 
 std::vector<double> SampleHandlerTutorial::ReturnKinematicParameterBinning(std::string KinematicParameterStr) {
+  if (IsSubEventVarString(KinematicParameterStr)) return ReturnKinematicVectorBinning(KinematicParameterStr);
+
   std::vector<double> binningVector;
   KinematicTypes KinematicParameter = static_cast<KinematicTypes>(ReturnKinematicParameterFromString(KinematicParameterStr));
 
@@ -300,3 +400,40 @@ std::vector<double> SampleHandlerTutorial::ReturnKinematicParameterBinning(std::
 
   return binningVector;
 }
+std::vector<double> SampleHandlerTutorial::ReturnKinematicVectorBinning(std::string KinematicVectorStr) {
+  std::vector<double> binningVector;
+  KinematicParticleVecs KinematicVector = static_cast<KinematicParticleVecs>(ReturnKinematicVectorFromString(KinematicVectorStr));
+
+  int nBins = 0;
+  double bin_start = 0;
+  double bin_end = 0;
+
+  switch(KinematicVector){
+    case(kParticleEnergy):
+      nBins = 20;
+      bin_start = 0;
+      bin_end = 6;
+      break;
+    case(kParticlePDG):
+      nBins = 100;
+      bin_start = -2000;
+      bin_end = 2000;
+      break;
+    case(kParticleBeamAngle):
+      nBins = 20;
+      bin_start = -M_PI;
+      bin_end = M_PI;
+      break;
+    default:
+      MACH3LOG_ERROR("Binning for {} not found", KinematicVectorStr);
+      throw MaCh3Exception(__FILE__,__LINE__);
+      break;
+  }
+  double step = (bin_end - bin_start)/nBins;
+
+  for(double bin_i = bin_start; bin_i <= bin_end; bin_i+=step){
+    binningVector.push_back(bin_i);
+  }
+
+  return binningVector;
+}

SamplesTutorial/SampleHandlerTutorial.h

@@ -3,14 +3,20 @@
 #include "Samples/SampleHandlerFD.h"
 #include "StructsTutorial.h"
 #include "SplinesTutorial/BinnedSplinesTutorial.h"
+#include <random>
 
 class SampleHandlerTutorial : public SampleHandlerFD
 {
  public:
   SampleHandlerTutorial(const std::string& config_name, ParameterHandlerGeneric* parameter_handler,
                         const std::shared_ptr<OscillationHandler>& Oscillator_ = nullptr);
   virtual ~SampleHandlerTutorial();
+
   enum KinematicTypes {kTrueNeutrinoEnergy, kTrueQ2, kM3Mode, kRecoNeutrinoEnergy};
+
+  // === JM enum for particle-level parameters ===
+  enum KinematicParticleVecs {kParticleEnergy, kParticlePDG, kParticleBeamAngle};
+  // =============================================
 
   std::vector<double> ReturnKinematicParameterBinning(std::string KinematicParameter) override;
 
@@ -24,11 +30,18 @@ class SampleHandlerTutorial : public SampleHandlerFD
 
   int SetupExperimentMC() override;
 
-  void CleanMemoryBeforeFit() override {};
+  void CleanMemoryBeforeFit() override;
 
   double ReturnKinematicParameter(KinematicTypes KinPar, int iEvent);
   double ReturnKinematicParameter(int KinematicVariable, int iEvent);
   double ReturnKinematicParameter(std::string KinematicParameter, int iEvent);
+
+  // === JM ReturnKinematicVector declarations for particle-level parameters ===
+  std::vector<double> ReturnKinematicVector(KinematicParticleVecs KinVec, int iEvent);
+  std::vector<double> ReturnKinematicVector(int KinematicVector, int iEvent);
+  std::vector<double> ReturnKinematicVector(std::string KinematicVector, int iEvent);
+  std::vector<double> ReturnKinematicVectorBinning(std::string KinematicParameter);
+  // ===========================================================================
 
   const double* GetPointerToKinematicParameter(KinematicTypes KinPar, int iEvent);
   const double* GetPointerToKinematicParameter(std::string KinematicParameter, int iEvent) override;
@@ -38,6 +51,7 @@ class SampleHandlerTutorial : public SampleHandlerFD
   void CalcWeightFunc(int iEvent) override {return; (void)iEvent;}
 
   std::vector<TutorialMCInfo> TutorialSamples;
+  std::vector<TutorialMCPlottingInfo> TutorialPlottingSamples;
 
   const std::unordered_map<std::string, int> KinematicParametersTutorial = {
     {"TrueNeutrinoEnergy", kTrueNeutrinoEnergy},
@@ -52,6 +66,20 @@ class SampleHandlerTutorial : public SampleHandlerFD
     {kM3Mode,"Mode"},
     {kRecoNeutrinoEnergy, "RecoNeutrinoEnergy"},
   };
+
+  // === JM maps for particle-level parameters ===
+  const std::unordered_map<std::string, int> KinematicVectorsTutorial = {
+    {"ParticleEnergy", kParticleEnergy},
+    {"ParticlePDG", kParticlePDG},
+    {"ParticleBeamAngle", kParticleBeamAngle},
+  };
+
+  const std::unordered_map<int, std::string> ReversedKinematicVectorsTutorial = {
+    {kParticleEnergy, "ParticleEnergy"},
+    {kParticlePDG, "ParticlePDG"},
+    {kParticleBeamAngle, "kParticleBeamAngle"},
+  };
+  // =============================================
 
   double pot;
   bool isATM;

SamplesTutorial/StructsTutorial.h

@@ -31,3 +31,14 @@ struct TutorialMCInfo {
   /// Lepton energy
   double ELep = M3::_BAD_DOUBLE_;
 };
+
+struct TutorialMCPlottingInfo {
+  // === JM Particle-level kinematic parameters ===
+  /// particle energy
+  std::vector<double> particle_energy = {};
+  /// particle pdg
+  std::vector<double> particle_pdg = {};
+  /// particle angle to beam
+  std::vector<double> particle_beamangle = {};
+  // ==============================================
+};