PPSProtonTransport/src/OpticalFunctionsTransport.cc

0001 #include "SimTransport/PPSProtonTransport/interface/OpticalFunctionsTransport.h"
0002 #include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"
0003
0004 OpticalFunctionsTransport::OpticalFunctionsTransport(const edm::ParameterSet& iConfig, edm::ConsumesCollector iC)
0005     : BaseProtonTransport(iConfig),
0006       lhcInfoToken_(iC.esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("lhcInfoLabel")))),
0007       beamParametersToken_(iC.esConsumes()),
0008       opticsToken_(iC.esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("opticsLabel")))),
0009       beamspotToken_(iC.esConsumes()),
0010       useEmpiricalApertures_(iConfig.getParameter<bool>("useEmpiricalApertures")),
0011       empiricalAperture45_xi0_int_(iConfig.getParameter<double>("empiricalAperture45_xi0_int")),
0012       empiricalAperture45_xi0_slp_(iConfig.getParameter<double>("empiricalAperture45_xi0_slp")),
0013       empiricalAperture45_a_int_(iConfig.getParameter<double>("empiricalAperture45_a_int")),
0014       empiricalAperture45_a_slp_(iConfig.getParameter<double>("empiricalAperture45_a_slp")),
0015       empiricalAperture56_xi0_int_(iConfig.getParameter<double>("empiricalAperture56_xi0_int")),
0016       empiricalAperture56_xi0_slp_(iConfig.getParameter<double>("empiricalAperture56_xi0_slp")),
0017       empiricalAperture56_a_int_(iConfig.getParameter<double>("empiricalAperture56_a_int")),
0018       empiricalAperture56_a_slp_(iConfig.getParameter<double>("empiricalAperture56_a_slp")) {
0019   MODE = TransportMode::OPTICALFUNCTIONS;
0020 }
0021 OpticalFunctionsTransport::~OpticalFunctionsTransport() {}
0022
0023 void OpticalFunctionsTransport::process(const HepMC::GenEvent* evt,
0024                                         const edm::EventSetup& iSetup,
0025                                         CLHEP::HepRandomEngine* engine) {
0026   clear();
0027
0028   lhcInfo_ = &iSetup.getData(lhcInfoToken_);
0029   beamParameters_ = &iSetup.getData(beamParametersToken_);
0030   opticalFunctions_ = &iSetup.getData(opticsToken_);
0031   beamspot_ = &iSetup.getData(beamspotToken_);
0032
0033   // Choose the optical function corresponding to the first station ono each side (it is in lhc ref. frame)
0034   optFunctionId45_ = 0;
0035   optFunctionId56_ = 0;
0036   for (const auto& ofp : (*opticalFunctions_)) {
0037     if (ofp.second.getScoringPlaneZ() < 0) {
0038       if (optFunctionId45_ == 0)
0039         optFunctionId45_ = ofp.first;
0040       if (opticalFunctions_->at(optFunctionId45_).getScoringPlaneZ() < ofp.second.getScoringPlaneZ())
0041         optFunctionId45_ = ofp.first;
0042     }
0043     if (ofp.second.getScoringPlaneZ() > 0) {
0044       if (optFunctionId56_ == 0)
0045         optFunctionId56_ = ofp.first;
0046       if (opticalFunctions_->at(optFunctionId56_).getScoringPlaneZ() > ofp.second.getScoringPlaneZ())
0047         optFunctionId56_ = ofp.first;
0048     }
0049   }
0050   //
0051   engine_ = engine;  // the engine needs to be updated for each event
0052
0053   for (HepMC::GenEvent::particle_const_iterator eventParticle = evt->particles_begin();
0054        eventParticle != evt->particles_end();
0055        ++eventParticle) {
0056     if (!((*eventParticle)->status() == 1 && (*eventParticle)->pdg_id() == 2212))
0057       continue;
0058
0059     if (!(fabs((*eventParticle)->momentum().eta()) > etaCut_ && fabs((*eventParticle)->momentum().pz()) > momentumCut_))
0060       continue;  // discard protons outside kinematic acceptance
0061
0062     unsigned int line = (*eventParticle)->barcode();
0063     HepMC::GenParticle* gpart = (*eventParticle);
0064     if (gpart->pdg_id() != 2212)
0065       continue;  // only transport stable protons
0066     if (gpart->status() != 1)
0067       continue;
0068     if (m_beamPart.find(line) != m_beamPart.end())  // assures this protons has not been already propagated
0069       continue;
0070
0071     transportProton(gpart);
0072   }
0073 }
0074
0075 bool OpticalFunctionsTransport::transportProton(const HepMC::GenParticle* in_trk) {
0076   const HepMC::FourVector& vtx_cms = in_trk->production_vertex()->position();  // in mm
0077   const HepMC::FourVector& mom_cms = in_trk->momentum();
0078
0079   // transformation to LHC/TOTEM convention
0080   HepMC::FourVector vtx_lhc(-vtx_cms.x(), vtx_cms.y(), -vtx_cms.z(), vtx_cms.t());
0081   HepMC::FourVector mom_lhc(-mom_cms.x(), mom_cms.y(), -mom_cms.z(), mom_cms.t());
0082
0083   // determine the LHC arm and related parameters
0084   double urad = 1.e-6;
0085   double beamMomentum = 0.;
0086   double xangle = 0.;
0087   double empiricalAperture_xi0_int, empiricalAperture_xi0_slp;
0088   double empiricalAperture_a_int, empiricalAperture_a_slp;
0089   unsigned int optFunctionId;
0090   // get the beam position at the IP in mm and in the LHC ref. frame
0091   double vtxXoffset;
0092   double vtxYoffset;
0093   if (useBeamPositionFromLHCInfo_) {
0094     vtxXoffset = -beamParameters_->getVtxOffsetX45() * cm_to_mm;
0095     vtxYoffset = beamParameters_->getVtxOffsetY45() * cm_to_mm;
0096   } else {
0097     vtxXoffset = -beamspot_->x() * cm_to_mm;
0098     vtxYoffset = beamspot_->y() * cm_to_mm;
0099   }
0100
0101   if (mom_lhc.z() < 0)  // sector 45
0102   {
0103     optFunctionId = optFunctionId45_;
0104     beamMomentum = beamParameters_->getBeamMom45();
0105     xangle = beamParameters_->getHalfXangleX45();
0106     empiricalAperture_xi0_int = empiricalAperture45_xi0_int_;
0107     empiricalAperture_xi0_slp = empiricalAperture45_xi0_slp_;
0108     empiricalAperture_a_int = empiricalAperture45_a_int_;
0109     empiricalAperture_a_slp = empiricalAperture45_a_slp_;
0110   } else {  // sector 56
0111     optFunctionId = optFunctionId56_;
0112     beamMomentum = beamParameters_->getBeamMom56();
0113     xangle = beamParameters_->getHalfXangleX56();
0114     empiricalAperture_xi0_int = empiricalAperture56_xi0_int_;
0115     empiricalAperture_xi0_slp = empiricalAperture56_xi0_slp_;
0116     empiricalAperture_a_int = empiricalAperture56_a_int_;
0117     empiricalAperture_a_slp = empiricalAperture56_a_slp_;
0118   }
0119   if (xangle > 1.0)
0120     xangle *= urad;
0121   // calculate kinematics for optics parametrisation, avoid the aproximation for small angles xangle -> tan(xangle)
0122   const double p = mom_lhc.rho();
0123   const double xi = 1. - p / beamMomentum;
0124   const double th_x_phys = mom_lhc.x() / abs(mom_lhc.z()) - tan(xangle);  //"-" in the LHC ref. frame
0125   const double th_y_phys = mom_lhc.y() / abs(mom_lhc.z());
0126   const double vtx_lhc_eff_x = vtx_lhc.x() - vtx_lhc.z() * (mom_lhc.x() / mom_lhc.z() + tan(xangle)) - (vtxXoffset);
0127   const double vtx_lhc_eff_y = vtx_lhc.y() - vtx_lhc.z() * (mom_lhc.y() / mom_lhc.z()) - (vtxYoffset);
0128
0129   if (verbosity_) {
0130     LogDebug("OpticalFunctionsTransport")
0131         << "simu: xi = " << xi << ", th_x_phys = " << th_x_phys << ", th_y_phys = " << th_y_phys
0132         << ", vtx_lhc_eff_x = " << vtx_lhc_eff_x << ", vtx_lhc_eff_y = " << vtx_lhc_eff_y;
0133   }
0134
0135   // check empirical aperture
0136   if (useEmpiricalApertures_) {
0137     const auto& xangle =
0138         (lhcInfo_->crossingAngle() > 1.0) ? lhcInfo_->crossingAngle() * urad : lhcInfo_->crossingAngle();
0139     const double xi_th = (empiricalAperture_xi0_int + xangle * empiricalAperture_xi0_slp) +
0140                          (empiricalAperture_a_int + xangle * empiricalAperture_a_slp) * th_x_phys;
0141
0142     if (xi > xi_th) {
0143       if (verbosity_) {
0144         LogDebug("OpticalFunctionsTransport") << "stop because of empirical appertures";
0145       }
0146       return false;
0147     }
0148   }
0149
0150   // transport the proton into  pot/scoring plane
0151   auto ofp = opticalFunctions_->at(optFunctionId);
0152   CTPPSDetId rpId(optFunctionId);
0153   const unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0154
0155   if (verbosity_)
0156     LogDebug("OpticalFunctionsTransport") << "  RP " << rpDecId << std::endl;
0157
0158   // transport proton
0159   LHCInterpolatedOpticalFunctionsSet::Kinematics k_in = {
0160       vtx_lhc_eff_x * cm_to_mm, th_x_phys, vtx_lhc_eff_y * cm_to_mm, th_y_phys, xi};  // conversions: mm -> cm
0161
0162   LHCInterpolatedOpticalFunctionsSet::Kinematics k_out;
0163   ofp.transport(k_in, k_out, true);
0164
0165   // Original code uses mm, but CMS uses cm, so keep it in cm
0166   double b_x = k_out.x * cm_to_mm, b_y = k_out.y * cm_to_mm;  // conversions: cm -> mm
0167   double a_x = k_out.th_x, a_y = k_out.th_y;
0168
0169   // if needed, subtract beam position and angle
0170   if (produceHitsRelativeToBeam_) {
0171     // determine beam position
0172     LHCInterpolatedOpticalFunctionsSet::Kinematics k_be_in = {0., -tan(xangle), 0., 0., 0.};
0173     LHCInterpolatedOpticalFunctionsSet::Kinematics k_be_out;
0174     ofp.transport(k_be_in, k_be_out, true);
0175
0176     a_x -= k_be_out.th_x;
0177     a_y -= k_be_out.th_y;
0178     b_x -= k_be_out.x * cm_to_mm;
0179     b_y -= k_be_out.y * cm_to_mm;  // conversions: cm -> mm
0180   }
0181
0182   const double z_scoringPlane = ofp.getScoringPlaneZ() * cm_to_mm;  // conversion: cm --> mm
0183
0184   if (verbosity_) {
0185     LogDebug("OpticalFunctionsTransport")
0186         << "    proton transported: a_x = " << a_x << " rad, a_y = " << a_y << " rad, b_x = " << b_x
0187         << " mm, b_y = " << b_y << " mm, z = " << z_scoringPlane << " mm";
0188   }
0189   //
0190   // Project the track back to the starting of PPS region in mm
0191   b_x -= (abs(z_scoringPlane) - ((z_scoringPlane < 0) ? fPPSRegionStart_45_ : fPPSRegionStart_56_) * 1e3) * a_x;
0192   b_y -= (abs(z_scoringPlane) - ((z_scoringPlane < 0) ? fPPSRegionStart_45_ : fPPSRegionStart_56_) * 1e3) * a_y;
0193
0194   unsigned int line = in_trk->barcode();
0195   double px = -p * a_x;
0196   double py = p * a_y;
0197   double pz = std::copysign(sqrt(p * p - px * px - py * py), mom_cms.z());
0198   double e = sqrt(px * px + py * py + pz * pz + ProtonMassSQ);
0199   TLorentzVector p_out(px, py, pz, e);
0200   m_beamPart[line] = p_out;
0201   m_xAtTrPoint[line] = -b_x;
0202   m_yAtTrPoint[line] = b_y;
0203   return true;
0204 }