File indexing completed on 2024-05-10 02:21:29
0001 #include "SimTransport/PPSProtonTransport/interface/TotemTransport.h"
0002 #include "FWCore/ParameterSet/interface/FileInPath.h"
0003 #include <CLHEP/Units/SystemOfUnits.h>
0004 #include "TLorentzVector.h"
0005 #include "TFile.h"
0006
0007 #include <cmath>
0008
0009 TotemTransport::TotemTransport(const edm::ParameterSet& iConfig)
0010 : BaseProtonTransport(iConfig),
0011 m_model_ip_150_r_name_(iConfig.getParameter<std::string>("Model_IP_150_R_Name")),
0012 m_model_ip_150_l_name_(iConfig.getParameter<std::string>("Model_IP_150_L_Name")),
0013 m_beampipe_aperture_radius_(iConfig.getParameter<double>("BeampipeApertureRadius")) {
0014 MODE = TransportMode::TOTEM;
0015 std::string s1 = iConfig.getParameter<std::string>("Beam1Filename");
0016 std::string s2 = iConfig.getParameter<std::string>("Beam2Filename");
0017 setBeamFileNames(s1, s2);
0018 double cax45 = iConfig.getParameter<double>("halfCrossingAngleXSector45");
0019 double cax56 = iConfig.getParameter<double>("halfCrossingAngleXSector56");
0020 setCrossingAngles(cax45, cax56, 0.0, 0.0);
0021 double stx = iConfig.getParameter<double>("BeamDivergenceX");
0022 double sty = iConfig.getParameter<double>("BeamDivergenceY");
0023 double sx = iConfig.getParameter<double>("BeamSigmaX");
0024 double sy = iConfig.getParameter<double>("BeamSigmaY");
0025 double se = iConfig.getParameter<double>("BeamEnergyDispersion");
0026 setBeamParameters(stx, sty, sx, sy, se);
0027
0028 if (fPPSRegionStart_56_ > 0)
0029 fPPSRegionStart_56_ *= -1;
0030
0031 edm::LogVerbatim("TotemTransport")
0032 << "=============================================================================\n"
0033 << " Bulding LHC Proton transporter based on TOTEM model\n"
0034 << "=============================================================================\n";
0035
0036 m_aprox_ip_150_r_ = ReadParameterization(m_model_ip_150_r_name_, beam1Filename_);
0037 m_aprox_ip_150_l_ = ReadParameterization(m_model_ip_150_l_name_, beam2Filename_);
0038
0039 if (m_aprox_ip_150_r_ == nullptr || m_aprox_ip_150_l_ == nullptr) {
0040 edm::LogError("TotemTransport") << "Parameterisation " << m_model_ip_150_r_name_ << " or " << m_model_ip_150_l_name_
0041 << " missing in file. Cannot proceed. ";
0042 throw edm::Exception(edm::errors::Configuration) << "TotemTransport is not properly initialized";
0043 }
0044 edm::LogVerbatim("TotemTransport") << "Parameterizations read from file, pointers:" << m_aprox_ip_150_r_ << " "
0045 << m_aprox_ip_150_l_ << " ";
0046 }
0047 TotemTransport::~TotemTransport() {}
0048
0049
0050
0051
0052 void TotemTransport::process(const HepMC::GenEvent* ievt,
0053 const edm::EventSetup& iSetup,
0054 CLHEP::HepRandomEngine* engine) {
0055 clear();
0056 engine_ = engine;
0057
0058 HepMC::GenEvent evt(*ievt);
0059
0060 for (HepMC::GenEvent::particle_const_iterator eventParticle = evt.particles_begin();
0061 eventParticle != evt.particles_end();
0062 ++eventParticle) {
0063 if (!((*eventParticle)->status() == 1 && (*eventParticle)->pdg_id() == 2212))
0064 continue;
0065
0066 if (!(fabs((*eventParticle)->momentum().eta()) > etaCut_ && fabs((*eventParticle)->momentum().pz()) > momentumCut_))
0067 continue;
0068
0069 unsigned int line = (*eventParticle)->barcode();
0070 HepMC::GenParticle* gpart = (*eventParticle);
0071 if (gpart->pdg_id() != 2212)
0072 continue;
0073 if (gpart->status() != 1)
0074 continue;
0075 if (m_beamPart.find(line) != m_beamPart.end())
0076 continue;
0077
0078 transportProton(gpart);
0079 }
0080 }
0081
0082
0083
0084
0085
0086 bool TotemTransport::transportProton(HepMC::GenParticle* in_trk) {
0087
0088 edm::LogVerbatim("TotemTransport") << "Starting proton transport using TOTEM method\n";
0089
0090 ApplyBeamCorrection(in_trk);
0091
0092 const HepMC::GenVertex* in_pos = in_trk->production_vertex();
0093 const HepMC::FourVector in_mom = in_trk->momentum();
0094
0095
0096
0097 double in_position[3] = {in_pos->position().x(), in_pos->position().y(), in_pos->position().z()};
0098
0099 double crossingAngleX = (in_mom.z() > 0) ? fCrossingAngleX_45_ : fCrossingAngleX_56_;
0100
0101
0102 in_position[0] =
0103 in_position[0] - in_position[2] * (in_mom.x() / in_mom.z() - crossingAngleX * urad);
0104 in_position[1] = in_position[1] - in_position[2] * (in_mom.y() / (in_mom.z()));
0105 in_position[2] = 0.;
0106 double in_momentum[3] = {in_mom.x(), in_mom.y(), in_mom.z()};
0107 double out_position[3];
0108 double out_momentum[3];
0109 edm::LogVerbatim("TotemTransport") << "before transport ->"
0110 << " position: " << in_position[0] << ", " << in_position[1] << ", "
0111 << in_position[2] << " momentum: " << in_momentum[0] << ", " << in_momentum[1]
0112 << ", " << in_momentum[2];
0113
0114 LHCOpticsApproximator* approximator = nullptr;
0115 double zin;
0116 double zout;
0117 if (in_mom.z() > 0) {
0118 approximator = m_aprox_ip_150_l_;
0119 zin = 0.0;
0120 zout = fPPSRegionStart_45_;
0121 } else {
0122 approximator = m_aprox_ip_150_r_;
0123 zin = 0.0;
0124 zout = fPPSRegionStart_56_;
0125 }
0126
0127 bool invert_beam_coord_system =
0128 true;
0129
0130 bool tracked = approximator->Transport_m_GeV(
0131 in_position, in_momentum, out_position, out_momentum, invert_beam_coord_system, zout - zin);
0132
0133 if (!tracked)
0134 return false;
0135
0136 edm::LogVerbatim("TotemTransport") << "after transport -> "
0137 << "position: " << out_position[0] << ", " << out_position[1] << ", "
0138 << out_position[2] << "momentum: " << out_momentum[0] << ", " << out_momentum[1]
0139 << ", " << out_momentum[2];
0140
0141 if (out_position[0] * out_position[0] + out_position[1] * out_position[1] >
0142 m_beampipe_aperture_radius_ * m_beampipe_aperture_radius_) {
0143 edm::LogVerbatim("TotemTransport") << "Proton ouside beampipe\n"
0144 << "===== END Transport "
0145 << "====================";
0146 return false;
0147 }
0148
0149 using CLHEP::meter;
0150 TVector3 out_pos(out_position[0] * meter, out_position[1] * meter, out_position[2] * meter);
0151 TVector3 out_mom(out_momentum[0], out_momentum[1], out_momentum[2]);
0152
0153 if (verbosity_) {
0154 LogDebug("TotemTransport") << "output -> position: ";
0155 out_pos.Print();
0156 LogDebug("TotemTransport") << " momentum: ";
0157 out_mom.Print();
0158 }
0159
0160 double px = -out_momentum[0];
0161 double py = out_momentum[1];
0162 double pz =
0163 out_momentum[2];
0164 double e = sqrt(px * px + py * py + pz * pz + ProtonMassSQ);
0165 TLorentzVector p_out(px, py, pz, e);
0166 double x1_ctpps = -out_position[0] * meter;
0167 double y1_ctpps = out_position[1] * meter;
0168
0169 unsigned int line = in_trk->barcode();
0170
0171 if (verbosity_)
0172 LogDebug("TotemTransport") << "TotemTransport:filterPPS: barcode = " << line << " x= " << x1_ctpps
0173 << " y= " << y1_ctpps;
0174
0175 m_beamPart[line] = p_out;
0176 m_xAtTrPoint[line] = x1_ctpps;
0177 m_yAtTrPoint[line] = y1_ctpps;
0178 return true;
0179 }
0180
0181 LHCOpticsApproximator* TotemTransport::ReadParameterization(const std::string& m_model_name,
0182 const std::string& rootfile) {
0183 edm::FileInPath fileName(rootfile.c_str());
0184 TFile* f = TFile::Open(fileName.fullPath().c_str(), "read");
0185 if (!f) {
0186 edm::LogError("TotemTransport") << "File " << fileName << " not found. Exiting.";
0187 return nullptr;
0188 }
0189 edm::LogVerbatim("TotemTransport") << "Root file opened, pointer:" << f;
0190
0191
0192 LHCOpticsApproximator* aprox = (LHCOpticsApproximator*)f->Get(m_model_name.c_str());
0193 f->Close();
0194 return aprox;
0195 }