Commit 5963ec78 authored by Frank Krauss's avatar Frank Krauss

Ironing out last issues in AHADIC

parent 3a02a05f
......@@ -9,6 +9,7 @@ using namespace ATOOLS;
Cluster_Splitter::Cluster_Splitter(std::list<Cluster *> * cluster_list,
Soft_Cluster_Handler * softclusters) :
Splitter_Base(cluster_list,softclusters),
m_mode(0),
m_output(false)
{}
......@@ -27,15 +28,88 @@ void Cluster_Splitter::Init(const bool & isgluon) {
bool Cluster_Splitter::MakeLongitudinalMomenta() {
CalculateLimits();
FixCoefficients(p_part1->Flavour(),p_part2->Flavour());
m_z1 = m_zselector(m_z1min,m_z1max);
FixCoefficients(p_part2->Flavour(),p_part1->Flavour());
m_z2 = m_zselector(m_z2min,m_z2max);
m_R12 = m_z1*(1.-m_z2)*m_Q2-m_kt2;
m_R21 = (1.-m_z1)*m_z2*m_Q2-m_kt2;
if (m_mode==0) {
FixCoefficients(p_part1->Flavour(),p_part2->Flavour());
m_z1 = m_zselector(m_z1min,m_z1max);
FixCoefficients(p_part2->Flavour(),p_part1->Flavour());
m_z2 = m_zselector(m_z2min,m_z2max);
m_R12 = m_z1*(1.-m_z2)*m_Q2-m_kt2;
m_R21 = (1.-m_z1)*m_z2*m_Q2-m_kt2;
}
else if (m_mode==1) {
do {
FixCoefficients(p_part1->Flavour(),p_part2->Flavour());
m_R12 = m_m12min + DeltaM2();
FixCoefficients(p_part2->Flavour(),p_part1->Flavour());
m_R21 = m_m22min + DeltaM2();
} while (m_R12+m_R21>m_Q2);
double e12 = (m_R12+m_kt2)/m_Q2, e21 = (m_R21+m_kt2)/m_Q2;
double disc = sqr(1-e12-e21)-4.*e12*e21;
if (disc<0.) return false;
disc = sqrt(disc);
m_z1 = (1.+e12-e21+disc)/2.;
m_z2 = (1.-e12+e21+disc)/2.;
}
return CheckIfAllowed();
}
void Cluster_Splitter::
FixCoefficients(const Flavour & flav1,const Flavour & flav2) {
m_flcnt = 0;
if (false && !(p_part1->IsLeading() || p_part2->IsLeading())) {
m_a = m_b = m_c = 0.;
return;
}
if (flav1==Flavour(kf_b) || flav1==Flavour(kf_b).Bar()) { //||
// flav1==Flavour(kf_c) || flav1==Flavour(kf_c).Bar()) {
m_flcnt = 1;
}
else if (flav1.IsDiQuark()) {
m_flcnt = 2;
}
m_a = m_alpha[m_flcnt];
m_b = m_beta[m_flcnt];
m_c = m_gamma[m_flcnt];
m_masses2 = sqr(p_constituents->Mass(flav1) +
p_constituents->Mass(flav2));
}
void Cluster_Splitter::CalculateLimits() {
m_m12min =
sqr(Min(p_softclusters->MinSingleMass(m_flavs1.first,m_newflav1),
p_softclusters->MinDoubleMass(m_flavs1.first,m_newflav1)));
m_m22min =
sqr(Min(p_softclusters->MinSingleMass(m_newflav2,m_flavs2.second),
p_softclusters->MinDoubleMass(m_newflav2,m_flavs2.second)));
double lambda = sqrt(sqr(m_Q2-m_m12min-m_m22min)-
4.*(m_m12min+m_kt2)*(m_m22min+m_kt2));
double centre1 = m_Q2-m_m22min+m_m12min;
double centre2 = m_Q2-m_m12min+m_m22min;
m_z1min = (centre1)/(2.*m_Q2);
m_z1max = (centre1+lambda)/(2.*m_Q2);
m_z2min = (centre2)/(2.*m_Q2);
m_z2max = (centre2+lambda)/(2.*m_Q2);
}
double Cluster_Splitter::DeltaM2() {
double deltaMmax = sqrt(m_Q2-m_m12min-m_m22min);
//double wtmax = pow(m_a*m_c,m_a)*exp(-m_a);
double wtMmax = m_a*deltaMmax/(m_c*(m_a+m_b));
double wtmax = pow(wtMmax,m_a)*pow(1.-wtMmax,m_b);
double deltaM, weight;
do {
deltaM = ran->Get()*deltaMmax;
weight = pow(deltaM/m_c,m_a)*pow(1.-deltaM/m_c,m_b);
//weight = pow(deltaM,m_a)*exp(-deltaM/m_c);
} while (weight<ran->Get()*wtmax);
return sqr(deltaM);
}
bool Cluster_Splitter::CheckIfAllowed() {
bool allowed = (m_R12>sqr(m_minQ_1) && m_R21>sqr(m_minQ_2));
return allowed;
}
double Cluster_Splitter::
WeightFunction(const double & z,const double & zmin,const double & zmax) {
double arg = dabs(m_c)>1.e-2? m_c*(4.*m_kt2+m_masses2) : 0. ;
......@@ -44,7 +118,7 @@ WeightFunction(const double & z,const double & zmin,const double & zmax) {
norm *= Max(pow(zmin,m_a), pow(1.-zmax,m_b));
else {
if (m_a<=0.) norm *= pow(zmin,m_a);
if (m_b<=0.) norm *= pow(1.-zmax,m_b);
if (m_b<=0.) norm *= pow(1.-zmax,m_b);
}
double wt = pow(z,m_a) * pow(1.-z,m_b) * exp(-arg/z);
if (wt>norm) {
......@@ -56,11 +130,6 @@ WeightFunction(const double & z,const double & zmin,const double & zmax) {
return wt / norm;
}
bool Cluster_Splitter::CheckIfAllowed() {
bool allowed = ((m_R12>sqr(m_minQ_1) && m_R21>sqr(m_minQ_2)));
return allowed;
}
bool Cluster_Splitter::FillParticlesInLists() {
for (size_t i=0;i<2;i++) p_cluster_list->push_back(MakeCluster(i));
return true;
......@@ -114,40 +183,3 @@ Cluster * Cluster_Splitter::MakeCluster(size_t i) {
return cluster;
}
void Cluster_Splitter::CalculateLimits() {
double m12 =
sqr(Min(p_softclusters->MinSingleMass(m_flavs1.first,m_newflav1),
p_softclusters->MinDoubleMass(m_flavs1.first,m_newflav1)));
double m22 =
sqr(Min(p_softclusters->MinSingleMass(m_newflav2,m_flavs2.second),
p_softclusters->MinDoubleMass(m_newflav2,m_flavs2.second)));
double lambda = sqrt(sqr(m_Q2-m12-m22)-4.*(m12+m_kt2)*(m22+m_kt2));
double centre1 = m_Q2-m22+m12;
double centre2 = m_Q2-m12+m22;
m_z1min = (centre1)/(2.*m_Q2);
m_z1max = (centre1+lambda)/(2.*m_Q2);
m_z2min = (centre2)/(2.*m_Q2);
m_z2max = (centre2+lambda)/(2.*m_Q2);
}
void Cluster_Splitter::
FixCoefficients(const Flavour & flav1,const Flavour & flav2) {
m_flcnt = 0;
if (false && !(p_part1->IsLeading() || p_part2->IsLeading())) {
m_a = m_b = m_c = 0.;
return;
}
if (flav1==Flavour(kf_b) || flav1==Flavour(kf_b).Bar()) { //||
// flav1==Flavour(kf_c) || flav1==Flavour(kf_c).Bar()) {
m_flcnt = 1;
}
else if (flav1.IsDiQuark()) {
m_flcnt = 2;
}
m_a = m_alpha[m_flcnt];
m_b = m_beta[m_flcnt];
m_c = m_gamma[m_flcnt];
m_masses2 = sqr(p_constituents->Mass(flav1) +
p_constituents->Mass(flav2));
}
......@@ -6,6 +6,7 @@
namespace AHADIC {
class Cluster_Splitter : public Splitter_Base {
private:
int m_mode;
double m_alpha[3], m_beta[3], m_gamma[3], m_a, m_b, m_c;
size_t m_flcnt;
......@@ -15,8 +16,7 @@ namespace AHADIC {
bool m_output;
bool MakeLongitudinalMomenta();
double DeltaM2(const double & MC2,const double & mt2,
const ATOOLS::Flavour & flav);
double DeltaM2();
void FixCoefficients(const ATOOLS::Flavour & flav1,
const ATOOLS::Flavour & flav2);
void CalculateLimits();
......
......@@ -70,6 +70,7 @@ void Singlet_Checker::Reset() {
}
bool Singlet_Checker::operator()() {
Reset();
list<Singlet *>::iterator lsit(p_singlets->begin());
while (lsit!=p_singlets->end()) {
p_singlet = (*lsit);
......@@ -101,7 +102,7 @@ bool Singlet_Checker::operator()() {
if (msg_LevelIsTracking()) {
for (list<list<Singlet *>::iterator>::iterator bit=m_badones.begin();
bit!=m_badones.end();bit++) {
msg_Out()<<(***bit)<<"\n";
msg_Error()<<(***bit)<<"\n";
}
}
return false;
......@@ -115,9 +116,9 @@ bool Singlet_Checker::CheckSinglet() {
for (list<Proto_Particle *>::iterator plit=p_singlet->begin();
plit!=p_singlet->end();plit++) {
if ((*plit)->Momentum()[0]<0. || (*plit)->Momentum().Abs2()<-1.e-4) {
msg_Out()<<"Error in "<<METHOD<<":\n"
<<" negative energy or mass^2 particle in singlet:\n"
<<(*p_singlet)<<"n";
msg_Error()<<"Error in "<<METHOD<<":\n"
<<" negative energy or mass^2 particle in singlet:\n"
<<(*p_singlet)<<"n";
//return false;
}
}
......@@ -130,7 +131,7 @@ bool Singlet_Checker::CheckSinglet() {
(*plit2)->Flavour()),
p_softclusters->MinDoubleMass((*plit1)->Flavour(),
(*plit2)->Flavour()));
if (mass < minmass) return false;
return (mass > minmass);
}
while (plit2!=p_singlet->end()) {
// this is more or less a plausibility check driven by some external mass
......@@ -154,8 +155,8 @@ bool Singlet_Checker::CheckSinglet() {
bool Singlet_Checker::FusePartonsInLowMassSinglet() {
if (p_singlet->front()->Flavour().IsGluon() &&
sqrt(m_mass) > 2.*m_minQmass && m_splitter(p_part1,p_part2)) {
// gluons heavy enough to be replaced by two quarks, splits gluon ring
// and necessitates reordering the singlet to have a quark as first
// gluon system is heavy enough to be replaced by two quarks, splits gluon
// ring and necessitates reordering the singlet to have a quark as first
// particle.
p_singlet->Reorder();
return true;
......@@ -214,7 +215,7 @@ void Singlet_Checker::SortProblematicSinglets() {
p_singlet = (**bit);
Flavour flav1 = p_singlet->front()->Flavour();
Flavour flav2 = p_singlet->back()->Flavour();
if (!flav1.IsGluon()) {
if (!flav1.IsGluon() && !flav2.IsGluon()) {
Flavour had = p_softclusters->LowestTransition(flav1,flav2);
if (had.Mass()>sqrt(p_singlet->Mass2())) {
AddOrUpdateTransition(p_singlet, had);
......@@ -370,19 +371,17 @@ bool Singlet_Checker::BoostRecoilerInNewSystem(const Vec4D & newmom) {
void Singlet_Checker::ForcedDecays() {
//msg_Out()<<METHOD<<"("<<m_badones.size()<<" bad singlets):\n";
list<list<Singlet *>::iterator>::iterator bit=m_badones.begin();
while (bit!=m_badones.end()) {
p_singlet = (**bit);
if (ForcedDecayOfTwoPartonSinglet()) {
p_singlets->erase((*bit));
bit = m_badones.erase(bit);
//msg_Out()<<" forced decay of two parton singlet worked out.\n";
}
else {
Flavour flav1 = (**bit)->front()->Flavour();
Flavour flav2 = (**bit)->back()->Flavour();
Flavour had = p_softclusters->LowestTransition(flav1,flav2);
//Flavour flav1 = p_singlet->front()->Flavour();
//Flavour flav2 = p_singlet->back()->Flavour();
//Flavour had = p_softclusters->LowestTransition(flav1,flav2);
//msg_Out()<<METHOD<<": "<<flav1<<" + "<<flav2<<" --> "<<had<<" "
// <<"(mass = "<<sqrt((**bit)->Mass2())<<" from "
// <<(**bit)->front()->Momentum().Abs2()<<" and "
......
......@@ -34,15 +34,14 @@ void Double_Transitions::FillMap(Single_Transitions * singletransitions)
pair.second = pair2.second;
double weight = constituents->TotWeight(popped.Bar());
if (weight<1.e-6) continue;
if (2.*popped.HadMass()+0.1<pair.first.HadMass()+pair.second.HadMass())
if (2.*constituents->Mass(popped)+0.1<
constituents->Mass(pair.first)+constituents->Mass(pair.second))
weight = 1.;
if (popped.IsDiQuark()) {
if (pair.first.Kfcode()==4) weight *= m_charm_baryon_modifier;
if (pair.second.Kfcode()==4) weight *= m_charm_baryon_modifier;
if (pair.first.Kfcode()==5) weight *= m_beauty_baryon_modifier;
if (pair.second.Kfcode()==5) weight *= m_beauty_baryon_modifier;
//msg_Out()<<METHOD<<"["<<pair.first<<", "<<pair.second<<"] "
// <<"for popped diquark; weight = "<<weight<<".\n";
}
if (m_transitions.find(pair)==m_transitions.end())
m_transitions[pair] = new Double_Transition_List;
......@@ -83,7 +82,7 @@ void Double_Transitions::Print(const bool & full) {
<<METHOD<<":\n";
for (Double_Transition_Map::iterator dtmit=m_transitions.begin();
dtmit!=m_transitions.end();dtmit++) {
if (dtmit->first.first.Kfcode()!=4) continue;
if (dtmit->first.first.Kfcode()!=4 && dtmit->first.second.Kfcode()!=4) continue;
if (full)
msg_Out()<<"---------------------------------------------------------\n"
<<"*** ["<<dtmit->first.first<<" "<<dtmit->first.second<<"]:\n";
......
......@@ -120,7 +120,7 @@ ReadClusterDecayParameters(Data_Reader & dataread) {
// Probably irrelevant as long as they are small.
// We will probably not have to tune them.
m_parametermap[string("decay_threshold")] =
dataread.GetValue<double>("DECAY_THRESHOLD", 0.500);
dataread.GetValue<double>("DECAY_THRESHOLD", 0.000);
m_parametermap[string("piphoton_threshold")] =
dataread.GetValue<double>("PI_PHOTON_THRESHOLD", 0.150);
m_parametermap[string("dipion_threshold")] =
......
......@@ -21,7 +21,7 @@ double KT_Selector::operator()(const double & ktmax,const double & M2) {
m_sig2 = m_sigma2;
do {
kttest = ktmax*ran->Get();
} while (WeightFunction(kttest)<ran->Get());
} while (ran->Get() > WeightFunction(kttest));
return kttest;
}
......
......@@ -205,7 +205,7 @@ double Soft_Cluster_Handler::RadiationWeight() {
double m2(sit->first.Mass());
if (m2>m_mass) break;
// wave-function overlap * phase-space (units of 1 in total)
weight = sit->second * PhaseSpace(m2,0.);
weight = sit->second * PhaseSpace(m2,0.,false);
totweight += weights[sit->first] = weight;
}
double disc = totweight * ran->Get();
......@@ -244,7 +244,9 @@ double Soft_Cluster_Handler::DecayWeight() {
double m2(dit->first.first.Mass()), m3(dit->first.second.Mass());
if (m2+m3>m_mass) break;
// wave-function overlap * phase-space (units of 1 in total)
weight = dit->second * PhaseSpace(m2,m3);
bool heavy = (dit->first.first.IsB_Hadron() || dit->first.first.IsC_Hadron() ||
dit->first.second.IsB_Hadron() || dit->first.second.IsC_Hadron());
weight = dit->second * PhaseSpace(m2,m3,heavy);
totweight += weights[dit->first] = weight;
}
......@@ -349,7 +351,7 @@ DefineHadronsInAnnihilation(const Flavour_Pair & one,const Flavour_Pair & two) {
m3 = tit->first.Mass();
if (m2+m3>m_mass) break;
// wave-function overlap * phase-space (units of 1 in total)
weight = oit->second * tit->second * PhaseSpace(m2,m3);
weight = oit->second * tit->second * PhaseSpace(m2,m3,false);
Flavour_Pair flpair;
flpair.first = oit->first; flpair.second = tit->first;
totweight += weights[flpair] = weight;
......@@ -368,13 +370,13 @@ DefineHadronsInAnnihilation(const Flavour_Pair & one,const Flavour_Pair & two) {
return totweight;
}
double Soft_Cluster_Handler::PhaseSpace(const double & m2,const double & m3) {
double Soft_Cluster_Handler::
PhaseSpace(const double & m2,const double & m3,const bool heavyB) {
if (heavyB) return 1.;
double m22(m2*m2),m32(m3*m3);
double ps = sqrt(sqr(m_mass2-m22-m32)-4.*m22*m32)/(8.*M_PI*m_mass2);
// extra weight to possible steer away from phase space only ... may give
// preference to higher or lower mass pairs
//double mwt = (pow(m2/m_mass,m_chi) + pow(m3/m_mass,m_chi));
//double mwt = 1./(pow(m_mass/(m2+0.001),m_chi) + pow(m_mass/(m3+0.001),m_chi));
double mwt = pow(m2/m_mass,m_chi) + pow(m3/m_mass,m_chi);
return ps * mwt;
}
......
......@@ -29,7 +29,7 @@ namespace AHADIC {
void FillFlavours(Cluster * cluster);
int CheckOutsideRange();
int Decay();
double PhaseSpace(const double & m2,const double & m3);
double PhaseSpace(const double & m2,const double & m3,const bool heavyB);
double DecayWeight();
double RadiationWeight();
double Annihilation();
......
......@@ -178,16 +178,17 @@ bool Sherpa::InitializeTheEventHandler()
bool Sherpa::GenerateOneEvent(bool reset)
{
if (m_evt_output_start>0 && m_evt_output_start==rpa->gen.NumberOfGeneratedEvents()+1) {
if (m_evt_output_start>0 &&
m_evt_output_start==rpa->gen.NumberOfGeneratedEvents()+1) {
msg->SetLevel(m_evt_output);
}
if(m_debuginterval>0 && rpa->gen.NumberOfGeneratedEvents()%m_debuginterval==0){
if (p_inithandler->GetMatrixElementHandler()->SeedMode()!=3 ||
rpa->gen.NumberOfGeneratedEvents()==0) {
if(m_debuginterval>0 &&
rpa->gen.NumberOfGeneratedEvents()%m_debuginterval==0 &&
(p_inithandler->GetMatrixElementHandler()->SeedMode()!=3 ||
rpa->gen.NumberOfGeneratedEvents()==0)) {
std::string fname=ToString(rpa->gen.NumberOfGeneratedEvents())+".dat";
ran->WriteOutStatus(("random."+fname).c_str());
}
}
if (m_debugstep>=0) {
if (p_inithandler->GetMatrixElementHandler()->SeedMode()!=3)
......@@ -202,6 +203,7 @@ bool Sherpa::GenerateOneEvent(bool reset)
if (reset) p_eventhandler->Reset();
if (p_eventhandler->GenerateEvent(p_inithandler->Mode())) {
//msg_Out()<<"Generate event worked out.\n";
if(m_debuginterval>0 && rpa->gen.NumberOfGeneratedEvents()%m_debuginterval==0){
std::string fname=ToString(rpa->gen.NumberOfGeneratedEvents())+".dat";
std::ofstream eventout(("refevent."+fname).c_str());
......
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