GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Private Member Functions | Private Attributes | List of all members
genie::LeptoHadronization Class Reference

Provides access to the LEPTO hadronization models.
. More...

#include <LeptoHadronization.h>

Inheritance diagram for genie::LeptoHadronization:
Inheritance graph
[legend]
Collaboration diagram for genie::LeptoHadronization:
Collaboration graph
[legend]

Public Member Functions

 LeptoHadronization ()
 
 LeptoHadronization (string config)
 
virtual ~LeptoHadronization ()
 
void ProcessEventRecord (GHepRecord *event) const
 
void Configure (const Registry &config)
 
void Configure (string config)
 
- Public Member Functions inherited from genie::EventRecordVisitorI
virtual ~EventRecordVisitorI ()
 
- Public Member Functions inherited from genie::Algorithm
virtual ~Algorithm ()
 
virtual void FindConfig (void)
 
virtual const RegistryGetConfig (void) const
 
RegistryGetOwnedConfig (void)
 
virtual const AlgIdId (void) const
 Get algorithm ID. More...
 
virtual AlgStatus_t GetStatus (void) const
 Get algorithm status. More...
 
virtual bool AllowReconfig (void) const
 
virtual AlgCmp_t Compare (const Algorithm *alg) const
 Compare with input algorithm. More...
 
virtual void SetId (const AlgId &id)
 Set algorithm ID. More...
 
virtual void SetId (string name, string config)
 
const AlgorithmSubAlg (const RgKey &registry_key) const
 
void AdoptConfig (void)
 
void AdoptSubstructure (void)
 
virtual void Print (ostream &stream) const
 Print algorithm info. More...
 

Private Member Functions

bool Hadronize (GHepRecord *event) const
 
void Initialize (void) const
 
void LoadConfig (void)
 

Private Attributes

TPythia6 * fPythia
 PYTHIA6 wrapper class. More...
 
int fMaxIterHad
 
double fPrimordialKT
 
double fRemnantPT
 
double fMinESinglet
 
bool fPromptPythiaList
 
double fWmin
 

Additional Inherited Members

- Static Public Member Functions inherited from genie::Algorithm
static string BuildParamVectKey (const std::string &comm_name, unsigned int i)
 
static string BuildParamVectSizeKey (const std::string &comm_name)
 
static string BuildParamMatKey (const std::string &comm_name, unsigned int i, unsigned int j)
 
static string BuildParamMatRowSizeKey (const std::string &comm_name)
 
static string BuildParamMatColSizeKey (const std::string &comm_name)
 
- Protected Member Functions inherited from genie::EventRecordVisitorI
 EventRecordVisitorI ()
 
 EventRecordVisitorI (string name)
 
 EventRecordVisitorI (string name, string config)
 
- Protected Member Functions inherited from genie::Algorithm
 Algorithm ()
 
 Algorithm (string name)
 
 Algorithm (string name, string config)
 
void Initialize (void)
 
void DeleteConfig (void)
 
void DeleteSubstructure (void)
 
RegistryExtractLocalConfig (const Registry &in) const
 
RegistryExtractLowerConfig (const Registry &in, const string &alg_key) const
 Split an incoming configuration Registry into a block valid for the sub-algo identified by alg_key. More...
 
template<class T >
bool GetParam (const RgKey &name, T &p, bool is_top_call=true) const
 
template<class T >
bool GetParamDef (const RgKey &name, T &p, const T &def) const
 
template<class T >
int GetParamVect (const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
 Handle to load vectors of parameters. More...
 
int GetParamVectKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const
 
template<class T >
int GetParamMat (const std::string &comm_name, TMatrixT< T > &mat, bool is_top_call=true) const
 Handle to load matrix of parameters. More...
 
template<class T >
int GetParamMatSym (const std::string &comm_name, TMatrixTSym< T > &mat, bool is_top_call=true) const
 
int GetParamMatKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const
 
int AddTopRegistry (Registry *rp, bool owns=true)
 add registry with top priority, also update ownership More...
 
int AddLowRegistry (Registry *rp, bool owns=true)
 add registry with lowest priority, also update ownership More...
 
int MergeTopRegistry (const Registry &r)
 
int AddTopRegisties (const vector< Registry * > &rs, bool owns=false)
 Add registries with top priority, also udated Ownerships. More...
 
- Protected Attributes inherited from genie::Algorithm
bool fAllowReconfig
 
bool fOwnsSubstruc
 true if it owns its substructure (sub-algs,...) More...
 
AlgId fID
 algorithm name and configuration set More...
 
vector< Registry * > fConfVect
 
vector< bool > fOwnerships
 ownership for every registry in fConfVect More...
 
AlgStatus_t fStatus
 algorithm execution status More...
 
AlgMapfOwnedSubAlgMp
 local pool for owned sub-algs (taken out of the factory pool) More...
 

Detailed Description

Provides access to the LEPTO hadronization models.
.

Author
Alfonso Garcia <alfonsog nikhef.nl> NIKHEF (Amsterdam)
Created:
October 18, 2019
License:
Copyright (c) 2003-2024, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org or see $GENIE/LICENSE

Definition at line 38 of file LeptoHadronization.h.

Constructor & Destructor Documentation

LeptoHadronization::LeptoHadronization ( )

Definition at line 52 of file LeptoHadronization.cxx.

References Initialize().

52  :
53 EventRecordVisitorI("genie::LeptoHadronization")
54 {
55  this->Initialize();
56 }
LeptoHadronization::LeptoHadronization ( string  config)

Definition at line 58 of file LeptoHadronization.cxx.

References Initialize().

58  :
59 EventRecordVisitorI("genie::LeptoHadronization", config)
60 {
61  this->Initialize();
62 }
LeptoHadronization::~LeptoHadronization ( )
virtual

Definition at line 64 of file LeptoHadronization.cxx.

65 {
66 
67 }

Member Function Documentation

void LeptoHadronization::Configure ( const Registry config)
virtual

Configure the algorithm with an external registry The registry is merged with the top level registry if it is owned, Otherwise a copy of it is added with the highest priority

Reimplemented from genie::Algorithm.

Definition at line 479 of file LeptoHadronization.cxx.

References genie::Algorithm::Configure(), and LoadConfig().

480 {
481  Algorithm::Configure(config);
482  this->LoadConfig();
483 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void LeptoHadronization::Configure ( string  config)
virtual

Configure the algorithm from the AlgoConfigPool based on param_set string given in input An algorithm contains a vector of registries coming from different xml configuration files, which are loaded according a very precise prioriy This methods will load a number registries in order of priority: 1) "Tunable" parameter set from CommonParametes. This is loaded with the highest prioriry and it is designed to be used for tuning procedure Usage not expected from the user. 2) For every string defined in "CommonParame" the corresponding parameter set will be loaded from CommonParameter.xml 3) parameter set specified by the config string and defined in the xml file of the algorithm 4) if config is not "Default" also the Default parameter set from the same xml file will be loaded Effectively this avoids the repetion of a parameter when it is not changed in the requested configuration

Reimplemented from genie::Algorithm.

Definition at line 485 of file LeptoHadronization.cxx.

References genie::Algorithm::Configure(), and LoadConfig().

486 {
487  Algorithm::Configure(config);
488  this->LoadConfig();
489 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
bool LeptoHadronization::Hadronize ( GHepRecord event) const
private

Definition at line 107 of file LeptoHadronization.cxx.

References genie::utils::math::LongLorentzVector::BoostZ(), genie::Interaction::ExclTag(), genie::XclsTag::FinalQuarkPdg(), genie::PDGLibrary::Find(), fMaxIterHad, fMinESinglet, fPrimordialKT, fPythia, fRemnantPT, fWmin, genie::Target::HitNucPdg(), genie::Target::HitQrkIsSet(), genie::Target::HitQrkPdg(), genie::Interaction::InitState(), genie::RandomGen::Instance(), genie::PDGLibrary::Instance(), genie::pdg::IsDiQuark(), genie::pdg::IsDQuark(), genie::pdg::IsProton(), genie::pdg::IsQuark(), genie::pdg::IsTQuark(), genie::pdg::IsUQuark(), genie::Interaction::Kine(), genie::Interaction::KinePtr(), genie::kIStDISPreFragmHadronicState, genie::kIStStableFinalState, genie::constants::kMw, genie::kPdgDDDiquarkS1, genie::kPdgHadronicSyst, genie::kPdgUDDiquarkS0, genie::kPdgUDDiquarkS1, genie::kPdgUUDiquarkS1, genie::constants::kPi, LOG, genie::units::millimeter, pDEBUG, pERROR, pINFO, pWARN, pyangl_(), pydecy_(), pykfdi_(), pyrobo_(), pyzdis_(), genie::RandomGen::RndHadro(), genie::units::second, genie::Kinematics::SetHadSystP4(), genie::Kinematics::SetW(), genie::InitialState::Tgt(), genie::Kinematics::W(), and genie::utils::kinematics::W().

Referenced by ProcessEventRecord().

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 }
TPythia6 * fPythia
PYTHIA6 wrapper class.
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
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
int HitQrkPdg(void) const
Definition: Target.cxx:242
virtual int HitNucleonPosition(void) const
Definition: GHepRecord.cxx:418
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
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:284
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
Summary information for an interaction.
Definition: Interaction.h:56
static constexpr double millimeter
Definition: Units.h:41
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
virtual GHepParticle * FinalStatePrimaryLepton(void) const
Definition: GHepRecord.cxx:333
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
#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
bool IsQuark(int pdgc)
Definition: PDGUtils.cxx:250
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
virtual GHepParticle * HitNucleon(void) const
Definition: GHepRecord.cxx:313
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
const Target & Tgt(void) const
Definition: InitialState.h:66
double pyangl_(double *, double *)
const int kPdgHadronicSyst
Definition: PDGCodes.h:210
enum genie::EGHepStatus GHepStatus_t
#define pDEBUG
Definition: Messenger.h:63
const int kPdgDDDiquarkS1
Definition: PDGCodes.h:55
void LeptoHadronization::Initialize ( void  ) const
private

Definition at line 69 of file LeptoHadronization.cxx.

References fPythia, and genie::RandomGen::Instance().

Referenced by LeptoHadronization().

70 {
71 #ifdef __GENIE_PYTHIA6_ENABLED__
72  fPythia = TPythia6::Instance();
73 
74  // sync GENIE/PYTHIA6 seed number
76 #endif
77 }
TPythia6 * fPythia
PYTHIA6 wrapper class.
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
void LeptoHadronization::LoadConfig ( void  )
private

Definition at line 491 of file LeptoHadronization.cxx.

References fMaxIterHad, fMinESinglet, fPrimordialKT, fPythia, fRemnantPT, fWmin, genie::Algorithm::GetParam(), and genie::constants::kMw.

Referenced by Configure().

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.
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
void LeptoHadronization::ProcessEventRecord ( GHepRecord event) const
virtual

Implements genie::EventRecordVisitorI.

Definition at line 79 of file LeptoHadronization.cxx.

References genie::gAbortingInErr, Hadronize(), genie::kHadroSysGenErr, LOG, pFATAL, pWARN, genie::exceptions::EVGThreadException::SetReason(), and genie::exceptions::EVGThreadException::SwitchOnFastForward().

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 }
#define pFATAL
Definition: Messenger.h:56
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
#define pWARN
Definition: Messenger.h:60
bool gAbortingInErr
Definition: Messenger.cxx:34
bool Hadronize(GHepRecord *event) const

Member Data Documentation

int genie::LeptoHadronization::fMaxIterHad
private

Definition at line 65 of file LeptoHadronization.h.

Referenced by Hadronize(), and LoadConfig().

double genie::LeptoHadronization::fMinESinglet
private

Definition at line 68 of file LeptoHadronization.h.

Referenced by Hadronize(), and LoadConfig().

double genie::LeptoHadronization::fPrimordialKT
private

Definition at line 66 of file LeptoHadronization.h.

Referenced by Hadronize(), and LoadConfig().

bool genie::LeptoHadronization::fPromptPythiaList
private

Definition at line 69 of file LeptoHadronization.h.

TPythia6* genie::LeptoHadronization::fPythia
mutableprivate

PYTHIA6 wrapper class.

Definition at line 61 of file LeptoHadronization.h.

Referenced by Hadronize(), Initialize(), and LoadConfig().

double genie::LeptoHadronization::fRemnantPT
private

Definition at line 67 of file LeptoHadronization.h.

Referenced by Hadronize(), and LoadConfig().

double genie::LeptoHadronization::fWmin
private

Definition at line 70 of file LeptoHadronization.h.

Referenced by Hadronize(), and LoadConfig().


The documentation for this class was generated from the following files: