PPSProtonTransport/src/BaseProtonTransport.cc

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/SystemOfUnits.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     // Totem uses negative Z for sector 56 while Hector uses always positive distance
0074     double ddd = (direction > 0) ? fPPSRegionStart_45_ : fabs(fPPSRegionStart_56_);
0075
0076     double time = (ddd * CLHEP::meter - gpart->production_vertex()->position().z() * CLHEP::mm);  // mm
0077
0078     //
0079     // ATTENTION: at this point, the vertex at PPS is already in mm
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;  //do not attach the incoming proton to this vertex to avoid duplicating data
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 }