tracking_validation.cpp 48 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871
  1. #include <iostream>
  2. #include <vector>
  3. #include <map>
  4. #include <utility>
  5. #include <numeric>
  6. #include <limits>
  7. #include <TSystem.h>
  8. #include "filval/filval.hpp"
  9. #include "filval/root/filval.hpp"
  10. #include <boost/range/combine.hpp>
  11. #include <boost/format.hpp>
  12. #include "analysis/TrackingNtuple.h"
  13. #include "analysis/TrackingNtupleObjs.hpp"
  14. #include "analysis/common.hpp"
  15. using namespace std;
  16. using namespace fv;
  17. using namespace fv::root;
  18. typedef std::pair<std::vector<float>,std::vector<float>> pair_vec;
  19. #define HIT_TYPE_PIXEL 0
  20. #define HIT_TYPE_GLUED 1
  21. #define HIT_TYPE_STRIP 2
  22. enum PixSubDet{
  23. Void = 0,
  24. Barrel = 1,
  25. EndCap = 2
  26. };
  27. struct PixDetElem{
  28. PixSubDet subdet;
  29. unsigned char layer;
  30. PixDetElem():subdet(PixSubDet::Void), layer(0) {}
  31. PixDetElem(const PixSubDet& subdet, const unsigned char& layer){
  32. this->subdet = subdet;
  33. this->layer = layer;
  34. }
  35. template<typename T>
  36. bool contains(const T& obj) const{
  37. return obj.subdet() == subdet && obj.layer() == layer;
  38. }
  39. bool contains(const PixSubDet& subdet, const unsigned char& layer) const{
  40. return subdet == this->subdet && layer == this->layer;
  41. }
  42. };
  43. std::map<PixSubDet,string> subdet_names = {{PixSubDet::Barrel, "BPIX"}, {PixSubDet::EndCap, "FPIX"}};
  44. typedef std::tuple<PixRecHit, SimHit> HitPair;
  45. vector<string> seedTypes =
  46. {"initialStepSeeds",
  47. "highPtTripletStepSeeds",
  48. "mixedTripletStepSeeds",
  49. "pixelLessStepSeeds",
  50. "tripletElectronSeeds",
  51. "pixelPairElectronSeeds",
  52. "stripPairElectronSeeds"};
  53. SuperClusterCollection super_clusters;
  54. SeedCollection seeds;
  55. PixRecHitCollection pixrec_hits;
  56. TrackCollection tracks;
  57. SimHitCollection sim_hits;
  58. SimTrackCollection sim_tracks;
  59. void setup_hit_pair_functions(){
  60. func<pair_vec(vector<HitPair>)>("calc_dphi_v_eta",
  61. FUNC(([](const vector<HitPair>& hit_pairs){
  62. vector<float> dphis;
  63. vector<float> etas;
  64. for(auto hit_pair : hit_pairs){
  65. auto& pixrec_hit = std::get<PixRecHit>(hit_pair);
  66. auto& sim_hit = std::get<SimHit>(hit_pair);
  67. dphis.push_back(atan2(sim_hit.x(), sim_hit.y()) - atan2(pixrec_hit.x(), pixrec_hit.y()));
  68. etas.push_back(pseudorapidity(pixrec_hit));
  69. }
  70. return std::make_pair(etas, dphis);
  71. })));
  72. func<pair_vec(vector<HitPair>)>("calc_dr_v_eta",
  73. FUNC(([](const vector<HitPair>& hit_pairs){
  74. vector<float> drs;
  75. vector<float> etas;
  76. for(auto hit_pair : hit_pairs){
  77. auto& pixrec_hit = std::get<PixRecHit>(hit_pair);
  78. auto& sim_hit = std::get<SimHit>(hit_pair);
  79. drs.push_back(displacement(pixrec_hit, sim_hit));
  80. etas.push_back(pseudorapidity(pixrec_hit));
  81. }
  82. return std::make_pair(etas, drs);
  83. })));
  84. func<pair_vec(vector<HitPair>)>("calc_dz_v_eta",
  85. FUNC(([](const vector<HitPair>& hit_pairs){
  86. vector<float> dzs;
  87. vector<float> etas;
  88. for(auto hit_pair : hit_pairs){
  89. auto& pixrec_hit = std::get<PixRecHit>(hit_pair);
  90. auto& sim_hit = std::get<SimHit>(hit_pair);
  91. dzs.push_back(sim_hit.z() - pixrec_hit.z());
  92. etas.push_back(pseudorapidity(pixrec_hit));
  93. }
  94. return std::make_pair(etas, dzs);
  95. })));
  96. func<pair_vec(vector<HitPair>)>("calc_drho_v_eta",
  97. FUNC(([](const vector<HitPair>& hit_pairs){
  98. vector<float> drhos;
  99. vector<float> etas;
  100. for(auto hit_pair : hit_pairs){
  101. auto& pixrec_hit = std::get<PixRecHit>(hit_pair);
  102. auto& sim_hit = std::get<SimHit>(hit_pair);
  103. drhos.push_back(rho(sim_hit) - rho(pixrec_hit));
  104. etas.push_back(pseudorapidity(pixrec_hit));
  105. }
  106. return std::make_pair(etas, drhos);
  107. })));
  108. }
  109. void register_objects(TrackingDataSet& tds){
  110. super_clusters.init(tds);
  111. seeds.init(tds);
  112. pixrec_hits.init(tds);
  113. tracks.init(tds);
  114. sim_hits.init(tds);
  115. sim_tracks.init(tds);
  116. }
  117. void setup_skipped_layer_hit_pairs(TrackingDataSet& tds){
  118. /* Finds SimHit/RecHit pairs
  119. */
  120. auto find_matched_nth_hit_in_layer_with_skip =
  121. func<vector<HitPair>(PixDetElem, int, bool)>("find_matched_nth_hit_in_layer_with_skip",
  122. FUNC(([](const PixDetElem&& pixdet,
  123. const unsigned int&& skip,
  124. const bool&& first){
  125. vector<HitPair> matched_hits;
  126. for(auto trk : tracks){ // loop over all tracks
  127. auto seed = seeds[trk.seedIdx()];
  128. // Require seed's source algorithm is one of those in seedTypes
  129. if(seed.algoOriginal() < 0 || seed.algoOriginal() >= seedTypes.size()) continue;
  130. // Require seed w/ at least 2 hits
  131. if(seed.hitIdx().size() < 2) continue;
  132. // take only pixel hits for now
  133. if(seed.hitType()[0] != HIT_TYPE_PIXEL || seed.hitType()[1] != HIT_TYPE_PIXEL) continue;
  134. auto rec_hit1 = pixrec_hits[seed.hitIdx()[0]];
  135. auto rec_hit2 = pixrec_hits[seed.hitIdx()[1]];
  136. if(first){ // Looking at first hit in the pair
  137. if(pixdet.contains(rec_hit1) && pixdet.contains((PixSubDet)rec_hit2.subdet(), rec_hit2.layer()-1-skip)){
  138. // We have the RecHit we want to work with, now find a properly* matched SimHit
  139. if(rec_hit1.simHitIdx().size() == 0) continue; // if rechit is matched to no simhits, give up.
  140. if(trk.simTrkIdx().size() == 0) continue; // if track is matched to no simtracks, give up.
  141. const int &simTrkIdx = trk.simTrkIdx()[0]; // get first matched track
  142. for(const int& simHitIdx : sim_tracks[simTrkIdx].simHitIdx()){ // loop over SimHits from SimTrack
  143. if(simHitIdx == rec_hit1.simHitIdx()[0]) // take first matched simhit (should be the closest one)
  144. matched_hits.push_back({rec_hit1, sim_hits[rec_hit1.simHitIdx()[0]]});
  145. }
  146. }
  147. }else{ // Looking at second hit in the pair
  148. if(pixdet.contains(rec_hit2) && pixdet.contains((PixSubDet)rec_hit1.subdet(), rec_hit2.layer()+1+skip)){
  149. // We have the RecHit we want to work with, now find a properly* matched SimHit
  150. if(rec_hit2.simHitIdx().size() == 0) continue; // if rechit is matched to no simhits, give up.
  151. if(trk.simTrkIdx().size() == 0) continue; // if track is matched to no simtracks, give up.
  152. const int &simTrkIdx = trk.simTrkIdx()[0]; // get first matched track
  153. for(const int& simHitIdx : sim_tracks[simTrkIdx].simHitIdx()){ // loop over SimHits from SimTrack
  154. if(simHitIdx == rec_hit1.simHitIdx()[0]) // take first matched simhit (should be the closest one)
  155. matched_hits.push_back({rec_hit2, sim_hits[rec_hit2.simHitIdx()[0]]});
  156. }
  157. }
  158. }
  159. }
  160. return matched_hits;
  161. })));
  162. auto BPIX_L1 = constant("BPIX_L1", PixDetElem(PixSubDet::Barrel, 1));
  163. auto BPIX_L2 = constant("BPIX_L2", PixDetElem(PixSubDet::Barrel, 2));
  164. auto FPIX_L1 = constant("FPIX_L1", PixDetElem(PixSubDet::EndCap, 1));
  165. auto FPIX_L2 = constant("FPIX_L2", PixDetElem(PixSubDet::EndCap, 2));
  166. auto skip_zero = constant("skip_zero", 0);
  167. auto skip_one = constant("skip_one", 1);
  168. auto first = constant("first", true);
  169. auto second = constant("second", false);
  170. // First hit on 1st/2nd bpix layers and second hit in 2nd/3rd
  171. auto first_hits_in_B1_skip_zero = fv::tup_apply(find_matched_nth_hit_in_layer_with_skip,
  172. fv::tuple(BPIX_L1, skip_zero, first), "first_hits_in_B1_skip_zero");
  173. auto first_hits_in_B2_skip_zero = fv::tup_apply(find_matched_nth_hit_in_layer_with_skip,
  174. fv::tuple(BPIX_L2, skip_zero, first), "first_hits_in_B2_skip_zero");
  175. // First hit on 1st/2nd fpix layers and second hit in 2nd/3rd
  176. auto first_hits_in_F1_skip_zero = fv::tup_apply(find_matched_nth_hit_in_layer_with_skip,
  177. fv::tuple(FPIX_L1, skip_zero, first), "first_hits_in_F1_skip_zero");
  178. auto first_hits_in_F2_skip_zero = fv::tup_apply(find_matched_nth_hit_in_layer_with_skip,
  179. fv::tuple(FPIX_L2, skip_zero, first), "first_hits_in_F2_skip_zero");
  180. // First hit on 1st/2nd fpix layers and second hit in 3rd/4th
  181. auto first_hits_in_B1_skip_one = fv::tup_apply(find_matched_nth_hit_in_layer_with_skip,
  182. fv::tuple(BPIX_L1, skip_one, first), "first_hits_in_B1_skip_one");
  183. auto first_hits_in_B2_skip_one = fv::tup_apply(find_matched_nth_hit_in_layer_with_skip,
  184. fv::tuple(BPIX_L2, skip_one, first), "first_hits_in_B2_skip_one");
  185. // First hit on 1st fpix layer and second hit in 3rd layer
  186. auto first_hits_in_F1_skip_one = fv::tup_apply(find_matched_nth_hit_in_layer_with_skip,
  187. fv::tuple(FPIX_L1, skip_one, first), "first_hits_in_F1_skip_one");
  188. TH2Params params_dz = {"$\\eta$", 100, -4, 4,
  189. "$\\Delta z$(rad)", 50, -.01, .01};
  190. auto calc_dz_v_eta = lookup_function<pair_vec(vector<HitPair>)>("calc_dz_v_eta");
  191. tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B1_skip_zero",
  192. fv::apply(calc_dz_v_eta, first_hits_in_B1_skip_zero),
  193. "First Hit in BPIX-L1 - Skip 0", params_dz);
  194. tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B2_skip_zero",
  195. fv::apply(calc_dz_v_eta, first_hits_in_B2_skip_zero),
  196. "First Hit in BPIX-L2 - Skip 0", params_dz);
  197. tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_F1_skip_zero",
  198. fv::apply(calc_dz_v_eta, first_hits_in_F1_skip_zero),
  199. "First Hit in FPIX-L1 - Skip 0", params_dz);
  200. tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_F2_skip_zero",
  201. fv::apply(calc_dz_v_eta, first_hits_in_F2_skip_zero),
  202. "First Hit in FPIX-L2 - Skip 0", params_dz);
  203. tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B1_skip_one",
  204. fv::apply(calc_dz_v_eta, first_hits_in_B1_skip_one),
  205. "First Hit in BPIX-L1 - Skip 1", params_dz);
  206. tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B2_skip_one",
  207. fv::apply(calc_dz_v_eta, first_hits_in_B2_skip_one),
  208. "First Hit in BPIX-L2 - Skip 1", params_dz);
  209. tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_F1_skip_one",
  210. fv::apply(calc_dz_v_eta, first_hits_in_F1_skip_one),
  211. "First Hit in FPIX-L1 - Skip 1", params_dz);
  212. }
  213. void setup_first_hit_pairs(TrackingDataSet& tds){
  214. auto find_matched_nth_hit_in_layer =
  215. func<vector<HitPair>(PixDetElem, int)>("find_matched_nth_hit_in_layer",
  216. FUNC(([](const PixDetElem&& pixdet,
  217. const unsigned int&& hit_number){
  218. vector<HitPair> matched_hits;
  219. for(auto track : tracks){ // loop over all tracks
  220. auto seed = seeds[track.seedIdx()];
  221. if(seed.hitIdx().size() <= hit_number) continue; // looking for hit_number'th hit
  222. if(seed.algoOriginal() < 0 || seed.algoOriginal() >= seedTypes.size()) continue;
  223. if(seed.hitType()[hit_number] != HIT_TYPE_PIXEL) continue; // take only pixel hits for now
  224. auto rec_hit = pixrec_hits[seed.hitIdx()[hit_number]];
  225. if(pixdet.contains(rec_hit)){
  226. // We have the RecHit we want to work with, now find a properly* matched SimHit
  227. if(rec_hit.simHitIdx().size() == 0) continue; // if rechit is matched to no simhits, give up.
  228. if(track.simTrkIdx().size() == 0) continue; // if rechit is matched to no simhits, give up.
  229. const int &simTrkIdx = track.simTrkIdx()[0]; // get first matched track
  230. for(const int& simHitIdx : sim_tracks[simTrkIdx].simHitIdx()){ // loop over SimHits from SimTrack
  231. if(simHitIdx == rec_hit.simHitIdx()[0]) // take first matched simhit (should be the closest one)
  232. matched_hits.push_back({rec_hit, sim_hits[rec_hit.simHitIdx()[0]]});
  233. }
  234. }
  235. }
  236. return matched_hits;
  237. })));
  238. //**//**
  239. auto BPIX_L1 = constant("BPIX_L1", PixDetElem(PixSubDet::Barrel, 1));
  240. auto BPIX_L2 = constant("BPIX_L2", PixDetElem(PixSubDet::Barrel, 2));
  241. auto BPIX_L3 = constant("BPIX_L3", PixDetElem(PixSubDet::Barrel, 3));
  242. auto BPIX_L4 = constant("BPIX_L4", PixDetElem(PixSubDet::Barrel, 4));
  243. auto FPIX_L1 = constant("FPIX_L1", PixDetElem(PixSubDet::EndCap, 1));
  244. auto FPIX_L2 = constant("FPIX_L2", PixDetElem(PixSubDet::EndCap, 2));
  245. auto FPIX_L3 = constant("FPIX_L3", PixDetElem(PixSubDet::EndCap, 3));
  246. auto first_hit = constant("1st", 0);
  247. auto second_hit = constant("2nd", 1);
  248. // First hits on inner two bpix layers
  249. auto first_hits_in_B1 = fv::tup_apply(find_matched_nth_hit_in_layer,
  250. fv::tuple(BPIX_L1, first_hit), "first_hits_in_B1");
  251. auto first_hits_in_B2 = fv::tup_apply(find_matched_nth_hit_in_layer,
  252. fv::tuple(BPIX_L2, first_hit), "first_hits_in_B2");
  253. // Second hits on outer three bpix layers
  254. auto second_hits_in_B2 = fv::tup_apply(find_matched_nth_hit_in_layer,
  255. fv::tuple(BPIX_L2, second_hit), "second_hits_in_B2");
  256. auto second_hits_in_B3 = fv::tup_apply(find_matched_nth_hit_in_layer,
  257. fv::tuple(BPIX_L3, second_hit), "second_hits_in_B3");
  258. auto second_hits_in_B4 = fv::tup_apply(find_matched_nth_hit_in_layer,
  259. fv::tuple(BPIX_L4, second_hit), "second_hits_in_B4");
  260. auto first_hits_in_F1 = fv::tup_apply(find_matched_nth_hit_in_layer,
  261. fv::tuple(FPIX_L1, first_hit), "first_hits_in_F1");
  262. auto first_hits_in_F2 = fv::tup_apply(find_matched_nth_hit_in_layer,
  263. fv::tuple(FPIX_L2, first_hit), "first_hits_in_F2");
  264. auto second_hits_in_F2 = fv::tup_apply(find_matched_nth_hit_in_layer,
  265. fv::tuple(FPIX_L2, second_hit), "second_hits_in_F2");
  266. auto second_hits_in_F3 = fv::tup_apply(find_matched_nth_hit_in_layer,
  267. fv::tuple(FPIX_L3, second_hit), "second_hits_in_F3");
  268. // Even vs Odd Ladders
  269. auto even = constant("even", false);
  270. auto odd = constant("odd", true);
  271. auto select_even_odd_ladder_hit_pairs =
  272. func<vector<HitPair>(vector<HitPair>, bool)>("select_even_odd_ladder_hit_pairs",
  273. FUNC(([](const vector<HitPair>& hit_pairs,
  274. const bool &&odd){
  275. vector<HitPair> even_pairs;
  276. for(const HitPair &hit_pair : hit_pairs){ // loop over all seeds
  277. if(std::get<PixRecHit>(hit_pair).ladder() % 2 == odd){
  278. even_pairs.push_back(hit_pair);
  279. }
  280. }
  281. return even_pairs;
  282. })));
  283. auto first_hits_in_B1_even_ladder = fv::tup_apply(select_even_odd_ladder_hit_pairs,
  284. fv::tuple(first_hits_in_B1, even), "first_hits_in_B1_even_ladder");
  285. auto first_hits_in_B1_odd_ladder = fv::tup_apply(select_even_odd_ladder_hit_pairs,
  286. fv::tuple(first_hits_in_B1, odd), "first_hits_in_B1_odd_ladder");
  287. auto first_hits_in_B2_even_ladder = fv::tup_apply(select_even_odd_ladder_hit_pairs,
  288. fv::tuple(first_hits_in_B2, even), "first_hits_in_B2_even_ladder");
  289. auto first_hits_in_B2_odd_ladder = fv::tup_apply(select_even_odd_ladder_hit_pairs,
  290. fv::tuple(first_hits_in_B2, odd), "first_hits_in_B2_odd_ladder");
  291. //Plots for dPhi of collections defined above
  292. auto calc_dphi_v_eta = lookup_function<pair_vec(vector<HitPair>)>("calc_dphi_v_eta");
  293. TH2Params params_dphi = {"$\\eta$", 100, -4, 4,
  294. "$\\Delta \\phi$(rad)", 150, -.0015, .0015};
  295. tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_first_hits_in_B1",
  296. fv::apply(calc_dphi_v_eta, first_hits_in_B1),
  297. "First Hit in BPIX-L1", params_dphi);
  298. tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_first_hits_in_B2",
  299. fv::apply(calc_dphi_v_eta, first_hits_in_B2),
  300. "First Hit in BPIX-L2", params_dphi);
  301. tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_first_hits_in_B1_even_ladder",
  302. fv::apply(calc_dphi_v_eta, first_hits_in_B1_even_ladder),
  303. "First Hit in BPIX-L1 - Even Ladders", params_dphi);
  304. tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_first_hits_in_B1_odd_ladder",
  305. fv::apply(calc_dphi_v_eta, first_hits_in_B1_odd_ladder),
  306. "First Hit in BPIX-L1 - Odd Ladders", params_dphi);
  307. tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_first_hits_in_B2_even_ladder",
  308. fv::apply(calc_dphi_v_eta, first_hits_in_B2_even_ladder),
  309. "First Hit in BPIX-L2 - Even Ladders", params_dphi);
  310. tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_first_hits_in_B2_odd_ladder",
  311. fv::apply(calc_dphi_v_eta, first_hits_in_B2_odd_ladder),
  312. "First Hit in BPIX-L2 - Odd Ladders", params_dphi);
  313. tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_second_hits_in_B2",
  314. fv::apply(calc_dphi_v_eta, second_hits_in_B2),
  315. "Second Hit in BPIX-L2", params_dphi);
  316. tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_second_hits_in_B3",
  317. fv::apply(calc_dphi_v_eta, second_hits_in_B3),
  318. "Second Hit in BPIX-L3", params_dphi);
  319. tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_second_hits_in_B4",
  320. fv::apply(calc_dphi_v_eta, second_hits_in_B4),
  321. "Second Hit in BPIX-L4", params_dphi);
  322. tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_first_hits_in_F1",
  323. fv::apply(calc_dphi_v_eta, first_hits_in_F1),
  324. "First Hit in FPIX-L1", params_dphi);
  325. tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_first_hits_in_F2",
  326. fv::apply(calc_dphi_v_eta, first_hits_in_F2),
  327. "First Hit in FPIX-L2", params_dphi);
  328. tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_second_hits_in_F2",
  329. fv::apply(calc_dphi_v_eta, second_hits_in_F2),
  330. "Second Hit in FPIX-L2", params_dphi);
  331. tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_second_hits_in_F3",
  332. fv::apply(calc_dphi_v_eta, second_hits_in_F3),
  333. "Second Hit in FPIX-L3", params_dphi);
  334. //Plots for dr of collections defined above
  335. auto calc_dr_v_eta = lookup_function<pair_vec(vector<HitPair>)>("calc_dr_v_eta");
  336. TH2Params params_dr = {"$\\eta$", 100, -4, 4,
  337. "$\\Delta \\r$(rad)", 50, 0, .006};
  338. tds.register_container<ContainerTH2Many<float>>("dr_v_eta_first_hits_in_B1",
  339. fv::apply(calc_dr_v_eta, first_hits_in_B1),
  340. "First Hit in BPIX-L1", params_dr);
  341. tds.register_container<ContainerTH2Many<float>>("dr_v_eta_first_hits_in_B2",
  342. fv::apply(calc_dr_v_eta, first_hits_in_B2),
  343. "First Hit in BPIX-L2", params_dr);
  344. tds.register_container<ContainerTH2Many<float>>("dr_v_eta_first_hits_in_B1_even_ladder",
  345. fv::apply(calc_dr_v_eta, first_hits_in_B1_even_ladder),
  346. "First Hit in BPIX-L1 - Even Ladders", params_dr);
  347. tds.register_container<ContainerTH2Many<float>>("dr_v_eta_first_hits_in_B1_odd_ladder",
  348. fv::apply(calc_dr_v_eta, first_hits_in_B1_odd_ladder),
  349. "First Hit in BPIX-L1 - Odd Ladders", params_dr);
  350. tds.register_container<ContainerTH2Many<float>>("dr_v_eta_first_hits_in_B2_even_ladder",
  351. fv::apply(calc_dr_v_eta, first_hits_in_B2_even_ladder),
  352. "First Hit in BPIX-L2 - Even Ladders", params_dr);
  353. tds.register_container<ContainerTH2Many<float>>("dr_v_eta_first_hits_in_B2_odd_ladder",
  354. fv::apply(calc_dr_v_eta, first_hits_in_B2_odd_ladder),
  355. "First Hit in BPIX-L2 - Odd Ladders", params_dr);
  356. tds.register_container<ContainerTH2Many<float>>("dr_v_eta_second_hits_in_B2",
  357. fv::apply(calc_dr_v_eta, second_hits_in_B2),
  358. "Second Hit in BPIX-L2", params_dr);
  359. tds.register_container<ContainerTH2Many<float>>("dr_v_eta_second_hits_in_B3",
  360. fv::apply(calc_dr_v_eta, second_hits_in_B3),
  361. "Second Hit in BPIX-L3", params_dr);
  362. tds.register_container<ContainerTH2Many<float>>("dr_v_eta_second_hits_in_B4",
  363. fv::apply(calc_dr_v_eta, second_hits_in_B4),
  364. "Second Hit in BPIX-L4", params_dr);
  365. tds.register_container<ContainerTH2Many<float>>("dr_v_eta_first_hits_in_F1",
  366. fv::apply(calc_dr_v_eta, first_hits_in_F1),
  367. "First Hit in FPIX-L1", params_dr);
  368. tds.register_container<ContainerTH2Many<float>>("dr_v_eta_first_hits_in_F2",
  369. fv::apply(calc_dr_v_eta, first_hits_in_F2),
  370. "First Hit in FPIX-L2", params_dr);
  371. tds.register_container<ContainerTH2Many<float>>("dr_v_eta_second_hits_in_F2",
  372. fv::apply(calc_dr_v_eta, second_hits_in_F2),
  373. "Second Hit in FPIX-L2", params_dr);
  374. tds.register_container<ContainerTH2Many<float>>("dr_v_eta_second_hits_in_F3",
  375. fv::apply(calc_dr_v_eta, second_hits_in_F3),
  376. "Second Hit in FPIX-L3", params_dr);
  377. //Plots for dz of collections defined above
  378. auto calc_dz_v_eta = lookup_function<pair_vec(vector<HitPair>)>("calc_dz_v_eta");
  379. TH2Params params_dz = {"$\\eta$", 100, -4, 4,
  380. "$\\Delta z$(rad)", 50, -.01, .01};
  381. tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B1",
  382. fv::apply(calc_dz_v_eta, first_hits_in_B1),
  383. "First Hit in BPIX-L1", params_dz);
  384. tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B2",
  385. fv::apply(calc_dz_v_eta, first_hits_in_B2),
  386. "First Hit in BPIX-L2", params_dz);
  387. tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B1_even_ladder",
  388. fv::apply(calc_dz_v_eta, first_hits_in_B1_even_ladder),
  389. "First Hit in BPIX-L1 - Even Ladders", params_dz);
  390. tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B1_odd_ladder",
  391. fv::apply(calc_dz_v_eta, first_hits_in_B1_odd_ladder),
  392. "First Hit in BPIX-L1 - Odd Ladders", params_dz);
  393. tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B2_even_ladder",
  394. fv::apply(calc_dz_v_eta, first_hits_in_B2_even_ladder),
  395. "First Hit in BPIX-L2 - Even Ladders", params_dz);
  396. tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B2_odd_ladder",
  397. fv::apply(calc_dz_v_eta, first_hits_in_B2_odd_ladder),
  398. "First Hit in BPIX-L2 - Odd Ladders", params_dz);
  399. tds.register_container<ContainerTH2Many<float>>("dz_v_eta_second_hits_in_B2",
  400. fv::apply(calc_dz_v_eta, second_hits_in_B2),
  401. "Second Hit in BPIX-L2", params_dz);
  402. tds.register_container<ContainerTH2Many<float>>("dz_v_eta_second_hits_in_B3",
  403. fv::apply(calc_dz_v_eta, second_hits_in_B3),
  404. "Second Hit in BPIX-L3", params_dz);
  405. tds.register_container<ContainerTH2Many<float>>("dz_v_eta_second_hits_in_B4",
  406. fv::apply(calc_dz_v_eta, second_hits_in_B4),
  407. "Second Hit in BPIX-L4", params_dz);
  408. //Plots for drho of collections defined above
  409. auto calc_drho_v_eta = lookup_function<pair_vec(vector<HitPair>)>("calc_drho_v_eta");
  410. TH2Params params_drho = {"$\\eta$", 100, -4, 4,
  411. "$\\Delta \\rho$(cm)", 50, -.01, .01};
  412. tds.register_container<ContainerTH2Many<float>>("drho_v_eta_first_hits_in_F1",
  413. fv::apply(calc_drho_v_eta, first_hits_in_F1),
  414. "First Hit in FPIX-L1", params_drho);
  415. tds.register_container<ContainerTH2Many<float>>("drho_v_eta_first_hits_in_F2",
  416. fv::apply(calc_drho_v_eta, first_hits_in_F2),
  417. "First Hit in FPIX-L2", params_drho);
  418. tds.register_container<ContainerTH2Many<float>>("drho_v_eta_second_hits_in_F2",
  419. fv::apply(calc_drho_v_eta, second_hits_in_F2),
  420. "Second Hit in FPIX-L2", params_drho);
  421. tds.register_container<ContainerTH2Many<float>>("drho_v_eta_second_hits_in_F3",
  422. fv::apply(calc_drho_v_eta, second_hits_in_F3),
  423. "Second Hit in FPIX-L3", params_drho);
  424. }
  425. void sc_value_hists_for_lb(TrackingDataSet& tds, std::function<bool(int)>& parity_check,
  426. const std::string& label){
  427. auto relabel = [label](const std::string&& base){
  428. return base + "_" + label;
  429. };
  430. auto BPIX_L1 = constant("BPIX_L1", PixDetElem(PixSubDet::Barrel, 1));
  431. auto FPIX_L1 = constant("FPIX_L1", PixDetElem(PixSubDet::EndCap, 1));
  432. auto BPIX_L2 = constant("BPIX_L2", PixDetElem(PixSubDet::Barrel, 2));
  433. auto FPIX_L2 = constant("FPIX_L2", PixDetElem(PixSubDet::EndCap, 2));
  434. auto BPIX_L3 = constant("BPIX_L3", PixDetElem(PixSubDet::Barrel, 3));
  435. auto FPIX_L3 = constant("FPIX_L2", PixDetElem(PixSubDet::EndCap, 3));
  436. auto BPIX_L4 = constant("BPIX_L4", PixDetElem(PixSubDet::Barrel, 4));
  437. auto first_hit = constant("1st", 0);
  438. auto second_hit = constant("2nd", 1);
  439. auto dRz = constant("dRz", true);
  440. auto dPhi = constant("dPhi", false);
  441. auto pcheck = constant(label, parity_check);
  442. auto sc_hit_residuals = lookup_function<
  443. pair_vec(PixDetElem, int, bool, std::function<bool(int)>)>("sc_hit_residuals");
  444. // First hits on inner three bpix layers
  445. auto sc_first_hits_in_B1_dz = fv::tup_apply(sc_hit_residuals,
  446. fv::tuple(BPIX_L1, first_hit, dRz, pcheck), relabel("sc_first_hits_in_B1_dz"));
  447. auto sc_first_hits_in_B2_dz = fv::tup_apply(sc_hit_residuals,
  448. fv::tuple(BPIX_L2, first_hit, dRz, pcheck), relabel("sc_first_hits_in_B2_dz"));
  449. auto sc_first_hits_in_B3_dz = fv::tup_apply(sc_hit_residuals,
  450. fv::tuple(BPIX_L3, first_hit, dRz, pcheck), relabel("sc_first_hits_in_B3_dz"));
  451. auto sc_first_hits_in_B1_dphi = fv::tup_apply(sc_hit_residuals,
  452. fv::tuple(BPIX_L1, first_hit, dPhi, pcheck), relabel("sc_first_hits_in_B1_dphi"));
  453. auto sc_first_hits_in_B2_dphi = fv::tup_apply(sc_hit_residuals,
  454. fv::tuple(BPIX_L2, first_hit, dPhi, pcheck), relabel("sc_first_hits_in_B2_dphi"));
  455. auto sc_first_hits_in_B3_dphi = fv::tup_apply(sc_hit_residuals,
  456. fv::tuple(BPIX_L3, first_hit, dPhi, pcheck), relabel("sc_first_hits_in_B3_dphi"));
  457. // Second hits on outer three bpix layers
  458. auto sc_second_hits_in_B2_dz = fv::tup_apply(sc_hit_residuals,
  459. fv::tuple(BPIX_L2, second_hit, dRz, pcheck), relabel("sc_second_hits_in_B2_dz"));
  460. auto sc_second_hits_in_B3_dz = fv::tup_apply(sc_hit_residuals,
  461. fv::tuple(BPIX_L3, second_hit, dRz, pcheck), relabel("sc_second_hits_in_B3_dz"));
  462. auto sc_second_hits_in_B4_dz = fv::tup_apply(sc_hit_residuals,
  463. fv::tuple(BPIX_L4, second_hit, dRz, pcheck), relabel("sc_second_hits_in_B4_dz"));
  464. auto sc_second_hits_in_B2_dphi = fv::tup_apply(sc_hit_residuals,
  465. fv::tuple(BPIX_L2, second_hit, dPhi, pcheck), relabel("sc_second_hits_in_B2_dphi"));
  466. auto sc_second_hits_in_B3_dphi = fv::tup_apply(sc_hit_residuals,
  467. fv::tuple(BPIX_L3, second_hit, dPhi, pcheck), relabel("sc_second_hits_in_B3_dphi"));
  468. auto sc_second_hits_in_B4_dphi = fv::tup_apply(sc_hit_residuals,
  469. fv::tuple(BPIX_L4, second_hit, dPhi, pcheck), relabel("sc_second_hits_in_B4_dphi"));
  470. TH2Params params_dz = {"$\\eta$", 100, -4, 4,
  471. "$\\delta z$(cm)", 100, -0.7, 0.7};
  472. tds.register_container<ContainerTH2Many<float>>(relabel("sc_first_hits_in_B1_dz"),
  473. sc_first_hits_in_B1_dz,
  474. "First Hit in BPIX-L1", params_dz);
  475. tds.register_container<ContainerTH2Many<float>>(relabel("sc_first_hits_in_B2_dz"),
  476. sc_first_hits_in_B2_dz,
  477. "First Hit in BPIX-L2", params_dz);
  478. tds.register_container<ContainerTH2Many<float>>(relabel("sc_first_hits_in_B3_dz"),
  479. sc_first_hits_in_B3_dz,
  480. "First Hit in BPIX-L3", params_dz);
  481. tds.register_container<ContainerTH2Many<float>>(relabel("sc_second_hits_in_B2_dz"),
  482. sc_second_hits_in_B2_dz,
  483. "Second Hit in BPIX-L2", params_dz);
  484. tds.register_container<ContainerTH2Many<float>>(relabel("sc_second_hits_in_B3_dz"),
  485. sc_second_hits_in_B3_dz,
  486. "Second Hit in BPIX-L3", params_dz);
  487. tds.register_container<ContainerTH2Many<float>>(relabel("sc_second_hits_in_B4_dz"),
  488. sc_second_hits_in_B4_dz,
  489. "Second Hit in BPIX-L4", params_dz);
  490. TH2Params params_dphi = {"$\\eta$", 100, -4, 4,
  491. "$\\delta \\phi$(rad)", 1500, -0.7, 0.7};
  492. tds.register_container<ContainerTH2Many<float>>(relabel("sc_first_hits_in_B1_dphi"),
  493. sc_first_hits_in_B1_dphi,
  494. "First Hit in BPIX-L1", params_dphi);
  495. tds.register_container<ContainerTH2Many<float>>(relabel("sc_first_hits_in_B2_dphi"),
  496. sc_first_hits_in_B2_dphi,
  497. "First Hit in BPIX-L2", params_dphi);
  498. tds.register_container<ContainerTH2Many<float>>(relabel("sc_first_hits_in_B3_dphi"),
  499. sc_first_hits_in_B3_dphi,
  500. "First Hit in BPIX-L3", params_dphi);
  501. tds.register_container<ContainerTH2Many<float>>(relabel("sc_second_hits_in_B2_dphi"),
  502. sc_second_hits_in_B2_dphi,
  503. "Second Hit in BPIX-L2", params_dphi);
  504. tds.register_container<ContainerTH2Many<float>>(relabel("sc_second_hits_in_B3_dphi"),
  505. sc_second_hits_in_B3_dphi,
  506. "Second Hit in BPIX-L3", params_dphi);
  507. tds.register_container<ContainerTH2Many<float>>(relabel("sc_second_hits_in_B4_dphi"),
  508. sc_second_hits_in_B4_dphi,
  509. "Second Hit in BPIX-L4", params_dphi);
  510. }
  511. void setup_sc_residuals(TrackingDataSet& tds){
  512. /* sc_hit_residuals looks through all the super clusters in the event and
  513. * examins the matched pixel hits. It calculates a residual based on the
  514. * extrapolated position of the SC track and the position of the matched
  515. * hit.
  516. * subdet: 1 for barrel, 2 for endcap
  517. * pix_layer: 1-4 for barrel, 1-3 for endcap
  518. * hit_number: examin either the first or second matched hit
  519. * sel_dRz: True: calculate dR, False: calculate dz
  520. * pcheck: a boolean function that filters residuals based on
  521. * ladder/blade. Normally used to select either even or odd ladders.
  522. */
  523. func<pair_vec(PixDetElem, int, bool, std::function<bool(int)>)>("sc_hit_residuals",
  524. FUNC(([](const PixDetElem&& pixdet,
  525. const int&& hit_number,
  526. const bool&& sel_dRz,
  527. const std::function<bool(int)>& pcheck){
  528. vector<float> residuals;
  529. vector<float> etas;
  530. for(auto super_cluster : super_clusters){ // loop over all super-clusters
  531. float dRz, dPhi;
  532. int ladder;
  533. unsigned int nSeeds = super_cluster.charge().size();
  534. for(unsigned int seedIdx=0; seedIdx<nSeeds; seedIdx++){ // loop over all seeds associated w/ sc
  535. if(hit_number == 0){
  536. if(!pixdet.contains((PixSubDet)super_cluster.subDet1()[seedIdx], super_cluster.lay1()[seedIdx])) continue;
  537. dRz = super_cluster.dRz1()[seedIdx];
  538. dPhi = super_cluster.dPhi1()[seedIdx];
  539. /* ladder = super_cluster.ladder1()[seedIdx]; */
  540. }else{
  541. if(!pixdet.contains((PixSubDet)super_cluster.subDet2()[seedIdx], super_cluster.lay2()[seedIdx])) continue;
  542. dRz = super_cluster.dRz2()[seedIdx];
  543. dPhi = super_cluster.dPhi2()[seedIdx];
  544. /* ladder = super_cluster.ladder2()[seedIdx]; */
  545. }
  546. ladder = 0; // TODO: Need to re-add ladder to tree for super-clusters
  547. DEBUG("New Seed, Idx: " << seedIdx << endl
  548. << "\tdRz: " << dRz << endl
  549. << "\tdPhi: " << dPhi << endl);
  550. if(!pcheck(ladder)) continue;
  551. if(sel_dRz)
  552. residuals.push_back(dRz);
  553. else
  554. residuals.push_back(dPhi);
  555. etas.push_back(pseudorapidityP(super_cluster));
  556. }
  557. }
  558. return std::make_pair(etas,residuals);
  559. })));
  560. std::function<bool(int)> even = [](const int& t){ return !(t%2); };
  561. std::function<bool(int)> odd = [](const int& t){ return (t%2); };
  562. sc_value_hists_for_lb(tds, even, "even");
  563. sc_value_hists_for_lb(tds, odd, "odd");
  564. }
  565. void setup_residual_cross_corrolations(TrackingDataSet& tds){
  566. /*
  567. * Interesting cross-corrolations
  568. * for hit1 in L1 or L2
  569. * dphi2 v dz1
  570. * dz2 v dz1
  571. * dphi2 v dphi1
  572. * dz2 v dphi1
  573. */
  574. enum Res{
  575. dz1,
  576. dz2,
  577. dPhi1,
  578. dPhi2
  579. };
  580. auto sc_cross_correlations = func<pair_vec(PixDetElem, PixDetElem, Res, Res)>("sc_cross_correlations",
  581. FUNC(([](const PixDetElem&& subdet1,
  582. const PixDetElem&& subdet2,
  583. const Res&& res1,
  584. const Res&& res2){
  585. vector<float> residuals1;
  586. vector<float> residuals2;
  587. for(auto super_cluster : super_clusters){ // loop over all supser-clusters
  588. unsigned int nSeeds = super_cluster.charge().size();
  589. for(unsigned int seedIdx=0; seedIdx<nSeeds; seedIdx++){ // loop over all seeds associated w/ sc
  590. if (!subdet1.contains((PixSubDet)super_cluster.subDet1()[seedIdx],super_cluster.lay1()[seedIdx]) ||
  591. !subdet2.contains((PixSubDet)super_cluster.subDet2()[seedIdx],super_cluster.lay2()[seedIdx]))
  592. continue;
  593. switch(res1){
  594. case Res::dz1:
  595. residuals1.push_back(super_cluster.dRz1()[seedIdx]);
  596. break;
  597. case Res::dz2:
  598. residuals1.push_back(super_cluster.dRz2()[seedIdx]);
  599. break;
  600. case Res::dPhi1:
  601. residuals1.push_back(super_cluster.dPhi1()[seedIdx]);
  602. break;
  603. case Res::dPhi2:
  604. residuals1.push_back(super_cluster.dPhi2()[seedIdx]);
  605. break;
  606. }
  607. switch(res2){
  608. case Res::dz1:
  609. residuals2.push_back(super_cluster.dRz1()[seedIdx]);
  610. break;
  611. case Res::dz2:
  612. residuals2.push_back(super_cluster.dRz2()[seedIdx]);
  613. break;
  614. case Res::dPhi1:
  615. residuals2.push_back(super_cluster.dPhi1()[seedIdx]);
  616. break;
  617. case Res::dPhi2:
  618. residuals2.push_back(super_cluster.dPhi2()[seedIdx]);
  619. break;
  620. }
  621. }
  622. }
  623. return std::make_pair(residuals1,residuals2);
  624. })));
  625. auto BPIX_L1 = constant("BPIX_L1", PixDetElem(PixSubDet::Barrel, 1));
  626. auto BPIX_L2 = constant("BPIX_L2", PixDetElem(PixSubDet::Barrel, 2));
  627. auto BPIX_L3 = constant("BPIX_L3", PixDetElem(PixSubDet::Barrel, 3));
  628. auto BPIX_L4 = constant("BPIX_L4", PixDetElem(PixSubDet::Barrel, 4));
  629. auto v_dz1 = constant("dz1", Res::dz1);
  630. auto v_dz2 = constant("dz2", Res::dz2);
  631. auto v_dPhi1 = constant("dPhi1", Res::dPhi1);
  632. auto v_dPhi2 = constant("dPhi2", Res::dPhi2);
  633. auto sc_dphi2_v_dz1_L1_L2 = fv::tup_apply(sc_cross_correlations,
  634. fv::tuple(BPIX_L1, BPIX_L2, v_dPhi2, v_dz1), "sc_dphi1_v_dz1_L1_L2");
  635. auto sc_dz2_v_dz1_L1_L2 = fv::tup_apply(sc_cross_correlations,
  636. fv::tuple(BPIX_L1, BPIX_L2, v_dz2, v_dz1), "sc_dz2_v_dz1_L1_L2");
  637. auto sc_dphi2_v_dphi1_L1_L2 = fv::tup_apply(sc_cross_correlations,
  638. fv::tuple(BPIX_L1, BPIX_L2, v_dPhi2, v_dPhi1), "sc_dphi1_v_dphi1_L1_L2");
  639. auto sc_dz2_v_dphi1_L1_L2 = fv::tup_apply(sc_cross_correlations,
  640. fv::tuple(BPIX_L1, BPIX_L2, v_dz2, v_dPhi1), "sc_dz2_v_dphi1_L1_L2");
  641. TH2Params hist_params = TH2Params::lookup("sc_dphi2_v_dz1_L1_L2");
  642. tds.register_container<ContainerTH2Many<float>>("sc_dphi2_v_dz1_L1_L2",
  643. sc_dphi2_v_dz1_L1_L2,
  644. "sc_dphi2_v_dz1_L1_L2", hist_params);
  645. hist_params = TH2Params::lookup("sc_dz2_v_dz1_L1_L2");
  646. tds.register_container<ContainerTH2Many<float>>("sc_dz2_v_dz1_L1_L2",
  647. sc_dz2_v_dz1_L1_L2,
  648. "sc_dz2_v_dz1_L1_L2", hist_params);
  649. hist_params = TH2Params::lookup("sc_dphi2_v_dphi1_L1_L2");
  650. tds.register_container<ContainerTH2Many<float>>("sc_dphi2_v_dphi1_L1_L2",
  651. sc_dphi2_v_dphi1_L1_L2,
  652. "sc_dphi2_v_dphi1_L1_L2", hist_params);
  653. hist_params = TH2Params::lookup("sc_dz2_v_dphi1_L1_L2");
  654. tds.register_container<ContainerTH2Many<float>>("sc_dz2_v_dphi1_L1_L2",
  655. sc_dz2_v_dphi1_L1_L2,
  656. "sc_dz2_v_dphi1_L1_L2", hist_params);
  657. hist_params = TH2Params::lookup("sc_dz2_v_dz1_L1_L2_zoom");
  658. tds.register_container<ContainerTH2Many<float>>("sc_dz2_v_dz1_L1_L2_zoom",
  659. sc_dz2_v_dz1_L1_L2,
  660. "sc_dz2_v_dz1_L1_L2_zoom", hist_params);
  661. hist_params = TH2Params::lookup("sc_dz2_v_dphi1_L1_L2_zoom");
  662. tds.register_container<ContainerTH2Many<float>>("sc_dz2_v_dphi1_L1_L2_zoom",
  663. sc_dz2_v_dphi1_L1_L2,
  664. "sc_dz2_v_dphi1_L1_L2_zoom", hist_params);
  665. hist_params = TH2Params::lookup("sc_dphi2_v_dphi1_L1_L2_zoom");
  666. tds.register_container<ContainerTH2Many<float>>("sc_dphi2_v_dphi1_L1_L2_zoom",
  667. sc_dphi2_v_dphi1_L1_L2,
  668. "sc_dphi2_v_dphi1_L1_L2_zoom", hist_params);
  669. auto sc_dphi2_v_dz1_L1_L3 = fv::tup_apply(sc_cross_correlations,
  670. fv::tuple(BPIX_L1, BPIX_L3, v_dPhi2, v_dz1), "sc_dphi1_v_dz1_L1_L3");
  671. auto sc_dz2_v_dz1_L1_L3 = fv::tup_apply(sc_cross_correlations,
  672. fv::tuple(BPIX_L1, BPIX_L3, v_dz2, v_dz1), "sc_dz2_v_dz1_L1_L3");
  673. auto sc_dphi2_v_dphi1_L1_L3 = fv::tup_apply(sc_cross_correlations,
  674. fv::tuple(BPIX_L1, BPIX_L3, v_dPhi2, v_dPhi1), "sc_dphi1_v_dphi1_L1_L3");
  675. auto sc_dz2_v_dphi1_L1_L3 = fv::tup_apply(sc_cross_correlations,
  676. fv::tuple(BPIX_L1, BPIX_L3, v_dz2, v_dPhi1), "sc_dz2_v_dphi1_L1_L3");
  677. hist_params = TH2Params::lookup("sc_dphi2_v_dz1_L1_L2");
  678. tds.register_container<ContainerTH2Many<float>>("sc_dphi2_v_dz1_L1_L3",
  679. sc_dphi2_v_dz1_L1_L3,
  680. "sc_dphi2_v_dz1_L1_L3", hist_params);
  681. hist_params = TH2Params::lookup("sc_dz2_v_dz1_L1_L2");
  682. tds.register_container<ContainerTH2Many<float>>("sc_dz2_v_dz1_L1_L3",
  683. sc_dz2_v_dz1_L1_L3,
  684. "sc_dz2_v_dz1_L1_L3", hist_params);
  685. hist_params = TH2Params::lookup("sc_dphi2_v_dphi1_L1_L2");
  686. tds.register_container<ContainerTH2Many<float>>("sc_dphi2_v_dphi1_L1_L3",
  687. sc_dphi2_v_dphi1_L1_L3,
  688. "sc_dphi2_v_dphi1_L1_L3", hist_params);
  689. hist_params = TH2Params::lookup("sc_dz2_v_dphi1_L1_L2");
  690. tds.register_container<ContainerTH2Many<float>>("sc_dz2_v_dphi1_L1_L3",
  691. sc_dz2_v_dphi1_L1_L3,
  692. "sc_dz2_v_dphi1_L1_L3", hist_params);
  693. hist_params = TH2Params::lookup("sc_dz2_v_dz1_L1_L2_zoom");
  694. tds.register_container<ContainerTH2Many<float>>("sc_dz2_v_dz1_L1_L3_zoom",
  695. sc_dz2_v_dz1_L1_L3,
  696. "sc_dz2_v_dz1_L1_L3_zoom", hist_params);
  697. hist_params = TH2Params::lookup("sc_dz2_v_dphi1_L1_L2_zoom");
  698. tds.register_container<ContainerTH2Many<float>>("sc_dz2_v_dphi1_L1_L3_zoom",
  699. sc_dz2_v_dphi1_L1_L3,
  700. "sc_dz2_v_dphi1_L1_L3_zoom", hist_params);
  701. hist_params = TH2Params::lookup("sc_dphi2_v_dphi1_L1_L2_zoom");
  702. tds.register_container<ContainerTH2Many<float>>("sc_dphi2_v_dphi1_L1_L3_zoom",
  703. sc_dphi2_v_dphi1_L1_L3,
  704. "sc_dphi2_v_dphi1_L1_L3_zoom", hist_params);
  705. }
  706. void run_analysis(const vector<fv::util::DataFileDescriptor>& dfds, const string output_filename, bool silent){
  707. gSystem->Load("libfilval.so");
  708. auto replace_suffix = [](const std::string& input, const std::string& new_suffix){
  709. return input.substr(0, input.find_last_of(".")) + new_suffix;
  710. };
  711. string log_filename = replace_suffix(output_filename, ".log");
  712. fv::util::Log::init_logger(log_filename, fv::util::LogPriority::kLogInfo);
  713. TrackingDataSet tds(output_filename, dfds, "trackingNtuple/tree");
  714. /* tds.set_max_events(10); */
  715. register_objects(tds);
  716. setup_hit_pair_functions();
  717. setup_first_hit_pairs(tds);
  718. setup_skipped_layer_hit_pairs(tds);
  719. /* setup_sc_residuals(tds); */
  720. /* setup_residual_cross_corrolations(tds); */
  721. tds.process(silent);
  722. tds.save_all();
  723. }
  724. int main(int argc, char * argv[]){
  725. fv::util::ArgParser args(argc, argv);
  726. bool silent = args.cmdOptionExists("-s");
  727. string output_filename = args.cmdOptionExists("-o") ? args.getCmdOption("-o") : "output.root";
  728. if(args.cmdOptionExists("-c")){
  729. fv::util::init_config(args.getCmdOption("-c"));
  730. }
  731. if(args.cmdOptionExists("-F")){
  732. auto file_list = fv::util::read_input_list(args.getCmdOption("-F"));
  733. run_analysis(file_list, output_filename, silent);
  734. }else if(args.cmdOptionExists("-f")){
  735. string input_filename = args.getCmdOption("-f");
  736. string data_label = args.cmdOptionExists("-l") ? args.getCmdOption("-l") : "";
  737. string data_category = args.cmdOptionExists("-c") ? args.getCmdOption("-c") : "";
  738. fv::util::DataFileDescriptor dfd(input_filename, data_label, data_category);
  739. run_analysis(std::vector<fv::util::DataFileDescriptor>({dfd}), output_filename, silent);
  740. }else {
  741. cout << "Usage: ./tracking_validation (-s) {-o output_filename} -F datafiles.txt" << endl;
  742. cout << " ./tracking_validation (-s) {-l DATA_LABEL} {-c DATA_CATEGORY} {-o output_filename} -f treefile.root" << endl;
  743. }
  744. return 0;
  745. }
  746. /* int main(int argc, char * argv[]){ */
  747. /* using namespace fv::util; */
  748. /* init_config(argv[1]); */
  749. /* int max_events = the_config.get("max-events").as<int>(); */
  750. /* cout << max_events << endl; */
  751. /* } */