Skip to content

Commit

Permalink
Merge pull request #20931 from civanch/exclude_gluons_from_g4tracking
Browse files Browse the repository at this point in the history 
Excluded gluons and quarks from Geant4 tracking
  • Loading branch information
@cmsbuild
cmsbuild committed Jan 12, 2018
2 parents 2be4c4b + 54a4f2b commit 173184a
Show file tree
Hide file tree
Showing 3 changed files with 103 additions and 76 deletions.
9 changes: 6 additions & 3 deletions SimG4Core/Application/python/g4SimHits_cfi.py
View file
Expand Up @@ -155,7 +155,6 @@
),
Generator = cms.PSet(
HectorEtaCut,
# string HepMCProductLabel = "generatorSmeared"
HepMCProductLabel = cms.InputTag('generatorSmeared'),
ApplyPCuts = cms.bool(True),
ApplyPtransCut = cms.bool(False),
Expand All @@ -168,9 +167,13 @@
LDecLenCut = cms.double(30.0), ## (cm) decay volume length
ApplyPhiCuts = cms.bool(False),
MinPhiCut = cms.double(-3.14159265359), ## (radians)
MaxPhiCut = cms.double(3.14159265359), ## according to CMS conventions
MaxPhiCut = cms.double(3.14159265359), ## according to CMS conventions
ApplyLumiMonitorCuts = cms.bool(False), ## primary for lumi monitors
Verbosity = cms.untracked.int32(0)
Verbosity = cms.untracked.int32(0),
PDGselection = cms.PSet(
PDGfilterSel = cms.bool(False), ## filter out unwanted particles
PDGfilter = cms.vint32(21,1,2,3,4,5,6) ## list of unwanted particles (gluons and quarks)
)
),
RunAction = cms.PSet(
StopFile = cms.string('StopRun')
Expand Down
5 changes: 3 additions & 2 deletions SimG4Core/Generators/interface/Generator.h
View file
Expand Up @@ -33,8 +33,9 @@ class Generator
private:

bool particlePassesPrimaryCuts(const G4ThreeVector& p) const;
bool isExotic(HepMC::GenParticle* p) const;
bool isExoticNonDetectable(HepMC::GenParticle* p) const;
bool isExotic(int pdgcode) const;
bool isExoticNonDetectable(int pdgcode) const;
bool IsInTheFilterList(int pdgcode) const;
void particleAssignDaughters(G4PrimaryParticle * p, HepMC::GenParticle * hp,
double length);
void setGenId(G4PrimaryParticle* p, int id) const
Expand Down
165 changes: 94 additions & 71 deletions SimG4Core/Generators/src/Generator.cc
View file
Expand Up @@ -17,9 +17,9 @@
#include "G4PhysicalConstants.hh"
#include "G4Log.hh"

#include <sstream>

using namespace edm;
//using std::cout;
//using std::endl;

Generator::Generator(const ParameterSet & p) :
fPCuts(p.getParameter<bool>("ApplyPCuts")),
Expand All @@ -42,7 +42,7 @@ Generator::Generator(const ParameterSet & p) :
Z_lmax(0),
Z_hector(0),
pdgFilterSel(false),
fPDGFilter(false)
fPDGFilter(false)
{
bool lumi = p.getParameter<bool>("ApplyLumiMonitorCuts");
if(lumi) { fLumiFilter = new LumiMonitorFilter(); }
Expand All @@ -57,20 +57,22 @@ Generator::Generator(const ParameterSet & p) :
pdgFilter.resize(0);
if ( p.exists("PDGselection") ) {
pdgFilterSel = (p.getParameter<edm::ParameterSet>("PDGselection")).
getParameter<bool>("PDGfilterSel");
getParameter<bool>("PDGfilterSel");
pdgFilter = (p.getParameter<edm::ParameterSet>("PDGselection")).
getParameter<std::vector< int > >("PDGfilter");
getParameter<std::vector< int > >("PDGfilter");
if(!pdgFilter.empty()) {
fPDGFilter = true;
std::stringstream ss;
ss << "SimG4Core/Generator: ";
if (pdgFilterSel) {
ss << " Selecting only PDG ID = ";
} else {
ss << " Filtering out PDG ID = ";
}
for ( unsigned int ii = 0; ii < pdgFilter.size(); ++ii) {
if (pdgFilterSel) {
edm::LogWarning("SimG4CoreGenerator")
<< " *** Selecting only PDG ID = " << pdgFilter[ii];
} else {
edm::LogWarning("SimG4CoreGenerator")
<< " *** Filtering out PDG ID = " << pdgFilter[ii];
}
ss << pdgFilter[ii] << " ";
}
edm::LogInfo("SimG4CoreGenerator") << ss.str();
}
}

Expand Down Expand Up @@ -105,7 +107,10 @@ void Generator::HepMC2G4(const HepMC::GenEvent * evt_orig, G4Event * g4evt)
HepMC::GenEvent *evt=new HepMC::GenEvent(*evt_orig);

if ( *(evt->vertices_begin()) == nullptr ) {
throw SimG4Exception("SimG4CoreGenerator: Corrupted Event - GenEvent with no vertex");
std::stringstream ss;
ss << "SimG4Core/Generator: in event " << g4evt->GetEventID()
<< " Corrupted Event - GenEvent with no vertex \n";
throw SimG4Exception(ss.str());
}

if (!evt->weights().empty()) {
Expand Down Expand Up @@ -152,7 +157,8 @@ void Generator::HepMC2G4(const HepMC::GenEvent * evt_orig, G4Event * g4evt)
// 2: particles are decayed by generator but need to be propagated by GEANT
// 3: particles are decayed by generator but do not need to be propagated by GEANT
int status = (*pitr)->status();
if (status > 3 && isExotic(*pitr) && (!(isExoticNonDetectable(*pitr)))) {
int pdg = (*pitr)->pdg_id();
if (status > 3 && isExotic(pdg) && (!(isExoticNonDetectable(pdg)))) {
// In Pythia 8, there are many status codes besides 1, 2, 3.
// By setting the status to 2 for exotic particles, they will be checked:
// if its decay vertex is outside the beampipe, it will be propagated by GEANT.
Expand All @@ -163,6 +169,12 @@ void Generator::HepMC2G4(const HepMC::GenEvent * evt_orig, G4Event * g4evt)

// Particles which are not decayed by generator
if (status == 1) {

// filter out unwanted particles and vertices
if(fPDGFilter && !pdgFilterSel && IsInTheFilterList(pdg)) {
continue;
}

qvtx = true;
if (verbose > 2) LogDebug("SimG4CoreGenerator")
<< "GenVertex barcode = " << (*vitr)->barcode()
Expand All @@ -189,6 +201,7 @@ void Generator::HepMC2G4(const HepMC::GenEvent * evt_orig, G4Event * g4evt)
// particles with status 2 without end_vertex are
// equivalent to stable
qvtx = true;
break;
}
}
}
Expand All @@ -209,43 +222,57 @@ void Generator::HepMC2G4(const HepMC::GenEvent * evt_orig, G4Event * g4evt)
for (pitr= (*vitr)->particles_begin(HepMC::children);
pitr != (*vitr)->particles_end(HepMC::children); ++pitr){

int status = (*pitr)->status();
int pdg = (*pitr)->pdg_id();
bool hasDecayVertex = (nullptr != (*pitr)->end_vertex()) ? true : false;

// Filter on allowed particle species if required
if (fPDGFilter) {
std::vector<int>::iterator it = find(pdgFilter.begin(),pdgFilter.end(),
(*pitr)->pdg_id());
if ( (it != pdgFilter.end() && !pdgFilterSel)
|| ( it == pdgFilter.end() && pdgFilterSel ) ) {
if (verbose > 2) LogDebug("SimG4CoreGenerator")
<< "Skip GenParticle barcode = " << (*pitr)->barcode()
<< " PDG Id = " << (*pitr)->pdg_id();
if (fPDGFilter) {
bool isInTheList = IsInTheFilterList(pdg);
if((!pdgFilterSel && isInTheList) || (pdgFilterSel && !isInTheList)) {
edm::LogInfo("SimG4CoreGenerator")
<< " Skiped GenParticle barcode= " << (*pitr)->barcode()
<< " PDGid= " << pdg
<< " status= " << status
<< " isExotic: " << isExotic(pdg)
<< " isExoticNotDet: " << isExoticNonDetectable(pdg)
<< " isInTheList: " << isInTheList
<< " hasDecayVertex: " << hasDecayVertex;
continue;
}
}

double x2 = 0.0;
double y2 = 0.0;
double z2 = 0.0;
double decay_length = 0.0;
int status = (*pitr)->status();
edm::LogInfo("SimG4CoreGenerator")
<< " pdg= " << pdg
<< " status= " << status
<< " hasPreDefinedDecay: " << hasDecayVertex
<< " isExotic: " << isExotic(pdg)
<< " isExoticNotDet: " << isExoticNonDetectable(pdg)
<< " isInTheList: " << IsInTheFilterList(pdg);

if (status > 3 && isExotic(*pitr) && (!(isExoticNonDetectable(*pitr))) ) {
status = 2;
if (status > 3 && isExotic(pdg) && (!(isExoticNonDetectable(pdg))) ) {
status = hasDecayVertex ? 2 : 1;
}

// this particle has predefined decay but has no vertex
if (2 == status && !hasDecayVertex) {
edm::LogWarning("SimG4CoreGenerator: in event ")
<< g4evt->GetEventID()
<< " a particle "
<< " pdgid= " << pdg
<< " has status=2 but has no decay vertex, so will be fully tracked by Geant4";
status = 1;
}

// check the status, 2 has end point with decay defined by generator
if (1 == status || 2 == status) {

// this particle has decayed but have no vertex
// it is an exotic case
if ( !(*pitr)->end_vertex() ) {
status = 1;
} else {
x2 = (*pitr)->end_vertex()->position().x();
y2 = (*pitr)->end_vertex()->position().y();
z2 = (*pitr)->end_vertex()->position().z();
decay_length =
std::sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1));
}
double x2 = x1;
double y2 = y1;
double z2 = z1;
double decay_length = 0.0;
if(2 == status) {
x2 = (*pitr)->end_vertex()->position().x();
y2 = (*pitr)->end_vertex()->position().y();
z2 = (*pitr)->end_vertex()->position().z();
decay_length = std::sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1));
}

bool toBeAdded = false;
Expand All @@ -262,7 +289,7 @@ void Generator::HepMC2G4(const HepMC::GenEvent * evt_orig, G4Event * g4evt)
// protection against numerical problems for extremely low momenta
// compute impact point at transition to Hector
const double minTan = 1.e-20;
if(fabs(z1) < Z_hector && fabs(pz) >= minTan*ptot) {
if(std::abs(z1) < Z_hector && std::abs(pz) >= minTan*ptot) {
if(pz > 0.0) { zimpact = Z_hector; }
else { zimpact = -Z_hector; }
double del = (zimpact - z1)/pz;
Expand All @@ -273,7 +300,7 @@ void Generator::HepMC2G4(const HepMC::GenEvent * evt_orig, G4Event * g4evt)

if (verbose > 2) LogDebug("SimG4CoreGenerator")
<< "Processing GenParticle barcode= " << (*pitr)->barcode()
<< " pdg= " << (*pitr)-> pdg_id()
<< " pdg= " << pdg
<< " status= " << (*pitr)->status()
<< " st= " << status
<< " rimpact(cm)= " << std::sqrt(rimpact2)/cm
Expand All @@ -284,7 +311,7 @@ void Generator::HepMC2G4(const HepMC::GenEvent * evt_orig, G4Event * g4evt)
// Particles of status 1 trasnported along the beam pipe for forward
// detectors (HECTOR) always pass to Geant4 without cuts
if( 1 == status &&
fabs(zimpact) >= Z_hector && rimpact2 <= theDecRCut2) {
std::abs(zimpact) >= Z_hector && rimpact2 <= theDecRCut2) {
toBeAdded = true;
if ( verbose > 2 ) LogDebug("SimG4CoreGenerator")
<< "GenParticle barcode = " << (*pitr)->barcode()
Expand All @@ -293,7 +320,7 @@ void Generator::HepMC2G4(const HepMC::GenEvent * evt_orig, G4Event * g4evt)

// Standard case: particles not decayed by the generator
if(1 == status &&
(fabs(zimpact) < Z_hector || rimpact2 > theDecRCut2)) {
(std::abs(zimpact) < Z_hector || rimpact2 > theDecRCut2)) {

// Ptot cut for all particles
if (fPCuts && (ptot < theMinPCut || ptot > theMaxPCut)) {
Expand All @@ -316,7 +343,7 @@ void Generator::HepMC2G4(const HepMC::GenEvent * evt_orig, G4Event * g4evt)
double zi = z1;

// can be propagated along Z
if(fabs(pz) >= minTan*ptot) {
if(std::abs(pz) >= minTan*ptot) {
if((zi >= Z_lmax) & (pz < 0.0)) {
zi = Z_lmax;
} else if((zi <= Z_lmin) & (pz > 0.0)) {
Expand Down Expand Up @@ -346,7 +373,7 @@ void Generator::HepMC2G4(const HepMC::GenEvent * evt_orig, G4Event * g4evt)
// are used for Geant4 tracking with predefined decay channel
// In the case of decay in vacuum particle is not tracked by Geant4
} else if(2 == status &&
x2*x2 + y2*y2 >= theDecRCut2 && fabs(z2) < Z_hector) {
x2*x2 + y2*y2 >= theDecRCut2 && std::abs(z2) < Z_hector) {

toBeAdded = true;
if ( verbose > 2 ) LogDebug("SimG4CoreGenerator")
Expand All @@ -356,9 +383,8 @@ void Generator::HepMC2G4(const HepMC::GenEvent * evt_orig, G4Event * g4evt)
}
}
if(toBeAdded){
G4int pdgcode= (*pitr)-> pdg_id();
G4PrimaryParticle* g4prim=
new G4PrimaryParticle(pdgcode, px*GeV, py*GeV, pz*GeV);
new G4PrimaryParticle(pdg, px*GeV, py*GeV, pz*GeV);

if ( g4prim->GetG4code() != nullptr ){
g4prim->SetMass( g4prim->GetG4code()->GetPDGMass() );
Expand Down Expand Up @@ -513,35 +539,32 @@ bool Generator::particlePassesPrimaryCuts(const G4ThreeVector& p) const
return flag;
}

bool Generator::isExotic(HepMC::GenParticle* p) const
bool Generator::isExotic(int pdgcode) const
{
int pdgid = abs(p->pdg_id());
if ((pdgid >= 1000000 && pdgid < 4000000) || // SUSY, R-hadron, and technicolor particles
pdgid == 17 || // 4th generation lepton
pdgid == 34 || // W-prime
pdgid == 37) // charged Higgs
{
return true;
}

return false;
int pdgid = std::abs(pdgcode);
return ((pdgid >= 1000000 && pdgid < 4000000) || // SUSY, R-hadron, and technicolor particles
pdgid == 17 || // 4th generation lepton
pdgid == 34 || // W-prime
pdgid == 37) // charged Higgs
? true : false;
}

bool Generator::isExoticNonDetectable(HepMC::GenParticle* p) const
bool Generator::isExoticNonDetectable(int pdgcode) const
{
int pdgid = abs(p->pdg_id());
HepPDT::ParticleID pid(p->pdg_id());
int pdgid = std::abs(pdgcode);
HepPDT::ParticleID pid(pdgcode);
int charge = pid.threeCharge();
if ((charge==0) && (pdgid >= 1000000 && pdgid < 1000040)) // SUSY
{
return true;
}
return ((charge==0) && (pdgid >= 1000000 && pdgid < 1000040)) // SUSY
? true : false;
}

bool Generator::IsInTheFilterList(int pdgcode) const
{
int pdgid = std::abs(pdgcode);
for(auto & pdg : pdgFilter) { if(std::abs(pdg) == pdgid) { return true; } }
return false;
}



void Generator::nonBeamEvent2G4(const HepMC::GenEvent * evt, G4Event * g4evt)
{
int i = 0;
Expand Down

0 comments on commit 173184a

Please sign in to comment.