Skip to content

Commit

Permalink
Add parameter for Z boost
Browse files Browse the repository at this point in the history
  • Loading branch information
@echapon
echapon committed Nov 2, 2016
1 parent 51edd2c commit bf9d278
Show file tree
Hide file tree
Showing 2 changed files with 50 additions and 30 deletions.
5 changes: 5 additions & 0 deletions GeneratorInterface/GenFilters/interface/MCParticlePairFilter.h
View file
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
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));
}

0 comments on commit bf9d278

Please sign in to comment.