77 if ( *(evt_orig->vertices_begin()) == 0 )
79 throw SimG4Exception(
"SimG4CoreGenerator: Corrupted Event - GenEvent with no vertex" ) ;
83 HepMC::GenEvent* evt =
new HepMC::GenEvent(*evt_orig) ;
85 if ( evt->weights().size() > 0 )
88 for (
unsigned int iw=1; iw<evt->weights().size(); iw++ )
91 if ( evt->weights()[iw] <= 0 )
break;
98 (*(evt->vertices_begin()))->position().y(),
99 (*(evt->vertices_begin()))->position().z(),
100 (*(evt->vertices_begin()))->position().t());
104 LogDebug(
"SimG4CoreGenerator") <<
"Primary Vertex = (" <<
vtx_->x() <<
","
112 unsigned int ng4vtx = 0;
114 for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
115 vitr != evt->vertices_end(); ++vitr ) {
120 for (HepMC::GenVertex::particle_iterator pitr= (*vitr)->particles_begin(HepMC::children);
121 pitr != (*vitr)->particles_end(HepMC::children); ++pitr) {
123 if ((*pitr)->status()==1) {
125 if (
verbose > 2 )
LogDebug(
"SimG4CoreGenerator") <<
"GenVertex barcode = " << (*vitr)->barcode()
126 <<
" selected for GenParticle barcode = " << (*pitr)->barcode() << std::endl;
132 else if ( (*pitr)->status()== 2 ) {
133 if ( (*pitr)->end_vertex() != 0 ) {
136 double xx = (*pitr)->end_vertex()->position().x();
137 double yy = (*pitr)->end_vertex()->position().y();
141 if (
verbose > 2 )
LogDebug(
"SimG4CoreGenerator") <<
"GenVertex barcode = " << (*vitr)->barcode()
142 <<
" selected for GenParticle barcode = " << (*pitr)->barcode() <<
" radius = " << r_dd << std::endl;
154 double x1 = (*vitr)->position().x();
155 double y1 = (*vitr)->position().y();
156 double z1 = (*vitr)->position().z();
157 double t1 = (*vitr)->position().t();
158 G4PrimaryVertex* g4vtx=
159 new G4PrimaryVertex(x1*mm, y1*mm, z1*mm, t1*mm/c_light);
161 for (HepMC::GenVertex::particle_iterator vpitr= (*vitr)->particles_begin(HepMC::children);
162 vpitr != (*vitr)->particles_end(HepMC::children); ++vpitr){
172 double r_decay_length=-1;
173 double decay_length=-1;
175 if ( (*vpitr)->status() == 1 || (*vpitr)->status() == 2 ) {
177 if ( (*vpitr)->end_vertex() != 0 ) {
182 double x2 = (*vpitr)->end_vertex()->position().x();
183 double y2 = (*vpitr)->end_vertex()->position().y();
184 double z2 = (*vpitr)->end_vertex()->position().z();
186 decay_length=
std::sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2));
192 bool toBeAdded =
false;
194 (*vpitr)->momentum().py(),
195 (*vpitr)->momentum().pz(),
196 (*vpitr)->momentum().e());
199 const double minTan = 1.e-20;
201 if ( fabs(z1) <
Z_hector && fabs(
p.pt()/
p.pz()) >= minTan ) {
206 if (
verbose > 2 )
LogDebug(
"SimG4CoreGenerator") <<
"Processing GenParticle barcode = " << (*vpitr)->barcode()
207 <<
" status = " << (*vpitr)->status()
208 <<
" zimpact = " << zimpact;
216 if (
verbose > 2 )
LogDebug(
"SimG4CoreGenerator") <<
"Skip GenParticle barcode = " << (*vpitr)->barcode()
217 <<
" PDG Id = " << (*vpitr)->pdg_id() << std::endl;
224 if( (*vpitr)->status() == 1 && fabs(zimpact) <
Z_hector ) {
229 if (
verbose > 2 )
LogDebug(
"SimG4CoreGenerator") <<
"GenParticle barcode = " << (*vpitr)->barcode() <<
" passed case 1" << std::endl;
233 else if((*vpitr)->status() == 2 && r_decay_length >
theRDecLenCut &&
236 if (
verbose > 2 )
LogDebug(
"SimG4CoreGenerator") <<
"GenParticle barcode = " << (*vpitr)->barcode() <<
" passed case 2" << std::endl;
241 else if( (*vpitr)->status() == 1 && fabs(zimpact) >=
Z_hector && fabs(z1) >=
Z_hector) {
243 if (
verbose > 2 )
LogDebug(
"SimG4CoreGenerator") <<
"GenParticle barcode = " << (*vpitr)->barcode() <<
" passed case 3" << std::endl;
248 G4int pdgcode= (*vpitr)-> pdg_id();
249 G4PrimaryParticle* g4prim=
250 new G4PrimaryParticle(pdgcode,
p.Px()*GeV,
p.Py()*GeV,
p.Pz()*GeV);
252 if ( g4prim->GetG4code() != 0 ){
253 g4prim->SetMass( g4prim->GetG4code()->GetPDGMass() ) ;
254 g4prim->SetCharge( g4prim->GetG4code()->GetPDGCharge() ) ;
260 setGenId( g4prim, (*vpitr)->barcode() ) ;
262 if ( (*vpitr)->status() == 2 )
264 if (
verbose > 1 ) g4prim->Print();
265 g4vtx->SetPrimary(g4prim);
267 if ( (*vpitr)->status()==1 && (*vpitr)->end_vertex()!=0) {
268 double proper_time=decay_length/(
p.Beta()*
p.Gamma()*c_light);
269 if (
verbose > 1 )
LogDebug(
"SimG4CoreGenerator") <<
"Setting proper time for beta="<<
p.Beta()<<
" gamma="<<
p.Gamma()<<
" Proper time=" <<proper_time<<
" ns" ;
270 g4prim->SetProperTime(proper_time*ns);
275 if (
verbose > 1 ) g4vtx->Print();
276 g4evt->AddPrimaryVertex(g4vtx);
286 G4PrimaryVertex* g4vtx=
287 new G4PrimaryVertex(x1*mm, y1*mm, z1*mm, t1*mm/c_light);
288 if (
verbose > 1 ) g4vtx->Print();
289 g4evt->AddPrimaryVertex(g4vtx);
void setGenId(G4PrimaryParticle *p, int id) const
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
void particleAssignDaughters(G4PrimaryParticle *p, HepMC::GenParticle *hp, double length)
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
math::XYZTLorentzVector * vtx_
std::vector< int > pdgFilter
bool particlePassesPrimaryCuts(const G4PrimaryParticle *p) const