evgen::HepMCFileGen Class Reference

Inheritance diagram for evgen::HepMCFileGen:

Public Member Functions
	HepMCFileGen (fhicl::ParameterSet const &p)
void	produce (art::Event &e) override
void	beginJob () override
void	beginRun (art::Run &run) override
void	endSubRun (art::SubRun &sr) override

Private Member Functions
std::ifstream	open_file ()

Private Attributes
std::string	fInputFilePath
	Path to the HEPMC input file, relative to FW_SEARCH_PATH. More...
std::ifstream *	fInputFile
double	fEventsPerPOT
	Number of events per POT (to be set) More...
int	fEventsPerSubRun
	Keeps track of the number of processed events per subrun. More...

Detailed Description

This module assumes that the input file has the hepevt format for each event to be simulated. In brief each event contains at least two lines. The first line contains two entries, the event number (which is ignored in ART/LArSoft) and the number of particles in the event. Each following line containes 15 entries to describe each particle. The entries are:

1. status code (should be set to 1 for any particle to be tracked, others won't be tracked)
2. the pdg code for the particle
3. the entry of the first mother for this particle in the event, 0 means no mother
4. the entry of the second mother for this particle in the event, 0 means no mother
5. the entry of the first daughter for this particle in the event, 0 means no daughter
6. the entry of the second daughter for this particle in the event, 0 means no daughter
7. x component of the particle momentum
8. y component of the particle momentum
9. z component of the particle momentum
10. energy of the particle
11. mass of the particle
12. x position of the particle initial position
13. y position of the particle initial position
14. z position of the particle initial position
15. time of the particle production

For example, if you want to simulate a single muon with a 5 GeV energy moving only in the z direction, the entry would be (see onemuon.hepmc):

0 1
1 13 0 0 0 0 0. 0. 1.0 5.0011 0.105 1.0 1.0 1.0 0.0

Or if you want to simulate a muon neutrino event (see oneneutrino.hepmc):

0 3
0 14 0 0 0 0 0.00350383 0.002469 0.589751 0.589766 0 208.939 63.9671 10.9272 4026.32
1 13 1 0 0 0 -0.168856 -0.0498011 0.44465 0.489765 105.658 208.939 63.9671 10.9272 4026.32
1 2212 1 0 0 0 0.151902 -0.124578 0.0497377 0.959907 938.272 208.939 63.9671 10.9272 4026.32

here the first particle is the initial state neutrino (status code 0, meaning initial state, not to be prooagated by GEANT). The other particles are the initial state particles. For the neutrino, write ccnc in place of 1st daugther, and mode (qe, res, ...) in place of 2nd daugther.

There are some assumptions that go into using this format that may not be obvious. The first is that only particles with status code = 1 are tracked in the LArSoft/Geant4 combination making the mother daughter relations somewhat irrelevant. That also means that you should let Geant4 handle any decays.

The units in LArSoft are cm for distances and ns for time. The use of TLorentzVector below does not imply space and time have the same units (do not use TLorentzVector::Boost()).

Definition at line 101 of file HepMCFileGen_module.cc.

Constructor & Destructor Documentation

evgen::HepMCFileGen::HepMCFileGen ( fhicl::ParameterSet const &  p )

explicit

Definition at line 118 of file HepMCFileGen_module.cc.

  : EDProducer{p}
  , fInputFilePath(p.get<std::string>("InputFilePath"))
  , fInputFile(nullptr)
  , fEventsPerPOT{p.get<double>("EventsPerPOT", -1.)}
  , fEventsPerSubRun(0)
{
  produces< std::vector<simb::MCTruth>   >();
  produces< sumdata::RunData, art::InRun >();
  produces< sumdata::POTSummary, art::InSubRun >();
}

Member Function Documentation

void evgen::HepMCFileGen::beginJob ( )

override

Definition at line 178 of file HepMCFileGen_module.cc.

{
  fInputFile = new std::ifstream(open_file());
  
}

void evgen::HepMCFileGen::beginRun ( art::Run &  run )

override

Definition at line 184 of file HepMCFileGen_module.cc.

{
    fEventsPerSubRun = 0;
    art::ServiceHandle<geo::Geometry const> geo;
    run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()));
  }

void evgen::HepMCFileGen::endSubRun ( art::SubRun &  sr )

override

Definition at line 191 of file HepMCFileGen_module.cc.

  {
    auto p = std::make_unique<sumdata::POTSummary>();
    p->totpot     = fEventsPerSubRun * fEventsPerPOT;
    p->totgoodpot = fEventsPerSubRun * fEventsPerPOT;
    sr.put(std::move(p));
    return;
  }

std::ifstream evgen::HepMCFileGen::open_file ( )

private

Definition at line 131 of file HepMCFileGen_module.cc.

{
  /*
   * The plan:
   *  1. expand the path in FW_SEARCH_PATH (only if relative path)
   *  2. copy it into scratch area (only if starts with `/pnfs`)
   *  3. open the file (original or copy) and return the opened file
   * 
   * Throws a cet::exception if eventually file is not found.
   */
  
  std::string fullFileName = fInputFilePath;
  
  cet::search_path searchPath("FW_SEARCH_PATH");
  if (searchPath.find_file(fInputFilePath, fullFileName)) {
    mf::LogDebug("HepMCFileGen")
      << "Input file '" << fInputFilePath << "' found in FW_SEARCH_PATH:\n"
      << fullFileName
      ;
  }
  
  //
  // prepare the file with IFDH, if path starts with `/pnfs`
  //
  if (fullFileName.compare(0, 6, "/pnfs/") == 0) { 
    fullFileName = art::ServiceHandle<IFDH>()->fetchInput(fullFileName);
    mf::LogDebug("HepMCFileGen")
      << "IFDH fetch: '" << fInputFilePath << "' -> '" << fullFileName << "'";
  }
  
  //
  // attempt to open
  //
  mf::LogDebug("HepMCFileGen")
    << "Reading input file '" << fInputFilePath << "' as:\n" << fullFileName;
  std::ifstream inputFile(fullFileName);
  if (inputFile) return inputFile;
  
  // all attempts failed, give up:
  throw cet::exception("HepMCFileGen")
    << "HEPMC input file '" << fInputFilePath << "' can't be opened.\n";
  
} // evgen::HepMCFileGen::open_file()

void evgen::HepMCFileGen::produce ( art::Event &  e )

override

Definition at line 200 of file HepMCFileGen_module.cc.

{
  // check that the file is still good
  if( !fInputFile->good() || fInputFile->peek() == EOF) {
    open_file();
  }
  if( !fInputFile->good() || fInputFile->peek() == EOF) {
    throw cet::exception("HepMCFileGen") << "input text file "
                                        << " cannot be read in produce().\n";
  }
  std::unique_ptr< std::vector<simb::MCTruth> > truthcol(new std::vector<simb::MCTruth>);
  simb::MCTruth truth;
  // declare the variables for reading in the event record
  int            event          = 0;
  unsigned short nParticles           = 0;
  int            status         = 0;
  int                   pdg            = 0;
  int                   firstMother    = 0;
  int                   secondMother   = 0;
  int                   firstDaughter  = 0;
  int                   secondDaughter = 0;
  double          xMomentum      = 0.;
  double          yMomentum           = 0.;
  double          zMomentum           = 0.;
  double          energy              = 0.;
  double          mass                = 0.;
  double          xPosition           = 0.;
  double          yPosition           = 0.;
  double          zPosition           = 0.;
  double          time                = 0.;
  bool set_neutrino = false;
  // neutrino
  int ccnc = -1, mode = -1, itype = -1, target = -1, nucleon = -1, quark = -1;
  double w = -1, x = -1, y = -1, qsqr = -1;
  // read in line to get event number and number of particles
  std::string oneLine;
  std::getline(*fInputFile, oneLine);
  std::istringstream inputLine;
  inputLine.str(oneLine);
  inputLine >> event >> nParticles;
  // now read in all the lines for the particles
  // in this interaction. only particles with
  // status = 1 get tracked in Geant4. see GENIE GHepStatus
  for(unsigned short i = 0; i < nParticles; ++i){
    std::getline(*fInputFile, oneLine);
    inputLine.clear();
    inputLine.str(oneLine);
    inputLine >> status >> pdg
              >> firstMother >> secondMother >> firstDaughter >> secondDaughter
              >> xMomentum   >> yMomentum    >> zMomentum     >> energy >> mass
              >> xPosition   >> yPosition    >> zPosition     >> time;
    std::cout<<"pdg_all"<<pdg<<"xmomentum_all"<<xMomentum<<"ymomentum_all"<<yMomentum<<"zposition_all"<<zPosition<<std::endl;
    TLorentzVector pos(xPosition, yPosition, zPosition, time);
    TLorentzVector mom(xMomentum, yMomentum, zMomentum, energy);
    simb::MCParticle part(i, pdg, "primary", firstMother, mass, status);
    part.AddTrajectoryPoint(pos, mom);
    //if (abs(pdg) == 18 || abs(pdg) == 12) 
    if (abs(pdg) == 52 )  // Animesh made changes
        {
      set_neutrino = true;
      ccnc = firstDaughter; // for the neutrino we write ccnc in place of 1st daugther
      mode = secondDaughter; // for the neutrino we write mode in place of 2nd daugther
      itype = -1;
      target = nucleon = quark = w = x = y = qsqr = -1;
     std::cout<<"pdg_mother"<<pdg<<"xmomentum_mother"<<xMomentum<<"ymomentum_mother"<<yMomentum<<"zposition_mother"<<zPosition<<std::endl;
    } 
    truth.Add(part);
     std::cout << i << "  Particle added with Pdg " << part.PdgCode() << ", Mother " << part.Mother() << ", track id " << part.TrackId() << ", ene " << part.E() << std::endl;
  }
 
  if (set_neutrino) {
    truth.SetNeutrino(ccnc,
                      mode,
                      itype,
                      target,
                      nucleon,
                      quark,
                      w,
                      x,
                      y,
                      qsqr);
    // set the neutrino information in MCTruth
    truth.SetOrigin(simb::kBeamNeutrino);
    // truth.SetGeneratorInfo(simb::Generator_t::kGENIE,
    //                          __GENIE_RELEASE__,
    //                          {{"tune", fTuneName}});
  }
  //std::cout << " neutrino " << truth.GetNeutrino() << std::endl;
  //std::cout << " lepton pdg " << truth.GetNeutrino().Lepton().PdgCode() << " ene " << truth.GetNeutrino().Lepton().E() << std::endl;
  truthcol->push_back(truth);
  e.put(std::move(truthcol));
  fEventsPerSubRun++;
  return;
}

Member Data Documentation

double evgen::HepMCFileGen::fEventsPerPOT

private

Number of events per POT (to be set)

Definition at line 113 of file HepMCFileGen_module.cc.

int evgen::HepMCFileGen::fEventsPerSubRun

private

Keeps track of the number of processed events per subrun.

Definition at line 114 of file HepMCFileGen_module.cc.

std::ifstream* evgen::HepMCFileGen::fInputFile

private

Definition at line 111 of file HepMCFileGen_module.cc.

std::string evgen::HepMCFileGen::fInputFilePath

private

Path to the HEPMC input file, relative to FW_SEARCH_PATH.

Definition at line 110 of file HepMCFileGen_module.cc.

---

The documentation for this class was generated from the following file:

• HepMCFileGen_module.cc