JShowerPositionFit.cc

Go to the documentation of this file.

#include <iostream>
#include <string>

#include "Jeep/JParser.hh"

#include "evt/Head.hh"
#include "evt/Det.hh"

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

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

#include "JFit/JEvt.hh"
#include "JFit/JEvtToolkit.hh"
#include "JFit/JFitParameters.hh"
#include "JFit/JHistory.hh"
#include "JFit/JPoint4DRegressor.hh"

#include "JLang/JTypeList.hh"

#include "JSupport/JParallelFileScanner.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMeta.hh"
#include "JSupport/JSupport.hh"

#include "JTrigger/JHit.hh"
#include "JTrigger/JAlgorithm.hh"
#include "JTrigger/JBind2nd.hh"
#include "JTrigger/JTriggerParameters.hh"
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JBuildL2.hh"
#include "JTrigger/JMatch3G.hh"

#include "JGeometry3D/JTrack3D.hh"

#include <TFile.h>
#include <TH2D.h>

int main(int argc, char **argv){

  using namespace std;
  using namespace JSUPPORT;
  using namespace KM3NETDAQ;
  using namespace JFIT;
  using namespace JAANET;

  typedef JParallelFileScanner< JTypeList<JDAQEvent, JEvt> >         JParallelFileScanner_t;
  typedef JParallelFileScanner_t::multi_pointer_type                 multi_pointer_type;
  typedef JTYPELIST<JAllTypes_t, JEvt>::typelist                     typelist;

  JParallelFileScanner_t   inputFile;
  JFileRecorder<typelist>  outputFile;
  JLimit_t&                numberOfEvents = inputFile.getLimit();
  string                   detectorFile;
  string                   pdfFile;
  double                   Tmax_ns;
  double                   ctMin;
  double                   Dmax_m;           
  JRange<double>           pos_grid;
  JRange<double>           time_grid;
  double                   pos_step;
  double                   time_step;
  size_t                   numberOfPrefits;
  int                      debug;

  try{
    JParser<> zap;

    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)       = "PositionFit.root";
    zap['a'] = make_field(detectorFile);
    zap['F'] = make_field(pdfFile)          = "PDF_ShowerPositionFit.root";
    zap['n'] = make_field(numberOfEvents)   = JLimit::max();
    zap['N'] = make_field(Tmax_ns)          = 20.0;  // [ns] 
    zap['C'] = make_field(ctMin)            = -0.7;  // optimal value 
    zap['D'] = make_field(Dmax_m)           = 80.0;  // [m] 
    zap['P'] = make_field(pos_grid)         = JRange<double>(-28, 28);  // [m] position grid   
    zap['p'] = make_field(pos_step)         = 6;  // [m] density of the grid                  
    zap['T'] = make_field(time_grid)        = JRange<double>(-60, 60);  // [ns] time grid      
    zap['t'] = make_field(time_step)        = 15;  // [ns] density of the grid  
    zap['R'] = make_field(numberOfPrefits)  = 1;
    zap['d'] = make_field(debug)            = 2;
  
    if (zap.read(argc, argv) != 0) return 1;
  }
  catch(const exception& error){
    FATAL(error.what() << endl);
  }

  using namespace JTRIGGER;
  using namespace JDETECTOR;
  using namespace JGEOMETRY3D;
  using namespace JTOOLS;

  JDetector detector;

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

  TFile* pdf_file = new TFile(pdfFile.c_str());
  TH2D* hpdf = (TH2D*)pdf_file->Get("hPDF2Dist");

  const JModuleRouter moduleRouter(detector);

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

  const JBuildL0<JHitL0>  buildL0;
  const JBuildL2<JHitL1>  buildL2(3, Tmax_ns, ctMin);
  const JMatch3G<JHitL1>  match3G(Dmax_m, Tmax_ns); // CAUSALITY FOR SHOWERS 

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

  typedef JRegressor<JPoint4D, JSimplex> JRegressor_t;
  JRegressor_t::debug = debug;
  JRegressor_t::MAXIMUM_ITERATIONS = 10000;

  JRegressor_t fit0(true);
  fit0.estimator.reset(new JMEstimatorTukey(10.0));
  JRegressor_t fit(hpdf, false);
  
  for (JTreeScanner<Evt> in(inputFile); inputFile.hasNext(); ) {

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

    multi_pointer_type ps = inputFile.next();

    JDAQEvent* tev = ps;
    JEvt*      evt = ps;
    JEvt       out;

    JEvt::iterator __end = evt->end();
    if (numberOfPrefits > 0) {
      advance(__end = evt->begin(), min(numberOfPrefits, evt->size()));
    }
    partial_sort(evt->begin(), __end, evt->end(), qualitySorter);

    if(!evt->empty()) {
      
      copy(evt->begin(), __end, back_inserter(out));

      JDataL0_t dataL0;
      JDataL0_t dataL0_selected;
      JDataL1_t data;

      buildL0(*tev, moduleRouter, true, back_inserter(dataL0));
      buildL2(JDAQTimeslice(*tev, false), moduleRouter, back_inserter(data));

      for (JEvt::const_iterator shower = evt->begin(); shower != __end; ++shower) {
        
        const JPoint4D  vertex(getPosition(*shower), shower->getT());

        for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {     
          
          JHitL0 hit(*i);
          const JPosition3D hit_pos(hit.getPosition());
          const double D = hit_pos.getDistance(vertex);
          const double t_res = hit.getT() - vertex.getT(hit_pos);
          const JDirection3D photonDir(hit_pos - vertex.getPosition());
          const double cosT = photonDir.getDot(hit.getDirection());

          // hit selection
          if(data.size()>0){
            if(D < Dmax_m && (t_res > -40 && t_res < 40) && (cosT >= -1 && cosT <= 0.1)){
              dataL0_selected.push_back(hit);
            }
          } else {
            if(D < Dmax_m && (t_res > -50 && t_res < 50) && (cosT >= -1 && cosT <= 0.1)){
              dataL0_selected.push_back(hit);
            }
          }
          
        }

        if(dataL0_selected.size() != 0){

          double chi2 = 0.0;
          
          vector<double> LogChi2;
          vector<JPoint4D> start_positions;
          
          for(double x = pos_grid.getLowerLimit(); x <= pos_grid.getUpperLimit();
              x += pos_step){
            for(double y = pos_grid.getLowerLimit(); y <= pos_grid.getUpperLimit(); 
                y += pos_step){
              for(double z = pos_grid.getLowerLimit(); z <= pos_grid.getUpperLimit(); 
                  z += pos_step){
                for(double t = time_grid.getLowerLimit(); t <= time_grid.getUpperLimit(); 
                    t += time_step){
                  
                  chi2 = 0.0;
                  
                  JPosition3D pos_shifted(vertex.getX() + x, vertex.getY() + y, vertex.getZ() + z);
                  JPoint4D point_shifted(pos_shifted, vertex.getT() + t);                 
                  
                  for(JDataL0_t::const_iterator i = dataL0_selected.begin(); i != dataL0_selected.end(); ++i){
                    JHitL0 hit(*i);
                    chi2 += fit0(point_shifted, hit);
                  }     
                  
                  LogChi2.push_back(chi2);
                  start_positions.push_back(point_shifted);
                
                }
              }
            }
          }
          
          while(LogChi2.size() > 35){ 
            vector<double>::iterator it;
            it = max_element(LogChi2.begin(), LogChi2.end());
            int pos = it - LogChi2.begin();
            LogChi2.erase(LogChi2.begin() + pos);
            start_positions.erase(start_positions.begin() + pos);         
          }

          for(vector<JPoint4D>::iterator start_pos = start_positions.begin(); start_pos != start_positions.end(); ++start_pos){

            double simplex_step = 1;  // [m]        
            fit0.step.resize(4);
            fit0.step[0] = JPoint4D(JVector3D(simplex_step, 0.0, 0.0), 0.0);
            fit0.step[1] = JPoint4D(JVector3D(0.0, simplex_step, 0.0), 0.0);
            fit0.step[2] = JPoint4D(JVector3D(0.0, 0.0, simplex_step), 0.0);
            fit0.step[3] = JPoint4D(JVector3D(0.0, 0.0, 0.0), 20);

            chi2 = fit0(JPoint4D(*start_pos), dataL0_selected.begin(), dataL0_selected.end());

            JVector3D fitPos0(fit0.value);
            JPoint4D pt(fitPos0, fit0.value.getT());
            
            simplex_step = 0.3;  // [m]
            fit.step.resize(4);
            fit.step[0] = JPoint4D(JVector3D(simplex_step, 0.0, 0.0), 0.0);
            fit.step[1] = JPoint4D(JVector3D(0.0, simplex_step, 0.0), 0.0);
            fit.step[2] = JPoint4D(JVector3D(0.0, 0.0, simplex_step), 0.0);
            fit.step[3] = JPoint4D(JVector3D(0.0, 0.0, 0.0), 10);  
            
            chi2 = fit(pt, dataL0_selected.begin(), dataL0_selected.end());
            
            int NDF = dataL0_selected.size() - fit.step.size();
            
            JVector3D fitPos(fit.value);
            JTrack3D fitResult(fitPos, JVersor3Z(0,0), fit.value.getT());
            
            out.push_back(getFit(JHistory(shower->getHistory()).add(JSHOWERPOSITIONFIT), fitResult, getQuality(chi2, NDF), NDF));
            out.rbegin()->setW(13, chi2);
            out.rbegin()->setW(14, NDF);
          }
        } else{
          DEBUG("Too few hits " << endl);
        }
      }
      sort(out.begin(), out.end(), qualitySorter);
    }
    outputFile.put(out);
  }
  STATUS(endl);

  JMultipleFileScanner<JRemove<typelist, JEvt>::typelist> io(inputFile);
  io >> outputFile;
  outputFile.close();

}