#ifndef __JSUPERNOVA_JSUPERNOVA__ #define __JSUPERNOVA_JSUPERNOVA__ #include #include "km3net-dataformat/online/JDAQEvent.hh" #include "km3net-dataformat/online/JDAQTriggeredHit.hh" #include "km3net-dataformat/online/JDAQTimeslice.hh" #include "km3net-dataformat/online/JDAQUTCExtended.hh" #include "JDetector/JDAQHitRouter.hh" #include "JDetector/JModuleRouter.hh" #include "JGeometry3D/JAxis3D.hh" #include "JTools/JRange.hh" #include "JTrigger/JDAQHitSelector.hh" #include "JTrigger/JEventToolkit.hh" #include "JTrigger/JHitR0.hh" #include "JTrigger/JMatchL0.hh" #include "JTrigger/JPreprocessor.hh" #include "JTrigger/JSuperFrame1D.hh" #include "JTrigger/JSuperFrame2D.hh" #include "JTrigger/JTriggerParameters.hh" #include "JSupport/JTriggerParametersSupportkit.hh" #include "TH1D.h" /** * \author mlincett */ namespace JSUPERNOVA { using namespace std; using namespace JPP; using namespace KM3NETDAQ; typedef vector multiplicities_t; typedef set JModuleSet; /** * Auxiliary class to store reduced information of a coincidence on an optical module */ class JCoincidenceSN { private: double time; int multiplicity; int moduleID; public: JCoincidenceSN(double t, int m, int dom) : time(t), multiplicity(m), moduleID(dom) { } int getMultiplicity() const { return multiplicity; } int getModule() const { return moduleID; } double getTime() const { return time; } bool operator<(const JCoincidenceSN& rhs) const { return (time < rhs.time); } }; /** * Auxiliary class to define a veto time window on a set of optical modules */ class JVeto { private: JTimeRange timeRange; JModuleSet moduleSet; public: /** * Default constructor * \param event DAQ event * \param hitRouter hit router as source of hit time calibration * */ JVeto(const JDAQEvent& event, const JDAQHitRouter& hitRouter) { timeRange = JTimeRange::DEFAULT_RANGE; typedef JDAQTriggeredHit hit_type; for (JDAQEvent::const_iterator hit = event.begin(); hit != event.end(); ++hit) { moduleSet.insert(hit->getModuleID()); timeRange.include(getTime(*hit, hitRouter.getPMT(*hit))); } } /** * Get length of veto time range */ double getLength() { return timeRange.getLength(); } /** * Determines if a coincidence is vetoed * \param in coincidence under test */ bool operator() (const JCoincidenceSN& in) const { return timeRange(in.getTime()) && (moduleSet.find(in.getModule()) != moduleSet.end()); } }; /** * Auxiliary class to build the supernova trigger dataset */ class JDataSN : public vector { private: int TMax_ns; int min_M; // JDAQHitSelector& hitSelector;// = JDAQHitDefaultSelector(); public: int frameIndex; JDAQUTCExtended timeUTC; /** * Default constructor * Configures the trigger with a time window and the pretrigger threshold. */ JDataSN(double window, int threshold = 4) // JDAQHitSelector& selector = JDAQHitDefaultSelector()) : TMax_ns(window), min_M(threshold) // , hitSelector(selector) { }; /** * Builds coincidences from calibrated hit data and loads them in the internal vector. * \param in hit data * \param moduleID optical module ID for the hit data */ void operator() (vector in, int moduleID) { if (in.size() > 1) { for (vector::const_iterator p = in.begin(); p != in.end(); ) { vector::const_iterator q = p; while (++q != in.end() && q->getT() - p->getT() <= TMax_ns ) {} int M = distance(p, q); if (M >= min_M) { this->push_back(JCoincidenceSN(p->getT(), M, moduleID)); } p = q; } } } /** * Builds coincidences from a timeslice and loads them into the internal vector. * Double pulses are merged according to the coincidence time window. * * \param timeslice input timeslice * \param moduleRouter detector module router * \param selector hit selector */ void operator() (const JDAQTimeslice* timeslice, const JModuleRouter& moduleRouter, const JDAQHitSelector& selector = JDAQHitDefaultSelector()) { frameIndex = timeslice->getFrameIndex(); timeUTC = timeslice->getTimesliceStart(); typedef JHitR0 hit_t; typedef JSuperFrame2D JSuperFrame2D_t; typedef JSuperFrame1D JSuperFrame1D_t; const JMatchL0 match(TMax_ns); for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) { int moduleID = frame->getModuleID(); JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*frame, moduleRouter.getModule(moduleID), selector); buffer.preprocess(JPreprocessor::join_t, match); JSuperFrame1D_t& data = JSuperFrame1D_t::multiplex(buffer); (*this)(data, moduleID); } sort(this->begin(), this->end()); } }; /** * Auxiliary class-operator to match a JVeto with a JCoincidenceSN object * Provides an operator to test if a coincidence is vetoed */ class JMatchVeto { private: JCoincidenceSN dut; public: /** * Default constructor * \param in coincidence to be matched against veto */ JMatchVeto(const JCoincidenceSN& in) : dut(in) {} /** * Operator * \param in veto to be matched against inner coincidence */ bool operator() (const JVeto& in) { return in(dut); } }; /** * Auxiliary class to manage a set of vetoes */ class JVetoSet : public vector { public: /** * Applies the veto set to a coincidence * \param in coincidence to be tested */ bool operator() (const JCoincidenceSN& in) const { return any_of(this->begin(), this->end(), JMatchVeto(in)); } }; /** * Auxiliary class to manage a cluster of coincidences */ class JClusterSN : public vector { public: /** * Finds coincidence with maximum multiplicity. */ JCoincidenceSN getPeak() const { JClusterSN::const_iterator p = this->begin(); for (JClusterSN::const_iterator q = p + 1; q != this->end(); q++) { if (q->getMultiplicity() > p->getMultiplicity()) { p = q; } } return (*p); } /* * Returns the set of modules spanned over by the cluster. */ JModuleSet getModules() const { JModuleSet out; for (JClusterSN::const_iterator p = this->begin(); p != this->end(); p++) { out.insert(p->getModule()); } return out; } }; /** * Interface for SN filter operator. * This operator is meant to determine the SN trigger count from a vector of JClusterSN */ class JSNFilter { public: virtual bool operator() (const JCoincidenceSN& in) const = 0; virtual bool operator() (const JClusterSN& in) const = 0; }; /** * SN filter based on multiplicity selection * optional suppression of multi-module coincidences * WARNING: no minimum threshold for the veto */ class JSNFilterM : public JSNFilter { private: JRange A; bool mode; public: JSNFilterM(JRange R, int m = 0) : A(R), mode(m) {} // select coincidence if within multiplicity acceptance bool operator() (const JCoincidenceSN& in) const { return A(in.getMultiplicity()); } // select cluster if any coincidence is accepted // optionally veto if more of one module is interested bool operator() (const JClusterSN& in) const { bool out = (*this)(in.getPeak()); if (mode == 1) { out = out && (in.getModules().size() == 1); } return out; } }; /** * SN filter based on veto window */ class JSNFilterMV : public JSNFilter { private: JRange A; JVetoSet V; public: JSNFilterMV(JRange& R, JVetoSet& S) : A(R), V(S) {} bool operator() (const JCoincidenceSN& in) const { return A(in.getMultiplicity()) && !V(in); } bool operator() (const JClusterSN& in) const { return (*this)(in.getPeak()); // return any_of(in.begin(), in.end(), *this); } }; /** * Auxiliary class to apply the supernova trigger to SN data */ class JTriggerSN : public vector { private: double TMaxCluster_ns = 1000.0; // JVetoSet veto; public: int frameIndex; JDAQUTCExtended timeUTC; /** * Default constructor * * \param TMax_ns time width for coincidence clustering (track / afterpulse) */ JTriggerSN(double TMax_ns) : TMaxCluster_ns(TMax_ns) {}; // void setVeto(JVetoSet &vt) { // veto = vt; // } /** * Builds trigger by clustering pretrigger data * * \param in pretrigger SN data */ void operator() (const JDataSN& in) { frameIndex = in.frameIndex; timeUTC = in.timeUTC; for (vector::const_iterator p = in.begin(); p != in.end(); ) { JClusterSN cluster; cluster.push_back(*p); vector::const_iterator q = p + 1; while(q != in.end() && (q->getTime() <= (p->getTime() + TMaxCluster_ns))) { cluster.push_back(*q); ++q; } p = q; this->push_back(cluster); } } /** * Get triggered modules after track veto * * \return std::set containing triggered modules IDs */ JModuleSet getModules(JRange A = JRange(6,10)) { JModuleSet triggeredModules; JSNFilterM F(A, 1); for (JTriggerSN::const_iterator p = this->begin(); p != this->end(); p++) { if ( F(*p) ) { // only clusters with one module pass the selection JSNFilterM(A, 1) triggeredModules.insert((*p)[0].getModule()); } } return triggeredModules; } /** * Get triggered modules according to given filter */ JModuleSet getModules(const JSNFilter& F) { JModuleSet out; for (JTriggerSN::const_iterator p = this->begin(); p != this->end(); p++) { if ( F(*p) ) { out.insert(p->getPeak().getModule()); } } return out; } /** * Fill histogram with multiplicity spectrum resulting from given filter */ void fill(TH1D* out, const JSNFilter& F) { for (JTriggerSN::const_iterator p = this->begin(); p != this->end(); p++) { if (F(*p)) { out->Fill( p->getPeak().getMultiplicity() ); } } } vector getMultiplicities(const JSNFilter& F) { vector out; for (JTriggerSN::const_iterator p = this->begin(); p != this->end(); p++) { if (F(*p)) { out.push_back( p->getPeak().getMultiplicity() ); } } return out; } /** * > operator * used by (reverse) std:priority_queue, in absence of this operator the container will behave unpredictably */ bool operator>(const JTriggerSN& rhs) const { return (frameIndex > rhs.frameIndex); } /** * < operator */ bool operator<(const JTriggerSN& rhs) const { return (frameIndex < rhs.frameIndex); } }; /** * SN trigger statistics, the information is stored in the form of a count as a function of the trigger level. * the livetime needs to be set manually to compute the rates for the printing. */ class JTriggerSNStats : public vector { private: double livetime; public: /** * default constructor * * \param nModules number of modules of the detector */ JTriggerSNStats(const int nModules) : vector(nModules) {} void setLiveTime(const double lt) { livetime = lt; } /** * put statistics into printable form * outputs trigger level - rate - error */ string toString() { ostringstream os; os << "=== SUPERNOVA TRIGGER STATISTICS ==" << endl; os << "=> livetime [s] = " << livetime << endl; os << "Level" << '\t' << "Rate (error)" << endl; if (livetime > 0) { for (unsigned i = 0; i < this->size(); i++) { double count = (*this)[i]; double rate = count / livetime; double error = 100 * (sqrt(count) / count); if (rate > 0) { os << i << '\t'; os << scientific << setprecision(2) << rate << "Hz "; os << fixed << setprecision(0) << "(" << setw(2) << error << "%)"; os << endl; } } } return os.str(); } /** * put statistics into printable form * outputs trigger level - rate - error */ string toSummaryFile() { ostringstream os; if (livetime > 0) { os << livetime << endl; for (unsigned i = 0; i < this->size(); i++) { double count = (*this)[i]; double rate = count / livetime; double error = (sqrt(count) / count); if (rate > 0) { os << i << ","; os << scientific; os << setprecision(2); os << rate << ","; os << error; os << endl; } } } return os.str(); } }; } #endif