81 if ( *(evt_orig->vertices_begin()) == 0 )
83 throw SimG4Exception(
"SimG4CoreGenerator: Corrupted Event - GenEvent with no vertex" ) ;
87 HepMC::GenEvent* evt =
new HepMC::GenEvent(*evt_orig) ;
89 if ( evt->weights().size() > 0 )
92 for (
unsigned int iw=1; iw<evt->weights().size(); iw++ )
95 if ( evt->weights()[iw] <= 0 )
break;
102 (*(evt->vertices_begin()))->position().y(),
103 (*(evt->vertices_begin()))->position().z(),
104 (*(evt->vertices_begin()))->position().t());
108 LogDebug(
"SimG4CoreGenerator") <<
"Primary Vertex = (" <<
vtx_->x() <<
","
116 unsigned int ng4vtx = 0;
118 for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
119 vitr != evt->vertices_end(); ++vitr ) {
124 for (HepMC::GenVertex::particle_iterator pitr= (*vitr)->particles_begin(HepMC::children);
125 pitr != (*vitr)->particles_end(HepMC::children); ++pitr) {
127 if ((*pitr)->status()==1) {
130 <<
"GenVertex barcode = " << (*vitr)->barcode()
131 <<
" selected for GenParticle barcode = " << (*pitr)->barcode()
138 else if ( (*pitr)->status()== 2 ) {
139 if ( (*pitr)->end_vertex() != 0 ) {
142 double xx = (*pitr)->end_vertex()->position().x();
143 double yy = (*pitr)->end_vertex()->position().y();
148 <<
"GenVertex barcode = " << (*vitr)->barcode()
149 <<
" selected for GenParticle barcode = "
150 << (*pitr)->barcode() <<
" radius = " << r_dd << std::endl;
162 double x1 = (*vitr)->position().x();
163 double y1 = (*vitr)->position().y();
164 double z1 = (*vitr)->position().z();
165 double t1 = (*vitr)->position().t();
166 G4PrimaryVertex* g4vtx=
167 new G4PrimaryVertex(x1*mm, y1*mm, z1*mm, t1*mm/c_light);
169 for (HepMC::GenVertex::particle_iterator vpitr= (*vitr)->particles_begin(HepMC::children);
170 vpitr != (*vitr)->particles_end(HepMC::children); ++vpitr){
180 double r_decay_length=-1;
181 double decay_length=-1;
183 if ( (*vpitr)->status() == 1 || (*vpitr)->status() == 2 ) {
185 if ( (*vpitr)->end_vertex() != 0 ) {
190 double x2 = (*vpitr)->end_vertex()->position().x();
191 double y2 = (*vpitr)->end_vertex()->position().y();
192 double z2 = (*vpitr)->end_vertex()->position().z();
194 decay_length=
std::sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2));
200 bool toBeAdded =
false;
202 (*vpitr)->momentum().py(),
203 (*vpitr)->momentum().pz(),
204 (*vpitr)->momentum().e());
207 const double minTan = 1.e-20;
209 if ( fabs(z1) <
Z_hector && fabs(
p.pt()/
p.pz()) >= minTan ) {
215 <<
"Processing GenParticle barcode = " << (*vpitr)->barcode()
216 <<
" status = " << (*vpitr)->status()
217 <<
" zimpact = " << zimpact;
226 <<
"Skip GenParticle barcode = " << (*vpitr)->barcode()
227 <<
" PDG Id = " << (*vpitr)->pdg_id() << std::endl;
234 if( (*vpitr)->status() == 1 && fabs(zimpact) <
Z_hector ) {
240 <<
"GenParticle barcode = " << (*vpitr)->barcode()
241 <<
" passed case 1" << std::endl;
245 else if((*vpitr)->status() == 2 && r_decay_length >
theRDecLenCut &&
249 <<
"GenParticle barcode = " << (*vpitr)->barcode()
250 <<
" passed case 2" << std::endl;
255 else if( (*vpitr)->status() == 1 && fabs(zimpact) >=
Z_hector && fabs(z1) >=
Z_hector) {
258 <<
"GenParticle barcode = " << (*vpitr)->barcode()
259 <<
" passed case 3" << std::endl;
264 G4int pdgcode= (*vpitr)-> pdg_id();
265 G4PrimaryParticle* g4prim=
266 new G4PrimaryParticle(pdgcode,
p.Px()*GeV,
p.Py()*GeV,
p.Pz()*GeV);
268 if ( g4prim->GetG4code() != 0 ){
269 g4prim->SetMass( g4prim->GetG4code()->GetPDGMass() ) ;
270 g4prim->SetCharge( g4prim->GetG4code()->GetPDGCharge() ) ;
276 setGenId( g4prim, (*vpitr)->barcode() ) ;
278 if ( (*vpitr)->status() == 2 ) {
282 if (
verbose > 1 ) g4prim->Print();
283 g4vtx->SetPrimary(g4prim);
285 if ( (*vpitr)->status()==1 && (*vpitr)->end_vertex()!=0) {
286 double proper_time=decay_length/(
p.Beta()*
p.Gamma()*c_light);
288 <<
"Setting proper time for beta="<<
p.Beta()<<
" gamma="
289 <<
p.Gamma()<<
" Proper time=" <<proper_time/ns<<
" ns" ;
290 g4prim->SetProperTime(proper_time);
295 if (
verbose > 1 ) g4vtx->Print();
296 g4evt->AddPrimaryVertex(g4vtx);
306 G4PrimaryVertex* g4vtx=
307 new G4PrimaryVertex(x1*mm, y1*mm, z1*mm, t1*mm/c_light);
308 if (
verbose > 1 ) g4vtx->Print();
309 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_
return(e1-e2)*(e1-e2)+dp *dp
std::vector< int > pdgFilter
bool particlePassesPrimaryCuts(const G4PrimaryParticle *p) const