Commit 6aa9e2c9 authored by Frank Krauss's avatar Frank Krauss

Seems to work fine.

parent e9dad575
......@@ -63,9 +63,9 @@ Cluster_Decay_Analysis::~Cluster_Decay_Analysis()
void Cluster_Decay_Analysis::AnalyseThis(Blob * blob)
{
return;
int Npiplus=0,Npiminus=0,Npi0=0,NKplus=0,NKminus=0,NK0=0,NK0b=0,Neta=0,Netaprime=0,NPS=0;
int Nrhoplus=0,Nrhominus=0,Nrho0=0,NKstarplus=0,NKstarminus=0,NKstar0=0,NKstar0b=0,
Nomega=0,Nphi=0,NV=0;
int Npiplus(0),Npiminus(0),Npi0(0),NKplus(0),NKminus(0),NK0(0),NK0b(0),Neta(0),Netaprime(0),
Nrhoplus(0),Nrhominus(0),Nrho0(0),NKstarplus(0),NKstarminus(0),NKstar0(0),NKstar0b(0),
Nomega(0),Nphi(0);
Particle * part;
int kfc, LambdaCount(0), LambdaP(0), LambdaM(0);
Vec4D QLambda(0.,0.,0.,0.);
......@@ -83,7 +83,7 @@ void Cluster_Decay_Analysis::AnalyseThis(Blob * blob)
break;
}
}
if (max_x=max_xB) m_histograms[string("x_E_B-Quarks_L")]->Insert(max_x);
if (max_x==max_xB) m_histograms[string("x_E_B-Quarks_L")]->Insert(max_x);
max_x = max_xB = 0.;
......@@ -202,7 +202,7 @@ void Cluster_Decay_Analysis::AnalyseThis(Blob * blob)
}
if (LambdaCount==2 && LambdaP==1 && LambdaM==1)
m_histograms[string("Q(Lambda_Lambda)")]->Insert(sqrt(-QLambda.Abs2()));
if (max_x=max_xB) m_histograms[string("x_E_B-Mesons_L")]->Insert(max_x);
if (max_x==max_xB) m_histograms[string("x_E_B-Mesons_L")]->Insert(max_x);
m_histograms[string("pi+_Number")]->Insert(Npiplus);
m_histograms[string("pi-_Number")]->Insert(Npiminus);
......
......@@ -8,9 +8,9 @@ using namespace ATOOLS;
using namespace std;
Cluster_Part::Cluster_Part(Dipole_Splitter * splitter,bool ana) :
m_ana(ana),
m_pt2max_factor(sqr(hadpars.Get(std::string("ptmax_factor")))),
p_splitter(splitter),
m_ana(ana)
p_splitter(splitter)
{
if (m_ana) {
m_histograms[string("PT_Cluster")] = new Histogram(0,0.,1.5,150);
......
......@@ -33,9 +33,6 @@ Cluster_Formation_Handler::Cluster_Formation_Handler(Cluster_List* clulist,
m_histograms[string("Cluster_Mass_Transformed")] = new Histogram(0,0.,100.,200);
m_histograms[string("Cluster_Number_Formation")] = new Histogram(0,0.,20.,20);
m_histograms[string("Cluster_Number_Transformed")] = new Histogram(0,0.,20.,20);
m_histograms[string("XB_bquark_before_shift")] = new Histogram(0,0.,1.,100);
m_histograms[string("XB_bquark_after_shift")] = new Histogram(0,0.,1.,100);
m_histograms[string("XB_bquark_in_cluster")] = new Histogram(0,0.,1.,100);
}
}
......@@ -74,41 +71,10 @@ int Cluster_Formation_Handler::FormClusters(Blob * blob) {
<<"In "<<METHOD<<": hadronize "<<blob->NInP()<<" partons."<<std::endl
<<(*blob)<<std::endl;
if (!ExtractSinglets(blob)) { return -1; }
Histogram * histoxb;
LPPL_Iterator pplit;
PPL_Iterator pit;
for(pplit=m_partlists.begin(); pplit!=m_partlists.end(); ++pplit) {
for(pit=(*pplit)->begin(); pit!=(*pplit)->end(); ++pit) {
if (m_analyse && (*pit)->m_flav.Kfcode()==5) {
histoxb = (m_histograms.find(string("XB_bquark_before_shift")))->second;
histoxb->Insert((*pit)->m_mom[0]/45.625);
}
}
}
if (!ShiftOnMassShells()) { return -1; }
for(pplit=m_partlists.begin(); pplit!=m_partlists.end(); ++pplit) {
for(pit=(*pplit)->begin(); pit!=(*pplit)->end(); ++pit) {
if (m_analyse && (*pit)->m_flav.Kfcode()==5) {
histoxb = (m_histograms.find(string("XB_bquark_after_shift")))->second;
histoxb->Insert((*pit)->m_mom[0]/45.625);
}
}
}
if (!FormOriginalClusters()) { return -1; }
if (!ApplyColourReconnections()) { return 0; }
if (!MergeClusterListsIntoOne()) { return 0; }
for(Cluster_Iterator cit=p_clulist->begin(); cit!=p_clulist->end(); cit++) {
if (m_analyse && (*cit)->GetTrip()->m_flav.Kfcode()==5) {
histoxb = (m_histograms.find(string("XB_bquark_in_cluster")))->second;
histoxb->Insert((*cit)->GetTrip()->m_mom[0]/45.625);
msg_Debugging()<<"Clusters (b) : "<<(*cit)->GetTrip()->m_mom<<"."<<std::endl;
}
if (m_analyse && (*cit)->GetAnti()->m_flav.Kfcode()==5) {
histoxb = (m_histograms.find(string("XB_bquark_in_cluster")))->second;
histoxb->Insert((*cit)->GetAnti()->m_mom[0]/45.625);
msg_Debugging()<<"Clusters (bbar) : "<<(*cit)->GetAnti()->m_mom<<"."<<std::endl;
}
}
if (!ClustersToHadrons(blob)) { return -1; }
return 1;
......@@ -143,9 +109,10 @@ void Cluster_Formation_Handler::Reset() {
bool Cluster_Formation_Handler::ExtractSinglets(Blob * blob)
{
Proto_Particle_List * pli(NULL);
bool construct(false);
int col1, col2;
Particle * part;
bool construct(false);
unsigned int col1(0), col2(0);
leading::code leading(hadpars.GetSplitter()->Leading());
Particle * part(NULL);
for (int i=0;i<blob->NInP();i++) {
part = blob->InParticle(i);
if ((part->Status()!=part_status::active && part->Status()!=part_status::fragmented) ||
......@@ -153,6 +120,7 @@ bool Cluster_Formation_Handler::ExtractSinglets(Blob * blob)
if (construct) {
if (part->GetFlow(2)==col1) {
Proto_Particle * copy = new Proto_Particle(part->Flav(),part->Momentum(),'L');
SetInfoTagForPrimaryParticle(copy,leading);
#ifdef memchecker
std::cout<<"### New Proto_Particle ("
<<copy<<"/"<<part->Flav()<<") from "<<METHOD<<"."<<std::endl;
......@@ -175,6 +143,7 @@ bool Cluster_Formation_Handler::ExtractSinglets(Blob * blob)
col2 = part->GetFlow(2);
pli = new Proto_Particle_List;
Proto_Particle * copy = new Proto_Particle(part->Flav(),part->Momentum(),'L');
SetInfoTagForPrimaryParticle(copy,leading);
#ifdef memchecker
std::cout<<"### New Proto_Particle ("
<<copy<<"/"<<part->Flav()<<") from "<<METHOD<<"."<<std::endl;
......@@ -187,6 +156,25 @@ bool Cluster_Formation_Handler::ExtractSinglets(Blob * blob)
return true;
}
void Cluster_Formation_Handler::SetInfoTagForPrimaryParticle(Proto_Particle * proto,
leading::code lead) const {
switch (lead) {
case leading::quarks_and_gluons:
case leading::quarks_and_gluons2:
break;
case leading::only_quarks:
if (proto->m_flav.IsQuark() || proto->m_flav.IsDiQuark())
proto->m_info='L';
else
proto->m_info='l';
break;
case leading::none:
default:
proto->m_info='l';
break;
}
}
bool Cluster_Formation_Handler::ShiftOnMassShells() {
ListOfPPLs shiftables, nonshiftables;
LPPL_Iterator pplit;
......@@ -330,8 +318,7 @@ bool Cluster_Formation_Handler::FormOriginalClusters()
while (!m_partlists.empty()) {
pplit=m_partlists.begin();
msg_Tracking()<<"======= "<<METHOD<<" for :"<<std::endl
<<(**pplit)<<std::endl;
msg_Tracking()<<"======= "<<METHOD<<" for :"<<std::endl<<(**pplit)<<std::endl;
if(p_gludecayer->DecayList(*pplit)) {
clist = new Cluster_List;
p_cformer->ConstructClusters(*pplit,clist);
......@@ -443,8 +430,8 @@ bool Cluster_Formation_Handler::ClustersToHadrons(Blob * blob)
<<checkbef<<" --> "<<checkaft<<"."<<std::endl;
#endif
Histogram * histomass, * histonumb;
if (m_analyse) {
Histogram * histomass, * histonumb;
histomass = (m_histograms.find(string("Cluster_Mass_Transformed")))->second;
histonumb = (m_histograms.find(string("Cluster_Number_Transformed")))->second;
int numb = p_clulist->size();
......
......@@ -32,11 +32,11 @@ namespace AHADIC {
bool ShiftList(Proto_Particle_List *);
Proto_Particle_List * SelectFromList(ListOfPPLs * lppl,Proto_Particle_List * ppl=NULL);
bool ExtractSinglets(ATOOLS::Blob *);
void SetInfoTagForPrimaryParticle(Proto_Particle * proto,leading::code lead) const;
bool FormOriginalClusters();
bool ApplyColourReconnections();
bool ClustersToHadrons(ATOOLS::Blob *);
bool MergeClusterListsIntoOne();
public:
Cluster_Formation_Handler(Cluster_List*, bool=false);
~Cluster_Formation_Handler();
......
......@@ -12,7 +12,6 @@ Cluster_Former::~Cluster_Former() { }
void Cluster_Former::ConstructClusters(Proto_Particle_List * plin, Cluster_List * clout)
{
Cluster * cluster(NULL);
int lead(0);
PPL_Iterator pit1,pit2;
while (!plin->empty()) {
pit1 = plin->begin();
......
......@@ -22,7 +22,6 @@ void Colour_Reconnections::Singlet_CR(Cluster_List * clin)
bool direction = (ran.Get()>0.5);
if (direction) clin->reverse();
int lead1,lead2;
cit1 = cit2 = clin->begin(); cit2++;
do {
if (TestClusters((*cit1),(*cit2),gen)) {
......
......@@ -10,7 +10,7 @@ using namespace ATOOLS;
Gluon_Decayer::Gluon_Decayer(Dipole_Splitter * splitter,bool ana) :
p_splitter(splitter),
m_pt2max_factor(sqr(hadpars.Get(std::string("ptmax_factor")))),
m_analyse(ana)
m_analyse(ana),m_tot(0),m_s(0),m_u(0),m_d(0)
{
double norm(0.);
for (FlavCCMap_Iterator fdit=hadpars.GetConstituents()->CCMap.begin();
......@@ -63,6 +63,10 @@ Gluon_Decayer::~Gluon_Decayer() {
}
m_histograms.clear();
}
msg_Out()<<"Flavour production in "<<m_tot<<" gluon splittings: "
<<"s_rate = "<<double(m_s)/double(m_tot)<<", "
<<"u_rate = "<<double(m_u)/double(m_tot)<<", "
<<"d_rate = "<<double(m_d)/double(m_tot)<<"."<<std::endl;
for (FDIter fdit=m_options.begin();fdit!=m_options.end();++fdit) {
if (fdit->second!=NULL) { delete fdit->second; fdit->second=NULL; }
......@@ -189,7 +193,7 @@ bool Gluon_Decayer::DecayDipoles() {
msg_Out()<<"--- "<<METHOD<<" splits the following dipole :"<<(*dipiter)<<std::endl;
(*dipiter)->Output();
}
if (!p_splitter->SplitDipole((*dipiter),PT2Max((*dipiter)),true)) {
if (!p_splitter->SplitDipole((*dipiter),PT2Max((*dipiter)),false)) {
switch (Rescue(dipiter)) {
case -1:
msg_Tracking()<<"--- Rescue failed."<<std::endl;
......@@ -254,7 +258,7 @@ int Gluon_Decayer::Rescue(DipIter & dip) {
(*dip)->AntiTriplet()->m_flav.IsGluon()) {
if ((*dip)!=(*m_dipoles.rbegin())) {
partner++;
if (p_splitter->SplitDipole((*partner),PT2Max((*partner)),true)) {
if (p_splitter->SplitDipole((*partner),PT2Max((*partner)),false)) {
AfterSplit(partner);
dip = partner;
return 1;
......@@ -263,7 +267,7 @@ int Gluon_Decayer::Rescue(DipIter & dip) {
if (dip!=m_dipoles.begin()) {
partner = dip;
partner--;
if (p_splitter->SplitDipole((*partner),PT2Max((*partner)),true)) {
if (p_splitter->SplitDipole((*partner),PT2Max((*partner)),false)) {
AfterSplit(partner);
dip = partner;
return 1;
......@@ -284,7 +288,7 @@ int Gluon_Decayer::Rescue(DipIter & dip) {
else if (!(*dip)->Triplet()->m_flav.IsGluon() &&
(*dip)->AntiTriplet()->m_flav.IsGluon()) {
partner++;
if (p_splitter->SplitDipole((*partner),PT2Max((*partner)),true)) {
if (p_splitter->SplitDipole((*partner),PT2Max((*partner)),false)) {
AfterSplit(partner);
dip = partner;
return 1;
......@@ -294,7 +298,7 @@ int Gluon_Decayer::Rescue(DipIter & dip) {
else if (!(*dip)->AntiTriplet()->m_flav.IsGluon() &&
(*dip)->Triplet()->m_flav.IsGluon()) {
partner--;
if (p_splitter->SplitDipole((*partner),PT2Max((*partner)),true)) {
if (p_splitter->SplitDipole((*partner),PT2Max((*partner)),false)) {
AfterSplit(partner);
dip = partner;
return 1;
......@@ -406,11 +410,16 @@ void Gluon_Decayer::AfterSplit(DipIter dip) {
}
void Gluon_Decayer::SplitIt(DipIter dipiter,Vec4D checkbef) {
msg_Tracking()<<"--- "<<METHOD<<"(in): two new proto_particles:"<<std::endl;
//msg_Out()<<"--- "<<METHOD<<"(in): two new proto_particles:"<<std::endl;
Proto_Particle * new1, * new2;
p_splitter->GetNewParticles(new1,new2);
msg_Tracking()<<" "<<(*new1)<<std::endl
<<" "<<(*new2)<<std::endl;
//msg_Out()<<" new1: "<<(*new1)<<std::endl<<" new2: "<<(*new2)<<std::endl;
if (new2->m_flav==Flavour(kf_d)) m_d++;
else if (new2->m_flav==Flavour(kf_u)) m_u++;
else if (new2->m_flav==Flavour(kf_s)) m_s++;
m_tot++;
DipIter partner;
Dipole * dip((*dipiter));
#ifdef AHAmomcheck
......
......@@ -31,6 +31,8 @@ namespace AHADIC {
void MergeDipoles(DipIter & dip1,DipIter & dip2);
void SplitIt(DipIter dip,ATOOLS::Vec4D=ATOOLS::Vec4D(0.,0.,0.,0.));
double PT2Max(Dipole * dip) const;
long int m_tot, m_s, m_u, m_d, m_di0, m_di1;
public:
Gluon_Decayer(Dipole_Splitter *,bool=false);
~Gluon_Decayer();
......
......@@ -18,11 +18,11 @@ namespace AHADIC {
class Proto_Particle {
static std::list<Proto_Particle * > s_actives;
public:
Proto_Particle * p_partner;
ATOOLS::Flavour m_flav;
ATOOLS::Vec4D m_mom;
char m_info;
double m_mass, m_kt2max;
Proto_Particle * p_partner;
public:
Proto_Particle(const Proto_Particle &);
Proto_Particle(ATOOLS::Flavour flav,ATOOLS::Vec4D mom,char info);
......
......@@ -21,11 +21,11 @@ Constituents::Constituents(bool no_diquarks) :
double norm = 1./total;
Flavour flav(Flavour(kf_d));
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),1,udfrac*norm);
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),1,2.*udfrac*norm);
flav = Flavour(kf_u);
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),1,udfrac*norm);
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),1,2.*udfrac*norm);
flav = Flavour(kf_s);
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),1,sfrac*norm);
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),1,2.*sfrac*norm);
flav = Flavour(kf_c);
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),1,0.);
flav = Flavour(kf_b);
......@@ -35,30 +35,33 @@ Constituents::Constituents(bool no_diquarks) :
if (!no_diquarks && bfrac>0.) {
// Light Di-quarks, spin 0
flav = Flavour(kf_ud_0);
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),0,ud0*norm);
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),0,bfrac*ud0*norm);
flav = Flavour(kf_sd_0);
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),0,qssup*norm);
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),0,bfrac*qssup*norm);
flav = Flavour(kf_su_0);
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),0,qssup*norm);
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),0,bfrac*qssup*norm);
// Light Di-quarks, spin 1
flav = Flavour(kf_uu_1);
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),2,3.*sp1sup*norm);
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),2,3.*bfrac*sp1sup*norm);
flav = Flavour(kf_ud_1);
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),2,3.*sp1sup*norm);
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),2,3.*bfrac*sp1sup*norm);
flav = Flavour(kf_dd_1);
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),2,3.*sp1sup*norm);
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),2,3.*bfrac*sp1sup*norm);
flav = Flavour(kf_su_1);
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),2,3.*sp1sup*qssup*norm);
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),2,3.*bfrac*sp1sup*qssup*norm);
flav = Flavour(kf_sd_1);
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),2,3.*sp1sup*qssup*norm);
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),2,3.*bfrac*sp1sup*qssup*norm);
flav = Flavour(kf_ss_1);
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),2,3.*sp1sup*sssup*norm);
CCMap[flav] = new ConstituentCharacteristic(flav.HadMass(),2,3.*bfrac*sp1sup*sssup*norm);
}
//double totwt(0.);
for (FlavCCMap_Iterator cmit=CCMap.begin(); cmit!=CCMap.end();cmit++) {
if (cmit->first!=Flavour(kf_gluon) &&
cmit->second->Mass()<m_minmass) m_minmass = cmit->second->Mass();
//totwt += cmit->second->TotWeight();
//std::cout<<"Weight for "<<cmit->first<<" : "<<cmit->second->TotWeight()<<"."<<std::endl;
}
}
......
......@@ -8,13 +8,18 @@ using namespace ATOOLS;
Dipole_Splitter::Dipole_Splitter(Strong_Coupling * as,const bool & analyse) :
p_tools(new Splitting_Tools(as,analyse)),
m_zform(DefineZForm(int(hadpars.Get(std::string("Z_Form"))))),
Dipole_Splitter::Dipole_Splitter(Strong_Coupling * as,const leading::code & lead,
const bool & analyse) :
p_tools(new Splitting_Tools(as,lead,analyse)),
m_ptorder(DefinePTOrder(int(hadpars.Get(std::string("PT_Order"))))),
m_zform(DefineZForm(int(hadpars.Get(std::string("Z_Form"))))),
m_analyse(analyse)
{
if (m_analyse) {
m_histograms[std::string("PT_G")] = new Histogram(0,0.,25.,50);
m_histograms[std::string("PT_Q")] = new Histogram(0,0.,25.,50);
m_histograms[std::string("PT_primary")] = new Histogram(0,0.,25.,50);
m_histograms[std::string("ZQ")] = new Histogram(0,0.,2.,50);
m_histograms[std::string("ZQ")] = new Histogram(0,0.,2.,50);
m_histograms[std::string("XQ")] = new Histogram(0,0.,1.,50);
m_histograms[std::string("EQ")] = new Histogram(0,0.,10.,50);
......@@ -196,7 +201,7 @@ void Dipole_Splitter::AnalyseClusterSplitting(Cluster * cluster) {
void Dipole_Splitter::AnalyseKinematics(const ATOOLS::Vec4D q1,const ATOOLS::Vec4D q2,
const ATOOLS::Vec4D q3,const bool glusplit) {
double Ebeam = rpa.gen.Ecms()/2.;
//double Ebeam = rpa.gen.Ecms()/2.;
Histogram * histo;
double kt2(p_tools->PT2()), kt(sqrt(kt2)), z(p_tools->Z());
if (glusplit) {
......
......@@ -31,7 +31,8 @@ namespace AHADIC {
void AnalyseKinematics(const ATOOLS::Vec4D,const ATOOLS::Vec4D,
const ATOOLS::Vec4D,const bool=true);
public:
Dipole_Splitter(Strong_Coupling *,const bool & analyse);
Dipole_Splitter(Strong_Coupling *,const leading::code & lead,
const bool & analyse);
~Dipole_Splitter();
bool SplitCluster(Cluster * clu,const double pt2max=-1.);
......@@ -45,6 +46,8 @@ namespace AHADIC {
const double & PT2() const;
const double & Mass2() const;
const leading::code & Leading() const { return p_tools->Leading(); }
};
inline void Dipole_Splitter::SetOptions(FlavDecayMap * options) {
......
......@@ -278,14 +278,12 @@ void Hadron_Init::Init() {
s_kftable[3324]=new Particle_Info(3324,1.5318,0.0091,0,0,3,1,0,"Xi(1530)","Xi(1530)");
s_kftable[3334]=new Particle_Info(3334,1.67245,8.01e-15,-3,0,3,1,0,"Omega-","Omega-");
s_kftable[4114]=new Particle_Info(4114,2.5175,0.0150,0,0,3,1,0,"Sigma(c)(2520)","Sigma(c)(2520)");
s_kftable[4124]=new Particle_Info(4124,2.625,0.002,3,0,3,1,0,"Lambda(c)(2625)+","Lambda(c)(2625)+");
s_kftable[4214]=new Particle_Info(4214,2.5159,0.0150,3,0,3,1,0,"Sigma(c)(2520)+","Sigma(c)(2520)+");
s_kftable[4224]=new Particle_Info(4224,2.5194,0.0150,6,0,3,1,0,"Sigma(c)(2520)++","Sigma(c)(2520)++");
s_kftable[4314]=new Particle_Info(4314,2.63,0.003,0,0,3,1,0,"Xi(c)*","Xi(c)*");
s_kftable[4324]=new Particle_Info(4324,2.63,0.003,3,0,3,1,0,"Xi(c)*+","Xi(c)*+");
s_kftable[4334]=new Particle_Info(4334,2.8,0.001,0,0,3,1,0,"Omega(c)*","Omega(c)*");
s_kftable[5114]=new Particle_Info(5114,5.81,0.01,-3,0,3,1,0,"Sigma(b)*-","Sigma(b)*-");
s_kftable[5124]=new Particle_Info(5124,5.91,0.002,0,0,3,1,0,"Lambda(b)(5910)","Lambda(b)(5910)");
s_kftable[5214]=new Particle_Info(5214,5.81,0.01,0,0,3,1,0,"Sigma(b)*","Sigma(b)*");
s_kftable[5224]=new Particle_Info(5224,5.81,0.01,3,0,3,1,0,"Sigma(b)*+","Sigma(b)*+");
s_kftable[5314]=new Particle_Info(5314,5.98,0.001,-3,0,3,1,0,"Xi(b)*-","Xi(b)*-");
......@@ -322,6 +320,8 @@ void Hadron_Init::Init() {
s_kftable[1214]=new Particle_Info(1214,1.52,0.12,0,0,3,0,0,"N(1520)","N(1520)");
s_kftable[2124]=new Particle_Info(2124,1.52,0.12,3,0,3,0,0,"N(1520)+","N(1520)+");
s_kftable[3124]=new Particle_Info(3124,1.5195,0.0156,0,0,3,0,0,"Lambda(1520)","Lambda(1520)");
s_kftable[4124]=new Particle_Info(4124,2.625,0.002,3,0,3,1,0,"Lambda(c)(2625)+","Lambda(c)(2625)+");
s_kftable[5124]=new Particle_Info(5124,5.91,0.002,0,0,3,1,0,"Lambda(b)(5910)","Lambda(b)(5910)");
s_kftable[31214]=new Particle_Info(31214,1.72,0.2,0,0,3,0,0,"N(1720)","N(1720)");
s_kftable[32124]=new Particle_Info(32124,1.72,0.2,3,0,3,0,0,"N(1720)+","N(1720)+");
s_kftable[33122]=new Particle_Info(33122,1.67,0.06,0,0,1,0,0,"Lambda(1670)","Lambda(1670)");
......
......@@ -55,6 +55,8 @@ void All_Hadron_Multiplets::ConstructWaveFunctions()
fl2 = int(kfcode/100)-int(kfcode/1000)*10;
fl3 = int(kfcode/1000)-int(kfcode/10000)*10;
//std::cout<<"Construct wave function for "<<hadron<<" --> "
// <<fl1<<" "<<fl2<<" "<<fl3<<" spin = "<<kfcode-10*int(kfcode/10)<<"."<<std::endl;
if (fl3==0) {
if (int(fl2/2.)!=fl2/2.) { help = fl1; fl1 = fl2; fl2 = help; }
wavefunction = ConstructMesonWaveFunction(int(kfcode/9000000),
......@@ -190,8 +192,7 @@ All_Hadron_Multiplets::ConstructBaryonWaveFunction(int lp,int spin,
// Baryon wave functions according to Lichtenberg, Namgung, Wills & Predazzi
if (spin==0) return NULL;
int wf=-1;
int pos1,pos2,pos3,di;
int wf(-1),pos1(-1),pos2(-1),pos3(-1),di(-1);
if ((spin==2 || lp!=0) && (fl1<fl2 || fl2<fl3)) {
// Octet
......@@ -245,6 +246,15 @@ All_Hadron_Multiplets::ConstructBaryonWaveFunction(int lp,int spin,
}
}
if (pos1<0 || pos2<0 || pos3<0 || wf<0) {
msg_Tracking()<<"Warning in "<<METHOD<<"("<<lp<<","<<spin<<","<<fl1<<","<<fl2<<","<<fl3<<"):"
<<std::endl
<<" leads to unknown combination of wf type and positions."<<std::endl
<<" wf = "<<wf<<", pos123 = "<<pos1<<" "<<pos2<<" "<<pos3<<";"
<<" will ignore this hadron."<<std::endl;
return NULL;
}
Hadron_Wave_Function * wavefunction = new Hadron_Wave_Function;
Flavour_Pair * pair;
......@@ -406,7 +416,7 @@ void All_Hadron_Multiplets::LookUpAngles(const int angular,const int spin,double
void All_Hadron_Multiplets::CreateMultiplets()
{
p_multiplets = new Hadron_Multiplet_Map;
int kfcode,dkfc,spin,orbital,radial,iso0;
int kfcode,spin,orbital,radial,iso0;
Hadron_Multiplet_Miter mpletiter;
Hadron_Multiplet * multiplet;
std::string prefix;
......@@ -627,7 +637,6 @@ Hadron_Wave_Function * All_Hadron_Multiplets::GetWaveFunction(ATOOLS::Flavour fl
void All_Hadron_Multiplets::PrintWaveFunctions()
{
Hadron_WF_Miter wfm;
Hadron_Wave_Function * wf;
map<Flavour,double> checkit;
for (Hadron_Multiplet_Miter mplet=p_multiplets->begin();mplet!=p_multiplets->end();mplet++) {
if (mplet->second->Weight()<=0.) continue;
......@@ -636,10 +645,9 @@ void All_Hadron_Multiplets::PrintWaveFunctions()
for (FlSetIter fl=mplet->second->GetElements()->begin();
fl!=mplet->second->GetElements()->end();fl++) {
wfm = p_wavefunctions->find((*fl));
if (wfm!=p_wavefunctions->end()) { wf = wfm->second; msg_Out()<<(*wf); }
if (wfm!=p_wavefunctions->end()) msg_Out()<<(*wfm->second);
if (mplet->second->Name()==string("Nucleons (1/2)")) {
WFcomponent * wfc = wf->GetWaves();
WFcomponent * wfc = wfm->second->GetWaves();
for (WFcompiter cit = wfc->begin();cit!=wfc->end();cit++) {
if (checkit.find(cit->first->first)==checkit.end())
checkit[cit->first->first] = sqr(cit->second);
......
......@@ -18,7 +18,8 @@ Hadronisation_Parameters::Hadronisation_Parameters() :
p_constituents(NULL),p_multiplets(NULL),
p_singletransitions(NULL),p_doubletransitions(NULL),
p_softclusters(NULL),
m_asform(asform::constant), p_coupling(NULL), p_splitter(NULL)
m_leading(leading::none),m_asform(asform::constant),
p_coupling(NULL), p_splitter(NULL)
{ }
Hadronisation_Parameters::~Hadronisation_Parameters() {
......@@ -55,12 +56,24 @@ void Hadronisation_Parameters::Init(string dir,string file)
p_doubletransitions = new Double_Transitions();
if (msg_LevelIsTracking())
p_doubletransitions->PrintDoubleTransitions();
//abort();
//exit(1);
switch (int(m_parametermap[string("asform")])) {
case 64:
m_asform = asform::Exponential;
switch (int(m_parametermap[string("leading")])) {
case 3:
m_leading = leading::quarks_and_gluons2;
break;
case 2:
m_leading = leading::quarks_and_gluons;
break;
case 1:
m_leading = leading::only_quarks;
break;
case 0:
default:
m_leading = leading::none;
break;
}
switch (int(m_parametermap[string("asform")])) {
case 8:
m_asform = asform::GDH_inspired;
break;
......@@ -75,8 +88,8 @@ void Hadronisation_Parameters::Init(string dir,string file)
m_asform = asform::constant;
break;
}
p_coupling = new Strong_Coupling(m_asform,dabs(m_parametermap[string("pt02")]));
p_splitter = new Dipole_Splitter(p_coupling,ana);
p_coupling = new Strong_Coupling(m_asform);
p_splitter = new Dipole_Splitter(p_coupling,m_leading,ana);
p_softclusters = new Soft_Cluster_Handler(ana);
}
......@@ -88,6 +101,8 @@ void Hadronisation_Parameters::ReadParameters(string dir,string file)
dataread.SetInputPath(dir);
dataread.SetInputFile(file);
// General switches for operational modes
m_parametermap[string("leading")] =
dataread.GetValue<int>("LEADING",1);
m_parametermap[string("asform")] =
dataread.GetValue<int>("ALPHAS_FORM",1);
m_parametermap[string("Mass_treatment")] =
......@@ -98,12 +113,13 @@ void Hadronisation_Parameters::ReadParameters(string dir,string file)
dataread.GetValue<int>("DECAY_SELECTION",3);
m_parametermap[string("PT_Ordering")] =
dataread.GetValue<int>("PT_ORDERING",3);
m_parametermap[string("Z_Form")] =
m_parametermap[string("Z_Form")] =
dataread.GetValue<int>("Z_FORM",2);
// Parameters
m_parametermap[string("asfix")] =
dataread.GetValue<double>("AS_FIX",1.0);
m_parametermap[string("pt2min")] =
dataread.GetValue<double>("PT^2_MIN",0.09);
m_parametermap[string("pt02")] =
dataread.GetValue<double>("PT^2_0",9.7738e-01);
m_parametermap[string("ptmax_factor")] =
......@@ -116,15 +132,15 @@ void Hadronisation_Parameters::ReadParameters(string dir,string file)
dataread.GetValue<double>("TRANSITION_OFFSET",-1.4258e-03);
m_parametermap[string("Offset_C->HH")] =
dataread.GetValue<double>("DECAY_OFFSET",-1.9844e+00);
m_parametermap[string("MassExponent_C->H")] =
m_parametermap[string("MassExponent_C->H")] =
dataread.GetValue<double>("TRANSITION_EXPONENT",5.0236e-01);
m_parametermap[string("WidthExponent_C->H")] =
dataread.GetValue<double>("TRANSITION_EXPONENT2",1.3163e-02);
m_parametermap[string("MassExponent_C->HH")] =
dataread.GetValue<double>("DECAY_EXPONENT",3.7959e+00);
m_parametermap[string("Strange_fraction")] =
m_parametermap[string("Strange_fraction")] =
dataread.GetValue<double>("STRANGE_FRACTION",1.0);
m_parametermap[string("Baryon_fraction")] =
m_parametermap[string("Baryon_fraction")] =
dataread.GetValue<double>("BARYON_FRACTION",1.0);
m_parametermap[string("P_qs_by_P_qq")] =
dataread.GetValue<double>("P_{QS}/P_{QQ}",1.0);
......@@ -183,7 +199,7 @@ void Hadronisation_Parameters::ReadParameters(string dir,string file)
m_parametermap[string("Mass_glue")] =
dataread.GetValue<double>("M_GLUE",0.);
Flavour(kf_gluon).SetHadMass(m_parametermap["Mass_glue"]);
double mud = m_parametermap[string("Mass_down")] =
double mud = m_parametermap[string("Mass_updown")] =
dataread.GetValue<double>("M_UP_DOWN",0.01);
double ms = m_parametermap[string("Mass_strange")] =
dataread.GetValue<double>("M_STRANGE",0.2);
......
......@@ -28,6 +28,7 @@ namespace AHADIC {
Single_Transitions * p_singletransitions;
Double_Transitions * p_doubletransitions;
Soft_Cluster_Handler * p_softclusters;
leading::code m_leading;
asform::code m_asform;
Strong_Coupling * p_coupling;
Dipole_Splitter * p_splitter;
......
This diff is collapsed.
......@@ -37,7 +37,8 @@ namespace AHADIC {
TransWeight::code m_tweightmode;
DecayWeight::code m_dweightmode;
double m_transitionoffset,m_decayoffset,m_kappa,m_lambda,m_chi,m_pt2max;
double m_transitionoffset,m_decayoffset,m_kappa,m_lambda,m_chi;
double m_pt2min,m_pt2maxfac;
Cluster_List * m_clin;
long int m_transitions, m_dtransitions, m_decays, m_forceddecays, m_lists;
......
This diff is collapsed.
......@@ -15,6 +15,16 @@ namespace AHADIC {
typedef std::map<ATOOLS::Flavour,DecaySpecs *> FlavDecayMap;
typedef FlavDecayMap::iterator FDIter;
struct leading {
enum code {
none = 0,
only_quarks = 1,
quarks_and_gluons = 2,
quarks_and_gluons2 = 3
};
};
struct PTOrder {
enum code {
total = 3,
......@@ -25,6 +35,7 @@ namespace AHADIC {
};
PTOrder::code DefinePTOrder(const int & ptorder);
struct ZForm {
enum code {
fragmentation = 30,
......@@ -37,10 +48,12 @@ namespace AHADIC {
ZForm::code DefineZForm(const int & zform);
class Splitting_Functions {
private:
int m_masstreatment;
double m_alpha,m_beta,m_sigma,m_kappa;
double m_alpha,m_sigma,m_kappa;
double m_m1_2,m_m2_2,m_m3_2,m_Q2,m_y,m_viji,m_vijk,m_zp,m_zm,m_pipj;
......@@ -78,6 +91,7 @@ namespace AHADIC {