123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238 |
- #include <iostream>
- #include <vector>
- #include <cmath>
- #include <TSystem.h>
- #include "filval.hpp"
- #include "root_filval.hpp"
- #include <boost/range/combine.hpp>
- #include "analysis/TrackingNtupleObjs.hpp"
- SeedCollection seeds;
- SimTrackCollection sim_tracks;
- SimVertexCollection sim_vertices;
- TrackCollection gsf_tracks;
- void register_objects(TrackingDataSet& tds){
- seeds.init(tds);
- sim_tracks.init(tds);
- sim_vertices.init(tds);
- gsf_tracks.init(tds);
- }
- struct {
- Value<float>* x;
- Value<float>* y;
- Value<float>* z;
- Value<float>* sigma_x;
- Value<float>* sigma_y;
- Value<float>* sigma_z;
- } bs;
- bool in_lum_region(const SimVertex& vertex) {
- /* Here approximate the luminous region as a cylinder with length 5*\sigma_z and radius
- * sqrt(\simga_x^2 + \sigma_y^2) centered as beamspot position.
- */
- float dx = vertex.x() - bs.x->get();
- float dy = vertex.y() - bs.y->get();
- float dz = vertex.z() - bs.z->get();
- auto radius = static_cast<float>(5 * sqrt(pow(bs.sigma_x->get(), 2) + pow(bs.sigma_y->get(), 2)));
- float half_len = 5*bs.sigma_z->get();
- return sqrt(dx*dx + dy*dy) < radius and abs(dz) < half_len;
- };
- bool is_good_sim(const SimTrack& sim_track) {
- const auto& vertex = sim_vertices[sim_track.parentVtxIdx()];
- return abs(sim_track.pdgId()) == 11 and in_lum_region(vertex);
- };
- void run(bool silent){
- using namespace std;
- using namespace fv;
- using namespace fv_root;
- auto file_list = the_config->get_source_files();
- string output_filename = the_config->get_output_filename();
- TrackingDataSet tds(output_filename, file_list, "trackingNtuple/tree");
- register_objects(tds);
- bs = {
- tds.track_branch<float>("bsp_x"),
- tds.track_branch<float>("bsp_y"),
- tds.track_branch<float>("bsp_z"),
- tds.track_branch<float>("bsp_sigmax"),
- tds.track_branch<float>("bsp_sigmay"),
- tds.track_branch<float>("bsp_sigmaz")
- };
- std::map<std::pair<int, int>, ContainerTH2<float>*> BPIX_residuals_dRz;
- std::map<std::pair<int, int>, ContainerTH2<float>*> BPIX_residuals_dPhi;
- std::map<std::pair<int, int>, ContainerTH2<float>*> FPIX_residuals_dRz;
- std::map<std::pair<int, int>, ContainerTH2<float>*> FPIX_residuals_dPhi;
- stringstream name;
- TH2Params hist_params;
- for (int layer=1; layer<=4; layer++) {
- for (int hit=1; hit<=3; hit++) {
- hist_params = (hit==1) ? TH2Params::lookup("dRz_v_Et") : TH2Params::lookup("dRz_v_Et_outer_hits");
- name.str("");
- name << "dRz_BPIX_L" << layer << "_H" << hit << "_v_Et";
- BPIX_residuals_dRz[{layer,hit}] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
- name.str("");
- name << "dRz_FPIX_L" << layer << "_H" << hit << "_v_Et";
- FPIX_residuals_dRz[{layer,hit}] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
- hist_params = (hit==1) ? TH2Params::lookup("dPhi_v_Et") : TH2Params::lookup("dPhi_v_Et_outer_hits");
- name.str("");
- name << "dPhi_BPIX_L" << layer << "_H" << hit << "_v_Et";
- BPIX_residuals_dPhi[{layer,hit}] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
- name.str("");
- name << "dPhi_FPIX_L" << layer << "_H" << hit << "_v_Et";
- FPIX_residuals_dPhi[{layer,hit}] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
- }
- }
- auto& seed_eff_v_pt = *tds.register_container<EfficiencyContainer<float>>("seed_eff_v_pt", TH1Params::lookup("eff_v_pt"));
- auto& seed_eff_v_eta = *tds.register_container<EfficiencyContainer<float>>("seed_eff_v_eta", TH1Params::lookup("eff_v_eta"));
- auto& seed_eff_v_phi = *tds.register_container<EfficiencyContainer<float>>("seed_eff_v_phi", TH1Params::lookup("eff_v_phi"));
- auto& tracking_eff_v_pt = *tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_pt", TH1Params::lookup("eff_v_pt"));
- auto& tracking_eff_v_eta = *tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_eta", TH1Params::lookup("eff_v_eta"));
- auto& tracking_eff_v_phi = *tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_phi", TH1Params::lookup("eff_v_phi"));
- auto& seed_pur_v_pt = *tds.register_container<EfficiencyContainer<float>>("seed_pur_v_pt", TH1Params::lookup("pur_v_pt"));
- auto& seed_pur_v_eta = *tds.register_container<EfficiencyContainer<float>>("seed_pur_v_eta", TH1Params::lookup("pur_v_eta"));
- auto& seed_pur_v_phi = *tds.register_container<EfficiencyContainer<float>>("seed_pur_v_phi", TH1Params::lookup("pur_v_phi"));
- auto& tracking_pur_v_pt = *tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_pt", TH1Params::lookup("pur_v_pt"));
- auto& tracking_pur_v_eta = *tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_eta", TH1Params::lookup("pur_v_eta"));
- auto& tracking_pur_v_phi = *tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_phi", TH1Params::lookup("pur_v_phi"));
- auto& hit_vs_layer_barrel = *tds.register_container<ContainerTH2<int>>("hit_vs_layer_barrel", TH2Params::lookup("hit_vs_layer"));
- auto& hit_vs_layer_forward = *tds.register_container<ContainerTH2<int>>("hit_vs_layer_forward", TH2Params::lookup("hit_vs_layer"));
- while (tds.next(!silent)) {
- for (const auto& sim_track : sim_tracks) {
- if (!is_good_sim(sim_track)) continue;
- if (seeds.size() == 0) continue;
- bool matched = !sim_track.seedIdx().empty();
- if (abs(sim_track.eta()) < 2.4)
- seed_eff_v_pt.fill(sim_track.pt(), matched);
- if (sim_track.pt() > 20.0)
- seed_eff_v_eta.fill(sim_track.eta(), matched);
- if (abs(sim_track.eta()) < 2.4 and sim_track.pt() > 20)
- seed_eff_v_phi.fill(sim_track.phi(), matched);
- }
- for (const auto& sim_track : sim_tracks) {
- if (!is_good_sim(sim_track)) continue;
- if (gsf_tracks.size() == 0) continue;
- bool matched = !sim_track.trkIdx().empty();
- if (abs(sim_track.eta()) < 2.4)
- tracking_eff_v_pt.fill(sim_track.pt(), matched);
- if (sim_track.pt() > 20.0)
- tracking_eff_v_eta.fill(sim_track.eta(), matched);
- if (abs(sim_track.eta()) < 2.4 and sim_track.pt() > 20)
- tracking_eff_v_phi.fill(sim_track.phi(), matched);
- }
- for (const auto& seed : seeds) {
- if (!seed.isECALDriven()) continue;
- bool match = false;
- for (const auto& sim_track_idx : seed.simTrkIdx()) {
- if(is_good_sim(sim_tracks[sim_track_idx])) {
- match = true;
- break;
- }
- }
- if (abs(seed.eta()) < 2.4)
- seed_pur_v_pt.fill(seed.pt(), match);
- if (seed.pt() > 20)
- seed_pur_v_eta.fill(seed.eta(), match);
- if (abs(seed.eta()) < 2.4 and seed.pt() > 20)
- seed_pur_v_phi.fill(seed.phi(), match);
- }
- for (const auto& gsf_track : gsf_tracks) {
- bool match = false;
- for (const auto& sim_track_idx : gsf_track.simTrkIdx()) {
- if(is_good_sim(sim_tracks[sim_track_idx])) {
- match = true;
- break;
- }
- }
- if (abs(gsf_track.eta()) < 2.4)
- tracking_pur_v_pt.fill(gsf_track.pt(), match);
- if (gsf_track.pt() > 20)
- tracking_pur_v_eta.fill(gsf_track.eta(), match);
- if (abs(gsf_track.eta()) < 2.4 and gsf_track.pt() > 20)
- tracking_pur_v_phi.fill(gsf_track.phi(), match);
- }
- for (const auto& seed : seeds) {
- if (seed.trkIdx() < 0) continue;
- if (!seed.isECALDriven()) continue;
- const auto the_track = gsf_tracks[seed.trkIdx()];
- vector<int> hits_valid;
- vector<float> hits_dRz;
- vector<float> hits_dPhi;
- if (the_track.q() == 1) {
- hits_valid = seed.isValidPos();
- hits_dRz = seed.dRZPos();
- hits_dPhi = seed.dPhiPos();
- } else {
- hits_valid = seed.isValidNeg();
- hits_dRz = seed.dRZNeg();
- hits_dPhi = seed.dPhiNeg();
- }
- const vector<int>& hits_isBarrel = seed.isBarrel();
- const vector<int>& hits_layerOrDiskNr = seed.layerOrDiskNr();
- size_t nHits = hits_valid.size();
- for(size_t hit_idx=0; hit_idx<nHits; hit_idx++) {
- if (!hits_valid[hit_idx]) continue;
- bool isBarrel = hits_isBarrel[hit_idx] == 1;
- int layerOrDiskNr = hits_layerOrDiskNr[hit_idx];
- float dRz = abs(hits_dRz[hit_idx]);
- float dPhi = abs(hits_dPhi[hit_idx]);
- if (isBarrel) {
- BPIX_residuals_dRz[{layerOrDiskNr, hit_idx+1}]->fill(dRz, seed.Et());
- BPIX_residuals_dPhi[{layerOrDiskNr, hit_idx+1}]->fill(dPhi, seed.Et());
- } else {
- FPIX_residuals_dRz[{layerOrDiskNr, hit_idx+1}]->fill(dRz, seed.Et());
- FPIX_residuals_dPhi[{layerOrDiskNr, hit_idx+1}]->fill(dPhi, seed.Et());
- }
- if (isBarrel)
- hit_vs_layer_barrel.fill(layerOrDiskNr, static_cast<const int &>(hit_idx + 1));
- else
- hit_vs_layer_forward.fill(layerOrDiskNr, static_cast<const int &>(hit_idx + 1));
- }
- }
- }
- tds.save_all();
- }
- int main(int argc, char * argv[]){
- using namespace fv_util;
- ArgParser args(argc, argv);
- bool silent = args.cmdOptionExists("-s");
- if(args.cmdOptionExists("-c")) {
- init_config(args.getCmdOption("-c"));
- init_log(LogPriority::kLogInfo);
- // gSystem->Load("libfilval.so");
- run(silent);
- } else {
- cout << "Usage: ./" << argv[0] << " (-s) -c config_file.yaml" << endl;
- }
- return 0;
- }
|