106 using namespace KM3NETDAQ;
113 JMultipleFileScanner<JDAQTimesliceTypes_t> inputFile;
118 JRange<double> rateRange_Hz;
119 JRange<double> totRange_ns;
120 JRange<double> Tail_ns;
121 JRange<int> multiplicity;
123 JROOTClassSelector selector;
129 JParser<> zap(
"Auxiliary program to determine PMT parameters from K40 data.");
131 zap[
'f'] =
make_field(inputFile,
"input file.");
133 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
134 zap[
'a'] =
make_field(detectorFile,
"detector file.");
135 zap[
'T'] =
make_field(TMax_ns,
"time window [ns].") = 20.0;
136 zap[
'V'] =
make_field(rateRange_Hz,
"PMT rate range [Hz].") = JRange<double>(0.0, 20.0e3);
137 zap[
't'] =
make_field(totRange_ns,
"PMT time-over-threshold range [ns].") = JRange<double>(0.0, ToTmax_ns);
138 zap[
'b'] =
make_field(background,
"background estimation method.") = rates_t, count_t, tails_t, rndms_t;
139 zap[
'B'] =
make_field(Tail_ns,
"time window used for background estimation.") = JRange<double>(15.0, 20.0);
140 zap[
'M'] =
make_field(multiplicity,
"multiplicity range of hits on DOM.") = JRange<int>(2, 31);
141 zap[
'D'] =
make_field(deadTime_us,
"L1 dead time (us)") = 0.0;
142 zap[
'C'] =
make_field(selector,
"timeslice selector, e.g. JDAQTimesliceL1.") = getROOTClassSelection<JDAQTimesliceTypes_t>();
147 catch(
const exception &error) {
148 FATAL(error.what() << endl);
160 catch(
const JException& error) {
161 FATAL(
"No trigger parameters from input:" << error.what() << endl);
166 if (selector == JDAQTimesliceL2::Class_Name() ||
167 selector == JDAQTimesliceSN::Class_Name()) {
168 FATAL(
"Option -C <selector> " << selector <<
" not compatible with calibration method." << endl);
171 if ((selector == JDAQTimeslice ::Class_Name() && parameters.writeL1.prescale > 0) ||
172 (selector == JDAQTimesliceL1::Class_Name())) {
174 if (parameters.TMaxLocal_ns < TMax_ns) {
175 FATAL(
"Option -T <TMax_ns> = " << TMax_ns <<
" is larger than in the trigger " << parameters.TMaxLocal_ns << endl);
178 if (background == rndms_t ||
179 background == count_t) {
180 FATAL(
"Option -C <selector> " << selector <<
" incompatible with option -b <background> " << background << endl);
184 if (!multiplicity.is_valid() || multiplicity.getLowerLimit() < 2) {
185 FATAL(
"Invalid option -M <multiplicity> " << multiplicity << endl);
188 if (!totRange_ns.is_valid()) {
189 FATAL(
"Invalid option -t <totRange_ns> " << totRange_ns << endl);
192 if (background == tails_t) {
194 if (!Tail_ns.is_valid()) {
195 FATAL(
"Invalid option -B <Tail_ns> " << Tail_ns << endl);
198 if (Tail_ns.getUpperLimit() > TMax_ns) {
199 FATAL(
"Invalid option -B <Tail_ns> " << Tail_ns <<
"; upper limit larger than option -T <TMax_ns> " << TMax_ns << endl);
210 load(detectorFile, detector);
212 catch(
const JException& error) {
216 if (detector.empty()) {
217 FATAL(
"Empty detector." << endl);
220 const JModuleRouter router(detector);
221 JSummaryRouter summaryRouter;
230 t0.push_back(i * 2 * TMax_ns);
233 const JMatch_t match(TMax_ns);
239 const JPMTAnalogueSignalProcessor cpu;
242 const double RTU = (cpu.getIntegralOfProbability(cpu.getNPE(totRange_ns.getLowerLimit()), cpu.getNPE(totRange_ns.getUpperLimit()), NPE)
244 cpu.getIntegralOfProbability(cpu.getNPE(0.0), cpu.getNPE(ToTmax_ns), NPE));
259 weights_hist.GetXaxis()->SetBinLabel(2,
WS_t);
260 weights_hist.GetXaxis()->SetBinLabel(3,
WB_t);
262 const int nx = combinatorics.getNumberOfPairs();
263 const double xmin = -0.5;
264 const double xmax = nx - 0.5;
266 const double ymin = -floor(TMax_ns) + 0.5;
267 const double ymax = +floor(TMax_ns) - 0.5;
268 const int ny = (int) ((ymax - ymin) / 1.0);
271 for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
273 const JModuleAddress& address = router.getAddress(module->getID());
275 NOTICE(
"Booking histograms for module " << module->getID() << endl);
277 zmap[address.first] = JHistogram(
new TH2D(
MAKE_CSTRING(module->getID() <<
_2S), NULL, nx, xmin, xmax, ny, ymin, ymax),
283 JTreeScanner<JDAQSummaryslice, JDAQEvaluator> summaryFile(inputFile);
288 const double deadTime_ns = deadTime_us * 1e3;
290 JObjectMultiplexer<JDAQTimesliceTypes_t, JDAQTimeslice> in(inputFile, selector);
294 for ( ; in.hasNext() && counter != inputFile.getLimit(); ++counter) {
296 STATUS(
"event: " << setw(10) << counter <<
'\r');
DEBUG(endl);
300 if (background == rates_t) {
302 const Long64_t index = summaryFile.find(*timeslice);
308 summaryRouter.update(summaryFile.getEntry(index));
310 if (timeslice->
getFrameIndex() != summaryRouter.getFrameIndex()) {
312 ERROR(
"Frame indices do not match "
313 <<
"[counter=" << counter <<
"]"
315 <<
"[summaryslice=" << summaryRouter.getFrameIndex() <<
"]" << endl);
320 FATAL(
"ROOT TTrees not aligned; Run number " << timeslice->
getRunNumber() << endl);
325 for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
327 if (router.hasModule(frame->getModuleID())) {
338 if (background == count_t) {
346 }
else if (background == rates_t) {
348 if (!summaryRouter.hasSummaryFrame(frame->getModuleID())) {
349 FATAL(
"Summary frame not found " << timeslice->
getFrameIndex() <<
":" << frame->getModuleID() << endl);
352 const JDAQSummaryFrame& summary_frame = summaryRouter.getSummaryFrame(frame->getModuleID());
357 rate_Hz[i] = RTU * summary_frame.
getRate(i);
384 const JModuleAddress& address = router.getAddress(frame->getModuleID());
385 const JModule& module = detector.getModule(address);
387 combinatorics.configure(module.size());
389 combinatorics.sort(JPairwiseComparator(module));
391 TH2D* h2s = zmap[address.first].h2s;
392 TH1D* h1b = zmap[address.first].h1b;
393 TH1D* h1L = zmap[address.first].h1L;
399 for (
size_t i = 0; i < combinatorics.getNumberOfPairs(); ++i) {
401 const int pmt1 = combinatorics.getPair(i).first;
402 const int pmt2 = combinatorics.getPair(i).second;
404 if (!veto[pmt1] && rateRange_Hz(rate_Hz[pmt1]) &&
405 !veto[pmt2] && rateRange_Hz(rate_Hz[pmt2]) ) {
415 JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*frame, module);
417 for (JSuperFrame2D_t::iterator i = buffer.begin(); i != buffer.end(); ++i) {
419 const int pmt = i->getPMTAddress();
421 if ( veto[pmt] || !rateRange_Hz(rate_Hz[pmt]) ) {
433 DEBUG(
"Number of hits " << timeslice->
getFrameIndex() <<
":" << frame->getModuleID() <<
' ' << frame->size() <<
' ' << data.size() << endl);
442 double t1 = -numeric_limits<double>::max();
448 while (++q != data.end() && q->getT() - p->getT() <= TMax_ns) {}
450 const int N = distance(p,q);
452 if (multiplicity(N)) {
454 if (p->getT() > t1 + deadTime_ns) {
459 if (totRange_ns( __p->getToT()) && totRange_ns(__q->getToT())) {
460 h2s->Fill((
double) combinatorics.getIndex(__p->getPMT(),__q->getPMT()),
461 JCombinatorics::getSign(__p->getPMT(),__q->getPMT()) * (__q->getT() - __p->getT()));
476 if (background == rndms_t) {
479 *i = JHitR0(i->getPMT(), JHit(i->getT() + t0[i->getPMT()], i->getToT()));
482 sort(data.begin(), data.end());
488 while (++q != data.end() && q->getT() - p->getT() <= TMax_ns) {}
490 const int N = distance(p,q);
492 if (multiplicity(N)) {
497 if (totRange_ns(__p->getToT()) && totRange_ns(__q->getToT())) {
498 h1b->Fill((
double) combinatorics.getIndex(p->getPMT(),q->getPMT()), 1.0);
507 }
else if (background == rates_t ||
508 background == count_t) {
514 const Double_t R1 = rate_Hz[i];
516 if (!veto[i] && rateRange_Hz(R1)) {
524 const Double_t R1 = rate_Hz[i];
525 const Double_t R2 = rate_Hz[j];
527 if (!veto[i] && rateRange_Hz(R1) &&
528 !veto[j] && rateRange_Hz(R2)) {
531 const double ED = (Rs - rate_Hz[i] - rate_Hz[j]) * 2 * TMax_ns * 1e-9;
532 const double E1 = rate_Hz[i] * 2 * TMax_ns * 1e-9;
533 const double E2 = rate_Hz[j] * 2 * TMax_ns * 1e-9;
537 multiplicity.getLowerLimit(),
538 multiplicity.getUpperLimit()) *
getFrameTime() / (2*TMax_ns);
540 h1b->Fill((
double) combinatorics.getIndex(i,j), N);
550 if (background == tails_t ) {
554 TH2D* h2s = hist->h2s;
555 TH1D* h1b = hist->h1b;
563 int k = combinatorics.getIndex(i,j);
565 for (
int l = 1; l <= h2s->GetYaxis()->GetNbins(); ++l) {
567 const Double_t dt = h2s->GetYaxis()->GetBinCenter(l) ;
569 if (Tail_ns(fabs(dt))) {
570 Y += h2s->GetBinContent(k+1,l);
571 W += h2s->GetBinError(k+1,l) * h2s->GetBinError(k+1,l);
576 h1b->SetBinContent(k+1, Y / N * (2.0*TMax_ns) / ((ymax - ymin)/ny) );
577 h1b->SetBinError (k+1, sqrt(W/N) * (2.0*TMax_ns) / ((ymax - ymin)/ny) );
588 weights_hist.Fill(
WS_t, (ymax - ymin)/ny);
589 weights_hist.Fill(
WB_t, 2.0*TMax_ns);
596 weights_hist.Write() ;
Utility class to parse command line options.
Data structure for all trigger parameters.
double getRate(const int tdc, const double factor=1.0) const
Get count rate.
static const char *const WB_t
Named bin for value of TMax_ns in JCalibrateK40.
#define MAKE_CSTRING(A)
Make C-string.
Long64_t counter_type
Type definition for counter.
JSuperFrame2D< hit_type > JSuperFrame2D_t
int getRunNumber() const
Get run number.
int getFrameIndex() const
Get frame index.
JLimit JLimit_t
Type definition of limit.
static const char *const WS_t
Named bin for time residual bin width.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
static const char *const _1B
Name extension for 1D background.
double getFrameTime()
Get frame time duration.
Data storage class for rate measurements of all PMTs in one module.
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
double coincidenceP(Double_t E1, Double_t E2, Double_t ED, int M_min, int M_max)
Coincidence probability of two PMTs.
bool testHighRateVeto() const
Test high-rate veto status.
static const char *const weights_hist_t
Histogram naming.
static const char *const _1L
Name extension for 1D live time.
static const char *const W1_overall_t
Named bin for duration of the run.
const JLimit & getLimit() const
Get limit.
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
const JDAQFrameStatus & getDAQFrameStatus() const
Get DAQ frame status.
JTriggerParameters getTriggerParameters(const JMultipleFileScanner_t &file_list)
Get trigger parameters.
static const char *const _2S
Name extension for 2D counts.
#define DEBUG(A)
Message macros.
JSuperFrame1D< hit_type > JSuperFrame1D_t
bool testFIFOStatus() const
Test FIFO status.
bool putObject(TDirectory *dir, const T &object)
Write object to ROOT directory.