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