LightFlavorRatios/unbinnedMaxLikelihood_fitter/Kshort_ScaleWidthDifferential.C

0001 #ifndef KSHORT_SCALEWIDTHDIFFERENTIAL_C
0002 #define KSHORT_SCALEWIDTHDIFFERENTIAL_C
0003
0004 #include <algorithm>
0005 #include <cctype>
0006 #include <cstdint>
0007 #include <format>
0008 #include <functional>
0009 #include <iomanip>
0010 #include <iostream>
0011 #include <regex>
0012 #include <sstream>
0013 #include <string>
0014 #include <vector>
0015
0016 #include <TCut.h>
0017 #include <TSystem.h>
0018
0019 #include "./RDataframeToRoofit.C"
0020 #include "./RDataframeToRoofit_SOnly.C"
0021 #include "./fitutil.h"
0022
0023 namespace
0024 {
0025 inline void replace_all(std::string &s, const std::string &from, const std::string &to)
0026 {
0027     if (from.empty())
0028         return;
0029     size_t pos = 0;
0030     while ((pos = s.find(from, pos)) != std::string::npos)
0031     {
0032         s.replace(pos, from.size(), to);
0033         pos += to.size();
0034     }
0035 }
0036
0037 inline std::string to_hex_u64(std::uint64_t v)
0038 {
0039     std::ostringstream oss;
0040     oss << std::hex << std::setw(16) << std::setfill('0') << v;
0041     return oss.str();
0042 }
0043
0044 inline std::string shorten_tag_if_needed(const std::string &tag, std::size_t maxlen = 180)
0045 {
0046     if (tag.size() <= maxlen)
0047         return tag;
0048
0049     const std::uint64_t h = static_cast<std::uint64_t>(std::hash<std::string>{}(tag));
0050     const std::string suffix = "_h" + to_hex_u64(h);
0051
0052     const std::size_t keep = (maxlen > suffix.size()) ? (maxlen - suffix.size()) : 0;
0053     return tag.substr(0, keep) + suffix;
0054 }
0055
0056 inline std::string trim_number(const std::string &text)
0057 {
0058     std::ostringstream oss;
0059     oss << std::fixed << std::setprecision(6) << std::stod(text);
0060
0061     std::string out = oss.str();
0062     while (out.size() > 1 && out.back() == '0')
0063         out.pop_back();
0064     if (!out.empty() && out.back() == '.')
0065         out.push_back('0');
0066     return out;
0067 }
0068
0069 struct CutRangeText
0070 {
0071     bool has_low = false;
0072     bool has_high = false;
0073     std::string low;
0074     std::string high;
0075     std::string low_bracket = "[";
0076     std::string high_bracket = ")";
0077 };
0078
0079 inline CutRangeText parse_cut_range(const std::string &cut_string, const std::string &variable)
0080 {
0081     CutRangeText range;
0082     const std::regex re(variable + R"(\s*(>=|<=|==|>|<)\s*([-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?))");
0083
0084     for (std::sregex_iterator it(cut_string.begin(), cut_string.end(), re), end; it != end; ++it)
0085     {
0086         const std::string op = (*it)[1].str();
0087         const std::string value = trim_number((*it)[2].str());
0088
0089         if (op == ">=" || op == ">")
0090         {
0091             range.has_low = true;
0092             range.low = value;
0093             range.low_bracket = (op == ">=") ? "[" : "(";
0094         }
0095         else if (op == "<=" || op == "<")
0096         {
0097             range.has_high = true;
0098             range.high = value;
0099             range.high_bracket = (op == "<=") ? "]" : ")";
0100         }
0101         else if (op == "==")
0102         {
0103             range.has_low = true;
0104             range.has_high = true;
0105             range.low = value;
0106             range.high = value;
0107             range.low_bracket = "[";
0108             range.high_bracket = "]";
0109         }
0110     }
0111
0112     return range;
0113 }
0114
0115 inline void add_range_text(std::vector<std::string> &entries, const CutRangeText &range, const std::string &label, const std::string &separator, const std::string &unit)
0116 {
0117     if (range.has_low && range.has_high)
0118     {
0119         entries.push_back(label + " = " + range.low_bracket + range.low + separator + range.high + range.high_bracket + unit);
0120     }
0121     else if (range.has_low)
0122     {
0123         entries.push_back(label + " #geq " + range.low + unit);
0124     }
0125     else if (range.has_high)
0126     {
0127         entries.push_back(label + " < " + range.high + unit);
0128     }
0129 }
0130
0131 inline std::vector<std::string> make_cut_legend_entries(const std::string &cut_string)
0132 {
0133     std::vector<std::string> entries;
0134
0135     add_range_text(entries, parse_cut_range(cut_string, "K_S0_pT"), "p_{T}", "-", " GeV");
0136     add_range_text(entries, parse_cut_range(cut_string, "K_S0_pseudorapidity"), "#eta", ", ", "");
0137     add_range_text(entries, parse_cut_range(cut_string, "K_S0_phi"), "#phi", ", ", " rad");
0138
0139     const std::regex min_track_pt_re(R"(min\s*\(\s*track_1_pT\s*,\s*track_2_pT\s*\)\s*(>=|<=|==|>|<)\s*([-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?))");
0140     std::smatch min_track_pt_match;
0141     if (std::regex_search(cut_string, min_track_pt_match, min_track_pt_re))
0142     {
0143         const std::string op = min_track_pt_match[1].str();
0144         const std::string value = trim_number(min_track_pt_match[2].str());
0145         if (op == ">=")
0146             entries.push_back("min track p_{T} " + value + " GeV");
0147         else
0148             entries.push_back("min track p_{T} " + op + " " + value + " GeV");
0149     }
0150
0151     return entries;
0152 }
0153 } // namespace
0154
0155 std::string modify(std::string str)
0156 {
0157     replace_all(str, ">=", "geq");
0158     replace_all(str, "<=", "leq");
0159     replace_all(str, "==", "eq");
0160     replace_all(str, "&&", "AND");
0161     replace_all(str, "||", "OR");
0162     replace_all(str, ">", "gt");
0163     replace_all(str, "<", "lt");
0164
0165     std::replace(str.begin(), str.end(), '.', 'p'); // decimal point
0166     std::replace(str.begin(), str.end(), '-', 'm'); // minus sign (negative numbers)
0167     std::replace(str.begin(), str.end(), '+', 'p'); // plus sign
0168
0169     for (char &c : str)
0170     {
0171         const bool ok = (std::isalnum(static_cast<unsigned char>(c)) != 0) || (c == '_');
0172         if (!ok)
0173             c = '_';
0174     }
0175
0176     std::string out;
0177     out.reserve(str.size());
0178     bool prev_us = false;
0179     for (char c : str)
0180     {
0181         if (c == '_')
0182         {
0183             if (!prev_us)
0184                 out.push_back(c);
0185             prev_us = true;
0186         }
0187         else
0188         {
0189             out.push_back(c);
0190             prev_us = false;
0191         }
0192     }
0193
0194     // Trim leading/trailing underscores
0195     while (!out.empty() && out.front() == '_')
0196         out.erase(out.begin());
0197     while (!out.empty() && out.back() == '_')
0198         out.pop_back();
0199
0200     return out;
0201 }
0202
0203 void configure_fitparam(const std::string cutstring, fitparam_config &fit_conf)
0204 {
0205     if (cutstring == "K_S0_pT>=0.60&&K_S0_pT<0.70&&K_S0_phi>=1.00&&K_S0_phi<2.00")
0206     {
0207         // Mass window and binning
0208         fit_conf.minMass = 0.44;
0209         fit_conf.maxMass = 0.56;
0210         fit_conf.nBins = 60;
0211
0212         // Signal starting params / ranges
0213         fit_conf.sigmodel = "Gaussian";
0214         fit_conf.mean = 0.5;
0215         fit_conf.mean_low = 0.485;
0216         fit_conf.mean_high = 0.51;
0217
0218         fit_conf.sigma = 5e-3;
0219         fit_conf.sigma_low = 1e-3;
0220         fit_conf.sigma_high = 8e-3;
0221
0222         // Background model
0223         fit_conf.bkgmodel = "Polynomial3";
0224         fit_conf.p1 = -1.6724e+00;
0225         fit_conf.p1_low = -10.;
0226         fit_conf.p1_high = 10.;
0227
0228         // fit_conf.p2 = 0;
0229         fit_conf.p2_low = -10.;
0230         fit_conf.p2_high = 10.;
0231
0232         // fit_conf.p3 = 0;
0233         fit_conf.p3_low = -10.;
0234         fit_conf.p3_high = 10.;
0235
0236         fit_conf.nSig = 0.0;     // if > 0, will be used as fixed starting value
0237         fit_conf.nSig_low = 0.0; // if > 0, will be used as fixed starting value
0238         fit_conf.nBkg = 0.0;     // if > 0, will be used as fixed starting value
0239         fit_conf.nBkg_low = 0.0; // if > 0, will be used as fixed starting value
0240     }
0241     else if (cutstring == "K_S0_pT>=0.60&&K_S0_pT<0.70&&K_S0_phi>=2.00&&K_S0_phi<3.15")
0242     {
0243         // Mass window and binning
0244         fit_conf.minMass = 0.44;
0245         fit_conf.maxMass = 0.56;
0246         fit_conf.nBins = 60;
0247
0248         // Signal starting params / ranges
0249         fit_conf.sigmodel = "Gaussian";
0250         fit_conf.mean = 0.488;
0251         fit_conf.mean_low = 0.48;
0252         fit_conf.mean_high = 0.51;
0253
0254         fit_conf.sigma = 4e-3;
0255         fit_conf.sigma_low = 1e-3;
0256         fit_conf.sigma_high = 8e-3;
0257
0258         // Background model
0259         fit_conf.bkgmodel = "Polynomial3";
0260         fit_conf.p1 = -1.6724e+00;
0261         fit_conf.p1_low = -10.;
0262         fit_conf.p1_high = 10.;
0263
0264         // fit_conf.p2 = 0;
0265         fit_conf.p2_low = -10.;
0266         fit_conf.p2_high = 10.;
0267
0268         // fit_conf.p3 = 0;
0269         fit_conf.p3_low = -10.;
0270         fit_conf.p3_high = 10.;
0271
0272         fit_conf.nSig = 0.0;     // if > 0, will be used as fixed starting value
0273         fit_conf.nSig_low = 0.0; // if > 0, will be used as fixed starting value
0274         fit_conf.nBkg = 0.0;     // if > 0, will be used as fixed starting value
0275         fit_conf.nBkg_low = 0.0; // if > 0, will be used as fixed starting value
0276     }
0277     else if (cutstring == "K_S0_pT>=0.80&&K_S0_pT<0.90&&K_S0_phi>=-1.00&&K_S0_phi<0.00")
0278     {
0279         // Mass window and binning
0280         fit_conf.minMass = 0.44;
0281         fit_conf.maxMass = 0.56;
0282         fit_conf.nBins = 60;
0283
0284         // Signal starting params / ranges
0285         fit_conf.sigmodel = "Gaussian";
0286         fit_conf.mean = 0.488;
0287         fit_conf.mean_low = 0.48;
0288         fit_conf.mean_high = 0.51;
0289
0290         fit_conf.sigma = 4e-3;
0291         fit_conf.sigma_low = 1e-3;
0292         fit_conf.sigma_high = 8e-3;
0293
0294         // Background model
0295         fit_conf.bkgmodel = "Polynomial3";
0296         fit_conf.p1 = -1.6724e+00;
0297         fit_conf.p1_low = -10.;
0298         fit_conf.p1_high = 10.;
0299
0300         // fit_conf.p2 = 0;
0301         fit_conf.p2_low = -10.;
0302         fit_conf.p2_high = 10.;
0303
0304         // fit_conf.p3 = 0;
0305         fit_conf.p3_low = -10.;
0306         fit_conf.p3_high = 10.;
0307
0308         fit_conf.nSig = 0.0;     // if > 0, will be used as fixed starting value
0309         fit_conf.nSig_low = 0.0; // if > 0, will be used as fixed starting value
0310         fit_conf.nBkg = 0.0;     // if > 0, will be used as fixed starting value
0311         fit_conf.nBkg_low = 0.0; // if > 0, will be used as fixed starting value
0312     }
0313     else if (cutstring == "K_S0_pT>=0.90&&K_S0_pT<1.00&&K_S0_phi>=-2.00&&K_S0_phi<-1.00")
0314     {
0315         // Mass window and binning
0316         fit_conf.minMass = 0.44;
0317         fit_conf.maxMass = 0.56;
0318         fit_conf.nBins = 60;
0319
0320         // Signal starting params / ranges
0321         fit_conf.sigmodel = "Gaussian";
0322         fit_conf.mean = 0.488;
0323         fit_conf.mean_low = 0.48;
0324         fit_conf.mean_high = 0.51;
0325
0326         fit_conf.sigma = 4e-3;
0327         fit_conf.sigma_low = 1e-3;
0328         fit_conf.sigma_high = 8e-3;
0329
0330         // Background model
0331         fit_conf.bkgmodel = "Polynomial3";
0332         fit_conf.p1 = -1.6724e+00;
0333         fit_conf.p1_low = -10.;
0334         fit_conf.p1_high = 10.;
0335
0336         // fit_conf.p2 = 0;
0337         fit_conf.p2_low = -10.;
0338         fit_conf.p2_high = 10.;
0339
0340         // fit_conf.p3 = 0;
0341         fit_conf.p3_low = -10.;
0342         fit_conf.p3_high = 10.;
0343
0344         fit_conf.nSig = 0.0;     // if > 0, will be used as fixed starting value
0345         fit_conf.nSig_low = 0.0; // if > 0, will be used as fixed starting value
0346         fit_conf.nBkg = 0.0;     // if > 0, will be used as fixed starting value
0347         fit_conf.nBkg_low = 0.0; // if > 0, will be used as fixed starting value
0348     }
0349     else if (cutstring == "K_S0_pT>=1.50&&K_S0_pT<1.80&&K_S0_phi>=0.00&&K_S0_phi<1.00")
0350     {
0351         // Mass window and binning
0352         fit_conf.minMass = 0.44;
0353         fit_conf.maxMass = 0.56;
0354         fit_conf.nBins = 60;
0355
0356         // Signal starting params / ranges
0357         fit_conf.sigmodel = "Gaussian";
0358         fit_conf.mean = 0.488;
0359         fit_conf.mean_low = 0.48;
0360         fit_conf.mean_high = 0.51;
0361
0362         fit_conf.sigma = 4e-3;
0363         fit_conf.sigma_low = 1e-3;
0364         fit_conf.sigma_high = 15e-3;
0365
0366         // Background model
0367         fit_conf.bkgmodel = "Polynomial3";
0368         fit_conf.p1 = -1.6724e+00;
0369         fit_conf.p1_low = -10.;
0370         fit_conf.p1_high = 10.;
0371
0372         // fit_conf.p2 = 0;
0373         fit_conf.p2_low = -10.;
0374         fit_conf.p2_high = 10.;
0375
0376         // fit_conf.p3 = 0;
0377         fit_conf.p3_low = -10.;
0378         fit_conf.p3_high = 10.;
0379
0380         fit_conf.nSig = 0.0;     // if > 0, will be used as fixed starting value
0381         fit_conf.nSig_low = 0.0; // if > 0, will be used as fixed starting value
0382         fit_conf.nBkg = 0.0;     // if > 0, will be used as fixed starting value
0383         fit_conf.nBkg_low = 0.0; // if > 0, will be used as fixed starting value
0384     }
0385     else
0386     {
0387         // Mass window and binning
0388         fit_conf.minMass = 0.44;
0389         fit_conf.maxMass = 0.56;
0390         fit_conf.nBins = 60;
0391
0392         // Signal starting params / ranges
0393         fit_conf.sigmodel = "Gaussian";
0394         fit_conf.mean = 0.5;
0395         fit_conf.mean_low = 0.485;
0396         fit_conf.mean_high = 0.51;
0397
0398         fit_conf.sigma = 1e-2;
0399         fit_conf.sigma_low = 1e-3;
0400         fit_conf.sigma_high = 5e-1;
0401
0402         // Tail params (used for CB/DoubleCB only)
0403         fit_conf.alpha1 = 0.5;
0404         fit_conf.alpha1_low = 1E-3;
0405         fit_conf.alpha1_high = 5.0;
0406
0407         fit_conf.n1 = 10;
0408         fit_conf.n1_low = 5;
0409         fit_conf.n1_high = 1E3;
0410
0411         fit_conf.alpha2 = 0.5;
0412         fit_conf.alpha2_low = 1E-3;
0413         fit_conf.alpha2_high = 5.0;
0414
0415         fit_conf.n2 = 10;
0416         fit_conf.n2_low = 5;
0417         fit_conf.n2_high = 1E3;
0418
0419         fit_conf.frac = 0.5;
0420
0421         // Background model
0422         fit_conf.bkgmodel = "Polynomial3";
0423         fit_conf.p1 = -1.6724e+00;
0424         fit_conf.p1_low = -10.;
0425         fit_conf.p1_high = 10.;
0426
0427         // fit_conf.p2 = 0;
0428         fit_conf.p2_low = -10.;
0429         fit_conf.p2_high = 10.;
0430
0431         // fit_conf.p3 = 0;
0432         fit_conf.p3_low = -10.;
0433         fit_conf.p3_high = 10.;
0434
0435         // fitconfig.nSig_fracTol = 0.7;
0436         // fitconfig.nSig_fracTol_low = 0.5;
0437         // fitconfig.nSig_fracTol_high = 1.00;
0438         // fitconfig.nBkg_fracTol = 1.0 - fitconfig.nSig_fracTol;
0439         // fitconfig.nBkg_fracTol_low = 1.0 - fitconfig.nSig_fracTol_high;
0440         // fitconfig.nBkg_fracTol_high = 1.0;
0441         fit_conf.nSig = 0.0;     // if > 0, will be used as fixed starting value
0442         fit_conf.nSig_low = 0.0; // if > 0, will be used as fixed starting value
0443         fit_conf.nBkg = 0.0;     // if > 0, will be used as fixed starting value
0444         fit_conf.nBkg_low = 0.0; // if > 0, will be used as fixed starting value
0445     }
0446 }
0447
0448 void Kshort_ScaleWidthDifferential(                                                           //
0449     std::string additionalCuts = "K_S0_pT>=1.50&&K_S0_pT<1.80&&K_S0_phi>=0.00&&K_S0_phi<1.00" //
0450 )
0451 {
0452     const std::string figure_dir = "./figure_kshort_crosscheck/";
0453     system(("mkdir -p " + figure_dir).c_str());
0454
0455     const bool isMC = false;
0456
0457     // TCut KS0_optimizeCut = "K_S0_DIRA >= 0.985000 && K_S0_decayLength > 0.125000 && min(abs(track_1_IP_xy), abs(track_2_IP_xy)) >= 0.020000";
0458     const TCut KS0_optimizeCut = "min(track_1_pT, track_2_pT) >= 0.2";
0459     const TCut KSO_additionalCuts(additionalCuts.c_str());
0460     const std::vector<std::string> cut_legend_entries = make_cut_legend_entries(additionalCuts + "&&" + KS0_optimizeCut.GetTitle());
0461
0462     const TCut totalselection = KS0_optimizeCut && KSO_additionalCuts;
0463     std::cout << "Selection: " << totalselection.GetTitle() << std::endl;
0464
0465     std::string plotdir_prefix = modify(additionalCuts);
0466     plotdir_prefix = shorten_tag_if_needed(plotdir_prefix);
0467     std::cout << "Plotdir prefix: " << plotdir_prefix << std::endl;
0468
0469     const std::string plotdir = figure_dir + "/" + plotdir_prefix;
0470     system(("mkdir -p " + plotdir).c_str());
0471
0472     // Fit configuration defaults
0473     fitparam_config fitconfig;
0474     fitconfig.branch = "K_S0_mass";
0475     fitconfig.decaystring = "K_{S}^{0}#rightarrow#pi^{-}#pi^{+}";
0476
0477     // Signal model
0478     // fitconfig.sigmodel = "DoubleCrystalBall";
0479     // fitconfig.sigmodel = "Gaussian";
0480
0481     configure_fitparam(additionalCuts, fitconfig);
0482
0483     // // Mass window and binning
0484     // fitconfig.minMass = (isMC) ? 0.45 : 0.44;
0485     // fitconfig.maxMass = (isMC) ? 0.54 : 0.56;
0486     // fitconfig.nBins = (isMC) ? 90 : 60;
0487
0488     // // Signal starting params / ranges
0489     // fitconfig.mean = 0.5;
0490     // fitconfig.mean_low = 0.485;
0491     // fitconfig.mean_high = 0.51;
0492
0493     // fitconfig.sigma = 1e-2;
0494     // fitconfig.sigma_low = 1e-3;
0495     // fitconfig.sigma_high = 5e-1;
0496
0497     // // Tail params (used for CB/DoubleCB only)
0498     // fitconfig.alpha1 = 0.5;
0499     // fitconfig.alpha1_low = 1E-3;
0500     // fitconfig.alpha1_high = 5.0;
0501
0502     // fitconfig.n1 = 10;
0503     // fitconfig.n1_low = 5;
0504     // fitconfig.n1_high = 1E3;
0505
0506     // fitconfig.alpha2 = 0.5;
0507     // fitconfig.alpha2_low = 1E-3;
0508     // fitconfig.alpha2_high = 5.0;
0509
0510     // fitconfig.n2 = 10;
0511     // fitconfig.n2_low = 5;
0512     // fitconfig.n2_high = 1E3;
0513
0514     // fitconfig.frac = 0.5;
0515
0516     // // Background model
0517     // fitconfig.bkgmodel = "Polynomial3";
0518     // fitconfig.p1 = -1.6724e+00;
0519     // fitconfig.p1_low = -10.;
0520     // fitconfig.p1_high = 10.;
0521
0522     // // fitconfig.p2 = 0;
0523     // fitconfig.p2_low = -10.;
0524     // fitconfig.p2_high = 10.;
0525
0526     // // fitconfig.p3 = 0;
0527     // fitconfig.p3_low = -10.;
0528     // fitconfig.p3_high = 10.;
0529
0530     // // fitconfig.nSig_fracTol = 0.7;
0531     // // fitconfig.nSig_fracTol_low = 0.5;
0532     // // fitconfig.nSig_fracTol_high = 1.00;
0533     // // fitconfig.nBkg_fracTol = 1.0 - fitconfig.nSig_fracTol;
0534     // // fitconfig.nBkg_fracTol_low = 1.0 - fitconfig.nSig_fracTol_high;
0535     // // fitconfig.nBkg_fracTol_high = 1.0;
0536     // fitconfig.nSig = 0.0;     // if > 0, will be used as fixed starting value
0537     // fitconfig.nSig_low = 0.0; // if > 0, will be used as fixed starting value
0538     // fitconfig.nBkg = 0.0;     // if > 0, will be used as fixed starting value
0539     // fitconfig.nBkg_low = 0.0; // if > 0, will be used as fixed starting value
0540
0541     // fitconfig.useStagedFit = false;
0542
0543     if (!isMC)
0544     {
0545         const bool doSnapshot = true;
0546         const std::string snapshotName = std::format("./snapshot/filtered_tree_{}.root", plotdir_prefix);
0547         // const std::string inputFile = "/sphenix/user/hjheng/sPHENIXRepo/TrackingAnalysis/Physics_Val_TF/misc/KFParticle_merged/merged_pipi_reco_00079516.root";
0548         const std::string inputFile = "/sphenix/tg/tg01/hf/aopatton/SVLooseJun4/6RunsCombinedKShortSVLoose.root";
0549
0550         const double significance = RDataframeToRoofit(doSnapshot, snapshotName, inputFile, totalselection, fitconfig, plotdir, cut_legend_entries);
0551
0552         std::cout << "Returned significance: " << significance << std::endl;
0553     }
0554     else
0555     {
0556         const bool doSnapshot = false;
0557         const std::string snapshotName = "./snapshot/filtered_tree.root";
0558         const std::string inputFile = "/sphenix/user/hjheng/sPHENIXRepo/TrackingAnalysis/Physics_Val_TF/misc/KFParticle_merged/merged_pipi_reco_00079516.root";
0559
0560         RDataframeToRoofit_SOnly(doSnapshot, snapshotName, inputFile, totalselection, fitconfig, plotdir);
0561     }
0562 }
0563
0564 #endif // KSHORT_SCALEWIDTHDIFFERENTIAL_C