NuclearInteractionSimulator.cc

Go to the documentation of this file.

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

#include "FastSimulation/MaterialEffects/interface/NuclearInteractionSimulator.h"
#include "FastSimulation/Particle/interface/makeParticle.h"
#include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"

#include "FastSimDataFormats/NuclearInteractions/interface/NUEvent.h"

//#include "DQMServices/Core/interface/DaqMonitorBEInterface.h"
//#include "FWCore/ServiceRegistry/interface/Service.h"

#include <iostream>
#include <sys/stat.h>
#include <cmath>
#include "TFile.h"
#include "TTree.h"
#include "TROOT.h"

// Internal variable and vectors with name started frm "thePion" means
// vectors/variable not only for pions but for all type of hadrons
// treated inside this code

NuclearInteractionSimulator::NuclearInteractionSimulator(std::vector<double>& hadronEnergies,
                                                         std::vector<int>& hadronTypes,
                                                         std::vector<std::string>& hadronNames,
                                                         std::vector<double>& hadronMasses,
                                                         std::vector<double>& hadronPMin,
                                                         double pionEnergy,
                                                         std::vector<double>& lengthRatio,
                                                         std::vector<std::vector<double> >& ratios,
                                                         std::map<int, int>& idMap,
                                                         std::string inputFile,
                                                         unsigned int distAlgo,
                                                         double distCut)
    : MaterialEffectsSimulator(),
      thePionEN(hadronEnergies),
      thePionID(hadronTypes),
      thePionNA(hadronNames),
      thePionMA(hadronMasses),
      thePionPMin(hadronPMin),
      thePionEnergy(pionEnergy),
      theLengthRatio(lengthRatio),
      theRatios(ratios),
      theIDMap(idMap),
      theDistAlgo(distAlgo),
      theDistCut(distCut),
      currentValuesWereSet(false) {
  std::string fullPath;

  // Prepare the map of files
  // Loop over the particle names
  TFile* aVFile = nullptr;
  std::vector<TTree*> aVTree(thePionEN.size(), static_cast<TTree*>(nullptr));
  std::vector<TBranch*> aVBranch(thePionEN.size(), static_cast<TBranch*>(nullptr));
  std::vector<NUEvent*> aVNUEvents(thePionEN.size(), static_cast<NUEvent*>(nullptr));
  std::vector<unsigned> aVCurrentEntry(thePionEN.size(), static_cast<unsigned>(0));
  std::vector<unsigned> aVCurrentInteraction(thePionEN.size(), static_cast<unsigned>(0));
  std::vector<unsigned> aVNumberOfEntries(thePionEN.size(), static_cast<unsigned>(0));
  std::vector<unsigned> aVNumberOfInteractions(thePionEN.size(), static_cast<unsigned>(0));
  std::vector<std::string> aVFileName(thePionEN.size(), static_cast<std::string>(""));
  std::vector<double> aVPionCM(thePionEN.size(), static_cast<double>(0));

  theTrees.resize(thePionNA.size());
  theBranches.resize(thePionNA.size());
  theNUEvents.resize(thePionNA.size());
  theCurrentEntry.resize(thePionNA.size());
  theCurrentInteraction.resize(thePionNA.size());
  theNumberOfEntries.resize(thePionNA.size());
  theNumberOfInteractions.resize(thePionNA.size());
  theFileNames.resize(thePionNA.size());
  thePionCM.resize(thePionNA.size());
  theFile = aVFile;
  for (unsigned iname = 0; iname < thePionNA.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;
    thePionCM[iname] = aVPionCM;
  }

  // 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 files
  //  for ( unsigned file=0; file<theFileNames.size(); ++file ) {
  edm::FileInPath myDataFile("FastSimulation/MaterialEffects/data/NuclearInteractions.root");

  fullPath = myDataFile.fullPath();
  theFile = TFile::Open(fullPath.c_str());

  unsigned fileNb = 0;
  for (unsigned iname = 0; iname < thePionNA.size(); ++iname) {
    for (unsigned iene = 0; iene < thePionEN.size(); ++iene) {
      //std::cout << "iname/iene " << iname << " " << iene << std::endl;
      std::ostringstream filename;
      double theEne = thePionEN[iene];
      filename << "NuclearInteractionsVal_" << thePionNA[iname] << "_E" << theEne << ".root";
      theFileNames[iname][iene] = filename.str();
      //std::cout << "thePid/theEne " << thePionID[iname] << " " << theEne << std::endl;

      ++fileNb;
      std::string treeName = "NuclearInteractions_" + thePionNA[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");
      //std::cout << "The branch = " << theBranches[iname][iene] << std::endl;
      if (!theBranches[iname][iene])
        throw cms::Exception("FastSimulation/MaterialEffects")
            << "Branch with name nuEvent not found in " << theFileNames[iname][iene];
      //
      theNUEvents[iname][iene] = new NUEvent();
      //std::cout << "The branch = " << theBranches[iname][iene] << std::endl;
      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(thePionEN[iene] * thePionEN[iene] - thePionMA[iname] * thePionMA[iname]), thePionEN[iene]);
      thePionCM[iname][iene] = (Reference + Proton).M();
    }
  }

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

  gROOT->cd();

  // Information (Should be on LogInfo)
  //  std::cout << " ---> A total of " << fileNb
  //        << " nuclear-interaction files was sucessfully open" << std::endl;

  //  dbe = edm::Service<DaqMonitorBEInterface>().operator->();
  //  htot = dbe->book1D("Total", "All particles",150,0.,150.);
  //  helas = dbe->book1D("Elastic", "Elastic interactions",150,0.,150.);
  //  hinel = dbe->book1D("Inelastic", "Inelastic interactions",150,0.,150.);
  //  hscatter = dbe->book1D("Scattering","Elastic Scattering angle",200,0.,2.);
  //  hscatter2 = dbe->book2D("Scattering2","Elastic Scattering angle vs p",100,0.,10.,200,0.,2.);
  //  hAfter = dbe->book1D("eAfter","Energy after collision",200,0.,4.);
  //  hAfter2 = dbe->book2D("eAfter2","Energy after collision",100,-2.5,2.5,100,0.,4);
  //  hAfter3 = dbe->book2D("eAfter3","Energy after collision",100,0.,1000.,100,0.,4);
}

NuclearInteractionSimulator::~NuclearInteractionSimulator() {
  // 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;
    }
  }

  // Close the output file
  myOutputFile.close();

  //  dbe->save("test.root");
}

void NuclearInteractionSimulator::compute(ParticlePropagator& Particle, RandomEngineAndDistribution const* random) {
  if (!currentValuesWereSet) {
    currentValuesWereSet = true;
    for (unsigned iname = 0; iname < thePionNA.size(); ++iname) {
      for (unsigned iene = 0; iene < thePionEN.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());
      }
    }
  }

  // Read a Nuclear Interaction in a random manner

  double pHadron = std::sqrt(Particle.particle().Vect().Mag2());
  //  htot->Fill(pHadron);

  // The hadron has enough momentum to create some relevant final state
  if (pHadron > thePionEnergy) {
    // The particle type
    std::map<int, int>::const_iterator thePit = theIDMap.find(Particle.particle().pid());

    int thePid = thePit != theIDMap.end() ? thePit->second : Particle.particle().pid();

    // Is this particle type foreseen?
    unsigned fPid = abs(thePid);
    if (fPid != 211 && fPid != 130 && fPid != 321 && fPid != 2112 && fPid != 2212) {
      return;
      //std::cout << "Unknown particle type = " << thePid << std::endl;
      //thePid = 211;
    }

    // The inelastic interaction length at p(pion) = 5 GeV/c
    unsigned thePidIndex = index(thePid);
    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)
                              //    double theElasticLength = ( 0.15 + 0.195 / log(pHadron/0.4) )
                              //    double theElasticLength = ( 0.15 + 0.305 / log(pHadron/0.35) )
                              * theInelasticLength;

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

    // Probability to interact is dl/L0 (maximum for 4 GeV pion)
    double aNuclInteraction = -std::log(random->flatShoot());
    if (aNuclInteraction < theTotalInteractionLength) {
      // The elastic part
      double elastic = random->flatShoot();
      if (elastic < theElasticLength / theTotalInteractionLength) {
        //          helas->Fill(pHadron);

        // Characteristic scattering angle for the elastic part
        double theta0 = std::sqrt(3.) / std::pow(theA(), 1. / 3.) * Particle.particle().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. * 3.14159265358979323 * random->flatShoot();

        // Rotate the particle accordingly
        RawParticle::Rotation rotation1(orthogonal(Particle.particle().Vect()), theta);
        RawParticle::Rotation rotation2(Particle.particle().Vect(), phi);
        Particle.particle().rotate(rotation1);
        Particle.particle().rotate(rotation2);

        // Distance
        double distance = std::sin(theta);

        // Create a daughter if the kink is large engough
        if (distance > theDistCut) {
          _theUpdatedState.reserve(1);
          _theUpdatedState.clear();
          _theUpdatedState.emplace_back(Particle.particle());
        }

        //  hscatter->Fill(myTheta);
        //  hscatter2->Fill(pHadron,myTheta);

      }

      // The inelastic part
      else {
        // Avoid multiple map access
        const std::vector<double>& aPionCM = thePionCM[thePidIndex];
        const std::vector<double>& aRatios = theRatios[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;
        // The smallest momentum for inelastic interactions
        double pMin = thePionPMin[thePidIndex];
        // The correspong smallest four vector
        XYZTLorentzVector Hadron0(0., 0., pMin, std::sqrt(pMin * pMin + Particle.particle().M2()));

        // The current centre-of-mass energy
        double ecm = (Proton + Hadron).M();
        //std::cout << "Proton = " << Proton << std::endl;
        //std::cout << "Hadron = " << Hadron << std::endl;
        //std::cout << "ecm = " << ecm << std::endl;
        // Get the files of interest (closest c.m. energies)
        unsigned ene1 = 0;
        unsigned ene2 = 0;
        // The smallest centre-of-mass energy
        //  double ecm1=1.63;
        double ecm1 = (Proton + Hadron0).M();
        //std::cout << "ecm1 = " << ecm1 << std::endl;
        //std::cout << "ecm[0] = " << aPionCM[0] << std::endl;
        //std::cout << "ecm[11] = " << aPionCM [ aPionCM.size()-1 ] << std::endl;
        double ecm2 = aPionCM[0];
        double ratio1 = 0.;
        double ratio2 = aRatios[0];
        if (ecm > aPionCM[0] && ecm < aPionCM[aPionCM.size() - 1]) {
          for (unsigned ene = 1; ene < aPionCM.size() && ecm > aPionCM[ene - 1]; ++ene) {
            if (ecm < aPionCM[ene]) {
              ene2 = ene;
              ene1 = ene2 - 1;
              ecm1 = aPionCM[ene1];
              ecm2 = aPionCM[ene2];
              ratio1 = aRatios[ene1];
              ratio2 = aRatios[ene2];
            }
          }
        } else if (ecm > aPionCM[aPionCM.size() - 1]) {
          ene1 = aPionCM.size() - 1;
          ene2 = aPionCM.size() - 2;
          ecm1 = aPionCM[ene1];
          ecm2 = aPionCM[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.;

        //std::cout << "Inelastic = " << ratio1 << " " << ratio2 << " " << inelastic << std::endl;
        //      std::cout << "Energy/Inelastic : "
        //      << Hadron.e() << " " << inelastic << std::endl;

        // 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];
          //      hinel->Fill(pHadron);
          //      std::cout << "INELASTIC INTERACTION ! "
          //            << pHadron << " " << theInelasticLength << " "
          //            << inelastic * theInelasticLength << std::endl;
          // 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;

          //std::cout << "Ecm1/2 = " << ecm1 << " " << ecm2 << std::endl;
          //std::cout << "Ratio1/2 = " << ratio1 << " " << ratio2 << std::endl;
          //std::cout << "Ene = " << ene << " slope = " << slope << std::endl;

          //std::cout << "Pion energy = " << Hadron.E()
          //        << "File chosen " << theFileNames[thePidIndex][ene]
          //        << std::endl;

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

          // std::cout << "File chosen : " << thePid << "/" << ene
          //        << " Current interaction = " << aCurrentInteraction[ene]
          //        << " Total interactions = " << aNumberOfInteractions[ene]
          //        << std::endl;
          //      theFiles[thePidIndex][ene]->cd();
          //      gDirectory->ls();

          // 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::cout << "End of interaction bunch ! ";
            std::vector<unsigned>& aCurrentEntry = theCurrentEntry[thePidIndex];
            std::vector<unsigned>& aNumberOfEntries = theNumberOfEntries[thePidIndex];
            std::vector<TTree*>& aTrees = theTrees[thePidIndex];
            ++aCurrentEntry[ene];
            // std::cerr << "Read the next entry "
            //           << aCurrentEntry[ene] << std::endl;
            aCurrentInteraction[ene] = 0;
            if (aCurrentEntry[ene] == aNumberOfEntries[ene]) {
              aCurrentEntry[ene] = 0;
              //  std::cout << "End of file - Rewind! " << std::endl;
            }
            unsigned myEntry = aCurrentEntry[ene];
            // std::cout << "The new entry " << myEntry
            //           << " is read ... in TTree " << aTrees[ene] << " ";
            aTrees[ene]->GetEntry(myEntry);
            // std::cout << "The number of interactions in the new entry is ... ";
            aNumberOfInteractions[ene] = aNUEvents[ene]->nInteractions();
            // std::cout << aNumberOfInteractions[ene] << std::endl;
          }

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

          unsigned firstTrack = anInteraction.first;
          unsigned lastTrack = anInteraction.last;
          //      std::cout << "First and last tracks are " << firstTrack << " " << lastTrack << std::endl;

          _theUpdatedState.reserve(lastTrack - firstTrack + 1);
          _theUpdatedState.clear();

          double distMin = 1E99;

          // Some rotation around the boost axis, for more randomness
          XYZVector theAxis = theBoost.Vect().Unit();
          double theAngle = random->flatShoot() * 2. * 3.14159265358979323;
          RawParticle::Rotation axisRotation(theAxis, theAngle);
          RawParticle::Boost axisBoost(theBoost.x(), theBoost.y(), theBoost.z());

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

          // A few checks
          // double eAfter = 0.;

          // Loop on the nuclear interaction products
          for (unsigned iTrack = firstTrack; iTrack <= lastTrack; ++iTrack) {
            unsigned idaugh = iTrack - firstTrack;
            NUEvent::NUParticle aParticle = aNUEvents[ene]->theNUParticles()[iTrack];
            //  std::cout << "Track " << iTrack
            //        << " id/px/py/pz/mass "
            //        << aParticle.id << " "
            //        << aParticle.px << " "
            //        << aParticle.py << " "
            //        << aParticle.pz << " "
            //        << aParticle.mass << " " << endl;

            // 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));

            RawParticle& aDaughter = _theUpdatedState.emplace_back(makeParticle(
                Particle.particleDataTable(),
                aParticle.id,
                XYZTLorentzVector(aParticle.px * ecm, aParticle.py * ecm, aParticle.pz * ecm, energy * ecm)));
            // Rotate to the collision axis
            aDaughter.rotate(orthRotation);

            // Rotate around the boost axis for more randomness
            aDaughter.rotate(axisRotation);

            // Boost it in the lab frame
            aDaughter.boost(axisBoost);

            // Store the closest daughter index (for later tracking purposes, so charged particles only)
            double distance = distanceToPrimary(Particle.particle(), aDaughter);
            // Find the closest daughter, if closer than a given upper limit.
            if (distance < distMin && distance < theDistCut) {
              distMin = distance;
              theClosestChargedDaughterId = idaugh;
            }

            // eAfter += aDaughter.E();
          }

          /*
      double eBefore = Particle.E();
      double rapid = Particle.momentum().Eta();
      if ( eBefore > 0. ) { 
        hAfter->Fill(eAfter/eBefore);
        hAfter2->Fill(rapid,eAfter/eBefore);
        hAfter3->Fill(eBefore,eAfter/eBefore);
      }
      */

          // 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) {
          // A fake particle with 0 momentum as a daughter!
          _theUpdatedState.reserve(1);
          _theUpdatedState.clear();
          _theUpdatedState.emplace_back(
              makeParticle(Particle.particleDataTable(), 22, XYZTLorentzVector(0., 0., 0., 0.)));
        }
      }
    }
  }
}

double NuclearInteractionSimulator::distanceToPrimary(const RawParticle& Particle, const RawParticle& aDaughter) const {
  double distance = 2E99;

  // Compute the distance only for charged primaries
  if (fabs(Particle.charge()) > 1E-12) {
    // The secondary must have the same charge
    double chargeDiff = fabs(aDaughter.charge() - Particle.charge());
    if (fabs(chargeDiff) < 1E-12) {
      // Here are two distance definitions * to be tuned *
      switch (theDistAlgo) {
        case 1:
          // sin(theta12)
          distance = (aDaughter.Vect().Unit().Cross(Particle.Vect().Unit())).R();
          break;

        case 2:
          // sin(theta12) * p1/p2
          distance = (aDaughter.Vect().Cross(Particle.Vect())).R() / aDaughter.Vect().Mag2();
          break;

        default:
          // Should not happen
          distance = 2E99;
          break;
      }
    }
  }

  return distance;
}

void NuclearInteractionSimulator::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");
    //    myOutputFile.seekp(0); // No need to rewind in that case
  }
  //
  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];
  }
  //
  myOutputFile.write((const char*)(&theCurrentEntries.front()), size1);
  myOutputFile.write((const char*)(&theCurrentInteractions.front()), size2);
  myOutputFile.flush();
}

bool NuclearInteractionSimulator::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)
  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 NuclearInteractionSimulator::index(int thePid) {
  unsigned myIndex = 0;
  while (thePid != thePionID[myIndex])
    ++myIndex;
  //  std::cout << "pid/index = " << thePid << " " << myIndex << std::endl;
  return myIndex;
}