CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MixBoostEvtVtxGenerator.cc
Go to the documentation of this file.
1 
5 
12 
13 #include "CLHEP/Units/GlobalSystemOfUnits.h"
14 #include "CLHEP/Units/GlobalPhysicalConstants.h"
15 #include "HepMC/SimpleVector.h"
16 #include "TMatrixD.h"
17 
18 #include <iostream>
19 
20 using namespace edm;
21 using namespace std;
22 using namespace CLHEP;
23 
24 class RandGaussQ;
25 class FourVector ;
26 
28 public:
30  virtual ~MixBoostEvtVtxGenerator();
31 
32  virtual HepMC::FourVector* newVertex() ;
33  virtual void produce( edm::Event&, const edm::EventSetup& );
34  virtual TMatrixD* GetInvLorentzBoost();
35  virtual HepMC::FourVector* getVertex(edm::Event&);
36  virtual HepMC::FourVector* getRecVertex(edm::Event&);
37 
39  void sigmaZ(double s=1.0);
40 
42  void X0(double m=0) { fX0=m; }
44  void Y0(double m=0) { fY0=m; }
46  void Z0(double m=0) { fZ0=m; }
47 
49  void Phi(double m=0) { phi_=m; }
51  void Alpha(double m=0) { alpha_=m; }
52  void Beta(double m=0) { beta_=m; }
53 
55  void betastar(double m=0) { fbetastar=m; }
57  void emittance(double m=0) { femittance=m; }
58 
60  double BetaFunction(double z, double z0);
61  // CLHEP::HepRandomEngine& getEngine();
62 
63 private:
68 
69  double alpha_, phi_;
70  //TMatrixD boost_;
71  double beta_;
72  double fX0, fY0, fZ0;
73  double fSigmaZ;
74  //double fdxdz, fdydz;
75  double fbetastar, femittance;
76  double falpha;
77 
78  HepMC::FourVector* fVertex ;
79  TMatrixD *boost_;
80  double fTimeOffset;
81 
82  // CLHEP::HepRandomEngine* fEngine;
84 
85  // CLHEP::RandGaussQ* fRandom ;
86 
90  std::vector<double> vtxOffset;
91  bool verbosity_;
92 
93 };
94 
95 
97  fVertex(0), boost_(0), fTimeOffset(0),
98  signalLabel(pset.getParameter<edm::InputTag>("signalLabel")),
99  hiLabel(pset.getParameter<edm::InputTag>("heavyIonLabel")),
100  useRecVertex(pset.exists("useRecVertex")?pset.getParameter<bool>("useRecVertex"):false),
101  verbosity_(pset.getUntrackedParameter<bool>("verbosity",false))
102 {
103 
104  vtxOffset.resize(3);
105  if(pset.exists("vtxOffset")) vtxOffset=pset.getParameter< std::vector<double> >("vtxOffset");
106  beta_ = pset.getParameter<double>("Beta");
107 
108  alpha_ = 0;
109  phi_ = 0;
110  if(pset.exists("Alpha")){
111  alpha_ = pset.getParameter<double>("Alpha")*radian;
112  phi_ = pset.getParameter<double>("Phi")*radian;
113  }
114 
115  produces<bool>("matchedVertex");
116 
117 }
118 
120 {
121  delete fVertex ;
122  if (boost_ != 0 ) delete boost_;
123 }
124 
125 
127  return 0;
128 }
129 
130 double MixBoostEvtVtxGenerator::BetaFunction(double z, double z0)
131 {
132  return sqrt(femittance*(fbetastar+(((z-z0)*(z-z0))/fbetastar)));
133 
134 }
135 
136 
138 {
139  if (s>=0 ) {
140  fSigmaZ=s;
141  }
142  else {
143  throw cms::Exception("LogicError")
144  << "Error in MixBoostEvtVtxGenerator::sigmaZ: "
145  << "Illegal resolution in Z (negative)";
146  }
147 }
148 
150 
151  TMatrixD tmpboost(4,4);
152  TMatrixD tmpboostZ(4,4);
153  TMatrixD tmpboostXYZ(4,4);
154 
155  // Lorentz boost to frame where the collision is head-on
156  // phi is the half crossing angle in the plane ZS
157  // alpha is the angle to the S axis from the X axis in the XY plane
158 
159  tmpboost(0,0) = 1./cos(phi_);
160  tmpboost(0,1) = - cos(alpha_)*sin(phi_);
161  tmpboost(0,2) = - tan(phi_)*sin(phi_);
162  tmpboost(0,3) = - sin(alpha_)*sin(phi_);
163  tmpboost(1,0) = - cos(alpha_)*tan(phi_);
164  tmpboost(1,1) = 1.;
165  tmpboost(1,2) = cos(alpha_)*tan(phi_);
166  tmpboost(1,3) = 0.;
167  tmpboost(2,0) = 0.;
168  tmpboost(2,1) = - cos(alpha_)*sin(phi_);
169  tmpboost(2,2) = cos(phi_);
170  tmpboost(2,3) = - sin(alpha_)*sin(phi_);
171  tmpboost(3,0) = - sin(alpha_)*tan(phi_);
172  tmpboost(3,1) = 0.;
173  tmpboost(3,2) = sin(alpha_)*tan(phi_);
174  tmpboost(3,3) = 1.;
175  double gama=1.0/sqrt(1-beta_*beta_);
176  tmpboostZ(0,0)=gama;
177  tmpboostZ(0,1)=0.;
178  tmpboostZ(0,2)=-1.0*beta_*gama;
179  tmpboostZ(0,3)=0.;
180  tmpboostZ(1,0)=0.;
181  tmpboostZ(1,1) = 1.;
182  tmpboostZ(1,2)=0.;
183  tmpboostZ(1,3)=0.;
184  tmpboostZ(2,0)=-1.0*beta_*gama;
185  tmpboostZ(2,1) = 0.;
186  tmpboostZ(2,2)=gama;
187  tmpboostZ(2,3) = 0.;
188  tmpboostZ(3,0)=0.;
189  tmpboostZ(3,1)=0.;
190  tmpboostZ(3,2)=0.;
191  tmpboostZ(3,3) = 1.;
192 
193  tmpboostXYZ=tmpboost*tmpboostZ;
194  tmpboostXYZ.Invert();
195 
196  boost_ = new TMatrixD(tmpboostXYZ);
197  if ( verbosity_ )boost_->Print();
198 
199  return boost_;
200 }
201 
202 HepMC::FourVector* MixBoostEvtVtxGenerator::getVertex( Event& evt){
203 
205  evt.getByLabel(hiLabel,input);
206 
207  const HepMC::GenEvent* inev = input->GetEvent();
208  HepMC::GenVertex* genvtx = inev->signal_process_vertex();
209  if(!genvtx){
210  cout<<"No Signal Process Vertex!"<<endl;
211  HepMC::GenEvent::particle_const_iterator pt=inev->particles_begin();
212  HepMC::GenEvent::particle_const_iterator ptend=inev->particles_end();
213  while(!genvtx || ( genvtx->particles_in_size() == 1 && pt != ptend ) ){
214  if(!genvtx) cout<<"No Gen Vertex!"<<endl;
215  if(pt == ptend) cout<<"End reached!"<<endl;
216  genvtx = (*pt)->production_vertex();
217  ++pt;
218  }
219  }
220 
221  double aX,aY,aZ,aT;
222 
223  aX = genvtx->position().x();
224  aY = genvtx->position().y();
225  aZ = genvtx->position().z();
226  aT = genvtx->position().t();
227 
228  if(!fVertex) fVertex = new HepMC::FourVector();
229  fVertex->set(aX,aY,aZ,aT);
230 
231  return fVertex;
232 
233 }
234 
235 
237 
239  evt.getByLabel(hiLabel,input);
240 
241  double aX,aY,aZ;
242 
243  aX = input->begin()->position().x() + vtxOffset[0];
244  aY = input->begin()->position().y() + vtxOffset[1];
245  aZ = input->begin()->position().z() + vtxOffset[2];
246 
247  if(!fVertex) fVertex = new HepMC::FourVector();
248  fVertex->set(10.0*aX,10.0*aY,10.0*aZ,0.0); // HepMC positions in mm (RECO in cm)
249 
250  return fVertex;
251 
252 }
253 
254 
256 {
257 
258 
259  Handle<HepMCProduct> HepMCEvt ;
260 
261  evt.getByLabel( signalLabel, HepMCEvt ) ;
262 
263  // generate new vertex & apply the shift
264  //
265 
266  HepMCEvt->boostToLab( GetInvLorentzBoost(), "vertex" );
267  HepMCEvt->boostToLab( GetInvLorentzBoost(), "momentum" );
268 
269  HepMCEvt->applyVtxGen( useRecVertex ? getRecVertex(evt) : getVertex(evt) ) ;
270 
271  // OK, create a (pseudo)product and put in into edm::Event
272  //
273  auto_ptr<bool> NewProduct(new bool(true)) ;
274  evt.put( NewProduct ,"matchedVertex") ;
275 
276  return ;
277 
278 }
279 
nocap nocap const skelname & operator=(const skelname &)
T getParameter(std::string const &) const
void Z0(double m=0)
set mean in Z in cm
std::vector< double > vtxOffset
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
void emittance(double m=0)
emittance (no the normalized)
void X0(double m=0)
set mean in X in cm
bool exists(std::string const &parameterName) const
checks if a parameter exists
double double double z
void sigmaZ(double s=1.0)
set resolution in Z in cm
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:85
T sqrt(T t)
Definition: SSEVec.h:46
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
virtual HepMC::FourVector * newVertex()
virtual TMatrixD * GetInvLorentzBoost()
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
MixBoostEvtVtxGenerator(const edm::ParameterSet &p)
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:356
void Alpha(double m=0)
angle between crossing plane and horizontal plane
void Phi(double m=0)
set half crossing angle
void betastar(double m=0)
set beta_star
void Y0(double m=0)
set mean in Y in cm
virtual void produce(edm::Event &, const edm::EventSetup &)
tuple cout
Definition: gather_cfg.py:121
virtual HepMC::FourVector * getRecVertex(edm::Event &)
virtual HepMC::FourVector * getVertex(edm::Event &)
double BetaFunction(double z, double z0)
beta function