Jpp 19.3.0-rc.1
the software that should make you happy
Loading...
Searching...
No Matches
JOMGsim.cc
Go to the documentation of this file.
1#include <string>
2#include <iostream>
3#include <iomanip>
4#include <limits>
5
6#include "TROOT.h"
7#include "TFile.h"
8#include "TH1D.h"
9#include "TH2D.h"
10#include "TRandom3.h"
11
14
16
21
24#include "JSupport/JSupport.hh"
25#include "JSupport/JMeta.hh"
26
28#include "JTools/JRouter.hh"
29
30#include "JAAnet/JHead.hh"
31
33
34#include "JROOT/JRootToolkit.hh"
36
37#include "Jeep/JPrint.hh"
38#include "Jeep/JParser.hh"
39#include "Jeep/JMessage.hh"
40
41
42/**
43 * \file
44 *
45 * Auxiliary program to determine PMT parameters from OMGSim data.
46 *
47 * \author mdejong
48 */
49int main(int argc, char **argv)
50{
51 using namespace std;
52 using namespace JPP;
53 using namespace KM3NETDAQ;
54
56 JLimit_t& numberOfEvents = inputFile.getLimit();
57 string outputFile;
58 string detectorFile;
59 Double_t TMax_ns;
60 double R;
61 JPMTParameters parameters;
62 ULong_t seed;
63 int debug;
64
65 parameters.TTS_ns = 0.0;
66 parameters.gainSpread = 0.0;
67 parameters.PunderAmplified = 0.0;
68 parameters.slewing = false;
69
70 try {
71
72 JProperties properties = parameters.getProperties();
73
74 JParser<> zap("Auxiliary program to determine PMT parameters from OMGsim data.");
75
76 zap['f'] = make_field(inputFile, "input file.");
77 zap['o'] = make_field(outputFile, "output file.") = "omgsim.root";
78 zap['n'] = make_field(numberOfEvents) = JLimit::max();
79 zap['a'] = make_field(detectorFile, "detector file.");
80 zap['T'] = make_field(TMax_ns, "time window [ns].") = 20.0;
81 zap['R'] = make_field(R, "radioactivity [Bq]");
82 zap['P'] = make_field(properties, "PMT parameters") = JPARSER::initialised();
83 zap['S'] = make_field(seed, "seed") = 0;
84 zap['d'] = make_field(debug, "debug flag.") = 1;
85
86 zap(argc, argv);
87 }
88 catch(const exception &error) {
89 FATAL(error.what() << endl);
90 }
91
92
93 gRandom->SetSeed(seed);
94
95
97
98 try {
99 load(detectorFile, detector);
100 }
101 catch(const JException& error) {
102 FATAL(error);
103 }
104
105 if (detector.size() != 1u) {
106 FATAL("Wrong number of modules " << detector.size() << endl);
107 }
108
109
110 const JModule& module = detector[0];
111
112 JCombinatorics combinatorics(module.size());
113
114 combinatorics.sort(JPairwiseComparator(module));
115
116 const Int_t nx = combinatorics.getNumberOfPairs();
117 const Double_t xmin = -0.5;
118 const Double_t xmax = nx - 0.5;
119
120 const Double_t ymin = -floor(TMax_ns) + 0.125;
121 const Double_t ymax = +floor(TMax_ns) - 0.125;
122 const Int_t ny = (Int_t) ((ymax - ymin) / 0.25);
123
124 TH2D h2s(MAKE_CSTRING(module.getID() << _2S), NULL, nx, xmin, xmax, ny, ymin, ymax);
125
126 h2s.Sumw2();
127
128
129 JRouter<int> router;
130
131 for (size_t i = 0; i != module.size(); ++i) {
132 router.put(module[i].getID(), i);
133 }
134
135 double numberOfPrimaries = 0.0;
136
137 Head header;
138
139 try {
140
141 STATUS("Extracting header data... " << flush);
142
143 header = getHeader(inputFile);
144
145 const JHead buffer(header);
146
147 if (buffer.is_valid(&JHead::norma))
148 numberOfPrimaries = buffer.norma.numberOfPrimaries;
149 else
150 THROW(JException, "Missing header information.");
151
152 STATUS("OK" << endl);
153
154 } catch(const exception& error) {
155 FATAL(error.what() << endl);
156 }
157
158 if (numberOfPrimaries == 0.0) {
159 FATAL("Number of primaries " << numberOfPrimaries << endl);
160 }
161
162
163 const JPMTAnalogueSignalProcessor cpu(parameters);
164
166 JPMTData<JPMTPulse> output;
167
168 struct hit_type {
169 int pmt;
170 double t1;
171 };
172
173 vector<hit_type> buffer;
174
175
176 TFile out(outputFile.c_str(), "recreate");
177
178 putObject(out, JMeta(argc, argv));
179
180 while (inputFile.hasNext()) {
181
182 if (inputFile.getCounter()%10000== 0) {
183 STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
184 }
185
186 const Evt* evt = inputFile.next();
187
188 if (evt->mc_hits.size() >= 2u) {
189
190 buffer.clear();
191
192 for (const auto& hit : evt->mc_hits) {
193
194 if (router.has(hit.pmt_id)) {
195
196 const int pmt = router.get(hit.pmt_id);
197
198 input .clear();
199 output.clear();
200
201 input.push_back(JPMTSignal(hit.t, hit.a));
202
203 cpu(module[pmt].getCalibration(), input, output);
204
205 for (const auto& hit : output) {
206 buffer.push_back({pmt, hit.t_ns});
207 }
208 }
209 }
210
211
212 for (vector<hit_type>::const_iterator p = buffer.begin(); p != buffer.end(); ++p) {
213 for (vector<hit_type>::const_iterator q = buffer.begin(); p != q; ++q) {
214 h2s.Fill((double) combinatorics.getIndex(p->pmt,q->pmt),
215 JCombinatorics::getSign(p->pmt,q->pmt) * (q->t1 - p->t1));
216 }
217 }
218 }
219 }
220 STATUS(endl);
221
222
223 const double W = R / numberOfPrimaries;
224
225 h2s.Scale(W);
226
227 out << h2s;
228
229 out.Write();
230 out.Close();
231}
string outputFile
KM3NeT DAQ constants, bit handling, etc.
Data structure for detector geometry and calibration.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define STATUS(A)
Definition JMessage.hh:63
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
ROOT I/O of application specific meta data.
Scanning of objects from multiple files according a format that follows from the extension of each fi...
int main(int argc, char **argv)
Definition JOMGsim.cc:49
PMT analogue signal processor.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
I/O formatting auxiliaries.
#define MAKE_CSTRING(A)
Make C-string.
Definition JPrint.hh:72
ROOT TTree parameter settings of various packages.
Monte Carlo run header.
Definition JHead.hh:1236
JAANET::norma norma
Definition JHead.hh:1603
bool is_valid(T JHead::*pd) const
Check validity of given data member in JHead.
Definition JHead.hh:1319
Detector data structure.
Definition JDetector.hh:96
Data structure for a composite optical module.
Definition JModule.hh:75
Template data structure for PMT I/O.
Data structure for PMT parameters.
double gainSpread
gain spread [unit]
double TTS_ns
transition time spread [ns]
JProperties getProperties(const JEquationParameters &equation=JPMTParameters::getEquationParameters())
Get properties of this class.
double PunderAmplified
probability of underamplified hit
bool slewing
time slewing of analogue signal
Utility class to parse parameter values.
General exception.
Definition JException.hh:24
int getID() const
Get identifier.
Definition JObjectID.hh:50
Utility class to parse command line options.
Definition JParser.hh:1698
General purpose class for object reading from a list of file names.
virtual bool hasNext() override
Check availability of next element.
counter_type getCounter() const
Get counter.
virtual const pointer_type & next() override
Get next element.
Auxiliary class to convert pair of indices to unique index and back.
int getIndex(const int first, const int second) const
Get index of pair of indices.
static int getSign(const int first, const int second)
Sign of pair of indices.
void sort(JComparator_t comparator)
Sort address pairs.
size_t getNumberOfPairs() const
Get number of pairs.
Direct addressing of elements with unique identifiers.
Definition JRouter.hh:27
static const char *const _2S
Name extension for 2D counts.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
JCalibration getCalibration(const JCalibration &first, const JCalibration &second)
Get calibration to go from first to second calibration.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
bool putObject(TDirectory &dir, const TObject &object)
Write object to ROOT directory.
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition Evt.hh:21
std::vector< Hit > mc_hits
MC: list of MC truth hits.
Definition Evt.hh:48
The Head class reflects the header of Monte-Carlo event files, which consists of keys (also referred ...
Definition Head.hh:65
Detector file.
Definition JHead.hh:227
double numberOfPrimaries
Number of primaries.
Definition JHead.hh:826
Acoustic hit.
Definition JBillabong.cc:70
Auxiliary class to sort pairs of PMT addresses within optical module.
Data structure for PMT analogue signal.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
Auxiliary class for ROOT I/O of application specific meta data.
Definition JMeta.hh:72