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