JOmega3D.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <cmath>

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

#include "JMath/JConstants.hh"
#include "JTools/JQuantile.hh"
#include "JTools/JRange.hh"
#include "JROOT/JRootToolkit.hh"

#include "JGeometry3D/JOmega3D.hh"

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


/**
 * \file
 *
 * Example program to test JGEOMETRY3D::JOmega3D class.
 * \author mdejong
 */
int main(int argc, char**argv)
{
  using namespace std;
  using namespace JPP;

  string         outputFile;
  int            numberOfEvents;
  JRange<double> range_Deg;
  int            debug;

  try {

    JParser<> zap("Example program to test JOmega3D class.");

    zap['o'] = make_field(outputFile)          = "";
    zap['n'] = make_field(numberOfEvents)      = 100;
    zap['G'] = make_field(range_Deg)           = JRange<double>(0.5, 10.5);
    zap['d'] = make_field(debug)               = 3;

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


  const int N = 10;
  
  TH1D h0("h0[grid",       NULL, N, range_Deg.getLowerLimit(), range_Deg.getUpperLimit());
  TH1D h1("h1[RMS]",       NULL, N, range_Deg.getLowerLimit(), range_Deg.getUpperLimit());
  TH1D h2("h2[deviation]", NULL, N, range_Deg.getLowerLimit(), range_Deg.getUpperLimit());

  
  vector<JOmega3D>  omega;
  vector<JQuantile> quantile;

  for (int i = 1; i <= h0.GetNbinsX(); ++i) {

    const Double_t x = h0.GetBinCenter(i);

    omega   .push_back(JOmega3D(x * PI/180.0));
    quantile.push_back(JQuantile());    
  }

  
  for (int i = 0; i != numberOfEvents; ++i) {

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

    const double   theta = gRandom->Uniform(0.0, PI);
    const double   phi   = gRandom->Uniform(0.0, 2*PI);

    const JAngle3D angle(theta, phi);

    for (size_t j = 0; j != omega.size(); ++j) {
      
      const int      pos = omega[j].find(angle);
      const double   dot = angle.getDot(omega[j][pos]);

      quantile[j].put(acos(dot) * 180.0/PI); 
    }
  }
  STATUS(endl);

  
  for (int i = 1; i <= h0.GetNbinsX(); ++i) {
    h0.SetBinContent(i, omega   [i-1].size());
    h1.SetBinContent(i, quantile[i-1].getSTDev());
    h2.SetBinContent(i, quantile[i-1].getDeviation());
  }

  if (outputFile != "") {
    
    TFile out(outputFile.c_str(), "recreate");

    out << h0 << h1 << h2;
    
    out.Write();
    out.Close();
  }

  for (int i = 1; i <= h2.GetNbinsX(); ++i) {

    const Double_t x = h2.GetBinCenter (i);
    const Double_t y = h2.GetBinContent(i);

    NOTICE("Grid test " << x << " [deg]: " << y << " [deg]." << endl); 
    ASSERT(y < x / sqrt(2));
  }

  return 0;
}