#include <OpticalFunctionsTransport.h>

Inheritance diagram for OpticalFunctionsTransport:

Public Member Functions
	OpticalFunctionsTransport (const edm::ParameterSet &ps, edm::ConsumesCollector iC)

void	process (const HepMC::GenEvent ev, const edm::EventSetup &es, CLHEP::HepRandomEngine engine) override

	~OpticalFunctionsTransport () override

Public Member Functions inherited from BaseProtonTransport
void	addPartToHepMC (const HepMC::GenEvent , HepMC::GenEvent )

void	ApplyBeamCorrection (HepMC::GenParticle *p)

void	ApplyBeamCorrection (TLorentzVector &p)

	BaseProtonTransport (const edm::ParameterSet &iConfig)

double	beamEnergy ()

double	beamMomentum ()

void	clear ()

std::vector< LHCTransportLink > &	getCorrespondenceMap ()

void	setBeamEnergy (double e)

virtual	~BaseProtonTransport ()

Private Member Functions
bool	transportProton (const HepMC::GenParticle *)

Private Attributes
const CTPPSBeamParameters *	beamParameters_

edm::ESGetToken < CTPPSBeamParameters, CTPPSBeamParametersRcd >	beamParametersToken_

const BeamSpotObjects *	beamspot_

edm::ESGetToken < BeamSpotObjects, BeamSpotObjectsRcd >	beamspotToken_

double	empiricalAperture45_a_int_

double	empiricalAperture45_a_slp_

double	empiricalAperture45_xi0_int_

double	empiricalAperture45_xi0_slp_

double	empiricalAperture56_a_int_

double	empiricalAperture56_a_slp_

double	empiricalAperture56_xi0_int_

double	empiricalAperture56_xi0_slp_

const LHCInfo *	lhcInfo_

edm::ESGetToken< LHCInfo, LHCInfoRcd >	lhcInfoToken_

unsigned int	optFunctionId45_

unsigned int	optFunctionId56_

const LHCInterpolatedOpticalFunctionsSetCollection *	opticalFunctions_

edm::ESGetToken < LHCInterpolatedOpticalFunctionsSetCollection, CTPPSInterpolatedOpticsRcd >	opticsToken_

bool	produceHitsRelativeToBeam_

bool	useEmpiricalApertures_

Additional Inherited Members
Protected Types inherited from BaseProtonTransport
enum	TransportMode { TransportMode::HECTOR, TransportMode::TOTEM, TransportMode::OPTICALFUNCTIONS }

Protected Attributes inherited from BaseProtonTransport
bool	bApplyZShift

std::string	beam1Filename_

std::string	beam2Filename_

double	beamEnergy_

double	beamMomentum_

CLHEP::HepRandomEngine *	engine_

double	etaCut_

double	fCrossingAngleX_45

double	fCrossingAngleX_56

double	fCrossingAngleY_45

double	fCrossingAngleY_56

double	fPPSRegionStart_45

double	fPPSRegionStart_56

std::map< unsigned int, TLorentzVector >	m_beamPart

std::vector< LHCTransportLink >	m_CorrespondenceMap

double	m_sig_E

double	m_sigmaSTX

double	m_sigmaSTY

double	m_sigmaSX

double	m_sigmaSY

std::map< unsigned int, double >	m_xAtTrPoint

std::map< unsigned int, double >	m_yAtTrPoint

TransportMode	MODE

double	momentumCut_

int	NEvent

bool	produceHitsRelativeToBeam_

bool	useBeamPositionFromLHCInfo_

bool	verbosity_

Detailed Description

Definition at line 23 of file OpticalFunctionsTransport.h.

Constructor & Destructor Documentation

OpticalFunctionsTransport::OpticalFunctionsTransport	(	const edm::ParameterSet &	ps,
		edm::ConsumesCollector	iC
	)

Definition at line 4 of file OpticalFunctionsTransport.cc.

References BaseProtonTransport::MODE, and BaseProtonTransport::OPTICALFUNCTIONS.

     : BaseProtonTransport(iConfig),
       lhcInfoToken_(iC.esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("lhcInfoLabel")))),
       beamParametersToken_(iC.esConsumes()),
       opticsToken_(iC.esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("opticsLabel")))),
       beamspotToken_(iC.esConsumes()),
       useEmpiricalApertures_(iConfig.getParameter<bool>("useEmpiricalApertures")),
       empiricalAperture45_xi0_int_(iConfig.getParameter<double>("empiricalAperture45_xi0_int")),
       empiricalAperture45_xi0_slp_(iConfig.getParameter<double>("empiricalAperture45_xi0_slp")),
       empiricalAperture45_a_int_(iConfig.getParameter<double>("empiricalAperture45_a_int")),
       empiricalAperture45_a_slp_(iConfig.getParameter<double>("empiricalAperture45_a_slp")),
       empiricalAperture56_xi0_int_(iConfig.getParameter<double>("empiricalAperture56_xi0_int")),
       empiricalAperture56_xi0_slp_(iConfig.getParameter<double>("empiricalAperture56_xi0_slp")),
       empiricalAperture56_a_int_(iConfig.getParameter<double>("empiricalAperture56_a_int")),
       empiricalAperture56_a_slp_(iConfig.getParameter<double>("empiricalAperture56_a_slp")) {
   MODE = TransportMode::OPTICALFUNCTIONS;
 }

OpticalFunctionsTransport::~OpticalFunctionsTransport ( )

inlineoverride

Definition at line 26 of file OpticalFunctionsTransport.h.

26 {};

Member Function Documentation

void OpticalFunctionsTransport::process	(	const HepMC::GenEvent *	ev,
		const edm::EventSetup &	es,
		CLHEP::HepRandomEngine *	engine
	)

overridevirtual

Implements BaseProtonTransport.

Definition at line 21 of file OpticalFunctionsTransport.cc.

References beamParameters_, beamParametersToken_, beamspot_, beamspotToken_, BaseProtonTransport::clear(), BaseProtonTransport::engine_, BaseProtonTransport::etaCut_, GenParticle::GenParticle, edm::EventSetup::getData(), lhcInfo_, lhcInfoToken_, geometryCSVtoXML::line, BaseProtonTransport::m_beamPart, BaseProtonTransport::momentumCut_, optFunctionId45_, optFunctionId56_, opticalFunctions_, opticsToken_, and transportProton().

                                                                        {
   this->clear();
 
   lhcInfo_ = &iSetup.getData(lhcInfoToken_);
   beamParameters_ = &iSetup.getData(beamParametersToken_);
   opticalFunctions_ = &iSetup.getData(opticsToken_);
   beamspot_ = &iSetup.getData(beamspotToken_);
 
   // Choose the optical function corresponding to the first station ono each side (it is in lhc ref. frame)
   optFunctionId45_ = 0;
   optFunctionId56_ = 0;
   for (const auto& ofp : (*opticalFunctions_)) {
     if (ofp.second.getScoringPlaneZ() < 0) {
       if (optFunctionId45_ == 0)
         optFunctionId45_ = ofp.first;
       if (opticalFunctions_->at(optFunctionId45_).getScoringPlaneZ() < ofp.second.getScoringPlaneZ())
         optFunctionId45_ = ofp.first;
     }
     if (ofp.second.getScoringPlaneZ() > 0) {
       if (optFunctionId56_ == 0)
         optFunctionId56_ = ofp.first;
       if (opticalFunctions_->at(optFunctionId56_).getScoringPlaneZ() > ofp.second.getScoringPlaneZ())
         optFunctionId56_ = ofp.first;
     }
   }
   //
   engine_ = _engine;  // the engine needs to be updated for each event
 
   for (HepMC::GenEvent::particle_const_iterator eventParticle = evt->particles_begin();
        eventParticle != evt->particles_end();
        ++eventParticle) {
     if (!((*eventParticle)->status() == 1 && (*eventParticle)->pdg_id() == 2212))
       continue;
 
     if (!(fabs((*eventParticle)->momentum().eta()) > etaCut_ && fabs((*eventParticle)->momentum().pz()) > momentumCut_))
       continue;  // discard protons outside kinematic acceptance
 
     unsigned int line = (*eventParticle)->barcode();
     HepMC::GenParticle* gpart = (*eventParticle);
     if (gpart->pdg_id() != 2212)
       continue;  // only transport stable protons
     if (gpart->status() != 1)
       continue;
     if (m_beamPart.find(line) != m_beamPart.end())  // assures this protons has not been already propagated
       continue;
 
     transportProton(gpart);
   }
 }

bool OpticalFunctionsTransport::transportProton ( const HepMC::GenParticle * in_trk )

private

Definition at line 73 of file OpticalFunctionsTransport.cc.

References funct::abs(), CTPPSDetId::arm(), BaseProtonTransport::beamMomentum(), beamParameters_, beamspot_, cm_to_mm, LHCInfo::crossingAngle(), alignCSCRings::e, empiricalAperture45_a_int_, empiricalAperture45_a_slp_, empiricalAperture45_xi0_int_, empiricalAperture45_xi0_slp_, empiricalAperture56_a_int_, empiricalAperture56_a_slp_, empiricalAperture56_xi0_int_, empiricalAperture56_xi0_slp_, BaseProtonTransport::fPPSRegionStart_45, BaseProtonTransport::fPPSRegionStart_56, CTPPSBeamParameters::getBeamMom45(), CTPPSBeamParameters::getBeamMom56(), CTPPSBeamParameters::getHalfXangleX45(), CTPPSBeamParameters::getHalfXangleX56(), CTPPSBeamParameters::getVtxOffsetX45(), CTPPSBeamParameters::getVtxOffsetY45(), CTPPSBeamParameters::getVtxOffsetZ45(), BeamSpotObjects::GetX(), BeamSpotObjects::GetY(), BeamSpotObjects::GetZ(), lhcInfo_, geometryCSVtoXML::line, LogDebug, BaseProtonTransport::m_beamPart, BaseProtonTransport::m_xAtTrPoint, BaseProtonTransport::m_yAtTrPoint, optFunctionId45_, optFunctionId56_, opticalFunctions_, AlCaHLTBitMon_ParallelJobs::p, produceHitsRelativeToBeam_, ProtonMassSQ, CTPPSDetId::rp(), mathSSE::sqrt(), CTPPSDetId::station(), funct::tan(), LHCInterpolatedOpticalFunctionsSet::Kinematics::th_x, LHCInterpolatedOpticalFunctionsSet::Kinematics::th_y, urad, BaseProtonTransport::useBeamPositionFromLHCInfo_, useEmpiricalApertures_, BaseProtonTransport::verbosity_, LHCInterpolatedOpticalFunctionsSet::Kinematics::x, and LHCInterpolatedOpticalFunctionsSet::Kinematics::y.

Referenced by process().

                                                                               {
   const HepMC::FourVector& vtx_cms = in_trk->production_vertex()->position();  // in mm
   const HepMC::FourVector& mom_cms = in_trk->momentum();
 
   // transformation to LHC/TOTEM convention
   HepMC::FourVector vtx_lhc(-vtx_cms.x(), vtx_cms.y(), -vtx_cms.z(), vtx_cms.t());
   HepMC::FourVector mom_lhc(-mom_cms.x(), mom_cms.y(), -mom_cms.z(), mom_cms.t());
 
   // determine the LHC arm and related parameters
   double urad = 1.e-6;
   double beamMomentum = 0.;
   double xangle = 0.;
   double empiricalAperture_xi0_int, empiricalAperture_xi0_slp;
   double empiricalAperture_a_int, empiricalAperture_a_slp;
   unsigned int optFunctionId_;
   // get the beam position at the IP in mm and in the LHC ref. frame
   double vtxXoffset_;
   double vtxYoffset_;
   double vtxZoffset_;
   if (useBeamPositionFromLHCInfo_) {
     vtxXoffset_ = -beamParameters_->getVtxOffsetX45() * cm_to_mm;
     vtxYoffset_ = beamParameters_->getVtxOffsetY45() * cm_to_mm;
     vtxZoffset_ = -beamParameters_->getVtxOffsetZ45() * cm_to_mm;
   } else {
     vtxXoffset_ = -beamspot_->GetX() * cm_to_mm;
     vtxYoffset_ = beamspot_->GetY() * cm_to_mm;
     vtxZoffset_ = -beamspot_->GetZ() * cm_to_mm;
   }
   vtxZoffset_ *= 1.0;   // just to avoid compilation error, should be used in the future
                         //
   if (mom_lhc.z() < 0)  // sector 45
   {
     optFunctionId_ = optFunctionId45_;
     beamMomentum = beamParameters_->getBeamMom45();
     xangle = beamParameters_->getHalfXangleX45();
     empiricalAperture_xi0_int = empiricalAperture45_xi0_int_;
     empiricalAperture_xi0_slp = empiricalAperture45_xi0_slp_;
     empiricalAperture_a_int = empiricalAperture45_a_int_;
     empiricalAperture_a_slp = empiricalAperture45_a_slp_;
   } else {  // sector 56
     optFunctionId_ = optFunctionId56_;
     beamMomentum = beamParameters_->getBeamMom56();
     xangle = beamParameters_->getHalfXangleX56();
     empiricalAperture_xi0_int = empiricalAperture56_xi0_int_;
     empiricalAperture_xi0_slp = empiricalAperture56_xi0_slp_;
     empiricalAperture_a_int = empiricalAperture56_a_int_;
     empiricalAperture_a_slp = empiricalAperture56_a_slp_;
   }
   if (xangle > 1.0)
     xangle *= urad;
   // calculate kinematics for optics parametrisation, avoid the aproximation for small angles xangle -> tan(xangle)
   const double p = mom_lhc.rho();
   const double xi = 1. - p / beamMomentum;
   const double th_x_phys = mom_lhc.x() / abs(mom_lhc.z()) - tan(xangle);  //"-" in the LHC ref. frame
   const double th_y_phys = mom_lhc.y() / abs(mom_lhc.z());
   const double vtx_lhc_eff_x = vtx_lhc.x() - vtx_lhc.z() * (mom_lhc.x() / mom_lhc.z() + tan(xangle)) - (vtxXoffset_);
   const double vtx_lhc_eff_y = vtx_lhc.y() - vtx_lhc.z() * (mom_lhc.y() / mom_lhc.z()) - (vtxYoffset_);
 
   if (verbosity_) {
     LogDebug("OpticalFunctionsTransport")
         << "simu: xi = " << xi << ", th_x_phys = " << th_x_phys << ", th_y_phys = " << th_y_phys
         << ", vtx_lhc_eff_x = " << vtx_lhc_eff_x << ", vtx_lhc_eff_y = " << vtx_lhc_eff_y << std::endl;
   }
 
   // check empirical aperture
   if (useEmpiricalApertures_) {
     const auto& xangle =
         (lhcInfo_->crossingAngle() > 1.0) ? lhcInfo_->crossingAngle() * urad : lhcInfo_->crossingAngle();
     const double xi_th = (empiricalAperture_xi0_int + xangle * empiricalAperture_xi0_slp) +
                          (empiricalAperture_a_int + xangle * empiricalAperture_a_slp) * th_x_phys;
 
     if (xi > xi_th) {
       if (verbosity_) {
         LogDebug("OpticalFunctionsTransport") << "stop because of empirical appertures";
       }
       return false;
     }
   }
 
   // transport the proton into  pot/scoring plane
   auto ofp = opticalFunctions_->at(optFunctionId_);
   CTPPSDetId rpId(optFunctionId_);
   const unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
 
   if (verbosity_)
     LogDebug("OpticalFunctionsTransport") << "  RP " << rpDecId << std::endl;
 
   // transport proton
   LHCInterpolatedOpticalFunctionsSet::Kinematics k_in = {
       vtx_lhc_eff_x * 1E-1, th_x_phys, vtx_lhc_eff_y * 1E-1, th_y_phys, xi};  // conversions: mm -> cm
 
   LHCInterpolatedOpticalFunctionsSet::Kinematics k_out;
   ofp.transport(k_in, k_out, true);
 
   // Original code uses mm, but CMS uses cm, so keep it in cm
   double b_x = k_out.x * 1E1, b_y = k_out.y * 1E1;  // conversions: cm -> mm
   double a_x = k_out.th_x, a_y = k_out.th_y;
 
   // if needed, subtract beam position and angle
   if (produceHitsRelativeToBeam_) {
     // determine beam position
     LHCInterpolatedOpticalFunctionsSet::Kinematics k_be_in = {0., -tan(xangle), 0., 0., 0.};
     LHCInterpolatedOpticalFunctionsSet::Kinematics k_be_out;
     ofp.transport(k_be_in, k_be_out, true);
 
     a_x -= k_be_out.th_x;
     a_y -= k_be_out.th_y;
     b_x -= k_be_out.x * 1E1;
     b_y -= k_be_out.y * 1E1;  // conversions: cm -> mm
   }
 
   const double z_scoringPlane = ofp.getScoringPlaneZ() * 1E1;  // conversion: cm --> mm
 
   if (verbosity_) {
     LogDebug("OpticalFunctionsTransport")
         << "    proton transported: a_x = " << a_x << " rad, a_y = " << a_y << " rad, b_x = " << b_x
         << " mm, b_y = " << b_y << " mm, z = " << z_scoringPlane << " mm" << std::endl;
   }
   //
   // Project the track back to the starting of PPS region
   b_x -= (abs(z_scoringPlane) - (double)((z_scoringPlane < 0) ? fPPSRegionStart_45 : fPPSRegionStart_56) * 1e3) *
          a_x;  // z_scoringPlane is in mm
   b_y -= (abs(z_scoringPlane) - (double)((z_scoringPlane < 0) ? fPPSRegionStart_45 : fPPSRegionStart_56) * 1e3) * a_y;
 
   unsigned int line = in_trk->barcode();
   double px = -p * a_x;
   double py = p * a_y;
   double pz = std::copysign(sqrt(p * p - px * px - py * py), mom_cms.z());
   double e = sqrt(px * px + py * py + pz * pz + ProtonMassSQ);
   TLorentzVector p_out(px, py, pz, e);
   m_beamPart[line] = p_out;
   m_xAtTrPoint[line] = -b_x;
   m_yAtTrPoint[line] = b_y;
   return true;
 }