JFootprint.cc

Go to the documentation of this file.

#include <iostream>
#include <iomanip>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF2.h"
#include "TFitResult.h"

#include "JROOT/JMinimizer.hh"

#include "JSupport/JMultipleFileScanner.hh"

#include "JTools/JAbstractHistogram.hh"
#include "JTools/JRange.hh"

#include "JAcoustics/JEvt.hh"
#include "JAcoustics/JSuperEvt.hh"
#include "JAcoustics/JSupport.hh"

#include "JLang/JObjectMultiplexer.hh"

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

namespace {

  const char* const _F2 = ".f2";  //!< extension of name for fit function
}

/**
 * \file
 *
 * Auxiliary application to determine tilt angles of seabed based on measured tilt angles.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  typedef JAbstractHistogram<Double_t> JHistogram_t;
  typedef JRange            <Double_t> JRange_t;

  typedef JTYPELIST<JEvt, JSuperEvt>::typelist  typelist;
  
  JMultipleFileScanner<typelist> inputFile;
  JLimit_t&           numberOfEvents = inputFile.getLimit();
  string              outputFile;
  JHistogram_t        X;
  JHistogram_t        Y;
  double              Z;
  JRange_t            x;
  JRange_t            y;
  string              option;
  bool                writeFits;
  int                 debug;

  try {

    JParser<> zap("Auxiliary application to determine tilt angles of seabed based on measured tilt angles.");

    zap['f'] = make_field(inputFile,       "input file (output of JKatoomba[.sh]/JFremantle[.sh])");
    zap['n'] = make_field(numberOfEvents)   = JLimit::max();
    zap['o'] = make_field(outputFile)       = "footprint.root";
    zap['x'] = make_field(X,               "histogram x abscissa")                             = JHistogram_t(1000, -25.0, +25.0);
    zap['y'] = make_field(Y,               "histogram y abscissa")                             = JHistogram_t(1000, -25.0, +25.0);
    zap['X'] = make_field(x,               "fit x-range centered at top [mrad]")               = JRange_t(-2.0, +2.0);
    zap['Y'] = make_field(y,               "fit y-range centered at top [mrad]")               = JRange_t(-2.0, +2.0);
    zap['Z'] = make_field(Z,               "detector height (for 2nd order tilt correction)")  = 0.0;
    zap['O'] = make_field(option,          "ROOT fit option, see TH1::Fit.")                   = "L";
    zap['w'] = make_field(writeFits,       "write fit function(s) to ROOT file " << "(\"<name>" << _F2 << "\")");
    zap['d'] = make_field(debug)            = 2;

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

  Z *= 0.5; // half height

  if (option.find('O') == string::npos) { option += "O"; }
  if (option.find("R") == string::npos) { option += "R"; }
  if (option.find("S") == string::npos) { option += "S"; }
  //if (option.find('N') == string::npos) { option += "N"; }
  if (debug ==  0 && option.find('Q') == string::npos) { option += "Q"; }

  TH2D H2("%", NULL,
          X.getNumberOfBins(), X.getLowerLimit(), X.getUpperLimit(),
          Y.getNumberOfBins(), Y.getLowerLimit(), Y.getUpperLimit());

  JObjectMultiplexer<typelist, JEvt> in(inputFile);

  for (counter_type counter = 0; in.hasNext() && counter != inputFile.getLimit(); ++counter) {

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

    const JEvt* evt = in.next();

    if (!evt->empty()) {

      double Tx = 0.0;
      double Ty = 0.0;
    
      for (JEvt::const_iterator i = evt->begin(); i != evt->end(); ++i) {

        const double tx = (i->tx + i->tx2 * Z) * 1.0e3;  // [mrad]
        const double ty = (i->ty + i->ty2 * Z) * 1.0e3;  // [mrad]

        Tx += tx;
        Ty += ty;
      }

      Tx /= evt->size();
      Ty /= evt->size();

      H2.Fill(Tx, Ty);
    }
  }
  STATUS(endl);


  double x0 = 0.0;
  double y0 = 0.0;
  double z0 = 0.0;
  
  for (Int_t ix = 1; ix <= H2.GetXaxis()->GetNbins(); ++ix) {
    for (Int_t iy = 1; iy <= H2.GetYaxis()->GetNbins(); ++iy) {
      if (H2.GetBinContent(ix,iy) > z0) {
        x0 = H2.GetXaxis()->GetBinCenter(ix);
        y0 = H2.GetYaxis()->GetBinCenter(iy);
        z0 = H2.GetBinContent(ix,iy);
      }
    }
  }

  STATUS("top: " << FIXED(7,3) << x0 << ' ' << FIXED(7,3) << y0 << ' ' << FIXED(9,1) << z0 << endl);


  TF2 f2("f2", "[0] * exp(-0.5 * (x-[1])*(x-[1]) / ([2]*[2])) * exp(-0.5 * (y-[3])*(y-[3]) / ([4]*[4])) + [5] + [6]*x + [7]*y");

  f2.SetParameter(0, z0 * 0.7);
  f2.SetParameter(1, x0);
  f2.SetParameter(2, 0.5);
  f2.SetParameter(3, y0);
  f2.SetParameter(4, 0.5);
  f2.SetParameter(5, 0.0);
  f2.SetParameter(6, 0.0);
  f2.SetParameter(7, 0.0);

  f2.SetParError(0, 0.1);
  f2.SetParError(1, 0.02);
  f2.SetParError(2, 0.02);
  f2.SetParError(3, 0.02);
  f2.SetParError(4, 0.02);
  f2.SetParError(5, 0.001);
  f2.SetParError(6, 0.001);
  f2.SetParError(7, 0.001);

  f2.SetRange(x0 + x.getLowerLimit(), y0 + y.getLowerLimit(), x0 + x.getUpperLimit(), y0 + y.getUpperLimit());

  H2.GetXaxis()->SetRangeUser(x0 + x.getLowerLimit(), x0 + x.getUpperLimit());
  H2.GetYaxis()->SetRangeUser(y0 + y.getLowerLimit(), y0 + y.getUpperLimit());

  TFitResultPtr result = H2.Fit(&f2, option.c_str(), "same");
  
  H2.GetXaxis()->SetRangeUser(H2.GetXaxis()->GetXmin(), H2.GetXaxis()->GetXmax());
  H2.GetYaxis()->SetRangeUser(H2.GetYaxis()->GetXmin(), H2.GetYaxis()->GetXmax());

  if (result.Get() == NULL) {
    FATAL("Invalid TFitResultPtr " << H2.GetName() << endl);
  }

  if (debug >= status_t || !result->IsValid()) {
    cout << "chi2/NDF " << FIXED(7,3) << result->Chi2() << '/' << result->Ndf() << ' ' << (result->IsValid() ? "" : "failed") << endl;
  }

  cout << "tilt: " << FIXED(9,6) << f2.GetParameter(1)*1.0e-3 << ' ' << FIXED(9,6) << f2.GetParameter(3)*1.0e-3 << endl;

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

  out << H2;

  if (writeFits) {
    
    TH2D* h2 = (TH2D*) H2.Clone(MAKE_CSTRING(H2.GetName() << _F2));

    h2->Reset();

    for (Int_t ix = 1; ix <= h2->GetXaxis()->GetNbins(); ++ix) {
      for (Int_t iy = 1; iy <= h2->GetYaxis()->GetNbins(); ++iy) {

        const Double_t x = h2->GetXaxis()->GetBinCenter(ix);
        const Double_t y = h2->GetYaxis()->GetBinCenter(iy);

        h2->SetBinContent(ix, iy, f2.Eval(x,y));
        h2->SetBinError  (ix, iy, 0.0);
      }
    }
  }

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