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 "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
 
#include "JDAQ/JDAQEventIO.hh"
#include "JDAQ/JDAQTimesliceIO.hh"
#include "km3net-dataformat/online/JDAQClock.hh"
 
#include "JSirene/JPulse.hh"
 
#include "JTools/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 "JTrigger/JTimeConverter.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);
        
        JTimeConverter 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);
 
          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();
}