JRose.cc

Go to the documentation of this file.

#include <iostream>
#include <iomanip>
#include <vector>
#include <map>

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

#include "JSupport/JMultipleFileScanner.hh"
#include "JROOT/JManager.hh"
#include "JROOT/JGraph.hh"
#include "JLang/JTitle.hh"

#include "JCompass/JModel.hh"
#include "JCompass/JEvt.hh"
#include "JCompass/JEvtToolkit.hh"
#include "JCompass/JSupport.hh"

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


namespace {

  /**
   * Get rotation angle of given quaternion.
   *
   * \param  Q             quaternion
   * \return               angle [rad]
   */
  inline double getAngle(const JGEOMETRY3D::JQuaternion3D& Q)
  {
    using namespace std;
    using namespace JPP;

    return atan2(sqrt(Q.getB()*Q.getB() +
                      Q.getC()*Q.getC() +
                      Q.getD()*Q.getD()),
                 Q.getA()) * 2.0;
  }

  
  /**
   * Get signed rotation angle of given quaternion.
   *
   * \param  Q             quaternion
   * \param  v             direction
   * \return               angle [rad]
   */
  inline double getAngle(const JGEOMETRY3D::JQuaternion3D& Q, const JGEOMETRY3D::JVector3D& v)
  {
    using namespace std;
    using namespace JPP;

    return getAngle(Q) * (v.getDot(Q) >= 0.0 ? +1.0 : -1.0);
  }


  /**
   * Auxiliary data structure for organising TGraph's.
   */
  struct map_type :
    public std::map<int, JROOT::JGraph_t>,
    public JLANG::JTitle
  {
    /**
     * Constructor.
     *
     * \param  title         title
     */
    map_type(const std::string& title) :
      JTitle(title)
    {}
  };
}

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

  JMultipleFileScanner<JEvt> inputFile;
  JLimit_t&       numberOfEvents = inputFile.getLimit();
  string          outputFile;
  double          stdev          = numeric_limits<double>::max();      // maximal chi2 / NDF
  int             debug;

  try {

    JProperties properties;

    properties.insert(gmake_property(stdev));

    JParser<> zap("Example program to plot compass fit results.");

    zap['f'] = make_field(inputFile,    "output of JCompass");
    zap['n'] = make_field(numberOfEvents)                                    = JLimit::max();
    zap['@'] = make_field(properties)                                        = JPARSER::initialised();
    zap['o'] = make_field(outputFile)    = "rose.root";
    zap['d'] = make_field(debug)         = 2;

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


  JManager<int, TH1D> H0(new TH1D("H0[%]", NULL, 100, 0.0, +0.200));
  JManager<int, TH1D> H1(new TH1D("H1[%]", NULL, 100, 0.0, +0.005));
  JManager<int, TH1D> HC(new TH1D("HC[%]", NULL, 300,-2.0, +1.0));

  map_type G0("Q0.twist");
  map_type G1("Q1.twist");
  map_type GA("Q0.swing");
  map_type GB("Q0.atan2");


  for (JModel previous; inputFile.hasNext(); ) {

    STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);

    const JEvt*  evt   = inputFile.next();
    const JModel model = getModel(*evt);

    if (inputFile.getCounter() > 1) {
      H0[evt->id]->Fill(getAngle(model.Q0, previous.Q0));
      H1[evt->id]->Fill(getAngle(model.Q1, previous.Q1));
      HC[evt->id]->Fill(log10(evt->chi2 / evt->ndf));
    }

    previous = model;

    if (evt->chi2 / evt->ndf <= stdev) {

      const JQuaternion3D::decomposition q0(model.Q0, JVector3Z_t);
      const JQuaternion3D::decomposition q1(model.Q1, JVector3Z_t);

      const double t1 = 0.5e-3 * (evt->UNIXTimeStart + evt->UNIXTimeStop);

      G0[evt->id].put(t1, getAngle(q0.twist, JVector3Z_t));
      G1[evt->id].put(t1, getAngle(q1.twist, JVector3Z_t));
      GA[evt->id].put(t1, getAngle(q0.swing));
      GB[evt->id].put(t1, atan2(q0.swing.getB(), q0.swing.getC()));   // conform JAcoustics::JModel::getAngle
    }
  }
  STATUS(endl);


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

  for (JManager<int, TH1D>* h1 : { &H0, &H1, &HC }) {
    for (JManager<int, TH1D>::const_iterator i = h1->begin(); i != h1->end(); ++i) {
      out << *(i->second);
    }
  }

  for (map_type* g1 : { &G0, &G1, &GA, &GB }) {
    for (map_type::const_iterator i = g1->begin(); i != g1->end(); ++i) {
      out << JGraph(i->second, MAKE_CSTRING("[" << i->first << "]." << g1->getTitle()));
    }
  }

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