JAnglerFish.cc File Reference

Example application to search for sets of consecutive hits on a given PMT,
that are separated in time by a constant value. More...

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

Go to the source code of this file.

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

Detailed Description

Example application to search for sets of consecutive hits on a given PMT,
that are separated in time by a constant value.

The number of consecutive hits
as well as the time difference are selected from the command line.
The searches are done within the JDAQSnapshotHits from JDAQEvents\ as well as within the JDAQHits from the JDAQTimeslice data.

Author

rgruiz

Definition in file JAnglerFish.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 41 of file JAnglerFish.cc.

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