MdCornerClippingCorrecter.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 <utl/CoordinateSystemPtr.h>

#include <fwk/LocalCoordinateSystem.h>

#include "MdCornerClippingCorrecter.h"

#include <sstream>

using namespace fwk;
using namespace utl;
using namespace MdCornerClippingCorrecterAG;

namespace MdCornerClippingCorrecterAG {
  MdCornerClippingCorrecter::MdCornerClippingCorrecter()
    { }

  MdCornerClippingCorrecter::~MdCornerClippingCorrecter()
    { }

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

    std::ostringstream os;
    os << "\ncomputeCornerClippingProbability = " << fComputeCornerClippingProbability << '\n'
       << "par0 = " << fPar0 << '\n'
       << "par1 = " << fPar1 << '\n'
       << "par2 = " << fPar2 << '\n'
       << "par3 = " << fPar3;
    INFO(os.str());

    return eSuccess;
  }

  fwk::VModule::ResultFlag
  MdCornerClippingCorrecter::Run(evt::Event& event)
  {
    if (!fComputeCornerClippingProbability) {
      INFO("No corner clipping probability will be computed.");
      return eSuccess;
    }

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

    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 corner clipping probability for the rejected counter " << counterId << ".";
        INFO(os);
        continue;
      }

      if (counter.IsSilent()) {
        os.str("");
        os << "No corner clipping probability 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;
        INFO(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;
        }

        FillCornerClippingProbability(event, module, mdm);

      }
    }
    return eSuccess;
  }

  fwk::VModule::ResultFlag
  MdCornerClippingCorrecter::Finish()
  {
    return eSuccess;
  }

  // private methods
  void
  MdCornerClippingCorrecter::FillCornerClippingProbability(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 cosTheta = rAxis.GetCosTheta(coreLocalCS);
    const double phi = rAxis.GetPhi(coreLocalCS) / utl::rad;

    // Orientation of module
    const double phi0 = mdetModule.GetPhi0() / utl::rad;

    // pcc computation
    const double sindphi = std::abs(std::sin(phi-phi0));
    const double m = fPar0 + fPar1 / cosTheta;
    const double b = fPar2 + fPar3 / cosTheta;
    const double pcc = b + m*sindphi;

    // Filling mRecData
    mevt::ModuleRecData& mRecData = module.GetRecData();
    mRecData.SetCornerClippingProbability(pcc);

    std::ostringstream os;
    os << "phi = " << phi << " rad, "
          "phi = " << phi / utl::deg << " deg, "
          "phi0 = " << phi0 << " rad, "
          "phi0 = " << phi0 / utl::deg << " deg, "
          "sin(phi-phi0)  = " << sindphi << ", "
          "cosTheta = " << cosTheta  << ", "
          "m = " << m << ", "
          "b = " << b << ", "
          "pcc = " << pcc;
    INFO(os);
  }

}