PYTHIA8/pythia8130/src/UserHooks.cxx

   1 // UserHooks.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2008 Torbjorn Sjostrand.
   3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
   4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
   5
   6 // Function definitions (not found in the header) for the UserHooks class.
   7
   8 #include "UserHooks.h"
   9
  10 namespace Pythia8 {
  11
  12 //**************************************************************************
  13
  14 // The UserHooks class.
  15
  16 //*********
  17
  18 // multiplySigmaBy allows the user to introduce a multiplicative factor
  19 // that modifies the cross section of a hard process. Since it is called
  20 // from before the event record is generated in full, the normal analysis
  21 // does not work. The code here provides a rather extensive summary of
  22 // which methods actually do work. It is a convenient starting point for
  23 // writing your own derived routine.
  24
  25 double UserHooks::multiplySigmaBy( const SigmaProcess* sigmaProcessPtr,
  26   const PhaseSpace* phaseSpacePtr, bool inEvent) {
  27
  28   // Process code, necessary when some to be treated differently.
  29   //int code       = sigmaProcessPtr->code();
  30
  31   // Final multiplicity, i.e. whether 2 -> 1 or 2 -> 2.
  32   //int nFinal     = sigmaProcessPtr->nFinal();
  33
  34   // Incoming x1 and x2 to the hard collision, and factorization scale.
  35   //double x1      = phaseSpacePtr->x1();
  36   //double x2      = phaseSpacePtr->x2();
  37   //double Q2Fac   = sigmaProcessPtr->Q2Fac();
  38
  39   // Renormalization scale and assumed alpha_strong and alpha_EM.
  40   //double Q2Ren   = sigmaProcessPtr->Q2Ren();
  41   //double alphaS  = sigmaProcessPtr->alphaSRen();
  42   //double alphaEM = sigmaProcessPtr->alphaEMRen();
  43
  44   // Subprocess mass-square.
  45   //double sHat = phaseSpacePtr->sHat();
  46
  47   // Now methods only relevant for 2 -> 2.
  48   //if (nFinal == 2) {
  49
  50     // Mandelstam variables and hard-process pT.
  51     //double tHat  = phaseSpacePtr->tHat();
  52     //double uHat  = phaseSpacePtr->uHat();
  53     //double pTHat = phaseSpacePtr->pTHat();
  54
  55     // Masses of the final-state particles. (Here 0 for light quarks.)
  56     //double m3    = sigmaProcessPtr->m(3);
  57     //double m4    = sigmaProcessPtr->m(4);
  58   //}
  59
  60   // Dummy statement to avoid compiler warnings.
  61   return ((inEvent && sigmaProcessPtr->code() == 0
  62     && phaseSpacePtr->sHat() < 0.) ? 0. : 1.);
  63
  64 }
  65
  66 //*********
  67
  68 // subEvent extracts currently resolved partons in the hard process.
  69
  70 void UserHooks::subEvent(const Event& event, bool isHardest) {
  71
  72   // Reset work event to be empty.
  73   workEvent.clear();
  74
  75   // Find which subsystem to study.
  76   int iSys = 0;
  77   if (!isHardest) iSys = event.sizeSystems() - 1;
  78
  79   // Loop through all the final partons of the given subsystem.
  80   for (int i = 2; i < event.sizeSystem(iSys); ++i) {
  81     int iOld = event.getInSystem( iSys, i);
  82
  83     // Copy partons to work event.
  84     int iNew = workEvent.append( event[iOld]);
  85
  86     // No mothers. Position in full event as daughters.
  87     workEvent[iNew].mothers( 0, 0);
  88     workEvent[iNew].daughters( iOld, iOld);
  89   }
  90
  91 }
  92
  93 //**************************************************************************
  94
  95 // The SuppressSmallPT class, derived from UserHooks.
  96
  97 //*********
  98
  99 // Modify event weight at the trial level, before selection.
 100
 101 double SuppressSmallPT::multiplySigmaBy( const SigmaProcess* sigmaProcessPtr,
 102   const PhaseSpace* phaseSpacePtr, bool ) {
 103
 104   // Need to initialize first time this method is called.
 105   if (!isInit) {
 106
 107     // Calculate pT0 as for multiple interactions.
 108     // Fudge factor allows offset relative to MI framework.
 109     double eCM    = phaseSpacePtr->ecm();
 110     double pT0Ref = Settings::parm("MultipleInteractions:pT0Ref");
 111     double ecmRef = Settings::parm("MultipleInteractions:ecmRef");
 112     double ecmPow = Settings::parm("MultipleInteractions:ecmPow");
 113     double pT0    = pT0timesMI * pT0Ref * pow(eCM / ecmRef, ecmPow);
 114     pT20          = pT0 * pT0;
 115
 116     // Initialize alpha_strong object as for multiple interactions,
 117     // alternatively as for hard processes.
 118     double alphaSvalue;
 119     int    alphaSorder;
 120     if (useSameAlphaSasMI) {
 121       alphaSvalue = Settings::parm("MultipleInteractions:alphaSvalue");
 122       alphaSorder = Settings::mode("MultipleInteractions:alphaSorder");
 123     } else {
 124       alphaSvalue = Settings::parm("SigmaProcess:alphaSvalue");
 125       alphaSorder = Settings::mode("SigmaProcess:alphaSorder");
 126     }
 127     alphaS.init( alphaSvalue, alphaSorder);
 128
 129     // Initialization finished.
 130     isInit = true;
 131   }
 132
 133   // Only modify 2 -> 2 processes.
 134   int nFinal = sigmaProcessPtr->nFinal();
 135   if (nFinal != 2) return 1.;
 136
 137   // pT scale of process. Weight pT^4 / (pT^2 + pT0^2)^2
 138   double pTHat     = phaseSpacePtr->pTHat();
 139   double pT2       = pTHat * pTHat;
 140   double wt        = pow2( pT2 / (pT20 + pT2) );
 141
 142   if (numberAlphaS > 0) {
 143     // Renormalization scale and assumed alpha_strong.
 144     double Q2RenOld  = sigmaProcessPtr->Q2Ren();
 145     double alphaSOld = sigmaProcessPtr->alphaSRen();
 146
 147     // Reweight to new alpha_strong at new scale.
 148     double Q2RenNew  = pT20 + Q2RenOld;
 149     double alphaSNew = alphaS.alphaS(Q2RenNew);
 150     wt              *= pow( alphaSNew / alphaSOld, numberAlphaS);
 151   }
 152
 153   // End weight calculation.
 154   return wt;
 155
 156 }
 157
 158
 159 //**************************************************************************
 160
 161 } // end namespace Pythia8