Reference Manual

Go to the documentation of this file.
00001 #include <iostream>
00002 #include "ResidualRefitting.h"
00003 #include <iomanip>
00004 
00005 //framework includes
00006 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
00007 #include "DataFormats/MuonDetId/interface/DTWireId.h"
00008 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
00009 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
00010 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
00011 #include "Geometry/Records/interface/IdealGeometryRecord.h"
00012 
00013 #include "FWCore/Framework/interface/Event.h"
00014 #include "FWCore/Framework/interface/EventSetup.h"
00015 #include "FWCore/Framework/interface/ESHandle.h"
00016 #include "FWCore/Utilities/interface/EDMException.h"
00017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00018 #include "FWCore/Framework/interface/MakerMacros.h"
00019 
00020 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
00021 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
00022 #include "Geometry/CommonDetUnit/interface/GeomDetEnumerators.h"
00023 
00024 #include "RecoMuon/DetLayers/interface/MuonDetLayerGeometry.h"
00025 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBuilder.h"
00026 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
00027 
00028 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
00029 
00030 #include "TrackingTools/DetLayers/interface/DetLayer.h"
00031 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00032 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00033 
00034 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00035 
00036 
00037 
00038 ResidualRefitting::ResidualRefitting( const edm::ParameterSet & cfg ) :
00039   outputFileName_               ( cfg.getUntrackedParameter<std::string>("histoutputFile") ),
00040   PropagatorSource_             ( cfg.getParameter<std::string>("propagator")), 
00041   muons_                        ( cfg.getParameter<edm::InputTag>( "muons"              ) ),
00042   muonsRemake_                  ( cfg.getParameter<edm::InputTag>("muonsRemake" ) ),                    //This Feels Misalignment
00043   muonsNoStation1_              ( cfg.getParameter<edm::InputTag>("muonsNoStation1") ),
00044   muonsNoStation2_              ( cfg.getParameter<edm::InputTag>("muonsNoStation2") ),
00045   muonsNoStation3_              ( cfg.getParameter<edm::InputTag>("muonsNoStation3") ),
00046   muonsNoStation4_              ( cfg.getParameter<edm::InputTag>("muonsNoStation4") ),
00047 
00048 /*
00049   muonsNoPXBLayer1_             ( cfg.getParameter<edm::InputTag>("muonsNoPXBLayer1"    ) ),
00050   muonsNoPXBLayer2_             ( cfg.getParameter<edm::InputTag>("muonsNoPXBLayer1"    ) ),
00051   muonsNoPXBLayer3_             ( cfg.getParameter<edm::InputTag>("muonsNoPXBLayer1"    ) ),
00052 
00053   muonsNoTIBLayer1_                     ( cfg.getParameter<edm::InputTag>("muonsNoTIBLayer1"    ) ),
00054   muonsNoTIBLayer2_                     ( cfg.getParameter<edm::InputTag>("muonsNoTIBLayer2"    ) ),
00055   muonsNoTIBLayer3_                     ( cfg.getParameter<edm::InputTag>("muonsNoTIBLayer3"    ) ),
00056   muonsNoTIBLayer4_                     ( cfg.getParameter<edm::InputTag>("muonsNoTIBLayer4"    ) ),
00057 
00058   muonsNoTOBLayer1_                     ( cfg.getParameter<edm::InputTag>("muonsNoTOBLayer1"    ) ),
00059   muonsNoTOBLayer2_                     ( cfg.getParameter<edm::InputTag>("muonsNoTOBLayer2"    ) ),
00060   muonsNoTOBLayer3_                     ( cfg.getParameter<edm::InputTag>("muonsNoTOBLayer3"    ) ),
00061   muonsNoTOBLayer4_                     ( cfg.getParameter<edm::InputTag>("muonsNoTOBLayer4"    ) ),
00062   muonsNoTOBLayer5_                     ( cfg.getParameter<edm::InputTag>("muonsNoTOBLayer5"    ) ),
00063   muonsNoTOBLayer6_                     ( cfg.getParameter<edm::InputTag>("muonsNoTOBLayer6"    ) ),*/
00064   debug_                                ( cfg.getUntrackedParameter<bool>("doDebug"     ) ),
00065   outputFile_(0),
00066   outputTree_(0),
00067   outputBranch_(0),
00068   theField(0)
00069 {
00070         eventInfo_.evtNum_ = 0;
00071         eventInfo_.evtNum_ = 0;
00072 
00073 // service parameters
00074         edm::ParameterSet serviceParameters = cfg.getParameter<edm::ParameterSet>("ServiceParameters");
00075   
00076 // the services
00077         theService = new MuonServiceProxy(serviceParameters);
00078         
00079 }  //The constructor
00080 
00081 void ResidualRefitting::analyze(const edm::Event& event, const edm::EventSetup& eventSetup) {
00082         
00083         if(debug_) printf("STARTING EVENT\n");
00084 
00085         eventInfo_.evtNum_ = (int)event.id().run();
00086         eventInfo_.runNum_ = (int)event.id().event();
00087 
00088         // Generator Collection
00089 
00090 // The original muon collection that is sitting in memory
00091         edm::Handle<reco::MuonCollection> muons;
00092 
00093         edm::Handle<reco::TrackCollection> muonTracks;
00094         edm::Handle<reco::TrackCollection> muonsNoSt1;
00095         edm::Handle<reco::TrackCollection> muonsNoSt2;
00096         edm::Handle<reco::TrackCollection> muonsNoSt3;
00097         edm::Handle<reco::TrackCollection> muonsNoSt4;
00098 
00099 
00100         event.getByLabel(muons_                         , muons ); //set label to muons
00101         event.getByLabel(muonsRemake_           , muonTracks);  
00102         event.getByLabel(muonsNoStation1_       , muonsNoSt1);
00103         event.getByLabel(muonsNoStation2_       , muonsNoSt2);
00104         event.getByLabel(muonsNoStation3_       , muonsNoSt3);
00105         event.getByLabel(muonsNoStation4_       , muonsNoSt4);
00106         
00107         
00108 /*
00109 //      std::cout<<"Muon Collection No PXB "<<std::endl;
00110 //Tracker Barrel Pixel Refits
00111         edm::Handle<MuonCollection> muonsNoPXBLayer1Coll;
00112         event.getByLabel(muonsNoPXBLayer1_, muonsNoPXBLayer1Coll);
00113         edm::Handle<MuonCollection> muonsNoPXBLayer2Coll;
00114         event.getByLabel(muonsNoPXBLayer2_, muonsNoPXBLayer2Coll);
00115         edm::Handle<MuonCollection> muonsNoPXBLayer3Coll;
00116         event.getByLabel(muonsNoPXBLayer3_, muonsNoPXBLayer3Coll);
00117 //      std::cout<<"Muon Collection No TIB "<<std::endl;
00118 // Tracker Inner Barrel Refits
00119         edm::Handle<MuonCollection> muonsNoTIBLayer1Coll;
00120         event.getByLabel(muonsNoTIBLayer1_, muonsNoTIBLayer1Coll);
00121         edm::Handle<MuonCollection> muonsNoTIBLayer2Coll;
00122         event.getByLabel(muonsNoTIBLayer2_, muonsNoTIBLayer2Coll);
00123         edm::Handle<MuonCollection> muonsNoTIBLayer3Coll;
00124         event.getByLabel(muonsNoTIBLayer3_, muonsNoTIBLayer3Coll);
00125         edm::Handle<MuonCollection> muonsNoTIBLayer4Coll;
00126         event.getByLabel(muonsNoTIBLayer4_, muonsNoTIBLayer4Coll);
00127 
00128 //      std::cout<<"Muon Collection No TOB "<<std::endl;
00129 
00130 //Tracker outer barrel refits
00131         edm::Handle<MuonCollection> muonsNoTOBLayer1Coll;
00132         event.getByLabel(muonsNoTOBLayer1_, muonsNoTOBLayer1Coll);
00133         edm::Handle<MuonCollection> muonsNoTOBLayer2Coll;
00134         event.getByLabel(muonsNoTOBLayer2_, muonsNoTOBLayer2Coll);
00135         edm::Handle<MuonCollection> muonsNoTOBLayer3Coll;
00136         event.getByLabel(muonsNoTOBLayer3_, muonsNoTOBLayer3Coll);
00137         edm::Handle<MuonCollection> muonsNoTOBLayer4Coll;
00138         event.getByLabel(muonsNoTOBLayer4_, muonsNoTOBLayer4Coll);
00139         edm::Handle<MuonCollection> muonsNoTOBLayer5Coll;
00140         event.getByLabel(muonsNoTOBLayer5_, muonsNoTOBLayer5Coll);
00141         edm::Handle<MuonCollection> muonsNoTOBLayer6Coll;
00142         event.getByLabel(muonsNoTOBLayer6_, muonsNoTOBLayer6Coll);
00143 */
00144         //magnetic field information    
00145         edm::ESHandle<MagneticField> field;
00146         edm::ESHandle<GlobalTrackingGeometry> globalTrackingGeometry;
00147         eventSetup.get<IdealMagneticFieldRecord>().get(field);
00148         eventSetup.get<GlobalTrackingGeometryRecord>().get(globalTrackingGeometry);
00149         eventSetup.get<TrackingComponentsRecord>().get( PropagatorSource_, thePropagator );
00150         theField = &*field;
00151         
00152         
00153         theService->update(eventSetup);
00154 
00155 //Zero storage
00156     zero_storage();
00157 
00158   
00159 //Do the Gmr Muons from the unModified Collection
00160 
00161 /*
00162         int iGmr = 0;
00163         if ( (muons->end() - muons->begin()) > 0) printf("Data Dump:: Original GMR Muons\n");
00164         for ( MuonCollection::const_iterator muon = muons->begin(); muon!=muons->end(); muon++, iGmr++) {
00165                 if ( iGmr >= ResidualRefitting::N_MAX_STORED) break; // error checking
00166                 if (!debug
00167                 
00168                 dumpTrackRef(muon->combinedMuon(), "cmb"); 
00169                 dumpTrackRef(muon->standAloneMuon(), "sam");
00170                 dumpTrackRef(muon->track(), "trk");
00171                 
00172 
00173         }
00174         storageGmrOld_.n_ = iGmr;
00175         storageSamNew_.n_ = iGmr;
00176 */
00177 
00178 //Refitted muons
00179         if (debug_) printf("Data Dump:: Rebuilt GMR Muon Track With TeV refitter default\n");
00180         int iGmrRemake = 0;
00181         for ( reco::TrackCollection::const_iterator muon = muonTracks->begin(); muon!=muonTracks->end(); muon++, iGmrRemake++) {
00182                 if ( iGmrRemake >= ResidualRefitting::N_MAX_STORED) break; // error checking
00183                         // from TrackInfoProducer/test/TrackInfoAnalyzerExample.cc
00184                 reco::TrackRef trackref=reco::TrackRef(muonTracks,iGmrRemake);
00185 
00186                         if (debug_) dumpTrackRef(trackref, "gmr");
00187                         muonInfo(storageGmrNew_,trackref,iGmrRemake);
00188 
00189         }
00190         storageGmrNew_.n_ = iGmrRemake;
00191                 
00192 
00193 
00194         if (debug_) printf("muons Remake");
00195         if (debug_) printf("-----------------------------------------\n");
00196         CollectTrackHits(muonTracks, storageTrackExtrapRec_, eventSetup);
00197 
00198 
00199         if (true) {
00200                 printf("muons No Station 1");
00201                 printf("-----------------------------------------\n");
00202         }
00203         NewTrackMeasurements(muonTracks, muonsNoSt1, storageTrackExtrapRecNoSt1_);
00204 
00205         if (true) {
00206                 printf("muons No Station 2");
00207                 printf("-----------------------------------------\n");
00208         }
00209         NewTrackMeasurements(muonTracks, muonsNoSt2, storageTrackExtrapRecNoSt2_);
00210 
00211         if (true) {
00212                 printf("muons No Station 3");
00213                 printf("-----------------------------------------\n");
00214         }
00215         NewTrackMeasurements(muonTracks, muonsNoSt3, storageTrackExtrapRecNoSt3_);
00216 
00217         if (true) {
00218                 printf("muons No Station 4");
00219                 printf("-----------------------------------------\n");
00220         }
00221         NewTrackMeasurements(muonTracks, muonsNoSt4, storageTrackExtrapRecNoSt4_);
00222 
00223 
00224 //      dumpMuonRecHits(storageRecMuon_); 
00225                 
00226 /****************************************************************************************************************************************/
00227 
00228   
00229 /*
00230  *      extrapolates track to a cylinder.
00231  *  commented for cosmic runs with no tracker in reco muons!!
00232  *
00233 */
00234 
00235 
00236         int iGmrCyl = 0;
00237         for (reco::MuonCollection::const_iterator muon = muons->begin(); muon != muons->end(); muon++, iGmrCyl++) {
00238 
00239                 dumpTrackRef(muon->combinedMuon(), "cmb"); 
00240                 dumpTrackRef(muon->standAloneMuon(), "sam");
00241                 dumpTrackRef(muon->track(), "trk");
00242 
00243                 cylExtrapTrkSam(iGmrCyl, muon->standAloneMuon() , samExtrap120_         , 120.);
00244                 cylExtrapTrkSam(iGmrCyl, muon->track()                  , trackExtrap120_       , 120.);
00245 
00246         }
00247         samExtrap120_.n_         = iGmrCyl;
00248         trackExtrap120_.n_       = iGmrCyl; 
00249 
00250 
00251         if (iGmrRemake > 0 || iGmrCyl > 0) {
00252                 outputTree_ -> Fill();
00253                 std::cout << "FILLING NTUPLE!" << std::endl;
00254                 std::cout << "Entries Recorded: " << outputTree_ -> GetEntries() << " Branch :: " << outputBranch_ -> GetEntries() <<  std::endl<<std::endl;
00255         } else std::cout<<"no tracks -- no fill!\n"<<std::endl<<std::endl;
00256 
00257 //  /*************************************************************************************************************/
00258 //  //END OF ntuple dumper
00259 //  //END OF ntuple dumper
00260 //  /***********************************************************************************************************/
00261 }
00262 //end Analyze() main function
00263 
00264 //------------------------------------------------------------------------------
00265 //  
00266 // Destructor
00267 // 
00268 ResidualRefitting::~ResidualRefitting() {
00269   delete outputFile_;
00270   delete theService;
00271 }
00272 //
00273 // Track Collection Analysis
00274 //
00275 void ResidualRefitting::CollectTrackHits(edm::Handle<reco::TrackCollection> trackColl,
00276                                          ResidualRefitting::storage_trackExtrap& trackExtrap,
00277                                          const edm::EventSetup& eventSetup) {
00278 
00279         //Retrieve tracker topology from geometry
00280         edm::ESHandle<TrackerTopology> tTopoHandle;
00281         eventSetup.get<IdealGeometryRecord>().get(tTopoHandle);
00282         const TrackerTopology* const tTopo = tTopoHandle.product();
00283 
00284         int iMuonHit = 0;
00285         int iTrackHit= 0;
00286         int numTracks= 0;       
00287 
00288         for ( reco::TrackCollection::const_iterator muon = trackColl->begin(); muon!=trackColl->end(); muon++) {
00289 
00290                 int iTrack = muon - trackColl->begin();
00291             reco::TrackRef trackref=reco::TrackRef(trackColl,iTrack);
00292                 FreeTrajectoryState recoStart = ResidualRefitting::freeTrajStateMuon(trackref);
00293 
00294                 if (debug_) dumpTrackRef(trackref, "CollectTrackHits Track");
00295         
00296                 for (trackingRecHit_iterator rec = muon->recHitsBegin(); rec != muon->recHitsEnd(); rec++) {
00297 
00298                         int iRec = rec - muon->recHitsBegin();  
00299                         DetId detid = (*rec)->geographicalId(); 
00300 
00301                         if (detid.det() != DetId::Muon && detid.det() != DetId::Tracker) {
00302                                 if (debug_) printf("Rec Hit not from muon system or tracker... continuing...\n");
00303                                 continue;
00304                         }
00305 //                      numTracks++;
00306 // Get Local and Global Position of Hits        
00307 
00308                         LocalPoint lp = (*rec)->localPosition();
00309                         float lpX       = lp.x();
00310                         float lpY       = lp.y();
00311                         float lpZ       = lp.z();
00312 
00313                         MuonTransientTrackingRecHit::MuonRecHitPointer mrhp =   
00314                                 MuonTransientTrackingRecHit::specificBuild(theService->trackingGeometry()->idToDet((**rec).geographicalId())
00315                                 ,&(**rec)); 
00316 
00317                         GlobalPoint gp = mrhp->globalPosition();
00318                         float gpRecX    = gp.x();
00319                         float gpRecY    = gp.y();
00320                         float gpRecZ    = gp.z();
00321                         float gpRecEta  = gp.eta();
00322                         float gpRecPhi  = gp.phi();
00323 
00324                         if (detid.det() == DetId::Muon) {
00325                         
00326                                 int systemMuon  = detid.subdetId(); // 1 DT; 2 CSC; 3 RPC
00327                                 int endcap              = -999;
00328                                 int station             = -999;
00329                                 int ring                = -999;
00330                                 int chamber             = -999;
00331                                 int layer               = -999;
00332                                 int superLayer  = -999;
00333                                 int wheel               = -999;
00334                                 int sector              = -999;
00335                                 if ( systemMuon == MuonSubdetId::CSC) {
00336                                         CSCDetId id(detid.rawId());
00337                                         endcap          = id.endcap();
00338                                         station         = id.station();
00339                                         ring            = id.ring();
00340                                         chamber         = id.chamber();
00341                                         layer           = id.layer();
00342                                         if (debug_)printf("CSC\t[endcap][station][ringN][chamber][layer]:[%d][%d][%d][%d][%d]\t",
00343                                                  endcap, station, ring, chamber, layer);
00344 
00345                                 }
00346                                 else if ( systemMuon == MuonSubdetId::DT ) {
00347                                         DTWireId id(detid.rawId());
00348                                         station         = id.station();
00349                                         layer           = id.layer();
00350                                         superLayer      = id.superLayer();
00351                                         wheel           = id.wheel();
00352                                         sector          = id.sector();
00353                                         if (debug_) printf("DT \t[station][layer][superlayer]:[%d][%d][%d]\n", station,layer,superLayer);
00354                                 
00355                                 }
00356                                 else if ( systemMuon == MuonSubdetId::RPC) {
00357                                         RPCDetId id(detid.rawId());
00358                                         station         = id.station();
00359                                         if (debug_) printf("RPC\t[station]:[%d]\n", station);
00360                                 }
00361 
00362                         
00363                                 storageRecMuon_.muonLink_       [iMuonHit]      = iTrack;                       
00364                                 storageRecMuon_.system_         [iMuonHit]      = systemMuon;
00365                                 storageRecMuon_.endcap_         [iMuonHit]      = endcap;
00366                                 storageRecMuon_.station_        [iMuonHit]      = station;
00367                                 storageRecMuon_.ring_           [iMuonHit]      = ring;
00368                                 storageRecMuon_.chamber_        [iMuonHit]      = chamber;
00369                                 storageRecMuon_.layer_          [iMuonHit]      = layer;
00370                                 storageRecMuon_.superLayer_     [iMuonHit]      = superLayer;
00371                                 storageRecMuon_.wheel_          [iMuonHit]      = wheel;
00372                                 storageRecMuon_.sector_         [iMuonHit]      = sector;
00373                 
00374                                 storageRecMuon_.gpX_            [iMuonHit]      = gpRecX;
00375                                 storageRecMuon_.gpY_            [iMuonHit]      = gpRecY;
00376                                 storageRecMuon_.gpZ_            [iMuonHit]      = gpRecZ;
00377                                 storageRecMuon_.gpEta_          [iMuonHit]      = gpRecEta;
00378                                 storageRecMuon_.gpPhi_          [iMuonHit]      = gpRecPhi;
00379                                 storageRecMuon_.lpX_            [iMuonHit]      = lpX;
00380                                 storageRecMuon_.lpY_            [iMuonHit]      = lpY;
00381                                 storageRecMuon_.lpZ_            [iMuonHit]      = lpZ;
00382                                 iMuonHit++;
00383         
00384                         }
00385                         else if (detid.det() == DetId::Tracker) {
00386                                 
00387                                 if (debug_) printf("Tracker\n");
00388 
00389                                 StoreTrackerRecHits(detid, tTopo, iTrack, iTrackHit);
00390 
00391                                 storageTrackHit_.gpX_           [iTrackHit]     = gpRecX;
00392                                 storageTrackHit_.gpY_           [iTrackHit]     = gpRecY;
00393                                 storageTrackHit_.gpZ_           [iTrackHit]     = gpRecZ;
00394                                 storageTrackHit_.gpEta_         [iTrackHit]     = gpRecEta;
00395                                 storageTrackHit_.gpPhi_         [iTrackHit]     = gpRecPhi;
00396                                 storageTrackHit_.lpX_           [iTrackHit]     = lpX;
00397                                 storageTrackHit_.lpY_           [iTrackHit]     = lpY;
00398                                 storageTrackHit_.lpZ_           [iTrackHit]     = lpZ;
00399                                 iTrackHit++;
00400                          }              
00401                         else printf("THIS CAN NOT HAPPEN\n");           
00402                         
00403                         trkExtrap(detid, numTracks, iTrack, iRec, recoStart, lp, trackExtrap);
00404                         numTracks++;
00405 
00406                         if (debug_) printf("\tLocal Positon:  \tx = %2.2f\ty = %2.2f\tz = %2.2f\n",lpX, lpY, lpZ);
00407                         if (debug_) printf("\tGlobal Position: \tx = %6.2f\ty = %6.2f\tz = %6.2f\teta = %4.2f\tphi = %3.2f\n",
00408                                                         gpRecX,gpRecY,gpRecZ,gpRecEta,gpRecPhi);
00409                                 
00410 
00411                 }
00412 
00413         }
00414 
00415         storageRecMuon_ .n_     = iMuonHit; 
00416         storageTrackHit_.n_             = iTrackHit;
00417         trackExtrap             .n_     = numTracks;
00418 
00419 }
00420 //
00421 // Deal with Re-Fitted Track with some station omitted. 
00422 //
00423 // This should take the new track, match it to its track before refitting with the hits dumped, and extrapolate out to the
00424 //  rec hits that were removed from the fit.  
00425 //
00426 //
00427 
00428 void ResidualRefitting::NewTrackMeasurements(edm::Handle<reco::TrackCollection> trackCollOrig,
00429          edm::Handle<reco::TrackCollection> trackColl,  ResidualRefitting::storage_trackExtrap& trackExtrap) {
00430 
00431         int numTracks   = 0;    
00432         int recCounter  = 0;
00433 
00434         for ( reco::TrackCollection::const_iterator muon = trackColl->begin(); muon!=trackColl->end(); muon++) {
00435 
00436                 int iTrack = muon - trackColl->begin(); 
00437         
00438             reco::TrackRef trackref=reco::TrackRef(trackColl,iTrack);
00439                 FreeTrajectoryState recoStart = ResidualRefitting::freeTrajStateMuon(trackref);
00440 
00441                 int iTrackLink =  MatchTrackWithRecHits(muon,trackCollOrig);
00442                 reco::TrackRef ref = reco::TrackRef(trackCollOrig, iTrackLink);
00443 
00444                 for (trackingRecHit_iterator rec1 = ref->recHitsBegin(); rec1!= ref->recHitsEnd(); rec1++, recCounter++) {
00445 
00446                         //int iRec = rec1 - ref->recHitsBegin();
00447                         
00448                         bool unbiasedRec = true;
00449 
00450                         for ( trackingRecHit_iterator rec2 = muon->recHitsBegin(); rec2!=muon->recHitsEnd();rec2++) {
00451                                                 
00452                                 if (IsSameHit(rec1,rec2)) {
00453                                         unbiasedRec = false;
00454                                         break;
00455                                 }
00456                         }
00457                         if (!unbiasedRec) continue;
00458                                 
00459                         DetId detid = (*rec1)->geographicalId(); 
00460 
00461                         MuonTransientTrackingRecHit::MuonRecHitPointer mrhp =   
00462                                 MuonTransientTrackingRecHit::specificBuild(theService->trackingGeometry()->idToDet((**rec1).geographicalId())
00463                                 ,&(**rec1)); 
00464                 
00465                         trkExtrap(detid, numTracks, iTrackLink, recCounter, recoStart, (*rec1)->localPosition(), trackExtrap);
00466                         numTracks++;    
00467                 
00468                 }       
00469 
00470         }
00471 
00472         trackExtrap.n_ = numTracks;
00473 
00474 }
00475 //
00476 // Find the original track that corresponds to the re-fitted track
00477 //
00478 int ResidualRefitting::MatchTrackWithRecHits(reco::TrackCollection::const_iterator trackIt,
00479          edm::Handle<reco::TrackCollection> ref) {
00480 
00481         if (debug_) printf("Matching a re-fitted track to the original track.\n");
00482         
00483         int TrackMatch = -1;
00484 
00485         for (trackingRecHit_iterator rec = trackIt->recHitsBegin(); rec!=trackIt->recHitsEnd(); rec++) {
00486                 
00487                 bool foundMatch = false;
00488                 for (reco::TrackCollection::const_iterator refIt = ref->begin(); refIt!=ref->end(); refIt++) {
00489                 
00490                         int iTrackMatch = refIt - ref->begin();
00491                         if (foundMatch && TrackMatch !=iTrackMatch) break;
00492                         for (trackingRecHit_iterator recRef = refIt->recHitsBegin(); recRef!=refIt->recHitsEnd(); recRef++) {
00493 
00494                                 if (!IsSameHit(rec,recRef)) continue;
00495                                 
00496                                 foundMatch = true;
00497                                 TrackMatch = iTrackMatch;
00498                         //      printf("Rec hit match for original track %d\n", iTrackMatch);
00499                 
00500                         }
00501                 }
00502                 if (!foundMatch) {
00503                         printf("SOMETHING WENT WRONG! Could not match Track with original track!");
00504                         exit(1);
00505                 }
00506                 
00507         }
00508         if (debug_) printf("Rec hit match for original track %d\n", TrackMatch);
00509 
00510 //      reco::TrackRef trackref=reco::TrackRef(ref,TrackMatch);
00511         return TrackMatch;
00512 }
00513 /*
00514 //
00515 // Match two tracks to see if one is a subset of the other
00516 //
00517 
00518 bool ResidualRefitting::TrackSubset(reco::TrackRef trackSub, reco::TrackRef trackTop) {
00519 
00520         
00521         bool matchAll = true;
00522 
00523         for (trackingRecHit_iterator recSub = trackSub->recHitsBegin(); recSub!=trackSub->recHitsEnd(); recSub++) {
00524 
00525                 bool matchSub = false;
00526 
00527 
00528                 for (trackingRecHit_iterator recTop = trackTop->recHitsBegin(); recTop!=trackTop->recHitsEnd(); recTop++) {
00529                 
00530                         if ( recSub == recTop ) matchSub = true;
00531                         if (matchSub) break;
00532 
00533                 }
00534                 if (!matchSub) return false;
00535 
00536         }
00537 
00538         return matchAll;
00539 
00540 }
00541 */
00542 
00543 //
00544 // Check to see if the rec hits are the same
00545 //
00546 bool ResidualRefitting::IsSameHit(trackingRecHit_iterator hit1, trackingRecHit_iterator hit2) {
00547 
00548         
00549         double lpx1 = (*hit1)->localPosition().x();
00550         double lpy1 = (*hit1)->localPosition().y();
00551         double lpz1 = (*hit1)->localPosition().z();
00552 
00553         double lpx2 = (*hit2)->localPosition().x();
00554         double lpy2 = (*hit2)->localPosition().y();
00555         double lpz2 = (*hit2)->localPosition().z();
00556         if ( fabs( lpx1 - lpx2) > 1e-3) return false;
00557 //      printf("Match lpx...\n");
00558         if ( fabs( lpy1 - lpy2) > 1e-3) return false;
00559 //      printf("Match lpy...\n");
00560         if ( fabs( lpz1 - lpz2) > 1e-3) return false;
00561 //      printf("Match lpz...\n");
00562 
00563         return true;
00564 
00565 }
00566 
00567 //
00568 // Store Tracker Rec Hits
00569 //
00570 void ResidualRefitting::StoreTrackerRecHits(DetId detid, const TrackerTopology* tTopo, int iTrack, int iRec) {
00571 
00572 
00573                         int detector    = -1;
00574                         int subdetector = -1;
00575                         int blade               = -1;
00576                         int disk                = -1;
00577                         int ladder              = -1;
00578                         int layer               = -1;
00579                         int module              = -1;
00580                         int panel               = -1;
00581                         int ring                = -1;
00582                         int side                = -1;
00583                         int wheel               = -1;
00584                         
00585 //Detector Info
00586 
00587                         detector = detid.det();
00588                         subdetector = detid.subdetId();
00589                         
00590                         if (detector != DetId::Tracker) { 
00591                                 std::cout<<"OMFG NOT THE TRACKER\n"<<std::endl;
00592                                 return;
00593                         }
00594 
00595                         if (debug_) std::cout<<"Tracker:: ";
00596                         if (subdetector == ResidualRefitting::PXB) {
00597                                 layer   = tTopo->pxbLayer(detid.rawId());
00598                                 ladder  = tTopo->pxbLadder(detid.rawId());
00599                                 module  = tTopo->pxbModule(detid.rawId());
00600                                 if (debug_)     std::cout       <<      "PXB"
00601                                                                                 <<      "\tlayer = "    << layer
00602                                                                                 <<      "\tladder = "   << ladder
00603                                                                                 <<      "\tmodule = "   << module;
00604                                 
00605                         } 
00606                         else if (subdetector == ResidualRefitting::PXF) {
00607                                 side    = tTopo->pxfSide(detid.rawId());
00608                                 disk    = tTopo->pxfDisk(detid.rawId());
00609                                 blade   = tTopo->pxfBlade(detid.rawId());
00610                                 panel   = tTopo->pxfPanel(detid.rawId());
00611                                 module  = tTopo->pxfModule(detid.rawId());
00612                                 if (debug_)     std::cout       <<  "PXF"
00613                                                                                 <<      "\tside = "             << side
00614                                                                                 <<      "\tdisk = "             << disk
00615                                                                                 <<      "\tblade = "    << blade
00616                                                                                 <<      "\tpanel = "    << panel
00617                                                                                 <<      "\tmodule = "   << module;
00618                                                         
00619                         }
00620                         else if (subdetector == ResidualRefitting::TIB) {
00621                                 layer   = tTopo->tibLayer(detid.rawId());
00622                                 module  = tTopo->tibModule(detid.rawId());
00623                                 if (debug_)     std::cout       << "TIB"
00624                                                                                 << "\tlayer = " << layer
00625                                                                                 << "\tmodule = "<< module;
00626                         }
00627                         else if (subdetector == ResidualRefitting::TID) {
00628                                 side    = tTopo->tidSide(detid.rawId());
00629                                 wheel   = tTopo->tidWheel(detid.rawId());
00630                                 ring    = tTopo->tidRing(detid.rawId());
00631                                 if (debug_)     std::cout       <<"TID"
00632                                                                                 << "\tside = "  << side
00633                                                                                 << "\twheel = " << wheel
00634                                                                                 << "\tring = "  << ring;
00635                         
00636                         }
00637                         else if (subdetector == ResidualRefitting::TOB) {
00638                                 layer   = tTopo->tobLayer(detid.rawId());
00639                                 module  = tTopo->tobModule(detid.rawId());
00640                                 if (debug_)     std::cout       <<"TOB"
00641                                                                                 <<"\tlayer = "  << layer
00642                                                                                 <<"\tmodule = " << module;
00643                         
00644                         }
00645                         else if (subdetector == ResidualRefitting::TEC) {
00646                                 ring    = tTopo->tecRing(detid.rawId());
00647                                 module  = tTopo->tecModule(detid.rawId());
00648                                 if (debug_)     std::cout       <<"TEC"
00649                                                                                 << "\tring = "  << ring
00650                                                                                 << "\tmodule = "<< module;
00651                         }
00652 
00653 
00654 //Do Storage
00655 
00656                         storageTrackHit_.muonLink_              [iRec] =iTrack          ;
00657                         storageTrackHit_.detector_              [iRec] =detector        ;
00658                         storageTrackHit_.subdetector_           [iRec] =subdetector     ;
00659                         storageTrackHit_.blade_                 [iRec] =blade           ;
00660                         storageTrackHit_.disk_                  [iRec] =disk            ;
00661                         storageTrackHit_.ladder_                [iRec] =ladder          ;
00662                         storageTrackHit_.layer_                 [iRec] =layer           ;
00663                         storageTrackHit_.module_                [iRec] =module          ;
00664                         storageTrackHit_.panel_                 [iRec] =panel           ;
00665                         storageTrackHit_.ring_                  [iRec] =ring            ;
00666                         storageTrackHit_.side_                  [iRec] =side            ;
00667                         storageTrackHit_.wheel_                 [iRec] =wheel           ;
00668                         
00669 }
00670 
00671 //
00672 // Store Muon info on P, Pt, eta, phi
00673 //
00674 void ResidualRefitting::muonInfo(ResidualRefitting::storage_muon& storeMuon, reco::TrackRef muon, int val) {
00675 
00676 
00677         storeMuon.pt_ [val]                     = muon->pt();
00678     storeMuon.p_  [val]                 = muon->p();
00679     storeMuon.eta_[val]                 = muon->eta();
00680         storeMuon.phi_[val]                     = muon->phi();
00681         storeMuon.charge_[val]          = muon->charge();
00682         storeMuon.numRecHits_[val]      = muon->numberOfValidHits();
00683         storeMuon.chiSq_[val]           = muon->chi2();
00684         storeMuon.ndf_[val]                     = muon->ndof();
00685         storeMuon.chiSqOvrNdf_[val]     = muon->normalizedChi2();
00686 
00687 }
00688 // 
00689 // Fill a track extrapolation 
00690 // 
00691 void ResidualRefitting::trkExtrap(const DetId& detid, int iTrk, int iTrkLink, int iRec,
00692                                   const FreeTrajectoryState& freeTrajState,
00693                                   const LocalPoint& recPoint,
00694                                   storage_trackExtrap& storeTemp){
00695 
00696         bool dump_ = debug_;
00697         
00698         if (dump_) std::cout<< "In the trkExtrap function"<<std::endl;
00699 
00700         float gpExtrapX         = -99999;
00701         float gpExtrapY         = -99999;
00702         float gpExtrapZ         = -99999;
00703         float gpExtrapEta       = -99999;
00704         float gpExtrapPhi       = -99999;
00705 
00706         float lpX                       = -99999;
00707         float lpY                       = -99999;
00708         float lpZ                       = -99999;
00709 
00710         //
00711         // Get the local positions for the recHits
00712         //
00713 
00714         float recLpX = recPoint.x();
00715         float recLpY = recPoint.y();
00716         float recLpZ = recPoint.z();
00717 
00718         float resX = -9999;
00719         float resY = -9999;
00720         float resZ = -9999;
00721 
00722         const GeomDet* gdet = theService->trackingGeometry()->idToDet(detid);
00723         
00724 //      TrajectoryStateOnSurface surfTest =  prop.propagate(freeTrajState, gdet->surface());
00725         TrajectoryStateOnSurface surfTest =  thePropagator->propagate(freeTrajState, gdet->surface());
00726         
00727         if (surfTest.isValid()) {
00728         
00729                 GlobalPoint globTest    = surfTest.globalPosition();            
00730                 gpExtrapX                               = globTest.x();
00731                 gpExtrapY                               = globTest.y();
00732                 gpExtrapZ                               = globTest.z();
00733                 gpExtrapEta                             = globTest.eta();
00734                 gpExtrapPhi                             = globTest.phi();
00735                 LocalPoint loc                  = surfTest.localPosition();
00736                 if (detid.det() == DetId::Muon || detid.det() == DetId::Tracker) {
00737                         lpX                                             = loc.x();
00738                         lpY                                             = loc.y();
00739                         lpZ                                             = loc.z();
00740 
00741                         resX = lpX - recLpX;
00742                         resY = lpY - recLpY;
00743                         resZ = lpZ - recLpZ;
00744                         
00745                 }
00746 
00747         }
00748         storeTemp.muonLink_     [iTrk] = iTrkLink       ;
00749         storeTemp.recLink_      [iTrk] = iRec           ;
00750         storeTemp.gpX_          [iTrk] = gpExtrapX      ;
00751         storeTemp.gpY_          [iTrk] = gpExtrapY      ;
00752         storeTemp.gpZ_          [iTrk] = gpExtrapZ      ;
00753         storeTemp.gpEta_        [iTrk] = gpExtrapEta;
00754         storeTemp.gpPhi_        [iTrk] = gpExtrapPhi;   
00755         storeTemp.lpX_          [iTrk] = lpX            ;
00756         storeTemp.lpY_          [iTrk] = lpY            ;
00757         storeTemp.lpZ_          [iTrk] = lpZ            ;
00758         storeTemp.resX_         [iTrk] = resX           ;
00759         storeTemp.resY_         [iTrk] = resY           ;
00760         storeTemp.resZ_         [iTrk] = resZ           ;
00761         
00762         printf("station: %d\tsector: %d\tresX storage: %4.2f\n", ReturnStation(detid), ReturnSector(detid), resX);
00763 
00764 }
00765 //
00766 // Return the station 
00767 //
00768 int ResidualRefitting::ReturnStation(DetId detid) {
00769         
00770         int station             = -999;
00771 
00772         if (detid.det() == DetId::Muon) {
00773         
00774                 int systemMuon  = detid.subdetId(); // 1 DT; 2 CSC; 3 RPC
00775                 if ( systemMuon == MuonSubdetId::CSC) {
00776                         CSCDetId id(detid.rawId());
00777                         station         = id.station();
00778 
00779                 }
00780                 else if ( systemMuon == MuonSubdetId::DT ) {
00781                         DTWireId id(detid.rawId());
00782                         station         = id.station();
00783                 
00784                 }
00785                 else if ( systemMuon == MuonSubdetId::RPC) {
00786                         RPCDetId id(detid.rawId());
00787                         station         = id.station();
00788                 }
00789 
00790         }
00791         
00792         return station;         
00793 }
00794 //
00795 // Return the sector 
00796 //
00797 int ResidualRefitting::ReturnSector(DetId detid) {
00798         
00799         int sector              = -999;
00800 
00801         if (detid.det() == DetId::Muon) {
00802         
00803                 int systemMuon  = detid.subdetId(); // 1 DT; 2 CSC; 3 RPC
00804                 if ( systemMuon == MuonSubdetId::DT ) {
00805                         DTWireId id(detid.rawId());
00806                         sector          = id.sector();
00807                 
00808                 }
00809 
00810         }
00811         
00812         return sector;          
00813 }
00814 
00815 // 
00816 // Store the SAM and Track position info at a particular rho
00817 // 
00818 void ResidualRefitting::cylExtrapTrkSam(int recNum, reco::TrackRef track, ResidualRefitting::storage_trackExtrap& storage, double rho) {
00819 
00820         Cylinder::PositionType pos(0,0,0);
00821         Cylinder::RotationType rot;
00822 
00823         Cylinder::CylinderPointer myCylinder = Cylinder::build(pos, rot, rho);
00824 //      SteppingHelixPropagator inwardProp  ( theField, oppositeToMomentum );
00825 //      SteppingHelixPropagator outwardProp ( theField, alongMomentum );
00826         FreeTrajectoryState recoStart = freeTrajStateMuon(track);
00827 //      TrajectoryStateOnSurface recoProp = outwardProp.propagate(recoStart, *myCylinder);
00828         TrajectoryStateOnSurface recoProp = thePropagator->propagate(recoStart, *myCylinder);
00829 
00830         double xVal             = -9999;
00831         double yVal             = -9999;
00832         double zVal             = -9999;
00833         double phiVal   = -9999;
00834         double etaVal   = -9999;
00835 
00836         if(recoProp.isValid()) {
00837                 GlobalPoint recoPoint = recoProp.globalPosition();
00838                 xVal = recoPoint.x();
00839                 yVal = recoPoint.y();
00840                 zVal = recoPoint.z();
00841                 phiVal = recoPoint.phi();
00842                 etaVal = recoPoint.eta();               
00843         }
00844         storage.muonLink_[recNum]       = recNum;
00845         storage.gpX_    [recNum]        = xVal;
00846         storage.gpY_    [recNum]        = yVal;
00847         storage.gpZ_    [recNum]        = zVal;
00848         storage.gpEta_  [recNum]        = etaVal;
00849         storage.gpPhi_  [recNum]        = phiVal;
00850 
00851         float rhoVal = sqrt( xVal*xVal + yVal*yVal);
00852 
00853         printf("Cylinder: rho = %4.2f\tphi = %4.2f\teta = %4.2f\n", rhoVal, phiVal, etaVal);
00854         if (debug_) printf("Cylinder: rho = %4.2f\tphi = %4.2f\teta = %4.2f\n", rhoVal, phiVal, etaVal);
00855 
00856 
00857 }
00859 //Pre-Job junk
00861 
00862 //  
00863 // zero storage
00864 // 
00865 void ResidualRefitting::zero_storage() {
00866         if (debug_)     printf("zero_storage\n");
00867 
00868         zero_muon(&storageGmrOld_       );
00869         zero_muon(&storageGmrNew_       );
00870         zero_muon(&storageSamNew_       );
00871         zero_muon(&storageTrkNew_       );
00872         zero_muon(&storageGmrNoSt1_     );
00873         zero_muon(&storageSamNoSt1_     );
00874         zero_muon(&storageGmrNoSt2_     );
00875         zero_muon(&storageSamNoSt2_     );
00876         zero_muon(&storageGmrNoSt3_     );
00877         zero_muon(&storageSamNoSt3_     );
00878         zero_muon(&storageGmrNoSt4_     );
00879         zero_muon(&storageSamNoSt4_     );
00880 //zero out the tracker
00881         zero_muon(&storageGmrNoPXBLayer1        );
00882         zero_muon(&storageGmrNoPXBLayer2        );
00883         zero_muon(&storageGmrNoPXBLayer3        );
00884 
00885         zero_muon(&storageGmrNoPXF      );
00886 
00887         zero_muon(&storageGmrNoTIBLayer1        );
00888         zero_muon(&storageGmrNoTIBLayer2        );
00889         zero_muon(&storageGmrNoTIBLayer3        );
00890         zero_muon(&storageGmrNoTIBLayer4        );
00891 
00892         zero_muon(&storageGmrNoTID      );
00893 
00894         zero_muon(&storageGmrNoTOBLayer1        );
00895         zero_muon(&storageGmrNoTOBLayer2        );
00896         zero_muon(&storageGmrNoTOBLayer3        );
00897         zero_muon(&storageGmrNoTOBLayer4        );
00898         zero_muon(&storageGmrNoTOBLayer5        );
00899         zero_muon(&storageGmrNoTOBLayer6        );
00900 
00901         zero_muon(&storageGmrNoTEC      );
00902 
00903         zero_muon(&storageTrkNoPXBLayer1        );
00904         zero_muon(&storageTrkNoPXBLayer2        );
00905         zero_muon(&storageTrkNoPXBLayer3        );
00906 
00907         zero_muon(&storageTrkNoPXF      );
00908 
00909         zero_muon(&storageTrkNoTIBLayer1        );
00910         zero_muon(&storageTrkNoTIBLayer2        );
00911         zero_muon(&storageTrkNoTIBLayer3        );
00912         zero_muon(&storageTrkNoTIBLayer4        );
00913 
00914         zero_muon(&storageTrkNoTID      );
00915 
00916         zero_muon(&storageTrkNoTOBLayer1        );
00917         zero_muon(&storageTrkNoTOBLayer2        );
00918         zero_muon(&storageTrkNoTOBLayer3        );
00919         zero_muon(&storageTrkNoTOBLayer4        );
00920         zero_muon(&storageTrkNoTOBLayer5        );
00921         zero_muon(&storageTrkNoTOBLayer6        );
00922 
00923         zero_muon(&storageTrkNoTEC      );
00924 
00925         zero_trackExtrap(&storageTrackExtrapRec_                );
00926         zero_trackExtrap(&storageTrackExtrapTracker_    );
00927         zero_trackExtrap(&storageTrackExtrapRecNoSt1_   );
00928         zero_trackExtrap(&storageTrackExtrapRecNoSt2_   );
00929         zero_trackExtrap(&storageTrackExtrapRecNoSt3_   );
00930         zero_trackExtrap(&storageTrackExtrapRecNoSt4_   );
00931 
00932         zero_trackExtrap(&trackExtrap120_                               );
00933 
00934         zero_trackExtrap(&samExtrap120_                                 );
00935 
00936         zero_trackExtrap(&storageTrackNoPXBLayer1               );
00937         zero_trackExtrap(&storageTrackNoPXBLayer2               );
00938         zero_trackExtrap(&storageTrackNoPXBLayer3               );
00939 
00940         zero_trackExtrap(&storageTrackNoPXF                             );
00941 
00942         zero_trackExtrap(&storageTrackNoTIBLayer1               );
00943         zero_trackExtrap(&storageTrackNoTIBLayer2               );
00944         zero_trackExtrap(&storageTrackNoTIBLayer3               );
00945         zero_trackExtrap(&storageTrackNoTIBLayer4               );
00946 
00947         zero_trackExtrap(&storageTrackNoTOBLayer1               );
00948         zero_trackExtrap(&storageTrackNoTOBLayer2               );
00949         zero_trackExtrap(&storageTrackNoTOBLayer3               );
00950         zero_trackExtrap(&storageTrackNoTOBLayer4               );
00951         zero_trackExtrap(&storageTrackNoTOBLayer5               );
00952         zero_trackExtrap(&storageTrackNoTOBLayer6               );
00953 
00954         zero_trackExtrap(&storageTrackNoTEC                             );
00955 
00956         zero_trackExtrap(&storageTrackNoTID                             );
00957 
00958         storageRecMuon_         .n_ = 0;
00959         storageTrackHit_        .n_ = 0;
00960 
00961 }
00962 //
00963 // Zero out a muon reference
00964 //
00965 void ResidualRefitting::zero_muon(ResidualRefitting::storage_muon* str){
00966 
00967         str->n_ = 0;
00968         
00969         for (int i = 0; i < ResidualRefitting::N_MAX_STORED; i++) {
00970     
00971                 str->pt_                        [i] = -9999;
00972                 str->eta_                       [i] = -9999;
00973                 str->p_                         [i] = -9999;
00974                 str->phi_                       [i] = -9999;    
00975                 str->numRecHits_        [i] = -9999;
00976                 str->chiSq_             [i] = -9999;
00977                 str->ndf_                       [i] = -9999;
00978                 str->chiSqOvrNdf_       [i] = -9999;
00979         }
00980 
00981 }
00982 //
00983 // Zero track extrapolation
00984 //
00985 void ResidualRefitting::zero_trackExtrap(ResidualRefitting::storage_trackExtrap* str) {
00986         
00987         str->n_ = 0;
00988         for (int i = 0; i < ResidualRefitting::N_MAX_STORED_HIT; i++) {
00989 
00990                 str->muonLink_  [i] = -9999;
00991                 str->recLink_   [i] = -9999;
00992                 str->gpX_               [i] = -9999;
00993                 str->gpY_               [i] = -9999;
00994                 str->gpZ_               [i] = -9999;
00995                 str->gpEta_             [i] = -9999;
00996                 str->gpPhi_             [i] = -9999;
00997                 str->lpX_               [i] = -9999;
00998                 str->lpY_               [i] = -9999;
00999                 str->lpZ_               [i] = -9999;
01000                 str->resX_              [i] = -9999;
01001                 str->resY_              [i] = -9999;
01002                 str->resZ_              [i] = -9999;
01003         }
01004         
01005 }
01006 //  
01007 // Begin Job
01008 // 
01009 void ResidualRefitting::beginJob() {
01010 
01011         std::cout<<"Creating file "<< outputFileName_.c_str()<<std::endl;
01012 
01013         outputFile_  = new TFile( outputFileName_.c_str(), "RECREATE" ) ; 
01014 
01015         outputTree_ = new TTree("outputTree","outputTree");
01016 
01017         outputTree_->Branch("eventInfo",&eventInfo_,
01018                                                 "evtNum_/I:"
01019                                                 "runNum_/I"
01020                                 );
01021 
01022 
01023         ResidualRefitting::branchMuon(storageGmrOld_    , "gmrOld");
01024         ResidualRefitting::branchMuon(storageGmrNew_    , "gmrNew");
01025         ResidualRefitting::branchMuon(storageGmrNoSt1_  , "gmrNoSt1");
01026         ResidualRefitting::branchMuon(storageGmrNoSt2_  , "gmrNoSt2");
01027         ResidualRefitting::branchMuon(storageGmrNoSt3_  , "gmrNoSt3");
01028         ResidualRefitting::branchMuon(storageGmrNoSt4_  , "gmrNoSt4");
01029 
01030         ResidualRefitting::branchMuon(storageSamNew_    , "samNew");
01031         ResidualRefitting::branchMuon(storageSamNoSt1_  , "samNoSt1");
01032         ResidualRefitting::branchMuon(storageSamNoSt2_  , "samNoSt2");
01033         ResidualRefitting::branchMuon(storageSamNoSt3_  , "samNoSt3");
01034         ResidualRefitting::branchMuon(storageSamNoSt4_  , "samNoSt4");
01035 
01036         ResidualRefitting::branchMuon(storageTrkNew_    , "trkNew");
01037         ResidualRefitting::branchMuon(storageGmrNoPXBLayer1     , "gmrNoPXBLayer1");
01038         ResidualRefitting::branchMuon(storageGmrNoPXBLayer2     , "gmrNoPXBLayer2");
01039         ResidualRefitting::branchMuon(storageGmrNoPXBLayer3     , "gmrNoPXBLayer3");
01040         ResidualRefitting::branchMuon(storageGmrNoPXF   , "gmrNoPXF");
01041         ResidualRefitting::branchMuon(storageGmrNoTIBLayer1     , "gmrNoTIBLayer1");
01042         ResidualRefitting::branchMuon(storageGmrNoTIBLayer2     , "gmrNoTIBLayer2");
01043         ResidualRefitting::branchMuon(storageGmrNoTIBLayer3     , "gmrNoTIBLayer3");
01044         ResidualRefitting::branchMuon(storageGmrNoTIBLayer4     , "gmrNoTIBLayer4");
01045         ResidualRefitting::branchMuon(storageGmrNoTID   , "gmrNoTID");
01046         ResidualRefitting::branchMuon(storageGmrNoTOBLayer1     , "gmrNoTOBLayer1");
01047         ResidualRefitting::branchMuon(storageGmrNoTOBLayer2     , "gmrNoTOBLayer2");
01048         ResidualRefitting::branchMuon(storageGmrNoTOBLayer3     , "gmrNoTOBLayer3");
01049         ResidualRefitting::branchMuon(storageGmrNoTOBLayer4     , "gmrNoTOBLayer4");
01050         ResidualRefitting::branchMuon(storageGmrNoTOBLayer5     , "gmrNoTOBLayer5");
01051         ResidualRefitting::branchMuon(storageGmrNoTOBLayer6     , "gmrNoTOBLayer6");
01052         ResidualRefitting::branchMuon(storageGmrNoTEC   , "gmrNoTEC");
01053 
01054         ResidualRefitting::branchMuon(storageTrkNoPXBLayer1     , "trkNoPXBLayer1");
01055         ResidualRefitting::branchMuon(storageTrkNoPXBLayer2     , "trkNoPXBLayer2");
01056         ResidualRefitting::branchMuon(storageTrkNoPXBLayer3     , "trkNoPXBLayer3");
01057         ResidualRefitting::branchMuon(storageTrkNoPXF   , "trkNoPXF");
01058         ResidualRefitting::branchMuon(storageTrkNoTIBLayer1     , "trkNoTIBLayer1");
01059         ResidualRefitting::branchMuon(storageTrkNoTIBLayer2     , "trkNoTIBLayer2");
01060         ResidualRefitting::branchMuon(storageTrkNoTIBLayer3     , "trkNoTIBLayer3");
01061         ResidualRefitting::branchMuon(storageTrkNoTIBLayer4     , "trkNoTIBLayer4");
01062         ResidualRefitting::branchMuon(storageTrkNoTID   , "trkNoTID");
01063         ResidualRefitting::branchMuon(storageTrkNoTOBLayer1     , "trkNoTOBLayer1");
01064         ResidualRefitting::branchMuon(storageTrkNoTOBLayer2     , "trkNoTOBLayer2");
01065         ResidualRefitting::branchMuon(storageTrkNoTOBLayer3     , "trkNoTOBLayer3");
01066         ResidualRefitting::branchMuon(storageTrkNoTOBLayer4     , "trkNoTOBLayer4");
01067         ResidualRefitting::branchMuon(storageTrkNoTOBLayer5     , "trkNoTOBLayer5");
01068         ResidualRefitting::branchMuon(storageTrkNoTOBLayer6     , "trkNoTOBLayer6");
01069         ResidualRefitting::branchMuon(storageTrkNoTEC   , "trkNoTEC");
01070                                         
01071         outputBranch_ = outputTree_ -> Branch("recHitsNew", &storageRecMuon_, 
01072 
01073                 "n_/I:"
01074                 "muonLink_[1000]/I:"
01075                 
01076                 "system_[1000]/I:"
01077                 "endcap_[1000]/I:"
01078                 "station_[1000]/I:"
01079                 "ring_[1000]/I:"
01080                 "chamber_[1000]/I:"
01081                 "layer_[1000]/I:"
01082                 "superLayer_[1000]/I:"
01083                 "wheel_[1000]/I:"
01084                 "sector_[1000]/I:"
01085                 
01086         
01087                 "gpX_[1000]/F:"
01088                 "gpY_[1000]/F:"
01089                 "gpZ_[1000]/F:"
01090                 "gpEta_[1000]/F:"
01091                 "gpPhi_[1000]/F:"
01092                 "lpX_[1000]/F:"
01093                 "lpY_[1000]/F:"
01094                 "lpZ_[1000]/F"
01095                 );
01096                 
01097                 
01098         outputBranch_ = outputTree_ -> Branch("recHitsTracker", &storageTrackHit_,
01099         
01100                 "n_/I:"
01101                 
01102                 "muonLink_[1000]/I:"
01103                 "detector_[1000]/I:"
01104                 "subdetector_[1000]/I:"
01105                 "blade_[1000]/I:"
01106                 "disk_[1000]/I:"
01107                 "ladder_[1000]/I:"
01108                 "layer_[1000]/I:"
01109                 "module_[1000]/I:"
01110                 "panel_[1000]/I:"
01111                 "ring_[1000]/I:"
01112                 "side_[1000]/I:"
01113                 "wheel_[1000]/I:"
01114                                 
01115                 "gpX_[1000]/F:"
01116                 "gpY_[1000]/F:"
01117                 "gpZ_[1000]/F:"
01118                 "gpEta_[1000]/F:"
01119                 "gpPhi_[1000]/F:"
01120                 "lpX_[1000]/F:"
01121                 "lpY_[1000]/F:"
01122                 "lpZ_[1000]/F"
01123                 );      
01124                 
01125                 
01126         ResidualRefitting::branchTrackExtrap(storageTrackExtrapRec_             , "trkExtrap");
01127         ResidualRefitting::branchTrackExtrap(storageTrackExtrapRecNoSt1_        , "trkExtrapNoSt1");
01128         ResidualRefitting::branchTrackExtrap(storageTrackExtrapRecNoSt2_        , "trkExtrapNoSt2");
01129         ResidualRefitting::branchTrackExtrap(storageTrackExtrapRecNoSt3_        , "trkExtrapNoSt3");
01130         ResidualRefitting::branchTrackExtrap(storageTrackExtrapRecNoSt4_        , "trkExtrapNoSt4");
01131 
01132         ResidualRefitting::branchTrackExtrap(storageTrackExtrapTracker_, "trkExtrapTracker");
01133         ResidualRefitting::branchTrackExtrap(storageTrackNoPXF                  , "trkExtrapNoPXF");
01134         ResidualRefitting::branchTrackExtrap(storageTrackNoPXBLayer1                    , "trkExtrapNoPXBLayer1");
01135         ResidualRefitting::branchTrackExtrap(storageTrackNoPXBLayer2                    , "trkExtrapNoPXBLayer2");
01136         ResidualRefitting::branchTrackExtrap(storageTrackNoPXBLayer3                    , "trkExtrapNoPXBLayer3");
01137         ResidualRefitting::branchTrackExtrap(storageTrackNoTIBLayer1                    , "trkExtrapNoTIBLayer1");
01138         ResidualRefitting::branchTrackExtrap(storageTrackNoTIBLayer2                    , "trkExtrapNoTIBLayer2");
01139         ResidualRefitting::branchTrackExtrap(storageTrackNoTIBLayer3                    , "trkExtrapNoTIBLayer3");
01140         ResidualRefitting::branchTrackExtrap(storageTrackNoTIBLayer4                    , "trkExtrapNoTIBLayer4");
01141         ResidualRefitting::branchTrackExtrap(storageTrackNoTID                  , "trkExtrapNoTID");
01142         ResidualRefitting::branchTrackExtrap(storageTrackNoTOBLayer1                    , "trkExtrapNoTOBLayer1");
01143         ResidualRefitting::branchTrackExtrap(storageTrackNoTOBLayer2                    , "trkExtrapNoTOBLayer2");
01144         ResidualRefitting::branchTrackExtrap(storageTrackNoTOBLayer3                    , "trkExtrapNoTOBLayer3");
01145         ResidualRefitting::branchTrackExtrap(storageTrackNoTOBLayer4                    , "trkExtrapNoTOBLayer4");
01146         ResidualRefitting::branchTrackExtrap(storageTrackNoTOBLayer5                    , "trkExtrapNoTOBLayer5");
01147         ResidualRefitting::branchTrackExtrap(storageTrackNoTOBLayer6                    , "trkExtrapNoTOBLayer6");
01148         ResidualRefitting::branchTrackExtrap(storageTrackNoTEC                  , "trkExtrapNoTEC");
01149 
01150         ResidualRefitting::branchTrackExtrap(trackExtrap120_                    , "trackCyl120");
01151         ResidualRefitting::branchTrackExtrap(samExtrap120_                              , "samCyl120");
01152 
01153 }
01154 //
01155 // Set the Muon Branches
01156 //
01157 void ResidualRefitting::branchMuon(ResidualRefitting::storage_muon& storageTmp, std::string branchName){
01158 
01159         outputBranch_ = outputTree_ -> Branch(branchName.c_str(), &storageTmp, 
01160                                                                         "n_/I:"
01161                                                                         "charge_[10]/I:"
01162                                                                         "pt_[10]/F:"
01163                                                                         "eta_[10]/F:"
01164                                                                         "p_[10]/F:"
01165                                                                         "phi_[10]/F:"
01166                                                                         "numRecHits_[10]/I:"
01167                                                                         "chiSq_[10]/F:"
01168                                                                         "ndf_[10]/F:"
01169                                                                         "chiSqOvrNdf_[10]/F"
01170                                                                                 
01171                                                                                 );
01172 
01173 }
01174 //
01175 // Set the Branches for Track Extrapolations
01176 //
01177 void ResidualRefitting::branchTrackExtrap(ResidualRefitting::storage_trackExtrap& storageTmp, std::string branchName){
01178 
01179         outputBranch_ = outputTree_ -> Branch(branchName.c_str(), &storageTmp, 
01180                                                                         "n_/I:"
01181                                                                         "muonLink_[1000]/I:"
01182                                                                         "recLink_[1000]/I:"
01183                                                                         "gpX_[1000]/F:"
01184                                                                         "gpY_[1000]/F:"
01185                                                                         "gpZ_[1000]/F:"
01186                                                                         "gpEta_[1000]/F:"
01187                                                                         "gpPhi_[1000]/F:"
01188                                                                         "lpX_[1000]/F:"
01189                                                                         "lpY_[1000]/F:"
01190                                                                         "lpZ_[1000]/F:"
01191                                                                         "resX_[1000]/F:"
01192                                                                         "resY_[1000]/F:"
01193                                                                         "resZ_[1000]/F"
01194                                                                                 
01195                                                                         );
01196 
01197 }
01198 // 
01199 // End Job
01200 // 
01201 void ResidualRefitting::endJob() {
01202 
01203 
01204   outputFile_ -> Write();
01205 
01206   outputFile_ -> Close() ;
01207 
01208 }
01209 // 
01210 // Return a Free Trajectory state for a muon track
01211 // 
01212 FreeTrajectoryState ResidualRefitting::freeTrajStateMuon(reco::TrackRef muon){
01213                                         
01214                         math::XYZPoint  innerPos = muon -> referencePoint();
01215                         math::XYZVector innerMom = muon -> momentum();
01216                         if (debug_) std::cout   <<  "Inner Pos: "
01217                                                                         <<      "\tx = "        << innerPos.X()
01218                                                                         <<      "\ty = "        << innerPos.Y()
01219                                                                         <<      "\tz = "        << innerPos.Z()
01220                                                                         <<      std::endl;
01221                                                                 
01222                         GlobalPoint   innerPoint( innerPos.X(), innerPos.Y(), innerPos.Z());
01223                         GlobalVector  innerVec  ( innerMom.X(), innerMom.Y(), innerMom.Z());
01224 
01225                         FreeTrajectoryState recoStart( innerPoint, innerVec, muon ->charge(), theField ); 
01226                         return recoStart;
01227 
01228 }
01229 
01231 // nTuple value Dumps
01233 
01234 // 
01235 // dump Track Extrapolation
01236 //
01237 void ResidualRefitting::dumpTrackExtrap(ResidualRefitting::storage_trackExtrap track) {
01238         std::cout<<"\n\nExtrapolation Dump:\n";
01239         for (unsigned int i = 0; i < (unsigned int)track.n_; i++) {
01240                 
01241 //              double rho = sqrt( (float)track.gpX_[i] * (float)track.gpX_[i] + (float)track.gpY_[i] * (float)track.gpY_[i]  );
01242 
01243                 printf ("%d\tmuonLink= %d",i, (int)track.muonLink_[i]);         
01244                 printf ("\trecLink = %d",       (int)track.recLink_[i]  );
01245 //              printf ("\tGlobal\tx = %0.3f"           ,               (float)track.gpX_[i]    );
01246 //              printf ("\ty = %0.3f"           ,               (float)track.gpY_[i]    );
01247 //              printf ("\tz = %0.3f"           ,               (float)track.gpZ_[i]    );
01248 //              printf ("\trho =%0.3f"          ,               rho                                             );
01249 //              printf ("\teta = %0.3f"         ,               (float)track.gpEta_[i]  );
01250 //              printf ("\tphi = %0.3f"         ,               (float)track.gpPhi_[i]  );
01251                 printf ("\t\tLocal\tx = %0.3f"  ,               (float)track.lpX_[i]    );
01252                 printf ("\ty = %0.3f"           ,       (float)track.lpY_[i]            );
01253                 printf ("\tz = %0.3f\n" ,       (float)track.lpZ_[i]            );
01254                 
01255         }
01256         
01257 }
01258 //
01259 // dump Muon Rec Hits
01260 //
01261 void ResidualRefitting::dumpMuonRecHits(ResidualRefitting::storage_hit hit) {
01262         std::cout<<"Muon Rec Hits Dump:\n";
01263         for (unsigned int i = 0; i < (unsigned int)hit.n_; i++) {
01264                 
01265 //              double rho = sqrt( (float)hit.gpX_[i] * (float)hit.gpX_[i] + (float)hit.gpY_[i] * (float)hit.gpY_[i]  );
01266                 
01267                 printf ("%d\tsubdetector = %d\t superLayer =%d" , i,    (int)hit.system_[i], (int)hit.superLayer_[i]    );
01268 //              printf ("\tGlobal\tx = %0.3f"                   ,               (float)hit.gpX_[i]                      );
01269 //              printf ("\ty = %0.3f"                           ,               (float)hit.gpY_[i]                      );
01270 //              printf ("\tz = %0.3f"                           ,               (float)hit.gpZ_[i]                      );
01271 //              printf ("\trho =%0.3f"                          ,               rho                                                     );
01272 //              printf ("\teta = %0.3f"                         ,               (float)hit.gpEta_[i]            );
01273 //              printf ("\tphi = %0.3f\n"                       ,               (float)hit.gpPhi_[i]            );
01274                 printf ("\t\tLocal\tx = %0.3f"          ,               (float)hit.lpX_[i]                      );
01275                 printf ("\ty = %0.3f"                           ,               (float)hit.lpY_[i]                      );
01276                 printf ("\tz = %0.3f\n"                 ,               (float)hit.lpZ_[i]                      );
01277                 
01278         }
01279         
01280 }
01281 //
01282 // dump Tracker Rec Hits
01283 //
01284 void ResidualRefitting::dumpTrackHits(ResidualRefitting::storage_trackHit hit) {
01285         std::cout<<"Tracker Rec Hits Dump:\n";
01286         for (unsigned int i = 0; i < (unsigned int)hit.n_; i++) {
01287                 
01288 //              double rho = sqrt( (float)hit.gpX_[i] * (float)hit.gpX_[i] + (float)hit.gpY_[i] * (float)hit.gpY_[i]  );
01289                 
01290                 printf ("%d\tsubdetector = %d"          , i,    (int)hit.subdetector_[i]        );
01291                 printf ("\tlayer = %d"                          ,       (int)hit.layer_[i]      );
01292 //              printf ("\tGlobal\tx = %0.3f"                   ,               (float)hit.gpX_[i]                      );
01293 //              printf ("\ty = %0.3f"                           ,               (float)hit.gpY_[i]                      );
01294 //              printf ("\tz = %0.3f"                           ,               (float)hit.gpZ_[i]                      );
01295 //              printf ("\trho =%0.3f"                          ,               rho                                                     );
01296 //              printf ("\teta = %0.3f"                         ,               (float)hit.gpEta_[i]            );
01297 //              printf ("\tphi = %0.3f\n"                       ,               (float)hit.gpPhi_[i]            );
01298                 printf ("\t\tLocal\tx = %0.3f"          ,               (float)hit.lpX_[i]                      );
01299                 printf ("\ty = %0.3f"                           ,               (float)hit.lpY_[i]                      );
01300                 printf ("\tz = %0.3f\n"                         ,               (float)hit.lpZ_[i]                      );
01301                 
01302         }
01303         
01304 }
01305 //
01306 //Dump a TrackRef 
01307 //
01308 void ResidualRefitting::dumpTrackRef(reco::TrackRef muon, std::string str) {
01309 
01310         float pt = muon->pt();
01311         float p  = muon->p  ();
01312         float eta = muon->eta();
01313         float phi = muon->phi();
01314         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);
01315 
01316 }
01317 
01318 DEFINE_FWK_MODULE( ResidualRefitting );
01319 
01320 
01322 //Deprecated
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/Alignment/OfflineValidation/plugins/ResidualRefitting.cc
