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

◆ NuclearInteraction()

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

Constructor.

Definition at line 348 of file NuclearInteraction.cc.

References antineutronsID, antiprotonsID, looper::cfg, gather_cfg::cout, currentValuesWereSet, Exception, corrVsCorr::filename, contentValuesFiles::fullPath, edm::FileInPath::fullPath(), mps_fire::i, l1ctLayer2EG_cff::id, ien4, fastsim::initializeOnce, inputFile, dqmiolumiharvest::j, K0LsID, KminussesID, KplussesID, myOutputBuffer, myOutputFile, neutronsID, PiminussesID, PiplussesID, protonsID, read(), mathSSE::sqrt(), AlCaHLTBitMon_QueryRunRegistry::string, theBranches, theCurrentEntry, theCurrentInteraction, theDistCut, theFile, theFileNames, theHadronCM, theHadronEN, theHadronEnergy, theHadronID, theHadronMA, theHadronNA, theIDMap, theNUEvents, theNumberOfEntries, theNumberOfInteractions, theRatios, theRatiosMap, theTrees, and to_string().

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

◆ ~NuclearInteraction()

fastsim::NuclearInteraction::~NuclearInteraction ( )

override

Default destructor.

Definition at line 494 of file NuclearInteraction.cc.

References cuy::save, and makeListRunsInFiles::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

◆ index()

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

private

Return a hashed index for a given particle ID.

Definition at line 897 of file NuclearInteraction.cc.

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

◆ interact()

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 513 of file NuclearInteraction.cc.

References A, funct::abs(), fastsim::Particle::charge(), hcalRecHitTable_cff::energy, JetChargeProducer_cfi::exp, NUEvent::NUInteraction::first, RandomEngineAndDistribution::flatShoot(), NUEvent::NUParticle::id, NUEvent::NUInteraction::last, dqm-mbProfile::log, M_PI, NUEvent::NUParticle::mass, fastsim::Particle::momentum(), EgammaValidation_cff::pdgId, fastsim::Particle::pdgId(), ALCARECOTkAlMinBias_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 HLT_2023v12_cff::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.);
     }
   }
 }

◆ orthogonal()

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

private

Return an orthogonal vector.

Definition at line 905 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.);
 }

◆ read()

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 840 of file NuclearInteraction.cc.

References makeListRunsInFiles::inputFile, mysort::results, findQualityFiles::size, and edm_modernize_messagelogger::stat.

Referenced by edmIntegrityCheck.PublishToFileSystem::get(), and NuclearInteraction().

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

◆ save()

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 795 of file NuclearInteraction.cc.

                                      {
   // 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

◆ antineutronsID

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.

Public Member Functions

Private Member Functions

Private Attributes

Static Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ NuclearInteraction()

◆ ~NuclearInteraction()

Member Function Documentation

◆ index()

◆ interact()

◆ orthogonal()

◆ read()

◆ save()

Member Data Documentation

◆ antineutronsID

◆ antiprotonsID

◆ currentValuesWereSet

◆ ien4

◆ inputFile

◆ K0LsID

◆ KminussesID

◆ KplussesID

◆ myOutputBuffer

◆ myOutputFile

◆ neutronsID

◆ PiminussesID

◆ PiplussesID

◆ protonsID

◆ theBranches

◆ theCurrentEntry

◆ theCurrentInteraction

◆ theDistCut

◆ theFile

◆ theFileNames

◆ theHadronCM

◆ theHadronEN

◆ theHadronEnergy

◆ theHadronID

◆ theHadronMA

◆ theHadronNA

◆ theHadronPMin

◆ theIDMap

◆ theLengthRatio

◆ theNUEvents

◆ theNumberOfEntries

◆ theNumberOfInteractions

◆ theRatios

◆ theRatiosMap

◆ theTrees