diff --git a/MatrixElement/FxFx/FxFxEventHandler.cc b/MatrixElement/FxFx/FxFxEventHandler.cc --- a/MatrixElement/FxFx/FxFxEventHandler.cc +++ b/MatrixElement/FxFx/FxFxEventHandler.cc @@ -1,652 +1,673 @@ // -*- C++ -*- // // Based on: // FxFxEventHandler.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2019 Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the FxFxEventHandler class. // #include "FxFxEventHandler.h" #include "FxFxReader.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/RefVector.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Repository/Repository.h" #include "ThePEG/Handlers/LuminosityFunction.h" #include "ThePEG/Handlers/XComb.h" #include "ThePEG/Handlers/CascadeHandler.h" #include "ThePEG/Cuts/Cuts.h" #include "ThePEG/PDF/PartonExtractor.h" #include "ThePEG/Utilities/LoopGuard.h" #include "ThePEG/EventRecord/Event.h" #include "ThePEG/EventRecord/Collision.h" #include "ThePEG/EventRecord/Step.h" #include "ThePEG/EventRecord/SubProcess.h" #include "ThePEG/Utilities/EnumIO.h" #include "ThePEG/Utilities/Maths.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Utilities/Throw.h" #include "ThePEG/Utilities/Debug.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace ThePEG; FxFxEventHandler::~FxFxEventHandler() {} IBPtr FxFxEventHandler::clone() const { return new_ptr(*this); } IBPtr FxFxEventHandler::fullclone() const { return new_ptr(*this); } void FxFxEventHandler::doinit() { EventHandler::doinit(); for ( int i = 0, N = readers().size(); i < N; ++i ) { readers()[i]->init(); //FxFxReader & reader = *readers()[i]; //reader.initialize(*this); //weightnames = reader.optWeightsNamesFunc(); } /* XSecStat* initxsecs = new XSecStat[weightnames.size()]; for(int ww = 0; ww < weightnames.size(); ww++){ initxsecs[ww].reset(); optstats.insert(std::make_pair(weightnames[ww], initxsecs[ww])); opthistStats.insert(std::make_pair(weightnames[ww], initxsecs[ww])); CrossSection initxs = 0.*picobarn; optxs.insert(std::make_pair(weightnames[ww], initxs)); }*/ ntries = 0; } void FxFxEventHandler::initialize() { if ( lumiFnPtr() ) Repository::clog() << "The LuminosityFunction '" << lumiFnPtr()->name() << "' assigned to the FxFxEventHandler '" << name() << "' will not be active in this run. Instead the incoming " << "particles will be determined by the used FxFxReader objects.\n" << Exception::warning; if ( readers().empty() ) throw FxFxInitError() << "No readers were defined for the FxFxEventHandler '" << name() << "'" << Exception::warning; // Go through all the readers and collect information about cross // sections and processes. typedef map ProcessMap; ProcessMap processes; PDPair incoming; Energy MaxEA = ZERO; Energy MaxEB = ZERO; for ( int i = 0, N = readers().size(); i < N; ++i ) { FxFxReader & reader = *readers()[i]; reader.initialize(*this); weightnames = reader.optWeightsNamesFunc(); // Check that the incoming particles are consistent between the // readers. if ( !incoming.first ) { incoming.first = getParticleData(reader.heprup.IDBMUP.first); if ( !incoming.first ) Throw() << "Unknown beam PID " << reader.heprup.IDBMUP.first << ". Have you created a matching BeamParticle object?" << Exception::runerror; } if ( !incoming.second ) { incoming.second = getParticleData(reader.heprup.IDBMUP.second); if ( !incoming.second ) Throw() << "Unknown beam PID " << reader.heprup.IDBMUP.first << ". Have you created a matching BeamParticle object?" << Exception::runerror; } if ( incoming.first->id() != reader.heprup.IDBMUP.first || incoming.second->id() != reader.heprup.IDBMUP.second ) Repository::clog() << "The different FxFxReader objects in the " << "FxFxEventHandler '" << name() << "' have different " << "types of colliding particles." << Exception::warning; MaxEA = max(MaxEA, reader.heprup.EBMUP.first*GeV); MaxEB = max(MaxEB, reader.heprup.EBMUP.second*GeV); // Check that the weighting of the events in the different readers // is consistent with the ones requested for this event // handler. Also collect the sum of the maximum weights. if ( reader.negativeWeights() && weightOption() > 0 ) throw FxFxInitError() << "The reader '" << reader.name() << "' contains negatively weighted events, " << "which is not allowed for the FxFxEventHandler '" << name() << "'." << Exception::warning; // Check that we do not have the same process numbers in different // readers. for ( int ip = 0; ip < reader.heprup.NPRUP; ++ip ) { if ( reader.heprup.LPRUP[ip] ) { ProcessMap::iterator pit = processes.find(reader.heprup.LPRUP[ip]); if ( pit == processes.end() ) processes[reader.heprup.LPRUP[ip]] = readers()[i]; else if ( warnPNum ) { Throw() << "In the FxFxEventHandler '" << name() << "', both the '" << pit->second->name() << "' and '" << reader.name() << "' contains sub-process number " << pit->first << ". This process may be double-counted in this run." << Exception::warning; } } } selector().insert(reader.stats.maxXSec(), i); } stats.maxXSec(selector().sum()); histStats.maxXSec(selector().sum()); for (map::iterator it= optstats.begin(); it!=optstats.end(); ++it){ (it->second).maxXSec(selector().sum()); } // Check that we have any cross section at all. if ( stats.maxXSec() <= ZERO ) throw FxFxInitError() << "The sum of the cross sections of the readers in the " << "FxFxEventHandler '" << name() << "' was zero." << Exception::warning; // We now create a LuminosityFunction object to inform others about // the energy of the beam. theIncoming = incoming; lumiFn(new_ptr(LuminosityFunction(MaxEA, MaxEB))); } void FxFxEventHandler::doinitrun() { EventHandler::doinitrun(); stats.reset(); histStats.reset(); for ( int i = 0, N = readers().size(); i < N; ++i ) { readers()[i]->init(); FxFxReader & reader = *readers()[i]; reader.initialize(*this); weightnames = reader.optWeightsNamesFunc(); } // XSecStat initxsecs; XSecStat* initxsecs = new XSecStat[weightnames.size()]; for(size_t ww = 0; ww < weightnames.size(); ww++){ initxsecs[ww].reset(); // optstats.insert(std::make_pair(weightnames[ww], initxsecs[ww])); // opthistStats.insert(std::make_pair(weightnames[ww], initxsecs[ww])); CrossSection initxs = 0.*picobarn; //optxs.insert(std::make_pair(weightnames[ww], initxs)); optstats.insert(std::make_pair(weightnames[ww], initxsecs[ww])); opthistStats.insert(std::make_pair(weightnames[ww], initxsecs[ww])); optxs.insert(std::make_pair(weightnames[ww], initxs)); } ntries = 0; } EventPtr FxFxEventHandler::generateEvent() { LoopGuard loopGuard(*this, maxLoop()); while ( true ) { loopGuard(); currentReader(readers()[selector().select(UseRandom::current())]); skipEvents(); currentReader()->reset(); double weight = currentReader()->getEvent(); if ( weightOption() == unitweight && weight < 0.0 ) weight = 0.0; if ( weightOption() == unitweight || weightOption() == unitnegweight ) { CrossSection newmax = selector().reweight(weight); if ( newmax > CrossSection() ) increaseMaxXSec(newmax); } select(weight/currentReader()->preweight); histStats.select(weight); if ( !weighted() ) { if ( weightOption() == unitweight || weightOption() == unitnegweight ) { if ( !rndbool(abs(weight)) ) continue; weight = Math::sign(1.0, weight); } else if ( weight == 0.0 ) continue; } else if ( weight == 0.0 ) continue; accept(); // Divide by the bias introduced by the preweights in the reader. weight /= currentReader()->preweight; // fact for weight normalization double fact = theNormWeight ? double(selector().sum()/picobarn) : 1.; try { theLastXComb = currentReader()->getXComb(); - currentEvent(new_ptr(Event(lastParticles(), this, generator()->runName(), - generator()->currentEventNumber(), weight*fact))); + //whether to use the LHE event number or not for the event identification + if(UseLHEEvent==0 || currentReader()->LHEEventNum() == -1) { + currentEvent(new_ptr(Event(lastParticles(), this, generator()->runName(), + generator()->currentEventNumber(), weight*fact ))); + } + else if(UseLHEEvent==1 && currentReader()->LHEEventNum() != -1) { + currentEvent(new_ptr(Event(lastParticles(), this, generator()->runName(), + currentReader()->LHEEventNum(), weight*fact ))); + } currentEvent()->optionalWeights() = currentReader()->optionalEventWeights(); // normalize the optional weights for(map::iterator it = currentEvent()->optionalWeights().begin(); it!=currentEvent()->optionalWeights().end();++it) { if(it->first!="ecom"&& it->second!=-999 && it->second!=-111 && it->second!=-222 && it->second!=-333) { it->second *= fact; } } //print optional weights here // cout << "event weight = " << weight << " fact = " << fact << endl; /*for (map::const_iterator it= currentReader()->optionalEventWeights().begin(); it!=currentReader()->optionalEventWeights().end(); ++it){ std::cout << it->first << " => " << it->second << '\n'; } cout << endl;*/ //print npLO and npNLO // cout << currentReader()->optionalEventnpLO() << "\t" << currentReader()->optionalEventnpNLO() << endl; performCollision(); if ( !currentCollision() ) throw Veto(); return currentEvent(); } catch (Veto) { reject(weight); } catch (Stop) { } catch (Exception &) { reject(weight); throw; } } } void FxFxEventHandler::skipEvents() { if ( weightOption() == 2 || weightOption() == -2 ) return; //does it make sense to skip events if we are using varying weights? // Don't do this for readers which seem to generate events on the fly. if ( currentReader()->active() || currentReader()->NEvents() <= 0 ) return; // Estimate the fration of the total events available from // currentReader() which will be requested. double frac = currentReader()->stats.maxXSec()/stats.maxXSec(); if ( stats.accepted() > 0 ) frac *= double(stats.attempts())/double(stats.accepted()); else frac *= double(stats.attempts() + 1); double xscan = generator()->N()*frac/currentReader()->NEvents(); // Estimate the number of times we need to go through the events for // the currentReader(), and how many events on average we need to // skip for each attempted event to go through the file an integer // number of times. double nscan = ceil(xscan); double meanskip = nscan/xscan - 1.0; // Skip an average numer of steps with a Poissonian distribution. currentReader()-> skip(UseRandom::rndPoisson(meanskip)%currentReader()->NEvents()); } void FxFxEventHandler::select(double weight) { stats.select(weight); currentReader()->select(weight); vector w; for (map::const_iterator it= currentReader()->optionalEventWeights().begin(); it!=currentReader()->optionalEventWeights().end(); ++it){ w.push_back(it->second); } int ii = 0; for (map::iterator it= opthistStats.begin(); it!=opthistStats.end(); ++it){ (it->second).select(w[ii]); ii++; } ii = 0; for (map::iterator it= optstats.begin(); it!=optstats.end(); ++it){ (it->second).select(w[ii]); ii++; } } tCollPtr FxFxEventHandler::performCollision() { lastExtractor()->select(lastXCombPtr()); if ( CKKWHandler() ) CKKWHandler()->setXComb(lastXCombPtr()); currentCollision(new_ptr(Collision(lastParticles(), currentEvent(), this))); if ( currentEvent() ) currentEvent()->addCollision(currentCollision()); currentStep(new_ptr(Step(currentCollision(), this))); currentCollision()->addStep(currentStep()); currentStep()->addSubProcess(currentReader()->getSubProcess()); lastExtractor()->constructRemnants(lastXCombPtr()->partonBinInstances(), subProcess(), currentStep()); if ( !currentReader()->cuts().passCuts(*currentCollision()) ) throw Veto(); initGroups(); if ( ThePEG_DEBUG_ITEM(1) ) { if ( currentEvent() ) generator()->logfile() << *currentEvent(); else generator()->logfile() << *currentCollision(); } return continueCollision(); } EventPtr FxFxEventHandler::continueEvent() { try { continueCollision(); } catch (Veto) { const double fact = theNormWeight ? double(selector().sum()/picobarn) : 1.; reject(currentEvent()->weight()/fact); } catch (Stop) { } catch (Exception &) { const double fact = theNormWeight ? double(selector().sum()/picobarn) : 1.; reject(currentEvent()->weight()/fact); throw; } return currentEvent(); } void FxFxEventHandler::dofinish() { EventHandler::dofinish(); if ( selector().compensating() ) generator()->log() << "Warning: The run was ended while the FxFxEventHandler '" << name() << "' was still trying to compensate for weights larger than 1. " << "The cross section estimates may therefore be statistically " << "inaccurate." << endl; } void FxFxEventHandler::statistics(ostream & os) const { if ( statLevel() == 0 ) return; string line = "=======================================" "=======================================\n"; if ( stats.accepted() <= 0 ) { os << line << "No events generated by event handler '" << name() << "'." << endl; return; } os << line << "Statistics for Les Houches event handler \'" << name() << "\':\n" << " " << "generated number of Cross-section\n" << " " << " events attempts (nb)\n"; os << line << "Total:" << setw(42) << stats.accepted() << setw(13) << stats.attempts() << setw(17) << ouniterr(stats.xSec(), stats.xSecErr(), nanobarn) << endl << line; if ( statLevel() == 1 ) return; if ( statLevel() == 2 ) { os << "Per Les Houches Reader breakdown:\n"; for ( int i = 0, N = readers().size(); i < N; ++i ) { FxFxReader & reader = *readers()[i]; string n = reader.name(); n.resize(37, ' '); os << n << setw(11) << reader.stats.accepted() << setw(13) << reader.stats.attempts() << setw(17) << ouniterr(reader.stats.xSec(), reader.stats.xSecErr(), nanobarn) << endl; } os << line; } else { os << "Per Les Houches Reader (and process #) breakdown:\n"; for ( int i = 0, N = readers().size(); i < N; ++i ) { FxFxReader & reader = *readers()[i]; string n = reader.name() + " (all)"; n.resize(37, ' '); os << n << setw(11) << reader.stats.accepted() << setw(13) << reader.stats.attempts() << setw(17) << ouniterr(reader.stats.xSec(), reader.stats.xSecErr(), nanobarn) << endl; CrossSection xsectot = reader.stats.xSec(); if ( xsectot != ZERO ) xsectot /= reader.stats.sumWeights(); typedef FxFxReader::StatMap::const_iterator const_iterator; for ( const_iterator i = reader.statmap.begin(); i != reader.statmap.end(); ++i ) { ostringstream ss; ss << reader.name() << " (" << i->first << ")"; string n = ss.str(); n.resize(37, ' '); os << n << setw(11) << i->second.accepted() << setw(13) << i->second.attempts() << setw(17) << ouniterr(i->second.sumWeights()*xsectot, sqrt(i->second.sumWeights2())*xsectot, nanobarn) << endl; } os << line; } } string warn = "Warning: Result may be statistically incorrect since\n" " the following FxFxReaders were oversampled:\n"; for ( int i = 0, N = readers().size(); i < N; ++i ) { FxFxReader & reader = *readers()[i]; if ( reader.NEvents() > 0 && reader.stats.attempts() > reader.NEvents() ) { os << warn; warn = ""; os << "'" << reader.name() << "' (by a factor " << double(reader.stats.attempts())/double(reader.NEvents()) << ")" << endl; } } } void FxFxEventHandler::increaseMaxXSec(CrossSection maxxsec) { stats.maxXSec(selector().sum()); histStats.maxXSec(selector().sum()); currentReader()->increaseMaxXSec(maxxsec); } void FxFxEventHandler::accept() { ntries++; stats.accept(); histStats.accept(); currentReader()->accept(); for (map::iterator it= opthistStats.begin(); it!=opthistStats.end(); ++it){ (it->second).accept(); } for (map::iterator it= optstats.begin(); it!=optstats.end(); ++it){ (it->second).accept(); } } void FxFxEventHandler::reject(double w) { ntries++; stats.reject(w); histStats.reject(w); currentReader()->reject(w); vector wv; for (map::const_iterator it= currentReader()->optionalEventWeights().begin(); it!=currentReader()->optionalEventWeights().end(); ++it){ wv.push_back(it->second); } int ii = 0; for (map::iterator it= opthistStats.begin(); it!=opthistStats.end(); ++it){ (it->second).reject(wv[ii]); ii++; } ii = 0; for (map::iterator it= optstats.begin(); it!=optstats.end(); ++it){ (it->second).reject(wv[ii]); ii++; } } map FxFxEventHandler::optintegratedXSecMap() const { map result; for (map::const_iterator it= optstats.begin(); it!=optstats.end(); ++it){ result[it->first] = (it->second.sumWeights()/it->second.attempts()) * picobarn; } return result; } CrossSection FxFxEventHandler::histogramScale() const { return histStats.xSec()/histStats.sumWeights(); } CrossSection FxFxEventHandler::integratedXSec() const { return histStats.xSec(); } CrossSection FxFxEventHandler::integratedXSecErr() const { return histStats.xSecErr(); } int FxFxEventHandler::ntriesinternal() const { return stats.attempts(); } void FxFxEventHandler::persistentOutput(PersistentOStream & os) const { os << stats << histStats << theReaders << theSelector << oenum(theWeightOption) << theUnitTolerance << theCurrentReader << warnPNum - << theNormWeight; + << theNormWeight << UseLHEEvent; } void FxFxEventHandler::persistentInput(PersistentIStream & is, int) { is >> stats >> histStats >> theReaders >> theSelector >> ienum(theWeightOption) >> theUnitTolerance >> theCurrentReader >> warnPNum - >> theNormWeight; + >> theNormWeight >> UseLHEEvent;; } ClassDescription FxFxEventHandler::initFxFxEventHandler; // Definition of the static class description member. void FxFxEventHandler::setUnitTolerance(double x) { theUnitTolerance = x; selector().tolerance(unitTolerance()); } void FxFxEventHandler::Init() { static ClassDocumentation documentation ("This is the main class administrating the selection of hard " "subprocesses from a set of ThePEG::FxFxReader objects."); static RefVector interfaceFxFxReaders ("FxFxReaders", "Objects capable of reading events from an event file or an " "external matrix element generator.", &FxFxEventHandler::theReaders, -1, false, false, true, false, false); static Switch interfaceWeightOption ("WeightOption", "The different ways to weight events in the Les Houches event handler. " "Whether weighted or not and whether or not negative weights are allowed.", &FxFxEventHandler::theWeightOption, unitweight, true, false); static SwitchOption interfaceWeightOptionUnitWeight (interfaceWeightOption, "UnitWeight", "All events have unit weight.", unitweight); static SwitchOption interfaceWeightOptionNegUnitWeight (interfaceWeightOption, "NegUnitWeight", "All events have weight +1 or maybe -1.", unitnegweight); static SwitchOption interfaceWeightOptionVarWeight (interfaceWeightOption, "VarWeight", "Events may have varying but positive weights.", varweight); static SwitchOption interfaceWeightOptionVarNegWeight (interfaceWeightOption, "VarNegWeight", "Events may have varying weights, both positive and negative.", varnegweight); static Switch interfaceWarnPNum ("WarnPNum", "Warn if the same process number is used in more than one " "FxFxReader.", &FxFxEventHandler::warnPNum, true, true, false); static SwitchOption interfaceWarnPNumWarning (interfaceWarnPNum, "Warning", "Give a warning message.", true); static SwitchOption interfaceWarnPNumNoWarning (interfaceWarnPNum, "NoWarning", "Don't give a warning message.", false); static Parameter interfaceUnitTolerance ("UnitTolerance", "If the WeightOption is set to unit weight, do not start compensating unless the a weight is found to be this much larger than unity.", &FxFxEventHandler::theUnitTolerance, 1.0e-6, 0.0, 0, true, false, Interface::lowerlim, &FxFxEventHandler::setUnitTolerance, (double(FxFxEventHandler::*)()const)(0), (double(FxFxEventHandler::*)()const)(0), (double(FxFxEventHandler::*)()const)(0), (double(FxFxEventHandler::*)()const)(0)); static Switch interfaceWeightNormalization ("WeightNormalization", "How to normalize the output weights", &FxFxEventHandler::theNormWeight, 0, false, false); static SwitchOption interfaceWeightNormalizationUnit (interfaceWeightNormalization, "Normalized", "Standard normalization, i.e. +/- for unweighted events", 0); static SwitchOption interfaceWeightNormalizationCrossSection (interfaceWeightNormalization, "CrossSection", "Normalize the weights to the max cross section in pb", 1); + static Switch interfaceEventNumbering + ("EventNumbering", + "How to number the events", + &FxFxEventHandler::UseLHEEvent, 0, false, false); + static SwitchOption interfaceEventNumberingIncremental + (interfaceEventNumbering, + "Incremental", + "Standard incremental numbering (i.e. as they are generated)", + 0); + static SwitchOption interfaceEventNumberingLHE + (interfaceEventNumbering, + "LHE", + "Corresponding to the LHE event number", + 1); interfaceFxFxReaders.rank(10); interfaceWeightOption.rank(9); } diff --git a/MatrixElement/FxFx/FxFxEventHandler.h b/MatrixElement/FxFx/FxFxEventHandler.h --- a/MatrixElement/FxFx/FxFxEventHandler.h +++ b/MatrixElement/FxFx/FxFxEventHandler.h @@ -1,447 +1,452 @@ // -*- C++ -*- // // FxFxEventHandler.h is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2019 Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef THEPEG_FxFxEventHandler_H #define THEPEG_FxFxEventHandler_H // // This is the declaration of the FxFxEventHandler class. // #include "ThePEG/Handlers/EventHandler.h" #include "FxFxEventHandler.fh" #include "FxFxReader.fh" #include "ThePEG/Utilities/CompSelector.h" #include "ThePEG/Utilities/XSecStat.h" namespace ThePEG { /** * The FxFxEventHandler inherits from the general EventHandler * class and administers the reading of events generated by external * matrix element generator programs according to the Les Houches * accord. * * The class has a list of FxFxReaders which * typically are connected to files with event data produced by * external matrix element generator programs. When an event is * requested by FxFxEventHandler, one of the readers are chosen, * an event is read in and then passed to the different * StepHandler defined in the underlying * EventHandler class. * * @see \ref FxFxEventHandlerInterfaces "The interfaces" * defined for FxFxEventHandler. */ class FxFxEventHandler: public EventHandler { public: /** * A vector of FxFxReader objects. */ typedef vector ReaderVector; /** * A selector of readers. */ typedef CompSelector ReaderSelector; /** * Enumerate the weighting options. */ enum WeightOpt { unitweight = 1, /**< All events have unit weight. */ unitnegweight = -1, /**< All events have wight +/- 1. */ varweight = 2, /**< Varying positive weights. */ varnegweight = -2 /**< Varying positive or negative weights. */ }; friend class FxFxHandler; public: /** @name Standard constructors and destructors. */ //@{ /** * The default constructor. */ FxFxEventHandler() - : theWeightOption(unitweight), theUnitTolerance(1.0e-6), warnPNum(true), theNormWeight(0) + : theWeightOption(unitweight), theUnitTolerance(1.0e-6), warnPNum(true), theNormWeight(0), UseLHEEvent(0) { selector().tolerance(unitTolerance()); } /** * The destructor. */ virtual ~FxFxEventHandler(); //@} public: /** @name Initialization and finalization functions. */ //@{ /** * Initialize this event handler and all related objects needed to * generate events. */ virtual void initialize(); /** * Write out accumulated statistics about intergrated cross sections * and stuff. */ virtual void statistics(ostream &) const; /** * Histogram scale. A histogram bin which has been filled with the * weights associated with the Event objects should be scaled by * this factor to give the correct cross section. */ virtual CrossSection histogramScale() const; /** * The estimated total integrated cross section of the processes * generated in this run. * @return 0 if no integrated cross section could be estimated. */ virtual CrossSection integratedXSec() const; virtual int ntriesinternal() const; /** * The estimated error in the total integrated cross section of the * processes generated in this run. * @return 0 if no integrated cross section error could be estimated. */ virtual CrossSection integratedXSecErr() const; virtual map optintegratedXSecMap() const; //@} /** @name Functions used for the actual generation */ //@{ /** * Generate an event. */ virtual EventPtr generateEvent(); /** * Create the Event and Collision objects. Used by the * generateEvent() function. */ virtual tCollPtr performCollision(); /** * Continue generating an event if the generation has been stopped * before finishing. */ virtual EventPtr continueEvent(); //@} /** @name Functions to manipulate statistics. */ //@{ /** * An event has been selected. Signal that an event has been * selected with the given \a weight. If unit weights are requested, * the event will be accepted with that weight. This also takes care * of the statistics collection of the selected reader object. */ void select(double weight); /** * Accept the current event, taking care of the statistics * collection of the corresponding reader objects. */ void accept(); /** * Reject the current event, taking care of the statistics * collection of the corresponding reader objects. */ void reject(double weight); /** * Increase the overestimated cross section for the selected reader. */ void increaseMaxXSec(CrossSection maxxsec); /** * Skip some events. To ensure a reader file is scanned an even * number of times, skip a number of events for the selected reader. */ void skipEvents(); //@} /** @name Simple access functions. */ //@{ /** * The way weights are to be treated. */ WeightOpt weightOption() const { return theWeightOption; } /** * If the weight option is set to unit weight, do not start * compensating unless the weight is this much larger than unity. */ double unitTolerance() const { return theUnitTolerance; } /** * Access the list of readers. */ const ReaderVector & readers() const { return theReaders; } /** * The selector to choose readers according to their overestimated * cross section. */ const ReaderSelector & selector() const { return theSelector; } /** * The currently selected reader object. */ tFxFxReaderPtr currentReader() const { return theCurrentReader; } /** * Set the currently selected reader object. */ void currentReader(tFxFxReaderPtr x) { theCurrentReader = x; } //@} public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); /** * The currently selected reader object. */ tFxFxReaderPtr theCurrentReader; protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); /** * Initialize this object. Called in the run phase just before * a run begins. */ virtual void doinitrun(); /** * Finalize this object. Called in the run phase just after a * run has ended. Used eg. to write out statistics. */ virtual void dofinish(); //@} protected: /** * Access the list of readers. */ ReaderVector & readers() { return theReaders; } /** * The selector to choose readers according to their overestimated * cross section. */ ReaderSelector & selector() { return theSelector; } /** * Helper function for the interface; */ void setUnitTolerance(double); /** * Collect statistics for this event handler. */ XSecStat stats; map optstats; map optxs; int ntries; map OptStatsFunc() { return optstats; } map OptXsFunc() { return optxs; } /** * Collect statistics for this event handler. To be used for * histogram scaling. */ XSecStat histStats; map opthistStats; /* * The weight identifiers for the events */ vector weightnames; private: /** * The list of readers. */ ReaderVector theReaders; /** * The selector to choose readers according to their overestimated * cross section. */ ReaderSelector theSelector; /** * The way weights are to be treated. */ WeightOpt theWeightOption; /** * If the weight option is set to unit weight, do not start * compensating unless the weight is this much larger than unity. */ double theUnitTolerance; /** * Warn if the same process number is used in more than one * FxFxReader. */ bool warnPNum; /** * How to normalize the weights */ unsigned int theNormWeight; + /** + * How to number the events + */ + unsigned int UseLHEEvent; + public: /** @cond EXCEPTIONCLASSES */ /** * Exception class used if no readers were assigned. */ class FxFxInitError: public InitException {}; /** * Exception class used if the same process number is used by more * than ne reader. */ class FxFxPNumException: public InitException {}; /** @endcond */ private: /** * The static object used to initialize the description of this class. * Indicates that this is a concrete class with persistent data. */ static ClassDescription initFxFxEventHandler; /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ FxFxEventHandler & operator=(const FxFxEventHandler &) = delete; }; } // CLASSDOC OFF #include "ThePEG/Utilities/ClassTraits.h" namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ /** This template specialization informs ThePEG about the * base classes of FxFxEventHandler. */ template <> struct BaseClassTrait { /** Typedef of the first base class of FxFxEventHandler. */ typedef EventHandler NthBase; }; /** This template specialization informs ThePEG about the name of * the FxFxEventHandler class and the shared object where it is defined. */ template <> struct ClassTraits : public ClassTraitsBase { /** Return a platform-independent class name */ static string className() { return "Herwig::FxFxEventHandler"; } /** Return the name of the shared library be loaded to get access to * the FxFxEventHandler class and every other class it uses * (except the base class). */ static string library() { return "HwFxFx.so"; } }; /** @endcond */ } #endif /* THEPEG_FxFxEventHandler_H */ diff --git a/MatrixElement/FxFx/FxFxFileReader.cc b/MatrixElement/FxFx/FxFxFileReader.cc --- a/MatrixElement/FxFx/FxFxFileReader.cc +++ b/MatrixElement/FxFx/FxFxFileReader.cc @@ -1,929 +1,946 @@ // -*- C++ -*- // // FxFxFileReader.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2019 Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the FxFxFileReader class. // #include "FxFxFileReader.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Utilities/Throw.h" #include "ThePEG/PDT/DecayMode.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include #include using namespace ThePEG; FxFxFileReader:: FxFxFileReader(const FxFxFileReader & x) : FxFxReader(x), neve(x.neve), ieve(0), LHFVersion(x.LHFVersion), outsideBlock(x.outsideBlock), headerBlock(x.headerBlock), initComments(x.initComments), initAttributes(x.initAttributes), eventComments(x.eventComments), eventAttributes(x.eventAttributes), theFileName(x.theFileName), theQNumbers(x.theQNumbers), theIncludeFxFxTags(x.theIncludeFxFxTags), theIncludeCentral(x.theIncludeCentral), theDecayer(x.theDecayer) {} FxFxFileReader::~FxFxFileReader() {} IBPtr FxFxFileReader::clone() const { return new_ptr(*this); } IBPtr FxFxFileReader::fullclone() const { return new_ptr(*this); } bool FxFxFileReader::preInitialize() const { return true; } void FxFxFileReader::doinit() { FxFxReader::doinit(); // are we using QNUMBERS if(!theQNumbers) return; // parse the header block and create // any new particles needed in QNUMBERS blocks string block = headerBlock; string line = ""; bool readingSLHA = false; int (*pf)(int) = tolower; unsigned int newNumber(0); do { line = StringUtils::car(block,"\r\n"); block = StringUtils::cdr(block,"\r\n"); if(line[0]=='#') continue; // are we reading the SLHA block if(readingSLHA) { // reached the end of slha block ? if(line.find(" split = StringUtils::split(line,"#"); // check for a qnumbers block transform(split[0].begin(), split[0].end(), split[0].begin(), pf); // if not contine if(split[0].find("block qnumbers")==string::npos) continue; // get name from comment string name; if(split.size()>=2) { name = StringUtils::stripws(split[1]); } else { ++newNumber; ostringstream tname; tname << "NP" << newNumber; name = tname.str(); } // extract the PDG code split = StringUtils::split(split[0]," "); istringstream is(split[2]); long PDGCode(0); is >> PDGCode; // get the charge, spin, colour and whether an antiparticle int charge(0),spin(0),colour(0),anti(0); for(unsigned int ix=0;ix<4;++ix) { line = StringUtils::car(block,"\r\n"); block = StringUtils::cdr(block,"\r\n"); int dummy[2]; istringstream is(line); is >> dummy[0] >> dummy[1]; switch (dummy[0]) { case 1: charge = dummy[1]; break; case 2: spin = dummy[1]; break; case 3: colour = dummy[1]; break; case 4: anti = dummy[1]; break; default: assert(false); } } // check if particles already exist PDPair newParticle; newParticle.first = getParticleData(PDGCode); if(newParticle.first) Throw() << "Particle with PDG code " << PDGCode << " whose creation was requested in a QNUMBERS Block" << " already exists. Retaining the original particle" << Exception::warning; if(anti) { newParticle.second = getParticleData(-PDGCode); if(newParticle.second) Throw() << "Anti-particle with PDG code " << -PDGCode << " whose creation was requested in a QNUMBERS Block" << " already exists. Retaining the original particle" << Exception::warning; if(( newParticle.first && !newParticle.second ) || ( newParticle.second && !newParticle.first ) ) Throw() << "Either particle or anti-particle with PDG code " << PDGCode << " whose creation was requested in a QNUMBERS Block" << " already exists, but not both the particle and antiparticle. " << " Something dodgy here stopping" << Exception::runerror; } // already exists continue if(newParticle.first) continue; // create the particles // particle with no anti particle if( anti == 0 ) { // construct the name if(name=="") { ostringstream temp; temp << PDGCode; name = temp.str(); } // create the ParticleData object newParticle.first = ParticleData::Create(PDGCode,name); } // particle anti-particle pair else { // construct the names string nameAnti; if(name=="") { ostringstream temp; temp << PDGCode; name = temp.str(); ostringstream temp2; temp << -PDGCode; nameAnti = temp2.str(); } else { nameAnti=name; for(string::iterator it=nameAnti.begin();it!=nameAnti.end();++it) { if(*it=='+') nameAnti.replace(it,it+1,"-"); else if(*it=='-') nameAnti.replace(it,it+1,"+"); } if(nameAnti==name) nameAnti += "bar"; } // create the ParticleData objects newParticle = ParticleData::Create(PDGCode,name,nameAnti); } // set the particle properties if(colour==1) colour = 0; newParticle.first->iColour(PDT::Colour(colour)); newParticle.first->iSpin (PDT::Spin (spin )); newParticle.first->iCharge(PDT::Charge(charge)); // register it generator()->preinitRegister(newParticle.first, "/Herwig/Particles/"+newParticle.first->PDGName()); // set the antiparticle properties if(newParticle.second) { if(colour==3||colour==6) colour *= -1; charge = -charge; newParticle.second->iColour(PDT::Colour(colour)); newParticle.second->iSpin (PDT::Spin (spin )); newParticle.second->iCharge(PDT::Charge(charge)); // register it generator()->preinitRegister(newParticle.second, "/Herwig/Particles/"+newParticle.second->PDGName()); } } // start of SLHA block ? else if(line.find("> id >> mass; // skip resetting masses on SM particles // as it can cause problems later on in event generation if(abs(id)<=6 || (abs(id)>=11 && abs(id)<=16) || abs(id)==23 || abs(id)==24) { // Throw() << "Standard model mass for PID " // << id // << " will not be changed." // << Exception::warning; block = StringUtils::cdr(block,"\r\n"); line = StringUtils::car(block,"\r\n"); continue; } // magnitude of mass for susy models mass = abs(mass); // set the mass tPDPtr particle = getParticleData(id); if(!particle) throw SetupException() << "FxFxFileReader::doinit() - Particle with PDG code not" << id << " not found." << Exception::runerror; const InterfaceBase * ifb = BaseRepository::FindInterface(particle, "NominalMass"); ostringstream os; os << mass; ifb->exec(*particle, "set", os.str()); // read the next line block = StringUtils::cdr(block,"\r\n"); line = StringUtils::car(block,"\r\n"); }; } // found a decay block else if(line.find("decay") == 0) { // get PGD code and width istringstream iss(line); string dummy; long parent(0); Energy width(ZERO); iss >> dummy >> parent >> iunit(width, GeV); // get the ParticleData object PDPtr inpart = getParticleData(parent); if(!inpart) { throw SetupException() << "FxFxFileReader::doinit() - A ParticleData object with the PDG code " << parent << " does not exist. " << Exception::runerror; return; } if ( abs(inpart->id()) == 6 || abs(inpart->id()) == 15 || abs(inpart->id()) == 23 || abs(inpart->id()) == 24 || abs(inpart->id()) == 25 ) { Throw() << "\n" "************************************************************************\n" "* Your LHE file changes the width of " << inpart->PDGName() << ".\n" "* This can cause serious problems in the event generation!\n" "************************************************************************\n" "\n" << Exception::warning; } else if (inpart->width() > ZERO && width <= ZERO) { Throw() << "\n" "************************************************************************\n" "* Your LHE file zeroes the non-zero width of " << inpart->PDGName() << ".\n" "* If " << inpart->PDGName() << " is a decaying SM particle,\n" "* this can cause serious problems in the event generation!\n" "************************************************************************\n" "\n" << Exception::warning; } // set the width inpart->width(width); if( width > ZERO ) { inpart->cTau(hbarc/width); inpart->widthCut(5.*width); inpart->stable(false); } // construct prefix for DecayModes string prefix(inpart->name() + "->"), tag(prefix),line(""); unsigned int nmode(0); // read any decay modes line = StringUtils::car(block,"\r\n"); while(line[0] != 'D' && line[0] != 'B' && line[0] != 'd' && line[0] != 'b' && line[0] != '<' && line != "") { // skip comments if(line[0] == '#') { block = StringUtils::cdr(block,"\r\n"); line = StringUtils::car(block,"\r\n"); continue; } // read decay mode and construct the tag istringstream is(line); double brat(0.); unsigned int nda(0),npr(0); is >> brat >> nda; while( true ) { long t; is >> t; if( is.fail() ) break; if( t == abs(parent) ) throw SetupException() << "An error occurred while read a decay of the " << inpart->PDGName() << ". One of its products has the same PDG code " << "as the parent particle in FxFxFileReader::doinit()." << " Please check the Les Houches file.\n" << Exception::runerror; tcPDPtr p = getParticleData(t); if( !p ) throw SetupException() << "FxFxFileReader::doinit() -" << " An unknown PDG code has been encounterd " << "while reading a decay mode. ID: " << t << Exception::runerror; ++npr; tag += p->name() + ","; } if( npr != nda ) throw SetupException() << "FxFxFileReader::doinit() - While reading a decay of the " << inpart->PDGName() << " from an SLHA file, an inconsistency " << "between the number of decay products and the value in " << "the 'NDA' column was found. Please check if the spectrum " << "file is correct.\n" << Exception::warning; // create the DecayMode if( npr > 1 ) { if( nmode==0 ) { generator()->preinitInterface(inpart, "VariableRatio" , "set","false"); if(inpart->massGenerator()) { ok = false; Throw() << inpart->PDGName() << " already has a WidthGenerator set" << " this is incompatible with using QNUMBERS " << "Use\n" << "set " << inpart->fullName() << ":Width_generator NULL\n" << "to fix this." << Exception::warning; } unsigned int ntemp=0; for(DecaySet::const_iterator dit = inpart->decayModes().begin(); dit != inpart->decayModes().end(); ++dit ) { if((**dit).on()) ++ntemp; } if(ntemp!=0) { ok = false; Throw() << inpart->PDGName() << " already has DecayModes" << " this is incompatible with using QNUMBERS " << "Use\n" << "do " << inpart->fullName() << ":SelectDecayModes none\n" << " to fix this." << Exception::warning; } } inpart->stable(false); tag.replace(tag.size() - 1, 1, ";"); DMPtr dm = generator()->findDecayMode(tag); if(!theDecayer) Throw() << "FxFxFileReader::doinit() Decayer must be set using the " << "FxFxFileReader:Decayer" << " must be set to allow the creation of new" << " decay modes." << Exception::runerror; if(!dm) { dm = generator()->preinitCreateDecayMode(tag); if(!dm) Throw() << "FxFxFileReader::doinit() - Needed to create " << "new decaymode but one could not be created for the tag " << tag << Exception::warning; } generator()->preinitInterface(dm, "Decayer", "set", theDecayer->fullName()); ostringstream br; br << setprecision(13) << brat; generator()->preinitInterface(dm, "BranchingRatio", "set", br.str()); generator()->preinitInterface(dm, "Active", "set", "Yes"); if(dm->CC()) { generator()->preinitInterface(dm->CC(), "BranchingRatio", "set", br.str()); generator()->preinitInterface(dm->CC(), "Active", "set", "Yes"); } ++nmode; } tag=prefix; // read the next line block = StringUtils::cdr(block,"\r\n"); line = StringUtils::car(block,"\r\n"); }; if(nmode>0) { inpart->update(); if(inpart->CC()) inpart->CC()->update(); } } } // start of SLHA block ? else if(line.find("() << "The file associated with '" << name() << "' does not contain a " << "proper formatted Les Houches event file. The events may not be " << "properly sampled." << Exception::warning; } //vector FxFxFileReader::optWeightNamesFunc() { return optionalWeightsNames; } vector FxFxFileReader::optWeightsNamesFunc() { return optionalWeightsNames; } void FxFxFileReader::open() { if ( filename().empty() ) throw FxFxFileError() << "No Les Houches file name. " << "Use 'set " << name() << ":FileName'." << Exception::runerror; cfile.open(filename()); if ( !cfile ) throw FxFxFileError() << "The FxFxFileReader '" << name() << "' could not open the " << "event file called '" << theFileName << "'." << Exception::runerror; cfile.readline(); if ( !cfile.find(" attributes = StringUtils::xmlAttributes("LesHouchesEvents", cfile.getline()); LHFVersion = attributes["version"]; //cout << LHFVersion << endl; if ( LHFVersion.empty() ) return; bool readingHeader = false; bool readingInit = false; headerBlock = ""; // char (cwgtinfo_weights_info[250][15]); string hs; // int cwgtinfo_nn(0); // Loop over all lines until we hit the tag. bool readingInitWeights = false, readingInitWeights_sc = false; string weightinfo; while ( cfile.readline() ) { // found the init block for multiple weights if(cfile.find(""; erase_substr(sub, str_arrow); scalename = sub; } ++ws; } while (isc); /* now get the relevant information * e.g. scales or PDF sets used */ string startDEL = ">"; //starting delimiter string stopDEL = ""; //end delimiter unsigned firstLim = hs.find(startDEL); //find start of delimiter // unsigned lastLim = hs.find(stopDEL); //find end of delimitr string scinfo = hs.substr(firstLim); //define the information for the scale erase_substr(scinfo,stopDEL); erase_substr(scinfo,startDEL); scinfo = StringUtils::stripws(scinfo); //cout << "scinfo = " << scinfo << endl; /* fill in the map * indicating the information to be appended to each scale * i.e. scinfo for each scalname */ scalemap[scalename] = scinfo.c_str(); string str_id = "id="; string str_prime = "'"; erase_substr(scalename, str_id); erase_substr(scalename, str_prime); optionalWeightsNames.push_back(scalename); } } if ( cfile.find("") ) { //cout << "found init block" << endl; // We have hit the init block, so we should expect to find the // standard information in the following. But first check for // attributes. initAttributes = StringUtils::xmlAttributes("init", cfile.getline()); readingInit = true; cfile.readline(); if ( !( cfile >> heprup.IDBMUP.first >> heprup.IDBMUP.second >> heprup.EBMUP.first >> heprup.EBMUP.second >> heprup.PDFGUP.first >> heprup.PDFGUP.second >> heprup.PDFSUP.first >> heprup.PDFSUP.second >> heprup.IDWTUP >> heprup.NPRUP ) ) { heprup.NPRUP = -42; LHFVersion = ""; return; } heprup.resize(); for ( int i = 0; i < heprup.NPRUP; ++i ) { cfile.readline(); if ( !( cfile >> heprup.XSECUP[i] >> heprup.XERRUP[i] >> heprup.XMAXUP[i] >> heprup.LPRUP[i] ) ) { heprup.NPRUP = -42; LHFVersion = ""; return; } } } if ( cfile.find("> sub; if(we==2) { npLO = atoi(sub.c_str()); } if(we==5) { npNLO = atoi(sub.c_str()); } ++we; } while (ievat); optionalnpLO = npLO; optionalnpNLO = npNLO; std::stringstream npstringstream; npstringstream << "np " << npLO << " " << npNLO; std::string npstrings = npstringstream.str(); /* the FxFx merging information * becomes part of the optionalWeights, labelled -999 * for future reference */ if(theIncludeFxFxTags) optionalWeights[npstrings.c_str()] = -999; string ecomstring = "ecom"; optionalWeights[ecomstring.c_str()] = heprup.EBMUP.first+heprup.EBMUP.second; if ( !cfile.readline() ) return false; // The first line determines how many subsequent particle lines we // have. if ( !( cfile >> hepeup.NUP >> hepeup.IDPRUP >> hepeup.XWGTUP >> hepeup.SCALUP >> hepeup.AQEDUP >> hepeup.AQCDUP ) ) return false; hepeup.resize(); // Read all particle lines. for ( int i = 0; i < hepeup.NUP; ++i ) { if ( !cfile.readline() ) return false; if ( !( cfile >> hepeup.IDUP[i] >> hepeup.ISTUP[i] >> hepeup.MOTHUP[i].first >> hepeup.MOTHUP[i].second >> hepeup.ICOLUP[i].first >> hepeup.ICOLUP[i].second >> hepeup.PUP[i][0] >> hepeup.PUP[i][1] >> hepeup.PUP[i][2] >> hepeup.PUP[i][3] >> hepeup.PUP[i][4] >> hepeup.VTIMUP[i] >> hepeup.SPINUP[i] ) ) return false; if(std::isnan(hepeup.PUP[i][0])||std::isnan(hepeup.PUP[i][1])|| std::isnan(hepeup.PUP[i][2])||std::isnan(hepeup.PUP[i][3])|| std::isnan(hepeup.PUP[i][4])) throw Exception() << "nan's as momenta in Les Houches file " << Exception::eventerror; if(hepeup.MOTHUP[i].first -1==i || hepeup.MOTHUP[i].second-1==i) { throw Exception() << "Particle has itself as a mother in Les Houches " << "file, this is not allowed\n" << Exception::eventerror; } } + LHEeventnum = -1; // set the LHEeventnum to -1, this will be the default if the tag is not found + + // Now read any additional comments and named weights. // read until the end of rwgt is found bool readingWeights = false, readingaMCFast = false, readingMG5ClusInfo = false; while ( cfile.readline() && !cfile.find("")) { if(cfile.find("> sub; if(wi==1) { string str_arrow = ">" ; erase_substr(sub, str_arrow); weightName = sub; } if(wi==2) weightValue = atof(sub.c_str()); ++wi; } while (iss); // store the optional weights found in the temporary map //cout << "weightName, weightValue= " << weightName << ", " << weightValue << endl; optionalWeightsTemp[weightName] = weightValue; } /* reading of aMCFast weights */ if(readingaMCFast) { std::stringstream amcfstringstream; amcfstringstream << "aMCFast " << cfile.getline(); std::string amcfstrings = amcfstringstream.str(); string str_newline = "\n"; erase_substr(amcfstrings,str_newline); optionalWeights[amcfstrings.c_str()] = -111; //for the aMCFast we give them a weight -111 for future identification } /* read additional MG5 Clustering information * used in LO merging */ if(readingMG5ClusInfo) { string hs = cfile.getline(); string startDEL = ""; //end delimiter unsigned firstLim = hs.find(startDEL); //find start of delimiter // unsigned lastLim = hs.find(stopDEL); //find end of delimitr string mg5clusinfo = hs.substr(firstLim); //define the information for the scale erase_substr(mg5clusinfo,stopDEL); erase_substr(mg5clusinfo,startDEL); string str_arrow = ">"; erase_substr(mg5clusinfo,str_arrow); string str_quotation = "\""; erase_substr(mg5clusinfo,str_quotation); string str_newline= "\n"; erase_substr(mg5clusinfo,str_newline); optionalWeights[mg5clusinfo.c_str()] = -222; //for the mg5 scale info weights we give them a weight -222 for future identification } + + //the event num tag + if(cfile.find("::const_iterator it=optionalWeightsTemp.begin(); it!=optionalWeightsTemp.end(); ++it){ string itfirst = it->first; erase_substr(itfirst, "'"); erase_substr(itfirst, "\""); for (map::const_iterator it2=scalemap.begin(); it2!=scalemap.end(); ++it2){ //find the scale id in the scale information and add this information string it2first = it2->first; erase_substr(it2first, "'"); erase_substr(it2first, "\""); //cout << "itfirst, it2first = " << itfirst << "\t" << it2first << endl; if(itfirst==it2first) { string info = it2->second; string str_newline = "\n"; erase_substr(info, str_newline); //set the optional weights optionalWeights[info] = it->second; } } } /* additionally, we set the "central" scale * this is actually the default event weight */ string central = "central"; if (theIncludeCentral) optionalWeights[central] = hepeup.XWGTUP; if ( !cfile ) return false; return true; } void FxFxFileReader::close() { cfile.close(); } void FxFxFileReader::persistentOutput(PersistentOStream & os) const { os << neve << LHFVersion << outsideBlock << headerBlock << initComments << initAttributes << eventComments << eventAttributes << theFileName << theQNumbers << theIncludeFxFxTags << theIncludeCentral << theDecayer; } void FxFxFileReader::persistentInput(PersistentIStream & is, int) { is >> neve >> LHFVersion >> outsideBlock >> headerBlock >> initComments >> initAttributes >> eventComments >> eventAttributes >> theFileName >> theQNumbers >> theIncludeFxFxTags >> theIncludeCentral >> theDecayer; ieve = 0; } ClassDescription FxFxFileReader::initFxFxFileReader; // Definition of the static class description member. void FxFxFileReader::Init() { static ClassDocumentation documentation ("ThePEG::FxFxFileReader is an base class to be used for objects " "which reads event files from matrix element generators. This class is " "able to read plain event files conforming to the Les Houches Event File " "accord."); static Parameter interfaceFileName ("FileName", "The name of a file containing events conforming to the Les Houches " "protocol to be read into ThePEG. A file name ending in " ".gz will be read from a pipe which uses " "zcat. If a file name ends in | the " "preceeding string is interpreted as a command, the output of which " "will be read through a pipe.", &FxFxFileReader::theFileName, "", false, false); interfaceFileName.fileType(); interfaceFileName.rank(11); static Switch interfaceQNumbers ("QNumbers", "Whether or not to read search for and read a QNUMBERS" " block in the header of the file.", &FxFxFileReader::theQNumbers, false, false, false); static SwitchOption interfaceQNumbersYes (interfaceQNumbers, "Yes", "Use QNUMBERS", true); static SwitchOption interfaceQNumbersNo (interfaceQNumbers, "No", "Don't use QNUMBERS", false); static Switch interfaceIncludeFxFxTags ("IncludeFxFxTags", "Include FxFx tags", &FxFxFileReader::theIncludeFxFxTags, true, true, false); static SwitchOption interfaceIncludeFxFxTagsYes (interfaceIncludeFxFxTags, "Yes", "Use the FxFx tags", true); static SwitchOption interfaceIncludeFxFxTagsNo (interfaceIncludeFxFxTags, "No", "Don't use the FxFx tags", false); static Switch interfaceIncludeCentral ("IncludeCentral", "Include definition of central weight", &FxFxFileReader::theIncludeCentral, false, true, false); static SwitchOption interfaceIncludeCentralYes (interfaceIncludeCentral, "Yes", "include definition of central weight", true); static SwitchOption interfaceIncludeCentralNo (interfaceIncludeCentral, "No", "Don't include definition of central weight", false); static Reference interfaceDecayer ("Decayer", "Decayer to use for any decays read from the QNUMBERS Blocks", &FxFxFileReader::theDecayer, false, false, true, true, false); } void FxFxFileReader::erase_substr(std::string& subject, const std::string& search) { size_t pos = 0; while((pos = subject.find(search, pos)) != std::string::npos) { subject.erase( pos, search.length() ); } } diff --git a/MatrixElement/FxFx/FxFxReader.cc b/MatrixElement/FxFx/FxFxReader.cc --- a/MatrixElement/FxFx/FxFxReader.cc +++ b/MatrixElement/FxFx/FxFxReader.cc @@ -1,1601 +1,1605 @@ // -*- C++ -*- // // FxFxReader.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2019 Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the FxFxReader class. // #include "FxFxReader.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/RefVector.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Command.h" //#include "config.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/PDF/PartonExtractor.h" #include "ThePEG/PDF/NoPDF.h" #include "ThePEG/Cuts/Cuts.h" #include "ThePEG/EventRecord/TmpTransform.h" #include "ThePEG/Utilities/UtilityBase.h" #include "ThePEG/Handlers/XComb.h" #include "ThePEG/Handlers/CascadeHandler.h" #include "FxFxEventHandler.h" #include "ThePEG/Utilities/Throw.h" #include "ThePEG/Utilities/HoldFlag.h" #include "ThePEG/Utilities/Debug.h" #include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h" using namespace ThePEG; FxFxReader::FxFxReader(bool active) : theNEvents(0), position(0), reopened(0), theMaxScan(-1), scanning(false), isActive(active), theCacheFileName(""), doCutEarly(true), preweight(1.0), reweightPDF(false), doInitPDFs(false), theMaxMultCKKW(0), theMinMultCKKW(0), lastweight(1.0), maxFactor(1.0), optionalnpLO(0), optionalnpNLO(0), weightScale(1.0*picobarn), skipping(false), theMomentumTreatment(0), useWeightWarnings(true),theReOpenAllowed(true), theIncludeSpin(true) {} FxFxReader::FxFxReader(const FxFxReader & x) : HandlerBase(x), LastXCombInfo<>(x), heprup(x.heprup), hepeup(x.hepeup), inData(x.inData), inPDF(x.inPDF), outPDF(x.outPDF), thePartonExtractor(x.thePartonExtractor), thePartonBins(x.thePartonBins), theXCombs(x.theXCombs), theCuts(x.theCuts), theNEvents(x.theNEvents), position(x.position), reopened(x.reopened), theMaxScan(x.theMaxScan), scanning(false), isActive(x.isActive), theCacheFileName(x.theCacheFileName), doCutEarly(x.doCutEarly), stats(x.stats), statmap(x.statmap), thePartonBinInstances(x.thePartonBinInstances), reweights(x.reweights), preweights(x.preweights), preweight(x.preweight), reweightPDF(x.reweightPDF), doInitPDFs(x.doInitPDFs), theMaxMultCKKW(x.theMaxMultCKKW), theMinMultCKKW(x.theMinMultCKKW), lastweight(x.lastweight), maxFactor(x.maxFactor), weightScale(x.weightScale), xSecWeights(x.xSecWeights), maxWeights(x.maxWeights), skipping(x.skipping), theMomentumTreatment(x.theMomentumTreatment), useWeightWarnings(x.useWeightWarnings), theReOpenAllowed(x.theReOpenAllowed), theIncludeSpin(x.theIncludeSpin) {} FxFxReader::~FxFxReader() {} void FxFxReader::doinitrun() { HandlerBase::doinitrun(); stats.reset(); for ( StatMap::iterator i = statmap.begin(); i != statmap.end(); ++i ) i->second.reset(); open(); if ( cacheFileName().length() ) openReadCacheFile(); position = 0; reopened = 0; } bool FxFxReader::preInitialize() const { if ( HandlerBase::preInitialize() ) return true; if ( doInitPDFs && ! ( inPDF.first && inPDF.second ) ) return true; return false; } void FxFxReader::doinit() { HandlerBase::doinit(); open(); close(); if ( !heprup.IDBMUP.first || !heprup.IDBMUP.second ) Throw() << "No information about incoming particles were found in " << "FxFxReader '" << name() << "'." << Exception::warning; inData = make_pair(getParticleData(heprup.IDBMUP.first), getParticleData(heprup.IDBMUP.second)); if ( heprup.EBMUP.first <= 0.0 || heprup.EBMUP.second <= 0.0 ) Throw() << "No information about the energy of incoming particles were found in " << "FxFxReader '" << name() << "'." << Exception::warning; if ( doInitPDFs && ! ( inPDF.first && inPDF.second ) ) { initPDFs(); if ( ! ( inPDF.first && inPDF.second ) ) Throw() << "FxFxReader '" << name() << "' could not create PDFBase objects in pre-initialization." << Exception::warning; } else if ( !inPDF.first || !inPDF.second ) Throw() << "No information about the PDFs of incoming particles were found in " << "FxFxReader '" << name() << "'." << Exception::warning; } void FxFxReader::initPDFs() { if ( inPDF.first && inPDF.second ) return; string remhname; if ( heprup.PDFSUP.first && !inPDF.first) { inPDF.first = dynamic_ptr_cast (generator()->preinitCreate("ThePEG::LHAPDF", fullName() + "/PDFA", "ThePEGLHAPDF.so")); if ( !inPDF.first ) { Throw() << "FxFxReader '" << name() << "' could not use information " << "about the PDFs used because the LHAPDF library was not properly " "defined." << Exception::warning; return; } remhname = fullName() + "/DummyRemH"; generator()->preinitCreate("ThePEG::NoRemnants", remhname); generator()->preinitInterface(inPDF.first, "RemnantHandler", "set", remhname); if ( heprup.PDFGUP.first > 0 && heprup.PDFGUP.first < 10 ) { ostringstream os; os << heprup.PDFGUP.first << " " << heprup.PDFSUP.first; generator()->preinitInterface(inPDF.first, "PDFLIBNumbers", "set", os.str()); } else { ostringstream os; os << heprup.PDFGUP.first*1000 + heprup.PDFSUP.first; generator()->preinitInterface(inPDF.first, "PDFNumber", "set", os.str()); } generator()->preinitInterface(inPDF.first, "RangeException", "newdef", "Freeze"); } if ( heprup.PDFSUP.second && !inPDF.second) { inPDF.second = dynamic_ptr_cast (generator()->preinitCreate("ThePEG::LHAPDF", fullName() + "/PDFB", "ThePEGLHAPDF.so")); if ( !inPDF.second ) { Throw() << "FxFxReader '" << name() << "' could not use information " << "about the PDFs used because the LHAPDF library was not properly " "defined." << Exception::warning; return; } if ( remhname == "" ) { remhname = fullName() + "/DummyRemH"; generator()->preinitCreate("ThePEG::NoRemnants", remhname); } generator()->preinitInterface(inPDF.second, "RemnantHandler", "set", remhname); if ( heprup.PDFGUP.second > 0 && heprup.PDFGUP.second < 10 ) { ostringstream os; os << heprup.PDFGUP.second << " " << heprup.PDFSUP.second; generator()->preinitInterface(inPDF.second, "PDFLIBNumbers", "set", os.str()); } else { ostringstream os; os << heprup.PDFGUP.second*1000 + heprup.PDFSUP.second; generator()->preinitInterface(inPDF.second, "PDFNumber", "set", os.str()); } generator()->preinitInterface(inPDF.second, "RangeException", "newdef", "Freeze"); } if ( ! ( inPDF.first && inPDF.second ) ) Throw() << "FxFxReader '" << name() << "' could not find information about the PDFs used." << Exception::warning; } void FxFxReader::initialize(FxFxEventHandler & eh) { Energy2 Smax = ZERO; double Y = 0.0; if ( !theCuts ) { theCuts = eh.cuts(); if ( !theCuts ) Throw() << "No Cuts object was assigned to the FxFxReader '" << name() << "' nor was one\nassigned to the controlling " << "FxFxEventHandler '" << eh.name() << "'.\nAt least one of them " << "needs to have a Cuts object." << Exception::runerror; Smax = cuts().SMax(); Y = cuts().Y(); } theCKKW = eh.CKKWHandler(); if ( !partonExtractor() ) { thePartonExtractor = eh.partonExtractor(); if ( !partonExtractor() ) Throw() << "No PartonExtractor object was assigned to the FxFxReader '" << name() << "' nor was one\nassigned to the controlling " << "FxFxEventHandler '" << eh.name() << "'.\nAt least one of them " << "needs to have a PartonExtractor object." << Exception::runerror; } open(); Energy emax = 2.0*sqrt(heprup.EBMUP.first*heprup.EBMUP.second)*GeV; theCuts->initialize(sqr(emax), 0.5*log(heprup.EBMUP.first/heprup.EBMUP.second)); if ( Smax > ZERO && ( Smax != cuts().SMax() || Y != cuts().Y() ) ) Throw() << "The FxFxReader '" << name() << "' uses the same Cuts object " << "as another FxFxReader which has not got the same energies of " << "the colliding particles. For the generation to work properly " << "different FxFxReader object with different colliding particles " << "must be assigned different (although possibly identical) Cuts " << "objects." << Exception::warning; thePartonBins = partonExtractor()->getPartons(emax, inData, cuts()); for ( int i = 0, N = partonBins().size(); i < N; ++i ) { theXCombs[partonBins()[i]] = new_ptr(XComb(emax, inData, &eh, partonExtractor(), CKKWHandler(), partonBins()[i], theCuts)); partonExtractor()->nDims(partonBins()[i]); } outPDF = make_pair(partonExtractor()->getPDF(inData.first), partonExtractor()->getPDF(inData.second)); close(); if ( !heprup.IDWTUP && useWeightWarnings ) Throw() << "No information about the weighting scheme was found. The events " << "produced by FxFxReader " << name() << " may not be sampled correctly." << Exception::warning; if ( abs(heprup.IDWTUP) > 1 && !eh.weighted() && useWeightWarnings ) Throw() << "FxFxReader " << name() << " has the IDWTUP flag set to " << heprup.IDWTUP << " which is not supported by this reader, the " << "produced events may not be sampled correctly. It is up to " << "sub-classes of FxFxReader to correctly convert to match IDWTUP " << "+/- 1. Will try to make intelligent guesses to get " << "correct statistics.\nIn most cases this should be sufficient. " << "Unset WeightWarnings to avoid this message," << "or set Weighted to on." << Exception::warning; if ( heprup.IDWTUP != eh.weightOption() && abs(heprup.IDWTUP) < 3 && useWeightWarnings ) Throw() << "FxFxReader " << name() << " has the IDWTUP flag set to " << heprup.IDWTUP << ", which does not correspond\nto the weight option " << eh.weightOption() << " set in " << "the FxFxEventHandler " << eh.name() << ".\n\n" << "Use the following handler setting instead:\n" << " set " << eh.name() << ":WeightOption " << heprup.IDWTUP << "\nWill try to make intelligent guesses to get " << "correct statistics. In most cases this should be sufficient. " << "Unset WeightWarnings to avoid this message" << Exception::warning; scan(); initStat(); } long FxFxReader::scan() { open(); // Shall we write the events to a cache file for fast reading? If so // we write to a temporary file if the caches events should be // randomized. if ( cacheFileName().length() ) openWriteCacheFile(); // Keep track of the number of events scanned. long neve = 0; long cuteve = 0; bool negw = false; // If the open() has not already gotten information about subprocesses // and cross sections we have to scan through the events. if ( !heprup.NPRUP || cacheFile() || abs(heprup.IDWTUP) != 1 ) { // why scan if IDWTUP != 1? HoldFlag<> isScanning(scanning); double oldsum = 0.0; vector lprup; vector newmax; vector oldeve; vector neweve; vector sumlprup; vector sumsqlprup; vector nscanned; for ( int i = 0; ( maxScan() < 0 || i < maxScan() ) && readEvent(); ++i ) { if ( !checkPartonBin() ) Throw() << "Found event in LesHouchesReader '" << name() << "' which cannot be handeled by the assigned PartonExtractor '" << partonExtractor()->name() << "'." << Exception::runerror; vector::iterator idit = find(lprup.begin(), lprup.end(), hepeup.IDPRUP); int id = lprup.size(); if ( idit == lprup.end() ) { lprup.push_back(hepeup.IDPRUP); newmax.push_back(0.0); neweve.push_back(0); oldeve.push_back(0); sumlprup.push_back(0.); sumsqlprup.push_back(0.); nscanned.push_back(0); } else { id = idit - lprup.begin(); } ++neve; ++oldeve[id]; oldsum += hepeup.XWGTUP; sumlprup[id] += hepeup.XWGTUP; sumsqlprup[id] += sqr(hepeup.XWGTUP); ++nscanned[id]; if ( cacheFile() ) { if ( eventWeight() == 0.0 ) { ++cuteve; continue; } cacheEvent(); } ++neweve[id]; newmax[id] = max(newmax[id], abs(eventWeight())); if ( eventWeight() < 0.0 ) negw = true; } //end of scanning events xSecWeights.resize(oldeve.size(), 1.0); for ( int i = 0, N = oldeve.size(); i < N; ++i ) if ( oldeve[i] ) xSecWeights[i] = double(neweve[i])/double(oldeve[i]); if ( maxScan() < 0 || neve > NEvents() ) NEvents(neve - cuteve); if ( lprup.size() == heprup.LPRUP.size() ) { for ( int id = 0, N = lprup.size(); id < N; ++id ) { vector::iterator idit = find(heprup.LPRUP.begin(), heprup.LPRUP.end(), hepeup.IDPRUP); if ( idit == heprup.LPRUP.end() ) { Throw() << "When scanning events, the LesHouschesReader '" << name() << "' found undeclared processes." << Exception::warning; heprup.NPRUP = 0; break; } int idh = idit - heprup.LPRUP.begin(); heprup.XMAXUP[idh] = newmax[id]; } } if ( heprup.NPRUP == 0 ) { // No heprup block was supplied or something went wrong. heprup.NPRUP = lprup.size(); heprup.LPRUP.resize(lprup.size()); heprup.XMAXUP.resize(lprup.size()); for ( int id = 0, N = lprup.size(); id < N; ++id ) { heprup.LPRUP[id] = lprup[id]; heprup.XMAXUP[id] = newmax[id]; } } if ( abs(heprup.IDWTUP) != 1 ) { // Try to fix things if abs(heprup.IDWTUP) != 1. double sumxsec = 0.0; if(abs(heprup.IDWTUP)==3) { for ( int id = 0; id < heprup.NPRUP; ++id ) sumxsec += heprup.XSECUP[id]; } else { for ( int id = 0; id < heprup.NPRUP; ++id ) { //set the cross section directly from the event weights read heprup.XSECUP[id] = sumlprup[id]/nscanned[id]; heprup.XERRUP[id] = (sumsqlprup[id]/nscanned[id] - sqr(sumlprup[id]/nscanned[id])) / nscanned[id]; if(fabs(heprup.XERRUP[id]) < 1e-10) { heprup.XERRUP[id] = 0; } if(fabs(newmax[id]) < 1e-10) { newmax[id] = 0; } if(heprup.XERRUP[id] < 0.) { if( heprup.XERRUP[id]/(sumsqlprup[id]/nscanned[id])>-1e-10) heprup.XERRUP[id] = 0.; else { Throw() << "Negative error when scanning events in LesHouschesReader '" << name() << Exception::warning; heprup.XERRUP[id] = 0.; } } heprup.XERRUP[id] = sqrt( heprup.XERRUP[id] ); heprup.XMAXUP[id] = newmax[id]; //cout << "heprup.XMAXUP[id] = " << heprup.XMAXUP[id] << endl; sumxsec += heprup.XSECUP[id]; } } //cout << "sumxsec = " << sumxsec << endl; //weightScale = picobarn*neve*sumxsec/oldsum; weightScale = 1.0*picobarn; // temporary fix? // cout << "weightscale = " << weightScale/picobarn << endl; } } if ( cacheFile() ) closeCacheFile(); if ( negw ) heprup.IDWTUP = min(-abs(heprup.IDWTUP), -1); return neve; } void FxFxReader::setWeightScale(long) {} void FxFxReader::initStat() { stats.reset(); statmap.clear(); if ( heprup.NPRUP <= 0 ) return; double sumx = 0.0; xSecWeights.resize(heprup.NPRUP, 1.0); maxWeights.clear(); for ( int ip = 0; ip < heprup.NPRUP; ++ip ) { sumx = max(heprup.XMAXUP[ip]*xSecWeights[ip], sumx); statmap[heprup.LPRUP[ip]] = XSecStat(heprup.XMAXUP[ip]*weightScale*xSecWeights[ip]); maxWeights[heprup.LPRUP[ip]] = heprup.XMAXUP[ip]; } stats.maxXSec(sumx*weightScale); maxFactor = 1.0; } void FxFxReader::increaseMaxXSec(CrossSection maxxsec) { for ( int i = 0; i < heprup.NPRUP; ++i ) statmap[heprup.LPRUP[i]].maxXSec(statmap[heprup.LPRUP[i]].maxXSec()* maxxsec/stats.maxXSec()); maxFactor *= maxxsec/stats.maxXSec(); stats.maxXSec(maxxsec); } tXCombPtr FxFxReader::getXComb() { if ( lastXCombPtr() ) return lastXCombPtr(); fillEvent(); connectMothers(); tcPBPair sel = createPartonBinInstances(); tXCombPtr lastXC = xCombs()[sel]; // clean up the old XComb object before switching to a new one if ( theLastXComb && theLastXComb != lastXC ) theLastXComb->clean(); theLastXComb = lastXC; lastXCombPtr()->subProcess(SubProPtr()); lastXCombPtr()->setPartonBinInstances(partonBinInstances(), sqr(hepeup.SCALUP)*GeV2); lastXCombPtr()->lastAlphaS(hepeup.AQCDUP); lastXCombPtr()->lastAlphaEM(hepeup.AQEDUP); return lastXCombPtr(); } tSubProPtr FxFxReader::getSubProcess() { getXComb(); if ( subProcess() ) return subProcess(); lastXCombPtr()->subProcess(new_ptr(SubProcess(lastPartons(), tCollPtr(), this))); subProcess()->setOutgoing(outgoing().begin(), outgoing().end()); subProcess()->setIntermediates(intermediates().begin(), intermediates().end()); return subProcess(); } void FxFxReader::fillEvent() { if ( !particleIndex.empty() ) return; particleIndex.clear(); colourIndex.clear(); colourIndex(0, tColinePtr()); createParticles(); createBeams(); } void FxFxReader::reopen() { // If we didn't know how many events there were, we know now. if ( NEvents() <= 0 ) NEvents(position); ++reopened; // How large fraction of the events have we actually used? And how // large will we have used if we go through the file again? double frac = double(stats.attempts())/double(NEvents()); if ( frac*double(reopened + 1)/double(reopened) > 1.0 && NEvents() - stats.attempts() < generator()->N() - generator()->currentEventNumber() ) { if(theReOpenAllowed) generator()->logWarning(FxFxReopenWarning() << "Reopening FxFxReader '" << name() << "' after accessing " << stats.attempts() << " events out of " << NEvents() << Exception::warning); else throw FxFxReopenWarning() << "More events requested than available in FxFxReader " << name() << Exception::runerror; } if ( cacheFile() ) { closeCacheFile(); openReadCacheFile(); if ( !uncacheEvent() ) Throw() << "Could not reopen FxFxReader '" << name() << "'." << Exception::runerror; } else { close(); open(); if ( !readEvent() ) Throw() << "Could not reopen FxFxReader '" << name() << "'." << Exception::runerror; } } void FxFxReader::reset() { particleIndex.clear(); colourIndex.clear(); if ( theLastXComb ) theLastXComb->clean(); theLastXComb = tXCombPtr(); } bool FxFxReader::readEvent() { reset(); if ( !doReadEvent() ) return false; // If we are just skipping event we do not need to reweight or do // anything fancy. if ( skipping ) return true; if ( cacheFile() && !scanning ) return true; // Reweight according to the re- and pre-weights objects in the // FxFxReader base class. lastweight = reweight(); if ( !reweightPDF && !cutEarly() ) return true; // We should try to reweight the PDFs or make early cuts here. fillEvent(); double x1 = incoming().first->momentum().plus()/ beams().first->momentum().plus(); if ( reweightPDF && inPDF.first && outPDF.first && inPDF.first != outPDF.first ) { if ( hepeup.XPDWUP.first <= 0.0 ) hepeup.XPDWUP.first = inPDF.first->xfx(inData.first, incoming().first->dataPtr(), sqr(hepeup.SCALUP*GeV), x1); double xf = outPDF.first->xfx(inData.first, incoming().first->dataPtr(), sqr(hepeup.SCALUP*GeV), x1); lastweight *= xf/hepeup.XPDWUP.first; hepeup.XPDWUP.first = xf; } double x2 = incoming().second->momentum().minus()/ beams().second->momentum().minus(); if ( reweightPDF && inPDF.second && outPDF.second && inPDF.second != outPDF.second ) { if ( hepeup.XPDWUP.second <= 0.0 ) hepeup.XPDWUP.second = inPDF.second->xfx(inData.second, incoming().second->dataPtr(), sqr(hepeup.SCALUP*GeV), x2); double xf = outPDF.second->xfx(inData.second, incoming().second->dataPtr(), sqr(hepeup.SCALUP*GeV), x2); lastweight *= xf/hepeup.XPDWUP.second; hepeup.XPDWUP.second = xf; } if ( cutEarly() ) { if ( !cuts().initSubProcess((incoming().first->momentum() + incoming().second->momentum()).m2(), 0.5*log(x1/x2)) ) lastweight = 0.0; tSubProPtr sub = getSubProcess(); TmpTransform tmp(sub, Utilities::getBoostToCM(sub->incoming())); if ( !cuts().passCuts(*sub) ) lastweight = 0.0; } return true; } double FxFxReader::getEvent() { if ( cacheFile() ) { if ( !uncacheEvent() ) reopen(); } else { if ( !readEvent() ) reopen(); } ++position; double max = maxWeights[hepeup.IDPRUP]*maxFactor; // cout << "maxFactor = " << maxFactor << " maxWeights[hepeup.IDPRUP] = " << maxWeights[hepeup.IDPRUP] << endl; // normalize all the weights to the max weight for(map::iterator it=optionalWeights.begin(); it!=optionalWeights.end();++it) { if(it->first!="ecom" && it->second!=-999 && it->second!=-111 && it->second!=-222 && it->second!=-333) { it->second = (max != 0.0) ? it->second/max : 0.0; } } //cout << "hepeup.XWGTUP = " << hepeup.XWGTUP << endl; //cout << "max = " << max << " " << eventWeight() << endl; return max != 0.0? eventWeight()/max: 0.0; } void FxFxReader::skip(long n) { HoldFlag<> skipflag(skipping); while ( n-- ) getEvent(); } double FxFxReader::reweight() { preweight = 1.0; if ( reweights.empty() && preweights.empty() && !( CKKWHandler() && maxMultCKKW() > 0 && maxMultCKKW() > minMultCKKW() ) ) return 1.0; fillEvent(); getSubProcess(); for ( int i = 0, N = preweights.size(); i < N; ++i ) { preweights[i]->setXComb(lastXCombPtr()); preweight *= preweights[i]->weight(); } double weight = preweight; for ( int i = 0, N = reweights.size(); i < N; ++i ) { reweights[i]->setXComb(lastXCombPtr()); weight *= reweights[i]->weight(); } // If we are caching events we do not want to do CKKW reweighting. if ( cacheFile() ) return weight; if ( CKKWHandler() && maxMultCKKW() > 0 && maxMultCKKW() > minMultCKKW() ) { CKKWHandler()->setXComb(lastXCombPtr()); weight *= CKKWHandler()->reweightCKKW(minMultCKKW(), maxMultCKKW()); } return weight; } bool FxFxReader::checkPartonBin() { // First find the positions of the incoming partons. pair< vector, vector > inc; for ( int i = 0; i < hepeup.NUP; ++i ) { if ( hepeup.ISTUP[i] == -9 ) { if ( inc.first.empty() ) inc.first.push_back(i); else if ( inc.second.empty() ) inc.second.push_back(i); } else if ( hepeup.ISTUP[i] == -1 ) { if ( inc.first.size() && hepeup.MOTHUP[i].first == inc.first.back() + 1 ) inc.first.push_back(i); else if ( inc.second.size() && hepeup.MOTHUP[i].first == inc.second.back() + 1 ) inc.second.push_back(i); else if ( inc.first.empty() ) { inc.first.push_back(-1); inc.first.push_back(i); } else if ( inc.second.empty() ) { inc.second.push_back(-1); inc.second.push_back(i); } else if ( inc.first.size() <= inc.second.size() ) inc.first.push_back(i); else inc.second.push_back(i); } } // Now store the corresponding id numbers pair< vector, vector > ids; ids.first.push_back(inc.first[0] < 0? heprup.IDBMUP.first: hepeup.IDUP[inc.first[0]]); for ( int i = 1, N = inc.first.size(); i < N; ++i ) ids.first.push_back(hepeup.IDUP[inc.first[i]]); ids.second.push_back(inc.second[0] < 0? heprup.IDBMUP.second: hepeup.IDUP[inc.second[0]]); for ( int i = 1, N = inc.second.size(); i < N; ++i ) ids.second.push_back(hepeup.IDUP[inc.second[i]]); // Find the correct pair of parton bins. PBPair pbp; for ( int i = 0, N = partonBins().size(); i < N; ++i ) { tcPBPtr curr = partonBins()[i].first; int icurr = inc.first.size() - 1; while ( curr && icurr >= 0 ) { if ( curr->parton()->id () != ids.first[icurr] ) break; curr = curr->incoming(); --icurr; } if(!(!partonBins()[i].first->incoming() && !partonBins()[i].first->particle() && partonBins()[i].first->parton()->id () == ids.first[0] && ( inc.first.size()==1 || (inc.first.size()==2 && ids.first[0]==ids.first[1]))) && ( curr || icurr >= 0 ) ) continue; curr = partonBins()[i].second; icurr = inc.second.size() - 1; while ( curr && icurr >= 0 ) { if ( curr->parton()->id () != ids.second[icurr] ) break; curr = curr->incoming(); --icurr; } if(!(!partonBins()[i].second->incoming() && !partonBins()[i].second->particle() && partonBins()[i].second->parton()->id () == ids.second[0] && ( inc.second.size()==1 || (inc.second.size()==2 && ids.second[0]==ids.second[1]))) && ( curr || icurr >= 0 ) ) continue; pbp = partonBins()[i]; } // If we are only checking we return here. return ( pbp.first && pbp.second ); } namespace { bool recursionNotNull(tcPBPtr bin, tcPPtr p) { while ( bin && p ) { if ( p->dataPtr() != bin->parton() ) break; bin = bin->incoming(); p = p->parents().size()? p->parents()[0]: tPPtr(); } return bin || p; } } tcPBPair FxFxReader::createPartonBinInstances() { tcPBPair sel; for ( int i = 0, N = partonBins().size(); i < N; ++i ) { tcPBPtr bin = partonBins()[i].first; tcPPtr p = incoming().first; if ( recursionNotNull(bin,p) ) continue; bin = partonBins()[i].second; p = incoming().second; if ( recursionNotNull(bin,p) ) continue; sel = partonBins()[i]; break; } if ( !sel.first || !sel.second ) Throw() << "Could not find appropriate PartonBin objects for event produced by " << "FxFxReader '" << name() << "'." << Exception::runerror; Direction<0> dir(true); thePartonBinInstances.first = new_ptr(PartonBinInstance(incoming().first, sel.first, -sqr(hepeup.SCALUP*GeV))); if ( thePartonBinInstances.first->xi() > 1.00001 ) { Throw() << "Found an event with momentum fraction larger than unity (x1=" << thePartonBinInstances.first->xi() << "). The event will be skipped." << Exception::warning; throw Veto(); } dir.reverse(); thePartonBinInstances.second = new_ptr(PartonBinInstance(incoming().second, sel.second, -sqr(hepeup.SCALUP*GeV))); if ( thePartonBinInstances.second->xi() > 1.00001 ) { Throw() << "Found an event with momentum fraction larger than unity (x2=" << thePartonBinInstances.second->xi() << "). The event will be skipped." << Exception::warning; throw Veto(); } return sel; } void FxFxReader::createParticles() { theBeams = PPair(); theIncoming = PPair(); theOutgoing = PVector(); theIntermediates = PVector(); for ( int i = 0, N = hepeup.IDUP.size(); i < N; ++i ) { if ( !hepeup.IDUP[i] ) continue; Lorentz5Momentum mom(hepeup.PUP[i][0]*GeV, hepeup.PUP[i][1]*GeV, hepeup.PUP[i][2]*GeV, hepeup.PUP[i][3]*GeV, hepeup.PUP[i][4]*GeV); if(theMomentumTreatment == 1) mom.rescaleEnergy(); else if(theMomentumTreatment == 2) mom.rescaleMass(); PDPtr pd = getParticleData(hepeup.IDUP[i]); if (!pd) { Throw() << "FxFxReader '" << name() << "' found unknown particle ID " << hepeup.IDUP[i] << " in Les Houches common block structure.\n" << "You need to define the new particle in an input file.\n" << Exception::runerror; } if ( ! pd->coloured() && ( hepeup.ICOLUP[i].first != 0 || hepeup.ICOLUP[i].second != 0 ) ) { Throw() << "FxFxReader " << name() << ": " << pd->PDGName() << " is not a coloured particle.\nIt should not have " << "(anti-)colour lines " << hepeup.ICOLUP[i].first << ' ' << hepeup.ICOLUP[i].second << " set; the event file needs to be fixed." << Exception::runerror; } PPtr p = pd->produceParticle(mom); if(hepeup.ICOLUP[i].first>=0 && hepeup.ICOLUP[i].second >=0) { tColinePtr c = colourIndex(hepeup.ICOLUP[i].first); if ( c ) c->addColoured(p); c = colourIndex(hepeup.ICOLUP[i].second); if ( c ) c->addAntiColoured(p); } else { tColinePtr c1 = colourIndex(abs(hepeup.ICOLUP[i].first )); tColinePtr c2 = colourIndex(abs(hepeup.ICOLUP[i].second)); if(pd->hasColour()) { c1->addColouredIndexed(p,1); c2->addColouredIndexed(p,2); } else { c1->addAntiColouredIndexed(p,1); c2->addAntiColouredIndexed(p,2); } } particleIndex(i + 1, p); switch ( hepeup.ISTUP[i] ) { case -9: if ( !theBeams.first ) theBeams.first = p; else if ( !theBeams.second ) theBeams.second = p; else Throw() << "To many incoming beam particles in the FxFxReader '" << name() << "'." << Exception::runerror; break; case -1: if ( !theIncoming.first ) theIncoming.first = p; else if ( !theIncoming.second ) theIncoming.second = p; else if ( particleIndex(theIncoming.first) == hepeup.MOTHUP[i].first ) theIncoming.first = p; else if ( particleIndex(theIncoming.second) == hepeup.MOTHUP[i].first ) theIncoming.second = p; else Throw() << "To many incoming particles to hard subprocess in the " << "FxFxReader '" << name() << "'." << Exception::runerror; p->scale(sqr(hepeup.SCALUP*GeV)); break; case 1: theOutgoing.push_back(p); p->scale(sqr(hepeup.SCALUP*GeV)); break; case -2: case 2: case 3: theIntermediates.push_back(p); break; default: Throw() << "Unknown status code (" << hepeup.ISTUP[i] << ") in the FxFxReader '" << name() << "'." << Exception::runerror; } // value 9 is defined as "Unknown or unpolarized particles" double spinup = hepeup.SPINUP[i]; if ( abs(spinup - 9) < 1.0e-3 ) spinup = 0.; if ( spinup < -1. || spinup > 1. ) { Throw() << "Polarization must be between -1 and 1, not " << spinup << " as found in the " << "FxFx event file.\nThe event file needs to be fixed." << Exception::runerror; } if( theIncludeSpin && abs(pd->id()) == ParticleID::tauminus && spinup !=0) { if(pd->iSpin() == PDT::Spin1Half ) { vector wave; Helicity::SpinorWaveFunction(wave,p,Helicity::outgoing,true); RhoDMatrix rho(pd->iSpin(),true); rho(0,0) = 0.5*(1.-spinup); rho(1,1) = 0.5*(1.+spinup); p->spinInfo()->rhoMatrix() = rho; p->spinInfo()-> DMatrix() = rho; } } } // check the colour flows, and if necessary create any sources/sinks // hard process // get the particles in the hard process PVector external; for ( int i = 0, N = hepeup.IDUP.size(); i < N; ++i ) { unsigned int moth; switch ( hepeup.ISTUP[i] ) { case -1: external.push_back(particleIndex.find(i+1)); break; case 1: case 2: case 3: moth = hepeup.MOTHUP[i].first; if(moth!=0 && (hepeup.ISTUP[moth]==-1||hepeup.ISTUP[moth]==-2|| hepeup.ISTUP[moth]==-9)) external.push_back(particleIndex.find(i+1)); moth = hepeup.MOTHUP[i].second; if(moth!=0 && (hepeup.ISTUP[moth]==-1||hepeup.ISTUP[moth]==-2|| hepeup.ISTUP[moth]==-9)) external.push_back(particleIndex.find(i+1)); break; case -2: case -9: default: break; } } // check the incoming/outgoing lines match vector unMatchedColour,unMatchedAntiColour; for(unsigned int ix=0;ix col = external[ix]->colourInfo()-> colourLines(); vector anti = external[ix]->colourInfo()->antiColourLines(); if(hepeup.ISTUP[particleIndex(external[ix])-1]<0) swap(col,anti); if(!col.empty()) { for(unsigned int ic1=0;ic1 col2; if(hepeup.ISTUP[particleIndex(external[iy])-1]<0) { if(external[iy]->colourInfo()->colourLines().empty()) continue; col2 = external[iy]->colourInfo()->colourLines(); } else if(hepeup.ISTUP[particleIndex(external[iy])-1]>0) { if(external[iy]->colourInfo()->antiColourLines().empty()) continue; col2 = external[iy]->colourInfo()->antiColourLines(); } for(unsigned int ic2=0;ic2(col[ic1])); } } if(!anti.empty()) { - for(unsigned int ic1=0;ic1 anti2; if(hepeup.ISTUP[particleIndex(external[iy])-1]<0) { - if(external[iy]->colourInfo()->colourLines().empty()) continue; + if(external[iy]->colourInfo()->antiColourLines().empty()) continue; anti2 = external[iy]->colourInfo()->antiColourLines(); } else if(hepeup.ISTUP[particleIndex(external[iy])-1]>0) { - if(external[iy]->colourInfo()->antiColourLines().empty()) continue; + if(external[iy]->colourInfo()->colourLines().empty()) continue; anti2 = external[iy]->colourInfo()->colourLines(); } for(unsigned int ic2=0;ic2(anti[ic1])); } } } // might have source/sink if( unMatchedColour.size() + unMatchedAntiColour.size() != 0) { if(unMatchedColour.size() == 3 ) { unMatchedColour[0]->setSourceNeighbours(unMatchedColour[1], unMatchedColour[2]); } else if(unMatchedColour.size() != 0 && ThePEG_DEBUG_LEVEL) { Throw() << "FxFxReader '" << name() << "' found inconsistent colour " << "flow in Les Houches common block structure for hard process.\n" << hepeup << Exception::runerror; } if(unMatchedAntiColour.size() == 3 ) { unMatchedAntiColour[0]->setSinkNeighbours(unMatchedAntiColour[1], unMatchedAntiColour[2]); } else if(unMatchedAntiColour.size() != 0 && ThePEG_DEBUG_LEVEL) { Throw() << "FxFxReader '" << name() << "' found inconsistent colour " << "flow in Les Houches common block structure for hard process.\n" << hepeup << Exception::runerror; } } // any subsequent decays for ( int i = 0, N = hepeup.IDUP.size(); i < N; ++i ) { if(hepeup.ISTUP[i] !=2 && hepeup.ISTUP[i] !=3) continue; PVector external; external.push_back(particleIndex.find(i+1)); for ( int j = 0; j < N; ++j ) { if(hepeup.MOTHUP[j].first==i+1|| hepeup.MOTHUP[j].second==i+1) external.push_back(particleIndex.find(j+1)); } // check the incoming/outgoing lines match vector unMatchedColour,unMatchedAntiColour; for(unsigned int ix=0;ix col = external[ix]->colourInfo()-> colourLines(); vector anti = external[ix]->colourInfo()->antiColourLines(); if(ix==0) swap(col,anti); if(!col.empty()) { for(unsigned int ic1=0;ic1 col2; if(iy==0) { if(external[iy]->colourInfo()->colourLines().empty()) continue; col2 = external[iy]->colourInfo()->colourLines(); } else { if(external[iy]->colourInfo()->antiColourLines().empty()) continue; col2 = external[iy]->colourInfo()->antiColourLines(); } for(unsigned int ic2=0;ic2(col[ic1])); } } if(!anti.empty()) { for(unsigned int ic1=0;ic1 anti2; if(iy==0) { if(external[iy]->colourInfo()->antiColourLines().empty()) continue; anti2 = external[iy]->colourInfo()->antiColourLines(); } else { if(external[iy]->colourInfo()->colourLines().empty()) continue; anti2 = external[iy]->colourInfo()->colourLines(); } for(unsigned int ic2=0;ic2(anti[ic1])); } } } // might have source/sink if( unMatchedColour.size() + unMatchedAntiColour.size() != 0) { if(unMatchedColour.size() == 3 ) { unMatchedColour[0]->setSourceNeighbours(unMatchedColour[1], unMatchedColour[2]); } else if(unMatchedColour.size() != 0 && ThePEG_DEBUG_LEVEL) { Throw() << "FxFxReader '" << name() << "' found inconsistent colour " << "flow in Les Houches common block structure for decay of \n" << *external[0] << "\n" << hepeup << Exception::runerror; } if(unMatchedAntiColour.size() == 3 ) { unMatchedAntiColour[0]->setSinkNeighbours(unMatchedAntiColour[1], unMatchedAntiColour[2]); } else if(unMatchedAntiColour.size() != 0 && ThePEG_DEBUG_LEVEL) { Throw() << "FxFxReader '" << name() << "' found inconsistent colour " << "flow in Les Houches common block structure for decay of\n" << *external[0] << "\n" << hepeup << Exception::runerror; } } } } void FxFxReader::createBeams() { if ( !theBeams.first && dynamic_ptr_cast::tcp>(inPDF.first) ) { theBeams.first = theIncoming.first; } else if ( !theBeams.first ) { theBeams.first = getParticleData(heprup.IDBMUP.first)->produceParticle(); double m = theBeams.first->mass()/GeV; theBeams.first->set5Momentum (Lorentz5Momentum(ZERO, ZERO, sqrt(sqr(heprup.EBMUP.first) - sqr(m))*GeV, heprup.EBMUP.first*GeV, m*GeV)); hepeup.IDUP.push_back(heprup.IDBMUP.first); hepeup.ISTUP.push_back(-9); hepeup.MOTHUP.push_back(make_pair(0, 0)); hepeup.ICOLUP.push_back(make_pair(0, 0)); hepeup.VTIMUP.push_back(0.0); hepeup.SPINUP.push_back(0.0); particleIndex(hepeup.IDUP.size(), theBeams.first); hepeup.MOTHUP[particleIndex(theIncoming.first) - 1].first = hepeup.IDUP.size(); } if ( !theBeams.second && dynamic_ptr_cast::tcp>(inPDF.second) ) { theBeams.second = theIncoming.second; } else if ( !theBeams.second ) { theBeams.second = getParticleData(heprup.IDBMUP.second)->produceParticle(); double m = theBeams.second->mass()/GeV; theBeams.second->set5Momentum (Lorentz5Momentum(ZERO, ZERO, -sqrt(sqr(heprup.EBMUP.second) - sqr(m))*GeV, heprup.EBMUP.second*GeV, m*GeV)); hepeup.IDUP.push_back(heprup.IDBMUP.second); hepeup.ISTUP.push_back(-9); hepeup.MOTHUP.push_back(make_pair(0, 0)); hepeup.ICOLUP.push_back(make_pair(0, 0)); hepeup.VTIMUP.push_back(0.0); hepeup.SPINUP.push_back(0.0); particleIndex(hepeup.IDUP.size(), theBeams.second); hepeup.MOTHUP[particleIndex(theIncoming.second) - 1].first = hepeup.IDUP.size(); } } void FxFxReader::connectMothers() { const ObjectIndexer & pi = particleIndex; for ( int i = 0, N = hepeup.IDUP.size(); i < N; ++i ) { if ( pi(hepeup.MOTHUP[i].first) ) pi(hepeup.MOTHUP[i].first)->addChild(pi(i + 1)); if ( pi(hepeup.MOTHUP[i].second) && hepeup.MOTHUP[i].second != hepeup.MOTHUP[i].first ) pi(hepeup.MOTHUP[i].second)->addChild(pi(i + 1)); } } void FxFxReader::openReadCacheFile() { if ( cacheFile() ) closeCacheFile(); cacheFile().open(cacheFileName(), "r"); position = 0; } void FxFxReader::openWriteCacheFile() { if ( cacheFile() ) closeCacheFile(); cacheFile().open(cacheFileName(), "w"); } void FxFxReader::closeCacheFile() { cacheFile().close(); } void FxFxReader::cacheEvent() const { static vector buff; cacheFile().write(&hepeup.NUP, sizeof(hepeup.NUP)); buff.resize(eventSize(hepeup.NUP)); char * pos = &buff[0]; pos = mwrite(pos, hepeup.IDPRUP); pos = mwrite(pos, hepeup.XWGTUP); pos = mwrite(pos, hepeup.XPDWUP); pos = mwrite(pos, hepeup.SCALUP); pos = mwrite(pos, hepeup.AQEDUP); pos = mwrite(pos, hepeup.AQCDUP); pos = mwrite(pos, hepeup.IDUP[0], hepeup.NUP); pos = mwrite(pos, hepeup.ISTUP[0], hepeup.NUP); pos = mwrite(pos, hepeup.MOTHUP[0], hepeup.NUP); pos = mwrite(pos, hepeup.ICOLUP[0], hepeup.NUP); for ( int i = 0; i < hepeup.NUP; ++i ) pos = mwrite(pos, hepeup.PUP[i][0], 5); pos = mwrite(pos, hepeup.VTIMUP[0], hepeup.NUP); pos = mwrite(pos, hepeup.SPINUP[0], hepeup.NUP); pos = mwrite(pos, lastweight); pos = mwrite(pos, optionalWeights); + pos = mwrite(pos, LHEeventnum); for(size_t ff = 0; ff < optionalWeightsNames.size(); ff++) { pos = mwrite(pos, optionalWeightsNames[ff]); } /* for(int f = 0; f < optionalWeightsNames.size(); f++) { cout << "optionalWeightsNames = " << optionalWeightsNames[f] << endl; }*/ pos = mwrite(pos, optionalnpLO); pos = mwrite(pos, optionalnpNLO); pos = mwrite(pos, preweight); cacheFile().write(&buff[0], buff.size(), 1); } bool FxFxReader::uncacheEvent() { reset(); static vector buff; if ( cacheFile().read(&hepeup.NUP, sizeof(hepeup.NUP)) != 1 ) return false; buff.resize(eventSize(hepeup.NUP)); if ( cacheFile().read(&buff[0], buff.size()) != 1 ) return false; const char * pos = &buff[0]; pos = mread(pos, hepeup.IDPRUP); pos = mread(pos, hepeup.XWGTUP); pos = mread(pos, hepeup.XPDWUP); pos = mread(pos, hepeup.SCALUP); pos = mread(pos, hepeup.AQEDUP); pos = mread(pos, hepeup.AQCDUP); hepeup.IDUP.resize(hepeup.NUP); pos = mread(pos, hepeup.IDUP[0], hepeup.NUP); hepeup.ISTUP.resize(hepeup.NUP); pos = mread(pos, hepeup.ISTUP[0], hepeup.NUP); hepeup.MOTHUP.resize(hepeup.NUP); pos = mread(pos, hepeup.MOTHUP[0], hepeup.NUP); hepeup.ICOLUP.resize(hepeup.NUP); pos = mread(pos, hepeup.ICOLUP[0], hepeup.NUP); hepeup.PUP.resize(hepeup.NUP); for ( int i = 0; i < hepeup.NUP; ++i ) pos = mread(pos, hepeup.PUP[i][0], 5); hepeup.VTIMUP.resize(hepeup.NUP); pos = mread(pos, hepeup.VTIMUP[0], hepeup.NUP); hepeup.SPINUP.resize(hepeup.NUP); pos = mread(pos, hepeup.SPINUP[0], hepeup.NUP); pos = mread(pos, lastweight); pos = mread(pos, optionalWeights); for(size_t ff = 0; ff < optionalWeightsNames.size(); ff++) { pos = mread(pos, optionalWeightsNames[ff]); } pos = mread(pos, optionalnpLO); pos = mread(pos, optionalnpNLO); pos = mread(pos, preweight); + pos = mread(pos, LHEeventnum); // If we are skipping, we do not have to do anything else. if ( skipping ) return true; if ( CKKWHandler() && maxMultCKKW() > 0 && maxMultCKKW() > minMultCKKW() ) { // The cached event has not been submitted to CKKW reweighting, so // we do that now. fillEvent(); getSubProcess(); CKKWHandler()->setXComb(lastXCombPtr()); lastweight *= CKKWHandler()->reweightCKKW(minMultCKKW(), maxMultCKKW()); } return true; } void FxFxReader::persistentOutput(PersistentOStream & os) const { os << heprup.IDBMUP << heprup.EBMUP << heprup.PDFGUP << heprup.PDFSUP << heprup.IDWTUP << heprup.NPRUP << heprup.XSECUP << heprup.XERRUP << heprup.XMAXUP << heprup.LPRUP << hepeup.NUP << hepeup.IDPRUP << hepeup.XWGTUP << hepeup.XPDWUP << hepeup.SCALUP << hepeup.AQEDUP << hepeup.AQCDUP << hepeup.IDUP << hepeup.ISTUP << hepeup.MOTHUP << hepeup.ICOLUP << hepeup.PUP << hepeup.VTIMUP << hepeup.SPINUP << inData << inPDF << outPDF << thePartonExtractor << theCKKW << thePartonBins << theXCombs << theCuts << theNEvents << position << reopened << theMaxScan << isActive << theCacheFileName << doCutEarly << stats << statmap << thePartonBinInstances << theBeams << theIncoming << theOutgoing << theIntermediates << reweights << preweights << preweight << reweightPDF << doInitPDFs - << theLastXComb << theMaxMultCKKW << theMinMultCKKW << lastweight << optionalWeights << optionalnpLO << optionalnpNLO + << theLastXComb << theMaxMultCKKW << theMinMultCKKW << lastweight << optionalWeights << optionalnpLO << optionalnpNLO << LHEeventnum << maxFactor << ounit(weightScale, picobarn) << xSecWeights << maxWeights << theMomentumTreatment << useWeightWarnings << theReOpenAllowed << theIncludeSpin; } void FxFxReader::persistentInput(PersistentIStream & is, int) { if ( cacheFile() ) closeCacheFile(); is >> heprup.IDBMUP >> heprup.EBMUP >> heprup.PDFGUP >> heprup.PDFSUP >> heprup.IDWTUP >> heprup.NPRUP >> heprup.XSECUP >> heprup.XERRUP >> heprup.XMAXUP >> heprup.LPRUP >> hepeup.NUP >> hepeup.IDPRUP >> hepeup.XWGTUP >> hepeup.XPDWUP >> hepeup.SCALUP >> hepeup.AQEDUP >> hepeup.AQCDUP >> hepeup.IDUP >> hepeup.ISTUP >> hepeup.MOTHUP >> hepeup.ICOLUP >> hepeup.PUP >> hepeup.VTIMUP >> hepeup.SPINUP >> inData >> inPDF >> outPDF >> thePartonExtractor >> theCKKW >> thePartonBins >> theXCombs >> theCuts >> theNEvents >> position >> reopened >> theMaxScan >> isActive >> theCacheFileName >> doCutEarly >> stats >> statmap >> thePartonBinInstances >> theBeams >> theIncoming >> theOutgoing >> theIntermediates >> reweights >> preweights >> preweight >> reweightPDF >> doInitPDFs - >> theLastXComb >> theMaxMultCKKW >> theMinMultCKKW >> lastweight >> optionalWeights >> optionalnpLO >> optionalnpNLO + >> theLastXComb >> theMaxMultCKKW >> theMinMultCKKW >> lastweight >> optionalWeights >> optionalnpLO >> optionalnpNLO >> LHEeventnum >> maxFactor >> iunit(weightScale, picobarn) >> xSecWeights >> maxWeights >> theMomentumTreatment >> useWeightWarnings >> theReOpenAllowed >> theIncludeSpin; } AbstractClassDescription FxFxReader::initFxFxReader; // Definition of the static class description member. void FxFxReader::setBeamA(long id) { heprup.IDBMUP.first = id; } long FxFxReader::getBeamA() const { return heprup.IDBMUP.first; } void FxFxReader::setBeamB(long id) { heprup.IDBMUP.second = id; } long FxFxReader::getBeamB() const { return heprup.IDBMUP.second; } void FxFxReader::setEBeamA(Energy e) { heprup.EBMUP.first = e/GeV; } Energy FxFxReader::getEBeamA() const { return heprup.EBMUP.first*GeV; } void FxFxReader::setEBeamB(Energy e) { heprup.EBMUP.second = e/GeV; } Energy FxFxReader::getEBeamB() const { return heprup.EBMUP.second*GeV; } void FxFxReader::setPDFA(PDFPtr pdf) { inPDF.first = pdf; } PDFPtr FxFxReader::getPDFA() const { return inPDF.first; } void FxFxReader::setPDFB(PDFPtr pdf) { inPDF.second = pdf; } PDFPtr FxFxReader::getPDFB() const { return inPDF.second; } void FxFxReader::Init() { static ClassDocumentation documentation ("ThePEG::FxFxReader is an abstract base class to be used " "for objects which reads event files or streams from matrix element " "generators."); static Parameter interfaceBeamA ("BeamA", "The PDG id of the incoming particle along the positive z-axis. " "If zero the corresponding information is to be deduced from the " "event stream/file.", 0, 0, 0, 0, true, false, false, &FxFxReader::setBeamA, &FxFxReader::getBeamA, 0, 0, 0); static Parameter interfaceBeamB ("BeamB", "The PDG id of the incoming particle along the negative z-axis. " "If zero the corresponding information is to be deduced from the " "event stream/file.", 0, 0, 0, 0, true, false, false, &FxFxReader::setBeamB, &FxFxReader::getBeamB, 0, 0, 0); static Parameter interfaceEBeamA ("EBeamA", "The energy of the incoming particle along the positive z-axis. " "If zero the corresponding information is to be deduced from the " "event stream/file.", 0, GeV, ZERO, ZERO, 1000000000.0*GeV, true, false, true, &FxFxReader::setEBeamA, &FxFxReader::getEBeamA, 0, 0, 0); static Parameter interfaceEBeamB ("EBeamB", "The energy of the incoming particle along the negative z-axis. " "If zero the corresponding information is to be deduced from the " "event stream/file.", 0, GeV, ZERO, ZERO, 1000000000.0*GeV, true, false, true, &FxFxReader::setEBeamB, &FxFxReader::getEBeamB, 0, 0, 0); static Reference interfacePDFA ("PDFA", "The PDF used for incoming particle along the positive z-axis. " "If null the corresponding information is to be deduced from the " "event stream/file.", 0, true, false, true, true, false, &FxFxReader::setPDFA, &FxFxReader::getPDFA, 0); static Reference interfacePDFB ("PDFB", "The PDF used for incoming particle along the negative z-axis. " "If null the corresponding information is to be deduced from the " "event stream/file.", 0, true, false, true, true, false, &FxFxReader::setPDFB, &FxFxReader::getPDFB, 0); static Parameter interfaceMaxScan ("MaxScan", "The maximum number of events to scan to obtain information about " "processes and cross section in the intialization.", &FxFxReader::theMaxScan, -1, 0, 0, true, false, false); static Parameter interfaceCacheFileName ("CacheFileName", "Name of file used to cache the events form the reader in a fast-readable " "form. If empty, no cache file will be generated.", &FxFxReader::theCacheFileName, "", true, false); interfaceCacheFileName.fileType(); static Switch interfaceCutEarly ("CutEarly", "Determines whether to apply cuts to events before converting to " "ThePEG format.", &FxFxReader::doCutEarly, true, true, false); static SwitchOption interfaceCutEarlyYes (interfaceCutEarly, "Yes", "Event are cut before converted.", true); static SwitchOption interfaceCutEarlyNo (interfaceCutEarly, "No", "Events are not cut before converted.", false); static Reference interfacePartonExtractor ("PartonExtractor", "The PartonExtractor object used to construct remnants. If no object is " "provided the FxFxEventHandler object must provide one instead.", &FxFxReader::thePartonExtractor, true, false, true, true, false); static Reference interfaceCuts ("Cuts", "The Cuts object to be used for this reader. Note that these " "must not be looser cuts than those used in the actual generation. " "If no object is provided the FxFxEventHandler object must " "provide one instead.", &FxFxReader::theCuts, true, false, true, true, false); static RefVector interfaceReweights ("Reweights", "A list of ThePEG::ReweightBase objects to modify this the weight of " "this reader.", &FxFxReader::reweights, 0, false, false, true, false); static RefVector interfacePreweights ("Preweights", "A list of ThePEG::ReweightBase objects to bias the phase space for this " "reader without influencing the actual cross section.", &FxFxReader::preweights, 0, false, false, true, false); static Switch interfaceReweightPDF ("ReweightPDF", "If the PDFs used in the generation for this reader is different " "from the ones assumed by the associated PartonExtractor object, " "should the events be reweighted to fit the latter?", &FxFxReader::reweightPDF, false, true, false); static SwitchOption interfaceReweightPDFNo (interfaceReweightPDF, "No", "The event weights are kept as they are.", false); static SwitchOption interfaceReweightPDFYes (interfaceReweightPDF, "Yes", "The events are reweighted.", true); static Switch interfaceInitPDFs ("InitPDFs", "If no PDFs were specified in PDFA or " "PDFBfor this reader, try to extract the " "information from the event file and assign the relevant PDFBase" "objects when the reader is initialized.", &FxFxReader::doInitPDFs, false, true, false); static SwitchOption interfaceInitPDFsYes (interfaceInitPDFs, "Yes", "Extract PDFs during initialization.", true); static SwitchOption interfaceInitPDFsNo (interfaceInitPDFs, "No", "Do not extract PDFs during initialization.", false); static Parameter interfaceMaxMultCKKW ("MaxMultCKKW", "If this reader is to be used (possibly together with others) for CKKW-" "reweighting and veto, this should give the multiplicity of outgoing " "particles in the highest multiplicity matrix element in the group. " "If set to zero, no CKKW procedure should be applied.", &FxFxReader::theMaxMultCKKW, 0, 0, 0, true, false, Interface::lowerlim); static Parameter interfaceMinMultCKKW ("MinMultCKKW", "If this reader is to be used (possibly together with others) for CKKW-" "reweighting and veto, this should give the multiplicity of outgoing " "particles in the lowest multiplicity matrix element in the group. If " "larger or equal to MaxMultCKKW, no CKKW " "procedure should be applied.", &FxFxReader::theMinMultCKKW, 0, 0, 0, true, false, Interface::lowerlim); static Switch interfaceMomentumTreatment ("MomentumTreatment", "Treatment of the momenta supplied by the interface", &FxFxReader::theMomentumTreatment, 0, false, false); static SwitchOption interfaceMomentumTreatmentAccept (interfaceMomentumTreatment, "Accept", "Just accept the momenta given", 0); static SwitchOption interfaceMomentumTreatmentRescaleEnergy (interfaceMomentumTreatment, "RescaleEnergy", "Rescale the energy supplied so it is consistent with the mass", 1); static SwitchOption interfaceMomentumTreatmentRescaleMass (interfaceMomentumTreatment, "RescaleMass", "Rescale the mass supplied so it is consistent with the" " energy and momentum", 2); static Switch interfaceWeightWarnings ("WeightWarnings", "Determines if warnings about possible weight incompatibilities should " "be issued when this reader is initialized.", &FxFxReader::useWeightWarnings, true, true, false); static SwitchOption interfaceWeightWarningsWarnAboutWeights (interfaceWeightWarnings, "WarnAboutWeights", "Warn about possible incompatibilities with the weight option in the " "Les Houches common block and the requested weight treatment.", true); static SwitchOption interfaceWeightWarningsDontWarnAboutWeights (interfaceWeightWarnings, "DontWarnAboutWeights", "Do not warn about possible incompatibilities with the weight option " "in the Les Houches common block and the requested weight treatment.", false); static Switch interfaceAllowedTopReOpen ("AllowedToReOpen", "Can the file be reopened if more events are requested than the file contains?", &FxFxReader::theReOpenAllowed, true, false, false); static SwitchOption interfaceAllowedTopReOpenYes (interfaceAllowedTopReOpen, "Yes", "Allowed to reopen the file", true); static SwitchOption interfaceAllowedTopReOpenNo (interfaceAllowedTopReOpen, "No", "Not allowed to reopen the file", false); static Switch interfaceIncludeSpin ("IncludeSpin", "Use the spin information present in the event file, for tau leptons" " only as this is the only case which makes any sense", &FxFxReader::theIncludeSpin, true, false, false); static SwitchOption interfaceIncludeSpinYes (interfaceIncludeSpin, "Yes", "Use the spin information", true); static SwitchOption interfaceIncludeSpinNo (interfaceIncludeSpin, "No", "Don't use the spin information", false); interfaceCuts.rank(8); interfacePartonExtractor.rank(7); interfaceBeamA.rank(5); interfaceBeamB.rank(4); interfaceEBeamA.rank(3); interfaceEBeamB.rank(2); interfaceMaxMultCKKW.rank(1.5); interfaceMinMultCKKW.rank(1.0); interfaceBeamA.setHasDefault(false); interfaceBeamB.setHasDefault(false); interfaceEBeamA.setHasDefault(false); interfaceEBeamB.setHasDefault(false); interfaceMaxMultCKKW.setHasDefault(false); interfaceMinMultCKKW.setHasDefault(false); } namespace ThePEG { ostream & operator<<(ostream & os, const HEPEUP & h) { os << "\n" << " " << setw(4) << h.NUP << " " << setw(6) << h.IDPRUP << " " << setw(14) << h.XWGTUP << " " << setw(14) << h.SCALUP << " " << setw(14) << h.AQEDUP << " " << setw(14) << h.AQCDUP << "\n"; for ( int i = 0; i < h.NUP; ++i ) os << " " << setw(8) << h.IDUP[i] << " " << setw(2) << h.ISTUP[i] << " " << setw(4) << h.MOTHUP[i].first << " " << setw(4) << h.MOTHUP[i].second << " " << setw(4) << h.ICOLUP[i].first << " " << setw(4) << h.ICOLUP[i].second << " " << setw(14) << h.PUP[i][0] << " " << setw(14) << h.PUP[i][1] << " " << setw(14) << h.PUP[i][2] << " " << setw(14) << h.PUP[i][3] << " " << setw(14) << h.PUP[i][4] << " " << setw(1) << h.VTIMUP[i] << " " << setw(1) << h.SPINUP[i] << std::endl; os << "" << std::endl; return os; } } diff --git a/MatrixElement/FxFx/FxFxReader.h b/MatrixElement/FxFx/FxFxReader.h --- a/MatrixElement/FxFx/FxFxReader.h +++ b/MatrixElement/FxFx/FxFxReader.h @@ -1,1006 +1,1016 @@ // -*- C++ -*- // // FxFxReader.h is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2019 Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef THEPEG_FxFxReader_H #define THEPEG_FxFxReader_H // This is the declaration of the FxFxReader class. #include "FxFx.h" #include "ThePEG/Handlers/HandlerBase.h" #include "ThePEG/Utilities/ObjectIndexer.h" #include "ThePEG/Utilities/Exception.h" #include "ThePEG/Utilities/XSecStat.h" #include "ThePEG/PDF/PartonBinInstance.h" #include "ThePEG/PDF/PartonBin.fh" #include "ThePEG/MatrixElement/ReweightBase.h" #include "FxFxEventHandler.fh" #include "FxFxReader.fh" #include "ThePEG/Utilities/CFile.h" #include #include namespace ThePEG { /** * FxFxReader is an abstract base class to be used for objects * which reads event files or streams from matrix element * generators. Derived classes must at least implement the open() and * doReadEvent() methods to read in information about the whole run into * the HEPRUP variable and next event into the HEPEUP variable * respectively. Also the close() function to close the file or stream * read must be implemented. Although these functions are named as if * we are reading from event files, they could just as well implement * the actual generation of events. * * After filling the HEPRUP and HEPEUP variables, which are protected * and easily accesible from the sub-class, this base class will then * be responsible for transforming this data to the ThePEG Event * record in the getEvent() method. FxFxReaders can * only be used inside FxFxEventHandler objects. * * In the initialization the virtual open() and scan() functions are * called. Here the derived class must provide the information about * the processes in the variables corresponding to the HEPRUP common * block. Note that the IDWTUP is required to be +/- 1, and sub * classes are required to change the information accordingly to * ensure the correct corss section sampling. Note also that the * controlling FxFxEventHandler may choose to generate weighted * events even if IDWTUP is 1. * * Note that the information given per process in e.g. the XSECUP and * XMAXUP vectors is not used by the FxFxEventHandler and by * default the FxFxReader is not assumed to be able to actively * choose between the sub-processes. Instead, the * FxFxEventHandler can handle several FxFxReader objects * and choose between them. However, a sub-class of FxFxReader * may set the flag isActive, in which case it is assumed to be able * to select between its sub-processes itself. * * The FxFxReader may be assigned a number ReweightBase objects * which either completely reweights the events produced (in the * reweights vector), or only biases the selection without influencing * the cross section (in the preweights vector). Note that it is the * responsibility of a sub-class to call the reweight() function and * multiply the weight according to its return value (typically done * in the readEvent() function). * * @see \ref FxFxReaderInterfaces "The interfaces" * defined for FxFxReader. * @see Event * @see FxFxEventHandler */ class FxFxReader: public HandlerBase, public LastXCombInfo<> { /** * FxFxEventHandler should have access to our private parts. */ friend class FxFxEventHandler; /** * Map for accumulating statistics of cross sections per process * number. */ typedef map StatMap; /** * Map of XComb objects describing the incoming partons indexed by * the corresponding PartonBin pair. */ typedef map XCombMap; /** * A vector of pointers to ReweightBase objects. */ typedef vector ReweightVector; public: /** @name Standard constructors and destructors. */ //@{ /** * Default constructor. If the optional argument is true, the reader * is assumed to be able to produce events on demand for a given * process. */ FxFxReader(bool active = false); /** * Copy-constructor. */ FxFxReader(const FxFxReader &); /** * Destructor. */ virtual ~FxFxReader(); //@} public: /** @name Main virtual fuctions to be overridden in * sub-classes. They are named as if we are reading from event * files, but could equally well implement the actual generation of * events. */ //@{ /** * Open a file or stream with events and read in the run information * into the heprup variable. */ virtual void open() = 0; /** * Read the next event from the file or stream into the * corresponding protected variables. Return false if there is no * more events. */ virtual bool doReadEvent() = 0; /** * Close the file or stream from which events have been read. */ virtual void close() = 0; /** * return the weight names */ // virtual vector optWeightsNamesFunc(); virtual vector optWeightsNamesFunc() = 0; //virtual vector optWeightNamesFunc() = 0; vector optionalWeightsNames; /** * The ID (e.g. 100x, 2001) for the weight */ // vector optionalWeightsNames; //@} /** @name Other important function which may be overridden in * sub-classes which wants to bypass the basic HEPRUP or HEPEUP * variables or otherwise facilitate the conversion to ThePEG * objects. */ //@{ /** * Initialize. This function is called by the FxFxEventHandler * to which this object is assigned. */ virtual void initialize(FxFxEventHandler & eh); /** * Calls readEvent() or uncacheEvent() to read information into the * FxFx common block variables. This function is called by the * FxFxEventHandler if this reader has been selectod to * produce an event. * * @return the weight asociated with this event. If negative weights * are allowed it should be between -1 and 1, otherwise between 0 * and 1. If outside these limits the previously estimated maximum * is violated. Note that the estimated maximum then should be * updated from the outside. */ virtual double getEvent(); /** * Calls doReadEvent() and performs pre-defined reweightings. A * sub-class overrides this function it must make sure that the * corresponding reweightings are done. */ virtual bool readEvent(); /** * Skip \a n events. Used by FxFxEventHandler to make sure * that a file is scanned an even number of times in case the events * are not ramdomly distributed in the file. */ virtual void skip(long n); /** * Get an XComb object. Converts the information in the Les Houches * common block variables to an XComb object describing the sub * process. This is the way information is conveyed from the reader * to the controlling FxFxEventHandler. */ tXCombPtr getXComb(); /** * Get a SubProcess object corresponding to the information in the * Les Houches common block variables. */ tSubProPtr getSubProcess(); /** * Scan the file or stream to obtain information about cross section * weights and particles etc. This function should fill the * variables corresponding to the /HEPRUP/ common block. The * function returns the number of events scanned. */ virtual long scan(); /** * Take the information corresponding to the HEPRUP common block and * initialize the statistics for this reader. */ virtual void initStat(); /** * Reweights the current event using the reweights and preweights * vectors. It is the responsibility of the sub-class to call this * function after the HEPEUP information has been retrieved. */ double reweight(); /** * Converts the information in the Les Houches common block * variables into a Particle objects. */ virtual void fillEvent(); /** * Removes the particles created in the last generated event, * preparing to produce a new one. */ void reset(); /** * Possibility for subclasses to recover from non-conformant * settings of XMAXUP when an event file has been scanned with \a * neve events. Should set weightScale so that the average XMAXUP * times weightScale gives the cross section for a process. (This is * needed for MadEvent). */ virtual void setWeightScale(long neve); //@} /** @name Access information about the current event. */ //@{ /** * Return the size of this event in bytes. To be used for the cache * file. \a npart is the number of particles. If \a npart is 0, the * number is taken from NUP. */ static size_t eventSize(int N) { return (N + 1)*sizeof(int) + // IDPRUP, ISTUP (7*N + 4)*sizeof(double) + // XWGTUP, SCALUP, AQEDUP, AQCDUP, PUP, // VTIMUP, SPINUP N*sizeof(long) + // IDUP 2*N*sizeof(pair) + // MOTHUP, ICOLUP sizeof(pair) + // XPDWUP. 2*sizeof(double); // lastweight and preweight } /** * The current event weight given by XWGTUP times possible * reweighting. Note that this is not necessarily the same as what * is returned by getEvent(), which is scaled with the maximum * weight. */ double eventWeight() const { return hepeup.XWGTUP*lastweight; } /** * Return the optional named weights associated to the current event. */ const map& optionalEventWeights() const { return optionalWeights; } + /** + * Return the Les Houches event number associated with the current event + */ + const long& LHEEventNum() const { return LHEeventnum; } + /** * Return the optional npLO and npNLO */ const int& optionalEventnpLO() const { return optionalnpLO; } const int& optionalEventnpNLO() const { return optionalnpNLO; } /** * The pair of PartonBinInstance objects describing the current * incoming partons in the event. */ const PBIPair & partonBinInstances() const { return thePartonBinInstances; } /** * Return the instances of the beam particles for the current event. */ const PPair & beams() const { return theBeams; } /** * Return the instances of the incoming particles to the sub process * for the current event. */ const PPair & incoming() const { return theIncoming; } /** * Return the instances of the outgoing particles from the sub process * for the current event. */ const PVector & outgoing() const { return theOutgoing; } /** * Return the instances of the intermediate particles in the sub * process for the current event. */ const PVector & intermediates() const { return theIntermediates; } /** * If this reader is to be used (possibly together with others) for * CKKW reweighting and veto, this should give the multiplicity of * outgoing particles in the highest multiplicity matrix element in * the group. */ int maxMultCKKW() const { return theMaxMultCKKW; } /** * If this reader is to be used (possibly together with others) for * CKKW reweighting and veto, this should give the multiplicity of * outgoing particles in the lowest multiplicity matrix element in * the group. */ int minMultCKKW() const { return theMinMultCKKW; } //@} /** @name Other inlined access functions. */ //@{ /** * The number of events found in this reader. If less than zero the * number of events are unlimited. */ long NEvents() const { return theNEvents; } /** * The number of events produced so far. Is reset to zero if an * event file is reopened. */ long currentPosition() const { return position; } /** * The maximum number of events to scan to collect information about * processes and cross sections. If less than 0, all events will be * scanned. */ long maxScan() const { return theMaxScan; } /** * Return true if this reader is active. */ bool active() const { return isActive; } /** * True if negative weights may be produced. */ bool negativeWeights() const { return heprup.IDWTUP < 0; } /** * The collected cross section statistics for this reader. */ const XSecStat & xSecStats() const { return stats; } /** * Collected statistics about the individual processes. */ const StatMap & processStats() const { return statmap; } /** * Select the current event. It will later be rejected with a * probability given by \a weight. */ void select(double weight) { stats.select(weight); statmap[hepeup.IDPRUP].select(weight); } /** * Accept the current event assuming it was previously selcted. */ void accept() { stats.accept(); statmap[hepeup.IDPRUP].accept(); } /** * Reject the current event assuming it was previously accepted. */ void reject(double w) { stats.reject(w); statmap[hepeup.IDPRUP].reject(w); } /** * Increase the overestimated cross section for this reader. */ virtual void increaseMaxXSec(CrossSection maxxsec); /** * The PartonExtractor object used to construct remnants. */ tPExtrPtr partonExtractor() const { return thePartonExtractor; } /** * Return a possibly null pointer to a CascadeHandler to be used for * CKKW-reweighting. */ tCascHdlPtr CKKWHandler() const { return theCKKW; } /** * The pairs of PartonBin objects describing the partons which can * be extracted by the PartonExtractor object. */ const PartonPairVec & partonBins() const { return thePartonBins; } /** * The map of XComb objects indexed by the corresponding PartonBin * pair. */ const XCombMap & xCombs() const { return theXCombs; } /** * The Cuts object to be used for this reader. */ const Cuts & cuts() const { return *theCuts; } //@} protected: /** @name Functions for manipulating cache files. */ //@{ /** * Name of file used to cache the events form the reader in a * fast-readable form. If empty, no cache file will be generated. */ string cacheFileName() const { return theCacheFileName; } /** * Determines whether to apply cuts to events converting them to * ThePEG format. */ bool cutEarly() const { return doCutEarly; } /** * File stream for the cache. */ CFile cacheFile() const { return theCacheFile;} /** * Open the cache file for reading. */ void openReadCacheFile(); /** * Open the cache file for writing. */ void openWriteCacheFile(); /** * Close the cache file; */ void closeCacheFile(); /** * Write the current event to the cache file. */ void cacheEvent() const; /** * Read an event from the cache file. Return false if something went wrong. */ bool uncacheEvent(); /** * Reopen a reader. If we have reached the end of an event file, * reopen it and issue a warning if we have used up a large fraction * of it. */ void reopen(); /** * Helper function to write a variable to a memory location */ template static char * mwrite(char * pos, const T & t, size_t n = 1) { std::memcpy(pos, &t, n*sizeof(T)); return pos + n*sizeof(T); } /** * Helper function to read a variable from a memory location */ template static const char * mread(const char * pos, T & t, size_t n = 1) { std::memcpy(&t, pos, n*sizeof(T)); return pos + n*sizeof(T); } //@} /** @name Auxilliary virtual methods which may be verridden by sub-classes. */ //@{ /** * Check the existence of a pair of PartonBin objects corresponding * to the current event. * * @return false if no pair of suitable PartonBin objects was found. */ virtual bool checkPartonBin(); /** * Create instances of all particles in the event and store them * in particleIndex. */ virtual void createParticles(); /** * Using the already created particles create a pair of * PartonBinInstance objects corresponding to the incoming * partons. Return the corresponding PartonBin objects. */ virtual tcPBPair createPartonBinInstances(); /** * Create instances of the incoming beams in the event and store * them in particleIndex. If no beam particles are included in the * event they are created from the run info. */ virtual void createBeams(); /** * Go through the mother indices and connect up the Particles. */ virtual void connectMothers(); //@} public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * Standard Init function used to initialize the interfaces. */ static void Init(); protected: /** @name Set functions for some variables not in the Les Houches accord. */ //@{ /** * The number of events in this reader. If less than zero the number * of events is unlimited. */ void NEvents(long x) { theNEvents = x; } /** * The map of XComb objects indexed by the corresponding PartonBin * pair. */ XCombMap & xCombs() { return theXCombs; } //@} /** @name Standard (and non-standard) Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); /** * Initialize this object. Called in the run phase just before * a run begins. */ virtual void doinitrun(); /** * Finalize this object. Called in the run phase just after a * run has ended. Used eg. to write out statistics. */ virtual void dofinish() { close(); HandlerBase::dofinish(); } /** * Return true if this object needs to be initialized before all * other objects because it needs to extract PDFs from the event file. */ virtual bool preInitialize() const; /** * Called from doinit() to extract PDFs from the event file and add * the corresponding objects to the current EventGenerator. */ virtual void initPDFs(); //@} protected: /** * The HEPRUP common block. */ HEPRUP heprup; /** * The HEPEUP common block. */ HEPEUP hepeup; /** * The ParticleData objects corresponding to the incoming particles. */ tcPDPair inData; /** * The PDFBase objects which has been used for the beam particle * when generating the events being read. Specified in the interface * or derived from PDFGUP and PDFSUP. */ pair inPDF; /** * The PDFBase object to be used in the subsequent generation. */ pair outPDF; /** * The PartonExtractor object used to construct remnants. */ PExtrPtr thePartonExtractor; /** * A pointer to a CascadeHandler to be used for CKKW-reweighting. */ tCascHdlPtr theCKKW; /** * The pairs of PartonBin objects describing the partons which can * be extracted by the PartonExtractor object. */ PartonPairVec thePartonBins; /** * The map of XComb objects indexed by the corresponding PartonBin * pair. */ XCombMap theXCombs; /** * The Cuts object to be used for this reader. */ CutsPtr theCuts; /** * The number of events in this reader. If less than zero the number * of events is unlimited. */ long theNEvents; /** * The number of events produced by this reader so far. Is reset * every time an event file is reopened. */ long position; /** * The number of times this reader has been reopened. */ int reopened; /** * The maximum number of events to scan to collect information about * processes and cross sections. If less than 0, all events will be * scanned. */ long theMaxScan; /** * Flag to tell whether we are in the process of scanning. */ bool scanning; /** * True if this is an active reader. */ bool isActive; /** * Name of file used to cache the events form the reader in a * fast-readable form. If empty, no cache file will be generated. */ string theCacheFileName; /** * Determines whether to apply cuts to events before converting them * to ThePEG format. */ bool doCutEarly; /** * Collect statistics for this reader. */ XSecStat stats; /** * Collect statistics for each individual process. */ StatMap statmap; /** * The pair of PartonBinInstance objects describing the current * incoming partons in the event. */ PBIPair thePartonBinInstances; /** * Association between ColourLines and colour indices in the current * translation. */ ObjectIndexer colourIndex; /** * Association between Particles and indices in the current * translation. */ ObjectIndexer particleIndex; /** * The instances of the beam particles for the current event. */ PPair theBeams; /** * The instances of the incoming particles to the sub process for * the current event. */ PPair theIncoming; /** * The instances of the outgoing particles from the sub process for * the current event. */ PVector theOutgoing; /** * The instances of the intermediate particles in the sub process for * the current event. */ PVector theIntermediates; /** * File stream for the cache. */ CFile theCacheFile; /** * The reweight objects modifying the weights of this reader. */ ReweightVector reweights; /** * The preweight objects modifying the weights of this reader. */ ReweightVector preweights; /** * The factor with which this reader was last pre-weighted. */ double preweight; /** * Should the event be reweighted by PDFs used by the PartonExtractor? */ bool reweightPDF; /** * Should PDFBase objects be constructed from the information in the * event file in the initialization? */ bool doInitPDFs; /** * If this reader is to be used (possibly together with others) for * CKKW reweighting and veto, this should give the multiplicity of * outgoing particles in the highest multiplicity matrix element in * the group. */ int theMaxMultCKKW; /** * If this reader is to be used (possibly together with others) for * CKKW reweighting and veto, this should give the multiplicity of * outgoing particles in the lowest multiplicity matrix element in * the group. */ int theMinMultCKKW; /** * The weight multiplying the last read event due to PDF * reweighting, CKKW reweighting or assigned reweight and preweight * objects. */ double lastweight; /** * The optional weights associated to the last read events. */ map optionalWeights; + + /** + * The event number + */ + long LHEeventnum; /** * If the maximum cross section of this reader has been increased * with increaseMaxXSec(), this is the total factor with which it * has been increased. */ double maxFactor; /** * npLO for FxFx merging */ int optionalnpLO; /** * npNLO for FxFx merging */ int optionalnpNLO; /** * The (reweighted) XWGTUP value should be scaled with this cross * section when compared to the overestimated cross section. */ CrossSection weightScale; /** * Individual scales for different sub-processes if reweighted. */ vector xSecWeights; /** * Individual maximum weights for individual (possibly reweighted) * processes. */ map maxWeights; /** * Is set to true when getEvent() is called from skip(int). */ bool skipping; /** * Option for the treatment of the momenta supplied */ unsigned int theMomentumTreatment; /** * Set to true if warnings about possible weight incompatibilities * should be issued. */ bool useWeightWarnings; /** * Option to allow reopening of the file */ bool theReOpenAllowed; /** * Use the spin information */ bool theIncludeSpin; private: /** Access function for the interface. */ void setBeamA(long id); /** Access function for the interface. */ long getBeamA() const; /** Access function for the interface. */ void setBeamB(long id); /** Access function for the interface. */ long getBeamB() const; /** Access function for the interface. */ void setEBeamA(Energy e); /** Access function for the interface. */ Energy getEBeamA() const; /** Access function for the interface. */ void setEBeamB(Energy e); /** Access function for the interface. */ Energy getEBeamB() const; /** Access function for the interface. */ void setPDFA(PDFPtr); /** Access function for the interface. */ PDFPtr getPDFA() const; /** Access function for the interface. */ void setPDFB(PDFPtr); /** Access function for the interface. */ PDFPtr getPDFB() const; private: /** * Describe an abstract base class with persistent data. */ static AbstractClassDescription initFxFxReader; /** * Private and non-existent assignment operator. */ FxFxReader & operator=(const FxFxReader &) = delete; public: /** @cond EXCEPTIONCLASSES */ /** Exception class used by FxFxReader in case inconsistencies * are encountered. */ class FxFxInconsistencyError: public Exception {}; /** Exception class used by FxFxReader in case more events than available are requested. */ class FxFxReopenWarning: public Exception {}; /** Exception class used by FxFxReader in case reopening an event file fails. */ class FxFxReopenError: public Exception {}; /** Exception class used by FxFxReader in case there is information missing in the initialization phase. */ class FxFxInitError: public InitException {}; /** @endcond */ }; /// Stream output for HEPEUP ostream & operator<<(ostream & os, const HEPEUP & h); } #include "ThePEG/Utilities/ClassTraits.h" namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ /** * This template specialization informs ThePEG about the * base class of FxFxReader. */ template <> struct BaseClassTrait: public ClassTraitsType { /** Typedef of the base class of FxFxReader. */ typedef HandlerBase NthBase; }; /** * This template specialization informs ThePEG about the name of the * FxFxReader class and the shared object where it is * defined. */ template <> struct ClassTraits : public ClassTraitsBase { /** * Return the class name. */ static string className() { return "Herwig::FxFxReader"; } /** * Return the name of the shared library to be loaded to get access * to the FxFxReader class and every other class it uses * (except the base class). */ static string library() { return "HwFxFx.so"; } }; /** @endcond */ } #endif /* THEPEG_FxFxReader_H */ diff --git a/MatrixElement/Matchbox/External/VBFNLO/VBFNLOPhasespace.cc b/MatrixElement/Matchbox/External/VBFNLO/VBFNLOPhasespace.cc --- a/MatrixElement/Matchbox/External/VBFNLO/VBFNLOPhasespace.cc +++ b/MatrixElement/Matchbox/External/VBFNLO/VBFNLOPhasespace.cc @@ -1,275 +1,277 @@ // -*- C++ -*- // // VBFNLOPhasespace.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the VBFNLOPhasespace class. // #include "VBFNLOPhasespace.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "Herwig/Utilities/GSLBisection.h" #include "ThePEG/Utilities/DynamicLoader.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/MatrixElement/Matchbox/Base/MatchboxAmplitude.h" #include "VBFNLO/utilities/BLHAinterface.h" #define DEFSTR(s) CPPSTR(s) #define CPPSTR(s) #s using namespace Herwig; VBFNLOPhasespace::VBFNLOPhasespace() : lastSqrtS(0*GeV), needToReshuffle(false), VBFNLOlib_(DEFSTR(VBFNLOLIB)) {} void VBFNLOPhasespace::loadVBFNLO() { if ( ! DynamicLoader::load(VBFNLOlib_+"/libVBFNLO.so") ) { string error1 = DynamicLoader::lastErrorMessage; if ( ! DynamicLoader::load(VBFNLOlib_+"/libVBFNLO.dylib") ) { string error2 = DynamicLoader::lastErrorMessage; if ( ! DynamicLoader::load("libVBFNLO.so") ) { string error3 = DynamicLoader::lastErrorMessage; if ( ! DynamicLoader::load("libVBFNLO.dylib") ) { string error4 = DynamicLoader::lastErrorMessage; throw Exception() << "VBFNLOPhasespace: failed to load libVBFNLO.so/dylib\n" << "Error messages are:\n\n" << "* " << VBFNLOlib_ << "/libVBFNLO.so:\n" << error1 << "\n" << "* " << VBFNLOlib_ << "/libVBFNLO.dylib:\n" << error2 << "\n" << "* libVBFNLO.so:\n" << error3 << "\n" << "* libVBFNLO.dylib:\n" << error4 << "\n" << Exception::runerror; } } } } } VBFNLOPhasespace::~VBFNLOPhasespace() {} IBPtr VBFNLOPhasespace::clone() const { return new_ptr(*this); } IBPtr VBFNLOPhasespace::fullclone() const { return new_ptr(*this); } void VBFNLOPhasespace::setXComb(tStdXCombPtr xco) { MatchboxPhasespace::setXComb(xco); // test for resuffling needToReshuffle = false; if ( xco ) { for ( cPDVector::const_iterator d = mePartonData().begin(); d != mePartonData().end(); ++d ) { // Higgs is massive -> does not need reshuffling if ( ( (**d).id() != ParticleID::h0 ) && ( (**d).hardProcessMass() != ZERO ) ) { needToReshuffle = true; break; } } } // set CMS energy int pStatus = 0; double zero = 0.0; double value = sqrt(lastXCombPtr()->lastS())/GeV; if (value && (value != lastSqrtS/GeV)) { lastSqrtS = value*GeV; string name = "sqrtS"; OLP_SetParameter(const_cast(name.c_str()),&value,&zero,&pStatus); if ( !pStatus ) throw Exception() << "VBFNLOPhasespace::setXComb(): VBFNLO failed to set parameter '" << name << "' to " << value << "\n" << Exception::runerror; } } double VBFNLOPhasespace::generateTwoToNKinematics(const double* random, vector& momenta) { double weight; int id = olpId()[ProcessType::oneLoopInterference] ? olpId()[ProcessType::oneLoopInterference] : olpId()[ProcessType::treeME2]; double* p = new double[4*momenta.size()]; OLP_PhaseSpacePoint(&id, const_cast(random), const_cast(random+1), p, &weight); if (weight < 0) { throw Exception() << "VBFNLOPhasespace::generateTwoToNKinematics(): Negative weight in VBFNLOPhaseSpace\n" << Exception::runerror; } if (weight == 0) { delete[] p; return 0; } for ( size_t i = 0; i < momenta.size(); ++i ) { momenta[i].setT(p[4*i] *GeV); momenta[i].setX(p[4*i+1]*GeV); momenta[i].setY(p[4*i+2]*GeV); momenta[i].setZ(p[4*i+3]*GeV); momenta[i].rescaleMass(); } delete[] p; Energy beamenergy = sqrt(lastXCombPtr()->lastS())/2.; double x1 = momenta[0].e()/beamenergy; double x2 = momenta[1].e()/beamenergy; Energy2 thisSHat = (momenta[0] + momenta[1]).m2(); // reshuffle so that particles have correct mass if ( needToReshuffle ) { // boost final-state into partonic CMS Boost toCMS = (momenta[0]+momenta[1]).findBoostToCM(); for ( size_t i = 2; i < momenta.size(); ++i ) { momenta[i].boost(toCMS); } // copied from MatchboxRambo phasespace double xi; ReshuffleEquation solve(sqrt(thisSHat),mePartonData().begin()+2,mePartonData().end(), momenta.begin()+2,momenta.end()); GSLBisection solver(1e-10,1e-8,10000); try { xi = solver.value(solve,0.0,1.1); } catch (GSLBisection::GSLerror) { return 0.; } catch (GSLBisection::IntervalError) { return 0.; } weight *= pow(xi,3.*(momenta.size()-3.)); Energy num = ZERO; Energy den = ZERO; cPDVector::const_iterator d = mePartonData().begin()+2; for ( vector::iterator k = momenta.begin()+2; k != momenta.end(); ++k, ++d ) { num += (*k).vect().mag2()/(*k).t(); Energy q = (*k).t(); (*k).setT(sqrt(sqr((**d).hardProcessMass())+xi*xi*sqr((*k).t()))); (*k).setVect(xi*(*k).vect()); weight *= q/(*k).t(); den += (*k).vect().mag2()/(*k).t(); (*k).setMass((**d).hardProcessMass()); } // unboost for ( size_t i = 2; i < momenta.size(); ++i ) { momenta[i].boost(-toCMS); } } if ( !matchConstraints(momenta) ) return 0.; lastXCombPtr()->lastX1X2(make_pair(x1,x2)); lastXCombPtr()->lastSHat(thisSHat); weight /= pow(thisSHat/GeV2,momenta.size()-4); weight /= x1*x2; + fillDiagramWeights(); + return weight; } int VBFNLOPhasespace::nDimPhasespace(int nFinal) const { return 3*nFinal; //get this from within VBFNLO int pStatus = 0; double value, zero; string name = "PSdimension"; OLP_GetParameter(const_cast(name.c_str()),&value,&zero,&pStatus); if ( pStatus != 1) { throw Exception() << "VBFNLOPhasespace::nDimPhasespace(): Cannot get phasespace dimension in VBFNLOPhaseSpace\n" << "error code: " << pStatus << "\n" << Exception::runerror; } // one additional number (first) needed for channel selection // one additional number (last) needed for global phi integration return value+2; } Energy VBFNLOPhasespace::ReshuffleEquation::operator() (double xi) const { cPDVector::const_iterator d = dataBegin; vector::const_iterator p = momentaBegin; Energy res = -w; for ( ; d != dataEnd; ++d, ++p ) { res += sqrt(sqr((**d).hardProcessMass()) + xi*xi*sqr(p->t())); } return res; } void VBFNLOPhasespace::doinit() { loadVBFNLO(); MatchboxPhasespace::doinit(); } void VBFNLOPhasespace::doinitrun() { loadVBFNLO(); MatchboxPhasespace::doinitrun(); } void VBFNLOPhasespace::persistentOutput(PersistentOStream & os) const { os << needToReshuffle << theLastXComb; } void VBFNLOPhasespace::persistentInput(PersistentIStream & is, int) { is >> needToReshuffle >> theLastXComb; } // *** Attention *** The following static variable is needed for the type // description system in ThePEG. Please check that the template arguments // are correct (the class and its base class), and that the constructor // arguments are correct (the class name and the name of the dynamically // loadable library where the class implementation can be found). DescribeClass describeHerwigVBFNLOPhasespace("Herwig::VBFNLOPhasespace", "HwMatchboxVBFNLO.so"); void VBFNLOPhasespace::Init() { static ClassDocumentation documentation ("VBFNLOPhasespace is an interface to the internal phasespace generator " "of VBFNLO. It uses the information passed via the BLHA interface to " "obtain information on the required channels."); }