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 = 250;
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) {
450 FATAL(
"Summary frame of module " << frame->getModuleID() <<
" not found at frame index " << timeslice->
getFrameIndex() << endl);
462 !rateRange_Hz(rate_Hz[
i])) {
468 const JCombinatorics& combinatorics = histogram.getCombinatorics();
470 TH2D* h2s = histogram.h2s;
471 TH1D* h1b = histogram.h1b;
472 TH1D* h1L = histogram.h1L;
478 JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*frame, module);
480 buffer.preprocess(option, match);
482 for (JSuperFrame2D_t::iterator
i = buffer.begin();
i != buffer.end(); ++
i) {
483 if (veto[
i->getPMTAddress()]) {
488 JSuperFrame1D_t&
data = JSuperFrame1D_t::multiplex(buffer);
490 DEBUG(
"Number of hits " << timeslice->
getFrameIndex() <<
":" << frame->getModuleID() <<
' ' << frame->size() <<
' ' << data.size() << endl);
495 size_t numberOfSignalEvents = 0;
497 double t1 = numeric_limits<double>::lowest();
503 while (++q != data.end() && q->getT() - p->getT() <= TMax_ns) {}
507 if (multiplicity(N)) {
509 numberOfSignalEvents += 1;
511 if (p->getT() > t1 + deadTime_ns) {
516 if (totRange_ns( __p->getToT()) && totRange_ns(__q->getToT())) {
517 h2s->Fill((
double) combinatorics.getIndex(__p->getPMT(),__q->getPMT()),
518 JCombinatorics::getSign(__p->getPMT(),__q->getPMT()) * (__q->getT() - __p->getT()));
533 if (background == rndms_t) {
539 sort(data.begin(), data.end());
541 double t1 = numeric_limits<double>::lowest();
547 while (++q != data.end() && q->getT() - p->getT() <= TMax_ns) {}
551 if (multiplicity(N)) {
553 if (p->getT() > t1 + deadTime_ns) {
558 if (totRange_ns(__p->getToT()) && totRange_ns(__q->getToT())) {
559 h1b->Fill((
double) combinatorics.getIndex(p->getPMT(),q->getPMT()), 1.0);
571 }
else if (background == rates_t ||
572 background == count_t) {
585 if (!veto[
i] && !veto[
j]) {
587 const double R1 = rate_Hz[
i];
588 const double R2 = rate_Hz[
j];
595 const double N =
getP((R1) * 2 * TMax_ns * 1e-9,
596 (R2) * 2 * TMax_ns * 1e-9,
597 (Rs - R1 - R2) * 2 * TMax_ns * 1e-9,
598 multiplicity.getLowerLimit(),
599 multiplicity.getUpperLimit()) *
getFrameTime() / (2*TMax_ns);
601 h1b->Fill((
double) combinatorics.getIndex(
i,j),
N);
613 for (
size_t i = 0;
i < combinatorics.getNumberOfPairs(); ++
i) {
615 const int pmt1 = combinatorics.getPair(
i).first;
616 const int pmt2 = combinatorics.getPair(
i).second;
618 if (!veto[pmt1] && !veto[pmt2]) {
619 h1L->Fill(
i, livetime_s);
628 if (background == tails_t ) {
632 if (hist->is_valid()) {
636 TH2D* h2s = hist->h2s;
637 TH1D* h1b = hist->h1b;
647 for (
int l = 1; l <= h2s->GetYaxis()->GetNbins(); ++l) {
649 const Double_t dt = h2s->GetYaxis()->GetBinCenter(l) ;
651 if (Tail_ns(fabs(dt))) {
652 Y += h2s->GetBinContent(k+1,l);
653 W += h2s->GetBinError(k+1,l) * h2s->GetBinError(k+1,l);
658 h1b->SetBinContent(k+1, Y / N * (2.0*TMax_ns) / ((ymax - ymin)/ny) );
659 h1b->SetBinError (k+1, sqrt(W/N) * (2.0*TMax_ns) / ((ymax - ymin)/ny) );
673 weights_hist.GetXaxis()->SetBinLabel(2,
WS_t);
674 weights_hist.GetXaxis()->SetBinLabel(3,
WB_t);
677 weights_hist.Fill(
WS_t, (ymax - ymin)/ny);
678 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 warning Cannot perform comparison test for histogram
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.
double getValue(const int tdc) const
Get value.
static const char *const _2S
Name extension for 2D counts.
#define DEBUG(A)
Message macros.