JGenerator.cc

Go to the documentation of this file.

 
#include <string>
#include <iostream>
#include <iomanip>
 
#include "TROOT.h"
#include "TFile.h"
#include "TRandom3.h"
 
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
 
#include "JAAnet/JHead.hh"
#include "JAAnet/JAAnetToolkit.hh"
 
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
 
#include "JGeometry3D/JCylinder3D.hh"
#include "JGeometry3D/JGeometry3DTestkit.hh"
 
#include "JTools/JRange.hh"
 
#include "JROOT/JRandom.hh"
 
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
namespace {
 
  const char* const volume_t   =  "volume";
  const char* const surface_t  =  "surface";
}
 
/**
 * \file
 * Basic event generator.
 *
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  typedef JRange<double> JRange_t;
 
  JFileRecorder<JTYPELIST<JAAnetTypes_t, JMeta, JRootTypes_t>::typelist> outputFile;
  counter_type    numberOfEvents;
  JCylinder3D     cylinder;
  int             type;
  JRange_t        E_GeV;
  string          option;
  JROOT::JRandom  seed;
  int             debug;
 
  try {
 
    JParser<> zap("Basic event generator.");
 
    zap['o'] = make_field(outputFile)          = "generator.root";
    zap['n'] = make_field(numberOfEvents)      = 1;
    zap['C'] = make_field(cylinder,   "generation volume");
    zap['P'] = make_field(type,       "PDG particle code");
    zap['E'] = make_field(E_GeV,      "energy range [log10(E)]");
    zap['O'] = make_field(option)              = surface_t, volume_t;
    zap['S'] = make_field(seed)                = 0;
    zap['d'] = make_field(debug)               = 2;
 
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
  seed.set(gRandom);
    
  outputFile.open();
 
  outputFile.put(JMeta(argc, argv));
  outputFile.put(*gRandom);
 
  const double A[] = {
    PI * cylinder.getRadius()*cylinder.getRadius(),
    2*PI * cylinder.getRadius()*(cylinder.getZmax() - cylinder.getZmin()),
    PI * cylinder.getRadius()*cylinder.getRadius()
  };
  
  double W = 0.0;
 
  for (int i = 0; i != sizeof(A)/sizeof(A[0]); ++i) {
    W += A[i];
  }
  
  for (counter_type counter = 0; counter != numberOfEvents; ++counter) {
 
    STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl);
 
    JVector3D    pos;
    JDirection3D dir;
 
    for ( ; ; ) {
 
      if        (option == volume_t) {
 
        pos = getRandomPosition(cylinder);
 
      } else if (option == surface_t) {
 
        const double phi = gRandom->Uniform(-PI, +PI);
        const double u   = gRandom->Uniform(0.0, W);
 
        double R;
        double z;
 
        if        (u <= A[0]) {
 
          R = sqrt(gRandom->Uniform(0.0, cylinder.getRadius()*cylinder.getRadius()));
          z = cylinder.getZmin();
 
        } else if (u <= A[1]) {
 
          R = cylinder.getRadius();
          z = gRandom->Uniform(cylinder.getZmin(), cylinder.getZmax());
        
        } else if (u <= A[2]) {
 
          R = sqrt(gRandom->Uniform(0.0, cylinder.getRadius()*cylinder.getRadius()));
          z = cylinder.getZmax();
 
        } else {
 
          continue;
        }
 
        pos = JVector3D(cylinder.getX() + R*cos(phi),
                        cylinder.getY() + R*sin(phi),
                        z);
 
      } else {
 
        FATAL("Invalid option " << option << endl);
      }
 
      randomize(&dir);
 
      if (cylinder.getIntersection(JAxis3D(pos,dir)).second > 0.0) {
        break;
      }
    }
 
    const double x = gRandom->Uniform(E_GeV.getLowerLimit(), E_GeV.getUpperLimit());
    const double E = pow(10.0, x);
 
    Evt out;
 
    out.mc_t   = 0.0;
 
    Trk trk;
 
    trk.status = TRK_ST_FINALSTATE;
    trk.type   = type;
    trk.pos    = getPosition (pos);
    trk.dir    = getDirection(dir);
    trk.E      = E;
    trk.t      = 0.0;
 
    out.mc_trks.push_back(trk);
 
    outputFile.put(out);
  }
  STATUS(endl);
 
  Head  header;
 
  {
    JHead buffer(header);
 
    buffer.fixedcan.xcenter  = cylinder.getX();
    buffer.fixedcan.ycenter  = cylinder.getY();
    buffer.fixedcan.radius   = cylinder.getRadius();
    buffer.fixedcan.zmin     = cylinder.getZmin();
    buffer.fixedcan.zmax     = cylinder.getZmax();
  
    buffer.push(&JHead::fixedcan);
 
    copy(buffer, header);
  }
 
  outputFile.put(header);
  outputFile.put(*gRandom);
 
  outputFile.close();
}