File indexing completed on 2023-03-17 11:25:55
0001 #include "SimTransport/PPSProtonTransport/interface/BaseProtonTransport.h"
0002 #include "FWCore/Utilities/interface/Exception.h"
0003 #include "Utilities/PPS/interface/PPSUnitConversion.h"
0004 #include <CLHEP/Random/RandGauss.h>
0005 #include <CLHEP/Units/GlobalSystemOfUnits.h>
0006
0007 BaseProtonTransport::BaseProtonTransport(const edm::ParameterSet& iConfig)
0008 : verbosity_(iConfig.getParameter<bool>("Verbosity")),
0009 bApplyZShift_(iConfig.getParameter<bool>("ApplyZShift")),
0010 useBeamPositionFromLHCInfo_(iConfig.getParameter<bool>("useBeamPositionFromLHCInfo")),
0011 produceHitsRelativeToBeam_(iConfig.getParameter<bool>("produceHitsRelativeToBeam")),
0012 fPPSRegionStart_45_(iConfig.getParameter<double>("PPSRegionStart_45")),
0013 fPPSRegionStart_56_(iConfig.getParameter<double>("PPSRegionStart_56")),
0014 etaCut_(iConfig.getParameter<double>("EtaCut")),
0015 momentumCut_(iConfig.getParameter<double>("MomentumCut")),
0016 beamEnergy_(iConfig.getParameter<double>("BeamEnergy")) {
0017 beamMomentum_ = std::sqrt(beamEnergy_ * beamEnergy_ - ProtonMassSQ);
0018 }
0019
0020 BaseProtonTransport::~BaseProtonTransport() { clear(); }
0021
0022 void BaseProtonTransport::ApplyBeamCorrection(HepMC::GenParticle* p) {
0023 TLorentzVector p_out;
0024 p_out.SetPx(p->momentum().px());
0025 p_out.SetPy(p->momentum().py());
0026 p_out.SetPz(p->momentum().pz());
0027 p_out.SetE(p->momentum().e());
0028 ApplyBeamCorrection(p_out);
0029 p->set_momentum(HepMC::FourVector(p_out.Px(), p_out.Py(), p_out.Pz(), p_out.E()));
0030 }
0031
0032 void BaseProtonTransport::ApplyBeamCorrection(TLorentzVector& p_out) {
0033 double thetax = atan(p_out.Px() / fabs(p_out.Pz()));
0034 double thetay = atan(p_out.Py() / fabs(p_out.Pz()));
0035 double energy = p_out.E();
0036 double urad = 1e-6;
0037
0038 int direction = (p_out.Pz() > 0) ? 1 : -1;
0039
0040 if (MODE == TransportMode::TOTEM)
0041 thetax += (p_out.Pz() > 0) ? fCrossingAngleX_45_ * urad : fCrossingAngleX_56_ * urad;
0042
0043 double dtheta_x = (m_sigmaSTX <= 0.0) ? 0.0 : CLHEP::RandGauss::shoot(engine_, 0., m_sigmaSTX);
0044 double dtheta_y = (m_sigmaSTY <= 0.0) ? 0.0 : CLHEP::RandGauss::shoot(engine_, 0., m_sigmaSTY);
0045 double denergy = (m_sig_E <= 0.0) ? 0.0 : CLHEP::RandGauss::shoot(engine_, 0., m_sig_E);
0046
0047 double s_theta = std::sqrt(pow(thetax + dtheta_x * urad, 2) + std::pow(thetay + dtheta_y * urad, 2));
0048 double s_phi = std::atan2(thetay + dtheta_y * urad, thetax + dtheta_x * urad);
0049 energy += denergy;
0050 double p = std::sqrt(std::pow(energy, 2) - ProtonMassSQ);
0051 double sint = std::sin(s_theta);
0052
0053 p_out.SetPx(p * sint * std::cos(s_phi));
0054 p_out.SetPy(p * sint * std::sin(s_phi));
0055 p_out.SetPz(p * std::cos(s_theta) * direction);
0056 p_out.SetE(energy);
0057 }
0058
0059 void BaseProtonTransport::addPartToHepMC(const HepMC::GenEvent* in_evt, HepMC::GenEvent* evt) {
0060 m_CorrespondenceMap.clear();
0061
0062 int direction = 0;
0063 HepMC::GenParticle* gpart;
0064
0065 unsigned int line;
0066
0067 for (auto const& it : m_beamPart) {
0068 line = (it).first;
0069 gpart = in_evt->barcode_to_particle(line);
0070
0071 direction = (gpart->momentum().pz() > 0) ? 1 : -1;
0072
0073
0074 double ddd = (direction > 0) ? fPPSRegionStart_45_ : fabs(fPPSRegionStart_56_);
0075
0076 double time = (ddd * meter - gpart->production_vertex()->position().z() * mm);
0077
0078
0079
0080
0081 if (ddd == 0.)
0082 continue;
0083
0084 if (verbosity_) {
0085 LogDebug("BaseProtonTransportEventProcessing")
0086 << "BaseProtonTransport:: x= " << (*(m_xAtTrPoint.find(line))).second << "\n"
0087 << "BaseProtonTransport:: y= " << (*(m_yAtTrPoint.find(line))).second << "\n"
0088 << "BaseProtonTransport:: z= " << ddd * direction * m_to_mm << "\n"
0089 << "BaseProtonTransport:: t= " << time;
0090 }
0091 TLorentzVector const& p_out = (it).second;
0092
0093 HepMC::GenVertex* vert = new HepMC::GenVertex(HepMC::FourVector((*(m_xAtTrPoint.find(line))).second,
0094 (*(m_yAtTrPoint.find(line))).second,
0095 ddd * direction * m_to_mm,
0096 time + time * 0.001));
0097
0098 vert->add_particle_out(new HepMC::GenParticle(
0099 HepMC::FourVector(p_out.Px(), p_out.Py(), p_out.Pz(), p_out.E()), gpart->pdg_id(), 1, gpart->flow()));
0100 evt->add_vertex(vert);
0101
0102 int ingoing = 0;
0103 int outgoing = (*vert->particles_out_const_begin())->barcode();
0104
0105 LHCTransportLink theLink(ingoing, outgoing);
0106 if (verbosity_)
0107 LogDebug("BaseProtonTransportEventProcessing")
0108 << "BaseProtonTransport:addPartToHepMC: LHCTransportLink " << theLink;
0109 m_CorrespondenceMap.push_back(theLink);
0110 }
0111 }
0112 void BaseProtonTransport::clear() {
0113 m_beamPart.clear();
0114 m_xAtTrPoint.clear();
0115 m_yAtTrPoint.clear();
0116 m_CorrespondenceMap.clear();
0117 }