Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2025-01-11 03:38:36

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 
0040   if (isVtxGenApplied())
0041     return;
0042 
0043   for (HepMC::GenEvent::vertex_iterator vt = evt_->vertices_begin(); vt != evt_->vertices_end(); ++vt) {
0044     double x = (*vt)->position().x() + vtxShift.x();
0045     double y = (*vt)->position().y() + vtxShift.y();
0046     double z = (*vt)->position().z() + vtxShift.z();
0047     double t = (*vt)->position().t() + vtxShift.t();
0048     //std::cout << " vertex (x,y,z)= " << x <<" " << y << " " << z << std::endl;
0049     (*vt)->set_position(HepMC::FourVector(x, y, z, t));
0050   }
0051 
0052   isVtxGenApplied_ = true;
0053 
0054   return;
0055 }
0056 
0057 void HepMCProduct::boostToLab(TMatrixD const* lorentz, std::string const& type) {
0058   //std::cout << "from boostToLab:" << std::endl;
0059 
0060   if (lorentz == nullptr) {
0061     //std::cout << " lorentz = 0 " << std::endl;
0062     return;
0063   }
0064 
0065   //lorentz->Print();
0066 
0067   TMatrixD tmplorentz(*lorentz);
0068   //tmplorentz.Print();
0069 
0070   if (type == "vertex") {
0071     if (isVtxBoostApplied()) {
0072       //std::cout << " isVtxBoostApplied true " << std::endl;
0073       return;
0074     }
0075 
0076     for (HepMC::GenEvent::vertex_iterator vt = evt_->vertices_begin(); vt != evt_->vertices_end(); ++vt) {
0077       // change basis to lorentz boost definition: (t,x,z,y)
0078       TMatrixD p4(4, 1);
0079       p4(0, 0) = (*vt)->position().t();
0080       p4(1, 0) = (*vt)->position().x();
0081       p4(2, 0) = (*vt)->position().z();
0082       p4(3, 0) = (*vt)->position().y();
0083 
0084       TMatrixD p4lab(4, 1);
0085       p4lab = tmplorentz * p4;
0086       //std::cout << " vertex lorentz: " << p4lab(1,0) << " " << p4lab(3,0) << " " << p4lab(2,0) << std::endl;
0087       (*vt)->set_position(HepMC::FourVector(p4lab(1, 0), p4lab(3, 0), p4lab(2, 0), p4lab(0, 0)));
0088     }
0089 
0090     isVtxBoostApplied_ = true;
0091   } else if (type == "momentum") {
0092     if (isPBoostApplied()) {
0093       //std::cout << " isPBoostApplied true " << std::endl;
0094       return;
0095     }
0096 
0097     for (HepMC::GenEvent::particle_iterator part = evt_->particles_begin(); part != evt_->particles_end(); ++part) {
0098       // change basis to lorentz boost definition: (E,Px,Pz,Py)
0099       TMatrixD p4(4, 1);
0100       p4(0, 0) = (*part)->momentum().e();
0101       p4(1, 0) = (*part)->momentum().x();
0102       p4(2, 0) = (*part)->momentum().z();
0103       p4(3, 0) = (*part)->momentum().y();
0104 
0105       TMatrixD p4lab(4, 1);
0106       p4lab = tmplorentz * p4;
0107       //std::cout << " momentum lorentz: " << p4lab(1,0) << " " << p4lab(3,0) << " " << p4lab(2,0) << std::endl;
0108       (*part)->set_momentum(HepMC::FourVector(p4lab(1, 0), p4lab(3, 0), p4lab(2, 0), p4lab(0, 0)));
0109     }
0110 
0111     isPBoostApplied_ = true;
0112   } else {
0113     std::cout << " no type found for boostToLab(std::string), options are vertex or momentum" << std::endl;
0114   }
0115 
0116   return;
0117 }
0118 
0119 const HepMC::GenEvent& HepMCProduct::getHepMCData() const { return *evt_; }
0120 
0121 // copy constructor
0122 HepMCProduct::HepMCProduct(HepMCProduct const& other) : evt_(nullptr) {
0123   if (other.evt_)
0124     evt_ = new HepMC::GenEvent(*other.evt_);
0125   isVtxGenApplied_ = other.isVtxGenApplied_;
0126   isVtxBoostApplied_ = other.isVtxBoostApplied_;
0127   isPBoostApplied_ = other.isPBoostApplied_;
0128 }
0129 
0130 // swap
0131 void HepMCProduct::swap(HepMCProduct& other) {
0132   std::swap(evt_, other.evt_);
0133   std::swap(isVtxGenApplied_, other.isVtxGenApplied_);
0134   std::swap(isVtxBoostApplied_, other.isVtxBoostApplied_);
0135   std::swap(isPBoostApplied_, other.isPBoostApplied_);
0136 }
0137 
0138 // assignment: use copy/swap idiom for exception safety.
0139 HepMCProduct& HepMCProduct::operator=(HepMCProduct const& other) {
0140   HepMCProduct temp(other);
0141   swap(temp);
0142   return *this;
0143 }
0144 
0145 // move, needed explicitly as we have raw pointer...
0146 HepMCProduct::HepMCProduct(HepMCProduct&& other)
0147     : evt_(nullptr), isVtxGenApplied_(false), isVtxBoostApplied_(false), isPBoostApplied_(false) {
0148   swap(other);
0149 }
0150 HepMCProduct& HepMCProduct::operator=(HepMCProduct&& other) {
0151   swap(other);
0152   return *this;
0153 }