Skip to content

Commit

Permalink
Add more configuration options (backport)
Browse files Browse the repository at this point in the history
  • Loading branch information
Oz Amram committed Jul 30, 2024
1 parent b770785 commit 938b266
Show file tree
Hide file tree
Showing 2 changed files with 101 additions and 45 deletions.
9 changes: 6 additions & 3 deletions IOMC/ParticleGuns/interface/CloseByParticleGunProducer.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,11 +25,14 @@ namespace edm {

protected:
// data members
bool fControlledByEta;
double fEnMin, fEnMax, fEtaMin, fEtaMax, fRMin, fRMax, fZMin, fZMax, fDelta, fPhiMin, fPhiMax, fTMin, fTMax,
bool fControlledByEta, fControlledByREta;
double fVarMin, fVarMax, fEtaMin, fEtaMax, fRMin, fRMax, fZMin, fZMax, fDelta, fPhiMin, fPhiMax, fTMin, fTMax,
fOffsetFirst;
double log_fVarMin = 0., log_fVarMax = 0.;
int fNParticles;
bool fMaxEnSpread = false;
bool fLogSpacedVar = false;
bool fMaxVarSpread = false;
bool fFlatPtGeneration = false;
bool fPointing = false;
bool fOverlapping = false;
bool fRandomShoot = false;
Expand Down
137 changes: 95 additions & 42 deletions IOMC/ParticleGuns/src/CloseByParticleGunProducer.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,10 +14,10 @@
#include "FWCore/Utilities/interface/RandomNumberGenerator.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "CLHEP/Random/RandFlat.h"
#include "CLHEP/Units/GlobalSystemOfUnits.h"
#include "CLHEP/Units/GlobalPhysicalConstants.h"
#include "CLHEP/Random/RandFlat.h"
#include <CLHEP/Random/RandFlat.h>
#include <CLHEP/Units/SystemOfUnits.h>
#include <CLHEP/Units/GlobalPhysicalConstants.h>
#include <CLHEP/Random/RandFlat.h>

using namespace edm;
using namespace std;
Expand All @@ -26,31 +26,53 @@ CloseByParticleGunProducer::CloseByParticleGunProducer(const ParameterSet& pset)
: BaseFlatGunProducer(pset), m_fieldToken(esConsumes()) {
ParameterSet pgun_params = pset.getParameter<ParameterSet>("PGunParameters");
fControlledByEta = pgun_params.getParameter<bool>("ControlledByEta");
fEnMax = pgun_params.getParameter<double>("EnMax");
fEnMin = pgun_params.getParameter<double>("EnMin");
if (fEnMin < 1)
LogError("CloseByParticleGunProducer") << " Please choose a minimum energy greater than 1 GeV, otherwise time "
"information may be invalid or not reliable";

fMaxEnSpread = pgun_params.getParameter<bool>("MaxEnSpread");
if (fControlledByEta) {
fControlledByREta = pgun_params.getParameter<bool>("ControlledByREta");
if (fControlledByEta and fControlledByREta)
throw cms::Exception("CloseByParticleGunProducer")
<< " Conflicting configuration, cannot have both ControlledByEta and ControlledByREta ";

fVarMax = pgun_params.getParameter<double>("VarMax");
fVarMin = pgun_params.getParameter<double>("VarMin");
fMaxVarSpread = pgun_params.getParameter<bool>("MaxVarSpread");
fLogSpacedVar = pgun_params.getParameter<bool>("LogSpacedVar");
fFlatPtGeneration = pgun_params.getParameter<bool>("FlatPtGeneration");
if (fVarMin < 1 && !fFlatPtGeneration)
throw cms::Exception("CloseByParticleGunProducer")
<< " Please choose a minimum energy greater than 1 GeV, otherwise time "
"information may be invalid or not reliable";
if (fVarMin < 0 && fLogSpacedVar)
throw cms::Exception("CloseByParticleGunProducer") << " Minimum energy must be greater than zero for log spacing";
else {
log_fVarMin = std::log(fVarMin);
log_fVarMax = std::log(fVarMax);
}

if (fControlledByEta || fControlledByREta) {
fEtaMax = pgun_params.getParameter<double>("MaxEta");
fEtaMin = pgun_params.getParameter<double>("MinEta");
if (fEtaMax <= fEtaMin)
LogError("CloseByParticleGunProducer") << " Please fix MinEta and MaxEta values in the configuration";
throw cms::Exception("CloseByParticleGunProducer") << " Please fix MinEta and MaxEta values in the configuration";
} else {
fRMax = pgun_params.getParameter<double>("RMax");
fRMin = pgun_params.getParameter<double>("RMin");
if (fRMax <= fRMin)
LogError("CloseByParticleGunProducer") << " Please fix RMin and RMax values in the configuration";
throw cms::Exception("CloseByParticleGunProducer") << " Please fix RMin and RMax values in the configuration";
}
if (!fControlledByREta) {
fZMax = pgun_params.getParameter<double>("ZMax");
fZMin = pgun_params.getParameter<double>("ZMin");

if (fZMax <= fZMin)
throw cms::Exception("CloseByParticleGunProducer") << " Please fix ZMin and ZMax values in the configuration";
}
fZMax = pgun_params.getParameter<double>("ZMax");
fZMin = pgun_params.getParameter<double>("ZMin");
fDelta = pgun_params.getParameter<double>("Delta");
fPhiMin = pgun_params.getParameter<double>("MinPhi");
fPhiMax = pgun_params.getParameter<double>("MaxPhi");
fPointing = pgun_params.getParameter<bool>("Pointing");
fOverlapping = pgun_params.getParameter<bool>("Overlapping");
if (fFlatPtGeneration && !fPointing)
throw cms::Exception("CloseByParticleGunProducer")
<< " Can't generate non pointing FlatPt samples; please disable FlatPt generation or generate pointing sample";
fRandomShoot = pgun_params.getParameter<bool>("RandomShoot");
fNParticles = pgun_params.getParameter<int>("NParticles");
fPartIDs = pgun_params.getParameter<vector<int>>("PartID");
Expand All @@ -60,7 +82,7 @@ CloseByParticleGunProducer::CloseByParticleGunProducer(const ParameterSet& pset)
fTMax = pgun_params.getParameter<double>("TMax");
fTMin = pgun_params.getParameter<double>("TMin");
if (fTMax <= fTMin)
LogError("CloseByParticleGunProducer") << " Please fix TMin and TMax values in the configuration";
throw cms::Exception("CloseByParticleGunProducer") << " Please fix TMin and TMax values in the configuration";
// set a fixed time offset for the particles
fOffsetFirst = pgun_params.getParameter<double>("OffsetFirst");

Expand All @@ -78,10 +100,13 @@ void CloseByParticleGunProducer::fillDescriptions(ConfigurationDescriptions& des
{
edm::ParameterSetDescription psd0;
psd0.add<bool>("ControlledByEta", false);
psd0.add<bool>("ControlledByREta", false);
psd0.add<double>("Delta", 10);
psd0.add<double>("EnMax", 200.0);
psd0.add<double>("EnMin", 25.0);
psd0.add<bool>("MaxEnSpread", false);
psd0.add<double>("VarMax", 200.0);
psd0.add<double>("VarMin", 25.0);
psd0.add<bool>("MaxVarSpread", false);
psd0.add<bool>("LogSpacedVar", false);
psd0.add<bool>("FlatPtGeneration", false);
psd0.add<double>("MaxEta", 2.7);
psd0.add<double>("MaxPhi", 3.14159265359);
psd0.add<double>("MinEta", 1.7);
Expand Down Expand Up @@ -125,15 +150,24 @@ void CloseByParticleGunProducer::produce(Event& e, const EventSetup& es) {
unsigned int numParticles = fRandomShoot ? CLHEP::RandFlat::shoot(engine, 1, fNParticles) : fNParticles;

double phi = CLHEP::RandFlat::shoot(engine, fPhiMin, fPhiMax);
double fZ = CLHEP::RandFlat::shoot(engine, fZMin, fZMax);
double fR;
double fZ;
double fR, fEta;
double fT;

if (!fControlledByEta) {
fR = CLHEP::RandFlat::shoot(engine, fRMin, fRMax);
if (!fControlledByREta) {
fZ = CLHEP::RandFlat::shoot(engine, fZMin, fZMax);

if (!fControlledByEta) {
fR = CLHEP::RandFlat::shoot(engine, fRMin, fRMax);
fEta = asinh(fZ / fR);
} else {
fEta = CLHEP::RandFlat::shoot(engine, fEtaMin, fEtaMax);
fR = (fZ / sinh(fEta));
}
} else {
double fEta = CLHEP::RandFlat::shoot(engine, fEtaMin, fEtaMax);
fR = (fZ / sinh(fEta));
fR = CLHEP::RandFlat::shoot(engine, fRMin, fRMax);
fEta = CLHEP::RandFlat::shoot(engine, fEtaMin, fEtaMax);
fZ = sinh(fEta) / fR;
}

if (fUseDeltaT) {
Expand All @@ -153,26 +187,44 @@ void CloseByParticleGunProducer::produce(Event& e, const EventSetup& es) {
} else
phi += fDelta / fR;

double fEn;
if (numParticles > 1 && fMaxEnSpread)
fEn = fEnMin + ip * (fEnMax - fEnMin) / (numParticles - 1);
double fVar;
if (numParticles > 1 && fMaxVarSpread)
fVar = fVarMin + ip * (fVarMax - fVarMin) / (numParticles - 1);
else if (fLogSpacedVar) {
double fVar_log = CLHEP::RandFlat::shoot(engine, log_fVarMin, log_fVarMax);
fVar = std::exp(fVar_log);
}

else
fEn = CLHEP::RandFlat::shoot(engine, fEnMin, fEnMax);
fVar = CLHEP::RandFlat::shoot(engine, fVarMin, fVarMax);

int partIdx = CLHEP::RandFlat::shoot(engine, 0, fPartIDs.size());
int PartID = fPartIDs[partIdx];
const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID)));
double mass = PData->mass().value();
double mom2 = fEn * fEn - mass * mass;
double mom = 0.;
if (mom2 > 0.) {
mom = sqrt(mom2);
}
double px = 0.;
double py = 0.;
double pz = mom;
double energy = fEn;

double mom, px, py, pz;
double energy;

if (!fFlatPtGeneration) {
double mom2 = fVar * fVar - mass * mass;
mom = 0.;
if (mom2 > 0.) {
mom = sqrt(mom2);
}
px = 0.;
py = 0.;
pz = mom;
energy = fVar;
} else {
double theta = 2. * atan(exp(-fEta));
mom = fVar / sin(theta);
px = fVar * cos(phi);
py = fVar * sin(phi);
pz = mom * cos(theta);
double energy2 = mom * mom + mass * mass;
energy = sqrt(energy2);
}
// Compute Vertex Position
double x = fR * cos(phi);
double y = fR * sin(phi);
Expand All @@ -189,7 +241,7 @@ void CloseByParticleGunProducer::produce(Event& e, const EventSetup& es) {

// compute correct path assuming uniform magnetic field in CMS
double pathLength = 0.;
const double speed = p.pz() / p.e() * c_light / cm;
const double speed = p.pz() / p.e() * c_light / CLHEP::cm;
if (PData->charge()) {
// Radius [cm] = P[GeV/c] * 10^9 / (c[mm/ns] * 10^6 * q[C] * B[T]) * 100[cm/m]
const double radius = std::sqrt(p.px() * p.px() + p.py() * p.py()) * std::pow(10, 5) /
Expand All @@ -202,9 +254,10 @@ void CloseByParticleGunProducer::produce(Event& e, const EventSetup& es) {

// if not pointing time doesn't mean a lot, keep the old way
const double pathTime = fPointing ? (pathLength / speed) : (std::sqrt(x * x + y * y + fZ * fZ) / speed);
double timeOffset = fOffsetFirst + (pathTime + ip * fT) * ns * c_light;
double timeOffset = fOffsetFirst + (pathTime + ip * fT) * CLHEP::ns * c_light;

HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(x * cm, y * cm, fZ * cm, timeOffset));
HepMC::GenVertex* Vtx =
new HepMC::GenVertex(HepMC::FourVector(x * CLHEP::cm, y * CLHEP::cm, fZ * CLHEP::cm, timeOffset));

HepMC::GenParticle* Part = new HepMC::GenParticle(p, PartID, 1);
Part->suggest_barcode(barcode);
Expand Down

0 comments on commit 938b266

Please sign in to comment.