JDAQProfile.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>

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

#include "JDAQ/JDAQSummarysliceIO.hh"
#include "JDAQ/JDAQEventIO.hh"
#include "km3net-dataformat/online/JDAQClock.hh"
#include "JDAQ/JDAQEvaluator.hh"

#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JSupportToolkit.hh"
#include "JSupport/JSupport.hh"

#include "JROOT/JManager.hh"

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


/**
 * \file
 *
 * Example program to histogram various data profiles.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  string         inputFile;
  JLimit_t       numberOfEvents;
  string         outputFile;
  Long64_t       prescale;
  int            debug;

  try { 

    JParser<> zap("Example program to histogram various data profiles.");
    
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)          = "profile.root";
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['P'] = make_field(prescale)            = 1;
    zap['d'] = make_field(debug)               = 1;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }


  NOTICE("Determine frame index range... " << flush);

  JFrameIndexRange frame_index = getFrameIndexRange<JDAQSummaryslice>(inputFile);

  NOTICE("= (" << frame_index.first << "," << frame_index.second << ")" << endl);

  const Long64_t NX   = frame_index.second  -  frame_index.first  +  1;
  const Long64_t nx   = (NX + prescale - 1) / prescale;
  const Double_t xmin = (Double_t) frame_index.first  - 0.5;
  const Double_t xmax = (Double_t) frame_index.second + 0.5;


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

  typedef JManager<int, TH1D>  JManager_t;

  JManager_t map_summary(new TH1D("Summary[%]", NULL, NX, xmin, xmax));
  JManager_t map_trigger(new TH1D("Trigger[%]", NULL, nx, xmin, xmax));
  JManager_t map_status (new TH1D("Status[%]",  NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5));
 
  TH1D hn("N",       NULL, NX, xmin, xmax);
  TH1D h0("L0",      NULL, NX, xmin, xmax);
  TH1D h1("WR",      NULL, NX, xmin, xmax);
  TH1D h2("trigger", NULL, nx, xmin, xmax);

  h2.Sumw2();

  for (JMultipleFileScanner<JDAQSummaryslice> in(inputFile, numberOfEvents); in.hasNext(); ) {

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

    JDAQSummaryslice* summaryslice = in.next();

    const Double_t x = summaryslice->getFrameIndex();


    for (JDAQSummaryslice::const_iterator frame = summaryslice->begin(); frame != summaryslice->end(); ++frame) {

      if (frame->testHighRateVeto()) {

        for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {

          if (frame->testHighRateVeto(pmt)) {
            map_status[frame->getModuleID()]->Fill((Double_t) pmt);
          }
        }
      }

      Double_t y = 0.0;

      for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
        if (!frame->testHighRateVeto(pmt)) {
          y += frame->getRate(pmt) * 1e-3;  // kHz
        }
      }

      map_summary[frame->getModuleID()]->Fill(x, y);
    }


    Double_t y = 0.0;

    for (JDAQSummaryslice::const_iterator frame = summaryslice->begin(); frame != summaryslice->end(); ++frame) {

      for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
        y += frame->getRate(pmt) * 1e-3;  // kHz
      }

      h1.Fill(x, (frame->testWhiteRabbitStatus() ? 1.0 : 0.0));
    }

    hn.Fill(x, summaryslice->size());
    h0.Fill(x, y);
  }
  STATUS(endl);


  for (JMultipleFileScanner<JDAQEvent> in(inputFile, numberOfEvents); in.hasNext(); ) {

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

    JDAQEvent* event = in.next();

    const Double_t x = event->getFrameIndex();
    const Double_t y = 1.0e9 / getFrameTime() / prescale;

    h2.Fill(x, y);

    typedef JDAQTriggeredHit JHit_t;
    //typedef JDAQSnapshotHit JHit_t;

    for (JDAQEvent::const_iterator<JHit_t> hit = event->begin<JHit_t>(); hit != event->end<JHit_t>(); ++hit) {
      map_trigger[hit->getModuleID()]->Fill(x);
    }
  }
  STATUS(endl);

  map_summary.Write(out);
  map_trigger.Write(out);
  map_status .Write(out);

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