mergeTimeOffsets.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"

#include "JDAQ/JDAQTimesliceIO.hh"
#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 "JTrigger/JTriggerToolkit.hh"

#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


namespace {

  using JDETECTOR::JModule;
  using namespace JPP ;   
  /**
   * Compare the t0s of two optical modules.
   *
   * \param  module_a        module
   * \param  module_b        module
   * \return q       
   */
  inline JQuantile compareT0(const JModule& module_a, const JModule& module_b)
  {
    JQuantile q ;
    
    JModule::const_iterator pmt_a = module_a.begin();
    JModule::const_iterator pmt_b = module_b.begin();

    for ( ; pmt_a != module_a.end() && pmt_b != module_b.end(); ++pmt_a , ++pmt_b) {
      q.put (pmt_a->getT0()-pmt_b->getT0());
    }      
    return q;
  }
}


/**
 * \file
 *
 * Auxiliary program to merge the time offsets of two detector files.\n
 * The time offsets can be averaged (-T "average") or merge accordingly to the odd/even 
 * PMT identification (-T "odd"). If the keyword <tt>odd</tt> is used, the time offsets 
 * of the PMTs with an odd identity are extracted from the detector file corresponding 
 * to the -a option. The other time offsets are extracted from the detector file 
 * corresponding to the -b option. The parity is calculated as following: 
 * <pre>
 *      (floor id - 1)*31 + PMT id
 * </pre>
 * \author acreusot
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  string  detectorFile_a;
  string  detectorFile_b;
  string  outputFile;
  string  mergeType;
  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['T'] = make_field(mergeType)             = "average";
    zap['d'] = make_field(debug)                 = 1;
    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);

  JDetector detectorMerged;
  try {
    load(outputFile, detectorMerged);
  }
  catch(const JException& error) {
    FATAL(error);
  }
  const JModuleRouter module_routerMerged(detectorMerged);

  // 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;
  }

  // compare the modules that the files have in common
  DEBUG("Comparing module by module." << endl);
  for (JDetector::iterator module = detectorMerged.begin(); module != detectorMerged.end(); ++module) {
    if (!module_router_b.hasModule(module->getID())) {
      continue;
      is_equal = false;
    }
    const JModule* module_b = &module_router_b.getModule(module->getID());
    if (!module_router_a.hasModule(module->getID())) {
      continue;
      is_equal = false;
    }
    const JModule* module_a = &module_router_a.getModule(module->getID());
    DEBUG("  Module " << setw(10) << module->getID());

    // comparing by PMT
    // t0
    for (unsigned int pmt = 0; pmt < module->size(); ++pmt) {
      const JPMT& pmt_a = module_a->getPMT(pmt);
      const JPMT& pmt_b = module_b->getPMT(pmt);
      double newT0 = 0;
      if (mergeType == "average") {
        newT0 = 0.5*(pmt_a.getT0() + pmt_b.getT0());
      } else if (mergeType == "odd") {
        const int floorId = module_a->getFloor();
        const int parityId = (floorId - 1)*NUMBER_OF_PMTS + pmt;
        if (parityId % 2 != 0) {
          newT0 = pmt_a.getT0();
        } else {
          newT0 = pmt_b.getT0();
        }
      }
      module->getPMT(pmt).setT0(newT0);
    }
  }    
  DEBUG(endl);
  //NOTICE("Store calibration data on file " << outputFile << endl);
  //detectorMerged.comment.add(JMeta(argc, argv));
  store(outputFile, detectorMerged);
  ASSERT(is_equal);
  return 0;
}