EICInteractionVertexSmear.cxx

Go to the documentation of this file.

#include "DDG4/Factories.h"
 
// Framework include files
#include "DD4hep/Printout.h"
#include "DD4hep/InstanceCount.h"
#include "DDG4/Geant4Random.h"
#include "DDG4/Geant4Context.h"
#include "DDG4/Geant4InputHandling.h"
#include "npdet/EICInteractionVertexSmear.h"
 
// C/C++ include files
#include <cmath>
 
#include "Math/Vector3D.h"
#include "Math/Vector4D.h"
#include "Math/LorentzRotation.h"
#include "Math/AxisAngle.h"
#include "Math/Rotation3D.h"
#include "Math/RotationX.h"
#include "Math/RotationY.h"
#include "Math/Boost.h"
#include "Math/BoostX.h"
#include "Math/BoostX.h"
#include "Math/PxPyPzM4D.h"
 
namespace npdet::sim {
 
using namespace dd4hep::sim;
 
EICInteractionVertexSmear::EICInteractionVertexSmear(Geant4Context* ctxt, const std::string& nam)
 : Geant4GeneratorAction(ctxt, nam)
{
 dd4hep::InstanceCount::increment(this);
 declareProperty("Offset", m_offset);
 declareProperty("Sigma_ion", m_sigma_Ion);
 declareProperty("Sigma_e", m_sigma_Electron);
 declareProperty("Mask", m_mask = 0);
 m_needsControl = true;
}
 
EICInteractionVertexSmear::~EICInteractionVertexSmear() { dd4hep::InstanceCount::decrement(this); }
 
void EICInteractionVertexSmear::smear(Interaction* inter) const {
 using namespace ROOT::Math;
 Geant4Random& rndm = context()->event().random();
 if (inter) {
 double dx = rndm.gauss(m_offset.x(), m_sigma_Ion.x());
 double dy = rndm.gauss(m_offset.y(), m_sigma_Ion.y());
 double dz = rndm.gauss(m_offset.z(), m_sigma_Ion.z());
 double dt = rndm.gauss(m_offset.t(), m_sigma_Ion.t());
 double dxe = rndm.gauss(m_offset.x(), m_sigma_Electron.x());
 double dye = rndm.gauss(m_offset.y(), m_sigma_Electron.y());
 double dze = rndm.gauss(m_offset.z(), m_sigma_Electron.z());
 double dte = rndm.gauss(m_offset.t(), m_sigma_Electron.t());
 XYZVector ion_dir(0, 0, 1.0);
 RotationY rotx(-dx);
 RotationX roty(dy);
 auto new_ion_dir = rotx(roty(ion_dir));
 
 Geant4PrimaryEvent::Interaction::VertexMap::iterator iv;
 Geant4PrimaryEvent::Interaction::ParticleMap::iterator ip;
 
 double alpha = new_ion_dir.Theta();
 auto rotation_axis = ion_dir.Cross(new_ion_dir).Unit();
 
 if (inter->locked) {
 this->abortRun("Locked interactions may not be boosted!",
 "Cannot boost interactions with a native G4 primary record!");
 } else if (alpha != 0.0) {
 
 double tanalpha = std::tan(alpha/2.0);
 double gamma = std::sqrt(1 + tanalpha * tanalpha);
 double betagamma = tanalpha;
 double beta = betagamma / gamma;
 
 Polar3D boost_vec(beta, M_PI / 2.0, new_ion_dir.Phi());
 
 std::cout << " beta = " << beta << "\n";
 
 Boost bst(boost_vec);
 LorentzRotation roty(AxisAngle(rotation_axis, alpha / 2.0));
 
 // Now move begin and end-vertex of all primary vertices accordingly
 for (iv = inter->vertices.begin(); iv != inter->vertices.end(); ++iv) {
 for (Geant4Vertex* v : (*iv).second) {
 XYZTVector v0(v->x, v->y, v->z, v->time);
 auto v1 = roty(bst * v0);
 v->x = v1.x();
 v->y = v1.y();
 v->z = v1.z();
 v->time = v1.t();
 }
 }
 // Now move begin and end-vertex of all primary and generator particles accordingly
 for (ip = inter->particles.begin(); ip != inter->particles.end(); ++ip) {
 Geant4ParticleHandle p = (*ip).second;
 XYZTVector v0{p->vsx, p->vsy, p->vsz, p->time};
 auto v1 = roty(bst * v0);
 
 PxPyPzMVector p0{p->psx, p->psy, p->psz, p->mass};
 auto p1 = roty(bst * p0);
 
 p->vsx = v1.x();
 p->vsy = v1.y();
 p->vsz = v1.z();
 p->time = v1.t();
 
 p->psx = p1.x();
 p->psy = p1.y();
 p->psz = p1.z();
 }
 }
 
 // smearInteraction(this,inter,dx,dy,dz,dt);
 return;
 } else {
 //print("+++ Smearing primary vertex for interaction type %d (%d Vertices, %d particles) "
 // "by (%+.2e mm, %+.2e mm, %+.2e mm, %+.2e ns)",
 // m_mask, int(inter->vertices.size()), int(inter->particles.size()), dx, dy, dz, dt);
 print("+++ No interaction of type %d present.", m_mask);
 }
}
 
void EICInteractionVertexSmear::operator()(G4Event*) {
 typedef std::vector<Geant4PrimaryInteraction*> _I;
 Geant4PrimaryEvent* evt = context()->event().extension<Geant4PrimaryEvent>();
 
 if ( m_mask >= 0 ) {
 Interaction* inter = evt->get(m_mask);
 smear(inter);
 return;
 }
 _I interactions = evt->interactions();
 for(_I::iterator i=interactions.begin(); i != interactions.end(); ++i)
 smear(*i);
}
}
 
namespace dd4hep::sim {
using EICInteractionVertexSmear = npdet::sim::EICInteractionVertexSmear;
}
DECLARE_GEANT4ACTION(EICInteractionVertexSmear)