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