JSlewingK40.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <algorithm>
#include <random>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF1.h"
#include "TProfile.h"
#include "TProfile2D.h"

#include "km3net-dataformat/online/JDAQ.hh"
#include "JDAQ/JDAQTimesliceIO.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JMath/JMathToolkit.hh"
#include "JTrigger/JHitR0.hh"
#include "JTrigger/JHitToolkit.hh"
#include "JTools/JRange.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JSupport.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * \file
 *
 * Auxiliary program to determine time slewing from K40 data.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  JMultipleFileScanner<JDAQTimeslice> inputFile;
  JLimit_t& numberOfEvents =          inputFile.getLimit();
  string             outputFile;
  string             detectorFile;
  Double_t           TMax_ns;
  JRange<double>     totVeto_ns;
  JRange<int>        multiplicity;
  JRange<double>     ct;
  bool               slewing;
  int                debug;

  try { 

    JParser<> zap("Auxiliary program to determine time slewing from K40 data.");
    
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)          = "monitor.root";
    zap['a'] = make_field(detectorFile);
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['T'] = make_field(TMax_ns)             = 20.5;                      // [ns]
    zap['t'] = make_field(totVeto_ns)          = JRange<double>(0, 1e4);    // time over threshold of reference hits [ns]
    zap['M'] = make_field(multiplicity)        = JRange<int>(2, 2);
    zap['C'] = make_field(ct)                  = JRange<double>(0.0, 1.0);  // cosine space angle between PMTs
    zap['S'] = make_field(slewing);
    zap['d'] = make_field(debug)               = 1;
      
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  using namespace KM3NETDAQ;

  if (!multiplicity.is_valid())          { FATAL("Invalid multiplicity        " << multiplicity << endl); }
  if ( multiplicity.getLowerLimit() < 2) { FATAL("Invalid multiplicity        " << multiplicity << endl); }
  if (!totVeto_ns  .is_valid())          { FATAL("Invalid time over threshold " << totVeto_ns   << endl); }

  JHit::setSlewing(slewing);

  JDetector detector;

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

  if (detector.empty()) {
    FATAL("Empty detector." << endl);
  }

  const JModuleRouter router(detector);

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

  TH1D       h0("h0", NULL,  41, -TMax_ns, +TMax_ns);
  TH1D       h1("h1", NULL,  50, -0.5, 49.5);
  TProfile2D hx("hx", NULL,  50, -0.5, 49.5, 50, -0.5, 49.5, -2*TMax_ns, +2*TMax_ns);
  TProfile   h2("h2", NULL, 100,  0.5, 100.5, -2*TMax_ns, +2*TMax_ns);

  h0.Sumw2();
  h1.Sumw2();

  random_device rd;
  mt19937 g(rd());

  vector<JHitR0> buffer;

  while (inputFile.hasNext()) {

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

    const JDAQTimeslice* timeslice = inputFile.next();

    for (JDAQTimeslice::const_iterator super_frame = timeslice->begin(); super_frame != timeslice->end(); ++super_frame) {

      if (router.hasModule(super_frame->getModuleID()) && !super_frame->empty()) {

        const JModule& module = router.getModule(super_frame->getModuleID());

        // process data

        buffer.resize(super_frame->size());
          
        vector<JHitR0>::iterator out = buffer.begin();
          
        for (JDAQSuperFrame::const_iterator i = super_frame->begin(); i != super_frame->end(); ++i) {
          *out = JHitR0(i->getPMT(), JHitToolkit<JHit>::getHit(*i, module.getPMT(i->getPMT())));
          ++out;
        }

        sort(buffer.begin(), buffer.end());     
          
        for (vector<JHitR0>::iterator p = buffer.begin(); p != buffer.end(); ) {
            
          vector<JHitR0>::iterator q = p;
            
          while (++q != buffer.end() && q->getT() - p->getT() <= TMax_ns ) {}

          if (multiplicity(distance(p,q))) {

            shuffle(p,q,g);

            for (vector<JHitR0>::const_iterator __p = p; __p != q; ++__p) {
              for (vector<JHitR0>::const_iterator __q = __p; ++__q != q; ) {

                const double dot = JMATH::getDot(module.getPMT(__p->getPMT()),
                                                 module.getPMT(__q->getPMT()));

                if (ct(dot)) {

                  h0.Fill(__q->getT() - __p->getT());
                  h1.Fill(__q->getToT());
                  hx.Fill(__p->getToT(), __q->getToT(), __q->getT() - __p->getT());

                  if (totVeto_ns(__p->getToT())) {
                    h2.Fill(__q->getToT(), __q->getT() - __p->getT());
                  }
                }
              }
            }

            sort(p,q);
          }

          p = q;
        }
      }
    }
  }
  STATUS(endl);


  // Fit function

  TF1 f1("f1", "[0]*exp([1]*sqrt(x) + [2]*x) + [3]");

  f1.SetParameter(0, h2.GetMaximum());
  f1.SetParameter(1, -0.01);
  f1.SetParameter(2, -0.05);
  f1.SetParameter(3, h2.GetMinimum());

  h2.ProjectionX()->Fit(&f1, "", "same");

  for (int i = 0; i != f1.GetNpar(); ++i) {
    cout << "\tstatic double p" << i << "() { return " << setw(9) << setprecision(5) << f1.GetParameter(i) << "; }" << endl;
  }


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