JOMGsim.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <limits>
 
#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/JPMTParameters.hh"
#include "JDetector/JPMTAnalogueSignalProcessor.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 "JROOT/JRootToolkit.hh"
#include "JROOT/JRootFileWriter.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \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             R;
  JPMTParameters     parameters;
  double             QE = 1.0;
  ULong_t            seed;
  int                debug;
 
  parameters.TTS_ns          = 0.0;
  parameters.gainSpread      = 0.0;
  parameters.PunderAmplified = 0.0;
  parameters.slewing         = false;    
 
  try { 
 
    JProperties properties = parameters.getProperties();
 
    properties["QE"] = QE;
 
    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['R'] = make_field(R,              "radioactivity [Bq]");
    zap['P'] = make_field(properties,     "PMT parameters")                              = JPARSER::initialised();
    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 Double_t WS   = 0.25;
 
  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*WS;
  const Double_t ymax = +floor(TMax_ns) - 0.5*WS;
  const Int_t    ny   = (Int_t) ((ymax - ymin) / WS);
 
  TH2D h2r(MAKE_CSTRING(module.getID() << _2R), NULL, nx, xmin, xmax, ny, ymin, ymax);
 
  h2r.Sumw2();
 
 
  JRouter<int> router;
 
  for (size_t i = 0; i != module.size(); ++i) {
    router.put(module[i].getID(), i);
  }
 
  double numberOfPrimaries = 0.0;
 
  Head header;
 
  try {
 
    STATUS("Extracting header data... " << flush);
 
    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); 
  }
 
 
  parameters.QE *= QE;
 
  const JPMTAnalogueSignalProcessor cpu(parameters);
  
  JPMTData<JPMTSignal> input;
  JPMTData<JPMTPulse>  output;
 
  struct hit_type {
    int    pmt;
    double t1;
  };
  
  vector<hit_type> buffer;
 
 
  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) {
 
      buffer.clear();
 
      for (const auto& hit : evt->mc_hits) {
 
        if (router.has(hit.pmt_id)) {
 
          const int pmt = router.get(hit.pmt_id);
 
          input .clear();
          output.clear();
 
          input.push_back(JPMTSignal(hit.t, hit.a));
 
          cpu(module[pmt].getCalibration(), input, output);
 
          for (const auto& hit : output) {
            buffer.push_back({pmt, hit.t_ns});
          }
        }
      }
 
 
      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->t1 - p->t1));
        }
      }
    }
  }
  STATUS(endl);
 
 
  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);
      }
    }
  }
 
  const double W = R / (WS * numberOfPrimaries);
 
  h2r.Scale(W);
 
  out << h2r;
 
  out.Write();
  out.Close();
}