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 #16425

Merged
merged 1 commit into from
Nov 17, 2016
Merged

Allow a boost on particles in MCParticlePairFilter #16425

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 {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is not needed given the way that FourVector is used below - you must have included the class definition somewhere else

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry about the late reply. Would it be OK if instead of this declaration I moved the included header https://github.com/echapon/cmssw/blob/bf9d278f552c725aa16496ed0a02029303bc8b4e/GeneratorInterface/GenFilters/src/MCParticlePairFilter.cc#L6 from the .cc to the .h?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi-
class FourVector;

lets you reference pointers to your class. To return a full FourVector you need an include "FourVector.h" (or where-ever the class is defined). So, this header file should work the same without the class FourVector portion.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've just tried to remove this part, but then it does not compile... I could also explicitly include the header where this class is define: HepMC/SimpleVector.h

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@davidlange6 sorry for the frequent pings... Could you tell me explicitly how I should change the code?

  • removing the "class FourVector;" line is not an option, because it does no longer compile.
  • I could instead include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h" (it is currently in the .cc)
  • I could instead include "HepMC/SimpleVector.h"

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @echapon - i will try to look tomorrow, something is not right in what you see.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @davidlange6 did you have a chance to have a look?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i just did. don't really understand it. will merge and we see in the future how things progress as compilers change

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));
}