59 static const char*
const rates_t =
"rates";
60 static const char*
const count_t =
"counts";
61 static const char*
const tails_t =
"tails";
62 static const char*
const rndms_t =
"randoms";
67 static const double ToTmax_ns = 250;
103 const Int_t nx = this->getNumberOfPairs();
104 const Double_t
xmin = -0.5;
105 const Double_t
xmax = nx - 0.5;
123 return (h2s != NULL &&
144 inline double getP(
const double E1,
const double E2,
const double ED,
const int M_min,
const int M_max)
148 double P = E1 * E2 * exp(-(E1 + E2 + ED));
154 for ( ; MD + 2 <= M_min; ++MD) {
162 for ( ; MD + 2 <= M_max; ++MD) {
201 int main(
int argc,
char **argv)
224 JParser<> zap(
"Auxiliary program to determine PMT parameters from K40 data.");
226 zap[
'f'] =
make_field(inputFile,
"input file.");
228 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
229 zap[
'a'] =
make_field(detectorFile,
"detector file.");
230 zap[
'T'] =
make_field(TMax_ns,
"time window [ns].") = 20.0;
232 zap[
't'] =
make_field(totRange_ns,
"PMT time-over-threshold range [ns].") = { 4.0, 7.0, ToTmax_ns };
233 zap[
'b'] =
make_field(background,
"background estimation method.") = rates_t, count_t, tails_t, rndms_t;
236 zap[
'D'] =
make_field(deadTime_us,
"L1 dead time (us)") = 0.0;
237 zap[
'C'] =
make_field(selector,
"timeslice selector, e.g. JDAQTimesliceL1.") = getROOTClassSelection<JDAQTimesliceTypes_t>();
238 zap[
'O'] =
make_field(option,
"hit pre-processing option.") = JPreprocessor::getOptions(JPreprocessor::filter_t);
243 catch(
const exception &error) {
244 FATAL(error.what() << endl);
252 if (selector == JDAQTimesliceL2::Class_Name() ||
253 selector == JDAQTimesliceSN::Class_Name()) {
254 FATAL(
"Option -C <selector> " << selector <<
" not compatible with calibration method." << endl);
257 if (selector == JDAQTimeslice ::Class_Name() ||
258 selector == JDAQTimesliceL1::Class_Name()) {
266 FATAL(
"No trigger parameters from input:" << error.
what() << endl);
269 if ((selector == JDAQTimeslice ::Class_Name() && parameters.writeL1.prescale > 0) ||
270 (selector == JDAQTimesliceL1::Class_Name())) {
272 if (parameters.TMaxLocal_ns < TMax_ns) {
273 FATAL(
"Option -T <TMax_ns> = " << TMax_ns <<
" is larger than in the trigger " << parameters.TMaxLocal_ns << endl);
276 if (background == rndms_t ||
277 background == count_t) {
278 FATAL(
"Option -C <selector> " << selector <<
" incompatible with option -b <background> " << background << endl);
284 FATAL(
"Invalid option -M <multiplicity> " << multiplicity << endl);
287 if (totRange_ns[0] > totRange_ns[1] ||
288 totRange_ns[1] > totRange_ns[2]) {
289 FATAL(
"Invalid option -t <totRange_ns> " <<
JEEPZ() << totRange_ns << endl);
292 if (background == tails_t) {
295 FATAL(
"Invalid option -B <Tail_ns> " << Tail_ns << endl);
299 FATAL(
"Invalid option -B <Tail_ns> " << Tail_ns <<
"; upper limit larger than option -T <TMax_ns> " << TMax_ns << endl);
317 FATAL(
"Empty detector." << endl);
329 t0.push_back(i * 2 * TMax_ns);
350 const Double_t ymin = -floor(TMax_ns) + 0.5;
351 const Double_t ymax = +floor(TMax_ns) - 0.5;
352 const Int_t ny = (Int_t) ((ymax - ymin) / 1.0);
354 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ++module) {
358 if (!module->empty()) {
360 NOTICE(
"Booking histograms for module " << module->
getID() << endl);
373 const double deadTime_ns = deadTime_us * 1e3;
382 for (JMultipleFileScanner_t::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
384 STATUS(
"Processing: " << *i << endl);
393 STATUS(
"event: " << setw(10) << counter <<
'\r');
DEBUG(endl);
405 if (background == rates_t) {
411 ERROR(
"Frame indices do not match at "
412 <<
"[counter = " << counter <<
"]"
420 for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
422 if (router.
hasModule(frame->getModuleID())) {
433 if (background == count_t) {
439 }
else if (background == tails_t) {
443 }
else if (background == rndms_t) {
447 }
else if (background == rates_t) {
454 rate_Hz[i] = sum.
getRate (i) *RTU;
460 FATAL(
"Summary frame of module " << frame->getModuleID() <<
" not found at frame index " << timeslice->
getFrameIndex() << endl);
472 !rateRange_Hz(rate_Hz[i])) {
478 const JCombinatorics& combinatorics = histogram.getCombinatorics();
480 TH2D* h2s = histogram.h2s;
481 TH1D* h1b = histogram.h1b;
482 TH1D* h1L = histogram.h1L;
488 JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*frame, module);
490 buffer.preprocess(option, match);
492 for (JSuperFrame2D_t::iterator i = buffer.begin(); i != buffer.end(); ++i) {
493 if (veto[i->getPMTAddress()]) {
498 JSuperFrame1D_t&
data = JSuperFrame1D_t::multiplex(buffer);
500 DEBUG(
"Number of hits " << timeslice->
getFrameIndex() <<
":" << frame->getModuleID() <<
' ' << frame->size() <<
' ' <<
data.size() << endl);
505 size_t numberOfSignalEvents = 0;
507 double t1 = numeric_limits<double>::lowest();
513 while (++q !=
data.end() && q->getT() - p->getT() <= TMax_ns) {}
517 if (multiplicity(N)) {
519 numberOfSignalEvents += 1;
521 if (p->getT() > t1 + deadTime_ns) {
526 if ((__p->getToT() >= totRange_ns[0]) &&
527 (__q->getToT() >= totRange_ns[0]) &&
528 (__p->getToT() >= totRange_ns[1] ||
529 __q->getToT() >= totRange_ns[1]) &&
530 (__p->getToT() <= totRange_ns[2]) &&
531 (__q->getToT() <= totRange_ns[2])) {
533 h2s->Fill((
double) combinatorics.
getIndex(__p->getPMT(),__q->getPMT()),
534 JCombinatorics::getSign(__p->getPMT(),__q->getPMT()) * (__q->getT() - __p->getT()));
549 if (background == rndms_t) {
552 *i =
hit_type(i->getPMT(),
JHit(i->getT() + t0[i->getPMT()], i->getToT()));
557 double t1 = numeric_limits<double>::lowest();
563 while (++q !=
data.end() && q->getT() - p->getT() <= TMax_ns) {}
567 if (multiplicity(N)) {
569 if (p->getT() > t1 + deadTime_ns) {
574 if ((__p->getToT() >= totRange_ns[0]) &&
575 (__q->getToT() >= totRange_ns[0]) &&
576 (__p->getToT() >= totRange_ns[1] ||
577 __q->getToT() >= totRange_ns[1]) &&
578 (__p->getToT() <= totRange_ns[2]) &&
579 (__q->getToT() <= totRange_ns[2])) {
581 h1b->Fill((
double) combinatorics.
getIndex(p->getPMT(),q->getPMT()), 1.0);
593 }
else if (background == rates_t ||
594 background == count_t) {
607 if (!veto[i] && !veto[
j]) {
609 const double R1 = rate_Hz[i];
610 const double R2 = rate_Hz[
j];
617 const double N =
getP((
R1) * 2 * TMax_ns * 1e-9,
618 (R2) * 2 * TMax_ns * 1e-9,
619 (Rs -
R1 - R2) * 2 * TMax_ns * 1e-9,
623 h1b->Fill((
double) combinatorics.
getIndex(i,
j), N);
637 const int pmt1 = combinatorics.
getPair(i).first;
638 const int pmt2 = combinatorics.
getPair(i).second;
640 if (!veto[pmt1] && !veto[pmt2]) {
641 h1L->Fill(i, livetime_s);
650 if (background == tails_t ) {
652 const double R = (2.0*TMax_ns) / ((ymax - ymin)/ny);
656 if (hr->is_valid()) {
661 for (
int ix = 1; ix <= h1b->GetXaxis()->GetNbins(); ++ix) {
667 for (
int iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {
669 const double y = h2s->GetYaxis()->GetBinCenter(iy) ;
671 if (Tail_ns(fabs(
y))) {
672 Y += h2s->GetBinContent(ix,iy);
673 W += h2s->GetBinError (ix,iy) * h2s->GetBinError(ix,iy);
678 h1b->SetBinContent(ix, Y/N * R);
679 h1b->SetBinError (ix, sqrt(W/N) * R);
692 weights_hist.GetXaxis()->SetBinLabel(2,
WS_t);
693 weights_hist.GetXaxis()->SetBinLabel(3,
WB_t);
696 weights_hist.Fill(
WS_t, (ymax - ymin)/ny);
697 weights_hist.Fill(
WB_t, 2.0*TMax_ns);
int main(int argc, char **argv)
KM3NeT DAQ constants, bit handling, etc.
Data structure for detector geometry and calibration.
Basic data structure for L0 hit.
Match operator for consecutive hits.
General purpose messaging.
#define DEBUG(A)
Message macros.
Direct access to module in detector data structure.
Scanning of objects from multiple files according a format that follows from the extension of each fi...
PMT analogue signal processor.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Auxiliaries for pre-processing of hits.
I/O formatting auxiliaries.
#define MAKE_CSTRING(A)
Make C-string.
Auxiliary class to define a range between two values.
Scanning of objects from a single file according a format that follows from the extension of each fil...
ROOT TTree parameter settings of various packages.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Address of module in detector data structure.
int first
index of module in detector data structure
Router for direct addressing of module data in detector data structure.
const JModuleAddress & getAddress(const JObjectID &id) const
Get address of module.
const JModule & getModule(const JObjectID &id) const
Get module parameters.
bool hasModule(const JObjectID &id) const
Has module.
Data structure for a composite optical module.
virtual const char * what() const override
Get error message.
Auxiliary class for object identification.
int getID() const
Get identifier.
Auxiliary class for multiplexing object iterators.
virtual bool hasNext() override
Check availability of next element.
virtual const pointer_type & next() override
Get next element.
Utility class to parse command line options.
Object reading from a list of files.
File router for fast addressing of summary data.
void update(const JDAQHeader &header)
Update router.
bool hasSummaryFrame(const JDAQModuleIdentifier &module) const
Has summary frame.
const JDAQSummaryFrame & getSummaryFrame() const
Get default summary frame.
const JDAQSummaryslice * getSummaryslice() const
Get summary slice.
Reduced data structure for L0 hit.
1-dimensional frame with time calibrated data from one optical module.
2-dimensional frame with time calibrated data from one optical module.
int getDetectorID() const
Get detector identifier.
int getRunNumber() const
Get run number.
int getFrameIndex() const
Get frame index.
const JDAQFrameStatus & getDAQFrameStatus() const
Get DAQ frame status.
Data storage class for rate measurements of all PMTs in one module.
double getRate(const int tdc, const double factor=1.0) const
Get count rate.
JRate_t getValue(const int tdc) const
Get value.
static const char *const _2S
Name extension for 2D counts.
static const char *const WB_t
Named bin for value of TMax_ns in JCalibrateK40.
static const char *const weights_hist_t
Histogram naming.
static const char *const _1B
Name extension for 1D background.
static const char *const WS_t
Named bin for time residual bin width.
static const char *const W1_overall_t
Named bin for duration of the run.
static const char *const _1L
Name extension for 1D live time.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
double getP(const double expval, bool hit)
Get Poisson probability to observe a hit or not for given expectation value for the number of hits.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
bool putObject(TDirectory &dir, const TObject &object)
Write object to ROOT directory.
Long64_t counter_type
Type definition for counter.
bool is_valid(const json &js)
Check validity of JSon data.
JTriggerParameters getTriggerParameters(const JMultipleFileScanner_t &file_list)
Get trigger parameters.
bool getDAQStatus(const JDAQFrameStatus &frame, const JStatus &status)
Test status of DAQ.
bool getPMTStatus(const JStatus &status)
Test status of PMT.
KM3NeT DAQ data structures and auxiliaries.
double getFrameTime()
Get frame time duration.
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Auxiliary class to sort pairs of PMT addresses within optical module.
PMT analogue signal processor.
virtual double getNPE(const double tot_ns) const override
Get number of photo-electrons.
double getIntegralOfChargeProbability(const double xmin, const double xmax, const int NPE) const
Get integral of probability.
Auxiliary data structure for streaming of STL containers.
Auxiliary class to select ROOT class based on class name.
Auxiliary base class for list of file names.
Auxiliary class for specifying the way of pre-processing of hits.