RdScintPlaneFit.cc

Go to the documentation of this file.

#include "RdScintPlaneFit.h"

#include <cmath>

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

#include <det/Detector.h>
#include <rdet/RDetector.h>

#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/config.h>
#include <utl/PhysicalConstants.h>
#include <utl/AxialVector.h>
#include <utl/Math.h>
#include <utl/ErrorLogger.h>

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerRRecData.h>
#include <evt/ShowerRRecDataQuantities.h>

#include <revt/REvent.h>
#include <revt/Station.h>
#include <revt/Header.h>
#include <revt/StationRecData.h>
#include <revt/StationRRecDataQuantities.h>

#include <CLHEP/Matrix/Matrix.h>
#include <CLHEP/Matrix/Vector.h>

#include <TMinuit.h>

using namespace std;
using namespace revt;
using namespace evt;
using namespace fwk;
using namespace utl;

#define OUT(x) if ((x) <= fInfoLevel) cerr << "  "
#define NL "\n  "

using CLHEP::HepMatrix;
using CLHEP::HepVector;

namespace RdScintPlaneFit {
  // global variables
  REvent* RdScintPlaneFit::fgCurrentREvent;
  Point RdScintPlaneFit::fgBarycenter;
  TimeStamp RdScintPlaneFit::fgBaryTime;
  CoordinateSystemPtr RdScintPlaneFit::fgLocalCS;
  Point RdScintPlaneFit::fgCoordinateOrigin;
  double RdScintPlaneFit::fgMeanEventTime(0);

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

  VModule::ResultFlag RdScintPlaneFit::Init()
  {
    // Initialize your module here. This method
    // is called once at the beginning of the run.
    // The eSuccess flag indicates the method ended
    // successfully.  For other possible return types,
    // see the VModule documentation.

    INFO("RdScintPlaneFit::Init()");

    // Read in the configurations of the xml file
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdScintPlaneFit");

    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);

    fMinuitOutput = (fInfoLevel >= eMinuit) ? true : false;

    topBranch.GetChild("minNumberOfStations").GetData(fMinNumberOfStations);
    topBranch.GetChild("allowUnphysicalCosines").GetData(fAllowUnphysicalCosines);
    
    return eSuccess;
  }

// analytic linear fit in local CS; approximation by neglecting
// z-coordinate, since z ~ 0 in local bary coordinate system,
// therefore loosing the w component of the axis. it is funny
// though, that stations living on a plane can not distinguish
// between upwards or downwards-going showers.
  VModule::ResultFlag RdScintPlaneFit::LinearFit(FitParameters& fit, const bool reuseFit) const
  {
    
    const Vector approx_axis = (reuseFit ? Vector(fit.u, fit.v, fit.w, fgLocalCS) : Vector());
    
    const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();
    
    double s1(0);
    double sx(0);
    double sy(0);
    double st(0);
    double sxx(0);
    double sxy(0);
    double syy(0);
    double sxt(0);
    double syt(0);
    //double stt(0);
    
    //loop over all stations
    
    int loop(0);
    
    for(REvent::StationIterator sIt = fgCurrentREvent->StationsBegin(); sIt != fgCurrentREvent->StationsEnd(); ++sIt) {
        
      
      const Vector sPos = rDetector.GetStation(*sIt).GetPosition() - fgBarycenter;
      
      const double x = sPos.GetX(fgLocalCS);
      
      const double y = sPos.GetY(fgLocalCS);
      
      const StationRecData& sRec = sIt->GetRecData();
      
      if(sRec.HasParameter(eAERAScintStationSignalTime)) {
        const double t = sRec.GetParameter(eAERAScintStationSignalTime) - fgMeanEventTime;
        const double invSigma2 = 1. / (5 * 5);
        s1 += invSigma2;
        const double x_s = x * invSigma2;
        sx += x_s;
        sxx += x * x_s;
        sxy += y * x_s;
        sxt += t * x_s;
        const double y_s = y * invSigma2;
        sy += y_s;
        syy += y * y_s;
        syt += t * y_s;
        st += t * invSigma2;
      }
      ++loop;
    }
     
    st *= kSpeedOfLight;
    sxt *= kSpeedOfLight;
    syt *= kSpeedOfLight;

    HepMatrix mat(3, 3);

    mat[0][0] = sxx;
    mat[0][1] = sxy;
    mat[0][2] = -sx;
    mat[1][0] = sxy;
    mat[1][1] = syy;
    mat[1][2] = -sy;
    mat[2][0] = -sx;
    mat[2][1] = -sy;
    mat[2][2] = s1;

    int ierr(0);
    mat.invert(ierr);
    
    if(ierr) {
      INFO("Matrix inversion failed.");
      return eFailure;
    }

    HepVector k(3);
     
    k[0] = -sxt;
    k[1] = -syt;
    k[2] = st;

    const HepVector p = mat * k;
    double w2 = 1 - Sqr(p[0]) - Sqr(p[1]);
    
    if(w2 < 0) {
      WARNING("Unphysical direction cosines");
      if(fAllowUnphysicalCosines) {
        WARNING("Setting z-compenent of the shower axis to 0");
        w2 = 0;
      } else {
        return eFailure;
      }
    }
     
    fit.u = p[0];
    fit.v = p[1];
    fit.w = sqrt(w2);
    fit.ct0 = p[2];

    fit.sigmaU2 = mat[0][0];
    fit.sigmaV2 = mat[1][1];
    fit.sigmaUV = mat[0][1];
    fit.sigmact02 = mat[2][2];
     
    
    
    return eSuccess;
  }

  VModule::ResultFlag RdScintPlaneFit::Run(evt::Event& event)
  {
    INFO("Plane fit");

    // nothing to do if there is no REvent
    if(!event.HasREvent()) {
      OUT(eIntermediate)
      << "eContinueLoop: No REvent, so no reconstruction"
                        << endl;
      return eContinueLoop;
    }

    fgCurrentREvent = &event.GetREvent();
    //const REvent& rEvent = *fgCurrentREvent;

    const det::Detector& detector = det::Detector::GetInstance();
    //const rdet::RDetector& rDetector = detector.GetRDetector();
    const CoordinateSystemPtr siteCS = detector.GetSiteCoordinateSystem();

    evt::ShowerRecData& shower = event.GetRecShower();
    evt::ShowerRRecData& rShower = shower.GetRRecShower();
     
    int nStations = rShower.GetParameter(eNumberOfAERAScintStationWithPulseFound);
     
    // not worth it
    if(nStations < fMinNumberOfStations) {
      OUT(eFinal) << "eContinueLoop: Not enough stations."
                  << endl;
                 
     return eContinueLoop;
    }
     
    //ShowerRRecData& showerrrec = event.GetRecShower().GetRRecShower();
    //const TimeStamp& eventTime = rEvent.GetHeader().GetTime();

    //weighted sums
    TimeInterval timeSum(0);
    double weightSum(0);
    for(REvent::StationIterator sIt = fgCurrentREvent->StationsBegin(); sIt != fgCurrentREvent->StationsEnd(); ++sIt) {

      //const rdet::Station& dStation = rDetector.GetStation(*sIt);
      const StationRecData& sRec = sIt->GetRecData();
      if(sRec.HasParameter(eAERAScintStationSignalTime)) {

        const double weight = sqrt(sRec.GetParameter(eAERAScintStationSignalHeight));

        weightSum += weight;
        timeSum += weight * sRec.GetParameter(eAERAScintStationSignalTime);

      }
    }  //station loop
         
    timeSum /= weightSum;

    fgMeanEventTime= rShower.GetParameter(eBaryTimeAERAScint);
    fgBaryTime = rShower.GetParameter(eBaryTimeAERAScint) + timeSum;
    
    const CoordinateSystemPtr referenceCS = detector.GetReferenceCoordinateSystem();
    
    const Point bary(rShower.GetParameter(eBaryCenterAERAScintX),rShower.GetParameter(eBaryCenterAERAScintY),rShower.GetParameter(eBaryCenterAERAScintZ),referenceCS);
    
    fgBarycenter = bary;
    fgLocalCS = LocalCoordinateSystem::Create(fgBarycenter);

    OUT(eIntermediate) << "# candidate stations = "
                       << nStations << endl;

  
    if(!rShower.HasParameter(eBaryCenterAERAScintX) || !rShower.HasParameter(eBaryTimeAERAScint)){
      OUT(eFinal) << "eContinueLoop: Not Bary Data available!."
                  << endl;
        
      return eContinueLoop;
    }
    
    //site coordinate system, only you used for terminal output
    OUT(eIntermediate) << "barycenter = "
                       << ToString(fgBarycenter, siteCS) << " (site)" << NL
                       << "bary time = "
                       << timeSum / nanosecond << " [ns] (to event time)" << endl;

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

    // this will point to the final result
    FitParameters finalFitResult;
   
    // first approximation
    FitParameters fitLin;
  
     if(LinearFit(fitLin) != eSuccess) {
      OUT(eIntermediate)
      << "eContinueLoop: Linear fit was not succesfull"
                        << endl;
       return eContinueLoop;
     }
    
    // try again with previous values as approxiamtions for sigma
    LinearFit(fitLin, true);
   
    // with first approximation on u and v, the nonlinear fit is attempted
    FitParameters fitNonlin = fitLin;
    bool useNonLinearFit(false);
   
    if(PlaneFit3DDriver(fitNonlin) == eSuccess) {
      OUT(eIntermediate) << "Stage "
                         << fitLin.stage << ": linear plane fit" << NL
                         << "  axis = ("
                         << fitLin.u << ", " << fitLin.v << ", " << fitLin.w << ") (bary)" << NL
                         << "  theta = "
                         << acos(fitLin.w) / degree << NL
                         << "  phi   = "
                         << atan2(fitLin.v, fitLin.u) / degree << NL
                         << "  ct0 = "
                         << fitLin.ct0 / meter << " [m]" << endl;
      useNonLinearFit = true;
      finalFitResult = fitNonlin;
    
    } else {
      OUT(eIntermediate)
      << "  Nonlinear fit failed, using linear approximation."
                        << endl;
      finalFitResult = fitLin;

    }
    // Fill the event
    const Vector axis(finalFitResult.u, finalFitResult.v, finalFitResult.w, fgLocalCS);

    rShower.SetParameter(eAERAScintShowerAxisX, axis.GetX(fgLocalCS));
    rShower.SetParameter(eAERAScintShowerAxisY, axis.GetY(fgLocalCS));
    rShower.SetParameter(eAERAScintShowerAxisZ, axis.GetZ(fgLocalCS));

    double chi2(0);
    //calculates residuals and fill them in station rec
    CalculateTimeResidual(axis, finalFitResult.ct0, chi2);

    const int angleNDOF = nStations - 3;
    rShower.SetParameter(eAERAScintPlaneFitNDF, angleNDOF);
    rShower.SetParameter(eAERAScintPlaneFitChi2, chi2);
    
    OUT(eFinal) << "Stage "
                << finalFitResult.stage << ": " << (useNonLinearFit ? "3d" : "linear") << " plane fit" << NL
                << "  axis  = ("
                << finalFitResult.u << ", " << finalFitResult.v << ", " << finalFitResult.w << ") (bary)" << NL
                << "  theta = "
                << acos(finalFitResult.w) / degree << " +/- " << NL
                << "  phi   = " << atan2(finalFitResult.v, finalFitResult.u) / degree << " +/- " << NL
                << "  ct0   = " << finalFitResult.ct0 / meter << " [m]" << NL
                << endl;

    return eSuccess;
  }

  VModule::ResultFlag RdScintPlaneFit::PlaneFit3DDriver(FitParameters& fit) const
  {
    TMinuit minuit(3);
    int errFlag(0);
    double argList[10];
    argList[0] = -1;
    if(!fMinuitOutput) {
      minuit.mnexcm("SET PRINTOUT", argList, 1, errFlag);
      minuit.mnexcm("SET NOWARNINGS", argList, 0, errFlag);
      minuit.SetPrintLevel(-1);
    }

    if(fAllowUnphysicalCosines) {
      minuit.SetFCN(RdScintPlaneFit::PlaneFit3DHorizonFnc);
      argList[0] = 1; // chi2
      minuit.mnexcm("SET ERRORDEF", argList, 1, errFlag);
      minuit.mnparm(0, "theta", acos(fit.w), 0.01, 0, 0, errFlag);
      minuit.mnparm(1, "phi", atan2(fit.v, fit.u), 0.01, 0, 0, errFlag);
      minuit.mnparm(2, "ct0", fit.ct0, 10 * meter, 0, 0, errFlag);
    } else {
      minuit.SetFCN(RdScintPlaneFit::PlaneFit3DFnc);
      argList[0] = 1; // chi2
      minuit.mnexcm("SET ERRORDEF", argList, 1, errFlag);
      minuit.mnparm(0, "u", fit.u, 0.01, 0, 0, errFlag);
      minuit.mnparm(1, "v", fit.v, 0.01, 0, 0, errFlag);
      minuit.mnparm(2, "ct0", fit.ct0, 10 * meter, 0, 0, errFlag);
    }

    // init fnc
    argList[0] = 1;
    minuit.mnexcm("CALI", argList, 1, errFlag);

    argList[0] = 500;
    minuit.mnexcm("MINIMIZE", argList, 1, errFlag);
    if(errFlag) {
      OUT(eObscure)
      << "  minuit minimize error: "
                   << errFlag << '\n';
      minuit.mnexcm("MINOS", argList, 0, errFlag);
      if(errFlag) {
        OUT(eObscure)
        << "  minuit minos error: "
                     << errFlag << '\n';
        return eFailure;
      }
    }

    // Get results
    double foo(0);
    double par1(0);
    double par2(0);

    minuit.GetParameter(0, par1, foo);
    minuit.GetParameter(1, par2, foo);

    if(fAllowUnphysicalCosines) {
      fit.u = sin(par1) * cos(par2);
      fit.v = sin(par1) * sin(par2);
      fit.w = cos(par1);
    } else {
      fit.u = par1;
      fit.v = par2;
      double w2 = 1 - Sqr(fit.u) - Sqr(fit.v);
      if(w2 < 0) {
        OUT(eIntermediate) << " eFailure: Fit failed due to unphysical angle cosines."
                           << endl;
        return eFailure;
      }
      fit.w = sqrt(w2);
    }

    minuit.GetParameter(2, fit.ct0, foo);

    // error matrix
    double emat[3][3];
    minuit.mnemat(&emat[0][0], 3);
    fit.sigmaU2 = emat[0][0];
    fit.sigmaV2 = emat[1][1];
    fit.sigmaUV = emat[0][1];
    fit.sigmact02 = emat[2][2];

    return eSuccess;
  }

  VModule::ResultFlag RdScintPlaneFit::Finish()
  {
    // Put any termination or cleanup code here.
    // This method is called once at the end of the run.
    INFO("RdPlaneFit::Finish()");
    return eSuccess;
  }

  void RdScintPlaneFit::PlaneFit3DFnc(int& /*nPar*/, double* const /*grad*/, double& value, double* const par, const int flag)
  {
    const double& u = par[0];
    const double& v = par[1];
    const double& ct0 = par[2];
    static int nStations = 0;
    static vector<double> sTiming;     //station timing
    static vector<double> sSignal;     //station signal
    static vector<Vector> sPosition;   //station position
    static vector<double> sTimeSigma2; //Uncurtainty on the pulse timing

    if(flag == 1) {  // fnc init

      sTiming.clear();
      sSignal.clear();
      sPosition.clear();
      sTimeSigma2.clear();

      const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();

      for(REvent::StationIterator sIt = fgCurrentREvent->StationsBegin(); sIt != fgCurrentREvent->StationsEnd(); ++sIt) { //loop over stations

        const StationRecData& sRec = sIt->GetRecData();
        
        if(!sRec.HasParameter(eAERAScintStationSignalTime)){
            continue;
        }
   
          sTiming.push_back(kSpeedOfLight * (sRec.GetParameter(eAERAScintStationSignalTime) - fgMeanEventTime));
          sSignal.push_back(sRec.GetParameter(eAERAScintStationSignalHeight));
          sPosition.push_back(rDetector.GetStation(*sIt).GetPosition() - fgBarycenter);
          sTimeSigma2.push_back(sRec.GetParameterError(eAERAScintStationSignalTime) * sRec.GetParameterError(eAERAScintStationSignalTime));

      }

      nStations = sTiming.size();

    } // fnc init

    const double w2 = 1 - Sqr(u) - Sqr(v);

    if(w2 < 0) {
      value = 1e30;
      return;
    }

    const double w = sqrt(w2);

    const Vector axis(u, v, w, fgLocalCS);

    double chi2(0);

    for(int i = 0; i < nStations; ++i)
      chi2 += Sqr(sTiming[i] - ct0 + axis * sPosition[i]) / sTimeSigma2[i];

    value = chi2 / kSpeedOfLight2;
  }

  void RdScintPlaneFit::PlaneFit3DHorizonFnc(int& /*nPar*/, double* const /*grad*/, double& value, double* const par, const int flag)
  {
    const double& theta = par[0];
    const double& phi = par[1];
    const double& ct0 = par[2];
    static int nStations = 0;
    static vector<double> sTiming;     //station timing
    static vector<double> sSignal;     //station signal
    static vector<Vector> sPosition;   //station position
    static vector<double> sTimeSigma2; //Uncurtainty on the pulse timing
//  static const sdet::STimeVariance* timeVariance = 0;

    if(flag == 1) {  // fnc init

      sTiming.clear();
      sSignal.clear();
      sPosition.clear();
      sTimeSigma2.clear();

      const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();

      for(REvent::StationIterator sIt = fgCurrentREvent->StationsBegin(); sIt != fgCurrentREvent->StationsEnd(); ++sIt) { //loop over stations

        const StationRecData& sRec = sIt->GetRecData();
        if(sRec.HasParameter(eAERAScintStationSignalHeight)) {
          sTiming.push_back(kSpeedOfLight * (sRec.GetParameter(eAERAScintStationSignalTime) - fgMeanEventTime));
          sSignal.push_back(sRec.GetParameter(eAERAScintStationSignalHeight));
          sPosition.push_back(rDetector.GetStation(*sIt).GetPosition() - fgBarycenter);
          sTimeSigma2.push_back((sRec.GetParameterError(eAERAScintStationSignalTime) * sRec.GetParameterError(eAERAScintStationSignalTime)));
        }
      }

      nStations = sTiming.size();

    } // fnc init

    const Vector axis(cos(phi) * sin(theta), sin(phi) * sin(theta), cos(theta), fgLocalCS);

    double chi2(0);

    for(int i = 0; i < nStations; ++i)
      chi2 += Sqr(sTiming[i] - ct0 + axis * sPosition[i]) / sTimeSigma2[i];

    value = chi2 / kSpeedOfLight2;
  }

  void 
  RdScintPlaneFit::CalculateTimeResidual(const utl::Vector& axis, const double ct0, double& chi2) const
  {
      
    const double t0 = ct0 / kSpeedOfLight;
    const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();
    chi2 = 0;
    for(REvent::StationIterator sIt = fgCurrentREvent->StationsBegin(); sIt != fgCurrentREvent->StationsEnd(); ++sIt) { //loop over stations
        StationRecData& sRec = sIt->GetRecData();
        if(sRec.HasParameter(eAERAScintStationSignalHeight)) {
          const Vector stLoc_c = (rDetector.GetStation(*sIt).GetPosition() - fgBarycenter) / kSpeedOfLight;
          const double t = (sRec.GetParameter(eAERAScintStationSignalTime) - fgMeanEventTime);
          double sigma_t =sRec.GetParameterError(eAERAScintStationSignalTime);
          double timeResidual = t - t0 + (stLoc_c * axis); 
          sRec.SetParameter(eAERAScintStationTimeResidual, timeResidual);
          sRec.SetParameterError(eAERAScintStationTimeResidual, sigma_t);
          chi2 += pow(timeResidual / sigma_t,2);

        }
      }
  }

}