Back to home page

Project CMSSW displayed by LXR

 
 

    


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