CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
BasicGenParticleValidation.cc
Go to the documentation of this file.
1 /*class BasicGenParticleValidation
2  *
3  * Class to fill dqm monitor elements from existing EDM file
4  *
5  * $Date: 2010/05/25 16:50:50 $
6  * $Revision: 1.1 $
7  */
8 
10 
11 #include "CLHEP/Units/defs.h"
12 #include "CLHEP/Units/PhysicalConstants.h"
13 
14 using namespace edm;
15 
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))
22 {
23  dbe = 0;
25 }
26 
28 
30 {
31  if(dbe){
33  dbe->setCurrentFolder("Generator/GenParticles");
34 
36 
37  // Number of analyzed events
38  nEvt = dbe->book1D("nEvt", "n analyzed Events", 1, 0., 1.);
39 
41  genPMultiplicity = dbe->book1D("genPMultiplicty", "Log(No. all GenParticles)", 50, -1, 5); //Log
42  //difference in HepMC and reco multiplicity
43  genMatched = dbe->book1D("genMatched", "Difference reco - matched", 50, -25, 25);
44  //multiple matching
45  multipleMatching = dbe->book1D("multipleMatching", "multiple reco HepMC matching", 50, 0, 50);
46  //momentum difference of matched particles
47  matchedResolution = dbe->book1D("matchedResolution", "log10(momentum difference of matched particles)", 70, -10., -3.);
48 
49  // GenJet general distributions
50  genJetMult = dbe->book1D("genJetMult", "GenJet multiplicity", 50, 0, 50);
51  genJetEnergy = dbe->book1D("genJetEnergy", "Log10(GenJet energy)", 60, -1, 5);
52  genJetPt = dbe->book1D("genJetPt", "Log10(GenJet pt)", 60, -1, 5);
53  genJetEta = dbe->book1D("genJetEta", "GenJet eta", 220, -11, 11);
54  genJetPhi = dbe->book1D("genJetPhi", "GenJet phi", 360, -180, 180);
55  genJetDeltaEtaMin = dbe->book1D("genJetDeltaEtaMin", "GenJet minimum rapidity gap", 30, 0, 30);
56 
57  genJetPto1 = dbe->book1D("genJetPto1", "GenJet multiplicity above 1 GeV", 50, 0, 50);
58  genJetPto10 = dbe->book1D("genJetPto10", "GenJet multiplicity above 10 GeV", 50, 0, 50);
59  genJetPto100 = dbe->book1D("genJetPto100", "GenJet multiplicity above 100 GeV", 50, 0, 50);
60  genJetCentral = dbe->book1D("genJetCentral", "GenJet multiplicity |eta|.lt.2.5", 50, 0, 50);
61 
62  genJetTotPt = dbe->book1D("genJetTotPt", "Log10(GenJet total pt)", 100, -5, 5);
63 
64  }
65  return;
66 }
67 
70 {
72  iSetup.getData( fPDGTable );
73  return;
74 }
75 void BasicGenParticleValidation::endRun(const edm::Run& iRun,const edm::EventSetup& iSetup){return;}
77 {
78 
79  unsigned int initSize = 1000;
80 
83  iEvent.getByLabel(hepmcCollection_, evt);
84 
85  //Get HepMC EVENT
86  HepMC::GenEvent *myGenEvent = new HepMC::GenEvent(*(evt->GetEvent()));
87 
88  nEvt->Fill(0.5);
89 
90  std::vector<const HepMC::GenParticle*> hepmcGPCollection;
91  std::vector<int> barcodeList;
92  hepmcGPCollection.reserve(initSize);
93  barcodeList.reserve(initSize);
94 
95  //Looping through HepMC::GenParticle collection to search for status 1 particles
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());
100  if ( verbosity_ > 0 ) {
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;
103  }
104  }
105  }
106 
107 
108  // Gather information on the reco::GenParticle collection
110  iEvent.getByLabel(genparticleCollection_, genParticles );
111 
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);
117  if ( verbosity_ > 0 ) {
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;
120  }
121  }
122  }
123 
124  unsigned int nReco = particles.size();
125  unsigned int nHepMC = hepmcGPCollection.size();
126 
127  genPMultiplicity->Fill(std::log10(nReco));
128 
129  // Define vector containing index of hepmc corresponding to the reco::GenParticle
130  std::vector<int> hepmcMatchIndex;
131  hepmcMatchIndex.reserve(initSize);
132 
133  // Matching procedure
134 
135  // Check array size consistency
136 
137  if ( nReco != nHepMC ) {
138  edm::LogWarning("CollectionSizeInconsistency") << "Collection size inconsistency: HepMC::GenParticle = " << nHepMC << " reco::GenParticle = " << nReco;
139  }
140 
141  // Match each HepMC with a reco
142 
143  for ( unsigned int i = 0; i < nReco; ++i ){
144  for ( unsigned int j = 0; j < nHepMC; ++j ){
145  if ( matchParticles( hepmcGPCollection[j], particles[i] ) ) {
146  hepmcMatchIndex.push_back((int)j);
147  if ( hepmcGPCollection[j]->momentum().rho() != 0. ) {
148  double reso = 1.-particles[i]->p()/hepmcGPCollection[j]->momentum().rho();
149  if ( verbosity_ > 0 ) {
150  std::cout << "Matching momentum: reco = " << particles[i]->p() << " HepMC = "
151  << hepmcGPCollection[j]->momentum().rho() << " resoultion = " << reso << std::endl;
152  }
153  matchedResolution->Fill(std::log10(std::fabs(reso))); }
154  continue;
155  }
156  }
157  }
158 
159  // Check unicity of matching
160 
161  unsigned int nMatched = hepmcMatchIndex.size();
162 
163  if ( nMatched != nReco ) {
164  edm::LogWarning("IncorrectMatching") << "Incorrect number of matched indexes: GenParticle = " << nReco << " matched indexes = " << nMatched;
165  }
166  genMatched->Fill(int(nReco-nMatched));
167 
168  unsigned int nWrMatch = 0;
169 
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];
175  nWrMatch++;
176  }
177  }
178  }
179  multipleMatching->Fill(int(nWrMatch));
180 
181  // Gather information in the GenJet collection
183  iEvent.getByLabel(genjetCollection_, genJets );
184 
185  int nJets = 0;
186  int nJetso1 = 0;
187  int nJetso10 = 0;
188  int nJetso100 = 0;
189  int nJetsCentral = 0;
190  double totPt = 0.;
191 
192  std::vector<double> jetEta;
193  jetEta.reserve(initSize);
194 
195  for (reco::GenJetCollection::const_iterator iter=genJets->begin();iter!=genJets->end();++iter){
196  nJets++;
197  double pt = (*iter).pt();
198  totPt += 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);
205 
206  genJetEnergy->Fill(std::log10((*iter).energy()));
207  genJetPt->Fill(std::log10(pt));
208  genJetEta->Fill(eta);
209  genJetPhi->Fill((*iter).phi()/CLHEP::degree);
210  }
211 
212  genJetMult->Fill(nJets);
213  genJetPto1->Fill(nJetso1);
214  genJetPto10->Fill(nJetso10);
215  genJetPto100->Fill(nJetso100);
216  genJetCentral->Fill(nJetsCentral);
217 
218  genJetTotPt->Fill(std::log10(totPt));
219 
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]));
225  }
226  }
227  }
228 
229  genJetDeltaEtaMin->Fill(deltaEta);
230 
231  delete myGenEvent;
232 }//analyze
233 
235 
236  bool state = false;
237 
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())) {
242  state = true; }
243 
244  return state;
245 
246 }
virtual void beginRun(const edm::Run &, const edm::EventSetup &)
int i
Definition: DBlmapReader.cc:9
virtual int pdgId() const
PDG identifier.
edm::ESHandle< HepPDT::ParticleDataTable > fPDGTable
PDT table.
DQMStore * dbe
ME&#39;s &quot;container&quot;.
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
Definition: DQMStore.cc:514
virtual void analyze(const edm::Event &, const edm::EventSetup &)
Definition: DDAxes.h:10
#define min(a, b)
Definition: mlp_lapack.h:161
bool matchParticles(const HepMC::GenParticle *&, const reco::GenParticle *&)
T eta() const
void getData(T &iHolder) const
Definition: EventSetup.h:67
void Fill(long long x)
int iEvent
Definition: GenABIO.cc:243
int j
Definition: DBlmapReader.cc:9
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:355
virtual double px() const
x coordinate of momentum vector
virtual double pz() const
z coordinate of momentum vector
char state
Definition: procUtils.cc:75
tuple cout
Definition: gather_cfg.py:41
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)
Definition: DQMStore.cc:232
Definition: Run.h:32