Line data Source code
1 : // LesHouches.cc is a part of the PYTHIA event generator.
2 : // Copyright (C) 2015 Torbjorn Sjostrand.
3 : // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4 : // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 :
6 : // Function definitions (not found in the header) for the LHAup and
7 : // LHAupLHEF classes.
8 :
9 : #include "Pythia8/LesHouches.h"
10 :
11 : // Access time information.
12 : #include <ctime>
13 :
14 : namespace Pythia8 {
15 :
16 : //==========================================================================
17 :
18 : // LHAup class.
19 :
20 : //--------------------------------------------------------------------------
21 :
22 : // Constants: could be changed here if desired, but normally should not.
23 : // These are of technical nature, as described for each.
24 :
25 : // LHA convention with cross section in pb may require conversion from mb.
26 : const double LHAup::CONVERTMB2PB = 1e9;
27 :
28 : //--------------------------------------------------------------------------
29 :
30 : // Print the initialization info; to check it worked.
31 :
32 : void LHAup::listInit(ostream& os) {
33 :
34 : // Header.
35 0 : os << "\n -------- LHA initialization information ------------ \n";
36 :
37 : // Beam info.
38 0 : os << fixed << setprecision(3)
39 0 : << "\n beam kind energy pdfgrp pdfset \n"
40 0 : << " A " << setw(6) << idBeamASave
41 0 : << setw(12) << eBeamASave
42 0 : << setw(8) << pdfGroupBeamASave
43 0 : << setw(8) << pdfSetBeamASave << "\n"
44 0 : << " B " << setw(6) << idBeamBSave
45 0 : << setw(12) << eBeamBSave
46 0 : << setw(8) << pdfGroupBeamBSave
47 0 : << setw(8) << pdfSetBeamBSave << "\n";
48 :
49 : // Event weighting strategy.
50 0 : os << "\n Event weighting strategy = " << setw(2)
51 0 : << strategySave << "\n" ;
52 :
53 : // Process list.
54 0 : os << scientific << setprecision(4)
55 0 : << "\n Processes, with strategy-dependent cross section info \n"
56 0 : << " number xsec (pb) xerr (pb) xmax (pb) \n" ;
57 0 : for (int ip = 0; ip < int(processes.size()); ++ip) {
58 0 : os << setw(8) << processes[ip].idProc
59 0 : << setw(15) << processes[ip].xSecProc
60 0 : << setw(15) << processes[ip].xErrProc
61 0 : << setw(15) << processes[ip].xMaxProc << "\n";
62 : }
63 :
64 : // Finished.
65 0 : os << "\n -------- End LHA initialization information -------- \n";
66 :
67 0 : }
68 :
69 : //--------------------------------------------------------------------------
70 :
71 : // Print the event info; to check it worked.
72 :
73 : void LHAup::listEvent(ostream& os) {
74 :
75 : // Header.
76 0 : os << "\n -------- LHA event information and listing -------------"
77 0 : << "--------------------------------------------------------- \n";
78 :
79 : // Basic event info.
80 0 : os << scientific << setprecision(4)
81 0 : << "\n process = " << setw(8) << idProc
82 0 : << " weight = " << setw(12) << weightProc
83 0 : << " scale = " << setw(12) << scaleProc << " (GeV) \n"
84 0 : << " "
85 0 : << " alpha_em = " << setw(12) << alphaQEDProc
86 0 : << " alpha_strong = " << setw(12) << alphaQCDProc << "\n";
87 :
88 : // Particle list
89 0 : os << fixed << setprecision(3)
90 0 : << "\n Participating Particles \n"
91 0 : << " no id stat mothers colours p_x "
92 0 : << "p_y p_z e m tau spin \n" ;
93 0 : for (int ip = 1; ip < int(particles.size()); ++ip) {
94 0 : os << setw(6) << ip
95 0 : << setw(10) << particles[ip].idPart
96 0 : << setw(5) << particles[ip].statusPart
97 0 : << setw(6) << particles[ip].mother1Part
98 0 : << setw(6) << particles[ip].mother2Part
99 0 : << setw(6) << particles[ip].col1Part
100 0 : << setw(6) << particles[ip].col2Part
101 0 : << setw(11) << particles[ip].pxPart
102 0 : << setw(11) << particles[ip].pyPart
103 0 : << setw(11) << particles[ip].pzPart
104 0 : << setw(11) << particles[ip].ePart
105 0 : << setw(11) << particles[ip].mPart
106 0 : << setw(8) << particles[ip].tauPart
107 0 : << setw(8) << particles[ip].spinPart << "\n";
108 : }
109 :
110 : // PDF info - optional.
111 0 : if (pdfIsSetSave) os << "\n pdf: id1 =" << setw(5) << id1pdfSave
112 0 : << " id2 =" << setw(5) << id2pdfSave
113 0 : << " x1 =" << scientific << setw(10) << x1pdfSave
114 0 : << " x2 =" << setw(10) << x2pdfSave
115 0 : << " scalePDF =" << setw(10) << scalePDFSave
116 0 : << " pdf1 =" << setw(10) << pdf1Save
117 0 : << " pdf2 =" << setw(10) << pdf2Save << "\n";
118 :
119 : // Finished.
120 0 : os << "\n -------- End LHA event information and listing ---------"
121 0 : << "--------------------------------------------------------- \n";
122 :
123 0 : }
124 :
125 : //--------------------------------------------------------------------------
126 :
127 : // Open and write header to a Les Houches Event File.
128 :
129 : bool LHAup::openLHEF(string fileNameIn) {
130 :
131 : // Open file for writing. Reset it to be empty.
132 0 : fileName = fileNameIn;
133 0 : const char* cstring = fileName.c_str();
134 0 : osLHEF.open(cstring, ios::out | ios::trunc);
135 0 : if (!osLHEF) {
136 0 : infoPtr->errorMsg("Error in LHAup::openLHEF:"
137 0 : " could not open file", fileName);
138 0 : return false;
139 : }
140 :
141 : // Read out current date and time.
142 0 : time_t t = time(0);
143 0 : strftime(dateNow,12,"%d %b %Y",localtime(&t));
144 0 : strftime(timeNow,9,"%H:%M:%S",localtime(&t));
145 :
146 : // Write header.
147 0 : osLHEF << "<LesHouchesEvents version=\"1.0\">\n"
148 0 : << "<!--\n"
149 0 : << " File written by Pythia8::LHAup on "
150 0 : << dateNow << " at " << timeNow << "\n"
151 0 : << "-->" << endl;
152 :
153 : // Done.
154 : return true;
155 :
156 0 : }
157 :
158 : //--------------------------------------------------------------------------
159 :
160 : // Write initialization information to a Les Houches Event File.
161 :
162 : bool LHAup::initLHEF() {
163 :
164 : // Write information on beams.
165 0 : osLHEF << "<init>\n" << scientific << setprecision(6)
166 0 : << " " << idBeamASave << " " << idBeamBSave
167 0 : << " " << eBeamASave << " " << eBeamBSave
168 0 : << " " << pdfGroupBeamASave << " " << pdfGroupBeamBSave
169 0 : << " " << pdfSetBeamASave << " " << pdfSetBeamBSave
170 0 : << " " << strategySave << " " << processes.size() << "\n";
171 :
172 : // Write information on all the subprocesses.
173 0 : for (int ip = 0; ip < int(processes.size()); ++ip)
174 0 : osLHEF << " " << setw(13) << processes[ip].xSecProc
175 0 : << " " << setw(13) << processes[ip].xErrProc
176 0 : << " " << setw(13) << processes[ip].xMaxProc
177 0 : << " " << setw(6) << processes[ip].idProc << "\n";
178 :
179 : // Done.
180 0 : osLHEF << "</init>" << endl;
181 0 : return true;
182 :
183 : }
184 :
185 : //--------------------------------------------------------------------------
186 :
187 : // Write event information to a Les Houches Event File.
188 : // Normal mode is to line up event info in columns, but the non-verbose
189 : // altnernative saves space at the expense of human readability.
190 :
191 : bool LHAup::eventLHEF(bool verbose) {
192 :
193 : // Default verbose option.
194 0 : if (verbose) {
195 :
196 : // Write information on process as such.
197 0 : osLHEF << "<event>\n" << scientific << setprecision(6)
198 0 : << " " << setw(5) << particles.size() - 1
199 0 : << " " << setw(5) << idProc
200 0 : << " " << setw(13) << weightProc
201 0 : << " " << setw(13) << scaleProc
202 0 : << " " << setw(13) << alphaQEDProc
203 0 : << " " << setw(13) << alphaQCDProc << "\n";
204 :
205 : // Write information on the particles, excluding zeroth.
206 0 : for (int ip = 1; ip < int(particles.size()); ++ip) {
207 0 : LHAParticle& ptNow = particles[ip];
208 0 : osLHEF << " " << setw(8) << ptNow.idPart
209 0 : << " " << setw(5) << ptNow.statusPart
210 0 : << " " << setw(5) << ptNow.mother1Part
211 0 : << " " << setw(5) << ptNow.mother2Part
212 0 : << " " << setw(5) << ptNow.col1Part
213 0 : << " " << setw(5) << ptNow.col2Part << setprecision(10)
214 0 : << " " << setw(17) << ptNow.pxPart
215 0 : << " " << setw(17) << ptNow.pyPart
216 0 : << " " << setw(17) << ptNow.pzPart
217 0 : << " " << setw(17) << ptNow.ePart
218 0 : << " " << setw(17) << ptNow.mPart << setprecision(6);
219 0 : if (ptNow.tauPart == 0.) osLHEF << " 0.";
220 0 : else osLHEF << " " << setw(13) << ptNow.tauPart;
221 0 : if (ptNow.spinPart == 9.) osLHEF << " 9.";
222 0 : else osLHEF << " " << setw(13) << ptNow.spinPart;
223 0 : osLHEF << "\n";
224 : }
225 :
226 : // Optionally write information on PDF values at hard interaction.
227 0 : if (pdfIsSetSave) osLHEF << "#pdf"
228 0 : << " " << setw(4) << id1pdfSave
229 0 : << " " << setw(4) << id2pdfSave
230 0 : << " " << setw(13) << x1pdfSave
231 0 : << " " << setw(13) << x2pdfSave
232 0 : << " " << setw(13) << scalePDFSave
233 0 : << " " << setw(13) << pdf1Save
234 0 : << " " << setw(13) << pdf2Save << "\n";
235 :
236 : // Alternative non-verbose option.
237 : } else {
238 :
239 : // Write information on process as such.
240 0 : osLHEF << "<event>\n" << scientific << setprecision(6)
241 0 : << particles.size() - 1 << " " << idProc << " "
242 0 : << weightProc << " " << scaleProc << " "
243 0 : << alphaQEDProc << " " << alphaQCDProc << "\n";
244 :
245 : // Write information on the particles, excluding zeroth.
246 0 : for (int ip = 1; ip < int(particles.size()); ++ip) {
247 0 : LHAParticle& ptNow = particles[ip];
248 0 : osLHEF << ptNow.idPart << " " << ptNow.statusPart
249 0 : << " " << ptNow.mother1Part << " " << ptNow.mother2Part
250 0 : << " " << ptNow.col1Part << " " << ptNow.col2Part
251 0 : << setprecision(10) << " " << ptNow.pxPart
252 0 : << " " << ptNow.pyPart << " " << ptNow.pzPart
253 0 : << " " << ptNow.ePart << " " << ptNow.mPart
254 0 : << setprecision(6);
255 0 : if (ptNow.tauPart == 0.) osLHEF << " 0.";
256 0 : else osLHEF << " " << setw(13) << ptNow.tauPart;
257 0 : if (ptNow.spinPart == 9.) osLHEF << " 9.";
258 0 : else osLHEF << " " << setw(13) << ptNow.spinPart;
259 0 : osLHEF << "\n";
260 : }
261 :
262 : // Optionally write information on PDF values at hard interaction.
263 0 : if (pdfIsSetSave) osLHEF << "#pdf" << " " << id1pdfSave
264 0 : << " " << id2pdfSave << " " << x1pdfSave << " " << x2pdfSave
265 0 : << " " << scalePDFSave << " " << pdf1Save << " " << pdf2Save
266 0 : << "\n";
267 : }
268 :
269 : // Done.
270 0 : osLHEF << "</event>" << endl;
271 0 : return true;
272 :
273 : }
274 :
275 : //--------------------------------------------------------------------------
276 :
277 : // Write end of a Les Houches Event File and close it.
278 :
279 : bool LHAup::closeLHEF(bool updateInit) {
280 :
281 : // Write an end to the file.
282 0 : osLHEF << "</LesHouchesEvents>" << endl;
283 0 : osLHEF.close();
284 :
285 : // Optionally update the cross section information.
286 0 : if (updateInit) {
287 0 : const char* cstring = fileName.c_str();
288 0 : osLHEF.open(cstring, ios::in | ios::out);
289 :
290 : // Rewrite header; identically with what openLHEF did.
291 0 : osLHEF << "<LesHouchesEvents version=\"1.0\">\n"
292 0 : << "<!--\n"
293 0 : << " File written by Pythia8::LHAup on "
294 0 : << dateNow << " at " << timeNow << "\n"
295 0 : << "-->" << endl;
296 :
297 : // Redo initialization information.
298 0 : initLHEF();
299 0 : osLHEF.close();
300 0 : }
301 :
302 : // Done.
303 0 : return true;
304 :
305 : }
306 :
307 : //--------------------------------------------------------------------------
308 :
309 : // Read in initialization information from a Les Houches Event File.
310 :
311 : bool LHAup::setInitLHEF(istream& is, bool readHeaders) {
312 :
313 : // Check that first line is consistent with proper LHEF file.
314 0 : string line;
315 0 : if (!getline(is, line)) return false;
316 0 : if (line.find("<LesHouchesEvents") == string::npos) return false;
317 0 : if (line.find("version=\"1.0\"" ) == string::npos ) return false;
318 :
319 : // What to search for if reading headers; if not reading
320 : // headers then return to default behaviour
321 0 : string headerTag = (readHeaders) ? "<header>" : "<init";
322 :
323 : // Loop over lines until an <init (or optionally <header>) tag
324 : // is found first on a line.
325 0 : string tag = " ";
326 : do {
327 0 : if (!getline(is, line)) return false;
328 0 : if (line.find_first_not_of(" \n\t\v\b\r\f\a") != string::npos) {
329 0 : istringstream getfirst(line);
330 0 : getfirst >> tag;
331 0 : if (!getfirst) return false;
332 0 : }
333 0 : } while (tag != "<init>" && tag != "<init" && tag != headerTag);
334 :
335 : // If header tag found, process if required
336 0 : if (readHeaders == true && tag == headerTag) {
337 : // Temporary local storage
338 0 : map < string, string > headerMap;
339 :
340 : // Loop over lines until an <init> tag is found.
341 : bool read = true, newKey = false;
342 0 : string key = "base";
343 0 : vector < string > keyVec;
344 :
345 0 : while (true) {
346 0 : if (!getline(is, line)) return false;
347 :
348 : // Break lines containing multiple tags into two segments.
349 : // (Could be generalized to multiple segments but this is
350 : // sufficient to handle at least <tag>info</tag> on same line.
351 0 : size_t firstTagEnd = line.find_first_of(">");
352 0 : size_t secondTagBegin = line.find_first_of("<",firstTagEnd);
353 0 : vector<string> lineVec;
354 0 : if (firstTagEnd != string::npos && secondTagBegin != string::npos) {
355 0 : lineVec.push_back(line.substr(0,secondTagBegin));
356 0 : lineVec.push_back(line.substr(secondTagBegin,
357 0 : line.size()-secondTagBegin));
358 0 : }
359 : else {
360 0 : lineVec.push_back(line);
361 : }
362 :
363 : // Loop over segments of current line
364 0 : for (int iVec=0; iVec<int(lineVec.size()); ++iVec) {
365 0 : line = lineVec[iVec];
366 :
367 : // Clean line to contain only valid characters
368 0 : size_t posBeg = line.find_first_not_of(" \n\t\v\b\r\f\a");
369 0 : size_t posEnd = line.find_last_not_of(" \n\t\v\b\r\f\a");
370 0 : string lineClean = " ";
371 0 : if (posBeg != string::npos && posEnd != string::npos && posBeg
372 0 : < posEnd) {
373 0 : lineClean = line.substr(posBeg, posEnd - posBeg + 1);
374 : posBeg = 0;
375 0 : posEnd = lineClean.size();
376 0 : }
377 :
378 : // Check for empty line
379 0 : if (lineClean == " " || posBeg >= posEnd) continue;
380 :
381 : // PZS Jan 2015: Allow multiple open/close tags on a single line.
382 0 : size_t tagBeg = lineClean.find_first_of("<");
383 0 : size_t tagEnd = lineClean.find_first_of(">");
384 :
385 0 : while (tagBeg != string::npos && tagBeg < tagEnd) {
386 :
387 : // Update remainder (non-tag) part of line, for later storage
388 0 : posBeg = tagEnd+1;
389 :
390 : // Only take the first word of the tag,
391 0 : tag = lineClean.substr(tagBeg + 1, tagEnd - tagBeg - 1);
392 0 : istringstream getfirst(tag);
393 0 : getfirst >> tag;
394 :
395 : // Prepare for next while iteration:
396 : // Look for next tag on line and update posBeg and posEnd.
397 0 : tagBeg = lineClean.find_first_of("<",tagEnd);
398 0 : tagEnd = lineClean.find_first_of(">",tagBeg+1);
399 :
400 : // Tag present, so handle here
401 0 : if (getfirst) {
402 :
403 : // Exit condition
404 0 : if (tag == "init") break;
405 :
406 : // End of header block; keep reading until <init> tag,
407 : // but do not store any further information
408 0 : else if (tag == "/header") {
409 : read = false;
410 0 : continue;
411 :
412 : // Opening tag
413 0 : } else if (tag[0] != '/') {
414 0 : keyVec.push_back(tag);
415 : newKey = true;
416 0 : continue;
417 :
418 : // Closing tag that matches current key
419 0 : } else if (tag == "/" + keyVec.back()) {
420 0 : keyVec.pop_back();
421 : newKey = true;
422 0 : continue;
423 :
424 : // Also check for forgotten close tag: next-to-last element
425 0 : } else if (keyVec.size() >= 2
426 0 : && tag == "/" + keyVec[keyVec.size()-2]) {
427 0 : infoPtr->errorMsg("Warning in LHAup::setInitLHEF:"
428 0 : " corrupt LHEF end tag",keyVec.back());
429 0 : keyVec.pop_back();
430 0 : keyVec.pop_back();
431 : newKey = true;
432 0 : continue;
433 : }
434 :
435 : } // if (getfirst)
436 :
437 0 : } // Loop over tags
438 :
439 : // Exit condition
440 0 : if (tag == "init") break;
441 :
442 : // At this point the (rest of) the line is not a tag;
443 : // If no longer reading anything, skip.
444 0 : if (!read) continue;
445 :
446 : // Check for key change
447 0 : if (newKey) {
448 0 : if (keyVec.empty()) key = "base";
449 0 : else key = keyVec[0];
450 0 : for (size_t i = 1; i < keyVec.size(); i++)
451 0 : key += "." + keyVec[i];
452 : newKey = false;
453 0 : }
454 :
455 : // Check if anything remains to store of this line
456 0 : posBeg = line.find_first_not_of(" \n\t\v\b\r\f\a",posBeg);
457 0 : if (posBeg == string::npos || posBeg > posEnd) continue;
458 :
459 : // Append information to local storage
460 0 : headerMap[key] += line.substr(posBeg,posEnd - posBeg + 1) + "\n";
461 :
462 0 : } // Loop over line segments
463 :
464 : // Exit condition
465 0 : if (tag == "init") break;
466 :
467 0 : } // while (true)
468 :
469 : // Copy information to info using LHAup::setInfoHeader
470 0 : for (map < string, string >::iterator it = headerMap.begin();
471 0 : it != headerMap.end(); it++)
472 0 : setInfoHeader(it->first, it->second);
473 :
474 0 : } // if (readHeaders == true && tag == headerTag)
475 :
476 : // Read in first info line; done if empty.
477 0 : if (!getline(is, line)) return false;
478 0 : if (line.find("</init") != string::npos) return true;
479 :
480 : // Read in beam and strategy info, and store it.
481 0 : int idbmupA, idbmupB;
482 0 : double ebmupA, ebmupB;
483 0 : int pdfgupA, pdfgupB, pdfsupA, pdfsupB, idwtup, nprup;
484 0 : istringstream getbms(line);
485 0 : getbms >> idbmupA >> idbmupB >> ebmupA >> ebmupB >> pdfgupA
486 0 : >> pdfgupB >> pdfsupA >> pdfsupB >> idwtup >> nprup;
487 0 : if (!getbms) return false;
488 0 : setBeamA(idbmupA, ebmupA, pdfgupA, pdfsupA);
489 0 : setBeamB(idbmupB, ebmupB, pdfgupB, pdfsupB);
490 0 : setStrategy(idwtup);
491 :
492 : // Read in process info, one process at a time, and store it.
493 0 : double xsecup, xerrup, xmaxup;
494 0 : xSecSumSave = 0.;
495 0 : xErrSumSave = 0.;
496 0 : int lprup;
497 0 : for (int ip = 0; ip < nprup; ++ip) {
498 0 : if (!getline(is, line)) return false;
499 0 : istringstream getpro(line);
500 0 : getpro >> xsecup >> xerrup >> xmaxup >> lprup ;
501 0 : if (!getpro) return false;
502 0 : addProcess(lprup, xsecup, xerrup, xmaxup);
503 0 : xSecSumSave += xsecup;
504 0 : xErrSumSave += pow2(xerrup);
505 0 : }
506 0 : xErrSumSave = sqrt(xErrSumSave);
507 :
508 : // Reading worked.
509 0 : return true;
510 :
511 0 : }
512 :
513 : //--------------------------------------------------------------------------
514 :
515 : // Read in event information from a Les Houches Event File,
516 : // into a staging area where it can be reused by setOldEventLHEF.
517 :
518 : bool LHAup::setNewEventLHEF(istream& is) {
519 :
520 : // Loop over lines until an <event tag is found first on a line.
521 0 : string line, tag;
522 0 : do {
523 0 : if (!getline(is, line)) return false;
524 0 : if (line.find_first_not_of(" \n\t\v\b\r\f\a") != string::npos) {
525 0 : istringstream getfirst(line);
526 0 : getfirst >> tag;
527 0 : if (!getfirst) return false;
528 0 : }
529 0 : } while (tag != "<event>" && tag != "<event");
530 :
531 : // Read in process info and store it.
532 0 : if (!getline(is, line)) return false;
533 0 : istringstream getpro(line);
534 0 : getpro >> nupSave >> idprupSave >> xwgtupSave >> scalupSave
535 0 : >> aqedupSave >> aqcdupSave;
536 0 : if (!getpro) return false;
537 :
538 : // Reset particlesSave vector, add slot-0 empty particle.
539 0 : particlesSave.clear();
540 0 : particlesSave.push_back( LHAParticle() );
541 :
542 : // Read in particle info one by one, and store it.
543 : // Note unusual C++ loop range, to better reflect LHA/Fortran standard.
544 : // (Recall that process(...) above added empty particle at index 0.)
545 0 : int idup, istup, mothup1, mothup2, icolup1, icolup2;
546 0 : double pup1, pup2, pup3, pup4, pup5, vtimup, spinup;
547 0 : for (int ip = 1; ip <= nupSave; ++ip) {
548 0 : if (!getline(is, line)) return false;
549 0 : istringstream getall(line);
550 0 : getall >> idup >> istup >> mothup1 >> mothup2 >> icolup1 >> icolup2
551 0 : >> pup1 >> pup2 >> pup3 >> pup4 >> pup5 >> vtimup >> spinup;
552 0 : if (!getall) return false;
553 0 : particlesSave.push_back( LHAParticle( idup, istup, mothup1, mothup2,
554 0 : icolup1, icolup2, pup1, pup2, pup3, pup4, pup5, vtimup, spinup, -1.) );
555 0 : }
556 :
557 : // Flavour and x values of hard-process initiators.
558 0 : id1InSave = particlesSave[1].idPart;
559 0 : id2InSave = particlesSave[2].idPart;
560 0 : x1InSave = (eBeamASave > 0.) ? particlesSave[1].ePart / eBeamASave : 0.;
561 0 : x2InSave = (eBeamBSave > 0.) ? particlesSave[2].ePart / eBeamBSave : 0.;
562 :
563 : // Continue parsing till </event>. Look for optional info on the way.
564 0 : getPDFSave = false;
565 0 : getScale = false;
566 0 : do {
567 0 : if (!getline(is, line)) return false;
568 0 : istringstream getinfo(line);
569 0 : getinfo >> tag;
570 0 : if (!getinfo) return false;
571 :
572 : // Extract PDF info if present.
573 0 : if (tag == "#pdf" && !getPDFSave) {
574 0 : getinfo >> id1pdfInSave >> id2pdfInSave >> x1pdfInSave >> x2pdfInSave
575 0 : >> scalePDFInSave >> pdf1InSave >> pdf2InSave;
576 0 : if (!getinfo) return false;
577 0 : getPDFSave = true;
578 :
579 : // Extract scale info if present.
580 0 : } else if (tag == "#" && !getScale) {
581 0 : double scaleIn = 0;
582 0 : for (int i = 3; i < int(particlesSave.size()); ++i)
583 0 : if (particlesSave[i].statusPart == 1) {
584 0 : if ( !(getinfo >> scaleIn) ) return false;
585 0 : particlesSave[i].scalePart = scaleIn;
586 0 : }
587 0 : if (!getinfo) return false;
588 0 : getScale = true;
589 0 : }
590 0 : } while (tag != "</event>" && tag != "</event");
591 :
592 : // Need id and x values even when no PDF info. Rest empty.
593 0 : if (!getPDFSave) {
594 0 : id1pdfInSave = id1InSave;
595 0 : id2pdfInSave = id2InSave;
596 0 : x1pdfInSave = x1InSave;
597 0 : x2pdfInSave = x2InSave;
598 0 : scalePDFInSave = 0.;
599 0 : pdf1InSave = 0.;
600 0 : pdf2InSave = 0.;
601 0 : }
602 :
603 : // Reading worked.
604 0 : return true;
605 :
606 0 : }
607 :
608 : //--------------------------------------------------------------------------
609 :
610 : // Make current event information read in by setNewEventLHEF.
611 :
612 : bool LHAup::setOldEventLHEF() {
613 :
614 : // Store saved event, optionally also parton density information.
615 0 : setProcess( idprupSave, xwgtupSave, scalupSave, aqedupSave, aqcdupSave);
616 0 : for (int ip = 1; ip <= nupSave; ++ip) addParticle( particlesSave[ip] );
617 0 : setIdX( id1InSave, id2InSave, x1InSave, x2InSave);
618 0 : setPdf( id1pdfInSave, id2pdfInSave, x1pdfInSave, x2pdfInSave,
619 0 : scalePDFInSave, pdf1InSave, pdf2InSave, getPDFSave);
620 :
621 : // Done;
622 0 : return true;
623 :
624 : }
625 :
626 : //--------------------------------------------------------------------------
627 :
628 : // Open a file using provided ifstream and return a pointer to an istream
629 : // that can be used to process the file.
630 :
631 : istream* LHAup::openFile(const char *fn, ifstream &ifs) {
632 : // Open the file
633 0 : ifs.open(fn);
634 0 : return (istream *) &ifs;
635 : }
636 :
637 : //--------------------------------------------------------------------------
638 :
639 : // Companion method to 'openFile', above.
640 : // Correctly deallocates memory if required before closing the file.
641 :
642 : void LHAup::closeFile(istream *&is, ifstream &ifs) {
643 : // If the istream pointer is not NULL and is not the
644 : // same as the ifstream, then delete pointer.
645 0 : if (is && is != &ifs) delete is;
646 0 : is = NULL;
647 :
648 : // Close the file
649 0 : if (ifs.is_open()) ifs.close();
650 0 : }
651 :
652 : //==========================================================================
653 :
654 : // LHAupLHEF class.
655 :
656 : //--------------------------------------------------------------------------
657 :
658 : // Routine for doing the job of reading and setting initialization info.
659 :
660 : bool LHAupLHEF::setInitLHEF( istream & isIn, bool readHead ) {
661 :
662 : // Done if there was a problem with initialising the reader
663 0 : if (!reader.isGood) return false;
664 :
665 : // Construct header information (stored in comments strings or optional
666 : // header file), so that reading of headers is possible.
667 0 : string comments;
668 0 : comments+="<LesHouchesEvents version =\"3.0\">\n";
669 0 : comments+="<header>\n";
670 0 : comments+=reader.headerComments;
671 0 : comments+="</header>\n";
672 0 : comments+="<init>\n";
673 0 : comments+=reader.initComments;
674 0 : comments+="</init>\n";
675 0 : istringstream is1(comments);
676 0 : bool useComments = (headerfile == NULL);
677 0 : istream & iss((useComments ? is1 : isIn));
678 :
679 : // Check that first line is consistent with proper LHEF file.
680 0 : string line;
681 0 : if ( useComments && !getline(iss,line)) return false;
682 0 : if (!useComments && !getLine(line)) return false;
683 :
684 : // What to search for if reading headers; if not reading
685 : // headers then return to default behaviour
686 0 : string headerTag = (readHead) ? "<header>" : "<init";
687 :
688 : // Loop over lines until an <init (or optionally <header>) tag
689 : // is found first on a line.
690 0 : string tag = " ";
691 : do {
692 0 : if ( useComments && !getline(iss,line)) return false;
693 0 : if (!useComments && !getLine(line)) return false;
694 0 : if (line.find_first_not_of(" \n\t\v\b\r\f\a") != string::npos) {
695 0 : istringstream getfirst(line);
696 0 : getfirst >> tag;
697 0 : if (!getfirst) return false;
698 0 : }
699 0 : } while (tag != "<init>" && tag != "<init" && tag != headerTag);
700 :
701 : // If header tag found, process if required
702 0 : if (readHead == true && tag == headerTag) {
703 : // Temporary local storage
704 0 : map < string, string > headerMap;
705 :
706 : // Loop over lines until an <init> tag is found.
707 : bool read = true, newKey = false;
708 : int commentDepth = 0;
709 0 : string key = "base";
710 0 : vector < string > keyVec;
711 0 : while (true) {
712 0 : if ( useComments && !getline(iss,line)) return false;
713 0 : if (!useComments && !getLine(line)) return false;
714 :
715 : // Tell XML parser to ignore comment and CDATA blocks
716 : // If we are currently inside a comment block, check for block end
717 0 : if (commentDepth >= 1 && line.find("-->") != string::npos)
718 0 : commentDepth--;
719 0 : if (commentDepth >= 1 && line.find("]]>") != string::npos)
720 0 : commentDepth--;
721 : // If the comment block did not end on this line, skip to next line
722 0 : if (commentDepth >= 1) continue;
723 : // Check for beginning of comment blocks (parse until comment begins)
724 0 : if (line.find("<!--") != string::npos) {
725 0 : if (line.find("-->") == string::npos) commentDepth++;
726 0 : int comBeg = line.find("<!--");
727 0 : line = line.substr(0,comBeg);
728 0 : }
729 : // Check for beginning of CDATA statement (parse until CDATA begins)
730 0 : if (line.find("<![cdata[") != string::npos
731 0 : || line.find("<![CDATA[") != string::npos) {
732 0 : if (line.find("]]>") == string::npos) commentDepth++;
733 0 : int comBeg = line.find("<![");
734 0 : line = line.substr(0,comBeg);
735 0 : }
736 :
737 : // Break lines containing multiple tags into two segments.
738 : // (Could be generalized to multiple segments but this is
739 : // sufficient to handle at least <tag>info</tag> on same line.
740 0 : size_t firstTagEnd = line.find_first_of(">");
741 0 : size_t secondTagBegin = line.find_first_of("<",firstTagEnd);
742 0 : vector<string> lineVec;
743 0 : if (firstTagEnd != string::npos && secondTagBegin != string::npos) {
744 0 : lineVec.push_back(line.substr(0,secondTagBegin));
745 0 : lineVec.push_back(line.substr(secondTagBegin,
746 0 : line.size()-secondTagBegin));
747 0 : }
748 : else {
749 0 : lineVec.push_back(line);
750 : }
751 :
752 : // Loop over segments of current line
753 0 : for (int iVec=0; iVec<int(lineVec.size()); ++iVec) {
754 0 : line = lineVec[iVec];
755 :
756 : // Clean line to contain only valid characters
757 0 : size_t posBeg = line.find_first_not_of(" \n\t\v\b\r\f\a");
758 0 : size_t posEnd = line.find_last_not_of(" \n\t\v\b\r\f\a");
759 0 : string lineClean = " ";
760 0 : if (posBeg != string::npos && posEnd != string::npos && posBeg
761 0 : < posEnd) {
762 0 : lineClean = line.substr(posBeg, posEnd - posBeg + 1);
763 : posBeg = 0;
764 0 : posEnd = lineClean.size();
765 0 : }
766 :
767 : // Check for empty line
768 0 : if (lineClean == " " || posBeg >= posEnd) continue;
769 :
770 : // PZS Jan 2015: Allow multiple open/close tags on a single line.
771 0 : size_t tagBeg = lineClean.find_first_of("<");
772 0 : size_t tagEnd = lineClean.find_first_of(">");
773 :
774 0 : while (tagBeg != string::npos && tagBeg < tagEnd) {
775 :
776 : // Update remainder (non-tag) part of line, for later storage
777 0 : posBeg = tagEnd+1;
778 :
779 : // Only take the first word of the tag,
780 0 : tag = lineClean.substr(tagBeg + 1, tagEnd - tagBeg - 1);
781 0 : istringstream getfirst(tag);
782 0 : getfirst >> tag;
783 :
784 : // Prepare for next while iteration:
785 : // Look for next tag on line and update posBeg and posEnd.
786 0 : tagBeg = lineClean.find_first_of("<",tagEnd);
787 0 : tagEnd = lineClean.find_first_of(">",tagBeg+1);
788 :
789 : // Tag present, so handle here
790 0 : if (getfirst) {
791 :
792 : // Exit condition
793 0 : if (tag == "init") break;
794 :
795 : // End of header block; keep reading until <init> tag,
796 : // but do not store any further information
797 0 : else if (tag == "/header") {
798 : read = false;
799 0 : continue;
800 :
801 : // Opening tag
802 0 : } else if (tag[0] != '/') {
803 0 : keyVec.push_back(tag);
804 : newKey = true;
805 0 : continue;
806 :
807 : // Closing tag that matches current key
808 0 : } else if (tag == "/" + keyVec.back()) {
809 0 : keyVec.pop_back();
810 : newKey = true;
811 0 : continue;
812 :
813 : // Also check for forgotten close tag: next-to-last element
814 0 : } else if (keyVec.size() >= 2
815 0 : && tag == "/" + keyVec[keyVec.size()-2]) {
816 0 : infoPtr->errorMsg("Warning in LHAupLHEF::setInitLHEF:"
817 0 : " corrupt LHEF end tag",keyVec.back());
818 0 : keyVec.pop_back();
819 0 : keyVec.pop_back();
820 : newKey = true;
821 0 : continue;
822 : }
823 :
824 : } // if (getfirst)
825 :
826 0 : } // Loop over tags
827 :
828 : // Exit condition
829 0 : if (tag == "init") break;
830 :
831 : // At this point the (rest of) the line is not a tag;
832 : // If no longer reading anything, skip.
833 0 : if (!read) continue;
834 :
835 : // Check for key change
836 0 : if (newKey) {
837 0 : if (keyVec.empty()) key = "base";
838 0 : else key = keyVec[0];
839 0 : for (size_t i = 1; i < keyVec.size(); i++)
840 0 : key += "." + keyVec[i];
841 : newKey = false;
842 0 : }
843 :
844 : // Check if anything remains to store of this line
845 0 : posBeg = line.find_first_not_of(" \n\t\v\b\r\f\a",posBeg);
846 0 : if (posBeg == string::npos || posBeg > posEnd) continue;
847 :
848 : // Append information to local storage
849 0 : headerMap[key] += line.substr(posBeg,posEnd - posBeg + 1) + "\n";
850 :
851 0 : } // Loop over line segments
852 :
853 : // Exit condition
854 0 : if (tag == "init") break;
855 :
856 0 : } // while (true)
857 :
858 : // Copy information to info using LHAup::setInfoHeader
859 0 : for (map < string, string >::iterator it = headerMap.begin();
860 0 : it != headerMap.end(); it++)
861 0 : setInfoHeader(it->first, it->second);
862 :
863 0 : } // if (readHead == true && tag == headerTag)
864 :
865 : // Extract beam and strategy info, and store it.
866 : int idbmupA, idbmupB;
867 : double ebmupA, ebmupB;
868 : int pdfgupA, pdfgupB, pdfsupA, pdfsupB, idwtup, nprup;
869 :
870 0 : idbmupA = reader.heprup.IDBMUP.first;
871 0 : idbmupB = reader.heprup.IDBMUP.second;
872 0 : ebmupA = reader.heprup.EBMUP.first;
873 0 : ebmupB = reader.heprup.EBMUP.second;
874 0 : pdfgupA = reader.heprup.PDFGUP.first;
875 : pdfgupB = reader.heprup.PDFGUP.first;
876 0 : pdfsupA = reader.heprup.PDFSUP.first;
877 0 : pdfsupB = reader.heprup.PDFSUP.second;
878 0 : idwtup = reader.heprup.IDWTUP;
879 0 : nprup = reader.heprup.NPRUP;
880 :
881 0 : setBeamA(idbmupA, ebmupA, pdfgupA, pdfsupA);
882 0 : setBeamB(idbmupB, ebmupB, pdfgupB, pdfsupB);
883 0 : setStrategy(idwtup);
884 :
885 : // Read in process info, one process at a time, and store it.
886 0 : double xsecup, xerrup, xmaxup;
887 0 : xSecSumSave = 0.;
888 0 : xErrSumSave = 0.;
889 : int lprup;
890 0 : infoPtr->sigmaLHEFSave.resize(0);
891 0 : for (int ip = 0; ip < nprup; ++ip) {
892 0 : xsecup = reader.heprup.XSECUP[ip];
893 0 : xerrup = reader.heprup.XERRUP[ip];
894 0 : xmaxup = reader.heprup.XMAXUP[ip];
895 0 : lprup = reader.heprup.LPRUP[ip];
896 0 : addProcess(lprup, xsecup, xerrup, xmaxup);
897 0 : xSecSumSave += xsecup;
898 0 : xErrSumSave += pow(xerrup,2);
899 0 : infoPtr->sigmaLHEFSave.push_back(xsecup);
900 : }
901 0 : xErrSumSave = sqrt(xErrSumSave);
902 :
903 : // Now extract LHEF 2.0/3.0 novelties.
904 0 : infoPtr->setLHEF3InitInfo();
905 0 : if (reader.version > 1) {
906 0 : infoPtr->setLHEF3InitInfo( reader.version,
907 0 : &reader.heprup.initrwgt, &(reader.heprup.generators),
908 0 : &(reader.heprup.weightgroups), &(reader.heprup.weights));
909 : }
910 :
911 : // Reading worked.
912 : return true;
913 0 : }
914 :
915 : //--------------------------------------------------------------------------
916 :
917 : // Routine for doing the job of reading and setting info on next event.
918 :
919 : bool LHAupLHEF::setNewEventLHEF() {
920 :
921 : // Done if the reader finished preemptively.
922 0 : if(!reader.readEvent()) return false;
923 :
924 : // Extract process info and store it.
925 0 : nupSave = reader.hepeup.NUP;
926 0 : idprupSave = reader.hepeup.IDPRUP;
927 0 : xwgtupSave = reader.hepeup.XWGTUP;
928 0 : scalupSave = reader.hepeup.SCALUP;
929 0 : aqedupSave = reader.hepeup.AQEDUP;
930 0 : aqcdupSave = reader.hepeup.AQCDUP;
931 :
932 : // Reset particlesSave vector, add slot-0 empty particle.
933 0 : particlesSave.clear();
934 0 : particlesSave.push_back( Pythia8::LHAParticle() );
935 :
936 : // Read in particle info one by one, and store it.
937 : // Note unusual C++ loop range, to better reflect LHA/Fortran standard.
938 : // (Recall that process(...) above added empty particle at index 0.)
939 : int idup, istup, mothup1, mothup2, icolup1, icolup2;
940 : double pup1, pup2, pup3, pup4, pup5, vtimup, spinup;
941 0 : for ( int i = 0; i < reader.hepeup.NUP; ++i ) {
942 : // Extract information stored in reader.
943 0 : idup = reader.hepeup.IDUP[i];
944 0 : istup = reader.hepeup.ISTUP[i];
945 0 : mothup1 = reader.hepeup.MOTHUP[i].first;
946 0 : mothup2 = reader.hepeup.MOTHUP[i].second;
947 0 : icolup1 = reader.hepeup.ICOLUP[i].first;
948 0 : icolup2 = reader.hepeup.ICOLUP[i].second;
949 0 : pup1 = reader.hepeup.PUP[i][0];
950 0 : pup2 = reader.hepeup.PUP[i][1];
951 0 : pup3 = reader.hepeup.PUP[i][2];
952 0 : pup4 = reader.hepeup.PUP[i][3];
953 0 : pup5 = reader.hepeup.PUP[i][4];
954 0 : vtimup = reader.hepeup.VTIMUP[i];
955 0 : spinup = reader.hepeup.SPINUP[i];
956 0 : particlesSave.push_back( Pythia8::LHAParticle( idup,istup,mothup1,mothup2,
957 : icolup1, icolup2, pup1, pup2, pup3, pup4, pup5, vtimup, spinup, -1.) );
958 : }
959 :
960 : // Flavour and x values of hard-process initiators.
961 0 : id1InSave = particlesSave[1].idPart;
962 0 : id2InSave = particlesSave[2].idPart;
963 0 : x1InSave = (eBeamA() > 0.) ? particlesSave[1].ePart / eBeamA() : 0.;
964 0 : x2InSave = (eBeamB() > 0.) ? particlesSave[2].ePart / eBeamB() : 0.;
965 :
966 : // Parse event comments and look for optional info on the way.
967 0 : std::string line, tag;
968 0 : std::stringstream ss(reader.eventComments);
969 0 : getPDFSave = false;
970 0 : getScale = false;
971 0 : getScale = (setScalesFromLHEF && reader.version == 1) ? false : true;
972 0 : while (getline(ss, line)) {
973 0 : istringstream getinfo(line);
974 0 : getinfo >> tag;
975 0 : if (!getinfo) break;
976 : // Extract PDF info if present.
977 0 : if (tag == "#pdf" && !getPDFSave) {
978 0 : getinfo >> id1pdfInSave >> id2pdfInSave >> x1pdfInSave >> x2pdfInSave
979 0 : >> scalePDFInSave >> pdf1InSave >> pdf2InSave;
980 0 : if (!getinfo) return false;
981 0 : getPDFSave = true;
982 : // Extract scale info if present.
983 0 : } else if (tag == "#" && !getScale) {
984 0 : double scaleIn = 0;
985 0 : for (int i = 3; i < int(particlesSave.size()); ++i)
986 0 : if (particlesSave[i].statusPart == 1) {
987 0 : if ( !(getinfo >> scaleIn) ) return false;
988 0 : particlesSave[i].scalePart = scaleIn;
989 0 : }
990 0 : if (!getinfo) return false;
991 0 : getScale = true;
992 0 : }
993 0 : }
994 :
995 : // Set production scales from <scales> tag.
996 0 : if ( setScalesFromLHEF && reader.version > 1 ){
997 0 : for ( map<string,double>::const_iterator
998 0 : it = reader.hepeup.scales.attributes.begin();
999 0 : it != reader.hepeup.scales.attributes.end(); ++it ) {
1000 0 : if ( it->first.find_last_of("_") != string::npos) {
1001 0 : unsigned iFound = it->first.find_last_of("_") + 1;
1002 0 : int iPos = atoi(it->first.substr(iFound).c_str());
1003 : // Only set production scales of final particles.
1004 0 : if ( iPos < int(particlesSave.size())
1005 0 : && particlesSave[iPos].statusPart == 1)
1006 0 : particlesSave[iPos].scalePart = it->second;
1007 0 : }
1008 : }
1009 0 : }
1010 :
1011 : // Need id and x values even when no PDF info. Rest empty.
1012 0 : if (!getPDFSave) {
1013 0 : id1pdfInSave = id1InSave;
1014 0 : id2pdfInSave = id2InSave;
1015 0 : x1pdfInSave = x1InSave;
1016 0 : x2pdfInSave = x2InSave;
1017 0 : scalePDFInSave = 0.;
1018 0 : pdf1InSave = 0.;
1019 0 : pdf2InSave = 0.;
1020 0 : }
1021 :
1022 : // Now extract LHEF 2.0/3.0 novelties.
1023 0 : infoPtr->setLHEF3EventInfo();
1024 : // Set everything for 2.0 and 3.0
1025 0 : if (reader.version > 1) {
1026 0 : infoPtr->setLHEF3EventInfo( &reader.hepeup.attributes,
1027 0 : &reader.hepeup.weights_detailed, &reader.hepeup.weights_compressed,
1028 0 : &reader.hepeup.scales, &reader.hepeup.weights, &reader.hepeup.rwgt);
1029 : // Try to at least set the event attributes for 1.0
1030 : } else {
1031 0 : infoPtr->setLHEF3EventInfo( &reader.hepeup.attributes, 0, 0, 0, 0, 0);
1032 : }
1033 :
1034 : // Reading worked.
1035 0 : return true;
1036 :
1037 0 : }
1038 :
1039 : //==========================================================================
1040 :
1041 : // LHAupFromPYTHIA8 class.
1042 :
1043 : //--------------------------------------------------------------------------
1044 :
1045 : // Read in initialization information from PYTHIA 8.
1046 :
1047 : bool LHAupFromPYTHIA8::setInit() {
1048 :
1049 : // Read in beam from Info class. Parton density left empty.
1050 0 : int idbmupA = infoPtr->idA();
1051 0 : int idbmupB = infoPtr->idB();
1052 0 : double ebmupA = infoPtr->eA();
1053 0 : double ebmupB = infoPtr->eB();
1054 : int pdfgupA = 0;
1055 : int pdfgupB = 0;
1056 : int pdfsupA = 0;
1057 : int pdfsupB = 0;
1058 0 : setBeamA(idbmupA, ebmupA, pdfgupA, pdfsupA);
1059 0 : setBeamB(idbmupB, ebmupB, pdfgupB, pdfsupB);
1060 :
1061 : // Currently only one allowed strategy.
1062 : int idwtup = 3;
1063 0 : setStrategy(idwtup);
1064 :
1065 : // Only one process with dummy information. (Can overwrite at the end.)
1066 : int lprup = 9999;
1067 : double xsecup = 1.;
1068 : double xerrup = 0.;
1069 : double xmaxup = 1.;
1070 0 : addProcess(lprup, xsecup, xerrup, xmaxup);
1071 :
1072 : // Done.
1073 0 : return true;
1074 :
1075 : }
1076 :
1077 : //--------------------------------------------------------------------------
1078 :
1079 : // Read in event information from PYTHIA 8.
1080 :
1081 : bool LHAupFromPYTHIA8::setEvent( int) {
1082 :
1083 : // Read process information from Info class, and store it.
1084 : // Note: renormalization scale here, factorization further down.
1085 : // For now always convert to process 9999, instead of infoPtr->code().
1086 : int idprup = 9999;
1087 0 : double xwgtup = infoPtr->weight();
1088 0 : double scalup = infoPtr->QRen();
1089 0 : double aqedup = infoPtr->alphaEM();
1090 0 : double aqcdup = infoPtr->alphaS();
1091 0 : setProcess(idprup, xwgtup, scalup, aqedup, aqcdup);
1092 :
1093 : // Read in particle info one by one, excluding zero and beams, and store it.
1094 : // Note unusual C++ loop range, to better reflect LHA/Fortran standard.
1095 0 : int nup = processPtr->size() - 3;
1096 : int idup, statusup, istup, mothup1, mothup2, icolup1, icolup2;
1097 : double pup1, pup2, pup3, pup4, pup5, vtimup, spinup;
1098 0 : for (int ip = 1; ip <= nup; ++ip) {
1099 0 : Particle& particle = (*processPtr)[ip + 2];
1100 0 : idup = particle.id();
1101 : // Convert from PYTHIA8 to LHA status codes.
1102 0 : statusup = particle.status();
1103 0 : if (ip < 3) istup = -1;
1104 0 : else if (statusup < 0) istup = 2;
1105 : else istup = 1;
1106 0 : mothup1 = max(0, particle.mother1() - 2);
1107 0 : mothup2 = max(0, particle.mother2() - 2);
1108 0 : icolup1 = particle.col();
1109 0 : icolup2 = particle.acol();
1110 0 : pup1 = particle.px();
1111 0 : pup2 = particle.py();
1112 0 : pup3 = particle.pz();
1113 0 : pup4 = particle.e();
1114 0 : pup5 = particle.m();
1115 0 : vtimup = particle.tau();
1116 0 : spinup = particle.pol();
1117 0 : addParticle(idup, istup, mothup1, mothup2, icolup1, icolup2,
1118 : pup1, pup2, pup3, pup4, pup5, vtimup, spinup, -1.) ;
1119 : }
1120 :
1121 : // Extract hard-process initiator information from Info class, and store it.
1122 0 : int id1up = infoPtr->id1();
1123 0 : int id2up = infoPtr->id2();
1124 0 : double x1up = infoPtr->x1();
1125 0 : double x2up = infoPtr->x2();
1126 0 : setIdX( id1up, id2up, x1up, x2up);
1127 :
1128 : // Also extract pdf information from Info class, and store it.
1129 0 : int id1pdfup = infoPtr->id1pdf();
1130 0 : int id2pdfup = infoPtr->id2pdf();
1131 0 : double x1pdfup = infoPtr->x1pdf();
1132 0 : double x2pdfup = infoPtr->x2pdf();
1133 0 : double scalePDFup = infoPtr->QFac();
1134 0 : double pdf1up = infoPtr->pdf1();
1135 0 : double pdf2up = infoPtr->pdf2();
1136 0 : setPdf(id1pdfup, id2pdfup, x1pdfup, x2pdfup, scalePDFup, pdf1up,
1137 : pdf2up, true);
1138 :
1139 : // Done.
1140 0 : return true;
1141 :
1142 : }
1143 :
1144 : //--------------------------------------------------------------------------
1145 :
1146 : // Update cross-section information at the end of the run.
1147 :
1148 : bool LHAupFromPYTHIA8::updateSigma() {
1149 :
1150 : // Read out information from PYTHIA 8 and send it in to LHA.
1151 0 : double sigGen = CONVERTMB2PB * infoPtr->sigmaGen();
1152 0 : double sigErr = CONVERTMB2PB * infoPtr->sigmaErr();
1153 0 : setXSec(0, sigGen);
1154 0 : setXErr(0, sigErr);
1155 :
1156 : // Done.
1157 0 : return true;
1158 :
1159 : }
1160 :
1161 : //==========================================================================
1162 :
1163 : } // end namespace Pythia8
|
