19 #include <Geant4/G4Event.hh>
20 #include <Geant4/G4PrimaryParticle.hh>
21 #include <Geant4/G4PrimaryVertex.hh>
22 #include <Geant4/G4VUserPrimaryParticleInformation.hh>
25 #pragma GCC diagnostic push
26 #pragma GCC diagnostic ignored "-Wshadow"
27 #include <Eigen/Dense>
28 #pragma GCC diagnostic pop
45 : m_TruthInfoContainer(nullptr)
46 , m_LowerKeyPrevExist(0)
47 , m_UpperKeyPrevExist(0)
57 std::cout <<
"PHG4TruthEventAction::EndOfEventAction - unable to find G4TruthInfo node" << std::endl;
75 std::cout <<
"PHG4TruthEventAction::EndOfEventAction - unable to find G4TruthInfo node" << std::endl;
86 std::set<int> savelist;
87 std::set<int> savevtxlist;
89 for (std::set<int>::const_iterator write_iter =
m_WriteSet.begin();
93 std::vector<int> wrttracks;
94 std::vector<int> wrtvtx;
97 int mytrkid = *write_iter;
101 if (savelist.find(mytrkid) != savelist.end())
107 wrttracks.push_back(mytrkid);
114 while (savelist.find(parentid) == savelist.end() && parentid != 0)
117 wrttracks.push_back(parentid);
124 while (wrttracks.begin() != wrttracks.end())
126 savelist.insert(wrttracks.back());
127 wrttracks.pop_back();
130 while (wrtvtx.begin() != wrtvtx.end())
132 savevtxlist.insert(wrtvtx.back());
143 int removed[4] = {0};
146 while (truthiter != truth_range.second)
149 int trackid = truthiter->first;
150 if (savelist.find(trackid) == savelist.end())
169 savevtxlist.insert((truthiter->second)->get_vtx_id());
182 while (vtxiter != vtxrange.second)
185 if (savevtxlist.find(vtxiter->first) == savevtxlist.end())
198 G4PrimaryVertex* pvtx = evt->GetPrimaryVertex();
201 G4PrimaryParticle* part = pvtx->GetPrimary();
213 part = part->GetNext();
215 pvtx = pvtx->GetNext();
236 std::cout <<
"PHG4TruthEventAction::SetInterfacePointers - unable to find G4TruthInfo" << std::endl;
251 while ((thisNode = iter()))
253 if (thisNode->getType() ==
"PHCompositeNode")
255 SearchNode(static_cast<PHCompositeNode*>(thisNode));
257 else if (thisNode->getType() ==
"PHIODataNode")
259 if (thisNode->getName().find(
"G4HIT_") == 0)
285 iter != range.second;
290 std::set<int> remove_ids;
295 int g4particle_id = *jter;
299 remove_ids.insert(g4particle_id);
304 for (std::set<int>::iterator jter = remove_ids.begin();
305 jter != remove_ids.end();
311 std::set<int> remove_more_ids;
312 for (std::map<
int, std::set<PHG4HitDefs::keytype> >::iterator jter = shower->
begin_g4hit_id();
316 int g4hitmap_id = jter->first;
317 std::map<int, PHG4HitContainer*>::iterator mapiter =
m_HitContainerMap.find(g4hitmap_id);
324 for (std::set<PHG4HitDefs::keytype>::iterator kter = jter->second.begin();
325 kter != jter->second.end();)
329 PHG4Hit* g4hit = mapiter->second->findHit(g4hit_id);
333 jter->second.erase(kter++);
342 if (jter->second.empty())
344 remove_more_ids.insert(g4hitmap_id);
348 for (std::set<int>::iterator jter = remove_more_ids.begin();
349 jter != remove_more_ids.end();
358 iter != range.second;)
379 shwiter != range.second;
385 std::vector<std::vector<float> > points;
386 std::vector<float> weights;
390 for (std::map<
int, std::set<PHG4HitDefs::keytype> >::iterator iter = shower->
begin_g4hit_id();
394 int g4hitmap_id = iter->first;
395 std::map<int, PHG4HitContainer*>::iterator mapiter =
m_HitContainerMap.find(g4hitmap_id);
403 unsigned int nhits = 0;
406 float light_yield = 0.0;
411 for (std::set<PHG4HitDefs::keytype>::iterator kter = iter->second.begin();
412 kter != iter->second.end();
420 cout <<
PHWHERE <<
" missing g4hit" << endl;
427 cout <<
PHWHERE <<
" missing g4particle for track "
435 cout <<
PHWHERE <<
" missing g4vertex" << endl;
441 if (!isnan(g4hit->
get_x(0)) &&
442 !isnan(g4hit->
get_y(0)) &&
443 !isnan(g4hit->
get_z(0)))
445 std::vector<float>
entry(3);
446 entry[0] = g4hit->
get_x(0);
447 entry[1] = g4hit->
get_y(0);
448 entry[2] = g4hit->
get_z(0);
450 points.push_back(entry);
452 weights.push_back(w);
457 if (!isnan(g4hit->
get_x(1)) &&
458 !isnan(g4hit->
get_y(1)) &&
459 !isnan(g4hit->
get_z(1)))
461 std::vector<float>
entry(3);
462 entry[0] = g4hit->
get_x(1);
463 entry[1] = g4hit->
get_y(1);
464 entry[2] = g4hit->
get_z(1);
466 points.push_back(entry);
468 weights.push_back(w);
477 if (abs(particle->
get_pid()) == 11)
497 if (nhits) shower->
set_nhits(g4hitmap_id, nhits);
498 if (edep != 0.0) shower->
set_edep(g4hitmap_id, edep);
499 if (eion != 0.0) shower->
set_eion(g4hitmap_id, eion);
500 if (light_yield != 0.0) shower->
set_light_yield(g4hitmap_id, light_yield);
501 if (edep_h != 0.0) shower->
set_eh_ratio(g4hitmap_id, edep_e / edep_h);
507 Eigen::Matrix<double, Eigen::Dynamic, 3>
X(points.size(), 3);
508 Eigen::Matrix<double, Eigen::Dynamic, 1> W(weights.size(), 1);
510 for (
unsigned int i = 0;
i < points.size(); ++
i)
512 for (
unsigned int j = 0;
j < 3; ++
j)
514 X(
i,
j) = points[
i][
j];
516 W(
i, 0) = weights[
i];
520 double prefactor = 1.0 / sumw;
521 Eigen::Matrix<double, 1, 3>
mean = prefactor * W.transpose() *
X;
524 for (
unsigned int i = 0;
i < points.size(); ++
i)
526 for (
unsigned int j = 0;
j < 3; ++
j)
X(
i,
j) = points[
i][
j] -
mean(0,
j);
530 prefactor = sumw / (sumw*sumw - sumw2);
531 Eigen::Matrix<double, 3, 3> covar = prefactor * (
X.transpose() * W.asDiagonal() *
X);
537 for (
unsigned int i = 0;
i < 3; ++
i)
539 for (
unsigned int j = 0;
j <=
i; ++
j)