Jpp  18.6.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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"
22 #include "JAAnet/JHeadToolkit.hh"
23 #include "JAAnet/JAAnetToolkit.hh"
24 
25 #include "JDAQ/JDAQEventIO.hh"
26 #include "JDAQ/JDAQTimesliceIO.hh"
28 
29 #include "JSirene/JPulse.hh"
30 
31 #include "JPhysics/JConstants.hh"
32 
33 #include "JDetector/JDetector.hh"
36 #include "JDetector/JPMTRouter.hh"
37 #include "JDetector/JTimeRange.hh"
38 
39 #include "JTrigger/JHit.hh"
40 #include "JTrigger/JFrame.hh"
41 #include "JTrigger/JTimeslice.hh"
42 #include "JTrigger/JHitL0.hh"
43 #include "JTrigger/JHitL1.hh"
44 #include "JTrigger/JBuildL1.hh"
45 #include "JTrigger/JBuildL2.hh"
46 #include "JTrigger/JAlgorithm.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 
60 namespace {
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  */
212 int 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 
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 }
Match operator for Cherenkov light from muon in any direction.
Router for direct addressing of PMT data in detector data structure.
Definition: JPMTRouter.hh:35
Utility class to parse command line options.
Definition: JParser.hh:1711
General exception.
Definition: JException.hh:24
data_type w[N+1][M+1]
Definition: JPolint.hh:867
1D match criterion.
Definition: JMatch1D.hh:31
T getLowerLimit() const
Get lower limit.
Definition: JRange.hh:202
static struct JTRIGGER::clusterizeWeight clusterizeWeight
int main(int argc, char *argv[])
Definition: Main.cc:15
ROOT TTree parameter settings of various packages.
int getModuleID() const
Get module identifier.
Match operator for Cherenkov light from muon with given direction.
Data structure for L1 hit.
Definition: JHitL1.hh:34
const JModule & getModule(const JObjectID &id) const
Get module parameters.
then usage $script< input file >[option[primary[working directory]]] nWhere option can be E
Definition: JMuonPostfit.sh:40
Vec getOffset(const JHead &header)
Get offset.
Algorithms for hit clustering and sorting.
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.
#define STATUS(A)
Definition: JMessage.hh:63
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
Detector data structure.
Definition: JDetector.hh:89
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
Router for direct addressing of module data in detector data structure.
Rotation matrix.
Definition: JRotation3D.hh:111
Auxiliary class to synchronously read DAQ events and Monte Carlo events (and optionally other events)...
double getT(const unsigned int i) const
Get time of hit i.
Definition: JHitL1.hh:170
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
static counter_type max()
Get maximum counter value.
Definition: JLimit.hh:128
then warning Cannot perform comparison test for histogram
Auxiliary class to convert DAQ hit time to/from Monte Carlo hit time.
3D match criterion for acoustic signals.
Basic data structure for time and time over threshold information of hit.
string outputFile
double E
Energy [GeV] (either MC truth or reconstructed)
Definition: Trk.hh:20
Data structure for detector geometry and calibration.
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
bool is_noise(const Hit &hit)
Verify hit origin.
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
Basic data structure for L0 hit.
double getCosThetaC()
Get average cosine of Cherenkov angle of water corresponding to group velocity.
const int n
Definition: JPolint.hh:786
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
The Vec class is a straightforward 3-d vector, which also works in pyroot.
Definition: Vec.hh:12
Template L2 builder.
Definition: JBuildL2.hh:45
void transform(const JRotation3D &R, const JVector3D &pos)
Transform hit.
Definition: JHitL1.hh:335
Detector file.
Definition: JHead.hh:226
Match operator for Cherenkov light from muon in any direction.
JDirection3D getDirection(const Vec &dir)
Get direction.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
set_variable E_E log10(E_{fit}/E_{#mu})"
JPosition3D getPosition(const Vec &pos)
Get position.
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
Data time slice.
Physics constants.
Direct access to PMT in detector data structure.
General purpose messaging.
Auxiliary include file for time conversion between DAQ/trigger hit and Monte Carlo hit...
#define FATAL(A)
Definition: JMessage.hh:67
Direct access to module in detector data structure.
Time-over-threshold (ToT) pulse from a PMT.
Data structure for L2 parameters.
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
double getZ() const
Get z position.
Definition: JHitL1.hh:157
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Definition: JMuonPostfit.sh:40
Auxiliary class for a time-over-threshold pulse from a PMT.
Definition: JPulse.hh:33
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Utility class to parse command line options.
const double getInverseSpeedOfLight()
Get inverse speed of light.
no fit printf nominal n $STRING awk v X
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
int j
Definition: JPolint.hh:792
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:84
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
do set_variable DETECTOR_TXT $WORKDIR detector
Data structure for normalised vector in three dimensions.
Definition: JVersor3D.hh:26
General purpose class for multiple pointers.
double getTime() const
Get DAQ/trigger time minus Monte Carlo time.
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition: Trk.hh:14
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:186
double getX() const
Get x position.
Definition: JHitL1.hh:133
int getParentModuleID(const JObjectID &id) const
Get parent module identifier.
Definition: JPMTRouter.hh:176
static struct JTRIGGER::clusterize clusterize
Basic data structure for L1 hit.
int debug
debug level
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:20
JTOOLS::JRange< double > JTimeRange
Type definition for time range (unit [s]).
double getSinThetaC()
Get average sine of Cherenkov angle of water corresponding to group velocity.
3D match criterion with road width.
Definition: JMatch3B.hh:34
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62