Commit 7c2b1d2b authored by Frank Krauss's avatar Frank Krauss Committed by Enrico Bothmann

Implement new AMISIC and remnannt handling

 - Functional for simple LEP- and LHC-type events
 - Minor bug-fixes in AHADIC
 - Minor updates to ATOOLS to allow for better conservation-law checks
   in blobs
 - Pull remnants out of the ISR_Handler in CSS and out of PDF into a new
   remnant handler, adding the logic of colour extraction
 - Move mass mode to allow AMISIC to have access to the Mass_Selector
 - AMISIC and REMNANTS are ready to be tuned
parent 1bf9800f
......@@ -33,8 +33,6 @@ bool Cluster_Splitter::MakeLongitudinalMomenta() {
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;
//msg_Out()<<METHOD<<"("<<p_part1->Flavour()<<"/"<<m_newflav1<<", "<<m_R12<<", "<<m_minQ_1<<") "
// <<"("<<m_newflav2<<"/"<<p_part2->Flavour()<<", "<<m_R21<<", "<<m_minQ_2<<")\n";
return CheckIfAllowed();
}
......
......@@ -41,6 +41,7 @@ bool Gluon_Decayer::operator()(Singlet * singlet) {
msg_Error()<<"Couldn't deal with 2-parton singlet.\n"<<(*singlet)<<"\n";
exit(1);
}
delete p_singlet;
return flag;
}
m_origsize=p_singlet->size();
......@@ -77,6 +78,8 @@ bool Gluon_Decayer::operator()(Singlet * singlet) {
else
p_singlet->pop_back();
case 0:
//msg_Out()<<METHOD<<": Step("<<part1<<", "<<part2<<") yields 0, now "
// <<p_singlet->size()<<" partons in singlet.\n";
break;
case -1:
default:
......
......@@ -55,8 +55,8 @@ Singlet_Checker::Singlet_Checker(list<Singlet *> * singlets,
{}
Singlet_Checker::~Singlet_Checker() {
msg_Out()<<METHOD<<" with "<<m_direct_transitions
<<" direct enforced transitions in total.\n";
msg_Tracking()<<METHOD<<" with "<<m_direct_transitions
<<" direct enforced transitions in total.\n";
}
void Singlet_Checker::Init() {
......@@ -98,6 +98,12 @@ bool Singlet_Checker::operator()() {
if (m_badones.size()>0) {
if (!DealWithProblematicSinglets()) {
msg_Error()<<METHOD<<" throws error - no rescue possible.\n";
if (msg_LevelIsTracking()) {
for (list<list<Singlet *>::iterator>::iterator bit=m_badones.begin();
bit!=m_badones.end();bit++) {
msg_Out()<<(***bit)<<"\n";
}
}
return false;
}
}
......@@ -289,8 +295,8 @@ bool Singlet_Checker::TransitProblematicSinglets() {
if (totmom.Abs2()<sqr(totmass)) {
for (map<Singlet *,Flavour>::iterator tit=m_transitions.begin();
tit!=m_transitions.end();tit++,i++) {
msg_Out()<<"Singlet with "<<tit->first->Momentum()<<" --> "
<<tit->second<<" ("<<tit->second.Mass()<<")\n";
msg_Debugging()<<"Singlet with "<<tit->first->Momentum()<<" --> "
<<tit->second<<" ("<<tit->second.Mass()<<")\n";
}
delete[] moms;
delete[] masses;
......
......@@ -28,21 +28,21 @@ namespace AHADIC {
long int m_direct_transitions;
bool CheckSinglet();
bool FusePartonsInLowMassSinglet();
bool DealWithProblematicSinglets();
void SortProblematicSinglets();
bool FindOtherSingletToTransit();
bool FindRecoilerForTransit();
bool TestRecoiler();
bool TransitProblematicSinglets();
bool TransitProblematicSingletWithRecoiler();
bool BoostRecoilerInNewSystem(const ATOOLS::Vec4D & newmom);
void ForcedDecays();
bool ForcedDecayOfTwoPartonSinglet();
bool ExtractAndCheckFlavours();
bool TwoGluonSingletToHadrons();
bool TwoQuarkSingletToHadrons();
bool CheckSinglet();
bool FusePartonsInLowMassSinglet();
bool DealWithProblematicSinglets();
void SortProblematicSinglets();
bool FindOtherSingletToTransit();
bool FindRecoilerForTransit();
bool TestRecoiler();
bool TransitProblematicSinglets();
bool TransitProblematicSingletWithRecoiler();
bool BoostRecoilerInNewSystem(const ATOOLS::Vec4D & newmom);
void ForcedDecays();
bool ForcedDecayOfTwoPartonSinglet();
bool ExtractAndCheckFlavours();
bool TwoGluonSingletToHadrons();
bool TwoQuarkSingletToHadrons();
public:
Singlet_Checker(std::list<Singlet *> * p_singlets,
Soft_Cluster_Handler * softclusters);
......
......@@ -48,7 +48,7 @@ Return_Value::code Ahadic::Hadronize(Blob_List * blobs)
switch (Hadronize(blob)) {
case Return_Value::Success :
break;
case Return_Value::Retry_Event :
case Return_Value::Retry_Event :
blobs->ColorConservation();
msg_Tracking()<<"ERROR in "<<METHOD<<" :\n"
<<" Hadronization for blob "
......@@ -78,13 +78,19 @@ Return_Value::code Ahadic::Hadronize(Blob_List * blobs)
}
Return_Value::code Ahadic::Hadronize(Blob * blob, int retry) {
//msg_Out()<<"######################################################################\n"
// <<(*blob)<<"\n";
// <<"######################################################################\n"
// <<"######################################################################\n";
Reset();
m_totmom = blob->CheckMomentumConservation();
if (!ExtractSinglets(blob) || !ShiftBeamParticles() || !CheckSinglets() ||
!DecayGluons() ||!DecayClusters()) {
msg_Error()<<"ERROR in "<<METHOD<<": Will retry event!\n"
<<(*blob);
//exit(1);
Reset(blob);
Reset();
return Return_Value::Retry_Event;
}
blob->UnsetStatus(blob_status::needs_hadronization);
......@@ -113,7 +119,6 @@ bool Ahadic::ExtractSinglets(Blob * blob)
bool Ahadic::ShiftBeamParticles()
{
//msg_Out()<<METHOD<<"\n";
if (!m_beamparticles()) {
msg_Error()<<METHOD<<" could not shift beam particles on mass shells.\n";
return false;
......@@ -123,7 +128,6 @@ bool Ahadic::ShiftBeamParticles()
bool Ahadic::CheckSinglets()
{
//msg_Out()<<METHOD<<"\n";
if (!m_singletchecker()) {
msg_Error()<<METHOD<<" singlets did not check out.\n";
return false;
......@@ -132,7 +136,6 @@ bool Ahadic::CheckSinglets()
}
bool Ahadic::DecayGluons() {
//msg_Out()<<METHOD<<"\n";
while (!m_singlet_list.empty()) {
if (m_gluondecayer(m_singlet_list.front()))
m_singlet_list.pop_front();
......@@ -145,14 +148,12 @@ bool Ahadic::DecayGluons() {
}
bool Ahadic::DecayClusters() {
//msg_Out()<<METHOD<<":\n"<<m_cluster_list<<"\n";
bool success = m_clusterdecayer();
if (!success) msg_Error()<<METHOD<<" could not decay all clusters.\n";
return success;
}
void Ahadic::FillOutgoingParticles(Blob * blob) {
//msg_Out()<<METHOD<<"\n";
while (!m_hadron_list.empty()) {
Particle * part = (*m_hadron_list.front())();
part->SetNumber();
......@@ -197,5 +198,6 @@ bool Ahadic::SanityCheck(Blob * blob,double norm2) {
}
void Ahadic::CleanUp(Blob * blob) {
Reset();
if(blob) blob->DeleteOutParticles(0);
}
......@@ -834,11 +834,11 @@ void Hadron_Init::Init() {
"thomas","thomas");
s_kftable[5508]=new Particle_Info(5508,1000000.0,1000,0,0,0,0,
"tanju","tanju");
s_kftable[5508]=new Particle_Info(5509,1000000.0,1000,0,0,0,0,
s_kftable[5509]=new Particle_Info(5509,1000000.0,1000,0,0,0,0,
"jennifer","jennifer");
s_kftable[5508]=new Particle_Info(6505,1000000.0,1000,0,0,0,0,
s_kftable[6505]=new Particle_Info(6505,1000000.0,1000,0,0,0,0,
"hendrik","hendrik");
s_kftable[5508]=new Particle_Info(6506,1000000.0,1000,0,0,0,0,
s_kftable[6506]=new Particle_Info(6506,1000000.0,1000,0,0,0,0,
"jan","jan");
// ##########################################################################
// ##########################################################################
......
......@@ -6,10 +6,10 @@ using namespace AHADIC;
using namespace ATOOLS;
using namespace std;
Singlet::Singlet() { }
Singlet::Singlet() {}
Singlet::~Singlet() {
while (size()>0) {
while (!empty()) {
delete front();
pop_front();
}
......@@ -89,7 +89,7 @@ std::ostream& AHADIC::operator<<(std::ostream & str,const Singlet & sing) {
<<(*pliter)->Flavour()
<<" (lead = "<<(*pliter)->IsLeading()<<", "
<<"beam = "<<(*pliter)->IsBeam()<<")"
<<": "<<(*pliter)->Momentum()<<"\n";
<<": "<<(*pliter)->Momentum()<<"("<<(*pliter)->Momentum().Abs2()<<")\n";
str<<"***********************************************************\n";
return str;
}
......
This diff is collapsed.
#ifndef AMISIC_Main_Amisic_H
#define AMISIC_Main_Amisic_H
#include "AMISIC++/Perturbative/Single_Collision_Handler.H"
#include "AMISIC++/Perturbative/MI_Processes.H"
#include "AMISIC++/Tools/Over_Estimator.H"
#include "AMISIC++/Tools/Impact_Parameter.H"
#include "MODEL/Main/Model_Base.H"
#include "BEAM/Main/Beam_Spectra_Handler.H"
#include "PDF/Main/ISR_Handler.H"
#include "ATOOLS/Org/File_IO_Base.H"
#include "ATOOLS/Org/Default_Reader.H"
#include "ATOOLS/Math/Histogram.H"
#include <map>
#include <string>
/*!
\file Amisic.H
\brief Declares the class Amisic
*/
#include "AMISIC++/Main/MI_Base.H"
#include "MODEL/Main/Model_Base.H"
#include "BEAM/Main/Beam_Spectra_Handler.H"
#include "PDF/Main/ISR_Handler.H"
namespace ATOOLS {
class Cluster_Amplitude;
}
namespace AMISIC {
/*!
\namespace AMISIC
\brief Houses all classes dedicated to the generation of
the underlying event.
\brief The module for the generation of a perturbative underlying
event according to the Sjostrand-van der Zijl model.
AMISIC++ is an accronym for <em>A</em> <em>M</em>ultiple
<em>I</em>nteraction <em>S</em>imulation <em>I</em>n
<em>C++</em>. Accordingly the AMISIC namespace houses all
classes which are concerned with the generation of underlying
events.
*/
class Amisic: public ATOOLS::File_IO_Base {
private:
std::string m_hardmodel, m_softmodel;
MI_Base *p_hardbase, *p_softbase;
<em>C++</em>.
MODEL::Model_Base *p_model;
BEAM::Beam_Spectra_Handler *p_beam;
PDF::ISR_Handler *p_isr;
The AMISIC namespace contains the module for the generation of a
perturbative underlying event according to the Sjostrand-van der Zijl
model. The original model is based on perturbative QCD 2->2 scatters
modified through a simple IR regularisation and convoluted with regular
PDFs and and a matter-overlap function between the incident hadronic
states which gives rise to an interaction probability.
bool m_external;
In SHERPA we added/plan to add other 2->2 scatters, for example for photon
and quarkonia production.
*/
class Amisic: public ATOOLS::File_IO_Base, ATOOLS::Mass_Selector {
private:
double m_b, m_bfac, m_residualE1, m_residualE2, m_pt2;
public:
MI_Processes * p_processes;
Over_Estimator m_overestimator;
Impact_Parameter m_impact;
Single_Collision_Handler m_singlecollision;
// constructors
bool m_ana;
size_t m_nscatters;
std::map<std::string,ATOOLS::Histogram *> m_histos;
bool InitParameters(ATOOLS::Default_Reader *const defaultreader);
void CreateAmplitudeLegs(ATOOLS::Cluster_Amplitude * ampl,
ATOOLS::Blob * blob);
void FillAmplitudeSettings(ATOOLS::Cluster_Amplitude * ampl);
void InitAnalysis();
void FinishAnalysis();
void Analyse(const bool & last);
public:
Amisic();
/*!
\fn Amisic()
\brief The default constructor.
*/
Amisic(MODEL::Model_Base *const model,
BEAM::Beam_Spectra_Handler *const beam,
PDF::ISR_Handler *const isr);
/*!
\fn Amisic(MODEL::Model_Base *const model,
BEAM::Beam_Spectra_Handler *const beam,
PDF::ISR_Handler *const isr)
\brief The standard constructor.
Creates the underlying event handler employing
external ISR, beam spectra and model handlers.
*/
// destructor
~Amisic();
// member functions
bool Initialize();
/*!
\fn bool Initialize()
\brief Initializes the underlying event handling.
*/
void Reset();
/*!
\fn void Reset()
\brief Resets the underlying event handler.
This method is to be called before any event to be generated.
*/
void CleanUp();
/*!
\fn void CleanUp()
\brief Cleans the underlying event handler.
This method is to be called after any event generated.
*/
bool GenerateHardProcess(ATOOLS::Blob *blob);
/*!
\fn bool GenerateHardProcess(ATOOLS::Blob *blob)
\brief Generates one hard interaction, which is part of
the hard underlying event.
*/
bool GenerateHardEvent(ATOOLS::Blob_List *blobs);
/*!
\fn bool GenerateHardEvent(ATOOLS::Blob_List *blobs)
\brief Generates a multitude of hard interactions
which taken together form the hard underlying event.
*/
void SameHardProcess(ATOOLS::Blob *blob);
/*!
\fn void SameHardProcess(ATOOLS::Blob *blob)
\brief Returns last hard intreraction.
*/
bool GenerateSoftProcess(ATOOLS::Blob *blob);
/*!
\fn bool GenerateSoftProcess(ATOOLS::Blob *blob)
\brief Generates one soft interaction, which is part of
the soft underlying event.
*/
bool GenerateSoftEvent(ATOOLS::Blob_List *blobs);
/*!
\fn bool GenerateSoftEvent(ATOOLS::Blob_List *blobs)
\brief Generates a multitude of soft interactions
which taken together form the soft underlying event.
*/
void SameSoftProcess(ATOOLS::Blob *blob);
/*!
\fn void SameSoftProcess(ATOOLS::Blob *blob)
\brief Returns last soft intreraction.
*/
bool Initialize(ATOOLS::Default_Reader *const defaultreader,
MODEL::Model_Base *const model,
PDF::ISR_Handler *const isr);
bool VetoEvent(const double & scale=-1.);
ATOOLS::Blob * GenerateScatter();
bool VetoScatter(ATOOLS::Blob * blob);
int ShiftMasses(ATOOLS::Cluster_Amplitude * ampl);
ATOOLS::Cluster_Amplitude * ClusterConfiguration(ATOOLS::Blob * blob);
void SetMaxEnergies(const double & E1,const double & E2);
void SetMaxScale(const double & scale);
void SetB(const double & b=-1.);
const double ScaleMin() const;
const double ScaleMax() const;
inline void SetMassMode(const int & mode) {
p_processes->SetMassMode(mode);
}
inline double Mass(const ATOOLS::Flavour &fl) const {
return p_processes->Mass(fl);
}
bool VetoHardProcess(ATOOLS::Blob *blob);
/*!
\fn bool VetoHardProcess(ATOOLS::Blob *blob)
\brief Vetoes the hard interaction according to the
probability density employed in the hard underlying event model.
*/
bool GenerateEvent(ATOOLS::Blob_List *blobs);
/*!
\fn bool GenerateEvent(ATOOLS::Blob_List *blobs)
\brief Generates one underlying event, using the active
hard and soft underlying event models.
*/
bool SelectHardModel(const std::string &hardmodel);
/*!
\fn bool SelectHardModel(const std::string &hardmodel)
\brief Selects the hard underlying event model
according to the given tag.
*/
bool SelectSoftModel(const std::string &softmodel);
/*!
\fn bool SelectSoftModel(const std::string &softmodel)
\brief Selects the soft underlying event model
according to the given tag.
*/
// inline functions
inline const std::string &HardModel() const { return m_hardmodel; }
inline const std::string &SoftModel() const { return m_softmodel; }
inline MI_Base *HardBase() const { return p_hardbase; }
inline MI_Base *SoftBase() const { return p_softbase; }
inline PHASIC::Process_Base *HardXS() { return p_hardbase->XS(); }
inline PHASIC::Process_Base *SoftXS() { return p_softbase->XS(); }
void Reset();
void CleanUp();
}; // end of class Amisic
/*!
\class Amisic
\brief The overall steering of the AMISIC package
This class implements the overall steering of the
underlying event generation. It is interfaced by the Sherpa
framework but it may also act as a stand-alone underlying event
generator. It updates the MI_Base instances with the information
which is necessary for the next event generation and delivers
the created underlying event as a blob to the outside world.
*/
void Test();
};
} // end of namespace AMISIC
}
#endif
#include "AMISIC++/Main/MI_Base.H"
#include "ATOOLS/Phys/Particle.H"
#include "ATOOLS/Org/Exception.H"
using namespace AMISIC;
MI_Base::String_MI_Base_Map *MI_Base::s_bases;
bool MI_Base::s_stophard=true;
bool MI_Base::s_stopsoft=true;
bool MI_Base::s_cleaned=true;
MI_Base *MI_Base::s_hard=NULL;
MI_Base *MI_Base::s_soft=NULL;
MI_Base::MI_Base(std::string name,TypeID type,unsigned int nparameter,
unsigned int infiles,unsigned int outfiles):
File_IO_Base(infiles,outfiles),
m_name(name),
m_type(type),
m_start(NULL), m_stop(NULL), m_last(NULL),
m_nparameter(nparameter),
p_xs(NULL)
{
static bool initialized=false;
if (!initialized) {
s_bases = new String_MI_Base_Map();
initialized=true;
}
for (String_MI_Base_Map::iterator nbit=s_bases->begin();
nbit!=s_bases->end();++nbit) {
if (nbit->first==m_name) {
THROW(fatal_error,"MI_Base already exists!");
}
}
if (m_type==Unknown) {
THROW(fatal_error,"MI base has no type!");
}
if (m_nparameter>0) {
m_start = new double[m_nparameter];
m_stop = new double[m_nparameter];
m_last = new double[m_nparameter];
}
(*s_bases)[m_name]=this;
switch (m_type) {
case SoftEvent:
if (m_name!=TypeToString(type)+" None") s_soft=this;
break;
case HardEvent:
if (m_name!=TypeToString(type)+" None") s_hard=this;
break;
default:
break;
}
}
MI_Base::~MI_Base()
{
for (String_MI_Base_Map::iterator nbit=s_bases->begin();
nbit!=s_bases->end();++nbit) {
if (nbit->first==m_name) {
s_bases->erase(nbit);
break;
}
}
if (m_nparameter>0) {
delete [] m_start;
delete [] m_stop;
delete [] m_last;
}
}
void MI_Base::UpdateAll(const MI_Base *mibase)
{
for (String_MI_Base_Map::iterator nbit=s_bases->begin();
nbit!=s_bases->end();++nbit) {
nbit->second->Update(mibase);
}
}
void MI_Base::Update(const MI_Base *mibase)
{
msg_Error()<<"MI_Base::Update("<<mibase<<"): "
<<"Virtual method called!"<<std::endl;
return;
}
bool MI_Base::Initialize()
{
msg_Error()<<"MI_Base::Initialize(): "
<<"Virtual method called!"<<std::endl;
return false;
}
void MI_Base::Reset()
{
msg_Error()<<"MI_Base::Reset(): "
<<"Virtual method called!"<<std::endl;
return;
}
void MI_Base::CleanUp()
{
for (String_MI_Base_Map::iterator nbit=s_bases->begin();
nbit!=s_bases->end();++nbit) {
nbit->second->m_inparticles.Clear();
nbit->second->m_outparticles.Clear();
}
s_stophard=false;
s_stopsoft=false;
s_cleaned=true;
}
bool MI_Base::VetoProcess(ATOOLS::Blob *blob)
{
msg_Error()<<"MI_Base::VetoProcess(): "
<<"Virtual method called!"<<std::endl;
return false;
}
bool MI_Base::GenerateProcess()
{
msg_Error()<<"MI_Base::GenerateProcess(): "
<<"Virtual method called!"<<std::endl;
return false;
}
void MI_Base::ResetAll()
{
for (String_MI_Base_Map::iterator nbit=s_bases->begin();
nbit!=s_bases->end();++nbit) {
nbit->second->Reset();
}
}
bool MI_Base::FillBlob(ATOOLS::Blob *blob)
{
if (blob==NULL) {
msg_Error()<<"MI_Base::FillBlob(..): "
<<"Blob is not initialized!"<<std::endl
<<" Cannot proceed in filling."<<std::endl;
return false;
}
if (!m_generatedprocess) return false;
if (m_inparticles.empty()) {
msg_Error()<<"MI_Base::FillBlob(..): "
<<"Did not create any process yet!"<<std::endl
<<" Cannot proceed in filling."<<std::endl;
return false;
}
bool generatedprocess=m_generatedprocess;
m_generatedprocess=false;
if (m_type==HardEvent) {
blob->SetType(ATOOLS::btp::Hard_Collision);
blob->SetStatus(ATOOLS::blob_status::needs_showers &
ATOOLS::blob_status::needs_beams &
ATOOLS::blob_status::needs_hadronization);
}
else {
blob->SetType(ATOOLS::btp::Soft_Collision);
blob->SetStatus(ATOOLS::blob_status::needs_beams &
ATOOLS::blob_status::needs_hadronization);
}
ATOOLS::Particle *particle;
for (size_t i=0;i<m_inparticles.size();++i) {
particle = new ATOOLS::Particle(-1,m_inparticles[i]->Flav(),
m_inparticles[i]->Momentum());
particle->SetFlow(1,m_inparticles[i]->GetFlow(1));
particle->SetFlow(2,m_inparticles[i]->GetFlow(2));
particle->SetNumber(1);
particle->SetStatus(ATOOLS::part_status::active);
particle->SetInfo('G');
blob->AddToInParticles(particle);
}
for (size_t i=0;i<m_outparticles.size();++i) {
particle = new ATOOLS::Particle(-1,m_outparticles[i]->Flav(),
m_outparticles[i]->Momentum());
particle->SetFlow(1,m_outparticles[i]->GetFlow(1));
particle->SetFlow(2,m_outparticles[i]->GetFlow(2));
particle->SetNumber(1);
particle->SetStatus(ATOOLS::part_status::active);
particle->SetInfo('H');
blob->AddToOutParticles(particle);
}
return generatedprocess;
}
std::string MI_Base::TypeToString(TypeID type)
{
switch (type) {
case HardEvent: return std::string("Hard Event");
case SoftEvent: return std::string("Soft Event");
case Unknown : return std::string("Unknown");
}
return std::string("Unknown");
}
MI_Base::TypeID MI_Base::StringToType(std::string type)
{
if (type==std::string("Hard Event")) return HardEvent;