19 #include <Geant4/G4Event.hh>
20 #include <Geant4/G4PrimaryParticle.hh>
21 #include <Geant4/G4PrimaryVertex.hh>
22 #include <Geant4/G4VUserPrimaryParticleInformation.hh>
24 #include <Eigen/Dense>
41 : m_TruthInfoContainer(nullptr)
42 , m_LowerKeyPrevExist(0)
43 , m_UpperKeyPrevExist(0)
53 std::cout <<
"PHG4TruthEventAction::EndOfEventAction - unable to find G4TruthInfo node" << std::endl;
71 std::cout <<
"PHG4TruthEventAction::EndOfEventAction - unable to find G4TruthInfo node" << std::endl;
82 std::set<int> savelist;
83 std::set<int> savevtxlist;
85 for (std::set<int>::const_iterator write_iter =
m_WriteSet.begin();
89 std::vector<int> wrttracks;
90 std::vector<int> wrtvtx;
93 int mytrkid = *write_iter;
97 if (savelist.find(mytrkid) != savelist.end())
103 wrttracks.push_back(mytrkid);
110 while (savelist.find(parentid) == savelist.end() && parentid != 0)
113 wrttracks.push_back(parentid);
120 while (wrttracks.begin() != wrttracks.end())
122 savelist.insert(wrttracks.back());
123 wrttracks.pop_back();
126 while (wrtvtx.begin() != wrtvtx.end())
128 savevtxlist.insert(wrtvtx.back());
139 int removed[4] = {0};
142 while (truthiter != truth_range.second)
145 int trackid = truthiter->first;
146 if (savelist.find(trackid) == savelist.end())
165 savevtxlist.insert((truthiter->second)->get_vtx_id());
178 while (vtxiter != vtxrange.second)
181 if (savevtxlist.find(vtxiter->first) == savevtxlist.end())
232 std::cout <<
"PHG4TruthEventAction::SetInterfacePointers - unable to find G4TruthInfo" << std::endl;
247 while ((thisNode = iter()))
249 if (thisNode->getType() ==
"PHCompositeNode")
251 SearchNode(static_cast<PHCompositeNode*>(thisNode));
253 else if (thisNode->getType() ==
"PHIODataNode")
255 if (thisNode->getName().find(
"G4HIT_") == 0)
281 iter != range.second;
286 std::set<int> remove_ids;
291 int g4particle_id = *jter;
295 remove_ids.insert(g4particle_id);
300 for (std::set<int>::iterator jter = remove_ids.begin();
301 jter != remove_ids.end();
307 std::set<int> remove_more_ids;
308 for (std::map<
int, std::set<PHG4HitDefs::keytype> >::iterator jter = shower->
begin_g4hit_id();
312 int g4hitmap_id = jter->first;
313 std::map<int, PHG4HitContainer*>::iterator mapiter =
m_HitContainerMap.find(g4hitmap_id);
320 for (std::set<PHG4HitDefs::keytype>::iterator kter = jter->second.begin();
321 kter != jter->second.end();)
325 PHG4Hit* g4hit = mapiter->second->findHit(g4hit_id);
329 jter->second.erase(kter++);
338 if (jter->second.empty())
340 remove_more_ids.insert(g4hitmap_id);
344 for (std::set<int>::iterator jter = remove_more_ids.begin();
345 jter != remove_more_ids.end();
354 iter != range.second;)
375 shwiter != range.second;
381 std::vector<std::vector<float> > points;
382 std::vector<float> weights;
386 for (std::map<
int, std::set<PHG4HitDefs::keytype> >::iterator iter = shower->
begin_g4hit_id();
390 int g4hitmap_id = iter->first;
391 std::map<int, PHG4HitContainer*>::iterator mapiter =
m_HitContainerMap.find(g4hitmap_id);
399 unsigned int nhits = 0;
402 float light_yield = 0.0;
407 for (std::set<PHG4HitDefs::keytype>::iterator kter = iter->second.begin();
408 kter != iter->second.end();
416 cout <<
PHWHERE <<
" missing g4hit" << endl;
423 cout <<
PHWHERE <<
" missing g4particle for track "
431 cout <<
PHWHERE <<
" missing g4vertex" << endl;
437 if (!isnan(g4hit->
get_x(0)) &&
438 !isnan(g4hit->
get_y(0)) &&
439 !isnan(g4hit->
get_z(0)))
441 std::vector<float> entry(3);
442 entry[0] = g4hit->
get_x(0);
443 entry[1] = g4hit->
get_y(0);
444 entry[2] = g4hit->
get_z(0);
446 points.push_back(entry);
448 weights.push_back(w);
453 if (!isnan(g4hit->
get_x(1)) &&
454 !isnan(g4hit->
get_y(1)) &&
455 !isnan(g4hit->
get_z(1)))
457 std::vector<float> entry(3);
458 entry[0] = g4hit->
get_x(1);
459 entry[1] = g4hit->
get_y(1);
460 entry[2] = g4hit->
get_z(1);
462 points.push_back(entry);
464 weights.push_back(w);
493 if (nhits) shower->
set_nhits(g4hitmap_id, nhits);
494 if (edep != 0.0) shower->
set_edep(g4hitmap_id, edep);
495 if (eion != 0.0) shower->
set_eion(g4hitmap_id, eion);
496 if (light_yield != 0.0) shower->
set_light_yield(g4hitmap_id, light_yield);
497 if (edep_h != 0.0) shower->
set_eh_ratio(g4hitmap_id, edep_e / edep_h);
503 Eigen::Matrix<double, Eigen::Dynamic, 3>
X(points.size(), 3);
504 Eigen::Matrix<double, Eigen::Dynamic, 1> W(weights.size(), 1);
506 for (
unsigned int i = 0; i < points.size(); ++i)
508 for (
unsigned int j = 0; j < 3; ++j)
510 X(i, j) = points[i][j];
512 W(i, 0) = weights[i];
516 double prefactor = 1.0 / sumw;
517 Eigen::Matrix<double, 1, 3> mean = prefactor * W.transpose() *
X;
520 for (
unsigned int i = 0; i < points.size(); ++i)
522 for (
unsigned int j = 0; j < 3; ++j)
X(i, j) = points[i][j] - mean(0, j);
526 prefactor = sumw / (sumw*sumw - sumw2);
527 Eigen::Matrix<double, 3, 3> covar = prefactor * (
X.transpose() * W.asDiagonal() *
X);
529 shower->
set_x(mean(0, 0));
530 shower->
set_y(mean(0, 1));
531 shower->
set_z(mean(0, 2));
533 for (
unsigned int i = 0; i < 3; ++i)
535 for (
unsigned int j = 0; j <= i; ++j)