examples/JDynamics/JBallarat.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <limits>
#include <map>

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

#include "JROOT/JGraph.hh"
#include "JROOT/JRootToolkit.hh"
#include "JROOT/JManager.hh"

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JDetectorSupportkit.hh"
#include "JDetector/JModuleRouter.hh"

#include "JSupport/JMultipleFileScanner.hh"
#include "JCompass/JEvt.hh"
#include "JCompass/JSupport.hh"

#include "JDynamics/JDynamics.hh"

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


namespace {

  using JUTC::JUTCTimeRange;
  using JGEOMETRY3D::JQuaternion3D;

  /**
   * Auxiliary method to get UTC time range of calibration data.
   *
   * \param  __begin                   begin of calibration data
   * \param  __end                     end   of calibration data
   * \return                           UTC time range
   */
  template<class T>
  static inline JUTCTimeRange getUTCTimeRange(T __begin, T __end)
  {
    JUTCTimeRange range = JUTCTimeRange::DEFAULT_RANGE();

    for (T i = __begin; i != __end; ++i) {
      if (!i->second.empty()) {
        range.combine(JUTCTimeRange(i->second.getXmin(), i->second.getXmax()));
      }
    }
 
    return range;
  }
  
  /**
   * Get twist angle.
   *
   * \param  Q                         quaternion
   * \return                           twist angle [rad]
   */
  inline double getTwist(const JQuaternion3D& Q)
  {
    using namespace JPP;

    const JQuaternion3D::decomposition q(Q, JVector3Z_t);

    double phi = (q.twist.getD() > 0.0 ? +2.0 : -2.0) * acos(q.twist.getA());

    while (phi > +PI) { phi -= 2*PI; }
    while (phi < -PI) { phi += 2*PI; }

    return phi;
  }
}


/**
 * \file
 *
 * Example program to plot orientation data.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  JMultipleFileScanner<JCOMPASS::JOrientation> inputFile;
  JLimit_t&       numberOfEvents = inputFile.getLimit();
  string          detectorFile;
  string          outputFile;
  double          Tmax_s;
  int             debug;

  try {

    JParser<> zap("Example program to plot orientation data.");

    zap['f'] = make_field(inputFile,    "output of JBallarat[.sh]");
    zap['n'] = make_field(numberOfEvents)                                    = JLimit::max();
    zap['a'] = make_field(detectorFile);
    zap['o'] = make_field(outputFile)    = "gold.root";
    zap['T'] = make_field(Tmax_s);
    zap['d'] = make_field(debug)         = 2;

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


  JDetector detector;

  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }

  const JModuleRouter router(detector);

  if (!hasDetectorAddressMap(detector.getID())) {
    FATAL("No detector address map for detector identier " << detector.getID() << endl);
  }

  const JDetectorBuilder& demo = getDetectorBuilder(detector.getID());


  JDynamics dynamics(detector, Tmax_s);

  STATUS("loading input from file(s) " << inputFile << "... " << flush);

  dynamics.load(inputFile);

  STATUS("OK" << endl);


  map<int, JGraph_t> Z0;

  for (JDynamics::JOrientation::const_iterator module = dynamics.orientation.begin(); module != dynamics.orientation.end(); ++module) {

    if (router.getModule(module->first).getFloor() != 0) {

      STATUS("Creating graphics for input data " << setw(10) << module->first << flush << '\r');

      JQuaternion3D Q0 = getRotation(demo.getModule(module->first), router.getModule(module->first));

      Q0.conjugate();

      JGraph_t& z0 = Z0[module->first];

      for (JDynamics::JOrientation::function_type::const_iterator i = module->second.begin(); i != module->second.end(); ++i) {

        const JQuaternion3D Q1 = router.getModule(module->first).getQuaternion() * i->getY();

        z0.put(i->getX(), getTwist(Q1 * Q0));
      }
    }
  }
  STATUS(endl);

  map<int, TH1D*> H1;

  for (JDynamics::JOrientation::const_iterator module = dynamics.orientation.begin(); module != dynamics.orientation.end(); ++module) {

    if (router.getModule(module->first).getFloor() != 0 && module->second.size() > 1u) {

      STATUS("Creating graphics for input data " << setw(10) << module->first << flush << '\r');

      TH1D* h1 = H1[module->first] = new TH1D(MAKE_CSTRING("U[" << module->first << "].twist"), NULL, 500, -180.0, +180.0);

      for (JDynamics::JOrientation::function_type::const_iterator i1 = module->second.begin(), i0 = i1++; i1 != module->second.end(); ++i0, ++i1) {

        JQuaternion3D Q0 = i0->getY();
        JQuaternion3D Q1 = i1->getY();

        h1->Fill(getTwist(Q1 * Q0.conjugate()) * 180.0 / PI);
      }
    }
  }
  STATUS(endl);

  map<int, TH1D*> H0;

  for (JDynamics::JOrientation::const_iterator module = dynamics.orientation.begin(); module != dynamics.orientation.end(); ++module) {

    STATUS("Creating graphics for compiled data " << setw(10) << module->first << flush << '\r');

    if (router.getModule(module->first).getFloor() != 0 && !module->second.empty()) {

      JQuaternion3D Q0 = getRotation(demo.getModule(module->first), router.getModule(module->first));

      Q0.conjugate();

      const Double_t xmin = module->second.getXmin();
      const Double_t xmax = module->second.getXmax();
      const Int_t    nx   = (Int_t) ((xmax - xmin) / Tmax_s);

      TH1D* h0 = H0[module->first] = new TH1D(MAKE_CSTRING("H[" << module->first << "].twist"), NULL, nx, xmin, xmax);

      for (Int_t i = 1; i <= h0->GetXaxis()->GetNbins(); ++i) {

        const Double_t      x = h0->GetXaxis()->GetBinCenter(i);
        const JQuaternion3D Q = module->second(x) * Q0;

        h0->SetBinContent(i, getTwist(Q));
      }
    }
  }
  STATUS(endl);

  const Double_t xmin = getUTCTimeRange(dynamics.orientation.begin(), dynamics.orientation.end()).getLowerLimit();
  const Double_t xmax = getUTCTimeRange(dynamics.orientation.begin(), dynamics.orientation.end()).getUpperLimit();
  const Int_t    nx   = (Int_t) ((xmax - xmin) / Tmax_s);

  JManager<int, TH1D> X0(new TH1D("X[%].twist", NULL, nx, xmin, xmax));

  for (Int_t i = 1; i <= X0->GetXaxis()->GetNbins(); ++i) {

    STATUS("Creating graphics for detector data " << FIXED(5,1) << (double) (i * 100) / (double) X0->GetXaxis()->GetNbins() << "%" << flush << '\r');

    dynamics.update(X0->GetXaxis()->GetBinCenter(i));

    for (JDetector::const_iterator module = dynamics.begin(); module != dynamics.end(); ++module) {

      if (module->getFloor() != 0) {
      
        TH1D* x0 = X0[module->getID()];

        const JQuaternion3D Q = getRotation(router.getModule(module->getID()), *module);

        x0->SetBinContent(i, getTwist(Q));
      }
    }
  }
  STATUS(endl);


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

  for (map<int, TH1D*>::const_iterator i = H0.begin(); i != H0.end(); ++i) {
    out << *i->second;
  }

  for (map<int, TH1D*>::const_iterator i = H1.begin(); i != H1.end(); ++i) {
    out << *i->second;
  }

  out << X0;

  for (map<int, JGraph_t>::const_iterator i = Z0.begin(); i != Z0.end(); ++i) {
    out << JGraph(i->second, MAKE_CSTRING("G[" << i->first << "].twist"));
  }

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