Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:45:55

0001 #include <iostream>
0002 #include "ResidualRefitting.h"
0003 #include <iomanip>
0004 
0005 //framework includes
0006 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0007 #include "DataFormats/MuonDetId/interface/DTWireId.h"
0008 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0009 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0010 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0011 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0012 
0013 #include "FWCore/Framework/interface/ConsumesCollector.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/Utilities/interface/EDMException.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/Framework/interface/MakerMacros.h"
0019 
0020 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0021 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0022 #include "Geometry/CommonDetUnit/interface/GeomDetEnumerators.h"
0023 
0024 #include "RecoMuon/DetLayers/interface/MuonDetLayerGeometry.h"
0025 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBuilder.h"
0026 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
0027 
0028 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
0029 
0030 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0031 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0032 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0033 
0034 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0035 
0036 ResidualRefitting::ResidualRefitting(const edm::ParameterSet& cfg)
0037     : magFieldToken_(esConsumes()),
0038       topoToken_(esConsumes()),
0039       trackingGeometryToken_(esConsumes()),
0040       propagatorToken_(esConsumes(edm::ESInputTag("", cfg.getParameter<std::string>("propagator")))),
0041       outputFileName_(cfg.getUntrackedParameter<std::string>("histoutputFile")),
0042       muons_(cfg.getParameter<edm::InputTag>("muons")),
0043       muonsRemake_(cfg.getParameter<edm::InputTag>("muonsRemake")),  //This Feels Misalignment
0044       muonsNoStation1_(cfg.getParameter<edm::InputTag>("muonsNoStation1")),
0045       muonsNoStation2_(cfg.getParameter<edm::InputTag>("muonsNoStation2")),
0046       muonsNoStation3_(cfg.getParameter<edm::InputTag>("muonsNoStation3")),
0047       muonsNoStation4_(cfg.getParameter<edm::InputTag>("muonsNoStation4")),
0048 
0049       /*
0050   muonsNoPXBLayer1_     ( cfg.getParameter<edm::InputTag>("muonsNoPXBLayer1"    ) ),
0051   muonsNoPXBLayer2_     ( cfg.getParameter<edm::InputTag>("muonsNoPXBLayer1"    ) ),
0052   muonsNoPXBLayer3_     ( cfg.getParameter<edm::InputTag>("muonsNoPXBLayer1"    ) ),
0053 
0054   muonsNoTIBLayer1_         ( cfg.getParameter<edm::InputTag>("muonsNoTIBLayer1"    ) ),
0055   muonsNoTIBLayer2_         ( cfg.getParameter<edm::InputTag>("muonsNoTIBLayer2"    ) ),
0056   muonsNoTIBLayer3_         ( cfg.getParameter<edm::InputTag>("muonsNoTIBLayer3"    ) ),
0057   muonsNoTIBLayer4_         ( cfg.getParameter<edm::InputTag>("muonsNoTIBLayer4"    ) ),
0058 
0059   muonsNoTOBLayer1_         ( cfg.getParameter<edm::InputTag>("muonsNoTOBLayer1"    ) ),
0060   muonsNoTOBLayer2_         ( cfg.getParameter<edm::InputTag>("muonsNoTOBLayer2"    ) ),
0061   muonsNoTOBLayer3_         ( cfg.getParameter<edm::InputTag>("muonsNoTOBLayer3"    ) ),
0062   muonsNoTOBLayer4_         ( cfg.getParameter<edm::InputTag>("muonsNoTOBLayer4"    ) ),
0063   muonsNoTOBLayer5_         ( cfg.getParameter<edm::InputTag>("muonsNoTOBLayer5"    ) ),
0064   muonsNoTOBLayer6_         ( cfg.getParameter<edm::InputTag>("muonsNoTOBLayer6"    ) ),*/
0065       debug_(cfg.getUntrackedParameter<bool>("doDebug")),
0066       outputFile_(nullptr),
0067       outputTree_(nullptr),
0068       outputBranch_(nullptr),
0069       theField(nullptr) {
0070   eventInfo_.evtNum_ = 0;
0071   eventInfo_.evtNum_ = 0;
0072 
0073   // service parameters
0074   edm::ParameterSet serviceParameters = cfg.getParameter<edm::ParameterSet>("ServiceParameters");
0075 
0076   // the services
0077   theService = new MuonServiceProxy(serviceParameters, consumesCollector());
0078 
0079 }  //The constructor
0080 
0081 void ResidualRefitting::analyze(const edm::Event& event, const edm::EventSetup& eventSetup) {
0082   if (debug_)
0083     printf("STARTING EVENT\n");
0084 
0085   eventInfo_.evtNum_ = (int)event.id().run();
0086   eventInfo_.runNum_ = (int)event.id().event();
0087 
0088   // Generator Collection
0089 
0090   // The original muon collection that is sitting in memory
0091   edm::Handle<reco::MuonCollection> muons;
0092 
0093   edm::Handle<reco::TrackCollection> muonTracks;
0094   edm::Handle<reco::TrackCollection> muonsNoSt1;
0095   edm::Handle<reco::TrackCollection> muonsNoSt2;
0096   edm::Handle<reco::TrackCollection> muonsNoSt3;
0097   edm::Handle<reco::TrackCollection> muonsNoSt4;
0098 
0099   event.getByLabel(muons_, muons);  //set label to muons
0100   event.getByLabel(muonsRemake_, muonTracks);
0101   event.getByLabel(muonsNoStation1_, muonsNoSt1);
0102   event.getByLabel(muonsNoStation2_, muonsNoSt2);
0103   event.getByLabel(muonsNoStation3_, muonsNoSt3);
0104   event.getByLabel(muonsNoStation4_, muonsNoSt4);
0105 
0106   /*
0107 //  std::cout<<"Muon Collection No PXB "<<std::endl;
0108 //Tracker Barrel Pixel Refits
0109     edm::Handle<MuonCollection> muonsNoPXBLayer1Coll;
0110     event.getByLabel(muonsNoPXBLayer1_, muonsNoPXBLayer1Coll);
0111     edm::Handle<MuonCollection> muonsNoPXBLayer2Coll;
0112     event.getByLabel(muonsNoPXBLayer2_, muonsNoPXBLayer2Coll);
0113     edm::Handle<MuonCollection> muonsNoPXBLayer3Coll;
0114     event.getByLabel(muonsNoPXBLayer3_, muonsNoPXBLayer3Coll);
0115 //  std::cout<<"Muon Collection No TIB "<<std::endl;
0116 // Tracker Inner Barrel Refits
0117     edm::Handle<MuonCollection> muonsNoTIBLayer1Coll;
0118     event.getByLabel(muonsNoTIBLayer1_, muonsNoTIBLayer1Coll);
0119     edm::Handle<MuonCollection> muonsNoTIBLayer2Coll;
0120     event.getByLabel(muonsNoTIBLayer2_, muonsNoTIBLayer2Coll);
0121     edm::Handle<MuonCollection> muonsNoTIBLayer3Coll;
0122     event.getByLabel(muonsNoTIBLayer3_, muonsNoTIBLayer3Coll);
0123     edm::Handle<MuonCollection> muonsNoTIBLayer4Coll;
0124     event.getByLabel(muonsNoTIBLayer4_, muonsNoTIBLayer4Coll);
0125 
0126 //  std::cout<<"Muon Collection No TOB "<<std::endl;
0127 
0128 //Tracker outer barrel refits
0129     edm::Handle<MuonCollection> muonsNoTOBLayer1Coll;
0130     event.getByLabel(muonsNoTOBLayer1_, muonsNoTOBLayer1Coll);
0131     edm::Handle<MuonCollection> muonsNoTOBLayer2Coll;
0132     event.getByLabel(muonsNoTOBLayer2_, muonsNoTOBLayer2Coll);
0133     edm::Handle<MuonCollection> muonsNoTOBLayer3Coll;
0134     event.getByLabel(muonsNoTOBLayer3_, muonsNoTOBLayer3Coll);
0135     edm::Handle<MuonCollection> muonsNoTOBLayer4Coll;
0136     event.getByLabel(muonsNoTOBLayer4_, muonsNoTOBLayer4Coll);
0137     edm::Handle<MuonCollection> muonsNoTOBLayer5Coll;
0138     event.getByLabel(muonsNoTOBLayer5_, muonsNoTOBLayer5Coll);
0139     edm::Handle<MuonCollection> muonsNoTOBLayer6Coll;
0140     event.getByLabel(muonsNoTOBLayer6_, muonsNoTOBLayer6Coll);
0141 */
0142   //magnetic field information
0143   theField = &eventSetup.getData(magFieldToken_);
0144   edm::ESHandle<GlobalTrackingGeometry> globalTrackingGeometry = eventSetup.getHandle(trackingGeometryToken_);
0145   thePropagator = eventSetup.getHandle(propagatorToken_);
0146   theService->update(eventSetup);
0147 
0148   //Zero storage
0149   zero_storage();
0150 
0151   //Do the Gmr Muons from the unModified Collection
0152 
0153   /*
0154     int iGmr = 0;
0155     if ( (muons->end() - muons->begin()) > 0) printf("Data Dump:: Original GMR Muons\n");
0156     for ( MuonCollection::const_iterator muon = muons->begin(); muon!=muons->end(); muon++, iGmr++) {
0157         if ( iGmr >= ResidualRefitting::N_MAX_STORED) break; // error checking
0158         if (!debug
0159         
0160         dumpTrackRef(muon->combinedMuon(), "cmb"); 
0161         dumpTrackRef(muon->standAloneMuon(), "sam");
0162         dumpTrackRef(muon->track(), "trk");
0163         
0164 
0165     }
0166     storageGmrOld_.n_ = iGmr;
0167     storageSamNew_.n_ = iGmr;
0168 */
0169 
0170   //Refitted muons
0171   if (debug_)
0172     printf("Data Dump:: Rebuilt GMR Muon Track With TeV refitter default\n");
0173   int iGmrRemake = 0;
0174   for (reco::TrackCollection::const_iterator muon = muonTracks->begin(); muon != muonTracks->end();
0175        muon++, iGmrRemake++) {
0176     if (iGmrRemake >= ResidualRefitting::N_MAX_STORED)
0177       break;  // error checking
0178               // from TrackInfoProducer/test/TrackInfoAnalyzerExample.cc
0179     reco::TrackRef trackref = reco::TrackRef(muonTracks, iGmrRemake);
0180 
0181     if (debug_)
0182       dumpTrackRef(trackref, "gmr");
0183     muonInfo(storageGmrNew_, trackref, iGmrRemake);
0184   }
0185   storageGmrNew_.n_ = iGmrRemake;
0186 
0187   if (debug_)
0188     printf("muons Remake");
0189   if (debug_)
0190     printf("-----------------------------------------\n");
0191   CollectTrackHits(muonTracks, storageTrackExtrapRec_, eventSetup);
0192 
0193   if (true) {
0194     printf("muons No Station 1");
0195     printf("-----------------------------------------\n");
0196   }
0197   NewTrackMeasurements(muonTracks, muonsNoSt1, storageTrackExtrapRecNoSt1_);
0198 
0199   if (true) {
0200     printf("muons No Station 2");
0201     printf("-----------------------------------------\n");
0202   }
0203   NewTrackMeasurements(muonTracks, muonsNoSt2, storageTrackExtrapRecNoSt2_);
0204 
0205   if (true) {
0206     printf("muons No Station 3");
0207     printf("-----------------------------------------\n");
0208   }
0209   NewTrackMeasurements(muonTracks, muonsNoSt3, storageTrackExtrapRecNoSt3_);
0210 
0211   if (true) {
0212     printf("muons No Station 4");
0213     printf("-----------------------------------------\n");
0214   }
0215   NewTrackMeasurements(muonTracks, muonsNoSt4, storageTrackExtrapRecNoSt4_);
0216 
0217   //    dumpMuonRecHits(storageRecMuon_);
0218 
0219   /****************************************************************************************************************************************/
0220 
0221   /*
0222  *  extrapolates track to a cylinder.
0223  *  commented for cosmic runs with no tracker in reco muons!!
0224  *
0225 */
0226 
0227   int iGmrCyl = 0;
0228   for (reco::MuonCollection::const_iterator muon = muons->begin(); muon != muons->end(); muon++, iGmrCyl++) {
0229     dumpTrackRef(muon->combinedMuon(), "cmb");
0230     dumpTrackRef(muon->standAloneMuon(), "sam");
0231     dumpTrackRef(muon->track(), "trk");
0232 
0233     cylExtrapTrkSam(iGmrCyl, muon->standAloneMuon(), samExtrap120_, 120.);
0234     cylExtrapTrkSam(iGmrCyl, muon->track(), trackExtrap120_, 120.);
0235   }
0236   samExtrap120_.n_ = iGmrCyl;
0237   trackExtrap120_.n_ = iGmrCyl;
0238 
0239   if (iGmrRemake > 0 || iGmrCyl > 0) {
0240     outputTree_->Fill();
0241     std::cout << "FILLING NTUPLE!" << std::endl;
0242     std::cout << "Entries Recorded: " << outputTree_->GetEntries() << " Branch :: " << outputBranch_->GetEntries()
0243               << std::endl
0244               << std::endl;
0245   } else
0246     std::cout << "no tracks -- no fill!\n" << std::endl << std::endl;
0247 
0248   //  /*************************************************************************************************************/
0249   //  //END OF ntuple dumper
0250   //  //END OF ntuple dumper
0251   //  /***********************************************************************************************************/
0252 }
0253 //end Analyze() main function
0254 
0255 //------------------------------------------------------------------------------
0256 //
0257 // Destructor
0258 //
0259 ResidualRefitting::~ResidualRefitting() {
0260   delete outputFile_;
0261   delete theService;
0262 }
0263 //
0264 // Track Collection Analysis
0265 //
0266 void ResidualRefitting::CollectTrackHits(edm::Handle<reco::TrackCollection> trackColl,
0267                                          ResidualRefitting::storage_trackExtrap& trackExtrap,
0268                                          const edm::EventSetup& eventSetup) {
0269   //Retrieve tracker topology from geometry
0270   const TrackerTopology* const tTopo = &eventSetup.getData(topoToken_);
0271 
0272   int iMuonHit = 0;
0273   int iTrackHit = 0;
0274   int numTracks = 0;
0275 
0276   for (reco::TrackCollection::const_iterator muon = trackColl->begin(); muon != trackColl->end(); muon++) {
0277     int iTrack = muon - trackColl->begin();
0278     reco::TrackRef trackref = reco::TrackRef(trackColl, iTrack);
0279     FreeTrajectoryState recoStart = ResidualRefitting::freeTrajStateMuon(trackref);
0280 
0281     if (debug_)
0282       dumpTrackRef(trackref, "CollectTrackHits Track");
0283 
0284     int iRec = 0;
0285     for (auto const& rec : muon->recHits()) {
0286       DetId detid = rec->geographicalId();
0287 
0288       if (detid.det() != DetId::Muon && detid.det() != DetId::Tracker) {
0289         if (debug_)
0290           printf("Rec Hit not from muon system or tracker... continuing...\n");
0291         continue;
0292       }
0293       //            numTracks++;
0294       // Get Local and Global Position of Hits
0295 
0296       LocalPoint lp = rec->localPosition();
0297       float lpX = lp.x();
0298       float lpY = lp.y();
0299       float lpZ = lp.z();
0300 
0301       auto mrhp = MuonTransientTrackingRecHit::specificBuild(
0302           theService->trackingGeometry()->idToDet(rec->geographicalId()), rec);
0303 
0304       GlobalPoint gp = mrhp->globalPosition();
0305       float gpRecX = gp.x();
0306       float gpRecY = gp.y();
0307       float gpRecZ = gp.z();
0308       float gpRecEta = gp.eta();
0309       float gpRecPhi = gp.phi();
0310 
0311       if (detid.det() == DetId::Muon) {
0312         int systemMuon = detid.subdetId();  // 1 DT; 2 CSC; 3 RPC
0313         int endcap = -999;
0314         int station = -999;
0315         int ring = -999;
0316         int chamber = -999;
0317         int layer = -999;
0318         int superLayer = -999;
0319         int wheel = -999;
0320         int sector = -999;
0321         if (systemMuon == MuonSubdetId::CSC) {
0322           CSCDetId id(detid.rawId());
0323           endcap = id.endcap();
0324           station = id.station();
0325           ring = id.ring();
0326           chamber = id.chamber();
0327           layer = id.layer();
0328           if (debug_)
0329             printf("CSC\t[endcap][station][ringN][chamber][layer]:[%d][%d][%d][%d][%d]\t",
0330                    endcap,
0331                    station,
0332                    ring,
0333                    chamber,
0334                    layer);
0335 
0336         } else if (systemMuon == MuonSubdetId::DT) {
0337           DTWireId id(detid.rawId());
0338           station = id.station();
0339           layer = id.layer();
0340           superLayer = id.superLayer();
0341           wheel = id.wheel();
0342           sector = id.sector();
0343           if (debug_)
0344             printf("DT \t[station][layer][superlayer]:[%d][%d][%d]\n", station, layer, superLayer);
0345 
0346         } else if (systemMuon == MuonSubdetId::RPC) {
0347           RPCDetId id(detid.rawId());
0348           station = id.station();
0349           if (debug_)
0350             printf("RPC\t[station]:[%d]\n", station);
0351         }
0352 
0353         storageRecMuon_.muonLink_[iMuonHit] = iTrack;
0354         storageRecMuon_.system_[iMuonHit] = systemMuon;
0355         storageRecMuon_.endcap_[iMuonHit] = endcap;
0356         storageRecMuon_.station_[iMuonHit] = station;
0357         storageRecMuon_.ring_[iMuonHit] = ring;
0358         storageRecMuon_.chamber_[iMuonHit] = chamber;
0359         storageRecMuon_.layer_[iMuonHit] = layer;
0360         storageRecMuon_.superLayer_[iMuonHit] = superLayer;
0361         storageRecMuon_.wheel_[iMuonHit] = wheel;
0362         storageRecMuon_.sector_[iMuonHit] = sector;
0363 
0364         storageRecMuon_.gpX_[iMuonHit] = gpRecX;
0365         storageRecMuon_.gpY_[iMuonHit] = gpRecY;
0366         storageRecMuon_.gpZ_[iMuonHit] = gpRecZ;
0367         storageRecMuon_.gpEta_[iMuonHit] = gpRecEta;
0368         storageRecMuon_.gpPhi_[iMuonHit] = gpRecPhi;
0369         storageRecMuon_.lpX_[iMuonHit] = lpX;
0370         storageRecMuon_.lpY_[iMuonHit] = lpY;
0371         storageRecMuon_.lpZ_[iMuonHit] = lpZ;
0372         iMuonHit++;
0373 
0374       } else if (detid.det() == DetId::Tracker) {
0375         if (debug_)
0376           printf("Tracker\n");
0377 
0378         StoreTrackerRecHits(detid, tTopo, iTrack, iTrackHit);
0379 
0380         storageTrackHit_.gpX_[iTrackHit] = gpRecX;
0381         storageTrackHit_.gpY_[iTrackHit] = gpRecY;
0382         storageTrackHit_.gpZ_[iTrackHit] = gpRecZ;
0383         storageTrackHit_.gpEta_[iTrackHit] = gpRecEta;
0384         storageTrackHit_.gpPhi_[iTrackHit] = gpRecPhi;
0385         storageTrackHit_.lpX_[iTrackHit] = lpX;
0386         storageTrackHit_.lpY_[iTrackHit] = lpY;
0387         storageTrackHit_.lpZ_[iTrackHit] = lpZ;
0388         iTrackHit++;
0389       } else
0390         printf("THIS CAN NOT HAPPEN\n");
0391 
0392       trkExtrap(detid, numTracks, iTrack, iRec, recoStart, lp, trackExtrap);
0393       numTracks++;
0394 
0395       if (debug_)
0396         printf("\tLocal Positon:  \tx = %2.2f\ty = %2.2f\tz = %2.2f\n", lpX, lpY, lpZ);
0397       if (debug_)
0398         printf("\tGlobal Position: \tx = %6.2f\ty = %6.2f\tz = %6.2f\teta = %4.2f\tphi = %3.2f\n",
0399                gpRecX,
0400                gpRecY,
0401                gpRecZ,
0402                gpRecEta,
0403                gpRecPhi);
0404 
0405       ++iRec;
0406     }
0407   }
0408 
0409   storageRecMuon_.n_ = iMuonHit;
0410   storageTrackHit_.n_ = iTrackHit;
0411   trackExtrap.n_ = numTracks;
0412 }
0413 //
0414 // Deal with Re-Fitted Track with some station omitted.
0415 //
0416 // This should take the new track, match it to its track before refitting with the hits dumped, and extrapolate out to the
0417 //  rec hits that were removed from the fit.
0418 //
0419 //
0420 
0421 void ResidualRefitting::NewTrackMeasurements(edm::Handle<reco::TrackCollection> trackCollOrig,
0422                                              edm::Handle<reco::TrackCollection> trackColl,
0423                                              ResidualRefitting::storage_trackExtrap& trackExtrap) {
0424   int numTracks = 0;
0425   int recCounter = 0;
0426 
0427   for (reco::TrackCollection::const_iterator muon = trackColl->begin(); muon != trackColl->end(); muon++) {
0428     int iTrack = muon - trackColl->begin();
0429 
0430     reco::TrackRef trackref = reco::TrackRef(trackColl, iTrack);
0431     FreeTrajectoryState recoStart = ResidualRefitting::freeTrajStateMuon(trackref);
0432 
0433     int iTrackLink = MatchTrackWithRecHits(muon, trackCollOrig);
0434     reco::TrackRef ref = reco::TrackRef(trackCollOrig, iTrackLink);
0435 
0436     for (auto const& rec1 : ref->recHits()) {
0437       bool unbiasedRec = true;
0438 
0439       for (auto const& rec2 : muon->recHits()) {
0440         if (IsSameHit(*rec1, *rec2)) {
0441           unbiasedRec = false;
0442           break;
0443         }
0444       }
0445       if (!unbiasedRec)
0446         continue;
0447 
0448       DetId detid = rec1->geographicalId();
0449 
0450       auto mrhp = MuonTransientTrackingRecHit::specificBuild(
0451           theService->trackingGeometry()->idToDet(rec1->geographicalId()), rec1);
0452 
0453       trkExtrap(detid, numTracks, iTrackLink, recCounter, recoStart, rec1->localPosition(), trackExtrap);
0454       numTracks++;
0455     }
0456   }
0457 
0458   trackExtrap.n_ = numTracks;
0459 }
0460 //
0461 // Find the original track that corresponds to the re-fitted track
0462 //
0463 int ResidualRefitting::MatchTrackWithRecHits(reco::TrackCollection::const_iterator trackIt,
0464                                              edm::Handle<reco::TrackCollection> ref) {
0465   if (debug_)
0466     printf("Matching a re-fitted track to the original track.\n");
0467 
0468   int TrackMatch = -1;
0469 
0470   for (auto const& rec : trackIt->recHits()) {
0471     bool foundMatch = false;
0472     for (reco::TrackCollection::const_iterator refIt = ref->begin(); refIt != ref->end(); refIt++) {
0473       int iTrackMatch = refIt - ref->begin();
0474       if (foundMatch && TrackMatch != iTrackMatch)
0475         break;
0476       for (auto const& recRef : refIt->recHits()) {
0477         if (!IsSameHit(*rec, *recRef))
0478           continue;
0479 
0480         foundMatch = true;
0481         TrackMatch = iTrackMatch;
0482         //  printf("Rec hit match for original track %d\n", iTrackMatch);
0483       }
0484     }
0485     if (!foundMatch) {
0486       printf("SOMETHING WENT WRONG! Could not match Track with original track!");
0487       exit(1);
0488     }
0489   }
0490   if (debug_)
0491     printf("Rec hit match for original track %d\n", TrackMatch);
0492 
0493   //    reco::TrackRef trackref=reco::TrackRef(ref,TrackMatch);
0494   return TrackMatch;
0495 }
0496 /*
0497 //
0498 // Match two tracks to see if one is a subset of the other
0499 //
0500 
0501 bool ResidualRefitting::TrackSubset(reco::TrackRef trackSub, reco::TrackRef trackTop) {
0502 
0503     
0504     bool matchAll = true;
0505 
0506     for (trackingRecHit_iterator recSub = trackSub->recHits().begin(); recSub!=trackSub->recHits().end(); recSub++) {
0507 
0508         bool matchSub = false;
0509 
0510 
0511         for (trackingRecHit_iterator recTop = trackTop->recHits().begin(); recTop!=trackTop->recHits().end(); recTop++) {
0512         
0513             if ( recSub == recTop ) matchSub = true;
0514             if (matchSub) break;
0515 
0516         }
0517         if (!matchSub) return false;
0518 
0519     }
0520 
0521     return matchAll;
0522 
0523 }
0524 */
0525 
0526 //
0527 // Check to see if the rec hits are the same
0528 //
0529 bool ResidualRefitting::IsSameHit(TrackingRecHit const& hit1, TrackingRecHit const& hit2) {
0530   double lpx1 = hit1.localPosition().x();
0531   double lpy1 = hit1.localPosition().y();
0532   double lpz1 = hit1.localPosition().z();
0533 
0534   double lpx2 = hit2.localPosition().x();
0535   double lpy2 = hit2.localPosition().y();
0536   double lpz2 = hit2.localPosition().z();
0537   if (fabs(lpx1 - lpx2) > 1e-3)
0538     return false;
0539   //    printf("Match lpx...\n");
0540   if (fabs(lpy1 - lpy2) > 1e-3)
0541     return false;
0542   //    printf("Match lpy...\n");
0543   if (fabs(lpz1 - lpz2) > 1e-3)
0544     return false;
0545   //    printf("Match lpz...\n");
0546 
0547   return true;
0548 }
0549 
0550 //
0551 // Store Tracker Rec Hits
0552 //
0553 void ResidualRefitting::StoreTrackerRecHits(DetId detid, const TrackerTopology* tTopo, int iTrack, int iRec) {
0554   int detector = -1;
0555   int subdetector = -1;
0556   int blade = -1;
0557   int disk = -1;
0558   int ladder = -1;
0559   int layer = -1;
0560   int module = -1;
0561   int panel = -1;
0562   int ring = -1;
0563   int side = -1;
0564   int wheel = -1;
0565 
0566   //Detector Info
0567 
0568   detector = detid.det();
0569   subdetector = detid.subdetId();
0570 
0571   if (detector != DetId::Tracker) {
0572     std::cout << "OMFG NOT THE TRACKER\n" << std::endl;
0573     return;
0574   }
0575 
0576   if (debug_)
0577     std::cout << "Tracker:: ";
0578   if (subdetector == ResidualRefitting::PXB) {
0579     layer = tTopo->pxbLayer(detid.rawId());
0580     ladder = tTopo->pxbLadder(detid.rawId());
0581     module = tTopo->pxbModule(detid.rawId());
0582     if (debug_)
0583       std::cout << "PXB"
0584                 << "\tlayer = " << layer << "\tladder = " << ladder << "\tmodule = " << module;
0585 
0586   } else if (subdetector == ResidualRefitting::PXF) {
0587     side = tTopo->pxfSide(detid.rawId());
0588     disk = tTopo->pxfDisk(detid.rawId());
0589     blade = tTopo->pxfBlade(detid.rawId());
0590     panel = tTopo->pxfPanel(detid.rawId());
0591     module = tTopo->pxfModule(detid.rawId());
0592     if (debug_)
0593       std::cout << "PXF"
0594                 << "\tside = " << side << "\tdisk = " << disk << "\tblade = " << blade << "\tpanel = " << panel
0595                 << "\tmodule = " << module;
0596 
0597   } else if (subdetector == ResidualRefitting::TIB) {
0598     layer = tTopo->tibLayer(detid.rawId());
0599     module = tTopo->tibModule(detid.rawId());
0600     if (debug_)
0601       std::cout << "TIB"
0602                 << "\tlayer = " << layer << "\tmodule = " << module;
0603   } else if (subdetector == ResidualRefitting::TID) {
0604     side = tTopo->tidSide(detid.rawId());
0605     wheel = tTopo->tidWheel(detid.rawId());
0606     ring = tTopo->tidRing(detid.rawId());
0607     if (debug_)
0608       std::cout << "TID"
0609                 << "\tside = " << side << "\twheel = " << wheel << "\tring = " << ring;
0610 
0611   } else if (subdetector == ResidualRefitting::TOB) {
0612     layer = tTopo->tobLayer(detid.rawId());
0613     module = tTopo->tobModule(detid.rawId());
0614     if (debug_)
0615       std::cout << "TOB"
0616                 << "\tlayer = " << layer << "\tmodule = " << module;
0617 
0618   } else if (subdetector == ResidualRefitting::TEC) {
0619     ring = tTopo->tecRing(detid.rawId());
0620     module = tTopo->tecModule(detid.rawId());
0621     if (debug_)
0622       std::cout << "TEC"
0623                 << "\tring = " << ring << "\tmodule = " << module;
0624   }
0625 
0626   //Do Storage
0627 
0628   storageTrackHit_.muonLink_[iRec] = iTrack;
0629   storageTrackHit_.detector_[iRec] = detector;
0630   storageTrackHit_.subdetector_[iRec] = subdetector;
0631   storageTrackHit_.blade_[iRec] = blade;
0632   storageTrackHit_.disk_[iRec] = disk;
0633   storageTrackHit_.ladder_[iRec] = ladder;
0634   storageTrackHit_.layer_[iRec] = layer;
0635   storageTrackHit_.module_[iRec] = module;
0636   storageTrackHit_.panel_[iRec] = panel;
0637   storageTrackHit_.ring_[iRec] = ring;
0638   storageTrackHit_.side_[iRec] = side;
0639   storageTrackHit_.wheel_[iRec] = wheel;
0640 }
0641 
0642 //
0643 // Store Muon info on P, Pt, eta, phi
0644 //
0645 void ResidualRefitting::muonInfo(ResidualRefitting::storage_muon& storeMuon, reco::TrackRef muon, int val) {
0646   storeMuon.pt_[val] = muon->pt();
0647   storeMuon.p_[val] = muon->p();
0648   storeMuon.eta_[val] = muon->eta();
0649   storeMuon.phi_[val] = muon->phi();
0650   storeMuon.charge_[val] = muon->charge();
0651   storeMuon.numRecHits_[val] = muon->numberOfValidHits();
0652   storeMuon.chiSq_[val] = muon->chi2();
0653   storeMuon.ndf_[val] = muon->ndof();
0654   storeMuon.chiSqOvrNdf_[val] = muon->normalizedChi2();
0655 }
0656 //
0657 // Fill a track extrapolation
0658 //
0659 void ResidualRefitting::trkExtrap(const DetId& detid,
0660                                   int iTrk,
0661                                   int iTrkLink,
0662                                   int iRec,
0663                                   const FreeTrajectoryState& freeTrajState,
0664                                   const LocalPoint& recPoint,
0665                                   storage_trackExtrap& storeTemp) {
0666   bool dump_ = debug_;
0667 
0668   if (dump_)
0669     std::cout << "In the trkExtrap function" << std::endl;
0670 
0671   float gpExtrapX = -99999;
0672   float gpExtrapY = -99999;
0673   float gpExtrapZ = -99999;
0674   float gpExtrapEta = -99999;
0675   float gpExtrapPhi = -99999;
0676 
0677   float lpX = -99999;
0678   float lpY = -99999;
0679   float lpZ = -99999;
0680 
0681   //
0682   // Get the local positions for the recHits
0683   //
0684 
0685   float recLpX = recPoint.x();
0686   float recLpY = recPoint.y();
0687   float recLpZ = recPoint.z();
0688 
0689   float resX = -9999;
0690   float resY = -9999;
0691   float resZ = -9999;
0692 
0693   const GeomDet* gdet = theService->trackingGeometry()->idToDet(detid);
0694 
0695   //    TrajectoryStateOnSurface surfTest =  prop.propagate(freeTrajState, gdet->surface());
0696   TrajectoryStateOnSurface surfTest = thePropagator->propagate(freeTrajState, gdet->surface());
0697 
0698   if (surfTest.isValid()) {
0699     GlobalPoint globTest = surfTest.globalPosition();
0700     gpExtrapX = globTest.x();
0701     gpExtrapY = globTest.y();
0702     gpExtrapZ = globTest.z();
0703     gpExtrapEta = globTest.eta();
0704     gpExtrapPhi = globTest.phi();
0705     LocalPoint loc = surfTest.localPosition();
0706     if (detid.det() == DetId::Muon || detid.det() == DetId::Tracker) {
0707       lpX = loc.x();
0708       lpY = loc.y();
0709       lpZ = loc.z();
0710 
0711       resX = lpX - recLpX;
0712       resY = lpY - recLpY;
0713       resZ = lpZ - recLpZ;
0714     }
0715   }
0716   storeTemp.muonLink_[iTrk] = iTrkLink;
0717   storeTemp.recLink_[iTrk] = iRec;
0718   storeTemp.gpX_[iTrk] = gpExtrapX;
0719   storeTemp.gpY_[iTrk] = gpExtrapY;
0720   storeTemp.gpZ_[iTrk] = gpExtrapZ;
0721   storeTemp.gpEta_[iTrk] = gpExtrapEta;
0722   storeTemp.gpPhi_[iTrk] = gpExtrapPhi;
0723   storeTemp.lpX_[iTrk] = lpX;
0724   storeTemp.lpY_[iTrk] = lpY;
0725   storeTemp.lpZ_[iTrk] = lpZ;
0726   storeTemp.resX_[iTrk] = resX;
0727   storeTemp.resY_[iTrk] = resY;
0728   storeTemp.resZ_[iTrk] = resZ;
0729 
0730   printf("station: %d\tsector: %d\tresX storage: %4.2f\n", ReturnStation(detid), ReturnSector(detid), resX);
0731 }
0732 //
0733 // Return the station
0734 //
0735 int ResidualRefitting::ReturnStation(DetId detid) {
0736   int station = -999;
0737 
0738   if (detid.det() == DetId::Muon) {
0739     int systemMuon = detid.subdetId();  // 1 DT; 2 CSC; 3 RPC
0740     if (systemMuon == MuonSubdetId::CSC) {
0741       CSCDetId id(detid.rawId());
0742       station = id.station();
0743 
0744     } else if (systemMuon == MuonSubdetId::DT) {
0745       DTWireId id(detid.rawId());
0746       station = id.station();
0747 
0748     } else if (systemMuon == MuonSubdetId::RPC) {
0749       RPCDetId id(detid.rawId());
0750       station = id.station();
0751     }
0752   }
0753 
0754   return station;
0755 }
0756 //
0757 // Return the sector
0758 //
0759 int ResidualRefitting::ReturnSector(DetId detid) {
0760   int sector = -999;
0761 
0762   if (detid.det() == DetId::Muon) {
0763     int systemMuon = detid.subdetId();  // 1 DT; 2 CSC; 3 RPC
0764     if (systemMuon == MuonSubdetId::DT) {
0765       DTWireId id(detid.rawId());
0766       sector = id.sector();
0767     }
0768   }
0769 
0770   return sector;
0771 }
0772 
0773 //
0774 // Store the SAM and Track position info at a particular rho
0775 //
0776 void ResidualRefitting::cylExtrapTrkSam(int recNum,
0777                                         reco::TrackRef track,
0778                                         ResidualRefitting::storage_trackExtrap& storage,
0779                                         double rho) {
0780   Cylinder::PositionType pos(0, 0, 0);
0781   Cylinder::RotationType rot;
0782 
0783   Cylinder::CylinderPointer myCylinder = Cylinder::build(pos, rot, rho);
0784   //    SteppingHelixPropagator inwardProp  ( theField, oppositeToMomentum );
0785   //    SteppingHelixPropagator outwardProp ( theField, alongMomentum );
0786   FreeTrajectoryState recoStart = freeTrajStateMuon(track);
0787   //    TrajectoryStateOnSurface recoProp = outwardProp.propagate(recoStart, *myCylinder);
0788   TrajectoryStateOnSurface recoProp = thePropagator->propagate(recoStart, *myCylinder);
0789 
0790   double xVal = -9999;
0791   double yVal = -9999;
0792   double zVal = -9999;
0793   double phiVal = -9999;
0794   double etaVal = -9999;
0795 
0796   if (recoProp.isValid()) {
0797     GlobalPoint recoPoint = recoProp.globalPosition();
0798     xVal = recoPoint.x();
0799     yVal = recoPoint.y();
0800     zVal = recoPoint.z();
0801     phiVal = recoPoint.phi();
0802     etaVal = recoPoint.eta();
0803   }
0804   storage.muonLink_[recNum] = recNum;
0805   storage.gpX_[recNum] = xVal;
0806   storage.gpY_[recNum] = yVal;
0807   storage.gpZ_[recNum] = zVal;
0808   storage.gpEta_[recNum] = etaVal;
0809   storage.gpPhi_[recNum] = phiVal;
0810 
0811   float rhoVal = sqrt(xVal * xVal + yVal * yVal);
0812 
0813   printf("Cylinder: rho = %4.2f\tphi = %4.2f\teta = %4.2f\n", rhoVal, phiVal, etaVal);
0814   if (debug_)
0815     printf("Cylinder: rho = %4.2f\tphi = %4.2f\teta = %4.2f\n", rhoVal, phiVal, etaVal);
0816 }
0817 ///////////////////////////////////////////////////////////////////////////////
0818 //Pre-Job junk
0819 ///////////////////////////////////////////////////////////////////////////////
0820 
0821 //
0822 // zero storage
0823 //
0824 void ResidualRefitting::zero_storage() {
0825   if (debug_)
0826     printf("zero_storage\n");
0827 
0828   zero_muon(&storageGmrOld_);
0829   zero_muon(&storageGmrNew_);
0830   zero_muon(&storageSamNew_);
0831   zero_muon(&storageTrkNew_);
0832   zero_muon(&storageGmrNoSt1_);
0833   zero_muon(&storageSamNoSt1_);
0834   zero_muon(&storageGmrNoSt2_);
0835   zero_muon(&storageSamNoSt2_);
0836   zero_muon(&storageGmrNoSt3_);
0837   zero_muon(&storageSamNoSt3_);
0838   zero_muon(&storageGmrNoSt4_);
0839   zero_muon(&storageSamNoSt4_);
0840   //zero out the tracker
0841   zero_muon(&storageGmrNoPXBLayer1);
0842   zero_muon(&storageGmrNoPXBLayer2);
0843   zero_muon(&storageGmrNoPXBLayer3);
0844 
0845   zero_muon(&storageGmrNoPXF);
0846 
0847   zero_muon(&storageGmrNoTIBLayer1);
0848   zero_muon(&storageGmrNoTIBLayer2);
0849   zero_muon(&storageGmrNoTIBLayer3);
0850   zero_muon(&storageGmrNoTIBLayer4);
0851 
0852   zero_muon(&storageGmrNoTID);
0853 
0854   zero_muon(&storageGmrNoTOBLayer1);
0855   zero_muon(&storageGmrNoTOBLayer2);
0856   zero_muon(&storageGmrNoTOBLayer3);
0857   zero_muon(&storageGmrNoTOBLayer4);
0858   zero_muon(&storageGmrNoTOBLayer5);
0859   zero_muon(&storageGmrNoTOBLayer6);
0860 
0861   zero_muon(&storageGmrNoTEC);
0862 
0863   zero_muon(&storageTrkNoPXBLayer1);
0864   zero_muon(&storageTrkNoPXBLayer2);
0865   zero_muon(&storageTrkNoPXBLayer3);
0866 
0867   zero_muon(&storageTrkNoPXF);
0868 
0869   zero_muon(&storageTrkNoTIBLayer1);
0870   zero_muon(&storageTrkNoTIBLayer2);
0871   zero_muon(&storageTrkNoTIBLayer3);
0872   zero_muon(&storageTrkNoTIBLayer4);
0873 
0874   zero_muon(&storageTrkNoTID);
0875 
0876   zero_muon(&storageTrkNoTOBLayer1);
0877   zero_muon(&storageTrkNoTOBLayer2);
0878   zero_muon(&storageTrkNoTOBLayer3);
0879   zero_muon(&storageTrkNoTOBLayer4);
0880   zero_muon(&storageTrkNoTOBLayer5);
0881   zero_muon(&storageTrkNoTOBLayer6);
0882 
0883   zero_muon(&storageTrkNoTEC);
0884 
0885   zero_trackExtrap(&storageTrackExtrapRec_);
0886   zero_trackExtrap(&storageTrackExtrapTracker_);
0887   zero_trackExtrap(&storageTrackExtrapRecNoSt1_);
0888   zero_trackExtrap(&storageTrackExtrapRecNoSt2_);
0889   zero_trackExtrap(&storageTrackExtrapRecNoSt3_);
0890   zero_trackExtrap(&storageTrackExtrapRecNoSt4_);
0891 
0892   zero_trackExtrap(&trackExtrap120_);
0893 
0894   zero_trackExtrap(&samExtrap120_);
0895 
0896   zero_trackExtrap(&storageTrackNoPXBLayer1);
0897   zero_trackExtrap(&storageTrackNoPXBLayer2);
0898   zero_trackExtrap(&storageTrackNoPXBLayer3);
0899 
0900   zero_trackExtrap(&storageTrackNoPXF);
0901 
0902   zero_trackExtrap(&storageTrackNoTIBLayer1);
0903   zero_trackExtrap(&storageTrackNoTIBLayer2);
0904   zero_trackExtrap(&storageTrackNoTIBLayer3);
0905   zero_trackExtrap(&storageTrackNoTIBLayer4);
0906 
0907   zero_trackExtrap(&storageTrackNoTOBLayer1);
0908   zero_trackExtrap(&storageTrackNoTOBLayer2);
0909   zero_trackExtrap(&storageTrackNoTOBLayer3);
0910   zero_trackExtrap(&storageTrackNoTOBLayer4);
0911   zero_trackExtrap(&storageTrackNoTOBLayer5);
0912   zero_trackExtrap(&storageTrackNoTOBLayer6);
0913 
0914   zero_trackExtrap(&storageTrackNoTEC);
0915 
0916   zero_trackExtrap(&storageTrackNoTID);
0917 
0918   storageRecMuon_.n_ = 0;
0919   storageTrackHit_.n_ = 0;
0920 }
0921 //
0922 // Zero out a muon reference
0923 //
0924 void ResidualRefitting::zero_muon(ResidualRefitting::storage_muon* str) {
0925   str->n_ = 0;
0926 
0927   for (int i = 0; i < ResidualRefitting::N_MAX_STORED; i++) {
0928     str->pt_[i] = -9999;
0929     str->eta_[i] = -9999;
0930     str->p_[i] = -9999;
0931     str->phi_[i] = -9999;
0932     str->numRecHits_[i] = -9999;
0933     str->chiSq_[i] = -9999;
0934     str->ndf_[i] = -9999;
0935     str->chiSqOvrNdf_[i] = -9999;
0936   }
0937 }
0938 //
0939 // Zero track extrapolation
0940 //
0941 void ResidualRefitting::zero_trackExtrap(ResidualRefitting::storage_trackExtrap* str) {
0942   str->n_ = 0;
0943   for (int i = 0; i < ResidualRefitting::N_MAX_STORED_HIT; i++) {
0944     str->muonLink_[i] = -9999;
0945     str->recLink_[i] = -9999;
0946     str->gpX_[i] = -9999;
0947     str->gpY_[i] = -9999;
0948     str->gpZ_[i] = -9999;
0949     str->gpEta_[i] = -9999;
0950     str->gpPhi_[i] = -9999;
0951     str->lpX_[i] = -9999;
0952     str->lpY_[i] = -9999;
0953     str->lpZ_[i] = -9999;
0954     str->resX_[i] = -9999;
0955     str->resY_[i] = -9999;
0956     str->resZ_[i] = -9999;
0957   }
0958 }
0959 //
0960 // Begin Job
0961 //
0962 void ResidualRefitting::beginJob() {
0963   std::cout << "Creating file " << outputFileName_.c_str() << std::endl;
0964 
0965   outputFile_ = new TFile(outputFileName_.c_str(), "RECREATE");
0966 
0967   outputTree_ = new TTree("outputTree", "outputTree");
0968 
0969   outputTree_->Branch("eventInfo",
0970                       &eventInfo_,
0971                       "evtNum_/I:"
0972                       "runNum_/I");
0973 
0974   ResidualRefitting::branchMuon(storageGmrOld_, "gmrOld");
0975   ResidualRefitting::branchMuon(storageGmrNew_, "gmrNew");
0976   ResidualRefitting::branchMuon(storageGmrNoSt1_, "gmrNoSt1");
0977   ResidualRefitting::branchMuon(storageGmrNoSt2_, "gmrNoSt2");
0978   ResidualRefitting::branchMuon(storageGmrNoSt3_, "gmrNoSt3");
0979   ResidualRefitting::branchMuon(storageGmrNoSt4_, "gmrNoSt4");
0980 
0981   ResidualRefitting::branchMuon(storageSamNew_, "samNew");
0982   ResidualRefitting::branchMuon(storageSamNoSt1_, "samNoSt1");
0983   ResidualRefitting::branchMuon(storageSamNoSt2_, "samNoSt2");
0984   ResidualRefitting::branchMuon(storageSamNoSt3_, "samNoSt3");
0985   ResidualRefitting::branchMuon(storageSamNoSt4_, "samNoSt4");
0986 
0987   ResidualRefitting::branchMuon(storageTrkNew_, "trkNew");
0988   ResidualRefitting::branchMuon(storageGmrNoPXBLayer1, "gmrNoPXBLayer1");
0989   ResidualRefitting::branchMuon(storageGmrNoPXBLayer2, "gmrNoPXBLayer2");
0990   ResidualRefitting::branchMuon(storageGmrNoPXBLayer3, "gmrNoPXBLayer3");
0991   ResidualRefitting::branchMuon(storageGmrNoPXF, "gmrNoPXF");
0992   ResidualRefitting::branchMuon(storageGmrNoTIBLayer1, "gmrNoTIBLayer1");
0993   ResidualRefitting::branchMuon(storageGmrNoTIBLayer2, "gmrNoTIBLayer2");
0994   ResidualRefitting::branchMuon(storageGmrNoTIBLayer3, "gmrNoTIBLayer3");
0995   ResidualRefitting::branchMuon(storageGmrNoTIBLayer4, "gmrNoTIBLayer4");
0996   ResidualRefitting::branchMuon(storageGmrNoTID, "gmrNoTID");
0997   ResidualRefitting::branchMuon(storageGmrNoTOBLayer1, "gmrNoTOBLayer1");
0998   ResidualRefitting::branchMuon(storageGmrNoTOBLayer2, "gmrNoTOBLayer2");
0999   ResidualRefitting::branchMuon(storageGmrNoTOBLayer3, "gmrNoTOBLayer3");
1000   ResidualRefitting::branchMuon(storageGmrNoTOBLayer4, "gmrNoTOBLayer4");
1001   ResidualRefitting::branchMuon(storageGmrNoTOBLayer5, "gmrNoTOBLayer5");
1002   ResidualRefitting::branchMuon(storageGmrNoTOBLayer6, "gmrNoTOBLayer6");
1003   ResidualRefitting::branchMuon(storageGmrNoTEC, "gmrNoTEC");
1004 
1005   ResidualRefitting::branchMuon(storageTrkNoPXBLayer1, "trkNoPXBLayer1");
1006   ResidualRefitting::branchMuon(storageTrkNoPXBLayer2, "trkNoPXBLayer2");
1007   ResidualRefitting::branchMuon(storageTrkNoPXBLayer3, "trkNoPXBLayer3");
1008   ResidualRefitting::branchMuon(storageTrkNoPXF, "trkNoPXF");
1009   ResidualRefitting::branchMuon(storageTrkNoTIBLayer1, "trkNoTIBLayer1");
1010   ResidualRefitting::branchMuon(storageTrkNoTIBLayer2, "trkNoTIBLayer2");
1011   ResidualRefitting::branchMuon(storageTrkNoTIBLayer3, "trkNoTIBLayer3");
1012   ResidualRefitting::branchMuon(storageTrkNoTIBLayer4, "trkNoTIBLayer4");
1013   ResidualRefitting::branchMuon(storageTrkNoTID, "trkNoTID");
1014   ResidualRefitting::branchMuon(storageTrkNoTOBLayer1, "trkNoTOBLayer1");
1015   ResidualRefitting::branchMuon(storageTrkNoTOBLayer2, "trkNoTOBLayer2");
1016   ResidualRefitting::branchMuon(storageTrkNoTOBLayer3, "trkNoTOBLayer3");
1017   ResidualRefitting::branchMuon(storageTrkNoTOBLayer4, "trkNoTOBLayer4");
1018   ResidualRefitting::branchMuon(storageTrkNoTOBLayer5, "trkNoTOBLayer5");
1019   ResidualRefitting::branchMuon(storageTrkNoTOBLayer6, "trkNoTOBLayer6");
1020   ResidualRefitting::branchMuon(storageTrkNoTEC, "trkNoTEC");
1021 
1022   outputBranch_ = outputTree_->Branch("recHitsNew",
1023                                       &storageRecMuon_,
1024 
1025                                       "n_/I:"
1026                                       "muonLink_[1000]/I:"
1027 
1028                                       "system_[1000]/I:"
1029                                       "endcap_[1000]/I:"
1030                                       "station_[1000]/I:"
1031                                       "ring_[1000]/I:"
1032                                       "chamber_[1000]/I:"
1033                                       "layer_[1000]/I:"
1034                                       "superLayer_[1000]/I:"
1035                                       "wheel_[1000]/I:"
1036                                       "sector_[1000]/I:"
1037 
1038                                       "gpX_[1000]/F:"
1039                                       "gpY_[1000]/F:"
1040                                       "gpZ_[1000]/F:"
1041                                       "gpEta_[1000]/F:"
1042                                       "gpPhi_[1000]/F:"
1043                                       "lpX_[1000]/F:"
1044                                       "lpY_[1000]/F:"
1045                                       "lpZ_[1000]/F");
1046 
1047   outputBranch_ = outputTree_->Branch("recHitsTracker",
1048                                       &storageTrackHit_,
1049 
1050                                       "n_/I:"
1051 
1052                                       "muonLink_[1000]/I:"
1053                                       "detector_[1000]/I:"
1054                                       "subdetector_[1000]/I:"
1055                                       "blade_[1000]/I:"
1056                                       "disk_[1000]/I:"
1057                                       "ladder_[1000]/I:"
1058                                       "layer_[1000]/I:"
1059                                       "module_[1000]/I:"
1060                                       "panel_[1000]/I:"
1061                                       "ring_[1000]/I:"
1062                                       "side_[1000]/I:"
1063                                       "wheel_[1000]/I:"
1064 
1065                                       "gpX_[1000]/F:"
1066                                       "gpY_[1000]/F:"
1067                                       "gpZ_[1000]/F:"
1068                                       "gpEta_[1000]/F:"
1069                                       "gpPhi_[1000]/F:"
1070                                       "lpX_[1000]/F:"
1071                                       "lpY_[1000]/F:"
1072                                       "lpZ_[1000]/F");
1073 
1074   ResidualRefitting::branchTrackExtrap(storageTrackExtrapRec_, "trkExtrap");
1075   ResidualRefitting::branchTrackExtrap(storageTrackExtrapRecNoSt1_, "trkExtrapNoSt1");
1076   ResidualRefitting::branchTrackExtrap(storageTrackExtrapRecNoSt2_, "trkExtrapNoSt2");
1077   ResidualRefitting::branchTrackExtrap(storageTrackExtrapRecNoSt3_, "trkExtrapNoSt3");
1078   ResidualRefitting::branchTrackExtrap(storageTrackExtrapRecNoSt4_, "trkExtrapNoSt4");
1079 
1080   ResidualRefitting::branchTrackExtrap(storageTrackExtrapTracker_, "trkExtrapTracker");
1081   ResidualRefitting::branchTrackExtrap(storageTrackNoPXF, "trkExtrapNoPXF");
1082   ResidualRefitting::branchTrackExtrap(storageTrackNoPXBLayer1, "trkExtrapNoPXBLayer1");
1083   ResidualRefitting::branchTrackExtrap(storageTrackNoPXBLayer2, "trkExtrapNoPXBLayer2");
1084   ResidualRefitting::branchTrackExtrap(storageTrackNoPXBLayer3, "trkExtrapNoPXBLayer3");
1085   ResidualRefitting::branchTrackExtrap(storageTrackNoTIBLayer1, "trkExtrapNoTIBLayer1");
1086   ResidualRefitting::branchTrackExtrap(storageTrackNoTIBLayer2, "trkExtrapNoTIBLayer2");
1087   ResidualRefitting::branchTrackExtrap(storageTrackNoTIBLayer3, "trkExtrapNoTIBLayer3");
1088   ResidualRefitting::branchTrackExtrap(storageTrackNoTIBLayer4, "trkExtrapNoTIBLayer4");
1089   ResidualRefitting::branchTrackExtrap(storageTrackNoTID, "trkExtrapNoTID");
1090   ResidualRefitting::branchTrackExtrap(storageTrackNoTOBLayer1, "trkExtrapNoTOBLayer1");
1091   ResidualRefitting::branchTrackExtrap(storageTrackNoTOBLayer2, "trkExtrapNoTOBLayer2");
1092   ResidualRefitting::branchTrackExtrap(storageTrackNoTOBLayer3, "trkExtrapNoTOBLayer3");
1093   ResidualRefitting::branchTrackExtrap(storageTrackNoTOBLayer4, "trkExtrapNoTOBLayer4");
1094   ResidualRefitting::branchTrackExtrap(storageTrackNoTOBLayer5, "trkExtrapNoTOBLayer5");
1095   ResidualRefitting::branchTrackExtrap(storageTrackNoTOBLayer6, "trkExtrapNoTOBLayer6");
1096   ResidualRefitting::branchTrackExtrap(storageTrackNoTEC, "trkExtrapNoTEC");
1097 
1098   ResidualRefitting::branchTrackExtrap(trackExtrap120_, "trackCyl120");
1099   ResidualRefitting::branchTrackExtrap(samExtrap120_, "samCyl120");
1100 }
1101 //
1102 // Set the Muon Branches
1103 //
1104 void ResidualRefitting::branchMuon(ResidualRefitting::storage_muon& storageTmp, std::string branchName) {
1105   outputBranch_ = outputTree_->Branch(branchName.c_str(),
1106                                       &storageTmp,
1107                                       "n_/I:"
1108                                       "charge_[10]/I:"
1109                                       "pt_[10]/F:"
1110                                       "eta_[10]/F:"
1111                                       "p_[10]/F:"
1112                                       "phi_[10]/F:"
1113                                       "numRecHits_[10]/I:"
1114                                       "chiSq_[10]/F:"
1115                                       "ndf_[10]/F:"
1116                                       "chiSqOvrNdf_[10]/F"
1117 
1118   );
1119 }
1120 //
1121 // Set the Branches for Track Extrapolations
1122 //
1123 void ResidualRefitting::branchTrackExtrap(ResidualRefitting::storage_trackExtrap& storageTmp, std::string branchName) {
1124   outputBranch_ = outputTree_->Branch(branchName.c_str(),
1125                                       &storageTmp,
1126                                       "n_/I:"
1127                                       "muonLink_[1000]/I:"
1128                                       "recLink_[1000]/I:"
1129                                       "gpX_[1000]/F:"
1130                                       "gpY_[1000]/F:"
1131                                       "gpZ_[1000]/F:"
1132                                       "gpEta_[1000]/F:"
1133                                       "gpPhi_[1000]/F:"
1134                                       "lpX_[1000]/F:"
1135                                       "lpY_[1000]/F:"
1136                                       "lpZ_[1000]/F:"
1137                                       "resX_[1000]/F:"
1138                                       "resY_[1000]/F:"
1139                                       "resZ_[1000]/F"
1140 
1141   );
1142 }
1143 //
1144 // End Job
1145 //
1146 void ResidualRefitting::endJob() {
1147   outputFile_->Write();
1148 
1149   outputFile_->Close();
1150 }
1151 //
1152 // Return a Free Trajectory state for a muon track
1153 //
1154 FreeTrajectoryState ResidualRefitting::freeTrajStateMuon(reco::TrackRef muon) {
1155   math::XYZPoint innerPos = muon->referencePoint();
1156   math::XYZVector innerMom = muon->momentum();
1157   if (debug_)
1158     std::cout << "Inner Pos: "
1159               << "\tx = " << innerPos.X() << "\ty = " << innerPos.Y() << "\tz = " << innerPos.Z() << std::endl;
1160 
1161   GlobalPoint innerPoint(innerPos.X(), innerPos.Y(), innerPos.Z());
1162   GlobalVector innerVec(innerMom.X(), innerMom.Y(), innerMom.Z());
1163 
1164   FreeTrajectoryState recoStart(innerPoint, innerVec, muon->charge(), theField);
1165   return recoStart;
1166 }
1167 
1168 /////////////////////////////////////////////////////////////////////////////////
1169 // nTuple value Dumps
1170 /////////////////////////////////////////////////////////////////////////////////
1171 
1172 //
1173 // dump Track Extrapolation
1174 //
1175 void ResidualRefitting::dumpTrackExtrap(const ResidualRefitting::storage_trackExtrap& track) {
1176   std::cout << "\n\nExtrapolation Dump:\n";
1177   for (unsigned int i = 0; i < (unsigned int)track.n_; i++) {
1178     //      double rho = sqrt( (float)track.gpX_[i] * (float)track.gpX_[i] + (float)track.gpY_[i] * (float)track.gpY_[i]  );
1179 
1180     printf("%d\tmuonLink= %d", i, (int)track.muonLink_[i]);
1181     printf("\trecLink = %d", (int)track.recLink_[i]);
1182     //      printf ("\tGlobal\tx = %0.3f"       ,       (float)track.gpX_[i]    );
1183     //      printf ("\ty = %0.3f"       ,       (float)track.gpY_[i]    );
1184     //      printf ("\tz = %0.3f"       ,       (float)track.gpZ_[i]    );
1185     //      printf ("\trho =%0.3f"      ,       rho                     );
1186     //      printf ("\teta = %0.3f"     ,       (float)track.gpEta_[i]  );
1187     //      printf ("\tphi = %0.3f"     ,       (float)track.gpPhi_[i]  );
1188     printf("\t\tLocal\tx = %0.3f", (float)track.lpX_[i]);
1189     printf("\ty = %0.3f", (float)track.lpY_[i]);
1190     printf("\tz = %0.3f\n", (float)track.lpZ_[i]);
1191   }
1192 }
1193 //
1194 // dump Muon Rec Hits
1195 //
1196 void ResidualRefitting::dumpMuonRecHits(const ResidualRefitting::storage_hit& hit) {
1197   std::cout << "Muon Rec Hits Dump:\n";
1198   for (unsigned int i = 0; i < (unsigned int)hit.n_; i++) {
1199     //      double rho = sqrt( (float)hit.gpX_[i] * (float)hit.gpX_[i] + (float)hit.gpY_[i] * (float)hit.gpY_[i]  );
1200 
1201     printf("%d\tsubdetector = %d\t superLayer =%d", i, (int)hit.system_[i], (int)hit.superLayer_[i]);
1202     //      printf ("\tGlobal\tx = %0.3f"           ,       (float)hit.gpX_[i]          );
1203     //      printf ("\ty = %0.3f"               ,       (float)hit.gpY_[i]          );
1204     //      printf ("\tz = %0.3f"               ,       (float)hit.gpZ_[i]          );
1205     //      printf ("\trho =%0.3f"              ,       rho                         );
1206     //      printf ("\teta = %0.3f"             ,       (float)hit.gpEta_[i]        );
1207     //      printf ("\tphi = %0.3f\n"           ,       (float)hit.gpPhi_[i]        );
1208     printf("\t\tLocal\tx = %0.3f", (float)hit.lpX_[i]);
1209     printf("\ty = %0.3f", (float)hit.lpY_[i]);
1210     printf("\tz = %0.3f\n", (float)hit.lpZ_[i]);
1211   }
1212 }
1213 //
1214 // dump Tracker Rec Hits
1215 //
1216 void ResidualRefitting::dumpTrackHits(const ResidualRefitting::storage_trackHit& hit) {
1217   std::cout << "Tracker Rec Hits Dump:\n";
1218   for (unsigned int i = 0; i < (unsigned int)hit.n_; i++) {
1219     //      double rho = sqrt( (float)hit.gpX_[i] * (float)hit.gpX_[i] + (float)hit.gpY_[i] * (float)hit.gpY_[i]  );
1220 
1221     printf("%d\tsubdetector = %d", i, (int)hit.subdetector_[i]);
1222     printf("\tlayer = %d", (int)hit.layer_[i]);
1223     //      printf ("\tGlobal\tx = %0.3f"           ,       (float)hit.gpX_[i]          );
1224     //      printf ("\ty = %0.3f"               ,       (float)hit.gpY_[i]          );
1225     //      printf ("\tz = %0.3f"               ,       (float)hit.gpZ_[i]          );
1226     //      printf ("\trho =%0.3f"              ,       rho                         );
1227     //      printf ("\teta = %0.3f"             ,       (float)hit.gpEta_[i]        );
1228     //      printf ("\tphi = %0.3f\n"           ,       (float)hit.gpPhi_[i]        );
1229     printf("\t\tLocal\tx = %0.3f", (float)hit.lpX_[i]);
1230     printf("\ty = %0.3f", (float)hit.lpY_[i]);
1231     printf("\tz = %0.3f\n", (float)hit.lpZ_[i]);
1232   }
1233 }
1234 //
1235 //Dump a TrackRef
1236 //
1237 void ResidualRefitting::dumpTrackRef(reco::TrackRef muon, std::string str) {
1238   float pt = muon->pt();
1239   float p = muon->p();
1240   float eta = muon->eta();
1241   float phi = muon->phi();
1242   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);
1243 }
1244 
1245 DEFINE_FWK_MODULE(ResidualRefitting);
1246 
1247 ////////////////////////////////////////////////////////////////////////////////////////////
1248 //Deprecated
1249 ////////////////////////////////////////////////////////////////////////////////////////////