CMS 3D CMS Logo

MCSingleParticleYPt.cc
Go to the documentation of this file.
1 // Filter based on MCSingleParticleFilter.cc, but using rapidity instead of eta
2 
5 #include <iostream>
6 
7 using namespace edm;
8 using namespace std;
9 
10 
12 token_(consumes<edm::HepMCProduct>(edm::InputTag(iConfig.getUntrackedParameter("moduleLabel",std::string("generator")),"unsmeared")))
13 {
14  fVerbose = iConfig.getUntrackedParameter("verbose",0);
15  fchekantiparticle = iConfig.getUntrackedParameter("CheckAntiparticle",true);
16  //here do whatever other initialization is needed
17  vector<int> defpid ;
18  defpid.push_back(0) ;
19  particleID = iConfig.getUntrackedParameter< vector<int> >("ParticleID",defpid);
20  vector<double> defptmin ;
21  defptmin.push_back(0.);
22  ptMin = iConfig.getUntrackedParameter< vector<double> >("MinPt", defptmin);
23  vector<double> defrapmin ;
24  defrapmin.push_back(-10.);
25  rapMin = iConfig.getUntrackedParameter< vector<double> >("MinY", defrapmin);
26  vector<double> defrapmax ;
27  defrapmax.push_back(10.);
28  rapMax = iConfig.getUntrackedParameter< vector<double> >("MaxY", defrapmax);
29  vector<int> defstat ;
30  defstat.push_back(0);
31  status = iConfig.getUntrackedParameter< vector<int> >("Status", defstat);
32 
33  // check for same size
34  if ( (ptMin.size() > 1 && particleID.size() != ptMin.size())
35  || (rapMin.size() > 1 && particleID.size() != rapMin.size())
36  || (rapMax.size() > 1 && particleID.size() != rapMax.size())
37  || (status.size() > 1 && particleID.size() != status.size()) ) {
38  edm::LogWarning("MCSingleParticleYPt") << "WARNING: MCSingleParticleYPt : size of vector cuts do not match!!" << endl;
39  }
40 
41  // if ptMin size smaller than particleID , fill up further with defaults
42  if (particleID.size() > ptMin.size() ){
43  vector<double> defptmin2 ;
44  for (unsigned int i = 0; i < particleID.size(); i++){ defptmin2.push_back(0.);}
45  ptMin = defptmin2;
46  }
47  // if etaMin size smaller than particleID , fill up further with defaults
48  if (particleID.size() > rapMin.size() ){
49  vector<double> defrapmin2 ;
50  for (unsigned int i = 0; i < particleID.size(); i++){ defrapmin2.push_back(-10.);}
51  rapMin = defrapmin2;
52  }
53  // if etaMax size smaller than particleID , fill up further with defaults
54  if (particleID.size() > rapMax.size() ){
55  vector<double> defrapmax2 ;
56  for (unsigned int i = 0; i < particleID.size(); i++){ defrapmax2.push_back(10.);}
57  rapMax = defrapmax2;
58  }
59  // if status size smaller than particleID , fill up further with defaults
60  if (particleID.size() > status.size() ){
61  vector<int> defstat2 ;
62  for (unsigned int i = 0; i < particleID.size(); i++){ defstat2.push_back(0);}
63  status = defstat2;
64  }
65 
66  if (fVerbose > 0) {
67  edm::LogInfo("MCSingleParticleYPt") << "----------------------------------------------------------------------" << std::endl;
68  edm::LogInfo("MCSingleParticleYPt") << "----- MCSingleParticleYPt" << std::endl;
69  for (unsigned int i=0; i<particleID.size(); ++i) {
70  edm::LogInfo("MCSingleParticleYPt") << " ID: " << particleID[i] << " pT > " << ptMin[i] << ", " << rapMin[i] << " < y < " << rapMax[i] << ", status = " << status[i] << std::endl;
71  }
72  if (fchekantiparticle) edm::LogInfo("MCSingleParticleYPt") << " anti-particles will be tested as well." << std::endl;
73  edm::LogInfo("MCSingleParticleYPt") << "----------------------------------------------------------------------" << std::endl;
74  }
75 }
76 
77 
79 {
80  // do anything here that needs to be done at desctruction time
81  // (e.g. close files, deallocate resources etc.)
82 }
83 
84 
85 // ------------ method called to skim the data ------------
87 {
88  using namespace edm;
89  bool accepted = false;
91  iEvent.getByToken(token_, evt);
92 
93  const HepMC::GenEvent * myGenEvent = evt->GetEvent();
94  for ( HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin();
95  p != myGenEvent->particles_end(); ++p ) {
96  if (fVerbose > 3) edm::LogInfo("MCSingleParticleYPt") << "Looking at particle : " << (*p)->pdg_id() << " status : " << (*p)->status() << std::endl;
97 
98  for (unsigned int i = 0; i < particleID.size(); i++) {
99  if (particleID[i] == (*p)->pdg_id() || (fchekantiparticle && (-particleID[i] == (*p)->pdg_id())) || particleID[i] == 0) {
100  // calculate rapidity just for the desired particle and make sure, this particles has enough energy
101  rapidity = ((*p)->momentum().e()-(*p)->momentum().pz()) > 0. ? 0.5*log( ((*p)->momentum().e()+(*p)->momentum().pz()) / ((*p)->momentum().e()-(*p)->momentum().pz()) ) : rapMax[i]+.1;
102  if (fVerbose > 2) edm::LogInfo("MCSingleParticleYPt") << "Testing particle : " << (*p)->pdg_id() << " pT: " << (*p)->momentum().perp() << " y: " << rapidity << " status : " << (*p)->status() << endl;
103  if ( (*p)->momentum().perp() > ptMin[i]
104  && rapidity > rapMin[i]
105  && rapidity < rapMax[i]
106  && ((*p)->status() == status[i] || status[i] == 0) ) {
107  accepted = true;
108  if (fVerbose > 1)
109  edm::LogInfo("MCSingleParticleYPt") << "Accepted particle : " << (*p)->pdg_id() << " pT: " << (*p)->momentum().perp() << " y: " << rapidity << " status : " << (*p)->status() << endl;
110  break;
111  }
112 
113  }
114  }
115  if (accepted) break;
116  }
117 
118  if (accepted) { return true; }
119  else { return false; }
120 }
121 
T getUntrackedParameter(std::string const &, T const &) const
std::vector< int > particleID
std::vector< int > status
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:460
std::vector< double > rapMax
virtual bool filter(edm::Event &, const edm::EventSetup &)
int iEvent
Definition: GenABIO.cc:230
std::vector< double > ptMin
const HepMC::GenEvent * GetEvent() const
Definition: HepMCProduct.h:38
HLT enums.
MCSingleParticleYPt(const edm::ParameterSet &)
edm::EDGetTokenT< edm::HepMCProduct > token_
std::vector< double > rapMin