HectorTransport.cc

Go to the documentation of this file.

#include "SimTransport/PPSProtonTransport/interface/HectorTransport.h"
#include "Utilities/PPS/interface/PPSUtilities.h"
#include <CLHEP/Random/RandGauss.h>
#include "TLorentzVector.h" 
//Hector headers
#include "H_BeamLine.h"
#include "H_BeamParticle.h"


HectorTransport::HectorTransport(){};

HectorTransport::HectorTransport(const edm::ParameterSet & param, bool verbosity)
{
   // Create LHC beam line
   MODE=TransportMode::HECTOR; // defines the MODE for the transport
   m_verbosity=verbosity;
   edm::ParameterSet hector_par = param.getParameter<edm::ParameterSet>("PPSHector");
   
   m_f_ctpps_f          = (float) hector_par.getParameter<double>("PPSf");
   m_b_ctpps_b          = (float) hector_par.getParameter<double>("PPSb");
   fCrossingAngle_56    = hector_par.getParameter<double>("CrossingAngleBeam1");
   fCrossingAngle_45    = hector_par.getParameter<double>("CrossingAngleBeam2");
   fBeamEnergy          = hector_par.getParameter<double>("BeamEnergy"); // beam energy in GeV
   m_fEtacut              = hector_par.getParameter<double>("EtaCutForHector");
   m_fMomentumMin         = hector_par.getParameter<double>("MomentumMin");
   fBeamMomentum        = sqrt(fBeamEnergy*fBeamEnergy - PPSTools::ProtonMassSQ);

    // User definitons
    m_lengthctpps     = hector_par.getParameter<double>("BeamLineLengthPPS" );

    m_beam1filename   = hector_par.getParameter<string>("Beam1");
    m_beam2filename   = hector_par.getParameter<string>("Beam2");  
    fVtxMeanX       = param.getParameter<double>("VtxMeanX");
    fVtxMeanY       = param.getParameter<double>("VtxMeanY");
    fVtxMeanZ       = param.getParameter<double>("VtxMeanZ");
    m_sigmaSTX      = hector_par.getParameter<double>("sigmaSTX" );
    m_sigmaSTY      = hector_par.getParameter<double>("sigmaSTY" );
    m_sigmaSX       = hector_par.getParameter<double>("sigmaSX");
    m_sigmaSY       = hector_par.getParameter<double>("sigmaSY");
    m_sig_E         = hector_par.getParameter<double>("sigmaEnergy");
    fBeamXatIP      = hector_par.getUntrackedParameter<double>("BeamXatIP",fVtxMeanX);
    fBeamYatIP      = hector_par.getUntrackedParameter<double>("BeamYatIP",fVtxMeanY);
    bApplyZShift    = hector_par.getParameter<bool>("ApplyZShift");
    //PPS
    edm::LogInfo("ProtonTransport")
         <<"=============================================================================\n"
         <<"             Bulding LHC Proton transporter based on HECTOR model\n"
         <<"=============================================================================\n";
    setBeamLine();
    fPPSRegionStart_56=m_b_ctpps_b;
    fPPSRegionStart_45=m_f_ctpps_f;
}
HectorTransport::~HectorTransport()
{
        this->clear();

}
void HectorTransport::process( const HepMC::GenEvent * ev , const edm::EventSetup & iSetup, CLHEP::HepRandomEngine * _engine)
{
     engine = _engine;
     iSetup.getData( m_pdt );
     genProtonsLoop(ev,iSetup);
     addPartToHepMC(const_cast<HepMC::GenEvent*>(ev));
}
bool HectorTransport::transportProton(const HepMC::GenParticle* gpart)
{
     edm::LogInfo("ProtonTransport")<<"Starting proton transport using HECTOR method\n";

     double px,py,pz,e;
     unsigned int line = (gpart)->barcode();

     double mass = gpart->generatedMass();
     double charge = 1.;

     px = gpart->momentum().px();     
     py = gpart->momentum().py();     
     pz = gpart->momentum().pz();     
     e =  gpart->momentum().e();      

     e = sqrt(pow(mass,2)+pow(px,2)+pow(py,2)+pow(pz,2));

     int direction = (pz>0)?1:-1;

     // Apply Beam and Crossing Angle Corrections
     TLorentzVector p_out(px,py,pz,e);
     PPSTools::LorentzBoost(p_out,"LAB", {fCrossingAngle_56 ,//Beam1
           fCrossingAngle_45, //Beam2
           fBeamMomentum,fBeamEnergy});

     ApplyBeamCorrection(p_out);

     // from mm to cm        
     double XforPosition = gpart->production_vertex()->position().x()/cm;//cm
     double YforPosition = gpart->production_vertex()->position().y()/cm;//cm
     double ZforPosition = gpart->production_vertex()->position().z()/cm;//cm

     H_BeamParticle h_p(mass,charge);
     h_p.set4Momentum(-direction*p_out.Px(), p_out.Py(), fabs(p_out.Pz()), p_out.E());
// shift the beam position to the given beam position at IP (in cm)
     XforPosition=(XforPosition-fVtxMeanX)+fBeamXatIP*mm_to_cm;
     YforPosition=(YforPosition-fVtxMeanY)+fBeamYatIP*mm_to_cm;
     //ZforPosition stays the same, move the closed orbit only in the X,Y plane
//
// shift the starting position of the track to Z=0 if configured so (all the variables below are still in cm)
     if (bApplyZShift) {
        double fCrossingAngle = (p_out.Pz()>0)?fCrossingAngle_45:-fCrossingAngle_56;
        XforPosition = XforPosition+(tan((long double)fCrossingAngle*urad)-((long double)p_out.Px())/((long double)p_out.Pz()))*ZforPosition;
        YforPosition = YforPosition-((long double)p_out.Py())/((long double)p_out.Pz())*ZforPosition;
        ZforPosition = 0.;
     }

//
// set position, but do not invert the coordinate for the angles (TX,TY) because it is done by set4Momentum
     h_p.setPosition(-direction*XforPosition*cm_to_um,YforPosition*cm_to_um,h_p.getTX(),h_p.getTY(),-direction*ZforPosition*cm_to_m);
     bool is_stop;
     float x1_ctpps;
     float y1_ctpps;

     H_BeamLine* _beamline = nullptr;
     double _targetZ=0;
     switch (direction) { 
             case -1: _beamline = &*m_beamline56;// negative side propagation
                      _targetZ  = m_b_ctpps_b;
                      break;
             case 1:  _beamline = &*m_beamline45;
                      _targetZ  = m_f_ctpps_f;
                      break;
     }
     // insert protection for NULL beamlines here
     h_p.computePath(&*_beamline);
     is_stop = h_p.stopped(&*_beamline);
     if(m_verbosity) LogDebug("HectorTransportEventProcessing") << "HectorTransport:filterPPS: barcode = "
             << line << " is_stop=  "<< is_stop;
     if (is_stop) {
       return false;
     }
     //
     //propagating
     //
     h_p.propagate( _targetZ );

     p_out = PPSTools::HectorParticle2LorentzVector(h_p,direction);

     p_out.SetPx(direction*p_out.Px());
     x1_ctpps = direction*h_p.getX()*um_to_mm; 
     y1_ctpps = h_p.getY()*um_to_mm;

     if(m_verbosity) LogDebug("HectorTransportEventProcessing") <<
             "HectorTransport:filterPPS: barcode = " << line << " x=  "<< x1_ctpps <<" y= " << y1_ctpps;
         
     m_beamPart[line]    = p_out;
     m_xAtTrPoint[line]  = x1_ctpps;
     m_yAtTrPoint[line]  = y1_ctpps;
     return true;
}
void HectorTransport::genProtonsLoop( const HepMC::GenEvent * evt ,const edm::EventSetup & iSetup)
{
/*
   Loop over genVertex looking for transportable protons
*/
    unsigned int line;

    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())>m_fEtacut && fabs((*eventParticle)->momentum().pz())>m_fMomentumMin)) continue;
        line = (*eventParticle)->barcode();

        if ( m_beamPart.find(line) != m_beamPart.end() ) continue;

        HepMC::GenParticle * gpart = (*eventParticle);  

        transportProton(gpart);
    } 
}
bool HectorTransport::setBeamLine()
{

    m_beamline45=nullptr;
    m_beamline56=nullptr;
    edm::FileInPath b1(m_beam1filename.c_str());
    edm::FileInPath b2(m_beam2filename.c_str());
    if(m_verbosity) {
        edm::LogInfo("HectorTransportSetup") << "===================================================================\n"  
            << " * * * * * * * * * * * * * * * * * * * * * * * * * * * *           \n"  
            << " *                                                         *       \n"  
            << " *                   --<--<--  A fast simulator --<--<--     *     \n"  
            << " *                 | --<--<--     of particle   --<--<--     *     \n"  
            << " *  ----HECTOR----<                                          *     \n"  
            << " *                 | -->-->-- transport through-->-->--      *     \n"   
            << " *                   -->-->-- generic beamlines -->-->--     *     \n"  
            << " *                                                           *     \n"   
            << " * JINST 2:P09005 (2007)                                     *     \n"  
            << " *      X Rouby, J de Favereau, K Piotrzkowski (CP3)         *     \n"  
            << " *       http://www.fynu.ucl.ac.be/hector.html               *     \n"  
            << " *                                                           *     \n"  
            << " * Center for Cosmology, Particle Physics and Phenomenology  *     \n"  
            << " *              Universite catholique de Louvain             *     \n"  
            << " *                 Louvain-la-Neuve, Belgium                 *     \n"  
            << " *                                                         *       \n"  
            << " * * * * * * * * * * * * * * * * * * * * * * * * * * * *           \n"   
            << " HectorTransport configuration: \n" 
            << " lengthctpps      = " << m_lengthctpps << "\n"
            << " m_f_ctpps_f      =  " << m_f_ctpps_f << "\n"
            << " m_b_ctpps_b      =  " << m_b_ctpps_b << "\n"
            << "===================================================================\n";
    }  

    // construct beam line for PPS (forward 1 backward 2):                                                                                           
    if(m_lengthctpps>0. ) {
        m_beamline45 = std::unique_ptr<H_BeamLine>(new H_BeamLine(-1, m_lengthctpps + 0.1 )); // (direction, length)
        m_beamline45->fill( b2.fullPath(), 1, "IP5");
        m_beamline56 = std::unique_ptr<H_BeamLine>(new H_BeamLine( 1, m_lengthctpps + 0.1 )); //
        m_beamline56->fill( b1.fullPath(), 1, "IP5");
    }
    else {
        if ( m_verbosity ) LogDebug("HectorTransportSetup") << "HectorTransport: WARNING: lengthctpps=  " << m_lengthctpps;
        return false;
    }
    if (m_verbosity) {
          edm::LogInfo("HectorTransportSetup")
                     << "====================================================================\n"
                     << "                  Forward beam line elements \n";
          m_beamline45->showElements();
          edm::LogInfo("HectorTransportSetup")
                    << "====================================================================\n"
                    << "                 Backward beam line elements \n";
          m_beamline56->showElements();
    }

    return true;
}