include/Pythia8Plugins/SetLHEDecayProductHook.h

0001 // SetLHEDecayProductHook.h is part of the PYTHIA event generator.
0002 // Copyright (C) 2024 Stephen Mrenna, Torbjorn Sjostrand.
0003 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
0004 // Please respect the MCnet Guidelines, see GUIDELINES for details.
0005 // Author: Stephen Mrenna, December 2022.
0006
0007 // This class is used to modify the resonance decay products from an
0008 // input LHE file.
0009
0010 // The following settings are available:
0011 // flag SetLHEDecayProduct:filter = false
0012 //   Activate filter if true.
0013
0014 #ifndef Pythia8_SetLHEDecayProductHooks_H
0015 #define Pythia8_SetLHEDecayProductHooks_H
0016
0017 // Includes.
0018 #include "Pythia8/Pythia.h"
0019 #include "Pythia8/UserHooks.h"
0020 #include "Pythia8/Event.h"
0021
0022 namespace Pythia8 {
0023
0024 class SetLHEDecayProductHook : public UserHooks {
0025 public:
0026
0027   // Constructor.
0028   SetLHEDecayProductHook(Settings &settings, const ParticleData* pdtPtrIn);
0029
0030   // Override base class methods.
0031   bool canVetoProcessLevel() override {return true;}
0032   bool doVetoProcessLevel(Event& process) override {
0033     return checkVetoProcessLevel(process);}
0034   bool initAfterBeams() override;
0035
0036   // Class specific.
0037   bool checkVetoProcessLevel(Event& process);
0038   unsigned long int returnCounter() {return counter;};
0039
0040 private:
0041
0042   // Return the custom particle mass.
0043   double m0powheg(const int idIn) {
0044     double mReturn(-1.0);
0045     switch( idIn ) {
0046     case 3: mReturn=0.2; break;
0047     case 4: mReturn=1.5; break;
0048     case 1: mReturn=0.1; break;
0049     case 2: mReturn=0.1; break;
0050     default: mReturn = pdtPtr->m0(idIn); break;
0051     }
0052     return mReturn;
0053   }
0054
0055   // Data members.
0056   bool filter;
0057   const ParticleData* pdtPtr;
0058   unsigned long int counter;
0059
0060 };
0061
0062 //--------------------------------------------------------------------------
0063
0064 // Constructor.
0065
0066 SetLHEDecayProductHook::SetLHEDecayProductHook(Settings &settings,
0067   const ParticleData* pdtPtrIn) :
0068   pdtPtr(pdtPtrIn), counter(0) {
0069   settings.addFlag("SetLHEDecayProduct:filter", false);
0070 }
0071
0072 //--------------------------------------------------------------------------
0073
0074 // Intialize the user hook after the beams.
0075
0076 bool SetLHEDecayProductHook::initAfterBeams() {
0077   filter = settingsPtr->flag("SetLHEDecayProduct:filter");
0078   return true;
0079 }
0080
0081 //--------------------------------------------------------------------------
0082
0083 // Return true if the resonance decays are vetoed.
0084
0085 bool SetLHEDecayProductHook::checkVetoProcessLevel(Event& process) {
0086
0087   if (!filter) return false;
0088   counter++;
0089   // Determine which W decays hadronically
0090   int hadronW = ( rndmPtr->flat() < 0.5 ) ? 1 : 2;
0091   int idQuark = ( rndmPtr->flat() < 0.5 ) ? 1 : 3;
0092   int newCol  = process.nextColTag();
0093
0094   int countW(0);
0095   for (int i = 0; i < process.size(); ++i) {
0096     if (process[i].idAbs() == 24 && process[i].status() == -22) {
0097       countW++;
0098       if( countW > 2 ) break;
0099       int idFermion = ( hadronW == countW ) ? idQuark : 11 ;
0100       // Identify W+- daughters, properly ordered.
0101       int idV = process[i].id();
0102       int i1  = process[i].daughter1();
0103       int i2  = process[i].daughter2();
0104       int id1  = process[i1].idAbs();
0105       int id2  = process[i2].idAbs();
0106       if( id1 == idFermion || id2 == idFermion ) continue;
0107
0108       double mV = process[i].m();
0109       Vec4 p1 = process[i1].p();
0110       p1.bstback( process[i].p() );
0111       double pMag = p1.pAbs();
0112
0113       // The current distributions are for the particle with the
0114       // same charge sign as the mother, i.e. W- -> e-.
0115       if (process[i1].id() * idV > 0) swap( i1, i2);
0116       process[i1].id(idFermion * process[i1].id()/abs(process[i1].id()) );
0117       process[i2].id((idFermion+1) * process[i2].id()/abs(process[i2].id()) );
0118       if( idFermion != 11 ) {
0119         if( process[i1].id() > 0 ) {
0120           process[i1].col(newCol);
0121           process[i2].acol(newCol);
0122
0123         } else {
0124           process[i2].col(newCol);
0125           process[i1].acol(newCol);
0126         }
0127       }
0128
0129       double m1 = m0powheg(process[i1].idAbs());
0130       double m2 = m0powheg(process[i2].idAbs());
0131
0132       // Energy and absolute momentum of first decay product in W rest frame.
0133       double e1 = 0.5* (pow2(mV) + pow2(m1) - pow2(m2))/mV;
0134       double pA = sqrt(pow2(e1) - pow2(m1));
0135
0136       p1.rescale3( pA/ pMag );
0137       p1.e( e1 );
0138       Vec4 p2   = Vec4(0,0,0,mV) - p1;
0139
0140
0141       p1.bst( process[i].p() );
0142       p2.bst( process[i].p() );
0143       process[i1].p( p1 );
0144       process[i2].p( p2 );
0145       process[i1].m( m1 );
0146       process[i2].m( m2 );
0147     }
0148     // End of loop over W's. Do not veto any events.
0149
0150   }
0151   return false;
0152 }
0153
0154 //==========================================================================
0155
0156 } // end namespace Pythia8
0157
0158 #endif // end Pythia8_SetLHEDecayProductHooks_H