File indexing completed on 2021-02-14 13:26:41
0001
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
0012
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
0022
0023
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
0053
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
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
0094 myOutputFile.open("NuclearInteractionOutputFile.txt");
0095 myOutputBuffer = 0;
0096
0097
0098
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
0108 std::ostringstream filename;
0109 double theEne = thePionEN[iene];
0110 filename << "NuclearInteractionsVal_" << thePionNA[iname] << "_E" << theEne << ".root";
0111 theFileNames[iname][iene] = filename.str();
0112
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
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
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
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
0149 ien4 = 0;
0150 while (thePionEN[ien4] < 4.0)
0151 ++ien4;
0152
0153 gROOT->cd();
0154
0155
0156
0157
0158
0159
0160
0161
0162
0163
0164
0165
0166
0167
0168 }
0169
0170 NuclearInteractionSimulator::~NuclearInteractionSimulator() {
0171
0172
0173
0174 theFile->Close();
0175 delete theFile;
0176
0177 for (auto& vEvents : theNUEvents) {
0178 for (auto evtPtr : vEvents) {
0179 delete evtPtr;
0180 }
0181 }
0182
0183
0184 myOutputFile.close();
0185
0186
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
0206
0207 double pHadron = std::sqrt(Particle.particle().Vect().Mag2());
0208
0209
0210
0211 if (pHadron > thePionEnergy) {
0212
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
0218 unsigned fPid = abs(thePid);
0219 if (fPid != 211 && fPid != 130 && fPid != 321 && fPid != 2112 && fPid != 2212) {
0220 return;
0221
0222
0223 }
0224
0225
0226 unsigned thePidIndex = index(thePid);
0227 double theInelasticLength = radLengths * theLengthRatio[thePidIndex];
0228
0229
0230
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
0234
0235 * theInelasticLength;
0236
0237
0238 double theTotalInteractionLength = theInelasticLength + theElasticLength;
0239
0240
0241 double aNuclInteraction = -std::log(random->flatShoot());
0242 if (aNuclInteraction < theTotalInteractionLength) {
0243
0244 double elastic = random->flatShoot();
0245 if (elastic < theElasticLength / theTotalInteractionLength) {
0246
0247
0248
0249 double theta0 = std::sqrt(3.) / std::pow(theA(), 1. / 3.) * Particle.particle().mass() / pHadron;
0250
0251
0252 double theta = theta0 * std::sqrt(-2. * std::log(random->flatShoot()));
0253 double phi = 2. * 3.14159265358979323 * random->flatShoot();
0254
0255
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
0262 double distance = std::sin(theta);
0263
0264
0265 if (distance > theDistCut) {
0266 _theUpdatedState.reserve(1);
0267 _theUpdatedState.clear();
0268 _theUpdatedState.emplace_back(Particle.particle());
0269 }
0270
0271
0272
0273
0274 }
0275
0276
0277 else {
0278
0279 const std::vector<double>& aPionCM = thePionCM[thePidIndex];
0280 const std::vector<double>& aRatios = theRatios[thePidIndex];
0281
0282
0283 XYZTLorentzVector Proton(0., 0., 0., 0.939);
0284
0285 const XYZTLorentzVector& Hadron = (const XYZTLorentzVector&)Particle;
0286
0287 double pMin = thePionPMin[thePidIndex];
0288
0289 XYZTLorentzVector Hadron0(0., 0., pMin, std::sqrt(pMin * pMin + Particle.particle().M2()));
0290
0291
0292 double ecm = (Proton + Hadron).M();
0293
0294
0295
0296
0297 unsigned ene1 = 0;
0298 unsigned ene2 = 0;
0299
0300
0301 double ecm1 = (Proton + Hadron0).M();
0302
0303
0304
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
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
0334
0335
0336
0337
0338 if (elastic > 1. - (inelastic * theInelasticLength) / theTotalInteractionLength) {
0339
0340 std::vector<unsigned>& aCurrentInteraction = theCurrentInteraction[thePidIndex];
0341 std::vector<unsigned>& aNumberOfInteractions = theNumberOfInteractions[thePidIndex];
0342 std::vector<NUEvent*>& aNUEvents = theNUEvents[thePidIndex];
0343
0344
0345
0346
0347
0348
0349
0350 unsigned ene;
0351 if (random->flatShoot() < slope || aNumberOfInteractions[ene1] == 0)
0352 ene = ene2;
0353 else
0354 ene = ene1;
0355
0356
0357
0358
0359
0360
0361
0362
0363
0364
0365 XYZTLorentzVector theBoost = Proton + Hadron;
0366 theBoost /= theBoost.e();
0367
0368
0369
0370
0371
0372
0373
0374
0375
0376
0377 if (aCurrentInteraction[ene] == aNumberOfInteractions[ene]) {
0378
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
0384
0385 aCurrentInteraction[ene] = 0;
0386 if (aCurrentEntry[ene] == aNumberOfEntries[ene]) {
0387 aCurrentEntry[ene] = 0;
0388
0389 }
0390 unsigned myEntry = aCurrentEntry[ene];
0391
0392
0393 aTrees[ene]->GetEntry(myEntry);
0394
0395 aNumberOfInteractions[ene] = aNUEvents[ene]->nInteractions();
0396
0397 }
0398
0399
0400 NUEvent::NUInteraction anInteraction = aNUEvents[ene]->theNUInteractions()[aCurrentInteraction[ene]];
0401
0402 unsigned firstTrack = anInteraction.first;
0403 unsigned lastTrack = anInteraction.last;
0404
0405
0406 _theUpdatedState.reserve(lastTrack - firstTrack + 1);
0407 _theUpdatedState.clear();
0408
0409 double distMin = 1E99;
0410
0411
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
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
0424
0425
0426
0427 for (unsigned iTrack = firstTrack; iTrack <= lastTrack; ++iTrack) {
0428 unsigned idaugh = iTrack - firstTrack;
0429 NUEvent::NUParticle aParticle = aNUEvents[ene]->theNUParticles()[iTrack];
0430
0431
0432
0433
0434
0435
0436
0437
0438
0439
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
0448 aDaughter.rotate(orthRotation);
0449
0450
0451 aDaughter.rotate(axisRotation);
0452
0453
0454 aDaughter.boost(axisBoost);
0455
0456
0457 double distance = distanceToPrimary(Particle.particle(), aDaughter);
0458
0459 if (distance < distMin && distance < theDistCut) {
0460 distMin = distance;
0461 theClosestChargedDaughterId = idaugh;
0462 }
0463
0464
0465 }
0466
0467
0468
0469
0470
0471
0472
0473
0474
0475
0476
0477
0478
0479
0480
0481
0482 ++aCurrentInteraction[ene];
0483
0484
0485 } else if (pHadron < 4. && elastic > 1. - (inelastic4 * theInelasticLength) / theTotalInteractionLength) {
0486
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
0501 if (fabs(Particle.charge()) > 1E-12) {
0502
0503 double chargeDiff = fabs(aDaughter.charge() - Particle.charge());
0504 if (fabs(chargeDiff) < 1E-12) {
0505
0506 switch (theDistAlgo) {
0507 case 1:
0508
0509 distance = (aDaughter.Vect().Unit().Cross(Particle.Vect().Unit())).R();
0510 break;
0511
0512 case 2:
0513
0514 distance = (aDaughter.Vect().Cross(Particle.Vect())).R() / aDaughter.Vect().Mag2();
0515 break;
0516
0517 default:
0518
0519 distance = 2E99;
0520 break;
0521 }
0522 }
0523 }
0524
0525 return distance;
0526 }
0527
0528 void NuclearInteractionSimulator::save() {
0529
0530 ++myOutputBuffer;
0531
0532
0533 if (myOutputBuffer / 1000 * 1000 == myOutputBuffer) {
0534 myOutputFile.close();
0535 myOutputFile.open("NuclearInteractionOutputFile.txt");
0536
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
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
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
0591 myInputFile.open(inputFile.c_str());
0592 if (myInputFile.is_open()) {
0593
0594 if (stat(inputFile.c_str(), &results) == 0)
0595 size = results.st_size;
0596 else
0597 return false;
0598
0599
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
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
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
0636 return myIndex;
0637 }