PYTHIA8/pythia8145/src/UserHooks.cxx

   1 // UserHooks.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2010 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 // omitResonanceDecays omits resonance decay chains from process record.
  69
  70 void UserHooks::omitResonanceDecays(const Event& process) {
  71
  72   // Reset work event to be empty
  73   workEvent.clear();
  74
  75   // Loop through all partons. Beam particles should be copied.
  76   for (int i = 0; i < process.size(); ++i) {
  77     bool doCopy  = false;
  78     bool isFinal = false;
  79     if (i < 3) doCopy = true;
  80
  81     // Daughters of beams should be copied.
  82     else {
  83       int iMother = process[i].mother1();
  84       if (iMother == 1 || iMother == 2) doCopy = true;
  85
  86       // Granddaughters of beams should be copied and are final.
  87       else if (iMother > 2) {
  88         int iGrandMother =  process[iMother].mother1();
  89         if (iGrandMother == 1 || iGrandMother == 2) {
  90           doCopy  = true;
  91           isFinal = true;
  92         }
  93       }
  94     }
  95
  96     // Do copying and modify status/daughters of final.
  97     if (doCopy) {
  98       int iNew = workEvent.append( process[i]);
  99       if (isFinal) {
 100         workEvent[iNew].statusPos();
 101         workEvent[iNew].daughters( 0, 0);
 102       }
 103     }
 104   }
 105
 106 }
 107
 108 //--------------------------------------------------------------------------
 109
 110 // subEvent extracts currently resolved partons in the hard process.
 111
 112 void UserHooks::subEvent(const Event& event, bool isHardest) {
 113
 114   // Reset work event to be empty.
 115   workEvent.clear();
 116
 117   // Find which subsystem to study.
 118   int iSys = 0;
 119   if (!isHardest) iSys = partonSystemsPtr->sizeSys() - 1;
 120
 121   // Loop through all the final partons of the given subsystem.
 122   for (int i = 0; i < partonSystemsPtr->sizeOut(iSys); ++i) {
 123     int iOld = partonSystemsPtr->getOut( iSys, i);
 124
 125     // Copy partons to work event.
 126     int iNew = workEvent.append( event[iOld]);
 127
 128     // No mothers. Position in full event as daughters.
 129     workEvent[iNew].mothers( 0, 0);
 130     workEvent[iNew].daughters( iOld, iOld);
 131   }
 132
 133 }
 134
 135 //==========================================================================
 136
 137 // The SuppressSmallPT class, derived from UserHooks.
 138
 139 //--------------------------------------------------------------------------
 140
 141 // Modify event weight at the trial level, before selection.
 142
 143 double SuppressSmallPT::multiplySigmaBy( const SigmaProcess* sigmaProcessPtr,
 144   const PhaseSpace* phaseSpacePtr, bool ) {
 145
 146   // Need to initialize first time this method is called.
 147   if (!isInit) {
 148
 149     // Calculate pT0 as for multiple interactions.
 150     // Fudge factor allows offset relative to MI framework.
 151     double eCM    = phaseSpacePtr->ecm();
 152     double pT0Ref = settingsPtr->parm("MultipleInteractions:pT0Ref");
 153     double ecmRef = settingsPtr->parm("MultipleInteractions:ecmRef");
 154     double ecmPow = settingsPtr->parm("MultipleInteractions:ecmPow");
 155     double pT0    = pT0timesMI * pT0Ref * pow(eCM / ecmRef, ecmPow);
 156     pT20          = pT0 * pT0;
 157
 158     // Initialize alpha_strong object as for multiple interactions,
 159     // alternatively as for hard processes.
 160     double alphaSvalue;
 161     int    alphaSorder;
 162     if (useSameAlphaSasMI) {
 163       alphaSvalue = settingsPtr->parm("MultipleInteractions:alphaSvalue");
 164       alphaSorder = settingsPtr->mode("MultipleInteractions:alphaSorder");
 165     } else {
 166       alphaSvalue = settingsPtr->parm("SigmaProcess:alphaSvalue");
 167       alphaSorder = settingsPtr->mode("SigmaProcess:alphaSorder");
 168     }
 169     alphaS.init( alphaSvalue, alphaSorder);
 170
 171     // Initialization finished.
 172     isInit = true;
 173   }
 174
 175   // Only modify 2 -> 2 processes.
 176   int nFinal = sigmaProcessPtr->nFinal();
 177   if (nFinal != 2) return 1.;
 178
 179   // pT scale of process. Weight pT^4 / (pT^2 + pT0^2)^2
 180   double pTHat     = phaseSpacePtr->pTHat();
 181   double pT2       = pTHat * pTHat;
 182   double wt        = pow2( pT2 / (pT20 + pT2) );
 183
 184   if (numberAlphaS > 0) {
 185     // Renormalization scale and assumed alpha_strong.
 186     double Q2RenOld  = sigmaProcessPtr->Q2Ren();
 187     double alphaSOld = sigmaProcessPtr->alphaSRen();
 188
 189     // Reweight to new alpha_strong at new scale.
 190     double Q2RenNew  = pT20 + Q2RenOld;
 191     double alphaSNew = alphaS.alphaS(Q2RenNew);
 192     wt              *= pow( alphaSNew / alphaSOld, numberAlphaS);
 193   }
 194
 195   // End weight calculation.
 196   return wt;
 197
 198 }
 199
 200
 201 //==========================================================================
 202
 203 } // end namespace Pythia8