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,17 @@ public:
 
   StatusCode initialize() override {
     return Consumer::initialize().andThen( [&] {
-      /*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 = "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_pdiff_bins, m_mother, m_mass, m_sigma, m_minpt,
+                                 std::forward_as_tuple( this, m_bins, m_pdiff_bins, m_particles[0], 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++ ) {
+      for ( unsigned int i = 0; i < m_particles.size(); i++ ) {
         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 ) );
       }
     } );
   }
@@ -64,24 +54,22 @@ 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", 18 * Gaudi::Units::MeV};
-  Gaudi::Property<float>        m_minpt{this, "MinPt", 800 * 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};
+  Gaudi::Property<unsigned int>              m_bins{this, "Bins", 100};
+  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", 18 * Gaudi::Units::MeV};
+  Gaudi::Property<float>                     m_minpt{this, "MinPt", 800 * 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;
@@ -144,58 +132,57 @@ void ParticleMassMonitor::operator()( LHCb::Particle::Range const& particles ) c
     ++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
+    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
+      const auto id = p->particleID().pid();
       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();
-      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();
+      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( 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;
+  if ( particles.size() == 2 ) {
+    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;
+  } else {
+    return 0.;
+  }
 }
\ No newline at end of file