EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
sHEPGen.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file sHEPGen.C
1 #include "sHEPGen.h"
2 
3 #include <phhepmc/PHHepMCGenEvent.h>
4 
5 #include <fun4all/Fun4AllReturnCodes.h>
6 #include <phool/PHCompositeNode.h>
7 #include <phool/PHIODataNode.h>
8 #include <phool/PHNodeIterator.h>
9 #include <phool/PHRandomSeed.h>
10 
11 /* HEPGen includes */
12 #include <hgenmanager.h>
13 
14 /* HepMC includes */
15 #include <HepMC/GenEvent.h>
16 
17 using namespace std;
18 
19 typedef PHIODataNode<PHObject> PHObjectNode_t;
20 
21 sHEPGen::sHEPGen(const std::string &name)
22  : SubsysReco(name)
23  , _eventcount(0)
24  , _p_electron_lab(-20)
25  , _p_hadron_lab(250)
26  , _p4_electron_lab(nullptr)
27  , _p4_hadron_lab(nullptr)
28  , _p4_hadron_lab_invert(nullptr)
29  , _p4_electron_prest(nullptr)
30  , _p4_hadron_prest(nullptr)
31  , _hgenManager(NULL)
32  , _datacardFile("hepgen_dvcs.data")
33 {
34  PHHepMCGenHelper::set_embedding_id(1); // default embedding ID to 1
35 }
36 
37 sHEPGen::~sHEPGen()
38 {
39 }
40 
41 int sHEPGen::Init(PHCompositeNode *topNode)
42 {
43  printlogo();
44 
45  /* electron and proton mass */
46  double mass_e = 5.109989e-4;
47  double mass_p = 9.382720e-1;
48 
49  /* 4-Vectors of colliding electron and proton in laboratory frame */
50  _p4_electron_lab = new HLorentzVector(0.,
51  0.,
52  _p_electron_lab,
53  sqrt(_p_electron_lab * _p_electron_lab + mass_e * mass_e));
54 
55  _p4_hadron_lab = new HLorentzVector(0.,
56  0.,
57  _p_hadron_lab,
58  sqrt(_p_hadron_lab * _p_hadron_lab + mass_p * mass_p));
59 
60  /* The current version of HEPGen supports only a "Fixed Target" mode, i.e.
61  the target (proton) is assumed at rest. Therefore, need to boost collision
62  from the laboratory frame to the proton-at-rest frame. */
63  _p4_hadron_lab_invert = new HLorentzVector(0.,
64  0.,
65  -1 * _p_hadron_lab,
66  sqrt(_p_hadron_lab * _p_hadron_lab + mass_p * mass_p));
67 
68  /* 4-Vectors of colliding electron and proton in proton-at-rest frame */
69  _p4_electron_prest = new HLorentzVector(*_p4_electron_lab);
70  _p4_electron_prest->boost(sqrt(_p4_hadron_lab_invert->getQuare()), *_p4_hadron_lab_invert);
71 
72  _p4_hadron_prest = new HLorentzVector(*_p4_hadron_lab);
73  _p4_hadron_prest->boost(sqrt(_p4_hadron_lab_invert->getQuare()), *_p4_hadron_lab_invert);
74 
75  if (verbosity > 2)
76  {
77  cout << "Electron and proton in laboratory frame:" << endl;
78  _p4_electron_lab->print();
79  _p4_hadron_lab->print();
80 
81  cout << "Electron and proton in proton-at-rest frame:" << endl;
82  _p4_electron_prest->print();
83  _p4_hadron_prest->print();
84  }
85 
86  /* flip sign of electron momentum z-component for HEPGen */
87  HVector3 flip_pz(_p4_electron_prest->getVector());
88  flip_pz.setZ(flip_pz.Z() * -1);
89  _p4_electron_prest->setVector(flip_pz);
90 
91  /* get instance of HepGenManager */
92  _hgenManager = HGenManager::getInstance();
93  _hgenManager->loadSettings(_datacardFile);
94  _hgenManager->setupGenerator();
95 
96  /* set beam parameters */
97  cout << "Colliding " << _p4_electron_lab->getVector().Z() << " GeV electron with " << _p4_hadron_lab->getVector().Z() << " GeV proton (laboratory frame)" << endl;
98  cout << "----ELEPT (proton-at-rest): " << _p4_electron_prest->getVector().Z() << " GeV" << endl;
99  _hgenManager->getParamManager()->getStruct()->ELEPT = _p4_electron_prest->getVector().Z();
100  _hgenManager->getParamManager()->getStruct()->PARL.at(2) = _p4_electron_prest->getVector().Z();
101 
102  /* set random seed */
103  unsigned int seed = PHRandomSeed();
104  _hgenManager->setSeed(seed);
105 
106  /* enable detailed event record printput for debugging */
107  if (verbosity > 2)
108  _hgenManager->enableDebug();
109 
110  create_node_tree(topNode);
111 
112  _eventcount = 0;
113 
114  return Fun4AllReturnCodes::EVENT_OK;
115 }
116 
117 int sHEPGen::End(PHCompositeNode *topNode)
118 {
119  cout << "Reached the sHEPGen::End()" << endl;
120 
121  //-* dump out closing info (cross-sections, etc)
122 
123  if (verbosity > 1) cout << "sHEPGen::End - I'm here!" << endl;
124 
125  return Fun4AllReturnCodes::EVENT_OK;
126 }
127 
128 int sHEPGen::process_event(PHCompositeNode *topNode)
129 {
130  if (verbosity > 1) cout << "sHEPGen::process_event - event: " << _eventcount << endl;
131 
132  _hgenManager->oneShot();
133 
134  if (verbosity > 2)
135  _hgenManager->getEvent()->printDebug();
136 
137  HEvent *evt_mc = _hgenManager->getEvent();
138 
139  /* Create HepMC GenEvent */
140  HepMC::GenEvent *evt = new HepMC::GenEvent();
141 
142  /* define the units (Pythia uses GeV and mm) */
143  evt->use_units(HepMC::Units::GEV, HepMC::Units::MM);
144 
145  /* add global information to the event */
146  evt->set_event_number(_eventcount);
147 
148  /* Set the PDF information */
149  HepMC::PdfInfo pdfinfo;
150  pdfinfo.set_scalePDF(evt_mc->getQsq());
151  pdfinfo.set_x2(evt_mc->getXbj());
152  evt->set_pdf_info(pdfinfo);
153 
154  /* Set additional event parameters */
155  evt->set_event_scale(evt_mc->getQsq());
156 
157  /* Event kinematics */
158  /* @TODO How can we store this information in HepMC record? */
159  //evt_mc->getNu();
160  //evt_mc->getY();
161  //evt_mc->getWsq();
162  //evt_mc->getXbj();
163  //evt_mc->getT();
164  //evt_mc->getQsq();
165  //evt_mc->getS();
166  //evt_mc->getEmiss();
167 
168  /* Create single HepMC vertex for event */
169  /* @TODO: Implement multiple vertices e.g. for decay particles */
170  HepMC::GenVertex *hepmcvtx = new HepMC::GenVertex(HepMC::FourVector(0,
171  0,
172  0,
173  0));
174 
175  /* Create HepMC particle records */
176  HEventData *edata = evt_mc->getStruct();
177  for (unsigned p = 0; p < edata->listOfParticles.size(); p++)
178  {
179  if (verbosity > 4)
180  {
181  cout << "______new particle_______" << endl;
182  cout << "Index: " << p + 1 << endl;
183  cout << "PID: " << edata->listOfParticles.at(p)->getParticleType() << " -- " << (edata->listOfParticles.at(p) == &(edata->incBeamParticle)) << endl;
184  cout << "Particle aux flag: " << edata->listOfParticles.at(p)->getParticleAuxFlag() << endl;
185  cout << "Particle origin: " << edata->listOfParticles.at(p)->getParticleOrigin() << endl;
186  cout << "Particle daughter1: " << edata->listOfParticles.at(p)->getParticleDaughter1() << endl;
187  cout << "Particle daughter2: " << edata->listOfParticles.at(p)->getParticleDaughter2() << endl;
188  }
189 
190  HLorentzVector v4_particle_p = edata->listOfParticles.at(p)->getVector();
191 
192  /* flip z axis component of particle momentum: sPHENIX EIC detector coordinate system uses
193  electron flying in 'negative z' direction */
194  HVector3 flip_pz(v4_particle_p.getVector());
195  flip_pz.setZ(flip_pz.Z() * -1);
196  v4_particle_p.setVector(flip_pz);
197 
198  /* Boost particle from proton-at-rest frame to laboratory frame */
199  HLorentzVector v4_particle_p_lab(v4_particle_p);
200 
201  v4_particle_p_lab.boost(sqrt(_p4_hadron_lab->getQuare()), *_p4_hadron_lab);
202 
203  if (verbosity > 2)
204  {
205  cout << "EVENT RECORD particle: " << edata->listOfParticles.at(p)->getParticleType()
206  << " (status: " << edata->listOfParticles.at(p)->getParticleAuxFlag() << ")" << endl;
207  cout << " --> PROTON REST frame: ";
208  v4_particle_p.print();
209  cout << " --> LABORATORY frame: ";
210  v4_particle_p_lab.print();
211  }
212 
213  HepMC::GenParticle *particle_hepmc = new HepMC::GenParticle(HepMC::FourVector(v4_particle_p_lab.getVector().X(),
214  v4_particle_p_lab.getVector().Y(),
215  v4_particle_p_lab.getVector().Z(),
216  v4_particle_p_lab.getEnergy()),
217  edata->listOfParticles.at(p)->getParticleType());
218  particle_hepmc->set_status(edata->listOfParticles.at(p)->getParticleAuxFlag());
219  hepmcvtx->add_particle_out(particle_hepmc);
220  }
221 
222  /* Add vertex to event */
223  evt->add_vertex(hepmcvtx);
224 
225  /* pass HepMC to PHNode */
226  PHHepMCGenEvent *success = PHHepMCGenHelper::insert_event(evt);
227  if (!success)
228  {
229  cout << "sHEPGen::process_event - Failed to add event to HepMC record!" << endl;
230  return Fun4AllReturnCodes::ABORTRUN;
231  }
232 
233  /* print outs */
234  if (verbosity > 2) cout << "sHEPGen::process_event - FINISHED WHOLE EVENT" << endl;
235 
236  ++_eventcount;
237 
238  return Fun4AllReturnCodes::EVENT_OK;
239 }
240 
241 void sHEPGen::printlogo()
242 {
243  cout << endl
244  << endl;
245  for (int i = 0; i < hepconst::logolength; i++)
246  cout << hepconst::logo[i] << endl;
247  cout << endl
248  << endl;
249 }