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.

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

Close by gun generator #26065

Merged
merged 7 commits into from Mar 14, 2019
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
44 changes: 44 additions & 0 deletions IOMC/EventVertexGenerators/interface/PassThroughEvtVtxGenerator.h
@@ -0,0 +1,44 @@
#ifndef IOMC_EventVertexGenerators_PassThroughEvtVtxGenerator_H
#define IOMC_EventVertexGenerators_PassThroughEvtVtxGenerator_H
/*
*/

#include "IOMC/EventVertexGenerators/interface/BaseEvtVtxGenerator.h"
#include "FWCore/Framework/interface/stream/EDProducer.h"
#include "FWCore/Utilities/interface/EDGetToken.h"

#include "TMatrixD.h"

namespace HepMC {
class FourVector ;
}

namespace CLHEP {
class HepRandomEngine;
}

namespace edm {
class HepMCProduct;
}

class PassThroughEvtVtxGenerator : public BaseEvtVtxGenerator
{
public:

// ctor & dtor
explicit PassThroughEvtVtxGenerator( const edm::ParameterSet& );
~PassThroughEvtVtxGenerator() override;

void produce( edm::Event&, const edm::EventSetup&) override;

HepMC::FourVector newVertex(CLHEP::HepRandomEngine*) const override;

TMatrixD const* GetInvLorentzBoost() const override { return nullptr;};

private :

edm::EDGetTokenT<edm::HepMCProduct> sourceToken;

};

#endif
67 changes: 67 additions & 0 deletions IOMC/EventVertexGenerators/src/PassThroughEvtVtxGenerator.cc
@@ -0,0 +1,67 @@

/*
*/

#include "IOMC/EventVertexGenerators/interface/PassThroughEvtVtxGenerator.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"

#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/Utilities/interface/RandomNumberGenerator.h"

#include "FWCore/Utilities/interface/Exception.h"

#include "DataFormats/Provenance/interface/Provenance.h"
#include "FWCore/Utilities/interface/EDMException.h"

//#include "HepMC/GenEvent.h"
// #include "CLHEP/Vector/ThreeVector.h"
// #include "HepMC/SimpleVector.h"

using namespace edm;
using namespace CLHEP;
//using namespace HepMC;


PassThroughEvtVtxGenerator::PassThroughEvtVtxGenerator( const ParameterSet& pset )
: BaseEvtVtxGenerator(pset)
{
Service<RandomNumberGenerator> rng;
if ( ! rng.isAvailable()) {
throw cms::Exception("Configuration")
<< "The PassThroughEvtVtxGenerator requires the RandomNumberGeneratorService\n"
"which is not present in the configuration file. \n"
"You must add the service\n"
"in the configuration file or remove the modules that require it.";
}
sourceToken=consumes<edm::HepMCProduct>(pset.getParameter<edm::InputTag>("src"));
}

PassThroughEvtVtxGenerator::~PassThroughEvtVtxGenerator()
{
}

HepMC::FourVector PassThroughEvtVtxGenerator::newVertex(CLHEP::HepRandomEngine*) const {
return HepMC::FourVector(0.,0.,0.);
Copy link
Contributor

Choose a reason for hiding this comment

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

Here you assume the implicit time=0 of the class, one could even set it explicitly to 0

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Indeed, I could change it into

return HepMC::FourVector(0.,0.,0.,0.);

}

void PassThroughEvtVtxGenerator::produce( Event& evt, const EventSetup& )
Copy link
Contributor

Choose a reason for hiding this comment

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

what is the purpose of not using the base class as elsewhere? To avoid the dummy LorentzBoost with a copy-only method I understand

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yes, not sure if that will gain anything, but that's what I did.

{
edm::Service<edm::RandomNumberGenerator> rng;

Handle<HepMCProduct> HepUnsmearedMCEvt ;

evt.getByToken( sourceToken, HepUnsmearedMCEvt ) ;

// Copy the HepMC::GenEvent
HepMC::GenEvent* genevt = new HepMC::GenEvent(*HepUnsmearedMCEvt->GetEvent());
std::unique_ptr<edm::HepMCProduct> HepMCEvt(new edm::HepMCProduct(genevt));

evt.put(std::move(HepMCEvt)) ;

return ;
}
2 changes: 2 additions & 0 deletions IOMC/EventVertexGenerators/src/module.cc
Expand Up @@ -3,6 +3,7 @@
//#include "IOMC/EventVertexGenerators/interface/VertexGenerator.h"

#include "IOMC/EventVertexGenerators/interface/BaseEvtVtxGenerator.h"
#include "IOMC/EventVertexGenerators/interface/PassThroughEvtVtxGenerator.h"
#include "IOMC/EventVertexGenerators/interface/GaussEvtVtxGenerator.h"
#include "IOMC/EventVertexGenerators/interface/FlatEvtVtxGenerator.h"
#include "IOMC/EventVertexGenerators/interface/BeamProfileVtxGenerator.h"
Expand All @@ -18,6 +19,7 @@

//using edm::VertexGenerator;
//DEFINE_FWK_MODULE(VertexGenerator) ;
DEFINE_FWK_MODULE(PassThroughEvtVtxGenerator) ;
DEFINE_FWK_MODULE(GaussEvtVtxGenerator) ;
DEFINE_FWK_MODULE(FlatEvtVtxGenerator) ;
DEFINE_FWK_MODULE(BeamProfileVtxGenerator) ;
Expand Down
29 changes: 29 additions & 0 deletions IOMC/ParticleGuns/interface/CloseByParticleGunProducer.h
@@ -0,0 +1,29 @@
#ifndef IOMC_ParticleGun_CloseByParticleGunProducer_H
#define IOMC_ParticleGun_CloseByParticleGunProducer_H

#include "IOMC/ParticleGuns/interface/BaseFlatGunProducer.h"

namespace edm
{

class CloseByParticleGunProducer : public BaseFlatGunProducer
{

public:
CloseByParticleGunProducer(const ParameterSet &);
~CloseByParticleGunProducer() override;

private:

void produce(Event & e, const EventSetup& es) override;

protected :

// data members
double fEn,fR,fZ,fDelta;
bool fPointing = false;
std::vector<int> fPartIDs;
};
}

#endif
125 changes: 125 additions & 0 deletions IOMC/ParticleGuns/src/CloseByParticleGunProducer.cc
@@ -0,0 +1,125 @@
#include <ostream>

#include "IOMC/ParticleGuns/interface/CloseByParticleGunProducer.h"

#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"

#include "DataFormats/Math/interface/Vector3D.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#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"

using namespace edm;
using namespace std;

CloseByParticleGunProducer::CloseByParticleGunProducer(const ParameterSet& pset) :
BaseFlatGunProducer(pset)
{

ParameterSet defpset ;
ParameterSet pgun_params =
pset.getParameter<ParameterSet>("PGunParameters") ;

fEn = pgun_params.getParameter<double>("En");
fR = pgun_params.getParameter<double>("R");
fZ = pgun_params.getParameter<double>("Z");
fDelta = pgun_params.getParameter<double>("Delta");
fPartIDs = pgun_params.getParameter< vector<int> >("PartID");
fPointing = pgun_params.getParameter<bool>("Pointing");

Copy link
Contributor

Choose a reason for hiding this comment

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

do you have an example configuration fragment usable for tests?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I have, since I already privately used it extensively. Would like me to put the snippet somewhere in this package? If so, where exactly?

This is an example:

process.generator = cms.EDProducer("CloseByParticleGunProducer",
                                   PGunParameters = cms.PSet(PartID = cms.vint32(22, 22),
                                                             En = cms.double(50.),
                                                             R = cms.double(60),
                                                             Z = cms.double(320),
                                                             Delta = cms.double(2.5),
                                                             Pointing = cms.bool(True),
                                                             MaxEta = cms.double(2.9),
                                                             MaxPhi = cms.double(3.14159265359/6.),
                                                             MinEta = cms.double(1.5),
                                                             MinPhi = cms.double(-3.14159265359/6.),
                                                             ),
                                   Verbosity = cms.untracked.int32(10),
                                   psethack = cms.string('two close by particles'),
                                   AddAntiParticle = cms.bool(False),
                                   firstRun = cms.untracked.uint32(1)
                                   )

Copy link
Contributor

Choose a reason for hiding this comment

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

We normally store these default fragments in Configuration/Generator/python , where they may be used for instance by RelVals .

produces<HepMCProduct>("unsmeared");
produces<GenEventInfoProduct>();
}

CloseByParticleGunProducer::~CloseByParticleGunProducer()
{
// no need to cleanup GenEvent memory - done in HepMCProduct
}

void CloseByParticleGunProducer::produce(Event &e, const EventSetup& es)
{
edm::Service<edm::RandomNumberGenerator> rng;
CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());

if ( fVerbosity > 0 )
{
LogDebug("CloseByParticleGunProducer") << " CloseByParticleGunProducer : Begin New Event Generation" << endl ;
}
fEvt = new HepMC::GenEvent() ;

// loop over particles
//
int barcode = 1 ;
double phi = CLHEP::RandFlat::shoot(engine, -3.14159265358979323846, 3.14159265358979323846);
for (unsigned int ip=0; ip<fPartIDs.size(); ++ip, phi += fDelta/fR)
{

int PartID = fPartIDs[ip] ;
const HepPDT::ParticleData *PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID))) ;
double mass = PData->mass().value() ;
double mom = sqrt(fEn*fEn-mass*mass);
double px = 0.;
double py = 0.;
double pz = mom;
double energy = fEn;

// Compute Vertex Position
double x=fR*cos(phi);
double y=fR*sin(phi);
constexpr double c= 2.99792458e+1; // cm/ns
double timeOffset = sqrt(x*x + y*y + fZ*fZ)/c*ns*c_light;
Copy link
Contributor

Choose a reason for hiding this comment

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

@rovere , I would expect here CMS units (cm, GeV, ns)? In that case of CLHEP (mm, MeV, ns) should not be used.
Having in one expression both your constant "c" and CLHEP "c_light" seems to be a problem.

To resolve this I would put a comment line about fX, fX, fZ unit, and select one or another.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I used CMS units for the configuration parameters, namely cm, GeV and ns. Then I compute the timing in ns using CMS units and the "correct" c value and convert it back to whatever CLHEP/hepmc is expecting. Finally, I put back also x,y,fZ back to units that are digested by HepMC.
I have to admit that this mix of units is, to put it mildly, quite ugly and confusing, I agree, but it's not coming from this code but rather imposed from externals.
By the way, I think you will find the same kind of "units massaging" also in other Generators.

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

HepMC::FourVector p(px,py,pz,energy) ;
// If we are requested to be pointing to (0,0,0), correct the momentum direction
if (fPointing) {
math::XYZVector direction(x,y,fZ);
math::XYZVector momentum = direction.unit() * mom;
p.setX(momentum.x());
p.setY(momentum.y());
p.setZ(momentum.z());
}
HepMC::GenParticle* Part = new HepMC::GenParticle(p,PartID,1);
Part->suggest_barcode( barcode );
barcode++;

Vtx->add_particle_out(Part);

if (fVerbosity > 0) {
Vtx->print();
Part->print();
}
fEvt->add_vertex(Vtx);
}


fEvt->set_event_number(e.id().event());
fEvt->set_signal_process_id(20);

if ( fVerbosity > 0 )
{
fEvt->print();
}

unique_ptr<HepMCProduct> BProduct(new HepMCProduct());
BProduct->addHepMCData( fEvt );
e.put(std::move(BProduct), "unsmeared");

unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(fEvt));
e.put(std::move(genEventInfo));

if ( fVerbosity > 0 )
{
LogDebug("CloseByParticleGunProducer") << " CloseByParticleGunProducer : Event Generation Done " << endl;
}
}

7 changes: 5 additions & 2 deletions IOMC/ParticleGuns/src/SealModule.cc
Expand Up @@ -19,12 +19,13 @@
#include "IOMC/ParticleGuns/interface/RandomtXiGunProducer.h"
#include "IOMC/ParticleGuns/interface/FlatRandomPtAndDxyGunProducer.h"
#include "IOMC/ParticleGuns/interface/RandomMultiParticlePGunProducer.h"
#include "IOMC/ParticleGuns/interface/CloseByParticleGunProducer.h"


// particle gun prototypes
//


/*
using edm::FlatEGunASCIIWriter;
DEFINE_FWK_MODULE(FlatEGunASCIIWriter);
Expand Down Expand Up @@ -58,3 +59,5 @@ using edm::FlatRandomPtAndDxyGunProducer;
DEFINE_FWK_MODULE(FlatRandomPtAndDxyGunProducer);
using edm::RandomMultiParticlePGunProducer;
DEFINE_FWK_MODULE(RandomMultiParticlePGunProducer);
using edm::CloseByParticleGunProducer;
DEFINE_FWK_MODULE(CloseByParticleGunProducer);