Cluster_Splitter.C 7.55 KB
Newer Older
1 2 3 4 5 6 7
#include "AHADIC++/Decays/Cluster_Splitter.H"
#include "AHADIC++/Tools/Hadronisation_Parameters.H"
#include "ATOOLS/Org/Message.H"
#include "ATOOLS/Math/Random.H"

using namespace AHADIC;
using namespace ATOOLS;
8
using namespace std;
9

10
Cluster_Splitter::Cluster_Splitter(list<Cluster *> * cluster_list,
11 12
				   Soft_Cluster_Handler * softclusters) :
  Splitter_Base(cluster_list,softclusters),
13 14 15
  m_mode(1),
  m_output(false) 
{
16
  m_analyse = false;
17 18 19 20 21 22 23 24
  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);
  }
}
25 26 27 28 29 30 31 32 33 34 35 36

void Cluster_Splitter::Init(const bool & isgluon) {
  Splitter_Base::Init(false);
  m_alpha[0] = hadpars->Get("alphaL");
  m_beta[0]  = hadpars->Get("betaL");
  m_gamma[0] = hadpars->Get("gammaL");
  m_alpha[1] = hadpars->Get("alphaH");
  m_beta[1]  = hadpars->Get("betaH");
  m_gamma[1] = hadpars->Get("gammaH");
  m_alpha[2] = hadpars->Get("alphaD");
  m_beta[2]  = hadpars->Get("betaD");
  m_gamma[2] = hadpars->Get("gammaD");
37 38 39 40
  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"));
41 42 43 44
}

bool Cluster_Splitter::MakeLongitudinalMomenta() {
  CalculateLimits();
45 46 47 48 49 50 51 52 53 54 55
  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());
56
      m_R12   = m_m12min + DeltaM2(0); 
57
      FixCoefficients(p_part2->Flavour(),p_part1->Flavour());
58
      m_R21   = m_m22min + DeltaM2(1);
59 60 61 62 63 64 65 66
    } 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.;
  }
67 68 69
  return CheckIfAllowed();
}

70 71
void Cluster_Splitter::
FixCoefficients(const Flavour & flav1,const Flavour & flav2) {
72 73 74 75 76
  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()) 
77
    m_flcnt = 1;
78 79 80 81
  else if (p_part1->IsLeading() ||
	   ((m_mode==0) && p_part2->IsLeading())) {
    m_flcnt = 1;
    m_masses2 *= 4.;
82
  }
83
  else if (flav1.IsDiQuark())
84
    m_flcnt = 2;
85 86
  else if (p_part1->IsBeam() || p_part2->IsBeam())
    m_flcnt = 3;
87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102
  m_a = m_alpha[m_flcnt];
  m_b = m_beta[m_flcnt];
  m_c = m_gamma[m_flcnt];
}

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;
103 104 105 106 107 108 109 110
  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));
111 112
}

113 114 115 116
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);
117 118 119
  double deltaM, weight;
  do {
    deltaM = ran->Get()*deltaMmax;
120
    weight = pow(deltaM/sigmaM,m_a)*exp(-sqr(deltaM)/sigmaM);
121 122 123 124 125
  } while (weight<ran->Get()*wtmax);
  return sqr(deltaM);
}

bool Cluster_Splitter::CheckIfAllowed() {
126 127 128 129 130 131 132
  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()); 
133 134 135
  return allowed;
}

136 137
double Cluster_Splitter::
WeightFunction(const double & z,const double & zmin,const double & zmax) {
138 139
  double arg  = dabs(m_c)>1.e-2? m_c*(m_kt2+m_masses2)/m_kt02 : 0. ;
  double norm = exp(-arg/zmax); //*(1.-zmax)); 
140 141 142 143
  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);
144
    if (m_b<=0.) norm *= pow(1.-zmax,m_b);
145
  }
146
  double wt = pow(z,m_a) * pow(1.-z,m_b) * exp(-arg/z); //*(1.-z));
147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191
  if (wt>norm) {
    msg_Error()<<"Error in "<<METHOD<<": wt(z) = "<<wt<<"("<<z<<") "
	       <<"for wtmax = "<<norm<<" "
	       <<"[a, b, c = "<<m_a<<", "<<m_b<<", "<<m_c<<"].\n";
    exit(1);
  }
  return wt / norm;
}

bool Cluster_Splitter::FillParticlesInLists() {
  for (size_t i=0;i<2;i++) p_cluster_list->push_back(MakeCluster(i));
  return true;
}

Cluster * Cluster_Splitter::MakeCluster(size_t i) {
  double alpha = (i==0? m_z1  : 1.-m_z1 );
  double beta  = (i==0? m_z2  : 1.-m_z2 );
  double mass2 = (i==0? m_m12 : m_m22);
  double sign  = (i==0?    1. : -1.);
  double R2    = (i==0? m_R12 : m_R21);
  double R02   = mass2+(m_popped_mass2+m_kt2);
  double ab    = 4.*mass2*(m_popped_mass2+m_kt2);
  double x = 1;
  if (sqr(R2-R02)>ab) {
    double centre = (R2+mass2-(m_popped_mass2+m_kt2))/(2.*R2);
    double lambda = Lambda(R2,mass2,m_popped_mass2+m_kt2);
    x = (i==0)? centre+lambda : centre-lambda;
  }
  double y      = mass2/(x*R2);
  // This is the overall cluster momentum - we do not need it - and its
  // individual components, i.e. the momenta of the Proto_Particles
  // it is made of.
  Vec4D newmom11 = (m_E*(     x*alpha*s_AxisP+     y*(1.-beta)*s_AxisM));
  Vec4D newmom12 = (m_E*((1.-x)*alpha*s_AxisP+(1.-y)*(1.-beta)*s_AxisM) +
		    sign * m_ktvec);
  Vec4D clumom = m_E*(alpha*s_AxisP + (1.-beta)*s_AxisM) + sign * m_ktvec;

  // back into lab system
  m_rotat.RotateBack(newmom11);
  m_boost.BoostBack(newmom11);
  m_rotat.RotateBack(newmom12);
  m_boost.BoostBack(newmom12);
  (i==0?p_part1:p_part2)->SetMomentum(newmom11);

  Proto_Particle * newp =
192 193
    new Proto_Particle((i==0?m_newflav1:m_newflav2),newmom12,
		       false,p_part1->IsBeam()||p_part2->IsBeam());
194 195 196 197 198 199 200 201 202 203 204 205 206
  Cluster * cluster(i==0?new Cluster(p_part1,newp):new Cluster(newp,p_part2));
  /*if (m_output) {
  double mass = sqrt((newmom11+newmom12).Abs2());
  if (mass<m_popped_mass+(i==0?m_mass1:m_mass2)) {
    msg_Out()<<"Small masses(i="<<i<<") --> [";
    if (i==0) msg_Out()<<p_part1->Flavour()<<"/"<<m_newflav1;
    if (i==1) msg_Out()<<m_newflav2<<"/"<<p_part2->Flavour();
    msg_Out()<<"], m = "<<mass<<" vs. ";
    if (i==0) msg_Out()<<m_mass1<<"/"<<m_popped_mass;
    if (i==1) msg_Out()<<m_popped_mass<<"/"<<m_mass2;
    msg_Out()<<"\n"<<"    "<<(*cluster);
  }
  }*/
207 208 209 210 211 212 213 214 215
  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));
  }
216 217 218
  return cluster;
}