Commit d27bd8e9 authored by Frank Krauss's avatar Frank Krauss

more minor changes to accommodate for a tuning

parent 42e77e1a
......@@ -34,26 +34,24 @@ Cluster_Part::~Cluster_Part()
bool Cluster_Part::TestDecay(Cluster * const cluster)
{
if (cluster->GetTrip()->m_info=='B' || cluster->GetAnti()->m_info=='B') {
msg_Tracking()<<":::::::::::::::::::::::::::::::::::::::::::::::"<<std::endl
<<"::: "<<METHOD<<" : Try "<<cluster->Number()<<" ("
<<cluster->GetTrip()->m_flav<<" "<<cluster->GetAnti()->m_flav<<", "
<<"m = "<<cluster->Mass()<<") --> "<<std::endl;
}
msg_Tracking()<<":::::::::::::::::::::::::::::::::::::::::::::::"<<std::endl
<<"::: "<<METHOD<<" : Try "<<cluster->Number()<<" ("
<<cluster->GetTrip()->m_flav<<" "<<cluster->GetAnti()->m_flav<<", "
<<"m = "<<cluster->Mass()<<") --> "<<std::endl;
if (!p_splitter->SplitCluster(cluster)) {
if (cluster->GetTrip()->m_info=='B' || cluster->GetAnti()->m_info=='B') {
msg_Tracking()<<"::: Warning in "<<METHOD<<":"<<std::endl
<<"::: Could not split cluster ("<<cluster->Number()<<"): "
<<cluster->GetTrip()->m_flav<<"/"<<cluster->GetAnti()->m_flav<<", "
<<"mass = "<<cluster->Mass()<<","<<std::endl
<<"::: try to enforce splitting (not implemented)."<<std::endl;
}
if (cluster->GetTrip()->m_info=='B' || cluster->GetAnti()->m_info=='B') {
msg_Tracking()<<"::: "<<METHOD<<" : yields enforced decay of "
<<cluster->Number()<<" succeded."<<std::endl
msg_Tracking()<<"::: Warning in "<<METHOD<<":"<<std::endl
<<"::: Could not split cluster ("<<cluster->Number()<<"): "
<<cluster->GetTrip()->m_flav<<"/"<<cluster->GetAnti()->m_flav<<", "
<<"mass = "<<cluster->Mass()<<","<<std::endl
<<"::: try to enforce splitting (not implemented)."<<std::endl;
bool worked(p_splitter->EnforceSplit(cluster));
if (worked) {
msg_Tracking()<<"::: "<<METHOD<<" : forced decay of "<<cluster->Momentum()<<" succeded."<<std::endl
<<(*cluster->GetLeft())<<(*cluster->GetRight())
<<":::::::::::::::::::::::::::::::::::::::::::::::"<<std::endl;
}
return p_splitter->EnforceSplit(cluster);
return worked;
}
if (m_ana) {
Vec4D lmom(cluster->GetLeft()->Momentum());
......@@ -64,9 +62,8 @@ bool Cluster_Part::TestDecay(Cluster * const cluster)
#ifdef AHAmomcheck
cluster->CheckConsistency(msg_Error(),METHOD);
#endif
if (cluster->GetTrip()->m_info=='B' || cluster->GetAnti()->m_info=='B') {
msg_Tracking()<<"::: "<<METHOD<<" : decay of "<<cluster->Number()<<" succeded."<<std::endl
<<":::::::::::::::::::::::::::::::::::::::::::::::"<<std::endl;
}
msg_Tracking()<<"::: "<<METHOD<<" : decay of "<<cluster->Momentum()<<" succeded."<<std::endl
<<(*cluster->GetLeft())<<(*cluster->GetRight())
<<":::::::::::::::::::::::::::::::::::::::::::::::"<<std::endl;
return true;
}
......@@ -55,6 +55,7 @@ bool Dipole_Splitter::SplitCluster(Cluster * cluster) {
cluster->GetAnti()->m_flav.IsDiQuark());
if (!EmitGluon(dip1,dip2)) return false;
double kt2_emit(p_tools->PT2()), z_emit(p_tools->Z());
bool swapped(SelectOrder(dip1,dip2));
......@@ -71,6 +72,7 @@ bool Dipole_Splitter::SplitCluster(Cluster * cluster) {
return false;
}
}
double kt2_split(p_tools->PT2()), z_split(p_tools->Z());
Proto_Particle * out1, * out2;
p_tools->GetNewParticles(out1,out2);
......@@ -86,22 +88,29 @@ bool Dipole_Splitter::SplitCluster(Cluster * cluster) {
dip1->Update();
dip2->Update();
if (cluster->Mass()>80.) {
msg_Tracking()<<"::: split massive cluster ("<<cluster->GetTrip()->m_info<<cluster->GetAnti()->m_info<<", "
<<"mass = "<<cluster->Mass()<<") with "
<<"kt_emit = "<<sqrt(kt2_emit)<<" & z_emit = "<<z_emit<<", "
<<"kt_split = "<<sqrt(kt2_split)<<" & z_split = "<<z_split<<") "
<<" --> "<<sqrt(dip1->MassBar2())<<" + "<<sqrt(dip2->MassBar2())<<std::endl;
}
else {
msg_Tracking()<<"::: split cluster ("<<cluster->GetTrip()->m_info<<cluster->GetAnti()->m_info<<", "
<<"mass = "<<cluster->Mass()<<") with "
<<"kt_emit = "<<sqrt(kt2_emit)<<" & z_emit = "<<z_emit<<", "
<<"kt_split = "<<sqrt(kt2_split)<<" & z_split = "<<z_split<<") "
<<" --> "<<sqrt(dip1->MassBar2())<<" + "<<sqrt(dip2->MassBar2())<<std::endl;
}
return Produce2Clusters(cluster,dip1,dip2);
}
bool Dipole_Splitter::EmitGluon(Dipole * dip1,Dipole *& dip2) {
double pt2max(p_tools->SetSpectatorAndSplitter(dip1));
if (dip1->Triplet()->m_info=='B' || dip1->AntiTriplet()->m_info=='B') {
msg_Debugging()<<"==============================================================="<<std::endl
<<"=== "<<METHOD<<"(pt_max = "<<sqrt(pt2max)<<", "
<<"m^2 = "<<dip1->Mass2()<<") for "
<<dip1->Triplet()->m_flav<<"("<<dip1->Triplet()->m_info<<") "
<<dip1->AntiTriplet()->m_flav<<"("<<dip1->AntiTriplet()->m_info<<")"<<std::endl
<<" "<<dip1->Triplet()->m_mom<<" (kt^2_max = "<<dip1->Triplet()->m_kt2max<<") "
<<dip1->AntiTriplet()->m_mom<<" (kt^2_max = "<<dip1->AntiTriplet()->m_kt2max<<")"
<<"."<<std::endl;
}
if (!p_tools->PrepareKinematics(dip1) || !p_tools->DetermineSplitting(dip1)) {
bool first(false);
if (dip1->Triplet()->m_info=='L' && dip1->AntiTriplet()->m_info=='L' &&
dip1->Triplet()->m_kt2max==dip1->AntiTriplet()->m_kt2max) first = true;
if (!p_tools->PrepareKinematics(dip1,first) || !p_tools->DetermineSplitting(dip1)) {
pt2max = p_tools->SwapSpectatorAndSplitter(dip1);
if (dip1->Triplet()->m_info=='B' || dip1->AntiTriplet()->m_info=='B') {
msg_Debugging()<<"==============================================================="<<std::endl
......@@ -113,7 +122,7 @@ bool Dipole_Splitter::EmitGluon(Dipole * dip1,Dipole *& dip2) {
<<dip1->AntiTriplet()->m_mom<<" (kt^2_max = "<<dip1->AntiTriplet()->m_kt2max<<")"
<<"."<<std::endl;
}
if (!p_tools->PrepareKinematics(dip1) || !p_tools->DetermineSplitting(dip1)) {
if (!p_tools->PrepareKinematics(dip1,first) || !p_tools->DetermineSplitting(dip1)) {
delete dip1->Triplet();
delete dip1->AntiTriplet();
delete dip1;
......@@ -158,7 +167,7 @@ bool Dipole_Splitter::Produce2Clusters(Cluster * cluster,Dipole *& dip1,Dipole *
#endif
if (m_analyse) AnalyseClusterSplitting(cluster);
msg_Tracking()<<(*cluster)<<std::endl;
msg_Debugging()<<(*cluster)<<std::endl;
return true;
}
......@@ -175,12 +184,16 @@ bool Dipole_Splitter::SplitDipole(Dipole * dip,const bool & first,const bool & v
if (!p_tools->PrepareKinematics(dip,first) ||
!p_tools->DetermineSplitting(dip,vetodiquark)) {
msg_Tracking()<<"=== "<<METHOD<<"(out): could not split dipole."<<std::endl;
msg_Debugging()<<"=== "<<METHOD<<"(out): could not split dipole."<<std::endl;
return false;
}
p_tools->AftermathOfSplitting(dip);
if (first)
msg_Tracking()<<"::: split initial gluon with "
<<"kt = "<<sqrt(p_tools->PT2())<<", z = "<<p_tools->Z()<<", "
<<"dipole mass = "<<sqrt(dip->Mass2())<<"."<<std::endl;
msg_Tracking()<<"=== "<<METHOD<<"(out): succeded to split dipole."<<std::endl;
msg_Debugging()<<"=== "<<METHOD<<"(out): succeded to split dipole."<<std::endl;
return true;
}
......@@ -194,12 +207,12 @@ bool Dipole_Splitter::EnforceSplit(Dipole * dip1,Dipole * dip2) {
<<" "<<dipole->Triplet()->m_mom<<" ( "<<dipole->Triplet()->m_kt2max<<" ) "
<<dipole->AntiTriplet()->m_mom<<" ( "<<dipole->AntiTriplet()->m_kt2max<<" )."<<std::endl;
if (!p_tools->PrepareKinematics(dipole,false,true) || !p_tools->EnforcedSplitting(dipole)) {
msg_Tracking()<<"=== "<<METHOD<<"(out): could not split dipole."<<std::endl;
msg_Debugging()<<"=== "<<METHOD<<"(out): could not split dipole."<<std::endl;
return false;
}
p_tools->AftermathOfSplitting(dipole);
msg_Tracking()<<"=== "<<METHOD<<"(out): succeded to split dipole."<<std::endl;
msg_Debugging()<<"=== "<<METHOD<<"(out): succeded to split dipole."<<std::endl;
return true;
}
......
......@@ -127,7 +127,7 @@ All_Hadron_Multiplets::ConstructMesonWaveFunction(const int iso0,const int rp,co
wavefunction->AddToWaves(pair,weight);
}
else {
weight = (sinth/sqrt(6.)+costh/sqrt(3.))*costh*m_singletsuppression;
weight = (sinth/sqrt(6.)+costh/sqrt(3.))*(1.-m_singletsuppression);
if (dabs(weight)>1.e-3) {
wavefunction = new Hadron_Wave_Function;
wavefunction->AddToWaves(pair,weight);
......@@ -137,7 +137,7 @@ All_Hadron_Multiplets::ConstructMesonWaveFunction(const int iso0,const int rp,co
pair->second = flavs[0].Bar();
wavefunction->AddToWaves(pair,weight);
}
weight = (-2.*sinth/sqrt(6.)+costh/sqrt(3.))*costh*m_singletsuppression;
weight = (-2.*sinth/sqrt(6.)+costh/sqrt(3.))*(1.-m_singletsuppression);
if (dabs(weight)>1.e-3) {
flavs[0] = Flavour((kf_code)(3));
pair = new Flavour_Pair;
......@@ -149,12 +149,12 @@ All_Hadron_Multiplets::ConstructMesonWaveFunction(const int iso0,const int rp,co
}
}
else if (fl1==3) {
weight = (-2.*costh/sqrt(6.)-sinth/sqrt(3.))*sinth*m_singletsuppression;
weight = (-2.*costh/sqrt(6.)-sinth/sqrt(3.))*m_singletsuppression;
if (dabs(weight)>1.e-3) {
wavefunction = new Hadron_Wave_Function;
wavefunction->AddToWaves(pair,weight);
}
weight = (costh/sqrt(6.)-sinth/sqrt(3.))*sinth*m_singletsuppression;
weight = (costh/sqrt(6.)-sinth/sqrt(3.))*m_singletsuppression;
if (dabs(weight)>1.e-3) {
flavs[0] = Flavour((kf_code)(1));
pair = new Flavour_Pair;
......
This diff is collapsed.
This diff is collapsed.
......@@ -71,8 +71,8 @@ namespace AHADIC {
void SelectFlavour(const double & maxmass,const bool & vetodiquark);
bool KinCheck(const double & kt2,const double & z);
bool ConstructKinematics(const bool glusplit);
double KT2Max(const bool & first=false);
double KT2Min(const double & kt2max,const bool & glusplit);
double KT2Max(const double & kt2min,const bool & first=false);
double KT2Min(const bool & glusplit,const bool & first=false);
void SetInfoTagsForOutgoings() const;
void AnalyseKinematics(const ATOOLS::Vec4D & q1,const ATOOLS::Vec4D & q2,
const ATOOLS::Vec4D & q3);
......
......@@ -88,7 +88,8 @@ const double Strong_Coupling::operator()(double q2,bool reweight) const {
const double Strong_Coupling::SelectPT(const double & scale2max,const double & scale2min) {
double pt2(0.), pt2max(Min(scale2max,m_pt2max)), ran1;
double mini(m_pt02+scale2min),maxi(m_pt02+pt2max),expo(m_eta==1.?0.:1./1.-m_eta);;
double mini(m_pt02+scale2min),maxi(m_pt02+pt2max),expo(dabs(m_eta-1.)<0.01?0.:1.-m_eta), invexpo(1./expo);
double mini_pow(pow(mini,expo)), maxi_pow(pow(maxi,expo));
bool runit(true);
while (runit) {
ran1 = ran.Get();
......@@ -123,8 +124,10 @@ const double Strong_Coupling::SelectPT(const double & scale2max,const double & s
case asform::GDH_inspired:
case asform::constant:
default:
if (m_eta==1.) pt2 = -m_pt02+mini*pow(maxi/mini,ran1);
else pt2 = -m_pt02+pow(mini*(1.-ran1)+maxi*ran1,expo);
if (m_eta==1.)
pt2 = -m_pt02+mini*pow(maxi/mini,ran1);
else
pt2 = -m_pt02+pow(ran1*maxi_pow+(1.-ran1)*mini_pow,invexpo);
if ((*this)(pt2)/m_asmax>ran.Get()) runit = false;
break;
}
......
......@@ -189,14 +189,13 @@ Double_Transitions::Double_Transitions() :
anti = iter2->first;
if (anti.IsQuark()) anti = anti.Bar();
flpair.second = wf2.second = anti;
//std::cout<<" "<<trip<<" / "<<anti<<std::endl;
//if (trip==Flavour(kf_u) && anti==Flavour(kf_u).Bar())
// std::cout<<" "<<trip<<" / "<<anti<<std::endl;
for (FlavCCMap_Iterator iter3=constituents.begin();
iter3!=constituents.end();iter3++) {
norm = 0.;
popwt = iter3->second->TotWeight();
//std::cout<<" pop "<<iter3->first<<" ( wt = "<<popwt<<") "
// <<" between "<<trip<<" and "<<anti<<"."<<std::endl;
if (popwt==0.) continue;
popped = iter3->first;
if (popped.IsQuark()) popped = popped.Bar();
......@@ -217,10 +216,12 @@ Double_Transitions::Double_Transitions() :
haditer1!=hads1->end();haditer1++) {
for (Single_Transition_Siter haditer2=hads2->begin();
haditer2!=hads2->end();haditer2++) {
weight = haditer1->second * haditer2->second/norm * popwt;
//std::cout<<" add : "<<weight/popwt<<" for "
// <<haditer1->first<<"("<<haditer1->second<<") * "
// <<haditer2->first<<"("<<haditer2->second<<")."<<std::endl;
weight = haditer1->second * haditer2->second * popwt;
//if (trip==Flavour(kf_u) && anti==Flavour(kf_u).Bar())
// std::cout<<" add : "<<weight/popwt<<" "
// <<"("<<(haditer1->second*haditer2->second)<<") for "
// <<haditer1->first<<"("<<haditer1->second<<") * "
// <<haditer2->first<<"("<<haditer2->second<<")."<<std::endl;
hadpair.first = haditer1->first;
hadpair.second = haditer2->first;
if (weight>0.) {
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment