Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HepMCCompress.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file HepMCCompress.cc
1 #include "HepMCCompress.h"
2 #include "PHG4InEvent.h"
3 #include "PHG4Particle.h"
4 
8 
9 #include <half/half.h>
10 
12 
13 #include <phool/PHCompositeNode.h>
14 #include <phool/getClass.h>
15 
16 #include <HepMC/GenEvent.h>
17 
18 #include <gsl/gsl_const.h>
19 
20 #include <list>
21 
22 using namespace std;
23 const double mm_over_c_to_sec = 0.1 / GSL_CONST_CGS_SPEED_OF_LIGHT; // pythia vtx time seems to be in mm/c
25 
28 {
29 public:
31  bool operator()( const HepMC::GenParticle* p )
32  {
33  if ( !p->end_vertex() && p->status() == 1 ) return 1;
34  return 0;
35  }
36 };
37 
39 
41  SubsysReco(name)
42 {}
43 
44 int
46 {
47  VariableArrayContainer *vararraycontainer = findNode::getClass<VariableArrayContainer>(topNode, "HEPMC_VarArray");
48  if (!vararraycontainer)
49  {
50  PHNodeIterator iter( topNode );
51  PHCompositeNode *dstNode;
52  dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST" ));
53 
54  vararraycontainer = new VariableArrayContainer();
56  vararraycontainer->AddVarArray(vararray);
57  vararray = new VariableArray(varids::G4PARTICLEV1);
58  vararraycontainer->AddVarArray(vararray);
59 
60  PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(vararraycontainer, "HEPMC_VarArray", "PHObject");
61  dstNode->addNode(newNode);
62  }
63  return 0;
64 }
65 
66 int
68 {
69  HepMC::GenEvent *evt = findNode::getClass<HepMC::GenEvent>(topNode, "HEPMC");
70  if (!evt)
71  {
72  cout << PHWHERE << " no evt pointer under HEPMC Node found" << endl;
74  }
75  VariableArrayContainer *vararraycontainer = findNode::getClass<VariableArrayContainer>(topNode, "HEPMC_VarArray");
76  if (!vararraycontainer)
77  {
78  cout << PHWHERE << "no PHG4INEVENT node" << endl;
80  }
81 
82  vector<short> shepmcvtxvec;
83 
84  std::list<HepMC::GenParticle*> finalstateparticles;
85  // units in G4 interface are GeV and CM
86  // const double mom_factor = HepMC::Units::conversion_factor( evt->momentum_unit(), HepMC::Units::GEV );
87  const double length_factor = HepMC::Units::conversion_factor( evt->length_unit(), HepMC::Units::CM );
88  for ( HepMC::GenEvent::vertex_iterator v = evt->vertices_begin();
89  v != evt->vertices_end(); ++v )
90  {
91  finalstateparticles.clear();
92  for (HepMC::GenVertex::particle_iterator p = (*v)->particles_begin(HepMC::children); p != (*v)->particles_end(HepMC::children); ++p)
93  {
94  if (isfinal(*p))
95  {
96  if (!select_pid.empty())
97  {
98  if (select_pid.find((*p)->pdg_id()) != select_pid.end())
99  {
100  finalstateparticles.push_back(*p);
101  }
102  continue;
103  }
104  if (!exclude_pid.empty())
105  {
106  if (exclude_pid.find((*p)->pdg_id()) != exclude_pid.end())
107  {
108  continue;
109  }
110  }
111  finalstateparticles.push_back(*p);
112  }
113  }
114  if (!finalstateparticles.empty())
115  {
116  // cout << "Vertex : " << endl;
117  // (*v)->print();
118  // cout << "id: " << (*v)->barcode() << endl;
119  // cout << "x: " << (*v)->position().x() << endl;
120  // cout << "y: " << (*v)->position().y() << endl;
121  // cout << "z: " << (*v)->position().z() << endl;
122  // cout << "t: " << (*v)->position().t() << endl;
123  // cout << "Particles" << endl;
124  shepmcvtxvec.push_back((*v)->barcode());
125  shepmcvtxvec.push_back(FloatToInt((*v)->position().x()*length_factor));
126  shepmcvtxvec.push_back(FloatToInt((*v)->position().y()*length_factor));
127  shepmcvtxvec.push_back(FloatToInt((*v)->position().z()*length_factor));
128 
129  // for (fiter = finalstateparticles.begin(); fiter != finalstateparticles.end(); fiter++)
130  // {
131  // // (*fiter)->print();
132  // PHG4Particle *particle = new PHG4Particle();
133  // particle->set_pid((*fiter)->pdg_id());
134  // particle->set_px((*fiter)->momentum().px()*mom_factor);
135  // particle->set_py((*fiter)->momentum().py()*mom_factor);
136  // particle->set_pz((*fiter)->momentum().pz()*mom_factor);
137  // ineve->AddParticle((*v)->barcode(), particle);
138  // }
139  // }
140  }
141  }
142  // ineve->identify();
144 }
145 
146 short int
147 HepMCCompress::FloatToInt(const float rval) const
148 {
149  half ftoi(rval);
150  return ftoi.bits();
151 }