herwig is hosted by Hepforge, IPPP Durham
Herwig  7.2.1
MPIHandler.h
1 // -*- C++ -*-
2 //
3 // MPIHandler.h is a part of Herwig - A multi-purpose Monte Carlo event generator
4 // Copyright (C) 2002-2019 The Herwig Collaboration
5 //
6 // Herwig is licenced under version 3 of the GPL, see COPYING for details.
7 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
8 //
9 #ifndef HERWIG_MPIHandler_H
10 #define HERWIG_MPIHandler_H
11 //
12 // This is the declaration of the MPIHandler class.
13 //
14 #include "ThePEG/Interface/Interfaced.h"
15 #include "ThePEG/Handlers/StandardEventHandler.h"
16 #include "ThePEG/Repository/EventGenerator.h"
17 #include "Herwig/PDT/StandardMatchers.h"
18 #include "Herwig/Utilities/GSLBisection.h"
19 //#include "Herwig/Utilities/GSLMultiRoot.h"
20 #include "Herwig/Utilities/GSLIntegrator.h"
21 #include "Herwig/Shower/UEBase.h"
22 
23 #include <cassert>
24 #include "ProcessHandler.h"
25 #include "MPIHandler.fh"
26 
27 
28 namespace Herwig {
29 using namespace ThePEG;
30 
45 class MPIHandler: public UEBase {
46 
50  static const unsigned int maxScatters_ = 99;
51 
55  friend struct Eikonalization;
56  friend struct TotalXSecBisection;
57  friend struct slopeAndTotalXSec;
58  friend struct slopeInt;
59  friend struct slopeBisection;
60 
61 public:
62 
64  typedef vector<SubHdlPtr> SubHandlerList;
65 
67  typedef vector<CutsPtr> CutsList;
68 
70  typedef vector<ProHdlPtr> ProcessHandlerList;
71 
73  typedef vector<CrossSection> XSVector;
74 
76  typedef pair<unsigned int, unsigned int> MPair;
77 
83  MPIHandler(): softMult_(0), identicalToUE_(-1),
84  PtOfQCDProc_(-1.0*GeV), Ptmin_(-1.0*GeV),
85  hardXSec_(0*millibarn), softXSec_(0*millibarn),
86  totalXSecExp_(0*millibarn),
87  softMu2_(ZERO), beta_(100.0/GeV2),
88  algorithm_(2), numSubProcs_(0),
89  colourDisrupt_(0.0), softInt_(true), twoComp_(true),
90  DLmode_(2), avgNhard_(0.0), avgNsoft_(0.0),
91  energyExtrapolation_(3), EEparamA_(0.6*GeV),
92  EEparamB_(37.5*GeV), refScale_(7000.*GeV),
93  pT0_(2.875*GeV), b_(0.3101), offset_(622.204*GeV) {}
94 
98  virtual ~MPIHandler(){}
100 
101 public:
102 
105 
106  /*
107  * @return true if for this beam setup MPI can be generated
108  */
109  virtual bool beamOK() const;
110 
114  virtual bool softInt() const {return softInt_;}
115 
121  virtual unsigned int softMultiplicity() const {return softMult_;}
122 
127  virtual unsigned int multiplicity(unsigned int sel=0);
128 
135  virtual tStdXCombPtr generate(unsigned int sel=0);
137 
138 
145  void persistentOutput(PersistentOStream & os) const;
146 
152  void persistentInput(PersistentIStream & is, int version);
154 
161  static void Init();
162 
167  virtual void initialize();
168 
173  virtual void finalize();
174 
179  virtual void clean();
180 
184  void statistics() const;
185 
191  int statLevel() const {return eventHandler()->statLevel();}
192 
196  CrossSection hardXSec() const { return hardXSec_; }
197 
201  CrossSection softXSec() const { return softXSec_; }
202 
206  CrossSection inelasticXSec() const { return inelXSec_; }
207 
212  CrossSection nonDiffractiveXSec() const {
213  return (1.-diffratio_)*inelXSec_;
214  }
215 
222  CrossSection diffractiveXSec() const {
223  return diffratio_*inelXSec_;
224  }
225 
226 
229 
238  tEHPtr eventHandler() const {return theHandler;}
239 
243  static const MPIHandler * currentHandler() {
244  return currentHandler_;
245  }
246 
250  virtual int Algorithm() const {return algorithm_;}
251 
255  virtual Energy Ptmin() const {
256  if(Ptmin_ > ZERO)
257  return Ptmin_;
258  else
259  throw Exception() << "MPIHandler::Ptmin called without initialize before"
261  }
262 
266  virtual InvEnergy2 beta() const {
267  if(beta_ != 100.0/GeV2)
268  return beta_;
269  else
270  throw Exception() << "MPIHandler::beta called without initialization"
272  }
273 
278  virtual Energy PtForVeto() const {return PtOfQCDProc_;}
279 
284  virtual unsigned int additionalHardProcs() const {return numSubProcs_-1;}
285 
290  virtual double colourDisrupt() const {return colourDisrupt_;}
291 
292 protected:
293 
300  virtual IBPtr clone() const;
301 
306  virtual IBPtr fullclone() const;
308 
309 private:
310 
314  const SubHandlerList & subProcesses()
315  const {return theSubProcesses;}
316 
320  SubHandlerList & subProcesses() {return theSubProcesses;}
321 
325  const CutsList & cuts() const {return theCuts;}
326 
330  CutsList & cuts() {return theCuts;}
331 
335  const ProcessHandlerList & processHandlers()
336  const {return theProcessHandlers;}
337 
341  ProcessHandlerList & processHandlers() {return theProcessHandlers;}
342 
343 
348  void Probs(XSVector UEXSecs);
349 
354  void MultDistribution(string filename) const;
355 
364  InvArea OverlapFunction(Length b, Energy2 mu2=ZERO) const;
365 
370  double poisson(Length b, CrossSection sigma,
371  unsigned int N, Energy2 mu2=ZERO) const;
372 
378  CrossSection totalXSecExp() const;
379 
386  CrossSection totalXSecDiff(CrossSection softXSec,
387  Energy2 softMu2=ZERO) const;
388 
395  InvEnergy2 slopeDiff(CrossSection softXSec,
396  Energy2 softMu2=ZERO) const;
397 
403  InvEnergy2 slopeExp() const;
404 
405 
409  void overrideUECuts();
410 
411 
412 private:
413 
418  MPIHandler & operator=(const MPIHandler &) = delete;
419 
427 
431  SubHandlerList theSubProcesses;
432 
436  CutsList theCuts;
437 
441  ProcessHandlerList theProcessHandlers;
442 
448 
454  unsigned int softMult_;
455 
460 
469 
475  Energy PtOfQCDProc_;
476 
480  Energy Ptmin_;
481 
485  CrossSection hardXSec_;
486 
490  CrossSection softXSec_;
491 
495  CrossSection inelXSec_;
496 
503  CrossSection totalXSecExp_;
504 
509  Energy2 softMu2_;
510 
515  InvEnergy2 beta_;
521 
525  Energy2 invRadius_;
526 
530  unsigned int numSubProcs_;
531 
539 
544  bool softInt_;
545 
550  bool twoComp_;
551 
556  unsigned int DLmode_;
557 
561  double avgNhard_;
562 
566  double avgNsoft_;
567 
572 
577  unsigned int energyExtrapolation_;
578 
582  Energy EEparamA_;
583  Energy EEparamB_;
584  Energy refScale_;
585  Energy pT0_;
586  double b_;
587  Energy offset_;
588 
592  double diffratio_=0.2;
593 
594 protected:
595 
603  class InitError: public Exception {};
604 
607 };
608 
609 }
610 
611 namespace Herwig {
612 
619  struct slopeAndTotalXSec : public GSLHelper<CrossSection, CrossSection> {
620 
621  public:
622 
626  slopeAndTotalXSec(tcMPIHPtr handler): handler_(handler) {}
627 
629  typedef Energy2 ArgType2;
630 
632  typedef InvEnergy2 ValType2;
633 
638  CrossSection f1(ArgType softXSec, ArgType2 softMu2) const {
639  return handler_->totalXSecDiff(softXSec, softMu2);
640  }
641 
646  InvEnergy2 f2(ArgType softXSec, ArgType2 softMu2) const {
647  return handler_->slopeDiff(softXSec, softMu2);
648  }
649 
651  virtual ValType vUnit() const {return 1.0*millibarn;}
652 
654  virtual ArgType aUnit() const {return 1.0*millibarn;}
655 
657  ValType2 vUnit2() const {return 1.0/GeV2;}
658 
660  ArgType2 aUnit2() const {return GeV2;}
661 
662  private:
663 
667  tcMPIHPtr handler_;
668 
669  };
670 
675  struct betaBisection : public GSLHelper<Energy2, InvEnergy2>{
676  public:
684  betaBisection(CrossSection soft, DiffXSec dsig, Energy ptmin)
685  : softXSec_(soft), dsig_(dsig), ptmin_(ptmin) {}
686 
690  virtual Energy2 operator ()(InvEnergy2 beta) const
691  {
692  if( fabs(beta*GeV2) < 1.E-4 )
693  beta = (beta > ZERO) ? 1.E-4/GeV2 : -1.E-4/GeV2;
694 
695  return (exp(beta*sqr(ptmin_)) - 1.0)/beta - softXSec_/dsig_;
696  }
697 
699  virtual ValType vUnit() const {return 1.0*GeV2;}
700 
702  virtual ArgType aUnit() const {return 1.0/GeV2;}
703 
704  private:
705 
707  CrossSection softXSec_;
708 
710  DiffXSec dsig_;
711 
713  Energy ptmin_;
714  };
715 
721  struct slopeBisection : public GSLHelper<InvEnergy2, Energy2> {
722  public:
724  slopeBisection(tcMPIHPtr handler) : handler_(handler) {}
725 
731  InvEnergy2 operator ()(Energy2 arg) const;
732 
734  CrossSection softXSec() const {return softXSec_;}
735 
736  private:
738  tcMPIHPtr handler_;
740  mutable CrossSection softXSec_;
741  };
742 
748  struct TotalXSecBisection : public GSLHelper<CrossSection, CrossSection> {
749  public:
750 
756  TotalXSecBisection(tcMPIHPtr handler, Energy2 softMu2=ZERO):
757  handler_(handler), softMu2_(softMu2) {}
758 
763  CrossSection operator ()(CrossSection argument) const {
764  return handler_->totalXSecDiff(argument, softMu2_);
765  }
766 
768  virtual ValType vUnit() const {return 1.0*millibarn;}
769 
771  virtual ArgType aUnit() const {return 1.0*millibarn;}
772 
773  private:
774 
778  tcMPIHPtr handler_;
779 
783  Energy2 softMu2_;
784 
785  };
786 
790  typedef decltype(mm/GeV2) LengthDiff;
791 
795  struct slopeInt : public GSLHelper<LengthDiff, Length>{
796 
797  public:
804  slopeInt(tcMPIHPtr handler, CrossSection hard,
805  CrossSection soft=0*millibarn, Energy2 softMu2=ZERO)
806  : handler_(handler), hardXSec_(hard),
807  softXSec_(soft), softMu2_(softMu2) {}
808 
813  ValType operator ()(ArgType arg) const;
814 
815  private:
816 
820  tcMPIHPtr handler_;
821 
825  CrossSection hardXSec_;
826 
830  CrossSection softXSec_;
831 
835  Energy2 softMu2_;
836 
837  };
838 
842  struct Eikonalization : public GSLHelper<Length, Length>{
843 
858  Eikonalization(tcMPIHPtr handler, int option, CrossSection hard,
859  CrossSection soft=0*millibarn, Energy2 softMu2=ZERO)
860  : theHandler(handler), theoption(option), hardXSec_(hard),
861  softXSec_(soft), softMu2_(softMu2) {}
862 
866  Length operator ()(Length argument) const;
867 
868  private:
872  tcMPIHPtr theHandler;
873 
879 
883  CrossSection hardXSec_;
884 
888  CrossSection softXSec_;
889 
893  Energy2 softMu2_;
894 
895  };
896 }
897 
898 #endif /* HERWIG_MPIHandler_H */
virtual ValType vUnit() const
provide the actual units of use
Definition: MPIHandler.h:768
Eikonalization(tcMPIHPtr handler, int option, CrossSection hard, CrossSection soft=0 *millibarn, Energy2 softMu2=ZERO)
The constructor.
Definition: MPIHandler.h:858
static const MPIHandler * currentHandler()
Return the current handler.
Definition: MPIHandler.h:243
tEHPtr theHandler
A pointer to the EventHandler that calls us.
Definition: MPIHandler.h:426
bool softInt_
Flag to store whether soft interactions, i.e.
Definition: MPIHandler.h:544
decltype(mm/GeV2) typedef LengthDiff
Typedef for derivative of the length.
Definition: MPIHandler.h:790
unsigned int softMult_
Variable to store the soft multiplicity generated for a event.
Definition: MPIHandler.h:454
vector< int > additionalMultiplicities_
Variable to store the multiplicity of the second hard process.
Definition: MPIHandler.h:459
double avgNsoft_
Variable to store the average soft multiplicity.
Definition: MPIHandler.h:566
CutsList & cuts()
Access the list of cuts.
Definition: MPIHandler.h:330
ArgType2 aUnit2() const
otherwise rounding errors may get significant
Definition: MPIHandler.h:660
ProcessHandlerList theProcessHandlers
List of ProcessHandler used to sample different processes independently.
Definition: MPIHandler.h:441
Energy2 ArgType2
second argument type
Definition: MPIHandler.h:629
virtual ~MPIHandler()
The destructor.
Definition: MPIHandler.h:98
tcMPIHPtr theHandler
Pointer to the Handler that calls this integrand.
Definition: MPIHandler.h:872
CrossSection ValType
Typedef for Value type.
Definition: GSLHelper.h:54
int identicalToUE_
Variable to store the information, which process is identical to the UE one (QCD dijets).
Definition: MPIHandler.h:468
CrossSection softXSec_
soft cross section
Definition: MPIHandler.h:707
Abstract base class used to minimize the dependence between MPIHandler and all Shower classes...
Definition: UEBase.h:25
const CutsList & cuts() const
Access the list of cuts.
Definition: MPIHandler.h:325
SubHandlerList & subProcesses()
Access the list of sub-process handlers.
Definition: MPIHandler.h:320
betaBisection(CrossSection soft, DiffXSec dsig, Energy ptmin)
Constructor.
Definition: MPIHandler.h:684
Energy ptmin_
pt cutoff
Definition: MPIHandler.h:713
CrossSection hardXSec_
The hard cross section to be eikonalized.
Definition: MPIHandler.h:825
A struct for the eikonalization of the inclusive cross section.
Definition: MPIHandler.h:842
slopeBisection(tcMPIHPtr handler)
Constructor.
Definition: MPIHandler.h:724
slopeAndTotalXSec(tcMPIHPtr handler)
Constructor.
Definition: MPIHandler.h:626
virtual Energy Ptmin() const
Return the ptmin parameter of the model.
Definition: MPIHandler.h:255
Energy2 softMu2_
The inv radius^2 of the soft interactions.
Definition: MPIHandler.h:893
A struct for the root finding that is necessary to determine the slope of the soft pt spectrum to mat...
Definition: MPIHandler.h:675
ThePEG::Ptr< InterfacedBase >::pointer IBPtr
double avgNhard_
Variable to store the average hard multiplicity.
Definition: MPIHandler.h:561
from, which can then be used by the GSL algorithms that need a pointer to a function and don&#39;t know a...
Definition: GSLHelper.h:31
constexpr auto sqr(const T &x) -> decltype(x *x)
tcMPIHPtr handler_
Pointer to the Handler that calls this integrand.
Definition: MPIHandler.h:820
tcMPIHPtr handler_
Pointer to the handler.
Definition: MPIHandler.h:667
double colourDisrupt_
Variable to store the relative number of colour disrupted connections to additional subprocesses...
Definition: MPIHandler.h:538
CrossSection f1(ArgType softXSec, ArgType2 softMu2) const
first element of the vector like function to find root for
Definition: MPIHandler.h:638
Energy2 softMu2_
The inv radius^2 of the soft interactions.
Definition: MPIHandler.h:835
tEHPtr eventHandler() const
Return the ThePEG::EventHandler assigned to this handler.
Definition: MPIHandler.h:238
CrossSection totalXSecExp_
Variable to store the total pp cross section (assuming rho=0!) as measured at LHC.
Definition: MPIHandler.h:503
virtual ArgType aUnit() const
otherwise rounding errors may get significant
Definition: MPIHandler.h:654
virtual ArgType aUnit() const
otherwise rounding errors may get significant
Definition: MPIHandler.h:771
A struct for the root finding that is necessary to determine the soft cross section and soft mu2 that...
Definition: MPIHandler.h:721
A struct for the root finding that is necessary to determine the soft cross section that is needed to...
Definition: MPIHandler.h:748
slopeInt(tcMPIHPtr handler, CrossSection hard, CrossSection soft=0 *millibarn, Energy2 softMu2=ZERO)
Constructor.
Definition: MPIHandler.h:804
CrossSection softXSec() const
Return the soft cross section that has been calculated.
Definition: MPIHandler.h:734
ValType2 vUnit2() const
provide the actual units of use
Definition: MPIHandler.h:657
CrossSection softXSec() const
Return the soft cross section below ptmin.
Definition: MPIHandler.h:201
virtual ValType vUnit() const
provide the actual units of use
Definition: MPIHandler.h:699
A struct for the integrand for the slope.
Definition: MPIHandler.h:795
ThePEG::Ptr< EventHandler >::transient_pointer tEHPtr
DiffXSec dsig_
dsigma/dp_T^2 at ptmin
Definition: MPIHandler.h:710
vector< CrossSection > XSVector
A vector of cross sections.
Definition: MPIHandler.h:73
static MPIHandler * currentHandler_
The current handler.
Definition: MPIHandler.h:571
virtual int Algorithm() const
Return theAlgorithm.
Definition: MPIHandler.h:250
virtual ArgType aUnit() const
provide the actual units of use
Definition: MPIHandler.h:702
virtual unsigned int softMultiplicity() const
Get the soft multiplicity from the pretabulated multiplicity distribution.
Definition: MPIHandler.h:121
int theoption
A flag to switch between the calculation of total and inelastic cross section or calculations for the...
Definition: MPIHandler.h:878
Energy2 invRadius_
Inverse hadron Radius squared .
Definition: MPIHandler.h:525
CrossSection inelXSec_
Variable to store the inelastic cross section.
Definition: MPIHandler.h:495
CrossSection nonDiffractiveXSec() const
Return the non-diffractive cross section assumed by the model.
Definition: MPIHandler.h:212
ProcessHandlerList & processHandlers()
Access the list of sub-process handlers.
Definition: MPIHandler.h:341
pair< unsigned int, unsigned int > MPair
A pair of multiplicities: hard, soft.
Definition: MPIHandler.h:76
Energy Ptmin_
Variable to store the parameter ptmin.
Definition: MPIHandler.h:480
InvEnergy2 f2(ArgType softXSec, ArgType2 softMu2) const
second element of the vector like function to find root for
Definition: MPIHandler.h:646
virtual Energy PtForVeto() const
Return the pt Cutoff of the Interaction that is identical to the UE one.
Definition: MPIHandler.h:278
virtual double colourDisrupt() const
Return the fraction of colour disrupted connections to the suprocesses.
Definition: MPIHandler.h:290
int algorithm_
Switch to be set from outside to determine the algorithm used for UE activity.
Definition: MPIHandler.h:520
virtual unsigned int additionalHardProcs() const
Return the number of additional "hard" processes ( = multiple parton scattering)
Definition: MPIHandler.h:284
SubHandlerList theSubProcesses
The list of SubProcessHandlers.
Definition: MPIHandler.h:431
virtual ValType vUnit() const
provide the actual units of use
Definition: MPIHandler.h:651
pair< double, double > generate(const Generator< Density > &gen, double r)
Generate a random variable and return its weight.
const ProcessHandlerList & processHandlers() const
Access the list of sub-process handlers.
Definition: MPIHandler.h:335
unsigned int DLmode_
Switch to determine which Donnachie & Landshoff parametrization should be used.
Definition: MPIHandler.h:556
CrossSection diffractiveXSec() const
Return the diffractive cross section assumed by the model.
Definition: MPIHandler.h:222
CutsList theCuts
The kinematical cuts used for this collision handler.
Definition: MPIHandler.h:436
tcMPIHPtr handler_
The handler.
Definition: MPIHandler.h:778
-*- C++ -*-
int statLevel() const
The level of statistics.
Definition: MPIHandler.h:191
This class is responsible for generating additional semi hard partonic interactions.
Definition: MPIHandler.h:45
tcMPIHPtr handler_
const pointer to the MPIHandler to give access to member functions.
Definition: MPIHandler.h:738
virtual InvEnergy2 beta() const
Return the slope of the soft pt spectrum as calculated.
Definition: MPIHandler.h:266
unsigned int numSubProcs_
Member variable to store the actual number of separate SubProcesses.
Definition: MPIHandler.h:530
MPIHandler()
The default constructor.
Definition: MPIHandler.h:83
virtual bool softInt() const
Return true or false depending on whether soft interactions are enabled.
Definition: MPIHandler.h:114
CrossSection hardXSec() const
Return the hard cross section above ptmin.
Definition: MPIHandler.h:196
CrossSection ArgType
Typedef for Agrument type.
Definition: GSLHelper.h:49
CrossSection softXSec_
The soft cross section to be eikonalized.
Definition: MPIHandler.h:888
const SubHandlerList & subProcesses() const
Access the list of sub-process handlers.
Definition: MPIHandler.h:314
unsigned int energyExtrapolation_
Flag to store whether to calculate the minimal UE pt according to an extrapolation formula or whether...
Definition: MPIHandler.h:577
CrossSection softXSec_
soft cross section that is determined on the fly.
Definition: MPIHandler.h:740
Energy2 softMu2_
Variable to store the soft radius, that is calculated during initialization for the two-component mod...
Definition: MPIHandler.h:509
vector< CutsPtr > CutsList
A vector of Cuts.
Definition: MPIHandler.h:67
InvEnergy2 ValType2
second value type
Definition: MPIHandler.h:632
CrossSection inelasticXSec() const
Return the inelastic cross section.
Definition: MPIHandler.h:206
Energy EEparamA_
Parameters for the energy extrapolation formula.
Definition: MPIHandler.h:582
Energy PtOfQCDProc_
Variable to store the minimal pt of the process that is identical to the UE one.
Definition: MPIHandler.h:475
A struct for the 2D root finding that is necessary to determine the soft cross section and the soft r...
Definition: MPIHandler.h:619
TotalXSecBisection(tcMPIHPtr handler, Energy2 softMu2=ZERO)
Constructor.
Definition: MPIHandler.h:756
vector< SubHdlPtr > SubHandlerList
A vector of SubProcessHandlers.
Definition: MPIHandler.h:64
constexpr ZeroUnit ZERO
vector< ProHdlPtr > ProcessHandlerList
A vector of ProcessHandlers.
Definition: MPIHandler.h:70
Selector< MPair > theMultiplicities
A ThePEG::Selector where the individual Probabilities P_N are stored and the actual Multiplicities ca...
Definition: MPIHandler.h:447
InvEnergy2 beta_
slope to the non-perturbative pt spectrum: .
Definition: MPIHandler.h:515
CrossSection hardXSec_
Variable to store the hard cross section above ptmin.
Definition: MPIHandler.h:485
CrossSection softXSec_
Variable to store the final soft cross section below ptmin.
Definition: MPIHandler.h:490
CrossSection hardXSec_
The hard cross section to be eikonalized.
Definition: MPIHandler.h:883
bool twoComp_
Flag to steer wheather the soft part has a different radius, that will be dynamically fixed...
Definition: MPIHandler.h:550
CrossSection softXSec_
The soft cross section to be eikonalized.
Definition: MPIHandler.h:830