File indexing completed on 2023-03-17 11:10:13
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019 #include "IOMC/EventVertexGenerators/interface/BetafuncEvtVtxGenerator.h"
0020 #include "FWCore/Utilities/interface/Exception.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022 #include "FWCore/Framework/interface/EventSetup.h"
0023 #include "FWCore/Framework/interface/ESHandle.h"
0024
0025 #include "CLHEP/Random/RandGaussQ.h"
0026 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0027 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0028
0029 #include "HepMC/SimpleVector.h"
0030
0031 #include <iostream>
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;
0042 fTimeOffset = p.getParameter<double>("TimeOffset") * ns * c_light;
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 beamToken_ = esConsumes<SimBeamSpotObjects, SimBeamSpotObjectsRcd, edm::Transition::BeginLuminosityBlock>();
0052 }
0053 }
0054
0055 BetafuncEvtVtxGenerator::~BetafuncEvtVtxGenerator() {}
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
0065 fX0 = beamhandle->fX0;
0066 fY0 = beamhandle->fY0;
0067 fZ0 = beamhandle->fZ0;
0068
0069 fSigmaZ = beamhandle->fSigmaZ;
0070 fTimeOffset = beamhandle->fTimeOffset;
0071 fbetastar = beamhandle->fbetastar;
0072 femittance = beamhandle->femittance;
0073 setBoost(beamhandle->fAlpha, beamhandle->fPhi);
0074 }
0075 }
0076
0077 HepMC::FourVector BetafuncEvtVtxGenerator::newVertex(CLHEP::HepRandomEngine* engine) const {
0078 double X, Y, Z;
0079
0080 double tmp_sigz = CLHEP::RandGaussQ::shoot(engine, 0., fSigmaZ);
0081 Z = tmp_sigz + fZ0;
0082
0083 double tmp_sigx = BetaFunction(Z, fZ0);
0084
0085 tmp_sigx /= sqrt(2.0);
0086 X = CLHEP::RandGaussQ::shoot(engine, 0., tmp_sigx) + fX0;
0087
0088 double tmp_sigy = BetaFunction(Z, fZ0);
0089
0090 tmp_sigy /= sqrt(2.0);
0091 Y = CLHEP::RandGaussQ::shoot(engine, 0., tmp_sigy) + fY0;
0092
0093 double tmp_sigt = CLHEP::RandGaussQ::shoot(engine, 0., fSigmaZ);
0094 double T = tmp_sigt + fTimeOffset;
0095
0096 return HepMC::FourVector(X, Y, Z, T);
0097 }
0098
0099 double BetafuncEvtVtxGenerator::BetaFunction(double z, double z0) const {
0100 return sqrt(femittance * (fbetastar + (((z - z0) * (z - z0)) / fbetastar)));
0101 }
0102
0103 void BetafuncEvtVtxGenerator::setBoost(double alpha, double phi) {
0104
0105
0106 TMatrixD tmpboost(4, 4);
0107
0108
0109
0110
0111
0112
0113
0114 tmpboost(0, 0) = 1. / cos(phi);
0115 tmpboost(0, 1) = -cos(alpha) * sin(phi);
0116 tmpboost(0, 2) = -tan(phi) * sin(phi);
0117 tmpboost(0, 3) = -sin(alpha) * sin(phi);
0118 tmpboost(1, 0) = -cos(alpha) * tan(phi);
0119 tmpboost(1, 1) = 1.;
0120 tmpboost(1, 2) = cos(alpha) * tan(phi);
0121 tmpboost(1, 3) = 0.;
0122 tmpboost(2, 0) = 0.;
0123 tmpboost(2, 1) = -cos(alpha) * sin(phi);
0124 tmpboost(2, 2) = cos(phi);
0125 tmpboost(2, 3) = -sin(alpha) * sin(phi);
0126 tmpboost(3, 0) = -sin(alpha) * tan(phi);
0127 tmpboost(3, 1) = 0.;
0128 tmpboost(3, 2) = sin(alpha) * tan(phi);
0129 tmpboost(3, 3) = 1.;
0130
0131 tmpboost.Invert();
0132 boost_ = tmpboost;
0133
0134 }
0135
0136 void BetafuncEvtVtxGenerator::sigmaZ(double s) {
0137 if (s >= 0) {
0138 fSigmaZ = s;
0139 } else {
0140 throw cms::Exception("LogicError") << "Error in BetafuncEvtVtxGenerator::sigmaZ: "
0141 << "Illegal resolution in Z (negative)";
0142 }
0143 }
0144
0145 TMatrixD const* BetafuncEvtVtxGenerator::GetInvLorentzBoost() const { return &boost_; }