ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GlobalVertexReco.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file GlobalVertexReco.cc
1 #include "GlobalVertexReco.h"
2 
3 #include "GlobalVertexMapv1.h"
4 #include "GlobalVertexMap.h" // for GlobalVertexMap
5 #include "GlobalVertex.h" // for GlobalVertex, GlobalVe...
6 #include "GlobalVertexv1.h"
7 
8 #include <g4bbc/BbcVertex.h>
9 #include <g4bbc/BbcVertexMap.h>
10 
13 
15 #include <fun4all/SubsysReco.h> // for SubsysReco
16 
17 #include <phool/PHCompositeNode.h>
18 #include <phool/PHIODataNode.h>
19 #include <phool/PHNode.h> // for PHNode
20 #include <phool/PHNodeIterator.h>
21 #include <phool/PHObject.h> // for PHObject
22 #include <phool/getClass.h>
23 #include <phool/phool.h> // for PHWHERE
24 
25 #include <cfloat>
26 #include <cmath>
27 #include <cstdlib> // for exit
28 #include <set> // for _Rb_tree_const_iterator
29 #include <iostream>
30 #include <utility> // for pair
31 
32 using namespace std;
33 
35  : SubsysReco(name)
36  , _xdefault(0.0)
37  , _xerr(0.3)
38  , _ydefault(0.0)
39  , _yerr(0.3)
40  , _tdefault(0.0)
41  , _terr(0.2)
42 {
43 }
44 
46 
48 {
50 }
51 
53 {
54  if (Verbosity() > 0)
55  {
56  cout << "======================= GlobalVertexReco::InitRun() =======================" << endl;
57  cout << "===========================================================================" << endl;
58  }
59 
60  return CreateNodes(topNode);
61 }
62 
64 {
65  if (Verbosity() > 1) cout << "GlobalVertexReco::process_event -- entered" << endl;
66 
67  //---------------------------------
68  // Get Objects off of the Node Tree
69  //---------------------------------
70  GlobalVertexMap *globalmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
71  if (!globalmap)
72  {
73  cout << PHWHERE << "::ERROR - cannot find GlobalVertexMap" << endl;
74  exit(-1);
75  }
76 
77  SvtxVertexMap *svtxmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
78  BbcVertexMap *bbcmap = findNode::getClass<BbcVertexMap>(topNode, "BbcVertexMap");
79 
80  // we will make 3 different kinds of global vertexes
81  // (1) SVTX+BBC vertexes - we match SVTX vertex to the nearest BBC vertex within 3 sigma in zvertex
82  // the spatial point comes from the SVTX, the timing from the BBC
83  // number of SVTX+BBC vertexes <= number of SVTX vertexes
84 
85  // (2) SVTX only vertexes - for those cases where the BBC wasn't simulated,
86  // or all the BBC vertexes are outside the 3 sigma matching requirement
87  // we pass forward the 3d point from the SVTX with the default timing info
88 
89  // (3) BBC only vertexes - use the default x,y positions on this module and
90  // pull in the bbc z and bbc t
91 
92  // there may be some quirks as we get to large luminosity and the BBC becomes
93  // untrust worthy, I'm guessing analyzers would resort exclusively to (1) or (2)
94  // in those cases
95 
96  std::set<unsigned int> used_svtx_vtxids;
97  std::set<unsigned int> used_bbc_vtxids;
98 
99  if (svtxmap && bbcmap)
100  {
101  if (Verbosity()) cout << "GlobalVertexReco::process_event - svtxmap && bbcmap" << endl;
102 
103  for (SvtxVertexMap::ConstIter svtxiter = svtxmap->begin();
104  svtxiter != svtxmap->end();
105  ++svtxiter)
106  {
107  const SvtxVertex *svtx = svtxiter->second;
108 
109  const BbcVertex *bbc_best = nullptr;
110  float min_sigma = FLT_MAX;
111  for (BbcVertexMap::ConstIter bbciter = bbcmap->begin();
112  bbciter != bbcmap->end();
113  ++bbciter)
114  {
115  const BbcVertex *bbc = bbciter->second;
116 
117  float combined_error = sqrt(svtx->get_error(2, 2) + pow(bbc->get_z_err(), 2));
118  float sigma = fabs(svtx->get_z() - bbc->get_z()) / combined_error;
119  if (sigma < min_sigma)
120  {
121  min_sigma = sigma;
122  bbc_best = bbc;
123  }
124  }
125 
126  if (min_sigma > 3.0 || !bbc_best) continue;
127 
128  // we have a matching pair
129  GlobalVertex *vertex = new GlobalVertexv1();
130 
131  for (unsigned int i = 0; i < 3; ++i)
132  {
133  vertex->set_position(i, svtx->get_position(i));
134  for (unsigned int j = i; j < 3; ++j)
135  {
136  vertex->set_error(i, j, svtx->get_error(i, j));
137  }
138  }
139 
140  vertex->set_t(bbc_best->get_t());
141  vertex->set_t_err(bbc_best->get_t_err());
142 
143  vertex->set_chisq(svtx->get_chisq());
144  vertex->set_ndof(svtx->get_ndof());
145 
146  vertex->insert_vtxids(GlobalVertex::SVTX, svtx->get_id());
147  used_svtx_vtxids.insert(svtx->get_id());
148  vertex->insert_vtxids(GlobalVertex::BBC, bbc_best->get_id());
149  used_bbc_vtxids.insert(bbc_best->get_id());
150 
151  globalmap->insert(vertex);
152 
153  if (Verbosity()) vertex->identify();
154  }
155  }
156 
157  // okay now loop over all unused SVTX vertexes (2nd class)...
158  if (svtxmap)
159  {
160  if (Verbosity()) cout << "GlobalVertexReco::process_event - svtxmap " << endl;
161 
162  for (SvtxVertexMap::ConstIter svtxiter = svtxmap->begin();
163  svtxiter != svtxmap->end();
164  ++svtxiter)
165  {
166  const SvtxVertex *svtx = svtxiter->second;
167 
168  if (used_svtx_vtxids.find(svtx->get_id()) != used_svtx_vtxids.end()) continue;
169  if (isnan(svtx->get_z())) continue;
170 
171  // we have a standalone SVTX vertex
172  GlobalVertex *vertex = new GlobalVertexv1();
173 
174  for (unsigned int i = 0; i < 3; ++i)
175  {
176  vertex->set_position(i, svtx->get_position(i));
177  for (unsigned int j = i; j < 3; ++j)
178  {
179  vertex->set_error(i, j, svtx->get_error(i, j));
180  }
181  }
182 
183  // default time could also come from somewhere else at some point
184  vertex->set_t(_tdefault);
185  vertex->set_t_err(_terr);
186 
187  vertex->set_chisq(svtx->get_chisq());
188  vertex->set_ndof(svtx->get_ndof());
189 
190  vertex->insert_vtxids(GlobalVertex::SVTX, svtx->get_id());
191  used_svtx_vtxids.insert(svtx->get_id());
192 
193  globalmap->insert(vertex);
194 
195  if (Verbosity()) vertex->identify();
196  }
197  }
198 
199  // okay now loop over all unused BBC vertexes (3rd class)...
200  if (bbcmap)
201  {
202  if (Verbosity()) cout << "GlobalVertexReco::process_event - bbcmap" << endl;
203 
204  for (BbcVertexMap::ConstIter bbciter = bbcmap->begin();
205  bbciter != bbcmap->end();
206  ++bbciter)
207  {
208  const BbcVertex *bbc = bbciter->second;
209 
210  if (used_bbc_vtxids.find(bbc->get_id()) != used_bbc_vtxids.end()) continue;
211  if (isnan(bbc->get_z())) continue;
212 
213  GlobalVertex *vertex = new GlobalVertexv1();
214 
215  // nominal beam location
216  // could be replaced with a beam spot some day
217  vertex->set_x(_xdefault);
218  vertex->set_y(_ydefault);
219  vertex->set_z(bbc->get_z());
220 
221  vertex->set_t(bbc->get_t());
222  vertex->set_t_err(bbc->get_t_err());
223 
224  vertex->set_error(0, 0, pow(_xerr, 2));
225  vertex->set_error(0, 1, 0.0);
226  vertex->set_error(0, 2, 0.0);
227 
228  vertex->set_error(1, 1, 0.0);
229  vertex->set_error(1, 1, pow(_yerr, 2));
230  vertex->set_error(1, 1, 0.0);
231 
232  vertex->set_error(0, 2, 0.0);
233  vertex->set_error(1, 2, 0.0);
234  vertex->set_error(2, 2, pow(bbc->get_z_err(), 2));
235 
236  vertex->insert_vtxids(GlobalVertex::BBC, bbc->get_id());
237  used_bbc_vtxids.insert(bbc->get_id());
238 
239  globalmap->insert(vertex);
240 
241  if (Verbosity()) vertex->identify();
242  }
243  }
244 
246 }
247 
249 {
251 }
252 
254 {
255  PHNodeIterator iter(topNode);
256 
257  // Looking for the DST node
258  PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
259  if (!dstNode)
260  {
261  cout << PHWHERE << "DST Node missing, doing nothing." << endl;
263  }
264 
265  // store the GLOBAL stuff under a sub-node directory
266  PHCompositeNode *globalNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "GLOBAL"));
267  if (!globalNode)
268  {
269  globalNode = new PHCompositeNode("GLOBAL");
270  dstNode->addNode(globalNode);
271  }
272 
273  // create the GlobalVertexMap
274  GlobalVertexMap *vertexes = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
275  if (!vertexes)
276  {
277  vertexes = new GlobalVertexMapv1();
278  PHIODataNode<PHObject> *VertexMapNode = new PHIODataNode<PHObject>(vertexes, "GlobalVertexMap", "PHObject");
279  globalNode->addNode(VertexMapNode);
280  }
281  else
282  {
283  cout << PHWHERE << "::ERROR - GlobalVertexMap pre-exists, but should not, perhaps you are also running the FastSim" << endl;
284  exit(-1);
285  }
286 
288 }