ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GlobalVertexFastSimReco.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file GlobalVertexFastSimReco.cc
2 
3 #include "GlobalVertexMap.h" // for GlobalVertexMap
4 #include "GlobalVertexMapv1.h"
5 #include "GlobalVertex.h" // for GlobalVertex
6 #include "GlobalVertexv1.h"
7 
9 #include <g4main/PHG4VtxPoint.h>
10 
12 #include <fun4all/SubsysReco.h> // for SubsysReco
13 
14 #include <phool/PHCompositeNode.h>
15 #include <phool/PHIODataNode.h>
16 #include <phool/PHNode.h> // for PHNode
17 #include <phool/PHNodeIterator.h>
18 #include <phool/PHObject.h> // for PHObject
19 #include <phool/PHRandomSeed.h>
20 #include <phool/getClass.h>
21 #include <phool/phool.h> // for PHWHERE
22 
23 #include <gsl/gsl_randist.h>
24 #include <gsl/gsl_rng.h> // for gsl_rng_alloc, gsl_rng_free
25 
26 #include <cmath>
27 #include <cstdlib> // for exit
28 #include <iostream>
29 
30 using namespace std;
31 
33  : SubsysReco(name)
34  , _x_smear(NAN)
35  , _y_smear(NAN)
36  , _z_smear(NAN)
37  , _t_smear(NAN)
38 {
39  RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
40 }
41 
43 {
44  gsl_rng_free(RandomGenerator);
45 }
46 
48 {
50 }
51 
53 {
54  if (isnan(_x_smear) ||
55  isnan(_y_smear) ||
56  isnan(_z_smear) ||
57  isnan(_t_smear))
58  {
59  cout << PHWHERE << "::ERROR - smearing must be defined for (x,y,z,t) via set_?_smearing(float)" << endl;
60  exit(-1);
61  }
62 
63  unsigned int seed = PHRandomSeed(); // fixed seed handled in PHRandomSeed()
64  gsl_rng_set(RandomGenerator, seed);
65 
66  if (Verbosity() > 0)
67  {
68  cout << "=================== GlobalVertexFastSimReco::InitRun() ====================" << endl;
69  cout << " x smearing: " << _x_smear << " cm " << endl;
70  cout << " y smearing: " << _y_smear << " cm " << endl;
71  cout << " z smearing: " << _z_smear << " cm " << endl;
72  cout << " t smearing: " << _t_smear << " cm " << endl;
73  cout << " random seed: " << seed << endl;
74  cout << "===========================================================================" << endl;
75  }
76 
77  return CreateNodes(topNode);
78 }
79 
81 {
82  if (Verbosity() > 1) cout << "GlobalVertexFastSimReco::process_event -- entered" << endl;
83 
84  //---------------------------------
85  // Get Objects off of the Node Tree
86  //---------------------------------
87  PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
88  if (!truthinfo)
89  {
90  cout << PHWHERE << "::ERROR - cannot find G4TruthInfo" << endl;
91  exit(-1);
92  }
93 
94  GlobalVertexMap *vertexes = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
95  if (!vertexes)
96  {
97  cout << PHWHERE << "::ERROR - cannot find GlobalVertexMap" << endl;
98  exit(-1);
99  }
100 
101  //---------------------
102  // Smear Truth Vertexes (only one per crossing right now)
103  //---------------------
104 
105  PHG4VtxPoint *point = truthinfo->GetPrimaryVtx(truthinfo->GetPrimaryVertexIndex());
106 
107  GlobalVertex *vertex = new GlobalVertexv1();
108 
109  vertex->set_x(point->get_x() + gsl_ran_gaussian(RandomGenerator, _x_smear));
110  vertex->set_y(point->get_y() + gsl_ran_gaussian(RandomGenerator, _y_smear));
111  vertex->set_z(point->get_z() + gsl_ran_gaussian(RandomGenerator, _z_smear));
112 
113  vertex->set_error(0, 0, _x_smear * _x_smear);
114  vertex->set_error(0, 1, 0.0);
115  vertex->set_error(0, 2, 0.0);
116 
117  vertex->set_error(1, 0, 0.0);
118  vertex->set_error(1, 1, _y_smear * _y_smear);
119  vertex->set_error(1, 2, 0.0);
120 
121  vertex->set_error(2, 0, 0.0);
122  vertex->set_error(2, 1, 0.0);
123  vertex->set_error(2, 2, _z_smear * _z_smear);
124 
125  vertex->set_t(point->get_t() + gsl_ran_gaussian(RandomGenerator, _t_smear));
126  vertex->set_t_err(_t_smear);
127 
128  vertexes->insert(vertex);
129 
131 }
132 
134 {
136 }
137 
139 {
140  PHNodeIterator iter(topNode);
141 
142  // Looking for the DST node
143  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
144  if (!dstNode)
145  {
146  cout << PHWHERE << "DST Node missing, doing nothing." << endl;
148  }
149 
150  // store the GLOBAL stuff under a sub-node directory
151  PHCompositeNode *globalNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "GLOBAL"));
152  if (!globalNode)
153  {
154  globalNode = new PHCompositeNode("GLOBAL");
155  dstNode->addNode(globalNode);
156  }
157 
158  // create the GlobalVertexMap
159  GlobalVertexMap *vertexes = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
160  if (!vertexes)
161  {
162  vertexes = new GlobalVertexMapv1();
163  PHIODataNode<PHObject> *VertexMapNode = new PHIODataNode<PHObject>(vertexes, "GlobalVertexMap", "PHObject");
164  globalNode->addNode(VertexMapNode);
165  }
166  else
167  {
168  cout << PHWHERE << "::ERROR - GlobalVertexMap pre-exists, but should not if running FastSim" << endl;
169  exit(-1);
170  }
171 
173 }