MdGeometryFitter.cc

Go to the documentation of this file.

#include "MdGeometryFitter.h"

#include <fwk/CentralConfig.h>
#include <fwk/RunController.h>
#include <fwk/LocalCoordinateSystem.h>

#include <mdet/MDetector.h>
#include <mdet/MTimeVariance.h>

#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>

#include <sevt/SEvent.h>
#include <sevt/EventTrigger.h>
#include <mevt/Counter.h>

using namespace utl;
using namespace std;

#define OUT(x) if ((x) <= fInfoLevel) cerr

namespace MdGeometryFitterAG {

  template<typename G>
  inline
  string
  ToString(const G& g, const CoordinateSystemPtr cs)
  {
    ostringstream os;
    os << '('
       << g.GetX(cs)/meter << ", "
       << g.GetY(cs)/meter << ", "
       << g.GetZ(cs)/meter << ") [m]";
    return os.str();
  }


  // perpendicular distance from the axis
  inline
  double
  RPerp(const TVector3& axis, const TVector3& station)
  {
    const double scal = axis.Dot(station);
    double r2 = station.Dot(station) - scal*scal;

    return sqrt(r2);
  }


  fwk::VModule::ResultFlag
  MdGeometryFitter::Init()
  {
    fwk::CentralConfig* cc = fwk::CentralConfig::GetInstance();

    const Branch topB = cc->GetTopBranch("MdGeometryFitter");

    topB.GetChild("infoLevel").GetData(fInfoLevel);
    if (fInfoLevel < eNone || fInfoLevel > eMinuit) {
      INFO("<infoLevel> out of bounds");
      return eFailure;
    }
    fMinuitOutput = (fInfoLevel >= eMinuit) ? true : false;

    topB.GetChild("minNumberOfStations").GetData(fMinNumberOfStations);
    if (fMinNumberOfStations <= 2) {
      INFO("<minNumberOfStations> must be > 2");
      return eFailure;
    }

    topB.GetChild("minNumberForFullCurvatureFit").GetData(fMinNumberForFullCurvatureFit);
    if (fMinNumberForFullCurvatureFit <= 3) {
      INFO("<fMinNumberForFullCurvatureFit> must be > 3");
      return eFailure;
    }

    topB.GetChild("CurvatureRadiusParameter1").GetData(fCurvatureRadiusParameter1);
    topB.GetChild("CurvatureRadiusParameter2").GetData(fCurvatureRadiusParameter2);

    return eSuccess;
  }

  fwk::VModule::ResultFlag
  MdGeometryFitter::Run(evt::Event& event)
  {
    INFO("MD geometrical reconstruction");

    // nothing to do if there is neither a SDEvent nor a MEvent
    if ( !event.HasMEvent() || !event.HasSEvent() || !event.GetSEvent().HasTrigger() ) {
      return eContinueLoop;
    }

    // the SD reconstruction is a prerequisite for the md plane fit  
    if ( !event.HasRecShower() || !event.GetRecShower().HasSRecShower() ) {
      return eContinueLoop;    
    }

    mevt::MEvent* mEvent = &event.GetMEvent();

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

    // the absolute points in time and space
    const CoordinateSystemPtr siteCS = detector.GetSiteCoordinateSystem();
    const utl::Point siteOrigin(0,0,0, siteCS);

    const evt::ShowerSRecData& sdRecShower = event.GetRecShower().GetSRecShower();
    const utl::Point sdCorePosition = sdRecShower.GetCorePosition();
    const utl::TimeStamp sdCoreTime = sdRecShower.GetCoreTime();
    const utl::CoordinateSystemPtr localCS = fwk::LocalCoordinateSystem::Create(sdCorePosition);

    TVector3 sdAxis;
    const double sdTheta = sdRecShower.GetAxis().GetTheta(localCS);
    const double sdPhi = sdRecShower.GetAxis().GetPhi(localCS);
    sdAxis.SetMagThetaPhi( 1., sdTheta, sdPhi); // taking the seed axis from the SD reconstruction
//    std::cout << " sdAxis = (" << sdAxis.x() << ", " << sdAxis.y() << ", " << sdAxis.z() << ")" << std::endl;  
    
    const double curvatureRadius = fCurvatureRadiusParameter1 + fCurvatureRadiusParameter2*sdTheta*sdTheta;

    // Create the geometry fitter
    GeometryFitter geometryFitter;
    geometryFitter.SetMinuitOutput(fMinuitOutput);

    TVector3 sdCorePositionLCS;
    sdCorePositionLCS.SetX( sdCorePosition.GetX(localCS) );
    sdCorePositionLCS.SetY( sdCorePosition.GetY(localCS) );
    sdCorePositionLCS.SetZ( sdCorePosition.GetZ(localCS) );
    geometryFitter.SetCore(sdCorePositionLCS);   

    geometryFitter.SetRadius(curvatureRadius);  

    SetTimeData(mEvent, sdCorePosition, sdCoreTime, sdAxis, geometryFitter);
    int nCounters = geometryFitter.GetNDetectors();

    // not worth it
    if (nCounters < fMinNumberOfStations) {
      OUT(eFinal)
        << "Not enough counters to reconstruct the geometry." << endl;
      return eContinueLoop;
    }

    geometryFitter.SetCoreFix(); // fix the core position

    if ( nCounters < fMinNumberForFullCurvatureFit ) {
      geometryFitter.SetCurvatureFix();  // fix the radius of curvature
    } else {
      geometryFitter.SetCurvatureFree(); // release the radius of curvature
    }

//    OUT(eObscure)
//      << "sd core postion = " << ToString(sdCorePosition, siteCS) << " (site)\n" 

//    OUT(eObscure)
//      << "local/site zenith angle diff. = " << acos(Vector(0,0,1, localCS)*Vector(0,0,1, siteCS))/degree << " [deg]" << endl;

    // Fit geometry
    if (geometryFitter.Fit() != EXIT_SUCCESS) {
      OUT(eIntermediate) << "Failed geometrical reconstruction.\n";
      return eContinueLoop;
    }

    // Fill the event
    FillRecShower(event, geometryFitter, sdCorePosition, sdCoreTime, sdAxis);

    FillCounter(mEvent, geometryFitter);

    OUT(eFinal)
      << "MD geometrical reconstruction \n"
         "  selected counters = " << nCounters << "\n" 
         "  axis  = (" << geometryFitter.GetU() << ", " << geometryFitter.GetV() << ", " << geometryFitter.GetW() << ") (SD core)\n"
         "  theta = " << geometryFitter.GetTheta()/degree << " +/- "
                      << geometryFitter.GetThetaError()/degree << " [deg] (SD core)\n"
         "  phi   = " << geometryFitter.GetPhi()/degree << " +/- "
                      << geometryFitter.GetPhiError()/degree << " [deg]\n"
         "  radius of curvature = " << geometryFitter.GetRadius() << " +/- " 
                                    << geometryFitter.GetRadiusError() << " [m]\n"
         "  ct0   = " << geometryFitter.GetCt0()/meter << " [m]\n"
         "  time chi2     = " << geometryFitter.GetChi2() << " / " << geometryFitter.GetNdof() << "\n"
         "  time residual = " << geometryFitter.GetMeanTimeResidual() << " +/- " << geometryFitter.GetTimeResidualSpread() << "\n"
         "  sd axis diff = " << acos( sdAxis.Dot(geometryFitter.GetAxis()) ) / degree  << " [deg]" << endl << endl;    

    ++fwk::RunController::GetInstance().GetRunData().GetNamedCounters()["MdGeometryFitter/Geometry"];

    return eSuccess;
  }

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

  void 
  MdGeometryFitter::SetTimeData(const mevt::MEvent* mEvent, utl::Point sdCorePosition, utl::TimeStamp sdCoreTime, const TVector3& axis, GeometryFitter& fitter) 
  {

    utl::CoordinateSystemPtr localCS = fwk::LocalCoordinateSystem::Create(sdCorePosition);
    const mdet::MDetector& mDetector = det::Detector::GetInstance().GetMDetector();
    const mdet::MTimeVariance& timeVariance = mdet::MTimeVariance::GetInstance();

    fitter.RemoveDetectors(); // remove previous data

    // Fill the fitter with time data
    for (mevt::MEvent::CounterConstIterator cIt = mEvent->CountersBegin(); cIt != mEvent->CountersEnd(); ++cIt) {

      if ( !cIt->IsCandidate() || !cIt->HasT50() || cIt->IsEmpty() ) {  // skip counters without t50
        continue;
      }

      const int counterId = cIt->GetId();

      const Vector cPos = mDetector.GetCounter(counterId).GetPosition() - sdCorePosition;
      const double x = cPos.GetX(localCS);
      const double y = cPos.GetY(localCS);
      const double z = cPos.GetZ(localCS);
      const TVector3 pos(x, y, z);
      const double time = (cIt->GetSignalT50() - sdCoreTime).GetInterval();

      const int nChannelsOn = cIt->GetNumberOfChannelsOn();
      const double distance = RPerp(axis, pos);   // using the SD axis to measure distances
      const double sigma = sqrt( timeVariance.GetTimeSigma2(nChannelsOn, distance) );

      fitter.AddDetector (counterId, pos, time, sigma);
    
    }

  }

  void 
  MdGeometryFitter::FillRecShower(evt::Event& event, const GeometryFitter& geometryFitter, 
                                  utl::Point sdCorePosition, utl::TimeStamp sdCoreTime, TVector3 sdAxis) {

    if (!event.HasRecShower()) {
      event.MakeRecShower();
    }

    if (!event.GetRecShower().HasMRecShower()) {
      event.GetRecShower().MakeMRecShower();
    }

    evt::ShowerMRecData& showerMRec = event.GetRecShower().GetMRecShower();

    const utl::CoordinateSystemPtr localCS = fwk::LocalCoordinateSystem::Create(sdCorePosition);

    // Fill the axis 
    const Vector axis(geometryFitter.GetU(), geometryFitter.GetV(), geometryFitter.GetW(), localCS);
    showerMRec.SetAxis(axis);

    showerMRec.SetThetaError( geometryFitter.GetThetaError() );
    showerMRec.SetPhiError( geometryFitter.GetPhiError() );
    showerMRec.SetThetaPhiCorrelation( geometryFitter.GetPhiError() );

    // Fill the core time
    const utl::TimeStamp coreTime = sdCoreTime + TimeInterval( geometryFitter.GetCt0()/kSpeedOfLight );
    const double coreTimeError = TimeInterval( geometryFitter.GetCt0Error()/kSpeedOfLight );
    showerMRec.SetCoreTime(coreTime, coreTimeError);

    // Fill the radius of curvature
    const double curvatureRadius = geometryFitter.GetRadius();
    const double curvatureRadiusError = geometryFitter.GetRadiusError(); 
    showerMRec.SetCurvature(curvatureRadius, curvatureRadiusError);

    // Fill the core position
    const TVector3 corePosition1 = geometryFitter.GetCore();
    const double xcore = corePosition1.x();
    const double ycore = corePosition1.y();
    const double zcore = corePosition1.z();  
    const utl::Point corePosition2( xcore, ycore, zcore, localCS); 
    showerMRec.SetCorePosition(corePosition2);

    // Fill the chi2
    const double chi2 = geometryFitter.GetChi2();
    const int angleNDOF = geometryFitter.GetNdof();
    showerMRec.SetAngleChi2(chi2, angleNDOF);

    // Angle between the md and sd axis
    const double mdSdAngle = acos( sdAxis.Dot(geometryFitter.GetAxis())); 
    showerMRec.SetMdSdAngle(mdSdAngle); 

    // Fill the reconstruction flags
    showerMRec.SetGeometryReconstructed(true);
    showerMRec.SetCurvatureFixed( geometryFitter.IsCurvatureFix() );
    showerMRec.SetCoreFixedGeo( geometryFitter.IsCoreFix() ); 
     
    // Fill the mean time residual
    showerMRec.SetTimeResidualMean( geometryFitter.GetMeanTimeResidual() );
    showerMRec.SetTimeResidualSpread( geometryFitter.GetTimeResidualSpread() );

  }


  void 
  MdGeometryFitter::FillCounter(mevt::MEvent* mEvent, const GeometryFitter& fitter) 
  {

    for (mevt::MEvent::CounterIterator cIt = mEvent->CountersBegin(); cIt != mEvent->CountersEnd(); ++cIt) {

      if ( !cIt->IsCandidate() || !cIt->HasT50() || cIt->IsEmpty() ) {  // skip counters without t50
        continue;
      }

      const int counterId = cIt->GetId();

      if (fitter.HasDetector(counterId)) {

        const double residual = fitter.GetTimeResidual(counterId);
        cIt->SetTimeResidual(residual);

        const double planeFrontDelay = fitter.GetPlaneFrontDelay(counterId);
        cIt->SetPlaneFrontDelay(planeFrontDelay);        

        const double t50error = fitter.GetTimeError(counterId);
        cIt->SetT50Error(t50error);

      } // end if

    } // end for loop

  } // end FillCounter

} // end namespace MdGeometryFitterAG

// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// End: