Jpp 19.3.0-rc.2
the software that should make you happy
Loading...
Searching...
No Matches
JSlewingK40.cc
Go to the documentation of this file.
1
2#include <string>
3#include <iostream>
4#include <iomanip>
5#include <algorithm>
6#include <random>
7
8#include "TROOT.h"
9#include "TFile.h"
10#include "TH1D.h"
11#include "TH2D.h"
12#include "TF1.h"
13#include "TProfile.h"
14#include "TProfile2D.h"
15
16#include "JROOT/JMinimizer.hh"
17
23#include "JMath/JMathToolkit.hh"
24#include "JTrigger/JHitR0.hh"
26#include "JTools/JRange.hh"
28#include "JSupport/JSupport.hh"
29#include "Jeep/JParser.hh"
30#include "Jeep/JMessage.hh"
31
32
33/**
34 * \file
35 *
36 * Auxiliary program to determine time slewing from K40 data.
37 * \author mdejong
38 */
39int main(int argc, char **argv)
40{
41 using namespace std;
42 using namespace JPP;
43 using namespace KM3NETDAQ;
44
46 JLimit_t& numberOfEvents = inputFile.getLimit();
47 string outputFile;
48 string detectorFile;
49 Double_t TMax_ns;
50 JRange<double> totVeto_ns;
51 JRange<int> multiplicity;
53 bool slewing;
54 int debug;
55
56 try {
57
58 JParser<> zap("Auxiliary program to determine time slewing from K40 data.");
59
60 zap['f'] = make_field(inputFile);
61 zap['o'] = make_field(outputFile) = "monitor.root";
62 zap['a'] = make_field(detectorFile);
63 zap['n'] = make_field(numberOfEvents) = JLimit::max();
64 zap['T'] = make_field(TMax_ns) = 20.5; // [ns]
65 zap['t'] = make_field(totVeto_ns) = JRange<double>(0, 1e4); // time over threshold of reference hits [ns]
66 zap['M'] = make_field(multiplicity) = JRange<int>(2, 2);
67 zap['C'] = make_field(ct) = JRange<double>(0.0, 1.0); // cosine space angle between PMTs
68 zap['S'] = make_field(slewing);
69 zap['d'] = make_field(debug) = 1;
70
71 zap(argc, argv);
72 }
73 catch(const exception &error) {
74 FATAL(error.what() << endl);
75 }
76
77
78 using namespace KM3NETDAQ;
79
80 if (!multiplicity.is_valid()) { FATAL("Invalid multiplicity " << multiplicity << endl); }
81 if ( multiplicity.getLowerLimit() < 2) { FATAL("Invalid multiplicity " << multiplicity << endl); }
82 if (!totVeto_ns .is_valid()) { FATAL("Invalid time over threshold " << totVeto_ns << endl); }
83
84 JHit::setSlewing(slewing);
85
87
88 try {
89 load(detectorFile, detector);
90 }
91 catch(const JException& error) {
92 FATAL(error);
93 }
94
95 if (detector.empty()) {
96 FATAL("Empty detector." << endl);
97 }
98
99 const JModuleRouter router(detector);
100
101
102 TFile out(outputFile.c_str(), "recreate");
103
104 TH1D h0("h0", NULL, 41, -TMax_ns, +TMax_ns);
105 TH1D h1("h1", NULL, 50, -0.5, 49.5);
106 TProfile2D hx("hx", NULL, 50, -0.5, 49.5, 50, -0.5, 49.5, -2*TMax_ns, +2*TMax_ns);
107 TProfile h2("h2", NULL, 100, 0.5, 100.5, -2*TMax_ns, +2*TMax_ns);
108
109 h0.Sumw2();
110 h1.Sumw2();
111
112 random_device rd;
113 mt19937 g(rd());
114
115 vector<JHitR0> buffer;
116
117 while (inputFile.hasNext()) {
118
119 STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
120
121 const JDAQTimeslice* timeslice = inputFile.next();
122
123 for (JDAQTimeslice::const_iterator super_frame = timeslice->begin(); super_frame != timeslice->end(); ++super_frame) {
124
125 if (router.hasModule(super_frame->getModuleID()) && !super_frame->empty()) {
126
127 const JModule& module = router.getModule(super_frame->getModuleID());
128
129 // process data
130
131 buffer.resize(super_frame->size());
132
133 vector<JHitR0>::iterator out = buffer.begin();
134
135 for (JDAQSuperFrame::const_iterator i = super_frame->begin(); i != super_frame->end(); ++i) {
136 *out = JHitR0(i->getPMT(), JHitToolkit<JHit>::getHit(*i, module.getPMT(i->getPMT())));
137 ++out;
138 }
139
140 sort(buffer.begin(), buffer.end());
141
142 for (vector<JHitR0>::iterator p = buffer.begin(); p != buffer.end(); ) {
143
145
146 while (++q != buffer.end() && q->getT() - p->getT() <= TMax_ns ) {}
147
148 if (multiplicity(distance(p,q))) {
149
150 shuffle(p,q,g);
151
152 for (vector<JHitR0>::const_iterator __p = p; __p != q; ++__p) {
153 for (vector<JHitR0>::const_iterator __q = __p; ++__q != q; ) {
154
155 const double dot = JMATH::getDot(module.getPMT(__p->getPMT()),
156 module.getPMT(__q->getPMT()));
157
158 if (ct(dot)) {
159
160 h0.Fill(__q->getT() - __p->getT());
161 h1.Fill(__q->getToT());
162 hx.Fill(__p->getToT(), __q->getToT(), __q->getT() - __p->getT());
163
164 if (totVeto_ns(__p->getToT())) {
165 h2.Fill(__q->getToT(), __q->getT() - __p->getT());
166 }
167 }
168 }
169 }
170
171 sort(p,q);
172 }
173
174 p = q;
175 }
176 }
177 }
178 }
179 STATUS(endl);
180
181
182 // Fit function
183
184 TF1 f1("f1", "[0]*exp([1]*sqrt(x) + [2]*x) + [3]");
185
186 f1.SetParameter(0, h2.GetMaximum());
187 f1.SetParameter(1, -0.01);
188 f1.SetParameter(2, -0.05);
189 f1.SetParameter(3, h2.GetMinimum());
190
191 h2.ProjectionX()->Fit(&f1, "", "same");
192
193 for (int i = 0; i != f1.GetNpar(); ++i) {
194 cout << "\tstatic double p" << i << "() { return " << setw(9) << setprecision(5) << f1.GetParameter(i) << "; }" << endl;
195 }
196
197
198 out.Write();
199 out.Close();
200}
string outputFile
KM3NeT DAQ constants, bit handling, etc.
Data structure for detector geometry and calibration.
Basic data structure for L0 hit.
Tools for handling different hit types.
Auxiliary methods for geometrical methods.
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define STATUS(A)
Definition JMessage.hh:63
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
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
Definition JParser.hh:2142
Auxiliary class to define a range between two values.
int main(int argc, char **argv)
ROOT TTree parameter settings of various packages.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Detector data structure.
Definition JDetector.hh:96
Router for direct addressing of module data in detector data structure.
bool hasModule(const JObjectID &id) const
Has module.
const JModule & getModule(const JObjectID &id) const
Get module parameters.
Data structure for a composite optical module.
Definition JModule.hh:75
const JPMT & getPMT(const int index) const
Get PMT.
Definition JModule.hh:172
General exception.
Definition JException.hh:24
Utility class to parse command line options.
Definition JParser.hh:1698
General purpose class for object reading from a list of file names.
virtual bool hasNext() override
Check availability of next element.
counter_type getCounter() const
Get counter.
virtual const pointer_type & next() override
Get next element.
Range of values.
Definition JRange.hh:42
bool is_valid() const
Check validity of range.
Definition JRange.hh:311
T getLowerLimit() const
Get lower limit.
Definition JRange.hh:202
Reduced data structure for L0 hit.
Definition JHitR0.hh:27
static void setSlewing(const bool slewing)
Set slewing option.
Hit data structure.
Definition JDAQHit.hh:35
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
double getDot(const JFirst_t &first, const JSecond_t &second)
Get dot product of objects.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
bool is_valid(const json &js)
Check validity of JSon data.
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
Detector file.
Definition JHead.hh:227
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
Template definition of hit toolkit.