CMS 3D CMS Logo

EcalGenEvtSelectorFrag.cc
Go to the documentation of this file.
1 #include <iostream>
4 #include "HepMC/GenVertex.h"
5 using namespace std;
6 
8 
9  partonId_ = pset.getParameter<vector<int> >("partons");
10  partonStatus_ = pset.getParameter<vector<int> >("partonStatus");
11  partonPt_ = pset.getParameter<vector<double> >("partonPt");
12 
13  particleId_ = pset.getParameter<vector<int> >("particles");
14  particleStatus_ = pset.getParameter<vector<int> >("particleStatus");
15  particlePt_ = pset.getParameter<vector<double> >("particlePt");
16 
17  etaMax_ = pset.getParameter<double>("etaMax");
18 
19  int id = partonId_.size();
20  int st = partonStatus_.size();
21  int pt = partonPt_.size();
22 
23  if(partonId_.size() != partonStatus_.size() || partonId_.size() != partonPt_.size()){
24  throw edm::Exception(edm::errors::LogicError)<<id<<st<<pt<<endl;
25  }
26 
27  id = particleId_.size();
28  st = particleStatus_.size();
29  pt = particlePt_.size();
30 
31  if(particleId_.size() != particleStatus_.size() || particleId_.size() != particlePt_.size()){
32  throw edm::Exception(edm::errors::LogicError)<<id<<st<<pt<<endl;
33 
34  }
35 
36 }
37 
39 
40  HepMC::GenEvent::particle_const_iterator begin = evt->particles_begin();
41  HepMC::GenEvent::particle_const_iterator end = evt->particles_end();
42 
43  bool foundParticle = false;
44  bool foundParton = false;
45 
46  HepMC::GenEvent::particle_const_iterator it = begin;
47  while((!foundParton || !foundParticle) && it != end){
48  bool isFrag = false;
49  bool isThisPhoton = false;
50 
51  foundParton = true;
52  /*for(unsigned i = 0; i < partonId_.size(); ++i){
53  if(selectParticle(*it, partonStatus_[i], partonId_[i], partonPt_[i], etaMax_)) foundParton = true;
54  }*/
55 
56  for(unsigned i = 0; i < particleId_.size(); ++i){
57  if(selectParticle(*it, particleStatus_[i], particleId_[i], particlePt_[i], etaMax_)) isThisPhoton =true;
58  }
59 
60  // Now you have to requre the partcile is "prompt", meaning its mom is parton
61 
62  if ( !((*it)->production_vertex()) ) {
63  isFrag = false;
64  }
65  else {
66  const HepMC::GenVertex* productionVertex = (*it)->production_vertex();
67 
68  size_t numberOfMothers = productionVertex->particles_in_size();
69  if ( numberOfMothers <= 0 ) {
70  isFrag = false ;
71  // cout << "number of mothers = " << numberOfMothers << endl;
72  }
73  else {
74  // cout << "number of mothers = " << numberOfMothers << endl;
75  HepMC::GenVertex::particles_in_const_iterator motherIt = productionVertex->particles_in_const_begin();
76  for( ; motherIt != productionVertex->particles_in_const_end(); motherIt++) {
77  if ( fabs( (*motherIt)->pdg_id() ) <= 22 ) {
78  isFrag = true;
79  }
80  }
81  }
82  }
83 
84  if ( (isFrag == true) && ( isThisPhoton == true) ) {
85  //cout << "Found fragmentation photon!!" << endl ;
86  foundParticle = true;
87  }
88 
89  ++it;
90 
91 
92 
93 
94  }
95 
96  foundParton = true ; // We don't care of the parton
97  return (foundParton && foundParticle);
98 }
T getParameter(std::string const &) const
EcalGenEvtSelectorFrag(const edm::ParameterSet &pset)
bool selectParticle(HepMC::GenParticle *par, int status, int pdg, double ptMin, double etaMax)
std::vector< int > partonStatus_
bool filter(HepMC::GenEvent *) override
std::vector< double > particlePt_
std::vector< int > partonId_
#define end
Definition: vmac.h:39
std::vector< int > particleId_
std::vector< double > partonPt_
#define begin
Definition: vmac.h:32
std::vector< int > particleStatus_