JPrintRMS.cc

Go to the documentation of this file.

#include <iostream>
#include <iomanip>

#include "TROOT.h"
#include "TFile.h"

#include "JSupport/JMultipleFileScanner.hh"

#include "JAcoustics/JEvent.hh"
#include "JAcoustics/JSupport.hh"

#include "JTools/JQuantile.hh"

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


/**
 * \file
 *
 * Auxiliary program to print RMS of acoustic events.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  JMultipleFileScanner<JEvent> inputFile;
  JLimit_t&   numberOfEvents = inputFile.getLimit();
  size_t      numberOfOutliers;
  int         debug;

  try {

    JParser<> zap("Auxiliary program to print RMS of acoustic events.");

    zap['f'] = make_field(inputFile);
    zap['n'] = make_field(numberOfEvents)   = JLimit_t::max();
    zap['N'] = make_field(numberOfOutliers) = 0;
    zap['d'] = make_field(debug)            = 1;

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


  JQuantile Q;

  while (inputFile.hasNext()) {

    const JEvent* evt = inputFile.next();

    vector<double> buffer;

    for (JEvent::const_iterator i = evt->begin(); i != evt->end(); ++i) {
      buffer.push_back(i->getToE());
    }

    if (buffer.size() >= 2 * numberOfOutliers + 5) {

      sort(buffer.begin(), buffer.end());

      vector<double>::iterator         p = buffer.begin();
      vector<double>::reverse_iterator q = buffer.rbegin();

      for (size_t i = 0; i != numberOfOutliers; ++i) {
        ++p;
        ++q;
      }

      const double x0  = 0.5 * (*p + *q);

      for (vector<double>::iterator i = p; i != q.base(); ++i) {
        *i -= x0;
      }

      const JQuantile q1("", p, q.base());

      Q.put(q1.getSTDev());
    }
  }

  cout << Q.getMean();
}