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 #include "HepMC/SimpleVector.h"
0028
0029 using CLHEP::cm;
0030 using CLHEP::ns;
0031 using CLHEP::radian;
0032
0033 BetafuncEvtVtxGenerator::BetafuncEvtVtxGenerator(const edm::ParameterSet& p) : BaseEvtVtxGenerator(p), boost_(4, 4) {
0034   readDB_ = p.getParameter<bool>("readDB");
0035   if (!readDB_) {
0036     fX0 = p.getParameter<double>("X0") * cm;
0037     fY0 = p.getParameter<double>("Y0") * cm;
0038     fZ0 = p.getParameter<double>("Z0") * cm;
0039     fSigmaZ = p.getParameter<double>("SigmaZ") * cm;
0040     fbetastar = p.getParameter<double>("BetaStar") * cm;
0041     femittance = p.getParameter<double>("Emittance") * cm;              // this is not the normalized emittance
0042     fTimeOffset = p.getParameter<double>("TimeOffset") * ns * c_light;  // HepMC distance units are in mm
0043
0044     setBoost(p.getParameter<double>("Alpha") * radian, p.getParameter<double>("Phi") * radian);
0045     if (fSigmaZ <= 0) {
0046       throw cms::Exception("Configuration") << "Error in BetafuncEvtVtxGenerator: "
0047                                             << "Illegal resolution in Z (SigmaZ is negative)";
0048     }
0049   }
0050   if (readDB_) {
0051     // NOTE: this is currently watching LS transitions, while it should watch Run transitions,
0052     // even though in reality there is no Run Dependent MC (yet) in CMS
0053     beamToken_ = esConsumes<SimBeamSpotObjects, SimBeamSpotObjectsRcd, edm::Transition::BeginLuminosityBlock>();
0054   }
0055 }
0056
0057 void BetafuncEvtVtxGenerator::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const& iEventSetup) {
0058   update(iEventSetup);
0059 }
0060
0061 void BetafuncEvtVtxGenerator::update(const edm::EventSetup& iEventSetup) {
0062   if (readDB_ && parameterWatcher_.check(iEventSetup)) {
0063     edm::ESHandle<SimBeamSpotObjects> beamhandle = iEventSetup.getHandle(beamToken_);
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   }
0073 }
0074
0075 HepMC::FourVector BetafuncEvtVtxGenerator::newVertex(CLHEP::HepRandomEngine* engine) const {
0076   double X, Y, Z;
0077
0078   double tmp_sigz = CLHEP::RandGaussQ::shoot(engine, 0., fSigmaZ);
0079   Z = tmp_sigz + fZ0;
0080
0081   double tmp_sigx = BetaFunction(Z, fZ0);
0082   // need sqrt(2) for beamspot width relative to single beam width
0083   tmp_sigx /= sqrt(2.0);
0084   X = CLHEP::RandGaussQ::shoot(engine, 0., tmp_sigx) + fX0;  // + Z*fdxdz ;
0085
0086   double tmp_sigy = BetaFunction(Z, fZ0);
0087   // need sqrt(2) for beamspot width relative to single beam width
0088   tmp_sigy /= sqrt(2.0);
0089   Y = CLHEP::RandGaussQ::shoot(engine, 0., tmp_sigy) + fY0;  // + Z*fdydz;
0090
0091   double tmp_sigt = CLHEP::RandGaussQ::shoot(engine, 0., fSigmaZ);
0092   double T = tmp_sigt + fTimeOffset;
0093
0094   return HepMC::FourVector(X, Y, Z, T);
0095 }
0096
0097 double BetafuncEvtVtxGenerator::BetaFunction(double z, double z0) const {
0098   return sqrt(femittance * (fbetastar + (((z - z0) * (z - z0)) / fbetastar)));
0099 }
0100
0101 void BetafuncEvtVtxGenerator::setBoost(double alpha, double phi) {
0102   //boost_.ResizeTo(4,4);
0103   //boost_ = new TMatrixD(4,4);
0104   TMatrixD tmpboost(4, 4);
0105
0106   //if ( (alpha_ == 0) && (phi_==0) ) { boost_->Zero(); return boost_; }
0107
0108   // Lorentz boost to frame where the collision is head-on
0109   // phi is the half crossing angle in the plane ZS
0110   // alpha is the angle to the S axis from the X axis in the XY plane
0111
0112   tmpboost(0, 0) = 1. / cos(phi);
0113   tmpboost(0, 1) = -cos(alpha) * sin(phi);
0114   tmpboost(0, 2) = -tan(phi) * sin(phi);
0115   tmpboost(0, 3) = -sin(alpha) * sin(phi);
0116   tmpboost(1, 0) = -cos(alpha) * tan(phi);
0117   tmpboost(1, 1) = 1.;
0118   tmpboost(1, 2) = cos(alpha) * tan(phi);
0119   tmpboost(1, 3) = 0.;
0120   tmpboost(2, 0) = 0.;
0121   tmpboost(2, 1) = -cos(alpha) * sin(phi);
0122   tmpboost(2, 2) = cos(phi);
0123   tmpboost(2, 3) = -sin(alpha) * sin(phi);
0124   tmpboost(3, 0) = -sin(alpha) * tan(phi);
0125   tmpboost(3, 1) = 0.;
0126   tmpboost(3, 2) = sin(alpha) * tan(phi);
0127   tmpboost(3, 3) = 1.;
0128
0129   tmpboost.Invert();
0130   boost_ = tmpboost;
0131   //boost_->Print();
0132 }
0133
0134 void BetafuncEvtVtxGenerator::sigmaZ(double s) {
0135   if (s >= 0) {
0136     fSigmaZ = s;
0137   } else {
0138     throw cms::Exception("LogicError") << "Error in BetafuncEvtVtxGenerator::sigmaZ: "
0139                                        << "Illegal resolution in Z (negative)";
0140   }
0141 }
0142
0143 TMatrixD const* BetafuncEvtVtxGenerator::GetInvLorentzBoost() const { return &boost_; }
0144
0145 void BetafuncEvtVtxGenerator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0146   edm::ParameterSetDescription desc;
0147   desc.add<double>("X0", 0.0)->setComment("in cm");
0148   desc.add<double>("Y0", 0.0)->setComment("in cm");
0149   desc.add<double>("Z0", 0.0)->setComment("in cm");
0150   desc.add<double>("SigmaZ", 0.0)->setComment("in cm");
0151   desc.add<double>("BetaStar", 0.0)->setComment("in cm");
0152   desc.add<double>("Emittance", 0.0)->setComment("in cm");
0153   desc.add<double>("Alpha", 0.0)->setComment("in radians");
0154   desc.add<double>("Phi", 0.0)->setComment("in radians");
0155   desc.add<double>("TimeOffset", 0.0)->setComment("in ns");
0156   desc.add<edm::InputTag>("src");
0157   desc.add<bool>("readDB");
0158   descriptions.add("BetafuncEvtVtxGenerator", desc);
0159 }