11 #include "CLHEP/Units/defs.h"
12 #include "CLHEP/Units/PhysicalConstants.h"
17 hepmcCollection_(iPSet.getParameter<edm::
InputTag>(
"hepmcCollection")),
18 genparticleCollection_(iPSet.getParameter<edm::
InputTag>(
"genparticleCollection")),
19 genjetCollection_(iPSet.getParameter<edm::
InputTag>(
"genjetsCollection")),
20 matchPr_(iPSet.getParameter<double>(
"matchingPrecision")),
21 verbosity_(iPSet.getUntrackedParameter<unsigned int>(
"verbosity",0))
79 unsigned int initSize = 1000;
86 HepMC::GenEvent *myGenEvent =
new HepMC::GenEvent(*(evt->GetEvent()));
90 std::vector<const HepMC::GenParticle*> hepmcGPCollection;
91 std::vector<int> barcodeList;
92 hepmcGPCollection.reserve(initSize);
93 barcodeList.reserve(initSize);
96 for (HepMC::GenEvent::particle_const_iterator iter = myGenEvent->particles_begin(); iter != myGenEvent->particles_end(); ++iter){
97 if ( (*iter)->status() == 1) {
98 hepmcGPCollection.push_back(*iter);
99 barcodeList.push_back((*iter)->barcode());
101 std::cout <<
"HepMC " << std::setw(14) << std::fixed << (*iter)->pdg_id() << std::setw(14) << std::fixed << (*iter)->momentum().px() << std::setw(14) << std::fixed
102 << (*iter)->momentum().py() << std::setw(14) << std::fixed << (*iter)->momentum().pz() << std::endl;
112 std::vector<const reco::GenParticle*> particles;
113 particles.reserve(initSize);
114 for (reco::GenParticleCollection::const_iterator iter=genParticles->begin();iter!=genParticles->end();++iter){
115 if ( (*iter).status() == 1) {
116 particles.push_back(&*iter);
118 std::cout <<
"reco " << std::setw(14) << std::fixed << (*iter).pdgId() << std::setw(14) << std::fixed << (*iter).px()
119 << std::setw(14) << std::fixed << (*iter).py() << std::setw(14) << std::fixed << (*iter).pz() << std::endl;
124 unsigned int nReco = particles.size();
125 unsigned int nHepMC = hepmcGPCollection.size();
130 std::vector<int> hepmcMatchIndex;
131 hepmcMatchIndex.reserve(initSize);
137 if ( nReco != nHepMC ) {
138 edm::LogWarning(
"CollectionSizeInconsistency") <<
"Collection size inconsistency: HepMC::GenParticle = " << nHepMC <<
" reco::GenParticle = " << nReco;
143 for (
unsigned int i = 0;
i < nReco; ++
i ){
144 for (
unsigned int j = 0;
j < nHepMC; ++
j ){
146 hepmcMatchIndex.push_back((
int)j);
147 if ( hepmcGPCollection[j]->momentum().
rho() != 0. ) {
148 double reso = 1.-particles[
i]->p()/hepmcGPCollection[
j]->momentum().rho();
150 std::cout <<
"Matching momentum: reco = " << particles[
i]->p() <<
" HepMC = "
151 << hepmcGPCollection[
j]->momentum().rho() <<
" resoultion = " << reso << std::endl;
161 unsigned int nMatched = hepmcMatchIndex.size();
163 if ( nMatched != nReco ) {
164 edm::LogWarning(
"IncorrectMatching") <<
"Incorrect number of matched indexes: GenParticle = " << nReco <<
" matched indexes = " << nMatched;
168 unsigned int nWrMatch = 0;
170 for (
unsigned int i = 0;
i < nMatched; ++
i ){
171 for (
unsigned int j =
i+1;
j < nMatched; ++
j ){
172 if ( hepmcMatchIndex[
i] == hepmcMatchIndex[
j] ) {
173 int theIndex = hepmcMatchIndex[
i];
174 edm::LogWarning(
"DuplicatedMatching") <<
"Multiple matching occurencies for GenParticle barcode = " << barcodeList[theIndex];
189 int nJetsCentral = 0;
192 std::vector<double>
jetEta;
193 jetEta.reserve(initSize);
195 for (reco::GenJetCollection::const_iterator iter=genJets->begin();iter!=genJets->end();++iter){
197 double pt = (*iter).pt();
199 if (pt > 1.) nJetso1++;
200 if (pt > 10.) nJetso10++;
201 if (pt > 100.) nJetso100++;
202 double eta = (*iter).eta();
203 if ( std::fabs(eta) < 2.5 ) nJetsCentral++;
204 jetEta.push_back(eta);
220 double deltaEta = 999.;
221 if ( jetEta.size() > 1 ) {
222 for (
unsigned int i = 0;
i < jetEta.size();
i++){
223 for (
unsigned int j =
i+1;
j < jetEta.size();
j++){
224 deltaEta =
std::min(deltaEta,std::fabs(jetEta[
i]-jetEta[
j]));
238 if ( hepmcP->pdg_id() != recoP->
pdgId() )
return state;
239 if ( std::fabs(hepmcP->momentum().px()-recoP->
px()) < std::fabs(
matchPr_*hepmcP->momentum().px()) &&
240 std::fabs(hepmcP->momentum().py()-recoP->
py()) < std::fabs(
matchPr_*hepmcP->momentum().py()) &&
241 std::fabs(hepmcP->momentum().pz()-recoP->
pz()) < std::fabs(
matchPr_*hepmcP->momentum().pz())) {
virtual void beginRun(const edm::Run &, const edm::EventSetup &)
MonitorElement * genJetCentral
virtual int pdgId() const
PDG identifier.
edm::ESHandle< HepPDT::ParticleDataTable > fPDGTable
PDT table.
MonitorElement * genMatched
DQMStore * dbe
ME's "container".
MonitorElement * multipleMatching
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
virtual void analyze(const edm::Event &, const edm::EventSetup &)
MonitorElement * genJetEta
MonitorElement * genJetEnergy
MonitorElement * genJetTotPt
edm::InputTag genparticleCollection_
bool matchParticles(const HepMC::GenParticle *&, const reco::GenParticle *&)
virtual ~BasicGenParticleValidation()
void getData(T &iHolder) const
MonitorElement * genJetMult
edm::InputTag hepmcCollection_
MonitorElement * genJetDeltaEtaMin
MonitorElement * genJetPto1
MonitorElement * genJetPt
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
edm::InputTag genjetCollection_
virtual double px() const
x coordinate of momentum vector
MonitorElement * matchedResolution
MonitorElement * genJetPto10
virtual double pz() const
z coordinate of momentum vector
MonitorElement * genPMultiplicity
MonitorElement * genJetPhi
MonitorElement * genJetPto100
BasicGenParticleValidation(const edm::ParameterSet &)
virtual void endRun(const edm::Run &, const edm::EventSetup &)
virtual double py() const
y coordinate of momentum vector
void setCurrentFolder(const std::string &fullpath)