Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:57:07 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!