JSignalToAWriter.cc

Go to the documentation of this file.

#include <string>
#include <iostream>

#include "TRandom3.h"
#include "TH1D.h"

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"

#include "JAcoustics/JToA.hh"
#include "JAcoustics/JSupport.hh"
#include "JAcoustics/JReceiver.hh"
#include "JAcoustics/JSoundVelocity.hh"

#include "JGeometry3D/JCylinder3D.hh"
#include "JGeometry3D/JRotation3D.hh"

#include "JTools/JRange.hh"

#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMeta.hh"

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

namespace JROOT {
  template<>
  class JRootCreateFlatTree<JACOUSTICS::JToA> : public std::true_type
  {};
}

namespace {

  /**
   * Get range of sound.
   *
   * \param  E        neutrino energy [GeV]
   * \param  l_abs    absorption length [m]
   * \return          range [m]
   */
  inline double getRange(const double E, const double l_abs)
  {
    const double precision = 1.0e-6;
    const double y0        = 1.0e-3;

    double xmin = 0.0;
    double xmax = 2.0;

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

      const double x = 0.5 * (xmin + xmax);
      const double R = 1.0e3 * pow(10.0, x);

      const double y = E * 2.0e-9 * exp(-R/l_abs) / R;

      if (fabs(y - y0) <= precision) {
        return R;
      }

      if (y > y0)
        xmin = x;
      else
        xmax = x;
    }

    return pow(10.0, 0.5 * (xmin + xmax));
  }
}

/**
 * \file
 * Auxiliary program to write signal acoustic data.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  typedef JFileRecorder<JTYPELIST<JToA, TH1D, JMeta>::typelist>  JFileRecorder_t;

  typedef JRange<double>  JRange_t;

  JFileRecorder_t  outputFile;
  string           detectorFile;
  double           T_s;
  size_t           numberOfEvents;
  double           E_GeV;
  int              run;
  int              id;
  JRange_t         X;
  JRange_t         Y;
  JRange_t         Z;

  struct {
    double         l_abs  = 5.0e3;         // absorption lentgh [m]
    double         ZMax_m = 150.0;         // signal length [m]
  } parameters;

  JSoundVelocity   V = getSoundVelocity;   // default sound velocity
  UInt_t           seed;
  int              debug;

  try {

    JProperties properties;

    properties.insert(gmake_property(parameters.l_abs));
    properties.insert(gmake_property(parameters.ZMax_m));

    JParser<> zap(" Auxiliary program to write signal acoustic data.");

    zap['a'] = make_field(detectorFile,        "detector.");
    zap['@'] = make_field(properties)                                   = JPARSER::initialised();
    zap['o'] = make_field(outputFile,          "output file");
    zap['V'] = make_field(V,                   "sound velocity")        = JPARSER::initialised();
    zap['E'] = make_field(E_GeV,               "neutrino energy [GeV]");
    zap['R'] = make_field(run,                 "run number")            = -1;
    zap['T'] = make_field(T_s,                 "time interval between events [s]");
    zap['n'] = make_field(numberOfEvents)                               = 1;
    zap['X'] = make_field(X,                   "x-range")               = JRange_t(0.0, 0.0);
    zap['Y'] = make_field(Y,                   "y-range")               = JRange_t(0.0, 0.0);
    zap['Z'] = make_field(Z,                   "z-range")               = JRange_t(0.0, 0.0);
    zap['W'] = make_field(id,                  "waveform identifier");
    zap['S'] = make_field(seed,                "seed")        = 0;
    zap['d'] = make_field(debug,               "debug")       = 0;

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


  gRandom->SetSeed(seed);

  JDetector detector;

  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }

  const double RMax_m = getRange(E_GeV, parameters.l_abs);
  
  const JCylinder3D   cylinder(detector.begin(), detector.end());
  const JPosition3D   center(cylinder.getX(), cylinder.getY(), 0.5 * (cylinder.getZmin() + cylinder.getZmax()));

  TH1D h1("h1", NULL, 60, 6.0, 12.0);

  for (Int_t i = 1; i <= h1.GetXaxis()->GetNbins(); ++i) {

    const double x = h1.GetXaxis()->GetBinCenter(i);
    const double E = pow(10.0, x);

    h1.SetBinContent(i, getRange(E, parameters.l_abs));
  }
   
  outputFile.open();

  outputFile.put(JMeta(argc, argv));
  outputFile.put(h1);

  JToA toa;

  toa.DETID                 = detector.getID();
  toa.RUN                   = run;
  toa.WAVEFORMID            = id;
  toa.QUALITYFACTOR         = 2.0e3;
  toa.QUALITYNORMALISATION  = 0.0e3;

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

    JPosition3D pos(gRandom->Uniform(X.getLowerLimit(), X.getUpperLimit()), 
                    gRandom->Uniform(Y.getLowerLimit(), Y.getUpperLimit()),
                    gRandom->Uniform(Z.getLowerLimit(), Z.getUpperLimit()));

    pos.add(center);
    
    Double_t x, y, z;

    gRandom->Sphere(x, y, z, 1.0);

    const JVersor3D dir(x,y,z);

    const JRotation3D R(dir);

    const double t0 = T_s * i;             // event time

    size_t number_of_hits = 0;
    
    for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {

      if (module->getFloor() != 0) {

        JPosition3D P(module->getPosition());

        P.sub(pos);
        P.rotate(R);

        if (sqrt(P.getX()*P.getX() + P.getY()*P.getY()) <= RMax_m && fabs(P.getZ()) <= parameters.ZMax_m) {

          const JReceiver receiver(module->getID(), module->getPosition(), module->getT0() * 1.0e-9);

          const double t1 = V.getTime(pos.getDistance(receiver), pos.getZ(), receiver.getZ());

          toa.DOMID        = receiver.getID();
          toa.TOA_NS       = llrint(1E9*receiver.putT(t0 + t1));

          outputFile.put(toa);

          number_of_hits += 1;
        }
      }
    }

    STATUS("event: " << setw(6) << i << ' '
           << FIXED(12,2) <<  t0                          << ' '
           << FIXED(9,2)  << (pos.getX() - center.getX()) << ' '
           << FIXED(9,2)  << (pos.getY() - center.getY()) << ' '
           << FIXED(9,2)  << (pos.getZ() - center.getZ()) << ' '
           << FIXED(6,3)  <<  dir.getDX()                 << ' '
           << FIXED(6,3)  <<  dir.getDY()                 << ' '
           << FIXED(6,3)  <<  dir.getDZ()                 << ' '
           << setw(4)     <<  number_of_hits << endl);    
  }  
  outputFile.close();
  
  return 0;
}