JFilter.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include <map>
#include <algorithm>
#include <iterator>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TProfile.h"

#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Hit.hh"
#include "km3net-dataformat/tools/time_converter.hh"

#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"

#include "JDAQ/JDAQEventIO.hh"
#include "JDAQ/JDAQTimesliceIO.hh"
#include "km3net-dataformat/online/JDAQClock.hh"

#include "JSirene/JPulse.hh"

#include "JPhysics/JConstants.hh"

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

#include "JTrigger/JHit.hh"
#include "JTrigger/JFrame.hh"
#include "JTrigger/JTimeslice.hh"
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JHitL1.hh"
#include "JTrigger/JBuildL1.hh"
#include "JTrigger/JBuildL2.hh"
#include "JTrigger/JAlgorithm.hh"
#include "JTrigger/JMatch3B.hh"
#include "JTrigger/JMatch3D.hh"
#include "JTrigger/JMatch1D.hh"
#include "JTrigger/JBind2nd.hh"

#include "JSupport/JTriggeredFileScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"

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


namespace {

  using namespace JPP;


  /**
   * Data structure for Monte Carlo hits.
   */
  class JBuffer :
    public std::vector<JPulse>
  {
  public: 
    /**
     * Default constructor.
     */
    JBuffer() :
      std::vector<JPulse>()
    {}


    /**
     * Compile data to L1 hits.
     *
     * An L1 hit consists of two (or more) hits from different PMTs which are coincident in time.
     * The start(stop) time of the L1 hit is defined by the earliest(latest) time of the hits 
     * that constitute the L1 hit.
     * The PMT identifier of the earliest hit is maintained.
     * The compilation may reduce the number of hits.
     *
     * \param  Tmax       maximal time difference between hits [ns]
     * \param  L1         minimum number of L0 hits to trigger L1
     */
    void compile(const double       Tmax,
                 const unsigned int L1 = 2)
    {
      std::sort(this->begin(), this->end());

      iterator out = this->begin();

      for (const_iterator i = this->begin(); i != this->end(); ) {
        
        if        (L1 == 1) {

          *out = JPulse(*i, *i);
          ++out;
          
          ++i;

        } else if (L1 >  1) {

          std::set<int> buffer;

          buffer.insert(i->getID());
          
          const_iterator j = i;
          
          for (double t1 = i->getLowerLimit(); ++j != this->end() && j->getLowerLimit() - t1 < Tmax; t1 = j->getLowerLimit()) {
            buffer.insert(j->getID());
          }
          
          if (buffer.size() >= L1) {      
            *out = JPulse(*i, *((--j)++));
            ++out;
          }
          
          i = j;
        }
      }

      this->erase(out, end());
    }


    /**
     * Check if a hit with given time range is present in buffer.
     *
     * \param  t1         minimal time [ns]
     * \param  t2         maximal time [ns]
     * \return            true if hit present; else false
     */
    inline bool has(const JTimeRange& timerange) const
    {
      const_iterator i = std::lower_bound(this->begin(), this->end(), timerange.getLowerLimit());

      return i != this->end() && i->overlap(timerange);
    }
  };


  using namespace JPP;

  /**
   * Match all L1 hits.
   *
   * \param  first      first  hit
   * \param  second     second hit
   * \return            true
   */
  inline bool matchAll(const JHitL1& first, const JHitL1& second)
  {
    return true;
  }


  /**
   * Auxiliary class to select PMTs looking below a given horizon.
   */
  struct JHorizon {
    /**
     * Constructor.
     * 
     * \param  ct               cosine elevation angle of horizon
     * \param  N                mininum number of PMTs below horizon
     */
    JHorizon(const double ct,
             const int    N)
    {
      this->ct = ct;
      this->N  = N;
    }

    
    /**
     * Select L1 hit.
     * 
     * \param  hit              L1 hit
     * \return                  true if accepted; false if rejected
     */
    bool operator()(const JHitL1& hit) const
    {
      int n = 0;

      for (JHitL1::const_iterator i = hit.begin(); i != hit.end(); ++i) {
        if (i->getDZ() < ct)
          ++n;
      }

      return n >= N;
    }


    double ct;
    int    N;
  };
}


/**
 * \file
 * Example program to test performance of various hit filters based on JTRIGGER::clusterize method.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  JTriggeredFileScanner<> inputFile;
  JLimit_t&      numberOfEvents = inputFile.getLimit();
  string         outputFile;
  string         detectorFile;
  double         Tmax_ns;
  double         roadWidth_m;
  double         ctMin;
  double         thetaRad;
  char           cluster;
  int            histogram;
  int            debug;

  try { 

    JParser<> zap("Example program to test performance of various hit filters.");
    
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)          = "filter.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['Q'] = make_field(thetaRad)            =   0.2;               //   
    zap['c'] = make_field(cluster)             = 'A', 'B', 'C', 'D', 'E', 'F';
    zap['H'] = make_field(histogram)           = 1, 2;
    zap['d'] = make_field(debug)               = 1;
    
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  using namespace KM3NETDAQ;

  cout.tie(&cerr);


  JDetector detector;

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

  const JModuleRouter moduleRouter(detector);
  const JPMTRouter    pmtRouter   (detector);


  detector -= get<JPosition3D>(getHeader(inputFile));


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

  TProfile* he;
  TProfile* hp;

  if (histogram == 1) {
 
    vector<double> X;

    double x = -0.5;

    for ( ; x <  30.0; x +=  1.0)
      X.push_back(x);

    for ( ; x <  50.0; x +=  2.0)
      X.push_back(x);

    for ( ; x < 100.0; x +=  5.0)
      X.push_back(x);

    for ( ; x < 200.0; x += 10.0)
      X.push_back(x);

    he = new TProfile("he", NULL, X.size() - 1, X.data());
    hp = new TProfile("hp", NULL, X.size() - 1, X.data());

  } else {

    he = new TProfile("he", NULL, 28, 0.0, 7.0);
    hp = new TProfile("hp", NULL, 28, 0.0, 7.0);
  }

  TH1D ht1("ht1", NULL, 550, -50.0, +500.0);
  TH1D ht2("ht2", NULL, 550, -50.0, +500.0);
  TH1D hd1("hd1", NULL, 100,  -1.0,   +1.0);
  TH1D hd2("hd2", NULL, 100,  -1.0,   +1.0);
  TH1D hx1("hx1", NULL, 100,   0.0, +250.0);
  TH1D hx2("hx2", NULL, 100,   0.0, +250.0);
  TH1D hw1("hw1", NULL, 100,  -0.5,  +99.5);
  TH1D hw2("hw2", NULL, 100,  -0.5,  +99.5);



  typedef map<int, JBuffer> JMap_t;
  typedef vector<JHitL1>    JDataL1_t;

  JMap_t                 zmap;
  const JBuildL2<JHitL1> buildL2(JL2Parameters(2, Tmax_ns, ctMin));
  const JMatch3B<JHitL1> match3B(roadWidth_m, Tmax_ns);
  const JMatch3D<JHitL1> match3D(Tmax_ns);
  const JMatch1D<JHitL1> match1D(roadWidth_m, Tmax_ns);
  const JHorizon         horizon(0.2, 2);

  const JVersor3D        gui(getSinThetaC(), 0.0, getCosThetaC());   // photon emission direction


  while (inputFile.hasNext()) {

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

    JTriggeredFileScanner<>::multi_pointer_type ps = inputFile.next();

    const JDAQEvent* tev   = ps;
    const Evt*       event = ps;
    
    if (has_neutrino(*event)) {

      const Trk&   neutrino = get_neutrino(*event);
      const double E        = neutrino.E;              // GeV


      vector<Trk>::const_iterator muon = event->mc_trks.end();

      {
        double Emax = 0.0;

        for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
          
          // (anti) muon
          
          if (is_muon(*track) && track->E > Emax) {
            muon = track;
            Emax = track->E;
          }
        }
      }

      if (muon != event->mc_trks.end()) {

        JPosition3D  pos = getPosition (*muon);
        double       t0  = muon->t;
        
        const JRotation3D R(getDirection(*muon));
        
        pos.rotate(R);
        
        const time_converter converter(*event, *tev);

        // Monte Carlo hits

        zmap.clear();

        const vector<Hit>& hitlist = event->mc_hits;

        for (vector<Hit>::const_iterator i = hitlist.begin(); i != hitlist.end(); ++i) {
          if (!is_noise(*i)) {
            zmap[pmtRouter.getParentModuleID(i->pmt_id)].push_back(*i);
          }
        }


        // L1 Monte Carlo true data

        for (JMap_t::iterator i = zmap.begin(); i != zmap.end(); ) {

          i->second.compile(Tmax_ns, 2);

          if (i->second.empty())
            zmap.erase(i++);
          else
            ++i;
        }

        const int L1mc = zmap.size();


        // L1 triggered data

        JDataL1_t dataL1;
        JDataL1_t zbuf;

        JDAQTimeslice timeslice(*tev, true);

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

          zbuf.clear();

          buildL2(*i, moduleRouter.getModule(i->getModuleID()), back_inserter(zbuf));

          if (!zbuf.empty()) {

            sort(zbuf.begin(), zbuf.end(), timeSorter<JHitL1>);
            
            dataL1.push_back(*zbuf.begin());
          }
        }


        // cluster

        JDataL1_t::iterator __end = dataL1.begin();

        switch (cluster) {

        case 'A':
          //__end = dataL1.end();
          __end = clusterize(dataL1.begin(), dataL1.end(), make_match(matchAll));
          break;

        case 'B':
          __end = clusterize(dataL1.begin(), dataL1.end(), weightSorter<JHitL1>, match3D);
          break;

        case 'C':
          __end = clusterize(dataL1.begin(), dataL1.end(), match3B);
          break;

        case 'D':
          __end = clusterizeWeight(dataL1.begin(), dataL1.end(), match3B);
          break;

        case 'E':
        
          for (JDataL1_t::iterator i = dataL1.begin(); i != dataL1.end(); ++i)
            i->rotate(R);

          __end = clusterizeWeight(dataL1.begin(), dataL1.end(), match1D);
        
          for (JDataL1_t::iterator i = dataL1.begin(); i != dataL1.end(); ++i)
            i->rotate_back(R);

          break;

        case 'F':
        
          __end = clusterizeWeight(dataL1.begin(), dataL1.end(), match3B);

          for (JDataL1_t::iterator i = dataL1.begin(); i != dataL1.end(); ++i)
            i->rotate(R);

          __end = partition(dataL1.begin(), __end, horizon);
          __end = clusterizeWeight(dataL1.begin(), __end, match1D);

          for (JDataL1_t::iterator i = dataL1.begin(); i != dataL1.end(); ++i)
            i->rotate_back(R);

          break;
        
        default:
          __end = dataL1.end();
          break;
        }


        // L1 monitor

        JDataL1_t::iterator __q = __end;

        set<int> L1;
        set<int> L1ok;

        for (vector<JHitL1>::iterator hit = dataL1.begin(); hit != __q; ) {

          L1.insert(hit->getModuleID());

          const double t1 = converter.getTime(hit->getT());

          if (zmap[hit->getModuleID()].has(JTimeRange(t1 - 0.5*Tmax_ns, t1 + 0.5*Tmax_ns))) {

            L1ok.insert(hit->getModuleID());

            ++hit;

          } else {

            iter_swap(hit, --__q);
          }
        }


        if (L1mc != 0 && !L1.empty()) {

          Double_t x = numeric_limits<Double_t>::max();

          if (histogram == 1)
            x = L1mc;
          else
            x = log10(E);

          he->Fill(x, (Double_t) L1ok.size() / (Double_t) L1mc); 
          hp->Fill(x, (Double_t) L1ok.size() / (Double_t) L1.size()); 


          // transform hits to muon system
        
          for (JDataL1_t::iterator hit = dataL1.begin(); hit != __end; ++hit) {
            hit->transform(R, pos);
          }
        
        
          for (JDataL1_t::iterator hit = dataL1.begin(); hit != __end; ++hit) {
          
            const double t1 = t0  +  getInverseSpeedOfLight() * (hit->getZ() + hit->getX() * getTanThetaC());
          
            double dot = +1.0;

            for (JHitL1::const_iterator i = hit->begin(); i != hit->end(); ++i) {
              if (i->getDot(gui) < dot)
                dot = i->getDot(gui);
            }

            const double w = (hit->rbegin()->getT() - hit->begin()->getT()); // / hit->size();

            if (distance(hit, __q) > 0) {
              ht1.Fill(converter.getTime(hit->getT())  -  t1);
              hd1.Fill(dot);
              hx1.Fill(hit->getX());
              hw1.Fill(w);
            } else {
              ht2.Fill(converter.getTime(hit->getT())  -  t1);
              hd2.Fill(dot);
              hx2.Fill(hit->getX());
              hw2.Fill(w);
            }
          }
        }
      }
    }
  }
  STATUS(endl);

  out.Write();
  out.Close();
}