Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Allow a boost on particles in MCParticlePairFilter #16621

Merged
merged 1 commit into from
Jan 17, 2017
+50 −30
Merged

Allow a boost on particles in MCParticlePairFilter #16621

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions GeneratorInterface/GenFilters/interface/MCParticlePairFilter.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,9 @@
namespace edm {
class HepMCProduct;
}
namespace HepMC {
class FourVector;
}

class MCParticlePairFilter : public edm::EDFilter {
public:
Expand All @@ -49,6 +52,7 @@ class MCParticlePairFilter : public edm::EDFilter {
private:
// ----------memeber function----------------------
int charge(const int& Id);
HepMC::FourVector zboost(const HepMC::FourVector&);

// ----------member data ---------------------------

Expand All @@ -67,6 +71,7 @@ class MCParticlePairFilter : public edm::EDFilter {
double maxDeltaPhi;
double minDeltaR;
double maxDeltaR;
double betaBoost;

};
#endif
75 changes: 45 additions & 30 deletions GeneratorInterface/GenFilters/src/MCParticlePairFilter.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,8 @@ maxInvMass(iConfig.getUntrackedParameter("MaxInvMass", 14000.)),
minDeltaPhi(iConfig.getUntrackedParameter("MinDeltaPhi", 0.)),
maxDeltaPhi(iConfig.getUntrackedParameter("MaxDeltaPhi", 6.3)),
minDeltaR(iConfig.getUntrackedParameter("MinDeltaR",0.)),
maxDeltaR(iConfig.getUntrackedParameter("MaxDeltaR",10000.))
maxDeltaR(iConfig.getUntrackedParameter("MaxDeltaR",10000.)),
betaBoost(iConfig.getUntrackedParameter("BetaBoost",0.))
{
//here do whatever other initialization is needed
vector<int> defpid1;
Expand Down Expand Up @@ -125,16 +126,17 @@ bool MCParticlePairFilter::filter(edm::Event& iEvent, const edm::EventSetup& iSe
}
}
if(gottypeAID) {
if ( (*p)->momentum().perp() > ptMin[0] && (*p)->momentum().rho() > pMin[0] && (*p)->momentum().eta() > etaMin[0]
&& (*p)->momentum().eta() < etaMax[0] && ((*p)->status() == status[0] || status[0] == 0)) {
HepMC::FourVector mom = zboost((*p)->momentum());
if ( mom.perp() > ptMin[0] && mom.rho() > pMin[0] && mom.eta() > etaMin[0]
&& mom.eta() < etaMax[0] && ((*p)->status() == status[0] || status[0] == 0)) {
// passed A type conditions ...
// ... now check pair-conditions with B type passed particles
unsigned int i=0;
double deltaphi;
double phi1 = (*p)->momentum().phi();
double phi1 = mom.phi();
double phi2;
double deltaeta;
double eta1 = (*p)->momentum().eta();
double eta1 = mom.eta();
double eta2;
double deltaR;
//HepLorentzVector momentum1 = (*p)->momentum();
Expand All @@ -147,25 +149,26 @@ bool MCParticlePairFilter::filter(edm::Event& iEvent, const edm::EventSetup& iSe
int charge1 = 0;
int combcharge = 0;
while(!accepted && i<typeBpassed.size()) {
tot_x=(*p)->momentum().px();
tot_y=(*p)->momentum().py();
tot_z=(*p)->momentum().pz();
tot_e=(*p)->momentum().e();
tot_x=mom.px();
tot_y=mom.py();
tot_z=mom.pz();
tot_e=mom.e();
charge1 = charge((*p)->pdg_id());
//totmomentum = momentum1 + typeBpassed[i]->momentum();
//invmass = totmomentum.m();
tot_x += typeBpassed[i]->momentum().px();
tot_y += typeBpassed[i]->momentum().py();
tot_z += typeBpassed[i]->momentum().pz();
tot_e += typeBpassed[i]->momentum().e();
HepMC::FourVector mom_i = zboost(typeBpassed[i]->momentum());
tot_x += mom_i.px();
tot_y += mom_i.py();
tot_z += mom_i.pz();
tot_e += mom_i.e();
invmass=sqrt(tot_e*tot_e-tot_x*tot_x-tot_y*tot_y-tot_z*tot_z);
combcharge = charge1 * charge(typeBpassed[i]->pdg_id());
if(invmass > minInvMass && invmass < maxInvMass) {
phi2 = typeBpassed[i]->momentum().phi();
phi2 = mom_i.phi();
deltaphi = fabs(phi1-phi2);
if(deltaphi > pi) deltaphi = 2.*pi-deltaphi;
if(deltaphi > minDeltaPhi && deltaphi < maxDeltaPhi) {
eta2 = typeBpassed[i]->momentum().eta();
eta2 = mom_i.eta();
deltaeta=fabs(eta1-eta2);
deltaR = sqrt(deltaeta*deltaeta+deltaphi*deltaphi);
if(deltaR > minDeltaR && deltaR < maxDeltaR) {
Expand Down Expand Up @@ -195,16 +198,17 @@ bool MCParticlePairFilter::filter(edm::Event& iEvent, const edm::EventSetup& iSe
}
}
if(gottypeBID) {
if ( (*p)->momentum().perp() > ptMin[1] && (*p)->momentum().rho() > pMin[1] && (*p)->momentum().eta() > etaMin[1]
&& (*p)->momentum().eta() < etaMax[1] && ((*p)->status() == status[1] || status[1] == 0)) {
HepMC::FourVector mom = zboost((*p)->momentum());
if ( mom.perp() > ptMin[1] && mom.rho() > pMin[1] && mom.eta() > etaMin[1]
&& mom.eta() < etaMax[1] && ((*p)->status() == status[1] || status[1] == 0)) {
// passed B type conditions ...
// ... now check pair-conditions with A type passed particles vector
unsigned int i=0;
double deltaphi;
double phi1 = (*p)->momentum().phi();
double phi1 = mom.phi();
double phi2;
double deltaeta;
double eta1 = (*p)->momentum().eta();
double eta1 = mom.eta();
double eta2;
double deltaR;
//HepLorentzVector momentum1 = (*p)->momentum();
Expand All @@ -218,25 +222,26 @@ bool MCParticlePairFilter::filter(edm::Event& iEvent, const edm::EventSetup& iSe
int combcharge = 0;
while(!accepted && i<typeApassed.size()) {
if((*p) != typeApassed[i]) {
tot_x=(*p)->momentum().px();
tot_y=(*p)->momentum().py();
tot_z=(*p)->momentum().pz();
tot_e=(*p)->momentum().e();
tot_x=mom.px();
tot_y=mom.py();
tot_z=mom.pz();
tot_e=mom.e();
charge1 = charge((*p)->pdg_id());
//totmomentum = momentum1 + typeApassed[i]->momentum();
//invmass = totmomentum.m();
tot_x += typeApassed[i]->momentum().px();
tot_y += typeApassed[i]->momentum().py();
tot_z += typeApassed[i]->momentum().pz();
tot_e += typeApassed[i]->momentum().e();
//invmass = totmomentum.m();
HepMC::FourVector mom_i = zboost(mom_i);
tot_x += mom_i.px();
tot_y += mom_i.py();
tot_z += mom_i.pz();
tot_e += mom_i.e();
invmass=sqrt(tot_e*tot_e-tot_x*tot_x-tot_y*tot_y-tot_z*tot_z);
combcharge = charge1 * charge(typeApassed[i]->pdg_id());
if(invmass > minInvMass && invmass < maxInvMass) {
phi2 = typeApassed[i]->momentum().phi();
phi2 = mom_i.phi();
deltaphi = fabs(phi1-phi2);
if(deltaphi > pi) deltaphi = 2.*pi-deltaphi;
if(deltaphi > minDeltaPhi && deltaphi < maxDeltaPhi) {
eta2 = typeApassed[i]->momentum().eta();
eta2 = mom_i.eta();
deltaeta=fabs(eta1-eta2);
deltaR = sqrt(deltaeta*deltaeta+deltaphi*deltaphi);
if(deltaR > minDeltaR && deltaR < maxDeltaR) {
Expand Down Expand Up @@ -334,3 +339,13 @@ int MCParticlePairFilter::charge(const int& Id){
// cout << hepchg<< endl;
return hepchg;
}

HepMC::FourVector MCParticlePairFilter::zboost(const HepMC::FourVector& mom) {
//Boost this Lorentz vector (from TLorentzVector::Boost)
double b2 = betaBoost*betaBoost;
double gamma = 1.0 / sqrt(1.0 - b2);
double bp = betaBoost*mom.pz();
double gamma2 = b2 > 0 ? (gamma - 1.0)/b2 : 0.0;

return HepMC::FourVector(mom.px(), mom.py(), mom.pz() + gamma2*bp*betaBoost + gamma*betaBoost*mom.e(), gamma*(mom.e()+bp));
}