Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:26:47

0001 // system headers
0002 #include <cmath>
0003 #include <memory>
0004 #include <vector>
0005 #include <map>
0006 #include <string>
0007 #include <fstream>
0008 #include <iostream>
0009 #include <sys/stat.h>
0010 #include <cmath>
0011 #include <TFile.h>
0012 #include <TTree.h>
0013 #include <TROOT.h>
0014 #include <Math/AxisAngle.h>
0015 #include "Math/Boost.h"
0016 #include "DataFormats/Math/interface/Vector3D.h"
0017 
0018 // Framework Headers
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/Utilities/interface/Exception.h"
0022 #include "FWCore/ParameterSet/interface/FileInPath.h"
0023 
0024 #include "FWCore/Utilities/interface/thread_safety_macros.h"
0025 
0026 // Fast Sim headers
0027 #include "FastSimulation/SimplifiedGeometryPropagator/interface/Particle.h"
0028 #include "FastSimulation/SimplifiedGeometryPropagator/interface/SimplifiedGeometry.h"
0029 #include "FastSimulation/SimplifiedGeometryPropagator/interface/InteractionModelFactory.h"
0030 #include "FastSimulation/SimplifiedGeometryPropagator/interface/InteractionModel.h"
0031 #include "FastSimulation/SimplifiedGeometryPropagator/interface/Constants.h"
0032 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0033 #include "FastSimDataFormats/NuclearInteractions/interface/NUEvent.h"
0034 
0035 // Math
0036 #include "DataFormats/Math/interface/LorentzVector.h"
0037 
0038 ///////////////////////////////////////////////
0039 // Author: Patrick Janot
0040 // Date: 25-Jan-2007
0041 //
0042 // Revision: Class structure modified to match SimplifiedGeometryPropagator
0043 //           Unified the way, the closest charged daughter is defined
0044 //           S. Kurz, 29 May 2017
0045 //////////////////////////////////////////////////////////
0046 
0047 // TODO: save() function called in destructer. Does that actually make sense?
0048 
0049 typedef math::XYZVector XYZVector;
0050 typedef math::XYZTLorentzVector XYZTLorentzVector;
0051 
0052 namespace fastsim {
0053   //! Implementation of nuclear interactions of hadrons in the tracker layers (based on fully simulated interactions).
0054   /*!
0055         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.
0056         Also, another implementation of nuclear interactions can be used that is based on G4 (NuclearInteractionFTF).
0057     */
0058   class NuclearInteraction : public InteractionModel {
0059   public:
0060     //! Constructor.
0061     NuclearInteraction(const std::string& name, const edm::ParameterSet& cfg);
0062 
0063     //! Default destructor.
0064     ~NuclearInteraction() override;
0065 
0066     //! Perform the interaction.
0067     /*!
0068             \param particle The particle that interacts with the matter.
0069             \param layer The detector layer that interacts with the particle.
0070             \param secondaries Particles that are produced in the interaction (if any).
0071             \param random The Random Engine.
0072         */
0073     void interact(fastsim::Particle& particle,
0074                   const SimplifiedGeometry& layer,
0075                   std::vector<std::unique_ptr<fastsim::Particle> >& secondaries,
0076                   const RandomEngineAndDistribution& random) override;
0077 
0078   private:
0079     //! Return a hashed index for a given particle ID
0080     unsigned index(int thePid);
0081 
0082     //! Return an orthogonal vector.
0083     XYZVector orthogonal(const XYZVector& aVector) const;
0084 
0085     //! Save the nuclear interactions to a file, so you can reproduce the events (e.g. in case of a crash).
0086     void save();
0087 
0088     //! Read the nuclear interactions from a file, so you can reproduce the events (e.g. in case of a crash).
0089     bool read(std::string inputFile);
0090 
0091     double theDistCut;       //!< Cut on deltaR for the FastSim Tracking (ClosestChargedDaughter algorithm)
0092     double theHadronEnergy;  //!< Minimum energy for nuclear interaction
0093     std::string inputFile;   //!< Directory/Name of input file in case you want to read interactions from file
0094 
0095     //////////
0096     // Read/Save nuclear interactions from FullSim
0097     ///////////////
0098 
0099     TFile* theFile = nullptr;                                     //!< Necessary to read the FullSim interactions
0100     std::vector<std::vector<TTree*> > theTrees;                   //!< Necessary to read the FullSim interactions
0101     std::vector<std::vector<TBranch*> > theBranches;              //!< Necessary to read the FullSim interactions
0102     std::vector<std::vector<NUEvent*> > theNUEvents;              //!< Necessary to read the FullSim interactions
0103     std::vector<std::vector<unsigned> > theCurrentEntry;          //!< Necessary to read the FullSim interactions
0104     std::vector<std::vector<unsigned> > theCurrentInteraction;    //!< Necessary to read the FullSim interactions
0105     std::vector<std::vector<unsigned> > theNumberOfEntries;       //!< Necessary to read the FullSim interactions
0106     std::vector<std::vector<unsigned> > theNumberOfInteractions;  //!< Necessary to read the FullSim interactions
0107     std::vector<std::vector<std::string> > theFileNames;          //!< Necessary to read the FullSim interactions
0108     std::vector<std::vector<double> > theHadronCM;                //!< Necessary to read the FullSim interactions
0109 
0110     unsigned ien4;  //!< Find the index for which EN = 4
0111 
0112     std::ofstream myOutputFile;  //!< Output file to save interactions
0113     unsigned myOutputBuffer;     //!< Needed to save interactions to file
0114 
0115     bool currentValuesWereSet;  //!< Read data from file that was created in a previous run
0116 
0117     //////////
0118     // Properties of the Hadrons
0119     ///////////////
0120 
0121     //! The evolution of the interaction lengths with energy
0122     static std::vector<std::vector<double> > theRatiosMap;
0123     //! Build the ID map (i.e., what is to be considered as a proton, etc...)
0124     static std::map<int, int> theIDMap;
0125 
0126     //! Filled into 'theRatiosMap' (evolution of the interaction lengths with energy)
0127     const std::vector<double> theHadronEN = {
0128         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};
0129 
0130     // Use same order for all vectors below (Properties for "piplus", "piminus", "K0L", "Kplus", "Kminus", "p", "pbar", "n", "nbar")
0131 
0132     //! ID of the hadrons
0133     const std::vector<int> theHadronID = {211, -211, 130, 321, -321, 2212, -2212, 2112, -2112};
0134     //! Names of the hadrons
0135     const std::vector<std::string> theHadronNA = {
0136         "piplus", "piminus", "K0L", "Kplus", "Kminus", "p", "pbar", "n", "nbar"};
0137     //! Masses of the hadrons
0138     const std::vector<double> theHadronMA = {
0139         0.13957, 0.13957, 0.497648, 0.493677, 0.493677, 0.93827, 0.93827, 0.939565, 0.939565};
0140     //! Smallest momentum for inelastic interactions
0141     const std::vector<double> theHadronPMin = {0.7, 0.0, 1.0, 1.0, 0.0, 1.1, 0.0, 1.1, 0.0};
0142     //! Inelastic interaction length at p(Hadron) = 5 GeV/c (relativ to radionLength of material)
0143     const std::vector<double> theLengthRatio = {//pi+      pi-    K0L      K+      K-      p      pbar     n      nbar
0144                                                 0.2257,
0145                                                 0.2294,
0146                                                 0.3042,
0147                                                 0.2591,
0148                                                 0.2854,
0149                                                 0.3101,
0150                                                 0.5216,
0151                                                 0.3668,
0152                                                 0.4898};
0153     //! Filled into 'theRatiosMap' (evolution of the interaction lengths with energy)
0154     const std::vector<double> theRatios = {// pi+ (211)
0155                                            0.031390573,
0156                                            0.531842852,
0157                                            0.819614219,
0158                                            0.951251711,
0159                                            0.986382750,
0160                                            1.000000000,
0161                                            0.985087033,
0162                                            0.982996773,
0163                                            0.990832192,
0164                                            0.992237923,
0165                                            0.994841580,
0166                                            0.973816742,
0167                                            0.967264815,
0168                                            0.971714258,
0169                                            0.969122824,
0170                                            0.978681792,
0171                                            0.977312732,
0172                                            0.984255819,
0173                                            // pi- (-211)
0174                                            0.035326512,
0175                                            0.577356403,
0176                                            0.857118809,
0177                                            0.965683504,
0178                                            0.989659360,
0179                                            1.000000000,
0180                                            0.989599240,
0181                                            0.980665408,
0182                                            0.988384816,
0183                                            0.981038152,
0184                                            0.975002104,
0185                                            0.959996152,
0186                                            0.953310808,
0187                                            0.954705592,
0188                                            0.957615400,
0189                                            0.961150456,
0190                                            0.965022184,
0191                                            0.960573304,
0192                                            // K0L (130)
0193                                            0.000000000,
0194                                            0.370261189,
0195                                            0.649793096,
0196                                            0.734342408,
0197                                            0.749079499,
0198                                            0.753360057,
0199                                            0.755790543,
0200                                            0.755872164,
0201                                            0.751337674,
0202                                            0.746685288,
0203                                            0.747519634,
0204                                            0.739357554,
0205                                            0.735004444,
0206                                            0.803039922,
0207                                            0.832749896,
0208                                            0.890900187,
0209                                            0.936734805,
0210                                            1.000000000,
0211                                            // K+ (321)
0212                                            0.000000000,
0213                                            0.175571717,
0214                                            0.391683394,
0215                                            0.528946472,
0216                                            0.572818635,
0217                                            0.614210280,
0218                                            0.644125538,
0219                                            0.670304050,
0220                                            0.685144573,
0221                                            0.702870161,
0222                                            0.714708513,
0223                                            0.730805263,
0224                                            0.777711536,
0225                                            0.831090576,
0226                                            0.869267129,
0227                                            0.915747562,
0228                                            0.953370523,
0229                                            1.000000000,
0230                                            // K- (-321)
0231                                            0.000000000,
0232                                            0.365353210,
0233                                            0.611663677,
0234                                            0.715315908,
0235                                            0.733498956,
0236                                            0.738361302,
0237                                            0.745253654,
0238                                            0.751459671,
0239                                            0.750628335,
0240                                            0.746442657,
0241                                            0.750850669,
0242                                            0.744895986,
0243                                            0.735093960,
0244                                            0.791663444,
0245                                            0.828609543,
0246                                            0.889993040,
0247                                            0.940897842,
0248                                            1.000000000,
0249                                            // proton (2212)
0250                                            0.000000000,
0251                                            0.042849136,
0252                                            0.459103223,
0253                                            0.666165343,
0254                                            0.787930873,
0255                                            0.890397011,
0256                                            0.920999533,
0257                                            0.937832788,
0258                                            0.950920131,
0259                                            0.966595049,
0260                                            0.979542270,
0261                                            0.988061653,
0262                                            0.983260159,
0263                                            0.988958431,
0264                                            0.991723494,
0265                                            0.995273237,
0266                                            1.000000000,
0267                                            0.999962634,
0268                                            // anti-proton (-2212)
0269                                            1.000000000,
0270                                            0.849956907,
0271                                            0.775625988,
0272                                            0.802018230,
0273                                            0.816207485,
0274                                            0.785899785,
0275                                            0.754998487,
0276                                            0.728977244,
0277                                            0.710010673,
0278                                            0.670890339,
0279                                            0.665627872,
0280                                            0.652682888,
0281                                            0.613334247,
0282                                            0.647534574,
0283                                            0.667910938,
0284                                            0.689919693,
0285                                            0.709200185,
0286                                            0.724199928,
0287                                            // neutron (2112)
0288                                            0.000000000,
0289                                            0.059216484,
0290                                            0.437844536,
0291                                            0.610370629,
0292                                            0.702090648,
0293                                            0.780076890,
0294                                            0.802143073,
0295                                            0.819570432,
0296                                            0.825829666,
0297                                            0.840079750,
0298                                            0.838435509,
0299                                            0.837529986,
0300                                            0.835687165,
0301                                            0.885205014,
0302                                            0.912450156,
0303                                            0.951451221,
0304                                            0.973215562,
0305                                            1.000000000,
0306                                            // anti-neutron
0307                                            1.000000000,
0308                                            0.849573257,
0309                                            0.756479495,
0310                                            0.787147094,
0311                                            0.804572414,
0312                                            0.791806302,
0313                                            0.760234588,
0314                                            0.741109531,
0315                                            0.724118186,
0316                                            0.692829761,
0317                                            0.688465897,
0318                                            0.671806061,
0319                                            0.636461171,
0320                                            0.675314029,
0321                                            0.699134460,
0322                                            0.724305037,
0323                                            0.742556115,
0324                                            0.758504713};
0325 
0326     //////////
0327     // Used to build the ID map 'theIDMap' (i.e., what is to be considered as a proton, etc...)
0328     ///////////////
0329 
0330     const std::vector<int> protonsID = {2212, 3222, -101, -102, -103, -104};          //!< PdgID of protons
0331     const std::vector<int> antiprotonsID = {-2212, -3222};                            //!< PdgID of anti-protons
0332     const std::vector<int> neutronsID = {2112, 3122, 3112, 3312, 3322, 3334, -3334};  //!< PdgID of neutrons
0333     const std::vector<int> antineutronsID = {-2112, -3122, -3112, -3312, -3322};      //!< PdgID of anti-neutrons
0334     const std::vector<int> K0LsID = {130, 310};                                       //!< PdgID of K0
0335     const std::vector<int> KplussesID = {321};                                        //!< PdgID of K+
0336     const std::vector<int> KminussesID = {-321};                                      //!< PdgID of K-
0337     const std::vector<int> PiplussesID = {211};                                       //!< PdgID of pt+
0338     const std::vector<int> PiminussesID = {-211};                                     //!< PdgID of pi-
0339   };
0340 
0341   // TODO: Is this correct?
0342   // Thread safety
0343   static std::once_flag initializeOnce;
0344   CMS_THREAD_GUARD(initializeOnce) std::vector<std::vector<double> > NuclearInteraction::theRatiosMap;
0345   CMS_THREAD_GUARD(initializeOnce) std::map<int, int> NuclearInteraction::theIDMap;
0346 }  // namespace fastsim
0347 
0348 fastsim::NuclearInteraction::NuclearInteraction(const std::string& name, const edm::ParameterSet& cfg)
0349     : fastsim::InteractionModel(name), currentValuesWereSet(false) {
0350   // Full path to FullSim root file
0351   std::string fullPath;
0352 
0353   // Read from config file
0354   theDistCut = cfg.getParameter<double>("distCut");
0355   theHadronEnergy = cfg.getParameter<double>("hadronEnergy");
0356   inputFile = cfg.getUntrackedParameter<std::string>("inputFile", "");
0357 
0358   // The evolution of the interaction lengths with energy
0359   // initialize once for all possible instances
0360   std::call_once(initializeOnce, [this]() {
0361     theRatiosMap.resize(theHadronID.size());
0362     for (unsigned i = 0; i < theHadronID.size(); ++i) {
0363       for (unsigned j = 0; j < theHadronEN.size(); ++j) {
0364         theRatiosMap[i].push_back(theRatios[i * theHadronEN.size() + j]);
0365       }
0366     }
0367 
0368     // Build the ID map (i.e., what is to be considered as a proton, etc...)
0369     // Protons
0370     for (const auto& id : protonsID)
0371       theIDMap[id] = 2212;
0372     // Anti-Protons
0373     for (const auto& id : antiprotonsID)
0374       theIDMap[id] = -2212;
0375     // Neutrons
0376     for (const auto& id : neutronsID)
0377       theIDMap[id] = 2112;
0378     // Anti-Neutrons
0379     for (const auto& id : antineutronsID)
0380       theIDMap[id] = -2112;
0381     // K0L's
0382     for (const auto& id : K0LsID)
0383       theIDMap[id] = 130;
0384     // K+'s
0385     for (const auto& id : KplussesID)
0386       theIDMap[id] = 321;
0387     // K-'s
0388     for (const auto& id : KminussesID)
0389       theIDMap[id] = -321;
0390     // pi+'s
0391     for (const auto& id : PiplussesID)
0392       theIDMap[id] = 211;
0393     // pi-'s
0394     for (const auto& id : PiminussesID)
0395       theIDMap[id] = -211;
0396   });
0397 
0398   // Prepare the map of files
0399   // Loop over the particle names
0400   TFile* aVFile = nullptr;
0401   std::vector<TTree*> aVTree(theHadronEN.size());
0402   std::vector<TBranch*> aVBranch(theHadronEN.size());
0403   std::vector<NUEvent*> aVNUEvents(theHadronEN.size());
0404   std::vector<unsigned> aVCurrentEntry(theHadronEN.size());
0405   std::vector<unsigned> aVCurrentInteraction(theHadronEN.size());
0406   std::vector<unsigned> aVNumberOfEntries(theHadronEN.size());
0407   std::vector<unsigned> aVNumberOfInteractions(theHadronEN.size());
0408   std::vector<std::string> aVFileName(theHadronEN.size());
0409   std::vector<double> aVHadronCM(theHadronEN.size());
0410   theTrees.resize(theHadronNA.size());
0411   theBranches.resize(theHadronNA.size());
0412   theNUEvents.resize(theHadronNA.size());
0413   theCurrentEntry.resize(theHadronNA.size());
0414   theCurrentInteraction.resize(theHadronNA.size());
0415   theNumberOfEntries.resize(theHadronNA.size());
0416   theNumberOfInteractions.resize(theHadronNA.size());
0417   theFileNames.resize(theHadronNA.size());
0418   theHadronCM.resize(theHadronNA.size());
0419   theFile = aVFile;
0420   for (unsigned iname = 0; iname < theHadronNA.size(); ++iname) {
0421     theTrees[iname] = aVTree;
0422     theBranches[iname] = aVBranch;
0423     theNUEvents[iname] = aVNUEvents;
0424     theCurrentEntry[iname] = aVCurrentEntry;
0425     theCurrentInteraction[iname] = aVCurrentInteraction;
0426     theNumberOfEntries[iname] = aVNumberOfEntries;
0427     theNumberOfInteractions[iname] = aVNumberOfInteractions;
0428     theFileNames[iname] = aVFileName;
0429     theHadronCM[iname] = aVHadronCM;
0430   }
0431 
0432   // Read the information from a previous run (to keep reproducibility)
0433   currentValuesWereSet = this->read(inputFile);
0434   if (currentValuesWereSet)
0435     std::cout << "***WARNING*** You are reading nuclear-interaction information from the file " << inputFile
0436               << " created in an earlier run." << std::endl;
0437 
0438   // Open the file for saving the information of the current run
0439   myOutputFile.open("NuclearInteractionOutputFile.txt");
0440   myOutputBuffer = 0;
0441 
0442   // Open the root file
0443   edm::FileInPath myDataFile("FastSimulation/MaterialEffects/data/NuclearInteractions.root");
0444   fullPath = myDataFile.fullPath();
0445   theFile = TFile::Open(fullPath.c_str());
0446 
0447   // Open the trees
0448   unsigned fileNb = 0;
0449   for (unsigned iname = 0; iname < theHadronNA.size(); ++iname) {
0450     for (unsigned iene = 0; iene < theHadronEN.size(); ++iene) {
0451       std::ostringstream filename;
0452       double theEne = theHadronEN[iene];
0453       filename << "NuclearInteractionsVal_" << theHadronNA[iname] << "_E" << theEne << ".root";
0454       theFileNames[iname][iene] = filename.str();
0455 
0456       ++fileNb;
0457       std::string treeName = "NuclearInteractions_" + theHadronNA[iname] + "_E" + std::to_string(int(theEne));
0458       theTrees[iname][iene] = (TTree*)theFile->Get(treeName.c_str());
0459 
0460       if (!theTrees[iname][iene])
0461         throw cms::Exception("FastSimulation/MaterialEffects") << "Tree with name " << treeName << " not found ";
0462       theBranches[iname][iene] = theTrees[iname][iene]->GetBranch("nuEvent");
0463       if (!theBranches[iname][iene])
0464         throw cms::Exception("FastSimulation/MaterialEffects")
0465             << "Branch with name nuEvent not found in " << theFileNames[iname][iene];
0466 
0467       theNUEvents[iname][iene] = new NUEvent();
0468       theBranches[iname][iene]->SetAddress(&theNUEvents[iname][iene]);
0469       theNumberOfEntries[iname][iene] = theTrees[iname][iene]->GetEntries();
0470 
0471       if (currentValuesWereSet) {
0472         theTrees[iname][iene]->GetEntry(theCurrentEntry[iname][iene]);
0473         unsigned NInteractions = theNUEvents[iname][iene]->nInteractions();
0474         theNumberOfInteractions[iname][iene] = NInteractions;
0475       }
0476 
0477       // Compute the corresponding cm energies of the nuclear interactions
0478       XYZTLorentzVector Proton(0., 0., 0., 0.986);
0479       XYZTLorentzVector Reference(
0480           0.,
0481           0.,
0482           std::sqrt(theHadronEN[iene] * theHadronEN[iene] - theHadronMA[iname] * theHadronMA[iname]),
0483           theHadronEN[iene]);
0484       theHadronCM[iname][iene] = (Reference + Proton).M();
0485     }
0486   }
0487 
0488   // Find the index for which EN = 4. (or thereabout)
0489   ien4 = 0;
0490   while (theHadronEN[ien4] < 4.0)
0491     ++ien4;
0492 
0493   gROOT->cd();
0494 }
0495 
0496 fastsim::NuclearInteraction::~NuclearInteraction() {
0497   // Close all local files
0498   // Among other things, this allows the TROOT destructor to end up
0499   // without crashing, while trying to close these files from outside
0500   theFile->Close();
0501   delete theFile;
0502 
0503   for (auto& vEvents : theNUEvents) {
0504     for (auto evtPtr : vEvents) {
0505       delete evtPtr;
0506     }
0507   }
0508 
0509   // Save data
0510   save();
0511   // Close the output file
0512   myOutputFile.close();
0513 }
0514 
0515 void fastsim::NuclearInteraction::interact(fastsim::Particle& particle,
0516                                            const SimplifiedGeometry& layer,
0517                                            std::vector<std::unique_ptr<fastsim::Particle> >& secondaries,
0518                                            const RandomEngineAndDistribution& random) {
0519   int pdgId = particle.pdgId();
0520   //
0521   // no valid PDGid
0522   //
0523   if (abs(pdgId) <= 100 || abs(pdgId) >= 1000000) {
0524     return;
0525   }
0526 
0527   double radLengths = layer.getThickness(particle.position(), particle.momentum());
0528   // TEC layers have some fudge factor for nuclear interactions
0529   radLengths *= layer.getNuclearInteractionThicknessFactor();
0530   //
0531   // no material
0532   //
0533   if (radLengths < 1E-10) {
0534     return;
0535   }
0536 
0537   // In case the events are not read from (old) saved file, then pick a random event from FullSim file
0538   if (!currentValuesWereSet) {
0539     currentValuesWereSet = true;
0540     for (unsigned iname = 0; iname < theHadronNA.size(); ++iname) {
0541       for (unsigned iene = 0; iene < theHadronEN.size(); ++iene) {
0542         theCurrentEntry[iname][iene] = (unsigned)(theNumberOfEntries[iname][iene] * random.flatShoot());
0543 
0544         theTrees[iname][iene]->GetEntry(theCurrentEntry[iname][iene]);
0545         unsigned NInteractions = theNUEvents[iname][iene]->nInteractions();
0546         theNumberOfInteractions[iname][iene] = NInteractions;
0547 
0548         theCurrentInteraction[iname][iene] = (unsigned)(theNumberOfInteractions[iname][iene] * random.flatShoot());
0549       }
0550     }
0551   }
0552 
0553   // Momentum of interacting hadron
0554   double pHadron = std::sqrt(particle.momentum().Vect().Mag2());
0555 
0556   //
0557   // The hadron has not enough momentum to create some relevant final state
0558   //
0559   if (pHadron < theHadronEnergy) {
0560     return;
0561   }
0562 
0563   // The particle type
0564   std::map<int, int>::const_iterator thePit = theIDMap.find(pdgId);
0565   // Use map for unique ID (e.g. proton = {2212, 3222, -101, -102, -103, -104})
0566   int thePid = (thePit != theIDMap.end() ? thePit->second : pdgId);
0567 
0568   // Is this particle type foreseen?
0569   unsigned fPid = abs(thePid);
0570   if (fPid != 211 && fPid != 130 && fPid != 321 && fPid != 2112 && fPid != 2212) {
0571     return;
0572   }
0573 
0574   // The hashed ID
0575   unsigned thePidIndex = index(thePid);
0576   // The inelastic interaction length at p(Hadron) = 5 GeV/c
0577   double theInelasticLength = radLengths * theLengthRatio[thePidIndex];
0578 
0579   // The elastic interaction length
0580   // The baroque parameterization is a fit to Fig. 40.13 of the PDG
0581   double ee = pHadron > 0.6 ? exp(-std::sqrt((pHadron - 0.6) / 1.122)) : exp(std::sqrt((0.6 - pHadron) / 1.122));
0582   double theElasticLength = (0.8753 * ee + 0.15) * theInelasticLength;
0583 
0584   // The total interaction length
0585   double theTotalInteractionLength = theInelasticLength + theElasticLength;
0586 
0587   //
0588   // Probability to interact is dl/L0 (maximum for 4 GeV Hadron)
0589   //
0590   double aNuclInteraction = -std::log(random.flatShoot());
0591   if (aNuclInteraction > theTotalInteractionLength) {
0592     return;
0593   }
0594 
0595   // The elastic part
0596   double elastic = random.flatShoot();
0597   if (elastic < theElasticLength / theTotalInteractionLength) {
0598     // Characteristic scattering angle for the elastic part
0599     // A of silicon
0600     double A = 28.0855;
0601     double theta0 = std::sqrt(3.) / std::pow(A, 1. / 3.) * particle.momentum().mass() / pHadron;
0602 
0603     // Draw an angle with theta/theta0*exp[(-theta/2theta0)**2] shape
0604     double theta = theta0 * std::sqrt(-2. * std::log(random.flatShoot()));
0605     double phi = 2. * M_PI * random.flatShoot();
0606 
0607     // Rotate the particle accordingly
0608     ROOT::Math::AxisAngle rotation1(orthogonal(particle.momentum().Vect()), theta);
0609     ROOT::Math::AxisAngle rotation2(particle.momentum().Vect(), phi);
0610     // Rotate!
0611     XYZVector rotated = rotation2((rotation1(particle.momentum().Vect())));
0612 
0613     // Create a daughter if the kink is large enough
0614     if (std::sin(theta) > theDistCut) {
0615       secondaries.emplace_back(
0616           new fastsim::Particle(pdgId,
0617                                 particle.position(),
0618                                 XYZTLorentzVector(rotated.X(), rotated.Y(), rotated.Z(), particle.momentum().E())));
0619 
0620       // Set the ClosestCharged Daughter thing for tracking
0621       if (particle.charge() != 0) {
0622         secondaries.back()->setMotherDeltaR(particle.momentum());
0623         secondaries.back()->setMotherPdgId(pdgId);
0624         secondaries.back()->setMotherSimTrackIndex(particle.simTrackIndex());
0625       }
0626 
0627       // The particle is destroyed
0628       particle.momentum().SetXYZT(0., 0., 0., 0.);
0629     } else {
0630       // If kink is small enough just rotate particle
0631       particle.momentum().SetXYZT(rotated.X(), rotated.Y(), rotated.Z(), particle.momentum().E());
0632     }
0633     // The inelastic part
0634   } else {
0635     // Avoid multiple map access
0636     const std::vector<double>& aHadronCM = theHadronCM[thePidIndex];
0637     const std::vector<double>& aRatios = theRatiosMap[thePidIndex];
0638     // Find the file with the closest c.m energy
0639     // The target nucleon
0640     XYZTLorentzVector Proton(0., 0., 0., 0.939);
0641     // The current particle
0642     const XYZTLorentzVector& Hadron = (const XYZTLorentzVector&)particle.momentum();
0643     // The smallest momentum for inelastic interactions
0644     double pMin = theHadronPMin[thePidIndex];
0645     // The corresponding smallest four vector
0646     XYZTLorentzVector Hadron0(0., 0., pMin, std::sqrt(pMin * pMin + particle.momentum().M2()));
0647 
0648     // The current centre-of-mass energy
0649     double ecm = (Proton + Hadron).M();
0650 
0651     // Get the files of interest (closest c.m. energies)
0652     unsigned ene1 = 0;
0653     unsigned ene2 = 0;
0654     // The smallest centre-of-mass energy
0655     double ecm1 = (Proton + Hadron0).M();
0656 
0657     double ecm2 = aHadronCM[0];
0658     double ratio1 = 0.;
0659     double ratio2 = aRatios[0];
0660     if (ecm > aHadronCM[0] && ecm < aHadronCM[aHadronCM.size() - 1]) {
0661       for (unsigned ene = 1; ene < aHadronCM.size() && ecm > aHadronCM[ene - 1]; ++ene) {
0662         if (ecm < aHadronCM[ene]) {
0663           ene2 = ene;
0664           ene1 = ene2 - 1;
0665           ecm1 = aHadronCM[ene1];
0666           ecm2 = aHadronCM[ene2];
0667           ratio1 = aRatios[ene1];
0668           ratio2 = aRatios[ene2];
0669         }
0670       }
0671     } else if (ecm > aHadronCM[aHadronCM.size() - 1]) {
0672       ene1 = aHadronCM.size() - 1;
0673       ene2 = aHadronCM.size() - 2;
0674       ecm1 = aHadronCM[ene1];
0675       ecm2 = aHadronCM[ene2];
0676       ratio1 = aRatios[ene2];
0677       ratio2 = aRatios[ene2];
0678     }
0679 
0680     // The inelastic part of the cross section depends cm energy
0681     double slope = (std::log10(ecm) - std::log10(ecm1)) / (std::log10(ecm2) - std::log10(ecm1));
0682     double inelastic = ratio1 + (ratio2 - ratio1) * slope;
0683     double inelastic4 = pHadron < 4. ? aRatios[ien4] : 1.;
0684 
0685     // Simulate an inelastic interaction
0686     if (elastic > 1. - (inelastic * theInelasticLength) / theTotalInteractionLength) {
0687       // Avoid mutliple map access
0688       std::vector<unsigned>& aCurrentInteraction = theCurrentInteraction[thePidIndex];
0689       std::vector<unsigned>& aNumberOfInteractions = theNumberOfInteractions[thePidIndex];
0690       std::vector<NUEvent*>& aNUEvents = theNUEvents[thePidIndex];
0691 
0692       // Choice of the file to read according the the log10(ecm) distance
0693       // and protection against low momentum proton and neutron that never interacts
0694       // (i.e., empty files)
0695       unsigned ene;
0696       if (random.flatShoot() < slope || aNumberOfInteractions[ene1] == 0) {
0697         ene = ene2;
0698       } else {
0699         ene = ene1;
0700       }
0701 
0702       // The boost characteristics
0703       XYZTLorentzVector theBoost = Proton + Hadron;
0704       theBoost /= theBoost.e();
0705 
0706       // Check we are not either at the end of an interaction bunch
0707       // or at the end of a file
0708       if (aCurrentInteraction[ene] == aNumberOfInteractions[ene]) {
0709         std::vector<unsigned>& aCurrentEntry = theCurrentEntry[thePidIndex];
0710         std::vector<unsigned>& aNumberOfEntries = theNumberOfEntries[thePidIndex];
0711         std::vector<TTree*>& aTrees = theTrees[thePidIndex];
0712         ++aCurrentEntry[ene];
0713 
0714         aCurrentInteraction[ene] = 0;
0715         if (aCurrentEntry[ene] == aNumberOfEntries[ene]) {
0716           aCurrentEntry[ene] = 0;
0717         }
0718 
0719         unsigned myEntry = aCurrentEntry[ene];
0720         aTrees[ene]->GetEntry(myEntry);
0721         aNumberOfInteractions[ene] = aNUEvents[ene]->nInteractions();
0722       }
0723 
0724       // Read the interaction
0725       NUEvent::NUInteraction anInteraction = aNUEvents[ene]->theNUInteractions()[aCurrentInteraction[ene]];
0726 
0727       unsigned firstTrack = anInteraction.first;
0728       unsigned lastTrack = anInteraction.last;
0729 
0730       // Some rotation around the boost axis, for more randomness
0731       XYZVector theAxis = theBoost.Vect().Unit();
0732       double theAngle = random.flatShoot() * 2. * M_PI;
0733       ROOT::Math::AxisAngle axisRotation(theAxis, theAngle);
0734       ROOT::Math::Boost axisBoost(theBoost.x(), theBoost.y(), theBoost.z());
0735 
0736       // A rotation to bring the particles back to the Hadron direction
0737       XYZVector zAxis(0., 0., 1.);
0738       XYZVector orthAxis = (zAxis.Cross(theBoost.Vect())).Unit();
0739       double orthAngle = acos(theBoost.Vect().Unit().Z());
0740       ROOT::Math::AxisAngle orthRotation(orthAxis, orthAngle);
0741 
0742       // Loop on the nuclear interaction products
0743       for (unsigned iTrack = firstTrack; iTrack <= lastTrack; ++iTrack) {
0744         NUEvent::NUParticle aParticle = aNUEvents[ene]->theNUParticles()[iTrack];
0745 
0746         // Add a RawParticle with the proper energy in the c.m. frame of
0747         // the nuclear interaction
0748         double energy = std::sqrt(aParticle.px * aParticle.px + aParticle.py * aParticle.py +
0749                                   aParticle.pz * aParticle.pz + aParticle.mass * aParticle.mass / (ecm * ecm));
0750 
0751         XYZTLorentzVector daugtherMomentum(aParticle.px * ecm, aParticle.py * ecm, aParticle.pz * ecm, energy * ecm);
0752 
0753         // Rotate to the collision axis
0754         XYZVector rotated = orthRotation(daugtherMomentum.Vect());
0755         // Rotate around the boost axis for more randomness
0756         rotated = axisRotation(rotated);
0757 
0758         // Rotated the daughter
0759         daugtherMomentum.SetXYZT(rotated.X(), rotated.Y(), rotated.Z(), daugtherMomentum.E());
0760 
0761         // Boost it in the lab frame
0762         daugtherMomentum = axisBoost(daugtherMomentum);
0763 
0764         // Create secondary
0765         secondaries.emplace_back(new fastsim::Particle(aParticle.id, particle.position(), daugtherMomentum));
0766 
0767         // The closestCharged Daughter thing for tracking
0768         // BUT: 'aParticle' also has to be charged, only then the mother should be set
0769         // Unfortunately, NUEvent::NUParticle does not contain any info about the charge
0770         // Did some tests and effect is absolutely negligible!
0771         if (particle.charge() != 0) {
0772           secondaries.back()->setMotherDeltaR(particle.momentum());
0773           secondaries.back()->setMotherPdgId(pdgId);
0774           secondaries.back()->setMotherSimTrackIndex(particle.simTrackIndex());
0775         }
0776       }
0777 
0778       // The particle is destroyed
0779       particle.momentum().SetXYZT(0., 0., 0., 0.);
0780 
0781       // This is a note from previous version of code but I don't understand it:
0782       // ERROR The way this loops through the events breaks
0783       // replay. Which events are retrieved depends on
0784       // which previous events were processed.
0785 
0786       // Increment for next time
0787       ++aCurrentInteraction[ene];
0788 
0789       // Simulate a stopping hadron (low momentum)
0790     } else if (pHadron < 4. && elastic > 1. - (inelastic4 * theInelasticLength) / theTotalInteractionLength) {
0791       // The particle is destroyed
0792       particle.momentum().SetXYZT(0., 0., 0., 0.);
0793     }
0794   }
0795 }
0796 
0797 void fastsim::NuclearInteraction::save() {
0798   // Size of buffer
0799   ++myOutputBuffer;
0800 
0801   // Periodically close the current file and open a new one
0802   if (myOutputBuffer / 1000 * 1000 == myOutputBuffer) {
0803     myOutputFile.close();
0804     myOutputFile.open("NuclearInteractionOutputFile.txt");
0805   }
0806   //
0807   unsigned size1 = theCurrentEntry.size() * theCurrentEntry.begin()->size();
0808   std::vector<unsigned> theCurrentEntries;
0809   theCurrentEntries.resize(size1);
0810   size1 *= sizeof(unsigned);
0811   //
0812   unsigned size2 = theCurrentInteraction.size() * theCurrentInteraction.begin()->size();
0813   std::vector<unsigned> theCurrentInteractions;
0814   theCurrentInteractions.resize(size2);
0815   size2 *= sizeof(unsigned);
0816 
0817   // Save the current entries
0818   std::vector<std::vector<unsigned> >::const_iterator aCurrentEntry = theCurrentEntry.begin();
0819   std::vector<std::vector<unsigned> >::const_iterator lastCurrentEntry = theCurrentEntry.end();
0820   unsigned allEntries = 0;
0821   for (; aCurrentEntry != lastCurrentEntry; ++aCurrentEntry) {
0822     unsigned size = aCurrentEntry->size();
0823     for (unsigned iene = 0; iene < size; ++iene)
0824       theCurrentEntries[allEntries++] = (*aCurrentEntry)[iene];
0825   }
0826 
0827   // Save the current interactions
0828   std::vector<std::vector<unsigned> >::const_iterator aCurrentInteraction = theCurrentInteraction.begin();
0829   std::vector<std::vector<unsigned> >::const_iterator lastCurrentInteraction = theCurrentInteraction.end();
0830   unsigned allInteractions = 0;
0831   for (; aCurrentInteraction != lastCurrentInteraction; ++aCurrentInteraction) {
0832     unsigned size = aCurrentInteraction->size();
0833     for (unsigned iene = 0; iene < size; ++iene)
0834       theCurrentInteractions[allInteractions++] = (*aCurrentInteraction)[iene];
0835   }
0836   // Write to file
0837   myOutputFile.write((const char*)(&theCurrentEntries.front()), size1);
0838   myOutputFile.write((const char*)(&theCurrentInteractions.front()), size2);
0839   myOutputFile.flush();
0840 }
0841 
0842 bool fastsim::NuclearInteraction::read(std::string inputFile) {
0843   std::ifstream myInputFile;
0844   struct stat results;
0845   //
0846   unsigned size1 = theCurrentEntry.size() * theCurrentEntry.begin()->size();
0847   std::vector<unsigned> theCurrentEntries;
0848   theCurrentEntries.resize(size1);
0849   size1 *= sizeof(unsigned);
0850   //
0851   unsigned size2 = theCurrentInteraction.size() * theCurrentInteraction.begin()->size();
0852   std::vector<unsigned> theCurrentInteractions;
0853   theCurrentInteractions.resize(size2);
0854   size2 *= sizeof(unsigned);
0855   //
0856   unsigned size = 0;
0857 
0858   // Open the file (if any), otherwise return false
0859   myInputFile.open(inputFile.c_str());
0860   if (myInputFile.is_open()) {
0861     // Get the size of the file
0862     if (stat(inputFile.c_str(), &results) == 0)
0863       size = results.st_size;
0864     else
0865       return false;  // Something is wrong with that file !
0866 
0867     // Position the pointer just before the last record
0868     myInputFile.seekg(size - size1 - size2);
0869     myInputFile.read((char*)(&theCurrentEntries.front()), size1);
0870     myInputFile.read((char*)(&theCurrentInteractions.front()), size2);
0871     myInputFile.close();
0872 
0873     // Read the current entries
0874     std::vector<std::vector<unsigned> >::iterator aCurrentEntry = theCurrentEntry.begin();
0875     std::vector<std::vector<unsigned> >::iterator lastCurrentEntry = theCurrentEntry.end();
0876     unsigned allEntries = 0;
0877     for (; aCurrentEntry != lastCurrentEntry; ++aCurrentEntry) {
0878       unsigned size = aCurrentEntry->size();
0879       for (unsigned iene = 0; iene < size; ++iene)
0880         (*aCurrentEntry)[iene] = theCurrentEntries[allEntries++];
0881     }
0882 
0883     // Read the current interactions
0884     std::vector<std::vector<unsigned> >::iterator aCurrentInteraction = theCurrentInteraction.begin();
0885     std::vector<std::vector<unsigned> >::iterator lastCurrentInteraction = theCurrentInteraction.end();
0886     unsigned allInteractions = 0;
0887     for (; aCurrentInteraction != lastCurrentInteraction; ++aCurrentInteraction) {
0888       unsigned size = aCurrentInteraction->size();
0889       for (unsigned iene = 0; iene < size; ++iene)
0890         (*aCurrentInteraction)[iene] = theCurrentInteractions[allInteractions++];
0891     }
0892 
0893     return true;
0894   }
0895 
0896   return false;
0897 }
0898 
0899 unsigned fastsim::NuclearInteraction::index(int thePid) {
0900   // Find hashed particle ID
0901   unsigned myIndex = 0;
0902   while (thePid != theHadronID[myIndex])
0903     ++myIndex;
0904   return myIndex;
0905 }
0906 
0907 XYZVector fastsim::NuclearInteraction::orthogonal(const XYZVector& aVector) const {
0908   double x = fabs(aVector.X());
0909   double y = fabs(aVector.Y());
0910   double z = fabs(aVector.Z());
0911 
0912   if (x < y)
0913     return x < z ? XYZVector(0., aVector.Z(), -aVector.Y()) : XYZVector(aVector.Y(), -aVector.X(), 0.);
0914   else
0915     return y < z ? XYZVector(-aVector.Z(), 0., aVector.X()) : XYZVector(aVector.Y(), -aVector.X(), 0.);
0916 }
0917 
0918 DEFINE_EDM_PLUGIN(fastsim::InteractionModelFactory, fastsim::NuclearInteraction, "fastsim::NuclearInteraction");