herwig is hosted by Hepforge, IPPP Durham
Herwig  7.1.5
MPIHandler.h
1 // -*- C++ -*-
2 //
3 // MPIHandler.h is a part of Herwig - A multi-purpose Monte Carlo event generator
4 // Copyright (C) 2002-2017 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_(2), EEparamA_(0.6*GeV),
92  EEparamB_(37.5*GeV), refScale_(7000.*GeV),
93  pT0_(3.11*GeV), b_(0.21) {}
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 
210 
219  tEHPtr eventHandler() const {return theHandler;}
220 
224  static const MPIHandler * currentHandler() {
225  return currentHandler_;
226  }
227 
231  virtual int Algorithm() const {return algorithm_;}
232 
236  virtual Energy Ptmin() const {
237  if(Ptmin_ > ZERO)
238  return Ptmin_;
239  else
240  throw Exception() << "MPIHandler::Ptmin called without initialize before"
242  }
243 
247  virtual InvEnergy2 beta() const {
248  if(beta_ != 100.0/GeV2)
249  return beta_;
250  else
251  throw Exception() << "MPIHandler::beta called without initialization"
253  }
254 
259  virtual Energy PtForVeto() const {return PtOfQCDProc_;}
260 
265  virtual unsigned int additionalHardProcs() const {return numSubProcs_-1;}
266 
271  virtual double colourDisrupt() const {return colourDisrupt_;}
272 
273 protected:
274 
281  virtual IBPtr clone() const;
282 
287  virtual IBPtr fullclone() const;
289 
290 private:
291 
295  const SubHandlerList & subProcesses()
296  const {return theSubProcesses;}
297 
301  SubHandlerList & subProcesses() {return theSubProcesses;}
302 
306  const CutsList & cuts() const {return theCuts;}
307 
311  CutsList & cuts() {return theCuts;}
312 
316  const ProcessHandlerList & processHandlers()
317  const {return theProcessHandlers;}
318 
322  ProcessHandlerList & processHandlers() {return theProcessHandlers;}
323 
324 
329  void Probs(XSVector UEXSecs);
330 
335  void MultDistribution(string filename) const;
336 
345  InvArea OverlapFunction(Length b, Energy2 mu2=ZERO) const;
346 
351  double poisson(Length b, CrossSection sigma,
352  unsigned int N, Energy2 mu2=ZERO) const;
353 
357  double factorial (unsigned int n) const;
358 
364  CrossSection totalXSecExp() const;
365 
372  CrossSection totalXSecDiff(CrossSection softXSec,
373  Energy2 softMu2=ZERO) const;
374 
381  InvEnergy2 slopeDiff(CrossSection softXSec,
382  Energy2 softMu2=ZERO) const;
383 
389  InvEnergy2 slopeExp() const;
390 
391 
395  void overrideUECuts();
396 
397 
398 private:
399 
405 
410  MPIHandler & operator=(const MPIHandler &) = delete;
411 
419 
423  SubHandlerList theSubProcesses;
424 
428  CutsList theCuts;
429 
433  ProcessHandlerList theProcessHandlers;
434 
440 
446  unsigned int softMult_;
447 
452 
461 
467  Energy PtOfQCDProc_;
468 
472  Energy Ptmin_;
473 
477  CrossSection hardXSec_;
478 
482  CrossSection softXSec_;
483 
487  CrossSection inelXSec_;
488 
495  CrossSection totalXSecExp_;
496 
501  Energy2 softMu2_;
502 
507  InvEnergy2 beta_;
513 
517  Energy2 invRadius_;
518 
522  unsigned int numSubProcs_;
523 
531 
536  bool softInt_;
537 
542  bool twoComp_;
543 
548  unsigned int DLmode_;
549 
553  double avgNhard_;
554 
558  double avgNsoft_;
559 
564 
569  unsigned int energyExtrapolation_;
570 
574  Energy EEparamA_;
575  Energy EEparamB_;
576  Energy refScale_;
577  Energy pT0_;
578  double b_;
579 
580 protected:
581 
589  class InitError: public Exception {};
590 
593 };
594 
595 }
596 
597 #include "ThePEG/Utilities/ClassTraits.h"
598 
599 namespace ThePEG {
600 
605 template <>
606 struct BaseClassTrait<Herwig::MPIHandler,1> {
608  typedef Interfaced NthBase;
609 };
610 
613 template <>
614 struct ClassTraits<Herwig::MPIHandler>
615  : public ClassTraitsBase<Herwig::MPIHandler> {
617  static string className() { return "Herwig::MPIHandler"; }
621  static string library() { return "JetCuts.so SimpleKTCut.so HwMPI.so"; }
622 };
623 
626 }
627 
628 namespace Herwig {
629 
636  struct slopeAndTotalXSec : public GSLHelper<CrossSection, CrossSection> {
637 
638  public:
639 
643  slopeAndTotalXSec(tcMPIHPtr handler): handler_(handler) {}
644 
646  typedef Energy2 ArgType2;
647 
649  typedef InvEnergy2 ValType2;
650 
655  CrossSection f1(ArgType softXSec, ArgType2 softMu2) const {
656  return handler_->totalXSecDiff(softXSec, softMu2);
657  }
658 
663  InvEnergy2 f2(ArgType softXSec, ArgType2 softMu2) const {
664  return handler_->slopeDiff(softXSec, softMu2);
665  }
666 
668  virtual ValType vUnit() const {return 1.0*millibarn;}
669 
671  virtual ArgType aUnit() const {return 1.0*millibarn;}
672 
674  ValType2 vUnit2() const {return 1.0/GeV2;}
675 
677  ArgType2 aUnit2() const {return GeV2;}
678 
679  private:
680 
684  tcMPIHPtr handler_;
685 
686  };
687 
692  struct betaBisection : public GSLHelper<Energy2, InvEnergy2>{
693  public:
701  betaBisection(CrossSection soft, DiffXSec dsig, Energy ptmin)
702  : softXSec_(soft), dsig_(dsig), ptmin_(ptmin) {}
703 
707  virtual Energy2 operator ()(InvEnergy2 beta) const
708  {
709  if( fabs(beta*GeV2) < 1.E-4 )
710  beta = (beta > ZERO) ? 1.E-4/GeV2 : -1.E-4/GeV2;
711 
712  return (exp(beta*sqr(ptmin_)) - 1.0)/beta - softXSec_/dsig_;
713  }
714 
716  virtual ValType vUnit() const {return 1.0*GeV2;}
717 
719  virtual ArgType aUnit() const {return 1.0/GeV2;}
720 
721  private:
722 
724  CrossSection softXSec_;
725 
727  DiffXSec dsig_;
728 
730  Energy ptmin_;
731  };
732 
738  struct slopeBisection : public GSLHelper<InvEnergy2, Energy2> {
739  public:
741  slopeBisection(tcMPIHPtr handler) : handler_(handler) {}
742 
748  InvEnergy2 operator ()(Energy2 arg) const;
749 
751  CrossSection softXSec() const {return softXSec_;}
752 
753  private:
755  tcMPIHPtr handler_;
757  mutable CrossSection softXSec_;
758  };
759 
765  struct TotalXSecBisection : public GSLHelper<CrossSection, CrossSection> {
766  public:
767 
773  TotalXSecBisection(tcMPIHPtr handler, Energy2 softMu2=ZERO):
774  handler_(handler), softMu2_(softMu2) {}
775 
780  CrossSection operator ()(CrossSection argument) const {
781  return handler_->totalXSecDiff(argument, softMu2_);
782  }
783 
785  virtual ValType vUnit() const {return 1.0*millibarn;}
786 
788  virtual ArgType aUnit() const {return 1.0*millibarn;}
789 
790  private:
791 
795  tcMPIHPtr handler_;
796 
800  Energy2 softMu2_;
801 
802  };
803 
807  typedef Qty<1,-2,0> LengthDiff;
808 
812  struct slopeInt : public GSLHelper<LengthDiff, Length>{
813 
814  public:
821  slopeInt(tcMPIHPtr handler, CrossSection hard,
822  CrossSection soft=0*millibarn, Energy2 softMu2=ZERO)
823  : handler_(handler), hardXSec_(hard),
824  softXSec_(soft), softMu2_(softMu2) {}
825 
830  ValType operator ()(ArgType arg) const;
831 
832  private:
833 
837  tcMPIHPtr handler_;
838 
842  CrossSection hardXSec_;
843 
847  CrossSection softXSec_;
848 
852  Energy2 softMu2_;
853 
854  };
855 
859  struct Eikonalization : public GSLHelper<Length, Length>{
860 
875  Eikonalization(tcMPIHPtr handler, int option, CrossSection hard,
876  CrossSection soft=0*millibarn, Energy2 softMu2=ZERO)
877  : theHandler(handler), theoption(option), hardXSec_(hard),
878  softXSec_(soft), softMu2_(softMu2) {}
879 
883  Length operator ()(Length argument) const;
884 
885  private:
889  tcMPIHPtr theHandler;
890 
896 
900  CrossSection hardXSec_;
901 
905  CrossSection softXSec_;
906 
910  Energy2 softMu2_;
911 
912  };
913 }
914 
915 
916 #ifndef ThePEG_TEMPLATES_IN_CC_FILE
917 // #include "MPIHandler.tcc"
918 #endif
919 
920 #endif /* HERWIG_MPIHandler_H */
virtual ValType vUnit() const
provide the actual units of use
Definition: MPIHandler.h:785
Eikonalization(tcMPIHPtr handler, int option, CrossSection hard, CrossSection soft=0 *millibarn, Energy2 softMu2=ZERO)
The constructor.
Definition: MPIHandler.h:875
static const MPIHandler * currentHandler()
Return the current handler.
Definition: MPIHandler.h:224
tEHPtr theHandler
A pointer to the EventHandler that calls us.
Definition: MPIHandler.h:418
bool softInt_
Flag to store whether soft interactions, i.e.
Definition: MPIHandler.h:536
unsigned int softMult_
Variable to store the soft multiplicity generated for a event.
Definition: MPIHandler.h:446
vector< int > additionalMultiplicities_
Variable to store the multiplicity of the second hard process.
Definition: MPIHandler.h:451
double avgNsoft_
Variable to store the average soft multiplicity.
Definition: MPIHandler.h:558
CutsList & cuts()
Access the list of cuts.
Definition: MPIHandler.h:311
ArgType2 aUnit2() const
otherwise rounding errors may get significant
Definition: MPIHandler.h:677
ProcessHandlerList theProcessHandlers
List of ProcessHandler used to sample different processes independently.
Definition: MPIHandler.h:433
Energy2 ArgType2
second argument type
Definition: MPIHandler.h:646
virtual ~MPIHandler()
The destructor.
Definition: MPIHandler.h:98
tcMPIHPtr theHandler
Pointer to the Handler that calls this integrand.
Definition: MPIHandler.h:889
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:460
CrossSection softXSec_
soft cross section
Definition: MPIHandler.h:724
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:306
SubHandlerList & subProcesses()
Access the list of sub-process handlers.
Definition: MPIHandler.h:301
betaBisection(CrossSection soft, DiffXSec dsig, Energy ptmin)
Constructor.
Definition: MPIHandler.h:701
Energy ptmin_
pt cutoff
Definition: MPIHandler.h:730
CrossSection hardXSec_
The hard cross section to be eikonalized.
Definition: MPIHandler.h:842
A struct for the eikonalization of the inclusive cross section.
Definition: MPIHandler.h:859
slopeBisection(tcMPIHPtr handler)
Constructor.
Definition: MPIHandler.h:741
slopeAndTotalXSec(tcMPIHPtr handler)
Constructor.
Definition: MPIHandler.h:643
virtual Energy Ptmin() const
Return the ptmin parameter of the model.
Definition: MPIHandler.h:236
Energy2 softMu2_
The inv radius^2 of the soft interactions.
Definition: MPIHandler.h:910
A struct for the root finding that is necessary to determine the slope of the soft pt spectrum to mat...
Definition: MPIHandler.h:692
ThePEG::Ptr< InterfacedBase >::pointer IBPtr
double avgNhard_
Variable to store the average hard multiplicity.
Definition: MPIHandler.h:553
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
tcMPIHPtr handler_
Pointer to the Handler that calls this integrand.
Definition: MPIHandler.h:837
tcMPIHPtr handler_
Pointer to the handler.
Definition: MPIHandler.h:684
double colourDisrupt_
Variable to store the relative number of colour disrupted connections to additional subprocesses...
Definition: MPIHandler.h:530
CrossSection f1(ArgType softXSec, ArgType2 softMu2) const
first element of the vector like function to find root for
Definition: MPIHandler.h:655
Energy2 softMu2_
The inv radius^2 of the soft interactions.
Definition: MPIHandler.h:852
tEHPtr eventHandler() const
Return the ThePEG::EventHandler assigned to this handler.
Definition: MPIHandler.h:219
CrossSection totalXSecExp_
Variable to store the total pp cross section (assuming rho=0!) as measured at LHC.
Definition: MPIHandler.h:495
virtual ArgType aUnit() const
otherwise rounding errors may get significant
Definition: MPIHandler.h:671
virtual ArgType aUnit() const
otherwise rounding errors may get significant
Definition: MPIHandler.h:788
A struct for the root finding that is necessary to determine the soft cross section and soft mu2 that...
Definition: MPIHandler.h:738
A struct for the root finding that is necessary to determine the soft cross section that is needed to...
Definition: MPIHandler.h:765
slopeInt(tcMPIHPtr handler, CrossSection hard, CrossSection soft=0 *millibarn, Energy2 softMu2=ZERO)
Constructor.
Definition: MPIHandler.h:821
CrossSection softXSec() const
Return the soft cross section that has been calculated.
Definition: MPIHandler.h:751
ValType2 vUnit2() const
provide the actual units of use
Definition: MPIHandler.h:674
CrossSection softXSec() const
Return the soft cross section below ptmin.
Definition: MPIHandler.h:201
static ClassDescription< MPIHandler > initMPIHandler
The static object used to initialize the description of this class.
Definition: MPIHandler.h:404
constexpr BinaryOpTraits< T, T >::MulT sqr(const T &x)
virtual ValType vUnit() const
provide the actual units of use
Definition: MPIHandler.h:716
A struct for the integrand for the slope.
Definition: MPIHandler.h:812
ThePEG::Ptr< EventHandler >::transient_pointer tEHPtr
DiffXSec dsig_
dsigma/dp_T^2 at ptmin
Definition: MPIHandler.h:727
vector< CrossSection > XSVector
A vector of cross sections.
Definition: MPIHandler.h:73
static MPIHandler * currentHandler_
The current handler.
Definition: MPIHandler.h:563
virtual int Algorithm() const
Return theAlgorithm.
Definition: MPIHandler.h:231
virtual ArgType aUnit() const
provide the actual units of use
Definition: MPIHandler.h:719
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:895
Energy2 invRadius_
Inverse hadron Radius squared .
Definition: MPIHandler.h:517
Qty< 1,-2, 0 > LengthDiff
Typedef for derivative of the length.
Definition: MPIHandler.h:807
CrossSection inelXSec_
Variable to store the inelastic cross section.
Definition: MPIHandler.h:487
ProcessHandlerList & processHandlers()
Access the list of sub-process handlers.
Definition: MPIHandler.h:322
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:472
InvEnergy2 f2(ArgType softXSec, ArgType2 softMu2) const
second element of the vector like function to find root for
Definition: MPIHandler.h:663
virtual Energy PtForVeto() const
Return the pt Cutoff of the Interaction that is identical to the UE one.
Definition: MPIHandler.h:259
virtual double colourDisrupt() const
Return the fraction of colour disrupted connections to the suprocesses.
Definition: MPIHandler.h:271
int algorithm_
Switch to be set from outside to determine the algorithm used for UE activity.
Definition: MPIHandler.h:512
virtual unsigned int additionalHardProcs() const
Return the number of additional "hard" processes ( = multiple parton scattering)
Definition: MPIHandler.h:265
SubHandlerList theSubProcesses
The list of SubProcessHandlers.
Definition: MPIHandler.h:423
virtual ValType vUnit() const
provide the actual units of use
Definition: MPIHandler.h:668
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:316
unsigned int DLmode_
Switch to determine which Donnachie & Landshoff parametrization should be used.
Definition: MPIHandler.h:548
CutsList theCuts
The kinematical cuts used for this collision handler.
Definition: MPIHandler.h:428
tcMPIHPtr handler_
The handler.
Definition: MPIHandler.h:795
-*- 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:755
virtual InvEnergy2 beta() const
Return the slope of the soft pt spectrum as calculated.
Definition: MPIHandler.h:247
unsigned int numSubProcs_
Member variable to store the actual number of separate SubProcesses.
Definition: MPIHandler.h:522
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:905
const SubHandlerList & subProcesses() const
Access the list of sub-process handlers.
Definition: MPIHandler.h:295
unsigned int energyExtrapolation_
Flag to store whether to calculate the minimal UE pt according to an extrapolation formula or whether...
Definition: MPIHandler.h:569
CrossSection softXSec_
soft cross section that is determined on the fly.
Definition: MPIHandler.h:757
Energy2 softMu2_
Variable to store the soft radius, that is calculated during initialization for the two-component mod...
Definition: MPIHandler.h:501
vector< CutsPtr > CutsList
A vector of Cuts.
Definition: MPIHandler.h:67
constexpr long double factorial(unsigned int n)
InvEnergy2 ValType2
second value type
Definition: MPIHandler.h:649
CrossSection inelasticXSec() const
Return the inelastic cross section.
Definition: MPIHandler.h:206
Energy EEparamA_
Parameters for the energy extrapolation formula.
Definition: MPIHandler.h:574
Energy PtOfQCDProc_
Variable to store the minimal pt of the process that is identical to the UE one.
Definition: MPIHandler.h:467
A struct for the 2D root finding that is necessary to determine the soft cross section and the soft r...
Definition: MPIHandler.h:636
TotalXSecBisection(tcMPIHPtr handler, Energy2 softMu2=ZERO)
Constructor.
Definition: MPIHandler.h:773
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:439
InvEnergy2 beta_
slope to the non-perturbative pt spectrum: .
Definition: MPIHandler.h:507
CrossSection hardXSec_
Variable to store the hard cross section above ptmin.
Definition: MPIHandler.h:477
CrossSection softXSec_
Variable to store the final soft cross section below ptmin.
Definition: MPIHandler.h:482
CrossSection hardXSec_
The hard cross section to be eikonalized.
Definition: MPIHandler.h:900
bool twoComp_
Flag to steer wheather the soft part has a different radius, that will be dynamically fixed...
Definition: MPIHandler.h:542
CrossSection softXSec_
The soft cross section to be eikonalized.
Definition: MPIHandler.h:847