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

Jump to bottom

Allow a boost on particles in MCSingleParticle Filters [80X Backport] #18535

Merged
merged 3 commits into from
May 29, 2017
Merged

Allow a boost on particles in MCSingleParticle Filters [80X Backport] #18535

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
22ccc08
Included the betaboost in the MCSingleParticle filters
stahlleiton May 1, 2017
28286f1
return error message if input beta is >= 1
stahlleiton May 2, 2017
ff3f79e
Use std::abs and call the boost function only for particles that pass…
stahlleiton May 5, 2017
File filter

Filter by extension

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

class MCSingleParticleFilter : public edm::EDFilter {
public:
Expand All @@ -47,6 +50,9 @@ class MCSingleParticleFilter : public edm::EDFilter {

virtual bool filter(edm::Event&, const edm::EventSetup&);
private:
// ----------memeber function----------------------
HepMC::FourVector zboost(const HepMC::FourVector&);

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

edm::EDGetTokenT<edm::HepMCProduct> token_;
Expand All @@ -55,5 +61,6 @@ class MCSingleParticleFilter : public edm::EDFilter {
std::vector<double> etaMin;
std::vector<double> etaMax;
std::vector<int> status;
double betaBoost;
};
#endif
7 changes: 7 additions & 0 deletions GeneratorInterface/GenFilters/interface/MCSmartSingleParticleFilter.h
View file
Open in desktop
Expand Up @@ -35,6 +35,9 @@
namespace edm {
class HepMCProduct;
}
namespace HepMC {
class FourVector;
}

class MCSmartSingleParticleFilter : public edm::EDFilter {
public:
Expand All @@ -44,6 +47,9 @@ class MCSmartSingleParticleFilter : public edm::EDFilter {

virtual bool filter(edm::Event&, const edm::EventSetup&);
private:
// ----------memeber function----------------------
HepMC::FourVector zboost(const HepMC::FourVector&);

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

edm::EDGetTokenT<edm::HepMCProduct> token_;
Expand All @@ -57,5 +63,6 @@ class MCSmartSingleParticleFilter : public edm::EDFilter {
std::vector<double> decayRadiusMax;
std::vector<double> decayZMin;
std::vector<double> decayZMax;
double betaBoost;
};
#endif
32 changes: 26 additions & 6 deletions GeneratorInterface/GenFilters/src/MCSingleParticleFilter.cc
View file
Open in desktop
Expand Up @@ -9,7 +9,8 @@ using namespace std;


MCSingleParticleFilter::MCSingleParticleFilter(const edm::ParameterSet& iConfig) :
token_(consumes<edm::HepMCProduct>(edm::InputTag(iConfig.getUntrackedParameter("moduleLabel",std::string("generator")),"unsmeared")))
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 ;
Expand Down Expand Up @@ -63,7 +64,10 @@ token_(consumes<edm::HepMCProduct>(edm::InputTag(iConfig.getUntrackedParameter("
status = defstat2;
}


// check if beta is smaller than 1
if (std::abs(betaBoost) >= 1 ){
edm::LogError("MCSingleParticleFilter") << "Input beta boost is >= 1 !";
}

}

Expand Down Expand Up @@ -95,10 +99,15 @@ bool MCSingleParticleFilter::filter(edm::Event& iEvent, const edm::EventSetup& i
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)->momentum().eta() > etaMin[i]
&& (*p)->momentum().eta() < etaMax[i] && ((*p)->status() == status[i] || status[i] == 0)) {
accepted = true;
}
if ( (*p)->momentum().perp() > ptMin[i]
&& ((*p)->status() == status[i] || status[i] == 0)) {

HepMC::FourVector mom = zboost((*p)->momentum());
if ( mom.eta() > etaMin[i] && mom.eta() < etaMax[i] ) {
accepted = true;
}

}

}
}
Expand All @@ -110,3 +119,14 @@ bool MCSingleParticleFilter::filter(edm::Event& iEvent, const edm::EventSetup& i

}


HepMC::FourVector MCSingleParticleFilter::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));
}

57 changes: 39 additions & 18 deletions GeneratorInterface/GenFilters/src/MCSmartSingleParticleFilter.cc
View file
Open in desktop
Expand Up @@ -9,7 +9,8 @@ using namespace std;


MCSmartSingleParticleFilter::MCSmartSingleParticleFilter(const edm::ParameterSet& iConfig) :
token_(consumes<edm::HepMCProduct>(edm::InputTag(iConfig.getUntrackedParameter("moduleLabel",std::string("generator")),"unsmeared")))
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 ;
Expand Down Expand Up @@ -119,6 +120,10 @@ token_(consumes<edm::HepMCProduct>(edm::InputTag(iConfig.getUntrackedParameter("
decayZMax = decayZmax2;
}

// check if beta is smaller than 1
if (std::abs(betaBoost) >= 1 ){
edm::LogError("MCSmartSingleParticleFilter") << "Input beta boost is >= 1 !";
}

}

Expand Down Expand Up @@ -149,25 +154,30 @@ bool MCSmartSingleParticleFilter::filter(edm::Event& iEvent, const edm::EventSet

for (unsigned int i = 0; i < particleID.size(); i++){
if (particleID[i] == (*p)->pdg_id() || particleID[i] == 0) {

if ( (*p)->momentum().rho() > pMin[i] && (*p)->momentum().perp() > ptMin[i]
&& (*p)->momentum().eta() > etaMin[i] && (*p)->momentum().eta() < etaMax[i]
&& ((*p)->status() == status[i] || status[i] == 0)) {

if (!((*p)->production_vertex())) continue;
if ( (*p)->momentum().perp() > ptMin[i]
&& ((*p)->status() == status[i] || status[i] == 0)) {

HepMC::FourVector mom = zboost((*p)->momentum());
if ( mom.rho() > pMin[i]
&& mom.eta() > etaMin[i] && mom.eta() < etaMax[i]) {

if (!((*p)->production_vertex())) continue;

double decx = (*p)->production_vertex()->position().x();
double decy = (*p)->production_vertex()->position().y();
double decrad = sqrt(decx*decx+decy*decy);
if (decrad<decayRadiusMin[i] ) continue;
if (decrad>decayRadiusMax[i] ) continue;

double decz = (*p)->production_vertex()->position().z();
if (decz<decayZMin[i] ) continue;
if (decz>decayZMax[i] ) continue;

accepted = true;
}
double decx = (*p)->production_vertex()->position().x();
double decy = (*p)->production_vertex()->position().y();
double decrad = sqrt(decx*decx+decy*decy);
if (decrad<decayRadiusMin[i] ) continue;
if (decrad>decayRadiusMax[i] ) continue;

double decz = (*p)->production_vertex()->position().z();
if (decz<decayZMin[i] ) continue;
if (decz>decayZMax[i] ) continue;

accepted = true;
}

}

}
}
Expand All @@ -179,3 +189,14 @@ bool MCSmartSingleParticleFilter::filter(edm::Event& iEvent, const edm::EventSet

}


HepMC::FourVector MCSmartSingleParticleFilter::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));
}