EventVertexGenerators/src/BeamProfileVtxGenerator.cc

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 //Hep3Vector * BeamProfileVtxGenerator::newVertex() {
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  // for endcap testbeam, use 3rd angle psi
0114   {
0115     const HepGeom::Vector3D<double> av(aX, aY, fMeanZ);
0116
0117     /*
0118      static const double kRadToDeg ( 180./M_PI ) ;
0119      std::cout<<"theta = "<<std::setw(7) << std::setiosflags( std::ios::fixed ) << std::setprecision(5)
0120           <<fTheta*kRadToDeg<<", phi="<<std::setw(7) << std::setiosflags( std::ios::fixed ) << std::setprecision(5)
0121           << fPhi*kRadToDeg <<", PSI="<<std::setw(7) << std::setiosflags( std::ios::fixed ) << std::setprecision(5)
0122           <<fPsi*kRadToDeg<<std::endl ;
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      std::cout<<"\n\n$$$$$$$$$$$Transform="
0131           <<" thetaZ = "<<std::setw(7) << std::setiosflags( std::ios::fixed ) << std::setprecision(5)
0132           <<TRF.getRotation().thetaZ()*kRadToDeg
0133           <<", phiZ = "<<std::setw(7) << std::setiosflags( std::ios::fixed ) << std::setprecision(5)
0134           <<TRF.getRotation().phiZ()*kRadToDeg
0135           <<", thetaY = "<<std::setw(7) << std::setiosflags( std::ios::fixed ) << std::setprecision(5)
0136           <<TRF.getRotation().thetaY()*kRadToDeg
0137           <<", phiY = "<<std::setw(7) << std::setiosflags( std::ios::fixed ) << std::setprecision(5)
0138           <<TRF.getRotation().phiY()*kRadToDeg
0139           <<", thetaX = "<<std::setw(7) << std::setiosflags( std::ios::fixed ) << std::setprecision(5)
0140           <<TRF.getRotation().thetaX()*kRadToDeg
0141           <<", phiX = "<<std::setw(7) << std::setiosflags( std::ios::fixed ) << std::setprecision(5)
0142           <<TRF.getRotation().phiX()*kRadToDeg
0143           <<std::setw(7) << std::setiosflags( std::ios::fixed ) << std::setprecision(5)
0144           <<TRF.getTranslation()<<std::endl ;
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; }