PythiaHadronization.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab
9 
10  Changes required to implement the GENIE Boosted Dark Matter module
11  were installed by Josh Berger (Univ. of Wisconsin)
12 */
13 //____________________________________________________________________________
14 
15 #include <RVersion.h>
16 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,15,6)
17 #include <TMCParticle.h>
18 #else
19 #include <TMCParticle6.h>
20 #endif
21 #include <TClonesArray.h>
22 #include <TMath.h>
23 #include <TH1D.h>
24 
38 
39 using namespace genie;
40 using namespace genie::constants;
41 
42 // the actual PYTHIA call
43 extern "C" void py2ent_(int *, int *, int *, double *);
44 
45 //____________________________________________________________________________
47 HadronizationModelBase("genie::PythiaHadronization")
48 {
49  this->Initialize();
50 }
51 //____________________________________________________________________________
53 HadronizationModelBase("genie::PythiaHadronization", config)
54 {
55  this->Initialize();
56 }
57 //____________________________________________________________________________
59 {
60 
61 }
62 //____________________________________________________________________________
64 {
65  fPythia = TPythia6::Instance();
66 
67  // sync GENIE/PYTHIA6 seed number
69 }
70 //____________________________________________________________________________
71 TClonesArray *
73  const Interaction * interaction) const
74 {
75  LOG("PythiaHad", pNOTICE) << "Running PYTHIA hadronizer";
76 
77  if(!this->AssertValidity(interaction)) {
78  LOG("PythiaHad", pERROR) << "Returning a null particle list!";
79  return 0;
80  }
81 
82  // get kinematics / init-state / process-info
83 
84  const Kinematics & kinematics = interaction->Kine();
85  const InitialState & init_state = interaction->InitState();
86  const ProcessInfo & proc_info = interaction->ProcInfo();
87  const Target & target = init_state.Tgt();
88 
89  assert(target.HitQrkIsSet());
90 
91  double W = kinematics.W();
92 
93  int probe = init_state.ProbePdg();
94  int hit_nucleon = target.HitNucPdg();
95  int hit_quark = target.HitQrkPdg();
96  bool from_sea = target.HitSeaQrk();
97 
98  LOG("PythiaHad", pNOTICE)
99  << "Hit nucleon pdgc = " << hit_nucleon << ", W = " << W;
100  LOG("PythiaHad", pNOTICE)
101  << "Selected hit quark pdgc = " << hit_quark
102  << ((from_sea) ? "[sea]" : "[valence]");
103 
104  // check hit-nucleon assignment, input neutrino & interaction type
105  bool isp = pdg::IsProton (hit_nucleon);
106  bool isn = pdg::IsNeutron (hit_nucleon);
107  bool isv = pdg::IsNeutrino (probe);
108  bool isvb = pdg::IsAntiNeutrino (probe);
109 //bool isl = pdg::IsNegChargedLepton (probe);
110 //bool islb = pdg::IsPosChargedLepton (probe);
111  bool iscc = proc_info.IsWeakCC ();
112  bool isnc = proc_info.IsWeakNC ();
113  bool isdm = proc_info.IsDarkMatter ();
114  bool isem = proc_info.IsEM ();
115  bool isu = pdg::IsUQuark (hit_quark);
116  bool isd = pdg::IsDQuark (hit_quark);
117  bool iss = pdg::IsSQuark (hit_quark);
118  bool isub = pdg::IsAntiUQuark (hit_quark);
119  bool isdb = pdg::IsAntiDQuark (hit_quark);
120  bool issb = pdg::IsAntiSQuark (hit_quark);
121 
122  //
123  // Generate the quark system (q + qq) initiating the hadronization
124  //
125 
126  int final_quark = 0; // leading quark (hit quark after the interaction)
127  int diquark = 0; // remnant diquark (xF<0 at hadronic CMS)
128 
129  // Figure out the what happens to the hit quark after the interaction
130  if (isnc || isem || isdm) {
131  // NC, EM
132  final_quark = hit_quark;
133  } else {
134  // CC
135  if (isv && isd ) final_quark = kPdgUQuark;
136  else if (isv && iss ) final_quark = kPdgUQuark;
137  else if (isv && isub) final_quark = kPdgAntiDQuark;
138  else if (isvb && isu ) final_quark = kPdgDQuark;
139  else if (isvb && isdb) final_quark = kPdgAntiUQuark;
140  else if (isvb && issb) final_quark = kPdgAntiUQuark;
141  else {
142  LOG("PythiaHad", pERROR)
143  << "Not allowed mode. Refused to make a final quark assignment!";
144  return 0;
145  }
146  }//CC
147 
148  // Figure out what the remnant diquark is.
149  // Note from Hugh, following a conversation with his local HEP theorist
150  // (Gary Goldstein): "I am told that the probability of finding the diquark
151  // in the singlet vs. triplet states is 50-50."
152 
153  // hit quark = valence quark
154  if(!from_sea) {
155  if (isp && isu) diquark = kPdgUDDiquarkS1; /* u(->q) + ud */
156  if (isp && isd) diquark = kPdgUUDiquarkS1; /* d(->q) + uu */
157  if (isn && isu) diquark = kPdgDDDiquarkS1; /* u(->q) + dd */
158  if (isn && isd) diquark = kPdgUDDiquarkS1; /* d(->q) + ud */
159  }
160  // hit quark = sea quark
161  else {
162  if(isp && isu) diquark = kPdgUDDiquarkS1; /* u(->q) + bar{u} uud (=ud) */
163  if(isp && isd) diquark = kPdgUUDiquarkS1; /* d(->q) + bar{d} uud (=uu) */
164  if(isn && isu) diquark = kPdgDDDiquarkS1; /* u(->q) + bar{u} udd (=dd) */
165  if(isn && isd) diquark = kPdgUDDiquarkS1; /* d(->q) + bar{d} udd (=ud) */
166 
167  // The following section needs revisiting.
168 
169  // The lepton is scattered off a sea antiquark, materializing its quark
170  // partner and leaving me with a 5q system ( <qbar + q> + qqq(valence) )
171  // I will force few qbar+q annhilations below to get my quark/diquark system
172  // Probably it is best to leave the qqq system in the final state and then
173  // just do the fragmentation of the qbar q system? But how do I figure out
174  // how to split the available energy?
175 
176  /* bar{u} (-> bar{d}) + u uud => u + uu */
177  if(isp && isub && iscc) {final_quark = kPdgUQuark; diquark = kPdgUUDiquarkS1;}
178  /* bar{u} (-> bar{u}) + u uud => u + ud */
179  if(isp && isub && (isnc||isem||isdm)) {final_quark = kPdgUQuark; diquark = kPdgUDDiquarkS1;}
180  /* bar{d} (-> bar{u}) + d uud => d + ud */
181  if(isp && isdb && iscc) {final_quark = kPdgDQuark; diquark = kPdgUDDiquarkS1;}
182  /* bar{d} (-> bar{d}) + d uud => d + uu */
183  if(isp && isdb && (isnc||isem||isdm)) {final_quark = kPdgDQuark; diquark = kPdgUUDiquarkS1;}
184  /* bar{u} (-> bar{d}) + u udd => u + ud */
185  if(isn && isub && iscc) {final_quark = kPdgUQuark; diquark = kPdgUDDiquarkS1;}
186  /* bar{u} (-> bar{u}) + u udd => u + dd */
187  if(isn && isub && (isnc||isem||isdm)) {final_quark = kPdgUQuark; diquark = kPdgDDDiquarkS1;}
188  /* bar{d} (-> bar{u}) + d udd => d + dd */
189  if(isn && isdb && iscc) {final_quark = kPdgDQuark; diquark = kPdgDDDiquarkS1;}
190  /* bar{d} (-> bar{d}) + d udd => d + ud */
191  if(isn && isdb && (isnc||isem||isdm)) {final_quark = kPdgDQuark; diquark = kPdgUDDiquarkS1;}
192 
193  // The neutrino is scatterred off s or sbar sea quarks
194  // For the time being I will handle s like d and sbar like dbar (copy & paste
195  // from above) so that I conserve charge.
196 
197  if(iss || issb) {
198  LOG("PythiaHad", pNOTICE)
199  << "Can not really handle a hit s or sbar quark / Faking it";
200 
201  if(isp && iss) { diquark = kPdgUUDiquarkS1; }
202  if(isn && iss) { diquark = kPdgUDDiquarkS1; }
203 
204  if(isp && issb && iscc) {final_quark = kPdgDQuark; diquark = kPdgUDDiquarkS1;}
205  if(isp && issb && (isnc||isem||isdm)) {final_quark = kPdgDQuark; diquark = kPdgUUDiquarkS1;}
206  if(isn && issb && iscc) {final_quark = kPdgDQuark; diquark = kPdgDDDiquarkS1;}
207  if(isn && issb && (isnc||isem||isdm)) {final_quark = kPdgDQuark; diquark = kPdgUDDiquarkS1;}
208  }
209 
210  // if the diquark is a ud, switch it to the singlet state with 50% probability
211  if(diquark == kPdgUDDiquarkS1) {
213  double Rqq = rnd->RndHadro().Rndm();
214  if(Rqq<0.5) diquark = kPdgUDDiquarkS0;
215  }
216  }
217  assert(diquark!=0);
218 
219  //
220  // PYTHIA -> HADRONIZATION
221  //
222 
223  LOG("PythiaHad", pNOTICE)
224  << "Fragmentation / Init System: "
225  << "q = " << final_quark << ", qq = " << diquark;
226  int ip = 0;
227 
228  // Determine how jetset treats un-stable particles appearing in hadronization
229 
230  int pi0_decflag = fPythia->GetMDCY(fPythia->Pycomp(kPdgPi0), 1);
231  int K0_decflag = fPythia->GetMDCY(fPythia->Pycomp(kPdgK0), 1);
232  int K0b_decflag = fPythia->GetMDCY(fPythia->Pycomp(kPdgAntiK0), 1);
233  int L0_decflag = fPythia->GetMDCY(fPythia->Pycomp(kPdgLambda), 1);
234  int L0b_decflag = fPythia->GetMDCY(fPythia->Pycomp(kPdgAntiLambda), 1);
235  int Dm_decflag = fPythia->GetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaM), 1);
236  int D0_decflag = fPythia->GetMDCY(fPythia->Pycomp(kPdgP33m1232_Delta0), 1);
237  int Dp_decflag = fPythia->GetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaP), 1);
238  int Dpp_decflag = fPythia->GetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaPP), 1);
239 
240 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
241  LOG("PythiaHad", pDEBUG) << "Original decay flag for pi0 = " << pi0_decflag;
242  LOG("PythiaHad", pDEBUG) << "Original decay flag for K0 = " << K0_decflag;
243  LOG("PythiaHad", pDEBUG) << "Original decay flag for \bar{K0} = " << K0b_decflag;
244  LOG("PythiaHad", pDEBUG) << "Original decay flag for Lambda = " << L0_decflag;
245  LOG("PythiaHad", pDEBUG) << "Original decay flag for \bar{Lambda0} = " << L0b_decflag;
246  LOG("PythiaHad", pDEBUG) << "Original decay flag for D- = " << Dm_decflag;
247  LOG("PythiaHad", pDEBUG) << "Original decay flag for D0 = " << D0_decflag;
248  LOG("PythiaHad", pDEBUG) << "Original decay flag for D+ = " << Dp_decflag;
249  LOG("PythiaHad", pDEBUG) << "Original decay flag for D++ = " << Dpp_decflag;
250 #endif
251 
252  fPythia->SetMDCY(fPythia->Pycomp(kPdgPi0), 1,0); // don't decay pi0
253  fPythia->SetMDCY(fPythia->Pycomp(kPdgK0), 1,0); // don't decay K0
254  fPythia->SetMDCY(fPythia->Pycomp(kPdgAntiK0), 1,0); // don't decay \bar{K0}
255  fPythia->SetMDCY(fPythia->Pycomp(kPdgLambda), 1,0); // don't decay Lambda0
256  fPythia->SetMDCY(fPythia->Pycomp(kPdgAntiLambda), 1,0); // don't decay \bar{Lambda0}
257  fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaM), 1,1); // decay Delta-
258  fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_Delta0), 1,1); // decay Delta0
259  fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaP), 1,1); // decay Delta+
260  fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaPP), 1,1); // decay Delta++
261 
262  // -- hadronize --
263  py2ent_(&ip, &final_quark, &diquark, &W); // hadronizer
264 
265  // restore pythia decay settings so as not to interfere with decayer
266  fPythia->SetMDCY(fPythia->Pycomp(kPdgPi0), 1, pi0_decflag);
267  fPythia->SetMDCY(fPythia->Pycomp(kPdgK0), 1, K0_decflag);
268  fPythia->SetMDCY(fPythia->Pycomp(kPdgAntiK0), 1, K0b_decflag);
269  fPythia->SetMDCY(fPythia->Pycomp(kPdgLambda), 1, L0_decflag);
270  fPythia->SetMDCY(fPythia->Pycomp(kPdgAntiLambda), 1, L0b_decflag);
271  fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaM), 1, Dm_decflag);
272  fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_Delta0), 1, D0_decflag);
273  fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaP), 1, Dp_decflag);
274  fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaPP), 1, Dpp_decflag);
275 
276  // get LUJETS record
277  fPythia->GetPrimaries();
278  TClonesArray * pythia_particles =
279  (TClonesArray *) fPythia->ImportParticles("All");
280 
281  // copy PYTHIA container to a new TClonesArray so as to transfer ownership
282  // of the container and of its elements to the calling method
283 
284  int np = pythia_particles->GetEntries();
285  assert(np>0);
286  TClonesArray * particle_list = new TClonesArray("TMCParticle", np);
287  particle_list->SetOwner(true);
288 
289  unsigned int i = 0;
290  TMCParticle * particle = 0;
291  TIter particle_iter(pythia_particles);
292 
293  while( (particle = (TMCParticle *) particle_iter.Next()) ) {
294  LOG("PythiaHad", pDEBUG)
295  << "Adding final state particle pdgc = " << particle->GetKF()
296  << " with status = " << particle->GetKS();
297 
298  if(particle->GetKS() == 1) {
299  if( pdg::IsQuark (particle->GetKF()) ||
300  pdg::IsDiQuark(particle->GetKF()) ) {
301  LOG("PythiaHad", pERROR)
302  << "Hadronization failed! Bare quark/di-quarks appear in final state!";
303  particle_list->Delete();
304  delete particle_list;
305  return 0;
306  }
307  }
308 
309  // fix numbering scheme used for mother/daughter assignments
310  particle->SetParent (particle->GetParent() - 1);
311  particle->SetFirstChild (particle->GetFirstChild() - 1);
312  particle->SetLastChild (particle->GetLastChild() - 1);
313 
314  // insert the particle in the list
315  new ( (*particle_list)[i++] ) TMCParticle(*particle);
316  }
317 
318  utils::fragmrec::Print(particle_list);
319  return particle_list;
320 }
321 //____________________________________________________________________________
322 PDGCodeList *
324  const Interaction * interaction) const
325 {
326 // Works the opposite way (compared with the KNO hadronization model)
327 // Rather than having this method as one of the hadronization model components,
328 // we extract the list of particles from the fragmentation record after the
329 // hadronization has been completed.
330 
331  TClonesArray * particle_list = this->Hadronize(interaction);
332 
333  if(!particle_list) return 0;
334 
335  bool allowdup=true;
336  PDGCodeList * pdgcv = new PDGCodeList(allowdup);
337  pdgcv->reserve(particle_list->GetEntries());
338 
339  TMCParticle * particle = 0;
340  TIter particle_iter(particle_list);
341 
342  while ((particle = (TMCParticle *) particle_iter.Next()))
343  {
344  if (particle->GetKS()==1) pdgcv->push_back(particle->GetKF());
345  }
346  particle_list->Delete();
347  delete particle_list;
348 
349  return pdgcv;
350 }
351 //____________________________________________________________________________
353  const Interaction * interaction, Option_t * opt) const
354 {
355 // Similar comments apply as in SelectParticles()
356 
357  if(!this->AssertValidity(interaction)) {
358  LOG("PythiaHad", pWARN)
359  << "Returning a null multipicity probability distribution!";
360  return 0;
361  }
362  double maxmult = this->MaxMult(interaction);
363  TH1D * mult_prob = this->CreateMultProbHist(maxmult);
364 
365  const int nev=500;
366  TMCParticle * particle = 0;
367 
368  for(int iev=0; iev<nev; iev++) {
369 
370  TClonesArray * particle_list = this->Hadronize(interaction);
371  double weight = this->Weight();
372 
373  if(!particle_list) { iev--; continue; }
374 
375  int n = 0;
376  TIter particle_iter(particle_list);
377  while ((particle = (TMCParticle *) particle_iter.Next()))
378  {
379  if (particle->GetKS()==1) n++;
380  }
381  particle_list->Delete();
382  delete particle_list;
383  mult_prob->Fill( (double)n, weight);
384  }
385 
386  double integral = mult_prob->Integral("width");
387  if(integral>0) {
388  // Normalize the probability distribution
389  mult_prob->Scale(1.0/integral);
390  } else {
391  SLOG("PythiaHad", pWARN) << "probability distribution integral = 0";
392  return mult_prob;
393  }
394 
395  string option(opt);
396 
397  bool apply_neugen_Rijk = option.find("+LowMultSuppr") != string::npos;
398  bool renormalize = option.find("+Renormalize") != string::npos;
399 
400  // Apply the NeuGEN probability scaling factors -if requested-
401  if(apply_neugen_Rijk) {
402  SLOG("KNOHad", pINFO) << "Applying NeuGEN scaling factors";
403  // Only do so for W<Wcut
404  const Kinematics & kinematics = interaction->Kine();
405  double W = kinematics.W();
406  if(W<fWcut) {
407  this->ApplyRijk(interaction, renormalize, mult_prob);
408  } else {
409  SLOG("PythiaHad", pDEBUG)
410  << "W = " << W << " < Wcut = " << fWcut
411  << " - Will not apply scaling factors";
412  }//<wcut?
413  }//apply?
414 
415  return mult_prob;
416 }
417 //____________________________________________________________________________
418 double PythiaHadronization::Weight(void) const
419 {
420  return 1.; // does not generate weighted events
421 }
422 //____________________________________________________________________________
424 {
425  Algorithm::Configure(config);
426  this->LoadConfig();
427 }
428 //____________________________________________________________________________
430 {
431  Algorithm::Configure(config);
432  this->LoadConfig();
433 }
434 //____________________________________________________________________________
436 {
437  // the configurable PYTHIA parameters used here are the ones used by NUX
438  // (see A.Rubbia's talk @ NuINT-01)
439  // The defaults are the values used by PYTHIA
440  // Use the NUX config set to set the tuned values as used in NUX.
441 
442  GetParam( "PYTHIA-SSBarSuppression", fSSBarSuppression ) ;
443  GetParam( "PYTHIA-GaussianPt2", fGaussianPt2 ) ;
444  GetParam( "PYTHIA-NonGaussianPt2Tail", fNonGaussianPt2Tail ) ;
445  GetParam( "PYTHIA-RemainingEnergyCutoff", fRemainingECutoff ) ;
446 
447  fPythia->SetPARJ(2, fSSBarSuppression);
448  fPythia->SetPARJ(21, fGaussianPt2);
449  fPythia->SetPARJ(23, fNonGaussianPt2Tail);
450  fPythia->SetPARJ(33, fRemainingECutoff);
451 
452  // Load Wcut determining the phase space area where the multiplicity prob.
453  // scaling factors would be applied -if requested-
454  GetParam( "Wcut", fWcut ) ;
455 
456  // decayer
457  fDecayer = 0;
458  if( GetConfig().Exists("Decayer") ) {
459  fDecayer = dynamic_cast<const DecayModelI *> (this->SubAlg("Decayer"));
460  assert(fDecayer);
461  }
462 
463  // Load NEUGEN multiplicity probability scaling parameters Rijk
464  //neutrinos
465  GetParam( "DIS-HMultWgt-vp-CC-m2", fRvpCCm2 ) ;
466  GetParam( "DIS-HMultWgt-vp-CC-m3", fRvpCCm3 ) ;
467  GetParam( "DIS-HMultWgt-vp-NC-m2", fRvpNCm2 ) ;
468  GetParam( "DIS-HMultWgt-vp-NC-m3", fRvpNCm3 ) ;
469  GetParam( "DIS-HMultWgt-vn-CC-m2", fRvnCCm2 ) ;
470  GetParam( "DIS-HMultWgt-vn-CC-m3", fRvnCCm3 ) ;
471  GetParam( "DIS-HMultWgt-vn-NC-m2", fRvnNCm2 ) ;
472  GetParam( "DIS-HMultWgt-vn-NC-m3", fRvnNCm3 ) ;
473  //Anti-neutrinos
474  GetParam( "DIS-HMultWgt-vbp-CC-m2", fRvbpCCm2 ) ;
475  GetParam( "DIS-HMultWgt-vbp-CC-m3", fRvbpCCm3 ) ;
476  GetParam( "DIS-HMultWgt-vbp-NC-m2", fRvbpNCm2 ) ;
477  GetParam( "DIS-HMultWgt-vbp-NC-m3", fRvbpNCm3 ) ;
478  GetParam( "DIS-HMultWgt-vbn-CC-m2", fRvbnCCm2 ) ;
479  GetParam( "DIS-HMultWgt-vbn-CC-m3", fRvbnCCm3 ) ;
480  GetParam( "DIS-HMultWgt-vbn-NC-m2", fRvbnNCm2 ) ;
481  GetParam( "DIS-HMultWgt-vbn-NC-m3", fRvbnNCm3 ) ;
482 
483  LOG("PythiaHad", pDEBUG) << GetConfig() ;
484 }
485 //____________________________________________________________________________
486 bool PythiaHadronization::AssertValidity(const Interaction * interaction) const
487 {
488  // check that there is no charm production
489  // (GENIE uses a special model for these cases)
490  if(interaction->ExclTag().IsCharmEvent()) {
491  LOG("PythiaHad", pWARN) << "Can't hadronize charm events";
492  return false;
493  }
494  // check the available mass
495  double W = utils::kinematics::W(interaction);
496  if(W < this->Wmin()) {
497  LOG("PythiaHad", pWARN) << "Low invariant mass, W = " << W << " GeV!!";
498  return false;
499  }
500 
501  const InitialState & init_state = interaction->InitState();
502  const ProcessInfo & proc_info = interaction->ProcInfo();
503  const Target & target = init_state.Tgt();
504 
505  if( ! target.HitQrkIsSet() ) {
506  LOG("PythiaHad", pWARN) << "Hit quark was not set!";
507  return false;
508  }
509 
510  int probe = init_state.ProbePdg();
511  int hit_nucleon = target.HitNucPdg();
512  int hit_quark = target.HitQrkPdg();
513 //bool from_sea = target.HitSeaQrk();
514 
515  // check hit-nucleon assignment, input neutrino & weak current
516  bool isp = pdg::IsProton (hit_nucleon);
517  bool isn = pdg::IsNeutron (hit_nucleon);
518  bool isv = pdg::IsNeutrino (probe);
519  bool isvb = pdg::IsAntiNeutrino (probe);
520  bool isdm = pdg::IsDarkMatter (probe);
521  bool isl = pdg::IsNegChargedLepton (probe);
522  bool islb = pdg::IsPosChargedLepton (probe);
523  bool iscc = proc_info.IsWeakCC ();
524  bool isnc = proc_info.IsWeakNC ();
525  bool isdmi = proc_info.IsDarkMatter ();
526  bool isem = proc_info.IsEM ();
527  if( !(iscc||isnc||isem||isdmi) ) {
528  LOG("PythiaHad", pWARN)
529  << "Can only handle electro-weak interactions";
530  return false;
531  }
532  if( !(isp||isn) || !(isv||isvb||isl||islb||isdm) ) {
533  LOG("PythiaHad", pWARN)
534  << "Invalid initial state: probe = "
535  << probe << ", hit_nucleon = " << hit_nucleon;
536  return false;
537  }
538 
539  // assert that the interaction mode is allowed
540  bool isu = pdg::IsUQuark (hit_quark);
541  bool isd = pdg::IsDQuark (hit_quark);
542  bool iss = pdg::IsSQuark (hit_quark);
543  bool isub = pdg::IsAntiUQuark (hit_quark);
544  bool isdb = pdg::IsAntiDQuark (hit_quark);
545  bool issb = pdg::IsAntiSQuark (hit_quark);
546 
547  bool allowed = (iscc && isv && (isd||isub||iss)) ||
548  (iscc && isvb && (isu||isdb||issb)) ||
549  (isnc && (isv||isvb) && (isu||isd||isub||isdb||iss||issb)) ||
550  (isdmi && isdm && (isu||isd||isub||isdb||iss||issb)) ||
551  (isem && (isl||islb) && (isu||isd||isub||isdb||iss||issb));
552  if(!allowed) {
553  LOG("PythiaHad", pWARN)
554  << "Impossible interaction type / probe / hit quark combination!";
555  return false;
556  }
557 
558  return true;
559 }
560 //____________________________________________________________________________
561 /*
562 void PythiaHadronization::SwitchDecays(int pdgc, bool on_off) const
563 {
564  LOG("PythiaHad", pNOTICE)
565  << "Switching " << ((on_off) ? "ON" : "OFF")
566  << " all PYTHIA decay channels for particle = " << pdgc;
567 
568  int flag = (on_off) ? 1 : 0;
569  int kc = fPythia->Pycomp(pdgc);
570  int first_ch = fPythia->GetMDCY(kc,2);
571  int last_ch = fPythia->GetMDCY(kc,2) + fPythia->GetMDCY(kc,3) - 1;
572 
573  for(int ich = first_ch; ich < last_ch; ich++) fPythia->SetMDME(ich,1,flag);
574 }
575 */
576 //____________________________________________________________________________
577 /*
578 void PythiaHadronization::HandleDecays(TClonesArray * plist) const
579 {
580 // Handle decays of unstable particles if requested through the XML config.
581 // The default is not to decay the particles at this stage (during event
582 // generation, the UnstableParticleDecayer event record visitor decays what
583 // is needed to be decayed later on). But, when comparing various models
584 // (eg PYTHIA vs KNO) independently and not within the full MC simulation
585 // framework it might be necessary to force the decays at this point.
586 
587  if(!fDecayer) {
588  LOG("PythiaHad", pWARN) << "No decayer was specified!";
589  return;
590  }
591 
592  this->SwitchDecays(kPdgLambda, true); // decay Lambda
593  this->SwitchDecays(kPdgAntiLambda, true); // decay \bar{Lambda}
594  this->SwitchDecays(kPdgSigmaP, true); // decay Sigma+
595  this->SwitchDecays(kPdgSigma0, true); // decay Sigma0
596  this->SwitchDecays(kPdgSigmaM, true); // decay Sigma-
597  this->SwitchDecays(kPdgAntiSigmaP, true); // decay Sigma+
598  this->SwitchDecays(kPdgAntiSigma0, true); // decay Sigma0
599  this->SwitchDecays(kPdgAntiSigmaM, true); // decay Sigma-
600  this->SwitchDecays(kPdgXi0, true); // decay Xi0
601  this->SwitchDecays(kPdgXiM, true); // decay Xi-
602  this->SwitchDecays(kPdgAntiXi0, true); // decay \bar{Xi0}
603  this->SwitchDecays(kPdgAntiXiP, true); // decay \bar{Xi+}
604  this->SwitchDecays(kPdgOmegaM, true); // decay Omega-
605  this->SwitchDecays(kPdgAntiOmegaP, true); // decay \bar{Omega+}
606 
607  int mstj21 = fPythia->GetMSTJ(21);
608  fPythia->SetMSTJ(21,1);
609  fPythia->SetMSTJ(22,2);
610  fPythia->SetPARJ(71,100);
611 
612  //-- loop through the fragmentation event record & decay unstables
613  int idecaying = -1; // position of decaying particle
614  TMCParticle * p = 0; // current particle
615 
616  TIter piter(plist);
617  while ( (p = (TMCParticle *) piter.Next()) ) {
618  idecaying++;
619  int status = p->GetKS();
620  int pdg = p->GetKF();
621 
622  bool decay_it = (status<10) &&
623  ( pdg == kPdgLambda ||
624  pdg == kPdgAntiLambda ||
625  pdg == kPdgSigmaP ||
626  pdg == kPdgSigma0 ||
627  pdg == kPdgSigmaM ||
628  pdg == kPdgAntiSigmaP ||
629  pdg == kPdgAntiSigma0 ||
630  pdg == kPdgAntiSigmaM ||
631  pdg == kPdgXi0 ||
632  pdg == kPdgXiM ||
633  pdg == kPdgAntiXi0 ||
634  pdg == kPdgAntiXiP ||
635  pdg == kPdgOmegaM ||
636  pdg == kPdgAntiOmegaP );
637 
638  // bother for final state particle only
639  if(decay_it) {
640 
641  LOG("PythiaHad", pINFO)
642  << "Decaying particle with pdgc = " << p->GetKF();
643 
644  DecayerInputs_t dinp;
645 
646  TLorentzVector p4;
647  p4.SetPxPyPzE(p->GetPx(), p->GetPy(), p->GetPz(), p->GetEnergy());
648 
649  dinp.PdgCode = p->GetKF();
650  dinp.P4 = &p4;
651 
652  TClonesArray * decay_products = fDecayer->Decay(dinp);
653  if(decay_products) {
654  //-- mark the parent particle as decayed & set daughters
655  p->SetKS(11);
656 
657  int nfp = plist->GetEntries(); // n. fragm. products
658  int ndp = decay_products->GetEntries(); // n. decay products
659 
660  p->SetFirstChild ( nfp ); // decay products added at
661  p->SetLastChild ( nfp + ndp -1 ); // the end of the fragm.rec.
662 
663  //-- add decay products to the fragmentation record
664  TMCParticle * dp = 0;
665  TIter dpiter(decay_products);
666 
667  while ( (dp = (TMCParticle *) dpiter.Next()) ) {
668  if(dp->GetKS()>10) continue;
669  dp->SetParent(idecaying);
670  new ( (*plist)[plist->GetEntries()] ) TMCParticle(*dp);
671  }
672 
673  //-- clean up decay products
674  decay_products->Delete();
675  delete decay_products;
676  }
677 
678  } // KS < 10 : final state particle (as in PYTHIA LUJETS record)
679  } // particles in fragmentation record
680 
681  fPythia->SetMSTJ(21,mstj21); // restore mstj(21)
682 }
683 */
684 //____________________________________________________________________________
685 
int iev
Definition: runWimpSim.h:118
const int kPdgP33m1232_DeltaPP
Definition: PDGCodes.h:91
double fRvnCCm3
neugen&#39;s Rijk: vn, CC, multiplicity = 3
PDGCodeList * SelectParticles(const Interaction *) const
const int kPdgUUDiquarkS1
Definition: PDGCodes.h:58
double W(bool selected=false) const
Definition: Kinematics.cxx:167
double fRvbpCCm3
neugen&#39;s Rijk: vbp, CC, multiplicity = 3
Basic constants.
double fSSBarSuppression
ssbar suppression
bool HitSeaQrk(void) const
Definition: Target.cxx:316
bool IsWeakCC(void) const
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:108
const XML_Char * target
Definition: expat.h:268
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
const int kPdgLambda
Definition: PDGCodes.h:69
double MaxMult(const Interaction *i) const
bool IsUQuark(int pdgc)
Definition: PDGUtils.cxx:249
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
int HitNucPdg(void) const
Definition: Target.cxx:321
double fRvnNCm2
neugen&#39;s Rijk: vn, NC, multiplicity = 2
double fRvbpCCm2
neugen&#39;s Rijk: vbp, CC, multiplicity = 2
const Var weight
int HitQrkPdg(void) const
Definition: Target.cxx:259
double fRvpCCm2
neugen&#39;s Rijk: vp, CC, multiplicity = 2
const int kPdgUQuark
Definition: PDGCodes.h:42
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
const DecayModelI * fDecayer
bool IsDarkMatter(int pdgc)
Definition: PDGUtils.cxx:125
double fNonGaussianPt2Tail
non gaussian pt2 tail parameterization
Definition: config.py:1
TString ip
Definition: loadincs.C:5
bool IsSQuark(int pdgc)
Definition: PDGUtils.cxx:259
void py2ent_(int *, int *, int *, double *)
bool IsAntiSQuark(int pdgc)
Definition: PDGUtils.cxx:279
double fRvbnNCm3
neugen&#39;s Rijk: vbn, NC, multiplicity = 3
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
bool IsAntiDQuark(int pdgc)
Definition: PDGUtils.cxx:274
virtual const Registry & GetConfig(void) const
Definition: Algorithm.cxx:254
A list of PDG codes.
Definition: PDGCodeList.h:33
const int kPdgP33m1232_DeltaP
Definition: PDGCodes.h:90
double fRvpNCm2
neugen&#39;s Rijk: vp, NC, multiplicity = 2
const int kPdgP33m1232_DeltaM
Definition: PDGCodes.h:88
const int kPdgK0
Definition: PDGCodes.h:151
bool IsCharmEvent(void) const
Definition: XclsTag.h:48
double W(const Interaction *const i)
Definition: KineUtils.cxx:1015
TH1D * CreateMultProbHist(double maxmult) const
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:304
bool IsPosChargedLepton(int pdgc)
Definition: PDGUtils.cxx:140
Summary information for an interaction.
Definition: Interaction.h:56
const int kPdgAntiUQuark
Definition: PDGCodes.h:43
double fGaussianPt2
gaussian pt2 distribution width
int iscc
Pure abstract base class. Defines the DecayModelI interface to be implemented by any algorithmic clas...
Definition: DecayModelI.h:31
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
void Configure(const Registry &config)
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:97
double fRvpNCm3
neugen&#39;s Rijk: vp, NC, multiplicity = 3
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:44
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:116
const int kPdgAntiDQuark
Definition: PDGCodes.h:45
bool IsDiQuark(int pdgc)
Definition: PDGUtils.cxx:223
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:41
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
int ProbePdg(void) const
Definition: InitialState.h:65
double fRvbnNCm2
neugen&#39;s Rijk: vbn, NC, multiplicity = 2
const int kPdgPi0
Definition: PDGCodes.h:137
const int kPdgDQuark
Definition: PDGCodes.h:44
#define pINFO
Definition: Messenger.h:63
const int kPdgAntiK0
Definition: PDGCodes.h:152
const int kPdgUDDiquarkS0
Definition: PDGCodes.h:56
double fRvpCCm3
neugen&#39;s Rijk: vp, CC, multiplicity = 3
double fRvbpNCm2
neugen&#39;s Rijk: vbp, NC, multiplicity = 2
#define pWARN
Definition: Messenger.h:61
const int kPdgP33m1232_Delta0
Definition: PDGCodes.h:89
bool IsEM(void) const
double fRvnNCm3
neugen&#39;s Rijk: vn, NC, multiplicity = 3
void ApplyRijk(const Interaction *i, bool norm, TH1D *mp) const
TRandom3 & RndHadro(void) const
rnd number generator used by hadronization models
Definition: RandomGen.h:54
double Wmin(void) const
Various utility methods common to hadronization models.
double fRemainingECutoff
remaining E cutoff for stopping fragmentation
bool HitQrkIsSet(void) const
Definition: Target.cxx:309
An abstract class. It avoids implementing the HadronizationModelI interface, leaving it for the concr...
bool IsDarkMatter(void) const
TPythia6 * fPythia
PYTHIA6 wrapper class.
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
bool IsQuark(int pdgc)
Definition: PDGUtils.cxx:233
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
bool IsDQuark(int pdgc)
Definition: PDGUtils.cxx:254
TH1D * MultiplicityProb(const Interaction *, Option_t *opt="") const
void Print(const TClonesArray *const part_list)
assert(nhit_max >=nhit_nbins)
double fWcut
configuration data common to all hadronizers
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double fRvbnCCm3
neugen&#39;s Rijk: vbn, CC, multiplicity = 3
bool AssertValidity(const Interaction *i) const
#define pNOTICE
Definition: Messenger.h:62
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:67
const int kPdgAntiLambda
Definition: PDGCodes.h:70
TClonesArray * Hadronize(const Interaction *) const
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:85
void kinematics()
Definition: kinematics.C:10
#define W(x)
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
bool IsNegChargedLepton(int pdgc)
Definition: PDGUtils.cxx:131
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:67
double fRvbpNCm3
neugen&#39;s Rijk: vbp, NC, multiplicity = 3
bool IsAntiUQuark(int pdgc)
Definition: PDGUtils.cxx:269
double fRvnCCm2
neugen&#39;s Rijk: vn, CC, multiplicity = 2
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
double fRvbnCCm2
neugen&#39;s Rijk: vbn, CC, multiplicity = 2
const int kPdgDDDiquarkS1
Definition: PDGCodes.h:55
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353