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;
}