JPreSim_HTR.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <map>
#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/JMultipleFileScanner.hh"
#include "JSupport/JParallelFileScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"

#include "JTools/JConstants.hh"
#include "JTools/JPermutation.hh"

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

#include "JGeometry3D/JOmega3D.hh"
#include "JGeometry3D/JRotation3D.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 monitor hit time residuals from reconstructed tracks using JPrefit+JSimplex
 *
 * 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;
 * With the option <tt>-I</tt>, tracks are fitted using only DOMs hit on the same DU as the excluded DOM.
 *   Additionally, using this option will also reduce the number of scan directions in JPrefit to a zenith scan only. 
 * \author mkarel
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  JMultipleFileScanner<JDAQEvent> inputFile;
  JLimit_t&      numberOfEvents = inputFile.getLimit();
  string         outputFileName;
  string         detectorFile;
  double         Tmax_ns;
  double         roadWidth_m;
  double         ctMin;
  int            numberOfOutliers;
  double         gridAngle_deg;
  size_t         numberOfPrefits;
  double         sigma_ns;
  bool           useL0;
  bool           noInterDU;
  int            minHits;
  int            debug;

  try { 

    JParser<> zap("Program to monitor hit time residuals from reconstructed tracks using JPrefit+JSimplex");
    
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFileName)      = "HTR_monitor.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['N'] = make_field(numberOfPrefits)     = 1;
    zap['U'] = make_field(useL0);
    zap['N'] = make_field(numberOfPrefits)     = numeric_limits<size_t>::max();
    zap['L'] = make_field(minHits)             = 4;
    zap['I'] = make_field(noInterDU);
    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();
  const double TMAX_NS                = 50.0;       // [ns]

  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;

  JOmega3D_t                 omega ;
  if (noInterDU) {

    const double thetaMin = 0.75*PI ;
    const double thetaMax = PI ;
    const double rad  = thetaMax - thetaMin;
    const double bin  = rad / floor(rad/(gridAngle_deg*PI/180.0) + 0.5);        // polar angle step size

    for (double theta = thetaMin; theta < thetaMax + 0.5*bin; theta += bin) {
      omega.push_back(JAngle3D(theta,0.0));
    }

  } else {

    omega = JOmega3D(gridAngle_deg * PI/180.0) ;

  }

  NOTICE( "Scanning over " << omega.size() << "JPrefit directions" <<endl ) ;

  typedef JRegressor<JLine3Z, JSimplex>  JRegressor_t;

  JRegressor_t fit(sigma_ns);

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

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

  TFile outputFile(outputFileName.c_str(),"recreate");

  // exclude one DOM from the fit
  for (JModuleRouter::const_iterator excl_mod = detector.begin(); excl_mod != detector.end(); ++excl_mod) {

    // output tree (one per excluded DOM)
    double htr_dt ;
    JFit bestfit ;
    TTree htr_tree(Form("%d.htr",excl_mod->getID()),Form("%d.htr",excl_mod->getID())) ;
    htr_tree.Branch("dt",&htr_dt) ;
    htr_tree.Branch("bestfit",&bestfit) ;

    STATUS( "Excluding DOM "<<excl_mod->getID()<<endl ) ;

  while (inputFile.hasNext()) {

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

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


    // L1 and L0 data
                
    JDataL0_t dataL0;
    JDataL1_t dataL1;
    JDataL1_t dataL1_excluded;
    JDataR1_t data;
    
    JDAQTimeslice timeslice(*tev, true);

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

  JModule module = router.getModule(i->getModuleID()) ;
  buffer(*i, module);

  if (module.getID() != excl_mod->getID()) {

    if (!noInterDU || module.getString() == excl_mod->getString()) {

      buildL0(buffer, back_inserter(dataL0));
      buildL2(buffer, back_inserter(dataL1));

    }

  } else {

    buildL2(buffer, back_inserter(dataL1_excluded));

  }
    }


    // 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) >= minHits && dataL1_excluded.size() > 0) {

      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) < minHits) {
          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)) >= minHits) {
                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(JPRESIM, 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);

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

      } else {

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

// JSimplex

    if (!out.empty()) {

      JEvt::iterator __end = out.end();
 
      if (numberOfPrefits > 0) {
        advance(__end = out.begin(), min(numberOfPrefits, out.size()));
      }
 
      for (JEvt::const_iterator track = out.begin(); track != __end; ++track) {

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

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

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

              JHitR1 hit(*i);
              
              hit.rotate(R);
              
              if (match(hit)) {
                data.push_back(hit);
              }
            }
          }
        }


        fit.step.resize(5);
        
        fit.step[0] = JLine3Z(JLine1Z(JVector3D(0.5, 0.0, 0.0), 0.0));
        fit.step[1] = JLine3Z(JLine1Z(JVector3D(0.0, 0.5, 0.0), 0.0));
        fit.step[2] = JLine3Z(JLine1Z(JVector3D(0.0, 0.0, 0.0), 1.0));
        fit.step[3] = JLine3Z(JLine1Z(JVector3D(), 0.0), JVersor3Z(0.005, 0.0));
        fit.step[4] = JLine3Z(JLine1Z(JVector3D(), 0.0), JVersor3Z(0.0, 0.005));
        
        const int NDF = getCount(data.begin(), data.end()) - fit.step.size();

        if (NDF >= 0) {
              
          const double chi2 = fit(JLine3Z(tz), data.begin(), data.end());
              
          JTrack3D tb(fit.value);
              
          tb.rotate_back(R);

          out2.push_back(getFit(JPRESIM, tb, getQuality(chi2, NDF), NDF)); 
        }
      }

      // apply default sorter

      sort(out2.begin(), out2.end(), qualitySorter);

    }  // !out.empty()

    }  // distance(dataL1.begin(), __end) >= minHits

    if (!out2.empty()) {

      bestfit = *out2.begin() ;

      const double t_exp = getTrack(bestfit).getT(excl_mod->getPosition()) ;

      sort(dataL1_excluded.begin(), dataL1_excluded.end()) ;

      htr_dt = dataL1_excluded.begin()->getT()-t_exp ;

      htr_tree.Fill() ;
 
    }

  }  // inputFile.hasNext()

  inputFile.rewind() ;

  outputFile.cd() ;
  htr_tree.Write() ;

  } // excluded DOM

  STATUS(endl);

  outputFile.Write() ;

}