EventVertexGenerators/src/BetafuncEvtVtxGenerator.cc

0001 /*
0002 ________________________________________________________________________
0003
0004  BetafuncEvtVtxGenerator
0005
0006  Smear vertex according to the Beta function on the transverse plane
0007  and a Gaussian on the z axis. It allows the beam to have a crossing
0008  angle (slopes dxdz and dydz).
0009
0010  Based on GaussEvtVtxGenerator
0011  implemented by Francisco Yumiceva (yumiceva@fnal.gov)
0012
0013  FERMILAB
0014  2006
0015 ________________________________________________________________________
0016 */
0017
0018 #include "IOMC/EventVertexGenerators/interface/BetafuncEvtVtxGenerator.h"
0019 #include "FWCore/Utilities/interface/Exception.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/Framework/interface/EventSetup.h"
0022 #include "FWCore/Framework/interface/ESHandle.h"
0023
0024 #include <CLHEP/Random/RandGaussQ.h>
0025 #include <CLHEP/Units/SystemOfUnits.h>
0026 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0027
0028 using CLHEP::cm;
0029 using CLHEP::ns;
0030 using CLHEP::radian;
0031
0032 BetafuncEvtVtxGenerator::BetafuncEvtVtxGenerator(const edm::ParameterSet& p) : BaseEvtVtxGenerator(p), boost_(4, 4) {
0033   readDB_ = p.getParameter<bool>("readDB");
0034   if (!readDB_) {
0035     fX0 = p.getParameter<double>("X0") * cm;
0036     fY0 = p.getParameter<double>("Y0") * cm;
0037     fZ0 = p.getParameter<double>("Z0") * cm;
0038     fSigmaZ = p.getParameter<double>("SigmaZ") * cm;
0039     fbetastar = p.getParameter<double>("BetaStar") * cm;
0040     femittance = p.getParameter<double>("Emittance") * cm;              // this is not the normalized emittance
0041     fTimeOffset = p.getParameter<double>("TimeOffset") * ns * c_light;  // HepMC distance units are in mm
0042
0043     setBoost(p.getParameter<double>("Alpha") * radian, p.getParameter<double>("Phi") * radian);
0044     if (fSigmaZ <= 0) {
0045       throw cms::Exception("Configuration") << "Error in BetafuncEvtVtxGenerator: "
0046                                             << "Illegal resolution in Z (SigmaZ is negative)";
0047     }
0048   }
0049   if (readDB_) {
0050     // NOTE: this is currently watching LS transitions, while it should watch Run transitions,
0051     // even though in reality there is no Run Dependent MC (yet) in CMS
0052     beamToken_ = esConsumes<SimBeamSpotObjects, SimBeamSpotObjectsRcd, edm::Transition::BeginLuminosityBlock>();
0053   }
0054 }
0055
0056 void BetafuncEvtVtxGenerator::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const& iEventSetup) {
0057   update(iEventSetup);
0058 }
0059
0060 void BetafuncEvtVtxGenerator::update(const edm::EventSetup& iEventSetup) {
0061   if (readDB_ && parameterWatcher_.check(iEventSetup)) {
0062     edm::ESHandle<SimBeamSpotObjects> beamhandle = iEventSetup.getHandle(beamToken_);
0063     if (!beamhandle->isGaussian()) {
0064       fX0 = beamhandle->x() * cm;
0065       fY0 = beamhandle->y() * cm;
0066       fZ0 = beamhandle->z() * cm;
0067       fSigmaZ = beamhandle->sigmaZ() * cm;
0068       fTimeOffset = beamhandle->timeOffset() * ns * c_light;  // HepMC distance units are in mm
0069       fbetastar = beamhandle->betaStar() * cm;
0070       femittance = beamhandle->emittance() * cm;
0071       setBoost(beamhandle->alpha() * radian, beamhandle->phi() * radian);
0072     } else {
0073       throw cms::Exception("Configuration")
0074           << "Error in BetafuncEvtVtxGenerator::update: The provided SimBeamSpotObjects is Gaussian.\n"
0075           << "Please check the configuration and ensure that the beam spot parameters are appropriate for a Betafunc "
0076              "distribution.";
0077     }
0078   }
0079 }
0080
0081 ROOT::Math::XYZTVector BetafuncEvtVtxGenerator::vertexShift(CLHEP::HepRandomEngine* engine) const {
0082   double X, Y, Z;
0083
0084   double tmp_sigz = CLHEP::RandGaussQ::shoot(engine, 0., fSigmaZ);
0085   Z = tmp_sigz + fZ0;
0086
0087   double tmp_sigx = BetaFunction(Z, fZ0);
0088   // need sqrt(2) for beamspot width relative to single beam width
0089   tmp_sigx /= sqrt(2.0);
0090   X = CLHEP::RandGaussQ::shoot(engine, 0., tmp_sigx) + fX0;  // + Z*fdxdz ;
0091
0092   double tmp_sigy = BetaFunction(Z, fZ0);
0093   // need sqrt(2) for beamspot width relative to single beam width
0094   tmp_sigy /= sqrt(2.0);
0095   Y = CLHEP::RandGaussQ::shoot(engine, 0., tmp_sigy) + fY0;  // + Z*fdydz;
0096
0097   double tmp_sigt = CLHEP::RandGaussQ::shoot(engine, 0., fSigmaZ);
0098   double T = tmp_sigt + fTimeOffset;
0099
0100   return ROOT::Math::XYZTVector(X, Y, Z, T);
0101 }
0102
0103 double BetafuncEvtVtxGenerator::BetaFunction(double z, double z0) const {
0104   return sqrt(femittance * (fbetastar + (((z - z0) * (z - z0)) / fbetastar)));
0105 }
0106
0107 void BetafuncEvtVtxGenerator::setBoost(double alpha, double phi) {
0108   //boost_.ResizeTo(4,4);
0109   //boost_ = new TMatrixD(4,4);
0110   TMatrixD tmpboost(4, 4);
0111
0112   //if ( (alpha_ == 0) && (phi_==0) ) { boost_->Zero(); return boost_; }
0113
0114   // Lorentz boost to frame where the collision is head-on
0115   // phi is the half crossing angle in the plane ZS
0116   // alpha is the angle to the S axis from the X axis in the XY plane
0117
0118   tmpboost(0, 0) = 1. / cos(phi);
0119   tmpboost(0, 1) = -cos(alpha) * sin(phi);
0120   tmpboost(0, 2) = -tan(phi) * sin(phi);
0121   tmpboost(0, 3) = -sin(alpha) * sin(phi);
0122   tmpboost(1, 0) = -cos(alpha) * tan(phi);
0123   tmpboost(1, 1) = 1.;
0124   tmpboost(1, 2) = cos(alpha) * tan(phi);
0125   tmpboost(1, 3) = 0.;
0126   tmpboost(2, 0) = 0.;
0127   tmpboost(2, 1) = -cos(alpha) * sin(phi);
0128   tmpboost(2, 2) = cos(phi);
0129   tmpboost(2, 3) = -sin(alpha) * sin(phi);
0130   tmpboost(3, 0) = -sin(alpha) * tan(phi);
0131   tmpboost(3, 1) = 0.;
0132   tmpboost(3, 2) = sin(alpha) * tan(phi);
0133   tmpboost(3, 3) = 1.;
0134
0135   tmpboost.Invert();
0136   boost_ = tmpboost;
0137   //boost_->Print();
0138 }
0139
0140 void BetafuncEvtVtxGenerator::sigmaZ(double s) {
0141   if (s >= 0) {
0142     fSigmaZ = s;
0143   } else {
0144     throw cms::Exception("LogicError") << "Error in BetafuncEvtVtxGenerator::sigmaZ: "
0145                                        << "Illegal resolution in Z (negative)";
0146   }
0147 }
0148
0149 TMatrixD const* BetafuncEvtVtxGenerator::GetInvLorentzBoost() const { return &boost_; }
0150
0151 void BetafuncEvtVtxGenerator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0152   edm::ParameterSetDescription desc;
0153   desc.add<double>("X0", 0.0)->setComment("in cm");
0154   desc.add<double>("Y0", 0.0)->setComment("in cm");
0155   desc.add<double>("Z0", 0.0)->setComment("in cm");
0156   desc.add<double>("SigmaZ", 0.0)->setComment("in cm");
0157   desc.add<double>("BetaStar", 0.0)->setComment("in cm");
0158   desc.add<double>("Emittance", 0.0)->setComment("in cm");
0159   desc.add<double>("Alpha", 0.0)->setComment("in radians");
0160   desc.add<double>("Phi", 0.0)->setComment("in radians");
0161   desc.add<double>("TimeOffset", 0.0)->setComment("in ns");
0162   desc.add<edm::InputTag>("src");
0163   desc.add<bool>("readDB");
0164   descriptions.add("BetafuncEvtVtxGenerator", desc);
0165 }