8 #ifndef Pythia8_DireHistory_H
9 #define Pythia8_DireHistory_H
11 #include "Pythia8/Basics.h"
12 #include "Pythia8/BeamParticle.h"
13 #include "Pythia8/Event.h"
14 #include "Pythia8/Info.h"
15 #include "Pythia8/ParticleData.h"
16 #include "Pythia8/PartonLevel.h"
17 #include "Pythia8/PythiaStdlib.h"
18 #include "Pythia8/Settings.h"
19 #include "Pythia8/StandardModel.h"
20 #include "Pythia8/MergingHooks.h"
21 #include "Pythia8/SimpleWeakShowerMEs.h"
22 #include "Pythia8/DireWeightContainer.h"
39 int radPos()
const {
return emittor; }
40 int emtPos()
const {
return emitted; }
41 int recPos()
const {
return recoiler; }
42 const Particle* rad() {
return radSave; }
43 const Particle* emt() {
return emtSave; }
44 const Particle* rec() {
return recSave; }
72 string name()
const {
return splitName;}
96 emitted = inSystem.emitted;
97 emittor = inSystem.emittor;
98 recoiler = inSystem.recoiler;
99 partner = inSystem.partner;
100 pTscale = inSystem.pTscale;
101 flavRadBef = inSystem.flavRadBef;
102 spinRadBef = inSystem.spinRadBef;
103 radBef = inSystem.radBef;
104 recBef = inSystem.recBef;
105 radSave = inSystem.radSave;
106 emtSave = inSystem.emtSave;
107 recSave = inSystem.recSave;
108 splitName = inSystem.splitName;
113 { emitted = c.emitted; emittor = c.emittor; recoiler = c.recoiler;
114 partner = c.partner; pTscale = c.pTscale; flavRadBef = c.flavRadBef;
115 spinRadBef = c.spinRadBef; radBef = c.radBef; recBef = c.recBef;
116 radSave = c.radSave; emtSave = c.emtSave; recSave = c.recSave;
117 splitName = c.splitName;}
return *
this; }
120 DireClustering(
int emtPosIn,
int radPosIn,
int recPosIn,
int partnerIn,
122 const Particle* recIn,
string splitNameIn,
123 int flavRadBefIn = 0,
int spinRadBefIn = 9,
124 int radBefIn = 0,
int recBefIn = 0)
125 : emitted(emtPosIn), emittor(radPosIn), recoiler(recPosIn),
126 partner(partnerIn), pTscale(pTscaleIn), radSave(radIn), emtSave(emtIn),
127 recSave(recIn), flavRadBef(flavRadBefIn), spinRadBef(spinRadBefIn),
128 radBef(radBefIn), recBef(recBefIn), splitName(splitNameIn) {}
131 double pT()
const {
return pTscale; }
134 double mass()
const {
135 double sik = 2.*radSave->p()*recSave->p();
136 double sij = 2.*radSave->p()*emtSave->p();
137 double sjk = 2.*emtSave->p()*recSave->p();
140 if ( radSave->isFinal() && recSave->isFinal()) m2 = sik+sij+sjk;
141 else if ( radSave->isFinal() && !recSave->isFinal()) m2 =-sik+sij-sjk;
142 else if (!radSave->isFinal() && recSave->isFinal()) m2 =-sik-sij+sjk;
143 else if (!radSave->isFinal() && !recSave->isFinal()) m2 = sik-sij-sjk;
186 MergingHooksPtr mergingHooksPtrIn,
192 shared_ptr<DireTimes> fsrIn,
193 shared_ptr<DireSpace> isrIn,
198 double clusterProbIn,
199 double clusterCouplIn,
200 double prodOfProbsIn,
201 double prodOfProbsFullIn,
206 {
for (
int i = 0, N = children.size(); i < N; ++i )
delete children[i]; }
209 bool projectOntoDesiredHistories();
228 double weightMEC() {
return MECnum/MECden; }
259 double weight_UNLOPS_CORRECTION(
int order,
PartonLevel* trial,
264 bool foundAllowedHistories() {
265 return (children.size() > 0 && foundAllowedPath); }
267 bool foundOrderedHistories() {
268 return (children.size() > 0 && foundOrderedPath); }
270 bool foundCompleteHistories() {
271 return (children.size() > 0 && foundCompletePath); }
274 void getStartingConditions(
const double RN,
Event& outState );
276 bool getClusteredEvent(
const double RN,
int nSteps,
Event& outState);
278 bool getFirstClusteredEventAboveTMS(
const double RN,
int nDesired,
279 Event& process,
int & nPerformed,
bool updateProcess =
true );
285 Event lowestMultProc(
const double RN) {
287 return (select(RN)->state);
291 double getPDFratio(
int side,
bool forSudakov,
bool useHardPDF,
292 int flavNum,
double xNum,
double muNum,
293 int flavDen,
double xDen,
double muDen);
297 void printHistory(
const double RN );
313 static const int NTRIAL;
317 static const double MCWINDOW, MBWINDOW, PROBMAXFAC;
321 void setScalesInHistory();
330 void findPath(vector<int>& out);
350 void setScales( vector<int> index,
bool forward);
359 void scaleCopies(
int iPart,
const Event& refEvent,
double rho);
365 void setEventScales();
370 int type = (!mother) ? 0
371 : ( clusterIn.rad()->isFinal() && clusterIn.rec()->isFinal()) ? 2
372 : ( clusterIn.rad()->isFinal() && !clusterIn.rec()->isFinal()) ? 1
373 : (!clusterIn.rad()->isFinal() && clusterIn.rec()->isFinal()) ?
375 cout << scientific << setprecision(6);
376 cout <<
" size " << state.size() <<
" scale " << scale
377 <<
" clusterIn " << clusterIn.pT() <<
" state.scale() " << state.scale()
378 <<
" scaleEffective " << scaleEffective;
379 if (type==-2) cout <<
"\t\t" <<
" II splitting emt="
380 << clusterIn.emt()->id() << endl;
381 if (type==-1) cout <<
"\t\t" <<
" IF splitting emt="
382 << clusterIn.emt()->id() << endl;
383 if (type== 1) cout <<
"\t\t" <<
" FI splitting emt="
384 << clusterIn.emt()->id() << endl;
385 if (type== 2) cout <<
"\t\t" <<
" FF splitting emt="
386 << clusterIn.emt()->id() << endl;
387 if ( mother ) mother->printScales();
391 bool trimHistories();
393 void remove(){ doInclude =
false; }
395 bool keep(){
return doInclude; }
400 bool isOrderedPath(
double maxscale);
402 bool hasScalesAboveCutoff();
404 bool followsInputPath();
408 bool allIntermediateAboveRhoMS(
double rhoms,
bool good =
true );
410 bool foundAnyOrderedPaths();
416 Event clusteredState(
int nSteps);
436 double weight(
PartonLevel* trial,
double as0,
double aem0,
438 AlphaEM * aemFSR,
AlphaEM * aemISR,
double& asWeight,
double& aemWeight,
442 double weightALPHAS(
double as0,
AlphaStrong * asFSR,
443 AlphaStrong * asISR,
int njetMin = -1 ,
int njetMax = -1 );
445 double weightALPHAEM(
double aem0,
AlphaEM * aemFSR,
446 AlphaEM * aemISR,
int njetMin = -1,
int njetMax = -1 );
448 vector<double> weightCouplings();
449 vector<double> weightCouplingsDenominator();
451 double weightPDFs(
double maxscale,
double pdfScale,
int njetMin = -1,
454 double weightEmissions(
PartonLevel* trial,
int type,
int njetMin,
455 int njetMax,
double maxscale );
456 vector<double> weightEmissionsVec(
PartonLevel* trial,
int type,
457 int njetMin,
int njetMax,
double maxscale );
460 double weightFirst(
PartonLevel* trial,
double as0,
double muR,
465 double weightFirstALPHAS(
double as0,
double muR,
AlphaStrong * asFSR,
469 double weightFirstALPHAEM(
double aem0,
double muR,
AlphaEM * aemFSR,
473 double weightFirstPDFs(
double as0,
double maxscale,
double pdfScale,
477 double weightFirstEmissions(
PartonLevel* trial,
double as0,
double maxscale,
481 double hardFacScale(
const Event& event);
483 double hardRenScale(
const Event& event);
485 double hardProcessScale(
const Event& event);
487 double hardStartScale(
const Event& event);
497 vector<double> doTrialShower(
PartonLevel* trial,
int type,
double maxscale,
498 double minscale = 0.);
501 bool updateind(vector<int> & ind,
int i,
int N);
504 vector<double> countEmissions(
PartonLevel* trial,
double maxscale,
505 double minscale,
int showerType,
double as0,
AlphaStrong * asFSR,
506 AlphaStrong * asISR,
int N,
bool fixpdf,
bool fixas);
510 double monteCarloPDFratios(
int flav,
double x,
double maxScale,
511 double minScale,
double pdfScale,
double asME,
Rndm* rndmPtr);
516 bool onlyOrderedPaths();
521 bool onlyAllowedPaths();
530 bool isAllowed,
bool isComplete);
535 vector<DireClustering> getAllClusterings(
const Event& event);
536 vector<DireClustering> getClusterings(
int emt,
int rad,
const Event& event);
538 void attachClusterings (vector<DireClustering>& clus,
int iEmt,
int iRad,
539 int iRec,
int iPartner,
double pT,
string name,
const Event& event);
560 double pdfForSudakov();
564 double hardProcessME(
const Event& event);
565 double hardProcessCouplings(
const Event& event,
int order = 0,
566 double renormMultFac = 1.,
AlphaStrong* alphaS = NULL,
567 AlphaEM* alphaEM = NULL,
bool fillCouplCounters =
false,
568 bool with2pi =
true);
582 int getRadBeforeFlav(
const int radAfter,
const int emtAfter,
589 int getRadBeforeSpin(
const int radAfter,
const int emtAfter,
590 const int spinRadAfter,
const int spinEmtAfter,
599 int getRadBeforeCol(
const int radAfter,
const int emtAfter,
608 int getRadBeforeAcol(
const int radAfter,
const int emtAfter,
616 int getColPartner(
const int in,
const Event& event);
623 int getAcolPartner(
const int in,
const Event& event);
633 vector<int> getReclusteredPartners(
const int rad,
const int emt,
652 bool getColSinglet(
const int flavType,
const int iParton,
653 const Event& event, vector<int>& exclude,
654 vector<int>& colSinglet);
660 bool isColSinglet(
const Event& event, vector<int> system);
667 bool isFlavSinglet(
const Event& event,
668 vector<int> system,
int flav=0);
679 bool connectRadiator(
Particle& radiator,
const int radType,
680 const Particle& recoiler,
const int recType,
681 const Event& event );
697 int FindCol(
int col,
int iExclude1,
int iExclude2,
698 const Event& event,
int type,
bool isHardIn);
706 int FindParticle(
const Particle& particle,
const Event& event,
707 bool checkStatus =
true );
714 bool allowedClustering(
int rad,
int emt,
int rec,
int partner,
715 string name,
const Event& event );
717 bool hasConnections(
int arrSize,
int nIncIDs[],
int nOutIDs[]);
718 bool canConnectFlavs( map<int,int> nIncIDs, map<int,int> nOutIDs);
725 bool isSinglett(
int rad,
int emt,
int rec,
const Event& event );
733 bool validEvent(
const Event& event );
742 double choseHardScale(
const Event& event )
const;
746 int getCurrentFlav(
const int)
const;
750 double getCurrentX(
const int)
const;
752 double getCurrentZ(
const int rad,
const int rec,
const int emt,
753 int idRadBef = 0)
const;
756 double pTLund(
const Event& event,
int radAfterBranch,
int emtAfterBranch,
757 int recAfterBranch,
string name);
761 int posChangedIncoming(
const Event& event,
bool before);
763 vector<int> getSplittingPos(
const Event& event,
int type);
769 double pdfFactor(
const Event& process,
const Event& event,
const int type,
770 double pdfScale,
double mu );
775 double integrand(
int flav,
double x,
double scaleInt,
double z);
779 vector<int> posFlavCKM(
int flav);
784 bool checkFlavour(vector<int>& flavCounts,
int flavRad,
int flavRadBef,
789 bool isQCD2to2(
const Event& event);
793 bool isEW2to1(
const Event& event);
796 bool isMassless2to2(
const Event& event);
798 bool isDIS2to2(
const Event& event);
801 bool mayHaveEffectiveVertex(
string process, vector<int> in, vector<int> out);
804 void setSelectedChild();
806 void setGoodChildren();
807 void setGoodSisters();
808 void setProbabilities();
814 bool hasTag(
string key) {
815 if(find(tagSave.begin(), tagSave.end(), key) != tagSave.end())
820 void setEffectiveScales();
823 double getShowerPluginScale(
const Event& event,
int rad,
int emt,
int rec,
824 string name,
string key,
double scalePythia);
826 pair<int,double> getCoupling(
const Event& event,
int rad,
int emt,
827 int rec,
string name);
830 map<string,int> count = map<string,int>());
850 vector<DireHistory *> children;
851 vector<DireHistory *> goodSisters;
859 map<double,DireHistory *> paths;
870 map<double,DireHistory *> goodBranches, badBranches;
876 double sumGoodBranches, sumBadBranches;
880 bool foundOrderedPath;
884 bool foundAllowedPath;
888 bool foundCompletePath;
890 bool foundOrderedChildren;
895 double scale, scaleEffective, couplEffective;
897 bool allowedOrderingPath;
905 double clusterProb, clusterCoupl, prodOfProbs, prodOfProbsFull;
909 int iReclusteredOld, iReclusteredNew;
915 double MECnum, MECden, MECcontrib;
917 vector<int> goodChildren;
920 MergingHooksPtr mergingHooksPtr;
949 shared_ptr<DireTimes> fsr;
950 shared_ptr<DireSpace> isr;
958 bool doSingleLegSudakovs;
960 vector<string> tagSave;
964 if (mother)
return mother->probMax();
967 void updateProbMax(
double probIn,
bool isComplete =
false) {
968 if (mother) mother->updateProbMax(probIn, isComplete);
969 if (!isComplete && !foundCompletePath)
return;
970 if (abs(probIn) > probMaxSave) probMaxSave = probIn;
973 int depth, minDepthSave;
975 if ( mother )
return mother->minDepth();
978 void updateMinDepth(
int depthIn) {
979 if ( mother )
return mother->updateMinDepth(depthIn);
980 minDepthSave = (minDepthSave>0) ? min(minDepthSave,depthIn) : depthIn;
983 map<string,int> couplingPowCount;
991 #endif // end Pythia8_DireHistory_H