Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Run3-TB51 Correct the beam momentum generator for HGCal TB studies #29114

Merged
merged 3 commits into from Mar 13, 2020
Merged

Run3-TB51 Correct the beam momentum generator for HGCal TB studies #29114

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
b6eb179
Correct the beam momentum generator for HGCal TB studies
Mar 5, 2020
5fb8530
Code check
Mar 5, 2020
2f5fb2c
Take Kevin's comments
Mar 9, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
2 changes: 1 addition & 1 deletion IOMC/ParticleGuns/interface/BeamMomentumGunProducer.h
View file
Open in desktop
Expand Up @@ -37,7 +37,7 @@ namespace edm {
TBranch *b_parX_, *b_parY_, *b_parZ_;
TBranch *b_parPx_, *b_parPy_, *b_parPz_;

static constexpr double mm2cm_ = 0.1;
static constexpr double cm2mm_ = 10.0;
static constexpr double MeV2GeV_ = 0.001;
};
} // namespace edm
Expand Down
40 changes: 30 additions & 10 deletions IOMC/ParticleGuns/src/BeamMomentumGunProducer.cc
View file
Open in desktop
@@ -1,4 +1,3 @@
#include <ostream>
#include "IOMC/ParticleGuns/interface/BeamMomentumGunProducer.h"

#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
Expand All @@ -13,6 +12,8 @@
#include "CLHEP/Random/RandFlat.h"

#include "TFile.h"
#include <cmath>
#include <ostream>
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this include should not be needed, please remove it


namespace CLHEP {
class HepRandomEngine;
Expand Down Expand Up @@ -104,23 +105,42 @@ namespace edm {
double mass = pData->mass().value();
if (fVerbosity > 0)
edm::LogVerbatim("BeamMomentumGun") << "PDGId: " << partID << " mass: " << mass;
double xp = (xoff_ + mm2cm_ * parX_->at(ip));
double yp = (yoff_ + mm2cm_ * parY_->at(ip));
HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(xp, yp, zpos_));
double xp = (xoff_ * cm2mm_ + (-1) * parY_->at(ip)); // 90 degree rotation applied
double yp = (yoff_ * cm2mm_ + parX_->at(ip)); // 90 degree rotation applied
double zp = zpos_ * cm2mm_;
HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(xp, yp, zp));
double pxGeV = MeV2GeV_ * parPx_->at(ip);
double pyGeV = MeV2GeV_ * parPy_->at(ip);
double pzGeV = MeV2GeV_ * parPz_->at(ip);
double momRand2 = pxGeV * pxGeV + pyGeV * pyGeV + pzGeV * pzGeV;
double energy = std::sqrt(momRand2 + mass * mass);
double mom = std::sqrt(momRand2);
HepMC::FourVector pGeV(pxGeV, pyGeV, pzGeV, energy);
double ptheta = pGeV.theta();
double pphi = pGeV.phi();
double theta = CLHEP::RandFlat::shoot(engine, fMinTheta, fMaxTheta);
double phi = CLHEP::RandFlat::shoot(engine, fMinPhi, fMaxPhi);
double px = mom * sin(theta) * cos(phi);
double py = mom * sin(theta) * sin(phi);
double pz = mom * cos(theta);

if (fVerbosity > 0)
edm::LogVerbatim("BeamMomentumGun") << "px:py:pz " << px << ":" << py << ":" << pz;
if (phi > M_PI)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

there are functions in https://github.com/cms-sw/cmssw/blob/master/DataFormats/Math/interface/deltaPhi.h to do this automatically, better to use them

phi = -(2 * M_PI - phi);
double newtheta = ptheta + theta;
if (newtheta > M_PI && newtheta <= 2 * M_PI)
newtheta = 2 * M_PI - newtheta;
double newphi = pphi + phi;
if (newphi > M_PI && newphi <= 2 * M_PI)
newphi = -(2 * M_PI - newphi);
else if (newphi < -M_PI && newphi >= -2 * M_PI)
newphi = 2 * M_PI + newphi;
double px = mom * sin(newtheta) * cos(newphi);
double py = mom * sin(newtheta) * sin(newphi);
double pz = mom * cos(newtheta);

if (fVerbosity > 0) {
edm::LogVerbatim("BeamMomentumGun") << "ptheta:pphi " << ptheta << ":" << pphi << "\ntheta:phi " << theta << ":"
<< phi << "\nnewtheta:newphi " << newtheta << ":" << newphi;

edm::LogVerbatim("BeamMomentumGun")
<< "x:y:z [mm] " << xp << ":" << yp << ":" << zpos_ << "\npx:py:pz [GeV] " << px << ":" << py << ":" << pz;
}

HepMC::FourVector p(px, py, pz, energy);
HepMC::GenParticle* part = new HepMC::GenParticle(p, partID, 1);
Expand Down
90 changes: 71 additions & 19 deletions SimG4CMS/HGCalTestBeam/plugins/HGCalTBAnalyzer.cc
View file
Open in desktop
Expand Up @@ -80,7 +80,7 @@ class HGCalTBAnalyzer : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::on
const bool ifEE_, ifFH_, ifBH_, ifBeam_;
const bool doSimHits_, doDigis_, doRecHits_;
const bool doTree_, doTreeCell_;
const bool doPassive_, doPassiveEE_, doPassiveHE_, doPassiveBH_, addP_;
const bool doPassive_, doPassiveEE_, doPassiveHE_, doPassiveBH_, addP_, doBeam_;
const std::string detectorEE_, detectorFH_;
const std::string detectorBH_, detectorBeam_;
const double zFrontEE_, zFrontFH_, zFrontBH_;
Expand Down Expand Up @@ -131,6 +131,10 @@ class HGCalTBAnalyzer : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::on

double xBeam_, yBeam_, zBeam_, pBeam_;
double thetaBeam_, phiBeam_;
int nBeamMC_;
std::vector<int> pdgIdBeamMC_;
std::vector<float> xBeamMC_, yBeamMC_, zBeamMC_;
std::vector<float> pxBeamMC_, pyBeamMC_, pzBeamMC_, pBeamMC_;
};

HGCalTBAnalyzer::HGCalTBAnalyzer(const edm::ParameterSet& iConfig)
Expand All @@ -148,6 +152,7 @@ HGCalTBAnalyzer::HGCalTBAnalyzer(const edm::ParameterSet& iConfig)
doPassiveHE_(iConfig.getParameter<bool>("doPassiveHE")),
doPassiveBH_(iConfig.getParameter<bool>("doPassiveBH")),
addP_(iConfig.getParameter<bool>("addP")),
doBeam_(iConfig.getParameter<bool>("doBeam")),
detectorEE_(iConfig.getParameter<std::string>("detectorEE")),
detectorFH_(iConfig.getParameter<std::string>("detectorFH")),
detectorBH_(iConfig.getParameter<std::string>("detectorBH")),
Expand Down Expand Up @@ -213,10 +218,9 @@ HGCalTBAnalyzer::HGCalTBAnalyzer(const edm::ParameterSet& iConfig)
tmp3 = iConfig.getParameter<edm::InputTag>("recHitSrcFH");
tok_hitrFH_ = consumes<HGCRecHitCollection>(tmp3);
#ifdef EDM_ML_DEBUG
if (ifFH_) {
if (ifFH_)
edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer:: Detector " << detectorFH_ << " with tags " << tmp1 << ", " << tmp2
<< ", " << tmp3;
}
#endif
tmp1 = iConfig.getParameter<std::string>("caloHitSrcBH");
tok_hitsBH_ = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", tmp1));
Expand All @@ -242,17 +246,15 @@ HGCalTBAnalyzer::HGCalTBAnalyzer(const edm::ParameterSet& iConfig)
tok_hgcPHBeam_ = consumes<edm::PassiveHitContainer>(tmp);

#ifdef EDM_ML_DEBUG
if (ifBH_) {
if (ifBH_)
edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer:: Detector " << detectorBH_ << " with tags " << tmp1 << ", " << tmp2
<< ", " << tmp3;
}
#endif
tmp1 = iConfig.getParameter<std::string>("caloHitSrcBeam");
tok_hitsBeam_ = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", tmp1));
#ifdef EDM_ML_DEBUG
if (ifBeam_) {
if (ifBeam_)
edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer:: Detector " << detectorBeam_ << " with tags " << tmp1;
}
#endif
}

Expand Down Expand Up @@ -302,6 +304,7 @@ void HGCalTBAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descripti
desc.add<bool>("doPassiveHE", false);
desc.add<bool>("doPassiveBH", false);
desc.add<bool>("addP", false);
desc.add<bool>("doBeam", false);
desc.addUntracked<double>("gev2mip200", 57.0e-6);
desc.addUntracked<double>("gev2mip300", 85.5e-6);
desc.addUntracked<double>("stoc_smear_time_200", 10.24);
Expand Down Expand Up @@ -407,6 +410,17 @@ void HGCalTBAnalyzer::beginJob() {
tree_->Branch("pBeam", &pBeam_, "pBeam/D");
tree_->Branch("thetaBeam", &thetaBeam_, "thetaBeam/D");
tree_->Branch("phiBeam", &phiBeam_, "phiBeam/D");
if (doBeam_) {
tree_->Branch("nBeamMC", &nBeamMC_, "nBeamMC/I");
tree_->Branch("pdgIdBeamMC", &pdgIdBeamMC_);
tree_->Branch("xBeamMC", &xBeamMC_);
tree_->Branch("yBeamMC", &yBeamMC_);
tree_->Branch("zBeamMC", &zBeamMC_);
tree_->Branch("pxBeamMC", &pxBeamMC_);
tree_->Branch("pyBeamMC", &pyBeamMC_);
tree_->Branch("pzBeamMC", &pzBeamMC_);
tree_->Branch("pBeamMC", &pBeamMC_);
}
if (doTreeCell_) {
tree_->Branch("simHitCellIdEE", &simHitCellIdEE_);
tree_->Branch("simHitCellEnEE", &simHitCellEnEE_);
Expand Down Expand Up @@ -551,10 +565,12 @@ void HGCalTBAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& i
HepMC::FourVector pxyz(0, 0, 0, 0);
for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
++p, ++k) {
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HGCSim") << "Particle [" << k << "] with p " << (*p)->momentum().rho() << " theta "
<< (*p)->momentum().theta() << " phi " << (*p)->momentum().phi() << " pxyz ("
<< (*p)->momentum().px() << ", " << (*p)->momentum().py() << ", "
<< (*p)->momentum().pz() << ")";
#endif
if (addP_) {
pxyz.setPx(pxyz.px() + (*p)->momentum().px());
pxyz.setPy(pxyz.py() + (*p)->momentum().py());
Expand All @@ -565,8 +581,10 @@ void HGCalTBAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& i
}
}
hBeam_->Fill(pxyz.rho());
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HGCSim") << "Particle with p " << pxyz.rho() << " theta " << pxyz.theta() << " phi "
<< pxyz.phi();
#endif
}

// Now the Simhits
Expand Down Expand Up @@ -1122,46 +1140,80 @@ void HGCalTBAnalyzer::analyzeSimHits(int type, std::vector<PCaloHit>& hits, doub

void HGCalTBAnalyzer::analyzeSimTracks(edm::Handle<edm::SimTrackContainer> const& SimTk,
edm::Handle<edm::SimVertexContainer> const& SimVtx) {
xBeam_ = yBeam_ = zBeam_ = pBeam_ = -1000000;
thetaBeam_ = phiBeam_ = -100000;
xBeam_ = yBeam_ = zBeam_ = pBeam_ = -9999;
nBeamMC_ = thetaBeam_ = phiBeam_ = -9999;
int nParBeam = 0;
int vertIndex(-1);
if (doBeam_) {
pdgIdBeamMC_.clear();
xBeamMC_.clear();
yBeamMC_.clear();
zBeamMC_.clear();
pxBeamMC_.clear();
pyBeamMC_.clear();
pzBeamMC_.clear();
pBeamMC_.clear();
}
std::vector<float> verX, verY, verZ;
for (const auto& simVtxItr : *SimVtx) {
verX.push_back(simVtxItr.position().X());
verY.push_back(simVtxItr.position().Y());
verZ.push_back(simVtxItr.position().Z());
}
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HGCSim") << "Size of track " << SimTk->size();
#endif
HepMC::FourVector pxyz(0, 0, 0, 0);
for (const auto& simTrkItr : *SimTk) {
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HGCSim") << "Track " << simTrkItr.trackId() << " Vertex " << simTrkItr.vertIndex() << " Type "
<< simTrkItr.type() << " Charge " << simTrkItr.charge() << " momentum "
<< simTrkItr.momentum() << " " << simTrkItr.momentum().P() << " GenIndex "
<< simTrkItr.genpartIndex();
#endif
if (addP_ && !(simTrkItr.noGenpart())) {
pxyz.setPx(pxyz.px() + simTrkItr.momentum().px());
pxyz.setPy(pxyz.py() + simTrkItr.momentum().py());
pxyz.setPz(pxyz.pz() + simTrkItr.momentum().pz());
pxyz.setE(pxyz.e() + simTrkItr.momentum().e());
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HGCSim") << "Track " << simTrkItr.trackId() << " Vertex " << simTrkItr.vertIndex() << " Type "
<< simTrkItr.type() << " Charge " << simTrkItr.charge() << " px "
<< simTrkItr.momentum().px() << " py " << simTrkItr.momentum().py() << " pz "
<< simTrkItr.momentum().pz() << " P " << simTrkItr.momentum().P() << " GenIndex "
<< simTrkItr.genpartIndex();
edm::LogVerbatim("HGCSim") << "Vertex " << simTrkItr.vertIndex()
<< " position-> X: " << verX[simTrkItr.vertIndex()]
<< " Y: " << verY[simTrkItr.vertIndex()] << " Z: " << verZ[simTrkItr.vertIndex()];
#endif
}
if (doBeam_ && !(simTrkItr.noGenpart())) {
nParBeam++;
pdgIdBeamMC_.push_back(simTrkItr.type());
xBeamMC_.push_back(verX[simTrkItr.vertIndex()]);
yBeamMC_.push_back(verY[simTrkItr.vertIndex()]);
zBeamMC_.push_back(verZ[simTrkItr.vertIndex()]);
pxBeamMC_.push_back(simTrkItr.momentum().px());
pyBeamMC_.push_back(simTrkItr.momentum().py());
pzBeamMC_.push_back(simTrkItr.momentum().pz());
pBeamMC_.push_back(simTrkItr.momentum().P());
} else if (!addP_ && (vertIndex == -1)) {
pxyz = simTrkItr.momentum();
}
if (vertIndex == -1)
vertIndex = simTrkItr.vertIndex();
}
nBeamMC_ = nParBeam;
pBeam_ = pxyz.rho();
thetaBeam_ = pxyz.theta();
phiBeam_ = pxyz.phi();
if (phiBeam_ < 0)
phiBeam_ += (2 * M_PI);
if (vertIndex != -1 && vertIndex < (int)SimVtx->size()) {
edm::SimVertexContainer::const_iterator simVtxItr = SimVtx->begin();
for (int iv = 0; iv < vertIndex; iv++)
for (int iv = 0; iv < vertIndex; iv++) {
simVtxItr++;
}
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HGCSim") << "Vertex " << vertIndex << " position " << simVtxItr->position();
#endif
xBeam_ = simVtxItr->position().X();
yBeam_ = simVtxItr->position().Y();
zBeam_ = simVtxItr->position().Z();
xBeam_ = verX[0];
yBeam_ = verY[0];
zBeam_ = verZ[0];
}
}

Expand Down