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::ProducerBase &producer) 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 = 0
	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 = {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) 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 = {0.13957, 0.13957, 0.497648, 0.493677, 0.493677, 0.93827, 0.93827, 0.939565, 0.939565}
	Masses of the hadrons. More...
const std::vector< std::string >	theHadronNA = {"piplus", "piminus", "K0L", "Kplus", "Kminus", "p", "pbar", "n", "nbar"}
	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 61 of file NuclearInteraction.cc.

Constructor & Destructor Documentation

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

Constructor.

Definition at line 210 of file NuclearInteraction.cc.

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

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

References myOutputFile, save(), theFile, and theNUEvents.

                                               {
    // 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 770 of file NuclearInteraction.cc.

References theHadronID.

Referenced by BeautifulSoup.PageElement::insert(), and interact().

{   
    // 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 369 of file NuclearInteraction.cc.

References patCaloMETCorrections_cff::A, funct::abs(), fastsim::Particle::charge(), currentValuesWereSet, JetChargeProducer_cfi::exp, NUEvent::NUInteraction::first, RandomEngineAndDistribution::flatShoot(), fastsim::SimplifiedGeometry::getNuclearInteractionThicknessFactor(), fastsim::SimplifiedGeometry::getThickness(), NUEvent::NUParticle::id, ien4, index(), NUEvent::NUInteraction::last, cmsBatch::log, M_PI, NUEvent::NUParticle::mass, fastsim::Particle::momentum(), orthogonal(), common_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(), theCurrentEntry, theCurrentInteraction, theDistCut, theHadronCM, theHadronEN, theHadronEnergy, theHadronNA, theHadronPMin, theIDMap, theLengthRatio, theNUEvents, theNumberOfEntries, theNumberOfInteractions, theRatiosMap, theta(), theTrees, 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 779 of file NuclearInteraction.cc.

References DEFINE_EDM_PLUGIN, x, y, and z.

Referenced by interact().

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

References findQualityFiles::size, trackingPlots::stat, theCurrentEntry, and theCurrentInteraction.

Referenced by edmIntegrityCheck.PublishToFileSystem::get(), Vispa.Plugins.EdmBrowser.EdmDataAccessor.EdmDataAccessor::goto(), NuclearInteraction(), and Vispa.Plugins.EdmBrowser.EdmDataAccessor.EdmDataAccessor::setFilterBranches().

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

References myOutputBuffer, myOutputFile, findQualityFiles::size, theCurrentEntry, and theCurrentInteraction.

Referenced by Vispa.Main.TabController.TabController::allowClose(), Vispa.Main.TabController.TabController::checkModificationTimestamp(), and ~NuclearInteraction().

{
    // 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 194 of file NuclearInteraction.cc.

Referenced by NuclearInteraction().

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

private

PdgID of anti-protons.

Definition at line 192 of file NuclearInteraction.cc.

Referenced by NuclearInteraction().

bool fastsim::NuclearInteraction::currentValuesWereSet

private

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

Definition at line 117 of file NuclearInteraction.cc.

Referenced by interact(), and NuclearInteraction().

unsigned fastsim::NuclearInteraction::ien4

private

Find the index for which EN = 4.

Definition at line 112 of file NuclearInteraction.cc.

Referenced by interact(), and NuclearInteraction().

std::string fastsim::NuclearInteraction::inputFile

private

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

Definition at line 94 of file NuclearInteraction.cc.

Referenced by NuclearInteraction().

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

private

PdgID of K0.

Definition at line 195 of file NuclearInteraction.cc.

Referenced by NuclearInteraction().

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

private

PdgID of K-.

Definition at line 197 of file NuclearInteraction.cc.

Referenced by NuclearInteraction().

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

private

PdgID of K+.

Definition at line 196 of file NuclearInteraction.cc.

Referenced by NuclearInteraction().

unsigned fastsim::NuclearInteraction::myOutputBuffer

private

Needed to save interactions to file.

Definition at line 115 of file NuclearInteraction.cc.

Referenced by NuclearInteraction(), and save().

std::ofstream fastsim::NuclearInteraction::myOutputFile

private

Output file to save interactions.

Definition at line 114 of file NuclearInteraction.cc.

Referenced by NuclearInteraction(), save(), and ~NuclearInteraction().

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

private

PdgID of neutrons.

Definition at line 193 of file NuclearInteraction.cc.

Referenced by NuclearInteraction().

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

private

PdgID of pi-.

Definition at line 199 of file NuclearInteraction.cc.

Referenced by NuclearInteraction().

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

private

PdgID of pt+.

Definition at line 198 of file NuclearInteraction.cc.

Referenced by NuclearInteraction().

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

private

PdgID of protons.

Definition at line 191 of file NuclearInteraction.cc.

Referenced by NuclearInteraction().

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

private

Necessary to read the FullSim interactions.

Definition at line 103 of file NuclearInteraction.cc.

Referenced by NuclearInteraction().

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

private

Necessary to read the FullSim interactions.

Definition at line 105 of file NuclearInteraction.cc.

Referenced by interact(), NuclearInteraction(), read(), and save().

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

private

Necessary to read the FullSim interactions.

Definition at line 106 of file NuclearInteraction.cc.

Referenced by interact(), NuclearInteraction(), read(), and save().

double fastsim::NuclearInteraction::theDistCut

private

Cut on deltaR for the FastSim Tracking (ClosestChargedDaughter algorithm)

Definition at line 92 of file NuclearInteraction.cc.

Referenced by interact(), and NuclearInteraction().

TFile* fastsim::NuclearInteraction::theFile = 0

private

Necessary to read the FullSim interactions.

Definition at line 101 of file NuclearInteraction.cc.

Referenced by NuclearInteraction(), and ~NuclearInteraction().

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

private

Necessary to read the FullSim interactions.

Definition at line 109 of file NuclearInteraction.cc.

Referenced by NuclearInteraction().

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

private

Necessary to read the FullSim interactions.

Definition at line 110 of file NuclearInteraction.cc.

Referenced by interact(), and NuclearInteraction().

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

private

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

Definition at line 129 of file NuclearInteraction.cc.

Referenced by interact(), and NuclearInteraction().

double fastsim::NuclearInteraction::theHadronEnergy

private

Minimum energy for nuclear interaction.

Definition at line 93 of file NuclearInteraction.cc.

Referenced by interact(), and NuclearInteraction().

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

private

ID of the hadrons.

Definition at line 134 of file NuclearInteraction.cc.

Referenced by index(), and NuclearInteraction().

const std::vector<double> fastsim::NuclearInteraction::theHadronMA = {0.13957, 0.13957, 0.497648, 0.493677, 0.493677, 0.93827, 0.93827, 0.939565, 0.939565}

private

Masses of the hadrons.

Definition at line 138 of file NuclearInteraction.cc.

Referenced by NuclearInteraction().

const std::vector<std::string> fastsim::NuclearInteraction::theHadronNA = {"piplus", "piminus", "K0L", "Kplus", "Kminus", "p", "pbar", "n", "nbar"}

private

Names of the hadrons.

Definition at line 136 of file NuclearInteraction.cc.

Referenced by interact(), and NuclearInteraction().

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

Referenced by interact().

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

Referenced by interact().

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

Referenced by interact().

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

private

Necessary to read the FullSim interactions.

Definition at line 104 of file NuclearInteraction.cc.

Referenced by interact(), NuclearInteraction(), and ~NuclearInteraction().

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

private

Necessary to read the FullSim interactions.

Definition at line 107 of file NuclearInteraction.cc.

Referenced by interact(), and NuclearInteraction().

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

private

Necessary to read the FullSim interactions.

Definition at line 108 of file NuclearInteraction.cc.

Referenced by interact(), and NuclearInteraction().

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

private

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

Definition at line 148 of file NuclearInteraction.cc.

Referenced by NuclearInteraction().

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

staticprivate

The evolution of the interaction lengths with energy.

Definition at line 124 of file NuclearInteraction.cc.

Referenced by interact().

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

private

Necessary to read the FullSim interactions.

Definition at line 102 of file NuclearInteraction.cc.

Referenced by interact(), and NuclearInteraction().