Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JTriggerTester.cc
Go to the documentation of this file.
1 
2 #include <string>
3 #include <iostream>
4 #include <iomanip>
5 #include <utility>
6 
10 
11 #include "JAAnet/JHead.hh"
12 #include "JAAnet/JHeadToolkit.hh"
13 #include "JAAnet/JAAnetToolkit.hh"
14 
16 #include "JDetector/JDetector.hh"
18 #include "JDetector/JPMTRouter.hh"
20 #include "JDetector/JK40Rates.hh"
23 
25 
28 #include "JSupport/JSupport.hh"
29 
30 #include "JTools/JHashMap.hh"
31 #include "JTools/JWeight.hh"
32 #include "JTools/JQuantile.hh"
33 
34 #include "Jeep/JPrint.hh"
35 #include "Jeep/JParser.hh"
36 #include "Jeep/JMessage.hh"
37 
38 namespace {
39 
40  using namespace JPP;
41 
42  /**
43  * Auxiliary data structure for monitoring number of photo-electrons.
44  */
45  struct tuple_type :
46  public std::vector<JWeight>
47  {
48  /**
49  * Default constructor.
50  */
51  tuple_type() :
52  std::vector<JWeight> (3, JWeight())
53  {}
54  };
55 
56 
57  /**
58  * Hash for PMT identifier
59  */
60  struct JHashEvaluator_t {
61  /**
62  * Get hash value.
63  *
64  * \param pmt PMT identifier
65  * \return hash value
66  */
67  inline int operator()(const JPMTIdentifier& pmt) const
68  {
69  return pmt.getID() * KM3NETDAQ::NUMBER_OF_PMTS + pmt.getPMTAddress();
70  }
71  };
72 }
73 
74 
75 /**
76  * \file
77  * Auxiliary program to verify processing of Monte Carlo events.
78  *
79  * Note that the available options are identical to those of JTriggerEfficiency.cc.
80  *
81  * \author mdejong
82  */
83 int main(int argc, char **argv)
84 {
85  using namespace std;
86  using namespace JPP;
87 
88  JMultipleFileScanner<Evt> inputFile;
89  JLimit_t& numberOfEvents = inputFile.getLimit();
90  string detectorFileA;
91  string detectorFileB;
93  JPMTParametersMap pmtParameters;
94  JK40Rates rates_Hz;
95  int debug;
96 
97  try {
98 
99  JParser<> zap("Auxiliary program to verify processing of Monte Carlo events.");
100 
101  zap['f'] = make_field(inputFile, "input file (output of detector simulation)");
102  zap['n'] = make_field(numberOfEvents) = JLimit::max();
103  zap['a'] = make_field(detectorFileA, "detector used for converion from Monte Carlo truth to raw data.");
104  zap['b'] = make_field(detectorFileB, "detector used for conversion of raw data to calibrated data.") = "";
105  zap['@'] = make_field(parameters, "Trigger parameters (or corresponding file name)") = JPARSER::initialised();
106  zap['P'] = make_field(pmtParameters, "PMT simulation data (or corresponding file name)") = JPARSER::initialised();
107  zap['B'] = make_field(rates_Hz, "background rates [Hz]") = JPARSER::initialised();
108  zap['d'] = make_field(debug, "debug") = 2;
109 
110  zap(argc, argv);
111  }
112  catch(const exception &error) {
113  FATAL(error.what() << endl);
114  }
115 
116  cout.tie(&cerr);
117 
118  if (detectorFileB == "") {
119  detectorFileB = detectorFileA;
120  }
121 
122 
123  JDetector detectorA;
124  JDetector detectorB;
125 
126  try {
127  load(detectorFileA, detectorA);
128  load(detectorFileB, detectorB);
129  }
130  catch(const JException& error) {
131  FATAL(error);
132  }
133 
134  NOTICE("Number of modules in detector A <" << detectorFileA << ">: " << setw(4) << detectorA.size() << endl);
135  NOTICE("Number of modules in detector B <" << detectorFileB << ">: " << setw(4) << detectorB.size() << endl);
136 
137  JPMTParametersMap::Throw(true);
138 
139  if (!pmtParameters.is_valid()) {
140  FATAL("Invalid PMT parameters " << pmtParameters << endl);
141  }
142 
143  if (pmtParameters.getQE() != 1.0) {
144 
145  NOTICE("Correct background rates with global efficiency " << pmtParameters.getQE() << endl);
146 
147  rates_Hz.correct(pmtParameters.getQE());
148 
149  NOTICE("Back ground rates: " << rates_Hz << endl);
150  }
151 
152  const JPMTRouter pmtRouter (detectorA);
153  const JModuleRouter moduleRouter(detectorB);
154 
155  parameters.set(getMaximalDistance(detectorB));
156 
157  DEBUG("Trigger:" << endl << parameters << endl);
158  DEBUG("PMT parameters:" << endl << pmtParameters << endl);
159 
160  Head header;
161 
162  try {
163  header = getHeader(inputFile);
164  }
165  catch(const JException& error) {
166  FATAL(error);
167  }
168 
169  const JPosition3D center = get<JPosition3D>(header);
170 
171  NOTICE("Apply detector offset from Monte Carlo run header (" << center << ")" << endl);
172 
173  detectorA -= center;
174  detectorB -= center;
175 
180 
181  JQuantile QT;
182 
183  while (inputFile.hasNext()) {
184 
185  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r');
186 
187  Evt* event = inputFile.next();
188 
189  if (!event->mc_hits.empty()) {
190  QT.put(getTimeRange(*event).getLength());
191  }
192 
193  for (vector<Hit>::const_iterator hit = event->mc_hits.begin(); hit != event->mc_hits.end(); ++hit) {
194 
195  if (!pmtRouter.hasPMT(hit->pmt_id)) {
196 
197  miss[hit->pmt_id].put(hit->pure_a);
198 
199  continue;
200  }
201 
202  const JPMTIdentifier pmt = pmtRouter.getIdentifier(hit->pmt_id);
203 
204  if (!moduleRouter.hasModule(pmt.getID())) {
205 
206  lost[pmt].put(hit->pure_a);
207 
208  continue;
209  }
210 
211  const double t0 = hit->pure_t;
212  const double t1 = putTime(t0, pmtRouter .getPMT(hit->pmt_id));
213  const double t2 = getTime(t1, moduleRouter.getPMT(pmt));
214 
215  if (fabs(t2 - t0) > parameters.TMaxLocal_ns) {
216  slip[pmt].put(hit->pure_a);
217  }
218 
219  const JPMTParameters data = pmtParameters.getPMTParameters(pmt);
220  const double P = getSurvivalProbability(data);
221 
222  npe[pmt][0].put(hit->pure_a);
223  npe[pmt][1].put(hit->pure_a * P);
224  npe[pmt][2].put(hit->pure_a * P * data.QE);
225  }
226  }
227  STATUS(endl);
228 
229 
230  NOTICE("Number of PMTs absent in detector A: " << setw(6) << miss.size() << endl);
231 
232  for (JHashMap<int, JWeight>::const_iterator i = miss.begin(); i != miss.end(); ++i) {
233  DEBUG(setw(5) << i->first << ' ' << FIXED(8,0) << i->second.getTotal() << endl);
234  }
235 
236  NOTICE("Number of PMTs absent in detector B: " << setw(6) << lost.size() << endl);
237 
238  for (JHashMap<JPMTIdentifier, JWeight>::const_iterator i = lost.begin(); i != lost.end(); ++i) {
239  DEBUG(setw(8) << i->first.getModuleID() << "[" << setw(2) << i->first.getPMTAddress() << "] " << FIXED(8,0) << i->second.getTotal() << endl);
240  }
241 
242  NOTICE("Number of PMTs with t0 detector A - B > " << FIXED(4,1) << parameters.TMaxLocal_ns << " [ns]: " << setw(6) << slip.size() << endl);
243 
244  for (JHashMap<JPMTIdentifier, JWeight>::const_iterator i = slip.begin(); i != slip.end(); ++i) {
245  DEBUG(setw(8) << i->first.getModuleID() << "[" << setw(2) << i->first.getPMTAddress() << "] " << FIXED(8,0) << i->second.getTotal() << endl);
246  }
247 
248 
249  NOTICE("Number of true photo-electrons, passed threshold and survived QE." << endl);
250 
251  tuple_type total;
252 
253  for (JHashMap<JPMTIdentifier, tuple_type>::const_iterator p = npe.begin(); p != npe.end(); ++p) {
254 
255  DEBUG(setw(8) << p->first.getModuleID() << "[" << setw(2) << p->first.getPMTAddress() << "]");
256 
257  for (size_t i = 0; i != p->second.size(); ++i) {
258  DEBUG(' ' << FIXED(8,0) << p->second[i].getTotal());
259  }
260  DEBUG(endl);
261 
262  for (size_t i = 0; i != p->second.size(); ++i) {
263  total[i].put(p->second[i].getTotal());
264  }
265  }
266 
267  NOTICE(setw(12) << "total");
268 
269  for (size_t i = 0; i != total.size(); ++i) {
270  NOTICE(' ' << FIXED(8,0) << total[i].getTotal());
271  }
272 
273  NOTICE(endl);
274 
275  NOTICE("Time range of hits [ns]: " << QT.getMin() << " - " << QT.getMax() << endl);
276 }
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
ROOT TTree parameter settings.
General purpose class for hash map of unique keys.
Definition: JHashMap.hh:71
JPMTIdentifier getIdentifier(const JPMTAddress &address) const
Get identifier of PMT.
Definition: JPMTRouter.hh:126
Quantile calculator.
Definition: JQuantile.hh:83
General purpose class for hash map of unique elements.
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:80
Router for direct addressing of module data in detector data structure.
const JPMTParameters & getPMTParameters() const
Get PMT parameters.
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
Definition: diff-Tuna.sh:38
bool hasPMT(const JObjectID &id) const
Has PMT.
Definition: JPMTRouter.hh:114
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:63
Auxiliary data structure for floating point format specification.
Definition: JPrint.hh:481
double getTime(const Hit &hit)
Get true time of hit.
Data structure for detector geometry and calibration.
esac $JPP_DIR examples JDetector JTransitTime o $OUTPUT_FILE n N $NPE T $TTS_NS d $DEBUG for HISTOGRAM in tts tt2 pmt
Definition: JTransitTime.sh:36
JTimeRange getTimeRange(const Evt &event)
Get time range (i.e. time between earliest and latest hit) of Monte Carlo event.
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
double putTime(const T &t1, const JCalibration &cal)
Get de-calibrated time.
double getMaximalDistance(const JDetector &detector)
Get maximal distance between modules in detector.
I/O formatting auxiliaries.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1954
int getID() const
Get identifier.
Definition: JObjectID.hh:55
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
#define NOTICE(A)
Definition: JMessage.hh:64
int getPMTAddress() const
Get PMT identifier (= TDC).
Auxiliary class for map of PMT parameters.
Weight calculator.
Definition: JWeight.hh:23
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
General purpose messaging.
The Head class reflects the header of Monte-Carlo event files, which consists of keys (also referred ...
Definition: Head.hh:66
#define FATAL(A)
Definition: JMessage.hh:67
Scanning of objects from multiple files according a format that follows from the extension of each fi...
const JPMT & getPMT(const JPMTIdentifier &id) const
Get PMT parameters.
Direct access to module in detector data structure.
General purpose class for object reading from a list of file names.
Utility class to parse command line options.
double getSurvivalProbability(const JPMTParameters &parameters)
Get model dependent probability that a one photo-electron hit survives the simulation of the PMT assu...
bool hasModule(const JObjectID &id) const
Has module.
Data structure for position in three dimensions.
Definition: JPosition3D.hh:35
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
container_type::const_iterator const_iterator
Definition: JHashMap.hh:85
Data structure for PMT parameters.
KM3NeT DAQ constants, bit handling, etc.
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
double QE
relative quantum efficiency
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:19
then $DIR JPlotNPE PDG P
Definition: JPlotNPE-PDG.sh:60
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
Auxiliary class for K40 rates.
Definition: JK40Rates.hh:41
int main(int argc, char *argv[])
Definition: Main.cpp:15