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();
}