JFilter.cc File Reference

#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"

Go to the source code of this file.

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

Detailed Description

Example program to test performance of various hit filters based on JTRIGGER::clusterize method.

Author

mdejong

Definition in file JFilter.cc.

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 210 of file JFilter.cc.

{
  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();
}