CharmHadronization.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 
16 #include <RVersion.h>
17 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,15,6)
18 #include <TMCParticle.h>
19 #else
20 #include <TMCParticle6.h>
21 #endif
22 #include <TPythia6.h>
23 #include <TVector3.h>
24 #include <TF1.h>
25 #include <TROOT.h>
26 
42 
43 using namespace genie;
44 using namespace genie::constants;
45 using namespace genie::controls;
46 
47 extern "C" void py1ent_(int *, int *, double *, double *, double *);
48 extern "C" void py2ent_(int *, int *, int *, double *);
49 
50 //____________________________________________________________________________
52 HadronizationModelI("genie::CharmHadronization")
53 {
54  this->Initialize();
55 }
56 //____________________________________________________________________________
58 HadronizationModelI("genie::CharmHadronization", config)
59 {
60  this->Initialize();
61 }
62 //____________________________________________________________________________
64 {
65  delete fCharmPT2pdf;
66  fCharmPT2pdf = 0;
67 
68  delete fD0FracSpl;
69  fD0FracSpl = 0;
70 
71  delete fDpFracSpl;
72  fDpFracSpl = 0;
73 
74  delete fDsFracSpl;
75  fDsFracSpl = 0;
76 }
77 //____________________________________________________________________________
79 {
80  fPythia = TPythia6::Instance();
81 }
82 //____________________________________________________________________________
84  const Interaction * interaction) const
85 {
86  LOG("CharmHad", pNOTICE) << "** Running CHARM hadronizer";
87 
88  PDGLibrary * pdglib = PDGLibrary::Instance();
90 
91  // ....................................................................
92  // Get information on the input event
93  //
94  const InitialState & init_state = interaction -> InitState();
95  const Kinematics & kinematics = interaction -> Kine();
96  const Target & target = init_state.Tgt();
97 
98  const TLorentzVector & p4Had = kinematics.HadSystP4();
99 
100  double Ev = init_state.ProbeE(kRfLab);
101  double W = kinematics.W(true);
102 
103  TVector3 beta = -1 * p4Had.BoostVector(); // boost vector for LAB' -> HCM'
104  TLorentzVector p4H(0,0,0,W); // hadronic system 4p @ HCM'
105 
106  double Eh = p4Had.Energy();
107 
108  LOG("CharmHad", pNOTICE) << "Ehad (LAB) = " << Eh << ", W = " << W;
109 
110  int nu_pdg = init_state.ProbePdg();
111  int nuc_pdg = target.HitNucPdg();
112 //int qpdg = target.HitQrkPdg();
113 //bool sea = target.HitSeaQrk();
114  bool isp = pdg::IsProton (nuc_pdg);
115  bool isn = pdg::IsNeutron(nuc_pdg);
116  bool isnu = pdg::IsNeutrino(nu_pdg);
117  bool isnub = pdg::IsAntiNeutrino(nu_pdg);
118  bool isdm = pdg::IsDarkMatter(nu_pdg);
119 
120  // ....................................................................
121  // Attempt to generate a charmed hadron & its 4-momentum
122  //
123 
124  TLorentzVector p4C(0,0,0,0);
125  int ch_pdg = -1;
126 
127  bool got_charmed_hadron = false;
128  unsigned int itry=0;
129 
130  while(itry++ < kRjMaxIterations && !got_charmed_hadron) {
131 
132  // Generate a charmed hadron PDG code
133  int pdg = this->GenerateCharmHadron(nu_pdg,Ev); // generate hadron
134  double mc = pdglib->Find(pdg)->Mass(); // lookup mass
135 
136  LOG("CharmHad", pNOTICE)
137  << "Trying charm hadron = " << pdg << "(m = " << mc << ")";
138 
139  if(mc>=W) continue; // dont' accept
140 
141  // Generate the charmed hadron energy based on the input
142  // fragmentation function
143  double z = fFragmFunc->GenerateZ(); // generate z(=Eh/Ev)
144  double Ec = z*Eh; // @ LAB'
145  double mc2 = TMath::Power(mc,2);
146  double Ec2 = TMath::Power(Ec,2);
147  double pc2 = Ec2-mc2;
148 
149  LOG("CharmHad", pINFO)
150  << "Trying charm hadron z = " << z << ", E = " << Ec;
151 
152  if(pc2<=0) continue;
153 
154  // Generate the charm hadron pT^2 and pL^2 (with respect to the
155  // hadronic system direction @ the LAB)
156  double ptc2 = fCharmPT2pdf->GetRandom();
157  double plc2 = Ec2 - ptc2 - mc2;
158  LOG("CharmHad", pINFO)
159  << "Trying charm hadron pT^2 (tranv to pHad) = " << ptc2;
160  if(plc2<0) continue;
161 
162  // Generate the charm hadron momentum components (@ LAB', z:\vec{pHad})
163  double ptc = TMath::Sqrt(ptc2);
164  double plc = TMath::Sqrt(plc2);
165  double phi = (2*kPi) * rnd->RndHadro().Rndm();
166  double pxc = ptc * TMath::Cos(phi);
167  double pyc = ptc * TMath::Sin(phi);
168  double pzc = plc;
169 
170  p4C.SetPxPyPzE(pxc,pyc,pzc,Ec); // @ LAB'
171 
172  // Boost charm hadron 4-momentum from the LAB' to the HCM' frame
173  //
174  LOG("CharmHad", pDEBUG)
175  << "Charm hadron p4 (@LAB') = " << utils::print::P4AsString(&p4C);
176 
177  p4C.Boost(beta);
178 
179  LOG("CharmHad", pDEBUG)
180  << "Charm hadron p4 (@HCM') = " << utils::print::P4AsString(&p4C);
181 
182  // Hadronic non-charm remnant 4p at HCM'
183  TLorentzVector p4 = p4H - p4C;
184  double wr = p4.M();
185  LOG("CharmHad", pINFO)
186  << "Invariant mass of remnant hadronic system= " << wr;
187  if(wr < kNucleonMass + kPionMass + kASmallNum) {
188  LOG("CharmHad", pINFO) << "Too small hadronic remnant mass!";
189  continue;
190  }
191 
192  ch_pdg = pdg;
193  got_charmed_hadron = true;
194 
195  LOG("CharmHad", pNOTICE)
196  << "Generated charm hadron = " << pdg << "(m = " << mc << ")";
197  LOG("CharmHad", pNOTICE)
198  << "Generated charm hadron z = " << z << ", E = " << Ec;
199  }
200 
201  // ....................................................................
202  // Check whether the code above had difficulty generating the charmed
203  // hadron near the W - threshold.
204  // If yes, attempt a phase space decay of a low mass charm hadron + nucleon
205  // pair that maintains the charge.
206  // That's a desperate solution but don't want to quit too early as that
207  // would distort the generated dsigma/dW distribution near threshold.
208  //
209  bool used_lowW_strategy = false;
210  int fs_nucleon_pdg = -1;
211  if(ch_pdg==-1 && W < 3.){
212  LOG("CharmHad", pNOTICE)
213  << "Had difficulty generating charm hadronic system near W threshold";
214  LOG("CharmHad", pNOTICE)
215  << "Trying an alternative strategy";
216 
217  double qfsl = interaction->FSPrimLepton()->Charge() / 3.;
218  double qinit = pdglib->Find(nuc_pdg)->Charge() / 3.;
219  int qhad = (int) (qinit - qfsl);
220 
221  int remn_pdg = -1;
222  int chrm_pdg = -1;
223 
224  //cc-only: qhad(nu) = +1,+2, qhad(nubar)= -1,0
225  //
226  if(qhad == 2) {
227  chrm_pdg = kPdgDP; remn_pdg = kPdgProton;
228  } else if(qhad == 1) {
229  if(rnd->RndHadro().Rndm() > 0.5) {
230  chrm_pdg = kPdgD0; remn_pdg = kPdgProton;
231  } else {
232  chrm_pdg = kPdgDP; remn_pdg = kPdgNeutron;
233  }
234  } else if(qhad == 0) {
235  chrm_pdg = kPdgAntiD0; remn_pdg = kPdgNeutron;
236  } else if(qhad == -1) {
237  chrm_pdg = kPdgDM; remn_pdg = kPdgNeutron;
238  }
239 
240  double mc = pdglib->Find(chrm_pdg)->Mass();
241  double mn = pdglib->Find(remn_pdg)->Mass();
242 
243  if(mc+mn < W) {
244  // Set decay
245  double mass[2] = {mc, mn};
246  bool permitted = fPhaseSpaceGenerator.SetDecay(p4H, 2, mass);
247  assert(permitted);
248 
249  // Get the maximum weight
250  double wmax = -1;
251  for(int i=0; i<200; i++) {
252  double w = fPhaseSpaceGenerator.Generate();
253  wmax = TMath::Max(wmax,w);
254  }
255 
256  if(wmax>0) {
257  wmax *= 2;
258 
259  // Generate unweighted decay
260  bool accept_decay=false;
261  unsigned int idecay_try=0;
262  while(!accept_decay)
263  {
264  idecay_try++;
265 
266  if(idecay_try>kMaxUnweightDecayIterations) {
267  LOG("CharmHad", pWARN)
268  << "Couldn't generate an unweighted phase space decay after "
269  << idecay_try << " attempts";
270  }
271  double w = fPhaseSpaceGenerator.Generate();
272  if(w > wmax) {
273  LOG("CharmHad", pWARN)
274  << "Decay weight = " << w << " > max decay weight = " << wmax;
275  }
276  double gw = wmax * rnd->RndHadro().Rndm();
277  accept_decay = (gw<=w);
278 
279  if(accept_decay) {
280  used_lowW_strategy = true;
281  TLorentzVector * p4 = fPhaseSpaceGenerator.GetDecay(0);
282  p4C = *p4;
283  ch_pdg = chrm_pdg;
284  fs_nucleon_pdg = remn_pdg;
285  }
286  } // decay loop
287  }//wmax>0
288 
289  }// allowed decay
290  } // alt low-W strategy
291 
292  // ....................................................................
293  // Check success in generating the charm hadron & compute 4p for
294  // remnant system
295  //
296  if(ch_pdg==-1){
297  LOG("CharmHad", pWARN)
298  << "Couldn't generate charm hadron for: " << *interaction;
299  return 0;
300  }
301 
302  TLorentzVector p4R = p4H - p4C;
303  double WR = p4R.M();
304  double MC = pdglib->Find(ch_pdg)->Mass();
305 
306  LOG("CharmHad", pNOTICE) << "Remnant hadronic system mass = " << WR;
307 
308  // ....................................................................
309  // Handle case where the user doesn't want to remnant system to be
310  // hadronized (add as 'hadronic blob')
311  //
312  if(fCharmOnly) {
313  // Create particle list (fragmentation record)
314  TClonesArray * particle_list = new TClonesArray("TMCParticle", 2);
315  particle_list->SetOwner(true);
316 
317  // insert the generated particles
318  new ((*particle_list)[0]) TMCParticle (1,ch_pdg,
319  -1,-1,-1, p4C.Px(),p4C.Py(),p4C.Pz(),p4C.E(),MC, 0,0,0,0,0);
320  new ((*particle_list)[1]) TMCParticle (1,kPdgHadronicBlob,
321  -1,-1,-1, p4R.Px(),p4R.Py(),p4R.Pz(),p4R.E(),WR, 0,0,0,0,0);
322 
323  return particle_list;
324  }
325 
326  // ....................................................................
327  // Handle case where the remnant system is already known and doesn't
328  // have to be hadronized. That happens when (close to the W threshold)
329  // the hadronic system was generated by a simple 2-body decay
330  //
331  if(used_lowW_strategy) {
332  // Create particle list (fragmentation record)
333  TClonesArray * particle_list = new TClonesArray("TMCParticle", 3);
334  particle_list->SetOwner(true);
335 
336  // insert the generated particles
337  new ((*particle_list)[0]) TMCParticle (1,ch_pdg,
338  -1,-1,-1, p4C.Px(),p4C.Py(),p4C.Pz(),p4C.E(),MC, 0,0,0,0,0);
339  new ((*particle_list)[1]) TMCParticle (11,kPdgHadronicBlob,
340  -1,2,2, p4R.Px(),p4R.Py(),p4R.Pz(),p4R.E(),WR, 0,0,0,0,0);
341  new ((*particle_list)[2]) TMCParticle (1,fs_nucleon_pdg,
342  1,-1,-1, p4R.Px(),p4R.Py(),p4R.Pz(),p4R.E(),WR, 0,0,0,0,0);
343 
344  return particle_list;
345  }
346 
347  // ....................................................................
348  // --------------------------------------------------------------------
349  // Hadronize non-charm hadronic blob using PYTHIA/JETSET
350  // --------------------------------------------------------------------
351  // ....................................................................
352 
353  // Create output event record
354  // Insert the generated charm hadron & the hadronic (non-charm) blob.
355  // In this case the hadronic blob is entered as a pre-fragm. state.
356 
357  TClonesArray * particle_list = new TClonesArray("TMCParticle");
358  particle_list->SetOwner(true);
359 
360  new ((*particle_list)[0]) TMCParticle (1,ch_pdg,
361  -1,-1,-1, p4C.Px(),p4C.Py(),p4C.Pz(),p4C.E(),MC, 0,0,0,0,0);
362  new ((*particle_list)[1]) TMCParticle (11,kPdgHadronicBlob,
363  -1,2,3, p4R.Px(),p4R.Py(),p4R.Pz(),p4R.E(),WR, 0,0,0,0,0);
364 
365  unsigned int rpos =2; // offset in event record
366 
367  bool use_pythia = (WR>1.5);
368 
369 /*
370  // Determining quark systems to input to PYTHIA based on simple quark model
371  // arguments
372  //
373  // Neutrinos
374  // ------------------------------------------------------------------
375  // Scattering off valence q
376  // ..................................................................
377  // p: [uu]+d
378  // |--> c --> D0 <c+\bar(u)> : [u]
379  // --> D+ <c+\bar(d)> : [d]
380  // --> Ds+ <c+\bar(s)> : [s]
381  // --> Lamda_c+ <c+ud > : [\bar(ud)]
382  //
383  // (for n: [uu] -> 50%[ud]_{0} + 50%[ud]_{1})
384  //
385  // Scattering off sea q
386  // ..................................................................
387  // p: [uud] + [\bar(d)]d (or)
388  // [\bar(s)]s
389  // |--> c --> D0 <c+\bar(u)> : [u]
390  // --> D+ <c+\bar(d)> : [d]
391  // --> Ds+ <c+\bar(s)> : [s]
392  // --> Lamda_c+ <c+ud > : [\bar(ud)]
393  // Anti-Neutrinos
394  // ------------------------------------------------------------------
395  // Scattering off sea q
396  // ..................................................................
397  // p: [uud] + [d] \bar(d) (or)
398  // [s] \bar(s)
399  // |----> \bar(c) --> \bar(D0) <\bar(c)+u> : [\bar(u)]
400  // --> D- <\bar(c)+d> : [\bar(d)]
401  // --> Ds- <\bar(c)+s> : [\bar(s)]
402  // [Summary]
403  // Qq
404  // | v + p [val/d] --> D0 + { u uu }(+2) / u,uu
405  // | v + p [val/d] --> D+ + { d uu }(+1) / d,uu
406  // | v + p [val/d] --> Ds+ + { s uu }(+1) / s,uu
407  // | v + p [val/d] --> Lc+ + { \bar(ud) uu }(+1) / \bar(d),u
408  // | v + n [val/d] --> D0 + { u ud }(+1) / u,ud
409  // | v + n [val/d] --> D+ + { d ud }( 0) / d,ud
410  // | v + n [val/d] --> Ds+ + { s ud }( 0) / s,ud
411  // | v + n [val/d] --> Lc+ + { \bar(ud) ud }( 0) / \bar(d),d
412  // | v + p [sea/d] --> D0 + { uud \bar(d) u }(+2) / u,uu
413  // | v + p [sea/d] --> D+ + { uud \bar(d) d }(+1) / d,uu
414  // | v + p [sea/d] --> Ds+ + { uud \bar(d) s }(+1) / s,uu
415  // | v + p [sea/d] --> Lc+ + { uud \bar(d) \bar(ud) }(+1) / \bar(d),u
416  // | v + n [sea/d] --> D0 + { udd \bar(d) u }(+1) / u,ud
417  // | v + n [sea/d] --> D+ + { udd \bar(d) d }( 0) / d,ud
418  // | v + n [sea/d] --> Ds+ + { udd \bar(d) s }( 0) / s,ud
419  // | v + n [sea/d] --> Lc+ + { udd \bar(d) \bar(ud) }( 0) / \bar(d),d
420  // | v + p [sea/s] --> D0 + { uud \bar(s) u }(+2) / u,uu
421  // | v + p [sea/s] --> D+ + { uud \bar(s) d }(+1) / d,uu
422  // | v + p [sea/s] --> Ds+ + { uud \bar(s) s }(+1) / s,uu
423  // | v + p [sea/s] --> Lc+ + { uud \bar(s) \bar(ud) }(+1) / \bar(d),u
424  // | v + n [sea/s] --> D0 + { udd \bar(s) u }(+1) / u,ud
425  // | v + n [sea/s] --> D+ + { udd \bar(s) d }( 0) / d,ud
426  // | v + n [sea/s] --> Ds+ + { udd \bar(s) s }( 0) / s,ud
427  // | v + n [sea/s] --> Lc+ + { udd \bar(s) \bar(ud) }( 0) / \bar(d),d
428 
429  // | \bar(v) + p [sea/\bar(d)] --> \bar(D0) + { uud d \bar(u) }( 0) / d,ud
430  // | \bar(v) + p [sea/\bar(d)] --> D- + { uud d \bar(d) }(+1) / d,uu
431  // | \bar(v) + p [sea/\bar(d)] --> Ds- + { uud d \bar(s) }(+1) / d,uu
432  // | \bar(v) + n [sea/\bar(d)] --> \bar(D0) + { udd d \bar(u) }(-1) / d,dd
433  // | \bar(v) + n [sea/\bar(d)] --> D- + { udd d \bar(d) }( 0) / d,ud
434  // | \bar(v) + n [sea/\bar(d)] --> Ds- + { udd d \bar(s) }( 0) / d,ud
435  // | \bar(v) + p [sea/\bar(s)] --> \bar(D0) + { uud s \bar(u) }( 0) / d,ud
436  // | \bar(v) + p [sea/\bar(s)] --> D- + { uud s \bar(d) }(+1) / d,uu
437  // | \bar(v) + p [sea/\bar(s)] --> Ds- + { uud s \bar(s) }(+1) / d,uu
438  // | \bar(v) + n [sea/\bar(s)] --> \bar(D0) + { udd s \bar(u) }(-1) / d,dd
439  // | \bar(v) + n [sea/\bar(s)] --> D- + { udd s \bar(d) }( 0) / d,ud
440  // | \bar(v) + n [sea/\bar(s)] --> Ds- + { udd s \bar(s) }( 0) / d,ud
441  //
442  //
443  // Taking some short-cuts below :
444  // In the process of obtaining 2 q systems (while conserving the charge) I might tread
445  // d,s or \bar(d),\bar(s) as the same
446  // In the future I should perform the first steps of the multi-q hadronization manualy
447  // (to reduce the number of q's input to PYTHIA) or use py3ent_(), py4ent_() ...
448  //
449 */
450 
451  if(use_pythia) {
452  int qrkSyst1 = 0;
453  int qrkSyst2 = 0;
454  if(isnu||isdm) { // neutrinos
455  if(ch_pdg==kPdgLambdaPc) {
456  if(isp) { qrkSyst1 = kPdgAntiDQuark; qrkSyst2 = kPdgUQuark; }
457  if(isn) { qrkSyst1 = kPdgAntiDQuark; qrkSyst2 = kPdgDQuark; }
458  } else {
459  if(isp) { qrkSyst2 = kPdgUUDiquarkS1; }
460  if(isn) { qrkSyst2 = (rnd->RndHadro().Rndm()<0.5) ? kPdgUDDiquarkS0 : kPdgUDDiquarkS1; }
461  if (ch_pdg==kPdgD0 ) { qrkSyst1 = kPdgUQuark; }
462  if (ch_pdg==kPdgDP ) { qrkSyst1 = kPdgDQuark; }
463  if (ch_pdg==kPdgDPs ) { qrkSyst1 = kPdgSQuark; }
464  }
465  }
466  if(isnub) { // antineutrinos
467  qrkSyst1 = kPdgDQuark;
468  if (isp && ch_pdg==kPdgAntiD0) { qrkSyst2 = (rnd->RndHadro().Rndm()<0.5) ? kPdgUDDiquarkS0 : kPdgUDDiquarkS1; }
469  if (isp && ch_pdg==kPdgDM ) { qrkSyst2 = kPdgUUDiquarkS1; }
470  if (isp && ch_pdg==kPdgDMs ) { qrkSyst2 = kPdgUUDiquarkS1; }
471  if (isn && ch_pdg==kPdgAntiD0) { qrkSyst2 = kPdgDDDiquarkS1; }
472  if (isn && ch_pdg==kPdgDM ) { qrkSyst2 = (rnd->RndHadro().Rndm()<0.5) ? kPdgUDDiquarkS0 : kPdgUDDiquarkS1; }
473  if (isn && ch_pdg==kPdgDMs ) { qrkSyst2 = (rnd->RndHadro().Rndm()<0.5) ? kPdgUDDiquarkS0 : kPdgUDDiquarkS1; }
474  }
475  if(qrkSyst1 == 0 && qrkSyst2 == 0) {
476  LOG("CharmHad", pWARN)
477  << "Couldn't generate quark systems for PYTHIA in: " << *interaction;
478  return 0;
479  }
480 
481  //
482  // Run PYTHIA for the hadronization of remnant system
483  //
484  fPythia->SetMDCY(fPythia->Pycomp(kPdgPi0), 1,0); // don't decay pi0
485  fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaM), 1,1); // decay Delta+
486  fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_Delta0), 1,1); // decay Delta++
487  fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaP), 1,1); // decay Delta++
488  fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaPP), 1,1); // decay Delta++
489 // fPythia->SetMDCY(fPythia->Pycomp(kPdgDeltaP), 1,1); // decay Delta+
490 // fPythia->SetMDCY(fPythia->Pycomp(kPdgDeltaPP), 1,1); // decay Delta++
491 
492  int ip = 0;
493  py2ent_(&ip, &qrkSyst1, &qrkSyst2, &WR); // hadronize
494 
495  fPythia->SetMDCY(fPythia->Pycomp(kPdgPi0),1,1); // restore
496 
497  //-- Get PYTHIA's LUJETS event record
498  TClonesArray * remnants = 0;
499  fPythia->GetPrimaries();
500  remnants = dynamic_cast<TClonesArray *>(fPythia->ImportParticles("All"));
501  if(!remnants) {
502  LOG("CharmHad", pWARN) << "Couldn't hadronize (non-charm) remnants!";
503  return 0;
504  }
505 
506  // PYTHIA performs the hadronization at the *remnant hadrons* centre of mass
507  // frame (not the hadronic centre of mass frame).
508  // Boost all hadronic blob fragments to the HCM', fix their mother/daughter
509  // assignments and add them to the fragmentation record.
510 
511  TVector3 rmnbeta = +1 * p4R.BoostVector(); // boost velocity
512 
513  TMCParticle * remn = 0; // remnant
514  TMCParticle * bremn = 0; // boosted remnant
515  TIter remn_iter(remnants);
516  while( (remn = (TMCParticle *) remn_iter.Next()) ) {
517 
518  // insert and get a pointer to inserted object for mods
519  bremn = new ((*particle_list)[rpos++]) TMCParticle (*remn);
520 
521  // boost
522  TLorentzVector p4(remn->GetPx(),remn->GetPy(),remn->GetPz(),remn->GetEnergy());
523  p4.Boost(rmnbeta);
524  bremn -> SetPx (p4.Px());
525  bremn -> SetPy (p4.Py());
526  bremn -> SetPz (p4.Pz());
527  bremn -> SetEnergy (p4.E() );
528 
529  // handle insertion of charmed hadron
530  int jp = bremn->GetParent();
531  int ifc = bremn->GetFirstChild();
532  int ilc = bremn->GetLastChild();
533  bremn -> SetParent ( (jp == 0 ? 1 : jp +1) );
534  bremn -> SetFirstChild ( (ifc == 0 ? -1 : ifc+1) );
535  bremn -> SetLastChild ( (ilc == 0 ? -1 : ilc+1) );
536  }
537  } // use_pythia
538 
539  // ....................................................................
540  // Hadronizing low-W non-charm hadronic blob using a phase space decay
541  // ....................................................................
542 
543  else {
544  // Just a small fraction of events (low-W remnant syste) causing trouble
545  // to PYTHIA/JETSET
546  // Set a 2-body N+pi system that matches the remnant system charge and
547  // do a simple phase space decay
548  //
549  // q(remn) remn/syst
550  // +2 : (p pi+)
551  // +1 : 50%(p pi0) + 50% n pi+
552  // 0 : 50%(p pi-) + 50% n pi0
553  // -1 : (n pi-)
554  //
555  double qfsl = interaction->FSPrimLepton()->Charge() / 3.;
556  double qinit = pdglib->Find(nuc_pdg)->Charge() / 3.;
557  double qch = pdglib->Find(ch_pdg)->Charge() / 3.;
558  int Q = (int) (qinit - qfsl - qch); // remnant hadronic system charge
559 
560  bool allowdup=true;
561  PDGCodeList pd(allowdup);
562  if(Q==2) {
564  else if (Q==1) {
565  if(rnd->RndHadro().Rndm()<0.5) {
567  else {
569  }
570  else if (Q==0) {
571  if(rnd->RndHadro().Rndm()<0.5) {
573  else {
575  }
576  else if (Q==-1) {
578 
579  double mass[2] = {
580  pdglib->Find(pd[0])->Mass(), pdglib->Find(pd[1])->Mass()
581  };
582 
583  // Set the decay
584  bool permitted = fPhaseSpaceGenerator.SetDecay(p4R, 2, mass);
585  if(!permitted) {
586  LOG("CharmHad", pERROR) << " *** Phase space decay is not permitted";
587  return 0;
588  }
589  // Get the maximum weight
590  double wmax = -1;
591  for(int i=0; i<200; i++) {
592  double w = fPhaseSpaceGenerator.Generate();
593  wmax = TMath::Max(wmax,w);
594  }
595  if(wmax<=0) {
596  LOG("CharmHad", pERROR) << " *** Non-positive maximum weight";
597  LOG("CharmHad", pERROR) << " *** Can not generate an unweighted phase space decay";
598  return 0;
599  }
600 
601  LOG("CharmHad", pINFO)
602  << "Max phase space gen. weight @ current hadronic system: " << wmax;
603 
604  // *** generating an un-weighted decay ***
605  wmax *= 1.3;
606  bool accept_decay=false;
607  unsigned int idectry=0;
608  while(!accept_decay)
609  {
610  idectry++;
611  if(idectry>kMaxUnweightDecayIterations) {
612  // report, clean-up and return
613  LOG("Char,Had", pWARN)
614  << "Couldn't generate an unweighted phase space decay after "
615  << itry << " attempts";
616  return 0;
617  }
618  double w = fPhaseSpaceGenerator.Generate();
619  if(w > wmax) {
620  LOG("CharmHad", pWARN)
621  << "Decay weight = " << w << " > max decay weight = " << wmax;
622  }
623  double gw = wmax * rnd->RndHadro().Rndm();
624  accept_decay = (gw<=w);
625  LOG("CharmHad", pDEBUG)
626  << "Decay weight = " << w << " / R = " << gw << " - accepted: " << accept_decay;
627  }
628  for(unsigned int i=0; i<2; i++) {
629  int pdgc = pd[i];
630  TLorentzVector * p4d = fPhaseSpaceGenerator.GetDecay(i);
631  new ( (*particle_list)[rpos+i] ) TMCParticle(
632  1,pdgc,1,-1,-1,p4d->Px(),p4d->Py(),p4d->Pz(),p4d->Energy(),
633  mass[i],0,0,0,0,0);
634  }
635  }
636 
637  //-- Print & return the fragmentation record
638  //utils::fragmrec::Print(particle_list);
639  return particle_list;
640 }
641 //____________________________________________________________________________
642 double CharmHadronization::Weight(void) const
643 {
644  return 1.; // does not generate weighted events
645 }
646 //____________________________________________________________________________
648  const Interaction * /*interaction*/) const
649 {
650  return 0;
651 }
652 //____________________________________________________________________________
654  const Interaction * /*interaction*/, Option_t * /*opt*/) const
655 {
656  return 0;
657 }
658 //____________________________________________________________________________
659 int CharmHadronization::GenerateCharmHadron(int nu_pdg, double EvLab) const
660 {
661  // generate a charmed hadron pdg code using a charm fraction table
662 
664  double r = rnd->RndHadro().Rndm();
665 
666  if(pdg::IsNeutrino(nu_pdg)) {
667  double tf = 0;
668  if (r < (tf+=fD0FracSpl->Evaluate(EvLab))) return kPdgD0; // D^0
669  else if (r < (tf+=fDpFracSpl->Evaluate(EvLab))) return kPdgDP; // D^+
670  else if (r < (tf+=fDsFracSpl->Evaluate(EvLab))) return kPdgDPs; // Ds^+
671  else return kPdgLambdaPc; // Lamda_c^+
672 
673  } else if(pdg::IsAntiNeutrino(nu_pdg)) {
674  if (r < fD0BarFrac) return kPdgAntiD0;
675  else if (r < fD0BarFrac+fDmFrac) return kPdgDM;
676  else return kPdgDMs;
677  }
678 
679  LOG("CharmHad", pERROR) << "Could not generate a charm hadron!";
680  return 0;
681 }
682 //____________________________________________________________________________
684 {
685  Algorithm::Configure(config);
686  this->LoadConfig();
687 }
688 //____________________________________________________________________________
690 {
691  Algorithm::Configure(config);
692  this->LoadConfig();
693 }
694 //____________________________________________________________________________
696 {
697 
698  bool hadronize_remnants = true ;
699  GetParam( "HadronizeRemnants", hadronize_remnants, false ) ;
700 
701  fCharmOnly = ! hadronize_remnants ;
702 
703  //-- Get a fragmentation function
704  fFragmFunc = dynamic_cast<const FragmentationFunctionI *> (
705  this->SubAlg("FragmentationFunc"));
707 
708  fCharmPT2pdf = new TF1("fCharmPT2pdf", "exp(-0.213362-6.62464*x)",0,0.6);
709  // stop ROOT from deleting this object of its own volition
710  gROOT->GetListOfFunctions()->Remove(fCharmPT2pdf);
711 
712  // neutrino charm fractions: D^0, D^+, Ds^+ (remainder: Lamda_c^+)
713  //
714  const int nc = 15;
715  double ec[nc] = {0.,5.,10.,15.,20.,25.,30.,35.,40.,50.,60.,70.,80.,90.,100.};
716 
717  double d0frac[nc] = { .000, .320, .460, .500, .520, .530, .540, .540,
718  .540, .550, .550, .560, .570, .580, .600 };
719  double dpfrac[nc] = { .000, .120, .180, .200, .200, .210, .210, .210,
720  .210, .210, .220, .220, .220, .230, .230 };
721  double dsfrac[nc] = { .000, .054, .078, .130, .130, .140, .140, .140,
722  .140, .140, .140, .140, .140, .150, .150 };
723 
724  fD0FracSpl = new Spline(nc, ec, d0frac);
725  fDpFracSpl = new Spline(nc, ec, dpfrac);
726  fDsFracSpl = new Spline(nc, ec, dsfrac);
727 
728  // anti-neutrino charm fractions: bar(D^0), D^-, (remainder: Ds^-)
729  //
730  fD0BarFrac = 0.667;
731  fDmFrac = 0.222;
732 }
733 //____________________________________________________________________________
const int kPdgAntiD0
Definition: PDGCodes.h:161
const double kPi
const int kPdgP33m1232_DeltaPP
Definition: PDGCodes.h:91
const int kPdgDPs
Definition: PDGCodes.h:162
const int kPdgUUDiquarkS1
Definition: PDGCodes.h:58
double W(bool selected=false) const
Definition: Kinematics.cxx:167
double fDmFrac
nubar D- charm fraction
Basic constants.
TPythia6 * fPythia
remnant (non-charm) hadronizer
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
Spline * fDpFracSpl
nu charm fraction vs Ev: D+
#define pERROR
Definition: Messenger.h:60
static const double kNucleonMass
Definition: Constants.h:78
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
int HitNucPdg(void) const
Definition: Target.cxx:321
virtual double GenerateZ(void) const =0
const int kPdgHadronicBlob
Definition: PDGCodes.h:188
A numeric analysis tool class for interpolating 1-D functions.
Definition: Spline.h:47
void Initialize(void) const
define the HadronizationModelI interface
const int kPdgUQuark
Definition: PDGCodes.h:42
void py1ent_(int *, int *, double *, double *, double *)
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
string P4AsString(const TLorentzVector *p)
Definition: PrintUtils.cxx:34
static const unsigned int kMaxUnweightDecayIterations
Definition: Controls.h:61
bool IsDarkMatter(int pdgc)
Definition: PDGUtils.cxx:125
const TLorentzVector & HadSystP4(void) const
Definition: Kinematics.h:67
Definition: config.py:1
TString ip
Definition: loadincs.C:5
Double_t beta
Timing fit.
double Evaluate(double x) const
Definition: Spline.cxx:362
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
PDGCodeList * SelectParticles(const Interaction *) const
void py2ent_(int *, int *, int *, double *)
const int kPdgSQuark
Definition: PDGCodes.h:46
A list of PDG codes.
Definition: PDGCodeList.h:33
const int kPdgP33m1232_DeltaP
Definition: PDGCodes.h:90
const int kPdgP33m1232_DeltaM
Definition: PDGCodes.h:88
TClonesArray * Hadronize(const Interaction *) const
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:304
Summary information for an interaction.
Definition: Interaction.h:56
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
const int kPdgUDDiquarkS1
Definition: PDGCodes.h:57
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:116
const int kPdgAntiDQuark
Definition: PDGCodes.h:45
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:41
const FragmentationFunctionI * fFragmFunc
charm hadron fragmentation func
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
int ProbePdg(void) const
Definition: InitialState.h:65
double fD0BarFrac
nubar {D0} charm fraction
const int kPdgLambdaPc
Definition: PDGCodes.h:83
const int kPdgPiP
Definition: PDGCodes.h:135
const int kPdgPi0
Definition: PDGCodes.h:137
void Configure(const Registry &config)
const int kPdgDQuark
Definition: PDGCodes.h:44
static const double kASmallNum
Definition: Controls.h:40
#define pINFO
Definition: Messenger.h:63
const int kPdgUDDiquarkS0
Definition: PDGCodes.h:56
Misc GENIE control constants.
z
Definition: test.py:28
#define pWARN
Definition: Messenger.h:61
const int kPdgP33m1232_Delta0
Definition: PDGCodes.h:89
bool fCharmOnly
don&#39;t hadronize non-charm blob
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
Pure abstract base class. Defines the FragmentationFunctionI interface to be implemented by any algor...
const int kPdgDP
Definition: PDGCodes.h:158
TRandom3 & RndHadro(void) const
rnd number generator used by hadronization models
Definition: RandomGen.h:54
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
static const double kPionMass
Definition: Constants.h:74
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
int GenerateCharmHadron(int nupdg, double EvLab) const
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
const int kPdgDMs
Definition: PDGCodes.h:163
Pure abstract base class. Defines the HadronizationModelI interface to be implemented by any algorith...
const int kPdgDM
Definition: PDGCodes.h:159
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
const int kPdgPiM
Definition: PDGCodes.h:136
Spline * fDsFracSpl
nu charm fraction vs Ev: Ds+
const int kPdgProton
Definition: PDGCodes.h:65
#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
static const unsigned int kRjMaxIterations
Definition: Controls.h:27
void kinematics()
Definition: kinematics.C:10
Float_t w
Definition: plot.C:20
#define W(x)
double ProbeE(RefFrame_t rf) const
const int kPdgNeutron
Definition: PDGCodes.h:67
const int kPdgD0
Definition: PDGCodes.h:160
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
Spline * fD0FracSpl
nu charm fraction vs Ev: D0
TF1 * fCharmPT2pdf
charm hadron pT^2 pdf
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:67
TGenPhaseSpace fPhaseSpaceGenerator
a phase space generator
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
TH1D * MultiplicityProb(const Interaction *, Option_t *opt="") const
const int kPdgDDDiquarkS1
Definition: PDGCodes.h:55
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353