Jpp test-rotations-new
the software that should make you happy
Loading...
Searching...
No Matches
JSupernova.hh
Go to the documentation of this file.
1#ifndef __JSUPERNOVA_JSUPERNOVA__
2#define __JSUPERNOVA_JSUPERNOVA__
3
4#include <set>
5
10
13
15
16#include "JTools/JRange.hh"
17
20#include "JTrigger/JHitR0.hh"
21#include "JTrigger/JMatchL0.hh"
27
28#include "TH1D.h"
29
32
33/**
34 * \author mlincett,selhedri,iagoos
35 */
36
37namespace JSUPERNOVA {
38
39 using namespace std;
40 using namespace JPP;
41 using namespace KM3NETDAQ;
42
44
46
47 /**
48 * Auxiliary class to store reduced information of a coincidence on an optical module
49 * This class allows storing the observables associated with searches for low-energy
50 * neutrinos from Core-Collapse supernovae:
51 */
53
54 public:
55 double time;
59 // Observables computed using hits in 10ns window in same DOM
61 double mean_dir_ctheta = -2; // mean PMT direction cos(theta)
62 double mean_dir_norm = -1; // hit concentration |R|
63 double deltaT = -1; // mean time spread
64 double total_ToT = -1; // total ToT
65
66 void clear(){
67 this->vetoIndex = 0;
68 this->multiplicity = 0;
69 this->mean_dir_ctheta = -2;
70 this->mean_dir_norm = -1;
71 this->total_ToT = -1;
72 this->deltaT = -1;
73 }
74
75 /**
76 * Constructor:
77 *
78 * \param t Event timing within a timeslice
79 * \param m Number of hits within 10ns
80 * \param dom ID of the module
81 * \param vetoIndex Muon veto index: 0 if pass all vetoes, 1 if pass only trigger veto, 2 if pass only correlation veto, 3 if fails all vetoes
82 * \param ctheta Mean z position of hit cluster on DOM
83 * \param rnorm Hit concentration |R|
84 * \param deltaT Average time spread of hit cluster
85 * \param total_ToT Total ToT in DOM
86 * \param timeslice_time Absolute time of the timeslice
87 */
88 JCoincidenceSN(double t, int m, int dom, int vetoIndex=-1, double ctheta=-2, double rnorm=-1, double deltaT=-1, double total_ToT=-1, double timeslice_time=-1)
90 { }
91
92 /**
93 * Compute single-DOM observables |R|, cos(theta), total ToT, and Delta(t) and set the corresponding class attributes
94 *
95 * \param hits Hit vector from timeslice
96 * \param det Single-DOM detx file used to simulate CCSN signal
97 */
98 bool operator()(const vector<JHitR0>& hits, const JDetector& det){
99 if (hits.size() == 0) return 1;
100 JVector3D mean_pmt_dir;
101 double first_hit_time = 0;
102 this->total_ToT = 0;
103 this->multiplicity = 0;
104 this->deltaT = 0;
105 int hitcount=0;
106 for(auto &hit: hits){
107 this->multiplicity++;
108 this->total_ToT += hit.getToT();
109 if (hitcount > 0) {
110 this->deltaT += hit.getT() - first_hit_time;
111 }
112 else first_hit_time = hit.getT();
113 int channel_id = hit.getPMT();
114 JPMTAddress pmt_address(JModuleRouter(det).getAddress(101),channel_id);
115 auto pmt = det.getPMT(pmt_address).getDirection();
116 mean_pmt_dir += pmt.getDirection();
117 hitcount++;
118 }
119 if (this->multiplicity >= 2) {
120 this->mean_dir_norm = mean_pmt_dir.getLength()/this->multiplicity;
121 this->mean_dir_ctheta = mean_pmt_dir.getZ()/this->multiplicity;
122 this->deltaT /= (this->multiplicity-1);
123 }
124 return 0;
125 }
126
127 int getMultiplicity() const {
128 return multiplicity;
129 }
130
131 int getModule() const {
132 return moduleID;
133 }
134
135 double getTime() const {
136 return time;
137 }
138
139 bool operator<(const JCoincidenceSN& rhs) const {
140 return (time < rhs.time);
141 }
142
143 };
144
145 /**
146 * Auxiliary class to define a veto time window on a set of optical modules
147 */
148 class JVeto {
149
150 private:
153
154 public:
155
156 /**
157 * Default constructor
158 * \param event DAQ event
159 * \param hitRouter hit router as source of hit time calibration
160 *
161 */
162 JVeto(const JDAQEvent& event, const JDAQHitRouter& hitRouter) {
163
165
167
169 hit != event.end<hit_type>();
170 ++hit) {
171
172 moduleSet.insert(hit->getModuleID());
173
174 timeRange.include(getTime(*hit, hitRouter.getPMT(*hit)));
175
176 }
177 }
178
179
180 /**
181 * Get length of veto time range
182 */
183 double getLength() {
184 return timeRange.getLength();
185 }
186
187 /**
188 * Determines if a coincidence is vetoed
189 * \param in coincidence under test
190 */
191 bool operator() (const JCoincidenceSN& in) const {
192 return timeRange(in.getTime()) &&
193 (moduleSet.find(in.getModule()) != moduleSet.end());
194 }
195
196 };
197
198
199 /**
200 * Auxiliary class to build the supernova trigger dataset
201 */
203 : public vector<JCoincidenceSN> {
204
205 private:
207 int min_M;
208 // JDAQHitSelector& hitSelector;// = JDAQHitDefaultSelector();
209
210 public:
213
214 /**
215 * Default constructor
216 * Configures the trigger with a time window and the pretrigger threshold.
217 */
218 JDataSN(double window, int threshold = 4) // JDAQHitSelector& selector = JDAQHitDefaultSelector())
219 : TMax_ns(window), min_M(threshold) // , hitSelector(selector)
220 { };
221
222 /**
223 * Builds coincidences from calibrated hit data and loads them in the internal vector.
224 * \param in hit data
225 * \param moduleID optical module ID for the hit data
226 * \param det single-DOM detx file used to simulate CCSN neutrinos (default: empty detector)
227 */
228 void operator() (const vector<JHitR0>& in, const int& moduleID, const JDetector& det = JDetector()) {
229
230 if (in.size() > 1) {
231
232 for (vector<JHitR0>::const_iterator p = in.begin(); p != in.end(); ) {
233
234 vector<JHitR0>::const_iterator q = p;
235
236 while (++q != in.end() && q->getT() - p->getT() <= TMax_ns ) {}
237
238 int M = distance(p, q);
239
240 if (M >= min_M) {
241 JCoincidenceSN coincidence(p->getT(), M, moduleID);
242 if (det.size() > 0) {
243 vector<JHitR0> coincidence_hits;
244 coincidence_hits.assign(p,q);
245 coincidence(coincidence_hits, det);
246 }
247 this->push_back(coincidence);
248 }
249
250 p = q;
251
252 }
253 }
254 }
255
256
257 /**
258 * Builds coincidences from a timeslice and loads them into the internal vector.
259 * Double pulses are merged according to the coincidence time window.
260 *
261 * \param timeslice input timeslice
262 * \param moduleRouter detector module router
263 * \param selector hit selector
264 * \param det single-DOM detx file used to simulate signal events
265 */
266 void operator() (const JDAQTimeslice* timeslice, const JModuleRouter& moduleRouter, const JDAQHitSelector& selector = JDAQHitDefaultSelector(), const JDetector& det = JDetector()) {
267
268 frameIndex = timeslice->getFrameIndex();
269 timeUTC = timeslice->getTimesliceStart();
270
271 typedef JHitR0 hit_t;
272
273 typedef JSuperFrame2D<hit_t> JSuperFrame2D_t;
274
275 typedef JSuperFrame1D<hit_t> JSuperFrame1D_t;
276
277 const JMatchL0<hit_t> match(TMax_ns);
278
279 for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
280
281 int moduleID = frame->getModuleID();
282
283 JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*frame,
284 moduleRouter.getModule(moduleID),
285 selector);
286
287 buffer.preprocess(JPreprocessor::join_t, match);
288
289 JSuperFrame1D_t& data = JSuperFrame1D_t::multiplex(buffer);
290
291 if (det.size() < 0)(*this)(data, moduleID);
292 else (*this)(data, moduleID, det);
293
294 }
295
296 sort(this->begin(), this->end());
297
298 }
299
300 };
301
302
303 /**
304 * Auxiliary class-operator to match a JVeto with a JCoincidenceSN object
305 * Provides an operator to test if a coincidence is vetoed
306 */
307
309
310 private:
312
313 public:
314 /**
315 * Default constructor
316 * \param in coincidence to be matched against veto
317 */
318 JMatchVeto(const JCoincidenceSN& in) : dut(in) {}
319
320 /**
321 * Operator
322 * \param in veto to be matched against inner coincidence
323 */
324 bool operator() (const JVeto& in) {
325 return in(dut);
326 }
327
328 };
329
330
331 /**
332 * Auxiliary class to manage a set of vetoes
333 */
334 class JVetoSet : public vector<JVeto> {
335
336 public:
337 /**
338 * Applies the veto set to a coincidence
339 * \param in coincidence to be tested
340 */
341 bool operator() (const JCoincidenceSN& in) const {
342 return any_of(this->begin(), this->end(), JMatchVeto(in));
343 }
344
345 };
346
347
348 /**
349 * Auxiliary class to manage a cluster of coincidences
350 */
351 class JClusterSN : public vector<JCoincidenceSN> {
352
353 public:
354
355 /**
356 * Finds coincidence with maximum multiplicity.
357 */
359
360 JClusterSN::const_iterator p = this->begin();
361
362 for (JClusterSN::const_iterator q = p + 1; q != this->end(); q++) {
363 if (q->getMultiplicity() > p->getMultiplicity()) {
364 p = q;
365 }
366 }
367
368 return (*p);
369
370 }
371
372
373 /*
374 * Returns the set of modules spanned over by the cluster.
375 */
377
378 JModuleSet out;
379
380 for (JClusterSN::const_iterator p = this->begin(); p != this->end(); p++) {
381 out.insert(p->getModule());
382 }
383
384 return out;
385 }
386
387 };
388
389
390 /**
391 * Interface for SN filter operator.
392 * This operator is meant to determine the SN trigger count from a vector of JClusterSN
393 */
394 class JSNFilter {
395 public:
396 virtual bool operator() (const JCoincidenceSN& in) const = 0;
397 virtual bool operator() (const JClusterSN& in) const = 0;
398 };
399
400
401 /**
402 * SN filter based on multiplicity selection
403 * optional suppression of multi-module coincidences
404 * WARNING: no minimum threshold for the veto
405 */
406 class JSNFilterM : public JSNFilter {
407
408 private:
410 bool mode;
411
412 public:
414 : A(R), mode(m)
415 {}
416
417 // select coincidence if within multiplicity acceptance
418 bool operator() (const JCoincidenceSN& in) const {
419 return A(in.getMultiplicity());
420 }
421
422 // select cluster if any coincidence is accepted
423 // optionally veto if more of one module is interested
424 bool operator() (const JClusterSN& in) const {
425
426 bool out = (*this)(in.getPeak());
427
428 if (mode == 1) {
429 out = out && (in.getModules().size() == 1);
430 }
431
432 return out;
433 }
434
435 };
436
437
438 /**
439 * SN filter based on veto window
440 */
441 class JSNFilterMV : public JSNFilter {
442
443 private:
446
447 public:
449 : A(R), V(S)
450 {}
451
452 bool operator() (const JCoincidenceSN& in) const {
453 return A(in.getMultiplicity()) && !V(in);
454 }
455
456 bool operator() (const JClusterSN& in) const {
457 return (*this)(in.getPeak());
458 // return any_of(in.begin(), in.end(), *this);
459 }
460
461 };
462
463
464
465 /**
466 * Auxiliary class to apply the supernova trigger to SN data
467 */
468 class JTriggerSN : public vector<JClusterSN> {
469
470 private:
471
472 double TMaxCluster_ns = 1000.0;
473
474 // JVetoSet veto;
475
476 public:
477
480
481 /**
482 * Default constructor
483 *
484 * \param TMax_ns time width for coincidence clustering (track / afterpulse)
485 */
486 JTriggerSN(double TMax_ns) :
487 TMaxCluster_ns(TMax_ns)
488 {};
489
490 // void setVeto(JVetoSet &vt) {
491 // veto = vt;
492 // }
493
494
495 /**
496 * Builds trigger by clustering pretrigger data
497 *
498 * \param in pretrigger SN data
499 */
500 void operator() (const JDataSN& in) {
501
503 timeUTC = in.timeUTC;
504
505 for (vector<JCoincidenceSN>::const_iterator p = in.begin(); p != in.end(); ) {
506
507 JClusterSN cluster;
508
509 cluster.push_back(*p);
510
511 vector<JCoincidenceSN>::const_iterator q = p + 1;
512
513 while(q != in.end() && (q->getTime() <= (p->getTime() + TMaxCluster_ns))) {
514 cluster.push_back(*q);
515 ++q;
516 }
517
518 p = q;
519
520 this->push_back(cluster);
521
522 }
523 }
524
525
526 /**
527 * Get triggered modules after track veto
528 *
529 * \return std::set containing triggered modules IDs
530 */
532
533 JModuleSet triggeredModules;
534
535 JSNFilterM F(A, 1);
536
537 for (JTriggerSN::const_iterator p = this->begin(); p != this->end(); p++) {
538
539 if ( F(*p) ) {
540
541 // only clusters with one module pass the selection JSNFilterM(A, 1)
542 triggeredModules.insert((*p)[0].getModule());
543
544 }
545
546 }
547
548 return triggeredModules;
549
550 }
551
552
553 /**
554 * Get triggered modules according to given filter
555 */
557
558 JModuleSet out;
559
560 for (JTriggerSN::const_iterator p = this->begin(); p != this->end(); p++) {
561 if ( F(*p) ) {
562 out.insert(p->getPeak().getModule());
563 }
564 }
565 return out;
566 }
567
568 /**
569 * Fill histogram with multiplicity spectrum resulting from given filter
570 */
571 void fill(TH1D* out, const JSNFilter& F) {
572 for (JTriggerSN::const_iterator p = this->begin(); p != this->end(); p++) {
573 if (F(*p)) {
574 out->Fill( p->getPeak().getMultiplicity() );
575 }
576 }
577 }
578
580 vector<double> out;
581
582 for (JTriggerSN::const_iterator p = this->begin(); p != this->end(); p++) {
583 if (F(*p)) {
584 out.push_back( p->getPeak().getMultiplicity() );
585 }
586 }
587 return out;
588 }
589
590
591
592
593 /**
594 * > operator
595 * used by (reverse) std:priority_queue, in absence of this operator the container will behave unpredictably
596 */
597 bool operator>(const JTriggerSN& rhs) const {
598 return (frameIndex > rhs.frameIndex);
599 }
600
601 /**
602 * < operator
603 */
604 bool operator<(const JTriggerSN& rhs) const {
605 return (frameIndex < rhs.frameIndex);
606 }
607
608
609 };
610
611
612 /**
613 * SN trigger statistics, the information is stored in the form of a count as a function of the trigger level.
614 * the livetime needs to be set manually to compute the rates for the printing.
615 */
616
617 class JTriggerSNStats : public vector<double> {
618
619 private:
620 double livetime;
621
622 public:
623 /**
624 * default constructor
625 *
626 * \param nModules number of modules of the detector
627 */
628 JTriggerSNStats(const int nModules) : vector<double>(nModules) {}
629
630 void setLiveTime(const double lt) {
631 livetime = lt;
632 }
633
634 /**
635 * put statistics into printable form
636 * outputs trigger level - rate - error
637 */
638 string toString() {
639
640 ostringstream os;
641
642 os << "=== SUPERNOVA TRIGGER STATISTICS ==" << endl;
643 os << "=> livetime [s] = " << livetime << endl;
644 os << "Level" << '\t' << "Rate (error)" << endl;
645
646 if (livetime > 0) {
647 for (unsigned i = 0; i < this->size(); i++) {
648
649 double count = (*this)[i];
650 double rate = count / livetime;
651 double error = 100 * (sqrt(count) / count);
652
653 if (rate > 0) {
654 os << i << '\t';
655 os << scientific << setprecision(2) << rate << "Hz ";
656 os << fixed << setprecision(0) << "(" << setw(2) << error << "%)";
657 os << endl;
658 }
659 }
660 }
661
662 return os.str();
663
664 }
665
666
667 /**
668 * put statistics into printable form
669 * outputs trigger level - rate - error
670 */
671 string toSummaryFile() {
672
673 ostringstream os;
674
675 if (livetime > 0) {
676
677 os << livetime << endl;
678
679 for (unsigned i = 0; i < this->size(); i++) {
680
681 double count = (*this)[i];
682 double rate = count / livetime;
683 double error = (sqrt(count) / count);
684
685 if (rate > 0) {
686 os << i << ",";
687 os << scientific;
688 os << setprecision(2);
689 os << rate << ",";
690 os << error;
691 os << endl;
692 }
693 }
694 }
695
696
697 return os.str();
698
699 }
700
701 };
702
703}
704
705#endif
Direct access to PMT data in detector data structure for DAQ hits.
Data structure for detector geometry and calibration.
Basic data structure for L0 hit.
Match operator for consecutive hits.
Direct access to module in detector data structure.
Auxiliaries for pre-processing of hits.
Auxiliary class to define a range between two values.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Simple wrapper around JModuleRouter class for direct addressing of PMT data in detector data structur...
const JPMT & getPMT(const JDAQKeyHit &hit) const
Get PMT parameters.
Detector data structure.
Definition JDetector.hh:96
const JPMT & getPMT(const JPMTAddress &address) const
Get PMT parameters.
Definition JDetector.hh:348
Router for direct addressing of module data in detector data structure.
const JModule & getModule(const JObjectID &id) const
Get module parameters.
Address of PMT in detector data structure.
const JDirection3D & getDirection() const
Get direction.
Data structure for vector in three dimensions.
Definition JVector3D.hh:36
double getLength() const
Get length.
Definition JVector3D.hh:246
double getZ() const
Get z position.
Definition JVector3D.hh:115
Auxiliary class to manage a cluster of coincidences.
JCoincidenceSN getPeak() const
Finds coincidence with maximum multiplicity.
JModuleSet getModules() const
Auxiliary class to store reduced information of a coincidence on an optical module This class allows ...
Definition JSupernova.hh:52
JCoincidenceSN(double t, int m, int dom, int vetoIndex=-1, double ctheta=-2, double rnorm=-1, double deltaT=-1, double total_ToT=-1, double timeslice_time=-1)
Constructor:
Definition JSupernova.hh:88
bool operator()(const vector< JHitR0 > &hits, const JDetector &det)
Compute single-DOM observables |R|, cos(theta), total ToT, and Delta(t) and set the corresponding cla...
Definition JSupernova.hh:98
bool operator<(const JCoincidenceSN &rhs) const
Auxiliary class to build the supernova trigger dataset.
JDAQUTCExtended timeUTC
JDataSN(double window, int threshold=4)
Default constructor Configures the trigger with a time window and the pretrigger threshold.
void operator()(const vector< JHitR0 > &in, const int &moduleID, const JDetector &det=JDetector())
Builds coincidences from calibrated hit data and loads them in the internal vector.
Auxiliary class-operator to match a JVeto with a JCoincidenceSN object Provides an operator to test i...
bool operator()(const JVeto &in)
Operator.
JMatchVeto(const JCoincidenceSN &in)
Default constructor.
SN filter based on veto window.
JSNFilterMV(JRange< int > &R, JVetoSet &S)
bool operator()(const JCoincidenceSN &in) const
SN filter based on multiplicity selection optional suppression of multi-module coincidences WARNING: ...
JSNFilterM(JRange< int > R, int m=0)
bool operator()(const JCoincidenceSN &in) const
Interface for SN filter operator.
virtual bool operator()(const JCoincidenceSN &in) const =0
SN trigger statistics, the information is stored in the form of a count as a function of the trigger ...
string toString()
put statistics into printable form outputs trigger level - rate - error
void setLiveTime(const double lt)
JTriggerSNStats(const int nModules)
default constructor
string toSummaryFile()
put statistics into printable form outputs trigger level - rate - error
Auxiliary class to apply the supernova trigger to SN data.
void fill(TH1D *out, const JSNFilter &F)
Fill histogram with multiplicity spectrum resulting from given filter.
void operator()(const JDataSN &in)
Builds trigger by clustering pretrigger data.
JDAQUTCExtended timeUTC
JTriggerSN(double TMax_ns)
Default constructor.
JModuleSet getModules(const JSNFilter &F)
Get triggered modules according to given filter.
JModuleSet getModules(JRange< int > A=JRange< int >(6, 10))
Get triggered modules after track veto.
vector< double > getMultiplicities(const JSNFilter &F)
bool operator<(const JTriggerSN &rhs) const
< operator
bool operator>(const JTriggerSN &rhs) const
Auxiliary class to manage a set of vetoes.
bool operator()(const JCoincidenceSN &in) const
Applies the veto set to a coincidence.
Auxiliary class to define a veto time window on a set of optical modules.
double getLength()
Get length of veto time range.
bool operator()(const JCoincidenceSN &in) const
Determines if a coincidence is vetoed.
JTimeRange timeRange
JVeto(const JDAQEvent &event, const JDAQHitRouter &hitRouter)
Default constructor.
JModuleSet moduleSet
T getLength() const
Get length (difference between upper and lower limit).
Definition JRange.hh:289
static JRange< double, std::less< double > > DEFAULT_RANGE()
Definition JRange.hh:555
range_type & include(argument_type x)
Include given value to range.
Definition JRange.hh:397
Reduced data structure for L0 hit.
Definition JHitR0.hh:27
L0 match criterion.
Definition JMatchL0.hh:29
1-dimensional frame with time calibrated data from one optical module.
2-dimensional frame with time calibrated data from one optical module.
JDAQUTCExtended getTimesliceStart() const
Get start of timeslice.
int getFrameIndex() const
Get frame index.
Template const_iterator.
Definition JDAQEvent.hh:68
const_iterator< T > end() const
Get end of data.
const_iterator< T > begin() const
Get begin of data.
Data structure for UTC time.
const JModule & getModule(const JType< JDetector_t > type, const int id, const JLocation &location=JLocation())
Get module according given detector type.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
set< int > JModuleSet
Definition JSupernova.hh:45
vector< double > multiplicities_t
Definition JSupernova.hh:43
const char * getTime()
Get current local time conform ISO-8601 standard.
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
Normalisation of MUPAGE events.
Definition JHead.hh:835
Acoustic hit.
Definition JBillabong.cc:70
Default class to select DAQ hits.
Auxiliary class to select DAQ hits.
@ join_t
join consecutive hits according match criterion