CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
AlpgenSource.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <iomanip>
3 #include <sstream>
4 #include <fstream>
5 #include <string>
6 #include <memory>
7 #include <cmath>
8 
9 #include <boost/algorithm/string/classification.hpp>
10 #include <boost/algorithm/string/split.hpp>
11 
20 
25 
28 
30 public:
32  AlpgenSource(const edm::ParameterSet &params,
33  const edm::InputSourceDescription &desc);
34 
36  virtual ~AlpgenSource();
37 
38 private:
40  virtual void produce(edm::Event &event) override;
41  virtual void beginRun(edm::Run &run) override;
42 
44  template<typename T>
47  template<typename T>
48  T getParameter(AlpgenHeader::Parameter index, const T &defValue) const;
49 
53  const std::string &comment) const;
54 
58  unsigned int processID() const;
59 
61  bool readAlpgenEvent(lhef::HEPEUP &hepeup);
62 
65 
67  std::unique_ptr<std::ifstream> inputFile_;
68 
70  unsigned long skipEvents_;
71 
73  unsigned long nEvent_;
74 
79 
80  std::unique_ptr<lhef::HEPEUP> hepeup_;
81 
84 
87 
92 };
93 
95  const edm::InputSourceDescription &desc) :
96  edm::ProducerSourceFromFiles(params, desc, false),
97  skipEvents_(params.getUntrackedParameter<unsigned int>("skipEvents", 0)),
98  nEvent_(0), lheAlpgenUnwParHeader("AlpgenUnwParFile"),
99  extraHeaderFileName_(params.getUntrackedParameter<std::string>("extraHeaderFileName","")),
100  extraHeaderName_(params.getUntrackedParameter<std::string>("extraHeaderName","")),
101  writeAlpgenWgtFile(params.getUntrackedParameter<bool>("writeAlpgenWgtFile", true)),
102  writeAlpgenParFile(params.getUntrackedParameter<bool>("writeAlpgenParFile", true)),
103  writeExtraHeader(params.getUntrackedParameter<bool>("writeExtraHeader", false))
104 {
105  std::vector<std::string> allFileNames = fileNames();
106 
107  // Only one filename
108  if (allFileNames.size() != 1)
109  throw cms::Exception("Generator|AlpgenInterface")
110  << "AlpgenSource needs exactly one file specified "
111  "for now." << std::endl;
112 
113  fileName_ = allFileNames[0];
114 
115  // Strip the "file:" prefix
116  if (fileName_.find("file:") != 0)
117  throw cms::Exception("Generator|AlpgenInterface") << "AlpgenSource only supports the file: scheme "
118  "for now." << std::endl;
119  fileName_.erase(0, 5);
120 
121  // Open the _unw.par file to store additional
122  // informations in the LHERunInfoProduct
123  std::ifstream reader((fileName_ + "_unw.par").c_str());
124  if (!reader.good())
125  throw cms::Exception("Generator|AlpgenInterface")
126  << "AlpgenSource was unable to open the file \""
127  << fileName_ << "_unw.par\"." << std::endl;
128 
129  // A full copy of the _unw.par file in the LHE header.
130  char buffer[256];
131  while(reader.getline(buffer, sizeof buffer))
132  lheAlpgenUnwParHeader.addLine(std::string(buffer) + "\n");
133 
134  // Parse that header to setup an Alpgen header,
135  // which will be used in the production itself.
138  throw cms::Exception("Generator|AlpgenInterface")
139  << "AlpgenSource was unable to parse the Alpgen "
140  << "unweighted parameter file." << std::endl;
141 
142  // Declare the products.
143  produces<LHERunInfoProduct, edm::InRun>();
144  produces<LHEEventProduct>();
145 }
146 
148 {
149 }
150 
152  const std::string &comment) const
153 {
154  std::ostringstream ss;
155  ss << std::setw(9) << pdgId << " " << std::scientific
156  << std::setprecision(9) << header.masses[mass] << " # "
157  << comment << std::endl;
158  return ss.str();
159 }
160 
162 {
163  // At this point, the lheUnwParHeader has the full contents of the _unw.par
164  // file. So we can get the HEPRUP information from the LHE header itself.
165  lhef::HEPRUP heprup;
166 
167  // Get basic run information.
168  // Beam identity.
169  heprup.IDBMUP.first = 2212;
170  switch(getParameter<int>(AlpgenHeader::ih2)) {
171  case 1:
172  heprup.IDBMUP.second = 2212;
173  break;
174  case -1:
175  heprup.IDBMUP.second = -2212;
176  break;
177  default:
178  throw cms::Exception("Generator|AlpgenInterface")
179  << "AlpgenSource was unable to understand the ih2 "
180  << "parameter." << std::endl;
181  }
182 
183  // Beam energy.
184  heprup.EBMUP.second = heprup.EBMUP.first =
185  getParameter<double>(AlpgenHeader::ebeam);
186 
187  // PDF info. Initially, Alpgen doesn't fill it.
188  heprup.PDFGUP.first = -1;
189  heprup.PDFGUP.second = -1;
190  heprup.PDFSUP.first = -1;
191  heprup.PDFSUP.second = -1;
192 
193  // Unweighted events.
194  heprup.IDWTUP = 3;
195 
196  // Only one process.
197  heprup.resize(1);
198 
199  // Cross section and error.
200  heprup.XSECUP[0] = header.xsec;
201  heprup.XERRUP[0] = header.xsecErr;
202 
203  // Maximum weight.
204  heprup.XMAXUP[0] = header.xsec;
205 
206  // Process code for Pythia.
207  heprup.LPRUP[0] = processID();
208 
209  // Comments on top.
211  comments.addLine("\n");
212  comments.addLine("\tExtracted by AlpgenInterface\n");
213 
214  // Add SLHA header containing particle masses from Alpgen.
215  // Pythia6Hadronisation will feed the masses to Pythia automatically.
216  LHERunInfoProduct::Header slha("slha");
217  slha.addLine("\n# SLHA header containing masses from Alpgen\n");
218  slha.addLine("Block MASS # Mass spectrum (kinematic masses)\n");
219  slha.addLine("# PDG Mass\n");
220  slha.addLine(slhaMassLine(4, AlpgenHeader::mc, "charm pole mass"));
221  slha.addLine(slhaMassLine(5, AlpgenHeader::mb, "bottom pole mass"));
222  slha.addLine(slhaMassLine(6, AlpgenHeader::mt, "top pole mass"));
223  slha.addLine(slhaMassLine(23, AlpgenHeader::mz, "Z mass"));
224  slha.addLine(slhaMassLine(24, AlpgenHeader::mw, "W mass"));
225  slha.addLine(slhaMassLine(25, AlpgenHeader::mh, "H mass"));
226 
227  char buffer[512];
228 
229  // We also add the information on weighted events.
230  LHERunInfoProduct::Header lheAlpgenWgtHeader("AlpgenWgtFile");
231  if (writeAlpgenWgtFile) {
232  std::ifstream wgtascii((fileName_+".wgt").c_str());
233  while(wgtascii.getline(buffer,512)) {
234  lheAlpgenWgtHeader.addLine(std::string(buffer) + "\n");
235  }
236  }
237 
238  LHERunInfoProduct::Header lheAlpgenParHeader("AlpgenParFile");
239  if (writeAlpgenParFile) {
240  std::ifstream parascii((fileName_+".par").c_str());
241  while(parascii.getline(buffer,512)) {
242  lheAlpgenParHeader.addLine(std::string(buffer) + "\n");
243  }
244  }
245 
246  // If requested by the user, we also add any specific header provided.
247  // Nota bene: the header is put in the LHERunInfoProduct AS IT WAS GIVEN.
248  // That means NO CROSS-CHECKS WHATSOEVER. Use with care.
249  LHERunInfoProduct::Header extraHeader(extraHeaderName_.c_str());
250  if(writeExtraHeader) {
251  std::ifstream extraascii(extraHeaderFileName_.c_str());
252  while(extraascii.getline(buffer,512)) {
253  extraHeader.addLine(std::string(buffer) + "\n");
254  }
255  }
256 
257  // Build the final Run info object. Backwards-compatible order.
258  std::auto_ptr<LHERunInfoProduct> runInfo(new LHERunInfoProduct(heprup));
259  runInfo->addHeader(comments);
260  runInfo->addHeader(lheAlpgenUnwParHeader);
261  if (writeAlpgenWgtFile)
262  runInfo->addHeader(lheAlpgenWgtHeader);
263  if (writeAlpgenParFile)
264  runInfo->addHeader(lheAlpgenParHeader);
265  runInfo->addHeader(slha);
266  if(writeExtraHeader)
267  runInfo->addHeader(extraHeader);
268  run.put(runInfo);
269 
270  // Open the .unw file in the heap, and set the global pointer to it.
271  inputFile_.reset(new std::ifstream((fileName_ + ".unw").c_str()));
272  if (!inputFile_->good())
273  throw cms::Exception("Generator|AlpgenInterface")
274  << "AlpgenSource was unable to open the file \""
275  << fileName_ << ".unw\"." << std::endl;
276 
277 }
278 
279 template<typename T>
281 {
282  std::map<AlpgenHeader::Parameter, double>::const_iterator pos =
283  header.params.find(index);
284  if (pos == header.params.end())
285  throw cms::Exception("Generator|AlpgenInterface")
286  << "Requested Alpgen parameter \""
287  << AlpgenHeader::parameterName(index) << "\" "
288  "not found in Alpgen parameter file." << std::endl;
289 
290  return T(pos->second);
291 }
292 
293 template<typename T>
295  const T &defValue) const
296 {
297  std::map<AlpgenHeader::Parameter, double>::const_iterator pos =
298  header.params.find(index);
299  if (pos == header.params.end())
300  return defValue;
301  else
302  return T(pos->second);
303 }
304 
305 unsigned int AlpgenSource::processID() const
306 {
307  // return 661; // The original, old thing.
308  // digits #XYZ: X = ihrd, Y = ihvy, Z = njets
309  return header.ihrd * 100 +
310  getParameter<unsigned int>(AlpgenHeader::ihvy, 0) * 10 +
311  getParameter<unsigned int>(AlpgenHeader::njets, 0);
312 }
313 
315 {
316  char buffer[512];
317  double dummy;
318  int nPart;
319  double sWgtRes;
320  double sQ;
321 
322  inputFile_->getline(buffer, sizeof buffer);
323  if (!inputFile_->good())
324  return false;
325 
326  std::istringstream ls(buffer);
327 
328  // Event number and process don't matter (or do they?)
329  ls >> dummy >> dummy;
330 
331  // Number of particles in the record, sample's average weight and Q scale
332  ls >> nPart >> sWgtRes >> sQ;
333 
334  if (ls.bad() || nPart < 1 || nPart > 1000)
335  return false;
336 
337  // Make room for the particles listed in the Alpgen file
338  hepeup.resize(nPart);
339 
340  // Scales, weights and process ID.
341  hepeup.SCALUP = sQ;
342  hepeup.XWGTUP = sWgtRes;
343  hepeup.IDPRUP = processID();
344 
345  // Incoming lines
346  for(int i = 0; i != 2; ++i) {
347  inputFile_->getline(buffer, sizeof buffer);
348  std::istringstream ls(buffer);
349  int flavour; ls >> flavour;
350  int colour1; ls >> colour1;
351  int colour2; ls >> colour2;
352  double zMomentum; ls >> zMomentum;
353 
354  if (inputFile_->bad())
355  return false;
356 
357  // Setting the HEPEUP of the incoming lines.
358  hepeup.IDUP[i] = flavour;
359  hepeup.ISTUP[i] = -1;
360  hepeup.MOTHUP[i].first = 0;
361  hepeup.MOTHUP[i].second = 0;
362  hepeup.PUP[i][0] = 0.;
363  hepeup.PUP[i][1] = 0.;
364  hepeup.PUP[i][2] = zMomentum;
365  hepeup.PUP[i][3] = std::fabs(zMomentum);
366  hepeup.PUP[i][4] = 0.;
367  if (colour1) colour1 += 500;
368  if (colour2) colour2 += 500;
369  hepeup.ICOLUP[i].first = colour1;
370  hepeup.ICOLUP[i].second = colour2;
371  }
372 
373  // Outgoing lines
374  for(int i = 2; i != nPart; ++i) {
375  inputFile_->getline(buffer, sizeof buffer);
376  std::istringstream ls(buffer);
377  int flavour; ls >> flavour;
378  int colour1; ls >> colour1;
379  int colour2; ls >> colour2;
380  double px, py, pz, mass;
381  ls >> px >> py >> pz >> mass;
382  double energy = std::sqrt(px*px + py*py + pz*pz + mass*mass);
383 
384  if (inputFile_->bad())
385  return false;
386 
387  // Setting the HEPEUP of the outgoing lines.
388  hepeup.IDUP[i] = flavour;
389  hepeup.ISTUP[i] = 1;
390  hepeup.MOTHUP[i].first = 1;
391  hepeup.MOTHUP[i].second = 2;
392  hepeup.PUP[i][0] = px;
393  hepeup.PUP[i][1] = py;
394  hepeup.PUP[i][2] = pz;
395  hepeup.PUP[i][3] = energy;
396  hepeup.PUP[i][4] = mass;
397  if (colour1) colour1 += 500;
398  if (colour2) colour2 += 500;
399  hepeup.ICOLUP[i].first = colour1;
400  hepeup.ICOLUP[i].second = colour2;
401  }
402 
403  return true;
404 }
405 
407 {
408  // The LHE Event Record
409  hepeup_.reset(new lhef::HEPEUP);
410 
411  lhef::HEPEUP& hepeup = *hepeup_;
412 
413  // Read the .unw file until it is over.
414  for(;;) {
415  if (!readAlpgenEvent(hepeup)) {
416  if (inputFile_->eof())
417  return false;
418 
419  throw cms::Exception("Generator|AlpgenInterface")
420  << "AlpgenSource is not able read event no. "
421  << nEvent_ << std::endl;
422  }
423 
424  nEvent_++;
425  if (skipEvents_ > 0)
426  skipEvents_--;
427  else
428  break;
429  }
430  return true;
431 }
432 
434 {
435 
436  // Here are the Alpgen routines for filling up the rest
437  // of the LHE Event Record. The .unw file has the information
438  // in a compressed way, e.g. it doesn't list the W boson -
439  // one has to reconstruct it from the e nu pair.
440  lhef::HEPEUP& hepeup = *hepeup_;
441 
442  switch(header.ihrd) {
443  case 1:
444  case 2:
445  case 3:
446  case 4:
447  case 10:
448  case 14:
449  case 15:
450  alpgen::fixEventWZ(hepeup);
451  break;
452  case 5:
454  break;
455  case 6:
456  alpgen::fixEventTTbar(hepeup);
457  break;
458  case 8:
460  break;
461  case 13:
464  break;
465  case 7:
466  case 9:
467  case 11:
468  case 12:
469  case 16:
470  // No fixes needed.
471  break;
472 
473  default:
474  throw cms::Exception("Generator|AlpgenInterface")
475  << "Unrecognized IHRD process code" << std::endl;
476  }
477 
478  // Create the LHEEventProduct and put it into the Event.
479  std::auto_ptr<LHEEventProduct> lheEvent(new LHEEventProduct(hepeup));
480  event.put(lheEvent);
481 
482  hepeup_.reset();
483 }
484 
int i
Definition: DBlmapReader.cc:9
virtual ~AlpgenSource()
Destructor.
void resize(int nrup)
Definition: LesHouches.h:52
std::map< Parameter, double > params
Definition: AlpgenHeader.h:67
AlpgenSource(const edm::ParameterSet &params, const edm::InputSourceDescription &desc)
Constructor.
Definition: AlpgenSource.cc:94
bool parse(const std::vector< std::string >::const_iterator &begin, const std::vector< std::string >::const_iterator &end)
Definition: AlpgenHeader.cc:67
std::pair< double, double > EBMUP
Definition: LesHouches.h:78
AlpgenHeader header
Alpgen _unw.par file as an AlpgenHeader.
Definition: AlpgenSource.cc:78
def ls
Definition: eostools.py:346
std::pair< int, int > IDBMUP
Definition: LesHouches.h:73
RunNumber_t run() const
std::unique_ptr< std::ifstream > inputFile_
Pointer to the input file.
Definition: AlpgenSource.cc:67
static std::string parameterName(Parameter index)
const_iterator end() const
bool writeAlpgenParFile
Definition: AlpgenSource.cc:90
bool writeExtraHeader
Definition: AlpgenSource.cc:91
virtual bool setRunAndEventInfo(edm::EventID &, edm::TimeValue_t &, edm::EventAuxiliary::ExperimentType &) override
unsigned long nEvent_
Number of events.
Definition: AlpgenSource.cc:73
unsigned long skipEvents_
Number of events to skip.
Definition: AlpgenSource.cc:70
std::vector< std::string > const & fileNames() const
Definition: FromFiles.h:22
std::pair< int, int > PDFGUP
Definition: LesHouches.h:84
bool writeAlpgenWgtFile
configuration flags
Definition: AlpgenSource.cc:89
std::unique_ptr< lhef::HEPEUP > hepeup_
Definition: AlpgenSource.cc:80
void fixEventWZ(lhef::HEPEUP &hepeup)
Fixes Event Record for ihrd = 1,2,3,4,10,14,15.
std::vector< std::pair< int, int > > MOTHUP
Definition: LesHouches.h:236
LHERunInfoProduct::Header lheAlpgenUnwParHeader
Alpgen _unw.par file as a LHE header.
Definition: AlpgenSource.cc:76
void fixEventMultiBoson(lhef::HEPEUP &hepeup)
Fixes Event Record for ihrd = 5.
#define DEFINE_FWK_INPUT_SOURCE(type)
void fixEventHiggsTTbar(lhef::HEPEUP &hepeup)
Fixes Event Record for ihrd = 8.
void resize(int nup)
Definition: LesHouches.h:161
std::string extraHeaderFileName_
Name of the extra header file.
Definition: AlpgenSource.cc:83
T sqrt(T t)
Definition: SSEVec.h:48
std::vector< FiveVector > PUP
Definition: LesHouches.h:248
void addLine(const std::string &line)
std::vector< int > ISTUP
Definition: LesHouches.h:230
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
unsigned long long TimeValue_t
Definition: Timestamp.h:28
virtual void produce(edm::Event &event) override
std::vector< int > IDUP
Definition: LesHouches.h:225
std::vector< double > XERRUP
Definition: LesHouches.h:114
std::string fileName_
Name of the input file.
Definition: AlpgenSource.cc:64
std::vector< double > XMAXUP
Definition: LesHouches.h:119
TString nPart(Int_t part, TString string, TString delimit=";", Bool_t removerest=true)
std::string slhaMassLine(int pdgId, AlpgenHeader::Masses mass, const std::string &comment) const
std::string extraHeaderName_
Name given to the extra header.
Definition: AlpgenSource.cc:86
void fixEventSingleTop(lhef::HEPEUP &hepeup, double mb, int itopprc)
Fixes Event Record for ihrd = 13.
virtual void beginRun(edm::Run &run) override
std::pair< int, int > PDFSUP
Definition: LesHouches.h:90
unsigned int ihrd
Definition: AlpgenHeader.h:68
EventNumber_t event() const
T getParameter(AlpgenHeader::Parameter index) const
Function to get parameter by name from AlpgenHeader.
double xsecErr
Definition: AlpgenHeader.h:70
void put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Run.h:107
void fixEventTTbar(lhef::HEPEUP &hepeup)
Fixes Event Record for ihrd = 6.
volatile std::atomic< bool > shutdown_flag false
std::vector< double > XSECUP
Definition: LesHouches.h:108
long double T
int flavour(const Candidate &part)
Definition: pdgIdUtils.h:31
double XWGTUP
Definition: LesHouches.h:196
double masses[MASS_MAX]
Definition: AlpgenHeader.h:73
bool readAlpgenEvent(lhef::HEPEUP &hepeup)
Read an event and put it into the HEPEUP.
std::vector< std::pair< int, int > > ICOLUP
Definition: LesHouches.h:242
const_iterator begin() const
Definition: Run.h:41
std::vector< int > LPRUP
Definition: LesHouches.h:124
double SCALUP
Definition: LesHouches.h:210
unsigned int processID() const
#define comment(par)
Definition: vmac.h:161