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