54 static const char* rates_t =
"rates";
55 static const char* count_t =
"counts";
56 static const char* tails_t =
"tails";
57 static const char* rndms_t =
"randoms";
62 static const double ToTmax_ns = std::numeric_limits<KM3NETDAQ::JDAQHit::JTOT_t>::max();
104 return (h2s != NULL &&
140 int main(
int argc,
char **argv)
144 using namespace KM3NETDAQ;
163 JParser<> zap(
"Auxiliary program to determine PMT parameters from K40 data.");
165 zap[
'f'] =
make_field(inputFile,
"input file.");
167 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
168 zap[
'a'] =
make_field(detectorFile,
"detector file.");
169 zap[
'T'] =
make_field(TMax_ns,
"time window [ns].") = 20.0;
172 zap[
'b'] =
make_field(background,
"background estimation method.") = rates_t, count_t, tails_t, rndms_t;
175 zap[
'D'] =
make_field(deadTime_us,
"L1 dead time (us)") = 0.0;
176 zap[
'C'] =
make_field(selector,
"timeslice selector, e.g. JDAQTimesliceL1.") = getROOTClassSelection<JDAQTimesliceTypes_t>();
177 zap[
'O'] =
make_field(option,
"hit pre-processing option.") = JPreprocessor::getOptions();
182 catch(
const exception &error) {
183 FATAL(error.what() << endl);
192 if (selector == JDAQTimesliceL2::Class_Name() ||
193 selector == JDAQTimesliceSN::Class_Name()) {
194 FATAL(
"Option -C <selector> " << selector <<
" not compatible with calibration method." << endl);
197 if (selector == JDAQTimeslice ::Class_Name() ||
198 selector == JDAQTimesliceL1::Class_Name()) {
206 FATAL(
"No trigger parameters from input:" << error.
what() << endl);
209 if ((selector == JDAQTimeslice ::Class_Name() && parameters.writeL1.prescale > 0) ||
210 (selector == JDAQTimesliceL1::Class_Name())) {
212 if (parameters.TMaxLocal_ns < TMax_ns) {
213 FATAL(
"Option -T <TMax_ns> = " << TMax_ns <<
" is larger than in the trigger " << parameters.TMaxLocal_ns << endl);
216 if (background == rndms_t ||
217 background == count_t) {
218 FATAL(
"Option -C <selector> " << selector <<
" incompatible with option -b <background> " << background << endl);
223 if (!multiplicity.is_valid() || multiplicity.getLowerLimit() < 2) {
224 FATAL(
"Invalid option -M <multiplicity> " << multiplicity << endl);
227 if (!totRange_ns.is_valid()) {
228 FATAL(
"Invalid option -t <totRange_ns> " << totRange_ns << endl);
231 if (background == tails_t) {
233 if (!Tail_ns.is_valid()) {
234 FATAL(
"Invalid option -B <Tail_ns> " << Tail_ns << endl);
237 if (Tail_ns.getUpperLimit() > TMax_ns) {
238 FATAL(
"Invalid option -B <Tail_ns> " << Tail_ns <<
"; upper limit larger than option -T <TMax_ns> " << TMax_ns << endl);
256 FATAL(
"Empty detector." << endl);
268 t0.push_back(i * 2 * TMax_ns);
278 const double RTU = (cpu.getIntegralOfChargeProbability(cpu.getNPE(totRange_ns.getLowerLimit()), cpu.getNPE(totRange_ns.getUpperLimit()), NPE)
280 cpu.getIntegralOfChargeProbability(cpu.getNPE(0.0), cpu.getNPE(ToTmax_ns), NPE));
293 const double xmin = -0.5;
294 const double xmax = nx - 0.5;
296 const double ymin = -floor(TMax_ns) + 0.5;
297 const double ymax = +floor(TMax_ns) - 0.5;
298 const int ny = (int) ((ymax - ymin) / 1.0);
301 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ++module) {
305 if (!module->empty()) {
307 NOTICE(
"Booking histograms for module " << module->getID() << endl);
322 const double deadTime_ns = deadTime_us * 1e3;
331 for (JMultipleFileScanner_t::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
333 STATUS(
"Processing: " << *i << endl);
342 STATUS(
"event: " << setw(10) << counter <<
'\r');
DEBUG(endl);
354 if (background == rates_t) {
360 ERROR(
"Frame indices do not match at "
361 <<
"[counter = " << counter <<
"]"
368 FATAL(
"ROOT TTrees not aligned; Run number " << timeslice->
getRunNumber() << endl);
373 for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
375 if (router.
hasModule(frame->getModuleID())) {
386 if (background == count_t) {
392 }
else if (background == tails_t) {
396 }
else if (background == rndms_t) {
400 }
else if (background == rates_t) {
407 rate_Hz[i] = RTU * sum.
getRate(i);
412 FATAL(
"Summary frame of module " << frame->getModuleID() <<
" not found at frame index " << timeslice->
getFrameIndex() << endl);
438 TH2D* h2s = zmap[address.
first].h2s;
439 TH1D* h1b = zmap[address.
first].h1b;
440 TH1D* h1L = zmap[address.
first].h1L;
446 JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*frame, module);
448 buffer.preprocess(option, match);
450 for (JSuperFrame2D_t::iterator i = buffer.begin(); i != buffer.end(); ++i) {
452 const int pmt = i->getPMTAddress();
454 if (veto[pmt] || !rateRange_Hz(rate_Hz[pmt])) {
459 JSuperFrame1D_t& data = JSuperFrame1D_t::multiplex(buffer);
461 DEBUG(
"Number of hits " << timeslice->
getFrameIndex() <<
":" << frame->getModuleID() <<
' ' << frame->size() <<
' ' << data.size() << endl);
470 size_t numberOfSignalEvents = 0;
472 double t1 = numeric_limits<double>::lowest();
478 while (++q != data.end() && q->getT() - p->getT() <= TMax_ns) {}
482 if (multiplicity(N)) {
484 numberOfSignalEvents += 1;
486 if (p->getT() > t1 + deadTime_ns) {
491 if (totRange_ns( __p->getToT()) && totRange_ns(__q->getToT())) {
492 h2s->Fill((
double) combinatorics.
getIndex(__p->getPMT(),__q->getPMT()),
493 JCombinatorics::getSign(__p->getPMT(),__q->getPMT()) * (__q->getT() - __p->getT()));
508 if (background == rndms_t) {
511 *i =
hit_type(i->getPMT(),
JHit(i->getT() + t0[i->getPMT()], i->getToT()));
514 sort(data.begin(), data.end());
516 double t1 = numeric_limits<double>::lowest();
522 while (++q != data.end() && q->getT() - p->getT() <= TMax_ns) {}
526 if (multiplicity(N)) {
528 if (p->getT() > t1 + deadTime_ns) {
533 if (totRange_ns(__p->getToT()) && totRange_ns(__q->getToT())) {
534 h1b->Fill((
double) combinatorics.
getIndex(p->getPMT(),q->getPMT()), 1.0);
546 }
else if (background == rates_t ||
547 background == count_t) {
553 const Double_t
R1 = rate_Hz[i];
555 if (!veto[i] && rateRange_Hz(R1)) {
563 const Double_t
R1 = rate_Hz[i];
564 const Double_t R2 = rate_Hz[
j];
566 if (!veto[i] && rateRange_Hz(R1) &&
567 !veto[
j] && rateRange_Hz(R2)) {
570 const double ED = (Rs - rate_Hz[i] - rate_Hz[
j]) * 2 * TMax_ns * 1e-9;
571 const double E1 = rate_Hz[i] * 2 * TMax_ns * 1e-9;
572 const double E2 = rate_Hz[
j] * 2 * TMax_ns * 1e-9;
576 multiplicity.getLowerLimit(),
577 multiplicity.getUpperLimit()) *
getFrameTime() / (2*TMax_ns);
579 h1b->Fill((
double) combinatorics.
getIndex(i,
j),
N);
593 const int pmt1 = combinatorics.
getPair(i).first;
594 const int pmt2 = combinatorics.
getPair(i).second;
596 if (!veto[pmt1] && rateRange_Hz(rate_Hz[pmt1]) &&
597 !veto[pmt2] && rateRange_Hz(rate_Hz[pmt2]) ) {
599 h1L->Fill(i, livetime_s);
608 if (background == tails_t ) {
612 if (hist->is_valid()) {
614 TH2D* h2s = hist->h2s;
615 TH1D* h1b = hist->h1b;
625 for (
int l = 1; l <= h2s->GetYaxis()->GetNbins(); ++l) {
627 const Double_t dt = h2s->GetYaxis()->GetBinCenter(l) ;
629 if (Tail_ns(fabs(dt))) {
630 Y += h2s->GetBinContent(k+1,l);
631 W += h2s->GetBinError(k+1,l) * h2s->GetBinError(k+1,l);
636 h1b->SetBinContent(k+1, Y / N * (2.0*TMax_ns) / ((ymax - ymin)/ny) );
637 h1b->SetBinError (k+1, sqrt(W/N) * (2.0*TMax_ns) / ((ymax - ymin)/ny) );
651 weights_hist.GetXaxis()->SetBinLabel(2,
WS_t);
652 weights_hist.GetXaxis()->SetBinLabel(3,
WB_t);
655 weights_hist.Fill(
WS_t, (ymax - ymin)/ny);
656 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.
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.
then JShowerPostfit f $INPUT_FILE o $OUTPUT_FILE N
*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 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 JAcoustics.sh -- typeset -A TRIPODS get_tripods $WORKDIR/tripod.txt TRIPODS XMEAN
int getFrameIndex() const
Get frame index.
1-dimensional frame with time calibrated data from one optical module.
then JCookie sh JDataQuality D $DETECTOR_ID R $RUNS[*] o $QUALITY_TXT d $DEBUG!fi fi JDataQuality f $QUALITY_TXT Q livetime_s
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.
bool is_valid(const json &js)
Check validity of JSon data.
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.
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.
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.
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
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.
double coincidenceP(double E1, double E2, double ED, int M_min, int M_max)
Coincidence probability of two PMTs within one module.
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 source JAcoustics sh $DETECTOR_ID CHECK_EXIT_CODE typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Auxiliary class for specifying the way of pre-processing of hits.
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.