Commit 613cf074 authored by Frank Krauss's avatar Frank Krauss

moving new AHADIC into trunk - finished

parent f38b20a0
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -50,11 +50,6 @@ Cluster_Decay_Analysis::~Cluster_Decay_Analysis()
return;
Histogram * histo;
string name;
// Sync histos in case of MPI (even on one core, otherwise empty histos)
#ifdef USING__MPI
for (map<string,Histogram *>::iterator hit=m_histograms.begin();
hit!=m_histograms.end();hit++) hit->second->MPISync();
#endif
for (map<string,Histogram *>::iterator hit=m_histograms.begin();
hit!=m_histograms.end();hit++) {
histo = hit->second;
......
#include "AHADIC++/Decays/Cluster_Decay_Handler.H"
#include "AHADIC++/Decays/Cluster_Part.H"
#include "AHADIC++/Tools/Hadronisation_Parameters.H"
using namespace AHADIC;
using namespace ATOOLS;
using namespace std;
Cluster_Decay_Handler::Cluster_Decay_Handler(Cluster_List * clulist,bool ana) :
p_softclusters(hadpars->GetSoftClusterHandler()),
p_clus(new Cluster_Part(ana)),
p_clulist(clulist),
p_analysis(ana?new Cluster_Decay_Analysis():NULL)
{ }
Cluster_Decay_Handler::~Cluster_Decay_Handler()
{
if (p_clus) { delete p_clus; p_clus=NULL; }
if (p_analysis) { delete p_analysis; p_analysis=NULL; }
}
int Cluster_Decay_Handler::DecayClusters(Blob * blob)
{
Cluster * cluster;
Cluster_Iterator cit(p_clulist->begin());
//msg_Out()<<":::::: "<<METHOD<<" with "<<p_clulist->size()<<" clusters.\n";
while (!p_clulist->empty()) {
cluster = p_clulist->front();
if (!cluster->Active()) return -1;
if (p_clus->TestDecay(cluster)) {
//msg_Out()<<":::::: "<<METHOD<<": decay ok for\n"<<(*cluster);
Cluster_List * clist(cluster->GetClusters());
if (!p_softclusters->TreatClusterList(clist,blob)) {
msg_Error()<<"Error in "<<METHOD<<" : \n"
<<" Did not find a kinematically allowed "
<<"solution for the cluster list of length " << clist->size() <<".\n"
<<" Will trigger retrying the event.\n";
return -1;
}
while (!clist->empty()) {
p_clulist->push_back(clist->front());
clist->pop_front();
}
//msg_Out()<<":::::: "<<p_clulist->size()<<" clusters now.\n";
}
else {
//msg_Out()<<":::::: "<<METHOD<<":\n"
// <<" Enter soft cluster treatment for undecayed cluster.\n"
// <<(*cluster);
Cluster_List clist;
clist.push_back(cluster);
if (!p_softclusters->TreatClusterList(&clist,blob)) {
msg_Error()<<"Error in "<<METHOD<<".\n";
return -1;
}
}
delete (p_clulist->front()->GetTrip());
delete (p_clulist->front()->GetAnti());
delete (p_clulist->front());
p_clulist->pop_front();
}
if (p_analysis) p_analysis->AnalyseThis(blob);
return 1;
}
#ifndef AHADIC_Decays_Cluster_Decay_Handler_H
#define AHADIC_Decays_Cluster_Decay_Handler_H
#include "AHADIC++/Decays/Cluster_Part.H"
#include "AHADIC++/Tools/Soft_Cluster_Handler.H"
#include "ATOOLS/Phys/Blob_List.H"
#include "AHADIC++/Tools/Cluster.H"
#include "AHADIC++/Decays/Cluster_Decay_Analysis.H"
namespace AHADIC {
class Cluster_Decay_Handler {
private:
Soft_Cluster_Handler * p_softclusters;
Cluster_Part * p_clus;
Cluster_List * p_clulist;
Cluster_Decay_Analysis * p_analysis;
ATOOLS::Blob * DecayIt(Cluster *);
public:
Cluster_Decay_Handler(Cluster_List *,bool=false);
~Cluster_Decay_Handler();
int DecayClusters(ATOOLS::Blob *);
};
}
#endif
#include "AHADIC++/Decays/Cluster_Decayer.H"
#include "AHADIC++/Tools/Hadronisation_Parameters.H"
#include "ATOOLS/Math/Random.H"
#include "ATOOLS/Org/Message.H"
using namespace AHADIC;
using namespace ATOOLS;
using namespace std;
Cluster_Decayer::Cluster_Decayer(list<Cluster *> * cluster_list,
Soft_Cluster_Handler * softclusters) :
p_cluster_list(cluster_list), p_softclusters(softclusters),
m_splitter(Cluster_Splitter(cluster_list,softclusters))
{}
Cluster_Decayer::~Cluster_Decayer() {}
void Cluster_Decayer::Init() { m_splitter.Init(); }
void Cluster_Decayer::Reset() {}
bool Cluster_Decayer::operator()() {
while (!p_cluster_list->empty()) {
if (!Treat(p_cluster_list->front())) {
return false;
}
p_cluster_list->pop_front();
}
return true;
}
bool Cluster_Decayer::Treat(Cluster * cluster) {
bool mustdecay = p_softclusters->MustPromptDecay(cluster);
if (!mustdecay && m_splitter((*cluster)[0],(*cluster)[1])) {
delete cluster;
return true;
}
switch (p_softclusters->Treat(cluster,true)) {
case -1:
// cluster cannot decay into anything - return false (triggers new event)
cluster->Clear();
delete cluster;
return false;
case 1:
// cluster decayed into hadrons - delete it and carry on.
cluster->Clear();
delete cluster;
return true;
case 0:
default:
//cluster should have decayed into clusters - throw error
break;
}
msg_Error()<<METHOD<<" throws error for:\n"<<(*cluster)<<"\n";
return false;
}
#ifndef AHADIC_Decays_Gluon_Decayer_H
#define AHADIC_Decays_Gluon_Decayer_H
#include "AHADIC++/Decays/Cluster_Splitter.H"
#include "AHADIC++/Tools/Soft_Cluster_Handler.H"
#include "AHADIC++/Tools/Proto_Particle.H"
#include <list>
namespace AHADIC {
class Cluster_Decayer {
private:
std::list<Cluster *> * p_cluster_list;
Soft_Cluster_Handler * p_softclusters;
Cluster_Splitter m_splitter;
bool Treat(Cluster * cluster);
public:
Cluster_Decayer(std::list<Cluster *> * cluster_list,
Soft_Cluster_Handler * softclusters);
~Cluster_Decayer();
void Init();
void Reset();
bool operator()();
};
}
#endif
#include "AHADIC++/Decays/Cluster_Part.H"
#include "ATOOLS/Org/Message.H"
#include "ATOOLS/Math/Random.H"
using namespace AHADIC;
using namespace ATOOLS;
using namespace std;
Cluster_Part::Cluster_Part(bool ana) :
m_ana(ana), m_fails(0), m_att(0)
{
p_csplitter = new Cluster_Splitter();
}
Cluster_Part::~Cluster_Part()
{
if (p_csplitter) delete p_csplitter;
}
bool Cluster_Part::TestDecay(Cluster * const cluster)
{
m_att++;
if (p_csplitter && !(*p_csplitter)(cluster)) {
m_fails++;
msg_Out()<<"...... "<<METHOD<<" fails for\n"<<(*cluster);
return false;
}
return true;
}
#ifndef AHADIC_Decays_Cluster_Part_H
#define AHADIC_Decays_Cluster_Part_H
#include "AHADIC++/Tools/Cluster.H"
#include "AHADIC++/Tools/Cluster_Splitter.H"
#include "ATOOLS/Math/Histogram.H"
#include "AHADIC++/Tools/Hadronisation_Parameters.H"
#include "ATOOLS/Org/Return_Value.H"
#include "ATOOLS/Math/Histogram.H"
#include <map>
namespace AHADIC {
class Cluster_Part {
protected:
bool m_ana;
Cluster_Splitter * p_csplitter;
FlavDecayMap m_options;
long int m_fails, m_att;
std::map<std::string,ATOOLS::Histogram *> m_histograms;
public:
Cluster_Part(bool=false);
~Cluster_Part();
bool TestDecay(Cluster * const);
};
}
#endif
#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;
Cluster_Splitter::Cluster_Splitter(std::list<Cluster *> * cluster_list,
Soft_Cluster_Handler * softclusters) :
Splitter_Base(cluster_list,softclusters),
m_output(false)
{}
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");
}
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;
return CheckIfAllowed();
}
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);
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);
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::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;
}
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 =
new Proto_Particle((i==0?m_newflav1:m_newflav2),newmom12);
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);
}
}*/
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));
}
#ifndef AHADIC_Decays_Cluster_Splitter_H
#define AHADIC_Decays_Cluster_Splitter_H
#include "AHADIC++/Tools/Splitter_Base.H"
namespace AHADIC {
class Cluster_Splitter : public Splitter_Base {
private:
double m_alpha[3], m_beta[3], m_gamma[3], m_a, m_b, m_c;
size_t m_flcnt;
double m_R12, m_R21;
double m_m12min, m_m22min, m_masses2;
bool m_output;
bool MakeLongitudinalMomenta();
double DeltaM2(const double & MC2,const double & mt2,
const ATOOLS::Flavour & flav);
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);
double WeightFunction(const double & z,
const double & zmin,const double & zax);
void SetOutput(const bool & out) { m_output = out; }
};
}
#endif
......@@ -4,15 +4,12 @@ SVNTAG = AHADIC++/Decays
include ../../svn.make
libAhadicDecays_la_SOURCES = \
SVN_Info.C \
Cluster_Decay_Handler.C \
Cluster_Decay_Analysis.C \
Cluster_Part.C
SVN_Info.C \
Cluster_Decayer.C \
Cluster_Splitter.C
localincdir = $(pkgincludedir)/AHADIC++/Decays
localinc_HEADERS = \
Cluster_Decay_Handler.H \
Cluster_Decay_Analysis.H \
Cluster_Part.H
Cluster_Decayer.H \
Cluster_Splitter.H
#include "ATOOLS/Org/SVN_Info.H"
static ATOOLS::SVN_Info initializer
("AHADIC++/Decays","branches/rel-2-2-2_newhad","31037:31391M","a996e1db8b65efc861a3d529d779a078");
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
#ifndef AHADIC_Formations_Beam_Particles_Shifter_H
#define AHADIC_Formations_Beam_Particles_Shifter_H
#include "AHADIC++/Tools/Singlet_Tools.H"
#include "AHADIC++/Tools/Constituents.H"
#include <list>
namespace AHADIC {
class Beam_Particles_Shifter {
private:
std::list<Singlet *> * p_singlets;
std::list<Proto_Particle *> m_beamparts;
Constituents * p_constituents;
void ExtractBeamParticles();
bool ShiftBeamParticles();
public:
Beam_Particles_Shifter(std::list<Singlet *> * singlets);
~Beam_Particles_Shifter();
void Init();
void Reset();
bool operator()();
};
}
#endif
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
#include "ATOOLS/Org/SVN_Info.H"
static ATOOLS::SVN_Info initializer
("AHADIC++/Formation","branches/rel-2-2-2_newhad","31037:31391M","ab96bb66a9ab2455d9e33233eeee4ab7");
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
#include "ATOOLS/Org/SVN_Info.H"
static ATOOLS::SVN_Info initializer
("AHADIC++/Main","branches/rel-2-2-2_newhad","31037:31391M","6aec0d693c797d36f254e1137a42f6f9");
MAKE = make
SUBDIRS = Tools Formation Decays Main
SUBDIRS = Tools Formation Decays Main
DIST_SUBDIRS = Tools Formation Decays Main
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.