tracking_eff.cpp 35 KB

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