tracking_eff.cpp 35 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760
  1. #include <iostream>
  2. #include <vector>
  3. #include <cmath>
  4. #include <numeric>
  5. #include <TSystem.h>
  6. #include <Math/VectorUtil.h>
  7. #include "filval.hpp"
  8. #include "root_filval.hpp"
  9. #include "common.hpp"
  10. #include "analysis/TrackingNtupleObjs.hpp"
  11. typedef ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<float>> Vec4;
  12. SeedCollection seeds;
  13. SimTrackCollection sim_tracks;
  14. SimVertexCollection sim_vertices;
  15. TrackCollection gsf_tracks;
  16. SuperClusterCollection scls;
  17. void register_objects(TrackingDataSet& tds){
  18. seeds.init(&tds);
  19. sim_tracks.init(&tds);
  20. sim_vertices.init(&tds);
  21. gsf_tracks.init(&tds);
  22. scls.init(&tds);
  23. }
  24. struct {
  25. Value<float>* x;
  26. Value<float>* y;
  27. Value<float>* z;
  28. Value<float>* sigma_x;
  29. Value<float>* sigma_y;
  30. Value<float>* sigma_z;
  31. } bs;
  32. template <typename T>
  33. struct KinColl {
  34. T* pt;
  35. T* eta;
  36. T* phi;
  37. };
  38. bool in_lum_region(const SimVertex& vertex) {
  39. /* Here approximate the luminous region as a cylinder with length 5*\sigma_z and radius
  40. * sqrt(\simga_x^2 + \sigma_y^2) centered as beamspot position.
  41. */
  42. float dx = vertex.x() - bs.x->get();
  43. float dy = vertex.y() - bs.y->get();
  44. float dz = vertex.z() - bs.z->get();
  45. auto radius = static_cast<float>(5 * sqrt(pow(bs.sigma_x->get(), 2) + pow(bs.sigma_y->get(), 2)));
  46. float half_len = 5*bs.sigma_z->get();
  47. return sqrt(dx*dx + dy*dy) < radius and abs(dz) < half_len;
  48. };
  49. bool is_good_sim(const SimTrack& sim_track) {
  50. const auto& vertex = sim_vertices[sim_track.parentVtxIdx()];
  51. return abs(sim_track.pdgId()) == 11 and // electrons
  52. in_lum_region(vertex) and // from luminous region
  53. abs(sim_track.eta())<3.0; // inside ECAL acceptance
  54. };
  55. //bool is_good_fake(const SimTrack& sim_track) {
  56. // const auto& vertex = sim_vertices[sim_track.parentVtxIdx()];
  57. // return abs(sim_track.pdgId()) != 11 and sim_track.q() != 0 and in_lum_region(vertex);
  58. //};
  59. bool is_good_seed(const Seed& seed, float hoe_cut) {
  60. return seed.isECALDriven() and seed.hoe() < hoe_cut;
  61. }
  62. float reco_energy_rel_err(const SimTrack& sim_track, const Seed& seed) {
  63. return (sim_track.pt() - seed.Et()) / sim_track.pt() ;
  64. }
  65. float reco_energy_rel_err(const SimTrack& sim_track, const Track& track) {
  66. return (sim_track.pt() - track.pt()) / sim_track.pt() ;
  67. }
  68. template<typename TrackOrSeed>
  69. bool reco_energy_consistent(const SimTrack& sim_track, const TrackOrSeed& track, float consistency_cut=0.1) {
  70. return fabs(reco_energy_rel_err(sim_track, track)) < consistency_cut;
  71. }
  72. Vec4 scl_p4(const SuperCluster& scl) {
  73. return Vec4(scl.px(), scl.py(), scl.pz(), scl.e());
  74. }
  75. double deltaR(float scl_eta, float scl_phi, const SimTrack& sim_track) {
  76. return sqrt(pow(sim_track.eta() - scl_eta, 2.0) +
  77. pow(sim_track.phi() - scl_phi, 2.0));
  78. }
  79. double deltaR(const SimTrack& sim_track, const Track& gsf_track) {
  80. return sqrt(pow(sim_track.eta() - gsf_track.eta(), 2.0) +
  81. pow(sim_track.phi() - gsf_track.phi(), 2.0));
  82. }
  83. auto get_match_stats(bool dRMatch) {
  84. int nSimTrack = 0;
  85. int nGSFTrack = 0;
  86. int nMatched = 0; // 1-to-1
  87. int nMerged = 0; // n-to-1
  88. int nLost = 0; // 1-to-0
  89. int nSplit = 0; // 1-to-n
  90. int nFaked = 0; // 0-to-1
  91. int nFlipped = 0;
  92. vector<float> matched_dR;
  93. vector<float> matched_dpT;
  94. std::map<int, std::vector<int>> sim_matches;
  95. std::map<int, std::vector<int>> gsf_matches;
  96. for(const auto& sim_track : sim_tracks)
  97. if (is_good_sim(sim_track)) sim_matches[sim_track.idx] = {};
  98. for(const auto& gsf_track : gsf_tracks)
  99. gsf_matches[gsf_track.idx] = {};
  100. nSimTrack = sim_matches.size();
  101. nGSFTrack = gsf_matches.size();
  102. for (const auto& sim_track : sim_tracks) {
  103. if (!is_good_sim(sim_track)) continue;
  104. if (dRMatch) {
  105. for (const auto& gsf_track : gsf_tracks) {
  106. float dR = deltaR(sim_track, gsf_track);
  107. if (dR < 0.2) {
  108. sim_matches[sim_track.idx].push_back(gsf_track.idx);
  109. gsf_matches[gsf_track.idx].push_back(sim_track.idx);
  110. matched_dR.push_back(dR);
  111. matched_dpT.push_back((sim_track.pt() - gsf_track.pt())/sim_track.pt());
  112. if (gsf_track.q() != sim_track.q()) nFlipped++;
  113. }
  114. }
  115. } else {
  116. for (const auto& gsfIdx : sim_track.trkIdx()) {
  117. const Track& gsf_track = gsf_tracks[gsfIdx];
  118. float dR = deltaR(sim_track, gsf_track);
  119. sim_matches[sim_track.idx].push_back(gsf_track.idx);
  120. gsf_matches[gsf_track.idx].push_back(sim_track.idx);
  121. matched_dR.push_back(dR);
  122. matched_dpT.push_back((sim_track.pt() - gsf_track.pt())/sim_track.pt());
  123. if (gsf_track.q() != sim_track.q()) nFlipped++;
  124. }
  125. }
  126. }
  127. for (const auto& p : sim_matches) {
  128. const auto& matches = p.second;
  129. if (matches.size() == 0) nLost++;
  130. if (matches.size() == 1 and gsf_matches[matches[0]].size() == 1) nMatched++;
  131. if (matches.size() > 1) nSplit++;
  132. }
  133. for (const auto& p : gsf_matches) {
  134. const auto& matches = p.second;
  135. if (matches.size() == 0) nFaked++;
  136. if (matches.size() > 1) nMerged++;
  137. }
  138. return std::make_tuple(nSimTrack, nGSFTrack, nMatched, nMerged, nLost, nSplit, nFaked, nFlipped, matched_dR, matched_dpT);
  139. }
  140. void run(){
  141. using namespace std;
  142. using namespace fv;
  143. using namespace fv_root;
  144. using namespace fv_util;
  145. auto file_list = the_config->get_source_files();
  146. string output_filename = the_config->get_output_filename();
  147. TrackingDataSet tds(output_filename, file_list, "trackingNtuple/tree");
  148. register_objects(tds);
  149. auto hoe_cut = the_config->get("hoe_cut").as<float>(999);
  150. bs = {
  151. tds.track_branch<float>("bsp_x"),
  152. tds.track_branch<float>("bsp_y"),
  153. tds.track_branch<float>("bsp_z"),
  154. tds.track_branch<float>("bsp_sigmax"),
  155. tds.track_branch<float>("bsp_sigmay"),
  156. tds.track_branch<float>("bsp_sigmaz")
  157. };
  158. enum TMType {
  159. NoMatch = 0,
  160. SeedMatched = 1,
  161. TrackMatched = 2,
  162. SeedAndTrackMatched = 3,
  163. };
  164. std::map<std::tuple<int, int>, ContainerTH2<float>*> residuals_dRz;
  165. std::map<std::tuple<int, int>, ContainerTH2<float>*> residuals_dPhi;
  166. std::map<std::tuple<int, int, int>, ContainerTH2<float>*> BPIX_residuals_dRz;
  167. std::map<std::tuple<int, int, int>, ContainerTH2<float>*> BPIX_residuals_dPhi;
  168. std::map<std::tuple<int, int, int>, ContainerTH2<float>*> FPIX_residuals_dRz;
  169. std::map<std::tuple<int, int, int>, ContainerTH2<float>*> FPIX_residuals_dPhi;
  170. stringstream name;
  171. auto set_name = [&name](const std::string& var, const std::string& region,
  172. const int& layer, const int& hit, const int& tm_type) {
  173. name.str("");
  174. name << var << "_" << region;
  175. if (layer)
  176. name << "_L" << layer;
  177. name << "_H" << hit << "_v_Et";
  178. switch (tm_type) {
  179. case NoMatch:
  180. name << "_NoMatch";
  181. break;
  182. case SeedMatched:
  183. name << "_SeedMatched";
  184. break;
  185. case TrackMatched:
  186. name << "_TrackMatched";
  187. break;
  188. case SeedAndTrackMatched:
  189. name << "_SeedAndTrackMatched";
  190. break;
  191. default:
  192. break;
  193. }
  194. };
  195. THParams hist_params;
  196. for (int hit=1; hit<=3; hit++) {
  197. for (int tm_type = NoMatch; tm_type <= SeedAndTrackMatched; tm_type++) {
  198. for (int layer=1; layer<=4; layer++) {
  199. hist_params = (hit == 1) ? THParams::lookup("dRz_v_Et") : THParams::lookup("dRz_v_Et_outer_hits");
  200. set_name("dRz", "BPIX", layer, hit, tm_type);
  201. BPIX_residuals_dRz[make_tuple(layer, hit, tm_type)] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
  202. set_name("dRz", "FPIX", layer, hit, tm_type);
  203. FPIX_residuals_dRz[make_tuple(layer, hit, tm_type)] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
  204. hist_params = (hit == 1) ? THParams::lookup("dPhi_v_Et") : THParams::lookup("dPhi_v_Et_outer_hits");
  205. set_name("dPhi", "BPIX", layer, hit, tm_type);
  206. BPIX_residuals_dPhi[make_tuple(layer, hit, tm_type)] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
  207. set_name("dPhi", "FPIX", layer, hit, tm_type);
  208. FPIX_residuals_dPhi[make_tuple(layer, hit, tm_type)] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
  209. }
  210. hist_params = (hit == 1) ? THParams::lookup("dRz_v_Et") : THParams::lookup("dRz_v_Et_outer_hits");
  211. set_name("dRz", "ALL", 0, hit, tm_type);
  212. residuals_dRz[std::make_tuple(hit, tm_type)] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
  213. hist_params = (hit == 1) ? THParams::lookup("dPhi_v_Et") : THParams::lookup("dPhi_v_Et_outer_hits");
  214. set_name("dPhi", "ALL", 0, hit, tm_type);
  215. residuals_dPhi[make_tuple(hit, tm_type)] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
  216. }
  217. }
  218. std::map<int,vector<float>> eta_regions = {{1, {0.0, 1.1, 1.8, 3.0}},
  219. {2, {0.0, 1.4, 2.3, 3.0}},
  220. {3, {0.0, 1.0, 2.0, 3.0}}};
  221. auto get_region = [&eta_regions](int hit, float eta) {
  222. auto regions = eta_regions[hit];
  223. if (regions.size() <= 1) return 1;
  224. float abseta = abs(eta);
  225. for (int r_idx=1; r_idx < regions.size(); r_idx++) {
  226. if (abseta > regions[r_idx-1] and abseta < regions[r_idx]) {
  227. return r_idx;
  228. }
  229. }
  230. return static_cast<int>(regions.size());
  231. };
  232. std::map<int, ContainerTH2<float>*> dPhi_residuals_v_eta;
  233. std::map<int, ContainerTH2<float>*> dRz_residuals_v_eta;
  234. std::map<std::pair<int, int>, ContainerTH2<float>*> dPhi_residuals_v_region;
  235. std::map<std::pair<int, int>, ContainerTH2<float>*> dRz_residuals_v_region;
  236. for (int hit=1; hit<=3; hit++) {
  237. name.str("");
  238. name << "dPhi_residuals_v_eta_H" << hit;
  239. hist_params = (hit==1) ? THParams::lookup("dPhi_v_eta") : THParams::lookup("dPhi_v_eta_outer_hits");
  240. dPhi_residuals_v_eta[hit] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
  241. name.str("");
  242. name << "dRz_residuals_v_eta_H" << hit;
  243. hist_params = (hit==1) ? THParams::lookup("dRz_v_eta") : THParams::lookup("dRz_v_eta_outer_hits");
  244. dRz_residuals_v_eta[hit] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
  245. for (int region=1; region<=eta_regions[hit].size()-1; region++){
  246. hist_params = (hit==1) ? THParams::lookup("dRz_v_Et") : THParams::lookup("dRz_v_Et_outer_hits");
  247. name.str("");
  248. name << "dRz_residuals_H" << hit << "_R" << region;
  249. dRz_residuals_v_region[{hit, region}] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
  250. hist_params = (hit==1) ? THParams::lookup("dPhi_v_Et") : THParams::lookup("dPhi_v_Et_outer_hits");
  251. name.str("");
  252. name << "dPhi_residuals_H" << hit << "_R" << region;
  253. dPhi_residuals_v_region[{hit, region}] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
  254. }
  255. }
  256. KinColl<ContainerTH1<float>> good_sim_kinems = {
  257. tds.register_container<ContainerTH1<float>>("good_sim_v_pt", THParams::lookup("pt")),
  258. tds.register_container<ContainerTH1<float>>("good_sim_v_eta", THParams::lookup("eta")),
  259. tds.register_container<ContainerTH1<float>>("good_sim_v_phi", THParams::lookup("phi"))};
  260. KinColl<ContainerTH1<float>> gsf_kinems = {
  261. tds.register_container<ContainerTH1<float>>("gsf_track_v_pt", THParams::lookup("pt")),
  262. tds.register_container<ContainerTH1<float>>("gsf_track_v_eta", THParams::lookup("eta")),
  263. tds.register_container<ContainerTH1<float>>("gsf_track_v_phi", THParams::lookup("phi"))};
  264. KinColl<ContainerTH1<float>> seed_kinems = {
  265. tds.register_container<ContainerTH1<float>>("seed_v_pt", THParams::lookup("pt")),
  266. tds.register_container<ContainerTH1<float>>("seed_v_eta", THParams::lookup("eta")),
  267. tds.register_container<ContainerTH1<float>>("seed_v_phi", THParams::lookup("phi"))};
  268. KinColl<ContainerTH1<float>> scl_kinems = {
  269. tds.register_container<ContainerTH1<float>>("scl_v_pt", THParams::lookup("pt")),
  270. tds.register_container<ContainerTH1<float>>("scl_v_eta", THParams::lookup("eta")),
  271. tds.register_container<ContainerTH1<float>>("scl_v_phi", THParams::lookup("phi"))};
  272. KinColl<EfficiencyContainer<float>> seed_eff = {
  273. tds.register_container<EfficiencyContainer<float>>("seed_eff_v_pt", THParams::lookup("pt")),
  274. tds.register_container<EfficiencyContainer<float>>("seed_eff_v_eta", THParams::lookup("eta")),
  275. tds.register_container<EfficiencyContainer<float>>("seed_eff_v_phi", THParams::lookup("phi"))};
  276. KinColl<EfficiencyContainer<float>> seed_pur = {
  277. tds.register_container<EfficiencyContainer<float>>("seed_pur_v_pt", THParams::lookup("pt")),
  278. tds.register_container<EfficiencyContainer<float>>("seed_pur_v_eta", THParams::lookup("eta")),
  279. tds.register_container<EfficiencyContainer<float>>("seed_pur_v_phi", THParams::lookup("phi"))};
  280. KinColl<EfficiencyContainer<float>> tracking_eff = {
  281. tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_pt", THParams::lookup("pt")),
  282. tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_eta", THParams::lookup("eta")),
  283. tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_phi", THParams::lookup("phi"))};
  284. KinColl<EfficiencyContainer<float>> tracking_pur = {
  285. tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_pt", THParams::lookup("pt")),
  286. tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_eta", THParams::lookup("eta")),
  287. tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_phi", THParams::lookup("phi"))};
  288. KinColl<EfficiencyContainer<float>> tracking_eff_dR = {
  289. tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_pt_dR", THParams::lookup("pt")),
  290. tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_eta_dR", THParams::lookup("eta")),
  291. tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_phi_dR", THParams::lookup("phi"))};
  292. KinColl<EfficiencyContainer<float>> tracking_pur_dR = {
  293. tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_pt_dR", THParams::lookup("pt")),
  294. tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_eta_dR", THParams::lookup("eta")),
  295. tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_phi_dR", THParams::lookup("phi"))};
  296. KinColl<EfficiencyContainer<float>> partial_fake_rate = {
  297. tds.register_container<EfficiencyContainer<float>>("partial_fake_rate_v_pt", THParams::lookup("pt")),
  298. tds.register_container<EfficiencyContainer<float>>("partial_fake_rate_v_eta", THParams::lookup("eta")),
  299. tds.register_container<EfficiencyContainer<float>>("partial_fake_rate_v_phi", THParams::lookup("phi"))};
  300. KinColl<EfficiencyContainer<float>> full_fake_rate = {
  301. tds.register_container<EfficiencyContainer<float>>("full_fake_rate_v_pt", THParams::lookup("pt")),
  302. tds.register_container<EfficiencyContainer<float>>("full_fake_rate_v_eta", THParams::lookup("eta")),
  303. tds.register_container<EfficiencyContainer<float>>("full_fake_rate_v_phi", THParams::lookup("phi"))};
  304. KinColl<EfficiencyContainer<float>> clean_fake_rate = {
  305. tds.register_container<EfficiencyContainer<float>>("clean_fake_rate_v_pt", THParams::lookup("pt")),
  306. tds.register_container<EfficiencyContainer<float>>("clean_fake_rate_v_eta", THParams::lookup("eta")),
  307. tds.register_container<EfficiencyContainer<float>>("clean_fake_rate_v_phi", THParams::lookup("phi"))};
  308. auto& hit_vs_layer_barrel = *tds.register_container<ContainerTH2<int>>("hit_vs_layer_barrel", THParams::lookup("hit_vs_layer"));
  309. auto& hit_vs_layer_forward = *tds.register_container<ContainerTH2<int>>("hit_vs_layer_forward", THParams::lookup("hit_vs_layer"));
  310. auto& gsf_tracks_nmatch_sim_tracks = *tds.register_container<ContainerTH1<size_t>>("gsf_tracks_nmatch_sim_tracks ", THParams::lookup("gsf_tracks_nmatch_sim_tracks"));
  311. auto& ecal_energy_resolution = *tds.register_container<ContainerTH1<float>>("ecal_energy_resolution ", THParams::lookup("ecal_energy_resolution"));
  312. auto& n_seeds = *tds.register_container<ContainerTH1<size_t>>("n_seeds", THParams::lookup("n_seeds"));
  313. auto& n_good_seeds = *tds.register_container<ContainerTH1<size_t>>("n_good_seeds", THParams::lookup("n_seeds"));
  314. auto& n_good_sim = *tds.register_container<ContainerTH1<int>>("n_good_sim", THParams::lookup("n_seeds"));
  315. auto& n_gsf_tracks = *tds.register_container<ContainerTH1<int>>("n_gsf_track", THParams::lookup("n_seeds"));
  316. auto& n_matched = *tds.register_container<ContainerTH1<int>>("n_matched", THParams::lookup("n_seeds"));
  317. auto& n_merged = *tds.register_container<ContainerTH1<int>>("n_merged", THParams::lookup("n_seeds"));
  318. auto& n_lost = *tds.register_container<ContainerTH1<int>>("n_lost", THParams::lookup("n_seeds"));
  319. auto& n_split = *tds.register_container<ContainerTH1<int>>("n_split", THParams::lookup("n_seeds"));
  320. auto& n_faked = *tds.register_container<ContainerTH1<int>>("n_faked", THParams::lookup("n_seeds"));
  321. auto& n_flipped = *tds.register_container<ContainerTH1<int>>("n_flipped", THParams::lookup("n_seeds"));
  322. auto& matched_dR = *tds.register_container<ContainerTH1<float>>("matched_dR", THParams::lookup("matched_dR"));
  323. auto& matched_dpT = *tds.register_container<ContainerTH1<float>>("matched_dpT", THParams::lookup("matched_dpT"));
  324. auto& n_matched_dR = *tds.register_container<ContainerTH1<int>>("n_matched_dR", THParams::lookup("n_seeds"));
  325. auto& n_merged_dR = *tds.register_container<ContainerTH1<int>>("n_merged_dR", THParams::lookup("n_seeds"));
  326. auto& n_lost_dR = *tds.register_container<ContainerTH1<int>>("n_lost_dR", THParams::lookup("n_seeds"));
  327. auto& n_split_dR = *tds.register_container<ContainerTH1<int>>("n_split_dR", THParams::lookup("n_seeds"));
  328. auto& n_faked_dR = *tds.register_container<ContainerTH1<int>>("n_faked_dR", THParams::lookup("n_seeds"));
  329. auto& n_flipped_dR = *tds.register_container<ContainerTH1<int>>("n_flipped_dR", THParams::lookup("n_seeds"));
  330. auto& matched_dR_dR = *tds.register_container<ContainerTH1<float>>("matched_dR_dR", THParams::lookup("matched_dR"));
  331. auto& matched_dpT_dR = *tds.register_container<ContainerTH1<float>>("matched_dpT_dR", THParams::lookup("matched_dpT"));
  332. auto& sim_pt = *tds.register_container<ContainerTH1<float>>("sim_pt", THParams::lookup("pt_fine"));
  333. auto& sim_eta = *tds.register_container<ContainerTH1<float>>("sim_eta", THParams::lookup("eta_fine"));
  334. auto& sim_phi = *tds.register_container<ContainerTH1<float>>("sim_phi", THParams::lookup("phi_fine"));
  335. auto& reco_pt = *tds.register_container<ContainerTH1<float>>("reco_pt", THParams::lookup("pt_fine"));
  336. auto& reco_eta = *tds.register_container<ContainerTH1<float>>("reco_eta", THParams::lookup("eta_fine"));
  337. auto& reco_phi = *tds.register_container<ContainerTH1<float>>("reco_phi", THParams::lookup("phi_fine"));
  338. while (tds.next()) {
  339. n_seeds.fill(seeds.size());
  340. size_t _n_good_seeds = 0;
  341. for (const auto& seed : seeds) {
  342. if (is_good_seed(seed, hoe_cut)) {
  343. _n_good_seeds++;
  344. seed_kinems.pt->fill(seed.pt());
  345. seed_kinems.eta->fill(seed.eta());
  346. seed_kinems.phi->fill(seed.phi());
  347. }
  348. }
  349. n_good_seeds.fill(_n_good_seeds);
  350. for (const auto& sim_track : sim_tracks) {
  351. if (is_good_sim(sim_track)) {
  352. good_sim_kinems.pt->fill(sim_track.pt());
  353. good_sim_kinems.eta->fill(sim_track.eta());
  354. good_sim_kinems.phi->fill(sim_track.phi());
  355. }
  356. }
  357. for (const auto& gsf_track : gsf_tracks) {
  358. gsf_kinems.pt->fill(gsf_track.pt());
  359. gsf_kinems.eta->fill(gsf_track.eta());
  360. gsf_kinems.phi->fill(gsf_track.phi());
  361. }
  362. for (const auto& scl : scls) {
  363. scl_kinems.pt->fill(hypot(scl.px(), scl.py()));
  364. scl_kinems.eta->fill(pseudorapidity(scl.x(), scl.y(), scl.z()));
  365. scl_kinems.phi->fill(atan2(scl.py(), scl.px()));
  366. }
  367. auto match_stats = get_match_stats(false);
  368. n_good_sim.fill(std::get<0>(match_stats));
  369. n_gsf_tracks.fill(std::get<1>(match_stats));
  370. n_matched.fill(std::get<2>(match_stats));
  371. n_merged.fill(std::get<3>(match_stats));
  372. n_lost.fill(std::get<4>(match_stats));
  373. n_split.fill(std::get<5>(match_stats));
  374. n_faked.fill(std::get<6>(match_stats));
  375. n_flipped.fill(std::get<7>(match_stats));
  376. for(float dR : std::get<8>(match_stats)) matched_dR.fill(dR);
  377. for(float dpT : std::get<9>(match_stats)) matched_dpT.fill(dpT);
  378. match_stats = get_match_stats(true);
  379. n_matched_dR.fill(std::get<2>(match_stats));
  380. n_merged_dR.fill(std::get<3>(match_stats));
  381. n_lost_dR.fill(std::get<4>(match_stats));
  382. n_split_dR.fill(std::get<5>(match_stats));
  383. n_faked_dR.fill(std::get<6>(match_stats));
  384. n_flipped_dR.fill(std::get<7>(match_stats));
  385. for(float dR : std::get<8>(match_stats)) matched_dR_dR.fill(dR);
  386. for(float dpT : std::get<9>(match_stats)) matched_dpT_dR.fill(dpT);
  387. for (const auto& sim_track : sim_tracks) {
  388. if (!is_good_sim(sim_track)) continue;
  389. sim_pt.fill(sim_track.pt());
  390. sim_eta.fill(sim_track.eta());
  391. sim_phi.fill(sim_track.phi());
  392. }
  393. for (const auto& reco_track : gsf_tracks) {
  394. reco_pt.fill(reco_track.pt());
  395. reco_eta.fill(reco_track.eta());
  396. reco_phi.fill(reco_track.phi());
  397. }
  398. for (const auto& sim_track : sim_tracks) {
  399. if (!is_good_sim(sim_track)) continue;
  400. if (seeds.size() == 0) continue;
  401. for (const auto &seed_idx : sim_track.seedIdx()) {
  402. const Seed& seed = seeds[seed_idx];
  403. if (!is_good_seed(seed, hoe_cut)) continue;
  404. ecal_energy_resolution.fill(reco_energy_rel_err(sim_track, seed));
  405. }
  406. }
  407. // Seeding Efficiency
  408. for (const auto& sim_track : sim_tracks) {
  409. if (!is_good_sim(sim_track)) continue;
  410. if (seeds.size() == 0) continue;
  411. bool matched = false;
  412. for (const auto& seed_idx : sim_track.seedIdx()) {
  413. const Seed& seed = seeds[seed_idx];
  414. if (is_good_seed(seed, hoe_cut)) {
  415. matched=true;
  416. break;
  417. }
  418. }
  419. if (abs(sim_track.eta()) < 2.4)
  420. seed_eff.pt->fill(sim_track.pt(), matched);
  421. if (sim_track.pt() > 20.0)
  422. seed_eff.eta->fill(sim_track.eta(), matched);
  423. if (abs(sim_track.eta()) < 2.4 and sim_track.pt() > 20)
  424. seed_eff.phi->fill(sim_track.phi(), matched);
  425. }
  426. // Tracking Efficiency
  427. for (const auto& sim_track : sim_tracks) {
  428. if (!is_good_sim(sim_track)) continue;
  429. if (gsf_tracks.size() == 0) continue;
  430. bool matched = false;
  431. for (const auto& track_idx : sim_track.trkIdx()) {
  432. const Seed& seed = seeds[gsf_tracks[track_idx].seedIdx()];
  433. if (is_good_seed(seed, hoe_cut)) {
  434. matched=true;
  435. break;
  436. }
  437. }
  438. if (abs(sim_track.eta()) < 2.4)
  439. tracking_eff.pt->fill(sim_track.pt(), matched);
  440. if (sim_track.pt() > 20.0)
  441. tracking_eff.eta->fill(sim_track.eta(), matched);
  442. if (abs(sim_track.eta()) < 2.4 and sim_track.pt() > 20)
  443. tracking_eff.phi->fill(sim_track.phi(), matched);
  444. // dR-matching
  445. matched = false;
  446. for (const auto& gsf_track : gsf_tracks) {
  447. const Seed& seed = seeds[gsf_track.seedIdx()];
  448. double dR = deltaR(sim_track, gsf_track);
  449. if (is_good_seed(seed, hoe_cut) and dR < 0.2) {
  450. matched = true;
  451. break;
  452. }
  453. }
  454. if (abs(sim_track.eta()) < 2.4)
  455. tracking_eff_dR.pt->fill(sim_track.pt(), matched);
  456. if (sim_track.pt() > 20.0)
  457. tracking_eff_dR.eta->fill(sim_track.eta(), matched);
  458. if (abs(sim_track.eta()) < 2.4 and sim_track.pt() > 20)
  459. tracking_eff_dR.phi->fill(sim_track.phi(), matched);
  460. }
  461. // Seeding Purity
  462. for (const auto& seed : seeds) {
  463. if (!is_good_seed(seed, hoe_cut)) continue;
  464. bool match = false;
  465. for (const auto& sim_track_idx : seed.simTrkIdx()) {
  466. if(is_good_sim(sim_tracks[sim_track_idx])) {
  467. match = true;
  468. break;
  469. }
  470. }
  471. if (abs(seed.eta()) < 2.4)
  472. seed_pur.pt->fill(seed.pt(), match);
  473. if (seed.pt() > 20)
  474. seed_pur.eta->fill(seed.eta(), match);
  475. if (abs(seed.eta()) < 2.4 and seed.pt() > 20)
  476. seed_pur.phi->fill(seed.phi(), match);
  477. }
  478. // Tracking Purity
  479. for (const auto& gsf_track : gsf_tracks) {
  480. gsf_tracks_nmatch_sim_tracks.fill(gsf_track.simTrkIdx().size());
  481. const Seed& seed = seeds[gsf_track.seedIdx()];
  482. if (!is_good_seed(seed, hoe_cut)) continue;
  483. bool matched = false;
  484. for (const auto& sim_track_idx : gsf_track.simTrkIdx()) {
  485. if (is_good_sim(sim_tracks[sim_track_idx])) {
  486. matched = true;
  487. break;
  488. }
  489. }
  490. if (abs(gsf_track.eta()) < 2.4)
  491. tracking_pur.pt->fill(gsf_track.pt(), matched);
  492. if (gsf_track.pt() > 20)
  493. tracking_pur.eta->fill(gsf_track.eta(), matched);
  494. if (abs(gsf_track.eta()) < 2.4 and gsf_track.pt() > 20)
  495. tracking_pur.phi->fill(gsf_track.phi(), matched);
  496. // dR-matching
  497. matched = false;
  498. for (const auto& sim_track : sim_tracks) {
  499. double dR = deltaR(sim_track, gsf_track);
  500. if (is_good_sim(sim_track) and dR < 0.2) {
  501. matched = true;
  502. break;
  503. }
  504. }
  505. if (abs(gsf_track.eta()) < 2.4)
  506. tracking_pur_dR.pt->fill(gsf_track.pt(), matched);
  507. if (gsf_track.pt() > 20.0)
  508. tracking_pur_dR.eta->fill(gsf_track.eta(), matched);
  509. if (abs(gsf_track.eta()) < 2.4 and gsf_track.pt() > 20)
  510. tracking_pur_dR.phi->fill(gsf_track.phi(), matched);
  511. }
  512. // Fake-Rate
  513. std::map<size_t, std::vector<int>> scl2trk; // Maps superclusters to tm-status of matched gsf-tracks
  514. for (const auto &gsf_track : gsf_tracks) {
  515. bool gsf_truth_matched = false;
  516. for(const auto& sim_track_idx : gsf_track.simTrkIdx()) {
  517. if (is_good_sim(sim_tracks[sim_track_idx])) {
  518. gsf_truth_matched = true;
  519. break;
  520. }
  521. }
  522. int scl_idx = seeds[gsf_track.seedIdx()].sclIdx();
  523. scl2trk[scl_idx].push_back(gsf_truth_matched);
  524. }
  525. std::set<size_t> tm_scls; // Set of super-clusters with well matched sim-electrons
  526. for (const auto &scl : scls) {
  527. Vec4 p4 = scl_p4(scl);
  528. for(const auto& sim_track : sim_tracks) {
  529. if (!is_good_sim(sim_track)) continue;
  530. if (deltaR(p4.eta(), p4.phi(), sim_track) > 0.2) continue;
  531. if (((p4.Et() - sim_track.pt()) / p4.Et()) > 0.1) continue;
  532. tm_scls.insert(scl.idx);
  533. break;
  534. }
  535. }
  536. /* cout << "scls: "; */
  537. /* for (const auto& scl: scls) cout << scl.idx << ", "; */
  538. /* cout << endl << "tm_scls: "; */
  539. /* for (const auto& idx: tm_scls) cout << idx << ", "; */
  540. /* cout << endl; */
  541. for (const auto &scl : scls) {
  542. if (scl.hoe() > hoe_cut) continue;
  543. float scl_pt = hypot(scl.px(), scl.py());
  544. float scl_eta = pseudorapidity(scl.x(), scl.y(), scl.z());
  545. float scl_phi = atan2(scl.py(), scl.px());
  546. int ntracks = scl2trk[scl.idx].size();
  547. int ntmtracks = std::accumulate(scl2trk[scl.idx].begin(), scl2trk[scl.idx].end(), 0);
  548. bool partial_fake = ntmtracks > 0 and ntracks > ntmtracks;
  549. bool full_fake = ntmtracks == 0 and ntracks > 0;
  550. /* cout << "ntracks: " << ntracks << " "; */
  551. /* cout << "count: " << tm_scls.count(scl.idx); */
  552. if (ntracks > 0) {
  553. partial_fake_rate.pt->fill(scl_pt, partial_fake);
  554. partial_fake_rate.eta->fill(scl_eta, partial_fake);
  555. partial_fake_rate.phi->fill(scl_phi, partial_fake);
  556. full_fake_rate.pt->fill(scl_pt, full_fake);
  557. full_fake_rate.eta->fill(scl_eta, full_fake);
  558. full_fake_rate.phi->fill(scl_phi, full_fake);
  559. if (tm_scls.count(scl.idx) == 0) {
  560. /* cout << "Filling clean_fake_rate" << endl; */
  561. clean_fake_rate.pt->fill(scl_pt, full_fake);
  562. clean_fake_rate.eta->fill(scl_eta, full_fake);
  563. clean_fake_rate.phi->fill(scl_phi, full_fake);
  564. }
  565. }
  566. /* cout << endl; */
  567. }
  568. // Hit Residuals
  569. for (const auto& seed : seeds) {
  570. if (seed.trkIdx() < 0) continue; // require that seed produced gsf-track
  571. if (!is_good_seed(seed, hoe_cut)) continue;
  572. bool is_seed_sim_matched = false;
  573. for (const auto& sim_track_idx : seed.simTrkIdx()) {
  574. const SimTrack& sim_track = sim_tracks[sim_track_idx];
  575. if (is_good_sim(sim_track) and reco_energy_consistent(sim_track, seed)) {
  576. is_seed_sim_matched = true;
  577. break;
  578. }
  579. }
  580. bool is_track_sim_matched = false;
  581. const auto the_track = gsf_tracks[seed.trkIdx()];
  582. for (const auto& sim_track_idx : the_track.simTrkIdx()) {
  583. const SimTrack& sim_track = sim_tracks[sim_track_idx];
  584. if (is_good_sim(sim_track) and reco_energy_consistent(sim_track, the_track)) {
  585. is_track_sim_matched = true;
  586. break;
  587. }
  588. }
  589. std::vector<TMType> tm_types;
  590. if (is_seed_sim_matched and is_track_sim_matched) {
  591. tm_types = {SeedAndTrackMatched, SeedMatched, TrackMatched};
  592. } else if (is_seed_sim_matched) {
  593. tm_types = {SeedMatched};
  594. } else if (is_track_sim_matched) {
  595. tm_types = {TrackMatched};
  596. } else {
  597. tm_types = {NoMatch};
  598. }
  599. vector<int> hits_valid;
  600. vector<float> hits_dRz;
  601. vector<float> hits_dPhi;
  602. if (the_track.q() == 1) {
  603. hits_valid = seed.isValidPos();
  604. hits_dRz = seed.dRZPos();
  605. hits_dPhi = seed.dPhiPos();
  606. } else {
  607. hits_valid = seed.isValidNeg();
  608. hits_dRz = seed.dRZNeg();
  609. hits_dPhi = seed.dPhiNeg();
  610. }
  611. const vector<int>& hits_isBarrel = seed.isBarrel();
  612. const vector<int>& hits_layerOrDiskNr = seed.layerOrDiskNr();
  613. for (const auto& tm_type : tm_types) {
  614. size_t nHits = hits_valid.size();
  615. for (size_t hit_idx = 0; hit_idx < nHits; hit_idx++) {
  616. if (!hits_valid[hit_idx]) continue;
  617. bool isBarrel = hits_isBarrel[hit_idx] == 1;
  618. int layerOrDiskNr = hits_layerOrDiskNr[hit_idx];
  619. float dRz = abs(hits_dRz[hit_idx]);
  620. float dPhi = abs(hits_dPhi[hit_idx]);
  621. int eta_region = get_region(static_cast<int>(hit_idx + 1), seed.eta());
  622. if (is_seed_sim_matched) {
  623. dRz_residuals_v_eta[hit_idx + 1]->fill(dRz, abs(seed.eta()));
  624. dPhi_residuals_v_eta[hit_idx + 1]->fill(dPhi, abs(seed.eta()));
  625. dRz_residuals_v_region[make_pair(hit_idx + 1, eta_region)]->fill(dRz, seed.Et());
  626. dPhi_residuals_v_region[make_pair(hit_idx + 1, eta_region)]->fill(dPhi, seed.Et());
  627. }
  628. residuals_dRz[make_tuple(hit_idx + 1, tm_type)]->fill(dRz, seed.Et());
  629. residuals_dPhi[make_tuple(hit_idx + 1, tm_type)]->fill(dPhi, seed.Et());
  630. if (layerOrDiskNr == -1) continue; // subDet is not set w/ old seeding
  631. if (isBarrel)
  632. hit_vs_layer_barrel.fill(layerOrDiskNr, static_cast<const int &>(hit_idx + 1));
  633. else
  634. hit_vs_layer_forward.fill(layerOrDiskNr, static_cast<const int &>(hit_idx + 1));
  635. if (isBarrel) {
  636. BPIX_residuals_dRz[make_tuple(layerOrDiskNr, hit_idx + 1, tm_type)]->fill(dRz, seed.Et());
  637. BPIX_residuals_dPhi[make_tuple(layerOrDiskNr, hit_idx + 1, tm_type)]->fill(dPhi, seed.Et());
  638. } else {
  639. FPIX_residuals_dRz[make_tuple(layerOrDiskNr, hit_idx + 1, tm_type)]->fill(dRz, seed.Et());
  640. FPIX_residuals_dPhi[make_tuple(layerOrDiskNr, hit_idx + 1, tm_type)]->fill(dPhi, seed.Et());
  641. }
  642. }
  643. }
  644. }
  645. }
  646. tds.save_all();
  647. }
  648. int main(int argc, char * argv[]){
  649. using namespace fv_util;
  650. ArgParser args(argc, argv);
  651. if(args.cmd_option_exists("-c")) {
  652. init_config(args.get_cmd_option("-c"));
  653. args.update_config();
  654. if (args.cmd_option_exists("-s")) {
  655. the_config->update_key("silent", true);
  656. }
  657. if (args.cmd_option_exists("-b")) {
  658. the_config->update_key("batch", true);
  659. }
  660. init_log(LogPriority::kLogInfo);
  661. // gSystem->Load("libfilval.so");
  662. run();
  663. } else {
  664. cout << "Usage: ./" << argv[0] << " (-s) (-b) -c config_file.yaml" << endl;
  665. }
  666. return 0;
  667. }