JShowerFit.cc

Go to the documentation of this file.

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

#include "evt/Head.hh"
#include "evt/Det.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/JHitL0.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JTriggerParameters.hh"

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

#include "JFit/JHitW0.hh"
#include "JFit/JPMTW0.hh"
#include "JFit/JGandalf.hh"
#include "JFit/JFitToolkit.hh"
#include "JFit/JEvt.hh"
#include "JFit/JEvtToolkit.hh"
#include "JFit/JFitParameters.hh"
#include "JFit/JShower3EZ.hh"
#include "JFit/JRegressor.hh"
#include "JFit/JShower3EZRegressor.hh"

#include "JTools/JConstants.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"

#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JPDFToolkit.hh"

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

#include <JGeometry3D/JOmega3D.hh>

#include "JAAnet/JAAnetToolkit.hh"


/**
 * \file
 *
 * Program to perform EM Shower Fit for ORCA with I/O of JFIT::JEvt data.
 * The reconstruction is made at the PMT level.
 * \author adomi
 */

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

  using namespace std;
  using namespace JFIT;
  using namespace JSUPPORT;
  using namespace KM3NETDAQ;
  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         Dmax_m;
  double         R_Hz;
  size_t         numberOfPrefits;  
  int            debug;


   try { 

    JParser<> zap;
   
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)          = "gandalf.root";
    zap['a'] = make_field(detectorFile);
    zap['F'] = make_field(pdfFile)             = "quickPDF_elec-CC.root";
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['N'] = make_field(Tmax_ns)             = 20;  // [ns]
    zap['D'] = make_field(Dmax_m)              = 80;  // [m]
    zap['B'] = make_field(R_Hz)                = 10.0e3;
    zap['R'] = make_field(numberOfPrefits)     = 1; 
    zap['d'] = make_field(debug)               = 2;
   
    zap(argc, argv);
  }
  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);
  }
  
  const JModuleRouter moduleRouter(detector);
  
  typedef vector<JHitL0>  JDataL0_t;
  typedef vector<JPMTW0> JPMTW0_t;
  const JBuildL0<JHitL0>  buildL0;

  typedef JRegressor<JShower3EZ, JSimplex>   JRegressor_t;
  
  JRegressor_t::debug  = debug;
  JRegressor_t::T_ns.setRange(-Tmax_ns, Tmax_ns);
  JRegressor_t::Vmax_npe = 20.0;   
  JRegressor_t::MAXIMUM_ITERATIONS = 1000;
  JRegressor_t fit(pdfFile);

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

  JTreeScanner<Evt> in(inputFile);

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

    if(!evt->empty()) {

      JEvt::iterator __end = evt->end();

      copy(evt->begin(), __end, back_inserter(out));
      
      if (numberOfPrefits > 0) {
        advance(__end = evt->begin(), min(numberOfPrefits, evt->size()));
      }
 
      partial_sort(evt->begin(), __end, evt->end(), qualitySorter);
      
      JDataL0_t dataL0;
      JDAQTimeslice timeslice(*tev, true);
        
      buildL0(*tev, moduleRouter, true, back_inserter(dataL0));
       
      for (JEvt::const_iterator shower = evt->begin(); shower != __end; ++shower) {

        JVector3D pos(getPosition(*shower)); //.rotate(R));
        
        JPoint4D pt(pos, shower->getT());
        JVector3D start_dir(0,0,0);       

        JVector3D DetC(0.0, 0.0, 117);
        double radius = pos.getDistance(DetC);
        
        multiset<JDAQPMTIdentifier> top;

        double Nhits = 0;

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

          if(D < Dmax_m && (t_res >= -25 && t_res <= 25) && (cosT >= -1 && cosT <= 0.1)){  
          top.insert(hit.getPMTIdentifier());
          JVector3D d(photonDir.getDX(), photonDir.getDY(), photonDir.getDZ());
          start_dir.add(d);
          }
          if(D < Dmax_m && (t_res > -30 && t_res < 20)){
          Nhits++;
          }
        }

        double start_E = exp((Nhits + 26)/26);  // starting energy [GeV]
        if(start_E > 100) start_E = 100;
        
        vector<double> LogLik, Energies;
        vector<JVector3D> start_directions;
        vector<JPMTW0> buffer;
        double chi2 = 0.0;
        double max_theta = 0, max_phi = 0, scan_step = 0;
        double Nmin = 0;

        if(start_E > 10 && radius < 100){
          max_theta = 30;
          max_phi = 30;
          scan_step = 5;
          Nmin = 10;
        } else {
          max_theta = 35;
          max_phi = 35;
          scan_step = 5;
          Nmin = 20;
          }

        JAngle3D start_dir_angles(start_dir.getX(), start_dir.getY(), start_dir.getZ());

        for(double E = -15; E <= 5; E += 7.5){  // scan in energy
          for(double th = -max_theta; th <= max_theta; th += scan_step){
            for(double ph = -max_phi; ph <= max_phi; ph += scan_step){

              buffer.clear();
              chi2 = 0.0;
            
              double theta = th * 3.14 / 180;
              double phi = ph * 3.14 / 180;
              double scan_E = start_E + E;  // [GeV]
              
              JAngle3D dir_shifted_angles(start_dir_angles.getTheta() + theta, start_dir_angles.getPhi() + phi);
              JVector3D dir_shifted(dir_shifted_angles.getDX(), dir_shifted_angles.getDY(), dir_shifted_angles.getDZ()); 
              JDirection3D conversion(dir_shifted.getX(), dir_shifted.getY(), dir_shifted.getZ());
              JRotation3D R(conversion);
              pt.rotate(R);
            
              for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
                JPosition3D pos(module->getPosition());
                pos.rotate(R);
                if (pt.getDistance(pos) < Dmax_m) {
                  for (unsigned int i = 0; i != module->size(); ++i) {
                    JDAQPMTIdentifier id(module->getID(), i);
                    JPMT pmt(module->getPMT(i));
                    pmt.rotate(R);
                    buffer.push_back(JPMTW0(pmt, R_Hz, top.count(id)));
                  }
                }
              }
              
              for(JPMTW0_t::iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt ){
                chi2 += fit(JShower3EZ(pt, JVersor3Z(), JEnergy(log10(scan_E))), *pmt);
              }

              LogLik.push_back(chi2);
              start_directions.push_back(dir_shifted);
              Energies.push_back(scan_E);
              pt.rotate_back(R);
              
            } 
          } 
        } 
        
        while(LogLik.size() > Nmin){

          vector<double>::iterator it;
          it = max_element(LogLik.begin(), LogLik.end());
          int pos = it - LogLik.begin();
          LogLik.erase(LogLik.begin() + pos);
          start_directions.erase(start_directions.begin() + pos);
          Energies.erase(Energies.begin() + pos);

        }
        
        int n = 0;
        for(vector<JVector3D>::iterator dir = start_directions.begin(); dir != start_directions.end(); ++dir){    

          fit.step.resize(3);
          double dir_step = 0.05;
          fit.step[0] = JShower3EZ(JPoint4D(JVector3D(), 0.0), JVersor3Z(dir_step, 0.0), JEnergy());
          fit.step[1] = JShower3EZ(JPoint4D(JVector3D(), 0.0), JVersor3Z(0.0, dir_step), JEnergy());
          fit.step[2] = JShower3EZ(JPoint4D(JVector3D(), 0.0), JVersor3Z(), 0.05);

          double scan_E = Energies[n];
          n++;
          buffer.clear();
          
          JVersor3D startVersor(dir->getX(), dir->getY(), dir->getZ());
          
          JDirection3D conversion(startVersor.getDX(), startVersor.getDY(), startVersor.getDZ());
          JRotation3D R(conversion);
          pt.rotate(R);
          
            for(JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
              JPosition3D pos(module->getPosition());
              pos.rotate(R);
              if (pt.getDistance(pos) < Dmax_m) {
                for (unsigned int i = 0; i != module->size(); ++i) {
                  JDAQPMTIdentifier id(module->getID(), i);
                  JPMT pmt(module->getPMT(i));
                  pmt.rotate(R);
                  buffer.push_back(JPMTW0(pmt, R_Hz, top.count(id)));
                }
              }
            }
            
            int NDF = buffer.size() - fit.step.size();
            
            chi2 = fit(JShower3EZ(pt, JVersor3Z(), JEnergy(log10(scan_E))), buffer.begin(), buffer.end());  
            
            pt.rotate_back(R); 
            
            JVector3D resPos(fit.value);
            JVersor3D resDir(fit.value.getDX(), fit.value.getDY(), fit.value.getDZ());  
            
            JTrack3D td(resPos, resDir, shower->getT());          
            double E_reco_corrected = fit.value.getE() - 4;
            JTrack3E tb(td, E_reco_corrected); 
            tb.rotate_back(R);                                        
            
            out.push_back(getFit(JHistory(shower->getHistory()).add(JSHOWERCOMPLETEFIT), tb, -chi2, NDF));
            out.rbegin()->setE(E_reco_corrected);
        }
      }
      sort(out.begin(), out.end(), qualitySorter);
    }
    outputFile.put(out);
  }
  STATUS(endl);
  
  JMultipleFileScanner<JRemove<typelist, JEvt>::typelist> io(inputFile);
  io >> outputFile;
  outputFile.close();

}