You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

101 lines
3.0 KiB

10 months ago
  1. from parameterisations.utils.parse_regression_coef_to_array import (
  2. parse_regression_coef_to_array,
  3. )
  4. from parameterisations.utils.fit_linear_regression_model import (
  5. fit_linear_regression_model,
  6. )
  7. import uproot
  8. import argparse
  9. from pathlib import Path
  10. def parameterise_search_window(
  11. input_file: str = "data/param_data_selected_all_p.root",
  12. tree_name: str = "Selected",
  13. ) -> Path:
  14. """Function that calculates parameters for estimating the hit search window border.
  15. Args:
  16. input_file (str, optional): Defaults to "data/param_data_selected.root".
  17. tree_name (str, optional): Defaults to "Selected".
  18. Returns:
  19. Path: Path to cpp code file.
  20. """
  21. input_tree = uproot.open({input_file: tree_name})
  22. # this is an event list of dictionaries containing awkward arrays
  23. array = input_tree.arrays()
  24. array["x_straight_diff_ref"] = (
  25. array["x"] + array["tx"] * (array["z_ref"] - array["z"]) - array["x_ref"]
  26. )
  27. array["x_straight_diff_ref_abs"] = abs(array["x_straight_diff_ref"])
  28. array["inv_p_gev"] = 1000.0 / array["p"]
  29. array["pol_qop_gev"] = array["signed_rel_current"] * array["qop"] * 1000.0
  30. # now fit model for the reference plane
  31. model_ref, poly_features_ref = fit_linear_regression_model(
  32. array,
  33. target_feat="x_straight_diff_ref_abs",
  34. features=["ty", "tx", "inv_p_gev", "pol_qop_gev"],
  35. keep=[
  36. "inv_p_gev",
  37. "ty^2 inv_p_gev",
  38. "tx^2 inv_p_gev",
  39. "tx inv_p_gev pol_qop_gev",
  40. "inv_p_gev^3",
  41. "tx^3 pol_qop_gev",
  42. "tx^2 inv_p_gev^2",
  43. "tx inv_p_gev^2 pol_qop_gev",
  44. "inv_p_gev^4",
  45. "ty^2 tx^2 inv_p_gev",
  46. "ty^2 tx inv_p_gev pol_qop_gev",
  47. "ty^2 inv_p_gev^3",
  48. "tx^4 inv_p_gev",
  49. ],
  50. degree=5,
  51. fit_intercept=False,
  52. )
  53. cpp_ref = parse_regression_coef_to_array(
  54. model_ref,
  55. poly_features_ref,
  56. "momentumWindowParamsRef",
  57. )
  58. outpath = Path("parameterisations/result/search_window_params.hpp")
  59. outpath.parent.mkdir(parents=True, exist_ok=True)
  60. with open(outpath, "w") as result:
  61. result.writelines(cpp_ref)
  62. return outpath
  63. if __name__ == "__main__":
  64. parser = argparse.ArgumentParser()
  65. parser.add_argument(
  66. "--input-file",
  67. type=str,
  68. help="Path to the input file",
  69. required=False,
  70. )
  71. parser.add_argument(
  72. "--tree-name",
  73. type=str,
  74. help="Path to the input file",
  75. required=False,
  76. )
  77. args = parser.parse_args()
  78. args_dict = {arg: val for arg, val in vars(args).items() if val is not None}
  79. outfile = parameterise_search_window(**args_dict)
  80. try:
  81. import subprocess
  82. # run clang-format for nicer looking result
  83. subprocess.run(
  84. [
  85. "clang-format",
  86. "-i",
  87. f"{outfile}",
  88. ],
  89. check=True,
  90. )
  91. except:
  92. pass