G/FdProfileConstrainedGeometryFit.cc

Go to the documentation of this file.

#include <evt/Event.h>
//#include <evt/ShowerRecData.h>
#include <evt/ShowerFRecData.h>

#include <fevt/Eye.h>
#include <fevt/Telescope.h>
#include <fevt/Pixel.h>
#include <fevt/EyeHeader.h>
#include <fevt/FEvent.h>
#include <fevt/EyeRecData.h>
#include <fevt/TelescopeRecData.h>

#include <det/Detector.h>

#include <fdet/Eye.h>
#include <fdet/Pixel.h>
#include <fdet/FDetector.h>

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

#include <utl/PhysicalConstants.h>

#include "FdProfileConstrainedGeometryFit.h"

#include "ChiZeroRegression.h"
#include "ProfileChi2.h"

using namespace std;
using namespace fevt;
using namespace utl;


namespace FdProfileConstrainedGeometryFitPG {

fwk::VModule::ResultFlag
FdProfileConstrainedGeometryFit::Init()
{
  fPCGFitter.Init();

  Branch topBranch = fwk::CentralConfig::GetInstance()->GetTopBranch("FdProfileConstrainedGeometryFitPG");
  topBranch.GetChild("eyeCut").GetData(fEyeCut);

  return eSuccess;
}


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


fwk::VModule::ResultFlag
FdProfileConstrainedGeometryFit::Run(evt::Event& event)
{
  if (!event.HasFEvent())
    return eSuccess;
  FEvent& fevent = event.GetFEvent();

  bool success = false;

  for (FEvent::EyeIterator eye = fevent.EyesBegin(ComponentSelector::eHasData),
       end = fevent.EyesEnd(ComponentSelector::eHasData); eye != end; ++eye) {

    if ((!eye->HasRecData() ||
        !eye->GetRecData().HasFRecShower()) &&
        !det::Detector::GetInstance().GetFDetector().GetEye(*eye).IsVirtual()
        /*eye->GetRecData().GetTimeFitNDof() <=0*/)
      continue;

    unsigned int eyeId = eye->GetId();
    cout << "PCGF for eye " << eyeId << " ...\n";
    unsigned int eyeBit = (fEyeCut / int(pow(10., eyeId-1))) % 10;
    if (eyeBit != 1) {
      cout << "skipping due to eyeCut" << "\n";
      continue;
    }

    if (det::Detector::GetInstance().GetFDetector().GetEye(*eye).IsVirtual()) {
      const int prevEye = GetDataFromPreviousFit(event, eyeId);
      if (prevEye == 0) {
        INFO("Could not find results from previous fits --- jump to next eye!");
        continue;
      }
    }

    if (fPCGFitter.Run(*eye)) {
      success = true;
    } else {
      eye->SetStatus(ComponentSelector::eDeSelected);
    }

    cout << "PCGF for eye " << eye->GetId() << " DONE\n";
  }

  if (success) {
    ++fwk::RunController::GetInstance().GetRunData().GetNamedCounters()["FdProfileConstrainedGeometryFitPG/Good"];
    return eSuccess;
  }
  
  return eContinueLoop;
}

int
FdProfileConstrainedGeometryFit::GetDataFromPreviousFit(evt::Event& event, int fgCurEye)
{

    const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
    fevt::FEvent& fevent = event.GetFEvent();

    set<unsigned int> eyeIds;

    // if the eye is virtual, take the data from any physical eye
    if (detFD.GetEye(fgCurEye).IsVirtual()) {
      for (fdet::FDetector::EyeIterator eyeIter = detFD.EyesBegin(fdet::FDetComponentSelector::ePhysical);
          eyeIter != detFD.EyesEnd(fdet::FDetComponentSelector::ePhysical);
          ++eyeIter) {

        if (fevent.HasEye(eyeIter->GetId())) {
          eyeIds.insert(eyeIter->GetId());
        }

      }
    }
    // if the eye is physical, take the data from this eye
    // or in case no previous data is found take it from any eye
    else {
      bool hasData = false;
      if (fevent.GetEye(fgCurEye).HasRecData()) {
        if ((fevent.GetEye(fgCurEye).GetRecData().HasFRecShower()) && (fevent.GetEye(fgCurEye).GetRecData().GetRp() != 0.)) {
          hasData = true;
        }
      }

      if (hasData) {
        eyeIds.insert(fgCurEye);
      }
      else {
        for (fdet::FDetector::EyeIterator eyeIter = detFD.EyesBegin(fdet::FDetComponentSelector::ePhysical);
            eyeIter != detFD.EyesEnd(fdet::FDetComponentSelector::ePhysical);
            ++eyeIter) {

          if (fevent.HasEye(eyeIter->GetId())) {
            eyeIds.insert(eyeIter->GetId());
          }

        }
      }
    }

    // Find the most time fit pixels
    int maxNFitPixels = 0;
    unsigned int maxNFitPixelsId = 0;
    for (set<unsigned int>::const_iterator IdIter = eyeIds.begin();
        IdIter != eyeIds.end();
        ++IdIter) {

      // Check if the previous modules found something
      if (!fevent.GetEye(*IdIter).HasRecData()) {
        continue;
      }
      else if (!fevent.GetEye(*IdIter).GetRecData().HasFRecShower()) {
        continue;
      }

      const fevt::EyeRecData& eyeRecData = fevent.GetEye(*IdIter).GetRecData();
      double nFitPixels = eyeRecData.GetNTimeFitPixels();

      if (nFitPixels > maxNFitPixels) {
        maxNFitPixels = nFitPixels;
        maxNFitPixelsId = *IdIter;
      }

      cout << "Eye " << *IdIter << " has " << nFitPixels << " pixels in the time fit" << endl;

    } // for eyeIds

    int takeFromEye = 0;
    if (maxNFitPixelsId != 0) {
      takeFromEye = maxNFitPixelsId;
    }
    else {
      return takeFromEye;
    }
    cout<<"take from "<<takeFromEye<<endl;

    const fdet::Eye& detEye = det::Detector::GetInstance().GetFDetector().GetEye(takeFromEye);
    fevt::Eye& eventEye = event.GetFEvent().GetEye(takeFromEye);

    const fevt::EyeRecData& eyeRecData = eventEye.GetRecData();

    // Get the data
//    const Vector& fSDP = eyeRecData.GetSDP();
//    const double fChiZero = eyeRecData.GetChiZero();
    const double fRp = eyeRecData.GetRp();
    double fTZero = eyeRecData.GetTZero();
    const Vector& fAxis = eyeRecData.GetFRecShower().GetAxis();
    const Point& fCore = eyeRecData.GetFRecShower().GetCorePosition();
    const double fTZeroErr = eyeRecData.GetTZeroError();
    const double fRpErr = eyeRecData.GetRpError();
    const double fChiZeroErr = eyeRecData.GetChiZeroError();
    const Point telPos = detEye.GetPosition();
    const CoordinateSystemPtr telCS = fwk::LocalCoordinateSystem::Create(telPos);
    const double sdpTheta = eyeRecData.GetSDP().GetTheta(telCS);
    const double sdpPhi = eyeRecData.GetSDP().GetPhi(telCS);
    const Vector sdp = Vector(1, sdpTheta, sdpPhi, telCS, Vector::kSpherical);
    const Point P_in_tel = telPos + Normalized(Cross(sdp, fAxis))*fRp;

    // Readjust the TZero according to the trigger offset of the different eyes
    const utl::TimeInterval offset = eventEye.GetHeader().GetTimeStamp() - fevent.GetEye(fgCurEye).GetHeader().GetTimeStamp();
    fTZero += offset.GetInterval();

    fevt::EyeRecData& curEyeRecData = fevent.GetEye(fgCurEye).GetRecData();
    if (!curEyeRecData.HasFRecShower()) {
      curEyeRecData.MakeFRecShower();
    }

    const utl::TimeStamp eyeStamp = fevent.GetEye(fgCurEye).GetHeader().GetTimeStamp();
    const Point eyePos = det::Detector::GetInstance().GetFDetector().GetEye(fgCurEye).GetPosition();

    bool fillEye = true;
    for (fdet::Eye::TelescopeIterator detTelIt = detFD.GetEye(fgCurEye).TelescopesBegin();
       detTelIt != detFD.GetEye(fgCurEye).TelescopesEnd(); ++detTelIt)
    {
      const int curtelId = detTelIt->GetId();
      const Point curtelPos = detTelIt->GetPosition();
      const CoordinateSystemPtr curtelCS = fwk::LocalCoordinateSystem::Create(curtelPos);
      const Point core_in_curtel = fCore-fAxis*fCore.GetZ(curtelCS)/cos(fAxis.GetTheta(curtelCS));

      const Vector curtel_to_core = core_in_curtel-curtelPos;
      const double chi0_in_curtel = Angle(curtel_to_core, fAxis);
      const double rp_in_curtel = curtel_to_core.GetMag()*sin(chi0_in_curtel);
      const utl::AxialVector SDP_in_curtel = Normalized(Cross(fAxis, curtel_to_core));
      const Point P_in_curtel = curtelPos + Normalized(Cross(SDP_in_curtel, fAxis))*rp_in_curtel;
      const Vector Pt_to_Pct = P_in_curtel - P_in_tel;
      const double t0_in_curtel = fTZero + ((Pt_to_Pct*fAxis > 0.) ? -1. : 1.)*Pt_to_Pct.GetMag()/utl::kSpeedOfLight;

      if (fevent.GetEye(fgCurEye).HasTelescope(curtelId)) {
        fevt::Telescope& tel = fevent.GetEye(fgCurEye).GetTelescope(curtelId);

        if (!tel.HasRecData())
          tel.MakeRecData();

        fevt::TelescopeRecData& telRecData = tel.GetRecData();
        telRecData.SetSDP(SDP_in_curtel);
        telRecData.SetChiZero(chi0_in_curtel, fChiZeroErr);
        telRecData.SetTZero(t0_in_curtel, fTZeroErr);
        telRecData.SetRp(rp_in_curtel, fRpErr);
        telRecData.SetTimeFitChiSquare(0,1);
        telRecData.SetTimeFitCorrelations(0,0,0);
      }

      //fill eye - choose tel with the same position as eye (usually, more of them are available)
      if (fillEye && (curtelPos - eyePos).GetMag() < 1e-3) {
        curEyeRecData.SetChiZero(chi0_in_curtel, fChiZeroErr);
        curEyeRecData.SetTZero(t0_in_curtel, fTZeroErr);
        curEyeRecData.SetRp(rp_in_curtel, fRpErr);
        curEyeRecData.SetSDP(SDP_in_curtel);
        curEyeRecData.SetTimeFitChiSquare(0,1);
        curEyeRecData.SetTimeFitCorrelations(0,0,0);
        curEyeRecData.SetSDPFitChiSquare(0,1);

        evt::ShowerFRecData& frecShower = curEyeRecData.GetFRecShower();
        frecShower.SetAxis(fAxis);
        frecShower.SetCorePosition(core_in_curtel);
        try {
          double tCore = t0_in_curtel;
          tCore -= (rp_in_curtel / tan(chi0_in_curtel)) / kSpeedOfLight;
          utl::TimeStamp timeAtCore = eyeStamp - TimeInterval(1000*100*ns) + TimeInterval(tCore);
          frecShower.SetCoreTime(timeAtCore, 0); //, rp_in_tel, chi0_in_tel, t0_in_tel);
        } catch (OutOfBoundException& e) {
          ERROR(e.GetMessage());
          return false;
        }
        fillEye = false;
      }

      //SDP pixels
      const unsigned int parentEyeId = detTelIt->GetParentPhysicalEyeId();
      const unsigned int parentTelId = detTelIt->GetParentPhysicalId();
      if (fevent.HasEye(parentEyeId) && fevent.GetEye(parentEyeId).HasRecData() && fevent.GetEye(fgCurEye).HasTelescope(curtelId)) {
        fevt::EyeRecData& parentEyeRecData = fevent.GetEye(parentEyeId).GetRecData();
        for (fevt::EyeRecData::PixelIterator pixelIter = parentEyeRecData.SDPPixelsBegin();
          pixelIter != parentEyeRecData.SDPPixelsEnd(); ++pixelIter) {

          if (pixelIter->GetTelescopeId() == parentTelId) {
            curEyeRecData.AddSDPPixel(fevent.GetEye(fgCurEye).GetTelescope(curtelId).GetPixel(pixelIter->GetId()));
          }
        }

        //time fit pixels
        for (fevt::EyeRecData::PixelIterator pixelIter = parentEyeRecData.TimeFitPixelsBegin();
          pixelIter != parentEyeRecData.TimeFitPixelsEnd(); ++pixelIter) {
      
          if (pixelIter->GetTelescopeId() == parentTelId) {
            curEyeRecData.AddTimeFitPixel(fevent.GetEye(fgCurEye).GetTelescope(curtelId).GetPixel(pixelIter->GetId()));
          }
        }
      }
    }


    for (fevt::EyeRecData::PixelIterator pixelIter = curEyeRecData.PulsedPixelsBegin();
      pixelIter != curEyeRecData.PulsedPixelsEnd(); ++pixelIter) {

      unsigned int telId = pixelIter->GetTelescopeId();
      if (fevent.GetEye(fgCurEye).HasTelescope(telId)) {
        fevt::TelescopeRecData& telRecData = fevent.GetEye(fgCurEye).GetTelescope(telId).GetRecData();
        const Point telPos = det::Detector::GetInstance().GetFDetector().GetEye(fgCurEye).GetTelescope(telId).GetPosition();
        const CoordinateSystemPtr telCS = fwk::LocalCoordinateSystem::Create(telPos);
        const double sdpTheta = telRecData.GetSDP().GetTheta(telCS);
        const double sdpPhi = telRecData.GetSDP().GetPhi(telCS);

        const Vector sdp = Vector(1, sdpTheta, sdpPhi, telCS, Vector::kSpherical);
        const Vector vertical(0,0,1, telCS);
        const Vector horizontalInSDP = Normalized(Cross(sdp, vertical));

        const Vector pixelDirection = det::Detector::GetInstance().GetFDetector().GetPixel(*pixelIter).GetDirection();
        const Vector pixelDirInSDP = pixelDirection - (sdp * pixelDirection)*sdp;
        const double pixelChi = Angle(pixelDirInSDP, horizontalInSDP);

        pixelIter->GetRecData().SetChi_i(pixelChi);
      }
    }

    return takeFromEye;

  } // int HybridGeometryFinder::GetDataFromPreviousFit()

/*
void
FdProfileConstrainedGeometryFit::FitTimingOfGivenGeometry(const Eye& eye)
  const
{
  const EyeRecData& eyeRec = eye.GetRecData();
  ChiZeroRegression chi0reg(eyeRec);

  cout << "Parameters from Axis Finder:\n"
          "Chi0: " << eyeRec.GetChiZero() << "\n"
          "Chi^2: " << eyeRec.GetTimeFitChiSquare() << "\t"
          "Rp: " << eyeRec.GetRp()/m << "\t"
          "T0: " << eyeRec.GetTZero()/ns << "\n";

   double chi2 = 0;
   double rp = 0;
   double t0 = 0;
   chi0reg(eyeRec.GetChiZero(), chi2, rp, t0);

   cout << "scanning Chi_0\n";

   for (int i = 0; i <= 180; ++i) {
     cout << i << " degrees\n";

     chi0reg(i*degree, chi2, rp, t0);
     cout << "Rp: " << rp/m << "\t"
             "T0: " << t0/ns << "\t"
             "Chi^2_t: " << chi2 << "\t";
  }
}
*/
}