MdBiasCorrecter.cc

Go to the documentation of this file.

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

#include <det/Detector.h>
#include <mdet/MDetector.h>
#include <mdet/Counter.h>
#include <mdet/Module.h>

#include <evt/Event.h>
#include <mevt/MEvent.h>
#include <mevt/Counter.h>
#include <sevt/SEvent.h>

#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>
#include <evt/ShowerMRecData.h>
#include <mevt/Scintillator.h>

#include <utl/Branch.h>
#include <utl/ConfigUtils.h>

#include <fwk/LocalCoordinateSystem.h>

#include "MdBiasCorrecter.h"

#include <sstream>

using namespace fwk;
using namespace utl;
using namespace MdBiasCorrecterAG;


namespace MdBiasCorrecterAG {

  VModule::ResultFlag
  MdBiasCorrecter::Init()
  {
    // Configuration reading
    utl::Branch config = fwk::CentralConfig::GetInstance()->GetTopBranch("MdBiasCorrecter");
    LoadConfig(config, "biasCorrection", fBiasCorrection, false);
    LoadConfig(config, "par0", fPar0, 0);
    LoadConfig(config, "par1", fPar1, 0);
    LoadConfig(config, "par2", fPar2, 0);
    LoadConfig(config, "par3", fPar3, 0);
    LoadConfig(config, "par4", fPar4, 0);

    return eSuccess;
  }


  VModule::ResultFlag
  MdBiasCorrecter::Run(evt::Event& event)
  {
    if (!fBiasCorrection) {
      INFO("No bias correction applied.");
      return eSuccess;
    }

    if (!event.HasMEvent()) {
      WARNING("Event without MEvent: nothing to be done. Skipping to next event.");
      return eSuccess;
    }

    INFO("Starting");

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

    std::ostringstream os;
    for (auto cIt = me.CountersBegin(), cEnd = me.CountersEnd(); cIt != cEnd; ++cIt) {

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

      if (counter.IsRejected()) {
        os.str("");
        os << "No bias correction for the rejected counter " << counterId << ".";
        INFO(os);
        continue;
      }

      if (counter.IsSilent()) {
        os.str("");
        os << "No bias correction for the silent counter " << counterId << ".";
        INFO(os);
        continue;
      }

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

      const mdet::Counter& mdc = md.GetCounter(counterId);

      for (auto mIt = counter.ModulesBegin(), mEnd = counter.ModulesEnd(); mIt != mEnd; ++mIt) {

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

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

        const mdet::Module& mdm = mdc.GetModule(moduleId);

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

        CorrectBias(event, module, mdm);

      }

    }

    return eSuccess;
  }


  // Private methods

  double
  MdBiasCorrecter::CorrectionFactor(const double phi, const double theta)
  {
   const double a = fPar0 + fPar1*(1 + fPar2*std::cos(theta)) * std::sin(theta);
   const double b = fPar3*(1 + fPar4*std::cos(theta)) * std::sin(theta);
   return 1 / (1 + a + b*fabs(std::sin(phi)));
  }


  void
  MdBiasCorrecter::CorrectBias(const evt::Event& event, mevt::Module& module, const mdet::Module& mdetModule)
  {
    // Shower geometry
    const evt::ShowerSRecData& sRecShower = event.GetRecShower().GetSRecShower();
    const utl::Point& rCore = sRecShower.GetCorePosition();
    const utl::Vector& rAxis = sRecShower.GetAxis();
    const utl::CoordinateSystemPtr& coreLocalCS = LocalCoordinateSystem::Create(rCore);
    const double theta = rAxis.GetTheta(coreLocalCS);
    const double phi = rAxis.GetPhi(coreLocalCS);
    
    // Orientation of module
    const double phi0 = mdetModule.GetPhi0();

    // Apply correction to numberOfMuons and errors
    mevt::ModuleRecData& mRecData = module.GetRecData();
    const double numberOfMuons = mRecData.GetNumberOfEstimatedMuons();
    const double corrfac = CorrectionFactor(phi-phi0, theta);
    const double numberOfMuonsCor = numberOfMuons * corrfac;
    mRecData.SetNumberOfEstimatedMuons(numberOfMuonsCor);
    // currently no error estimate for bias corrected number of muons
    mRecData.SetNumberOfMuonsErrorHigh(0); 
    mRecData.SetNumberOfMuonsErrorLow(0);

    std::ostringstream os;
    os << "Uncorrected number of muons = " << numberOfMuons << ", "
          "corrected  = " << numberOfMuonsCor << ", "
          "correction factor = " << corrfac << ", "
          "for theta = " << theta << ", "
          "phi = " << phi << ", "
          "phi0 = " << phi0;
    INFO(os);
  }

}