Jpp test-rotations-old
the software that should make you happy
Loading...
Searching...
No Matches
JMonitorL1dt.cc
Go to the documentation of this file.
10#include "JTrigger/JHitL1.hh"
11#include "JTrigger/JBuildL2.hh"
15#include "JSupport/JSupport.hh"
18#include "Jeep/JParser.hh"
19#include "JROOT/JManager.hh"
21
22#include "TFile.h"
23#include "TH1D.h"
24#include "TH2D.h"
25#include "TMath.h"
26
27#include <string>
28#include <iostream>
29
30
31using namespace JTOOLS;
32namespace MONITORL1DT {
33
34 /**
35 * Data structure for hit time and DOM identifier.
36 */
37 class JElement {
38 public:
40 id(0),
41 t (0.0)
42 {}
43
44 JElement(const int __id,
45 const double __t) :
46 id(__id),
47 t (__t)
48 {}
49
50 int id;
51 double t;
52 };
53
54
55 /**
56 * Sort JElement in time.
57 *
58 * \param first first element
59 * \param second second element
60 * \return true if second later than first; else false
61 */
62 inline bool operator<(const JElement& first, const JElement& second)
63 {
64 return first.t < second.t;
65 }
66
67
68 /**
69 * Auxiliary data structure for histogram management.
70 */
71 struct JHistogram {
72 /**
73 * Default constructor.
74 */
76 h2s(NULL),
77 h1l(NULL)
78 {}
79
80
81 /**
82 * Constructor.
83 *
84 * \param __h2s 2D histogram for signal
85 * \param __h1l 1D histogram for background
86 */
87 JHistogram(TH2D* __h2s,
88 TH1D* __h1l) :
89 h2s(__h2s),
90 h1l(__h1l)
91 {
92 h2s->Sumw2();
93 h1l->Sumw2();
94 }
95
96 TH2D* h2s;
97 TH1D* h1l;
98 };
99
100 // Taken from C. Pieterse nanobeacon code
102 inline bool comparepair(const pair_type& A, const pair_type& B) {
103 return( abs(A.first - A.second) > abs(B.first - B.second) ) ;
104 }
105}
106
107
108/**
109 * \file
110 *
111 * Auxiliary program to plot L1 time difference correlations between DOMs.
112 *
113 * \author kmelis, lnauta
114 */
115int main(int argc, char **argv)
116{
117 using namespace std;
118 using namespace MONITORL1DT;
119 using namespace JPP;
120
122 string outputFile;
123 string detectorFile;
124 double Timewindow_ns;
125 int binwidth; // since the data is taken on a ns level, having bins of non-integer value is non-sensical
126 double TmaxL1_ns;
127 JROOTClassSelector selector;
128 unsigned int multiplicity;
129 bool correct_time;
130 double livetime_s;
131 int debug;
132
133 try {
134
135 JParser<> zap("Program to create L1 hit time difference histograms from raw data.");
136
137 zap['f'] = make_field(inputFile, "input file");
138 zap['o'] = make_field(outputFile, "output file") = "monitor.root";
139 zap['a'] = make_field(detectorFile, "detector file");
140 zap['t'] = make_field(TmaxL1_ns, "max time between L1 hits [ns]") = 1000.0;
141 zap['T'] = make_field(Timewindow_ns, "time window around t=0 [ns]") = 2400.0;
142 zap['w'] = make_field(binwidth, "binwidth [ns]") = 1;
143 zap['C'] = make_field(selector, "datastream selector") = getROOTClassSelection<JDAQTimesliceTypes_t>();
144 zap['m'] = make_field(multiplicity, "minimal multiplicity of the L1 hits") = 2;
145 zap['c'] = make_field(correct_time, "subtract expected arrival time from delta-t");
146 zap['L'] = make_field(livetime_s, "livetime of the data, set to positive value") = -1.0; // set to positive value if you want to set the livetime manually (i.e. MC)
147 zap['d'] = make_field(debug) = 1;
148
149 if (zap.read(argc, argv) != 0)
150 return 1;
151 }
152 catch(const exception &error) {
153 FATAL(error.what() << endl);
154 }
155
156
158
159 try {
160 load(detectorFile, detector);
161 }
162 catch(const JException& error) {
163 FATAL(error);
164 }
165
166 if (detector.empty()) {
167 FATAL("Empty detector." << endl);
168 }
169
170 const JModuleRouter router(detector);
171
172 typedef vector<JHitL1> JFrameL1_t;
173
174 const double ctMin = -1; // full angle over the DOM
175 const JBuildL2<JHitL1> buildL2(JL2Parameters(multiplicity, TmaxL1_ns, ctMin));
176
177 // Manager histograms for single-DOM output
178 vector<JHistogram> zmap(detector.size());
179
180 const int nx = detector.size();
181 const double xmin = -0.5;
182 const double xmax = nx - 0.5;
183
184 const double ymin = -floor(Timewindow_ns) + 0.5;
185 const double ymax = +floor(Timewindow_ns) + 0.5; // changed - to + to get an even number for Nbins, rebin needs numbers that have divisors.
186 const int ny = (int) ((ymax - ymin) / binwidth);
187
188 // Manage histograms for single-DU output
189 Int_t num_of_floors = getNumberOfFloors(detector);
190 JCombinatorics c(num_of_floors);
191 Int_t npairs = c.getNumberOfPairs();
193
194 JManager<int, TH2D>* manager;
195 manager = new JManager <int, TH2D>(new TH2D("h%", "", npairs, 0.5, npairs+0.5, ny, ymin, ymax));
196
197 NOTICE("Running JMonitorL1dt: Monitoring L1 time differences and creating histograms." << endl);
198 for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
199
200 const JModuleAddress& address = router.getAddress(module->getID());
201
202 STATUS("Booking histograms for module " << module->getID() << endl);
203
204 const JString title(module->getID());
205 JString titleString1D;
206 JString titleString2D;
207 titleString1D = title + ".1L";
208 titleString2D = title + ".2S";
209
210 zmap[address.first] = JHistogram(new TH2D((titleString2D).c_str(), NULL, nx, xmin, xmax, ny, ymin, ymax),
211 new TH1D((titleString1D).c_str(), NULL, nx, xmin, xmax));
212
213 for (JDetector::iterator mod = detector.begin(); mod != detector.end(); ++mod) {
214 zmap[address.first].h2s->GetXaxis()->SetBinLabel(distance(detector.begin(), mod)+1, Form("%i", mod->getID()));
215 zmap[address.first].h1l->GetXaxis()->SetBinLabel(distance(detector.begin(), mod)+1, Form("%i", mod->getID()));
216 }
217 }
218
219
221
223
224 int counter = 0;
225
226 for ( ; in.hasNext() && counter != inputFile.getLimit(); ++counter) {
227
228 STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl);
229
230 const JDAQTimeslice* timeslice = in.next();
231
232 JFrameL1_t frameL1;
234 vector<bool> DOM_OK(detector.size(), true);
235
236 for (JDAQTimeslice::const_iterator super_frame = timeslice->begin(); super_frame != timeslice->end(); ++super_frame) {
237 if (router.hasModule(super_frame->getModuleID()) && !super_frame->empty()) {
238
239 const JModuleAddress& address = router.getAddress(super_frame->getModuleID());
240 const JModule& module = detector.getModule(address);
241
242 frameL1.clear();
243
244 buildL2(*super_frame, module, back_inserter(frameL1));
245
246 for (JFrameL1_t::iterator L1hit = frameL1.begin(); L1hit != frameL1.end(); ++L1hit) {
247 buffer.push_back(JElement(address.first, L1hit->begin()->getT()));
248 }
249 }
250 }
251
252 for (vector<JHistogram>::iterator h1 = zmap.begin(); h1 != zmap.end(); ++h1) {
253 if (!DOM_OK[distance(zmap.begin(), h1)]) {
254 continue;
255 }
256 for (unsigned int i = 0; i < detector.size(); ++i) {
257 if (DOM_OK[i]) { h1->h1l->Fill(i, getFrameTime() * 1e-9); } // fill the th1 with the frametime (how long did it measure) in seconds, all weight will give total measured time
258 }
259 }
260
261 // fill histograms with correlations
262 sort(buffer.begin(), buffer.end());
263
264 for (vector<JElement>::const_iterator p = buffer.begin(); p != buffer.end(); ) {
265 vector<JElement>::const_iterator q = p;
266
267 const JModule& module_p = detector.getModule((JModuleAddress)p->id);
268
269 while (++q != buffer.end() && q->t - p->t <= Timewindow_ns ) {
270 // if (p->id == q->id) { continue; }
271
272 const JModule& module_q = detector.getModule((JModuleAddress)q->id);
273
274 double dom_distance = getDistance(module_p.getPosition(), module_q.getPosition());
275 double time_correction = (correct_time ? (dom_distance / getSpeedOfLight()) : 0);
276
277 zmap[p->id].h2s->Fill(q->id, q->t - p->t - time_correction);
278 zmap[q->id].h2s->Fill(p->id, p->t - q->t + time_correction);
279
280 if (module_p.getString() == module_q.getString()) {
281 // indices are pairs running from 0-17 between 0-152, while floors are 1-18 and bins are 1-153
282 int xbin = c.getIndex(module_p.getFloor() - 1, module_q.getFloor() - 1) + 1;
283 (*manager)[module_p.getString()]->Fill(xbin, q->t - p->t - time_correction);
284 }
285 }
286 p++; // move to the next L1 hit
287 }
288 buffer.clear();
289 }
290 STATUS(endl);
291
292 // livetime in case set manually
293 if (livetime_s > 0.0) {
294 for (vector<JHistogram>::iterator i = zmap.begin(); i != zmap.end(); ++i) {
295 TH1D* hl = i->h1l;
296 for (int ibin = 1; ibin <= hl->GetNbinsX(); ++ibin) {
297 hl->SetBinContent(ibin, livetime_s);
298 hl->SetBinError(ibin, 0.0000001);
299 }
300 }
301 }
302
303 TFile out(outputFile.c_str(), "recreate");
304
305 for (vector<JHistogram>::iterator i = zmap.begin(); i != zmap.end(); ++i) {
306 i->h2s->Write();
307 i->h1l->Write();
308 }
309
310 manager->Write(out);
311 out.Write();
312 out.Close();
313}
string outputFile
KM3NeT DAQ constants, bit handling, etc.
Data structure for detector geometry and calibration.
Basic data structure for L1 hit.
Dynamic ROOT object management.
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define STATUS(A)
Definition JMessage.hh:63
#define NOTICE(A)
Definition JMessage.hh:64
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
Direct access to module in detector data structure.
int main(int argc, char **argv)
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
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
int getFloor() const
Get floor number.
Definition JLocation.hh:146
int getString() const
Get string number.
Definition JLocation.hh:135
Address of module in detector data structure.
int first
index of module in detector data structure
Router for direct addressing of module data in detector data structure.
bool hasModule(const JObjectID &id) const
Has module.
const JModuleAddress & getAddress(const JObjectID &id) const
Get address of module.
Data structure for a composite optical module.
Definition JModule.hh:75
const JPosition3D & getPosition() const
Get position.
General exception.
Definition JException.hh:24
Auxiliary class for multiplexing object iterators.
virtual bool hasNext() override
Check availability of next element.
virtual const pointer_type & next() override
Get next element.
Wrapper class around STL string class.
Definition JString.hh:29
Utility class to parse command line options.
Definition JParser.hh:1698
int read(const int argc, const char *const argv[])
Parse the program's command line options.
Definition JParser.hh:1992
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys.
Definition JManager.hh:47
General purpose class for object reading from a list of file names.
Template definition for direct access of elements in ROOT TChain.
Auxiliary class to convert pair of indices to unique index and back.
int getIndex(const int first, const int second) const
Get index of pair of indices.
void sort(JComparator_t comparator)
Sort address pairs.
size_t getNumberOfPairs() const
Get number of pairs.
Template L2 builder.
Definition JBuildL2.hh:49
Data structure for hit time and DOM identifier.
JElement(const int __id, const double __t)
JCombinatorics::pair_type pair_type
int getNumberOfFloors(const JDetector &detector)
Get number of floors.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
double getDistance(const JFirst_t &first, const JSecond_t &second)
Get distance between objects.
const double getSpeedOfLight()
Get speed of light.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
std::set< JROOTClassSelector > getROOTClassSelection(const bool option=false)
Get ROOT class selection.
Auxiliary classes and methods for multi-dimensional interpolations and histograms.
double getFrameTime()
Get frame time duration.
Definition JDAQClock.hh:162
bool operator<(const JElement &first, const JElement &second)
Sort JElement in time.
bool comparepair(const pair_type &A, const pair_type &B)
Detector file.
Definition JHead.hh:227
Auxiliary class to select ROOT class based on class name.
Data structure for a pair of indices.
Auxiliary class to convert value to element.
Definition JMergeSort.cc:66
Data structure for L2 parameters.
Auxiliary data structure for histogram management.
JHistogram(TH2D *__h2s, TH1D *__h1l)
Constructor.
JHistogram()
Default constructor.