#include "AddOns/Analysis/Observables/Soft_Photon_Observables.H" #include "AddOns/Analysis/Main/Primitive_Analysis.H" #include "ATOOLS/Math/Poincare.H" #include "ATOOLS/Org/MyStrStream.H" #include "ATOOLS/Org/Run_Parameter.H" #include "ATOOLS/Org/Message.H" using namespace ANALYSIS; using namespace ATOOLS; template Primitive_Observable_Base *GetObservable(const Analysis_Key& key) { ATOOLS::Scoped_Settings s{ key.m_settings }; const auto parameters = s.SetDefault({}).GetVector(); if (parameters.size() < 7) THROW(missing_input, "Missing parameter values."); if (parameters.size() > 15) THROW(missing_input, "Too many parameter values."); std::vector f(parameters.size()-5); for (size_t i=0;i(parameters[i]); f[i]=ATOOLS::Flavour((kf_code)abs(kf)); if (kf<0) f[i]=f[i].Bar(); } std::string list=parameters[parameters.size()-1]; return new Class(f,HistogramType(parameters[parameters.size()-2]), s.Interprete(parameters[parameters.size()-5]), s.Interprete(parameters[parameters.size()-4]), s.Interprete(parameters[parameters.size()-3]),list); } #define DEFINE_GETTER_METHOD(CLASS,NAME) \ Primitive_Observable_Base * \ ATOOLS::Getter::operator()(const Analysis_Key& key) const \ { return GetObservable(key); } #define DEFINE_PRINT_METHOD(NAME) \ void ATOOLS::Getter::PrintInfo(std::ostream &str,const size_t width) const \ { str<<"[kf1, ..., kfn, min, max, bins, Lin|LinErr|Log|LogErr, list] ... 1& flavs, int type, double xmin, double xmax, int nbins, const std::string& listname, const std::string& name) : Primitive_Observable_Base(type,xmin,xmax,nbins), f_special(false) { if(flavs.size()<2) { msg_Error()<<"Error in Soft_Photon_Observable_Base:"<=0.0) f_special=true; } void Soft_Photon_Observable_Base::Evaluate(double value, double weight, double ncount) { p_histo->Insert(value,weight,ncount); } void Soft_Photon_Observable_Base::Evaluate(int nout, const Vec4D* moms, const Flavour* flavs, double weight, double ncount) { msg_Out()<<"Flavour array method"<Type()==btp::QED_Radiation && blobs[i]->TypeSpec()==m_blobtype) { QEDblob=blobs[i]; break; } } if (!QEDblob) { Evaluate(0.,weight,ncount); return; } Particle_Vector multipole,photons; for (int i=0;iNOutP();++i) { for (size_t j=0;jOutParticle(i)->Flav()==m_flavs[j]) { multipole.push_back(QEDblob->OutParticle(i)); continue; } } if (QEDblob->OutParticle(i)->Info()=='S') { photons.push_back(QEDblob->OutParticle(i)); } } Evaluate(multipole,photons,weight,ncount); } //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DEFINE_OBSERVABLE_GETTER(Soft_Photon_Energy, Soft_Photon_Energy_Getter, "SoftPhotonEnergy") void Soft_Photon_Energy::Evaluate(const ATOOLS::Particle_Vector& multipole, const ATOOLS::Particle_Vector& photons, double weight, double ncount) { Vec4D decayer(0.,0.,0.,0.); for (size_t i=0;iMomentum(); } Vec4D sumphot(0.,0.,0.,0.); for (size_t i=0;iMomentum(); } decayer+=sumphot; Poincare decframe(decayer); decframe.Boost(sumphot); p_histo->Insert(sumphot[0],weight,ncount); } Soft_Photon_Energy::Soft_Photon_Energy(const std::vector& flavs, int type,double xmin,double xmax, int nbins, const std::string & listname) : Soft_Photon_Observable_Base(flavs,type,xmin,xmax,nbins,listname, "SoftPhotonEnergy") {} Primitive_Observable_Base* Soft_Photon_Energy::Copy() const { return new Soft_Photon_Energy(m_flavs,m_type,m_xmin,m_xmax,m_nbins, m_listname); } //============================================================================= DEFINE_OBSERVABLE_GETTER(Soft_Photon_Angle,Soft_Photon_Angle_Getter, "SoftPhotonAngle") void Soft_Photon_Angle::Evaluate(const ATOOLS::Particle_Vector& multipole, const ATOOLS::Particle_Vector& photons, double weight, double ncount) { Vec4D multipolesum(0.,0.,0.,0.); Vec4D photonsum(0.,0.,0.,0.); Vec4D axpart(0.,0.,0.,1.); Vec4D axis(0.,0.,0.,1.); int chargesum(0); // FS charged momentum sum for (size_t i=0;iMomentum(); chargesum+=multipole[i]->Flav().Charge(); } // photon momentum sum for (size_t i=0;iMomentum(); } // neutral resonance -> nothing further // charged resonance -> add to multipole, change definition for theta=0 if (chargesum!=0) { // add charged resonance to multipole Vec4D resonance(multipolesum+photonsum); multipolesum+=resonance; // resonance at theta=0 axpart=resonance; } Poincare multipoleboost(multipolesum); multipoleboost.Boost(axpart); Poincare rotation(axpart,axis); for (size_t i=0;iMomentum(); multipoleboost.Boost(k); rotation.Rotate(k); double theta = acos((Vec3D(k)*Vec3D(axis))/(Vec3D(k).Abs())); p_histo->Insert(theta,weight,ncount); } } Soft_Photon_Angle::Soft_Photon_Angle(const std::vector& flavs, int type, double xmin, double xmax, int nbins, const std::string & listname) : Soft_Photon_Observable_Base(flavs,type,xmin,xmax,nbins,listname, "SoftPhotonAngle") {} Primitive_Observable_Base* Soft_Photon_Angle::Copy() const { return new Soft_Photon_Angle(m_flavs,m_type,m_xmin,m_xmax,m_nbins, m_listname); }