|
|
|
#include <Generator.h>
Public Member Functions | |
| virtual const double | eventWeight () const |
| Generator (const edm::ParameterSet &p) | |
| virtual const HepMC::GenEvent * | genEvent () const |
| virtual const math::XYZTLorentzVector * | genVertex () const |
| void | HepMC2G4 (const HepMC::GenEvent *g, G4Event *e) |
| void | nonBeamEvent2G4 (const HepMC::GenEvent *g, G4Event *e) |
| void | setGenEvent (const HepMC::GenEvent *inpevt) |
| virtual | ~Generator () |
Private Member Functions | |
| void | particleAssignDaughters (G4PrimaryParticle *p, HepMC::GenParticle *hp, double length) |
| bool | particlePassesPrimaryCuts (const math::XYZTLorentzVector &mom, const double zimp) const |
| bool | particlePassesPrimaryCuts (const G4PrimaryParticle *p) const |
| void | setGenId (G4PrimaryParticle *p, int id) const |
Private Attributes | |
| HepMC::GenEvent * | evt_ |
| bool | fEtaCuts |
| bool | fPCuts |
| bool | fPhiCuts |
| double | theEtaCutForHector |
| double | theMaxEtaCut |
| double | theMaxPCut |
| double | theMaxPhiCut |
| double | theMinEtaCut |
| double | theMinPCut |
| double | theMinPhiCut |
| double | theRDecLenCut |
| int | verbose |
| math::XYZTLorentzVector * | vtx_ |
| double | weight_ |
| double | Z_hector |
| double | Z_lmax |
| double | Z_lmin |
Definition at line 21 of file Generator.h.
| Generator::Generator | ( | const edm::ParameterSet & | p | ) |
Definition at line 22 of file Generator.cc.
References funct::exp(), fEtaCuts, LogDebug, theEtaCutForHector, theMaxEtaCut, theMinEtaCut, theRDecLenCut, Z_hector, Z_lmax, and Z_lmin.
: fPCuts(p.getParameter<bool>("ApplyPCuts")), fEtaCuts(p.getParameter<bool>("ApplyEtaCuts")), fPhiCuts(p.getParameter<bool>("ApplyPhiCuts")), theMinPhiCut(p.getParameter<double>("MinPhiCut")), // now operates in radians (CMS standard) theMaxPhiCut(p.getParameter<double>("MaxPhiCut")), theMinEtaCut(p.getParameter<double>("MinEtaCut")), theMaxEtaCut(p.getParameter<double>("MaxEtaCut")), theMinPCut(p.getParameter<double>("MinPCut")), // now operates in GeV (CMS standard) theMaxPCut(p.getParameter<double>("MaxPCut")), theRDecLenCut(p.getParameter<double>("RDecLenCut")*cm), theEtaCutForHector(p.getParameter<double>("EtaCutForHector")), verbose(p.getUntrackedParameter<int>("Verbosity",0)), evt_(0), vtx_(0), weight_(0), Z_lmin(0), Z_lmax(0), Z_hector(0) { if(fEtaCuts){ Z_lmax = theRDecLenCut*( ( 1 - exp(-2*theMaxEtaCut) ) / ( 2*exp(-theMaxEtaCut) ) ); Z_lmin = theRDecLenCut*( ( 1 - exp(-2*theMinEtaCut) ) / ( 2*exp(-theMinEtaCut) ) ); } Z_hector = theRDecLenCut*( ( 1 - exp(-2*theEtaCutForHector) ) / ( 2*exp(-theEtaCutForHector) ) ); if ( verbose > 0 ) LogDebug("SimG4CoreGenerator") << "Z_min = " << Z_lmin << " Z_max = " << Z_lmax << " Z_hector = " << Z_hector; }
| Generator::~Generator | ( | ) | [virtual] |
Definition at line 54 of file Generator.cc.
{
}
| virtual const double Generator::eventWeight | ( | ) | const [inline, virtual] |
Definition at line 32 of file Generator.h.
References weight_.
Referenced by RunManager::produce().
{ return weight_; }
| virtual const HepMC::GenEvent* Generator::genEvent | ( | ) | const [inline, virtual] |
Definition at line 30 of file Generator.h.
References evt_.
Referenced by RunManager::produce().
{ return evt_; }
| virtual const math::XYZTLorentzVector* Generator::genVertex | ( | ) | const [inline, virtual] |
Definition at line 31 of file Generator.h.
References vtx_.
Referenced by RunManager::produce().
{ return vtx_; }
| void Generator::HepMC2G4 | ( | const HepMC::GenEvent * | g, |
| G4Event * | e | ||
| ) |
Definition at line 58 of file Generator.cc.
References configurableAnalysis::GenParticle, LogDebug, AlCaHLTBitMon_ParallelJobs::p, particleAssignDaughters(), particlePassesPrimaryCuts(), setGenId(), mathSSE::sqrt(), theRDecLenCut, vtx_, weight_, and Z_hector.
Referenced by RunManager::generateEvent().
{
if ( *(evt_orig->vertices_begin()) == 0 )
{
throw SimG4Exception( "SimG4CoreGenerator: Corrupted Event - GenEvent with no vertex" ) ;
}
HepMC::GenEvent* evt = new HepMC::GenEvent(*evt_orig) ;
if ( evt->weights().size() > 0 )
{
weight_ = evt->weights()[0] ;
for ( unsigned int iw=1; iw<evt->weights().size(); iw++ )
{
// terminate if the versot of weights contains a zero-weight
if ( evt->weights()[iw] <= 0 ) break;
weight_ *= evt->weights()[iw] ;
}
}
if (vtx_ != 0) delete vtx_;
vtx_ = new math::XYZTLorentzVector((*(evt->vertices_begin()))->position().x(),
(*(evt->vertices_begin()))->position().y(),
(*(evt->vertices_begin()))->position().z(),
(*(evt->vertices_begin()))->position().t());
if(verbose >0){
evt->print();
LogDebug("SimG4CoreGenerator") << "Primary Vertex = (" << vtx_->x() << ","
<< vtx_->y() << ","
<< vtx_->z() << ")";
}
// double x0 = vtx_->x();
// double y0 = vtx_->y();
unsigned int ng4vtx = 0;
for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
vitr != evt->vertices_end(); ++vitr ) {
// loop for vertex ...
// real vertex?
G4bool qvtx=false;
for (HepMC::GenVertex::particle_iterator pitr= (*vitr)->particles_begin(HepMC::children);
pitr != (*vitr)->particles_end(HepMC::children); ++pitr) {
// Admit also status=1 && end_vertex for long vertex special decay treatment
if ((*pitr)->status()==1) {
qvtx=true;
if ( verbose > 2 ) LogDebug("SimG4CoreGenerator") << "GenVertex barcode = " << (*vitr)->barcode()
<< " selected for GenParticle barcode = " << (*pitr)->barcode() << std::endl;
break;
}
// The selection is made considering if the partcile with status = 2 have the end_vertex
// with a radius (R) greater then the theRDecLenCut that means: the end_vertex is outside
// the beampipe cilinder (no requirement on the Z of the vertex is applyed).
else if ( (*pitr)->status()== 2 ) {
if ( (*pitr)->end_vertex() != 0 ) {
//double xx = x0-(*pitr)->end_vertex()->position().x();
//double yy = y0-(*pitr)->end_vertex()->position().y();
double xx = (*pitr)->end_vertex()->position().x();
double yy = (*pitr)->end_vertex()->position().y();
double r_dd=std::sqrt(xx*xx+yy*yy);
if (r_dd>theRDecLenCut){
qvtx=true;
if ( verbose > 2 ) LogDebug("SimG4CoreGenerator") << "GenVertex barcode = " << (*vitr)->barcode()
<< " selected for GenParticle barcode = " << (*pitr)->barcode() << " radius = " << r_dd << std::endl;
break;
}
}
}
}
if (!qvtx) {
continue;
}
double x1 = (*vitr)->position().x();
double y1 = (*vitr)->position().y();
double z1 = (*vitr)->position().z();
double t1 = (*vitr)->position().t();
G4PrimaryVertex* g4vtx=
new G4PrimaryVertex(x1*mm, y1*mm, z1*mm, t1*mm/c_light);
for (HepMC::GenVertex::particle_iterator vpitr= (*vitr)->particles_begin(HepMC::children);
vpitr != (*vitr)->particles_end(HepMC::children); ++vpitr){
// Special cases:
// 1) import in Geant4 a full decay chain (e.g. also particles with status == 2)
// from the generator in case the decay radius is larger than theRDecLenCut
// In this case no cuts will be applied to select the particles hat has to be
// processed by geant
// 2) import in Geant4 particles with status == 1 but a final end vertex.
// The time of the vertex is used as the time of flight to be forced for the particle
double r_decay_length=-1;
double decay_length=-1;
if ( (*vpitr)->status() == 1 || (*vpitr)->status() == 2 ) {
// this particle has decayed
if ( (*vpitr)->end_vertex() != 0 ) {
// needed some particles have status 2 and no end_vertex
// Which are the particles with status 2 and not end_vertex, what are suppose to di such kind of particles
// when propagated to geant?
double x2 = (*vpitr)->end_vertex()->position().x();
double y2 = (*vpitr)->end_vertex()->position().y();
double z2 = (*vpitr)->end_vertex()->position().z();
r_decay_length=std::sqrt(x2*x2+y2*y2);
decay_length=std::sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2));
}
}
// end modification
bool toBeAdded = false;
math::XYZTLorentzVector p((*vpitr)->momentum().px(),
(*vpitr)->momentum().py(),
(*vpitr)->momentum().pz(),
(*vpitr)->momentum().e());
// protection against numerical problems for extremely low momenta
const double minTan = 1.e-20;
double zimpact = Z_hector+1.;
if ( fabs(z1) < Z_hector && fabs(p.pt()/p.pz()) >= minTan ) {
// write tan(p.Theta()) as p.Pt()/p.Pz()
zimpact = (theRDecLenCut-sqrt(x1*x1+y1*y1))*(p.pz()/p.pt())+z1;
}
if ( verbose > 2 ) LogDebug("SimG4CoreGenerator") << "Processing GenParticle barcode = " << (*vpitr)->barcode()
<< " status = " << (*vpitr)->status()
<< " zimpact = " << zimpact;
// Standard case: particles not decayed by the generator
if( (*vpitr)->status() == 1 && fabs(zimpact) < Z_hector ) {
if ( !particlePassesPrimaryCuts( p, zimpact ) ) {
continue ;
}
toBeAdded = true;
if ( verbose > 2 ) LogDebug("SimG4CoreGenerator") << "GenParticle barcode = " << (*vpitr)->barcode() << " passed case 1" << std::endl;
}
// Decay chain entering exiting the fiducial cylinder defined by theRDecLenCut
else if((*vpitr)->status() == 2 && r_decay_length > theRDecLenCut &&
fabs(zimpact) < Z_hector ) {
toBeAdded=true;
if ( verbose > 2 ) LogDebug("SimG4CoreGenerator") << "GenParticle barcode = " << (*vpitr)->barcode() << " passed case 2" << std::endl;
}
// Particles trasnported along the beam pipe for forward detectors (HECTOR)
// Always pass to Geant4 without cuts (to be checked)
else if( (*vpitr)->status() == 1 && fabs(zimpact) >= Z_hector && fabs(z1) >= Z_hector) {
toBeAdded = true;
if ( verbose > 2 ) LogDebug("SimG4CoreGenerator") << "GenParticle barcode = " << (*vpitr)->barcode() << " passed case 3" << std::endl;
}
if(toBeAdded){
G4int pdgcode= (*vpitr)-> pdg_id();
G4PrimaryParticle* g4prim=
new G4PrimaryParticle(pdgcode, p.Px()*GeV, p.Py()*GeV, p.Pz()*GeV);
if ( g4prim->GetG4code() != 0 ){
g4prim->SetMass( g4prim->GetG4code()->GetPDGMass() ) ;
g4prim->SetCharge( g4prim->GetG4code()->GetPDGCharge() ) ;
}
g4prim->SetWeight( 10000*(*vpitr)->barcode() ) ;
setGenId( g4prim, (*vpitr)->barcode() ) ;
if ( (*vpitr)->status() == 2 )
particleAssignDaughters(g4prim,(HepMC::GenParticle *) *vpitr, decay_length);
if ( verbose > 1 ) g4prim->Print();
g4vtx->SetPrimary(g4prim);
// impose also proper time for status=1 and available end_vertex
if ( (*vpitr)->status()==1 && (*vpitr)->end_vertex()!=0) {
double proper_time=decay_length/(p.Beta()*p.Gamma()*c_light);
if ( verbose > 1 ) LogDebug("SimG4CoreGenerator") <<"Setting proper time for beta="<<p.Beta()<<" gamma="<<p.Gamma()<<" Proper time=" <<proper_time<<" ns" ;
g4prim->SetProperTime(proper_time*ns);
}
}
}
if (verbose > 1 ) g4vtx->Print();
g4evt->AddPrimaryVertex(g4vtx);
ng4vtx++;
}
// Add a protection for completely empty events (produced by LHCTransport): add a dummy vertex with no particle attached to it
if ( ng4vtx == 0 ) {
double x1 = 0.;
double y1 = 0.;
double z1 = 0.;
double t1 = 0.;
G4PrimaryVertex* g4vtx=
new G4PrimaryVertex(x1*mm, y1*mm, z1*mm, t1*mm/c_light);
if (verbose > 1 ) g4vtx->Print();
g4evt->AddPrimaryVertex(g4vtx);
}
delete evt ;
return ;
}
| void Generator::nonBeamEvent2G4 | ( | const HepMC::GenEvent * | g, |
| G4Event * | e | ||
| ) |
Definition at line 368 of file Generator.cc.
References g, configurableAnalysis::GenParticle, i, particlePassesPrimaryCuts(), setGenId(), and v.
Referenced by RunManager::generateEvent().
{
int i = 0;
for(HepMC::GenEvent::particle_const_iterator it = evt->particles_begin();
it != evt->particles_end(); ++it )
{
i++;
HepMC::GenParticle * g = (*it);
int g_status = g->status();
// storing only particle with status == 1
if (g_status == 1)
{
int g_id = g->pdg_id();
G4PrimaryParticle * g4p =
new G4PrimaryParticle(g_id,g->momentum().px()*GeV,g->momentum().py()*GeV,g->momentum().pz()*GeV);
if (g4p->GetG4code() != 0)
{
g4p->SetMass(g4p->GetG4code()->GetPDGMass());
g4p->SetCharge(g4p->GetG4code()->GetPDGCharge()) ;
}
g4p->SetWeight(i*10000);
setGenId(g4p,i);
if (particlePassesPrimaryCuts(g4p))
{
G4PrimaryVertex * v = new
G4PrimaryVertex(g->production_vertex()->position().x()*mm,
g->production_vertex()->position().y()*mm,
g->production_vertex()->position().z()*mm,
g->production_vertex()->position().t()*mm/c_light);
v->SetPrimary(g4p);
g4evt->AddPrimaryVertex(v);
if(verbose >0) {
v->Print();
}
}
}
} // end loop on HepMC particles
}
| void Generator::particleAssignDaughters | ( | G4PrimaryParticle * | p, |
| HepMC::GenParticle * | hp, | ||
| double | length | ||
| ) | [private] |
Definition at line 266 of file Generator.cc.
References createTree::dd, LogDebug, AlCaHLTBitMon_ParallelJobs::p, setGenId(), and mathSSE::sqrt().
Referenced by HepMC2G4().
{
if ( !(vp->end_vertex()) ) return ;
if ( verbose > 1 )
LogDebug("SimG4CoreGenerator") << "Special case of long decay length \n"
<< "Assign daughters with to mother with decaylength=" << decaylength << "mm";
math::XYZTLorentzVector p(vp->momentum().px(), vp->momentum().py(), vp->momentum().pz(), vp->momentum().e());
double proper_time=decaylength/(p.Beta()*p.Gamma()*c_light);
if ( verbose > 2 ) {
LogDebug("SimG4CoreGenerator") <<" px="<<vp->momentum().px()
<<" py="<<vp->momentum().py()
<<" pz="<<vp->momentum().pz()
<<" e="<<vp->momentum().e()
<<" beta="<<p.Beta()
<<" gamma="<<p.Gamma()
<<" Proper time=" <<proper_time<<" ns" ;
}
g4p->SetProperTime(proper_time*ns); // the particle will decay after the same length if it has not interacted before
double x1 = vp->end_vertex()->position().x();
double y1 = vp->end_vertex()->position().y();
double z1 = vp->end_vertex()->position().z();
for (HepMC::GenVertex::particle_iterator
vpdec= vp->end_vertex()->particles_begin(HepMC::children);
vpdec != vp->end_vertex()->particles_end(HepMC::children); ++vpdec) {
//transform decay products such that in the rest frame of mother
math::XYZTLorentzVector pdec((*vpdec)->momentum().px(),
(*vpdec)->momentum().py(),
(*vpdec)->momentum().pz(),
(*vpdec)->momentum().e());
// children should only be taken into account once
G4PrimaryParticle * g4daught=
new G4PrimaryParticle((*vpdec)->pdg_id(), pdec.x()*GeV, pdec.y()*GeV, pdec.z()*GeV);
if ( g4daught->GetG4code() != 0 )
{
g4daught->SetMass( g4daught->GetG4code()->GetPDGMass() ) ;
g4daught->SetCharge( g4daught->GetG4code()->GetPDGCharge() ) ;
}
g4daught->SetWeight( 10000*(*vpdec)->barcode() ) ;
setGenId( g4daught, (*vpdec)->barcode() ) ;
if ( verbose > 1 ) LogDebug("SimG4CoreGenerator") <<"Assigning a "<<(*vpdec)->pdg_id()
<<" as daughter of a " <<vp->pdg_id() ;
if ( (*vpdec)->status() == 2 && (*vpdec)->end_vertex() != 0 )
{
double x2 = (*vpdec)->end_vertex()->position().x();
double y2 = (*vpdec)->end_vertex()->position().y();
double z2 = (*vpdec)->end_vertex()->position().z();
double dd = std::sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2));
particleAssignDaughters(g4daught,*vpdec,dd);
}
(*vpdec)->set_status(1000+(*vpdec)->status());
g4p->SetDaughter(g4daught);
}
return;
}
| bool Generator::particlePassesPrimaryCuts | ( | const math::XYZTLorentzVector & | mom, |
| const double | zimp | ||
| ) | const [private] |
Definition at line 324 of file Generator.cc.
References fEtaCuts, fPCuts, fPhiCuts, LogDebug, phi, theMaxPCut, theMaxPhiCut, and Z_lmax.
{
double phi = mom.Phi() ;
double nrg = mom.P() ;
bool flag = true;
if ( (fEtaCuts) && (zimp < Z_lmin || zimp > Z_lmax) ) {
flag = false;
} else if ( (fPCuts) && (nrg < theMinPCut || nrg > theMaxPCut) ) {
flag = false;
} else if ( (fPhiCuts) && (phi < theMinPhiCut || phi > theMaxPhiCut) ) {
flag = false;
}
if ( verbose > 2 ) LogDebug("SimG4CoreGenerator") << "Generator p = " << nrg << " z_imp = " << zimp << " phi = " << mom.Phi() << " Flag = " << flag;
return flag;
}
| bool Generator::particlePassesPrimaryCuts | ( | const G4PrimaryParticle * | p | ) | const [private] |
Definition at line 344 of file Generator.cc.
References eta(), fEtaCuts, fPCuts, fPhiCuts, funct::log(), LogDebug, phi, mathSSE::sqrt(), funct::tan(), theMaxEtaCut, theMaxPCut, and theMaxPhiCut.
Referenced by HepMC2G4(), and nonBeamEvent2G4().
{
G4ThreeVector mom = p->GetMomentum();
double phi = mom.phi() ;
double nrg = sqrt(p->GetPx()*p->GetPx() + p->GetPy()*p->GetPy() + p->GetPz()*p->GetPz());
nrg /= GeV ; // need to convert, since Geant4 operates in MeV
double eta = -log(tan(mom.theta()/2));
bool flag = true;
if ( (fEtaCuts) && (eta < theMinEtaCut || eta > theMaxEtaCut) ) {
flag = false;
} else if ( (fPCuts) && (nrg < theMinPCut || nrg > theMaxPCut) ) {
flag = false;
} else if ( (fPhiCuts) && (phi < theMinPhiCut || phi > theMaxPhiCut) ) {
flag = false;
}
if ( verbose > 2 ) LogDebug("SimG4CoreGenerator") << "Generator p = " << nrg << " eta = " << eta << " theta = " << mom.theta() << " phi = " << phi << " Flag = " << flag;
return flag;
}
| void Generator::setGenEvent | ( | const HepMC::GenEvent * | inpevt | ) | [inline] |
Definition at line 27 of file Generator.h.
References evt_.
Referenced by RunManager::generateEvent().
{ evt_ = (HepMC::GenEvent*)inpevt; return ; }
| void Generator::setGenId | ( | G4PrimaryParticle * | p, |
| int | id | ||
| ) | const [inline, private] |
Definition at line 37 of file Generator.h.
Referenced by HepMC2G4(), nonBeamEvent2G4(), and particleAssignDaughters().
{p->SetUserInformation(new GenParticleInfo(id));}
HepMC::GenEvent* Generator::evt_ [private] |
Definition at line 53 of file Generator.h.
Referenced by genEvent(), and setGenEvent().
bool Generator::fEtaCuts [private] |
Definition at line 42 of file Generator.h.
Referenced by Generator(), and particlePassesPrimaryCuts().
bool Generator::fPCuts [private] |
Definition at line 41 of file Generator.h.
Referenced by particlePassesPrimaryCuts().
bool Generator::fPhiCuts [private] |
Definition at line 43 of file Generator.h.
Referenced by particlePassesPrimaryCuts().
double Generator::theEtaCutForHector [private] |
Definition at line 51 of file Generator.h.
Referenced by Generator().
double Generator::theMaxEtaCut [private] |
Definition at line 47 of file Generator.h.
Referenced by Generator(), and particlePassesPrimaryCuts().
double Generator::theMaxPCut [private] |
Definition at line 49 of file Generator.h.
Referenced by particlePassesPrimaryCuts().
double Generator::theMaxPhiCut [private] |
Definition at line 45 of file Generator.h.
Referenced by particlePassesPrimaryCuts().
double Generator::theMinEtaCut [private] |
Definition at line 46 of file Generator.h.
Referenced by Generator().
double Generator::theMinPCut [private] |
Definition at line 48 of file Generator.h.
double Generator::theMinPhiCut [private] |
Definition at line 44 of file Generator.h.
double Generator::theRDecLenCut [private] |
Definition at line 50 of file Generator.h.
Referenced by Generator(), and HepMC2G4().
int Generator::verbose [private] |
Definition at line 52 of file Generator.h.
math::XYZTLorentzVector* Generator::vtx_ [private] |
Definition at line 54 of file Generator.h.
Referenced by genVertex(), and HepMC2G4().
double Generator::weight_ [private] |
Definition at line 55 of file Generator.h.
Referenced by eventWeight(), and HepMC2G4().
double Generator::Z_hector [private] |
Definition at line 56 of file Generator.h.
Referenced by Generator(), and HepMC2G4().
double Generator::Z_lmax [private] |
Definition at line 56 of file Generator.h.
Referenced by Generator(), and particlePassesPrimaryCuts().
double Generator::Z_lmin [private] |
Definition at line 56 of file Generator.h.
Referenced by Generator().
1.7.3