File indexing completed on 2023-10-25 09:53:37
0001
0002
0003 #include "IOMC/EventVertexGenerators/interface/BeamProfileVtxGenerator.h"
0004 #include "FWCore/Utilities/interface/Exception.h"
0005
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008
0009 #include "CLHEP/Geometry/Transform3D.h"
0010 #include "CLHEP/Random/RandFlat.h"
0011 #include "CLHEP/Random/RandGaussQ.h"
0012 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0013 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0014 #include "HepMC/SimpleVector.h"
0015
0016 #include <fstream>
0017 #include <string>
0018
0019 BeamProfileVtxGenerator::BeamProfileVtxGenerator(const edm::ParameterSet& p) : BaseEvtVtxGenerator(p) {
0020 meanX(p.getParameter<double>("BeamMeanX") * cm);
0021 meanY(p.getParameter<double>("BeamMeanY") * cm);
0022 beamPos(p.getParameter<double>("BeamPosition") * cm);
0023 sigmaX(p.getParameter<double>("BeamSigmaX") * cm);
0024 sigmaY(p.getParameter<double>("BeamSigmaY") * cm);
0025 double fMinEta = p.getParameter<double>("MinEta");
0026 double fMaxEta = p.getParameter<double>("MaxEta");
0027 double fMinPhi = p.getParameter<double>("MinPhi");
0028 double fMaxPhi = p.getParameter<double>("MaxPhi");
0029 eta(0.5 * (fMinEta + fMaxEta));
0030 phi(0.5 * (fMinPhi + fMaxPhi));
0031 psi(p.getParameter<double>("Psi"));
0032 nBinx = p.getParameter<int>("BinX");
0033 nBiny = p.getParameter<int>("BinY");
0034 ffile = p.getParameter<bool>("UseFile");
0035 fTimeOffset = p.getParameter<double>("TimeOffset") * ns * c_light;
0036
0037 if (ffile) {
0038 std::string file = p.getParameter<std::string>("File");
0039 std::ifstream is(file.c_str(), std::ios::in);
0040 if (is) {
0041 double elem, sum = 0;
0042 while (!is.eof()) {
0043 is >> elem;
0044 fdistn.push_back(elem);
0045 sum += elem;
0046 }
0047 if (((int)(fdistn.size())) >= nBinx * nBiny) {
0048 double last = 0;
0049 for (unsigned int i = 0; i < fdistn.size(); i++) {
0050 fdistn[i] /= sum;
0051 fdistn[i] += last;
0052 last = fdistn[i];
0053 }
0054 setType(false);
0055 } else {
0056 ffile = false;
0057 }
0058 } else {
0059 ffile = false;
0060 }
0061 }
0062 if (!ffile) {
0063 setType(p.getParameter<bool>("GaussianProfile"));
0064 }
0065
0066 edm::LogInfo("VertexGenerator") << "BeamProfileVtxGenerator: with "
0067 << "beam along eta = " << fEta << " (Theta = " << fTheta / deg
0068 << ") phi = " << fPhi / deg << ") psi = " << fPsi / deg << " centred at (" << fMeanX
0069 << ", " << fMeanY << ", " << fMeanZ << ") "
0070 << "and spread (" << fSigmaX << ", " << fSigmaY << ") of type Gaussian = " << fType
0071 << " use file " << ffile;
0072 if (ffile)
0073 edm::LogInfo("VertexGenerator") << "With " << nBinx << " bins "
0074 << " along X and " << nBiny << " bins along Y";
0075 }
0076
0077 BeamProfileVtxGenerator::~BeamProfileVtxGenerator() {}
0078
0079
0080 HepMC::FourVector BeamProfileVtxGenerator::newVertex(CLHEP::HepRandomEngine* engine) const {
0081 double aX, aY;
0082 if (ffile) {
0083 double r1 = engine->flat();
0084 int ixy = 0, ix, iy;
0085 for (unsigned int i = 0; i < fdistn.size(); i++) {
0086 if (r1 > fdistn[i])
0087 ixy = i + 1;
0088 }
0089 if (ixy >= nBinx * nBiny) {
0090 ix = nBinx - 1;
0091 iy = nBiny - 1;
0092 } else {
0093 ix = ixy % nBinx;
0094 iy = (ixy - ix) / nBinx;
0095 }
0096 aX = 0.5 * (2 * ix - nBinx + 2 * engine->flat()) * fSigmaX + fMeanX;
0097 aY = 0.5 * (2 * iy - nBiny + 2 * engine->flat()) * fSigmaY + fMeanY;
0098 } else {
0099 if (fType) {
0100 aX = fSigmaX * CLHEP::RandGaussQ::shoot(engine) + fMeanX;
0101 aY = fSigmaY * CLHEP::RandGaussQ::shoot(engine) + fMeanY;
0102 } else {
0103 aX = CLHEP::RandFlat::shoot(engine, -0.5 * fSigmaX, 0.5 * fSigmaX) + fMeanX;
0104 aY = CLHEP::RandFlat::shoot(engine, -0.5 * fSigmaY, 0.5 * fSigmaY) + fMeanY;
0105 }
0106 }
0107
0108 double xp, yp, zp;
0109 if (2. * M_PI < fabs(fPsi)) {
0110 xp = -aX * cos(fTheta) * cos(fPhi) + aY * sin(fPhi) + fMeanZ * sin(fTheta) * cos(fPhi);
0111 yp = -aX * cos(fTheta) * sin(fPhi) - aY * cos(fPhi) + fMeanZ * sin(fTheta) * sin(fPhi);
0112 zp = aX * sin(fTheta) + fMeanZ * cos(fTheta);
0113 } else
0114 {
0115 const HepGeom::Vector3D<double> av(aX, aY, fMeanZ);
0116
0117
0118
0119
0120
0121
0122
0123
0124 const HepGeom::RotateZ3D R1(fPhi - M_PI);
0125 const HepGeom::Point3D<double> xUnit(0, 1, 0);
0126 const HepGeom::Point3D<double> zUnit(0, 0, 1);
0127 const HepGeom::Transform3D RXRZ(HepGeom::Rotate3D(-fTheta, R1 * xUnit) * R1);
0128 const HepGeom::Transform3D TRF(HepGeom::Rotate3D(fPsi, RXRZ * zUnit) * RXRZ);
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138
0139
0140
0141
0142
0143
0144
0145
0146 const HepGeom::Vector3D<double> pv(TRF * av);
0147
0148 xp = pv.x();
0149 yp = pv.y();
0150 zp = pv.z();
0151 }
0152
0153 LogDebug("VertexGenerator") << "BeamProfileVtxGenerator: Vertex created "
0154 << "at (" << xp << ", " << yp << ", " << zp << ", " << fTimeOffset << ")";
0155 return HepMC::FourVector(xp, yp, zp, fTimeOffset);
0156 ;
0157 }
0158
0159 void BeamProfileVtxGenerator::sigmaX(double s) {
0160 if (s >= 0) {
0161 fSigmaX = s;
0162 } else {
0163 edm::LogWarning("VertexGenerator") << "Warning BeamProfileVtxGenerator:"
0164 << " Illegal resolution in X " << s << "- set to default value 0 cm";
0165 fSigmaX = 0;
0166 }
0167 }
0168
0169 void BeamProfileVtxGenerator::sigmaY(double s) {
0170 if (s >= 0) {
0171 fSigmaY = s;
0172 } else {
0173 edm::LogWarning("VertexGenerator") << "Warning BeamProfileVtxGenerator:"
0174 << " Illegal resolution in Y " << s << "- set to default value 0 cm";
0175 fSigmaY = 0;
0176 }
0177 }
0178
0179 void BeamProfileVtxGenerator::eta(double s) {
0180 fEta = s;
0181 fTheta = 2.0 * atan(exp(-fEta));
0182 }
0183
0184 void BeamProfileVtxGenerator::setType(bool s) { fType = s; }