#include <string>
#include <iostream>
#include <iomanip>
#include <limits>
#include <array>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "km3net-dataformat/online/JDAQ.hh"
#include "km3net-dataformat/online/JDAQHeader.hh"
#include "JDAQ/JDAQTimesliceIO.hh"
#include "JDAQ/JDAQSummarysliceIO.hh"
#include "JDAQ/JDAQEvaluator.hh"
#include "JTrigger/JSuperFrame1D.hh"
#include "JTrigger/JSuperFrame2D.hh"
#include "JTrigger/JHitR0.hh"
#include "JTrigger/JMatchL0.hh"
#include "JTrigger/JHitToolkit.hh"
#include "JTrigger/JTriggerParameters.hh"
#include "JTrigger/JTriggerToolkit.hh"
#include "JTrigger/JPreprocessor.hh"
#include "JTrigger/JTriggerToolkit.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JPMTAnalogueSignalProcessor.hh"
#include "JSupport/JSingleFileScanner.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JSummaryFileRouter.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
#include "JSupport/JTriggerParametersSupportkit.hh"
#include "JTools/JRange.hh"
#include "JTools/JAbstractHistogram.hh"
#include "JROOT/JROOTClassSelector.hh"
#include "JROOT/JRootToolkit.hh"
#include "JROOT/JRootFileWriter.hh"
#include "JLang/JObjectMultiplexer.hh"
#include "JCalibrate/JCalibrateK40.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
namespace {
using JLANG::JObjectID;
using JTOOLS::JCombinatorics;
using JTOOLS::JAbstractHistogram;
using JDETECTOR::JModule;
// Methods for background estimation
static const char* const rates_t = "rates"; //!< singles rates from summary data
static const char* const count_t = "counts"; //!< singles rates from L0 data
static const char* const tails_t = "tails"; //!< tails of time residual distribution
static const char* const rndms_t = "randoms"; //!< random time offsets (only L0 data)
/**
* Maximal time-over-threshold.
*/
static const double ToTmax_ns = 250;
/**
* Auxiliary data structure for histogram management.
*/
struct JHistogram :
public JObjectID,
public JCombinatorics
{
/**
* Default constructor.
*/
JHistogram() :
h2s(NULL),
h1b(NULL),
h1L(NULL)
{}
/**
* Constructor.
*
* \param module module
* \param axis axis
*/
JHistogram(const JModule& module, const JAbstractHistogram<Double_t>& axis) :
JObjectID(module.getID())
{
using namespace JPP;
// sort the PMT pairs by opening angle in the order largest angle -> smallest angle
this->configure(module.size());
this->sort(JPairwiseComparator(module));
const Int_t nx = this->getNumberOfPairs();
const Double_t xmin = -0.5;
const Double_t xmax = nx - 0.5;
h2s = new TH2D(MAKE_CSTRING(module.getID() << _2S), NULL, nx, xmin, xmax, axis.getNumberOfBins(), axis.getLowerLimit(), axis.getUpperLimit());
h1b = new TH1D(MAKE_CSTRING(module.getID() << _1B), NULL, nx, xmin, xmax);
h1L = new TH1D(MAKE_CSTRING(module.getID() << _1L), NULL, nx, xmin, xmax);
h2s->Sumw2();
h1b->Sumw2();
h1L->Sumw2();
}
/**
* Check validity.
*
* \return true if valid; else false
*/
bool is_valid() const
{
return (h2s != NULL &&
h1b != NULL &&
h1L != NULL);
}
TH2D* h2s; //!< time distribution per PMT pair
TH1D* h1b; //!< background rates per PMT pair
TH1D* h1L; //!< livetimes per PMT pair
};
/**
* Get coincidence probability of two PMTs within one module due to random background.
*
* \param E1 expected counts of 1st PMT
* \param E2 expected counts of 2nd PMT
* \param ED expected counts of module, excluding PMT 1 and PMT 2
* \param M_min multiplicity lower limit
* \param M_max multiplicity upper limit
* \return probability
*/
inline double getP(const double E1, const double E2, const double ED, const int M_min, const int M_max)
{
// 1st, always calculate the probability for twofold multiplicity
double P = E1 * E2 * exp(-(E1 + E2 + ED));
// 2nd, calculate the probability of multiplicity at lower limit
int MD = 1;
for ( ; MD + 2 <= M_min; ++MD) {
P *= ED/MD;
}
// 3rd, add the probabilities of multiplicities up to upper limit
double P_sum = P;
for ( ; MD + 2 <= M_max; ++MD) {
P *= ED/MD;
P_sum += P;
}
return P_sum;
}
}
/**
* \file
*
* Auxiliary program to determine PMT parameters from K40 data.
*
* By default, the combinatorial background is estimated from the singles rate.\n
* In case L0 data are taken (and no summary data are available),
* the singles rate can be determined from the amount of L0 data.\n
* One can also estimate the background from the tails in the time residual distribution.\n
* In that case, option -B can be used to specify the minimal and maximal time offset in ns.\n
* For example -B "10 20".\n
* Note that the minimal and maximal time should be larger that the width of
* the time residual distribution and less than the time window of the L1 coincidence, respectively.\n
* The time window is applied symmetrically around a time offset of zero.\n
* Finally, if L0 data are available, one can estimate the background using
* the same procedure but with a random (read wrong) time calibration.
*
* The option -t can be used to specify three time-over-threshold values.\n
* The time-over-threshold of each hit in the coincidence should be greater or equal the first value and less or equal the last value.\n
* Furthermore, the time-over-threshold of one of the two hits should be greater or equals the second value.\n
* If you change the (default) values, you may also want to change the corresponding values in JEditPMTParameters.cc.
*
* The option -M can be used to specify a range of multiplicities.\n
* For example -M "2 2" will strictly select two-fold coincidences.
* Note that such a selection can bias the result.
*
* Also consult <a href="https://common.pages.km3net.de/jpp/JMonitor.PDF">documentation</a>.
* \author mdejong
*/
int main(int argc, char **argv)
{
using namespace std;
using namespace JPP;
using namespace KM3NETDAQ;
JMultipleFileScanner_t inputFile;
counter_type numberOfEvents;
string outputFile;
string detectorFile;
Double_t TMax_ns;
JRange<double> rateRange_Hz;
array <double, 3> totRange_ns;
JRange<double> Tail_ns;
JRange<int> multiplicity;
double deadTime_us;
JROOTClassSelector selector;
string background;
JPreprocessor option;
int debug;
try {
JParser<> zap("Auxiliary program to determine PMT parameters from K40 data.");
zap['f'] = make_field(inputFile, "input file.");
zap['o'] = make_field(outputFile, "output file.") = "calibrate_k40.root";
zap['n'] = make_field(numberOfEvents) = JLimit::max();
zap['a'] = make_field(detectorFile, "detector file.");
zap['T'] = make_field(TMax_ns, "time window [ns].") = 20.0;
zap['V'] = make_field(rateRange_Hz, "PMT rate range [Hz].") = JRange<double>(0.0, 20.0e3);
zap['t'] = make_field(totRange_ns, "PMT time-over-threshold range [ns].") = { 4.0, 7.0, ToTmax_ns };
zap['b'] = make_field(background, "background estimation method.") = rates_t, count_t, tails_t, rndms_t;
zap['B'] = make_field(Tail_ns, "time window used for background estimation.") = JRange<double>(15.0, 20.0);
zap['M'] = make_field(multiplicity, "multiplicity range of hits on DOM.") = JRange<int>(2, 31);
zap['D'] = make_field(deadTime_us, "L1 dead time (us)") = 0.0;
zap['C'] = make_field(selector, "timeslice selector, e.g. JDAQTimesliceL1.") = getROOTClassSelection<JDAQTimesliceTypes_t>();
zap['O'] = make_field(option, "hit pre-processing option.") = JPreprocessor::getOptions(JPreprocessor::filter_t);
zap['d'] = make_field(debug, "debug flag.") = 1;
zap(argc, argv);
}
catch(const exception &error) {
FATAL(error.what() << endl);
}
//-----------------------------------------------------------
// check the input parameters
//-----------------------------------------------------------
if (selector == JDAQTimesliceL2::Class_Name() ||
selector == JDAQTimesliceSN::Class_Name()) {
FATAL("Option -C <selector> " << selector << " not compatible with calibration method." << endl);
}
if (selector == JDAQTimeslice ::Class_Name() ||
selector == JDAQTimesliceL1::Class_Name()) {
JTriggerParameters parameters;
try {
parameters = getTriggerParameters(inputFile);
}
catch(const JException& error) {
FATAL("No trigger parameters from input:" << error.what() << endl);
}
if ((selector == JDAQTimeslice ::Class_Name() && parameters.writeL1.prescale > 0) ||
(selector == JDAQTimesliceL1::Class_Name())) {
if (parameters.TMaxLocal_ns < TMax_ns) {
FATAL("Option -T <TMax_ns> = " << TMax_ns << " is larger than in the trigger " << parameters.TMaxLocal_ns << endl);
}
if (background == rndms_t ||
background == count_t) {
FATAL("Option -C <selector> " << selector << " incompatible with option -b <background> " << background << endl);
}
}
}
if (!multiplicity.is_valid() || multiplicity.getLowerLimit() < 2 || multiplicity.getUpperLimit() > NUMBER_OF_PMTS) {
FATAL("Invalid option -M <multiplicity> " << multiplicity << endl);
}
if (totRange_ns[0] > totRange_ns[1] ||
totRange_ns[1] > totRange_ns[2]) {
FATAL("Invalid option -t <totRange_ns> " << JEEPZ() << totRange_ns << endl);
}
if (background == tails_t) {
if (!Tail_ns.is_valid()) {
FATAL("Invalid option -B <Tail_ns> " << Tail_ns << endl);
}
if (Tail_ns.getUpperLimit() > TMax_ns) {
FATAL("Invalid option -B <Tail_ns> " << Tail_ns << "; upper limit larger than option -T <TMax_ns> " << TMax_ns << endl);
}
}
//-----------------------------------------------------------
// load detector file
//-----------------------------------------------------------
JDetector detector;
try {
load(detectorFile, detector);
}
catch(const JException& error) {
FATAL(error);
}
if (detector.empty()) {
FATAL("Empty detector." << endl);
}
const JModuleRouter router(detector);
//-----------------------------------------------------------
// determine time offsets for background method rndms_t
//-----------------------------------------------------------
vector<double> t0; // time offsets for random background evaluation [ns]
for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
t0.push_back(i * 2 * TMax_ns);
}
//-----------------------------------------------------------
// correction factor for rate due to specified time-over-threshold range
//-----------------------------------------------------------
const JPMTAnalogueSignalProcessor cpu;
const int NPE = 1;
const double RTU = (cpu.getIntegralOfChargeProbability(cpu.getNPE(totRange_ns[0]), cpu.getNPE(totRange_ns[2]), NPE)
/
cpu.getIntegralOfChargeProbability(cpu.getNPE(0.0), cpu.getNPE(ToTmax_ns), NPE));
//-----------------------------------------------------------
// initialise histograms
//-----------------------------------------------------------
vector<JHistogram> zmap(detector.size());
const Double_t ymin = -floor(TMax_ns) + 0.5;
const Double_t ymax = +floor(TMax_ns) - 0.5;
const Int_t ny = (Int_t) ((ymax - ymin) / 1.0);
for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
const JModuleAddress& address = router.getAddress(module->getID());
if (!module->empty()) {
NOTICE("Booking histograms for module " << module->getID() << endl);
zmap[address.first] = JHistogram(*module, JAbstractHistogram<Double_t>(ny, ymin, ymax));
}
}
typedef JHitR0 hit_type;
typedef JSuperFrame2D<hit_type> JSuperFrame2D_t;
typedef JSuperFrame1D<hit_type> JSuperFrame1D_t;
const JMatchL0<hit_type> match(TMax_ns); // time window self-coincidences [ns]
const double deadTime_ns = deadTime_us * 1e3;
TFile out(outputFile.c_str(), "recreate");
putObject(out, JMeta(argc, argv));
counter_type counter = 0;
for (JMultipleFileScanner_t::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
STATUS("Processing: " << *i << endl);
JSummaryFileRouter summary(*i);
JSingleFileScanner<JDAQTimesliceTypes_t> reader(*i);
JObjectMultiplexer<JDAQTimesliceTypes_t, JDAQTimeslice> in(reader, selector);
for (JDAQHeader header; in.hasNext() && counter != numberOfEvents; ++counter) {
STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl);
const JDAQTimeslice* timeslice = in.next();
if (header.getDetectorID() != timeslice->getDetectorID() ||
header.getRunNumber () != timeslice->getRunNumber ()) {
header = timeslice->getDAQHeader();
putObject(out, header);
}
if (background == rates_t) {
summary.update(timeslice->getDAQHeader());
if (timeslice->getFrameIndex() != summary.getSummaryslice()->getFrameIndex()) {
ERROR("Frame indices do not match at "
<< "[counter = " << counter << "]"
<< "[timeslice = " << timeslice->getFrameIndex() << "]"
<< "[summaryslice = " << summary.getSummaryslice()->getFrameIndex() << "]" << endl);
continue;
}
}
for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
if (router.hasModule(frame->getModuleID())) {
const JModule& module = router.getModule(frame->getModuleID());
//-----------------------------------------------------------
// determine background rates and veto per PMT
//-----------------------------------------------------------
vector<double> rate_Hz(NUMBER_OF_PMTS, 0.0);
vector<bool> veto (NUMBER_OF_PMTS, false);
if (background == count_t) {
for (JDAQSuperFrame::const_iterator i = frame->begin(); i != frame->end(); ++i) {
rate_Hz[i->getPMT()] += RTU * 1e9 / getFrameTime();
}
} else if (background == tails_t) {
// see below
} else if (background == rndms_t) {
// see below
} else if (background == rates_t) {
if (summary.hasSummaryFrame(frame->getModuleID())) {
const JDAQSummaryFrame& sum = summary.getSummaryFrame(frame->getModuleID());
for (int i = 0; i != NUMBER_OF_PMTS; ++i){
rate_Hz[i] = sum.getRate (i) *RTU;
veto [i] = sum.getValue(i) == 0;
}
} else {
FATAL("Summary frame of module " << frame->getModuleID() << " not found at frame index " << timeslice->getFrameIndex() << endl);
}
}
const JDAQFrameStatus status = frame->getDAQFrameStatus();
// Set veto according DAQ or PMT status and rate;
// Hits of PMTs with veto will not be counted in livetime nor coincidences nor background
for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
if (!getDAQStatus(status, module, i) ||
!getPMTStatus(status, module, i) ||
!rateRange_Hz(rate_Hz[i])) {
veto[i] = true;
}
}
const JHistogram& histogram = zmap[router.getAddress(frame->getModuleID()).first];
const JCombinatorics& combinatorics = histogram.getCombinatorics();
TH2D* h2s = histogram.h2s;
TH1D* h1b = histogram.h1b;
TH1D* h1L = histogram.h1L;
//-----------------------------------------------------------
// process data
//-----------------------------------------------------------
JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*frame, module);
buffer.preprocess(option, match);
for (JSuperFrame2D_t::iterator i = buffer.begin(); i != buffer.end(); ++i) {
if (veto[i->getPMTAddress()]) {
i->reset();
}
}
JSuperFrame1D_t& data = JSuperFrame1D_t::multiplex(buffer);
DEBUG("Number of hits " << timeslice->getFrameIndex() << ":" << frame->getModuleID() << ' ' << frame->size() << ' ' << data.size() << endl);
// Signal;
// Hits from PMTs that do not comply with rate cuts have been exluded above
size_t numberOfSignalEvents = 0;
double t1 = numeric_limits<double>::lowest();
for (vector<hit_type>::const_iterator p = data.begin(); p != data.end(); ) {
vector<hit_type>::const_iterator q = p;
while (++q != data.end() && q->getT() - p->getT() <= TMax_ns) {}
const int N = distance(p,q);
if (multiplicity(N)) {
numberOfSignalEvents += 1;
if (p->getT() > t1 + deadTime_ns) {
for (vector<hit_type>::const_iterator __p = p; __p != q; ++__p) {
for (vector<hit_type>::const_iterator __q = __p; ++__q != q; ) {
if ((__p->getToT() >= totRange_ns[0]) &&
(__q->getToT() >= totRange_ns[0]) &&
(__p->getToT() >= totRange_ns[1] ||
__q->getToT() >= totRange_ns[1]) &&
(__p->getToT() <= totRange_ns[2]) &&
(__q->getToT() <= totRange_ns[2])) {
h2s->Fill((double) combinatorics.getIndex(__p->getPMT(),__q->getPMT()),
JCombinatorics::getSign(__p->getPMT(),__q->getPMT()) * (__q->getT() - __p->getT()));
}
}
}
t1 = p->getT();
}
}
p = q;
}
// Background;
// Note that rate cuts and veto have already been accounted for when data buffer was filled
if (background == rndms_t) {
for (vector<hit_type>::iterator i = data.begin(); i != data.end(); ++i) {
*i = hit_type(i->getPMT(), JHit(i->getT() + t0[i->getPMT()], i->getToT()));
}
sort(data.begin(), data.end());
double t1 = numeric_limits<double>::lowest();
for (vector<hit_type>::const_iterator p = data.begin(); p != data.end(); ) {
vector<hit_type>::const_iterator q = p;
while (++q != data.end() && q->getT() - p->getT() <= TMax_ns) {}
const int N = distance(p,q);
if (multiplicity(N)) {
if (p->getT() > t1 + deadTime_ns) {
for (vector<hit_type>::const_iterator __p = p; __p != q; ++__p) {
for (vector<hit_type>::const_iterator __q = __p; ++__q != q; ) {
if ((__p->getToT() >= totRange_ns[0]) &&
(__q->getToT() >= totRange_ns[0]) &&
(__p->getToT() >= totRange_ns[1] ||
__q->getToT() >= totRange_ns[1]) &&
(__p->getToT() <= totRange_ns[2]) &&
(__q->getToT() <= totRange_ns[2])) {
h1b->Fill((double) combinatorics.getIndex(p->getPMT(),q->getPMT()), 1.0);
}
}
}
t1 = p->getT();
}
}
p = q;
}
} else if (background == rates_t ||
background == count_t) {
double Rs = 0.0; // total rate [Hz]
for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
if (!veto[i]) {
Rs += rate_Hz[i]; // [Hz]
}
}
for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
for (int j = i; ++j != NUMBER_OF_PMTS; ) {
if (!veto[i] && !veto[j]) {
const double R1 = rate_Hz[i]; // [Hz]
const double R2 = rate_Hz[j]; // [Hz]
// evaluate expected counts within a time window of 2 * TMax_ns for the two PMTs and the other PMTs inside the optical module;
// expected rate due to random coincidences then is product of the probability for the specific observation and the number of trials
// the observation is one hit in PMT 1, one hit in PMT 2 and a number of hits in the other PMTs inside the same module
// corresponding to the given multiplicity range
const double N = getP((R1) * 2 * TMax_ns * 1e-9,
(R2) * 2 * TMax_ns * 1e-9,
(Rs - R1 - R2) * 2 * TMax_ns * 1e-9,
multiplicity.getLowerLimit(),
multiplicity.getUpperLimit()) * getFrameTime() / (2*TMax_ns);
h1b->Fill((double) combinatorics.getIndex(i,j), N);
}
}
}
}
//-----------------------------------------------------------
// fill the livetime by PMT pairs
//-----------------------------------------------------------
const double livetime_s = getFrameTime()*1e-9 * exp(-deadTime_ns * numberOfSignalEvents/getFrameTime());
for (size_t i = 0; i < combinatorics.getNumberOfPairs(); ++i) {
const int pmt1 = combinatorics.getPair(i).first;
const int pmt2 = combinatorics.getPair(i).second;
if (!veto[pmt1] && !veto[pmt2]) {
h1L->Fill(i, livetime_s);
}
}
}
}
}
}
STATUS(endl);
if (background == tails_t ) {
const double R = (2.0*TMax_ns) / ((ymax - ymin)/ny);
for (vector<JHistogram>::iterator hr = zmap.begin() ; hr != zmap.end() ; ++hr) {
if (hr->is_valid()) {
TH2D* h2s = hr->h2s;
TH1D* h1b = hr->h1b;
for (int ix = 1; ix <= h1b->GetXaxis()->GetNbins(); ++ix) {
double Y = 0.0; // sum of counts in tail regions
double W = 0.0; // square of error of Y
int N = 0; // number of bins in tail regions
for (int iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {
const double y = h2s->GetYaxis()->GetBinCenter(iy) ;
if (Tail_ns(fabs(y))) {
Y += h2s->GetBinContent(ix,iy);
W += h2s->GetBinError (ix,iy) * h2s->GetBinError(ix,iy);
N += 1;
}
}
h1b->SetBinContent(ix, Y/N * R);
h1b->SetBinError (ix, sqrt(W/N) * R);
}
}
}
}
//---------------------------------------------
// save normalisation constants and store histograms
//---------------------------------------------
TH1D weights_hist(weights_hist_t, NULL, 3, 0.0, 3.0);
weights_hist.GetXaxis()->SetBinLabel(1, W1_overall_t); // [s]
weights_hist.GetXaxis()->SetBinLabel(2, WS_t); // [ns]
weights_hist.GetXaxis()->SetBinLabel(3, WB_t); // [ns]
weights_hist.Fill(W1_overall_t, counter*getFrameTime()*1e-9); // [s]
weights_hist.Fill(WS_t, (ymax - ymin)/ny); // [ns]
weights_hist.Fill(WB_t, 2.0*TMax_ns); // [ns]
out << weights_hist;
for (vector<JHistogram>::iterator i = zmap.begin(); i != zmap.end(); ++i) {
if (i->is_valid()) {
out << *(i->h2s);
out << *(i->h1b);
out << *(i->h1L);
}
}
out.Write();
out.Close();
}