12 #include "TProfile2D.h"
35 int main(
int argc,
char **argv)
41 JMultipleFileScanner<JDAQTimeslice> inputFile;
46 JRange<double> totVeto_ns;
47 JRange<int> multiplicity;
54 JParser<> zap(
"Auxiliary program to determine time slewing from K40 data.");
59 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
61 zap[
't'] =
make_field(totVeto_ns) = JRange<double>(0, 1e4);
62 zap[
'M'] =
make_field(multiplicity) = JRange<int>(2, 2);
63 zap[
'C'] =
make_field(ct) = JRange<double>(0.0, 1.0);
69 catch(
const exception &error) {
70 FATAL(error.what() << endl);
77 if (!multiplicity.is_valid()) {
FATAL(
"Invalid multiplicity " << multiplicity << endl); }
78 if ( multiplicity.getLowerLimit() < 2) {
FATAL(
"Invalid multiplicity " << multiplicity << endl); }
79 if (!totVeto_ns .
is_valid()) {
FATAL(
"Invalid time over threshold " << totVeto_ns << endl); }
81 JHit::setSlewing(slewing);
86 load(detectorFile, detector);
88 catch(
const JException& error) {
92 if (detector.empty()) {
93 FATAL(
"Empty detector." << endl);
96 const JModuleRouter router(detector);
101 TH1D h0(
"h0", NULL, 41, -TMax_ns, +TMax_ns);
102 TH1D h1(
"h1", NULL, 50, -0.5, 49.5);
103 TProfile2D hx(
"hx", NULL, 50, -0.5, 49.5, 50, -0.5, 49.5, -2*TMax_ns, +2*TMax_ns);
104 TProfile h2(
"h2", NULL, 100, 0.5, 100.5, -2*TMax_ns, +2*TMax_ns);
112 while (inputFile.hasNext()) {
114 STATUS(
"event: " << setw(10) << inputFile.getCounter() <<
'\r');
DEBUG(endl);
118 for (JDAQTimeslice::const_iterator super_frame = timeslice->begin(); super_frame != timeslice->end(); ++super_frame) {
120 if (router.hasModule(super_frame->getModuleID()) && !super_frame->empty()) {
122 const JModule& module = router.getModule(super_frame->getModuleID());
126 buffer.resize(super_frame->size());
135 sort(buffer.begin(), buffer.end());
141 while (++q != buffer.end() && q->getT() - p->getT() <= TMax_ns ) {}
150 const double dot =
JMATH::getDot(module.getPMT(__p->getPMT()),
151 module.getPMT(__q->getPMT()));
155 h0.Fill(__q->getT() - __p->getT());
156 h1.Fill(__q->getToT());
157 hx.Fill(__p->getToT(), __q->getToT(), __q->getT() - __p->getT());
159 if (totVeto_ns(__p->getToT())) {
160 h2.Fill(__q->getToT(), __q->getT() - __p->getT());
179 TF1 f1(
"f1",
"[0]*exp([1]*sqrt(x) + [2]*x) + [3]");
181 f1.SetParameter(0, h2.GetMaximum());
182 f1.SetParameter(1, -0.01);
183 f1.SetParameter(2, -0.05);
184 f1.SetParameter(3, h2.GetMinimum());
186 h2.ProjectionX()->Fit(&f1,
"",
"same");
188 for (
int i = 0; i != f1.GetNpar(); ++i) {
189 cout <<
"\tstatic double p" << i <<
"() { return " << setw(9) << setprecision(5) << f1.GetParameter(i) <<
"; }" << endl;