SupernovaGen_module.cc

Go to the documentation of this file.

///
/// \brief   SNOwGLoBES-based Event Generator for Supernovae
/// \author  Justin Vasel <jvasel@indiana.edu>
/// \author  Matt Tamsett <m.tamsett@sussex.ac.uk>
///

// C++ includes
#include <cmath>
#include <iostream>
#include <memory>
#include <string>
#include <vector>

// ROOT includes
#include "TDatabasePDG.h"
#include "TFile.h"
#include "TH2F.h"
#include "TRandom3.h"

// CLHEP includes
#include "CLHEP/Random/RandFlat.h"
#include "CLHEP/Random/RandGaussQ.h"

// Framework includes
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/SubRun.h"
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Framework/Services/Optional/RandomNumberGenerator.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// NOvASoft includes
#include "nugen/EventGeneratorBase/evgenbase.h"
#include "Geometry/Geometry.h"
#include "NovaDAQConventions/DAQConventions.h"
#include "nusimdata/SimulationBase/MCParticle.h"
#include "nusimdata/SimulationBase/MCTruth.h"
#include "SummaryData/RunData.h"
#include "SummaryData/SubRunData.h"


namespace evgen {
  class SupernovaGen : public art::EDProducer {
  public:
    explicit SupernovaGen(fhicl::ParameterSet const& p);
    virtual ~SupernovaGen();                        

    void produce(art::Event& e);  
    void reconfigure(const fhicl::ParameterSet& p);
    void beginRun(art::Run& r);
    void endSubRun(art::SubRun& sr);

  private:

    void Generate(std::unique_ptr< std::vector<simb::MCTruth>> &truthcol);

    /// Depending on configuration of fPMeaning
    /// calculate the momentum of the particle
    double GetMomentum(const double& pkf, const double& m) const;
    
    static const int    kUNIF = 0;
    static const int    kGAUS = 1;
    
    int                 fPDG;                ///< PDG code of particles to generate
    int                 fPMeaning;           ///< Meaning of P0, fSigmaP and fPDist. 
                                             ///< By default (fP0Meaning=0), P0 and sigmaP is momentum
                                             /// If fPMeaning=1, then P0 and sigmaP is kinetic energy in GeV
                                             /// If fPMeaning=2, then P0and sigmaP is total energy in GeV
    int                 fCycle;              ///< MC production cycle
    double              fXlow;               ///< Lower x position (cm)
    double              fXhigh;              ///< Upper x position (cm)
    double              fYlow;               ///< Lower y position (cm)
    double              fYhigh;              ///< Upper y position (cm)
    double              fZlow;               ///< Lower z position (cm)
    double              fZhigh;              ///< Upper z position (cm)
    double              fTlow;               ///< Lower t position (ns)
    double              fThigh;              ///< Upper t position (ns)
    double              fCosZ0;              ///< Cosine of central angle wrt z-axis
    double              fPhiXY0;             ///< Central angle in the x-y plane (degrees)
    double              fSigmaCosZ;          ///< Size of variation of cosz
    double              fSigmaPhiXY;         ///< Size of variation in phixy (degrees)
    int                 fAngleDist;          ///< How to distribute angles (0=uniform, 1=gaussian)
    std::string         fPositronPDFFile;    ///< ROOT file name where the positron PDFs are stored
    bool                fOnlyNoise;          ///< Generate no signal
                                             // Multiplicative scaling factor to change the distance of the SN 
                                             // the given rates are proportional to the distance^2 (10^2)
    double              fDistanceScaling;    ///< Scaling factor to PDF values
    art::EventNumber_t  fEventPadding;       ///< Number of empty events to place at beginning of file
    
    // Histograms
    TH2F* fPositronPDF;
    
    int fNtimePoints, fNenergyBins;
    double fMinPositronEnergy = 0.0007488;   // GeV
    double fMaxPositronEnergy = 0.09975;     // GeV

    double fSigmaP;
    int fTimeIndex = 0;
  };
}

////////////////////////////////////////////////////////////////////////

namespace evgen{

  //----------------------------------------------------------------------------
  SupernovaGen::SupernovaGen(fhicl::ParameterSet const& p)
  : fPMeaning(0)
  , fCycle   (p.get<int>("Cycle", 0))
  {
    this->reconfigure(p);

    // get the random number seed, use a random default if not specified    
    // in the configuration file.  
    int seed = p.get<unsigned int>("Seed", evgb::GetRandomNumberSeed());
  
    createEngine( seed );
    
    produces<std::vector<simb::MCTruth>>();
    produces<sumdata::SubRunData, art::InSubRun>();
    produces<sumdata::RunData,    art::InRun   >();
  }

  //----------------------------------------------------------------------------
  SupernovaGen::~SupernovaGen() { }

  //----------------------------------------------------------------------------
  void SupernovaGen::reconfigure(const fhicl::ParameterSet& p)
  {
    fPDG = p.get<int>("PDG");
    try{
      fPMeaning = p.get<int>("PMeaning");
      if (fPMeaning == 1) {
        std::cout << "SupernovaGen: Using Kinetic energy for the meaning of P0, SigmaP and PDist" << std::endl;
      }
      else if(fPMeaning == 2) {
        std::cout << "SupernovaGen: Using Total energy for the meaning of P0, SigmaP and PDist" << std::endl;
      }
    }
    catch(...){
      fPMeaning = 0;
    }
    
    fXlow               = p.get<double>      ("Xlow");
    fYlow               = p.get<double>      ("Ylow");
    fZlow               = p.get<double>      ("Zlow");
    fXhigh              = p.get<double>      ("Xhigh");
    fYhigh              = p.get<double>      ("Yhigh");
    fZhigh              = p.get<double>      ("Zhigh");
    fTlow               = p.get<double>      ("Tlow");
    fThigh              = p.get<double>      ("Thigh");
    fCosZ0              = p.get<double>      ("CosZ0");
    fPhiXY0             = p.get<double>      ("PhiXY0");
    fSigmaPhiXY         = p.get<double>      ("SigmaPhiXY");
    fSigmaCosZ          = p.get<double>      ("SigmaCosZ");
    fAngleDist          = p.get<int>         ("AngleDist");
    fPositronPDFFile    = p.get<std::string> ("PositronPDFFile");
    fOnlyNoise          = p.get<bool>        ("OnlyNoise");
    fDistanceScaling    = p.get<double>      ("DistanceScaling");
    
    fEventPadding       = p.get<art::EventNumber_t>("EventPadding");
    
    TFile *root_file = new TFile(fPositronPDFFile.c_str());
    root_file->GetObject("ibd", fPositronPDF);
    
    fNtimePoints = fPositronPDF->GetNbinsX();
    fNenergyBins = fPositronPDF->GetNbinsY();
    
    fSigmaP = (fMaxPositronEnergy - fMinPositronEnergy) / (2. * double(fNenergyBins));
  }

  //----------------------------------------------------------------------------
  void SupernovaGen::beginRun(art::Run& r)
  {
    // Get the detector information and push it into the event as RunData
    art::ServiceHandle<geo::Geometry> geo;
    std::unique_ptr<sumdata::RunData> runcol(new sumdata::RunData(geo->DetId(),
                                                                  geo->FileBaseName(),
                                                                  geo->ExtractGDML()));
    r.put(std::move(runcol));
    return;
  }
  
  //----------------------------------------------------------------------------
  void SupernovaGen::endSubRun(art::SubRun& sr)
  {
    // store the cycle information
    art::ServiceHandle<geo::Geometry> geo;
    std::unique_ptr< sumdata::SubRunData > sd(new sumdata::SubRunData(fCycle));
    sr.put(std::move(sd));
  }

  //----------------------------------------------------------------------------
  void SupernovaGen::produce(art::Event& e)
  {    
    std::unique_ptr<std::vector<simb::MCTruth>> truthcol(new std::vector<simb::MCTruth>);
        
    // Only generate particles if we're done padding with empty events
    if (e.event() > fEventPadding){
      fTimeIndex++;
      
      // check that time index is within the available bounds
      // TODO exit gracefully rather than crashing
      if (fTimeIndex > fNtimePoints) return;
      
      Generate(truthcol);
    }
      
    e.put(std::move(truthcol));
    return;
  }

  //----------------------------------------------------------------------------
  void SupernovaGen::Generate(std::unique_ptr<std::vector<simb::MCTruth>> &truthcol) 
  {
    // Set up the generators
    art::ServiceHandle<art::RandomNumberGenerator> rng;
    
    CLHEP::HepRandomEngine &engine = rng->getEngine();
    CLHEP::RandFlat         flat(engine);
    CLHEP::RandGaussQ       gauss(engine);
    
    //
    // every event will have the positron distribution from the fPositronPDF array
    //
    int generatedPositrons = 0;

    for (int energyBin = 0; energyBin < fNenergyBins; ++energyBin) {
      // the central momentum is given by the bin centre
      double P0 = fMinPositronEnergy + (double(energyBin + 1) * 2. * fSigmaP);
      double positronPDF = 0.;
      
      if (!fOnlyNoise) 
        positronPDF = fPositronPDF->GetBinContent(fTimeIndex, energyBin + 1);
        
      TRandom3 rnd(0); // zero seeds randomly
        
      int nPositrons = rnd.Poisson(positronPDF * fDistanceScaling);
      
      // Silence this output if no positrons produced
      if (nPositrons > 0) {
        std::cout << "\tEnergy bin["         << energyBin 
                  << "] ("                   << P0
                  << " GeV), positron PDF: " << positronPDF
                  << std::endl;
                  
        // Generate a Poisson-fluctated number of positrons based on this
        std::cout << "\tGenerating "         << nPositrons
                  << " positrons"            << std::endl; 
      }
      
      for (int positronIdx = 0; positronIdx < nPositrons; ++positronIdx) {
        generatedPositrons++;
        simb::MCTruth mct;
        mct.SetOrigin(simb::kSingleParticle);
      
        // Momentum, kinetic energy or full energy
        double pkf = 0.0;
        pkf = flat.fire(P0 - fSigmaP, P0 + fSigmaP);
        
        const TDatabasePDG* pdgt = TDatabasePDG::Instance();
        const TParticlePDG* pdgp = pdgt->GetParticle(fPDG);
        
        // Mass in GeV
        double    m = 0.0;
        if (pdgp) m = pdgp->Mass();
  
        // Momentum in GeV/c
        const double p = GetMomentum(pkf, m);
  
        // Choose position
        double x[4];
        x[0] = flat.fire(fXlow, fXhigh);
        x[1] = flat.fire(fYlow, fYhigh);
        x[2] = flat.fire(fZlow, fZhigh);
        x[3] = flat.fire(fTlow, fThigh);
        
        std::cout << "\t\tpositron[ " << generatedPositrons-1
                  << "]: p: "         << p
                  << " GeV, x: "      << x[0]
                  << " cm, y: "       << x[1]
                  << " cm, z: "       << x[2]
                  << " cm, t: "       << x[3]
                  << " ns"            << std::endl;
                  
        const TLorentzVector pos(x[0], x[1], x[2], x[3]);
        
        // Choose angles
        double cosz, phi;
        unsigned int tryIdx;
        for (tryIdx = 0; tryIdx < 1000000; ++tryIdx) {
          if (fAngleDist == kGAUS) {
            cosz = gauss.fire(fCosZ0,  fSigmaCosZ);
            phi  = gauss.fire(fPhiXY0, fSigmaPhiXY);
          } else {
            cosz = flat.fire(fCosZ0 - fSigmaCosZ,
                             fCosZ0 + fSigmaCosZ);
      
            phi  = flat.fire(fPhiXY0 - fSigmaPhiXY,
                             fPhiXY0 + fSigmaPhiXY);
          }
          if (cosz >= -1.0 && cosz <= 1.0) break;
        }
        if (cosz < -1.0 || cosz > 1.0) {
          mf::LogError("SupernovaGen") << __FILE__ << ":" << __LINE__
                        << " Failed to find an acceptable cos(theta_z)"
                        << " after many tries.\n"
                        << " Please adjust CosZ0 and SigmaCosZ0 in your"
                        << " SupernovaGen.fcl file and rerun";
          abort();
        }
        
        const double sinz    = sqrt(1.0 - cosz * cosz);
        const double sinphi  = sin(phi * M_PI / 180.0);
        const double cosphi  = cos(phi * M_PI / 180.0);
        
        // Set track id to -i as these are all primary particles and have id <= 0
        const int trackid = -1 * (generatedPositrons + 1);
        
        const std::string primary("primary");
        simb::MCParticle part(trackid, fPDG, primary);
        
        const TLorentzVector pvec(cosphi * sinz * p,
                                  sinphi * sinz * p,
                                  cosz * p,
                                  sqrt(p * p + m * m));
  
        part.AddTrajectoryPoint(pos, pvec);
      
        mct.Add(part);
        truthcol->push_back(mct);
      } // end of loop over positrons
    } // end loop over energy bins particles
    
    std::cout << "Generated a total of " << generatedPositrons
              << " positrons for event " << fTimeIndex << std::endl;
  }

  //------------------------------------------------------------------------------
  double SupernovaGen::GetMomentum(const double& pkf, const double& m) const{

      if (fPMeaning == 0) return pkf;

      double total_energy = 0.0;
      if (fPMeaning == 1){
        total_energy = pkf + m;
      }
      else if(fPMeaning == 2){
        total_energy = pkf;
      }
      else{
        total_energy = sqrt(pkf * pkf + m * m);
      }

      return sqrt(total_energy * total_energy - m * m);
  }

  DEFINE_ART_MODULE(SupernovaGen)

} //end namespace