JCompareDetector.cc File Reference

#include <string>
#include <iostream>
#include <iomanip>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "JGeometry3D/JQuaternion3D.hh"
#include "JTools/JRange.hh"
#include "JTools/JConstants.hh"
#include "JTools/JQuantile.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 outout file with histograms is optionally produced.

Author

mjongen

Definition in file JCompareDetector.cc.

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 56 of file JCompareDetector.cc.

{
  using namespace std;
  using namespace JPP;
 
  string  detectorFile_a;
  string  detectorFile_b;
  string  outputFile;
  double  precision;
  int     debug;
 
  try {
 
    JParser<> zap("Auxiliary program to find differences between two detector files.");
 
    zap['a'] = make_field(detectorFile_a);
    zap['b'] = make_field(detectorFile_b);
    zap['o'] = make_field(outputFile)            = "";
    zap['p'] = make_field(precision)             = 0.001; // [ns,m]
    zap['d'] = make_field(debug)                 = 3;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
  
  JDetector detector_a;
  JDetector detector_b;
 
  try {
    load(detectorFile_a, detector_a);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
  try {
    load(detectorFile_b, detector_b);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
 
  bool is_equal = true;
 
  const JModuleRouter module_router_a(detector_a);
  const JModuleRouter module_router_b(detector_b);
 
 
  // compare detector IDs
 
  if (detector_a.getID() != detector_b.getID()) {
 
    DEBUG("* Different detector identifiers " 
          << setw(5) << detector_a.getID() << " (A) and " << endl 
          << setw(5) << detector_b.getID() << " (B)."     << endl
          << endl);
 
    is_equal = false;
  }
 
 
  // check whether the same modules are present in the file
 
  for (JDetector::iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
 
    if (!module_router_b.hasModule(module->getID())) {
 
      DEBUG("* Module " << setw(10) << module->getID() << " is in A " << getModuleLabel(*module) << " but not in B" << endl);
 
      is_equal = false;
    }
  }
 
  for (JDetector::iterator module = detector_b.begin(); module != detector_b.end(); ++module) {
 
    if (!module_router_a.hasModule(module->getID())) {
 
      DEBUG("* Module " << setw(10) << module->getID() << " is in B " << getModuleLabel(*module) << " but not in A" << endl);
 
      is_equal = false;
    }
  }
  DEBUG(endl);
 
 
  // compare the modules that the files have in common
 
  DEBUG("Comparing module by module." << endl);
 
  for (JDetector::iterator module_a = detector_a.begin(); module_a != detector_a.end(); ++module_a) {
 
    if (!module_router_b.hasModule(module_a->getID())) {
 
      continue;
 
      is_equal = false;
    }
    
    const JModule* module_b = &module_router_b.getModule(module_a->getID());
    
    DEBUG("  Module " << setw(10) << module_a->getID());
           
    // string and floor numbers
 
    if (module_a->getLocation() == module_b->getLocation()) {
 
      DEBUG(" " << getModuleLabel(*module_a) << endl);
 
    } else {
 
      DEBUG(endl);
      DEBUG("  * different location: "
            << getModuleLabel(*module_a) << " (A), "
            << getModuleLabel(*module_b) << " (B)" << endl);
 
      is_equal = false;
    }
 
    // average DOM positions
 
    if (module_a->getPosition().getDistance(module_b->getPosition()) > precision) {
 
      DEBUG("  * different position: "
            << module_a->getPosition() << " (A), " 
            << module_b->getPosition() << " (B)" 
            << ", B - A " << JPosition3D(module_b->getPosition()-module_a->getPosition()) << endl);
 
      is_equal = false;
    }
 
    // number of PMTs
 
    if (module_a->size() != module_b->size()) {
 
      DEBUG("  * different number of PMTs: "
            << module_a->size() << " (A), " 
            << module_b->size() << " (B)" << endl);
 
      is_equal = false;
    }
           
    // average t0
 
    const JQuantile q = compareT0(*module_a, *module_b);
 
    if (fabs(q.getMean()) > precision) {
 
      DEBUG("  * different average t0: "
            << ", B - A " << q.getMean()
            << endl);
      
      is_equal = false;
    }
 
    // comparing by PMT
 
    // t0
 
    for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
 
      const JPMT& pmt_a = module_a->getPMT(pmt);
      const JPMT& pmt_b = module_b->getPMT(pmt);
 
      if (fabs(pmt_a.getT0() - pmt_b.getT0()) > precision) {
 
        DEBUG("  * different T0 PMT " << pmt << ": "
              << pmt_a.getT0() << " (A" << FILL(2,'0') << pmt << "), " 
              << pmt_b.getT0() << " (B" << FILL(2,'0') << pmt << ")" 
              << ", B - A " << pmt_b.getT0() - pmt_a.getT0()  
              << endl);
        
        is_equal = false;
      }
    }
 
    // positions
 
    for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
 
      const JPMT& pmt_a = module_a->getPMT(pmt);
      const JPMT& pmt_b = module_b->getPMT(pmt);
 
      // PMT positions
 
      if (pmt_a.getPosition().getDistance(pmt_b.getPosition()) > precision) {
 
        DEBUG("  * different PMT position: "
              << pmt_a.getPosition() << " (A" << FILL(2,'0') << pmt << "), " 
              << pmt_b.getPosition() << " (B" << FILL(2,'0') << pmt << ")" 
              << ", B - A " << JPosition3D(pmt_b.getPosition() - pmt_a.getPosition()) << endl);
        
        is_equal = false;
      }
    }
 
    // status
 
    for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
 
      const JPMT& pmt_a = module_a->getPMT(pmt);
      const JPMT& pmt_b = module_b->getPMT(pmt);
 
      if (pmt_a.getStatus() != pmt_b.getStatus()) {
 
        DEBUG("  * different status PMT " << pmt << ": "
              << pmt_a.getStatus() << " (A" << FILL(2,'0') << pmt << "), " 
              << pmt_b.getStatus() << " (B" << FILL(2,'0') << pmt << ")"  
              << endl);
        
        is_equal = false;
      }
    }
  }    
  DEBUG(endl);
 
 
  if (outputFile != "") {
 
    typedef JRange<int>   JRange_t;
    
    const JRange_t string = combine(JRange_t(detector_a.begin(), detector_a.end(), &JModule::getString),
                                    JRange_t(detector_b.begin(), detector_b.end(), &JModule::getString));
    const JRange_t floor  = combine(JRange_t(detector_a.begin(), detector_a.end(), &JModule::getFloor),
                                    JRange_t(detector_b.begin(), detector_b.end(), &JModule::getFloor));
 
    TFile out(outputFile.c_str(), "recreate");    
    
    TH2D M2("M2", NULL,
            string.getLength() + 1,
            string.getLowerLimit() - 0.5,
            string.getUpperLimit() + 0.5,
            floor.getLength() + 1,
            floor.getLowerLimit() - 0.5,
            floor.getUpperLimit() + 0.5);
 
    TH2D* X2  = (TH2D*) M2.Clone("X2" );
    TH2D* Y2  = (TH2D*) M2.Clone("Y2" );
    TH2D* Z2  = (TH2D*) M2.Clone("Z2" );
    TH2D* T2  = (TH2D*) M2.Clone("T2" );
    TH2D* RM2 = (TH2D*) M2.Clone("RM2");
    TH2D* R2  = (TH2D*) M2.Clone("R2" );
    TH2D* QA  = (TH2D*) M2.Clone("QA" );
    TH2D* QB  = (TH2D*) M2.Clone("QB" );
    TH2D* QC  = (TH2D*) M2.Clone("QC" );
    TH2D* QD  = (TH2D*) M2.Clone("QD" );
 
    for( JDetector::iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
      if( !module_router_b.hasModule(module->getID()) ) {
        M2.Fill(module->getString(), module->getFloor(), -1.0);
      }
    }
    
    for( JDetector::iterator module = detector_b.begin(); module != detector_b.end(); ++module) {
      if( !module_router_a.hasModule(module->getID()) ) {
        M2.Fill(module->getString(), module->getFloor(), +1.0);
      }
    }
 
    for( JDetector::iterator module_a = detector_a.begin(); module_a != detector_a.end(); ++module_a) {
 
      if (!module_router_b.hasModule(module_a->getID())) {
        continue;
      }
      
      const JModule* module_b = &module_router_b.getModule(module_a->getID());
 
      double phi = 0.0;
 
      for (unsigned int pmt = 0; pmt != module_a->size() && pmt != module_b->size(); ++pmt) {
 
        double x = module_a->getPMT(pmt).getPhi()  -  module_b->getPMT(pmt).getPhi();
 
        while (x > +PI) { x -= PI; }
        while (x < -PI) { x += PI; }
      
        phi += x; 
      }
 
      phi /= min(module_a->size(),
                 module_b->size());
 
      JQuaternion3D Q = getRotation(*module_a, *module_b);
 
      JQuantile q = compareT0(*module_a, *module_b);
 
      X2 ->Fill(module_a->getString(), module_a->getFloor(), module_a->getX() - module_b->getX());
      Y2 ->Fill(module_a->getString(), module_a->getFloor(), module_a->getY() - module_b->getY());
      Z2 ->Fill(module_a->getString(), module_a->getFloor(), module_a->getZ() - module_b->getZ());
      T2 ->Fill(module_a->getString(), module_a->getFloor(), q.getMean());
      RM2->Fill(module_a->getString(), module_a->getFloor(), q.getRMS());
      R2 ->Fill(module_a->getString(), module_a->getFloor(), phi * 180.0/PI);
      QA ->Fill(module_a->getString(), module_a->getFloor(), Q.getA());
      QB ->Fill(module_a->getString(), module_a->getFloor(), Q.getB());
      QC ->Fill(module_a->getString(), module_a->getFloor(), Q.getC());
      QD ->Fill(module_a->getString(), module_a->getFloor(), Q.getD());
    }
 
    out.Write();
    out.Close();
  }
 
  ASSERT(is_equal);
 
  return 0;
}