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