#include"SHERPA/PerturbativePhysics/Hard_Decay_Handler.H" #include "ATOOLS/Org/Message.H" #include "ATOOLS/Org/Data_Reader.H" #include "ATOOLS/Org/Return_Value.H" #include "ATOOLS/Org/Shell_Tools.H" #include "ATOOLS/Phys/Blob_List.H" #include "ATOOLS/Phys/Cluster_Amplitude.H" #include "ATOOLS/Phys/NLO_Subevt.H" #include "ATOOLS/Phys/Weight_Info.H" #include "ATOOLS/Math/Random.H" #include "PHASIC++/Decays/Decay_Map.H" #include "PHASIC++/Decays/Decay_Table.H" #include "PHASIC++/Decays/Decay_Channel.H" #include "MODEL/Main/Model_Base.H" #include "MODEL/Main/Single_Vertex.H" #include "MODEL/Main/Color_Function.H" #include "PHASIC++/Channels/Multi_Channel.H" #include "PHASIC++/Channels/Rambo.H" #include "PHASIC++/Channels/Decay_Dalitz.H" #include "METOOLS/Main/Spin_Structure.H" #include "METOOLS/SpinCorrelations/Polarized_CrossSections_Handler.H" #include "ATOOLS/Org/Run_Parameter.H" #include "SHERPA/SoftPhysics/Soft_Photon_Handler.H" #include "ATOOLS/Phys/Variations.H" #include "ATOOLS/Phys/KF_Table.H" #include "EXTRA_XS/One2Two/Comix1to2.H" #include "EXTRA_XS/One2Three/Comix1to3.H" #include #include #include #include #include #include using namespace SHERPA; using namespace ATOOLS; using namespace MODEL; using namespace EXTRAXS; using namespace PHASIC; using namespace METOOLS; using namespace std; class ParticlePairFirstEnergySort { public: bool operator()(const ParticlePair& a,const ParticlePair& b) { return (a.first->Momentum()[0]Momentum()[0]); } }; class ParticlePairPairFirstEnergySort { public: bool operator()(const ParticlePairPair& a,const ParticlePairPair& b) { return (a.first.first->Momentum()[0]+a.first.second->Momentum()[0] Momentum()[0]+b.first.second->Momentum()[0]); } }; Hard_Decay_Handler::Hard_Decay_Handler() : p_newsublist(NULL), m_resultdir(""), m_offshell(""), m_decay_tau(false), m_set_widths(false), m_br_weights(true), m_usemass(true), m_min_prop_width(0.0) { auto& s = Settings::GetMainSettings(); auto ds = s["HARD_DECAYS"]; m_mass_smearing = ds["Mass_Smearing"].SetDefault(1).Get(); /* TODO: Writing out a 1->3 channel which might have a large width in one resolved configuration, and a small one in another? */ m_store_results = ds["Store_Results"].SetDefault(0).Get(); m_br_weights = ds["Apply_Branching_Ratios"].SetDefault(true).Get(); m_decay_tau = ds["Decay_Tau"].SetDefault(false).Get(); m_set_widths = ds["Set_Widths"].SetDefault(false).Get(); m_min_prop_width = ds["Min_Prop_Width"].SetDefault(0.0).Get(); m_int_accuracy = ds["Int_Accuracy"].SetDefault(0.01).Get(); m_int_niter = ds["Int_NIter"].SetDefault(2500).Get(); m_int_target_mode = ds["Int_Target_Mode"].SetDefault(0).Get(); m_offshell = ds["Resolve_Decays"] .SetDefault("Threshold") .UseNoneReplacements() .Get(); m_resultdir = ds["Result_Directory"] .SetDefault(s["RESULT_DIRECTORY"].Get() + "/Decays") .Get(); if (m_store_results) MakeDir(m_resultdir, true); m_spincorr=rpa->gen.HardSC(); m_polcrosssec = ds["Pol_Cross_Section"]["Enabled"].SetDefault(false).Get(); // if also final state polarization should be enabled, initialize polarization handler in a separate method in the // Initialization_Handler analogous to the Hard_Decay_Handler if (m_polcrosssec){ if (!m_spincorr){ THROW(fatal_error, "Calculation of polarized cross sections only possible together with spin correlations") } p_polarization_handler = new METOOLS::Polarized_CrossSections_Handler(); } // also need to tell shower whats massive now // TODO: need to use the same mass-selector // for now, implement alike SetDecayMasses(); m_qedmode = ds["QED_Corrections"].SetDefault(1).Get(); if (m_qedmode && m_spincorr) m_qedmode=2; if (m_qedmode>0 && !m_usemass) { THROW(fatal_error,std::string("QED corrections to hard decays only ") +std::string("available in massive mode.")); } DEBUG_FUNC(""); p_decaymap = new Decay_Map(this); KFCode_ParticleInfo_Map::const_iterator it; for (it=s_kftable.begin();it!=s_kftable.end();it++) { Flavour flav(it->first); if (Decays(flav)) { Decay_Table* dt=new Decay_Table(flav, this); p_decaymap->insert(make_pair(flav,dt)); ReadDecayTable(flav); } if (flav!=flav.Bar() && Decays(flav.Bar())) { Decay_Table* dt=new Decay_Table(flav.Bar(), this); p_decaymap->insert(make_pair(flav.Bar(),dt)); ReadDecayTable(flav.Bar()); } } // initialize them sorted by masses: Decay_Map::iterator dmit; msg_Debugging()<<"Initialising hard decay tables: two-body decays.\n"; for (dmit=p_decaymap->begin(); dmit!=p_decaymap->end(); ++dmit) { InitializeDirectDecays(dmit->second); } msg_Debugging()<<"Initialising hard decay tables: three-body decays.\n"; for (dmit=p_decaymap->begin(); dmit!=p_decaymap->end(); ++dmit) { InitializeOffshellDecays(dmit->second); dmit->second->UpdateChannelStatuses(); } // overwrite partial widths from HO SM or user input if (ds["Use_HO_SM_Widths"].SetDefault(true).Get()) { if (Flavour(kf_h0).Mass()<125.07 || Flavour(kf_h0).Mass()>125.11 || Flavour(kf_t).Mass()<172.4 || Flavour(kf_t).Mass()>172.6) THROW(fatal_error, "Use_HO_SM_Widths specified, but particle masses not in SM range."); SetHOSMWidths(ds); } for (dmit=p_decaymap->begin(); dmit!=p_decaymap->end(); ++dmit) { for (auto dc: *(dmit->second)) { auto s = ds["Channels"][dc->IDCode()]["Width"]; if (!s.HasDefault()) s.SetDefault(dc->IWidth()); dc->SetWidth(s.Get()); } dmit->second->UpdateWidth(); } if (m_set_widths) for (dmit=p_decaymap->begin(); dmit!=p_decaymap->end(); ++dmit) { dmit->second->Flav().SetWidth(dmit->second->TotalWidth()); } if (p_decaymap->size()) msg_Info()<size(); ++i) delete (*p_newsublist)[i]; p_newsublist->clear(); delete p_newsublist; } } void Hard_Decay_Handler::SetDecayMasses() { Settings& s = Settings::GetMainSettings(); ATOOLS::Flavour_Vector allflavs(MODEL::s_model->IncludedFlavours()); std::vector defpsmassive,defpsmassless; const std::vector psmassive { s["MASSIVE_PS"].GetVector() }; const std::vector psmassless { s["MASSLESS_PS"].GetVector() }; const bool respect{ s["RESPECT_MASSIVE_FLAG"].Get() }; // check consistency for (size_t i(0);i def: dusgy massless, rest massive // respect=1 -> def: only massive massive, rest massless // TODO: need to fill in those that are massive already? if (!respect) { defpsmassless.push_back(kf_d); defpsmassless.push_back(kf_u); defpsmassless.push_back(kf_s); defpsmassless.push_back(kf_gluon); defpsmassless.push_back(kf_photon); for (size_t i(0);iMass(true)!=fit->Mass(false)) mf.push_back(*fit); msg_Info()<Flav()); Flavour inflav=dt->Flav(); Vertex_Table::const_iterator vlit=s_model->VertexTable()->find(inflav); const Vertex_List& vlist=(vlit!=s_model->VertexTable()->end()) ? (vlit->second) : Vertex_List(); // temporary hack: // prepare reduced vertex list, not containing duplicates // in the sense of having identical external flavours Vertex_List reduced_vlist; for (Vertex_List::const_iterator vit=vlist.begin();vit!=vlist.end();++vit){ Vertex_List::const_iterator vjt=reduced_vlist.begin(); for(;vjt!=reduced_vlist.end(); ++vjt) if((*vjt)->in == (*vit)->in) break; if(vjt==reduced_vlist.end()) reduced_vlist.push_back(*vit); } msg_Debugging()<<"Vertices:"<NLegs(); ++j) dc->AddDecayProduct(sv->in[j]); Comix1to2* diagram=new Comix1to2(dc->Flavs()); dc->AddDiagram(diagram); dc->SetChannels(new Multi_Channel("")); dc->Channels()->SetNin(1); dc->Channels()->SetNout(dc->NOut()); Rambo* rambo = new Rambo(1,dc->NOut(),&dc->Flavs().front(),this); dc->Channels()->Add(rambo); dc->Channels()->Reset(); auto s = Settings::GetMainSettings()["HARD_DECAYS"]["Channels"][dc->IDCode()]; dc->SetActive(s["Status"].SetDefault(dc->Active()).GetVector()); if (CalculateWidth(dc)) dt->AddDecayChannel(dc); else delete dc; } dt->UpdateWidth(); if (m_set_widths) dt->Flav().SetWidth(dt->TotalWidth()); } void Hard_Decay_Handler::InitializeOffshellDecays(Decay_Table* dt) { DEBUG_FUNC(dt->Flav()<<" "<size()); size_t dtsize=dt->size(); for (size_t i=0;iat(i); vector new_dcs=ResolveDecay(dc); if (TriggerOffshell(dc, new_dcs)) { dc->SetActiveAll(-1); for (size_t j=0; jGetDecayChannel(new_dcs[j]->Flavs()); if (dup && dup->Active(0)>=0) { DEBUG_INFO("Adding new diagram to "<<*dup); for (size_t k=0; kGetDiagrams().size(); ++k) { dup->AddDiagram(new_dcs[j]->GetDiagrams()[k]); } for (size_t k=0; kChannels()->Channels().size(); ++k) { dup->AddChannel(new_dcs[j]->Channels()->Channels()[k]); } new_dcs[j]->ResetDiagrams(); new_dcs[j]->ResetChannels(); delete new_dcs[j]; new_dcs[j]=NULL; CalculateWidth(dup); } else { DEBUG_INFO("Adding "<Name()); auto s = Settings::GetMainSettings()["HARD_DECAYS"]["Channels"][new_dcs[j]->IDCode()]; new_dcs[j]->SetActive(s["Status"].SetDefault(new_dcs[j]->Active()).GetVector()); dt->AddDecayChannel(new_dcs[j]); } } } else { DEBUG_INFO("Keeping factorised."); for (size_t j=0; jSetActiveAll(-1); } } } dt->UpdateWidth(); if (m_set_widths) dt->Flav().SetWidth(dt->TotalWidth()); } void Hard_Decay_Handler::SetHOSMWidths(ATOOLS::Scoped_Settings& s) { // Higgs WG BRs 2022 s["Channels"]["25,5,-5"] ["Width"].SetDefault(2.382E-03); s["Channels"]["25,15,-15"] ["Width"].SetDefault(2.565E-04); s["Channels"]["25,13,-13"] ["Width"].SetDefault(8.901E-07); s["Channels"]["25,4,-4"] ["Width"].SetDefault(1.182E-04); s["Channels"]["25,3,-3"] ["Width"].SetDefault(1E-06); s["Channels"]["25,21,21"] ["Width"].SetDefault(3.354E-04); s["Channels"]["25,22,22"] ["Width"].SetDefault(9.307E-06); s["Channels"]["25,23,22"] ["Width"].SetDefault(6.318E-06); s["Channels"]["24,2,-1"] ["Width"].SetDefault(0.7041); s["Channels"]["24,4,-3"] ["Width"].SetDefault(0.7041); s["Channels"]["24,12,-11"] ["Width"].SetDefault(0.2256); s["Channels"]["24,14,-13"] ["Width"].SetDefault(0.2256); s["Channels"]["24,16,-15"] ["Width"].SetDefault(0.2256); s["Channels"]["-24,-2,1"] ["Width"].SetDefault(0.7041); s["Channels"]["-24,-4,3"] ["Width"].SetDefault(0.7041); s["Channels"]["-24,-12,11"]["Width"].SetDefault(0.2256); s["Channels"]["-24,-14,13"]["Width"].SetDefault(0.2256); s["Channels"]["-24,-16,15"]["Width"].SetDefault(0.2256); s["Channels"]["23,1,-1"] ["Width"].SetDefault(0.3828); s["Channels"]["23,2,-2"] ["Width"].SetDefault(0.2980); s["Channels"]["23,3,-3"] ["Width"].SetDefault(0.3828); s["Channels"]["23,4,-4"] ["Width"].SetDefault(0.2980); s["Channels"]["23,5,-5"] ["Width"].SetDefault(0.3828); s["Channels"]["23,11,-11"] ["Width"].SetDefault(0.0840); s["Channels"]["23,12,-12"] ["Width"].SetDefault(0.1663); s["Channels"]["23,13,-13"] ["Width"].SetDefault(0.0840); s["Channels"]["23,14,-14"] ["Width"].SetDefault(0.1663); s["Channels"]["23,15,-15"] ["Width"].SetDefault(0.0840); s["Channels"]["23,16,-16"] ["Width"].SetDefault(0.1663); s["Channels"]["6,24,5"] ["Width"].SetDefault(1.32); s["Channels"]["-6,-24,-5"] ["Width"].SetDefault(1.32); } bool Hard_Decay_Handler::TriggerOffshell(Decay_Channel* dc, vector new_dcs) { DEBUG_FUNC(dc->Name()<<"... "<Flavs().size(); ++j) outmass+=this->Mass(dc->Flavs()[j]); DEBUG_INFO(this->Mass(dc->Flavs()[0])<<" vs "<Mass(dc->Flavs()[0])Width(); } return (sum_resolved_widths>dc->Width()); } else if (m_offshell=="None") { return false; } else { THROW(fatal_error, "Parameter HARD_DECAYS:Resolve_Decays set to wrong value."); } return false; } vector Hard_Decay_Handler::ResolveDecay(Decay_Channel* dc1) { DEBUG_FUNC(dc1->Name()); vector new_dcs; const std::vector flavs1(dc1->Flavs()); for (size_t j=1;j Z Z? // Further iterations: W+ -> b t -> b W b -> b b .. ? if (flavs1[j]==flavs1[k]) ignore=true; } if (ignore) continue; Vertex_Table::const_iterator it=s_model->VertexTable()->find(flavs1[j]); const Vertex_List& vlist(it->second); // temporary hack: // prepare reduced vertex list, not containing duplicates // in the sense of having identical external flavours Vertex_List reduced_vlist; for (Vertex_List::const_iterator vit=vlist.begin();vit!=vlist.end();++vit){ Vertex_List::const_iterator vjt=reduced_vlist.begin(); for(;vjt!=reduced_vlist.end(); ++vjt) if((*vjt)->in == (*vit)->in) break; if(vjt==reduced_vlist.end()) reduced_vlist.push_back(*vit); } for (size_t k=0;k3 only Decay_Channel* dc=new Decay_Channel(flavs1[0], this); size_t nonprop(0), propi(0), propj(0); dc->AddDecayProduct(flavs1[3-j]); dc->AddDecayProduct(sv->in[1]); dc->AddDecayProduct(sv->in[2]); DEBUG_FUNC("trying "<Name()); // TODO what about W+ -> b t -> b W b -> b b ... two diagrams, factor 2, ...? // TODO what about identical particles like W' -> b t -> b W b ... two diagrams, factor 2, ...? // TODO what about W -> W gamma -> l v gamma? for (size_t l=1; l<4; ++l) { if (dc->Flavs()[l]==flavs1[3-j]) nonprop=l; } for (size_t l=1; l<4; ++l) { if (l!=nonprop && propi>0) propj=l; else if (l!=nonprop && propi==0) propi=l; } assert(dc1->GetDiagrams().size()==1); DEBUG_VAR(dc->Flavs()); DEBUG_VAR(flavs1[j]); Comix1to3* diagram=new Comix1to3(dc->Flavs(),flavs1[j], nonprop, propi, propj); dc->AddDiagram(diagram); dc->SetChannels(new Multi_Channel("")); dc->Channels()->SetNin(1); dc->Channels()->SetNout(dc->NOut()); dc->Channels()->Add(new Rambo(1,dc->NOut(),&dc->Flavs().front(),this)); if (flavs1[j].Width()>0.0) { dc->Channels()->Add(new Decay_Dalitz(&dc->Flavs().front(), flavs1[j].Mass(), flavs1[j].Width(), nonprop, propi, propj, this)); } dc->Channels()->Reset(); if (CalculateWidth(dc)) new_dcs.push_back(dc); else delete dc; } } return new_dcs; /* TODO: What about cases where a 1->3 was resolved from one 1->2, but would have been possible from another 1->2 where it was decided not to be resolved? */ } bool Hard_Decay_Handler::CalculateWidth(Decay_Channel* dc) { // Integrate or use results read from decay table file double outmass=0.0; for (size_t l=1; lFlavs().size(); ++l) outmass+=this->Mass(dc->Flavs()[l]); if (this->Mass(dc->Flavs()[0])>outmass) { DEBUG_FUNC("Starting calculation of "<Name()<<" width now."); if (m_store_results && m_read.find(dc->Flavs()[0])!=m_read.end()) { if (m_read[dc->Flavs()[0]].find(dc->IDCode())!= m_read[dc->Flavs()[0]].end()) { const vector& results(m_read[dc->Flavs()[0]][dc->IDCode()]); dc->SetIWidth(results[0]); dc->SetIDeltaWidth(results[1]); dc->SetMax(results[2]); } else { msg_Tracking()<<" Integrating "<Name()<CalculateWidth(m_int_accuracy, m_int_target_mode==0 ? dc->Flavs()[0].Width() : 0.0, m_int_niter); } } else { msg_Tracking()<<" Integrating "<Name()<CalculateWidth(m_int_accuracy, m_int_target_mode==0 ? dc->Flavs()[0].Width() : 0.0, m_int_niter); } } else { dc->SetActiveAll(-1); dc->SetIWidth(0.0); dc->SetIDeltaWidth(0.0); dc->SetMax(0.0); } auto s = Settings::GetMainSettings()["HARD_DECAYS"]["Channels"][dc->IDCode()]; dc->SetWidth(dc->IWidth()); dc->SetDeltaWidth(dc->IDeltaWidth()); return true; } bool Hard_Decay_Handler::ProperVertex(MODEL::Single_Vertex* sv) { if (sv->dec) return false; for (int i(0); iNLegs(); ++i) if (sv->in[i].IsDummy()) return false; if (sv->NLegs()!=3) return false; // TODO // TODO: ignore radiation graphs. should we? for (int i=1; iNLegs(); ++i) { if (sv->in[i].Kfcode()==sv->in[0].Kfcode()) { return false; } } // what about extra particles like Z4 if Z stable? return true; } void Hard_Decay_Handler::CreateDecayBlob(ATOOLS::Particle* inpart) { DEBUG_FUNC(inpart->Flav()); if(inpart->DecayBlob()) THROW(fatal_error,"Decay blob already exists."); if(!Decays(inpart->Flav())) return; Blob* blob = p_bloblist->AddBlob(btp::Hard_Decay); blob->SetStatus(blob_status::needs_showers); blob->AddStatus(blob_status::needs_extraQED); blob->AddToInParticles(inpart); blob->SetTypeSpec("Sherpa"); Decay_Table* table=p_decaymap->FindDecay(blob->InParticle(0)->Flav()); if (table==NULL) { msg_Error()<AddData("dc",new Blob_Data(table->Select())); DEBUG_INFO("p_onshell="<Momentum()); blob->AddData("p_onshell",new Blob_Data(inpart->Momentum())); DEBUG_INFO("succeeded."); } void Hard_Decay_Handler::FindDecayProducts(Particle* decayer, list& decayprods) { if (decayer->DecayBlob()==NULL) { decayprods.push_back(decayer); } else { for (size_t i=0; iDecayBlob()->NOutP(); ++i) { FindDecayProducts(decayer->DecayBlob()->OutParticle(i), decayprods); } } } double Hard_Decay_Handler::BRFactor(ATOOLS::Blob* blob) const { double brfactor=1.0; for (size_t i=0; iNOutP(); ++i) { Particle* part=blob->OutParticle(i); Decay_Table* dt=p_decaymap->FindDecay(part->RefFlav()); if (dt) { brfactor*=dt->ActiveWidth(0)/dt->TotalWidth(); if (part->DecayBlob() && part->DecayBlob()->Type()==btp::Hard_Decay) brfactor*=BRFactor(part->DecayBlob()); } } return brfactor; } void Hard_Decay_Handler::TreatInitialBlob(ATOOLS::Blob* blob, METOOLS::Amplitude2_Tensor* amps, const Particle_Vector& origparts) { Decay_Handler_Base::TreatInitialBlob(blob, amps, origparts); double brfactor=m_br_weights ? BRFactor(blob) : 1.0; DEBUG_VAR(brfactor); Blob_Data_Base * bdbmeweight((*blob)["MEWeight"]); if (bdbmeweight) { // msg_Out()<Get()<<", " // <<"BR = "<InParticle(0)->Flav()<<" " // <InParticle(1)->Flav()<<" -->"; // for (size_t i=0;iNOutP();i++) // msg_Out()<<" "<OutParticle(i)->Flav(); // msg_Out()<<".\n"; bdbmeweight->Set(brfactor*bdbmeweight->Get()); } Blob_Data_Base * wgtinfo((*blob)["MEWeightInfo"]); if (wgtinfo) *wgtinfo->Get()*=brfactor; Blob_Data_Base * wgtmap_bdb((*blob)["WeightsMap"]); if (wgtmap_bdb) wgtmap_bdb->Get()["BR"]=brfactor; NLO_subevtlist* sublist(NULL); Blob_Data_Base * bdb((*blob)["NLO_subeventlist"]); if (bdb) sublist=bdb->Get(); if (sublist) { // If the blob contains a NLO_subeventlist, we have to attach decays // in the sub-events as well. The decay has to be identical for infrared // cancellations, so we simply take each decay and boost it to the sub // kinematics to replace the particle in the subevent DEBUG_FUNC(""); vector > decayprods(blob->NOutP()); size_t newn(2); for (size_t i=0; iNOutP()-1; ++i) { // iterate over out-particles excluding real emission parton list decayprods_i; FindDecayProducts(blob->OutParticle(i), decayprods_i); DEBUG_VAR(blob->OutParticle(i)->Flav()); list::const_iterator it; for (it=decayprods_i.begin(); it!=decayprods_i.end(); ++it) { DEBUG_VAR((*it)->Flav()); } decayprods[i]=decayprods_i; newn+=decayprods_i.size(); } DEBUG_VAR(newn); if (p_newsublist) { for (size_t i=0; isize(); ++i) delete (*p_newsublist)[i]; p_newsublist->clear(); } else p_newsublist=new NLO_subevtlist(); for (size_t i=0; isize(); ++i) { // iterate over sub events and replace decayed particles NLO_subevt* sub((*sublist)[i]); DEBUG_VAR(*(*sublist)[i]); if (sub->IsReal()) newn+=1; Flavour* newfls = new Flavour[newn]; Vec4D* newmoms = new Vec4D[newn]; size_t* newid = new size_t[newn]; for (size_t n=0; nm_n=newn; newsub->p_id=newid; newsub->p_fl=newfls; newsub->p_mom=newmoms; newsub->m_delete=true; p_newsublist->push_back(newsub); int nin=blob->NInP(); for (size_t j=0, jnew=0; jNOutP()-1; ++j, ++jnew) { if (j<2 || decayprods[j-nin].size()==1) { newfls[jnew]=sub->p_fl[j]; newmoms[jnew]=sub->p_mom[j]; continue; } if (sub->p_fl[j]!=blob->OutParticle(j-nin)->Flav()) { THROW(fatal_error, "Internal Error 1"); } Vec4D oldmom=blob->OutParticle(j-nin)->Momentum(); Vec4D newmom=sub->p_mom[j]; Poincare cms(oldmom); Poincare newframe(newmom); newframe.Invert(); list::const_iterator it; for (it=decayprods[j-nin].begin(); it!=decayprods[j-nin].end(); ++it) { newfls[jnew]=(*it)->Flav(); newmoms[jnew]=Vec4D(newframe*(cms*(*it)->Momentum())); jnew++; } jnew--; // because we replaced one particle } if (sub->IsReal()) { // Add remaining parton for real event newfls[newn-1]=sub->p_fl[sub->m_n-1]; newmoms[newn-1]=sub->p_mom[sub->m_n-1]; } } bdb->Set(p_newsublist); DEBUG_INFO("New subevts:"); for (size_t i=0;isize();++i) { (*p_newsublist)[i]->m_result*=brfactor; (*p_newsublist)[i]->m_results["BR"]*=brfactor; (*p_newsublist)[i]->m_me*=brfactor; (*p_newsublist)[i]->m_mewgt*=brfactor; DEBUG_VAR(*(*p_newsublist)[i]); } } } bool Hard_Decay_Handler::DefineInitialConditions(Cluster_Amplitude* ampl, Blob* initial_blob) { DEBUG_FUNC(this); DEBUG_VAR(*ampl); for (int i=0; iNOutP(); ++i) { ampl->Leg(initial_blob->NInP()+i)->SetMom (initial_blob->OutParticle(i)->Momentum()); } if (ampl->NIn()==2) { for (Cluster_Amplitude *campl(ampl); campl;campl=campl->Next()) { if (-campl->Leg(0)->Mom()[0]>rpa->gen.PBeam(0)[0] || -campl->Leg(1)->Mom()[0]>rpa->gen.PBeam(1)[0]) return false; } } size_t imax=ampl->Legs().size()-1; for (int i=0; iNOutP(); ++i) { if (initial_blob->OutParticle(i)->DecayBlob()) { AddDecayClustering(ampl, initial_blob->OutParticle(i)->DecayBlob(), imax, 1<<(initial_blob->NInP()+i)); } } return true; } void Hard_Decay_Handler::AddDecayClustering(ATOOLS::Cluster_Amplitude*& ampl, Blob* blob, size_t& imax, size_t idmother) { DEBUG_FUNC("blob->Id()="<Id()<<" idmother="<GetOutParticles().size();++i) { Particle * p(blob->OutParticle(i)); if (p->Info()=='S') photons.push_back(make_pair(p,p)); else if (p->Info()=='s') splitphotonproducts.push_back(p); else daughters.push_back(p); } msg_Debugging()<<"daughters: "; for (size_t i(0);iFlav().IDName()<<" "; msg_Debugging()<<" + "<InitPrev(); copy->CopyFrom(ampl); copy->SetNLO(0); copy->SetFlag(1); copy->SetMS(ampl->MS()); Cluster_Leg *lij(ampl->IdLeg(idmother)); copy->SetKT2(lij->Mom().Abs2()); for (size_t i=0; iLegs().size(); ++i) ampl->Leg(i)->SetStat(ampl->Leg(i)->Stat()|1); lij->SetStat(1|2|4); size_t idk(0); for (size_t i=0; iLegs().size(); ++i) { copy->Leg(i)->SetK(0); if (copy->Leg(i)->Id()!=idmother && (copy->Leg(i)->Col().m_i==lij->Col().m_j || copy->Leg(i)->Col().m_j==lij->Col().m_i)) idk=copy->Leg(i)->Id(); } if (lij->Col().m_i==0 && lij->Col().m_j==0) { // Ad hoc EW partner size_t ampl_nout=ampl->Legs().size()-ampl->NIn(); if (ampl_nout==1) idk=ampl->Leg(0)->Id(); else { size_t select(0); do { select=ampl->NIn()+floor(ran->Get()*ampl_nout); } while (ampl->Leg(select)->Id()&idmother || select>ampl->Legs().size()-1); idk=ampl->Leg(select)->Id(); } } if (idk==0) THROW(fatal_error,"Colour partner not found"); lij->SetK(idk); Cluster_Leg *d1(copy->IdLeg(idmother)); size_t stat1(0), stat2(0); d1->SetMom(RecombinedMomentum(daughters[0],photons,splitphotons,stat1)); d1->SetStat(stat1); d1->SetFlav(daughters[0]->Flav()); d1->SetFromDec(true); copy->CreateLeg(RecombinedMomentum(daughters[1],photons,splitphotons,stat2), daughters[1]->RefFlav()); copy->Legs().back()->SetFromDec(true); size_t idnew=1<<(++imax); copy->Legs().back()->SetId(idnew); copy->Legs().back()->SetStat(stat2); Cluster_Amplitude::SetColours(ampl->IdLeg(idmother), copy->IdLeg(idmother), copy->Legs().back()); copy->SetIdNew(idnew); DEBUG_VAR(*copy); Cluster_Amplitude* tmp=copy; while (tmp->Next()) { tmp=tmp->Next(); if (tmp->IdNew()&idmother) tmp->SetIdNew(tmp->IdNew()|idnew); for (size_t i=0; iLegs().size(); ++i) { if (tmp->Leg(i)->Id()&idmother) { tmp->Leg(i)->SetId(tmp->Leg(i)->Id()|idnew); } if (tmp->Leg(i)->K()&idmother) { tmp->Leg(i)->SetK(tmp->Leg(i)->K()|idnew); } } DEBUG_VAR(*tmp); } std::vector ids; ids.push_back(idmother); ids.push_back(idnew); // always combine radiated (split) photons with // identified primary charged decay particle while (splitphotons.size()) AddSplitPhotonsClustering(copy, daughters, splitphotons, imax, ids); while (photons.size()) AddPhotonsClustering(copy, daughters, photons, imax, ids); if (daughters[0]->DecayBlob()) AddDecayClustering(copy, daughters[0]->DecayBlob(), imax, idmother); if (daughters[1]->DecayBlob()) AddDecayClustering(copy, daughters[1]->DecayBlob(), imax, idnew); ampl=copy; } else if (daughters.size()==3) { msg_Debugging()<<"1 to 3 case"< 0 P[->1 2] // propagator always combines daughters 1+2 Cluster_Amplitude* step1=ampl->InitPrev(); step1->CopyFrom(ampl); step1->SetNLO(0); step1->SetFlag(1); step1->SetMS(ampl->MS()); Cluster_Leg *lij(ampl->IdLeg(idmother)); step1->SetKT2(lij->Mom().Abs2()); if (!lij) THROW(fatal_error,"Cluster leg of id "+ToString(idmother) +" not found."); for (size_t i=0; iLegs().size(); ++i) ampl->Leg(i)->SetStat(ampl->Leg(i)->Stat()|1); lij->SetStat(1|2|4); size_t idk(0); for (size_t i=0; iLegs().size(); ++i) { step1->Leg(i)->SetK(0); if (step1->Leg(i)->Id()!=idmother) if (step1->Leg(i)->Col().m_i==lij->Col().m_j || step1->Leg(i)->Col().m_j==lij->Col().m_i) idk=step1->Leg(i)->Id(); } if (lij->Col().m_i==0 && lij->Col().m_j==0) { // Ad hoc EW partner size_t ampl_nout=ampl->Legs().size()-ampl->NIn(); if (ampl_nout==1) idk=ampl->Leg(0)->Id(); else { size_t select(0); do { select=ampl->NIn()+floor(ran->Get()*ampl_nout); } while (ampl->Leg(select)->Id()&idmother || select>ampl->Legs().size()-1); idk=ampl->Leg(select)->Id(); } } if (idk==0) THROW(fatal_error,"Colour partner not found"); lij->SetK(idk); Cluster_Leg *d1(step1->IdLeg(idmother)); size_t stat1(0),stat2(0),stat3(0); d1->SetMom(RecombinedMomentum(daughters[0],photons,splitphotons,stat1)); d1->SetStat(stat1); d1->SetFlav(daughters[0]->Flav()); d1->SetFromDec(true); // todo: 1->2 qcd shower with ew fs recoil partner // d1->SetK(idmother);// not that simple: w->qq' has color connection in fs Decay_Channel* dc(NULL); Blob_Data_Base* data = (*blob)["dc"]; if (data) { dc=data->Get(); DEBUG_VAR(*dc); } else THROW(fatal_error, "Internal error."); Comix1to3* amp=dynamic_cast(dc->GetDiagrams()[0]); if (!amp) THROW(fatal_error, "Internal error."); Flavour prop_flav=amp->Prop(); Vec4D momd2=RecombinedMomentum(daughters[1],photons,splitphotons,stat2); Vec4D momd3=RecombinedMomentum(daughters[2],photons,splitphotons,stat3); Vec4D prop_mom=momd2+momd3; step1->CreateLeg(prop_mom, prop_flav); size_t idnew1=1<<(++imax); step1->Legs().back()->SetId(idnew1); step1->Legs().back()->SetStat(0); step1->Legs().back()->SetFromDec(true); Cluster_Amplitude::SetColours(ampl->IdLeg(idmother), step1->IdLeg(idmother), step1->Legs().back()); step1->SetIdNew(idnew1); DEBUG_VAR(*step1); Cluster_Amplitude* tmp=step1; while (tmp->Next()) { tmp=tmp->Next(); if (tmp->IdNew()&idmother) tmp->SetIdNew(tmp->IdNew()|idnew1); for (size_t i=0; iLegs().size(); ++i) { if (tmp->Leg(i)->Id()&idmother) { tmp->Leg(i)->SetId(tmp->Leg(i)->Id()|idnew1); } if (tmp->Leg(i)->K()&idmother) { tmp->Leg(i)->SetK(tmp->Leg(i)->K()|idnew1); } } DEBUG_VAR(*tmp); } Cluster_Amplitude* step2=step1->InitPrev(); step2->CopyFrom(step1); step2->SetNLO(0); step2->SetFlag(1); step2->SetMS(step1->MS()); for (size_t i=0; iLegs().size(); ++i) step1->Leg(i)->SetStat(step1->Leg(i)->Stat()|1); step1->IdLeg(idnew1)->SetStat(1|4); step1->IdLeg(idnew1)->SetK(idk); for (size_t i=0; iLegs().size(); ++i) step2->Leg(i)->SetK(0); Cluster_Leg *d2(step2->IdLeg(idnew1)); d2->SetMom(momd2); d2->SetStat(stat2); d2->SetFlav(daughters[1]->Flav()); step2->CreateLeg(momd3, daughters[2]->Flav()); size_t idnew2=1<<(++imax); step2->Legs().back()->SetId(idnew2); step2->Legs().back()->SetStat(stat3); step2->Legs().back()->SetFromDec(true); Cluster_Amplitude::SetColours(step1->IdLeg(idnew1), step2->IdLeg(idnew1), step2->Legs().back()); step2->SetIdNew(idnew2); DEBUG_VAR(*step2); tmp=step2; while (tmp->Next()) { tmp=tmp->Next(); if (tmp->IdNew()&(idmother|idnew1)) tmp->SetIdNew(tmp->IdNew()|idnew2); for (size_t i=0; iLegs().size(); ++i) { if (tmp->Leg(i)->Id()&idnew1) { tmp->Leg(i)->SetId(tmp->Leg(i)->Id()|idnew2); } if (tmp->Leg(i)->K()&idnew1) { tmp->Leg(i)->SetK(tmp->Leg(i)->K()|idnew2); } } DEBUG_VAR(*tmp); } std::vector ids; ids.push_back(idmother); ids.push_back(idnew1); ids.push_back(idnew2); // always combine radiated (split) photons with // identified primary charged decay particle while (splitphotons.size()) AddSplitPhotonsClustering(step2,daughters,splitphotons,imax,ids); while (photons.size()) AddPhotonsClustering(step2,daughters,photons,imax,ids); if (daughters[0]->DecayBlob()) AddDecayClustering(step2, daughters[0]->DecayBlob(), imax, idmother); if (daughters[1]->DecayBlob()) AddDecayClustering(step2, daughters[1]->DecayBlob(), imax, idnew1); if (daughters[2]->DecayBlob()) AddDecayClustering(step2, daughters[2]->DecayBlob(), imax, idnew2); ampl=step2; } else { PRINT_VAR(*blob); THROW(fatal_error, "1 -> n not implemented yet."); } } void Hard_Decay_Handler::AddSplitPhotonsClustering(Cluster_Amplitude*& ampl, const Particle_Vector daughters, ParticlePairPair_Vector& splitphotons, size_t& imax, const std::vector& ids) { DEBUG_FUNC(splitphotons.size()<<" split photons to be clustered"); // will need to construct two cluster steps // 1) the (offshell) photon is radiated // 2) the (offshell) photon splits into the identified pair Particle * splitphoton1(splitphotons.back().first.first); Particle * splitphoton2(splitphotons.back().first.second); Particle * daughter(splitphotons.back().second); splitphotons.pop_back(); size_t idmother(0),idphoton(0); if (daughter==daughters[0]) idmother=ids[0]; else if (daughter==daughters[1]) idmother=ids[1]; else if (daughter==daughters[2]) idmother=ids[2]; else THROW(fatal_error,"Did not find id for "+daughter->Flav().IDName()); // construct recombined (off-shell) photon momentum Vec4D pmom(splitphoton1->Momentum()+splitphoton2->Momentum()); msg_Debugging()<<"Cluster recombined photon with "<Flav()<<" "<InitPrev(); copy->CopyFrom(ampl); copy->SetNLO(0); copy->SetFlag(1); copy->SetMS(ampl->MS()); Cluster_Leg *lij(ampl->IdLeg(idmother)); for (size_t i=0; iLegs().size(); ++i) ampl->Leg(i)->SetStat(ampl->Leg(i)->Stat()|1); lij->SetStat(1|2|4); size_t idk(0); for (size_t i=0; iLegs().size(); ++i) copy->Leg(i)->SetK(0); if (lij->Col().m_i!=0 || lij->Col().m_j!=0) THROW(fatal_error,"Adding QED to coloured particle."); // Ad hoc QED partner, must not be another soft photon size_t ampl_nout=ampl->Legs().size()-ampl->NIn(); if (ampl_nout==1) idk=ampl->Leg(0)->Id(); else { size_t select(0); size_t nvalid(0); for (size_t i(ampl->NIn());iLegs().size();++i) { if (!(ampl->Leg(i)->Id()&idmother || i>ampl->Legs().size()-1 || ampl->Leg(i)->Flav().Kfcode()==kf_photon)) { nvalid++; } } if (nvalid==0) select=0; else { do { select=ampl->NIn()+floor(ran->Get()*ampl_nout); } while (ampl->Leg(select)->Id()&idmother || select>ampl->Legs().size()-1 || ampl->Leg(select)->Flav().Kfcode()==kf_photon); } msg_Debugging()<<"choose ("<Leg(select)->Id())<<") " <Leg(select)->Flav()<Leg(select)->Id(); } if (idk==0) THROW(fatal_error,"Colour partner not found"); lij->SetK(idk); Cluster_Leg *d1(copy->IdLeg(idmother)); size_t stat1(0), stat2(0); d1->SetMom(RecombinedMomentum(daughter,ParticlePair_Vector(),splitphotons,stat1)); d1->SetStat(stat1); d1->SetFlav(daughter->Flav()); copy->CreateLeg(pmom,Flavour(kf_photon)); idphoton=1<<(++imax); copy->Legs().back()->SetId(idphoton); copy->Legs().back()->SetStat(stat2); Cluster_Amplitude::SetColours(ampl->IdLeg(idmother), copy->IdLeg(idmother), copy->Legs().back()); copy->SetIdNew(idphoton); DEBUG_VAR(*copy); // update IDs Cluster_Amplitude* tmp=copy; while (tmp->Next()) { tmp=tmp->Next(); if (tmp->IdNew()&idmother) tmp->SetIdNew(tmp->IdNew()|idphoton); for (size_t i=0; iLegs().size(); ++i) { if (tmp->Leg(i)->Id()&idmother) { tmp->Leg(i)->SetId(tmp->Leg(i)->Id()|idphoton); } if (tmp->Leg(i)->K()&idmother) { tmp->Leg(i)->SetK(tmp->Leg(i)->K()|idphoton); } } DEBUG_VAR(*tmp); } ampl=copy; // construct the photon splitting cluster step copy=ampl->InitPrev(); copy->CopyFrom(ampl); copy->SetNLO(0); copy->SetFlag(1); copy->SetMS(ampl->MS()); lij=ampl->IdLeg(idphoton); for (size_t i=0; iLegs().size(); ++i) ampl->Leg(i)->SetStat(ampl->Leg(i)->Stat()|1); lij->SetStat(1|2|4); idk=0; for (size_t i=0; iLegs().size(); ++i) copy->Leg(i)->SetK(0); if (lij->Col().m_i!=0 || lij->Col().m_j!=0) THROW(fatal_error,"Adding QED to coloured particle."); // Ad hoc QED partner, must not be another soft photon ampl_nout=ampl->Legs().size()-ampl->NIn(); if (ampl_nout==1) idk=ampl->Leg(0)->Id(); else { size_t select(0); size_t nvalid(0); for (size_t i(ampl->NIn());iLegs().size();++i) { if (!(ampl->Leg(i)->Id()&idmother || i>ampl->Legs().size()-1 || ampl->Leg(i)->Flav().Kfcode()==kf_photon)) { nvalid++; } } if (nvalid==0) select=0; else { do { select=ampl->NIn()+floor(ran->Get()*ampl_nout); } while (ampl->Leg(select)->Id()&idmother || select>ampl->Legs().size()-1 || ampl->Leg(select)->Flav().Kfcode()==kf_photon); } msg_Debugging()<<"choose ("<Leg(select)->Id())<<") " <Leg(select)->Flav()<Leg(select)->Id(); } if (idk==0) THROW(fatal_error,"Colour partner not found"); lij->SetK(idk); d1=copy->IdLeg(idphoton); stat1=0; stat2=0; // we know the splitting products are not evolved further yet d1->SetMom(splitphoton1->Momentum()); d1->SetStat(stat1); d1->SetFlav(splitphoton1->Flav()); copy->CreateLeg(splitphoton2->Momentum(),splitphoton2->Flav()); size_t idnew=1<<(++imax); copy->Legs().back()->SetId(idnew); copy->Legs().back()->SetStat(stat2); Cluster_Amplitude::SetColours(ampl->IdLeg(idmother), copy->IdLeg(idmother), copy->Legs().back()); copy->SetIdNew(idnew); DEBUG_VAR(*copy); // update IDs tmp=copy; while (tmp->Next()) { tmp=tmp->Next(); if (tmp->IdNew()&idphoton) tmp->SetIdNew(tmp->IdNew()|idnew); for (size_t i=0; iLegs().size(); ++i) { if (tmp->Leg(i)->Id()&idphoton) { tmp->Leg(i)->SetId(tmp->Leg(i)->Id()|idnew); } if (tmp->Leg(i)->K()&idphoton) { tmp->Leg(i)->SetK(tmp->Leg(i)->K()|idnew); } } DEBUG_VAR(*tmp); } ampl=copy; } void Hard_Decay_Handler::AddPhotonsClustering(Cluster_Amplitude*& ampl, const Particle_Vector daughters, ParticlePair_Vector& photons, size_t& imax, const std::vector& ids) { DEBUG_FUNC(photons.size()<<" photons to be clustered"); Particle * photon(photons.back().first); Particle * daughter(photons.back().second); photons.pop_back(); size_t idmother(0); if (daughter==daughters[0]) idmother=ids[0]; else if (daughter==daughters[1]) idmother=ids[1]; else if (daughter==daughters[2]) idmother=ids[2]; else THROW(fatal_error,"Did not find id for "+daughter->Flav().IDName()); msg_Debugging()<<"Cluster "<Flav()<<" "<Momentum() <<" with "<Flav()<<" "<InitPrev(); copy->CopyFrom(ampl); copy->SetNLO(0); copy->SetFlag(1); copy->SetMS(ampl->MS()); Cluster_Leg *lij(ampl->IdLeg(idmother)); for (size_t i=0; iLegs().size(); ++i) ampl->Leg(i)->SetStat(ampl->Leg(i)->Stat()|1); lij->SetStat(1|2|4); size_t idk(0); for (size_t i=0; iLegs().size(); ++i) { copy->Leg(i)->SetK(0); } if (lij->Col().m_i!=0 || lij->Col().m_j!=0) { THROW(fatal_error,"Adding QED to coloured particle."); } // Ad hoc QED partner, must not be another soft photon size_t ampl_nout=ampl->Legs().size()-ampl->NIn(); if (ampl_nout==1) idk=ampl->Leg(0)->Id(); else { size_t select(0); size_t nvalid(0); for (size_t i(ampl->NIn());iLegs().size();++i) { if (!(ampl->Leg(i)->Id()&idmother || i>ampl->Legs().size()-1 || ampl->Leg(i)->Flav().Kfcode()==kf_photon)) { nvalid++; } } if (nvalid==0) select=0; else { do { select=ampl->NIn()+floor(ran->Get()*ampl_nout); } while (ampl->Leg(select)->Id()&idmother || select>ampl->Legs().size()-1 || ampl->Leg(select)->Flav().Kfcode()==kf_photon); } msg_Debugging()<<"choose ("<Leg(select)->Id())<<") " <Leg(select)->Flav()<Leg(select)->Id(); } if (idk==0) THROW(fatal_error,"Colour partner not found"); lij->SetK(idk); Cluster_Leg *d1(copy->IdLeg(idmother)); size_t stat1(0), stat2(0); d1->SetMom(RecombinedMomentum(daughter,photons,ParticlePairPair_Vector(),stat1)); d1->SetStat(stat1); d1->SetFlav(daughter->Flav()); copy->CreateLeg(photon->Momentum(),photon->RefFlav()); size_t idnew=1<<(++imax); copy->Legs().back()->SetId(idnew); copy->Legs().back()->SetStat(stat2); Cluster_Amplitude::SetColours(ampl->IdLeg(idmother), copy->IdLeg(idmother), copy->Legs().back()); copy->SetIdNew(idnew); DEBUG_VAR(*copy); // update IDs Cluster_Amplitude* tmp=copy; while (tmp->Next()) { tmp=tmp->Next(); if (tmp->IdNew()&idmother) tmp->SetIdNew(tmp->IdNew()|idnew); for (size_t i=0; iLegs().size(); ++i) { if (tmp->Leg(i)->Id()&idmother) { tmp->Leg(i)->SetId(tmp->Leg(i)->Id()|idnew); } if (tmp->Leg(i)->K()&idmother) { tmp->Leg(i)->SetK(tmp->Leg(i)->K()|idnew); } } DEBUG_VAR(*tmp); } ampl=copy; } void Hard_Decay_Handler::UnsplitPhotons(const ATOOLS::Particle_Vector& spp, ParticlePairPair_Vector& splitphotons) { // undo the splitting of a photon, recombine particle-anti-particle pair // into (off-shell) photon if (spp.size()) { if (spp.size()%2!=0) THROW(fatal_error,"Uneven number of photon-splitting products, cannot proceed."); if (spp.size()==2) { splitphotons.push_back(make_pair(make_pair(spp[0],spp[1]),spp[0])); } else { // find (particle,anti-particle)-pairs in spp, // order by invariant mass, starting with the smallest std::map > pairtable; for (size_t i(0);iFlav()==spp[j]->Flav().Bar()) { double m2((spp[i]->Momentum()+spp[j]->Momentum()).Mass()); pairtable[m2]=make_pair(i,j); } } } if (pairtable.empty()) { msg_Debugging()<<"no pairs found"< >::const_iterator it=pairtable.begin();it!=pairtable.end();++it) msg_Debugging()<second.first<<" "<second.second<<": " <<", m2="<first< usedindices; for (std::map >::const_iterator it=pairtable.begin(); it!=pairtable.end();++it) { bool valid(true); for (size_t i(0);isecond.first==usedindices[i] || it->second.second==usedindices[i]) { valid=false; break; } if (!valid) continue; usedindices.push_back(it->second.first); usedindices.push_back(it->second.second); msg_Debugging()<<"constructing split pair: P -> " <second.first]->Flav()<<" " <second.second]->Flav()<<" , m2 = " <first<second.first], spp[it->second.second]), spp[it->second.first])); } if (2*splitphotons.size()!=spp.size()) { msg_Error()<Flav().Charge()) cdaughters.push_back(daughters[i]); if (cdaughters.size()==1) { for (size_t i(0);iMomentum()); cmom+=cmoms[i]; } Poincare ccms(cmom); for (size_t i(0);iMomentum()); ccms.Boost(pmom); size_t id(0); double dR(pmom.DR(cmoms[0])); for (size_t j(1);jFlav()<<" " <Momentum() <<" assigned to "<Flav()<Flav().Charge()) cdaughters.push_back(daughters[i]); if (cdaughters.size()==1) { for (size_t i(0);iMomentum()); cmom+=cmoms[i]; } Poincare ccms(cmom); for (size_t i(0);iMomentum() +splitphotons[i].first.second->Momentum()); ccms.Boost(pmom); size_t id(0); double dR(pmom.DR(cmoms[0])); for (size_t j(1);jFlav()<<" " <Momentum()<<" and " <Flav()<<" " <Momentum() <<" assigned to "<Flav()<Momentum() +splitphotons[i].first.second->Momentum(); stat|=2|4; } } for (size_t i(0);iMomentum(); stat|=2|4; } } msg_Debugging()<Flav()<<": "<Momentum(); } void Hard_Decay_Handler::ReadDecayTable(Flavour decayer) { DEBUG_FUNC(decayer); if (!m_store_results) return; Data_Reader reader = Data_Reader("|",";","!"); reader.AddComment("#"); reader.AddComment("//"); reader.AddWordSeparator("\t"); reader.SetInputPath(m_resultdir); reader.SetInputFile(decayer.ShellName()); vector > file; if(reader.MatrixFromFile(file)) { for (size_t iline=0; iline results(3); for (size_t i=0; i<3; ++i) results[i]=ToType(file[iline][i+1]); m_read[decayer].insert(make_pair(decaychannel, results)); } else { PRINT_INFO("Wrong format in decay table in "<begin(); dmit!=p_decaymap->end(); ++dmit) { ofstream ostr((m_resultdir+"/"+dmit->first.ShellName()).c_str()); ostr<<"# Decay table for "<first<second->begin(); dtit!=dmit->second->end(); ++dtit) { ostr<IDCode()<<"\t" <IWidth()<<"\t" <IDeltaWidth()<<"\t" <Max()<