Jpp
 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 
19 #include "JAAnet/JHead.hh"
20 #include "JAAnet/JHeadToolkit.hh"
21 
22 #include "JDAQ/JDAQEventIO.hh"
23 #include "JDAQ/JDAQTimesliceIO.hh"
25 
26 #include "JSirene/JPulse.hh"
27 
28 #include "JTools/JConstants.hh"
29 
30 #include "JDetector/JDetector.hh"
33 #include "JDetector/JPMTRouter.hh"
34 #include "JDetector/JTimeRange.hh"
35 
36 #include "JTrigger/JHit.hh"
37 #include "JTrigger/JFrame.hh"
38 #include "JTrigger/JTimeslice.hh"
39 #include "JTrigger/JHitL0.hh"
40 #include "JTrigger/JHitL1.hh"
41 #include "JTrigger/JBuildL1.hh"
42 #include "JTrigger/JBuildL2.hh"
43 #include "JTrigger/JAlgorithm.hh"
44 #include "JTrigger/JMatch3B.hh"
45 #include "JTrigger/JMatch3D.hh"
46 #include "JTrigger/JMatch1D.hh"
47 #include "JTrigger/JBind2nd.hh"
49 
52 #include "JSupport/JSupport.hh"
53 
54 #include "Jeep/JParser.hh"
55 #include "Jeep/JMessage.hh"
56 
57 
58 namespace {
59 
60  using namespace JPP;
61 
62 
63  /**
64  * Data structure for Monte Carlo hits.
65  */
66  class JBuffer :
67  public std::vector<JPulse>
68  {
69  public:
70  /**
71  * Default constructor.
72  */
73  JBuffer() :
74  std::vector<JPulse>()
75  {}
76 
77 
78  /**
79  * Compile data to L1 hits.
80  *
81  * An L1 hit consists of two (or more) hits from different PMTs which are coincident in time.
82  * The start(stop) time of the L1 hit is defined by the earliest(latest) time of the hits
83  * that constitute the L1 hit.
84  * The PMT identifier of the earliest hit is maintained.
85  * The compilation may reduce the number of hits.
86  *
87  * \param Tmax maximal time difference between hits [ns]
88  * \param L1 minimum number of L0 hits to trigger L1
89  */
90  void compile(const double Tmax,
91  const unsigned int L1 = 2)
92  {
93  std::sort(this->begin(), this->end());
94 
95  iterator out = this->begin();
96 
97  for (const_iterator i = this->begin(); i != this->end(); ) {
98 
99  if (L1 == 1) {
100 
101  *out = JPulse(*i, *i);
102  ++out;
103 
104  ++i;
105 
106  } else if (L1 > 1) {
107 
108  std::set<int> buffer;
109 
110  buffer.insert(i->getID());
111 
112  const_iterator j = i;
113 
114  for (double t1 = i->getLowerLimit(); ++j != this->end() && j->getLowerLimit() - t1 < Tmax; t1 = j->getLowerLimit()) {
115  buffer.insert(j->getID());
116  }
117 
118  if (buffer.size() >= L1) {
119  *out = JPulse(*i, *((--j)++));
120  ++out;
121  }
122 
123  i = j;
124  }
125  }
126 
127  this->erase(out, end());
128  }
129 
130 
131  /**
132  * Check if a hit with given time range is present in buffer.
133  *
134  * \param t1 minimal time [ns]
135  * \param t2 maximal time [ns]
136  * \return true if hit present; else false
137  */
138  inline bool has(const JTimeRange& timerange) const
139  {
140  const_iterator i = std::lower_bound(this->begin(), this->end(), timerange.getLowerLimit());
141 
142  return i != this->end() && i->overlap(timerange);
143  }
144  };
145 
146 
147  using namespace JPP;
148 
149  /**
150  * Match all L1 hits.
151  *
152  * \param first first hit
153  * \param second second hit
154  * \return true
155  */
156  inline bool matchAll(const JHitL1& first, const JHitL1& second)
157  {
158  return true;
159  }
160 
161 
162  /**
163  * Auxiliary class to select PMTs looking below a given horizon.
164  */
165  struct JHorizon {
166  /**
167  * Constructor.
168  *
169  * \param ct cosine elevation angle of horizon
170  * \param N mininum number of PMTs below horizon
171  */
172  JHorizon(const double ct,
173  const int N)
174  {
175  this->ct = ct;
176  this->N = N;
177  }
178 
179 
180  /**
181  * Select L1 hit.
182  *
183  * \param hit L1 hit
184  * \return true if accepted; false if rejected
185  */
186  bool operator()(const JHitL1& hit) const
187  {
188  int n = 0;
189 
190  for (JHitL1::const_iterator i = hit.begin(); i != hit.end(); ++i) {
191  if (i->getDZ() < ct)
192  ++n;
193  }
194 
195  return n >= N;
196  }
197 
198 
199  double ct;
200  int N;
201  };
202 }
203 
204 
205 /**
206  * \file
207  * Example program to test performance of various hit filters based on JTRIGGER::clusterize method.
208  * \author mdejong
209  */
210 int main(int argc, char **argv)
211 {
212  using namespace std;
213  using namespace JPP;
214  using namespace KM3NETDAQ;
215 
216  JTriggeredFileScanner<> inputFile;
217  JLimit_t& numberOfEvents = inputFile.getLimit();
218  string outputFile;
219  string detectorFile;
220  double Tmax_ns;
221  double roadWidth_m;
222  double ctMin;
223  double thetaRad;
224  char cluster;
225  int histogram;
226  int debug;
227 
228  try {
229 
230  JParser<> zap("Example program to test performance of various hit filters.");
231 
232  zap['f'] = make_field(inputFile);
233  zap['o'] = make_field(outputFile) = "filter.root";
234  zap['a'] = make_field(detectorFile);
235  zap['n'] = make_field(numberOfEvents) = JLimit::max();
236  zap['T'] = make_field(Tmax_ns) = 20.0; // [ns]
237  zap['R'] = make_field(roadWidth_m) = 150.0; // [m]
238  zap['C'] = make_field(ctMin) = 0.0; //
239  zap['Q'] = make_field(thetaRad) = 0.2; //
240  zap['c'] = make_field(cluster) = 'A', 'B', 'C', 'D', 'E', 'F';
241  zap['H'] = make_field(histogram) = 1, 2;
242  zap['d'] = make_field(debug) = 1;
243 
244  zap(argc, argv);
245  }
246  catch(const exception &error) {
247  FATAL(error.what() << endl);
248  }
249 
250 
251  using namespace KM3NETDAQ;
252 
253  cout.tie(&cerr);
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 
269  detector -= get<JPosition3D>(getHeader(inputFile));
270 
271 
272  TFile out(outputFile.c_str(), "recreate");
273 
274  TProfile* he;
275  TProfile* hp;
276 
277  if (histogram == 1) {
278 
279  vector<double> X;
280 
281  double x = -0.5;
282 
283  for ( ; x < 30.0; x += 1.0)
284  X.push_back(x);
285 
286  for ( ; x < 50.0; x += 2.0)
287  X.push_back(x);
288 
289  for ( ; x < 100.0; x += 5.0)
290  X.push_back(x);
291 
292  for ( ; x < 200.0; x += 10.0)
293  X.push_back(x);
294 
295  he = new TProfile("he", NULL, X.size() - 1, X.data());
296  hp = new TProfile("hp", NULL, X.size() - 1, X.data());
297 
298  } else {
299 
300  he = new TProfile("he", NULL, 28, 0.0, 7.0);
301  hp = new TProfile("hp", NULL, 28, 0.0, 7.0);
302  }
303 
304  TH1D ht1("ht1", NULL, 550, -50.0, +500.0);
305  TH1D ht2("ht2", NULL, 550, -50.0, +500.0);
306  TH1D hd1("hd1", NULL, 100, -1.0, +1.0);
307  TH1D hd2("hd2", NULL, 100, -1.0, +1.0);
308  TH1D hx1("hx1", NULL, 100, 0.0, +250.0);
309  TH1D hx2("hx2", NULL, 100, 0.0, +250.0);
310  TH1D hw1("hw1", NULL, 100, -0.5, +99.5);
311  TH1D hw2("hw2", NULL, 100, -0.5, +99.5);
312 
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  JTimeConverter 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);
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:33
Utility class to parse command line options.
Definition: JParser.hh:1493
General exception.
Definition: JException.hh:23
data_type w[N+1][M+1]
Definition: JPolint.hh:708
1D match criterion.
Definition: JMatch1D.hh:31
T getLowerLimit() const
Get lower limit.
Definition: JRange.hh:215
ROOT TTree parameter settings.
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.
static struct JTRIGGER::@69 clusterize
Anonymous structure for clustering of hits.
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:80
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:181
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
double getTime() const
Get DAQ/trigger minus Monte Carlo hit time.
static counter_type max()
Get maximum counter value.
Definition: JLimit.hh:117
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:18
Data structure for detector geometry and calibration.
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
JRange< double > JTimeRange
Type definition for time range.
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.
static struct JTRIGGER::@71 clusterizeWeight
Anonymous struct for weighed clustering of hits.
double getSinThetaC()
Get average sine of Cherenkov angle of water.
Definition: JConstants.hh:155
const double getInverseSpeedOfLight()
Get inverse speed of light.
Definition: JConstants.hh:100
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
Constants.
Template L2 builder.
Definition: JBuildL2.hh:45
void transform(const JRotation3D &R, const JVector3D &pos)
Transform hit.
Definition: JHitL1.hh:346
Detector file.
Definition: JHead.hh:130
Match operator for Cherenkov light from muon in any direction.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1954
double getCosThetaC()
Get average cosine of Cherenkov angle of water.
Definition: JConstants.hh:144
Data time slice.
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
Direct access to PMT in detector data structure.
int debug
debug level
Definition: JSirene.cc:61
then usage $script[distance] fi case set_variable R
Definition: JDrawLED.sh:40
General purpose messaging.
#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.
double getZ() const
Get z position.
Definition: JHitL1.hh:168
Auxiliary class for a time-over-threshold pulse from a PMT.
Definition: JPulse.hh:33
Utility class to parse command line options.
alias put_queue eval echo n
Definition: qlib.csh:19
JDirection3D getDirection(const Vec &v)
Get direction.
Auxiliary class to convert DAQ/trigger hit time to/from Monte Carlo hit time.
int j
Definition: JPolint.hh:634
Data structure for position in three dimensions.
Definition: JPosition3D.hh:35
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
Data structure for normalised vector in three dimensions.
Definition: JVersor3D.hh:23
double getTanThetaC()
Get average tangent of Cherenkov angle of water.
Definition: JConstants.hh:133
General purpose class for multiple pointers.
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition: Trk.hh:12
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:185
double getX() const
Get x position.
Definition: JHitL1.hh:144
int getParentModuleID(const JObjectID &id) const
Get parent module identifier.
Definition: JPMTRouter.hh:174
then usage $script[input file[working directory[option]]] nWhere option can be N
Definition: JMuonPostfit.sh:37
3D match criterion.
Definition: JMatch3D.hh:29
Basic data structure for L1 hit.
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:19
then usage $script[input file[working directory[option]]] nWhere option can be E
Definition: JMuonPostfit.sh:37
3D match criterion with road width.
Definition: JMatch3B.hh:34
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
JPosition3D getPosition(const Vec &v)
Get position.
int main(int argc, char *argv[])
Definition: Main.cpp:15