2 #include "HepMC/GenEvent.h"
3 #include "HepMC/GenVertex.h"
4 #include "HepMC/GenParticle.h"
24 using namespace HepPDT;
33 : nSimTracks(0), nSimVertices(0), nGenParticles(0), nChargedParticleTracks(0), initialSize(5000) {
90 void FBaseSimEvent::fill(
const std::vector<SimTrack>& simTracks,
const std::vector<SimVertex>& simVertices) {
96 unsigned nVtx = simVertices.size();
97 unsigned nTks = simTracks.size();
104 std::vector<int> myVertices(nVtx, -1);
105 std::vector<int> myTracks(nTks, -1);
110 std::map<unsigned, unsigned> geantToIndex;
111 for (
unsigned it = 0; it < simTracks.size(); ++it) {
112 geantToIndex[simTracks[it].trackId()] = it;
140 for (
unsigned trackId = 0; trackId < nTks; ++trackId) {
155 std::map<unsigned, unsigned>::iterator
association = geantToIndex.find(motherGeantId);
156 if (association != geantToIndex.end())
157 motherId = association->second;
159 int originId = motherId == -1 ? -1 : myTracks[motherId];
173 if (myVertices[vertexId] == -1)
181 int motherType = motherId == -1 ? 0 : simTracks[motherId].type();
183 bool notBremInDetector = (
abs(motherType) != 11 &&
std::abs(motherType) != 13) || motherType != track.
type() ||
186 if (notBremInDetector) {
200 if (myTracks[trackId] >= 0) {
205 myTracks[trackId] = myTracks[motherId];
206 if (myTracks[trackId] >= 0) {
214 for (
unsigned vertexId = 0; vertexId < nVtx; ++vertexId) {
216 if (myVertices[vertexId] != -1)
228 std::map<unsigned, unsigned>::iterator
association = geantToIndex.find(motherGeantId);
229 if (association != geantToIndex.end())
230 motherId = association->second;
232 int originId = motherId == -1 ? -1 : myTracks[motherId];
250 for (
int fsimi = 0; fsimi < (int)
nTracks(); ++fsimi) {
252 double trackerSurfaceTime =
313 int genEventSize = myGenEvent.particles_size();
314 std::vector<int> myGenVertices(genEventSize, static_cast<int>(0));
317 if (myGenEvent.particles_empty())
324 HepMC::GenVertex*
primaryVertex = *(myGenEvent.vertices_begin());
328 primaryVertex->position().y() / 10.,
329 primaryVertex->position().z() / 10.,
330 primaryVertex->position().t() / 10.);
335 int initialBarcode = 0;
336 if (myGenEvent.particles_begin() != myGenEvent.particles_end()) {
337 initialBarcode = (*myGenEvent.particles_begin())->barcode();
341 for (
auto piter = myGenEvent.particles_begin(); piter != myGenEvent.particles_end(); ++piter) {
357 HepMC::GenVertex* productionVertex = p->production_vertex();
358 if (productionVertex) {
359 unsigned productionMother = productionVertex->particles_in_size();
360 if (productionMother) {
361 unsigned motherId = (*(productionVertex->particles_in_const_begin()))->pdg_id();
362 if (motherId < 1000000)
363 productionVertexPosition =
XYZTLorentzVector(productionVertex->position().x() / 10.,
364 productionVertex->position().y() / 10.,
365 productionVertex->position().z() / 10.,
366 productionVertex->position().t() / 10.);
372 int abspdgId =
std::abs(p->pdg_id());
373 HepMC::GenVertex* endVertex = p->end_vertex();
377 bool testStable = p->status() % 1000 == 1;
380 if (p->status() == 2 && abspdgId < 1000000) {
383 endVertex->position().y() / 10.,
384 endVertex->position().z() / 10.,
385 endVertex->position().t() / 10.);
393 bool testDaugh =
false;
394 if (!testStable && p->status() == 2 && endVertex && endVertex->particles_out_size()) {
395 HepMC::GenVertex::particles_out_const_iterator firstDaughterIt = endVertex->particles_out_const_begin();
396 HepMC::GenVertex::particles_out_const_iterator lastDaughterIt = endVertex->particles_out_const_end();
397 for (; firstDaughterIt != lastDaughterIt; ++firstDaughterIt) {
399 if (daugh->status() % 1000 == 1) {
401 if (abspdgId == 11 || abspdgId == 13) {
404 endVertex->position().y() / 10.,
405 endVertex->position().z() / 10.,
406 endVertex->position().t() / 10.);
421 if (!testStable && !testDaugh && p->production_vertex()) {
422 XYZTLorentzVector productionVertexPosition(p->production_vertex()->position().x() / 10.,
423 p->production_vertex()->position().y() / 10.,
424 p->production_vertex()->position().z() / 10.,
425 p->production_vertex()->position().t() / 10.);
426 dist = (primaryVertexPosition - productionVertexPosition).Vect().Mag2();
428 bool testDecay = (dist > 1
e-8) ?
true :
false;
431 if (testStable || testDaugh || testDecay) {
437 int motherBarcode = p->production_vertex() && p->production_vertex()->particles_in_const_begin() !=
438 p->production_vertex()->particles_in_const_end()
439 ? (*(p->production_vertex()->particles_in_const_begin()))->barcode()
442 int originVertex = motherBarcode && myGenVertices[motherBarcode - initialBarcode]
443 ? myGenVertices[motherBarcode - initialBarcode]
446 XYZTLorentzVector momentum(p->momentum().px(), p->momentum().py(), p->momentum().pz(), p->momentum().e());
451 int theTrack = testStable && p->end_vertex() ?
463 (testStable && p->end_vertex() && !p->end_vertex()->particles_out_size())
470 p->end_vertex()->position().y() / 10.,
471 p->end_vertex()->position().z() / 10.,
472 p->end_vertex()->position().t() / 10.);
477 myGenVertices[p->barcode() - initialBarcode] = theVertex;
499 if (!
vertex(iv).noParent()) {
510 (*theSimTracks)[trackId] =
521 FSimTrack(p, iv, ig, trackId,
this);
544 (*theSimVertices)[vertexId] =
FSimVertex(v, im, vertexId,
this);
552 std::cout <<
"Id Gen Name eta phi pT E Vtx1 "
554 <<
"Moth Vtx2 eta phi R Z Da1 Da2 Ecal?" << std::endl;
556 for (HepMC::GenEvent::particle_const_iterator piter = myGenEvent.particles_begin();
557 piter != myGenEvent.particles_end();
561 int partId = p->pdg_id();
564 if (
pdt->particle(ParticleID(partId)) !=
nullptr) {
565 name = (
pdt->particle(ParticleID(partId)))->
name();
570 XYZTLorentzVector momentum1(p->momentum().px(), p->momentum().py(), p->momentum().pz(), p->momentum().e());
574 if (!p->production_vertex())
577 XYZVector vertex1(p->production_vertex()->position().x() / 10.,
578 p->production_vertex()->position().y() / 10.,
579 p->production_vertex()->position().z() / 10.);
580 vertexId1 = p->production_vertex()->barcode();
582 std::cout.setf(std::ios::fixed, std::ios::floatfield);
583 std::cout.setf(std::ios::right, std::ios::adjustfield);
587 for (
unsigned int k = 0;
k < 11 - name.length() &&
k < 12;
k++)
590 double eta = momentum1.eta();
595 std::cout << std::setw(6) << std::setprecision(2) << eta <<
" " << std::setw(6) << std::setprecision(2)
596 << momentum1.phi() <<
" " << std::setw(7) << std::setprecision(2) << momentum1.pt() <<
" " << std::setw(7)
597 << std::setprecision(2) << momentum1.e() <<
" " << std::setw(4) << vertexId1 <<
" " << std::setw(6)
598 << std::setprecision(1) << vertex1.x() <<
" " << std::setw(6) << std::setprecision(1) << vertex1.y()
599 <<
" " << std::setw(6) << std::setprecision(1) << vertex1.z() <<
" ";
601 const HepMC::GenParticle* mother = *(p->production_vertex()->particles_in_const_begin());
604 std::cout << std::setw(4) << mother->barcode() <<
" ";
608 if (p->end_vertex()) {
610 p->end_vertex()->position().y() / 10.,
611 p->end_vertex()->position().z() / 10.,
612 p->end_vertex()->position().t() / 10.);
613 int vertexId2 = p->end_vertex()->barcode();
615 std::vector<const HepMC::GenParticle*> children;
616 HepMC::GenVertex::particles_out_const_iterator firstDaughterIt = p->end_vertex()->particles_out_const_begin();
617 HepMC::GenVertex::particles_out_const_iterator lastDaughterIt = p->end_vertex()->particles_out_const_end();
618 for (; firstDaughterIt != lastDaughterIt; ++firstDaughterIt) {
619 children.push_back(*firstDaughterIt);
622 std::cout << std::setw(4) << vertexId2 <<
" " << std::setw(6) << std::setprecision(2) << vertex2.eta() <<
" "
623 << std::setw(6) << std::setprecision(2) << vertex2.phi() <<
" " << std::setw(5) << std::setprecision(1)
624 << vertex2.pt() <<
" " << std::setw(6) << std::setprecision(1) << vertex2.z() <<
" ";
625 for (
unsigned id = 0;
id < children.size(); ++
id)
626 std::cout << std::setw(4) << children[
id]->barcode() <<
" ";
633 std::cout <<
" Id Gen Name eta phi pT E Vtx1 "
635 <<
"Moth Vtx2 eta phi R Z Daughters Ecal?" << std::endl;
const ParticleDataTable * pdt
double lateVertexPosition
const math::XYZVectorD & trackerSurfacePosition() const
int addSimVertex(const XYZTLorentzVector &decayVertex, int im=-1, FSimVertexType::VertexType type=FSimVertexType::ANY)
Add a new vertex to the Event and to the various lists.
float charge() const
charge
bool acceptParticle(const RawParticle &p) const
bool propagateToPreshowerLayer1(bool first=true)
const HepMC::GenParticle * embdGenpart(int i) const
return MC track with a given id
uint16_t *__restrict__ id
HepPDT::ParticleDataTable ParticleDataTable
std::vector< FSimTrack > * theSimTracks
std::vector< unsigned > * theChargedTracks
FBaseSimEvent(const edm::ParameterSet &kine)
Default constructor.
const XYZTLorentzVector & momentum() const
Temporary (until move of SimTrack to Mathcore) - No! Actually very useful.
unsigned int theTrackSize
unsigned int theVertexSize
KineParticleFilter * myFilter
The particle filter.
int addSimTrack(const RawParticle *p, int iv, int ig=-1, const HepMC::GenVertex *ev=nullptr)
Add a new track to the Event and to the various lists.
const HepPDT::ParticleDataTable * theTable() const
Get the pointer to the particle data table.
void setEcal(const RawParticle &pp, int success)
Set the ecal variables.
int getSuccess() const
Has propagation been performed and was barrel or endcap reached ?
void addParticles(const HepMC::GenEvent &hev)
Add the particles and their vertices to the list.
int pid() const
get the HEP particle ID number
void setLayer2(const RawParticle &pp, int success)
Set the preshower layer2 variables.
void clear()
clear the FBaseSimEvent content before the next event
void setMagneticField(double b)
Set the magnetic field.
RawParticle const & particle() const
The particle being propagated.
bool acceptVertex(const math::XYZTLorentzVector &p) const
bool notYetToEndVertex(const XYZTLorentzVector &pos) const
Compare the end vertex position with another position.
FSimVertex & vertex(int id) const
Return vertex with given Id.
unsigned int nVertices() const
Number of vertices.
RawParticle makeParticle(HepPDT::ParticleDataTable const *, int id, const math::XYZTLorentzVector &p)
FSimVertexType & vertexType(int id) const
Return vertex with given Id.
bool propagateToVFcalEntrance(bool first=true)
void setLayer1(const RawParticle &pp, int success)
Set the preshower layer1 variables.
const XYZTLorentzVector & momentum() const
the momentum fourvector
double t() const
vertex time
std::tuple< layerClusterToCaloParticle, caloParticleToLayerCluster > association
~FBaseSimEvent()
usual virtual destructor
const math::XYZTLorentzVector & position() const
Temporary (until CMSSW moves to Mathcore) - No ! Actually very useful.
unsigned int nTracks() const
Number of tracks.
bool propagateToHcalExit(bool first=true)
unsigned int nSimVertices
int genpartIndex() const
index of the corresponding Generator particle in the Event container (-1 if no Genpart) ...
double cos2ThetaV() const
A FSimVertexType hold the information on the vertex origine.
Abs< T >::type abs(const T &t)
const math::XYZTLorentzVectorD & position() const
std::vector< FSimVertexType > FSimVertexTypeCollection
collection of FSimVertexType objects
unsigned int nGenParticles
std::vector< FSimVertex > * theSimVertices
bool propagateToEcalEntrance(bool first=true)
void fill(const HepMC::GenEvent &hev)
fill the FBaseSimEvent from the current HepMC::GenEvent
int vertIndex() const
index of the vertex in the Event container (-1 if no vertex)
const XYZTLorentzVector & vertex() const
the vertex fourvector
unsigned int theChargedSize
const math::XYZTLorentzVectorD & trackerSurfaceMomentum() const
void addChargedTrack(int id)
Add an id in the vector of charged tracks id's.
void setHcalExit(const RawParticle &pp, int success)
Set the hcal exit variables.
unsigned int nChargedParticleTracks
int type() const
particle type (HEP PDT convension)
bool propagateToHcalEntrance(bool first=true)
void initializePdt(const HepPDT::ParticleDataTable *aPdt)
Initialize the particle data table.
void setEndVertex(int endv)
Set the end vertex.
unsigned int nGenParts() const
Number of generator particles.
const math::XYZTLorentzVectorD & momentum() const
static int position[264][3]
double mass(int pdgID, const HepPDT::ParticleDataTable *pdt)
void addDaughter(int i)
Add a RecHit for a track on a layer.
void setHO(const RawParticle &pp, int success)
Set the ho variables.
int chargedTrack(int id) const
return "reconstructed" charged tracks index.
void setVFcal(const RawParticle &pp, int success)
Set the hcal variables.
void setHcal(const RawParticle &pp, int success)
Set the hcal variables.
const FSimVertex vertex() const
Origin vertex.
void printMCTruth(const HepMC::GenEvent &hev)
print the original MCTruth event
std::vector< HepMC::GenParticle * > * theGenParticles
math::XYZVector XYZVector
bool propagateToHOLayer(bool first=true)
bool propagateToPreshowerLayer2(bool first=true)
math::XYZTLorentzVector XYZTLorentzVector
void print() const
print the FBaseSimEvent in an intelligible way
FSimTrack & track(int id) const
Return track with given Id.
FSimVertexTypeCollection * theFSimVerticesType