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