ProtonTransport.cc

Go to the documentation of this file.

#include "SimTransport/PPSProtonTransport/interface/ProtonTransport.h"
#include "Utilities/PPS/interface/PPSUnitConversion.h"
#include <CLHEP/Random/RandGauss.h>
#include <CLHEP/Units/GlobalSystemOfUnits.h>

ProtonTransport::ProtonTransport(){};
ProtonTransport::~ProtonTransport(){};
void ProtonTransport::clear() {
  m_beamPart.clear();
  m_xAtTrPoint.clear();
  m_yAtTrPoint.clear();
};

void ProtonTransport::addPartToHepMC(HepMC::GenEvent* evt) {
  NEvent++;
  m_CorrespondenceMap.clear();

  int direction = 0;
  HepMC::GenParticle* gpart;

  unsigned int line;

  for (auto const& it : m_beamPart) {
    line = (it).first;
    gpart = evt->barcode_to_particle(line);

    direction = (gpart->momentum().pz() > 0) ? 1 : -1;

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

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

    //
    // ATTENTION: at this point, the vertex at PPS is already in mm
    //
    if (ddd == 0.)
      continue;
    if (m_verbosity) {
      LogDebug("ProtonTransportEventProcessing")
          << "ProtonTransport:: x= " << (*(m_xAtTrPoint.find(line))).second << "\n"
          << "ProtonTransport:: y= " << (*(m_yAtTrPoint.find(line))).second << "\n"
          << "ProtonTransport:: z= " << ddd * direction * m_to_mm << "\n"
          << "ProtonTransport:: t= " << time;
    }
    TLorentzVector const& p_out = (it).second;

    HepMC::GenVertex* vert = new HepMC::GenVertex(HepMC::FourVector((*(m_xAtTrPoint.find(line))).second,
                                                                    (*(m_yAtTrPoint.find(line))).second,
                                                                    ddd * direction * m_to_mm,
                                                                    time + time * 0.001));

    gpart->set_status(2);
    vert->add_particle_in(gpart);
    vert->add_particle_out(new HepMC::GenParticle(
        HepMC::FourVector(p_out.Px(), p_out.Py(), p_out.Pz(), p_out.E()), gpart->pdg_id(), 1, gpart->flow()));
    evt->add_vertex(vert);

    int ingoing = (*vert->particles_in_const_begin())->barcode();
    int outgoing = (*vert->particles_out_const_begin())->barcode();

    LHCTransportLink theLink(ingoing, outgoing);
    if (m_verbosity)
      LogDebug("ProtonTransportEventProcessing") << "ProtonTransport:addPartToHepMC: LHCTransportLink " << theLink;
    m_CorrespondenceMap.push_back(theLink);
  }
}
void ProtonTransport::ApplyBeamCorrection(HepMC::GenParticle* p) {
  TLorentzVector p_out;
  p_out.SetPx(p->momentum().px());
  p_out.SetPy(p->momentum().py());
  p_out.SetPz(p->momentum().pz());
  p_out.SetE(p->momentum().e());
  ApplyBeamCorrection(p_out);
  p->set_momentum(HepMC::FourVector(p_out.Px(), p_out.Py(), p_out.Pz(), p_out.E()));
}
void ProtonTransport::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) ? fCrossingAngle_45 * urad : fCrossingAngle_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 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);
  energy += denergy;
  double p = sqrt(pow(energy, 2) - ProtonMassSQ);

  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.SetE(energy);
}