File indexing completed on 2024-08-13 05:00:11
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023 #include <memory>
0024 #include <string>
0025 #include <map>
0026 #include <vector>
0027 #include <fstream>
0028 #include <typeinfo>
0029
0030
0031
0032
0033 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0034 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0035 #include "FWCore/Framework/interface/ESWatcher.h"
0036 #include "FWCore/Framework/interface/Frameworkfwd.h"
0037 #include "FWCore/Framework/interface/Event.h"
0038 #include "FWCore/Framework/interface/EventSetup.h"
0039 #include "FWCore/Framework/interface/MakerMacros.h"
0040 #include "FWCore/Utilities/interface/InputTag.h"
0041 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0042 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0043 #include "FWCore/ServiceRegistry/interface/Service.h"
0044
0045
0046
0047 #include "DataFormats/DetId/interface/DetId.h"
0048 #include "DataFormats/MuonReco/interface/Muon.h"
0049 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0050 #include "DataFormats/TrackReco/interface/Track.h"
0051 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0052 #include "DataFormats/GeometrySurface/interface/TangentPlane.h"
0053 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0054 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0055 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0056 #include "MagneticField/Engine/interface/MagneticField.h"
0057 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0058
0059 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBuilder.h"
0060 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0061 #include "TrackingTools/TrackRefitter/interface/TrackTransformer.h"
0062 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0063 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0064 #include "TrackingTools/GeomPropagators/interface/SmartPropagator.h"
0065 #include "TrackingTools/GeomPropagators/interface/PropagationDirectionChooser.h"
0066 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0067 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
0068 #include "TrackPropagation/RungeKutta/interface/defaultRKPropagator.h"
0069 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0070 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0071 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCartesianToLocal.h"
0072
0073
0074 #include "FWCore/Framework/interface/ESHandle.h"
0075 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0076
0077
0078 #include "CondFormats/Alignment/interface/Alignments.h"
0079 #include "CondFormats/Alignment/interface/AlignTransform.h"
0080 #include "CondFormats/AlignmentRecord/interface/GlobalPositionRcd.h"
0081
0082
0083 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0084 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0085 #include "TrackingTools/TrackFitters/interface/KFTrajectoryFitter.h"
0086 #include "TrackingTools/TrackFitters/interface/KFTrajectorySmoother.h"
0087 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0088 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0089
0090 #include <TH1.h>
0091 #include <TH2.h>
0092 #include <TFile.h>
0093
0094 #include <iostream>
0095 #include <cmath>
0096 #include <cassert>
0097 #include <fstream>
0098 #include <iomanip>
0099
0100 using namespace edm;
0101 using namespace std;
0102 using namespace reco;
0103
0104
0105
0106
0107
0108 class GlobalTrackerMuonAlignment : public edm::one::EDAnalyzer<> {
0109 public:
0110 explicit GlobalTrackerMuonAlignment(const edm::ParameterSet&);
0111 ~GlobalTrackerMuonAlignment() override = default;
0112 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0113
0114 void analyzeTrackTrack(const edm::Event&, const edm::EventSetup&);
0115 void analyzeTrackTrajectory(const edm::Event&, const edm::EventSetup&);
0116 void bookHist();
0117 void fitHist();
0118 void trackFitter(reco::TrackRef, reco::TransientTrack&, PropagationDirection, TrajectoryStateOnSurface&);
0119 void muonFitter(reco::TrackRef, reco::TransientTrack&, PropagationDirection, TrajectoryStateOnSurface&);
0120 void debugTrackHit(const std::string, reco::TrackRef);
0121 void debugTrackHit(const std::string, reco::TransientTrack&);
0122 void debugTrajectorySOS(const std::string, TrajectoryStateOnSurface&);
0123 void debugTrajectorySOSv(const std::string, TrajectoryStateOnSurface);
0124 void debugTrajectory(const std::string, Trajectory&);
0125
0126 void gradientGlobal(GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&, AlgebraicSymMatrix66&);
0127 void gradientLocal(GlobalVector&,
0128 GlobalVector&,
0129 GlobalVector&,
0130 GlobalVector&,
0131 GlobalVector&,
0132 CLHEP::HepSymMatrix&,
0133 CLHEP::HepMatrix&,
0134 CLHEP::HepVector&,
0135 AlgebraicVector4&);
0136 void gradientGlobalAlg(GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&, AlgebraicSymMatrix66&);
0137 void misalignMuon(GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&);
0138 void misalignMuonL(GlobalVector&,
0139 GlobalVector&,
0140 GlobalVector&,
0141 GlobalVector&,
0142 GlobalVector&,
0143 GlobalVector&,
0144 AlgebraicVector4&,
0145 TrajectoryStateOnSurface&,
0146 TrajectoryStateOnSurface&);
0147 void writeGlPosRcd(CLHEP::HepVector& d3);
0148 inline double CLHEP_dot(const CLHEP::HepVector& a, const CLHEP::HepVector& b) {
0149 return a(1) * b(1) + a(2) * b(2) + a(3) * b(3);
0150 }
0151
0152 private:
0153 void beginJob() override;
0154 void analyze(const edm::Event&, const edm::EventSetup&) override;
0155 void endJob() override;
0156
0157
0158 const edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> m_TkGeometryToken;
0159 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> m_MagFieldToken;
0160 const edm::ESGetToken<Alignments, GlobalPositionRcd> m_globalPosToken;
0161 const edm::ESGetToken<Propagator, TrackingComponentsRecord> m_propToken;
0162 const edm::ESGetToken<TransientTrackingRecHitBuilder, TransientRecHitRecord> m_ttrhBuilderToken;
0163
0164 const edm::InputTag trackTags_;
0165 const edm::InputTag muonTags_;
0166 const edm::InputTag gmuonTags_;
0167 const edm::InputTag smuonTags_;
0168 const std::string propagator_;
0169 bool cosmicMuonMode_;
0170 bool isolatedMuonMode_;
0171 bool collectionCosmic;
0172 bool collectionIsolated;
0173 bool refitMuon_;
0174 bool refitTrack_;
0175 const std::string rootOutFile_;
0176 const std::string txtOutFile_;
0177 const bool writeDB_;
0178 const bool debug_;
0179
0180 const edm::EDGetTokenT<reco::TrackCollection> trackIsolateToken_;
0181 const edm::EDGetTokenT<reco::TrackCollection> muonIsolateToken_;
0182 const edm::EDGetTokenT<reco::TrackCollection> gmuonIsolateToken_;
0183 const edm::EDGetTokenT<reco::MuonCollection> smuonIsolateToken_;
0184 const edm::EDGetTokenT<reco::TrackCollection> trackCosmicToken_;
0185 const edm::EDGetTokenT<reco::TrackCollection> muonCosmicToken_;
0186 const edm::EDGetTokenT<reco::TrackCollection> gmuonCosmicToken_;
0187 const edm::EDGetTokenT<reco::MuonCollection> smuonCosmicToken_;
0188 const edm::EDGetTokenT<reco::TrackCollection> trackToken_;
0189 const edm::EDGetTokenT<reco::TrackCollection> muonToken_;
0190 const edm::EDGetTokenT<reco::TrackCollection> gmuonToken_;
0191 const edm::EDGetTokenT<reco::MuonCollection> smuonToken_;
0192
0193 edm::ESWatcher<GlobalTrackingGeometryRecord> watchTrackingGeometry_;
0194 edm::ESHandle<GlobalTrackingGeometry> theTrackingGeometry;
0195 const GlobalTrackingGeometry* trackingGeometry_;
0196
0197 edm::ESWatcher<IdealMagneticFieldRecord> watchMagneticFieldRecord_;
0198 const MagneticField* magneticField_;
0199
0200 edm::ESWatcher<TrackingComponentsRecord> watchTrackingComponentsRecord_;
0201
0202 edm::ESWatcher<GlobalPositionRcd> watchGlobalPositionRcd_;
0203 const Alignments* globalPositionRcd_;
0204
0205 KFTrajectoryFitter* theFitter;
0206 KFTrajectorySmoother* theSmoother;
0207 KFTrajectoryFitter* theFitterOp;
0208 KFTrajectorySmoother* theSmootherOp;
0209 bool defineFitter;
0210 MuonTransientTrackingRecHitBuilder* MuRHBuilder;
0211 const TransientTrackingRecHitBuilder* TTRHBuilder;
0212
0213
0214 AlgebraicVector3 Gfr;
0215 AlgebraicSymMatrix33 Inf;
0216
0217 CLHEP::HepVector Grad;
0218 CLHEP::HepMatrix Hess;
0219
0220 CLHEP::HepVector GradL;
0221 CLHEP::HepMatrix HessL;
0222
0223 int N_event;
0224 int N_track;
0225
0226 std::vector<AlignTransform>::const_iterator
0227 iteratorTrackerRcd;
0228 std::vector<AlignTransform>::const_iterator iteratorMuonRcd;
0229 std::vector<AlignTransform>::const_iterator iteratorEcalRcd;
0230 std::vector<AlignTransform>::const_iterator iteratorHcalRcd;
0231
0232 CLHEP::HepVector MuGlShift;
0233 CLHEP::HepVector MuGlAngle;
0234
0235 std::string MuSelect;
0236
0237 ofstream OutGlobalTxt;
0238
0239 TFile* file;
0240 TH1F* histo;
0241 TH1F* histo2;
0242 TH1F* histo3;
0243 TH1F* histo4;
0244 TH1F* histo5;
0245 TH1F* histo6;
0246 TH1F* histo7;
0247 TH1F* histo8;
0248
0249 TH1F* histo11;
0250 TH1F* histo12;
0251 TH1F* histo13;
0252 TH1F* histo14;
0253 TH1F* histo15;
0254 TH1F* histo16;
0255 TH1F* histo17;
0256 TH1F* histo18;
0257 TH1F* histo19;
0258 TH1F* histo20;
0259 TH1F* histo21;
0260 TH1F* histo22;
0261 TH1F* histo23;
0262 TH1F* histo24;
0263 TH1F* histo25;
0264 TH1F* histo26;
0265 TH1F* histo27;
0266 TH1F* histo28;
0267 TH1F* histo29;
0268 TH1F* histo30;
0269 TH1F* histo31;
0270 TH1F* histo32;
0271 TH1F* histo33;
0272 TH1F* histo34;
0273 TH1F* histo35;
0274
0275 TH2F* histo101;
0276 TH2F* histo102;
0277 };
0278
0279
0280
0281
0282
0283
0284
0285
0286
0287
0288
0289
0290 GlobalTrackerMuonAlignment::GlobalTrackerMuonAlignment(const edm::ParameterSet& iConfig)
0291 : m_TkGeometryToken(esConsumes()),
0292 m_MagFieldToken(esConsumes()),
0293 m_globalPosToken(esConsumes()),
0294 m_propToken(esConsumes(edm::ESInputTag("", iConfig.getParameter<string>("Propagator")))),
0295 m_ttrhBuilderToken(esConsumes(edm::ESInputTag("", "witTrackAngle"))),
0296 trackTags_(iConfig.getParameter<edm::InputTag>("tracks")),
0297 muonTags_(iConfig.getParameter<edm::InputTag>("muons")),
0298 gmuonTags_(iConfig.getParameter<edm::InputTag>("gmuons")),
0299 smuonTags_(iConfig.getParameter<edm::InputTag>("smuons")),
0300 cosmicMuonMode_(iConfig.getParameter<bool>("cosmics")),
0301 isolatedMuonMode_(iConfig.getParameter<bool>("isolated")),
0302 refitMuon_(iConfig.getParameter<bool>("refitmuon")),
0303 refitTrack_(iConfig.getParameter<bool>("refittrack")),
0304 rootOutFile_(iConfig.getUntrackedParameter<string>("rootOutFile")),
0305 txtOutFile_(iConfig.getUntrackedParameter<string>("txtOutFile")),
0306 writeDB_(iConfig.getUntrackedParameter<bool>("writeDB")),
0307 debug_(iConfig.getUntrackedParameter<bool>("debug")),
0308 trackIsolateToken_(consumes<reco::TrackCollection>(edm::InputTag("ALCARECOMuAlCalIsolatedMu", "TrackerOnly"))),
0309 muonIsolateToken_(consumes<reco::TrackCollection>(edm::InputTag("ALCARECOMuAlCalIsolatedMu", "StandAlone"))),
0310 gmuonIsolateToken_(consumes<reco::TrackCollection>(edm::InputTag("ALCARECOMuAlCalIsolatedMu", "GlobalMuon"))),
0311 smuonIsolateToken_(consumes<reco::MuonCollection>(edm::InputTag("ALCARECOMuAlCalIsolatedMu", "SelectedMuons"))),
0312 trackCosmicToken_(
0313 consumes<reco::TrackCollection>(edm::InputTag("ALCARECOMuAlCOMMuGlobalCosmicsMu", "TrackerOnly"))),
0314 muonCosmicToken_(
0315 consumes<reco::TrackCollection>(edm::InputTag("ALCARECOMuAlCOMMuGlobalCosmicsMu", "StandAlone"))),
0316 gmuonCosmicToken_(
0317 consumes<reco::TrackCollection>(edm::InputTag("ALCARECOMuAlCOMMuGlobalCosmicsMu", "GlobalMuon"))),
0318 smuonCosmicToken_(
0319 consumes<reco::MuonCollection>(edm::InputTag("ALCARECOMuAlCOMMuGlobalCosmicsMu", "SelectedMuons"))),
0320 trackToken_(consumes<reco::TrackCollection>(trackTags_)),
0321 muonToken_(consumes<reco::TrackCollection>(muonTags_)),
0322 gmuonToken_(consumes<reco::TrackCollection>(gmuonTags_)),
0323 smuonToken_(consumes<reco::MuonCollection>(smuonTags_)),
0324 defineFitter(true) {
0325 #ifdef NO_FALSE_FALSE_MODE
0326 if (cosmicMuonMode_ == false && isolatedMuonMode_ == false) {
0327 throw cms::Exception("BadConfig") << "Muon collection not set! "
0328 << "Use from GlobalTrackerMuonAlignment_test_cfg.py.";
0329 }
0330 #endif
0331 if (cosmicMuonMode_ == true && isolatedMuonMode_ == true) {
0332 throw cms::Exception("BadConfig") << "Muon collection can be either cosmic or isolated! "
0333 << "Please set only one to true.";
0334 }
0335 }
0336
0337 void GlobalTrackerMuonAlignment::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0338 edm::ParameterSetDescription desc;
0339 desc.add<std::string>("propagator", "SteppingHelixPropagator");
0340 desc.add<edm::InputTag>("tracks", edm::InputTag("ALCARECOMuAlGlobalCosmics:TrackerOnly"));
0341 desc.add<edm::InputTag>("muons", edm::InputTag("ALCARECOMuAlGlobalCosmics:StandAlone"));
0342 desc.add<edm::InputTag>("gmuons", edm::InputTag("ALCARECOMuAlGlobalCosmics:GlobalMuon"));
0343 desc.add<edm::InputTag>("smuons", edm::InputTag("ALCARECOMuAlGlobalCosmics:SelectedMuons"));
0344 desc.add<bool>("isolated", false);
0345 desc.add<bool>("cosmics", false);
0346 desc.add<bool>("refitmuon", false);
0347 desc.add<bool>("refittrack", false);
0348 desc.addUntracked<std::string>("rootOutFile", "outfile.root");
0349 desc.addUntracked<std::string>("txtOutFile", "outglobal.txt");
0350 desc.addUntracked<bool>("writeDB", false);
0351 desc.addUntracked<bool>("debug", false);
0352 descriptions.add("globalTrackerMuonAlignment", desc);
0353 }
0354
0355
0356
0357
0358
0359 void GlobalTrackerMuonAlignment::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0360
0361 GlobalTrackerMuonAlignment::analyzeTrackTrajectory(iEvent, iSetup);
0362 }
0363
0364
0365 void GlobalTrackerMuonAlignment::beginJob() {
0366 N_event = 0;
0367 N_track = 0;
0368
0369 if (cosmicMuonMode_ == true && isolatedMuonMode_ == false) {
0370 collectionCosmic = true;
0371 collectionIsolated = false;
0372 } else if (cosmicMuonMode_ == false && isolatedMuonMode_ == true) {
0373 collectionCosmic = false;
0374 collectionIsolated = true;
0375 } else {
0376 collectionCosmic = false;
0377 collectionIsolated = false;
0378 }
0379
0380 for (int i = 0; i <= 2; i++) {
0381 Gfr(i) = 0.;
0382 for (int j = 0; j <= 2; j++) {
0383 Inf(i, j) = 0.;
0384 }
0385 }
0386
0387 Grad = CLHEP::HepVector(6, 0);
0388 Hess = CLHEP::HepSymMatrix(6, 0);
0389
0390 GradL = CLHEP::HepVector(6, 0);
0391 HessL = CLHEP::HepSymMatrix(6, 0);
0392
0393
0394 TDirectory* dirsave = gDirectory;
0395
0396 file = new TFile(rootOutFile_.c_str(), "recreate");
0397 const bool oldAddDir = TH1::AddDirectoryStatus();
0398
0399 TH1::AddDirectory(true);
0400
0401 this->bookHist();
0402
0403 TH1::AddDirectory(oldAddDir);
0404 dirsave->cd();
0405 }
0406
0407
0408 void GlobalTrackerMuonAlignment::endJob() {
0409 bool alarm = false;
0410
0411 this->fitHist();
0412
0413 AlgebraicVector3 d(0., 0., 0.);
0414
0415 AlgebraicSymMatrix33 InfI;
0416 for (int i = 0; i <= 2; i++)
0417 for (int j = 0; j <= 2; j++) {
0418 if (j < i)
0419 continue;
0420 InfI(i, j) += Inf(i, j);
0421 }
0422 bool ierr = !InfI.Invert();
0423 if (ierr) {
0424 if (alarm)
0425 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Error inverse Inf matrix !!!!!!!!!!!";
0426 }
0427
0428 for (int i = 0; i <= 2; i++)
0429 for (int k = 0; k <= 2; k++)
0430 d(i) -= InfI(i, k) * Gfr(k);
0431
0432
0433
0434 CLHEP::HepVector d3 = CLHEP::solve(Hess, -Grad);
0435 int iEr3;
0436 CLHEP::HepMatrix Errd3 = Hess.inverse(iEr3);
0437 if (iEr3 != 0) {
0438 if (alarm)
0439 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " endJob Error inverse Hess matrix !!!!!!!!!!!";
0440 }
0441
0442
0443
0444 CLHEP::HepVector dLI = CLHEP::solve(HessL, -GradL);
0445 int iErI;
0446 CLHEP::HepMatrix ErrdLI = HessL.inverse(iErI);
0447 if (iErI != 0) {
0448 if (alarm)
0449 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " endJob Error inverse HessL matrix !!!!!!!!!!!";
0450 }
0451
0452
0453
0454 edm::LogVerbatim("GlobalTrackerMuonAlignment")
0455 << " ---- " << N_event << " event " << N_track << " tracks " << MuSelect << " ---- ";
0456 if (collectionIsolated)
0457 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ALCARECOMuAlCalIsolatedMu";
0458 else if (collectionCosmic)
0459 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ALCARECOMuAlGlobalCosmics";
0460 else
0461 edm::LogVerbatim("GlobalTrackerMuonAlignment") << smuonTags_;
0462 edm::LogVerbatim("GlobalTrackerMuonAlignment")
0463 << " Similated shifts[cm] " << MuGlShift(1) << " " << MuGlShift(2) << " " << MuGlShift(3) << " "
0464 << " angles[rad] " << MuGlAngle(1) << " " << MuGlAngle(2) << " " << MuGlAngle(3) << " ";
0465 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " d " << -d;
0466
0467 edm::LogVerbatim("GlobalTrackerMuonAlignment")
0468 << " +- " << std::sqrt(InfI(0, 0)) << " " << std::sqrt(InfI(1, 1)) << " " << std::sqrt(InfI(2, 2));
0469 edm::LogVerbatim("GlobalTrackerMuonAlignment")
0470 << " dG " << d3(1) << " " << d3(2) << " " << d3(3) << " " << d3(4) << " " << d3(5) << " " << d3(6);
0471
0472 edm::LogVerbatim("GlobalTrackerMuonAlignment")
0473 << " +- " << std::sqrt(Errd3(1, 1)) << " " << std::sqrt(Errd3(2, 2)) << " " << std::sqrt(Errd3(3, 3)) << " "
0474 << std::sqrt(Errd3(4, 4)) << " " << std::sqrt(Errd3(5, 5)) << " " << std::sqrt(Errd3(6, 6));
0475 edm::LogVerbatim("GlobalTrackerMuonAlignment")
0476 << " dL " << dLI(1) << " " << dLI(2) << " " << dLI(3) << " " << dLI(4) << " " << dLI(5) << " " << dLI(6);
0477
0478 edm::LogVerbatim("GlobalTrackerMuonAlignment")
0479 << " +- " << std::sqrt(ErrdLI(1, 1)) << " " << std::sqrt(ErrdLI(2, 2)) << " " << std::sqrt(ErrdLI(3, 3)) << " "
0480 << std::sqrt(ErrdLI(4, 4)) << " " << std::sqrt(ErrdLI(5, 5)) << " " << std::sqrt(ErrdLI(6, 6));
0481
0482
0483 CLHEP::HepVector vectorToDb(6, 0), vectorErrToDb(6, 0);
0484
0485
0486 vectorToDb = -dLI;
0487 for (unsigned int i = 1; i <= 6; i++)
0488 vectorErrToDb(i) = std::sqrt(ErrdLI(i, i));
0489
0490
0491 file->Write();
0492 file->Close();
0493
0494
0495 OutGlobalTxt.open(txtOutFile_.c_str(), ios::out);
0496 if (!OutGlobalTxt.is_open())
0497 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " outglobal.txt is not open !!!!!";
0498 else {
0499 OutGlobalTxt << vectorToDb(1) << " " << vectorToDb(2) << " " << vectorToDb(3) << " " << vectorToDb(4) << " "
0500 << vectorToDb(5) << " " << vectorToDb(6) << " muon Global.\n";
0501 OutGlobalTxt << vectorErrToDb(1) << " " << vectorErrToDb(1) << " " << vectorErrToDb(1) << " " << vectorErrToDb(1)
0502 << " " << vectorErrToDb(1) << " " << vectorErrToDb(1) << " errors.\n";
0503 OutGlobalTxt << N_event << " events are processed.\n";
0504
0505 if (collectionIsolated)
0506 OutGlobalTxt << "ALCARECOMuAlCalIsolatedMu.\n";
0507 else if (collectionCosmic)
0508 OutGlobalTxt << " ALCARECOMuAlGlobalCosmics.\n";
0509 else
0510 OutGlobalTxt << smuonTags_ << ".\n";
0511 OutGlobalTxt.close();
0512 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Write to the file outglobal.txt done ";
0513 }
0514
0515
0516 if (debug_)
0517 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " writeBD_ " << writeDB_;
0518 if (writeDB_)
0519 GlobalTrackerMuonAlignment::writeGlPosRcd(vectorToDb);
0520 }
0521
0522
0523 void GlobalTrackerMuonAlignment::bookHist() {
0524 double PI = 3.1415927;
0525 histo = new TH1F("Pt", "pt", 1000, 0, 100);
0526 histo2 = new TH1F("P", "P [GeV/c]", 400, 0., 400.);
0527 histo2->GetXaxis()->SetTitle("momentum [GeV/c]");
0528 histo3 = new TH1F("outerLambda", "#lambda outer", 100, -PI / 2., PI / 2.);
0529 histo3->GetXaxis()->SetTitle("#lambda outer");
0530 histo4 = new TH1F("phi", "#phi [rad]", 100, -PI, PI);
0531 histo4->GetXaxis()->SetTitle("#phi [rad]");
0532 histo5 = new TH1F("Rmuon", "inner muon hit R [cm]", 100, 0., 800.);
0533 histo5->GetXaxis()->SetTitle("R of muon [cm]");
0534 histo6 = new TH1F("Zmuon", "inner muon hit Z[cm]", 100, -1000., 1000.);
0535 histo6->GetXaxis()->SetTitle("Z of muon [cm]");
0536 histo7 = new TH1F("(Pm-Pt)/Pt", " (Pmuon-Ptrack)/Ptrack", 100, -2., 2.);
0537 histo7->GetXaxis()->SetTitle("(Pmuon-Ptrack)/Ptrack");
0538 histo8 = new TH1F("chi muon-track", "#chi^{2}(muon-track)", 1000, 0., 1000.);
0539 histo8->GetXaxis()->SetTitle("#chi^{2} of muon w.r.t. propagated track");
0540 histo11 = new TH1F("distance muon-track", "distance muon w.r.t track [cm]", 100, 0., 30.);
0541 histo11->GetXaxis()->SetTitle("distance of muon w.r.t. track [cm]");
0542 histo12 = new TH1F("Xmuon-Xtrack", "Xmuon-Xtrack [cm]", 200, -20., 20.);
0543 histo12->GetXaxis()->SetTitle("Xmuon - Xtrack [cm]");
0544 histo13 = new TH1F("Ymuon-Ytrack", "Ymuon-Ytrack [cm]", 200, -20., 20.);
0545 histo13->GetXaxis()->SetTitle("Ymuon - Ytrack [cm]");
0546 histo14 = new TH1F("Zmuon-Ztrack", "Zmuon-Ztrack [cm]", 200, -20., 20.);
0547 histo14->GetXaxis()->SetTitle("Zmuon-Ztrack [cm]");
0548 histo15 = new TH1F("NXmuon-NXtrack", "NXmuon-NXtrack [rad]", 200, -.1, .1);
0549 histo15->GetXaxis()->SetTitle("N_{X}(muon)-N_{X}(track) [rad]");
0550 histo16 = new TH1F("NYmuon-NYtrack", "NYmuon-NYtrack [rad]", 200, -.1, .1);
0551 histo16->GetXaxis()->SetTitle("N_{Y}(muon)-N_{Y}(track) [rad]");
0552 histo17 = new TH1F("NZmuon-NZtrack", "NZmuon-NZtrack [rad]", 200, -.1, .1);
0553 histo17->GetXaxis()->SetTitle("N_{Z}(muon)-N_{Z}(track) [rad]");
0554 histo18 = new TH1F("expected error of Xinner", "outer hit of inner tracker", 100, 0, .01);
0555 histo18->GetXaxis()->SetTitle("expected error of Xinner [cm]");
0556 histo19 = new TH1F("expected error of Xmuon", "inner hit of muon", 100, 0, .1);
0557 histo19->GetXaxis()->SetTitle("expected error of Xmuon [cm]");
0558 histo20 = new TH1F("expected error of Xmuon-Xtrack", "muon w.r.t. propagated track", 100, 0., 10.);
0559 histo20->GetXaxis()->SetTitle("expected error of Xmuon-Xtrack [cm]");
0560 histo21 = new TH1F("pull of Xmuon-Xtrack", "pull of Xmuon-Xtrack", 100, -10., 10.);
0561 histo21->GetXaxis()->SetTitle("(Xmuon-Xtrack)/expected error");
0562 histo22 = new TH1F("pull of Ymuon-Ytrack", "pull of Ymuon-Ytrack", 100, -10., 10.);
0563 histo22->GetXaxis()->SetTitle("(Ymuon-Ytrack)/expected error");
0564 histo23 = new TH1F("pull of Zmuon-Ztrack", "pull of Zmuon-Ztrack", 100, -10., 10.);
0565 histo23->GetXaxis()->SetTitle("(Zmuon-Ztrack)/expected error");
0566 histo24 = new TH1F("pull of PXmuon-PXtrack", "pull of PXmuon-PXtrack", 100, -10., 10.);
0567 histo24->GetXaxis()->SetTitle("(P_{X}(muon)-P_{X}(track))/expected error");
0568 histo25 = new TH1F("pull of PYmuon-PYtrack", "pull of PYmuon-PYtrack", 100, -10., 10.);
0569 histo25->GetXaxis()->SetTitle("(P_{Y}(muon)-P_{Y}(track))/expected error");
0570 histo26 = new TH1F("pull of PZmuon-PZtrack", "pull of PZmuon-PZtrack", 100, -10., 10.);
0571 histo26->GetXaxis()->SetTitle("(P_{Z}(muon)-P_{Z}(track))/expected error");
0572 histo27 = new TH1F("N_x", "Nx of tangent plane", 120, -1.1, 1.1);
0573 histo27->GetXaxis()->SetTitle("normal vector projection N_{X}");
0574 histo28 = new TH1F("N_y", "Ny of tangent plane", 120, -1.1, 1.1);
0575 histo28->GetXaxis()->SetTitle("normal vector projection N_{Y}");
0576 histo29 = new TH1F("lenght of track", "lenght of track", 200, 0., 400);
0577 histo29->GetXaxis()->SetTitle("lenght of track [cm]");
0578 histo30 = new TH1F("lenght of muon", "lenght of muon", 200, 0., 800);
0579 histo30->GetXaxis()->SetTitle("lenght of muon [cm]");
0580
0581 histo31 = new TH1F("local chi muon-track", "#local chi^{2}(muon-track)", 1000, 0., 1000.);
0582 histo31->GetXaxis()->SetTitle("#local chi^{2} of muon w.r.t. propagated track");
0583 histo32 = new TH1F("pull of Px/Pz local", "pull of Px/Pz local", 100, -10., 10.);
0584 histo32->GetXaxis()->SetTitle("local (Px/Pz(muon) - Px/Pz(track))/expected error");
0585 histo33 = new TH1F("pull of Py/Pz local", "pull of Py/Pz local", 100, -10., 10.);
0586 histo33->GetXaxis()->SetTitle("local (Py/Pz(muon) - Py/Pz(track))/expected error");
0587 histo34 = new TH1F("pull of X local", "pull of X local", 100, -10., 10.);
0588 histo34->GetXaxis()->SetTitle("local (Xmuon - Xtrack)/expected error");
0589 histo35 = new TH1F("pull of Y local", "pull of Y local", 100, -10., 10.);
0590 histo35->GetXaxis()->SetTitle("local (Ymuon - Ytrack)/expected error");
0591
0592 histo101 = new TH2F("Rtr/mu vs Ztr/mu", "hit of track/muon", 100, -800., 800., 100, 0., 600.);
0593 histo101->GetXaxis()->SetTitle("Z of track/muon [cm]");
0594 histo101->GetYaxis()->SetTitle("R of track/muon [cm]");
0595 histo102 = new TH2F("Ytr/mu vs Xtr/mu", "hit of track/muon", 100, -600., 600., 100, -600., 600.);
0596 histo102->GetXaxis()->SetTitle("X of track/muon [cm]");
0597 histo102->GetYaxis()->SetTitle("Y of track/muon [cm]");
0598 }
0599
0600
0601 void GlobalTrackerMuonAlignment::fitHist() {
0602 histo7->Fit("gaus", "Q");
0603
0604 histo12->Fit("gaus", "Q");
0605 histo13->Fit("gaus", "Q");
0606 histo14->Fit("gaus", "Q");
0607 histo15->Fit("gaus", "Q");
0608 histo16->Fit("gaus", "Q");
0609 histo17->Fit("gaus", "Q");
0610
0611 histo21->Fit("gaus", "Q");
0612 histo22->Fit("gaus", "Q");
0613 histo23->Fit("gaus", "Q");
0614 histo24->Fit("gaus", "Q");
0615 histo25->Fit("gaus", "Q");
0616 histo26->Fit("gaus", "Q");
0617
0618 histo32->Fit("gaus", "Q");
0619 histo33->Fit("gaus", "Q");
0620 histo34->Fit("gaus", "Q");
0621 histo35->Fit("gaus", "Q");
0622 }
0623
0624
0625 void GlobalTrackerMuonAlignment::analyzeTrackTrack(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0626 using namespace edm;
0627 using namespace reco;
0628
0629 bool info = false;
0630 bool alarm = false;
0631
0632 double PI = 3.1415927;
0633
0634 cosmicMuonMode_ = true;
0635 isolatedMuonMode_ = !cosmicMuonMode_;
0636
0637
0638
0639
0640 N_event++;
0641 if (info || debug_) {
0642 edm::LogVerbatim("GlobalTrackerMuonAlignment") << "----- event " << N_event << " -- tracks " << N_track << " ---";
0643 if (cosmicMuonMode_)
0644 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " treated as CosmicMu ";
0645 if (isolatedMuonMode_)
0646 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " treated as IsolatedMu ";
0647 edm::LogVerbatim("GlobalTrackerMuonAlignment");
0648 }
0649
0650 const edm::Handle<reco::TrackCollection>& tracks =
0651 ((collectionIsolated)
0652 ? iEvent.getHandle(trackIsolateToken_)
0653 : ((collectionCosmic) ? iEvent.getHandle(trackCosmicToken_) : iEvent.getHandle(trackToken_)));
0654 const edm::Handle<reco::TrackCollection>& muons =
0655 ((collectionIsolated) ? iEvent.getHandle(muonIsolateToken_)
0656 : ((collectionCosmic) ? iEvent.getHandle(muonCosmicToken_) : iEvent.getHandle(muonToken_)));
0657 const edm::Handle<reco::TrackCollection>& gmuons =
0658 ((collectionIsolated)
0659 ? iEvent.getHandle(gmuonIsolateToken_)
0660 : ((collectionCosmic) ? iEvent.getHandle(gmuonCosmicToken_) : iEvent.getHandle(gmuonToken_)));
0661 const edm::Handle<reco::MuonCollection>& smuons =
0662 ((collectionIsolated)
0663 ? iEvent.getHandle(smuonIsolateToken_)
0664 : ((collectionCosmic) ? iEvent.getHandle(smuonCosmicToken_) : iEvent.getHandle(smuonToken_)));
0665
0666 if (debug_) {
0667 edm::LogVerbatim("GlobalTrackerMuonAlignment")
0668 << " ievBunch " << iEvent.bunchCrossing() << " runN " << (int)iEvent.run();
0669 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " N tracks s/amu gmu selmu " << tracks->size() << " "
0670 << muons->size() << " " << gmuons->size() << " " << smuons->size();
0671 for (MuonCollection::const_iterator itMuon = smuons->begin(); itMuon != smuons->end(); ++itMuon) {
0672 edm::LogVerbatim("GlobalTrackerMuonAlignment")
0673 << " is isolatValid Matches " << itMuon->isIsolationValid() << " " << itMuon->isMatchesValid();
0674 }
0675 }
0676
0677 if (isolatedMuonMode_) {
0678 if (tracks->size() != 1)
0679 return;
0680 if (muons->size() != 1)
0681 return;
0682 if (gmuons->size() != 1)
0683 return;
0684 if (smuons->size() != 1)
0685 return;
0686 }
0687
0688 if (cosmicMuonMode_) {
0689 if (smuons->size() > 2)
0690 return;
0691 if (tracks->size() != smuons->size())
0692 return;
0693 if (muons->size() != smuons->size())
0694 return;
0695 if (gmuons->size() != smuons->size())
0696 return;
0697 }
0698
0699
0700
0701
0702
0703
0704
0705
0706
0707
0708
0709 if (watchTrackingGeometry_.check(iSetup) || !trackingGeometry_) {
0710 edm::ESHandle<GlobalTrackingGeometry> trackingGeometry = iSetup.getHandle(m_TkGeometryToken);
0711 trackingGeometry_ = &*trackingGeometry;
0712 }
0713
0714 if (watchMagneticFieldRecord_.check(iSetup) || !magneticField_) {
0715 edm::ESHandle<MagneticField> magneticField = iSetup.getHandle(m_MagFieldToken);
0716 magneticField_ = &*magneticField;
0717 }
0718
0719 if (watchGlobalPositionRcd_.check(iSetup) || !globalPositionRcd_) {
0720 edm::ESHandle<Alignments> globalPositionRcd = iSetup.getHandle(m_globalPosToken);
0721 globalPositionRcd_ = &*globalPositionRcd;
0722 for (std::vector<AlignTransform>::const_iterator i = globalPositionRcd_->m_align.begin();
0723 i != globalPositionRcd_->m_align.end();
0724 ++i) {
0725 if (DetId(DetId::Tracker).rawId() == i->rawId())
0726 iteratorTrackerRcd = i;
0727 if (DetId(DetId::Muon).rawId() == i->rawId())
0728 iteratorMuonRcd = i;
0729 if (DetId(DetId::Ecal).rawId() == i->rawId())
0730 iteratorEcalRcd = i;
0731 if (DetId(DetId::Hcal).rawId() == i->rawId())
0732 iteratorHcalRcd = i;
0733 }
0734 if (true || debug_) {
0735 edm::LogVerbatim("GlobalTrackerMuonAlignment") << "=== iteratorMuonRcd " << iteratorMuonRcd->rawId() << " ====\n"
0736 << " translation " << iteratorMuonRcd->translation() << "\n"
0737 << " angles " << iteratorMuonRcd->rotation().eulerAngles();
0738 edm::LogVerbatim("GlobalTrackerMuonAlignment") << iteratorMuonRcd->rotation();
0739 }
0740 }
0741
0742 edm::ESHandle<Propagator> propagator = iSetup.getHandle(m_propToken);
0743
0744 SteppingHelixPropagator alongStHePr = SteppingHelixPropagator(magneticField_, alongMomentum);
0745 SteppingHelixPropagator oppositeStHePr = SteppingHelixPropagator(magneticField_, oppositeToMomentum);
0746
0747
0748 defaultRKPropagator::Product aprod(magneticField_, alongMomentum, 5.e-5);
0749 auto& alongRKPr = aprod.propagator;
0750 defaultRKPropagator::Product oprod(magneticField_, oppositeToMomentum, 5.e-5);
0751 auto& oppositeRKPr = oprod.propagator;
0752
0753 float epsilon = 5.;
0754 SmartPropagator alongSmPr = SmartPropagator(alongRKPr, alongStHePr, magneticField_, alongMomentum, epsilon);
0755 SmartPropagator oppositeSmPr =
0756 SmartPropagator(oppositeRKPr, oppositeStHePr, magneticField_, oppositeToMomentum, epsilon);
0757
0758
0759
0760 for (MuonCollection::const_iterator itMuon = smuons->begin(); itMuon != smuons->end(); ++itMuon) {
0761 if (debug_) {
0762 edm::LogVerbatim("GlobalTrackerMuonAlignment")
0763 << " mu gM is GM Mu SaM tM " << itMuon->isGlobalMuon() << " " << itMuon->isMuon() << " "
0764 << itMuon->isStandAloneMuon() << " " << itMuon->isTrackerMuon() << " ";
0765 }
0766
0767
0768 TransientTrack muTT(itMuon->outerTrack(), magneticField_, trackingGeometry_);
0769 TrajectoryStateOnSurface innerMuTSOS = muTT.innermostMeasurementState();
0770 TrajectoryStateOnSurface outerMuTSOS = muTT.outermostMeasurementState();
0771
0772
0773 TransientTrack trackTT(itMuon->track(), magneticField_, trackingGeometry_);
0774 TrajectoryStateOnSurface outerTrackTSOS = trackTT.outermostMeasurementState();
0775 TrajectoryStateOnSurface innerTrackTSOS = trackTT.innermostMeasurementState();
0776
0777 GlobalPoint pointTrackIn = innerTrackTSOS.globalPosition();
0778 GlobalPoint pointTrackOut = outerTrackTSOS.globalPosition();
0779 float lenghtTrack = (pointTrackOut - pointTrackIn).mag();
0780 GlobalPoint pointMuonIn = innerMuTSOS.globalPosition();
0781 GlobalPoint pointMuonOut = outerMuTSOS.globalPosition();
0782 float lenghtMuon = (pointMuonOut - pointMuonIn).mag();
0783 GlobalVector momentumTrackOut = outerTrackTSOS.globalMomentum();
0784 GlobalVector momentumTrackIn = innerTrackTSOS.globalMomentum();
0785 float distanceInIn = (pointTrackIn - pointMuonIn).mag();
0786 float distanceInOut = (pointTrackIn - pointMuonOut).mag();
0787 float distanceOutIn = (pointTrackOut - pointMuonIn).mag();
0788 float distanceOutOut = (pointTrackOut - pointMuonOut).mag();
0789
0790 if (debug_) {
0791 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " pointTrackIn " << pointTrackIn;
0792 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Out " << pointTrackOut << " lenght " << lenghtTrack;
0793 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " point MuonIn " << pointMuonIn;
0794 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Out " << pointMuonOut << " lenght " << lenghtMuon;
0795 edm::LogVerbatim("GlobalTrackerMuonAlignment")
0796 << " momeTrackIn Out " << momentumTrackIn << " " << momentumTrackOut;
0797 edm::LogVerbatim("GlobalTrackerMuonAlignment")
0798 << " doi io ii oo " << distanceOutIn << " " << distanceInOut << " " << distanceInIn << " " << distanceOutOut;
0799 }
0800
0801 if (lenghtTrack < 90.)
0802 continue;
0803 if (lenghtMuon < 300.)
0804 continue;
0805 if (momentumTrackIn.mag() < 15. || momentumTrackOut.mag() < 15.)
0806 continue;
0807 if (trackTT.charge() != muTT.charge())
0808 continue;
0809
0810 if (debug_)
0811 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " passed lenght momentum cuts";
0812
0813 GlobalVector GRm, GPm, Nl, Rm, Pm, Rt, Pt, Rt0;
0814 AlgebraicSymMatrix66 Cm, C0, Ce, C1;
0815
0816 TrajectoryStateOnSurface extrapolationT;
0817 int tsosMuonIf = 0;
0818
0819 if (isolatedMuonMode_) {
0820 const Surface& refSurface = innerMuTSOS.surface();
0821 ConstReferenceCountingPointer<TangentPlane> tpMuLocal(refSurface.tangentPlane(innerMuTSOS.localPosition()));
0822 Nl = tpMuLocal->normalVector();
0823 ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
0824 GlobalVector Ng = tpMuGlobal->normalVector();
0825 Surface* surf = (Surface*)&refSurface;
0826 const Plane* refPlane = dynamic_cast<Plane*>(surf);
0827 GlobalVector Nlp = refPlane->normalVector();
0828
0829 if (debug_) {
0830 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
0831 << " alfa " << int(asin(Nl.x()) * 57.29578);
0832 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " global " << Ng.x() << " " << Ng.y() << " " << Ng.z();
0833 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " lp " << Nlp.x() << " " << Nlp.y() << " " << Nlp.z();
0834 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug")
0835 << " Nlocal Nglobal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " " << Ng.x() << " " << Ng.y() << " "
0836 << Ng.z() << " alfa " << static_cast<int>(asin(Nl.x()) * 57.29578);
0837 }
0838
0839
0840 extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
0841
0842
0843 if (!extrapolationT.isValid()) {
0844 if (false & alarm)
0845 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
0846
0847 ;
0848 continue;
0849 }
0850 tsosMuonIf = 1;
0851 Rt = GlobalVector((extrapolationT.globalPosition()).x(),
0852 (extrapolationT.globalPosition()).y(),
0853 (extrapolationT.globalPosition()).z());
0854
0855 Pt = extrapolationT.globalMomentum();
0856
0857 GRm = GlobalVector(
0858 (innerMuTSOS.globalPosition()).x(), (innerMuTSOS.globalPosition()).y(), (innerMuTSOS.globalPosition()).z());
0859 GPm = innerMuTSOS.globalMomentum();
0860 Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
0861 (outerTrackTSOS.globalPosition()).y(),
0862 (outerTrackTSOS.globalPosition()).z());
0863 Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + innerMuTSOS.cartesianError().matrix());
0864 C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
0865 Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
0866 C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
0867
0868 }
0869
0870 if (cosmicMuonMode_) {
0871
0872 if ((distanceOutIn <= distanceInOut) & (distanceOutIn <= distanceInIn) &
0873 (distanceOutIn <= distanceOutOut)) {
0874 if (debug_)
0875 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ----- Out - In -----";
0876
0877 const Surface& refSurface = innerMuTSOS.surface();
0878 ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
0879 Nl = tpMuGlobal->normalVector();
0880
0881
0882 extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
0883
0884
0885 if (!extrapolationT.isValid()) {
0886 if (false & alarm)
0887 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
0888
0889 ;
0890 continue;
0891 }
0892 if (debug_)
0893 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " extrapolationT.isValid " << extrapolationT.isValid();
0894
0895 tsosMuonIf = 1;
0896 Rt = GlobalVector((extrapolationT.globalPosition()).x(),
0897 (extrapolationT.globalPosition()).y(),
0898 (extrapolationT.globalPosition()).z());
0899
0900 Pt = extrapolationT.globalMomentum();
0901
0902 GRm = GlobalVector(
0903 (innerMuTSOS.globalPosition()).x(), (innerMuTSOS.globalPosition()).y(), (innerMuTSOS.globalPosition()).z());
0904 GPm = innerMuTSOS.globalMomentum();
0905 Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
0906 (outerTrackTSOS.globalPosition()).y(),
0907 (outerTrackTSOS.globalPosition()).z());
0908 Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + innerMuTSOS.cartesianError().matrix());
0909 C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
0910 Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
0911 C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
0912
0913 if (false & debug_) {
0914 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ->propDir " << propagator->propagationDirection();
0915 PropagationDirectionChooser Chooser;
0916 edm::LogVerbatim("GlobalTrackerMuonAlignment")
0917 << " propDirCh " << Chooser.operator()(*outerTrackTSOS.freeState(), refSurface) << " Ch == along "
0918 << (alongMomentum == Chooser.operator()(*outerTrackTSOS.freeState(), refSurface));
0919 edm::LogVerbatim("GlobalTrackerMuonAlignment")
0920 << " --- Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
0921 << " alfa " << int(asin(Nl.x()) * 57.29578);
0922 }
0923 }
0924
0925 else if ((distanceInOut <= distanceInIn) & (distanceInOut <= distanceOutIn) &
0926 (distanceInOut <= distanceOutOut)) {
0927 if (debug_)
0928 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ----- In - Out -----";
0929
0930 const Surface& refSurface = outerMuTSOS.surface();
0931 ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
0932 Nl = tpMuGlobal->normalVector();
0933
0934
0935 extrapolationT = oppositeSmPr.propagate(innerTrackTSOS, refSurface);
0936
0937 if (!extrapolationT.isValid()) {
0938 if (false & alarm)
0939 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
0940 << "\a\a\a\a\a\a\a\a" << extrapolationT.isValid();
0941 continue;
0942 }
0943 if (debug_)
0944 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " extrapolationT.isValid " << extrapolationT.isValid();
0945
0946 tsosMuonIf = 2;
0947 Rt = GlobalVector((extrapolationT.globalPosition()).x(),
0948 (extrapolationT.globalPosition()).y(),
0949 (extrapolationT.globalPosition()).z());
0950
0951 Pt = extrapolationT.globalMomentum();
0952
0953 GRm = GlobalVector(
0954 (outerMuTSOS.globalPosition()).x(), (outerMuTSOS.globalPosition()).y(), (outerMuTSOS.globalPosition()).z());
0955 GPm = outerMuTSOS.globalMomentum();
0956 Rt0 = GlobalVector((innerTrackTSOS.globalPosition()).x(),
0957 (innerTrackTSOS.globalPosition()).y(),
0958 (innerTrackTSOS.globalPosition()).z());
0959 Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + outerMuTSOS.cartesianError().matrix());
0960 C0 = AlgebraicSymMatrix66(innerTrackTSOS.cartesianError().matrix());
0961 Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
0962 C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
0963
0964 if (false & debug_) {
0965 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ->propDir " << propagator->propagationDirection();
0966 PropagationDirectionChooser Chooser;
0967 edm::LogVerbatim("GlobalTrackerMuonAlignment")
0968 << " propDirCh " << Chooser.operator()(*innerTrackTSOS.freeState(), refSurface) << " Ch == oppisite "
0969 << (oppositeToMomentum == Chooser.operator()(*innerTrackTSOS.freeState(), refSurface));
0970 edm::LogVerbatim("GlobalTrackerMuonAlignment")
0971 << " --- Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
0972 << " alfa " << int(asin(Nl.x()) * 57.29578);
0973 }
0974 }
0975
0976 else if ((distanceOutOut <= distanceInOut) & (distanceOutOut <= distanceInIn) &
0977 (distanceOutOut <= distanceOutIn)) {
0978 if (debug_)
0979 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ----- Out - Out -----";
0980
0981
0982 continue;
0983
0984 const Surface& refSurface = outerMuTSOS.surface();
0985 ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
0986 Nl = tpMuGlobal->normalVector();
0987
0988
0989 extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
0990
0991 if (!extrapolationT.isValid()) {
0992 if (alarm)
0993 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " !!!!!!!!! Catastrophe Out-Out extrapolationT.isValid "
0994 << "\a\a\a\a\a\a\a\a" << extrapolationT.isValid();
0995 continue;
0996 }
0997 if (debug_)
0998 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " extrapolationT.isValid " << extrapolationT.isValid();
0999
1000 tsosMuonIf = 2;
1001 Rt = GlobalVector((extrapolationT.globalPosition()).x(),
1002 (extrapolationT.globalPosition()).y(),
1003 (extrapolationT.globalPosition()).z());
1004
1005 Pt = extrapolationT.globalMomentum();
1006
1007 GRm = GlobalVector(
1008 (outerMuTSOS.globalPosition()).x(), (outerMuTSOS.globalPosition()).y(), (outerMuTSOS.globalPosition()).z());
1009 GPm = outerMuTSOS.globalMomentum();
1010 Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
1011 (outerTrackTSOS.globalPosition()).y(),
1012 (outerTrackTSOS.globalPosition()).z());
1013 Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + outerMuTSOS.cartesianError().matrix());
1014 C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
1015 Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
1016 C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
1017
1018 if (debug_) {
1019 PropagationDirectionChooser Chooser;
1020 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1021 << " propDirCh " << Chooser.operator()(*outerTrackTSOS.freeState(), refSurface) << " Ch == along "
1022 << (alongMomentum == Chooser.operator()(*outerTrackTSOS.freeState(), refSurface));
1023 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1024 << " --- Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1025 << " alfa " << int(asin(Nl.x()) * 57.29578);
1026 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1027 << " Nornal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1028 << " alfa " << int(asin(Nl.x()) * 57.29578);
1029 }
1030 }
1031 else {
1032 if (alarm)
1033 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1034 << " ------- !!!!!!!!!!!!!!! No proper Out - In \a\a\a\a\a\a\a";
1035 continue;
1036 }
1037
1038 }
1039
1040 if (tsosMuonIf == 0) {
1041 if (info) {
1042 edm::LogVerbatim("GlobalTrackerMuonAlignment") << "No tsosMuon !!!!!!";
1043 continue;
1044 }
1045 }
1046 TrajectoryStateOnSurface tsosMuon;
1047 if (tsosMuonIf == 1)
1048 tsosMuon = muTT.innermostMeasurementState();
1049 else
1050 tsosMuon = muTT.outermostMeasurementState();
1051
1052
1053 AlgebraicVector4 LPRm;
1054 GlobalTrackerMuonAlignment::misalignMuonL(GRm, GPm, Nl, Rt, Rm, Pm, LPRm, extrapolationT, tsosMuon);
1055
1056 GlobalVector resR = Rm - Rt;
1057 GlobalVector resP0 = Pm - Pt;
1058 GlobalVector resP = Pm / Pm.mag() - Pt / Pt.mag();
1059 float RelMomResidual = (Pm.mag() - Pt.mag()) / (Pt.mag() + 1.e-6);
1060 ;
1061
1062 AlgebraicVector6 Vm;
1063 Vm(0) = resR.x();
1064 Vm(1) = resR.y();
1065 Vm(2) = resR.z();
1066 Vm(3) = resP0.x();
1067 Vm(4) = resP0.y();
1068 Vm(5) = resP0.z();
1069 float Rmuon = Rm.perp();
1070 float Zmuon = Rm.z();
1071 float alfa_x = atan2(Nl.x(), Nl.y()) * 57.29578;
1072
1073 if (debug_) {
1074 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1075 << " Nx Ny Nz alfa_x " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " " << alfa_x;
1076 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug")
1077 << " Rm " << Rm << "\n Rt " << Rt << "\n resR " << resR << "\n resP " << resP << " dp/p " << RelMomResidual;
1078 }
1079
1080 double chi_d = 0;
1081 for (int i = 0; i <= 5; i++)
1082 chi_d += Vm(i) * Vm(i) / Cm(i, i);
1083
1084 AlgebraicVector5 Vml(tsosMuon.localParameters().vector() - extrapolationT.localParameters().vector());
1085 AlgebraicSymMatrix55 m(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
1086 AlgebraicSymMatrix55 Cml(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
1087 bool ierrLoc = !m.Invert();
1088 if (ierrLoc && debug_ && info) {
1089 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ==== Error inverting Local covariance matrix ==== ";
1090 continue;
1091 }
1092 double chi_Loc = ROOT::Math::Similarity(Vml, m);
1093 if (debug_)
1094 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1095 << " chi_Loc px/pz/err " << chi_Loc << " " << Vml(1) / std::sqrt(Cml(1, 1));
1096
1097 if (Pt.mag() < 15.)
1098 continue;
1099 if (Pm.mag() < 5.)
1100 continue;
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144 if (Rmuon < 120. || Rmuon > 450.)
1145 continue;
1146 if (Zmuon < -720. || Zmuon > 720.)
1147 continue;
1148 MuSelect = " Barrel+EndCaps";
1149
1150 if (debug_)
1151 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " .............. passed all cuts";
1152
1153 N_track++;
1154
1155
1156 GlobalTrackerMuonAlignment::gradientGlobalAlg(Rt, Pt, Rm, Nl, Cm);
1157 GlobalTrackerMuonAlignment::gradientGlobal(Rt, Pt, Rm, Pm, Nl, Cm);
1158
1159 CLHEP::HepSymMatrix covLoc(4, 0);
1160 for (int i = 1; i <= 4; i++)
1161 for (int j = 1; j <= i; j++) {
1162 covLoc(i, j) = (tsosMuon.localError().matrix() + extrapolationT.localError().matrix())(i, j);
1163
1164 }
1165
1166 const Surface& refSurface = tsosMuon.surface();
1167 CLHEP::HepMatrix rotLoc(3, 3, 0);
1168 rotLoc(1, 1) = refSurface.rotation().xx();
1169 rotLoc(1, 2) = refSurface.rotation().xy();
1170 rotLoc(1, 3) = refSurface.rotation().xz();
1171
1172 rotLoc(2, 1) = refSurface.rotation().yx();
1173 rotLoc(2, 2) = refSurface.rotation().yy();
1174 rotLoc(2, 3) = refSurface.rotation().yz();
1175
1176 rotLoc(3, 1) = refSurface.rotation().zx();
1177 rotLoc(3, 2) = refSurface.rotation().zy();
1178 rotLoc(3, 3) = refSurface.rotation().zz();
1179
1180 CLHEP::HepVector posLoc(3);
1181 posLoc(1) = refSurface.position().x();
1182 posLoc(2) = refSurface.position().y();
1183 posLoc(3) = refSurface.position().z();
1184
1185 GlobalTrackerMuonAlignment::gradientLocal(Rt, Pt, Rm, Pm, Nl, covLoc, rotLoc, posLoc, LPRm);
1186
1187 if (debug_) {
1188 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Norm " << Nl;
1189 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " posLoc " << posLoc.T();
1190 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " rotLoc " << rotLoc;
1191 }
1192
1193
1194 histo->Fill(itMuon->track()->pt());
1195
1196
1197 histo2->Fill(Pt.mag());
1198 histo3->Fill((PI / 2. - itMuon->track()->outerTheta()));
1199 histo4->Fill(itMuon->track()->phi());
1200 histo5->Fill(Rmuon);
1201 histo6->Fill(Zmuon);
1202 histo7->Fill(RelMomResidual);
1203
1204 histo8->Fill(chi_d);
1205
1206 histo101->Fill(Zmuon, Rmuon);
1207 histo101->Fill(Rt0.z(), Rt0.perp());
1208 histo102->Fill(Rt0.x(), Rt0.y());
1209 histo102->Fill(Rm.x(), Rm.y());
1210
1211 histo11->Fill(resR.mag());
1212 if (fabs(Nl.x()) < 0.98)
1213 histo12->Fill(resR.x());
1214 if (fabs(Nl.y()) < 0.98)
1215 histo13->Fill(resR.y());
1216 if (fabs(Nl.z()) < 0.98)
1217 histo14->Fill(resR.z());
1218 histo15->Fill(resP.x());
1219 histo16->Fill(resP.y());
1220 histo17->Fill(resP.z());
1221
1222 if ((fabs(Nl.x()) < 0.98) && (fabs(Nl.y()) < 0.98) && (fabs(Nl.z()) < 0.98)) {
1223 histo18->Fill(std::sqrt(C0(0, 0)));
1224 histo19->Fill(std::sqrt(C1(0, 0)));
1225 histo20->Fill(std::sqrt(C1(0, 0) + Ce(0, 0)));
1226 }
1227 if (fabs(Nl.x()) < 0.98)
1228 histo21->Fill(Vm(0) / std::sqrt(Cm(0, 0)));
1229 if (fabs(Nl.y()) < 0.98)
1230 histo22->Fill(Vm(1) / std::sqrt(Cm(1, 1)));
1231 if (fabs(Nl.z()) < 0.98)
1232 histo23->Fill(Vm(2) / std::sqrt(Cm(2, 2)));
1233 histo24->Fill(Vm(3) / std::sqrt(C1(3, 3) + Ce(3, 3)));
1234 histo25->Fill(Vm(4) / std::sqrt(C1(4, 4) + Ce(4, 4)));
1235 histo26->Fill(Vm(5) / std::sqrt(C1(5, 5) + Ce(5, 5)));
1236 histo27->Fill(Nl.x());
1237 histo28->Fill(Nl.y());
1238 histo29->Fill(lenghtTrack);
1239 histo30->Fill(lenghtMuon);
1240 histo31->Fill(chi_Loc);
1241 histo32->Fill(Vml(1) / std::sqrt(Cml(1, 1)));
1242 histo33->Fill(Vml(2) / std::sqrt(Cml(2, 2)));
1243 histo34->Fill(Vml(3) / std::sqrt(Cml(3, 3)));
1244 histo35->Fill(Vml(4) / std::sqrt(Cml(4, 4)));
1245
1246 if (debug_) {
1247
1248 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " diag 0[ " << C0(0, 0) << " " << C0(1, 1) << " " << C0(2, 2)
1249 << " " << C0(3, 3) << " " << C0(4, 4) << " " << C0(5, 5) << " ]";
1250 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " diag e[ " << Ce(0, 0) << " " << Ce(1, 1) << " " << Ce(2, 2)
1251 << " " << Ce(3, 3) << " " << Ce(4, 4) << " " << Ce(5, 5) << " ]";
1252 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " diag 1[ " << C1(0, 0) << " " << C1(1, 1) << " " << C1(2, 2)
1253 << " " << C1(3, 3) << " " << C1(4, 4) << " " << C1(5, 5) << " ]";
1254 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Rm " << Rm.x() << " " << Rm.y() << " " << Rm.z() << " Pm "
1255 << Pm.x() << " " << Pm.y() << " " << Pm.z();
1256 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Rt " << Rt.x() << " " << Rt.y() << " " << Rt.z() << " Pt "
1257 << Pt.x() << " " << Pt.y() << " " << Pt.z();
1258 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Nl*(Rm-Rt) " << Nl.dot(Rm - Rt);
1259 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " resR " << resR.x() << " " << resR.y() << " " << resR.z()
1260 << " resP " << resP.x() << " " << resP.y() << " " << resP.z();
1261 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1262 << " Rm-t " << (Rm - Rt).x() << " " << (Rm - Rt).y() << " " << (Rm - Rt).z() << " Pm-t " << (Pm - Pt).x()
1263 << " " << (Pm - Pt).y() << " " << (Pm - Pt).z();
1264 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Vm " << Vm;
1265 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1266 << " +- " << std::sqrt(Cm(0, 0)) << " " << std::sqrt(Cm(1, 1)) << " " << std::sqrt(Cm(2, 2)) << " "
1267 << std::sqrt(Cm(3, 3)) << " " << std::sqrt(Cm(4, 4)) << " " << std::sqrt(Cm(5, 5));
1268 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1269 << " Pmuon Ptrack dP/Ptrack " << itMuon->outerTrack()->p() << " " << itMuon->track()->outerP() << " "
1270 << (itMuon->outerTrack()->p() - itMuon->track()->outerP()) / itMuon->track()->outerP();
1271 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " cov matrix ";
1272 edm::LogVerbatim("GlobalTrackerMuonAlignment") << Cm;
1273 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " diag [ " << Cm(0, 0) << " " << Cm(1, 1) << " " << Cm(2, 2)
1274 << " " << Cm(3, 3) << " " << Cm(4, 4) << " " << Cm(5, 5) << " ]";
1275
1276 AlgebraicSymMatrix66 Ro;
1277 double Diag[6];
1278 for (int i = 0; i <= 5; i++)
1279 Diag[i] = std::sqrt(Cm(i, i));
1280 for (int i = 0; i <= 5; i++)
1281 for (int j = 0; j <= 5; j++)
1282 Ro(i, j) = Cm(i, j) / Diag[i] / Diag[j];
1283 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " correlation matrix ";
1284 edm::LogVerbatim("GlobalTrackerMuonAlignment") << Ro;
1285
1286 AlgebraicSymMatrix66 CmI;
1287 for (int i = 0; i <= 5; i++)
1288 for (int j = 0; j <= 5; j++)
1289 CmI(i, j) = Cm(i, j);
1290
1291 bool ierr = !CmI.Invert();
1292 if (ierr) {
1293 if (alarm || debug_)
1294 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Error inverse covariance matrix !!!!!!!!!!!";
1295 continue;
1296 }
1297 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " inverse cov matrix ";
1298 edm::LogVerbatim("GlobalTrackerMuonAlignment") << Cm;
1299
1300 double chi = ROOT::Math::Similarity(Vm, CmI);
1301 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " chi chi_d " << chi << " " << chi_d;
1302 }
1303
1304 }
1305
1306 }
1307
1308
1309 void GlobalTrackerMuonAlignment::analyzeTrackTrajectory(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
1310 using namespace edm;
1311 using namespace reco;
1312
1313 bool info = false;
1314 bool alarm = false;
1315
1316 double PI = 3.1415927;
1317
1318
1319 cosmicMuonMode_ = true;
1320
1321
1322 isolatedMuonMode_ = !cosmicMuonMode_;
1323
1324 N_event++;
1325 if (info || debug_) {
1326 edm::LogVerbatim("GlobalTrackerMuonAlignment") << "----- event " << N_event << " -- tracks " << N_track << " ---";
1327 if (cosmicMuonMode_)
1328 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " treated as CosmicMu ";
1329 if (isolatedMuonMode_)
1330 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " treated as IsolatedMu ";
1331 edm::LogVerbatim("GlobalTrackerMuonAlignment");
1332 }
1333
1334 const edm::Handle<reco::TrackCollection>& tracks =
1335 ((collectionIsolated)
1336 ? iEvent.getHandle(trackIsolateToken_)
1337 : ((collectionCosmic) ? iEvent.getHandle(trackCosmicToken_) : iEvent.getHandle(trackToken_)));
1338 const edm::Handle<reco::TrackCollection>& muons =
1339 ((collectionIsolated) ? iEvent.getHandle(muonIsolateToken_)
1340 : ((collectionCosmic) ? iEvent.getHandle(muonCosmicToken_) : iEvent.getHandle(muonToken_)));
1341 const edm::Handle<reco::TrackCollection>& gmuons =
1342 ((collectionIsolated)
1343 ? iEvent.getHandle(gmuonIsolateToken_)
1344 : ((collectionCosmic) ? iEvent.getHandle(gmuonCosmicToken_) : iEvent.getHandle(gmuonToken_)));
1345 const edm::Handle<reco::MuonCollection>& smuons =
1346 ((collectionIsolated)
1347 ? iEvent.getHandle(smuonIsolateToken_)
1348 : ((collectionCosmic) ? iEvent.getHandle(smuonCosmicToken_) : iEvent.getHandle(smuonToken_)));
1349
1350 if (debug_) {
1351 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1352 << " ievBunch " << iEvent.bunchCrossing() << " runN " << (int)iEvent.run();
1353 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " N tracks s/amu gmu selmu " << tracks->size() << " "
1354 << muons->size() << " " << gmuons->size() << " " << smuons->size();
1355 for (MuonCollection::const_iterator itMuon = smuons->begin(); itMuon != smuons->end(); ++itMuon) {
1356 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1357 << " is isolatValid Matches " << itMuon->isIsolationValid() << " " << itMuon->isMatchesValid();
1358 }
1359 }
1360
1361 if (isolatedMuonMode_) {
1362 if (tracks->size() != 1)
1363 return;
1364 if (muons->size() != 1)
1365 return;
1366 if (gmuons->size() != 1)
1367 return;
1368 if (smuons->size() != 1)
1369 return;
1370 }
1371
1372 if (cosmicMuonMode_) {
1373 if (smuons->size() > 2)
1374 return;
1375 if (tracks->size() != smuons->size())
1376 return;
1377 if (muons->size() != smuons->size())
1378 return;
1379 if (gmuons->size() != smuons->size())
1380 return;
1381 }
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393 if (watchTrackingGeometry_.check(iSetup) || !trackingGeometry_) {
1394 edm::ESHandle<GlobalTrackingGeometry> trackingGeometry = iSetup.getHandle(m_TkGeometryToken);
1395 trackingGeometry_ = &*trackingGeometry;
1396 theTrackingGeometry = trackingGeometry;
1397 }
1398
1399 if (watchMagneticFieldRecord_.check(iSetup) || !magneticField_) {
1400 edm::ESHandle<MagneticField> magneticField = iSetup.getHandle(m_MagFieldToken);
1401 magneticField_ = &*magneticField;
1402 }
1403
1404 if (watchGlobalPositionRcd_.check(iSetup) || !globalPositionRcd_) {
1405 edm::ESHandle<Alignments> globalPositionRcd = iSetup.getHandle(m_globalPosToken);
1406 globalPositionRcd_ = &*globalPositionRcd;
1407 for (std::vector<AlignTransform>::const_iterator i = globalPositionRcd_->m_align.begin();
1408 i != globalPositionRcd_->m_align.end();
1409 ++i) {
1410 if (DetId(DetId::Tracker).rawId() == i->rawId())
1411 iteratorTrackerRcd = i;
1412 if (DetId(DetId::Muon).rawId() == i->rawId())
1413 iteratorMuonRcd = i;
1414 if (DetId(DetId::Ecal).rawId() == i->rawId())
1415 iteratorEcalRcd = i;
1416 if (DetId(DetId::Hcal).rawId() == i->rawId())
1417 iteratorHcalRcd = i;
1418 }
1419 if (debug_) {
1420 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1421 << "=== iteratorTrackerRcd " << iteratorTrackerRcd->rawId() << " ====\n"
1422 << " translation " << iteratorTrackerRcd->translation() << "\n"
1423 << " angles " << iteratorTrackerRcd->rotation().eulerAngles();
1424 edm::LogVerbatim("GlobalTrackerMuonAlignment") << iteratorTrackerRcd->rotation();
1425 edm::LogVerbatim("GlobalTrackerMuonAlignment") << "=== iteratorMuonRcd " << iteratorMuonRcd->rawId() << " ====\n"
1426 << " translation " << iteratorMuonRcd->translation() << "\n"
1427 << " angles " << iteratorMuonRcd->rotation().eulerAngles();
1428 edm::LogVerbatim("GlobalTrackerMuonAlignment") << iteratorMuonRcd->rotation();
1429 }
1430 }
1431
1432 edm::ESHandle<Propagator> propagator = iSetup.getHandle(m_propToken);
1433
1434 SteppingHelixPropagator alongStHePr = SteppingHelixPropagator(magneticField_, alongMomentum);
1435 SteppingHelixPropagator oppositeStHePr = SteppingHelixPropagator(magneticField_, oppositeToMomentum);
1436
1437
1438 defaultRKPropagator::Product aprod(magneticField_, alongMomentum, 5.e-5);
1439 auto& alongRKPr = aprod.propagator;
1440 defaultRKPropagator::Product oprod(magneticField_, oppositeToMomentum, 5.e-5);
1441 auto& oppositeRKPr = oprod.propagator;
1442
1443 float epsilon = 5.;
1444 SmartPropagator alongSmPr = SmartPropagator(alongRKPr, alongStHePr, magneticField_, alongMomentum, epsilon);
1445 SmartPropagator oppositeSmPr =
1446 SmartPropagator(oppositeRKPr, oppositeStHePr, magneticField_, oppositeToMomentum, epsilon);
1447
1448 if (defineFitter) {
1449 if (debug_)
1450 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ............... DEFINE FITTER ...................";
1451 KFUpdator* theUpdator = new KFUpdator();
1452
1453 Chi2MeasurementEstimator* theEstimator = new Chi2MeasurementEstimator(100000, 100000);
1454 theFitter = new KFTrajectoryFitter(alongSmPr, *theUpdator, *theEstimator);
1455 theSmoother = new KFTrajectorySmoother(alongSmPr, *theUpdator, *theEstimator);
1456 theFitterOp = new KFTrajectoryFitter(oppositeSmPr, *theUpdator, *theEstimator);
1457 theSmootherOp = new KFTrajectorySmoother(oppositeSmPr, *theUpdator, *theEstimator);
1458
1459 edm::ESHandle<TransientTrackingRecHitBuilder> builder = iSetup.getHandle(m_ttrhBuilderToken);
1460 this->TTRHBuilder = &(*builder);
1461 this->MuRHBuilder = new MuonTransientTrackingRecHitBuilder(theTrackingGeometry);
1462 if (debug_) {
1463 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1464 << " theTrackingGeometry.isValid() " << theTrackingGeometry.isValid();
1465 edm::LogVerbatim("GlobalTrackerMuonAlignment") << "get also the MuonTransientTrackingRecHitBuilder"
1466 << "\n";
1467 edm::LogVerbatim("GlobalTrackerMuonAlignment") << "get also the TransientTrackingRecHitBuilder"
1468 << "\n";
1469 }
1470 defineFitter = false;
1471 }
1472
1473
1474
1475 for (MuonCollection::const_iterator itMuon = smuons->begin(); itMuon != smuons->end(); ++itMuon) {
1476 if (debug_) {
1477 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1478 << " mu gM is GM Mu SaM tM " << itMuon->isGlobalMuon() << " " << itMuon->isMuon() << " "
1479 << itMuon->isStandAloneMuon() << " " << itMuon->isTrackerMuon() << " ";
1480 }
1481
1482
1483 TransientTrack muTT(itMuon->outerTrack(), magneticField_, trackingGeometry_);
1484 TrajectoryStateOnSurface innerMuTSOS = muTT.innermostMeasurementState();
1485 TrajectoryStateOnSurface outerMuTSOS = muTT.outermostMeasurementState();
1486
1487
1488 TransientTrack trackTT(itMuon->track(), magneticField_, trackingGeometry_);
1489 TrajectoryStateOnSurface outerTrackTSOS = trackTT.outermostMeasurementState();
1490 TrajectoryStateOnSurface innerTrackTSOS = trackTT.innermostMeasurementState();
1491
1492 GlobalPoint pointTrackIn = innerTrackTSOS.globalPosition();
1493 GlobalPoint pointTrackOut = outerTrackTSOS.globalPosition();
1494 float lenghtTrack = (pointTrackOut - pointTrackIn).mag();
1495 GlobalPoint pointMuonIn = innerMuTSOS.globalPosition();
1496 GlobalPoint pointMuonOut = outerMuTSOS.globalPosition();
1497 float lenghtMuon = (pointMuonOut - pointMuonIn).mag();
1498 GlobalVector momentumTrackOut = outerTrackTSOS.globalMomentum();
1499 GlobalVector momentumTrackIn = innerTrackTSOS.globalMomentum();
1500 float distanceInIn = (pointTrackIn - pointMuonIn).mag();
1501 float distanceInOut = (pointTrackIn - pointMuonOut).mag();
1502 float distanceOutIn = (pointTrackOut - pointMuonIn).mag();
1503 float distanceOutOut = (pointTrackOut - pointMuonOut).mag();
1504
1505 if (debug_) {
1506 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " pointTrackIn " << pointTrackIn;
1507 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Out " << pointTrackOut << " lenght " << lenghtTrack;
1508 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " point MuonIn " << pointMuonIn;
1509 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Out " << pointMuonOut << " lenght " << lenghtMuon;
1510 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1511 << " momeTrackIn Out " << momentumTrackIn << " " << momentumTrackOut;
1512 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1513 << " doi io ii oo " << distanceOutIn << " " << distanceInOut << " " << distanceInIn << " " << distanceOutOut;
1514 }
1515
1516 if (lenghtTrack < 90.)
1517 continue;
1518 if (lenghtMuon < 300.)
1519 continue;
1520 if (innerMuTSOS.globalMomentum().mag() < 5. || outerMuTSOS.globalMomentum().mag() < 5.)
1521 continue;
1522 if (momentumTrackIn.mag() < 15. || momentumTrackOut.mag() < 15.)
1523 continue;
1524 if (trackTT.charge() != muTT.charge())
1525 continue;
1526
1527 if (debug_)
1528 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " passed lenght momentum cuts";
1529
1530 GlobalVector GRm, GPm, Nl, Rm, Pm, Rt, Pt, Rt0;
1531 AlgebraicSymMatrix66 Cm, C0, Ce, C1;
1532
1533
1534 TrajectoryStateOnSurface extrapolationT;
1535 int tsosMuonIf = 0;
1536
1537 TrajectoryStateOnSurface muonFittedTSOS;
1538 TrajectoryStateOnSurface trackFittedTSOS;
1539
1540 if (isolatedMuonMode_) {
1541 if (debug_)
1542 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ------ Isolated (out-in) ----- ";
1543 const Surface& refSurface = innerMuTSOS.surface();
1544 ConstReferenceCountingPointer<TangentPlane> tpMuLocal(refSurface.tangentPlane(innerMuTSOS.localPosition()));
1545 Nl = tpMuLocal->normalVector();
1546 ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
1547 GlobalVector Ng = tpMuGlobal->normalVector();
1548 Surface* surf = (Surface*)&refSurface;
1549 const Plane* refPlane = dynamic_cast<Plane*>(surf);
1550 GlobalVector Nlp = refPlane->normalVector();
1551
1552 if (debug_) {
1553 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1554 << " alfa " << int(asin(Nl.x()) * 57.29578);
1555 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " global " << Ng.x() << " " << Ng.y() << " " << Ng.z();
1556 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " lp " << Nlp.x() << " " << Nlp.y() << " " << Nlp.z();
1557 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug")
1558 << " Nlocal Nglobal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " " << Ng.x() << " " << Ng.y() << " "
1559 << Ng.z() << " alfa " << static_cast<int>(asin(Nl.x()) * 57.29578);
1560 }
1561
1562
1563
1564
1565 if (!refitTrack_)
1566 extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
1567 else {
1568 GlobalTrackerMuonAlignment::trackFitter(itMuon->track(), trackTT, alongMomentum, trackFittedTSOS);
1569 if (trackFittedTSOS.isValid())
1570 extrapolationT = alongSmPr.propagate(trackFittedTSOS, refSurface);
1571 }
1572
1573 if (!extrapolationT.isValid()) {
1574 if (false & alarm)
1575 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
1576
1577 ;
1578 continue;
1579 }
1580 tsosMuonIf = 1;
1581 Rt = GlobalVector((extrapolationT.globalPosition()).x(),
1582 (extrapolationT.globalPosition()).y(),
1583 (extrapolationT.globalPosition()).z());
1584
1585 Pt = extrapolationT.globalMomentum();
1586
1587 GRm = GlobalVector(
1588 (innerMuTSOS.globalPosition()).x(), (innerMuTSOS.globalPosition()).y(), (innerMuTSOS.globalPosition()).z());
1589 GPm = innerMuTSOS.globalMomentum();
1590
1591 Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
1592 (outerTrackTSOS.globalPosition()).y(),
1593 (outerTrackTSOS.globalPosition()).z());
1594 Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + innerMuTSOS.cartesianError().matrix());
1595 C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
1596 Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
1597 C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
1598
1599 if (refitMuon_)
1600 GlobalTrackerMuonAlignment::muonFitter(itMuon->outerTrack(), muTT, oppositeToMomentum, muonFittedTSOS);
1601
1602 }
1603
1604 if (cosmicMuonMode_) {
1605
1606 if ((distanceOutIn <= distanceInOut) & (distanceOutIn <= distanceInIn) &
1607 (distanceOutIn <= distanceOutOut)) {
1608 if (debug_)
1609 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ----- Out - In -----";
1610
1611 const Surface& refSurface = innerMuTSOS.surface();
1612 ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
1613 Nl = tpMuGlobal->normalVector();
1614
1615
1616
1617 if (!refitTrack_)
1618 extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
1619 else {
1620 GlobalTrackerMuonAlignment::trackFitter(itMuon->track(), trackTT, alongMomentum, trackFittedTSOS);
1621 if (trackFittedTSOS.isValid())
1622 extrapolationT = alongSmPr.propagate(trackFittedTSOS, refSurface);
1623 }
1624
1625 if (!extrapolationT.isValid()) {
1626 if (false & alarm)
1627 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
1628
1629 ;
1630 continue;
1631 }
1632 if (debug_)
1633 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " extrapolationT.isValid " << extrapolationT.isValid();
1634
1635 tsosMuonIf = 1;
1636 Rt = GlobalVector((extrapolationT.globalPosition()).x(),
1637 (extrapolationT.globalPosition()).y(),
1638 (extrapolationT.globalPosition()).z());
1639
1640 Pt = extrapolationT.globalMomentum();
1641
1642 GRm = GlobalVector(
1643 (innerMuTSOS.globalPosition()).x(), (innerMuTSOS.globalPosition()).y(), (innerMuTSOS.globalPosition()).z());
1644 GPm = innerMuTSOS.globalMomentum();
1645 Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
1646 (outerTrackTSOS.globalPosition()).y(),
1647 (outerTrackTSOS.globalPosition()).z());
1648 Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + innerMuTSOS.cartesianError().matrix());
1649 C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
1650 Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
1651 C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
1652
1653 if (refitMuon_)
1654 GlobalTrackerMuonAlignment::muonFitter(itMuon->outerTrack(), muTT, oppositeToMomentum, muonFittedTSOS);
1655
1656 if (false & debug_) {
1657 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ->propDir " << propagator->propagationDirection();
1658 PropagationDirectionChooser Chooser;
1659 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1660 << " propDirCh " << Chooser.operator()(*outerTrackTSOS.freeState(), refSurface) << " Ch == along "
1661 << (alongMomentum == Chooser.operator()(*outerTrackTSOS.freeState(), refSurface));
1662 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1663 << " --- Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1664 << " alfa " << int(asin(Nl.x()) * 57.29578);
1665 }
1666 }
1667
1668 else if ((distanceInOut <= distanceInIn) & (distanceInOut <= distanceOutIn) &
1669 (distanceInOut <= distanceOutOut)) {
1670 if (debug_)
1671 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ----- In - Out -----";
1672
1673 const Surface& refSurface = outerMuTSOS.surface();
1674 ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
1675 Nl = tpMuGlobal->normalVector();
1676
1677
1678
1679 if (!refitTrack_)
1680 extrapolationT = oppositeSmPr.propagate(innerTrackTSOS, refSurface);
1681 else {
1682 GlobalTrackerMuonAlignment::trackFitter(itMuon->track(), trackTT, oppositeToMomentum, trackFittedTSOS);
1683 if (trackFittedTSOS.isValid())
1684 extrapolationT = oppositeSmPr.propagate(trackFittedTSOS, refSurface);
1685 }
1686
1687 if (!extrapolationT.isValid()) {
1688 if (false & alarm)
1689 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
1690 << "\a\a\a\a\a\a\a\a" << extrapolationT.isValid();
1691 continue;
1692 }
1693 if (debug_)
1694 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " extrapolationT.isValid " << extrapolationT.isValid();
1695
1696 tsosMuonIf = 2;
1697 Rt = GlobalVector((extrapolationT.globalPosition()).x(),
1698 (extrapolationT.globalPosition()).y(),
1699 (extrapolationT.globalPosition()).z());
1700
1701 Pt = extrapolationT.globalMomentum();
1702
1703 GRm = GlobalVector(
1704 (outerMuTSOS.globalPosition()).x(), (outerMuTSOS.globalPosition()).y(), (outerMuTSOS.globalPosition()).z());
1705 GPm = outerMuTSOS.globalMomentum();
1706 Rt0 = GlobalVector((innerTrackTSOS.globalPosition()).x(),
1707 (innerTrackTSOS.globalPosition()).y(),
1708 (innerTrackTSOS.globalPosition()).z());
1709 Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + outerMuTSOS.cartesianError().matrix());
1710 C0 = AlgebraicSymMatrix66(innerTrackTSOS.cartesianError().matrix());
1711 Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
1712 C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
1713
1714 if (refitMuon_)
1715 GlobalTrackerMuonAlignment::muonFitter(itMuon->outerTrack(), muTT, alongMomentum, muonFittedTSOS);
1716
1717 if (false & debug_) {
1718 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ->propDir " << propagator->propagationDirection();
1719 PropagationDirectionChooser Chooser;
1720 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1721 << " propDirCh " << Chooser.operator()(*innerTrackTSOS.freeState(), refSurface) << " Ch == oppisite "
1722 << (oppositeToMomentum == Chooser.operator()(*innerTrackTSOS.freeState(), refSurface));
1723 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1724 << " --- Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1725 << " alfa " << int(asin(Nl.x()) * 57.29578);
1726 }
1727 }
1728
1729 else if ((distanceOutOut <= distanceInOut) & (distanceOutOut <= distanceInIn) &
1730 (distanceOutOut <= distanceOutIn)) {
1731 if (debug_)
1732 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ----- Out - Out -----";
1733
1734
1735 continue;
1736
1737 const Surface& refSurface = outerMuTSOS.surface();
1738 ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
1739 Nl = tpMuGlobal->normalVector();
1740
1741
1742 extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
1743
1744 if (!extrapolationT.isValid()) {
1745 if (alarm)
1746 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " !!!!!!!!! Catastrophe Out-Out extrapolationT.isValid "
1747 << "\a\a\a\a\a\a\a\a" << extrapolationT.isValid();
1748 continue;
1749 }
1750 if (debug_)
1751 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " extrapolationT.isValid " << extrapolationT.isValid();
1752
1753 tsosMuonIf = 2;
1754 Rt = GlobalVector((extrapolationT.globalPosition()).x(),
1755 (extrapolationT.globalPosition()).y(),
1756 (extrapolationT.globalPosition()).z());
1757
1758 Pt = extrapolationT.globalMomentum();
1759
1760 GRm = GlobalVector(
1761 (outerMuTSOS.globalPosition()).x(), (outerMuTSOS.globalPosition()).y(), (outerMuTSOS.globalPosition()).z());
1762 GPm = outerMuTSOS.globalMomentum();
1763 Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
1764 (outerTrackTSOS.globalPosition()).y(),
1765 (outerTrackTSOS.globalPosition()).z());
1766 Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + outerMuTSOS.cartesianError().matrix());
1767 C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
1768 Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
1769 C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
1770
1771 if (debug_) {
1772 PropagationDirectionChooser Chooser;
1773 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1774 << " propDirCh " << Chooser.operator()(*outerTrackTSOS.freeState(), refSurface) << " Ch == along "
1775 << (alongMomentum == Chooser.operator()(*outerTrackTSOS.freeState(), refSurface));
1776 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1777 << " --- Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1778 << " alfa " << int(asin(Nl.x()) * 57.29578);
1779 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1780 << " Nornal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1781 << " alfa " << int(asin(Nl.x()) * 57.29578);
1782 }
1783 }
1784 else {
1785 if (alarm)
1786 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1787 << " ------- !!!!!!!!!!!!!!! No proper Out - In \a\a\a\a\a\a\a";
1788 continue;
1789 }
1790
1791 }
1792
1793 if (tsosMuonIf == 0) {
1794 if (info) {
1795 edm::LogVerbatim("GlobalTrackerMuonAlignment") << "No tsosMuon !!!!!!";
1796 continue;
1797 }
1798 }
1799 TrajectoryStateOnSurface tsosMuon;
1800 if (tsosMuonIf == 1)
1801 tsosMuon = muTT.innermostMeasurementState();
1802 else
1803 tsosMuon = muTT.outermostMeasurementState();
1804
1805
1806 AlgebraicVector4 LPRm;
1807 GlobalTrackerMuonAlignment::misalignMuonL(GRm, GPm, Nl, Rt, Rm, Pm, LPRm, extrapolationT, tsosMuon);
1808
1809 if (refitTrack_) {
1810 if (!trackFittedTSOS.isValid()) {
1811 if (info)
1812 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ================= trackFittedTSOS notValid !!!!!!!! ";
1813 continue;
1814 }
1815 if (debug_)
1816 this->debugTrajectorySOS(" trackFittedTSOS ", trackFittedTSOS);
1817 }
1818
1819 if (refitMuon_) {
1820 if (!muonFittedTSOS.isValid()) {
1821 if (info)
1822 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ================= muonFittedTSOS notValid !!!!!!!! ";
1823 continue;
1824 }
1825 if (debug_)
1826 this->debugTrajectorySOS(" muonFittedTSOS ", muonFittedTSOS);
1827 Rm = GlobalVector((muonFittedTSOS.globalPosition()).x(),
1828 (muonFittedTSOS.globalPosition()).y(),
1829 (muonFittedTSOS.globalPosition()).z());
1830 Pm = muonFittedTSOS.globalMomentum();
1831 LPRm = AlgebraicVector4(muonFittedTSOS.localParameters().vector()(1),
1832 muonFittedTSOS.localParameters().vector()(2),
1833 muonFittedTSOS.localParameters().vector()(3),
1834 muonFittedTSOS.localParameters().vector()(4));
1835 }
1836 GlobalVector resR = Rm - Rt;
1837 GlobalVector resP0 = Pm - Pt;
1838 GlobalVector resP = Pm / Pm.mag() - Pt / Pt.mag();
1839 float RelMomResidual = (Pm.mag() - Pt.mag()) / (Pt.mag() + 1.e-6);
1840 ;
1841
1842 AlgebraicVector6 Vm;
1843 Vm(0) = resR.x();
1844 Vm(1) = resR.y();
1845 Vm(2) = resR.z();
1846 Vm(3) = resP0.x();
1847 Vm(4) = resP0.y();
1848 Vm(5) = resP0.z();
1849 float Rmuon = Rm.perp();
1850 float Zmuon = Rm.z();
1851 float alfa_x = atan2(Nl.x(), Nl.y()) * 57.29578;
1852
1853 if (debug_) {
1854 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1855 << " Nx Ny Nz alfa_x " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " " << alfa_x;
1856 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug")
1857 << " Rm " << Rm << "\n Rt " << Rt << "\n resR " << resR << "\n resP " << resP << " dp/p " << RelMomResidual;
1858 }
1859
1860 double chi_d = 0;
1861 for (int i = 0; i <= 5; i++)
1862 chi_d += Vm(i) * Vm(i) / Cm(i, i);
1863
1864 AlgebraicVector5 Vml(tsosMuon.localParameters().vector() - extrapolationT.localParameters().vector());
1865 AlgebraicSymMatrix55 m(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
1866 AlgebraicSymMatrix55 Cml(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
1867 bool ierrLoc = !m.Invert();
1868 if (ierrLoc && debug_ && info) {
1869 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ==== Error inverting Local covariance matrix ==== ";
1870 continue;
1871 }
1872 double chi_Loc = ROOT::Math::Similarity(Vml, m);
1873 if (debug_)
1874 edm::LogVerbatim("GlobalTrackerMuonAlignment")
1875 << " chi_Loc px/pz/err " << chi_Loc << " " << Vml(1) / std::sqrt(Cml(1, 1));
1876
1877 if (Pt.mag() < 15.)
1878 continue;
1879 if (Pm.mag() < 5.)
1880 continue;
1881
1882
1883
1884
1885
1886
1887
1888
1889
1890
1891
1892
1893
1894
1895 if (fabs(resR.x()) > 20.)
1896 continue;
1897 if (fabs(resR.y()) > 20.)
1898 continue;
1899 if (fabs(resR.z()) > 20.)
1900 continue;
1901 if (fabs(resR.mag()) > 30.)
1902 continue;
1903 if (fabs(resP.x()) > 0.06)
1904 continue;
1905 if (fabs(resP.y()) > 0.06)
1906 continue;
1907 if (fabs(resP.z()) > 0.06)
1908 continue;
1909 if (chi_d > 40.)
1910 continue;
1911
1912
1913
1914
1915
1916
1917
1918
1919
1920
1921
1922
1923
1924
1925
1926
1927
1928
1929
1930 if (Rmuon < 120. || Rmuon > 450.)
1931 continue;
1932 if (Zmuon < -720. || Zmuon > 720.)
1933 continue;
1934 MuSelect = " Barrel+EndCaps";
1935
1936 if (debug_)
1937 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " .............. passed all cuts";
1938
1939 N_track++;
1940
1941
1942 GlobalTrackerMuonAlignment::gradientGlobalAlg(Rt, Pt, Rm, Nl, Cm);
1943 GlobalTrackerMuonAlignment::gradientGlobal(Rt, Pt, Rm, Pm, Nl, Cm);
1944
1945 CLHEP::HepSymMatrix covLoc(4, 0);
1946 for (int i = 1; i <= 4; i++)
1947 for (int j = 1; j <= i; j++) {
1948 covLoc(i, j) = (tsosMuon.localError().matrix() + extrapolationT.localError().matrix())(i, j);
1949
1950 }
1951
1952 const Surface& refSurface = tsosMuon.surface();
1953 CLHEP::HepMatrix rotLoc(3, 3, 0);
1954 rotLoc(1, 1) = refSurface.rotation().xx();
1955 rotLoc(1, 2) = refSurface.rotation().xy();
1956 rotLoc(1, 3) = refSurface.rotation().xz();
1957
1958 rotLoc(2, 1) = refSurface.rotation().yx();
1959 rotLoc(2, 2) = refSurface.rotation().yy();
1960 rotLoc(2, 3) = refSurface.rotation().yz();
1961
1962 rotLoc(3, 1) = refSurface.rotation().zx();
1963 rotLoc(3, 2) = refSurface.rotation().zy();
1964 rotLoc(3, 3) = refSurface.rotation().zz();
1965
1966 CLHEP::HepVector posLoc(3);
1967 posLoc(1) = refSurface.position().x();
1968 posLoc(2) = refSurface.position().y();
1969 posLoc(3) = refSurface.position().z();
1970
1971 GlobalTrackerMuonAlignment::gradientLocal(Rt, Pt, Rm, Pm, Nl, covLoc, rotLoc, posLoc, LPRm);
1972
1973 if (debug_) {
1974 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Norm " << Nl;
1975 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " posLoc " << posLoc.T();
1976 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " rotLoc " << rotLoc;
1977 }
1978
1979
1980 histo->Fill(itMuon->track()->pt());
1981
1982
1983 histo2->Fill(Pt.mag());
1984 histo3->Fill((PI / 2. - itMuon->track()->outerTheta()));
1985 histo4->Fill(itMuon->track()->phi());
1986 histo5->Fill(Rmuon);
1987 histo6->Fill(Zmuon);
1988 histo7->Fill(RelMomResidual);
1989
1990 histo8->Fill(chi_d);
1991
1992 histo101->Fill(Zmuon, Rmuon);
1993 histo101->Fill(Rt0.z(), Rt0.perp());
1994 histo102->Fill(Rt0.x(), Rt0.y());
1995 histo102->Fill(Rm.x(), Rm.y());
1996
1997 histo11->Fill(resR.mag());
1998 if (fabs(Nl.x()) < 0.98)
1999 histo12->Fill(resR.x());
2000 if (fabs(Nl.y()) < 0.98)
2001 histo13->Fill(resR.y());
2002 if (fabs(Nl.z()) < 0.98)
2003 histo14->Fill(resR.z());
2004 histo15->Fill(resP.x());
2005 histo16->Fill(resP.y());
2006 histo17->Fill(resP.z());
2007
2008 if ((fabs(Nl.x()) < 0.98) && (fabs(Nl.y()) < 0.98) && (fabs(Nl.z()) < 0.98)) {
2009 histo18->Fill(std::sqrt(C0(0, 0)));
2010 histo19->Fill(std::sqrt(C1(0, 0)));
2011 histo20->Fill(std::sqrt(C1(0, 0) + Ce(0, 0)));
2012 }
2013 if (fabs(Nl.x()) < 0.98)
2014 histo21->Fill(Vm(0) / std::sqrt(Cm(0, 0)));
2015 if (fabs(Nl.y()) < 0.98)
2016 histo22->Fill(Vm(1) / std::sqrt(Cm(1, 1)));
2017 if (fabs(Nl.z()) < 0.98)
2018 histo23->Fill(Vm(2) / std::sqrt(Cm(2, 2)));
2019 histo24->Fill(Vm(3) / std::sqrt(C1(3, 3) + Ce(3, 3)));
2020 histo25->Fill(Vm(4) / std::sqrt(C1(4, 4) + Ce(4, 4)));
2021 histo26->Fill(Vm(5) / std::sqrt(C1(5, 5) + Ce(5, 5)));
2022 histo27->Fill(Nl.x());
2023 histo28->Fill(Nl.y());
2024 histo29->Fill(lenghtTrack);
2025 histo30->Fill(lenghtMuon);
2026 histo31->Fill(chi_Loc);
2027 histo32->Fill(Vml(1) / std::sqrt(Cml(1, 1)));
2028 histo33->Fill(Vml(2) / std::sqrt(Cml(2, 2)));
2029 histo34->Fill(Vml(3) / std::sqrt(Cml(3, 3)));
2030 histo35->Fill(Vml(4) / std::sqrt(Cml(4, 4)));
2031
2032 if (debug_) {
2033
2034 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " diag 0[ " << C0(0, 0) << " " << C0(1, 1) << " " << C0(2, 2)
2035 << " " << C0(3, 3) << " " << C0(4, 4) << " " << C0(5, 5) << " ]";
2036 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " diag e[ " << Ce(0, 0) << " " << Ce(1, 1) << " " << Ce(2, 2)
2037 << " " << Ce(3, 3) << " " << Ce(4, 4) << " " << Ce(5, 5) << " ]";
2038 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " diag 1[ " << C1(0, 0) << " " << C1(1, 1) << " " << C1(2, 2)
2039 << " " << C1(3, 3) << " " << C1(4, 4) << " " << C1(5, 5) << " ]";
2040 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Rm " << Rm.x() << " " << Rm.y() << " " << Rm.z() << " Pm "
2041 << Pm.x() << " " << Pm.y() << " " << Pm.z();
2042 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Rt " << Rt.x() << " " << Rt.y() << " " << Rt.z() << " Pt "
2043 << Pt.x() << " " << Pt.y() << " " << Pt.z();
2044 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Nl*(Rm-Rt) " << Nl.dot(Rm - Rt);
2045 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " resR " << resR.x() << " " << resR.y() << " " << resR.z()
2046 << " resP " << resP.x() << " " << resP.y() << " " << resP.z();
2047 edm::LogVerbatim("GlobalTrackerMuonAlignment")
2048 << " Rm-t " << (Rm - Rt).x() << " " << (Rm - Rt).y() << " " << (Rm - Rt).z() << " Pm-t " << (Pm - Pt).x()
2049 << " " << (Pm - Pt).y() << " " << (Pm - Pt).z();
2050 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Vm " << Vm;
2051 edm::LogVerbatim("GlobalTrackerMuonAlignment")
2052 << " +- " << std::sqrt(Cm(0, 0)) << " " << std::sqrt(Cm(1, 1)) << " " << std::sqrt(Cm(2, 2)) << " "
2053 << std::sqrt(Cm(3, 3)) << " " << std::sqrt(Cm(4, 4)) << " " << std::sqrt(Cm(5, 5));
2054 edm::LogVerbatim("GlobalTrackerMuonAlignment")
2055 << " Pmuon Ptrack dP/Ptrack " << itMuon->outerTrack()->p() << " " << itMuon->track()->outerP() << " "
2056 << (itMuon->outerTrack()->p() - itMuon->track()->outerP()) / itMuon->track()->outerP();
2057 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " cov matrix ";
2058 edm::LogVerbatim("GlobalTrackerMuonAlignment") << Cm;
2059 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " diag [ " << Cm(0, 0) << " " << Cm(1, 1) << " " << Cm(2, 2)
2060 << " " << Cm(3, 3) << " " << Cm(4, 4) << " " << Cm(5, 5) << " ]";
2061
2062 AlgebraicSymMatrix66 Ro;
2063 double Diag[6];
2064 for (int i = 0; i <= 5; i++)
2065 Diag[i] = std::sqrt(Cm(i, i));
2066 for (int i = 0; i <= 5; i++)
2067 for (int j = 0; j <= 5; j++)
2068 Ro(i, j) = Cm(i, j) / Diag[i] / Diag[j];
2069 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " correlation matrix ";
2070 edm::LogVerbatim("GlobalTrackerMuonAlignment") << Ro;
2071
2072 AlgebraicSymMatrix66 CmI;
2073 for (int i = 0; i <= 5; i++)
2074 for (int j = 0; j <= 5; j++)
2075 CmI(i, j) = Cm(i, j);
2076
2077 bool ierr = !CmI.Invert();
2078 if (ierr) {
2079 if (alarm || debug_)
2080 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Error inverse covariance matrix !!!!!!!!!!!";
2081 continue;
2082 }
2083 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " inverse cov matrix ";
2084 edm::LogVerbatim("GlobalTrackerMuonAlignment") << Cm;
2085
2086 double chi = ROOT::Math::Similarity(Vm, CmI);
2087 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " chi chi_d " << chi << " " << chi_d;
2088 }
2089
2090 }
2091
2092 }
2093
2094
2095 void GlobalTrackerMuonAlignment::gradientGlobalAlg(
2096 GlobalVector& Rt, GlobalVector& Pt, GlobalVector& Rm, GlobalVector& Nl, AlgebraicSymMatrix66& Cm) {
2097
2098
2099
2100 AlgebraicMatrix33 Jac;
2101 AlgebraicVector3 Wi, R_m, R_t, P_t, Norm, dR;
2102
2103 R_m(0) = Rm.x();
2104 R_m(1) = Rm.y();
2105 R_m(2) = Rm.z();
2106 R_t(0) = Rt.x();
2107 R_t(1) = Rt.y();
2108 R_t(2) = Rt.z();
2109 P_t(0) = Pt.x();
2110 P_t(1) = Pt.y();
2111 P_t(2) = Pt.z();
2112 Norm(0) = Nl.x();
2113 Norm(1) = Nl.y();
2114 Norm(2) = Nl.z();
2115
2116 for (int i = 0; i <= 2; i++) {
2117 if (Cm(i, i) > 1.e-20)
2118 Wi(i) = 1. / Cm(i, i);
2119 else
2120 Wi(i) = 1.e-10;
2121 dR(i) = R_m(i) - R_t(i);
2122 }
2123
2124 float PtN = P_t(0) * Norm(0) + P_t(1) * Norm(1) + P_t(2) * Norm(2);
2125
2126 Jac(0, 0) = 1. - P_t(0) * Norm(0) / PtN;
2127 Jac(0, 1) = -P_t(0) * Norm(1) / PtN;
2128 Jac(0, 2) = -P_t(0) * Norm(2) / PtN;
2129
2130 Jac(1, 0) = -P_t(1) * Norm(0) / PtN;
2131 Jac(1, 1) = 1. - P_t(1) * Norm(1) / PtN;
2132 Jac(1, 2) = -P_t(1) * Norm(2) / PtN;
2133
2134 Jac(2, 0) = -P_t(2) * Norm(0) / PtN;
2135 Jac(2, 1) = -P_t(2) * Norm(1) / PtN;
2136 Jac(2, 2) = 1. - P_t(2) * Norm(2) / PtN;
2137
2138 AlgebraicSymMatrix33 Itr;
2139
2140 for (int i = 0; i <= 2; i++)
2141 for (int j = 0; j <= 2; j++) {
2142 if (j < i)
2143 continue;
2144 Itr(i, j) = 0.;
2145 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ij " << i << " " << j;
2146 for (int k = 0; k <= 2; k++) {
2147 Itr(i, j) += Jac(k, i) * Wi(k) * Jac(k, j);
2148 }
2149 }
2150
2151 for (int i = 0; i <= 2; i++)
2152 for (int j = 0; j <= 2; j++) {
2153 if (j < i)
2154 continue;
2155 Inf(i, j) += Itr(i, j);
2156 }
2157
2158 AlgebraicVector3 Gtr(0., 0., 0.);
2159 for (int i = 0; i <= 2; i++)
2160 for (int k = 0; k <= 2; k++)
2161 Gtr(i) += dR(k) * Wi(k) * Jac(k, i);
2162 for (int i = 0; i <= 2; i++)
2163 Gfr(i) += Gtr(i);
2164
2165 if (debug_) {
2166 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Wi " << Wi;
2167 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " N " << Norm;
2168 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " P_t " << P_t;
2169 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " (Pt*N) " << PtN;
2170 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " dR " << dR;
2171 edm::LogVerbatim("GlobalTrackerMuonAlignment")
2172 << " +/- " << 1. / std::sqrt(Wi(0)) << " " << 1. / std::sqrt(Wi(1)) << " " << 1. / std::sqrt(Wi(2)) << " ";
2173 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Jacobian dr/ddx ";
2174 edm::LogVerbatim("GlobalTrackerMuonAlignment") << Jac;
2175 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " G-- " << Gtr;
2176 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Itrack ";
2177 edm::LogVerbatim("GlobalTrackerMuonAlignment") << Itr;
2178 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Gfr " << Gfr;
2179 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -- Inf --";
2180 edm::LogVerbatim("GlobalTrackerMuonAlignment") << Inf;
2181 }
2182
2183 return;
2184 }
2185
2186
2187 void GlobalTrackerMuonAlignment::gradientGlobal(GlobalVector& GRt,
2188 GlobalVector& GPt,
2189 GlobalVector& GRm,
2190 GlobalVector& GPm,
2191 GlobalVector& GNorm,
2192 AlgebraicSymMatrix66& GCov) {
2193
2194
2195
2196
2197
2198 bool alarm = false;
2199 bool info = false;
2200
2201 int Nd = 6;
2202
2203
2204 CLHEP::HepSymMatrix w(Nd, 0);
2205 for (int i = 1; i <= Nd; i++)
2206 for (int j = 1; j <= Nd; j++) {
2207 if (j <= i)
2208 w(i, j) = GCov(i - 1, j - 1);
2209
2210
2211 if ((i == j) && (i <= 3) && (GCov(i - 1, j - 1) < 1.e-20))
2212 w(i, j) = 1.e20;
2213 if (i != j)
2214 w(i, j) = 0.;
2215 }
2216
2217
2218
2219
2220 CLHEP::HepVector V(Nd), Rt(3), Pt(3), Rm(3), Pm(3), Norm(3);
2221 Rt(1) = GRt.x();
2222 Rt(2) = GRt.y();
2223 Rt(3) = GRt.z();
2224 Pt(1) = GPt.x();
2225 Pt(2) = GPt.y();
2226 Pt(3) = GPt.z();
2227 Rm(1) = GRm.x();
2228 Rm(2) = GRm.y();
2229 Rm(3) = GRm.z();
2230 Pm(1) = GPm.x();
2231 Pm(2) = GPm.y();
2232 Pm(3) = GPm.z();
2233 Norm(1) = GNorm.x();
2234 Norm(2) = GNorm.y();
2235 Norm(3) = GNorm.z();
2236
2237 V = dsum(Rm - Rt, Pm - Pt);
2238 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " V " << V.T();
2239
2240
2241 double PmN = CLHEP_dot(Pm, Norm);
2242
2243 CLHEP::HepMatrix Jac(Nd, Nd, 0);
2244 for (int i = 1; i <= 3; i++)
2245 for (int j = 1; j <= 3; j++) {
2246 Jac(i, j) = Pm(i) * Norm(j) / PmN;
2247 if (i == j)
2248 Jac(i, j) -= 1.;
2249 }
2250
2251
2252 Jac(4, 4) = 0.;
2253 Jac(5, 4) = -Pm(3);
2254 Jac(6, 4) = Pm(2);
2255 Jac(4, 5) = Pm(3);
2256 Jac(5, 5) = 0.;
2257 Jac(6, 5) = -Pm(1);
2258 Jac(4, 6) = -Pm(2);
2259 Jac(5, 6) = Pm(1);
2260 Jac(6, 6) = 0.;
2261
2262 CLHEP::HepVector dsda(3);
2263 dsda(1) = (Norm(2) * Rm(3) - Norm(3) * Rm(2)) / PmN;
2264 dsda(2) = (Norm(3) * Rm(1) - Norm(1) * Rm(3)) / PmN;
2265 dsda(3) = (Norm(1) * Rm(2) - Norm(2) * Rm(1)) / PmN;
2266
2267
2268 Jac(1, 4) = Pm(1) * dsda(1);
2269 Jac(2, 4) = -Rm(3) + Pm(2) * dsda(1);
2270 Jac(3, 4) = Rm(2) + Pm(3) * dsda(1);
2271
2272 Jac(1, 5) = Rm(3) + Pm(1) * dsda(2);
2273 Jac(2, 5) = Pm(2) * dsda(2);
2274 Jac(3, 5) = -Rm(1) + Pm(3) * dsda(2);
2275
2276 Jac(1, 6) = -Rm(2) + Pm(1) * dsda(3);
2277 Jac(2, 6) = Rm(1) + Pm(2) * dsda(3);
2278 Jac(3, 6) = Pm(3) * dsda(3);
2279
2280 CLHEP::HepSymMatrix W(Nd, 0);
2281 int ierr;
2282 W = w.inverse(ierr);
2283 if (ierr != 0) {
2284 if (alarm)
2285 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " gradientGlobal: inversion of matrix w fail ";
2286 return;
2287 }
2288
2289 CLHEP::HepMatrix W_Jac(Nd, Nd, 0);
2290 W_Jac = Jac.T() * W;
2291
2292 CLHEP::HepVector grad3(Nd);
2293 grad3 = W_Jac * V;
2294
2295 CLHEP::HepMatrix hess3(Nd, Nd);
2296 hess3 = Jac.T() * W * Jac;
2297
2298
2299 Grad += grad3;
2300 Hess += hess3;
2301
2302 CLHEP::HepVector d3I = CLHEP::solve(Hess, -Grad);
2303 int iEr3I;
2304 CLHEP::HepMatrix Errd3I = Hess.inverse(iEr3I);
2305 if (iEr3I != 0) {
2306 if (alarm)
2307 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " gradientGlobal error inverse Hess matrix !!!!!!!!!!!";
2308 }
2309
2310 if (info || debug_) {
2311 edm::LogVerbatim("GlobalTrackerMuonAlignment")
2312 << " dG " << d3I(1) << " " << d3I(2) << " " << d3I(3) << " " << d3I(4) << " " << d3I(5) << " " << d3I(6);
2313 ;
2314 edm::LogVerbatim("GlobalTrackerMuonAlignment")
2315 << " +- " << std::sqrt(Errd3I(1, 1)) << " " << std::sqrt(Errd3I(2, 2)) << " " << std::sqrt(Errd3I(3, 3))
2316 << " " << std::sqrt(Errd3I(4, 4)) << " " << std::sqrt(Errd3I(5, 5)) << " " << std::sqrt(Errd3I(6, 6));
2317 }
2318
2319 #ifdef CHECK_OF_DERIVATIVES
2320
2321
2322 CLHEP::HepVector d(3, 0), a(3, 0);
2323 CLHEP::HepMatrix T(3, 3, 1);
2324
2325 CLHEP::HepMatrix Ti = T.T();
2326
2327
2328 double A = CLHEP_dot(Ti * Pm, Norm);
2329 double B = CLHEP_dot((Rt - Ti * Rm + Ti * d), Norm);
2330 double s0 = B / A;
2331
2332 CLHEP::HepVector r0(3, 0), p0(3, 0);
2333 r0 = Ti * Rm - Ti * d + s0 * (Ti * Pm) - Rt;
2334 p0 = Ti * Pm - Pt;
2335
2336 double delta = 0.0001;
2337
2338 int ii = 3;
2339 d(ii) += delta;
2340
2341
2342
2343 Ti = T.T();
2344
2345
2346 A = CLHEP_dot(Ti * Pm, Norm);
2347 B = CLHEP_dot((Rt - Ti * Rm + Ti * d), Norm);
2348 double s = B / A;
2349
2350 CLHEP::HepVector r(3, 0), p(3, 0);
2351 r = Ti * Rm - Ti * d + s * (Ti * Pm) - Rt;
2352 p = Ti * Pm - Pt;
2353
2354 edm::LogVerbatim("GlobalTrackerMuonAlignment")
2355 << " s0 s num dsda(" << ii << ") " << s0 << " " << s << " " << (s - s0) / delta << " " << dsda(ii);
2356
2357 edm::LogVerbatim("GlobalTrackerMuonAlignment")
2358 << " -- An d(r,p)/d(" << ii << ") " << Jac(1, ii) << " " << Jac(2, ii) << " " << Jac(3, ii) << " " << Jac(4, ii)
2359 << " " << Jac(5, ii) << " " << Jac(6, ii);
2360 edm::LogVerbatim("GlobalTrackerMuonAlignment")
2361 << " Nu d(r,p)/d(" << ii << ") " << (r(1) - r0(1)) / delta << " " << (r(2) - r0(2)) / delta << " "
2362 << (r(3) - r0(3)) / delta << " " << (p(1) - p0(1)) / delta << " " << (p(2) - p0(2)) / delta << " "
2363 << (p(3) - p0(3)) / delta;
2364
2365 edm::LogVerbatim("GlobalTrackerMuonAlignment")
2366 << " -- An d(r,p)/a(" << ii << ") " << Jac(1, ii + 3) << " " << Jac(2, ii + 3) << " " << Jac(3, ii + 3) << " "
2367 << Jac(4, ii + 3) << " " << Jac(5, ii + 3) << " " << Jac(6, ii + 3);
2368 edm::LogVerbatim("GlobalTrackerMuonAlignment")
2369 << " Nu d(r,p)/a(" << ii << ") " << (r(1) - r0(1)) / delta << " " << (r(2) - r0(2)) / delta << " "
2370 << (r(3) - r0(3)) / delta << " " << (p(1) - p0(1)) / delta << " " << (p(2) - p0(2)) / delta << " "
2371 << (p(3) - p0(3)) / delta;
2372
2373 #endif
2374
2375 return;
2376 }
2377
2378
2379 void GlobalTrackerMuonAlignment::gradientLocal(GlobalVector& GRt,
2380 GlobalVector& GPt,
2381 GlobalVector& GRm,
2382 GlobalVector& GPm,
2383 GlobalVector& GNorm,
2384 CLHEP::HepSymMatrix& covLoc,
2385 CLHEP::HepMatrix& rotLoc,
2386 CLHEP::HepVector& R0,
2387 AlgebraicVector4& LPRm) {
2388
2389
2390
2391
2392 bool alarm = true;
2393
2394 bool info = false;
2395
2396 if (debug_)
2397 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " gradientLocal ";
2398
2399
2400
2401
2402
2403
2404
2405
2406
2407
2408
2409
2410
2411
2412
2413
2414
2415
2416 CLHEP::HepVector Rt(3), Pt(3), Rm(3), Pm(3), Norm(3);
2417 Rt(1) = GRt.x();
2418 Rt(2) = GRt.y();
2419 Rt(3) = GRt.z();
2420 Pt(1) = GPt.x();
2421 Pt(2) = GPt.y();
2422 Pt(3) = GPt.z();
2423 Rm(1) = GRm.x();
2424 Rm(2) = GRm.y();
2425 Rm(3) = GRm.z();
2426 Pm(1) = GPm.x();
2427 Pm(2) = GPm.y();
2428 Pm(3) = GPm.z();
2429 Norm(1) = GNorm.x();
2430 Norm(2) = GNorm.y();
2431 Norm(3) = GNorm.z();
2432
2433 CLHEP::HepVector V(4), Rml(3), Pml(3), Rtl(3), Ptl(3);
2434
2435
2436
2437
2438
2439
2440
2441 Rml = rotLoc * (Rm - R0);
2442 Rtl = rotLoc * (Rt - R0);
2443 Pml = rotLoc * Pm;
2444 Ptl = rotLoc * Pt;
2445
2446 V(1) = LPRm(0) - Ptl(1) / Ptl(3);
2447 V(2) = LPRm(1) - Ptl(2) / Ptl(3);
2448 V(3) = LPRm(2) - Rtl(1);
2449 V(4) = LPRm(3) - Rtl(2);
2450
2451
2452
2453
2454
2455
2456
2457
2458
2459
2460
2461
2462
2463
2464 CLHEP::HepSymMatrix W = covLoc;
2465
2466 int ierr;
2467 W.invert(ierr);
2468 if (ierr != 0) {
2469 if (alarm)
2470 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " gradientLocal: inversion of matrix W fail ";
2471 return;
2472 }
2473
2474
2475
2476
2477
2478
2479
2480
2481
2482 CLHEP::HepMatrix JacToLoc(4, 6, 0);
2483 for (int i = 1; i <= 2; i++)
2484 for (int j = 1; j <= 3; j++) {
2485 JacToLoc(i, j + 3) = (rotLoc(i, j) - rotLoc(3, j) * Pml(i) / Pml(3)) / Pml(3);
2486 JacToLoc(i + 2, j) = rotLoc(i, j);
2487 }
2488
2489
2490
2491 double PmN = CLHEP_dot(Pm, Norm);
2492
2493 CLHEP::HepMatrix Jac(6, 6, 0);
2494 for (int i = 1; i <= 3; i++)
2495 for (int j = 1; j <= 3; j++) {
2496 Jac(i, j) = Pm(i) * Norm(j) / PmN;
2497 if (i == j)
2498 Jac(i, j) -= 1.;
2499 }
2500
2501
2502 Jac(4, 4) = 0.;
2503 Jac(5, 4) = -Pm(3);
2504 Jac(6, 4) = Pm(2);
2505 Jac(4, 5) = Pm(3);
2506 Jac(5, 5) = 0.;
2507 Jac(6, 5) = -Pm(1);
2508 Jac(4, 6) = -Pm(2);
2509 Jac(5, 6) = Pm(1);
2510 Jac(6, 6) = 0.;
2511
2512 CLHEP::HepVector dsda(3);
2513 dsda(1) = (Norm(2) * Rm(3) - Norm(3) * Rm(2)) / PmN;
2514 dsda(2) = (Norm(3) * Rm(1) - Norm(1) * Rm(3)) / PmN;
2515 dsda(3) = (Norm(1) * Rm(2) - Norm(2) * Rm(1)) / PmN;
2516
2517
2518 Jac(1, 4) = Pm(1) * dsda(1);
2519 Jac(2, 4) = -Rm(3) + Pm(2) * dsda(1);
2520 Jac(3, 4) = Rm(2) + Pm(3) * dsda(1);
2521
2522 Jac(1, 5) = Rm(3) + Pm(1) * dsda(2);
2523 Jac(2, 5) = Pm(2) * dsda(2);
2524 Jac(3, 5) = -Rm(1) + Pm(3) * dsda(2);
2525
2526 Jac(1, 6) = -Rm(2) + Pm(1) * dsda(3);
2527 Jac(2, 6) = Rm(1) + Pm(2) * dsda(3);
2528 Jac(3, 6) = Pm(3) * dsda(3);
2529
2530
2531 CLHEP::HepMatrix JacCorLoc(4, 6, 0);
2532 JacCorLoc = JacToLoc * Jac;
2533
2534
2535 CLHEP::HepMatrix W_Jac(6, 4, 0);
2536 W_Jac = JacCorLoc.T() * W;
2537
2538 CLHEP::HepVector gradL(6);
2539 gradL = W_Jac * V;
2540
2541 CLHEP::HepMatrix hessL(6, 6);
2542 hessL = JacCorLoc.T() * W * JacCorLoc;
2543
2544 GradL += gradL;
2545 HessL += hessL;
2546
2547 CLHEP::HepVector dLI = CLHEP::solve(HessL, -GradL);
2548 int iErI;
2549 CLHEP::HepMatrix ErrdLI = HessL.inverse(iErI);
2550 if (iErI != 0) {
2551 if (alarm)
2552 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " gradLocal Error inverse Hess matrix !!!!!!!!!!!";
2553 }
2554
2555 if (info || debug_) {
2556 edm::LogVerbatim("GlobalTrackerMuonAlignment")
2557 << " dL " << dLI(1) << " " << dLI(2) << " " << dLI(3) << " " << dLI(4) << " " << dLI(5) << " " << dLI(6);
2558 ;
2559 edm::LogVerbatim("GlobalTrackerMuonAlignment")
2560 << " +- " << std::sqrt(ErrdLI(1, 1)) << " " << std::sqrt(ErrdLI(2, 2)) << " " << std::sqrt(ErrdLI(3, 3))
2561 << " " << std::sqrt(ErrdLI(4, 4)) << " " << std::sqrt(ErrdLI(5, 5)) << " " << std::sqrt(ErrdLI(6, 6));
2562 }
2563
2564 if (debug_) {
2565 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug")
2566 << " dV(da3) {" << -JacCorLoc(1, 6) * 0.002 << " " << -JacCorLoc(2, 6) * 0.002 << " "
2567 << -JacCorLoc(3, 6) * 0.002 << " " << -JacCorLoc(4, 6) * 0.002 << "}";
2568 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " JacCLo {" << JacCorLoc(1, 6) << " " << JacCorLoc(2, 6)
2569 << " " << JacCorLoc(3, 6) << " " << JacCorLoc(4, 6) << "}";
2570 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug")
2571 << "Jpx/yx {" << Jac(4, 6) / Pm(3) << " " << Jac(5, 6) / Pm(3) << " " << Jac(1, 6) << " " << Jac(2, 6) << "}";
2572 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug")
2573 << "Jac(,a3){" << Jac(1, 6) << " " << Jac(2, 6) << " " << Jac(3, 6) << " " << Jac(4, 6) << " " << Jac(5, 6)
2574 << " " << Jac(6, 6);
2575 int i = 5;
2576 if (GNorm.z() > 0.95)
2577 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Ecap1 N " << GNorm;
2578 else if (GNorm.z() < -0.95)
2579 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Ecap2 N " << GNorm;
2580 else
2581 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Barrel N " << GNorm;
2582 edm::LogVerbatim("GlobalTrackerMuonAlignment")
2583 << " JacCLo(i," << i << ") = {" << JacCorLoc(1, i) << " " << JacCorLoc(2, i) << " " << JacCorLoc(3, i) << " "
2584 << JacCorLoc(4, i) << "}";
2585 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " rotLoc " << rotLoc;
2586 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " position " << R0;
2587 edm::LogVerbatim("GlobalTrackerMuonAlignment")
2588 << " Pm,l " << Pm.T() << " " << Pml(1) / Pml(3) << " " << Pml(2) / Pml(3);
2589 edm::LogVerbatim("GlobalTrackerMuonAlignment")
2590 << " Pt,l " << Pt.T() << " " << Ptl(1) / Ptl(3) << " " << Ptl(2) / Ptl(3);
2591 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " V " << V.T();
2592 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Cov \n" << covLoc;
2593 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " W*Cov " << W * covLoc;
2594
2595 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " JacToLoc " << JacToLoc;
2596 }
2597
2598 #ifdef CHECK_OF_JACOBIAN_CARTESIAN_TO_LOCAL
2599
2600 CLHEP::HepVector V0(4, 0);
2601 V0(1) = Pml(1) / Pml(3) - Ptl(1) / Ptl(3);
2602 V0(2) = Pml(2) / Pml(3) - Ptl(2) / Ptl(3);
2603 V0(3) = Rml(1) - Rtl(1);
2604 V0(4) = Rml(2) - Rtl(2);
2605 int ii = 3;
2606 float delta = 0.01;
2607 CLHEP::HepVector V1(4, 0);
2608 if (ii <= 3) {
2609 Rm(ii) += delta;
2610 Rml = rotLoc * (Rm - R0);
2611 } else {
2612 Pm(ii - 3) += delta;
2613 Pml = rotLoc * Pm;
2614 }
2615
2616
2617 V1(1) = Pml(1) / Pml(3) - Ptl(1) / Ptl(3);
2618 V1(2) = Pml(2) / Pml(3) - Ptl(2) / Ptl(3);
2619 V1(3) = Rml(1) - Rtl(1);
2620 V1(4) = Rml(2) - Rtl(2);
2621
2622 if (GNorm.z() > 0.95)
2623 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Ecap1 N " << GNorm;
2624 else if (GNorm.z() < -0.95)
2625 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Ecap2 N " << GNorm;
2626 else
2627 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Barrel N " << GNorm;
2628 edm::LogVerbatim("GlobalTrackerMuonAlignment")
2629 << " dldc Num (i," << ii << ") " << (V1(1) - V0(1)) / delta << " " << (V1(2) - V0(2)) / delta << " "
2630 << (V1(3) - V0(3)) / delta << " " << (V1(4) - V0(4)) / delta;
2631 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " dldc Ana (i," << ii << ") " << JacToLoc(1, ii) << " "
2632 << JacToLoc(2, ii) << " " << JacToLoc(3, ii) << " " << JacToLoc(4, ii);
2633 float dtxdpx = (rotLoc(1, 1) - rotLoc(3, 1) * Pml(1) / Pml(3)) / Pml(3);
2634 float dtydpx = (rotLoc(2, 1) - rotLoc(3, 2) * Pml(2) / Pml(3)) / Pml(3);
2635 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " dtx/dpx dty/ " << dtxdpx << " " << dtydpx;
2636 #endif
2637
2638 return;
2639 }
2640
2641
2642 void GlobalTrackerMuonAlignment::misalignMuon(
2643 GlobalVector& GRm, GlobalVector& GPm, GlobalVector& Nl, GlobalVector& Rt, GlobalVector& Rm, GlobalVector& Pm) {
2644 CLHEP::HepVector d(3, 0), a(3, 0), Norm(3), Rm0(3), Pm0(3), Rm1(3), Pm1(3), Rmi(3), Pmi(3);
2645
2646 d(1) = 0.0;
2647 d(2) = 0.0;
2648 d(3) = 0.0;
2649 a(1) = 0.0000;
2650 a(2) = 0.0000;
2651 a(3) = 0.0000;
2652
2653
2654
2655
2656
2657
2658
2659
2660
2661
2662
2663 MuGlShift = d;
2664 MuGlAngle = a;
2665
2666 if ((d(1) == 0.) && (d(2) == 0.) && (d(3) == 0.) && (a(1) == 0.) && (a(2) == 0.) && (a(3) == 0.)) {
2667 Rm = GRm;
2668 Pm = GPm;
2669 if (debug_)
2670 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ...... without misalignment ";
2671 return;
2672 }
2673
2674 Rm0(1) = GRm.x();
2675 Rm0(2) = GRm.y();
2676 Rm0(3) = GRm.z();
2677 Pm0(1) = GPm.x();
2678 Pm0(2) = GPm.y();
2679 Pm0(3) = GPm.z();
2680 Norm(1) = Nl.x();
2681 Norm(2) = Nl.y();
2682 Norm(3) = Nl.z();
2683
2684 CLHEP::HepMatrix T(3, 3, 1);
2685
2686
2687
2688
2689
2690 double s1 = std::sin(a(1)), c1 = std::cos(a(1));
2691 double s2 = std::sin(a(2)), c2 = std::cos(a(2));
2692 double s3 = std::sin(a(3)), c3 = std::cos(a(3));
2693 T(1, 1) = c2 * c3;
2694 T(1, 2) = c1 * s3 + s1 * s2 * c3;
2695 T(1, 3) = s1 * s3 - c1 * s2 * c3;
2696 T(2, 1) = -c2 * s3;
2697 T(2, 2) = c1 * c3 - s1 * s2 * s3;
2698 T(2, 3) = s1 * c3 + c1 * s2 * s3;
2699 T(3, 1) = s2;
2700 T(3, 2) = -s1 * c2;
2701 T(3, 3) = c1 * c2;
2702
2703
2704
2705 CLHEP::HepMatrix Ti = T.T();
2706 CLHEP::HepVector di = -Ti * d;
2707
2708 CLHEP::HepVector ex0(3, 0), ey0(3, 0), ez0(3, 0), ex(3, 0), ey(3, 0), ez(3, 0);
2709 ex0(1) = 1.;
2710 ey0(2) = 1.;
2711 ez0(3) = 1.;
2712 ex = Ti * ex0;
2713 ey = Ti * ey0;
2714 ez = Ti * ez0;
2715
2716 CLHEP::HepVector TiN = Ti * Norm;
2717
2718
2719
2720
2721 double si = CLHEP_dot((Ti * Rm0 - Rm0 + di), TiN) / CLHEP_dot(Pm0, TiN);
2722 Rm1(1) = CLHEP_dot(ex, Rm0 + si * Pm0 - di);
2723 Rm1(2) = CLHEP_dot(ey, Rm0 + si * Pm0 - di);
2724 Rm1(3) = CLHEP_dot(ez, Rm0 + si * Pm0 - di);
2725 Pm1 = T * Pm0;
2726
2727 Rm = GlobalVector(Rm1(1), Rm1(2), Rm1(3));
2728
2729 Pm = GlobalVector(CLHEP_dot(Pm0, ex), CLHEP_dot(Pm0, ey), CLHEP_dot(Pm0, ez));
2730
2731 if (debug_) {
2732
2733 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ----- Pm " << Pm;
2734 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " T*Pm0 " << Pm1.T();
2735 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " Ti*T*Pm0 " << (Ti * (T * Pm0)).T();
2736
2737
2738 CLHEP::HepVector Rml = Rm1(1) * ex + Rm1(2) * ey + Rm1(3) * ez + di;
2739
2740 CLHEP::HepVector Rt0(3);
2741 Rt0(1) = Rt.x();
2742 Rt0(2) = Rt.y();
2743 Rt0(3) = Rt.z();
2744
2745
2746 double sl = CLHEP_dot(T * (Rt0 - Rml), T * Norm) / CLHEP_dot(Ti * Pm1, Norm);
2747 CLHEP::HepVector Rl = Rml + sl * (Ti * Pm1);
2748
2749
2750
2751
2752 double A = CLHEP_dot(Ti * Pm1, Norm);
2753 double B = CLHEP_dot(Rt0, Norm) - CLHEP_dot(Ti * Rm1, Norm) + CLHEP_dot(Ti * d, Norm);
2754 double s = B / A;
2755 CLHEP::HepVector Dr = Ti * Rm1 - Ti * d + s * (Ti * Pm1) - Rm0;
2756 CLHEP::HepVector Dp = Ti * Pm1 - Pm0;
2757
2758 CLHEP::HepVector TiR = Ti * Rm0 + di;
2759 CLHEP::HepVector Ri = T * TiR + d;
2760
2761 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " --TTi-0 " << (Ri - Rm0).T();
2762 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -- Dr " << Dr.T();
2763 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -- Dp " << Dp.T();
2764 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ex " << ex.T();
2765 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ey " << ey.T();
2766 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ez " << ez.T();
2767 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ---- T ---- " << T;
2768 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ---- Ti ---- " << Ti;
2769 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ---- d ---- " << d.T();
2770 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ---- di ---- " << di.T();
2771 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -- si sl s " << si << " " << sl << " " << s;
2772 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- si sl " << si << " " << sl;
2773 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- si s " << si << " " << s;
2774 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -- Rml-(Rm0+si*Pm0) " << (Rml - Rm0 - si * Pm0).T();
2775 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- Rm0 " << Rm0.T();
2776 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- Rm1 " << Rm1.T();
2777 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- Rmi " << Rmi.T();
2778 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --siPm " << (si * Pm0).T();
2779 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --s2Pm " << (s2 * (T * Pm1)).T();
2780 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --R1-0 " << (Rm1 - Rm0).T();
2781 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --Ri-0 " << (Rmi - Rm0).T();
2782 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " --Rl-0 " << (Rl - Rm0).T();
2783 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --Pi-0 " << (Pmi - Pm0).T();
2784 }
2785
2786 return;
2787
2788 }
2789
2790
2791 void GlobalTrackerMuonAlignment::misalignMuonL(GlobalVector& GRm,
2792 GlobalVector& GPm,
2793 GlobalVector& Nl,
2794 GlobalVector& Rt,
2795 GlobalVector& Rm,
2796 GlobalVector& Pm,
2797 AlgebraicVector4& Vm,
2798 TrajectoryStateOnSurface& tsosTrack,
2799 TrajectoryStateOnSurface& tsosMuon) {
2800 CLHEP::HepVector d(3, 0), a(3, 0), Norm(3), Rm0(3), Pm0(3), Rm1(3), Pm1(3), Rmi(3), Pmi(3);
2801
2802 d(1) = 0.0;
2803 d(2) = 0.0;
2804 d(3) = 0.0;
2805 a(1) = 0.0000;
2806 a(2) = 0.0000;
2807 a(3) = 0.0000;
2808
2809
2810
2811
2812
2813
2814
2815
2816
2817
2818
2819
2820
2821 MuGlShift = d;
2822 MuGlAngle = a;
2823
2824 if ((d(1) == 0.) && (d(2) == 0.) && (d(3) == 0.) && (a(1) == 0.) && (a(2) == 0.) && (a(3) == 0.)) {
2825 Rm = GRm;
2826 Pm = GPm;
2827 AlgebraicVector4 Vm0;
2828 Vm = AlgebraicVector4(tsosMuon.localParameters().vector()(1),
2829 tsosMuon.localParameters().vector()(2),
2830 tsosMuon.localParameters().vector()(3),
2831 tsosMuon.localParameters().vector()(4));
2832 if (debug_)
2833 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ...... without misalignment ";
2834 return;
2835 }
2836
2837 Rm0(1) = GRm.x();
2838 Rm0(2) = GRm.y();
2839 Rm0(3) = GRm.z();
2840 Pm0(1) = GPm.x();
2841 Pm0(2) = GPm.y();
2842 Pm0(3) = GPm.z();
2843 Norm(1) = Nl.x();
2844 Norm(2) = Nl.y();
2845 Norm(3) = Nl.z();
2846
2847 CLHEP::HepMatrix T(3, 3, 1);
2848
2849
2850
2851
2852
2853 double s1 = std::sin(a(1)), c1 = std::cos(a(1));
2854 double s2 = std::sin(a(2)), c2 = std::cos(a(2));
2855 double s3 = std::sin(a(3)), c3 = std::cos(a(3));
2856 T(1, 1) = c2 * c3;
2857 T(1, 2) = c1 * s3 + s1 * s2 * c3;
2858 T(1, 3) = s1 * s3 - c1 * s2 * c3;
2859 T(2, 1) = -c2 * s3;
2860 T(2, 2) = c1 * c3 - s1 * s2 * s3;
2861 T(2, 3) = s1 * c3 + c1 * s2 * s3;
2862 T(3, 1) = s2;
2863 T(3, 2) = -s1 * c2;
2864 T(3, 3) = c1 * c2;
2865
2866
2867
2868
2869 CLHEP::HepMatrix Ti = T.T();
2870 CLHEP::HepVector di = -Ti * d;
2871
2872 CLHEP::HepVector ex0(3, 0), ey0(3, 0), ez0(3, 0), ex(3, 0), ey(3, 0), ez(3, 0);
2873 ex0(1) = 1.;
2874 ey0(2) = 1.;
2875 ez0(3) = 1.;
2876 ex = Ti * ex0;
2877 ey = Ti * ey0;
2878 ez = Ti * ez0;
2879
2880 CLHEP::HepVector TiN = Ti * Norm;
2881
2882
2883
2884
2885 double si = CLHEP_dot((Ti * Rm0 - Rm0 + di), TiN) / CLHEP_dot(Pm0, TiN);
2886 Rm1(1) = CLHEP_dot(ex, Rm0 + si * Pm0 - di);
2887 Rm1(2) = CLHEP_dot(ey, Rm0 + si * Pm0 - di);
2888 Rm1(3) = CLHEP_dot(ez, Rm0 + si * Pm0 - di);
2889 Pm1 = T * Pm0;
2890
2891
2892 Rm = GlobalVector(Rm1(1), Rm1(2), Rm1(3));
2893
2894 Pm = GlobalVector(CLHEP_dot(Pm0, ex), CLHEP_dot(Pm0, ey), CLHEP_dot(Pm0, ez));
2895
2896
2897 const Surface& refSurface = tsosMuon.surface();
2898 CLHEP::HepMatrix Tl(3, 3, 0);
2899
2900 Tl(1, 1) = refSurface.rotation().xx();
2901 Tl(1, 2) = refSurface.rotation().xy();
2902 Tl(1, 3) = refSurface.rotation().xz();
2903
2904 Tl(2, 1) = refSurface.rotation().yx();
2905 Tl(2, 2) = refSurface.rotation().yy();
2906 Tl(2, 3) = refSurface.rotation().yz();
2907
2908 Tl(3, 1) = refSurface.rotation().zx();
2909 Tl(3, 2) = refSurface.rotation().zy();
2910 Tl(3, 3) = refSurface.rotation().zz();
2911
2912 CLHEP::HepVector R0(3, 0), newR0(3, 0), newRl(3, 0), newPl(3, 0);
2913 R0(1) = refSurface.position().x();
2914 R0(2) = refSurface.position().y();
2915 R0(3) = refSurface.position().z();
2916
2917 newRl = Tl * (Rm1 - R0);
2918 newPl = Tl * Pm1;
2919
2920 Vm(0) = newPl(1) / newPl(3);
2921 Vm(1) = newPl(2) / newPl(3);
2922 Vm(2) = newRl(1);
2923 Vm(3) = newRl(2);
2924
2925 #ifdef CHECK_DERIVATIVES_LOCAL_VS_ANGLES
2926 float del = 0.0001;
2927
2928 int ii = 2;
2929 T(3, 1) -= del;
2930 T(1, 3) += del;
2931
2932 AlgebraicVector4 Vm0 = Vm;
2933 CLHEP::HepVector Rm10 = Rm1, Pm10 = Pm1;
2934
2935 CLHEP::HepMatrix newTl = Tl * T;
2936 Ti = T.T();
2937 di = -Ti * d;
2938 ex = Ti * ex0;
2939 ey = Ti * ey0;
2940 ez = Ti * ez0;
2941 TiN = Ti * Norm;
2942 si = CLHEP_dot((Ti * Rm0 - Rm0 + di), TiN) / CLHEP_dot(Pm0, TiN);
2943 Rm1(1) = CLHEP_dot(ex, Rm0 + si * Pm0 - di);
2944 Rm1(2) = CLHEP_dot(ey, Rm0 + si * Pm0 - di);
2945 Rm1(3) = CLHEP_dot(ez, Rm0 + si * Pm0 - di);
2946 Pm1 = T * Pm0;
2947
2948 newRl = Tl * (Rm1 - R0);
2949 newPl = Tl * Pm1;
2950
2951 Vm(0) = newPl(1) / newPl(3);
2952 Vm(1) = newPl(2) / newPl(3);
2953 Vm(2) = newRl(1);
2954 Vm(3) = newRl(2);
2955 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ========= dVm/da" << ii << " " << (Vm - Vm0) * (1. / del);
2956 if (debug_)
2957 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug")
2958 << " dRm/da3 " << ((Rm1 - Rm10) * (1. / del)).T() << " " << ((Pm1 - Pm10) * (1. / del)).T();
2959 #endif
2960
2961 if (debug_) {
2962 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " Norm " << Norm.T();
2963 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ex " << (Tl.T() * ex0).T();
2964 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ey " << (Tl.T() * ey0).T();
2965 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ez " << (Tl.T() * ez0).T();
2966
2967 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " pxyz/p gl 0 " << (Pm0 / Pm0.norm()).T();
2968 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " pxyz/p loc0 " << (Tl * Pm0 / (Tl * Pm0).norm()).T();
2969 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " pxyz/p glob " << (Pm1 / Pm1.norm()).T();
2970 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " pxyz/p loc " << (newPl / newPl.norm()).T();
2971
2972 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " Rot " << refSurface.rotation();
2973 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " Tl " << Tl;
2974 edm::LogVerbatim("GlobalTrackerMuonAlignment")
2975 << " ---- phi g0 g1 l0 l1 " << atan2(Pm0(2), Pm0(1)) << " " << atan2(Pm1(2), Pm1(1)) << " "
2976 << atan2((Tl * Pm0)(2), (Tl * Pm0)(1)) << " " << atan2(newPl(2), newPl(1));
2977 edm::LogVerbatim("GlobalTrackerMuonAlignment")
2978 << " ---- angl Gl Loc " << atan2(Pm1(2), Pm1(1)) - atan2(Pm0(2), Pm0(1)) << " "
2979 << atan2(newPl(2), newPl(1)) - atan2((Tl * Pm0)(2), (Tl * Pm0)(1)) << " "
2980 << atan2(newPl(3), newPl(2)) - atan2((Tl * Pm0)(3), (Tl * Pm0)(2)) << " "
2981 << atan2(newPl(1), newPl(3)) - atan2((Tl * Pm0)(1), (Tl * Pm0)(3)) << " ";
2982 }
2983
2984 if (debug_) {
2985 CLHEP::HepMatrix newTl(3, 3, 0);
2986 CLHEP::HepVector R0(3, 0), newRl(3, 0), newPl(3, 0);
2987 newTl = Tl * Ti.T();
2988 newR0 = Ti * R0 + di;
2989
2990 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " N " << Norm.T();
2991 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " dtxl yl " << Vm(0) - tsosMuon.localParameters().vector()(1)
2992 << " " << Vm(1) - tsosMuon.localParameters().vector()(2);
2993 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " dXl dYl " << Vm(2) - tsosMuon.localParameters().vector()(3)
2994 << " " << Vm(3) - tsosMuon.localParameters().vector()(4);
2995 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " localM " << tsosMuon.localParameters().vector();
2996 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Vm " << Vm;
2997
2998 CLHEP::HepVector Rvc(3, 0), Pvc(3, 0), Rvg(3, 0), Pvg(3, 0);
2999 Rvc(1) = Vm(2);
3000 Rvc(2) = Vm(3);
3001 Pvc(3) = tsosMuon.localParameters().pzSign() * Pm0.norm() / std::sqrt(1 + Vm(0) * Vm(0) + Vm(1) * Vm(1));
3002 Pvc(1) = Pvc(3) * Vm(0);
3003 Pvc(2) = Pvc(3) * Vm(1);
3004
3005 Rvg = newTl.T() * Rvc + newR0;
3006 Pvg = newTl.T() * Pvc;
3007
3008 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " newPl " << newPl.T();
3009 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Pvc " << Pvc.T();
3010 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Rvg " << Rvg.T();
3011 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Rm1 " << Rm1.T();
3012 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Pvg " << Pvg.T();
3013 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Pm1 " << Pm1.T();
3014 }
3015
3016 if (debug_) {
3017
3018 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ----- Pm " << Pm;
3019 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " T*Pm0 " << Pm1.T();
3020 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " Ti*T*Pm0 " << (Ti * (T * Pm0)).T();
3021
3022
3023 CLHEP::HepVector Rml = Rm1(1) * ex + Rm1(2) * ey + Rm1(3) * ez + di;
3024
3025 CLHEP::HepVector Rt0(3);
3026 Rt0(1) = Rt.x();
3027 Rt0(2) = Rt.y();
3028 Rt0(3) = Rt.z();
3029
3030
3031 double sl = CLHEP_dot(T * (Rt0 - Rml), T * Norm) / CLHEP_dot(Ti * Pm1, Norm);
3032 CLHEP::HepVector Rl = Rml + sl * (Ti * Pm1);
3033
3034
3035
3036
3037 double A = CLHEP_dot(Ti * Pm1, Norm);
3038 double B = CLHEP_dot(Rt0, Norm) - CLHEP_dot(Ti * Rm1, Norm) + CLHEP_dot(Ti * d, Norm);
3039 double s = B / A;
3040 CLHEP::HepVector Dr = Ti * Rm1 - Ti * d + s * (Ti * Pm1) - Rm0;
3041 CLHEP::HepVector Dp = Ti * Pm1 - Pm0;
3042
3043 CLHEP::HepVector TiR = Ti * Rm0 + di;
3044 CLHEP::HepVector Ri = T * TiR + d;
3045
3046 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " --TTi-0 " << (Ri - Rm0).T();
3047 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -- Dr " << Dr.T();
3048 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -- Dp " << Dp.T();
3049 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ex " << ex.T();
3050 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ey " << ey.T();
3051 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ez " << ez.T();
3052 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ---- T ---- " << T;
3053 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ---- Ti ---- " << Ti;
3054 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ---- d ---- " << d.T();
3055 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ---- di ---- " << di.T();
3056 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -- si sl s " << si << " " << sl << " " << s;
3057 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- si sl " << si << " " << sl;
3058 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- si s " << si << " " << s;
3059 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -- Rml-(Rm0+si*Pm0) " << (Rml - Rm0 - si * Pm0).T();
3060 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- Rm0 " << Rm0.T();
3061 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- Rm1 " << Rm1.T();
3062 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- Rmi " << Rmi.T();
3063 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --siPm " << (si * Pm0).T();
3064 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --s2Pm " << (s2 * (T * Pm1)).T();
3065 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --R1-0 " << (Rm1 - Rm0).T();
3066 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --Ri-0 " << (Rmi - Rm0).T();
3067 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " --Rl-0 " << (Rl - Rm0).T();
3068 edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --Pi-0 " << (Pmi - Pm0).T();
3069 }
3070
3071 return;
3072
3073 }
3074
3075
3076 void GlobalTrackerMuonAlignment::trackFitter(reco::TrackRef alongTr,
3077 reco::TransientTrack& alongTTr,
3078 PropagationDirection direction,
3079 TrajectoryStateOnSurface& trackFittedTSOS) {
3080 bool info = false;
3081 bool smooth = false;
3082
3083 if (info)
3084 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ......... start of trackFitter ......... ";
3085 if (false && info) {
3086 this->debugTrackHit(" Tracker track hits ", alongTr);
3087 this->debugTrackHit(" Tracker TransientTrack hits ", alongTTr);
3088 }
3089
3090 edm::OwnVector<TrackingRecHit> recHit;
3091 DetId trackDetId = DetId(alongTr->innerDetId());
3092 for (auto const& hit : alongTr->recHits())
3093 recHit.push_back(hit->clone());
3094 if (direction != alongMomentum) {
3095 edm::OwnVector<TrackingRecHit> recHitAlong = recHit;
3096 recHit.clear();
3097 for (unsigned int ihit = recHitAlong.size() - 1; ihit + 1 > 0; ihit--) {
3098 recHit.push_back(recHitAlong[ihit]);
3099 }
3100 recHitAlong.clear();
3101 }
3102 if (info)
3103 edm::LogVerbatim("GlobalTrackerMuonAlignment")
3104 << " ... Own recHit.size() " << recHit.size() << " " << trackDetId.rawId();
3105
3106
3107 TransientTrackingRecHit::RecHitContainer recHitTrack;
3108 for (auto const& hit : alongTr->recHits())
3109 recHitTrack.push_back(TTRHBuilder->build(hit));
3110
3111 if (info)
3112 edm::LogVerbatim("GlobalTrackerMuonAlignment")
3113 << " ... recHitTrack.size() " << recHit.size() << " " << trackDetId.rawId() << " ======> recHitMu ";
3114
3115 MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu = recHitTrack;
3116
3117
3118
3119
3120 if (info)
3121 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ...along.... recHitMu.size() " << recHitMu.size();
3122 if (direction != alongMomentum) {
3123 MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMuAlong = recHitMu;
3124 recHitMu.clear();
3125 for (unsigned int ihit = recHitMuAlong.size() - 1; ihit + 1 > 0; ihit--) {
3126 recHitMu.push_back(recHitMuAlong[ihit]);
3127 }
3128 recHitMuAlong.clear();
3129 }
3130 if (info)
3131 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ...opposite... recHitMu.size() " << recHitMu.size();
3132
3133 TrajectoryStateOnSurface firstTSOS;
3134 if (direction == alongMomentum)
3135 firstTSOS = alongTTr.innermostMeasurementState();
3136 else
3137 firstTSOS = alongTTr.outermostMeasurementState();
3138
3139 AlgebraicSymMatrix55 CovLoc;
3140 CovLoc(0, 0) = 1. * firstTSOS.localError().matrix()(0, 0);
3141 CovLoc(1, 1) = 10. * firstTSOS.localError().matrix()(1, 1);
3142 CovLoc(2, 2) = 10. * firstTSOS.localError().matrix()(2, 2);
3143 CovLoc(3, 3) = 100. * firstTSOS.localError().matrix()(3, 3);
3144 CovLoc(4, 4) = 100. * firstTSOS.localError().matrix()(4, 4);
3145 TrajectoryStateOnSurface initialTSOS(
3146 firstTSOS.localParameters(), LocalTrajectoryError(CovLoc), firstTSOS.surface(), &*magneticField_);
3147
3148 PTrajectoryStateOnDet PTraj = trajectoryStateTransform::persistentState(initialTSOS, trackDetId.rawId());
3149
3150 const TrajectorySeed seedT(PTraj, recHit, direction);
3151
3152 std::vector<Trajectory> trajVec;
3153 if (direction == alongMomentum)
3154 trajVec = theFitter->fit(seedT, recHitMu, initialTSOS);
3155 else
3156 trajVec = theFitterOp->fit(seedT, recHitMu, initialTSOS);
3157
3158 if (info) {
3159 this->debugTrajectorySOSv(" innermostTSOS of TransientTrack", alongTTr.innermostMeasurementState());
3160 this->debugTrajectorySOSv(" outermostTSOS of TransientTrack", alongTTr.outermostMeasurementState());
3161 this->debugTrajectorySOS(" initialTSOS for theFitter ", initialTSOS);
3162 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " . . . . . trajVec.size() " << trajVec.size();
3163 if (!trajVec.empty())
3164 this->debugTrajectory(" theFitter trajectory ", trajVec[0]);
3165 }
3166
3167 if (!smooth) {
3168 if (!trajVec.empty())
3169 trackFittedTSOS = trajVec[0].lastMeasurement().updatedState();
3170 } else {
3171 std::vector<Trajectory> trajSVec;
3172 if (!trajVec.empty()) {
3173 if (direction == alongMomentum)
3174 trajSVec = theSmoother->trajectories(*(trajVec.begin()));
3175 else
3176 trajSVec = theSmootherOp->trajectories(*(trajVec.begin()));
3177 }
3178 if (info)
3179 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " . . . . trajSVec.size() " << trajSVec.size();
3180 if (!trajSVec.empty())
3181 this->debugTrajectorySOSv("smoothed geom InnermostState", trajSVec[0].geometricalInnermostState());
3182 if (!trajSVec.empty())
3183 trackFittedTSOS = trajSVec[0].geometricalInnermostState();
3184 }
3185
3186 if (info)
3187 this->debugTrajectorySOSv(" trackFittedTSOS ", trackFittedTSOS);
3188
3189 }
3190
3191
3192 void GlobalTrackerMuonAlignment::muonFitter(reco::TrackRef alongTr,
3193 reco::TransientTrack& alongTTr,
3194 PropagationDirection direction,
3195 TrajectoryStateOnSurface& trackFittedTSOS) {
3196 bool info = false;
3197 bool smooth = false;
3198
3199 if (info)
3200 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ......... start of muonFitter ........";
3201 if (false && info) {
3202 this->debugTrackHit(" Muon track hits ", alongTr);
3203 this->debugTrackHit(" Muon TransientTrack hits ", alongTTr);
3204 }
3205
3206 edm::OwnVector<TrackingRecHit> recHit;
3207 DetId trackDetId = DetId(alongTr->innerDetId());
3208 for (auto const& hit : alongTr->recHits())
3209 recHit.push_back(hit->clone());
3210 if (direction != alongMomentum) {
3211 edm::OwnVector<TrackingRecHit> recHitAlong = recHit;
3212 recHit.clear();
3213 for (unsigned int ihit = recHitAlong.size() - 1; ihit + 1 > 0; ihit--) {
3214 recHit.push_back(recHitAlong[ihit]);
3215 }
3216 recHitAlong.clear();
3217 }
3218 if (info)
3219 edm::LogVerbatim("GlobalTrackerMuonAlignment")
3220 << " ... Own recHit.size() DetId==Muon " << recHit.size() << " " << trackDetId.rawId();
3221
3222 MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu =
3223 MuRHBuilder->build(alongTTr.recHitsBegin(), alongTTr.recHitsEnd());
3224 if (info)
3225 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ...along.... recHitMu.size() " << recHitMu.size();
3226 if (direction != alongMomentum) {
3227 MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMuAlong = recHitMu;
3228 recHitMu.clear();
3229 for (unsigned int ihit = recHitMuAlong.size() - 1; ihit + 1 > 0; ihit--) {
3230 recHitMu.push_back(recHitMuAlong[ihit]);
3231 }
3232 recHitMuAlong.clear();
3233 }
3234 if (info)
3235 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ...opposite... recHitMu.size() " << recHitMu.size();
3236
3237 TrajectoryStateOnSurface firstTSOS;
3238 if (direction == alongMomentum)
3239 firstTSOS = alongTTr.innermostMeasurementState();
3240 else
3241 firstTSOS = alongTTr.outermostMeasurementState();
3242
3243 AlgebraicSymMatrix55 CovLoc;
3244 CovLoc(0, 0) = 1. * firstTSOS.localError().matrix()(0, 0);
3245 CovLoc(1, 1) = 10. * firstTSOS.localError().matrix()(1, 1);
3246 CovLoc(2, 2) = 10. * firstTSOS.localError().matrix()(2, 2);
3247 CovLoc(3, 3) = 100. * firstTSOS.localError().matrix()(3, 3);
3248 CovLoc(4, 4) = 100. * firstTSOS.localError().matrix()(4, 4);
3249 TrajectoryStateOnSurface initialTSOS(
3250 firstTSOS.localParameters(), LocalTrajectoryError(CovLoc), firstTSOS.surface(), &*magneticField_);
3251
3252 PTrajectoryStateOnDet PTraj = trajectoryStateTransform::persistentState(initialTSOS, trackDetId.rawId());
3253
3254 const TrajectorySeed seedT(PTraj, recHit, direction);
3255
3256 std::vector<Trajectory> trajVec;
3257 if (direction == alongMomentum)
3258 trajVec = theFitter->fit(seedT, recHitMu, initialTSOS);
3259 else
3260 trajVec = theFitterOp->fit(seedT, recHitMu, initialTSOS);
3261
3262 if (info) {
3263 this->debugTrajectorySOSv(" innermostTSOS of TransientTrack", alongTTr.innermostMeasurementState());
3264 this->debugTrajectorySOSv(" outermostTSOS of TransientTrack", alongTTr.outermostMeasurementState());
3265 this->debugTrajectorySOS(" initialTSOS for theFitter ", initialTSOS);
3266 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " . . . . . trajVec.size() " << trajVec.size();
3267 if (!trajVec.empty())
3268 this->debugTrajectory(" theFitter trajectory ", trajVec[0]);
3269 }
3270
3271 if (!smooth) {
3272 if (!trajVec.empty())
3273 trackFittedTSOS = trajVec[0].lastMeasurement().updatedState();
3274 } else {
3275 std::vector<Trajectory> trajSVec;
3276 if (!trajVec.empty()) {
3277 if (direction == alongMomentum)
3278 trajSVec = theSmoother->trajectories(*(trajVec.begin()));
3279 else
3280 trajSVec = theSmootherOp->trajectories(*(trajVec.begin()));
3281 }
3282 if (info)
3283 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " . . . . trajSVec.size() " << trajSVec.size();
3284 if (!trajSVec.empty())
3285 this->debugTrajectorySOSv("smoothed geom InnermostState", trajSVec[0].geometricalInnermostState());
3286 if (!trajSVec.empty())
3287 trackFittedTSOS = trajSVec[0].geometricalInnermostState();
3288 }
3289
3290 }
3291
3292
3293 void GlobalTrackerMuonAlignment::debugTrackHit(const std::string title, TransientTrack& alongTr) {
3294 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ------- " << title << " --------";
3295 int nHit = 1;
3296 for (trackingRecHit_iterator i = alongTr.recHitsBegin(); i != alongTr.recHitsEnd(); i++) {
3297 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Hit " << nHit++ << " DetId " << (*i)->geographicalId().det();
3298 if ((*i)->geographicalId().det() == DetId::Tracker)
3299 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Tracker ";
3300 else if ((*i)->geographicalId().det() == DetId::Muon)
3301 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Muon ";
3302 else
3303 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Unknown ";
3304 if (!(*i)->isValid())
3305 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " not valid ";
3306 else
3307 edm::LogVerbatim("GlobalTrackerMuonAlignment");
3308 }
3309 }
3310
3311
3312 void GlobalTrackerMuonAlignment::debugTrackHit(const std::string title, reco::TrackRef alongTr) {
3313 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ------- " << title << " --------";
3314 int nHit = 1;
3315 for (auto const& hit : alongTr->recHits()) {
3316 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Hit " << nHit++ << " DetId " << hit->geographicalId().det();
3317 if (hit->geographicalId().det() == DetId::Tracker)
3318 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Tracker ";
3319 else if (hit->geographicalId().det() == DetId::Muon)
3320 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Muon ";
3321 else
3322 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Unknown ";
3323 if (!hit->isValid())
3324 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " not valid ";
3325 else
3326 edm::LogVerbatim("GlobalTrackerMuonAlignment");
3327 }
3328 }
3329
3330
3331 void GlobalTrackerMuonAlignment::debugTrajectorySOS(const std::string title, TrajectoryStateOnSurface& trajSOS) {
3332 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " --- " << title << " --- ";
3333 if (!trajSOS.isValid()) {
3334 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Not valid !!!! ";
3335 return;
3336 }
3337 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " R |p| " << trajSOS.globalPosition().perp() << " "
3338 << trajSOS.globalMomentum().mag() << " charge " << trajSOS.charge();
3339 edm::LogVerbatim("GlobalTrackerMuonAlignment")
3340 << " x p " << trajSOS.globalParameters().vector()(0) << " " << trajSOS.globalParameters().vector()(1) << " "
3341 << trajSOS.globalParameters().vector()(2) << " " << trajSOS.globalParameters().vector()(3) << " "
3342 << trajSOS.globalParameters().vector()(4) << " " << trajSOS.globalParameters().vector()(5);
3343 edm::LogVerbatim("GlobalTrackerMuonAlignment")
3344 << " +/- " << std::sqrt(trajSOS.cartesianError().matrix()(0, 0)) << " "
3345 << std::sqrt(trajSOS.cartesianError().matrix()(1, 1)) << " " << std::sqrt(trajSOS.cartesianError().matrix()(2, 2))
3346 << " " << std::sqrt(trajSOS.cartesianError().matrix()(3, 3)) << " "
3347 << std::sqrt(trajSOS.cartesianError().matrix()(4, 4)) << " "
3348 << std::sqrt(trajSOS.cartesianError().matrix()(5, 5));
3349 edm::LogVerbatim("GlobalTrackerMuonAlignment")
3350 << " q/p dxy/dz xy " << trajSOS.localParameters().vector()(0) << " " << trajSOS.localParameters().vector()(1)
3351 << " " << trajSOS.localParameters().vector()(2) << " " << trajSOS.localParameters().vector()(3) << " "
3352 << trajSOS.localParameters().vector()(4);
3353 edm::LogVerbatim("GlobalTrackerMuonAlignment")
3354 << " +/- error " << std::sqrt(trajSOS.localError().matrix()(0, 0)) << " "
3355 << std::sqrt(trajSOS.localError().matrix()(1, 1)) << " " << std::sqrt(trajSOS.localError().matrix()(2, 2)) << " "
3356 << std::sqrt(trajSOS.localError().matrix()(3, 3)) << " " << std::sqrt(trajSOS.localError().matrix()(4, 4)) << " ";
3357 }
3358
3359
3360 void GlobalTrackerMuonAlignment::debugTrajectorySOSv(const std::string title, TrajectoryStateOnSurface trajSOS) {
3361 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " --- " << title << " --- ";
3362 if (!trajSOS.isValid()) {
3363 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Not valid !!!! ";
3364 return;
3365 }
3366 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " R |p| " << trajSOS.globalPosition().perp() << " "
3367 << trajSOS.globalMomentum().mag() << " charge " << trajSOS.charge();
3368 edm::LogVerbatim("GlobalTrackerMuonAlignment")
3369 << " x p " << trajSOS.globalParameters().vector()(0) << " " << trajSOS.globalParameters().vector()(1) << " "
3370 << trajSOS.globalParameters().vector()(2) << " " << trajSOS.globalParameters().vector()(3) << " "
3371 << trajSOS.globalParameters().vector()(4) << " " << trajSOS.globalParameters().vector()(5);
3372 edm::LogVerbatim("GlobalTrackerMuonAlignment")
3373 << " +/- " << std::sqrt(trajSOS.cartesianError().matrix()(0, 0)) << " "
3374 << std::sqrt(trajSOS.cartesianError().matrix()(1, 1)) << " " << std::sqrt(trajSOS.cartesianError().matrix()(2, 2))
3375 << " " << std::sqrt(trajSOS.cartesianError().matrix()(3, 3)) << " "
3376 << std::sqrt(trajSOS.cartesianError().matrix()(4, 4)) << " "
3377 << std::sqrt(trajSOS.cartesianError().matrix()(5, 5));
3378 edm::LogVerbatim("GlobalTrackerMuonAlignment")
3379 << " q/p dxy/dz xy " << trajSOS.localParameters().vector()(0) << " " << trajSOS.localParameters().vector()(1)
3380 << " " << trajSOS.localParameters().vector()(2) << " " << trajSOS.localParameters().vector()(3) << " "
3381 << trajSOS.localParameters().vector()(4);
3382 edm::LogVerbatim("GlobalTrackerMuonAlignment")
3383 << " +/- error " << std::sqrt(trajSOS.localError().matrix()(0, 0)) << " "
3384 << std::sqrt(trajSOS.localError().matrix()(1, 1)) << " " << std::sqrt(trajSOS.localError().matrix()(2, 2)) << " "
3385 << std::sqrt(trajSOS.localError().matrix()(3, 3)) << " " << std::sqrt(trajSOS.localError().matrix()(4, 4)) << " ";
3386 }
3387
3388
3389 void GlobalTrackerMuonAlignment::debugTrajectory(const std::string title, Trajectory& traj) {
3390 edm::LogVerbatim("GlobalTrackerMuonAlignment") << "\n"
3391 << " ...... " << title << " ...... ";
3392 if (!traj.isValid()) {
3393 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Not valid !!!!!!!! ";
3394 return;
3395 }
3396 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " chi2/Nhit " << traj.chiSquared() << " / " << traj.foundHits();
3397 if (traj.direction() == alongMomentum)
3398 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " alongMomentum >>>>";
3399 else
3400 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " oppositeToMomentum <<<<";
3401 this->debugTrajectorySOSv(" firstMeasurementTSOS ", traj.firstMeasurement().updatedState());
3402
3403 this->debugTrajectorySOSv(" lastMeasurementTSOS ", traj.lastMeasurement().updatedState());
3404
3405 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " . . . . . . . . . . . . . . . . . . . . . . . . . . . . \n";
3406 }
3407
3408
3409 void GlobalTrackerMuonAlignment::writeGlPosRcd(CLHEP::HepVector& paramVec) {
3410
3411
3412
3413
3414
3415
3416
3417
3418
3419
3420
3421 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " paramVector " << paramVec.T();
3422
3423 CLHEP::Hep3Vector colX, colY, colZ;
3424
3425 #ifdef NOT_EXACT_ROTATION_MATRIX
3426 colX = CLHEP::Hep3Vector(1., -paramVec(6), paramVec(5));
3427 colY = CLHEP::Hep3Vector(paramVec(6), 1., -paramVec(4));
3428 colZ = CLHEP::Hep3Vector(-paramVec(5), paramVec(4), 1.);
3429 #else
3430 double s1 = std::sin(paramVec(4)), c1 = std::cos(paramVec(4));
3431 double s2 = std::sin(paramVec(5)), c2 = std::cos(paramVec(5));
3432 double s3 = std::sin(paramVec(6)), c3 = std::cos(paramVec(6));
3433 colX = CLHEP::Hep3Vector(c2 * c3, -c2 * s3, s2);
3434 colY = CLHEP::Hep3Vector(c1 * s3 + s1 * s2 * c3, c1 * c3 - s1 * s2 * s3, -s1 * c2);
3435 colZ = CLHEP::Hep3Vector(s1 * s3 - c1 * s2 * c3, s1 * c3 + c1 * s2 * s3, c1 * c2);
3436 #endif
3437
3438 CLHEP::HepVector param0(6, 0);
3439
3440 Alignments globalPositions{};
3441
3442
3443 AlignTransform tracker(
3444 iteratorTrackerRcd->translation(), iteratorTrackerRcd->rotation(), DetId(DetId::Tracker).rawId());
3445
3446 CLHEP::Hep3Vector posMuGlRcd = iteratorMuonRcd->translation();
3447 CLHEP::HepRotation rotMuGlRcd = iteratorMuonRcd->rotation();
3448 CLHEP::HepEulerAngles angMuGlRcd = iteratorMuonRcd->rotation().eulerAngles();
3449 if (debug_)
3450 edm::LogVerbatim("GlobalTrackerMuonAlignment")
3451 << " Old muon Rcd Euler angles " << angMuGlRcd.phi() << " " << angMuGlRcd.theta() << " " << angMuGlRcd.psi();
3452 AlignTransform muon;
3453 if ((angMuGlRcd.phi() == 0.) && (angMuGlRcd.theta() == 0.) && (angMuGlRcd.psi() == 0.) && (posMuGlRcd.x() == 0.) &&
3454 (posMuGlRcd.y() == 0.) && (posMuGlRcd.z() == 0.)) {
3455 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " New muon parameters are stored in Rcd";
3456
3457 AlignTransform muonNew(AlignTransform::Translation(paramVec(1), paramVec(2), paramVec(3)),
3458 AlignTransform::Rotation(colX, colY, colZ),
3459 DetId(DetId::Muon).rawId());
3460 muon = muonNew;
3461 } else if ((paramVec(1) == 0.) && (paramVec(2) == 0.) && (paramVec(3) == 0.) && (paramVec(4) == 0.) &&
3462 (paramVec(5) == 0.) && (paramVec(6) == 0.)) {
3463 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Old muon parameters are stored in Rcd";
3464
3465 AlignTransform muonNew(iteratorMuonRcd->translation(), iteratorMuonRcd->rotation(), DetId(DetId::Muon).rawId());
3466 muon = muonNew;
3467 } else {
3468 edm::LogVerbatim("GlobalTrackerMuonAlignment") << " New + Old muon parameters are stored in Rcd";
3469 CLHEP::Hep3Vector posMuGlRcdThis = CLHEP::Hep3Vector(paramVec(1), paramVec(2), paramVec(3));
3470 CLHEP::HepRotation rotMuGlRcdThis = CLHEP::HepRotation(colX, colY, colZ);
3471 CLHEP::Hep3Vector posMuGlRcdNew = rotMuGlRcdThis * posMuGlRcd + posMuGlRcdThis;
3472 CLHEP::HepRotation rotMuGlRcdNew = rotMuGlRcdThis * rotMuGlRcd;
3473
3474 AlignTransform muonNew(posMuGlRcdNew, rotMuGlRcdNew, DetId(DetId::Muon).rawId());
3475 muon = muonNew;
3476 }
3477
3478
3479 AlignTransform ecal(iteratorEcalRcd->translation(), iteratorEcalRcd->rotation(), DetId(DetId::Ecal).rawId());
3480
3481 AlignTransform hcal(iteratorHcalRcd->translation(), iteratorHcalRcd->rotation(), DetId(DetId::Hcal).rawId());
3482
3483 AlignTransform calo(AlignTransform::Translation(param0(1), param0(2), param0(3)),
3484 AlignTransform::EulerAngles(param0(4), param0(5), param0(6)),
3485 DetId(DetId::Calo).rawId());
3486
3487 edm::LogVerbatim("GlobalTrackerMuonAlignment")
3488 << "Tracker (" << tracker.rawId() << ") at " << tracker.translation() << " " << tracker.rotation().eulerAngles();
3489 edm::LogVerbatim("GlobalTrackerMuonAlignment") << tracker.rotation();
3490
3491 edm::LogVerbatim("GlobalTrackerMuonAlignment")
3492 << "Muon (" << muon.rawId() << ") at " << muon.translation() << " " << muon.rotation().eulerAngles();
3493 edm::LogVerbatim("GlobalTrackerMuonAlignment")
3494 << " rotations angles around x,y,z "
3495 << " ( " << -muon.rotation().zy() << " " << muon.rotation().zx() << " " << -muon.rotation().yx() << " )";
3496 edm::LogVerbatim("GlobalTrackerMuonAlignment") << muon.rotation();
3497
3498 edm::LogVerbatim("GlobalTrackerMuonAlignment")
3499 << "Ecal (" << ecal.rawId() << ") at " << ecal.translation() << " " << ecal.rotation().eulerAngles();
3500 edm::LogVerbatim("GlobalTrackerMuonAlignment") << ecal.rotation();
3501
3502 edm::LogVerbatim("GlobalTrackerMuonAlignment")
3503 << "Hcal (" << hcal.rawId() << ") at " << hcal.translation() << " " << hcal.rotation().eulerAngles();
3504 edm::LogVerbatim("GlobalTrackerMuonAlignment") << hcal.rotation();
3505
3506 edm::LogVerbatim("GlobalTrackerMuonAlignment")
3507 << "Calo (" << calo.rawId() << ") at " << calo.translation() << " " << calo.rotation().eulerAngles();
3508 edm::LogVerbatim("GlobalTrackerMuonAlignment") << calo.rotation();
3509
3510 globalPositions.m_align.push_back(tracker);
3511 globalPositions.m_align.push_back(muon);
3512 globalPositions.m_align.push_back(ecal);
3513 globalPositions.m_align.push_back(hcal);
3514 globalPositions.m_align.push_back(calo);
3515
3516 edm::LogVerbatim("GlobalTrackerMuonAlignment") << "Uploading to the database...";
3517
3518 edm::Service<cond::service::PoolDBOutputService> poolDbService;
3519
3520 if (!poolDbService.isAvailable())
3521 throw cms::Exception("NotAvailable") << "PoolDBOutputService not available";
3522
3523
3524
3525
3526
3527
3528 poolDbService->writeOneIOV<Alignments>(globalPositions, poolDbService->currentTime(), "GlobalPositionRcd");
3529 edm::LogVerbatim("GlobalTrackerMuonAlignment") << "done!";
3530
3531 return;
3532 }
3533
3534
3535 DEFINE_FWK_MODULE(GlobalTrackerMuonAlignment);