New plot to monitor alignment weak modes added to ParticleMassMonitor

Some updates to the ParticleMassMonitor:

• Addition of a new profile histogram showing the difference between the average of the reconstructed mass distribution and the central value of the mass histogram in bins of the momentum difference between the positively and negatively charged daughters. Useful to monitor the curvature bias effect on alignment
• Generalisation of the monitor to remove the hard-coded particle decays to solve the issue: #506 (closed)

Goes with Moore!2818 (merged) and Alignment!436 (merged)

Phys/ParticleMaker/src/ParticleMassMonitor.cpp

@@ -35,20 +35,20 @@ public:
 
   StatusCode initialize() override {
     return Consumer::initialize().andThen( [&] {
-      /*Assign the names of the daughter particles of the decay. So far only considering D0 -> KPi, Z -> mumu and JPs ->
+      /*Assign the names of the daughter particles of the decay. So far only considering D0 -> PiK, Z -> mumu and JPs ->
        mumu used for alignment mass constraints*/
       if ( m_mother == "D0" ) {
-        m_daughter1 = "K";
-        m_daughter2 = "Pi";
+        m_daughter1 = "Pi";
+        m_daughter2 = "K";
       } else {
         m_daughter1 = "Mu1";
         m_daughter2 = "Mu2";
       }
 
       // Mass histograms
-      m_histograms_mass.emplace(
-          std::piecewise_construct, std::forward_as_tuple( 0 ),
-          std::forward_as_tuple( this, m_bins, m_mother, m_mass, m_sigma, m_minpt, m_maxpt, m_minp, m_maxp ) );
+      m_histograms_mass.emplace( std::piecewise_construct, std::forward_as_tuple( 0 ),
+                                 std::forward_as_tuple( this, m_bins, m_pdiff_bins, m_mother, m_mass, m_sigma, m_minpt,
+                                                        m_maxpt, m_minp, m_maxp, m_pdiff_range / Gaudi::Units::GeV ) );
 
       // Momentum histograms
       std::vector<std::string> particles = {m_mother, m_daughter1, m_daughter2};
@@ -65,13 +65,15 @@ public:
 private:
   // Histogram settings. Default values are the standard ones for D0 -> KPi candidates
   Gaudi::Property<unsigned int> m_bins{this, "Bins", 100};
+  Gaudi::Property<unsigned int> m_pdiff_bins{this, "PdiffBins", 20};
   Gaudi::Property<std::string>  m_mother{this, "ParticleType", "D0"};
   Gaudi::Property<float>        m_mass{this, "MassMean", 1864.84 * Gaudi::Units::MeV};
-  Gaudi::Property<float>        m_sigma{this, "MassSigma", 90 * Gaudi::Units::MeV};
+  Gaudi::Property<float>        m_sigma{this, "MassSigma", 18 * Gaudi::Units::MeV};
   Gaudi::Property<float>        m_minpt{this, "MinPt", 800 * Gaudi::Units::MeV};
-  Gaudi::Property<float>        m_maxpt{this, "MaxPt", 6000 * Gaudi::Units::MeV};
+  Gaudi::Property<float>        m_maxpt{this, "MaxPt", 5000 * Gaudi::Units::MeV};
   Gaudi::Property<float>        m_minp{this, "MinP", 2000 * Gaudi::Units::MeV};
   Gaudi::Property<float>        m_maxp{this, "MaxP", 150000 * Gaudi::Units::MeV};
+  Gaudi::Property<float>        m_pdiff_range{this, "PdiffDelta", 40000 * Gaudi::Units::MeV};
 
   // Flag for tuple generation
   Gaudi::Property<bool> m_generateTuples{this, "GenerateTuples", false};
@@ -92,12 +94,14 @@ private:
   };
 
   struct HistoMass {
-    mutable Gaudi::Accumulators::Histogram<1> CandidateMass;
-    mutable Gaudi::Accumulators::Histogram<2> MassvsPt;
-    mutable Gaudi::Accumulators::Histogram<2> MassvsP;
-
-    HistoMass( const ParticleMassMonitor* owner, unsigned int const& bins, std::string const& name, float const& mass,
-               float const& sigma, float const& min_pt, float const& max_pt, float const& min_p, float const& max_p )
+    mutable Gaudi::Accumulators::Histogram<1>        CandidateMass;
+    mutable Gaudi::Accumulators::Histogram<2>        MassvsPt;
+    mutable Gaudi::Accumulators::Histogram<2>        MassvsP;
+    mutable Gaudi::Accumulators::ProfileHistogram<1> MassvsPdiff;
+
+    HistoMass( const ParticleMassMonitor* owner, unsigned int const& bins, unsigned int const& bins_pdiff,
+               std::string const& name, float const& mass, float const& sigma, float const& min_pt, float const& max_pt,
+               float const& min_p, float const& max_p, float const& range_pdiff )
         : CandidateMass{owner,
                         name + "_mass",
                         name + " mass",
@@ -109,12 +113,19 @@ private:
         , MassvsP{owner,
                   name + "_MassvsP",
                   name + " mass vs p",
-                  {{bins, min_p, max_p, "p [MeV/c]"}, {bins, mass - 5 * sigma, mass + 5 * sigma, "m [MeV/c^{2}]"}}} {}
+                  {{bins, min_p, max_p, "p [MeV/c]"}, {bins, mass - 5 * sigma, mass + 5 * sigma, "m [MeV/c^{2}]"}}}
+        , MassvsPdiff{owner,
+                      name + "_MassvsPdiff",
+                      fmt::format( "m({}) - {:.2f} MeV/c^{{2}} vs daughter p difference", name, mass ),
+                      {bins_pdiff, -range_pdiff, range_pdiff, "Daughter p_{+} - p_{-} [GeV/c]"}} {}
   };
 
   // Histogram objects are stored inside of maps
   std::map<int, HistoMomentum> m_histograms_momentum;
   std::map<int, HistoMass>     m_histograms_mass;
+
+  // Function to calculate the momentum difference of the daughter particles
+  double calculate_momentum_difference( LHCb::Particle::ConstVector const& particles ) const;
 };
 
 DECLARE_COMPONENT( ParticleMassMonitor )
@@ -122,16 +133,18 @@ DECLARE_COMPONENT( ParticleMassMonitor )
 // Algorithm
 void ParticleMassMonitor::operator()( LHCb::Particle::Range const& particles ) const {
   for ( const auto& p : particles ) {
-    const double m         = p->momentum().M();
-    const auto   daughters = p->daughtersVector();
-    const double pt        = p->pt();
-    const double mom       = p->p();
+    const double m              = p->momentum().M();
+    const auto   daughters      = p->daughtersVector();
+    const double pt             = p->pt();
+    const double mom            = p->p();
+    const double mom_difference = calculate_momentum_difference( daughters );
 
     // Mass plots
     auto& mass_histo = m_histograms_mass.at( 0 );
     ++mass_histo.CandidateMass[m];
     ++mass_histo.MassvsPt[{pt, m}];
     ++mass_histo.MassvsP[{mom, m}];
+    mass_histo.MassvsPdiff[mom_difference / Gaudi::Units::GeV] += m - m_mass; // Momentum in GeV for this plot
 
     // Momentum plots
     int daughters_counter = 0;
@@ -159,6 +172,7 @@ void ParticleMassMonitor::operator()( LHCb::Particle::Range const& particles ) c
       // Creation of the tuple
       Tuple particleTuple = nTuple( "Particletuples", "" );
       particleTuple->column( m_mother + "_M", m ).ignore();
+      particleTuple->column( "Daughter_p_plus_p_minus", mom_difference ).ignore();
       particleTuple->column( m_mother + "_Pt", pt ).ignore();
       particleTuple->column( m_daughter1 + "_Pt", daughters[0]->pt() ).ignore();
       particleTuple->column( m_daughter2 + "_Pt", daughters[1]->pt() ).ignore();
@@ -175,4 +189,13 @@ void ParticleMassMonitor::operator()( LHCb::Particle::Range const& particles ) c
       particleTuple->write().ignore();
     }
   }
+}
+
+double ParticleMassMonitor::calculate_momentum_difference( LHCb::Particle::ConstVector const& particles ) const {
+  const int    charge_p1 = particles[0]->charge();
+  const int    charge_p2 = particles[1]->charge();
+  const double p1        = particles[0]->p();
+  const double p2        = particles[1]->p();
+  const double p_diff    = ( charge_p1 > 0 && charge_p2 < 0 ) ? p1 - p2 : p2 - p1;
+  return p_diff;
 }
\ No newline at end of file