Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TruthClusterizerBase.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TruthClusterizerBase.cc
1 #include "TruthClusterizerBase.h"
2 
3 #include <g4main/PHG4Hit.h>
9 
10 /* #include <trackbase/InttDefs.h> */
11 #include <trackbase/TpcDefs.h>
14 #include <trackbase/TrkrDefs.h>
15 #include <trackbase/TrkrHit.h>
16 #include <trackbase/TrkrHitSet.h>
18 #include <trackbase/TrkrHitv2.h> // for TrkrHit
19 
22 
24 #include <fun4all/SubsysReco.h>
25 
27 #include <mvtx/MvtxHitPruner.h>
28 #include <mvtx/MvtxClusterizer.h>
29 
30 #include <phool/PHCompositeNode.h>
31 #include <phool/PHIODataNode.h>
32 #include <phool/PHNode.h>
33 #include <phool/PHNodeIterator.h>
34 #include <phool/PHObject.h> // for PHObject
35 #include <phool/getClass.h>
36 #include <phool/phool.h>
37 
38 #include <TMatrixFfwd.h> // for TMatrixF
39 #include <TMatrixT.h> // for TMatrixT, operator*
40 #include <TMatrixTUtils.h> // for TMatrixTRow
41 
42 
43 #include <array>
44 #include <cmath>
45 #include <iostream>
46 #include <set>
47 #include <vector>
48 
50  : m_hits { new TrkrHitSetContainerv1 }
51 { }
52 
53 void TruthClusterizerBase::init_clusterizer_base( PHCompositeNode*& _topNode, int _verbosity ) {
54  m_topNode = _topNode;
55  m_verbosity = _verbosity;
57  auto dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
58  assert(dstNode);
59  PHNodeIterator dstiter(dstNode);
60  auto DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
61  if (!DetNode)
62  {
63  DetNode = new PHCompositeNode("TRKR");
64  dstNode->addNode(DetNode);
65  }
66 
67  m_truthtracks = findNode::getClass<TrkrTruthTrackContainer>(m_topNode, "TRKR_TRUTHTRACKCONTAINER");
68  if (!m_truthtracks)
69  {
70  PHNodeIterator dstiter(dstNode);
72  auto newNode = new PHIODataNode<PHObject>(m_truthtracks, "TRKR_TRUTHTRACKCONTAINER", "PHObject");
73  DetNode->addNode(newNode);
74  }
75 
76  m_clusters = findNode::getClass<TrkrClusterContainer>(m_topNode, "TRKR_TRUTHCLUSTERCONTAINER");
77  if (!m_clusters)
78  {
80  auto newNode = new PHIODataNode<PHObject>(m_clusters, "TRKR_TRUTHCLUSTERCONTAINER", "PHObject");
81  DetNode->addNode(newNode);
82  }
83 
84  m_truthinfo = findNode::getClass<PHG4TruthInfoContainer>(m_topNode, "G4TruthInfo");
85  if (!m_truthinfo)
86  {
87  std::cout << PHWHERE << " PHG4TruthInfoContainer node not found on node tree" << std::endl;
89  }
90 }
91 
93  delete m_hits;
94 }
95 
97  int new_trkid = (hit==nullptr) ? -1 : hit->get_trkid();
98  m_is_new_track = (new_trkid != m_trkid);
99  if (m_verbosity>5) std::cout << PHWHERE << std::endl << " -> Checking status of PHG4Hit. Track id("<<new_trkid<<")" << std::endl;
100  if (!m_is_new_track) return;
101 
102  m_trkid = new_trkid; // although m_current_track won't catch up until update_track is called
105 }
106 
107 // to call if m_was_emb=true and after clustering
109  m_hits->Reset(); // clear out the old hits
110  for (auto hitsetkey : pass_clusters->getHitSetKeys()) {
111  m_hitsetkey_cnt.try_emplace(hitsetkey,0);
112  unsigned int& cnt = m_hitsetkey_cnt[hitsetkey];
113  auto range = pass_clusters->getClusters(hitsetkey);
114  for (auto cluster = range.first; cluster != range.second; ++cluster) {
115  auto ckey = TrkrDefs::genClusKey(hitsetkey, cnt);
116  m_clusters->addClusterSpecifyKey(ckey, cluster->second);
118  ++cnt;
119  }
120  }
121  m_was_emb = false;
122 }
123 
124 // if m_is_new_track
126  m_current_track = m_is_emb ? m_truthtracks->getTruthTrack(m_trkid, m_truthinfo) : nullptr;
127 }
128 
132  float neffelectrons)
133 {
134  if (!m_is_emb) return;
136  // See if this hit already exists
137  TrkrHit *hit = nullptr;
138  hit = hitsetit->second->getHit(hitkey);
139  if (!hit)
140  {
141  // create a new one
142  hit = new TrkrHitv2();
143  hitsetit->second->addHitSpecificKey(hitkey, hit);
144  }
145  // Either way, add the energy to it -- adc values will be added at digitization
146  hit->addEnergy(neffelectrons);
147 }
148 
149 void TruthClusterizerBase::print_clusters(int nclusprint) {
150  std::cout << PHWHERE << ": content of clusters " << std::endl;
151  auto& tmap = m_truthtracks->getMap();
152  std::cout << " Number of tracks: " << tmap.size() << std::endl;
153  for (auto& _pair : tmap) {
154  auto& track = _pair.second;
155 
156  printf("id(%2i) phi:eta:pt(", (int)track->getTrackid());
157  std::cout << "phi:eta:pt(";
158  printf("%5.2f:%5.2f:%5.2f", track->getPhi(), track->getPseudoRapidity(), track->getPt());
159  /* Form("%5.2:%5.2:%5.2", track->getPhi(), track->getPseudoRapidity(), track->getPt()) */
160  //<<track->getPhi()<<":"<<track->getPseudoRapidity()<<":"<<track->getPt()
161  std::cout << ") nclusters(" << track->getClusters().size() <<") ";
162  if (m_verbosity <= 10) { std::cout << std::endl; }
163  else {
164  int nclus = 0;
165  for (auto cluskey : track->getClusters()) {
166  std::cout << " "
167  << ((int) TrkrDefs::getHitSetKeyFromClusKey(cluskey)) <<":index(" <<
168  ((int) TrkrDefs::getClusIndex(cluskey)) << ")";
169  ++nclus;
170  if (nclusprint > 0 && nclus >= nclusprint) {
171  std::cout << " ... ";
172  break;
173  }
174  }
175  }
176  }
177  std::cout << PHWHERE << " ----- end of clusters " << std::endl;
178 }
179