GeneratorProducts/src/HepMC3Product.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 <HepMC3/GenEvent.h>
0011 #include <HepMC3/GenVertex.h>
0012 #include <HepMC3/GenParticle.h>
0013 #include "SimDataFormats/GeneratorProducts/interface/HepMC3Product.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015
0016 using namespace edm;
0017 using namespace std;
0018
0019 HepMC3Product::HepMC3Product(HepMC3::GenEvent* evt)
0020     : isVtxGenApplied_(false), isVtxBoostApplied_(false), isPBoostApplied_(false) {
0021   addHepMCData(evt);
0022 }
0023
0024 HepMC3Product::~HepMC3Product() = default;
0025
0026 void HepMC3Product::addHepMCData(HepMC3::GenEvent* evt) { evt->write_data(evt_); }
0027
0028 void HepMC3Product::applyVtxGen(HepMC3::FourVector const& vtxShift) {
0029   //std::cout<< " applyVtxGen called " << isVtxGenApplied_ << endl;
0030   //fTimeOffset = 0;
0031   if (isVtxGenApplied())
0032     return;
0033   HepMC3::GenEvent evt;
0034   evt.read_data(evt_);
0035   evt.shift_position_by(vtxShift);
0036   evt.write_data(evt_);
0037   isVtxGenApplied_ = true;
0038   return;
0039 }
0040
0041 void HepMC3Product::boostToLab(TMatrixD const* lorentz, std::string const& type) {
0042   //std::cout << "from boostToLab:" << std::endl;
0043   HepMC3::GenEvent evt;
0044   evt.read_data(evt_);
0045
0046   if (lorentz == nullptr) {
0047     //std::cout << " lorentz = 0 " << std::endl;
0048     return;
0049   }
0050
0051   //lorentz->Print();
0052
0053   TMatrixD tmplorentz(*lorentz);
0054   //tmplorentz.Print();
0055
0056   if (type == "vertex") {
0057     if (isVtxBoostApplied()) {
0058       //std::cout << " isVtxBoostApplied true " << std::endl;
0059       return;
0060     }
0061
0062     for (const HepMC3::GenVertexPtr& vt : evt.vertices()) {
0063       // change basis to lorentz boost definition: (t,x,z,y)
0064       TMatrixD p4(4, 1);
0065       p4(0, 0) = vt->position().t();
0066       p4(1, 0) = vt->position().x();
0067       p4(2, 0) = vt->position().z();
0068       p4(3, 0) = vt->position().y();
0069
0070       TMatrixD p4lab(4, 1);
0071       p4lab = tmplorentz * p4;
0072       //std::cout << " vertex lorentz: " << p4lab(1,0) << " " << p4lab(3,0) << " " << p4lab(2,0) << std::endl;
0073       vt->set_position(HepMC3::FourVector(p4lab(1, 0), p4lab(3, 0), p4lab(2, 0), p4lab(0, 0)));
0074     }
0075
0076     isVtxBoostApplied_ = true;
0077   } else if (type == "momentum") {
0078     if (isPBoostApplied()) {
0079       //std::cout << " isPBoostApplied true " << std::endl;
0080       return;
0081     }
0082
0083     for (const HepMC3::GenParticlePtr& part : evt.particles()) {
0084       // change basis to lorentz boost definition: (E,Px,Pz,Py)
0085       TMatrixD p4(4, 1);
0086       p4(0, 0) = part->momentum().e();
0087       p4(1, 0) = part->momentum().x();
0088       p4(2, 0) = part->momentum().z();
0089       p4(3, 0) = part->momentum().y();
0090
0091       TMatrixD p4lab(4, 1);
0092       p4lab = tmplorentz * p4;
0093       //std::cout << " momentum lorentz: " << p4lab(1,0) << " " << p4lab(3,0) << " " << p4lab(2,0) << std::endl;
0094       part->set_momentum(HepMC3::FourVector(p4lab(1, 0), p4lab(3, 0), p4lab(2, 0), p4lab(0, 0)));
0095     }
0096
0097     isPBoostApplied_ = true;
0098   } else {
0099     edm::LogWarning("GeneratorProducts") << "no type found for boostToLab(std::string), options are vertex or momentum";
0100     //throw edm::Exception(edm::errors::Configuration, "GeneratorProducts")
0101     //  << "no type found for boostToLab(std::string), options are vertex or momentum \n";
0102   }
0103
0104   evt.write_data(evt_);
0105   return;
0106 }