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