CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
BetafuncEvtVtxGenerator.cc
Go to the documentation of this file.
1 
2 // $Id: BetafuncEvtVtxGenerator.cc,v 1.11 2010/12/06 10:57:25 yumiceva Exp $
3 /*
4 ________________________________________________________________________
5 
6  BetafuncEvtVtxGenerator
7 
8  Smear vertex according to the Beta function on the transverse plane
9  and a Gaussian on the z axis. It allows the beam to have a crossing
10  angle (slopes dxdz and dydz).
11 
12  Based on GaussEvtVtxGenerator
13  implemented by Francisco Yumiceva (yumiceva@fnal.gov)
14 
15  FERMILAB
16  2006
17 ________________________________________________________________________
18 */
19 
20 
21 
25 
26 #include "CLHEP/Random/RandGaussQ.h"
27 #include "CLHEP/Units/GlobalSystemOfUnits.h"
28 #include "CLHEP/Units/GlobalPhysicalConstants.h"
29 //#include "CLHEP/Vector/ThreeVector.h"
30 #include "HepMC/SimpleVector.h"
31 
32 #include <iostream>
33 
34 
35 
38 {
39 
40  fRandom = new CLHEP::RandGaussQ(getEngine());
41 
42  fX0 = p.getParameter<double>("X0")*cm;
43  fY0 = p.getParameter<double>("Y0")*cm;
44  fZ0 = p.getParameter<double>("Z0")*cm;
45  fSigmaZ = p.getParameter<double>("SigmaZ")*cm;
46  alpha_ = p.getParameter<double>("Alpha")*radian;
47  phi_ = p.getParameter<double>("Phi")*radian;
48  fbetastar = p.getParameter<double>("BetaStar")*cm;
49  femittance = p.getParameter<double>("Emittance")*cm; // this is not the normalized emittance
50  fTimeOffset = p.getParameter<double>("TimeOffset")*ns*c_light; // HepMC time units are mm
51 
52  if (fSigmaZ <= 0) {
53  throw cms::Exception("Configuration")
54  << "Error in BetafuncEvtVtxGenerator: "
55  << "Illegal resolution in Z (SigmaZ is negative)";
56  }
57 
58 
59 }
60 
62 {
63  delete fRandom;
64 }
65 
66 
67 //Hep3Vector* BetafuncEvtVtxGenerator::newVertex() {
68 HepMC::FourVector* BetafuncEvtVtxGenerator::newVertex() {
69 
70 
71  double X,Y,Z;
72 
73  double tmp_sigz = fRandom->fire(0., fSigmaZ);
74  Z = tmp_sigz + fZ0;
75 
76  double tmp_sigx = BetaFunction(Z,fZ0);
77  // need sqrt(2) for beamspot width relative to single beam width
78  tmp_sigx /= sqrt(2.0);
79  X = fRandom->fire(0.,tmp_sigx) + fX0; // + Z*fdxdz ;
80 
81  double tmp_sigy = BetaFunction(Z,fZ0);
82  // need sqrt(2) for beamspot width relative to single beam width
83  tmp_sigy /= sqrt(2.0);
84  Y = fRandom->fire(0.,tmp_sigy) + fY0; // + Z*fdydz;
85 
86  double tmp_sigt = fRandom->fire(0., fSigmaZ);
87  double T = tmp_sigt + fTimeOffset;
88 
89  if ( fVertex == 0 ) fVertex = new HepMC::FourVector();
90  fVertex->set(X,Y,Z,T);
91 
92  return fVertex;
93 }
94 
95 double BetafuncEvtVtxGenerator::BetaFunction(double z, double z0)
96 {
97  return sqrt(femittance*(fbetastar+(((z-z0)*(z-z0))/fbetastar)));
98 
99 }
100 
101 
103 {
104  if (s>=0 ) {
105  fSigmaZ=s;
106  }
107  else {
108  throw cms::Exception("LogicError")
109  << "Error in BetafuncEvtVtxGenerator::sigmaZ: "
110  << "Illegal resolution in Z (negative)";
111  }
112 }
113 
115 
116  //alpha_ = 0;
117  //phi_ = 142.e-6;
118 
119  if (boost_ != 0 ) return boost_;
120 
121  //boost_.ResizeTo(4,4);
122  //boost_ = new TMatrixD(4,4);
123  TMatrixD tmpboost(4,4);
124 
125  //if ( (alpha_ == 0) && (phi_==0) ) { boost_->Zero(); return boost_; }
126 
127  // Lorentz boost to frame where the collision is head-on
128  // phi is the half crossing angle in the plane ZS
129  // alpha is the angle to the S axis from the X axis in the XY plane
130 
131  tmpboost(0,0) = 1./cos(phi_);
132  tmpboost(0,1) = - cos(alpha_)*sin(phi_);
133  tmpboost(0,2) = - tan(phi_)*sin(phi_);
134  tmpboost(0,3) = - sin(alpha_)*sin(phi_);
135  tmpboost(1,0) = - cos(alpha_)*tan(phi_);
136  tmpboost(1,1) = 1.;
137  tmpboost(1,2) = cos(alpha_)*tan(phi_);
138  tmpboost(1,3) = 0.;
139  tmpboost(2,0) = 0.;
140  tmpboost(2,1) = - cos(alpha_)*sin(phi_);
141  tmpboost(2,2) = cos(phi_);
142  tmpboost(2,3) = - sin(alpha_)*sin(phi_);
143  tmpboost(3,0) = - sin(alpha_)*tan(phi_);
144  tmpboost(3,1) = 0.;
145  tmpboost(3,2) = sin(alpha_)*tan(phi_);
146  tmpboost(3,3) = 1.;
147 
148  tmpboost.Invert();
149  boost_ = new TMatrixD(tmpboost);
150  //boost_->Print();
151 
152  return boost_;
153 }
const double Z[kNumberCalorimeter]
T getParameter(std::string const &) const
void sigmaZ(double s=1.0)
set resolution in Z in cm
CLHEP::HepRandomEngine & getEngine()
virtual HepMC::FourVector * newVertex()
return a new event vertex
HepMC::FourVector * fVertex
BetafuncEvtVtxGenerator(const edm::ParameterSet &p)
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
#define X(str)
Definition: MuonsGrabber.cc:49
double BetaFunction(double z, double z0)
beta function
virtual TMatrixD * GetInvLorentzBoost()
Definition: DDAxes.h:10
T sqrt(T t)
Definition: SSEVec.h:28
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
string s
Definition: asciidump.py:422