PHASIC%2B%2B/Decays/Decay_Channel.H

0001 #ifndef Decay_Channel_h
0002 #define Decay_Channel_h
0003
0004 #include "ATOOLS/Phys/Flavour.H"
0005 #include "ATOOLS/Org/Message.H"
0006 #include "ATOOLS/Org/Exception.H"
0007 #include "ATOOLS/Org/MyStrStream.H"
0008 #include "ATOOLS/Math/MyComplex.H"
0009 #include "ATOOLS/Math/Vector.H"
0010 #include "ATOOLS/Phys/Particle.H"
0011 #include <set>
0012 #include "ATOOLS/Org/Scoped_Settings.H"
0013
0014 namespace METOOLS {
0015   class Spin_Amplitudes;
0016   class Amplitude2_Tensor;
0017   class Spin_Density;
0018 }
0019
0020 namespace PHASIC {
0021   class Multi_Channel;
0022   class Single_Channel;
0023
0024   class Decay_Channel {
0025   private:
0026     static std::map<std::string, double> s_kinmaxfails;
0027
0028     /**
0029      * This method makes the phase space generator choose a phase space point
0030      * and calculate the corresponding ME and PS weights.
0031      * The ME weights (\b not squared) for each helicity saved in
0032      * Decay_Channel::m_diagrams are then summed for all available diagrams and
0033      * helicities taking into account the symmetry factor and
0034      * averaging over the incoming helicity states.
0035      * So far, this Differential method does not care about spin correlations
0036      * it merely sums up all absolute squared entries of those
0037      * in order to obtain a value for \f$d\Gamma\f$, and averages over incoming
0038      * spin.
0039      *
0040      * @param momenta Reference to phase space point, which has size=NOut()+1
0041      * and the incoming momentum already set.
0042      * @param anti Whether to consider the charge conjugated matrix elements.
0043      * @param sigma If not NULL, this spin density matrix is considered for the
0044      * incoming particle.
0045      *
0046      * @return \f$d\Gamma(\f$
0047      */
0048     double Differential(ATOOLS::Vec4D_Vector& momenta, bool anti,
0049                         METOOLS::Spin_Density* sigma=NULL,
0050                         const std::vector<ATOOLS::Particle*>& p=
0051                           std::vector<ATOOLS::Particle*>());
0052     double ME2(const ATOOLS::Vec4D_Vector& momenta, bool anti,
0053                METOOLS::Spin_Density* sigma=NULL,
0054                const std::vector<ATOOLS::Particle*>& p=
0055                std::vector<ATOOLS::Particle*>());
0056
0057   protected :
0058     double      m_width, m_deltawidth, m_minmass, m_max, m_symfac;
0059
0060     /// integrated width (as opposed to the one given by BR)
0061     double m_iwidth;
0062     /// integrated deltawidth (as opposed to the one given by BR)
0063     double m_ideltawidth;
0064
0065     /**
0066      * m_active[n] specifies the status of the decay channel for the decay
0067      * of the @noop <n>'th particle of this type in the event.
0068      * It can take the following values: m_active[n]=
0069      *   0: Disabled decay channel, but contributes to total width
0070      *   1: Enabled decay channel
0071      *   2: Forced decay channel
0072      *   -1: Disabled decay channel, and does not contribute to width
0073      */
0074     std::vector<int> m_active;
0075
0076     ATOOLS::Flavour_Vector m_flavours;
0077     std::vector<METOOLS::Spin_Amplitudes*> m_diagrams;
0078     PHASIC::Multi_Channel* p_channels;
0079     METOOLS::Amplitude2_Tensor* p_amps;
0080     /**
0081      * The Mass_Selector knows which type of mass is to be used for all decays
0082      * associated with this decay table. This can be the HadMass() if the decays
0083      * belong to the non-perturbative regime of the event, or the Mass() if
0084      * they belong to the perturbative part.
0085      */
0086     const ATOOLS::Mass_Selector* p_ms;
0087
0088     double Lambda(const double& a, const double& b, const double& c) const;
0089     double MassWeight(const double& s, const double& sp,
0090                       const double& b, const double& c) const;
0091
0092   public :
0093     Decay_Channel(const ATOOLS::Flavour &, const ATOOLS::Mass_Selector* ms);
0094     virtual ~Decay_Channel();
0095     void AddDecayProduct(const ATOOLS::Flavour & flout,const bool & sort=true);
0096     inline void SetWidth(const double & width) { m_width =width; }
0097     inline void SetDeltaWidth(const double & delta) { m_deltawidth=delta; }
0098
0099     std::string Name() const;
0100     std::string IDCode() const;
0101     static std::string IDCode(const ATOOLS::Flavour& decayer,
0102                               const ATOOLS::Flavour_Vector& daughters);
0103     std::string FSIDCode();
0104     inline const double& Width() const { return m_width; }
0105     inline const double& DeltaWidth() const { return m_deltawidth; }
0106     inline const double& IWidth() const { return m_iwidth; }
0107     inline void SetIWidth(const double& width) { m_iwidth=width; }
0108     inline const double& IDeltaWidth() const { return m_ideltawidth; }
0109     inline void SetIDeltaWidth(const double& delta) { m_ideltawidth=delta; }
0110     inline const double& Max() const { return m_max; }
0111     inline void SetMax(const double& max) { m_max=max; }
0112     inline const double& MinimalMass() const { return m_minmass; }
0113     inline int Active(const size_t& count) const { return count<m_active.size() ? m_active[count] : m_active[0]; }
0114     inline const std::vector<int>& Active() const { return m_active; }
0115     inline void SetActive(std::vector<int> active) { m_active=active; }
0116     inline void SetActive(size_t count, int active) {
0117       if (count>=m_active.size()) m_active.resize(count+1,m_active[0]);
0118       m_active[count]=active;
0119     }
0120     inline void SetActiveAll(int active) { for (size_t i=0; i<m_active.size(); ++i) m_active[i]=active; }
0121
0122     double GenerateMass(const double& max, const double &width) const;
0123     double SymmetryFactor();
0124
0125     inline int NOut() const { return m_flavours.size()-1; }
0126     inline const std::vector<ATOOLS::Flavour>& Flavs() const{return m_flavours;}
0127     inline const ATOOLS::Flavour& GetDecaying() const { return m_flavours[0]; }
0128     inline const ATOOLS::Flavour& GetDecayProduct(size_t i) const {
0129       return m_flavours[i+1];
0130     }
0131     inline const std::vector<METOOLS::Spin_Amplitudes*>& GetDiagrams() const
0132     { return m_diagrams; }
0133     void AddDiagram(METOOLS::Spin_Amplitudes* amp);
0134     inline PHASIC::Multi_Channel* Channels() const { return p_channels; }
0135     inline void SetChannels(PHASIC::Multi_Channel* chans) { p_channels=chans; }
0136     void AddChannel(PHASIC::Single_Channel* chan);
0137
0138     inline void ResetDiagrams() { m_diagrams.clear(); }
0139     void ResetChannels();
0140
0141     void Output() const;
0142
0143     void CalculateWidth(double acc=0.01, double ref=0.0, int iter=2500);
0144     void GenerateKinematics(ATOOLS::Vec4D_Vector& momenta, bool anti,
0145                 METOOLS::Spin_Density* sigma=NULL,
0146                 const std::vector<ATOOLS::Particle*>& p=
0147                 std::vector<ATOOLS::Particle*>());
0148     inline METOOLS::Amplitude2_Tensor* Amps() { return p_amps; }
0149
0150     friend std::ostream &operator<<(std::ostream &os, const Decay_Channel &dt);
0151
0152     static void SortFlavours(ATOOLS::Flavour_Vector& flavs);
0153     static bool FlavourSort(const ATOOLS::Flavour &fl1,
0154                             const ATOOLS::Flavour &fl2);
0155     static void PrintMaxKinFailStatistics(std::ostream &str);
0156   };
0157 }
0158 #endif