PYTHIA  8.312
SetLHEDecayProductHook.h
1 // SetLHEDecayProductHook.h is part of the PYTHIA event generator.
2 // Copyright (C) 2024 Stephen Mrenna, Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 // Author: Stephen Mrenna, December 2022.
6 
7 // This class is used to modify the resonance decay products from an
8 // input LHE file.
9 
10 // The following settings are available:
11 // flag SetLHEDecayProduct:filter = false
12 // Activate filter if true.
13 
14 #ifndef Pythia8_SetLHEDecayProductHooks_H
15 #define Pythia8_SetLHEDecayProductHooks_H
16 
17 // Includes.
18 #include "Pythia8/Pythia.h"
19 #include "Pythia8/UserHooks.h"
20 #include "Pythia8/Event.h"
21 
22 namespace Pythia8 {
23 
25 public:
26 
27  // Constructor.
28  SetLHEDecayProductHook(Settings &settings, const ParticleData* pdtPtrIn);
29 
30  // Override base class methods.
31  bool canVetoProcessLevel() override {return true;}
32  bool doVetoProcessLevel(Event& process) override {
33  return checkVetoProcessLevel(process);}
34  bool initAfterBeams() override;
35 
36  // Class specific.
37  bool checkVetoProcessLevel(Event& process);
38  unsigned long int returnCounter() {return counter;};
39 
40 private:
41 
42  // Return the custom particle mass.
43  double m0powheg(const int idIn) {
44  double mReturn(-1.0);
45  switch( idIn ) {
46  case 3: mReturn=0.2; break;
47  case 4: mReturn=1.5; break;
48  case 1: mReturn=0.1; break;
49  case 2: mReturn=0.1; break;
50  default: mReturn = pdtPtr->m0(idIn); break;
51  }
52  return mReturn;
53  }
54 
55  // Data members.
56  bool filter;
57  const ParticleData* pdtPtr;
58  unsigned long int counter;
59 
60 };
61 
62 //--------------------------------------------------------------------------
63 
64 // Constructor.
65 
67  const ParticleData* pdtPtrIn) :
68  pdtPtr(pdtPtrIn), counter(0) {
69  settings.addFlag("SetLHEDecayProduct:filter", false);
70 }
71 
72 //--------------------------------------------------------------------------
73 
74 // Intialize the user hook after the beams.
75 
77  filter = settingsPtr->flag("SetLHEDecayProduct:filter");
78  return true;
79 }
80 
81 //--------------------------------------------------------------------------
82 
83 // Return true if the resonance decays are vetoed.
84 
86 
87  if (!filter) return false;
88  counter++;
89  // Determine which W decays hadronically
90  int hadronW = ( rndmPtr->flat() < 0.5 ) ? 1 : 2;
91  int idQuark = ( rndmPtr->flat() < 0.5 ) ? 1 : 3;
92  int newCol = process.nextColTag();
93 
94  int countW(0);
95  for (int i = 0; i < process.size(); ++i) {
96  if (process[i].idAbs() == 24 && process[i].status() == -22) {
97  countW++;
98  if( countW > 2 ) break;
99  int idFermion = ( hadronW == countW ) ? idQuark : 11 ;
100  // Identify W+- daughters, properly ordered.
101  int idV = process[i].id();
102  int i1 = process[i].daughter1();
103  int i2 = process[i].daughter2();
104  int id1 = process[i1].idAbs();
105  int id2 = process[i2].idAbs();
106  if( id1 == idFermion || id2 == idFermion ) continue;
107 
108  double mV = process[i].m();
109  Vec4 p1 = process[i1].p();
110  p1.bstback( process[i].p() );
111  double pMag = p1.pAbs();
112 
113  // The current distributions are for the particle with the
114  // same charge sign as the mother, i.e. W- -> e-.
115  if (process[i1].id() * idV > 0) swap( i1, i2);
116  process[i1].id(idFermion * process[i1].id()/abs(process[i1].id()) );
117  process[i2].id((idFermion+1) * process[i2].id()/abs(process[i2].id()) );
118  if( idFermion != 11 ) {
119  if( process[i1].id() > 0 ) {
120  process[i1].col(newCol);
121  process[i2].acol(newCol);
122 
123  } else {
124  process[i2].col(newCol);
125  process[i1].acol(newCol);
126  }
127  }
128 
129  double m1 = m0powheg(process[i1].idAbs());
130  double m2 = m0powheg(process[i2].idAbs());
131 
132  // Energy and absolute momentum of first decay product in W rest frame.
133  double e1 = 0.5* (pow2(mV) + pow2(m1) - pow2(m2))/mV;
134  double pA = sqrt(pow2(e1) - pow2(m1));
135 
136  p1.rescale3( pA/ pMag );
137  p1.e( e1 );
138  Vec4 p2 = Vec4(0,0,0,mV) - p1;
139 
140 
141  p1.bst( process[i].p() );
142  p2.bst( process[i].p() );
143  process[i1].p( p1 );
144  process[i2].p( p2 );
145  process[i1].m( m1 );
146  process[i2].m( m2 );
147  }
148  // End of loop over W's. Do not veto any events.
149 
150  }
151  return false;
152 }
153 
154 //==========================================================================
155 
156 } // end namespace Pythia8
157 
158 #endif // end Pythia8_SetLHEDecayProductHooks_H
constexpr double pow2(const double &x)
Powers of small integers - for balance speed/code clarity.
Definition: PythiaStdlib.h:182
bool flag(string keyIn)
Give back current value, with check that key exists.
Definition: Settings.cc:1478
Settings * settingsPtr
Pointer to the settings database.
Definition: PhysicsBase.h:81
void bstback(const Vec4 &pIn)
Boost given by a Vec4 p; boost in opposite direction.
Definition: Basics.cc:509
void rescale3(double fac)
Member functions that perform operations.
Definition: Basics.h:95
Definition: SetLHEDecayProductHook.h:24
The Event class holds all info on the generated event.
Definition: Event.h:453
bool canVetoProcessLevel() override
Override base class methods.
Definition: SetLHEDecayProductHook.h:31
bool initAfterBeams() override
Intialize the user hook after the beams.
Definition: SetLHEDecayProductHook.h:76
Rndm * rndmPtr
Pointer to the random number generator.
Definition: PhysicsBase.h:93
SetLHEDecayProductHook(Settings &settings, const ParticleData *pdtPtrIn)
Constructor.
Definition: SetLHEDecayProductHook.h:66
void bst(double betaX, double betaY, double betaZ)
Boost (simple).
Definition: Basics.cc:434
int size() const
Event record size.
Definition: Event.h:504
double m2(const Vec4 &v1)
The squared invariant mass of one or more four-vectors.
Definition: Basics.cc:605
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
bool checkVetoProcessLevel(Event &process)
Class specific.
Definition: SetLHEDecayProductHook.h:85
UserHooks is base class for user access to program execution.
Definition: UserHooks.h:32
bool doVetoProcessLevel(Event &process) override
Definition: SetLHEDecayProductHook.h:32
double flat()
Generate next random number uniformly between 0 and 1.
Definition: Basics.cc:189
This class holds a map of all ParticleDataEntries.
Definition: ParticleData.h:422
Definition: Basics.h:32
Definition: Settings.h:195
void addFlag(string keyIn, bool defaultIn)
Add new entry.
Definition: Settings.h:267