14 #include <Math/AxisAngle.h> 15 #include "Math/Boost.h" 81 unsigned index(
int thePid);
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};
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};
145 0.2257, 0.2294, 0.3042, 0.2591, 0.2854, 0.3101, 0.5216, 0.3668, 0.4898
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,
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,
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,
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,
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,
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,
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,
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,
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
191 const std::vector<int>
protonsID = {2212, 3222, -101, -102, -103, -104};
193 const std::vector<int>
neutronsID = {2112, 3122, 3112, 3312, 3322, 3334, -3334};
195 const std::vector<int>
K0LsID = {130, 310};
224 std::call_once(initializeOnce, [
this] () {
234 for(
const auto &
id :
protonsID) theIDMap[
id] = 2212;
242 for(
const auto &
id :
K0LsID) theIDMap[
id] = 130;
255 TFile* aVFile=
nullptr;
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());
275 for(
unsigned iname=0; iname<
theHadronNA.size(); ++iname){
290 std::cout <<
"***WARNING*** You are reading nuclear-interaction information from the file " 291 << inputFile <<
" created in an earlier run." 299 edm::FileInPath myDataFile(
"FastSimulation/MaterialEffects/data/NuclearInteractions.root");
301 theFile = TFile::Open(fullPath.c_str());
305 for(
unsigned iname=0; iname<
theHadronNA.size(); ++iname){
306 for(
unsigned iene=0; iene<
theHadronEN.size(); ++iene){
309 filename <<
"NuclearInteractionsVal_" <<
theHadronNA[iname] <<
"_E"<< theEne <<
".root";
317 <<
"Tree with name " << treeName <<
" not found ";
320 <<
"Branch with name nuEvent not found in " <<
theFileNames[iname][iene];
328 unsigned NInteractions =
theNUEvents[iname][iene]->nInteractions();
357 for(
auto evtPtr: vEvents){
375 if(
abs(pdgId) <= 100 ||
abs(pdgId) >= 1000000)
386 if(radLengths < 1E-10)
394 for (
unsigned iname=0; iname<
theHadronNA.size(); ++iname ) {
395 for (
unsigned iene=0; iene<
theHadronEN.size(); ++iene ) {
399 unsigned NInteractions =
theNUEvents[iname][iene]->nInteractions();
418 std::map<int,int>::const_iterator thePit =
theIDMap.find(pdgId);
420 int thePid = (thePit !=
theIDMap.end() ? thePit->second :
pdgId);
423 unsigned fPid =
abs(thePid);
424 if(fPid != 211 && fPid != 130 && fPid != 321 && fPid != 2112 && fPid != 2212)
430 unsigned thePidIndex =
index(thePid);
432 double theInelasticLength = radLengths *
theLengthRatio[thePidIndex];
437 double theElasticLength = (0.8753 * ee + 0.15) * theInelasticLength;
440 double theTotalInteractionLength = theInelasticLength + theElasticLength;
446 if(aNuclInteraction > theTotalInteractionLength)
453 if(elastic < theElasticLength/theTotalInteractionLength){
465 ROOT::Math::AxisAngle rotation2(particle.
momentum().Vect(),phi);
479 if(particle.
charge() != 0){
480 secondaries.back()->setMotherDeltaR(particle.
momentum());
481 secondaries.back()->setMotherPdgId(pdgId);
482 secondaries.back()->setMotherSimTrackIndex(particle.
simTrackIndex());
486 particle.
momentum().SetXYZT(0., 0., 0., 0.);
489 particle.
momentum().SetXYZT(rotated.X(),
497 const std::vector<double>& aHadronCM =
theHadronCM[thePidIndex];
498 const std::vector<double>& aRatios =
theRatiosMap[thePidIndex];
513 double ecm = (Proton+Hadron).M();
519 double ecm1= (Proton+Hadron0).M();
521 double ecm2=aHadronCM[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]){
529 ecm1 = aHadronCM[ene1];
530 ecm2 = aHadronCM[ene2];
531 ratio1 = aRatios[ene1];
532 ratio2 = aRatios[ene2];
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];
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.;
550 if(elastic > 1.- (inelastic*theInelasticLength) / theTotalInteractionLength){
554 std::vector<NUEvent*>& aNUEvents =
theNUEvents[thePidIndex];
560 if(random.
flatShoot() < slope || aNumberOfInteractions[ene1] == 0){
568 theBoost /= theBoost.e();
572 if(aCurrentInteraction[ene] == aNumberOfInteractions[ene]){
575 std::vector<TTree*>& aTrees =
theTrees[thePidIndex];
576 ++aCurrentEntry[ene];
578 aCurrentInteraction[ene] = 0;
579 if(aCurrentEntry[ene] == aNumberOfEntries[ene]){
580 aCurrentEntry[ene] = 0;
583 unsigned myEntry = aCurrentEntry[ene];
584 aTrees[ene]->GetEntry(myEntry);
585 aNumberOfInteractions[ene] = aNUEvents[ene]->nInteractions();
591 unsigned firstTrack = anInteraction.
first;
592 unsigned lastTrack = anInteraction.
last;
595 XYZVector theAxis = theBoost.Vect().Unit();
597 ROOT::Math::AxisAngle axisRotation(theAxis,theAngle);
598 ROOT::Math::Boost axisBoost(theBoost.x(),theBoost.y(),theBoost.z());
602 XYZVector orthAxis = (zAxis.Cross(theBoost.Vect())).Unit();
603 double orthAngle = acos(theBoost.Vect().Unit().Z());
604 ROOT::Math::AxisAngle orthRotation(orthAxis,orthAngle);
607 for(
unsigned iTrack=firstTrack; iTrack<=lastTrack; ++iTrack){
613 + aParticle.
py*aParticle.
py 614 + aParticle.
pz*aParticle.
pz 615 + aParticle.
mass*aParticle.
mass / (ecm * ecm));
620 XYZVector rotated = orthRotation(daugtherMomentum.Vect());
622 rotated = axisRotation(rotated);
625 daugtherMomentum.SetXYZT(rotated.X(), rotated.Y(), rotated.Z(), daugtherMomentum.E());
628 daugtherMomentum = axisBoost(daugtherMomentum);
637 if(particle.
charge() != 0){
638 secondaries.back()->setMotherDeltaR(particle.
momentum());
639 secondaries.back()->setMotherPdgId(pdgId);
640 secondaries.back()->setMotherSimTrackIndex(particle.
simTrackIndex());
645 particle.
momentum().SetXYZT(0., 0., 0., 0.);
653 ++aCurrentInteraction[ene];
656 }
else if(pHadron < 4. && elastic > 1.- (inelastic4*theInelasticLength) / theTotalInteractionLength){
658 particle.
momentum().SetXYZT(0., 0., 0., 0.);
676 std::vector<unsigned> theCurrentEntries;
677 theCurrentEntries.resize(size1);
678 size1*=
sizeof(unsigned);
681 std::vector<unsigned> theCurrentInteractions;
682 theCurrentInteractions.resize(size2);
683 size2 *=
sizeof(unsigned);
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];
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];
705 myOutputFile.write((
const char*)(&theCurrentEntries.front()), size1);
706 myOutputFile.write((
const char*)(&theCurrentInteractions.front()), size2);
713 std::ifstream myInputFile;
717 std::vector<unsigned> theCurrentEntries;
718 theCurrentEntries.resize(size1);
719 size1*=
sizeof(unsigned);
722 std::vector<unsigned> theCurrentInteractions;
723 theCurrentInteractions.resize(size2);
724 size2 *=
sizeof(unsigned);
730 myInputFile.open (inputFile.c_str());
731 if(myInputFile.is_open()){
734 if(
stat(inputFile.c_str(), &results) == 0) size = results.st_size;
738 myInputFile.seekg(size-size1-size2);
739 myInputFile.read((
char*)(&theCurrentEntries.front()),size1);
740 myInputFile.read((
char*)(&theCurrentInteractions.front()),size2);
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++];
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++];
781 double x = fabs(aVector.X());
782 double y = fabs(aVector.Y());
783 double z = fabs(aVector.Z());
799 "fastsim::NuclearInteraction" const std::vector< int > antineutronsID
PdgID of anti-neutrons.
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.
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)
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.
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.
Base class for any interaction model between a particle and a tracker layer.
const std::vector< double > theRatios
Filled into 'theRatiosMap' (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.
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.
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)
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.
const std::vector< int > neutronsID
PdgID of neutrons.
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
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.
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-.
const std::vector< double > theHadronEN
Filled into 'theRatiosMap' (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
const math::XYZTLorentzVector & momentum() const
Return momentum of the particle.
#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...
math::XYZVector XYZVector
const std::vector< std::string > theHadronNA
Names of the hadrons.
Power< A, B >::type pow(const A &a, const B &b)
math::XYZTLorentzVector XYZTLorentzVector
unsigned myOutputBuffer
Needed to save interactions to file.