ResidualRefitting.cc

Go to the documentation of this file.

#include <iostream>
#include "ResidualRefitting.h"
#include <iomanip>

//framework includes
#include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
#include "DataFormats/MuonDetId/interface/DTWireId.h"
#include "DataFormats/MuonDetId/interface/CSCDetId.h"
#include "DataFormats/MuonDetId/interface/RPCDetId.h"
#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
#include "Geometry/Records/interface/TrackerTopologyRcd.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Utilities/interface/EDMException.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
#include "Geometry/CommonDetUnit/interface/GeomDet.h"
#include "Geometry/CommonDetUnit/interface/GeomDetEnumerators.h"

#include "RecoMuon/DetLayers/interface/MuonDetLayerGeometry.h"
#include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBuilder.h"
#include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"

#include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"

#include "TrackingTools/DetLayers/interface/DetLayer.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"

#include "TrackingTools/Records/interface/TrackingComponentsRecord.h"



ResidualRefitting::ResidualRefitting( const edm::ParameterSet & cfg ) :
  outputFileName_       ( cfg.getUntrackedParameter<std::string>("histoutputFile") ),
  PropagatorSource_     ( cfg.getParameter<std::string>("propagator")), 
  muons_                ( cfg.getParameter<edm::InputTag>( "muons"      ) ),
  muonsRemake_          ( cfg.getParameter<edm::InputTag>("muonsRemake" ) ),            //This Feels Misalignment
  muonsNoStation1_      ( cfg.getParameter<edm::InputTag>("muonsNoStation1") ),
  muonsNoStation2_      ( cfg.getParameter<edm::InputTag>("muonsNoStation2") ),
  muonsNoStation3_      ( cfg.getParameter<edm::InputTag>("muonsNoStation3") ),
  muonsNoStation4_      ( cfg.getParameter<edm::InputTag>("muonsNoStation4") ),

/*
  muonsNoPXBLayer1_     ( cfg.getParameter<edm::InputTag>("muonsNoPXBLayer1"    ) ),
  muonsNoPXBLayer2_     ( cfg.getParameter<edm::InputTag>("muonsNoPXBLayer1"    ) ),
  muonsNoPXBLayer3_     ( cfg.getParameter<edm::InputTag>("muonsNoPXBLayer1"    ) ),

  muonsNoTIBLayer1_         ( cfg.getParameter<edm::InputTag>("muonsNoTIBLayer1"    ) ),
  muonsNoTIBLayer2_         ( cfg.getParameter<edm::InputTag>("muonsNoTIBLayer2"    ) ),
  muonsNoTIBLayer3_         ( cfg.getParameter<edm::InputTag>("muonsNoTIBLayer3"    ) ),
  muonsNoTIBLayer4_         ( cfg.getParameter<edm::InputTag>("muonsNoTIBLayer4"    ) ),

  muonsNoTOBLayer1_         ( cfg.getParameter<edm::InputTag>("muonsNoTOBLayer1"    ) ),
  muonsNoTOBLayer2_         ( cfg.getParameter<edm::InputTag>("muonsNoTOBLayer2"    ) ),
  muonsNoTOBLayer3_         ( cfg.getParameter<edm::InputTag>("muonsNoTOBLayer3"    ) ),
  muonsNoTOBLayer4_         ( cfg.getParameter<edm::InputTag>("muonsNoTOBLayer4"    ) ),
  muonsNoTOBLayer5_         ( cfg.getParameter<edm::InputTag>("muonsNoTOBLayer5"    ) ),
  muonsNoTOBLayer6_         ( cfg.getParameter<edm::InputTag>("muonsNoTOBLayer6"    ) ),*/
  debug_                ( cfg.getUntrackedParameter<bool>("doDebug" ) ),
  outputFile_(0),
  outputTree_(0),
  outputBranch_(0),
  theField(0)
{
        eventInfo_.evtNum_ = 0;
        eventInfo_.evtNum_ = 0;

// service parameters
    edm::ParameterSet serviceParameters = cfg.getParameter<edm::ParameterSet>("ServiceParameters");
  
// the services
    theService = new MuonServiceProxy(serviceParameters);
    
}  //The constructor

void ResidualRefitting::analyze(const edm::Event& event, const edm::EventSetup& eventSetup) {
    
    if(debug_) printf("STARTING EVENT\n");

    eventInfo_.evtNum_ = (int)event.id().run();
    eventInfo_.runNum_ = (int)event.id().event();

    // Generator Collection

// The original muon collection that is sitting in memory
    edm::Handle<reco::MuonCollection> muons;

    edm::Handle<reco::TrackCollection> muonTracks;
    edm::Handle<reco::TrackCollection> muonsNoSt1;
    edm::Handle<reco::TrackCollection> muonsNoSt2;
    edm::Handle<reco::TrackCollection> muonsNoSt3;
    edm::Handle<reco::TrackCollection> muonsNoSt4;


    event.getByLabel(muons_             , muons ); //set label to muons
    event.getByLabel(muonsRemake_       , muonTracks);  
    event.getByLabel(muonsNoStation1_   , muonsNoSt1);
    event.getByLabel(muonsNoStation2_   , muonsNoSt2);
    event.getByLabel(muonsNoStation3_   , muonsNoSt3);
    event.getByLabel(muonsNoStation4_   , muonsNoSt4);
    
    
/*
//  std::cout<<"Muon Collection No PXB "<<std::endl;
//Tracker Barrel Pixel Refits
    edm::Handle<MuonCollection> muonsNoPXBLayer1Coll;
    event.getByLabel(muonsNoPXBLayer1_, muonsNoPXBLayer1Coll);
    edm::Handle<MuonCollection> muonsNoPXBLayer2Coll;
    event.getByLabel(muonsNoPXBLayer2_, muonsNoPXBLayer2Coll);
    edm::Handle<MuonCollection> muonsNoPXBLayer3Coll;
    event.getByLabel(muonsNoPXBLayer3_, muonsNoPXBLayer3Coll);
//  std::cout<<"Muon Collection No TIB "<<std::endl;
// Tracker Inner Barrel Refits
    edm::Handle<MuonCollection> muonsNoTIBLayer1Coll;
    event.getByLabel(muonsNoTIBLayer1_, muonsNoTIBLayer1Coll);
    edm::Handle<MuonCollection> muonsNoTIBLayer2Coll;
    event.getByLabel(muonsNoTIBLayer2_, muonsNoTIBLayer2Coll);
    edm::Handle<MuonCollection> muonsNoTIBLayer3Coll;
    event.getByLabel(muonsNoTIBLayer3_, muonsNoTIBLayer3Coll);
    edm::Handle<MuonCollection> muonsNoTIBLayer4Coll;
    event.getByLabel(muonsNoTIBLayer4_, muonsNoTIBLayer4Coll);

//  std::cout<<"Muon Collection No TOB "<<std::endl;

//Tracker outer barrel refits
    edm::Handle<MuonCollection> muonsNoTOBLayer1Coll;
    event.getByLabel(muonsNoTOBLayer1_, muonsNoTOBLayer1Coll);
    edm::Handle<MuonCollection> muonsNoTOBLayer2Coll;
    event.getByLabel(muonsNoTOBLayer2_, muonsNoTOBLayer2Coll);
    edm::Handle<MuonCollection> muonsNoTOBLayer3Coll;
    event.getByLabel(muonsNoTOBLayer3_, muonsNoTOBLayer3Coll);
    edm::Handle<MuonCollection> muonsNoTOBLayer4Coll;
    event.getByLabel(muonsNoTOBLayer4_, muonsNoTOBLayer4Coll);
    edm::Handle<MuonCollection> muonsNoTOBLayer5Coll;
    event.getByLabel(muonsNoTOBLayer5_, muonsNoTOBLayer5Coll);
    edm::Handle<MuonCollection> muonsNoTOBLayer6Coll;
    event.getByLabel(muonsNoTOBLayer6_, muonsNoTOBLayer6Coll);
*/
    //magnetic field information    
    edm::ESHandle<MagneticField> field;
    edm::ESHandle<GlobalTrackingGeometry> globalTrackingGeometry;
    eventSetup.get<IdealMagneticFieldRecord>().get(field);
    eventSetup.get<GlobalTrackingGeometryRecord>().get(globalTrackingGeometry);
    eventSetup.get<TrackingComponentsRecord>().get( PropagatorSource_, thePropagator );
    theField = &*field;
    
    
    theService->update(eventSetup);

//Zero storage
    zero_storage();

  
//Do the Gmr Muons from the unModified Collection

/*
    int iGmr = 0;
    if ( (muons->end() - muons->begin()) > 0) printf("Data Dump:: Original GMR Muons\n");
    for ( MuonCollection::const_iterator muon = muons->begin(); muon!=muons->end(); muon++, iGmr++) {
        if ( iGmr >= ResidualRefitting::N_MAX_STORED) break; // error checking
        if (!debug
        
        dumpTrackRef(muon->combinedMuon(), "cmb"); 
        dumpTrackRef(muon->standAloneMuon(), "sam");
        dumpTrackRef(muon->track(), "trk");
        

    }
    storageGmrOld_.n_ = iGmr;
    storageSamNew_.n_ = iGmr;
*/

//Refitted muons
    if (debug_) printf("Data Dump:: Rebuilt GMR Muon Track With TeV refitter default\n");
    int iGmrRemake = 0;
    for ( reco::TrackCollection::const_iterator muon = muonTracks->begin(); muon!=muonTracks->end(); muon++, iGmrRemake++) {
        if ( iGmrRemake >= ResidualRefitting::N_MAX_STORED) break; // error checking
            // from TrackInfoProducer/test/TrackInfoAnalyzerExample.cc
            reco::TrackRef trackref=reco::TrackRef(muonTracks,iGmrRemake);

            if (debug_) dumpTrackRef(trackref, "gmr");
            muonInfo(storageGmrNew_,trackref,iGmrRemake);

    }
    storageGmrNew_.n_ = iGmrRemake;
        


    if (debug_) printf("muons Remake");
    if (debug_) printf("-----------------------------------------\n");
    CollectTrackHits(muonTracks, storageTrackExtrapRec_, eventSetup);


    if (true) {
        printf("muons No Station 1");
        printf("-----------------------------------------\n");
    }
    NewTrackMeasurements(muonTracks, muonsNoSt1, storageTrackExtrapRecNoSt1_);

    if (true) {
        printf("muons No Station 2");
        printf("-----------------------------------------\n");
    }
    NewTrackMeasurements(muonTracks, muonsNoSt2, storageTrackExtrapRecNoSt2_);

    if (true) {
        printf("muons No Station 3");
        printf("-----------------------------------------\n");
    }
    NewTrackMeasurements(muonTracks, muonsNoSt3, storageTrackExtrapRecNoSt3_);

    if (true) {
        printf("muons No Station 4");
        printf("-----------------------------------------\n");
    }
    NewTrackMeasurements(muonTracks, muonsNoSt4, storageTrackExtrapRecNoSt4_);


//  dumpMuonRecHits(storageRecMuon_); 
        
/****************************************************************************************************************************************/

  
/*
 *  extrapolates track to a cylinder.
 *  commented for cosmic runs with no tracker in reco muons!!
 *
*/


    int iGmrCyl = 0;
    for (reco::MuonCollection::const_iterator muon = muons->begin(); muon != muons->end(); muon++, iGmrCyl++) {

        dumpTrackRef(muon->combinedMuon(), "cmb"); 
        dumpTrackRef(muon->standAloneMuon(), "sam");
        dumpTrackRef(muon->track(), "trk");

        cylExtrapTrkSam(iGmrCyl, muon->standAloneMuon() , samExtrap120_     , 120.);
        cylExtrapTrkSam(iGmrCyl, muon->track()          , trackExtrap120_   , 120.);

    }
    samExtrap120_.n_     = iGmrCyl;
    trackExtrap120_.n_   = iGmrCyl; 


    if (iGmrRemake > 0 || iGmrCyl > 0) {
        outputTree_ -> Fill();
        std::cout << "FILLING NTUPLE!" << std::endl;
        std::cout << "Entries Recorded: " << outputTree_ -> GetEntries() << " Branch :: " << outputBranch_ -> GetEntries() <<  std::endl<<std::endl;
    } else std::cout<<"no tracks -- no fill!\n"<<std::endl<<std::endl;

//  /*************************************************************************************************************/
//  //END OF ntuple dumper
//  //END OF ntuple dumper
//  /***********************************************************************************************************/
}
//end Analyze() main function

//------------------------------------------------------------------------------
//  
// Destructor
// 
ResidualRefitting::~ResidualRefitting() {
  delete outputFile_;
  delete theService;
}
//
// Track Collection Analysis
//
void ResidualRefitting::CollectTrackHits(edm::Handle<reco::TrackCollection> trackColl,
                                         ResidualRefitting::storage_trackExtrap& trackExtrap,
                                         const edm::EventSetup& eventSetup) {

    //Retrieve tracker topology from geometry
    edm::ESHandle<TrackerTopology> tTopoHandle;
    eventSetup.get<TrackerTopologyRcd>().get(tTopoHandle);
    const TrackerTopology* const tTopo = tTopoHandle.product();

    int iMuonHit = 0;
    int iTrackHit= 0;
    int numTracks= 0;   

    for ( reco::TrackCollection::const_iterator muon = trackColl->begin(); muon!=trackColl->end(); muon++) {

        int iTrack = muon - trackColl->begin();
        reco::TrackRef trackref=reco::TrackRef(trackColl,iTrack);
        FreeTrajectoryState recoStart = ResidualRefitting::freeTrajStateMuon(trackref);

        if (debug_) dumpTrackRef(trackref, "CollectTrackHits Track");
    
        for (trackingRecHit_iterator rec = muon->recHitsBegin(); rec != muon->recHitsEnd(); rec++) {

            int iRec = rec - muon->recHitsBegin();  
            DetId detid = (*rec)->geographicalId(); 

            if (detid.det() != DetId::Muon && detid.det() != DetId::Tracker) {
                if (debug_) printf("Rec Hit not from muon system or tracker... continuing...\n");
                continue;
            }
//          numTracks++;
// Get Local and Global Position of Hits    

            LocalPoint lp = (*rec)->localPosition();
            float lpX   = lp.x();
            float lpY   = lp.y();
            float lpZ   = lp.z();

            MuonTransientTrackingRecHit::MuonRecHitPointer mrhp =   
                MuonTransientTrackingRecHit::specificBuild(theService->trackingGeometry()->idToDet((**rec).geographicalId())
                ,&(**rec)); 

            GlobalPoint gp = mrhp->globalPosition();
            float gpRecX    = gp.x();
            float gpRecY    = gp.y();
            float gpRecZ    = gp.z();
            float gpRecEta  = gp.eta();
            float gpRecPhi  = gp.phi();

            if (detid.det() == DetId::Muon) {
            
                int systemMuon  = detid.subdetId(); // 1 DT; 2 CSC; 3 RPC
                int endcap      = -999;
                int station     = -999;
                int ring        = -999;
                int chamber     = -999;
                int layer       = -999;
                int superLayer  = -999;
                int wheel       = -999;
                int sector      = -999;
                if ( systemMuon == MuonSubdetId::CSC) {
                    CSCDetId id(detid.rawId());
                    endcap      = id.endcap();
                    station     = id.station();
                    ring        = id.ring();
                    chamber     = id.chamber();
                    layer       = id.layer();
                    if (debug_)printf("CSC\t[endcap][station][ringN][chamber][layer]:[%d][%d][%d][%d][%d]\t",
                         endcap, station, ring, chamber, layer);

                }
                else if ( systemMuon == MuonSubdetId::DT ) {
                    DTWireId id(detid.rawId());
                    station     = id.station();
                    layer       = id.layer();
                    superLayer  = id.superLayer();
                    wheel       = id.wheel();
                    sector      = id.sector();
                    if (debug_) printf("DT \t[station][layer][superlayer]:[%d][%d][%d]\n", station,layer,superLayer);
                
                }
                else if ( systemMuon == MuonSubdetId::RPC) {
                    RPCDetId id(detid.rawId());
                    station     = id.station();
                    if (debug_) printf("RPC\t[station]:[%d]\n", station);
                }

            
                storageRecMuon_.muonLink_   [iMuonHit]  = iTrack;           
                storageRecMuon_.system_     [iMuonHit]  = systemMuon;
                storageRecMuon_.endcap_     [iMuonHit]  = endcap;
                storageRecMuon_.station_    [iMuonHit]  = station;
                storageRecMuon_.ring_       [iMuonHit]  = ring;
                storageRecMuon_.chamber_    [iMuonHit]  = chamber;
                storageRecMuon_.layer_      [iMuonHit]  = layer;
                storageRecMuon_.superLayer_ [iMuonHit]  = superLayer;
                storageRecMuon_.wheel_      [iMuonHit]  = wheel;
                storageRecMuon_.sector_     [iMuonHit]  = sector;
        
                storageRecMuon_.gpX_        [iMuonHit]  = gpRecX;
                storageRecMuon_.gpY_        [iMuonHit]  = gpRecY;
                storageRecMuon_.gpZ_        [iMuonHit]  = gpRecZ;
                storageRecMuon_.gpEta_      [iMuonHit]  = gpRecEta;
                storageRecMuon_.gpPhi_      [iMuonHit]  = gpRecPhi;
                storageRecMuon_.lpX_        [iMuonHit]  = lpX;
                storageRecMuon_.lpY_        [iMuonHit]  = lpY;
                storageRecMuon_.lpZ_        [iMuonHit]  = lpZ;
                iMuonHit++;
    
            }
            else if (detid.det() == DetId::Tracker) {
                
                if (debug_) printf("Tracker\n");

                StoreTrackerRecHits(detid, tTopo, iTrack, iTrackHit);

                storageTrackHit_.gpX_       [iTrackHit] = gpRecX;
                storageTrackHit_.gpY_       [iTrackHit] = gpRecY;
                storageTrackHit_.gpZ_       [iTrackHit] = gpRecZ;
                storageTrackHit_.gpEta_     [iTrackHit] = gpRecEta;
                storageTrackHit_.gpPhi_     [iTrackHit] = gpRecPhi;
                storageTrackHit_.lpX_       [iTrackHit] = lpX;
                storageTrackHit_.lpY_       [iTrackHit] = lpY;
                storageTrackHit_.lpZ_       [iTrackHit] = lpZ;
                iTrackHit++;
             }      
            else printf("THIS CAN NOT HAPPEN\n");       
            
            trkExtrap(detid, numTracks, iTrack, iRec, recoStart, lp, trackExtrap);
            numTracks++;

            if (debug_) printf("\tLocal Positon:  \tx = %2.2f\ty = %2.2f\tz = %2.2f\n",lpX, lpY, lpZ);
            if (debug_) printf("\tGlobal Position: \tx = %6.2f\ty = %6.2f\tz = %6.2f\teta = %4.2f\tphi = %3.2f\n",
                            gpRecX,gpRecY,gpRecZ,gpRecEta,gpRecPhi);
                

        }

    }

    storageRecMuon_ .n_     = iMuonHit; 
    storageTrackHit_.n_     = iTrackHit;
    trackExtrap     .n_     = numTracks;

}
//
// Deal with Re-Fitted Track with some station omitted. 
//
// This should take the new track, match it to its track before refitting with the hits dumped, and extrapolate out to the
//  rec hits that were removed from the fit.  
//
//

void ResidualRefitting::NewTrackMeasurements(edm::Handle<reco::TrackCollection> trackCollOrig,
     edm::Handle<reco::TrackCollection> trackColl,  ResidualRefitting::storage_trackExtrap& trackExtrap) {

    int numTracks   = 0;    
    int recCounter  = 0;

    for ( reco::TrackCollection::const_iterator muon = trackColl->begin(); muon!=trackColl->end(); muon++) {

        int iTrack = muon - trackColl->begin(); 
    
        reco::TrackRef trackref=reco::TrackRef(trackColl,iTrack);
        FreeTrajectoryState recoStart = ResidualRefitting::freeTrajStateMuon(trackref);

        int iTrackLink =  MatchTrackWithRecHits(muon,trackCollOrig);
        reco::TrackRef ref = reco::TrackRef(trackCollOrig, iTrackLink);

        for (trackingRecHit_iterator rec1 = ref->recHitsBegin(); rec1!= ref->recHitsEnd(); rec1++, recCounter++) {

            //int iRec = rec1 - ref->recHitsBegin();
            
            bool unbiasedRec = true;

            for ( trackingRecHit_iterator rec2 = muon->recHitsBegin(); rec2!=muon->recHitsEnd();rec2++) {
                        
                if (IsSameHit(rec1,rec2)) {
                    unbiasedRec = false;
                    break;
                }
            }
            if (!unbiasedRec) continue;
                
            DetId detid = (*rec1)->geographicalId(); 

            MuonTransientTrackingRecHit::MuonRecHitPointer mrhp =   
                MuonTransientTrackingRecHit::specificBuild(theService->trackingGeometry()->idToDet((**rec1).geographicalId())
                ,&(**rec1)); 
        
            trkExtrap(detid, numTracks, iTrackLink, recCounter, recoStart, (*rec1)->localPosition(), trackExtrap);
            numTracks++;    
        
        }   

    }

    trackExtrap.n_ = numTracks;

}
//
// Find the original track that corresponds to the re-fitted track
//
int ResidualRefitting::MatchTrackWithRecHits(reco::TrackCollection::const_iterator trackIt,
     edm::Handle<reco::TrackCollection> ref) {

    if (debug_) printf("Matching a re-fitted track to the original track.\n");
    
    int TrackMatch = -1;

    for (trackingRecHit_iterator rec = trackIt->recHitsBegin(); rec!=trackIt->recHitsEnd(); rec++) {
        
        bool foundMatch = false;
        for (reco::TrackCollection::const_iterator refIt = ref->begin(); refIt!=ref->end(); refIt++) {
        
            int iTrackMatch = refIt - ref->begin();
            if (foundMatch && TrackMatch !=iTrackMatch) break;
            for (trackingRecHit_iterator recRef = refIt->recHitsBegin(); recRef!=refIt->recHitsEnd(); recRef++) {

                if (!IsSameHit(rec,recRef)) continue;
                
                foundMatch = true;
                TrackMatch = iTrackMatch;
            //  printf("Rec hit match for original track %d\n", iTrackMatch);
        
            }
        }
        if (!foundMatch) {
            printf("SOMETHING WENT WRONG! Could not match Track with original track!");
            exit(1);
        }
        
    }
    if (debug_) printf("Rec hit match for original track %d\n", TrackMatch);

//  reco::TrackRef trackref=reco::TrackRef(ref,TrackMatch);
    return TrackMatch;
}
/*
//
// Match two tracks to see if one is a subset of the other
//

bool ResidualRefitting::TrackSubset(reco::TrackRef trackSub, reco::TrackRef trackTop) {

    
    bool matchAll = true;

    for (trackingRecHit_iterator recSub = trackSub->recHitsBegin(); recSub!=trackSub->recHitsEnd(); recSub++) {

        bool matchSub = false;


        for (trackingRecHit_iterator recTop = trackTop->recHitsBegin(); recTop!=trackTop->recHitsEnd(); recTop++) {
        
            if ( recSub == recTop ) matchSub = true;
            if (matchSub) break;

        }
        if (!matchSub) return false;

    }

    return matchAll;

}
*/

//
// Check to see if the rec hits are the same
//
bool ResidualRefitting::IsSameHit(trackingRecHit_iterator hit1, trackingRecHit_iterator hit2) {

    
    double lpx1 = (*hit1)->localPosition().x();
    double lpy1 = (*hit1)->localPosition().y();
    double lpz1 = (*hit1)->localPosition().z();

    double lpx2 = (*hit2)->localPosition().x();
    double lpy2 = (*hit2)->localPosition().y();
    double lpz2 = (*hit2)->localPosition().z();
    if ( fabs( lpx1 - lpx2) > 1e-3) return false;
//  printf("Match lpx...\n");
    if ( fabs( lpy1 - lpy2) > 1e-3) return false;
//  printf("Match lpy...\n");
    if ( fabs( lpz1 - lpz2) > 1e-3) return false;
//  printf("Match lpz...\n");

    return true;

}

//
// Store Tracker Rec Hits
//
void ResidualRefitting::StoreTrackerRecHits(DetId detid, const TrackerTopology* tTopo, int iTrack, int iRec) {


            int detector    = -1;
            int subdetector = -1;
            int blade       = -1;
            int disk        = -1;
            int ladder      = -1;
            int layer       = -1;
            int module      = -1;
            int panel       = -1;
            int ring        = -1;
            int side        = -1;
            int wheel       = -1;
            
//Detector Info

            detector = detid.det();
            subdetector = detid.subdetId();
            
            if (detector != DetId::Tracker) { 
                std::cout<<"OMFG NOT THE TRACKER\n"<<std::endl;
                return;
            }

            if (debug_) std::cout<<"Tracker:: ";
            if (subdetector == ResidualRefitting::PXB) {
                layer   = tTopo->pxbLayer(detid.rawId());
                ladder  = tTopo->pxbLadder(detid.rawId());
                module  = tTopo->pxbModule(detid.rawId());
                if (debug_) std::cout   <<  "PXB"
                                        <<  "\tlayer = "    << layer
                                        <<  "\tladder = "   << ladder
                                        <<  "\tmodule = "   << module;
                
            } 
            else if (subdetector == ResidualRefitting::PXF) {
                side    = tTopo->pxfSide(detid.rawId());
                disk    = tTopo->pxfDisk(detid.rawId());
                blade   = tTopo->pxfBlade(detid.rawId());
                panel   = tTopo->pxfPanel(detid.rawId());
                module  = tTopo->pxfModule(detid.rawId());
                if (debug_) std::cout   <<  "PXF"
                                        <<  "\tside = "     << side
                                        <<  "\tdisk = "     << disk
                                        <<  "\tblade = "    << blade
                                        <<  "\tpanel = "    << panel
                                        <<  "\tmodule = "   << module;
                            
            }
            else if (subdetector == ResidualRefitting::TIB) {
                layer   = tTopo->tibLayer(detid.rawId());
                module  = tTopo->tibModule(detid.rawId());
                if (debug_) std::cout   << "TIB"
                                        << "\tlayer = " << layer
                                        << "\tmodule = "<< module;
            }
            else if (subdetector == ResidualRefitting::TID) {
                side    = tTopo->tidSide(detid.rawId());
                wheel   = tTopo->tidWheel(detid.rawId());
                ring    = tTopo->tidRing(detid.rawId());
                if (debug_) std::cout   <<"TID"
                                        << "\tside = "  << side
                                        << "\twheel = " << wheel
                                        << "\tring = "  << ring;
            
            }
            else if (subdetector == ResidualRefitting::TOB) {
                layer   = tTopo->tobLayer(detid.rawId());
                module  = tTopo->tobModule(detid.rawId());
                if (debug_) std::cout   <<"TOB"
                                        <<"\tlayer = "  << layer
                                        <<"\tmodule = " << module;
            
            }
            else if (subdetector == ResidualRefitting::TEC) {
                ring    = tTopo->tecRing(detid.rawId());
                module  = tTopo->tecModule(detid.rawId());
                if (debug_) std::cout   <<"TEC"
                                        << "\tring = "  << ring
                                        << "\tmodule = "<< module;
            }


//Do Storage

            storageTrackHit_.muonLink_      [iRec] =iTrack      ;
            storageTrackHit_.detector_      [iRec] =detector    ;
            storageTrackHit_.subdetector_       [iRec] =subdetector ;
            storageTrackHit_.blade_         [iRec] =blade       ;
            storageTrackHit_.disk_          [iRec] =disk        ;
            storageTrackHit_.ladder_        [iRec] =ladder      ;
            storageTrackHit_.layer_         [iRec] =layer       ;
            storageTrackHit_.module_        [iRec] =module      ;
            storageTrackHit_.panel_         [iRec] =panel       ;
            storageTrackHit_.ring_          [iRec] =ring        ;
            storageTrackHit_.side_          [iRec] =side        ;
            storageTrackHit_.wheel_         [iRec] =wheel       ;
            
}

//
// Store Muon info on P, Pt, eta, phi
//
void ResidualRefitting::muonInfo(ResidualRefitting::storage_muon& storeMuon, reco::TrackRef muon, int val) {


    storeMuon.pt_ [val]         = muon->pt();
    storeMuon.p_  [val]         = muon->p();
    storeMuon.eta_[val]         = muon->eta();
    storeMuon.phi_[val]         = muon->phi();
    storeMuon.charge_[val]      = muon->charge();
    storeMuon.numRecHits_[val]  = muon->numberOfValidHits();
    storeMuon.chiSq_[val]       = muon->chi2();
    storeMuon.ndf_[val]         = muon->ndof();
    storeMuon.chiSqOvrNdf_[val] = muon->normalizedChi2();

}
// 
// Fill a track extrapolation 
// 
void ResidualRefitting::trkExtrap(const DetId& detid, int iTrk, int iTrkLink, int iRec,
                  const FreeTrajectoryState& freeTrajState,
                  const LocalPoint& recPoint,
                  storage_trackExtrap& storeTemp){

    bool dump_ = debug_;
    
    if (dump_) std::cout<< "In the trkExtrap function"<<std::endl;

    float gpExtrapX     = -99999;
    float gpExtrapY     = -99999;
    float gpExtrapZ     = -99999;
    float gpExtrapEta   = -99999;
    float gpExtrapPhi   = -99999;

    float lpX           = -99999;
    float lpY           = -99999;
    float lpZ           = -99999;

    //
    // Get the local positions for the recHits
    //

    float recLpX = recPoint.x();
    float recLpY = recPoint.y();
    float recLpZ = recPoint.z();

    float resX = -9999;
    float resY = -9999;
    float resZ = -9999;

    const GeomDet* gdet = theService->trackingGeometry()->idToDet(detid);
    
//  TrajectoryStateOnSurface surfTest =  prop.propagate(freeTrajState, gdet->surface());
    TrajectoryStateOnSurface surfTest =  thePropagator->propagate(freeTrajState, gdet->surface());
    
    if (surfTest.isValid()) {
    
        GlobalPoint globTest    = surfTest.globalPosition();        
        gpExtrapX               = globTest.x();
        gpExtrapY               = globTest.y();
        gpExtrapZ               = globTest.z();
        gpExtrapEta             = globTest.eta();
        gpExtrapPhi             = globTest.phi();
        LocalPoint loc          = surfTest.localPosition();
        if (detid.det() == DetId::Muon || detid.det() == DetId::Tracker) {
            lpX                     = loc.x();
            lpY                     = loc.y();
            lpZ                     = loc.z();

            resX = lpX - recLpX;
            resY = lpY - recLpY;
            resZ = lpZ - recLpZ;
            
        }

    }
    storeTemp.muonLink_ [iTrk] = iTrkLink   ;
    storeTemp.recLink_  [iTrk] = iRec       ;
    storeTemp.gpX_      [iTrk] = gpExtrapX  ;
    storeTemp.gpY_      [iTrk] = gpExtrapY  ;
    storeTemp.gpZ_      [iTrk] = gpExtrapZ  ;
    storeTemp.gpEta_    [iTrk] = gpExtrapEta;
    storeTemp.gpPhi_    [iTrk] = gpExtrapPhi;   
    storeTemp.lpX_      [iTrk] = lpX        ;
    storeTemp.lpY_      [iTrk] = lpY        ;
    storeTemp.lpZ_      [iTrk] = lpZ        ;
    storeTemp.resX_     [iTrk] = resX       ;
    storeTemp.resY_     [iTrk] = resY       ;
    storeTemp.resZ_     [iTrk] = resZ       ;
    
    printf("station: %d\tsector: %d\tresX storage: %4.2f\n", ReturnStation(detid), ReturnSector(detid), resX);

}
//
// Return the station 
//
int ResidualRefitting::ReturnStation(DetId detid) {
    
    int station     = -999;

    if (detid.det() == DetId::Muon) {
    
        int systemMuon  = detid.subdetId(); // 1 DT; 2 CSC; 3 RPC
        if ( systemMuon == MuonSubdetId::CSC) {
            CSCDetId id(detid.rawId());
            station     = id.station();

        }
        else if ( systemMuon == MuonSubdetId::DT ) {
            DTWireId id(detid.rawId());
            station     = id.station();
        
        }
        else if ( systemMuon == MuonSubdetId::RPC) {
            RPCDetId id(detid.rawId());
            station     = id.station();
        }

    }
    
    return station;     
}
//
// Return the sector 
//
int ResidualRefitting::ReturnSector(DetId detid) {
    
    int sector      = -999;

    if (detid.det() == DetId::Muon) {
    
        int systemMuon  = detid.subdetId(); // 1 DT; 2 CSC; 3 RPC
        if ( systemMuon == MuonSubdetId::DT ) {
            DTWireId id(detid.rawId());
            sector      = id.sector();
        
        }

    }
    
    return sector;      
}

// 
// Store the SAM and Track position info at a particular rho
// 
void ResidualRefitting::cylExtrapTrkSam(int recNum, reco::TrackRef track, ResidualRefitting::storage_trackExtrap& storage, double rho) {

    Cylinder::PositionType pos(0,0,0);
    Cylinder::RotationType rot;

    Cylinder::CylinderPointer myCylinder = Cylinder::build(pos, rot, rho);
//  SteppingHelixPropagator inwardProp  ( theField, oppositeToMomentum );
//  SteppingHelixPropagator outwardProp ( theField, alongMomentum );
    FreeTrajectoryState recoStart = freeTrajStateMuon(track);
//  TrajectoryStateOnSurface recoProp = outwardProp.propagate(recoStart, *myCylinder);
    TrajectoryStateOnSurface recoProp = thePropagator->propagate(recoStart, *myCylinder);

    double xVal     = -9999;
    double yVal     = -9999;
    double zVal     = -9999;
    double phiVal   = -9999;
    double etaVal   = -9999;

    if(recoProp.isValid()) {
        GlobalPoint recoPoint = recoProp.globalPosition();
        xVal = recoPoint.x();
        yVal = recoPoint.y();
        zVal = recoPoint.z();
        phiVal = recoPoint.phi();
        etaVal = recoPoint.eta();       
    }
    storage.muonLink_[recNum]   = recNum;
    storage.gpX_    [recNum]    = xVal;
    storage.gpY_    [recNum]    = yVal;
    storage.gpZ_    [recNum]    = zVal;
    storage.gpEta_  [recNum]    = etaVal;
    storage.gpPhi_  [recNum]    = phiVal;

    float rhoVal = sqrt( xVal*xVal + yVal*yVal);

    printf("Cylinder: rho = %4.2f\tphi = %4.2f\teta = %4.2f\n", rhoVal, phiVal, etaVal);
    if (debug_) printf("Cylinder: rho = %4.2f\tphi = %4.2f\teta = %4.2f\n", rhoVal, phiVal, etaVal);


}
//Pre-Job junk

//  
// zero storage
// 
void ResidualRefitting::zero_storage() {
    if (debug_) printf("zero_storage\n");

    zero_muon(&storageGmrOld_   );
    zero_muon(&storageGmrNew_   );
    zero_muon(&storageSamNew_   );
    zero_muon(&storageTrkNew_   );
    zero_muon(&storageGmrNoSt1_ );
    zero_muon(&storageSamNoSt1_ );
    zero_muon(&storageGmrNoSt2_ );
    zero_muon(&storageSamNoSt2_ );
    zero_muon(&storageGmrNoSt3_ );
    zero_muon(&storageSamNoSt3_ );
    zero_muon(&storageGmrNoSt4_ );
    zero_muon(&storageSamNoSt4_ );
//zero out the tracker
    zero_muon(&storageGmrNoPXBLayer1    );
    zero_muon(&storageGmrNoPXBLayer2    );
    zero_muon(&storageGmrNoPXBLayer3    );

    zero_muon(&storageGmrNoPXF  );

    zero_muon(&storageGmrNoTIBLayer1    );
    zero_muon(&storageGmrNoTIBLayer2    );
    zero_muon(&storageGmrNoTIBLayer3    );
    zero_muon(&storageGmrNoTIBLayer4    );

    zero_muon(&storageGmrNoTID  );

    zero_muon(&storageGmrNoTOBLayer1    );
    zero_muon(&storageGmrNoTOBLayer2    );
    zero_muon(&storageGmrNoTOBLayer3    );
    zero_muon(&storageGmrNoTOBLayer4    );
    zero_muon(&storageGmrNoTOBLayer5    );
    zero_muon(&storageGmrNoTOBLayer6    );

    zero_muon(&storageGmrNoTEC  );

    zero_muon(&storageTrkNoPXBLayer1    );
    zero_muon(&storageTrkNoPXBLayer2    );
    zero_muon(&storageTrkNoPXBLayer3    );

    zero_muon(&storageTrkNoPXF  );

    zero_muon(&storageTrkNoTIBLayer1    );
    zero_muon(&storageTrkNoTIBLayer2    );
    zero_muon(&storageTrkNoTIBLayer3    );
    zero_muon(&storageTrkNoTIBLayer4    );

    zero_muon(&storageTrkNoTID  );

    zero_muon(&storageTrkNoTOBLayer1    );
    zero_muon(&storageTrkNoTOBLayer2    );
    zero_muon(&storageTrkNoTOBLayer3    );
    zero_muon(&storageTrkNoTOBLayer4    );
    zero_muon(&storageTrkNoTOBLayer5    );
    zero_muon(&storageTrkNoTOBLayer6    );

    zero_muon(&storageTrkNoTEC  );

    zero_trackExtrap(&storageTrackExtrapRec_        );
    zero_trackExtrap(&storageTrackExtrapTracker_    );
    zero_trackExtrap(&storageTrackExtrapRecNoSt1_   );
    zero_trackExtrap(&storageTrackExtrapRecNoSt2_   );
    zero_trackExtrap(&storageTrackExtrapRecNoSt3_   );
    zero_trackExtrap(&storageTrackExtrapRecNoSt4_   );

    zero_trackExtrap(&trackExtrap120_               );

    zero_trackExtrap(&samExtrap120_                 );

    zero_trackExtrap(&storageTrackNoPXBLayer1       );
    zero_trackExtrap(&storageTrackNoPXBLayer2       );
    zero_trackExtrap(&storageTrackNoPXBLayer3       );

    zero_trackExtrap(&storageTrackNoPXF             );

    zero_trackExtrap(&storageTrackNoTIBLayer1       );
    zero_trackExtrap(&storageTrackNoTIBLayer2       );
    zero_trackExtrap(&storageTrackNoTIBLayer3       );
    zero_trackExtrap(&storageTrackNoTIBLayer4       );

    zero_trackExtrap(&storageTrackNoTOBLayer1       );
    zero_trackExtrap(&storageTrackNoTOBLayer2       );
    zero_trackExtrap(&storageTrackNoTOBLayer3       );
    zero_trackExtrap(&storageTrackNoTOBLayer4       );
    zero_trackExtrap(&storageTrackNoTOBLayer5       );
    zero_trackExtrap(&storageTrackNoTOBLayer6       );

    zero_trackExtrap(&storageTrackNoTEC             );

    zero_trackExtrap(&storageTrackNoTID             );

    storageRecMuon_     .n_ = 0;
    storageTrackHit_    .n_ = 0;

}
//
// Zero out a muon reference
//
void ResidualRefitting::zero_muon(ResidualRefitting::storage_muon* str){

    str->n_ = 0;
    
    for (int i = 0; i < ResidualRefitting::N_MAX_STORED; i++) {
    
        str->pt_            [i] = -9999;
        str->eta_           [i] = -9999;
        str->p_             [i] = -9999;
        str->phi_           [i] = -9999;    
        str->numRecHits_    [i] = -9999;
        str->chiSq_         [i] = -9999;
        str->ndf_           [i] = -9999;
        str->chiSqOvrNdf_   [i] = -9999;
    }

}
//
// Zero track extrapolation
//
void ResidualRefitting::zero_trackExtrap(ResidualRefitting::storage_trackExtrap* str) {
    
    str->n_ = 0;
    for (int i = 0; i < ResidualRefitting::N_MAX_STORED_HIT; i++) {

        str->muonLink_  [i] = -9999;
        str->recLink_   [i] = -9999;
        str->gpX_       [i] = -9999;
        str->gpY_       [i] = -9999;
        str->gpZ_       [i] = -9999;
        str->gpEta_     [i] = -9999;
        str->gpPhi_     [i] = -9999;
        str->lpX_       [i] = -9999;
        str->lpY_       [i] = -9999;
        str->lpZ_       [i] = -9999;
        str->resX_      [i] = -9999;
        str->resY_      [i] = -9999;
        str->resZ_      [i] = -9999;
    }
    
}
//  
// Begin Job
// 
void ResidualRefitting::beginJob() {

    std::cout<<"Creating file "<< outputFileName_.c_str()<<std::endl;

    outputFile_  = new TFile( outputFileName_.c_str(), "RECREATE" ) ; 

    outputTree_ = new TTree("outputTree","outputTree");

    outputTree_->Branch("eventInfo",&eventInfo_,
                        "evtNum_/I:"
                        "runNum_/I"
                );


    ResidualRefitting::branchMuon(storageGmrOld_    , "gmrOld");
    ResidualRefitting::branchMuon(storageGmrNew_    , "gmrNew");
    ResidualRefitting::branchMuon(storageGmrNoSt1_  , "gmrNoSt1");
    ResidualRefitting::branchMuon(storageGmrNoSt2_  , "gmrNoSt2");
    ResidualRefitting::branchMuon(storageGmrNoSt3_  , "gmrNoSt3");
    ResidualRefitting::branchMuon(storageGmrNoSt4_  , "gmrNoSt4");

    ResidualRefitting::branchMuon(storageSamNew_    , "samNew");
    ResidualRefitting::branchMuon(storageSamNoSt1_  , "samNoSt1");
    ResidualRefitting::branchMuon(storageSamNoSt2_  , "samNoSt2");
    ResidualRefitting::branchMuon(storageSamNoSt3_  , "samNoSt3");
    ResidualRefitting::branchMuon(storageSamNoSt4_  , "samNoSt4");

    ResidualRefitting::branchMuon(storageTrkNew_    , "trkNew");
    ResidualRefitting::branchMuon(storageGmrNoPXBLayer1 , "gmrNoPXBLayer1");
    ResidualRefitting::branchMuon(storageGmrNoPXBLayer2 , "gmrNoPXBLayer2");
    ResidualRefitting::branchMuon(storageGmrNoPXBLayer3 , "gmrNoPXBLayer3");
    ResidualRefitting::branchMuon(storageGmrNoPXF   , "gmrNoPXF");
    ResidualRefitting::branchMuon(storageGmrNoTIBLayer1 , "gmrNoTIBLayer1");
    ResidualRefitting::branchMuon(storageGmrNoTIBLayer2 , "gmrNoTIBLayer2");
    ResidualRefitting::branchMuon(storageGmrNoTIBLayer3 , "gmrNoTIBLayer3");
    ResidualRefitting::branchMuon(storageGmrNoTIBLayer4 , "gmrNoTIBLayer4");
    ResidualRefitting::branchMuon(storageGmrNoTID   , "gmrNoTID");
    ResidualRefitting::branchMuon(storageGmrNoTOBLayer1 , "gmrNoTOBLayer1");
    ResidualRefitting::branchMuon(storageGmrNoTOBLayer2 , "gmrNoTOBLayer2");
    ResidualRefitting::branchMuon(storageGmrNoTOBLayer3 , "gmrNoTOBLayer3");
    ResidualRefitting::branchMuon(storageGmrNoTOBLayer4 , "gmrNoTOBLayer4");
    ResidualRefitting::branchMuon(storageGmrNoTOBLayer5 , "gmrNoTOBLayer5");
    ResidualRefitting::branchMuon(storageGmrNoTOBLayer6 , "gmrNoTOBLayer6");
    ResidualRefitting::branchMuon(storageGmrNoTEC   , "gmrNoTEC");

    ResidualRefitting::branchMuon(storageTrkNoPXBLayer1 , "trkNoPXBLayer1");
    ResidualRefitting::branchMuon(storageTrkNoPXBLayer2 , "trkNoPXBLayer2");
    ResidualRefitting::branchMuon(storageTrkNoPXBLayer3 , "trkNoPXBLayer3");
    ResidualRefitting::branchMuon(storageTrkNoPXF   , "trkNoPXF");
    ResidualRefitting::branchMuon(storageTrkNoTIBLayer1 , "trkNoTIBLayer1");
    ResidualRefitting::branchMuon(storageTrkNoTIBLayer2 , "trkNoTIBLayer2");
    ResidualRefitting::branchMuon(storageTrkNoTIBLayer3 , "trkNoTIBLayer3");
    ResidualRefitting::branchMuon(storageTrkNoTIBLayer4 , "trkNoTIBLayer4");
    ResidualRefitting::branchMuon(storageTrkNoTID   , "trkNoTID");
    ResidualRefitting::branchMuon(storageTrkNoTOBLayer1 , "trkNoTOBLayer1");
    ResidualRefitting::branchMuon(storageTrkNoTOBLayer2 , "trkNoTOBLayer2");
    ResidualRefitting::branchMuon(storageTrkNoTOBLayer3 , "trkNoTOBLayer3");
    ResidualRefitting::branchMuon(storageTrkNoTOBLayer4 , "trkNoTOBLayer4");
    ResidualRefitting::branchMuon(storageTrkNoTOBLayer5 , "trkNoTOBLayer5");
    ResidualRefitting::branchMuon(storageTrkNoTOBLayer6 , "trkNoTOBLayer6");
    ResidualRefitting::branchMuon(storageTrkNoTEC   , "trkNoTEC");
                    
    outputBranch_ = outputTree_ -> Branch("recHitsNew", &storageRecMuon_, 

        "n_/I:"
        "muonLink_[1000]/I:"
        
        "system_[1000]/I:"
        "endcap_[1000]/I:"
        "station_[1000]/I:"
        "ring_[1000]/I:"
        "chamber_[1000]/I:"
        "layer_[1000]/I:"
        "superLayer_[1000]/I:"
        "wheel_[1000]/I:"
        "sector_[1000]/I:"
        
    
        "gpX_[1000]/F:"
        "gpY_[1000]/F:"
        "gpZ_[1000]/F:"
        "gpEta_[1000]/F:"
        "gpPhi_[1000]/F:"
        "lpX_[1000]/F:"
        "lpY_[1000]/F:"
        "lpZ_[1000]/F"
        );
        
        
    outputBranch_ = outputTree_ -> Branch("recHitsTracker", &storageTrackHit_,
    
        "n_/I:"
        
        "muonLink_[1000]/I:"
        "detector_[1000]/I:"
        "subdetector_[1000]/I:"
        "blade_[1000]/I:"
        "disk_[1000]/I:"
        "ladder_[1000]/I:"
        "layer_[1000]/I:"
        "module_[1000]/I:"
        "panel_[1000]/I:"
        "ring_[1000]/I:"
        "side_[1000]/I:"
        "wheel_[1000]/I:"
                
        "gpX_[1000]/F:"
        "gpY_[1000]/F:"
        "gpZ_[1000]/F:"
        "gpEta_[1000]/F:"
        "gpPhi_[1000]/F:"
        "lpX_[1000]/F:"
        "lpY_[1000]/F:"
        "lpZ_[1000]/F"
        );  
        
        
    ResidualRefitting::branchTrackExtrap(storageTrackExtrapRec_     , "trkExtrap");
    ResidualRefitting::branchTrackExtrap(storageTrackExtrapRecNoSt1_    , "trkExtrapNoSt1");
    ResidualRefitting::branchTrackExtrap(storageTrackExtrapRecNoSt2_    , "trkExtrapNoSt2");
    ResidualRefitting::branchTrackExtrap(storageTrackExtrapRecNoSt3_    , "trkExtrapNoSt3");
    ResidualRefitting::branchTrackExtrap(storageTrackExtrapRecNoSt4_    , "trkExtrapNoSt4");

    ResidualRefitting::branchTrackExtrap(storageTrackExtrapTracker_, "trkExtrapTracker");
    ResidualRefitting::branchTrackExtrap(storageTrackNoPXF          , "trkExtrapNoPXF");
    ResidualRefitting::branchTrackExtrap(storageTrackNoPXBLayer1            , "trkExtrapNoPXBLayer1");
    ResidualRefitting::branchTrackExtrap(storageTrackNoPXBLayer2            , "trkExtrapNoPXBLayer2");
    ResidualRefitting::branchTrackExtrap(storageTrackNoPXBLayer3            , "trkExtrapNoPXBLayer3");
    ResidualRefitting::branchTrackExtrap(storageTrackNoTIBLayer1            , "trkExtrapNoTIBLayer1");
    ResidualRefitting::branchTrackExtrap(storageTrackNoTIBLayer2            , "trkExtrapNoTIBLayer2");
    ResidualRefitting::branchTrackExtrap(storageTrackNoTIBLayer3            , "trkExtrapNoTIBLayer3");
    ResidualRefitting::branchTrackExtrap(storageTrackNoTIBLayer4            , "trkExtrapNoTIBLayer4");
    ResidualRefitting::branchTrackExtrap(storageTrackNoTID          , "trkExtrapNoTID");
    ResidualRefitting::branchTrackExtrap(storageTrackNoTOBLayer1            , "trkExtrapNoTOBLayer1");
    ResidualRefitting::branchTrackExtrap(storageTrackNoTOBLayer2            , "trkExtrapNoTOBLayer2");
    ResidualRefitting::branchTrackExtrap(storageTrackNoTOBLayer3            , "trkExtrapNoTOBLayer3");
    ResidualRefitting::branchTrackExtrap(storageTrackNoTOBLayer4            , "trkExtrapNoTOBLayer4");
    ResidualRefitting::branchTrackExtrap(storageTrackNoTOBLayer5            , "trkExtrapNoTOBLayer5");
    ResidualRefitting::branchTrackExtrap(storageTrackNoTOBLayer6            , "trkExtrapNoTOBLayer6");
    ResidualRefitting::branchTrackExtrap(storageTrackNoTEC          , "trkExtrapNoTEC");

    ResidualRefitting::branchTrackExtrap(trackExtrap120_            , "trackCyl120");
    ResidualRefitting::branchTrackExtrap(samExtrap120_              , "samCyl120");

}
//
// Set the Muon Branches
//
void ResidualRefitting::branchMuon(ResidualRefitting::storage_muon& storageTmp, std::string branchName){

    outputBranch_ = outputTree_ -> Branch(branchName.c_str(), &storageTmp, 
                                    "n_/I:"
                                    "charge_[10]/I:"
                                    "pt_[10]/F:"
                                    "eta_[10]/F:"
                                    "p_[10]/F:"
                                    "phi_[10]/F:"
                                    "numRecHits_[10]/I:"
                                    "chiSq_[10]/F:"
                                    "ndf_[10]/F:"
                                    "chiSqOvrNdf_[10]/F"
                                        
                                        );

}
//
// Set the Branches for Track Extrapolations
//
void ResidualRefitting::branchTrackExtrap(ResidualRefitting::storage_trackExtrap& storageTmp, std::string branchName){

    outputBranch_ = outputTree_ -> Branch(branchName.c_str(), &storageTmp, 
                                    "n_/I:"
                                    "muonLink_[1000]/I:"
                                    "recLink_[1000]/I:"
                                    "gpX_[1000]/F:"
                                    "gpY_[1000]/F:"
                                    "gpZ_[1000]/F:"
                                    "gpEta_[1000]/F:"
                                    "gpPhi_[1000]/F:"
                                    "lpX_[1000]/F:"
                                    "lpY_[1000]/F:"
                                    "lpZ_[1000]/F:"
                                    "resX_[1000]/F:"
                                    "resY_[1000]/F:"
                                    "resZ_[1000]/F"
                                        
                                    );

}
// 
// End Job
// 
void ResidualRefitting::endJob() {


  outputFile_ -> Write();

  outputFile_ -> Close() ;

}
// 
// Return a Free Trajectory state for a muon track
// 
FreeTrajectoryState ResidualRefitting::freeTrajStateMuon(reco::TrackRef muon){
                    
            math::XYZPoint  innerPos = muon -> referencePoint();
            math::XYZVector innerMom = muon -> momentum();
            if (debug_) std::cout   <<  "Inner Pos: "
                                    <<  "\tx = "    << innerPos.X()
                                    <<  "\ty = "    << innerPos.Y()
                                    <<  "\tz = "    << innerPos.Z()
                                    <<  std::endl;
                                
            GlobalPoint   innerPoint( innerPos.X(), innerPos.Y(), innerPos.Z());
            GlobalVector  innerVec  ( innerMom.X(), innerMom.Y(), innerMom.Z());

            FreeTrajectoryState recoStart( innerPoint, innerVec, muon ->charge(), theField ); 
            return recoStart;

}

// nTuple value Dumps

// 
// dump Track Extrapolation
//
void ResidualRefitting::dumpTrackExtrap(const ResidualRefitting::storage_trackExtrap& track) {
    std::cout<<"\n\nExtrapolation Dump:\n";
    for (unsigned int i = 0; i < (unsigned int)track.n_; i++) {
        
//      double rho = sqrt( (float)track.gpX_[i] * (float)track.gpX_[i] + (float)track.gpY_[i] * (float)track.gpY_[i]  );

        printf ("%d\tmuonLink= %d",i, (int)track.muonLink_[i]);     
        printf ("\trecLink = %d",   (int)track.recLink_[i]  );
//      printf ("\tGlobal\tx = %0.3f"       ,       (float)track.gpX_[i]    );
//      printf ("\ty = %0.3f"       ,       (float)track.gpY_[i]    );
//      printf ("\tz = %0.3f"       ,       (float)track.gpZ_[i]    );
//      printf ("\trho =%0.3f"      ,       rho                     );
//      printf ("\teta = %0.3f"     ,       (float)track.gpEta_[i]  );
//      printf ("\tphi = %0.3f"     ,       (float)track.gpPhi_[i]  );
        printf ("\t\tLocal\tx = %0.3f"  ,       (float)track.lpX_[i]    );
        printf ("\ty = %0.3f"       ,   (float)track.lpY_[i]        );
        printf ("\tz = %0.3f\n" ,   (float)track.lpZ_[i]        );
        
    }
    
}
//
// dump Muon Rec Hits
//
void ResidualRefitting::dumpMuonRecHits(const ResidualRefitting::storage_hit& hit) {
    std::cout<<"Muon Rec Hits Dump:\n";
    for (unsigned int i = 0; i < (unsigned int)hit.n_; i++) {
        
//      double rho = sqrt( (float)hit.gpX_[i] * (float)hit.gpX_[i] + (float)hit.gpY_[i] * (float)hit.gpY_[i]  );
        
        printf ("%d\tsubdetector = %d\t superLayer =%d" , i,    (int)hit.system_[i], (int)hit.superLayer_[i]    );
//      printf ("\tGlobal\tx = %0.3f"           ,       (float)hit.gpX_[i]          );
//      printf ("\ty = %0.3f"               ,       (float)hit.gpY_[i]          );
//      printf ("\tz = %0.3f"               ,       (float)hit.gpZ_[i]          );
//      printf ("\trho =%0.3f"              ,       rho                         );
//      printf ("\teta = %0.3f"             ,       (float)hit.gpEta_[i]        );
//      printf ("\tphi = %0.3f\n"           ,       (float)hit.gpPhi_[i]        );
        printf ("\t\tLocal\tx = %0.3f"      ,       (float)hit.lpX_[i]          );
        printf ("\ty = %0.3f"               ,       (float)hit.lpY_[i]          );
        printf ("\tz = %0.3f\n"         ,       (float)hit.lpZ_[i]          );
        
    }
    
}
//
// dump Tracker Rec Hits
//
void ResidualRefitting::dumpTrackHits(const ResidualRefitting::storage_trackHit& hit) {
    std::cout<<"Tracker Rec Hits Dump:\n";
    for (unsigned int i = 0; i < (unsigned int)hit.n_; i++) {
        
//      double rho = sqrt( (float)hit.gpX_[i] * (float)hit.gpX_[i] + (float)hit.gpY_[i] * (float)hit.gpY_[i]  );
        
        printf ("%d\tsubdetector = %d"      , i,    (int)hit.subdetector_[i]    );
        printf ("\tlayer = %d"              ,   (int)hit.layer_[i]  );
//      printf ("\tGlobal\tx = %0.3f"           ,       (float)hit.gpX_[i]          );
//      printf ("\ty = %0.3f"               ,       (float)hit.gpY_[i]          );
//      printf ("\tz = %0.3f"               ,       (float)hit.gpZ_[i]          );
//      printf ("\trho =%0.3f"              ,       rho                         );
//      printf ("\teta = %0.3f"             ,       (float)hit.gpEta_[i]        );
//      printf ("\tphi = %0.3f\n"           ,       (float)hit.gpPhi_[i]        );
        printf ("\t\tLocal\tx = %0.3f"      ,       (float)hit.lpX_[i]          );
        printf ("\ty = %0.3f"               ,       (float)hit.lpY_[i]          );
        printf ("\tz = %0.3f\n"             ,       (float)hit.lpZ_[i]          );
        
    }
    
}
//
//Dump a TrackRef 
//
void ResidualRefitting::dumpTrackRef(reco::TrackRef muon, std::string str) {

    float pt = muon->pt();
    float p  = muon->p  ();
    float eta = muon->eta();
    float phi = muon->phi();
    printf("\t%s: \tp = %4.2f \t pt = %4.2f \t eta = %4.2f \t phi = %4.2f\n",str.c_str(), p, pt, eta, phi);

}

DEFINE_FWK_MODULE( ResidualRefitting );


//Deprecated