generators/sphstarlight/StarLight_To_HepMc.cc

0001 #include <HepMC/GenEvent.h>
0002 #include <HepMC/GenParticle.h>
0003 #include <HepMC/GenVertex.h>
0004 #include <HepMC/IO_BaseClass.h>  // for IO_BaseClass
0005 #include <HepMC/IO_GenEvent.h>
0006 #include <HepMC/SimpleVector.h>  // for FourVector
0007 #include <HepMC/Units.h>         // for GEV, MM
0008
0009 #include <CLHEP/Vector/LorentzVector.h>
0010
0011 #include <TDatabasePDG.h>
0012 #include <TParticlePDG.h>  // for TParticlePDG
0013
0014 #include <cmath>
0015 #include <fstream>
0016 #include <iostream>
0017 #include <sstream>
0018 #include <string>
0019
0020 TDatabasePDG* PDGdb;
0021
0022 // Function to parse the STARlight output file and return a vector of HepMC::GenParticle pointers
0023 // std::vector<HepMC::GenParticle*> parseStarlightOutput(const std::string& filename)
0024 int fillEvent(HepMC::GenEvent* evt, std::ifstream& file)
0025 {
0026   std::string line;
0027   std::string label;
0028
0029   int found_event = 0;
0030   while ( getline( file, line ) )
0031   {
0032     // keep going until EVENT line is found
0033     if ( line.starts_with("EVENT") )
0034     {
0035       found_event = 1;
0036       break;
0037     }
0038   }
0039
0040   if ( found_event == 1 )
0041   {
0042     int nevt, ntrk, nvtx;
0043
0044     std::stringstream evtline( line );
0045
0046     evtline >> label >> nevt >> ntrk >> nvtx;
0047
0048     // first line should be the event
0049     // EVENT: n ntracks nvertices ,
0050     if (label != "EVENT:")
0051     {
0052       std::cout << "label is " << label << std::endl;
0053       return -1;
0054     }
0055
0056     static int nprint = 100;
0057     if (nevt % nprint == 0)
0058     {
0059       std::cout << nevt << std::endl;
0060       if ( nevt>10000 )
0061       {
0062         nprint = 10000;
0063       }
0064       else if ( nevt>1000 )
0065       {
0066         nprint = 1000;
0067       }
0068     }
0069
0070     evt->set_event_number(nevt);
0071
0072     // Next line should be the vertex
0073     // Note: starlight currently only has one vertex, but a future version could have more
0074     // VERTEX: x y z t nv nproc nparent ndaughters
0075     double x, y, z, t;
0076     int nv, nproc, npar, ndau;
0077     file >> label >> x >> y >> z >> t >> nv >> nproc >> npar >> ndau;
0078     if (label != "VERTEX:")
0079     {
0080       return -1;
0081     }
0082     CLHEP::HepLorentzVector newVertex;
0083     newVertex = CLHEP::HepLorentzVector(x, y, z, t);
0084     HepMC::GenVertex* v0 = new HepMC::GenVertex(newVertex);
0085     evt->add_vertex(v0);
0086
0087     // TRACK: GPID px py py nev ntr stopv PDGPID
0088     double px, py, pz;  // three vector components of the track's momentum
0089     int gpid, nev, ntr, stopv, pdgpid;
0090     for (int itrk = 0; itrk < ndau; itrk++)
0091     {
0092       file >> label >> gpid >> px >> py >> pz >> nev >> ntr >> stopv >> pdgpid;
0093       if (label != "TRACK:")
0094       {
0095         return -1;
0096       }
0097
0098       // Get mass of particle
0099       TParticlePDG* partinfo = PDGdb->GetParticle(pdgpid);
0100       double mass = partinfo->Mass();
0101
0102       // Calculate energy
0103       double E = sqrt(mass * mass + px * px + py * py + pz * pz);
0104
0105       CLHEP::HepLorentzVector p4(px, py, pz, E);
0106
0107       const int particle_status = 1;
0108
0109       // Create a HepMC::GenParticle using status code 1 for stable particle
0110       v0->add_particle_out(new HepMC::GenParticle(p4, pdgpid, particle_status));
0111     }
0112   }
0113   else
0114   {
0115     // end of file, or failure to read file
0116     return -2;
0117   }
0118
0119   return 1;
0120 }
0121
0122 // Main function to convert STARlight output to HEPMC format
0123 int main(int argc, char* argv[])
0124 {
0125   if (argc != 3)
0126   {
0127     std::cerr << "Usage: " << argv[0] << "<starlight_file> <HepMC_output_file>" << std::endl;
0128     return 1;
0129   }
0130
0131   // std::string starlightfname = "slight.out"; // Slight Input file Hardcoded.
0132   // This way the starlight input file is not hardcoded which is bad.
0133   std::string starlightfname = argv[1];
0134   std::string hepmcfname = argv[2];
0135
0136   // Open the STARlight output file
0137   std::ifstream starlightfile(starlightfname);
0138   if (!starlightfile.is_open())
0139   {
0140     std::cerr << " STARlight file not found: " << starlightfname << std::endl;
0141     return -1;
0142   }
0143
0144   // Open pdg database
0145   PDGdb = TDatabasePDG::Instance();
0146
0147   // Open the HepMC output file
0148   std::ofstream outputFile(hepmcfname);
0149   HepMC::IO_GenEvent ascii_io(hepmcfname.c_str(), std::ios::out);
0150
0151 //  unsigned int events = 0;
0152
0153   while (true)
0154   {
0155 //    events++;
0156
0157     HepMC::GenEvent* evt = new HepMC::GenEvent();
0158     evt->use_units(HepMC::Units::GEV, HepMC::Units::MM);
0159
0160     int status = fillEvent(evt, starlightfile);
0161     if (status < 0)
0162     {
0163       delete evt;
0164       break;
0165     }
0166
0167     ascii_io << evt;
0168
0169     delete evt;
0170   }
0171
0172   std::cout << "Conversion completed successfully. HEPMC file saved to " << hepmcfname << std::endl;
0173
0174   starlightfile.close();
0175
0176   return 0;
0177 }