LHEInterface/src/lheh5.cc

0001 #include "GeneratorInterface/LHEInterface/interface/lheh5.h"
0002
0003 namespace lheh5 {
0004
0005   inline std::ostream& operator<<(std::ostream& os, Particle const& p) {
0006     os << "\tpx: " << p.px << " py: " << p.py << " pz: " << p.pz << " e: " << p.e << "\n";
0007     return os;
0008   }
0009
0010   inline std::ostream& operator<<(std::ostream& os, EventHeader const& eh) {
0011     os << "\tnparticles: " << eh.nparticles << " procid: " << eh.pid << " weight: " << eh.weight
0012        << " trials: " << eh.trials << "\n";
0013     os << "\tscale: " << eh.scale << " rscale: " << eh.rscale << " fscale: " << eh.fscale << " aqed: " << eh.aqed
0014        << " aqcd: " << eh.aqcd << "\n";
0015     os << "\tnpLO: " << eh.npLO << " npNLO: " << eh.npNLO << "\n";
0016     return os;
0017   }
0018
0019   Particle Events::mkParticle(size_t idx) const {
0020     return {_vid[idx],
0021             _vstatus[idx],
0022             _vmother1[idx],
0023             _vmother2[idx],
0024             _vcolor1[idx],
0025             _vcolor2[idx],
0026             _vpx[idx],
0027             _vpy[idx],
0028             _vpz[idx],
0029             _ve[idx],
0030             _vm[idx],
0031             _vlifetime[idx],
0032             _vspin[idx]};
0033   }
0034
0035   std::vector<Particle> Events::mkEvent(size_t ievent) const {
0036     std::vector<Particle> _E;
0037     // NOTE we need to subtract the particle offset here as the
0038     // particle properties run from 0 and not the global index when using batches
0039     for (size_t id = _vstart[ievent] - _particle_offset; id < _vend[ievent] - _particle_offset; ++id) {
0040       _E.push_back(mkParticle(id));
0041     }
0042
0043     // Make sure beam particles are ordered according to convention i.e. first particle has positive z-momentum
0044     if (_E[0].pz < 0) {
0045       std::swap<Particle>(_E[0], _E[1]);
0046     }
0047
0048     return _E;
0049   }
0050
0051   EventHeader Events::mkEventHeader(int iev) const {
0052     return {
0053         _vnparticles[iev],
0054         _vpid[iev],
0055         _vweight[iev],
0056         _vtrials[iev],
0057         _vscale[iev],
0058         _vrscale[iev],
0059         _vfscale[iev],
0060         _vaqed[iev],
0061         _vaqcd[iev],
0062         _vnpLO[iev],
0063         _vnpNLO[iev],
0064     };
0065   }
0066
0067   // Read function, returns an Events struct
0068   Events readEvents(HighFive::Group& g_index,
0069                     HighFive::Group& g_particle,
0070                     HighFive::Group& g_event,
0071                     size_t first_event,
0072                     size_t n_events) {
0073     // Lookup
0074     std::vector<size_t> _vstart, _vend;
0075     // Particles
0076     std::vector<int> _vid, _vstatus, _vmother1, _vmother2, _vcolor1, _vcolor2;
0077     std::vector<double> _vpx, _vpy, _vpz, _ve, _vm, _vlifetime, _vspin;
0078     // Event info
0079     std::vector<int> _vnparticles, _vpid, _vnpLO, _vnpNLO;
0080     std::vector<size_t> _vtrials;
0081     std::vector<double> _vweight, _vscale, _vrscale, _vfscale, _vaqed, _vaqcd;
0082
0083     //double starttime, endtime;
0084     //starttime = MPI_Wtime();
0085     // Lookup
0086     HighFive::DataSet _start = g_index.getDataSet("start");
0087     HighFive::DataSet _end = g_index.getDataSet("end");
0088     // Particles
0089     HighFive::DataSet _id = g_particle.getDataSet("id");
0090     HighFive::DataSet _status = g_particle.getDataSet("status");
0091     HighFive::DataSet _mother1 = g_particle.getDataSet("mother1");
0092     HighFive::DataSet _mother2 = g_particle.getDataSet("mother2");
0093     HighFive::DataSet _color1 = g_particle.getDataSet("color1");
0094     HighFive::DataSet _color2 = g_particle.getDataSet("color2");
0095     HighFive::DataSet _px = g_particle.getDataSet("px");
0096     HighFive::DataSet _py = g_particle.getDataSet("py");
0097     HighFive::DataSet _pz = g_particle.getDataSet("pz");
0098     HighFive::DataSet _e = g_particle.getDataSet("e");
0099     HighFive::DataSet _m = g_particle.getDataSet("m");
0100     HighFive::DataSet _lifetime = g_particle.getDataSet("lifetime");
0101     HighFive::DataSet _spin = g_particle.getDataSet("spin");
0102     // Event info
0103     HighFive::DataSet _nparticles = g_event.getDataSet("nparticles");
0104     HighFive::DataSet _pid = g_event.getDataSet("pid");
0105     HighFive::DataSet _weight = g_event.getDataSet("weight");
0106     HighFive::DataSet _trials = g_event.getDataSet("trials");
0107     HighFive::DataSet _scale = g_event.getDataSet("scale");
0108     HighFive::DataSet _rscale = g_event.getDataSet("rscale");
0109     HighFive::DataSet _fscale = g_event.getDataSet("fscale");
0110     HighFive::DataSet _aqed = g_event.getDataSet("aqed");
0111     HighFive::DataSet _aqcd = g_event.getDataSet("aqcd");
0112     HighFive::DataSet _npLO = g_event.getDataSet("npLO");
0113     HighFive::DataSet _npNLO = g_event.getDataSet("npNLO");
0114
0115     //endtime   = MPI_Wtime();
0116     //printf("DS took %f seconds\n", endtime-starttime);
0117     std::vector<size_t> offset_e = {first_event};
0118     std::vector<size_t> readsize_e = {n_events};
0119     //_vstart.reserve(n_events);
0120     //_vend.reserve(n_events);
0121
0122     _start.select(offset_e, readsize_e).read(_vstart);
0123     _end.select(offset_e, readsize_e).read(_vend);
0124     std::vector<size_t> offset_p = {_vstart.front()};
0125     std::vector<size_t> readsize_p = {_vend.back() - _vstart.front()};
0126
0127     int RESP = _vend.back() - _vstart.front();
0128     _vid.reserve(RESP);
0129     _vstatus.reserve(RESP);
0130     _vmother1.reserve(RESP);
0131     _vmother2.reserve(RESP);
0132     _vcolor1.reserve(RESP);
0133     _vcolor2.reserve(RESP);
0134     _vpx.reserve(RESP);
0135     _vpy.reserve(RESP);
0136     _vpz.reserve(RESP);
0137     _ve.reserve(RESP);
0138     _vm.reserve(RESP);
0139     _vlifetime.reserve(RESP);
0140     _vspin.reserve(RESP);
0141
0142     _vnparticles.reserve(n_events);
0143     _vpid.reserve(n_events);
0144     _vweight.reserve(n_events);
0145     _vtrials.reserve(n_events);
0146     _vscale.reserve(n_events);
0147     _vrscale.reserve(n_events);
0148     _vfscale.reserve(n_events);
0149     _vaqed.reserve(n_events);
0150     _vaqcd.reserve(n_events);
0151     _vnpLO.reserve(n_events);
0152     _vnpNLO.reserve(n_events);
0153
0154     //starttime = MPI_Wtime();
0155     // This is using HighFive's read
0156     _id.select(offset_p, readsize_p).read(_vid);
0157     _status.select(offset_p, readsize_p).read(_vstatus);
0158     _mother1.select(offset_p, readsize_p).read(_vmother1);
0159     _mother2.select(offset_p, readsize_p).read(_vmother2);
0160     _color1.select(offset_p, readsize_p).read(_vcolor1);
0161     _color2.select(offset_p, readsize_p).read(_vcolor2);
0162     _px.select(offset_p, readsize_p).read(_vpx);
0163     _py.select(offset_p, readsize_p).read(_vpy);
0164     _pz.select(offset_p, readsize_p).read(_vpz);
0165     _e.select(offset_p, readsize_p).read(_ve);
0166     _m.select(offset_p, readsize_p).read(_vm);
0167     _lifetime.select(offset_p, readsize_p).read(_vlifetime);
0168     _spin.select(offset_p, readsize_p).read(_vspin);
0169
0170     //endtime   = MPI_Wtime();
0171     //printf("SELP took %f seconds\n", endtime-starttime);
0172     //starttime = MPI_Wtime();
0173     _nparticles.select(offset_e, readsize_e).read(_vnparticles);
0174     _pid.select(offset_e, readsize_e).read(_vpid);
0175     _weight.select(offset_e, readsize_e).read(_vweight);
0176     _trials.select(offset_e, readsize_e).read(_vtrials);
0177     _scale.select(offset_e, readsize_e).read(_vscale);
0178     _rscale.select(offset_e, readsize_e).read(_vrscale);
0179     _fscale.select(offset_e, readsize_e).read(_vfscale);
0180     _aqed.select(offset_e, readsize_e).read(_vaqed);
0181     _aqcd.select(offset_e, readsize_e).read(_vaqcd);
0182     _npLO.select(offset_e, readsize_e).read(_vnpLO);
0183     _npNLO.select(offset_e, readsize_e).read(_vnpNLO);
0184     //endtime   = MPI_Wtime();
0185     //printf("SELE took %f seconds\n", endtime-starttime);
0186
0187     return {
0188         std::move(_vstart),      std::move(_vend),    std::move(_vid),     std::move(_vstatus),   std::move(_vmother1),
0189         std::move(_vmother2),    std::move(_vcolor1), std::move(_vcolor2), std::move(_vpx),       std::move(_vpy),
0190         std::move(_vpz),         std::move(_ve),      std::move(_vm),      std::move(_vlifetime), std::move(_vspin),
0191         std::move(_vnparticles), std::move(_vpid),    std::move(_vweight), std::move(_vtrials),   std::move(_vscale),
0192         std::move(_vrscale),     std::move(_vfscale), std::move(_vaqed),   std::move(_vaqcd),     std::move(_vnpLO),
0193         std::move(_vnpNLO),      offset_p[0],
0194     };
0195   }
0196
0197   Particle Events2::mkParticle(size_t idx) const {
0198     return {_vid[idx],
0199             _vstatus[idx],
0200             _vmother1[idx],
0201             _vmother2[idx],
0202             _vcolor1[idx],
0203             _vcolor2[idx],
0204             _vpx[idx],
0205             _vpy[idx],
0206             _vpz[idx],
0207             _ve[idx],
0208             _vm[idx],
0209             _vlifetime[idx],
0210             _vspin[idx]};
0211   }
0212
0213   std::vector<Particle> Events2::mkEvent(size_t ievent) const {
0214     std::vector<Particle> _E;
0215     // NOTE we need to subtract the particle offset here as the
0216     // particle properties run from 0 and not the global index when using batches
0217     size_t partno = _vstart[ievent] - _particle_offset;
0218     for (int id = 0; id < _vnparticles[ievent]; ++id) {
0219       _E.push_back(mkParticle(partno));
0220       partno++;
0221     }
0222
0223     // Make sure beam particles are ordered according to convention i.e. first particle has positive z-momentum
0224     if (_E[0].pz < 0) {
0225       std::swap<Particle>(_E[0], _E[1]);
0226     }
0227
0228     return _E;
0229   }
0230
0231   EventHeader Events2::mkEventHeader(int iev) const {
0232     return {
0233         _vnparticles[iev],
0234         _vpid[iev],
0235         _vweight[iev],
0236         _vtrials[iev],
0237         _vscale[iev],
0238         _vrscale[iev],
0239         _vfscale[iev],
0240         _vaqed[iev],
0241         _vaqcd[iev],
0242         npLO,
0243         npNLO,
0244     };
0245   }
0246
0247   // Read function, returns an Events struct --- this is for the new structure
0248   Events2 readEvents(
0249       HighFive::Group& g_particle, HighFive::Group& g_event, size_t first_event, size_t n_events, int npLO, int npNLO) {
0250     // Lookup
0251     std::vector<size_t> _vstart;
0252     // Particles
0253     std::vector<int> _vid, _vstatus, _vmother1, _vmother2, _vcolor1, _vcolor2;
0254     std::vector<double> _vpx, _vpy, _vpz, _ve, _vm, _vlifetime, _vspin;
0255     // Event info
0256     std::vector<int> _vnparticles, _vpid;
0257     std::vector<size_t> _vtrials;
0258     std::vector<double> _vweight, _vscale, _vrscale, _vfscale, _vaqed, _vaqcd;
0259
0260     //double starttime, endtime;
0261     // Lookup
0262     HighFive::DataSet _start = g_event.getDataSet("start");
0263     // Particles
0264     HighFive::DataSet _id = g_particle.getDataSet("id");
0265     HighFive::DataSet _status = g_particle.getDataSet("status");
0266     HighFive::DataSet _mother1 = g_particle.getDataSet("mother1");
0267     HighFive::DataSet _mother2 = g_particle.getDataSet("mother2");
0268     HighFive::DataSet _color1 = g_particle.getDataSet("color1");
0269     HighFive::DataSet _color2 = g_particle.getDataSet("color2");
0270     HighFive::DataSet _px = g_particle.getDataSet("px");
0271     HighFive::DataSet _py = g_particle.getDataSet("py");
0272     HighFive::DataSet _pz = g_particle.getDataSet("pz");
0273     HighFive::DataSet _e = g_particle.getDataSet("e");
0274     HighFive::DataSet _m = g_particle.getDataSet("m");
0275     HighFive::DataSet _lifetime = g_particle.getDataSet("lifetime");
0276     HighFive::DataSet _spin = g_particle.getDataSet("spin");
0277     // Event info
0278     HighFive::DataSet _nparticles = g_event.getDataSet("nparticles");
0279     HighFive::DataSet _pid = g_event.getDataSet("pid");
0280     HighFive::DataSet _weight = g_event.getDataSet("weight");
0281     HighFive::DataSet _trials = g_event.getDataSet("trials");
0282     HighFive::DataSet _scale = g_event.getDataSet("scale");
0283     HighFive::DataSet _rscale = g_event.getDataSet("rscale");
0284     HighFive::DataSet _fscale = g_event.getDataSet("fscale");
0285     HighFive::DataSet _aqed = g_event.getDataSet("aqed");
0286     HighFive::DataSet _aqcd = g_event.getDataSet("aqcd");
0287
0288     std::vector<size_t> offset_e = {first_event};
0289     std::vector<size_t> readsize_e = {n_events};
0290
0291     // We now know the first event to read
0292     _start.select(offset_e, readsize_e).read(_vstart);
0293
0294     // That's the first particle
0295     std::vector<size_t> offset_p = {_vstart.front()};
0296     // The last particle is last entry in start + nparticles of that event
0297     _vnparticles.reserve(n_events);
0298     _nparticles.select(offset_e, readsize_e).read(_vnparticles);
0299
0300     size_t RESP = _vstart.back() - _vstart.front() + _vnparticles.back();
0301     std::vector<size_t> readsize_p = {RESP};
0302
0303     _vid.reserve(RESP);
0304     _vstatus.reserve(RESP);
0305     _vmother1.reserve(RESP);
0306     _vmother2.reserve(RESP);
0307     _vcolor1.reserve(RESP);
0308     _vcolor2.reserve(RESP);
0309     _vpx.reserve(RESP);
0310     _vpy.reserve(RESP);
0311     _vpz.reserve(RESP);
0312     _ve.reserve(RESP);
0313     _vm.reserve(RESP);
0314     _vlifetime.reserve(RESP);
0315     _vspin.reserve(RESP);
0316
0317     _vpid.reserve(n_events);
0318     _vweight.reserve(n_events);
0319     _vtrials.reserve(n_events);
0320     _vscale.reserve(n_events);
0321     _vrscale.reserve(n_events);
0322     _vfscale.reserve(n_events);
0323     _vaqed.reserve(n_events);
0324     _vaqcd.reserve(n_events);
0325
0326     // This is using HighFive's read
0327     _id.select(offset_p, readsize_p).read(_vid);
0328     _status.select(offset_p, readsize_p).read(_vstatus);
0329     _mother1.select(offset_p, readsize_p).read(_vmother1);
0330     _mother2.select(offset_p, readsize_p).read(_vmother2);
0331     _color1.select(offset_p, readsize_p).read(_vcolor1);
0332     _color2.select(offset_p, readsize_p).read(_vcolor2);
0333     _px.select(offset_p, readsize_p).read(_vpx);
0334     _py.select(offset_p, readsize_p).read(_vpy);
0335     _pz.select(offset_p, readsize_p).read(_vpz);
0336     _e.select(offset_p, readsize_p).read(_ve);
0337     _m.select(offset_p, readsize_p).read(_vm);
0338     _lifetime.select(offset_p, readsize_p).read(_vlifetime);
0339     _spin.select(offset_p, readsize_p).read(_vspin);
0340
0341     _pid.select(offset_e, readsize_e).read(_vpid);
0342     _weight.select(offset_e, readsize_e).read(_vweight);
0343     _trials.select(offset_e, readsize_e).read(_vtrials);
0344     _scale.select(offset_e, readsize_e).read(_vscale);
0345     _rscale.select(offset_e, readsize_e).read(_vrscale);
0346     _fscale.select(offset_e, readsize_e).read(_vfscale);
0347     _aqed.select(offset_e, readsize_e).read(_vaqed);
0348     _aqcd.select(offset_e, readsize_e).read(_vaqcd);
0349
0350     return {
0351         std::move(_vstart),
0352         std::move(_vid),
0353         std::move(_vstatus),
0354         std::move(_vmother1),
0355         std::move(_vmother2),
0356         std::move(_vcolor1),
0357         std::move(_vcolor2),
0358         std::move(_vpx),
0359         std::move(_vpy),
0360         std::move(_vpz),
0361         std::move(_ve),
0362         std::move(_vm),
0363         std::move(_vlifetime),
0364         std::move(_vspin),
0365         std::move(_vnparticles),
0366         std::move(_vpid),
0367         std::move(_vweight),
0368         std::move(_vtrials),
0369         std::move(_vscale),
0370         std::move(_vrscale),
0371         std::move(_vfscale),
0372         std::move(_vaqed),
0373         std::move(_vaqcd),
0374         npLO,
0375         npNLO,
0376         offset_p[0],
0377     };
0378   }
0379
0380 }  // namespace lheh5