JOMGsim.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <limits>
#include <vector>
#include <algorithm>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TRandom3.h"
 
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
 
#include "km3net-dataformat/online/JDAQ.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JPMTParametersMap.hh"
#include "JDetector/JPMTAnalogueSignalProcessor.hh"
 
#include "JTrigger/JHit.hh"
 
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
 
#include "JTools/JCombinatorics.hh"
#include "JTools/JRouter.hh"
 
#include "JAAnet/JHead.hh"
 
#include "JCalibrate/JCalibrateK40.hh"
 
#include "JLang/JComparator.hh"
#include "JLang/JComparison.hh"
 
#include "JROOT/JRootToolkit.hh"
#include "JROOT/JRootFileWriter.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
namespace {
 
  using JTRIGGER::JHit;
  using JDETECTOR::JPMTPulse;
 
  /**
   * Type definition to convert PMT pulse to hit.
   */
  struct hit_type :
    public JHit
  {
    /**
     * Constructor.
     *
     * \param  pmt           PMT identifier
     * \param  hit           PMT pulse
     */
    hit_type(const size_t pmt, const JPMTPulse& hit) :
      JHit(round(hit.t_ns), round(hit.tot_ns)),                  // minic TDC in CLB
      pmt(pmt)
    {}
 
    size_t pmt;
  };
}
 
/**
 * \file
 *
 * Auxiliary program to determine PMT parameters from OMGSim data. 
 *
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
 
  JMultipleFileScanner<Evt> inputFile;
  JLimit_t&          numberOfEvents = inputFile.getLimit();
  string             outputFile;
  string             detectorFile;
  Double_t           TMax_ns;
  Double_t           WBin_ns;
  double             R;
  JPMTParametersMap  parameters;
  size_t             recycling;
  ULong_t            seed;
  int                debug;
 
  try { 
 
    JParser<> zap("Auxiliary program to determine PMT parameters from OMGsim data.");
    
    zap['f'] = make_field(inputFile,      "input file.");
    zap['o'] = make_field(outputFile,     "output file.")                                = "omgsim.root";
    zap['n'] = make_field(numberOfEvents)                                                = JLimit::max();
    zap['a'] = make_field(detectorFile,   "detector file.");
    zap['T'] = make_field(TMax_ns,        "time window [ns].")                           = 20.0;
    zap['W'] = make_field(WBin_ns,        "bin width [ns].")                             = 1.0;
    zap['R'] = make_field(R,              "radioactivity [Bq]");
    zap['P'] = make_field(parameters,     "PMT parameters")                              = JPARSER::initialised();
    zap['N'] = make_field(recycling,      "recycling of events")                         = 1;
    zap['S'] = make_field(seed,           "seed")                                        = 0;
    zap['d'] = make_field(debug,          "debug flag.")                                 = 1;
     
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  gRandom->SetSeed(seed);
 
 
  JDetector detector;
 
  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
  if (detector.size() != 1u) {
    FATAL("Wrong number of modules " << detector.size() << endl);
  }
 
 
  const JModule& module = detector[0];
 
  JCombinatorics combinatorics(module.size());
 
  combinatorics.sort(JPairwiseComparator(module));
 
  const Int_t    nx   = combinatorics.getNumberOfPairs();
  const Double_t xmin = -0.5;
  const Double_t xmax = nx - 0.5;
 
  const Double_t ymin = -floor(TMax_ns) - 0.5*WBin_ns;
  const Double_t ymax = +floor(TMax_ns) + 0.5*WBin_ns;
  const Int_t    ny   = (Int_t) ((ymax - ymin) / WBin_ns);
 
  TH2D h2r(MAKE_CSTRING(module.getID() << _2R), NULL, nx, xmin, xmax, ny, ymin, ymax);
 
  h2r.Sumw2();
 
  TH1D h1("M", NULL, 101, -0.5, 100.5);
 
  JRouter<int> router;
 
  size_t slewing = module.size();
 
  //size_t slewing = 0;
 
  vector<JPMTAnalogueSignalProcessor> cpu;
 
  for (size_t i = 0; i != module.size(); ++i) {
 
    router.put(module[i].getID(), i);
 
    const JPMTParameters& wip = parameters.getPMTParameters(JPMTIdentifier(module.getID(),i));
 
    if (wip.slewing) {
      //slewing += 1;
    }
 
    cpu.push_back(JPMTAnalogueSignalProcessor(wip));
  }
 
  if      (slewing == 0)
    JHit::setSlewing(false);
  else if (slewing == module.size())
    JHit::setSlewing(true);
  else
    FATAL("Inconsistent PMT slewing options " << slewing << endl);
 
  NOTICE("PMT slewing option " << JHit::getSlewing() << endl);
 
 
  double numberOfPrimaries = 0.0;
 
  try {
 
    STATUS("Extracting header data... " << flush);
 
    const Head  header = getHeader(inputFile);
 
    const JHead buffer(header);
 
    if (buffer.is_valid(&JHead::norma)) 
      numberOfPrimaries = buffer.norma.numberOfPrimaries;
    else
      THROW(JException, "Missing header information.");
 
    STATUS("OK" << endl);
 
  } catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
  if (numberOfPrimaries == 0.0) {
    FATAL("Number of primaries " << numberOfPrimaries << endl); 
  }
  
 
  TFile out(outputFile.c_str(), "recreate");
 
  putObject(out, JMeta(argc, argv));
 
  while (inputFile.hasNext()) {
 
    if (inputFile.getCounter()%10000== 0) {
      STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
    }
 
    const Evt* evt = inputFile.next();
 
    if (evt->mc_hits.size() >= 2u) {
 
      JModuleData input(module.size());
 
      for (const auto& hit : evt->mc_hits) {
 
        if (router.has(hit.pmt_id)) {
 
          const int pmt = router.get(hit.pmt_id);
 
          input[pmt].push_back(JPMTSignal(hit.t, hit.a));
        }
      }
 
      for (size_t i = 0; i != recycling; ++i) {
 
        vector<hit_type> buffer;
 
        for (size_t pmt = 0; pmt != module.size(); ++pmt) {
 
          if (!input[pmt].empty()) {
 
            JPMTData<JPMTPulse> output;
 
            cpu[pmt](module[pmt].getCalibration(), input[pmt], output);
 
            for (const auto& hit : output) {
              buffer.push_back(hit_type(pmt, hit));
            }
          }
        }
 
        if (buffer.size() >= 2u) {
 
          sort(buffer.begin(), buffer.end());
 
          for (vector<hit_type>::iterator p = buffer.begin(); p != buffer.end(); ) {
 
            vector<hit_type>::iterator q = p;
 
            for (++q; q != buffer.end() && q->getT() - p->getT() <= TMax_ns; ++q) {}                           // time window
 
            if (distance(p,q) >= 2) {
 
              sort(p, q, make_comparator(&hit_type::pmt, JComparison::lt()));                                  //
 
              const int M = distance(p, unique(p, q, make_comparator(&hit_type::pmt, JComparison::eq())));     // unique PMTs
 
              if (M >= 2) {
                h1.Fill(M);
              }
 
              sort(p, q);
            }
 
            p = q;
          }
 
          if (debug >= debug_t) {
 
            cout << "hit:";
 
            for (const auto& hit : buffer) {
              cout << ' ' << setw(2) << hit.pmt << ' ' << FIXED(10,1) << hit.getT() << ' ' << FIXED(3,0) << hit.getToT();
            }
 
            cout << endl;
          }
 
          for (vector<hit_type>::const_iterator p = buffer.begin(); p != buffer.end(); ++p) {
            for (vector<hit_type>::const_iterator q = buffer.begin(); p != q; ++q) {
              h2r.Fill((double) combinatorics.getIndex(p->pmt,q->pmt), 
                       JCombinatorics::getSign(p->pmt,q->pmt) * (q->getT() - p->getT()));
            }
          }
        }
      }
    }
  }
  STATUS(endl);
 
  // reset under and overflows.
 
  for (int ix = 0; ix <= h2r.GetXaxis()->GetNbins() + 1; ++ix) {
    h2r.SetBinContent(ix, 0, 0.0);  h2r.SetBinContent(ix, h2r.GetYaxis()->GetNbins() + 1, 0.0);
    h2r.SetBinError  (ix, 0, 0.0);  h2r.SetBinError  (ix, h2r.GetYaxis()->GetNbins() + 1, 0.0);
  }
  
  for (int iy = 0; iy <= h2r.GetYaxis()->GetNbins() + 1; ++iy) {
    h2r.SetBinContent(0, iy, 0.0);  h2r.SetBinContent(h2r.GetXaxis()->GetNbins() + 1, iy, 0.0);
    h2r.SetBinError  (0, iy, 0.0);  h2r.SetBinError  (h2r.GetXaxis()->GetNbins() + 1, iy, 0.0);
  }
  
 
  for (Int_t ix = 1; ix <= h2r.GetXaxis()->GetNbins(); ++ix) {
    for (Int_t iy = 1; iy <= h2r.GetYaxis()->GetNbins(); ++iy) {
      if (h2r.GetBinError(ix,iy) == 0.0) {
        h2r.SetBinError(ix, iy, 1.0);
      }
    }
  }
 
  h2r.Scale(R / (numberOfPrimaries * recycling * WBin_ns));
  h1 .Scale(R / (numberOfPrimaries * recycling));
 
  out << h1 << h2r;
 
  out.Write();
  out.Close();
}