JSupernovaL0.cc

Go to the documentation of this file.

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

#include "TROOT.h"
#include "TFile.h"
#include "TH2D.h"

#include "JGizmo/JManager.hh"

#include "JDAQ/JDAQ.hh"
#include "JDAQ/JDAQTimeslice.hh"
#include "JDAQ/JDAQSummaryslice.hh"
#include "JDAQ/JDAQEvaluator.hh"
#include "JDAQ/JDAQSummaryFrame.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"

#include "JTools/JRange.hh"

#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JSupport.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

/**
 * \author mlincett
 * \file 
 *
 * Preliminary study of summary rates for supernova detection
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  JMultipleFileScanner<JDAQSummaryslice> inputFile;
  string             outputFile;
  int                debug;
  int                nLines;

  try { 

    JParser<> zap("Monitoring of PMT time over threshold distributions.");
    
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)        = "summaryrates.root";
    zap['d'] = make_field(debug)             = 1;
    zap['n'] = make_field(nLines)            = 2;
    zap(argc, argv);

  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }

  JTreeScanner<JDAQSummaryslice, JDAQEvaluator> in(inputFile);

  cout.tie(&cerr);

  using namespace KM3NETDAQ;

  const double factor = 1.0/1000 ;
  const int Nbins = 90 ; // reasonable range for singles rates from 0-20 kHz

  set<int> modules;


  const double maxRate_kHz = 20;

  JManager<int, TH2D> SinglesRatedistr(new TH2D("%.2S", NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS-0.5, Nbins, JDAQRate::getData(factor))) ;

  JManager<int, TH1D> PMT(new TH1D("pmt_S%", NULL, Nbins, JDAQRate::getData(factor)));
  JManager<int, TH1D> DOM(new TH1D("dom_S%", NULL, 1000, 0, maxRate_kHz * NUMBER_OF_PMTS ));
  JManager<int, TH1D> DET(new TH1D("det_S%", NULL, 1000, 0, nLines * 18 * maxRate_kHz * NUMBER_OF_PMTS));

  JManager<int, TH1D> Delta(new TH1D("delta_S%", NULL, 2000, -5000, +5000));
  
  JManager<int, TH1D> History(new TH1D("history_S%", NULL, 300000, 0, 300000)); 
  
  const int nStages = 3;

  vector<double> deltaRate(nStages, 0);
  vector<double> detRateBuf(nStages, 0);

  while (in.hasNext()) {

    vector<double> detRate(nStages, 0);

    int counter = in.getCounter();

    STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl);

    JDAQSummaryslice* summary = in.next();

    JDAQSummaryslice* previous_summary;

    int nDetectorActivePMTs = 0;

    int nDetectorTotalPMTs = 0;

    // loop over modules

    for (KM3NETDAQ::JDAQSummaryslice::const_iterator summary_frame = summary->begin(); summary_frame != summary->end(); ++summary_frame) {
      
      vector<double> domRate(nStages, 0); 
     
      int DOMID = summary_frame->getModuleID() ;
      
      TH2D* h_i = SinglesRatedistr[DOMID];

      int nModuleActivePMTs = 0;

      for (int ipmt = 0 ; ipmt < NUMBER_OF_PMTS ; ipmt++) {     

        double pmtRate = summary_frame->getRate(ipmt, factor);

        h_i->Fill(ipmt, pmtRate, summary_frame->getWeight(ipmt, factor));
        
        PMT[0]->Fill(pmtRate);

        domRate[0] += pmtRate;

        nDetectorTotalPMTs++;

        if ( !( summary_frame->testHighRateVeto(ipmt) || summary_frame->testFIFOStatus(ipmt) ) ) {

          nModuleActivePMTs++;

          PMT[1]->Fill(pmtRate);

          domRate[1] += pmtRate;

          /*
          if (counter != 0) {

            KM3NETDAQ::JDAQSummaryslice::const_iterator previous_frame = previous_summary->begin();

            for ( ; previous_frame != previous_summary->end(); ++previous_frame) {
              if (DOMID == previous_frame->getModuleID()) {
                break;
              }
            }

            if (!( previous_frame->testHighRateVeto(ipmt) || previous_frame->testFIFOStatus(ipmt) ) ) {
              PMT[2]->Fill(pmtRate);
              domRate[2] += pmtRate;
            }
          }     
        
          */      
        }
      }


      // DOM summary

      nDetectorActivePMTs += nModuleActivePMTs;
      nDetectorTotalPMTs += NUMBER_OF_PMTS;

      for (int i = 0; i < nStages; i++) {
        if (domRate[i] > 0) {
          detRate[i] += domRate[i];
          if (i == 1) {
            domRate[i] /= ((double)nModuleActivePMTs) / NUMBER_OF_PMTS;
          }
          DOM[i]->Fill(domRate[i]);

          if (i == 0) { modules.insert(DOMID); }
        }
      }
    }

    for (int i = 0; i < nStages; i++) {

      DET[i]->Fill(detRate[i]);
      if (i == 1) {
        detRate[i] /= ((double)nDetectorActivePMTs)/((double)nDetectorTotalPMTs);
      }
      History[i]->Fill(counter, detRate[i]);

      if (counter != 0) {
        double deltaRate = detRate[i] - detRateBuf[i];
        Delta[i]->Fill(deltaRate);
      }

      detRateBuf[i] = detRate[i];

    }

    previous_summary = summary;
    
  }
  
  STATUS(endl);

  STATUS(modules.size() << " modules active in the run." << endl);

  // Store histogram(s).
  
  TFile out(outputFile.c_str(), "recreate");
  
  PMT.Write(out);
  DOM.Write(out);
  DET.Write(out);
  Delta.Write(out);
  History.Write(out);
  
  
  // SinglesRatedistr.Write(out) ;
  
  out.Write();
  out.Close();

}