CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/GeneratorInterface/HiGenCommon/plugins/MixBoostEvtVtxGenerator.cc

Go to the documentation of this file.
00001 
00002 // $Id: MixBoostEvtVtxGenerator.cc,v 1.1 2012/06/08 22:19:37 yilmaz Exp $
00003 /*
00004 ________________________________________________________________________
00005 
00006  MixBoostEvtVtxGenerator
00007 
00008  Smear vertex according to the Beta function on the transverse plane
00009  and a Gaussian on the z axis. It allows the beam to have a crossing
00010  angle (slopes dxdz and dydz).
00011 
00012  Based on GaussEvtVtxGenerator
00013  implemented by Francisco Yumiceva (yumiceva@fnal.gov)
00014 
00015  FERMILAB
00016  2006
00017 ________________________________________________________________________
00018 */
00019 
00020 //lingshan: add beta for z-axis boost
00021 
00022 //#include "IOMC/EventVertexGenerators/interface/BetafuncEvtVtxGenerator.h"
00023 
00024 #include "FWCore/Framework/interface/Event.h"
00025 #include "FWCore/Utilities/interface/Exception.h"
00026 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00027 
00028 #include "FWCore/Framework/interface/EDProducer.h"
00029 #include "FWCore/Utilities/interface/InputTag.h"
00030 #include "FWCore/ServiceRegistry/interface/Service.h"
00031 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
00032 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
00033 #include "DataFormats/VertexReco/interface/Vertex.h"
00034 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00035 
00036 #include "CLHEP/Random/RandGaussQ.h"
00037 #include "CLHEP/Units/GlobalSystemOfUnits.h"
00038 #include "CLHEP/Units/GlobalPhysicalConstants.h"
00039 //#include "CLHEP/Vector/ThreeVector.h"
00040 #include "HepMC/SimpleVector.h"
00041 #include "TMatrixD.h"
00042 
00043 #include <iostream>
00044 
00045 using namespace edm;
00046 using namespace std;
00047 using namespace CLHEP;
00048 
00049 class RandGaussQ;
00050 class FourVector ;
00051 
00052 class MixBoostEvtVtxGenerator : public edm::EDProducer{
00053 public:
00054   MixBoostEvtVtxGenerator(const edm::ParameterSet & p);
00055   virtual ~MixBoostEvtVtxGenerator();
00056 
00058   //virtual CLHEP::Hep3Vector * newVertex();
00059   virtual HepMC::FourVector* newVertex() ;
00060   virtual void produce( edm::Event&, const edm::EventSetup& );
00061   virtual TMatrixD* GetInvLorentzBoost();
00062   virtual HepMC::FourVector* getVertex(edm::Event&);
00063   virtual HepMC::FourVector* getRecVertex(edm::Event&);
00064 
00066   void sigmaZ(double s=1.0);
00067 
00069   void X0(double m=0) { fX0=m; }
00071   void Y0(double m=0) { fY0=m; }
00073   void Z0(double m=0) { fZ0=m; }
00074 
00076   void Phi(double m=0) { phi_=m; }
00078   void Alpha(double m=0) { alpha_=m; }
00079   void Beta(double m=0) { beta_=m; }
00080 
00082   void betastar(double m=0) { fbetastar=m; }
00084   void emittance(double m=0) { femittance=m; }
00085 
00087   double BetaFunction(double z, double z0);
00088   CLHEP::HepRandomEngine& getEngine();
00089 
00090 private:
00092   MixBoostEvtVtxGenerator(const MixBoostEvtVtxGenerator &p);
00094   MixBoostEvtVtxGenerator&  operator = (const MixBoostEvtVtxGenerator & rhs );
00095 
00096   double alpha_, phi_;
00097   //TMatrixD boost_;
00098   double beta_;
00099   double fX0, fY0, fZ0;
00100   double fSigmaZ;
00101   //double fdxdz, fdydz;
00102   double fbetastar, femittance;
00103   double falpha;
00104 
00105   HepMC::FourVector*       fVertex ;
00106   TMatrixD *boost_;
00107   double fTimeOffset;
00108   
00109   CLHEP::HepRandomEngine*  fEngine;
00110   edm::InputTag            sourceLabel;
00111 
00112   CLHEP::RandGaussQ*  fRandom ;
00113 
00114   edm::InputTag            signalLabel;
00115   edm::InputTag            hiLabel;
00116   bool                     useRecVertex;
00117   std::vector<double>      vtxOffset;
00118 
00119 };
00120 
00121 
00122 MixBoostEvtVtxGenerator::MixBoostEvtVtxGenerator(const edm::ParameterSet & pset ):
00123   fVertex(0), boost_(0), fTimeOffset(0), fEngine(0),
00124   signalLabel(pset.getParameter<edm::InputTag>("signalLabel")),
00125   hiLabel(pset.getParameter<edm::InputTag>("heavyIonLabel")),
00126   useRecVertex(pset.exists("useRecVertex")?pset.getParameter<bool>("useRecVertex"):false)
00127 { 
00128 
00129   vtxOffset.resize(3);
00130   if(pset.exists("vtxOffset")) vtxOffset=pset.getParameter< std::vector<double> >("vtxOffset"); 
00131   edm::Service<edm::RandomNumberGenerator> rng;
00132 
00133   if ( ! rng.isAvailable()) {
00134     
00135     throw cms::Exception("Configuration")
00136       << "The BaseEvtVtxGenerator requires the RandomNumberGeneratorService\n"
00137       "which is not present in the configuration file.  You must add the service\n"
00138       "in the configuration file or remove the modules that require it.";
00139   }
00140 
00141   CLHEP::HepRandomEngine& engine = rng->getEngine();
00142   fEngine = &engine;
00143 
00144   produces<bool>("matchedVertex"); 
00145   
00146 }
00147 
00148 MixBoostEvtVtxGenerator::~MixBoostEvtVtxGenerator() 
00149 {
00150   delete fVertex ;
00151   if (boost_ != 0 ) delete boost_;
00152   delete fRandom; 
00153 }
00154 
00155 CLHEP::HepRandomEngine& MixBoostEvtVtxGenerator::getEngine(){
00156   return *fEngine;
00157 }
00158 
00159 
00160 //Hep3Vector* MixBoostEvtVtxGenerator::newVertex() {
00161 HepMC::FourVector* MixBoostEvtVtxGenerator::newVertex() {
00162 
00163         
00164         double X,Y,Z;
00165         
00166         double tmp_sigz = fRandom->fire(0., fSigmaZ);
00167         Z = tmp_sigz + fZ0;
00168 
00169         double tmp_sigx = BetaFunction(Z,fZ0); 
00170         // need sqrt(2) for beamspot width relative to single beam width
00171         tmp_sigx /= sqrt(2.0);
00172         X = fRandom->fire(0.,tmp_sigx) + fX0; // + Z*fdxdz ;
00173 
00174         double tmp_sigy = BetaFunction(Z,fZ0);
00175         // need sqrt(2) for beamspot width relative to single beam width
00176         tmp_sigy /= sqrt(2.0);
00177         Y = fRandom->fire(0.,tmp_sigy) + fY0; // + Z*fdydz;
00178 
00179         double tmp_sigt = fRandom->fire(0., fSigmaZ);
00180         double T = tmp_sigt + fTimeOffset; 
00181 
00182         if ( fVertex == 0 ) fVertex = new HepMC::FourVector();
00183         fVertex->set(X,Y,Z,T);
00184                 
00185         return fVertex;
00186 }
00187 
00188 double MixBoostEvtVtxGenerator::BetaFunction(double z, double z0)
00189 {
00190         return sqrt(femittance*(fbetastar+(((z-z0)*(z-z0))/fbetastar)));
00191 
00192 }
00193 
00194 
00195 void MixBoostEvtVtxGenerator::sigmaZ(double s) 
00196 { 
00197         if (s>=0 ) {
00198                 fSigmaZ=s; 
00199         }
00200         else {
00201                 throw cms::Exception("LogicError")
00202                         << "Error in MixBoostEvtVtxGenerator::sigmaZ: "
00203                         << "Illegal resolution in Z (negative)";
00204         }
00205 }
00206 
00207 TMatrixD* MixBoostEvtVtxGenerator::GetInvLorentzBoost() {
00208 
00209         //alpha_ = 0;
00210         //phi_ = 142.e-6;
00211 //      if (boost_ != 0 ) return boost_;
00212         
00213         //boost_.ResizeTo(4,4);
00214         //boost_ = new TMatrixD(4,4);
00215         TMatrixD tmpboost(4,4);
00216         TMatrixD tmpboostZ(4,4);
00217         TMatrixD tmpboostXYZ(4,4);
00218 
00219         //if ( (alpha_ == 0) && (phi_==0) ) { boost_->Zero(); return boost_; }
00220         
00221         // Lorentz boost to frame where the collision is head-on
00222         // phi is the half crossing angle in the plane ZS
00223         // alpha is the angle to the S axis from the X axis in the XY plane
00224         
00225         tmpboost(0,0) = 1./cos(phi_);
00226         tmpboost(0,1) = - cos(alpha_)*sin(phi_);
00227         tmpboost(0,2) = - tan(phi_)*sin(phi_);
00228         tmpboost(0,3) = - sin(alpha_)*sin(phi_);
00229         tmpboost(1,0) = - cos(alpha_)*tan(phi_);
00230         tmpboost(1,1) = 1.;
00231         tmpboost(1,2) = cos(alpha_)*tan(phi_);
00232         tmpboost(1,3) = 0.;
00233         tmpboost(2,0) = 0.;
00234         tmpboost(2,1) = - cos(alpha_)*sin(phi_);
00235         tmpboost(2,2) = cos(phi_);
00236         tmpboost(2,3) = - sin(alpha_)*sin(phi_);
00237         tmpboost(3,0) = - sin(alpha_)*tan(phi_);
00238         tmpboost(3,1) = 0.;
00239         tmpboost(3,2) = sin(alpha_)*tan(phi_);
00240         tmpboost(3,3) = 1.;
00241        //cout<<"beta "<<beta_;
00242        double gama=1.0/sqrt(1-beta_*beta_);
00243        tmpboostZ(0,0)=gama;
00244        tmpboostZ(0,1)=0.;
00245        tmpboostZ(0,2)=-1.0*beta_*gama;
00246        tmpboostZ(0,3)=0.;
00247        tmpboostZ(1,0)=0.;
00248        tmpboostZ(1,1) = 1.;
00249        tmpboostZ(1,2)=0.;
00250        tmpboostZ(1,3)=0.;
00251        tmpboostZ(2,0)=-1.0*beta_*gama;
00252        tmpboostZ(2,1) = 0.;
00253        tmpboostZ(2,2)=gama;
00254        tmpboostZ(2,3) = 0.;
00255        tmpboostZ(3,0)=0.;
00256        tmpboostZ(3,1)=0.;
00257        tmpboostZ(3,2)=0.;
00258        tmpboostZ(3,3) = 1.;
00259 
00260        tmpboostXYZ=tmpboost*tmpboostZ;
00261        tmpboost.Invert();
00262 
00263 
00264 
00265        boost_ = new TMatrixD(tmpboostXYZ);
00266        boost_->Print();
00267         
00268         return boost_;
00269 }
00270 
00271 HepMC::FourVector* MixBoostEvtVtxGenerator::getVertex( Event& evt){
00272   
00273   Handle<HepMCProduct> input;
00274   evt.getByLabel(hiLabel,input);
00275 
00276   const HepMC::GenEvent* inev = input->GetEvent();
00277   HepMC::GenVertex* genvtx = inev->signal_process_vertex();
00278   if(!genvtx){
00279     cout<<"No Signal Process Vertex!"<<endl;
00280     HepMC::GenEvent::particle_const_iterator pt=inev->particles_begin();
00281     HepMC::GenEvent::particle_const_iterator ptend=inev->particles_end();
00282     while(!genvtx || ( genvtx->particles_in_size() == 1 && pt != ptend ) ){
00283       if(!genvtx) cout<<"No Gen Vertex!"<<endl;
00284       if(pt == ptend) cout<<"End reached!"<<endl;
00285       genvtx = (*pt)->production_vertex();
00286       ++pt;
00287     }
00288   }
00289   double aX,aY,aZ,aT;
00290   
00291   aX = genvtx->position().x();
00292   aY = genvtx->position().y();
00293   aZ = genvtx->position().z();
00294   aT = genvtx->position().t();
00295    
00296   if(!fVertex) fVertex = new HepMC::FourVector();
00297   fVertex->set(aX,aY,aZ,aT);
00298   
00299   return fVertex;
00300   
00301 }
00302  
00303  
00304 HepMC::FourVector* MixBoostEvtVtxGenerator::getRecVertex( Event& evt){
00305  
00306   Handle<reco::VertexCollection> input;
00307   evt.getByLabel(hiLabel,input);
00308 
00309   double aX,aY,aZ;
00310  
00311   aX = input->begin()->position().x() + vtxOffset[0];
00312   aY = input->begin()->position().y() + vtxOffset[1];
00313   aZ = input->begin()->position().z() + vtxOffset[2];
00314  
00315   if(!fVertex) fVertex = new HepMC::FourVector();
00316   fVertex->set(10.0*aX,10.0*aY,10.0*aZ,0.0); // HepMC positions in mm (RECO in cm)
00317    
00318   return fVertex;
00319  
00320 }
00321 
00322 
00323 void MixBoostEvtVtxGenerator::produce( Event& evt, const EventSetup& )
00324 {
00325     
00326     
00327   Handle<HepMCProduct> HepMCEvt ;
00328   
00329   evt.getByLabel( signalLabel, HepMCEvt ) ;
00330     
00331   // generate new vertex & apply the shift 
00332   //
00333  
00334   HepMCEvt->applyVtxGen( useRecVertex ? getRecVertex(evt) : getVertex(evt) ) ;
00335  
00336   //   HepMCEvt->boostToLab( GetInvLorentzBoost(), "vertex" );
00337   //   HepMCEvt->boostToLab( GetInvLorentzBoost(), "momentum" );
00338   
00339   // OK, create a (pseudo)product and put in into edm::Event
00340   //
00341   auto_ptr<bool> NewProduct(new bool(true)) ;      
00342   evt.put( NewProduct ,"matchedVertex") ;
00343        
00344   return ;
00345   
00346 }
00347 
00348 #include "FWCore/Framework/interface/MakerMacros.h"
00349 DEFINE_FWK_MODULE(MixBoostEvtVtxGenerator);