Commit 759970b3 authored by Stefan Hoeche's avatar Stefan Hoeche

svn merge -r22957:22959 branches/nlomerging_comix

svn merge -r22984:23105 branches/nlomerging_comix
parent 82318aa6
......@@ -9,7 +9,7 @@ using namespace std;
Cluster_Decay_Handler::Cluster_Decay_Handler(Cluster_List * clulist,bool ana) :
p_softclusters(hadpars->GetSoftClusterHandler()),
p_clus(new Cluster_Part(hadpars->GetSplitter(),ana)),
p_clus(new Cluster_Part(ana)),
p_clulist(clulist),
p_analysis(ana?new Cluster_Decay_Analysis():NULL)
{ }
......@@ -25,82 +25,33 @@ Cluster_Decay_Handler::~Cluster_Decay_Handler()
int Cluster_Decay_Handler::DecayClusters(Blob * blob)
{
Cluster * cluster;
Cluster_List clist;
msg_Tracking()<<"+++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl
<<"+++ "<<METHOD<<" for "<<p_clulist->size()<<" clusters."<<std::endl;
Cluster_Iterator cit(p_clulist->begin());
while (!p_clulist->empty()) {
cluster = p_clulist->front();
if (!cluster->Active()) {
msg_Error()<<"Error in "<<METHOD<<" inactive cluster:"<<std::endl;
cluster->Print();
msg_Error()<<" Will continue and hope for the best."<<std::endl;
return -1;
}
msg_Tracking()
<<"+++ Test cluster ["
<<cluster->GetTrip()->m_flav<<"("<<cluster->GetTrip()->m_info<<"), "
<<cluster->GetAnti()->m_flav<<"("<<cluster->GetAnti()->m_info<<")].\n";
if (!cluster->Active()) return -1;
if (p_clus->TestDecay(cluster)) {
clist.empty();
clist.push_back(cluster->GetLeft());
clist.push_back(cluster->GetRight());
//if (cluster->GetTrip()->m_info=='B' ||
// cluster->GetAnti()->m_info=='B')
// msg_Out()<<"++++ From "<<cluster->Number()
// <<"("<<cluster->Mass()<<", "<<cluster->Momentum()<<") "
// <<cluster->GetLeft()->Number()
// <<" ("<<cluster->GetLeft()->Mass()<<") + "
// <<cluster->GetRight()->Number()
// <<" ("<<cluster->GetRight()->Mass()<<"), "
// <<"popped "<<cluster->GetLeft()->GetAnti()->m_flav<<"\n"
// <<(*cluster);
if (!p_softclusters->TreatClusterDecay(&clist,blob)) {
msg_Error()<<"Error in "<<METHOD<<" : "<<std::endl
Cluster_List * clist(cluster->GetClusters());
if (!p_softclusters->TreatClusterDecay(clist,blob)) {
msg_Error()<<"Error in "<<METHOD<<" : \n"
<<" Did not find a kinematically allowed "
<<"solution for the cluster list."<<std::endl
<<" Will trigger retrying the event."<<std::endl;
<<"solution for the cluster list.\n"
<<" Will trigger retrying the event.\n";
return -1;
}
while (!clist.empty()) {
p_clulist->push_back(clist.front());
clist.pop_front();
while (!clist->empty()) {
p_clulist->push_back(clist->front());
clist->pop_front();
}
}
else {
msg_Tracking()<<"+++ TestDecay did not work out - enforce decay.\n";
if (!p_softclusters->EnforcedDecay(cluster,blob,true)) {
msg_Error()<<"+++ EnforcedDecay did not work out, will return -1.\n";
return -1;
}
msg_Tracking()<<"+++ EnforcedDecay did work out, will continue.\n";
if (!p_softclusters->EnforcedDecay(cluster,blob,true)) return -1;
}
delete (p_clulist->front()->GetTrip());
delete (p_clulist->front()->GetAnti());
delete (p_clulist->front());
p_clulist->pop_front();
msg_Tracking()<<"++++ After popping front cluster list has now size = "
<<p_clulist->size()<<"."<<std::endl;
}
if (p_analysis) p_analysis->AnalyseThis(blob);
return 1;
}
ATOOLS::Blob * Cluster_Decay_Handler::
ClusterDecayBlob(Cluster * cluster,Cluster_List * p_clulist) {
Blob * decblob(cluster->ConstructDecayBlob());
if (cluster->GetLeft()!=NULL &&
cluster->GetLeft()->GetFlav()==Flavour(kf_cluster)) {
p_clulist->push_back(cluster->GetLeft());
}
if (cluster->GetRight()!=NULL &&
cluster->GetRight()->GetFlav()==Flavour(kf_cluster)) {
p_clulist->push_back(cluster->GetRight());
}
if (cluster) cluster->SetActive(false);
return decblob;
}
......@@ -2,7 +2,6 @@
#define AHADIC_Decays_Cluster_Decay_Handler_H
#include "AHADIC++/Decays/Cluster_Part.H"
#include "AHADIC++/Tools/Dipole_Splitter.H"
#include "AHADIC++/Tools/Soft_Cluster_Handler.H"
#include "ATOOLS/Phys/Blob_List.H"
#include "AHADIC++/Tools/Cluster.H"
......@@ -19,91 +18,11 @@ namespace AHADIC {
Cluster_Decay_Analysis * p_analysis;
ATOOLS::Blob * DecayIt(Cluster *);
ATOOLS::Blob * ClusterDecayBlob(Cluster *,Cluster_List *);
public:
Cluster_Decay_Handler(Cluster_List *,bool=false);
~Cluster_Decay_Handler();
int DecayClusters(ATOOLS::Blob *);
};
/*!
\file Cluster_Decay_Handler.H
\brief Contains the class AHADIC::Cluster_Formeration_Handler.
*/
/*!
\class Cluster_Decay_Handler
\brief The Cluster_Decay_Handler organises the decay of clusters at the end of
the parton shower phase.
The central (and also the only public) method of this class is
Cluster_Decay_Handler::DecayClusters(Cluster_List *,ATOOLS::Blob_List *),
which is being called, once the previous cluster formation from partons,
managed by the class Cluster_Formation_Handler, is finished.
*/
/*!
\var Soft_Cluster_Handler * p_softclusters
\brief A pointer to the common Soft_Cluster_Handler, used in both the formation
and decay phase to decide on forced cluster to hadron transitions and decays
and to perform them.
*/
/*!
\var Cluster_Part * p_clus
\brief Pointer to the Cluster_Part, responsible for the definition of the actual
cluster decay kinematics.
*/
/*!
\var Cluster_Decay_Analysis * p_analysis
\brief Pointer to the Cluster_Decay_Analysis.
*/
/*!
\fn ATOOLS::Blob * Cluster_Decay_Handler::DecayIt(Cluster *)
\brief Managing for the decay of one single cluster.
The following steps are taken:
- a decay blob (status: blob_status::needs_hadrondecays) is initialised
and the cluster is added as incoming particle;
- the method Cluster_Part::TestDecay(Cluster*) decides upon a flavour to be
popped out of the vacuum and instantiates two new clusters with according
flavour content. They are then linked as left and right pointers to the
decaying one;
- finally, the cluster decay blob (without the two outgoing clusters added)
is returned.
*/
/*!
\fn Cluster_Decay_Handler::Cluster_Decay_Handler(Soft_Cluster_Handler *,bool=false)
\brief Constructor of the class. Initialises a new Cluster_Part and, eventually,
a Cluster_Decay_Analysis object.
The boolean steers, whether the Cluster_Decay_Analysis object is instantiated.
*/
/*!
\fn Cluster_Decay_Handler::~Cluster_Decay_Handler()
\brief Deletes the Cluster_Part and the Cluster_Decay_Analysis object.
*/
/*!
\fn int Cluster_Decay_Handler::DecayClusters(Cluster_List *,ATOOLS::Blob_List * = NULL)
\brief Managing the decays of the clusters in a Cluster_List into secondary
clusters and, ultimately, the primordial hadrons.
The idea is the following: each decaying cluster is connected with a cluster list,
filled by its offsprings. These emerge from the call to the method
Cluster_Decay_Handler::DecayIt(Cluster *) and, at that point, are still clusters.
Then the Soft_Cluster_Handler tests the two offspring clusters, whether they must
decay or transit into hadrons. This is done through the method
Soft_Cluster_Handler::TreatClusterList(Cluster_List *,ATOOLS::Blob *). Now the
offspring list is iterated over. Each of them is treated in the following way:
- First, the ATOOLS::Particle related to it is added as outgoing particle to the decay
ATOOLS::Blob of the original cluster.
- If the offspring does not have any filled left and right pointers, it did not decay into
hadrons and hence it can safely be added to the original cluster list for further treatment.
- If, on the other hand, the offspring does have filled left and right pointers, it decayed
into hadrons or clusters. Then the offsprings of the offspring are tested: If they are
not hadrons, they are added to the original cluster list for further treatment.
Finally, the original cluster is deleted and taken out of the original cluster list.
This procedure is repeated until the original cluster list is empty, i.e., until all
clusters are decayed into hadrons.
*/
}
#endif
......@@ -7,46 +7,23 @@ using namespace AHADIC;
using namespace ATOOLS;
using namespace std;
Cluster_Part::Cluster_Part(Dipole_Splitter * splitter,bool ana) :
m_ana(ana), p_splitter(splitter), m_fails(0), m_att(0)
Cluster_Part::Cluster_Part(bool ana) :
m_ana(ana), m_fails(0), m_att(0)
{
if (m_ana) {
m_histograms[string("PT_Cluster")] = new Histogram(0,0.,1.5,150);
}
p_csplitter = new Cluster_Splitter();
}
Cluster_Part::~Cluster_Part()
{
if (m_ana) {
msg_Tracking()<<METHOD<<" (had "<<m_fails<<" fails in "<<m_att<<" attempts).\n";
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();
}
if (p_csplitter) delete p_csplitter;
}
bool Cluster_Part::TestDecay(Cluster * const cluster)
{
m_att++;
if (!p_splitter->SplitCluster(cluster)) {
if (p_csplitter && !(*p_csplitter)(cluster)) {
m_fails++;
return false;
}
if (m_ana) {
Vec4D lmom(cluster->GetLeft()->Momentum());
double pt = sqrt(sqr(lmom[1]) + sqr(lmom[2]));
Histogram* histo((m_histograms.find(std::string("PT_Cluster")))->second);
histo->Insert(pt);
}
#ifdef AHAmomcheck
cluster->CheckConsistency(msg_Error(),METHOD);
#endif
return true;
}
......@@ -3,7 +3,7 @@
#include "AHADIC++/Tools/Cluster.H"
#include "AHADIC++/Tools/Dipole_Splitter.H"
#include "AHADIC++/Tools/Cluster_Splitter.H"
#include "ATOOLS/Math/Histogram.H"
#include "AHADIC++/Tools/Hadronisation_Parameters.H"
#include "ATOOLS/Org/Return_Value.H"
......@@ -14,120 +14,16 @@
namespace AHADIC {
class Cluster_Part {
protected:
bool m_ana;
Dipole_Splitter * p_splitter;
long int m_fails, m_att;
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(Dipole_Splitter * =NULL,bool=false);
Cluster_Part(bool=false);
~Cluster_Part();
bool TestDecay(Cluster * const);
};
/*!
\file Cluster_Part.H
\brief This file contains the class Cluster_Part.
*/
/*!
\class Cluster_Part
\brief The class cluster part is responsible for the decays of clusters into test clusters.
These are either kept as clusters for further decays or they are transformed into hadrons by
the Soft_Cluster_Handler.
*/
/*!
\var int Cluster_Part::m_4Qmode
\brief Up to now: 0 by default.
\todo: Think about improvement.
*/
/*!
\var int Cluster_Part::m_ybar_mode
\brief Up to now: 0 by default. For details, cf.
Cluster_Part::GetYBar(const double,const double)
*/
/*!
\var double Cluster_Part::m_fraction
\brief The variable \f$f\f$ in Cluster_Part::BuildKinematics(Cluster * const,ATOOLS::Flavour &).
It is set by the parameter "CLUSTER_MASS_FRACTION", default is 0.5.
*/
/*!
\var double Cluster_Part::m_ystar_expvalue
\brief The expectation value of the \f$y^*\f$ shift in terms of \f$y^{*,max}\f$.
For more details see the description of
Cluster_Part::GetYStar(const double,const double), there it is named
\f$\langle\delta y^*\rangle\f$.
It is set by the parameter "<Y*_SHIFT>", default is 0.5.
*/
/*!
\var double Cluster_Part::m_ystar_sigma
\brief The width of the \f$y^*\f$ distribution in a Gaussian.
For more details see the description of
Cluster_Part::GetYStar(const double,const double), there it is named \f$\sigma_y\f$.
It is set by the parameter "<Y*_WIDTH>", default is 0.5.
*/
/*!
\var ATOOLS::Vec4D Cluster_Part::m_momenta[4]
\brief The four momenta of the cluster constituents and the popped quark pair.
*/
/*!
\var bool Cluster_Part::m_ana
\brief Switches the analysis on and off. Default is true.
*/
/*!
\var std::map<std::string,ATOOLS::Histogram *> Cluster_Part::m_histograms
\brief A number of histograms monitoring the decay kinematics.
In particular, the \f$p_\perp\f$ distribution and the centre-of-mass energy squared of
the popped quarks is monitored, as well as the flavour composition of the popped quarks.
In addition, \f$y^*\f$, \f$y^{*}/y^{*, max}\f$, \f$\bar y\f$, \f$\bar y/\bar y^{max}\f$
are monitored.
*/
/*!
\fn Cluster_Part::Cluster_Part()
\brief Initialises the constants and, if neccessary, the histograms.
*/
/*!
\fn Cluster_Part::~Cluster_Part()
\brief Outputs the histograms from the analysis.
*/
/*!
\fn bool Cluster_Part::TestDecay(Cluster *)
\brief Decides the popped flavour and instantiates two clusters accordingly.
The procedure is as follows:
- First the decaying cluster is boosted in its c.m. frame and its constitutents
are oriented along the z-axis;
- then Cluster_Part::BuildKinematics(Cluster * const,ATOOLS::Flavour &) is invoked
to pop a flavour and fix the decay kinematics;
- finally, offspring clusters are constructed with flavour content and momenta
(eventually these clusters may become hadrons through the Soft_Cluster_Handler);
- the system is boosted and rotated back, the offsprings are linked to the decaying
cluster, corresponding particles are instantiated for them.
*/
/*!
\fn bool Cluster_Part::UpdateDecay(Cluster *,const int)
\brief Updates decay kinematics after one or both outgoing particles have been
shifted to a new mass-shell.
The basic idea is that the original decay is boosted in its c.m. system, with orientation
of the original clusters constituents along the z-axis. The energies and longitudinal (z-)
momenta of the outgoing particles are then recalculated with the new shifted masses of
the offsrpings. Their \f$p_\perp\f$ is kept and eventually rotated around the z-axis.
Finally the full system is boosted and rotated back in its original frame.
*/
/*!
\fn bool Cluster_Part::BuildKinematics(Cluster * const,ATOOLS::Flavour &)
\brief Builds the kinematics of a cluster->cluster+cluster decay in its own
rest frame.
The algorithm is as follows:
The procedure is repeated, until physical energies for the original cluster
constituents (larger than their mass) emerge.
*/
}
#endif
......@@ -22,7 +22,7 @@ Cluster_Formation_Handler::Cluster_Formation_Handler(Cluster_List* clulist,
bool ana) :
m_single_cr(hadpars->Get(string("colour_reconnections"))==1),
m_double_cr(false),
p_gludecayer(new Gluon_Decayer(hadpars->GetSplitter(),ana)),
p_gludecayer(new Gluon_Decayer(ana)),
p_cformer(new Cluster_Former()),
p_recons(new Colour_Reconnections(2,1,hadpars->Get(string("pt02")))),
p_softclusters(hadpars->GetSoftClusterHandler()),
......@@ -69,11 +69,14 @@ int Cluster_Formation_Handler::FormClusters(Blob * blob) {
m_partlists.clear();
m_clulists.clear();
}
//msg_Out()<<"\n\n\n\n"
// <<"====================================================\n"
// <<"====================================================\n"
// <<"====================================================\n"
// <<"In "<<METHOD<<": hadronize "<<blob->NInP()<<" partons.\n";
// Vec4D cms(0.,0.,0.,0.);
// for (size_t i=0;i<blob->NInP();i++) cms += blob->InParticle(i)->Momentum();
// msg_Out()<<"\n\n\n\n"
// <<"====================================================\n"
// <<"====================================================\n"
// <<"====================================================\n"
// <<"In "<<METHOD<<": hadronize "<<blob->NInP()<<" partons "
// <<"with E = "<<sqrt(cms.Abs2())<<".\n"<<(*blob)<<"\n";
if (!ExtractSinglets(blob)) { Reset(); return -1; }
if (!ShiftOnMassShells()) { Reset(); return -1; }
if (!FormOriginalClusters()) { Reset(); return -1; }
......@@ -82,8 +85,21 @@ int Cluster_Formation_Handler::FormClusters(Blob * blob) {
if (!ClustersToHadrons(blob)) {
Reset(); return -1;
}
msg_Tracking()<<METHOD<<":\n"<<(*p_clulist)<<"\n";
//msg_Out()<<(*blob)<<"\n";
// Vec4D blobmom(0.,0.,0.,0.);
// for (size_t i(0);i<blob->NOutP();i++)
// blobmom+=blob->OutParticle(i)->Momentum();
// msg_Out()<<"____________________________________________________\n"
// <<"____________________________________________________\n"
// <<"Cluster list after all merging etc.:\n"<<(*p_clulist)
// <<"____________________________________________________\n"
// <<"Blob momentum: "<<blobmom<<";\n"<<(*blob)<<"\n"
// <<"____________________________________________________\n"
// <<"____________________________________________________\n";
// msg_Out()<<"##################################################\n"
// <<METHOD<<" was successful: "
// <<p_clulist->size()<<" clusters left to decay.\n"
// <<(*p_clulist)<<"\n"
// <<"##################################################\n";
return 1;
}
......@@ -114,7 +130,6 @@ bool Cluster_Formation_Handler::ExtractSinglets(Blob * blob)
Proto_Particle_List * pli(NULL);
bool construct(false);
unsigned int col1(0), col2(0);
leading::code leading(hadpars->GetSplitter()->Leading());
Particle * part(NULL);
for (int i=0;i<blob->NInP();i++) {
part = blob->InParticle(i);
......@@ -126,7 +141,7 @@ bool Cluster_Formation_Handler::ExtractSinglets(Blob * blob)
Proto_Particle * copy =
new Proto_Particle(part->Flav(),part->Momentum(),
(part->Info()=='B')?'B':'L');
SetInfoTagForPrimaryParticle(copy,leading);
SetInfoTagForPrimaryParticle(copy);
pli->push_back(copy);
col1 = part->GetFlow(1);
if (col1==col2) construct = false;
......@@ -147,7 +162,7 @@ bool Cluster_Formation_Handler::ExtractSinglets(Blob * blob)
Proto_Particle * copy =
new Proto_Particle(part->Flav(),part->Momentum(),
part->Info()=='B'?'B':'L');
SetInfoTagForPrimaryParticle(copy,leading);
SetInfoTagForPrimaryParticle(copy);
pli->push_back(copy);
m_partlists.push_back(pli);
construct = true;
......@@ -157,25 +172,9 @@ bool Cluster_Formation_Handler::ExtractSinglets(Blob * blob)
}
void Cluster_Formation_Handler::
SetInfoTagForPrimaryParticle(Proto_Particle * proto,
leading::code lead) const {
SetInfoTagForPrimaryParticle(Proto_Particle * proto) const {
if (proto->m_info=='B') return;
switch (lead) {
case leading::quarks_and_gluons:
break;
case leading::quarks_and_gluons2:
break;
case leading::only_quarks:
if (proto->m_flav.IsQuark() || proto->m_flav.IsDiQuark())
proto->m_info='L';
else
proto->m_info='l';
break;
case leading::none:
default:
proto->m_info='l';
break;
}
proto->m_info=(proto->m_flav.IsQuark() || proto->m_flav.IsDiQuark())?'L':'l';
}
bool Cluster_Formation_Handler::ShiftOnMassShells() {
......@@ -401,40 +400,7 @@ bool Cluster_Formation_Handler::MergeClusterListsIntoOne() {
bool Cluster_Formation_Handler::ClustersToHadrons(Blob * blob)
{
#ifdef AHAmomcheck
Vec4D checkbef(0.,0.,0.,0.);
for (Cluster_Iterator cit=p_clulist->begin();cit!=p_clulist->end();cit++) {
checkbef += (*cit)->Momentum();
}
#endif
if (!p_softclusters->TreatClusterList(p_clulist,blob)) {
msg_Error()<<"Error in "<<METHOD<<" : \n"
<<" Did not find a kinematically allowed solution for the "
<<"cluster list.\n"
<<" Will trigger a new event.\n";
msg_Out()<<(*blob)<<"\n";
exit(1);
return false;
}
#ifdef AHAmomcheck
Vec4D checkaft(0.,0.,0.,0.);
for (Cluster_Iterator cit=p_clulist->begin();cit!=p_clulist->end();cit++) {
checkaft += (*cit)->Momentum();
}
for (short int i=0;i<blob->NOutP();i++) {
checkaft += blob->OutParticle(i)->Momentum();
}
double Q2(dabs((checkbef-checkaft).Abs2()));
if (Q2>1.e-12 || IsNan(Q2)) {
msg_Error()<<METHOD<<" yields a momentum violation : \n"
<<" "<<checkbef<<" - "<<checkaft<<" --> "
<<(checkbef-checkaft).Abs2()<<".\n";
}
else msg_Tracking()<<METHOD<<" conserves momentum:"
<<checkbef<<" --> "<<checkaft<<".\n";
#endif
if (!p_softclusters->TreatClusterList(p_clulist,blob)) return false;
if (m_analyse) {
Histogram * histomass, * histonumb;
......
This diff is collapsed.
#ifndef AHADIC_Formation_Gluon_Decayer_H
#define AHADIC_Formation_Gluon_Decayer_H
#include "AHADIC++/Tools/Dipole_Splitter.H"
#include "AHADIC++/Tools/Cluster.H"
#include "AHADIC++/Tools/Gluon_Splitter.H"
#include "AHADIC++/Tools/Dipole.H"
#include "ATOOLS/Phys/Blob.H"
#include "ATOOLS/Math/Histogram.H"
#include <map>
......@@ -10,15 +10,15 @@
namespace AHADIC {
class Gluon_Decayer {
private:
Dipole_Splitter * p_splitter;
double m_pt2max, m_pt2max_factor, m_pt02;
FlavDecayMap m_options;
DipoleList m_dipoles;
Gluon_Splitter * p_gsplitter;
double m_pt2max, m_pt2max_factor, m_pt02;
FlavDecayMap m_options;
DipoleList m_dipoles;
bool m_analyse;
std::map<std::string,ATOOLS::Histogram *> m_histograms;
void Init();
bool FillDipoleList(Proto_Particle_List * plin=NULL);
void UpdatePPList(Proto_Particle_List * plin=NULL);
void PrintDipoleList();
......@@ -33,7 +33,7 @@ namespace AHADIC {
long int m_tot, m_s, m_u, m_d, m_di0, m_di1;
public:
Gluon_Decayer(Dipole_Splitter *,bool=false);
Gluon_Decayer(bool ana=false);
~Gluon_Decayer();
bool DecayList(Proto_Particle_List * plin=NULL);
};
......
#include <cassert>
#include "AHADIC++/Main/Ahadic.H"
#include "AHADIC++/Tools/Soft_Cluster_Handler.H"
#include "AHADIC++/Tools/Dipole.H"
#include "AHADIC++/Tools/Cluster.H"
#include "ATOOLS/Org/Message.H"
#include "ATOOLS/Org/Run_Parameter.H"
......@@ -44,19 +45,13 @@ Ahadic::~Ahadic()
Return_Value::code Ahadic::Hadronize(Blob_List * blobs)
{
// msg_Out()<<"=======================================================\n"
// <<"=======================================================\n"
// <<"=======================================================\n"
// <<"=======================================================\n"
// <<"=======================================================\n";
static std::string mname(METHOD);
Return_Value::IncCall(mname);
if (msg->LevelIsDebugging()) {
msg_Out()<<"#######################################################\n"
<<"################### IN ################################\n";
for (Blob_List::iterator blit=blobs->begin();blit!=blobs->end();++blit) {
if ((*blit)->Has(blob_status::needs_hadronization) &&
(*blit)->Type()==btp::Fragmentation)
msg_Out()<<"#################################################\n"
<<(**blit)<<"\n"
<<"#################################################\n";
}
}
Blob * blob(NULL);
bool moveon(false);
double norm2(sqr(rpa->gen.Ecms()));
......@@ -65,35 +60,18 @@ Return_Value::code Ahadic::Hadronize(Blob_List * blobs)
if ((*blit)->Has(blob_status::needs_hadronization) &&
(*blit)->Type()==btp::Fragmentation) {
blob = (*blit);
//Particle * part;
//double DeltaR, DeltaY, DeltaPhi,S12;
//Vec4D mom1, mom2;
//for (int i=0;i<blob->NInP()-1;i++) {
// part = blob->InParticle(i);
// mom1 = blob->InParticle(i)->Momentum();
// mom2 = blob->InParticle(i+1)->Momentum();
// DeltaPhi = dabs(mom1.Phi()-mom2.Phi());
// DeltaY = (mom1.Eta()-mom2.Eta());
// DeltaR = sqrt(DeltaY*DeltaY+DeltaPhi*DeltaPhi);
// S12 = (mom1+mom2).Abs2();
//}
//msg_Out()<<" {"<<mom1.Y()<<", "<<mom1.Phi()<<", "<<mom1.PPerp()<<"} "
// <<" ["<<part->GetFlow(1)<<", "<<part->GetFlow(2)<<"] "
// <<" from "<<part->ProductionBlob()->Id()<<"\n"
// <<"-----------------------------------------------\n";
//part = blob->InParticle(blob->NInP()-1);
moveon = false;
Reset();
for (short int i=0;i<m_maxtrials;i++) {
try {
result = Hadronize(blob,i);
} catch (Return_Value::code ret) {
msg_Error()<<"ERROR in "<<METHOD<<" : "<<std::endl
<<" Hadronization for blob "
msg_Error()<<"ERROR in "<<METHOD<<" : \n"
<<" Hadronization for blob "<<(*blob)<<"\n"
<<"("<<blob<<"; "<<blob->NInP()<<" -> "
<<blob->NOutP()<<") "
<<"did not work out,"<<std::endl
<<" will trigger retrying the event."<<std::endl;
<<"did not work out,\n"
<<" will trigger retrying the event.\n";
CleanUp(blob);
return Return_Value::Retry_Event;
}
......@@ -106,16 +84,17 @@ Return_Value::code Ahadic::Hadronize(Blob_List * blobs)
case Return_Value::Retry_Event :
{
blobs->ColorConservation();
msg_Tracking()<<"ERROR in "<<METHOD<<" : "<<std::endl
msg_Tracking()<<"ERROR in "<<METHOD<<" :\n"
<<" Hadronization for blob "
<<"("<<blob<<"; "<<blob->NInP()<<" -> "
<<blob->NOutP()<<") "
<<"did not work out,"<<std::endl
<<" will trigger retrying the event."<<std::endl;
<<" will trigger retrying the event.\n";
CleanUp(blob);
if (rpa->gen.Beam1().IsLepton() ||
rpa->gen.Beam2().IsLepton()) {
msg_Tracking()<<METHOD<<"(): Non-hh collision. Request new event instead.\n";
msg_Tracking()<<METHOD<<": "
<<"Non-hh collision. Request new event instead.\n";
return Return_Value::New_Event;
}
return result;
......@@ -227,7 +206,8 @@ void Ahadic::Reset() {
}
bool Ahadic::SanityCheck(Blob * blob,double norm2) {
if (dabs(blob->CheckMomentumConservation().Abs2())/norm2>1.e-12 ||
Vec4D checkmom(blob->CheckMomentumConservation());
if (dabs(checkmom.Abs2())/norm2>1.e-12 ||
(norm2<0. && norm2>0.) ||
Cluster::RemainingClusters()!=0 ||
control::s_AHAparticles!=blob->NOutP()
......@@ -238,8 +218,9 @@ bool Ahadic::SanityCheck(Blob * blob,double norm2) {
<<" remaining clusters (norm2 = "<<norm2<<")."<<endl
<<" Protoparticles = "<<control::s_AHAprotoparticles
<<"/ parts = "<<control::s_AHAparticles<<" vs. "<<blob->NOutP()
<<" : "<<blob->CheckMomentumConservation()<<endl
<<" : "<<checkmom<<" ("<<sqrt(Max(0.,checkmom.Abs2()))<<")\n"
<<(*blob)<<endl;
exit(0);
return false;
}
msg_Tracking()<<"Passed "<<METHOD<<" with "
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
#ifndef AHADIC_Tools_Cluster_Splitter_H
#define AHADIC_Tools_Cluster_Splitter_H
#include "AHADIC++/Tools/Splitter_Base.H"
#include "AHADIC++/Tools/Cluster.H"
namespace AHADIC {
class Cluster_Splitter : public Splitter_Base {
private:
size_t m_nmax;
double m_etax, m_etax_lead, m_etay, m_etay_lead;
size_t m_pairs;