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;
  UInt_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();
 
    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 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.125;
  const Double_t ymax = +floor(TMax_ns) - 0.125;
  const Int_t    ny   = (Int_t) ((ymax - ymin) / 0.25);

  TH2D h2s(MAKE_CSTRING(module.getID() << _2S), NULL, nx, xmin, xmax, ny, ymin, ymax);

  h2s.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); 
  }


  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) {
          h2s.Fill((double) combinatorics.getIndex(p->pmt,q->pmt), 
                   JCombinatorics::getSign(p->pmt,q->pmt) * (q->t1 - p->t1));
        }
      }
    }
  }
  STATUS(endl);


  const double W = R / numberOfPrimaries;

  h2s.Scale(W);

  out << h2s;

  out.Write();
  out.Close();
}