JAlignDetector.cc File Reference

Auxiliary program to align two detectors. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <vector>
#include <set>
#include "JMath/JMath.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JTripod.hh"
#include "JFit/JSimplex.hh"
#include "JLang/JException.hh"
#include "JSupport/JMeta.hh"
#include "JTools/JRange.hh"
#include "Jeep/JContainer.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 align two detectors.

Author

mdejong

Definition in file JAlignDetector.cc.

Function Documentation

main()

int main	(	int	argc,
		char **	argv )

Definition at line 300 of file JAlignDetector.cc.

{
  using namespace std;
  using namespace JPP;
 
  const set<string>  PARAMETERS = { X_t, Y_t, Z_t, PHI_t, TX_t, TY_t };
 
  string             detectorFile_a;
  string             detectorFile_b;
  bool               overwriteDetector;
  string             tripodFile;
  double             sigma_m;
  bool               option;
  JRange_t           range;
  set<string>        parameters;
  int                debug;
 
  try {
 
    JParser<> zap("Auxiliary program to align two detectors.");
 
    zap['a'] = make_field(detectorFile_a,      "detector - subject to alignment (option -A)");
    zap['b'] = make_field(detectorFile_b,      "detector - reference for alignment");
    zap['A'] = make_field(overwriteDetector,   "overwrite detector file provided through '-a' with modified positions.");
    zap['T'] = make_field(tripodFile,          "tripods")                               = "";
    zap['s'] = make_field(sigma_m)               = 0.2;
    zap['O'] = make_field(option,              "keep strings vertical");
    zap['r'] = make_field(range,               "range of floors used in fit")           = JRange_t();
    zap['P'] = make_field(parameters,          "fit parameters (empty = all), possible values: " << JEEPZ() << PARAMETERS) = JPARSER::initialised();
    zap['d'] = make_field(debug)                 = 2;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  if (overwriteDetector) {
    if (tripodFile == "") {
      FATAL("No tripod file.");
    }
  }
  
  for (const auto& i : parameters) {
    if (PARAMETERS.count(i) == 0) {
      FATAL("Invalid fit paramater " << i << ", possible values: " << JEEPZ() << PARAMETERS << 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);
  }
 
 
  const JFit_t fit(detector_b, sigma_m, option, range);
 
  vector<JModule> data;
 
  for (JDetector::const_iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
    if (fit.hasModule(module->getID())) {
      data.push_back(*module);
    }
  }
 
 
  JSimplex<JModel_t> simplex;
  
  JSimplex<JModel_t>::debug              = debug;
  JSimplex<JModel_t>::MAXIMUM_ITERATIONS = 100000;
 
  simplex.value = JModel_t(0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
 
 
  if (parameters.empty()) {
    parameters = PARAMETERS;
  }
 
  if (parameters.count(X_t))   { simplex.step.push_back(JModel_t(0.01, 0.00, 0.00, 0.0,    0.0, 0.0)); }
  if (parameters.count(Y_t))   { simplex.step.push_back(JModel_t(0.00, 0.01, 0.00, 0.0,    0.0, 0.0)); }
  if (parameters.count(Z_t))   { simplex.step.push_back(JModel_t(0.00, 0.00, 0.01, 0.0,    0.0, 0.0)); }
  if (parameters.count(PHI_t)) { simplex.step.push_back(JModel_t(0.00, 0.00, 0.00, 5.0e-4, 0.0, 0.0)); }
  if (parameters.count(TX_t))  { simplex.step.push_back(JModel_t(0.00, 0.00, 0.00, 0.0, 1.0e-4, 0.0)); }
  if (parameters.count(TY_t))  { simplex.step.push_back(JModel_t(0.00, 0.00, 0.00, 0.0, 0.0, 1.0e-4)); }
 
 
  const double chi2 = simplex(fit, data.begin(), data.end());
 
  cout << "Number of iterations " << simplex.numberOfIterations << endl;
  cout << "chi2/NDF " << FIXED(7,3) << chi2 << '/' << (detector_a.size() - simplex.step.size()) << endl;
 
  cout << "model:" << endl;
 
  cout << "x    " << FIXED(7,3) << simplex.value.x   << endl;
  cout << "y    " << FIXED(7,3) << simplex.value.y   << endl;
  cout << "z    " << FIXED(7,3) << simplex.value.z   << endl;
  cout << "phi  " << FIXED(9,5) << simplex.value.phi << endl;
  cout << "Tx   " << FIXED(9,5) << simplex.value.tx  << endl;
  cout << "Ty   " << FIXED(9,5) << simplex.value.ty  << endl;
 
 
  if (overwriteDetector) {
 
    NOTICE("Store alignment data on files " << detectorFile_a << " and " << tripodFile << endl);
 
    detector_a.comment.add(JMeta(argc,argv));
 
    for (JDetector::iterator module = detector_a.begin(); module != detector_a.end(); ++module) {
 
      const JPosition3D p1 = module->getPosition();
      const JPosition3D p2 = fit.getPosition(simplex.value, p1);
 
      module->add(p2 - p1);
    }
 
    try {
      store(detectorFile_a, detector_a);
    }
    catch(const JException& error) {
      FATAL(error);
    }
 
    JContainer< vector<JTripod> > tripods;
 
    tripods.load(tripodFile.c_str());
 
    tripods.comment.add(JMeta(argc,argv));
    
    for (vector<JTripod>::iterator tripod = tripods.begin(); tripod != tripods.end(); ++tripod) {
      
      const JPosition3D p1 = tripod->getUTMPosition() - detector_a.getUTMPosition();
      const JPosition3D p2 = fit.getPosition(simplex.value, p1);
 
      tripod->add(p2 - p1);
    }
    
    tripods.store(tripodFile.c_str());
  }
}