GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AGCharm2019.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
8 
9  Changes required to implement the GENIE Boosted Dark Matter module
10  were installed by Josh Berger (Univ. of Wisconsin)
11 */
12 //____________________________________________________________________________
13 
14 #include <RVersion.h>
15 #include <TPythia6.h>
16 #include <TVector3.h>
17 #include <TF1.h>
18 #include <TROOT.h>
19 
20 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,15,6)
21 #include <TMCParticle.h>
22 #else
23 #include <TMCParticle6.h>
24 #endif
25 
46 
47 using namespace genie;
48 using namespace genie::constants;
49 using namespace genie::controls;
50 
51 extern "C" void py1ent_(int *, int *, double *, double *, double *);
52 extern "C" void py2ent_(int *, int *, int *, double *);
53 
54 //____________________________________________________________________________
56 EventRecordVisitorI("genie::AGCharm2019")
57 {
58  this->Initialize();
59 }
60 //____________________________________________________________________________
61 AGCharm2019::AGCharm2019(string config) :
62 EventRecordVisitorI("genie::AGCharm2019", config)
63 {
64  this->Initialize();
65 }
66 //____________________________________________________________________________
68 {
69  delete fCharmPT2pdf;
70  fCharmPT2pdf = 0;
71 
72  delete fD0FracSpl;
73  fD0FracSpl = 0;
74 
75  delete fDpFracSpl;
76  fDpFracSpl = 0;
77 
78  delete fDsFracSpl;
79  fDsFracSpl = 0;
80 }
81 //____________________________________________________________________________
82 void AGCharm2019::Initialize(void) const
83 {
84  fPythia = TPythia6::Instance();
85 }
86 //____________________________________________________________________________
88 {
89  Interaction * interaction = event->Summary();
90  TClonesArray * particle_list = this->Hadronize(interaction);
91 
92  if(! particle_list ) {
93  LOG("AGCharm2019", pWARN) << "Got an empty particle list. Hadronizer failed!";
94  LOG("AGCharm2019", pWARN) << "Quitting the current event generation thread";
95 
96  event->EventFlags()->SetBitNumber(kHadroSysGenErr, true);
97 
99  exception.SetReason("Could not simulate the hadronic system");
100  exception.SwitchOnFastForward();
101  throw exception;
102 
103  return;
104  }
105 
106  int mom = event->FinalStateHadronicSystemPosition();
107  assert(mom!=-1);
108 
109  // find the proper status for the particles we are going to put in event record
110  bool is_nucleus = interaction->InitState().Tgt().IsNucleus();
111  GHepStatus_t istfin = (is_nucleus) ?
113 
114  // retrieve the hadronic blob lorentz boost
115  // Because Hadronize() returned particles not in the LAB reference frame
116  const TLorentzVector * had_syst = event -> Particle(mom) -> P4() ;
117  TVector3 beta( 0., 0., had_syst -> P()/ had_syst -> Energy() ) ;
118 
119  // Vector defining rotation from LAB to LAB' (z:= \vec{phad})
120  TVector3 unitvq = had_syst -> Vect().Unit();
121 
122  GHepParticle * neutrino = event->Probe();
123  const TLorentzVector & vtx = *(neutrino->X4());
124 
125  GHepParticle * particle = 0;
126  TIter particle_iter(particle_list);
127  while ((particle = (GHepParticle *) particle_iter.Next())) {
128 
129  int pdgc = particle -> Pdg() ;
130 
131  // bring the particle in the LAB reference frame
132  particle -> P4() -> Boost( beta ) ;
133  particle -> P4() -> RotateUz( unitvq ) ;
134 
135  // set the proper status according to a number of things:
136  // interaction on a nucleaus or nucleon, particle type
137  GHepStatus_t ist = ( particle -> Status() ==1 ) ? istfin : kIStDISPreFragmHadronicState;
138 
139  // handle gammas, and leptons that might come from internal pythia decays
140  // mark them as final state particles
141  bool not_hadron = ( pdgc == kPdgGamma ||
142  pdg::IsNeutralLepton(pdgc) ||
143  pdg::IsChargedLepton(pdgc) ) ;
144 
145  if(not_hadron) { ist = kIStStableFinalState; }
146  particle -> SetStatus( ist ) ;
147 
148  int im = mom + 1 + particle -> FirstMother() ;
149  int ifc = ( particle -> FirstDaughter() == -1) ? -1 : mom + 1 + particle -> FirstDaughter();
150  int ilc = ( particle -> LastDaughter() == -1) ? -1 : mom + 1 + particle -> LastDaughter();
151 
152  particle -> SetFirstMother( im ) ;
153  if ( ifc > -1 ) {
154  particle -> SetFirstDaughter( ifc ) ;
155  particle -> SetLastDaughter( ilc ) ;
156  }
157 
158  // the Pythia particle position is overridden
159  particle -> SetPosition( vtx ) ;
160 
161  event->AddParticle(*particle);
162  }
163 
164  particle_list -> Delete() ;
165  delete particle_list ;
166 
167  // update the weight of the event
168  event -> SetWeight ( Weight() * event->Weight() );
169 
170 }
171 //____________________________________________________________________________
172 TClonesArray * AGCharm2019::Hadronize(
173  const Interaction * interaction) const
174 {
175  LOG("CharmHad", pNOTICE) << "** Running CHARM hadronizer";
176 
177  PDGLibrary * pdglib = PDGLibrary::Instance();
178  RandomGen * rnd = RandomGen::Instance();
179 
180  // ....................................................................
181  // Get information on the input event
182  //
183  const InitialState & init_state = interaction -> InitState();
184  const Kinematics & kinematics = interaction -> Kine();
185  const Target & target = init_state.Tgt();
186 
187  const TLorentzVector & p4Had = kinematics.HadSystP4();
188 
189  double Ev = init_state.ProbeE(kRfLab);
190  double W = kinematics.W(true);
191 
192  TVector3 beta = -1 * p4Had.BoostVector(); // boost vector for LAB' -> HCM'
193  TLorentzVector p4H(0,0,0,W); // hadronic system 4p @ HCM'
194 
195  double Eh = p4Had.Energy();
196 
197  LOG("CharmHad", pNOTICE) << "Ehad (LAB) = " << Eh << ", W = " << W;
198 
199  int nu_pdg = init_state.ProbePdg();
200  int nuc_pdg = target.HitNucPdg();
201 //int qpdg = target.HitQrkPdg();
202 //bool sea = target.HitSeaQrk();
203  bool isp = pdg::IsProton (nuc_pdg);
204  bool isn = pdg::IsNeutron(nuc_pdg);
205  bool isnu = pdg::IsNeutrino(nu_pdg);
206  bool isnub = pdg::IsAntiNeutrino(nu_pdg);
207  bool isdm = pdg::IsDarkMatter(nu_pdg);
208 
209  // ....................................................................
210  // Attempt to generate a charmed hadron & its 4-momentum
211  //
212 
213  TLorentzVector p4C(0,0,0,0);
214  int ch_pdg = -1;
215 
216  bool got_charmed_hadron = false;
217  unsigned int itry=0;
218 
219  while(itry++ < kRjMaxIterations && !got_charmed_hadron) {
220 
221  // Generate a charmed hadron PDG code
222  int pdg = this->GenerateCharmHadron(nu_pdg,Ev); // generate hadron
223  double mc = pdglib->Find(pdg)->Mass(); // lookup mass
224 
225  LOG("CharmHad", pNOTICE)
226  << "Trying charm hadron = " << pdg << "(m = " << mc << ")";
227 
228  if(mc>=W) continue; // dont' accept
229 
230  // Generate the charmed hadron energy based on the input
231  // fragmentation function
232  double z = fFragmFunc->GenerateZ(); // generate z(=Eh/Ev)
233  double Ec = z*Eh; // @ LAB'
234  double mc2 = TMath::Power(mc,2);
235  double Ec2 = TMath::Power(Ec,2);
236  double pc2 = Ec2-mc2;
237 
238  LOG("CharmHad", pINFO)
239  << "Trying charm hadron z = " << z << ", E = " << Ec;
240 
241  if(pc2<=0) continue;
242 
243  // Generate the charm hadron pT^2 and pL^2 (with respect to the
244  // hadronic system direction @ the LAB)
245  double ptc2 = fCharmPT2pdf->GetRandom();
246  double plc2 = Ec2 - ptc2 - mc2;
247  LOG("CharmHad", pINFO)
248  << "Trying charm hadron pT^2 (tranv to pHad) = " << ptc2;
249  if(plc2<0) continue;
250 
251  // Generate the charm hadron momentum components (@ LAB', z:\vec{pHad})
252  double ptc = TMath::Sqrt(ptc2);
253  double plc = TMath::Sqrt(plc2);
254  double phi = (2*kPi) * rnd->RndHadro().Rndm();
255  double pxc = ptc * TMath::Cos(phi);
256  double pyc = ptc * TMath::Sin(phi);
257  double pzc = plc;
258 
259  p4C.SetPxPyPzE(pxc,pyc,pzc,Ec); // @ LAB'
260 
261  // Boost charm hadron 4-momentum from the LAB' to the HCM' frame
262  //
263  LOG("CharmHad", pDEBUG)
264  << "Charm hadron p4 (@LAB') = " << utils::print::P4AsString(&p4C);
265 
266  p4C.Boost(beta);
267 
268  LOG("CharmHad", pDEBUG)
269  << "Charm hadron p4 (@HCM') = " << utils::print::P4AsString(&p4C);
270 
271  // Hadronic non-charm remnant 4p at HCM'
272  TLorentzVector p4 = p4H - p4C;
273  double wr = p4.M();
274  LOG("CharmHad", pINFO)
275  << "Invariant mass of remnant hadronic system= " << wr;
276  if(wr < kNucleonMass + kPionMass + kASmallNum) {
277  LOG("CharmHad", pINFO) << "Too small hadronic remnant mass!";
278  continue;
279  }
280 
281  ch_pdg = pdg;
282  got_charmed_hadron = true;
283 
284  LOG("CharmHad", pNOTICE)
285  << "Generated charm hadron = " << pdg << "(m = " << mc << ")";
286  LOG("CharmHad", pNOTICE)
287  << "Generated charm hadron z = " << z << ", E = " << Ec;
288  }
289 
290  // ....................................................................
291  // Check whether the code above had difficulty generating the charmed
292  // hadron near the W - threshold.
293  // If yes, attempt a phase space decay of a low mass charm hadron + nucleon
294  // pair that maintains the charge.
295  // That's a desperate solution but don't want to quit too early as that
296  // would distort the generated dsigma/dW distribution near threshold.
297  //
298  bool used_lowW_strategy = false;
299  int fs_nucleon_pdg = -1;
300  if(ch_pdg==-1 && W < 3.){
301  LOG("CharmHad", pNOTICE)
302  << "Had difficulty generating charm hadronic system near W threshold";
303  LOG("CharmHad", pNOTICE)
304  << "Trying an alternative strategy";
305 
306  double qfsl = interaction->FSPrimLepton()->Charge() / 3.;
307  double qinit = pdglib->Find(nuc_pdg)->Charge() / 3.;
308  int qhad = (int) (qinit - qfsl);
309 
310  int remn_pdg = -1;
311  int chrm_pdg = -1;
312 
313  //cc-only: qhad(nu) = +1,+2, qhad(nubar)= -1,0
314  //
315  if(qhad == 2) {
316  chrm_pdg = kPdgDP; remn_pdg = kPdgProton;
317  } else if(qhad == 1) {
318  if(rnd->RndHadro().Rndm() > 0.5) {
319  chrm_pdg = kPdgD0; remn_pdg = kPdgProton;
320  } else {
321  chrm_pdg = kPdgDP; remn_pdg = kPdgNeutron;
322  }
323  } else if(qhad == 0) {
324  chrm_pdg = kPdgAntiD0; remn_pdg = kPdgNeutron;
325  } else if(qhad == -1) {
326  chrm_pdg = kPdgDM; remn_pdg = kPdgNeutron;
327  }
328 
329  double mc = pdglib->Find(chrm_pdg)->Mass();
330  double mn = pdglib->Find(remn_pdg)->Mass();
331 
332  if(mc+mn < W) {
333  // Set decay
334  double mass[2] = {mc, mn};
335  bool permitted = fPhaseSpaceGenerator.SetDecay(p4H, 2, mass);
336  assert(permitted);
337 
338  // Get the maximum weight
339  double wmax = -1;
340  for(int i=0; i<200; i++) {
341  double w = fPhaseSpaceGenerator.Generate();
342  wmax = TMath::Max(wmax,w);
343  }
344 
345  if(wmax>0) {
346  wmax *= 2;
347 
348  // Generate unweighted decay
349  bool accept_decay=false;
350  unsigned int idecay_try=0;
351  while(!accept_decay)
352  {
353  idecay_try++;
354 
355  if(idecay_try>kMaxUnweightDecayIterations) {
356  LOG("CharmHad", pWARN)
357  << "Couldn't generate an unweighted phase space decay after "
358  << idecay_try << " attempts";
359  }
360  double w = fPhaseSpaceGenerator.Generate();
361  if(w > wmax) {
362  LOG("CharmHad", pWARN)
363  << "Decay weight = " << w << " > max decay weight = " << wmax;
364  }
365  double gw = wmax * rnd->RndHadro().Rndm();
366  accept_decay = (gw<=w);
367 
368  if(accept_decay) {
369  used_lowW_strategy = true;
370  TLorentzVector * p4 = fPhaseSpaceGenerator.GetDecay(0);
371  p4C = *p4;
372  ch_pdg = chrm_pdg;
373  fs_nucleon_pdg = remn_pdg;
374  }
375  } // decay loop
376  }//wmax>0
377 
378  }// allowed decay
379  } // alt low-W strategy
380 
381  // ....................................................................
382  // Check success in generating the charm hadron & compute 4p for
383  // remnant system
384  //
385  if(ch_pdg==-1){
386  LOG("CharmHad", pWARN)
387  << "Couldn't generate charm hadron for: " << *interaction;
388  return 0;
389  }
390 
391  TLorentzVector p4R = p4H - p4C;
392  double WR = p4R.M();
393  //double MC = pdglib->Find(ch_pdg)->Mass();
394 
395  LOG("CharmHad", pNOTICE) << "Remnant hadronic system mass = " << WR;
396 
397  // ....................................................................
398  // Handle case where the user doesn't want to remnant system to be
399  // hadronized (add as 'hadronic blob')
400  //
401  if(fCharmOnly) {
402  // Create particle list (fragmentation record)
403  TClonesArray * particle_list = new TClonesArray("genie::GHepParticle", 2);
404  particle_list->SetOwner(true);
405 
406  // insert the generated particles
407  new ((*particle_list)[0]) GHepParticle (ch_pdg,kIStStableFinalState,
408  -1,-1,-1,-1, p4C.Px(),p4C.Py(),p4C.Pz(),p4C.E(), 0,0,0,0);
409  new ((*particle_list)[1]) GHepParticle (kPdgHadronicBlob,kIStStableFinalState,
410  -1,-1,-1,-1, p4R.Px(),p4R.Py(),p4R.Pz(),p4R.E(), 0,0,0,0);
411 
412  return particle_list;
413  }
414 
415  // ....................................................................
416  // Handle case where the remnant system is already known and doesn't
417  // have to be hadronized. That happens when (close to the W threshold)
418  // the hadronic system was generated by a simple 2-body decay
419  //
420  if(used_lowW_strategy) {
421  // Create particle list (fragmentation record)
422  TClonesArray * particle_list = new TClonesArray("genie::GHepParticle", 3);
423  particle_list->SetOwner(true);
424 
425  // insert the generated particles
426  new ((*particle_list)[0]) GHepParticle (ch_pdg,kIStStableFinalState,
427  -1,-1,-1,-1, p4C.Px(),p4C.Py(),p4C.Pz(),p4C.E(), 0,0,0,0);
428  new ((*particle_list)[1]) GHepParticle (kPdgHadronicBlob,kIStNucleonTarget,
429  -1,-1,2,2, p4R.Px(),p4R.Py(),p4R.Pz(),p4R.E(), 0,0,0,0);
430  new ((*particle_list)[2]) GHepParticle (fs_nucleon_pdg,kIStStableFinalState,
431  1,1,-1,-1, p4R.Px(),p4R.Py(),p4R.Pz(),p4R.E(), 0,0,0,0);
432 
433  return particle_list;
434  }
435 
436  // ....................................................................
437  // --------------------------------------------------------------------
438  // Hadronize non-charm hadronic blob using PYTHIA/JETSET
439  // --------------------------------------------------------------------
440  // ....................................................................
441 
442  // Create output event record
443  // Insert the generated charm hadron & the hadronic (non-charm) blob.
444  // In this case the hadronic blob is entered as a pre-fragm. state.
445 
446  TClonesArray * particle_list = new TClonesArray("genie::GHepParticle");
447  particle_list->SetOwner(true);
448 
449  new ((*particle_list)[0]) GHepParticle (ch_pdg,kIStStableFinalState,
450  -1,-1,-1,-1, p4C.Px(),p4C.Py(),p4C.Pz(),p4C.E(), 0,0,0,0);
451  new ((*particle_list)[1]) GHepParticle (kPdgHadronicBlob,kIStNucleonTarget,
452  -1,-1,2,3, p4R.Px(),p4R.Py(),p4R.Pz(),p4R.E(), 0,0,0,0);
453 
454  unsigned int rpos =2; // offset in event record
455 
456  bool use_pythia = (WR>1.5);
457 
458 /*
459  // Determining quark systems to input to PYTHIA based on simple quark model
460  // arguments
461  //
462  // Neutrinos
463  // ------------------------------------------------------------------
464  // Scattering off valence q
465  // ..................................................................
466  // p: [uu]+d
467  // |--> c --> D0 <c+\bar(u)> : [u]
468  // --> D+ <c+\bar(d)> : [d]
469  // --> Ds+ <c+\bar(s)> : [s]
470  // --> Lamda_c+ <c+ud > : [\bar(ud)]
471  //
472  // (for n: [uu] -> 50%[ud]_{0} + 50%[ud]_{1})
473  //
474  // Scattering off sea q
475  // ..................................................................
476  // p: [uud] + [\bar(d)]d (or)
477  // [\bar(s)]s
478  // |--> c --> D0 <c+\bar(u)> : [u]
479  // --> D+ <c+\bar(d)> : [d]
480  // --> Ds+ <c+\bar(s)> : [s]
481  // --> Lamda_c+ <c+ud > : [\bar(ud)]
482  // Anti-Neutrinos
483  // ------------------------------------------------------------------
484  // Scattering off sea q
485  // ..................................................................
486  // p: [uud] + [d] \bar(d) (or)
487  // [s] \bar(s)
488  // |----> \bar(c) --> \bar(D0) <\bar(c)+u> : [\bar(u)]
489  // --> D- <\bar(c)+d> : [\bar(d)]
490  // --> Ds- <\bar(c)+s> : [\bar(s)]
491  // [Summary]
492  // Qq
493  // | v + p [val/d] --> D0 + { u uu }(+2) / u,uu
494  // | v + p [val/d] --> D+ + { d uu }(+1) / d,uu
495  // | v + p [val/d] --> Ds+ + { s uu }(+1) / s,uu
496  // | v + p [val/d] --> Lc+ + { \bar(ud) uu }(+1) / \bar(d),u
497  // | v + n [val/d] --> D0 + { u ud }(+1) / u,ud
498  // | v + n [val/d] --> D+ + { d ud }( 0) / d,ud
499  // | v + n [val/d] --> Ds+ + { s ud }( 0) / s,ud
500  // | v + n [val/d] --> Lc+ + { \bar(ud) ud }( 0) / \bar(d),d
501  // | v + p [sea/d] --> D0 + { uud \bar(d) u }(+2) / u,uu
502  // | v + p [sea/d] --> D+ + { uud \bar(d) d }(+1) / d,uu
503  // | v + p [sea/d] --> Ds+ + { uud \bar(d) s }(+1) / s,uu
504  // | v + p [sea/d] --> Lc+ + { uud \bar(d) \bar(ud) }(+1) / \bar(d),u
505  // | v + n [sea/d] --> D0 + { udd \bar(d) u }(+1) / u,ud
506  // | v + n [sea/d] --> D+ + { udd \bar(d) d }( 0) / d,ud
507  // | v + n [sea/d] --> Ds+ + { udd \bar(d) s }( 0) / s,ud
508  // | v + n [sea/d] --> Lc+ + { udd \bar(d) \bar(ud) }( 0) / \bar(d),d
509  // | v + p [sea/s] --> D0 + { uud \bar(s) u }(+2) / u,uu
510  // | v + p [sea/s] --> D+ + { uud \bar(s) d }(+1) / d,uu
511  // | v + p [sea/s] --> Ds+ + { uud \bar(s) s }(+1) / s,uu
512  // | v + p [sea/s] --> Lc+ + { uud \bar(s) \bar(ud) }(+1) / \bar(d),u
513  // | v + n [sea/s] --> D0 + { udd \bar(s) u }(+1) / u,ud
514  // | v + n [sea/s] --> D+ + { udd \bar(s) d }( 0) / d,ud
515  // | v + n [sea/s] --> Ds+ + { udd \bar(s) s }( 0) / s,ud
516  // | v + n [sea/s] --> Lc+ + { udd \bar(s) \bar(ud) }( 0) / \bar(d),d
517 
518  // | \bar(v) + p [sea/\bar(d)] --> \bar(D0) + { uud d \bar(u) }( 0) / d,ud
519  // | \bar(v) + p [sea/\bar(d)] --> D- + { uud d \bar(d) }(+1) / d,uu
520  // | \bar(v) + p [sea/\bar(d)] --> Ds- + { uud d \bar(s) }(+1) / d,uu
521  // | \bar(v) + n [sea/\bar(d)] --> \bar(D0) + { udd d \bar(u) }(-1) / d,dd
522  // | \bar(v) + n [sea/\bar(d)] --> D- + { udd d \bar(d) }( 0) / d,ud
523  // | \bar(v) + n [sea/\bar(d)] --> Ds- + { udd d \bar(s) }( 0) / d,ud
524  // | \bar(v) + p [sea/\bar(s)] --> \bar(D0) + { uud s \bar(u) }( 0) / d,ud
525  // | \bar(v) + p [sea/\bar(s)] --> D- + { uud s \bar(d) }(+1) / d,uu
526  // | \bar(v) + p [sea/\bar(s)] --> Ds- + { uud s \bar(s) }(+1) / d,uu
527  // | \bar(v) + n [sea/\bar(s)] --> \bar(D0) + { udd s \bar(u) }(-1) / d,dd
528  // | \bar(v) + n [sea/\bar(s)] --> D- + { udd s \bar(d) }( 0) / d,ud
529  // | \bar(v) + n [sea/\bar(s)] --> Ds- + { udd s \bar(s) }( 0) / d,ud
530  //
531  //
532  // Taking some short-cuts below :
533  // In the process of obtaining 2 q systems (while conserving the charge) I might tread
534  // d,s or \bar(d),\bar(s) as the same
535  // In the future I should perform the first steps of the multi-q hadronization manualy
536  // (to reduce the number of q's input to PYTHIA) or use py3ent_(), py4ent_() ...
537  //
538 */
539 
540  if(use_pythia) {
541  int qrkSyst1 = 0;
542  int qrkSyst2 = 0;
543  if(isnu||isdm) { // neutrinos
544  if(ch_pdg==kPdgLambdaPc) {
545  if(isp) { qrkSyst1 = kPdgAntiDQuark; qrkSyst2 = kPdgUQuark; }
546  if(isn) { qrkSyst1 = kPdgAntiDQuark; qrkSyst2 = kPdgDQuark; }
547  } else {
548  if(isp) { qrkSyst2 = kPdgUUDiquarkS1; }
549  if(isn) { qrkSyst2 = (rnd->RndHadro().Rndm()<0.5) ? kPdgUDDiquarkS0 : kPdgUDDiquarkS1; }
550  if (ch_pdg==kPdgD0 ) { qrkSyst1 = kPdgUQuark; }
551  if (ch_pdg==kPdgDP ) { qrkSyst1 = kPdgDQuark; }
552  if (ch_pdg==kPdgDPs ) { qrkSyst1 = kPdgSQuark; }
553  }
554  }
555  if(isnub) { // antineutrinos
556  qrkSyst1 = kPdgDQuark;
557  if (isp && ch_pdg==kPdgAntiD0) { qrkSyst2 = (rnd->RndHadro().Rndm()<0.5) ? kPdgUDDiquarkS0 : kPdgUDDiquarkS1; }
558  if (isp && ch_pdg==kPdgDM ) { qrkSyst2 = kPdgUUDiquarkS1; }
559  if (isp && ch_pdg==kPdgDMs ) { qrkSyst2 = kPdgUUDiquarkS1; }
560  if (isn && ch_pdg==kPdgAntiD0) { qrkSyst2 = kPdgDDDiquarkS1; }
561  if (isn && ch_pdg==kPdgDM ) { qrkSyst2 = (rnd->RndHadro().Rndm()<0.5) ? kPdgUDDiquarkS0 : kPdgUDDiquarkS1; }
562  if (isn && ch_pdg==kPdgDMs ) { qrkSyst2 = (rnd->RndHadro().Rndm()<0.5) ? kPdgUDDiquarkS0 : kPdgUDDiquarkS1; }
563  }
564  if(qrkSyst1 == 0 && qrkSyst2 == 0) {
565  LOG("CharmHad", pWARN)
566  << "Couldn't generate quark systems for PYTHIA in: " << *interaction;
567  return 0;
568  }
569 
570  //
571  // Run PYTHIA for the hadronization of remnant system
572  //
573  fPythia->SetMDCY(fPythia->Pycomp(kPdgPi0), 1,0); // don't decay pi0
574  fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaM), 1,1); // decay Delta+
575  fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_Delta0), 1,1); // decay Delta++
576  fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaP), 1,1); // decay Delta++
577  fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaPP), 1,1); // decay Delta++
578 // fPythia->SetMDCY(fPythia->Pycomp(kPdgDeltaP), 1,1); // decay Delta+
579 // fPythia->SetMDCY(fPythia->Pycomp(kPdgDeltaPP), 1,1); // decay Delta++
580 
581  int ip = 0;
582  py2ent_(&ip, &qrkSyst1, &qrkSyst2, &WR); // hadronize
583 
584  fPythia->SetMDCY(fPythia->Pycomp(kPdgPi0),1,1); // restore
585 
586  //-- Get PYTHIA's LUJETS event record
587  TClonesArray * pythia_remnants = 0;
588  fPythia->GetPrimaries();
589  pythia_remnants = dynamic_cast<TClonesArray *>(fPythia->ImportParticles("All"));
590  if(!pythia_remnants) {
591  LOG("CharmHad", pWARN) << "Couldn't hadronize (non-charm) remnants!";
592  return 0;
593  }
594 
595  int np = pythia_remnants->GetEntries();
596  assert(np>0);
597 
598  // PYTHIA performs the hadronization at the *remnant hadrons* centre of mass
599  // frame (not the hadronic centre of mass frame).
600  // Boost all hadronic blob fragments to the HCM', fix their mother/daughter
601  // assignments and add them to the fragmentation record.
602 
603  TVector3 rmnbeta = +1 * p4R.BoostVector(); // boost velocity
604 
605  TMCParticle * pythia_remn = 0; // remnant
606  GHepParticle * bremn = 0; // boosted remnant
607  TIter remn_iter(pythia_remnants);
608  while( (pythia_remn = (TMCParticle *) remn_iter.Next()) ) {
609 
610  // insert and get a pointer to inserted object for mods
611  bremn = new ((*particle_list)[rpos++]) GHepParticle ( pythia_remn->GetKF(), // pdg
612  GHepStatus_t(pythia_remn->GetKS()), // status
613  pythia_remn->GetParent(), // first parent
614  -1, // second parent
615  pythia_remn->GetFirstChild(), // first daughter
616  pythia_remn->GetLastChild(), // second daughter
617  pythia_remn -> GetPx(), // px
618  pythia_remn -> GetPy(), // py
619  pythia_remn -> GetPz(), // pz
620  pythia_remn -> GetEnergy(), // e
621  pythia_remn->GetVx(), // x
622  pythia_remn->GetVy(), // y
623  pythia_remn->GetVz(), // z
624  pythia_remn->GetTime() // t
625  );
626 
627  // boost
628  bremn -> P4() -> Boost( rmnbeta ) ;
629 
630  // handle insertion of charmed hadron
631  int jp = bremn->FirstMother();
632  int ifc = bremn->FirstDaughter();
633  int ilc = bremn->LastDaughter();
634 
635  bremn -> SetFirstMother( (jp == 0 ? 1 : jp +1) );
636  bremn -> SetFirstDaughter ( (ifc == 0 ? -1 : ifc+1) );
637  bremn -> SetLastDaughter ( (ilc == 0 ? -1 : ilc+1) );
638  }
639  } // use_pythia
640 
641  // ....................................................................
642  // Hadronizing low-W non-charm hadronic blob using a phase space decay
643  // ....................................................................
644 
645  else {
646  // Just a small fraction of events (low-W remnant syste) causing trouble
647  // to PYTHIA/JETSET
648  // Set a 2-body N+pi system that matches the remnant system charge and
649  // do a simple phase space decay
650  //
651  // q(remn) remn/syst
652  // +2 : (p pi+)
653  // +1 : 50%(p pi0) + 50% n pi+
654  // 0 : 50%(p pi-) + 50% n pi0
655  // -1 : (n pi-)
656  //
657  double qfsl = interaction->FSPrimLepton()->Charge() / 3.;
658  double qinit = pdglib->Find(nuc_pdg)->Charge() / 3.;
659  double qch = pdglib->Find(ch_pdg)->Charge() / 3.;
660  int Q = (int) (qinit - qfsl - qch); // remnant hadronic system charge
661 
662  bool allowdup=true;
663  PDGCodeList pd(allowdup);
664  if(Q==2) {
666  else if (Q==1) {
667  if(rnd->RndHadro().Rndm()<0.5) {
669  else {
671  }
672  else if (Q==0) {
673  if(rnd->RndHadro().Rndm()<0.5) {
675  else {
677  }
678  else if (Q==-1) {
680 
681  double mass[2] = {
682  pdglib->Find(pd[0])->Mass(), pdglib->Find(pd[1])->Mass()
683  };
684 
685  // Set the decay
686  bool permitted = fPhaseSpaceGenerator.SetDecay(p4R, 2, mass);
687  if(!permitted) {
688  LOG("CharmHad", pERROR) << " *** Phase space decay is not permitted";
689  return 0;
690  }
691  // Get the maximum weight
692  double wmax = -1;
693  for(int i=0; i<200; i++) {
694  double w = fPhaseSpaceGenerator.Generate();
695  wmax = TMath::Max(wmax,w);
696  }
697  if(wmax<=0) {
698  LOG("CharmHad", pERROR) << " *** Non-positive maximum weight";
699  LOG("CharmHad", pERROR) << " *** Can not generate an unweighted phase space decay";
700  return 0;
701  }
702 
703  LOG("CharmHad", pINFO)
704  << "Max phase space gen. weight @ current hadronic system: " << wmax;
705 
706  // *** generating an un-weighted decay ***
707  wmax *= 1.3;
708  bool accept_decay=false;
709  unsigned int idectry=0;
710  while(!accept_decay)
711  {
712  idectry++;
713  if(idectry>kMaxUnweightDecayIterations) {
714  // report, clean-up and return
715  LOG("Char,Had", pWARN)
716  << "Couldn't generate an unweighted phase space decay after "
717  << itry << " attempts";
718  return 0;
719  }
720  double w = fPhaseSpaceGenerator.Generate();
721  if(w > wmax) {
722  LOG("CharmHad", pWARN)
723  << "Decay weight = " << w << " > max decay weight = " << wmax;
724  }
725  double gw = wmax * rnd->RndHadro().Rndm();
726  accept_decay = (gw<=w);
727  LOG("CharmHad", pDEBUG)
728  << "Decay weight = " << w << " / R = " << gw << " - accepted: " << accept_decay;
729  }
730  for(unsigned int i=0; i<2; i++) {
731  int pdgc = pd[i];
732  TLorentzVector * p4d = fPhaseSpaceGenerator.GetDecay(i);
733  new ( (*particle_list)[rpos+i] ) GHepParticle(
734  pdgc,kIStStableFinalState,1,1,-1,-1,p4d->Px(),p4d->Py(),p4d->Pz(),p4d->Energy(),
735  0,0,0,0);
736  }
737  }
738 
739  //-- Print & return the fragmentation record
740  //utils::fragmrec::Print(particle_list);
741  return particle_list;
742 }
743 //____________________________________________________________________________
744 int AGCharm2019::GenerateCharmHadron(int nu_pdg, double EvLab) const
745 {
746  // generate a charmed hadron pdg code using a charm fraction table
747 
748  RandomGen * rnd = RandomGen::Instance();
749  double r = rnd->RndHadro().Rndm();
750 
751  // the ratios are giving up to a certain value.
752  // The rations saturates at high energies, so the values used above that enery
753  // are the evaluated at
754 
755  if(pdg::IsNeutrino(nu_pdg)) {
756 
757  // the ratios are giving up to a certain value.
758  // The rations saturates at high energies, so the values used above that enery
759  // are the evaluated at the maximum energies avaiable for the ratios
760 
761  EvLab = TMath::Min( EvLab, fFracMaxEnergy ) ;
762 
763  double tf = 0;
764  if (r < (tf+=fD0FracSpl->Evaluate(EvLab))) return kPdgD0; // D^0
765  else if (r < (tf+=fDpFracSpl->Evaluate(EvLab))) return kPdgDP; // D^+
766  else if (r < (tf+=fDsFracSpl->Evaluate(EvLab))) return kPdgDPs; // Ds^+
767  else return kPdgLambdaPc; // Lamda_c^+
768 
769  } else if(pdg::IsAntiNeutrino(nu_pdg)) {
770  if (r < fD0BarFrac) return kPdgAntiD0;
771  else if (r < fD0BarFrac+fDmFrac) return kPdgDM;
772  else return kPdgDMs;
773  }
774 
775  LOG("CharmHad", pERROR) << "Could not generate a charm hadron!";
776  return 0;
777 }
778 //____________________________________________________________________________
779 double AGCharm2019::Weight(void) const
780 {
781  return 1. ;
782 }
783 
784 //____________________________________________________________________________
785 void AGCharm2019::Configure(const Registry & config)
786 {
787  Algorithm::Configure(config);
788  this->LoadConfig();
789 }
790 //____________________________________________________________________________
791 void AGCharm2019::Configure(string config)
792 {
793  Algorithm::Configure(config);
794  this->LoadConfig();
795 }
796 //____________________________________________________________________________
798 {
799 
800  bool hadronize_remnants ;
801  GetParamDef( "HadronizeRemnants", hadronize_remnants, true ) ;
802 
803  fCharmOnly = ! hadronize_remnants ;
804 
805  //-- Get a fragmentation function
806  fFragmFunc = dynamic_cast<const FragmentationFunctionI *> (
807  this->SubAlg("FragmentationFunc"));
808  assert(fFragmFunc);
809 
810  string pt_function ;
811  this -> GetParam( "PTFunction", pt_function ) ;
812 
813  fCharmPT2pdf = new TF1("fCharmPT2pdf", pt_function.c_str(),0,0.6);
814 
815  // stop ROOT from deleting this object of its own volition
816  gROOT->GetListOfFunctions()->Remove(fCharmPT2pdf);
817 
818  // neutrino charm fractions: D^0, D^+, Ds^+ (remainder: Lamda_c^+)
819  std::vector<double> ec, d0frac, dpfrac, dsfrac ;
820 
821  std::string raw ;
822  std::vector<std::string> bits ;
823 
824  bool invalid_configuration = false ;
825 
826  // load energy points
827  this -> GetParamVect( "CharmFrac-E", ec ) ;
828  fFracMaxEnergy = ec.back() - controls::kASmallNum ;
829 
830  // load D0 fractions
831  this -> GetParamVect( "CharmFrac-D0", d0frac ) ;
832 
833  // check the size
834  if ( d0frac.size() != ec.size() ) {
835  LOG("AGCharm2019", pFATAL) << "E entries don't match D0 fraction entries";
836  LOG("AGCharm2019", pFATAL) << "E: " << ec.size() ;
837  LOG("AGCharm2019", pFATAL) << "D0: " << d0frac.size() ;
838  invalid_configuration = true ;
839  }
840 
841  // load D+ fractions
842  this -> GetParamVect( "CharmFrac-D+", dpfrac ) ;
843 
844  // check the size
845  if ( dpfrac.size() != ec.size() ) {
846  LOG("AGCharm2019", pFATAL) << "E entries don't match D+ fraction entries";
847  LOG("AGCharm2019", pFATAL) << "E: " << ec.size() ;
848  LOG("AGCharm2019", pFATAL) << "D+: " << dpfrac.size() ;
849  invalid_configuration = true ;
850  }
851 
852  // load D_s fractions
853  this -> GetParamVect( "CharmFrac-Ds", dsfrac ) ;
854 
855  // check the size
856  if ( dsfrac.size() != ec.size() ) {
857  LOG("AGCharm2019", pFATAL) << "E entries don't match Ds fraction entries";
858  LOG("AGCharm2019", pFATAL) << "E: " << ec.size() ;
859  LOG("AGCharm2019", pFATAL) << "Ds: " << dsfrac.size() ;
860  invalid_configuration = true ;
861  }
862 
863  fD0FracSpl = new Spline( ec.size(), & ec[0], & d0frac[0] );
864  fDpFracSpl = new Spline( ec.size(), & ec[0], & dpfrac[0] );
865  fDsFracSpl = new Spline( ec.size(), & ec[0], & dsfrac[0] );
866 
867  // anti-neutrino charm fractions: bar(D^0), D^-, (remainder: Ds^-)
868 
869  this -> GetParam( "CharmFrac-D0bar", fD0BarFrac ) ;
870  this -> GetParam( "CharmFrac-D-", fDmFrac ) ;
871 
872  if ( invalid_configuration ) {
873 
874  LOG("AGCharm2019", pFATAL)
875  << "Invalid configuration: Exiting" ;
876 
877  // From the FreeBSD Library Functions Manual
878  //
879  // EX_CONFIG (78) Something was found in an unconfigured or miscon-
880  // figured state.
881 
882  exit( 78 ) ;
883  }
884 }
885 //____________________________________________________________________________
const int kPdgAntiD0
Definition: PDGCodes.h:184
const int kPdgP33m1232_DeltaPP
Definition: PDGCodes.h:107
const int kPdgDPs
Definition: PDGCodes.h:185
const int kPdgUUDiquarkS1
Definition: PDGCodes.h:58
double W(bool selected=false) const
Definition: Kinematics.cxx:157
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
#define pERROR
Definition: Messenger.h:59
static const double kNucleonMass
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
int HitNucPdg(void) const
Definition: Target.cxx:304
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
Spline * fDpFracSpl
nu charm fraction vs Ev: D+
Definition: AGCharm2019.h:82
int FirstDaughter(void) const
Definition: GHepParticle.h:68
const int kPdgHadronicBlob
Definition: PDGCodes.h:211
void Initialize(void) const
Definition: AGCharm2019.cxx:82
int GenerateCharmHadron(int nupdg, double EvLab) const
A numeric analysis tool class for interpolating 1-D functions.
Definition: Spline.h:58
#define pFATAL
Definition: Messenger.h:56
bool IsNucleus(void) const
Definition: Target.cxx:272
const int kPdgUQuark
Definition: PDGCodes.h:42
void Configure(const Registry &config)
int GetParamVect(const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
Handle to load vectors of parameters.
TF1 * fCharmPT2pdf
charm hadron pT^2 pdf
Definition: AGCharm2019.h:76
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
TGenPhaseSpace fPhaseSpaceGenerator
a phase space generator
Definition: AGCharm2019.h:71
string P4AsString(const TLorentzVector *p)
Definition: PrintUtils.cxx:27
static const unsigned int kMaxUnweightDecayIterations
Definition: Controls.h:61
bool IsDarkMatter(int pdgc)
Definition: PDGUtils.cxx:127
bool IsChargedLepton(int pdgc)
Definition: PDGUtils.cxx:101
virtual double Weight(void) const
Definition: GHepRecord.h:124
const TLorentzVector & HadSystP4(void) const
Definition: Kinematics.h:66
double Evaluate(double x) const
Definition: Spline.cxx:363
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
virtual ~AGCharm2019()
Definition: AGCharm2019.cxx:67
const int kPdgSQuark
Definition: PDGCodes.h:46
Spline * fD0FracSpl
nu charm fraction vs Ev: D0
Definition: AGCharm2019.h:81
A list of PDG codes.
Definition: PDGCodeList.h:32
const int kPdgP33m1232_DeltaP
Definition: PDGCodes.h:106
bool fCharmOnly
don&#39;t hadronize non-charm blob
Definition: AGCharm2019.h:75
const int kPdgP33m1232_DeltaM
Definition: PDGCodes.h:104
double fFracMaxEnergy
Maximum energy available for the Meson fractions.
Definition: AGCharm2019.h:79
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
int FirstMother(void) const
Definition: GHepParticle.h:66
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
Summary information for an interaction.
Definition: Interaction.h:56
TClonesArray * Hadronize(const Interaction *) const
int LastDaughter(void) const
Definition: GHepParticle.h:69
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
const int kPdgUDDiquarkS1
Definition: PDGCodes.h:57
double fD0BarFrac
nubar {D0} charm fraction
Definition: AGCharm2019.h:84
const int kPdgGamma
Definition: PDGCodes.h:189
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:118
const int kPdgAntiDQuark
Definition: PDGCodes.h:45
void ProcessEventRecord(GHepRecord *event) const
Definition: AGCharm2019.cxx:87
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
int ProbePdg(void) const
Definition: InitialState.h:64
const int kPdgLambdaPc
Definition: PDGCodes.h:99
void py2ent_(int *, int *, int *, double *)
const int kPdgPiP
Definition: PDGCodes.h:158
const int kPdgPi0
Definition: PDGCodes.h:160
const int kPdgDQuark
Definition: PDGCodes.h:44
static const double kASmallNum
Definition: Controls.h:40
#define pINFO
Definition: Messenger.h:62
double Weight(void) const
const int kPdgUDDiquarkS0
Definition: PDGCodes.h:56
#define pWARN
Definition: Messenger.h:60
const int kPdgP33m1232_Delta0
Definition: PDGCodes.h:105
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
Pure abstract base class. Defines the FragmentationFunctionI interface to be implemented by any algor...
bool IsNeutralLepton(int pdgc)
Definition: PDGUtils.cxx:95
void py1ent_(int *, int *, double *, double *, double *)
const int kPdgDP
Definition: PDGCodes.h:181
void LoadConfig(void)
TRandom3 & RndHadro(void) const
rnd number generator used by hadronization models
Definition: RandomGen.h:53
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
Singleton class to load &amp; serve a TDatabasePDG.
Definition: PDGLibrary.h:35
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
const int kPdgDMs
Definition: PDGCodes.h:186
const FragmentationFunctionI * fFragmFunc
charm hadron fragmentation func
Definition: AGCharm2019.h:77
const int kPdgDM
Definition: PDGCodes.h:182
virtual double GenerateZ(void) const =0
const int kPdgPiM
Definition: PDGCodes.h:159
const InitialState & InitState(void) const
Definition: Interaction.h:69
const int kPdgProton
Definition: PDGCodes.h:81
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
#define pNOTICE
Definition: Messenger.h:61
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:66
static const unsigned int kRjMaxIterations
Definition: Controls.h:26
double ProbeE(RefFrame_t rf) const
const int kPdgNeutron
Definition: PDGCodes.h:83
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
const int kPdgD0
Definition: PDGCodes.h:183
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
TPythia6 * fPythia
remnant (non-charm) hadronizer
Definition: AGCharm2019.h:86
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:58
double fDmFrac
nubar D- charm fraction
Definition: AGCharm2019.h:85
enum genie::EGHepStatus GHepStatus_t
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63
Spline * fDsFracSpl
nu charm fraction vs Ev: Ds+
Definition: AGCharm2019.h:83
const int kPdgDDDiquarkS1
Definition: PDGCodes.h:55
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345