Corral/src/CQResonances.h

0001 #ifndef CQRESONANCES_H
0002 #define CQRESONANCES_H
0003 //---- CQResonances.h -- the two-body decays seen as peaks in C(Q), and the variable Q binning of the ir2=2 block
0004 //----   (Q_inv, pT_pair) built from them (README_meanpT_from_CQ.md sec 11). One list for both, so the C(Q)
0005 //----   signposts (DrawCQSignposts, corral_loop.cxx) and the fine Q bins always agree.
0006 //---- Needs fluct_common.h (particle IDs) included before it.
0007 #include <vector>
0008 #include <cmath>
0009 #include <algorithm>
0010 #include <string>
0011 #include <cstdio>
0012
0013 //---- Q_inv (= dq = 2k*, CalcRm::PairInfo) at which a two-body decay R -> 1+2 lands:
0014 //----   dq = sqrt((M^2-(m1+m2)^2)(M^2-(m1-m2)^2))/M; -1 below threshold. Order of the two legs doesn't matter.
0015 inline double CQResonanceQ(double M, double m1, double m2){
0016     double a    = (M*M-(m1+m2)*(m1+m2))*(M*M-(m1-m2)*(m1-m2));
0017     return (a>0.) ? sqrt(a)/M : -1.;
0018 }
0019
0020 //---- species a, b (either order), PDG mass M and width G (GeV; 0 = weak decay), leg masses m1, m2, label
0021 struct CQResonance { int a, b; double M, G, m1, m2; const char* lab; };
0022 namespace CQRes {
0023     const double mpi = 0.13957, mK0 = 0.497611, mL = 1.115683, mp = 0.938272, mK = 0.493677;        // PDG
0024     const int PIP = kParticleIDPionPlus, PIM = kParticleIDPionMinus;
0025     const int PP = kParticleIDProton, PM = kParticleIDAntiProton, KP = kParticleIDKaonPlus, KM = kParticleIDKaonMinus;
0026     const int KS = kParticleIDKshort, LA = kParticleIDLambda, AL = kParticleIDAntiLambda;
0027 }
0028 const CQResonance CQ_RESONANCES[]   = {
0029     { CQRes::PIP, CQRes::PIM, 0.497611, 0.,       CQRes::mpi, CQRes::mpi, "K^{0}_{S}"        },
0030     { CQRes::PIP, CQRes::PIM, 0.77526,  0.1474,   CQRes::mpi, CQRes::mpi, "#rho^{0}"         },
0031     { CQRes::PIP, CQRes::PIM, 0.990,    0.055,    CQRes::mpi, CQRes::mpi, "f_{0}(980)"       },
0032     { CQRes::PIP, CQRes::PIM, 1.2754,   0.1866,   CQRes::mpi, CQRes::mpi, "f_{2}(1270)"      },
0033     { CQRes::KS,  CQRes::PIP, 0.89167,  0.0514,   CQRes::mK0, CQRes::mpi, "K*^{+}"           },
0034     { CQRes::KS,  CQRes::PIM, 0.89167,  0.0514,   CQRes::mK0, CQRes::mpi, "K*^{-}"           },
0035     { CQRes::LA,  CQRes::PIP, 1.38280,  0.0360,   CQRes::mL,  CQRes::mpi, "#Sigma*^{+}"      },
0036     { CQRes::LA,  CQRes::PIM, 1.3872,   0.0394,   CQRes::mL,  CQRes::mpi, "#Sigma*^{-}"      },
0037     { CQRes::LA,  CQRes::PIM, 1.32171,  0.,       CQRes::mL,  CQRes::mpi, "#Xi^{-}"          },
0038     { CQRes::AL,  CQRes::PIM, 1.38280,  0.0360,   CQRes::mL,  CQRes::mpi, "#bar{#Sigma*}^{-}"},
0039     { CQRes::AL,  CQRes::PIP, 1.3872,   0.0394,   CQRes::mL,  CQRes::mpi, "#bar{#Sigma*}^{+}"},
0040     { CQRes::AL,  CQRes::PIP, 1.32171,  0.,       CQRes::mL,  CQRes::mpi, "#bar{#Xi}^{+}"    },
0041     { CQRes::PP,  CQRes::PIM, 1.115683, 0.,       CQRes::mp,  CQRes::mpi, "#Lambda"          },     // Lambdas not reconstructed as V0s (README_PID.md)
0042     { CQRes::PP,  CQRes::PIM, 1.232,    0.117,    CQRes::mp,  CQRes::mpi, "#Delta^{0}"       },
0043     { CQRes::PP,  CQRes::PIP, 1.232,    0.117,    CQRes::mp,  CQRes::mpi, "#Delta^{++}"      },
0044     { CQRes::PM,  CQRes::PIP, 1.115683, 0.,       CQRes::mp,  CQRes::mpi, "#bar{#Lambda}"    },
0045     { CQRes::PM,  CQRes::PIP, 1.232,    0.117,    CQRes::mp,  CQRes::mpi, "#bar{#Delta}^{0}" },
0046     { CQRes::PM,  CQRes::PIM, 1.232,    0.117,    CQRes::mp,  CQRes::mpi, "#bar{#Delta}^{--}"},
0047     { CQRes::KP,  CQRes::KM,  1.019461, 0.004249, CQRes::mK,  CQRes::mK,  "#phi"             },
0048     //---- 2026-10-06: the 46-pairtype table (PDG masses; a0(980)->KK0s and f0(980)->K0sK0s sit below threshold)
0049     { CQRes::KP,  CQRes::PIM, 0.89555,  0.0473,   CQRes::mK,  CQRes::mpi, "K*^{0}"           },
0050     { CQRes::KM,  CQRes::PIP, 0.89555,  0.0473,   CQRes::mK,  CQRes::mpi, "#bar{K*}^{0}"     },
0051     { CQRes::KP,  CQRes::PIM, 1.4324,   0.109,    CQRes::mK,  CQRes::mpi, "K_{2}*^{0}(1430)" },
0052     { CQRes::KM,  CQRes::PIP, 1.4324,   0.109,    CQRes::mK,  CQRes::mpi, "#bar{K}_{2}*^{0}(1430)"},
0053     { CQRes::KM,  CQRes::PIP, 1.86484,  0.,       CQRes::mK,  CQRes::mpi, "D^{0}"            },
0054     { CQRes::KP,  CQRes::PIM, 1.86484,  0.,       CQRes::mK,  CQRes::mpi, "#bar{D}^{0}"      },
0055     { CQRes::KP,  CQRes::KM,  1.86484,  0.,       CQRes::mK,  CQRes::mK,  "D^{0}"            },
0056     { CQRes::KM,  CQRes::PP,  1.5194,   0.0157,   CQRes::mK,  CQRes::mp,  "#Lambda(1520)"    },
0057     { CQRes::KP,  CQRes::PM,  1.5194,   0.0157,   CQRes::mK,  CQRes::mp,  "#bar{#Lambda}(1520)"},
0058     { CQRes::KM,  CQRes::PP,  1.820,    0.080,    CQRes::mK,  CQRes::mp,  "#Lambda(1820)"    },
0059     { CQRes::KP,  CQRes::PM,  1.820,    0.080,    CQRes::mK,  CQRes::mp,  "#bar{#Lambda}(1820)"},
0060     { CQRes::KS,  CQRes::PIP, 1.86966,  0.,       CQRes::mK0, CQRes::mpi, "D^{+}"            },
0061     { CQRes::KS,  CQRes::PIM, 1.86966,  0.,       CQRes::mK0, CQRes::mpi, "D^{-}"            },
0062     { CQRes::KS,  CQRes::KP,  1.86966,  0.,       CQRes::mK0, CQRes::mK,  "D^{+}"            },
0063     { CQRes::KS,  CQRes::KM,  1.86966,  0.,       CQRes::mK0, CQRes::mK,  "D^{-}"            },
0064     { CQRes::KS,  CQRes::KP,  1.96835,  0.,       CQRes::mK0, CQRes::mK,  "D_{s}^{+}"        },
0065     { CQRes::KS,  CQRes::KM,  1.96835,  0.,       CQRes::mK0, CQRes::mK,  "D_{s}^{-}"        },
0066     { CQRes::KS,  CQRes::PP,  2.28646,  0.,       CQRes::mK0, CQRes::mp,  "#Lambda_{c}^{+}"  },
0067     { CQRes::KS,  CQRes::PM,  2.28646,  0.,       CQRes::mK0, CQRes::mp,  "#bar{#Lambda}_{c}^{-}"},
0068     { CQRes::LA,  CQRes::PIP, 2.28646,  0.,       CQRes::mL,  CQRes::mpi, "#Lambda_{c}^{+}"  },
0069     { CQRes::AL,  CQRes::PIM, 2.28646,  0.,       CQRes::mL,  CQRes::mpi, "#bar{#Lambda}_{c}^{-}"},
0070     { CQRes::LA,  CQRes::KM,  1.67245,  0.,       CQRes::mL,  CQRes::mK,  "#Omega^{-}"       },
0071     { CQRes::AL,  CQRes::KP,  1.67245,  0.,       CQRes::mL,  CQRes::mK,  "#bar{#Omega}^{+}" },
0072     { CQRes::LA,  CQRes::KM,  1.690,    0.020,    CQRes::mL,  CQRes::mK,  "#Xi(1690)^{-}"    },
0073     { CQRes::AL,  CQRes::KP,  1.690,    0.020,    CQRes::mL,  CQRes::mK,  "#bar{#Xi}(1690)^{+}"},
0074     { CQRes::LA,  CQRes::KM,  1.823,    0.024,    CQRes::mL,  CQRes::mK,  "#Xi(1820)^{-}"    },
0075     { CQRes::AL,  CQRes::KP,  1.823,    0.024,    CQRes::mL,  CQRes::mK,  "#bar{#Xi}(1820)^{+}"},
0076     { CQRes::LA,  CQRes::KS,  1.690,    0.020,    CQRes::mL,  CQRes::mK0, "#Xi(1690)^{0}"    },
0077     { CQRes::AL,  CQRes::KS,  1.690,    0.020,    CQRes::mL,  CQRes::mK0, "#bar{#Xi}(1690)^{0}"},
0078     { CQRes::LA,  CQRes::KS,  1.823,    0.024,    CQRes::mL,  CQRes::mK0, "#Xi(1820)^{0}"    },
0079     { CQRes::AL,  CQRes::KS,  1.823,    0.024,    CQRes::mL,  CQRes::mK0, "#bar{#Xi}(1820)^{0}"},
0080     { CQRes::KS,  CQRes::KS,  1.5174,   0.086,    CQRes::mK0, CQRes::mK0, "f'_{2}(1525)"     },
0081     { CQRes::KS,  CQRes::KS,  1.733,    0.123,    CQRes::mK0, CQRes::mK0, "f_{0}(1710)"      },
0082 };
0083 const int NCQ_RESONANCES    = sizeof(CQ_RESONANCES)/sizeof(CQ_RESONANCES[0]);
0084 //---- h+, h- (2026-10-11) carry the pion mass, so they take the pion entries: h+h- shows the pi+pi- signposts and bins
0085 inline int CQResonancePid(int pid){ return (pid==kParticleIDhPlus) ? kParticleIDPionPlus : ((pid==kParticleIDhMinus) ? kParticleIDPionMinus : pid); }
0086 inline bool CQResonanceMatches(const CQResonance &r, int pid1, int pid2){
0087     pid1 = CQResonancePid(pid1); pid2 = CQResonancePid(pid2);
0088     return (r.a==pid1 && r.b==pid2) || (r.a==pid2 && r.b==pid1);
0089 }
0090
0091 //---- ir2=2 axes (README_meanpT_from_CQ.md sec 11, user 2026-10-10)
0092 //---- Q_inv: 0-2 GeV (as hQsib), VARIABLE WIDTH, per species pair. Built on a 20 MeV grid: every 20 MeV cell is
0093 //----   split into the finest bin width any of these rules asks for, so every bin nests inside the uniform
0094 //----   20 MeV bins and every 1 MeV bin of hQsib nests inside it:
0095 //----   - Q < 0.10 GeV: 2 MeV; 0.10-0.20 GeV: 5 MeV (Coulomb and short-range correlations)
0096 //----   - each CQ_RESONANCES entry of this species pair at Q0: width_Q = max(Gamma_Q/2.355, sigma_Q) with the
0097 //----     PLACEHOLDER resolution sigma_Q = max(2 MeV, 0.03 Q0) (no final tracking yet; written on the pages);
0098 //----     window Q0 +- (5 width_Q + 0.03 Q0) (the last term for the momentum scale); bin = the largest of
0099 //----     1, 2, 4, 5, 10, 20 MeV that is <= width_Q/4 (never underbin a sharp peak)
0100 //----   - elsewhere 20 MeV
0101 //---- pT_pair = |pT1+pT2| (= 2 kT): 0-5 GeV/c, 50 uniform bins.
0102 const double QPT_QMAX       = 2.0;
0103 const int    QPT_NPT        = 50;
0104 const double QPT_PTMAX      = 5.0;
0105 const char* const QPT_NOTE  = "variable width binning in Q";
0106 inline std::vector<double> QBinEdges(int pid1, int pid2){
0107     const int NCELL = (int)lround(QPT_QMAX/0.020);
0108     std::vector<int> w(NCELL,20);                       // bin width (MeV) per 20 MeV cell
0109     for (int c=0;c<NCELL;c++){ double q = 0.020*c; if (q<0.0999) w[c] = 2; else if (q<0.1999) w[c] = 5; }
0110     const int allowed[6]    = {1,2,4,5,10,20};
0111     for (int i=0;i<NCQ_RESONANCES;i++){
0112         const CQResonance &r    = CQ_RESONANCES[i];
0113         if (!CQResonanceMatches(r,pid1,pid2)) continue;
0114         double q0   = CQResonanceQ(r.M,r.m1,r.m2);
0115         if (q0<=0. || q0>=QPT_QMAX) continue;
0116         double gq   = 0.;
0117         if (r.G>0.){ double qh = CQResonanceQ(r.M+r.G/2.,r.m1,r.m2), ql = CQResonanceQ(r.M-r.G/2.,r.m1,r.m2); gq = qh - std::max(ql,0.); }
0118         double sig  = std::max(0.002,0.03*q0);
0119         double wq   = std::max(gq/2.355,sig);
0120         double hw   = 5.*wq + 0.03*q0;
0121         int bw      = 1; for (int a : allowed) if (a<=1000.*wq/4.+1e-9) bw = a;
0122         int c1      = std::max(0,(int)floor((q0-hw)/0.020)), c2 = std::min(NCELL-1,(int)floor((q0+hw)/0.020));
0123         for (int c=c1;c<=c2;c++) w[c] = std::min(w[c],bw);
0124     }
0125     std::vector<double> e;
0126     for (int c=0;c<NCELL;c++) for (int k=0;k<20/w[c];k++) e.push_back(0.001*(20*c+k*w[c]));
0127     e.push_back(QPT_QMAX);
0128     return e;
0129 }
0130 //---- the Q edges as segments "0.000-0.100 @2" ... for the log
0131 inline std::string QBinSummary(const std::vector<double> &e){
0132     std::string s; char buf[64];
0133     size_t i = 0;
0134     while (i+1<e.size()){
0135         int w   = (int)lround(1000.*(e[i+1]-e[i]));
0136         size_t j = i+1; while (j+1<e.size() && (int)lround(1000.*(e[j+1]-e[j]))==w) j++;
0137         snprintf(buf,sizeof(buf),"%s%.3f-%.3f @%d",s.empty()?"":", ",e[i],e[j],w); s += buf;
0138         i = j;
0139     }
0140     return s;
0141 }
0142 #endif