JPrefit.cc

Go to the documentation of this file.

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

#include "evt/Head.hh"
#include "evt/Evt.hh"
#include "JDAQ/JDAQEvent.hh"
#include "JDAQ/JDAQTimeslice.hh"
#include "JDAQ/JDAQSummaryslice.hh"

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"

#include "JTrigger/JHit.hh"
#include "JTrigger/JFrame.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/JBuildL1.hh"
#include "JTrigger/JBuildL2.hh"
#include "JTrigger/JAlgorithm.hh"
#include "JTrigger/JMatch1D.hh"
#include "JTrigger/JMatch3B.hh"
#include "JTrigger/JBind2nd.hh"
#include "JTrigger/JTriggerParameters.hh"

#include "JSupport/JSingleFileScanner.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.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 "JTools/JConstants.hh"
#include "JTools/JPermutation.hh"

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

#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


namespace {

  using namespace JPP;


  /**
   * Sort hits according times corrected for position along z-axis.
   *
   * \param  first             first  hit
   * \param  second            second hit
   * \return                   true if first hit earlier than second hit; else false
   */
  inline bool cmz(const JHitR1& first, const JHitR1& second)
  {
    return (first .getT() * getSpeedOfLight()  -  first .getZ()  <
            second.getT() * getSpeedOfLight()  -  second.getZ());
  }
}
               

/**
 * \file
 *
 * Program to perform linear fits JFIT::JEstimator<JLine1Z> covering the solid angle producing JFIT::JEvt data.
 *
 * The option <tt>-U</tt> has the following consequences:
 *    - L0 hits are included in the fit;
 *    - L1 hits are obtained from the list of triggered hits;
 *    - quality of the fit is determined counting all L0 hits;
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  typedef JTYPELIST<JRemove<JAllTypes_t, JDAQTimesliceTypes_t>::typelist, JEvt>::typelist   typelist;

  JSingleFileScanner<JDAQEvent> inputFile;
  JFileRecorder     <typelist>  outputFile;
  JLimit_t&      numberOfEvents = inputFile.getLimit();
  string         detectorFile;
  double         Tmax_ns;
  double         roadWidth_m;
  double         ctMin;
  double         sigma_ns;
  int            numberOfOutliers;
  double         gridAngle_deg;
  bool           useL0;
  size_t         numberOfPrefits;
  int            debug;

  try { 

    JParser<> zap("Program to perform pre-fit of muon trajectory to data.");
    
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)          = "prefit.root";
    zap['a'] = make_field(detectorFile);
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['T'] = make_field(Tmax_ns)             =  20.0;               // [ns]
    zap['R'] = make_field(roadWidth_m)         = 150.0;               // [m]    
    zap['C'] = make_field(ctMin)               =   0.0;               //   
    zap['S'] = make_field(sigma_ns);
    zap['O'] = make_field(numberOfOutliers)    = 0;
    zap['g'] = make_field(gridAngle_deg);                             // [deg]
    zap['U'] = make_field(useL0);
    zap['N'] = make_field(numberOfPrefits)     = numeric_limits<size_t>::max();
    zap['d'] = make_field(debug)               = 1;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }



  cout.tie(&cerr);

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


  JDetector detector;

  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }

  const JModuleRouter router(detector);


  typedef vector<JHitL0>     JDataL0_t;
  typedef vector<JHitL1>     JDataL1_t;
  typedef vector<JHitR1>     JDataR1_t;

  JSuperFrame2D <JHit>       buffer;
  const JBuildL0<JHitL0>     buildL0;
  const JBuildL2<JHitL1>     buildL2(2, Tmax_ns, ctMin);
  const JMatch1D<JHitR1>     match1D(roadWidth_m, Tmax_ns);
  const JMatch3B<JHitL1>     match3B(roadWidth_m, Tmax_ns);
  const JHitL1Comparator     compare;
  const JOmega3D             omega(gridAngle_deg * PI/180.0);


  outputFile.open();

  outputFile.put(JMeta(argc, argv));

  while (inputFile.hasNext()) {

    STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);

    JDAQEvent* tev   = inputFile.next();
    JEvt       out;


    // L1 and L0 data
                
    JDataL0_t dataL0;
    JDataL1_t dataL1;
    
    buildL0(JDAQTimeslice(*tev,  true),  router, back_inserter(dataL0));
    buildL2(JDAQTimeslice(*tev, !useL0), router, back_inserter(dataL1));


    // 3D cluster of unique optical modules

    JDataL1_t::iterator __end = unique(dataL1.begin(), dataL1.end(), equal_to<JDAQModuleIdentifier>());
    
    __end = clusterizeWeight(dataL1.begin(), __end, match3B);


    if (distance(dataL1.begin(), __end) >= NUMBER_OF_L1HITS) {


      JDataR1_t data;


      for (JOmega3D_t::const_iterator dir = omega.begin(); dir != omega.end(); ++dir) {

        const JRotation3D R(*dir);

        data.clear();

        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 (distance(data.begin(), __end1) > MAXIMUM_NUMBER_OF_HITS) {

          advance(__end1 = data.begin(), MAXIMUM_NUMBER_OF_HITS);

          partial_sort(data.begin(), __end1, data.end(), cmz);
        }


        // 1D cluster

        __end1 = clusterizeWeight(data.begin(), __end1, match1D);

        if (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 (find_if(data.begin(), __end1, bind2nd(equal_to<JDAQModuleIdentifier>(), i->getModuleID())) == __end1) {

              JHitR1 hit(*i);

              hit.rotate(R);

              if (count_if(data.begin(), __end1, JBind2nd(match1D,hit)) >= NUMBER_OF_L1HITS) {
                data.push_back(hit);
              }
            }
          }

          __end1 = clusterize(__end1, data.end(), match1D);
        }

          
        if (distance(data.begin(), __end1) <= JEstimator<JLine1Z>::NUMBER_OF_PARAMETERS) {
          continue;
        }

        // 1D fit
          
        JLine1Z  tz;
        double   chi2 = numeric_limits<double>::max();
        int      NDF  = 0;
        int      N    = 0;
          
        JMatrixNZ V;
        JVectorNZ Y;
            
        if (distance(data.begin(), __end1) <= FACTORY_LIMIT) {
            
          double ymin = numeric_limits<double>::max();
            
          JDataR1_t::iterator __end2 = __end1;
              
          for (int n = 0; n <= numberOfOutliers && distance(data.begin(), __end2) > JEstimator<JLine1Z>::NUMBER_OF_PARAMETERS; ++n, --__end2) {
              
            sort(data.begin(), __end1, compare);
              
            do {

              try {
                
                JEstimator<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) {
                  WARNING(endl << "chi2(1) " << y << endl);
                } else if (y < ymin) {
                  ymin = y;
                  tz   = fit;
                  NDF  = distance(data.begin(), __end2) - JEstimator<JLine1Z>::NUMBER_OF_PARAMETERS;
                  N    = getCount(data.begin(), __end2);
                  chi2 = y;
                }
              }
              catch(const JException& error) {}

            } while (next_permutation(data.begin(), __end2, __end1, compare));
              
            ymin -= STANDARD_DEVIATIONS * STANDARD_DEVIATIONS;
          }
            
        } else {
            
          try {
              
            const int number_of_outliers = distance(data.begin(), __end1) - JEstimator<JLine1Z>::NUMBER_OF_PARAMETERS - 1;
              
            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  = distance(data.begin(), __end1) - JEstimator<JLine1Z>::NUMBER_OF_PARAMETERS;
            N    = 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 = getChi2(Y, V, i);
                  
                if        (y <= 0.0) {
                  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 -= getCount(data[kill]);
            }

            chi2 = getChi2(Y, V);
            tz   = fit;
          }
          catch(const JException& error) {}
        }

        if (NDF != 0) {

          tz.rotate_back(R);

          out.push_back(getFit(JHistory(JMUONPREFIT), tz, *dir, getQuality(chi2, (useL0 ? N : NDF)), NDF));
        }
      }

      // apply default sorter

      if (numberOfPrefits > 0) {

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

        advance(__end = out.begin(), min(numberOfPrefits, out.size()));

        partial_sort(out.begin(), __end, out.end(), qualitySorter);

        // add downward pointing solutions if available but absent so far

        if (count_if(out.begin(), __end, make_predicate(&JFit::getDZ, 0.0, JComparison::le())) == 0) {

          JEvt::iterator __begin = __end;
          JEvt::iterator __q     = partition(__begin, out.end(), make_predicate(&JFit::getDZ, 0.0, JComparison::le()));

          advance(__end = __begin, min(numberOfPrefits, (size_t) distance(__begin, __q)));

          partial_sort(__begin, __end, __q, qualitySorter);
        }

        out.erase(__end, out.end());

      } else {

        sort(out.begin(), out.end(), qualitySorter);
      }
    }

    outputFile.put(out);
  }
  STATUS(endl);

  JSingleFileScanner<JRemove<typelist, JEvt>::typelist> io(inputFile);

  io >> outputFile;

  outputFile.close();
}