MdMuonCounter.cc

Go to the documentation of this file.

#include "MdMuonCounter.h"

#ifndef DEBUG
#  define DEBUG 1
#endif

#include "VHalfMuonCounterFunctor.h"
#include "GapStrategy.h"
#include "AmountInWindowStrategy.h"
#include "ConsecutiveInWindowStrategy.h"
#include "AmountGlobalStrategy.h"
#include "../MdLDFFinderAG/ProfLike.h"

#include <fwk/CentralConfig.h>
#include <det/VManager.h>
#include <mdet/MDetector.h>
#include <sdet/SDetector.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>
#include <sevt/SEvent.h>
#include <mevt/Module.h>

#include <utl/ConsecutiveEnumFactory.h>
#include <utl/Trace.h>
#include <utl/AugerUnits.h>
#include <utl/Branch.h>
#include <utl/RandomEngine.h>
#include <utl/ConfigUtils.h>
#include <utl/TimeStamp.h>
#include <utl/PhysicalConstants.h>
#include <utl/TimeInterval.h>

#include <utl/ErrorLogger.h>
#include <utl/Verbosity.h>

#include <sstream>
//#include <cassert>

#include <TH1I.h>

using namespace fwk;
using std::cout;
using std::endl;


namespace MdMuonCounterAG {

  const char* const MdMuonCounter::kStrategyTags[] = { "Gap", "GapRelax", "AmountInWindow", "Consecutive", "AmountGlobal" };


  VModule::ResultFlag
  MdMuonCounter::Init()
  {
    // Configuration reading.
    utl::Branch config = fwk::CentralConfig::GetInstance()->GetTopBranch("MdMuonCounter");
    fUnits.SetTimeDefault(utl::ns, "ns");
    fUnits.Configure(config);
    fLog.Configure(config);
    std::string tag;
    LoadConfig(config, "strategy", tag, kStrategyTags[0]);
    fStrategy = utl::ConsecutiveEnumFactory<Strategy, eAmountGlobal, kStrategyTags>::Create(tag);

    // Init according to selected strategy.
    switch (fStrategy) {
    case eAmountInWindow:
      LoadConfig(config, "windowSize", fWindowSize, 8);
      LoadConfig(config, "nOnes", fNOnes, 2);
      fCounter.reset(new AmountInWindowStrategy(fWindowSize, fNOnes));
      fLog("#Window size", 1)
        (fWindowSize)("-")
        ("Number of ones")(fNOnes)("-")
        ("Strategy")(kStrategyTags[ fStrategy ])(".");
      break;
    case eConsecutive:
      LoadConfig(config, "windowSize", fWindowSize, 8);
      LoadConfig(config, "nOnes", fNOnes, 2);
      fCounter.reset(new ConsecutiveInWindowStrategy(fWindowSize, fNOnes));
      fLog("#Window size", 1)
        (fWindowSize)("-")
        ("Number of ones")(fNOnes)("-")
        ("Strategy")(kStrategyTags[ fStrategy ])(".");
      break;
    case eGap:
      LoadConfig(config, "windowSize", fWindowSize, 8);
      LoadConfig(config, "nOnes", fNOnes, 2);
      LoadConfig(config, "nGaps", fNGaps, 1);
      fCounter.reset(new GapStrategy(fWindowSize, fNOnes, fNGaps, "strict"));
      fLog("#Window size", 1)
        (fWindowSize)("-")
        ("Number of ones")(fNOnes)("-")
        ("Number of gaps")(fNGaps)("-")
        ("Strategy")(kStrategyTags[ fStrategy ])(".");
      break;
    case eGapRelax:
      LoadConfig(config, "windowSize", fWindowSize, 8);
      LoadConfig(config, "nOnes", fNOnes, 2);//minumum number of 1s
      LoadConfig(config, "nGaps", fNGaps, 6);//maximum number of gaps (0s) allowed between two 1s
      fCounter.reset(new GapStrategy(fWindowSize, fNOnes, fNGaps, "relax"));
      fLog("#Window size", 1)
        (fWindowSize)("-")
        ("Number of ones")(fNOnes)("-")
        ("Number of gaps")(fNGaps)("-")
        ("Strategy")(kStrategyTags[ fStrategy ])(".");
      break;
    case eAmountGlobal:
      LoadConfig(config, "nOnesPerPatternMatch", fNOnesPerPatternMatch, 2.0);
      LoadConfig(config, "nMinOnes", fNMinOnes, 0);
      fCounter.reset(new AmountGlobalStrategy(fNOnesPerPatternMatch, fNMinOnes));
      fWindowSize = 0;
      fLog("#Number of ones", 1)
        (fNOnesPerPatternMatch)("-")
        ("Minimum Number of ones in a scintillator")(fNMinOnes)("-")
        ("Strategy")(kStrategyTags[ fStrategy ])(".");
      break;
    }

    LoadConfig(config, "setTimeStamps", fSetTimeStamps, false);
    if (fSetTimeStamps)
      INFO("Setting the MD timestamps");
    else
      INFO("Not setting the MD timestamps");

    //utl::Branch topB_SD = CentralConfig::GetInstance()->GetTopBranch("SdCalibrator");
    //topB_SD.GetChild("latchBin").GetData(fSdLatchNanos);

    //LoadConfig( topB, "sdLatchNanos", fSdLatchNanos, 6400.);

    // Fill the fModuleDelay map from the configuration file
    if (fSetTimeStamps) {
      utl::Branch moduleDelays = config.GetChild("moduleDelays");
      for (utl::Branch b = moduleDelays.GetFirstChild(); b; b = b.GetNextSibling()) {
        utl::AttributeMap attibutes = b.GetAttributes();
        const int counterId = det::VManager::FindComponent<unsigned int>("counterId", attibutes);
        const int moduleId = det::VManager::FindComponent<unsigned int>("moduleId", attibutes);
        const int globalModuleId = counterId*100 + moduleId;
        const double delay = b.Get<double>();
        fModuleDelay[globalModuleId] = delay;
      }
    }

    LoadConfig(config, "syncCountingHisto", fSyncCountingHisto, false);
    return eSuccess;
  }


  VModule::ResultFlag
  MdMuonCounter::Run(evt::Event& event)
  {
    //INFO("Starting");

    if (!event.HasMEvent()) {
      WARNING("\n+++++++++\nEvent without MEvent: nothing to be done. Skipping to next event.\n++++++++++\n");
      return eContinueLoop;
    }

    mevt::MEvent& me = event.GetMEvent();
    const mdet::MDetector& theMdDetector = det::Detector::GetInstance().GetMDetector();

    std::ostringstream os;
    for (mevt::MEvent::CounterIterator ic = me.CountersBegin(), ec = me.CountersEnd(); ic != ec; ++ic) {

      mevt::Counter& counter = *ic;
      const int counterId = counter.GetId();

      if (counter.IsRejected()) {
         os.str("");
         os << "Skipping the rejected counter " << counterId;
         INFO(os);
         continue;
      }

      os.str("");
      os << "Processing counter " << counterId;
      INFO(os);

      const mdet::Counter& mdetCounter = theMdDetector.GetCounter(counterId);
      int stationId = mdetCounter.GetAssociatedTankId();

      // GPS time of the partner SD station
      const utl::TimeStamp sdTraceStartTime = GetSdTraceStartTime(event, stationId);
      utl::TimeStamp mdTraceStartTime;
      const bool hasSdTimeStamps = bool(sdTraceStartTime);
      //cout << "hasSdTimeStamps is " << std::boolalpha << hasSdTimeStamps << endl;

      for (mevt::Counter::ModuleIterator im = counter.ModulesBegin(), em = counter.ModulesEnd(); im != em; ++im) {

        DEBUGLOG("Starting module loop");
        mevt::Module& module = *im;
        const int moduleId = module.GetId();
        const int moduleGlobalId = counterId*100 + moduleId;

        os.str("");
        os << "Processing module " << moduleId;
        DEBUGLOG(os);

        const mdet::Module& mdetModule = mdetCounter.GetModule(moduleId);

        // Check if module was flagged as rejected (by a selector) or, in the case of real data, if T1 from SD was not matched
        // this should imply an empty (all "0s" ) module traces
        if (module.IsRejected()) {
           os.str("");
           os << "Module " << module.GetId() << " of counter " << counterId;
           os << " discarded because flagged as: " << module.GetRejectionReason();
           WARNING(os);
           continue;
        }

        const bool hasModuleDelay = fModuleDelay.count(moduleGlobalId);
        //cout << "hasModuleDelay is " << std::boolalpha << hasModuleDelay << endl;

        // Set timestamps only if the following three conditions are met
        // 1 - the configuration variable fSetTimeStamps is true
        // 2 - the partner sd station has a timestamp
        // 3 - the module is in the fModuleDelay list
        const bool setTimeStamps = fSetTimeStamps && hasSdTimeStamps && hasModuleDelay;
        //cout << "setTimeStamps is " << std::boolalpha << setTimeStamps << endl;

        // Get the start time of the traces
        if (setTimeStamps)
          mdTraceStartTime = sdTraceStartTime + GetTraceOffset(module, counterId);

        //INFO("After setting mdTraceStartTime");

        for (mevt::Module::ChannelIterator ch = module.ChannelsBegin(), ech = module.ChannelsEnd(); ch != ech; ++ch) {
          //DEBUGLOG("Starting channel loop");
          mevt::Channel& channel = *ch;

          if (setTimeStamps)
            channel.SetTraceStartTime(mdTraceStartTime);

          //DEBUGLOG("After setting the channel trace start time");

          const unsigned int nRecPatternMatches = FillChannelRecData(channel);
          //DEBUGLOG("After filling the channel reconstructed data");

          os.str("");
          os << "counter " << counterId << ": "
             << "module "  << module.GetId() << "  "
             << "channel " << channel.GetId() << "  "
             << "were detected "<< nRecPatternMatches <<  " muons ";
          DEBUGLOG(os);
          if ( nRecPatternMatches ) fLog( os.str(), 3 );

        }

        FillModuleRecData(module, mdetModule);
        DEBUGLOG("After filling the module reconstructed data");

        // Set the start time of the channelsOn histogram
        if (setTimeStamps) {
          const double samplingTime = mdetModule.IsSiPM() ?
            mdetModule.GetFrontEndSiPM().GetMeanSampleRatePeriod() :
            mdetModule.GetFrontEnd().GetMeanSampleRatePeriod();
          // adjusting the start time by half sample period. ChannelsOn is a histogram, not a trace!
          const utl::TimeStamp channelsOnStartTime = mdTraceStartTime - utl::TimeInterval(0.5*samplingTime);
          module.GetRecData().SetChannelsOnStartTime(channelsOnStartTime);
        }
        DEBUGLOG("After setting the module histogram start time");

      } // end module loop

      // Setting the counter T50

      // Skip not candidate and empty counters
      if (!counter.IsCandidate() || counter.IsEmpty())
        continue;

      // The time correction requires the SD geometrical reconstruction
      if (!event.HasRecShower() ||  !event.GetRecShower().HasSRecShower())
        continue;

      // Skip counter if the SD gps is unavailable
      if (!sdTraceStartTime)
        continue;

      const evt::ShowerSRecData& sRecShower = event.GetRecShower().GetSRecShower();
      const utl::Vector& rAxis = sRecShower.GetAxis();

      // Compute the t50 in nanoseconds from the SD trace start
      // Apply the md syncronization and correct modules times by position
      const double t50 = ComputeSignalT50(counter, rAxis);

      // Refer the t50 to the GPS timestamp
      const utl::TimeStamp t50GpsTime = sdTraceStartTime + utl::TimeInterval(t50);

      counter.SetSignalT50(t50GpsTime);

      // Set the offset between the t50 and the start time of the SD signal

      // Checking that the partner SD station exists
      if (!event.HasSEvent() || !event.GetSEvent().HasStation(stationId)) {
        os.str("");
        os << "SD station " << stationId << " not found.";
        WARNING(os);
        continue;
      }

      const sevt::Station& station = event.GetSEvent().GetStation(stationId);

      if (!station.HasRecData()) {
        os.str("");
        os << "Reconstructed data of SD station " << stationId << " not found.";
        WARNING(os);
        continue;
      }

      const utl::TimeStamp sdSignalStartTime = station.GetRecData().GetSignalStartTime();

      utl::TimeInterval t50ToSdStart = t50GpsTime - sdSignalStartTime;
      double t50ToSdStartNanos = t50ToSdStart.GetInterval() / utl::nanosecond;
      counter.SetT502SdStart(t50ToSdStartNanos);

      os.str("");
      os << "MD t50 - SD t0 = " << t50ToSdStartNanos << " ns";
      INFO(os);

    } // end counter loop

    return eSuccess;
  }


  VModule::ResultFlag
  MdMuonCounter::Finish()
  {
    return eSuccess;
  }


  // Private methods

  unsigned int
  MdMuonCounter::FillChannelRecData(mevt::Channel& channel)
  {
    if (!channel.HasTrace())
      return 0;

    const utl::TraceB trace = channel.GetTrace();

    // Set the trace start time
    //channel.SetTraceStartTime(mdTraceStartTime);

    if (channel.HasRecData()) {
      if (channel.GetRecData().GetNumberOfPatternMatchs()) {
        fLog("Clearing", 2)
          (channel.GetRecData().GetNumberOfPatternMatchs())
          ("preexisting muons for channel=")(channel.GetId())(':')
          ("in rec info.");
        channel.GetRecData().ClearPatternMatchBins();
      }
    } else
      channel.MakeRecData();

    mevt::ChannelRecData& rData = channel.GetRecData();
    // reconstructed muons
    const unsigned int nRecPatternMatches = (*fCounter)(trace, rData);

    return nRecPatternMatches;
  }


  void
  MdMuonCounter::FillModuleRecData(mevt::Module& module, const mdet::Module& mdetModule)
    //const
  {
    DEBUGLOG("Starting FillModuleRecData");

    std::ostringstream os;

    if (!module.HasRecData())
      module.MakeRecData();
    mevt::ModuleRecData& mRecData = module.GetRecData();
    DEBUGLOG("After making mRecData");

    // Get a vector with the bin numbers of the pattern matches
    const std::vector<unsigned int>& vPatternMatches = GetPatternMatchBins(module);
    DEBUGLOG("After getting vPatternMatches");

    // Fill the channelsOn histogram
    FillChannelsOn(module, vPatternMatches, mdetModule);
    DEBUGLOG("After filling channelsOn");

    // Fill the pattern match times vector
    const double samplingTime = mdetModule.IsSiPM() ?
      mdetModule.GetFrontEndSiPM().GetMeanSampleRatePeriod() :
      mdetModule.GetFrontEnd().GetMeanSampleRatePeriod();

    for (std::vector<unsigned int>::const_iterator it = vPatternMatches.begin();
         it != vPatternMatches.end(); ++it) {
      const double matchTime = (*it) * samplingTime;
      mRecData.AddPatternMatchTime(matchTime);
    }
    DEBUGLOG("After filling the pattern match times vector");

    const size_t segmentation = module.GetNumberOfActiveChannels();
    mRecData.SetSegmentation(segmentation);
    DEBUGLOG("After setting the segmentation");

    const mdet::Scintillator& scintillator = *(mdetModule.ScintillatorsBegin());
    const double scintillatorArea = scintillator.GetArea();
    const double moduleArea = scintillatorArea * segmentation;
    mRecData.SetActiveArea(moduleArea);
    DEBUGLOG("After setting the active area");

    //os.str("");
    //os << "\nchannelsOn ";
    //for (utl::TraceUI::ConstIterator it=mRecData.GetChannelsOn().Begin(); it!=mRecData.GetChannelsOn().End(); ++it)
    //  os << (*it) << ' ';
    //os << " (total: " << mRecData.GetNumberOfChannelsOn() << ')';
    //INFO(os);

    mRecData.SetWindowSize(fWindowSize);
    DEBUGLOG("After setting the window size");

    MdLDFFinderAG::ProfLike modLike(segmentation, mRecData.GetChannelsOn());

    // Set the number of estimated muons
    if (!mRecData.IsSaturated()) {    // muons can be estimated for unsaturated modules only

      const double nunmberOfMuons = EstimateNumberOfMuons(mRecData.GetChannelsOn(), segmentation);
      mRecData.SetNumberOfEstimatedMuons(nunmberOfMuons);
      DEBUGLOG("After setting the number of muons");
      // calculate errors with the profile likelihood

      const double errorHigh = modLike.GetErrorHigh();
      mRecData.SetNumberOfMuonsErrorHigh(errorHigh);
      const double errorLow = modLike.GetErrorLow();
      mRecData.SetNumberOfMuonsErrorLow(errorLow);

      os.str("");
      os << "module " << module.GetId() << ", "
         << std::fixed << std::setprecision(1)
         << "area = " << moduleArea << " m2, "
         << "bars = " << segmentation << ", "
         << "muons = " << nunmberOfMuons << " +" << errorHigh << " -" << errorLow;
      INFO(os);

    } else {
      os.str("");
      os << "module " << module.GetId() << " is saturated";
      INFO(os);
    }

    // Set the lower limit to the number of muons for both saturated and unsaturated modules
    const double lowLimit = modLike.GetLowLimit();
    mRecData.SetNumberOfMuonsLowLimit(lowLimit);

    DEBUGLOG("Finishing FillModuleRecData");
  }


  std::vector<unsigned int>
  MdMuonCounter::GetPatternMatchBins(mevt::Module& module)
    //const
  {
    // Get pattern match bins
    std::vector<unsigned int> vMatchBins;
    for (mevt::Module::ChannelIterator channel = module.ChannelsBegin(), ech = module.ChannelsEnd();
         channel != ech; ++channel) {

      if (channel->HasRecData() && !channel->IsMasked()) {  // use unmasked channels only
        for (mevt::ChannelRecData::PatternMatchBinIterator mut = channel->GetRecData().PatternMatchBinsBegin(),
             end = channel->GetRecData().PatternMatchBinsEnd(); mut != end; ++mut) {
          vMatchBins.push_back(*mut);

          fLog("Matched pattern in module", 50)(module.GetId())("- channel")(channel->GetId())
              ("bin")(*mut);
        }
      }

    }

    return vMatchBins;
  }


  void
  MdMuonCounter::FillChannelsOn(mevt::Module& module, const std::vector<unsigned int>& vPatternMatches, const mdet::Module& mdetModule)
    const
  {
    // Get the size and binning of the module (they are the same for all channels)

    const unsigned int nTraceSamples = mdetModule.IsSiPM() ?
      mdetModule.GetFrontEndSiPM().GetBufferLength() :
      mdetModule.GetFrontEnd().GetBufferLength();

    // Intermediate histogram to group pattern matches in time windows
    unsigned int histogramStartTime = 1;
    unsigned int histogramSize = nTraceSamples;

    if (fSyncCountingHisto && vPatternMatches.size() > 0) {
        histogramStartTime = *min_element(vPatternMatches.begin(), vPatternMatches.end());
        histogramSize = nTraceSamples-histogramStartTime+1;
    }

    TH1I hChannelsOn("h1", "Pattern matches", histogramSize, histogramStartTime, nTraceSamples);
    for (std::vector<unsigned int>::const_iterator it = vPatternMatches.begin(); it != vPatternMatches.end(); ++it) {
      const unsigned int bin = *it + 1;     // beware shift between vector and TH1 indexes
      hChannelsOn.Fill(bin);  // load the histogram
    }
    if (fWindowSize <= histogramSize)
        hChannelsOn.Rebin(fWindowSize);  // group bins according to the window size
    else
        hChannelsOn.Rebin(histogramSize);

    //Create the ChannelsOn histogram
    mevt::ModuleRecData& mRecData = module.GetRecData();
    if (!mRecData.HasChannelsOn()) {
      mRecData.MakeChannelsOn();
    }
    utl::TraceUI& channelsOnTrace = mRecData.GetChannelsOn();

    // Fill the ChannelsOn histogram
    const unsigned int nChannelsOnBins = hChannelsOn.GetNbinsX();
    for (unsigned int bin = 1; bin <= nChannelsOnBins; ++bin) {
      const unsigned int nChannelsOn =  hChannelsOn.GetBinContent(bin);
      channelsOnTrace.PushBack(nChannelsOn);
    }

    // Set the ChannelsOn binning
    const double samplingTime = mdetModule.IsSiPM() ?
      mdetModule.GetFrontEndSiPM().GetMeanSampleRatePeriod() :
      mdetModule.GetFrontEnd().GetMeanSampleRatePeriod();
    const double channelsOnBinning = samplingTime * fWindowSize;
    channelsOnTrace.SetBinning(channelsOnBinning);

    std::ostringstream info;
    info << "Counting channels on\n"
            "Number of bins = " << nChannelsOnBins << "\n"
            "Binning (ns) = " << channelsOnBinning;
    DEBUGLOG(info);
  }


  double
  MdMuonCounter::EstimateNumberOfMuons(const utl::TraceUI& channelsOn, const size_t segmentation)
    const
  {
    double estimatedMuons = 0;

    for (utl::TraceUI::ConstIterator it = channelsOn.Begin(); it != channelsOn.End(); ++it)
      if (*it)
        estimatedMuons -= segmentation * std::log(1 - double(*it)/segmentation);

    return estimatedMuons;
  }


  utl::TimeInterval
  MdMuonCounter::GetTraceOffset(const mevt::Module& module, const int counterId)
    const
  {
    const int moduleId = module.GetId();
    const int moduleGlobalId = counterId*100 + moduleId;
    const double moduleDelay = fModuleDelay.find(moduleGlobalId)->second;

    const sdet::SDetector& sDetector = det::Detector::GetInstance().GetSDetector();

    const mdet::MDetector& theMdDetector = det::Detector::GetInstance().GetMDetector();
    const mdet::Counter& mdetCounter = theMdDetector.GetCounter(counterId);
    const mdet::Module& mdetModule = mdetCounter.GetModule(moduleId);

    const unsigned int nLatchBin = mdetModule.IsSiPM() ?
      mdetModule.GetFrontEndSiPM().GetPreT1BufferLength() :
      mdetModule.GetFrontEnd().GetPreT1BufferLength();

    const double samplingTime = mdetModule.IsSiPM() ?
      mdetModule.GetFrontEndSiPM().GetMeanSampleRatePeriod() :
      mdetModule.GetFrontEnd().GetMeanSampleRatePeriod();

    const double traceLatchTime = (nLatchBin - 1) * samplingTime;
    const int stationId = mdetCounter.GetAssociatedTankId();
    double sdLatchNanos = sDetector.GetStation(stationId).GetLatchBin() * sDetector.GetStation(stationId).GetFADCBinSize();
    utl::TimeInterval traceOffset(moduleDelay + sdLatchNanos - traceLatchTime);

    return traceOffset;
  }


  utl::TimeStamp
  MdMuonCounter::GetSdTraceStartTime(const evt::Event& event, const int stationId)
    const
  {
    // Checking that the partner SD station exists
    if (!event.HasSEvent() || !event.GetSEvent().HasStation(stationId)) {
      std::ostringstream os;
      os << "Not setting the timestamps for the MD counter " << stationId << ". "
            "Cannot find the associated SD station.";
      WARNING(os);
      return utl::TimeStamp();
    }
    const sevt::Station& station = event.GetSEvent().GetStation(stationId);

    // Check that the station has GPS data
    if (!station.HasGPSData()) {
      std::ostringstream os;
      os << "Not setting the timestamps for the MD counter " << stationId << ". "
            "Cannot find the gps time of associated SD station.";
      WARNING(os);
      return utl::TimeStamp();
    }

    // Get the GPS time
    return utl::TimeStamp(station.GetTraceStartTime());
  }


  double
  MdMuonCounter::ComputeSignalT50(mevt::Counter& counter, const utl::Vector& rAxis)
    const
  {
    const int counterId = counter.GetId();

    const det::Detector& detector = det::Detector::GetInstance();

    const mdet::MDetector& theMdDetector = detector.GetMDetector();
    const mdet::Counter& mdetCounter = theMdDetector.GetCounter(counterId);

    const int stationId = mdetCounter.GetAssociatedTankId();

    const sdet::SDetector& sDetector = detector.GetSDetector();
    const utl::Point& sdPosition = sDetector.GetStation(stationId).GetPosition();

    std::vector<double> counterPatternMatches;

    for (mevt::Counter::ModuleIterator im = counter.ModulesBegin(), em = counter.ModulesEnd(); im != em; ++im) {

      DEBUGLOG("Starting module loop");
      mevt::Module& module = *im;
      const int moduleId = module.GetId();

      // Skip rejected and silent modules
      if (!module.IsCandidate())
        continue;

      const double moduleTraceOffset = GetTraceOffset(module, counterId);

      // Compute the correction to the muon times due to the position of each module relative to the SD (height included)
      const utl::Point& modulePosition = mdetCounter.GetModule(moduleId).GetPosition();
      const utl::Vector& dPosition =  modulePosition - sdPosition;
      const double positionCorrection = GetTimeCorrection(dPosition, rAxis);

      const std::vector<double> modulePatternMatches = module.GetRecData().GetPatternMatchTimes();

      for (std::vector<double>::const_iterator ip = modulePatternMatches.begin(), ep = modulePatternMatches.end(); ip != ep; ++ip) {
        const double correctedPatternMatchTime = *ip + moduleTraceOffset + positionCorrection;
        counterPatternMatches.push_back(correctedPatternMatchTime);
      }

    } // end module loop

    const double t50 = GetMedian(counterPatternMatches);

    return t50;
  }


  double
  MdMuonCounter::GetMedian(const std::vector<double> v)
    const
  {

    unsigned int ndata = v.size();
    assert(ndata>0);  // at least one datum needed to estimate the median

    // temporary vector to get rid of the ctness of vdata as required by sort
    std::vector<double> vtemp(v);
    std::sort(vtemp.begin(), vtemp.end());

    double median = 0;

    if (!(ndata % 2)) {   // if number of data is even
      const unsigned int i = ndata / 2;
      median = 0.5*(vtemp[i-1] + vtemp[i]); // median is the average of central values
    } else {  // if number of data is odd
      const unsigned int i = (ndata - 1) / 2;
      median = vtemp[i];  // usual median
    }

    return median;
  }


  double
  MdMuonCounter::GetTimeCorrection(const utl::Vector& dposition, const utl::Vector& axis)
    const
  {
    const double s = dposition * axis;
    const double timeCorrection = s / utl::kSpeedOfLight;;

    return timeCorrection;
  }

}