14 label_(iConfig.getUntrackedParameter(
"moduleLabel",std::string(
"generator"))),
15 particleCharge(iConfig.getUntrackedParameter(
"ParticleCharge",0)),
16 minInvMass(iConfig.getUntrackedParameter(
"MinInvMass", 0.)),
17 maxInvMass(iConfig.getUntrackedParameter(
"MaxInvMass", 14000.)),
18 minDeltaPhi(iConfig.getUntrackedParameter(
"MinDeltaPhi", 0.)),
19 maxDeltaPhi(iConfig.getUntrackedParameter(
"MaxDeltaPhi", 6.3)),
20 minDeltaR(iConfig.getUntrackedParameter(
"MinDeltaR",0.)),
21 maxDeltaR(iConfig.getUntrackedParameter(
"MaxDeltaR",10000.))
30 vector<double> defptmin;
31 defptmin.push_back(0.);
33 vector<double> defpmin;
34 defpmin.push_back(0.);
37 vector<double> defetamin;
38 defetamin.push_back(-10.);
40 vector<double> defetamax ;
41 defetamax.push_back(10.);
54 cout <<
"WARNING: MCParticlePairFilter : size of some vectors not matching with 2!!" << endl;
58 if (2 >
ptMin.size() ){
59 vector<double> defptmin2 ;
60 for (
unsigned int i = 0;
i < 2;
i++){ defptmin2.push_back(0.);}
64 if (2 >
pMin.size() ){
65 vector<double> defpmin2 ;
66 for (
unsigned int i = 0;
i < 2;
i++){ defpmin2.push_back(0.);}
71 vector<double> defetamin2 ;
72 for (
unsigned int i = 0;
i < 2;
i++){ defetamin2.push_back(-10.);}
77 vector<double> defetamax2 ;
78 for (
unsigned int i = 0;
i < 2;
i++){ defetamax2.push_back(10.);}
83 vector<int> defstat2 ;
84 for (
unsigned int i = 0;
i < 2;
i++){ defstat2.push_back(0);}
104 bool accepted =
false;
107 const double pi = 3.14159;
109 vector<HepMC::GenParticle*> typeApassed;
110 vector<HepMC::GenParticle*> typeBpassed;
113 const HepMC::GenEvent * myGenEvent = evt->GetEvent();
116 for ( HepMC::GenEvent::particle_const_iterator
p = myGenEvent->particles_begin();
117 p != myGenEvent->particles_end(); ++
p ) {
120 bool gottypeAID =
false;
128 if ( (*p)->momentum().perp() >
ptMin[0] && (*p)->momentum().rho() >
pMin[0] && (*p)->momentum().eta() >
etaMin[0]
129 && (*p)->momentum().eta() <
etaMax[0] && ((*p)->status() ==
status[0] ||
status[0] == 0)) {
134 double phi1 = (*p)->momentum().phi();
137 double eta1 = (*p)->momentum().eta();
149 while(!accepted && i<typeBpassed.size()) {
150 tot_x=(*p)->momentum().px();
151 tot_y=(*p)->momentum().py();
152 tot_z=(*p)->momentum().pz();
153 tot_e=(*p)->momentum().e();
154 charge1 =
charge((*p)->pdg_id());
157 tot_x += typeBpassed[
i]->momentum().px();
158 tot_y += typeBpassed[
i]->momentum().py();
159 tot_z += typeBpassed[
i]->momentum().pz();
160 tot_e += typeBpassed[
i]->momentum().e();
161 invmass=
sqrt(tot_e*tot_e-tot_x*tot_x-tot_y*tot_y-tot_z*tot_z);
162 combcharge = charge1 *
charge(typeBpassed[i]->pdg_id());
164 phi2 = typeBpassed[
i]->momentum().phi();
165 deltaphi = fabs(phi1-phi2);
166 if(deltaphi >
pi) deltaphi = 2.*
pi-deltaphi;
168 eta2 = typeBpassed[
i]->momentum().eta();
169 deltaeta=fabs(eta1-eta2);
170 deltaR =
sqrt(deltaeta*deltaeta+deltaphi*deltaphi);
183 typeApassed.push_back(*
p);
190 bool gottypeBID =
false;
198 if ( (*p)->momentum().perp() >
ptMin[1] && (*p)->momentum().rho() >
pMin[1] && (*p)->momentum().eta() >
etaMin[1]
199 && (*p)->momentum().eta() <
etaMax[1] && ((*p)->status() ==
status[1] ||
status[1] == 0)) {
204 double phi1 = (*p)->momentum().phi();
207 double eta1 = (*p)->momentum().eta();
219 while(!accepted && i<typeApassed.size()) {
220 if((*
p) != typeApassed[
i]) {
221 tot_x=(*p)->momentum().px();
222 tot_y=(*p)->momentum().py();
223 tot_z=(*p)->momentum().pz();
224 tot_e=(*p)->momentum().e();
225 charge1 =
charge((*p)->pdg_id());
228 tot_x += typeApassed[
i]->momentum().px();
229 tot_y += typeApassed[
i]->momentum().py();
230 tot_z += typeApassed[
i]->momentum().pz();
231 tot_e += typeApassed[
i]->momentum().e();
232 invmass=
sqrt(tot_e*tot_e-tot_x*tot_x-tot_y*tot_y-tot_z*tot_z);
233 combcharge = charge1 *
charge(typeApassed[i]->pdg_id());
235 phi2 = typeApassed[
i]->momentum().phi();
236 deltaphi = fabs(phi1-phi2);
237 if(deltaphi >
pi) deltaphi = 2.*
pi-deltaphi;
239 eta2 = typeApassed[
i]->momentum().eta();
240 deltaeta=fabs(eta1-eta2);
241 deltaR =
sqrt(deltaeta*deltaeta+deltaphi*deltaphi);
255 typeBpassed.push_back(*
p);
261 if (accepted){
return true; }
else {
return false;}
273 int kqa,kq1,kq2,kq3,kqj,irt,kqx,kqn;
277 int ichg[109]={-1,2,-1,2,-1,2,-1,2,0,0,-3,0,-3,0,-3,0,
278 -3,0,0,0,0,0,0,3,0,0,0,0,0,0,3,0,3,6,0,0,3,6,0,0,-1,2,-1,2,-1,2,0,0,0,0,
279 -3,0,-3,0,-3,0,0,0,0,0,-1,2,-1,2,-1,2,0,0,0,0,
280 -3,0,-3,0,-3,0,3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
286 kqn=kqa/1000000000%10;
296 if(kqa==0 || kqa >= 10000000) {
298 if(kqn==1) {hepchg=0;}
301 else if(kqa<=100) {hepchg = ichg[kqa-1];}
303 else if(kqj == 0) {hepchg = 0;}
305 else if(kqx>0 && irt<100)
307 hepchg = ichg[irt-1];
308 if(kqa==1000017 || kqa==1000018) {hepchg = 0;}
309 if(kqa==1000034 || kqa==1000052) {hepchg = 0;}
310 if(kqa==1000053 || kqa==1000054) {hepchg = 0;}
311 if(kqa==5100061 || kqa==5100062) {hepchg = 6;}
317 hepchg = ichg[kq2-1]-ichg[kq1-1];
319 if((kq2==3) || (kq2==5)) {hepchg = ichg[kq1-1]-ichg[kq2-1];}
323 hepchg = ichg[kq3-1] + ichg[kq2-1];
328 hepchg = ichg[kq3-1]+ichg[kq2-1]+ichg[kq1-1];
332 if(Id<0 && hepchg!=0) {hepchg = -1*hepchg;}
T getUntrackedParameter(std::string const &, T const &) const
std::vector< double > ptMin
std::vector< double > pMin
std::vector< double > etaMin
virtual bool filter(edm::Event &, const edm::EventSetup &)
MCParticlePairFilter(const edm::ParameterSet &)
std::vector< int > particleID1
std::vector< double > etaMax
std::vector< int > particleID2
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
double deltaR(double eta1, double eta2, double phi1, double phi2)
int charge(const int &Id)
std::vector< int > status