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