Commit d62c6637 authored by Frank Krauss's avatar Frank Krauss

another bug ... numerical precision and some physics in the transition to h+gamma

parent 4c758efb
...@@ -26,12 +26,12 @@ int Cluster_Decay_Handler::DecayClusters(Blob * blob) ...@@ -26,12 +26,12 @@ int Cluster_Decay_Handler::DecayClusters(Blob * blob)
{ {
Cluster * cluster; Cluster * cluster;
Cluster_Iterator cit(p_clulist->begin()); Cluster_Iterator cit(p_clulist->begin());
//msg_Out()<<METHOD<<" with "<<p_clulist->size()<<" clusters.\n"; //msg_Out()<<":::::: "<<METHOD<<" with "<<p_clulist->size()<<" clusters.\n";
while (!p_clulist->empty()) { while (!p_clulist->empty()) {
cluster = p_clulist->front(); cluster = p_clulist->front();
if (!cluster->Active()) return -1; if (!cluster->Active()) return -1;
if (p_clus->TestDecay(cluster)) { if (p_clus->TestDecay(cluster)) {
//msg_Out()<<METHOD<<": decay ok for\n"<<(*cluster); //msg_Out()<<":::::: "<<METHOD<<": decay ok for\n"<<(*cluster);
Cluster_List * clist(cluster->GetClusters()); Cluster_List * clist(cluster->GetClusters());
if (!p_softclusters->TreatClusterList(clist,blob)) { if (!p_softclusters->TreatClusterList(clist,blob)) {
msg_Error()<<"Error in "<<METHOD<<" : \n" msg_Error()<<"Error in "<<METHOD<<" : \n"
...@@ -44,15 +44,17 @@ int Cluster_Decay_Handler::DecayClusters(Blob * blob) ...@@ -44,15 +44,17 @@ int Cluster_Decay_Handler::DecayClusters(Blob * blob)
p_clulist->push_back(clist->front()); p_clulist->push_back(clist->front());
clist->pop_front(); clist->pop_front();
} }
//msg_Out()<<":::::: "<<p_clulist->size()<<" clusters now.\n";
} }
else { else {
//if (cluster==NULL) msg_Out()<<"Aha.\n"; //msg_Out()<<":::::: "<<METHOD<<":\n"
//msg_Out()<<"Enter soft cluster treatment with \n"<<(*cluster); // <<" Enter soft cluster treatment for undecayed cluster.\n"
// <<(*cluster);
Cluster_List clist; Cluster_List clist;
clist.push_back(cluster); clist.push_back(cluster);
if (!p_softclusters->TreatClusterList(&clist,blob)) { if (!p_softclusters->TreatClusterList(&clist,blob)) {
msg_Error()<<"Error in "<<METHOD<<".\n"; msg_Error()<<"Error in "<<METHOD<<".\n";
exit(1); return -1;
} }
} }
delete (p_clulist->front()->GetTrip()); delete (p_clulist->front()->GetTrip());
......
...@@ -23,7 +23,7 @@ bool Cluster_Part::TestDecay(Cluster * const cluster) ...@@ -23,7 +23,7 @@ bool Cluster_Part::TestDecay(Cluster * const cluster)
m_att++; m_att++;
if (p_csplitter && !(*p_csplitter)(cluster)) { if (p_csplitter && !(*p_csplitter)(cluster)) {
m_fails++; m_fails++;
//msg_Out()<<METHOD<<" fails for\n"<<(*cluster)<<"\n"; msg_Out()<<"...... "<<METHOD<<" fails for\n"<<(*cluster);
return false; return false;
} }
return true; return true;
......
...@@ -91,18 +91,18 @@ int Cluster_Formation_Handler::FormClusters(Blob * blob) { ...@@ -91,18 +91,18 @@ int Cluster_Formation_Handler::FormClusters(Blob * blob) {
if (!ClustersToHadrons(blob)) { if (!ClustersToHadrons(blob)) {
Reset(); return -1; Reset(); return -1;
} }
/*
Vec4D blobmom(0.,0.,0.,0.); // Vec4D blobmom(0.,0.,0.,0.);
for (size_t i(0);i<blob->NOutP();i++) // for (size_t i(0);i<blob->NOutP();i++)
blobmom+=blob->OutParticle(i)->Momentum(); // blobmom+=blob->OutParticle(i)->Momentum();
msg_Out()<<"____________________________________________________\n" // msg_Out()<<"____________________________________________________\n"
<<"____________________________________________________\n" // <<"____________________________________________________\n"
<<"Cluster list after all merging etc.:\n"<<(*p_clulist) // <<"Cluster list after all merging etc.:\n"<<(*p_clulist)
<<"____________________________________________________\n" // <<"____________________________________________________\n"
<<"Blob momentum: "<<blobmom<<";\n"<<(*blob)<<"\n" // <<"Blob momentum: "<<blobmom<<";\n"<<(*blob)<<"\n"
<<"____________________________________________________\n" // <<"____________________________________________________\n"
<<"____________________________________________________\n"; // <<"____________________________________________________\n";
*/
return 1; return 1;
} }
......
...@@ -55,10 +55,10 @@ Return_Value::code Ahadic::Hadronize(Blob_List * blobs) ...@@ -55,10 +55,10 @@ Return_Value::code Ahadic::Hadronize(Blob_List * blobs)
if ((*blit)->Has(blob_status::needs_hadronization) && if ((*blit)->Has(blob_status::needs_hadronization) &&
(*blit)->Type()==btp::Fragmentation) { (*blit)->Type()==btp::Fragmentation) {
blob = (*blit); blob = (*blit);
//msg_Out()<<"====================================================\n" // msg_Out()<<"====================================================\n"
// <<"====================================================\n" // <<"====================================================\n"
// <<"====================================================\n" // <<"====================================================\n"
// <<(*blob)<<"\n"; // <<(*blob)<<"\n";
moveon = false; moveon = false;
Reset(); Reset();
for (short int i=0;i<m_maxtrials;i++) { for (short int i=0;i<m_maxtrials;i++) {
...@@ -175,7 +175,6 @@ Return_Value::code Ahadic::Hadronize(Blob * blob,int retry) { ...@@ -175,7 +175,6 @@ Return_Value::code Ahadic::Hadronize(Blob * blob,int retry) {
break; break;
} }
//msg_Out()<<METHOD<<" formed clusters:\n"<<(*blob)<<".\n";
switch (p_cdechandler->DecayClusters(blob)) { switch (p_cdechandler->DecayClusters(blob)) {
case -1 : case -1 :
msg_Tracking()<<"ERROR in "<<METHOD<<" :"<<std::endl msg_Tracking()<<"ERROR in "<<METHOD<<" :"<<std::endl
......
...@@ -317,15 +317,16 @@ void Cluster::Delete() { ...@@ -317,15 +317,16 @@ void Cluster::Delete() {
bool Cluster::EnsureMomentum() { bool Cluster::EnsureMomentum() {
Vec4D check(Momentum()); Vec4D check(Momentum());
double E(dabs(check[0]));
for (Cluster_Iterator cit(m_clusters.begin());cit!=m_clusters.end();cit++) for (Cluster_Iterator cit(m_clusters.begin());cit!=m_clusters.end();cit++)
check -= (*cit)->Momentum(); check -= (*cit)->Momentum();
if (dabs(check.Abs2())>1.e-10 || dabs(check[0])>1.e-10 || if (dabs(check.Abs2()/(E*E))>1.e-6 || dabs(check[0]/E)>1.e-6 ||
dabs(check.PSpat())>1.e-10) { dabs(check.PSpat()/E)>1.e-6) {
// msg_Out()<<"*** "<<METHOD<<" yields "<<check<<" ("<<check.Abs2()<<") " msg_Out()<<" --> "<<METHOD<<": simple construction not ok:\n"
// <<"for cluster "<<Number()<<"\n"<<(*this) <<" --> "<<check<<" ("<<check.Abs2()<<").\n";//<<(*this);
// <<"********************************************************\n";
return false; return false;
} }
//msg_Out()<<" --> "<<METHOD<<": four momentum ok now.\n";
return true; return true;
} }
......
...@@ -32,8 +32,9 @@ bool Cluster_Splitter::operator()(Cluster * cluster) { ...@@ -32,8 +32,9 @@ bool Cluster_Splitter::operator()(Cluster * cluster) {
if (ConstructLightC() && ConstructSystem(cluster)) { if (ConstructLightC() && ConstructSystem(cluster)) {
if (m_ana) Analysis(); if (m_ana) Analysis();
Reset(); Reset();
if (!cluster->EnsureMomentum()) return EnforceMomentum(cluster); if (!cluster->EnsureMomentum() && !EnforceMomentum(cluster)) {
//msg_Out()<<(*cluster)<<"\n"; return false;
}
return true; return true;
} }
UndoTrafos(); UndoTrafos();
...@@ -67,7 +68,6 @@ bool Cluster_Splitter::ConstructSystem(Cluster * cluster) { ...@@ -67,7 +68,6 @@ bool Cluster_Splitter::ConstructSystem(Cluster * cluster) {
m_cms = p_split->m_mom + p_spect->m_mom; m_cms = p_split->m_mom + p_spect->m_mom;
pair<double,double> exponents(FixExponents()); pair<double,double> exponents(FixExponents());
bool hit; bool hit;
//msg_Out()<<METHOD<<": pt^2_max = "<<m_pt2max<<".\n";
double pt2max(m_pt2max); double pt2max(m_pt2max);
if (m_leadsplit) pt2max *= m_pt2max/Max(m_pt2max,m_LC.m_msplit2); if (m_leadsplit) pt2max *= m_pt2max/Max(m_pt2max,m_LC.m_msplit2);
if (m_leadspect) pt2max *= m_pt2max/Max(m_pt2max,m_LC.m_mspect2); if (m_leadspect) pt2max *= m_pt2max/Max(m_pt2max,m_LC.m_mspect2);
...@@ -88,7 +88,11 @@ bool Cluster_Splitter::ConstructSystem(Cluster * cluster) { ...@@ -88,7 +88,11 @@ bool Cluster_Splitter::ConstructSystem(Cluster * cluster) {
} }
if (hit) break; if (hit) break;
} }
if (!hit) return false; if (!hit) {
//msg_Out()<<" --> "<<METHOD<<" failed for m_min = "
// <<sqrt(m_mmin2)<<".\n";
return false;
}
MakeKinematics(); MakeKinematics();
MakeClusters(cluster); MakeClusters(cluster);
return true; return true;
...@@ -120,7 +124,7 @@ ConstructKinematics(const double & etax,const double & etay) { ...@@ -120,7 +124,7 @@ ConstructKinematics(const double & etax,const double & etay) {
weight = 0.; weight = 0.;
} }
else { else {
z = SelectZ(m_mmin2/sqq,m_leadspect || m_leadsplit); z = SelectZ(4.*m_mmin2/sqq,m_leadspect || m_leadsplit);
weight = exp(-(sqq-sqqmin)/(4.*m_pt02)); weight = exp(-(sqq-sqqmin)/(4.*m_pt02));
//weight = 1.; //weight = 1.;
} }
...@@ -140,8 +144,11 @@ ConstructKinematics(const double & etax,const double & etay) { ...@@ -140,8 +144,11 @@ ConstructKinematics(const double & etax,const double & etay) {
bool Cluster_Splitter::AcceptSystem(const double & pt2max) { bool Cluster_Splitter::AcceptSystem(const double & pt2max) {
PoppedPair * pop(m_popped.back()); PoppedPair * pop(m_popped.back());
pop->m_kt2 = pop->m_z*(1.-pop->m_z)*pop->m_sqq-pop->m_mpop2; pop->m_kt2 = pop->m_z*(1.-pop->m_z)*pop->m_sqq-pop->m_mpop2;
if (pop->m_kt2 < 0.) return false; if (pop->m_kt2 < 0. || pop->m_kt2 > pt2max) {
if (pop->m_kt2 > pt2max) return false; //msg_Out()<<" --> "<<METHOD<<" kt2 veto: "<<pop->m_kt2
// <<" from s = "<<pop->m_sqq<<", z = "<<pop->m_z<<".\n";
return false;
}
return (((*p_as)(pop->m_sqq,false)*(*p_as)(pop->m_kt2,false))/ return (((*p_as)(pop->m_sqq,false)*(*p_as)(pop->m_kt2,false))/
sqr(p_as->MaxValue())) > ran->Get(); sqr(p_as->MaxValue())) > ran->Get();
} }
...@@ -184,8 +191,6 @@ MakePairKinematics(PoppedPair * pop,Vec4D & test,Vec4D & test1) { ...@@ -184,8 +191,6 @@ MakePairKinematics(PoppedPair * pop,Vec4D & test,Vec4D & test1) {
double kt(sqrt(pop->m_kt2)); double kt(sqrt(pop->m_kt2));
double phi(2.*M_PI*ran->Get()); double phi(2.*M_PI*ran->Get());
Vec4D kperp(0.,kt*cos(phi),kt*sin(phi),0.); Vec4D kperp(0.,kt*cos(phi),kt*sin(phi),0.);
//msg_Out()<<METHOD<<": sqq = "<<pop->m_sqq<<", z = "<<pop->m_z<<", "
// <<"kt = "<<kt<<", phi = "<<phi<<"\n";
test += pop->m_outmom[0] = test += pop->m_outmom[0] =
pop->m_x*pop->m_z*m_LC.m_pA + pop->m_x*pop->m_z*m_LC.m_pA +
pop->m_y*(1.-pop->m_z)*m_LC.m_pB + kperp; pop->m_y*(1.-pop->m_z)*m_LC.m_pB + kperp;
...@@ -194,7 +199,6 @@ MakePairKinematics(PoppedPair * pop,Vec4D & test,Vec4D & test1) { ...@@ -194,7 +199,6 @@ MakePairKinematics(PoppedPair * pop,Vec4D & test,Vec4D & test1) {
pop->m_y*pop->m_z*m_LC.m_pB - kperp; pop->m_y*pop->m_z*m_LC.m_pB - kperp;
m_sumx += pop->m_x; m_sumx += pop->m_x;
m_sumy += pop->m_y; m_sumy += pop->m_y;
//msg_Out()<<pop->m_outmom[0]<<" + "<<pop->m_outmom[1]<<"\n";
for (size_t i(0);i<2;i++) { for (size_t i(0);i<2;i++) {
m_rotat.RotateBack(pop->m_outmom[i]); m_rotat.RotateBack(pop->m_outmom[i]);
m_boost.BoostBack(pop->m_outmom[i]); m_boost.BoostBack(pop->m_outmom[i]);
...@@ -368,19 +372,34 @@ void Cluster_Splitter::SelectPartners() { ...@@ -368,19 +372,34 @@ void Cluster_Splitter::SelectPartners() {
bool Cluster_Splitter::PoppedMassPossible(const double & m2) { bool Cluster_Splitter::PoppedMassPossible(const double & m2) {
PoppedPair * pop(m_popped.back()); PoppedPair * pop(m_popped.back());
pop->m_kt2 = pop->m_z*(1.-pop->m_z)*pop->m_sqq-m2; pop->m_kt2 = pop->m_z*(1.-pop->m_z)*pop->m_sqq-m2;
if (pop->m_kt2<0.) return false; if (pop->m_kt2<0.) {
//if (m2<=1.001*m_mmin2)
//msg_Out()<<" --> "<<METHOD<<"(sqq = "<<pop->m_sqq<<", "
// <<"z = "<<pop->m_z<<" --> "<<pop->m_kt2<<").\n";
return false;
}
double sumx = m_sumx+pop->m_x; double sumx = m_sumx+pop->m_x;
double sumy = m_sumy+pop->m_y; double sumy = m_sumy+pop->m_y;
double snew((1.-sumx)*(1.-sumy)*m_LC.m_smandel); double snew((1.-sumx)*(1.-sumy)*m_LC.m_smandel);
if (snew<sqr(m_LC.m_mspect+m_LC.m_msplit)) return false; if (snew<sqr(m_LC.m_mspect+m_LC.m_msplit) ||
if (sumx>1. || sumy>1.) return false; sumx>1. || sumy>1.||
if (!AlphaBeta(snew,m_alphanew,m_betanew)) return false; !AlphaBeta(snew,m_alphanew,m_betanew)) {
//msg_Out()<<" --> "<<METHOD<<" testing LC's failed.\n";
return false;
}
Vec4D tsplit((1.-sumx)*(1.-m_alphanew)*m_LC.m_pA+ Vec4D tsplit((1.-sumx)*(1.-m_alphanew)*m_LC.m_pA+
(1.-sumy)*m_betanew*m_LC.m_pB); (1.-sumy)*m_betanew*m_LC.m_pB);
Vec4D tspect((1.-sumx)*m_alphanew*m_LC.m_pA+ Vec4D tspect((1.-sumx)*m_alphanew*m_LC.m_pA+
(1.-sumy)*(1.-m_betanew)*m_LC.m_pB); (1.-sumy)*(1.-m_betanew)*m_LC.m_pB);
if (dabs(tsplit.Abs2()/m_LC.m_msplit2-1.)>1.e-8) return false; if (dabs(tsplit.Abs2()/m_LC.m_msplit2-1.)>1.e-6 ||
if (dabs(tspect.Abs2()/m_LC.m_mspect2-1.)>1.e-8) return false; dabs(tspect.Abs2()/m_LC.m_mspect2-1.)>1.e-6) {
//msg_Out()<<" --> "<<METHOD<<" failed: "
// <<dabs(tsplit.Abs2())<<"/"<<m_LC.m_msplit2
// <<" -> "<<dabs(tsplit.Abs2()/m_LC.m_msplit2-1.)<<" and "
// <<dabs(tspect.Abs2())<<"/"<<m_LC.m_mspect2
// <<" -> "<<dabs(tspect.Abs2()/m_LC.m_mspect2-1.)<<".\n";
return false;
}
return true; return true;
} }
...@@ -399,6 +418,13 @@ bool Cluster_Splitter::EnforceMomentum(Cluster * cluster) { ...@@ -399,6 +418,13 @@ bool Cluster_Splitter::EnforceMomentum(Cluster * cluster) {
cit!=cluster->GetClusters()->end();cit++) { cit!=cluster->GetClusters()->end();cit++) {
summom += (*cit)->Momentum(); summom += (*cit)->Momentum();
} }
//msg_Out()<<" --> "<<METHOD<<" for \n"
// <<summom<<" vs. "<<cluster->Momentum()<<" --> "
// <<(summom-cluster->Momentum())<<".\n";
//for (Cluster_Iterator cit(cluster->GetClusters()->begin());
// cit!=cluster->GetClusters()->end();cit++) {
//msg_Out()<<(**cit);
//}
Poincare rest(summom); Poincare rest(summom);
Poincare back(cluster->Momentum()); Poincare back(cluster->Momentum());
for (Cluster_Iterator cit(cluster->GetClusters()->begin()); for (Cluster_Iterator cit(cluster->GetClusters()->begin());
...@@ -406,5 +432,6 @@ bool Cluster_Splitter::EnforceMomentum(Cluster * cluster) { ...@@ -406,5 +432,6 @@ bool Cluster_Splitter::EnforceMomentum(Cluster * cluster) {
(*cit)->Boost(rest); (*cit)->Boost(rest);
(*cit)->BoostBack(back); (*cit)->BoostBack(back);
} }
return cluster->EnsureMomentum(); if (!cluster->EnsureMomentum()) exit(1);
return true;
} }
...@@ -18,6 +18,8 @@ namespace AHADIC { ...@@ -18,6 +18,8 @@ namespace AHADIC {
size_t m_popspliti, m_popspecti; size_t m_popspliti, m_popspecti;
Proto_Particle * p_trip, * p_anti; Proto_Particle * p_trip, * p_anti;
long int m_fourmom_err, m_kt2_err, m_LC_err, m_flav_err;
size_t SelectNumberOfPairs(const size_t & nmax); size_t SelectNumberOfPairs(const size_t & nmax);
bool SelectSplitter(Proto_Particle * part1,Proto_Particle * part2); bool SelectSplitter(Proto_Particle * part1,Proto_Particle * part2);
std::pair<double,double> FixExponents(); std::pair<double,double> FixExponents();
......
...@@ -25,7 +25,7 @@ Soft_Cluster_Handler::Soft_Cluster_Handler(bool ana) : ...@@ -25,7 +25,7 @@ Soft_Cluster_Handler::Soft_Cluster_Handler(bool ana) :
m_pt02(hadpars->Get(std::string("pt02"))), m_pt02(hadpars->Get(std::string("pt02"))),
m_transitions(0), m_dtransitions(0), m_decays(0), m_transitions(0), m_dtransitions(0), m_decays(0),
m_forceddecays(0), m_lists(0), m_update(0), m_forceddecays(0), m_lists(0), m_update(0),
m_ana(ana) m_ana(ana), m_out(false)
{ {
if (m_ana) { if (m_ana) {
m_histograms[string("PT_HH")] = new Histogram(0,0.,10.,100); m_histograms[string("PT_HH")] = new Histogram(0,0.,10.,100);
...@@ -53,17 +53,19 @@ Soft_Cluster_Handler::~Soft_Cluster_Handler() ...@@ -53,17 +53,19 @@ Soft_Cluster_Handler::~Soft_Cluster_Handler()
bool Soft_Cluster_Handler::TreatClusterList(Cluster_List * clin, Blob * blob) bool Soft_Cluster_Handler::TreatClusterList(Cluster_List * clin, Blob * blob)
{ {
//msg_Out()<<"#######################################################" // if (m_out)
// <<"#######################################################\n" // msg_Out()<<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++"
// <<"###### "<<METHOD<<"("<<clin->size()<<" clusters):\n"; // <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
// <<"++++++ "<<METHOD<<"("<<clin->size()<<" clusters):\n";
// checks for transitions to hadrons and attaches them, if neccessary // checks for transitions to hadrons and attaches them, if neccessary
if (!CheckListForTreatment(clin)) { if (!CheckListForTreatment(clin) && m_out) {
//msg_Out()<<"###### No hadrons produced. Just continue.\n" // msg_Out()<<"++++++ No hadrons produced. Just continue.\n"
// <<"#######################################################" // <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++"
// <<"#######################################################\n"; // <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++\n";
return true; return false;
} }
//msg_Out()<<"###### Hadrons produced, attach to blob.\n"; // if (m_out)
// msg_Out()<<"++++++ Hadrons produced, will attach to blob.\n";
return AttachHadronsToBlob(clin,blob); return AttachHadronsToBlob(clin,blob);
} }
...@@ -82,16 +84,21 @@ bool Soft_Cluster_Handler::CheckListForTreatment(Cluster_List * clin) { ...@@ -82,16 +84,21 @@ bool Soft_Cluster_Handler::CheckListForTreatment(Cluster_List * clin) {
} }
int Soft_Cluster_Handler::CheckCluster(Cluster * cluster) { int Soft_Cluster_Handler::CheckCluster(Cluster * cluster) {
//msg_Out()<<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++" // if (m_out)
// <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++\n" // msg_Out()<<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++"
// <<"++++++ "<<METHOD<<" for:\n"<<(*cluster); // <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
// <<"++++++ "<<METHOD<<" for:\n"<<(*cluster);
cluster->clear(); cluster->clear();
Flavour haddec1(Flavour(kf_none)), haddec2(Flavour(kf_none)); Flavour haddec1(Flavour(kf_none)), haddec2(Flavour(kf_none));
Flavour hadtrans(Flavour(kf_none)); Flavour hadtrans(Flavour(kf_none));
double decayweight(DecayWeight(cluster,haddec1,haddec2)); double decayweight(DecayWeight(cluster,haddec1,haddec2));
double transweight(TransformWeight(cluster,hadtrans)); double transweight(TransformWeight(cluster,hadtrans));
//msg_Out()<<"++++++ "<<METHOD<<"(dec = "<<decayweight<<", " // if (m_out)
// <<"trans = "<<transweight<<").\n"; // msg_Out()<<"++++++ "<<METHOD<<"["<<cluster->Mass()<<" "
// <<"("<<cluster->GetTrip()->m_flav<<" + "
// <<cluster->GetAnti()->m_flav<<"] --> "
// <<"(dec = "<<decayweight<<", "
// <<"trans = "<<transweight<<").\n";
if (decayweight>0.) { if (decayweight>0.) {
if (transweight>0.) { if (transweight>0.) {
double totweight(decayweight+transweight); double totweight(decayweight+transweight);
...@@ -100,20 +107,22 @@ int Soft_Cluster_Handler::CheckCluster(Cluster * cluster) { ...@@ -100,20 +107,22 @@ int Soft_Cluster_Handler::CheckCluster(Cluster * cluster) {
cluster->push_back(haddec1); cluster->push_back(haddec1);
cluster->push_back(haddec2); cluster->push_back(haddec2);
m_decays += 1; m_decays += 1;
//msg_Out()<<"++++++ decays to "<<haddec1<<" + "<<haddec2<<".\n" // if (m_out)
// <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++" // msg_Out()<<"++++++ decays to "<<haddec1<<" + "<<haddec2<<".\n"
// <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++" // <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++"
// <<"\n\n"; // <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++"
// <<"\n\n";
return 2; return 2;
} }
else { else {
// competition between decay and transition - transition wins // competition between decay and transition - transition wins
cluster->push_back(hadtrans); cluster->push_back(hadtrans);
cluster->push_back(Flavour(kf_photon)); cluster->push_back(Flavour(kf_photon));
//msg_Out()<<"++++++ decays to "<<hadtrans<<" + photon.\n" // if (m_out)
// <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++" // msg_Out()<<"++++++ decays to "<<hadtrans<<" + photon.\n"
// <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++" // <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++"
// <<"\n\n"; // <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++"
// <<"\n\n";
m_transitions += 1; m_transitions += 1;
return 2; return 2;
} }
...@@ -122,9 +131,11 @@ int Soft_Cluster_Handler::CheckCluster(Cluster * cluster) { ...@@ -122,9 +131,11 @@ int Soft_Cluster_Handler::CheckCluster(Cluster * cluster) {
// regular decay, and no simple transition open - so no competition // regular decay, and no simple transition open - so no competition
cluster->push_back(haddec1); cluster->push_back(haddec1);
cluster->push_back(haddec2); cluster->push_back(haddec2);
//msg_Out()<<"++++++ decays to "<<haddec1<<" + "<<haddec2<<".\n" // if (m_out)
// <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++" // msg_Out()<<"++++++ decays to "<<haddec1<<" + "<<haddec2<<".\n"
// <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++\n\n"; // <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++"
// <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++"
// <<"\n\n";
m_decays += 1; m_decays += 1;
return 2; return 2;
} }
...@@ -134,16 +145,18 @@ int Soft_Cluster_Handler::CheckCluster(Cluster * cluster) { ...@@ -134,16 +145,18 @@ int Soft_Cluster_Handler::CheckCluster(Cluster * cluster) {
if (transweight<=0.) transweight = TransformWeight(cluster,hadtrans,true); if (transweight<=0.) transweight = TransformWeight(cluster,hadtrans,true);
cluster->push_back(hadtrans); cluster->push_back(hadtrans);
cluster->push_back(Flavour(kf_photon)); cluster->push_back(Flavour(kf_photon));
//msg_Out()<<"++++++ decays to "<<hadtrans<<" + photon.\n" // if (m_out)
// <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++" // msg_Out()<<"++++++ decays to "<<hadtrans<<" + photon.\n"
// <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++\n\n"; // <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++"
// <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++\n\n";
m_transitions += 1; m_transitions += 1;
return 2; return 2;
} }
// no decay and no transition: both weight equal 0. // no decay and no transition: both weight equal 0.
//msg_Out()<<"++++++ no decay.\n" // if (m_out)
// <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++" // msg_Out()<<"++++++ no decay.\n"
// <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++\n\n"; // <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++"
// <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++\n\n";
cluster->clear(); cluster->clear();
return 0; return 0;
} }
...@@ -172,9 +185,11 @@ AttachHadronsToBlob(Cluster_List * clin,Blob * blob) ...@@ -172,9 +185,11 @@ AttachHadronsToBlob(Cluster_List * clin,Blob * blob)
break; break;
} }
} }
//msg_Out()<<"###### "<<METHOD<<" was successful.\n" // if (m_out)
// <<"#######################################################" // msg_Out()<<"++++++ "<<METHOD<<" was successful:"
// <<"#######################################################\n"; // <<blob->CheckMomentumConservation()<<"\n"
// <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++"
// <<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++\n";
return true; return true;
} }
...@@ -230,16 +245,16 @@ TransformWeight(Cluster * cluster,Flavour & hadron,const bool & enforce) ...@@ -230,16 +245,16 @@ TransformWeight(Cluster * cluster,Flavour & hadron,const bool & enforce)
break; break;
} }
} }
return wt/(16.*M_PI*MC); return totweight/(16.*M_PI*MC)/137.;
} }
double Soft_Cluster_Handler:: double Soft_Cluster_Handler::
TransformKin(const double MC,const Flavour & flav,const bool & enforce) { TransformKin(const double MC,const Flavour & flav,const bool & enforce) {
double mass2(sqr(flav.HadMass())); double mass2(sqr(flav.HadMass()));
double width2(sqr(Max(flav.Width(),1.e-6))); double width2(sqr(Max(flav.Width(),1.e-8)));
return return
pow(sqr(mass2)/(sqr(MC*MC-mass2) + mass2*width2),m_kappa) * pow(sqr(mass2)/(sqr(MC*MC-mass2) + mass2*width2),m_kappa) *
pow(mass2*width2/(sqr(MC*MC-mass2) + mass2*width2),m_lambda); pow(mass2*width2/(sqr(MC*MC-mass2) + mass2*width2),1.+m_lambda);
} }
...@@ -326,7 +341,7 @@ DecayWeight(Cluster * cluster,Flavour & had1,Flavour & had2) ...@@ -326,7 +341,7 @@ DecayWeight(Cluster * cluster,Flavour & had1,Flavour & had2)
} }
} }
} }
return wt/(16.*M_PI*MC*MC*MC); return totweight/(16.*M_PI*MC*MC*MC);
} }
...@@ -408,13 +423,13 @@ void Soft_Cluster_Handler::FixHHDecay(Cluster * cluster,Blob * blob, ...@@ -408,13 +423,13 @@ void Soft_Cluster_Handler::FixHHDecay(Cluster * cluster,Blob * blob,
blob->AddToOutParticles(left); blob->AddToOutParticles(left);
blob->AddToOutParticles(right); blob->AddToOutParticles(right);
} }
// if (cluster->GetTrip()->m_info=='B' || cluster->GetAnti()->m_info=='B') { //if (cluster->GetTrip()->m_info=='B' || cluster->GetAnti()->m_info=='B') {
// msg_Out()<<"==========================================================\n" // msg_Out()<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
// <<"Cluster decay (pt = "<<pt<<") for cluster \n"<<(*cluster)<<"\n" // <<"Cluster decay (pt + "<<pt<<") for cluster \n"<<(*cluster)<<"\n"
// <<"==> "<<left->Momentum()<<" + "<<right->Momentum()<<" for " // <<"++> "<<left->Momentum()<<" + "<<right->Momentum()<<" for "
// <<left->Flav()<<" + "<<right->Flav()<<"\n" // <<left->Flav()<<" + "<<right->Flav()<<"\n"
// <<"==========================================================\n"; // <<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n";
// } //}
if (m_ana) { if (m_ana) {
Histogram* histo((m_histograms.find(std::string("PT_HH")))->second); Histogram* histo((m_histograms.find(std::string("PT_HH")))->second);
histo->Insert(pt); histo->Insert(pt);
......
...@@ -25,7 +25,7 @@ namespace AHADIC { ...@@ -25,7 +25,7 @@ namespace AHADIC {
long int m_transitions, m_dtransitions; long int m_transitions, m_dtransitions;
long int m_decays, m_forceddecays, m_lists, m_update; long int m_decays, m_forceddecays, m_lists, m_update;
bool m_ana; bool m_ana, m_out;
std::map<std::string,ATOOLS::Histogram *> m_histograms; std::map<std::string,ATOOLS::Histogram *> m_histograms;
bool CheckListForTreatment(Cluster_List *); bool CheckListForTreatment(Cluster_List *);
......
...@@ -173,29 +173,40 @@ double Splitter_Base::SelectZ(const double & delta,const bool & lead) { ...@@ -173,29 +173,40 @@ double Splitter_Base::SelectZ(const double & delta,const bool & lead) {
bool Splitter_Base::SelectFlavour(const double & sqq,const bool & vetodi) { bool Splitter_Base::SelectFlavour(const double & sqq,const bool & vetodi) {
Flavour flav(kf_none); Flavour flav(kf_none);
double sumwt(CalculateSumWT(sqrt(sqq/4.),vetodi)); double mmax(sqrt(sqq/4.)), m2, sumwt;
double mmax(1.e12), m2;
long int calls(0); long int calls(0);
while (calls<100) { while (calls<100) {
flav = Flavour(kf_none); sumwt = CalculateSumWT(mmax,vetodi);
calls++; if (sumwt>1.e-12) {
sumwt *= ran->Get(); flav = Flavour(kf_none);
for (FDIter fdit=m_options.begin();fdit!=m_options.end();fdit++) { calls++;
if (fdit->second->popweight>0. && fdit->second->massmin<mmax && sumwt *= ran->Get();
!(vetodi && fdit->first.IsDiQuark())) for (FDIter fdit=m_options.begin();fdit!=m_options.end();fdit++) {
sumwt -= fdit->second->popweight; if (fdit->second->popweight>0. && fdit->second->massmin<mmax &&
if (sumwt<0.) { !(vetodi && fdit->first.IsDiQuark()))
flav = fdit->first; sumwt -= fdit->second->popweight;
m2 = sqr(fdit->second->massmin); if (sumwt<0.) {
break; flav = fdit->first;
m2 = sqr(fdit->second->massmin);
break;
}
} }
if (PoppedMassPossible(m2)) break;
mmax = sqrt(m2);
} }
if (PoppedMassPossible(m2)) break; else {
else mmax = sqrt(m2); for (FDIter fdit=m_options.begin();fdit!=m_options.end();fdit++) {
sumwt = CalculateSumWT(mmax); if (vetodi && fdit->first.IsDiQuark()) continue;
if (sumwt<=1.e-12) calls=100; }
calls = 100;
}
}
if (calls>=100) {
//msg_Out()<<" --> "<<METHOD<<" failed, too many calls.\n"
// <<" --> sumwt = "<<sumwt<<" for minimal mass = "
// <<sqrt(sqq/4.)<<".\n";
return false;
} }
if (calls>=100) return false;
m_popped.back()->m_flav = flav.IsDiQuark()?flav.Bar():flav; m_popped.back()->m_flav = flav.IsDiQuark()?flav.Bar():flav;
m_popped.back()->m_mpop2 = m_popped.back()->m_mpop2 =
sqr(hadpars->GetConstituents()->Mass(m_popped.back()->m_flav)); sqr(hadpars->GetConstituents()->Mass(m_popped.back()->m_flav));
......
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