ParticleGuns/src/BeamMomentumGunProducer.cc

0001 #include "IOMC/ParticleGuns/interface/BeamMomentumGunProducer.h"
0002
0003 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0004 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0005
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/ServiceRegistry/interface/Service.h"
0010 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0011
0012 #include "CLHEP/Random/RandFlat.h"
0013
0014 #include "TFile.h"
0015
0016 namespace CLHEP {
0017   class HepRandomEngine;
0018 }
0019
0020 namespace edm {
0021   BeamMomentumGunProducer::BeamMomentumGunProducer(const edm::ParameterSet& pset)
0022       : FlatBaseThetaGunProducer(pset),
0023         parPDGId_(nullptr),
0024         parX_(nullptr),
0025         parY_(nullptr),
0026         parZ_(nullptr),
0027         parPx_(nullptr),
0028         parPy_(nullptr),
0029         parPz_(nullptr),
0030         b_npar_(nullptr),
0031         b_eventId_(nullptr),
0032         b_parPDGId_(nullptr),
0033         b_parX_(nullptr),
0034         b_parY_(nullptr),
0035         b_parZ_(nullptr),
0036         b_parPx_(nullptr),
0037         b_parPy_(nullptr),
0038         b_parPz_(nullptr) {
0039     edm::ParameterSet pgun_params = pset.getParameter<edm::ParameterSet>("PGunParameters");
0040
0041     // doesn't seem necessary to check if pset is empty
0042     xoff_ = pgun_params.getParameter<double>("XOffset");
0043     yoff_ = pgun_params.getParameter<double>("YOffset");
0044     zpos_ = pgun_params.getParameter<double>("ZPosition");
0045     if (fVerbosity > 0)
0046       edm::LogVerbatim("BeamMomentumGun")
0047           << "Beam vertex offset (cm) " << xoff_ << ":" << yoff_ << " and z position " << zpos_;
0048
0049     edm::FileInPath fp = pgun_params.getParameter<edm::FileInPath>("FileName");
0050     const std::string& infileName = fp.fullPath();
0051
0052     fFile_ = new TFile(infileName.c_str());
0053     fFile_->GetObject("EventTree", fTree_);
0054     nentries_ = fTree_->GetEntriesFast();
0055     if (fVerbosity > 0)
0056       edm::LogVerbatim("BeamMomentumGun") << "Total Events: " << nentries_ << " in " << infileName;
0057
0058     // Set branch addresses and branch pointers
0059     int npart = fTree_->SetBranchAddress("npar", &npar_, &b_npar_);
0060     int event = fTree_->SetBranchAddress("eventId", &eventId_, &b_eventId_);
0061     int pdgid = fTree_->SetBranchAddress("parPDGId", &parPDGId_, &b_parPDGId_);
0062     int parxx = fTree_->SetBranchAddress("parX", &parX_, &b_parX_);
0063     int paryy = fTree_->SetBranchAddress("parY", &parY_, &b_parY_);
0064     int parzz = fTree_->SetBranchAddress("parZ", &parZ_, &b_parZ_);
0065     int parpx = fTree_->SetBranchAddress("parPx", &parPx_, &b_parPx_);
0066     int parpy = fTree_->SetBranchAddress("parPy", &parPy_, &b_parPy_);
0067     int parpz = fTree_->SetBranchAddress("parPz", &parPz_, &b_parPz_);
0068     if ((npart != 0) || (event != 0) || (pdgid != 0) || (parxx != 0) || (paryy != 0) || (parzz != 0) || (parpx != 0) ||
0069         (parpy != 0) || (parpz != 0))
0070       throw cms::Exception("GenException") << "Branch address wrong in i/p file\n";
0071
0072     produces<HepMCProduct>("unsmeared");
0073     produces<GenEventInfoProduct>();
0074
0075     if (fVerbosity > 0)
0076       edm::LogVerbatim("BeamMomentumGun") << "BeamMonetumGun is initialzed";
0077   }
0078
0079   void BeamMomentumGunProducer::produce(edm::Event& e, const edm::EventSetup& es) {
0080     if (fVerbosity > 0)
0081       edm::LogVerbatim("BeamMomentumGun") << "BeamMomentumGunProducer : Begin New Event Generation";
0082
0083     edm::Service<edm::RandomNumberGenerator> rng;
0084     CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
0085
0086     // event loop (well, another step in it...)
0087     // no need to clean up GenEvent memory - done in HepMCProduct
0088     // here re-create fEvt (memory)
0089     //
0090     fEvt = new HepMC::GenEvent();
0091
0092     // random entry generation for peaking event randomly from tree
0093     long int rjentry = static_cast<long int>(CLHEP::RandFlat::shoot(engine, 0, nentries_ - 1));
0094     fTree_->GetEntry(rjentry);
0095     if (fVerbosity > 0)
0096       edm::LogVerbatim("BeamMomentumGun") << "Entry " << rjentry << " : " << eventId_ << "  :  " << npar_;
0097
0098     // loop over particles
0099     int barcode = 1;
0100     for (unsigned int ip = 0; ip < parPDGId_->size(); ip++) {
0101       int partID = parPDGId_->at(ip);
0102       const HepPDT::ParticleData* pData = fPDGTable->particle(HepPDT::ParticleID(std::abs(partID)));
0103       double mass = pData->mass().value();
0104       if (fVerbosity > 0)
0105         edm::LogVerbatim("BeamMomentumGun") << "PDGId: " << partID << "   mass: " << mass;
0106       double xp = (xoff_ * cm2mm_ + (-1) * parY_->at(ip));  // 90 degree rotation applied
0107       double yp = (yoff_ * cm2mm_ + parX_->at(ip));         // 90 degree rotation applied
0108       double zp = zpos_ * cm2mm_;
0109       HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(xp, yp, zp));
0110       double pxGeV = MeV2GeV_ * (-1) * parPy_->at(ip);  // 90 degree rotation applied
0111       double pyGeV = MeV2GeV_ * parPx_->at(ip);         // 90 degree rotation applied
0112       double pzGeV = MeV2GeV_ * parPz_->at(ip);
0113       double theta = CLHEP::RandFlat::shoot(engine, fMinTheta, fMaxTheta);
0114       double phi = CLHEP::RandFlat::shoot(engine, fMinPhi, fMaxPhi);
0115       // rotation about Z axis
0116       double px1 = pxGeV * cos(phi) - pyGeV * sin(phi);
0117       double py1 = pxGeV * sin(phi) + pyGeV * cos(phi);
0118       double pz1 = pzGeV;
0119       // rotation about Y axis
0120       double px = px1 * cos(theta) + pz1 * sin(theta);
0121       double py = py1;
0122       double pz = -px1 * sin(theta) + pz1 * cos(theta);
0123       double energy = std::sqrt(px * px + py * py + pz * pz + mass * mass);
0124
0125       if (fVerbosity > 0) {
0126         edm::LogVerbatim("BeamMomentumGun") << "x:y:z [mm] " << xp << ":" << yp << ":" << zpos_;
0127         edm::LogVerbatim("BeamMomentumGun") << "px:py:pz [GeV] " << px << ":" << py << ":" << pz;
0128       }
0129
0130       HepMC::FourVector p(px, py, pz, energy);
0131       HepMC::GenParticle* part = new HepMC::GenParticle(p, partID, 1);
0132       part->suggest_barcode(barcode);
0133       barcode++;
0134       Vtx->add_particle_out(part);
0135
0136       if (fAddAntiParticle) {
0137         HepMC::FourVector ap(-px, -py, -pz, energy);
0138         int apartID = (partID == 22 || partID == 23) ? partID : -partID;
0139         HepMC::GenParticle* apart = new HepMC::GenParticle(ap, apartID, 1);
0140         apart->suggest_barcode(barcode);
0141         if (fVerbosity > 0)
0142           edm::LogVerbatim("BeamMomentumGun")
0143               << "Add anti-particle " << apartID << ":" << -px << ":" << -py << ":" << -pz;
0144         barcode++;
0145         Vtx->add_particle_out(apart);
0146       }
0147
0148       fEvt->add_vertex(Vtx);
0149     }
0150
0151     fEvt->set_event_number(e.id().event());
0152     fEvt->set_signal_process_id(20);
0153
0154     if (fVerbosity > 0)
0155       fEvt->print();
0156
0157     std::unique_ptr<HepMCProduct> BProduct(new HepMCProduct());
0158     BProduct->addHepMCData(fEvt);
0159     e.put(std::move(BProduct), "unsmeared");
0160
0161     std::unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(fEvt));
0162     e.put(std::move(genEventInfo));
0163
0164     if (fVerbosity > 0)
0165       edm::LogVerbatim("BeamMomentumGun") << "BeamMomentumGunProducer : Event Generation Done";
0166   }
0167 }  // namespace edm