Jpp test-rotations-new
the software that should make you happy
Loading...
Searching...
No Matches
JFilter.cc
Go to the documentation of this file.
1
2#include <string>
3#include <iostream>
4#include <iomanip>
5#include <vector>
6#include <set>
7#include <map>
8#include <algorithm>
9#include <iterator>
10
11#include "TROOT.h"
12#include "TFile.h"
13#include "TH1D.h"
14#include "TProfile.h"
15
20
21#include "JAAnet/JHead.hh"
24
25#include "JDAQ/JDAQEventIO.hh"
28
29#include "JSirene/JPulse.hh"
30
32
38
39#include "JTrigger/JHit.hh"
40#include "JTrigger/JFrame.hh"
42#include "JTrigger/JHitL0.hh"
43#include "JTrigger/JHitL1.hh"
44#include "JTrigger/JBuildL1.hh"
45#include "JTrigger/JBuildL2.hh"
47#include "JTrigger/JMatch3B.hh"
48#include "JTrigger/JMatch3D.hh"
49#include "JTrigger/JMatch1D.hh"
50#include "JTrigger/JBind2nd.hh"
51
54#include "JSupport/JSupport.hh"
55
56#include "Jeep/JParser.hh"
57#include "Jeep/JMessage.hh"
58
59
60namespace {
61
62 using namespace JPP;
63
64
65 /**
66 * Data structure for Monte Carlo hits.
67 */
68 class JBuffer :
69 public std::vector<JPulse>
70 {
71 public:
72 /**
73 * Default constructor.
74 */
75 JBuffer() :
76 std::vector<JPulse>()
77 {}
78
79
80 /**
81 * Compile data to L1 hits.
82 *
83 * An L1 hit consists of two (or more) hits from different PMTs which are coincident in time.
84 * The start(stop) time of the L1 hit is defined by the earliest(latest) time of the hits
85 * that constitute the L1 hit.
86 * The PMT identifier of the earliest hit is maintained.
87 * The compilation may reduce the number of hits.
88 *
89 * \param Tmax maximal time difference between hits [ns]
90 * \param L1 minimum number of L0 hits to trigger L1
91 */
92 void compile(const double Tmax,
93 const unsigned int L1 = 2)
94 {
95 std::sort(this->begin(), this->end());
96
97 iterator out = this->begin();
98
99 for (const_iterator i = this->begin(); i != this->end(); ) {
100
101 if (L1 == 1) {
102
103 *out = JPulse(*i, *i);
104 ++out;
105
106 ++i;
107
108 } else if (L1 > 1) {
109
110 std::set<int> buffer;
111
112 buffer.insert(i->getID());
113
114 const_iterator j = i;
115
116 for (double t1 = i->getLowerLimit(); ++j != this->end() && j->getLowerLimit() - t1 < Tmax; t1 = j->getLowerLimit()) {
117 buffer.insert(j->getID());
118 }
119
120 if (buffer.size() >= L1) {
121 *out = JPulse(*i, *((--j)++));
122 ++out;
123 }
124
125 i = j;
126 }
127 }
128
129 this->erase(out, end());
130 }
131
132
133 /**
134 * Check if a hit with given time range is present in buffer.
135 *
136 * \param t1 minimal time [ns]
137 * \param t2 maximal time [ns]
138 * \return true if hit present; else false
139 */
140 inline bool has(const JTimeRange& timerange) const
141 {
142 const_iterator i = std::lower_bound(this->begin(), this->end(), timerange.getLowerLimit());
143
144 return i != this->end() && i->overlap(timerange);
145 }
146 };
147
148
149 using namespace JPP;
150
151 /**
152 * Match all L1 hits.
153 *
154 * \param first first hit
155 * \param second second hit
156 * \return true
157 */
158 inline bool matchAll(const JHitL1& first, const JHitL1& second)
159 {
160 return true;
161 }
162
163
164 /**
165 * Auxiliary class to select PMTs looking below a given horizon.
166 */
167 struct JHorizon {
168 /**
169 * Constructor.
170 *
171 * \param ct cosine elevation angle of horizon
172 * \param N mininum number of PMTs below horizon
173 */
174 JHorizon(const double ct,
175 const int N)
176 {
177 this->ct = ct;
178 this->N = N;
179 }
180
181
182 /**
183 * Select L1 hit.
184 *
185 * \param hit L1 hit
186 * \return true if accepted; false if rejected
187 */
188 bool operator()(const JHitL1& hit) const
189 {
190 int n = 0;
191
192 for (JHitL1::const_iterator i = hit.begin(); i != hit.end(); ++i) {
193 if (i->getDZ() < ct)
194 ++n;
195 }
196
197 return n >= N;
198 }
199
200
201 double ct;
202 int N;
203 };
204}
205
206
207/**
208 * \file
209 * Example program to test performance of various hit filters based on JTRIGGER::clusterize method.
210 * \author mdejong
211 */
212int main(int argc, char **argv)
213{
214 using namespace std;
215 using namespace JPP;
216 using namespace KM3NETDAQ;
217
218 JTriggeredFileScanner<> inputFile;
219 JLimit_t& numberOfEvents = inputFile.getLimit();
220 string outputFile;
221 string detectorFile;
222 double Tmax_ns;
223 double roadWidth_m;
224 double ctMin;
225 double thetaRad;
226 char cluster;
227 int histogram;
228 int debug;
229
230 try {
231
232 JParser<> zap("Example program to test performance of various hit filters.");
233
234 zap['f'] = make_field(inputFile);
235 zap['o'] = make_field(outputFile) = "filter.root";
236 zap['a'] = make_field(detectorFile);
237 zap['n'] = make_field(numberOfEvents) = JLimit::max();
238 zap['T'] = make_field(Tmax_ns) = 20.0; // [ns]
239 zap['R'] = make_field(roadWidth_m) = 150.0; // [m]
240 zap['C'] = make_field(ctMin) = 0.0; //
241 zap['Q'] = make_field(thetaRad) = 0.2; //
242 zap['c'] = make_field(cluster) = 'A', 'B', 'C', 'D', 'E', 'F';
243 zap['H'] = make_field(histogram) = 1, 2;
244 zap['d'] = make_field(debug) = 1;
245
246 zap(argc, argv);
247 }
248 catch(const exception &error) {
249 FATAL(error.what() << endl);
250 }
251
252
253 using namespace KM3NETDAQ;
254
255
257
258 try {
259 load(detectorFile, detector);
260 }
261 catch(const JException& error) {
262 FATAL(error);
263 }
264
265 const JModuleRouter moduleRouter(detector);
266 const JPMTRouter pmtRouter (detector);
267
268 const Vec offset = getOffset(getHeader(inputFile));
269
270 detector -= getPosition(offset);
271
272
273 TFile out(outputFile.c_str(), "recreate");
274
275 TProfile* he;
276 TProfile* hp;
277
278 if (histogram == 1) {
279
281
282 double x = -0.5;
283
284 for ( ; x < 30.0; x += 1.0)
285 X.push_back(x);
286
287 for ( ; x < 50.0; x += 2.0)
288 X.push_back(x);
289
290 for ( ; x < 100.0; x += 5.0)
291 X.push_back(x);
292
293 for ( ; x < 200.0; x += 10.0)
294 X.push_back(x);
295
296 he = new TProfile("he", NULL, X.size() - 1, X.data());
297 hp = new TProfile("hp", NULL, X.size() - 1, X.data());
298
299 } else {
300
301 he = new TProfile("he", NULL, 28, 0.0, 7.0);
302 hp = new TProfile("hp", NULL, 28, 0.0, 7.0);
303 }
304
305 TH1D ht1("ht1", NULL, 550, -50.0, +500.0);
306 TH1D ht2("ht2", NULL, 550, -50.0, +500.0);
307 TH1D hd1("hd1", NULL, 100, -1.0, +1.0);
308 TH1D hd2("hd2", NULL, 100, -1.0, +1.0);
309 TH1D hx1("hx1", NULL, 100, 0.0, +250.0);
310 TH1D hx2("hx2", NULL, 100, 0.0, +250.0);
311 TH1D hw1("hw1", NULL, 100, -0.5, +99.5);
312 TH1D hw2("hw2", NULL, 100, -0.5, +99.5);
313
314
315 typedef map<int, JBuffer> JMap_t;
316 typedef vector<JHitL1> JDataL1_t;
317
318 JMap_t zmap;
319 const JBuildL2<JHitL1> buildL2(JL2Parameters(2, Tmax_ns, ctMin));
320 const JMatch3B<JHitL1> match3B(roadWidth_m, Tmax_ns);
321 const JMatch3D<JHitL1> match3D(Tmax_ns);
322 const JMatch1D<JHitL1> match1D(roadWidth_m, Tmax_ns);
323 const JHorizon horizon(0.2, 2);
324
325 const JVersor3D gui(getSinThetaC(), 0.0, getCosThetaC()); // photon emission direction
326
327
328 while (inputFile.hasNext()) {
329
330 STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
331
333
334 const JDAQEvent* tev = ps;
335 const Evt* event = ps;
336
337 if (has_neutrino(*event)) {
338
339 const Trk& neutrino = get_neutrino(*event);
340 const double E = neutrino.E; // GeV
341
342
343 vector<Trk>::const_iterator muon = event->mc_trks.end();
344
345 {
346 double Emax = 0.0;
347
348 for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
349
350 // (anti) muon
351
352 if (is_muon(*track) && track->E > Emax) {
353 muon = track;
354 Emax = track->E;
355 }
356 }
357 }
358
359 if (muon != event->mc_trks.end()) {
360
361 JPosition3D pos = getPosition (*muon);
362 double t0 = muon->t;
363
364 const JRotation3D R(getDirection(*muon));
365
366 pos.rotate(R);
367
368 const time_converter converter(*event, *tev);
369
370 // Monte Carlo hits
371
372 zmap.clear();
373
374 const vector<Hit>& hitlist = event->mc_hits;
375
376 for (vector<Hit>::const_iterator i = hitlist.begin(); i != hitlist.end(); ++i) {
377 if (!is_noise(*i)) {
378 zmap[pmtRouter.getParentModuleID(i->pmt_id)].push_back(*i);
379 }
380 }
381
382
383 // L1 Monte Carlo true data
384
385 for (JMap_t::iterator i = zmap.begin(); i != zmap.end(); ) {
386
387 i->second.compile(Tmax_ns, 2);
388
389 if (i->second.empty())
390 zmap.erase(i++);
391 else
392 ++i;
393 }
394
395 const int L1mc = zmap.size();
396
397
398 // L1 triggered data
399
400 JDataL1_t dataL1;
401 JDataL1_t zbuf;
402
403 JDAQTimeslice timeslice(*tev, true);
404
405 for (JDAQTimeslice::const_iterator i = timeslice.begin(); i != timeslice.end(); ++i) {
406
407 zbuf.clear();
408
409 buildL2(*i, moduleRouter.getModule(i->getModuleID()), back_inserter(zbuf));
410
411 if (!zbuf.empty()) {
412
413 sort(zbuf.begin(), zbuf.end(), timeSorter<JHitL1>);
414
415 dataL1.push_back(*zbuf.begin());
416 }
417 }
418
419
420 // cluster
421
422 JDataL1_t::iterator __end = dataL1.begin();
423
424 switch (cluster) {
425
426 case 'A':
427 //__end = dataL1.end();
428 __end = clusterize(dataL1.begin(), dataL1.end(), make_match(matchAll));
429 break;
430
431 case 'B':
432 __end = clusterize(dataL1.begin(), dataL1.end(), weightSorter<JHitL1>, match3D);
433 break;
434
435 case 'C':
436 __end = clusterize(dataL1.begin(), dataL1.end(), match3B);
437 break;
438
439 case 'D':
440 __end = clusterizeWeight(dataL1.begin(), dataL1.end(), match3B);
441 break;
442
443 case 'E':
444
445 for (JDataL1_t::iterator i = dataL1.begin(); i != dataL1.end(); ++i)
446 i->rotate(R);
447
448 __end = clusterizeWeight(dataL1.begin(), dataL1.end(), match1D);
449
450 for (JDataL1_t::iterator i = dataL1.begin(); i != dataL1.end(); ++i)
451 i->rotate_back(R);
452
453 break;
454
455 case 'F':
456
457 __end = clusterizeWeight(dataL1.begin(), dataL1.end(), match3B);
458
459 for (JDataL1_t::iterator i = dataL1.begin(); i != dataL1.end(); ++i)
460 i->rotate(R);
461
462 __end = partition(dataL1.begin(), __end, horizon);
463 __end = clusterizeWeight(dataL1.begin(), __end, match1D);
464
465 for (JDataL1_t::iterator i = dataL1.begin(); i != dataL1.end(); ++i)
466 i->rotate_back(R);
467
468 break;
469
470 default:
471 __end = dataL1.end();
472 break;
473 }
474
475
476 // L1 monitor
477
478 JDataL1_t::iterator __q = __end;
479
480 set<int> L1;
481 set<int> L1ok;
482
483 for (vector<JHitL1>::iterator hit = dataL1.begin(); hit != __q; ) {
484
485 L1.insert(hit->getModuleID());
486
487 const double t1 = converter.getTime(hit->getT());
488
489 if (zmap[hit->getModuleID()].has(JTimeRange(t1 - 0.5*Tmax_ns, t1 + 0.5*Tmax_ns))) {
490
491 L1ok.insert(hit->getModuleID());
492
493 ++hit;
494
495 } else {
496
497 iter_swap(hit, --__q);
498 }
499 }
500
501
502 if (L1mc != 0 && !L1.empty()) {
503
504 Double_t x = numeric_limits<Double_t>::max();
505
506 if (histogram == 1)
507 x = L1mc;
508 else
509 x = log10(E);
510
511 he->Fill(x, (Double_t) L1ok.size() / (Double_t) L1mc);
512 hp->Fill(x, (Double_t) L1ok.size() / (Double_t) L1.size());
513
514
515 // transform hits to muon system
516
517 for (JDataL1_t::iterator hit = dataL1.begin(); hit != __end; ++hit) {
518 hit->transform(R, pos);
519 }
520
521
522 for (JDataL1_t::iterator hit = dataL1.begin(); hit != __end; ++hit) {
523
524 const double t1 = t0 + getInverseSpeedOfLight() * (hit->getZ() + hit->getX() * getTanThetaC());
525
526 double dot = +1.0;
527
528 for (JHitL1::const_iterator i = hit->begin(); i != hit->end(); ++i) {
529 if (i->getDot(gui) < dot)
530 dot = i->getDot(gui);
531 }
532
533 const double w = (hit->rbegin()->getT() - hit->begin()->getT()); // / hit->size();
534
535 if (distance(hit, __q) > 0) {
536 ht1.Fill(converter.getTime(hit->getT()) - t1);
537 hd1.Fill(dot);
538 hx1.Fill(hit->getX());
539 hw1.Fill(w);
540 } else {
541 ht2.Fill(converter.getTime(hit->getT()) - t1);
542 hd2.Fill(dot);
543 hx2.Fill(hit->getX());
544 hw2.Fill(w);
545 }
546 }
547 }
548 }
549 }
550 }
551 STATUS(endl);
552
553 out.Write();
554 out.Close();
555}
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
Algorithms for hit clustering and sorting.
string outputFile
Data structure for detector geometry and calibration.
int main(int argc, char **argv)
Definition JFilter.cc:212
Basic data structure for L0 hit.
Basic data structure for L1 hit.
Match operator for Cherenkov light from muon with given direction.
Match operator for Cherenkov light from muon in any direction.
Match operator for Cherenkov light from muon in any direction.
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.
Direct access to PMT in detector data structure.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
Physics constants.
Time-over-threshold (ToT) pulse from a PMT.
ROOT TTree parameter settings of various packages.
Basic data structure for time and time over threshold information of hit.
Synchronously read DAQ events and Monte Carlo events (and optionally other events).
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
3D match criterion for acoustic signals.
Definition JBillabong.cc:96
Detector data structure.
Definition JDetector.hh:96
Router for direct addressing of module data in detector data structure.
const JModule & getModule(const JObjectID &id) const
Get module parameters.
Router for direct addressing of PMT data in detector data structure.
Definition JPMTRouter.hh:37
int getParentModuleID(const JObjectID &id) const
Get parent module identifier.
Data structure for position in three dimensions.
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Data structure for normalised vector in three dimensions.
Definition JVersor3D.hh:28
General exception.
Definition JException.hh:24
Utility class to parse command line options.
Definition JParser.hh:1698
Auxiliary class for a time-over-threshold pulse from a PMT.
Definition JPulse.hh:36
counter_type getCounter() const
Get counter.
T getLowerLimit() const
Get lower limit.
Definition JRange.hh:202
Template L2 builder.
Definition JBuildL2.hh:49
Data structure for L1 hit.
Definition JHitL1.hh:37
1D match criterion.
Definition JMatch1D.hh:33
3D match criterion with road width.
Definition JMatch3B.hh:36
Auxiliary class to convert DAQ hit time to/from Monte Carlo hit time.
double getTime() const
Get DAQ/trigger time minus Monte Carlo time.
JDirection3D getDirection(const Vec &dir)
Get direction.
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
JPosition3D getPosition(const Vec &pos)
Get position.
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
bool is_noise(const Hit &hit)
Verify hit origin.
Vec getOffset(const JHead &header)
Get offset.
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
double getCosThetaC()
Get average cosine of Cherenkov angle of water corresponding to group velocity.
const double getInverseSpeedOfLight()
Get inverse speed of light.
double getSinThetaC()
Get average sine of Cherenkov angle of water corresponding to group velocity.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
const int n
Definition JPolint.hh:791
int j
Definition JPolint.hh:801
JHitIterator_t clusterize(JHitIterator_t __begin, JHitIterator_t __end, const JMatch_t &match, const int Nmin=1)
Partition data according given binary match operator.
Definition JAlgorithm.hh:38
JHitIterator_t clusterizeWeight(JHitIterator_t __begin, JHitIterator_t __end, const JMatch_t &match)
Partition data according given binary match operator.
bool weightSorter(const JHit_t &first, const JHit_t &second)
Compare two hits by weight.
bool timeSorter(const JHit_t &first, const JHit_t &second)
Compare two hits by weight.
JMatchHelper< JHit_t > make_match(bool(*match)(const JHit_t &, const JHit_t &))
Auxiliary method to make JMatch object based on pointer to match function.
Definition JMatch.hh:115
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
Definition gui.cpp:136
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition Evt.hh:21
Detector file.
Definition JHead.hh:227
General purpose class for multiple pointers.
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
Auxiliary class to synchronously read DAQ events and Monte Carlo events (and optionally other events)...
virtual bool hasNext() override
Check availability of next element.
virtual const multi_pointer_type & next() override
Get next element.
Data structure for L2 parameters.
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition Trk.hh:15
double E
Energy [GeV] (either MC truth or reconstructed)
Definition Trk.hh:20
The Vec class is a straightforward 3-d vector, which also works in pyroot.
Definition Vec.hh:13
Auxiliary include file for time conversion between DAQ/trigger hit and Monte Carlo hit.