CRYGen.cxx

Go to the documentation of this file.

#include "EventGenerator/Cosmics/CRYGen.h"
#include <cmath>
#include <iostream>
// CRY include files
#include "CRYSetup.h"
#include "CRYParticle.h"
#include "CRYGenerator.h"
// ROOT include files
#include "TRandom.h"
// NOvA include files
#include "EventGenerator/evgen.h"
#include "Geometry/geo.h"
#include "RawData/DAQHeader.h"
#include "Simulation/MCTruth.h"
using namespace evgen;

//......................................................................

static double gsUniformRandom() 
{
  return gRandom->Uniform();
}

//......................................................................

CRYGen::CRYGen(int idet, int runno) :
  fGeom(geo::Geometry::Instance(runno, idet))
{
  // Apply an energy threshold to avoid extremely long lists of
  // irrelevant particles
  fEthresh = 50.0E-3; // GeV

  // By default, sample for 30 usec, 10 usec before spill, 20 usec
  // after
  fSampleTime =  30.0E-6;
  fToffset    = -10.0E-6;

  // Construct the CRY generator
  std::string config("date 1-1-2014 "
                     "returnGammas 1 "
                     "returnElectrons 1 "
                     "returnMuons 1 "
                     "returnPions 1 "
                     "returnNeutrons 1 "
                     "returnProtons 1 ");
  
  if (idet == rawdata::kIPND) {
    config += 
      "latitude 41.8 "
      "altitude 0 "
      "subboxLength 50 ";
  }
  if (idet == rawdata::kNear)  {
    config += 
      "latitude 41.8 "
      "altitude 0 "
      "subboxLength 300 ";
  }
  if (idet == rawdata::kFar) {
    config += 
      "latitude 48.3 "
      "altitude 0 "
      "subboxLength 300 ";
  }

  // Find the pointer to the CRY data tables
  std::string crydatadir(getenv("CRYDATADIR"));

  // Construct the event generator object
  fSetup = new CRYSetup(config, crydatadir);
  fSetup->setRandomFunction(gsUniformRandom);
  
  fGen = new CRYGenerator(fSetup);
  fEvt = new std::vector<CRYParticle*>;
}

//......................................................................

void CRYGen::Sample(sim::MCTruth& mctruth, double* w)
{
  double tstart;        // Generator time at start of sample
  TDatabasePDG pdgdata; // Database of PDG information

  tstart = fGen->timeSimulated();
  while (1) {
    fEvt->clear();
    fGen->genEvent(fEvt);
    
    mctruth.Clear();
    for (unsigned int i=0; i<fEvt->size(); ++i) {
      CRYParticle* cryp = (*fEvt)[i];
      
      // Pull out the PDG code
      int pdg = cryp->PDGid();

      // Get the energies of the particles
      double ke = cryp->ke()*1.0E-3; // MeV to GeV conversion
      if (ke<fEthresh) continue;

      double m    = pdgdata.GetParticle(pdg)->Mass(); // In GeV
      double etot = ke + m;
      double ptot = etot*etot-m*m;
      if (ptot>0.0) ptot = sqrt(ptot);
      else          ptot = 0.0;
      
      // Sort out the momentum components. Remember that the NOvA
      // frame has y up and z along the beam. So uvw -> zxy
      double px = ptot * cryp->v();
      double py = ptot * cryp->w();
      double pz = ptot * cryp->u();
      
      // Particle start position. CRY distributes uniformly in x-y
      // plane at fixed z, where z is the vertical direction. This
      // requires some offsets and rotations to put the particles at
      // the surface in the NOvA geometry as well as some rotations
      // since the NOvA coordinate frame has y up and z along the
      // beam.
      double vx = cryp->y()*100.0;
      double vy = cryp->z()*100.0 + fGeom.SurfaceY();
      double vz = cryp->x()*100.0 + 0.5*fGeom.DetLength();
      double t  = cryp->t()-tstart + fToffset; // seconds

      // Project backward to edge of world volume
      double xyz[3]  = { vx,  vy,  vz};
      double xyzo[3];
      double dxyz[3] = {-px, -py, -pz};
      double x1, x2;
      double y1, y2;
      double z1, z2;
      fGeom.WorldBox(&x1, &x2, &y1, &y2, &z1, &z2);
      geo::ProjectToBoxEdge(xyz, dxyz, x1, x2, y1, y2, z1, z2, xyzo);
      vx = xyzo[0];
      vy = xyzo[1];
      vz = xyzo[2];
      
      // Boiler plate...
      int istatus    =  1;
      int imother1   = kCosmicRayGenerator;
      int imother2   = kCosmicRayGenerator;
      int idaughter1 = -1;
      int idaughter2 = -1;
      
      // Push the particle onto the stack
      TParticle p(pdg,
                  istatus,
                  imother1, imother2,
                  idaughter1, idaughter2,
                  px, py, pz,
                  etot,
                  vx, vy, vz, t);
      mctruth.Add(p);
    } // Loop on particles in event
    
    // Check if we're done with this time sample
    if (fGen->timeSimulated()-tstart > fSampleTime) break;
    
  } // Loop on events simulated

  // Check if this time slice passes selection criteria
  // TODO
  if (w) *w = 1.0;
}

---