HepMCProduct.cc

Go to the documentation of this file.

/*************************
 *
 *  Date: 2005/08 $
 *  \author J Weng - F. Moortgat'
 */

#include <iostream>
#include <algorithm>  // because we use std::swap

#include <HepMC/GenEvent.h>
#include <HepMC/SimpleVector.h>
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"

using namespace edm;
using namespace std;

HepMCProduct::HepMCProduct(HepMC::GenEvent* evt)
    : evt_(evt), isVtxGenApplied_(false), isVtxBoostApplied_(false), isPBoostApplied_(false) {}

HepMCProduct::~HepMCProduct() {
  delete evt_;
  evt_ = nullptr;
  isVtxGenApplied_ = false;
  isVtxBoostApplied_ = false;
  isPBoostApplied_ = false;
}

void HepMCProduct::addHepMCData(HepMC::GenEvent* evt) {
  //  evt->print();
  // cout <<sizeof (evt)  <<"  " << sizeof ( HepMC::GenEvent)   << endl;
  evt_ = evt;

  // same story about vertex smearing - GenEvent won't know it...
  // in fact, would be better to implement CmsGenEvent...
}

void HepMCProduct::applyVtxGen(HepMC::FourVector const& vtxShift) {
  //std::cout<< " applyVtxGen called " << isVtxGenApplied_ << endl;
  //fTimeOffset = 0;

  if (isVtxGenApplied())
    return;

  for (HepMC::GenEvent::vertex_iterator vt = evt_->vertices_begin(); vt != evt_->vertices_end(); ++vt) {
    double x = (*vt)->position().x() + vtxShift.x();
    double y = (*vt)->position().y() + vtxShift.y();
    double z = (*vt)->position().z() + vtxShift.z();
    double t = (*vt)->position().t() + vtxShift.t();
    //std::cout << " vertex (x,y,z)= " << x <<" " << y << " " << z << std::endl;
    (*vt)->set_position(HepMC::FourVector(x, y, z, t));
  }

  isVtxGenApplied_ = true;

  return;
}

void HepMCProduct::boostToLab(TMatrixD const* lorentz, std::string const& type) {
  //std::cout << "from boostToLab:" << std::endl;

  if (lorentz == nullptr) {
    //std::cout << " lorentz = 0 " << std::endl;
    return;
  }

  //lorentz->Print();

  TMatrixD tmplorentz(*lorentz);
  //tmplorentz.Print();

  if (type == "vertex") {
    if (isVtxBoostApplied()) {
      //std::cout << " isVtxBoostApplied true " << std::endl;
      return;
    }

    for (HepMC::GenEvent::vertex_iterator vt = evt_->vertices_begin(); vt != evt_->vertices_end(); ++vt) {
      // change basis to lorentz boost definition: (t,x,z,y)
      TMatrixD p4(4, 1);
      p4(0, 0) = (*vt)->position().t();
      p4(1, 0) = (*vt)->position().x();
      p4(2, 0) = (*vt)->position().z();
      p4(3, 0) = (*vt)->position().y();

      TMatrixD p4lab(4, 1);
      p4lab = tmplorentz * p4;
      //std::cout << " vertex lorentz: " << p4lab(1,0) << " " << p4lab(3,0) << " " << p4lab(2,0) << std::endl;
      (*vt)->set_position(HepMC::FourVector(p4lab(1, 0), p4lab(3, 0), p4lab(2, 0), p4lab(0, 0)));
    }

    isVtxBoostApplied_ = true;
  } else if (type == "momentum") {
    if (isPBoostApplied()) {
      //std::cout << " isPBoostApplied true " << std::endl;
      return;
    }

    for (HepMC::GenEvent::particle_iterator part = evt_->particles_begin(); part != evt_->particles_end(); ++part) {
      // change basis to lorentz boost definition: (E,Px,Pz,Py)
      TMatrixD p4(4, 1);
      p4(0, 0) = (*part)->momentum().e();
      p4(1, 0) = (*part)->momentum().x();
      p4(2, 0) = (*part)->momentum().z();
      p4(3, 0) = (*part)->momentum().y();

      TMatrixD p4lab(4, 1);
      p4lab = tmplorentz * p4;
      //std::cout << " momentum lorentz: " << p4lab(1,0) << " " << p4lab(3,0) << " " << p4lab(2,0) << std::endl;
      (*part)->set_momentum(HepMC::FourVector(p4lab(1, 0), p4lab(3, 0), p4lab(2, 0), p4lab(0, 0)));
    }

    isPBoostApplied_ = true;
  } else {
    std::cout << " no type found for boostToLab(std::string), options are vertex or momentum" << std::endl;
  }

  return;
}

const HepMC::GenEvent& HepMCProduct::getHepMCData() const { return *evt_; }

// copy constructor
HepMCProduct::HepMCProduct(HepMCProduct const& other) : evt_(nullptr) {
  if (other.evt_)
    evt_ = new HepMC::GenEvent(*other.evt_);
  isVtxGenApplied_ = other.isVtxGenApplied_;
  isVtxBoostApplied_ = other.isVtxBoostApplied_;
  isPBoostApplied_ = other.isPBoostApplied_;
  //fTimeOffset = other.fTimeOffset;
}

// swap
void HepMCProduct::swap(HepMCProduct& other) {
  std::swap(evt_, other.evt_);
  std::swap(isVtxGenApplied_, other.isVtxGenApplied_);
  std::swap(isVtxBoostApplied_, other.isVtxBoostApplied_);
  std::swap(isPBoostApplied_, other.isPBoostApplied_);
  //std::swap(fTimeOffset, other.fTimeOffset);
}

// assignment: use copy/swap idiom for exception safety.
HepMCProduct& HepMCProduct::operator=(HepMCProduct const& other) {
  HepMCProduct temp(other);
  swap(temp);
  return *this;
}

// move, needed explicitly as we have raw pointer...
HepMCProduct::HepMCProduct(HepMCProduct&& other) : evt_(nullptr) { swap(other); }
HepMCProduct& HepMCProduct::operator=(HepMCProduct&& other) {
  swap(other);
  return *this;
}