|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
|
#include "GeneratorInterface/GenFilters/plugins/MCSingleParticleFilter.h"
#include "GeneratorInterface/GenFilters/plugins/MCFilterZboostHelper.h"
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include <iostream>
using namespace edm;
using namespace std;
MCSingleParticleFilter::MCSingleParticleFilter(const edm::ParameterSet& iConfig)
: token_(consumes<edm::HepMCProduct>(
edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))),
betaBoost(iConfig.getUntrackedParameter("BetaBoost", 0.)) {
//here do whatever other initialization is needed
vector<int> defpid;
defpid.push_back(0);
particleID = iConfig.getUntrackedParameter<vector<int> >("ParticleID", defpid);
vector<double> defptmin;
defptmin.push_back(0.);
ptMin = iConfig.getUntrackedParameter<vector<double> >("MinPt", defptmin);
vector<double> defetamin;
defetamin.push_back(-10.);
etaMin = iConfig.getUntrackedParameter<vector<double> >("MinEta", defetamin);
vector<double> defetamax;
defetamax.push_back(10.);
etaMax = iConfig.getUntrackedParameter<vector<double> >("MaxEta", defetamax);
vector<int> defstat;
defstat.push_back(0);
status = iConfig.getUntrackedParameter<vector<int> >("Status", defstat);
// check for same size
if ((ptMin.size() > 1 && particleID.size() != ptMin.size()) ||
(etaMin.size() > 1 && particleID.size() != etaMin.size()) ||
(etaMax.size() > 1 && particleID.size() != etaMax.size()) ||
(status.size() > 1 && particleID.size() != status.size())) {
edm::LogInfo("MCSingleParticleFilter")
<< "WARNING: size of MinPthat and/or MaxPthat not matching with ProcessID size!!" << endl;
}
// if ptMin size smaller than particleID , fill up further with defaults
if (particleID.size() > ptMin.size()) {
vector<double> defptmin2;
for (unsigned int i = 0; i < particleID.size(); i++) {
defptmin2.push_back(0.);
}
ptMin = defptmin2;
}
// if etaMin size smaller than particleID , fill up further with defaults
if (particleID.size() > etaMin.size()) {
vector<double> defetamin2;
for (unsigned int i = 0; i < particleID.size(); i++) {
defetamin2.push_back(-10.);
}
etaMin = defetamin2;
}
// if etaMax size smaller than particleID , fill up further with defaults
if (particleID.size() > etaMax.size()) {
vector<double> defetamax2;
for (unsigned int i = 0; i < particleID.size(); i++) {
defetamax2.push_back(10.);
}
etaMax = defetamax2;
}
// if status size smaller than particleID , fill up further with defaults
if (particleID.size() > status.size()) {
vector<int> defstat2;
for (unsigned int i = 0; i < particleID.size(); i++) {
defstat2.push_back(0);
}
status = defstat2;
}
// check if beta is smaller than 1
if (std::abs(betaBoost) >= 1) {
edm::LogError("MCSingleParticleFilter") << "Input beta boost is >= 1 !";
}
}
MCSingleParticleFilter::~MCSingleParticleFilter() {
// do anything here that needs to be done at desctruction time
// (e.g. close files, deallocate resources etc.)
}
// ------------ method called to skim the data ------------
bool MCSingleParticleFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
using namespace edm;
bool accepted = false;
Handle<HepMCProduct> evt;
iEvent.getByToken(token_, evt);
const HepMC::GenEvent* myGenEvent = evt->GetEvent();
for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
++p) {
for (unsigned int i = 0; i < particleID.size(); i++) {
if (particleID[i] == (*p)->pdg_id() || particleID[i] == 0) {
if ((*p)->momentum().perp() > ptMin[i] && ((*p)->status() == status[i] || status[i] == 0)) {
HepMC::FourVector mom = MCFilterZboostHelper::zboost((*p)->momentum(), betaBoost);
if (mom.eta() > etaMin[i] && mom.eta() < etaMax[i]) {
accepted = true;
}
}
}
}
}
return accepted;
}
|