Plume timing hists

Plume/PlumeReco/src/PlumeDigitMonitor.cpp

@@ -30,6 +30,12 @@
 #include "AIDA/IHistogram2D.h"
 #include "AIDA/IProfile1D.h"
 
+// ROOT
+#include <TF1.h>
+#include <TFitResult.h>
+#include <TGraph.h>
+#include <TGraphErrors.h>
+
 #include <TProfile.h>
 #include <sstream>
 
@@ -38,7 +44,6 @@ using Input   = LHCb::PlumeAdcs;
 using BXTypes = LHCb::ODIN::BXTypes;
 
 namespace {
-
   template <typename T>
   using axis1D = Gaudi::Accumulators::Axis<T>;
 
@@ -48,6 +53,13 @@ namespace {
     t.emplace( std::piecewise_construct, std::forward_as_tuple( key ),
                std::forward_as_tuple( owner, name, title, axis1 ) );
   }
+  template <typename T, typename K, typename OWNER>
+  void map_emplace2D( T& t, K&& key, OWNER* owner, std::string const& name, std::string const& title,
+                      Gaudi::Accumulators::Axis<typename T::mapped_type::AxisArithmeticType> axis1,
+                      Gaudi::Accumulators::Axis<typename T::mapped_type::AxisArithmeticType> axis2 ) {
+    t.emplace( std::piecewise_construct, std::forward_as_tuple( key ),
+               std::forward_as_tuple( owner, name, title, axis1, axis2 ) );
+  }
 
   template <typename T>
   std::string to_string( const T& streamable ) {
@@ -58,6 +70,10 @@ namespace {
 
   const auto AxisBCID     = axis1D<double>( 3565, -0.5, 3564.5 );
   const auto AxisADC      = axis1D<double>( 4096 + 256, -256, 4096 );
+  const auto AxisTIME     = axis1D<double>( 100, 0, 8 );
+  const auto AxisTIMEERR  = axis1D<double>( 100, 0, 0.4 );
+  const auto AxisSAMP     = axis1D<double>( 1000, 0, 1000 );
+  const auto AxisSAMPERR  = axis1D<double>( 100, 0, 10 );
   const auto AxisINDEXTAE = axis1D<double>( 9, 0.5, 9.5 );
 } // namespace
 
@@ -84,6 +100,14 @@ private:
 
   void monitor( const Input& digits, const LHCb::ODIN& odin ) const;
 
+  // timing s-shape
+  void fit_sshapes( const Input& adcs, std::map<unsigned int, std::vector<std::pair<float, float>>>& result,
+                    std::map<unsigned int, int>& status ) const;
+  void fit_single_sshape( std::vector<int> sshape, std::vector<std::pair<float, float>>& result, int& status ) const;
+  void initialize_function();
+  std::unique_ptr<TF1>   m_fit_function_tf1;
+  Gaudi::Property<float> m_threshold_sshape{this, "ThresholdSShape", 150.5};
+
   mutable std::mutex m_lazy_lock;
 
   mutable Gaudi::Accumulators::Histogram<1> m_over_thld_lumi{
@@ -157,6 +181,16 @@ private:
       this, "hits_back_h", "Number of hits in back layer horizontal PMTs",
       axis1D<double>{10, -5., 5., "hits_back_h_axis", {"46", "45", "44", "43", "42", "30", "31", "32", "33", "34"}}};
 
+  mutable std::map<unsigned int, Gaudi::Accumulators::Histogram<1>> m_time_channel_tot;
+  mutable std::map<unsigned int, Gaudi::Accumulators::Histogram<1>> m_time_err_channel_tot;
+  mutable std::map<unsigned int, Gaudi::Accumulators::Histogram<1>> m_sshape_amp_tot;
+  mutable std::map<unsigned int, Gaudi::Accumulators::Histogram<1>> m_sshape_amp_err_tot;
+
+  mutable std::map<unsigned int, Gaudi::Accumulators::Histogram<2>> m_time_channel;
+  mutable std::map<unsigned int, Gaudi::Accumulators::Histogram<2>> m_time_err_channel;
+  mutable std::map<unsigned int, Gaudi::Accumulators::Histogram<2>> m_sshape_amp;
+  mutable std::map<unsigned int, Gaudi::Accumulators::Histogram<2>> m_sshape_amp_err;
+
   mutable Gaudi::Accumulators::StatCounter<> m_numAdcOTLumi = {this, "Number of ADC over threshold for LUMI"};
   mutable Gaudi::Accumulators::StatCounter<> m_numCoinc     = {this, "Number of coincidences"};
 
@@ -236,6 +270,28 @@ StatusCode PlumeDigitMonitor::initialize() {
       map_emplace( m_mon_bx_adc, bx, this, "adc_MON_" + to_string( bx ), "ADC for MON/" + to_string( bx ), AxisADC );
       map_emplace( m_pin_bx_adc, bx, this, "adc_PIN_" + to_string( bx ), "ADC for PIN/" + to_string( bx ), AxisADC );
     }
+
+    for ( auto tchannel : {5, 11, 29, 35} ) {
+      map_emplace( m_time_channel_tot, tchannel, this, "time_sshape_tot_" + to_string( tchannel ),
+                   "Inflection point for " + to_string( tchannel ), AxisTIME );
+      map_emplace( m_time_err_channel_tot, tchannel, this, "time_sshape_err_tot_" + to_string( tchannel ),
+                   "Uncertainty per track for " + to_string( tchannel ), AxisTIMEERR );
+      map_emplace( m_sshape_amp_tot, tchannel, this, "amp_sshape_tot_" + to_string( tchannel ),
+                   "S-shape amplitude for " + to_string( tchannel ), AxisSAMP );
+      map_emplace( m_sshape_amp_err_tot, tchannel, this, "amp_sshape_err_tot_" + to_string( tchannel ),
+                   "S-shape amplitude error for " + to_string( tchannel ), AxisSAMPERR );
+
+      map_emplace2D( m_time_channel, tchannel, this, "time_sshape_" + to_string( tchannel ),
+                     "Inflection point for " + to_string( tchannel ), AxisBCID, AxisTIME );
+      map_emplace2D( m_time_err_channel, tchannel, this, "time_sshape_err_" + to_string( tchannel ),
+                     "Uncertainty per track for " + to_string( tchannel ), AxisBCID, AxisTIMEERR );
+      map_emplace2D( m_sshape_amp, tchannel, this, "amp_sshape_" + to_string( tchannel ),
+                     "S-shape amplitude for " + to_string( tchannel ), AxisBCID, AxisSAMP );
+      map_emplace2D( m_sshape_amp_err, tchannel, this, "amp_sshape_err_" + to_string( tchannel ),
+                     "S-shape amplitude error for " + to_string( tchannel ), AxisBCID, AxisSAMPERR );
+    }
+
+    initialize_function();
   } );
 }
 
@@ -271,12 +327,17 @@ void PlumeDigitMonitor::monitor( const Input& adcs, const LHCb::ODIN& odin ) con
 
   ++m_calibType[calibType];
 
+  std::map<unsigned int, std::vector<std::pair<float, float>>> sshape_result{};
+  std::map<unsigned int, int>                                  sshape_status{};
+  fit_sshapes( adcs, sshape_result, sshape_status );
+
   for ( const auto& digit : adcs ) {
-    const auto adc          = digit->adc();
-    const auto type         = digit->channelID().channelType();
-    const auto channel_id   = digit->channelID().channelID();
-    const auto over_thld    = overThreshold( *digit );
-    bool       complementOT = false;
+    const auto adc            = digit->adc();
+    const auto type           = digit->channelID().channelType();
+    const auto channel_id     = digit->channelID().channelID();
+    const auto channel_sub_id = digit->channelID().channelSubID();
+    const auto over_thld      = overThreshold( *digit );
+    bool       complementOT   = false;
     if ( type == ChannelType::LUMI || type == ChannelType::TIME ) {
       const auto* complement = adcs.object( ChannelID(
           type, channel_id < 24 ? ( channel_id + 24 ) : ( channel_id - 24 ), digit->channelID().channelSubID() ) );
@@ -290,6 +351,20 @@ void PlumeDigitMonitor::monitor( const Input& adcs, const LHCb::ODIN& odin ) con
       ++m_channel_bx_adc.at( key )[adc];
     }
 
+    if ( type == ChannelType::TIME && channel_sub_id == 0 ) {
+      if ( sshape_status.at( channel_id ) == true ) {
+        ++m_time_channel.at( channel_id )[{bcid, sshape_result.at( channel_id )[3].first}];
+        ++m_time_err_channel.at( channel_id )[{bcid, sshape_result.at( channel_id )[3].second}];
+        ++m_sshape_amp.at( channel_id )[{bcid, sshape_result.at( channel_id )[1].first}];
+        ++m_sshape_amp_err.at( channel_id )[{bcid, sshape_result.at( channel_id )[1].second}];
+
+        ++m_time_channel_tot.at( channel_id )[sshape_result.at( channel_id )[3].first];
+        ++m_time_err_channel_tot.at( channel_id )[sshape_result.at( channel_id )[3].second];
+        ++m_sshape_amp_tot.at( channel_id )[sshape_result.at( channel_id )[1].first];
+        ++m_sshape_amp_err_tot.at( channel_id )[sshape_result.at( channel_id )[1].second];
+      }
+    }
+
     switch ( type ) {
     case ChannelType::LUMI:
       nCoincLumi += over_thld && complementOT;
@@ -405,3 +480,83 @@ void PlumeDigitMonitor::monitor( const Input& adcs, const LHCb::ODIN& odin ) con
   m_bcid_anycoinc_lumi[bcid] += nCoincLumi > 0;
   m_numAdcOTLumi += nOTLumi;
 }
+
+void PlumeDigitMonitor::initialize_function() {
+  m_fit_function_tf1 = std::make_unique<TF1>( "sshape_tf1", "[0] + [1]*TMath::Erf( [2]*(x - [3]) )", 0., 8. );
+  m_fit_function_tf1->SetParNames( "A", "B", "a", "t" );
+  m_fit_function_tf1->SetParLimits( 0, -1000, 1000. );
+  m_fit_function_tf1->SetParLimits( 1, 0., 1000. );
+  m_fit_function_tf1->SetParLimits( 2, 0., 10. );
+  m_fit_function_tf1->SetParLimits( 3, 0., 10. );
+  m_fit_function_tf1->SetParameters( 30., 300., 0.4, 4. );
+}
+
+void PlumeDigitMonitor::fit_single_sshape( std::vector<int> sshape, std::vector<std::pair<float, float>>& result,
+                                           int& status ) const {
+  TGraphErrors sshape_graph{};
+  for ( const auto point : sshape ) {
+    sshape_graph.SetPoint( sshape_graph.GetN(), sshape_graph.GetN(), point );
+    sshape_graph.SetPointError( sshape_graph.GetN() - 1, 0, std::sqrt( std::abs( point ) ) );
+  }
+  sshape_graph.Fit( m_fit_function_tf1.get(), "NQ" );
+  auto fitResult = sshape_graph.Fit( m_fit_function_tf1.get(), "SNQ" );
+
+  status = fitResult->IsValid();
+
+  result.emplace_back( m_fit_function_tf1->GetParameter( 0 ), m_fit_function_tf1->GetParError( 0 ) );
+  result.emplace_back( m_fit_function_tf1->GetParameter( 1 ), m_fit_function_tf1->GetParError( 1 ) );
+  result.emplace_back( m_fit_function_tf1->GetParameter( 2 ), m_fit_function_tf1->GetParError( 2 ) );
+  result.emplace_back( m_fit_function_tf1->GetParameter( 3 ), m_fit_function_tf1->GetParError( 3 ) );
+}
+
+void PlumeDigitMonitor::fit_sshapes( const Input&                                                  adcs,
+                                     std::map<unsigned int, std::vector<std::pair<float, float>>>& result,
+                                     std::map<unsigned int, int>&                                  status ) const {
+  // map<<channelID, map<subchannelID, ADC>>
+  auto timing_map = std::map<unsigned int, std::map<unsigned int, int>>{};
+  // map<<channelID, vector<ADC>>
+  // where vector is ordered by s-shape order
+  auto sshapes = std::map<unsigned int, std::vector<int>>{};
+
+  // fill the timing_map with the ADCs
+  for ( const auto& digit : adcs ) {
+    if ( digit->channelID().channelType() != ChannelID::ChannelType::TIME ) continue;
+    const auto channelID    = digit->channelID().channelID();
+    const auto channelSubID = digit->channelID().channelSubID();
+
+    if ( timing_map.find( channelID ) == timing_map.end() ) timing_map[channelID] = std::map<unsigned int, int>{};
+    timing_map[channelID][channelSubID] = digit->adc();
+  }
+
+  std::map<unsigned int, bool> coincidence{};
+
+  // use the last s-shape point for checking for coincidence
+  coincidence[5] = coincidence[29] =
+      ( timing_map.at( 5 ).at( 15 ) - timing_map.at( 5 ).at( 7 ) ) > m_threshold_sshape &&
+      ( timing_map.at( 29 ).at( 15 ) - timing_map.at( 29 ).at( 7 ) ) > m_threshold_sshape;
+  coincidence[11] = coincidence[35] =
+      ( timing_map.at( 11 ).at( 15 ) - timing_map.at( 11 ).at( 7 ) ) > m_threshold_sshape &&
+      ( timing_map.at( 35 ).at( 15 ) - timing_map.at( 35 ).at( 7 ) ) > m_threshold_sshape;
+
+  // fill the s-shape map with differences
+  for ( const auto& tchannel : timing_map ) {
+    sshapes[tchannel.first] = std::vector<int>{};
+    for ( unsigned int tsubchannel = 0; tsubchannel < 8; ++tsubchannel ) {
+      sshapes.at( tchannel.first )
+          .push_back( tchannel.second.at( tsubchannel + 8 ) - tchannel.second.at( tsubchannel ) );
+    }
+  }
+
+  // fit the s-shapes
+  for ( const auto& sshape : sshapes ) {
+    const auto                           channelID = sshape.first;
+    std::vector<std::pair<float, float>> fit_result;
+    int                                  status_sshape = -1;
+    if ( coincidence[channelID] ) {
+      fit_single_sshape( sshape.second, fit_result, status_sshape );
+      result[channelID] = fit_result;
+    } else
+      result[channelID] = std::vector<std::pair<float, float>>( 4, std::pair<float, float>{-9999, -9999} );
+    status[channelID] = status_sshape;
+  }
+}

Plume/PlumeReco/src/PlumeRawToDigits.cpp

@@ -50,7 +50,9 @@ private:
   mutable Gaudi::Accumulators::MsgCounter<MSG::WARNING> m_noMonitoringBank{this, "No monitoring bank found for Plume"};
   mutable Gaudi::Accumulators::MsgCounter<MSG::WARNING> m_noTimingBank{this, "No timing bank found for Plume"};
   mutable Gaudi::Accumulators::MsgCounter<MSG::WARNING> m_tooManyBanks{this, "Too many banks for Plume"};
-  mutable Gaudi::Accumulators::WeightedHistogram<1>     m_sumADCPerChannel{
+  mutable Gaudi::Accumulators::MsgCounter<MSG::ERROR>   m_unknown_version{
+      this, "Plume raw bank version not known. Assuming version = 1"};
+  mutable Gaudi::Accumulators::WeightedHistogram<1> m_sumADCPerChannel{
       this, "SumADCPerChannel", "Sum of ADCs per channel", {4 * 32, 0, 4 * 32}};
   mutable Gaudi::Accumulators::Histogram<1, Gaudi::Accumulators::atomicity::full, uint64_t> m_nOverThresholdPerChannel{
       this, "NOverThresholdPerChannel", "Number of set over-threshold bits per channel", {4 * 32, 0, 4 * 32}};
@@ -83,9 +85,7 @@ std::tuple<LHCb::PlumeAdcs, LHCb::PlumeCoincidences> PlumeRawToDigits::
   bool foundTimingBank     = false;
 
   for ( const auto& bank : banks ) {
-    if ( bank->version() > 1 ) {
-      error() << fmt::format( "Plume raw data version {} not known. Assuming version = 1", bank->version() ) << endmsg;
-    }
+    if ( bank->version() > 1 ) { ++m_unknown_version; }
     const int sourceID = bank->sourceID();
 
     std::bitset<64> bits_ovt{__builtin_bswap32( bank->range<uint32_t>()[0] ) +