Jpp 19.3.0-rc.2
the software that should make you happy
Loading...
Searching...
No Matches
Functions
JHalibut.cc File Reference

Example program to determine N-fold coincidence rates. More...

#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <vector>
#include <map>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TF1.h"
#include "JROOT/JMinimizer.hh"
#include "JDAQ/JDAQTimesliceIO.hh"
#include "km3net-dataformat/online/JDAQClock.hh"
#include "JDAQ/JDAQEvaluator.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JTrigger/JHitR0.hh"
#include "JTrigger/JMatchL0.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JSuperFrame1D.hh"
#include "JTrigger/JSuperFrame2D.hh"
#include "JROOT/JManager.hh"
#include "JCalibrate/JCalibrateK40.hh"
#include "JROOT/JROOTClassSelector.hh"
#include "JTools/JQuantile.hh"
#include "JTools/JRange.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JAutoTreeScanner.hh"
#include "JSupport/JSupport.hh"
#include "JLang/JObjectMultiplexer.hh"
#include "JMath/JMathToolkit.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Example program to determine N-fold coincidence rates.

Author
mdejong

Definition in file JHalibut.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 54 of file JHalibut.cc.

55{
56 using namespace std;
57 using namespace JPP;
58 using namespace KM3NETDAQ;
59
60 typedef JRange<int> JRange_t;
61
62 JMultipleFileScanner<> inputFile;
63 JLimit_t& numberOfEvents = inputFile.getLimit();
64 string outputFile;
65 string detectorFile;
66 double TMax_ns;
67 JROOTClassSelector selector;
68 JRange_t M;
69 string summaryFile;
70 string option;
71 int debug;
72
73 try {
74
75 JParser<> zap("Example program to determine N-fold coincidence rates.");
76
77 zap['f'] = make_field(inputFile);
78 zap['o'] = make_field(outputFile) = "halibut.root";
79 zap['n'] = make_field(numberOfEvents) = JLimit::max();
80 zap['a'] = make_field(detectorFile);
81 zap['T'] = make_field(TMax_ns) = 20.0;
82 zap['C'] = make_field(selector) = getROOTClassSelection<JDAQTimesliceTypes_t>();
83 zap['M'] = make_field(M) = JRange_t(2, 31);
84 zap['s'] = make_field(summaryFile) = "halibut.txt";
85 zap['O'] = make_field(option) = "";
86 zap['d'] = make_field(debug) = 1;
87
88 zap(argc, argv);
89 }
90 catch(const exception &error) {
91 FATAL(error.what() << endl);
92 }
93
94
95 JDetector detector;
96
97 try {
98 load(detectorFile, detector);
99 }
100 catch(const JException& error) {
101 FATAL(error);
102 }
103
104 const JModuleRouter router(detector);
105
106 map<int, double> TMax_s; // histogram upper limit
107
108 TMax_s[2] = 5.0e-2;
109 TMax_s[3] = 0.2e+0;
110 TMax_s[4] = 1.0e+0;
111 TMax_s[5] = 1.0e+1;
112 TMax_s[6] = 1.0e+2;
113
114 typedef JSuperFrame2D<JHitR0> JSuperFrame2D_t;
115 typedef JSuperFrame1D<JHitR0> JSuperFrame1D_t;
116
117 const JMatchL0<JHitR0> match(TMax_ns); // time window self-coincidences [ns]
118
119 typedef JManager<int, TH1D> JManager_t;
120
121 JManager_t H1(new TH1D("H1[%]", NULL, 100, -TMax_ns, +TMax_ns));
122 JManager_t T1(new TH1D("T1[%]", NULL, 100, 0.0, TMax_s[M.getLowerLimit()]));
123
124 H1->Sumw2();
125 T1->Sumw2();
126
127 JTreeScanner<>::debug = debug;
128
129 JAutoTreeScanner<JDAQTimeslice, JDAQEvaluator> zmap = JType<JDAQTimesliceTypes_t>();
130
131 JTreeScannerInterface<JDAQTimeslice>* ps = zmap[selector];
132
133 counter_type counter = 0;
134
135 // process file-by-file to speed up JTreeScanner::configure();
136
137 for (JMultipleFileScanner<>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
138
139 STATUS("Processing " << *i << endl);
140
141 ps->configure(*i);
142
143 vector<double> t0(detector.size(), 0.0);
144
145 for ( ; ps->hasNext() && counter != inputFile.getLimit(); ++counter) {
146
147 STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl);
148
149 const JDAQTimeslice* timeslice = ps->next();
150
151 for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
152
153 JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*frame, router.getModule(frame->getModuleID()));
154
155 for (JSuperFrame2D_t::iterator i = buffer.begin(); i != buffer.end(); ++i) {
156 i->join(match);
157 }
158
159 JSuperFrame1D_t& data = JSuperFrame1D_t::multiplex(buffer);
160
161 if (data.size() > 1) {
162
163 TH1D* h1 = H1[frame->getModuleID()];
164 TH1D* t1 = T1[frame->getModuleID()];
165
166 for (vector<JHitR0>::const_iterator p = data.begin(); p != data.end(); ) {
167
168 vector<JHitR0>::const_iterator q = p;
169
170 while (++q != data.end() && q->getT() - p->getT() <= TMax_ns ) {}
171
172 const int N = distance(p,q);
173
174 if (M(N)) {
175
176 const int i = router.getIndex(frame->getModuleID());
177 const double ts = getTimeOfRTS(frame->getFrameIndex()) + p->getT();
178
179 t1->Fill((ts - t0[i]) * 1.0e-9); // [s]
180
181 t0[i] = ts;
182
183 const double W = factorial(N, 2);
184
185 for (vector<JHitR0>::const_iterator __p = p; __p != q; ++__p) {
186 for (vector<JHitR0>::const_iterator __q = __p; ++__q != q; ) {
187 h1->Fill(JCombinatorics::getSign(__p->getPMT(),__q->getPMT()) * (__q->getT() - __p->getT()), 1.0/W);
188 }
189 }
190 }
191
192 p = q;
193 }
194 }
195 }
196 }
197 STATUS(endl);
198 }
199
200 if (counter != 0) {
201
202 const double V = (H1->GetXaxis()->GetXmax() - H1->GetXaxis()->GetXmin()) / (double) H1->GetXaxis()->GetNbins(); // [ns]
203 const double W = counter * getFrameTime() * 1.0e-9; // [s]
204
205 for (JManager_t::iterator i = H1.begin(); i != H1.end(); ++i) {
206 i->second->Scale(1.0/(V*W));
207 }
208
209 for (JManager_t::iterator i = T1.begin(); i != T1.end(); ++i) {
210 i->second->Scale(1.0/i->second->GetMaximum());
211 }
212 }
213
214 if (summaryFile != "") {
215
216 const double V = (H1->GetXaxis()->GetXmax() - H1->GetXaxis()->GetXmin()) / (double) H1->GetXaxis()->GetNbins(); // [ns]
217
218 const int number_of_strings = getNumberOfStrings(detector);
219 const int number_of_floors = getNumberOfFloors (detector);
220 const int PRECISION = (M.getLowerLimit() > 2 ? 4 : 3);
221
222 ofstream out(summaryFile.c_str());
223
224 out << "Multiplicity " << M << endl;
225 out << "-------------------------------------------------------" << endl;
226 out << " location | Gauss | S - B | Total | slope " << endl;
227 out << " | [Hz] | [Hz] | [Hz] | [Hz] " << endl;
228 out << "-------------------------------------------------------" << endl;
229
230 JQuantile Q[4];
231
232 for (int string = 1; string <= number_of_strings; ++string) {
233 for (int floor = number_of_floors; floor >= 1; --floor) {
234
235 const int id = detector.getModule(JLocation(string,floor)).getID();
236
237 out << " " << setw(3) << string << ' ' << setw(2) << floor << " ";
238
239 TH1D* h1 = (H1.find(id) != H1.end() ? H1[id] : NULL);
240 TH1D* t1 = (T1.find(id) != T1.end() ? T1[id] : NULL);
241
242 if (h1 != NULL) {
243
244 TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))/(TMath::Sqrt(2*TMath::Pi())*[2]) + [3]");
245
246 f1.SetParameter(0, h1->GetMaximum());
247 f1.SetParameter(1, 0.0);
248 f1.SetParameter(2, h1->GetRMS() * 0.25);
249 f1.SetParameter(3, h1->GetMinimum());
250
251 h1->Fit(&f1, option.c_str(), "same");
252
253 out << " | " << FIXED(8,PRECISION) << f1.GetParameter(0);
254 out << " | " << FIXED(8,PRECISION) << (h1->GetSumOfWeights() - f1.GetParameter(3) * h1->GetNbinsX()) * V;
255 out << " | " << FIXED(8,PRECISION) << h1->GetSumOfWeights() * V;
256
257 Q[0].put( f1.GetParameter(0));
258 Q[1].put((h1->GetSumOfWeights() - f1.GetParameter(3) * h1->GetNbinsX()) * V);
259 Q[2].put( h1->GetSumOfWeights() * V);
260 }
261
262 if (t1 != NULL) {
263
264 TF1 f1("f1", "[0]*exp(-[1]*x)");
265
266 f1.SetParameter(0, t1->GetMaximum());
267 f1.SetParameter(1, 1.0 / t1->GetRMS());
268
269 t1->Fit(&f1, option.c_str(), "same");
270
271 out << " | " << FIXED(8,PRECISION) << f1.GetParameter(1);
272
273 Q[3].put(f1.GetParameter(1));
274 }
275
276 out << endl;
277 }
278
279 out << endl;
280 }
281
282 if (Q[0].getCount() != 0) {
283
284 out << "-------------------------------------------------------" << endl;
285 out << setw(10) << left << " average";
286
287 for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
288 out << " | " << FIXED(8,PRECISION) << Q[i].getMean();
289 }
290
291 out << endl;
292 }
293
294 out.close();
295 }
296
297 if (outputFile != "") {
298
299 TFile out(outputFile.c_str(), "RECREATE");
300
301 H1.Write(out);
302 T1.Write(out);
303
304 out.Close();
305 }
306}
outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37
DEBUG
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
STATUS
#define STATUS(A)
Definition JMessage.hh:63
FATAL
#define FATAL(A)
Definition JMessage.hh:67
debug
int debug
debug level
Definition JSirene.cc:72
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
distance
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Definition PhysicsEvent.hh:434
JDETECTOR::JDetector
Detector data structure.
Definition JDetector.hh:96
JDETECTOR::JLocation
Logical location of module.
Definition JLocation.hh:40
JDETECTOR::JModuleRouter
Router for direct addressing of module data in detector data structure.
Definition JModuleRouter.hh:38
JLANG::JException
General exception.
Definition JException.hh:24
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
JROOT::JManager
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys.
Definition JManager.hh:47
JSUPPORT::JTreeScannerInterface
Auxiliary interface for direct access of elements in ROOT TChain.
Definition JTreeScannerInterface.hh:412
JSUPPORT::JTreeScanner
Template definition for direct access of elements in ROOT TChain.
Definition JTreeScanner.hh:698
JTOOLS::JCombinatorics::getSign
static int getSign(const int first, const int second)
Sign of pair of indices.
Definition JCombinatorics.hh:215
JTOOLS::JRange
Range of values.
Definition JRange.hh:42
JTOOLS::JRange::getLowerLimit
T getLowerLimit() const
Get lower limit.
Definition JRange.hh:202
JTRIGGER::JMatchL0
L0 match criterion.
Definition JMatchL0.hh:29
JTRIGGER::JSuperFrame1D
1-dimensional frame with time calibrated data from one optical module.
Definition JSuperFrame1D.hh:40
JTRIGGER::JSuperFrame2D
2-dimensional frame with time calibrated data from one optical module.
Definition JSuperFrame2D.hh:46
GAUSS_LEGENDRE::f1
const JPolynome f1(1.0, 2.0, 3.0)
Function.
JDETECTOR::getNumberOfFloors
int getNumberOfFloors(const JDetector &detector)
Get number of floors.
Definition JDetectorToolkit.hh:243
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition JDetectorToolkit.hh:508
JDETECTOR::getNumberOfStrings
static const JStringCounter getNumberOfStrings
Function object to count unique strings.
Definition JStringCounter.hh:95
JLANG::getCount
size_t getCount(const array_type< T > &buffer, const JCompare_t &compare)
Count number of unique values.
Definition JVectorize.hh:262
JMATH::factorial
long long int factorial(const long long int n)
Determine factorial.
Definition JMathToolkit.hh:43
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:43
JROOT::getROOTClassSelection
std::set< JROOTClassSelector > getROOTClassSelection(const bool option=false)
Get ROOT class selection.
Definition JROOTClassSelector.hh:318
JROOT::counter_type
Long64_t counter_type
Type definition for counter.
Definition JROOT/JCounter.hh:24
KM3NETDAQ
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
KM3NETDAQ::getFrameTime
double getFrameTime()
Get frame time duration.
Definition JDAQClock.hh:162
KM3NETDAQ::getTimeOfRTS
double getTimeOfRTS(const JDAQChronometer &chronometer)
Get time of last RTS in ns since start of run for a given chronometer.
Definition JDAQChronometer.hh:223
FIXED
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
JAANET::JRange_t
Type definition of range.
Definition JHead.hh:43
JAANET::detector
Detector file.
Definition JHead.hh:227
JLANG::JType
Auxiliary class for a type holder.
Definition JType.hh:19
JROOT::JROOTClassSelector
Auxiliary class to select ROOT class based on class name.
Definition JROOTClassSelector.hh:38
JSUPPORT::JAutoTreeScanner
Auxiliary class to select JTreeScanner based on ROOT class name.
Definition JAutoTreeScanner.hh:36
JSUPPORT::JLimit
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45
JSUPPORT::JLimit::max
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
JTOOLS::JQuantile
Auxiliary data structure for running average, standard deviation and quantiles.
Definition JQuantile.hh:46
JTOOLS::JQuantile::put
void put(const double x, const double w=1.0)
Put value.
Definition JQuantile.hh:133
JTOOLS::JQuantile::getMean
double getMean() const
Get mean value.
Definition JQuantile.hh:252
Generated by doxygen 1.12.0