plugins/src/EICInteractionVertexSmear.cxx

0001 #include "DDG4/Factories.h"
0002
0003 // Framework include files
0004 #include "DD4hep/Printout.h"
0005 #include "DD4hep/InstanceCount.h"
0006 #include "DDG4/Geant4Random.h"
0007 #include "DDG4/Geant4Context.h"
0008 #include "DDG4/Geant4InputHandling.h"
0009 #include "npdet/EICInteractionVertexSmear.h"
0010
0011 // C/C++ include files
0012 #include <cmath>
0013
0014 #include "Math/Vector3D.h"
0015 #include "Math/Vector4D.h"
0016 #include "Math/LorentzRotation.h"
0017 #include "Math/AxisAngle.h"
0018 #include "Math/Rotation3D.h"
0019 #include "Math/RotationX.h"
0020 #include "Math/RotationY.h"
0021 #include "Math/Boost.h"
0022 #include "Math/BoostX.h"
0023 #include "Math/BoostX.h"
0024 #include "Math/PxPyPzM4D.h"
0025
0026 namespace npdet::sim {
0027
0028 using namespace dd4hep::sim;
0029
0030 /// Standard constructor
0031 EICInteractionVertexSmear::EICInteractionVertexSmear(Geant4Context* ctxt, const std::string& nam)
0032   : Geant4GeneratorAction(ctxt, nam)
0033 {
0034   dd4hep::InstanceCount::increment(this);
0035   declareProperty("Offset", m_offset);
0036   declareProperty("Sigma_ion",  m_sigma_Ion);
0037   declareProperty("Sigma_e",  m_sigma_Electron);
0038   declareProperty("Mask",   m_mask = 0);
0039   m_needsControl = true;
0040 }
0041
0042 /// Default destructor
0043 EICInteractionVertexSmear::~EICInteractionVertexSmear() { dd4hep::InstanceCount::decrement(this); }
0044
0045 /// Action to smear one single interaction according to the properties
0046 void EICInteractionVertexSmear::smear(Interaction* inter) const {
0047   using namespace ROOT::Math;
0048   Geant4Random& rndm = context()->event().random();
0049   if (inter) {
0050     double    dx  = rndm.gauss(m_offset.x(), m_sigma_Ion.x());
0051     double    dy  = rndm.gauss(m_offset.y(), m_sigma_Ion.y());
0052     //double    dz  = rndm.gauss(m_offset.z(), m_sigma_Ion.z());
0053     //double    dt  = rndm.gauss(m_offset.t(), m_sigma_Ion.t());
0054     //double    dxe = rndm.gauss(m_offset.x(), m_sigma_Electron.x());
0055     //double    dye = rndm.gauss(m_offset.y(), m_sigma_Electron.y());
0056     //double    dze = rndm.gauss(m_offset.z(), m_sigma_Electron.z());
0057     //double    dte = rndm.gauss(m_offset.t(), m_sigma_Electron.t());
0058     XYZVector ion_dir(0, 0, 1.0);
0059     RotationY rotx_ion(-dx);
0060     RotationX roty_ion(dy);
0061     auto      new_ion_dir = rotx_ion(roty_ion(ion_dir));
0062
0063     Geant4PrimaryEvent::Interaction::VertexMap::iterator   iv;
0064     Geant4PrimaryEvent::Interaction::ParticleMap::iterator ip;
0065
0066     double alpha         = new_ion_dir.Theta();
0067     auto   rotation_axis = ion_dir.Cross(new_ion_dir).Unit();
0068
0069     if (inter->locked) {
0070       this->abortRun("Locked interactions may not be boosted!",
0071                      "Cannot boost interactions with a native G4 primary record!");
0072     } else if (alpha != 0.0) {
0073
0074       double tanalpha  = std::tan(alpha/2.0);
0075       double gamma     = std::sqrt(1 + tanalpha * tanalpha);
0076       double betagamma = tanalpha;
0077       double beta      = betagamma / gamma;
0078
0079       Polar3D boost_vec(beta, M_PI / 2.0, new_ion_dir.Phi());
0080
0081       std::cout << " beta = " << beta << "\n";
0082
0083       Boost           bst(boost_vec);
0084       LorentzRotation roty(AxisAngle(rotation_axis, alpha / 2.0));
0085
0086       // Now move begin and end-vertex of all primary vertices accordingly
0087       for (iv = inter->vertices.begin(); iv != inter->vertices.end(); ++iv) {
0088         for (Geant4Vertex* v : (*iv).second) {
0089           XYZTVector v0(v->x, v->y, v->z, v->time);
0090           auto       v1 = roty(bst * v0);
0091           v->x          = v1.x();
0092           v->y          = v1.y();
0093           v->z          = v1.z();
0094           v->time       = v1.t();
0095         }
0096       }
0097       // Now move begin and end-vertex of all primary and generator particles accordingly
0098       for (ip = inter->particles.begin(); ip != inter->particles.end(); ++ip) {
0099         Geant4ParticleHandle p = (*ip).second;
0100         XYZTVector           v0{p->vsx, p->vsy, p->vsz, p->time};
0101         auto                 v1 = roty(bst * v0);
0102
0103         PxPyPzMVector p0{p->psx, p->psy, p->psz, p->mass};
0104         auto          p1 = roty(bst * p0);
0105
0106         p->vsx  = v1.x();
0107         p->vsy  = v1.y();
0108         p->vsz  = v1.z();
0109         p->time = v1.t();
0110
0111         p->psx = p1.x();
0112         p->psy = p1.y();
0113         p->psz = p1.z();
0114       }
0115     }
0116
0117     // smearInteraction(this,inter,dx,dy,dz,dt);
0118     return;
0119   } else {
0120     //print("+++ Smearing primary vertex for interaction type %d (%d Vertices, %d particles) "
0121     //      "by (%+.2e mm, %+.2e mm, %+.2e mm, %+.2e ns)",
0122     //      m_mask, int(inter->vertices.size()), int(inter->particles.size()), dx, dy, dz, dt);
0123     print("+++ No interaction of type %d present.", m_mask);
0124   }
0125 }
0126
0127 /// Callback to generate primary particles
0128 void EICInteractionVertexSmear::operator()(G4Event*) {
0129   typedef std::vector<Geant4PrimaryInteraction*> _I;
0130   Geant4PrimaryEvent* evt = context()->event().extension<Geant4PrimaryEvent>();
0131
0132   if ( m_mask >= 0 )  {
0133     Interaction* inter = evt->get(m_mask);
0134     smear(inter);
0135     return;
0136   }
0137   _I interactions = evt->interactions();
0138   for(_I::iterator i=interactions.begin(); i != interactions.end(); ++i)
0139     smear(*i);
0140 }
0141 }
0142
0143 namespace dd4hep::sim {
0144 using EICInteractionVertexSmear = npdet::sim::EICInteractionVertexSmear;
0145 }
0146 DECLARE_GEANT4ACTION(EICInteractionVertexSmear)
0147