17 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,15,6) 18 #include <TMCParticle.h> 20 #include <TMCParticle6.h> 43 using namespace genie;
47 extern "C" void py1ent_(
int *,
int *,
double *,
double *,
double *);
48 extern "C" void py2ent_(
int *,
int *,
int *,
double *);
86 LOG(
"CharmHad",
pNOTICE) <<
"** Running CHARM hadronizer";
94 const InitialState & init_state = interaction -> InitState();
98 const TLorentzVector & p4Had = kinematics.
HadSystP4();
101 double W = kinematics.
W(
true);
103 TVector3
beta = -1 * p4Had.BoostVector();
104 TLorentzVector p4H(0,0,0,W);
106 double Eh = p4Had.Energy();
108 LOG(
"CharmHad",
pNOTICE) <<
"Ehad (LAB) = " << Eh <<
", W = " <<
W;
124 TLorentzVector p4C(0,0,0,0);
127 bool got_charmed_hadron =
false;
134 double mc = pdglib->
Find(pdg)->Mass();
137 <<
"Trying charm hadron = " << pdg <<
"(m = " << mc <<
")";
145 double mc2 = TMath::Power(mc,2);
146 double Ec2 = TMath::Power(Ec,2);
147 double pc2 = Ec2-mc2;
150 <<
"Trying charm hadron z = " << z <<
", E = " << Ec;
157 double plc2 = Ec2 - ptc2 - mc2;
159 <<
"Trying charm hadron pT^2 (tranv to pHad) = " << ptc2;
166 double pxc = ptc * TMath::Cos(phi);
167 double pyc = ptc * TMath::Sin(phi);
170 p4C.SetPxPyPzE(pxc,pyc,pzc,Ec);
183 TLorentzVector
p4 = p4H - p4C;
186 <<
"Invariant mass of remnant hadronic system= " << wr;
188 LOG(
"CharmHad",
pINFO) <<
"Too small hadronic remnant mass!";
193 got_charmed_hadron =
true;
196 <<
"Generated charm hadron = " << pdg <<
"(m = " << mc <<
")";
198 <<
"Generated charm hadron z = " << z <<
", E = " << Ec;
209 bool used_lowW_strategy =
false;
210 int fs_nucleon_pdg = -1;
211 if(ch_pdg==-1 && W < 3.){
213 <<
"Had difficulty generating charm hadronic system near W threshold";
215 <<
"Trying an alternative strategy";
217 double qfsl = interaction->
FSPrimLepton()->Charge() / 3.;
218 double qinit = pdglib->
Find(nuc_pdg)->Charge() / 3.;
219 int qhad = (
int) (qinit - qfsl);
228 }
else if(qhad == 1) {
234 }
else if(qhad == 0) {
236 }
else if(qhad == -1) {
240 double mc = pdglib->
Find(chrm_pdg)->Mass();
241 double mn = pdglib->
Find(remn_pdg)->Mass();
245 double mass[2] = {
mc, mn};
251 for(
int i=0;
i<200;
i++) {
253 wmax = TMath::Max(wmax,w);
260 bool accept_decay=
false;
261 unsigned int idecay_try=0;
268 <<
"Couldn't generate an unweighted phase space decay after " 269 << idecay_try <<
" attempts";
274 <<
"Decay weight = " << w <<
" > max decay weight = " << wmax;
276 double gw = wmax * rnd->
RndHadro().Rndm();
277 accept_decay = (gw<=
w);
280 used_lowW_strategy =
true;
284 fs_nucleon_pdg = remn_pdg;
298 <<
"Couldn't generate charm hadron for: " << *interaction;
302 TLorentzVector p4R = p4H - p4C;
304 double MC = pdglib->
Find(ch_pdg)->Mass();
306 LOG(
"CharmHad",
pNOTICE) <<
"Remnant hadronic system mass = " << WR;
314 TClonesArray * particle_list =
new TClonesArray(
"TMCParticle", 2);
315 particle_list->SetOwner(
true);
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);
321 -1,-1,-1, p4R.Px(),p4R.Py(),p4R.Pz(),p4R.E(),WR, 0,0,0,0,0);
323 return particle_list;
331 if(used_lowW_strategy) {
333 TClonesArray * particle_list =
new TClonesArray(
"TMCParticle", 3);
334 particle_list->SetOwner(
true);
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);
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);
344 return particle_list;
357 TClonesArray * particle_list =
new TClonesArray(
"TMCParticle");
358 particle_list->SetOwner(
true);
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);
363 -1,2,3, p4R.Px(),p4R.Py(),p4R.Pz(),p4R.E(),WR, 0,0,0,0,0);
365 unsigned int rpos =2;
367 bool use_pythia = (WR>1.5);
475 if(qrkSyst1 == 0 && qrkSyst2 == 0) {
477 <<
"Couldn't generate quark systems for PYTHIA in: " << *interaction;
493 py2ent_(&ip, &qrkSyst1, &qrkSyst2, &WR);
498 TClonesArray * remnants = 0;
500 remnants =
dynamic_cast<TClonesArray *
>(
fPythia->ImportParticles(
"All"));
502 LOG(
"CharmHad",
pWARN) <<
"Couldn't hadronize (non-charm) remnants!";
511 TVector3 rmnbeta = +1 * p4R.BoostVector();
513 TMCParticle * remn = 0;
514 TMCParticle * bremn = 0;
515 TIter remn_iter(remnants);
516 while( (remn = (TMCParticle *) remn_iter.Next()) ) {
519 bremn =
new ((*particle_list)[rpos++]) TMCParticle (*remn);
522 TLorentzVector
p4(remn->GetPx(),remn->GetPy(),remn->GetPz(),remn->GetEnergy());
524 bremn -> SetPx (
p4.Px());
525 bremn -> SetPy (
p4.Py());
526 bremn -> SetPz (
p4.Pz());
527 bremn -> SetEnergy (
p4.E() );
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) );
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);
580 pdglib->
Find(pd[0])->Mass(), pdglib->
Find(pd[1])->Mass()
586 LOG(
"CharmHad",
pERROR) <<
" *** Phase space decay is not permitted";
591 for(
int i=0;
i<200;
i++) {
593 wmax = TMath::Max(wmax,w);
596 LOG(
"CharmHad",
pERROR) <<
" *** Non-positive maximum weight";
597 LOG(
"CharmHad",
pERROR) <<
" *** Can not generate an unweighted phase space decay";
602 <<
"Max phase space gen. weight @ current hadronic system: " << wmax;
606 bool accept_decay=
false;
607 unsigned int idectry=0;
614 <<
"Couldn't generate an unweighted phase space decay after " 615 << itry <<
" attempts";
621 <<
"Decay weight = " << w <<
" > max decay weight = " << wmax;
623 double gw = wmax * rnd->
RndHadro().Rndm();
624 accept_decay = (gw<=
w);
626 <<
"Decay weight = " << w <<
" / R = " << gw <<
" - accepted: " << accept_decay;
628 for(
unsigned int i=0;
i<2;
i++) {
631 new ( (*particle_list)[rpos+
i] ) TMCParticle(
632 1,pdgc,1,-1,-1,p4d->Px(),p4d->Py(),p4d->Pz(),p4d->Energy(),
639 return particle_list;
679 LOG(
"CharmHad",
pERROR) <<
"Could not generate a charm hadron!";
698 bool hadronize_remnants = true ;
699 GetParam(
"HadronizeRemnants", hadronize_remnants,
false ) ;
705 this->
SubAlg(
"FragmentationFunc"));
708 fCharmPT2pdf =
new TF1(
"fCharmPT2pdf",
"exp(-0.213362-6.62464*x)",0,0.6);
715 double ec[
nc] = {0.,5.,10.,15.,20.,25.,30.,35.,40.,50.,60.,70.,80.,90.,100.};
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 };
const int kPdgP33m1232_DeltaPP
const int kPdgUUDiquarkS1
double W(bool selected=false) const
double fDmFrac
nubar D- charm fraction
TPythia6 * fPythia
remnant (non-charm) hadronizer
bool IsNeutrino(int pdgc)
THE MAIN GENIE PROJECT NAMESPACE
Spline * fDpFracSpl
nu charm fraction vs Ev: D+
static const double kNucleonMass
static RandomGen * Instance()
Access instance.
int HitNucPdg(void) const
virtual double GenerateZ(void) const =0
const int kPdgHadronicBlob
virtual ~CharmHadronization()
A numeric analysis tool class for interpolating 1-D functions.
void Initialize(void) const
define the HadronizationModelI interface
void py1ent_(int *, int *, double *, double *, double *)
Generated/set kinematical variables for an event.
string P4AsString(const TLorentzVector *p)
static const unsigned int kMaxUnweightDecayIterations
bool IsDarkMatter(int pdgc)
const TLorentzVector & HadSystP4(void) const
double Evaluate(double x) const
A singleton holding random number generator classes. All random number generation in GENIE should tak...
PDGCodeList * SelectParticles(const Interaction *) const
void py2ent_(int *, int *, int *, double *)
const int kPdgP33m1232_DeltaP
const int kPdgP33m1232_DeltaM
TClonesArray * Hadronize(const Interaction *) const
Summary information for an interaction.
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
const int kPdgUDDiquarkS1
bool IsAntiNeutrino(int pdgc)
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
const FragmentationFunctionI * fFragmFunc
charm hadron fragmentation func
virtual void Configure(const Registry &config)
double fD0BarFrac
nubar {D0} charm fraction
void Configure(const Registry &config)
static const double kASmallNum
const int kPdgUDDiquarkS0
Misc GENIE control constants.
const int kPdgP33m1232_Delta0
bool fCharmOnly
don'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...
TRandom3 & RndHadro(void) const
rnd number generator used by hadronization models
static PDGLibrary * Instance(void)
static const double kPionMass
Singleton class to load & serve a TDatabasePDG.
Var Sqrt(const Var &v)
Use to take sqrt of a var.
int GenerateCharmHadron(int nupdg, double EvLab) const
A registry. Provides the container for algorithm configuration parameters.
Pure abstract base class. Defines the HadronizationModelI interface to be implemented by any algorith...
TParticlePDG * Find(int pdgc)
double Weight(void) const
assert(nhit_max >=nhit_nbins)
Spline * fDsFracSpl
nu charm fraction vs Ev: Ds+
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
static const unsigned int kRjMaxIterations
double ProbeE(RefFrame_t rf) const
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)
TGenPhaseSpace fPhaseSpaceGenerator
a phase space generator
Initial State information.
TH1D * MultiplicityProb(const Interaction *, Option_t *opt="") const
const int kPdgDDDiquarkS1
const Algorithm * SubAlg(const RgKey ®istry_key) const