GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PythiaBaseHadro2019.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
8 
9  Changes required to implement the GENIE Boosted Dark Matter module
10  were installed by Josh Berger (Univ. of Wisconsin)
11 */
12 //____________________________________________________________________________
13 
27 
28 using namespace genie;
29 using namespace genie::constants;
30 
31 //____________________________________________________________________________
34 {
35 
36 }
37 //____________________________________________________________________________
40 {
41 
42 }
43 //____________________________________________________________________________
44 PythiaBaseHadro2019::PythiaBaseHadro2019(string name, string config) :
45 EventRecordVisitorI(name, config)
46 {
47 
48 }
49 //____________________________________________________________________________
51 {
52 
53 }
54 //____________________________________________________________________________
56 {
57  Interaction * interaction = event->Summary();
58 
59  bool valid_process = this->AssertValidity(interaction);
60  if(!valid_process) {
61  LOG("PythiaHad", pFATAL)
62  << "Input interaction type is not allowed!!!";
63  LOG("PythiaHad", pFATAL)
64  << *event;
65  gAbortingInErr = true;
66  std::exit(1);
67  }
68 
69  // Decide the leading quark and remnant diquark PDG codes for this event
70  this->MakeQuarkDiquarkAssignments(interaction);
71 
72  // Copy original and set required PYTHIA decay flags
73  this->CopyOriginalDecayFlags();
74  this->SetDesiredDecayFlags();
75 
76  // Call PYTHIA6 or PYTHIA8 to obtain the fragmentation products
77  //TClonesArray * particle_list = this->Hadronize(interaction);
78  bool hadronized = this->Hadronize(event);
79 
80  // Restore decay flags
82 
83  if(!hadronized) {
84  LOG("PythiaHad", pWARN) << "Hadronization failed!";
85  event->EventFlags()->SetBitNumber(kHadroSysGenErr, true);
87  exception.SetReason("Could not simulate the hadronic system");
88  exception.SwitchOnFastForward();
89  throw exception;
90  return;
91  }
92 }
93 //____________________________________________________________________________
95  const Interaction * interaction) const
96 {
97  LOG("PythiaHad", pNOTICE)
98  << "Making leading quark / remnant di-quark assignments";
99 
100  // get kinematics / init-state / process-info
101 
102  const Kinematics & kinematics = interaction->Kine();
103  const InitialState & init_state = interaction->InitState();
104  const ProcessInfo & proc_info = interaction->ProcInfo();
105  const Target & target = init_state.Tgt();
106 
107  assert(target.HitQrkIsSet());
108 
109  double W = kinematics.W();
110 
111  int probe = init_state.ProbePdg();
112  int hit_nucleon = target.HitNucPdg();
113  int hit_quark = target.HitQrkPdg();
114  bool from_sea = target.HitSeaQrk();
115 
116  LOG("PythiaHad", pNOTICE)
117  << "Hit nucleon pdgc = " << hit_nucleon << ", W = " << W;
118  LOG("PythiaHad", pNOTICE)
119  << "Selected hit quark pdgc = " << hit_quark
120  << ((from_sea) ? "[sea]" : "[valence]");
121 
122  // check hit-nucleon assignment, input neutrino & interaction type
123  bool isp = pdg::IsProton (hit_nucleon);
124  bool isn = pdg::IsNeutron (hit_nucleon);
125  bool isv = pdg::IsNeutrino (probe);
126  bool isvb = pdg::IsAntiNeutrino (probe);
127 //bool isl = pdg::IsNegChargedLepton (probe);
128 //bool islb = pdg::IsPosChargedLepton (probe);
129  bool iscc = proc_info.IsWeakCC ();
130  bool isnc = proc_info.IsWeakNC ();
131  bool isdm = proc_info.IsDarkMatter ();
132  bool isem = proc_info.IsEM ();
133  bool isu = pdg::IsUQuark (hit_quark);
134  bool isd = pdg::IsDQuark (hit_quark);
135  bool iss = pdg::IsSQuark (hit_quark);
136  bool isub = pdg::IsAntiUQuark (hit_quark);
137  bool isdb = pdg::IsAntiDQuark (hit_quark);
138  bool issb = pdg::IsAntiSQuark (hit_quark);
139 
140  //
141  // Generate the quark system (q + qq) initiating the hadronization
142  //
143 
144  int leading_quark = 0; // leading quark (hit quark after the interaction)
145  int remnant_diquark = 0; // remnant diquark (xF<0 at hadronic CMS)
146 
147  // Figure out the what happens to the hit quark after the interaction
148  if (isnc || isem || isdm) {
149  // NC, EM
150  leading_quark = hit_quark;
151  } else {
152  // CC
153  if (isv && isd ) leading_quark = kPdgUQuark;
154  else if (isv && iss ) leading_quark = kPdgUQuark;
155  else if (isv && isub) leading_quark = kPdgAntiDQuark;
156  else if (isvb && isu ) leading_quark = kPdgDQuark;
157  else if (isvb && isdb) leading_quark = kPdgAntiUQuark;
158  else if (isvb && issb) leading_quark = kPdgAntiUQuark;
159  else {
160  LOG("PythiaHad", pERROR)
161  << "Not allowed mode. Refused to make a leading quark assignment!";
162  return;
163  }
164  }//CC
165 
166  // Figure out what the remnant diquark is.
167  // Note from Hugh, following a conversation with his local HEP theorist
168  // (Gary Goldstein): "I am told that the probability of finding the diquark
169  // in the singlet vs. triplet states is 50-50."
170 
171  // hit quark = valence quark
172  if(!from_sea) {
173  if (isp && isu) remnant_diquark = kPdgUDDiquarkS1; /* u(->q) + ud */
174  if (isp && isd) remnant_diquark = kPdgUUDiquarkS1; /* d(->q) + uu */
175  if (isn && isu) remnant_diquark = kPdgDDDiquarkS1; /* u(->q) + dd */
176  if (isn && isd) remnant_diquark = kPdgUDDiquarkS1; /* d(->q) + ud */
177  }
178  // hit quark = sea quark
179  else {
180  if(isp && isu) remnant_diquark = kPdgUDDiquarkS1; /* u(->q) + bar{u} uud (=ud) */
181  if(isp && isd) remnant_diquark = kPdgUUDiquarkS1; /* d(->q) + bar{d} uud (=uu) */
182  if(isn && isu) remnant_diquark = kPdgDDDiquarkS1; /* u(->q) + bar{u} udd (=dd) */
183  if(isn && isd) remnant_diquark = kPdgUDDiquarkS1; /* d(->q) + bar{d} udd (=ud) */
184 
185  // The following section needs revisiting.
186 
187  // The lepton is scattered off a sea antiquark, materializing its quark
188  // partner and leaving me with a 5q system ( <qbar + q> + qqq(valence) )
189  // I will force few qbar+q annhilations below to get my quark/diquark system
190  // Probably it is best to leave the qqq system in the final state and then
191  // just do the fragmentation of the qbar q system? But how do I figure out
192  // how to split the available energy?
193 
194  /* bar{u} (-> bar{d}) + u uud => u + uu */
195  if(isp && isub && iscc) {
196  leading_quark = kPdgUQuark;
197  remnant_diquark = kPdgUUDiquarkS1;
198  }
199  /* bar{u} (-> bar{u}) + u uud => u + ud */
200  if(isp && isub && (isnc||isem||isdm)) {
201  leading_quark = kPdgUQuark;
202  remnant_diquark = kPdgUDDiquarkS1;
203  }
204  /* bar{d} (-> bar{u}) + d uud => d + ud */
205  if(isp && isdb && iscc) {
206  leading_quark = kPdgDQuark;
207  remnant_diquark = kPdgUDDiquarkS1;
208  }
209  /* bar{d} (-> bar{d}) + d uud => d + uu */
210  if(isp && isdb && (isnc||isem||isdm)) {
211  leading_quark = kPdgDQuark;
212  remnant_diquark = kPdgUUDiquarkS1;
213  }
214  /* bar{u} (-> bar{d}) + u udd => u + ud */
215  if(isn && isub && iscc) {
216  leading_quark = kPdgUQuark;
217  remnant_diquark = kPdgUDDiquarkS1;
218  }
219  /* bar{u} (-> bar{u}) + u udd => u + dd */
220  if(isn && isub && (isnc||isem||isdm)) {
221  leading_quark = kPdgUQuark;
222  remnant_diquark = kPdgDDDiquarkS1;
223  }
224  /* bar{d} (-> bar{u}) + d udd => d + dd */
225  if(isn && isdb && iscc) {
226  leading_quark = kPdgDQuark;
227  remnant_diquark = kPdgDDDiquarkS1;
228  }
229  /* bar{d} (-> bar{d}) + d udd => d + ud */
230  if(isn && isdb && (isnc||isem||isdm)) {
231  leading_quark = kPdgDQuark;
232  remnant_diquark = kPdgUDDiquarkS1;
233  }
234 
235  // The neutrino is scatterred off s or sbar sea quarks
236  // For the time being I will handle s like d and sbar like dbar (copy & paste
237  // from above) so that I conserve charge.
238 
239  if(iss || issb) {
240  LOG("PythiaHad", pNOTICE)
241  << "Can not really handle a hit s or sbar quark / Faking it";
242 
243  if(isp && iss) { remnant_diquark = kPdgUUDiquarkS1; }
244  if(isn && iss) { remnant_diquark = kPdgUDDiquarkS1; }
245 
246  if(isp && issb && iscc) {
247  leading_quark = kPdgDQuark;
248  remnant_diquark = kPdgUDDiquarkS1;
249  }
250  if(isp && issb && (isnc||isem||isdm)) {
251  leading_quark = kPdgDQuark;
252  remnant_diquark = kPdgUUDiquarkS1;
253  }
254  if(isn && issb && iscc) {
255  leading_quark = kPdgDQuark;
256  remnant_diquark = kPdgDDDiquarkS1;
257  }
258  if(isn && issb && (isnc||isem||isdm)) {
259  leading_quark = kPdgDQuark;
260  remnant_diquark = kPdgUDDiquarkS1;
261  }
262  }
263 
264  // if the diquark is a ud, switch it to the singlet state with 50% probability
265  if(remnant_diquark == kPdgUDDiquarkS1) {
266  RandomGen * rnd = RandomGen::Instance();
267  double Rqq = rnd->RndHadro().Rndm();
268  if(Rqq<0.5) remnant_diquark = kPdgUDDiquarkS0;
269  }
270  }
271 
272  fLeadingQuark = leading_quark;
273  fRemnantDiquark = remnant_diquark;
274 }
275 //____________________________________________________________________________
276 bool PythiaBaseHadro2019::AssertValidity(const Interaction * interaction) const {
277 
278  // check that there is no charm production
279  // (GENIE uses a special model for these cases)
280  if(interaction->ExclTag().IsCharmEvent()) {
281  LOG("PythiaHad", pWARN) << "Can't hadronize charm events";
282  return false;
283  }
284  // check the available mass
285  double W = utils::kinematics::W(interaction);
286  double Wmin = kNucleonMass+kPionMass;
287  if(W < Wmin) {
288  LOG("PythiaHad", pWARN) << "Low invariant mass, W = "
289  << W << " GeV!!";
290  return false;
291  }
292 
293  const InitialState & init_state = interaction->InitState();
294  const ProcessInfo & proc_info = interaction->ProcInfo();
295  const Target & target = init_state.Tgt();
296 
297  if( ! target.HitQrkIsSet() ) {
298  LOG("PythiaHad", pWARN) << "Hit quark was not set!";
299  return false;
300  }
301 
302  int probe = init_state.ProbePdg();
303  int hit_nucleon = target.HitNucPdg();
304  int hit_quark = target.HitQrkPdg();
305  //bool from_sea = target.HitSeaQrk();
306 
307  // check hit-nucleon assignment, input neutrino & weak current
308  bool isp = pdg::IsProton (hit_nucleon);
309  bool isn = pdg::IsNeutron (hit_nucleon);
310  bool isv = pdg::IsNeutrino (probe);
311  bool isvb = pdg::IsAntiNeutrino (probe);
312  bool isdm = pdg::IsDarkMatter (probe);
313  bool isl = pdg::IsNegChargedLepton (probe);
314  bool islb = pdg::IsPosChargedLepton (probe);
315  bool iscc = proc_info.IsWeakCC ();
316  bool isnc = proc_info.IsWeakNC ();
317  bool isdmi = proc_info.IsDarkMatter ();
318  bool isem = proc_info.IsEM ();
319  if( !(iscc||isnc||isem||isdmi) ) {
320  LOG("PythiaHad", pWARN)
321  << "Can only handle electro-weak interactions";
322  return false;
323  }
324  if( !(isp||isn) || !(isv||isvb||isl||islb||isdm) ) {
325  LOG("PythiaHad", pWARN)
326  << "Invalid initial state: probe = "
327  << probe << ", hit_nucleon = " << hit_nucleon;
328  return false;
329  }
330 
331  // assert that the interaction mode is allowed
332  bool isu = pdg::IsUQuark (hit_quark);
333  bool isd = pdg::IsDQuark (hit_quark);
334  bool iss = pdg::IsSQuark (hit_quark);
335  bool isub = pdg::IsAntiUQuark (hit_quark);
336  bool isdb = pdg::IsAntiDQuark (hit_quark);
337  bool issb = pdg::IsAntiSQuark (hit_quark);
338 
339  bool allowed = (iscc && isv && (isd||isub||iss)) ||
340  (iscc && isvb && (isu||isdb||issb)) ||
341  (isnc && (isv||isvb) && (isu||isd||isub||isdb||iss||issb)) ||
342  (isdmi && isdm && (isu||isd||isub||isdb||iss||issb)) ||
343  (isem && (isl||islb) && (isu||isd||isub||isdb||iss||issb));
344  if(!allowed) {
345  LOG("PythiaHad", pWARN)
346  << "Impossible interaction type / probe / hit quark combination!";
347  return false;
348  }
349 
350  return true;
351 }
352 //____________________________________________________________________________
354 {
355  fLeadingQuark = 0;
356  fRemnantDiquark = 0;
357  fSSBarSuppression = 0.;
358  fGaussianPt2 = 0.;
359  fNonGaussianPt2Tail = 0.;
360  fRemainingECutoff = 0.;
361  fDiQuarkSuppression = 0.;
363  fSVMesonSuppression = 0.;
364  fLunda = 0.;
365  fLundb = 0.;
366  fLundaDiq = 0.;
367  fOriDecayFlag_pi0 = false;
368  fOriDecayFlag_K0 = false;
369  fOriDecayFlag_K0b = false;
370  fOriDecayFlag_L0 = false;
371  fOriDecayFlag_L0b = false;
372  fOriDecayFlag_Dm = false;
373  fOriDecayFlag_D0 = false;
374  fOriDecayFlag_Dp = false;
375  fOriDecayFlag_Dpp = false;
376  fReqDecayFlag_pi0 = false;
377  fReqDecayFlag_K0 = false;
378  fReqDecayFlag_K0b = false;
379  fReqDecayFlag_L0 = false;
380  fReqDecayFlag_L0b = false;
381  fReqDecayFlag_Dm = false;
382  fReqDecayFlag_D0 = false;
383  fReqDecayFlag_Dp = false;
384  fReqDecayFlag_Dpp = false;
385 }
386 //____________________________________________________________________________
388 {
389  // Get PYTHIA physics configuration parameters specified by GENIE
390  this->GetParam( "PYTHIA-SSBarSuppression", fSSBarSuppression );
391  this->GetParam( "PYTHIA-GaussianPt2", fGaussianPt2 );
392  this->GetParam( "PYTHIA-NonGaussianPt2Tail", fNonGaussianPt2Tail );
393  this->GetParam( "PYTHIA-RemainingEnergyCutoff", fRemainingECutoff );
394  this->GetParam( "PYTHIA-DiQuarkSuppression", fDiQuarkSuppression );
395  this->GetParam( "PYTHIA-LightVMesonSuppression", fLightVMesonSuppression );
396  this->GetParam( "PYTHIA-SVMesonSuppression", fSVMesonSuppression );
397  this->GetParam( "PYTHIA-Lunda", fLunda );
398  this->GetParam( "PYTHIA-Lundb", fLundb );
399  this->GetParam( "PYTHIA-LundaDiq", fLundaDiq );
400 
401  // Set required PYTHIA decay flags
402  fReqDecayFlag_pi0 = false; // don't decay pi0
403  fReqDecayFlag_K0 = false; // don't decay K0
404  fReqDecayFlag_K0b = false; // don't decay \bar{K0}
405  fReqDecayFlag_L0 = false; // don't decay Lambda0
406  fReqDecayFlag_L0b = false; // don't decay \bar{Lambda0}
407  fReqDecayFlag_Dm = true; // decay Delta-
408  fReqDecayFlag_D0 = true; // decay Delta0
409  fReqDecayFlag_Dp = true; // decay Delta+
410  fReqDecayFlag_Dpp = true; // decay Delta++
411 
412  // Load Wcut determining the phase space area where the multiplicity prob.
413  // scaling factors would be applied -if requested-
414  double Wcut, Wmin ;
415  this->GetParam( "Wcut", Wcut );
416  this->GetParam( "KNO2PYTHIA-Wmin", Wmin );
417 
418  if ( Wcut > Wmin ) {
419  LOG("PythiaHad", pERROR)
420  << "Wcut value too high and in conflict with the KNO2PYTHIA-Wmin!"
421  << "\n Wcut = " << Wcut
422  << "\n KNO2PYTHIA-Wmin = " << Wmin;
423  }
424 }
425 //____________________________________________________________________________
double fSSBarSuppression
ssbar suppression
const int kPdgUUDiquarkS1
Definition: PDGCodes.h:58
double W(bool selected=false) const
Definition: Kinematics.cxx:157
bool HitSeaQrk(void) const
Definition: Target.cxx:299
bool IsWeakCC(void) const
virtual void MakeQuarkDiquarkAssignments(const Interaction *in) const
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
#define pERROR
Definition: Messenger.h:59
static const double kNucleonMass
bool IsUQuark(int pdgc)
Definition: PDGUtils.cxx:266
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
int HitNucPdg(void) const
Definition: Target.cxx:304
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
virtual void SetDesiredDecayFlags(void) const =0
int HitQrkPdg(void) const
Definition: Target.cxx:242
virtual void RestoreOriginalDecayFlags(void) const =0
#define pFATAL
Definition: Messenger.h:56
const int kPdgUQuark
Definition: PDGCodes.h:42
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
virtual void ProcessEventRecord(GHepRecord *event) const
double fGaussianPt2
gaussian pt2 distribution width
bool IsDarkMatter(int pdgc)
Definition: PDGUtils.cxx:127
bool IsSQuark(int pdgc)
Definition: PDGUtils.cxx:276
double fDiQuarkSuppression
di-quark suppression parameter
bool IsAntiSQuark(int pdgc)
Definition: PDGUtils.cxx:306
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
bool IsAntiDQuark(int pdgc)
Definition: PDGUtils.cxx:301
double fSVMesonSuppression
strange vector meson suppression
bool IsCharmEvent(void) const
Definition: XclsTag.h:50
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
bool IsPosChargedLepton(int pdgc)
Definition: PDGUtils.cxx:148
Summary information for an interaction.
Definition: Interaction.h:56
const int kPdgAntiUQuark
Definition: PDGCodes.h:43
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
bool IsWeakNC(void) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
double fLightVMesonSuppression
light vector meson suppression
const int kPdgUDDiquarkS1
Definition: PDGCodes.h:57
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
double fLundaDiq
adjustment of Lund a for di-quark
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:118
const int kPdgAntiDQuark
Definition: PDGCodes.h:45
const Kinematics & Kine(void) const
Definition: Interaction.h:71
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
int ProbePdg(void) const
Definition: InitialState.h:64
double fNonGaussianPt2Tail
non gaussian pt2 tail parameterization
double fLunda
Lund a parameter.
virtual void CopyOriginalDecayFlags(void) const =0
const int kPdgDQuark
Definition: PDGCodes.h:44
const int kPdgUDDiquarkS0
Definition: PDGCodes.h:56
virtual bool Hadronize(GHepRecord *event) const =0
#define pWARN
Definition: Messenger.h:60
bool IsEM(void) const
TRandom3 & RndHadro(void) const
rnd number generator used by hadronization models
Definition: RandomGen.h:53
bool HitQrkIsSet(void) const
Definition: Target.cxx:292
bool IsDarkMatter(void) const
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
bool IsDQuark(int pdgc)
Definition: PDGUtils.cxx:271
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double fLundb
Lund b parameter.
#define pNOTICE
Definition: Messenger.h:61
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:66
virtual bool AssertValidity(const Interaction *in) const
bool gAbortingInErr
Definition: Messenger.cxx:34
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
double fRemainingECutoff
remaining E cutoff stopping fragmentation
bool IsNegChargedLepton(int pdgc)
Definition: PDGUtils.cxx:139
bool IsAntiUQuark(int pdgc)
Definition: PDGUtils.cxx:296
Initial State information.
Definition: InitialState.h:48
const int kPdgDDDiquarkS1
Definition: PDGCodes.h:55