Jpp 20.0.0-27-g39925593c-D
the software that should make you happy
Loading...
Searching...
No Matches
Functions
JFilter.cc File Reference

Example program to test performance of various hit filters based on JTRIGGER::clusterize method. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include <map>
#include <algorithm>
#include <iterator>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TProfile.h"
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Hit.hh"
#include "km3net-dataformat/tools/time_converter.hh"
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JDAQ/JDAQEventIO.hh"
#include "JDAQ/JDAQTimesliceIO.hh"
#include "km3net-dataformat/online/JDAQClock.hh"
#include "JSirene/JPulse.hh"
#include "JPhysics/JConstants.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JPMTRouter.hh"
#include "JDetector/JTimeRange.hh"
#include "JTrigger/JHit.hh"
#include "JTrigger/JFrame.hh"
#include "JTrigger/JTimeslice.hh"
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JHitL1.hh"
#include "JTrigger/JBuildL1.hh"
#include "JTrigger/JBuildL2.hh"
#include "JTrigger/JAlgorithm.hh"
#include "JTrigger/JMatch3B.hh"
#include "JTrigger/JMatch3D.hh"
#include "JTrigger/JMatch1D.hh"
#include "JTrigger/JBind2nd.hh"
#include "JSupport/JTriggeredFileScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Example program to test performance of various hit filters based on JTRIGGER::clusterize method.

Author
mdejong

Definition in file JFilter.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 212 of file JFilter.cc.

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
256 JDetector detector;
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
280 vector<double> X;
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
332 JTriggeredFileScanner<>::multi_pointer_type ps = inputFile.next();
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}
outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37
DEBUG
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
STATUS
#define STATUS(A)
Definition JMessage.hh:63
FATAL
#define FATAL(A)
Definition JMessage.hh:67
debug
int debug
debug level
Definition JSirene.cc:72
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
distance
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Definition build/antares-dataformat/PhysicsEvent.hh:434
JDETECTOR::JDetector
Detector data structure.
Definition JDetector.hh:96
JDETECTOR::JModuleRouter
Router for direct addressing of module data in detector data structure.
Definition JModuleRouter.hh:38
JDETECTOR::JPMTRouter
Router for direct addressing of PMT data in detector data structure.
Definition JPMTRouter.hh:37
JGEOMETRY3D::JPosition3D
Data structure for position in three dimensions.
Definition JPosition3D.hh:38
JGEOMETRY3D::JPosition3D::rotate
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition JPosition3D.hh:186
JGEOMETRY3D::JVersor3D
Data structure for normalised vector in three dimensions.
Definition JVersor3D.hh:28
JLANG::JException
General exception.
Definition JException.hh:24
JOSCPROB::JOscProbInterpolator
Template definition of a multi-dimensional oscillation probability interpolation table.
Definition JOscProbInterpolator.hh:71
JTOOLS::JMultiFunction::insert
void insert(const JMultiFunction< __JFunction_t, __JMaplist_t, __JDistance_t > &input)
Insert multidimensional input.
Definition JMultiFunction.hh:121
time_converter
Auxiliary class to convert DAQ hit time to/from Monte Carlo hit time.
Definition time_converter.hh:20
JAANET::getDirection
JDirection3D getDirection(const Vec &dir)
Get direction.
Definition JAAnetToolkit.hh:218
JAANET::has_neutrino
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
Definition JAAnetToolkit.hh:520
JAANET::getPosition
JPosition3D getPosition(const Vec &pos)
Get position.
Definition JAAnetToolkit.hh:182
JAANET::is_muon
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
Definition JAAnetToolkit.hh:366
JAANET::is_noise
bool is_noise(const Hit &hit)
Verify hit origin.
Definition JAAnetToolkit.hh:124
JAANET::getOffset
Vec getOffset(const JHead &header)
Get offset.
Definition JHeadToolkit.hh:291
JAANET::get_neutrino
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
Definition JAAnetToolkit.hh:542
JACOUSTICS::JTimeRange
JTOOLS::JRange< double > JTimeRange
Type definition for time range (unit [s]).
Definition JAcoustics/JTimeRange.hh:19
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition JDetectorToolkit.hh:520
JPHYSICS::getCosThetaC
double getCosThetaC()
Get average cosine of Cherenkov angle of water corresponding to group velocity.
Definition JPhysics/JConstants.hh:167
JPHYSICS::getInverseSpeedOfLight
const double getInverseSpeedOfLight()
Get inverse speed of light.
Definition JPhysics/JConstants.hh:123
JPHYSICS::getTanThetaC
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
Definition JPhysics/JConstants.hh:156
JPHYSICS::getSinThetaC
double getSinThetaC()
Get average sine of Cherenkov angle of water corresponding to group velocity.
Definition JPhysics/JConstants.hh:178
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:45
JSUPPORT::getHeader
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
Definition JMonteCarloFileSupportkit.hh:547
JTRIGGER::clusterize
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
JTRIGGER::clusterizeWeight
JHitIterator_t clusterizeWeight(JHitIterator_t __begin, JHitIterator_t __end, const JMatch_t &match)
Partition data according given binary match operator.
Definition JAlgorithm.hh:257
JTRIGGER::make_match
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
KM3NETDAQ
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
gui
Definition gui.cpp:136
Evt
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition Evt.hh:21
JAANET::detector
Detector file.
Definition JHead.hh:227
JLANG::JMultiPointer
General purpose class for multiple pointers.
Definition JMultiPointer.hh:24
JSUPPORT::JLimit
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45
JSUPPORT::JLimit::max
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
JTRIGGER::JL2Parameters
Data structure for L2 parameters.
Definition JTrigger/JTriggerParameters.hh:34
Trk
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition Trk.hh:15
Trk::E
double E
Energy [GeV] (either MC truth or reconstructed)
Definition Trk.hh:20
Vec
The Vec class is a straightforward 3-d vector, which also works in pyroot.
Definition Vec.hh:13
Generated by doxygen 1.10.0