Adds resolution plots for innermost hits on seeds

analysis/tracking_validation.cpp

@@ -21,12 +21,16 @@ using namespace fv::root;
 
 using namespace std;
 
+typedef std::tuple<Seed, PixRecHit, Track, SimHit, SimTrack> MatchedTrack;
+typedef std::pair<std::vector<float>,std::vector<float>> pair_vec;
+
 #define HIT_TYPE_PIXEL 0
 #define HIT_TYPE_GLUED 1
 #define HIT_TYPE_STRIP 2
 
 #define PIXEL_BARREL 1
 #define PIXEL_ENDCAP 2
+vector<string> subdet_names = {"", "BPIX", "FPIX"};
 
 template<typename A, typename B>
 float displacement(const A& a, const B& b){
@@ -35,8 +39,12 @@ float displacement(const A& a, const B& b){
                      pow(a.x-b.x, 2));
 }
 
-typedef std::tuple<Seed, PixRecHit, Track, SimHit, SimTrack> MatchedTrack;
-typedef std::pair<std::vector<float>,std::vector<float>> pair_vec;
+bool in_det(const MatchedTrack &mt, int &&det, int &&layer){
+    auto& hit = std::get<PixRecHit>(mt);
+    return hit.det == det && hit.lay == layer;
+};
+
+
 
 vector<string> seedTypes =
     {"initialStepSeeds",
@@ -46,7 +54,6 @@ vector<string> seedTypes =
      "tripletElectronSeeds",
      "pixelPairElectronSeeds",
      "stripPairElectronSeeds"};
-int nSeedTypes = seedTypes.size();
 
 Value<vector<Seed>>* seeds;
 Value<vector<PixRecHit>>* pixrec_hits;
@@ -64,8 +71,57 @@ void register_objects(TrackingDataSet& tds){
     sim_tracks = register_sim_tracks(tds);
 }
 
+void setup_first_hit_pairs(TrackingDataSet& tds){
+    // Finds pairs of (rechit, simhit) where the rechit is the innermost hit on
+    // a seed in a particular barrel layer.
+    typedef std::tuple<PixRecHit, SimHit> HitPair;
+    auto& matched_hits_in_layer =
+    func<vector<HitPair>(vector<Seed>,
+                         vector<PixRecHit>,
+                         vector<SimHit>, int)>("find_matched_innermost_hits_in_layer",
+        FUNC(([](const vector<Seed>& seeds,
+                 const vector<PixRecHit>& pixrec_hits,
+                 const vector<SimHit>& sim_hits,
+                 const int&& bpix_layer){
+            vector<HitPair> matched_hits;
+            for(const Seed &seed : seeds){
+                if( seed.algoOriginal < 0 || seed.algoOriginal >= seedTypes.size()) continue;
+                for(auto tup : boost::combine(seed.hitIdx, seed.hitType)){ //for hits in track's seed
+                    int hitIdx, hitType;
+                    boost::tie(hitIdx, hitType) = tup;
+                    if(hitType != HIT_TYPE_PIXEL) continue; // take only pixel hits for now
+                    const PixRecHit &rec_hit = pixrec_hits[hitIdx];
+                    if(rec_hit.det == PIXEL_BARREL && rec_hit.lay == bpix_layer){
+                        if(rec_hit.simHitIdx.size() > 0){
+                            matched_hits.push_back({rec_hit, sim_hits[rec_hit.simHitIdx[0]]});
+                        }
+                    }
+                    break;
+                }
+            }
+            return matched_hits;
+        })));
+
+    auto first_hits_in_B1 = fv::apply(matched_hits_in_layer, fv::tuple(seeds, pixrec_hits, sim_hits, constant("1", 1)), "first_hits_in_B1");
+    auto first_hits_in_B2 = fv::apply(matched_hits_in_layer, fv::tuple(seeds, pixrec_hits, sim_hits, constant("2", 2)), "first_hits_in_B2");
+
+    auto& calc_dphi = func<float(HitPair)>("calc_dphi",
+        FUNC(([](const HitPair& hit_pair){
+            const auto &rec_hit = std::get<PixRecHit>(hit_pair);
+            const auto &sim_hit = std::get<SimHit>(hit_pair);
+            return atan2(rec_hit.x, rec_hit.y) - atan2(sim_hit.x, sim_hit.y);
+        })));
+
+    tds.register_container<ContainerTH1Many<float>>("matched_hits_from_B1", "Matched Hits dPhi - B1",
+                                                    fv::map(calc_dphi, first_hits_in_B1),
+                                                    100, -0.001, 0.001,"dPhi(rad)");
+    tds.register_container<ContainerTH1Many<float>>("matched_hits_from_B2", "Matched Hits dPhi - B2",
+                                                    fv::map(calc_dphi, first_hits_in_B2),
+                                                    100, -0.001, 0.001,"dPhi(rad)");
+}
+
 void setup_matched_tracks(TrackingDataSet& tds){
-    // Finds pairs of (track,sim_track) where the cycle
+    // Finds sets of (seed, rechit, track, simhit sim_track) where the cycle
     // track->seed->rec_hit->sim_hit->sim_track->track
     // is closed
     auto& find_matched_tracks =
@@ -85,7 +141,7 @@ void setup_matched_tracks(TrackingDataSet& tds){
             vector<MatchedTrack> matched_tracks;
             for(const Track &track : tracks){
                 const Seed& seed = seeds[track.seedIdx];
-                if( seed.algoOriginal < 0 || seed.algoOriginal >= nSeedTypes) continue;
+                if( seed.algoOriginal < 0 || seed.algoOriginal >= seedTypes.size()) continue;
                 INFO(boost::format("       seed indx= %d   algorithm= %d   type %s\n")
                         % track.seedIdx % seed.algoOriginal % seedTypes[seed.algoOriginal]);
                 for(auto tup : boost::combine(seed.hitIdx, seed.hitType)){ //for hits in track's seed
@@ -97,7 +153,8 @@ void setup_matched_tracks(TrackingDataSet& tds){
                     TVector3 recoHitPoint;
                     recoHitPoint.SetXYZ( rec_hit.x, rec_hit.y, rec_hit.z);
                     float Rreco = recoHitPoint.Perp();
-                    INFO(boost::format("            hit= %d is at radius= %f\n") % hitIdx % Rreco);
+                    INFO(boost::format("            hit= %d is in %s-L%d at radius= %f\n")
+                                       % hitIdx % subdet_names[rec_hit.det] % rec_hit.lay % Rreco);
 
                     for(const int& simHitIdx : rec_hit.simHitIdx){ // for sim-hits matching rec-hit
                         const SimHit &sim_hit = sim_hits[simHitIdx];
@@ -117,11 +174,6 @@ void setup_matched_tracks(TrackingDataSet& tds){
     fv::apply(find_matched_tracks, fv::tuple(seeds, pixrec_hits, tracks, sim_hits, sim_tracks), "matched_tracks");
 }
 
-bool in_det(const MatchedTrack &mt, int &&det, int &&layer){
-    auto& hit = std::get<PixRecHit>(mt);
-    return hit.det == det && hit.lay == layer;
-};
-
 void setup_residuals(TrackingDataSet &tds){
 
 
@@ -264,15 +316,29 @@ void setup_residuals(TrackingDataSet &tds){
     tds.register_container<ContainerTH2Many<float>>("matched_track_dzs_v_eta_B4", "Matched Track dz - B4",
                                                     fv::apply(calc_dz_v_eta, fv::tuple(matched_tracks_B4)),
                                                     100, -4, 4, 100, -.002, .002,"Eta","dz(cm)");
-    tds.register_container<ContainerTH2Many<float>>("matched_track_dzs_v_eta_F1", "Matched Track dz - F1",
-                                                    fv::apply(calc_dz_v_eta, fv::tuple(matched_tracks_F1)),
-                                                    100, -4, 4, 100, -.002, .002,"Eta","dz(cm)");
-    tds.register_container<ContainerTH2Many<float>>("matched_track_dzs_v_eta_F2", "Matched Track dz - F2",
-                                                    fv::apply(calc_dz_v_eta, fv::tuple(matched_tracks_F2)),
-                                                    100, -4, 4, 100, -.002, .002,"Eta","dz(cm)");
-    tds.register_container<ContainerTH2Many<float>>("matched_track_dzs_v_eta_F3", "Matched Track dz - F3",
-                                                    fv::apply(calc_dz_v_eta, fv::tuple(matched_tracks_F3)),
-                                                    100, -4, 4, 100, -.002, .002,"Eta","dz(cm)");
+    auto& calc_dr_v_eta = func<pair_vec(vector<MatchedTrack>)>("matched_track_hit_drs",
+        FUNC(([](const vector<MatchedTrack>& matched_tracks){
+            vector<float> drs;
+            vector<float> etas;
+            for(auto matched_track : matched_tracks){
+                auto& rec_hit = std::get<PixRecHit>(matched_track);
+                auto& sim_hit = std::get<SimHit>(matched_track);
+                float sim_hit_r = sqrt(pow(sim_hit.x, 2) + pow(sim_hit.y, 2));
+                float rec_hit_r = sqrt(pow(rec_hit.x, 2) + pow(rec_hit.y, 2));
+                drs.push_back(sim_hit_r - rec_hit_r);
+                etas.push_back(std::get<SimTrack>(matched_track).eta);
+            }
+            return std::make_pair(etas, drs);
+        })));
+    tds.register_container<ContainerTH2Many<float>>("matched_track_drs_v_eta_F1", "Matched Track dz - F1",
+                                                    fv::apply(calc_dr_v_eta, fv::tuple(matched_tracks_F1)),
+                                                    100, -4, 4, 100, -.002, .002,"Eta","dr(cm)");
+    tds.register_container<ContainerTH2Many<float>>("matched_track_drs_v_eta_F2", "Matched Track dz - F2",
+                                                    fv::apply(calc_dr_v_eta, fv::tuple(matched_tracks_F2)),
+                                                    100, -4, 4, 100, -.002, .002,"Eta","dr(cm)");
+    tds.register_container<ContainerTH2Many<float>>("matched_track_drs_v_eta_F3", "Matched Track dz - F3",
+                                                    fv::apply(calc_dr_v_eta, fv::tuple(matched_tracks_F3)),
+                                                    100, -4, 4, 100, -.002, .002,"Eta","dr(cm)");
 }
 
 void run_analysis(const vector<fv::util::DataFileDescriptor>& dfds, const string output_filename, bool silent){
@@ -290,13 +356,13 @@ void run_analysis(const vector<fv::util::DataFileDescriptor>& dfds, const string
 
     setup_matched_tracks(tds);
     setup_residuals(tds);
+    setup_first_hit_pairs(tds);
 
     tds.process(silent);
     tds.save_all();
 }
 
-int main(int argc, char * argv[])
-{
+int main(int argc, char * argv[]){
     fv::util::ArgParser args(argc, argv);
     bool silent = args.cmdOptionExists("-s");
     string output_filename = args.cmdOptionExists("-o") ? args.getCmdOption("-o") : "output.root";

filval

@@ -1 +1 @@
-Subproject commit d7237e5c01868713d6b8d7491a28ffbebe2bbd8d
+Subproject commit f455a5d61f96309e0ea8cf526916a06a62141f20