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