Page MenuHomeHEPForge
Files F7878561

No OneTemporary
Actions

View Options

diff --git a/DipoleShower/Merging2/MFactory.cc b/DipoleShower/Merging2/MFactory.cc
--- a/DipoleShower/Merging2/MFactory.cc
+++ b/DipoleShower/Merging2/MFactory.cc
@@ -1,1130 +1,730 @@
- // -*- C++ -*-
- //
- // MergeboxFactory.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
- // Copyright (C) 2002-2012 The Herwig Collaboration
- //
- // Herwig is licenced under version 2 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 MergeboxFactory class.
- //
- #include "MFactory.h"
- #include "Node.h"
+ // -*- C++ -*-
+ //
+ // MergeboxFactory.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
+ // Copyright (C) 2002-2012 The Herwig Collaboration
+ //
+ // Herwig is licenced under version 2 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 MergeboxFactory class.
+ //
+#include "MFactory.h"
+#include "Node.h"
- #include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
- #include "ThePEG/Interface/ClassDocumentation.h"
- #include "ThePEG/Utilities/DescribeClass.h"
- #include "ThePEG/Interface/Reference.h"
- #include "ThePEG/Interface/RefVector.h"
- #include "ThePEG/Interface/Switch.h"
- #include "ThePEG/Interface/Parameter.h"
- #include "ThePEG/Interface/Command.h"
- #include "ThePEG/Utilities/StringUtils.h"
- #include "ThePEG/Repository/Repository.h"
- #include "ThePEG/Repository/EventGenerator.h"
- #include "Herwig/MatrixElement/Matchbox/Utility/DiagramDrawer.h"
- #include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
+#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Utilities/DescribeClass.h"
+#include "ThePEG/Interface/Reference.h"
+#include "ThePEG/Interface/RefVector.h"
+#include "ThePEG/Interface/Switch.h"
+#include "ThePEG/Interface/Parameter.h"
+#include "ThePEG/Interface/Command.h"
+#include "ThePEG/Utilities/StringUtils.h"
+#include "ThePEG/Repository/Repository.h"
+#include "ThePEG/Repository/EventGenerator.h"
+#include "Herwig/MatrixElement/Matchbox/Utility/DiagramDrawer.h"
+#include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
- #include "ThePEG/MatrixElement/MEBase.h"
+#include "ThePEG/MatrixElement/MEBase.h"
- #include "ThePEG/Persistency/PersistentOStream.h"
- #include "ThePEG/Persistency/PersistentIStream.h"
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
- #include "Herwig/MatrixElement/Matchbox/Base/DipoleRepository.h"
- #include "Herwig/MatrixElement/Matchbox/Utility/SU2Helper.h"
+#include "Herwig/MatrixElement/Matchbox/Base/DipoleRepository.h"
+#include "Herwig/MatrixElement/Matchbox/Utility/SU2Helper.h"
- #include "ThePEG/Cuts/JetFinder.h"
- #include "fastjet/ClusterSequence.hh"
+#include "ThePEG/Cuts/JetFinder.h"
+#include "fastjet/ClusterSequence.hh"
- #include <boost/progress.hpp>
+#include <boost/progress.hpp>
- #include <iterator>
- using std::ostream_iterator;
+#include <iterator>
+using std::ostream_iterator;
- using namespace Herwig;
- using std::ostream_iterator;
+using namespace Herwig;
+using std::ostream_iterator;
- MFactory::MFactory() :MatchboxFactory(),
- theonlyNLOParts(false),
- theonlyvirtualNLOParts(false),
- theonlyrealNLOParts(false),
- theunitarizeNLOParts(true),
- calc_born(true),
- calc_virtual(true),
- calc_real(true),
- calc_projected_virtual(true),
- calc_inclusive_real(true),
- calc_double_inclusive(true),
- needFOH(false),
- unitarized(true),
- NLOunitarized(true),
- theonlyUnlopsweights(false),
- theonlyk(-1),
- theonlymulti(-1),
- divideSub(-1),
- divideSubNumber(-1),
- theonlyabove(0),
- theonlysub(-1),
- theChooseHistory(3),
- theN(1),
- theM(-1),
- theMergePTsmearing(0.),
- theIRSafeRATIO(100),
- theStairFactor(0.),
- theMergePT(2.*GeV),
- theNLOMergePT(2.*GeV),
- theIRSafePT(1000000.0 * GeV),ransetup(false),defMERegionByJetAlg(false),
- theGamma(1.){
+MFactory::MFactory() :MatchboxFactory(),
+calc_born(true), calc_virtual(true),
+calc_real(true), theonlyUnlopsweights(false),
+theonlyk(-1), theonlymulti(-1),
+divideSub(-1), divideSubNumber(-1),
+theonlyabove(0), theonlysub(-1),
+ransetup(false){}
+
+MFactory::~MFactory(){}
+
+IBPtr MFactory::clone() const {
+ return new_ptr(*this);
+}
+
+IBPtr MFactory::fullclone() const {
+ return new_ptr(*this);
+}
+
+
+
+
+void MFactory::fill_amplitudes() {
+
+ olpProcesses().clear();
+ if ( particleGroups().find("j") == particleGroups().end() ) throw InitException() << "Could not find a jet particle group named 'j'";
+
+ // rebind the particle data objects
+ for ( map<string, PDVector>::iterator g = particleGroups().begin() ; g != particleGroups().end() ; ++g ) {
+ for ( PDVector::iterator p = g->second.begin() ; p != g->second.end() ; ++p ) {
+ *p = getParticleData((**p).id());
+ }
}
+
+ const PDVector& partons = particleGroups()["j"];
+
+ unsigned int nl = 0;
+
+ for ( PDVector::const_iterator p = partons.begin() ; p != partons.end() ; ++p )
+ if ( abs((**p).id()) < 6 ) ++nl;
+
+ nLight(nl / 2);
+ processMap[0] = getProcesses()[0];
+ if(MH()->M()>=0)
+ setHighestVirt(processMap[0].size()+MH()->M());
+
+ for ( int i = 1 ; i <= MH()->N() ; ++i ) {
+ processMap[i] = processMap[i - 1];
+ processMap[i].push_back("j");
+ }
+
+ for ( int i = 0 ; i <= MH()->N() ; ++i ) {
+ vector<Ptr<MatchboxMEBase>::ptr> ames = makeMEs(processMap[i], orderInAlphaS() + i,i<MH()->M()+1);
+ copy(ames.begin(), ames.end(), back_inserter(pureMEsMap()[i]));
+ }
+}
- MFactory::~MFactory() {
- }
-
- IBPtr MFactory::clone() const {
- return new_ptr(*this);
- }
-
- IBPtr MFactory::fullclone() const {
- return new_ptr(*this);
- }
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- void MFactory::fill_amplitudes() {
+void MFactory::prepare_BV(int i) {
+ // check if we have virtual contributions
+ bool haveVirtuals = true;
+ bool virtualsAreDR = false;
+ bool virtualsAreDRbar = false;
+ bool virtualsAreCDR = false;
+ bool virtualsAreCS = false;
+ bool virtualsAreBDK = false;
+ bool virtualsAreExpanded = false;
+ for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator born = pureMEsMap()[i].begin() ; born != pureMEsMap()[i].end() ; ++born ) {
- olpProcesses().clear();
+ prepareME(*born);
- if ( particleGroups().find("j") == particleGroups().end() ) throw InitException() << "Could not find a jet particle group named 'j'";
- // rebind the particle data objects
- for ( map<string, PDVector>::iterator g = particleGroups().begin() ; g != particleGroups().end() ; ++g ) {
- for ( PDVector::iterator p = g->second.begin() ; p != g->second.end() ; ++p ) {
- #ifndef NDEBUG
- long checkid = (**p).id();
- #endif
- *p = getParticleData((**p).id());
- assert((**p).id() == checkid);
+ if ( (*born)->isOLPTree() ) {
+ int id = orderOLPProcess((*born)->subProcess(),
+ (*born)->matchboxAmplitude(),
+ ProcessType::treeME2);
+ (*born)->olpProcess(ProcessType::treeME2,id);
+ id = orderOLPProcess((*born)->subProcess(),
+ (*born)->matchboxAmplitude(),
+ ProcessType::colourCorrelatedME2);
+ (*born)->olpProcess(ProcessType::colourCorrelatedME2,id);
+
+ bool haveGluon = false;
+ for ( PDVector::const_iterator p = (*born)->subProcess().legs.begin();
+ p != (*born)->subProcess().legs.end(); ++p )
+ if ( (**p).id() == 21 ) {
+ haveGluon = true;
+ break;
+ }
+ if ( haveGluon ) {
+ id = orderOLPProcess((*born)->subProcess(),
+ (*born)->matchboxAmplitude(),
+ ProcessType::spinColourCorrelatedME2);
+ (*born)->olpProcess(ProcessType::spinColourCorrelatedME2,id);
}
}
-
- const PDVector& partons = particleGroups()["j"];
-
- unsigned int nl = 0;
-
- for ( PDVector::const_iterator p = partons.begin() ; p != partons.end() ; ++p )
- if ( abs((**p).id()) < 6 ) ++nl;
-
- nLight(nl / 2);
- processMap[0] = getProcesses()[0];
- if(theM>=0&&!theonlyrealNLOParts)
- setHighestVirt(processMap[0].size()+theM);
-
- for ( int i = 1 ; i <= theN ; ++i ) {
- processMap[i] = processMap[i - 1];
- if (particleGroups().find("jm") != particleGroups().end()){
- for(vector<string>::iterator it= processMap[i].begin();it!= processMap[i].end();it++){
- if(*it=="j")*it="jm";
- if(*it=="p")*it="pm";
- }
- processMap[i].push_back("jm");
- }else{
- processMap[i].push_back("j");
- }
-
- }
-
- for ( int i = 0 ; i <= theN ; ++i ) {
- vector<Ptr<MatchboxMEBase>::ptr> ames = makeMEs(processMap[i], orderInAlphaS() + i,i<theM+1);
- copy(ames.begin(), ames.end(), back_inserter(pureMEsMap()[i]));
- }
- }
-
- void MFactory::prepare_BV(int i) {
- // check if we have virtual contributions
- bool haveVirtuals = true;
- bool virtualsAreDR = false;
- bool virtualsAreDRbar = false;
- bool virtualsAreCDR = false;
- bool virtualsAreCS = false;
- bool virtualsAreBDK = false;
- bool virtualsAreExpanded = false;
- for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator born = pureMEsMap()[i].begin() ; born != pureMEsMap()[i].end() ; ++born ) {
-
- prepareME(*born);
-
-
- if ( (*born)->isOLPTree() ) {
- int id = orderOLPProcess((*born)->subProcess(),
- (*born)->matchboxAmplitude(),
- ProcessType::treeME2);
- (*born)->olpProcess(ProcessType::treeME2,id);
+ if ( (*born)->isOLPLoop() && i <= MH()->M()) {
+ int id = orderOLPProcess((*born)->subProcess(),
+ (*born)->matchboxAmplitude(),
+ ProcessType::oneLoopInterference);
+ (*born)->olpProcess(ProcessType::oneLoopInterference,id);
+ if ( !(*born)->onlyOneLoop() && (*born)->needsOLPCorrelators() ) {
id = orderOLPProcess((*born)->subProcess(),
(*born)->matchboxAmplitude(),
ProcessType::colourCorrelatedME2);
(*born)->olpProcess(ProcessType::colourCorrelatedME2,id);
-
- bool haveGluon = false;
- for ( PDVector::const_iterator p = (*born)->subProcess().legs.begin();
- p != (*born)->subProcess().legs.end(); ++p )
- if ( (**p).id() == 21 ) {
- haveGluon = true;
- break;
- }
- if ( haveGluon ) {
- id = orderOLPProcess((*born)->subProcess(),
- (*born)->matchboxAmplitude(),
- ProcessType::spinColourCorrelatedME2);
- (*born)->olpProcess(ProcessType::spinColourCorrelatedME2,id);
- }
- }
- if ( (*born)->isOLPLoop() && i <= theM) {
- int id = orderOLPProcess((*born)->subProcess(),
- (*born)->matchboxAmplitude(),
- ProcessType::oneLoopInterference);
- (*born)->olpProcess(ProcessType::oneLoopInterference,id);
- if ( !(*born)->onlyOneLoop() && (*born)->needsOLPCorrelators() ) {
- id = orderOLPProcess((*born)->subProcess(),
- (*born)->matchboxAmplitude(),
- ProcessType::colourCorrelatedME2);
- (*born)->olpProcess(ProcessType::colourCorrelatedME2,id);
- }
- }
-
-
-
-
-
-
- haveVirtuals &= (**born).haveOneLoop();
- if ( (**born).haveOneLoop() ) {
- virtualsAreDR |= (**born).isDR();
- virtualsAreCDR |= !(**born).isDR();
- virtualsAreDRbar |= (**born).isDRbar();
- virtualsAreCS |= (**born).isCS();
- virtualsAreBDK |= (**born).isBDK();
- virtualsAreExpanded |= (**born).isExpanded();
}
}
- // check the additional insertion operators
- if ( !theVirtualsMap[i].empty() ) haveVirtuals = true;
- for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt = theVirtualsMap[i].begin() ; virt != theVirtualsMap[i].end() ; ++virt ) {
- (**virt).factory(this);
- virtualsAreDR |= (**virt).isDR();
- virtualsAreDRbar |= (**virt).isDRbar();
- virtualsAreCDR |= !(**virt).isDR();
- virtualsAreCS |= (**virt).isCS();
- virtualsAreBDK |= (**virt).isBDK();
- virtualsAreExpanded |= (**virt).isExpanded();
- }
- // check for consistent conventions on virtuals, if we are to include them
- if ( i <= theM ) {
- if ( virtualsAreDR && virtualsAreCDR ) {
- throw InitException() << "Virtual corrections use inconsistent regularization schemes.\n";
- }
- if ( (virtualsAreCS && virtualsAreBDK) || (virtualsAreCS && virtualsAreExpanded) || (virtualsAreBDK && virtualsAreExpanded)
- || (!virtualsAreCS && !virtualsAreBDK && !virtualsAreExpanded) ) {
- throw InitException() << "Virtual corrections use inconsistent conventions on finite terms.\n";
- }
- if ( !haveVirtuals ) {
- throw InitException() << "Could not find amplitudes for all virtual contributions needed.\n";
- }
- }
- // prepare dipole insertion operators
- // Need them for alpha-improvement
- //if ( i <= theM ) {
- for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt
- = DipoleRepository::insertionIOperators(dipoleSet()).begin();
- virt != DipoleRepository::insertionIOperators(dipoleSet()).end(); ++virt ) {
- (**virt).factory(this);
- if ( virtualsAreDRbar )
- (**virt).useDRbar();
- if ( virtualsAreDR )
- (**virt).useDR();
- else
- (**virt).useCDR();
- if ( virtualsAreCS )
- (**virt).useCS();
- if ( virtualsAreBDK )
- (**virt).useBDK();
- if ( virtualsAreExpanded )
- (**virt).useExpanded();
- }
- for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt
- = DipoleRepository::insertionPKOperators(dipoleSet()).begin();
- virt != DipoleRepository::insertionPKOperators(dipoleSet()).end(); ++virt ) {
-
- (**virt).factory(this);
- if ( virtualsAreDRbar )
- (**virt).useDRbar();
- if ( virtualsAreDR )
- (**virt).useDR();
- else
- (**virt).useCDR();
- if ( virtualsAreCS )
- (**virt).useCS();
- if ( virtualsAreBDK )
- (**virt).useBDK();
- if ( virtualsAreExpanded )
- (**virt).useExpanded();
- }
-
- //}
- thevirtualsAreExpandedMap[i] = virtualsAreExpanded;
- }
-
- void MFactory::prepare_R(int i) {
- for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator real = thePureMEsMap[i].begin() ; real != thePureMEsMap[i].end() ; ++real ) {
- prepareME(*real);
+
+ haveVirtuals &= (**born).haveOneLoop();
+ if ( (**born).haveOneLoop() ) {
+ virtualsAreDR |= (**born).isDR();
+ virtualsAreCDR |= !(**born).isDR();
+ virtualsAreDRbar |= (**born).isDRbar();
+ virtualsAreCS |= (**born).isCS();
+ virtualsAreBDK |= (**born).isBDK();
+ virtualsAreExpanded |= (**born).isExpanded();
}
}
+
+ // check the additional insertion operators
+ if ( !theVirtualsMap[i].empty() ) haveVirtuals = true;
+ for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt = theVirtualsMap[i].begin() ; virt != theVirtualsMap[i].end() ; ++virt ) {
+ (**virt).factory(this);
+ virtualsAreDR |= (**virt).isDR();
+ virtualsAreDRbar |= (**virt).isDRbar();
+ virtualsAreCDR |= !(**virt).isDR();
+ virtualsAreCS |= (**virt).isCS();
+ virtualsAreBDK |= (**virt).isBDK();
+ virtualsAreExpanded |= (**virt).isExpanded();
+ }
+
+
+ // check for consistent conventions on virtuals, if we are to include MH()->M()
+ if ( i <= MH()->M() ) {
+ if ( virtualsAreDR && virtualsAreCDR ) {
+ throw InitException() << "Virtual corrections use inconsistent regularization schemes.\n";
+ }
+ if ( (virtualsAreCS && virtualsAreBDK) || (virtualsAreCS && virtualsAreExpanded) || (virtualsAreBDK && virtualsAreExpanded)
+ || (!virtualsAreCS && !virtualsAreBDK && !virtualsAreExpanded) ) {
+ throw InitException() << "Virtual corrections use inconsistent conventions on finite terms.\n";
+ }
+ if ( !haveVirtuals ) {
+ throw InitException() << "Could not find amplitudes for all virtual contributions needed.\n";
+ }
+ }
+
+ // prepare dipole insertion operators
+ // Need MH()->M() for alpha-improvement
+ //if ( i <= MH()->M() ) {
+ for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt
+ = DipoleRepository::insertionIOperators(dipoleSet()).begin();
+ virt != DipoleRepository::insertionIOperators(dipoleSet()).end(); ++virt ) {
+ (**virt).factory(this);
+ if ( virtualsAreDRbar )
+ (**virt).useDRbar();
+ if ( virtualsAreDR )
+ (**virt).useDR();
+ else
+ (**virt).useCDR();
+ if ( virtualsAreCS )
+ (**virt).useCS();
+ if ( virtualsAreBDK )
+ (**virt).useBDK();
+ if ( virtualsAreExpanded )
+ (**virt).useExpanded();
+ }
+ for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt
+ = DipoleRepository::insertionPKOperators(dipoleSet()).begin();
+ virt != DipoleRepository::insertionPKOperators(dipoleSet()).end(); ++virt ) {
+
+ (**virt).factory(this);
+ if ( virtualsAreDRbar )
+ (**virt).useDRbar();
+ if ( virtualsAreDR )
+ (**virt).useDR();
+ else
+ (**virt).useCDR();
+ if ( virtualsAreCS )
+ (**virt).useCS();
+ if ( virtualsAreBDK )
+ (**virt).useBDK();
+ if ( virtualsAreExpanded )
+ (**virt).useExpanded();
+ }
+
+ //}
+ thevirtualsAreExpandedMap[i] = virtualsAreExpanded;
+
+}
- void MFactory::pushB(Ptr<MatchboxMEBase>::ptr born, int i, bool projected) {
- Ptr<MatchboxMEBase>::ptr bornme = (*born).cloneMe();
- bornme->maxMultCKKW(1);
- bornme->minMultCKKW(0);
+void MFactory::prepare_R(int i) {
+ for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator
+ real = thePureMEsMap[i].begin() ;
+ real!= thePureMEsMap[i].end() ; ++real )
+ prepareME(*real);
+}
-
- string pname = fullName() + "/" + bornme->name() + ".Born";
- if ( projected ) pname += "pro";
- if ( !(generator()->preinitRegister(bornme, pname)) ) throw InitException() << "Born ME " << pname << " already existing.";
-
-
-
- bornme->virtuals().clear();
- for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt = theVirtualsMap[i].begin() ; virt != theVirtualsMap[i].end() ; ++virt ) {
- if ( (**virt).apply((*born).diagrams().front()->partons()) ) bornme->virtuals().push_back(*virt);
- }
- for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt = DipoleRepository::insertionIOperators(dipoleSet()).begin() ;
- virt != DipoleRepository::insertionIOperators(dipoleSet()).end() ; ++virt ) {
- if ( (**virt).apply((*born).diagrams().front()->partons()) ) bornme->virtuals().push_back(*virt);
- }
- for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt = DipoleRepository::insertionPKOperators(dipoleSet()).begin() ;
- virt != DipoleRepository::insertionPKOperators(dipoleSet()).end() ; ++virt ) {
- if ( (**virt).apply((*born).diagrams().front()->partons()) ) bornme->virtuals().push_back(*virt);
- }
-
-
-
-
-
-
-
-
-
- Ptr<Node>::ptr clusternode = new_ptr(Node(bornme, i, 0, needFOH));
- clusternode->mergePt(theMergePT);
- clusternode->centralMergePt(theMergePT);
- clusternode->N(theN + getProcesses()[0].size());clusternode->N0( getProcesses()[0].size());
- clusternode->M(theM + getProcesses()[0].size());
- if(theonlyk!=-1)clusternode->onlyN(theonlyk + getProcesses()[0].size());
- clusternode->unitarized(unitarized);
- clusternode->NLOunitarized(NLOunitarized);
- clusternode->calculateInNode(true);
- clusternode->treefactory(this);
-
- theMergingHelper->firstNodeMap(bornme,clusternode);
-
- bornme->factory(this);
- bornme->merger(theMergingHelper);
- if ( projected ) bornme->projectorStage(1);
- else bornme->projectorStage(0);
- clusternode->deepHead(clusternode);
- clusternode->chooseHistory(theChooseHistory);
-
- vector<Ptr<Node>::ptr> temp;
- vector<Ptr<Node>::ptr> temp1;
- temp.push_back(clusternode);
- unsigned int k = 1;
- while (thePureMEsMap[i - k].size() != 0) {
- for ( unsigned int j = 0 ; j < temp.size() ; j++ ) {
- //cout<<"\ntemp[j]->children().size() "<<temp[j]->children().size();
- temp[j]->birth(thePureMEsMap[i - k]);
- for ( unsigned int m = 0 ; m < temp[j]->children().size() ; m++ ) {
- //int numofsplit=temp[j]->children()[m]->nodeME()->numberOfSplittings(DipoleRepository::dipoles(dipoleSet()),thePureMEsMap[i - k+1]);
- //cout<<"\nnumofsplit "<<numofsplit<<temp[j]->children()[m]->nodeME()->name()<<" k "<<k<<" j "<<j<<" m "<<m<<" "<<" "<<i;
- //temp[j]->children()[m]->numberOfSplittings(numofsplit);
-
- if ( i <= theM ) {
+void MFactory::pushB(Ptr<MatchboxMEBase>::ptr born, int i) {
+ Ptr<MatchboxMEBase>::ptr bornme = (*born).cloneMe();
+ bornme->maxMultCKKW(1);
+ bornme->minMultCKKW(0);
+
+
+ string pname = fullName() + "/" + bornme->name() + ".Born";
+ if ( !(generator()->preinitRegister(bornme, pname)) ) throw InitException() << "Born ME " << pname << " already existing.";
+
+ bornme->virtuals().clear();
+ for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt = theVirtualsMap[i].begin() ; virt != theVirtualsMap[i].end() ; ++virt ) {
+ if ( (**virt).apply((*born).diagrams().front()->partons()) ) bornme->virtuals().push_back(*virt);
+ }
+ for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt = DipoleRepository::insertionIOperators(dipoleSet()).begin() ;
+ virt != DipoleRepository::insertionIOperators(dipoleSet()).end() ; ++virt ) {
+ if ( (**virt).apply((*born).diagrams().front()->partons()) ) bornme->virtuals().push_back(*virt);
+ }
+ for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt = DipoleRepository::insertionPKOperators(dipoleSet()).begin() ;
+ virt != DipoleRepository::insertionPKOperators(dipoleSet()).end() ; ++virt ) {
+ if ( (**virt).apply((*born).diagrams().front()->partons()) ) bornme->virtuals().push_back(*virt);
+ }
+
+
+ Ptr<Node>::ptr clusternode = new_ptr(Node(bornme, i, 0));
+ theMergingHelper->firstNodeMap(bornme,clusternode);
+ bornme->factory(this);
+ bornme->merger(theMergingHelper);
+
+ vector<Ptr<Node>::ptr> temp;
+ vector<Ptr<Node>::ptr> temp1;
+ temp.push_back(clusternode);
+ unsigned int k = 1;
+ while (thePureMEsMap[i - k].size() != 0) {
+ for ( unsigned int j = 0 ; j < temp.size() ; j++ ) {
+ temp[j]->birth(thePureMEsMap[i - k]);
+ for ( unsigned int m = 0 ; m < temp[j]->children().size() ; m++ ) {
+ if ( i <= MH()->M() ) {
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt = DipoleRepository::insertionIOperators(dipoleSet()).begin() ;
virt != DipoleRepository::insertionIOperators(dipoleSet()).end() ; ++virt ) {
if ( (**virt).apply((*(temp[j]->children()[m]->nodeME())).diagrams().front()->partons()) ){
- cout<<"\nMatchboxInsertionOperator"<<flush;
Ptr<MatchboxInsertionOperator>::ptr myIOP = (**virt).cloneMe();
- //ostringstream pname;
- //pname << temp[j]->children()[m]->nodeME()->fullName() << "/" << (**virt).name();
- //if ( ! (generator()->preinitRegister(myIOP,pname.str()) ) )
- // throw Exception() << "MatchboxMEBase::cloneDependencies(): Insertion operator " << pname.str() << " already existing." << Exception::runerror;
temp[j]->children()[m]->nodeME()->virtuals().push_back(myIOP);
}
}
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt = DipoleRepository::insertionPKOperators(dipoleSet()).begin() ;
virt != DipoleRepository::insertionPKOperators(dipoleSet()).end() ; ++virt ) {
if ( (**virt).apply((*(temp[j]->children()[m]->nodeME())).diagrams().front()->partons()) ){
Ptr<MatchboxInsertionOperator>::ptr myIOP = (**virt).cloneMe();
- //ostringstream pname;
- //pname << temp[j]->children()[m]->nodeME()->fullName() << "/" << (**virt).name();
- //if ( ! (generator()->preinitRegister(myIOP,pname.str()) ) )
- // throw Exception() << "MatchboxMEBase::cloneDependencies(): Insertion operator " << pname.str() << " already existing." << Exception::runerror;
temp[j]->children()[m]->nodeME()->virtuals().push_back(myIOP);
}
}
-
-
temp[j]->children()[m]->nodeME()->noOneLoop();
- }
- temp1.push_back(temp[j]->children()[m]);
}
+ temp1.push_back(temp[j]->children()[m]);
}
- temp = temp1;
- temp1.clear();
- k++;
}
-
- if(theN>i)
- bornme->needsCorrelations();
- else bornme->needsNoCorrelations();
-
-
- bornme->cloneDependencies();
- if ( theonlyk == -1 || (theonlyk == i-1 &&unitarized) || theonlyk == i) {
- MEs().push_back(bornme);
+ temp = temp1;
+ temp1.clear();
+ k++;
+ }
+
+ if(MH()->N()>i)
+ bornme->needsCorrelations();
+ else bornme->needsNoCorrelations();
+
+ bornme->cloneDependencies();
+ if ( theonlyk == -1 || (theonlyk == i-1 &&unitarized) || theonlyk == i) {
+ MEs().push_back(bornme);
+ }
+}
+
+
+
+
+
+void MFactory::pushV(Ptr<MatchboxMEBase>::ptr born, int i) {
+
+ Ptr<MatchboxMEBase>::ptr nlo = (*born).cloneMe();
+
+ string pname = fullName() + "/" + nlo->name() + ".Virtual";
+ if ( !(generator()->preinitRegister(nlo, pname)) ) throw InitException() << "Virtual ME " << pname << " already existing.";
+ ////////////////////////////////////NLO///////////////////////////
+ nlo->virtuals().clear();
+ if ( !nlo->onlyOneLoop() ) {
+ for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt = theVirtualsMap[i].begin() ; virt != theVirtualsMap[i].end() ; ++virt ) {
+ if ( (**virt).apply((*born).diagrams().front()->partons()) ) nlo->virtuals().push_back(*virt);
+ }
+ for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt = DipoleRepository::insertionIOperators(dipoleSet()).begin() ;
+ virt != DipoleRepository::insertionIOperators(dipoleSet()).end() ; ++virt ) {
+ if ( (**virt).apply((*born).diagrams().front()->partons()) ) nlo->virtuals().push_back(*virt);
+ }
+ for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt = DipoleRepository::insertionPKOperators(dipoleSet()).begin() ;
+ virt != DipoleRepository::insertionPKOperators(dipoleSet()).end() ; ++virt ) {
+ if ( (**virt).apply((*born).diagrams().front()->partons()) ) nlo->virtuals().push_back(*virt);
+ }
+ if ( nlo->virtuals().empty() ) throw InitException() << "No insertion operators have been found for " << (*born).name() << "\n";
+ if ( checkPoles() ) {
+ if ( !thevirtualsAreExpandedMap[i] ) {
+ throw InitException() << "Cannot check epsilon poles if virtuals are not in `expanded' convention.\n";
+ }
}
}
-
-
-
-
-
- void MFactory::pushV(Ptr<MatchboxMEBase>::ptr born, int i, bool projected) {
-
- Ptr<MatchboxMEBase>::ptr nlo = (*born).cloneMe();
- nlo->maxMultCKKW(1);
- nlo->minMultCKKW(0);
- int pro = projected ? 1 : 0;
-
- string pname = fullName() + "/" + nlo->name() + ".Virtual";
- if ( projected ) pname += "pro";
- if ( !(generator()->preinitRegister(nlo, pname)) ) throw InitException() << "Virtual ME " << pname << " already existing.";
- ////////////////////////////////////NLO///////////////////////////
- nlo->virtuals().clear();
- if ( !nlo->onlyOneLoop() ) {
- for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt = theVirtualsMap[i].begin() ; virt != theVirtualsMap[i].end() ; ++virt ) {
- if ( (**virt).apply((*born).diagrams().front()->partons()) ) nlo->virtuals().push_back(*virt);
- }
- for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt = DipoleRepository::insertionIOperators(dipoleSet()).begin() ;
- virt != DipoleRepository::insertionIOperators(dipoleSet()).end() ; ++virt ) {
- if ( (**virt).apply((*born).diagrams().front()->partons()) ) nlo->virtuals().push_back(*virt);
- }
- for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt = DipoleRepository::insertionPKOperators(dipoleSet()).begin() ;
- virt != DipoleRepository::insertionPKOperators(dipoleSet()).end() ; ++virt ) {
- if ( (**virt).apply((*born).diagrams().front()->partons()) ) nlo->virtuals().push_back(*virt);
- }
- if ( nlo->virtuals().empty() ) throw InitException() << "No insertion operators have been found for " << (*born).name() << "\n";
- if ( checkPoles() ) {
- if ( !thevirtualsAreExpandedMap[i] ) {
- throw InitException() << "Cannot check epsilon poles if virtuals are not in `expanded' convention.\n";
- }
+ nlo->doOneLoopNoBorn();
+ ////////////////////////////////////NLO///////////////////////////
+ Ptr<Node>::ptr clusternode = new_ptr(Node(nlo, i, 0));
+ clusternode->virtualContribution(true);
+ theMergingHelper->firstNodeMap(nlo,clusternode);
+ nlo->merger(theMergingHelper);
+
+ vector<Ptr<Node>::ptr> temp;
+ vector<Ptr<Node>::ptr> temp1;
+ temp.push_back(clusternode);
+ unsigned int k = 1;
+ while (thePureMEsMap[i - k].size() != 0) {
+ for ( unsigned int j = 0 ; j < temp.size() ; j++ ) {
+ temp[j]->birth(thePureMEsMap[i - k]);
+ for ( unsigned int m = 0 ; m < temp[j]->children().size() ; ++m ) {
+ temp1.push_back(temp[j]->children()[m]);
}
}
- nlo->doOneLoopNoBorn();
- ////////////////////////////////////NLO///////////////////////////
- Ptr<Node>::ptr clusternode = new_ptr(Node(nlo, i, 0, needFOH));
- clusternode->mergePt(theMergePT);
- clusternode->centralMergePt(theMergePT);
- clusternode->unitarized(unitarized);
- clusternode->NLOunitarized(NLOunitarized);
- clusternode->N(theN + getProcesses()[0].size());clusternode->N0( getProcesses()[0].size());
- clusternode->M(theM + getProcesses()[0].size());
- if(theonlyk!=-1)clusternode->onlyN(theonlyk + getProcesses()[0].size());
- clusternode->calculateInNode(true);
- clusternode->virtualContribution(true);
- clusternode->treefactory(this);
- theMergingHelper->firstNodeMap(nlo,clusternode);
-
- nlo->merger(theMergingHelper);
- //nlo->firstNode(clusternode);
-
- if ( projected ) nlo->projectorStage(1);
- else nlo->projectorStage(0);
- clusternode->deepHead(clusternode);
- clusternode->chooseHistory(theChooseHistory);
- vector<Ptr<Node>::ptr> temp;
- vector<Ptr<Node>::ptr> temp1;
- temp.push_back(clusternode);
- unsigned int k = 1;
- while (thePureMEsMap[i - k].size() != 0) {
- for ( unsigned int j = 0 ; j < temp.size() ; j++ ) {
- temp[j]->birth(thePureMEsMap[i - k]);
- for ( unsigned int m = 0 ; m < temp[j]->children().size() ; ++m ) {
- temp1.push_back(temp[j]->children()[m]);
- }
- }
- temp = temp1;
- k++;
+ temp = temp1;
+ k++;
+ }
+
+ if ( nlo->isOLPLoop() ) {
+ int id = orderOLPProcess(nlo->subProcess(),
+ (*born).matchboxAmplitude(),
+ ProcessType::oneLoopInterference);
+ nlo->olpProcess(ProcessType::oneLoopInterference,id);
+ if ( !nlo->onlyOneLoop() && nlo->needsOLPCorrelators() ) {
+ id = orderOLPProcess(nlo->subProcess(),
+ (*born).matchboxAmplitude(),
+ ProcessType::colourCorrelatedME2);
+ nlo->olpProcess(ProcessType::colourCorrelatedME2,id);
}
-
- if ( nlo->isOLPLoop() ) {
- int id = orderOLPProcess(nlo->subProcess(),
- (*born).matchboxAmplitude(),
- ProcessType::oneLoopInterference);
- nlo->olpProcess(ProcessType::oneLoopInterference,id);
- if ( !nlo->onlyOneLoop() && nlo->needsOLPCorrelators() ) {
- id = orderOLPProcess(nlo->subProcess(),
- (*born).matchboxAmplitude(),
- ProcessType::colourCorrelatedME2);
- nlo->olpProcess(ProcessType::colourCorrelatedME2,id);
+ }
+
+ nlo->needsCorrelations();
+ nlo->cloneDependencies();
+ if ( theonlyk == -1 || theonlyk == i - 1|| theonlyk == i - 2 ) {
+ MEs().push_back(nlo);
+ }
+}
+
+void MFactory::pushProR(Ptr<MatchboxMEBase>::ptr born, int i) {
+ Ptr<MatchboxMEBase>::ptr bornme = (*born).cloneMe();
+
+ string pname = fullName() + "/" + bornme->name() + ".Real";
+ if ( !(generator()->preinitRegister(bornme, pname)) ) throw InitException() << "Subtracted ME " << pname << " already existing.";
+
+
+ Ptr<Node>::ptr clusternode = new_ptr(Node(bornme, i - 2, 1));
+ clusternode->subtractedReal(true);
+ theMergingHelper->firstNodeMap(bornme,clusternode);
+ bornme->merger(theMergingHelper);
+
+ vector<Ptr<Node>::ptr> temp;
+ vector<Ptr<Node>::ptr> temp1;
+ temp.push_back(clusternode);
+
+ unsigned int k = 1;
+ while (thePureMEsMap[i - k].size() != 0) {
+ for ( unsigned int j = 0 ; j < temp.size() ; j++ ) {
+ temp[j]->birth(thePureMEsMap[i - k]);
+ for ( unsigned int m = 0 ; m < temp[j]->children().size() ; ++m ) {
+ temp1.push_back(temp[j]->children()[m]);
}
}
+ temp = temp1;
+ k++;
+ }
+
+ if(MH()->N()>i)
+ bornme->needsCorrelations();
+ else bornme->needsNoCorrelations();
+
+ bornme->cloneDependencies(pname);
+
+ if ( theonlyk == -1 || theonlyk == i - 1 || theonlyk == i - 2 ) {
+ MEs().push_back(bornme);
+ }
+}
+
+void MFactory::orderOLPs() {
+
+}
+
+void MFactory::setup() {
+
+ useMe();
+
+ if(!ransetup){
-
-
-
- nlo->needsCorrelations();
- nlo->cloneDependencies();
- if ( theonlyk == -1 || theonlyk == i - pro|| theonlyk == i - pro-1 ) {
- MEs().push_back(nlo);
- }
- }
-
- void MFactory::pushProR(Ptr<MatchboxMEBase>::ptr born, int i, bool projected,bool finiteDipoles) {
- Ptr<MatchboxMEBase>::ptr bornme = (*born).cloneMe();
- bornme->maxMultCKKW(1);
- bornme->minMultCKKW(0);
- //Ptr<SubtractedME>::ptr sub = new_ptr(SubtractedME());
- int pro = projected ? 2 : 1;
-
- string pname = fullName() + "/" + bornme->name() + (!finiteDipoles?".proReal":"finDI");
- if ( projected ) pname += "pro";
-
- if ( !(generator()->preinitRegister(bornme, pname)) ) throw InitException() << "Subtracted ME " << pname << " already existing.";
-
- if ( projected ) {
- bornme->projectorStage(2);
- } else {
- bornme->projectorStage(1);
- }
-
- Ptr<Node>::ptr clusternode = new_ptr(Node(bornme, i - 2, 1, needFOH));
- clusternode->mergePt(theMergePT);
- clusternode->centralMergePt(theMergePT);
- clusternode->N(theN + getProcesses()[0].size());clusternode->N0( getProcesses()[0].size());
- clusternode->M(theM + getProcesses()[0].size());
- if(theonlyk!=-1)clusternode->onlyN(theonlyk + getProcesses()[0].size());
- clusternode->calculateInNode(true);
- clusternode->subtractedReal(true);
- clusternode->finiteDipoles(finiteDipoles);
- clusternode->unitarized(unitarized);
- clusternode->NLOunitarized(NLOunitarized);
- clusternode->treefactory(this);
- //bornme->firstNode(clusternode);
- theMergingHelper->firstNodeMap(bornme,clusternode);
- bornme->merger(theMergingHelper);
- if ( projected ) bornme->projectorStage(2);
- else bornme->projectorStage(1);
- //bornme->firstNode()->deepHead(clusternode);
- clusternode->chooseHistory(theChooseHistory);
- vector<Ptr<Node>::ptr> temp;
- vector<Ptr<Node>::ptr> temp1;
- temp.push_back(clusternode);
- unsigned int k = 1;
- while (thePureMEsMap[i - k].size() != 0) {
- for ( unsigned int j = 0 ; j < temp.size() ; j++ ) {
- temp[j]->birth(thePureMEsMap[i - k]);
- for ( unsigned int m = 0 ; m < temp[j]->children().size() ; ++m ) {
-
- // if ( temp[j]->children()[m]->nodeME()->isOLPTree() ) {
- // int id = orderOLPProcess(temp[j]->children()[m]->nodeME()->subProcess(),
- // temp[j]->children()[m]->nodeME()->matchboxAmplitude(),
- // ProcessType::treeME2);
- // temp[j]->children()[m]->nodeME()->olpProcess(ProcessType::treeME2,id);
- // id = orderOLPProcess(temp[j]->children()[m]->nodeME()->subProcess(),
- // temp[j]->children()[m]->nodeME()->matchboxAmplitude(),
- // ProcessType::colourCorrelatedME2);
- // temp[j]->children()[m]->nodeME()->olpProcess(ProcessType::colourCorrelatedME2,id);
- //
- // bool haveGluon = false;
- // for ( PDVector::const_iterator p = temp[j]->children()[m]->nodeME()->subProcess().legs.begin();
- // p != temp[j]->children()[m]->nodeME()->subProcess().legs.end(); ++p )
- // if ( (**p).id() == 21 ) {
- // haveGluon = true;
- // break;
- // }
- // if ( haveGluon ) {
- // id = orderOLPProcess(temp[j]->children()[m]->nodeME()->subProcess(),
- // temp[j]->children()[m]->nodeME()->matchboxAmplitude(),
- // ProcessType::spinColourCorrelatedME2);
- // temp[j]->children()[m]->nodeME()->olpProcess(ProcessType::spinColourCorrelatedME2,id);
- // }
- //
- //
- // }
-
-
- temp1.push_back(temp[j]->children()[m]);
- }
- }
- temp = temp1;
- k++;
- }
-
-
- // if ( bornme->isOLPTree() ) {
- // int id = orderOLPProcess(bornme->subProcess(),
- // (*born).matchboxAmplitude(),
- // ProcessType::treeME2);
- // bornme->olpProcess(ProcessType::treeME2,id);
- // id = orderOLPProcess(bornme->subProcess(),
- // (*born).matchboxAmplitude(),
- // ProcessType::colourCorrelatedME2);
- // bornme->olpProcess(ProcessType::colourCorrelatedME2,id);
- //
- // bool haveGluon = false;
- // for ( PDVector::const_iterator p = (*born).subProcess().legs.begin();
- // p != (*born).subProcess().legs.end(); ++p )
- // if ( (**p).id() == 21 ) {
- // haveGluon = true;
- // break;
- // }
- // if ( haveGluon ) {
- // id = orderOLPProcess((*born).subProcess(),
- // (*born).matchboxAmplitude(),
- // ProcessType::spinColourCorrelatedME2);
- // (*born).olpProcess(ProcessType::spinColourCorrelatedME2,id);
- // }
- //
- //
- // }
-
-
-
-
-
-
-
- if(theN>i)
- bornme->needsCorrelations();
- else bornme->needsNoCorrelations();
-
- bornme->cloneDependencies(pname);
-
- ///////////////////////////////REAL/////////////////////////
-
- //sub->setBorns(thePureMEsMap[i - 1]);
- //sub->head(bornme);
- //sub->factory(this);
- //sub->dependent().clear();
- //sub->getDipoles();
- //if ( sub->dependent().empty() ) {
- // cout << "???TreeMerge????" << flush;
- // abort();
- //}
-
- ///////////////////////////////REAL/////////////////////////
-
- //if ( projected ) {
-
- // Ptr<ReweightNumber>::ptr rw1 = new_ptr(ReweightNumber(1.));
- //sub->head()->addReweighter(rw1);
- //bornme->addReweighter(rw1);
- //} else {
- //Ptr<ReweightNumber>::ptr rw1 = new_ptr(ReweightNumber(1.));
- //sub->head()->addReweighter(rw1);
- //bornme->addReweighter(rw1);
- // }
- if ( theonlyk == -1 || theonlyk == i - pro || theonlyk == i - pro-1 ) {
- MEs().push_back(bornme);
- }
- }
-
- void MFactory::orderOLPs() {
-
- }
-
- void MFactory::setup() {
-
- useMe();
-
- if(!ransetup){
-
olpProcesses().clear();
externalAmplitudes().clear();
setFixedCouplings(true);
setFixedQEDCouplings(true);
const PDVector& partons = particleGroups()["j"];
unsigned int nl = 0;
for ( PDVector::const_iterator p = partons.begin();
p != partons.end(); ++p ) {
if ( abs((**p).id()) < 7 && (**p).mass() == ZERO )
++nl;
if ( (**p).id() > 0 && (**p).id() < 7 && (**p).mass() == ZERO )
nLightJetVec( (**p).id() );
if ( (**p).id() > 0 && (**p).id() < 7 && (**p).mass() != ZERO )
nHeavyJetVec( (**p).id() );
}
nLight(nl/2);
const PDVector& partonsInP = particleGroups()["p"];
for ( PDVector::const_iterator pip = partonsInP.begin();
pip != partonsInP.end(); ++pip ) {
if ( (**pip).id() > 0 && (**pip).id() < 7 && (**pip).mass() == ZERO )
nLightProtonVec( (**pip).id() );
}
for ( vector<Ptr<MatchboxAmplitude>::ptr>::iterator amp
= amplitudes().begin(); amp != amplitudes().end(); ++amp )
(**amp).factory(this);
- largeNBasis()->factory(this);
- // assert(!(theonlyk!=-1&&theonlysub!=-1));
- // assert(!(theonlyk!=-1&&divideSub!=-1));
+ theMergingHelper->largeNBasis()->factory(this);
+
assert(!(divideSub!=-1&&divideSubNumber==-1)||!(divideSub==-1&&divideSubNumber!=-1));
assert(!subProcessGroups());
//fill the amplitudes
if ( !amplitudes().empty() ) {
fill_amplitudes();
}
// prepare the Born and virtual matrix elements
- for ( int i = 0 ; i <= max(0, theN) ; ++i ) {
- prepare_BV(i);
- }
+ for ( int i = 0 ; i <= max(0, MH()->N()) ; ++i ) prepare_BV(i);
+
// prepare the real emission matrix elements
- for ( int i = 0 ; i <= theN ; ++i ) {
- prepare_R(i);
- }
+ for ( int i = 0 ; i <= MH()->N() ; ++i ) prepare_R(i);
+
orderOLPs();
// start creating matrix elements
MEs().clear();
+
int onlysubcounter=0;
int i = theonlyabove ;
- cout<<"\nstart filling LO "<<i<<" "<<theonlymulti<<flush;
- if (true) {
-
-
- for (; i <= max(0, theN) ; ++i ) {
- //if(i==theonlymulti||theonlymulti==-1)
+ if (calc_born) {
+ for (; i <= max(0, MH()->N()) ; ++i ) {
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator born = thePureMEsMap[i].begin() ; born != thePureMEsMap[i].end() ; ++born ) {
- if (!theonlyNLOParts&&!theonlyUnlopsweights&&( theonlyk == -1 || theonlyk == i -1|| theonlyk == i )){
+ if (!theonlyUnlopsweights&&( theonlyk == -1 || theonlyk == i -1|| theonlyk == i )){
if(((theonlysub==-1||theonlysub==onlysubcounter)&&divideSub==-1)||(divideSub!=-1&&onlysubcounter%divideSub==divideSubNumber))
- pushB(*born, i, false);
+ pushB(*born, i);
onlysubcounter++;
}
}
+ }
}
- }
- cout<<"\nstart filling NLO"<<flush;
-
-
-
- i = theonlyabove ;
- for (; i <=max(0, theN); ++i ) {
- // if(i==theonlymulti||theonlymulti==-1)
+
+
+ if (calc_virtual) {
+ i = theonlyabove ;
+ for (; i <=max(0, MH()->N()); ++i ) {
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator born = thePureMEsMap[i].begin() ; born != thePureMEsMap[i].end() ; ++born ) {
-
- if ( i <= theM && !theonlyrealNLOParts&&( theonlyk == -1 || theonlyk == i -1|| theonlyk == i )){
+ if ( i <= MH()->M() &&( theonlyk == -1 || theonlyk == i -1|| theonlyk == i )){
if(((theonlysub==-1||theonlysub==onlysubcounter)&&divideSub==-1)||(divideSub!=-1&&onlysubcounter%divideSub==divideSubNumber))
- pushV(*born, i, false);
+ pushV(*born, i);
onlysubcounter++;
}
}}
- if (false) {
- i = theonlyabove ;
- for (; i <= max(0, theN) ; ++i ) {
- //if((i)==theonlymulti+1||theonlymulti==-1)
+ }
+ if (calc_real) {
+ i = theonlyabove ;
+ for (; i <= max(0, MH()->N()) ; ++i ) {
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator born = thePureMEsMap[i].begin() ; born != thePureMEsMap[i].end() ; ++born ) {
- if ( i <= theM + 1 && i > 0 && !theonlyvirtualNLOParts&&!theonlyUnlopsweights&&( theonlyk == -1 || theonlyk == i -2|| theonlyk == i-1 )){
+ if ( i <= MH()->M() + 1 && i > 0 && !theonlyvirtualNLOParts&&!theonlyUnlopsweights&&( theonlyk == -1 || theonlyk == i -2|| theonlyk == i-1 )){
if(((theonlysub==-1||theonlysub==onlysubcounter)&&divideSub==-1)||(divideSub!=-1&&onlysubcounter%divideSub==divideSubNumber))
- pushProR(*born, i, false,false);
+ pushProR(*born, i);
onlysubcounter++;
}
}
+ }
}
- }
-
+
if ( !externalAmplitudes().empty() ) {
generator()->log() << "Initializing external amplitudes.\n" << flush;
boost::progress_display * progressBar =
new boost::progress_display(externalAmplitudes().size(),generator()->log());
for ( set<Ptr<MatchboxAmplitude>::tptr>::const_iterator ext =
externalAmplitudes().begin(); ext != externalAmplitudes().end(); ++ext ) {
if ( !(**ext).initializeExternal() ) {
throw InitException()
<< "error: failed to initialize amplitude '" << (**ext).name() << "'\n";
}
++(*progressBar);
}
delete progressBar;
generator()->log() << "--------------------------------------------------------------------------------\n"
<< flush;
}
if ( !olpProcesses().empty() ) {
generator()->log() << "Initializing one-loop provider(s).\n" << flush;
map<Ptr<MatchboxAmplitude>::tptr, map<pair<Process, int>, int> > olps;
for ( map<Ptr<MatchboxAmplitude>::tptr, map<pair<Process, int>, int> >::const_iterator oit = olpProcesses().begin() ; oit != olpProcesses().end() ;
++oit ) {
olps[oit->first] = oit->second;
}
boost::progress_display * progressBar = new boost::progress_display(olps.size(), generator()->log());
for ( map<Ptr<MatchboxAmplitude>::tptr, map<pair<Process, int>, int> >::const_iterator olpit = olps.begin() ; olpit != olps.end() ; ++olpit ) {
if ( !olpit->first->startOLP(olpit->second) ) {
throw InitException() << "error: failed to start OLP for amplitude '" << olpit->first->name() << "'\n";
}
++(*progressBar);
}
delete progressBar;
generator()->log() << "--------------------------------------------------------------------------------\n" << flush;
}
cout<<"\n Generated "<<MEs().size()<<" Subprocesses."<<flush;
if(theonlysub!=-1)cout<<" ( "<<theonlysub<<"/"<<onlysubcounter<<" )"<<flush;
cout<<"\n"<<flush;
generator()->log() << "process setup finished.\n" << flush;
-
-
-
- ransetup=true;
-
- }
-
-
- }
-
- void MFactory::persistentOutput(PersistentOStream & os) const {
- os<< theonlyNLOParts
- << theonlyvirtualNLOParts
- << theonlyrealNLOParts
- << theunitarizeNLOParts
- << calc_born
- << calc_virtual
- << calc_real
- << calc_projected_virtual
- << calc_inclusive_real
- << calc_double_inclusive
- << needFOH
- << unitarized
- << NLOunitarized
- << theonlyUnlopsweights
- << theonlyk
- << theonlymulti
- << divideSub
- << divideSubNumber
- << theonlyabove
- << theonlysub
- << theChooseHistory
- << theN
- << theM
- << theMergePTsmearing
- << theIRSafeRATIO
- << theStairFactor
- << ounit(theMergePT,GeV)
- << ounit(theNLOMergePT,GeV)
- << ounit(theIRSafePT,GeV)
- << theSubtractionData
- << theLargeNBasis
- << ransetup
- << processMap
- <<defMERegionByJetAlg
- <<theCuts
- <<theMergingJetFinder
- <<pp_dcut
- <<ee_ycut
- <<theGamma<<theMergingHelper;
- }
-
- void MFactory::persistentInput(PersistentIStream & is, int) {
- is>> theonlyNLOParts
- >> theonlyvirtualNLOParts
- >> theonlyrealNLOParts
- >> theunitarizeNLOParts
- >> calc_born
- >> calc_virtual
- >> calc_real
- >> calc_projected_virtual
- >> calc_inclusive_real
- >> calc_double_inclusive
- >> needFOH
- >> unitarized
- >> NLOunitarized
- >> theonlyUnlopsweights
- >> theonlyk
- >> theonlymulti
- >> divideSub
- >> divideSubNumber
- >> theonlyabove
- >> theonlysub
- >> theChooseHistory
- >> theN
- >> theM
- >> theMergePTsmearing
- >> theIRSafeRATIO
- >> theStairFactor
- >> iunit(theMergePT,GeV)
- >> iunit(theNLOMergePT,GeV)
- >> iunit(theIRSafePT,GeV)
- >> theSubtractionData
- >> theLargeNBasis
- >> ransetup
- >> processMap
- >>defMERegionByJetAlg
- >> theCuts
- >>theMergingJetFinder
- >>pp_dcut
- >>ee_ycut
- >>theGamma>>theMergingHelper;
- }
-
- void MFactory::Init() {
- static Parameter<MFactory, int> interfaceadditionalN("additionalN", "Set the number of additional jets to consider.", &MFactory::theN, 0, 0,
- 0, false, false, Interface::lowerlim);
-
- static Parameter<MFactory, int> interfacevirtualM("virtualM",
- "Set the number of virtual corrections to consider. -1 is default for no virtual correction.", &MFactory::theM, -1, -1, 0, false, false,
- Interface::lowerlim);
-
- static Parameter<MFactory, int> interfaceonlyk("onlyk", "calculate only the ME with k additional partons.", &MFactory::theonlyk, -1, -1, 0,
- false, false, Interface::lowerlim);
- static Parameter<MFactory, int> interfaceonlymulti("onlymulti", "calculate only the ME with k additional partons.", &MFactory::theonlymulti, -1, -1, 0,
- false, false, Interface::lowerlim);
- static Parameter<MFactory, int> interfaceonlyabove("onlyabove", "calculate only the ME with more than k additional partons.", &MFactory::theonlyabove, -1, -1, 0,
- false, false, Interface::lowerlim);
-
- static Parameter<MFactory, int> interfaceonlysub("onlysub", "calculate only one subProcess. this is for building grids.", &MFactory::theonlysub, -1, -1, 0,
- false, false, Interface::lowerlim);
-
-
-
-
-
-
-
-
-
- static Parameter<MFactory, int> interfacedivideSub("divideSub", "calculate only one subProcess. this is for building grids.", &MFactory::divideSub, -1, -1, 0,
- false, false, Interface::lowerlim);
-
-
- static Parameter<MFactory, int> interfacedivideSubNumber("divideSubNumber", "calculate only one subProcess. this is for building grids.", &MFactory::divideSubNumber, -1, -1, 0,
- false, false, Interface::lowerlim);
-
-
-
-
-
- static Parameter<MFactory, int> interfacechooseHistory("chooseHistory", "different ways of choosing the history", &MFactory::theChooseHistory, 3, -1, 0,
- false, false, Interface::lowerlim);
-
-
- static Switch<MFactory, bool> interface_calc_born("calc_born", "[debug] Switch on or off the born contribution.", &MFactory::calc_born, true,
- false, false);
- static SwitchOption interface_calc_bornOn(interface_calc_born, "On", "Switch on calculation of born.", true);
- static SwitchOption interface_calc_bornOff(interface_calc_born, "Off", "Switch off calculation of born.", false);
-
- static Switch<MFactory, bool> interface_calc_virtual("calc_virtual", "[debug] Switch on or off the virtual contribution.",
- &MFactory::calc_virtual, true, false, false);
- static SwitchOption interface_calc_virtualOn(interface_calc_virtual, "On", "Switch on calculation of virtual.", true);
- static SwitchOption interface_calc_virtualOff(interface_calc_virtual, "Off", "Switch off calculation of virtual.", false);
-
- static Switch<MFactory, bool> interface_calc_real("calc_real", "[debug] Switch on or off the real contribution.", &MFactory::calc_real, true,
- false, false);
- static SwitchOption interface_calc_realOn(interface_calc_real, "On", "Switch on calculation of real.", true);
- static SwitchOption interface_calc_realOff(interface_calc_real, "Off", "Switch off calculation of real.", false);
-
- static Switch<MFactory, bool> interface_calc_projected_virtual("calc_projected_virtual",
- "[debug] Switch on or off the projected_virtual contribution.", &MFactory::calc_projected_virtual, true, false, false);
- static SwitchOption interface_calc_projected_virtualOn(interface_calc_projected_virtual, "On", "Switch on calculation of projected_virtual.", true);
- static SwitchOption interface_calc_projected_virtualOff(interface_calc_projected_virtual, "Off", "Switch off calculation of projected_virtual.", false);
-
- static Switch<MFactory, bool> interface_calc_inclusive_real("calc_inclusive_real", "[debug] Switch on or off the real contribution.",
- &MFactory::calc_inclusive_real, true, false, false);
- static SwitchOption interface_calc_inclusive_realOn(interface_calc_inclusive_real, "On", "Switch on calculation of inclusive_real.", true);
- static SwitchOption interface_calc_inclusive_realOff(interface_calc_inclusive_real, "Off", "Switch off calculation of inclusive_real.", false);
-
- static Switch<MFactory, bool> interface_calc_double_inclusive("calc_double_inclusive", "[debug] Switch on or off the real contribution.",
- &MFactory::calc_double_inclusive, true, false, false);
- static SwitchOption interface_calc_double_inclusiveOn(interface_calc_double_inclusive, "On", "Switch on calculation of double_inclusive.", true);
- static SwitchOption interface_calc_double_inclusiveOff(interface_calc_double_inclusive, "Off", "Switch off calculation of double_inclusive.", false);
-
- static Switch<MFactory, bool> interface_needFOH("needFOH", "Switch on or off the full ordered histories.", &MFactory::needFOH, true, false,
- false);
- static SwitchOption interface_needFOHOn(interface_needFOH, "On", "Switch on the full ordered histories.", true);
- static SwitchOption interface_needFOHOff(interface_needFOH, "Off", "Switch off the full ordered histories.", false);
-
- static Switch<MFactory, bool> interface_theonlyNLOParts("onlyNLOParts", "Switch on or off the onlyNLOParts.", &MFactory::theonlyNLOParts, true, false,
- false);
- static SwitchOption interface_theonlyNLOPartsOn(interface_theonlyNLOParts, "On", "Switch on the theonlyNLOParts.", true);
- static SwitchOption interface_theonlyNLOPartsOff(interface_theonlyNLOParts, "Off", "Switch off the theonlyNLOParts.", false);
-
-
- // theonlyvirtualNLOParts;
- // theonlyrealNLOParts;
-
- static Switch<MFactory, bool> interface_theonlyvirtualNLOParts("onlyvirtualNLOParts", "Switch on or off the onlyvirtualNLOParts.", &MFactory::theonlyvirtualNLOParts, true, false,
- false);
- static SwitchOption interface_theonlyvirtualNLOPartsOn(interface_theonlyvirtualNLOParts, "On", "Switch on the theonlyvirtualNLOParts.", true);
- static SwitchOption interface_theonlyvirtualNLOPartsOff(interface_theonlyvirtualNLOParts, "Off", "Switch off the theonlyvirtualNLOParts.", false);
-
- static Switch<MFactory, bool> interface_theonlyrealNLOParts("onlyrealNLOParts", "Switch on or off the onlyrealNLOParts.", &MFactory::theonlyrealNLOParts, true, false,
- false);
- static SwitchOption interface_theonlyrealNLOPartsOn(interface_theonlyrealNLOParts, "On", "Switch on the theonlyrealNLOParts.", true);
- static SwitchOption interface_theonlyrealNLOPartsOff(interface_theonlyrealNLOParts, "Off", "Switch off the theonlyrealNLOParts.", false);
-
- static Switch<MFactory, bool> interface_theunitarizeNLOParts("unitarizeNLOParts", "Switch on or off the unitarizeNLOParts.", &MFactory::theunitarizeNLOParts, true, false,
- false);
- static SwitchOption interface_theunitarizeNLOPartsOn(interface_theunitarizeNLOParts, "On", "Switch on the unitarizeNLOParts.", true);
- static SwitchOption interface_theunitarizeNLOPartsOff(interface_theunitarizeNLOParts, "Off", "Switch off the unitarizeNLOParts.", false);
-
-
- static Switch<MFactory, bool> interface_theonlyUnlopsweights("onlyUnlopsweights", "Switch on or off the onlyUnlopsweights.", &MFactory::theonlyUnlopsweights, true, false,
- false);
- static SwitchOption interface_theonlyUnlopsweightsOn(interface_theonlyUnlopsweights, "On", "Switch on the onlyUnlopsweights.", true);
- static SwitchOption interface_theonlyUnlopsweightsOff(interface_theonlyUnlopsweights, "Off", "Switch off the onlyUnlopsweights.", false);
-
-
-
- static Switch<MFactory, bool> interface_Unitarized("Unitarized", "Unitarize the cross section.", &MFactory::unitarized, true, false, false);
- static SwitchOption interface_UnitarizedOn(interface_Unitarized, "On", "Switch on the unitarized cross section.", true);
- static SwitchOption interface_UnitarizedOff(interface_Unitarized, "Off", "Switch off the unitarized cross section.", false);
-
- static Switch<MFactory, bool> interface_NLOUnitarized("NLOUnitarized", "Unitarize the cross section.", &MFactory::NLOunitarized, true, false, false);
- static SwitchOption interface_NLOUnitarizedOn(interface_NLOUnitarized, "On", "Switch on the unitarized NLO cross section.", true);
- static SwitchOption interface_NLOUnitarizedOff(interface_NLOUnitarized, "Off", "Switch off the unitarized NLO cross section.", false);
-
- static Parameter<MFactory, Energy> interfacemergept("mergept", "Set the pt to be merged to consider.", &MFactory::theMergePT, GeV, 2.0 * GeV,
- ZERO, Constants::MaxEnergy, true, false, Interface::limited);
-
-
-
- static Parameter<MFactory, double> interfacemergeptsmearing("mergeptsmearing", "Set the percentage the merging pt should be smeared.",
- &MFactory::theMergePTsmearing, 0., 0.,
- 0.0, 0.5, true, false, Interface::limited);
-
- static Parameter<MFactory, double> interfacestairfactor("stairfactor", "Set the pt to be merged to consider.", &MFactory::theStairFactor, 1, 1,
- 1, 4, true, false, Interface::limited);
- interfacestairfactor.setHasDefault(false);
-
-
-
- interfacemergept.setHasDefault(false);
- static Parameter<MFactory, Energy> interfaceNLOmergept("NLOmergept", "NLO:Set the pt to be merged to consider.", &MFactory::theNLOMergePT, GeV, 2.0 * GeV,
- ZERO, Constants::MaxEnergy, true, false, Interface::limited);
- interfaceNLOmergept.setHasDefault(false);
-
- static Parameter<MFactory, Energy> interfaceIRSafePT("IRSafePT", "Set the pt for which a matrix element should be treated as IR-safe.",
- &MFactory::theIRSafePT, GeV, 0.0 * GeV, ZERO, Constants::MaxEnergy, true, false, Interface::limited);
- interfaceIRSafePT.setHasDefault(false);
-
- static Parameter<MFactory, double> interfaceIRSafeRATIO("IRSafeRATIO", "Set the RATIO for which a matrix element should be treated as IR-safe.",
- &MFactory::theIRSafeRATIO, 100.0, 0.0, ZERO, false, false, Interface::lowerlim);
- interfaceIRSafeRATIO.setHasDefault(false);
-
- /* static Command<MFactory> interfaceProcess("Process", "Set the process to consider.", &MFactory::doProcess, false);*/
-
- static Parameter<MFactory, string> interfaceSubtractionData("SubtractionData", "Prefix for subtraction check data.",
- &MFactory::theSubtractionData, "", false, false);
-
-
- static Reference<MFactory,ColourBasis> interfaceLargeNBasis
- ("LargeNBasis",
- "Set the large-N colour basis implementation.",
- &MFactory::theLargeNBasis, false, false, true, true, false);
-
-
-
-
-
-
- static Switch<MFactory,bool>
- interfacedefMERegionByJetAlg ("MERegionByJetAlg","",&MFactory::defMERegionByJetAlg, false, false, false);
-
- static SwitchOption interfacedefMERegionByJetAlgYes
- (interfacedefMERegionByJetAlg,"Yes","",true);
- static SwitchOption interfacedefMERegionByJetAlgNo
- (interfacedefMERegionByJetAlg,"No","",false);
-
-
- static Reference<MFactory,JetFinder> interfaceMergingJetFinder
- ("MergingJetFinder",
- "A reference to the jet finder used in Merging.",
- &MFactory::theMergingJetFinder, false, false, true, false, false);
-
- static Reference<MFactory,Cuts> interfaceCuts
- ("Cuts",
- "A reference to the Cuts.",
- &MFactory::theCuts, false, false, true, false, false);
-
- static Parameter<MFactory,double> interfacedcut
- ("pp_dcut",
- "The MergingScale.",
- &MFactory::pp_dcut, 5.0, 0.0, 0,
- false, false, Interface::lowerlim);
-
- static Parameter<MFactory,double> interfaceycut
- ("ee_ycut",
- "The MergingScale.",
- &MFactory::ee_ycut, 0.2, 0.0, 0,
- false, false, Interface::lowerlim);
-
- static Parameter<MFactory,double> interfacegamma
- ("gamma",
- "gamma parameter.",
- &MFactory::theGamma, 1.0, 0.0, 0,
- false, false, Interface::lowerlim);
-
-
-
-
- static Reference<MFactory,Merger> interfaceMergingHelper3
- ("MergingHelper",
- "",
- &MFactory::theMergingHelper, false, false, true, true, false);
-
-
-
+ ransetup=true;
}
+
+
+}
- // *** 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<MFactory, MatchboxFactory> describeHerwigMFactory("Herwig::MFactory", "HwDipoleShower.so");
+void MFactory::persistentOutput(PersistentOStream & os) const {
+ os
+ << calc_born << calc_virtual << calc_real
+ << theonlyUnlopsweights << theonlyk << theonlymulti
+ << divideSub << divideSubNumber << theonlyabove
+ << theonlysub << theSubtractionData << ransetup
+ << processMap << theMergingHelper;
+}
+
+void MFactory::persistentInput(PersistentIStream & is, int) {
+ is
+ >> calc_born >> calc_virtual >> calc_real
+ >> theonlyUnlopsweights >> theonlyk >> theonlymulti
+ >> divideSub >> divideSubNumber >> theonlyabove
+ >> theonlysub >> theSubtractionData >> ransetup
+ >> processMap >> theMergingHelper;
+}
+
+void MFactory::Init() {
+
+
+
+ static Parameter<MFactory, int> interfaceonlyk("onlyk", "calculate only the ME with k additional partons.", &MFactory::theonlyk, -1, -1, 0,
+ false, false, Interface::lowerlim);
+ static Parameter<MFactory, int> interfaceonlymulti("onlymulti", "calculate only the ME with k additional partons.", &MFactory::theonlymulti, -1, -1, 0,
+ false, false, Interface::lowerlim);
+ static Parameter<MFactory, int> interfaceonlyabove("onlyabove", "calculate only the ME with more than k additional partons.", &MFactory::theonlyabove, -1, -1, 0,
+ false, false, Interface::lowerlim);
+
+ static Parameter<MFactory, int> interfaceonlysub("onlysub", "calculate only one subProcess. this is for building grids.", &MFactory::theonlysub, -1, -1, 0,
+ false, false, Interface::lowerlim);
+
+
+
+
+
+
+
+
+
+ static Parameter<MFactory, int> interfacedivideSub("divideSub", "calculate only one subProcess. this is for building grids.", &MFactory::divideSub, -1, -1, 0,
+ false, false, Interface::lowerlim);
+
+
+ static Parameter<MFactory, int> interfacedivideSubNumber("divideSubNumber", "calculate only one subProcess. this is for building grids.", &MFactory::divideSubNumber, -1, -1, 0,
+ false, false, Interface::lowerlim);
+
+
+
+
+
+ static Switch<MFactory, bool> interface_calc_born("calc_born", "[debug] Switch on or off the born contribution.", &MFactory::calc_born, true,
+ false, false);
+ static SwitchOption interface_calc_bornOn(interface_calc_born, "On", "Switch on calculation of born.", true);
+ static SwitchOption interface_calc_bornOff(interface_calc_born, "Off", "Switch off calculation of born.", false);
+
+ static Switch<MFactory, bool> interface_calc_virtual("calc_virtual", "[debug] Switch on or off the virtual contribution.",
+ &MFactory::calc_virtual, true, false, false);
+ static SwitchOption interface_calc_virtualOn(interface_calc_virtual, "On", "Switch on calculation of virtual.", true);
+ static SwitchOption interface_calc_virtualOff(interface_calc_virtual, "Off", "Switch off calculation of virtual.", false);
+
+ static Switch<MFactory, bool> interface_calc_real("calc_real", "[debug] Switch on or off the real contribution.", &MFactory::calc_real, true,
+ false, false);
+ static SwitchOption interface_calc_realOn(interface_calc_real, "On", "Switch on calculation of real.", true);
+ static SwitchOption interface_calc_realOff(interface_calc_real, "Off", "Switch off calculation of real.", false);
+
+
+
+
+ static Switch<MFactory, bool> interface_theonlyNLOParts("onlyNLOParts", "Switch on or off the onlyNLOParts.", &MFactory::theonlyNLOParts, true, false,
+ false);
+ static SwitchOption interface_theonlyNLOPartsOn(interface_theonlyNLOParts, "On", "Switch on the theonlyNLOParts.", true);
+ static SwitchOption interface_theonlyNLOPartsOff(interface_theonlyNLOParts, "Off", "Switch off the theonlyNLOParts.", false);
+
+
+ // theonlyvirtualNLOParts;
+ // theonlyrealNLOParts;
+
+ static Switch<MFactory, bool> interface_theonlyvirtualNLOParts("onlyvirtualNLOParts", "Switch on or off the onlyvirtualNLOParts.", &MFactory::theonlyvirtualNLOParts, true, false,
+ false);
+ static SwitchOption interface_theonlyvirtualNLOPartsOn(interface_theonlyvirtualNLOParts, "On", "Switch on the theonlyvirtualNLOParts.", true);
+ static SwitchOption interface_theonlyvirtualNLOPartsOff(interface_theonlyvirtualNLOParts, "Off", "Switch off the theonlyvirtualNLOParts.", false);
+
+ static Switch<MFactory, bool> interface_theonlyrealNLOParts("onlyrealNLOParts", "Switch on or off the onlyrealNLOParts.", &MFactory::theonlyrealNLOParts, true, false,
+ false);
+ static SwitchOption interface_theonlyrealNLOPartsOn(interface_theonlyrealNLOParts, "On", "Switch on the theonlyrealNLOParts.", true);
+ static SwitchOption interface_theonlyrealNLOPartsOff(interface_theonlyrealNLOParts, "Off", "Switch off the theonlyrealNLOParts.", false);
+
+ static Switch<MFactory, bool> interface_theunitarizeNLOParts("unitarizeNLOParts", "Switch on or off the unitarizeNLOParts.", &MFactory::theunitarizeNLOParts, true, false,
+ false);
+ static SwitchOption interface_theunitarizeNLOPartsOn(interface_theunitarizeNLOParts, "On", "Switch on the unitarizeNLOParts.", true);
+ static SwitchOption interface_theunitarizeNLOPartsOff(interface_theunitarizeNLOParts, "Off", "Switch off the unitarizeNLOParts.", false);
+
+
+ static Switch<MFactory, bool> interface_theonlyUnlopsweights("onlyUnlopsweights", "Switch on or off the onlyUnlopsweights.", &MFactory::theonlyUnlopsweights, true, false,
+ false);
+ static SwitchOption interface_theonlyUnlopsweightsOn(interface_theonlyUnlopsweights, "On", "Switch on the onlyUnlopsweights.", true);
+ static SwitchOption interface_theonlyUnlopsweightsOff(interface_theonlyUnlopsweights, "Off", "Switch off the onlyUnlopsweights.", false);
+
+
+
+ static Switch<MFactory, bool> interface_Unitarized("Unitarized", "Unitarize the cross section.", &MFactory::unitarized, true, false, false);
+ static SwitchOption interface_UnitarizedOn(interface_Unitarized, "On", "Switch on the unitarized cross section.", true);
+ static SwitchOption interface_UnitarizedOff(interface_Unitarized, "Off", "Switch off the unitarized cross section.", false);
+
+
+
+
+ static Switch<MFactory, bool> interface_NLOUnitarized("NLOUnitarized", "Unitarize the cross section.", &MFactory::NLOunitarized, true, false, false);
+ static SwitchOption interface_NLOUnitarizedOn(interface_NLOUnitarized, "On", "Switch on the unitarized NLO cross section.", true);
+ static SwitchOption interface_NLOUnitarizedOff(interface_NLOUnitarized, "Off", "Switch off the unitarized NLO cross section.", false);
+
+
+
+
+ /* static Command<MFactory> interfaceProcess("Process", "Set the process to consider.", &MFactory::doProcess, false);*/
+
+ static Parameter<MFactory, string> interfaceSubtractionData("SubtractionData", "Prefix for subtraction check data.",
+ &MFactory::theSubtractionData, "", false, false);
+
+
+
+
+ static Reference<MFactory,Merger> interfaceMergingHelper
+ ("MergingHelper",
+ "",
+ &MFactory::theMergingHelper, false, false, true, true, false);
+
+
+
+
+}
+
+ // *** 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<MFactory, MatchboxFactory> describeHerwigMFactory("Herwig::MFactory", "HwDipoleShower.so");
diff --git a/DipoleShower/Merging2/MFactory.h b/DipoleShower/Merging2/MFactory.h
--- a/DipoleShower/Merging2/MFactory.h
+++ b/DipoleShower/Merging2/MFactory.h
@@ -1,338 +1,217 @@
-// -*- C++ -*-
-//
-// MFactory.h is a part of Herwig - A multi-purpose Monte Carlo event generator
-// Copyright (C) 2002-2012 The Herwig Collaboration
-//
-// Herwig is licenced under version 2 of the GPL, see COPYING for details.
-// Please respect the MCnet academic guidelines, see GUIDELINES for details.
-//
+ // -*- C++ -*-
+ //
+ // MFactory.h is a part of Herwig - A multi-purpose Monte Carlo event generator
+ // Copyright (C) 2002-2012 The Herwig Collaboration
+ //
+ // Herwig is licenced under version 2 of the GPL, see COPYING for details.
+ // Please respect the MCnet academic guidelines, see GUIDELINES for details.
+ //
#ifndef HERWIG_MFactory_H
#define HERWIG_MFactory_H
-//
-// This is the declaration of the MergeboxFactory class.
-//
+ //
+ // This is the declaration of the MergeboxFactory class.
+ //
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
#include "Node.fh"
#include "Merger.h"
-//#include "ThePEG/MatrixElement/ReweightConstant.h"
+ //#include "ThePEG/MatrixElement/ReweightConstant.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
namespace Herwig {
-
+
using namespace ThePEG;
-
+
/**
* \ingroup Matchbox
- * \author Simon Platzer,Johannes Bellm
+ * \author Johannes Bellm
*
* \brief MFactory automatically sets up a NLO
* QCD calculation carried out in dipole subtraction.
*
* @see \ref MFactoryInterfaces "The interfaces"
* defined for MergeboxFactory.
*/
class MFactory : public MatchboxFactory {
public:
-
+
/** @name Standard constructors and destructors. */
- //@{
+ //@{
/**
* The default constructor.
*/
MFactory();
-
+
/**
* The destructor.
*/
virtual ~MFactory();
- //@}
-
-
+ //@}
+
virtual void setup();
-
+
/**fill all amplitudes, stored in pureMEsMap
*/
void fill_amplitudes();
-
+
/** prepare the Born and virtual matrix elements.
*/
void prepare_BV(int i);
-
+
void prepare_R(int i);
-
- void pushB(Ptr<MatchboxMEBase>::ptr,int,bool projected);
-
- void pushV(Ptr<MatchboxMEBase>::ptr,int,bool projected);
-
- void pushProR(Ptr<MatchboxMEBase>::ptr,int,bool projected,bool);
-
-
+
+ void pushB(Ptr<MatchboxMEBase>::ptr,int);
+
+ void pushV(Ptr<MatchboxMEBase>::ptr,int);
+
+ void pushProR(Ptr<MatchboxMEBase>::ptr,int);
+
void orderOLPs();
- Ptr<ColourBasis>::ptr largeNBasis(){return theLargeNBasis;}
-
- void largeNBasis(Ptr<ColourBasis>::ptr x){theLargeNBasis=x;}
-
-
- double stairfactor() const {return theStairFactor;}
-
- double smear(){return theMergePTsmearing;}
-
int onlymulti()const {return theonlymulti==-1?-1:(theonlymulti+processMap.find(0)->second.size());}
bool onlyUnlopsweights() const {return theonlyUnlopsweights;}
- Ptr<Merger>::ptr mergingHelper(){return theMergingHelper;};
+ Ptr<Merger>::ptr MH(){return theMergingHelper;}
-
+
/**
* Return the Map of matrix elements to be considered
* (the Key is the number of additional jets)
*/
const map<int, vector<Ptr<MatchboxMEBase>::ptr> >& pureMEsMap() const {
return thePureMEsMap;
}
-
+
/**
* Access the Map of matrix elements to be considered
* (the Key is the number of additional jets)
*/
map<int, vector<Ptr<MatchboxMEBase>::ptr> >& pureMEsMap() {
return thePureMEsMap;
}
+ public:
-
-
- bool MERegionByJetAlg(){return defMERegionByJetAlg;}
-
- double gamma()const{return theGamma;}
-
- /**
- * Return the produced subtracted matrix elements
- */
- // const vector<Ptr<MatchboxMEBase>::ptr>& projectedMEs() const {
- // return theprojectedMEs;
- // }
- //
- // vector<Ptr<MatchboxMEBase>::ptr>& projectedMEs() {
- // return theprojectedMEs;
- // }
- 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);
- //@}
-
-
-
-
-
-
-
-
+ //@}
+
+
static void Init();
-
-
-
-
-
+
+
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;
- //@}
-
-
-
+ //@}
+
+
+
private:
-
+
bool theonlyNLOParts;
bool theonlyvirtualNLOParts;
bool theonlyrealNLOParts;
+ bool theonlyUnlopsweights;
bool theunitarizeNLOParts;
bool calc_born;
bool calc_virtual;
bool calc_real;
- bool calc_projected_virtual;
- bool calc_inclusive_real;
- bool calc_double_inclusive;
- bool needFOH;//Full ordered history
bool unitarized;
bool NLOunitarized;
- bool theonlyUnlopsweights;
- bool defMERegionByJetAlg;
int theonlyk;
int theonlymulti;
+ int theonlyabove;
+ int theonlysub;
int divideSub;
int divideSubNumber;
- int theonlyabove;
- int theonlysub;
- int theChooseHistory;
- int theN;
- int theM;
- double theMergePTsmearing;
- double theIRSafeRATIO;
- double theStairFactor;
- Energy theMergePT;
- Energy theNLOMergePT;
- Energy theIRSafePT;
-
- /**
- * Reference to the jet finder
- */
- Ptr<JetFinder>::ptr theMergingJetFinder;
- CutsPtr theCuts;
- double ee_ycut;
-
- double pp_dcut;
-
/**
* Prefix for subtraction data
*/
string theSubtractionData;
-
-
-
+
+
map< int, vector<string> > processMap;
+
map< int, bool > thevirtualsAreExpandedMap;
-
-
- /**
- * True, if the integral over the unresolved emission should be
- * calculated.
- */
- map<int, bool > theInclusiveMap;
-
-
- /**
- * True, if veto scales should be set
- * for the real emission
- */
- map<int, bool> theVetoScalesMap;
-
+
+
/**
* The matrix elements: int = number of additional jets
*/
-
+
map< int, vector<Ptr<MatchboxMEBase>::ptr> > thePureMEsMap;
-
+
/**
* The Born matrix elements to be considered
*/
vector<Ptr<MatchboxMEBase>::ptr> theBornMEs;
-
+
/**
* The virtual corrections to be considered
*/
map<int, vector<Ptr<MatchboxInsertionOperator>::ptr> > theVirtualsMap;
-
+
/**
* The produced NLO matrix elements
*/
vector<Ptr<MatchboxMEBase>::ptr> theBornVirtualMEs;
-
+
/**
* The real emission matrix elements to be considered
*/
-
+
vector<Ptr<MatchboxMEBase>::ptr> theRealEmissionMEs;
-
-
+
/**
- * The produced subtracted matrix elements
+ * A large-N colour basis to be used when reproducing the shower
+ * kernels.
*/
- vector<Ptr<SubtractedME>::ptr> theSubtractedMEs;
-
- /**
- * The produced finite real emission matrix elements
- */
- vector<Ptr<MatchboxMEBase>::ptr> theFiniteRealMEs;
-
-
-
-
- /**
- * The produced finite real emission matrix elements
- */
- map<int, vector<Ptr<MatchboxMEBase>::ptr> > the_Fin_R_MEsMap;
-
-
-
-
- /**
- * A large-N colour basis to be used when reproducing the shower
- * kernels.
- */
- Ptr<ColourBasis>::ptr theLargeNBasis;
-
Ptr<Merger>::ptr theMergingHelper;
-
-
- bool ransetup;
+ bool ransetup;
- double theGamma;
-
-/*
-
- *
- * Generate subprocesses.
-
- set<PDVector> makeSubProcesses(const vector<string>&) const;
-
- *
- * Generate matrix element objects for the given process.
-
- vector<Ptr<MatchboxMEBase>::ptr> makeMEs(const vector<string>&,
- unsigned int orderas) const;
-
-*/
-
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MFactory & operator=(const MFactory &);
-
-
-
-
-
+
};
-
+
}
#endif /* HERWIG_MFactory_H */
diff --git a/DipoleShower/Merging2/Merger.cc b/DipoleShower/Merging2/Merger.cc
--- a/DipoleShower/Merging2/Merger.cc
+++ b/DipoleShower/Merging2/Merger.cc
@@ -1,1421 +1,1494 @@
// -*- C++ -*-
//
// Merger.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig is licenced under version 2 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 Merger class.
//
#include "Merger.h"
#include "Node.h"
#include "MFactory.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "Herwig/DipoleShower/DipoleShowerHandler.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Config/Constants.h"
#include "Herwig/Shower/ShowerHandler.h"
#include "Herwig/PDF/HwRemDecayer.h"
#include "Herwig/MatrixElement/Matchbox/Phasespace/RandomHelpers.h"
#include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
#include "fastjet/ClusterSequence.hh"
using namespace Herwig;
Merger::Merger()
: MergerBase() {
StartingBorn0=NPtr();
StartingBorn1=NPtr();
StartingBorn2=NPtr();
NPtr CalcBorn0=NPtr();
NPtr CalcBorn1=NPtr();
NPtr CalcBorn2=NPtr();
theNf=5;
minusL=false;
Unlopsweights=true;
theKImproved=true;
MergingScale=20.*GeV;
theDipoleShowerHandler=Ptr<DipoleShowerHandler>::ptr();
}
Merger::~Merger() {}
IBPtr Merger::clone() const {
return new_ptr(*this);
}
IBPtr Merger::fullclone() const {
return new_ptr(*this);
}
double Merger::reweightCKKWBornGamma(NPtr Node,bool fast){
double weightCL=0.;
weight=1.;
if (!Node->children().empty()) {
PVector particles;
for (size_t i=0;i<Node->nodeME()->mePartonData().size();i++){
Ptr<ThePEG::Particle>::ptr p =Node->nodeME()->mePartonData()[i]->produceParticle(Node->nodeME()->lastMEMomenta()[i]);
particles.push_back(p);
}
if(!matrixElementRegion(particles, Node->children()[0]->dipol()->lastPt(), MergingScale))weight*= 0.;
}
NPtr Born= Node-> getLongestHistory_simple(true,xiQSh);//Longest irreführend
NPtr CLNode;
NPtr BornCL;
if( !Node->children().empty()){
if (UseRandom::rnd()<.5){
CLNode= Node->randomChild();
bool inhist=CLNode->isInHistoryOf(Born);
weight*=inhist?(-2.*int(Node->children().size())):0.;
projected=true;Node->nodeME()->projectorStage(1);
weightCL=2.*int(Node->children().size());
BornCL=CLNode-> getLongestHistory_simple(false,xiQSh);
}else{
weight=2.;projected=false;Node->nodeME()->projectorStage(0);
}
}else{
weight=1.;projected=false;Node->nodeME()->projectorStage(0);
}
- if (Node->treefactory()->onlymulti()!=-1&&
- Node->treefactory()->onlymulti()!=int(Node->xcomb()->mePartonData().size()-Node->nodeME()->projectorStage()))
+ if (treefactory()->onlymulti()!=-1&&
+ treefactory()->onlymulti()!=int(Node->xcomb()->mePartonData().size()-Node->nodeME()->projectorStage()))
return 0.;
- if(!Node->allAbove(Node->mergePt()-0.1*GeV))weight=0.;
- if(CLNode&&!CLNode->allAbove(Node->mergePt()-0.1*GeV))weightCL=0.;
+ if(!Node->allAbove(mergePt()-0.1*GeV))weight=0.;
+ if(CLNode&&!CLNode->allAbove(mergePt()-0.1*GeV))weightCL=0.;
if (!Born->xcomb()->willPassCuts()){
if (Born->parent()) {
//cout<<"\n parent not pass";
}
return 0.;
}
double res=0.;
bool maxMulti=Node->xcomb()->meMomenta().size()-2 == theMaxLegsLO;
- double gamma=Node->treefactory()->gamma();
+
if(weight!=0.){
Energy startscale=CKKW_StartScale(Born);
Energy projectedscale=startscale;
fillHistory( startscale, Born, Node,fast);
if (!fillProjector(projectedscale))weight=0.;
Node->runningPt(projectedscale);
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.&&weightCL==0.)return 0.;
- //Node->vetoPt((projected&&maxMulti)?Node->mergePt():history.back().scale);
+ //Node->vetoPt((projected&&maxMulti)?mergePt():history.back().scale);
- res+=weight*matrixElementWeight(startscale,Node,(!maxMulti&&!projected)?gamma:1.);
+ res+=weight*matrixElementWeight(startscale,Node,(!maxMulti&&!projected)?theGamma:1.);
}
//cout<<"\n"<<((!maxMulti&&!projected)?gamma:1.)
//<<history.back().weight<<" "<<alphaReweight()<<" "<<pdfReweight();;
- if(CLNode&&gamma!=1.){
+ if(CLNode&&theGamma!=1.){
//assert(false);
Energy startscale=CKKW_StartScale(BornCL);
Energy projectedscale=startscale;
fillHistory( startscale, BornCL, CLNode,fast);
if (!fillProjector(projectedscale))return 0.;
Node->runningPt(projectedscale);
weightCL*=history.back().weight*alphaReweight()*pdfReweight();
double diff=0;
Node->nodeME()->factory()->setAlphaParameter(1.);
diff-=weightCL*CLNode->dipol()->dSigHatDR(sqr(startscale*xiFacME))/nanobarn;
- Node->nodeME()->factory()->setAlphaParameter(gamma);
+ Node->nodeME()->factory()->setAlphaParameter(theGamma);
string alp=(CLNode->dipol()->aboveAlpha()?"Above":"Below");
diff+=weightCL*CLNode->dipol()->dSigHatDR(sqr(startscale*xiFacME))/nanobarn;
Node->nodeME()->factory()->setAlphaParameter(1.);
//cout<<"\n diff "<<diff;
//cout<<"\ndiff "<<diff<<" "<<history.back().weight<<" "<<alphaReweight()<<" "<<alp;
res+=diff;
}
return res;
}
-
double Merger::reweightCKKWBornStandard(NPtr Node,bool fast){
if (!Node->children().empty()) {
PVector particles;
for (size_t i=0;i<Node->nodeME()->mePartonData().size();i++){
Ptr<ThePEG::Particle>::ptr p =Node->nodeME()->mePartonData()[i]->produceParticle(Node->nodeME()->lastMEMomenta()[i]);
particles.push_back(p);
}
if(!matrixElementRegion(particles, Node->children()[0]->dipol()->lastPt(), MergingScale))return 0.;
}
NPtr Born= Node-> getLongestHistory_simple(true,xiQSh);//Longest irreführend
if( Born!= Node){
if (UseRandom::rnd()<.5){
weight=-2.;projected=true;Node->nodeME()->projectorStage(1);
}else{
weight=2.;projected=false;Node->nodeME()->projectorStage(0);
}
}else{
weight=1.;projected=false;Node->nodeME()->projectorStage(0);
}
- if (Node->treefactory()->onlymulti()!=-1&&
- Node->treefactory()->onlymulti()!=int(Node->xcomb()->mePartonData().size()-Node->nodeME()->projectorStage()))
+ if (treefactory()->onlymulti()!=-1&&
+ treefactory()->onlymulti()!=int(Node->xcomb()->mePartonData().size()-Node->nodeME()->projectorStage()))
return 0.;
- assert(Node->allAbove(Node->mergePt()-0.1*GeV));
+ assert(Node->allAbove(mergePt()-0.1*GeV));
if (!Born->xcomb()->willPassCuts()){
if (Born->parent()) {
//cout<<"\n parent not pass";
}
return 0.;
}
Energy startscale=CKKW_StartScale(Born);
Energy projectedscale=startscale;
fillHistory( startscale, Born, Node,fast);
if (!fillProjector(projectedscale))return 0.;
Node->runningPt(projectedscale);
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.)return 0.;
bool maxMulti=Node->xcomb()->meMomenta().size()-2 == theMaxLegsLO;
- Node->vetoPt((projected&&maxMulti)?Node->mergePt():history.back().scale);
+ Node->vetoPt((projected&&maxMulti)?mergePt():history.back().scale);
return weight*matrixElementWeight(startscale,Node,1.);
}
double Merger::reweightCKKWVirtualStandard(NPtr Node,bool fast){
NPtr Born= Node-> getLongestHistory_simple(true,xiQSh);
if( Born!= Node){
if (UseRandom::rnd()<.5){
weight=-2.;projected=true;Node->nodeME()->projectorStage(1);
}else{
weight=2.;projected=false;Node->nodeME()->projectorStage(0);
}
}else{
weight=1.;projected=false;Node->nodeME()->projectorStage(0);
}
if (!Born->xcomb()->willPassCuts())return 0.;
Energy startscale=CKKW_StartScale(Born);
Energy projectedscale=startscale;
fillHistory( startscale, Born, Node,fast);
if (history.size()==1&&Node->children().size()!=0) {
cout<<"\n1-->"<<startscale/GeV<<" "<<weight;
}
if (!fillProjector(projectedscale))return 0.;
Node->runningPt(projectedscale);
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.)return 0.;
bool maxMulti=Node->xcomb()->meMomenta().size()-2 == theMaxLegsNLO;
- Node->vetoPt((projected&&maxMulti)?Node->mergePt():history.back().scale);
+ Node->vetoPt((projected&&maxMulti)?mergePt():history.back().scale);
double matrixElement=matrixElementWeight(startscale,Node);
double Bornweight=Node->nodeME()->lastBorndSigHatDR();
double unlopsweight =(-sumpdfReweightUnlops()
-sumalphaReweightUnlops()
-sumfillHistoryUnlops())
*Bornweight
*SM().alphaS()/(2.*ThePEG::Constants::pi);
//assert(unlopsweight==0.||history.size()!=1);
if (history.size()==1&&Node->children().size()!=0) {
cout<<"\n unlopsweight "<<unlopsweight<<" "<<matrixElement;
}
return weight*
as(startscale*xiRenSh)/SM().alphaS()*
(matrixElement+unlopsweight);
}
double Merger::reweightCKKWRealStandard(NPtr Node,bool fast){
- bool allAbove=Node->allAbove(Node->mergePt());
+ bool allAbove=Node->allAbove(mergePt());
if (allAbove)return reweightCKKWRealAllAbove(Node, fast);
if (UseRandom::rnd()<.5)
return 2.*reweightCKKWRealBelowSubReal( Node, fast);
return 2.*reweightCKKWRealBelowSubInt( Node, fast);
}
double Merger::reweightCKKWRealAllAbove(NPtr Node,bool fast){
NPtr Born= Node-> getLongestHistory_simple(true,xiQSh);
NPtr CLNode= Node->randomChild();
bool inhist=CLNode->isInHistoryOf(Born);
Born=CLNode-> getLongestHistory_simple(false,xiQSh);
if( Born!= CLNode){
if (UseRandom::rnd()<.5){
weight=-2.; projected=true; Node->nodeME()->projectorStage(2);}
else{
weight= 2.; projected=false; Node->nodeME()->projectorStage(1);}
}else{
weight=1.; projected=false; Node->nodeME()->projectorStage(1);
}
- if (!CLNode->allAbove(Node->mergePt()))return 0.;
+ if (!CLNode->allAbove(mergePt()))return 0.;
if (!Born->xcomb()->willPassCuts())return 0.;
Energy startscale=CKKW_StartScale(Born);
Energy projectedscale=startscale;
fillHistory( startscale, Born, CLNode,fast);
if (!fillProjector(projectedscale))return 0.;
Node->runningPt(projectedscale);
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.)return 0.;
bool maxMulti=CLNode->xcomb()->meMomenta().size()-2 == theMaxLegsNLO;
- Node->vetoPt((projected&&maxMulti)?Node->mergePt():history.back().scale);
+ Node->vetoPt((projected&&maxMulti)?mergePt():history.back().scale);
double res= weight*as(startscale*xiRenSh)/SM().alphaS()*
(double)Node->children().size()*
((inhist?matrixElementWeight(startscale,Node):0.)
+CLNode->dipol()->dSigHatDR(sqr(startscale*xiFacME))/nanobarn);
// cout<<"\nCLNode->dipol()->dSigHatDR(sqr(startscale*xiFacME))/nanobarn "<<CLNode->dipol()->dSigHatDR(sqr(startscale*xiFacME))/nanobarn;
return res;
}
double Merger::reweightCKKWRealBelowSubReal(NPtr Node,bool fast){
NPtrVec children=Node->children();
Selector<NPtr> HistNodeSel;
Energy minScale=generator()->maximumCMEnergy();
for (NPtrVec::iterator child = children.begin();
child != children.end(); child++){
if ((*child)->dipol()->lastPt()<minScale) {
minScale=(*child)->dipol()->lastPt();
HistNodeSel.insert(1.,*child);
}
}
NPtr HistNode=HistNodeSel.select(UseRandom::rnd());
NPtr Born=HistNode-> getLongestHistory_simple(false,xiQSh);
if( Born!= HistNode){
if (UseRandom::rnd()<.5){
weight=-2.; projected=true; Node->nodeME()->projectorStage(1);}
else{
weight= 2.; projected=false; Node->nodeME()->projectorStage(0);}
}else{
weight=1.; projected=false; Node->nodeME()->projectorStage(0);
}
if (!Born->xcomb()->willPassCuts())return 0.;
Energy startscale=CKKW_StartScale(Born);
Energy projectedscale=startscale;
fillHistory( startscale, Born, HistNode,fast);
if (!fillProjector(projectedscale))return 0.;
Node->runningPt(projectedscale);
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.)return 0.;
bool maxMulti=HistNode->xcomb()->meMomenta().size()-2 == theMaxLegsNLO;
- Node->vetoPt((projected&&maxMulti)?Node->mergePt():history.back().scale);
+ Node->vetoPt((projected&&maxMulti)?mergePt():history.back().scale);
double sumPS=0;
for (NPtrVec::iterator child = children.begin();
child != children.end(); child++){
- if ((*child)->allAbove(Node->mergePt()))
+ if ((*child)->allAbove(mergePt()))
sumPS-=(*child)->calcPs(startscale*xiFacME);
}
//double me=matrixElementWeight(startscale,Node);
//cout<<"\nreweightCKKWRealBelowSubReal "<<me<<" "<<sumPS<<" "<<me/sumPS;
return weight*as(startscale*xiRenSh)/SM().alphaS()*
(matrixElementWeight(startscale,Node)-sumPS);
}
double Merger::reweightCKKWRealBelowSubInt(NPtr Node,bool fast){
NPtr CLNode= Node->randomChild();
NPtr Born=CLNode-> getLongestHistory_simple(false,xiQSh);
if( Born!= CLNode){
if (UseRandom::rnd()<.5){
weight=-2.; projected=true; Node->nodeME()->projectorStage(2);}
else{
weight= 2.; projected=false; Node->nodeME()->projectorStage(1);}
}else{
weight=1.; projected=false; Node->nodeME()->projectorStage(1);
}
- if (!CLNode->allAbove(Node->mergePt()))return 0.;
+ if (!CLNode->allAbove(mergePt()))return 0.;
if (!Born->xcomb()->willPassCuts())return 0.;
Energy startscale=CKKW_StartScale(Born);
Energy projectedscale=startscale;
fillHistory( startscale, Born, CLNode,fast);
if (!fillProjector(projectedscale))return 0.;
Node->runningPt(projectedscale);
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.)return 0.;
bool maxMulti=CLNode->xcomb()->meMomenta().size()-2 == theMaxLegsNLO;
- Node->vetoPt((projected&&maxMulti)?Node->mergePt():history.back().scale);
+ Node->vetoPt((projected&&maxMulti)?mergePt():history.back().scale);
double res=0.;
if (Born==CLNode&&!CLNode->children().empty())
res=-1.*CLNode->dipol()->dSigHatDR(sqr(startscale*xiFacME))/nanobarn;
else
res=CLNode->calcPsMinusDip(startscale*xiFacME);
return weight*as(startscale*xiRenSh)/SM().alphaS()*
(double)Node->children().size()*res;
}
double Merger::matrixElementWeight(Energy startscale,NPtr Node,double diffAlpha){
double res;
NPtr DeepHead=Node;//->deepHead();
renormscale(startscale);
DeepHead->nodeME()->factory()->scaleChoice()->setXComb(DeepHead->xcomb());
DeepHead->nodeME()->setScale(sqr(startscale),sqr(startscale));
theCalculateInNode=false;
res=DeepHead->nodeME()->dSigHatDR(false,diffAlpha)/nanobarn;
renormscale(0.0*GeV);
theCalculateInNode=true;
return res;
}
double Merger::matrixElementWeightWithLoops(Energy startscale,NPtr Node,bool fast){
double res=0.;
return res;
// The deephead should be calculated here.
NPtr DeepHead=Node;//->deepHead();
renormscale(startscale);
DeepHead->nodeME()->factory()->scaleChoice()->setXComb(DeepHead->xcomb());
DeepHead->nodeME()->setScale(sqr(startscale),sqr(startscale));
theCalculateInNode=false;
DeepHead->nodeME()->doOneLoopNoBorn();
res=DeepHead->nodeME()->dSigHatDR(fast)/nanobarn;
DeepHead->nodeME()->noOneLoopNoBorn();
renormscale(0.0*GeV);
theCalculateInNode=true;
return res;
}
bool Merger::fillProjector(Energy& prerunning){
if(history.begin()->node->deepHead()->nodeME()->projectorStage() == 0){
prerunning=(history.size()==1?xiQSh:1.)*history.back().scale;
return true;
}
for (History::iterator it=history.begin();it!=history.end();it++){
if (projectorStage((*it).node)&&history.begin()->node->deepHead()->nodeME()->projectorStage() != 0){
history.begin()->node->deepHead()->xcomb()->lastProjector((*it).node->xcomb());
prerunning=(it==history.begin()?xiQSh:1.)*(*it).scale;
return true;
}
}
return false;
}
double Merger::pdfReweight(){
double res=1.;
for(int side=0;side!=2;side++){
if(history[0].node->xcomb()->mePartonData()[side]->coloured()){
History::iterator it=history.begin();
for (;it+1!=history.end();it++){
res *= pdfratio((*it).node, (*it).scale,xiFacSh*((*it).node->dipol()->lastPt()), side);
}
res*=pdfratio(history.back().node,history.back().scale ,history[0].scale*xiFacME, side);
}
}
return res;
}
double Merger::alphaReweight(){
double res=1.;
Energy Q_R=xiRenME*history[0].scale;
res *= pow(as(Q_R) / SM().alphaS(), history[0].node->nodeME()->orderInAlphaS());
res *= pow(history[0].node->deepHead()->xcomb()->eventHandler().SM().alphaEMPtr()->value(history[0].node->nodeME()->factory()->scaleChoice()->renormalizationScaleQED())/ SM().alphaEMMZ(), history[0].node->nodeME()->orderInAlphaEW());
for (History::iterator it=history.begin();(it+1)!=history.end();it++){
if ((*it).node->parent()){
Energy q_i=xiRenSh* (*it).node->dipol()->lastPt();
res *= as(q_i)/ SM().alphaS()
*(theKImproved?(1.+((3.*(67./18.-1./6.*Constants::pi*Constants::pi)-5./9.*5.)*as(q_i))/2./Constants::pi):1.);
}
}
return res;
}
void Merger::fillHistory(Energy scale, NPtr Begin, NPtr EndNode,bool fast){
history.clear();
double sudakov0_n=1.;
history.push_back(HistoryStep(Begin,sudakov0_n,scale));
scale*=xiQSh;
- if (Begin->parent()||!Begin->deepHead()->unitarized()) {
+ if (Begin->parent()||!isUnitarized) {
while (Begin->parent() && (Begin != EndNode)) {
if (!dosudakov(Begin,scale, Begin->dipol()->lastPt(),sudakov0_n,fast)){
history.push_back(HistoryStep(Begin->parent(),0.,scale));
}
scale=Begin->dipol()->lastPt();
history.push_back(HistoryStep(Begin->parent(),sudakov0_n,Begin->dipol()->lastPt()));
Begin = Begin->parent();
}
Energy notunirunning=scale;
- if (!Begin->deepHead()->unitarized()&&Begin->deepHead()->N() > Begin->deepHead()->nodeME()->lastMEMomenta().size()) {
- if (!dosudakov(Begin,notunirunning,Begin->deepHead()->mergePt(),sudakov0_n,fast)){
+ if (!isUnitarized&&N() > Begin->deepHead()->nodeME()->lastMEMomenta().size()) {
+ if (!dosudakov(Begin,notunirunning,mergePt(),sudakov0_n,fast)){
history.back().weight=0.;
}else{
history.back().weight=sudakov0_n;
}
}
}
if( history.size()==1)scale/=xiQSh;
}
double Merger::sumpdfReweightUnlops(){
double res=0.;
Energy beam1Scale=history[0].scale*xiFacME;
Energy beam2Scale=history[0].scale*xiFacME;
History::iterator it=history.begin();
for (;it!=history.end();it++){
History::iterator ittmp=it;
ittmp++;
if((*it).node->xcomb()->mePartonData()[0]->coloured()&&(*it).node->nodeME()->lastX1()>0.99)return 0.;
if((*it).node->xcomb()->mePartonData()[1]->coloured()&&(*it).node->nodeME()->lastX2()>0.99)return 0.;
if (ittmp!=history.end()){
if((*it).node->nodeME()->lastX1()<0.00001)return 0.;
if((*it).node->nodeME()->lastX2()<0.00001)return 0.;
if ((*it).node->dipol()->bornEmitter() == 0 ){
res +=pdfUnlops((*it).node->nodeME()->lastParticles().first->dataPtr(),
(*it).node->nodeME()->lastPartons().first->dataPtr(),
(*it).node->xcomb()->partonBins().first->pdf(),
beam1Scale,
((*it).node->dipol()->lastPt()),
(*it).node->nodeME()->lastX1(),
theNf,
history[0].scale);
beam1Scale=((*it).node->dipol()->lastPt())*xiFacSh;
}
if ((*it).node->dipol()->bornEmitter() == 1 ){
res +=pdfUnlops((*it).node->nodeME()->lastParticles().second->dataPtr(),
(*it).node->nodeME()->lastPartons().second->dataPtr(),
(*it).node->xcomb()->partonBins().second->pdf(),
beam2Scale,
((*it).node->dipol()->lastPt()),
(*it).node->nodeME()->lastX2(),
theNf,
history[0].scale);
beam2Scale=((*it).node->dipol()->lastPt())*xiFacSh;
}
if ((*it).node->dipol()->bornSpectator() == 0 &&(*it).node->dipol()->bornEmitter() >1){//
res +=pdfUnlops((*it).node->nodeME()->lastParticles().first->dataPtr(),
(*it).node->nodeME()->lastPartons().first->dataPtr(),
(*it).node->xcomb()->partonBins().first->pdf(),
beam1Scale,
((*it).node->dipol()->lastPt()),
(*it).node->nodeME()->lastX1(),
theNf,
history[0].scale);
//pdfratio((*it).node, beam1Scale, sqrt((*it).node->dipol()->lastPt()), 1);
beam1Scale=((*it).node->dipol()->lastPt())*xiFacSh;
}
if ((*it).node->dipol()->bornSpectator() == 1 &&(*it).node->dipol()->bornEmitter() >1){//
res +=pdfUnlops((*it).node->nodeME()->lastParticles().second->dataPtr(),
(*it).node->nodeME()->lastPartons().second->dataPtr(),
(*it).node->xcomb()->partonBins().second->pdf(),
beam2Scale,
((*it).node->dipol()->lastPt()),
(*it).node->nodeME()->lastX2(),
theNf,
history[0].scale);
//pdfratio((*it).node, beam2Scale , sqrt((*it).node->dipol()->lastPt()), 2);
beam2Scale=((*it).node->dipol()->lastPt())*xiFacSh;
}
}
}
if (history[0].node->deepHead()->xcomb()->mePartonData()[0]->coloured()){
res +=pdfUnlops(history.back().node->nodeME()->lastParticles().first->dataPtr(),
history.back().node->nodeME()->lastPartons().first->dataPtr(),
history.back().node->xcomb()->partonBins().first->pdf(),
beam1Scale,
history[0].scale*xiFacME,
(history.back()).node->nodeME()->lastX1(),
theNf,
history[0].scale);
}
if (history[0].node->deepHead()->xcomb()->mePartonData()[1]->coloured()) {
res +=pdfUnlops(history.back().node->nodeME()->lastParticles().second->dataPtr(),
history.back().node->nodeME()->lastPartons().second->dataPtr(),
history.back().node->xcomb()->partonBins().second->pdf(),
beam2Scale,
history[0].scale*xiFacME,
(history.back()).node->nodeME()->lastX2(),
theNf,
history[0].scale);
//history[0].node->deepHead()->nodeME()->pdf2(sqr(beam2Scale))/history[0].node->deepHead()->nodeME()->pdf2(sqr(history[0].scale));
}
return res;
}
double Merger::pdfUnlops(tcPDPtr particle,tcPDPtr parton,tcPDFPtr pdf,Energy running, Energy next,double x,int nlp,Energy fixedScale) {
//copied from PKOperator
double res=0.;
int number=10;//*int((max(running/next,next/running)));
for (int nr=0;nr<number;nr++){
double restmp=0;
using namespace RandomHelpers;
double r = UseRandom::rnd();
double eps = 1e-3;
pair<double,double> zw =
generate((piecewise(),
flat(0.0,x),
match(inverse(0.0,x,1.0) + inverse(1.0+eps,x,1.0))),r);
double z = zw.first;
double mapz = zw.second;
double PDFxparton=pdf->xfx(particle,parton,sqr(fixedScale),x)/x;
double CA = SM().Nc();
double CF = (SM().Nc()*SM().Nc()-1.0)/(2.*SM().Nc());
if (abs(parton->id()) < 7) {
double PDFxByzgluon=pdf->xfx(particle,getParticleData(ParticleID::g),sqr(fixedScale),x/z)*z/x;
double PDFxByzparton=pdf->xfx(particle,parton,sqr(fixedScale),x/z)*z/x;
assert(abs(parton->id()) < 7);
restmp += CF*(3./2.+2.*log(1.-x)) * PDFxparton;
if ( z > x ) {
restmp+= 0.5 * ( sqr(z) + sqr(1.-z) ) * PDFxByzgluon / z;
restmp += CF*2.*(PDFxByzparton - z*PDFxparton)/(z*(1.-z));
restmp -= CF*PDFxByzparton * (1.+z)/z;
}
}else{
assert(parton->id() == ParticleID::g);
double PDFxByzgluon=pdf->xfx(particle,getParticleData(ParticleID::g),sqr(fixedScale),x/z)*z/x;
if ( z > x ){
double factor = CF * ( 1. + sqr(1.-z) ) / sqr(z);
for ( int f = -nlp; f <= nlp; ++f ) {
if ( f == 0 )
continue;
restmp += pdf->xfx(particle,getParticleData(f),sqr(fixedScale),x/z)*z/x*factor;
}
}
restmp += ( (11./6.) * CA - (1./3.)*theNf + 2.*CA*log(1.-x) ) *PDFxparton;
if ( z > x ) {
restmp += 2. * CA * ( PDFxByzgluon - z*PDFxparton ) / (z*(1.-z));
restmp += 2.* CA *( (1.-z)/z - 1. + z*(1.-z) ) * PDFxByzgluon / z;
}
}
if (PDFxparton<1e-8)restmp= 0.;
res+=1*restmp*log(sqr(running/next))/PDFxparton*mapz;
}
return res/number;
}
double Merger::sumalphaReweightUnlops(){
double res=0.;
if (!(history[0].node->children().empty())){
res +=alphasUnlops(history[0].scale,
history[0].scale);
}
// dsig is calculated with fixed alpha_s
for (History::iterator it=history.begin();(it+1)!=history.end();it++){
assert((*it).node->parent());
res +=alphasUnlops((*it).node->dipol()->lastPt()*xiRenSh ,history[0].scale*xiRenME);
}
return res;
}
double Merger::sumfillHistoryUnlops(){
double res=0.;
for (History::iterator it = history.begin(); (it+1) != history.end();it++){
doUNLOPS((*it).node,(it == history.begin()?xiQSh:1.)*(*it).scale , (*it).node->dipol()->lastPt() , history[0].scale, res);
}
return res;
}
Ptr<MFactory>::ptr Merger::treefactory(){return theTreeFactory;}
bool Merger::reweightCKKWSingle(Ptr<MatchboxXComb>::ptr SX, double & res,bool fast) {
assert(!fast);
Ptr<StandardEventHandler>::ptr eH =
dynamic_ptr_cast<Ptr<StandardEventHandler>::ptr>(generator()->eventHandler());
//cout<<"\nfast ";//<<fast<<" HS "<<history.size()<<" stnode "<<StartingBorn<<" "<<eH->didEstimate();
if (!eH->didEstimate()||fast) {
history.clear();
StartingBorn0=NPtr();
StartingBorn1=NPtr();
StartingBorn2=NPtr();
}
theNf=5;//TODO
if (!SX) return true;
assert(SX->eventHandlerPtr());
Ptr<MatchboxMEBase>::ptr ME = dynamic_ptr_cast<Ptr<MatchboxMEBase>::ptr>(SX->matchboxME());
if (!ME) return true;
if (theFirstNodeMap.find(ME)==theFirstNodeMap.end()) {
cout<<"\nnot in map:"<<theFirstNodeMap.size();
}
NPtr Node = theFirstNodeMap[SX->matchboxME()];
Ptr<MFactory>::ptr factory = dynamic_ptr_cast<Ptr<MFactory>::ptr>(Node->nodeME()->factory());
assert(factory);
if (!Node) return true;
NPtr MENode = Node;
Ptr<AlphaEMBase>::transient_pointer alphaEM = SX->eventHandler().SM().alphaEMPtr();
assert(DSH()->hardScaleIsMuF());
xiRenME=ME->renormalizationScaleFactor();
xiFacME=ME->factorizationScaleFactor();
xiRenSh=DSH()->renormalizationScaleFactor();
xiFacSh=DSH()->factorizationScaleFactor();
xiQSh=DSH()->hardScaleFactor();
if(Node->deepHead()->subtractedReal()){
res*=reweightCKKWRealStandard(Node);
}else if(ME->oneLoopNoBorn()){
res*=reweightCKKWVirtualStandard(Node);
- }else if(Node->treefactory()->gamma()!=1.){
+ }else if(theGamma!=1.){
res*=reweightCKKWBornGamma(Node,fast);
}else{
res*=reweightCKKWBornStandard(Node,fast);
}
SX->lastCentralScale(sqr(Node->runningPt()));
SX->lastShowerScale(sqr(Node->runningPt()));
if(SX->lastProjector()){
SX->lastProjector()->lastCentralScale(sqr(Node->runningPt()));
SX->lastProjector()->lastShowerScale(sqr(Node->runningPt()));
}
renormscale(0.0*GeV);
if (res == 0.){
history.clear();
StartingBorn0=NPtr();
StartingBorn1=NPtr();
StartingBorn2=NPtr();
return false;
}
cleanup(MENode);
cleanup(Node);
DSH()->eventRecord().clear();
SX->subProcess(SubProPtr());
- Node->VetoedShower(true);
if (!fast) {
history.clear();
StartingBorn0=NPtr();
StartingBorn1=NPtr();
StartingBorn2=NPtr();
}
return true;
}
//----------------------------------Reviewed--------------------------------------//
void Merger::CKKW_PrepareSudakov(NPtr Born,Energy running){
//cleanup(Born);
tSubProPtr sub = Born->xcomb()->construct();
DSH()->resetPDFs(make_pair(Born->xcomb()->partonBins().first->pdf(),
Born->xcomb()->partonBins().second->pdf()),
Born->xcomb()->partonBins());
DSH()->setCurrentHandler();
DSH()->currentHandler()->generator()->currentEventHandler(Born->deepHead()->xcomb()->eventHandlerPtr());
DSH()->currentHandler()->remnantDecayer()->setHadronContent(Born->deepHead()->xcomb()->lastParticles());
DSH()->eventRecord().clear();
DSH()->eventRecord().prepare(sub, dynamic_ptr_cast<tStdXCombPtr>(Born->xcomb()), DSH()->pdfs(), Born->deepHead()->xcomb()->lastParticles());
DSH()->hardScales(sqr(running));
}
Energy Merger::CKKW_StartScale(NPtr Born){
Energy res=generator()->maximumCMEnergy();;
if(!Born->children().empty()){
for (size_t i=0;i<Born->nodeME()->mePartonData().size();i++){
if (!Born->nodeME()->mePartonData()[i]->coloured())continue;
for (size_t j=i+1;j<Born->nodeME()->mePartonData().size();j++){
if (i==j||!Born->nodeME()->mePartonData()[j]->coloured())continue;
res= min(res,sqrt(2.*Born->nodeME()->lastMEMomenta()[i]*Born->nodeME()->lastMEMomenta()[j]));
}
}
}else{
Born->nodeME()->factory()->scaleChoice()->setXComb(Born->xcomb());
res= sqrt(Born->nodeME()->factory()->scaleChoice()->renormalizationScale());
}
Born->nodeME()->factory()->scaleChoice()->setXComb(Born->xcomb());
res=max(res, sqrt(Born->nodeME()->factory()->scaleChoice()->renormalizationScale()));
return res;
}
double Merger::alphasUnlops( Energy next,Energy fixedScale) {
double betaZero = (11./6.)*SM().Nc() - (1./3.)*theNf;
return (betaZero*log(sqr(fixedScale/next)))+
(theKImproved?(((3.*(67./18.-1./6.*Constants::pi*Constants::pi)-5./9.*theNf))):0.);
}
double Merger::pdfratio(NPtr Born,Energy & nominator_scale, Energy denominator_scale,int side){
StdXCombPtr bXc = Born->xcomb();
assert (bXc->mePartonData()[side]->coloured());
double from,to;
if (side==0){
if (denominator_scale==nominator_scale) {
return 1.;
}
from=Born->nodeME()->pdf1(sqr(denominator_scale));
to =Born->nodeME()->pdf1(sqr( nominator_scale ));
if ( (to < 1e-8||from < 1e-8)&&(to/from>10000000.)){cout<<"\npdfratio to="<<to<<" from="<<from;return 0.;}
}
else{
if (denominator_scale==nominator_scale) {
return 1.;
}
from=Born->nodeME()->pdf2(sqr(denominator_scale));
to=Born->nodeME()->pdf2( sqr( nominator_scale ));
if ( (to < 1e-8||from < 1e-8)&&(to/from>10000000.)){cout<<"\npdfratio to="<<to<<" from="<<from;return 0.;}
}
return to/from;
}
bool Merger::dosudakov(NPtr Born,Energy running, Energy next, double& sudakov0_n,bool fast) {
CKKW_PrepareSudakov(Born, running);
for ( list<DipoleChain>::iterator chain = DSH()->eventRecord().chains().begin() ;
chain != DSH()->eventRecord().chains().end() ; chain++ ) {
for ( list<Dipole>::iterator dip = (*chain).dipoles().begin() ; dip != (*chain).dipoles().end() ; ++dip ) {
sudakov0_n*=singlesudakov( dip, next,running,make_pair(true,false), fast );
sudakov0_n*=singlesudakov( dip, next,running,make_pair(false,true), fast );
if (sudakov0_n==0.0){
cleanup(Born);
return false;
}
}
}
cleanup(Born);
return true;
}
bool Merger::doUNLOPS(NPtr Born,Energy running, Energy next,Energy fixedScale, double& UNLOPS) {
CKKW_PrepareSudakov(Born, running);
for ( list<DipoleChain>::iterator chain = DSH()->eventRecord().chains().begin() ;
chain != DSH()->eventRecord().chains().end() ; chain++ ) {
for ( list<Dipole>::iterator dip = (*chain).dipoles().begin() ; dip != (*chain).dipoles().end() ; ++dip ) {
UNLOPS+=singleUNLOPS( dip, next,running,fixedScale,make_pair(true,false) );;
UNLOPS+=singleUNLOPS( dip, next,running,fixedScale,make_pair(false,true) );
}
}
cleanup(Born);
return true;
}
bool Merger::projectorStage(NPtr Born){
Born->deepHead()->nodeME()->mePartonData().size();
return (Born->deepHead()->nodeME()->projectorStage() ==
int((Born->deepHead()->nodeME()->mePartonData().size()
- Born->nodeME()->mePartonData().size())));
}
void Merger::cleanup(NPtr Born) {
DSH()->eventRecord().clear();
if(!Born->xcomb()->subProcess())return;
ParticleVector vecfirst = Born->xcomb()->subProcess()->incoming().first->children();
for ( ParticleVector::const_iterator it = vecfirst.begin() ; it != vecfirst.end() ; it++ )
Born->xcomb()->subProcess()->incoming().first->abandonChild(*it);
ParticleVector vecsecond = Born->xcomb()->subProcess()->incoming().second->children();
for ( ParticleVector::const_iterator it = vecsecond.begin() ; it != vecsecond.end() ; it++ )
Born->xcomb()->subProcess()->incoming().second->abandonChild(*it);
Born->xcomb()->subProcess(SubProPtr());
}
double Merger::singlesudakov(list<Dipole>::iterator dip ,Energy next,Energy running,pair<bool,bool> conf ,bool fast){
double res=1.;
tPPtr emitter = dip->emitter(conf);
tPPtr spectator = dip->spectator(conf);
DipoleSplittingInfo candidate((*dip).index(conf),conf,(*dip).emitterX(conf),(*dip).spectatorX(conf),emitter,spectator);
if ( DSH()->generators().find(candidate.index()) == DSH()->generators().end() ) DSH()->getGenerators(candidate.index());
pair<GeneratorMap2::iterator,GeneratorMap2::iterator> gens = DSH()->generators().equal_range(candidate.index());
for ( GeneratorMap2::iterator gen = gens.first; gen != gens.second; ++gen ) {
if ( !(gen->first == candidate.index()) )
continue;
Energy dScale = gen->second->splittingKinematics()->dipoleScale(emitter->momentum(),spectator->momentum());
assert (! isnan(dScale/GeV ) );
candidate.scale(dScale);
candidate.continuesEvolving();
Energy ptMax=(*gen).second->splittingKinematics()->ptMax(candidate.scale(),candidate.emitterX(), candidate.spectatorX(),
candidate.index(),*gen->second->splittingKernel());
candidate.hardPt(min(running,ptMax));
if (candidate.hardPt()>next){
res*=gen->second->sudakov(candidate,next,fast);
}
}
return res;
}
double Merger::singleUNLOPS(list<Dipole>::iterator dip ,Energy next,Energy running,Energy fixedScale,pair<bool,bool> conf ){
double res=0.;
tPPtr emitter = dip->emitter(conf);
tPPtr spectator = dip->spectator(conf);
DipoleSplittingInfo candidate((*dip).index(conf),conf,(*dip).emitterX(conf),(*dip).spectatorX(conf),emitter,spectator);
if ( DSH()->generators().find(candidate.index()) == DSH()->generators().end() ) DSH()->getGenerators(candidate.index());
pair<GeneratorMap2::iterator,GeneratorMap2::iterator> gens = DSH()->generators().equal_range(candidate.index());
for ( GeneratorMap2::iterator gen = gens.first; gen != gens.second; ++gen ) {
if ( !(gen->first == candidate.index()) )
continue;
Energy dScale = gen->second->splittingKinematics()->dipoleScale(emitter->momentum(),spectator->momentum());
assert (! isnan(dScale/GeV ) );
candidate.scale(dScale);
candidate.continuesEvolving();
Energy ptMax=gen->second->splittingKinematics()->ptMax(candidate.scale(),candidate.emitterX(), candidate.spectatorX(),
candidate.index(),*gen->second->splittingKernel());
candidate.hardPt(min(running,ptMax));
if (candidate.hardPt()>next){
res+=gen->second->unlops(candidate,next,fixedScale);
}
}
return res;
}
void Merger::firstNodeMap(Ptr<MatchboxMEBase>::ptr a,NPtr b){theFirstNodeMap.insert(make_pair(a,b));};
map<Ptr<MatchboxMEBase>::ptr,NPtr> Merger::firstNodeMap(){return theFirstNodeMap;}
void Merger::setXComb(Ptr<MatchboxMEBase>::ptr me,tStdXCombPtr xc,int st){
theFirstNodeMap[me]->setXComb(xc, st);
}
void Merger::setKinematics(Ptr<MatchboxMEBase>::ptr me){
theFirstNodeMap[me]->setKinematics();
}
void Merger::clearKinematics(Ptr<MatchboxMEBase>::ptr me){
theFirstNodeMap[me]->clearKinematics();
}
bool Merger::generateKinematics(Ptr<MatchboxMEBase>::ptr me,const double * r){
theFirstNodeMap[me]->firstgenerateKinematics(r, 0,me->lastXCombPtr()->lastSHat());
if (theFirstNodeMap[me]->cutStage()==0 ){
bool inAlphaPS=false;
NPtrVec children=theFirstNodeMap[me]->children();
for (NPtrVec::iterator child = children.begin();
child != children.end(); child++){
- double gamma=treefactory()->gamma();;
- treefactory()->setAlphaParameter(gamma);
+
+ treefactory()->setAlphaParameter(theGamma);
inAlphaPS|=(*child)->dipol()->aboveAlpha();
treefactory()->setAlphaParameter(1.);
}
SafeClusterMap temp=theFirstNodeMap[me]->clusterSafe();
for(SafeClusterMap::iterator
it=temp.begin();
it!=temp.end();++it){
if (!it->second.first&&!inAlphaPS)return false;
}
}
if (theFirstNodeMap[me]->cutStage()==1 ){
SafeClusterMap temp=theFirstNodeMap[me]->clusterSafe();
for(SafeClusterMap::iterator
it=temp.begin();
it!=temp.end();++it){
if (!it->second.first && !it->second.second)return false;
}
}
return true;
}
bool Merger::calculateInNode() const{
return theCalculateInNode;
}
void Merger::fillProjectors(Ptr<MatchboxMEBase>::ptr me){
for (unsigned int i = 0; i < (theFirstNodeMap[me]->Projector()).size(); ++i) {
me->lastXCombPtr()->projectors().insert(
(theFirstNodeMap[me]->Projector())[i].first,
(theFirstNodeMap[me]->Projector())[i].second->xcomb());
}
}
pair<bool,bool> Merger::clusterSafe(Ptr<MatchboxMEBase>::ptr me,int emit,int emis,int spec){
return theFirstNodeMap[me]->clusterSafe().find(make_pair(make_pair(emit,emis),spec))->second;
}
bool Merger::matrixElementRegion(PVector particles,Energy winnerScale,Energy cutscale){
//cout<<"\nparticles s"<<particles.size()<<" "<<particles[0]<<" "<<particles[1]<<flush;
/*
if (defMERegionByJetAlg && !particles[0]->coloured()&& !particles[1]->coloured()) {
assert(false);
vector<fastjet::PseudoJet> input_particles;
for(size_t em=2; em < particles.size();em++){
input_particles.push_back(fastjet::PseudoJet(particles[em]->momentum().x()/GeV,
particles[em]->momentum().y()/GeV,
particles[em]->momentum().z()/GeV,
particles[em]->momentum().e()/GeV));
}
fastjet::JetDefinition jet_def(fastjet::ee_kt_algorithm);
fastjet::ClusterSequence clust_seq(input_particles, jet_def);
size_t n = particles.size()-2;
vector<fastjet::PseudoJet> exclusive_jets = clust_seq.exclusive_jets_ycut(ee_ycut);
return n==exclusive_jets.size();
}
if (defMERegionByJetAlg ) {
assert(false);
size_t noncol=0;
vector<fastjet::PseudoJet> input_particles;
for(size_t em=2; em < particles.size();em++){
if (particles[em]->coloured())
input_particles.push_back(fastjet::PseudoJet(particles[em]->momentum().x()/GeV,
particles[em]->momentum().y()/GeV,
particles[em]->momentum().z()/GeV,
particles[em]->momentum().e()/GeV));
else
noncol++;
}
double Rparam = 1.0;
fastjet::Strategy strategy = fastjet::Best;
fastjet::RecombinationScheme recomb_scheme = fastjet::E_scheme;
fastjet::JetDefinition jet_def(fastjet::kt_algorithm, Rparam, recomb_scheme, strategy);
fastjet::ClusterSequence clust_seq(input_particles, jet_def);
size_t n = particles.size()-2-noncol;
vector<fastjet::PseudoJet> exclusive_jets = clust_seq.exclusive_jets(pp_dcut);
// cout<<"\nn= "<<n<<" jets= "<<exclusive_jets.size();
//for (size_t i=0; i<exclusive_jets.size(); i++) {
//cout<<"\nn= "<<n<<" jetspt= "<<exclusive_jets[i].perp();
//}
return n==exclusive_jets.size();
}
*/
//assert(false);
Energy ptx=1000000.*GeV;
bool foundwinnerpt=false;
//FF
for(size_t em=2; em < particles.size();em++){
if (!particles[em]->coloured()) continue;
for(size_t emm=2; emm < particles.size();emm++){
if (!particles[emm]->coloured()) continue;
if (em==emm) continue;
for(size_t spe=2; spe < particles.size();spe++){
if (!particles[spe]->coloured()) continue;
if (em==spe||emm==spe) continue;
if (!(particles[em]->id()==-particles[emm]->id()||particles[emm]->id()>6))continue;
// assert(false);
Lorentz5Momentum emittermom = particles[em]->momentum();
Lorentz5Momentum emissionmom = particles[emm]->momentum();
Lorentz5Momentum spectatormom = particles[spe]->momentum();
Energy scale = sqrt(2.*(emissionmom*emittermom + emissionmom*spectatormom + emittermom*spectatormom));
double y = emissionmom*emittermom / (emissionmom*emittermom + emissionmom*spectatormom + emittermom*spectatormom);
double z = emittermom*spectatormom / (emittermom*spectatormom + emissionmom*spectatormom);
if (abs(scale * sqrt(y*z*(1.-z))-winnerScale)<0.0001*GeV) {
foundwinnerpt=true;
}
// if(scale * sqrt(y*z*(1.-z))<optVeto&&winnerScale>optVeto)cout<<"\nFF "<<(scale * sqrt(y*z*(1.-z))/GeV);
ptx =min(ptx,scale * sqrt(y*z*(1.-z)));
}
}
}
//FI
for(size_t spe=0; spe < 2;spe++){
if (!particles[spe]->coloured()) continue;
for(size_t em=2; em < particles.size();em++){
if (!particles[em]->coloured()) continue;
for(size_t emm=2; emm < particles.size();emm++){
if (!particles[emm]->coloured()) continue;
if (em==emm) continue;
if (!(particles[em]->id()==-particles[emm]->id()||particles[emm]->id()>6))continue;
// assert(false);
Lorentz5Momentum emittermom = particles[em]->momentum();
Lorentz5Momentum emissionmom = particles[emm]->momentum();
Lorentz5Momentum spectatormom = particles[spe]->momentum();
Lorentz5Momentum Emsp=emittermom+emissionmom-spectatormom;
Energy scale = sqrt(-1.*Emsp*Emsp);
double x = (- emissionmom*emittermom + emissionmom*spectatormom + emittermom*spectatormom) /(emittermom*spectatormom + emissionmom*spectatormom);
double z = emittermom*spectatormom / (emittermom*spectatormom + emissionmom*spectatormom);
// if(scale * sqrt(z*(1.-z)*(1.-x)/x)<optVeto&&winnerScale>optVeto)cout<<"\nFI "<<(scale * sqrt(z*(1.-z)*(1.-x)/x)/GeV);
if (abs(scale *sqrt(z*(1.-z)*(1.-x)/x)-winnerScale)<0.0001*GeV) {
foundwinnerpt=true;
}
if((z*(1.-z)*(1.-x)/x)>0.)
ptx =min(ptx,scale * sqrt(z*(1.-z)*(1.-x)/x));
}
}
}
//IF
for(size_t em=0; em < 2;em++){
if (!particles[em]->coloured()) continue;
for(size_t emm=2; emm < particles.size();emm++){
if (!particles[emm]->coloured()) continue;
for(size_t spe=2; spe < particles.size();spe++){
if (!particles[spe]->coloured()) continue;
if (emm==spe) continue;
if (!(particles[em]->id()>6|| particles[em]->id()==particles[emm]->id() ||particles[emm]->id()>6))continue;
// assert(false);
Lorentz5Momentum emittermom = particles[em]->momentum();
Lorentz5Momentum emissionmom = particles[emm]->momentum();
Lorentz5Momentum spectatormom = particles[spe]->momentum();
Lorentz5Momentum Emsp=-emittermom+emissionmom+spectatormom;
Energy scale = sqrt(-1.*Emsp*Emsp);
double x = (-emissionmom*spectatormom + emittermom*spectatormom + emittermom*emissionmom) /(emittermom*emissionmom + emittermom*spectatormom);
double u = emittermom*emissionmom / (emittermom*emissionmom + emittermom*spectatormom);
if (abs(scale *sqrt(u*(1.-u)*(1.-x)/x)-winnerScale)<0.0001*GeV) {
foundwinnerpt=true;
}
if((u*(1.-u)*(1.-x)/x)>0.)
ptx =min(ptx,scale * sqrt(u*(1.-u)*(1.-x)/x));
}
}
}
//II
for(size_t em=0; em < 2;em++){
if (!particles[em]->coloured()) continue;
for(size_t spe=0; spe < 2;spe++){
if (!particles[spe]->coloured()) continue;
for(size_t emm=2; emm < particles.size();emm++){
if (!particles[emm]->coloured()) continue;
if (em==spe) continue;
if (!(particles[em]->id()>6||
particles[em]->id()==particles[emm]->id() ||
particles[emm]->id()>6))continue;
//assert(false);
Lorentz5Momentum emittermom = particles[em]->momentum();
Lorentz5Momentum emissionmom = particles[emm]->momentum();
Lorentz5Momentum spectatormom = particles[spe]->momentum();
Lorentz5Momentum Emsp=emittermom-emissionmom+spectatormom;
Energy scale = sqrt(Emsp*Emsp);
double x = (emittermom*spectatormom - emittermom*emissionmom - spectatormom*emissionmom)/(emittermom*spectatormom);
double v = (emittermom*emissionmom)/(emittermom*spectatormom);
if (abs(scale * sqrt(v*(1.-x-v)/x)-winnerScale)<0.0001*GeV) {
foundwinnerpt=true;
}
if ((v*(1.-x-v)/x)>0.)
ptx =min(ptx, scale * sqrt(v*(1.-x-v)/x));
}
}
}
// cout<<"\n"<<cutscale/GeV<< " "<<foundwinnerpt<<" "<<ptx/GeV;
assert(foundwinnerpt);
return (ptx>cutscale) ;
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void Merger::persistentOutput(PersistentOStream & os) const {
os <<xiRenME <<xiFacME <<xiRenSh
<<xiFacSh <<xiQSh <<theNf
<<minusL <<Unlopsweights
<< theKImproved << ounit(MergingScale,GeV)
<<theMaxLegsLO <<theMaxLegsNLO
<<theDipoleShowerHandler<<theTreeFactory<<theFirstNodeMap ;
}
void Merger::persistentInput(PersistentIStream & is, int) {
is >>xiRenME >>xiFacME >>xiRenSh
>>xiFacSh >>xiQSh >>theNf
>>minusL >>Unlopsweights
>> theKImproved >> iunit(MergingScale,GeV)
>>theMaxLegsLO >>theMaxLegsNLO
>>theDipoleShowerHandler>>theTreeFactory>>theFirstNodeMap;
}
// *** 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<Merger, MergerBase> describeHerwigMerger("Herwig::Merger", "HwDipoleShower.so");
//ClassDescription<Merger> Merger::initMerger;
// Definition of the static class description member.
void Merger::Init() {
static ClassDocumentation<Merger> documentation
("Merger.");
static Reference<Merger,DipoleShowerHandler> interfaceShowerHandler
("DipoleShowerHandler",
"",
&Merger::theDipoleShowerHandler, false, false, true, true, false);
static Switch<Merger,bool>
interfaceminusL ("minusL","",&Merger::minusL, false, false, false);
static SwitchOption interfaceminusLYes
(interfaceminusL,"Yes","",true);
static SwitchOption interfaceminusLNo
(interfaceminusL,"No","",false);
static Switch<Merger,bool>
interfaceUnlopsweights ("Unlopsweights","",&Merger::Unlopsweights, false, false, false);
static SwitchOption interfaceUnlopsweightsYes
(interfaceUnlopsweights,"Yes","",true);
static SwitchOption interfaceUnlopsweightsNo
(interfaceUnlopsweights,"No","",false);
static Switch<Merger,bool>
interfacetheKImproved ("KImproved","",&Merger::theKImproved, false, false, false);
static SwitchOption interfacetheKImprovedYes
(interfacetheKImproved,"Yes","",true);
static SwitchOption interfacetheKImprovedNo
(interfacetheKImproved,"No","",false);
static Parameter<Merger,Energy> interfaceMergerScale
("MergingScale",
"The MergingScale.",
&Merger::MergingScale, GeV, 20.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Parameter<Merger,unsigned int> interfaceMaxLegsLO
("MaxLegsLO",
".",
&Merger::theMaxLegsLO, 0, 0, 0,
false, false, Interface::lowerlim);
static Parameter<Merger,unsigned int> interfaceMaxLegsNLO
("MaxLegsNLO",
".",
&Merger::theMaxLegsNLO, 0, 0, 0,
false, false, Interface::lowerlim);
static Reference<Merger,MFactory> interfaceMergerHelper
("MFactory",
"",
&Merger::theTreeFactory, false, false, true, true, false);
+
+
+ static Parameter<Merger,double> interfacedcut
+ ("pp_dcut",
+ "The MergingScale.",
+ &Merger::pp_dcut, 5.0, 0.0, 0,
+ false, false, Interface::lowerlim);
+
+ static Parameter<Merger,double> interfaceycut
+ ("ee_ycut",
+ "The MergingScale.",
+ &Merger::ee_ycut, 0.2, 0.0, 0,
+ false, false, Interface::lowerlim);
+
+ static Parameter<Merger,double> interfacegamma
+ ("gamma",
+ "gamma parameter.",
+ &Merger::theGamma, 1.0, 0.0, 0,
+ false, false, Interface::lowerlim);
+
+
+ static Reference<Merger,JetFinder> interfaceMergingJetFinder
+ ("MergingJetFinder",
+ "A reference to the jet finder used in Merging.",
+ &Merger::theMergingJetFinder, false, false, true, false, false);
+
+
+
+ static Reference<Merger,ColourBasis> interfaceLargeNBasis
+ ("LargeNBasis",
+ "Set the large-N colour basis implementation.",
+ &Merger::theLargeNBasis, false, false, true, true, false);
+
+
+
+
+
+
+ static Switch<Merger,bool>
+ interfacedefMERegionByJetAlg ("MERegionByJetAlg","",&Merger::defMERegionByJetAlg, false, false, false);
+
+ static SwitchOption interfacedefMERegionByJetAlgYes
+ (interfacedefMERegionByJetAlg,"Yes","",true);
+ static SwitchOption interfacedefMERegionByJetAlgNo
+ (interfacedefMERegionByJetAlg,"No","",false);
+
+
+
+
+
+ static Parameter<Merger, Energy> interfaceIRSafePT("IRSafePT", "Set the pt for which a matrix element should be treated as IR-safe.",
+ &Merger::theIRSafePT, GeV, 0.0 * GeV, ZERO, Constants::MaxEnergy, true, false, Interface::limited);
+ interfaceIRSafePT.setHasDefault(false);
+
+
+
+ static Parameter<Merger, double> interfacemergePtsmearing("MergingScaleSmearing", "Set the percentage the merging pt should be smeared.",
+ &Merger::theSmearing, 0., 0.,
+ 0.0, 0.5, true, false, Interface::limited);
+
+
+
+ static Parameter<Merger, int> interfacechooseHistory("chooseHistory", "different ways of choosing the history", &Merger::theChooseHistory, 3, -1, 0,
+ false, false, Interface::lowerlim);
+
+
+
+ static Parameter<Merger, int> interfaceadditionalN("additionalN", "Set the number of additional jets to consider.", &Merger::theN, 0, 0,
+ 0, false, false, Interface::lowerlim);
+
+ static Parameter<Merger, int> interfacevirtualM("virtualM",
+ "Set the number of virtual corrections to consider. -1 is default for no virtual correction.", &Merger::theM, -1, -1, 0, false, false,
+ Interface::lowerlim);
+
+
}
-
+
diff --git a/DipoleShower/Merging2/Merger.h b/DipoleShower/Merging2/Merger.h
--- a/DipoleShower/Merging2/Merger.h
+++ b/DipoleShower/Merging2/Merger.h
@@ -1,328 +1,419 @@
// -*- C++ -*-
//
// Merger.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_Merger_H
#define HERWIG_Merger_H
//
// This is the declaration of the Merger class.
//
#include "MFactory.fh"
-#include "Node.h"
+#include "Node.fh"
#include "ThePEG/Handlers/HandlerBase.h"
#include "ThePEG/Vectors/SpinOneLorentzRotation.h"
#include "Herwig/DipoleShower/DipoleShowerHandler.h"
#include "Herwig/DipoleShower/Base/DipoleSplittingGenerator.h"
#include "Herwig/MatrixElement/Matchbox/Base/MergerBase.h"
#include "ThePEG/Cuts/JetFinder.h"
#include "ThePEG/Cuts/Cuts.h"
namespace Herwig {
using namespace ThePEG;
typedef Ptr<Node>::ptr NPtr;
typedef vector<Ptr<Node>::ptr> NPtrVec;
struct HistoryStep {
HistoryStep(){}
HistoryStep(NPtr cn,double w ,Energy sc){
node=cn;
weight=w;
scale=sc;
}
NPtr node;
double weight;
Energy scale;
};
typedef vector< HistoryStep > History;
typedef multimap<DipoleIndex,Ptr<DipoleSplittingGenerator>::ptr> GeneratorMap2;
/**
* \ingroup DipoleShower
* \author Johannes Bellm
*
* \brief Merger handles the Merger ....... //TODO .
*
* @see \ref MergerInterfaces "The interfaces"
* defined for Merger.
*/
class Merger: public MergerBase {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
Merger();
/**
* The destructor.
*/
virtual ~Merger();
//@}
public:
double singlesudakov(list<Dipole>::iterator,Energy,Energy,pair<bool,bool>,bool fast=false);
bool dosudakov(NPtr Born,Energy running, Energy next, double& sudakov0_n,bool fast=false);
//bool dosudakovold(NPtr, Energy, Energy, double&);
//bool sudakov(NPtr Born, Energy running, Energy next);
void cleanup(NPtr);
bool projectorStage(NPtr);
Energy CKKW_StartScale(NPtr);
void CKKW_PrepareSudakov(NPtr,Energy);
double matrixElementWeight(Energy startscale,NPtr,double diffalpha=1.);
double matrixElementWeightWithLoops(Energy startscale,NPtr,bool);
bool fillProjector(Energy&);
void fillHistory(Energy, NPtr, NPtr ,bool fast=false);
double sumpdfReweightUnlops();
double sumalphaReweightUnlops();
double pdfratio(NPtr,Energy&,Energy,int);
double pdfReweight();
double alphaReweight();
double sumfillHistoryUnlops();
double alphasUnlops( Energy next,Energy fixedScale);
double pdfUnlops(tcPDPtr,tcPDPtr,tcPDFPtr,Energy,Energy,double,int,Energy);
double singleUNLOPS(list<Dipole>::iterator,Energy,Energy,Energy,pair<bool,bool>);
bool doUNLOPS(NPtr Born,Energy running, Energy next,Energy fixedScale, double& UNLOPS);
bool reweightCKKWSingle(Ptr<MatchboxXComb>::ptr SX, double & res,bool fast=false) ;
// double reweightCKKWBorn(NPtr Node,bool fast=false);
//double reweightCKKWBorn2(NPtr Node,bool fast=false);
// double reweightCKKWBorn3(NPtr Node,bool fast=false);
double reweightCKKWBornStandard(NPtr Node,bool fast=false);
double reweightCKKWBornGamma(NPtr Node,bool fast=false);
double reweightCKKWVirtualStandard(NPtr Node,bool fast=false);
double reweightCKKWRealStandard(NPtr Node,bool fast=false);
double reweightCKKWRealAllAbove(NPtr Node,bool fast=false);
double reweightCKKWRealBelowSubReal(NPtr Node,bool fast=false);
double reweightCKKWRealBelowSubInt(NPtr Node,bool fast=false);
//double reweightCKKWVirt(NPtr Node);
//double reweightCKKWReal(NPtr Node);
double as(Energy q){return DSH()->as(q);}
Ptr<DipoleShowerHandler>::ptr DSH(){return theDipoleShowerHandler;}
Energy mergingScale()const {return MergingScale;}
size_t maxLegsLO() const {return theMaxLegsLO;}
size_t maxLegsNLO()const {return theMaxLegsNLO;}
Ptr<MFactory>::ptr treefactory();
map<Ptr<MatchboxMEBase>::ptr,NPtr> firstNodeMap();
void setXComb(Ptr<MatchboxMEBase>::ptr,tStdXCombPtr,int);
void setKinematics(Ptr<MatchboxMEBase>::ptr);
void clearKinematics(Ptr<MatchboxMEBase>::ptr);
bool generateKinematics(Ptr<MatchboxMEBase>::ptr,const double *);
bool calculateInNode() const;
void fillProjectors(Ptr<MatchboxMEBase>::ptr);
pair<bool,bool> clusterSafe(Ptr<MatchboxMEBase>::ptr,int,int,int);
bool theCalculateInNode;
bool matrixElementRegion(PVector particles,Energy winnerScale=0.*GeV,Energy cutscale=0.*GeV);
void renormscale(Energy x) {
therenormscale = x;
}
Energy renormscale() {
return therenormscale;
}
+ Energy mergePt() {
+ return theMergePt;
+ }
+
+ void mergePt(Energy x) {
+ theMergePt = x;
+ }
+
+ Energy centralMergePt() {
+ return theCentralMergePt;
+ }
+
+ void centralMergePt(Energy x) {
+ theCentralMergePt = x;
+ }
+
+
void firstNodeMap(Ptr<MatchboxMEBase>::ptr,NPtr);
+ void smeareMergePt(){theMergePt=centralMergePt()*(1.+0.*(-1.+2.*UseRandom::rnd())*smear());}
+ double gamma()const{return theGamma;}
+
+
+
+
+ double smear()const{return theSmearing;}
+
+
+ bool MERegionByJetAlg(){return defMERegionByJetAlg;}
+
+
+ Ptr<ColourBasis>::ptr largeNBasis(){return theLargeNBasis;}
+
+
+ void largeNBasis(Ptr<ColourBasis>::ptr x){theLargeNBasis=x;}
+
+ int M()const{return theM;}
+
+ int N()const{return theM;}
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();
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;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* Variations
*/
double xiRenME;
double xiFacME;
double xiRenSh;
double xiFacSh;
double xiQSh;
/**
*
**/
double theNf;
bool minusL;
bool Unlopsweights;
bool theKImproved;
Energy MergingScale;
unsigned int theMaxLegsLO;
unsigned int theMaxLegsNLO;
double projectorWeight2Born0;
double projectorWeight1Born1;
double projectorWeight2Born1;
double projectorWeight0Born2;
double projectorWeight1Born2;
double projectorWeight1Real1;
double projectorWeight2Real1;
double projectorWeight0Real2;//Theta<
double projectorWeight1Real2;
double projectorWeight2Real2;
double projectorWeight0;
double projectorWeight1;
NPtr CalcBorn;
bool projected;
double weight,weightCB;
History history;
NPtr StartingBorn0;
NPtr StartingBorn1;
NPtr StartingBorn2;
NPtr CalcBorn0;
NPtr CalcBorn1;
NPtr CalcBorn2;
Energy therenormscale;
+ /**
+ * If any clustering below the CutStage has a lower pT than the MergePt
+ * the phase space point has to be thrown away.
+ */
+
+ Energy theIRSafePT;
+
+
+ double theGamma;
+
+
+
+
+
+ int theChooseHistory;
+
+
+
+ Energy theMergePt;
+
+ double ee_ycut;
+
+ double pp_dcut;
+
+
+ Energy theCentralMergePt;
+
+ double theSmearing;
+
+ int theN0;
+
+ int theOnlyN;
+
+ int theN;
+
+ int theM;
+
+ bool isUnitarized;
+
+ bool isNLOUnitarized;
+
+
+ bool defMERegionByJetAlg;
+
+
+
+ Ptr<JetFinder>::ptr theMergingJetFinder;
+
+ Ptr<ColourBasis>::ptr theLargeNBasis;
+
+
+
/**
* The mean of the Gaussian distribution for
* the intrinsic pt of sea partons.
*/
Ptr<DipoleShowerHandler>::ptr theDipoleShowerHandler;
Ptr<MFactory>::ptr theTreeFactory;
map<Ptr<MatchboxMEBase>::ptr,NPtr> theFirstNodeMap;
+
+
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<Merger> initMerger;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
Merger & operator=(const Merger &);
};
}
#endif /* HERWIG_Merger_H */
diff --git a/DipoleShower/Merging2/Node.cc b/DipoleShower/Merging2/Node.cc
--- a/DipoleShower/Merging2/Node.cc
+++ b/DipoleShower/Merging2/Node.cc
@@ -1,1212 +1,603 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the Node class.
//
#include "Node.h"
#include "MFactory.h"
#include "Merger.h"
#include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
#include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
#include "Herwig/Shower/ShowerHandler.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Handlers/StdXCombGroup.h"
#include "Herwig/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig/MatrixElement/Matchbox/Phasespace/TildeKinematics.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
Node::Node() {
- theclusteredto=0;
}
bool NodeDebug=false;
-Node::Node(Ptr<MatchboxMEBase>::ptr nodeME, int deepprostage, int cutstage, bool nFOH) {
+Node::Node(Ptr<MatchboxMEBase>::ptr nodeME, int deepprostage, int cutstage) {
+ nodeME->maxMultCKKW(1);
+ nodeME->minMultCKKW(0);
+ theDeepHead = this;
thenodeMEPtr = nodeME;
theDeepProStage = deepprostage;
theCutStage = cutstage;
- theNeedFullOrderedHistory = nFOH;
- needsVetoedShower=false;
theSubtractedReal=false;
theVirtualContribution=false;
- thefiniteDipoles=false;
- thesubCorrection=false;
- theclusteredto=0;
- theOnlyN=-1;
- theNumberOfSplittings=0;
}
Node::Node(Ptr<Node>::ptr deephead, Ptr<Node>::ptr head, Ptr<SubtractionDipole>::ptr dipol, Ptr<MatchboxMEBase>::ptr nodeME,
int deepprostage, int cutstage) {
theDeepHead = deephead;
theparent = head;
thedipol = dipol;
thenodeMEPtr = nodeME;
theDeepProStage = deepprostage;
theCutStage = cutstage;
- theNeedFullOrderedHistory = deephead->needFullOrderedHistory();
- needsVetoedShower=false;
theSubtractedReal=false;
theVirtualContribution=false;
- thefiniteDipoles=false;
- thesubCorrection=false;
- theclusteredto=0;
- theOnlyN=-1;
- theNumberOfSplittings=0;
}
Node::~Node() { }
-
-
-Energy Node::theMergePt(2.*GeV);
-Energy Node::theCentralMergePt(2.*GeV);
-double Node::smearing(0.) ;
-Energy Node::theIRsafePt(1000000.*GeV);
-double Node::theIRCsafeRatio(1000);
-bool Node::isUnitarized(true);
-bool Node::isNLOUnitarized(true);
-int Node::theChooseHistory(3);
-
-unsigned int Node::theN0(0);
-int Node::theOnlyN(-1);
-unsigned int Node::theN(0);
-unsigned int Node::theM(0);
-
-
-
Ptr<SubtractionDipole>::ptr Node::dipol() const {
return thedipol;
}
/** returns the matrix element pointer */
const Ptr<MatchboxMEBase>::ptr Node::nodeME() const {
return thenodeMEPtr;
}
/** access the matrix element pointer */
Ptr<MatchboxMEBase>::ptr Node::nodeME() {
return thenodeMEPtr;
}
Ptr<Node>::ptr Node::randomChild() {
return thechildren[(int)(UseRandom::rnd() * thechildren.size())];
}
bool Node::allAbove(Energy pt){
for (vector<Ptr<Node>::ptr>::iterator it = thechildren.begin(); it != thechildren.end(); it++) {
if((*it)->dipol()->lastPt()<pt)return false;
}
return true;
}
bool Node::isInHistoryOf(Ptr<Node>::ptr other){
while (other->parent()) {
if(other==this)return true;
other=other->parent();
}
return false;
}
void Node::flushCaches() {
this->theProjectors.clear();
for ( unsigned int i = 0 ; i < thechildren.size() ; ++i ) {
thechildren[i]->dipol()->underlyingBornME()->flushCaches();
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator r = thechildren[i]->dipol()->reweights().begin() ; r != thechildren[i]->dipol()->reweights().end() ;
++r ) {
(**r).flushCaches();
}
if ( thechildren[i]->xcomb() ) thechildren[i]->xcomb()->clean();
thechildren[i]->nodeME()->flushCaches();
thechildren[i]->flushCaches();
}
}
-
-double Node::smear(){
- return treefactory()->smear();
-}
-
-
-
void Node::setKinematics() {
for ( unsigned int i = 0 ; i < thechildren.size() ; ++i ) {
thechildren[i]->dipol()->setXComb(thechildren[i]->xcomb());
thechildren[i]->dipol()->setKinematics();
thechildren[i]->nodeME()->setKinematics();
thechildren[i]->setKinematics();
}
}
void Node::clearKinematics() {
for ( unsigned int i = 0 ; i < thechildren.size() ; ++i ) {
thechildren[i]->dipol()->setXComb(thechildren[i]->xcomb());
thechildren[i]->nodeME()->clearKinematics();
thechildren[i]->dipol()->clearKinematics();
thechildren[i]->clearKinematics();
}
}
bool Node::generateKinematics(const double *r, int stage, Energy2 shat) {
isOrdered = false;
isthissafe=true;
for ( unsigned int i = 0 ; i < thechildren.size() ; ++i ) {
thechildren[i]->dipol()->setXComb(thechildren[i]->xcomb());
if ( !thechildren[i]->dipol()->generateKinematics(r) ) cout << "stop";
-
- // if(!thechildren[i]->xcomb()->willPassCuts()){
- // return false;
- // }
if ( dipol()->lastPt() < thechildren[i]->dipol()->lastPt() && parent()->ordered() ) {
isOrdered = true;
deepHead()->setOrderedSteps(stage + 1);
}
-
thechildren[i]->generateKinematics(r, stage + 1, thechildren[i]->xcomb()->lastSHat());
- isthissafe = (isthissafe && thechildren[i]->dipol()->lastPt() >=(1+stage*deepHead()->treefactory()->stairfactor())*theDeepHead->mergePt());
+ isthissafe = (isthissafe && thechildren[i]->dipol()->lastPt() >=MH()->mergePt());
}
return isthissafe;
}
void Node::firstgenerateKinematics(const double *r, int stage, Energy2 shat) {
-
- //TODO: Always good?
- // if(thechildren.size()>0)
flushCaches();
//Set here the new merge Pt for the next phase space point.( Smearing!!!)
- mergePt(centralMergePt()*(1.+0.*(-1.+2.*UseRandom::rnd())*smear()));
+ MH()->smeareMergePt();
isOrdered = true;
// if we consider the hard scale to play a role in the steps we must change here to 0.
setOrderedSteps(1);
this->orderedSteps();
clustersafer.clear();
for ( unsigned int i = 0 ; i < thechildren.size() ; ++i ) {
bool ifirst = true;
bool isecond = true;
thechildren[i]->dipol()->setXComb(thechildren[i]->xcomb());
if ( !thechildren[i]->dipol()->generateKinematics(r) ) cout << "stop";
- //cout<<"\n child shat "<<thechildren[i]->xcomb()->lastSHat()/GeV2<<" "<<shat/GeV2;
- //thechildren[i]->xcomb()->lastSHat(shat);
+
isecond = thechildren[i]->generateKinematics(r, stage + 1, thechildren[i]->xcomb()->lastSHat());
-
- //TODO rethink for NLO
- ifirst = (thechildren[i]->dipol()->lastPt() >= theDeepHead->mergePt());
- // if(!thechildren[i]->xcomb()->willPassCuts()){
- // ifirst = false;
- // isecond = false;
- // }
+ ifirst = (thechildren[i]->dipol()->lastPt() >= MH()->mergePt());
+
pair<pair<int, int>, int> EmitEmisSpec = make_pair(make_pair(thechildren[i]->dipol()->realEmitter(), thechildren[i]->dipol()->realEmission()),
thechildren[i]->dipol()->realSpectator());
clustersafer.insert(make_pair(EmitEmisSpec, make_pair(ifirst, isecond)));
}
-
}
+
+
void Node::setXComb(tStdXCombPtr xc, int proStage) {
if ( !parent() ) this->xcomb(xc);
for ( unsigned int i = 0 ; i < thechildren.size() ; ++i ) {
if ( !thechildren[i]->xcomb() ) {
thechildren[i]->xcomb(thechildren[i]->dipol()->makeBornXComb(xc));
assert(thechildren[i]->xcomb());
thechildren[i]->xcomb()->head(xc);
if ( !thechildren[i]->dipol()->lastXCombPtr() ) {
thechildren[i]->dipol()->setXComb(thechildren[i]->xcomb());
}
thechildren[i]->setXComb(thechildren[i]->xcomb(), (proStage - 1));
} else {
if ( !(thechildren[i]->dipol()->lastXCombPtr()->lastScale() == thechildren[i]->xcomb()->lastScale()) ) {
thechildren[i]->dipol()->setXComb(thechildren[i]->xcomb());
}
if ( thechildren[i]->xcomb()->head() != xc ) thechildren[i]->xcomb()->head(xc);
- //Here maybe problems.
thechildren[i]->setXComb(thechildren[i]->xcomb(), (proStage - 1));
}
}
}
void Node::birth(vector<Ptr<MatchboxMEBase>::ptr> vec) {
- vector<Ptr<SubtractionDipole>::ptr> dipoles = thenodeMEPtr->getDipoles(DipoleRepository::dipoles(thenodeMEPtr->factory()->dipoleSet()), vec,true);
+ vector<Ptr<SubtractionDipole>::ptr> dipoles =
+ nodeME()->getDipoles(DipoleRepository::dipoles(
+ nodeME()->factory()->dipoleSet()), vec,true);
for ( unsigned int j = 0 ; j < dipoles.size() ; ++j ) {
dipoles[j]->doSubtraction();
- Ptr<Node>::ptr node = new_ptr(
- Node(theDeepHead, this, dipoles[j], dipoles[j]->underlyingBornME(), theDeepHead->DeepProStage(), theDeepHead->cutStage()));
+ Ptr<Node>::ptr node = new_ptr(Node(theDeepHead,this, dipoles[j],
+ dipoles[j]->underlyingBornME(),
+ theDeepHead->DeepProStage(),
+ theDeepHead->cutStage()));
thechildren.push_back(node);
}
}
vector<Ptr<Node>::ptr> Node::getNextOrderedNodes(bool normal,double hardScaleFactor) {
- //cout<<"\ngetNextOrderedNodes";
+
vector<Ptr<Node>::ptr> temp = children();
vector<Ptr<Node>::ptr> res;
for ( vector<Ptr<Node>::ptr>::const_iterator it = temp.begin() ; it != temp.end() ; ++it ) {
- if((*it)->deepHead()->mergePt()>(*it)->dipol()->lastPt()){
+ if(MH()->mergePt()>(*it)->dipol()->lastPt()){
res.clear();
// if any of the nodes is below the merging scale return empty vector
return res;
continue;
}
if (parent()&& normal){
if ( (*it)->dipol()->lastPt() < dipol()->lastPt() ){
continue;
}
}
if ( (*it)->children().size() != 0 ){
vector<Ptr<Node>::ptr> tempdown = (*it)->children();
for ( vector<Ptr<Node>::ptr>::const_iterator itd = tempdown.begin() ; itd != tempdown.end() ; ++itd ) {
if( (*itd)->dipol()->lastPt() > (*it)->dipol()->lastPt()&&(*it)->inShowerPS((*itd)->dipol()->lastPt()) ){
res.push_back(*it);
break;
}
}
}
else {
(*it)->nodeME()->factory()->scaleChoice()->setXComb((*it)->xcomb());
- //cout<<"\n"<<sqrt((*it)->nodeME()->factory()->scaleChoice()->renormalizationScale())/GeV<<" "<<(*it)->dipol()->lastPt()/GeV<<" "<<(*it)->dipol()->name();
if ( sqr(hardScaleFactor)*(*it)->nodeME()->factory()->scaleChoice()->renormalizationScale()
>= (*it)->dipol()->lastPt() * (*it)->dipol()->lastPt()&&
(*it)->inShowerPS(hardScaleFactor*sqrt((*it)->nodeME()->factory()->scaleChoice()->renormalizationScale()))) {
res.push_back(*it);
}
}
}
return res;
}
bool Node::inShowerPS(Energy hardpT){
- //return true;
- assert(deepHead()->treefactory()->largeNBasis());
+
+ assert(MH()->largeNBasis());
double x=0.;
double z_=dipol()->lastZ();
// if (dipol()->lastPt()>50*GeV) {
// return false;
// }
string type;
// II
if( dipol()->bornEmitter()<2&&dipol()->bornSpectator()<2){
type="II";
x =dipol()->bornEmitter()==0?xcomb()->lastX1():xcomb()->lastX2();
double ratio = sqr(dipol()->lastPt()/dipol()->lastDipoleScale());
double x__ = z_*(1.-z_)/(1.-z_+ratio);
double v = ratio*z_ /(1.-z_+ratio);
if (dipol()->lastPt()>(1.-x) * dipol()->lastDipoleScale()/ (2.*sqrt(x)))return false;
assert(v< 1.-x__&&x > 0. && x < 1. && v > 0.);
- //return true;
-
}
// IF
if( dipol()->bornEmitter()<2&&dipol()->bornSpectator()>=2){
type="IF";
x =dipol()->bornEmitter()==0?xcomb()->lastX1():xcomb()->lastX2();
if (dipol()->lastPt()>dipol()->lastDipoleScale()* sqrt((1.- x)/x) /2.)return false;
- //return true;
}
// FI
if( dipol()->bornEmitter()>=2&&dipol()->bornSpectator()<2){
type="FI";
double lastx=dipol()->bornSpectator()==0?xcomb()->lastX1():xcomb()->lastX2();
if (dipol()->lastPt()>dipol()->lastDipoleScale()*sqrt((1.-lastx)/lastx) /2.)return false;
}
// FF
if( dipol()->bornEmitter()>=2&&dipol()->bornSpectator()>=2){
type="FF";
if (dipol()->lastPt()>dipol()->lastDipoleScale()/2.)return false;
- // cout<<"\n "<<dipol()->bornEmitter()<<" "<<dipol()->bornSpectator()<<" "<<dipol()->lastPt()/GeV<<" "<<hardpT/GeV<<" " << sqrt(1-sqr(dipol()->lastPt()/hardpT))<<" "<<z_;
-
-
}
- //return true;
double kappa=sqr(dipol()->lastPt()/hardpT);
- //kappa=sqr(dipol()->lastPt()/dipol()->lastDipoleScale());
+
double s = sqrt(1.-kappa);
pair<double,double> zbounds= make_pair(0.5*(1.+x-(1.-x)*s),0.5*(1.+x+(1.-x)*s));
- //if(zbounds.first<z_&&z_<zbounds.second)cout<<"\n "<<type<<" "<<dipol()->lastPt()/GeV<<" "<<hardpT/GeV;
- //else cout<<"\n XX "<<type<<" "<<dipol()->lastPt()/GeV<<" "<<hardpT/GeV;;
-
- if (false&&dipol()->lastBornR()<dipol()->lastRealR()&&(zbounds.first<z_&&z_<zbounds.second)) {
- cout<<"\n"<<dipol()->name()<<" pt "<<dipol()->lastPt()/GeV<<" bR "<<dipol()->lastBornR()<<" rR "<<dipol()->lastRealR();
- }
-
- if (theChooseHistory==-1) {
-
- return (zbounds.first<z_&&z_<zbounds.second)&&dipol()->lastBornR()>dipol()->lastRealR();
- }
- //dipol()->lastBornR()>dipol()->lastRealR()&&
+
return (zbounds.first<z_&&z_<zbounds.second);
}
-Energy Node::newHardPt(){
- assert(deepHead()->treefactory()->largeNBasis());
-
- if(dipol()->underlyingBornME()->largeNColourCorrelatedME2(
- make_pair(dipol()->bornEmitter(),dipol()->bornSpectator()),
- deepHead()->treefactory()->largeNBasis())==0.)return 1000000000.*GeV;
-
- double x=0.;
- double z_=dipol()->lastZ();
- // II
- if( dipol()->bornEmitter()<2&&dipol()->bornSpectator()<2){
- x =xcomb()->lastX1()*xcomb()->lastX2();
- }
- // IF
- if( dipol()->bornEmitter()<2&&dipol()->bornSpectator()>=2){
- x =dipol()->bornEmitter()==0?xcomb()->lastX1():xcomb()->lastX2();
- }
-
- double chi=1-sqr(2*z_-1-x)/sqr(1-x);
-
- return dipol()->lastPt()/sqrt(abs(chi));
-}
-
-
Ptr<Node>::ptr Node::getLongestHistory_simple(bool normal,double hardScaleFactor) {
Ptr<Node>::ptr res=this;
vector<Ptr<Node>::ptr> temp = getNextOrderedNodes(normal,hardScaleFactor);
Energy minpt=100000.*GeV;
Selector<Ptr<Node>::ptr> subprosel;
while (temp.size()!=0){
minpt=100000.*GeV;
subprosel.clear();
for (vector<Ptr<Node>::ptr>::iterator it=temp.begin();it!=temp.end();it++){
-
-
- if( (*it)->dipol()->underlyingBornME()->largeNColourCorrelatedME2(
- make_pair((*it)->dipol()->bornEmitter(),(*it)->dipol()->bornSpectator()),
- (*it)->deepHead()->treefactory()->largeNBasis())!=0.
+ if( (*it)->dipol()->underlyingBornME()->largeNColourCorrelatedME2(
+ make_pair((*it)->dipol()->bornEmitter(),(*it)->dipol()->bornSpectator()),
+ MH()->largeNBasis())!=0.
){
- Ptr<AlphaSBase>::transient_pointer alphaS = (*it)->xcomb()->eventHandler().SM().alphaSPtr();
-
- if((*it)->nodeME()->dSigHatDR()/nanobarn!=0.){
- subprosel.insert((abs((*it)->dipol()->dSigHatDR() /
- (*it)->nodeME()->dSigHatDR()*alphaS->value(
- (*it)->dipol()->lastPt()*(*it)->dipol()->lastPt()))), (*it));
+
+ if((*it)->nodeME()->dSigHatDR()/nanobarn!=0.){
+ subprosel.insert((abs((*it)->dipol()->dSigHatDR() /
+ (*it)->nodeME()->dSigHatDR()*MH()->as((*it)->dipol()->lastPt()))), (*it));
minpt=min(minpt,(*it)->dipol()->lastPt());
}
-
- // subprosel.insert(1., (*it));
-
-
-
+ //TODO choosehistories
}
}
- // cout<<"\n"<<subprosel.size();
if (subprosel.empty())
return res;
res = subprosel.select(UseRandom::rnd());
-
temp = res->getNextOrderedNodes(true,hardScaleFactor);
}
-
- if (res->parent()&&minpt>50.*GeV) {
- Energy mmmm=minpt;
- Ptr<Node>::ptr xxxx;
- vector<Ptr<Node>::ptr> tempch =children();
-
- for (vector<Ptr<Node>::ptr>::iterator it=tempch.begin();it!=tempch.end();it++){
- if (mmmm>(*it)->dipol()->lastPt()) {
-
-
- mmmm=min(mmmm,(*it)->dipol()->lastPt());
- xxxx=*it;
- }
- }
- if(false&&mmmm*3.<minpt){
- // return res->parent();
- cout<<"\n\ncluster to "<<minpt/GeV<<" dipole: "<<res->dipol()->name()<<" "<<res->dipol()->lastPt()/GeV;
- cout<<"\n"<<res->xcomb()->meMomenta()[2].perp()/GeV<<" mom "<<res->xcomb()->meMomenta()[2]/GeV;
- cout<<"\n"<<res->xcomb()->meMomenta()[3].perp()/GeV<<" mom "<<res->xcomb()->meMomenta()[3]/GeV;
-
- cout<<"\n but "<<mmmm/GeV<<" with "<<xxxx->dipol()->name()<<" "<<xxxx->dipol()->lastPt()/GeV;
- cout<<"\n"<<xxxx->xcomb()->meMomenta()[2].perp()/GeV<<" mom "<<xxxx->xcomb()->meMomenta()[2]/GeV;
- cout<<"\n"<<xxxx->xcomb()->meMomenta()[3].perp()/GeV<<" mom "<<xxxx->xcomb()->meMomenta()[3]/GeV;
-
- cout<<"\norig";
- cout<<"\n"<<xcomb()->meMomenta()[2].perp()/GeV<<" mom "<<xcomb()->meMomenta()[2]/GeV;
- cout<<"\n"<<xcomb()->meMomenta()[3].perp()/GeV<<" mom "<<xcomb()->meMomenta()[3]/GeV;
- cout<<"\n"<<xcomb()->meMomenta()[4].perp()/GeV<<" mom "<<xcomb()->meMomenta()[4]/GeV;
-
- cout<<"\n cl bornR "<<res->dipol()->lastBornR()<<" realR "<<res->dipol()->lastRealR();
- cout<<"\n sm bornR "<<xxxx->dipol()->lastBornR()<<" realR "<<xxxx->dipol()->lastRealR();
-
-
- double deta2 = sqr(xcomb()->meMomenta()[2].eta() - xcomb()->meMomenta()[3].eta());
- double dphi = abs(xcomb()->meMomenta()[2].phi() - xcomb()->meMomenta()[3].phi());
- if ( dphi > Constants::pi ) dphi = 2.0*Constants::pi - dphi;
- double dr = sqrt(deta2 + sqr(dphi));
- cout<<"\nR_23 = "<<dr<<" "<<sqrt(deta2)<<" "<<dphi;
-
- deta2 = sqr(xcomb()->meMomenta()[2].eta() - xcomb()->meMomenta()[4].eta());
- dphi = abs(xcomb()->meMomenta()[2].phi() - xcomb()->meMomenta()[4].phi());
- if ( dphi > Constants::pi ) dphi = 2.0*Constants::pi - dphi;
- dr = sqrt(deta2 + sqr(dphi));
- cout<<"\nR_24 = "<<dr<<" "<<sqrt(deta2)<<" "<<dphi;
-
- deta2 = sqr(xcomb()->meMomenta()[3].eta() - xcomb()->meMomenta()[4].eta());
- dphi = abs(xcomb()->meMomenta()[3].phi() - xcomb()->meMomenta()[4].phi());
- if ( dphi > Constants::pi ) dphi = 2.0*Constants::pi - dphi;
- dr = sqrt(deta2 + sqr(dphi));
- cout<<"\nR_34 = "<<dr<<" "<<sqrt(deta2)<<" "<<dphi;
-
-
-
- }
- }
-
-
- return res;
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-Ptr<Node>::ptr Node::getLongestHistory(bool normal) {
- Ptr<Node>::ptr res;
- vector<Ptr<Node>::ptr> temp = getNextOrderedNodes(normal);
- vector< vector <Ptr<Node>::ptr> > Historys;
- vector< vector <Ptr<Node>::ptr> > newHistorys;
- for ( vector<Ptr<Node>::ptr>::iterator it = temp.begin() ; it != temp.end() ; it++ ) {
- vector<Ptr<Node>::ptr> x;
- x.push_back(this);
- x.push_back(*it);
- Historys.push_back(x);
- newHistorys.push_back(x);
- }
- while(newHistorys.size()!=0){
- vector< vector <Ptr<Node>::ptr> > tempnewhist=newHistorys;
- newHistorys.clear();
- for (vector< vector <Ptr<Node>::ptr> >::iterator newhist =
- tempnewhist.begin() ; newhist != tempnewhist.end() ; newhist++ ) {
- temp =newhist->back()->getNextOrderedNodes(normal);
- // cout<<"\nts "<<temp.size();
- for ( vector<Ptr<Node>::ptr>::iterator it = temp.begin() ; it != temp.end() ; it++ ) {
- vector<Ptr<Node>::ptr> x=*newhist;
- x.push_back(*it);
- Historys.push_back(x);
- newHistorys.push_back(x);
- }
- }
- }
- if(Historys.empty())return this;
-
- size_t longestsize=Historys.back().size();
-
-
- map<Ptr<Node>::ptr, vector<Ptr<Node>::ptr> > Tree;
- // cout<<"\n----";
- // parent considered children
- for (vector< vector <Ptr<Node>::ptr> >::iterator newhist =
- Historys.end() ; newhist != Historys.begin() ; ) {
- newhist--;
- if (newhist->size()<longestsize){
- if(!Tree.empty())break;
- }
-
- longestsize=newhist->size();
-
-
- Ptr<Node>::ptr Born=newhist->back();
- if (Born->dipol()->lastPt()<Born->deepHead()->mergePt())continue;
- Energy Scale=100000000.*GeV;
- if(!Born->children().empty()){
- vector<Ptr<Node>::ptr> temp2 = Born->children();
- for (size_t i=0;i<Born->nodeME()->mePartonData().size();i++){
- if (!Born->nodeME()->mePartonData()[i]->coloured())continue;
- for (size_t j=i+1;j<Born->nodeME()->mePartonData().size();j++){
- if (i==j||!Born->nodeME()->mePartonData()[j]->coloured())continue;
-
-
- assert(false);//This needs to be consistent with the DipoleShower Startscale if the history is not fully ordered!
- Scale=min(Scale,sqrt(2.*Born->nodeME()->lastMEMomenta()[i]*Born->nodeME()->lastMEMomenta()[j]));
- }
- }
-
- }else{
- if(NodeDebug)cout<<"\nlong hist children";
- Born->nodeME()->factory()->scaleChoice()->setXComb(Born->xcomb());
-
- Scale = sqrt(Born->nodeME()->factory()->scaleChoice()->renormalizationScale());
- }
- if (Scale==100000000.*GeV)continue;
-
- while(Born!=this){
-
- // cout<<"\n"<<Born->dipol()->bornEmitter()<<" "<<Born->dipol()->bornSpectator()<<" "<<Born->dipol()->underlyingBornME()->largeNColourCorrelatedME2(
- // make_pair(Born->dipol()->bornEmitter(),Born->dipol()->bornSpectator()),
- // Born->deepHead()->treefactory()->largeNBasis());
- if(Born->dipol()->underlyingBornME()->largeNColourCorrelatedME2(
- make_pair(Born->dipol()->bornEmitter(),Born->dipol()->bornSpectator()),
- Born->deepHead()->treefactory()->largeNBasis())==0.)break;
- if(Born->dipol()->lastPt()>Scale) break;
- if (Born->dipol()->dSigHatDR() ==0.*nanobarn) break;
- if(!Born->inShowerPS(Scale)) break;
- if (Born->dipol()->lastPt()<Born->deepHead()->mergePt())break;
- if(Born->children().empty()&& !Born->xcomb()->willPassCuts() )break;
- Scale=Born->dipol()->lastPt();
- Born=Born->parent();
-
- }
- if (Born!=this)continue;
-
-
- //This history works:
- Born=newhist->back();
-
- while (Born!=this) {
- Tree[Born->parent()].push_back(Born);
- Born = Born->parent();
- }
-
-
- }
-
- res = this;
-
-
- while (Tree.find(res) != Tree.end()) {
- Ptr<Node>::ptr subpro;
- Selector<Ptr<Node>::ptr> subprosel;
- // cout<<"\n\ntree "<< Tree[res].size();
- for ( vector<Ptr<Node>::ptr>::iterator it = Tree[res].begin() ; it != Tree[res].end() ; it++ ) {
- assert((*it)->deepHead()->mergePt()<(*it)->dipol()->lastPt());
- // assert((*it)->xcomb()->willPassCuts() );
- // assert((*it)->dipol()->dSigHatDR() !=0.*nanobarn);
-
-
- //case0 : flat
- if((*it)->deepHead()->chooseHistory()==0){
- subprosel.insert(1., (*it));
- }
-
- //case1 : always take one of the smallest pt
- else if((*it)->deepHead()->chooseHistory()==1){
- if(!subpro) subpro=(*it);
- if((*it)->dipol()->lastPt()<subpro->dipol()->lastPt()){
- if(NodeDebug)cout<<"\n "<<(*it)->dipol()->lastPt()/GeV<<" < "<<subpro->dipol()->lastPt()/GeV;
- subpro=(*it);
-
- subprosel.clear();
- subprosel.insert(1., (*it));
- }else if ((*it)->dipol()->lastPt()==subpro->dipol()->lastPt() ) subprosel.insert(1., (*it));
- }
-
- //case2 : take correspondingly to 1/pt^2
- else if((*it)->deepHead()->chooseHistory()==2){
- subprosel.insert(1./(*it)->dipol()->lastPt()/(*it)->dipol()->lastPt()*GeV2, (*it));
- }
-
- //case3 : correspondingly to the dipoleweight
- else if((*it)->deepHead()->chooseHistory()==3){
- subprosel.insert(abs((*it)->dipol()->dSigHatDR() / nanobarn), (*it));
- }
- //case3 : correspondingly to the dipoleweight
- else if((*it)->deepHead()->chooseHistory()==5){
-
- Ptr<AlphaSBase>::transient_pointer alphaS = (*it)->xcomb()->eventHandler().SM().alphaSPtr();
- subprosel.insert((abs((*it)->dipol()->dSigHatDR() / nanobarn*alphaS->value((*it)->dipol()->lastPt()*(*it)->dipol()->lastPt())/((*it)->parent()->parent()?alphaS->value((*it)->parent()->dipol()->lastPt()*(*it)->parent()->dipol()->lastPt()):0.2))), (*it));
- }
- else if((*it)->deepHead()->chooseHistory()==8){
-
- Ptr<AlphaSBase>::transient_pointer alphaS = (*it)->xcomb()->eventHandler().SM().alphaSPtr();
- // cout<<"\n"<<(*it)->dipol()->dSigHatDR()/nanobarn<<" "<<(*it)->nodeME()->dSigHatDR()/nanobarn<<abs((*it)->dipol()->dSigHatDR() / (*it)->nodeME()->dSigHatDR()*alphaS->value((*it)->dipol()->lastPt()*(*it)->dipol()->lastPt())/((*it)->parent()->parent()?alphaS->value((*it)->parent()->dipol()->lastPt()*(*it)->parent()->dipol()->lastPt()):0.2));
- if((*it)->nodeME()->dSigHatDR()/nanobarn!=0.)
- subprosel.insert((abs((*it)->dipol()->dSigHatDR() / (*it)->nodeME()->dSigHatDR()*alphaS->value((*it)->dipol()->lastPt()*(*it)->dipol()->lastPt())/((*it)->parent()->parent()?alphaS->value((*it)->parent()->dipol()->lastPt()*(*it)->parent()->dipol()->lastPt()):0.2))), (*it));
- }
-
- else if((*it)->deepHead()->chooseHistory()==9){
-
- Ptr<AlphaSBase>::transient_pointer alphaS = (*it)->xcomb()->eventHandler().SM().alphaSPtr();
-
- subprosel.insert((abs((*it)->dipol()->dSigHatDR() / (*it)->nodeME()->dSigHatDR()*alphaS->value((*it)->dipol()->lastPt()*(*it)->dipol()->lastPt())/((*it)->parent()->parent()?alphaS->value((*it)->parent()->dipol()->lastPt()*(*it)->parent()->dipol()->lastPt()):0.2))), (*it));
- // cout<<"\n"<<(abs((*it)->dipol()->dSigHatDR() / nanobarn))<<" -> "
- // <<(abs((*it)->dipol()->dSigHatDR() / nanobarn*alphaS->value((*it)->dipol()->lastPt()*(*it)->dipol()->lastPt())/((*it)->parent()->parent()?alphaS->value((*it)->parent()->dipol()->lastPt()*(*it)->parent()->dipol()->lastPt()):0.2)))<<" ";
-
- }
-
- //case3 : correspondingly to the dipoleweight
- else if((*it)->deepHead()->chooseHistory()==6){
- Ptr<AlphaSBase>::transient_pointer alphaS = (*it)->xcomb()->eventHandler().SM().alphaSPtr();
- subprosel.insert((alphaS->value((*it)->dipol()->lastPt()*(*it)->dipol()->lastPt())/((*it)->parent()->parent()?alphaS->value((*it)->parent()->dipol()->lastPt()*(*it)->parent()->dipol()->lastPt()):0.2)), (*it));
- // cout<<"\n"<<(abs((*it)->dipol()->dSigHatDR() / nanobarn))<<" -> "
- // <<(abs((*it)->dipol()->dSigHatDR() / nanobarn*alphaS->value((*it)->dipol()->lastPt()*(*it)->dipol()->lastPt())/((*it)->parent()->parent()?alphaS->value((*it)->parent()->dipol()->lastPt()*(*it)->parent()->dipol()->lastPt()):0.2)))<<" ";
-
- }
-
- else if((*it)->deepHead()->chooseHistory()==7){
- Ptr<Node>::ptr xx= (*it);
- Energy sum_pt=0*GeV;
- while (xx->parent()){
- sum_pt=xx->dipol()->lastPt();
- xx=xx->parent();
- }
-
- subprosel.insert(1/sum_pt*GeV, (*it));
- // cout<<"\n"<<(abs((*it)->dipol()->dSigHatDR() / nanobarn))<<" -> "
- // <<(abs((*it)->dipol()->dSigHatDR() / nanobarn*alphaS->value((*it)->dipol()->lastPt()*(*it)->dipol()->lastPt())/((*it)->parent()->parent()?alphaS->value((*it)->parent()->dipol()->lastPt()*(*it)->parent()->dipol()->lastPt()):0.2)))<<" ";
-
- }
-
-
- //case4 : correspondingly to the dipoleweight*alpha_s of the last splitting
- else if((*it)->deepHead()->chooseHistory()==4){
- if(!subpro) subpro=(*it);
- if(
- sqrt(2.*(*it)->parent()->nodeME()->lastMEMomenta()[(*it)->dipol()->realEmitter()]*
- (*it)->parent()->nodeME()->lastMEMomenta()[(*it)->dipol()->realEmission()])
- <
- sqrt(2.*subpro->parent()->nodeME()->lastMEMomenta()[subpro->dipol()->realEmitter()]*
- subpro->parent()->nodeME()->lastMEMomenta()[subpro->dipol()->realEmission()])
- ){
- subpro=(*it);
- subprosel.clear();
- subprosel.insert(1., (*it));
- }else if (sqrt(2.*(*it)->parent()->nodeME()->lastMEMomenta()[(*it)->dipol()->realEmitter()]*
- (*it)->parent()->nodeME()->lastMEMomenta()[(*it)->dipol()->realEmission()])
- ==
- sqrt(2.*subpro->parent()->nodeME()->lastMEMomenta()[subpro->dipol()->realEmitter()]*
- subpro->parent()->nodeME()->lastMEMomenta()[subpro->dipol()->realEmission()]) ) subprosel.insert(1., (*it));
- }
- else{
- assert(false);
- }
-
-
- }
-
- if ( !subprosel.empty() ) {
- // cout<<"\nsubprosel.select "<<subprosel.size()<<flush;
- subpro = subprosel.select(UseRandom::rnd());
- // cout<<"\nc"<<flush;
- subprosel.clear();
- }
- // if(alllongestsize!=longestsize)cout<<"\n.."<<subpro->dipol()->name()<<" ";
- if ( subpro ){
- res = subpro;
- // if(res->parent())if(res->parent()->parent())
- // cout<<" "<<(1.*res->clusternumber())/(1.*res->parent()->clusternumber());
- res->addclusternumber();
- }
-
- else return res;
- }
-
- return res;
-}
-
-
-
-double Node::setProjectorStage(bool fast){
-
- nodeME()->projectorStage(0);
-
- if(deepHead()->onlyN()!=-1){
-
- if(!deepHead()->subtractedReal()){
- if(nodeME()->oneLoopNoBorn()&&!NLOunitarized()){
- assert ( onlyN()==int(nodeME()->lastMEMomenta().size()));
- return 1.;
- }else{
- // if(children().empty()&&unitarized())return 1.;
- if(int(nodeME()->lastMEMomenta().size())-onlyN()==0){
- nodeME()->projectorStage(0);
- return 1.;
- }
- if(int(nodeME()->lastMEMomenta().size())-onlyN()==1){
- nodeME()->projectorStage(1);
- return -1.;
- }
- assert(false);
- }
-
-
-
- }else{
- if((N0()+1) == nodeME()->lastMEMomenta().size() || !NLOunitarized()){
- nodeME()->projectorStage(1);
- return 1.;
- }else{
- assert((N0()+1) < nodeME()->lastMEMomenta().size());
- if(int(nodeME()->lastMEMomenta().size())-onlyN()==1){
- nodeME()->projectorStage(1);
- return 1.;
- }
- if(int(nodeME()->lastMEMomenta().size())-onlyN()==2){
- nodeME()->projectorStage(2);
- return -1.;
- }
- }
- assert(false);
- }
- assert(false);
- }
-
-
-
-
-
- double res=1.;
- nodeME()->projectorStage(0);
- if(!deepHead()->subtractedReal()){
- if(nodeME()->oneLoopNoBorn()&&!NLOunitarized())return res;
- if(!children().empty()&&unitarized()){
- res*=2.;
- if (UseRandom::rnd()<0.5){
- nodeME()->projectorStage(1);
- res*=-1.;
- }else{
- nodeME()->projectorStage(0);
- }
- }
- }else{
- if(!NLOunitarized()){
- nodeME()->projectorStage(1);
- return res;
- }
- if((N0()+1) < nodeME()->lastMEMomenta().size()){
- res*=2.;
- if (UseRandom::rnd()<0.5){
- nodeME()->projectorStage(2);
- res*=-1.;
- }else{
- nodeME()->projectorStage(1);
- }
- }else{
- nodeME()->projectorStage(1);
- }
- }
return res;
}
bool Node::headContribution(double hardScaleFactor){
bool allabove=true;
vector<Ptr<Node>::ptr> temp2 = children();
- // set<int> onebelow;
for (vector<Ptr<Node>::ptr>::iterator it = temp2.begin(); it != temp2.end(); it++) {
- // if ((*it)->dipol()->lastPt()<mergePt()){
- // if(onebelow.empty()){
- // onebelow.insert((*it)->dipol()->realEmitter());
- // onebelow.insert((*it)->dipol()->realEmission());
- // }
- // else{
- // if(onebelow.find((*it)->dipol()->realEmitter())==onebelow.end()&&onebelow.find((*it)->dipol()->realEmission())==onebelow.end()){
- // // return false;
- // }
- // }
- // }
- allabove&=(*it)->dipol()->lastPt()>mergePt();
+ allabove&=(*it)->dipol()->lastPt()>MH()->mergePt();
}
if(allabove){
Ptr<Node>::ptr tmpBorn = getLongestHistory_simple(true,hardScaleFactor);
if(!tmpBorn->parent())return false;
}
return true;
}
bool Node::DipolesAboveMergeingScale(Ptr<Node>::ptr& selectedNode,double & sum,Energy& minpt,int& number){
sum=0.;
Selector<Ptr<Node>::ptr> first_subpro;
vector<Ptr<Node>::ptr> tmp=children();
for (vector<Ptr<Node>::ptr>::iterator it = tmp.begin(); it != tmp.end(); it++) {
double Di=-1.* (*it)->dipol()->dSigHatDR(sqr(10.*GeV))/nanobarn;
if ((Di!=0.)&&(*it)->xcomb()->willPassCuts()) {
-
- //vector<Ptr<Node>::ptr> tmp2=(*it)->children();
- //for (vector<Ptr<Node>::ptr>::iterator it2 = tmp2.begin(); it2 != tmp2.end(); it2++) assert(((*it2)->dipol()->lastPt()>deepHead()->mergePt()));
-
assert((*it)->dipol()->clustersafe());
minpt=min(minpt,(*it)->dipol()->lastPt());
assert((*it)->dipol()->clustersafe());
- //assert (((*it)->dipol()->lastPt()>deepHead()->mergePt()));
sum+=Di;
first_subpro.insert(1., (*it));
}
}
number=int(first_subpro.size());
if (number!=0) {
selectedNode=first_subpro.select(UseRandom::rnd());
return true;
}
return false;
}
double Node::calcPsMinusDip(Energy scale){
- return -1.* dipol()->dipMinusPs(sqr(scale),deepHead()->treefactory()->largeNBasis())/nanobarn;
+ return -1.* dipol()->dipMinusPs(sqr(scale),MH()->largeNBasis())/nanobarn;
}
double Node::calcPs(Energy scale){
- return dipol()->ps(sqr(scale),deepHead()->treefactory()->largeNBasis())/nanobarn;
+ return dipol()->ps(sqr(scale),MH()->largeNBasis())/nanobarn;
}
bool Node::diffPsDipBelowMergeingScale(Ptr<Node>::ptr& selectedNode,double & sum,Energy& minpt,int& number){
Selector<Ptr<Node>::ptr> first_subpro;
vector<Ptr<Node>::ptr> tmp=children();
for (vector<Ptr<Node>::ptr>::iterator it = tmp.begin(); it != tmp.end(); it++) {
vector<Ptr<Node>::ptr> tmp2=(*it)->children();
bool isInSomePS=false;
Energy maxx=0.*GeV;
Energy minn=100000*GeV;
for (vector<Ptr<Node>::ptr>::iterator it2 = tmp2.begin(); it2 != tmp2.end(); it2++){
isInSomePS|=(*it)->inShowerPS((*it2)->dipol()->lastPt());
maxx=max(maxx,(*it)->dipol()->lastPt());
minn=min(minn,(*it)->dipol()->lastPt());
}
double Di=0.;
if(isInSomePS||(tmp2.empty())){
-
- //cout<<"\n "<<(*it)->dipol()->lastPt()/GeV;
-
- Di=-1.* (*it)->dipol()->dipMinusPs(sqr(10.*GeV),deepHead()->treefactory()->largeNBasis())/nanobarn;
+
+ Di=-1.* (*it)->dipol()->dipMinusPs(sqr(10.*GeV),MH()->largeNBasis())/nanobarn;
}else{
Di=-1.* (*it)->dipol()->dSigHatDR(sqr(10.*GeV))/nanobarn;
- // -1.*-1.* (*it)->dipol()->dipMinusPs(sqr(10.*GeV),deepHead()->treefactory()->largeNBasis())/nanobarn;
}
- if ((Di!=0.)&&(*it)->xcomb()->willPassCuts()){//&&((*it)->dipol()->lastPt()<deepHead()->mergePt())) {
+ if ((Di!=0.)&&(*it)->xcomb()->willPassCuts()){//&&((*it)->dipol()->lastPt()<MH()->mergePt())) {
vector<Ptr<Node>::ptr> tmp2=(*it)->children();
- for (vector<Ptr<Node>::ptr>::iterator it2 = tmp2.begin(); it2 != tmp2.end(); it2++) assert(((*it2)->dipol()->lastPt()>deepHead()->mergePt()));
+ for (vector<Ptr<Node>::ptr>::iterator it2 = tmp2.begin(); it2 != tmp2.end(); it2++) assert(((*it2)->dipol()->lastPt()>MH()->mergePt()));
if (Di!=0.) {
first_subpro.insert(1., (*it));
minpt=min(minpt,(*it)->dipol()->lastPt());
}
}
}
number=int(first_subpro.size());
if (number!=0) {
selectedNode=first_subpro.select(UseRandom::rnd());
vector<Ptr<Node>::ptr> tmp2=selectedNode->children();
bool isInSomePS=false;
for (vector<Ptr<Node>::ptr>::iterator it2 = tmp2.begin(); it2 != tmp2.end(); it2++){
isInSomePS|=selectedNode->inShowerPS((*it2)->dipol()->lastPt());
}
- if((isInSomePS||(tmp2.empty()))){//selectedNode->dipol()->lastPt()<deepHead()->mergePt()&&
- sum=-1.* selectedNode->dipol()->dipMinusPs(sqr(10.*GeV),deepHead()->treefactory()->largeNBasis())/nanobarn;
+ if((isInSomePS||(tmp2.empty()))){//selectedNode->dipol()->lastPt()<MH()->mergePt()&&
+ sum=-1.* selectedNode->dipol()->dipMinusPs(sqr(10.*GeV),MH()->largeNBasis())/nanobarn;
}else{
sum=-1.* selectedNode->dipol()->dSigHatDR(sqr(10.*GeV))/nanobarn;
}
-
- // sum=-1.*-1.* selectedNode->dipol()->dipMinusPs(sqr(10.*GeV),deepHead()->treefactory()->largeNBasis())/nanobarn;
-
-
-
return true;
}
return false;
}
bool Node::psBelowMergeingScale(Ptr<Node>::ptr& selectedNode,double & sum,Energy& minpt,int& number){
sum=0.;
Selector<Ptr<Node>::ptr> first_subpro;
vector<Ptr<Node>::ptr> tmp=children();
for (vector<Ptr<Node>::ptr>::iterator it = tmp.begin(); it != tmp.end(); it++) {
- double Di=-1.* (*it)->dipol()->ps(sqr(10.*GeV),deepHead()->treefactory()->largeNBasis())/nanobarn;
- if ((Di!=0)&&(*it)->xcomb()->willPassCuts()){//&&((*it)->dipol()->lastPt()<deepHead()->mergePt())) {
+ double Di=-1.* (*it)->dipol()->ps(sqr(10.*GeV),MH()->largeNBasis())/nanobarn;
+ if ((Di!=0)&&(*it)->xcomb()->willPassCuts()){//&&((*it)->dipol()->lastPt()<MH()->mergePt())) {
vector<Ptr<Node>::ptr> tmp2=(*it)->children();
bool isInSomePS=false;
for (vector<Ptr<Node>::ptr>::iterator it2 = tmp2.begin(); it2 != tmp2.end(); it2++){
- assert(((*it2)->dipol()->lastPt()>deepHead()->mergePt())||((*it)->dipol()->lastPt()>deepHead()->mergePt()));
+ assert(((*it2)->dipol()->lastPt()>MH()->mergePt())||((*it)->dipol()->lastPt()>MH()->mergePt()));
isInSomePS|=(*it)->inShowerPS((*it2)->dipol()->lastPt());
}
- //cout<<"\nisInSomePS "<< isInSomePS<< " "<< tmp2.size()<<" "<<maxx/GeV;
+
if(!isInSomePS&&!(tmp2.empty()))continue;
sum+=Di;
if (Di!=0) {
first_subpro.insert(1., (*it));
minpt=min(minpt,(*it)->dipol()->lastPt());
}
}
}
number=int(first_subpro.size());
if (number!=0) {
selectedNode=first_subpro.select(UseRandom::rnd());
return true;
}
return false;
}
bool Node::dipBelowMergeingScale(Ptr<Node>::ptr& selectedNode,double & sum,Energy& minpt,int& number){
sum=0.;
Selector<Ptr<Node>::ptr> first_subpro;
vector<Ptr<Node>::ptr> tmp=children();
for (vector<Ptr<Node>::ptr>::iterator it = tmp.begin(); it != tmp.end(); it++) {
- if (
- true
+ if (true
//(*it)->dipol()->clustersafe() &&
// (*it)->xcomb()->willPassCuts()
) {
bool calcdip=true;
vector<Ptr<Node>::ptr> tmp2=(*it)->children();
for (vector<Ptr<Node>::ptr>::iterator it2 = tmp2.begin(); it2 != tmp2.end(); it2++){
- assert(((*it2)->dipol()->lastPt()>deepHead()->mergePt())||((*it)->dipol()->lastPt()>deepHead()->mergePt()));
- if(((*it2)->dipol()->lastPt()<deepHead()->mergePt())){
- if((*it)->dipol()->lastPt()<deepHead()->mergePt()){
+ assert(((*it2)->dipol()->lastPt()>MH()->mergePt())||((*it)->dipol()->lastPt()>MH()->mergePt()));
+ if(((*it2)->dipol()->lastPt()<MH()->mergePt())){
+ if((*it)->dipol()->lastPt()<MH()->mergePt()){
sum=0;
return false;
}
calcdip=false;
}
}
if(calcdip){
double Di=-1.* (*it)->dipol()->dSigHatDR(sqr(10.*GeV))/nanobarn;
sum+=Di;
minpt=min(minpt,(*it)->dipol()->lastPt());
if (Di!=0) {
first_subpro.insert(1., (*it));
}
}
}
}
number=int(first_subpro.size());
if (number!=0) {
selectedNode=first_subpro.select(UseRandom::rnd());
return true;
}
return false;
}
-
-
-
-
-
-
-
-
-
-
-
-
-Ptr<MFactory>::ptr Node::treefactory(){return theTreeFactory;}
-
-void Node::treefactory(Ptr<MFactory>::ptr x){theTreeFactory=x;}
-
IBPtr Node::clone() const {
return new_ptr(*this);
}
IBPtr Node::fullclone() const {
return new_ptr(*this);
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void Node::persistentOutput(PersistentOStream & os) const {
- os << theheadxcomb <<
- thexcomb <<
- thenodeMEPtr <<
- thedipol <<
- thechildren <<
- theparent <<
- theProjectors <<
- theDeepHead <<
- theCutStage <<
- ounit(theMergePt, GeV) <<
- ounit(theCentralMergePt, GeV) <<
- smearing <<
- isthissafe <<
- ounit(theIRsafePt, GeV) <<
- theIRCsafeRatio <<
- theDeepProStage <<
- clustersafer <<
- isOrdered <<
- theOrderdSteps <<
- theNeedFullOrderedHistory <<
- theN0 <<
- theOnlyN <<
- theN <<
- theM <<
- theSudakovSteps <<
- isUnitarized <<
- isNLOUnitarized <<
- ounit(theVetoPt, GeV) <<
- ounit(theRunningPt, GeV) <<
- needsVetoedShower <<
- theCalculateInNode <<
- theNumberOfSplittings <<
- theSubtractedReal <<
- theVirtualContribution <<
- thefiniteDipoles <<
- thesubCorrection <<
- theclusteredto <<
- theTreeFactory <<
- theChooseHistory;
+ os <<
+ theheadxcomb << thexcomb << thenodeMEPtr <<
+ thedipol << thechildren << theparent <<
+ theProjectors << theDeepHead << theCutStage <<
+ theDeepProStage << clustersafer << ounit(theVetoPt, GeV) <<
+ ounit(theRunningPt, GeV) << theSubtractedReal << theVirtualContribution
+ <<theMergingHelper;
// *** ATTENTION *** os << ; // Add all member variable which should be written persistently here.
}
void Node::persistentInput(PersistentIStream & is, int) {
- is >> theheadxcomb>>
- thexcomb>>
- thenodeMEPtr>>
- thedipol>>
- thechildren>>
- theparent>>
- theProjectors>>
- theDeepHead>>
- theCutStage>>
- iunit(theMergePt, GeV)>>
- iunit(theCentralMergePt, GeV)>>
- smearing>>
- isthissafe>>
- iunit(theIRsafePt, GeV)>>
- theIRCsafeRatio>>
- theDeepProStage>>
- clustersafer>>
- isOrdered>>
- theOrderdSteps>>
- theNeedFullOrderedHistory>>
- theN0>>
- theOnlyN>>
- theN>>
- theM>>
- theSudakovSteps>>
- isUnitarized>>
- isNLOUnitarized>>
- iunit(theVetoPt, GeV)>>
- iunit(theRunningPt, GeV)>>
- needsVetoedShower>>
- theCalculateInNode>>
- theNumberOfSplittings>>
- theSubtractedReal>>
- theVirtualContribution>>
- thefiniteDipoles>>
- thesubCorrection>>
- theclusteredto>> theTreeFactory>>
- theChooseHistory;
+ is >>
+ theheadxcomb >> thexcomb >> thenodeMEPtr >>
+ thedipol >> thechildren >> theparent >>
+ theProjectors >> theDeepHead >> theCutStage >>
+ theDeepProStage >> clustersafer >> iunit(theVetoPt, GeV) >>
+ iunit(theRunningPt, GeV) >> theSubtractedReal >> theVirtualContribution
+ >>theMergingHelper;
// *** ATTENTION *** is >> ; // Add all member variable which should be read persistently here.
}
// *** 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<Node,Interfaced> describeHerwigNode("Herwig::Node", "HwDipoleShower.so");
void Node::Init() {
static ClassDocumentation<Node> documentation("There is no documentation for the Node class");
}
diff --git a/DipoleShower/Merging2/Node.h b/DipoleShower/Merging2/Node.h
--- a/DipoleShower/Merging2/Node.h
+++ b/DipoleShower/Merging2/Node.h
@@ -1,652 +1,446 @@
// -*- C++ -*-
#ifndef Herwig_Node_H
#define Herwig_Node_H
//
// This is the declaration of the Node class.
//
//#include "Node.fh"
#include "MFactory.fh"
+#include "Merger.h"
#include "ThePEG/Interface/Interfaced.h"
#include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.fh"
#include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.fh"
#include "Herwig/DipoleShower/Base/DipoleEventRecord.h"
#include "ThePEG/MatrixElement/MEBase.h"
#include "ThePEG/Config/Pointers.h"
#include <vector>
namespace Herwig {
using namespace ThePEG;
/**
* Here is the documentation of the Node class.
*
* @see \ref NodeInterfaces "The interfaces"
* defined for Node.
*/
/**
* Define the SafeClusterMap type map<pair<pair<emitter,emmision>,spectator >
* ,pair<first-clustering,second-clustering> >
*/
typedef map<pair<pair<int,int>,int >,pair<bool,bool> > SafeClusterMap;
template <typename T>
struct Nodecounter
{
Nodecounter()
{
objects_created++;
objects_alive++;
}
virtual ~Nodecounter()
{
--objects_alive;
}
static int objects_created;
static int objects_alive;
};
template <typename T> int Nodecounter<T>::objects_created( 0 );
template <typename T> int Nodecounter<T>::objects_alive( 0 );
class Node : public Interfaced,Nodecounter<Node> {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor. Do not use!
*/
Node();
- Node(Ptr<MatchboxMEBase>::ptr nodeME, int deepprostage,int cutstage, bool nFOH);
+ Node(Ptr<MatchboxMEBase>::ptr nodeME, int deepprostage,int cutstage);
- Node(Ptr<Node>::ptr deephead, Ptr<Node>::ptr head, Ptr<SubtractionDipole>::ptr dipol, Ptr<MatchboxMEBase>::ptr nodeME,
- int deepprostage, int cutstage);
+ Node(Ptr<Node>::ptr deephead,
+ Ptr<Node>::ptr head,
+ Ptr<SubtractionDipole>::ptr dipol,
+ Ptr<MatchboxMEBase>::ptr nodeME,
+ int deepprostage, int cutstage);
/**
* The destructor.
*/
virtual ~Node();
//@}
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();
public:
/** get children from vector<Ptr<MatchboxMEBase>::ptr> */
void birth(vector<Ptr<MatchboxMEBase>::ptr> vec);
/** recursive setXComb. proStage is the number of clusterings
* before the projectors get filled. */
void setXComb(tStdXCombPtr xc, int proStage);
-
-
- double setProjectorStage(bool fast=false);
-
bool headContribution(double hardscalefactor);
bool DipolesAboveMergeingScale(Ptr<Node>::ptr& selectedNode,double& sum,Energy& minpt,int& number);
bool diffPsDipBelowMergeingScale(Ptr<Node>::ptr& selectedNode,double & sum,Energy& minpt,int& number);
bool psBelowMergeingScale(Ptr<Node>::ptr& selectedNode,double & sum,Energy& minpt,int& number);
bool dipBelowMergeingScale(Ptr<Node>::ptr& selectedNode,double & sum,Energy& minpt,int& number);
double calcPsMinusDip(Energy scale);
double calcPs(Energy scale);
/** recursive flush caches and clean up XCombs. */
void flushCaches();
/** recursive clearKinematics */
void clearKinematics();
/** recursive setKinematics */
void setKinematics();
/** recursive generateKinematics using tilde kinematics of the dipoles*/
bool generateKinematics(const double *r, int stage,Energy2 shat);
void firstgenerateKinematics(const double *r, int stage,Energy2 shat);
/** returns the matrix element pointer */
const Ptr<MatchboxMEBase>::ptr nodeME()const;
/** access the matrix element pointer */
Ptr<MatchboxMEBase>::ptr nodeME();
/** the parent node */
Ptr<Node>::ptr parent()const {
return theparent;
}
/** map of children nodes*/
vector< Ptr<Node>::ptr > children() {
return thechildren;
}
Ptr<Node>::ptr randomChild();
bool allAbove(Energy pt);
bool isInHistoryOf(Ptr<Node>::ptr other);
/** set the first node (godfather). only use in factory*/
void deepHead(Ptr<Node>::ptr deephead) {
theDeepHead = deephead;
}
/** return the first node*/
Ptr<Node>::ptr deepHead() const {
return theDeepHead;
}
/** insert nodes to projector vector */
void Projector(double a, Ptr<Node>::ptr pro) {
- pair<double,Ptr<Node>::ptr> p;
- p.first = a;
- p.second = pro;
- theProjectors.push_back(p);
+ theProjectors.push_back(make_pair(a,pro));
}
/** insert nodes to projector vector */
vector< pair <double , Ptr<Node>::ptr > > Projector() {
return theProjectors;
}
/** returns the dipol of the node. */
Ptr<SubtractionDipole>::ptr dipol() const;
/** set the xcomb of the node */
void xcomb(StdXCombPtr xc) {
thexcomb = xc;
}
/** return the xcomb */
StdXCombPtr xcomb() const {
return thexcomb;
}
/** set the headxcomb of the node */
void headxcomb(StdXCombPtr xc) {
thexcomb = xc;
}
/** return the xcomb */
StdXCombPtr headxcomb() const {
return theheadxcomb;
}
- Energy mergePt() {
- return theMergePt;
- }
-
- void mergePt(Energy x) {
- theMergePt = x;
- }
-
- Energy centralMergePt() {
- return theCentralMergePt;
- }
-
- void centralMergePt(Energy x) {
- theCentralMergePt = x;
- }
-
- double smear();
+
Energy vetoPt() const {
return theVetoPt;
}
void vetoPt(Energy x) {
theVetoPt=x;
}
Energy runningPt(){return theRunningPt;}
void runningPt(Energy x){theRunningPt=x;}
int cutStage() const {
return theCutStage;
}
-
- void N(unsigned int N) {
- theN = N;
- }
-
- void N0(unsigned int N) {
- theN0 = N;
- }
-
-
- void onlyN( int N) {
- theOnlyN = N;
- }
- int onlyN( ) {
- return theOnlyN ;
- }
-
-
-
- unsigned int N() const {
- return theN;
- }
-
- unsigned int N0() const {
- return theN0;
- }
-
- void M(unsigned int M) {
- theM = M;
- }
- unsigned int M() const {
- return theM;
- }
-
- unsigned int sudakovSteps() const {
- return theSudakovSteps;
- }
-
+
unsigned int DeepProStage() const {
return theDeepProStage;
}
- void irsavePt(Energy x) {
- theIRsafePt = x;
- }
-
- Energy irsavePt() {
- return theIRsafePt;
- }
-
- void irsaveRatio(double x) {
- theIRCsafeRatio = x;
- }
-
- double irsaveRatio() {
- return theIRCsafeRatio;
- }
+
SafeClusterMap clusterSafe() const {
return clustersafer;
}
bool ordered() {
return isOrdered;
}
int orderedSteps() const {
return theOrderdSteps;
}
void setOrderedSteps(int x) {
theOrderdSteps = x;
}
bool isThisSafe() const {
return isthissafe;
}
vector<Ptr<Node>::ptr> getNextOrderedNodes(bool normal=true,double hardscalefactor=1.);
bool inShowerPS(Energy hardpt);
- Energy newHardPt();
vector<Ptr<Node>::ptr> getNextFullOrderedNodes();
bool hasFullHistory();
- bool unitarized() const {
- return isUnitarized;
- }
- void unitarized(bool is) {
- isUnitarized = is;
- }
-
- bool NLOunitarized() const {
- return isNLOUnitarized;
- }
- void NLOunitarized(bool is) {
- isNLOUnitarized = is;
- }
-
- bool needFullOrderedHistory() {
- return theNeedFullOrderedHistory;
- }
-
- Ptr<Node>::ptr getLongestHistory(bool normal=true);
Ptr<Node>::ptr getLongestHistory_simple(bool normal=true,double hardscalefactor=1.);
-
- tSubProPtr showeredSub(){return theShoweredSub;}
-
- void showeredSub( tSubProPtr x){theShoweredSub=x;}
-
- DipoleEventRecord& eventRec(){return theEventRec;}
- void eventRec(DipoleEventRecord& x){theEventRec=x;}
-
- bool VetoedShower(){return needsVetoedShower;}
- void VetoedShower(bool x){needsVetoedShower = x;}
- bool calculateInNode(){return theCalculateInNode;}
- void calculateInNode(bool x) {
- theCalculateInNode = x;
- }
bool subtractedReal() {
return theSubtractedReal;
}
void subtractedReal(bool x) {
theSubtractedReal = x;
}
bool virtualContribution() {
return theVirtualContribution ;
}
void virtualContribution(bool x) {
theVirtualContribution = x;
}
+ Ptr<Merger>::ptr MH(){return theMergingHelper;};
- bool finiteDipoles() {
- return thefiniteDipoles;
- }
- void finiteDipoles(bool x) {
- thefiniteDipoles = x;
- }
-
-
-
- bool subCorrection() {
- return thesubCorrection;
- }
- void subCorrection(bool x) {
- thesubCorrection = x;
- }
- void chooseHistory(int x){
- theChooseHistory=x;
- }
-
- int chooseHistory(){
- return theChooseHistory;
- }
-
- int numberOfSplittings(){
- return theNumberOfSplittings;
- }
- void numberOfSplittings(int n){
- theNumberOfSplittings=n;
- }
-
-
- Ptr<MFactory>::ptr treefactory();
-
- void treefactory(Ptr<MFactory>::ptr x);
-
-
- void addclusternumber(){theclusteredto++;}
- int clusternumber(){return theclusteredto;}
-
-
private:
StdXCombPtr theheadxcomb;
StdXCombPtr thexcomb;
/** the Matrixelement representing this node. */
Ptr<MatchboxMEBase>::ptr thenodeMEPtr;
/** the dipol used to substract
* and generate kinematics using tilde kinematics */
Ptr<SubtractionDipole>::ptr thedipol;
/** vector of the children node*/
vector< Ptr<Node>::ptr > thechildren;
/** the parent node*/
Ptr<Node>::ptr theparent;
/** The nodes of the projection stage. */
vector< pair <double,Ptr<Node>::ptr> > theProjectors;
/** The godfather node of whole tree.(Firstnode) */
Ptr<Node>::ptr theDeepHead;
/**
* The CutStage is number of clusterings which are possible without
* introducing a merging scale to cut away singularities.
* -> subtracted MEs have the CutStage 1.
* -> virtual and normal tree level ME get 0.
*/
int theCutStage;
/**
- * If any clustering below the CutStage has a lower pT than the MergePt
- * the phase space point has to be thrown away.
- */
-
- static Energy theMergePt;
-
-
- static Energy theCentralMergePt;
-
- static double smearing;
-
- /**
* all nodes downstream have pt over merged pt.
*/
bool isthissafe;
-
- /**
- * If any clustering below the CutStage has a lower pT than the MergePt
- * the phase space point has to be thrown away.
- */
-
- static Energy theIRsafePt;
-
- /**
- * If any clustering below the CutStage has a lower (pT/shat)^2 than this value
- * the phase space point has to be thrown away.
- */
-
- static double theIRCsafeRatio;
+
/**
* The DeepProStage is the number of additional partons of the DeepHead
* compared to the lowest multiplicity.
*/
int theDeepProStage;
/**
* For [[Emitter,Emission],Spectator] the mapped pair gives
* information if the first and the second cluster is safe.
*/
SafeClusterMap clustersafer;
bool isOrdered;
int theOrderdSteps;
bool theNeedFullOrderedHistory;
- static unsigned int theN0;
-
- static int theOnlyN;
-
- static unsigned int theN;
-
- static unsigned int theM;
-
unsigned int theSudakovSteps;
- static bool isUnitarized;
-
- static bool isNLOUnitarized;
-
Energy theVetoPt;
Energy theRunningPt;
- tSubProPtr theShoweredSub;
-
- DipoleEventRecord theEventRec;
-
- bool needsVetoedShower;
-
- bool theCalculateInNode;
-
- int theNumberOfSplittings;
-
bool theSubtractedReal;
bool theVirtualContribution;
- bool thefiniteDipoles;
- bool thesubCorrection;
-
- int theclusteredto;
-
-
-
-
- Ptr<MFactory>::ptr theTreeFactory;
-
-
- // 1 == always take one of the smallest Pts.
- // 2 == always take one of the highest Pts.
- // 3 == choose correspondingly to the dipolweight.
-
- static int theChooseHistory;
-
+ Ptr<Merger>::ptr theMergingHelper;
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;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
Node & operator=(const Node &);
};
}
#endif /* Herwig_Node_H */

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 6:18 PM (1 d, 21 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3800601
Default Alt Text
(202 KB)

Event Timeline