KNOHadronization.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  Hugh Gallagher <gallag \at minos.phy.tufts.edu>
11  Tufts University
12 
13  Tinjun Yang <tjyang \at stanford.edu>
14  Stanford University
15 
16  Strange baryon production, and adjusted hadronic shower production
17  to conserve strangeness, and to continue balancing charge and
18  maintaining correct multiplicity was implemented by Keith Hofmann
19  and Hugh Gallagher (Tufts)
20 
21  Production of etas was added by Ji Liu (W&M)
22 
23  Changes required to implement the GENIE Boosted Dark Matter module
24  were installed by Josh Berger (Univ. of Wisconsin)
25 */
26 //____________________________________________________________________________
27 
28 #include <cstdlib>
29 
30 #include <RVersion.h>
31 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,15,6)
32 #include <TMCParticle.h>
33 #else
34 #include <TMCParticle6.h>
35 #endif
36 #include <TSystem.h>
37 #include <TLorentzVector.h>
38 #include <TClonesArray.h>
39 #include <TH1D.h>
40 #include <TMath.h>
41 #include <TF1.h>
42 #include <TROOT.h>
43 
53 //#include "Framework/Numerical/Spline.h"
60 
61 using namespace genie;
62 using namespace genie::constants;
63 using namespace genie::controls;
64 using namespace genie::utils::print;
65 
66 //____________________________________________________________________________
68 HadronizationModelBase("genie::KNOHadronization")
69 {
70  fBaryonXFpdf = 0;
71  fBaryonPT2pdf = 0;
72 //fKNO = 0;
73 }
74 //____________________________________________________________________________
76 HadronizationModelBase("genie::KNOHadronization", config)
77 {
78  fBaryonXFpdf = 0;
79  fBaryonPT2pdf = 0;
80 //fKNO = 0;
81 }
82 //____________________________________________________________________________
84 {
85  if (fBaryonXFpdf ) delete fBaryonXFpdf;
86  if (fBaryonPT2pdf) delete fBaryonPT2pdf;
87 //if (fKNO ) delete fKNO;
88 }
89 //____________________________________________________________________________
90 // HadronizationModelI interface implementation:
91 //____________________________________________________________________________
93 {
94 
95 }
96 //____________________________________________________________________________
98  const Interaction * interaction) const
99 {
100 // Generate the hadronic system in a neutrino interaction using a KNO-based
101 // model.
102 
103  if(!this->AssertValidity(interaction)) {
104  LOG("KNOHad", pWARN) << "Returning a null particle list!";
105  return 0;
106  }
107  fWeight=1;
108 
109  double W = utils::kinematics::W(interaction);
110  LOG("KNOHad", pINFO) << "W = " << W << " GeV";
111 
112  //-- Select hadronic shower particles
113  PDGCodeList * pdgcv = this->SelectParticles(interaction);
114 
115  if(!pdgcv) {
116  LOG("KNOHad", pNOTICE)
117  << "Failed selecting particles for " << *interaction;
118  return 0;
119  }
120 
121  //-- Decay the hadronic final state
122  // Two strategies are considered (for N particles):
123  // 1- N (>=2) particles get passed to the phase space decayer. This is the
124  // old NeuGEN strategy.
125  // 2- decay strategy adopted at the July-2006 hadronization model mini-workshop
126  // (C.Andreopoulos, H.Gallagher, P.Kehayias, T.Yang)
127  // The generated baryon P4 gets selected from from experimental xF and pT^2
128  // distributions and the remaining N-1 particles are passed to the phase space
129  // decayer, with P4 = P4(Sum_Hadronic) - P4(Baryon).
130  // For N=2, generate a phase space decay and keep the solution according to its
131  // likelihood calculated based on the baryon xF and pT pdfs. Especially for N=2
132  // keep the option of using simple phase space decay with reweighting switched
133  // off (for consistency with the neugen/daikon version).
134  //
135  TClonesArray * particle_list = 0;
136  bool reweight_decays = fReWeightDecays;
138  bool use_isotropic_decay = (pdgcv->size()==2 && fUseIsotropic2BDecays);
139  if(use_isotropic_decay) {
140  particle_list = this->DecayMethod1(W,*pdgcv,false);
141  } else {
142  particle_list = this->DecayMethod2(W,*pdgcv,reweight_decays);
143  }
144  } else {
145  particle_list = this->DecayMethod1(W,*pdgcv,reweight_decays);
146  }
147 
148  if(!particle_list) {
149  LOG("KNOHad", pNOTICE)
150  << "Failed decaying a hadronic system @ W=" << W
151  << "with multiplicity=" << pdgcv->size();
152 
153  // clean-up and exit
154  delete pdgcv;
155  return 0;
156  }
157 
158  //-- Handle unstable particle decays (if requested)
159  this->HandleDecays(particle_list);
160 
161  //-- The container 'owns' its elements
162  particle_list->SetOwner(true);
163 
164  delete pdgcv;
165 
166  return particle_list;
167 }
168 //____________________________________________________________________________
170  const Interaction * interaction) const
171 {
172  if(!this->AssertValidity(interaction)) {
173  LOG("KNOHad", pWARN) << "Returning a null particle list!";
174  return 0;
175  }
176 
177  unsigned int min_mult = 2;
178  unsigned int mult = 0;
179  PDGCodeList * pdgcv = 0;
180 
181  double W = utils::kinematics::W(interaction);
182 
183  //-- Get the charge that the hadron shower needs to have so as to
184  // conserve charge in the interaction
185  int maxQ = this->HadronShowerCharge(interaction);
186  LOG("KNOHad", pINFO) << "Hadron Shower Charge = " << maxQ;
187 
188  //-- Build the multiplicity probabilities for the input interaction
189  LOG("KNOHad", pDEBUG) << "Building Multiplicity Probability distribution";
190  LOG("KNOHad", pDEBUG) << *interaction;
191  Option_t * opt = "+LowMultSuppr+Renormalize";
192  TH1D * mprob = this->MultiplicityProb(interaction,opt);
193 
194  if(!mprob) {
195  LOG("KNOHad", pWARN) << "Null multiplicity probability distribution!";
196  return 0;
197  }
198  if(mprob->Integral("width")<=0) {
199  LOG("KNOHad", pWARN) << "Empty multiplicity probability distribution!";
200  delete mprob;
201  return 0;
202  }
203 
204  //----- FIND AN ALLOWED SOLUTION FOR THE HADRONIC FINAL STATE
205 
206  bool allowed_state=false;
207  unsigned int itry = 0;
208 
209  while(!allowed_state)
210  {
211  itry++;
212 
213  //-- Go in error if a solution has not been found after many attempts
214  if(itry>kMaxKNOHadSystIterations) {
215  LOG("KNOHad", pERROR)
216  << "Couldn't select hadronic shower particles after: "
217  << itry << " attempts!";
218  delete mprob;
219  return 0;
220  }
221 
222  //-- Generate a hadronic multiplicity
223  mult = TMath::Nint( mprob->GetRandom() );
224 
225  LOG("KNOHad", pINFO) << "Hadron multiplicity = " << mult;
226 
227  //-- Check that the generated multiplicity is consistent with the charge
228  // that the hadronic shower is required to have - else retry
229  if(mult < (unsigned int) TMath::Abs(maxQ)) {
230  LOG("KNOHad", pWARN)
231  << "Multiplicity not enough to generate hadronic charge! Retrying.";
232  allowed_state = false;
233  continue;
234  }
235 
236  //-- Force a min multiplicity
237  // This should never happen if the multiplicity probability distribution
238  // was properly built
239  if(mult < min_mult) {
240  if(fForceMinMult) {
241  LOG("KNOHad", pWARN)
242  << "Low generated multiplicity: " << mult
243  << ". Forcing to minimum accepted multiplicity: " << min_mult;
244  mult = min_mult;
245  } else {
246  LOG("KNOHad", pWARN)
247  << "Generated multiplicity: " << mult << " is too low! Quitting";
248  delete mprob;
249  return 0;
250  }
251  }
252 
253  //-- Determine what kind of particles we have in the final state
254  pdgcv = this->GenerateHadronCodes(mult, maxQ, W);
255 
256  LOG("KNOHad", pNOTICE)
257  << "Generated multiplicity (@ W = " << W << "): " << pdgcv->size();
258 
259  // muliplicity might have been forced to smaller value if the invariant
260  // mass of the hadronic system was not sufficient
261  mult = pdgcv->size(); // update for potential change
262 
263  // is it an allowed decay?
264  double msum=0;
265  vector<int>::const_iterator pdg_iter;
266  for(pdg_iter = pdgcv->begin(); pdg_iter != pdgcv->end(); ++pdg_iter) {
267  int pdgc = *pdg_iter;
268  double m = PDGLibrary::Instance()->Find(pdgc)->Mass();
269 
270  msum += m;
271  LOG("KNOHad", pDEBUG) << "- PDGC=" << pdgc << ", m=" << m << " GeV";
272  }
273  bool permitted = (W > msum);
274 
275  if(!permitted) {
276  LOG("KNOHad", pWARN) << "*** Decay forbidden by kinematics! ***";
277  LOG("KNOHad", pWARN) << "sum{mass} = " << msum << ", W = " << W;
278  LOG("KNOHad", pWARN) << "Discarding hadronic system & re-trying!";
279  delete pdgcv;
280  allowed_state = false;
281  continue;
282  }
283 
284  allowed_state = true;
285 
286  LOG("KNOHad", pNOTICE)
287  << "Found an allowed hadronic state @ W=" << W
288  << " multiplicity=" << mult;
289 
290  } // attempts
291 
292  delete mprob;
293 
294  return pdgcv;
295 }
296 //____________________________________________________________________________
298  const Interaction * interaction, Option_t * opt) const
299 {
300 // Returns a multiplicity probability distribution for the input interaction.
301 // The input option (Default: "") can contain (combinations) of these strings:
302 // - "+LowMultSuppr": applies NeuGEN Rijk factors suppresing the low multipl.
303 // (1-pion and 2-pion) states as part of the DIS/RES joining scheme.
304 // - "+Renormalize": renormalizes the probability distribution after applying
305 // the NeuGEN scaling factors: Eg, when used as a hadronic multiplicity pdf
306 // the output hadronic multiplicity probability histogram needs to be re-
307 // normalized. But, when this method is called from a DIS cross section
308 // algorithm using the integrated probability reduction as a cross section
309 // section reduction factor then the output histogram should not be re-
310 // normalized after applying the scaling factors.
311 
312  if(!this->AssertValidity(interaction)) {
313  LOG("KNOHad", pWARN)
314  << "Returning a null multiplicity probability distribution!";
315  return 0;
316  }
317 
318  const InitialState & init_state = interaction->InitState();
319  int nu_pdg = init_state.ProbePdg();
320  int nuc_pdg = init_state.Tgt().HitNucPdg();
321 
322  // Compute the average charged hadron multiplicity as: <n> = a + b*ln(W^2)
323  // Calculate avergage hadron multiplicity (= 1.5 x charged hadron mult.)
324 
325  double W = utils::kinematics::W(interaction);
326  double avnch = this->AverageChMult(nu_pdg, nuc_pdg, W);
327  double avn = 1.5*avnch;
328 
329  SLOG("KNOHad", pINFO)
330  << "Average hadronic multiplicity (W=" << W << ") = " << avn;
331 
332  // Find the max possible multiplicity as W = Mneutron + (maxmult-1)*Mpion
333  double maxmult = this->MaxMult(interaction);
334 
335  // If required force the NeuGEN maximum multiplicity limit (10)
336  // Note: use for NEUGEN/GENIE comparisons, not physics MC production
337  if(fForceNeuGenLimit && maxmult>10) maxmult=10;
338 
339  // Set maximum multiplicity so that it does not exceed the max number of
340  // particles accepted by the ROOT phase space decayer (18)
341  // Change this if ROOT authors remove the TGenPhaseSpace limitation.
342  if(maxmult>18) maxmult=18;
343 
344  SLOG("KNOHad", pDEBUG) << "Computed maximum multiplicity = " << maxmult;
345 
346  if(maxmult<2) {
347  LOG("KNOHad", pWARN) << "Low maximum multiplicity! Quiting.";
348  return 0;
349  }
350 
351  // Create multiplicity probability histogram
352  TH1D * mult_prob = this->CreateMultProbHist(maxmult);
353 
354  // Compute the multiplicity probabilities values up to the bin corresponding
355  // to the computed maximum multiplicity
356 
357  if(maxmult>2) {
358  int nbins = mult_prob->FindBin(maxmult);
359 
360  for(int i = 1; i <= nbins; i++) {
361  // KNO distribution is <n>*P(n) vs n/<n>
362  double n = mult_prob->GetBinCenter(i); // bin centre
363  double z = n/avn; // z=n/<n>
364  double avnP = this->KNO(nu_pdg,nuc_pdg,z); // <n>*P(n)
365  double P = avnP / avn; // P(n)
366 
367  SLOG("KNOHad", pDEBUG)
368  << "n = " << n << " (n/<n> = " << z
369  << ", <n>*P = " << avnP << ") => P = " << P;
370 
371  mult_prob->Fill(n,P);
372  }
373  } else {
374  SLOG("KNOHad", pDEBUG) << "Fixing multiplicity to 2";
375  mult_prob->Fill(2,1.);
376  }
377 
378  double integral = mult_prob->Integral("width");
379  if(integral>0) {
380  // Normalize the probability distribution
381  mult_prob->Scale(1.0/integral);
382  } else {
383  SLOG("KNOHad", pWARN) << "probability distribution integral = 0";
384  return mult_prob;
385  }
386 
387  string option(opt);
388 
389  bool apply_neugen_Rijk = option.find("+LowMultSuppr") != string::npos;
390  bool renormalize = option.find("+Renormalize") != string::npos;
391 
392  // Apply the NeuGEN probability scaling factors -if requested-
393  if(apply_neugen_Rijk) {
394  SLOG("KNOHad", pINFO) << "Applying NeuGEN scaling factors";
395  // Only do so for W<Wcut
396  if(W<fWcut) {
397  this->ApplyRijk(interaction, renormalize, mult_prob);
398  } else {
399  SLOG("KNOHad", pDEBUG)
400  << "W = " << W << " < Wcut = " << fWcut
401  << " - Will not apply scaling factors";
402  }//<wcut?
403  }//apply?
404 
405  return mult_prob;
406 }
407 //____________________________________________________________________________
408 double KNOHadronization::Weight(void) const
409 {
410  return fWeight;
411 }
412 //____________________________________________________________________________
413 // methods overloading the default Algorithm interface implementation:
414 //____________________________________________________________________________
416 {
417  Algorithm::Configure(config);
418  this->LoadConfig();
419 }
420 //____________________________________________________________________________
422 {
423  Algorithm::Configure(config);
424  this->LoadConfig();
425 }
426 //____________________________________________________________________________
427 // private methods:
428 //____________________________________________________________________________
430 {
431  // Force decays of unstable hadronization products?
432  GetParamDef( "ForceDecays", fForceDecays, false ) ;
433 
434  // Force minimum multiplicity (if generated less than that) or abort?
435  GetParamDef( "ForceMinMultiplicity", fForceMinMult, true ) ;
436 
437  // Generate the baryon xF and pT^2 using experimental data as PDFs?
438  // In this case, only the N-1 other particles would be fed into the phase
439  // space decayer. This seems to improve hadronic system features such as
440  // bkw/fwd xF hemisphere average multiplicities.
441  // Note: not in the legacy KNO model (NeuGEN). Switch this feature off for
442  // comparisons or for reproducing old simulations.
443  GetParam( "KNO-UseBaryonPdfs-xFpT2", fUseBaryonXfPt2Param ) ;
444 
445  // Reweight the phase space decayer events to reproduce the experimentally
446  // measured pT^2 distributions.
447  // Note: not in the legacy KNO model (NeuGEN). Switch this feature off for
448  // comparisons or for reproducing old simulations.
449  GetParam( "KNO-PhaseSpDec-Reweight", fReWeightDecays ) ;
450 
451  // Parameter for phase space re-weighting. See ReWeightPt2()
452  GetParam( "KNO-PhaseSpDec-ReweightParm", fPhSpRwA ) ;
453 
454  // use isotropic non-reweighted 2-body phase space decays for consistency
455  // with neugen/daikon
456  GetParam( "KNO-UseIsotropic2BodyDec", fUseIsotropic2BDecays ) ;
457 
458  // Generated weighted or un-weighted hadronic systems?
459  GetParamDef( "GenerateWeighted", fGenerateWeighted, false ) ;
460 
461 
462  // Probabilities for producing hadron pairs
463 
464  //-- pi0 pi0
465  GetParam( "KNO-ProbPi0Pi0", fPpi0 ) ;
466  //-- pi+ pi-
467  GetParam( "KNO-ProbPiplusPiminus", fPpic ) ;
468  //-- K+ K-
469  GetParam( "KNO-ProbKplusKminus", fPKc ) ;
470  //-- K0 K0bar
471  GetParam( "KNO-ProbK0K0bar", fPK0 ) ;
472  //-- pi0 eta
473  GetParam( "KNO-ProbPi0Eta", fPpi0eta ) ;
474  //-- eta eta
475  GetParam( "KNO-ProbEtaEta", fPeta ) ;
476 
477  double fsum = fPeta + fPpi0eta + fPK0 + fPKc + fPpic + fPpi0;
478  double diff = TMath::Abs(1.-fsum);
479  if(diff>0.001) {
480  LOG("KNOHad", pWARN) << "KNO Probabilities do not sum to unity! Renormalizing..." ;
481  fPpi0 = fPpi0/fsum;
482  fPpic = fPpic/fsum;
483  fPKc = fPKc/fsum;
484  fPK0 = fPK0/fsum;
485  fPpi0eta = fPpi0eta/fsum;
486  fPeta = fPeta/fsum;
487  }
488 
489  // Decay unstable particles now or leave it for later? Which decayer to use?
490  fDecayer = 0;
491  if(fForceDecays) {
492  fDecayer = dynamic_cast<const DecayModelI *> (this->SubAlg("Decayer"));
493  assert(fDecayer);
494  }
495 
496  // Baryon pT^2 and xF parameterizations used as PDFs
497 
498  if (fBaryonXFpdf ) delete fBaryonXFpdf;
499  if (fBaryonPT2pdf) delete fBaryonPT2pdf;
500 
501  fBaryonXFpdf = new TF1("fBaryonXFpdf",
502  "0.083*exp(-0.5*pow(x+0.385,2.)/0.131)",-1,0.5);
503  fBaryonPT2pdf = new TF1("fBaryonPT2pdf",
504  "exp(-0.214-6.625*x)",0,0.6);
505  // stop ROOT from deleting these object of its own volition
506  gROOT->GetListOfFunctions()->Remove(fBaryonXFpdf);
507  gROOT->GetListOfFunctions()->Remove(fBaryonPT2pdf);
508 
509 /*
510  // load legacy KNO spline
511  fUseLegacyKNOSpline = fConfig->GetBoolDef("UseLegacyKNOSpl", false);
512 
513  // Delete the KNO spline from previous configuration of this instance
514  if(fKNO) delete fKNO;
515 
516  if(fUseLegacyKNOSpline) {
517  assert(gSystem->Getenv("GENIE"));
518  string basedir = gSystem->Getenv("GENIE");
519  string defknodata = basedir + "/data/kno/KNO.dat";
520  string knodata = fConfig->GetStringDef("kno-data",defknodata);
521 
522  LOG("KNOHad", pNOTICE) << "Loading KNO data from: " << knodata;
523  fKNO = new Spline(knodata);
524  }
525 */
526 
527  // Load parameters determining the average charged hadron multiplicity
528  GetParam( "KNO-Alpha-vp", fAvp ) ;
529  GetParam( "KNO-Alpha-vn", fAvn ) ;
530  GetParam( "KNO-Alpha-vbp", fAvbp ) ;
531  GetParam( "KNO-Alpha-vbn", fAvbn ) ;
532  GetParam( "KNO-Beta-vp", fBvp ) ;
533  GetParam( "KNO-Beta-vn", fBvn ) ;
534  GetParam( "KNO-Beta-vbp", fBvbp ) ;
535  GetParam( "KNO-Beta-vbn", fBvbn ) ;
536 
537  // Load parameters determining the prob of producing a strange baryon
538  // via associated production
539  GetParam( "KNO-Alpha-Hyperon", fAhyperon ) ;
540  GetParam( "KNO-Beta-Hyperon", fBhyperon ) ;
541 
542  // Load the Levy function parameter
543  GetParam( "KNO-LevyC-vp", fCvp ) ;
544  GetParam( "KNO-LevyC-vn", fCvn ) ;
545  GetParam( "KNO-LevyC-vbp", fCvbp ) ;
546  GetParam( "KNO-LevyC-vbn", fCvbn ) ;
547 
548  // Force NEUGEN upper limit in hadronic multiplicity (to be used only
549  // NEUGEN/GENIE comparisons)
550  GetParamDef( "ForceNeugenMultLimit", fForceNeuGenLimit, false ) ;
551 
552  // Load Wcut determining the phase space area where the multiplicity prob.
553  // scaling factors would be applied -if requested-
554  GetParam( "Wcut", fWcut ) ;
555 
556  // Load NEUGEN multiplicity probability scaling parameters Rijk
557  //neutrinos
558  GetParam( "DIS-HMultWgt-vp-CC-m2", fRvpCCm2 ) ;
559  GetParam( "DIS-HMultWgt-vp-CC-m3", fRvpCCm3 ) ;
560  GetParam( "DIS-HMultWgt-vp-NC-m2", fRvpNCm2 ) ;
561  GetParam( "DIS-HMultWgt-vp-NC-m3", fRvpNCm3 ) ;
562  GetParam( "DIS-HMultWgt-vn-CC-m2", fRvnCCm2 ) ;
563  GetParam( "DIS-HMultWgt-vn-CC-m3", fRvnCCm3 ) ;
564  GetParam( "DIS-HMultWgt-vn-NC-m2", fRvnNCm2 ) ;
565  GetParam( "DIS-HMultWgt-vn-NC-m3", fRvnNCm3 ) ;
566  //Anti-neutrinos
567  GetParam( "DIS-HMultWgt-vbp-CC-m2", fRvbpCCm2 ) ;
568  GetParam( "DIS-HMultWgt-vbp-CC-m3", fRvbpCCm3 ) ;
569  GetParam( "DIS-HMultWgt-vbp-NC-m2", fRvbpNCm2 ) ;
570  GetParam( "DIS-HMultWgt-vbp-NC-m3", fRvbpNCm3 ) ;
571  GetParam( "DIS-HMultWgt-vbn-CC-m2", fRvbnCCm2 ) ;
572  GetParam( "DIS-HMultWgt-vbn-CC-m3", fRvbnCCm3 ) ;
573  GetParam( "DIS-HMultWgt-vbn-NC-m2", fRvbnNCm2 ) ;
574  GetParam( "DIS-HMultWgt-vbn-NC-m3", fRvbnNCm3 ) ;
575 
576 }
577 //____________________________________________________________________________
578 double KNOHadronization::KNO(int probe_pdg, int nuc_pdg, double z) const
579 {
580 // Computes <n>P(n) for the input reduced multiplicity z=n/<n>
581 
582 /*
583  if(fUseLegacyKNOSpline) {
584  bool inrange = z > fKNO->XMin() && z < fKNO->XMax();
585  return (inrange) ? fKNO->Evaluate(z) : 0.;
586  }
587 */
588 
589  bool is_p = pdg::IsProton (nuc_pdg);
590  bool is_n = pdg::IsNeutron (nuc_pdg);
591  bool is_nu = pdg::IsNeutrino (probe_pdg);
592  bool is_nubar = pdg::IsAntiNeutrino (probe_pdg);
593  bool is_l = pdg::IsNegChargedLepton (probe_pdg);
594  bool is_lbar = pdg::IsPosChargedLepton (probe_pdg);
595  // EDIT
596  bool is_dm = pdg::IsDarkMatter (probe_pdg);
597 
598  double c=0; // Levy function parameter
599 
600  if ( is_p && (is_nu || is_l ) ) c=fCvp;
601  else if ( is_n && (is_nu || is_l ) ) c=fCvn;
602  else if ( is_p && (is_nubar || is_lbar) ) c=fCvbp;
603  else if ( is_n && (is_nubar || is_lbar) ) c=fCvbn;
604  // EDIT: assume it's neutrino-like for now...
605  else if ( is_p && is_dm ) c=fCvp;
606  else if ( is_n && is_dm ) c=fCvn;
607  else {
608  LOG("KNOHad", pERROR)
609  << "Invalid initial state (probe = " << probe_pdg << ", "
610  << "hit nucleon = " << nuc_pdg << ")";
611  return 0;
612  }
613 
614  double x = c*z+1;
615  double kno = 2*TMath::Exp(-c)*TMath::Power(c,x)/TMath::Gamma(x);
616 
617  return kno;
618 }
619 //____________________________________________________________________________
621  int probe_pdg,int nuc_pdg, double W) const
622 {
623 // computes the average charged multiplicity
624 //
625  bool is_p = pdg::IsProton (nuc_pdg);
626  bool is_n = pdg::IsNeutron (nuc_pdg);
627  bool is_nu = pdg::IsNeutrino (probe_pdg);
628  bool is_nubar = pdg::IsAntiNeutrino (probe_pdg);
629  bool is_l = pdg::IsNegChargedLepton (probe_pdg);
630  bool is_lbar = pdg::IsPosChargedLepton (probe_pdg);
631  // EDIT
632  bool is_dm = pdg::IsDarkMatter (probe_pdg);
633 
634  double a=0, b=0; // params controlling average multiplicity
635 
636  if ( is_p && (is_nu || is_l ) ) { a=fAvp; b=fBvp; }
637  else if ( is_n && (is_nu || is_l ) ) { a=fAvn; b=fBvn; }
638  else if ( is_p && (is_nubar || is_lbar) ) { a=fAvbp; b=fBvbp; }
639  else if ( is_n && (is_nubar || is_lbar) ) { a=fAvbn; b=fBvbn; }
640  // EDIT: assume it's neutrino-like for now...
641  else if ( is_p && is_dm ) { a=fAvp; b=fBvp; }
642  else if ( is_n && is_dm ) { a=fAvn; b=fBvn; }
643  else {
644  LOG("KNOHad", pERROR)
645  << "Invalid initial state (probe = " << probe_pdg << ", "
646  << "hit nucleon = " << nuc_pdg << ")";
647  return 0;
648  }
649 
650  double av_nch = a + b * 2*TMath::Log(W);
651  return av_nch;
652 }
653 //____________________________________________________________________________
655 {
656 // Returns the hadron shower charge in units of +e
657 // HadronShowerCharge = Q{initial} - Q{final state primary lepton}
658 // eg in v p -> l- X the hadron shower charge is +2
659 // in v n -> l- X the hadron shower charge is +1
660 // in v n -> v X the hadron shower charge is 0
661 //
662  int hadronShowerCharge = 0;
663 
664  // find out the charge of the final state lepton
665  double ql = interaction->FSPrimLepton()->Charge() / 3.;
666 
667  // find out the charge of the probe
668  double qp = interaction->InitState().Probe()->Charge() / 3.;
669 
670  // get the initial state, ask for the hit-nucleon and get
671  // its charge ( = initial state charge for vN interactions)
672  const InitialState & init_state = interaction->InitState();
673  int hit_nucleon = init_state.Tgt().HitNucPdg();
674 
675  assert( pdg::IsProton(hit_nucleon) || pdg::IsNeutron(hit_nucleon) );
676 
677  // Ask PDGLibrary for the nucleon charge
678  double qnuc = PDGLibrary::Instance()->Find(hit_nucleon)->Charge() / 3.;
679 
680  // calculate the hadron shower charge
681  hadronShowerCharge = (int) ( qp + qnuc - ql );
682 
683  return hadronShowerCharge;
684 }
685 //____________________________________________________________________________
687  double W, const PDGCodeList & pdgv, bool reweight_decays) const
688 {
689 // Simple phase space decay including all generated particles.
690 // The old NeuGEN decay strategy.
691 
692  LOG("KNOHad", pINFO) << "** Using Hadronic System Decay method 1";
693 
694  TLorentzVector p4had(0,0,0,W);
695  TClonesArray * plist = new TClonesArray("TMCParticle", pdgv.size());
696 
697  // do the decay
698  bool ok = this->PhaseSpaceDecay(*plist, p4had, pdgv, 0, reweight_decays);
699 
700  // clean-up and return NULL
701  if(!ok) {
702  plist->Delete();
703  delete plist;
704  return 0;
705  }
706  return plist;
707 }
708 //____________________________________________________________________________
710  double W, const PDGCodeList & pdgv, bool reweight_decays) const
711 {
712 // Generate the baryon based on experimental pT^2 and xF distributions
713 // Then pass the remaining system of N-1 particles to a phase space decayer.
714 // The strategy adopted at the July-2006 hadronization model mini-workshop.
715 
716  LOG("KNOHad", pINFO) << "** Using Hadronic System Decay method 2";
717 
718  // If only 2 particles are input then don't call the phase space decayer
719  if(pdgv.size() == 2) return this->DecayBackToBack(W,pdgv);
720 
721  // Now handle the more general case:
722 
723  // Take the baryon
724  int baryon = pdgv[0];
725  double MN = PDGLibrary::Instance()->Find(baryon)->Mass();
726  double MN2 = TMath::Power(MN, 2);
727 
728  // Check baryon code
729  // ...
730 
731  // Strip the PDG list from the baryon
732  bool allowdup = true;
733  PDGCodeList pdgv_strip(pdgv.size()-1, allowdup);
734  for(unsigned int i=1; i<pdgv.size(); i++) pdgv_strip[i-1] = pdgv[i];
735 
736  // Get the sum of all masses for the particles in the stripped list
737  double mass_sum = 0;
738  vector<int>::const_iterator pdg_iter = pdgv_strip.begin();
739  for( ; pdg_iter != pdgv_strip.end(); ++pdg_iter) {
740  int pdgc = *pdg_iter;
741  mass_sum += PDGLibrary::Instance()->Find(pdgc)->Mass();
742  }
743 
744  // Create the particle list
745  TClonesArray * plist = new TClonesArray("TMCParticle", pdgv.size());
746 
748  TLorentzVector p4had(0,0,0,W);
749  TLorentzVector p4N (0,0,0,0);
750  TLorentzVector p4d;
751 
752  // generate the N 4-p independently
753 
754  bool got_baryon_4p = false;
755  bool got_hadsyst_4p = false;
756 
757  while(!got_hadsyst_4p) {
758 
759  LOG("KNOHad", pINFO) << "Generating p4 for baryon with pdg = " << baryon;
760 
761  while(!got_baryon_4p) {
762 
763  //-- generate baryon xF and pT2
764  double xf = fBaryonXFpdf ->GetRandom();
765  double pt2 = fBaryonPT2pdf->GetRandom();
766 
767  //-- generate baryon px,py,pz
768  double pt = TMath::Sqrt(pt2);
769  double phi = (2*kPi) * rnd->RndHadro().Rndm();
770  double px = pt * TMath::Cos(phi);
771  double py = pt * TMath::Sin(phi);
772  double pz = xf*W/2;
773  double p2 = TMath::Power(pz,2) + pt2;
774  double E = TMath::Sqrt(p2+MN2);
775 
776  p4N.SetPxPyPzE(px,py,pz,E);
777 
778  LOG("KNOHad", pDEBUG) << "Trying nucleon xF= "<< xf<< ", pT2= "<< pt2;
779 
780  //-- check whether there is phase space for the remnant N-1 system
781  p4d = p4had-p4N; // 4-momentum vector for phase space decayer
782  double Mav = p4d.Mag();
783 
784  got_baryon_4p = (Mav > mass_sum);
785 
786  } // baryon xf,pt2 seletion
787 
788  LOG("KNOHad", pINFO)
789  << "Generated baryon with P4 = " << utils::print::P4AsString(&p4N);
790 
791  // Insert the baryon at the event record
792  new ((*plist)[0]) TMCParticle(
793  1,baryon,-1,-1,-1, p4N.Px(),p4N.Py(),p4N.Pz(),p4N.Energy(),MN, 0,0,0,0,0);
794 
795  // Do a phase space decay for the N-1 particles and add them to the list
796  LOG("KNOHad", pINFO)
797  << "Generating p4 for the remaining hadronic system";
798  LOG("KNOHad", pINFO)
799  << "Remaining system: Available mass = " << p4d.Mag()
800  << ", Particle masses = " << mass_sum;
801 
802  bool is_ok = this->PhaseSpaceDecay(
803  *plist, p4d, pdgv_strip, 1, reweight_decays);
804 
805  got_hadsyst_4p = is_ok;
806 
807  if(!got_hadsyst_4p) {
808  got_baryon_4p = false;
809  plist->Delete();
810  }
811  }
812 
813  // clean-up and return NULL
814  if(0) {
815  LOG("KNOHad", pERROR) << "*** Decay forbidden by kinematics! ***";
816  plist->Delete();
817  delete plist;
818  return 0;
819  }
820  return plist;
821 }
822 //____________________________________________________________________________
824  double W, const PDGCodeList & pdgv) const
825 {
826 // Handles a special case (only two particles) of the 2nd decay method
827 //
828 
829  LOG("KNOHad", pINFO) << "Generating two particles back-to-back";
830 
831  assert(pdgv.size()==2);
832 
834 
835  // Create the particle list
836  TClonesArray * plist = new TClonesArray("TMCParticle", pdgv.size());
837 
838  // Get xF,pT2 distribution (y-) maxima for the rejection method
839  double xFo = 1.1 * fBaryonXFpdf ->GetMaximum(-1,1);
840  double pT2o = 1.1 * fBaryonPT2pdf->GetMaximum( 0,1);
841 
842  TLorentzVector p4(0,0,0,W); // 2-body hadronic system 4p
843 
844  // Do the 2-body decay
845  bool accepted = false;
846  while(!accepted) {
847 
848  // Find an allowed (unweighted) phase space decay for the 2 particles
849  // and add them to the list
850  bool ok = this->PhaseSpaceDecay(*plist, p4, pdgv, 0, false);
851 
852  // If the decay isn't allowed clean-up and return NULL
853  if(!ok) {
854  LOG("KNOHad", pERROR) << "*** Decay forbidden by kinematics! ***";
855  plist->Delete();
856  delete plist;
857  return 0;
858  }
859 
860  // If the decay was allowed, then compute the baryon xF,pT2 and accept/
861  // reject the phase space decays so as to reproduce the xF,pT2 PDFs
862 
863  TMCParticle * baryon = (TMCParticle *) (*plist)[0];
864  assert(pdg::IsNeutronOrProton(baryon->GetKF()));
865 
866  double px = baryon->GetPx();
867  double py = baryon->GetPy();
868  double pz = baryon->GetPz();
869 
870  double pT2 = px*px + py*py;
871  double pL = pz;
872  double xF = pL/(W/2);
873 
874  double pT2rnd = pT2o * rnd->RndHadro().Rndm();
875  double xFrnd = xFo * rnd->RndHadro().Rndm();
876 
877  double pT2pdf = fBaryonPT2pdf->Eval(pT2);
878  double xFpdf = fBaryonXFpdf ->Eval(xF );
879 
880  LOG("KNOHad", pINFO) << "baryon xF = " << xF << ", pT2 = " << pT2;
881 
882  accepted = (xFrnd < xFpdf && pT2rnd < pT2pdf);
883 
884  LOG("KNOHad", pINFO) << ((accepted) ? "Decay accepted":"Decay rejected");
885  }
886  return plist;
887 }
888 //____________________________________________________________________________
890  TClonesArray & plist, TLorentzVector & pd,
891  const PDGCodeList & pdgv, int offset, bool reweight) const
892 {
893 // General method decaying the input particle system 'pdgv' with available 4-p
894 // given by 'pd'. The decayed system is used to populate the input TMCParticle
895 // array starting from the slot 'offset'.
896 //
897  LOG("KNOHad", pINFO) << "*** Performing a Phase Space Decay";
898  LOG("KNOHad", pINFO) << "pT reweighting is " << (reweight ? "on" : "off");
899 
900  assert ( offset >= 0);
901  assert ( pdgv.size() > 1);
902 
903  // Get the decay product masses
904 
905  vector<int>::const_iterator pdg_iter;
906  int i = 0;
907  double * mass = new double[pdgv.size()];
908  double sum = 0;
909  for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
910  int pdgc = *pdg_iter;
911  double m = PDGLibrary::Instance()->Find(pdgc)->Mass();
912  mass[i++] = m;
913  sum += m;
914  }
915 
916  LOG("KNOHad", pINFO)
917  << "Decaying N = " << pdgv.size() << " particles / total mass = " << sum;
918  LOG("KNOHad", pINFO)
919  << "Decaying system p4 = " << utils::print::P4AsString(&pd);
920 
921  // Set the decay
922  bool permitted = fPhaseSpaceGenerator.SetDecay(pd, pdgv.size(), mass);
923  if(!permitted) {
924  LOG("KNOHad", pERROR)
925  << " *** Phase space decay is not permitted \n"
926  << " Total particle mass = " << sum << "\n"
927  << " Decaying system p4 = " << utils::print::P4AsString(&pd);
928 
929  // clean-up and return
930  delete [] mass;
931  return false;
932  }
933 
934  // Get the maximum weight
935  //double wmax = fPhaseSpaceGenerator.GetWtMax();
936  double wmax = -1;
937  for(int idec=0; idec<200; idec++) {
938  double w = fPhaseSpaceGenerator.Generate();
939  if(reweight) { w *= this->ReWeightPt2(pdgv); }
940  wmax = TMath::Max(wmax,w);
941  }
942  assert(wmax>0);
943 
944  LOG("KNOHad", pNOTICE)
945  << "Max phase space gen. weight @ current hadronic system: " << wmax;
946 
947  // Generate a weighted or unweighted decay
948 
950 
951  if(fGenerateWeighted)
952  {
953  // *** generating weighted decays ***
954  double w = fPhaseSpaceGenerator.Generate();
955  if(reweight) { w *= this->ReWeightPt2(pdgv); }
956  fWeight *= TMath::Max(w/wmax, 1.);
957  }
958  else
959  {
960  // *** generating un-weighted decays ***
961  wmax *= 2.3;
962  bool accept_decay=false;
963  unsigned int itry=0;
964 
965  while(!accept_decay)
966  {
967  itry++;
968 
969  if(itry>kMaxUnweightDecayIterations) {
970  // report, clean-up and return
971  LOG("KNOHad", pWARN)
972  << "Couldn't generate an unweighted phase space decay after "
973  << itry << " attempts";
974  delete [] mass;
975  return false;
976  }
977 
978  double w = fPhaseSpaceGenerator.Generate();
979  if(reweight) { w *= this->ReWeightPt2(pdgv); }
980  if(w > wmax) {
981  LOG("KNOHad", pWARN)
982  << "Decay weight = " << w << " > max decay weight = " << wmax;
983  }
984  double gw = wmax * rnd->RndHadro().Rndm();
985  accept_decay = (gw<=w);
986 
987  LOG("KNOHad", pINFO)
988  << "Decay weight = " << w << " / R = " << gw
989  << " - accepted: " << accept_decay;
990 
991  bool return_after_not_accepted_decay = false;
992  if(return_after_not_accepted_decay && !accept_decay) {
993  LOG("KNOHad", pWARN)
994  << "Was instructed to return after a not-accepted decay";
995  delete [] mass;
996  return false;
997  }
998  }
999  }
1000 
1001  // Insert final state products into a TClonesArray of TMCParticles
1002 
1003  i=0;
1004  for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
1005 
1006  //-- current PDG code
1007  int pdgc = *pdg_iter;
1008 
1009  //-- get the 4-momentum of the i-th final state particle
1010  TLorentzVector * p4fin = fPhaseSpaceGenerator.GetDecay(i);
1011 
1012  new ( plist[offset+i] ) TMCParticle(
1013  1, /* KS Code */
1014  pdgc, /* PDG Code */
1015  -1, /* parent particle */
1016  -1, /* first child particle */
1017  -1, /* last child particle */
1018  p4fin->Px(), /* 4-momentum: px component */
1019  p4fin->Py(), /* 4-momentum: py component */
1020  p4fin->Pz(), /* 4-momentum: pz component */
1021  p4fin->Energy(), /* 4-momentum: E component */
1022  mass[i], /* particle mass */
1023  0, /* production vertex 4-vector: vx */
1024  0, /* production vertex 4-vector: vy */
1025  0, /* production vertex 4-vector: vz */
1026  0, /* production vertex 4-vector: time */
1027  0 /* lifetime */
1028  );
1029  i++;
1030  }
1031 
1032  // Clean-up
1033  delete [] mass;
1034 
1035  return true;
1036 }
1037 //____________________________________________________________________________
1038 double KNOHadronization::ReWeightPt2(const PDGCodeList & pdgcv) const
1039 {
1040 // Phase Space Decay re-weighting to reproduce exp(-pT2/<pT2>) pion pT2
1041 // distributions.
1042 // See: A.B.Clegg, A.Donnachie, A Description of Jet Structure by pT-limited
1043 // Phase Space.
1044 
1045  double w = 1;
1046 
1047  for(unsigned int i = 0; i < pdgcv.size(); i++) {
1048 
1049  //int pdgc = pdgcv[i];
1050  //if(pdgc!=kPdgPiP&&pdgc!=kPdgPiM) continue;
1051 
1052  TLorentzVector * p4 = fPhaseSpaceGenerator.GetDecay(i);
1053  double pt2 = TMath::Power(p4->Px(),2) + TMath::Power(p4->Py(),2);
1054  double wi = TMath::Exp(-fPhSpRwA*TMath::Sqrt(pt2));
1055  //double wi = (9.41 * TMath::Landau(pt2,0.24,0.12));
1056 
1057  w *= wi;
1058  }
1059  return w;
1060 }
1061 //____________________________________________________________________________
1063  int multiplicity, int maxQ, double W) const
1064 {
1065 // Selection of fragments (identical as in NeuGEN).
1066 
1067  // Get PDG library and rnd num generator
1070 
1071  // Create vector to add final state hadron PDG codes
1072  bool allowdup=true;
1073  PDGCodeList * pdgc = new PDGCodeList(allowdup);
1074  //pdgc->reserve(multiplicity);
1075  int hadrons_to_add = multiplicity;
1076 
1077  //
1078  // Assign baryon as p, n, Sigma+, Sigma- or Lambda
1079  //
1080 
1081  int baryon_code = this->GenerateBaryonPdgCode(multiplicity, maxQ, W);
1082  pdgc->push_back(baryon_code);
1083 
1084  bool baryon_is_strange = (baryon_code == kPdgSigmaP ||
1085  baryon_code == kPdgLambda ||
1086  baryon_code == kPdgSigmaM);
1087  bool baryon_chg_is_pos = (baryon_code == kPdgProton ||
1088  baryon_code == kPdgSigmaP);
1089  bool baryon_chg_is_neg = (baryon_code == kPdgSigmaM);
1090 
1091  // Update number of hadrons to add, available shower charge & invariant mass
1092  if(baryon_chg_is_pos) maxQ -= 1;
1093  if(baryon_chg_is_neg) maxQ += 1;
1094  hadrons_to_add--;
1095  W -= pdg->Find( (*pdgc)[0] )->Mass();
1096 
1097  //
1098  // Assign remaining hadrons up to n = multiplicity
1099  //
1100 
1101  // Conserve strangeness
1102  if(baryon_is_strange) {
1103  LOG("KNOHad", pDEBUG)
1104  << " Remnant baryon is strange. Conserving strangeness...";
1105 
1106  //conserve strangeness and handle charge imbalance with one particle
1107  if(multiplicity == 2) {
1108  if(maxQ == 1) {
1109  LOG("KNOHad", pDEBUG) << " -> Adding a K+";
1110  pdgc->push_back( kPdgKP );
1111 
1112  // update n-of-hadrons to add, avail. shower charge & invariant mass
1113  maxQ -= 1;
1114  hadrons_to_add--;
1115  W -= pdg->Find(kPdgKP)->Mass();
1116  }
1117  else if(maxQ == 0) {
1118  LOG("KNOHad", pDEBUG) << " -> Adding a K0";
1119  pdgc->push_back( kPdgK0 );
1120 
1121  // update n-of-hadrons to add, avail. shower charge & invariant mass
1122  hadrons_to_add--;
1123  W -= pdg->Find(kPdgK0)->Mass();
1124  }
1125  }
1126 
1127  //only two particles left to balance charge
1128  else if(multiplicity == 3 && maxQ == 2) {
1129  LOG("KNOHad", pDEBUG) << " -> Adding a K+";
1130  pdgc->push_back( kPdgKP );
1131 
1132  // update n-of-hadrons to add, avail. shower charge & invariant mass
1133  maxQ -= 1;
1134  hadrons_to_add--;
1135  W -= pdg->Find(kPdgKP)->Mass();
1136  }
1137  else if(multiplicity == 3 && maxQ == -1) { //adding K+ makes it impossible to balance charge
1138  LOG("KNOHad", pDEBUG) << " -> Adding a K0";
1139  pdgc->push_back( kPdgK0 );
1140 
1141  // update n-of-hadrons to add, avail. shower charge & invariant mass
1142  hadrons_to_add--;
1143  W -= pdg->Find(kPdgK0)->Mass();
1144  }
1145 
1146  //simply conserve strangeness, without regard to charge
1147  else {
1148  double y = rnd->RndHadro().Rndm();
1149  if(y < 0.5) {
1150  LOG("KNOHad", pDEBUG) <<" -> Adding a K+";
1151  pdgc->push_back( kPdgKP );
1152 
1153  // update n-of-hadrons to add, avail. shower charge & invariant mass
1154  maxQ -= 1;
1155  hadrons_to_add--;
1156  W -= pdg->Find(kPdgKP)->Mass();
1157  }
1158  else {
1159  LOG("KNOHad", pDEBUG) <<" -> Adding a K0";
1160  pdgc->push_back( kPdgK0 );
1161 
1162  // update n-of-hadrons to add, avail. shower charge & invariant mass
1163  hadrons_to_add--;
1164  W -= pdg->Find(kPdgK0)->Mass();
1165  }
1166  }
1167  }//if the baryon is strange
1168 
1169  // Handle charge imbalance
1170  while(maxQ != 0) {
1171 
1172  if (maxQ < 0) {
1173  // Need more negative charge
1174  LOG("KNOHad", pDEBUG) << "Need more negative charge -> Adding a pi-";
1175  pdgc->push_back( kPdgPiM );
1176 
1177  // update n-of-hadrons to add, avail. shower charge & invariant mass
1178  maxQ += 1;
1179  hadrons_to_add--;
1180 
1181  W -= pdg->Find(kPdgPiM)->Mass();
1182 
1183  } else if (maxQ > 0) {
1184  // Need more positive charge
1185  LOG("KNOHad", pDEBUG) << "Need more positive charge -> Adding a pi+";
1186  pdgc->push_back( kPdgPiP );
1187 
1188  // update n-of-hadrons to add, avail. shower charge & invariant mass
1189  maxQ -= 1;
1190  hadrons_to_add--;
1191 
1192  W -= pdg->Find(kPdgPiP)->Mass();
1193  }
1194  }
1195 
1196  // Add remaining neutrals or pairs up to the generated multiplicity
1197  if(maxQ == 0) {
1198 
1199  LOG("KNOHad", pDEBUG)
1200  << "Hadronic charge balanced. Now adding only neutrals or +- pairs";
1201 
1202  // Final state has correct charge.
1203  // Now add pi0 or pairs (pi0 pi0 / pi+ pi- / K+ K- / K0 K0bar) only
1204 
1205  // Masses of particle pairs
1206  double M2pi0 = 2 * pdg -> Find (kPdgPi0) -> Mass();
1207  double M2pic = pdg -> Find (kPdgPiP) -> Mass() +
1208  pdg -> Find (kPdgPiM) -> Mass();
1209  double M2Kc = pdg -> Find (kPdgKP ) -> Mass() +
1210  pdg -> Find (kPdgKM ) -> Mass();
1211  double M2K0 = 2 * pdg -> Find (kPdgK0 ) -> Mass();
1212  double M2Eta = 2 * pdg -> Find (kPdgEta) -> Mass();
1213  double Mpi0eta = pdg -> Find (kPdgPi0) -> Mass() +
1214  pdg -> Find (kPdgEta) -> Mass();
1215 
1216  // Prevent multiplicity overflow.
1217  // Check if we have an odd number of hadrons to add.
1218  // If yes, add a single pi0 and then go on and add pairs
1219 
1220  if( hadrons_to_add > 0 && hadrons_to_add % 2 == 1 ) {
1221 
1222  LOG("KNOHad", pDEBUG)
1223  << "Odd number of hadrons left to add -> Adding a pi0";
1224  pdgc->push_back( kPdgPi0 );
1225 
1226  // update n-of-hadrons to add & available invariant mass
1227  hadrons_to_add--;
1228  W -= pdg->Find(kPdgPi0)->Mass();
1229  }
1230 
1231  // Now add pairs (pi0 pi0 / pi+ pi- / K+ K- / K0 K0bar)
1232  assert( hadrons_to_add % 2 == 0 ); // even number
1233  LOG("KNOHad", pDEBUG)
1234  <<" hadrons_to_add = "<<hadrons_to_add<<" W= "<<W<<" M2pi0 = "<<M2pi0<<" M2pic = "<<M2pic<<" M2Kc = "<<M2Kc<<" M2K0= "<<M2K0<<" M2Eta= "<<M2Eta;
1235 
1236  while(hadrons_to_add > 0 && W >= M2pi0) {
1237 
1238  double x = rnd->RndHadro().Rndm();
1239  LOG("KNOHad", pDEBUG) << "rndm = " << x;
1240  // Add a pi0 pair
1241  if (x >= 0 && x < fPpi0) {
1242  LOG("KNOHad", pDEBUG) << " -> Adding a pi0pi0 pair";
1243  pdgc->push_back( kPdgPi0 );
1244  pdgc->push_back( kPdgPi0 );
1245  hadrons_to_add -= 2; // update the number of hadrons to add
1246  W -= M2pi0; // update the available invariant mass
1247  }
1248 
1249  // Add a pi+ pi- pair
1250  else if (x < fPpi0 + fPpic) {
1251  if(W >= M2pic) {
1252  LOG("KNOHad", pDEBUG) << " -> Adding a pi+pi- pair";
1253  pdgc->push_back( kPdgPiP );
1254  pdgc->push_back( kPdgPiM );
1255  hadrons_to_add -= 2; // update the number of hadrons to add
1256  W -= M2pic; // update the available invariant mass
1257  } else {
1258  LOG("KNOHad", pDEBUG)
1259  << "Not enough mass for a pi+pi-: trying something else";
1260  }
1261  }
1262 
1263  // Add a K+ K- pair
1264  else if (x < fPpi0 + fPpic + fPKc) {
1265  if(W >= M2Kc) {
1266  LOG("KNOHad", pDEBUG) << " -> Adding a K+K- pair";
1267  pdgc->push_back( kPdgKP );
1268  pdgc->push_back( kPdgKM );
1269  hadrons_to_add -= 2; // update the number of hadrons to add
1270  W -= M2Kc; // update the available invariant mass
1271  } else {
1272  LOG("KNOHad", pDEBUG)
1273  << "Not enough mass for a K+K-: trying something else";
1274  }
1275  }
1276 
1277  // Add a K0 - \bar{K0} pair
1278  else if (x <= fPpi0 + fPpic + fPKc + fPK0) {
1279  if( W >= M2K0 ) {
1280  LOG("KNOHad", pDEBUG) << " -> Adding a K0 K0bar pair";
1281  pdgc->push_back( kPdgK0 );
1282  pdgc->push_back( kPdgAntiK0 );
1283  hadrons_to_add -= 2; // update the number of hadrons to add
1284  W -= M2K0; // update the available invariant mass
1285  } else {
1286  LOG("KNOHad", pDEBUG)
1287  << "Not enough mass for a K0 K0bar: trying something else";
1288  }
1289  }
1290 
1291  // Add a Pi0-Eta pair
1292  else if (x <= fPpi0 + fPpic + fPKc + fPK0 + fPpi0eta) {
1293  if( W >= Mpi0eta ) {
1294  LOG("KNOHad", pDEBUG) << " -> Adding a Pi0-Eta pair";
1295  pdgc->push_back( kPdgPi0 );
1296  pdgc->push_back( kPdgEta );
1297  hadrons_to_add -= 2; // update the number of hadrons to add
1298  W -= Mpi0eta; // update the available invariant mass
1299  } else {
1300  LOG("KNOHad", pDEBUG)
1301  << "Not enough mass for a Pi0-Eta pair: trying something else";
1302  }
1303  }
1304 
1305  //Add a Eta pair
1306  else if(x <= fPpi0 + fPpic + fPKc + fPK0 + fPpi0eta + fPeta) {
1307  if( W >= M2Eta ){
1308  LOG("KNOHad", pDEBUG) << " -> Adding a eta-eta pair";
1309  pdgc->push_back( kPdgEta );
1310  pdgc->push_back( kPdgEta );
1311  hadrons_to_add -= 2; // update the number of hadrons to add
1312  W -= M2Eta; // update the available invariant mass
1313  } else {
1314  LOG("KNOHad", pDEBUG)
1315  << "Not enough mass for a Eta-Eta pair: trying something else";
1316  }
1317 
1318  } else {
1319  LOG("KNOHad", pERROR)
1320  << "Hadron Assignment Probabilities do not add up to 1!!";
1321  exit(1);
1322  }
1323 
1324  // make sure it has enough invariant mass to reach the
1325  // given multiplicity, even by adding only the lightest
1326  // hadron pairs (pi0's)
1327  // Otherwise force a lower multiplicity.
1328  if(W < M2pi0) hadrons_to_add = 0;
1329 
1330  } // while there are more hadrons to add
1331  } // if charge is balanced (maxQ == 0)
1332 
1333  return pdgc;
1334 }
1335 //____________________________________________________________________________
1337  int multiplicity, int maxQ, double W) const
1338 {
1339 // Selection of main target fragment (identical as in NeuGEN).
1340 // Assign baryon as p or n. Force it for ++ and - I=3/2 at mult. = 2
1341 
1343  double x = rnd->RndHadro().Rndm();
1344  double y = rnd->RndHadro().Rndm();
1345 
1346  // initialize to neutron & then change it to proton if you must
1347  int pdgc = kPdgNeutron;
1348 
1349  // Assign a probability for the given W for the baryon to become strange
1350  // using a function derived from a fit to the data in Jones et al. (1993)
1351  // Don't let the probability be larger than 1.
1352  double Pstr = fAhyperon + fBhyperon * TMath::Log(W*W);
1353  Pstr = TMath::Min(1.,Pstr);
1354  Pstr = TMath::Max(0.,Pstr);
1355 
1356  // Available hadronic system charge = 2
1357  if(maxQ == 2) {
1358  //for multiplicity ==2, force it to p
1359  if(multiplicity ==2 ) pdgc = kPdgProton;
1360  else {
1361  if(x < 0.66667) pdgc = kPdgProton;
1362  }
1363  }
1364  // Available hadronic system charge = 1
1365  if(maxQ == 1) {
1366  if(multiplicity == 2) {
1367  if(x < 0.33333) pdgc = kPdgProton;
1368  } else {
1369  if(x < 0.50000) pdgc = kPdgProton;
1370  }
1371  }
1372 
1373  // Available hadronic system charge = 0
1374  if(maxQ == 0) {
1375  if(multiplicity == 2) {
1376  if(x < 0.66667) pdgc = kPdgProton;
1377  } else {
1378  if(x < 0.50000) pdgc = kPdgProton;
1379  }
1380  }
1381  // Available hadronic system charge = -1
1382  if(maxQ == -1) {
1383  // for multiplicity == 2, force it to n
1384  if(multiplicity != 2) {
1385  if(x < 0.33333) pdgc = kPdgProton;
1386  }
1387  }
1388 
1389  // For neutrino interactions turn protons and neutrons to Sigma+ and
1390  // Lambda respectively (Lambda and Sigma- respectively for anti-neutrino
1391  // interactions).
1392  if(pdgc == kPdgProton && y < Pstr && maxQ > 0) {
1393  pdgc = kPdgSigmaP;
1394  }
1395  else if(pdgc == kPdgProton && y < Pstr && maxQ <= 0) {
1396  pdgc = kPdgLambda;
1397  }
1398  else if(pdgc == kPdgNeutron && y < Pstr && maxQ > 0) {
1399  pdgc = kPdgLambda;
1400  }
1401  else if(pdgc == kPdgNeutron && y < Pstr && maxQ <= 0) {
1402  pdgc = kPdgSigmaM;
1403  }
1404 
1405  if(pdgc == kPdgProton)
1406  LOG("KNOHad", pDEBUG) << " -> Adding a proton";
1407  if(pdgc == kPdgNeutron)
1408  LOG("KNOHad", pDEBUG) << " -> Adding a neutron";
1409  if(pdgc == kPdgSigmaP)
1410  LOG("KNOHad", pDEBUG) << " -> Adding a sigma+";
1411  if(pdgc == kPdgLambda)
1412  LOG("KNOHad", pDEBUG) << " -> Adding a lambda";
1413  if(pdgc == kPdgSigmaM)
1414  LOG("KNOHad", pDEBUG) << " -> Adding a sigma-";
1415 
1416  return pdgc;
1417 }
1418 //____________________________________________________________________________
1419 void KNOHadronization::HandleDecays(TClonesArray * plist) const
1420 {
1421 // Handle decays of unstable particles if requested through the XML config.
1422 // The default is not to decay the particles at this stage (during event
1423 // generation, the UnstableParticleDecayer event record visitor decays what
1424 // is needed to be decayed later on). But, when comparing various models
1425 // (eg PYTHIA vs KNO) independently and not within the full MC simulation
1426 // framework it might be necessary to force the decays at this point.
1427 
1428  if (fForceDecays) {
1429  assert(fDecayer);
1430 
1431  //-- loop through the fragmentation event record & decay unstables
1432  int idecaying = -1; // position of decaying particle
1433  TMCParticle * p = 0; // current particle
1434 
1435  TIter piter(plist);
1436  while ( (p = (TMCParticle *) piter.Next()) ) {
1437 
1438  idecaying++;
1439  int status = p->GetKS();
1440 
1441  // bother for final state particle only
1442  if(status < 10) {
1443 
1444  // until ROOT's T(MC)Particle(PDG) Lifetime() is fixed, decay only
1445  // pi^0's
1446  if ( p->GetKF() == kPdgPi0 ) {
1447 
1448  DecayerInputs_t dinp;
1449 
1450  TLorentzVector p4;
1451  p4.SetPxPyPzE(p->GetPx(), p->GetPy(), p->GetPz(), p->GetEnergy());
1452 
1453  dinp.PdgCode = p->GetKF();
1454  dinp.P4 = &p4;
1455 
1456  TClonesArray * decay_products = fDecayer->Decay(dinp);
1457 
1458  if(decay_products) {
1459  //-- mark the parent particle as decayed & set daughters
1460  p->SetKS(11);
1461 
1462  int nfp = plist->GetEntries(); // n. fragm. products
1463  int ndp = decay_products->GetEntries(); // n. decay products
1464 
1465  p->SetFirstChild ( nfp ); // decay products added at
1466  p->SetLastChild ( nfp + ndp -1 ); // the end of the fragm.rec.
1467 
1468  //-- add decay products to the fragmentation record
1469  TMCParticle * dp = 0;
1470  TIter dpiter(decay_products);
1471 
1472  while ( (dp = (TMCParticle *) dpiter.Next()) ) {
1473 
1474  dp->SetParent(idecaying);
1475  new ( (*plist)[plist->GetEntries()] ) TMCParticle(*dp);
1476  }
1477 
1478  //-- clean up decay products
1479  decay_products->Delete();
1480  delete decay_products;
1481  }
1482 
1483  } // particle is to be decayed
1484  } // KS < 10 : final state particle (as in PYTHIA LUJETS record)
1485  } // particles in fragmentation record
1486  } // force decay
1487 }
1488 //____________________________________________________________________________
1489 bool KNOHadronization::AssertValidity(const Interaction * interaction) const
1490 {
1491  if(interaction->ExclTag().IsCharmEvent()) {
1492  LOG("KNOHad", pWARN) << "Can't hadronize charm events";
1493  return false;
1494  }
1495  double W = utils::kinematics::W(interaction);
1496  if(W < this->Wmin()) {
1497  LOG("KNOHad", pWARN) << "Low invariant mass, W = " << W << " GeV!!";
1498  return false;
1499  }
1500  return true;
1501 }
1502 //____________________________________________________________________________
1503 
1504 
const double kPi
double fRvnCCm3
neugen&#39;s Rijk: vn, CC, multiplicity = 3
double fRvbpCCm3
neugen&#39;s Rijk: vbp, CC, multiplicity = 3
Basic constants.
double fPpi0eta
{Pi0 eta} production probability
int status
Definition: fabricate.py:1613
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:108
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
double ReWeightPt2(const PDGCodeList &pdgcv) const
#define pERROR
Definition: Messenger.h:60
const int kPdgLambda
Definition: PDGCodes.h:69
double MaxMult(const Interaction *i) const
int HadronShowerCharge(const Interaction *) const
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
double fAvn
offset in average charged hadron multiplicity = f(W) relation for vn
const char * p
Definition: xmltok.h:285
double fRvpCCm2
neugen&#39;s Rijk: vp, CC, multiplicity = 2
double fPeta
{eta eta} production probability
TClonesArray * Hadronize(const Interaction *) const
double fCvp
Levy function parameter for vp.
double fCvbp
Levy function parameter for vbp.
string P4AsString(const TLorentzVector *p)
Definition: PrintUtils.cxx:34
bool fReWeightDecays
Reweight phase space decays?
static const unsigned int kMaxUnweightDecayIterations
Definition: Controls.h:61
bool IsDarkMatter(int pdgc)
Definition: PDGUtils.cxx:125
TF1 * fBaryonPT2pdf
baryon pT^2 PDF
TParticlePDG * Probe(void) const
double Mass(Resonance_t res)
resonance mass (GeV)
TClonesArray * DecayMethod2(double W, const PDGCodeList &pdgv, bool reweight_decays) const
Definition: config.py:1
PDGCodeList * SelectParticles(const Interaction *) const
float Gamma() const
double fRvbnNCm3
neugen&#39;s Rijk: vbn, NC, multiplicity = 3
double fPpi0
{pi0 pi0 } production probability
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
double fPK0
{K0 K0bar} production probability
A list of PDG codes.
Definition: PDGCodeList.h:33
double fRvpNCm2
neugen&#39;s Rijk: vp, NC, multiplicity = 2
void HandleDecays(TClonesArray *particle_list) const
bool fForceNeuGenLimit
force upper hadronic multiplicity to NeuGEN limit
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
#define P(a, b, c, d, e, x)
bool PhaseSpaceDecay(TClonesArray &pl, TLorentzVector &pd, const PDGCodeList &pdgv, int offset=0, bool reweight=false) 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 nbins
Definition: cellShifts.C:15
int GenerateBaryonPdgCode(int mult, int maxQ, double W) const
Pure abstract base class. Defines the DecayModelI interface to be implemented by any algorithmic clas...
Definition: DecayModelI.h:31
const char * Find(const char *fname)
Definition: Icons.cxx:12
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
bool fForceDecays
force decays of unstable hadrons produced?
double fBvbp
slope in average charged hadron multiplicity = f(W) relation for vbp
PDGCodeList * GenerateHadronCodes(int mult, int maxQ, double W) 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
double fAvp
offset in average charged hadron multiplicity = f(W) relation for vp
Float_t E
Definition: plot.C:20
double fCvbn
Levy function parameter for vbn.
const int kPdgKM
Definition: PDGCodes.h:150
const double a
double fBhyperon
see above
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:116
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 kPdgKP
Definition: PDGCodes.h:149
double KNO(int nu, int nuc, double z) const
bool fGenerateWeighted
generate weighted events?
const int kPdgEta
Definition: PDGCodes.h:138
const int kPdgPiP
Definition: PDGCodes.h:135
const int kPdgPi0
Definition: PDGCodes.h:137
double fAvbn
offset in average charged hadron multiplicity = f(W) relation for vbn
TH1D * MultiplicityProb(const Interaction *, Option_t *opt="") const
#define pINFO
Definition: Messenger.h:63
bool AssertValidity(const Interaction *i) const
double fPKc
{K+ K- } production probability
const int kPdgAntiK0
Definition: PDGCodes.h:152
double fRvpCCm3
neugen&#39;s Rijk: vp, CC, multiplicity = 3
double fWeight
weight for generated event
double fRvbpNCm2
neugen&#39;s Rijk: vbp, NC, multiplicity = 2
Misc GENIE control constants.
z
Definition: test.py:28
TClonesArray * DecayBackToBack(double W, const PDGCodeList &pdgv) const
#define pWARN
Definition: Messenger.h:61
void Initialize(void) const
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
Simple printing utilities.
const int kPdgSigmaM
Definition: PDGCodes.h:73
double fRvnNCm3
neugen&#39;s Rijk: vn, NC, multiplicity = 3
double fCvn
Levy function parameter for vn.
TGenPhaseSpace fPhaseSpaceGenerator
a phase space generator
double fPhSpRwA
parameter for phase space decay reweighting
bool fUseIsotropic2BDecays
force isotropic, non-reweighted 2-body decays for consistency with neugen/daikon
double fBvn
slope in average charged hadron multiplicity = f(W) relation for vn
double fBvp
slope in average charged hadron multiplicity = f(W) relation for vp
const DecayModelI * fDecayer
decay algorithm
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 fAhyperon
parameter controlling strange baryon production probability via associated production (P=a+b*lnW^2) ...
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
const TLorentzVector * P4
4-momentum
Definition: DecayerInputs.h:36
double Wmin(void) const
Various utility methods common to hadronization models.
Singleton class to load & serve a TDatabasePDG.
Definition: PDGLibrary.h:30
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:324
double Weight(void) const
An abstract class. It avoids implementing the HadronizationModelI interface, leaving it for the concr...
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
double AverageChMult(int nu, int nuc, double W) const
bool IsNeutronOrProton(int pdgc)
Definition: PDGUtils.cxx:314
double fBvbn
slope in average charged hadron multiplicity = f(W) relation for vbn
bool fUseBaryonXfPt2Param
Generate baryon xF,pT2 from experimental parameterization?
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
exit(0)
const hit & b
Definition: hits.cxx:21
TClonesArray * DecayMethod1(double W, const PDGCodeList &pdgv, bool reweight_decays) const
assert(nhit_max >=nhit_nbins)
const int kPdgPiM
Definition: PDGCodes.h:136
double fWcut
configuration data common to all hadronizers
const int kPdgSigmaP
Definition: PDGCodes.h:71
const InitialState & InitState(void) const
Definition: Interaction.h:69
double fRvbnCCm3
neugen&#39;s Rijk: vbn, CC, multiplicity = 3
std::initializer_list< hep_hpc::hdf5::PropertyList > plist
bool fForceMinMult
force minimum multiplicity if (at low W) generated less?
static const unsigned int kMaxKNOHadSystIterations
Definition: Controls.h:57
const int kPdgProton
Definition: PDGCodes.h:65
#define pNOTICE
Definition: Messenger.h:62
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:67
double fPpic
{pi+ pi- } production probability
virtual TClonesArray * Decay(const DecayerInputs_t &inp) const =0
return a TClonesArray of TMCParticle objects (NOTE: all TMCParticle units in GeV^n [hbar=c=1]) ...
Double_t sum
Definition: plot.C:31
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:85
Float_t w
Definition: plot.C:20
#define W(x)
const int kPdgNeutron
Definition: PDGCodes.h:67
void Configure(const Registry &config)
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
TF1 * fBaryonXFpdf
baryon xF PDF
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
double fAvbp
offset in average charged hadron multiplicity = f(W) relation for vbp
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353