fastsim::NuclearInteraction Class Reference

Implementation of nuclear interactions of hadrons in the tracker layers (based on fully simulated interactions). More...

Inheritance diagram for fastsim::NuclearInteraction:

Public Member Functions
void	interact (fastsim::Particle &particle, const SimplifiedGeometry &layer, std::vector< std::unique_ptr< fastsim::Particle > > &secondaries, const RandomEngineAndDistribution &random) override
	Perform the interaction. More...
	NuclearInteraction (const std::string &name, const edm::ParameterSet &cfg)
	Constructor. More...
	~NuclearInteraction () override
	Default destructor. More...
Public Member Functions inherited from fastsim::InteractionModel
const std::string	getName ()
	Return (unique) name of this interaction. More...
	InteractionModel (std::string name)
	Constructor. More...
virtual void	registerProducts (edm::ProducesCollector) const
	In case interaction produces and stores content in the event (e.g. TrackerSimHits). More...
virtual void	storeProducts (edm::Event &iEvent)
	In case interaction produces and stores content in the event (e.g. TrackerSimHits). More...
virtual	~InteractionModel ()
	Default destructor. More...

Private Member Functions
unsigned	index (int thePid)
	Return a hashed index for a given particle ID. More...
XYZVector	orthogonal (const XYZVector &aVector) const
	Return an orthogonal vector. More...
bool	read (std::string inputFile)
	Read the nuclear interactions from a file, so you can reproduce the events (e.g. in case of a crash). More...
void	save ()
	Save the nuclear interactions to a file, so you can reproduce the events (e.g. in case of a crash). More...

Private Attributes
const std::vector< int >	antineutronsID = {-2112, -3122, -3112, -3312, -3322}
	PdgID of anti-neutrons. More...
const std::vector< int >	antiprotonsID = {-2212, -3222}
	PdgID of anti-protons. More...
bool	currentValuesWereSet
	Read data from file that was created in a previous run. More...
unsigned	ien4
	Find the index for which EN = 4. More...
std::string	inputFile
	Directory/Name of input file in case you want to read interactions from file. More...
const std::vector< int >	K0LsID = {130, 310}
	PdgID of K0. More...
const std::vector< int >	KminussesID = {-321}
	PdgID of K-. More...
const std::vector< int >	KplussesID = {321}
	PdgID of K+. More...
unsigned	myOutputBuffer
	Needed to save interactions to file. More...
std::ofstream	myOutputFile
	Output file to save interactions. More...
const std::vector< int >	neutronsID = {2112, 3122, 3112, 3312, 3322, 3334, -3334}
	PdgID of neutrons. More...
const std::vector< int >	PiminussesID = {-211}
	PdgID of pi-. More...
const std::vector< int >	PiplussesID = {211}
	PdgID of pt+. More...
const std::vector< int >	protonsID = {2212, 3222, -101, -102, -103, -104}
	PdgID of protons. More...
std::vector< std::vector < TBranch * > >	theBranches
	Necessary to read the FullSim interactions. More...
std::vector< std::vector < unsigned > >	theCurrentEntry
	Necessary to read the FullSim interactions. More...
std::vector< std::vector < unsigned > >	theCurrentInteraction
	Necessary to read the FullSim interactions. More...
double	theDistCut
	Cut on deltaR for the FastSim Tracking (ClosestChargedDaughter algorithm) More...
TFile *	theFile = nullptr
	Necessary to read the FullSim interactions. More...
std::vector< std::vector < std::string > >	theFileNames
	Necessary to read the FullSim interactions. More...
std::vector< std::vector < double > >	theHadronCM
	Necessary to read the FullSim interactions. More...
const std::vector< double >	theHadronEN
	Filled into 'theRatiosMap' (evolution of the interaction lengths with energy) More...
double	theHadronEnergy
	Minimum energy for nuclear interaction. More...
const std::vector< int >	theHadronID = {211, -211, 130, 321, -321, 2212, -2212, 2112, -2112}
	ID of the hadrons. More...
const std::vector< double >	theHadronMA
	Masses of the hadrons. More...
const std::vector< std::string >	theHadronNA
	Names of the hadrons. More...
const std::vector< double >	theHadronPMin = {0.7, 0.0, 1.0, 1.0, 0.0, 1.1, 0.0, 1.1, 0.0}
	Smallest momentum for inelastic interactions. More...
const std::vector< double >	theLengthRatio
	Inelastic interaction length at p(Hadron) = 5 GeV/c (relativ to radionLength of material) More...
std::vector< std::vector < NUEvent * > >	theNUEvents
	Necessary to read the FullSim interactions. More...
std::vector< std::vector < unsigned > >	theNumberOfEntries
	Necessary to read the FullSim interactions. More...
std::vector< std::vector < unsigned > >	theNumberOfInteractions
	Necessary to read the FullSim interactions. More...
const std::vector< double >	theRatios
	Filled into 'theRatiosMap' (evolution of the interaction lengths with energy) More...
std::vector< std::vector < TTree * > >	theTrees
	Necessary to read the FullSim interactions. More...

Static Private Attributes
static std::map< int, int >	theIDMap
	Build the ID map (i.e., what is to be considered as a proton, etc...) More...
static std::vector < std::vector< double > >	theRatiosMap
	The evolution of the interaction lengths with energy. More...

Detailed Description

Implementation of nuclear interactions of hadrons in the tracker layers (based on fully simulated interactions).

Computes the probability for hadrons to interact with a nucleon of the tracker material (inelastically) and then reads a nuclear interaction randomly from multiple fully simulated files. Also, another implementation of nuclear interactions can be used that is based on G4 (NuclearInteractionFTF).

Definition at line 58 of file NuclearInteraction.cc.

Constructor & Destructor Documentation

fastsim::NuclearInteraction::NuclearInteraction	(	const std::string &	name,
		const edm::ParameterSet &	cfg
	)

Constructor.

Definition at line 348 of file NuclearInteraction.cc.

References looper::cfg, gather_cfg::cout, Exception, lut2db_cfg::filename, contentValuesFiles::fullPath, edm::FileInPath::fullPath(), mps_fire::i, gpuClustering::id, fastsim::initializeOnce, writeEcalDQMStatus::inputFile, dqmiolumiharvest::j, SiPixelLorentzAngle_cfi::read, mathSSE::sqrt(), AlCaHLTBitMon_QueryRunRegistry::string, and interactiveExample::theFile.

    : 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
  unsigned fileNb = 0;
  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();

      ++fileNb;
      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 ( )

override

Default destructor.

Definition at line 496 of file NuclearInteraction.cc.

References cuy::save, and interactiveExample::theFile.

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

Member Function Documentation

unsigned fastsim::NuclearInteraction::index ( int  thePid )

private

Return a hashed index for a given particle ID.

Definition at line 899 of file NuclearInteraction.cc.

Referenced by BeautifulSoup.PageElement::_invert().

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

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

overridevirtual

Perform the interaction.

Parameters

particle	The particle that interacts with the matter.
layer	The detector layer that interacts with the particle.
secondaries	Particles that are produced in the interaction (if any).
random	The Random Engine.

Implements fastsim::InteractionModel.

Definition at line 515 of file NuclearInteraction.cc.

References funct::A, funct::abs(), fastsim::Particle::charge(), relval_parameters_module::energy, funct::exp(), NUEvent::NUInteraction::first, RandomEngineAndDistribution::flatShoot(), fastsim::SimplifiedGeometry::getNuclearInteractionThicknessFactor(), fastsim::SimplifiedGeometry::getThickness(), NUEvent::NUParticle::id, NUEvent::NUInteraction::last, log, M_PI, NUEvent::NUParticle::mass, fastsim::Particle::momentum(), fastsim::Particle::pdgId(), HLT_FULL_cff::pMin, fastsim::Particle::position(), funct::pow(), NUEvent::NUParticle::px, NUEvent::NUParticle::py, NUEvent::NUParticle::pz, fastsim::Particle::simTrackIndex(), funct::sin(), slope, mathSSE::sqrt(), theta(), and MetAnalyzer::zAxis.

                                                                                      {
  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.);
    }
  }
}

XYZVector fastsim::NuclearInteraction::orthogonal ( const XYZVector &  aVector ) const

private

Return an orthogonal vector.

Definition at line 907 of file NuclearInteraction.cc.

References x, y, and z.

                                                                              {
  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.);
}

bool fastsim::NuclearInteraction::read ( std::string  inputFile )

private

Read the nuclear interactions from a file, so you can reproduce the events (e.g. in case of a crash).

Definition at line 842 of file NuclearInteraction.cc.

References bookConverter::results, findQualityFiles::size, and edm_modernize_messagelogger::stat.

Referenced by edmIntegrityCheck.PublishToFileSystem::get().

                                                        {
  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;
}

void fastsim::NuclearInteraction::save ( )

private

Save the nuclear interactions to a file, so you can reproduce the events (e.g. in case of a crash).

Definition at line 797 of file NuclearInteraction.cc.

References findQualityFiles::size.

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

Member Data Documentation

const std::vector<int> fastsim::NuclearInteraction::antineutronsID = {-2112, -3122, -3112, -3312, -3322}

private

PdgID of anti-neutrons.

Definition at line 333 of file NuclearInteraction.cc.

const std::vector<int> fastsim::NuclearInteraction::antiprotonsID = {-2212, -3222}

private

PdgID of anti-protons.

Definition at line 331 of file NuclearInteraction.cc.

bool fastsim::NuclearInteraction::currentValuesWereSet

private

Read data from file that was created in a previous run.

Definition at line 115 of file NuclearInteraction.cc.

unsigned fastsim::NuclearInteraction::ien4

private

Find the index for which EN = 4.

Definition at line 110 of file NuclearInteraction.cc.

std::string fastsim::NuclearInteraction::inputFile

private

Directory/Name of input file in case you want to read interactions from file.

Definition at line 93 of file NuclearInteraction.cc.

const std::vector<int> fastsim::NuclearInteraction::K0LsID = {130, 310}

private

PdgID of K0.

Definition at line 334 of file NuclearInteraction.cc.

const std::vector<int> fastsim::NuclearInteraction::KminussesID = {-321}

private

PdgID of K-.

Definition at line 336 of file NuclearInteraction.cc.

const std::vector<int> fastsim::NuclearInteraction::KplussesID = {321}

private

PdgID of K+.

Definition at line 335 of file NuclearInteraction.cc.

unsigned fastsim::NuclearInteraction::myOutputBuffer

private

Needed to save interactions to file.

Definition at line 113 of file NuclearInteraction.cc.

std::ofstream fastsim::NuclearInteraction::myOutputFile

private

Output file to save interactions.

Definition at line 112 of file NuclearInteraction.cc.

const std::vector<int> fastsim::NuclearInteraction::neutronsID = {2112, 3122, 3112, 3312, 3322, 3334, -3334}

private

PdgID of neutrons.

Definition at line 332 of file NuclearInteraction.cc.

const std::vector<int> fastsim::NuclearInteraction::PiminussesID = {-211}

private

PdgID of pi-.

Definition at line 338 of file NuclearInteraction.cc.

const std::vector<int> fastsim::NuclearInteraction::PiplussesID = {211}

private

PdgID of pt+.

Definition at line 337 of file NuclearInteraction.cc.

const std::vector<int> fastsim::NuclearInteraction::protonsID = {2212, 3222, -101, -102, -103, -104}

private

PdgID of protons.

Definition at line 330 of file NuclearInteraction.cc.

std::vector<std::vector<TBranch*> > fastsim::NuclearInteraction::theBranches

private

Necessary to read the FullSim interactions.

Definition at line 101 of file NuclearInteraction.cc.

std::vector<std::vector<unsigned> > fastsim::NuclearInteraction::theCurrentEntry

private

Necessary to read the FullSim interactions.

Definition at line 103 of file NuclearInteraction.cc.

std::vector<std::vector<unsigned> > fastsim::NuclearInteraction::theCurrentInteraction

private

Necessary to read the FullSim interactions.

Definition at line 104 of file NuclearInteraction.cc.

double fastsim::NuclearInteraction::theDistCut

private

Cut on deltaR for the FastSim Tracking (ClosestChargedDaughter algorithm)

Definition at line 91 of file NuclearInteraction.cc.

TFile* fastsim::NuclearInteraction::theFile = nullptr

private

Necessary to read the FullSim interactions.

Definition at line 99 of file NuclearInteraction.cc.

std::vector<std::vector<std::string> > fastsim::NuclearInteraction::theFileNames

private

Necessary to read the FullSim interactions.

Definition at line 107 of file NuclearInteraction.cc.

std::vector<std::vector<double> > fastsim::NuclearInteraction::theHadronCM

private

Necessary to read the FullSim interactions.

Definition at line 108 of file NuclearInteraction.cc.

const std::vector<double> fastsim::NuclearInteraction::theHadronEN

private

Initial value:

= {
        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}

Filled into 'theRatiosMap' (evolution of the interaction lengths with energy)

Definition at line 127 of file NuclearInteraction.cc.

double fastsim::NuclearInteraction::theHadronEnergy

private

Minimum energy for nuclear interaction.

Definition at line 92 of file NuclearInteraction.cc.

const std::vector<int> fastsim::NuclearInteraction::theHadronID = {211, -211, 130, 321, -321, 2212, -2212, 2112, -2112}

private

ID of the hadrons.

Definition at line 133 of file NuclearInteraction.cc.

const std::vector<double> fastsim::NuclearInteraction::theHadronMA

private

Initial value:

= {
        0.13957, 0.13957, 0.497648, 0.493677, 0.493677, 0.93827, 0.93827, 0.939565, 0.939565}

Masses of the hadrons.

Definition at line 138 of file NuclearInteraction.cc.

const std::vector<std::string> fastsim::NuclearInteraction::theHadronNA

private

Initial value:

= {
        "piplus", "piminus", "K0L", "Kplus", "Kminus", "p", "pbar", "n", "nbar"}

Names of the hadrons.

Definition at line 135 of file NuclearInteraction.cc.

const std::vector<double> fastsim::NuclearInteraction::theHadronPMin = {0.7, 0.0, 1.0, 1.0, 0.0, 1.1, 0.0, 1.1, 0.0}

private

Smallest momentum for inelastic interactions.

Definition at line 141 of file NuclearInteraction.cc.

std::map< int, int > fastsim::NuclearInteraction::theIDMap

staticprivate

Build the ID map (i.e., what is to be considered as a proton, etc...)

Definition at line 124 of file NuclearInteraction.cc.

const std::vector<double> fastsim::NuclearInteraction::theLengthRatio

private

Initial value:

= {
                                                0.2257,
                                                0.2294,
                                                0.3042,
                                                0.2591,
                                                0.2854,
                                                0.3101,
                                                0.5216,
                                                0.3668,
                                                0.4898}

Inelastic interaction length at p(Hadron) = 5 GeV/c (relativ to radionLength of material)

Definition at line 143 of file NuclearInteraction.cc.

std::vector<std::vector<NUEvent*> > fastsim::NuclearInteraction::theNUEvents

private

Necessary to read the FullSim interactions.

Definition at line 102 of file NuclearInteraction.cc.

std::vector<std::vector<unsigned> > fastsim::NuclearInteraction::theNumberOfEntries

private

Necessary to read the FullSim interactions.

Definition at line 105 of file NuclearInteraction.cc.

std::vector<std::vector<unsigned> > fastsim::NuclearInteraction::theNumberOfInteractions

private

Necessary to read the FullSim interactions.

Definition at line 106 of file NuclearInteraction.cc.

const std::vector<double> fastsim::NuclearInteraction::theRatios

private

Filled into 'theRatiosMap' (evolution of the interaction lengths with energy)

Definition at line 154 of file NuclearInteraction.cc.

std::vector< std::vector< double > > fastsim::NuclearInteraction::theRatiosMap

staticprivate

The evolution of the interaction lengths with energy.

Definition at line 122 of file NuclearInteraction.cc.

std::vector<std::vector<TTree*> > fastsim::NuclearInteraction::theTrees

private

Necessary to read the FullSim interactions.

Definition at line 100 of file NuclearInteraction.cc.