generators/sHEPGen/sHEPGen.C

0001 #include "sHEPGen.h"
0002
0003 #include <phhepmc/PHHepMCGenEvent.h>
0004
0005 #include <fun4all/Fun4AllReturnCodes.h>
0006 #include <phool/PHCompositeNode.h>
0007 #include <phool/PHIODataNode.h>
0008 #include <phool/PHNodeIterator.h>
0009 #include <phool/PHRandomSeed.h>
0010
0011 /* HEPGen includes */
0012 #include <hgenmanager.h>
0013
0014 /* HepMC includes */
0015 #include <HepMC/GenEvent.h>
0016
0017 using namespace std;
0018
0019 typedef PHIODataNode<PHObject> PHObjectNode_t;
0020
0021 sHEPGen::sHEPGen(const std::string &name)
0022   : SubsysReco(name)
0023   , _eventcount(0)
0024   , _p_electron_lab(-20)
0025   , _p_hadron_lab(250)
0026   , _p4_electron_lab(nullptr)
0027   , _p4_hadron_lab(nullptr)
0028   , _p4_hadron_lab_invert(nullptr)
0029   , _p4_electron_prest(nullptr)
0030   , _p4_hadron_prest(nullptr)
0031   , _hgenManager(NULL)
0032   , _datacardFile("hepgen_dvcs.data")
0033 {
0034   PHHepMCGenHelper::set_embedding_id(1);  // default embedding ID to 1
0035 }
0036
0037 sHEPGen::~sHEPGen()
0038 {
0039 }
0040
0041 int sHEPGen::Init(PHCompositeNode *topNode)
0042 {
0043   printlogo();
0044
0045   /* electron and proton mass */
0046   double mass_e = 5.109989e-4;
0047   double mass_p = 9.382720e-1;
0048
0049   /* 4-Vectors of colliding electron and proton in laboratory frame */
0050   _p4_electron_lab = new HLorentzVector(0.,
0051                                         0.,
0052                                         _p_electron_lab,
0053                                         sqrt(_p_electron_lab * _p_electron_lab + mass_e * mass_e));
0054
0055   _p4_hadron_lab = new HLorentzVector(0.,
0056                                       0.,
0057                                       _p_hadron_lab,
0058                                       sqrt(_p_hadron_lab * _p_hadron_lab + mass_p * mass_p));
0059
0060   /* The current version of HEPGen supports only a "Fixed Target" mode, i.e.
0061      the target (proton) is assumed at rest. Therefore, need to boost collision
0062      from the laboratory frame to the proton-at-rest frame. */
0063   _p4_hadron_lab_invert = new HLorentzVector(0.,
0064                                              0.,
0065                                              -1 * _p_hadron_lab,
0066                                              sqrt(_p_hadron_lab * _p_hadron_lab + mass_p * mass_p));
0067
0068   /* 4-Vectors of colliding electron and proton in proton-at-rest frame */
0069   _p4_electron_prest = new HLorentzVector(*_p4_electron_lab);
0070   _p4_electron_prest->boost(sqrt(_p4_hadron_lab_invert->getQuare()), *_p4_hadron_lab_invert);
0071
0072   _p4_hadron_prest = new HLorentzVector(*_p4_hadron_lab);
0073   _p4_hadron_prest->boost(sqrt(_p4_hadron_lab_invert->getQuare()), *_p4_hadron_lab_invert);
0074
0075   if (verbosity > 2)
0076   {
0077     cout << "Electron and proton in laboratory frame:" << endl;
0078     _p4_electron_lab->print();
0079     _p4_hadron_lab->print();
0080
0081     cout << "Electron and proton in proton-at-rest frame:" << endl;
0082     _p4_electron_prest->print();
0083     _p4_hadron_prest->print();
0084   }
0085
0086   /* flip sign of electron momentum z-component for HEPGen */
0087   HVector3 flip_pz(_p4_electron_prest->getVector());
0088   flip_pz.setZ(flip_pz.Z() * -1);
0089   _p4_electron_prest->setVector(flip_pz);
0090
0091   /* get instance of HepGenManager */
0092   _hgenManager = HGenManager::getInstance();
0093   _hgenManager->loadSettings(_datacardFile);
0094   _hgenManager->setupGenerator();
0095
0096   /* set beam parameters */
0097   cout << "Colliding " << _p4_electron_lab->getVector().Z() << " GeV electron with " << _p4_hadron_lab->getVector().Z() << " GeV proton (laboratory frame)" << endl;
0098   cout << "----ELEPT (proton-at-rest): " << _p4_electron_prest->getVector().Z() << " GeV" << endl;
0099   _hgenManager->getParamManager()->getStruct()->ELEPT = _p4_electron_prest->getVector().Z();
0100   _hgenManager->getParamManager()->getStruct()->PARL.at(2) = _p4_electron_prest->getVector().Z();
0101
0102   /* set random seed */
0103   unsigned int seed = PHRandomSeed();
0104   _hgenManager->setSeed(seed);
0105
0106   /* enable detailed event record printput for debugging */
0107   if (verbosity > 2)
0108     _hgenManager->enableDebug();
0109
0110   create_node_tree(topNode);
0111
0112   _eventcount = 0;
0113
0114   return Fun4AllReturnCodes::EVENT_OK;
0115 }
0116
0117 int sHEPGen::End(PHCompositeNode *topNode)
0118 {
0119   cout << "Reached the sHEPGen::End()" << endl;
0120
0121   //-* dump out closing info (cross-sections, etc)
0122
0123   if (verbosity > 1) cout << "sHEPGen::End - I'm here!" << endl;
0124
0125   return Fun4AllReturnCodes::EVENT_OK;
0126 }
0127
0128 int sHEPGen::process_event(PHCompositeNode *topNode)
0129 {
0130   if (verbosity > 1) cout << "sHEPGen::process_event - event: " << _eventcount << endl;
0131
0132   _hgenManager->oneShot();
0133
0134   if (verbosity > 2)
0135     _hgenManager->getEvent()->printDebug();
0136
0137   HEvent *evt_mc = _hgenManager->getEvent();
0138
0139   /* Create HepMC GenEvent */
0140   HepMC::GenEvent *evt = new HepMC::GenEvent();
0141
0142   /* define the units (Pythia uses GeV and mm) */
0143   evt->use_units(HepMC::Units::GEV, HepMC::Units::MM);
0144
0145   /* add global information to the event */
0146   evt->set_event_number(_eventcount);
0147
0148   /* Set the PDF information */
0149   HepMC::PdfInfo pdfinfo;
0150   pdfinfo.set_scalePDF(evt_mc->getQsq());
0151   pdfinfo.set_x2(evt_mc->getXbj());
0152   evt->set_pdf_info(pdfinfo);
0153
0154   /* Set additional event parameters */
0155   evt->set_event_scale(evt_mc->getQsq());
0156
0157   /* Event kinematics */
0158   /* @TODO How can we store this information in HepMC record? */
0159   //evt_mc->getNu();
0160   //evt_mc->getY();
0161   //evt_mc->getWsq();
0162   //evt_mc->getXbj();
0163   //evt_mc->getT();
0164   //evt_mc->getQsq();
0165   //evt_mc->getS();
0166   //evt_mc->getEmiss();
0167
0168   /* Create single HepMC vertex for event */
0169   /* @TODO: Implement multiple vertices e.g. for decay particles */
0170   HepMC::GenVertex *hepmcvtx = new HepMC::GenVertex(HepMC::FourVector(0,
0171                                                                       0,
0172                                                                       0,
0173                                                                       0));
0174
0175   /* Create HepMC particle records */
0176   HEventData *edata = evt_mc->getStruct();
0177   for (unsigned p = 0; p < edata->listOfParticles.size(); p++)
0178   {
0179     if (verbosity > 4)
0180     {
0181       cout << "______new particle_______" << endl;
0182       cout << "Index:  " << p + 1 << endl;
0183       cout << "PID: " << edata->listOfParticles.at(p)->getParticleType() << " -- " << (edata->listOfParticles.at(p) == &(edata->incBeamParticle)) << endl;
0184       cout << "Particle aux flag: " << edata->listOfParticles.at(p)->getParticleAuxFlag() << endl;
0185       cout << "Particle origin: " << edata->listOfParticles.at(p)->getParticleOrigin() << endl;
0186       cout << "Particle daughter1: " << edata->listOfParticles.at(p)->getParticleDaughter1() << endl;
0187       cout << "Particle daughter2: " << edata->listOfParticles.at(p)->getParticleDaughter2() << endl;
0188     }
0189
0190     HLorentzVector v4_particle_p = edata->listOfParticles.at(p)->getVector();
0191
0192     /* flip z axis component of particle momentum: sPHENIX EIC detector coordinate system uses
0193          electron flying in 'negative z' direction */
0194     HVector3 flip_pz(v4_particle_p.getVector());
0195     flip_pz.setZ(flip_pz.Z() * -1);
0196     v4_particle_p.setVector(flip_pz);
0197
0198     /* Boost particle from proton-at-rest frame to laboratory frame */
0199     HLorentzVector v4_particle_p_lab(v4_particle_p);
0200
0201     v4_particle_p_lab.boost(sqrt(_p4_hadron_lab->getQuare()), *_p4_hadron_lab);
0202
0203     if (verbosity > 2)
0204     {
0205       cout << "EVENT RECORD particle: " << edata->listOfParticles.at(p)->getParticleType()
0206            << " (status: " << edata->listOfParticles.at(p)->getParticleAuxFlag() << ")" << endl;
0207       cout << " --> PROTON REST frame: ";
0208       v4_particle_p.print();
0209       cout << " --> LABORATORY  frame: ";
0210       v4_particle_p_lab.print();
0211     }
0212
0213     HepMC::GenParticle *particle_hepmc = new HepMC::GenParticle(HepMC::FourVector(v4_particle_p_lab.getVector().X(),
0214                                                                                   v4_particle_p_lab.getVector().Y(),
0215                                                                                   v4_particle_p_lab.getVector().Z(),
0216                                                                                   v4_particle_p_lab.getEnergy()),
0217                                                                 edata->listOfParticles.at(p)->getParticleType());
0218     particle_hepmc->set_status(edata->listOfParticles.at(p)->getParticleAuxFlag());
0219     hepmcvtx->add_particle_out(particle_hepmc);
0220   }
0221
0222   /* Add vertex to event */
0223   evt->add_vertex(hepmcvtx);
0224
0225   /* pass HepMC to PHNode */
0226   PHHepMCGenEvent *success = PHHepMCGenHelper::insert_event(evt);
0227   if (!success)
0228   {
0229     cout << "sHEPGen::process_event - Failed to add event to HepMC record!" << endl;
0230     return Fun4AllReturnCodes::ABORTRUN;
0231   }
0232
0233   /* print outs */
0234   if (verbosity > 2) cout << "sHEPGen::process_event - FINISHED WHOLE EVENT" << endl;
0235
0236   ++_eventcount;
0237
0238   return Fun4AllReturnCodes::EVENT_OK;
0239 }
0240
0241 void sHEPGen::printlogo()
0242 {
0243   cout << endl
0244        << endl;
0245   for (int i = 0; i < hepconst::logolength; i++)
0246     cout << hepconst::logo[i] << endl;
0247   cout << endl
0248        << endl;
0249 }