57 static const char* rates_t =
"rates";
58 static const char* count_t =
"counts";
59 static const char* tails_t =
"tails";
60 static const char* rndms_t =
"randoms";
65 static const double ToTmax_ns = std::numeric_limits<KM3NETDAQ::JDAQHit::JTOT_t>::max();
99 const Int_t nx = this->getNumberOfPairs();
100 const Double_t
xmin = -0.5;
101 const Double_t
xmax = nx - 0.5;
119 return (h2s != NULL &&
140 inline double getP(
const double E1,
const double E2,
const double ED,
const int M_min,
const int M_max)
144 double P = E1 * E2 *
exp(-(E1 + E2 + ED));
150 for ( ; MD + 2 <= M_min; ++MD) {
158 for ( ; MD + 2 <= M_max; ++MD) {
192 int main(
int argc,
char **argv)
196 using namespace KM3NETDAQ;
215 JParser<> zap(
"Auxiliary program to determine PMT parameters from K40 data.");
217 zap[
'f'] =
make_field(inputFile,
"input file.");
219 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
220 zap[
'a'] =
make_field(detectorFile,
"detector file.");
221 zap[
'T'] =
make_field(TMax_ns,
"time window [ns].") = 20.0;
224 zap[
'b'] =
make_field(background,
"background estimation method.") = rates_t, count_t, tails_t, rndms_t;
227 zap[
'D'] =
make_field(deadTime_us,
"L1 dead time (us)") = 0.0;
228 zap[
'C'] =
make_field(selector,
"timeslice selector, e.g. JDAQTimesliceL1.") = getROOTClassSelection<JDAQTimesliceTypes_t>();
229 zap[
'O'] =
make_field(option,
"hit pre-processing option.") = JPreprocessor::getOptions();
234 catch(
const exception &error) {
235 FATAL(error.what() << endl);
243 if (selector == JDAQTimesliceL2::Class_Name() ||
244 selector == JDAQTimesliceSN::Class_Name()) {
245 FATAL(
"Option -C <selector> " << selector <<
" not compatible with calibration method." << endl);
248 if (selector == JDAQTimeslice ::Class_Name() ||
249 selector == JDAQTimesliceL1::Class_Name()) {
257 FATAL(
"No trigger parameters from input:" << error.
what() << endl);
260 if ((selector == JDAQTimeslice ::Class_Name() && parameters.writeL1.prescale > 0) ||
261 (selector == JDAQTimesliceL1::Class_Name())) {
263 if (parameters.TMaxLocal_ns < TMax_ns) {
264 FATAL(
"Option -T <TMax_ns> = " << TMax_ns <<
" is larger than in the trigger " << parameters.TMaxLocal_ns << endl);
267 if (background == rndms_t ||
268 background == count_t) {
269 FATAL(
"Option -C <selector> " << selector <<
" incompatible with option -b <background> " << background << endl);
274 if (!multiplicity.is_valid() || multiplicity.getLowerLimit() < 2 || multiplicity.getUpperLimit() >
NUMBER_OF_PMTS) {
275 FATAL(
"Invalid option -M <multiplicity> " << multiplicity << endl);
278 if (!totRange_ns.is_valid()) {
279 FATAL(
"Invalid option -t <totRange_ns> " << totRange_ns << endl);
282 if (background == tails_t) {
284 if (!Tail_ns.is_valid()) {
285 FATAL(
"Invalid option -B <Tail_ns> " << Tail_ns << endl);
288 if (Tail_ns.getUpperLimit() > TMax_ns) {
289 FATAL(
"Invalid option -B <Tail_ns> " << Tail_ns <<
"; upper limit larger than option -T <TMax_ns> " << TMax_ns << endl);
307 FATAL(
"Empty detector." << endl);
319 t0.push_back(i * 2 * TMax_ns);
329 const double RTU = (cpu.getIntegralOfChargeProbability(cpu.getNPE(totRange_ns.getLowerLimit()), cpu.getNPE(totRange_ns.getUpperLimit()), NPE)
331 cpu.getIntegralOfChargeProbability(cpu.getNPE(0.0), cpu.getNPE(ToTmax_ns), NPE));
340 const Double_t ymin = -floor(TMax_ns) + 0.5;
341 const Double_t ymax = +floor(TMax_ns) - 0.5;
342 const Int_t ny = (Int_t) ((ymax - ymin) / 1.0);
344 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ++module) {
348 if (!module->empty()) {
350 NOTICE(
"Booking histograms for module " << module->
getID() << endl);
363 const double deadTime_ns = deadTime_us * 1e3;
372 for (JMultipleFileScanner_t::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
374 STATUS(
"Processing: " << *i << endl);
383 STATUS(
"event: " << setw(10) << counter <<
'\r');
DEBUG(endl);
395 if (background == rates_t) {
401 ERROR(
"Frame indices do not match at "
402 <<
"[counter = " << counter <<
"]"
410 for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
412 if (router.
hasModule(frame->getModuleID())) {
423 if (background == count_t) {
429 }
else if (background == tails_t) {
433 }
else if (background == rndms_t) {
437 }
else if (background == rates_t) {
444 rate_Hz[i] = RTU * sum.
getRate(i);
449 FATAL(
"Summary frame of module " << frame->getModuleID() <<
" not found at frame index " << timeslice->
getFrameIndex() << endl);
461 !rateRange_Hz(rate_Hz[i])) {
467 const JCombinatorics& combinatorics = histogram.getCombinatorics();
469 TH2D* h2s = histogram.h2s;
470 TH1D* h1b = histogram.h1b;
471 TH1D* h1L = histogram.h1L;
477 JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*frame, module);
479 buffer.preprocess(option, match);
481 for (JSuperFrame2D_t::iterator i = buffer.begin(); i != buffer.end(); ++i) {
482 if (veto[i->getPMTAddress()]) {
487 JSuperFrame1D_t& data = JSuperFrame1D_t::multiplex(buffer);
489 DEBUG(
"Number of hits " << timeslice->
getFrameIndex() <<
":" << frame->getModuleID() <<
' ' << frame->size() <<
' ' << data.size() << endl);
494 size_t numberOfSignalEvents = 0;
496 double t1 = numeric_limits<double>::lowest();
502 while (++q != data.end() && q->getT() - p->getT() <= TMax_ns) {}
506 if (multiplicity(N)) {
508 numberOfSignalEvents += 1;
510 if (p->getT() > t1 + deadTime_ns) {
515 if (totRange_ns( __p->getToT()) && totRange_ns(__q->getToT())) {
516 h2s->Fill((
double) combinatorics.getIndex(__p->getPMT(),__q->getPMT()),
517 JCombinatorics::getSign(__p->getPMT(),__q->getPMT()) * (__q->getT() - __p->getT()));
532 if (background == rndms_t) {
535 *i =
hit_type(i->getPMT(),
JHit(i->getT() + t0[i->getPMT()], i->getToT()));
538 sort(data.begin(), data.end());
540 double t1 = numeric_limits<double>::lowest();
546 while (++q != data.end() && q->getT() - p->getT() <= TMax_ns) {}
550 if (multiplicity(N)) {
552 if (p->getT() > t1 + deadTime_ns) {
557 if (totRange_ns(__p->getToT()) && totRange_ns(__q->getToT())) {
558 h1b->Fill((
double) combinatorics.getIndex(p->getPMT(),q->getPMT()), 1.0);
570 }
else if (background == rates_t ||
571 background == count_t) {
584 if (!veto[i] && !veto[
j]) {
586 const double R1 = rate_Hz[i];
587 const double R2 = rate_Hz[
j];
594 const double N =
getP((R1) * 2 * TMax_ns * 1e-9,
595 (R2) * 2 * TMax_ns * 1e-9,
596 (Rs - R1 - R2) * 2 * TMax_ns * 1e-9,
597 multiplicity.getLowerLimit(),
598 multiplicity.getUpperLimit()) *
getFrameTime() / (2*TMax_ns);
600 h1b->Fill((
double) combinatorics.getIndex(i,j),
N);
612 for (
size_t i = 0; i < combinatorics.getNumberOfPairs(); ++i) {
614 const int pmt1 = combinatorics.getPair(i).first;
615 const int pmt2 = combinatorics.getPair(i).second;
617 if (!veto[pmt1] && !veto[pmt2]) {
618 h1L->Fill(i, livetime_s);
627 if (background == tails_t ) {
631 if (hist->is_valid()) {
635 TH2D* h2s = hist->h2s;
636 TH1D* h1b = hist->h1b;
646 for (
int l = 1; l <= h2s->GetYaxis()->GetNbins(); ++l) {
648 const Double_t dt = h2s->GetYaxis()->GetBinCenter(l) ;
650 if (Tail_ns(fabs(dt))) {
651 Y += h2s->GetBinContent(k+1,l);
652 W += h2s->GetBinError(k+1,l) * h2s->GetBinError(k+1,l);
657 h1b->SetBinContent(k+1, Y / N * (2.0*TMax_ns) / ((ymax - ymin)/ny) );
658 h1b->SetBinError (k+1, sqrt(W/N) * (2.0*TMax_ns) / ((ymax - ymin)/ny) );
672 weights_hist.GetXaxis()->SetBinLabel(2,
WS_t);
673 weights_hist.GetXaxis()->SetBinLabel(3,
WB_t);
676 weights_hist.Fill(
WS_t, (ymax - ymin)/ny);
677 weights_hist.Fill(
WB_t, 2.0*TMax_ns);
Utility class to parse command line options.
Auxiliaries for pre-processing of hits.
then fatal No hydrophone data file $HYDROPHONE_TXT fi sort gr k
int main(int argc, char *argv[])
ROOT TTree parameter settings of various packages.
double getRate(const int tdc, const double factor=1.0) const
Get count rate.
Basic data structure for L0 hit.
const JModule & getModule(const JObjectID &id) const
Get module parameters.
Data structure for a composite optical module.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
int getDetectorID() const
Get detector identifier.
static const char *const WB_t
Named bin for value of TMax_ns in JCalibrateK40.
bool is_valid(const json &js)
Check validity of JSon data.
void update(const JDAQHeader &header)
Update router.
Auxiliary class to select ROOT class based on class name.
Router for direct addressing of module data in detector data structure.
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
#define MAKE_CSTRING(A)
Make C-string.
Long64_t counter_type
Type definition for counter.
then fatal Wrong number of arguments fi set_variable STRING $argv[1] set_variable DETECTORXY_TXT $WORKDIR $DETECTORXY_TXT tail read X Y CHI2 RMS printf optimum n $X $Y $CHI2 $RMS awk v Y
const JDAQSummaryFrame & getSummaryFrame() const
Get default summary frame.
Auxiliary class for multiplexing object iterators.
int getRunNumber() const
Get run number.
Data structure for detector geometry and calibration.
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
int getFrameIndex() const
Get frame index.
1-dimensional frame with time calibrated data from one optical module.
int first
index of module in detector data structure
Scanning of objects from a single file according a format that follows from the extension of each fil...
static const char *const WS_t
Named bin for time residual bin width.
I/O formatting auxiliaries.
Auxiliary class to sort pairs of PMT addresses within optical module.
Reduced data structure for L0 hit.
const JDAQSummaryslice * getSummaryslice() const
Get summary slice.
#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.
int getID() const
Get identifier.
File router for fast addressing of summary data.
virtual const pointer_type & next() override
Get next element.
Address of module in detector data structure.
Match operator for consecutive hits.
virtual bool hasNext() override
Check availability of next element.
bool putObject(TDirectory &dir, const TObject &object)
Write object to ROOT directory.
General purpose messaging.
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...
Scanning of objects from multiple files according a format that follows from the extension of each fi...
Direct access to module in detector data structure.
PMT analogue signal processor.
static const char *const weights_hist_t
Histogram naming.
static const char *const _1L
Name extension for 1D live time.
Auxiliary base class for list of file names.
virtual const char * what() const override
Get error message.
bool getPMTStatus(const JStatus &status)
Test status of PMT.
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Auxiliary class to define a range between two values.
Utility class to parse command line options.
bool hasSummaryFrame(const JDAQModuleIdentifier &module) const
Has summary frame.
static const char *const W1_overall_t
Named bin for duration of the run.
const JModuleAddress & getAddress(const JObjectID &id) const
Get address of module.
bool hasModule(const JObjectID &id) const
Has module.
2-dimensional frame with time calibrated data from one optical module.
PMT analogue signal processor.
Object reading from a list of files.
do set_variable DETECTOR_TXT $WORKDIR detector
bool getDAQStatus(const JDAQFrameStatus &frame, const JStatus &status)
Test status of DAQ.
KM3NeT DAQ constants, bit handling, etc.
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
Auxiliary class for specifying the way of pre-processing of hits.
const JDAQFrameStatus & getDAQFrameStatus() const
Get DAQ frame status.
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable STRING $argv[2] set_array QUANTILES set_variable FORMULA *[0] exp(-0.5 *(x-[1])*(x-[1])/([2]*[2]))" set_variable MODULE `getModule -a $DETECTOR -L "$STRING 0"` source JAcousticsToolkit.sh typeset -A TRIPODS get_tripods $WORKDIR/tripod.txt TRIPODS XMEAN
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.