Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-12-10 02:31:53

0001 /*
0002 */
0003 
0004 #include "IOMC/EventVertexGenerators/interface/BaseEvtVtxGenerator.h"
0005 
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008 
0009 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0010 #include "SimDataFormats/GeneratorProducts/interface/HepMC3Product.h"
0011 
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0014 
0015 #include "FWCore/Utilities/interface/Exception.h"
0016 
0017 #include "DataFormats/Provenance/interface/Provenance.h"
0018 #include "FWCore/Utilities/interface/EDMException.h"
0019 
0020 using namespace edm;
0021 using namespace CLHEP;
0022 
0023 BaseEvtVtxGenerator::BaseEvtVtxGenerator(const ParameterSet& pset) {
0024   Service<RandomNumberGenerator> rng;
0025   if (!rng.isAvailable()) {
0026     throw cms::Exception("Configuration") << "The BaseEvtVtxGenerator requires the RandomNumberGeneratorService\n"
0027                                              "which is not present in the configuration file. \n"
0028                                              "You must add the service\n"
0029                                              "in the configuration file or remove the modules that require it.";
0030   }
0031 
0032   sourceToken3 = consumes<edm::HepMC3Product>(pset.getParameter<edm::InputTag>("src"));
0033   sourceToken = consumes<edm::HepMCProduct>(pset.getParameter<edm::InputTag>("src"));
0034   produces<edm::HepMC3Product>();
0035   produces<edm::HepMCProduct>();
0036 }
0037 
0038 BaseEvtVtxGenerator::~BaseEvtVtxGenerator() {}
0039 
0040 void BaseEvtVtxGenerator::produce(Event& evt, const EventSetup&) {
0041   edm::Service<edm::RandomNumberGenerator> rng;
0042   CLHEP::HepRandomEngine* engine = &rng->getEngine(evt.streamID());
0043 
0044   Handle<HepMCProduct> HepUnsmearedMCEvt;
0045 
0046   bool found = evt.getByToken(sourceToken, HepUnsmearedMCEvt);
0047 
0048   if (found) {  // HepMC event exists
0049 
0050     // Make a copy
0051     HepMC::GenEvent* genevt = new HepMC::GenEvent(*HepUnsmearedMCEvt->GetEvent());
0052 
0053     std::unique_ptr<edm::HepMCProduct> HepMCEvt(new edm::HepMCProduct(genevt));
0054     // generate new vertex & apply the shift
0055     //
0056     ROOT::Math::XYZTVector VertexShift = vertexShift(engine);
0057     HepMCEvt->applyVtxGen(HepMC::FourVector(VertexShift.x(), VertexShift.y(), VertexShift.z(), VertexShift.t()));
0058 
0059     HepMCEvt->boostToLab(GetInvLorentzBoost(), "vertex");
0060     HepMCEvt->boostToLab(GetInvLorentzBoost(), "momentum");
0061 
0062     evt.put(std::move(HepMCEvt));
0063 
0064   } else {  // no HepMC event, try to get HepMC3 event
0065 
0066     Handle<HepMC3Product> HepUnsmearedMCEvt3;
0067     found = evt.getByToken(sourceToken3, HepUnsmearedMCEvt3);
0068 
0069     if (!found)
0070       throw cms::Exception("ProductAbsent") << "No HepMCProduct, tried to get HepMC3Product, but it is also absent.";
0071 
0072     HepMC3::GenEvent* genevt3 = new HepMC3::GenEvent();
0073     genevt3->read_data(*HepUnsmearedMCEvt3->GetEvent());
0074     HepMC3Product* productcopy3 = new HepMC3Product(genevt3);
0075     ROOT::Math::XYZTVector VertexShift = vertexShift(engine);
0076     productcopy3->applyVtxGen(HepMC3::FourVector(VertexShift.x(), VertexShift.y(), VertexShift.z(), VertexShift.t()));
0077 
0078     if (GetInvLorentzBoost() != nullptr) {
0079       TMatrixD tmplorentz(*GetInvLorentzBoost());
0080       TMatrixD p4(4, 1);
0081       p4(0, 0) = 1.;
0082       p4(1, 0) = 1.;
0083       p4(2, 0) = 1.;
0084       p4(3, 0) = 1.;  // Check if the boost matrix is not trivial
0085       TMatrixD p4lab(4, 1);
0086       p4lab = tmplorentz * p4;
0087       if (p4lab(0, 0) - p4(0, 0) != 0. || p4lab(1, 0) - p4(1, 0) != 0. || p4lab(2, 0) - p4(2, 0) != 0. ||
0088           p4lab(3, 0) - p4(3, 0) != 0.) {  // not trivial:
0089         productcopy3->boostToLab(GetInvLorentzBoost(), "vertex");
0090         productcopy3->boostToLab(GetInvLorentzBoost(), "momentum");
0091       }
0092     }
0093 
0094     std::unique_ptr<edm::HepMC3Product> HepMC3Evt(productcopy3);
0095     evt.put(std::move(HepMC3Evt));
0096   }
0097 
0098   return;
0099 }