Commit c4aa4cd6 authored by Frank Krauss's avatar Frank Krauss

updated hadronisation, remnants, and MIs

parent f68a9da8
......@@ -5,13 +5,23 @@
using namespace AHADIC;
using namespace ATOOLS;
using namespace std;
Cluster_Splitter::Cluster_Splitter(std::list<Cluster *> * cluster_list,
Cluster_Splitter::Cluster_Splitter(list<Cluster *> * cluster_list,
Soft_Cluster_Handler * softclusters) :
Splitter_Base(cluster_list,softclusters),
m_mode(0),
m_output(false)
{}
m_mode(1),
m_output(false)
{
m_analyse = true;
if (m_analyse) {
m_histograms[string("kt")] = new Histogram(0,0.,5.,100);
m_histograms[string("z1")] = new Histogram(0,0.,1.,100);
m_histograms[string("z2")] = new Histogram(0,0.,1.,100);
m_histograms[string("mass")] = new Histogram(0,0.,100.,200);
m_histograms[string("Rmass")] = new Histogram(0,0.,2.,100);
}
}
void Cluster_Splitter::Init(const bool & isgluon) {
Splitter_Base::Init(false);
......@@ -24,6 +34,10 @@ void Cluster_Splitter::Init(const bool & isgluon) {
m_alpha[2] = hadpars->Get("alphaD");
m_beta[2] = hadpars->Get("betaD");
m_gamma[2] = hadpars->Get("gammaD");
m_alpha[3] = hadpars->Get("alphaB");
m_beta[3] = hadpars->Get("betaB");
m_gamma[3] = hadpars->Get("gammaB");
m_kt02 = sqr(hadpars->Get("kt_o"));
}
bool Cluster_Splitter::MakeLongitudinalMomenta() {
......@@ -39,9 +53,9 @@ bool Cluster_Splitter::MakeLongitudinalMomenta() {
else if (m_mode==1) {
do {
FixCoefficients(p_part1->Flavour(),p_part2->Flavour());
m_R12 = m_m12min + DeltaM2();
m_R12 = m_m12min + DeltaM2(0);
FixCoefficients(p_part2->Flavour(),p_part1->Flavour());
m_R21 = m_m22min + DeltaM2();
m_R21 = m_m22min + DeltaM2(1);
} 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;
......@@ -55,23 +69,24 @@ bool Cluster_Splitter::MakeLongitudinalMomenta() {
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 = 0;
m_masses2 = sqr(Max(1.,(p_constituents->Mass(flav1) +
p_constituents->Mass(flav2))));
if (flav1==Flavour(kf_b) || flav1==Flavour(kf_b).Bar() ||
flav1==Flavour(kf_c) || flav1==Flavour(kf_c).Bar())
m_flcnt = 1;
else if (p_part1->IsLeading() ||
((m_mode==0) && p_part2->IsLeading())) {
m_flcnt = 1;
m_masses2 *= 4.;
}
else if (flav1.IsDiQuark()) {
else if (flav1.IsDiQuark())
m_flcnt = 2;
}
else if (p_part1->IsBeam() || p_part2->IsBeam())
m_flcnt = 3;
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() {
......@@ -85,42 +100,50 @@ void Cluster_Splitter::CalculateLimits() {
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);
m_z1min = (centre1-lambda)/(2.*m_Q2);
m_z1max = (centre1+lambda)/(2.*m_Q2);
m_z2min = (centre2-lambda)/(2.*m_Q2);
m_z2max = (centre2+lambda)/(2.*m_Q2);
m_mean12 = (m_minQ_1+m_maxQ_1)/2.;
m_mean21 = (m_minQ_2+m_maxQ_2)/2.;
m_sigma12 = sqr(Max(1.,m_maxQ_1-m_minQ_1));
m_sigma21 = sqr(Max(1.,m_maxQ_2-m_minQ_2));
}
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 Cluster_Splitter::DeltaM2(const size_t & cl) {
double deltaMmax = Min(sqrt(m_Q2-m_m12min-m_m22min),10.+(cl==0?m_maxQ_1:m_maxQ_2));
double sigmaM = m_c * (cl==0?m_sigma12:m_sigma21);
double wtmax = pow(deltaMmax/sigmaM,m_a);
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);
weight = pow(deltaM/sigmaM,m_a)*exp(-sqr(deltaM)/sigmaM);
} 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));
bool allowed = (m_R12>m_minQ_12 && m_R21>m_minQ_22);
//if (p_part1->IsLeading() && m_R12>sqr(m_maxQ_1))
// allowed = allowed &&
// (exp(-sqr(sqrt(m_R12)-m_maxQ_1)/(2.*m_sigma12))>ran->Get());
//if (p_part2->IsLeading() && m_R21>sqr(m_maxQ_2))
// allowed = allowed &&
// (exp(-sqr(sqrt(m_R21)-m_maxQ_2)/(2.*m_sigma21))>ran->Get());
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. ;
double norm = exp(-arg/zmax);
double arg = dabs(m_c)>1.e-2? m_c*(m_kt2+m_masses2)/m_kt02 : 0. ;
double norm = exp(-arg/zmax); //*(1.-zmax));
if (m_a<=0. && m_b<=0.)
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);
}
double wt = pow(z,m_a) * pow(1.-z,m_b) * exp(-arg/z);
double wt = pow(z,m_a) * pow(1.-z,m_b) * exp(-arg/z); //*(1.-z));
if (wt>norm) {
msg_Error()<<"Error in "<<METHOD<<": wt(z) = "<<wt<<"("<<z<<") "
<<"for wtmax = "<<norm<<" "
......@@ -166,7 +189,8 @@ Cluster * Cluster_Splitter::MakeCluster(size_t i) {
(i==0?p_part1:p_part2)->SetMomentum(newmom11);
Proto_Particle * newp =
new Proto_Particle((i==0?m_newflav1:m_newflav2),newmom12);
new Proto_Particle((i==0?m_newflav1:m_newflav2),newmom12,
false,p_part1->IsBeam()||p_part2->IsBeam());
Cluster * cluster(i==0?new Cluster(p_part1,newp):new Cluster(newp,p_part2));
/*if (m_output) {
double mass = sqrt((newmom11+newmom12).Abs2());
......@@ -180,6 +204,15 @@ Cluster * Cluster_Splitter::MakeCluster(size_t i) {
msg_Out()<<"\n"<<" "<<(*cluster);
}
}*/
if (m_analyse && m_Q>91.) {
if (i==1) {
m_histograms[string("kt")]->Insert(sqrt(m_kt2));
m_histograms[string("z1")]->Insert(m_z1);
m_histograms[string("z2")]->Insert(m_z2);
}
m_histograms[string("mass")]->Insert(sqrt(R2));
m_histograms[string("Rmass")]->Insert(2.*sqrt(R2/m_Q2));
}
return cluster;
}
......@@ -7,31 +7,30 @@ 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;
double m_alpha[4], m_beta[4], m_gamma[4], m_a, m_b, m_c, m_kt02;
size_t m_flcnt;
double m_R12, m_R21;
double m_m12min, m_m22min, m_masses2;
double m_mean12, m_mean21, m_sigma12, m_sigma21, m_sigma;
bool m_output;
bool m_output;
bool MakeLongitudinalMomenta();
double DeltaM2();
void FixCoefficients(const ATOOLS::Flavour & flav1,
const ATOOLS::Flavour & flav2);
void CalculateLimits();
bool FillParticlesInLists();
bool CheckIfAllowed();
bool CheckKinematics() { return true; }
bool MakeLongitudinalMomenta();
double DeltaM2(const size_t & cl);
void FixCoefficients(const ATOOLS::Flavour & flav1,
const ATOOLS::Flavour & flav2);
void CalculateLimits();
bool FillParticlesInLists();
bool CheckIfAllowed();
bool CheckKinematics() { return true; }
Cluster * MakeCluster(size_t i);
public:
Cluster_Splitter(std::list<Cluster *> * cluster_list,
Soft_Cluster_Handler * softclusters);
void Init(const bool & isgluon=false);
void Init(const bool & isgluon=false);
double WeightFunction(const double & z,
const double & zmin,const double & zax);
void SetOutput(const bool & out) { m_output = out; }
void SetOutput(const bool & out) { m_output = out; }
};
}
......
......@@ -101,7 +101,7 @@ bool Gluon_Decayer::SplitGluonRing() {
Proto_Particle * Gluon_Decayer::FirstGluon() {
double minm2(1.e12), m2thres(sqr(2.*m_breaker.MinMass())), m2;
list<Proto_Particle *>::iterator ppiter1, ppiter2, winner;
list<Proto_Particle *>::iterator ppiter1, ppiter2, winner(p_singlet->end());
for (list<Proto_Particle *>::iterator ppiter1=p_singlet->begin();
ppiter1!=p_singlet->end();ppiter1++) {
Proto_Particle * part1(*ppiter1);
......@@ -115,7 +115,7 @@ Proto_Particle * Gluon_Decayer::FirstGluon() {
winner = ppiter1;
}
}
return (*winner);
return (winner==p_singlet->end()?(*p_singlet->begin()):(*winner));
}
int Gluon_Decayer::Step(Proto_Particle * part1,Proto_Particle * part2,
......
......@@ -164,7 +164,8 @@ Cluster * Gluon_Splitter::MakeCluster() {
// Update momentum of original (anti-)(di-) quark after gluon splitting
p_part1->SetMomentum(newmom11);
// Make new particle
Proto_Particle * newp12 = new Proto_Particle(m_newflav1,newmom12,'l');
Proto_Particle * newp12 = new Proto_Particle(m_newflav1,newmom12,false,
p_part1->IsBeam() || p_part2->IsBeam());
// Take care of sequence in cluster = triplet + anti-triplet
Cluster * cluster(m_barrd?new Cluster(newp12,p_part1):new Cluster(p_part1,newp12));
return cluster;
......
......@@ -115,7 +115,7 @@ bool Singlet_Checker::CheckSinglet() {
// Checking the mass for pairs of colour-connected particles
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) {
if ((*plit)->Momentum()[0]<0. || (*plit)->Momentum().Abs2()<-1.e-2) {
msg_Error()<<"Error in "<<METHOD<<":\n"
<<" negative energy or mass^2 particle in singlet:\n"
<<(*p_singlet)<<"n";
......
......@@ -48,6 +48,7 @@ Singlet * Singlet_Former::MakeAnother() {
Particle * part = FindStart();
partlist->push_back(new Proto_Particle(*part));
if (part->Flav().IsQuark()) partlist->back()->SetLeading(true);
if (part->Beam()>-1) partlist->back()->SetBeam(true);
unsigned int col1 = part->GetFlow(1);
unsigned int col2 = part->GetFlow(2);
while (col2!=col1) {
......@@ -59,6 +60,7 @@ Singlet * Singlet_Former::MakeAnother() {
col1 = part->GetFlow(1);
partlist->push_back(new Proto_Particle(*part));
if (part->Flav().IsQuark()) partlist->back()->SetLeading(true);
if (part->Beam()>-1) partlist->back()->SetBeam(true);
break;
}
}
......
......@@ -51,7 +51,7 @@ bool Trivial_Splitter::operator()(Singlet * singlet) {
p_part1->SetMomentum(m_q1mom);
p_part2->SetMomentum(m_glumom);
p_singlet->push_back(new Proto_Particle(m_newflav.Bar(),m_q2mom,'l'));
p_singlet->push_back(new Proto_Particle(m_newflav.Bar(),m_q2mom));
return true;
}
......
......@@ -13,7 +13,7 @@ using namespace ATOOLS;
using namespace std;
Ahadic::Ahadic(string path, string file) :
Ahadic::Ahadic(string path, string file,string shower) :
m_softclusters(Soft_Cluster_Handler(&m_hadron_list)),
m_beamparticles(Beam_Particles_Shifter(&m_singlet_list, &m_softclusters)),
m_sformer(Singlet_Former(&m_singlet_list)),
......@@ -22,7 +22,7 @@ Ahadic::Ahadic(string path, string file) :
m_clusterdecayer(Cluster_Decayer(&m_cluster_list, &m_softclusters))
{
hadpars = new Hadronisation_Parameters();
hadpars->Init(path,file);
hadpars->Init(path,file,shower);
m_beamparticles.Init();
m_softclusters.Init();
......
......@@ -44,7 +44,8 @@ namespace AHADIC {
void CleanUp(ATOOLS::Blob * blob);
public:
Ahadic(std::string path=std::string("./"),
std::string file=std::string("Cluster.dat"));
std::string file=std::string("Cluster.dat"),
std::string shower="Dire");
~Ahadic();
ATOOLS::Return_Value::code Hadronize(ATOOLS::Blob_List *);
......
......@@ -38,10 +38,10 @@ void Double_Transitions::FillMap(Single_Transitions * singletransitions)
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;
if (int(pair.first.Kfcode())==4) weight *= m_charm_baryon_modifier;
if (int(pair.second.Kfcode())==4) weight *= m_charm_baryon_modifier;
if (int(pair.first.Kfcode())==5) weight *= m_beauty_baryon_modifier;
if (int(pair.second.Kfcode())==5) weight *= m_beauty_baryon_modifier;
}
if (m_transitions.find(pair)==m_transitions.end())
m_transitions[pair] = new Double_Transition_List;
......
This diff is collapsed.
......@@ -26,6 +26,7 @@ namespace AHADIC {
class Hadronisation_Parameters {
private:
int m_shower;
double m_offset;
Constituents * p_constituents;
Single_Transitions * p_stransitions;
......@@ -42,13 +43,12 @@ namespace AHADIC {
void ReadMesonWeights(ATOOLS::Data_Reader & dataread);
void ReadGluonSplittingParameters(ATOOLS::Data_Reader & dataread);
void ReadClusterDecayParameters(ATOOLS::Data_Reader & dataread);
void ReadClusterToMesonParameters(ATOOLS::Data_Reader & dataread);
void ReadDeprecatedParameters(ATOOLS::Data_Reader & dataread);
void ReadClusterToMesonPSParameters(ATOOLS::Data_Reader & dataread);
public:
Hadronisation_Parameters();
~Hadronisation_Parameters();
void Init(std::string,std::string);
void Init(std::string,std::string,std::string);
double Get(std::string keyword);
const bool & AnaOn() const { return m_ana; }
......
......@@ -18,10 +18,10 @@ void KT_Selector::Init(const bool & isgluon) {
double KT_Selector::operator()(const double & ktmax,const double & M2) {
double kttest(-1.);
m_sig2 = m_sigma2;
m_sig2 = M2*m_sigma2;
do {
kttest = ktmax*ran->Get();
} while (ran->Get() > WeightFunction(kttest));
} while (ran->Get() > WeightFunction(kttest));
return kttest;
}
......
......@@ -14,7 +14,7 @@ namespace AHADIC {
~KT_Selector();
void Init(const bool & isgluon);
double operator()(const double & ktmax,const double & M2=-1.);
double operator()(const double & ktmax,const double & M2=1.);
};
}
......
......@@ -9,7 +9,7 @@ using namespace ATOOLS;
using namespace std;
Soft_Cluster_Handler::Soft_Cluster_Handler(list<Proto_Particle *> * hadrons) :
p_hadrons(hadrons)
p_hadrons(hadrons), m_ktfac(1.)
{ }
Soft_Cluster_Handler::~Soft_Cluster_Handler()
......@@ -120,7 +120,7 @@ bool Soft_Cluster_Handler::TreatSingletCluster() {
}
// below two-pion threshold
else {
if (ran->Get()>0.5) {
if (ran->Get()>0.66) {
m_hads[0] = m_hads[1] = Flavour(kf_pi);
}
else {
......@@ -169,14 +169,15 @@ bool Soft_Cluster_Handler::FixKinematics() {
double M2(m_mass*m_mass);
double m12(sqr(m_hads[0].Mass())),m22(sqr(m_hads[1].Mass()));
double E1((M2+m12-m22)/(2.*m_mass)), p1(sqrt(sqr(E1)-m12));
//double cth = 1.-ran->Get(); // only between 0 and 1.
//double sth = sqrt(1.-cth*cth);
double pt = m_ktselector(p1);
double pl = sqrt(p1*p1-pt*pt);
double phi = 2.*M_PI*ran->Get();
m_moms[0] = Vec4D( E1, pt*cos(phi), pt*sin(phi), pl);
m_moms[1] = Vec4D(m_mass-E1,-pt*cos(phi),-pt*sin(phi),-pl);
double E1((M2+m12-m22)/(2.*m_mass)), p1(sqrt(sqr(E1)-m12)), pt;
bool lead = (*p_cluster)[0]->IsLeading() || (*p_cluster)[0]->IsLeading();
//if (lead)
pt = m_ktselector(p1,lead?1.:m_ktfac);
//else pt = p1 * sqrt(1.-sqr(1.-2.*ran->Get()));
double pl = sqrt(p1*p1-pt*pt);
double phi = 2.*M_PI*ran->Get();
m_moms[0] = Vec4D( E1, pt*cos(phi), pt*sin(phi), pl);
m_moms[1] = Vec4D(m_mass-E1,-pt*cos(phi),-pt*sin(phi),-pl);
bool print(false);
for (size_t i=0;i<2;i++) {
rotat.RotateBack(m_moms[i]);
......@@ -372,7 +373,7 @@ DefineHadronsInAnnihilation(const Flavour_Pair & one,const Flavour_Pair & two) {
double Soft_Cluster_Handler::
PhaseSpace(const double & m2,const double & m3,const bool heavyB) {
if (heavyB) return 1.;
if (m_chi<0. || 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
......
......@@ -14,7 +14,7 @@ namespace AHADIC {
Single_Transitions * p_singletransitions;
Double_Transitions * p_doubletransitions;
KT_Selector m_ktselector;
double m_light, m_chi;
double m_light, m_chi, m_ktfac;
double m_piphoton_threshold, m_dipion_threshold, m_open_threshold;
Cluster * p_cluster;
......
......@@ -10,9 +10,23 @@ using namespace std;
Splitter_Base::Splitter_Base(list<Cluster *> * cluster_list,
Soft_Cluster_Handler * softclusters) :
p_cluster_list(cluster_list), p_softclusters(softclusters),
m_attempts(100)
m_ktfac(1.),
m_attempts(100), m_analyse(false)
{}
Splitter_Base::~Splitter_Base() {
Histogram * histo;
string name;
for (map<string,Histogram *>::iterator hit=m_histograms.begin();
hit!=m_histograms.end();hit++) {
histo = hit->second;
name = string("Fragmentation_Analysis/")+hit->first+string(".dat");
histo->Output(name);
delete histo;
}
m_histograms.clear();
}
void Splitter_Base::Init(const bool & isgluon) {
p_singletransitions = hadpars->GetSingleTransitions();
p_doubletransitions = hadpars->GetDoubleTransitions();
......@@ -117,9 +131,10 @@ void Splitter_Base::PopFlavours() {
void Splitter_Base::DetermineMinimalMasses() {
if (!m_flavs1.first.IsGluon() && !m_flavs1.second.IsGluon()) {
m_minQ_1 = p_doubletransitions->GetLightestMass(m_flavs1);
if (!m_flavs1.first.IsDiQuark() && !m_flavs1.second.IsDiQuark())
m_minQ_1 = Min(m_minQ_1,
Max(0.,p_singletransitions->GetLightestMass(m_flavs1)));
m_maxQ_1 = p_doubletransitions->GetHeaviestMass(m_flavs1);
if (!(m_flavs1.first.IsDiQuark() && m_flavs1.second.IsDiQuark()))
m_minQ_1 = Min(m_minQ_1,
Max(0.,p_singletransitions->GetLightestMass(m_flavs1)));
}
else {
m_minQ_1 = (p_constituents->Mass(m_flavs1.first)+
......@@ -127,18 +142,17 @@ void Splitter_Base::DetermineMinimalMasses() {
}
if (!m_flavs2.first.IsGluon() && !m_flavs2.second.IsGluon()) {
m_minQ_2 = p_doubletransitions->GetLightestMass(m_flavs2);
if (!m_flavs2.first.IsDiQuark() && !m_flavs2.second.IsDiQuark())
m_minQ_2 = Min(m_minQ_2,
Max(0.,p_singletransitions->GetLightestMass(m_flavs2)));
m_maxQ_2 = p_doubletransitions->GetHeaviestMass(m_flavs2);
if (!(m_flavs2.first.IsDiQuark() && m_flavs2.second.IsDiQuark()))
m_minQ_2 = Min(m_minQ_2,
Max(0.,p_singletransitions->GetLightestMass(m_flavs2)));
}
else {
m_minQ_2 = (p_constituents->Mass(m_flavs2.first)+
p_constituents->Mass(m_flavs2.second));
}
//msg_Out()<<METHOD<<"("<<m_flavs1.first<<"+"<<m_flavs1.second<<" --> "<<m_minQ_1<<", "
// <<m_flavs2.first<<"+"<<m_flavs2.second<<" --> "<<m_minQ_2<<").\n";
m_minQ_12 = sqr(m_minQ_1);
m_minQ_22 = sqr(m_minQ_2);
m_minQ_12 = sqr(m_minQ_1); m_maxQ_12 = sqr(m_minQ_1);
m_minQ_22 = sqr(m_minQ_2); m_maxQ_22 = sqr(m_maxQ_2);
}
bool Splitter_Base::MakeKinematics() {
......@@ -148,13 +162,16 @@ bool Splitter_Base::MakeKinematics() {
void Splitter_Base::MakeTransverseMomentum() {
m_ktmax = (m_Emax-2.*m_popped_mass)/2.;
// have to make this a parameter for the beam breakup?
if (p_part1->IsBeam() || p_part2->IsBeam()) m_ktmax = Min(2.0,m_ktmax);
if (m_ktmax<0.) {
msg_Error()<<METHOD<<" yields error ktmax = "<<m_ktmax
<<" from "<<m_Emax<<", "<<m_popped_mass<<" vs. "
<<" min = "<<m_minmass<<".\n";
abort();
}
m_kt = m_ktselector(m_ktmax,sqr(m_Q-m_mass1-m_mass2));
bool islead = p_part1->IsLeading() || p_part2->IsLeading();
m_kt = m_ktselector(m_ktmax,islead || p_part2->Flavour().IsGluon()?1.:m_ktfac);
m_kt2 = m_kt*m_kt;
m_phi = 2.*M_PI*ran->Get();
m_ktvec = m_kt * Vec4D(0.,cos(m_phi),sin(m_phi),0.);
......
......@@ -6,6 +6,8 @@
#include "AHADIC++/Tools/Z_Selector.H"
#include "AHADIC++/Tools/Soft_Cluster_Handler.H"
#include "ATOOLS/Math/Poincare.H"
#include "ATOOLS/Math/Histogram.H"
#include <map>
namespace AHADIC {
class Splitter_Base {
......@@ -32,13 +34,17 @@ namespace AHADIC {
double m_popped_mass, m_popped_mass2;
ATOOLS::Flavour m_newflav1, m_newflav2;
double m_minQ_1, m_minQ_2, m_minQ_12, m_minQ_22;
double m_maxQ_1, m_maxQ_2, m_maxQ_12, m_maxQ_22;
double m_z1min, m_z1max, m_z2min, m_z2max, m_z1, m_z2;
double m_ktmax, m_kt, m_kt2, m_phi;
double m_ktmax, m_ktfac, m_kt, m_kt2, m_phi;
ATOOLS::Vec4D m_ktvec;
ATOOLS::Poincare m_boost, m_rotat;
bool m_analyse;
std::map<std::string,ATOOLS::Histogram *> m_histograms;
virtual bool InitSplitting(Proto_Particle * part1,Proto_Particle * part2,
Proto_Particle * part3);
virtual void FillMasses();
......@@ -56,6 +62,7 @@ namespace AHADIC {
public:
Splitter_Base(std::list<Cluster *> * cluster_list,
Soft_Cluster_Handler * softclusters);
~Splitter_Base();
virtual void Init(const bool & isgluon);
virtual bool operator()(Proto_Particle * part1,Proto_Particle * part2,
......
......@@ -170,6 +170,8 @@ void Amisic::Reset() {}
void Amisic::InitAnalysis() {
m_nscatters = 0;
m_histos[string("N_scatters")] = new Histogram(0,0,20,20);
m_histos[string("B")] = new Histogram(0,0,10,100);
m_histos[string("Bfac")] = new Histogram(0,0,10,100);
m_histos[string("P_T(1)")] = new Histogram(0,0,100,100);
m_histos[string("Y(1)")] = new Histogram(0,-10,10,10);
m_histos[string("Delta_Y(1)")] = new Histogram(0,0,10,10);
......@@ -212,6 +214,8 @@ void Amisic::Analyse(const bool & last) {
m_nscatters++;
if (last) {
m_histos[string("N_scatters")]->Insert(double(m_nscatters)+0.5);
m_histos[string("B")]->Insert(m_b);
m_histos[string("Bfac")]->Insert(m_bfac);
m_nscatters = 0;
}
}
......
......@@ -18,7 +18,10 @@ using namespace std;
MI_Process_Group::MI_Process_Group(const std::string & name) :
m_name(name), m_lastxs(0.), m_pref(M_PI)
{
m_pt02 = sqr((*mipars)("pt_0"));
m_muR_scheme = mipars->GetScaleScheme();
m_muR_fac = sqr((*mipars)("RenScale_Factor"));
m_muF_fac = sqr((*mipars)("FacScale_Factor"));
m_pt02 = sqr((*mipars)("pt_0"));
}
MI_Process_Group::~MI_Process_Group() {
......@@ -54,7 +57,13 @@ double MI_Process_Group::SoftCorrection(const double & pt2) const {
double MI_Process_Group::Scale(const double & pt2) const {
// Default scale, including an IR regularisation - maybe we should get more choices.
return pt2+m_pt02;
switch (m_muR_scheme) {
case scale_scheme::PT_with_Raps:
exit(1);
case scale_scheme::PT:
default:
return (pt2+m_pt02);
}
}
void MI_Process_Group::Output() const {
......@@ -116,7 +125,7 @@ operator()(const double & shat,const double & that,const double & uhat) {
double pref = m_pref/sqr(shat);
Flavour gluon(kf_gluon);
double pdf = p_pdf[0]->GetXPDF(gluon)*p_pdf[1]->GetXPDF(gluon);
double cpl = sqr((*p_alphaS)(Scale(m_scale)));
double cpl = sqr((*p_alphaS)(m_muR_fac*Scale(m_scale)));
double tot = 0.;
for (list<MI_Process * >::iterator mit=m_processes.begin();
mit!=m_processes.end();mit++) tot += (**mit)();
......@@ -198,7 +207,7 @@ double MI_QQB_Processes::
operator()(const double & shat,const double & that,const double & uhat) {
double tot = PreCalculate(shat,that,uhat);
double pref = m_pref/sqr(shat);
double cpl = sqr((*p_alphaS)(Scale(m_scale)));
double cpl = sqr((*p_alphaS)(m_muR_fac*Scale(m_scale)));
double pdf = 0;
for (size_t i=1;i<6;i++) {
Flavour flav = Flavour(i);
......@@ -236,7 +245,7 @@ double MI_QQ_Processes::
operator()(const double & shat,const double & that,const double & uhat) {
double tot = PreCalculate(shat,that,uhat);
double pref = m_pref/sqr(shat);
double cpl = sqr((*p_alphaS)(Scale(m_scale)));
double cpl = sqr((*p_alphaS)(m_muR_fac*Scale(m_scale)));
double pdf = 0;
for (size_t i=1;i<6;i++) {
Flavour flav = Flavour(i);
......@@ -300,7 +309,7 @@ double MI_QG_Processes::
operator()(const double & shat,const double & that,const double & uhat) {
double tot = PreCalculate(shat,that,uhat);
double pref = m_pref/sqr(shat);
double cpl = sqr((*p_alphaS)(Scale(m_scale)));
double cpl = sqr((*p_alphaS)(m_muR_fac*Scale(m_scale)));
double pdf = 0;
Flavour gluon(kf_gluon);
for (size_t i=1;i<6;i++) {
......@@ -370,7 +379,7 @@ double MI_Q1Q2_Processes::
operator()(const double & shat,const double & that,const double & uhat) {
double tot = PreCalculate(shat,that,uhat);
double pref = m_pref/sqr(shat);
double cpl = sqr((*p_alphaS)(Scale(m_scale)));
double cpl = sqr((*p_alphaS)(m_muR_fac*Scale(m_scale)));
double pdf = 0;
for (size_t i=1;i<6;i++) {
Flavour flav1 = Flavour(i);
......@@ -446,7 +455,7 @@ double MI_QG_QGamma_Processes::
operator()(const double & shat,const double & that,const double & uhat) {
double me2 = PreCalculate(shat,that,uhat);
double pref = m_pref/sqr(shat);
double cpl = (*p_alphaS)(Scale(m_scale)) * (*p_alpha)(Scale(m_scale));
double cpl = (*p_alphaS)(m_muR_fac*Scale(m_scale)) * (*p_alpha)(Scale(m_scale));
double pdf = 0;
Flavour gluon(kf_gluon);
for (size_t i=1;i<6;i++) {
......@@ -510,7 +519,7 @@ double MI_QQ_GGamma_Processes::
operator()(const double & shat,const double & that,const double & uhat) {
double me2 = PreCalculate(shat,that,uhat);
double pref = m_pref/sqr(shat);
double cpl = (*p_alphaS)(Scale(m_scale)) * (*p_alpha)(Scale(m_scale));
double cpl = (*p_alphaS)(m_muR_fac*Scale(m_scale)) * (*p_alpha)(Scale(m_scale));
double pdf = 0;
for (size_t i=1;i<6;i++) {
Flavour quark = Flavour(i);
......
......@@ -2,6 +2,7 @@
#define AMISIC_Perturbative_MI_Process_Group_H
#include "AMISIC++/Perturbative/MI_Process.H"
#include "AMISIC++/Tools/MI_Parameters.H"
#include "MODEL/Main/Model_Base.H"
#include "MODEL/Main/Running_AlphaQED.H"
#include "MODEL/Main/Running_AlphaS.H"
......@@ -11,8 +12,10 @@
namespace AMISIC {