diff --git a/Handlers/StandardEventHandler.cc b/Handlers/StandardEventHandler.cc --- a/Handlers/StandardEventHandler.cc +++ b/Handlers/StandardEventHandler.cc @@ -1,587 +1,592 @@ // -*- C++ -*- // // StandardEventHandler.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 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 StandardEventHandler class. // #include "StandardEventHandler.h" #include "ThePEG/Handlers/StandardXComb.h" #include "ThePEG/Handlers/StdXCombGroup.h" #include "ThePEG/Handlers/SubProcessHandler.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/RefVector.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Utilities/SimplePhaseSpace.h" #include "ThePEG/Utilities/LoopGuard.h" #include "ThePEG/Utilities/Debug.h" #include "ThePEG/PDF/PartonExtractor.h" #include "ThePEG/MatrixElement/MEBase.h" #include "ThePEG/MatrixElement/MEGroup.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Handlers/LuminosityFunction.h" #include "ThePEG/Handlers/SamplerBase.h" #include "ThePEG/Handlers/CascadeHandler.h" #include "ThePEG/Cuts/Cuts.h" #include "ThePEG/Config/algorithm.h" #include #include using namespace ThePEG; StandardEventHandler::StandardEventHandler() : EventHandler(false), collisionCuts(true), theLumiDim(0) { setupGroups(); } StandardEventHandler::~StandardEventHandler() {} void StandardEventHandler::reject(double w) { tStdXCombPtr last = dynamic_ptr_cast(lastXCombPtr()); if ( !last ) return; last->reject(w); xSecStats.reject(w); } void StandardEventHandler::reweight(double factor) const { tStdXCombPtr last = dynamic_ptr_cast(lastXCombPtr()); if ( !currentEvent() || !last ) return; double weight = currentEvent()->weight(); last->reweight(weight,factor*weight); xSecStats.reweight(weight,factor*weight); currentEvent()->weight(factor*weight); for ( map::iterator w = currentEvent()->optionalWeights().begin(); w != currentEvent()->optionalWeights().end(); ++w ) w->second *= factor; } void StandardEventHandler::doupdate() { EventHandler::doupdate(); bool redo = touched(); UpdateChecker::check(theIncomingA, redo); UpdateChecker::check(theIncomingB, redo); for_each(subProcesses(), UpdateChecker(redo)); if ( !redo ) return; theIncoming.first = theIncomingA; theIncoming.second = theIncomingB; for ( SubHandlerList::iterator sit = subProcesses().begin(); sit != subProcesses().end(); ++sit ) if ( !(**sit).pExtractor()->canHandle(incoming()) ) throw StandardEventHandlerUpdateException() << "Cannot use the parton extractor '" << (**sit).pExtractor()->name() << "' in the SubProcessHandler '" << (**sit).name() << "' in the " << "StandardEventHandler '" << name() << "' since it cannot handle " << "the requested types of incoming particles (" << theIncomingA->name() << "," << theIncomingB->name() << ")."; } void StandardEventHandler::doinit() { EventHandler::doinit(); if ( !lumiFnPtr() ) throw StandardEventHandlerUpdateException() << "The StandardEventHandler '" << name() << "' does not have any " << "LuminosityFunction object assigned to it, which it needs to be " << "able to generate events." << Exception::warning; } IBPtr StandardEventHandler::clone() const { return new_ptr(*this); } IBPtr StandardEventHandler::fullclone() const { return new_ptr(*this); } void StandardEventHandler:: addME(Energy maxEnergy, tSubHdlPtr sub, tPExtrPtr extractor, tCutsPtr cuts, tCascHdlPtr ckkw, tMEPtr me, const PBPair & pBins, const PartonPairVec& allPBins) { typedef MEBase::DiagramVector DiagramVector; typedef map DiagramMap; cPDPair pin(pBins.first->parton(), pBins.second->parton()); DiagramVector diag = me->diagrams(); DiagramMap tdiag; DiagramMap tmdiag; for ( int i = 0, N = diag.size(); i < N; ++i ) { cPDPair din(diag[i]->partons()[0], diag[i]->partons()[1]); if (!me->noMirror()) if ( din.first->id() < din.second->id() ) swap(din.first, din.second); if ( din == pin ) tdiag[diag[i]->getTag()].push_back(diag[i]); if (!me->noMirror()) if ( din.first == pin.second && din.second == pin.first ) tmdiag[diag[i]->getTag()].push_back(diag[i]); } if ( tdiag.empty() ) tdiag = tmdiag; for ( DiagramMap::iterator dit = tdiag.begin(); dit != tdiag.end(); ++dit ) { cPDPair din(dit->second.back()->partons()[0], dit->second.back()->partons()[1]); // check assert(me->noMirror() ? din == pin : true); StdXCombPtr xcomb = me->makeXComb(maxEnergy, incoming(), this, sub, extractor, ckkw, pBins, cuts, dit->second, din != pin, allPBins); if ( xcomb->checkInit() ) xCombs().push_back(xcomb); else generator()->logWarning( StandardEventHandlerInitError() << "The matrix element '" << xcomb->matrixElement()->name() << "' cannot generate the diagram '" << dit->first << "' when used together with the parton extractor '" << xcomb->pExtractor()->name() << "'.\nThe corresponding diagram is switched off, " << "check the collision energy and/or the cuts." << Exception::warning); } } void StandardEventHandler::initGroups() { tStdXCombPtr lastXC = dynamic_ptr_cast(lastXCombPtr()); if ( lastXC ) optGroups = lastXC->subProcessHandler()->groups(); EventHandler::initGroups(); } tCollPtr StandardEventHandler::performCollision() { tStdXCombPtr lastXC = dynamic_ptr_cast(lastXCombPtr()); if ( CKKWHandler() ) CKKWHandler()->setXComb(lastXCombPtr()); lastExtractor()->select(lastXC); currentCollision(new_ptr(Collision(lastParticles(), currentEvent(), this))); if ( currentEvent() ) currentEvent()->addCollision(currentCollision()); currentStep(new_ptr(Step(currentCollision()))); currentCollision()->addStep(currentStep()); currentStep()->addSubProcess(lastXC->construct()); + if ( currentEvent() ) { + map optionalWeights = lastXC->generateOptionalWeights(); + for ( const auto& weight : optionalWeights ) + currentEvent()->optionalWeights().insert(weight); + } lastExtractor()->construct(lastXC->partonBinInstances(), currentStep()); if ( collisionCuts ) if ( !lastCuts().passCuts(*currentCollision()) ) throw Veto(); initGroups(); if ( ThePEG_DEBUG_ITEM(1) ) { if ( currentEvent() ) generator()->logfile() << *currentEvent(); else generator()->logfile() << *currentCollision(); } return continueCollision(); } void StandardEventHandler::setScale(Energy2 scale) { lastXCombPtr()->lastScale(scale); } void StandardEventHandler::initialize() { theLumiDim = lumiFn().nDim(incoming()); Energy maxEnergy = lumiFn().maximumCMEnergy(); xCombs().clear(); cuts()->initialize(sqr(maxEnergy), lumiFn().Y()); for ( SubHandlerList::const_iterator sit = subProcesses().begin(); sit != subProcesses().end(); ++sit ) { CutsPtr kincuts = (**sit).cuts()? (**sit).cuts(): cuts(); if ( (**sit).cuts() ) kincuts->initialize(sqr(maxEnergy), lumiFn().Y()); PExtrPtr pextract = (**sit).pExtractor(); tCascHdlPtr ckkw = (**sit).CKKWHandler(); if ( !ckkw ) ckkw = CKKWHandler(); PartonPairVec vpc = pextract->getPartons(maxEnergy, incoming(), *kincuts); // The last parton bin pair was in fact the bins corresponding to // the incoming particles, so we remove them, but save them to // keep them referenced. PBPair orig = vpc.back(); vpc.pop_back(); for ( PartonPairVec::iterator ppit = vpc.begin(); ppit != vpc.end(); ++ppit ) for ( MEVector::const_iterator meit = (**sit).MEs().begin(); meit != (**sit).MEs().end(); ++meit ) { addME(maxEnergy, *sit, pextract, kincuts, ckkw, *meit, *ppit,vpc); } } theMaxDims.clear(); for ( int i = 0, N = xCombs().size(); i < N; ++i ) theMaxDims.push_back(xCombs()[i]->nDim()); sampler()->setEventHandler(this); sampler()->initialize(); } CrossSection StandardEventHandler:: dSigDR(const pair ll, Energy2 maxS, int ibin, int nr, const double * r) { PPair inc = make_pair(incoming().first->produceParticle(), incoming().second->produceParticle()); SimplePhaseSpace::CMS(inc, maxS); xCombs()[ibin]->prepare(inc); return xCombs()[ibin]->dSigDR(ll, nr, r); } tStdXCombPtr StandardEventHandler::select(int bin, double & weight) { tStdXCombPtr lastXC = xCombs()[bin]; // clean up the old XComb object before switching to a new one if ( theLastXComb && theLastXComb != lastXC ) theLastXComb->clean(); theLastXComb = lastXC; weight /= lastXC->matrixElement()->preWeight(); lastXC->select(weight); xSecStats.select(weight); lastXC->accept(); xSecStats.accept(); return lastXC; } int StandardEventHandler::nBins() const { return xCombs().size(); } void StandardEventHandler::statistics(ostream & os) const { if ( statLevel() == 0 ) return; map partonMap; map meMap; map extractMap; XSecStat tot(sampler()->maxXSec()); for ( int i = 0, N = xCombs().size(); i < N; ++i ) { const StandardXComb & x = *xCombs()[i]; if ( partonMap.find(x.partons()) == partonMap.end() ) partonMap[x.partons()] = XSecStat(sampler()->maxXSec()); partonMap[x.partons()] += x.stats(); if ( meMap.find(x.matrixElement()) == meMap.end() ) meMap[x.matrixElement()] = XSecStat(sampler()->maxXSec()); meMap[x.matrixElement()] += x.stats(); if ( extractMap.find(x.pExtractor()) == extractMap.end() ) extractMap[x.pExtractor()] = XSecStat(sampler()->maxXSec()); extractMap[x.pExtractor()] += x.stats(); tot += x.stats(); } string line = string(78, '=') + "\n"; if ( tot.accepted() <= 0 ) { os << line << "No events generated by event handler '" << name() << "'." << endl; return; } os << line << "Statistics for event handler \'" << name() << "\':\n" << " " << "generated number of Cross-section\n" << " " << " events attempts (nb)\n"; os << line << "Total (from attempted events): including vetoed events" << setw(23) << ouniterr(sampler()->integratedXSec(), sampler()->integratedXSecErr(), nanobarn) << endl; os << line << "Total (from generated events):" << setw(17) << tot.accepted() << setw(13) << tot.attempts() << setw(17) << ouniterr(tot.xSec(sampler()->attempts()),tot.xSecErr(sampler()->attempts()) , nanobarn) << "\n"; os << "Events carry "; if ( weighted() ) os << "varying weights."; else if ( sampler()->almostUnweighted() ) os << "varying weights, most of which are unit weights."; else os << "unit weights."; os << endl << line; if ( statLevel() == 1 ) return; os << "Per matrix element breakdown:\n"; for ( map::iterator i = meMap.begin(); i != meMap.end(); ++i ) { string n = i->first->name(); n.resize(37, ' '); os << n << setw(11) << i->second.accepted() << setw(13) << i->second.attempts() << setw(17) << ouniterr(i->second.xSec(sampler()->attempts()), i->second.xSecErr(sampler()->attempts()), nanobarn) << endl; } os << line; if ( statLevel() == 2 ) return; os << "Per parton extractor breakdown:\n"; for ( map::iterator i = extractMap.begin(); i != extractMap.end(); ++i ) { string n = i->first->name(); n.resize(37, ' '); os << n << setw(11) << i->second.accepted() << setw(13) << i->second.attempts() << setw(17) << ouniterr(i->second.xSec(sampler()->attempts()), i->second.xSecErr(sampler()->attempts()), nanobarn) << endl; } os << line; os << "Per incoming partons breakdown:\n"; for ( map::iterator i = partonMap.begin(); i != partonMap.end(); ++i ) { string n = i->first.first->PDGName() + " " + i->first.second->PDGName(); n.resize(37, ' '); os << n << setw(11) << i->second.accepted() << setw(13) << i->second.attempts() << setw(17) << ouniterr(i->second.xSec(sampler()->attempts()), i->second.xSecErr(sampler()->attempts()), nanobarn) << endl; } os << line; if ( statLevel() == 3 ) return; os << "Detailed breakdown:\n"; for ( int i = 0, N = xCombs().size(); i < N; ++i ) { const StandardXComb & x = *xCombs()[i]; XSecStat xstat(sampler()->maxXSec()); xstat += x.stats(); os << "(" << x.pExtractor()->name() << ") " << x.partons().first->PDGName() << " " << x.partons().second->PDGName() << " (" << x.matrixElement()->name() << " " << x.lastDiagram()->getTag() << ") " << endl << setw(48) << xstat.accepted() << setw(13) << xstat.attempts() << setw(17) << ouniterr(xstat.xSec(sampler()->attempts()), xstat.xSecErr(sampler()->attempts()), nanobarn) << endl; } os << line; } CrossSection StandardEventHandler::histogramScale() const { xSecStats.maxXSec(sampler()->maxXSec()); return xSecStats.xSec(sampler()->attempts())/xSecStats.sumWeights(); } CrossSection StandardEventHandler::integratedXSec() const { xSecStats.maxXSec(sampler()->maxXSec()); return xSecStats.xSec(sampler()->attempts()); } CrossSection StandardEventHandler::integratedXSecErr() const { xSecStats.maxXSec(sampler()->maxXSec()); return xSecStats.xSecErr(sampler()->attempts()); } CrossSection StandardEventHandler::integratedXSecNoReweight() const { xSecStats.maxXSec(sampler()->maxXSec()); return xSecStats.xSecNoReweight(sampler()->attempts()); } CrossSection StandardEventHandler::integratedXSecErrNoReweight() const { xSecStats.maxXSec(sampler()->maxXSec()); return xSecStats.xSecErrNoReweight(sampler()->attempts()); } void StandardEventHandler::doinitrun() { EventHandler::doinitrun(); for ( SubHandlerList::iterator sit = subProcesses().begin(); sit != subProcesses().end(); ++sit ) (**sit).initrun(); sampler()->initrun(); for ( int i = 0, N = xCombs().size(); i < N; ++i ) xCombs()[i]->reset(); xSecStats.reset(); } CrossSection StandardEventHandler::dSigDR(const vector & r) { double jac = 1.0; pair ll = lumiFn().generateLL(&r[0], jac); Energy2 maxS = sqr(lumiFn().maximumCMEnergy())/exp(ll.first + ll.second); int bin = sampler()->lastBin(); CrossSection x = jac*lumiFn().value(incoming(), ll.first, ll.second) *dSigDR(ll, maxS, bin, nDim(bin) - lumiDim(), &r[lumiDim()]); return x; } EventPtr StandardEventHandler::generateEvent() { LoopGuard loopGuard(*this, maxLoop()); while (1) { loopGuard(); EventHandler::clean(); double weight = sampler()->generate(); tStdXCombPtr lastXC = select(sampler()->lastBin(), weight); try { lumiFn().select(lastXC); currentEventBoost() = lumiFn().getBoost(); currentEvent(new_ptr(Event(lastParticles(), this, generator()->runName(), generator()->currentEventNumber(), weight))); performCollision(); if ( !currentCollision() ) throw Veto(); currentEvent()->transform(currentEventBoost()); return currentEvent(); } catch (Veto) { reject(currentEvent()->weight()); } catch (Stop) { break; } catch (Exception &) { reject(currentEvent()->weight()); throw; } } return currentEvent(); } EventPtr StandardEventHandler::continueEvent() { if ( !generator() ) throw StandardEventHandlerInitError() << "The event handler '" << name() << "' had not been isolated " << "from the setup phase before it was used." << Exception::maybeabort; try { continueCollision(); } catch (Veto) { reject(currentEvent()->weight()); } catch (Stop) { } catch (Exception &) { reject(currentEvent()->weight()); throw; } return currentEvent(); } void StandardEventHandler::select(tXCombPtr newXComb) { EventHandler::select(newXComb); lastExtractor()->select(newXComb); } void StandardEventHandler::clean() { if ( theLastXComb ) theLastXComb->clean(); for (size_t i=0; i < theXCombs.size(); ++i ) if ( theXCombs[i] ) theXCombs[i]->clean(); EventHandler::clean(); } void StandardEventHandler::dofinish() { clean(); EventHandler::dofinish(); } ClassDescription StandardEventHandler::initStandardEventHandler; // void StandardEventHandler::setIncomingA(PDPtr p) { theIncomingA = p; theIncoming.first = p; } void StandardEventHandler::setIncomingB(PDPtr p) { theIncomingB = p; theIncoming.second = p; } void StandardEventHandler::Init() { static ClassDocumentation documentation ("This is the standard event handler to generate hard sub-processes " "within ThePEG. It must specify a pair of incoming particle beams " "in BeamA and BeamB " "and a suiteable LuminosityFunction. In " "addition at least one object describing the sub-processes to be " "generated must be specified in " "SubProcessHandlers."); static Reference interfaceIncomingA ("BeamA", "The type of particles in first beam", &StandardEventHandler::theIncomingA, false, false, true, false, &StandardEventHandler::setIncomingA); static Reference interfaceIncomingB ("BeamB", "The type of particles in second beam", &StandardEventHandler::theIncomingB, false, false, true, false, &StandardEventHandler::setIncomingB); static RefVector interfaceSubhandlers ("SubProcessHandlers", "The list of sub-process handlers used in this StandardEventHandler. ", &StandardEventHandler::theSubProcesses, 0, false, false, true, false); static Reference interfaceCuts ("Cuts", "Common kinematical cuts for this StandardEventHandler. These cuts " "may be overidden in individual sub-process handlers.", &StandardEventHandler::theCuts, false, false, true, false); static Switch interfaceCollisionCuts ("CollisionCuts", "Switch on or off cuts on collision objects", &StandardEventHandler::collisionCuts, true, false, false); static SwitchOption interfaceCollisionCutsOn (interfaceCollisionCuts, "On", "Switch on cuts on collision objects", true); static SwitchOption interfaceCollisionCutsOff (interfaceCollisionCuts, "Off", "Switch off cuts on collision cuts", false); static SwitchOption interfaceCollisionCutsYes (interfaceCollisionCuts, "Yes", "Switch on cuts on collision objects", true); static SwitchOption interfaceCollisionCutsNo (interfaceCollisionCuts, "No", "Switch off cuts on collision cuts", false); static Reference interfaceSampler ("Sampler", "The phase space sampler responsible for generating phase space" "points according to the cross section given by this event handler", &StandardEventHandler::theSampler, false, false, true, true); interfaceSubhandlers.rank(11); interfaceIncomingA.rank(3); interfaceIncomingB.rank(2); } void StandardEventHandler::persistentOutput(PersistentOStream & os) const { os << theIncomingA << theIncomingB << theSubProcesses << theCuts << collisionCuts << theXCombs << theMaxDims << theSampler << theLumiDim << xSecStats; } void StandardEventHandler::persistentInput(PersistentIStream & is, int) { is >> theIncomingA >> theIncomingB >> theSubProcesses >> theCuts >> collisionCuts >> theXCombs >> theMaxDims >> theSampler >> theLumiDim >> xSecStats; } diff --git a/Handlers/StandardXComb.cc b/Handlers/StandardXComb.cc --- a/Handlers/StandardXComb.cc +++ b/Handlers/StandardXComb.cc @@ -1,828 +1,833 @@ // -*- C++ -*- // // StandardXComb.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 Leif Lonnblad // Copyright (C) 2009-2017 Simon Platzer // // 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 StandardXComb class. // #include "StandardXComb.h" #include "StdXCombGroup.h" #include "ThePEG/Handlers/StandardEventHandler.h" #include "ThePEG/Handlers/SubProcessHandler.h" #include "ThePEG/Cuts/Cuts.h" #include "ThePEG/PDF/PartonExtractor.h" #include "ThePEG/Utilities/Debug.h" #include "ThePEG/Utilities/Maths.h" #include "ThePEG/PDT/ParticleData.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Utilities/SimplePhaseSpace.h" #include "ThePEG/Utilities/UtilityBase.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/EventRecord/SubProcessGroup.h" #include "ThePEG/Vectors/LorentzRotation.h" #include "ThePEG/MatrixElement/MEBase.h" #include "ThePEG/MatrixElement/ColourLines.h" #include "ThePEG/Handlers/LuminosityFunction.h" #include "ThePEG/Handlers/CascadeHandler.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/EventRecord/TmpTransform.h" #ifdef ThePEG_TEMPLATES_IN_CC_FILE #include "StandardXComb.tcc" #endif using namespace ThePEG; StandardXComb::StandardXComb() : XComb(), isMirror(false), theNDim(0), partonDims(make_pair(0, 0)), theKinematicsGenerated(false), theLastDiagramIndex(0), theLastPDFWeight(0.0), theLastCrossSection(ZERO), theLastJacobian(1.0), theLastME2(-1.0), theLastPreweight(1.0), theLastMECrossSection(ZERO), theLastMEPDFWeight(1.0), theLastMECouplings(1.0), checkedCuts(false), passedCuts(false), theCutWeight(1.0), theNeedsReshuffling(false) {} StandardXComb:: StandardXComb(Energy newMaxEnergy, const cPDPair & inc, tEHPtr newEventHandler, tSubHdlPtr newSubProcessHandler, tPExtrPtr newExtractor, tCascHdlPtr newCKKW, const PBPair & newPartonBins, tCutsPtr newCuts, tMEPtr newME, const DiagramVector & newDiagrams, bool mir, tStdXCombPtr newHead) : XComb(newMaxEnergy, inc, newEventHandler, newExtractor, newCKKW, newPartonBins, newCuts), theSubProcessHandler(newSubProcessHandler), theME(newME), theDiagrams(newDiagrams), isMirror(mir), theNDim(0), partonDims(0,0), theKinematicsGenerated(false), theLastDiagramIndex(0), theLastPDFWeight(0.0), theLastCrossSection(ZERO), theLastJacobian(0.0), theLastME2(-1.0), theLastPreweight(1.0), theLastMECrossSection(ZERO), theLastMEPDFWeight(1.0), theLastMECouplings(1.0), theHead(newHead), checkedCuts(false), passedCuts(false), theCutWeight(1.0), theNeedsReshuffling(false) { partonDims = pExtractor()->nDims(partonBins()); if ( matrixElement()->haveX1X2() ) { partonDims.first = 0; partonDims.second = 0; } theNDim = matrixElement()->nDim() + partonDims.first + partonDims.second; mePartonData() = lastDiagram()->partons(); checkReshufflingNeeds(); } StandardXComb::StandardXComb(tMEPtr me, const tPVector & parts, DiagramIndex indx) : theME(me), isMirror(false), theNDim(0), partonDims(make_pair(0, 0)), theKinematicsGenerated(false), theLastDiagramIndex(0), theLastPDFWeight(0.0), theLastCrossSection(ZERO), theLastME2(-1.0), theLastPreweight(1.0), theLastMECrossSection(ZERO), theLastMEPDFWeight(1.0), theLastMECouplings(1.0), checkedCuts(false), passedCuts(false), theCutWeight(1.0), theNeedsReshuffling(false) { subProcess(new_ptr(SubProcess(make_pair(parts[0], parts[1]), tCollPtr(), me))); for ( int i = 0, N = parts.size(); i < N; ++i ) { subProcess()->addOutgoing(parts[i], false); theMEPartonData.push_back(parts[i]->dataPtr()); theMEMomenta.push_back(parts[i]->momentum()); } lastSHat((meMomenta()[0] + meMomenta()[1]).m2()); string tag = me->diagrams()[indx]->getTag(); for ( int i = 0, N = me->diagrams().size(); i < N; ++i ) if ( me->diagrams()[i]->getTag() == tag ) theDiagrams.push_back(me->diagrams()[i]); checkReshufflingNeeds(); } StandardXComb::StandardXComb(tStdXCombPtr newHead, const PBPair & newPartonBins, tMEPtr newME, const DiagramVector & newDiagrams) : XComb(newHead->maxEnergy(), newHead->particles(), newHead->eventHandlerPtr(), newHead->pExtractor(), newHead->CKKWHandler(), newPartonBins, newHead->cuts()), theSubProcessHandler(const_ptr_cast(newHead->subProcessHandler())), theME(newME), theDiagrams(newDiagrams), isMirror(newHead->mirror()), theNDim(0), partonDims(0,0), theKinematicsGenerated(false), theLastDiagramIndex(0), theLastPDFWeight(0.0), theLastCrossSection(ZERO), theLastJacobian(0.0), theLastME2(-1.0), theLastPreweight(1.0), theLastMECrossSection(ZERO), theLastMEPDFWeight(1.0), theLastMECouplings(1.0), theHead(newHead), checkedCuts(false), passedCuts(false), theCutWeight(1.0), theNeedsReshuffling(false) { partonDims = pExtractor()->nDims(partonBins()); if ( matrixElement()->haveX1X2() ) { partonDims.first = 0; partonDims.second = 0; } theNDim = matrixElement()->nDim() + partonDims.first + partonDims.second; mePartonData() = lastDiagram()->partons(); checkReshufflingNeeds(); /* // wait for C++11 StandardXComb(newHead->maxEnergy(),newHead->particles(), newHead->eventHandlerPtr(), const_ptr_cast(newHead->subProcessHandler()), newHead->pExtractor(),newHead->CKKWHandler(), newPartonBins,newHead->cuts(),newME,newDiagrams,newHead->mirror(), newHead); */ } StandardXComb::~StandardXComb() {} void StandardXComb::recreatePartonBinInstances(Energy2 scale) { PBIPair newBins; Direction<0> dir(true); newBins.first = new_ptr(PartonBinInstance(lastPartons().first,partonBins().first,scale)); dir.reverse(); newBins.second = new_ptr(PartonBinInstance(lastPartons().second,partonBins().second,scale)); resetPartonBinInstances(newBins); setPartonBinInfo(); lastPartons().first->scale(partonBinInstances().first->scale()); lastPartons().second->scale(partonBinInstances().second->scale()); } void StandardXComb::refillPartonBinInstances(const double* r) { pExtractor()->select(this); pExtractor()->updatePartonBinInstances(partonBinInstances()); pExtractor()->generateSHat(lastS(), partonBinInstances(), r, r + nDim() - partonDims.second,true); } bool StandardXComb::setIncomingPartons(tStdXCombPtr labHead) { if ( lastPartons().first ) return true; // If no labHead was given, we need to search the labhead if ( !labHead ){ labHead = head(); while ( labHead->head() && labHead->head() != labHead ) { labHead =labHead->head(); } } createPartonBinInstances(); lastParticles(labHead->lastParticles()); setPartonBinInfo(); lastPartons(make_pair(mePartonData()[0]->produceParticle(Lorentz5Momentum()), mePartonData()[1]->produceParticle(Lorentz5Momentum()))); Lorentz5Momentum pFirst = meMomenta()[0]; Lorentz5Momentum pSecond = meMomenta()[1]; if ( labHead->matrixElement()->wantCMS() ) { Boost toLab = (labHead->lastPartons().first->momentum() + labHead->lastPartons().second->momentum()).boostVector(); if ( toLab.mag2() > Constants::epsilon ) { pFirst.boost(toLab); pSecond.boost(toLab); } } lastPartons().first->set5Momentum(pFirst); lastPartons().second->set5Momentum(pSecond); partonBinInstances().first->parton(lastPartons().first); partonBinInstances().second->parton(lastPartons().second); lastS((lastParticles().first->momentum() + lastParticles().second->momentum()).m2()); lastSHat((lastPartons().first->momentum() + lastPartons().second->momentum()).m2()); lastP1P2(make_pair(labHead->lastP1(),labHead->lastP2())); double x1 = lastPartons().first->momentum().plus()/ lastParticles().first->momentum().plus(); double x2 = lastPartons().second->momentum().minus()/ lastParticles().second->momentum().minus(); if(x1<=0. || x2 <= 0. ) return false; lastX1X2(make_pair(x1,x2)); lastY((lastPartons().first->momentum()+ lastPartons().second->momentum()).rapidity()); return true; } void StandardXComb::fill(const PPair& newParticles, const PPair& newPartons, const vector& newMEMomenta, const DVector& newLastRandomNumbers) { lastParticles(newParticles); lastP1P2(make_pair(0.0, 0.0)); lastS((lastParticles().first->momentum() + lastParticles().second->momentum()).m2()); lastPartons(newPartons); lastSHat((lastPartons().first->momentum() + lastPartons().second->momentum()).m2()); lastX1X2(make_pair(lastPartons().first->momentum().plus()/ lastParticles().first->momentum().plus(), lastPartons().second->momentum().minus()/ lastParticles().second->momentum().minus())); lastY((lastPartons().first->momentum() + lastPartons().second->momentum()).rapidity()); meMomenta().resize(newMEMomenta.size()); copy(newMEMomenta.begin(),newMEMomenta.end(), meMomenta().begin()); lastRandomNumbers().resize(newLastRandomNumbers.size()); copy(newLastRandomNumbers.begin(),newLastRandomNumbers.end(), lastRandomNumbers().begin()); } bool StandardXComb::checkInit() { Energy summin = ZERO; for ( int i = 2, N = mePartonData().size(); i < N; ++i ) { summin += mePartonData()[i]->massMin(); } return ( summin < min(maxEnergy(), cuts()->mHatMax()) ); } bool StandardXComb::willPassCuts() { if ( checkedCuts ) return passedCuts; checkedCuts = true; if ( !head() ) { if ( !cuts()->initSubProcess(lastSHat(), lastY(), mirror()) ) { passedCuts = false; theCutWeight = 0.0; return false; } } else { cuts()->initSubProcess(lastSHat(), lastY(), mirror()); } // check for exceptional configurations which may happen in NLO // dijets subtraction with an extremely soft incoming parton giving // rise to about lightlike CM momentum if ( (meMomenta()[0]+meMomenta()[1]).m2() <= ZERO ) { passedCuts = false; return false; } tcPDVector outdata(mePartonData().begin()+2,mePartonData().end()); vector outmomenta(meMomenta().begin()+2,meMomenta().end()); Boost tocm = (meMomenta()[0]+meMomenta()[1]).findBoostToCM(); if ( tocm.mag2() > Constants::epsilon ) { for ( vector::iterator p = outmomenta.begin(); p != outmomenta.end(); ++p ) { p->boost(tocm); } } for ( vector::iterator p = outmomenta.begin(); p != outmomenta.end(); ++p ) if ( !std::isfinite(double(p->x()/GeV)) || !std::isfinite(double(p->y()/GeV)) || !std::isfinite(double(p->z()/GeV)) || !std::isfinite(double(p->t()/GeV)) ) throw Exception() << "Event momenta contain an invalid entry: " << (*p)/GeV << Exception::runerror; if ( !cuts()->passCuts(outdata,outmomenta,mePartonData()[0],mePartonData()[1]) ) { theCutWeight = cuts()->cutWeight(); passedCuts = false; return false; } theCutWeight = cuts()->cutWeight(); passedCuts = true; return true; } void StandardXComb::clean() { XComb::clean(); theLastPDFWeight = 0.0; theLastCrossSection = ZERO; theLastJacobian = 0.0; theLastME2 = 0.0; theLastPreweight = 1.0; theLastMECrossSection = ZERO; theLastMEPDFWeight = 0.0; theLastMECouplings = 0.0; theProjectors.clear(); theProjector = StdXCombPtr(); theKinematicsGenerated = false; checkedCuts = false; passedCuts = false; theCutWeight = 1.0; theExternalDiagram = tcDiagPtr(); matrixElement()->flushCaches(); } CrossSection StandardXComb::dSigDR(const double * r) { matrixElement()->setXComb(this); if ( !matrixElement()->apply() ) { subProcess(SubProPtr()); lastCrossSection(ZERO); return ZERO; } meMomenta().resize(mePartonData().size()); if ( !matrixElement()->generateKinematics(r) ) { subProcess(SubProPtr()); lastCrossSection(ZERO); return ZERO; } setIncomingPartons(); lastScale(matrixElement()->scale()); lastAlphaS(matrixElement()->alphaS()); lastAlphaEM(matrixElement()->alphaEM()); partonBinInstances().first->scale(lastScale()); partonBinInstances().second->scale(lastScale()); if ( (!willPassCuts() && !matrixElement()->headCuts() && !matrixElement()->ignoreCuts()) || !matrixElement()->apply() ) { subProcess(SubProPtr()); lastCrossSection(ZERO); return ZERO; } lastPDFWeight(head()->lastPDFWeight()); matrixElement()->setKinematics(); CrossSection xsec = matrixElement()->dSigHatDR() * lastPDFWeight(); xsec *= cutWeight(); subProcess(SubProPtr()); lastCrossSection(xsec); return xsec; } CrossSection StandardXComb:: dSigDR(const pair ll, int nr, const double * r) { if ( matrixElement()->keepRandomNumbers() ) { lastRandomNumbers().resize(nDim()); copy(r,r+nDim(),lastRandomNumbers().begin()); } pExtractor()->select(this); setPartonBinInfo(); lastP1P2(ll); lastS(sqr(maxEnergy())/exp(lastP1() + lastP2())); meMomenta().resize(mePartonData().size()); matrixElement()->setXComb(this); PPair partons; if ( !matrixElement()->haveX1X2() ) { if ( !pExtractor()->generateL(partonBinInstances(), r, r + nr - partonDims.second) ) { lastCrossSection(ZERO); return ZERO; } partons = make_pair(partonBinInstances().first->parton(), partonBinInstances().second->parton()); lastSHat(lastS()/exp(partonBinInstances().first->l() + partonBinInstances().second->l())); meMomenta()[0] = partons.first->momentum(); meMomenta()[1] = partons.second->momentum(); } else { if ( !matrixElement()->generateKinematics(r + partonDims.first) ) { lastCrossSection(ZERO); return ZERO; } lastSHat((meMomenta()[0]+meMomenta()[1]).m2()); matrixElement()->setKinematics(); lastScale(matrixElement()->scale()); partons.first = mePartonData()[0]->produceParticle(meMomenta()[0]); partons.second = mePartonData()[1]->produceParticle(meMomenta()[1]); Direction<0> dir(true); partonBinInstances().first = new_ptr(PartonBinInstance(lastParticles().first,partons.first, partonBins().first,lastScale())); dir.reverse(); partonBinInstances().second = new_ptr(PartonBinInstance(lastParticles().second,partons.second, partonBins().second,lastScale())); } lastPartons(partons); if ( lastSHat() < cuts()->sHatMin() ) { lastCrossSection(ZERO); return ZERO; } lastY(0.5*(partonBinInstances().second->l() - partonBinInstances().first->l())); if ( !cuts()->initSubProcess(lastSHat(), lastY(), mirror()) ) { lastCrossSection(ZERO); return ZERO; } if ( mirror() ) swap(meMomenta()[0], meMomenta()[1]); if ( matrixElement()->wantCMS() && !matrixElement()->haveX1X2() ) SimplePhaseSpace::CMS(meMomenta()[0], meMomenta()[1], lastSHat()); Energy summ = ZERO; if ( meMomenta().size() == 3 ) { if ( !matrixElement()->haveX1X2() ) meMomenta()[2] = Lorentz5Momentum(sqrt(lastSHat())); } else { for ( int i = 2, N = meMomenta().size(); i < N; ++i ) { if ( !matrixElement()->haveX1X2() ) meMomenta()[i] = Lorentz5Momentum(mePartonData()[i]->mass()); summ += mePartonData()[i]->massMin(); } if ( sqr(summ) >= lastSHat() ) { lastCrossSection(ZERO); return ZERO; } } if ( !matrixElement()->haveX1X2() ) lastScale(max(lastSHat()/4.0, cuts()->scaleMin())); lastSHat(pExtractor()->generateSHat(lastS(), partonBinInstances(), r, r + nr - partonDims.second, matrixElement()->haveX1X2())); if ( !cuts()->sHat(lastSHat()) ) { lastCrossSection(ZERO); return ZERO; } r += partonDims.first; lastX1X2(make_pair(lastPartons().first->momentum().plus()/ lastParticles().first->momentum().plus(), lastPartons().second->momentum().minus()/ lastParticles().second->momentum().minus())); if ( !cuts()->x1(lastX1()) || !cuts()->x2(lastX2()) ) { lastCrossSection(ZERO); return ZERO; } lastY((lastPartons().first->momentum() + lastPartons().second->momentum()).rapidity()); if ( !cuts()->yHat(lastY()) ) { lastCrossSection(ZERO); return ZERO; } if ( !cuts()->initSubProcess(lastSHat(), lastY(), mirror()) ) { lastCrossSection(ZERO); return ZERO; } meMomenta()[0] = lastPartons().first->momentum(); meMomenta()[1] = lastPartons().second->momentum(); if ( mirror() ) swap(meMomenta()[0], meMomenta()[1]); if ( matrixElement()->wantCMS() && !matrixElement()->haveX1X2() ) SimplePhaseSpace::CMS(meMomenta()[0], meMomenta()[1], lastSHat()); if ( meMomenta().size() == 3 ) { if ( !matrixElement()->haveX1X2() ) meMomenta()[2] = Lorentz5Momentum(sqrt(lastSHat())); } else { if ( sqr(summ) >= lastSHat() ) { lastCrossSection(ZERO); return ZERO; } } if ( !matrixElement()->haveX1X2() ) { if ( !matrixElement()->generateKinematics(r) ) { lastCrossSection(ZERO); return ZERO; } } lastScale(matrixElement()->scale()); if ( !cuts()->scale(lastScale()) ) { lastCrossSection(ZERO); return ZERO; } // get information on cuts; we don't take this into account here for // reasons of backward compatibility but this will change eventually willPassCuts(); pair evalPDFS = make_pair(matrixElement()->havePDFWeight1(), matrixElement()->havePDFWeight2()); if ( mirror() ) swap(evalPDFS.first,evalPDFS.second); lastPDFWeight(pExtractor()->fullFn(partonBinInstances(), lastScale(), evalPDFS)); if ( lastPDFWeight() == 0.0 ) { lastCrossSection(ZERO); return ZERO; } matrixElement()->setKinematics(); CrossSection xsec = matrixElement()->dSigHatDR() * lastPDFWeight(); if ( xsec == ZERO ) { lastCrossSection(ZERO); return ZERO; } const bool isCKKW=CKKWHandler() && matrixElement()->maxMultCKKW() > 0 && matrixElement()->maxMultCKKW() > matrixElement()->minMultCKKW(); if(!isCKKW) xsec *= cutWeight(); lastAlphaS (matrixElement()->orderInAlphaS () >0 ? matrixElement()->alphaS() : -1.); lastAlphaEM(matrixElement()->orderInAlphaEW() >0 ? matrixElement()->alphaEM() : -1.); matrixElement()->fillProjectors(); if ( !projectors().empty() ) { lastProjector(projectors().select(UseRandom::rnd())); } subProcess(SubProPtr()); if ( isCKKW ) { newSubProcess(); CKKWHandler()->setXComb(this); xsec *= CKKWHandler()->reweightCKKW(matrixElement()->minMultCKKW(), matrixElement()->maxMultCKKW()); } if ( matrixElement()->reweighted() ) { newSubProcess(); xsec *= matrixElement()->reWeight() * matrixElement()->preWeight(); } lastCrossSection(xsec); return xsec; } +map StandardXComb::generateOptionalWeights() { + matrixElement()->setXComb(this); + return matrixElement()->generateOptionalWeights(); +} + void StandardXComb::newSubProcess(bool group) { if ( subProcess() ) return; if ( head() && matrixElement()->wantCMS() ) { // first get the meMomenta in their CMS, as this may // not be the case Boost cm = (meMomenta()[0] + meMomenta()[1]).findBoostToCM(); if ( cm.mag2() > Constants::epsilon ) { for ( vector::iterator m = meMomenta().begin(); m != meMomenta().end(); ++m ) { *m = m->boost(cm); } } } if ( !lastProjector() ) { if ( !group ) subProcess(new_ptr(SubProcess(lastPartons(), tCollPtr(), matrixElement()))); else subProcess(new_ptr(SubProcessGroup(lastPartons(), tCollPtr(), matrixElement()))); if ( !theExternalDiagram ) lastDiagramIndex(matrixElement()->diagram(diagrams())); const ColourLines & cl = matrixElement()->selectColourGeometry(lastDiagram()); Lorentz5Momentum p1 = lastPartons().first->momentum(); Lorentz5Momentum p2 = lastPartons().second->momentum(); tPPair inc = lastPartons(); if ( mirror() ) swap(inc.first, inc.second); if ( matrixElement()->wantCMS() && !matrixElement()->haveX1X2() ) { LorentzRotation r = Utilities::boostToCM(inc); lastDiagram()->construct(subProcess(), *this, cl); subProcess()->transform(r.inverse()); lastPartons().first->set5Momentum(p1); lastPartons().second->set5Momentum(p2); } else { lastDiagram()->construct(subProcess(), *this, cl); } lastPartons().first ->scale(partonBinInstances().first ->scale()); lastPartons().second->scale(partonBinInstances().second->scale()); for ( int i = 0, N = subProcess()->outgoing().size(); i < N; ++i ) subProcess()->outgoing()[i]->scale(lastScale()); // construct the spin information for the interaction matrixElement()->constructVertex(subProcess(),&cl); // set veto scales matrixElement()->setVetoScales(subProcess()); } else { lastProjector()->newSubProcess(); subProcess(lastProjector()->subProcess()); lastPartons(lastProjector()->lastPartons()); lastSHat((lastPartons().first->momentum() + lastPartons().second->momentum()).m2()); lastX1X2(make_pair(lastPartons().first->momentum().plus()/ lastParticles().first->momentum().plus(), lastPartons().second->momentum().minus()/ lastParticles().second->momentum().minus())); lastY(log(lastX1()/lastX2())*0.5); lastCentralScale(lastProjector()->lastCentralScale()); lastShowerScale(lastProjector()->lastShowerScale()); partonBinInstances().first->parton(lastPartons().first); partonBinInstances().second->parton(lastPartons().second); if ( !matrixElement()->keepRandomNumbers() ) throw Exception() << "Matrix element needs to request random number storage " << "for creating projected subprocesses." << Exception::runerror; const double * r = &lastRandomNumbers()[0]; pExtractor()->generateSHat(lastS(), partonBinInstances(), r, r + nDim() - partonDims.second,true); pExtractor()->updatePartonBinInstances(partonBinInstances()); } } void StandardXComb::checkReshufflingNeeds() { theNeedsReshuffling = false; if ( eventHandler().cascadeHandler() ) { if ( eventHandler().cascadeHandler()->isReshuffling() ) return; } for ( cPDVector::const_iterator p = mePartonData().begin() + 2; p != mePartonData().end(); ++p ) { theNeedsReshuffling |= ( (**p).mass() != (**p).hardProcessMass() ); } } double StandardXComb::reshuffleEquation(double k, const vector >& coeffs, Energy2 Q2) const { double res = 0.; for ( vector >::const_iterator c = coeffs.begin(); c != coeffs.end(); ++c ) { double x = sqr(k)*c->first/Q2+c->second/Q2; if ( x < 0.0 ) throw Veto(); res += sqrt(x); } return res - 1.; } double StandardXComb::solveReshuffleEquation(const vector >& coeffs, Energy2 Q2) const { // shamelessly adapted from Herwig++'s QTildeReconstructor // in principle we could have included the exact solution for two // outgoing particles, but for this to be simple we would require // boosting to the CM so the gain is questionable double k1 = 0.,k2 = 1.,k = 0.; if ( reshuffleEquation(k1,coeffs,Q2) < 0.0 ) { while ( reshuffleEquation(k2, coeffs, Q2) < 0.0 ) { k1 = k2; k2 *= 2; } while ( abs( (k1 - k2)/(k1 + k2) ) > 1.e-10 ) { if( reshuffleEquation(k2,coeffs,Q2) == 0.0 ) { return k2; } else { k = (k1+k2)/2.; if ( reshuffleEquation(k,coeffs,Q2) > 0.0 ) { k2 = k; } else { k1 = k; } } } return k1; } throw Veto(); return 1.; } void StandardXComb::reshuffle(vector& pout) const { Lorentz5Momentum Q(ZERO,ZERO,ZERO,ZERO); for ( vector::const_iterator p = pout.begin(); p != pout.end(); ++p ) Q += *p; if ( Q.m2() <= ZERO ) throw Veto(); LorentzVector Qhat = Q/sqrt(Q.m2()); vector pQhat; vector > coeffs; vector::const_iterator p = pout.begin(); cPDVector::const_iterator pd = mePartonData().begin() + 2; for ( ; p != pout.end(); ++p, ++pd ) { pQhat.push_back((*p)*Qhat); coeffs.push_back(make_pair(sqr(pQhat.back())-sqr((**pd).hardProcessMass()), sqr((**pd).mass()))); } double k = solveReshuffleEquation(coeffs,Q.m2()); vector::iterator px = pout.begin(); vector::const_iterator pQ = pQhat.begin(); vector >::const_iterator coeff = coeffs.begin(); for ( ; px != pout.end(); ++px, ++pQ, ++coeff ) { *px = k*((*px) - (*pQ) * Qhat) + sqrt(sqr(k)*coeff->first + coeff->second)*Qhat; px->rescaleMass(); } } tSubProPtr StandardXComb::construct() { matrixElement()->setXComb(this); if ( !head() ) { if ( !cuts()->initSubProcess(lastSHat(), lastY()) ) throw Veto(); } else { cuts()->initSubProcess(lastSHat(), lastY()); } if ( head() ) if ( !matrixElement()->apply() ) return tSubProPtr(); setPartonBinInfo(); matrixElement()->setKinematics(); newSubProcess(); TmpTransform tmp(subProcess(), Utilities::getBoostToCM(subProcess()->incoming())); if ( !cuts()->passCuts(*subProcess()) ) throw Veto(); if ( head() ) { subProcess()->head(head()->subProcess()); if(lastCrossSection()!=ZERO&&head()->lastCrossSection()!=ZERO) subProcess()->groupWeight(lastCrossSection()/head()->lastCrossSection()); else subProcess()->groupWeight(0.); } return subProcess(); } void StandardXComb::Init() {} void StandardXComb::persistentOutput(PersistentOStream & os) const { os << theSubProcessHandler << theME << theStats << theDiagrams << isMirror << theNDim << partonDims << theLastDiagramIndex << theExternalDiagram << theMEInfo << theLastRandomNumbers << theMEPartonData << theLastPDFWeight << ounit(theLastCrossSection,nanobarn) << theLastJacobian << theLastME2 << theLastPreweight << ounit(theLastMECrossSection,nanobarn) << theLastMEPDFWeight << theLastMECouplings << theHead << theProjectors << theProjector << theKinematicsGenerated << checkedCuts << passedCuts << theCutWeight << theNeedsReshuffling; } void StandardXComb::persistentInput(PersistentIStream & is, int) { is >> theSubProcessHandler >> theME >> theStats >> theDiagrams >> isMirror >> theNDim >> partonDims >> theLastDiagramIndex >> theExternalDiagram >> theMEInfo >> theLastRandomNumbers >> theMEPartonData >> theLastPDFWeight >> iunit(theLastCrossSection,nanobarn) >> theLastJacobian >> theLastME2 >> theLastPreweight >> iunit(theLastMECrossSection,nanobarn) >> theLastMEPDFWeight >> theLastMECouplings >> theHead >> theProjectors >> theProjector >> theKinematicsGenerated >> checkedCuts >> passedCuts >> theCutWeight >> theNeedsReshuffling; } ClassDescription StandardXComb::initStandardXComb; diff --git a/Handlers/StandardXComb.h b/Handlers/StandardXComb.h --- a/Handlers/StandardXComb.h +++ b/Handlers/StandardXComb.h @@ -1,759 +1,767 @@ // -*- C++ -*- // // StandardXComb.h is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 Leif Lonnblad // Copyright (C) 2009-2017 Simon Platzer // // 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_StandardXComb_H #define ThePEG_StandardXComb_H // This is the declaration of the StandardXComb class. #include "ThePEG/Config/ThePEG.h" #include "SubProcessHandler.fh" #include "ThePEG/PDF/PartonExtractor.fh" #include "ThePEG/PDF/PartonBin.h" #include "ThePEG/PDF/PartonBinInstance.h" #include "ThePEG/Utilities/VSelector.h" #include "ThePEG/Utilities/ClassDescription.h" #include "ThePEG/Utilities/Maths.h" #include "ThePEG/Utilities/XSecStat.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/MatrixElement/MEBase.h" #include "ThePEG/Handlers/XComb.h" #include "ThePEG/Handlers/StandardEventHandler.h" #include "ThePEG/Handlers/SubProcessHandler.fh" #include "StandardXComb.fh" namespace ThePEG { /** * The StandardXComb class inherits from the more general XComb class * which stores all information about the generation of a hard * sub-proces for a given pair of incoming particles, a pair of * extracted partons, etc. This class stores more information related * to thestandard process generation scheme in ThePEG, such as the * PartonExtractor and MEBase object used. It also does some of the * administration of the process generation. * * The main function is dSigDR() which returns the differential cross * section w.r.t. a given vector of random numbers in the interval * ]0,1[. In the initialization this is used to pre-sample the phase * space. In the generation phase it is used to give the cross section * for a phase space point, and if this StandardXComb is chosen the * construct() function is called to generate the actual sub-process. * * @see ParonExtractor * @see MEBase * @see Cuts * @see StdXCombGroup */ class StandardXComb: public XComb { public: /** A vector of DiagramBase objects. */ typedef MEBase::DiagramVector DiagramVector; /** A vector of indices. */ typedef MEBase::DiagramIndex DiagramIndex; /** MEBase needs to be a friend. */ friend class MEBase; public: /** @name Standard constructors and destructors. */ //@{ /** * Standard constructor. */ StandardXComb(Energy newMaxEnergy, const cPDPair & inc, tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler, tPExtrPtr newExtractor, tCascHdlPtr newCKKW, const PBPair & newPartonBins, tCutsPtr newCuts, tMEPtr newME, const DiagramVector & newDiagrams, bool mir, tStdXCombPtr newHead = tStdXCombPtr()); /** * Constructor given a head xcomb. */ StandardXComb(tStdXCombPtr newHead, const PBPair & newPartonBins, tMEPtr newME, const DiagramVector & newDiagrams); /** * Default constructor. */ StandardXComb(); /** * Destructor. */ virtual ~StandardXComb(); /** * Constructor used by MEBase to create a temporary object to store info. */ StandardXComb(tMEPtr me, const tPVector & parts, DiagramIndex i); //@} /** @name Utilities for incoming partons. */ //@{ /** * Properly setup the PartonBinInstance objects provided a sub * process has been constructed using this XComb. */ void recreatePartonBinInstances(Energy2 scale); /** * Fill the variables needed to generate remnants; momenta will be * used from the partons set in this xcomb, but random numbers need * to be provided to (re)generate variables not fixed by the * incoming partons. */ void refillPartonBinInstances(const double* r); /** * Setup information on incoming partons depending * on the information previously supplied through the * choice of diagram and incoming momenta in the first * two entries of meMomenta(). Partons are not actually * extracted from the incoming particles, though a subprocess * detached from the current Event may be created. */ bool setIncomingPartons(tStdXCombPtr labHead = tStdXCombPtr()); /** * Fill phase space information as far as possible */ void fill(const PPair& newParticles, const PPair& newPartons, const vector& newMEMomenta, const DVector& newLastRandomNumbers = DVector()); //@} /** @name Access the assigned objects used in the generation. */ //@{ /** * Return a pointer to the corresponding sub-process handler. May be * null if the standard process generation in ThePEG was not used. */ tcSubHdlPtr subProcessHandler() const { return theSubProcessHandler; } /** * The matrix element to be used. */ tMEPtr matrixElement() const { return theME; } /** * Return a pointer to the head XComb this XComb * depends on. May return NULL, if this is not a * member of a XComb group. */ tStdXCombPtr head() const { return theHead; } /** * Set the head XComb pointer. */ void head(tStdXCombPtr headXC) { theHead = headXC; } /** * Return a selector object of xcombs to choose subprocesses * different than the one currently integrated. */ Selector& projectors() { return theProjectors; } /** * Return a selector object of xcombs to choose subprocesses * different than the one currently integrated. */ const Selector& projectors() const { return theProjectors; } /** * Return a pointer to a projector xcomb which will generate a subprocess * different from the one just integrated. */ tStdXCombPtr lastProjector() const { return theProjector; } /** * Set a pointer to a projector xcomb which will generate a subprocess * different from the one just integrated. */ void lastProjector(tStdXCombPtr pxc) { theProjector = pxc; } //@} /** @name Main functions used for the generation. */ //@{ /** * Try to determine if this subprocess is at all possible. */ virtual bool checkInit(); /** * The number of dimensions of the phase space used to generate this * process. */ virtual int nDim() const { return theNDim; } /** * Return the parton extraction dimensions */ const pair& partonDimensions() const { return partonDims; } /** * Return true, if the current configuration will pass the cuts */ bool willPassCuts(); /** * Return the cut weight encountered from fuzzy cuts */ double cutWeight() const { return theCutWeight; } /** * Reset all saved data about last generated phasespace point; */ virtual void clean(); /** * Return true, if kinematics have already been generated */ bool kinematicsGenerated() const { return theKinematicsGenerated; } /** * Indicate that kinematics have been generated */ void didGenerateKinematics() { theKinematicsGenerated = true; } /** * Generate a phase space point from a vector \a r of \a nr numbers * in the interval ]0,1[ and return the corresponding differential * cross section. */ virtual CrossSection dSigDR(const pair ll, int nr, const double * r); /** * If this XComb has a head XComb, return the cross section * differential in the variables previously supplied. The PDF weight * is taken from the lastPDFWeight supplied by the head XComb * object. */ CrossSection dSigDR(const double * r); /** + * If variations are available for the subprocess handled, generate + * and return a map of optional weights to be included for the + * event; this version defaults to an implementation in MEBase but + * can be overloaded by inheriting XComb objects. + */ + virtual map generateOptionalWeights(); + + /** * Return the PDF weight used in the last call to dSigDR */ double lastPDFWeight() const { return theLastPDFWeight; } /** * Return the cross section calculated in the last call to dSigDR */ CrossSection lastCrossSection() const { return theLastCrossSection; } /** * Check if a reshuffling is required when constructing the hard * subprocess. */ void checkReshufflingNeeds(); /** * Return true if a reshuffling is required when constructing the hard * subprocess. */ bool needsReshuffling() const { return theNeedsReshuffling; } /** * Perform the reshuffling from hardProcessMass to mass values, * given outgoing momenta */ void reshuffle(vector&) const; /** * Construct a sub-process object from the information available. */ virtual tSubProPtr construct(); //@} /** @name Functions used for collecting statistics. */ //@{ /** * The statistics object for this XComb. */ virtual const XSecStat & stats() const { return theStats; } /** * Select the current event. It will later be rejected with a * probability given by \a weight. */ virtual void select(double weight) { theStats.select(weight); } /** * Accept the current event assuming it was previously selcted. */ virtual void accept() { theStats.accept(); } /** * Reweight a selected and accepted event. */ void reweight(double oldWeight, double newWeight) { theStats.reweight(oldWeight,newWeight); } /** * Reject the current event assuming it was previously accepted. If * weighted events are produced, the \a weight should be the same as * the previous call to select(double). */ virtual void reject(double weight = 1.0) { theStats.reject(weight); } /** * Reset statistics. */ virtual void reset() { theStats.reset(); } //@} /** @name Access information used by the MEBase object. */ //@{ /** * The diagrams used by the matrix element. */ const DiagramVector & diagrams() const { return theDiagrams; } /** * True if the TreeDiagram's for this matrix element should in fact * be mirrored before used to create an actual sub-rocess. */ bool mirror() const { return isMirror; } /** * Return the momenta of the partons to be used by the matrix * element object, in the order specified by the TreeDiagram objects * given by the matrix element. */ const vector & meMomenta() const { return theMEMomenta; } /** * Return the last selected diagram. */ tcDiagPtr lastDiagram() const { if ( !theExternalDiagram ) return diagrams()[lastDiagramIndex()]; return theExternalDiagram; } /** * Return the parton types to be used by the matrix element object, * in the order specified by the TreeDiagram objects given by the * matrix element. */ const cPDVector & mePartonData() const { return theMEPartonData; } /** * Return the index of the last selected diagram. */ DiagramIndex lastDiagramIndex() const { return theLastDiagramIndex; } /** * Get information saved by the matrix element in the calculation of * the cross section to be used later when selecting diagrams and * colour flow. */ const DVector & meInfo() const { return theMEInfo; } /** * Set information saved by the matrix element in the calculation of * the cross section to be used later when selecting diagrams and * colour flow. */ void meInfo(const DVector & info) { theMEInfo = info; } /** * Return the random numbers used to generate the * last phase space point, if the matrix element * requested so. */ const DVector& lastRandomNumbers() const { return theLastRandomNumbers; } /** * Get the last jacobian obtained when generating the kinematics * for the call to dSigHatDR. */ double jacobian() const { return theLastJacobian; } /** * Return the matrix element squared as calculated * for the last phase space point. This may optionally * be used by a matrix element for caching. */ double lastME2() const { return theLastME2; } /** * Return the last preweight factor */ double lastPreweight() const { return theLastPreweight; } /** * Return the partonic cross section as calculated * for the last phase space point. This may optionally * be used by a matrix element for caching. */ CrossSection lastMECrossSection() const { return theLastMECrossSection; } /** * Return the PDF weight as calculated * for the last phase space point, if the matrix * element does supply PDF weights. This may optionally * be used by a matrix element for caching. */ double lastMEPDFWeight() const { return theLastMEPDFWeight; } /** * Return the coupling factor as calculated for the lats phase space * point. */ double lastMECouplings() const { return theLastMECouplings; } //@} /** * Construct the corresponding SubProcess object if it hasn't been * done before. */ virtual void newSubProcess(bool group = false); /** * Return the momenta of the partons to be used by the matrix * element object, in the order specified by the TreeDiagram objects * given by the matrix element. */ vector & meMomenta() { return theMEMomenta; } /** * Access the random numbers used to generate the * last phase space point, if the matrix element * requested so. */ DVector& lastRandomNumbers() { return theLastRandomNumbers; } /** * Return the parton types to be used by the matrix element object, * in the order specified by the TreeDiagram objects given by the * matrix element. */ cPDVector & mePartonData() { return theMEPartonData; } /** * Set a diagram to be used instead of the one selected by the matrix * element. */ void externalDiagram(tcDiagPtr diag) { theExternalDiagram = diag; } /** * Set the last selected diagram. */ void lastDiagramIndex(DiagramIndex i) { theLastDiagramIndex = i; } /** * Set the PDF weight used in the last call to dSigDR */ void lastPDFWeight(double w) { theLastPDFWeight = w; } /** * Set the cross section calculated in the last call to dSigDR */ void lastCrossSection(CrossSection s) { theLastCrossSection = s; } /** * Set the last jacobian obtained when generating the kinematics for * the call to dSigHatDR. */ void jacobian(double j) { theLastJacobian = j; } /** * Set the matrix element squared as calculated * for the last phase space point. This may optionally * be used by a matrix element for caching. */ void lastME2(double v) { theLastME2 = v; } /** * Set the last preweight factor */ void lastPreweight(double w) { theLastPreweight = w; } /** * Set the partonic cross section as calculated * for the last phase space point. This may optionally * be used by a matrix element for caching. */ void lastMECrossSection(CrossSection v) { theLastMECrossSection = v; } /** * Set the PDF weight as calculated * for the last phase space point, if the matrix * element does supply PDF weights. This may optionally * be used by a matrix element for caching. */ void lastMEPDFWeight(double v) { theLastMEPDFWeight = v; } /** * Set the coupling factor */ void lastMECouplings(double v) { theLastMECouplings = v; } 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 interface. */ static void Init(); private: /** * The corresponding sub-process handler */ tSubHdlPtr theSubProcessHandler; /** * The matrix element to be used. */ tMEPtr theME; /** * Statistics gathering for this XComb. */ XSecStat theStats; /** * The diagrams used by the matrix element. */ DiagramVector theDiagrams; /** * True if the TreeDiagram's for this matrix element should in fact * be mirrored before used to create an actual sub-rocess. */ bool isMirror; /** * The number of dimensions of the phase space used to generate this * process. */ int theNDim; protected: /** * The number of dimensions of the phase space used for each of the * incoming partons. */ pair partonDims; private: /** * True, if kinematics have already been generated */ bool theKinematicsGenerated; /** * The momenta of the partons to be used by the matrix element * object, in the order specified by the TreeDiagram objects given * by the matrix element. */ vector theMEMomenta; /** * The parton types to be used by the matrix element object, in the * order specified by the TreeDiagram objects given by the matrix * element. */ cPDVector theMEPartonData; /** * A diagram to be used instead of the one selected by the matrix element. */ tcDiagPtr theExternalDiagram; /** * The last selected tree diagram. */ DiagramIndex theLastDiagramIndex; /** * Information saved by the matrix element in the calculation of the * cross section to be used later when selecting diagrams and colour * flow. */ DVector theMEInfo; /** * The random numbers used to generate the * last phase space point, if the matrix element * requested so. */ DVector theLastRandomNumbers; /** * The PDF weight used in the last call to dSigDR */ double theLastPDFWeight; /** * The cross section calculated in the last call to dSigDR */ CrossSection theLastCrossSection; /** * Save the last jacobian obtained when generating the kinematics for * the call to dSigHatDR. */ double theLastJacobian; /** * The matrix element squared as calculated * for the last phase space point. This may optionally * be used by a matrix element for caching. */ double theLastME2; /** * The last preweight factor */ double theLastPreweight; /** * The partonic cross section as calculated * for the last phase space point. This may optionally * be used by a matrix element for caching. */ CrossSection theLastMECrossSection; /** * The PDF weight as calculated * for the last phase space point, if the matrix * element does supply PDF weights. This may optionally * be used by a matrix element for caching. */ double theLastMEPDFWeight; /** * The coupling factor */ double theLastMECouplings; /** * A pointer to the head XComb this XComb * depends on. May return NULL, if this is not a * member of a XComb group. */ tStdXCombPtr theHead; /** * A selector object of xcombs to choose subprocesses * different than the one currently integrated. */ Selector theProjectors; /** * A pointer to a projector xcomb which will generate a subprocess * different from the one just integrated. */ tStdXCombPtr theProjector; /** * True, if cuts have already been checked */ bool checkedCuts; /** * The result of the last call to willPassCuts */ bool passedCuts; /** * The cut weight encountered from fuzzy cuts */ double theCutWeight; /** * True if a reshuffling is required when constructing the hard * subprocess. */ bool theNeedsReshuffling; /** * Calculate the reshuffling equation given the coefficients */ double reshuffleEquation(double, const vector >&, Energy2) const; /** * Solve the reshuffling equation given the coefficients */ double solveReshuffleEquation(const vector >&, Energy2) const; private: /** * Describe a concrete class with persistent data. */ static ClassDescription initStandardXComb; /** * Private and non-existent assignment operator. */ StandardXComb & operator=(const StandardXComb &); }; /** @cond TRAITSPECIALIZATIONS */ /** * This template specialization informs ThePEG about the base class of * StandardXComb. */ template <> struct BaseClassTrait: public ClassTraitsType { /** Typedef of the base class of StandardXComb. */ typedef XComb NthBase; }; /** * This template specialization informs ThePEG about the name of the * StandardXComb class. */ template <> struct ClassTraits: public ClassTraitsBase { /** Return the class name. */ static string className() { return "ThePEG::StandardXComb"; } }; /** @endcond */ } #endif /* ThePEG_StandardXComb_H */ diff --git a/MatrixElement/MEBase.h b/MatrixElement/MEBase.h --- a/MatrixElement/MEBase.h +++ b/MatrixElement/MEBase.h @@ -1,666 +1,675 @@ // -*- C++ -*- // // MEBase.h is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 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_MEBase_H #define ThePEG_MEBase_H // This is the declaration of the MEBase class. #include "ThePEG/Handlers/HandlerBase.h" #include "ThePEG/EventRecord/SubProcess.h" #include "ThePEG/MatrixElement/DiagramBase.h" #include "ThePEG/MatrixElement/ColourLines.h" #include "ThePEG/MatrixElement/Amplitude.h" #include "ThePEG/Handlers/LastXCombInfo.h" #include "ThePEG/Handlers/StandardXComb.fh" #include "ReweightBase.h" #include "ThePEG/Handlers/EventHandler.fh" #include "ThePEG/Handlers/StandardEventHandler.fh" #include "ThePEG/Handlers/SubProcessHandler.fh" #include "ThePEG/PDF/PartonBin.fh" #include "MEBase.fh" namespace ThePEG { /** * The MEBase class is the base class of all objects * representing hard matrix elements in ThePEG. There are three * methods which must be overridden by a concrete subclass:
* * includedDiagrams(tcPDPair) should return a vector of DiagramBase * objects describing the diagrams used for this matrix element for * the given pair of incoming parton types. These DiagramBases are * used to identify the incoming and outgoing partons which can be * handled by the process generation scheme, and is also used to * cnstruct a corresponding SubProcess object. * * scale() should return the scale associated with the phase space * point set with the last call to setKinematics(...) or * generateKinematics(...). * * me() should return the the matrix element squared using the the * type and momentum of the incoming and outgoing partons, previously * set by the setKinematics(...) or generateKinematics(...) member * functions, accessible through the methods meMomenta() and * mePartonData() inherited from LastXCombInfo, and/or from * information stored by sub classes. The returned value should be * dimensionless suitable scaled by the total invariant mass squared * (accessible through the sHat() member function). Any user of this * method must make sure that the setKinematics(...) member function * has been appropriately called before. * * colourGeometries() should return a Selector with the possible * ColourLines objects weighted by their relative probabilities given * the information set by the last call to setKinematics(...) or * generateKinematics(...). * * There are other virtula functions which may be overridden as listed * below. * * @see \ref MEBaseInterfaces "The interfaces" * defined for MEBase. * @see DiagramBase * @see ColourLines * */ class MEBase: public HandlerBase, public LastXCombInfo { public: /** A vector of pointers to DiagramBase objects. */ typedef vector DiagramVector; /** The size_type used in the DiagramVector. */ typedef DiagramVector::size_type DiagramIndex; /** A vector of pointers to ReweightBase objects. */ typedef vector ReweightVector; public: /** @name Standard constructors and destructors. */ //@{ /** * Default constructor. */ MEBase(); /** * Destructor. */ virtual ~MEBase(); //@} public: /** @name Virtual functions to be overridden by sub-classes.. */ //@{ /** * Return the order in \f$\alpha_S\f$ in which this matrix element * is given. */ virtual unsigned int orderInAlphaS() const = 0; /** * Return the order in \f$\alpha_{EM}\f$ in which this matrix * element is given. Returns 0. */ virtual unsigned int orderInAlphaEW() const = 0; /** * Return the matrix element for the kinematical configuation * previously provided by the last call to setKinematics(), suitably * scaled by sHat() to give a dimension-less number. */ virtual double me2() const = 0; /** * Return the scale associated with the phase space point provided * by the last call to setKinematics(). */ virtual Energy2 scale() const = 0; /** * Return the value of \f$\alpha_S\f$ associated with the phase * space point provided by the last call to setKinematics(). This * versions returns SM().alphaS(scale()). */ virtual double alphaS() const; /** * Return the value of \f$\alpha_EM\f$ associated with the phase * space point provided by the last call to setKinematics(). This * versions returns SM().alphaEM(scale()). */ virtual double alphaEM() const; /** * Set the typed and momenta of the incoming and outgoing partons to * be used in subsequent calls to me() and colourGeometries(). */ void setKinematics(tPPair in, const PVector & out); /** * Set the typed and momenta of the incoming and outgoing partons to * be used in subsequent calls to me() and colourGeometries() * according to the associated XComb object. If the function is * overridden in a sub class the new function must call the base * class one first. */ virtual void setKinematics() {} /** * construct the spin information for the interaction */ virtual void constructVertex(tSubProPtr sub); /** * construct the spin information for the interaction */ virtual void constructVertex(tSubProPtr sub, const ColourLines* cl); /** * The number of internal degreed of freedom used in the matrix * element. This default version returns 0; */ virtual int nDim() const; /** * Generate internal degrees of freedom given nDim() uniform random * numbers in the interval ]0,1[. To help the phase space generator, * the 'dSigHatDR' should be a smooth function of these numbers, * although this is not strictly necessary. The return value should * be true of the generation succeeded. If so the generated momenta * should be stored in the meMomenta() vector. */ virtual bool generateKinematics(const double * r) = 0; /** * Return true, if this matrix element expects * the incoming partons in their center-of-mass system */ virtual bool wantCMS() const { return true; } /** * If this is a dependent matrix element in a ME group, return true, * if cuts should be inherited from the head matrix element, i.e. no * cut is being applied to the dependent matrix element if the head * configuration has passed the cuts. */ virtual bool headCuts() const { return false; } /** * If this is a dependent matrix element in a ME group, return true, * if cuts should be ignored. */ virtual bool ignoreCuts() const { return false; } /** * If this is a dependent matrix element in a ME group, return true, * if it applies to the process set in lastXComb() */ virtual bool apply() const { return true; } /** * Return the matrix element squared differential in the variables * given by the last call to generateKinematics(). */ virtual CrossSection dSigHatDR() const = 0; /** + * If variations are available for the subprocess handled, generate + * and return a map of optional weights to be included for the + * event. + */ + virtual map generateOptionalWeights() { + return map(); + } + + /** * Return true, if this matrix element will generate momenta for the * incoming partons itself. The matrix element is required to store * the incoming parton momenta in meMomenta()[0,1]. No mapping in * tau and y is performed by the PartonExtractor object, if a * derived class returns true here. The phase space jacobian is to * include a factor 1/(x1 x2). */ virtual bool haveX1X2() const { return false; } /** * Return true, if this matrix element provides the PDF * weight for the first incoming parton itself. */ virtual bool havePDFWeight1() const { return false; } /** * Return true, if this matrix element provides the PDF * weight for the second incoming parton itself. */ virtual bool havePDFWeight2() const { return false; } /** * Return true, if the XComb steering this matrix element * should keep track of the random numbers used to generate * the last phase space point */ virtual bool keepRandomNumbers() const { return false; } /** * Comlete a SubProcess object using the internal degrees of freedom * generated in the last generateKinematics() (and possible other * degrees of freedom which was intergated over in dSigHatDR(). This * default version does nothing. Will be made purely virtual in the * future. */ virtual void generateSubCollision(SubProcess &); /** * Clear the information previously provided by a call to * setKinematics(...). */ virtual void clearKinematics(); /** * Add all possible diagrams with the add() function. */ virtual void getDiagrams() const = 0; /** * Return true, if this matrix element does not want to * make use of mirroring processes; in this case all * possible partonic subprocesses with a fixed assignment * of incoming particles need to be provided through the diagrams * added with the add(...) method. */ virtual bool noMirror () const { return false; } /** * Return all possible diagrams. */ const DiagramVector & diagrams() const { if ( theDiagrams.empty() ) getDiagrams(); return theDiagrams; } /** * Return a Selector with possible colour geometries for the selected * diagram weighted by their relative probabilities. */ virtual Selector colourGeometries(tcDiagPtr diag) const = 0; /** * Select a ColpurLines geometry. The default version returns a * colour geometry selected among the ones returned from * colourGeometries(tcDiagPtr). */ virtual const ColourLines & selectColourGeometry(tcDiagPtr diag) const; /** * With the information previously supplied with the * setKinematics(...) method, a derived class may optionally * override this method to weight the given diagrams with their * (although certainly not physical) relative probabilities. */ virtual Selector diagrams(const DiagramVector &) const { return Selector(); } /** * Select a diagram. Default version uses diagrams(const * DiagramVector &) to select a diagram according to the * weights. This is the only method used that should be outside of * MEBase. */ virtual DiagramIndex diagram(const DiagramVector &) const; /** * Return true if this matrix element has associated (p)reWeight * objects assigned. */ inline bool reweighted() const { return reweights.size() > 0 || preweights.size() > 0; } /** * With the information previously supplied with the * setKinematics(...) methods, return the combined effects of the * reweighters. */ double reWeight() const; /** * With the information previously supplied with the * setKinematics(...) methods, return the comined effects of the * peweighters. */ double preWeight() const; /** * Add objects to the list of reweighters. */ void addReweighter(tReweightPtr rw); /** * Add objects to the list of preweighters. */ void addPreweighter(tReweightPtr rw); /** * Return the amplitude associated with this matrix element. This * function is allowed to return the null pointer if the amplitude * is not available. */ Ptr::pointer amplitude() const { return theAmplitude; } /** * Set the amplitude associated with this matrix element. */ void amplitude(Ptr::pointer amp) { theAmplitude = amp; } //@} public: /** @name Acces information about the last generated phase space point. */ //@{ /** * Return the last set invariant mass squared. */ Energy2 sHat() const { return lastSHat(); } /** * Return the factor with which this matrix element was last * pre-weighted. */ double preweight() const { return lastPreweight(); } /** * Inform this matrix element that a new phase space * point is about to be generated, so all caches should * be flushed. */ virtual void flushCaches() {} /** * For the given event generation setup return a xcomb object * appropriate to this matrix element. */ virtual StdXCombPtr makeXComb(Energy newMaxEnergy, const cPDPair & inc, tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler, tPExtrPtr newExtractor, tCascHdlPtr newCKKW, const PBPair & newPartonBins, tCutsPtr newCuts, const DiagramVector & newDiagrams, bool mir, const PartonPairVec& allPBins, tStdXCombPtr newHead = tStdXCombPtr(), tMEPtr newME = tMEPtr()); /** * For the given event generation setup return a dependent xcomb object * appropriate to this matrix element. */ virtual StdXCombPtr makeXComb(tStdXCombPtr newHead, const PBPair & newPartonBins, const DiagramVector & newDiagrams, tMEPtr newME = tMEPtr()); /** * Fill the projectors object of xcombs to choose subprocesses * different than the one currently integrated. */ virtual void fillProjectors() { } /** * Set the XComb object to be used in the next call to * generateKinematics() and dSigHatDR(). */ virtual void setXComb(tStdXCombPtr); /** * Retrieve information obtained in the calculation of the cross * section to be used later when selecting diagrams and colour flow. */ const DVector & meInfo() const; /** * Save information obtained in the calculation of the cross * section to be used later when selecting diagrams and colour flow. */ void meInfo(const DVector & info) const; /** * If this matrix element is to be used 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. */ virtual int maxMultCKKW() const { return theMaxMultCKKW; } /** * If this matrix element is to be used 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. */ virtual int minMultCKKW() const { return theMinMultCKKW; } /** * If this matrix element is to be used together with others for * CKKW reweighting and veto, this will set the multiplicity of * outgoing particles in the highest multiplicity matrix element in * the group. */ virtual void maxMultCKKW(int mult) { theMaxMultCKKW = mult; } /** * If this matrix element is to be used together with others for * CKKW reweighting and veto, this will set the multiplicity of * outgoing particles in the lowest multiplicity matrix element in * the group. */ virtual void minMultCKKW(int mult) { theMinMultCKKW = mult; } /** * Set veto scales on the particles at the given * SubProcess which has been generated using this * matrix element. */ virtual void setVetoScales(tSubProPtr) const {} //@} 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: /** * To be used by sub classes in the getDiagrams() method to add * included diagrams. */ void add(DiagPtr dp) const { theDiagrams.push_back(dp); } /** * Access the momenta set by the last call to generateKinematics(). */ vector & meMomenta(); using LastXCombInfo::meMomenta; /** * Set the matrix element squared as calculated * for the last phase space point. This may optionally * be used by a matrix element for caching. */ void lastME2(double v) const; using LastXCombInfo::lastME2; /** * Set the last preweight factor */ void lastPreweight(double w) const; using LastXCombInfo::lastPreweight; /** * Set the partonic cross section as calculated * for the last phase space point. This may optionally * be used by a matrix element for caching. */ void lastMECrossSection(CrossSection v) const; using LastXCombInfo::lastMECrossSection; /** * Set the PDF weight as calculated * for the last phase space point, if the matrix * element does supply PDF weights. This may optionally * be used by a matrix element for caching. */ void lastMEPDFWeight(double v) const; using LastXCombInfo::lastMEPDFWeight; /** * Set the coupling weight as calculated * for the last phase space point */ void lastMECouplings(double v) const; using LastXCombInfo::lastMECouplings; /** * Set the last jacobian obtained when generating the kinematics for * the call to dSigHatDR. */ void jacobian(double j); using LastXCombInfo::jacobian; /** * Initialize all member variables from another * MEBase object. * * @TODO remove? */ void use(tcMEPtr other); /** * Initialize the diagrams from another MEBase object. */ void useDiagrams(tcMEPtr other) 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(); //@} private: /** * The diagrams included for this matrix element. */ mutable DiagramVector theDiagrams; /** * The reweight objects modifying this matrix element. */ ReweightVector reweights; /** * The preweight objects modifying this matrix element. */ ReweightVector preweights; /** * The amplitude associated with this matrix element. */ Ptr::pointer theAmplitude; /** * If this matrix element is to be used 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 matrix element is to be used 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; private: /** * Describe an abstract base class with persistent data. */ static AbstractClassDescription initMEBase; /** * Private and non-existent assignment operator. */ MEBase & operator=(const MEBase &); }; } namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ /** * This template specialization informs ThePEG about the base class of * MEBase. */ template <> struct BaseClassTrait: public ClassTraitsType { /** Typedef of the base class of MEBase. */ typedef HandlerBase NthBase; }; /** * This template specialization informs ThePEG about the name of the * MEBase class. */ template <> struct ClassTraits: public ClassTraitsBase { /** Return the class name. */ static string className() { return "ThePEG::MEBase"; } }; /** @endcond */ } #include "ThePEG/Handlers/StandardXComb.h" #endif /* ThePEG_MEBase_H */