CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
BasicHepMCValidation.cc
Go to the documentation of this file.
1 /*class BasicHepMCValidation
2  *
3  * Class to fill dqm monitor elements from existing EDM file
4  *
5  */
6 
7 #include "Validation/EventGenerator/interface/BasicHepMCValidation.h"
8 
9 #include "CLHEP/Units/defs.h"
10 #include "CLHEP/Units/PhysicalConstants.h"
11 #include "Validation/EventGenerator/interface/DQMHelper.h"
12 using namespace edm;
13 
14 BasicHepMCValidation::BasicHepMCValidation(const edm::ParameterSet& iPSet):
15  wmanager_(iPSet,consumesCollector()),
16  hepmcCollection_(iPSet.getParameter<edm::InputTag>("hepmcCollection"))
17 {
18  hepmcCollectionToken_=consumes<HepMCProduct>(hepmcCollection_);
19 }
20 
21 BasicHepMCValidation::~BasicHepMCValidation() {}
22 
23 void BasicHepMCValidation::dqmBeginRun(const edm::Run& r, const edm::EventSetup& c) {
24  c.getData( fPDGTable );
25 }
26 
27 void BasicHepMCValidation::bookHistograms(DQMStore::IBooker &i, edm::Run const &, edm::EventSetup const &){
29  DQMHelper dqm(&i); i.setCurrentFolder("Generator/Particles");
30 
31  // Number of analyzed events
32  nEvt = dqm.book1dHisto("nEvt", "n analyzed Events", 1, 0., 1.);
33 
36  // quarks
37  particles.push_back(ParticleMonitor("u",1,i));
38  particles.push_back(ParticleMonitor("ubar",-1,i));
39  particles.push_back(ParticleMonitor("d",2,i));
40  particles.push_back(ParticleMonitor("dbar",-2,i));
41  particles.push_back(ParticleMonitor("s",3,i));
42  particles.push_back(ParticleMonitor("sbar",-3,i));
43  particles.push_back(ParticleMonitor("c",4,i));
44  particles.push_back(ParticleMonitor("cbar",-4,i));
45  particles.push_back(ParticleMonitor("b",5,i));
46  particles.push_back(ParticleMonitor("bbar",-5,i));
47  particles.push_back(ParticleMonitor("t",6,i));
48  particles.push_back(ParticleMonitor("tbar",-6,i));
49 
50  //leptons
51  particles.push_back(ParticleMonitor("eminus",11,i));
52  particles.push_back(ParticleMonitor("eplus",-11,i));
53  particles.push_back(ParticleMonitor("nue",12,i));
54  particles.push_back(ParticleMonitor("nuebar",-12,i));
55  particles.push_back(ParticleMonitor("muminus",13,i));
56  particles.push_back(ParticleMonitor("muplus",-13,i));
57  particles.push_back(ParticleMonitor("numu",14,i));
58  particles.push_back(ParticleMonitor("numubar",-14,i));
59  particles.push_back(ParticleMonitor("tauminus",15,i));
60  particles.push_back(ParticleMonitor("tauplus",-15,i));
61  particles.push_back(ParticleMonitor("nutau",16,i));
62  particles.push_back(ParticleMonitor("nutaubar",-16,i));
63 
64  //bosons
65  particles.push_back(ParticleMonitor("Wplus",24,i));
66  particles.push_back(ParticleMonitor("Wminus",-24,i));
67  particles.push_back(ParticleMonitor("Z",23,i));
68  particles.push_back(ParticleMonitor("gamma",22,i));
69  particles.push_back(ParticleMonitor("gluon",21,i));
70 
71  //mesons
72  particles.push_back(ParticleMonitor("piplus",211,i,true)); //log
73  particles.push_back(ParticleMonitor("piminus",-211,i,true)); //log
74  particles.push_back(ParticleMonitor("pizero",111,i,true)); //log
75  particles.push_back(ParticleMonitor("Kplus",321,i));
76  particles.push_back(ParticleMonitor("Kminus",-321,i));
77  particles.push_back(ParticleMonitor("Klzero",130,i));
78  particles.push_back(ParticleMonitor("Kszero",310,i));
79 
80  //baryons
81  particles.push_back(ParticleMonitor("p",2212,i,true)); //log
82  particles.push_back(ParticleMonitor("pbar",-2212,i,true)); //log
83  particles.push_back(ParticleMonitor("n",2112,i,true)); //log
84  particles.push_back(ParticleMonitor("nbar",-2112,i,true)); //log
85  particles.push_back(ParticleMonitor("lambda0",3122,i));
86  particles.push_back(ParticleMonitor("lambda0bar",-3122,i));
87 
88  //D mesons
89  particles.push_back(ParticleMonitor("Dplus",411,i));
90  particles.push_back(ParticleMonitor("Dminus",-411,i));
91  particles.push_back(ParticleMonitor("Dzero",421,i));
92  particles.push_back(ParticleMonitor("Dzerobar",-421,i));
93 
94  //B mesons
95  particles.push_back(ParticleMonitor("Bplus",521,i));
96  particles.push_back(ParticleMonitor("Bminus",-521,i));
97  particles.push_back(ParticleMonitor("Bzero",511,i));
98  particles.push_back(ParticleMonitor("Bzerobar",-511,i));
99  particles.push_back(ParticleMonitor("Bszero",531,i));
100  particles.push_back(ParticleMonitor("Bszerobar",-531,i));
101 
102  //
103  otherPtclNumber = dqm.book1dHisto("otherPtclNumber", "Log10(No. other ptcls)", 60, -1, 5,"log_{10}(No. other ptcls)","Number of Events"); //Log
104  otherPtclMomentum = dqm.book1dHisto("otherPtclMomentum", "Log10(p) other ptcls", 60, -2, 4,"log10(P^{other ptcls}) (log_{10}(GeV))","Number of Events");
105 
107  genPtclNumber = dqm.book1dHisto("genPtclNumber", "Log10(No. all particles)", 60, -1, 5,"log10(No. all particles)","Number of Events"); //Log
108  genVrtxNumber = dqm.book1dHisto("genVrtxNumber", "Log10(No. all vertexs)", 60, -1, 5,"log10(No. all vertexs)","Number of Events"); //Log
109  //
110  stablePtclNumber= dqm.book1dHisto("stablePtclNumber", "Log10(No. stable particles)", 50, 0, 5,"log10(No. stable particles)","Number of Events"); //Log
111  stablePtclPhi = dqm.book1dHisto("stablePtclPhi", "stable Ptcl Phi", 360, -180, 180,"#phi^{stable Ptcl} (rad)","Number of Events");
112  stablePtclEta = dqm.book1dHisto("stablePtclEta", "stable Ptcl Eta (pseudo rapidity)", 220, -11, 11,"#eta^{stable Ptcl}","Number of Events");
113  stablePtclCharge = dqm.book1dHisto("stablePtclCharge", "stablePtclCharge", 5, -2, 2,"Charge^{stable ptcls}","Number of Events");
114  stableChaNumber= dqm.book1dHisto("stableChaNumber", "Log10(No. stable charged particles)", 50, 0, 5,"log_{10}(No. stable charged particles)","Number of Events"); //Log
115  stablePtclp = dqm.book1dHisto("stablePtclp", "Log10(p) stable ptcl p", 80, -4, 4,"log_{10}(P^{stable ptcl}) (log_{10}(GeV))","Number of Events"); //Log
116  stablePtclpT = dqm.book1dHisto("stablePtclpT", "Log10(pT) stable ptcl pT", 80, -4, 4,"log_{10}(P_{t}^{stable ptcl}) (log_{10}(GeV))","Number of Events"); //Log
117  partonNumber = dqm.book1dHisto("partonNumber", "number of partons", 100, 0, 100,"number of partons","Number of Events");
118  partonpT = dqm.book1dHisto("partonpT", "Log10(pT) parton pT", 80, -4, 4,"Log10(P_{t}^{parton})","Number of Events"); //Log
119  outVrtxStablePtclNumber = dqm.book1dHisto("outVrtxStablePtclNumber", "No. outgoing stable ptcls from vrtx", 10, 0, 10,"No. outgoing stable ptcls from vrtx","Number of Events");
120  //
121  outVrtxPtclNumber = dqm.book1dHisto("outVrtxPtclNumber", "No. outgoing ptcls from vrtx", 30, 0, 30,"No. outgoing ptcls from vrtx","Number of Events");
122  vrtxZ = dqm.book1dHisto("VrtxZ", "VrtxZ", 50 , -250, 250,"Z_{Vtx}","Number of Events");
123  vrtxRadius = dqm.book1dHisto("vrtxRadius", "vrtxRadius", 50, 0, 50,"R_{vtx}","Number of Events");
124  //
125  unknownPDTNumber = dqm.book1dHisto("unknownPDTNumber", "Log10(No. unknown ptcls PDT)", 60, -1, 5,"log_{10}(No. unknown ptcls PDT)","Number of Events"); //Log
126  genPtclStatus = dqm.book1dHisto("genPtclStatus", "Status of genParticle", 200,0,200.,"","Number of Events");
127  //
128  Bjorken_x = dqm.book1dHisto("Bjorken_x", "Bjorken_x", 1000, 0.0, 1.0,"Bjorken_{x}","Number of Events");
129  //
130  status1ShortLived = dqm.book1dHisto("status1ShortLived","Status 1 short lived", 11, 0, 11,"","Number of Events");
131  status1ShortLived->setBinLabel(1,"d/dbar");
132  status1ShortLived->setBinLabel(2,"u/ubar");
133  status1ShortLived->setBinLabel(3,"s/sbar");
134  status1ShortLived->setBinLabel(4,"c/cbar");
135  status1ShortLived->setBinLabel(5,"b/bbar");
136  status1ShortLived->setBinLabel(6,"t/tbar");
137  status1ShortLived->setBinLabel(7,"g");
138  status1ShortLived->setBinLabel(8,"tau-/tau+");
139  status1ShortLived->setBinLabel(9,"Z0");
140  status1ShortLived->setBinLabel(10,"W-/W+");
141  status1ShortLived->setBinLabel(11,"PDG = 7,8,17,25-99");
142 
143  log10DeltaEcms =dqm.book1dHisto("DeltaEcms1log10","log_{10} of deviation from nominal Ecms", 200,-1., 5.,"log_{10}(#DeltaE) (log_{10}(GeV))","Number of Events");
144  DeltaEcms = dqm.book1dHisto("DeltaEcms1","deviation from nominal Ecms", 200,-1., 1.,"#DeltaE (GeV)","Number of Events");
145  DeltaPx = dqm.book1dHisto("DeltaPx1","deviation from nominal Px", 200,-1., 1.,"#DeltaP_{x} (GeV)","Number of Events");
146  DeltaPy = dqm.book1dHisto("DeltaPy1","deviation from nominal Py", 200,-1., 1.,"#DeltaP_{y} (GeV)","Number of Events");
147  DeltaPz = dqm.book1dHisto("DeltaPz1","deviation from nominal Pz", 200,-1., 1.,"#DeltaP_{z} (GeV)","Number of Events");
148  return;
149 }
150 
151 void BasicHepMCValidation::analyze(const edm::Event& iEvent,const edm::EventSetup& iSetup)
152 {
154  int partonNum = 0;
155  //
156  int outVrtxStablePtclNum = 0; int stablePtclNum = 0; int otherPtclNum = 0; int unknownPDTNum = 0; int stableChaNum = 0;
157  //
158  double bjorken = 0.;
159  //
160  double etotal = 0. ; double pxtotal = 0.; double pytotal = 0.; double pztotal = 0.;
161 
163  edm::Handle<HepMCProduct> evt;
164  iEvent.getByToken(hepmcCollectionToken_, evt);
165 
166  //Get EVENT
167  HepMC::GenEvent const *myGenEvent = evt->GetEvent();
168 
169  double weight = wmanager_.weight(iEvent);
170 
171  nEvt->Fill(0.5,weight);
172 
173  genPtclNumber->Fill(log10(myGenEvent->particles_size()),weight);
174  genVrtxNumber->Fill(log10(myGenEvent->vertices_size()),weight);
175 
177  HepMC::PdfInfo const *pdf = myGenEvent->pdf_info();
178  if(pdf){
179  bjorken = ((pdf->x1())/((pdf->x1())+(pdf->x2())));
180  }
181  Bjorken_x->Fill(bjorken,weight);
182  //Looping through the VERTICES in the event
183  HepMC::GenEvent::vertex_const_iterator vrtxBegin = myGenEvent->vertices_begin();
184  HepMC::GenEvent::vertex_const_iterator vrtxEnd = myGenEvent->vertices_end();
185  unsigned int nvtx(0);
186  for(HepMC::GenEvent::vertex_const_iterator vrtxIt = vrtxBegin; vrtxIt!=vrtxEnd; ++vrtxIt)
187  {
189  HepMC::GenVertex const *vrtx = *vrtxIt;
190  outVrtxPtclNumber->Fill(vrtx->particles_out_size(),weight); //std::cout << "all " << vrtx->particles_out_size() << '\n';
191 
192  if(nvtx==0){
193  vrtxZ->Fill(vrtx->point3d().z(),weight);
194  vrtxRadius->Fill(vrtx->point3d().perp(),weight);
195  }
197  HepMC::GenVertex::particles_out_const_iterator vrtxPtclBegin = vrtx->particles_out_const_begin();
198  HepMC::GenVertex::particles_out_const_iterator vrtxPtclEnd = vrtx->particles_out_const_end();
199  outVrtxStablePtclNum = 0;
200  for(HepMC::GenVertex::particles_out_const_iterator vrtxPtclIt = vrtxPtclBegin; vrtxPtclIt != vrtxPtclEnd; ++vrtxPtclIt)
201  {
202  HepMC::GenParticle const *vrtxPtcl = *vrtxPtclIt;
203  if (vrtxPtcl->status() == 1){
204  ++outVrtxStablePtclNum; //std::cout << "stable " << outVrtxStablePtclNum << '\n';
205  }
206  }
207  outVrtxStablePtclNumber->Fill(outVrtxStablePtclNum,weight);
208  nvtx++;
209  }//vertices
210 
212  HepMC::GenEvent::particle_const_iterator ptclBegin = myGenEvent->particles_begin();
213  HepMC::GenEvent::particle_const_iterator ptclEnd = myGenEvent->particles_end();
214  for(HepMC::GenEvent::particle_const_iterator ptclIt = ptclBegin; ptclIt!=ptclEnd; ++ptclIt)
215  {
216 
218  HepMC::GenParticle const *ptcl = *ptclIt;
219  int Id = ptcl->pdg_id(); // std::cout << Id << '\n';
220  float Log_p = log10( ptcl->momentum().rho() );
221  double charge = 999.; // for the charge it's needed a HepPDT method
222  int status = ptcl->status();
223  const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID(Id));
224  if(PData==0) {
225  // std::cout << "Unknown id = " << Id << '\n';
226  ++unknownPDTNum;
227  }
228  else
229  charge = PData->charge();
230 
232  genPtclStatus->Fill((float)status,weight);
233 
235  if(ptcl->status() == 1){
236  ++stablePtclNum;
237  stablePtclPhi->Fill(ptcl->momentum().phi()/CLHEP::degree,weight); //std::cout << ptcl->polarization().phi() << '\n';
238  stablePtclEta->Fill(ptcl->momentum().pseudoRapidity(),weight);
239  stablePtclCharge->Fill(charge,weight); // std::cout << ptclData.charge() << '\n';
240  stablePtclp->Fill(Log_p,weight);
241  stablePtclpT->Fill(log10(ptcl->momentum().perp()),weight);
242  if (charge != 0. && charge != 999.) ++stableChaNum;
243  if ( std::abs(Id) == 1 ) status1ShortLived->Fill(1,weight);
244  if ( std::abs(Id) == 2 ) status1ShortLived->Fill(2,weight);
245  if ( std::abs(Id) == 3 ) status1ShortLived->Fill(3,weight);
246  if ( std::abs(Id) == 4 ) status1ShortLived->Fill(4,weight);
247  if ( std::abs(Id) == 5 ) status1ShortLived->Fill(5,weight);
248  if ( std::abs(Id) == 6 ) status1ShortLived->Fill(6,weight);
249  if ( Id == 21 ) status1ShortLived->Fill(7,weight);
250  if ( std::abs(Id) == 15 ) status1ShortLived->Fill(8,weight);
251  if ( Id == 23 ) status1ShortLived->Fill(9,weight);
252  if ( std::abs(Id) == 24 ) status1ShortLived->Fill(10,weight);
253  if ( std::abs(Id) == 7 || std::abs(Id) == 8 || std::abs(Id) == 17 || (std::abs(Id) >= 25 && std::abs(Id) <= 99) ) status1ShortLived->Fill(11,weight);
254  etotal += ptcl->momentum().e();
255  pxtotal += ptcl->momentum().px();
256  pytotal += ptcl->momentum().py();
257  pztotal += ptcl->momentum().pz();
258  }
259 
260  if (abs(Id) < 6 || abs(Id) == 22){
261  ++partonNum; partonpT->Fill(Log_p,weight);
262  }
263 
264  bool indentified=false;
265  for(unsigned int i=0;i<particles.size();i++){if(particles.at(i).Fill(ptcl,weight)){indentified=true; break;}}
266  if(!indentified){ ++otherPtclNum; otherPtclMomentum->Fill(Log_p,weight);}
267  }//event particles
268 
269  // set a default sqrt(s) and then check in the event
270  double ecms = 13000.;
271  if ( myGenEvent->valid_beam_particles() ) {
272  ecms = myGenEvent->beam_particles().first->momentum().e()+myGenEvent->beam_particles().second->momentum().e();
273  }
274  log10DeltaEcms->Fill(log10(fabs(etotal-ecms)),weight);
275  DeltaEcms->Fill(etotal-ecms,weight);
276  DeltaPx->Fill(pxtotal,weight);
277  DeltaPy->Fill(pytotal,weight);
278  DeltaPz->Fill(pztotal,weight);
279 
280 
282  stablePtclNumber->Fill(log10(stablePtclNum+0.1),weight);
283  stableChaNumber->Fill(log10(stableChaNum+0.1),weight);
284  otherPtclNumber->Fill(log10(otherPtclNum+0.1),weight);
285  unknownPDTNumber->Fill(log10(unknownPDTNum+0.1),weight);
286  //
287  partonNumber->Fill(partonNum,weight);
288  for(unsigned int i=0;i<particles.size();i++){particles.at(i).FillCount(weight);};
289 
290 }//analyze
BasicHepMCValidation::genPtclStatus
MonitorElement * genPtclStatus
Definition: BasicHepMCValidation.h:169
BasicHepMCValidation::DeltaEcms
MonitorElement * DeltaEcms
Definition: BasicHepMCValidation.h:190
i
int i
Definition: DBlmapReader.cc:9
BasicHepMCValidation::genVrtxNumber
MonitorElement * genVrtxNumber
Definition: BasicHepMCValidation.h:166
BasicHepMCValidation::unknownPDTNumber
MonitorElement * unknownPDTNumber
Definition: BasicHepMCValidation.h:167
BasicHepMCValidation::outVrtxPtclNumber
MonitorElement * outVrtxPtclNumber
Definition: BasicHepMCValidation.h:168
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:457
BasicHepMCValidation::Bjorken_x
MonitorElement * Bjorken_x
Definition: BasicHepMCValidation.h:185
BasicHepMCValidation::otherPtclNumber
MonitorElement * otherPtclNumber
other ME&#39;s
Definition: BasicHepMCValidation.h:163
BasicHepMCValidation::dqmBeginRun
virtual void dqmBeginRun(const edm::Run &r, const edm::EventSetup &c) override
Definition: BasicHepMCValidation.cc:23
BasicHepMCValidation::BasicHepMCValidation
BasicHepMCValidation(const edm::ParameterSet &)
Definition: BasicHepMCValidation.cc:14
BasicHepMCValidation::otherPtclMomentum
MonitorElement * otherPtclMomentum
Definition: BasicHepMCValidation.h:164
MonitorElement::setBinLabel
void setBinLabel(int bin, const std::string &label, int axis=1)
set bin label for x, y or z axis (axis=1, 2, 3 respectively)
Definition: MonitorElement.cc:920
BasicHepMCValidation::stablePtclCharge
MonitorElement * stablePtclCharge
Definition: BasicHepMCValidation.h:175
BasicHepMCValidation::partonNumber
MonitorElement * partonNumber
Definition: BasicHepMCValidation.h:178
BasicHepMCValidation::hepmcCollectionToken_
edm::EDGetTokenT< edm::HepMCProduct > hepmcCollectionToken_
Definition: BasicHepMCValidation.h:195
DQMHelper::book1dHisto
MonitorElement * book1dHisto(std::string name, std::string title, int n, double xmin, double xmax, std::string xaxis, std::string yaxis)
Definition: DQMHelper.cc:12
edm::EventSetup::getData
void getData(T &iHolder) const
Definition: EventSetup.h:79
MonitorElement::Fill
void Fill(long long x)
Definition: MonitorElement.h:170
iEvent
int iEvent
Definition: GenABIO.cc:230
BasicHepMCValidation::vrtxRadius
MonitorElement * vrtxRadius
Definition: BasicHepMCValidation.h:183
BasicHepMCValidation::fPDGTable
edm::ESHandle< HepPDT::ParticleDataTable > fPDGTable
PDT table.
Definition: BasicHepMCValidation.h:49
BasicHepMCValidation::bookHistograms
virtual void bookHistograms(DQMStore::IBooker &i, edm::Run const &, edm::EventSetup const &) override
Definition: BasicHepMCValidation.cc:27
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
BasicHepMCValidation::log10DeltaEcms
MonitorElement * log10DeltaEcms
Definition: BasicHepMCValidation.h:189
BasicHepMCValidation::stablePtclPhi
MonitorElement * stablePtclPhi
Definition: BasicHepMCValidation.h:173
BasicHepMCValidation::stablePtclEta
MonitorElement * stablePtclEta
Definition: BasicHepMCValidation.h:174
ParticleData
HepPDT::ParticleData ParticleData
Definition: ParticleDataTable.h:9
BasicHepMCValidation::analyze
virtual void analyze(edm::Event const &, edm::EventSetup const &) override
Definition: BasicHepMCValidation.cc:151
BasicHepMCValidation::stablePtclp
MonitorElement * stablePtclp
Definition: BasicHepMCValidation.h:176
BasicHepMCValidation::status1ShortLived
MonitorElement * status1ShortLived
Definition: BasicHepMCValidation.h:187
BasicHepMCValidation::stablePtclpT
MonitorElement * stablePtclpT
Definition: BasicHepMCValidation.h:177
DQMStore::IBooker::setCurrentFolder
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:277
BasicHepMCValidation::stableChaNumber
MonitorElement * stableChaNumber
Definition: BasicHepMCValidation.h:172
BasicHepMCValidation::stablePtclNumber
MonitorElement * stablePtclNumber
Definition: BasicHepMCValidation.h:171
BasicHepMCValidation::genPtclNumber
MonitorElement * genPtclNumber
Definition: BasicHepMCValidation.h:165
BasicHepMCValidation::hepmcCollection_
edm::InputTag hepmcCollection_
Definition: BasicHepMCValidation.h:46
BasicHepMCValidation::outVrtxStablePtclNumber
MonitorElement * outVrtxStablePtclNumber
Definition: BasicHepMCValidation.h:180
histoStyle.weight
int weight
Definition: histoStyle.py:50
WeightManager::weight
double weight(const edm::Event &)
Definition: WeightManager.cc:25
BasicHepMCValidation::particles
std::vector< ParticleMonitor > particles
Definition: BasicHepMCValidation.h:160
edm::Run
Definition: Run.h:42
mps_update.status
tuple status
Definition: mps_update.py:57
Generated for CMSSW Reference Manual by doxygen 1.8.5