JAlignDetector.cc File Reference

Auxiliary program to align two detectors. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <vector>
#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

mjongen

Definition in file JAlignDetector.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 288 of file JAlignDetector.cc.

{
  using namespace std;
  using namespace JPP;

  string             detectorFile_a;
  string             detectorFile_b;
  bool               overwriteDetector;
  string             tripodFile;
  double             sigma_m;
  bool               option;
  JRange_t           range;
  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['d'] = make_field(debug)                 = 2;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  if (overwriteDetector) {
    if (tripodFile == "") {
      FATAL("No tripod file.");
    }
  }
  
  
  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()) && module->getFloor() != 0) {
      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);

  simplex.step.resize(6);

  simplex.step[0] = JModel_t(0.01, 0.00, 0.00, 0.0,    0.0, 0.0);
  simplex.step[1] = JModel_t(0.00, 0.01, 0.00, 0.0,    0.0, 0.0);
  simplex.step[2] = JModel_t(0.00, 0.00, 0.01, 0.0,    0.0, 0.0);
  simplex.step[3] = JModel_t(0.00, 0.00, 0.00, 5.0e-4, 0.0, 0.0);
  simplex.step[4] = JModel_t(0.00, 0.00, 0.00, 0.0, 1.0e-4, 0.0);
  simplex.step[5] = 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());
  }
}