Commit 8cdadc61 authored by Frank Krauss's avatar Frank Krauss

debugged the RescaleCluster bug

parent c0eafef5
......@@ -44,7 +44,7 @@ int Cluster_Decay_Handler::DecayClusters(Blob * blob)
clist.empty();
clist.push_back(cluster->GetLeft());
clist.push_back(cluster->GetRight());
msg_Tracking()<<"++++ From "<<cluster->Number()<<": "
msg_Tracking()<<"++++ From "<<cluster->Number()<<"("<<cluster->Mass()<<"): "
<<cluster->GetLeft()->Number()<<" ("<<cluster->GetLeft()->Mass()<<") "
<<cluster->GetRight()->Number()<<" ("<<cluster->GetRight()->Mass()<<"), "
<<"popped "<<cluster->GetLeft()->GetAnti()->m_flav<<"."<<std::endl;
......@@ -61,11 +61,12 @@ int Cluster_Decay_Handler::DecayClusters(Blob * blob)
}
}
else {
msg_Tracking()<<"+++ TestDecay did not work out - try to enforce decay."<<std::endl;
if (!p_softclusters->EnforcedDecay(cluster,blob,true)) {
msg_Tracking()<<" no decay."<<std::endl
<<"++++ TestDecay did not work out, will return -1."<<std::endl;
msg_Error()<<"+++ EnforcedDecay did not work out, will return -1."<<std::endl;
return -1;
}
msg_Tracking()<<"+++ EnforcedDecay did work out, will continue."<<std::endl;
}
delete (p_clulist->front()->GetTrip());
delete (p_clulist->front()->GetAnti());
......@@ -92,3 +93,4 @@ ATOOLS::Blob * Cluster_Decay_Handler::ClusterDecayBlob(Cluster * cluster,Cluster
}
......@@ -37,21 +37,14 @@ bool Cluster_Part::TestDecay(Cluster * const cluster)
msg_Tracking()<<":::::::::::::::::::::::::::::::::::::::::::::::"<<std::endl
<<"::: "<<METHOD<<" : Try "<<cluster->Number()<<" ("
<<cluster->GetTrip()->m_flav<<" "<<cluster->GetAnti()->m_flav<<", "
<<"m = "<<cluster->Mass()<<") --> "<<std::endl;
<<"mass = "<<cluster->Mass()<<") --> "<<std::endl;
if (!p_splitter->SplitCluster(cluster)) {
msg_Tracking()<<"::: Warning in "<<METHOD<<":"<<std::endl
<<"::: Could not split cluster ("<<cluster->Number()<<"): "
<<cluster->GetTrip()->m_flav<<"/"<<cluster->GetAnti()->m_flav<<", "
<<"mass = "<<cluster->Mass()<<","<<std::endl
<<"::: try to enforce splitting (not implemented)."<<std::endl;
bool worked(p_splitter->EnforceSplit(cluster));
if (worked) {
msg_Tracking()<<"::: "<<METHOD<<" : forced decay of "<<cluster->Momentum()<<" succeded."<<std::endl
<<(*cluster->GetLeft())<<(*cluster->GetRight())
<<":::::::::::::::::::::::::::::::::::::::::::::::"<<std::endl;
}
return worked;
<<"::: try to enforce cluster splitting into hadrons."<<std::endl;
return false;
}
if (m_ana) {
Vec4D lmom(cluster->GetLeft()->Momentum());
......
......@@ -216,11 +216,11 @@ bool Ahadic::SanityCheck(Blob * blob,double norm2) {
<<(*blob)<<endl;
return false;
}
msg_Debugging()<<"Passed "<<METHOD<<" with "
<<"protoparticles = "<<control::s_AHAprotoparticles
<<"/ parts = "<<control::s_AHAparticles<<" vs. "<<blob->NOutP()
<<" : "<<blob->CheckMomentumConservation()<<endl
<<(*blob)<<endl;
msg_Tracking()<<"Passed "<<METHOD<<" with "
<<"protoparticles = "<<control::s_AHAprotoparticles
<<"/ parts = "<<control::s_AHAparticles<<" vs. "<<blob->NOutP()
<<" : "<<blob->CheckMomentumConservation()<<endl;
msg_Debugging()<<(*blob)<<endl;
return true;
}
......
......@@ -199,17 +199,17 @@ void Cluster::RescaleMomentum(ATOOLS::Vec4D newmom)
Vec4D_Vector save(3+int(p_left!=NULL)+int(p_right!=NULL));
save[0] = m_momentum;
if (p_trip!=NULL) save[1] = p_trip->m_mom;
if (p_anti!=NULL) save[2] = p_anti->m_mom;
if (p_left!=NULL) save[3] = p_trip->m_mom;
if (p_right!=NULL) save[4] = p_anti->m_mom;
if (p_trip!=NULL) rest.Boost(p_trip->m_mom);
if (p_trip!=NULL) back.BoostBack(p_trip->m_mom);
if (p_anti!=NULL) rest.Boost(p_anti->m_mom);
if (p_anti!=NULL) back.BoostBack(p_anti->m_mom);
if (p_left!=NULL) p_left->Boost(rest);
if (p_left!=NULL) p_left->BoostBack(back);
if (p_trip!=NULL) save[1] = p_trip->m_mom;
if (p_anti!=NULL) save[2] = p_anti->m_mom;
if (p_left!=NULL) save[3] = p_trip->m_mom;
if (p_right!=NULL) save[4] = p_anti->m_mom;
if (p_trip!=NULL) rest.Boost(p_trip->m_mom);
if (p_trip!=NULL) back.BoostBack(p_trip->m_mom);
if (p_anti!=NULL) rest.Boost(p_anti->m_mom);
if (p_anti!=NULL) back.BoostBack(p_anti->m_mom);
if (p_left!=NULL) p_left->Boost(rest);
if (p_left!=NULL) p_left->BoostBack(back);
if (p_right!=NULL) p_right->Boost(rest);
if (p_right!=NULL) p_right->BoostBack(back);
m_momentum = newmom;
......@@ -221,19 +221,23 @@ void Cluster::RescaleMomentum(ATOOLS::Vec4D newmom)
Vec4D testmom = m_momentum-p_trip->m_mom-p_anti->m_mom;
if (dabs(testmom.Abs2()/save[0][0])>1.e-6 || testmom[0]/save[0][0]>1.e-6) {
msg_Error()<<"Error in "<<METHOD<<":"<<std::endl
<<" From "<<save[0]<<" ("<<sqrt(Max(0.,save[0].Abs2()))<<") to "<<m_momentum<<" with "<<std::endl
<<" From "<<save[0]<<" ("<<sqrt(Max(0.,save[0].Abs2()))<<") to "
<<m_momentum<<" with "<<std::endl
<<" "<<save[1]<<" ("<<save[1].Abs2()<<") + "
<<save[2]<<" ("<<save[2].Abs2()<<")"<<std::endl;
if (p_trip!=NULL) msg_Error()<<" Trip: "<<p_trip->m_mom<<" ("<<p_trip->m_mom.Abs2()<<")";
else msg_Error()<<"No triplet: "<<p_trip<<" ";
if (p_anti!=NULL) msg_Error()<<" Anti: "<<p_anti->m_mom<<" ("<<p_anti->m_mom.Abs2()<<")"<<std::endl;
else msg_Error()<<"No antitriplet: "<<p_anti<<" ";
if (p_trip!=NULL)
msg_Error()<<" Trip: "<<p_trip->m_mom<<" ("<<p_trip->m_mom.Abs2()<<")";
else
msg_Error()<<"No triplet: "<<p_trip<<" ";
if (p_anti!=NULL)
msg_Error()<<" Anti: "<<p_anti->m_mom<<" ("<<p_anti->m_mom.Abs2()<<")"<<std::endl;
else
msg_Error()<<"No antitriplet: "<<p_anti<<" ";
msg_Error()<<" diff: "<<testmom;
rest.Boost(m_momentum);
back.BoostBack(m_momentum);
//PRINT_VAR(m_momentum);
PRINT_VAR(m_momentum);
msg_Error()<<" from "<<newmom<<" --> "<<m_momentum<<"."<<std::endl;
//exit(1);
}
if (p_left!=NULL) {
msg_Error()<<"Maybe error in RescaleMomentum("<<save[0]<<" -> "<<m_momentum<<")"<<std::endl
......
......@@ -115,16 +115,7 @@ bool Dipole_Splitter::EmitGluon(Dipole * dip1,Dipole *& dip2) {
dip1->Triplet()->m_kt2max==dip1->AntiTriplet()->m_kt2max) first = true;
if (!p_tools->PrepareKinematics(dip1,first) || !p_tools->DetermineSplitting(dip1,first)) {
p_tools->SwapSpectatorAndSplitter(dip1);
if (dip1->Triplet()->m_info=='B' || dip1->AntiTriplet()->m_info=='B') {
msg_Debugging()<<"==============================================================="<<std::endl
<<"=== swapped splitter and spectator (now: m^2 = "<<dip1->Mass2()<<") for "
<<dip1->Triplet()->m_flav<<"("<<dip1->Triplet()->m_info<<") "
<<dip1->AntiTriplet()->m_flav<<"("<<dip1->AntiTriplet()->m_info<<")"<<std::endl
<<" "<<dip1->Triplet()->m_mom<<" (kt^2_max = "<<dip1->Triplet()->m_kt2max<<") "
<<dip1->AntiTriplet()->m_mom<<" (kt^2_max = "<<dip1->AntiTriplet()->m_kt2max<<")"
<<"."<<std::endl;
}
if (!p_tools->PrepareKinematics(dip1,false) || !p_tools->DetermineSplitting(dip1,false,false)) {
if (!p_tools->PrepareKinematics(dip1,first) || !p_tools->DetermineSplitting(dip1,first,false)) {
delete dip1->Triplet();
delete dip1->AntiTriplet();
delete dip1;
......
......@@ -56,7 +56,6 @@ void Hadronisation_Parameters::Init(string dir,string file)
p_doubletransitions = new Double_Transitions();
if (msg_LevelIsTracking())
p_doubletransitions->PrintDoubleTransitions();
//exit(1);
switch (int(m_parametermap[string("leading")])) {
case 3:
......
......@@ -82,6 +82,7 @@ Soft_Cluster_Handler::~Soft_Cluster_Handler()
bool Soft_Cluster_Handler::TreatClusterList(Cluster_List * clin, Blob * blob,bool cludec)
{
msg_Debugging()<<METHOD<<"-----------------------------------------------------"<<std::endl;
if (clin->size()==1) return TreatSingleCluster(*clin->begin(),blob);
if (cludec && clin->size()==2) return TreatClusterDecay(clin,blob);
......@@ -94,9 +95,9 @@ bool Soft_Cluster_Handler::TreatClusterList(Cluster_List * clin, Blob * blob,boo
<<METHOD<<" for "<<clin->size()<<" clusters."<<std::endl;
do {
if (!UpdateTransitions(clin)) {
msg_Debugging()<<"Error in "<<METHOD<<", E = "<<E<<" : "<<std::endl
<<" Cannot update any further. Leave and hope for the best."<<std::endl
<<(*clin)<<std::endl;
msg_Error()<<"Error in "<<METHOD<<", E = "<<E<<" : "<<std::endl
<<" Cannot update any further. Leave and hope for the best."<<std::endl
<<(*clin)<<std::endl;
return false;
}
} while (!CheckIfAllowed(clin,E));
......@@ -146,17 +147,22 @@ bool Soft_Cluster_Handler::TreatClusterDecay(Cluster_List * clin, Blob * blob)
{
if (!CheckListForTreatment(clin)) return true;
Cluster_Iterator cit=clin->begin();
Cluster * left((*cit++)), * right((*cit));
Cluster * cluster(right->GetPrev());
msg_Debugging()<<METHOD<<"------------------------------------------------------"<<std::endl
<<" for "<<clin->size()<<" clusters with masses "
<<left->Mass()<<" + "<<right->Mass()<<", original cluster: mass = "
<<cluster->Mass()<<"."<<std::endl;
double E(-1.);
while (!CheckIfAllowed(clin,E)) {
if (UpdateTransitions(clin)) continue;
Cluster_Iterator cit=clin->begin();
Cluster * left((*cit++)), * right((*cit));
if (left->GetPrev()!=right->GetPrev()) {
msg_Error()<<"Error in "<<METHOD<<" ("<<clin->size()<<" clusters) : "<<std::endl
<<" No common previous cluster."<<std::endl;
return false;
}
Cluster * cluster(left->GetPrev());
clin->clear();
if (!EnforcedDecay(cluster,blob,true,clin)) {
msg_Error()<<"Error in "<<METHOD<<" ("<<clin->size()<<" clusters) : "<<std::endl
......@@ -172,7 +178,7 @@ bool Soft_Cluster_Handler::TreatClusterDecay(Cluster_List * clin, Blob * blob)
bool Soft_Cluster_Handler::AttachHadronsToBlob(Cluster_List * clin,ATOOLS::Blob * blob)
{
msg_Tracking()<<"$$$$$$ "<<METHOD<<" for "<<clin->size()<<" clusters."<<std::endl;
msg_Debugging()<<"$$$$$$ "<<METHOD<<" for "<<clin->size()<<" clusters."<<std::endl;
Cluster_Iterator cit(clin->begin());
Particle * part;
Cluster * cluster;
......@@ -183,8 +189,8 @@ bool Soft_Cluster_Handler::AttachHadronsToBlob(Cluster_List * clin,ATOOLS::Blob
part = cluster->GetSelf();
part->SetFinalMass();
blob->AddToOutParticles(part);
msg_Tracking()<<"$$ attach one hadron ("<<part->Flav()<<", "<<part->Momentum()<<") "
<<"from cluster "<<cluster->Number()<<"."<<std::endl;
msg_Debugging()<<"$$ attach one hadron ("<<part->Flav()<<", "<<part->Momentum()<<") "
<<"from cluster "<<cluster->Number()<<"."<<std::endl;
delete cluster->GetTrip();
delete cluster->GetAnti();
delete cluster;
......@@ -272,22 +278,20 @@ int Soft_Cluster_Handler::CheckCluster(Cluster * cluster,bool lighter)
double decayweight(DecayWeight(cluster,haddec1,haddec2,false));
double transformweight(TransformWeight(cluster,hadtrans,lighter,false));
if (decayweight>0. || transformweight>0.) {
msg_Tracking()<<"@@@@ "<<METHOD<<" for cluster ("<<cluster->Number()<<":"
<<cluster->GetTrip()->m_flav<<"("<<cluster->GetTrip()->m_info<<"), "
<<cluster->GetAnti()->m_flav<<"("<<cluster->GetAnti()->m_info<<"), "
<<"m = "<<cluster->Mass()<<") --> ";
msg_Debugging()<<"@@@@ "<<METHOD<<" for cluster ("<<cluster->Number()<<":"
<<cluster->GetTrip()->m_flav<<"("<<cluster->GetTrip()->m_info<<"), "
<<cluster->GetAnti()->m_flav<<"("<<cluster->GetAnti()->m_info<<"), "
<<"m = "<<cluster->Mass()<<") --> ";
double totweight((decayweight+transformweight)*0.9999999);
if (totweight<=0. || decayweight/totweight>ran.Get()) {
m_decays += 1;
cluster->push_back(haddec1);
cluster->push_back(haddec2);
msg_Tracking()<<haddec1<<" + "<<haddec2<<"."<<std::endl;
return 2;
}
if (transformweight>0.) {
m_transitions += 1;
cluster->push_back(hadtrans);
msg_Tracking()<<hadtrans<<"."<<std::endl;
return 1;
}
}
......@@ -299,18 +303,17 @@ bool Soft_Cluster_Handler::CheckIfAllowed(Cluster_List * clin,double & E) {
double totmass(0.);
Vec4D totmom(0.,0.,0.,0.);
Cluster * cluster;
std::vector<Flavour> flavs;
for (Cluster_Iterator cit=clin->begin();cit!=clin->end();cit++) {
cluster = (*cit);
switch (cluster->size()) {
case 1:
totmass += (*cluster)[0].HadMass();
flavs.push_back((*cluster)[0]);
break;
case 2:
//if (cluster->Mass()<(*cluster)[0].HadMass()+(*cluster)[1].HadMass()) return false;
//break;
case 0:
default:
flavs.push_back(Flavour(kf_cluster));
totmass += cluster->Mass();
break;
}
......@@ -321,32 +324,50 @@ bool Soft_Cluster_Handler::CheckIfAllowed(Cluster_List * clin,double & E) {
}
bool Soft_Cluster_Handler::UpdateTransitions(Cluster_List * clin) {
//std::cout<<"@@@@ "<<METHOD<<" for "<<clin->size()<<" clusters. "<<std::endl;
Cluster * cluster, * winner(NULL);
Flavour hadron,winhad;
double tfwt, maxwt(0.);
Flavour hadron1,hadron2,winhad1,winhad2;
double wt, maxwt(0.);
int winno;
bool found(false);
for (Cluster_Iterator cit=clin->begin();cit!=clin->end();cit++) {
cluster = (*cit);
wt = -1.;
if (cluster->size()==1) {
hadron = (*cluster)[0];
tfwt = TransformWeight(cluster,hadron,true,true);
if (tfwt>maxwt) {
winner = cluster;
winhad = hadron;
maxwt = tfwt;
found = true;
hadron1 = (*cluster)[0];
wt = TransformWeight(cluster,hadron1,true,true);
if (wt>maxwt) {
winner = cluster;
winhad1 = hadron1;
winno = 1;
maxwt = wt;
found = true;
}
}
else if (wt<0. || cluster->size()==2) {
wt = TransformWeight(cluster,hadron1,hadron2,true);
if (wt>maxwt) {
winner = cluster;
winhad1 = hadron1;
winhad2 = hadron2;
winno = 2;
maxwt = wt;
found = true;
}
}
}
if (found) {
//std::cout<<"@@@@ "<<METHOD<<" found a winner: "<<winner<<"."<<std::endl;
//std::cout<<"@@@@ Replace "<<(*winner)[0]<<" with "<<winhad<<"."<<std::endl;
winner->clear();
winner->push_back(winhad);
return true;
if (winno==1) {
winner->clear();
winner->push_back(winhad1);
return true;
}
else if (winno==2) {
winner->clear();
winner->push_back(winhad1);
winner->push_back(winhad2);
return true;
}
}
//std::cout<<"@@@@ "<<METHOD<<": no winner found."<<std::endl;
return false;
}
......@@ -647,7 +668,7 @@ double Soft_Cluster_Handler::TransformKin(const double MC,const ATOOLS::Flavour
}
double Soft_Cluster_Handler::DecayWeight(Cluster * cluster,Flavour & had1,Flavour & had2,
const bool enforce)
const bool & enforce)
{
if (m_dweightmode==DecayWeight::off && !enforce) return 0.;
Flavour_Pair flpair;
......@@ -776,7 +797,7 @@ void Soft_Cluster_Handler::FixHHDecay(Cluster * cluster,Blob * blob,
if (constrained && m_pt2max<p2max) pt2max = m_pt2max;
if (cluster->GetTrip()->m_info=='L' || cluster->GetAnti()->m_info=='L') {
pt2max *= m_pt2maxfac * sqr(m_pt2max)/(Max(m_pt2max,m12)*Max(m_pt2max,m22));
pt2max *= sqr(m_pt2max)/(Max(m_pt2max,m12)*Max(m_pt2max,m22));
pt2min = m_pt2min/m12*m_pt2min/m22*m_pt2min;
if (pt2min>pt2max) pt2min = pt2max/4.;
pt2 = p_as->SelectPT(pt2max,pt2min);
......@@ -893,14 +914,14 @@ bool Soft_Cluster_Handler::TryLocalCompensation(Cluster_List * clin)
}
}
if (!partner) {
//msg_Out()<<"@@@ "<<METHOD<<": no partner for cluster "<<cluster->Number()<<"."<<std::endl;
//msg_Tracking()<<"@@@ "<<METHOD<<": no partner for cluster "<<cluster->Number()<<"."<<std::endl;
return false;
}
mass1 = (*cluster)[0].HadMass();
mass2 = partner->size()==1?(*partner)[0].HadMass():partner->Mass();
if (sqr(mass1+mass2)>(cluster->Momentum()+partner->Momentum()).Abs2()) {
if (direx==0) {
//msg_Out()<<"@@@ "<<METHOD<<": Tried "<<partner->Number()
//msg_Tracking()<<"@@@ "<<METHOD<<": Tried "<<partner->Number()
// <<" for cluster "<<cluster->Number()<<", not enough energy."<<std::endl;
return false;
}
......@@ -908,7 +929,7 @@ bool Soft_Cluster_Handler::TryLocalCompensation(Cluster_List * clin)
if (direx== 1) partner = cluster->GetNBAnti();
mass2 = partner->size()==1?(*partner)[0].HadMass():partner->Mass();
if (sqr(mass1+mass2)>(cluster->Momentum()+partner->Momentum()).Abs2()) {
//msg_Out()<<"@@@ "<<METHOD<<": Tried "<<partner->Number()
//msg_Tracking()<<"@@@ "<<METHOD<<": Tried "<<partner->Number()
// <<" for cluster "<<cluster->Number()<<", not enough energy."<<std::endl;
return false;
}
......@@ -920,7 +941,7 @@ bool Soft_Cluster_Handler::TryLocalCompensation(Cluster_List * clin)
momenta.push_back(cluster->Momentum());
momenta.push_back(partner->Momentum());
if (!hadpars.AdjustMomenta(2,&momenta.front(),&masses.front())) {
//msg_Out()<<"@@@ "<<METHOD<<": Tried "<<partner->Number()
//msg_Tracking()<<"@@@ "<<METHOD<<": Tried "<<partner->Number()
// <<" for cluster "<<cluster->Number()<<", adjusting didn't work."<<std::endl;
return false;
}
......
......@@ -51,13 +51,12 @@ namespace AHADIC {
int CheckCluster(Cluster *,bool);
double TransformWeight(Cluster *,ATOOLS::Flavour &,
const bool & lighter=false,
const bool & enforce=false);
const bool & lighter=false,const bool & enforce=false);
double TransformKin(const double MC,const ATOOLS::Flavour & flav,
const bool & enforce=false);
double DecayWeight(Cluster *,ATOOLS::Flavour &,ATOOLS::Flavour &,
const bool enforce=false);
const bool & enforce=false);
bool ClusterAnnihilation(Cluster * cluster,ATOOLS::Flavour &,ATOOLS::Flavour &);
void FixHHDecay(Cluster * cluster,ATOOLS::Blob * blob,
const ATOOLS::Flavour had1,const ATOOLS::Flavour had2,
......
......@@ -152,14 +152,15 @@ void Splitting_Tools::SetSpectatorAndSplitter(Dipole * dip) {
void Splitting_Tools::SwapSpectatorAndSplitter(Dipole * dip) {
Proto_Particle * help(p_split); p_split = p_spect; p_spect = help;
if (dip->IsSwitched()) dip->SetSwitched(false);
else dip->SetSwitched(true);
else dip->SetSwitched(true);
if (p_split->m_kt2max<=0.)
p_split->m_kt2max = p_spect->m_kt2max = ATOOLS::Max(p_split->m_kt2max,p_spect->m_kt2max);
}
bool Splitting_Tools::PrepareKinematics(Dipole * dip,const bool & first,const bool & enforce) {
if (dip->MassBar2()<4.*m_mmin_2) {
msg_Debugging()<<METHOD<<"(massmin^2 = "<<m_mmin_2<<", red.mass^2 = "<<dip->MassBar2()<<") :"<<std::endl
msg_Debugging()<<"~~ "<<METHOD<<"(massmin^2 = "<<m_mmin_2<<", red.mass^2 = "<<dip->MassBar2()<<") :"
<<std::endl
<<" --> Dipole cannot split."<<std::endl;
//(*dip).Output();
return false;
......@@ -250,32 +251,35 @@ bool Splitting_Tools::SelectFlavour(const bool & vetodiquark)
bool Splitting_Tools::DetermineSplitting(Dipole * dip1,const bool & first,const bool & vetodiquark) {
int trials(0);
msg_Debugging()<<"In "<<METHOD<<"(first = "<<first<<", veto = "<<vetodiquark<<")."<<std::endl;
//dip1->Output();
msg_Debugging()<<"~~ "<<METHOD<<"(first = "<<first<<", veto = "<<vetodiquark<<")."<<std::endl;
while (trials<1000) {
trials++;
if (ProduceKinematics(first,vetodiquark) &&
ConstructKinematics()) break;
}
if (m_analyse) {
Histogram * histo = (m_histograms.find(std::string("Splitting_Trials")))->second;
histo->Insert(trials);
ConstructKinematics()) {
if (m_analyse) {
Histogram * histo = (m_histograms.find(std::string("Splitting_Trials")))->second;
histo->Insert(trials);
}
msg_Debugging()<<"~~ "<<METHOD<<"(first = "<<first<<", "
<<"veto = "<<vetodiquark<<") -> succeeded."<<std::endl;
return true;
}
}
msg_Debugging()<<"In "<<METHOD<<"(first = "<<first<<", veto = "<<vetodiquark<<") -> succeeded."<<std::endl;
return true;
msg_Debugging()<<"~~ "<<METHOD<<"(first = "<<first<<", "
<<"veto = "<<vetodiquark<<") -> failed."<<std::endl;
return false;
}
bool Splitting_Tools::ProduceKinematics(const bool & first,const bool & vetodiquark) {
msg_Debugging()<<"In "<<METHOD<<"."<<std::endl;
msg_Debugging()<<"~~ "<<METHOD<<"."<<std::endl;
if (!SelectFlavour(vetodiquark)) {
msg_Error()<<"no flavour selected, return false."<<std::endl;
return false;
}
msg_Debugging()<<" flav = "<<m_flav<<"."<<std::endl;
if (!FixRanges(first)) return false;
msg_Debugging()<<" kt2min = "<<m_kt2_min<<", "
<<"z in ["<<m_zmin<<", "<<m_zmax<<"]."<<std::endl;
msg_Debugging()<<" kt2min = "<<m_kt2_min<<", z in ["<<m_zmin<<", "<<m_zmax<<"]."<<std::endl;
m_kt2 = m_z = -1.;
int trials(0);
......@@ -298,8 +302,8 @@ bool Splitting_Tools::ProduceKinematics(const bool & first,const bool & vetodiqu
}
}
m_phi = 2.*M_PI*ran.Get();
msg_Debugging()<<"___ "<<METHOD<<" : exit inner loop with kt2 = "<<m_kt2<<" "
<<"and z = "<<m_z<<" for "<<m_flav<<"."<<std::endl;
msg_Debugging()<<"~~ "<<METHOD<<" : exit inner loop with kt2 = "<<m_kt2<<" "
<<"and z = "<<m_z<<" for "<<m_flav<<"."<<std::endl;
if (m_analyse) {
Histogram * histo = (m_histograms.find(std::string("Kinematics_Trials")))->second;
......@@ -359,18 +363,19 @@ bool Splitting_Tools::UpdateRanges(const bool & first) {
else {
if ((p_spect->m_info=='B')) {
if (p_split->m_info=='B')
m_kt2_max = Min(m_kt2_max,m_pt2max_factor*sqrt(4.*p_spect->m_mom.PPerp2()*p_split->m_mom.PPerp2()));
m_kt2_max = Min(m_kt2_max,sqrt(m_pt2max_factor*p_spect->m_mom.PPerp2()*p_split->m_mom.PPerp2()));
else
m_kt2_max = Min(m_kt2_max,m_pt2max_factor*sqrt(4.*p_split->m_mom.PPerp2()*m_pt2min));
m_kt2_max = Min(m_kt2_max,sqrt(m_pt2max_factor*p_split->m_mom.PPerp2()*m_pt2min));
}
else {
double pt2(m_pt2max_factor*p_spect->m_mom.PPerp2(p_split->m_mom));
double pt2((m_glusplit?m_pt2max_factor:1.)*p_spect->m_mom.PPerp2(p_split->m_mom));
if (m_kt2_max>pt2) m_kt2_max = sqrt(m_kt2_max*pt2);
}
}
msg_Tracking()<<"... first: "<<p_split->m_info<<p_spect->m_info
<<"("<<p_split->m_flav<<" "<<p_spect->m_flav<<"): "
<<"kt_max = "<<sqrt(m_kt2_max)<<" from pt(kin) = "<<p_spect->m_mom.PPerp2(p_split->m_mom)<<", "
<<"kt_max = "<<sqrt(m_kt2_max)<<" from pt(kin) = "
<<p_spect->m_mom.PPerp2(p_split->m_mom)<<", "
<<"p_min = "<<Min(p_split->m_mom.PSpat(),p_spect->m_mom.PSpat())<<" and "
<<"cos = "<<p_split->m_mom.CosTheta(p_spect->m_mom)
<<", red. mass = "<<sqrt(m_Qt2)<<"."<<std::endl;
......@@ -441,7 +446,7 @@ void Splitting_Tools::AftermathOfSplitting(Dipole * dip1) {
}
}
msg_Debugging()<<"___ "<<METHOD<<" yields new particles : "<<std::endl
<<(*p_out1)<<(*p_out2)<<std::endl;
<<(*p_out1)<<(*p_out2)<<std::endl;
}
else {
p_split->m_mom = m_mom3;
......@@ -575,14 +580,14 @@ bool Splitting_Tools::ConstructKinematics() {
Vec4D q23 = m_mom0-q1;
if (sqr(q23*q1)-m_s23*m_m1_2<0.) {
//msg_Error()<<"Warning in "<<METHOD<<": gamma explodes: "<<(sqr(q23*q1)-m_s23*m_m1_2)
// <<" for kt^2 = "<<m_kt2<<"."<<std::endl;
msg_Debugging()<<"Warning in "<<METHOD<<": gamma explodes: "<<(sqr(q23*q1)-m_s23*m_m1_2)
<<" for kt^2 = "<<m_kt2<<"."<<std::endl;
return false;
}
double gam = q23*q1+sqrt(sqr(q23*q1)-m_s23*m_m1_2);
if (IsZero(m_s23*m_m1_2-gam*gam) || IsZero(gam)) {
msg_Error()<<"Warning in "<<METHOD<<": beta explodes: "<<(m_s23*m_m1_2-gam*gam)
<<" for kt^2 = "<<m_kt2<<" & gamma = "<<gam<<"."<<std::endl;
msg_Debugging()<<"Warning in "<<METHOD<<": beta explodes: "<<(m_s23*m_m1_2-gam*gam)
<<" for kt^2 = "<<m_kt2<<" & gamma = "<<gam<<"."<<std::endl;
return false;
}
double a23 = m_s23/gam, a1=m_m1_2/gam, beta=1.0/(1.0-a23*a1);
......@@ -590,9 +595,9 @@ bool Splitting_Tools::ConstructKinematics() {
Vec4D n = beta*(q1-a1*q23);
if (IsZero(a1) || IsZero(1.-m_y) && IsZero(gam) || IsZero((1.+a23*a1)) ||
IsZero((1./a1-(m_m2_2+m_m3_2)/gam))) {
msg_Error()<<"Warning in "<<METHOD<<": z explodes: "<<std::endl
<<" a23 = "<<a23<<", a1 = "<<a1<<", y "<<m_y<<" gamma = "<<gam
<<" & masses = "<<m_m2_2<<"/"<<m_m3_2<<"."<<std::endl;
msg_Debugging()<<"Warning in "<<METHOD<<": z explodes: "<<std::endl
<<" a23 = "<<a23<<", a1 = "<<a1<<", y "<<m_y<<" gamma = "<<gam
<<" & masses = "<<m_m2_2<<"/"<<m_m3_2<<"."<<std::endl;
return false;
}
m_z =
......@@ -611,16 +616,16 @@ bool Splitting_Tools::ConstructKinematics() {
}
q2 = m_mom0-q3-q1;
msg_Debugging()<<"___ "<<METHOD<<": before boosting back."<<std::endl
<<"___ "<<m_mom0<<" = "<<m_mom1<<" ("<<sqrt(Max(0.,m_mom1.Abs2()))<<")"
msg_Debugging()<<"~~ "<<METHOD<<": before boosting back."<<std::endl
<<"~~ "<<m_mom0<<" = "<<m_mom1<<" ("<<sqrt(Max(0.,m_mom1.Abs2()))<<")"
<<" + "<<m_mom3<<" ("<<sqrt(Max(0.,m_mom3.Abs2()))<<")"<<std::endl
<<"___ "<<q1<<" ("<<sqrt(Max(0.,q1.Abs2()))<<")"
<<"~~ "<<q1<<" ("<<sqrt(Max(0.,q1.Abs2()))<<")"
<<" + "<<q2<<" ("<<sqrt(Max(0.,q2.Abs2()))<<")"
<<" + "<<q3<<" ("<<sqrt(Max(0.,q3.Abs2()))<<")."<<std::endl;
if (m_m1_2==0.) {
msg_Debugging()<<"___ "<<METHOD<<" summary: "<<std::endl
<<" Masses: m_m2 = m_m3 = "<<m_m2<<" = "<<m_m3<<" --> s23 = "<<m_m23_2<<"."<<std::endl;
<<" Masses: m_m2 = m_m3 = "<<m_m2<<" = "<<m_m3<<" --> s23 = "<<m_m23_2<<"."<<std::endl;
}
if (q1[0]<0. || q2[0]<0. || q3[0]<0.) {
......@@ -648,7 +653,7 @@ bool Splitting_Tools::ConstructKinematics() {
m_cms.BoostBack(m_mom0);
msg_Debugging()<<"___ "<<METHOD<<"(out): succeeded to build kinematics "
<<"for kt^2 = "<<m_kt2<<", z = "<<m_z<<" for "<<m_flav<<"."<<std::endl;
<<"for kt^2 = "<<m_kt2<<", z = "<<m_z<<" for "<<m_flav<<"."<<std::endl;
return true;
}
......
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