JEqualizer.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 "TProfile.h"
 
#include "JDAQ/JDAQSummarysliceIO.hh"
#include "JDAQ/JSupport.hh"
 
#include "JROOT/JManager.hh"
 
#include "JDetector/JDetectorAddressMap.hh"
#include "JDetector/JDetectorSupportkit.hh"
#include "JDetector/JPMTParametersMap.hh"
#include "JDetector/JPMTParametersToolkit.hh"
 
#include "JGizmo/JGizmoToolkit.hh"
 
#include "JSupport/JMultipleFileScanner.hh"
 
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 *
 * Example program to equalize QE of PMTs based on summary data.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
 
  JMultipleFileScanner<JDAQSummaryslice> inputFile;
  JLimit_t&            numberOfEvents = inputFile.getLimit();
  string               outputFile;
  string               pmtFile;
  double               rate_kHz;
  string               rings;
  int                  debug;
 
  try { 
 
    JParser<> zap("Example program to equalize QE of PMTs based on summary data.");
    
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)          = "summaryslice.root";
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['P'] = make_field(pmtFile,           "specify PMT file name that can be given to JTriggerEfficiency.") = "";
    zap['R'] = make_field(rate_kHz,          "specify expected singles rate [Hz]")                             = 0.0;
    zap['r'] = make_field(rings,             "rings, e.g. \"AB\"")                                             = "ABCDEF";
    zap['d'] = make_field(debug)               = 2;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
 
  const double factor = 1.0e-3;   // [kHz]
 
  rings = to_upper(rings);
 
 
  JPMTParametersMap parameters;
 
  if (pmtFile != "") {
    parameters.load(pmtFile.c_str());
  }
 
 
  map<JPMTIdentifier, double> QE;
 
  for (JPMTParametersMap::const_iterator i = parameters.begin(); i != parameters.end(); ++i) {
    QE[i->first] = getSurvivalProbability(i->second) * i->second.QE;
  }
 
 
  JManager<int, TH1D>     H0(new TH1D    ("H0[%]", NULL, JDAQRate::getN(), JDAQRate::getData(factor)));
  JManager<int, TProfile> H1(new TProfile("H1[%]", NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5));
 
  int detector = -1;
 
  while (inputFile.hasNext()) {
    
    STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
    
    const JDAQSummaryslice* summary = inputFile.next();
 
    detector = summary->getDetectorID();
 
    const JDetectorAddressMap& demo = getDetectorAddressMap(summary->getDetectorID());
 
    for (JDAQSummaryslice::const_iterator frame = summary->begin(); frame != summary->end(); ++frame) {
 
      const JModuleAddressMap& memo = demo.get(frame->getModuleID());
 
      TH1D*     h0 = H0[frame->getModuleID()];
      TProfile* h1 = H1[frame->getModuleID()];
 
      for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
 
        if (rings.find(memo.getPMTPhysicalAddress(pmt).ring) != string::npos) {
 
          const JPMTIdentifier id(frame->getModuleID(), pmt);
 
          map<JPMTIdentifier, double>::const_iterator p = QE.find(id);
 
          double R = frame->getRate  (pmt, factor);
          double W = frame->getWeight(pmt, factor);
          double P = 1.0;
 
          if (p != QE.end()) {
            P = p->second;
          }
 
          if (P > 0.0) {
 
            h0->Fill(R / P, W);
            H0->Fill(R / P, W);
            
            h1->Fill(pmt, R / P);
            H1->Fill(pmt, R / P);
          }
        }
      }
    }
  }
  STATUS(endl);
 
 
  if (pmtFile != "" && rate_kHz > 0.0) {
 
    NOTICE("Write QEs in PMT file " << pmtFile << endl);
 
    for (JManager<int, TProfile>::const_iterator i = H1.begin(); i != H1.end(); ++i) {
 
      for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
 
        const JPMTIdentifier id(i->first, pmt);
 
        parameters[id].QE *= i->second->GetBinContent(pmt + 1)  /  rate_kHz;
      }
    }
 
    parameters.store(pmtFile.c_str());
  }
 
 
  const JDetectorAddressMap& demo = getDetectorAddressMap(detector);
 
  setAxisLabels(*H1, "X", demo.getDefaultModuleAddressMap());
 
  for (JManager<int, TProfile>::const_iterator i = H1.begin(); i != H1.end(); ++i) {
    setAxisLabels(*(i->second), "X", demo.get(i->first));
  }
 
 
  TFile out(outputFile.c_str(), "recreate");
 
  H0.Write(out, true);
  H1.Write(out, true);
 
  out.Write();
  out.Close();
}