MdMuonIntegrator.cc

Go to the documentation of this file.

#include <utl/ErrorLogger.h>
#include <fwk/CentralConfig.h>

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>
#include <mevt/MEvent.h>
#include <mevt/Counter.h>
#include <utl/ConfigUtils.h>
#include <fwk/LocalCoordinateSystem.h>
#include <utl/CoordinateSystem.h>
#include <utl/Branch.h>

#include "MdMuonIntegrator.h"

#include <sstream>
#include <string>
#include <iomanip>

using namespace fwk;
using namespace utl;


namespace MdMuonIntegratorAG {

  VModule::ResultFlag
  MdMuonIntegrator::Init()
  {
    utl::Branch config = fwk::CentralConfig::GetInstance()->GetTopBranch("MdMuonIntegrator");

    LoadConfig(config, "baseLineLowLimit", fBaseLineLowLimit, 10);
    LoadConfig(config, "baseLineHighLimit", fBaseLineHighLimit, 210);

    LoadConfig(config, "lowThreshold", fLowThreshold, 2);
    LoadConfig(config, "highThreshold", fHighThreshold, 8);

    LoadConfig(config, "showerMaxWidth", fShowerMaxWidth, 200);
    LoadConfig(config, "delayBinaryADC", fDelayBinaryADC, 24);

    return eSuccess;
  }


  VModule::ResultFlag
  MdMuonIntegrator::Run(evt::Event& event)
  {
    if (!event.HasMEvent()) {
      WARNING("\n+++++++++\nEvent without MEvent: nothing to be done. Skipping to next event.\n++++++++++\n");
      return eContinueLoop;
    }

    INFO("Starting");

    mevt::MEvent& me = event.GetMEvent();

    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);

      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();

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

        if (module.IsRejected()) {
          os.str("");
          os << "Module " << module.GetId() << " of counter " << counterId;
          os << " discarded because flagged as: " << module.GetRejectionReason();
          WARNING(os);
          continue;
        }

        GetMuonsWithADC(event, module);

      }

    }

    return eSuccess;
  }


  // Private methods

  void
  MdMuonIntegrator::GetMuonsWithADC(const evt::Event& event, mevt::Module& module)
  {
    mevt::ModuleRecData& mRecData = module.GetRecData();

    if (!mRecData.IsEmpty()) {
      // Look for the first muon in the module (needed to optimize the integration window
      int leadingMuonBin = 0;
      int lastMuonBin = 0;
      GetModuleFirstMuon(leadingMuonBin, lastMuonBin, module);

      double cosZenith = 1;
      if (event.HasRecShower() && event.GetRecShower().HasSRecShower()) {
        const evt::ShowerSRecData& sRecShower = event.GetRecShower().GetSRecShower();
        const utl::Point& rCore = sRecShower.GetCorePosition();
        const utl::Vector& rAxis = sRecShower.GetAxis();
        const utl::CoordinateSystemPtr coreCS = LocalCoordinateSystem::Create(rCore);
        cosZenith = rAxis.GetCosTheta(coreCS);
      }

      if (mRecData.IsADCCalibratedLG() && module.HasIntegratorATrace()) {
        // Get ADC signal charge
        double charge = 0;
        const unsigned int windowSize = mRecData.GetWindowSize();
        const utl::TraceUSI& adcTraceLG =  module.GetIntegratorATrace();

        GetSignalCharge(charge, leadingMuonBin, lastMuonBin, windowSize, adcTraceLG);

        double muons = charge * cosZenith / mRecData.GetMeanChargeMuonLG();
        double resolution = mRecData.GetStdDevChargeMuonLG()/mRecData.GetMeanChargeMuonLG();

        // Get uncertainties
        double errorLow = 0;
        double errorHigh = 0;
        if (muons > 0)
          GetUncertainties(errorLow, errorHigh, muons, resolution);

        // Set results in event
        mRecData.SetNumberOfEstimatedMuonsADCLG(muons);
        mRecData.SetNumberOfMuonsADCErrorHighLG(errorLow);
        mRecData.SetNumberOfMuonsADCErrorLowLG(errorHigh);

        mRecData.SetTotalChargeLG(charge);

        std::ostringstream os;
        os << "Number of muons ADC LG = " << muons << ", "
              "Error Low  = " << errorLow << ", Error high = " << errorHigh << " "
              "Calibration " << mRecData.GetMeanChargeMuonLG() << " Charge " << charge;
        INFO(os);
      }

      if (mRecData.IsADCCalibratedHG() && module.HasIntegratorBTrace()) {
        // Get ADC signal charge
        double charge = 0;
        const unsigned int windowSize = mRecData.GetWindowSize();
        const utl::TraceUSI& adcTraceHG =  module.GetIntegratorBTrace();

        GetSignalCharge(charge, leadingMuonBin, lastMuonBin, windowSize, adcTraceHG);

        double muons = charge * cosZenith / mRecData.GetMeanChargeMuonHG();
        double resolution = mRecData.GetStdDevChargeMuonHG()/mRecData.GetMeanChargeMuonHG();

        // Get uncertainties
        double errorLow = 0;
        double errorHigh = 0;
        if (muons > 0)
          GetUncertainties(errorLow, errorHigh, muons, resolution);

        // Set results in event
        mRecData.SetNumberOfEstimatedMuonsADCHG(muons);
        mRecData.SetNumberOfMuonsADCErrorHighHG(errorLow);
        mRecData.SetNumberOfMuonsADCErrorLowHG(errorHigh);

        mRecData.SetTotalChargeHG(charge);

        std::ostringstream os;
        os << "Number of muons ADC HG = " << muons << ", "
              "Error Low  = " << errorLow << ", Error high = " << errorHigh << " "
              "Calibration " << mRecData.GetMeanChargeMuonHG() << " Charge " << charge;
        INFO(os);
      }
    }
  }


  void
  MdMuonIntegrator::GetModuleFirstMuon(int& leadingMuonBin, int& lastMuonBin, mevt::Module& module)
  {
    const utl::TraceUI& channelsOn = module.GetRecData().GetChannelsOn();

    for (unsigned int i = 0; i < channelsOn.GetSize(); ++i) {
      if (channelsOn.At(i) > 0) {
        if (leadingMuonBin == 0) {
          leadingMuonBin = i;
        }
        lastMuonBin = i;
      }
    }
  }


  void
  MdMuonIntegrator::GetSignalCharge(double& charge,
                                    const double leadingMuonBin, const double lastMuonBin,
                                    const unsigned int windowSize, const utl::TraceUSI& adcTrace)
  {
    unsigned int minBin = 0;
    unsigned int maxBin = 0;
    int start = leadingMuonBin * windowSize / 2;
    int end = (lastMuonBin + 1) * windowSize / 2;

    // Calculate base line for th HG trace and the LG trace
    double base = 0;
    double baseStdDev = 0;
    int j = 0;

    for (j = fBaseLineLowLimit; j < fBaseLineHighLimit; ++j)
      base += adcTrace.At(j);

    base /= fBaseLineHighLimit - fBaseLineLowLimit;

    for (j = fBaseLineLowLimit; j < fBaseLineHighLimit; ++j)
      baseStdDev += pow(adcTrace.At(j) - base, 2);

    baseStdDev /= fBaseLineHighLimit - fBaseLineLowLimit;
    baseStdDev = sqrt(baseStdDev);

    if (leadingMuonBin > 0) {
      unsigned int i = 0;
      // We set a low and high threshold levels based on the baseline fluctuations
      double thrhigh = baseStdDev * fHighThreshold;
      double thrlow = baseStdDev * fLowThreshold;
      unsigned int thrBin = 0;

      // This looks for the first bin with signal over the high threshold
      // after the leading muon time
      for (i = start; i < adcTrace.GetSize(); ++i) {
        if (thrBin == 0 && adcTrace.At(i) - base >= thrhigh) {
          thrBin = i;
          break;
        }
      }

      // This looks backwards till the signal falls below the low threshold

      minBin = thrBin;

      for (int i = thrBin; i > start+10; --i) {
        if (adcTrace.At(i) - base > thrlow) {
          minBin = i;
        } else if (adcTrace.At(i-1) - base < thrlow && adcTrace.At(i-2) - base < thrlow)
          break;
      }

      maxBin = minBin + fShowerMaxWidth;
      if (maxBin > adcTrace.GetSize() - 1)
        maxBin = adcTrace.GetSize() - 10; // avoid out of range

      // Now we integrate between the first bin with signal and the last bin with signal above the low threshold
      for (i = minBin; i <= maxBin; ++i) {
        if (i > thrBin) {
          // I exit the loop when the signal goes below the baseline, but I want to make sure it does not go up some bins later
          bool exitLoop = true;
          for (unsigned int j = i+1; j < i+6; ++j) {
            if (adcTrace.At(j) - base > thrlow) {
              exitLoop = false;
            } else {
              for (unsigned int k = j+1; k < maxBin; ++k) {
                if (adcTrace.At(k) - base > thrhigh) {
                  exitLoop = false;
                }
              }
            }
          }
          if (exitLoop)
            break;
        }
        charge += adcTrace.At(i) - base;
      }

      if (thrBin == 0) {
        // did not reach threshold, integrate noise in the region of the muon
        minBin = start + fDelayBinaryADC;
        maxBin = end + windowSize / 2. + fDelayBinaryADC;

        if (maxBin > adcTrace.GetSize())
          maxBin = adcTrace.GetSize() - 10; // avoid out of range
        for (i = minBin; i <= maxBin; ++i) {
          if (adcTrace.At(i) - base) {
            charge += adcTrace.At(i) - base;
          }
        }
      }
    }
  }


  void
  MdMuonIntegrator::GetUncertainties(double& errorLow, double& errorHigh, const double nMu, const double resolution)
  {
    // This is a fast implementation to obtain the poisson uncertainty.
    // This is not an implementation for production...
    double s = 1; // 1 sigma interval

    double xmin = (nMu - sqrt(nMu)) / 2;
    double xmax = nMu * 2;

    double stepLow = (nMu - xmin) / 100;
    double stepHigh = (xmax - nMu) / 100;

    double poissonWeight = 0;
    double pmin = 0;
    double pmax = 0;
    int i = 0;

    while (poissonWeight < s) {
      pmin = nMu - i*stepLow;
      poissonWeight = 2*((pmin - nMu) + nMu * log(nMu / pmin));
      ++i;
    }
    poissonWeight = 0;
    i = 0;
    while (poissonWeight < s) {
      pmax = nMu + i * stepHigh;
      poissonWeight = 2*((pmax - nMu) + nMu * log(nMu / pmax));
      ++i;
    }

    errorLow = sqrt(pow(nMu - pmin, 2) + nMu * pow(resolution, 2));
    errorHigh = sqrt(pow(pmax - nMu, 2) + nMu * pow(resolution, 2));
  }

}