Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:31:08 

0001 #include "SimTransport/PPSProtonTransport/interface/HectorTransport.h"
0002 #include "Utilities/PPS/interface/PPSUtilities.h"
0003 #include "TLorentzVector.h"
0004 //Hector headers
0005 #include "H_BeamLine.h"
0006 #include "H_BeamParticle.h"
0007 #include <memory>
0008 
0009 #include <string>
0010 
0011 HectorTransport::HectorTransport(const edm::ParameterSet& iConfig, edm::ConsumesCollector iC)
0012     : BaseProtonTransport(iConfig),
0013       m_beamline45(nullptr),
0014       m_beamline56(nullptr),
0015       beamParametersToken_(iC.esConsumes()),
0016       beamspotToken_(iC.esConsumes()) {
0017   MODE = TransportMode::HECTOR;
0018   std::string s1 = iConfig.getParameter<std::string>("Beam1Filename");
0019   std::string s2 = iConfig.getParameter<std::string>("Beam2Filename");
0020   setBeamFileNames(s1, s2);
0021   double cax45 = iConfig.getParameter<double>("halfCrossingAngleXSector45");
0022   double cax56 = iConfig.getParameter<double>("halfCrossingAngleXSector56");
0023   double cay45 = iConfig.getParameter<double>("halfCrossingAngleYSector45");
0024   double cay56 = iConfig.getParameter<double>("halfCrossingAngleYSector56");
0025   setCrossingAngles(cax45, cax56, cay45, cay56);
0026   double stx = iConfig.getParameter<double>("BeamDivergenceX");
0027   double sty = iConfig.getParameter<double>("BeamDivergenceY");
0028   double sx = iConfig.getParameter<double>("BeamSigmaX");
0029   double sy = iConfig.getParameter<double>("BeamSigmaY");
0030   double se = iConfig.getParameter<double>("BeamEnergyDispersion");
0031   setBeamParameters(stx, sty, sx, sy, se);
0032   //PPS
0033   edm::LogVerbatim("ProtonTransport")
0034       << "=============================================================================\n"
0035       << "             Bulding LHC Proton transporter based on HECTOR model\n"
0036       << "=============================================================================\n";
0037   setBeamLine();
0038 }
0039 HectorTransport::~HectorTransport() {}
0040 //
0041 // this method is the same for all propagator, but since transportProton is different for each derived class
0042 // it needes to be overriden
0043 //
0044 void HectorTransport::process(const HepMC::GenEvent* evt,
0045                               const edm::EventSetup& iSetup,
0046                               CLHEP::HepRandomEngine* engine) {
0047   clear();
0048   engine_ = engine;  // the engine needs to be updated for each event
0049 
0050   beamspot_ = &iSetup.getData(beamspotToken_);
0051   beamParameters_ = &iSetup.getData(beamParametersToken_);
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 
0065     if (gpart->pdg_id() != 2212)
0066       continue;  // only transport stable protons
0067     if (gpart->status() != 1)
0068       continue;
0069     if (m_beamPart.find(line) != m_beamPart.end())
0070       continue;  // assures this protons has not been already propagated
0071 
0072     transportProton(gpart);
0073   }
0074 }
0075 bool HectorTransport::transportProton(const HepMC::GenParticle* gpart) {
0076   edm::LogVerbatim("ProtonTransport") << "Starting proton transport using HECTOR method\n";
0077 
0078   double px, py, pz, e;
0079   unsigned int line = (gpart)->barcode();
0080 
0081   double mass = gpart->generatedMass();
0082   double charge = 1.;
0083 
0084   // remove the CMS vertex shift
0085   double vtxXoffset;
0086   double vtxYoffset;
0087   if (useBeamPositionFromLHCInfo_) {
0088     vtxXoffset = beamParameters_->getVtxOffsetX45() * cm_to_mm;
0089     vtxYoffset = beamParameters_->getVtxOffsetY45() * cm_to_mm;
0090   } else {
0091     vtxXoffset = beamspot_->x() * cm_to_mm;
0092     vtxYoffset = beamspot_->y() * cm_to_mm;
0093   }
0094 
0095   // Momentum in LHC ref. frame
0096   px = -gpart->momentum().px();
0097   py = gpart->momentum().py();
0098   pz = -gpart->momentum().pz();
0099   e = gpart->momentum().e();
0100 
0101   int direction = (pz > 0) ? 1 : -1;  // in relation to LHC ref frame
0102 
0103   double XforPosition = -gpart->production_vertex()->position().x();  //mm in the ref. frame of LHC
0104   double YforPosition = gpart->production_vertex()->position().y();   //mm
0105   double ZforPosition = -gpart->production_vertex()->position().z();  //mm
0106   double fCrossingAngleX = (pz < 0) ? fCrossingAngleX_45_ : fCrossingAngleX_56_;
0107   double fCrossingAngleY = (pz < 0) ? fCrossingAngleY_45_ : fCrossingAngleY_56_;
0108   //
0109   H_BeamParticle h_p(mass, charge);
0110   h_p.set4Momentum(px, py, pz, e);
0111   //
0112   // shift the starting position of the track to Z=0 if configured so (all the variables below are still in cm)
0113   XforPosition = XforPosition - ZforPosition * (px / pz + fCrossingAngleX * urad);  // theta ~=tan(theta)
0114   YforPosition = YforPosition - ZforPosition * (py / pz + fCrossingAngleY * urad);
0115   XforPosition -= (-vtxXoffset);  // X was already in the LHC ref. frame
0116   YforPosition -= vtxYoffset;
0117 
0118   // set position, but do not invert the coordinate for the angles (TX,TY) because it is done by set4Momentum
0119   h_p.setPosition(XforPosition * mm_to_um, YforPosition * mm_to_um, h_p.getTX(), h_p.getTY(), 0.);
0120 
0121   bool is_stop;
0122   float x1_ctpps;
0123   float y1_ctpps;
0124 
0125   H_BeamLine* beamline = nullptr;
0126   double targetZ = 0;
0127   switch (direction) {
0128     case 1:
0129       beamline = &*m_beamline56;  // negative side propagation
0130       targetZ = fPPSRegionStart_56_;
0131       break;
0132     case -1:
0133       beamline = &*m_beamline45;
0134       targetZ = fPPSRegionStart_45_;
0135       break;
0136   }
0137   // insert protection for NULL beamlines here
0138   h_p.computePath(&*beamline);
0139   is_stop = h_p.stopped(&*beamline);
0140   if (verbosity_)
0141     LogDebug("HectorTransportEventProcessing")
0142         << "HectorTransport:filterPPS: barcode = " << line << " is_stop=  " << is_stop;
0143   if (is_stop) {
0144     return false;
0145   }
0146   //
0147   //propagating
0148   //
0149   h_p.propagate(targetZ);
0150   x1_ctpps = h_p.getX();
0151   y1_ctpps = h_p.getY();
0152 
0153   double thx = h_p.getTX();
0154   double thy = h_p.getTY();
0155 
0156   if (produceHitsRelativeToBeam_) {
0157     H_BeamParticle p_beam(mass, charge);
0158     p_beam.set4Momentum(0., 0., beamMomentum(), beamEnergy());
0159     p_beam.setPosition(0., 0., fCrossingAngleX * urad, fCrossingAngleY * urad, 0.);
0160     p_beam.computePath(&*beamline);
0161     thx -= p_beam.getTX();
0162     thy -= p_beam.getTY();
0163     x1_ctpps -= p_beam.getX();
0164     y1_ctpps -= p_beam.getY();
0165     edm::LogVerbatim("HectorTransportEventProcessing")
0166         << "Shifting the hit relative to beam :  X = " << p_beam.getX() << "(mm)      Y = " << p_beam.getY() << "(mm)";
0167   }
0168 
0169   double partP = sqrt(pow(h_p.getE(), 2) - ProtonMassSQ);
0170   double theta = sqrt(thx * thx + thy * thy) * urad;
0171 
0172   // copy the kinematic changing to CMS ref. frame, only the negative Pz needs to be changed
0173   TLorentzVector p_out(-tan(thx * urad) * partP * cos(theta),
0174                        tan(thy * urad) * partP * cos(theta),
0175                        -direction * partP * cos(theta),
0176                        h_p.getE());
0177 
0178   m_beamPart[line] = p_out;
0179   m_xAtTrPoint[line] = -x1_ctpps * um_to_mm;  // move to CMS ref. frame
0180   m_yAtTrPoint[line] = y1_ctpps * um_to_mm;
0181   if (verbosity_)
0182     LogDebug("HectorTransportEventProcessing")
0183         << "HectorTransport:filterPPS: barcode = " << line << " x=  " << x1_ctpps << " y= " << y1_ctpps;
0184 
0185   return true;
0186 }
0187 bool HectorTransport::setBeamLine() {
0188   edm::FileInPath b1(beam1Filename_.c_str());
0189   edm::FileInPath b2(beam2Filename_.c_str());
0190 
0191   // construct beam line for PPS (forward 1 backward 2):
0192   if (fPPSBeamLineLength_ > 0.) {
0193     m_beamline45 = std::make_unique<H_BeamLine>(-1, fPPSBeamLineLength_ + 0.1);
0194     m_beamline45->fill(b2.fullPath(), 1, "IP5");
0195     m_beamline56 = std::make_unique<H_BeamLine>(1, fPPSBeamLineLength_ + 0.1);
0196     m_beamline56->fill(b1.fullPath(), 1, "IP5");
0197   } else {
0198     if (verbosity_)
0199       LogDebug("HectorTransportSetup") << "HectorTransport: WARNING: lengthctpps=  " << fPPSBeamLineLength_;
0200     return false;
0201   }
0202   if (verbosity_) {
0203     edm::LogInfo("HectorTransportSetup") << "===================================================================\n"
0204                                          << " * * * * * * * * * * * * * * * * * * * * * * * * * * * *           \n"
0205                                          << " *                                                         *       \n"
0206                                          << " *                   --<--<--  A fast simulator --<--<--     *     \n"
0207                                          << " *                 | --<--<--     of particle   --<--<--     *     \n"
0208                                          << " *  ----HECTOR----<                                          *     \n"
0209                                          << " *                 | -->-->-- transport through-->-->--      *     \n"
0210                                          << " *                   -->-->-- generic beamlines -->-->--     *     \n"
0211                                          << " *                                                           *     \n"
0212                                          << " * JINST 2:P09005 (2007)                                     *     \n"
0213                                          << " *      X Rouby, J de Favereau, K Piotrzkowski (CP3)         *     \n"
0214                                          << " *       http://www.fynu.ucl.ac.be/hector.html               *     \n"
0215                                          << " *                                                           *     \n"
0216                                          << " * Center for Cosmology, Particle Physics and Phenomenology  *     \n"
0217                                          << " *              Universite catholique de Louvain             *     \n"
0218                                          << " *                 Louvain-la-Neuve, Belgium                 *     \n"
0219                                          << " *                                                         *       \n"
0220                                          << " * * * * * * * * * * * * * * * * * * * * * * * * * * * *           \n"
0221                                          << " HectorTransport configuration: \n"
0222                                          << " Beam line length      = " << fPPSBeamLineLength_ << "\n"
0223                                          << " PPS Region Start 44   =  " << fPPSRegionStart_45_ << "\n"
0224                                          << " PPS Region Start 56   =  " << fPPSRegionStart_56_ << "\n"
0225                                          << "===================================================================\n";
0226     edm::LogVerbatim("HectorTransportSetup") << "===================================================================\n"
0227                                              << "                  Forward beam line elements \n";
0228     m_beamline45->showElements();
0229     edm::LogVerbatim("HectorTransportSetup") << "===================================================================\n"
0230                                              << "                 Backward beam line elements \n";
0231     m_beamline56->showElements();
0232   }
0233   return true;
0234 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!