JTwist.cc

Go to the documentation of this file.

#include <iostream>
#include <iomanip>
#include <cmath>
#include <vector>

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

#include "JLang/JException.hh"
#include "JMath/JConstants.hh"
#include "JGeometry3D/JQuaternion3D.hh"
#include "JGeometry3D/JEigen3D.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JCompass/JHit.hh"
#include "JCompass/JModel.hh"
#include "JFit/JSimplex.hh"

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


/**
 * Test program for fit of quaternions.
 */
int main(int argc, char**argv)
{
  using namespace std;
  using namespace JPP;

  string       outputFile;
  int          numberOfEvents;
  double       X;
  double       Y;
  double       Z;
  double       alpha;
  JPosition3D  sigma;
  UInt_t       seed;
  bool         fit;
  int          debug;

  try {

    JParser<> zap("Test program for fit of quaternions.");

    zap['o'] = make_field(outputFile)          = "twist.root";
    zap['n'] = make_field(numberOfEvents)      = 1000;
    zap['x'] = make_field(X,               "tilt angle around x-axis [deg]") = 0.0;
    zap['y'] = make_field(Y,               "tilt angle around y-axis [deg]") = 0.0;
    zap['z'] = make_field(Z,               "tilt angle around z-axis [deg]") = 0.0;
    zap['a'] = make_field(alpha,           "twist angle [deg/m]");
    zap['s'] = make_field(sigma,           "(x,y,z) resolution [deg]");
    zap['S'] = make_field(seed)                = 0;
    zap['F'] = make_field(fit,             "enable fit");
    zap['d'] = make_field(debug)               = 2;

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


  gRandom->SetSeed(seed);


  TFile out(outputFile.c_str(), "recreate");

  TH1D hc("chi2/NDF", NULL, 100, 0.0,  2.0);
  TH1D h0("tilt",     NULL, 100, 0.0,  2.0);
  TH1D h1("twist",    NULL, 100, 0.0,  0.1);


  // model

  const JModel model(JQuaternion3D(JQuaternion3X(X     * PI / 180.0),
                                   JQuaternion3Y(Y     * PI / 180.0),
                                   JQuaternion3Z(Z     * PI / 180.0)),
                     JQuaternion3D(JQuaternion3Z(alpha * PI / 180.0)));

  STATUS("Q0   " << FIXED(12,8) << model.Q0 << endl);
  STATUS("Q1   " << FIXED(12,8) << model.Q1 << endl);
  

  JSimplex<JModel> simplex;

  JSimplex<JModel>::MAXIMUM_ITERATIONS = 10000;

  simplex.debug = debug;

  const JChi2 getChi2(EM_NORMAL);


  for (int counter = 0; counter != numberOfEvents; ++counter) {

    STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl);

    // generate data

    const int    N   = 18;                         // number of floors
    const double z1  = 30.0;                       // height of first floor
    const double dz  = 10.0;                       // distance between floors
    const double eps = 1.0e-6;                     // minimal sigma

    vector<JHit> data;

    for (int i = 1; i <= N; ++i) {

      const double z = z1 + (i - 1) * dz;          // height of floor i

      data.push_back(JHit(i, z, model(z), max(eps, sigma.getLength())));
    }

    // smearing

    for (vector<JHit>::iterator i = data.begin(); i != data.end(); ++i) {
      i->mul(JQuaternion3D(JQuaternion3X(gRandom->Gaus(0.0, sigma.getX()) * PI / 180.0),
                           JQuaternion3Y(gRandom->Gaus(0.0, sigma.getY()) * PI / 180.0),
                           JQuaternion3Z(gRandom->Gaus(0.0, sigma.getZ()) * PI / 180.0)));
    }

    for (vector<JHit>::const_iterator i = data.begin(); i != data.end(); ++i) {
      DEBUG(FIXED(4,1) << i->getZ() << ' ' << FIXED(12,8) << i->getQuaternion() << endl);
    }


    JModel result(data.begin(), data.end());      // prefit

    DEBUG("R0   " << FIXED(12,8) << result.Q0 << ' ' << getAngle(model.Q0, result.Q0) << endl);
    DEBUG("R1   " << FIXED(12,8) << result.Q1 << ' ' << getAngle(model.Q1, result.Q1) << endl);

    if (fit) {

      simplex.value = result;                     // start value

      simplex.step.resize(4);

      simplex.step[0] = JModel(JQuaternion3X(0.50 * PI / 180.0), JQuaternion3D::getIdentity());
      simplex.step[1] = JModel(JQuaternion3Y(0.50 * PI / 180.0), JQuaternion3D::getIdentity());
      simplex.step[2] = JModel(JQuaternion3Z(0.50 * PI / 180.0), JQuaternion3D::getIdentity());
      simplex.step[3] = JModel(JQuaternion3D::getIdentity(), JQuaternion3Z(0.05 * PI / 180.0));

      const double chi2 = simplex(getChi2, data.begin(), data.end());  
      const int    ndf  = data.size() * 4 - simplex.step.size();

      DEBUG("number of iterations " << simplex.numberOfIterations << endl);

      result = simplex.value;                     // final value

      DEBUG("chi2 " << chi2 << '/' << ndf << endl);
      DEBUG("R0   " << FIXED(12,8) << simplex.value.Q0 << ' ' << getAngle(model.Q0, result.Q0) << endl);
      DEBUG("R1   " << FIXED(12,8) << simplex.value.Q1 << ' ' << getAngle(model.Q1, result.Q1) << endl);

      hc.Fill(chi2 / ndf);
    }

    //h0.Fill(getAngle(model.Q0, result.Q0));
    h0.Fill(getAngle(JDirection3D(0.0, 0.0, 1.0).rotate(model .Q0), 
                     JDirection3D(0.0, 0.0, 1.0).rotate(result.Q0)));
    h1.Fill(getAngle(model.Q1, result.Q1));
  }
  STATUS(endl);

  out.Write();
  out.Close();
}