Back to home page

Project CMSSW displayed by LXR
 
 

    

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

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!