EICInteractionVertexBoost.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/Geant4InputHandling.h"
#include "npdet/EICInteractionVertexBoost.h"
 
#include "Math/Vector4D.h"
#include "Math/LorentzRotation.h"
#include "Math/Boost.h"
#include "Math/BoostX.h"
#include "Math/BoostX.h"
#include "Math/PxPyPzM4D.h"
 
 
namespace npdet::sim {
 
using namespace dd4hep::sim;
 
EICInteractionVertexBoost::EICInteractionVertexBoost(Geant4Context* ctxt, const std::string& nam)
 : Geant4GeneratorAction(ctxt, nam)
{
 dd4hep::InstanceCount::increment(this);
 //declareProperty("Angle", m_angle = 0);
 declareProperty("IonCrossingAngle", m_ionCrossingAngle = 0.0166667);
 declareProperty("ElectronCrossingAngle", m_eCrossingAngle = 0.00833333);
 declareProperty("Mask", m_mask = 1);
 m_needsControl = true;
}
 
EICInteractionVertexBoost::~EICInteractionVertexBoost() {
 dd4hep::InstanceCount::decrement(this);
}
 
void EICInteractionVertexBoost::boost(Interaction* inter) const {
 using namespace ROOT::Math;
 // boostInteraction(this, inter, m_angle);
 // int dd4hep::sim::boostInteraction(const Geant4Action* caller,
 // Geant4PrimaryEvent::Interaction* inter,
 // double alpha)
 if (inter) {
 Geant4PrimaryEvent::Interaction::VertexMap::iterator iv;
 Geant4PrimaryEvent::Interaction::ParticleMap::iterator ip;
 
 double alpha = m_ionCrossingAngle+m_eCrossingAngle;
 
 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;
 
 std::cout << " beta = " << beta << "\n";
 std::cout << " gamma = " << gamma << "\n";
 std::cout << " alpha = " << alpha*180.0/M_PI << " deg\n";
 
 BoostX bst(beta);
 LorentzRotation roty(RotationY(alpha/2.0-m_eCrossingAngle));
 
 // 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);
 //double t = gamma * v->time + betagamma * v->x / CLHEP::c_light;
 //double x = gamma * v->x + betagamma * CLHEP::c_light * v->time;
 //double y = v->y;
 //double z = v->z;
 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;
 //double t = gamma * p->time + betagamma * p->vsx / CLHEP::c_light;
 //double x = gamma * p->vsx + betagamma * CLHEP::c_light * p->time;
 //double y = p->vsy;
 //double z = p->vsz;
 
 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);
 //double m = p->mass;
 //double e2 = SQR(p->psx) + SQR(p->psy) + SQR(p->psz) + SQR(m);
 //double px = betagamma * std::sqrt(e2) + gamma * p->psx;
 //double py = p->psy;
 //double pz = p->psz;
 
 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();
 
 std::cout << p1 << "\n";
 std::cout << " theta = " << p1.Theta() << " \n";
 }
 }
 }
 else {
 print("+++ No interaction of mask/type %d present.", m_mask);
 }
}
 
void EICInteractionVertexBoost::operator()(G4Event*) {
 typedef std::vector<Geant4PrimaryInteraction*> _I;
 Geant4PrimaryEvent* evt = context()->event().extension<Geant4PrimaryEvent>();
 
 if ( m_mask >= 0 ) {
 std::cout << " mask is not zero\n";
 Interaction* inter = evt->get(m_mask);
 boost(inter);
 return;
 }
 _I interactions = evt->interactions();
 for(_I::iterator i=interactions.begin(); i != interactions.end(); ++i)
 boost(*i);
}
 
} // namespace npdet::sim
 
namespace dd4hep::sim {
using EICInteractionVertexBoost = npdet::sim::EICInteractionVertexBoost;
}
 
DECLARE_GEANT4ACTION(EICInteractionVertexBoost)