JMuonSimplex.hh

Go to the documentation of this file.

#ifndef JMUONSIMPLEX_INCLUDE
#define JMUONSIMPLEX_INCLUDE

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <algorithm>

#include "JDAQ/JDAQTimeslice.hh"
#include "JDAQ/JDAQSummaryslice.hh"

#include "JTrigger/JHit.hh"
#include "JTrigger/JTimeslice.hh"
#include "JTrigger/JSuperFrame2D.hh"
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JHitL1.hh"
#include "JTrigger/JHitR1.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JBuildL2.hh"

#include "JFit/JLine1Z.hh"
#include "JFit/JLine1ZEstimator.hh"
#include "JFit/JFitToolkit.hh"
#include "JFit/JEvt.hh"
#include "JFit/JEvtToolkit.hh"
#include "JFit/JFitParameters.hh"
#include "JFit/JFitStatus.hh"
#include "JFit/JLine3Z.hh"
#include "JFit/JModel.hh"
#include "JFit/JSimplex.hh"
#include "JFit/JLine3ZRegressor.hh"
#include "JFit/JMuonSimplexParameters_t.hh"

#include "JGeometry3D/JOmega3D.hh"
#include "JGeometry3D/JRotation3D.hh"

#include "JLang/JSharedPointer.hh"

/**
 * \author mdejong, gmaggi
 */

namespace JFIT
{
  /**
   * class to handle the Simplex angular reconstruction.
   * this should be combined with JMuonPrefit first, and subsequently with JMuonGandalf to achieve an optimal angular resolution.
   */
  class JMuonSimplex
  {
  private:
    JLANG::JSharedPointer<const JDETECTOR::JModuleRouter> router_;
    JFIT::JMuonSimplexParameters_t parameters_;

  public:

    /**
     * Default constructor
     */
    JMuonSimplex()
    {}    

    /**
     * Parameterized constructor
     *
     * \param router JSharedPointer of JModuleRouter, this is built via detector file.
     * \param parameters struct that holds default-optimized parameters for the reconstruction, for ARCA and ORCA configuration.
     */
    JMuonSimplex(const JLANG::JSharedPointer<const JDETECTOR::JModuleRouter> &router, 
                 const JFIT::JMuonSimplexParameters_t &parameters):
      router_(router),
      parameters_(parameters)
    {}
    
    ~JMuonSimplex()
    {}  
    /**
     * Declaration of Member function that actually performs the reconstruction
     *
     * \param timeSlice which is contructed via a JDAQEvent
     * \param InPreFits input JEvt which is a container of prefits, normally JMuonPrefit
     * \param OutFits non const JEvt.
     */
    void
    getJEvt(const KM3NETDAQ::JDAQTimeslice &timeSlice,
            JFIT::JEvt &InPreFits, // <- this parameter is declared as non const on purpose, see the algorthm below. 
            JFIT::JEvt &OutFits) const;
        
  };
}


/**
 * Member function definition
 */
void
JFIT::JMuonSimplex::getJEvt(const KM3NETDAQ::JDAQTimeslice &timeSlice,
                            JFIT::JEvt &InPreFits,
                            JFIT::JEvt &OutFits) const
{
  if ( InPreFits.empty() ) return;

  const bool useL0              =  parameters_.useL0;
  const bool reprocess          =  parameters_.reprocess; 
  const double sigma_ns         =  parameters_.sigma_ns;
  const size_t numberOfPrefits  =  parameters_.numberOfPrefits;
  const double Tmax_ns          =  parameters_.Tmax_ns;
  const double roadWidth_m      =  parameters_.roadWidth_m; 
  const double ctMin            =  parameters_.ctMin; 

  typedef std::vector<JTRIGGER::JHitL0>  JDataL0_t;
  typedef std::vector<JTRIGGER::JHitL1>  JDataL1_t;
  typedef std::vector<JTRIGGER::JHitR1>  JDataR1_t;

  JTRIGGER::JSuperFrame2D<JTRIGGER::JHit>        buffer;
  const JTRIGGER::JBuildL0<JTRIGGER::JHitL0>     buildL0;
  const JTRIGGER::JBuildL2<JTRIGGER::JHitL1>     buildL2(2, Tmax_ns, ctMin);

  typedef JFIT::JRegressor<JFIT::JLine3Z, JFIT::JSimplex>  JRegressor_t;

  JRegressor_t fit(sigma_ns);

  fit.estimator.reset(new JMEstimatorLorentzian());

  //JRegressor_t::debug              = debug; 
  JRegressor_t::MAXIMUM_ITERATIONS = 10000;

  const double TMAX_NS = 50.0; //[ns]

  JEvt::iterator __end = InPreFits.end(); 

  if (reprocess) {
    __end = std::partition(InPreFits.begin(), __end, JHistory::is_not_event(JFIT::JMUONSIMPLEX));
  }
  
  if (InPreFits.begin() != __end) {

    std::copy(InPreFits.begin(), __end, std::back_inserter(OutFits));

    if (numberOfPrefits > 0) {
      std::advance(__end = InPreFits.begin(), std::min(numberOfPrefits, OutFits.size()));
    }

    std::partial_sort(InPreFits.begin(), __end, InPreFits.end(), qualitySorter);

    __end = std::partition( InPreFits.begin(), __end, 
                            JFIT::JHistory::is_event(InPreFits.begin()->getHistory())
                            );

    JDataL0_t dataL0;
    JDataL1_t dataL1;

    for (KM3NETDAQ::JDAQTimeslice::const_iterator i = timeSlice.begin(); 
         i != timeSlice.end(); ++i) {

      if (router_->hasModule(i->getModuleID())) {
      
        buffer(*i, router_->getModule(i->getModuleID()));
        
        buildL0(buffer, std::back_inserter(dataL0));
        buildL2(buffer, std::back_inserter(dataL1));
      }
    }

    JDataR1_t data;

    for (JFIT::JEvt::const_iterator track = InPreFits.begin(); track != __end; ++track) {

      const JGEOMETRY3D::JRotation3D     R (JFIT::getDirection(*track));
      const JFIT::JLine1Z         tz(getPosition (*track).rotate(R), track->getT());
      const JFIT::JModel<JFIT::JLine1Z> match(tz, roadWidth_m, JTimeRange(-TMAX_NS, +TMAX_NS));

      data.clear();

      for (JDataL1_t::const_iterator i = dataL1.begin(); i != dataL1.end(); ++i) {

        JTRIGGER::JHitR1 hit(*i);

        hit.rotate(R);

        if (match(hit)) {
          data.push_back(hit);
        }
      }

      if (useL0) {

        data.reserve(data.size() + dataL0.size());

        JDataR1_t::iterator __end1 = data.end();

        for (JDataL0_t::iterator i = dataL0.begin(); i != dataL0.end(); ++i) {

          if ( std::find_if( data.begin(), __end1, 
                             std::bind2nd(std::equal_to<KM3NETDAQ::JDAQModuleIdentifier>(), i->getModuleID())) 
               == __end1) {

            JTRIGGER::JHitR1 hit(*i);

            hit.rotate(R);

            if (match(hit)) {
              data.push_back(hit);
            }
          }
        }
      }

      fit.step.resize(5);

      fit.step[0] = JFIT::JLine3Z( JFIT::JLine1Z( JGEOMETRY3D::JVector3D(0.5, 0.0, 0.0), 0.0) );
      fit.step[1] = JFIT::JLine3Z( JFIT::JLine1Z( JGEOMETRY3D::JVector3D(0.0, 0.5, 0.0), 0.0) );
      fit.step[2] = JFIT::JLine3Z( JFIT::JLine1Z( JGEOMETRY3D::JVector3D(0.0, 0.0, 0.0), 1.0) );
      fit.step[3] = JFIT::JLine3Z( JFIT::JLine1Z( JGEOMETRY3D::JVector3D(), 0.0), JGEOMETRY3D::JVersor3Z(0.005, 0.0) );
      fit.step[4] = JFIT::JLine3Z( JFIT::JLine1Z( JGEOMETRY3D::JVector3D(), 0.0), JGEOMETRY3D::JVersor3Z(0.0, 0.005) );

      const int NDF = JFIT::getCount(data.begin(), data.end()) - fit.step.size();

      if (NDF >= 0) {

        const double chi2 = fit( JFIT::JLine3Z(tz), data.begin(), data.end() );

        JGEOMETRY3D::JTrack3D track3D(fit.value);

        track3D.rotate_back(R);

        const double energy(0);
        OutFits.push_back( JFIT::getFit(JFIT::JHistory(track->getHistory()).add(JFIT::JFitApplication_t::JMUONSIMPLEX), track3D, 
                                    JFIT::getQuality(chi2, NDF), NDF,energy,JFIT::JFitStatus_t::OKAY) );
      }
    }

    // apply default sorter 

    std::sort(OutFits.begin(), OutFits.end(), qualitySorter);
  }
  
}


#endif