g4simulation/g4main/ReadEICFiles.cc

0001 #include "ReadEICFiles.h"
0002
0003 #include "EicEventHeader.h"
0004 #include "EicEventHeaderv1.h"
0005
0006 #include <fun4all/Fun4AllReturnCodes.h>
0007
0008 #include <phool/PHCompositeNode.h>
0009 #include <phool/PHIODataNode.h>
0010 #include <phool/PHNode.h>  // for PHNode
0011 #include <phool/PHNodeIterator.h>
0012 #include <phool/PHObject.h>  // for PHObject
0013 #include <phool/getClass.h>
0014 #include <phool/phool.h>  // for PHWHERE
0015
0016 #pragma GCC diagnostic push
0017 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0018 #include <HepMC/GenEvent.h>
0019 #include <HepMC/GenVertex.h>
0020 #pragma GCC diagnostic pop
0021
0022 #include <HepMC/GenParticle.h>   // for GenParticle
0023 #include <HepMC/PdfInfo.h>       // for PdfInfo
0024 #include <HepMC/SimpleVector.h>  // for FourVector
0025 #include <HepMC/Units.h>         // for GEV, MM
0026
0027 // eicsmear classes
0028 #include <eicsmear/erhic/EventMC.h>
0029 #include <eicsmear/erhic/EventMilou.h>
0030 #include <eicsmear/erhic/ParticleMC.h>  // for ParticleMC
0031 #include <eicsmear/erhic/Pid.h>         // for Pid
0032
0033 #include <eicsmear/erhic/EventDEMP.h>
0034
0035 // General Root and C++ classes
0036 #include <TBranch.h>  // for TBranch
0037 #include <TChain.h>
0038 #include <TVector3.h>  // for TVector3
0039
0040 #include <iostream>  // for operator<<, endl, basic_ostream
0041 #include <vector>    // for vector
0042
0043 class PHHepMCGenEvent;
0044
0045 using namespace std;
0046
0047 ///////////////////////////////////////////////////////////////////
0048
0049 ReadEICFiles::ReadEICFiles(const string &name)
0050   : SubsysReco(name)
0051   , filename("")
0052   , Tin(nullptr)
0053   , nEntries(0)
0054   , entry(0)
0055   , m_EvtGenId(EvtGen::Unknown)
0056   , GenEvent(nullptr)
0057   , _node_name("PHHepMCGenEvent")
0058 {
0059   PHHepMCGenHelper::set_embedding_id(1);  // default embedding ID to 1
0060   return;
0061 }
0062
0063 ///////////////////////////////////////////////////////////////////
0064
0065 ReadEICFiles::~ReadEICFiles()
0066 {
0067   delete Tin;
0068 }
0069
0070 ///////////////////////////////////////////////////////////////////
0071
0072 bool ReadEICFiles::OpenInputFile(const string &name)
0073 {
0074   filename = name;
0075   Tin = new TChain("EICTree", "EICTree");
0076   Tin->Add(name.c_str());
0077   GetTree();
0078   return true;
0079 }
0080
0081 ///////////////////////////////////////////////////////////////////
0082
0083 void ReadEICFiles::GetTree()
0084 {
0085   /* Print the actual class of the event branch record,
0086      i.e. erhic::EventMilou or other */
0087   string EventClass(Tin->GetBranch("event")->GetClassName());
0088   if (Verbosity() > 0)
0089   {
0090     cout << "ReadEICFiles: Input Branch Event Class = "
0091          << EventClass << endl;
0092   }
0093   if (EventClass.find("Milou") != string::npos)
0094   {
0095     m_EvtGenId = EvtGen::Milou;
0096   }
0097
0098   if (EventClass.find("DEMP") != string::npos)
0099   {
0100     m_EvtGenId = EvtGen::DEMP;
0101   }
0102
0103   Tin->SetBranchAddress("event", &GenEvent);
0104   nEntries = Tin->GetEntries();
0105 }
0106
0107 int ReadEICFiles::Init(PHCompositeNode *topNode)
0108 {
0109   /* Create node tree */
0110   CreateNodeTree(topNode);
0111   return 0;
0112 }
0113
0114 int ReadEICFiles::process_event(PHCompositeNode *topNode)
0115 {
0116   /* Check if there is an unused event left in input file */
0117   if (entry >= nEntries)
0118   {
0119     if (filename.size() > 0)
0120     {
0121       cout << "input file " << filename << " exhausted" << endl;
0122     }
0123     else
0124     {
0125       cout << "no file opened" << endl;
0126     }
0127     return Fun4AllReturnCodes::ABORTRUN;
0128   }
0129
0130   /* Get event record from input file */
0131   Tin->GetEntry(entry);
0132
0133   EicEventHeader *evthead = findNode::getClass<EicEventHeader>(topNode, "EicEventHeader");
0134   switch (m_EvtGenId)
0135   {
0136   case EvtGen::Milou:
0137   {
0138     erhic::EventMilou *gen = dynamic_cast<erhic::EventMilou *>(GenEvent);
0139     evthead->set_eventgenerator_type(EicEventHeader::EvtGen::Milou);
0140     evthead->set_milou_weight(gen->weight);
0141     evthead->set_milou_trueX(gen->trueX);
0142     evthead->set_milou_trueQ2(gen->trueQ2);
0143   }
0144   break;
0145   case EvtGen::DEMP:
0146   {
0147     erhic::EventDEMP *gen = dynamic_cast<erhic::EventDEMP *>(GenEvent);
0148     evthead->set_eventgenerator_type(EicEventHeader::EvtGen::DEMP);
0149     evthead->set_demp_weight(gen->weight);
0150   }
0151   break;
0152   case EvtGen::Unknown:
0153     cout << "unknown event generator" << endl;
0154     break;
0155   default:
0156     cout << "what is this " << m_EvtGenId << " ????" << endl;
0157     break;
0158   }
0159
0160   /* Create GenEvent */
0161   HepMC::GenEvent *evt = new HepMC::GenEvent();
0162
0163   /* define the units (Pythia uses GeV and mm) */
0164   evt->use_units(HepMC::Units::GEV, HepMC::Units::MM);
0165
0166   /* add global information to the event */
0167   evt->set_event_number(entry);
0168
0169   /* process ID from pythia */
0170   evt->set_signal_process_id(GenEvent->GetProcess());
0171
0172   /* Set the PDF information */
0173   HepMC::PdfInfo pdfinfo;
0174   pdfinfo.set_x1(1);
0175   pdfinfo.set_x2(GenEvent->GetX());
0176   pdfinfo.set_scalePDF(GenEvent->GetQ2());
0177   evt->set_pdf_info(pdfinfo);
0178
0179   /* Loop over all particles for this event in input file and fill
0180    * vector with HepMC particles */
0181   vector<HepMC::GenParticle *> hepmc_particles;
0182   vector<unsigned> origin_index;
0183
0184   /* save pointers to beam particles */
0185   HepMC::GenParticle *hepmc_beam1 = nullptr;
0186   HepMC::GenParticle *hepmc_beam2 = nullptr;
0187
0188   for (unsigned ii = 0; ii < GenEvent->GetNTracks(); ii++)
0189   {
0190     /* Get particle / track from event records.
0191      * Class documentation for erhic::VirtualParticle at
0192      * http://www.star.bnl.gov/~tpb/eic-smear/classerhic_1_1_virtual_particle.html */
0193     erhic::ParticleMC *track_ii = GenEvent->GetTrack(ii);
0194
0195     if (Verbosity() > 1)
0196     {
0197       cout << __PRETTY_FUNCTION__ << " : " << __LINE__ << endl;
0198       cout << "\ttrack " << ii << endl;
0199       cout << "\t4mom\t" << track_ii->GetPx() << "\t" << track_ii->GetPy() << "\t" << track_ii->GetPz() << "\t" << track_ii->GetE() << endl;
0200       cout << "\tstatus= " << track_ii->GetStatus() << "\tindex= " << track_ii->GetIndex() << "\tmass= " << track_ii->GetM() << endl;
0201     }
0202
0203     /* Create HepMC particle record */
0204     HepMC::GenParticle *hepmcpart = new HepMC::GenParticle(HepMC::FourVector(track_ii->GetPx(),
0205                                                                              track_ii->GetPy(),
0206                                                                              track_ii->GetPz(),
0207                                                                              track_ii->GetE()),
0208                                                            track_ii->Id());
0209
0210     /* translate eic-smear status codes to HepMC status codes */
0211     switch (track_ii->GetStatus())
0212     {
0213     case 1:
0214       hepmcpart->set_status(1);
0215       break;  // 'stable particle'
0216
0217     case 21:
0218       hepmcpart->set_status(3);
0219       break;  // 'documentation line'
0220
0221     default:
0222       hepmcpart->set_status(0);
0223       break;  // 'null entry'
0224     }
0225
0226     /* assume the first two particles are the beam particles (which getsHepMC status 4)*/
0227     if (ii < 2)
0228     {
0229       hepmcpart->set_status(4);
0230     }
0231
0232     /* add particle information */
0233     hepmcpart->setGeneratedMass(track_ii->GetM());
0234
0235     /* append particle to vector */
0236     hepmc_particles.push_back(hepmcpart);
0237     origin_index.push_back(track_ii->GetIndex());
0238
0239     /* if first particle, call this the first beam particle */
0240     if (ii == 0)
0241     {
0242       hepmc_beam1 = hepmcpart;
0243     }
0244
0245     /* if second particle, call this the second beam particle */
0246     if (ii == 1)
0247     {
0248       hepmc_beam2 = hepmcpart;
0249     }
0250   }
0251
0252   /* Check if hepmc_particles and origin_index vectors are the same size */
0253   if (hepmc_particles.size() != origin_index.size())
0254   {
0255     cout << "ReadEICFiles::process_event - Lengths of HepMC particles and Origin index vectors do not match!" << endl;
0256
0257     delete evt;
0258     return Fun4AllReturnCodes::ABORTRUN;
0259   }
0260
0261   /* add HepMC particles to Hep MC vertices; skip first two particles
0262    * in loop, assuming that they are the beam particles */
0263   vector<HepMC::GenVertex *> hepmc_vertices;
0264
0265   for (unsigned p = 2; p < hepmc_particles.size(); p++)
0266   {
0267     HepMC::GenParticle *pp = hepmc_particles.at(p);
0268
0269     /* continue if vertices for particle are already set */
0270     if (pp->production_vertex() && pp->end_vertex())
0271     {
0272       continue;
0273     }
0274
0275     /* access mother particle vertex */
0276     erhic::ParticleMC *track_pp = GenEvent->GetTrack(p);
0277
0278     unsigned parent_index = track_pp->GetParentIndex();
0279
0280     HepMC::GenParticle *pmother = nullptr;
0281     for (unsigned m = 0; m < hepmc_particles.size(); m++)
0282     {
0283       if (origin_index.at(m) == parent_index)
0284       {
0285         pmother = hepmc_particles.at(m);
0286         break;
0287       }
0288     }
0289
0290     /* if mother does not exist: create new vertex and add this particle as outgoing to vertex */
0291     if (!pmother)
0292     {
0293       HepMC::GenVertex *hepmcvtx = new HepMC::GenVertex(HepMC::FourVector(track_pp->GetVertex()[0],
0294                                                                           track_pp->GetVertex()[1],
0295                                                                           track_pp->GetVertex()[2],
0296                                                                           0));
0297       hepmc_vertices.push_back(hepmcvtx);
0298       hepmcvtx->add_particle_out(pp);
0299       continue;
0300     }
0301     /* if mother exists and has end vertex: add this particle as outgoing to the mother's end vertex */
0302     else if (pmother->end_vertex())
0303     {
0304       pmother->end_vertex()->add_particle_out(pp);
0305     }
0306     /* if mother exists and has no end vertex: create new vertex */
0307     else
0308     {
0309       HepMC::GenVertex *hepmcvtx = new HepMC::GenVertex(HepMC::FourVector(track_pp->GetVertex()[0],
0310                                                                           track_pp->GetVertex()[1],
0311                                                                           track_pp->GetVertex()[2],
0312                                                                           0));
0313       hepmc_vertices.push_back(hepmcvtx);
0314       hepmcvtx->add_particle_in(pmother);
0315       pmother->end_vertex()->add_particle_out(pp);
0316     }
0317   }
0318
0319   /* Add end vertex to beam particles if they don't have one yet */
0320   for (unsigned p = 0; p < 2; p++)
0321   {
0322     HepMC::GenParticle *pp = hepmc_particles.at(p);
0323
0324     if (!pp->end_vertex())
0325     {
0326       /* create collision vertex */
0327       HepMC::GenVertex *hepmcvtx = new HepMC::GenVertex(HepMC::FourVector(0,
0328                                                                           0,
0329                                                                           0,
0330                                                                           0));
0331       hepmc_vertices.push_back(hepmcvtx);
0332       hepmcvtx->add_particle_in(pp);
0333     }
0334   }
0335
0336   /* Check that all particles (except beam particles) have a production vertex */
0337   for (unsigned p = 2; p < hepmc_particles.size(); p++)
0338   {
0339     if (!hepmc_particles.at(p)->production_vertex())
0340     {
0341       cout << "ReadEICFiles::process_event - Missing production vertex for one or more non-beam particles!" << endl;
0342       delete evt;
0343       return Fun4AllReturnCodes::ABORTRUN;
0344     }
0345   }
0346
0347   /* Add HepMC vertices to event */
0348   for (auto &hepmc_vertice : hepmc_vertices)
0349   {
0350     evt->add_vertex(hepmc_vertice);
0351   }
0352
0353   /* set beam particles */
0354   evt->set_beam_particles(hepmc_beam1, hepmc_beam2);
0355
0356   /* pass HepMC to PHNode*/
0357   PHHepMCGenEvent *success = PHHepMCGenHelper::insert_event(evt);
0358   if (Verbosity() > 1)
0359   {
0360     cout << __PRETTY_FUNCTION__ << " : " << __LINE__ << endl;
0361     evt->print();
0362     cout << endl
0363          << endl;
0364   }
0365
0366   if (!success)
0367   {
0368     cout << "ReadEICFiles::process_event - Failed to add event to HepMC record!" << endl;
0369     return Fun4AllReturnCodes::ABORTRUN;
0370   }
0371   /* Count up number of 'used' events from input file */
0372   entry++;
0373
0374   /* Done */
0375   return 0;
0376 }
0377
0378 int ReadEICFiles::CreateNodeTree(PHCompositeNode *topNode)
0379 {
0380   PHHepMCGenHelper::create_node_tree(topNode);
0381
0382   PHNodeIterator iter(topNode);
0383   // Looking for the DST node
0384   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0385   if (!dstNode)
0386   {
0387     cout << PHWHERE << "DST Node missing, doing nothing." << endl;
0388     return Fun4AllReturnCodes::ABORTRUN;
0389   }
0390   EicEventHeader *evthead = findNode::getClass<EicEventHeader>(topNode, "EicEventHeader");
0391   if (!evthead)
0392   {
0393     evthead = new EicEventHeaderv1();
0394     PHIODataNode<PHObject> *HeaderNode = new PHIODataNode<PHObject>(evthead, "EicEventHeader", "PHObject");
0395     dstNode->addNode(HeaderNode);
0396   }
0397   return Fun4AllReturnCodes::EVENT_OK;
0398 }