JAnglerFish.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>

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

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

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

#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"

#include "JROOT/JManager.hh"

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

#include "JTrigger/JSuperFrame2D.hh"
/**
 * \file
 *
 * Example application to search for sets of consecutive hits on a given PMT,\n
 * that are separated in time by a constant value. The number of consecutive hits\n
 * as well as the time difference are selected from the command line.\n 
 * The searches are done within the JDAQSnapshotHits from JDAQEvents\
 * as well as within the JDAQHits from the JDAQTimeslice data.
 *
 * \author rgruiz
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  typedef JTYPELIST<JMeta , TH1D, TH2D>::typelist  typelist;
 
  string                  inputFile;
  JLimit_t                numberOfEvents;
  JFileRecorder<typelist> outputFile;
  string                  detectorFile;
  int                     deltaT;
  int                     multiplicity;
  int                     debug;

  try { 

    JParser<> zap("Program to search for sets of consecutive hits on a given PMT that are equidistant in time.");
    
    zap['f'] = make_field(inputFile);
    zap['a'] = make_field(detectorFile);
    zap['t'] = make_field(deltaT)         = 255;
    zap['m'] = make_field(multiplicity)   = 3;
    zap['o'] = make_field(outputFile)     = "anglerfish.root";
    zap['n'] = make_field(numberOfEvents) = JLimit::max();
    zap['d'] = make_field(debug)          = 1;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }

  outputFile.open();

  if (!outputFile.is_open()) FATAL("Error opening file " << outputFile << endl);

  outputFile.put(JMeta(argc, argv));
  
  JDetector detector;

  try {
    load (detectorFile, detector);
  }
  catch(const JException & error) {
    cerr << "FATAL ERROR. Could not open detector file '" << detectorFile << "'." << endl;
    exit(1);
  }

  const floor_range   range = getRangeOfFloors(detector);

  const JStringRouter stringRouter(detector);

  TH1D h1 ("MEvt", "", 50, 0.5, 50.5);

  TH2D h2 ("FEvt", NULL, 
           getNumberOfStrings(detector) + 0, -0.5, getNumberOfStrings(detector) - 0.5,
           range.getLength(), range.getLowerLimit() - 0.5, range.getUpperLimit() + 0.5);
  
  TH1D h3 ("ML0" , "", 50, 0.5, 50.5);

  TH2D h4 ("FL0", NULL, 
           getNumberOfStrings(detector) + 0, -0.5, getNumberOfStrings(detector) - 0.5,
           range.getLength(), range.getLowerLimit() - 0.5, range.getUpperLimit() + 0.5);

  for (Int_t i = 1; i <= h2.GetXaxis()->GetNbins(); ++i) {
    h2.GetXaxis()->SetBinLabel(i, MAKE_CSTRING(stringRouter.at(i-1)));
  }
 
  for (Int_t i = 1; i <= h2.GetXaxis()->GetNbins(); ++i) {
    h4.GetXaxis()->SetBinLabel(i, MAKE_CSTRING(stringRouter.at(i-1)));
  }

  const JDAQHitRouter hitRouter(detector);

  {
    typedef JDAQSnapshotHit JHit_t;
    typedef vector<JHit_t> JHitBuffer;

    for (JMultipleFileScanner<JDAQEvent> in(inputFile, numberOfEvents); in.hasNext(); ){

      JDAQEvent*  event  = in.next();  
      JHitBuffer& buffer = event->getHits<JHit_t>();

      std::sort(buffer.begin(),buffer.end(),less<JDAQKeyHit>());

      for (JHitBuffer::const_iterator p = buffer.begin() ; p != buffer.end() ; ) {

        JHitBuffer::const_iterator q = p;
      
        for (JDAQHit::JTDC_t t1 = p->getT(); 
             (++q              != buffer.end()     && 
              q->getModuleID() == p->getModuleID() &&
              q->getPMT()      == p->getPMT()      &&
              q->getT()        == t1 + deltaT); 
           t1 = q->getT()) {}
        
        int d = distance(p,q);
      
        const JPMTChannel& channel = hitRouter.getPMTChannel(*p);

        if (d >= multiplicity) 
          h2.Fill((double) stringRouter.getIndex(channel.getString()), (double) channel.getFloor());
        
        if (d >= 2) 
          h1.Fill(d);           
        
        p=q;
      }   
    }
  }

  {
    typedef JDAQHit                  JHit_t;
    typedef vector<JHit_t>           JHitBuffer;
    typedef JSuperFrame2D<JHit_t> JSuperFrame2D_t;

    const JModuleRouter moduleRouter(detector);

    for (JMultipleFileScanner<JDAQTimeslice> in(inputFile, numberOfEvents); in.hasNext(); ) {

      JDAQTimeslice* slice = in.next();
    
      for(JDAQTimeslice::iterator frame = slice->begin(); frame != slice->end() ; ++frame){
      
        const JModule&   module = moduleRouter.getModule(frame->getModuleID()); 
        JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*frame, module);
        
        for (JSuperFrame2D_t::iterator pmt = buffer.begin() ; pmt != buffer.end() ; ++pmt) {
          
          for (JHitBuffer::const_iterator p = pmt->begin() ; p != pmt->end() ;) {

            JHitBuffer::const_iterator q = p;
            
            for (JDAQHit::JTDC_t t1 = p->getT(); 
                 (++q              != pmt->end()     && 
                  q->getT()        == t1 + deltaT); 
                 t1 = q->getT()) {}
        
            int d = distance(p,q);

            const JPMTChannel& channel = hitRouter.getPMTChannel(JDAQKeyHit(frame->getModuleID(),*p));

            if (d >= multiplicity) 
              h4.Fill((double) stringRouter.getIndex(channel.getString()), (double) channel.getFloor());
              
            if (d >= 2) 
              h3.Fill(d);               

            p=q;
          }
        }
      }
    }
  }

  outputFile.put(h1);
  outputFile.put(h2);
  outputFile.put(h3);
  outputFile.put(h4);

  outputFile.close();
}