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