GeneratorProducts/src/HepMCProduct.cc

0001 /*************************
0002  *
0003  *  Date: 2005/08 $
0004  *  \author J Weng - F. Moortgat'
0005  */
0006
0007 #include <iostream>
0008 #include <algorithm>  // because we use std::swap
0009
0010 #include <HepMC/GenEvent.h>
0011 #include <HepMC/SimpleVector.h>
0012 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0013
0014 using namespace edm;
0015 using namespace std;
0016
0017 HepMCProduct::HepMCProduct(HepMC::GenEvent* evt)
0018     : evt_(evt), isVtxGenApplied_(false), isVtxBoostApplied_(false), isPBoostApplied_(false) {}
0019
0020 HepMCProduct::~HepMCProduct() {
0021   delete evt_;
0022   evt_ = nullptr;
0023   isVtxGenApplied_ = false;
0024   isVtxBoostApplied_ = false;
0025   isPBoostApplied_ = false;
0026 }
0027
0028 void HepMCProduct::addHepMCData(HepMC::GenEvent* evt) {
0029   //  evt->print();
0030   // cout <<sizeof (evt)  <<"  " << sizeof ( HepMC::GenEvent)   << endl;
0031   evt_ = evt;
0032
0033   // same story about vertex smearing - GenEvent won't know it...
0034   // in fact, would be better to implement CmsGenEvent...
0035 }
0036
0037 void HepMCProduct::applyVtxGen(HepMC::FourVector const& vtxShift) {
0038   //std::cout<< " applyVtxGen called " << isVtxGenApplied_ << endl;
0039   //fTimeOffset = 0;
0040
0041   if (isVtxGenApplied())
0042     return;
0043
0044   for (HepMC::GenEvent::vertex_iterator vt = evt_->vertices_begin(); vt != evt_->vertices_end(); ++vt) {
0045     double x = (*vt)->position().x() + vtxShift.x();
0046     double y = (*vt)->position().y() + vtxShift.y();
0047     double z = (*vt)->position().z() + vtxShift.z();
0048     double t = (*vt)->position().t() + vtxShift.t();
0049     //std::cout << " vertex (x,y,z)= " << x <<" " << y << " " << z << std::endl;
0050     (*vt)->set_position(HepMC::FourVector(x, y, z, t));
0051   }
0052
0053   isVtxGenApplied_ = true;
0054
0055   return;
0056 }
0057
0058 void HepMCProduct::boostToLab(TMatrixD const* lorentz, std::string const& type) {
0059   //std::cout << "from boostToLab:" << std::endl;
0060
0061   if (lorentz == nullptr) {
0062     //std::cout << " lorentz = 0 " << std::endl;
0063     return;
0064   }
0065
0066   //lorentz->Print();
0067
0068   TMatrixD tmplorentz(*lorentz);
0069   //tmplorentz.Print();
0070
0071   if (type == "vertex") {
0072     if (isVtxBoostApplied()) {
0073       //std::cout << " isVtxBoostApplied true " << std::endl;
0074       return;
0075     }
0076
0077     for (HepMC::GenEvent::vertex_iterator vt = evt_->vertices_begin(); vt != evt_->vertices_end(); ++vt) {
0078       // change basis to lorentz boost definition: (t,x,z,y)
0079       TMatrixD p4(4, 1);
0080       p4(0, 0) = (*vt)->position().t();
0081       p4(1, 0) = (*vt)->position().x();
0082       p4(2, 0) = (*vt)->position().z();
0083       p4(3, 0) = (*vt)->position().y();
0084
0085       TMatrixD p4lab(4, 1);
0086       p4lab = tmplorentz * p4;
0087       //std::cout << " vertex lorentz: " << p4lab(1,0) << " " << p4lab(3,0) << " " << p4lab(2,0) << std::endl;
0088       (*vt)->set_position(HepMC::FourVector(p4lab(1, 0), p4lab(3, 0), p4lab(2, 0), p4lab(0, 0)));
0089     }
0090
0091     isVtxBoostApplied_ = true;
0092   } else if (type == "momentum") {
0093     if (isPBoostApplied()) {
0094       //std::cout << " isPBoostApplied true " << std::endl;
0095       return;
0096     }
0097
0098     for (HepMC::GenEvent::particle_iterator part = evt_->particles_begin(); part != evt_->particles_end(); ++part) {
0099       // change basis to lorentz boost definition: (E,Px,Pz,Py)
0100       TMatrixD p4(4, 1);
0101       p4(0, 0) = (*part)->momentum().e();
0102       p4(1, 0) = (*part)->momentum().x();
0103       p4(2, 0) = (*part)->momentum().z();
0104       p4(3, 0) = (*part)->momentum().y();
0105
0106       TMatrixD p4lab(4, 1);
0107       p4lab = tmplorentz * p4;
0108       //std::cout << " momentum lorentz: " << p4lab(1,0) << " " << p4lab(3,0) << " " << p4lab(2,0) << std::endl;
0109       (*part)->set_momentum(HepMC::FourVector(p4lab(1, 0), p4lab(3, 0), p4lab(2, 0), p4lab(0, 0)));
0110     }
0111
0112     isPBoostApplied_ = true;
0113   } else {
0114     std::cout << " no type found for boostToLab(std::string), options are vertex or momentum" << std::endl;
0115   }
0116
0117   return;
0118 }
0119
0120 const HepMC::GenEvent& HepMCProduct::getHepMCData() const { return *evt_; }
0121
0122 // copy constructor
0123 HepMCProduct::HepMCProduct(HepMCProduct const& other) : evt_(nullptr) {
0124   if (other.evt_)
0125     evt_ = new HepMC::GenEvent(*other.evt_);
0126   isVtxGenApplied_ = other.isVtxGenApplied_;
0127   isVtxBoostApplied_ = other.isVtxBoostApplied_;
0128   isPBoostApplied_ = other.isPBoostApplied_;
0129   //fTimeOffset = other.fTimeOffset;
0130 }
0131
0132 // swap
0133 void HepMCProduct::swap(HepMCProduct& other) {
0134   std::swap(evt_, other.evt_);
0135   std::swap(isVtxGenApplied_, other.isVtxGenApplied_);
0136   std::swap(isVtxBoostApplied_, other.isVtxBoostApplied_);
0137   std::swap(isPBoostApplied_, other.isPBoostApplied_);
0138   //std::swap(fTimeOffset, other.fTimeOffset);
0139 }
0140
0141 // assignment: use copy/swap idiom for exception safety.
0142 HepMCProduct& HepMCProduct::operator=(HepMCProduct const& other) {
0143   HepMCProduct temp(other);
0144   swap(temp);
0145   return *this;
0146 }
0147
0148 // move, needed explicitly as we have raw pointer...
0149 HepMCProduct::HepMCProduct(HepMCProduct&& other) : evt_(nullptr) { swap(other); }
0150 HepMCProduct& HepMCProduct::operator=(HepMCProduct&& other) {
0151   swap(other);
0152   return *this;
0153 }