Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

PPS transport update #36589

Merged
merged 9 commits into from
Jan 4, 2022
Merged

PPS transport update #36589

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
a4a5fbe
Make PPS transport thread safe
civanch Dec 24, 2021
d154254
code format
civanch Dec 24, 2021
5184385
fixed static analysis warnings
civanch Dec 26, 2021
5134ac3
implement CMS code convention, removed unused variables and headers, …
civanch Dec 28, 2021
836573f
code format
civanch Dec 28, 2021
43afc3b
removed commented lines
civanch Dec 28, 2021
6528f03
move initialisation of common parameters to the base class
civanch Dec 29, 2021
369828a
reorganise initialisation
civanch Jan 3, 2022
d0f9ffe
code format
civanch Jan 3, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
64 changes: 40 additions & 24 deletions SimTransport/PPSProtonTransport/interface/BaseProtonTransport.h
View file
Open in desktop
Expand Up @@ -17,7 +17,7 @@ namespace CLHEP {
class BaseProtonTransport {
public:
BaseProtonTransport(const edm::ParameterSet& iConfig);
virtual ~BaseProtonTransport() { this->clear(); };
virtual ~BaseProtonTransport();

std::vector<LHCTransportLink>& getCorrespondenceMap() { return m_CorrespondenceMap; }
virtual void process(const HepMC::GenEvent* ev, const edm::EventSetup& es, CLHEP::HepRandomEngine* engine) = 0;
Expand All @@ -29,49 +29,65 @@ class BaseProtonTransport {
void ApplyBeamCorrection(HepMC::GenParticle* p);
void ApplyBeamCorrection(TLorentzVector& p);

void setBeamEnergy(double e) {
beamEnergy_ = e;
beamMomentum_ = sqrt(beamEnergy_ * beamEnergy_ - ProtonMassSQ);
}

double beamEnergy() { return beamEnergy_; };
double beamMomentum() { return beamMomentum_; };

protected:
void setBeamFileNames(const std::string& nam1, const std::string& nam2) {
beam1Filename_ = nam1;
beam2Filename_ = nam2;
};

void setBeamParameters(double stx, double sty, double sx, double sy, double se) {
m_sigmaSTX = stx;
m_sigmaSTY = sty;
m_sigmaSX = sx;
m_sigmaSY = sy;
m_sig_E = se;
};

void setCrossingAngles(double cx45, double cx56, double cy45, double cy56) {
fCrossingAngleX_45_ = cx45;
fCrossingAngleX_56_ = cx56;
fCrossingAngleY_45_ = cy45;
fCrossingAngleY_56_ = cy56;
};

enum class TransportMode { HECTOR, TOTEM, OPTICALFUNCTIONS };
TransportMode MODE;

int NEvent;
CLHEP::HepRandomEngine* engine_;
CLHEP::HepRandomEngine* engine_{nullptr};
std::vector<LHCTransportLink> m_CorrespondenceMap;
std::map<unsigned int, TLorentzVector> m_beamPart;
std::map<unsigned int, double> m_xAtTrPoint;
std::map<unsigned int, double> m_yAtTrPoint;

bool verbosity_;
bool bApplyZShift;
bool bApplyZShift_;
bool useBeamPositionFromLHCInfo_;
bool produceHitsRelativeToBeam_;

std::string beam1Filename_;
std::string beam2Filename_;
std::string beam1Filename_{""};
std::string beam2Filename_{""};

double fPPSRegionStart_45;
double fPPSRegionStart_56;
double fCrossingAngleX_45;
double fCrossingAngleX_56;
double fCrossingAngleY_45;
double fCrossingAngleY_56;
double fPPSRegionStart_45_;
double fPPSRegionStart_56_;
double fCrossingAngleX_45_{0.0};
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These fCrossingAngle* are needed inside BaseProtonTransport::ApplyBeamCorrection (that method is only used in the derived class TotemTransport) and in HectorTransport::transportProton

I think you should either initialize these class members in the constructor of BaseProtonTransport (via iConfig.getParameter etc.) or move the ApplyBeamCorrection method in the derived class TotemTransport, removing it from the base class, and let these fCrossingAngle* as class members of the two derived classes TotemTransport and HectorTransport

double fCrossingAngleX_56_{0.0};
double fCrossingAngleY_45_{0.0};
double fCrossingAngleY_56_{0.0};

double beamMomentum_;
double beamEnergy_;
double etaCut_;
double momentumCut_;

double m_sigmaSTX;
double m_sigmaSTY;
double m_sigmaSX;
double m_sigmaSY;
double m_sig_E;
private:
double beamMomentum_;
double beamEnergy_;

double m_sigmaSTX{0.0};
civanch marked this conversation as resolved.
Show resolved Hide resolved
civanch marked this conversation as resolved.
Show resolved Hide resolved
double m_sigmaSTY{0.0};
double m_sigmaSX{0.0};
double m_sigmaSY{0.0};
double m_sig_E{0.0};
};
#endif
15 changes: 4 additions & 11 deletions SimTransport/PPSProtonTransport/interface/HectorTransport.h
View file
Open in desktop
Expand Up @@ -7,7 +7,7 @@
#include "CondFormats/PPSObjects/interface/CTPPSBeamParameters.h"

#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Utilities/interface/EDGetToken.h"
#include "FWCore/Framework/interface/ConsumesCollector.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
Expand Down Expand Up @@ -44,20 +44,13 @@ class HectorTransport : public BaseProtonTransport {
void process(const HepMC::GenEvent* ev, const edm::EventSetup& es, CLHEP::HepRandomEngine* engine) override;

private:
bool m_verbosity;

static constexpr double fPPSBeamLineLength_ = 250.; // default beam line length

//!propagate the particles through a beamline to PPS
bool transportProton(const HepMC::GenParticle*);

// New function to calculate the LorentzBoost

// function to calculate the LorentzBoost
bool setBeamLine();
// Defaults

double m_fEtacut;
double m_fMomentumMin;

// PPSHector
std::unique_ptr<H_BeamLine> m_beamline45;
Expand All @@ -66,7 +59,7 @@ class HectorTransport : public BaseProtonTransport {
edm::ESGetToken<CTPPSBeamParameters, CTPPSBeamParametersRcd> beamParametersToken_;
edm::ESGetToken<BeamSpotObjects, BeamSpotObjectsRcd> beamspotToken_;

const CTPPSBeamParameters* beamParameters_;
const BeamSpotObjects* beamspot_;
const CTPPSBeamParameters* beamParameters_{nullptr};
const BeamSpotObjects* beamspot_{nullptr};
};
#endif
20 changes: 7 additions & 13 deletions SimTransport/PPSProtonTransport/interface/OpticalFunctionsTransport.h
View file
Open in desktop
Expand Up @@ -16,14 +16,10 @@
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include "SimTransport/PPSProtonTransport/interface/BaseProtonTransport.h"

#include <array>
#include <unordered_map>
#include <cmath>

class OpticalFunctionsTransport : public BaseProtonTransport {
public:
OpticalFunctionsTransport(const edm::ParameterSet& ps, edm::ConsumesCollector iC);
~OpticalFunctionsTransport() override{};
~OpticalFunctionsTransport() override;

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

Expand All @@ -35,20 +31,18 @@ class OpticalFunctionsTransport : public BaseProtonTransport {
edm::ESGetToken<LHCInterpolatedOpticalFunctionsSetCollection, CTPPSInterpolatedOpticsRcd> opticsToken_;
edm::ESGetToken<BeamSpotObjects, BeamSpotObjectsRcd> beamspotToken_;

const LHCInfo* lhcInfo_;
const CTPPSBeamParameters* beamParameters_;
const LHCInterpolatedOpticalFunctionsSetCollection* opticalFunctions_;
const BeamSpotObjects* beamspot_;
const LHCInfo* lhcInfo_{nullptr};
const CTPPSBeamParameters* beamParameters_{nullptr};
const LHCInterpolatedOpticalFunctionsSetCollection* opticalFunctions_{nullptr};
const BeamSpotObjects* beamspot_{nullptr};

unsigned int optFunctionId45_;
unsigned int optFunctionId56_;
unsigned int optFunctionId45_{0};
unsigned int optFunctionId56_{0};

bool useEmpiricalApertures_;
double empiricalAperture45_xi0_int_, empiricalAperture45_xi0_slp_, empiricalAperture45_a_int_,
empiricalAperture45_a_slp_;
double empiricalAperture56_xi0_int_, empiricalAperture56_xi0_slp_, empiricalAperture56_a_int_,
empiricalAperture56_a_slp_;

bool produceHitsRelativeToBeam_;
};
#endif
20 changes: 7 additions & 13 deletions SimTransport/PPSProtonTransport/interface/TotemTransport.h
View file
Open in desktop
Expand Up @@ -11,10 +11,6 @@
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"

#include <unordered_map>
#include <array>

namespace CLHEP {
class HepRandomEngine;
Expand All @@ -23,22 +19,20 @@ namespace CLHEP {
class TotemTransport : public BaseProtonTransport {
public:
TotemTransport(const edm::ParameterSet& ps);
~TotemTransport() override{};
~TotemTransport() override;
// look for scattered protons, propagates them, add them to the event
void process(const HepMC::GenEvent* ev, const edm::EventSetup& es, CLHEP::HepRandomEngine* engine) override;

private:
bool transportProton(const HepMC::GenParticle*);
bool transportProton(HepMC::GenParticle*);
LHCOpticsApproximator* ReadParameterization(const std::string&, const std::string&);

LHCOpticsApproximator* m_aprox_ip_150_r = nullptr;
LHCOpticsApproximator* m_aprox_ip_150_l = nullptr;
LHCOpticsApproximator* m_aprox_ip_150_r_;
LHCOpticsApproximator* m_aprox_ip_150_l_;

std::string m_model_ip_150_r_name;
std::string m_model_ip_150_l_name;
std::string m_model_ip_150_r_name_;
std::string m_model_ip_150_l_name_;

double m_beampipe_aperture_radius;
double m_fEtacut;
double m_fMomentumMin;
double m_beampipe_aperture_radius_;
};
#endif
46 changes: 23 additions & 23 deletions SimTransport/PPSProtonTransport/src/BaseProtonTransport.cc
View file
Open in desktop
Expand Up @@ -3,20 +3,22 @@
#include "Utilities/PPS/interface/PPSUnitConversion.h"
#include <CLHEP/Random/RandGauss.h>
#include <CLHEP/Units/GlobalSystemOfUnits.h>
#include <cctype>

BaseProtonTransport::BaseProtonTransport(const edm::ParameterSet& iConfig)
: verbosity_(iConfig.getParameter<bool>("Verbosity")),
civanch marked this conversation as resolved.
Show resolved Hide resolved
bApplyZShift(iConfig.getParameter<bool>("ApplyZShift")),
bApplyZShift_(iConfig.getParameter<bool>("ApplyZShift")),
useBeamPositionFromLHCInfo_(iConfig.getParameter<bool>("useBeamPositionFromLHCInfo")),
produceHitsRelativeToBeam_(iConfig.getParameter<bool>("produceHitsRelativeToBeam")),
fPPSRegionStart_45(iConfig.getParameter<double>("PPSRegionStart_45")),
fPPSRegionStart_56(iConfig.getParameter<double>("PPSRegionStart_56")),
beamEnergy_(iConfig.getParameter<double>("BeamEnergy")),
fPPSRegionStart_45_(iConfig.getParameter<double>("PPSRegionStart_45")),
fPPSRegionStart_56_(iConfig.getParameter<double>("PPSRegionStart_56")),
etaCut_(iConfig.getParameter<double>("EtaCut")),
momentumCut_(iConfig.getParameter<double>("MomentumCut")) {
setBeamEnergy(beamEnergy_);
momentumCut_(iConfig.getParameter<double>("MomentumCut")),
beamEnergy_(iConfig.getParameter<double>("BeamEnergy")) {
beamMomentum_ = std::sqrt(beamEnergy_ * beamEnergy_ - ProtonMassSQ);
}

BaseProtonTransport::~BaseProtonTransport() { clear(); }

void BaseProtonTransport::ApplyBeamCorrection(HepMC::GenParticle* p) {
TLorentzVector p_out;
p_out.SetPx(p->momentum().px());
Expand All @@ -26,37 +28,35 @@ void BaseProtonTransport::ApplyBeamCorrection(HepMC::GenParticle* p) {
ApplyBeamCorrection(p_out);
p->set_momentum(HepMC::FourVector(p_out.Px(), p_out.Py(), p_out.Pz(), p_out.E()));
}

void BaseProtonTransport::ApplyBeamCorrection(TLorentzVector& p_out) {
double theta = p_out.Theta();
double thetax = atan(p_out.Px() / fabs(p_out.Pz()));
double thetay = atan(p_out.Py() / fabs(p_out.Pz()));
double energy = p_out.E();
double urad = 1e-6;

int direction = (p_out.Pz() > 0) ? 1 : -1;

if (p_out.Pz() < 0)
theta = TMath::Pi() - theta;

if (MODE == TransportMode::TOTEM)
thetax += (p_out.Pz() > 0) ? fCrossingAngleX_45 * urad : fCrossingAngleX_56 * urad;
thetax += (p_out.Pz() > 0) ? fCrossingAngleX_45_ * urad : fCrossingAngleX_56_ * urad;

double dtheta_x = (double)CLHEP::RandGauss::shoot(engine_, 0., m_sigmaSTX);
double dtheta_y = (double)CLHEP::RandGauss::shoot(engine_, 0., m_sigmaSTY);
double denergy = (double)CLHEP::RandGauss::shoot(engine_, 0., m_sig_E);
double dtheta_x = (m_sigmaSTX <= 0.0) ? 0.0 : CLHEP::RandGauss::shoot(engine_, 0., m_sigmaSTX);
double dtheta_y = (m_sigmaSTY <= 0.0) ? 0.0 : CLHEP::RandGauss::shoot(engine_, 0., m_sigmaSTY);
double denergy = (m_sig_E <= 0.0) ? 0.0 : CLHEP::RandGauss::shoot(engine_, 0., m_sig_E);

double s_theta = sqrt(pow(thetax + dtheta_x * urad, 2) + pow(thetay + dtheta_y * urad, 2));
double s_phi = atan2(thetay + dtheta_y * urad, thetax + dtheta_x * urad);
double s_theta = std::sqrt(pow(thetax + dtheta_x * urad, 2) + std::pow(thetay + dtheta_y * urad, 2));
double s_phi = std::atan2(thetay + dtheta_y * urad, thetax + dtheta_x * urad);
energy += denergy;
double p = sqrt(pow(energy, 2) - ProtonMassSQ);
double p = std::sqrt(std::pow(energy, 2) - ProtonMassSQ);
double sint = std::sin(s_theta);

p_out.SetPx((double)p * sin(s_theta) * cos(s_phi));
p_out.SetPy((double)p * sin(s_theta) * sin(s_phi));
p_out.SetPz((double)p * (cos(s_theta)) * direction);
p_out.SetPx(p * sint * std::cos(s_phi));
p_out.SetPy(p * sint * std::sin(s_phi));
p_out.SetPz(p * std::cos(s_theta) * direction);
p_out.SetE(energy);
}

void BaseProtonTransport::addPartToHepMC(const HepMC::GenEvent* in_evt, HepMC::GenEvent* evt) {
NEvent++;
m_CorrespondenceMap.clear();

int direction = 0;
Expand All @@ -71,7 +71,7 @@ void BaseProtonTransport::addPartToHepMC(const HepMC::GenEvent* in_evt, HepMC::G
direction = (gpart->momentum().pz() > 0) ? 1 : -1;

// Totem uses negative Z for sector 56 while Hector uses always positive distance
double ddd = (direction > 0) ? fPPSRegionStart_45 : fabs(fPPSRegionStart_56);
double ddd = (direction > 0) ? fPPSRegionStart_45_ : fabs(fPPSRegionStart_56_);

double time = (ddd * meter - gpart->production_vertex()->position().z() * mm); // mm

Expand Down