Adds ability to specify multiple input root files via a text file

containing line delimited paths. Plots are work-in-progress.

.gitignore

@@ -2,3 +2,5 @@
 legacy/
 data/
 build/
+
+tags

analysis/obj_types.cpp

@@ -15,6 +15,7 @@ typedef TreeDataSet<TrackingNtuple> TrackingDataSet;
 //==============PIXEL REC HIT=================================================
 //============================================================================
 struct PixRecHit {
+    unsigned int   idx;
     unsigned short det;
     unsigned short lay;
     unsigned int   detId;
@@ -78,26 +79,27 @@ register_pixrec_hits(TrackingDataSet &tds){
                  const vector<vector<float>>& chargeFraction){
             std::vector<PixRecHit> pixrec_hits;
             for(int i=0; i<det.size(); i++)
-                pixrec_hits.push_back({det[i],
-                                   lay[i],
-                                   detId[i],
-                                   simType[i],
-                                   x[i],
-                                   y[i],
-                                   z[i],
-                                   xx[i],
-                                   xy[i],
-                                   yy[i],
-                                   yz[i],
-                                   zz[i],
-                                   zx[i],
-                                   radL[i],
-                                   bbxi[i],
-                                   trkIdx[i],
-                                   seeIdx[i],
-                                   simHitIdx[i],
-                                   chargeFraction[i]
-                                   });
+                pixrec_hits.push_back({i,
+                                       det[i],
+                                       lay[i],
+                                       detId[i],
+                                       simType[i],
+                                       x[i],
+                                       y[i],
+                                       z[i],
+                                       xx[i],
+                                       xy[i],
+                                       yy[i],
+                                       yz[i],
+                                       zz[i],
+                                       zx[i],
+                                       radL[i],
+                                       bbxi[i],
+                                       trkIdx[i],
+                                       seeIdx[i],
+                                       simHitIdx[i],
+                                       chargeFraction[i]
+                                       });
             return pixrec_hits;
         })));
 
@@ -129,19 +131,20 @@ register_pixrec_hits(TrackingDataSet &tds){
 //===================SIM HIT==================================================
 //============================================================================
 struct SimHit {
-   unsigned short det;
-   unsigned short lay;
-   unsigned int   detId;
-   float          x;
-   float          y;
-   float          z;
-   int            particle;
-   short          process;
-   float          eloss;
-   float          tof;
-   int            simTrkIdx;
-   vector<int>    hitIdx;
-   vector<int>    hitType;
+    unsigned int   idx;
+    unsigned short det;
+    unsigned short lay;
+    unsigned int   detId;
+    float          x;
+    float          y;
+    float          z;
+    int            particle;
+    short          process;
+    float          eloss;
+    float          tof;
+    int            simTrkIdx;
+    vector<int>    hitIdx;
+    vector<int>    hitType;
 };
 
 Value<vector<SimHit>>*
@@ -177,7 +180,8 @@ register_sim_hits(TrackingDataSet &tds){
                  ){
             std::vector<SimHit> sim_hits;
             for(int i=0; i<det.size(); i++)
-                sim_hits.push_back({det[i],
+                sim_hits.push_back({i,
+                                    det[i],
                                     lay[i],
                                     detId[i],
                                     x[i],
@@ -215,54 +219,55 @@ register_sim_hits(TrackingDataSet &tds){
 //=======================SEED=================================================
 //============================================================================
 struct Seed {
-   short          fitok;
-   float          px;
-   float          py;
-   float          pz;
-   float          pt;
-   float          eta;
-   float          phi;
-   float          dxy;
-   float          dz;
-   float          ptErr;
-   float          etaErr;
-   float          phiErr;
-   float          dxyErr;
-   float          dzErr;
-   float          chi2;
-   int            q;
-   unsigned int   nValid;
-   unsigned int   nPixel;
-   unsigned int   nGlued;
-   unsigned int   nStrip;
-   unsigned int   algo;
-   unsigned int   algoOriginal;
-   int            trkIdx;
-   vector<float>  shareFrac;
-   vector<int>    simTrkIdx;
-   vector<int>    hitIdx;
-   vector<int>    hitType;
-   unsigned int   offset;
-   unsigned char  hitsMask;
-   int            subDet2;
-   float          dRz2;
-   float          dPhi2;
-   float          dRz2Pos;
-   float          dPhi2Pos;
-   int            subDet1;
-   float          dRz1;
-   float          dPhi1;
-   float          dRz1Pos;
-   float          dPhi1Pos;
-   float          hoe1;
-   float          hoe2;
-   float          superClusterEnergy;
-   float          superClusterEta;
-   float          superClusterPhi;
-   float          superClusterEt;
-   int            superClusterIdx;
-   unsigned int   ecalDriven;
-   unsigned int   trkDriven;
+    unsigned int   idx;
+    short          fitok;
+    float          px;
+    float          py;
+    float          pz;
+    float          pt;
+    float          eta;
+    float          phi;
+    float          dxy;
+    float          dz;
+    float          ptErr;
+    float          etaErr;
+    float          phiErr;
+    float          dxyErr;
+    float          dzErr;
+    float          chi2;
+    int            q;
+    unsigned int   nValid;
+    unsigned int   nPixel;
+    unsigned int   nGlued;
+    unsigned int   nStrip;
+    unsigned int   algo;
+    unsigned int   algoOriginal;
+    int            trkIdx;
+    vector<float>  shareFrac;
+    vector<int>    simTrkIdx;
+    vector<int>    hitIdx;
+    vector<int>    hitType;
+    unsigned int   offset;
+    unsigned char  hitsMask;
+    int            subDet2;
+    float          dRz2;
+    float          dPhi2;
+    float          dRz2Pos;
+    float          dPhi2Pos;
+    int            subDet1;
+    float          dRz1;
+    float          dPhi1;
+    float          dRz1Pos;
+    float          dPhi1Pos;
+    float          hoe1;
+    float          hoe2;
+    float          superClusterEnergy;
+    float          superClusterEta;
+    float          superClusterPhi;
+    float          superClusterEt;
+    int            superClusterIdx;
+    unsigned int   ecalDriven;
+    unsigned int   trkDriven;
 };
 
 
@@ -367,7 +372,8 @@ register_seeds(TrackingDataSet &tds){
                const vector<unsigned int>& trkDriven){
             std::vector<Seed> seeds;
             for(int i=0; i<fitok.size(); i++)
-                seeds.push_back({fitok[i],
+                seeds.push_back({i,
+                                 fitok[i],
                                  px[i],
                                  py[i],
                                  pz[i],
@@ -476,59 +482,60 @@ register_seeds(TrackingDataSet &tds){
 //=======================TRACK================================================
 //============================================================================
 struct Track{
-   float         px;
-   float         py;
-   float         pz;
-   float         pt;
-   float         inner_px;
-   float         inner_py;
-   float         inner_pz;
-   float         inner_pt;
-   float         outer_px;
-   float         outer_py;
-   float         outer_pz;
-   float         outer_pt;
-   float         eta;
-   float         lambda;
-   float         cotTheta;
-   float         phi;
-   float         dxy;
-   float         dz;
-   float         ptErr;
-   float         etaErr;
-   float         lambdaErr;
-   float         phiErr;
-   float         dxyErr;
-   float         dzErr;
-   float         refpoint_x;
-   float         refpoint_y;
-   float         refpoint_z;
-   float         nChi2;
-   int           q;
-   unsigned int  nValid;
-   unsigned int  nInvalid;
-   unsigned int  nPixel;
-   unsigned int  nStrip;
-   unsigned int  nPixelLay;
-   unsigned int  nStripLay;
-   unsigned int  n3DLay;
-   unsigned int  nOuterLost;
-   unsigned int  nInnerLost;
-   unsigned int  algo;
-   unsigned int  originalAlgo;
-   ULong64_t     algoMask;
-   unsigned int  stopReason;
-   short         isHP;
-   int           seedIdx;
-   float         vtxx;
-   float         vtxy;
-   float         vtxz;
-   vector<float> shareFrac;
-   vector<int>   simTrkIdx;
-   vector<int>   hitIdx;
-   vector<int>   hitType;
-   int           genIdx;
-   float         genDR;
+    unsigned int   idx;
+    float         px;
+    float         py;
+    float         pz;
+    float         pt;
+    float         inner_px;
+    float         inner_py;
+    float         inner_pz;
+    float         inner_pt;
+    float         outer_px;
+    float         outer_py;
+    float         outer_pz;
+    float         outer_pt;
+    float         eta;
+    float         lambda;
+    float         cotTheta;
+    float         phi;
+    float         dxy;
+    float         dz;
+    float         ptErr;
+    float         etaErr;
+    float         lambdaErr;
+    float         phiErr;
+    float         dxyErr;
+    float         dzErr;
+    float         refpoint_x;
+    float         refpoint_y;
+    float         refpoint_z;
+    float         nChi2;
+    int           q;
+    unsigned int  nValid;
+    unsigned int  nInvalid;
+    unsigned int  nPixel;
+    unsigned int  nStrip;
+    unsigned int  nPixelLay;
+    unsigned int  nStripLay;
+    unsigned int  n3DLay;
+    unsigned int  nOuterLost;
+    unsigned int  nInnerLost;
+    unsigned int  algo;
+    unsigned int  originalAlgo;
+    ULong64_t     algoMask;
+    unsigned int  stopReason;
+    short         isHP;
+    int           seedIdx;
+    float         vtxx;
+    float         vtxy;
+    float         vtxz;
+    vector<float> shareFrac;
+    vector<int>   simTrkIdx;
+    vector<int>   hitIdx;
+    vector<int>   hitType;
+    int           genIdx;
+    float         genDR;
 };
 
 Value<vector<Track>>*
@@ -641,7 +648,8 @@ register_tracks(TrackingDataSet &tds){
                                     const vector<float>&         genDR){
             std::vector<Track> tracks;
             for(int i=0; i<px.size(); i++)
-                tracks.push_back({px[i],
+                tracks.push_back({i,
+                                  px[i],
                                   py[i],
                                   pz[i],
                                   pt[i],
@@ -758,34 +766,35 @@ register_tracks(TrackingDataSet &tds){
 //===================SIM TRACK================================================
 //============================================================================
 struct SimTrack{
-   int           event;
-   int           bunchCrossing;
-   int           pdgId;
-   float         px;
-   float         py;
-   float         pz;
-   float         pt;
-   float         eta;
-   float         phi;
-   float         pca_pt;
-   float         pca_eta;
-   float         pca_lambda;
-   float         pca_cotTheta;
-   float         pca_phi;
-   float         pca_dxy;
-   float         pca_dz;
-   int           q;
-   unsigned int  nValid;
-   unsigned int  nPixel;
-   unsigned int  nStrip;
-   unsigned int  nLay;
-   unsigned int  nPixelLay;
-   unsigned int  n3DLay;
-   vector<int>   trkIdx;
-   vector<float> shareFrac;
-   int           parentVtxIdx;
-   vector<int>   decayVtxIdx;
-   vector<int>   simHitIdx;
+    unsigned int   idx;
+    int           event;
+    int           bunchCrossing;
+    int           pdgId;
+    float         px;
+    float         py;
+    float         pz;
+    float         pt;
+    float         eta;
+    float         phi;
+    float         pca_pt;
+    float         pca_eta;
+    float         pca_lambda;
+    float         pca_cotTheta;
+    float         pca_phi;
+    float         pca_dxy;
+    float         pca_dz;
+    int           q;
+    unsigned int  nValid;
+    unsigned int  nPixel;
+    unsigned int  nStrip;
+    unsigned int  nLay;
+    unsigned int  nPixelLay;
+    unsigned int  n3DLay;
+    vector<int>   trkIdx;
+    vector<float> shareFrac;
+    int           parentVtxIdx;
+    vector<int>   decayVtxIdx;
+    vector<int>   simHitIdx;
 };
 
 
@@ -850,36 +859,35 @@ register_sim_tracks(TrackingDataSet &tds){
                const vector<vector<int>>&   simHitIdx){
             std::vector<SimTrack> sim_tracks;
             for(int i=0; i<event.size(); i++)
-                sim_tracks.push_back({
-                   event[i],
-                   bunchCrossing[i],
-                   pdgId[i],
-                   px[i],
-                   py[i],
-                   pz[i],
-                   pt[i],
-                   eta[i],
-                   phi[i],
-                   pca_pt[i],
-                   pca_eta[i],
-                   pca_lambda[i],
-                   pca_cotTheta[i],
-                   pca_phi[i],
-                   pca_dxy[i],
-                   pca_dz[i],
-                   q[i],
-                   nValid[i],
-                   nPixel[i],
-                   nStrip[i],
-                   nLay[i],
-                   nPixelLay[i],
-                   n3DLay[i],
-                   trkIdx[i],
-                   shareFrac[i],
-                   parentVtxIdx[i],
-                   decayVtxIdx[i],
-                   simHitIdx[i]
-                });
+                sim_tracks.push_back({i,
+                                      event[i],
+                                      bunchCrossing[i],
+                                      pdgId[i],
+                                      px[i],
+                                      py[i],
+                                      pz[i],
+                                      pt[i],
+                                      eta[i],
+                                      phi[i],
+                                      pca_pt[i],
+                                      pca_eta[i],
+                                      pca_lambda[i],
+                                      pca_cotTheta[i],
+                                      pca_phi[i],
+                                      pca_dxy[i],
+                                      pca_dz[i],
+                                      q[i],
+                                      nValid[i],
+                                      nPixel[i],
+                                      nStrip[i],
+                                      nLay[i],
+                                      nPixelLay[i],
+                                      n3DLay[i],
+                                      trkIdx[i],
+                                      shareFrac[i],
+                                      parentVtxIdx[i],
+                                      decayVtxIdx[i],
+                                      simHitIdx[i]});
             return sim_tracks;
         })));
     return apply(builder,

analysis/tracking_validation.cpp

@@ -11,6 +11,7 @@
 #include "filval/root/filval.hpp"
 
 #include <boost/range/combine.hpp>
+#include <boost/format.hpp>
 
 #include "TrackingNtuple.h"
 #include "analysis/obj_types.cpp"
@@ -25,6 +26,16 @@ using namespace std;
 #define HIT_TYPE_GLUED 1
 #define HIT_TYPE_STRIP 2
 
+vector<string> seedTypes =
+    {"initialStepSeeds",
+     "highPtTripletStepSeeds",
+     "mixedTripletStepSeeds",
+     "pixelLessStepSeeds",
+     "tripletElectronSeeds",
+     "pixelPairElectronSeeds",
+     "stripPairElectronSeeds"};
+int nSeedTypes = seedTypes.size();
+
 Value<vector<Seed>>* seeds;
 Value<vector<PixRecHit>>* pixrec_hits;
 Value<vector<Track>>* tracks;
@@ -41,128 +52,165 @@ void register_objects(TrackingDataSet& tds){
     sim_tracks = register_sim_tracks(tds);
 }
 
-void setup_dphi_vs_eta(TrackingDataSet& tds){
-    auto& calc_dphi_vs_eta =
-    func<std::pair<vector<float>,vector<float>>(vector<Seed>,
-                                                vector<PixRecHit>,
-                                                vector<SimHit>, int)>("calc_dphi_vs_eta",
-            FUNC(([](const vector<Seed>& seeds,
-                     const vector<PixRecHit> pixrec_hits,
-                     const vector<SimHit> sim_hits,
-                     int layer){
-                vector<float> dphi;
-                vector<float> eta;
-                for(const Seed &seed : seeds){
+void setup_matched_tracks(TrackingDataSet& tds){
+    // Finds pairs of (track,sim_track) where the cycle
+    // track->seed->rec_hit->sim_hit->sim_track->track
+    // is closed
+    typedef std::tuple<Seed, PixRecHit, Track, SimHit, SimTrack> MatchedTrack;
+
+    auto& find_matched_tracks =
+    func<vector<MatchedTrack>(vector<Seed>,
+                              vector<PixRecHit>,
+                              vector<Track>,
+                              vector<SimHit>,
+                              vector<SimTrack>)>("find_matched_tracks",
+        FUNC(([](const vector<Seed>& seeds,
+                 const vector<PixRecHit> pixrec_hits,
+                 const vector<Track>& tracks,
+                 const vector<SimHit> sim_hits,
+                 const vector<SimTrack> sim_tracks){
+            INFO("New event");
+            INFO(boost::format("Number of tracks is %d, number of seeds is %d, number of hits %d, number of simhits %d")
+                    % tracks.size() % seeds.size() % pixrec_hits.size() % sim_hits.size());
+            vector<MatchedTrack> matched_tracks;
+            for(const Track &track : tracks){
+                const Seed& seed = seeds[track.seedIdx];
+                if( seed.algoOriginal < 0 || seed.algoOriginal >= nSeedTypes) 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
                     int hitIdx, hitType;
-                    for(auto tup : boost::combine(seed.hitIdx, seed.hitType)){
-                        boost::tie(hitIdx, hitType) = tup;
-                        if(hitType == HIT_TYPE_PIXEL){
-                            const PixRecHit &hit = pixrec_hits[hitIdx];
-                            if(layer && hit.lay != layer) continue;
-                            float hit_phi = atan2(hit.y, hit.x);
-                            for(const int& simHitIdx : hit.simHitIdx){
-                                const SimHit &sim_hit = sim_hits[simHitIdx];
-                                dphi.push_back(hit_phi-atan2(sim_hit.y, sim_hit.x));
-                                eta.push_back(seed.eta);
+                    boost::tie(hitIdx, hitType) = tup;
+                    if(hitType != HIT_TYPE_PIXEL) continue; // take only pixel hits for now
+                    const PixRecHit &rec_hit = pixrec_hits[hitIdx];
+
+                    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);
+
+                    for(const int& simHitIdx : rec_hit.simHitIdx){ // for sim-hits matching rec-hit
+                        const SimHit &sim_hit = sim_hits[simHitIdx];
+                        const SimTrack &sim_track = sim_tracks[sim_hit.simTrkIdx];
+                        for(const int &trkIdx : sim_track.trkIdx){ // for tracks matched to this sim_track
+                            if(trkIdx == track.idx){  // sim_track matches with our original track
+                                matched_tracks.push_back({seed, rec_hit, track, sim_hit, sim_track});
+                                break;
                             }
                         }
                     }
                 }
-                return std::make_pair(dphi,eta);
-            })));
-
-    auto dphi_vs_eta_B1 =
-    fv::apply(calc_dphi_vs_eta, fv::tuple(seeds, pixrec_hits, sim_hits, constant<int>("1", 1)), "dphi_vs_eta_B1");
-    auto dphi_vs_eta_B2 =
-    fv::apply(calc_dphi_vs_eta, fv::tuple(seeds, pixrec_hits, sim_hits, constant<int>("2", 2)), "dphi_vs_eta_B2");
-    auto dphi_vs_eta_B3 =
-    fv::apply(calc_dphi_vs_eta, fv::tuple(seeds, pixrec_hits, sim_hits, constant<int>("3", 3)), "dphi_vs_eta_B3");
-    auto dphi_vs_eta_B4 =
-    fv::apply(calc_dphi_vs_eta, fv::tuple(seeds, pixrec_hits, sim_hits, constant<int>("4", 4)), "dphi_vs_eta_B4");
-
-
-    tds.register_container<ContainerTH2Many<float>>("dphi_vs_eta_B1", "dphi vs eta (BPIX L1)",
-                                                      dphi_vs_eta_B1, 100, -0.002, 0.002, 20, -3, 3,
-                                                      "dphi",  "eta");
-    tds.register_container<ContainerTH2Many<float>>("dphi_vs_eta_B2", "dphi vs eta (BPIX L2)",
-                                                      dphi_vs_eta_B2, 100, -0.002, 0.002, 20, -3, 3,
-                                                      "dphi",  "eta");
-    tds.register_container<ContainerTH2Many<float>>("dphi_vs_eta_B3", "dphi vs eta (BPIX L3)",
-                                                      dphi_vs_eta_B3, 100, -0.002, 0.002, 20, -3, 3,
-                                                      "dphi",  "eta");
-    tds.register_container<ContainerTH2Many<float>>("dphi_vs_eta_B4", "dphi vs eta (BPIX L4)",
-                                                      dphi_vs_eta_B4, 100, -0.002, 0.002, 20, -3, 3,
-                                                      "dphi",  "eta");
+            }
+            return matched_tracks;
+        })));
+
+    auto matched_tracks =
+        fv::apply(find_matched_tracks, fv::tuple(seeds, pixrec_hits, tracks, sim_hits, sim_tracks), "matched_tracks");
+
+    auto delta_pt = fv::map(func<float(MatchedTrack)>("matched_track_delta_pt",
+        FUNC(([](const MatchedTrack& matched_track){
+            const Track& track = std::get<Track>(matched_track);
+            const SimTrack& sim_track = std::get<SimTrack>(matched_track);
+            return track.pt - sim_track.pt;
+        }))), matched_tracks, "matched_track_delta_pt");
+
+    tds.register_container<ContainerTH1Many<float>>("matched_track_delta_pt", "Matched Track delta P_t",
+                                                      delta_pt, 20, -10, 10, "P_t");
+    /* auto dphi_vs_eta_B1 = */
+    /*     fv::apply(calc_dphi_vs_eta, fv::tuple(seeds, pixrec_hits, sim_hits, constant<int>("1", 1)), "dphi_vs_eta_B1"); */
+    /* auto dphi_vs_eta_B2 = */
+    /*     fv::apply(calc_dphi_vs_eta, fv::tuple(seeds, pixrec_hits, sim_hits, constant<int>("2", 2)), "dphi_vs_eta_B2"); */
+    /* auto dphi_vs_eta_B3 = */
+    /*     fv::apply(calc_dphi_vs_eta, fv::tuple(seeds, pixrec_hits, sim_hits, constant<int>("3", 3)), "dphi_vs_eta_B3"); */
+    /* auto dphi_vs_eta_B4 = */
+    /*     fv::apply(calc_dphi_vs_eta, fv::tuple(seeds, pixrec_hits, sim_hits, constant<int>("4", 4)), "dphi_vs_eta_B4"); */
+
+
+    /* tds.register_container<ContainerTH2Many<float>>("dphi_vs_eta_B1", "dphi vs eta (BPIX L1)", */
+    /*                                                   dphi_vs_eta_B1, 20, -3, 3, 100, -0.002, 0.002, */
+    /*                                                   "eta",  "dphi"); */
+    /* tds.register_container<ContainerTH2Many<float>>("dphi_vs_eta_B2", "dphi vs eta (BPIX L2)", */
+    /*                                                   dphi_vs_eta_B2, 20, -3, 3, 100, -0.002, 0.002, */
+    /*                                                   "eta",  "dphi"); */
+    /* tds.register_container<ContainerTH2Many<float>>("dphi_vs_eta_B3", "dphi vs eta (BPIX L3)", */
+    /*                                                   dphi_vs_eta_B3, 20, -3, 3, 100, -0.002, 0.002, */
+    /*                                                   "eta",  "dphi"); */
+    /* tds.register_container<ContainerTH2Many<float>>("dphi_vs_eta_B4", "dphi vs eta (BPIX L4)", */
+    /*                                                   dphi_vs_eta_B4, 20, -3, 3, 100, -0.002, 0.002, */
+    /*                                                   "eta",  "dphi"); */
 
 }
 
-void setup_dz_vs_eta(TrackingDataSet& tds){
-    auto& calc_dz_vs_eta =
-    func<std::pair<vector<float>,vector<float>>(vector<Seed>,
-                                                vector<PixRecHit>,
-                                                vector<SimHit>, int)>("calc_dz_vs_eta",
-            FUNC(([](const vector<Seed>& seeds,
-                     const vector<PixRecHit> pixrec_hits,
-                     const vector<SimHit> sim_hits,
-                     int layer){
-                vector<float> dz;
-                vector<float> eta;
-                for(const Seed &seed : seeds){
-                    int hitIdx, hitType;
-                    for(auto tup : boost::combine(seed.hitIdx, seed.hitType)){
-                        boost::tie(hitIdx, hitType) = tup;
-                        if(hitType == HIT_TYPE_PIXEL){
-                            const PixRecHit &hit = pixrec_hits[hitIdx];
-                            if(layer && hit.lay != layer) continue;
-                            for(const int& simHitIdx : hit.simHitIdx){
-                                dz.push_back(hit.z-sim_hits[simHitIdx].z);
-                                eta.push_back(seed.eta);
-                            }
-                        }
-                    }
-                }
-                return std::make_pair(dz,eta);
-            })));
-
-    auto dz_vs_eta_B1 =
-    fv::apply(calc_dz_vs_eta, fv::tuple(seeds, pixrec_hits, sim_hits, constant<int>("1", 1)), "dz_vs_eta_B1");
-    auto dz_vs_eta_B2 =
-    fv::apply(calc_dz_vs_eta, fv::tuple(seeds, pixrec_hits, sim_hits, constant<int>("2", 2)), "dz_vs_eta_B2");
-    auto dz_vs_eta_B3 =
-    fv::apply(calc_dz_vs_eta, fv::tuple(seeds, pixrec_hits, sim_hits, constant<int>("3", 3)), "dz_vs_eta_B3");
-    auto dz_vs_eta_B4 =
-    fv::apply(calc_dz_vs_eta, fv::tuple(seeds, pixrec_hits, sim_hits, constant<int>("4", 4)), "dz_vs_eta_B4");
-
-
-    tds.register_container<ContainerTH2Many<float>>("dz_vs_eta_B1", "dz vs eta (BPIX L1)",
-                                                      dz_vs_eta_B1, 20, -0.01, .01, 20, -3, 3,
-                                                      "dz",  "eta");
-    tds.register_container<ContainerTH2Many<float>>("dz_vs_eta_B2", "dz vs eta (BPIX L2)",
-                                                      dz_vs_eta_B2, 20, -0.01, .01, 20, -3, 3,
-                                                      "dz",  "eta");
-    tds.register_container<ContainerTH2Many<float>>("dz_vs_eta_B3", "dz vs eta (BPIX L3)",
-                                                      dz_vs_eta_B3, 20, -0.01, .01, 20, -3, 3,
-                                                      "dz",  "eta");
-    tds.register_container<ContainerTH2Many<float>>("dz_vs_eta_B4", "dz vs eta (BPIX L4)",
-                                                      dz_vs_eta_B4, 20, -0.01, .01, 20, -3, 3,
-                                                      "dz",  "eta");
 
-}
-
-void run_analysis(const std::string& input_filename, const std::string& data_label, bool silent){
+/* void setup_dz_vs_eta(TrackingDataSet& tds){ */
+/*     auto& calc_dz_vs_eta = */
+/*     func<std::pair<vector<float>,vector<float>>(vector<Seed>, */
+/*                                                 vector<PixRecHit>, */
+/*                                                 vector<SimHit>, int)>("calc_dz_vs_eta", */
+/*             FUNC(([](const vector<Seed>& seeds, */
+/*                      const vector<PixRecHit> pixrec_hits, */
+/*                      const vector<SimHit> sim_hits, */
+/*                      int layer){ */
+/*                 vector<float> dz; */
+/*                 vector<float> eta; */
+/*                 for(const Seed &seed : seeds){ */
+/*                     int hitIdx, hitType; */
+/*                     for(auto tup : boost::combine(seed.hitIdx, seed.hitType)){ */
+/*                         boost::tie(hitIdx, hitType) = tup; */
+/*                         if(hitType == HIT_TYPE_PIXEL){ */
+/*                             const PixRecHit &hit = pixrec_hits[hitIdx]; */
+/*                             if(layer && hit.lay != layer) continue; */
+/*                             for(const int& simHitIdx : hit.simHitIdx){ */
+/*                                 dz.push_back(hit.z-sim_hits[simHitIdx].z); */
+/*                                 eta.push_back(seed.eta); */
+/*                             } */
+/*                         } */
+/*                     } */
+/*                 } */
+/*                 return std::make_pair(eta,dz); */
+/*             }))); */
+
+/*     auto dz_vs_eta_B1 = */
+/*     fv::apply(calc_dz_vs_eta, fv::tuple(seeds, pixrec_hits, sim_hits, constant<int>("1", 1)), "dz_vs_eta_B1"); */
+/*     auto dz_vs_eta_B2 = */
+/*     fv::apply(calc_dz_vs_eta, fv::tuple(seeds, pixrec_hits, sim_hits, constant<int>("2", 2)), "dz_vs_eta_B2"); */
+/*     auto dz_vs_eta_B3 = */
+/*     fv::apply(calc_dz_vs_eta, fv::tuple(seeds, pixrec_hits, sim_hits, constant<int>("3", 3)), "dz_vs_eta_B3"); */
+/*     auto dz_vs_eta_B4 = */
+/*     fv::apply(calc_dz_vs_eta, fv::tuple(seeds, pixrec_hits, sim_hits, constant<int>("4", 4)), "dz_vs_eta_B4"); */
+
+
+/*     tds.register_container<ContainerTH2Many<float>>("dz_vs_eta_B1", "dz vs eta (BPIX L1)", */
+/*                                                       dz_vs_eta_B1, 20, -3, 3, 20, -.01, .01, */
+/*                                                       "eta",  "dz"); */
+/*     tds.register_container<ContainerTH2Many<float>>("dz_vs_eta_B2", "dz vs eta (BPIX L2)", */
+/*                                                       dz_vs_eta_B2, 20, -3, 3, 20, -.01, .01, */
+/*                                                       "eta",  "dz"); */
+/*     tds.register_container<ContainerTH2Many<float>>("dz_vs_eta_B3", "dz vs eta (BPIX L3)", */
+/*                                                       dz_vs_eta_B3, 20, -3, 3, 20, -.01, .01, */
+/*                                                       "eta",  "dz"); */
+/*     tds.register_container<ContainerTH2Many<float>>("dz_vs_eta_B4", "dz vs eta (BPIX L4)", */
+/*                                                       dz_vs_eta_B4, 20, -3, 3, 20, -.01, .01, */
+/*                                                       "eta",  "dz"); */
+
+/* } */
+
+void run_analysis(const vector<fv::util::DataFileDescriptor>& dfds, const string output_filename, bool silent){
     gSystem->Load("libfilval.so");
     auto replace_suffix = [](const std::string& input, const std::string& new_suffix){
         return input.substr(0, input.find_last_of(".")) + new_suffix;
     };
-    string log_filename = replace_suffix(input_filename, "_result.log");
+
+    string log_filename = replace_suffix(output_filename, ".log");
     fv::util::Log::init_logger(log_filename, fv::util::LogPriority::kLogDebug);
 
-    string output_filename = replace_suffix(input_filename, "_result.root");
-    TrackingDataSet tds(output_filename, input_filename, data_label, "trackingNtuple/tree");
+    TrackingDataSet tds(output_filename, dfds, "trackingNtuple/tree");
+    /* tds.set_max_events(10); */
     register_objects(tds);
 
-    setup_dz_vs_eta(tds);
-    setup_dphi_vs_eta(tds);
+    /* setup_dz_vs_eta(tds); */
+    /* setup_dphi_vs_eta(tds); */
+    setup_matched_tracks(tds);
 
 
     tds.process(silent);
@@ -172,12 +220,18 @@ void run_analysis(const std::string& input_filename, const std::string& data_lab
 int main(int argc, char * argv[])
 {
     fv::util::ArgParser args(argc, argv);
-    if(!args.cmdOptionExists("-l") || !args.cmdOptionExists("-f")) {
-        cout << "Usage: ./main (-s) -l DATA_LABEL -f treefile.root" << endl;
-        return -1;
-    }
     bool silent = args.cmdOptionExists("-s");
-    string input_filename = args.getCmdOption("-f");
-    string data_label = args.getCmdOption("-l");
-    run_analysis(input_filename, data_label, silent);
+    string output_filename = args.cmdOptionExists("-o") ? args.getCmdOption("-o") : "output.root";
+    if(args.cmdOptionExists("-F")){
+        auto file_list = fv::util::read_input_list(args.getCmdOption("-F"));
+        run_analysis(file_list, output_filename, silent);
+    }else if(args.cmdOptionExists("-f")) {
+        string input_filename = args.getCmdOption("-f");
+        string data_label = args.cmdOptionExists("-l") ? args.getCmdOption("-l") : "";
+        fv::util::DataFileDescriptor dfd(input_filename, data_label);
+        run_analysis(std::vector<fv::util::DataFileDescriptor>({dfd}), output_filename, silent);
+    }else {
+        cout << "Usage: ./tracking_validation (-s) {-o output_filename} -F datafiles.txt" << endl;
+        cout << "       ./tracking_validation (-s) {-l DATA_LABEL} {-o output_filename} -f treefile.root" << endl;
+    }
 }

filval

@@ -1 +1 @@
-Subproject commit 7c5a06c426dbfff8a46a305d5b084b7f59053da9
+Subproject commit 7f7e4233bec36e5ba9794ea2ada4dc74bef8eb38

util/DetId.hpp

@@ -0,0 +1,66 @@
+#ifndef UTIL_DETID_H
+#define UTIL_DETID_H
+
+
+#include <stdint.h>
+/** \class DetId
+
+Parent class for all detector ids in CMS.  The DetId is a 32-bit
+unsigned integer.  The four most significant bits ([31:28]) identify
+the large-scale detector (e.g. Tracker or Ecal) while the next three
+bits ([27:25]) identify a part of the detector (such as HcalBarrel
+(HB) for Hcal).
+
+// Taken from cmssw_9_0_X:DataFormats/DetId/interface/DetId.h
+*/
+class DetId {
+public:
+  static const int kDetOffset          = 28;
+  static const int kSubdetOffset       = 25;
+
+
+  enum Detector { Tracker=1,Muon=2,Ecal=3,Hcal=4,Calo=5,Forward=6,VeryForward=7 };
+  /// Create an empty or null id (also for persistence)
+  DetId()  : id_(0) { }
+  /// Create an id from a raw number
+  DetId(uint32_t id) : id_(id) { }
+  /// Create an id, filling the detector and subdetector fields as specified
+  DetId(Detector det, int subdet)  {
+    id_=((det&0xF)<<28)|((subdet&0x7)<<25);
+  }
+
+  /// get the detector field from this detid
+  Detector det() const { return Detector((id_>>kDetOffset)&0xF); }
+  /// get the contents of the subdetector field (not cast into any detector's numbering enum)
+  int subdetId() const { return ((id_>>kSubdetOffset)&0x7); }
+
+  uint32_t operator()() const { return id_; }
+  operator uint32_t() const { return id_; }
+
+  /// get the raw id
+  uint32_t rawId() const { return id_; }
+  /// is this a null id ?
+  bool null() const { return id_==0; }
+
+  /// equality
+  bool operator==(DetId id) const { return id_==id.id_; }
+  /// inequality
+  bool operator!=(DetId id) const { return id_!=id.id_; }
+  /// comparison
+  bool operator<(DetId id) const { return id_<id.id_; }
+
+protected:
+  uint32_t id_;
+};
+
+/// equality
+inline bool operator==(uint32_t i, DetId id)  { return i==id(); }
+inline bool operator==(DetId id, uint32_t i)  { return i==id(); }
+/// inequality
+inline bool operator!=(uint32_t i, DetId id)  { return i!=id(); }
+inline bool operator!=(DetId id, uint32_t i) { return i!=id(); }
+/// comparison
+inline bool operator<(uint32_t i, DetId id) { return i<id(); }
+inline bool operator<(DetId id, uint32_t i) { return id()<i; }
+
+#endif

util/PXBDetId.hpp

@@ -0,0 +1,57 @@
+#ifndef DataFormats_SiStripDetId_PXBDetId_H
+#define DataFormats_SiStripDetId_PXBDetId_H
+
+#include <ostream>
+#include "util/DetId.hpp"
+#include "util/PixelSubdetector.hpp"
+
+/**
+ *  Det identifier class for the PixelBarrel
+ */
+
+class PXBDetId;
+
+std::ostream& operator<<(std::ostream& os,const PXBDetId& id);
+
+class PXBDetId : public DetId {
+ public:
+  /** Constructor of a null id */
+  PXBDetId();
+  /** Constructor from a raw value */
+  PXBDetId(uint32_t rawid);
+  /**Construct from generic DetId */
+  PXBDetId(const DetId& id);
+
+  PXBDetId(uint32_t layer, uint32_t ladder, uint32_t module)
+    :DetId(DetId::Tracker,PixelSubdetector::PixelBarrel){
+      id_ |= (layer& layerMask_)   << layerStartBit_     |
+             (ladder& ladderMask_) << ladderStartBit_  |
+             (module& moduleMask_) << moduleStartBit_;
+  }
+
+
+  /// layer id
+  unsigned int layer() const{
+    return int((id_>>layerStartBit_) & layerMask_);}
+
+  /// ladder  id
+  unsigned int ladder() const
+    { return ((id_>>ladderStartBit_) & ladderMask_) ;}
+
+  /// det id
+  unsigned int module() const
+    { return ((id_>>moduleStartBit_)& moduleMask_) ;}
+
+ private:
+  /// two bits would be enough, but  we could use the number "0" as a wildcard
+  static const unsigned int layerStartBit_=   16;
+  static const unsigned int ladderStartBit_=   8;
+  static const unsigned int moduleStartBit_=      2;
+  /// two bits would be enough, but  we could use the number "0" as a wildcard
+  static const unsigned int layerMask_=       0xF;
+  static const unsigned int ladderMask_=      0xFF;
+  static const unsigned int moduleMask_=         0x3F;
+};
+
+
+#endif

util/PXFDetId.hpp

@@ -0,0 +1,75 @@
+#ifndef DataFormats_SiStripDetId_PXFDetId_H
+#define DataFormats_SiStripDetId_PXFDetId_H
+
+#include <ostream>
+#include "DataFormats/DetId/interface/DetId.h"
+#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
+
+/** 
+ *  Det identifier class for the PixelEndcap
+ */
+class PXFDetId;
+
+std::ostream& operator<<(std::ostream& os,const PXFDetId& id);
+
+class PXFDetId : public DetId {
+ public:
+  /** Constructor of a null id */
+  PXFDetId();
+  /** Constructor from a raw value */
+  PXFDetId(uint32_t rawid);
+  /**Construct from generic DetId */
+  PXFDetId(const DetId& id); 
+  
+  PXFDetId(uint32_t side,
+	   uint32_t disk,
+	   uint32_t blade,
+	   uint32_t panel,
+	   uint32_t module) : DetId(DetId::Tracker,PixelSubdetector::PixelEndcap){
+    id_ |= (side& sideMask_)  << sideStartBit_   |
+      (disk& diskMask_)        << diskStartBit_      |
+      (blade& bladeMask_)      << bladeStartBit_     |
+      (panel& panelMask_)      << panelStartBit_     |
+      (module& moduleMask_)    << moduleStartBit_  ;
+  }
+  
+  
+  /// positive or negative id
+  unsigned int side() const{
+    return int((id_>>sideStartBit_) & sideMask_);
+  }
+  
+  /// disk id
+  unsigned int disk() const{
+    return int((id_>>diskStartBit_) & diskMask_);
+  }
+  
+  /// blade id
+  unsigned int blade() const
+    { return ((id_>>bladeStartBit_) & bladeMask_) ;}
+  
+ /// panel id
+  unsigned int panel() const
+    { return ((id_>>panelStartBit_) & panelMask_) ;}
+
+  /// det id
+  unsigned int module() const
+    { return ((id_>>moduleStartBit_) & moduleMask_) ;}
+  
+ private:
+  /// two bits would be enough, but  we could use the number "0" as a wildcard
+  static const unsigned int sideStartBit_=   23;
+  static const unsigned int diskStartBit_=   16;
+  static const unsigned int bladeStartBit_=  10;
+  static const unsigned int panelStartBit_=  8;
+  static const unsigned int moduleStartBit_= 2;
+  /// two bits would be enough, but  we could use the number "0" as a wildcard
+  
+  static const unsigned int sideMask_=     0x3;
+  static const unsigned int diskMask_=     0xF;
+  static const unsigned int bladeMask_=    0x3F;
+  static const unsigned int panelMask_=    0x3;
+  static const unsigned int moduleMask_=   0x3F;
+};
+
+#endif

util/PixelSubdetector.h

@@ -0,0 +1,15 @@
+#ifndef DataFormats_SiPixelDetId_PixelSubdetector_H
+#define DataFormats_SiPixelDetId_PixelSubdetector_H
+
+/** 
+ *  Enumeration for Pixel Tracker Subdetectors
+ *
+ */
+
+class PixelSubdetector { 
+ public:
+  enum SubDetector {PixelBarrel=1,PixelEndcap=2};
+
+};
+
+#endif