JCoelacanth.cc

Go to the documentation of this file.

 
#include <string>
#include <iostream>
#include <iomanip>
#include <string>
#include <map>
 
#include "TROOT.h"
#include "TFile.h"
 
#include "JDAQ/JDAQEventIO.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
 
#include "JSupport/JSingleFileScanner.hh"
#include "JSupport/JSupport.hh"
 
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 *
 * Program to correct detector file for inconsistencies in configuration of data filter.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
 
  vector<string> inputFile;
  JLimit_t     numberOfEvents;
  string       outputFile;
  string       detectorFile;
  double       stdev;
  bool         overwriteDetector;
  int          qaqc;
  int          debug;
 
  try { 
 
    JParser<> zap("Program to correct detector file for inconsistencies in configuration of data filter.");
    
    zap['f'] = make_field(inputFile,         "original raw data file and reprocessed raw data file");
    zap['a'] = make_field(detectorFile);
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['S'] = make_field(stdev,             "number of standard deviations for comparison of counts per optical module.")  =  10.0;
    zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a'.");
    zap['Q'] = make_field(qaqc)                =  0;
    zap['d'] = make_field(debug)               = 2;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
 
  if (inputFile.size() != 2) {
    FATAL("Wrong number of input files " << inputFile.size() << " != " << 2 << endl);
  }
  
  JDetector detector;
 
  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
  const JModuleRouter router(detector);
 
  typedef map<int, size_t>  map_type;
 
  map_type zmap[2];
 
  for (int i = 0; i != 2; ++i) {
 
    STATUS("Processing " << inputFile[i] << endl);
 
    for (JSingleFileScanner<JDAQEvent> in(inputFile[i], numberOfEvents); in.hasNext(); ) {
 
      STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
 
      const JDAQEvent* event = in.next();
 
      for (JDAQEvent::const_iterator<JDAQTriggeredHit> hit = event->begin<JDAQTriggeredHit>(); hit != event->end<JDAQTriggeredHit>(); ++hit) {
        zmap[i][hit->getModuleID()] += 1;
      }
    }
    STATUS(endl);
  }
 
  size_t ns = 0;
 
  for (map_type::const_iterator
         p = zmap[0].begin(),
         q = zmap[1].begin(); (p != zmap[0].end() &&
                               q != zmap[1].end()); ) {
 
    if      (p->first < q->first) { ++p; }
    else if (q->first < p->first) { ++q; }
    else {
 
      if (max ((double) p->second,  (double) q->second)  >  stdev * sqrt(p->second + q->second)  &&
          fabs((double) p->second - (double) q->second)  >  stdev * sqrt(p->second + q->second)) {
 
        ns += 1;
 
        NOTICE("Module " << setw(8) << p->first << ' ' << setw(6) << p->second << " <> " << setw(6) << q->second << ' ' << (overwriteDetector ? "-> set position to (0,0,0)" : "") << endl);
 
        if (overwriteDetector) {
          detector[router.getIndex(p->first)].set(JVector3D(0.0,0.0,0.0));
        }
      }
 
      ++p;
      ++q;
    }
  }
 
  if (overwriteDetector) {
 
    NOTICE("Store calibration data on file " << detectorFile << endl);
 
    store(detectorFile, detector);
  }
 
  QAQC(ns << endl);
 
  return 0;
}