10 #include "TFitResult.h"
53 inline double get_time(
double t1,
const double T)
55 for (
int buffer[] = { 1000, 100, 10, 1, 0 }, *
i = buffer; *
i != 0; ++
i) {
57 const double xmin = -0.5 * (*i) *
T;
58 const double xmax = +0.5 * (*i) *
T;
60 while (t1 > xmax) { t1 -= (*i) *
T; }
61 while (t1 < xmin) { t1 += (*i) *
T; }
86 int main(
int argc,
char **argv)
90 using namespace KM3NETDAQ;
99 double laserFrequency_Hz;
100 bool overwriteDetector;
114 JParser<> zap(
"Application for dark room time calibration.");
116 zap[
'f'] =
make_field(inputFile,
"input file (time slice data from laser calibration).");
118 zap[
'a'] =
make_field(detectorFile,
"detector file.");
120 "Set reference PMTs, e.g."
121 "\n-! \"808969848 0 808982077 23\" sets PMT 0 of module 808969848 and PMT 23 of module 808982077 as references.");
122 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
123 zap[
'l'] =
make_field(laserFrequency_Hz,
"laser frequency [Hz]") = 10000;
124 zap[
'A'] =
make_field(overwriteDetector,
"overwrite detector file provided through '-a' with correct time offsets.");
125 zap[
'O'] =
make_field(option,
"ROOT fit option, see TH1::Fit.") =
"LS";
127 zap[
'C'] =
make_field(selector,
"timeslice selector, e.g. JDAQTimesliceL1.") = getROOTClassSelection<JDAQTimesliceTypes_t>();
128 zap[
'T'] =
make_field(T_ns,
"time window for time-over-threshold monitor") =
JRange_t(-10.0, +10.0);
134 catch(
const exception& error) {
135 FATAL(error.what() << endl);
141 if (laserFrequency_Hz <= 0.0) {
142 FATAL(
"Invalid laser frequency " << laserFrequency_Hz << endl);
145 const double laserPeriod_ns = 1.0e9 / laserFrequency_Hz;
147 if (option.find(
'R') == string::npos) { option +=
'R'; }
148 if (option.find(
'S') == string::npos) { option +=
'S'; }
174 const double xmin = (
X ==
JRange_t() ? -0.5 * laserPeriod_ns :
X.getLowerLimit());
175 const double xmax = (
X ==
JRange_t() ? +0.5 * laserPeriod_ns :
X.getUpperLimit());
176 const int nx = 2 * (int) (
xmax -
xmin);
178 TH1D h0(
"h0", NULL, nx,
xmin,
xmax);
179 TH1D h1(
"h1", NULL, 256, -0.5, +255.5);
183 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ++module) {
187 for (JTDC_t::const_iterator
i = range.first;
i != range.second; ++
i) {
190 <<
" string " << setw(3) << module->getString()
191 <<
" floor " << setw(2) << module->getFloor()
192 <<
" module " << setw(8) << module->getID()
193 <<
" channel " << setw(2) <<
i->second << endl);
201 zmap.insert(make_pair(
id,
new TH1D(os.str().c_str(), NULL, nx,
xmin,
xmax)));
215 for ( ; in.
hasNext() && counter != inputFile.getLimit(); ++counter) {
217 STATUS(
"event: " << setw(10) << counter <<
'\r');
DEBUG(endl);
221 for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
225 if (range.first != range.second) {
227 const double t0 = get_time(
getTimeOfFrame(frame->getFrameIndex()), laserPeriod_ns);
231 buildL0(*frame, moduleRouter.
getModule(frame->getModuleID()), back_inserter(dataL0));
233 for (JDataL0_t::const_iterator hit = dataL0.begin(); hit != dataL0.end(); ++hit) {
237 map_type::const_iterator p = zmap.find(
id);
239 if (p != zmap.end()) {
241 const double t1 = get_time(t0 + hit->getT(), laserPeriod_ns);
249 h1.Fill(hit->getToT());
264 TF1
f1(
"f1",
"[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))/(TMath::Sqrt(2.0*TMath::Pi())*[2]) + [3]");
266 for (map_type::iterator it = zmap.begin(); it != zmap.end(); ++it) {
269 TH1D* h1 = it->second;
271 if (h1->GetEntries() == 0) {
272 WARNING(
"Histogram " << h1->GetName() <<
" empty" << endl);
276 STATUS(
"--- PMT = " << pmt <<
"; histogram " << h1->GetName() << endl);
282 Double_t
sigma = 2.0;
285 for (
int i = 1;
i != h1->GetNbinsX(); ++
i) {
287 const Double_t
x = h1->GetBinCenter (
i);
288 const Double_t
y = h1->GetBinContent(
i);
296 f1.SetParameter(0, ymax/sigma);
297 f1.SetParameter(1, x0);
298 f1.SetParameter(2, sigma);
299 f1.SetParameter(3, ymin);
301 for (Int_t
i = 0;
i != f1.GetNpar(); ++
i) {
302 f1.SetParError(
i, 0.0);
308 TFitResultPtr
result = h1->Fit(&f1, option.c_str(),
"same", x0 - 3 *
sigma, x0 + 3 *
sigma);
310 if (result.Get() == NULL) {
311 FATAL(
"Invalid TFitResultPtr " << h1->GetName() << endl);
314 if (
debug >=
notice_t || !result->IsValid() || f1.GetParameter(0) < Wmin) {
315 cout <<
"Histogram " << h1->GetName() <<
" fit " << (result->IsValid() ?
"ok" :
"failed") << endl;
316 cout <<
"\tw = " <<
FIXED(12,3) << f1.GetParameter(0) << endl;
317 cout <<
"\tx0 = " <<
FIXED(12,3) << x0 << endl;
318 cout <<
"\tt0 = " <<
FIXED(12,3) << f1.GetParameter(1) << endl;
324 int number_of_peaks = 0;
326 Double_t dx = 2.0 * f1.GetParameter(2);
328 Double_t
Y = f1.GetParameter(3);
334 for (Int_t il = 1, ir = il; ir <= nx; ) {
336 for ( ; ir <= nx && h1->GetBinCenter(ir) <= h1->GetBinCenter(il) + dx; ++ir) {
337 W += h1->GetBinContent(ir) -
Y;
342 number_of_peaks += 1;
350 W -= h1->GetBinContent(il) -
Y;
355 if (number_of_peaks != 1) {
356 WARNING(
"Number of peaks " << h1->GetName() <<
' ' << number_of_peaks <<
" != 1" << endl);
359 if (result->IsValid() && f1.GetParameter(0) >= Wmin) {
360 t0[pmt] = f1.GetParameter(1);
370 const double W = laserFrequency_Hz * counter *
getFrameTime() * 1.0e-9;
372 NOTICE(
"Expection values [npe]" << endl);
375 NOTICE(
i->first <<
' ' <<
FIXED(7,3) <<
i->second / W << endl);
380 if (overwriteDetector) {
384 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ++module) {
388 if (range.first != range.second) {
394 if (p0 != t0.end()) {
396 const double t1 = p0->second;
403 module->getPMT(pmt).subT0(p1->second);
405 module->getPMT(pmt).subT0(t1);
410 if (!module->getPMT(
id.getPMTAddress()).has(
PMT_DISABLE)) {
415 << setw(10) << module->getID() <<
"/" <<
FILL(2,
'0') <<
id.getPMTAddress() <<
FILL() <<
' '
416 <<
"missing or insufficient signal." << endl);
423 NOTICE(
"Store calibration data on file " << detectorFile << endl);
431 FATAL(
"Number of errors " << errors <<
" != 0" << endl);
Utility class to parse command line options.
int main(int argc, char *argv[])
ROOT TTree parameter settings of various packages.
Basic data structure for L0 hit.
const JModule & getModule(const JObjectID &id) const
Get module parameters.
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
std::string getLabel(const JLocation &location)
Get module label for monitoring and other applications.
Auxiliary class to select ROOT class based on class name.
Router for direct addressing of module data in detector data structure.
Auxiliary class for TDC constraints.
Utility class to parse parameter values.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
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
Auxiliary data structure for floating point format specification.
Auxiliary class for multiplexing object iterators.
Basic data structure for time and time over threshold information of hit.
Data structure for detector geometry and calibration.
Utility class to parse parameter values.
double getTimeOfFrame(const int frame_index)
Get start time of frame in ns since start of run for a given frame index.
const JPolynome f1(1.0, 2.0, 3.0)
Function.
Auxiliary class for defining the range of iterations of objects.
Type definition of range.
Template specialisation of L0 builder for JHitR0 data type.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
double getFrameTime()
Get frame time duration.
do set_variable OUTPUT_DIRECTORY $WORKDIR T
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
static const int PMT_DISABLE
KM3NeT Data Definitions v3.3.0-2-g5cc95cf https://git.km3net.de/common/km3net-dataformat.
virtual const pointer_type & next() override
Get next element.
virtual bool hasNext() override
Check availability of next element.
bool putObject(TDirectory &dir, const TObject &object)
Write object to ROOT directory.
General purpose messaging.
Auxiliary data structure for sequence of same character.
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.
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Auxiliary class to define a range between two values.
General purpose class for object reading from a list of file names.
then fatal The output file must have the wildcard in the e g root fi 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
Utility class to parse command line options.
no fit printf nominal n $STRING awk v X
const JLimit & getLimit() const
Get limit.
std::map< int, range_type > map_type
do set_variable DETECTOR_TXT $WORKDIR detector
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
#define DEBUG(A)
Message macros.