Jpp  pmt_effective_area_update_2
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
JCompareDetector.cc File Reference

Auxiliary program to find differences between two detector files. More...

#include <string>
#include <iostream>
#include <iomanip>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "JGeometry3D/JQuaternion3D.hh"
#include "JMath/JConstants.hh"
#include "JTools/JRange.hh"
#include "JTools/JQuantile.hh"
#include "JROOT/JManager.hh"
#include "JROOT/JRootToolkit.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.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

Auxiliary program to find differences between two detector files.


A ROOT output file with histograms is optionally produced.

Author
mjongen

Definition in file JCompareDetector.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 71 of file JCompareDetector.cc.

72 {
73  using namespace std;
74  using namespace JPP;
75 
76  string detectorFile_a;
77  string detectorFile_b;
78  string outputFile;
79  double precision;
80  int debug;
81 
82  try {
83 
84  JParser<> zap("Auxiliary program to find differences between two detector files.");
85 
86  zap['a'] = make_field(detectorFile_a);
87  zap['b'] = make_field(detectorFile_b);
88  zap['o'] = make_field(outputFile) = "";
89  zap['p'] = make_field(precision) = 0.001; // [ns,m]
90  zap['d'] = make_field(debug) = 3;
91 
92  zap(argc, argv);
93  }
94  catch(const exception &error) {
95  FATAL(error.what() << endl);
96  }
97 
98 
99  JDetector detector_a;
100  JDetector detector_b;
101 
102  try {
103  load(detectorFile_a, detector_a);
104  }
105  catch(const JException& error) {
106  FATAL(error);
107  }
108 
109  try {
110  load(detectorFile_b, detector_b);
111  }
112  catch(const JException& error) {
113  FATAL(error);
114  }
115 
116  size_t numberOfPMTs = 0;
117 
118  bool is_equal = true;
119 
120  const JModuleRouter module_router_a(detector_a);
121  const JModuleRouter module_router_b(detector_b);
122 
123  setFormat<JPosition3D> (JFormat_t(15, 9, std::ios::fixed | std::ios::showpos));
124 
125  // compare detector identifiers
126 
127  if (detector_a.getID() != detector_b.getID()) {
128 
129  DEBUG("* Different detector identifiers "
130  << setw(5) << detector_a.getID() << " (A) and " << endl
131  << setw(5) << detector_b.getID() << " (B)." << endl
132  << endl);
133 
134  is_equal = false;
135  }
136 
137  for (JDetector::iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
138  if (module->size() > numberOfPMTs) {
139  numberOfPMTs = module->size();
140  }
141  }
142 
143  // check whether the same modules are present in the file
144 
145  for (JDetector::iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
146 
147  if (!module_router_b.hasModule(module->getID())) {
148 
149  DEBUG("* Module " << setw(10) << module->getID() << " is in A " << getLabel(*module) << " but not in B" << endl);
150 
151  is_equal = false;
152  }
153  }
154 
155  for (JDetector::iterator module = detector_b.begin(); module != detector_b.end(); ++module) {
156 
157  if (!module_router_a.hasModule(module->getID())) {
158 
159  DEBUG("* Module " << setw(10) << module->getID() << " is in B " << getLabel(*module) << " but not in A" << endl);
160 
161  is_equal = false;
162  }
163  }
164  DEBUG(endl);
165 
166 
167  // compare the modules that the files have in common
168 
169  DEBUG("Comparing module by module." << endl);
170 
171  for (JDetector::iterator module_a = detector_a.begin(); module_a != detector_a.end(); ++module_a) {
172 
173  if (!module_router_b.hasModule(module_a->getID())) {
174 
175  continue;
176 
177  is_equal = false;
178  }
179 
180  const JModule* module_b = &module_router_b.getModule(module_a->getID());
181 
182  DEBUG(" Module " << setw(10) << module_a->getID());
183 
184  // string and floor numbers
185 
186  if (module_a->getLocation() == module_b->getLocation()) {
187 
188  DEBUG(" " << getLabel(*module_a) << endl);
189 
190  } else {
191 
192  DEBUG(endl);
193  DEBUG(" * different location: "
194  << getLabel(*module_a) << " (A), "
195  << getLabel(*module_b) << " (B)" << endl);
196 
197  is_equal = false;
198  }
199 
200  // time offset
201 
202  if (fabs(module_a->getT0() - module_b->getT0()) > precision) {
203 
204  DEBUG(" * different T0: "
205  << module_a->getT0() << " (A), "
206  << module_b->getT0() << " (B) "
207  << ", B - A " << module_b->getT0() - module_a->getT0() << endl);
208 
209  is_equal = false;
210  }
211 
212  // quaternion calibration
213 
214  if (getAngle(module_a->getQuaternion(), module_b->getQuaternion()) > precision) {
215 
216  DEBUG(" * different quaternion calibration: "
217  << module_a->getQuaternion() << " (A), "
218  << module_b->getQuaternion() << " (B) "
219  << ", B - A " << getAngle(module_b->getQuaternion(), module_a->getQuaternion()) << endl);
220 
221  is_equal = false;
222  }
223 
224  // average DOM positions
225 
226  if (getDistance(module_a->getPosition(), module_b->getPosition()) > precision) {
227 
228  DEBUG(" * different position: "
229  << module_a->getPosition() << " (A), "
230  << module_b->getPosition() << " (B)"
231  << ", B - A " << JPosition3D(module_b->getPosition() - module_a->getPosition()) << endl);
232 
233  is_equal = false;
234  }
235 
236  // number of PMTs
237 
238  if (module_a->size() != module_b->size()) {
239 
240  DEBUG(" * different number of PMTs: "
241  << module_a->size() << " (A), "
242  << module_b->size() << " (B)" << endl);
243 
244  is_equal = false;
245  }
246 
247  // average t0
248 
249  if (!module_a->empty() &&
250  !module_b->empty()) {
251 
252  const JQuantile q = getQuantile(*module_a, *module_b);
253 
254  if (fabs(q.getMean()) > precision) {
255 
256  DEBUG(" * different average t0: "
257  << ", B - A " << q.getMean()
258  << endl);
259 
260  is_equal = false;
261  }
262  }
263 
264  // global t0
265 
266  if (fabs(module_a->getT0() - module_b->getT0()) > precision) {
267 
268  DEBUG(" * different global t0: "
269  << module_a->getT0() << " (A), "
270  << module_b->getT0() << " (B)"
271  << ", B - A " << module_b->getT0() - module_a->getT0()
272  << endl);
273 
274  is_equal = false;
275  }
276 
277  // comparing by PMT
278 
279  // t0
280 
281  for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
282 
283  const JPMT& pmt_a = module_a->getPMT(pmt);
284  const JPMT& pmt_b = module_b->getPMT(pmt);
285 
286  if (fabs(pmt_a.getT0() - pmt_b.getT0()) > precision) {
287 
288  DEBUG(" * different T0 PMT " << pmt << ": "
289  << pmt_a.getT0() << " (A" << FILL(2,'0') << pmt << "), " << FILL()
290  << pmt_b.getT0() << " (B" << FILL(2,'0') << pmt << ")" << FILL()
291  << ", B - A " << pmt_b.getT0() - pmt_a.getT0()
292  << endl);
293 
294  is_equal = false;
295  }
296  }
297 
298  // positions
299 
300  for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
301 
302  const JPMT& pmt_a = module_a->getPMT(pmt);
303  const JPMT& pmt_b = module_b->getPMT(pmt);
304 
305  // PMT positions
306 
307  if (pmt_a.getPosition().getDistance(pmt_b.getPosition()) > precision) {
308 
309  DEBUG(" * different PMT position: "
310  << pmt_a.getPosition() << " (A" << FILL(2,'0') << pmt << "), " << FILL()
311  << pmt_b.getPosition() << " (B" << FILL(2,'0') << pmt << ")" << FILL()
312  << ", B - A " << JPosition3D(pmt_b.getPosition() - pmt_a.getPosition()) << endl);
313 
314  is_equal = false;
315  }
316  }
317 
318  // status
319 
320  for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
321 
322  const JPMT& pmt_a = module_a->getPMT(pmt);
323  const JPMT& pmt_b = module_b->getPMT(pmt);
324 
325  if (pmt_a.getStatus() != pmt_b.getStatus()) {
326 
327  DEBUG(" * different status PMT " << pmt << ": "
328  << pmt_a.getStatus() << " (A" << FILL(2,'0') << pmt << "), " << FILL()
329  << pmt_b.getStatus() << " (B" << FILL(2,'0') << pmt << ")" << FILL()
330  << endl);
331 
332  is_equal = false;
333  }
334  }
335  }
336  DEBUG(endl);
337 
338 
339  if (outputFile != "") {
340 
341  set<int> string;
342  set<int> floor;
343 
344  for (JDetector::iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
345  string.insert(module->getString());
346  floor .insert(module->getFloor ());
347  }
348 
349  for (JDetector::iterator module = detector_b.begin(); module != detector_b.end(); ++module) {
350  string.insert(module->getString());
351  floor .insert(module->getFloor ());
352  }
353 
354  JManager<int, TH1D> H1(new TH1D("M[%]", NULL, numberOfPMTs, -0.5, numberOfPMTs - 0.5));
355 
356  TH2D M2("M2", NULL,
357  string.size(), -0.5, string.size() - 0.5,
358  floor .size(), -0.5, floor .size() - 0.5);
359 
360  for (set<int>::const_iterator i = string.begin(); i != string.end(); ++i) {
361  M2.GetXaxis()->SetBinLabel(distance(string.begin(), i) + 1, MAKE_CSTRING(*i));
362  }
363 
364  for (set<int>::const_iterator i = floor.begin(); i != floor.end(); ++i) {
365  M2.GetYaxis()->SetBinLabel(distance(floor.begin(), i) + 1, MAKE_CSTRING(*i));
366  }
367 
368  TH2D* X2 = (TH2D*) M2.Clone("X2" );
369  TH2D* Y2 = (TH2D*) M2.Clone("Y2" );
370  TH2D* Z2 = (TH2D*) M2.Clone("Z2" );
371  TH2D* T2 = (TH2D*) M2.Clone("T2" );
372  TH2D* RMS = (TH2D*) M2.Clone("RMS");
373  TH2D* R2 = (TH2D*) M2.Clone("R2" );
374  TH2D* QA = (TH2D*) M2.Clone("QA" );
375  TH2D* QB = (TH2D*) M2.Clone("QB" );
376  TH2D* QC = (TH2D*) M2.Clone("QC" );
377  TH2D* QD = (TH2D*) M2.Clone("QD" );
378  TH2D* Q2 = (TH2D*) M2.Clone("Q2" );
379 
380  for (JDetector::iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
381  if (!module_router_b.hasModule(module->getID()) ) {
382  M2.Fill(module->getString(), module->getFloor(), -1.0);
383  }
384  }
385 
386  for (JDetector::iterator module = detector_b.begin(); module != detector_b.end(); ++module) {
387  if (!module_router_a.hasModule(module->getID()) ) {
388  M2.Fill(module->getString(), module->getFloor(), +1.0);
389  }
390  }
391 
392 
393  for (JDetector::iterator module_a = detector_a.begin(); module_a != detector_a.end(); ++module_a) {
394 
395  if (!module_router_b.hasModule(module_a->getID())) {
396  continue;
397  }
398 
399  const JModule* module_b = &module_router_b.getModule(module_a->getID());
400 
401  for (size_t i = 0; i != numberOfPMTs; ++i) {
402 
403  if (module_a->getFloor() != 0) {
404  H1[module_a->getID()]->SetBinContent(i + 1, module_a->getPMT(i).getT0() - module_b->getPMT(i).getT0());
405  }
406  }
407 
408  X2 ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), module_a->getX() - module_b->getX());
409  Y2 ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), module_a->getY() - module_b->getY());
410  Z2 ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), module_a->getZ() - module_b->getZ());
411 
412  if (module_a->getFloor() != 0 &&
413  module_b->getFloor() != 0) {
414 
415  const JQuaternion3D Q = getRotation(*module_a, *module_b);
416  const JQuantile q = getQuantile(*module_a, *module_b);
417 
419 
420  const double phi = (JVector3Z_t.getDot(q1.twist) >= 0.0 ? +1.0 : -1.0) * q1.twist.getAngle();
421 
422  R2 ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), phi);
423  QA ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), Q.getA());
424  QB ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), Q.getB());
425  QC ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), Q.getC());
426  QD ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), Q.getD());
427  Q2 ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), Q.getAngle());
428  T2 ->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), q.getMean());
429  RMS->Fill(getBin(string, module_a->getString()), getBin(floor, module_a->getFloor()), q.getSTDev());
430  }
431  }
432 
433 
434  TFile out(outputFile.c_str(), "recreate");
435 
436  for (TH2D* h2 : { &M2, X2, Y2, Z2, T2, RMS, R2, QA, QB, QC, QD, Q2 }) {
437  out << *h2;
438  }
439 
440  out << H1;
441 
442  out.Write();
443  out.Close();
444  }
445 
446  ASSERT(is_equal);
447 
448  return 0;
449 }
Utility class to parse command line options.
Definition: JParser.hh:1500
double getAngle(const JQuaternion3D &first, const JQuaternion3D &second)
Get space angle between quanternions.
General exception.
Definition: JException.hh:23
Q(UTCMax_s-UTCMin_s)-livetime_s
const JStatus & getStatus() const
Get status.
Definition: JStatus.hh:67
double getB() const
Get b value.
int getFloor() const
Get floor number.
Definition: JLocation.hh:145
Data structure for a composite optical module.
Definition: JModule.hh:66
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
double getSTDev() const
Get standard deviation.
Definition: JQuantile.hh:266
Auxiliary data structure for running average, standard deviation and quantiles.
Definition: JQuantile.hh:43
do rm f tmp H1
std::string getLabel(const JLocation &location)
Get module label for monitoring and other applications.
Definition: JLocation.hh:246
Detector data structure.
Definition: JDetector.hh:81
Router for direct addressing of module data in detector data structure.
static JRotation getRotation
Function object to get rotation matrix to go from first to second module.
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:151
string outputFile
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition: JVector3D.hh:270
double getMean() const
Get mean value.
Definition: JQuantile.hh:252
double getDistance(const JFirst_t &first, const JSecond_t &second)
Get distance between objects.
const JQuaternion3D & getQuaternion() const
Get quaternion.
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
Definition: JManager.hh:43
#define ASSERT(A,...)
Assert macro.
Definition: JMessage.hh:90
const JLocation & getLocation() const
Get location.
Definition: JLocation.hh:69
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
Data structure for PMT geometry, calibration and status.
Definition: JPMT.hh:42
double getY() const
Get y position.
Definition: JVector3D.hh:104
int debug
debug level
Definition: JSirene.cc:63
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
const JPMT & getPMT(const int index) const
Get PMT.
Definition: JModule.hh:163
Auxiliary data structure for sequence of same character.
Definition: JManip.hh:328
#define FATAL(A)
Definition: JMessage.hh:67
double getD() const
Get d value.
Data structure for unit quaternion in three dimensions.
double getAngle() const
Get rotation angle.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
static const JVector3D JVector3Z_t(0, 0, 1)
unit z-vector
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
double getC() const
Get c value.
double getX() const
Get x position.
Definition: JVector3D.hh:94
double getA() const
Get a value.
double getDot(const JVector3D &vector) const
Get dot product.
Definition: JVector3D.hh:282
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
Auxiliary data structure for decomposition of quaternion in twist and swing quaternions.
Data structure for format specifications.
Definition: JManip.hh:522
double getZ() const
Get z position.
Definition: JVector3D.hh:115
double getT0() const
Get time offset.