JMuonPrefit.hh

Go to the documentation of this file.

#ifndef JMUONPREFIT_INCLUDE
#define JMUONPREFIT_INCLUDE
 
#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <algorithm>
 
#include "JDAQ/JDAQTimeslice.hh"
#include "JDAQ/JDAQSummaryslice.hh"
 
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JHitL1.hh"
#include "JTrigger/JHitR1.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JBuildL2.hh"
#include "JTrigger/JAlgorithm.hh"
#include "JTrigger/JMatch1D.hh"
#include "JTrigger/JMatch3B.hh"
#include "JTrigger/JBind2nd.hh"
 
#include "JFit/JLine1Z.hh"
#include "JFit/JLine1ZEstimator.hh"
#include "JFit/JMatrixNZ.hh"
#include "JFit/JVectorNZ.hh"
#include "JFit/JFitToolkit.hh"
#include "JFit/JEvt.hh"
#include "JFit/JEvtToolkit.hh"
#include "JFit/JFitParameters.hh"
#include "JFit/JFitStatus.hh"
#include "JFit/JMuonPrefitParameters_t.hh"
 
#include "JTools/JConstants.hh"
#include "JTools/JPermutation.hh"
 
#include "JGeometry3D/JOmega3D.hh"
#include "JGeometry3D/JRotation3D.hh"
 
#include "JLang/JPredicate.hh"
#include "JLang/JSharedPointer.hh"
 
/**
 * \author mdejong, gmaggi
 */
 
namespace JFIT
{
  /**
   * class to handle Prefit angular reconstruction. 
   * this should be combined with JMuonSimplex and JMuonGandalf to achieve an optimal angular resolution.
   */
  class JMuonPrefit
  {
  private:
    JLANG::JSharedPointer<const JDETECTOR::JModuleRouter> router_;
    JFIT::JMuonPrefitParameters_t parameters_;
    JGEOMETRY3D::JOmega3D omega_;    
 
    
  public:
 
    /**
     * Default constructor
     */
    JMuonPrefit()
    {}
 
    /**
     * 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.
     */
    JMuonPrefit(const JLANG::JSharedPointer<const JDETECTOR::JModuleRouter> &router,
                const JFIT::JMuonPrefitParameters_t &parameters):
      router_(router),
      parameters_(parameters),
      omega_(parameters.gridAngle_deg* JTOOLS::PI/180.0)
    {}
 
 
    /**
     * 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.
     * \param omega which is constructed by gridAngle_deg (part of JMuonPrefitParameters_t)
     */
    JMuonPrefit(const JLANG::JSharedPointer<const JDETECTOR::JModuleRouter> &router,
                const JFIT::JMuonPrefitParameters_t &parameters,
                const JGEOMETRY3D::JOmega3D &omega): 
      router_(router),
      parameters_(parameters),
      omega_(omega)
    {}
       
    ~JMuonPrefit()
    {}  
 
    /**
     * Declaration Member function that actually performs the reconstruction 
     *
     * \param timeSliceBuildL0 which is contructed via a JDAQEvent
     * \param timeSliceBuildL2 which is contructed via a JDAQEvent
     * \param out non const JEvt.
     */
    void
    getJEvt(const KM3NETDAQ::JDAQTimeslice &timeSliceBuildL0,
            const KM3NETDAQ::JDAQTimeslice &timeSliceBuildL2,
            JFIT::JEvt &out) const;
    
    /**
     * Sort hits according times corrected for position along z-axis.
     * used as a predicate in std::partial_sort
     * 
     * \param  first             first  hit 
     * \param  second            second hit
     * \return                   true if first hit earlier than second hit; else false
     */
    static bool 
    cmz(const JTRIGGER::JHitR1& first, 
        const JTRIGGER::JHitR1& second)
    {
      return (first .getT() * JTOOLS::getSpeedOfLight()  -  first .getZ()  <
              second.getT() * JTOOLS::getSpeedOfLight()  -  second.getZ());
    }
    
  };
}
 
 
/**
 * Member function definition
 */
 
void
JFIT::JMuonPrefit::getJEvt(const KM3NETDAQ::JDAQTimeslice &timeSliceBuildL0,
                           const KM3NETDAQ::JDAQTimeslice &timeSliceBuildL2,
                           JFIT::JEvt &out) const
{
  const double sigma_ns        = parameters_.sigma_ns;
  const double gridAngle_deg   = parameters_.gridAngle_deg;
  const bool useL0             = parameters_.useL0;
  const int numberOfOutliers   = parameters_.numberOfOutliers;
  const std::size_t numberOfPrefits = parameters_.numberOfPrefits;
  const double Tmax_ns         = parameters_.Tmax_ns;
  const double ctMin           = parameters_.ctMin;
  const double roadWidth_m     = parameters_.roadWidth_m;
 
  const int    FACTORY_LIMIT          =   8;                                // [unit]                           
  const double STANDARD_DEVIATIONS    =   3.0;                              // [unit]                           
  const double HIT_OFF                =   1.0e3 * sigma_ns * sigma_ns;      // [ns^2]                           
  const int    NUMBER_OF_L1HITS       =  (useL0 ? 1 : 4);
  const int    MAXIMUM_NUMBER_OF_HITS =  50; //numeric_limits<int>::max();                
 
  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(JTRIGGER::JL2Parameters(2, Tmax_ns, ctMin));
  const JTRIGGER::JMatch1D<JTRIGGER::JHitR1>     match1D(roadWidth_m, Tmax_ns);
  const JTRIGGER::JMatch3B<JTRIGGER::JHitL1>     match3B(roadWidth_m, Tmax_ns);
  const JFIT::JHitL1Comparator     compare;
 
  // L1 and L0 data
  JDataL0_t dataL0;
  JDataL1_t dataL1;
  
  buildL0( timeSliceBuildL0, *router_, back_inserter(dataL0) );
  buildL2( timeSliceBuildL2, *router_, back_inserter(dataL1) );
  
  // 3D cluster of unique optical modules
  
  JDataL1_t::iterator __end = std::unique(dataL1.begin(), dataL1.end(), 
                                          std::equal_to<KM3NETDAQ::JDAQModuleIdentifier>()
                                          );
 
  __end = JTRIGGER::clusterizeWeight(dataL1.begin(), __end, match3B);
 
  if (std::distance(dataL1.begin(), __end) >= NUMBER_OF_L1HITS) {
        
    JDataR1_t data;
        
    for (JGEOMETRY3D::JOmega3D_t::const_iterator dir = omega_.begin(); dir != omega_.end(); ++dir) {
      
      const JGEOMETRY3D::JRotation3D R(*dir);
      
      data.clear();
      
      std::copy(dataL1.begin(), __end, back_inserter(data));
      
      for (JDataR1_t::iterator i = data.begin(); i != data.end(); ++i) {
        i->rotate(R);
      }
      
      
      JDataR1_t::iterator __end1 = data.end();
      
      
      // reduce data
      
      if ( std::distance(data.begin(), __end1) > MAXIMUM_NUMBER_OF_HITS ) {
        
        std::advance(__end1 = data.begin(), MAXIMUM_NUMBER_OF_HITS);
        
        std::partial_sort(data.begin(), __end1, data.end(), cmz);
      }
      
      
      // 1D cluster
      
      __end1 = JTRIGGER::clusterizeWeight(data.begin(), __end1, match1D);
      
      if ( std::distance(data.begin(), __end1) < NUMBER_OF_L1HITS ) {
        continue;
      }
      
      if (useL0) {
        
        data.erase(__end1, data.end());
        
        data.reserve(data.size() + dataL0.size());
        
        __end1 = data.end();
        
        for (JDataL0_t::const_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 ){
            JHitR1 hit(*i);
            
            hit.rotate(R);
            
            if ( std::count_if(data.begin(), __end1, JBind2nd(match1D,hit)) >= NUMBER_OF_L1HITS ) {
              data.push_back(hit);
            }
          }
        }
        
        __end1 = JTRIGGER::clusterize(__end1, data.end(), match1D);
      }
      
            
      if ( std::distance(data.begin(), __end1) <= JFIT::JEstimator<JLine1Z>::NUMBER_OF_PARAMETERS ) {
        continue;
      }
     
      // 1D fit
           
      JFIT::JLine1Z  tz;
      double   chi2 = std::numeric_limits<double>::max();
      int      NDF  = 0;
      int      N    = 0;
      
      JFIT::JMatrixNZ V;
      JFIT::JVectorNZ Y;
      
      if ( std::distance(data.begin(), __end1) <= FACTORY_LIMIT ) {
        
        double ymin = std::numeric_limits<double>::max();
        
        JDataR1_t::iterator __end2 = __end1;
        
        for (int n = 0; 
             n <= numberOfOutliers && std::distance(data.begin(), __end2) > JFIT::JEstimator<JLine1Z>::NUMBER_OF_PARAMETERS; 
             ++n, --__end2) {
          
          std::sort(data.begin(), __end1, compare);
        
          do {
        
            try {
              
              JFIT::JEstimator<JFIT::JLine1Z> fit(data.begin(), __end2);
              
              V.set(fit, data.begin(), __end2, gridAngle_deg, sigma_ns);
              Y.set(fit, data.begin(), __end2);
              
              V.invert();
              
              double y = getChi2(Y, V);
              
              if        (y <= 0.0) {
                //std::cout<<"THIS IS A WARNING chi2(1) " << y <<std::endl;
                //WARNING(endl << "chi2(1) " << y << endl); this gives: text+0x4320): undefined reference to `debug'
              }else if (y < ymin) {
                ymin = y;
                tz   = fit;
                NDF  = std::distance(data.begin(), __end2) - JFIT::JEstimator<JLine1Z>::NUMBER_OF_PARAMETERS;
                N    = JFIT::getCount(data.begin(), __end2);
                chi2 = y;
              }
            }
            catch(const JLANG::JException& error) {}
        
          } while ( JTOOLS::next_permutation(data.begin(), __end2, __end1, compare) );
          
          ymin -= STANDARD_DEVIATIONS * STANDARD_DEVIATIONS;
        }
        
      }else {
        
        try {
          
          const int number_of_outliers = std::distance(data.begin(), __end1) - JFIT::JEstimator<JLine1Z>::NUMBER_OF_PARAMETERS - 1;
          
          JFIT::JEstimator<JLine1Z> fit(data.begin(), __end1);
          
          V.set(fit, data.begin(), __end1, gridAngle_deg, sigma_ns);
          Y.set(fit, data.begin(), __end1);
          
          V.invert();
          
          NDF  = std::distance(data.begin(), __end1) - JFIT::JEstimator<JFIT::JLine1Z>::NUMBER_OF_PARAMETERS;
          N    = JFIT::getCount(data.begin(), __end1);
          
          for (int n = 0; n <= number_of_outliers && NDF > 0; ++n, --NDF) {
            
            double ymax =  0.0;
            int    kill = -1;
            
            for (unsigned int i = 0; i != Y.size(); ++i) {
              
              double y = JFIT::getChi2(Y, V, i);
              
              if        (y <= 0.0) {
                //std::cout<<"THIS IS A WARNING chi2(1) " << y <<std::endl;
                //WARNING(endl << "chi2(2) " << y << endl);
              } else if (y > ymax) {
                ymax = y;
                kill = i;
              }
            }
            
            if (ymax < STANDARD_DEVIATIONS * STANDARD_DEVIATIONS) {
              break;
            }
            
            V.update(kill, HIT_OFF);
            
            fit.update(data.begin(), __end1, V);
            
            Y.set(fit, data.begin(), __end1);
            
            N -= JFIT::getCount(data[kill]);
          }
          
          chi2 = JFIT::getChi2(Y, V);
          tz   = fit;
        }
        catch(const JLANG::JException& error) {}
      }
             
      if (NDF != 0) {
        
        tz.rotate_back(R);
        
        const double energy(0);
        out.push_back( JFIT::getFit(JHistory(JMUONPREFIT), tz, *dir, getQuality(chi2, (useL0 ? N : NDF)), NDF,
                                    energy,JFIT::JFitStatus_t::OKAY) );
      }
     
    }
      
    // apply default sorter
    
    if (numberOfPrefits > 0) {
      
      JEvt::iterator __end = out.end();
      
      std::advance(__end = out.begin(), std::min(numberOfPrefits, out.size()));
      
      std::partial_sort(out.begin(), __end, out.end(), qualitySorter);
      
      // add downward pointing solutions if available but absent so far
      
      if ( std::count_if(out.begin(), __end, 
                         JLANG::make_predicate(&JFit::getDZ, 0.0, JLANG::JComparison::le()) 
                         ) == 0) {
        
        JEvt::iterator __begin = __end;
        JEvt::iterator __q     = std::partition(__begin, out.end(), 
                                                JLANG::make_predicate(&JFit::getDZ, 0.0, JLANG::JComparison::le() )
                                                );
        
        std::advance(__end = __begin, std::min(numberOfPrefits, (size_t) std::distance(__begin, __q)));
        
        std::partial_sort(__begin, __end, __q, qualitySorter);
      }
      
      out.erase(__end, out.end());
      
    } else {
      
      std::sort(out.begin(), out.end(), qualitySorter);
    }
    
  }
}
 
#endif