NuclearInteraction.cc

Go to the documentation of this file.

// system headers
#include <cmath>
#include <memory>
#include <vector>
#include <map>
#include <string>
#include <fstream>
#include <iostream>
#include <sys/stat.h>
#include <cmath>
#include <TFile.h>
#include <TTree.h>
#include <TROOT.h>
#include <Math/AxisAngle.h>
#include "Math/Boost.h"
#include "DataFormats/Math/interface/Vector3D.h"

// Framework Headers
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "FWCore/ParameterSet/interface/FileInPath.h"

#include "FWCore/Utilities/interface/thread_safety_macros.h"

// Fast Sim headers
#include "FastSimulation/SimplifiedGeometryPropagator/interface/Particle.h"
#include "FastSimulation/SimplifiedGeometryPropagator/interface/SimplifiedGeometry.h"
#include "FastSimulation/SimplifiedGeometryPropagator/interface/InteractionModelFactory.h"
#include "FastSimulation/SimplifiedGeometryPropagator/interface/InteractionModel.h"
#include "FastSimulation/SimplifiedGeometryPropagator/interface/Constants.h"
#include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
#include "FastSimDataFormats/NuclearInteractions/interface/NUEvent.h"

// Math
#include "DataFormats/Math/interface/LorentzVector.h"

// Author: Patrick Janot
// Date: 25-Jan-2007
//
// Revision: Class structure modified to match SimplifiedGeometryPropagator
//           Unified the way, the closest charged daughter is defined
//           S. Kurz, 29 May 2017

// TODO: save() function called in destructer. Does that actually make sense?

typedef math::XYZVector XYZVector;
typedef math::XYZTLorentzVector XYZTLorentzVector;

namespace fastsim {

  class NuclearInteraction : public InteractionModel {
  public:
    NuclearInteraction(const std::string& name, const edm::ParameterSet& cfg);

    ~NuclearInteraction() override;


    void interact(fastsim::Particle& particle,
                  const SimplifiedGeometry& layer,
                  std::vector<std::unique_ptr<fastsim::Particle> >& secondaries,
                  const RandomEngineAndDistribution& random) override;

  private:
    unsigned index(int thePid);

    XYZVector orthogonal(const XYZVector& aVector) const;

    void save();

    bool read(std::string inputFile);

    double theDistCut;       
    double theHadronEnergy;  
    std::string inputFile;   

    // Read/Save nuclear interactions from FullSim

    TFile* theFile = nullptr;                                     
    std::vector<std::vector<TTree*> > theTrees;                   
    std::vector<std::vector<TBranch*> > theBranches;              
    std::vector<std::vector<NUEvent*> > theNUEvents;              
    std::vector<std::vector<unsigned> > theCurrentEntry;          
    std::vector<std::vector<unsigned> > theCurrentInteraction;    
    std::vector<std::vector<unsigned> > theNumberOfEntries;       
    std::vector<std::vector<unsigned> > theNumberOfInteractions;  
    std::vector<std::vector<std::string> > theFileNames;          
    std::vector<std::vector<double> > theHadronCM;                

    unsigned ien4;  

    std::ofstream myOutputFile;  
    unsigned myOutputBuffer;     

    bool currentValuesWereSet;  

    // Properties of the Hadrons

    static std::vector<std::vector<double> > theRatiosMap;
    static std::map<int, int> theIDMap;

    const std::vector<double> theHadronEN = {
        1.0, 2.0, 3.0, 4.0, 5.0, 7.0, 9.0, 12.0, 15.0, 20.0, 30.0, 50.0, 100.0, 200.0, 300.0, 500.0, 700.0, 1000.0};

    // Use same order for all vectors below (Properties for "piplus", "piminus", "K0L", "Kplus", "Kminus", "p", "pbar", "n", "nbar")

    const std::vector<int> theHadronID = {211, -211, 130, 321, -321, 2212, -2212, 2112, -2112};
    const std::vector<std::string> theHadronNA = {
        "piplus", "piminus", "K0L", "Kplus", "Kminus", "p", "pbar", "n", "nbar"};
    const std::vector<double> theHadronMA = {
        0.13957, 0.13957, 0.497648, 0.493677, 0.493677, 0.93827, 0.93827, 0.939565, 0.939565};
    const std::vector<double> theHadronPMin = {0.7, 0.0, 1.0, 1.0, 0.0, 1.1, 0.0, 1.1, 0.0};
    const std::vector<double> theLengthRatio = {//pi+      pi-    K0L      K+      K-      p      pbar     n      nbar
                                                0.2257,
                                                0.2294,
                                                0.3042,
                                                0.2591,
                                                0.2854,
                                                0.3101,
                                                0.5216,
                                                0.3668,
                                                0.4898};
    const std::vector<double> theRatios = {// pi+ (211)
                                           0.031390573,
                                           0.531842852,
                                           0.819614219,
                                           0.951251711,
                                           0.986382750,
                                           1.000000000,
                                           0.985087033,
                                           0.982996773,
                                           0.990832192,
                                           0.992237923,
                                           0.994841580,
                                           0.973816742,
                                           0.967264815,
                                           0.971714258,
                                           0.969122824,
                                           0.978681792,
                                           0.977312732,
                                           0.984255819,
                                           // pi- (-211)
                                           0.035326512,
                                           0.577356403,
                                           0.857118809,
                                           0.965683504,
                                           0.989659360,
                                           1.000000000,
                                           0.989599240,
                                           0.980665408,
                                           0.988384816,
                                           0.981038152,
                                           0.975002104,
                                           0.959996152,
                                           0.953310808,
                                           0.954705592,
                                           0.957615400,
                                           0.961150456,
                                           0.965022184,
                                           0.960573304,
                                           // K0L (130)
                                           0.000000000,
                                           0.370261189,
                                           0.649793096,
                                           0.734342408,
                                           0.749079499,
                                           0.753360057,
                                           0.755790543,
                                           0.755872164,
                                           0.751337674,
                                           0.746685288,
                                           0.747519634,
                                           0.739357554,
                                           0.735004444,
                                           0.803039922,
                                           0.832749896,
                                           0.890900187,
                                           0.936734805,
                                           1.000000000,
                                           // K+ (321)
                                           0.000000000,
                                           0.175571717,
                                           0.391683394,
                                           0.528946472,
                                           0.572818635,
                                           0.614210280,
                                           0.644125538,
                                           0.670304050,
                                           0.685144573,
                                           0.702870161,
                                           0.714708513,
                                           0.730805263,
                                           0.777711536,
                                           0.831090576,
                                           0.869267129,
                                           0.915747562,
                                           0.953370523,
                                           1.000000000,
                                           // K- (-321)
                                           0.000000000,
                                           0.365353210,
                                           0.611663677,
                                           0.715315908,
                                           0.733498956,
                                           0.738361302,
                                           0.745253654,
                                           0.751459671,
                                           0.750628335,
                                           0.746442657,
                                           0.750850669,
                                           0.744895986,
                                           0.735093960,
                                           0.791663444,
                                           0.828609543,
                                           0.889993040,
                                           0.940897842,
                                           1.000000000,
                                           // proton (2212)
                                           0.000000000,
                                           0.042849136,
                                           0.459103223,
                                           0.666165343,
                                           0.787930873,
                                           0.890397011,
                                           0.920999533,
                                           0.937832788,
                                           0.950920131,
                                           0.966595049,
                                           0.979542270,
                                           0.988061653,
                                           0.983260159,
                                           0.988958431,
                                           0.991723494,
                                           0.995273237,
                                           1.000000000,
                                           0.999962634,
                                           // anti-proton (-2212)
                                           1.000000000,
                                           0.849956907,
                                           0.775625988,
                                           0.802018230,
                                           0.816207485,
                                           0.785899785,
                                           0.754998487,
                                           0.728977244,
                                           0.710010673,
                                           0.670890339,
                                           0.665627872,
                                           0.652682888,
                                           0.613334247,
                                           0.647534574,
                                           0.667910938,
                                           0.689919693,
                                           0.709200185,
                                           0.724199928,
                                           // neutron (2112)
                                           0.000000000,
                                           0.059216484,
                                           0.437844536,
                                           0.610370629,
                                           0.702090648,
                                           0.780076890,
                                           0.802143073,
                                           0.819570432,
                                           0.825829666,
                                           0.840079750,
                                           0.838435509,
                                           0.837529986,
                                           0.835687165,
                                           0.885205014,
                                           0.912450156,
                                           0.951451221,
                                           0.973215562,
                                           1.000000000,
                                           // anti-neutron
                                           1.000000000,
                                           0.849573257,
                                           0.756479495,
                                           0.787147094,
                                           0.804572414,
                                           0.791806302,
                                           0.760234588,
                                           0.741109531,
                                           0.724118186,
                                           0.692829761,
                                           0.688465897,
                                           0.671806061,
                                           0.636461171,
                                           0.675314029,
                                           0.699134460,
                                           0.724305037,
                                           0.742556115,
                                           0.758504713};

    // Used to build the ID map 'theIDMap' (i.e., what is to be considered as a proton, etc...)

    const std::vector<int> protonsID = {2212, 3222, -101, -102, -103, -104};          
    const std::vector<int> antiprotonsID = {-2212, -3222};                            
    const std::vector<int> neutronsID = {2112, 3122, 3112, 3312, 3322, 3334, -3334};  
    const std::vector<int> antineutronsID = {-2112, -3122, -3112, -3312, -3322};      
    const std::vector<int> K0LsID = {130, 310};                                       
    const std::vector<int> KplussesID = {321};                                        
    const std::vector<int> KminussesID = {-321};                                      
    const std::vector<int> PiplussesID = {211};                                       
    const std::vector<int> PiminussesID = {-211};                                     
  };

  // TODO: Is this correct?
  // Thread safety
  static std::once_flag initializeOnce;
  CMS_THREAD_GUARD(initializeOnce) std::vector<std::vector<double> > NuclearInteraction::theRatiosMap;
  CMS_THREAD_GUARD(initializeOnce) std::map<int, int> NuclearInteraction::theIDMap;
}  // namespace fastsim

fastsim::NuclearInteraction::NuclearInteraction(const std::string& name, const edm::ParameterSet& cfg)
    : fastsim::InteractionModel(name), currentValuesWereSet(false) {
  // Full path to FullSim root file
  std::string fullPath;

  // Read from config file
  theDistCut = cfg.getParameter<double>("distCut");
  theHadronEnergy = cfg.getParameter<double>("hadronEnergy");
  inputFile = cfg.getUntrackedParameter<std::string>("inputFile", "");

  // The evolution of the interaction lengths with energy
  // initialize once for all possible instances
  std::call_once(initializeOnce, [this]() {
    theRatiosMap.resize(theHadronID.size());
    for (unsigned i = 0; i < theHadronID.size(); ++i) {
      for (unsigned j = 0; j < theHadronEN.size(); ++j) {
        theRatiosMap[i].push_back(theRatios[i * theHadronEN.size() + j]);
      }
    }

    // Build the ID map (i.e., what is to be considered as a proton, etc...)
    // Protons
    for (const auto& id : protonsID)
      theIDMap[id] = 2212;
    // Anti-Protons
    for (const auto& id : antiprotonsID)
      theIDMap[id] = -2212;
    // Neutrons
    for (const auto& id : neutronsID)
      theIDMap[id] = 2112;
    // Anti-Neutrons
    for (const auto& id : antineutronsID)
      theIDMap[id] = -2112;
    // K0L's
    for (const auto& id : K0LsID)
      theIDMap[id] = 130;
    // K+'s
    for (const auto& id : KplussesID)
      theIDMap[id] = 321;
    // K-'s
    for (const auto& id : KminussesID)
      theIDMap[id] = -321;
    // pi+'s
    for (const auto& id : PiplussesID)
      theIDMap[id] = 211;
    // pi-'s
    for (const auto& id : PiminussesID)
      theIDMap[id] = -211;
  });

  // Prepare the map of files
  // Loop over the particle names
  TFile* aVFile = nullptr;
  std::vector<TTree*> aVTree(theHadronEN.size());
  std::vector<TBranch*> aVBranch(theHadronEN.size());
  std::vector<NUEvent*> aVNUEvents(theHadronEN.size());
  std::vector<unsigned> aVCurrentEntry(theHadronEN.size());
  std::vector<unsigned> aVCurrentInteraction(theHadronEN.size());
  std::vector<unsigned> aVNumberOfEntries(theHadronEN.size());
  std::vector<unsigned> aVNumberOfInteractions(theHadronEN.size());
  std::vector<std::string> aVFileName(theHadronEN.size());
  std::vector<double> aVHadronCM(theHadronEN.size());
  theTrees.resize(theHadronNA.size());
  theBranches.resize(theHadronNA.size());
  theNUEvents.resize(theHadronNA.size());
  theCurrentEntry.resize(theHadronNA.size());
  theCurrentInteraction.resize(theHadronNA.size());
  theNumberOfEntries.resize(theHadronNA.size());
  theNumberOfInteractions.resize(theHadronNA.size());
  theFileNames.resize(theHadronNA.size());
  theHadronCM.resize(theHadronNA.size());
  theFile = aVFile;
  for (unsigned iname = 0; iname < theHadronNA.size(); ++iname) {
    theTrees[iname] = aVTree;
    theBranches[iname] = aVBranch;
    theNUEvents[iname] = aVNUEvents;
    theCurrentEntry[iname] = aVCurrentEntry;
    theCurrentInteraction[iname] = aVCurrentInteraction;
    theNumberOfEntries[iname] = aVNumberOfEntries;
    theNumberOfInteractions[iname] = aVNumberOfInteractions;
    theFileNames[iname] = aVFileName;
    theHadronCM[iname] = aVHadronCM;
  }

  // Read the information from a previous run (to keep reproducibility)
  currentValuesWereSet = this->read(inputFile);
  if (currentValuesWereSet)
    std::cout << "***WARNING*** You are reading nuclear-interaction information from the file " << inputFile
              << " created in an earlier run." << std::endl;

  // Open the file for saving the information of the current run
  myOutputFile.open("NuclearInteractionOutputFile.txt");
  myOutputBuffer = 0;

  // Open the root file
  edm::FileInPath myDataFile("FastSimulation/MaterialEffects/data/NuclearInteractions.root");
  fullPath = myDataFile.fullPath();
  theFile = TFile::Open(fullPath.c_str());

  // Open the trees
  for (unsigned iname = 0; iname < theHadronNA.size(); ++iname) {
    for (unsigned iene = 0; iene < theHadronEN.size(); ++iene) {
      std::ostringstream filename;
      double theEne = theHadronEN[iene];
      filename << "NuclearInteractionsVal_" << theHadronNA[iname] << "_E" << theEne << ".root";
      theFileNames[iname][iene] = filename.str();

      std::string treeName = "NuclearInteractions_" + theHadronNA[iname] + "_E" + std::to_string(int(theEne));
      theTrees[iname][iene] = (TTree*)theFile->Get(treeName.c_str());

      if (!theTrees[iname][iene])
        throw cms::Exception("FastSimulation/MaterialEffects") << "Tree with name " << treeName << " not found ";
      theBranches[iname][iene] = theTrees[iname][iene]->GetBranch("nuEvent");
      if (!theBranches[iname][iene])
        throw cms::Exception("FastSimulation/MaterialEffects")
            << "Branch with name nuEvent not found in " << theFileNames[iname][iene];

      theNUEvents[iname][iene] = new NUEvent();
      theBranches[iname][iene]->SetAddress(&theNUEvents[iname][iene]);
      theNumberOfEntries[iname][iene] = theTrees[iname][iene]->GetEntries();

      if (currentValuesWereSet) {
        theTrees[iname][iene]->GetEntry(theCurrentEntry[iname][iene]);
        unsigned NInteractions = theNUEvents[iname][iene]->nInteractions();
        theNumberOfInteractions[iname][iene] = NInteractions;
      }

      // Compute the corresponding cm energies of the nuclear interactions
      XYZTLorentzVector Proton(0., 0., 0., 0.986);
      XYZTLorentzVector Reference(
          0.,
          0.,
          std::sqrt(theHadronEN[iene] * theHadronEN[iene] - theHadronMA[iname] * theHadronMA[iname]),
          theHadronEN[iene]);
      theHadronCM[iname][iene] = (Reference + Proton).M();
    }
  }

  // Find the index for which EN = 4. (or thereabout)
  ien4 = 0;
  while (theHadronEN[ien4] < 4.0)
    ++ien4;

  gROOT->cd();
}

fastsim::NuclearInteraction::~NuclearInteraction() {
  // Close all local files
  // Among other things, this allows the TROOT destructor to end up
  // without crashing, while trying to close these files from outside
  theFile->Close();
  delete theFile;

  for (auto& vEvents : theNUEvents) {
    for (auto evtPtr : vEvents) {
      delete evtPtr;
    }
  }

  // Save data
  save();
  // Close the output file
  myOutputFile.close();
}

void fastsim::NuclearInteraction::interact(fastsim::Particle& particle,
                                           const SimplifiedGeometry& layer,
                                           std::vector<std::unique_ptr<fastsim::Particle> >& secondaries,
                                           const RandomEngineAndDistribution& random) {
  int pdgId = particle.pdgId();
  //
  // no valid PDGid
  //
  if (abs(pdgId) <= 100 || abs(pdgId) >= 1000000) {
    return;
  }

  double radLengths = layer.getThickness(particle.position(), particle.momentum());
  // TEC layers have some fudge factor for nuclear interactions
  radLengths *= layer.getNuclearInteractionThicknessFactor();
  //
  // no material
  //
  if (radLengths < 1E-10) {
    return;
  }

  // In case the events are not read from (old) saved file, then pick a random event from FullSim file
  if (!currentValuesWereSet) {
    currentValuesWereSet = true;
    for (unsigned iname = 0; iname < theHadronNA.size(); ++iname) {
      for (unsigned iene = 0; iene < theHadronEN.size(); ++iene) {
        theCurrentEntry[iname][iene] = (unsigned)(theNumberOfEntries[iname][iene] * random.flatShoot());

        theTrees[iname][iene]->GetEntry(theCurrentEntry[iname][iene]);
        unsigned NInteractions = theNUEvents[iname][iene]->nInteractions();
        theNumberOfInteractions[iname][iene] = NInteractions;

        theCurrentInteraction[iname][iene] = (unsigned)(theNumberOfInteractions[iname][iene] * random.flatShoot());
      }
    }
  }

  // Momentum of interacting hadron
  double pHadron = std::sqrt(particle.momentum().Vect().Mag2());

  //
  // The hadron has not enough momentum to create some relevant final state
  //
  if (pHadron < theHadronEnergy) {
    return;
  }

  // The particle type
  std::map<int, int>::const_iterator thePit = theIDMap.find(pdgId);
  // Use map for unique ID (e.g. proton = {2212, 3222, -101, -102, -103, -104})
  int thePid = (thePit != theIDMap.end() ? thePit->second : pdgId);

  // Is this particle type foreseen?
  unsigned fPid = abs(thePid);
  if (fPid != 211 && fPid != 130 && fPid != 321 && fPid != 2112 && fPid != 2212) {
    return;
  }

  // The hashed ID
  unsigned thePidIndex = index(thePid);
  // The inelastic interaction length at p(Hadron) = 5 GeV/c
  double theInelasticLength = radLengths * theLengthRatio[thePidIndex];

  // The elastic interaction length
  // The baroque parameterization is a fit to Fig. 40.13 of the PDG
  double ee = pHadron > 0.6 ? exp(-std::sqrt((pHadron - 0.6) / 1.122)) : exp(std::sqrt((0.6 - pHadron) / 1.122));
  double theElasticLength = (0.8753 * ee + 0.15) * theInelasticLength;

  // The total interaction length
  double theTotalInteractionLength = theInelasticLength + theElasticLength;

  //
  // Probability to interact is dl/L0 (maximum for 4 GeV Hadron)
  //
  double aNuclInteraction = -std::log(random.flatShoot());
  if (aNuclInteraction > theTotalInteractionLength) {
    return;
  }

  // The elastic part
  double elastic = random.flatShoot();
  if (elastic < theElasticLength / theTotalInteractionLength) {
    // Characteristic scattering angle for the elastic part
    // A of silicon
    double A = 28.0855;
    double theta0 = std::sqrt(3.) / std::pow(A, 1. / 3.) * particle.momentum().mass() / pHadron;

    // Draw an angle with theta/theta0*exp[(-theta/2theta0)**2] shape
    double theta = theta0 * std::sqrt(-2. * std::log(random.flatShoot()));
    double phi = 2. * M_PI * random.flatShoot();

    // Rotate the particle accordingly
    ROOT::Math::AxisAngle rotation1(orthogonal(particle.momentum().Vect()), theta);
    ROOT::Math::AxisAngle rotation2(particle.momentum().Vect(), phi);
    // Rotate!
    XYZVector rotated = rotation2((rotation1(particle.momentum().Vect())));

    // Create a daughter if the kink is large enough
    if (std::sin(theta) > theDistCut) {
      secondaries.emplace_back(
          new fastsim::Particle(pdgId,
                                particle.position(),
                                XYZTLorentzVector(rotated.X(), rotated.Y(), rotated.Z(), particle.momentum().E())));

      // Set the ClosestCharged Daughter thing for tracking
      if (particle.charge() != 0) {
        secondaries.back()->setMotherDeltaR(particle.momentum());
        secondaries.back()->setMotherPdgId(pdgId);
        secondaries.back()->setMotherSimTrackIndex(particle.simTrackIndex());
      }

      // The particle is destroyed
      particle.momentum().SetXYZT(0., 0., 0., 0.);
    } else {
      // If kink is small enough just rotate particle
      particle.momentum().SetXYZT(rotated.X(), rotated.Y(), rotated.Z(), particle.momentum().E());
    }
    // The inelastic part
  } else {
    // Avoid multiple map access
    const std::vector<double>& aHadronCM = theHadronCM[thePidIndex];
    const std::vector<double>& aRatios = theRatiosMap[thePidIndex];
    // Find the file with the closest c.m energy
    // The target nucleon
    XYZTLorentzVector Proton(0., 0., 0., 0.939);
    // The current particle
    const XYZTLorentzVector& Hadron = (const XYZTLorentzVector&)particle.momentum();
    // The smallest momentum for inelastic interactions
    double pMin = theHadronPMin[thePidIndex];
    // The corresponding smallest four vector
    XYZTLorentzVector Hadron0(0., 0., pMin, std::sqrt(pMin * pMin + particle.momentum().M2()));

    // The current centre-of-mass energy
    double ecm = (Proton + Hadron).M();

    // Get the files of interest (closest c.m. energies)
    unsigned ene1 = 0;
    unsigned ene2 = 0;
    // The smallest centre-of-mass energy
    double ecm1 = (Proton + Hadron0).M();

    double ecm2 = aHadronCM[0];
    double ratio1 = 0.;
    double ratio2 = aRatios[0];
    if (ecm > aHadronCM[0] && ecm < aHadronCM[aHadronCM.size() - 1]) {
      for (unsigned ene = 1; ene < aHadronCM.size() && ecm > aHadronCM[ene - 1]; ++ene) {
        if (ecm < aHadronCM[ene]) {
          ene2 = ene;
          ene1 = ene2 - 1;
          ecm1 = aHadronCM[ene1];
          ecm2 = aHadronCM[ene2];
          ratio1 = aRatios[ene1];
          ratio2 = aRatios[ene2];
        }
      }
    } else if (ecm > aHadronCM[aHadronCM.size() - 1]) {
      ene1 = aHadronCM.size() - 1;
      ene2 = aHadronCM.size() - 2;
      ecm1 = aHadronCM[ene1];
      ecm2 = aHadronCM[ene2];
      ratio1 = aRatios[ene2];
      ratio2 = aRatios[ene2];
    }

    // The inelastic part of the cross section depends cm energy
    double slope = (std::log10(ecm) - std::log10(ecm1)) / (std::log10(ecm2) - std::log10(ecm1));
    double inelastic = ratio1 + (ratio2 - ratio1) * slope;
    double inelastic4 = pHadron < 4. ? aRatios[ien4] : 1.;

    // Simulate an inelastic interaction
    if (elastic > 1. - (inelastic * theInelasticLength) / theTotalInteractionLength) {
      // Avoid mutliple map access
      std::vector<unsigned>& aCurrentInteraction = theCurrentInteraction[thePidIndex];
      std::vector<unsigned>& aNumberOfInteractions = theNumberOfInteractions[thePidIndex];
      std::vector<NUEvent*>& aNUEvents = theNUEvents[thePidIndex];

      // Choice of the file to read according the the log10(ecm) distance
      // and protection against low momentum proton and neutron that never interacts
      // (i.e., empty files)
      unsigned ene;
      if (random.flatShoot() < slope || aNumberOfInteractions[ene1] == 0) {
        ene = ene2;
      } else {
        ene = ene1;
      }

      // The boost characteristics
      XYZTLorentzVector theBoost = Proton + Hadron;
      theBoost /= theBoost.e();

      // Check we are not either at the end of an interaction bunch
      // or at the end of a file
      if (aCurrentInteraction[ene] == aNumberOfInteractions[ene]) {
        std::vector<unsigned>& aCurrentEntry = theCurrentEntry[thePidIndex];
        std::vector<unsigned>& aNumberOfEntries = theNumberOfEntries[thePidIndex];
        std::vector<TTree*>& aTrees = theTrees[thePidIndex];
        ++aCurrentEntry[ene];

        aCurrentInteraction[ene] = 0;
        if (aCurrentEntry[ene] == aNumberOfEntries[ene]) {
          aCurrentEntry[ene] = 0;
        }

        unsigned myEntry = aCurrentEntry[ene];
        aTrees[ene]->GetEntry(myEntry);
        aNumberOfInteractions[ene] = aNUEvents[ene]->nInteractions();
      }

      // Read the interaction
      NUEvent::NUInteraction anInteraction = aNUEvents[ene]->theNUInteractions()[aCurrentInteraction[ene]];

      unsigned firstTrack = anInteraction.first;
      unsigned lastTrack = anInteraction.last;

      // Some rotation around the boost axis, for more randomness
      XYZVector theAxis = theBoost.Vect().Unit();
      double theAngle = random.flatShoot() * 2. * M_PI;
      ROOT::Math::AxisAngle axisRotation(theAxis, theAngle);
      ROOT::Math::Boost axisBoost(theBoost.x(), theBoost.y(), theBoost.z());

      // A rotation to bring the particles back to the Hadron direction
      XYZVector zAxis(0., 0., 1.);
      XYZVector orthAxis = (zAxis.Cross(theBoost.Vect())).Unit();
      double orthAngle = acos(theBoost.Vect().Unit().Z());
      ROOT::Math::AxisAngle orthRotation(orthAxis, orthAngle);

      // Loop on the nuclear interaction products
      for (unsigned iTrack = firstTrack; iTrack <= lastTrack; ++iTrack) {
        NUEvent::NUParticle aParticle = aNUEvents[ene]->theNUParticles()[iTrack];

        // Add a RawParticle with the proper energy in the c.m. frame of
        // the nuclear interaction
        double energy = std::sqrt(aParticle.px * aParticle.px + aParticle.py * aParticle.py +
                                  aParticle.pz * aParticle.pz + aParticle.mass * aParticle.mass / (ecm * ecm));

        XYZTLorentzVector daugtherMomentum(aParticle.px * ecm, aParticle.py * ecm, aParticle.pz * ecm, energy * ecm);

        // Rotate to the collision axis
        XYZVector rotated = orthRotation(daugtherMomentum.Vect());
        // Rotate around the boost axis for more randomness
        rotated = axisRotation(rotated);

        // Rotated the daughter
        daugtherMomentum.SetXYZT(rotated.X(), rotated.Y(), rotated.Z(), daugtherMomentum.E());

        // Boost it in the lab frame
        daugtherMomentum = axisBoost(daugtherMomentum);

        // Create secondary
        secondaries.emplace_back(new fastsim::Particle(aParticle.id, particle.position(), daugtherMomentum));

        // The closestCharged Daughter thing for tracking
        // BUT: 'aParticle' also has to be charged, only then the mother should be set
        // Unfortunately, NUEvent::NUParticle does not contain any info about the charge
        // Did some tests and effect is absolutely negligible!
        if (particle.charge() != 0) {
          secondaries.back()->setMotherDeltaR(particle.momentum());
          secondaries.back()->setMotherPdgId(pdgId);
          secondaries.back()->setMotherSimTrackIndex(particle.simTrackIndex());
        }
      }

      // The particle is destroyed
      particle.momentum().SetXYZT(0., 0., 0., 0.);

      // This is a note from previous version of code but I don't understand it:
      // ERROR The way this loops through the events breaks
      // replay. Which events are retrieved depends on
      // which previous events were processed.

      // Increment for next time
      ++aCurrentInteraction[ene];

      // Simulate a stopping hadron (low momentum)
    } else if (pHadron < 4. && elastic > 1. - (inelastic4 * theInelasticLength) / theTotalInteractionLength) {
      // The particle is destroyed
      particle.momentum().SetXYZT(0., 0., 0., 0.);
    }
  }
}

void fastsim::NuclearInteraction::save() {
  // Size of buffer
  ++myOutputBuffer;

  // Periodically close the current file and open a new one
  if (myOutputBuffer / 1000 * 1000 == myOutputBuffer) {
    myOutputFile.close();
    myOutputFile.open("NuclearInteractionOutputFile.txt");
  }
  //
  unsigned size1 = theCurrentEntry.size() * theCurrentEntry.begin()->size();
  std::vector<unsigned> theCurrentEntries;
  theCurrentEntries.resize(size1);
  size1 *= sizeof(unsigned);
  //
  unsigned size2 = theCurrentInteraction.size() * theCurrentInteraction.begin()->size();
  std::vector<unsigned> theCurrentInteractions;
  theCurrentInteractions.resize(size2);
  size2 *= sizeof(unsigned);

  // Save the current entries
  std::vector<std::vector<unsigned> >::const_iterator aCurrentEntry = theCurrentEntry.begin();
  std::vector<std::vector<unsigned> >::const_iterator lastCurrentEntry = theCurrentEntry.end();
  unsigned allEntries = 0;
  for (; aCurrentEntry != lastCurrentEntry; ++aCurrentEntry) {
    unsigned size = aCurrentEntry->size();
    for (unsigned iene = 0; iene < size; ++iene)
      theCurrentEntries[allEntries++] = (*aCurrentEntry)[iene];
  }

  // Save the current interactions
  std::vector<std::vector<unsigned> >::const_iterator aCurrentInteraction = theCurrentInteraction.begin();
  std::vector<std::vector<unsigned> >::const_iterator lastCurrentInteraction = theCurrentInteraction.end();
  unsigned allInteractions = 0;
  for (; aCurrentInteraction != lastCurrentInteraction; ++aCurrentInteraction) {
    unsigned size = aCurrentInteraction->size();
    for (unsigned iene = 0; iene < size; ++iene)
      theCurrentInteractions[allInteractions++] = (*aCurrentInteraction)[iene];
  }
  // Write to file
  myOutputFile.write((const char*)(&theCurrentEntries.front()), size1);
  myOutputFile.write((const char*)(&theCurrentInteractions.front()), size2);
  myOutputFile.flush();
}

bool fastsim::NuclearInteraction::read(std::string inputFile) {
  std::ifstream myInputFile;
  struct stat results;
  //
  unsigned size1 = theCurrentEntry.size() * theCurrentEntry.begin()->size();
  std::vector<unsigned> theCurrentEntries;
  theCurrentEntries.resize(size1);
  size1 *= sizeof(unsigned);
  //
  unsigned size2 = theCurrentInteraction.size() * theCurrentInteraction.begin()->size();
  std::vector<unsigned> theCurrentInteractions;
  theCurrentInteractions.resize(size2);
  size2 *= sizeof(unsigned);
  //
  unsigned size = 0;

  // Open the file (if any), otherwise return false
  myInputFile.open(inputFile.c_str());
  if (myInputFile.is_open()) {
    // Get the size of the file
    if (stat(inputFile.c_str(), &results) == 0)
      size = results.st_size;
    else
      return false;  // Something is wrong with that file !

    // Position the pointer just before the last record
    myInputFile.seekg(size - size1 - size2);
    myInputFile.read((char*)(&theCurrentEntries.front()), size1);
    myInputFile.read((char*)(&theCurrentInteractions.front()), size2);
    myInputFile.close();

    // Read the current entries
    std::vector<std::vector<unsigned> >::iterator aCurrentEntry = theCurrentEntry.begin();
    std::vector<std::vector<unsigned> >::iterator lastCurrentEntry = theCurrentEntry.end();
    unsigned allEntries = 0;
    for (; aCurrentEntry != lastCurrentEntry; ++aCurrentEntry) {
      unsigned size = aCurrentEntry->size();
      for (unsigned iene = 0; iene < size; ++iene)
        (*aCurrentEntry)[iene] = theCurrentEntries[allEntries++];
    }

    // Read the current interactions
    std::vector<std::vector<unsigned> >::iterator aCurrentInteraction = theCurrentInteraction.begin();
    std::vector<std::vector<unsigned> >::iterator lastCurrentInteraction = theCurrentInteraction.end();
    unsigned allInteractions = 0;
    for (; aCurrentInteraction != lastCurrentInteraction; ++aCurrentInteraction) {
      unsigned size = aCurrentInteraction->size();
      for (unsigned iene = 0; iene < size; ++iene)
        (*aCurrentInteraction)[iene] = theCurrentInteractions[allInteractions++];
    }

    return true;
  }

  return false;
}

unsigned fastsim::NuclearInteraction::index(int thePid) {
  // Find hashed particle ID
  unsigned myIndex = 0;
  while (thePid != theHadronID[myIndex])
    ++myIndex;
  return myIndex;
}

XYZVector fastsim::NuclearInteraction::orthogonal(const XYZVector& aVector) const {
  double x = fabs(aVector.X());
  double y = fabs(aVector.Y());
  double z = fabs(aVector.Z());

  if (x < y)
    return x < z ? XYZVector(0., aVector.Z(), -aVector.Y()) : XYZVector(aVector.Y(), -aVector.X(), 0.);
  else
    return y < z ? XYZVector(-aVector.Z(), 0., aVector.X()) : XYZVector(aVector.Y(), -aVector.X(), 0.);
}

DEFINE_EDM_PLUGIN(fastsim::InteractionModelFactory, fastsim::NuclearInteraction, "fastsim::NuclearInteraction");