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,27 +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 ->
-       mumu used for alignment mass constraints*/
-      if ( m_mother == "D0" ) {
-        m_daughter1 = "K";
-        m_daughter2 = "Pi";
-      } else {
-        m_daughter1 = "Mu1";
-        m_daughter2 = "Mu2";
-      }
+      for ( unsigned int i = 0; i < m_particles.size(); i++ ) {
 
-      // 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 ) );
+        // Mass histogram
+        if ( i == 0 ) {
+          m_histograms_mass.emplace( std::piecewise_construct, std::forward_as_tuple( 0 ),
+                                     std::forward_as_tuple( this, m_bins, m_pdiff_bins, m_particles[i], 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};
-      for ( int i = 0; i < 3; i++ ) {
+        // Momentum histograms
         m_histograms_momentum.emplace(
             std::piecewise_construct, std::forward_as_tuple( i ),
-            std::forward_as_tuple( this, m_bins, particles[i], m_minpt, m_maxpt, m_minp, m_maxp ) );
+            std::forward_as_tuple( this, m_bins, m_particles[i], m_minpt, m_maxpt, m_minp, m_maxp ) );
       }
     } );
   }
@@ -65,21 +58,21 @@ 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<std::string>  m_mother{this, "ParticleType", "D0"};
+  Gaudi::Property<unsigned int> m_pdiff_bins{this, "PdiffBins", 20};
   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};
+
+  // Vector with the name of the particles in the same order as in the particle combiner
+  Gaudi::Property<std::vector<std::string>> m_particles{this, "ParticleType", {"D0", "Pi", "K"}};
 
   // Flag for tuple generation
   Gaudi::Property<bool> m_generateTuples{this, "GenerateTuples", false};
 
-  // Daughter particles
-  std::string m_daughter1;
-  std::string m_daughter2;
-
   // Definition of histograms
   struct HistoMomentum {
     mutable Gaudi::Accumulators::Histogram<1> CandidatePt;
@@ -92,12 +85,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 +104,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( const LHCb::Particle& particle1, const LHCb::Particle& particle2 ) const;
 };
 
 DECLARE_COMPONENT( ParticleMassMonitor )
@@ -126,53 +128,59 @@ void ParticleMassMonitor::operator()( LHCb::Particle::Range const& particles ) c
     const auto   daughters = p->daughtersVector();
     const double pt        = p->pt();
     const double mom       = p->p();
+    const double mom_difference =
+        daughters.size() > 2 ? 0. : calculate_momentum_difference( *daughters[0], *daughters[1] );
 
     // Mass plots
     auto& mass_histo = m_histograms_mass.at( 0 );
     ++mass_histo.CandidateMass[m];
     ++mass_histo.MassvsPt[{pt, m}];
     ++mass_histo.MassvsP[{mom, m}];
+    if ( m_particles.size() == 3 ) {
+      mass_histo.MassvsPdiff[mom_difference / Gaudi::Units::GeV] += m - m_mass; // Momentum in GeV for this plot
+    }
 
     // Momentum plots
-    int daughters_counter = 0;
-    for ( int i = 0; i < 3; i++ ) {
+    for ( unsigned int i = 0; i < m_particles.size(); i++ ) {
       auto& momentum_histo = m_histograms_momentum.at( i );
       if ( i == 0 ) {
         ++momentum_histo.CandidatePt[pt];
         ++momentum_histo.CandidateMom[mom];
       } else {
-        ++momentum_histo.CandidatePt[daughters[daughters_counter]->pt()];
-        ++momentum_histo.CandidateMom[daughters[daughters_counter]->p()];
-        daughters_counter++;
+        ++momentum_histo.CandidatePt[daughters[i - 1]->pt()];
+        ++momentum_histo.CandidateMom[daughters[i - 1]->p()];
       }
     }
 
     // Code for tuple generation
     if ( m_generateTuples ) {
-      // Extra variables to store in the tuples
-      const auto slopes1 = daughters[0]->slopes();
-      const auto slopes2 = daughters[1]->slopes();
-      const auto id      = p->particleID().pid();
-      const auto id1     = daughters[0]->particleID().pid();
-      const auto id2     = daughters[1]->particleID().pid();
-
-      // Creation of the tuple
-      Tuple particleTuple = nTuple( "Particletuples", "" );
-      particleTuple->column( m_mother + "_M", m ).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();
-      particleTuple->column( m_mother + "_ID", id ).ignore();
-      particleTuple->column( m_daughter1 + "_ID", id1 ).ignore();
-      particleTuple->column( m_daughter2 + "_ID", id2 ).ignore();
-      particleTuple->column( m_mother + "_P", mom ).ignore();
-      particleTuple->column( m_daughter1 + "_P", daughters[0]->p() ).ignore();
-      particleTuple->column( m_daughter2 + "_P", daughters[1]->p() ).ignore();
-      particleTuple->column( m_daughter1 + "_PxoverPz", slopes1.X() ).ignore();
-      particleTuple->column( m_daughter2 + "_PxoverPz", slopes2.X() ).ignore();
-      particleTuple->column( m_daughter1 + "_PyoverPz", slopes1.Y() ).ignore();
-      particleTuple->column( m_daughter2 + "_PyoverPz", slopes2.Y() ).ignore();
+      const auto id            = p->particleID().pid();
+      Tuple      particleTuple = nTuple( "Particletuples", "" );
+      particleTuple->column( m_particles[0] + "_M", m ).ignore();
+      particleTuple->column( m_particles[0] + "_Pt", pt ).ignore();
+      particleTuple->column( m_particles[0] + "_ID", id ).ignore();
+      particleTuple->column( m_particles[0] + "_P", mom ).ignore();
+
+      for ( unsigned int i = 1; i < m_particles.size(); i++ ) {
+        particleTuple->column( m_particles[i] + "_Pt", daughters[i - 1]->pt() ).ignore();
+        particleTuple->column( m_particles[i] + "_ID", daughters[i - 1]->particleID().pid() ).ignore();
+        particleTuple->column( m_particles[i] + "_P", daughters[i - 1]->p() ).ignore();
+        particleTuple->column( m_particles[i] + "_PxoverPz", daughters[i - 1]->slopes().X() ).ignore();
+        particleTuple->column( m_particles[i] + "_PyoverPz", daughters[i - 1]->slopes().Y() ).ignore();
+      }
+
+      if ( m_particles.size() == 3 ) { particleTuple->column( "Daughter_p_plus_p_minus", mom_difference ).ignore(); }
       particleTuple->write().ignore();
     }
   }
+}
+
+double ParticleMassMonitor::calculate_momentum_difference( const LHCb::Particle& particle1,
+                                                           const LHCb::Particle& particle2 ) const {
+  const int    charge_p1 = particle1.charge();
+  const int    charge_p2 = particle2.charge();
+  const double p1        = particle1.p();
+  const double p2        = particle2.p();
+  const double p_diff    = ( charge_p1 > 0 && charge_p2 < 0 ) ? p1 - p2 : p2 - p1;
+  return p_diff;
 }
\ No newline at end of file