MaterialEffects/src/NuclearInteractionSimulator.cc

0001 //Framework Headers
0002 #include "FWCore/Utilities/interface/Exception.h"
0003 #include "FWCore/ParameterSet/interface/FileInPath.h"
0004
0005 #include "FastSimulation/MaterialEffects/interface/NuclearInteractionSimulator.h"
0006 #include "FastSimulation/Particle/interface/makeParticle.h"
0007 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0008
0009 #include "FastSimDataFormats/NuclearInteractions/interface/NUEvent.h"
0010
0011 //#include "DQMServices/Core/interface/DaqMonitorBEInterface.h"
0012 //#include "FWCore/ServiceRegistry/interface/Service.h"
0013
0014 #include <iostream>
0015 #include <sys/stat.h>
0016 #include <cmath>
0017 #include "TFile.h"
0018 #include "TTree.h"
0019 #include "TROOT.h"
0020
0021 // Internal variable and vectors with name started frm "thePion" means
0022 // vectors/variable not only for pions but for all type of hadrons
0023 // treated inside this code
0024
0025 NuclearInteractionSimulator::NuclearInteractionSimulator(std::vector<double>& hadronEnergies,
0026                                                          std::vector<int>& hadronTypes,
0027                                                          std::vector<std::string>& hadronNames,
0028                                                          std::vector<double>& hadronMasses,
0029                                                          std::vector<double>& hadronPMin,
0030                                                          double pionEnergy,
0031                                                          std::vector<double>& lengthRatio,
0032                                                          std::vector<std::vector<double> >& ratios,
0033                                                          std::map<int, int>& idMap,
0034                                                          std::string inputFile,
0035                                                          unsigned int distAlgo,
0036                                                          double distCut)
0037     : MaterialEffectsSimulator(),
0038       thePionEN(hadronEnergies),
0039       thePionID(hadronTypes),
0040       thePionNA(hadronNames),
0041       thePionMA(hadronMasses),
0042       thePionPMin(hadronPMin),
0043       thePionEnergy(pionEnergy),
0044       theLengthRatio(lengthRatio),
0045       theRatios(ratios),
0046       theIDMap(idMap),
0047       theDistAlgo(distAlgo),
0048       theDistCut(distCut),
0049       currentValuesWereSet(false) {
0050   std::string fullPath;
0051
0052   // Prepare the map of files
0053   // Loop over the particle names
0054   TFile* aVFile = nullptr;
0055   std::vector<TTree*> aVTree(thePionEN.size(), static_cast<TTree*>(nullptr));
0056   std::vector<TBranch*> aVBranch(thePionEN.size(), static_cast<TBranch*>(nullptr));
0057   std::vector<NUEvent*> aVNUEvents(thePionEN.size(), static_cast<NUEvent*>(nullptr));
0058   std::vector<unsigned> aVCurrentEntry(thePionEN.size(), static_cast<unsigned>(0));
0059   std::vector<unsigned> aVCurrentInteraction(thePionEN.size(), static_cast<unsigned>(0));
0060   std::vector<unsigned> aVNumberOfEntries(thePionEN.size(), static_cast<unsigned>(0));
0061   std::vector<unsigned> aVNumberOfInteractions(thePionEN.size(), static_cast<unsigned>(0));
0062   std::vector<std::string> aVFileName(thePionEN.size(), static_cast<std::string>(""));
0063   std::vector<double> aVPionCM(thePionEN.size(), static_cast<double>(0));
0064
0065   theTrees.resize(thePionNA.size());
0066   theBranches.resize(thePionNA.size());
0067   theNUEvents.resize(thePionNA.size());
0068   theCurrentEntry.resize(thePionNA.size());
0069   theCurrentInteraction.resize(thePionNA.size());
0070   theNumberOfEntries.resize(thePionNA.size());
0071   theNumberOfInteractions.resize(thePionNA.size());
0072   theFileNames.resize(thePionNA.size());
0073   thePionCM.resize(thePionNA.size());
0074   theFile = aVFile;
0075   for (unsigned iname = 0; iname < thePionNA.size(); ++iname) {
0076     theTrees[iname] = aVTree;
0077     theBranches[iname] = aVBranch;
0078     theNUEvents[iname] = aVNUEvents;
0079     theCurrentEntry[iname] = aVCurrentEntry;
0080     theCurrentInteraction[iname] = aVCurrentInteraction;
0081     theNumberOfEntries[iname] = aVNumberOfEntries;
0082     theNumberOfInteractions[iname] = aVNumberOfInteractions;
0083     theFileNames[iname] = aVFileName;
0084     thePionCM[iname] = aVPionCM;
0085   }
0086
0087   // Read the information from a previous run (to keep reproducibility)
0088   currentValuesWereSet = this->read(inputFile);
0089   if (currentValuesWereSet)
0090     std::cout << "***WARNING*** You are reading nuclear-interaction information from the file " << inputFile
0091               << " created in an earlier run." << std::endl;
0092
0093   // Open the file for saving the information of the current run
0094   myOutputFile.open("NuclearInteractionOutputFile.txt");
0095   myOutputBuffer = 0;
0096
0097   // Open the root files
0098   //  for ( unsigned file=0; file<theFileNames.size(); ++file ) {
0099   edm::FileInPath myDataFile("FastSimulation/MaterialEffects/data/NuclearInteractions.root");
0100
0101   fullPath = myDataFile.fullPath();
0102   theFile = TFile::Open(fullPath.c_str());
0103
0104   for (unsigned iname = 0; iname < thePionNA.size(); ++iname) {
0105     for (unsigned iene = 0; iene < thePionEN.size(); ++iene) {
0106       //std::cout << "iname/iene " << iname << " " << iene << std::endl;
0107       std::ostringstream filename;
0108       double theEne = thePionEN[iene];
0109       filename << "NuclearInteractionsVal_" << thePionNA[iname] << "_E" << theEne << ".root";
0110       theFileNames[iname][iene] = filename.str();
0111       //std::cout << "thePid/theEne " << thePionID[iname] << " " << theEne << std::endl;
0112
0113       std::string treeName = "NuclearInteractions_" + thePionNA[iname] + "_E" + std::to_string(int(theEne));
0114       //
0115       theTrees[iname][iene] = (TTree*)theFile->Get(treeName.c_str());
0116       if (!theTrees[iname][iene])
0117         throw cms::Exception("FastSimulation/MaterialEffects") << "Tree with name " << treeName << " not found ";
0118       //
0119       theBranches[iname][iene] = theTrees[iname][iene]->GetBranch("nuEvent");
0120       //std::cout << "The branch = " << theBranches[iname][iene] << std::endl;
0121       if (!theBranches[iname][iene])
0122         throw cms::Exception("FastSimulation/MaterialEffects")
0123             << "Branch with name nuEvent not found in " << theFileNames[iname][iene];
0124       //
0125       theNUEvents[iname][iene] = new NUEvent();
0126       //std::cout << "The branch = " << theBranches[iname][iene] << std::endl;
0127       theBranches[iname][iene]->SetAddress(&theNUEvents[iname][iene]);
0128       //
0129       theNumberOfEntries[iname][iene] = theTrees[iname][iene]->GetEntries();
0130
0131       if (currentValuesWereSet) {
0132         theTrees[iname][iene]->GetEntry(theCurrentEntry[iname][iene]);
0133         unsigned NInteractions = theNUEvents[iname][iene]->nInteractions();
0134         theNumberOfInteractions[iname][iene] = NInteractions;
0135       }
0136
0137       //
0138       // Compute the corresponding cm energies of the nuclear interactions
0139       XYZTLorentzVector Proton(0., 0., 0., 0.986);
0140       XYZTLorentzVector Reference(
0141           0., 0., std::sqrt(thePionEN[iene] * thePionEN[iene] - thePionMA[iname] * thePionMA[iname]), thePionEN[iene]);
0142       thePionCM[iname][iene] = (Reference + Proton).M();
0143     }
0144   }
0145
0146   // Find the index for which EN = 4. (or thereabout)
0147   ien4 = 0;
0148   while (thePionEN[ien4] < 4.0)
0149     ++ien4;
0150
0151   gROOT->cd();
0152
0153   //  dbe = edm::Service<DaqMonitorBEInterface>().operator->();
0154   //  htot = dbe->book1D("Total", "All particles",150,0.,150.);
0155   //  helas = dbe->book1D("Elastic", "Elastic interactions",150,0.,150.);
0156   //  hinel = dbe->book1D("Inelastic", "Inelastic interactions",150,0.,150.);
0157   //  hscatter = dbe->book1D("Scattering","Elastic Scattering angle",200,0.,2.);
0158   //  hscatter2 = dbe->book2D("Scattering2","Elastic Scattering angle vs p",100,0.,10.,200,0.,2.);
0159   //  hAfter = dbe->book1D("eAfter","Energy after collision",200,0.,4.);
0160   //  hAfter2 = dbe->book2D("eAfter2","Energy after collision",100,-2.5,2.5,100,0.,4);
0161   //  hAfter3 = dbe->book2D("eAfter3","Energy after collision",100,0.,1000.,100,0.,4);
0162 }
0163
0164 NuclearInteractionSimulator::~NuclearInteractionSimulator() {
0165   // Close all local files
0166   // Among other things, this allows the TROOT destructor to end up
0167   // without crashing, while trying to close these files from outside
0168   theFile->Close();
0169   delete theFile;
0170
0171   for (auto& vEvents : theNUEvents) {
0172     for (auto evtPtr : vEvents) {
0173       delete evtPtr;
0174     }
0175   }
0176
0177   // Close the output file
0178   myOutputFile.close();
0179
0180   //  dbe->save("test.root");
0181 }
0182
0183 void NuclearInteractionSimulator::compute(ParticlePropagator& Particle, RandomEngineAndDistribution const* random) {
0184   if (!currentValuesWereSet) {
0185     currentValuesWereSet = true;
0186     for (unsigned iname = 0; iname < thePionNA.size(); ++iname) {
0187       for (unsigned iene = 0; iene < thePionEN.size(); ++iene) {
0188         theCurrentEntry[iname][iene] = (unsigned)(theNumberOfEntries[iname][iene] * random->flatShoot());
0189
0190         theTrees[iname][iene]->GetEntry(theCurrentEntry[iname][iene]);
0191         unsigned NInteractions = theNUEvents[iname][iene]->nInteractions();
0192         theNumberOfInteractions[iname][iene] = NInteractions;
0193
0194         theCurrentInteraction[iname][iene] = (unsigned)(theNumberOfInteractions[iname][iene] * random->flatShoot());
0195       }
0196     }
0197   }
0198
0199   // Read a Nuclear Interaction in a random manner
0200
0201   double pHadron = std::sqrt(Particle.particle().Vect().Mag2());
0202   //  htot->Fill(pHadron);
0203
0204   // The hadron has enough momentum to create some relevant final state
0205   if (pHadron > thePionEnergy) {
0206     // The particle type
0207     std::map<int, int>::const_iterator thePit = theIDMap.find(Particle.particle().pid());
0208
0209     int thePid = thePit != theIDMap.end() ? thePit->second : Particle.particle().pid();
0210
0211     // Is this particle type foreseen?
0212     unsigned fPid = abs(thePid);
0213     if (fPid != 211 && fPid != 130 && fPid != 321 && fPid != 2112 && fPid != 2212) {
0214       return;
0215       //std::cout << "Unknown particle type = " << thePid << std::endl;
0216       //thePid = 211;
0217     }
0218
0219     // The inelastic interaction length at p(pion) = 5 GeV/c
0220     unsigned thePidIndex = index(thePid);
0221     double theInelasticLength = radLengths * theLengthRatio[thePidIndex];
0222
0223     // The elastic interaction length
0224     // The baroque parameterization is a fit to Fig. 40.13 of the PDG
0225     double ee = pHadron > 0.6 ? exp(-std::sqrt((pHadron - 0.6) / 1.122)) : exp(std::sqrt((0.6 - pHadron) / 1.122));
0226     double theElasticLength = (0.8753 * ee + 0.15)
0227                               //    double theElasticLength = ( 0.15 + 0.195 / log(pHadron/0.4) )
0228                               //    double theElasticLength = ( 0.15 + 0.305 / log(pHadron/0.35) )
0229                               * theInelasticLength;
0230
0231     // The total interaction length
0232     double theTotalInteractionLength = theInelasticLength + theElasticLength;
0233
0234     // Probability to interact is dl/L0 (maximum for 4 GeV pion)
0235     double aNuclInteraction = -std::log(random->flatShoot());
0236     if (aNuclInteraction < theTotalInteractionLength) {
0237       // The elastic part
0238       double elastic = random->flatShoot();
0239       if (elastic < theElasticLength / theTotalInteractionLength) {
0240         //          helas->Fill(pHadron);
0241
0242         // Characteristic scattering angle for the elastic part
0243         double theta0 = std::sqrt(3.) / std::pow(theA(), 1. / 3.) * Particle.particle().mass() / pHadron;
0244
0245         // Draw an angle with theta/theta0*exp[(-theta/2theta0)**2] shape
0246         double theta = theta0 * std::sqrt(-2. * std::log(random->flatShoot()));
0247         double phi = 2. * 3.14159265358979323 * random->flatShoot();
0248
0249         // Rotate the particle accordingly
0250         RawParticle::Rotation rotation1(orthogonal(Particle.particle().Vect()), theta);
0251         RawParticle::Rotation rotation2(Particle.particle().Vect(), phi);
0252         Particle.particle().rotate(rotation1);
0253         Particle.particle().rotate(rotation2);
0254
0255         // Distance
0256         double distance = std::sin(theta);
0257
0258         // Create a daughter if the kink is large engough
0259         if (distance > theDistCut) {
0260           _theUpdatedState.reserve(1);
0261           _theUpdatedState.clear();
0262           _theUpdatedState.emplace_back(Particle.particle());
0263         }
0264
0265         //  hscatter->Fill(myTheta);
0266         //  hscatter2->Fill(pHadron,myTheta);
0267
0268       }
0269
0270       // The inelastic part
0271       else {
0272         // Avoid multiple map access
0273         const std::vector<double>& aPionCM = thePionCM[thePidIndex];
0274         const std::vector<double>& aRatios = theRatios[thePidIndex];
0275         // Find the file with the closest c.m energy
0276         // The target nucleon
0277         XYZTLorentzVector Proton(0., 0., 0., 0.939);
0278         // The current particle
0279         const XYZTLorentzVector& Hadron = (const XYZTLorentzVector&)Particle;
0280         // The smallest momentum for inelastic interactions
0281         double pMin = thePionPMin[thePidIndex];
0282         // The correspong smallest four vector
0283         XYZTLorentzVector Hadron0(0., 0., pMin, std::sqrt(pMin * pMin + Particle.particle().M2()));
0284
0285         // The current centre-of-mass energy
0286         double ecm = (Proton + Hadron).M();
0287         //std::cout << "Proton = " << Proton << std::endl;
0288         //std::cout << "Hadron = " << Hadron << std::endl;
0289         //std::cout << "ecm = " << ecm << std::endl;
0290         // Get the files of interest (closest c.m. energies)
0291         unsigned ene1 = 0;
0292         unsigned ene2 = 0;
0293         // The smallest centre-of-mass energy
0294         //  double ecm1=1.63;
0295         double ecm1 = (Proton + Hadron0).M();
0296         //std::cout << "ecm1 = " << ecm1 << std::endl;
0297         //std::cout << "ecm[0] = " << aPionCM[0] << std::endl;
0298         //std::cout << "ecm[11] = " << aPionCM [ aPionCM.size()-1 ] << std::endl;
0299         double ecm2 = aPionCM[0];
0300         double ratio1 = 0.;
0301         double ratio2 = aRatios[0];
0302         if (ecm > aPionCM[0] && ecm < aPionCM[aPionCM.size() - 1]) {
0303           for (unsigned ene = 1; ene < aPionCM.size() && ecm > aPionCM[ene - 1]; ++ene) {
0304             if (ecm < aPionCM[ene]) {
0305               ene2 = ene;
0306               ene1 = ene2 - 1;
0307               ecm1 = aPionCM[ene1];
0308               ecm2 = aPionCM[ene2];
0309               ratio1 = aRatios[ene1];
0310               ratio2 = aRatios[ene2];
0311             }
0312           }
0313         } else if (ecm > aPionCM[aPionCM.size() - 1]) {
0314           ene1 = aPionCM.size() - 1;
0315           ene2 = aPionCM.size() - 2;
0316           ecm1 = aPionCM[ene1];
0317           ecm2 = aPionCM[ene2];
0318           ratio1 = aRatios[ene2];
0319           ratio2 = aRatios[ene2];
0320         }
0321
0322         // The inelastic part of the cross section depends cm energy
0323         double slope = (std::log10(ecm) - std::log10(ecm1)) / (std::log10(ecm2) - std::log10(ecm1));
0324         double inelastic = ratio1 + (ratio2 - ratio1) * slope;
0325         double inelastic4 = pHadron < 4. ? aRatios[ien4] : 1.;
0326
0327         //std::cout << "Inelastic = " << ratio1 << " " << ratio2 << " " << inelastic << std::endl;
0328         //      std::cout << "Energy/Inelastic : "
0329         //      << Hadron.e() << " " << inelastic << std::endl;
0330
0331         // Simulate an inelastic interaction
0332         if (elastic > 1. - (inelastic * theInelasticLength) / theTotalInteractionLength) {
0333           // Avoid mutliple map access
0334           std::vector<unsigned>& aCurrentInteraction = theCurrentInteraction[thePidIndex];
0335           std::vector<unsigned>& aNumberOfInteractions = theNumberOfInteractions[thePidIndex];
0336           std::vector<NUEvent*>& aNUEvents = theNUEvents[thePidIndex];
0337           //      hinel->Fill(pHadron);
0338           //      std::cout << "INELASTIC INTERACTION ! "
0339           //            << pHadron << " " << theInelasticLength << " "
0340           //            << inelastic * theInelasticLength << std::endl;
0341           // Choice of the file to read according the the log10(ecm) distance
0342           // and protection against low momentum proton and neutron that never interacts
0343           // (i.e., empty files)
0344           unsigned ene;
0345           if (random->flatShoot() < slope || aNumberOfInteractions[ene1] == 0)
0346             ene = ene2;
0347           else
0348             ene = ene1;
0349
0350           //std::cout << "Ecm1/2 = " << ecm1 << " " << ecm2 << std::endl;
0351           //std::cout << "Ratio1/2 = " << ratio1 << " " << ratio2 << std::endl;
0352           //std::cout << "Ene = " << ene << " slope = " << slope << std::endl;
0353
0354           //std::cout << "Pion energy = " << Hadron.E()
0355           //        << "File chosen " << theFileNames[thePidIndex][ene]
0356           //        << std::endl;
0357
0358           // The boost characteristics
0359           XYZTLorentzVector theBoost = Proton + Hadron;
0360           theBoost /= theBoost.e();
0361
0362           // std::cout << "File chosen : " << thePid << "/" << ene
0363           //        << " Current interaction = " << aCurrentInteraction[ene]
0364           //        << " Total interactions = " << aNumberOfInteractions[ene]
0365           //        << std::endl;
0366           //      theFiles[thePidIndex][ene]->cd();
0367           //      gDirectory->ls();
0368
0369           // Check we are not either at the end of an interaction bunch
0370           // or at the end of a file
0371           if (aCurrentInteraction[ene] == aNumberOfInteractions[ene]) {
0372             // std::cout << "End of interaction bunch ! ";
0373             std::vector<unsigned>& aCurrentEntry = theCurrentEntry[thePidIndex];
0374             std::vector<unsigned>& aNumberOfEntries = theNumberOfEntries[thePidIndex];
0375             std::vector<TTree*>& aTrees = theTrees[thePidIndex];
0376             ++aCurrentEntry[ene];
0377             // std::cerr << "Read the next entry "
0378             //           << aCurrentEntry[ene] << std::endl;
0379             aCurrentInteraction[ene] = 0;
0380             if (aCurrentEntry[ene] == aNumberOfEntries[ene]) {
0381               aCurrentEntry[ene] = 0;
0382               //  std::cout << "End of file - Rewind! " << std::endl;
0383             }
0384             unsigned myEntry = aCurrentEntry[ene];
0385             // std::cout << "The new entry " << myEntry
0386             //           << " is read ... in TTree " << aTrees[ene] << " ";
0387             aTrees[ene]->GetEntry(myEntry);
0388             // std::cout << "The number of interactions in the new entry is ... ";
0389             aNumberOfInteractions[ene] = aNUEvents[ene]->nInteractions();
0390             // std::cout << aNumberOfInteractions[ene] << std::endl;
0391           }
0392
0393           // Read the interaction
0394           NUEvent::NUInteraction anInteraction = aNUEvents[ene]->theNUInteractions()[aCurrentInteraction[ene]];
0395
0396           unsigned firstTrack = anInteraction.first;
0397           unsigned lastTrack = anInteraction.last;
0398           //      std::cout << "First and last tracks are " << firstTrack << " " << lastTrack << std::endl;
0399
0400           _theUpdatedState.reserve(lastTrack - firstTrack + 1);
0401           _theUpdatedState.clear();
0402
0403           double distMin = 1E99;
0404
0405           // Some rotation around the boost axis, for more randomness
0406           XYZVector theAxis = theBoost.Vect().Unit();
0407           double theAngle = random->flatShoot() * 2. * 3.14159265358979323;
0408           RawParticle::Rotation axisRotation(theAxis, theAngle);
0409           RawParticle::Boost axisBoost(theBoost.x(), theBoost.y(), theBoost.z());
0410
0411           // A rotation to bring the particles back to the pion direction
0412           XYZVector zAxis(0., 0., 1.);
0413           XYZVector orthAxis = (zAxis.Cross(theBoost.Vect())).Unit();
0414           double orthAngle = acos(theBoost.Vect().Unit().Z());
0415           RawParticle::Rotation orthRotation(orthAxis, orthAngle);
0416
0417           // A few checks
0418           // double eAfter = 0.;
0419
0420           // Loop on the nuclear interaction products
0421           for (unsigned iTrack = firstTrack; iTrack <= lastTrack; ++iTrack) {
0422             unsigned idaugh = iTrack - firstTrack;
0423             NUEvent::NUParticle aParticle = aNUEvents[ene]->theNUParticles()[iTrack];
0424             //  std::cout << "Track " << iTrack
0425             //        << " id/px/py/pz/mass "
0426             //        << aParticle.id << " "
0427             //        << aParticle.px << " "
0428             //        << aParticle.py << " "
0429             //        << aParticle.pz << " "
0430             //        << aParticle.mass << " " << endl;
0431
0432             // Add a RawParticle with the proper energy in the c.m frame of
0433             // the nuclear interaction
0434             double energy = std::sqrt(aParticle.px * aParticle.px + aParticle.py * aParticle.py +
0435                                       aParticle.pz * aParticle.pz + aParticle.mass * aParticle.mass / (ecm * ecm));
0436
0437             RawParticle& aDaughter = _theUpdatedState.emplace_back(makeParticle(
0438                 Particle.particleDataTable(),
0439                 aParticle.id,
0440                 XYZTLorentzVector(aParticle.px * ecm, aParticle.py * ecm, aParticle.pz * ecm, energy * ecm)));
0441             // Rotate to the collision axis
0442             aDaughter.rotate(orthRotation);
0443
0444             // Rotate around the boost axis for more randomness
0445             aDaughter.rotate(axisRotation);
0446
0447             // Boost it in the lab frame
0448             aDaughter.boost(axisBoost);
0449
0450             // Store the closest daughter index (for later tracking purposes, so charged particles only)
0451             double distance = distanceToPrimary(Particle.particle(), aDaughter);
0452             // Find the closest daughter, if closer than a given upper limit.
0453             if (distance < distMin && distance < theDistCut) {
0454               distMin = distance;
0455               theClosestChargedDaughterId = idaugh;
0456             }
0457
0458             // eAfter += aDaughter.E();
0459           }
0460
0461           /*
0462       double eBefore = Particle.E();
0463       double rapid = Particle.momentum().Eta();
0464       if ( eBefore > 0. ) {
0465         hAfter->Fill(eAfter/eBefore);
0466         hAfter2->Fill(rapid,eAfter/eBefore);
0467         hAfter3->Fill(eBefore,eAfter/eBefore);
0468       }
0469       */
0470
0471           // ERROR The way this loops through the events breaks
0472           // replay. Which events are retrieved depends on
0473           // which previous events were processed.
0474
0475           // Increment for next time
0476           ++aCurrentInteraction[ene];
0477
0478           // Simulate a stopping hadron (low momentum)
0479         } else if (pHadron < 4. && elastic > 1. - (inelastic4 * theInelasticLength) / theTotalInteractionLength) {
0480           // A fake particle with 0 momentum as a daughter!
0481           _theUpdatedState.reserve(1);
0482           _theUpdatedState.clear();
0483           _theUpdatedState.emplace_back(
0484               makeParticle(Particle.particleDataTable(), 22, XYZTLorentzVector(0., 0., 0., 0.)));
0485         }
0486       }
0487     }
0488   }
0489 }
0490
0491 double NuclearInteractionSimulator::distanceToPrimary(const RawParticle& Particle, const RawParticle& aDaughter) const {
0492   double distance = 2E99;
0493
0494   // Compute the distance only for charged primaries
0495   if (fabs(Particle.charge()) > 1E-12) {
0496     // The secondary must have the same charge
0497     double chargeDiff = fabs(aDaughter.charge() - Particle.charge());
0498     if (fabs(chargeDiff) < 1E-12) {
0499       // Here are two distance definitions * to be tuned *
0500       switch (theDistAlgo) {
0501         case 1:
0502           // sin(theta12)
0503           distance = (aDaughter.Vect().Unit().Cross(Particle.Vect().Unit())).R();
0504           break;
0505
0506         case 2:
0507           // sin(theta12) * p1/p2
0508           distance = (aDaughter.Vect().Cross(Particle.Vect())).R() / aDaughter.Vect().Mag2();
0509           break;
0510
0511         default:
0512           // Should not happen
0513           distance = 2E99;
0514           break;
0515       }
0516     }
0517   }
0518
0519   return distance;
0520 }
0521
0522 void NuclearInteractionSimulator::save() {
0523   // Size of buffer
0524   ++myOutputBuffer;
0525
0526   // Periodically close the current file and open a new one
0527   if (myOutputBuffer / 1000 * 1000 == myOutputBuffer) {
0528     myOutputFile.close();
0529     myOutputFile.open("NuclearInteractionOutputFile.txt");
0530     //    myOutputFile.seekp(0); // No need to rewind in that case
0531   }
0532   //
0533   unsigned size1 = theCurrentEntry.size() * theCurrentEntry.begin()->size();
0534   std::vector<unsigned> theCurrentEntries;
0535   theCurrentEntries.resize(size1);
0536   size1 *= sizeof(unsigned);
0537   //
0538   unsigned size2 = theCurrentInteraction.size() * theCurrentInteraction.begin()->size();
0539   std::vector<unsigned> theCurrentInteractions;
0540   theCurrentInteractions.resize(size2);
0541   size2 *= sizeof(unsigned);
0542
0543   // Save the current entries
0544   std::vector<std::vector<unsigned> >::const_iterator aCurrentEntry = theCurrentEntry.begin();
0545   std::vector<std::vector<unsigned> >::const_iterator lastCurrentEntry = theCurrentEntry.end();
0546   unsigned allEntries = 0;
0547   for (; aCurrentEntry != lastCurrentEntry; ++aCurrentEntry) {
0548     unsigned size = aCurrentEntry->size();
0549     for (unsigned iene = 0; iene < size; ++iene)
0550       theCurrentEntries[allEntries++] = (*aCurrentEntry)[iene];
0551   }
0552
0553   // Save the current interactions
0554   std::vector<std::vector<unsigned> >::const_iterator aCurrentInteraction = theCurrentInteraction.begin();
0555   std::vector<std::vector<unsigned> >::const_iterator lastCurrentInteraction = theCurrentInteraction.end();
0556   unsigned allInteractions = 0;
0557   for (; aCurrentInteraction != lastCurrentInteraction; ++aCurrentInteraction) {
0558     unsigned size = aCurrentInteraction->size();
0559     for (unsigned iene = 0; iene < size; ++iene)
0560       theCurrentInteractions[allInteractions++] = (*aCurrentInteraction)[iene];
0561   }
0562   //
0563   myOutputFile.write((const char*)(&theCurrentEntries.front()), size1);
0564   myOutputFile.write((const char*)(&theCurrentInteractions.front()), size2);
0565   myOutputFile.flush();
0566 }
0567
0568 bool NuclearInteractionSimulator::read(std::string inputFile) {
0569   std::ifstream myInputFile;
0570   struct stat results;
0571   //
0572   unsigned size1 = theCurrentEntry.size() * theCurrentEntry.begin()->size();
0573   std::vector<unsigned> theCurrentEntries;
0574   theCurrentEntries.resize(size1);
0575   size1 *= sizeof(unsigned);
0576   //
0577   unsigned size2 = theCurrentInteraction.size() * theCurrentInteraction.begin()->size();
0578   std::vector<unsigned> theCurrentInteractions;
0579   theCurrentInteractions.resize(size2);
0580   size2 *= sizeof(unsigned);
0581   //
0582   unsigned size = 0;
0583
0584   // Open the file (if any)
0585   myInputFile.open(inputFile.c_str());
0586   if (myInputFile.is_open()) {
0587     // Get the size of the file
0588     if (stat(inputFile.c_str(), &results) == 0)
0589       size = results.st_size;
0590     else
0591       return false;  // Something is wrong with that file !
0592
0593     // Position the pointer just before the last record
0594     myInputFile.seekg(size - size1 - size2);
0595     myInputFile.read((char*)(&theCurrentEntries.front()), size1);
0596     myInputFile.read((char*)(&theCurrentInteractions.front()), size2);
0597     myInputFile.close();
0598
0599     // Read the current entries
0600     std::vector<std::vector<unsigned> >::iterator aCurrentEntry = theCurrentEntry.begin();
0601     std::vector<std::vector<unsigned> >::iterator lastCurrentEntry = theCurrentEntry.end();
0602     unsigned allEntries = 0;
0603     for (; aCurrentEntry != lastCurrentEntry; ++aCurrentEntry) {
0604       unsigned size = aCurrentEntry->size();
0605       for (unsigned iene = 0; iene < size; ++iene)
0606         (*aCurrentEntry)[iene] = theCurrentEntries[allEntries++];
0607     }
0608
0609     // Read the current interactions
0610     std::vector<std::vector<unsigned> >::iterator aCurrentInteraction = theCurrentInteraction.begin();
0611     std::vector<std::vector<unsigned> >::iterator lastCurrentInteraction = theCurrentInteraction.end();
0612     unsigned allInteractions = 0;
0613     for (; aCurrentInteraction != lastCurrentInteraction; ++aCurrentInteraction) {
0614       unsigned size = aCurrentInteraction->size();
0615       for (unsigned iene = 0; iene < size; ++iene)
0616         (*aCurrentInteraction)[iene] = theCurrentInteractions[allInteractions++];
0617     }
0618
0619     return true;
0620   }
0621
0622   return false;
0623 }
0624
0625 unsigned NuclearInteractionSimulator::index(int thePid) {
0626   unsigned myIndex = 0;
0627   while (thePid != thePionID[myIndex])
0628     ++myIndex;
0629   //  std::cout << "pid/index = " << thePid << " " << myIndex << std::endl;
0630   return myIndex;
0631 }