JMorphology.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <memory>
#include <map>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TRandom3.h"
#include "TString.h"
#include "TF1.h"
 
#include "JGeometry3D/JAngle3D.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JGeometry3D/JRotation3D.hh"
 
#include "JMath/JMathToolkit.hh"
 
#include "JLang/JVectorize.hh"
 
#include "JROOT/JGraph.hh"
#include "JROOT/JRootToolkit.hh"
 
#include "JSystem/JDateAndTime.hh"
 
#include "JTools/JAbstractHistogram.hh"
#include "JTools/JRange.hh"
 
#include "JAstronomy/JAstronomy.hh"
#include "JAstronomy/JMorphology.hh"
#include "JAstronomy/JResolution.hh"
 
#include "Jeep/JProperties.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * Toy simulation to study morphology of source and detector resolution.
 */
int main(int argc, char**argv)
{
  using namespace std;
  using namespace JPP;
  
  typedef JAbstractHistogram<double> JHistogram_t;
  typedef JRange<JDateAndTime>       UTC_t;
  typedef JRange<double>             range_type;
 
  struct {
    range_type ct;      //!< cosine zenith angle of track on in detector
  } parameters;
  
  size_t           numberOfEvents;
  string           outputFile;
  JHistogram_t     X;
  JHistogram_t     Y;
  UTC_t            UTC;
  morphology_type  morphology;
  resolution_type  resolution;
  double           background;
  double           omega;
  UInt_t           seed;
  int              debug;
 
  try {
 
    JProperties properties;
 
    properties.insert(gmake_property(parameters.ct));
    
    JParser<> zap;
 
    zap['@'] = make_field(properties)                                  = JPARSER::initialised(); 
    zap['o'] = make_field(outputFile)                                  = "source.root";
    zap['n'] = make_field(numberOfEvents);
    zap['x'] = make_field(X,              "rigth ascension [deg]")     =  JHistogram_t(500, -5.0, +5.0);
    zap['y'] = make_field(Y,              "declination [deg]")         =  JHistogram_t(500, -5.0, +5.0);
    zap['U'] = make_field(UTC,            "UTC time range");
    zap['M'] = make_field(morphology,     "morphology: \"<type> (parameters)+\"\n\twhere <type> can be:" << get_keys(morphology_helper));
    zap['r'] = make_field(resolution,     "resolution: \"<type> (parameters)+\"\n\twhere <type> can be:" << get_keys(resolution_helper));
    zap['B'] = make_field(background,     "number of events / s / sr") = 0.0;
    zap['O'] = make_field(omega,          "effective bin width [sr]")  = 1.0e-4;
    zap['S'] = make_field(seed)                = 0;
    zap['d'] = make_field(debug)               = 1;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
  gRandom->SetSeed(seed);
 
 
  const JAstronomy astronomy(ARCA);
    
  const time_t  Tmin_s  =  UTC.getLowerLimit().getTime();
  const time_t  Tmax_s  =  UTC.getUpperLimit().getTime();
 
  double x0 = -1.0;
  {
    const double xmin  = getRadians(X.getLowerLimit());
    const double xmax  = getRadians(X.getUpperLimit());
    const double ymin  = sin(getRadians(Y.getLowerLimit()));
    const double ymax  = sin(getRadians(Y.getUpperLimit()));
 
    const double V  = (ymax - ymin) * (xmax - xmin);
 
    x0 =  1.0 - V / (sqrt(2.0)*2.0*PI);
  }
  const double dx = omega / (2.0*PI);
  const int    nx = (int) ((1.0 - x0) / dx);
 
  TH1D ha("ha", NULL,  70, -3.0,         +2.3);    // log10(alpha)
  TH1D h1("h1", NULL,  nx,  1.0 - nx*dx, +1.0);    // cos(alpha)
  TH2D h2("h2", NULL,
          X.getNumberOfBins(), X.getLowerLimit(), X.getUpperLimit(),
          Y.getNumberOfBins(), Y.getLowerLimit(), Y.getUpperLimit());
 
  JGraph_t g1;                             // background
  JGraph_t g2;                             // signal
 
  TH2D hs("hs", NULL, 100, 0.0, NUMBER_OF_SECONDS_PER_SEDERIAL_DAY / NUMBER_OF_SECONDS_PER_HOUR, 100, -1.0, +1.0);
  
  for (size_t counter = 0; counter != numberOfEvents; ++counter) {
 
    STATUS("counter: "<< setw(8) << counter << '\r'); DEBUG(endl);
 
 
    // generate neutrino at random time according morphology of source
 
    const double             t1  =  getTimeMJD(Tmin_s + gRandom->Rndm() * (Tmax_s - Tmin_s));
    const JSourceLocation    u1  =  morphology->get();
    
 
    // neutrino in detector coordinate system
 
    const JNeutrinoDirection u2  =  astronomy.getNeutrinoDirection(t1, u1);
 
 
    // track in detector coordinate system according experimental resolution
 
    JDirection3D u3 = resolution->get();
 
 
    // align track with respect to incident neutrino    
    {
      const JRotation3D R(u2);
 
      u3.rotate_back(R);
    }
 
 
    // apply cuts
 
    hs.Fill(fmod(t1, NUMBER_OF_SECONDS_PER_SEDERIAL_DAY) / NUMBER_OF_SECONDS_PER_HOUR, u3.getDZ());
 
    if (!parameters.ct(u3.getDZ())) {
      continue;
    }
 
 
    // track in astrophysical coordinates
    
    const JNeutrinoDirection u4(u3);
 
    const JSourceLocation    u5  =  astronomy.getSourceLocation(t1, u4);
    
 
    // evaluate 2D-residuals in source coordinates
 
    const JRotation3D Rs(morphology->getSourceLocation());
 
    JDirection3D us(morphology->getSourceLocation()); 
    JDirection3D vs(u5);
 
    us.rotate(Rs);
    vs.rotate(Rs);
 
    const double x = getDegrees(asin(vs.getDX() - us.getDX()));
    const double y = getDegrees(asin(vs.getDY() - us.getDY()));
 
    h2.Fill(x, y);
    g2.put (x, y);
 
    const double dot = getDot(us, vs);
 
    ha.Fill(log10(getDegrees(acos(dot))));
    h1.Fill(dot);
  }
  STATUS(endl);
 
 
  if (background > 0.0) {
 
    NOTICE("Background: " << FIXED(9,5) << background * (Tmax_s - Tmin_s) * omega << " events / bin" << endl);
 
    for (Int_t ix = 1; ix <= h2.GetXaxis()->GetNbins(); ++ix) {
      for (Int_t iy = 1; iy <= h2.GetYaxis()->GetNbins(); ++iy) {
 
        const double xmin  = getRadians(h2.GetXaxis()->GetBinLowEdge(ix));
        const double xmax  = getRadians(h2.GetXaxis()->GetBinUpEdge (ix));
        const double ymin  = sin(getRadians(h2.GetYaxis()->GetBinLowEdge(iy)));
        const double ymax  = sin(getRadians(h2.GetYaxis()->GetBinUpEdge (iy)));
 
        const double omega = (ymax - ymin) * (xmax - xmin);
 
        const double mu = background * omega * (Tmax_s - Tmin_s);
 
        const double x  = h2.GetXaxis()->GetBinCenter(ix);
        const double y  = h2.GetYaxis()->GetBinCenter(iy);
        const double z  = gRandom->Poisson(mu);
 
        h2.Fill(x, y, z);
      }
    }
 
    const double xmin  = getRadians(X.getLowerLimit());
    const double xmax  = getRadians(X.getUpperLimit());
    const double ymin  = sin(getRadians(Y.getLowerLimit()));
    const double ymax  = sin(getRadians(Y.getUpperLimit()));
 
    const double omega = (ymax - ymin) * (xmax - xmin);
 
    const double mu = background * omega * (Tmax_s - Tmin_s);
 
    const size_t N  = gRandom->Poisson(mu);
 
    for (size_t i = 0; i != N; ++i) {
 
      const double x = gRandom->Uniform(xmin, xmax);
      const double y = gRandom->Uniform(ymin, ymax);
 
      g1.put(getDegrees(x), getDegrees(asin(y)));
 
      const double dot = sqrt(1.0 - y*y - sin(x)*sin(x));
 
      h1.Fill(dot);
    }
  }
 
 
  TFile out(outputFile.c_str(), "recreate");
 
  out << ha << h1 << h2 << hs;
  out << JGraph(g1, "g1");
  out << JGraph(g2, "g2");
  
  out.Write();
  out.Close();
}