PPSProtonTransport/src/TotemTransport.cc

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;  // make sure sector 56 has negative position, as TOTEM convention
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 // this method is the same for all propagator, but since transportProton is different for each derived class
0050 // it needes to be overriden
0051 //
0052 void TotemTransport::process(const HepMC::GenEvent* ievt,
0053                              const edm::EventSetup& iSetup,
0054                              CLHEP::HepRandomEngine* engine) {
0055   clear();
0056   engine_ = engine;  // the engine needs to be updated for each event
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;  // discard protons outside kinematic acceptance
0068
0069     unsigned int line = (*eventParticle)->barcode();
0070     HepMC::GenParticle* gpart = (*eventParticle);
0071     if (gpart->pdg_id() != 2212)
0072       continue;  // only transport stable protons
0073     if (gpart->status() != 1)
0074       continue;
0075     if (m_beamPart.find(line) != m_beamPart.end())  // assures this protons has not been already propagated
0076       continue;
0077
0078     transportProton(gpart);
0079   }
0080 }
0081 //
0082 //
0083 // here comes the real thing
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   // ATTENTION: HepMC uses mm, vertex config of CMS uses cm and SimTransport uses mm
0096   //
0097   double in_position[3] = {in_pos->position().x(), in_pos->position().y(), in_pos->position().z()};  //in LHC ref. frame
0098
0099   double crossingAngleX = (in_mom.z() > 0) ? fCrossingAngleX_45_ : fCrossingAngleX_56_;
0100
0101   // Move the position to z=0. Do it in the CMS ref frame. Totem parameterization does the rotation internatlly
0102   in_position[0] =
0103       in_position[0] - in_position[2] * (in_mom.x() / in_mom.z() - crossingAngleX * urad);  // in CMS ref. frame
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;  // Totem propagations assumes the starting point at 0 (zero)
0120     zout = fPPSRegionStart_45_;
0121   } else {
0122     approximator = m_aprox_ip_150_r_;
0123     zin = 0.0;  // Totem propagations assumes the starting point at 0 (zero)
0124     zout = fPPSRegionStart_56_;
0125   }
0126
0127   bool invert_beam_coord_system =
0128       true;  // it doesn't matter the option here, it is hard coded as TRUE inside LHCOpticsApproximator!
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];  // calculates px by means of TH_X, which is in the LHC ref. frame.
0161   double py = out_momentum[1];   // this need to be checked again, since it seems an invertion is occuring in the prop.
0162   double pz =
0163       out_momentum[2];  // totem calculates output pz already in the CMS ref. frame, it doesn't need to be converted
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;  // Totem parameterization uses meter, one need it in millimeter
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   // read parametrization
0192   LHCOpticsApproximator* aprox = (LHCOpticsApproximator*)f->Get(m_model_name.c_str());
0193   f->Close();
0194   return aprox;
0195 }