examples/JSirene/JDomino.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <map>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"

#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh" 

#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JPMTRouter.hh"

#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"

#include "JLang/JPredicate.hh"
#include "JTools/JRange.hh"

#include "JROOT/JManager.hh"
#include "JROOT/JRootToolkit.hh"
#include "JGizmo/JGizmoToolkit.hh"

#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * \file
 *
 * Example program to verify Monte Carlo data.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  typedef JRange<size_t>  JRange_t;
  
  JMultipleFileScanner<Evt> inputFile;
  JLimit_t&      numberOfEvents = inputFile.getLimit();
  string         outputFile;
  string         detectorFile;
  JRange_t       M;
  unsigned int   L_cut;
  int            debug;

  try { 

    JParser<> zap("Example program to verify Monte Carlo data.");
    
    zap['f'] = make_field(inputFile);
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['o'] = make_field(outputFile)          = "domino.root";
    zap['a'] = make_field(detectorFile);
    zap['M'] = make_field(M)                   = JRange_t();
    zap['m'] = make_field(L_cut)               = 1;
    zap['d'] = make_field(debug)               = 2;
    
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  JDetector detector;

  if (detectorFile != "") {

    try {
      load(detectorFile, detector);
    }
    catch(const JException& error) {
      FATAL(error);
    }
  }

  const JPMTRouter router(detector);

  const JHead header = getHeader(inputFile);

  JPosition3D center = get<JPosition3D>(header);
  
  NOTICE("Apply detector offset " << center << endl); 

  detector -= center;

  vector<Double_t> X;

  double x = -100.0;

  for ( ; x <  -10.0; x +=  5.0) { X.push_back(x); }
  for ( ; x <  +20.0; x +=  1.0) { X.push_back(x); }
  for ( ; x <  +50.0; x +=  2.0) { X.push_back(x); }
  for ( ; x < +100.0; x +=  5.0) { X.push_back(x); }
  for ( ; x < +250.0; x += 10.0) { X.push_back(x); }
  for ( ; x < +500.0; x += 25.0) { X.push_back(x); }
  for ( ; x < +900.0; x += 50.0) { X.push_back(x); }

  JManager<int, TH1D> H1(new TH1D("h[%]", NULL, X.size() - 1, X.data()));

  map<int, int> npe;

  size_t number_of_events = 0;

  while (inputFile.hasNext()) {

    STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);

    const Evt* event = inputFile.next();

    if (!event->mc_hits.empty()) {
      ++number_of_events;
    }

    size_t number_of_muons = 0;
    
    for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
      if (is_muon(*track)) {
        number_of_muons += 1;
      }
    }

    if (!M(number_of_muons)) {
      continue;
    }

    if (event->mc_hits.size() < L_cut) continue; //simple multiplicity cut
    
    for (vector<Hit>::const_iterator hit = event->mc_hits.begin(); hit != event->mc_hits.end(); ++hit) {

      npe[hit->type] += hit->a;

      const double t1 = getTime(*hit);
              
      vector<Trk>::const_iterator track = find_if(event->mc_trks.begin(),
                                                  event->mc_trks.end(),
                                                  make_predicate(&Trk::id, hit->origin));

      if (track != event->mc_trks.end() && is_muon(*track)) {
      
        const JTrack3E muon = getTrack(*track);

        if (router.hasPMT(hit->pmt_id)) {
              
          const JPMT&  pmt = router.getPMT(hit->pmt_id);
          const double t0  = muon.getT(pmt);

          H1[pmt.getID()]->Fill(t1 - t0, hit->a);

          H1->Fill(t1 - t0, hit->a);
        }       
      }
    }
  }
  STATUS(endl); 

  if (number_of_events != 0) {

    const double W = 1.0 / (double) number_of_events;

    convertToPDF(*H1, "WE", W);

    for (auto& i : H1) {
      convertToPDF(*i.second, "WE", W);
    }
  }

  if (debug >= status_t) {

    cout << "photon-electron statistics" << endl;

    for (const auto& i : npe) {
      cout << setw(3) << i.first << ' ' << setw(6) << i.second << endl;
    }
  }

  TFile out(outputFile.c_str(), "recreate");

  out << *H1 << H1;

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