fastjet/contrib/MassFlav.hh

0001 ///
0002 /// \file MassFlav.hh
0003 ///
0004 /// This file is mainly intended for IRC tests of the IFNPlugin
0005 #ifndef __FJCONTRIB_MASSFLAV_HH__
0006 #define __FJCONTRIB_MASSFLAV_HH__
0007
0008 #ifdef __FJC_FLAVINFO_USEFJCORE__
0009 #include "FlavInfo.hh"
0010 #else
0011 #include "fastjet/PseudoJet.hh"
0012 #include "fastjet/JetDefinition.hh"
0013 #include "fastjet/contrib/FlavNeutraliser.hh"
0014 #endif
0015
0016 #include <limits>
0017 ///
0018 ///
0019
0020 FASTJET_BEGIN_NAMESPACE  // defined in fastjet/internal/base.hh
0021 namespace contrib{
0022
0023   //class MassFlavHistory : public FlavInfo {
0024   class MassFlavHistory : public FlavHistory {
0025    public:
0026     // MassFlavHistory(FlavInfo flavinfo, double m_in = 0.0) : FlavInfo(flavinfo),
0027     // _m(m_in) {}
0028     MassFlavHistory(FlavHistory flavhist, double m_in = 0.0)
0029         : FlavHistory(flavhist), _m(m_in) {}
0030     void set_mass(double m_in) { _m = m_in; }
0031     double mass() const { return _m; }
0032
0033    private:
0034     double _m;
0035   };
0036
0037   class MassFlavRecombiner : public FlavRecombiner {
0038    public:
0039     MassFlavRecombiner() : FlavRecombiner() {}
0040
0041     MassFlavRecombiner(const FlavRecombiner & flav_reco) : FlavRecombiner(flav_reco) {}
0042
0043
0044     /// Perform a recombination taking flavour into account
0045     /// as well as the known mass of the particles, so
0046     /// as to be able to correctly reconstruct the rapidity
0047     /// of the newly formed particle, even when pa and pb
0048     /// are at quite extreme rapidities.
0049     void recombine(const PseudoJet &pa, const PseudoJet &pb,
0050                    PseudoJet &pab) const {
0051
0052       // old code for the case eta_a << 0, eta_b >> 0
0053       // not sure whether it's actually needed
0054
0055       //PseudoJet pab_tilde;
0056       //FlavRecombiner::recombine(pa, pb, pab_tilde);
0057
0058       //double pab_m2 = ma*ma + mb*mb + 2*dot_product(pa,pb);
0059
0060       //// and set the correct rapidity
0061       //double rap = 0.5*log((max(0.0,pab_m2)+pab_tilde.kt2())/
0062       //                      pow(pab_tilde.E() + fabs(pab_tilde.pz()),2));
0063       //if (pab_tilde.pz() > 0) rap = -rap;
0064
0065       //pab = PtYPhiM(pab_tilde.pt(), rap, pab_tilde.phi(), sqrt(pab_m2));
0066       //pab.set_user_info(new MassFlavHistory(pab_tilde.user_info<FlavHistory>(),
0067       //                  sqrt(pab_m2)));
0068
0069       // special treatment if both rapidities are small, or if one
0070       // is small and it looks like the other will be small too
0071       constexpr double rap_lim = 0.1;
0072
0073       // if both particles (or one particle and its result) are at
0074       // small rapidity, just use 4-vector sum, which preserves
0075       // accurate small x,y,z components more reliably
0076       bool pa_small_rap = std::fabs(pa.rap()) < rap_lim;
0077       bool pb_small_rap = std::fabs(pb.rap()) < rap_lim;
0078       double pab_approx_rap = (pa.pt()*pa.rap() + pb.pt()*pb.rap())/(pa.pt()+pb.pt());
0079       bool pab_small_rap = std::fabs(pab_approx_rap) < rap_lim;
0080       int n_small_rap = int(pa_small_rap ) + int(pb_small_rap) + int(pab_small_rap);
0081       // at least two of the above at small rapidity: revert to 4-vector sum
0082       if (n_small_rap >= 2) {
0083         // open question as to what we do about the mass here...
0084         // (i.e. we could try to make it more accurate, but will not for now)
0085         FlavRecombiner::recombine(pa, pb, pab);
0086         pab.set_user_info(new MassFlavHistory(pab.user_info<FlavHistory>(),
0087                                            pab.m()));
0088         return;
0089       }
0090
0091       double ma = mass_of_particle(pa);
0092       double mb = mass_of_particle(pb);
0093
0094       double avrap = (pa.rap() + pb.rap()) / 2;
0095
0096       PseudoJet shifted_pa = PtYPhiM(pa.pt(), pa.rap() - avrap, pa.phi(), ma);
0097       PseudoJet shifted_pb = PtYPhiM(pb.pt(), pb.rap() - avrap, pb.phi(), mb);
0098
0099       shifted_pa.set_user_info(
0100           new MassFlavHistory(pa.user_info<FlavHistory>(), ma));
0101       shifted_pb.set_user_info(
0102           new MassFlavHistory(pb.user_info<FlavHistory>(), mb));
0103
0104       PseudoJet shifted_pab;
0105
0106       // Recombine using the default recombiner
0107       FlavRecombiner::recombine(shifted_pa, shifted_pb, shifted_pab);
0108
0109       pab = PtYPhiM(shifted_pab.pt(), shifted_pab.rap() + avrap,
0110                     shifted_pab.phi(), shifted_pab.m());
0111       pab.set_user_info(new MassFlavHistory(shifted_pab.user_info<FlavHistory>(),
0112                                         shifted_pab.m()));
0113     }
0114
0115     double mass_of_particle(const PseudoJet &p) const {
0116       if (p.has_user_info<MassFlavHistory>())
0117         return p.user_info<MassFlavHistory>().mass();
0118       // otherwise we expect the mass to be zero to within rounding errors
0119       // so check this
0120       double m2 = p.m2();
0121       const double safety_factor = 10.0;
0122       if (std::fabs(m2) > safety_factor *
0123                               std::numeric_limits<double>::epsilon() *
0124                               pow(p.E(), 2)) {
0125         throw Error(
0126             "MassFlavRecombiner found particle without MassFlavHistory, whose "
0127             "mass is inconsistent with zero");
0128       }
0129       // if things are OK, then just return zero.
0130       return 0.0;
0131     }
0132   };
0133
0134 } // namespace contrib
0135
0136 FASTJET_END_NAMESPACE  // defined in fastjet/internal/base.hh
0137
0138 #endif // __FJCONTRIB_MASSFLAV_HH__