CMS 3D CMS Logo

NuclearInteraction.cc
Go to the documentation of this file.
1 // system headers
2 #include <cmath>
3 #include <memory>
4 #include <vector>
5 #include <map>
6 #include <string>
7 #include <fstream>
8 #include <iostream>
9 #include <sys/stat.h>
10 #include <cmath>
11 #include <TFile.h>
12 #include <TTree.h>
13 #include <TROOT.h>
14 #include <Math/AxisAngle.h>
15 #include "Math/Boost.h"
17 
18 // Framework Headers
23 
25 
26 // Fast Sim headers
34 
35 // Math
37 
38 
40 // Author: Patrick Janot
41 // Date: 25-Jan-2007
42 //
43 // Revision: Class structure modified to match SimplifiedGeometryPropagator
44 // Unified the way, the closest charged daughter is defined
45 // S. Kurz, 29 May 2017
47 
48 
49 // TODO: save() function called in destructer. Does that actually make sense?
50 
53 
54 namespace fastsim
55 {
57 
62  {
63  public:
66 
68  ~NuclearInteraction() override;
69 
71 
77  void interact(fastsim::Particle & particle, const SimplifiedGeometry & layer,std::vector<std::unique_ptr<fastsim::Particle> > & secondaries,const RandomEngineAndDistribution & random) override;
78 
79  private:
81  unsigned index(int thePid);
82 
84  XYZVector orthogonal(const XYZVector& aVector) const;
85 
87  void save();
88 
91 
92  double theDistCut;
93  double theHadronEnergy;
95 
96 
98  // Read/Save nuclear interactions from FullSim
100 
101  TFile* theFile = nullptr;
102  std::vector< std::vector<TTree*> > theTrees;
103  std::vector< std::vector<TBranch*> > theBranches;
104  std::vector< std::vector<NUEvent*> > theNUEvents;
105  std::vector< std::vector<unsigned> > theCurrentEntry;
106  std::vector< std::vector<unsigned> > theCurrentInteraction;
107  std::vector< std::vector<unsigned> > theNumberOfEntries;
108  std::vector< std::vector<unsigned> > theNumberOfInteractions;
109  std::vector< std::vector<std::string> > theFileNames;
110  std::vector< std::vector<double> > theHadronCM;
111 
112  unsigned ien4;
113 
114  std::ofstream myOutputFile;
115  unsigned myOutputBuffer;
116 
118 
120  // Properties of the Hadrons
122 
124  static std::vector< std::vector<double> > theRatiosMap;
126  static std::map<int, int> theIDMap;
127 
129  const std::vector<double> theHadronEN = {1.0, 2.0, 3.0, 4.0, 5.0, 7.0, 9.0, 12.0, 15.0, 20.0, 30.0, 50.0, 100.0, 200.0, 300.0, 500.0, 700.0, 1000.0};
130 
131  // Use same order for all vectors below (Properties for "piplus", "piminus", "K0L", "Kplus", "Kminus", "p", "pbar", "n", "nbar")
132 
134  const std::vector<int> theHadronID = {211, -211, 130, 321, -321, 2212, -2212, 2112, -2112};
136  const std::vector<std::string> theHadronNA = {"piplus", "piminus", "K0L", "Kplus", "Kminus", "p", "pbar", "n", "nbar"};
138  const std::vector<double> theHadronMA = {0.13957, 0.13957, 0.497648, 0.493677, 0.493677, 0.93827, 0.93827, 0.939565, 0.939565};
140  const std::vector<double> theHadronPMin = {0.7, 0.0, 1.0, 1.0, 0.0, 1.1, 0.0, 1.1, 0.0};
142  const std::vector<double> theLengthRatio =
143  {
144  //pi+ pi- K0L K+ K- p pbar n nbar
145  0.2257, 0.2294, 0.3042, 0.2591, 0.2854, 0.3101, 0.5216, 0.3668, 0.4898
146  };
148  const std::vector<double> theRatios =
149  { // pi+ (211)
150  0.031390573,0.531842852,0.819614219,0.951251711,0.986382750,1.000000000,0.985087033,0.982996773,
151  0.990832192,0.992237923,0.994841580,0.973816742,0.967264815,0.971714258,0.969122824,0.978681792,
152  0.977312732,0.984255819,
153  // pi- (-211)
154  0.035326512,0.577356403,0.857118809,0.965683504,0.989659360,1.000000000,0.989599240,0.980665408,
155  0.988384816,0.981038152,0.975002104,0.959996152,0.953310808,0.954705592,0.957615400,0.961150456,
156  0.965022184,0.960573304,
157  // K0L (130)
158  0.000000000,0.370261189,0.649793096,0.734342408,0.749079499,0.753360057,0.755790543,0.755872164,
159  0.751337674,0.746685288,0.747519634,0.739357554,0.735004444,0.803039922,0.832749896,0.890900187,
160  0.936734805,1.000000000,
161  // K+ (321)
162  0.000000000,0.175571717,0.391683394,0.528946472,0.572818635,0.614210280,0.644125538,0.670304050,
163  0.685144573,0.702870161,0.714708513,0.730805263,0.777711536,0.831090576,0.869267129,0.915747562,
164  0.953370523,1.000000000,
165  // K- (-321)
166  0.000000000,0.365353210,0.611663677,0.715315908,0.733498956,0.738361302,0.745253654,0.751459671,
167  0.750628335,0.746442657,0.750850669,0.744895986,0.735093960,0.791663444,0.828609543,0.889993040,
168  0.940897842,1.000000000,
169  // proton (2212)
170  0.000000000,0.042849136,0.459103223,0.666165343,0.787930873,0.890397011,0.920999533,0.937832788,
171  0.950920131,0.966595049,0.979542270,0.988061653,0.983260159,0.988958431,0.991723494,0.995273237,
172  1.000000000,0.999962634,
173  // anti-proton (-2212)
174  1.000000000,0.849956907,0.775625988,0.802018230,0.816207485,0.785899785,0.754998487,0.728977244,
175  0.710010673,0.670890339,0.665627872,0.652682888,0.613334247,0.647534574,0.667910938,0.689919693,
176  0.709200185,0.724199928,
177  // neutron (2112)
178  0.000000000,0.059216484,0.437844536,0.610370629,0.702090648,0.780076890,0.802143073,0.819570432,
179  0.825829666,0.840079750,0.838435509,0.837529986,0.835687165,0.885205014,0.912450156,0.951451221,
180  0.973215562,1.000000000,
181  // anti-neutron
182  1.000000000,0.849573257,0.756479495,0.787147094,0.804572414,0.791806302,0.760234588,0.741109531,
183  0.724118186,0.692829761,0.688465897,0.671806061,0.636461171,0.675314029,0.699134460,0.724305037,
184  0.742556115,0.758504713
185  };
186 
188  // Used to build the ID map 'theIDMap' (i.e., what is to be considered as a proton, etc...)
190 
191  const std::vector<int> protonsID = {2212, 3222, -101, -102, -103, -104};
192  const std::vector<int> antiprotonsID = {-2212, -3222};
193  const std::vector<int> neutronsID = {2112, 3122, 3112, 3312, 3322, 3334, -3334};
194  const std::vector<int> antineutronsID = {-2112, -3122, -3112, -3312, -3322};
195  const std::vector<int> K0LsID = {130, 310};
196  const std::vector<int> KplussesID = {321};
197  const std::vector<int> KminussesID = {-321};
198  const std::vector<int> PiplussesID = {211};
199  const std::vector<int> PiminussesID = {-211};
200  };
201 
202  // TODO: Is this correct?
203  // Thread safety
204  static std::once_flag initializeOnce;
205  CMS_THREAD_GUARD(initializeOnce) std::vector< std::vector<double> > NuclearInteraction::theRatiosMap;
206  CMS_THREAD_GUARD(initializeOnce) std::map<int, int> NuclearInteraction::theIDMap;
207 }
208 
209 
211  : fastsim::InteractionModel(name)
213 {
214  // Full path to FullSim root file
215  std::string fullPath;
216 
217  // Read from config file
218  theDistCut = cfg.getParameter<double>("distCut");
219  theHadronEnergy = cfg.getParameter<double>("hadronEnergy");
220  inputFile = cfg.getUntrackedParameter<std::string>("inputFile","");
221 
222  // The evolution of the interaction lengths with energy
223  // initialize once for all possible instances
224  std::call_once(initializeOnce, [this] () {
225  theRatiosMap.resize(theHadronID.size());
226  for(unsigned i=0; i<theHadronID.size(); ++i){
227  for(unsigned j=0; j<theHadronEN.size(); ++j){
228  theRatiosMap[i].push_back(theRatios[i*theHadronEN.size() + j]);
229  }
230  }
231 
232  // Build the ID map (i.e., what is to be considered as a proton, etc...)
233  // Protons
234  for(const auto & id : protonsID) theIDMap[id] = 2212;
235  // Anti-Protons
236  for(const auto & id : antiprotonsID) theIDMap[id] = -2212;
237  // Neutrons
238  for(const auto & id : neutronsID) theIDMap[id] = 2112;
239  // Anti-Neutrons
240  for(const auto & id : antineutronsID) theIDMap[id] = -2112;
241  // K0L's
242  for(const auto & id : K0LsID) theIDMap[id] = 130;
243  // K+'s
244  for(const auto & id : KplussesID) theIDMap[id] = 321;
245  // K-'s
246  for(const auto & id : KminussesID) theIDMap[id] = -321;
247  // pi+'s
248  for(const auto & id : PiplussesID) theIDMap[id] = 211;
249  // pi-'s
250  for(const auto & id : PiminussesID) theIDMap[id] = -211;
251  });
252 
253  // Prepare the map of files
254  // Loop over the particle names
255  TFile* aVFile=nullptr;
256  std::vector<TTree*> aVTree(theHadronEN.size());
257  std::vector<TBranch*> aVBranch(theHadronEN.size());
258  std::vector<NUEvent*> aVNUEvents(theHadronEN.size());
259  std::vector<unsigned> aVCurrentEntry(theHadronEN.size());
260  std::vector<unsigned> aVCurrentInteraction(theHadronEN.size());
261  std::vector<unsigned> aVNumberOfEntries(theHadronEN.size());
262  std::vector<unsigned> aVNumberOfInteractions(theHadronEN.size());
263  std::vector<std::string> aVFileName(theHadronEN.size());
264  std::vector<double> aVHadronCM(theHadronEN.size());
265  theTrees.resize(theHadronNA.size());
266  theBranches.resize(theHadronNA.size());
267  theNUEvents.resize(theHadronNA.size());
268  theCurrentEntry.resize(theHadronNA.size());
269  theCurrentInteraction.resize(theHadronNA.size());
270  theNumberOfEntries.resize(theHadronNA.size());
271  theNumberOfInteractions.resize(theHadronNA.size());
272  theFileNames.resize(theHadronNA.size());
273  theHadronCM.resize(theHadronNA.size());
274  theFile = aVFile;
275  for(unsigned iname=0; iname<theHadronNA.size(); ++iname){
276  theTrees[iname] = aVTree;
277  theBranches[iname] = aVBranch;
278  theNUEvents[iname] = aVNUEvents;
279  theCurrentEntry[iname] = aVCurrentEntry;
280  theCurrentInteraction[iname] = aVCurrentInteraction;
281  theNumberOfEntries[iname] = aVNumberOfEntries;
282  theNumberOfInteractions[iname] = aVNumberOfInteractions;
283  theFileNames[iname] = aVFileName;
284  theHadronCM[iname] = aVHadronCM;
285  }
286 
287  // Read the information from a previous run (to keep reproducibility)
288  currentValuesWereSet = this->read(inputFile);
290  std::cout << "***WARNING*** You are reading nuclear-interaction information from the file "
291  << inputFile << " created in an earlier run."
292  << std::endl;
293 
294  // Open the file for saving the information of the current run
295  myOutputFile.open("NuclearInteractionOutputFile.txt");
296  myOutputBuffer = 0;
297 
298  // Open the root file
299  edm::FileInPath myDataFile("FastSimulation/MaterialEffects/data/NuclearInteractions.root");
300  fullPath = myDataFile.fullPath();
301  theFile = TFile::Open(fullPath.c_str());
302 
303  // Open the trees
304  unsigned fileNb = 0;
305  for(unsigned iname=0; iname<theHadronNA.size(); ++iname){
306  for(unsigned iene=0; iene<theHadronEN.size(); ++iene){
307  std::ostringstream filename;
308  double theEne = theHadronEN[iene];
309  filename << "NuclearInteractionsVal_" << theHadronNA[iname] << "_E"<< theEne << ".root";
310  theFileNames[iname][iene] = filename.str();
311 
312  ++fileNb;
313  std::string treeName="NuclearInteractions_"+theHadronNA[iname]+"_E"+std::to_string(int(theEne));
314  theTrees[iname][iene] = (TTree*) theFile->Get(treeName.c_str());
315 
316  if ( !theTrees[iname][iene] ) throw cms::Exception("FastSimulation/MaterialEffects")
317  << "Tree with name " << treeName << " not found ";
318  theBranches[iname][iene] = theTrees[iname][iene]->GetBranch("nuEvent");
319  if ( !theBranches[iname][iene] ) throw cms::Exception("FastSimulation/MaterialEffects")
320  << "Branch with name nuEvent not found in " << theFileNames[iname][iene];
321 
322  theNUEvents[iname][iene] = new NUEvent();
323  theBranches[iname][iene]->SetAddress(&theNUEvents[iname][iene]);
324  theNumberOfEntries[iname][iene] = theTrees[iname][iene]->GetEntries();
325 
327  theTrees[iname][iene]->GetEntry(theCurrentEntry[iname][iene]);
328  unsigned NInteractions = theNUEvents[iname][iene]->nInteractions();
329  theNumberOfInteractions[iname][iene] = NInteractions;
330  }
331 
332  // Compute the corresponding cm energies of the nuclear interactions
333  XYZTLorentzVector Proton(0.,0.,0.,0.986);
335  0.,
336  std::sqrt(theHadronEN[iene]*theHadronEN[iene]-theHadronMA[iname]*theHadronMA[iname]),
337  theHadronEN[iene]);
338  theHadronCM[iname][iene] = (Reference+Proton).M();
339  }
340  }
341 
342  // Find the index for which EN = 4. (or thereabout)
343  ien4 = 0;
344  while(theHadronEN[ien4] < 4.0) ++ien4;
345 
346  gROOT->cd();
347 }
348 
350  // Close all local files
351  // Among other things, this allows the TROOT destructor to end up
352  // without crashing, while trying to close these files from outside
353  theFile->Close();
354  delete theFile;
355 
356  for(auto& vEvents: theNUEvents){
357  for(auto evtPtr: vEvents){
358  delete evtPtr;
359  }
360  }
361 
362  // Save data
363  save();
364  // Close the output file
365  myOutputFile.close();
366 }
367 
368 
369 void fastsim::NuclearInteraction::interact(fastsim::Particle & particle, const SimplifiedGeometry & layer,std::vector<std::unique_ptr<fastsim::Particle> > & secondaries,const RandomEngineAndDistribution & random)
370 {
371  int pdgId = particle.pdgId();
372  //
373  // no valid PDGid
374  //
375  if(abs(pdgId) <= 100 || abs(pdgId) >= 1000000)
376  {
377  return;
378  }
379 
380  double radLengths = layer.getThickness(particle.position(),particle.momentum());
381  // TEC layers have some fudge factor for nuclear interactions
382  radLengths *= layer.getNuclearInteractionThicknessFactor();
383  //
384  // no material
385  //
386  if(radLengths < 1E-10)
387  {
388  return;
389  }
390 
391  // In case the events are not read from (old) saved file, then pick a random event from FullSim file
392  if(!currentValuesWereSet) {
393  currentValuesWereSet = true;
394  for ( unsigned iname=0; iname<theHadronNA.size(); ++iname ) {
395  for ( unsigned iene=0; iene<theHadronEN.size(); ++iene ) {
396  theCurrentEntry[iname][iene] = (unsigned) (theNumberOfEntries[iname][iene] * random.flatShoot());
397 
398  theTrees[iname][iene]->GetEntry(theCurrentEntry[iname][iene]);
399  unsigned NInteractions = theNUEvents[iname][iene]->nInteractions();
400  theNumberOfInteractions[iname][iene] = NInteractions;
401 
402  theCurrentInteraction[iname][iene] = (unsigned) (theNumberOfInteractions[iname][iene] * random.flatShoot());
403  }
404  }
405  }
406 
407  // Momentum of interacting hadron
408  double pHadron = std::sqrt(particle.momentum().Vect().Mag2());
409 
410  //
411  // The hadron has not enough momentum to create some relevant final state
412  //
413  if(pHadron < theHadronEnergy){
414  return;
415  }
416 
417  // The particle type
418  std::map<int,int>::const_iterator thePit = theIDMap.find(pdgId);
419  // Use map for unique ID (e.g. proton = {2212, 3222, -101, -102, -103, -104})
420  int thePid = (thePit != theIDMap.end() ? thePit->second : pdgId);
421 
422  // Is this particle type foreseen?
423  unsigned fPid = abs(thePid);
424  if(fPid != 211 && fPid != 130 && fPid != 321 && fPid != 2112 && fPid != 2212)
425  {
426  return;
427  }
428 
429  // The hashed ID
430  unsigned thePidIndex = index(thePid);
431  // The inelastic interaction length at p(Hadron) = 5 GeV/c
432  double theInelasticLength = radLengths * theLengthRatio[thePidIndex];
433 
434  // The elastic interaction length
435  // The baroque parameterization is a fit to Fig. 40.13 of the PDG
436  double ee = pHadron > 0.6 ? exp(-std::sqrt((pHadron-0.6)/1.122)) : exp(std::sqrt((0.6-pHadron)/1.122));
437  double theElasticLength = (0.8753 * ee + 0.15) * theInelasticLength;
438 
439  // The total interaction length
440  double theTotalInteractionLength = theInelasticLength + theElasticLength;
441 
442  //
443  // Probability to interact is dl/L0 (maximum for 4 GeV Hadron)
444  //
445  double aNuclInteraction = -std::log(random.flatShoot());
446  if(aNuclInteraction > theTotalInteractionLength)
447  {
448  return;
449  }
450 
451  // The elastic part
452  double elastic = random.flatShoot();
453  if(elastic < theElasticLength/theTotalInteractionLength){
454  // Characteristic scattering angle for the elastic part
455  // A of silicon
456  double A = 28.0855;
457  double theta0 = std::sqrt(3.)/ std::pow(A,1./3.) * particle.momentum().mass()/pHadron;
458 
459  // Draw an angle with theta/theta0*exp[(-theta/2theta0)**2] shape
460  double theta = theta0 * std::sqrt(-2.*std::log(random.flatShoot()));
461  double phi = 2. * M_PI * random.flatShoot();
462 
463  // Rotate the particle accordingly
464  ROOT::Math::AxisAngle rotation1(orthogonal(particle.momentum().Vect()),theta);
465  ROOT::Math::AxisAngle rotation2(particle.momentum().Vect(),phi);
466  // Rotate!
467  XYZVector rotated = rotation2((rotation1(particle.momentum().Vect())));
468 
469  // Create a daughter if the kink is large enough
470  if (std::sin(theta) > theDistCut) {
471  secondaries.emplace_back(new fastsim::Particle(pdgId
472  ,particle.position()
473  ,XYZTLorentzVector(rotated.X(),
474  rotated.Y(),
475  rotated.Z(),
476  particle.momentum().E())));
477 
478  // Set the ClosestCharged Daughter thing for tracking
479  if(particle.charge() != 0){
480  secondaries.back()->setMotherDeltaR(particle.momentum());
481  secondaries.back()->setMotherPdgId(pdgId);
482  secondaries.back()->setMotherSimTrackIndex(particle.simTrackIndex());
483  }
484 
485  // The particle is destroyed
486  particle.momentum().SetXYZT(0., 0., 0., 0.);
487  }else{
488  // If kink is small enough just rotate particle
489  particle.momentum().SetXYZT(rotated.X(),
490  rotated.Y(),
491  rotated.Z(),
492  particle.momentum().E());
493  }
494  // The inelastic part
495  }else{
496  // Avoid multiple map access
497  const std::vector<double>& aHadronCM = theHadronCM[thePidIndex];
498  const std::vector<double>& aRatios = theRatiosMap[thePidIndex];
499  // Find the file with the closest c.m energy
500  // The target nucleon
501  XYZTLorentzVector Proton(0.,0.,0.,0.939);
502  // The current particle
503  const XYZTLorentzVector& Hadron = (const XYZTLorentzVector&)particle.momentum();
504  // The smallest momentum for inelastic interactions
505  double pMin = theHadronPMin[thePidIndex];
506  // The corresponding smallest four vector
507  XYZTLorentzVector Hadron0(0.,
508  0.,
509  pMin,
510  std::sqrt(pMin*pMin+particle.momentum().M2()));
511 
512  // The current centre-of-mass energy
513  double ecm = (Proton+Hadron).M();
514 
515  // Get the files of interest (closest c.m. energies)
516  unsigned ene1=0;
517  unsigned ene2=0;
518  // The smallest centre-of-mass energy
519  double ecm1= (Proton+Hadron0).M();
520 
521  double ecm2=aHadronCM[0];
522  double ratio1=0.;
523  double ratio2=aRatios[0];
524  if(ecm > aHadronCM[0] && ecm < aHadronCM[aHadronCM.size()-1]){
525  for(unsigned ene=1; ene < aHadronCM.size() && ecm > aHadronCM[ene-1]; ++ene){
526  if(ecm < aHadronCM[ene]){
527  ene2 = ene;
528  ene1 = ene2-1;
529  ecm1 = aHadronCM[ene1];
530  ecm2 = aHadronCM[ene2];
531  ratio1 = aRatios[ene1];
532  ratio2 = aRatios[ene2];
533  }
534  }
535  }else if(ecm > aHadronCM[aHadronCM.size()-1]){
536  ene1 = aHadronCM.size()-1;
537  ene2 = aHadronCM.size()-2;
538  ecm1 = aHadronCM[ene1];
539  ecm2 = aHadronCM[ene2];
540  ratio1 = aRatios[ene2];
541  ratio2 = aRatios[ene2];
542  }
543 
544  // The inelastic part of the cross section depends cm energy
545  double slope = (std::log10(ecm) - std::log10(ecm1)) / (std::log10(ecm2) - std::log10(ecm1));
546  double inelastic = ratio1 + (ratio2 - ratio1) * slope;
547  double inelastic4 = pHadron < 4. ? aRatios[ien4] : 1.;
548 
549  // Simulate an inelastic interaction
550  if(elastic > 1.- (inelastic*theInelasticLength) / theTotalInteractionLength){
551  // Avoid mutliple map access
552  std::vector<unsigned>& aCurrentInteraction = theCurrentInteraction[thePidIndex];
553  std::vector<unsigned>& aNumberOfInteractions = theNumberOfInteractions[thePidIndex];
554  std::vector<NUEvent*>& aNUEvents = theNUEvents[thePidIndex];
555 
556  // Choice of the file to read according the the log10(ecm) distance
557  // and protection against low momentum proton and neutron that never interacts
558  // (i.e., empty files)
559  unsigned ene;
560  if(random.flatShoot() < slope || aNumberOfInteractions[ene1] == 0){
561  ene = ene2;
562  }else{
563  ene = ene1;
564  }
565 
566  // The boost characteristics
567  XYZTLorentzVector theBoost = Proton + Hadron;
568  theBoost /= theBoost.e();
569 
570  // Check we are not either at the end of an interaction bunch
571  // or at the end of a file
572  if(aCurrentInteraction[ene] == aNumberOfInteractions[ene]){
573  std::vector<unsigned>& aCurrentEntry = theCurrentEntry[thePidIndex];
574  std::vector<unsigned>& aNumberOfEntries = theNumberOfEntries[thePidIndex];
575  std::vector<TTree*>& aTrees = theTrees[thePidIndex];
576  ++aCurrentEntry[ene];
577 
578  aCurrentInteraction[ene] = 0;
579  if(aCurrentEntry[ene] == aNumberOfEntries[ene]){
580  aCurrentEntry[ene] = 0;
581  }
582 
583  unsigned myEntry = aCurrentEntry[ene];
584  aTrees[ene]->GetEntry(myEntry);
585  aNumberOfInteractions[ene] = aNUEvents[ene]->nInteractions();
586  }
587 
588  // Read the interaction
589  NUEvent::NUInteraction anInteraction = aNUEvents[ene]->theNUInteractions()[aCurrentInteraction[ene]];
590 
591  unsigned firstTrack = anInteraction.first;
592  unsigned lastTrack = anInteraction.last;
593 
594  // Some rotation around the boost axis, for more randomness
595  XYZVector theAxis = theBoost.Vect().Unit();
596  double theAngle = random.flatShoot() * 2. * M_PI;
597  ROOT::Math::AxisAngle axisRotation(theAxis,theAngle);
598  ROOT::Math::Boost axisBoost(theBoost.x(),theBoost.y(),theBoost.z());
599 
600  // A rotation to bring the particles back to the Hadron direction
601  XYZVector zAxis(0.,0.,1.);
602  XYZVector orthAxis = (zAxis.Cross(theBoost.Vect())).Unit();
603  double orthAngle = acos(theBoost.Vect().Unit().Z());
604  ROOT::Math::AxisAngle orthRotation(orthAxis,orthAngle);
605 
606  // Loop on the nuclear interaction products
607  for(unsigned iTrack=firstTrack; iTrack<=lastTrack; ++iTrack){
608  NUEvent::NUParticle aParticle = aNUEvents[ene]->theNUParticles()[iTrack];
609 
610  // Add a RawParticle with the proper energy in the c.m. frame of
611  // the nuclear interaction
612  double energy = std::sqrt(aParticle.px*aParticle.px
613  + aParticle.py*aParticle.py
614  + aParticle.pz*aParticle.pz
615  + aParticle.mass*aParticle.mass / (ecm * ecm));
616 
617  XYZTLorentzVector daugtherMomentum(aParticle.px*ecm, aParticle.py*ecm, aParticle.pz*ecm, energy*ecm);
618 
619  // Rotate to the collision axis
620  XYZVector rotated = orthRotation(daugtherMomentum.Vect());
621  // Rotate around the boost axis for more randomness
622  rotated = axisRotation(rotated);
623 
624  // Rotated the daughter
625  daugtherMomentum.SetXYZT(rotated.X(), rotated.Y(), rotated.Z(), daugtherMomentum.E());
626 
627  // Boost it in the lab frame
628  daugtherMomentum = axisBoost(daugtherMomentum);
629 
630  // Create secondary
631  secondaries.emplace_back(new fastsim::Particle(aParticle.id, particle.position(), daugtherMomentum));
632 
633  // The closestCharged Daughter thing for tracking
634  // BUT: 'aParticle' also has to be charged, only then the mother should be set
635  // Unfortunately, NUEvent::NUParticle does not contain any info about the charge
636  // Did some tests and effect is absolutely negligible!
637  if(particle.charge() != 0){
638  secondaries.back()->setMotherDeltaR(particle.momentum());
639  secondaries.back()->setMotherPdgId(pdgId);
640  secondaries.back()->setMotherSimTrackIndex(particle.simTrackIndex());
641  }
642  }
643 
644  // The particle is destroyed
645  particle.momentum().SetXYZT(0., 0., 0., 0.);
646 
647  // This is a note from previous version of code but I don't understand it:
648  // ERROR The way this loops through the events breaks
649  // replay. Which events are retrieved depends on
650  // which previous events were processed.
651 
652  // Increment for next time
653  ++aCurrentInteraction[ene];
654 
655  // Simulate a stopping hadron (low momentum)
656  }else if(pHadron < 4. && elastic > 1.- (inelastic4*theInelasticLength) / theTotalInteractionLength){
657  // The particle is destroyed
658  particle.momentum().SetXYZT(0., 0., 0., 0.);
659  }
660  }
661 }
662 
663 void
665 {
666  // Size of buffer
667  ++myOutputBuffer;
668 
669  // Periodically close the current file and open a new one
670  if(myOutputBuffer/1000*1000 == myOutputBuffer){
671  myOutputFile.close();
672  myOutputFile.open ("NuclearInteractionOutputFile.txt");
673  }
674  //
675  unsigned size1 = theCurrentEntry.size() * theCurrentEntry.begin()->size();
676  std::vector<unsigned> theCurrentEntries;
677  theCurrentEntries.resize(size1);
678  size1*=sizeof(unsigned);
679  //
680  unsigned size2 = theCurrentInteraction.size() * theCurrentInteraction.begin()->size();
681  std::vector<unsigned> theCurrentInteractions;
682  theCurrentInteractions.resize(size2);
683  size2 *= sizeof(unsigned);
684 
685  // Save the current entries
686  std::vector< std::vector<unsigned> >::const_iterator aCurrentEntry = theCurrentEntry.begin();
687  std::vector< std::vector<unsigned> >::const_iterator lastCurrentEntry = theCurrentEntry.end();
688  unsigned allEntries=0;
689  for(; aCurrentEntry!=lastCurrentEntry; ++aCurrentEntry){
690  unsigned size = aCurrentEntry->size();
691  for(unsigned iene=0; iene<size; ++iene)
692  theCurrentEntries[allEntries++] = (*aCurrentEntry)[iene];
693  }
694 
695  // Save the current interactions
696  std::vector< std::vector<unsigned> >::const_iterator aCurrentInteraction = theCurrentInteraction.begin();
697  std::vector< std::vector<unsigned> >::const_iterator lastCurrentInteraction = theCurrentInteraction.end();
698  unsigned allInteractions=0;
699  for (; aCurrentInteraction!=lastCurrentInteraction; ++aCurrentInteraction){
700  unsigned size = aCurrentInteraction->size();
701  for(unsigned iene=0; iene<size; ++iene)
702  theCurrentInteractions[allInteractions++] = (*aCurrentInteraction)[iene];
703  }
704  // Write to file
705  myOutputFile.write((const char*)(&theCurrentEntries.front()), size1);
706  myOutputFile.write((const char*)(&theCurrentInteractions.front()), size2);
707  myOutputFile.flush();
708 }
709 
710 bool
712 {
713  std::ifstream myInputFile;
714  struct stat results;
715  //
716  unsigned size1 = theCurrentEntry.size() * theCurrentEntry.begin()->size();
717  std::vector<unsigned> theCurrentEntries;
718  theCurrentEntries.resize(size1);
719  size1*=sizeof(unsigned);
720  //
721  unsigned size2 = theCurrentInteraction.size() * theCurrentInteraction.begin()->size();
722  std::vector<unsigned> theCurrentInteractions;
723  theCurrentInteractions.resize(size2);
724  size2 *= sizeof(unsigned);
725  //
726  unsigned size = 0;
727 
728 
729  // Open the file (if any), otherwise return false
730  myInputFile.open (inputFile.c_str());
731  if(myInputFile.is_open()){
732 
733  // Get the size of the file
734  if(stat(inputFile.c_str(), &results) == 0) size = results.st_size;
735  else return false; // Something is wrong with that file !
736 
737  // Position the pointer just before the last record
738  myInputFile.seekg(size-size1-size2);
739  myInputFile.read((char*)(&theCurrentEntries.front()),size1);
740  myInputFile.read((char*)(&theCurrentInteractions.front()),size2);
741  myInputFile.close();
742 
743  // Read the current entries
744  std::vector< std::vector<unsigned> >::iterator aCurrentEntry = theCurrentEntry.begin();
745  std::vector< std::vector<unsigned> >::iterator lastCurrentEntry = theCurrentEntry.end();
746  unsigned allEntries=0;
747  for(; aCurrentEntry!=lastCurrentEntry; ++aCurrentEntry){
748  unsigned size = aCurrentEntry->size();
749  for(unsigned iene=0; iene<size; ++iene)
750  (*aCurrentEntry)[iene] = theCurrentEntries[allEntries++];
751  }
752 
753  // Read the current interactions
754  std::vector< std::vector<unsigned> >::iterator aCurrentInteraction = theCurrentInteraction.begin();
755  std::vector< std::vector<unsigned> >::iterator lastCurrentInteraction = theCurrentInteraction.end();
756  unsigned allInteractions=0;
757  for(; aCurrentInteraction!=lastCurrentInteraction; ++aCurrentInteraction){
758  unsigned size = aCurrentInteraction->size();
759  for(unsigned iene=0; iene<size; ++iene)
760  (*aCurrentInteraction)[iene] = theCurrentInteractions[allInteractions++];
761  }
762 
763  return true;
764  }
765 
766  return false;
767 }
768 
769 unsigned
771 {
772  // Find hashed particle ID
773  unsigned myIndex=0;
774  while(thePid != theHadronID[myIndex]) ++myIndex;
775  return myIndex;
776 }
777 
778 XYZVector
780 {
781  double x = fabs(aVector.X());
782  double y = fabs(aVector.Y());
783  double z = fabs(aVector.Z());
784 
785  if (x < y)
786  return x < z ?
787  XYZVector(0.,aVector.Z(),-aVector.Y()) :
788  XYZVector(aVector.Y(),-aVector.X(),0.);
789  else
790  return y < z ?
791  XYZVector(-aVector.Z(),0.,aVector.X()) :
792  XYZVector(aVector.Y(),-aVector.X(),0.);
793 }
794 
795 
799  "fastsim::NuclearInteraction"
800  );
const std::vector< int > antineutronsID
PdgID of anti-neutrons.
size
Write out results.
const std::vector< int > PiplussesID
PdgID of pt+.
Implementation of a generic detector layer (base class for forward/barrel layers).
const math::XYZTLorentzVector & position() const
Return position of the particle.
Definition: Particle.h:142
double theDistCut
Cut on deltaR for the FastSim Tracking (ClosestChargedDaughter algorithm)
const std::vector< int > KminussesID
PdgID of K-.
double flatShoot(double xmin=0.0, double xmax=1.0) const
std::vector< std::vector< TTree * > > theTrees
Necessary to read the FullSim interactions.
static const double slope[3]
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
const double getNuclearInteractionThicknessFactor() const
Some layers have a different thickness for nuclear interactions.
Geom::Theta< T > theta() const
bool currentValuesWereSet
Read data from file that was created in a previous run.
std::vector< std::vector< unsigned > > theNumberOfEntries
Necessary to read the FullSim interactions.
const std::vector< int > KplussesID
PdgID of K+.
TFile * theFile
Necessary to read the FullSim interactions.
int pdgId() const
Return pdgId of the particle.
Definition: Particle.h:136
TRandom random
Definition: MVATrainer.cc:138
static std::vector< std::vector< double > > theRatiosMap
The evolution of the interaction lengths with energy.
Implementation of nuclear interactions of hadrons in the tracker layers (based on fully simulated int...
std::vector< std::vector< NUEvent * > > theNUEvents
Necessary to read the FullSim interactions.
const std::vector< int > protonsID
PdgID of protons.
const std::vector< double > theLengthRatio
Inelastic interaction length at p(Hadron) = 5 GeV/c (relativ to radionLength of material) ...
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
Base class for any interaction model between a particle and a tracker layer.
const std::vector< double > theRatios
Filled into &#39;theRatiosMap&#39; (evolution of the interaction lengths with energy)
void interact(fastsim::Particle &particle, const SimplifiedGeometry &layer, std::vector< std::unique_ptr< fastsim::Particle > > &secondaries, const RandomEngineAndDistribution &random) override
Perform the interaction.
const std::vector< int > K0LsID
PdgID of K0.
Definition: NUEvent.h:6
std::vector< std::vector< std::string > > theFileNames
Necessary to read the FullSim interactions.
const std::vector< int > antiprotonsID
PdgID of anti-protons.
static std::map< int, int > theIDMap
Build the ID map (i.e., what is to be considered as a proton, etc...)
std::vector< std::vector< unsigned > > theNumberOfInteractions
Necessary to read the FullSim interactions.
T sqrt(T t)
Definition: SSEVec.h:18
std::vector< std::vector< unsigned > > theCurrentInteraction
Necessary to read the FullSim interactions.
~NuclearInteraction() override
Default destructor.
unsigned index(int thePid)
Return a hashed index for a given particle ID.
#define CMS_THREAD_GUARD(_var_)
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
XYZVector orthogonal(const XYZVector &aVector) const
Return an orthogonal vector.
std::ofstream myOutputFile
Output file to save interactions.
math::XYZTLorentzVector XYZTLorentzVector
int simTrackIndex() const
Return index of the SimTrack.
Definition: Particle.h:155
const std::vector< int > neutronsID
PdgID of neutrons.
#define M_PI
void save()
Save the nuclear interactions to a file, so you can reproduce the events (e.g. in case of a crash)...
std::vector< std::vector< double > > theHadronCM
Necessary to read the FullSim interactions.
double theHadronEnergy
Minimum energy for nuclear interaction.
math::XYZVector XYZVector
std::vector< std::vector< unsigned > > theCurrentEntry
Necessary to read the FullSim interactions.
NuclearInteraction(const std::string &name, const edm::ParameterSet &cfg)
Constructor.
std::vector< std::vector< TBranch * > > theBranches
Necessary to read the FullSim interactions.
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:30
bool read(std::string inputFile)
Read the nuclear interactions from a file, so you can reproduce the events (e.g. in case of a crash)...
const std::vector< double > theHadronPMin
Smallest momentum for inelastic interactions.
double charge() const
Return charge of the particle.
Definition: Particle.h:139
virtual const double getThickness(const math::XYZTLorentzVector &position) const =0
Return thickness of the layer at a given position.
unsigned ien4
Find the index for which EN = 4.
const std::vector< double > theHadronMA
Masses of the hadrons.
const std::vector< int > PiminussesID
PdgID of pi-.
HLT enums.
const std::vector< double > theHadronEN
Filled into &#39;theRatiosMap&#39; (evolution of the interaction lengths with energy)
static std::once_flag initializeOnce
const std::vector< int > theHadronID
ID of the hadrons.
std::string fullPath() const
Definition: FileInPath.cc:163
const math::XYZTLorentzVector & momentum() const
Return momentum of the particle.
Definition: Particle.h:145
#define DEFINE_EDM_PLUGIN(factory, type, name)
std::string inputFile
Directory/Name of input file in case you want to read interactions from file.
Definition of a generic FastSim Particle which can be propagated through the detector (formerly Parti...
Definition: Particle.h:19
math::XYZVector XYZVector
Definition: RawParticle.h:28
const std::vector< std::string > theHadronNA
Names of the hadrons.
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
math::XYZTLorentzVector XYZTLorentzVector
Definition: RawParticle.h:27
unsigned myOutputBuffer
Needed to save interactions to file.