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)
95 JLimit_t& numberOfEvents = inputFile.getLimit();
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'; }
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);
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);
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);
Data structure for detector geometry and calibration.
Basic data structure for L0 hit.
General purpose messaging.
#define DEBUG(A)
Message macros.
Direct access to module in detector data structure.
Scanning of objects from multiple files according a format that follows from the extension of each fi...
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Utility class to parse parameter values.
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
int main(int argc, char **argv)
Auxiliary class to define a range between two values.
ROOT TTree parameter settings of various packages.
Basic data structure for time and time over threshold information of hit.
Router for direct addressing of module data in detector data structure.
const JModule & getModule(const JObjectID &id) const
Get module parameters.
Utility class to parse parameter values.
Auxiliary class for multiplexing object iterators.
virtual bool hasNext() override
Check availability of next element.
virtual const pointer_type & next() override
Get next element.
Utility class to parse command line options.
General purpose class for object reading from a list of file names.
Template specialisation of L0 builder for JHitR0 data type.
const JPolynome f1(1.0, 2.0, 3.0)
Function.
std::string getLabel(const JLocation &location)
Get module label for monitoring and other applications.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
bool putObject(TDirectory &dir, const TObject &object)
Write object to ROOT directory.
Long64_t counter_type
Type definition for counter.
KM3NeT DAQ data structures and auxiliaries.
double getFrameTime()
Get frame time duration.
double getTimeOfFrame(const int frame_index)
Get start time of frame in ns since start of run for a given frame index.
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
std::map< int, range_type > map_type
static const int PMT_DISABLE
KM3NeT Data Definitions v3.4.0-8-ge14cb17 https://git.km3net.de/common/km3net-dataformat.
Auxiliary data structure for sequence of same character.
Auxiliary data structure for floating point format specification.
Type definition of range.
Auxiliary class for TDC constraints.
range_type equal_range(const int id)
Get range of constraints for given module.
bool is_valid(const bool option=false) const
Check validity of TDC constrants.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Auxiliary class to select ROOT class based on class name.
Auxiliary class for defining the range of iterations of objects.