HiGenCommon/plugins/BetaBoostEvtVtxGenerator.cc

0001
0002 /*
0003   ________________________________________________________________________
0004
0005   BetaBoostEvtVtxGenerator
0006
0007   Smear vertex according to the Beta function on the transverse plane
0008   and a Gaussian on the z axis. It allows the beam to have a crossing
0009   angle (slopes dxdz and dydz).
0010
0011   Based on GaussEvtVtxGenerator
0012   implemented by Francisco Yumiceva (yumiceva@fnal.gov)
0013
0014   FERMILAB
0015   2006
0016   ________________________________________________________________________
0017 */
0018
0019 //lingshan: add beta for z-axis boost
0020
0021 //#include "IOMC/EventVertexGenerators/interface/BetafuncEvtVtxGenerator.h"
0022
0023 #include "FWCore/Framework/interface/Event.h"
0024 #include "FWCore/Utilities/interface/Exception.h"
0025 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0026
0027 #include "FWCore/Framework/interface/global/EDProducer.h"
0028 #include "FWCore/Utilities/interface/InputTag.h"
0029 #include "FWCore/ServiceRegistry/interface/Service.h"
0030 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0031 #include "FWCore/Utilities/interface/EDGetToken.h"
0032 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0033
0034 #include <CLHEP/Random/RandGaussQ.h>
0035 #include <CLHEP/Units/SystemOfUnits.h>
0036 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0037 //#include "CLHEP/Vector/ThreeVector.h"
0038 #include "HepMC/SimpleVector.h"
0039 #include "TMatrixD.h"
0040
0041 #include <iostream>
0042
0043 using namespace edm;
0044 using namespace std;
0045 using namespace CLHEP;
0046
0047 namespace CLHEP {
0048   class HepRandomEngine;
0049 }
0050
0051 class BetaBoostEvtVtxGenerator : public edm::global::EDProducer<> {
0052 public:
0053   BetaBoostEvtVtxGenerator(const edm::ParameterSet& p);
0054   /** Copy constructor */
0055   BetaBoostEvtVtxGenerator(const BetaBoostEvtVtxGenerator& p) = delete;
0056   /** Copy assignment operator */
0057   BetaBoostEvtVtxGenerator& operator=(const BetaBoostEvtVtxGenerator& rhs) = delete;
0058   ~BetaBoostEvtVtxGenerator() override = default;
0059
0060   /// return a new event vertex
0061   HepMC::FourVector newVertex(CLHEP::HepRandomEngine*) const;
0062   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0063   TMatrixD GetInvLorentzBoost() const;
0064
0065   /// beta function
0066   double BetaFunction(double z, double z0) const;
0067
0068 private:
0069   const double alpha_;  // angle between crossing plane and horizontal plane
0070   const double phi_;    // half crossing angle
0071   const double beta_;
0072   const double fX0;      // mean in X in cm
0073   const double fY0;      // mean in Y in cm
0074   const double fZ0;      // mean in Z in cm
0075   const double fSigmaZ;  // resolution in Z in cm
0076   //double fdxdz, fdydz;
0077   const double fbetastar;
0078   const double femittance;  // emittance (no the normalized)a
0079
0080   const double fTimeOffset;
0081
0082   const TMatrixD boost_;
0083
0084   const edm::EDGetTokenT<HepMCProduct> sourceLabel;
0085
0086   const bool verbosity_;
0087 };
0088
0089 BetaBoostEvtVtxGenerator::BetaBoostEvtVtxGenerator(const edm::ParameterSet& p)
0090     : alpha_(p.getParameter<double>("Alpha") * radian),
0091       phi_(p.getParameter<double>("Phi") * radian),
0092       beta_(p.getParameter<double>("Beta")),
0093       fX0(p.getParameter<double>("X0") * cm),
0094       fY0(p.getParameter<double>("Y0") * cm),
0095       fZ0(p.getParameter<double>("Z0") * cm),
0096       fSigmaZ(p.getParameter<double>("SigmaZ") * cm),
0097       fbetastar(p.getParameter<double>("BetaStar") * cm),
0098       femittance(p.getParameter<double>("Emittance") * cm),              // this is not the normalized emittance
0099       fTimeOffset(p.getParameter<double>("TimeOffset") * ns * c_light),  // HepMC time units are mm
0100       boost_(GetInvLorentzBoost()),
0101       sourceLabel(consumes<HepMCProduct>(p.getParameter<edm::InputTag>("src"))),
0102       verbosity_(p.getUntrackedParameter<bool>("verbosity", false)) {
0103   if (fSigmaZ <= 0) {
0104     throw cms::Exception("Configuration") << "Error in BetaBoostEvtVtxGenerator: "
0105                                           << "Illegal resolution in Z (SigmaZ is negative)";
0106   }
0107
0108   produces<edm::HepMCProduct>();
0109 }
0110
0111 HepMC::FourVector BetaBoostEvtVtxGenerator::newVertex(CLHEP::HepRandomEngine* engine) const {
0112   double X, Y, Z;
0113
0114   double tmp_sigz = CLHEP::RandGaussQ::shoot(engine, 0.0, fSigmaZ);
0115   Z = tmp_sigz + fZ0;
0116
0117   double tmp_sigx = BetaFunction(Z, fZ0);
0118   // need sqrt(2) for beamspot width relative to single beam width
0119   tmp_sigx /= sqrt(2.0);
0120   X = CLHEP::RandGaussQ::shoot(engine, 0.0, tmp_sigx) + fX0;  // + Z*fdxdz ;
0121
0122   double tmp_sigy = BetaFunction(Z, fZ0);
0123   // need sqrt(2) for beamspot width relative to single beam width
0124   tmp_sigy /= sqrt(2.0);
0125   Y = CLHEP::RandGaussQ::shoot(engine, 0.0, tmp_sigy) + fY0;  // + Z*fdydz;
0126
0127   double tmp_sigt = CLHEP::RandGaussQ::shoot(engine, 0.0, fSigmaZ);
0128   double T = tmp_sigt + fTimeOffset;
0129
0130   return HepMC::FourVector(X, Y, Z, T);
0131 }
0132
0133 double BetaBoostEvtVtxGenerator::BetaFunction(double z, double z0) const {
0134   return sqrt(femittance * (fbetastar + (((z - z0) * (z - z0)) / fbetastar)));
0135 }
0136
0137 TMatrixD BetaBoostEvtVtxGenerator::GetInvLorentzBoost() const {
0138   //alpha_ = 0;
0139   //phi_ = 142.e-6;
0140   //    if (boost_ != 0 ) return boost_;
0141
0142   //boost_.ResizeTo(4,4);
0143   //boost_ = new TMatrixD(4,4);
0144   TMatrixD tmpboost(4, 4);
0145   TMatrixD tmpboostZ(4, 4);
0146   TMatrixD tmpboostXYZ(4, 4);
0147
0148   //if ( (alpha_ == 0) && (phi_==0) ) { boost_->Zero(); return boost_; }
0149
0150   // Lorentz boost to frame where the collision is head-on
0151   // phi is the half crossing angle in the plane ZS
0152   // alpha is the angle to the S axis from the X axis in the XY plane
0153
0154   tmpboost(0, 0) = 1. / cos(phi_);
0155   tmpboost(0, 1) = -cos(alpha_) * sin(phi_);
0156   tmpboost(0, 2) = -tan(phi_) * sin(phi_);
0157   tmpboost(0, 3) = -sin(alpha_) * sin(phi_);
0158   tmpboost(1, 0) = -cos(alpha_) * tan(phi_);
0159   tmpboost(1, 1) = 1.;
0160   tmpboost(1, 2) = cos(alpha_) * tan(phi_);
0161   tmpboost(1, 3) = 0.;
0162   tmpboost(2, 0) = 0.;
0163   tmpboost(2, 1) = -cos(alpha_) * sin(phi_);
0164   tmpboost(2, 2) = cos(phi_);
0165   tmpboost(2, 3) = -sin(alpha_) * sin(phi_);
0166   tmpboost(3, 0) = -sin(alpha_) * tan(phi_);
0167   tmpboost(3, 1) = 0.;
0168   tmpboost(3, 2) = sin(alpha_) * tan(phi_);
0169   tmpboost(3, 3) = 1.;
0170   //cout<<"beta "<<beta_;
0171   double gama = 1.0 / sqrt(1 - beta_ * beta_);
0172   tmpboostZ(0, 0) = gama;
0173   tmpboostZ(0, 1) = 0.;
0174   tmpboostZ(0, 2) = -1.0 * beta_ * gama;
0175   tmpboostZ(0, 3) = 0.;
0176   tmpboostZ(1, 0) = 0.;
0177   tmpboostZ(1, 1) = 1.;
0178   tmpboostZ(1, 2) = 0.;
0179   tmpboostZ(1, 3) = 0.;
0180   tmpboostZ(2, 0) = -1.0 * beta_ * gama;
0181   tmpboostZ(2, 1) = 0.;
0182   tmpboostZ(2, 2) = gama;
0183   tmpboostZ(2, 3) = 0.;
0184   tmpboostZ(3, 0) = 0.;
0185   tmpboostZ(3, 1) = 0.;
0186   tmpboostZ(3, 2) = 0.;
0187   tmpboostZ(3, 3) = 1.;
0188
0189   tmpboostXYZ = tmpboostZ * tmpboost;
0190   tmpboostXYZ.Invert();
0191
0192   if (verbosity_) {
0193     tmpboostXYZ.Print();
0194   }
0195
0196   return tmpboostXYZ;
0197 }
0198
0199 void BetaBoostEvtVtxGenerator::produce(edm::StreamID, Event& evt, const EventSetup&) const {
0200   edm::Service<RandomNumberGenerator> rng;
0201   if (!rng.isAvailable()) {
0202     throw cms::Exception("Configuration")
0203         << "Attempt to get a random engine when the RandomNumberGeneratorService is not configured.\n"
0204            "You must configure the service if you want an engine.\n";
0205   }
0206   CLHEP::HepRandomEngine* engine = &rng->getEngine(evt.streamID());
0207
0208   const auto& HepUnsmearedMCEvt = evt.get(sourceLabel);
0209
0210   // Copy the HepMC::GenEvent
0211   auto HepMCEvt = std::make_unique<edm::HepMCProduct>(new HepMC::GenEvent(*HepUnsmearedMCEvt.GetEvent()));
0212
0213   // generate new vertex & apply the shift
0214   //
0215   auto vertex = newVertex(engine);
0216   HepMCEvt->applyVtxGen(&vertex);
0217
0218   //HepMCEvt->LorentzBoost( 0., 142.e-6 );
0219   HepMCEvt->boostToLab(&boost_, "vertex");
0220   HepMCEvt->boostToLab(&boost_, "momentum");
0221   evt.put(std::move(HepMCEvt));
0222 }
0223
0224 #include "FWCore/Framework/interface/MakerMacros.h"
0225 DEFINE_FWK_MODULE(BetaBoostEvtVtxGenerator);