#include #include #include #include #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, 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 inputFile; JLimit_t& numberOfEvents = inputFile.getLimit(); string outputFile; double stdev = numeric_limits::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 H0(new TH1D("H0[%]", NULL, 100, 0.0, +0.200)); JManager H1(new TH1D("H1[%]", NULL, 100, 0.0, +0.005)); JManager 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* h1 : { &H0, &H1, &HC }) { for (JManager::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(); }