BetafuncEvtVtxGenerator.cc

Go to the documentation of this file.

/*
________________________________________________________________________

 BetafuncEvtVtxGenerator

 Smear vertex according to the Beta function on the transverse plane
 and a Gaussian on the z axis. It allows the beam to have a crossing
 angle (slopes dxdz and dydz).

 Based on GaussEvtVtxGenerator
 implemented by Francisco Yumiceva (yumiceva@fnal.gov)

 FERMILAB
 2006
________________________________________________________________________
*/



#include "IOMC/EventVertexGenerators/interface/BetafuncEvtVtxGenerator.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"

#include "CLHEP/Random/RandGaussQ.h"
#include "CLHEP/Units/GlobalSystemOfUnits.h"
#include "CLHEP/Units/GlobalPhysicalConstants.h"
//#include "CLHEP/Vector/ThreeVector.h"
#include "HepMC/SimpleVector.h"

#include "CondFormats/DataRecord/interface/SimBeamSpotObjectsRcd.h"
#include "CondFormats/BeamSpotObjects/interface/SimBeamSpotObjects.h"

#include <iostream>



BetafuncEvtVtxGenerator::BetafuncEvtVtxGenerator(const edm::ParameterSet & p )
: BaseEvtVtxGenerator(p),
  boost_(4,4)
{ 
  readDB_=p.getParameter<bool>("readDB");
  if (!readDB_){
    fX0 =        p.getParameter<double>("X0")*cm;
    fY0 =        p.getParameter<double>("Y0")*cm;
    fZ0 =        p.getParameter<double>("Z0")*cm;
    fSigmaZ =    p.getParameter<double>("SigmaZ")*cm;
    fbetastar =  p.getParameter<double>("BetaStar")*cm;
    femittance = p.getParameter<double>("Emittance")*cm; // this is not the normalized emittance
    fTimeOffset = p.getParameter<double>("TimeOffset")*ns*c_light; // HepMC time units are mm

    setBoost(p.getParameter<double>("Alpha")*radian,
             p.getParameter<double>("Phi")*radian);
    if (fSigmaZ <= 0) {
      throw cms::Exception("Configuration")
    << "Error in BetafuncEvtVtxGenerator: "
    << "Illegal resolution in Z (SigmaZ is negative)";
    }
  }
}

BetafuncEvtVtxGenerator::~BetafuncEvtVtxGenerator() 
{
}

void BetafuncEvtVtxGenerator::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const& iEventSetup){
  update(iEventSetup);
}
void BetafuncEvtVtxGenerator::beginRun(const edm::Run & , const edm::EventSetup& iEventSetup){
  update(iEventSetup);
}

void BetafuncEvtVtxGenerator::update(const edm::EventSetup& iEventSetup){
  if (readDB_ &&  parameterWatcher_.check(iEventSetup)){
    edm::ESHandle< SimBeamSpotObjects > beamhandle;
    iEventSetup.get<SimBeamSpotObjectsRcd>().get(beamhandle);

    fX0=beamhandle->fX0;
    fY0=beamhandle->fY0;
    fZ0=beamhandle->fZ0;
    //    falpha=beamhandle->fAlpha;
    fSigmaZ=beamhandle->fSigmaZ;
    fTimeOffset=beamhandle->fTimeOffset;
    fbetastar=beamhandle->fbetastar;
    femittance=beamhandle->femittance;
    setBoost(beamhandle->fAlpha,beamhandle->fPhi);

  }
}

HepMC::FourVector BetafuncEvtVtxGenerator::newVertex(CLHEP::HepRandomEngine* engine) const {

    
    double X,Y,Z;
    
    double tmp_sigz = CLHEP::RandGaussQ::shoot(engine, 0., fSigmaZ);
    Z = tmp_sigz + fZ0;

    double tmp_sigx = BetaFunction(Z,fZ0); 
    // need sqrt(2) for beamspot width relative to single beam width
    tmp_sigx /= sqrt(2.0);
    X = CLHEP::RandGaussQ::shoot(engine, 0., tmp_sigx) + fX0; // + Z*fdxdz ;

    double tmp_sigy = BetaFunction(Z,fZ0);
    // need sqrt(2) for beamspot width relative to single beam width
    tmp_sigy /= sqrt(2.0);
    Y = CLHEP::RandGaussQ::shoot(engine, 0., tmp_sigy) + fY0; // + Z*fdydz;

    double tmp_sigt = CLHEP::RandGaussQ::shoot(engine, 0., fSigmaZ);
    double T = tmp_sigt + fTimeOffset; 

    return HepMC::FourVector(X,Y,Z,T);
}

double BetafuncEvtVtxGenerator::BetaFunction(double z, double z0) const
{
    return sqrt(femittance*(fbetastar+(((z-z0)*(z-z0))/fbetastar)));

}

void BetafuncEvtVtxGenerator::setBoost(double alpha, double phi) 
{
  //boost_.ResizeTo(4,4);
  //boost_ = new TMatrixD(4,4);
  TMatrixD tmpboost(4,4);
  
  //if ( (alpha_ == 0) && (phi_==0) ) { boost_->Zero(); return boost_; }
  
  // Lorentz boost to frame where the collision is head-on
  // phi is the half crossing angle in the plane ZS
  // alpha is the angle to the S axis from the X axis in the XY plane
  
  tmpboost(0,0) = 1./cos(phi);
  tmpboost(0,1) = - cos(alpha)*sin(phi);
  tmpboost(0,2) = - tan(phi)*sin(phi);
  tmpboost(0,3) = - sin(alpha)*sin(phi);
  tmpboost(1,0) = - cos(alpha)*tan(phi);
  tmpboost(1,1) = 1.;
  tmpboost(1,2) = cos(alpha)*tan(phi);
  tmpboost(1,3) = 0.;
  tmpboost(2,0) = 0.;
  tmpboost(2,1) = - cos(alpha)*sin(phi);
  tmpboost(2,2) = cos(phi);
  tmpboost(2,3) = - sin(alpha)*sin(phi);
  tmpboost(3,0) = - sin(alpha)*tan(phi);
  tmpboost(3,1) = 0.;
  tmpboost(3,2) = sin(alpha)*tan(phi);
  tmpboost(3,3) = 1.;
  
  tmpboost.Invert();
  boost_ = tmpboost;
  //boost_->Print();
}

void BetafuncEvtVtxGenerator::sigmaZ(double s) 
{ 
    if (s>=0 ) {
        fSigmaZ=s; 
    }
    else {
        throw cms::Exception("LogicError")
            << "Error in BetafuncEvtVtxGenerator::sigmaZ: "
            << "Illegal resolution in Z (negative)";
    }
}

TMatrixD const* BetafuncEvtVtxGenerator::GetInvLorentzBoost() const {

  return &boost_;
}