File indexing completed on 2024-04-06 12:11:21
0001
0002 #include <cmath>
0003 #include <memory>
0004 #include <vector>
0005 #include <map>
0006 #include <string>
0007 #include <fstream>
0008 #include <iostream>
0009 #include <sys/stat.h>
0010 #include <cmath>
0011 #include <TFile.h>
0012 #include <TTree.h>
0013 #include <TROOT.h>
0014 #include <Math/AxisAngle.h>
0015 #include "Math/Boost.h"
0016 #include "DataFormats/Math/interface/Vector3D.h"
0017
0018
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/Utilities/interface/Exception.h"
0022 #include "FWCore/ParameterSet/interface/FileInPath.h"
0023
0024 #include "FWCore/Utilities/interface/thread_safety_macros.h"
0025
0026
0027 #include "FastSimulation/SimplifiedGeometryPropagator/interface/Particle.h"
0028 #include "FastSimulation/SimplifiedGeometryPropagator/interface/SimplifiedGeometry.h"
0029 #include "FastSimulation/SimplifiedGeometryPropagator/interface/InteractionModelFactory.h"
0030 #include "FastSimulation/SimplifiedGeometryPropagator/interface/InteractionModel.h"
0031 #include "FastSimulation/SimplifiedGeometryPropagator/interface/Constants.h"
0032 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0033 #include "FastSimDataFormats/NuclearInteractions/interface/NUEvent.h"
0034
0035
0036 #include "DataFormats/Math/interface/LorentzVector.h"
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049 typedef math::XYZVector XYZVector;
0050 typedef math::XYZTLorentzVector XYZTLorentzVector;
0051
0052 namespace fastsim {
0053
0054
0055
0056
0057
0058 class NuclearInteraction : public InteractionModel {
0059 public:
0060
0061 NuclearInteraction(const std::string& name, const edm::ParameterSet& cfg);
0062
0063
0064 ~NuclearInteraction() override;
0065
0066
0067
0068
0069
0070
0071
0072
0073 void interact(fastsim::Particle& particle,
0074 const SimplifiedGeometry& layer,
0075 std::vector<std::unique_ptr<fastsim::Particle> >& secondaries,
0076 const RandomEngineAndDistribution& random) override;
0077
0078 private:
0079
0080 unsigned index(int thePid);
0081
0082
0083 XYZVector orthogonal(const XYZVector& aVector) const;
0084
0085
0086 void save();
0087
0088
0089 bool read(std::string inputFile);
0090
0091 double theDistCut;
0092 double theHadronEnergy;
0093 std::string inputFile;
0094
0095
0096
0097
0098
0099 TFile* theFile = nullptr;
0100 std::vector<std::vector<TTree*> > theTrees;
0101 std::vector<std::vector<TBranch*> > theBranches;
0102 std::vector<std::vector<NUEvent*> > theNUEvents;
0103 std::vector<std::vector<unsigned> > theCurrentEntry;
0104 std::vector<std::vector<unsigned> > theCurrentInteraction;
0105 std::vector<std::vector<unsigned> > theNumberOfEntries;
0106 std::vector<std::vector<unsigned> > theNumberOfInteractions;
0107 std::vector<std::vector<std::string> > theFileNames;
0108 std::vector<std::vector<double> > theHadronCM;
0109
0110 unsigned ien4;
0111
0112 std::ofstream myOutputFile;
0113 unsigned myOutputBuffer;
0114
0115 bool currentValuesWereSet;
0116
0117
0118
0119
0120
0121
0122 static std::vector<std::vector<double> > theRatiosMap;
0123
0124 static std::map<int, int> theIDMap;
0125
0126
0127 const std::vector<double> theHadronEN = {
0128 1.0, 2.0, 3.0, 4.0, 5.0, 7.0, 9.0, 12.0, 15.0, 20.0, 30.0, 50.0, 100.0, 200.0, 300.0, 500.0, 700.0, 1000.0};
0129
0130
0131
0132
0133 const std::vector<int> theHadronID = {211, -211, 130, 321, -321, 2212, -2212, 2112, -2112};
0134
0135 const std::vector<std::string> theHadronNA = {
0136 "piplus", "piminus", "K0L", "Kplus", "Kminus", "p", "pbar", "n", "nbar"};
0137
0138 const std::vector<double> theHadronMA = {
0139 0.13957, 0.13957, 0.497648, 0.493677, 0.493677, 0.93827, 0.93827, 0.939565, 0.939565};
0140
0141 const std::vector<double> theHadronPMin = {0.7, 0.0, 1.0, 1.0, 0.0, 1.1, 0.0, 1.1, 0.0};
0142
0143 const std::vector<double> theLengthRatio = {
0144 0.2257,
0145 0.2294,
0146 0.3042,
0147 0.2591,
0148 0.2854,
0149 0.3101,
0150 0.5216,
0151 0.3668,
0152 0.4898};
0153
0154 const std::vector<double> theRatios = {
0155 0.031390573,
0156 0.531842852,
0157 0.819614219,
0158 0.951251711,
0159 0.986382750,
0160 1.000000000,
0161 0.985087033,
0162 0.982996773,
0163 0.990832192,
0164 0.992237923,
0165 0.994841580,
0166 0.973816742,
0167 0.967264815,
0168 0.971714258,
0169 0.969122824,
0170 0.978681792,
0171 0.977312732,
0172 0.984255819,
0173
0174 0.035326512,
0175 0.577356403,
0176 0.857118809,
0177 0.965683504,
0178 0.989659360,
0179 1.000000000,
0180 0.989599240,
0181 0.980665408,
0182 0.988384816,
0183 0.981038152,
0184 0.975002104,
0185 0.959996152,
0186 0.953310808,
0187 0.954705592,
0188 0.957615400,
0189 0.961150456,
0190 0.965022184,
0191 0.960573304,
0192
0193 0.000000000,
0194 0.370261189,
0195 0.649793096,
0196 0.734342408,
0197 0.749079499,
0198 0.753360057,
0199 0.755790543,
0200 0.755872164,
0201 0.751337674,
0202 0.746685288,
0203 0.747519634,
0204 0.739357554,
0205 0.735004444,
0206 0.803039922,
0207 0.832749896,
0208 0.890900187,
0209 0.936734805,
0210 1.000000000,
0211
0212 0.000000000,
0213 0.175571717,
0214 0.391683394,
0215 0.528946472,
0216 0.572818635,
0217 0.614210280,
0218 0.644125538,
0219 0.670304050,
0220 0.685144573,
0221 0.702870161,
0222 0.714708513,
0223 0.730805263,
0224 0.777711536,
0225 0.831090576,
0226 0.869267129,
0227 0.915747562,
0228 0.953370523,
0229 1.000000000,
0230
0231 0.000000000,
0232 0.365353210,
0233 0.611663677,
0234 0.715315908,
0235 0.733498956,
0236 0.738361302,
0237 0.745253654,
0238 0.751459671,
0239 0.750628335,
0240 0.746442657,
0241 0.750850669,
0242 0.744895986,
0243 0.735093960,
0244 0.791663444,
0245 0.828609543,
0246 0.889993040,
0247 0.940897842,
0248 1.000000000,
0249
0250 0.000000000,
0251 0.042849136,
0252 0.459103223,
0253 0.666165343,
0254 0.787930873,
0255 0.890397011,
0256 0.920999533,
0257 0.937832788,
0258 0.950920131,
0259 0.966595049,
0260 0.979542270,
0261 0.988061653,
0262 0.983260159,
0263 0.988958431,
0264 0.991723494,
0265 0.995273237,
0266 1.000000000,
0267 0.999962634,
0268
0269 1.000000000,
0270 0.849956907,
0271 0.775625988,
0272 0.802018230,
0273 0.816207485,
0274 0.785899785,
0275 0.754998487,
0276 0.728977244,
0277 0.710010673,
0278 0.670890339,
0279 0.665627872,
0280 0.652682888,
0281 0.613334247,
0282 0.647534574,
0283 0.667910938,
0284 0.689919693,
0285 0.709200185,
0286 0.724199928,
0287
0288 0.000000000,
0289 0.059216484,
0290 0.437844536,
0291 0.610370629,
0292 0.702090648,
0293 0.780076890,
0294 0.802143073,
0295 0.819570432,
0296 0.825829666,
0297 0.840079750,
0298 0.838435509,
0299 0.837529986,
0300 0.835687165,
0301 0.885205014,
0302 0.912450156,
0303 0.951451221,
0304 0.973215562,
0305 1.000000000,
0306
0307 1.000000000,
0308 0.849573257,
0309 0.756479495,
0310 0.787147094,
0311 0.804572414,
0312 0.791806302,
0313 0.760234588,
0314 0.741109531,
0315 0.724118186,
0316 0.692829761,
0317 0.688465897,
0318 0.671806061,
0319 0.636461171,
0320 0.675314029,
0321 0.699134460,
0322 0.724305037,
0323 0.742556115,
0324 0.758504713};
0325
0326
0327
0328
0329
0330 const std::vector<int> protonsID = {2212, 3222, -101, -102, -103, -104};
0331 const std::vector<int> antiprotonsID = {-2212, -3222};
0332 const std::vector<int> neutronsID = {2112, 3122, 3112, 3312, 3322, 3334, -3334};
0333 const std::vector<int> antineutronsID = {-2112, -3122, -3112, -3312, -3322};
0334 const std::vector<int> K0LsID = {130, 310};
0335 const std::vector<int> KplussesID = {321};
0336 const std::vector<int> KminussesID = {-321};
0337 const std::vector<int> PiplussesID = {211};
0338 const std::vector<int> PiminussesID = {-211};
0339 };
0340
0341
0342
0343 static std::once_flag initializeOnce;
0344 CMS_THREAD_GUARD(initializeOnce) std::vector<std::vector<double> > NuclearInteraction::theRatiosMap;
0345 CMS_THREAD_GUARD(initializeOnce) std::map<int, int> NuclearInteraction::theIDMap;
0346 }
0347
0348 fastsim::NuclearInteraction::NuclearInteraction(const std::string& name, const edm::ParameterSet& cfg)
0349 : fastsim::InteractionModel(name), currentValuesWereSet(false) {
0350
0351 std::string fullPath;
0352
0353
0354 theDistCut = cfg.getParameter<double>("distCut");
0355 theHadronEnergy = cfg.getParameter<double>("hadronEnergy");
0356 inputFile = cfg.getUntrackedParameter<std::string>("inputFile", "");
0357
0358
0359
0360 std::call_once(initializeOnce, [this]() {
0361 theRatiosMap.resize(theHadronID.size());
0362 for (unsigned i = 0; i < theHadronID.size(); ++i) {
0363 for (unsigned j = 0; j < theHadronEN.size(); ++j) {
0364 theRatiosMap[i].push_back(theRatios[i * theHadronEN.size() + j]);
0365 }
0366 }
0367
0368
0369
0370 for (const auto& id : protonsID)
0371 theIDMap[id] = 2212;
0372
0373 for (const auto& id : antiprotonsID)
0374 theIDMap[id] = -2212;
0375
0376 for (const auto& id : neutronsID)
0377 theIDMap[id] = 2112;
0378
0379 for (const auto& id : antineutronsID)
0380 theIDMap[id] = -2112;
0381
0382 for (const auto& id : K0LsID)
0383 theIDMap[id] = 130;
0384
0385 for (const auto& id : KplussesID)
0386 theIDMap[id] = 321;
0387
0388 for (const auto& id : KminussesID)
0389 theIDMap[id] = -321;
0390
0391 for (const auto& id : PiplussesID)
0392 theIDMap[id] = 211;
0393
0394 for (const auto& id : PiminussesID)
0395 theIDMap[id] = -211;
0396 });
0397
0398
0399
0400 TFile* aVFile = nullptr;
0401 std::vector<TTree*> aVTree(theHadronEN.size());
0402 std::vector<TBranch*> aVBranch(theHadronEN.size());
0403 std::vector<NUEvent*> aVNUEvents(theHadronEN.size());
0404 std::vector<unsigned> aVCurrentEntry(theHadronEN.size());
0405 std::vector<unsigned> aVCurrentInteraction(theHadronEN.size());
0406 std::vector<unsigned> aVNumberOfEntries(theHadronEN.size());
0407 std::vector<unsigned> aVNumberOfInteractions(theHadronEN.size());
0408 std::vector<std::string> aVFileName(theHadronEN.size());
0409 std::vector<double> aVHadronCM(theHadronEN.size());
0410 theTrees.resize(theHadronNA.size());
0411 theBranches.resize(theHadronNA.size());
0412 theNUEvents.resize(theHadronNA.size());
0413 theCurrentEntry.resize(theHadronNA.size());
0414 theCurrentInteraction.resize(theHadronNA.size());
0415 theNumberOfEntries.resize(theHadronNA.size());
0416 theNumberOfInteractions.resize(theHadronNA.size());
0417 theFileNames.resize(theHadronNA.size());
0418 theHadronCM.resize(theHadronNA.size());
0419 theFile = aVFile;
0420 for (unsigned iname = 0; iname < theHadronNA.size(); ++iname) {
0421 theTrees[iname] = aVTree;
0422 theBranches[iname] = aVBranch;
0423 theNUEvents[iname] = aVNUEvents;
0424 theCurrentEntry[iname] = aVCurrentEntry;
0425 theCurrentInteraction[iname] = aVCurrentInteraction;
0426 theNumberOfEntries[iname] = aVNumberOfEntries;
0427 theNumberOfInteractions[iname] = aVNumberOfInteractions;
0428 theFileNames[iname] = aVFileName;
0429 theHadronCM[iname] = aVHadronCM;
0430 }
0431
0432
0433 currentValuesWereSet = this->read(inputFile);
0434 if (currentValuesWereSet)
0435 std::cout << "***WARNING*** You are reading nuclear-interaction information from the file " << inputFile
0436 << " created in an earlier run." << std::endl;
0437
0438
0439 myOutputFile.open("NuclearInteractionOutputFile.txt");
0440 myOutputBuffer = 0;
0441
0442
0443 edm::FileInPath myDataFile("FastSimulation/MaterialEffects/data/NuclearInteractions.root");
0444 fullPath = myDataFile.fullPath();
0445 theFile = TFile::Open(fullPath.c_str());
0446
0447
0448 for (unsigned iname = 0; iname < theHadronNA.size(); ++iname) {
0449 for (unsigned iene = 0; iene < theHadronEN.size(); ++iene) {
0450 std::ostringstream filename;
0451 double theEne = theHadronEN[iene];
0452 filename << "NuclearInteractionsVal_" << theHadronNA[iname] << "_E" << theEne << ".root";
0453 theFileNames[iname][iene] = filename.str();
0454
0455 std::string treeName = "NuclearInteractions_" + theHadronNA[iname] + "_E" + std::to_string(int(theEne));
0456 theTrees[iname][iene] = (TTree*)theFile->Get(treeName.c_str());
0457
0458 if (!theTrees[iname][iene])
0459 throw cms::Exception("FastSimulation/MaterialEffects") << "Tree with name " << treeName << " not found ";
0460 theBranches[iname][iene] = theTrees[iname][iene]->GetBranch("nuEvent");
0461 if (!theBranches[iname][iene])
0462 throw cms::Exception("FastSimulation/MaterialEffects")
0463 << "Branch with name nuEvent not found in " << theFileNames[iname][iene];
0464
0465 theNUEvents[iname][iene] = new NUEvent();
0466 theBranches[iname][iene]->SetAddress(&theNUEvents[iname][iene]);
0467 theNumberOfEntries[iname][iene] = theTrees[iname][iene]->GetEntries();
0468
0469 if (currentValuesWereSet) {
0470 theTrees[iname][iene]->GetEntry(theCurrentEntry[iname][iene]);
0471 unsigned NInteractions = theNUEvents[iname][iene]->nInteractions();
0472 theNumberOfInteractions[iname][iene] = NInteractions;
0473 }
0474
0475
0476 XYZTLorentzVector Proton(0., 0., 0., 0.986);
0477 XYZTLorentzVector Reference(
0478 0.,
0479 0.,
0480 std::sqrt(theHadronEN[iene] * theHadronEN[iene] - theHadronMA[iname] * theHadronMA[iname]),
0481 theHadronEN[iene]);
0482 theHadronCM[iname][iene] = (Reference + Proton).M();
0483 }
0484 }
0485
0486
0487 ien4 = 0;
0488 while (theHadronEN[ien4] < 4.0)
0489 ++ien4;
0490
0491 gROOT->cd();
0492 }
0493
0494 fastsim::NuclearInteraction::~NuclearInteraction() {
0495
0496
0497
0498 theFile->Close();
0499 delete theFile;
0500
0501 for (auto& vEvents : theNUEvents) {
0502 for (auto evtPtr : vEvents) {
0503 delete evtPtr;
0504 }
0505 }
0506
0507
0508 save();
0509
0510 myOutputFile.close();
0511 }
0512
0513 void fastsim::NuclearInteraction::interact(fastsim::Particle& particle,
0514 const SimplifiedGeometry& layer,
0515 std::vector<std::unique_ptr<fastsim::Particle> >& secondaries,
0516 const RandomEngineAndDistribution& random) {
0517 int pdgId = particle.pdgId();
0518
0519
0520
0521 if (abs(pdgId) <= 100 || abs(pdgId) >= 1000000) {
0522 return;
0523 }
0524
0525 double radLengths = layer.getThickness(particle.position(), particle.momentum());
0526
0527 radLengths *= layer.getNuclearInteractionThicknessFactor();
0528
0529
0530
0531 if (radLengths < 1E-10) {
0532 return;
0533 }
0534
0535
0536 if (!currentValuesWereSet) {
0537 currentValuesWereSet = true;
0538 for (unsigned iname = 0; iname < theHadronNA.size(); ++iname) {
0539 for (unsigned iene = 0; iene < theHadronEN.size(); ++iene) {
0540 theCurrentEntry[iname][iene] = (unsigned)(theNumberOfEntries[iname][iene] * random.flatShoot());
0541
0542 theTrees[iname][iene]->GetEntry(theCurrentEntry[iname][iene]);
0543 unsigned NInteractions = theNUEvents[iname][iene]->nInteractions();
0544 theNumberOfInteractions[iname][iene] = NInteractions;
0545
0546 theCurrentInteraction[iname][iene] = (unsigned)(theNumberOfInteractions[iname][iene] * random.flatShoot());
0547 }
0548 }
0549 }
0550
0551
0552 double pHadron = std::sqrt(particle.momentum().Vect().Mag2());
0553
0554
0555
0556
0557 if (pHadron < theHadronEnergy) {
0558 return;
0559 }
0560
0561
0562 std::map<int, int>::const_iterator thePit = theIDMap.find(pdgId);
0563
0564 int thePid = (thePit != theIDMap.end() ? thePit->second : pdgId);
0565
0566
0567 unsigned fPid = abs(thePid);
0568 if (fPid != 211 && fPid != 130 && fPid != 321 && fPid != 2112 && fPid != 2212) {
0569 return;
0570 }
0571
0572
0573 unsigned thePidIndex = index(thePid);
0574
0575 double theInelasticLength = radLengths * theLengthRatio[thePidIndex];
0576
0577
0578
0579 double ee = pHadron > 0.6 ? exp(-std::sqrt((pHadron - 0.6) / 1.122)) : exp(std::sqrt((0.6 - pHadron) / 1.122));
0580 double theElasticLength = (0.8753 * ee + 0.15) * theInelasticLength;
0581
0582
0583 double theTotalInteractionLength = theInelasticLength + theElasticLength;
0584
0585
0586
0587
0588 double aNuclInteraction = -std::log(random.flatShoot());
0589 if (aNuclInteraction > theTotalInteractionLength) {
0590 return;
0591 }
0592
0593
0594 double elastic = random.flatShoot();
0595 if (elastic < theElasticLength / theTotalInteractionLength) {
0596
0597
0598 double A = 28.0855;
0599 double theta0 = std::sqrt(3.) / std::pow(A, 1. / 3.) * particle.momentum().mass() / pHadron;
0600
0601
0602 double theta = theta0 * std::sqrt(-2. * std::log(random.flatShoot()));
0603 double phi = 2. * M_PI * random.flatShoot();
0604
0605
0606 ROOT::Math::AxisAngle rotation1(orthogonal(particle.momentum().Vect()), theta);
0607 ROOT::Math::AxisAngle rotation2(particle.momentum().Vect(), phi);
0608
0609 XYZVector rotated = rotation2((rotation1(particle.momentum().Vect())));
0610
0611
0612 if (std::sin(theta) > theDistCut) {
0613 secondaries.emplace_back(
0614 new fastsim::Particle(pdgId,
0615 particle.position(),
0616 XYZTLorentzVector(rotated.X(), rotated.Y(), rotated.Z(), particle.momentum().E())));
0617
0618
0619 if (particle.charge() != 0) {
0620 secondaries.back()->setMotherDeltaR(particle.momentum());
0621 secondaries.back()->setMotherPdgId(pdgId);
0622 secondaries.back()->setMotherSimTrackIndex(particle.simTrackIndex());
0623 }
0624
0625
0626 particle.momentum().SetXYZT(0., 0., 0., 0.);
0627 } else {
0628
0629 particle.momentum().SetXYZT(rotated.X(), rotated.Y(), rotated.Z(), particle.momentum().E());
0630 }
0631
0632 } else {
0633
0634 const std::vector<double>& aHadronCM = theHadronCM[thePidIndex];
0635 const std::vector<double>& aRatios = theRatiosMap[thePidIndex];
0636
0637
0638 XYZTLorentzVector Proton(0., 0., 0., 0.939);
0639
0640 const XYZTLorentzVector& Hadron = (const XYZTLorentzVector&)particle.momentum();
0641
0642 double pMin = theHadronPMin[thePidIndex];
0643
0644 XYZTLorentzVector Hadron0(0., 0., pMin, std::sqrt(pMin * pMin + particle.momentum().M2()));
0645
0646
0647 double ecm = (Proton + Hadron).M();
0648
0649
0650 unsigned ene1 = 0;
0651 unsigned ene2 = 0;
0652
0653 double ecm1 = (Proton + Hadron0).M();
0654
0655 double ecm2 = aHadronCM[0];
0656 double ratio1 = 0.;
0657 double ratio2 = aRatios[0];
0658 if (ecm > aHadronCM[0] && ecm < aHadronCM[aHadronCM.size() - 1]) {
0659 for (unsigned ene = 1; ene < aHadronCM.size() && ecm > aHadronCM[ene - 1]; ++ene) {
0660 if (ecm < aHadronCM[ene]) {
0661 ene2 = ene;
0662 ene1 = ene2 - 1;
0663 ecm1 = aHadronCM[ene1];
0664 ecm2 = aHadronCM[ene2];
0665 ratio1 = aRatios[ene1];
0666 ratio2 = aRatios[ene2];
0667 }
0668 }
0669 } else if (ecm > aHadronCM[aHadronCM.size() - 1]) {
0670 ene1 = aHadronCM.size() - 1;
0671 ene2 = aHadronCM.size() - 2;
0672 ecm1 = aHadronCM[ene1];
0673 ecm2 = aHadronCM[ene2];
0674 ratio1 = aRatios[ene2];
0675 ratio2 = aRatios[ene2];
0676 }
0677
0678
0679 double slope = (std::log10(ecm) - std::log10(ecm1)) / (std::log10(ecm2) - std::log10(ecm1));
0680 double inelastic = ratio1 + (ratio2 - ratio1) * slope;
0681 double inelastic4 = pHadron < 4. ? aRatios[ien4] : 1.;
0682
0683
0684 if (elastic > 1. - (inelastic * theInelasticLength) / theTotalInteractionLength) {
0685
0686 std::vector<unsigned>& aCurrentInteraction = theCurrentInteraction[thePidIndex];
0687 std::vector<unsigned>& aNumberOfInteractions = theNumberOfInteractions[thePidIndex];
0688 std::vector<NUEvent*>& aNUEvents = theNUEvents[thePidIndex];
0689
0690
0691
0692
0693 unsigned ene;
0694 if (random.flatShoot() < slope || aNumberOfInteractions[ene1] == 0) {
0695 ene = ene2;
0696 } else {
0697 ene = ene1;
0698 }
0699
0700
0701 XYZTLorentzVector theBoost = Proton + Hadron;
0702 theBoost /= theBoost.e();
0703
0704
0705
0706 if (aCurrentInteraction[ene] == aNumberOfInteractions[ene]) {
0707 std::vector<unsigned>& aCurrentEntry = theCurrentEntry[thePidIndex];
0708 std::vector<unsigned>& aNumberOfEntries = theNumberOfEntries[thePidIndex];
0709 std::vector<TTree*>& aTrees = theTrees[thePidIndex];
0710 ++aCurrentEntry[ene];
0711
0712 aCurrentInteraction[ene] = 0;
0713 if (aCurrentEntry[ene] == aNumberOfEntries[ene]) {
0714 aCurrentEntry[ene] = 0;
0715 }
0716
0717 unsigned myEntry = aCurrentEntry[ene];
0718 aTrees[ene]->GetEntry(myEntry);
0719 aNumberOfInteractions[ene] = aNUEvents[ene]->nInteractions();
0720 }
0721
0722
0723 NUEvent::NUInteraction anInteraction = aNUEvents[ene]->theNUInteractions()[aCurrentInteraction[ene]];
0724
0725 unsigned firstTrack = anInteraction.first;
0726 unsigned lastTrack = anInteraction.last;
0727
0728
0729 XYZVector theAxis = theBoost.Vect().Unit();
0730 double theAngle = random.flatShoot() * 2. * M_PI;
0731 ROOT::Math::AxisAngle axisRotation(theAxis, theAngle);
0732 ROOT::Math::Boost axisBoost(theBoost.x(), theBoost.y(), theBoost.z());
0733
0734
0735 XYZVector zAxis(0., 0., 1.);
0736 XYZVector orthAxis = (zAxis.Cross(theBoost.Vect())).Unit();
0737 double orthAngle = acos(theBoost.Vect().Unit().Z());
0738 ROOT::Math::AxisAngle orthRotation(orthAxis, orthAngle);
0739
0740
0741 for (unsigned iTrack = firstTrack; iTrack <= lastTrack; ++iTrack) {
0742 NUEvent::NUParticle aParticle = aNUEvents[ene]->theNUParticles()[iTrack];
0743
0744
0745
0746 double energy = std::sqrt(aParticle.px * aParticle.px + aParticle.py * aParticle.py +
0747 aParticle.pz * aParticle.pz + aParticle.mass * aParticle.mass / (ecm * ecm));
0748
0749 XYZTLorentzVector daugtherMomentum(aParticle.px * ecm, aParticle.py * ecm, aParticle.pz * ecm, energy * ecm);
0750
0751
0752 XYZVector rotated = orthRotation(daugtherMomentum.Vect());
0753
0754 rotated = axisRotation(rotated);
0755
0756
0757 daugtherMomentum.SetXYZT(rotated.X(), rotated.Y(), rotated.Z(), daugtherMomentum.E());
0758
0759
0760 daugtherMomentum = axisBoost(daugtherMomentum);
0761
0762
0763 secondaries.emplace_back(new fastsim::Particle(aParticle.id, particle.position(), daugtherMomentum));
0764
0765
0766
0767
0768
0769 if (particle.charge() != 0) {
0770 secondaries.back()->setMotherDeltaR(particle.momentum());
0771 secondaries.back()->setMotherPdgId(pdgId);
0772 secondaries.back()->setMotherSimTrackIndex(particle.simTrackIndex());
0773 }
0774 }
0775
0776
0777 particle.momentum().SetXYZT(0., 0., 0., 0.);
0778
0779
0780
0781
0782
0783
0784
0785 ++aCurrentInteraction[ene];
0786
0787
0788 } else if (pHadron < 4. && elastic > 1. - (inelastic4 * theInelasticLength) / theTotalInteractionLength) {
0789
0790 particle.momentum().SetXYZT(0., 0., 0., 0.);
0791 }
0792 }
0793 }
0794
0795 void fastsim::NuclearInteraction::save() {
0796
0797 ++myOutputBuffer;
0798
0799
0800 if (myOutputBuffer / 1000 * 1000 == myOutputBuffer) {
0801 myOutputFile.close();
0802 myOutputFile.open("NuclearInteractionOutputFile.txt");
0803 }
0804
0805 unsigned size1 = theCurrentEntry.size() * theCurrentEntry.begin()->size();
0806 std::vector<unsigned> theCurrentEntries;
0807 theCurrentEntries.resize(size1);
0808 size1 *= sizeof(unsigned);
0809
0810 unsigned size2 = theCurrentInteraction.size() * theCurrentInteraction.begin()->size();
0811 std::vector<unsigned> theCurrentInteractions;
0812 theCurrentInteractions.resize(size2);
0813 size2 *= sizeof(unsigned);
0814
0815
0816 std::vector<std::vector<unsigned> >::const_iterator aCurrentEntry = theCurrentEntry.begin();
0817 std::vector<std::vector<unsigned> >::const_iterator lastCurrentEntry = theCurrentEntry.end();
0818 unsigned allEntries = 0;
0819 for (; aCurrentEntry != lastCurrentEntry; ++aCurrentEntry) {
0820 unsigned size = aCurrentEntry->size();
0821 for (unsigned iene = 0; iene < size; ++iene)
0822 theCurrentEntries[allEntries++] = (*aCurrentEntry)[iene];
0823 }
0824
0825
0826 std::vector<std::vector<unsigned> >::const_iterator aCurrentInteraction = theCurrentInteraction.begin();
0827 std::vector<std::vector<unsigned> >::const_iterator lastCurrentInteraction = theCurrentInteraction.end();
0828 unsigned allInteractions = 0;
0829 for (; aCurrentInteraction != lastCurrentInteraction; ++aCurrentInteraction) {
0830 unsigned size = aCurrentInteraction->size();
0831 for (unsigned iene = 0; iene < size; ++iene)
0832 theCurrentInteractions[allInteractions++] = (*aCurrentInteraction)[iene];
0833 }
0834
0835 myOutputFile.write((const char*)(&theCurrentEntries.front()), size1);
0836 myOutputFile.write((const char*)(&theCurrentInteractions.front()), size2);
0837 myOutputFile.flush();
0838 }
0839
0840 bool fastsim::NuclearInteraction::read(std::string inputFile) {
0841 std::ifstream myInputFile;
0842 struct stat results;
0843
0844 unsigned size1 = theCurrentEntry.size() * theCurrentEntry.begin()->size();
0845 std::vector<unsigned> theCurrentEntries;
0846 theCurrentEntries.resize(size1);
0847 size1 *= sizeof(unsigned);
0848
0849 unsigned size2 = theCurrentInteraction.size() * theCurrentInteraction.begin()->size();
0850 std::vector<unsigned> theCurrentInteractions;
0851 theCurrentInteractions.resize(size2);
0852 size2 *= sizeof(unsigned);
0853
0854 unsigned size = 0;
0855
0856
0857 myInputFile.open(inputFile.c_str());
0858 if (myInputFile.is_open()) {
0859
0860 if (stat(inputFile.c_str(), &results) == 0)
0861 size = results.st_size;
0862 else
0863 return false;
0864
0865
0866 myInputFile.seekg(size - size1 - size2);
0867 myInputFile.read((char*)(&theCurrentEntries.front()), size1);
0868 myInputFile.read((char*)(&theCurrentInteractions.front()), size2);
0869 myInputFile.close();
0870
0871
0872 std::vector<std::vector<unsigned> >::iterator aCurrentEntry = theCurrentEntry.begin();
0873 std::vector<std::vector<unsigned> >::iterator lastCurrentEntry = theCurrentEntry.end();
0874 unsigned allEntries = 0;
0875 for (; aCurrentEntry != lastCurrentEntry; ++aCurrentEntry) {
0876 unsigned size = aCurrentEntry->size();
0877 for (unsigned iene = 0; iene < size; ++iene)
0878 (*aCurrentEntry)[iene] = theCurrentEntries[allEntries++];
0879 }
0880
0881
0882 std::vector<std::vector<unsigned> >::iterator aCurrentInteraction = theCurrentInteraction.begin();
0883 std::vector<std::vector<unsigned> >::iterator lastCurrentInteraction = theCurrentInteraction.end();
0884 unsigned allInteractions = 0;
0885 for (; aCurrentInteraction != lastCurrentInteraction; ++aCurrentInteraction) {
0886 unsigned size = aCurrentInteraction->size();
0887 for (unsigned iene = 0; iene < size; ++iene)
0888 (*aCurrentInteraction)[iene] = theCurrentInteractions[allInteractions++];
0889 }
0890
0891 return true;
0892 }
0893
0894 return false;
0895 }
0896
0897 unsigned fastsim::NuclearInteraction::index(int thePid) {
0898
0899 unsigned myIndex = 0;
0900 while (thePid != theHadronID[myIndex])
0901 ++myIndex;
0902 return myIndex;
0903 }
0904
0905 XYZVector fastsim::NuclearInteraction::orthogonal(const XYZVector& aVector) const {
0906 double x = fabs(aVector.X());
0907 double y = fabs(aVector.Y());
0908 double z = fabs(aVector.Z());
0909
0910 if (x < y)
0911 return x < z ? XYZVector(0., aVector.Z(), -aVector.Y()) : XYZVector(aVector.Y(), -aVector.X(), 0.);
0912 else
0913 return y < z ? XYZVector(-aVector.Z(), 0., aVector.X()) : XYZVector(aVector.Y(), -aVector.X(), 0.);
0914 }
0915
0916 DEFINE_EDM_PLUGIN(fastsim::InteractionModelFactory, fastsim::NuclearInteraction, "fastsim::NuclearInteraction");