CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
GenParticles2HepMCConverter.cc
Go to the documentation of this file.
3 
9 
13 
21 
22 #include <iostream>
23 #include <map>
24 
25 using namespace std;
26 
28 public:
31 
32  void beginRun(edm::Run const& iRun, edm::EventSetup const&) override;
33  void produce(edm::Event& event, const edm::EventSetup& eventSetup) override;
34 
35 private:
40 
41  std::vector<int> signalParticlePdgIds_;
42  const double cmEnergy_;
43  HepMC::GenCrossSection xsec_;
44 
45 private:
46  inline HepMC::FourVector FourVector(const reco::Candidate::Point& point) {
47  return HepMC::FourVector(10 * point.x(), 10 * point.y(), 10 * point.z(), 0);
48  };
49 
50  inline HepMC::FourVector FourVector(const reco::Candidate::LorentzVector& lvec) {
51  // Avoid negative mass, set minimum m^2 = 0
52  return HepMC::FourVector(lvec.px(), lvec.py(), lvec.pz(), std::hypot(lvec.P(), std::max(0., lvec.mass())));
53  };
54 };
55 
57  // dummy value to set incident proton pz for particle gun samples
58  : cmEnergy_(pset.getUntrackedParameter<double>("cmEnergy", 13000)) {
59  genParticlesToken_ = consumes<reco::CandidateView>(pset.getParameter<edm::InputTag>("genParticles"));
60  genEventInfoToken_ = consumes<GenEventInfoProduct>(pset.getParameter<edm::InputTag>("genEventInfo"));
61  genRunInfoToken_ = consumes<GenRunInfoProduct, edm::InRun>(pset.getParameter<edm::InputTag>("genEventInfo"));
62  signalParticlePdgIds_ = pset.getParameter<std::vector<int>>("signalParticlePdgIds");
63 
64  produces<edm::HepMCProduct>("unsmeared");
65 }
66 
68  edm::Handle<GenRunInfoProduct> genRunInfoHandle;
69  iRun.getByToken(genRunInfoToken_, genRunInfoHandle);
70 
71  xsec_.set_cross_section(genRunInfoHandle->internalXSec().value(), genRunInfoHandle->internalXSec().error());
72 }
73 
75  edm::Handle<reco::CandidateView> genParticlesHandle;
76  event.getByToken(genParticlesToken_, genParticlesHandle);
77 
78  edm::Handle<GenEventInfoProduct> genEventInfoHandle;
79  event.getByToken(genEventInfoToken_, genEventInfoHandle);
80 
81  eventSetup.getData(pTable_);
82 
83  HepMC::GenEvent* hepmc_event = new HepMC::GenEvent();
84  hepmc_event->set_event_number(event.id().event());
85  hepmc_event->set_signal_process_id(genEventInfoHandle->signalProcessID());
86  hepmc_event->set_event_scale(genEventInfoHandle->qScale());
87  hepmc_event->set_alphaQED(genEventInfoHandle->alphaQED());
88  hepmc_event->set_alphaQCD(genEventInfoHandle->alphaQCD());
89 
90  hepmc_event->weights() = genEventInfoHandle->weights();
91 
92  hepmc_event->set_cross_section(xsec_);
93 
94  // Set PDF
95  const gen::PdfInfo* pdf = genEventInfoHandle->pdf();
96  if (pdf != nullptr) {
97  const int pdf_id1 = pdf->id.first, pdf_id2 = pdf->id.second;
98  const double pdf_x1 = pdf->x.first, pdf_x2 = pdf->x.second;
99  const double pdf_scalePDF = pdf->scalePDF;
100  const double pdf_xPDF1 = pdf->xPDF.first, pdf_xPDF2 = pdf->xPDF.second;
101  HepMC::PdfInfo hepmc_pdfInfo(pdf_id1, pdf_id2, pdf_x1, pdf_x2, pdf_scalePDF, pdf_xPDF1, pdf_xPDF2);
102  hepmc_event->set_pdf_info(hepmc_pdfInfo);
103  }
104 
105  // Prepare list of HepMC::GenParticles
106  std::map<const reco::Candidate*, HepMC::GenParticle*> genCandToHepMCMap;
107  HepMC::GenParticle *hepmc_parton1 = nullptr, *hepmc_parton2 = nullptr;
108  std::vector<HepMC::GenParticle*> hepmc_particles;
109  const reco::Candidate *parton1 = nullptr, *parton2 = nullptr;
110  for (unsigned int i = 0, n = genParticlesHandle->size(); i < n; ++i) {
111  const reco::Candidate* p = &genParticlesHandle->at(i);
112  HepMC::GenParticle* hepmc_particle = new HepMC::GenParticle(FourVector(p->p4()), p->pdgId(), p->status());
113  hepmc_particle->suggest_barcode(i + 1);
114 
115  // Assign particle's generated mass from the standard particle data table
116  double particleMass;
117  if (pTable_->particle(p->pdgId()))
118  particleMass = pTable_->particle(p->pdgId())->mass();
119  else
120  particleMass = p->mass();
121 
122  hepmc_particle->set_generated_mass(particleMass);
123 
124  hepmc_particles.push_back(hepmc_particle);
125  genCandToHepMCMap[p] = hepmc_particle;
126 
127  // Find incident proton pair
128  if (p->mother() == nullptr and std::abs(p->eta()) > 5 and std::abs(p->pz()) > 1000) {
129  if (!parton1 and p->pz() > 0) {
130  parton1 = p;
131  hepmc_parton1 = hepmc_particle;
132  } else if (!parton2 and p->pz() < 0) {
133  parton2 = p;
134  hepmc_parton2 = hepmc_particle;
135  }
136  }
137  }
138 
139  HepMC::GenVertex* vertex1 = nullptr;
140  HepMC::GenVertex* vertex2 = nullptr;
141  if (parton1 == nullptr and parton2 == nullptr) {
142  // Particle gun samples do not have incident partons. Put dummy incident particle and prod vertex
143  // Note: leave parton1 and parton2 as nullptr since it is not used anymore after creating hepmc_parton1 and 2
144  const reco::Candidate::LorentzVector nullP4(0, 0, 0, 0);
145  const reco::Candidate::LorentzVector beamP4(0, 0, cmEnergy_ / 2, cmEnergy_ / 2);
146  vertex1 = new HepMC::GenVertex(FourVector(nullP4));
147  vertex2 = new HepMC::GenVertex(FourVector(nullP4));
148  hepmc_parton1 = new HepMC::GenParticle(FourVector(+beamP4), 2212, 4);
149  hepmc_parton2 = new HepMC::GenParticle(FourVector(-beamP4), 2212, 4);
150  } else {
151  // Put incident beam particles : proton -> parton vertex
152  vertex1 = new HepMC::GenVertex(FourVector(parton1->vertex()));
153  vertex2 = new HepMC::GenVertex(FourVector(parton2->vertex()));
154  }
155  hepmc_event->add_vertex(vertex1);
156  hepmc_event->add_vertex(vertex2);
157  vertex1->add_particle_in(hepmc_parton1);
158  vertex2->add_particle_in(hepmc_parton2);
159  //hepmc_event->set_beam_particles(hepmc_parton1, hepmc_parton2);
160 
161  // Prepare vertex list
162  typedef std::map<const reco::Candidate*, HepMC::GenVertex*> ParticleToVertexMap;
163  ParticleToVertexMap particleToVertexMap;
164  particleToVertexMap[parton1] = vertex1;
165  particleToVertexMap[parton2] = vertex2;
166  for (unsigned int i = 0, n = genParticlesHandle->size(); i < n; ++i) {
167  const reco::Candidate* p = &genParticlesHandle->at(i);
168  if (p == parton1 or p == parton2)
169  continue;
170 
171  // Connect mother-daughters for the other cases
172  for (unsigned int j = 0, nMothers = p->numberOfMothers(); j < nMothers; ++j) {
173  // Mother-daughter hierarchy defines vertex
174  const reco::Candidate* elder = p->mother(j)->daughter(0);
175  HepMC::GenVertex* vertex;
176  if (particleToVertexMap.find(elder) == particleToVertexMap.end()) {
177  vertex = new HepMC::GenVertex(FourVector(elder->vertex()));
178  hepmc_event->add_vertex(vertex);
179  particleToVertexMap[elder] = vertex;
180  } else {
181  vertex = particleToVertexMap[elder];
182  }
183 
184  // Vertex is found. Now connect each other
185  const reco::Candidate* mother = p->mother(j);
186  vertex->add_particle_in(genCandToHepMCMap[mother]);
187  vertex->add_particle_out(hepmc_particles[i]);
188  }
189  }
190 
191  // Finalize HepMC event record
192  bool hasSignalVertex = false;
193  if (!signalParticlePdgIds_.empty()) {
194  // Loop over all vertices to assign the signal vertex, decaying to a signal particle
195  for (auto v = hepmc_event->vertices_begin(); v != hepmc_event->vertices_end(); ++v) {
196  for (auto p = (*v)->particles_begin(HepMC::children); p != (*v)->particles_end(HepMC::children); ++p) {
197  const int pdgId = (*p)->pdg_id();
199  signalParticlePdgIds_.end()) {
200  hepmc_event->set_signal_process_vertex(*v);
201  hasSignalVertex = true;
202  break;
203  }
204  }
205  if (hasSignalVertex)
206  break;
207  }
208  }
209  // Set the default signal vertex if still not set
210  if (!hasSignalVertex)
211  hepmc_event->set_signal_process_vertex(*(vertex1->vertices_begin()));
212 
213  std::unique_ptr<edm::HepMCProduct> hepmc_product(new edm::HepMCProduct());
214  hepmc_product->addHepMCData(hepmc_event);
215  event.put(std::move(hepmc_product), "unsmeared");
216 }
217 
void produce(edm::Event &event, const edm::EventSetup &eventSetup) override
T getParameter(std::string const &) const
EventNumber_t event() const
Definition: EventID.h:40
edm::EDGetTokenT< GenRunInfoProduct > genRunInfoToken_
virtual double pz() const =0
z coordinate of momentum vector
edm::EDGetTokenT< GenEventInfoProduct > genEventInfoToken_
const PDF * pdf() const
double alphaQCD() const
void beginRun(edm::Run const &iRun, edm::EventSetup const &) override
virtual const Candidate * daughter(size_type i) const =0
return daughter at a given position, i = 0, ... numberOfDaughters() - 1 (read only mode) ...
std::pair< double, double > x
Definition: PdfInfo.h:13
edm::ESHandle< ParticleDataTable > pTable_
size_type size() const
virtual const Candidate * mother(size_type i=0) const =0
return pointer to mother
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
HepMC::FourVector FourVector(const reco::Candidate::Point &point)
edm::EDGetTokenT< reco::CandidateView > genParticlesToken_
std::pair< double, double > xPDF
Definition: PdfInfo.h:14
bool getData(T &iHolder) const
Definition: EventSetup.h:113
virtual int status() const =0
status word
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
double qScale() const
virtual const LorentzVector & p4() const =0
four-momentum Lorentz vector
HepMC::FourVector FourVector(const reco::Candidate::LorentzVector &lvec)
virtual int pdgId() const =0
PDG identifier.
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Run.h:315
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e< void, edm::EventID const &, edm::Timestamp const & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
Definition: Activities.doc:12
virtual size_type numberOfMothers() const =0
number of mothers (zero or one in most of but not all the cases)
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double alphaQED() const
std::pair< int, int > id
Definition: PdfInfo.h:12
virtual double eta() const =0
momentum pseudorapidity
virtual double mass() const =0
mass
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Candidate.h:37
const XSec & internalXSec() const
edm::EventID id() const
Definition: EventBase.h:59
std::vector< double > & weights()
const_reference at(size_type pos) const
math::XYZPoint Point
point in the space
Definition: Candidate.h:41
virtual const Point & vertex() const =0
vertex position
unsigned int signalProcessID() const
GenParticles2HepMCConverter(const edm::ParameterSet &pset)
double scalePDF
Definition: PdfInfo.h:15
def move(src, dest)
Definition: eostools.py:511
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
Definition: event.py:1
Definition: Run.h:45