Implement/macro/PtCalcMLTutorial.C

0001 #include <TSystem.h>
0002 #include <TStopwatch.h>
0003 #include <iostream>
0004 #include <vector>
0005
0006 // #include "PtCalculator.h"  // SiCaloPt::PtCalculator & friends
0007 #include "/sphenix/user/jzhang1/testcode4all/INTT-EMCAL/InttSeedingTrackDev/ML4Reco/Implement/src/PtCalculator.h"  // SiCaloPt::PtCalculator & friends
0008 R__LOAD_LIBRARY(/sphenix/user/jzhang1/testcode4all/INTT-EMCAL/InttSeedingTrackDev/ML4Reco/Implement/src/libPtCalc.so)
0009
0010 // ---- Weights(onnx) and Scalers(json) Path ---------------------------
0011 struct DemoPaths
0012 {
0013     std::string emd_onnx          = "/sphenix/user/jzhang1/testcode4all/INTT-EMCAL/InttSeedingTrackDev/ML4Reco/Implement/ML_Weight_Scaler/model_MLEMD.onnx";
0014     std::string emd_scaler_json   = "";
0015
0016     std::string eproj_onnx        = "/sphenix/user/jzhang1/testcode4all/INTT-EMCAL/InttSeedingTrackDev/ML4Reco/Implement/ML_Weight_Scaler/model_MLEproj.onnx";
0017     std::string eproj_scaler_json = "/sphenix/user/jzhang1/testcode4all/INTT-EMCAL/InttSeedingTrackDev/ML4Reco/Implement/ML_Weight_Scaler/scaler_MLEproj.json";
0018
0019     std::string combined_onnx         = "/sphenix/user/jzhang1/testcode4all/INTT-EMCAL/InttSeedingTrackDev/ML4Reco/Implement/ML_Weight_Scaler/model_MLCombined.onnx";
0020     std::string combined_scaler_json  = "";
0021 };
0022
0023 // ---- turn string into optional<string> for Config ------------------------
0024 template<typename Opt>
0025 Opt make_opt(const std::string& s)
0026 {
0027     if (s.empty()) return std::nullopt;
0028     return s;
0029 }
0030
0031 // ---- PtCalc Tutorial -------------------------------------
0032 void PtCalcMLTutorial()
0033 {
0034     // Load PtCalc shared library
0035     // gSystem->Load("libPtCalc.so");
0036     // gSystem->Load("/sphenix/user/jzhang1/testcode4all/INTT-EMCAL/InttSeedingTrackDev/ML4Reco/Implement/src/libPtCalc.so");
0037
0038     // Use appropriate paths in your environment, the default "DemoPaths" setup is correct here for mine
0039     DemoPaths WS_Path;
0040
0041     // SiCaloPt::PtCalculatorConfig setting(if you want to use the ML models, Formula-method donnot need these)
0042     SiCaloPt::PtCalculatorConfig cfg;
0043     cfg.mlEMD_model_path        = make_opt<decltype(cfg.mlEMD_model_path)>(WS_Path.emd_onnx);
0044     cfg.mlEMD_scaler_json       = make_opt<decltype(cfg.mlEMD_scaler_json)>(WS_Path.emd_scaler_json);
0045     cfg.mlEproj_model_path      = make_opt<decltype(cfg.mlEproj_model_path)>(WS_Path.eproj_onnx);
0046     cfg.mlEproj_scaler_json     = make_opt<decltype(cfg.mlEproj_scaler_json)>(WS_Path.eproj_scaler_json);
0047     cfg.mlCombined_model_path   = make_opt<decltype(cfg.mlCombined_model_path)>(WS_Path.combined_onnx);
0048     cfg.mlCombined_scaler_json  = make_opt<decltype(cfg.mlCombined_scaler_json)>(WS_Path.combined_scaler_json);
0049
0050     // PtCalculator instance
0051     SiCaloPt::PtCalculator calcTutorial(cfg);
0052
0053     // initialize (load models and scalers for ML methods)
0054     std::string err;
0055     if (!calcTutorial.init(&err))
0056     {
0057         std::cout << "[init] failed: " << err << std::endl;
0058         return;
0059     }
0060     std::cout << "[init] OK\n";
0061
0062     // ======================== EMD Formula ==========================
0063     {
0064         SiCaloPt::InputEMD in;
0065         in.EMD_Angle  = 0.025;  // delta_phi - EM Deflection angle in rad
0066         in.EMD_Eta    = 0.00;   // EMCal Cluster eta
0067         in.EMD_Radius = 93.5;   // EMCal Cluster radius in cm
0068
0069         calcTutorial.setParCeta(0.2);   // set C_eta parameter if needed
0070         calcTutorial.setParPower(-1.0); // set Power parameter if needed
0071
0072         auto r = calcTutorial.ComputePt(SiCaloPt::Method::MethodEMD, SiCaloPt::AnyInput{in});
0073         std::cout << "[EMD-analytic] ok=" << r.ok
0074                   << "  pt=" << r.pt_reco
0075                   << "  err=\"" << r.err << "\"\n";
0076     }
0077
0078     // ======================== Eproj Formula ==========================
0079     {
0080         SiCaloPt::InputEproj in;
0081         in.Energy_Calo   = 1.8;   // EMCal Cluster energy in GeV
0082         in.Radius_Calo   = 93.5;  // EMCal Cluster radius in cm
0083         in.Z_Calo        = 0.0;   // EMCal Cluster z in cm
0084         in.Radius_vertex = 0.0;   // Vertex radius in cm
0085         in.Z_vertex      = 0.0;   // Vertex z in cm
0086
0087         auto r = calcTutorial.ComputePt(SiCaloPt::Method::MethodEproj, SiCaloPt::AnyInput{in});
0088         std::cout << "[Eproj-analytic] ok=" << r.ok
0089                   << "  pt=" << r.pt_reco
0090                   << "  err=\"" << r.err << "\"\n";
0091     }
0092
0093     // ============= ML：MLEMD (2-d input: 1/dphi_EMD, eta_track = 0 now) ===============
0094     {
0095         // 2-d input features:{ dphi_EMD, eta_track }
0096         std::vector<float> featsMLEMD = {15, 0};
0097
0098         SiCaloPt::InputMLEMD in{featsMLEMD};
0099         auto r = calcTutorial.ComputePt(SiCaloPt::Method::MethodMLEMD, SiCaloPt::AnyInput{in});
0100         std::cout << "[MLEMD-2D] ok=" << r.ok
0101                 << "  pt=" << r.pt_reco
0102                 << "  err=\"" << r.err << "\"\n";
0103     }
0104
0105     // ====== ML：MLEproj (7-d input: INTT 3/4 layer R,Z, INTT 5/4 layer R,Z, Calo R,Z,Energy) ======
0106     {
0107         // 7-d input features:{ INTT 3/4 layer R, INTT 3/4 layer Z, INTT 5/4 layer R, INTT 5/4 layer Z, Calo R, Calo Z, Calo Energy }
0108         std::vector<float> featsMLEproj = { 10.0,  5.0,   // INTT 3/4 layer
0109                                             15.0,  7.5,   // INTT 5/4 layer
0110                                             100.0, 50.0, 8.0 }; // Calo R,Z,Energy
0111
0112         SiCaloPt::InputMLEproj in{featsMLEproj};
0113         auto r = calcTutorial.ComputePt(SiCaloPt::Method::MethodMLEproj, SiCaloPt::AnyInput{in});
0114         std::cout << "[MLEproj-7D] ok=" << r.ok
0115                 << "  pt=" << r.pt_reco
0116                 << "  err=\"" << r.err << "\"\n";
0117     }
0118
0119     // ============ ML：Combined/Gate (2-d input: pt_from_MLEMD, pt_from_MLEproj) ===================
0120     {
0121         // 2-d input features:{ pt_from_MLEMD, pt_from_MLEproj }
0122         std::vector<float> featsMLCombined = {8.0, 9.5};
0123
0124         SiCaloPt::InputMLCombined in{featsMLCombined};
0125         auto r = calcTutorial.ComputePt(SiCaloPt::Method::MethodMLCombined, SiCaloPt::AnyInput{in});
0126         std::cout << "[MLCombined] ok=" << r.ok
0127                 << "  pt=" << r.pt_reco
0128                 << "  err=\"" << r.err << "\"\n";
0129     }
0130
0131     // // ======================== Batch example ========================
0132     // {
0133     //     std::vector<SiCaloPt::InputEproj> batch(1000);
0134     //     for (size_t i = 0; i < batch.size(); ++i) {
0135     //     auto& x = batch[i];
0136     //     x.Energy_Calo   = 1.0 + 0.001 * i;
0137     //     x.Radius_Calo   = 93.5;
0138     //     x.Z_Calo        = 0.0;
0139     //     x.Radius_vertex = 0.0;
0140     //     x.Z_vertex      = 0.0;
0141     //     }
0142
0143     //     TStopwatch sw; sw.Start();
0144     //     double sum_pt = 0.0; size_t ok_cnt = 0;
0145     //     for (auto& x : batch) {
0146     //     auto r = calcTutorial.ComputePt(SiCaloPt::Method::MethodEproj, SiCaloPt::AnyInput{x});
0147     //     if (r.ok) { sum_pt += r.pt_reco; ++ok_cnt; }
0148     //     }
0149     //     sw.Stop();
0150
0151     //     std::cout << "[Batch Eproj] ok=" << ok_cnt << "/" << batch.size()
0152     //             << "  avg_pt=" << (ok_cnt ? sum_pt/ok_cnt : 0.0)
0153     //             << "  time=" << sw.RealTime() << " s\n";
0154     // }
0155 }