Jpp  17.2.1-pre0
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  cout.tie(&cerr);
256 
257 
259 
260  try {
261  load(detectorFile, detector);
262  }
263  catch(const JException& error) {
264  FATAL(error);
265  }
266 
267  const JModuleRouter moduleRouter(detector);
268  const JPMTRouter pmtRouter (detector);
269 
270 
271  detector -= get<JPosition3D>(getHeader(inputFile));
272 
273 
274  TFile out(outputFile.c_str(), "recreate");
275 
276  TProfile* he;
277  TProfile* hp;
278 
279  if (histogram == 1) {
280 
282 
283  double x = -0.5;
284 
285  for ( ; x < 30.0; x += 1.0)
286  X.push_back(x);
287 
288  for ( ; x < 50.0; x += 2.0)
289  X.push_back(x);
290 
291  for ( ; x < 100.0; x += 5.0)
292  X.push_back(x);
293 
294  for ( ; x < 200.0; x += 10.0)
295  X.push_back(x);
296 
297  he = new TProfile("he", NULL, X.size() - 1, X.data());
298  hp = new TProfile("hp", NULL, X.size() - 1, X.data());
299 
300  } else {
301 
302  he = new TProfile("he", NULL, 28, 0.0, 7.0);
303  hp = new TProfile("hp", NULL, 28, 0.0, 7.0);
304  }
305 
306  TH1D ht1("ht1", NULL, 550, -50.0, +500.0);
307  TH1D ht2("ht2", NULL, 550, -50.0, +500.0);
308  TH1D hd1("hd1", NULL, 100, -1.0, +1.0);
309  TH1D hd2("hd2", NULL, 100, -1.0, +1.0);
310  TH1D hx1("hx1", NULL, 100, 0.0, +250.0);
311  TH1D hx2("hx2", NULL, 100, 0.0, +250.0);
312  TH1D hw1("hw1", NULL, 100, -0.5, +99.5);
313  TH1D hw2("hw2", NULL, 100, -0.5, +99.5);
314 
315 
316 
317  typedef map<int, JBuffer> JMap_t;
318  typedef vector<JHitL1> JDataL1_t;
319 
320  JMap_t zmap;
321  const JBuildL2<JHitL1> buildL2(JL2Parameters(2, Tmax_ns, ctMin));
322  const JMatch3B<JHitL1> match3B(roadWidth_m, Tmax_ns);
323  const JMatch3D<JHitL1> match3D(Tmax_ns);
324  const JMatch1D<JHitL1> match1D(roadWidth_m, Tmax_ns);
325  const JHorizon horizon(0.2, 2);
326 
327  const JVersor3D gui(getSinThetaC(), 0.0, getCosThetaC()); // photon emission direction
328 
329 
330  while (inputFile.hasNext()) {
331 
332  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
333 
334  JTriggeredFileScanner<>::multi_pointer_type ps = inputFile.next();
335 
336  const JDAQEvent* tev = ps;
337  const Evt* event = ps;
338 
339  if (has_neutrino(*event)) {
340 
341  const Trk& neutrino = get_neutrino(*event);
342  const double E = neutrino.E; // GeV
343 
344 
345  vector<Trk>::const_iterator muon = event->mc_trks.end();
346 
347  {
348  double Emax = 0.0;
349 
350  for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
351 
352  // (anti) muon
353 
354  if (is_muon(*track) && track->E > Emax) {
355  muon = track;
356  Emax = track->E;
357  }
358  }
359  }
360 
361  if (muon != event->mc_trks.end()) {
362 
363  JPosition3D pos = getPosition (*muon);
364  double t0 = muon->t;
365 
366  const JRotation3D R(getDirection(*muon));
367 
368  pos.rotate(R);
369 
370  const time_converter converter(*event, *tev);
371 
372  // Monte Carlo hits
373 
374  zmap.clear();
375 
376  const vector<Hit>& hitlist = event->mc_hits;
377 
378  for (vector<Hit>::const_iterator i = hitlist.begin(); i != hitlist.end(); ++i) {
379  if (!is_noise(*i)) {
380  zmap[pmtRouter.getParentModuleID(i->pmt_id)].push_back(*i);
381  }
382  }
383 
384 
385  // L1 Monte Carlo true data
386 
387  for (JMap_t::iterator i = zmap.begin(); i != zmap.end(); ) {
388 
389  i->second.compile(Tmax_ns, 2);
390 
391  if (i->second.empty())
392  zmap.erase(i++);
393  else
394  ++i;
395  }
396 
397  const int L1mc = zmap.size();
398 
399 
400  // L1 triggered data
401 
402  JDataL1_t dataL1;
403  JDataL1_t zbuf;
404 
405  JDAQTimeslice timeslice(*tev, true);
406 
407  for (JDAQTimeslice::const_iterator i = timeslice.begin(); i != timeslice.end(); ++i) {
408 
409  zbuf.clear();
410 
411  buildL2(*i, moduleRouter.getModule(i->getModuleID()), back_inserter(zbuf));
412 
413  if (!zbuf.empty()) {
414 
415  sort(zbuf.begin(), zbuf.end(), timeSorter<JHitL1>);
416 
417  dataL1.push_back(*zbuf.begin());
418  }
419  }
420 
421 
422  // cluster
423 
424  JDataL1_t::iterator __end = dataL1.begin();
425 
426  switch (cluster) {
427 
428  case 'A':
429  //__end = dataL1.end();
430  __end = clusterize(dataL1.begin(), dataL1.end(), make_match(matchAll));
431  break;
432 
433  case 'B':
434  __end = clusterize(dataL1.begin(), dataL1.end(), weightSorter<JHitL1>, match3D);
435  break;
436 
437  case 'C':
438  __end = clusterize(dataL1.begin(), dataL1.end(), match3B);
439  break;
440 
441  case 'D':
442  __end = clusterizeWeight(dataL1.begin(), dataL1.end(), match3B);
443  break;
444 
445  case 'E':
446 
447  for (JDataL1_t::iterator i = dataL1.begin(); i != dataL1.end(); ++i)
448  i->rotate(R);
449 
450  __end = clusterizeWeight(dataL1.begin(), dataL1.end(), match1D);
451 
452  for (JDataL1_t::iterator i = dataL1.begin(); i != dataL1.end(); ++i)
453  i->rotate_back(R);
454 
455  break;
456 
457  case 'F':
458 
459  __end = clusterizeWeight(dataL1.begin(), dataL1.end(), match3B);
460 
461  for (JDataL1_t::iterator i = dataL1.begin(); i != dataL1.end(); ++i)
462  i->rotate(R);
463 
464  __end = partition(dataL1.begin(), __end, horizon);
465  __end = clusterizeWeight(dataL1.begin(), __end, match1D);
466 
467  for (JDataL1_t::iterator i = dataL1.begin(); i != dataL1.end(); ++i)
468  i->rotate_back(R);
469 
470  break;
471 
472  default:
473  __end = dataL1.end();
474  break;
475  }
476 
477 
478  // L1 monitor
479 
480  JDataL1_t::iterator __q = __end;
481 
482  set<int> L1;
483  set<int> L1ok;
484 
485  for (vector<JHitL1>::iterator hit = dataL1.begin(); hit != __q; ) {
486 
487  L1.insert(hit->getModuleID());
488 
489  const double t1 = converter.getTime(hit->getT());
490 
491  if (zmap[hit->getModuleID()].has(JTimeRange(t1 - 0.5*Tmax_ns, t1 + 0.5*Tmax_ns))) {
492 
493  L1ok.insert(hit->getModuleID());
494 
495  ++hit;
496 
497  } else {
498 
499  iter_swap(hit, --__q);
500  }
501  }
502 
503 
504  if (L1mc != 0 && !L1.empty()) {
505 
506  Double_t x = numeric_limits<Double_t>::max();
507 
508  if (histogram == 1)
509  x = L1mc;
510  else
511  x = log10(E);
512 
513  he->Fill(x, (Double_t) L1ok.size() / (Double_t) L1mc);
514  hp->Fill(x, (Double_t) L1ok.size() / (Double_t) L1.size());
515 
516 
517  // transform hits to muon system
518 
519  for (JDataL1_t::iterator hit = dataL1.begin(); hit != __end; ++hit) {
520  hit->transform(R, pos);
521  }
522 
523 
524  for (JDataL1_t::iterator hit = dataL1.begin(); hit != __end; ++hit) {
525 
526  const double t1 = t0 + getInverseSpeedOfLight() * (hit->getZ() + hit->getX() * getTanThetaC());
527 
528  double dot = +1.0;
529 
530  for (JHitL1::const_iterator i = hit->begin(); i != hit->end(); ++i) {
531  if (i->getDot(gui) < dot)
532  dot = i->getDot(gui);
533  }
534 
535  const double w = (hit->rbegin()->getT() - hit->begin()->getT()); // / hit->size();
536 
537  if (distance(hit, __q) > 0) {
538  ht1.Fill(converter.getTime(hit->getT()) - t1);
539  hd1.Fill(dot);
540  hx1.Fill(hit->getX());
541  hw1.Fill(w);
542  } else {
543  ht2.Fill(converter.getTime(hit->getT()) - t1);
544  hd2.Fill(dot);
545  hx2.Fill(hit->getX());
546  hw2.Fill(w);
547  }
548  }
549  }
550  }
551  }
552  }
553  STATUS(endl);
554 
555  out.Write();
556  out.Close();
557 }
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:1517
General exception.
Definition: JException.hh:23
data_type w[N+1][M+1]
Definition: JPolint.hh:778
1D match criterion.
Definition: JMatch1D.hh:31
T getLowerLimit() const
Get lower limit.
Definition: JRange.hh:202
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 E
Definition: JMuonPostfit.sh:35
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.
JTOOLS::JRange< double > JTimeRange
Type definition for time range (unit [ns]).
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
then JShowerPostfit f $INPUT_FILE o $OUTPUT_FILE N
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:117
Auxiliary class to convert DAQ hit time to/from Monte Carlo hit time.
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:697
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
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:1993
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.
then usage $script[distance] fi case set_variable R
Definition: JDrawLED.sh:43
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.
static struct JTRIGGER::@81 clusterizeWeight
Anonymous struct for weighed clustering of hits.
double getZ() const
Get z position.
Definition: JHitL1.hh:157
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:703
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
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
static struct JTRIGGER::@79 clusterize
Anonymous structure for clustering of hits.
int getParentModuleID(const JObjectID &id) const
Get parent module identifier.
Definition: JPMTRouter.hh:176
3D match criterion.
Definition: JMatch3D.hh:29
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
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