GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
LeptoHadronization.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  Alfonso Garcia <alfonsog \at nikhef.nl>
7  NIKHEF (Amsterdam)
8 */
9 //____________________________________________________________________________
10 
11 #include <RVersion.h>
12 #include <TClonesArray.h>
13 #include <TMath.h>
14 
17 #include "Framework/Conventions/GBuild.h"
29 
30 #ifdef __GENIE_PYTHIA6_ENABLED__
31 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,15,6)
32 #include <TMCParticle.h>
33 #else
34 #include <TMCParticle6.h>
35 #endif
36 #endif // __GENIE_PYTHIA6_ENABLED__
37 
38 using namespace genie;
39 using namespace genie::constants;
40 
41 // the actual PYTHIA call
42 extern "C" {
43  double pyangl_( double *, double * );
44  void pykfdi_( int *, int *, int *, int * );
45  void pyzdis_( int *, int *, double *, double * );
46  void pyrobo_( int *, int *, double *, double *, double *, double *, double * );
47  void pydecy_( int * );
48  void py2ent_( int *, int *, int *, double * );
49 }
50 
51 //____________________________________________________________________________
53 EventRecordVisitorI("genie::LeptoHadronization")
54 {
55  this->Initialize();
56 }
57 //____________________________________________________________________________
59 EventRecordVisitorI("genie::LeptoHadronization", config)
60 {
61  this->Initialize();
62 }
63 //____________________________________________________________________________
65 {
66 
67 }
68 //____________________________________________________________________________
70 {
71 #ifdef __GENIE_PYTHIA6_ENABLED__
72  fPythia = TPythia6::Instance();
73 
74  // sync GENIE/PYTHIA6 seed number
76 #endif
77 }
78 //____________________________________________________________________________
81  event // avoid unused variable warning if PYTHIA6 is not enabled
82  #endif
83 ) const
84 {
85 
86 #ifdef __GENIE_PYTHIA6_ENABLED__
87 
88  if(!this->Hadronize(event)) {
89  LOG("LeptoHad", pWARN) << "Hadronization failed!";
90  event->EventFlags()->SetBitNumber(kHadroSysGenErr, true);
92  exception.SetReason("Could not simulate the hadronic system");
93  exception.SwitchOnFastForward();
94  throw exception;
95  return;
96  }
97 
98 #else
99  LOG("LeptoHad", pFATAL)
100  << "Calling GENIE/PYTHIA6 without enabling PYTHIA6";
101  gAbortingInErr = true;
102  std::exit(1);
103 #endif
104 
105 }
106 //____________________________________________________________________________
109  event // avoid unused variable warning if PYTHIA6 is not enabled
110 #endif
111 ) const
112 {
113 
114 #ifdef __GENIE_PYTHIA6_ENABLED__
115  // Compute kinematics of hadronic system with energy/momentum conservation
116  LongLorentzVector p4v( * event->Probe()->P4() );
117  LongLorentzVector p4N( * event->HitNucleon()->P4() );
118  LongLorentzVector p4l( * event->FinalStatePrimaryLepton()->P4() );
119  LongLorentzVector p4Hadlong( p4v.Px()+p4N.Px()-p4l.Px(), p4v.Py()+p4N.Py()-p4l.Py(), p4v.Pz()+p4N.Pz()-p4l.Pz(), p4v.E()+p4N.E()-p4l.E() );
120 
121  LOG("LeptoHad", pDEBUG) << "v [LAB']: " << p4v.E() << " // " << p4v.M2() << " // [ " << p4v.Dx() << " , " << p4v.Dy() << " , " << p4v.Dz() << " ]";
122  LOG("LeptoHad", pDEBUG) << "N [LAB']: " << p4N.E() << " // " << p4N.M2() << " // [ " << p4N.Dx() << " , " << p4N.Dy() << " , " << p4N.Dz() << " ]";
123  LOG("LeptoHad", pDEBUG) << "l [LAB']: " << p4l.E() << " // " << p4l.M2() << " // [ " << p4l.Dx() << " , " << p4l.Dy() << " , " << p4l.Dz() << " ]";
124  LOG("LeptoHad", pDEBUG) << "H [LAB']: " << p4Hadlong.E() << " // " << p4Hadlong.M2() << " // [ " << p4Hadlong.Dx() << " , " << p4Hadlong.Dy() << " , " << p4Hadlong.Dz() << " ]";
125 
126  // Translate from long double to double
127  const TLorentzVector & vtx = *( event->Probe()->X4());
128  TLorentzVector p4Had( (double)p4Hadlong.Px(), (double)p4Hadlong.Py(), (double)p4Hadlong.Pz(), (double)p4Hadlong.E() );
129  event->AddParticle(kPdgHadronicSyst, kIStDISPreFragmHadronicState, event->HitNucleonPosition(),-1,-1,-1, p4Had, vtx);
130 
131  const Interaction * interaction = event->Summary();
132  interaction->KinePtr()->SetHadSystP4(p4Had);
133  interaction->KinePtr()->SetW(p4Hadlong.M());
134 
135  double W = interaction->Kine().W();
136  if(W < fWmin) {
137  LOG("LeptoHad", pWARN) << "Low invariant mass, W = " << W << " GeV!!";
138  return 0;
139  }
140 
141  const XclsTag & xclstag = interaction->ExclTag();
142  const Target & target = interaction->InitState().Tgt();
143 
144  assert(target.HitQrkIsSet());
145 
146  bool isp = pdg::IsProton(target.HitNucPdg());
147  int hit_quark = target.HitQrkPdg();
148  int frag_quark = xclstag.FinalQuarkPdg();
149 
150  LOG("LeptoHad", pDEBUG) << "Hit nucleon pdgc = " << target.HitNucPdg() << ", W = " << W;
151  LOG("LeptoHad", pDEBUG) << "Selected hit quark pdgc = " << hit_quark << " // Fragmentation quark = " << frag_quark;
152 
153  RandomGen * rnd = RandomGen::Instance();
154 
155  //
156  // Generate the hadron combination to input PYTHIA
157  //
158 
159  double pmas1_W = fPythia->GetPMAS(24,1);
160  double pmas2_W = fPythia->GetPMAS(24,2);
161  double pmas2_t = fPythia->GetPMAS(6,2);
162  fPythia->SetPMAS(24,1,kMw); //mass of the W boson (pythia=80.450 // genie=80.385)
163  fPythia->SetPMAS(24,2,0.); //set to 0 the width of the W boson to avoid problems with energy conservation
164  fPythia->SetPMAS(6,2,0.); //set to 0 the width of the top to avoid problems with energy conservation
165 
166  //If the hit quark is a d we have these options:
167  /* uud(->q) => uu + q */
168  /* uud d(->q)db => uu + q (d valence and db sea annihilates)*/
169  /* udd(->q) => ud + q */
170  /* udd d(->q)db => ud + q (d valence and db sea annihilates)*/
171  if ( pdg::IsDQuark(hit_quark) ) {
172  // choose diquark system depending on proton or neutron
173  int diquark = 0;
174  if (isp) diquark = kPdgUUDiquarkS1;
175  else diquark = rnd->RndHadro().Rndm()>0.75 ? kPdgUDDiquarkS1 : kPdgUDDiquarkS0;
176  // Check that the trasnferred energy is higher than the mass of the produced quarks
177  double m_frag = PDGLibrary::Instance()->Find(frag_quark)->Mass();
178  double m_diquark = PDGLibrary::Instance()->Find(diquark)->Mass();
179  if( W <= m_frag + m_diquark + fMinESinglet ) {
180  LOG("LeptoHad", pWARN) << "Low invariant mass, W = " << W << " GeV! Returning a null list";
181  LOG("LeptoHad", pWARN) << "frag_quark = " << frag_quark << " -> m = " << m_frag;
182  LOG("LeptoHad", pWARN) << "diquark = " << diquark << " -> m = " << m_diquark;
183  return 0;
184  }
185  // Input the two particles to PYTHIA back to back in the CM frame
186  // If a top quark is produced we decay it because it does not hadronize
187  fPythia->Py1ent( -1, frag_quark, (W*W - m_diquark*m_diquark + m_frag*m_frag)/2./W, 0., 0. ); //k(1,2) = 2
188  if ( pdg::IsTQuark(frag_quark) ) {
189  int ip = 1;
190  pydecy_(&ip);
191  }
192  fPythia->Py1ent( fPythia->GetN()+1, diquark, (W*W + m_diquark*m_diquark - m_frag*m_frag)/2./W, fPythia->GetPARU(1), 0. ); //k(2,2) = 1
193  }
194 
195  //If the hit quark is a u we have these options:
196  /* u(->q)ud => ud + q */
197  /* uud u(->q)ub => ud + q (u valence and ub sea annihilates)*/
198  /* u(->q)dd => dd + q */
199  /* udd u(->q)ub => dd + q (u valence and ub sea annihilates)*/
200  else if ( pdg::IsUQuark(hit_quark) ) {
201  // choose diquark system depending on proton or neutron
202  int diquark = 0;
203  if (isp) diquark = rnd->RndHadro().Rndm()>0.75 ? kPdgUDDiquarkS1 : kPdgUDDiquarkS0;
204  else diquark = kPdgDDDiquarkS1;
205  // Check that the trasnferred energy is higher than the mass of the produced quarks.
206  double m_frag = PDGLibrary::Instance()->Find(frag_quark)->Mass();
207  double m_diquark = PDGLibrary::Instance()->Find(diquark)->Mass();
208  if( W <= m_frag + m_diquark + fMinESinglet ) {
209  LOG("LeptoHad", pWARN) << "Low invariant mass, W = " << W << " GeV! Returning a null list";
210  LOG("LeptoHad", pWARN) << "frag_quark = " << frag_quark << " -> m = " << m_frag;
211  LOG("LeptoHad", pWARN) << "diquark = " << diquark << " -> m = " << m_diquark;
212  return 0;
213  }
214  // Input the two particles to PYTHIA back to back in the CM frame
215  fPythia->Py1ent( -1, frag_quark, (W*W - m_diquark*m_diquark + m_frag*m_frag)/2./W, 0., 0. ); //k(1,2) = 2
216  fPythia->Py1ent( fPythia->GetN()+1, diquark, (W*W + m_diquark*m_diquark - m_frag*m_frag)/2./W, fPythia->GetPARU(1), 0. ); //k(2,2) = 1
217  }
218  else {
219 
220  // If the hit quark is not u or d then is more complicated.
221  // We are using the same procedure use in LEPTO (see lqev.F)
222  // Our initial systemt will look like this -> qqq + hit_q(->frag_q) + rema_q
223  // And we have to input PYTHIA something like this -> frag_q + rema + hadron
224  // These are the posible combinations -> frag_q[q] + meson [qqb] + diquark [qq]
225  // -> frag_q[qb] + baryon [qqq] + quark [q]
226 
227  // Remnant of the hit quark (which is from the sea) will be of opposite charge
228  int rema_hit_quark = -hit_quark;
229 
230  // Check that the trasnfered energy is higher than the mass of the produce quarks plus remnant quark and nucleon
231  double m_frag = PDGLibrary::Instance()->Find(frag_quark)->Mass();
232  double m_rema_hit = PDGLibrary::Instance()->Find(rema_hit_quark)->Mass();
233  if (W <= m_frag + m_rema_hit + 0.9 + fMinESinglet ) {
234  LOG("LeptoHad", pWARN) << "Low invariant mass, W = " << W << " GeV! Returning a null list";
235  LOG("LeptoHad", pWARN) << " frag_quark = " << frag_quark << " -> m = " << m_frag;
236  LOG("LeptoHad", pWARN) << " rema_hit_quark = " << rema_hit_quark << " -> m = " << m_rema_hit;
237  return 0;
238  }
239 
240  //PDG of the two hadronic particles for the final state
241  int hadron = 0;
242  int rema = 0;
243 
244  int ntwoq = isp ? 2 : 1; //proton two ups & neutron one up
245  int counter = 0;
246 
247  // Here we select the id and kinematics of the hadron and rema particles
248  // Some combinations can be kinematically forbiden so we repeat this process
249  // up to 100 times before the event is discarded.
250  while( counter<fMaxIterHad ) {
251 
252  // Loop to create a combination of hadron + rema. Two options are possible:
253  // 1) diquark [qq] + meson [qqb]
254  // 2) quark [q] + baryon [qqq]
255  while(hadron==0) {
256  //choose a valence quark and the remaining will be a diquark system
257  int valquark = int(1.+ntwoq/3.+rnd->RndHadro().Rndm());
258  int diquark = 0;
259  if ( valquark==ntwoq ) diquark = rnd->RndHadro().Rndm()>0.75 ? kPdgUDDiquarkS1 : kPdgUDDiquarkS0;
260  else diquark = 1000*ntwoq+100*ntwoq+3;
261 
262  // Choose flavours using PYTHIA tool
263  int idum;
264  if ( rema_hit_quark>0 ) { //create a baryon (qqq)
265  pykfdi_(&diquark,&rema_hit_quark,&idum,&hadron);
266  rema = valquark;
267  }
268  else { //create a meson (qqbar)
269  pykfdi_(&valquark,&rema_hit_quark,&idum,&hadron);
270  rema = diquark;
271  }
272  }
273 
274  double m_hadron = PDGLibrary::Instance()->Find(hadron)->Mass();
275  double m_rema = PDGLibrary::Instance()->Find(rema)->Mass();
276 
277  // Give balancing pT to hadron and rema particles
278  double pT = fRemnantPT * TMath::Sqrt( -1*TMath::Log( rnd->RndHadro().Rndm() ) );
279  double pT2 = TMath::Power(pT,2);
280  double pr = TMath::Power(m_hadron,2)+pT2;
281  int kfl1 = 1;
282  int kfl3 = 0;
283  double z;
284  //to generate the longitudinal scaling variable z in jet fragmentation using PYTHIA function
285  // Split energy-momentum of remnant using PYTHIA function
286  // z=E-pz fraction for rema forming jet-system with frag_q
287  // 1-z=E-pz fraction for hadron
288  pyzdis_(&kfl1,&kfl3,&pr,&z);
289 
290  // Energy of trasnfered to the hadron
291  double tm_hadron = pr / z / W;
292  double E_hadron = 0.5 * ( z*W + tm_hadron ); //E_hadron - pz = zW
293  double E_pz = W - tm_hadron;
294  double WT = (1-z) * W * E_pz - pT2;
295 
296  // Check if energy in jet system is enough for fragmentation.
297  if ( WT > TMath::Power(m_frag+m_rema+fMinESinglet,2) ) {
298 
299  // Energy of transfered to the fragmented quark and rema system
300  // Applying energy conservation
301  WT = TMath::Sqrt( WT + pT2 );
302  double tm_rema = TMath::Power(m_rema,2) + pT2;
303  double E_frag = 0.5 * ( WT + ( TMath::Power(m_frag,2) - tm_rema)/WT ); //E_frag + E_rema = WT
304  double E_rema = 0.5 * ( WT + (-TMath::Power(m_frag,2) + tm_rema)/WT );
305  double x_rema = -1 * TMath::Sqrt( TMath::Power(E_rema,2) - tm_rema );
306  double theta_rema = pyangl_(&x_rema,&pT);
307 
308  // Select a phi angle between between particles randomly
309  double phi = 2*kPi*rnd->RndHadro().Rndm();
310 
311  // Input the three particles to PYTHIA in the CM frame
312  // If a top quark is produced we decay it because it does not hadronize
313  fPythia->Py1ent( -1, frag_quark, E_frag, 0., 0. ); //k(1,2) = 2
314  if (TMath::Abs(frag_quark) > 5 ) {
315  int ip = 1;
316  pydecy_(&ip);
317  }
318  fPythia->Py1ent( fPythia->GetN()+1, rema, E_rema, theta_rema, phi ); //k(2,2) = 1
319 
320  int imin = 0;
321  int imax = 0;
322  double the = 0.; double ph = 0.;
323  double dbex = 0.; double dbey = 0.; double dbez = (E_pz-(1-z)*W)/(E_pz+(1-z)*W);
324  pyrobo_( &imin , &imax, &the, &ph, &dbex, &dbey , &dbez );
325 
326  double pz_hadron = -0.5 * ( z*W - tm_hadron );
327  double theta_hadron = pyangl_(&pz_hadron,&pT);
328  fPythia->SetMSTU( 10, 1 ); //keep the mass value stored in P(I,5), whatever it is.
329  fPythia->SetP( fPythia->GetN()+1, 5, m_hadron );
330  fPythia->Py1ent( fPythia->GetN()+1, hadron, E_hadron, theta_hadron, phi + kPi );
331  fPythia->SetMSTU( 10, 2 ); //find masses according to mass tables as usual.
332 
333  // Target remnants required to go backwards in hadronic cms
334  if ( fPythia->GetP(fPythia->GetN()-1,3)<0 && fPythia->GetP(fPythia->GetN(),3)<0 ) break; //quit the while from line 368
335 
336  LOG("LeptoHad", pINFO) << "Not backward hadron or rema";
337  LOG("LeptoHad", pINFO) << "hadron = " << hadron << " -> Pz = " << fPythia->GetP(fPythia->GetN(),3) ;
338  LOG("LeptoHad", pINFO) << "rema = " << rema << " -> Pz = " << fPythia->GetP(fPythia->GetN()-1,3) ;
339 
340  }
341  else {
342  LOG("LeptoHad", pINFO) << "Low WT value ... ";
343  LOG("LeptoHad", pINFO) << "WT = " << TMath::Sqrt(WT) << " // m_frag = " << m_frag << " // m_rema = " << m_rema;
344  }
345 
346  LOG("LeptoHad", pINFO) << "Hadronization paricles not suitable. Trying again... " << counter;
347  counter++;
348  if (counter==100) {
349  LOG("LeptoHad", pWARN) << "Hadronization particles failed after " << counter << " iterations! Returning a null list";
350  return 0;
351  }
352 
353  }
354 
355  }
356 
357  // Introduce a primordial kT system
358  double pT = fPrimordialKT * TMath::Sqrt( -1*TMath::Log( rnd->RndHadro().Rndm() ) );
359  double phi = -2*kPi*rnd->RndHadro().Rndm();
360  double theta = 0.;
361  int imin = 0;
362  int imax = 0;
363  double dbex = 0.; double dbey = 0.; double dbez = 0;
364  pyrobo_( &imin , &imax, &theta, &phi, &dbex, &dbey , &dbez );
365  phi = -1 * phi;
366  theta = TMath::ATan(2.*pT/W);
367  pyrobo_( &imin , &imax, &theta, &phi, &dbex, &dbey , &dbez );
368 
369 
370  // Run PYTHIA with the input particles
371  fPythia->Pyexec();
372 
373  fPythia->SetPMAS(24,1,pmas1_W);
374  fPythia->SetPMAS(24,2,pmas2_W);
375  fPythia->SetPMAS(6,2,pmas2_t);
376 
377  // Use for debugging purposes
378  //fPythia->Pylist(3);
379 
380  // get LUJETS record
381  fPythia->GetPrimaries();
382  TClonesArray * pythia_particles = (TClonesArray *) fPythia->ImportParticles("All");
383 
384  // copy PYTHIA container to a new TClonesArray so as to transfer ownership
385  // of the container and of its elements to the calling method
386 
387  int np = pythia_particles->GetEntries();
388  assert(np>0);
389  TClonesArray * particle_list = new TClonesArray("genie::GHepParticle", np);
390  particle_list->SetOwner(true);
391 
392  // Boost velocity HCM -> LAB
393  long double beta = p4Hadlong.P()/p4Hadlong.E();
394 
395  //fix numbering for events with top
396  bool isTop = false;
397 
398  //-- Translate the fragmentation products from TMCParticles to
399  // GHepParticles and copy them to the event record.
400  int mom = event->FinalStateHadronicSystemPosition();
401  assert(mom!=-1);
402 
403  TMCParticle * p = 0;
404  TIter particle_iter(pythia_particles);
405  while( (p = (TMCParticle *) particle_iter.Next()) ) {
406 
407  int pdgc = p->GetKF();
408  int ks = p->GetKS();
409 
410  // Final state particles can not be quarks or diquarks but colorless
411  if(ks == 1) {
412  if( pdg::IsQuark(pdgc) || pdg::IsDiQuark(pdgc) ) {
413  LOG("LeptoHad", pERROR) << "Hadronization failed! Bare quark/di-quarks appear in final state!";
414  return false;
415  }
416  }
417 
418  // When top quark is produced, it is immidiately decay before hadronization. Then the decayed
419  // products are hadronized with the hadron remnants. Therefore, we remove the top quark from
420  // the list of particles so that the mother/daugher assigments is at the same level for decayed
421  // products and hadron remnants.
422  if ( pdg::IsTQuark( TMath::Abs(pdgc) ) ) { isTop=true; continue; }
423 
424  // fix numbering scheme used for mother/daughter assignments
425  if ( isTop ) {
426  (p->GetParent()==0) ? p->SetParent(p->GetParent() - 1) : p->SetParent(p->GetParent() - 2);
427  p->SetFirstChild (p->GetFirstChild() - 2);
428  p->SetLastChild (p->GetLastChild() - 2);
429  }
430  else {
431  p->SetParent(p->GetParent() - 1);
432  p->SetFirstChild (p->GetFirstChild() - 1);
433  p->SetLastChild (p->GetLastChild() - 1);
434  }
435 
436  LongLorentzVector p4long( p->GetPx(), p->GetPy(), p->GetPz(), p->GetEnergy() );
437  p4long.BoostZ(beta);
438  p4long.Rotate(p4Hadlong);
439 
440  // Translate from long double to double
441  TLorentzVector p4( (double)p4long.Px(), (double)p4long.Py(), (double)p4long.Pz(), (double)p4long.E() );
442 
443  // Somtimes PYTHIA output particles with E smaller than its mass. This is wrong,
444  // so we assume that the are at rest.
445  double massPDG = PDGLibrary::Instance()->Find(pdgc)->Mass();
446  if ( (ks==1 || ks==4) && p4.E() < massPDG ) {
447  LOG("LeptoHad", pINFO) << "Putting at rest one stable particle generated by PYTHIA because E < m";
448  LOG("LeptoHad", pINFO) << "PDG = " << pdgc << " // State = " << ks;
449  LOG("LeptoHad", pINFO) << "E = " << p4.E() << " // |p| = " << p4.P();
450  LOG("LeptoHad", pINFO) << "p = [ " << p4.Px() << " , " << p4.Py() << " , " << p4.Pz() << " ]";
451  LOG("LeptoHad", pINFO) << "m = " << p4.M() << " // mpdg = " << massPDG;
452  p4.SetXYZT(0,0,0,massPDG);
453  }
454 
455  // copy final state particles to the event record
457 
458  int im = mom + 1 + p->GetParent();
459  int ifc = (p->GetFirstChild() <= -1) ? -1 : mom + 1 + p->GetFirstChild();
460  int ilc = (p->GetLastChild() <= -1) ? -1 : mom + 1 + p->GetLastChild();
461 
462  double vx = vtx.X() + p->GetVx()*1e12; //pythia gives position in [mm] while genie uses [fm]
463  double vy = vtx.Y() + p->GetVy()*1e12;
464  double vz = vtx.Z() + p->GetVz()*1e12;
465  double vt = vtx.T() + p->GetTime()*(units::millimeter/units::second);
466  TLorentzVector pos( vx, vy, vz, vt );
467 
468  event->AddParticle( pdgc, ist, im,-1, ifc, ilc, p4, pos );
469 
470  }
471  return true;
472 
473 #else
474  return false;
475 #endif // __GENIE_PYTHIA6_ENABLED__
476 
477 }
478 //____________________________________________________________________________
480 {
481  Algorithm::Configure(config);
482  this->LoadConfig();
483 }
484 //____________________________________________________________________________
485 void LeptoHadronization::Configure(string config)
486 {
487  Algorithm::Configure(config);
488  this->LoadConfig();
489 }
490 //____________________________________________________________________________
492 {
493 
494 
495  GetParam("MaxIter-Had", fMaxIterHad ) ;
496 
497  // Width of Gaussian distribution for transverse momentums
498  // Define in LEPTO with PARL(3) and PARL(14)
499  GetParam("Primordial-kT", fPrimordialKT ) ;
500  GetParam("Remnant-pT", fRemnantPT ) ;
501 
502  // It is, with quark masses added, used to define the minimum allowable energy of a colour-singlet parton system.
503  GetParam( "Energy-Singlet", fMinESinglet ) ;
504 
505 #ifdef __GENIE_PYTHIA6_ENABLED__
506  // PYTHIA parameters only valid for HEDIS
507  GetParam( "Xsec-Wmin", fWmin ) ;
508  fPythia->SetPARP(2, fWmin); //(D = 10. GeV) lowest c.m. energy for the event as a whole that the program will accept to simulate. (bellow 2GeV pythia crashes)
509 
510  int warnings; GetParam( "Warnings", warnings ) ;
511  int errors; GetParam( "Errors", errors ) ;
512  int qrk_mass; GetParam( "QuarkMass", qrk_mass ) ;
513 
514  fPythia->SetMSTU(26, warnings); // (Default=10) maximum number of warnings that are printed
515  fPythia->SetMSTU(22, errors); // (Default=10) maximum number of errors that are printed
516  fPythia->SetMSTJ(93, qrk_mass); // light (d, u, s, c, b) quark masses are taken from PARF(101) - PARF(105) rather than PMAS(1,1) - PMAS(5,1). Diquark masses are given as sum of quark masses, without spin splitting term.
517 
518  fPythia->SetPMAS(24,1,kMw); //mass of the W boson (pythia=80.450 // genie=80.385)
519  fPythia->SetPMAS(24,2,0.); //set to 0 the width of the W boson to avoid problems with energy conservation
520  fPythia->SetPMAS(6,2,0.); //set to 0 the width of the top to avoid problems with energy conservation
521 #endif
522 
523 }
TPythia6 * fPythia
PYTHIA6 wrapper class.
void Configure(const Registry &config)
const int kPdgUUDiquarkS1
Definition: PDGCodes.h:58
double W(bool selected=false) const
Definition: Kinematics.cxx:157
#define pERROR
Definition: Messenger.h:59
void pyrobo_(int *, int *, double *, double *, double *, double *, double *)
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
#define __GENIE_PYTHIA6_ENABLED__
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
int HitQrkPdg(void) const
Definition: Target.cxx:242
#define pFATAL
Definition: Messenger.h:56
bool IsTQuark(int pdgc)
Definition: PDGUtils.cxx:291
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
Contains minimal information for tagging exclusive processes.
Definition: XclsTag.h:39
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
void ProcessEventRecord(GHepRecord *event) const
Summary information for an interaction.
Definition: Interaction.h:56
static constexpr double millimeter
Definition: Units.h:41
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
static constexpr double second
Definition: Units.h:37
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
const int kPdgUDDiquarkS1
Definition: PDGCodes.h:57
bool IsDiQuark(int pdgc)
Definition: PDGUtils.cxx:231
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
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void py2ent_(int *, int *, int *, double *)
#define pINFO
Definition: Messenger.h:62
const int kPdgUDDiquarkS0
Definition: PDGCodes.h:56
#define pWARN
Definition: Messenger.h:60
void pydecy_(int *)
int FinalQuarkPdg(void) const
Definition: XclsTag.h:72
void pykfdi_(int *, int *, int *, int *)
TRandom3 & RndHadro(void) const
rnd number generator used by hadronization models
Definition: RandomGen.h:53
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:279
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
bool HitQrkIsSet(void) const
Definition: Target.cxx:292
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
bool IsQuark(int pdgc)
Definition: PDGUtils.cxx:250
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
bool IsDQuark(int pdgc)
Definition: PDGUtils.cxx:271
void SetHadSystP4(const TLorentzVector &p4)
Definition: Kinematics.cxx:307
const InitialState & InitState(void) const
Definition: Interaction.h:69
void pyzdis_(int *, int *, double *, double *)
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:66
double pyangl_(double *, double *)
bool gAbortingInErr
Definition: Messenger.cxx:34
const int kPdgHadronicSyst
Definition: PDGCodes.h:210
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...
bool Hadronize(GHepRecord *event) const
enum genie::EGHepStatus GHepStatus_t
#define pDEBUG
Definition: Messenger.h:63
const int kPdgDDDiquarkS1
Definition: PDGCodes.h:55