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