CMS 3D CMS Logo
GlobalTrackerMuonAlignment.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: GlobalTrackerMuonAlignment
4 // Class: GlobalTrackerMuonAlignment
5 //
15 //
16 // Original Author: Alexandre Spiridonov
17 // Created: Fri Oct 16 15:59:05 CEST 2009
18 //
19 // $Id: GlobalTrackerMuonAlignment.cc,v 1.11 2012/07/16 12:17:53 eulisse Exp $
20 //
21 
22 // system include files
23 #include <memory>
24 #include <string>
25 #include <map>
26 #include <vector>
27 #include <fstream>
28 #include <typeinfo>
29 
30 // user include files
31 
32 // Framework
33 #include "FWCore/MessageLogger/interface/MessageLogger.h"
34 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
35 #include "FWCore/Framework/interface/ESWatcher.h"
36 #include "FWCore/Framework/interface/Frameworkfwd.h"
37 #include "FWCore/Framework/interface/Event.h"
38 #include "FWCore/Framework/interface/EventSetup.h"
39 #include "FWCore/Framework/interface/MakerMacros.h"
40 #include "FWCore/Utilities/interface/InputTag.h"
41 #include "FWCore/ParameterSet/interface/ParameterSet.h"
42 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
43 #include "FWCore/ServiceRegistry/interface/Service.h"
44 //#include "FWCore/Framework/interface/EventSetupRecordImplementation.h"
45 
46 // references
47 #include "DataFormats/DetId/interface/DetId.h"
48 #include "DataFormats/MuonReco/interface/Muon.h"
49 #include "DataFormats/MuonReco/interface/MuonFwd.h"
50 #include "DataFormats/TrackReco/interface/Track.h"
51 #include "DataFormats/TrackReco/interface/TrackFwd.h"
52 #include "DataFormats/GeometrySurface/interface/TangentPlane.h"
53 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
54 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
55 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
56 #include "MagneticField/Engine/interface/MagneticField.h"
57 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
58 
59 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBuilder.h"
60 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
61 #include "TrackingTools/TrackRefitter/interface/TrackTransformer.h"
62 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
63 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
64 #include "TrackingTools/GeomPropagators/interface/SmartPropagator.h"
65 #include "TrackingTools/GeomPropagators/interface/PropagationDirectionChooser.h"
66 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
67 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
68 #include "TrackPropagation/RungeKutta/interface/defaultRKPropagator.h"
69 #include "TrackingTools/PatternTools/interface/Trajectory.h"
70 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
71 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCartesianToLocal.h"
72 
73 // Database
74 #include "FWCore/Framework/interface/ESHandle.h"
75 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
76 
77 // Alignment
78 #include "CondFormats/Alignment/interface/Alignments.h"
79 #include "CondFormats/Alignment/interface/AlignTransform.h"
80 #include "CondFormats/AlignmentRecord/interface/GlobalPositionRcd.h"
81 
82 // Refit
83 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
84 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
85 #include "TrackingTools/TrackFitters/interface/KFTrajectoryFitter.h"
86 #include "TrackingTools/TrackFitters/interface/KFTrajectorySmoother.h"
87 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
88 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
89 
90 #include <TH1.h>
91 #include <TH2.h>
92 #include <TFile.h>
93 
94 #include <iostream>
95 #include <cmath>
96 #include <cassert>
97 #include <fstream>
98 #include <iomanip>
99 
100 using namespace edm;
101 using namespace std;
102 using namespace reco;
103 
104 //
105 // class declaration
106 //
107 
108 class GlobalTrackerMuonAlignment : public edm::one::EDAnalyzer<> {
109 public:
110  explicit GlobalTrackerMuonAlignment(const edm::ParameterSet&);
111  ~GlobalTrackerMuonAlignment() override = default;
112  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
113 
114  void analyzeTrackTrack(const edm::Event&, const edm::EventSetup&);
115  void analyzeTrackTrajectory(const edm::Event&, const edm::EventSetup&);
116  void bookHist();
117  void fitHist();
118  void trackFitter(reco::TrackRef, reco::TransientTrack&, PropagationDirection, TrajectoryStateOnSurface&);
119  void muonFitter(reco::TrackRef, reco::TransientTrack&, PropagationDirection, TrajectoryStateOnSurface&);
120  void debugTrackHit(const std::string, reco::TrackRef);
121  void debugTrackHit(const std::string, reco::TransientTrack&);
122  void debugTrajectorySOS(const std::string, TrajectoryStateOnSurface&);
123  void debugTrajectorySOSv(const std::string, TrajectoryStateOnSurface);
124  void debugTrajectory(const std::string, Trajectory&);
125 
126  void gradientGlobal(GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&, AlgebraicSymMatrix66&);
127  void gradientLocal(GlobalVector&,
128  GlobalVector&,
129  GlobalVector&,
130  GlobalVector&,
131  GlobalVector&,
132  CLHEP::HepSymMatrix&,
133  CLHEP::HepMatrix&,
134  CLHEP::HepVector&,
135  AlgebraicVector4&);
136  void gradientGlobalAlg(GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&, AlgebraicSymMatrix66&);
137  void misalignMuon(GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&);
138  void misalignMuonL(GlobalVector&,
139  GlobalVector&,
140  GlobalVector&,
141  GlobalVector&,
142  GlobalVector&,
143  GlobalVector&,
144  AlgebraicVector4&,
145  TrajectoryStateOnSurface&,
146  TrajectoryStateOnSurface&);
147  void writeGlPosRcd(CLHEP::HepVector& d3);
148  inline double CLHEP_dot(const CLHEP::HepVector& a, const CLHEP::HepVector& b) {
149  return a(1) * b(1) + a(2) * b(2) + a(3) * b(3);
150  }
151 
152 private:
153  void beginJob() override;
154  void analyze(const edm::Event&, const edm::EventSetup&) override;
155  void endJob() override;
156 
157  // ----------member data ---------------------------
158  const edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> m_TkGeometryToken;
159  const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> m_MagFieldToken;
160  const edm::ESGetToken<Alignments, GlobalPositionRcd> m_globalPosToken;
161  const edm::ESGetToken<Propagator, TrackingComponentsRecord> m_propToken;
162  const edm::ESGetToken<TransientTrackingRecHitBuilder, TransientRecHitRecord> m_ttrhBuilderToken;
163 
164  const edm::InputTag trackTags_; // used to select what tracks to read from configuration file
165  const edm::InputTag muonTags_; // used to select what standalone muons
166  const edm::InputTag gmuonTags_; // used to select what global muons
167  const edm::InputTag smuonTags_; // used to select what selected muons
168  const std::string propagator_; // name of the propagator
169  bool cosmicMuonMode_;
170  bool isolatedMuonMode_;
171  bool collectionCosmic; // if true, read muon collection expected for CosmicMu
172  bool collectionIsolated; // IsolateMu
173  bool refitMuon_; // if true, refit stand alone muon track
174  bool refitTrack_; // tracker track
175  const std::string rootOutFile_;
176  const std::string txtOutFile_;
177  const bool writeDB_; // write results to DB
178  const bool debug_; // run in debugging mode
179 
180  const edm::EDGetTokenT<reco::TrackCollection> trackIsolateToken_;
181  const edm::EDGetTokenT<reco::TrackCollection> muonIsolateToken_;
182  const edm::EDGetTokenT<reco::TrackCollection> gmuonIsolateToken_;
183  const edm::EDGetTokenT<reco::MuonCollection> smuonIsolateToken_;
184  const edm::EDGetTokenT<reco::TrackCollection> trackCosmicToken_;
185  const edm::EDGetTokenT<reco::TrackCollection> muonCosmicToken_;
186  const edm::EDGetTokenT<reco::TrackCollection> gmuonCosmicToken_;
187  const edm::EDGetTokenT<reco::MuonCollection> smuonCosmicToken_;
188  const edm::EDGetTokenT<reco::TrackCollection> trackToken_;
189  const edm::EDGetTokenT<reco::TrackCollection> muonToken_;
190  const edm::EDGetTokenT<reco::TrackCollection> gmuonToken_;
191  const edm::EDGetTokenT<reco::MuonCollection> smuonToken_;
192 
193  edm::ESWatcher<GlobalTrackingGeometryRecord> watchTrackingGeometry_;
194  edm::ESHandle<GlobalTrackingGeometry> theTrackingGeometry;
195  const GlobalTrackingGeometry* trackingGeometry_;
196 
197  edm::ESWatcher<IdealMagneticFieldRecord> watchMagneticFieldRecord_;
198  const MagneticField* magneticField_;
199 
200  edm::ESWatcher<TrackingComponentsRecord> watchTrackingComponentsRecord_;
201 
202  edm::ESWatcher<GlobalPositionRcd> watchGlobalPositionRcd_;
203  const Alignments* globalPositionRcd_;
204 
205  KFTrajectoryFitter* theFitter;
206  KFTrajectorySmoother* theSmoother;
207  KFTrajectoryFitter* theFitterOp;
208  KFTrajectorySmoother* theSmootherOp;
209  bool defineFitter;
210  MuonTransientTrackingRecHitBuilder* MuRHBuilder;
211  const TransientTrackingRecHitBuilder* TTRHBuilder;
212 
213  // // LSF for global d(3): Inf * d = Gfr
214  AlgebraicVector3 Gfr; // free terms
215  AlgebraicSymMatrix33 Inf; // information matrix
216  // // LSF for global d(6): Hess * d = Grad
217  CLHEP::HepVector Grad; // gradient of the objective function in global parameters
218  CLHEP::HepMatrix Hess; // Hessian -- // ---
219 
220  CLHEP::HepVector GradL; // gradient of the objective function in local parameters
221  CLHEP::HepMatrix HessL; // Hessian -- // ---
222 
223  int N_event; // processed events
224  int N_track; // selected tracks
225 
226  std::vector<AlignTransform>::const_iterator
227  iteratorTrackerRcd; // location of Tracker in container globalPositionRcd_->m_aligm
228  std::vector<AlignTransform>::const_iterator iteratorMuonRcd; // Muon
229  std::vector<AlignTransform>::const_iterator iteratorEcalRcd; // Ecal
230  std::vector<AlignTransform>::const_iterator iteratorHcalRcd; // Hcal
231 
232  CLHEP::HepVector MuGlShift; // evaluated global muon shifts
233  CLHEP::HepVector MuGlAngle; // evaluated global muon angles
234 
235  std::string MuSelect; // what part of muon system is selected for 1st hit
236 
237  ofstream OutGlobalTxt; // output the vector of global alignment as text
238 
239  TFile* file;
240  TH1F* histo;
241  TH1F* histo2; // outerP
242  TH1F* histo3; // outerLambda = PI/2-outerTheta
243  TH1F* histo4; // phi
244  TH1F* histo5; // outerR
245  TH1F* histo6; // outerZ
246  TH1F* histo7; // s
247  TH1F* histo8; // chi^2 of muon-track
248 
249  TH1F* histo11; // |Rm-Rt|
250  TH1F* histo12; // Xm-Xt
251  TH1F* histo13; // Ym-Yt
252  TH1F* histo14; // Zm-Zt
253  TH1F* histo15; // Nxm-Nxt
254  TH1F* histo16; // Nym-Nyt
255  TH1F* histo17; // Nzm-Nzt
256  TH1F* histo18; // Error X of inner track
257  TH1F* histo19; // Error X of muon
258  TH1F* histo20; // Error of Xm-Xt
259  TH1F* histo21; // pull(Xm-Xt)
260  TH1F* histo22; // pull(Ym-Yt)
261  TH1F* histo23; // pull(Zm-Zt)
262  TH1F* histo24; // pull(PXm-PXt)
263  TH1F* histo25; // pull(PYm-Pyt)
264  TH1F* histo26; // pull(PZm-PZt)
265  TH1F* histo27; // Nx of tangent plane
266  TH1F* histo28; // Ny of tangent plane
267  TH1F* histo29; // lenght of inner track
268  TH1F* histo30; // lenght of muon track
269  TH1F* histo31; // chi2 local
270  TH1F* histo32; // pull(Pxm/Pzm - Pxt/Pzt) local
271  TH1F* histo33; // pull(Pym/Pzm - Pyt/Pzt) local
272  TH1F* histo34; // pull(Xm - Xt) local
273  TH1F* histo35; // pull(Ym - Yt) local
274 
275  TH2F* histo101; // Rtrack/muon vs Ztrack/muon
276  TH2F* histo102; // Ytrack/muon vs Xtrack/muon
277 };
278 
279 //
280 // constants, enums and typedefs
281 //
282 
283 //
284 // static data member definitions
285 //
286 
287 //
288 // constructors and destructor
289 //
290 GlobalTrackerMuonAlignment::GlobalTrackerMuonAlignment(const edm::ParameterSet& iConfig)
291  : m_TkGeometryToken(esConsumes()),
292  m_MagFieldToken(esConsumes()),
293  m_globalPosToken(esConsumes()),
294  m_propToken(esConsumes(edm::ESInputTag("", iConfig.getParameter<string>("Propagator")))),
295  m_ttrhBuilderToken(esConsumes(edm::ESInputTag("", "witTrackAngle"))),
296  trackTags_(iConfig.getParameter<edm::InputTag>("tracks")),
297  muonTags_(iConfig.getParameter<edm::InputTag>("muons")),
298  gmuonTags_(iConfig.getParameter<edm::InputTag>("gmuons")),
299  smuonTags_(iConfig.getParameter<edm::InputTag>("smuons")),
300  cosmicMuonMode_(iConfig.getParameter<bool>("cosmics")),
301  isolatedMuonMode_(iConfig.getParameter<bool>("isolated")),
302  refitMuon_(iConfig.getParameter<bool>("refitmuon")),
303  refitTrack_(iConfig.getParameter<bool>("refittrack")),
304  rootOutFile_(iConfig.getUntrackedParameter<string>("rootOutFile")),
305  txtOutFile_(iConfig.getUntrackedParameter<string>("txtOutFile")),
306  writeDB_(iConfig.getUntrackedParameter<bool>("writeDB")),
307  debug_(iConfig.getUntrackedParameter<bool>("debug")),
308  trackIsolateToken_(consumes<reco::TrackCollection>(edm::InputTag("ALCARECOMuAlCalIsolatedMu", "TrackerOnly"))),
309  muonIsolateToken_(consumes<reco::TrackCollection>(edm::InputTag("ALCARECOMuAlCalIsolatedMu", "StandAlone"))),
310  gmuonIsolateToken_(consumes<reco::TrackCollection>(edm::InputTag("ALCARECOMuAlCalIsolatedMu", "GlobalMuon"))),
311  smuonIsolateToken_(consumes<reco::MuonCollection>(edm::InputTag("ALCARECOMuAlCalIsolatedMu", "SelectedMuons"))),
312  trackCosmicToken_(
313  consumes<reco::TrackCollection>(edm::InputTag("ALCARECOMuAlCOMMuGlobalCosmicsMu", "TrackerOnly"))),
314  muonCosmicToken_(
315  consumes<reco::TrackCollection>(edm::InputTag("ALCARECOMuAlCOMMuGlobalCosmicsMu", "StandAlone"))),
316  gmuonCosmicToken_(
317  consumes<reco::TrackCollection>(edm::InputTag("ALCARECOMuAlCOMMuGlobalCosmicsMu", "GlobalMuon"))),
318  smuonCosmicToken_(
319  consumes<reco::MuonCollection>(edm::InputTag("ALCARECOMuAlCOMMuGlobalCosmicsMu", "SelectedMuons"))),
320  trackToken_(consumes<reco::TrackCollection>(trackTags_)),
321  muonToken_(consumes<reco::TrackCollection>(muonTags_)),
322  gmuonToken_(consumes<reco::TrackCollection>(gmuonTags_)),
323  smuonToken_(consumes<reco::MuonCollection>(smuonTags_)),
324  defineFitter(true) {
325 #ifdef NO_FALSE_FALSE_MODE
326  if (cosmicMuonMode_ == false && isolatedMuonMode_ == false) {
327  throw cms::Exception("BadConfig") << "Muon collection not set! "
328  << "Use from GlobalTrackerMuonAlignment_test_cfg.py.";
329  }
330 #endif
331  if (cosmicMuonMode_ == true && isolatedMuonMode_ == true) {
332  throw cms::Exception("BadConfig") << "Muon collection can be either cosmic or isolated! "
333  << "Please set only one to true.";
334  }
335 }
336 
337 void GlobalTrackerMuonAlignment::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
338  edm::ParameterSetDescription desc;
339  desc.add<std::string>("propagator", "SteppingHelixPropagator");
340  desc.add<edm::InputTag>("tracks", edm::InputTag("ALCARECOMuAlGlobalCosmics:TrackerOnly"));
341  desc.add<edm::InputTag>("muons", edm::InputTag("ALCARECOMuAlGlobalCosmics:StandAlone"));
342  desc.add<edm::InputTag>("gmuons", edm::InputTag("ALCARECOMuAlGlobalCosmics:GlobalMuon"));
343  desc.add<edm::InputTag>("smuons", edm::InputTag("ALCARECOMuAlGlobalCosmics:SelectedMuons"));
344  desc.add<bool>("isolated", false);
345  desc.add<bool>("cosmics", false);
346  desc.add<bool>("refitmuon", false);
347  desc.add<bool>("refittrack", false);
348  desc.addUntracked<std::string>("rootOutFile", "outfile.root");
349  desc.addUntracked<std::string>("txtOutFile", "outglobal.txt");
350  desc.addUntracked<bool>("writeDB", false);
351  desc.addUntracked<bool>("debug", false);
352  descriptions.add("globalTrackerMuonAlignment", desc);
353 }
354 //
355 // member functions
356 //
357 
358 // ------------ method called to for each event ------------
359 void GlobalTrackerMuonAlignment::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
360  //GlobalTrackerMuonAlignment::analyzeTrackTrack(iEvent, iSetup);
361  GlobalTrackerMuonAlignment::analyzeTrackTrajectory(iEvent, iSetup);
362 }
363 
364 // ------------ method called once for job just before starting event loop ------------
365 void GlobalTrackerMuonAlignment::beginJob() {
366  N_event = 0;
367  N_track = 0;
368 
369  if (cosmicMuonMode_ == true && isolatedMuonMode_ == false) {
370  collectionCosmic = true;
371  collectionIsolated = false;
372  } else if (cosmicMuonMode_ == false && isolatedMuonMode_ == true) {
373  collectionCosmic = false;
374  collectionIsolated = true;
375  } else {
376  collectionCosmic = false;
377  collectionIsolated = false;
378  }
379 
380  for (int i = 0; i <= 2; i++) {
381  Gfr(i) = 0.;
382  for (int j = 0; j <= 2; j++) {
383  Inf(i, j) = 0.;
384  }
385  }
386 
387  Grad = CLHEP::HepVector(6, 0);
388  Hess = CLHEP::HepSymMatrix(6, 0);
389 
390  GradL = CLHEP::HepVector(6, 0);
391  HessL = CLHEP::HepSymMatrix(6, 0);
392 
393  // histograms
394  TDirectory* dirsave = gDirectory;
395 
396  file = new TFile(rootOutFile_.c_str(), "recreate");
397  const bool oldAddDir = TH1::AddDirectoryStatus();
398 
399  TH1::AddDirectory(true);
400 
401  this->bookHist();
402 
403  TH1::AddDirectory(oldAddDir);
404  dirsave->cd();
405 }
406 
407 // ------------ method called once each job just after ending the event loop ------------
408 void GlobalTrackerMuonAlignment::endJob() {
409  bool alarm = false;
410 
411  this->fitHist();
412 
413  AlgebraicVector3 d(0., 0., 0.); // ------------ alignmnet Global Algebraic
414 
415  AlgebraicSymMatrix33 InfI; // inverse it
416  for (int i = 0; i <= 2; i++)
417  for (int j = 0; j <= 2; j++) {
418  if (j < i)
419  continue;
420  InfI(i, j) += Inf(i, j);
421  }
422  bool ierr = !InfI.Invert();
423  if (ierr) {
424  if (alarm)
425  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Error inverse Inf matrix !!!!!!!!!!!";
426  }
427 
428  for (int i = 0; i <= 2; i++)
429  for (int k = 0; k <= 2; k++)
430  d(i) -= InfI(i, k) * Gfr(k);
431  // end of Global Algebraic
432 
433  // --------------- alignment Global CLHEP
434  CLHEP::HepVector d3 = CLHEP::solve(Hess, -Grad);
435  int iEr3;
436  CLHEP::HepMatrix Errd3 = Hess.inverse(iEr3);
437  if (iEr3 != 0) {
438  if (alarm)
439  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " endJob Error inverse Hess matrix !!!!!!!!!!!";
440  }
441  // end of Global CLHEP
442 
443  // ----------------- alignment Local CLHEP
444  CLHEP::HepVector dLI = CLHEP::solve(HessL, -GradL);
445  int iErI;
446  CLHEP::HepMatrix ErrdLI = HessL.inverse(iErI);
447  if (iErI != 0) {
448  if (alarm)
449  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " endJob Error inverse HessL matrix !!!!!!!!!!!";
450  }
451  // end of Local CLHEP
452 
453  // printout of final parameters
454  edm::LogVerbatim("GlobalTrackerMuonAlignment")
455  << " ---- " << N_event << " event " << N_track << " tracks " << MuSelect << " ---- ";
456  if (collectionIsolated)
457  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ALCARECOMuAlCalIsolatedMu";
458  else if (collectionCosmic)
459  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ALCARECOMuAlGlobalCosmics";
460  else
461  edm::LogVerbatim("GlobalTrackerMuonAlignment") << smuonTags_;
462  edm::LogVerbatim("GlobalTrackerMuonAlignment")
463  << " Similated shifts[cm] " << MuGlShift(1) << " " << MuGlShift(2) << " " << MuGlShift(3) << " "
464  << " angles[rad] " << MuGlAngle(1) << " " << MuGlAngle(2) << " " << MuGlAngle(3) << " ";
465  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " d " << -d;
466 
467  edm::LogVerbatim("GlobalTrackerMuonAlignment")
468  << " +- " << std::sqrt(InfI(0, 0)) << " " << std::sqrt(InfI(1, 1)) << " " << std::sqrt(InfI(2, 2));
469  edm::LogVerbatim("GlobalTrackerMuonAlignment")
470  << " dG " << d3(1) << " " << d3(2) << " " << d3(3) << " " << d3(4) << " " << d3(5) << " " << d3(6);
471 
472  edm::LogVerbatim("GlobalTrackerMuonAlignment")
473  << " +- " << std::sqrt(Errd3(1, 1)) << " " << std::sqrt(Errd3(2, 2)) << " " << std::sqrt(Errd3(3, 3)) << " "
474  << std::sqrt(Errd3(4, 4)) << " " << std::sqrt(Errd3(5, 5)) << " " << std::sqrt(Errd3(6, 6));
475  edm::LogVerbatim("GlobalTrackerMuonAlignment")
476  << " dL " << dLI(1) << " " << dLI(2) << " " << dLI(3) << " " << dLI(4) << " " << dLI(5) << " " << dLI(6);
477 
478  edm::LogVerbatim("GlobalTrackerMuonAlignment")
479  << " +- " << std::sqrt(ErrdLI(1, 1)) << " " << std::sqrt(ErrdLI(2, 2)) << " " << std::sqrt(ErrdLI(3, 3)) << " "
480  << std::sqrt(ErrdLI(4, 4)) << " " << std::sqrt(ErrdLI(5, 5)) << " " << std::sqrt(ErrdLI(6, 6));
481 
482  // what do we write to DB
483  CLHEP::HepVector vectorToDb(6, 0), vectorErrToDb(6, 0);
484  //vectorToDb = d3;
485  //for(unsigned int i=1; i<=6; i++) vectorErrToDb(i) = std::sqrt(Errd3(i,i));
486  vectorToDb = -dLI;
487  for (unsigned int i = 1; i <= 6; i++)
488  vectorErrToDb(i) = std::sqrt(ErrdLI(i, i));
489 
490  // write histograms to root file
491  file->Write();
492  file->Close();
493 
494  // write global parameters to text file
495  OutGlobalTxt.open(txtOutFile_.c_str(), ios::out);
496  if (!OutGlobalTxt.is_open())
497  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " outglobal.txt is not open !!!!!";
498  else {
499  OutGlobalTxt << vectorToDb(1) << " " << vectorToDb(2) << " " << vectorToDb(3) << " " << vectorToDb(4) << " "
500  << vectorToDb(5) << " " << vectorToDb(6) << " muon Global.\n";
501  OutGlobalTxt << vectorErrToDb(1) << " " << vectorErrToDb(1) << " " << vectorErrToDb(1) << " " << vectorErrToDb(1)
502  << " " << vectorErrToDb(1) << " " << vectorErrToDb(1) << " errors.\n";
503  OutGlobalTxt << N_event << " events are processed.\n";
504 
505  if (collectionIsolated)
506  OutGlobalTxt << "ALCARECOMuAlCalIsolatedMu.\n";
507  else if (collectionCosmic)
508  OutGlobalTxt << " ALCARECOMuAlGlobalCosmics.\n";
509  else
510  OutGlobalTxt << smuonTags_ << ".\n";
511  OutGlobalTxt.close();
512  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Write to the file outglobal.txt done ";
513  }
514 
515  // write new GlobalPositionRcd to DB
516  if (debug_)
517  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " writeBD_ " << writeDB_;
518  if (writeDB_)
519  GlobalTrackerMuonAlignment::writeGlPosRcd(vectorToDb);
520 }
521 
522 // ------------ book histogram ------------
523 void GlobalTrackerMuonAlignment::bookHist() {
524  double PI = 3.1415927;
525  histo = new TH1F("Pt", "pt", 1000, 0, 100);
526  histo2 = new TH1F("P", "P [GeV/c]", 400, 0., 400.);
527  histo2->GetXaxis()->SetTitle("momentum [GeV/c]");
528  histo3 = new TH1F("outerLambda", "#lambda outer", 100, -PI / 2., PI / 2.);
529  histo3->GetXaxis()->SetTitle("#lambda outer");
530  histo4 = new TH1F("phi", "#phi [rad]", 100, -PI, PI);
531  histo4->GetXaxis()->SetTitle("#phi [rad]");
532  histo5 = new TH1F("Rmuon", "inner muon hit R [cm]", 100, 0., 800.);
533  histo5->GetXaxis()->SetTitle("R of muon [cm]");
534  histo6 = new TH1F("Zmuon", "inner muon hit Z[cm]", 100, -1000., 1000.);
535  histo6->GetXaxis()->SetTitle("Z of muon [cm]");
536  histo7 = new TH1F("(Pm-Pt)/Pt", " (Pmuon-Ptrack)/Ptrack", 100, -2., 2.);
537  histo7->GetXaxis()->SetTitle("(Pmuon-Ptrack)/Ptrack");
538  histo8 = new TH1F("chi muon-track", "#chi^{2}(muon-track)", 1000, 0., 1000.);
539  histo8->GetXaxis()->SetTitle("#chi^{2} of muon w.r.t. propagated track");
540  histo11 = new TH1F("distance muon-track", "distance muon w.r.t track [cm]", 100, 0., 30.);
541  histo11->GetXaxis()->SetTitle("distance of muon w.r.t. track [cm]");
542  histo12 = new TH1F("Xmuon-Xtrack", "Xmuon-Xtrack [cm]", 200, -20., 20.);
543  histo12->GetXaxis()->SetTitle("Xmuon - Xtrack [cm]");
544  histo13 = new TH1F("Ymuon-Ytrack", "Ymuon-Ytrack [cm]", 200, -20., 20.);
545  histo13->GetXaxis()->SetTitle("Ymuon - Ytrack [cm]");
546  histo14 = new TH1F("Zmuon-Ztrack", "Zmuon-Ztrack [cm]", 200, -20., 20.);
547  histo14->GetXaxis()->SetTitle("Zmuon-Ztrack [cm]");
548  histo15 = new TH1F("NXmuon-NXtrack", "NXmuon-NXtrack [rad]", 200, -.1, .1);
549  histo15->GetXaxis()->SetTitle("N_{X}(muon)-N_{X}(track) [rad]");
550  histo16 = new TH1F("NYmuon-NYtrack", "NYmuon-NYtrack [rad]", 200, -.1, .1);
551  histo16->GetXaxis()->SetTitle("N_{Y}(muon)-N_{Y}(track) [rad]");
552  histo17 = new TH1F("NZmuon-NZtrack", "NZmuon-NZtrack [rad]", 200, -.1, .1);
553  histo17->GetXaxis()->SetTitle("N_{Z}(muon)-N_{Z}(track) [rad]");
554  histo18 = new TH1F("expected error of Xinner", "outer hit of inner tracker", 100, 0, .01);
555  histo18->GetXaxis()->SetTitle("expected error of Xinner [cm]");
556  histo19 = new TH1F("expected error of Xmuon", "inner hit of muon", 100, 0, .1);
557  histo19->GetXaxis()->SetTitle("expected error of Xmuon [cm]");
558  histo20 = new TH1F("expected error of Xmuon-Xtrack", "muon w.r.t. propagated track", 100, 0., 10.);
559  histo20->GetXaxis()->SetTitle("expected error of Xmuon-Xtrack [cm]");
560  histo21 = new TH1F("pull of Xmuon-Xtrack", "pull of Xmuon-Xtrack", 100, -10., 10.);
561  histo21->GetXaxis()->SetTitle("(Xmuon-Xtrack)/expected error");
562  histo22 = new TH1F("pull of Ymuon-Ytrack", "pull of Ymuon-Ytrack", 100, -10., 10.);
563  histo22->GetXaxis()->SetTitle("(Ymuon-Ytrack)/expected error");
564  histo23 = new TH1F("pull of Zmuon-Ztrack", "pull of Zmuon-Ztrack", 100, -10., 10.);
565  histo23->GetXaxis()->SetTitle("(Zmuon-Ztrack)/expected error");
566  histo24 = new TH1F("pull of PXmuon-PXtrack", "pull of PXmuon-PXtrack", 100, -10., 10.);
567  histo24->GetXaxis()->SetTitle("(P_{X}(muon)-P_{X}(track))/expected error");
568  histo25 = new TH1F("pull of PYmuon-PYtrack", "pull of PYmuon-PYtrack", 100, -10., 10.);
569  histo25->GetXaxis()->SetTitle("(P_{Y}(muon)-P_{Y}(track))/expected error");
570  histo26 = new TH1F("pull of PZmuon-PZtrack", "pull of PZmuon-PZtrack", 100, -10., 10.);
571  histo26->GetXaxis()->SetTitle("(P_{Z}(muon)-P_{Z}(track))/expected error");
572  histo27 = new TH1F("N_x", "Nx of tangent plane", 120, -1.1, 1.1);
573  histo27->GetXaxis()->SetTitle("normal vector projection N_{X}");
574  histo28 = new TH1F("N_y", "Ny of tangent plane", 120, -1.1, 1.1);
575  histo28->GetXaxis()->SetTitle("normal vector projection N_{Y}");
576  histo29 = new TH1F("lenght of track", "lenght of track", 200, 0., 400);
577  histo29->GetXaxis()->SetTitle("lenght of track [cm]");
578  histo30 = new TH1F("lenght of muon", "lenght of muon", 200, 0., 800);
579  histo30->GetXaxis()->SetTitle("lenght of muon [cm]");
580 
581  histo31 = new TH1F("local chi muon-track", "#local chi^{2}(muon-track)", 1000, 0., 1000.);
582  histo31->GetXaxis()->SetTitle("#local chi^{2} of muon w.r.t. propagated track");
583  histo32 = new TH1F("pull of Px/Pz local", "pull of Px/Pz local", 100, -10., 10.);
584  histo32->GetXaxis()->SetTitle("local (Px/Pz(muon) - Px/Pz(track))/expected error");
585  histo33 = new TH1F("pull of Py/Pz local", "pull of Py/Pz local", 100, -10., 10.);
586  histo33->GetXaxis()->SetTitle("local (Py/Pz(muon) - Py/Pz(track))/expected error");
587  histo34 = new TH1F("pull of X local", "pull of X local", 100, -10., 10.);
588  histo34->GetXaxis()->SetTitle("local (Xmuon - Xtrack)/expected error");
589  histo35 = new TH1F("pull of Y local", "pull of Y local", 100, -10., 10.);
590  histo35->GetXaxis()->SetTitle("local (Ymuon - Ytrack)/expected error");
591 
592  histo101 = new TH2F("Rtr/mu vs Ztr/mu", "hit of track/muon", 100, -800., 800., 100, 0., 600.);
593  histo101->GetXaxis()->SetTitle("Z of track/muon [cm]");
594  histo101->GetYaxis()->SetTitle("R of track/muon [cm]");
595  histo102 = new TH2F("Ytr/mu vs Xtr/mu", "hit of track/muon", 100, -600., 600., 100, -600., 600.);
596  histo102->GetXaxis()->SetTitle("X of track/muon [cm]");
597  histo102->GetYaxis()->SetTitle("Y of track/muon [cm]");
598 }
599 
600 // ------------ fit histogram ------------
601 void GlobalTrackerMuonAlignment::fitHist() {
602  histo7->Fit("gaus", "Q");
603 
604  histo12->Fit("gaus", "Q");
605  histo13->Fit("gaus", "Q");
606  histo14->Fit("gaus", "Q");
607  histo15->Fit("gaus", "Q");
608  histo16->Fit("gaus", "Q");
609  histo17->Fit("gaus", "Q");
610 
611  histo21->Fit("gaus", "Q");
612  histo22->Fit("gaus", "Q");
613  histo23->Fit("gaus", "Q");
614  histo24->Fit("gaus", "Q");
615  histo25->Fit("gaus", "Q");
616  histo26->Fit("gaus", "Q");
617 
618  histo32->Fit("gaus", "Q");
619  histo33->Fit("gaus", "Q");
620  histo34->Fit("gaus", "Q");
621  histo35->Fit("gaus", "Q");
622 }
623 
624 // ------------ method to analyze recoTrack & recoMuon of Global Muon ------------
625 void GlobalTrackerMuonAlignment::analyzeTrackTrack(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
626  using namespace edm;
627  using namespace reco;
628  //debug_ = true;
629  bool info = false;
630  bool alarm = false;
631  //bool alarm = true;
632  double PI = 3.1415927;
633 
634  cosmicMuonMode_ = true; // true: both Cosmic and IsolatedMuon are treated with 1,2 tracks
635  isolatedMuonMode_ = !cosmicMuonMode_; //true: only IsolatedMuon are treated with 1 track
636 
637  //if(N_event == 1)
638  //GlobalPositionRcdRead1::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup);
639 
640  N_event++;
641  if (info || debug_) {
642  edm::LogVerbatim("GlobalTrackerMuonAlignment") << "----- event " << N_event << " -- tracks " << N_track << " ---";
643  if (cosmicMuonMode_)
644  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " treated as CosmicMu ";
645  if (isolatedMuonMode_)
646  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " treated as IsolatedMu ";
647  edm::LogVerbatim("GlobalTrackerMuonAlignment");
648  }
649 
650  const edm::Handle<reco::TrackCollection>& tracks =
651  ((collectionIsolated)
652  ? iEvent.getHandle(trackIsolateToken_)
653  : ((collectionCosmic) ? iEvent.getHandle(trackCosmicToken_) : iEvent.getHandle(trackToken_)));
654  const edm::Handle<reco::TrackCollection>& muons =
655  ((collectionIsolated) ? iEvent.getHandle(muonIsolateToken_)
656  : ((collectionCosmic) ? iEvent.getHandle(muonCosmicToken_) : iEvent.getHandle(muonToken_)));
657  const edm::Handle<reco::TrackCollection>& gmuons =
658  ((collectionIsolated)
659  ? iEvent.getHandle(gmuonIsolateToken_)
660  : ((collectionCosmic) ? iEvent.getHandle(gmuonCosmicToken_) : iEvent.getHandle(gmuonToken_)));
661  const edm::Handle<reco::MuonCollection>& smuons =
662  ((collectionIsolated)
663  ? iEvent.getHandle(smuonIsolateToken_)
664  : ((collectionCosmic) ? iEvent.getHandle(smuonCosmicToken_) : iEvent.getHandle(smuonToken_)));
665 
666  if (debug_) {
667  edm::LogVerbatim("GlobalTrackerMuonAlignment")
668  << " ievBunch " << iEvent.bunchCrossing() << " runN " << (int)iEvent.run();
669  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " N tracks s/amu gmu selmu " << tracks->size() << " "
670  << muons->size() << " " << gmuons->size() << " " << smuons->size();
671  for (MuonCollection::const_iterator itMuon = smuons->begin(); itMuon != smuons->end(); ++itMuon) {
672  edm::LogVerbatim("GlobalTrackerMuonAlignment")
673  << " is isolatValid Matches " << itMuon->isIsolationValid() << " " << itMuon->isMatchesValid();
674  }
675  }
676 
677  if (isolatedMuonMode_) { // ---- Only 1 Isolated Muon --------
678  if (tracks->size() != 1)
679  return;
680  if (muons->size() != 1)
681  return;
682  if (gmuons->size() != 1)
683  return;
684  if (smuons->size() != 1)
685  return;
686  }
687 
688  if (cosmicMuonMode_) { // ---- 1,2 Cosmic Muon --------
689  if (smuons->size() > 2)
690  return;
691  if (tracks->size() != smuons->size())
692  return;
693  if (muons->size() != smuons->size())
694  return;
695  if (gmuons->size() != smuons->size())
696  return;
697  }
698 
699  // ok mc_isolated_mu
700  //edm::ESHandle<TrackerGeometry> trackerGeometry;
701  //iSetup.get<TrackerDigiGeometryRecord>().get(trackerGeometry);
702  // ok mc_isolated_mu
703  //edm::ESHandle<DTGeometry> dtGeometry;
704  //iSetup.get<MuonGeometryRecord>().get(dtGeometry);
705  // don't work
706  //edm::ESHandle<CSCGeometry> cscGeometry;
707  //iSetup.get<MuonGeometryRecord>().get(cscGeometry);
708 
709  if (watchTrackingGeometry_.check(iSetup) || !trackingGeometry_) {
710  edm::ESHandle<GlobalTrackingGeometry> trackingGeometry = iSetup.getHandle(m_TkGeometryToken);
711  trackingGeometry_ = &*trackingGeometry;
712  }
713 
714  if (watchMagneticFieldRecord_.check(iSetup) || !magneticField_) {
715  edm::ESHandle<MagneticField> magneticField = iSetup.getHandle(m_MagFieldToken);
716  magneticField_ = &*magneticField;
717  }
718 
719  if (watchGlobalPositionRcd_.check(iSetup) || !globalPositionRcd_) {
720  edm::ESHandle<Alignments> globalPositionRcd = iSetup.getHandle(m_globalPosToken);
721  globalPositionRcd_ = &*globalPositionRcd;
722  for (std::vector<AlignTransform>::const_iterator i = globalPositionRcd_->m_align.begin();
723  i != globalPositionRcd_->m_align.end();
724  ++i) {
725  if (DetId(DetId::Tracker).rawId() == i->rawId())
726  iteratorTrackerRcd = i;
727  if (DetId(DetId::Muon).rawId() == i->rawId())
728  iteratorMuonRcd = i;
729  if (DetId(DetId::Ecal).rawId() == i->rawId())
730  iteratorEcalRcd = i;
731  if (DetId(DetId::Hcal).rawId() == i->rawId())
732  iteratorHcalRcd = i;
733  }
734  if (true || debug_) {
735  edm::LogVerbatim("GlobalTrackerMuonAlignment") << "=== iteratorMuonRcd " << iteratorMuonRcd->rawId() << " ====\n"
736  << " translation " << iteratorMuonRcd->translation() << "\n"
737  << " angles " << iteratorMuonRcd->rotation().eulerAngles();
738  edm::LogVerbatim("GlobalTrackerMuonAlignment") << iteratorMuonRcd->rotation();
739  }
740  } // end of GlobalPositionRcd
741 
742  edm::ESHandle<Propagator> propagator = iSetup.getHandle(m_propToken);
743 
744  SteppingHelixPropagator alongStHePr = SteppingHelixPropagator(magneticField_, alongMomentum);
745  SteppingHelixPropagator oppositeStHePr = SteppingHelixPropagator(magneticField_, oppositeToMomentum);
746 
747  //double tolerance = 5.e-5;
748  defaultRKPropagator::Product aprod(magneticField_, alongMomentum, 5.e-5);
749  auto& alongRKPr = aprod.propagator;
750  defaultRKPropagator::Product oprod(magneticField_, oppositeToMomentum, 5.e-5);
751  auto& oppositeRKPr = oprod.propagator;
752 
753  float epsilon = 5.;
754  SmartPropagator alongSmPr = SmartPropagator(alongRKPr, alongStHePr, magneticField_, alongMomentum, epsilon);
755  SmartPropagator oppositeSmPr =
756  SmartPropagator(oppositeRKPr, oppositeStHePr, magneticField_, oppositeToMomentum, epsilon);
757 
758  // ................................................ selected/global muon
759  //itMuon --> Jim's globalMuon
760  for (MuonCollection::const_iterator itMuon = smuons->begin(); itMuon != smuons->end(); ++itMuon) {
761  if (debug_) {
762  edm::LogVerbatim("GlobalTrackerMuonAlignment")
763  << " mu gM is GM Mu SaM tM " << itMuon->isGlobalMuon() << " " << itMuon->isMuon() << " "
764  << itMuon->isStandAloneMuon() << " " << itMuon->isTrackerMuon() << " ";
765  }
766 
767  // get information about the innermost muon hit -------------------------
768  TransientTrack muTT(itMuon->outerTrack(), magneticField_, trackingGeometry_);
769  TrajectoryStateOnSurface innerMuTSOS = muTT.innermostMeasurementState();
770  TrajectoryStateOnSurface outerMuTSOS = muTT.outermostMeasurementState();
771 
772  // get information about the outermost tracker hit -----------------------
773  TransientTrack trackTT(itMuon->track(), magneticField_, trackingGeometry_);
774  TrajectoryStateOnSurface outerTrackTSOS = trackTT.outermostMeasurementState();
775  TrajectoryStateOnSurface innerTrackTSOS = trackTT.innermostMeasurementState();
776 
777  GlobalPoint pointTrackIn = innerTrackTSOS.globalPosition();
778  GlobalPoint pointTrackOut = outerTrackTSOS.globalPosition();
779  float lenghtTrack = (pointTrackOut - pointTrackIn).mag();
780  GlobalPoint pointMuonIn = innerMuTSOS.globalPosition();
781  GlobalPoint pointMuonOut = outerMuTSOS.globalPosition();
782  float lenghtMuon = (pointMuonOut - pointMuonIn).mag();
783  GlobalVector momentumTrackOut = outerTrackTSOS.globalMomentum();
784  GlobalVector momentumTrackIn = innerTrackTSOS.globalMomentum();
785  float distanceInIn = (pointTrackIn - pointMuonIn).mag();
786  float distanceInOut = (pointTrackIn - pointMuonOut).mag();
787  float distanceOutIn = (pointTrackOut - pointMuonIn).mag();
788  float distanceOutOut = (pointTrackOut - pointMuonOut).mag();
789 
790  if (debug_) {
791  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " pointTrackIn " << pointTrackIn;
792  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Out " << pointTrackOut << " lenght " << lenghtTrack;
793  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " point MuonIn " << pointMuonIn;
794  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Out " << pointMuonOut << " lenght " << lenghtMuon;
795  edm::LogVerbatim("GlobalTrackerMuonAlignment")
796  << " momeTrackIn Out " << momentumTrackIn << " " << momentumTrackOut;
797  edm::LogVerbatim("GlobalTrackerMuonAlignment")
798  << " doi io ii oo " << distanceOutIn << " " << distanceInOut << " " << distanceInIn << " " << distanceOutOut;
799  }
800 
801  if (lenghtTrack < 90.)
802  continue;
803  if (lenghtMuon < 300.)
804  continue;
805  if (momentumTrackIn.mag() < 15. || momentumTrackOut.mag() < 15.)
806  continue;
807  if (trackTT.charge() != muTT.charge())
808  continue;
809 
810  if (debug_)
811  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " passed lenght momentum cuts";
812 
813  GlobalVector GRm, GPm, Nl, Rm, Pm, Rt, Pt, Rt0;
814  AlgebraicSymMatrix66 Cm, C0, Ce, C1;
815 
816  TrajectoryStateOnSurface extrapolationT;
817  int tsosMuonIf = 0;
818 
819  if (isolatedMuonMode_) { //------------------------------- Isolated Muon -----
820  const Surface& refSurface = innerMuTSOS.surface();
821  ConstReferenceCountingPointer<TangentPlane> tpMuLocal(refSurface.tangentPlane(innerMuTSOS.localPosition()));
822  Nl = tpMuLocal->normalVector();
823  ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
824  GlobalVector Ng = tpMuGlobal->normalVector();
825  Surface* surf = (Surface*)&refSurface;
826  const Plane* refPlane = dynamic_cast<Plane*>(surf);
827  GlobalVector Nlp = refPlane->normalVector();
828 
829  if (debug_) {
830  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
831  << " alfa " << int(asin(Nl.x()) * 57.29578);
832  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " global " << Ng.x() << " " << Ng.y() << " " << Ng.z();
833  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " lp " << Nlp.x() << " " << Nlp.y() << " " << Nlp.z();
834  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug")
835  << " Nlocal Nglobal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " " << Ng.x() << " " << Ng.y() << " "
836  << Ng.z() << " alfa " << static_cast<int>(asin(Nl.x()) * 57.29578);
837  }
838 
839  // extrapolation to innermost muon hit
840  extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
841  //extrapolationT = propagator->propagate(outerTrackTSOS, refSurface);
842 
843  if (!extrapolationT.isValid()) {
844  if (false & alarm)
845  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
846  //<<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
847  ;
848  continue;
849  }
850  tsosMuonIf = 1;
851  Rt = GlobalVector((extrapolationT.globalPosition()).x(),
852  (extrapolationT.globalPosition()).y(),
853  (extrapolationT.globalPosition()).z());
854 
855  Pt = extrapolationT.globalMomentum();
856  // global parameters of muon
857  GRm = GlobalVector(
858  (innerMuTSOS.globalPosition()).x(), (innerMuTSOS.globalPosition()).y(), (innerMuTSOS.globalPosition()).z());
859  GPm = innerMuTSOS.globalMomentum();
860  Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
861  (outerTrackTSOS.globalPosition()).y(),
862  (outerTrackTSOS.globalPosition()).z());
863  Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + innerMuTSOS.cartesianError().matrix());
864  C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
865  Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
866  C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
867 
868  } // ------------------------------- end Isolated Muon -----
869 
870  if (cosmicMuonMode_) { //------------------------------- Cosmic Muon -----
871 
872  if ((distanceOutIn <= distanceInOut) & (distanceOutIn <= distanceInIn) &
873  (distanceOutIn <= distanceOutOut)) { // ----- Out - In ------
874  if (debug_)
875  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ----- Out - In -----";
876 
877  const Surface& refSurface = innerMuTSOS.surface();
878  ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
879  Nl = tpMuGlobal->normalVector();
880 
881  // extrapolation to innermost muon hit
882  extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
883  //extrapolationT = propagator->propagate(outerTrackTSOS, refSurface);
884 
885  if (!extrapolationT.isValid()) {
886  if (false & alarm)
887  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
888  //<<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
889  ;
890  continue;
891  }
892  if (debug_)
893  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " extrapolationT.isValid " << extrapolationT.isValid();
894 
895  tsosMuonIf = 1;
896  Rt = GlobalVector((extrapolationT.globalPosition()).x(),
897  (extrapolationT.globalPosition()).y(),
898  (extrapolationT.globalPosition()).z());
899 
900  Pt = extrapolationT.globalMomentum();
901  // global parameters of muon
902  GRm = GlobalVector(
903  (innerMuTSOS.globalPosition()).x(), (innerMuTSOS.globalPosition()).y(), (innerMuTSOS.globalPosition()).z());
904  GPm = innerMuTSOS.globalMomentum();
905  Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
906  (outerTrackTSOS.globalPosition()).y(),
907  (outerTrackTSOS.globalPosition()).z());
908  Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + innerMuTSOS.cartesianError().matrix());
909  C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
910  Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
911  C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
912 
913  if (false & debug_) {
914  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ->propDir " << propagator->propagationDirection();
915  PropagationDirectionChooser Chooser;
916  edm::LogVerbatim("GlobalTrackerMuonAlignment")
917  << " propDirCh " << Chooser.operator()(*outerTrackTSOS.freeState(), refSurface) << " Ch == along "
918  << (alongMomentum == Chooser.operator()(*outerTrackTSOS.freeState(), refSurface));
919  edm::LogVerbatim("GlobalTrackerMuonAlignment")
920  << " --- Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
921  << " alfa " << int(asin(Nl.x()) * 57.29578);
922  }
923  } // enf of ---- Out - In -----
924 
925  else if ((distanceInOut <= distanceInIn) & (distanceInOut <= distanceOutIn) &
926  (distanceInOut <= distanceOutOut)) { // ----- In - Out ------
927  if (debug_)
928  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ----- In - Out -----";
929 
930  const Surface& refSurface = outerMuTSOS.surface();
931  ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
932  Nl = tpMuGlobal->normalVector();
933 
934  // extrapolation to outermost muon hit
935  extrapolationT = oppositeSmPr.propagate(innerTrackTSOS, refSurface);
936 
937  if (!extrapolationT.isValid()) {
938  if (false & alarm)
939  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
940  << "\a\a\a\a\a\a\a\a" << extrapolationT.isValid();
941  continue;
942  }
943  if (debug_)
944  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " extrapolationT.isValid " << extrapolationT.isValid();
945 
946  tsosMuonIf = 2;
947  Rt = GlobalVector((extrapolationT.globalPosition()).x(),
948  (extrapolationT.globalPosition()).y(),
949  (extrapolationT.globalPosition()).z());
950 
951  Pt = extrapolationT.globalMomentum();
952  // global parameters of muon
953  GRm = GlobalVector(
954  (outerMuTSOS.globalPosition()).x(), (outerMuTSOS.globalPosition()).y(), (outerMuTSOS.globalPosition()).z());
955  GPm = outerMuTSOS.globalMomentum();
956  Rt0 = GlobalVector((innerTrackTSOS.globalPosition()).x(),
957  (innerTrackTSOS.globalPosition()).y(),
958  (innerTrackTSOS.globalPosition()).z());
959  Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + outerMuTSOS.cartesianError().matrix());
960  C0 = AlgebraicSymMatrix66(innerTrackTSOS.cartesianError().matrix());
961  Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
962  C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
963 
964  if (false & debug_) {
965  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ->propDir " << propagator->propagationDirection();
966  PropagationDirectionChooser Chooser;
967  edm::LogVerbatim("GlobalTrackerMuonAlignment")
968  << " propDirCh " << Chooser.operator()(*innerTrackTSOS.freeState(), refSurface) << " Ch == oppisite "
969  << (oppositeToMomentum == Chooser.operator()(*innerTrackTSOS.freeState(), refSurface));
970  edm::LogVerbatim("GlobalTrackerMuonAlignment")
971  << " --- Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
972  << " alfa " << int(asin(Nl.x()) * 57.29578);
973  }
974  } // enf of ---- In - Out -----
975 
976  else if ((distanceOutOut <= distanceInOut) & (distanceOutOut <= distanceInIn) &
977  (distanceOutOut <= distanceOutIn)) { // ----- Out - Out ------
978  if (debug_)
979  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ----- Out - Out -----";
980 
981  // reject: momentum of track has opposite direction to muon track
982  continue;
983 
984  const Surface& refSurface = outerMuTSOS.surface();
985  ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
986  Nl = tpMuGlobal->normalVector();
987 
988  // extrapolation to outermost muon hit
989  extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
990 
991  if (!extrapolationT.isValid()) {
992  if (alarm)
993  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " !!!!!!!!! Catastrophe Out-Out extrapolationT.isValid "
994  << "\a\a\a\a\a\a\a\a" << extrapolationT.isValid();
995  continue;
996  }
997  if (debug_)
998  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " extrapolationT.isValid " << extrapolationT.isValid();
999 
1000  tsosMuonIf = 2;
1001  Rt = GlobalVector((extrapolationT.globalPosition()).x(),
1002  (extrapolationT.globalPosition()).y(),
1003  (extrapolationT.globalPosition()).z());
1004 
1005  Pt = extrapolationT.globalMomentum();
1006  // global parameters of muon
1007  GRm = GlobalVector(
1008  (outerMuTSOS.globalPosition()).x(), (outerMuTSOS.globalPosition()).y(), (outerMuTSOS.globalPosition()).z());
1009  GPm = outerMuTSOS.globalMomentum();
1010  Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
1011  (outerTrackTSOS.globalPosition()).y(),
1012  (outerTrackTSOS.globalPosition()).z());
1013  Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + outerMuTSOS.cartesianError().matrix());
1014  C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
1015  Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
1016  C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
1017 
1018  if (debug_) {
1019  PropagationDirectionChooser Chooser;
1020  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1021  << " propDirCh " << Chooser.operator()(*outerTrackTSOS.freeState(), refSurface) << " Ch == along "
1022  << (alongMomentum == Chooser.operator()(*outerTrackTSOS.freeState(), refSurface));
1023  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1024  << " --- Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1025  << " alfa " << int(asin(Nl.x()) * 57.29578);
1026  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1027  << " Nornal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1028  << " alfa " << int(asin(Nl.x()) * 57.29578);
1029  }
1030  } // enf of ---- Out - Out -----
1031  else {
1032  if (alarm)
1033  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1034  << " ------- !!!!!!!!!!!!!!! No proper Out - In \a\a\a\a\a\a\a";
1035  continue;
1036  }
1037 
1038  } // ------------------------------- end Cosmic Muon -----
1039 
1040  if (tsosMuonIf == 0) {
1041  if (info) {
1042  edm::LogVerbatim("GlobalTrackerMuonAlignment") << "No tsosMuon !!!!!!";
1043  continue;
1044  }
1045  }
1046  TrajectoryStateOnSurface tsosMuon;
1047  if (tsosMuonIf == 1)
1048  tsosMuon = muTT.innermostMeasurementState();
1049  else
1050  tsosMuon = muTT.outermostMeasurementState();
1051 
1052  //GlobalTrackerMuonAlignment::misalignMuon(GRm, GPm, Nl, Rt, Rm, Pm);
1053  AlgebraicVector4 LPRm; // muon local (dx/dz, dy/dz, x, y)
1054  GlobalTrackerMuonAlignment::misalignMuonL(GRm, GPm, Nl, Rt, Rm, Pm, LPRm, extrapolationT, tsosMuon);
1055 
1056  GlobalVector resR = Rm - Rt;
1057  GlobalVector resP0 = Pm - Pt;
1058  GlobalVector resP = Pm / Pm.mag() - Pt / Pt.mag();
1059  float RelMomResidual = (Pm.mag() - Pt.mag()) / (Pt.mag() + 1.e-6);
1060  ;
1061 
1062  AlgebraicVector6 Vm;
1063  Vm(0) = resR.x();
1064  Vm(1) = resR.y();
1065  Vm(2) = resR.z();
1066  Vm(3) = resP0.x();
1067  Vm(4) = resP0.y();
1068  Vm(5) = resP0.z();
1069  float Rmuon = Rm.perp();
1070  float Zmuon = Rm.z();
1071  float alfa_x = atan2(Nl.x(), Nl.y()) * 57.29578;
1072 
1073  if (debug_) {
1074  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1075  << " Nx Ny Nz alfa_x " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " " << alfa_x;
1076  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug")
1077  << " Rm " << Rm << "\n Rt " << Rt << "\n resR " << resR << "\n resP " << resP << " dp/p " << RelMomResidual;
1078  }
1079 
1080  double chi_d = 0;
1081  for (int i = 0; i <= 5; i++)
1082  chi_d += Vm(i) * Vm(i) / Cm(i, i);
1083 
1084  AlgebraicVector5 Vml(tsosMuon.localParameters().vector() - extrapolationT.localParameters().vector());
1085  AlgebraicSymMatrix55 m(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
1086  AlgebraicSymMatrix55 Cml(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
1087  bool ierrLoc = !m.Invert();
1088  if (ierrLoc && debug_ && info) {
1089  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ==== Error inverting Local covariance matrix ==== ";
1090  continue;
1091  }
1092  double chi_Loc = ROOT::Math::Similarity(Vml, m);
1093  if (debug_)
1094  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1095  << " chi_Loc px/pz/err " << chi_Loc << " " << Vml(1) / std::sqrt(Cml(1, 1));
1096 
1097  if (Pt.mag() < 15.)
1098  continue;
1099  if (Pm.mag() < 5.)
1100  continue;
1101 
1102  //if(Pt.mag() < 30.) continue; // momenum cut < 30GeV
1103  //if(Pt.mag() < 60.) continue; // momenum cut < 30GeV
1104  //if(Pt.mag() > 50.) continue; // momenum cut > 50GeV
1105  //if(Pt.mag() > 100.) continue; // momenum cut > 100GeV
1106  //if(trackTT.charge() < 0) continue; // select positive charge
1107  //if(trackTT.charge() > 0) continue; // select negative charge
1108 
1109  //if(fabs(resR.x()) > 5.) continue; // strong cut X
1110  //if(fabs(resR.y()) > 5.) continue; // Y
1111  //if(fabs(resR.z()) > 5.) continue; // Z
1112  //if(fabs(resR.mag()) > 7.5) continue; // dR
1113 
1114  /*
1115  //if(fabs(RelMomResidual) > 0.5) continue;
1116  if(fabs(resR.x()) > 20.) continue;
1117  if(fabs(resR.y()) > 20.) continue;
1118  if(fabs(resR.z()) > 20.) continue;
1119  if(fabs(resR.mag()) > 30.) continue;
1120  if(fabs(resP.x()) > 0.06) continue;
1121  if(fabs(resP.y()) > 0.06) continue;
1122  if(fabs(resP.z()) > 0.06) continue;
1123  if(chi_d > 40.) continue;
1124  */
1125 
1126  // select Barrel
1127  //if(Rmuon < 400. || Rmuon > 450.) continue;
1128  //if(Zmuon < -600. || Zmuon > 600.) continue;
1129  //if(fabs(Nl.z()) > 0.95) continue;
1130  //MuSelect = " Barrel";
1131  // EndCap1
1132  //if(Rmuon < 120. || Rmuon > 450.) continue;
1133  //if(Zmuon < -720.) continue;
1134  //if(Zmuon > -580.) continue;
1135  //if(fabs(Nl.z()) < 0.95) continue;
1136  //MuSelect = " EndCap1";
1137  // EndCap2
1138  //if(Rmuon < 120. || Rmuon > 450.) continue;
1139  //if(Zmuon > 720.) continue;
1140  //if(Zmuon < 580.) continue;
1141  //if(fabs(Nl.z()) < 0.95) continue;
1142  //MuSelect = " EndCap2";
1143  // select All
1144  if (Rmuon < 120. || Rmuon > 450.)
1145  continue;
1146  if (Zmuon < -720. || Zmuon > 720.)
1147  continue;
1148  MuSelect = " Barrel+EndCaps";
1149 
1150  if (debug_)
1151  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " .............. passed all cuts";
1152 
1153  N_track++;
1154  // gradient and Hessian for each track
1155 
1156  GlobalTrackerMuonAlignment::gradientGlobalAlg(Rt, Pt, Rm, Nl, Cm);
1157  GlobalTrackerMuonAlignment::gradientGlobal(Rt, Pt, Rm, Pm, Nl, Cm);
1158 
1159  CLHEP::HepSymMatrix covLoc(4, 0);
1160  for (int i = 1; i <= 4; i++)
1161  for (int j = 1; j <= i; j++) {
1162  covLoc(i, j) = (tsosMuon.localError().matrix() + extrapolationT.localError().matrix())(i, j);
1163  //if(i != j) Cov(i,j) = 0.;
1164  }
1165 
1166  const Surface& refSurface = tsosMuon.surface();
1167  CLHEP::HepMatrix rotLoc(3, 3, 0);
1168  rotLoc(1, 1) = refSurface.rotation().xx();
1169  rotLoc(1, 2) = refSurface.rotation().xy();
1170  rotLoc(1, 3) = refSurface.rotation().xz();
1171 
1172  rotLoc(2, 1) = refSurface.rotation().yx();
1173  rotLoc(2, 2) = refSurface.rotation().yy();
1174  rotLoc(2, 3) = refSurface.rotation().yz();
1175 
1176  rotLoc(3, 1) = refSurface.rotation().zx();
1177  rotLoc(3, 2) = refSurface.rotation().zy();
1178  rotLoc(3, 3) = refSurface.rotation().zz();
1179 
1180  CLHEP::HepVector posLoc(3);
1181  posLoc(1) = refSurface.position().x();
1182  posLoc(2) = refSurface.position().y();
1183  posLoc(3) = refSurface.position().z();
1184 
1185  GlobalTrackerMuonAlignment::gradientLocal(Rt, Pt, Rm, Pm, Nl, covLoc, rotLoc, posLoc, LPRm);
1186 
1187  if (debug_) {
1188  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Norm " << Nl;
1189  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " posLoc " << posLoc.T();
1190  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " rotLoc " << rotLoc;
1191  }
1192 
1193  // ----------------------------------------------------- fill histogram
1194  histo->Fill(itMuon->track()->pt());
1195 
1196  //histo2->Fill(itMuon->track()->outerP());
1197  histo2->Fill(Pt.mag());
1198  histo3->Fill((PI / 2. - itMuon->track()->outerTheta()));
1199  histo4->Fill(itMuon->track()->phi());
1200  histo5->Fill(Rmuon);
1201  histo6->Fill(Zmuon);
1202  histo7->Fill(RelMomResidual);
1203  //histo8->Fill(chi);
1204  histo8->Fill(chi_d);
1205 
1206  histo101->Fill(Zmuon, Rmuon);
1207  histo101->Fill(Rt0.z(), Rt0.perp());
1208  histo102->Fill(Rt0.x(), Rt0.y());
1209  histo102->Fill(Rm.x(), Rm.y());
1210 
1211  histo11->Fill(resR.mag());
1212  if (fabs(Nl.x()) < 0.98)
1213  histo12->Fill(resR.x());
1214  if (fabs(Nl.y()) < 0.98)
1215  histo13->Fill(resR.y());
1216  if (fabs(Nl.z()) < 0.98)
1217  histo14->Fill(resR.z());
1218  histo15->Fill(resP.x());
1219  histo16->Fill(resP.y());
1220  histo17->Fill(resP.z());
1221 
1222  if ((fabs(Nl.x()) < 0.98) && (fabs(Nl.y()) < 0.98) && (fabs(Nl.z()) < 0.98)) {
1223  histo18->Fill(std::sqrt(C0(0, 0)));
1224  histo19->Fill(std::sqrt(C1(0, 0)));
1225  histo20->Fill(std::sqrt(C1(0, 0) + Ce(0, 0)));
1226  }
1227  if (fabs(Nl.x()) < 0.98)
1228  histo21->Fill(Vm(0) / std::sqrt(Cm(0, 0)));
1229  if (fabs(Nl.y()) < 0.98)
1230  histo22->Fill(Vm(1) / std::sqrt(Cm(1, 1)));
1231  if (fabs(Nl.z()) < 0.98)
1232  histo23->Fill(Vm(2) / std::sqrt(Cm(2, 2)));
1233  histo24->Fill(Vm(3) / std::sqrt(C1(3, 3) + Ce(3, 3)));
1234  histo25->Fill(Vm(4) / std::sqrt(C1(4, 4) + Ce(4, 4)));
1235  histo26->Fill(Vm(5) / std::sqrt(C1(5, 5) + Ce(5, 5)));
1236  histo27->Fill(Nl.x());
1237  histo28->Fill(Nl.y());
1238  histo29->Fill(lenghtTrack);
1239  histo30->Fill(lenghtMuon);
1240  histo31->Fill(chi_Loc);
1241  histo32->Fill(Vml(1) / std::sqrt(Cml(1, 1)));
1242  histo33->Fill(Vml(2) / std::sqrt(Cml(2, 2)));
1243  histo34->Fill(Vml(3) / std::sqrt(Cml(3, 3)));
1244  histo35->Fill(Vml(4) / std::sqrt(Cml(4, 4)));
1245 
1246  if (debug_) { //--------------------------------- debug print ----------
1247 
1248  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " diag 0[ " << C0(0, 0) << " " << C0(1, 1) << " " << C0(2, 2)
1249  << " " << C0(3, 3) << " " << C0(4, 4) << " " << C0(5, 5) << " ]";
1250  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " diag e[ " << Ce(0, 0) << " " << Ce(1, 1) << " " << Ce(2, 2)
1251  << " " << Ce(3, 3) << " " << Ce(4, 4) << " " << Ce(5, 5) << " ]";
1252  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " diag 1[ " << C1(0, 0) << " " << C1(1, 1) << " " << C1(2, 2)
1253  << " " << C1(3, 3) << " " << C1(4, 4) << " " << C1(5, 5) << " ]";
1254  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Rm " << Rm.x() << " " << Rm.y() << " " << Rm.z() << " Pm "
1255  << Pm.x() << " " << Pm.y() << " " << Pm.z();
1256  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Rt " << Rt.x() << " " << Rt.y() << " " << Rt.z() << " Pt "
1257  << Pt.x() << " " << Pt.y() << " " << Pt.z();
1258  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Nl*(Rm-Rt) " << Nl.dot(Rm - Rt);
1259  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " resR " << resR.x() << " " << resR.y() << " " << resR.z()
1260  << " resP " << resP.x() << " " << resP.y() << " " << resP.z();
1261  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1262  << " Rm-t " << (Rm - Rt).x() << " " << (Rm - Rt).y() << " " << (Rm - Rt).z() << " Pm-t " << (Pm - Pt).x()
1263  << " " << (Pm - Pt).y() << " " << (Pm - Pt).z();
1264  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Vm " << Vm;
1265  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1266  << " +- " << std::sqrt(Cm(0, 0)) << " " << std::sqrt(Cm(1, 1)) << " " << std::sqrt(Cm(2, 2)) << " "
1267  << std::sqrt(Cm(3, 3)) << " " << std::sqrt(Cm(4, 4)) << " " << std::sqrt(Cm(5, 5));
1268  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1269  << " Pmuon Ptrack dP/Ptrack " << itMuon->outerTrack()->p() << " " << itMuon->track()->outerP() << " "
1270  << (itMuon->outerTrack()->p() - itMuon->track()->outerP()) / itMuon->track()->outerP();
1271  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " cov matrix ";
1272  edm::LogVerbatim("GlobalTrackerMuonAlignment") << Cm;
1273  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " diag [ " << Cm(0, 0) << " " << Cm(1, 1) << " " << Cm(2, 2)
1274  << " " << Cm(3, 3) << " " << Cm(4, 4) << " " << Cm(5, 5) << " ]";
1275 
1276  AlgebraicSymMatrix66 Ro;
1277  double Diag[6];
1278  for (int i = 0; i <= 5; i++)
1279  Diag[i] = std::sqrt(Cm(i, i));
1280  for (int i = 0; i <= 5; i++)
1281  for (int j = 0; j <= 5; j++)
1282  Ro(i, j) = Cm(i, j) / Diag[i] / Diag[j];
1283  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " correlation matrix ";
1284  edm::LogVerbatim("GlobalTrackerMuonAlignment") << Ro;
1285 
1286  AlgebraicSymMatrix66 CmI;
1287  for (int i = 0; i <= 5; i++)
1288  for (int j = 0; j <= 5; j++)
1289  CmI(i, j) = Cm(i, j);
1290 
1291  bool ierr = !CmI.Invert();
1292  if (ierr) {
1293  if (alarm || debug_)
1294  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Error inverse covariance matrix !!!!!!!!!!!";
1295  continue;
1296  }
1297  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " inverse cov matrix ";
1298  edm::LogVerbatim("GlobalTrackerMuonAlignment") << Cm;
1299 
1300  double chi = ROOT::Math::Similarity(Vm, CmI);
1301  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " chi chi_d " << chi << " " << chi_d;
1302  } // end of debug_ printout --------------------------------------------
1303 
1304  } // end loop on selected muons, i.e. Jim's globalMuon
1305 
1306 } //end of analyzeTrackTrack
1307 
1308 // ------- method to analyze recoTrack & trajectoryMuon of Global Muon ------
1309 void GlobalTrackerMuonAlignment::analyzeTrackTrajectory(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
1310  using namespace edm;
1311  using namespace reco;
1312  //debug_ = true;
1313  bool info = false;
1314  bool alarm = false;
1315  //bool alarm = true;
1316  double PI = 3.1415927;
1317 
1318  //-*-*-*-*-*-*-
1319  cosmicMuonMode_ = true; // true: both Cosmic and IsolatedMuon are treated with 1,2 tracks
1320  //cosmicMuonMode_ = false; // true: both Cosmic and IsolatedMuon are treated with 1,2 tracks
1321  //-*-*-*-*-*-*-
1322  isolatedMuonMode_ = !cosmicMuonMode_; //true: only IsolatedMuon are treated with 1 track
1323 
1324  N_event++;
1325  if (info || debug_) {
1326  edm::LogVerbatim("GlobalTrackerMuonAlignment") << "----- event " << N_event << " -- tracks " << N_track << " ---";
1327  if (cosmicMuonMode_)
1328  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " treated as CosmicMu ";
1329  if (isolatedMuonMode_)
1330  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " treated as IsolatedMu ";
1331  edm::LogVerbatim("GlobalTrackerMuonAlignment");
1332  }
1333 
1334  const edm::Handle<reco::TrackCollection>& tracks =
1335  ((collectionIsolated)
1336  ? iEvent.getHandle(trackIsolateToken_)
1337  : ((collectionCosmic) ? iEvent.getHandle(trackCosmicToken_) : iEvent.getHandle(trackToken_)));
1338  const edm::Handle<reco::TrackCollection>& muons =
1339  ((collectionIsolated) ? iEvent.getHandle(muonIsolateToken_)
1340  : ((collectionCosmic) ? iEvent.getHandle(muonCosmicToken_) : iEvent.getHandle(muonToken_)));
1341  const edm::Handle<reco::TrackCollection>& gmuons =
1342  ((collectionIsolated)
1343  ? iEvent.getHandle(gmuonIsolateToken_)
1344  : ((collectionCosmic) ? iEvent.getHandle(gmuonCosmicToken_) : iEvent.getHandle(gmuonToken_)));
1345  const edm::Handle<reco::MuonCollection>& smuons =
1346  ((collectionIsolated)
1347  ? iEvent.getHandle(smuonIsolateToken_)
1348  : ((collectionCosmic) ? iEvent.getHandle(smuonCosmicToken_) : iEvent.getHandle(smuonToken_)));
1349 
1350  if (debug_) {
1351  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1352  << " ievBunch " << iEvent.bunchCrossing() << " runN " << (int)iEvent.run();
1353  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " N tracks s/amu gmu selmu " << tracks->size() << " "
1354  << muons->size() << " " << gmuons->size() << " " << smuons->size();
1355  for (MuonCollection::const_iterator itMuon = smuons->begin(); itMuon != smuons->end(); ++itMuon) {
1356  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1357  << " is isolatValid Matches " << itMuon->isIsolationValid() << " " << itMuon->isMatchesValid();
1358  }
1359  }
1360 
1361  if (isolatedMuonMode_) { // ---- Only 1 Isolated Muon --------
1362  if (tracks->size() != 1)
1363  return;
1364  if (muons->size() != 1)
1365  return;
1366  if (gmuons->size() != 1)
1367  return;
1368  if (smuons->size() != 1)
1369  return;
1370  }
1371 
1372  if (cosmicMuonMode_) { // ---- 1,2 Cosmic Muon --------
1373  if (smuons->size() > 2)
1374  return;
1375  if (tracks->size() != smuons->size())
1376  return;
1377  if (muons->size() != smuons->size())
1378  return;
1379  if (gmuons->size() != smuons->size())
1380  return;
1381  }
1382 
1383  // ok mc_isolated_mu
1384  //edm::ESHandle<TrackerGeometry> trackerGeometry;
1385  //iSetup.get<TrackerDigiGeometryRecord>().get(trackerGeometry);
1386  // ok mc_isolated_mu
1387  //edm::ESHandle<DTGeometry> dtGeometry;
1388  //iSetup.get<MuonGeometryRecord>().get(dtGeometry);
1389  // don't work
1390  //edm::ESHandle<CSCGeometry> cscGeometry;
1391  //iSetup.get<MuonGeometryRecord>().get(cscGeometry);
1392 
1393  if (watchTrackingGeometry_.check(iSetup) || !trackingGeometry_) {
1394  edm::ESHandle<GlobalTrackingGeometry> trackingGeometry = iSetup.getHandle(m_TkGeometryToken);
1395  trackingGeometry_ = &*trackingGeometry;
1396  theTrackingGeometry = trackingGeometry;
1397  }
1398 
1399  if (watchMagneticFieldRecord_.check(iSetup) || !magneticField_) {
1400  edm::ESHandle<MagneticField> magneticField = iSetup.getHandle(m_MagFieldToken);
1401  magneticField_ = &*magneticField;
1402  }
1403 
1404  if (watchGlobalPositionRcd_.check(iSetup) || !globalPositionRcd_) {
1405  edm::ESHandle<Alignments> globalPositionRcd = iSetup.getHandle(m_globalPosToken);
1406  globalPositionRcd_ = &*globalPositionRcd;
1407  for (std::vector<AlignTransform>::const_iterator i = globalPositionRcd_->m_align.begin();
1408  i != globalPositionRcd_->m_align.end();
1409  ++i) {
1410  if (DetId(DetId::Tracker).rawId() == i->rawId())
1411  iteratorTrackerRcd = i;
1412  if (DetId(DetId::Muon).rawId() == i->rawId())
1413  iteratorMuonRcd = i;
1414  if (DetId(DetId::Ecal).rawId() == i->rawId())
1415  iteratorEcalRcd = i;
1416  if (DetId(DetId::Hcal).rawId() == i->rawId())
1417  iteratorHcalRcd = i;
1418  }
1419  if (debug_) {
1420  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1421  << "=== iteratorTrackerRcd " << iteratorTrackerRcd->rawId() << " ====\n"
1422  << " translation " << iteratorTrackerRcd->translation() << "\n"
1423  << " angles " << iteratorTrackerRcd->rotation().eulerAngles();
1424  edm::LogVerbatim("GlobalTrackerMuonAlignment") << iteratorTrackerRcd->rotation();
1425  edm::LogVerbatim("GlobalTrackerMuonAlignment") << "=== iteratorMuonRcd " << iteratorMuonRcd->rawId() << " ====\n"
1426  << " translation " << iteratorMuonRcd->translation() << "\n"
1427  << " angles " << iteratorMuonRcd->rotation().eulerAngles();
1428  edm::LogVerbatim("GlobalTrackerMuonAlignment") << iteratorMuonRcd->rotation();
1429  }
1430  } // end of GlobalPositionRcd
1431 
1432  edm::ESHandle<Propagator> propagator = iSetup.getHandle(m_propToken);
1433 
1434  SteppingHelixPropagator alongStHePr = SteppingHelixPropagator(magneticField_, alongMomentum);
1435  SteppingHelixPropagator oppositeStHePr = SteppingHelixPropagator(magneticField_, oppositeToMomentum);
1436 
1437  //double tolerance = 5.e-5;
1438  defaultRKPropagator::Product aprod(magneticField_, alongMomentum, 5.e-5);
1439  auto& alongRKPr = aprod.propagator;
1440  defaultRKPropagator::Product oprod(magneticField_, oppositeToMomentum, 5.e-5);
1441  auto& oppositeRKPr = oprod.propagator;
1442 
1443  float epsilon = 5.;
1444  SmartPropagator alongSmPr = SmartPropagator(alongRKPr, alongStHePr, magneticField_, alongMomentum, epsilon);
1445  SmartPropagator oppositeSmPr =
1446  SmartPropagator(oppositeRKPr, oppositeStHePr, magneticField_, oppositeToMomentum, epsilon);
1447 
1448  if (defineFitter) {
1449  if (debug_)
1450  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ............... DEFINE FITTER ...................";
1451  KFUpdator* theUpdator = new KFUpdator();
1452  //Chi2MeasurementEstimator* theEstimator = new Chi2MeasurementEstimator(30);
1453  Chi2MeasurementEstimator* theEstimator = new Chi2MeasurementEstimator(100000, 100000);
1454  theFitter = new KFTrajectoryFitter(alongSmPr, *theUpdator, *theEstimator);
1455  theSmoother = new KFTrajectorySmoother(alongSmPr, *theUpdator, *theEstimator);
1456  theFitterOp = new KFTrajectoryFitter(oppositeSmPr, *theUpdator, *theEstimator);
1457  theSmootherOp = new KFTrajectorySmoother(oppositeSmPr, *theUpdator, *theEstimator);
1458 
1459  edm::ESHandle<TransientTrackingRecHitBuilder> builder = iSetup.getHandle(m_ttrhBuilderToken);
1460  this->TTRHBuilder = &(*builder);
1461  this->MuRHBuilder = new MuonTransientTrackingRecHitBuilder(theTrackingGeometry);
1462  if (debug_) {
1463  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1464  << " theTrackingGeometry.isValid() " << theTrackingGeometry.isValid();
1465  edm::LogVerbatim("GlobalTrackerMuonAlignment") << "get also the MuonTransientTrackingRecHitBuilder"
1466  << "\n";
1467  edm::LogVerbatim("GlobalTrackerMuonAlignment") << "get also the TransientTrackingRecHitBuilder"
1468  << "\n";
1469  }
1470  defineFitter = false;
1471  }
1472 
1473  // ................................................ selected/global muon
1474  //itMuon --> Jim's globalMuon
1475  for (MuonCollection::const_iterator itMuon = smuons->begin(); itMuon != smuons->end(); ++itMuon) {
1476  if (debug_) {
1477  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1478  << " mu gM is GM Mu SaM tM " << itMuon->isGlobalMuon() << " " << itMuon->isMuon() << " "
1479  << itMuon->isStandAloneMuon() << " " << itMuon->isTrackerMuon() << " ";
1480  }
1481 
1482  // get information about the innermost muon hit -------------------------
1483  TransientTrack muTT(itMuon->outerTrack(), magneticField_, trackingGeometry_);
1484  TrajectoryStateOnSurface innerMuTSOS = muTT.innermostMeasurementState();
1485  TrajectoryStateOnSurface outerMuTSOS = muTT.outermostMeasurementState();
1486 
1487  // get information about the outermost tracker hit -----------------------
1488  TransientTrack trackTT(itMuon->track(), magneticField_, trackingGeometry_);
1489  TrajectoryStateOnSurface outerTrackTSOS = trackTT.outermostMeasurementState();
1490  TrajectoryStateOnSurface innerTrackTSOS = trackTT.innermostMeasurementState();
1491 
1492  GlobalPoint pointTrackIn = innerTrackTSOS.globalPosition();
1493  GlobalPoint pointTrackOut = outerTrackTSOS.globalPosition();
1494  float lenghtTrack = (pointTrackOut - pointTrackIn).mag();
1495  GlobalPoint pointMuonIn = innerMuTSOS.globalPosition();
1496  GlobalPoint pointMuonOut = outerMuTSOS.globalPosition();
1497  float lenghtMuon = (pointMuonOut - pointMuonIn).mag();
1498  GlobalVector momentumTrackOut = outerTrackTSOS.globalMomentum();
1499  GlobalVector momentumTrackIn = innerTrackTSOS.globalMomentum();
1500  float distanceInIn = (pointTrackIn - pointMuonIn).mag();
1501  float distanceInOut = (pointTrackIn - pointMuonOut).mag();
1502  float distanceOutIn = (pointTrackOut - pointMuonIn).mag();
1503  float distanceOutOut = (pointTrackOut - pointMuonOut).mag();
1504 
1505  if (debug_) {
1506  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " pointTrackIn " << pointTrackIn;
1507  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Out " << pointTrackOut << " lenght " << lenghtTrack;
1508  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " point MuonIn " << pointMuonIn;
1509  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Out " << pointMuonOut << " lenght " << lenghtMuon;
1510  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1511  << " momeTrackIn Out " << momentumTrackIn << " " << momentumTrackOut;
1512  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1513  << " doi io ii oo " << distanceOutIn << " " << distanceInOut << " " << distanceInIn << " " << distanceOutOut;
1514  }
1515 
1516  if (lenghtTrack < 90.)
1517  continue;
1518  if (lenghtMuon < 300.)
1519  continue;
1520  if (innerMuTSOS.globalMomentum().mag() < 5. || outerMuTSOS.globalMomentum().mag() < 5.)
1521  continue;
1522  if (momentumTrackIn.mag() < 15. || momentumTrackOut.mag() < 15.)
1523  continue;
1524  if (trackTT.charge() != muTT.charge())
1525  continue;
1526 
1527  if (debug_)
1528  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " passed lenght momentum cuts";
1529 
1530  GlobalVector GRm, GPm, Nl, Rm, Pm, Rt, Pt, Rt0;
1531  AlgebraicSymMatrix66 Cm, C0, Ce, C1;
1532 
1533  //extrapolationT = propagator->propagate(outerTrackTSOS, refSurface);
1534  TrajectoryStateOnSurface extrapolationT;
1535  int tsosMuonIf = 0;
1536 
1537  TrajectoryStateOnSurface muonFittedTSOS;
1538  TrajectoryStateOnSurface trackFittedTSOS;
1539 
1540  if (isolatedMuonMode_) { //------------------------------- Isolated Muon --- Out-In --
1541  if (debug_)
1542  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ------ Isolated (out-in) ----- ";
1543  const Surface& refSurface = innerMuTSOS.surface();
1544  ConstReferenceCountingPointer<TangentPlane> tpMuLocal(refSurface.tangentPlane(innerMuTSOS.localPosition()));
1545  Nl = tpMuLocal->normalVector();
1546  ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
1547  GlobalVector Ng = tpMuGlobal->normalVector();
1548  Surface* surf = (Surface*)&refSurface;
1549  const Plane* refPlane = dynamic_cast<Plane*>(surf);
1550  GlobalVector Nlp = refPlane->normalVector();
1551 
1552  if (debug_) {
1553  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1554  << " alfa " << int(asin(Nl.x()) * 57.29578);
1555  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " global " << Ng.x() << " " << Ng.y() << " " << Ng.z();
1556  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " lp " << Nlp.x() << " " << Nlp.y() << " " << Nlp.z();
1557  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug")
1558  << " Nlocal Nglobal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " " << Ng.x() << " " << Ng.y() << " "
1559  << Ng.z() << " alfa " << static_cast<int>(asin(Nl.x()) * 57.29578);
1560  }
1561 
1562  // extrapolation to innermost muon hit
1563 
1564  //extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
1565  if (!refitTrack_)
1566  extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
1567  else {
1568  GlobalTrackerMuonAlignment::trackFitter(itMuon->track(), trackTT, alongMomentum, trackFittedTSOS);
1569  if (trackFittedTSOS.isValid())
1570  extrapolationT = alongSmPr.propagate(trackFittedTSOS, refSurface);
1571  }
1572 
1573  if (!extrapolationT.isValid()) {
1574  if (false & alarm)
1575  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
1576  //<<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
1577  ;
1578  continue;
1579  }
1580  tsosMuonIf = 1;
1581  Rt = GlobalVector((extrapolationT.globalPosition()).x(),
1582  (extrapolationT.globalPosition()).y(),
1583  (extrapolationT.globalPosition()).z());
1584 
1585  Pt = extrapolationT.globalMomentum();
1586  // global parameters of muon
1587  GRm = GlobalVector(
1588  (innerMuTSOS.globalPosition()).x(), (innerMuTSOS.globalPosition()).y(), (innerMuTSOS.globalPosition()).z());
1589  GPm = innerMuTSOS.globalMomentum();
1590 
1591  Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
1592  (outerTrackTSOS.globalPosition()).y(),
1593  (outerTrackTSOS.globalPosition()).z());
1594  Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + innerMuTSOS.cartesianError().matrix());
1595  C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
1596  Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
1597  C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
1598 
1599  if (refitMuon_)
1600  GlobalTrackerMuonAlignment::muonFitter(itMuon->outerTrack(), muTT, oppositeToMomentum, muonFittedTSOS);
1601 
1602  } // ------------------------------- end Isolated Muon -- Out - In ---
1603 
1604  if (cosmicMuonMode_) { //------------------------------- Cosmic Muon -----
1605 
1606  if ((distanceOutIn <= distanceInOut) & (distanceOutIn <= distanceInIn) &
1607  (distanceOutIn <= distanceOutOut)) { // ----- Out - In ------
1608  if (debug_)
1609  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ----- Out - In -----";
1610 
1611  const Surface& refSurface = innerMuTSOS.surface();
1612  ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
1613  Nl = tpMuGlobal->normalVector();
1614 
1615  // extrapolation to innermost muon hit
1616  //extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
1617  if (!refitTrack_)
1618  extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
1619  else {
1620  GlobalTrackerMuonAlignment::trackFitter(itMuon->track(), trackTT, alongMomentum, trackFittedTSOS);
1621  if (trackFittedTSOS.isValid())
1622  extrapolationT = alongSmPr.propagate(trackFittedTSOS, refSurface);
1623  }
1624 
1625  if (!extrapolationT.isValid()) {
1626  if (false & alarm)
1627  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
1628  //<<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
1629  ;
1630  continue;
1631  }
1632  if (debug_)
1633  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " extrapolationT.isValid " << extrapolationT.isValid();
1634 
1635  tsosMuonIf = 1;
1636  Rt = GlobalVector((extrapolationT.globalPosition()).x(),
1637  (extrapolationT.globalPosition()).y(),
1638  (extrapolationT.globalPosition()).z());
1639 
1640  Pt = extrapolationT.globalMomentum();
1641  // global parameters of muon
1642  GRm = GlobalVector(
1643  (innerMuTSOS.globalPosition()).x(), (innerMuTSOS.globalPosition()).y(), (innerMuTSOS.globalPosition()).z());
1644  GPm = innerMuTSOS.globalMomentum();
1645  Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
1646  (outerTrackTSOS.globalPosition()).y(),
1647  (outerTrackTSOS.globalPosition()).z());
1648  Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + innerMuTSOS.cartesianError().matrix());
1649  C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
1650  Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
1651  C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
1652 
1653  if (refitMuon_)
1654  GlobalTrackerMuonAlignment::muonFitter(itMuon->outerTrack(), muTT, oppositeToMomentum, muonFittedTSOS);
1655 
1656  if (false & debug_) {
1657  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ->propDir " << propagator->propagationDirection();
1658  PropagationDirectionChooser Chooser;
1659  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1660  << " propDirCh " << Chooser.operator()(*outerTrackTSOS.freeState(), refSurface) << " Ch == along "
1661  << (alongMomentum == Chooser.operator()(*outerTrackTSOS.freeState(), refSurface));
1662  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1663  << " --- Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1664  << " alfa " << int(asin(Nl.x()) * 57.29578);
1665  }
1666  } // enf of ---- Out - In -----
1667 
1668  else if ((distanceInOut <= distanceInIn) & (distanceInOut <= distanceOutIn) &
1669  (distanceInOut <= distanceOutOut)) { // ----- In - Out ------
1670  if (debug_)
1671  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ----- In - Out -----";
1672 
1673  const Surface& refSurface = outerMuTSOS.surface();
1674  ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
1675  Nl = tpMuGlobal->normalVector();
1676 
1677  // extrapolation to outermost muon hit
1678  //extrapolationT = oppositeSmPr.propagate(innerTrackTSOS, refSurface);
1679  if (!refitTrack_)
1680  extrapolationT = oppositeSmPr.propagate(innerTrackTSOS, refSurface);
1681  else {
1682  GlobalTrackerMuonAlignment::trackFitter(itMuon->track(), trackTT, oppositeToMomentum, trackFittedTSOS);
1683  if (trackFittedTSOS.isValid())
1684  extrapolationT = oppositeSmPr.propagate(trackFittedTSOS, refSurface);
1685  }
1686 
1687  if (!extrapolationT.isValid()) {
1688  if (false & alarm)
1689  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
1690  << "\a\a\a\a\a\a\a\a" << extrapolationT.isValid();
1691  continue;
1692  }
1693  if (debug_)
1694  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " extrapolationT.isValid " << extrapolationT.isValid();
1695 
1696  tsosMuonIf = 2;
1697  Rt = GlobalVector((extrapolationT.globalPosition()).x(),
1698  (extrapolationT.globalPosition()).y(),
1699  (extrapolationT.globalPosition()).z());
1700 
1701  Pt = extrapolationT.globalMomentum();
1702  // global parameters of muon
1703  GRm = GlobalVector(
1704  (outerMuTSOS.globalPosition()).x(), (outerMuTSOS.globalPosition()).y(), (outerMuTSOS.globalPosition()).z());
1705  GPm = outerMuTSOS.globalMomentum();
1706  Rt0 = GlobalVector((innerTrackTSOS.globalPosition()).x(),
1707  (innerTrackTSOS.globalPosition()).y(),
1708  (innerTrackTSOS.globalPosition()).z());
1709  Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + outerMuTSOS.cartesianError().matrix());
1710  C0 = AlgebraicSymMatrix66(innerTrackTSOS.cartesianError().matrix());
1711  Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
1712  C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
1713 
1714  if (refitMuon_)
1715  GlobalTrackerMuonAlignment::muonFitter(itMuon->outerTrack(), muTT, alongMomentum, muonFittedTSOS);
1716 
1717  if (false & debug_) {
1718  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ->propDir " << propagator->propagationDirection();
1719  PropagationDirectionChooser Chooser;
1720  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1721  << " propDirCh " << Chooser.operator()(*innerTrackTSOS.freeState(), refSurface) << " Ch == oppisite "
1722  << (oppositeToMomentum == Chooser.operator()(*innerTrackTSOS.freeState(), refSurface));
1723  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1724  << " --- Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1725  << " alfa " << int(asin(Nl.x()) * 57.29578);
1726  }
1727  } // enf of ---- In - Out -----
1728 
1729  else if ((distanceOutOut <= distanceInOut) & (distanceOutOut <= distanceInIn) &
1730  (distanceOutOut <= distanceOutIn)) { // ----- Out - Out ------
1731  if (debug_)
1732  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ----- Out - Out -----";
1733 
1734  // reject: momentum of track has opposite direction to muon track
1735  continue;
1736 
1737  const Surface& refSurface = outerMuTSOS.surface();
1738  ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
1739  Nl = tpMuGlobal->normalVector();
1740 
1741  // extrapolation to outermost muon hit
1742  extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
1743 
1744  if (!extrapolationT.isValid()) {
1745  if (alarm)
1746  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " !!!!!!!!! Catastrophe Out-Out extrapolationT.isValid "
1747  << "\a\a\a\a\a\a\a\a" << extrapolationT.isValid();
1748  continue;
1749  }
1750  if (debug_)
1751  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " extrapolationT.isValid " << extrapolationT.isValid();
1752 
1753  tsosMuonIf = 2;
1754  Rt = GlobalVector((extrapolationT.globalPosition()).x(),
1755  (extrapolationT.globalPosition()).y(),
1756  (extrapolationT.globalPosition()).z());
1757 
1758  Pt = extrapolationT.globalMomentum();
1759  // global parameters of muon
1760  GRm = GlobalVector(
1761  (outerMuTSOS.globalPosition()).x(), (outerMuTSOS.globalPosition()).y(), (outerMuTSOS.globalPosition()).z());
1762  GPm = outerMuTSOS.globalMomentum();
1763  Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
1764  (outerTrackTSOS.globalPosition()).y(),
1765  (outerTrackTSOS.globalPosition()).z());
1766  Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + outerMuTSOS.cartesianError().matrix());
1767  C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
1768  Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
1769  C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
1770 
1771  if (debug_) {
1772  PropagationDirectionChooser Chooser;
1773  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1774  << " propDirCh " << Chooser.operator()(*outerTrackTSOS.freeState(), refSurface) << " Ch == along "
1775  << (alongMomentum == Chooser.operator()(*outerTrackTSOS.freeState(), refSurface));
1776  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1777  << " --- Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1778  << " alfa " << int(asin(Nl.x()) * 57.29578);
1779  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1780  << " Nornal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1781  << " alfa " << int(asin(Nl.x()) * 57.29578);
1782  }
1783  } // enf of ---- Out - Out -----
1784  else {
1785  if (alarm)
1786  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1787  << " ------- !!!!!!!!!!!!!!! No proper Out - In \a\a\a\a\a\a\a";
1788  continue;
1789  }
1790 
1791  } // ------------------------------- end Cosmic Muon -----
1792 
1793  if (tsosMuonIf == 0) {
1794  if (info) {
1795  edm::LogVerbatim("GlobalTrackerMuonAlignment") << "No tsosMuon !!!!!!";
1796  continue;
1797  }
1798  }
1799  TrajectoryStateOnSurface tsosMuon;
1800  if (tsosMuonIf == 1)
1801  tsosMuon = muTT.innermostMeasurementState();
1802  else
1803  tsosMuon = muTT.outermostMeasurementState();
1804 
1805  //GlobalTrackerMuonAlignment::misalignMuon(GRm, GPm, Nl, Rt, Rm, Pm);
1806  AlgebraicVector4 LPRm; // muon local (dx/dz, dy/dz, x, y)
1807  GlobalTrackerMuonAlignment::misalignMuonL(GRm, GPm, Nl, Rt, Rm, Pm, LPRm, extrapolationT, tsosMuon);
1808 
1809  if (refitTrack_) {
1810  if (!trackFittedTSOS.isValid()) {
1811  if (info)
1812  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ================= trackFittedTSOS notValid !!!!!!!! ";
1813  continue;
1814  }
1815  if (debug_)
1816  this->debugTrajectorySOS(" trackFittedTSOS ", trackFittedTSOS);
1817  }
1818 
1819  if (refitMuon_) {
1820  if (!muonFittedTSOS.isValid()) {
1821  if (info)
1822  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ================= muonFittedTSOS notValid !!!!!!!! ";
1823  continue;
1824  }
1825  if (debug_)
1826  this->debugTrajectorySOS(" muonFittedTSOS ", muonFittedTSOS);
1827  Rm = GlobalVector((muonFittedTSOS.globalPosition()).x(),
1828  (muonFittedTSOS.globalPosition()).y(),
1829  (muonFittedTSOS.globalPosition()).z());
1830  Pm = muonFittedTSOS.globalMomentum();
1831  LPRm = AlgebraicVector4(muonFittedTSOS.localParameters().vector()(1),
1832  muonFittedTSOS.localParameters().vector()(2),
1833  muonFittedTSOS.localParameters().vector()(3),
1834  muonFittedTSOS.localParameters().vector()(4));
1835  }
1836  GlobalVector resR = Rm - Rt;
1837  GlobalVector resP0 = Pm - Pt;
1838  GlobalVector resP = Pm / Pm.mag() - Pt / Pt.mag();
1839  float RelMomResidual = (Pm.mag() - Pt.mag()) / (Pt.mag() + 1.e-6);
1840  ;
1841 
1842  AlgebraicVector6 Vm;
1843  Vm(0) = resR.x();
1844  Vm(1) = resR.y();
1845  Vm(2) = resR.z();
1846  Vm(3) = resP0.x();
1847  Vm(4) = resP0.y();
1848  Vm(5) = resP0.z();
1849  float Rmuon = Rm.perp();
1850  float Zmuon = Rm.z();
1851  float alfa_x = atan2(Nl.x(), Nl.y()) * 57.29578;
1852 
1853  if (debug_) {
1854  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1855  << " Nx Ny Nz alfa_x " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " " << alfa_x;
1856  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug")
1857  << " Rm " << Rm << "\n Rt " << Rt << "\n resR " << resR << "\n resP " << resP << " dp/p " << RelMomResidual;
1858  }
1859 
1860  double chi_d = 0;
1861  for (int i = 0; i <= 5; i++)
1862  chi_d += Vm(i) * Vm(i) / Cm(i, i);
1863 
1864  AlgebraicVector5 Vml(tsosMuon.localParameters().vector() - extrapolationT.localParameters().vector());
1865  AlgebraicSymMatrix55 m(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
1866  AlgebraicSymMatrix55 Cml(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
1867  bool ierrLoc = !m.Invert();
1868  if (ierrLoc && debug_ && info) {
1869  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ==== Error inverting Local covariance matrix ==== ";
1870  continue;
1871  }
1872  double chi_Loc = ROOT::Math::Similarity(Vml, m);
1873  if (debug_)
1874  edm::LogVerbatim("GlobalTrackerMuonAlignment")
1875  << " chi_Loc px/pz/err " << chi_Loc << " " << Vml(1) / std::sqrt(Cml(1, 1));
1876 
1877  if (Pt.mag() < 15.)
1878  continue;
1879  if (Pm.mag() < 5.)
1880  continue;
1881 
1882  //if(Pt.mag() < 30.) continue; // momenum cut < 30GeV
1883  //if(Pt.mag() < 60.) continue; // momenum cut < 30GeV
1884  //if(Pt.mag() > 50.) continue; // momenum cut > 50GeV
1885  //if(Pt.mag() > 100.) continue; // momenum cut > 100GeV
1886  //if(trackTT.charge() < 0) continue; // select positive charge
1887  //if(trackTT.charge() > 0) continue; // select negative charge
1888 
1889  //if(fabs(resR.x()) > 5.) continue; // strong cut X
1890  //if(fabs(resR.y()) > 5.) continue; // Y
1891  //if(fabs(resR.z()) > 5.) continue; // Z
1892  //if(fabs(resR.mag()) > 7.5) continue; // dR
1893 
1894  //if(fabs(RelMomResidual) > 0.5) continue;
1895  if (fabs(resR.x()) > 20.)
1896  continue;
1897  if (fabs(resR.y()) > 20.)
1898  continue;
1899  if (fabs(resR.z()) > 20.)
1900  continue;
1901  if (fabs(resR.mag()) > 30.)
1902  continue;
1903  if (fabs(resP.x()) > 0.06)
1904  continue;
1905  if (fabs(resP.y()) > 0.06)
1906  continue;
1907  if (fabs(resP.z()) > 0.06)
1908  continue;
1909  if (chi_d > 40.)
1910  continue;
1911 
1912  // select Barrel
1913  //if(Rmuon < 400. || Rmuon > 450.) continue;
1914  //if(Zmuon < -600. || Zmuon > 600.) continue;
1915  //if(fabs(Nl.z()) > 0.95) continue;
1916  //MuSelect = " Barrel";
1917  // EndCap1
1918  //if(Rmuon < 120. || Rmuon > 450.) continue;
1919  //if(Zmuon < -720.) continue;
1920  //if(Zmuon > -580.) continue;
1921  //if(fabs(Nl.z()) < 0.95) continue;
1922  //MuSelect = " EndCap1";
1923  // EndCap2
1924  //if(Rmuon < 120. || Rmuon > 450.) continue;
1925  //if(Zmuon > 720.) continue;
1926  //if(Zmuon < 580.) continue;
1927  //if(fabs(Nl.z()) < 0.95) continue;
1928  //MuSelect = " EndCap2";
1929  // select All
1930  if (Rmuon < 120. || Rmuon > 450.)
1931  continue;
1932  if (Zmuon < -720. || Zmuon > 720.)
1933  continue;
1934  MuSelect = " Barrel+EndCaps";
1935 
1936  if (debug_)
1937  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " .............. passed all cuts";
1938 
1939  N_track++;
1940  // gradient and Hessian for each track
1941 
1942  GlobalTrackerMuonAlignment::gradientGlobalAlg(Rt, Pt, Rm, Nl, Cm);
1943  GlobalTrackerMuonAlignment::gradientGlobal(Rt, Pt, Rm, Pm, Nl, Cm);
1944 
1945  CLHEP::HepSymMatrix covLoc(4, 0);
1946  for (int i = 1; i <= 4; i++)
1947  for (int j = 1; j <= i; j++) {
1948  covLoc(i, j) = (tsosMuon.localError().matrix() + extrapolationT.localError().matrix())(i, j);
1949  //if(i != j) Cov(i,j) = 0.;
1950  }
1951 
1952  const Surface& refSurface = tsosMuon.surface();
1953  CLHEP::HepMatrix rotLoc(3, 3, 0);
1954  rotLoc(1, 1) = refSurface.rotation().xx();
1955  rotLoc(1, 2) = refSurface.rotation().xy();
1956  rotLoc(1, 3) = refSurface.rotation().xz();
1957 
1958  rotLoc(2, 1) = refSurface.rotation().yx();
1959  rotLoc(2, 2) = refSurface.rotation().yy();
1960  rotLoc(2, 3) = refSurface.rotation().yz();
1961 
1962  rotLoc(3, 1) = refSurface.rotation().zx();
1963  rotLoc(3, 2) = refSurface.rotation().zy();
1964  rotLoc(3, 3) = refSurface.rotation().zz();
1965 
1966  CLHEP::HepVector posLoc(3);
1967  posLoc(1) = refSurface.position().x();
1968  posLoc(2) = refSurface.position().y();
1969  posLoc(3) = refSurface.position().z();
1970 
1971  GlobalTrackerMuonAlignment::gradientLocal(Rt, Pt, Rm, Pm, Nl, covLoc, rotLoc, posLoc, LPRm);
1972 
1973  if (debug_) {
1974  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Norm " << Nl;
1975  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " posLoc " << posLoc.T();
1976  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " rotLoc " << rotLoc;
1977  }
1978 
1979  // ----------------------------------------------------- fill histogram
1980  histo->Fill(itMuon->track()->pt());
1981 
1982  //histo2->Fill(itMuon->track()->outerP());
1983  histo2->Fill(Pt.mag());
1984  histo3->Fill((PI / 2. - itMuon->track()->outerTheta()));
1985  histo4->Fill(itMuon->track()->phi());
1986  histo5->Fill(Rmuon);
1987  histo6->Fill(Zmuon);
1988  histo7->Fill(RelMomResidual);
1989  //histo8->Fill(chi);
1990  histo8->Fill(chi_d);
1991 
1992  histo101->Fill(Zmuon, Rmuon);
1993  histo101->Fill(Rt0.z(), Rt0.perp());
1994  histo102->Fill(Rt0.x(), Rt0.y());
1995  histo102->Fill(Rm.x(), Rm.y());
1996 
1997  histo11->Fill(resR.mag());
1998  if (fabs(Nl.x()) < 0.98)
1999  histo12->Fill(resR.x());
2000  if (fabs(Nl.y()) < 0.98)
2001  histo13->Fill(resR.y());
2002  if (fabs(Nl.z()) < 0.98)
2003  histo14->Fill(resR.z());
2004  histo15->Fill(resP.x());
2005  histo16->Fill(resP.y());
2006  histo17->Fill(resP.z());
2007 
2008  if ((fabs(Nl.x()) < 0.98) && (fabs(Nl.y()) < 0.98) && (fabs(Nl.z()) < 0.98)) {
2009  histo18->Fill(std::sqrt(C0(0, 0)));
2010  histo19->Fill(std::sqrt(C1(0, 0)));
2011  histo20->Fill(std::sqrt(C1(0, 0) + Ce(0, 0)));
2012  }
2013  if (fabs(Nl.x()) < 0.98)
2014  histo21->Fill(Vm(0) / std::sqrt(Cm(0, 0)));
2015  if (fabs(Nl.y()) < 0.98)
2016  histo22->Fill(Vm(1) / std::sqrt(Cm(1, 1)));
2017  if (fabs(Nl.z()) < 0.98)
2018  histo23->Fill(Vm(2) / std::sqrt(Cm(2, 2)));
2019  histo24->Fill(Vm(3) / std::sqrt(C1(3, 3) + Ce(3, 3)));
2020  histo25->Fill(Vm(4) / std::sqrt(C1(4, 4) + Ce(4, 4)));
2021  histo26->Fill(Vm(5) / std::sqrt(C1(5, 5) + Ce(5, 5)));
2022  histo27->Fill(Nl.x());
2023  histo28->Fill(Nl.y());
2024  histo29->Fill(lenghtTrack);
2025  histo30->Fill(lenghtMuon);
2026  histo31->Fill(chi_Loc);
2027  histo32->Fill(Vml(1) / std::sqrt(Cml(1, 1)));
2028  histo33->Fill(Vml(2) / std::sqrt(Cml(2, 2)));
2029  histo34->Fill(Vml(3) / std::sqrt(Cml(3, 3)));
2030  histo35->Fill(Vml(4) / std::sqrt(Cml(4, 4)));
2031 
2032  if (debug_) { //--------------------------------- debug print ----------
2033 
2034  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " diag 0[ " << C0(0, 0) << " " << C0(1, 1) << " " << C0(2, 2)
2035  << " " << C0(3, 3) << " " << C0(4, 4) << " " << C0(5, 5) << " ]";
2036  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " diag e[ " << Ce(0, 0) << " " << Ce(1, 1) << " " << Ce(2, 2)
2037  << " " << Ce(3, 3) << " " << Ce(4, 4) << " " << Ce(5, 5) << " ]";
2038  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " diag 1[ " << C1(0, 0) << " " << C1(1, 1) << " " << C1(2, 2)
2039  << " " << C1(3, 3) << " " << C1(4, 4) << " " << C1(5, 5) << " ]";
2040  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Rm " << Rm.x() << " " << Rm.y() << " " << Rm.z() << " Pm "
2041  << Pm.x() << " " << Pm.y() << " " << Pm.z();
2042  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Rt " << Rt.x() << " " << Rt.y() << " " << Rt.z() << " Pt "
2043  << Pt.x() << " " << Pt.y() << " " << Pt.z();
2044  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Nl*(Rm-Rt) " << Nl.dot(Rm - Rt);
2045  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " resR " << resR.x() << " " << resR.y() << " " << resR.z()
2046  << " resP " << resP.x() << " " << resP.y() << " " << resP.z();
2047  edm::LogVerbatim("GlobalTrackerMuonAlignment")
2048  << " Rm-t " << (Rm - Rt).x() << " " << (Rm - Rt).y() << " " << (Rm - Rt).z() << " Pm-t " << (Pm - Pt).x()
2049  << " " << (Pm - Pt).y() << " " << (Pm - Pt).z();
2050  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Vm " << Vm;
2051  edm::LogVerbatim("GlobalTrackerMuonAlignment")
2052  << " +- " << std::sqrt(Cm(0, 0)) << " " << std::sqrt(Cm(1, 1)) << " " << std::sqrt(Cm(2, 2)) << " "
2053  << std::sqrt(Cm(3, 3)) << " " << std::sqrt(Cm(4, 4)) << " " << std::sqrt(Cm(5, 5));
2054  edm::LogVerbatim("GlobalTrackerMuonAlignment")
2055  << " Pmuon Ptrack dP/Ptrack " << itMuon->outerTrack()->p() << " " << itMuon->track()->outerP() << " "
2056  << (itMuon->outerTrack()->p() - itMuon->track()->outerP()) / itMuon->track()->outerP();
2057  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " cov matrix ";
2058  edm::LogVerbatim("GlobalTrackerMuonAlignment") << Cm;
2059  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " diag [ " << Cm(0, 0) << " " << Cm(1, 1) << " " << Cm(2, 2)
2060  << " " << Cm(3, 3) << " " << Cm(4, 4) << " " << Cm(5, 5) << " ]";
2061 
2062  AlgebraicSymMatrix66 Ro;
2063  double Diag[6];
2064  for (int i = 0; i <= 5; i++)
2065  Diag[i] = std::sqrt(Cm(i, i));
2066  for (int i = 0; i <= 5; i++)
2067  for (int j = 0; j <= 5; j++)
2068  Ro(i, j) = Cm(i, j) / Diag[i] / Diag[j];
2069  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " correlation matrix ";
2070  edm::LogVerbatim("GlobalTrackerMuonAlignment") << Ro;
2071 
2072  AlgebraicSymMatrix66 CmI;
2073  for (int i = 0; i <= 5; i++)
2074  for (int j = 0; j <= 5; j++)
2075  CmI(i, j) = Cm(i, j);
2076 
2077  bool ierr = !CmI.Invert();
2078  if (ierr) {
2079  if (alarm || debug_)
2080  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Error inverse covariance matrix !!!!!!!!!!!";
2081  continue;
2082  }
2083  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " inverse cov matrix ";
2084  edm::LogVerbatim("GlobalTrackerMuonAlignment") << Cm;
2085 
2086  double chi = ROOT::Math::Similarity(Vm, CmI);
2087  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " chi chi_d " << chi << " " << chi_d;
2088  } // end of debug_ printout --------------------------------------------
2089 
2090  } // end loop on selected muons, i.e. Jim's globalMuon
2091 
2092 } //end of analyzeTrackTrajectory
2093 
2094 // ---- calculate gradient and Hessian matrix (algebraic) to search global shifts ------
2095 void GlobalTrackerMuonAlignment::gradientGlobalAlg(
2096  GlobalVector& Rt, GlobalVector& Pt, GlobalVector& Rm, GlobalVector& Nl, AlgebraicSymMatrix66& Cm) {
2097  // ---------------------------- Calculate Information matrix and Gfree vector
2098  // Information == Hessian , Gfree == gradient of the objective function
2099 
2100  AlgebraicMatrix33 Jac;
2101  AlgebraicVector3 Wi, R_m, R_t, P_t, Norm, dR;
2102 
2103  R_m(0) = Rm.x();
2104  R_m(1) = Rm.y();
2105  R_m(2) = Rm.z();
2106  R_t(0) = Rt.x();
2107  R_t(1) = Rt.y();
2108  R_t(2) = Rt.z();
2109  P_t(0) = Pt.x();
2110  P_t(1) = Pt.y();
2111  P_t(2) = Pt.z();
2112  Norm(0) = Nl.x();
2113  Norm(1) = Nl.y();
2114  Norm(2) = Nl.z();
2115 
2116  for (int i = 0; i <= 2; i++) {
2117  if (Cm(i, i) > 1.e-20)
2118  Wi(i) = 1. / Cm(i, i);
2119  else
2120  Wi(i) = 1.e-10;
2121  dR(i) = R_m(i) - R_t(i);
2122  }
2123 
2124  float PtN = P_t(0) * Norm(0) + P_t(1) * Norm(1) + P_t(2) * Norm(2);
2125 
2126  Jac(0, 0) = 1. - P_t(0) * Norm(0) / PtN;
2127  Jac(0, 1) = -P_t(0) * Norm(1) / PtN;
2128  Jac(0, 2) = -P_t(0) * Norm(2) / PtN;
2129 
2130  Jac(1, 0) = -P_t(1) * Norm(0) / PtN;
2131  Jac(1, 1) = 1. - P_t(1) * Norm(1) / PtN;
2132  Jac(1, 2) = -P_t(1) * Norm(2) / PtN;
2133 
2134  Jac(2, 0) = -P_t(2) * Norm(0) / PtN;
2135  Jac(2, 1) = -P_t(2) * Norm(1) / PtN;
2136  Jac(2, 2) = 1. - P_t(2) * Norm(2) / PtN;
2137 
2138  AlgebraicSymMatrix33 Itr;
2139 
2140  for (int i = 0; i <= 2; i++)
2141  for (int j = 0; j <= 2; j++) {
2142  if (j < i)
2143  continue;
2144  Itr(i, j) = 0.;
2145  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ij " << i << " " << j;
2146  for (int k = 0; k <= 2; k++) {
2147  Itr(i, j) += Jac(k, i) * Wi(k) * Jac(k, j);
2148  }
2149  }
2150 
2151  for (int i = 0; i <= 2; i++)
2152  for (int j = 0; j <= 2; j++) {
2153  if (j < i)
2154  continue;
2155  Inf(i, j) += Itr(i, j);
2156  }
2157 
2158  AlgebraicVector3 Gtr(0., 0., 0.);
2159  for (int i = 0; i <= 2; i++)
2160  for (int k = 0; k <= 2; k++)
2161  Gtr(i) += dR(k) * Wi(k) * Jac(k, i);
2162  for (int i = 0; i <= 2; i++)
2163  Gfr(i) += Gtr(i);
2164 
2165  if (debug_) {
2166  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Wi " << Wi;
2167  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " N " << Norm;
2168  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " P_t " << P_t;
2169  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " (Pt*N) " << PtN;
2170  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " dR " << dR;
2171  edm::LogVerbatim("GlobalTrackerMuonAlignment")
2172  << " +/- " << 1. / std::sqrt(Wi(0)) << " " << 1. / std::sqrt(Wi(1)) << " " << 1. / std::sqrt(Wi(2)) << " ";
2173  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Jacobian dr/ddx ";
2174  edm::LogVerbatim("GlobalTrackerMuonAlignment") << Jac;
2175  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " G-- " << Gtr;
2176  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Itrack ";
2177  edm::LogVerbatim("GlobalTrackerMuonAlignment") << Itr;
2178  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Gfr " << Gfr;
2179  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -- Inf --";
2180  edm::LogVerbatim("GlobalTrackerMuonAlignment") << Inf;
2181  }
2182 
2183  return;
2184 }
2185 
2186 // ---- calculate gradient and Hessian matrix in global parameters ------
2187 void GlobalTrackerMuonAlignment::gradientGlobal(GlobalVector& GRt,
2188  GlobalVector& GPt,
2189  GlobalVector& GRm,
2190  GlobalVector& GPm,
2191  GlobalVector& GNorm,
2192  AlgebraicSymMatrix66& GCov) {
2193  // we search for 6D global correction vector (d, a), where
2194  // 3D vector of shihts d
2195  // 3D vector of rotation angles a
2196 
2197  //bool alarm = true;
2198  bool alarm = false;
2199  bool info = false;
2200 
2201  int Nd = 6; // dimension of vector of alignment pararmeters, d
2202 
2203  //double PtMom = GPt.mag();
2204  CLHEP::HepSymMatrix w(Nd, 0);
2205  for (int i = 1; i <= Nd; i++)
2206  for (int j = 1; j <= Nd; j++) {
2207  if (j <= i)
2208  w(i, j) = GCov(i - 1, j - 1);
2209  //if(i >= 3) w(i,j) /= PtMom;
2210  //if(j >= 3) w(i,j) /= PtMom;
2211  if ((i == j) && (i <= 3) && (GCov(i - 1, j - 1) < 1.e-20))
2212  w(i, j) = 1.e20; // w=0
2213  if (i != j)
2214  w(i, j) = 0.; // use diaginal elements
2215  }
2216 
2217  //GPt /= GPt.mag();
2218  //GPm /= GPm.mag(); // end of transform
2219 
2220  CLHEP::HepVector V(Nd), Rt(3), Pt(3), Rm(3), Pm(3), Norm(3);
2221  Rt(1) = GRt.x();
2222  Rt(2) = GRt.y();
2223  Rt(3) = GRt.z();
2224  Pt(1) = GPt.x();
2225  Pt(2) = GPt.y();
2226  Pt(3) = GPt.z();
2227  Rm(1) = GRm.x();
2228  Rm(2) = GRm.y();
2229  Rm(3) = GRm.z();
2230  Pm(1) = GPm.x();
2231  Pm(2) = GPm.y();
2232  Pm(3) = GPm.z();
2233  Norm(1) = GNorm.x();
2234  Norm(2) = GNorm.y();
2235  Norm(3) = GNorm.z();
2236 
2237  V = dsum(Rm - Rt, Pm - Pt);
2238  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " V " << V.T();
2239 
2240  //double PmN = CLHEP::dot(Pm, Norm);
2241  double PmN = CLHEP_dot(Pm, Norm);
2242 
2243  CLHEP::HepMatrix Jac(Nd, Nd, 0);
2244  for (int i = 1; i <= 3; i++)
2245  for (int j = 1; j <= 3; j++) {
2246  Jac(i, j) = Pm(i) * Norm(j) / PmN;
2247  if (i == j)
2248  Jac(i, j) -= 1.;
2249  }
2250 
2251  // dp/da
2252  Jac(4, 4) = 0.; // dpx/dax
2253  Jac(5, 4) = -Pm(3); // dpy/dax
2254  Jac(6, 4) = Pm(2); // dpz/dax
2255  Jac(4, 5) = Pm(3); // dpx/day
2256  Jac(5, 5) = 0.; // dpy/day
2257  Jac(6, 5) = -Pm(1); // dpz/day
2258  Jac(4, 6) = -Pm(2); // dpx/daz
2259  Jac(5, 6) = Pm(1); // dpy/daz
2260  Jac(6, 6) = 0.; // dpz/daz
2261 
2262  CLHEP::HepVector dsda(3);
2263  dsda(1) = (Norm(2) * Rm(3) - Norm(3) * Rm(2)) / PmN;
2264  dsda(2) = (Norm(3) * Rm(1) - Norm(1) * Rm(3)) / PmN;
2265  dsda(3) = (Norm(1) * Rm(2) - Norm(2) * Rm(1)) / PmN;
2266 
2267  // dr/da
2268  Jac(1, 4) = Pm(1) * dsda(1); // drx/dax
2269  Jac(2, 4) = -Rm(3) + Pm(2) * dsda(1); // dry/dax
2270  Jac(3, 4) = Rm(2) + Pm(3) * dsda(1); // drz/dax
2271 
2272  Jac(1, 5) = Rm(3) + Pm(1) * dsda(2); // drx/day
2273  Jac(2, 5) = Pm(2) * dsda(2); // dry/day
2274  Jac(3, 5) = -Rm(1) + Pm(3) * dsda(2); // drz/day
2275 
2276  Jac(1, 6) = -Rm(2) + Pm(1) * dsda(3); // drx/daz
2277  Jac(2, 6) = Rm(1) + Pm(2) * dsda(3); // dry/daz
2278  Jac(3, 6) = Pm(3) * dsda(3); // drz/daz
2279 
2280  CLHEP::HepSymMatrix W(Nd, 0);
2281  int ierr;
2282  W = w.inverse(ierr);
2283  if (ierr != 0) {
2284  if (alarm)
2285  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " gradientGlobal: inversion of matrix w fail ";
2286  return;
2287  }
2288 
2289  CLHEP::HepMatrix W_Jac(Nd, Nd, 0);
2290  W_Jac = Jac.T() * W;
2291 
2292  CLHEP::HepVector grad3(Nd);
2293  grad3 = W_Jac * V;
2294 
2295  CLHEP::HepMatrix hess3(Nd, Nd);
2296  hess3 = Jac.T() * W * Jac;
2297  //hess3(4,4) = 1.e-10; hess3(5,5) = 1.e-10; hess3(6,6) = 1.e-10; //????????????????
2298 
2299  Grad += grad3;
2300  Hess += hess3;
2301 
2302  CLHEP::HepVector d3I = CLHEP::solve(Hess, -Grad);
2303  int iEr3I;
2304  CLHEP::HepMatrix Errd3I = Hess.inverse(iEr3I);
2305  if (iEr3I != 0) {
2306  if (alarm)
2307  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " gradientGlobal error inverse Hess matrix !!!!!!!!!!!";
2308  }
2309 
2310  if (info || debug_) {
2311  edm::LogVerbatim("GlobalTrackerMuonAlignment")
2312  << " dG " << d3I(1) << " " << d3I(2) << " " << d3I(3) << " " << d3I(4) << " " << d3I(5) << " " << d3I(6);
2313  ;
2314  edm::LogVerbatim("GlobalTrackerMuonAlignment")
2315  << " +- " << std::sqrt(Errd3I(1, 1)) << " " << std::sqrt(Errd3I(2, 2)) << " " << std::sqrt(Errd3I(3, 3))
2316  << " " << std::sqrt(Errd3I(4, 4)) << " " << std::sqrt(Errd3I(5, 5)) << " " << std::sqrt(Errd3I(6, 6));
2317  }
2318 
2319 #ifdef CHECK_OF_DERIVATIVES
2320  // -------------------- check of derivatives
2321 
2322  CLHEP::HepVector d(3, 0), a(3, 0);
2323  CLHEP::HepMatrix T(3, 3, 1);
2324 
2325  CLHEP::HepMatrix Ti = T.T();
2326  //double A = CLHEP::dot(Ti*Pm, Norm);
2327  //double B = CLHEP::dot((Rt -Ti*Rm + Ti*d), Norm);
2328  double A = CLHEP_dot(Ti * Pm, Norm);
2329  double B = CLHEP_dot((Rt - Ti * Rm + Ti * d), Norm);
2330  double s0 = B / A;
2331 
2332  CLHEP::HepVector r0(3, 0), p0(3, 0);
2333  r0 = Ti * Rm - Ti * d + s0 * (Ti * Pm) - Rt;
2334  p0 = Ti * Pm - Pt;
2335 
2336  double delta = 0.0001;
2337 
2338  int ii = 3;
2339  d(ii) += delta; // d
2340  //T(2,3) += delta; T(3,2) -= delta; int ii = 1; // a1
2341  //T(3,1) += delta; T(1,3) -= delta; int ii = 2; // a2
2342  //T(1,2) += delta; T(2,1) -= delta; int ii = 3; // a2
2343  Ti = T.T();
2344  //A = CLHEP::dot(Ti*Pm, Norm);
2345  //B = CLHEP::dot((Rt -Ti*Rm + Ti*d), Norm);
2346  A = CLHEP_dot(Ti * Pm, Norm);
2347  B = CLHEP_dot((Rt - Ti * Rm + Ti * d), Norm);
2348  double s = B / A;
2349 
2350  CLHEP::HepVector r(3, 0), p(3, 0);
2351  r = Ti * Rm - Ti * d + s * (Ti * Pm) - Rt;
2352  p = Ti * Pm - Pt;
2353 
2354  edm::LogVerbatim("GlobalTrackerMuonAlignment")
2355  << " s0 s num dsda(" << ii << ") " << s0 << " " << s << " " << (s - s0) / delta << " " << dsda(ii);
2356  // d(r,p) / d shift
2357  edm::LogVerbatim("GlobalTrackerMuonAlignment")
2358  << " -- An d(r,p)/d(" << ii << ") " << Jac(1, ii) << " " << Jac(2, ii) << " " << Jac(3, ii) << " " << Jac(4, ii)
2359  << " " << Jac(5, ii) << " " << Jac(6, ii);
2360  edm::LogVerbatim("GlobalTrackerMuonAlignment")
2361  << " Nu d(r,p)/d(" << ii << ") " << (r(1) - r0(1)) / delta << " " << (r(2) - r0(2)) / delta << " "
2362  << (r(3) - r0(3)) / delta << " " << (p(1) - p0(1)) / delta << " " << (p(2) - p0(2)) / delta << " "
2363  << (p(3) - p0(3)) / delta;
2364  // d(r,p) / d angle
2365  edm::LogVerbatim("GlobalTrackerMuonAlignment")
2366  << " -- An d(r,p)/a(" << ii << ") " << Jac(1, ii + 3) << " " << Jac(2, ii + 3) << " " << Jac(3, ii + 3) << " "
2367  << Jac(4, ii + 3) << " " << Jac(5, ii + 3) << " " << Jac(6, ii + 3);
2368  edm::LogVerbatim("GlobalTrackerMuonAlignment")
2369  << " Nu d(r,p)/a(" << ii << ") " << (r(1) - r0(1)) / delta << " " << (r(2) - r0(2)) / delta << " "
2370  << (r(3) - r0(3)) / delta << " " << (p(1) - p0(1)) / delta << " " << (p(2) - p0(2)) / delta << " "
2371  << (p(3) - p0(3)) / delta;
2372  // ----------------------------- end of check
2373 #endif
2374 
2375  return;
2376 } // end gradientGlobal
2377 
2378 // ---- calculate gradient and Hessian matrix in local parameters ------
2379 void GlobalTrackerMuonAlignment::gradientLocal(GlobalVector& GRt,
2380  GlobalVector& GPt,
2381  GlobalVector& GRm,
2382  GlobalVector& GPm,
2383  GlobalVector& GNorm,
2384  CLHEP::HepSymMatrix& covLoc,
2385  CLHEP::HepMatrix& rotLoc,
2386  CLHEP::HepVector& R0,
2387  AlgebraicVector4& LPRm) {
2388  // we search for 6D global correction vector (d, a), where
2389  // 3D vector of shihts d
2390  // 3D vector of rotation angles a
2391 
2392  bool alarm = true;
2393  //bool alarm = false;
2394  bool info = false;
2395 
2396  if (debug_)
2397  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " gradientLocal ";
2398 
2399  /*
2400  const Surface& refSurface = tsosMuon.surface();
2401 
2402  CLHEP::HepMatrix rotLoc (3,3,0);
2403  rotLoc(1,1) = refSurface.rotation().xx();
2404  rotLoc(1,2) = refSurface.rotation().xy();
2405  rotLoc(1,3) = refSurface.rotation().xz();
2406 
2407  rotLoc(2,1) = refSurface.rotation().yx();
2408  rotLoc(2,2) = refSurface.rotation().yy();
2409  rotLoc(2,3) = refSurface.rotation().yz();
2410 
2411  rotLoc(3,1) = refSurface.rotation().zx();
2412  rotLoc(3,2) = refSurface.rotation().zy();
2413  rotLoc(3,3) = refSurface.rotation().zz();
2414  */
2415 
2416  CLHEP::HepVector Rt(3), Pt(3), Rm(3), Pm(3), Norm(3);
2417  Rt(1) = GRt.x();
2418  Rt(2) = GRt.y();
2419  Rt(3) = GRt.z();
2420  Pt(1) = GPt.x();
2421  Pt(2) = GPt.y();
2422  Pt(3) = GPt.z();
2423  Rm(1) = GRm.x();
2424  Rm(2) = GRm.y();
2425  Rm(3) = GRm.z();
2426  Pm(1) = GPm.x();
2427  Pm(2) = GPm.y();
2428  Pm(3) = GPm.z();
2429  Norm(1) = GNorm.x();
2430  Norm(2) = GNorm.y();
2431  Norm(3) = GNorm.z();
2432 
2433  CLHEP::HepVector V(4), Rml(3), Pml(3), Rtl(3), Ptl(3);
2434 
2435  /*
2436  R0(1) = refSurface.position().x();
2437  R0(2) = refSurface.position().y();
2438  R0(3) = refSurface.position().z();
2439  */
2440 
2441  Rml = rotLoc * (Rm - R0);
2442  Rtl = rotLoc * (Rt - R0);
2443  Pml = rotLoc * Pm;
2444  Ptl = rotLoc * Pt;
2445 
2446  V(1) = LPRm(0) - Ptl(1) / Ptl(3);
2447  V(2) = LPRm(1) - Ptl(2) / Ptl(3);
2448  V(3) = LPRm(2) - Rtl(1);
2449  V(4) = LPRm(3) - Rtl(2);
2450 
2451  /*
2452  CLHEP::HepSymMatrix Cov(4,0), W(4,0);
2453  for(int i=1; i<=4; i++)
2454  for(int j=1; j<=i; j++){
2455  Cov(i,j) = (tsosMuon.localError().matrix()
2456  + tsosTrack.localError().matrix())(i,j);
2457  //if(i != j) Cov(i,j) = 0.;
2458  //if((i == j) && ((i==1) || (i==2))) Cov(i,j) = 100.;
2459  //if((i == j) && ((i==3) || (i==4))) Cov(i,j) = 10000.;
2460  }
2461  W = Cov;
2462  */
2463 
2464  CLHEP::HepSymMatrix W = covLoc;
2465 
2466  int ierr;
2467  W.invert(ierr);
2468  if (ierr != 0) {
2469  if (alarm)
2470  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " gradientLocal: inversion of matrix W fail ";
2471  return;
2472  }
2473 
2474  // JacobianCartesianToLocal
2475 
2476  //AlgebraicMatrix56 jacobian // differ from calculation above
2477  //= JacobianCartesianToLocal::JacobianCartesianToLocal
2478  //(refSurface, tsosTrack.localParameters()).jacobian();
2479  //for(int i=1; i<=4; i++) for(int j=1; j<=6; j++){
2480  //int j1 = j - 1; JacToLoc(i,j) = jacobian(i, j1);}
2481 
2482  CLHEP::HepMatrix JacToLoc(4, 6, 0);
2483  for (int i = 1; i <= 2; i++)
2484  for (int j = 1; j <= 3; j++) {
2485  JacToLoc(i, j + 3) = (rotLoc(i, j) - rotLoc(3, j) * Pml(i) / Pml(3)) / Pml(3);
2486  JacToLoc(i + 2, j) = rotLoc(i, j);
2487  }
2488 
2489  // JacobianCorrectionsToCartesian
2490  //double PmN = CLHEP::dot(Pm, Norm);
2491  double PmN = CLHEP_dot(Pm, Norm);
2492 
2493  CLHEP::HepMatrix Jac(6, 6, 0);
2494  for (int i = 1; i <= 3; i++)
2495  for (int j = 1; j <= 3; j++) {
2496  Jac(i, j) = Pm(i) * Norm(j) / PmN;
2497  if (i == j)
2498  Jac(i, j) -= 1.;
2499  }
2500 
2501  // dp/da
2502  Jac(4, 4) = 0.; // dpx/dax
2503  Jac(5, 4) = -Pm(3); // dpy/dax
2504  Jac(6, 4) = Pm(2); // dpz/dax
2505  Jac(4, 5) = Pm(3); // dpx/day
2506  Jac(5, 5) = 0.; // dpy/day
2507  Jac(6, 5) = -Pm(1); // dpz/day
2508  Jac(4, 6) = -Pm(2); // dpx/daz
2509  Jac(5, 6) = Pm(1); // dpy/daz
2510  Jac(6, 6) = 0.; // dpz/daz
2511 
2512  CLHEP::HepVector dsda(3);
2513  dsda(1) = (Norm(2) * Rm(3) - Norm(3) * Rm(2)) / PmN;
2514  dsda(2) = (Norm(3) * Rm(1) - Norm(1) * Rm(3)) / PmN;
2515  dsda(3) = (Norm(1) * Rm(2) - Norm(2) * Rm(1)) / PmN;
2516 
2517  // dr/da
2518  Jac(1, 4) = Pm(1) * dsda(1); // drx/dax
2519  Jac(2, 4) = -Rm(3) + Pm(2) * dsda(1); // dry/dax
2520  Jac(3, 4) = Rm(2) + Pm(3) * dsda(1); // drz/dax
2521 
2522  Jac(1, 5) = Rm(3) + Pm(1) * dsda(2); // drx/day
2523  Jac(2, 5) = Pm(2) * dsda(2); // dry/day
2524  Jac(3, 5) = -Rm(1) + Pm(3) * dsda(2); // drz/day
2525 
2526  Jac(1, 6) = -Rm(2) + Pm(1) * dsda(3); // drx/daz
2527  Jac(2, 6) = Rm(1) + Pm(2) * dsda(3); // dry/daz
2528  Jac(3, 6) = Pm(3) * dsda(3); // drz/daz
2529 
2530  // JacobianCorrectionToLocal
2531  CLHEP::HepMatrix JacCorLoc(4, 6, 0);
2532  JacCorLoc = JacToLoc * Jac;
2533 
2534  // gradient and Hessian
2535  CLHEP::HepMatrix W_Jac(6, 4, 0);
2536  W_Jac = JacCorLoc.T() * W;
2537 
2538  CLHEP::HepVector gradL(6);
2539  gradL = W_Jac * V;
2540 
2541  CLHEP::HepMatrix hessL(6, 6);
2542  hessL = JacCorLoc.T() * W * JacCorLoc;
2543 
2544  GradL += gradL;
2545  HessL += hessL;
2546 
2547  CLHEP::HepVector dLI = CLHEP::solve(HessL, -GradL);
2548  int iErI;
2549  CLHEP::HepMatrix ErrdLI = HessL.inverse(iErI);
2550  if (iErI != 0) {
2551  if (alarm)
2552  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " gradLocal Error inverse Hess matrix !!!!!!!!!!!";
2553  }
2554 
2555  if (info || debug_) {
2556  edm::LogVerbatim("GlobalTrackerMuonAlignment")
2557  << " dL " << dLI(1) << " " << dLI(2) << " " << dLI(3) << " " << dLI(4) << " " << dLI(5) << " " << dLI(6);
2558  ;
2559  edm::LogVerbatim("GlobalTrackerMuonAlignment")
2560  << " +- " << std::sqrt(ErrdLI(1, 1)) << " " << std::sqrt(ErrdLI(2, 2)) << " " << std::sqrt(ErrdLI(3, 3))
2561  << " " << std::sqrt(ErrdLI(4, 4)) << " " << std::sqrt(ErrdLI(5, 5)) << " " << std::sqrt(ErrdLI(6, 6));
2562  }
2563 
2564  if (debug_) {
2565  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug")
2566  << " dV(da3) {" << -JacCorLoc(1, 6) * 0.002 << " " << -JacCorLoc(2, 6) * 0.002 << " "
2567  << -JacCorLoc(3, 6) * 0.002 << " " << -JacCorLoc(4, 6) * 0.002 << "}";
2568  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " JacCLo {" << JacCorLoc(1, 6) << " " << JacCorLoc(2, 6)
2569  << " " << JacCorLoc(3, 6) << " " << JacCorLoc(4, 6) << "}";
2570  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug")
2571  << "Jpx/yx {" << Jac(4, 6) / Pm(3) << " " << Jac(5, 6) / Pm(3) << " " << Jac(1, 6) << " " << Jac(2, 6) << "}";
2572  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug")
2573  << "Jac(,a3){" << Jac(1, 6) << " " << Jac(2, 6) << " " << Jac(3, 6) << " " << Jac(4, 6) << " " << Jac(5, 6)
2574  << " " << Jac(6, 6);
2575  int i = 5;
2576  if (GNorm.z() > 0.95)
2577  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Ecap1 N " << GNorm;
2578  else if (GNorm.z() < -0.95)
2579  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Ecap2 N " << GNorm;
2580  else
2581  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Barrel N " << GNorm;
2582  edm::LogVerbatim("GlobalTrackerMuonAlignment")
2583  << " JacCLo(i," << i << ") = {" << JacCorLoc(1, i) << " " << JacCorLoc(2, i) << " " << JacCorLoc(3, i) << " "
2584  << JacCorLoc(4, i) << "}";
2585  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " rotLoc " << rotLoc;
2586  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " position " << R0;
2587  edm::LogVerbatim("GlobalTrackerMuonAlignment")
2588  << " Pm,l " << Pm.T() << " " << Pml(1) / Pml(3) << " " << Pml(2) / Pml(3);
2589  edm::LogVerbatim("GlobalTrackerMuonAlignment")
2590  << " Pt,l " << Pt.T() << " " << Ptl(1) / Ptl(3) << " " << Ptl(2) / Ptl(3);
2591  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " V " << V.T();
2592  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Cov \n" << covLoc;
2593  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " W*Cov " << W * covLoc;
2594  //edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " JacCarToLoc = drldrc \n" << JacobianCartesianToLocal::JacobianCartesianToLocal(refSurface, tsosTrack.localParameters()).jacobian();
2595  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " JacToLoc " << JacToLoc;
2596  }
2597 
2598 #ifdef CHECK_OF_JACOBIAN_CARTESIAN_TO_LOCAL
2599  //---------------------- check of derivatives
2600  CLHEP::HepVector V0(4, 0);
2601  V0(1) = Pml(1) / Pml(3) - Ptl(1) / Ptl(3);
2602  V0(2) = Pml(2) / Pml(3) - Ptl(2) / Ptl(3);
2603  V0(3) = Rml(1) - Rtl(1);
2604  V0(4) = Rml(2) - Rtl(2);
2605  int ii = 3;
2606  float delta = 0.01;
2607  CLHEP::HepVector V1(4, 0);
2608  if (ii <= 3) {
2609  Rm(ii) += delta;
2610  Rml = rotLoc * (Rm - R0);
2611  } else {
2612  Pm(ii - 3) += delta;
2613  Pml = rotLoc * Pm;
2614  }
2615  //if(ii <= 3) {Rt(ii) += delta; Rtl = rotLoc * (Rt - R0);}
2616  //else {Pt(ii-3) += delta; Ptl = rotLoc * Pt;}
2617  V1(1) = Pml(1) / Pml(3) - Ptl(1) / Ptl(3);
2618  V1(2) = Pml(2) / Pml(3) - Ptl(2) / Ptl(3);
2619  V1(3) = Rml(1) - Rtl(1);
2620  V1(4) = Rml(2) - Rtl(2);
2621 
2622  if (GNorm.z() > 0.95)
2623  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Ecap1 N " << GNorm;
2624  else if (GNorm.z() < -0.95)
2625  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Ecap2 N " << GNorm;
2626  else
2627  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Barrel N " << GNorm;
2628  edm::LogVerbatim("GlobalTrackerMuonAlignment")
2629  << " dldc Num (i," << ii << ") " << (V1(1) - V0(1)) / delta << " " << (V1(2) - V0(2)) / delta << " "
2630  << (V1(3) - V0(3)) / delta << " " << (V1(4) - V0(4)) / delta;
2631  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " dldc Ana (i," << ii << ") " << JacToLoc(1, ii) << " "
2632  << JacToLoc(2, ii) << " " << JacToLoc(3, ii) << " " << JacToLoc(4, ii);
2633  float dtxdpx = (rotLoc(1, 1) - rotLoc(3, 1) * Pml(1) / Pml(3)) / Pml(3);
2634  float dtydpx = (rotLoc(2, 1) - rotLoc(3, 2) * Pml(2) / Pml(3)) / Pml(3);
2635  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " dtx/dpx dty/ " << dtxdpx << " " << dtydpx;
2636 #endif
2637 
2638  return;
2639 } // end gradientLocal
2640 
2641 // ---- simulate a misalignment of muon system ------
2642 void GlobalTrackerMuonAlignment::misalignMuon(
2643  GlobalVector& GRm, GlobalVector& GPm, GlobalVector& Nl, GlobalVector& Rt, GlobalVector& Rm, GlobalVector& Pm) {
2644  CLHEP::HepVector d(3, 0), a(3, 0), Norm(3), Rm0(3), Pm0(3), Rm1(3), Pm1(3), Rmi(3), Pmi(3);
2645 
2646  d(1) = 0.0;
2647  d(2) = 0.0;
2648  d(3) = 0.0;
2649  a(1) = 0.0000;
2650  a(2) = 0.0000;
2651  a(3) = 0.0000; // zero
2652  //d(1)=0.0; d(2)=0.0; d(3)=0.0; a(1)=0.0020; a(2)=0.0000; a(3)=0.0000; // a1=.002
2653  //d(1)=0.0; d(2)=0.0; d(3)=0.0; a(1)=0.0000; a(2)=0.0020; a(3)=0.0000; // a2=.002
2654  //d(1)=0.0; d(2)=0.0; d(3)=0.0; a(1)=0.0000; a(2)=0.0000; a(3)=0.0020; // a3=.002
2655  //d(1)=1.0; d(2)=2.0; d(3)=3.0; a(1)=0.0000; a(2)=0.0000; a(3)=0.0000; // 12300
2656  //d(1)=0.0; d(2)=0.0; d(3)=0.0; a(1)=0.0000; a(2)=0.0000; a(3)=0.0020; // a3=0.002
2657  //d(1)=10.; d(2)=20.; d(3)=30.; a(1)=0.0100; a(2)=0.0200; a(3)=0.0300; // huge
2658  //d(1)=10.; d(2)=20.; d(3)=30.; a(1)=0.2000; a(2)=0.2500; a(3)=0.3000; // huge angles
2659  //d(1)=0.3; d(2)=0.4; d(3)=0.5; a(1)=0.0005; a(2)=0.0006; a(3)=0.0007; // 0345 0567
2660  //d(1)=0.3; d(2)=0.4; d(3)=0.5; a(1)=0.0001; a(2)=0.0002; a(3)=0.0003; // 0345 0123
2661  //d(1)=0.2; d(2)=0.2; d(3)=0.2; a(1)=0.0002; a(2)=0.0002; a(3)=0.0002; // 0111 0111
2662 
2663  MuGlShift = d;
2664  MuGlAngle = a;
2665 
2666  if ((d(1) == 0.) && (d(2) == 0.) && (d(3) == 0.) && (a(1) == 0.) && (a(2) == 0.) && (a(3) == 0.)) {
2667  Rm = GRm;
2668  Pm = GPm;
2669  if (debug_)
2670  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ...... without misalignment ";
2671  return;
2672  }
2673 
2674  Rm0(1) = GRm.x();
2675  Rm0(2) = GRm.y();
2676  Rm0(3) = GRm.z();
2677  Pm0(1) = GPm.x();
2678  Pm0(2) = GPm.y();
2679  Pm0(3) = GPm.z();
2680  Norm(1) = Nl.x();
2681  Norm(2) = Nl.y();
2682  Norm(3) = Nl.z();
2683 
2684  CLHEP::HepMatrix T(3, 3, 1);
2685 
2686  //T(1,2) = a(3); T(1,3) = -a(2);
2687  //T(2,1) = -a(3); T(2,3) = a(1);
2688  //T(3,1) = a(2); T(3,2) = -a(1);
2689 
2690  double s1 = std::sin(a(1)), c1 = std::cos(a(1));
2691  double s2 = std::sin(a(2)), c2 = std::cos(a(2));
2692  double s3 = std::sin(a(3)), c3 = std::cos(a(3));
2693  T(1, 1) = c2 * c3;
2694  T(1, 2) = c1 * s3 + s1 * s2 * c3;
2695  T(1, 3) = s1 * s3 - c1 * s2 * c3;
2696  T(2, 1) = -c2 * s3;
2697  T(2, 2) = c1 * c3 - s1 * s2 * s3;
2698  T(2, 3) = s1 * c3 + c1 * s2 * s3;
2699  T(3, 1) = s2;
2700  T(3, 2) = -s1 * c2;
2701  T(3, 3) = c1 * c2;
2702 
2703  //int terr;
2704  //CLHEP::HepMatrix Ti = T.inverse(terr);
2705  CLHEP::HepMatrix Ti = T.T();
2706  CLHEP::HepVector di = -Ti * d;
2707 
2708  CLHEP::HepVector ex0(3, 0), ey0(3, 0), ez0(3, 0), ex(3, 0), ey(3, 0), ez(3, 0);
2709  ex0(1) = 1.;
2710  ey0(2) = 1.;
2711  ez0(3) = 1.;
2712  ex = Ti * ex0;
2713  ey = Ti * ey0;
2714  ez = Ti * ez0;
2715 
2716  CLHEP::HepVector TiN = Ti * Norm;
2717  //double si = CLHEP::dot((Ti*Rm0 - Rm0 + di), TiN) / CLHEP::dot(Pm0, TiN);
2718  //Rm1(1) = CLHEP::dot(ex, Rm0 + si*Pm0 - di);
2719  //Rm1(2) = CLHEP::dot(ey, Rm0 + si*Pm0 - di);
2720  //Rm1(3) = CLHEP::dot(ez, Rm0 + si*Pm0 - di);
2721  double si = CLHEP_dot((Ti * Rm0 - Rm0 + di), TiN) / CLHEP_dot(Pm0, TiN);
2722  Rm1(1) = CLHEP_dot(ex, Rm0 + si * Pm0 - di);
2723  Rm1(2) = CLHEP_dot(ey, Rm0 + si * Pm0 - di);
2724  Rm1(3) = CLHEP_dot(ez, Rm0 + si * Pm0 - di);
2725  Pm1 = T * Pm0;
2726 
2727  Rm = GlobalVector(Rm1(1), Rm1(2), Rm1(3));
2728  //Pm = GlobalVector(CLHEP::dot(Pm0,ex), CLHEP::dot(Pm0, ey), CLHEP::dot(Pm0,ez));// =T*Pm0
2729  Pm = GlobalVector(CLHEP_dot(Pm0, ex), CLHEP_dot(Pm0, ey), CLHEP_dot(Pm0, ez)); // =T*Pm0
2730 
2731  if (debug_) { // ------------- debug tranformation
2732 
2733  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ----- Pm " << Pm;
2734  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " T*Pm0 " << Pm1.T();
2735  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " Ti*T*Pm0 " << (Ti * (T * Pm0)).T();
2736 
2737  //CLHEP::HepVector Rml = (Rm0 + si*Pm0 - di) + di;
2738  CLHEP::HepVector Rml = Rm1(1) * ex + Rm1(2) * ey + Rm1(3) * ez + di;
2739 
2740  CLHEP::HepVector Rt0(3);
2741  Rt0(1) = Rt.x();
2742  Rt0(2) = Rt.y();
2743  Rt0(3) = Rt.z();
2744 
2745  //double sl = CLHEP::dot(T*(Rt0 - Rml), T*Norm) / CLHEP::dot(Ti * Pm1, Norm);
2746  double sl = CLHEP_dot(T * (Rt0 - Rml), T * Norm) / CLHEP_dot(Ti * Pm1, Norm);
2747  CLHEP::HepVector Rl = Rml + sl * (Ti * Pm1);
2748 
2749  //double A = CLHEP::dot(Ti*Pm1, Norm);
2750  //double B = CLHEP::dot(Rt0, Norm) - CLHEP::dot(Ti*Rm1, Norm)
2751  //+ CLHEP::dot(Ti*d, Norm);
2752  double A = CLHEP_dot(Ti * Pm1, Norm);
2753  double B = CLHEP_dot(Rt0, Norm) - CLHEP_dot(Ti * Rm1, Norm) + CLHEP_dot(Ti * d, Norm);
2754  double s = B / A;
2755  CLHEP::HepVector Dr = Ti * Rm1 - Ti * d + s * (Ti * Pm1) - Rm0;
2756  CLHEP::HepVector Dp = Ti * Pm1 - Pm0;
2757 
2758  CLHEP::HepVector TiR = Ti * Rm0 + di;
2759  CLHEP::HepVector Ri = T * TiR + d;
2760 
2761  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " --TTi-0 " << (Ri - Rm0).T();
2762  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -- Dr " << Dr.T();
2763  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -- Dp " << Dp.T();
2764  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ex " << ex.T();
2765  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ey " << ey.T();
2766  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ez " << ez.T();
2767  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ---- T ---- " << T;
2768  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ---- Ti ---- " << Ti;
2769  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ---- d ---- " << d.T();
2770  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ---- di ---- " << di.T();
2771  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -- si sl s " << si << " " << sl << " " << s;
2772  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- si sl " << si << " " << sl;
2773  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- si s " << si << " " << s;
2774  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -- Rml-(Rm0+si*Pm0) " << (Rml - Rm0 - si * Pm0).T();
2775  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- Rm0 " << Rm0.T();
2776  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- Rm1 " << Rm1.T();
2777  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- Rmi " << Rmi.T();
2778  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --siPm " << (si * Pm0).T();
2779  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --s2Pm " << (s2 * (T * Pm1)).T();
2780  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --R1-0 " << (Rm1 - Rm0).T();
2781  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --Ri-0 " << (Rmi - Rm0).T();
2782  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " --Rl-0 " << (Rl - Rm0).T();
2783  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --Pi-0 " << (Pmi - Pm0).T();
2784  } // -------- end of debug
2785 
2786  return;
2787 
2788 } // end of misalignMuon
2789 
2790 // ---- simulate a misalignment of muon system with Local ------
2791 void GlobalTrackerMuonAlignment::misalignMuonL(GlobalVector& GRm,
2792  GlobalVector& GPm,
2793  GlobalVector& Nl,
2794  GlobalVector& Rt,
2795  GlobalVector& Rm,
2796  GlobalVector& Pm,
2797  AlgebraicVector4& Vm,
2798  TrajectoryStateOnSurface& tsosTrack,
2799  TrajectoryStateOnSurface& tsosMuon) {
2800  CLHEP::HepVector d(3, 0), a(3, 0), Norm(3), Rm0(3), Pm0(3), Rm1(3), Pm1(3), Rmi(3), Pmi(3);
2801 
2802  d(1) = 0.0;
2803  d(2) = 0.0;
2804  d(3) = 0.0;
2805  a(1) = 0.0000;
2806  a(2) = 0.0000;
2807  a(3) = 0.0000; // zero
2808  //d(1)=0.0000001; d(2)=0.0; d(3)=0.0; a(1)=0.0000; a(2)=0.0000; a(3)=0.0000; // nonzero
2809  //d(1)=0.0; d(2)=0.0; d(3)=0.0; a(1)=0.0020; a(2)=0.0000; a(3)=0.0000; // a1=.002
2810  //d(1)=0.0; d(2)=0.0; d(3)=0.0; a(1)=0.0000; a(2)=0.0020; a(3)=0.0000; // a2=.002
2811  //d(1)=0.0; d(2)=0.0; d(3)=0.0; a(1)=0.0000; a(2)=0.0000; a(3)=0.0020; // a3=.002
2812  //d(1)=0.0; d(2)=0.0; d(3)=0.0; a(1)=0.0000; a(2)=0.0000; a(3)=0.0200; // a3=.020
2813  //d(1)=1.0; d(2)=2.0; d(3)=3.0; a(1)=0.0000; a(2)=0.0000; a(3)=0.0000; // 12300
2814  //d(1)=0.0; d(2)=0.0; d(3)=0.0; a(1)=0.0000; a(2)=0.0000; a(3)=0.0020; // a3=0.002
2815  //d(1)=10.; d(2)=20.; d(3)=30.; a(1)=0.0100; a(2)=0.0200; a(3)=0.0300; // huge
2816  //d(1)=10.; d(2)=20.; d(3)=30.; a(1)=0.2000; a(2)=0.2500; a(3)=0.3000; // huge angles
2817  //d(1)=0.3; d(2)=0.4; d(3)=0.5; a(1)=0.0005; a(2)=0.0006; a(3)=0.0007; // 0345 0567
2818  //d(1)=0.3; d(2)=0.4; d(3)=0.5; a(1)=0.0001; a(2)=0.0002; a(3)=0.0003; // 0345 0123
2819  //d(1)=0.1; d(2)=0.2; d(3)=0.3; a(1)=0.0003; a(2)=0.0004; a(3)=0.0005; // 0123 0345
2820 
2821  MuGlShift = d;
2822  MuGlAngle = a;
2823 
2824  if ((d(1) == 0.) && (d(2) == 0.) && (d(3) == 0.) && (a(1) == 0.) && (a(2) == 0.) && (a(3) == 0.)) {
2825  Rm = GRm;
2826  Pm = GPm;
2827  AlgebraicVector4 Vm0;
2828  Vm = AlgebraicVector4(tsosMuon.localParameters().vector()(1),
2829  tsosMuon.localParameters().vector()(2),
2830  tsosMuon.localParameters().vector()(3),
2831  tsosMuon.localParameters().vector()(4));
2832  if (debug_)
2833  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ...... without misalignment ";
2834  return;
2835  }
2836 
2837  Rm0(1) = GRm.x();
2838  Rm0(2) = GRm.y();
2839  Rm0(3) = GRm.z();
2840  Pm0(1) = GPm.x();
2841  Pm0(2) = GPm.y();
2842  Pm0(3) = GPm.z();
2843  Norm(1) = Nl.x();
2844  Norm(2) = Nl.y();
2845  Norm(3) = Nl.z();
2846 
2847  CLHEP::HepMatrix T(3, 3, 1);
2848 
2849  //T(1,2) = a(3); T(1,3) = -a(2);
2850  //T(2,1) = -a(3); T(2,3) = a(1);
2851  //T(3,1) = a(2); T(3,2) = -a(1);
2852 
2853  double s1 = std::sin(a(1)), c1 = std::cos(a(1));
2854  double s2 = std::sin(a(2)), c2 = std::cos(a(2));
2855  double s3 = std::sin(a(3)), c3 = std::cos(a(3));
2856  T(1, 1) = c2 * c3;
2857  T(1, 2) = c1 * s3 + s1 * s2 * c3;
2858  T(1, 3) = s1 * s3 - c1 * s2 * c3;
2859  T(2, 1) = -c2 * s3;
2860  T(2, 2) = c1 * c3 - s1 * s2 * s3;
2861  T(2, 3) = s1 * c3 + c1 * s2 * s3;
2862  T(3, 1) = s2;
2863  T(3, 2) = -s1 * c2;
2864  T(3, 3) = c1 * c2;
2865 
2866  // Ti, di what we have to apply for misalignment
2867  //int terr;
2868  //CLHEP::HepMatrix Ti = T.inverse(terr);
2869  CLHEP::HepMatrix Ti = T.T();
2870  CLHEP::HepVector di = -Ti * d;
2871 
2872  CLHEP::HepVector ex0(3, 0), ey0(3, 0), ez0(3, 0), ex(3, 0), ey(3, 0), ez(3, 0);
2873  ex0(1) = 1.;
2874  ey0(2) = 1.;
2875  ez0(3) = 1.;
2876  ex = Ti * ex0;
2877  ey = Ti * ey0;
2878  ez = Ti * ez0;
2879 
2880  CLHEP::HepVector TiN = Ti * Norm;
2881  //double si = CLHEP::dot((Ti*Rm0 - Rm0 + di), TiN) / CLHEP::dot(Pm0, TiN);
2882  //Rm1(1) = CLHEP::dot(ex, Rm0 + si*Pm0 - di);
2883  //Rm1(2) = CLHEP::dot(ey, Rm0 + si*Pm0 - di);
2884  //Rm1(3) = CLHEP::dot(ez, Rm0 + si*Pm0 - di);
2885  double si = CLHEP_dot((Ti * Rm0 - Rm0 + di), TiN) / CLHEP_dot(Pm0, TiN);
2886  Rm1(1) = CLHEP_dot(ex, Rm0 + si * Pm0 - di);
2887  Rm1(2) = CLHEP_dot(ey, Rm0 + si * Pm0 - di);
2888  Rm1(3) = CLHEP_dot(ez, Rm0 + si * Pm0 - di);
2889  Pm1 = T * Pm0;
2890 
2891  // Global muon parameters after misalignment
2892  Rm = GlobalVector(Rm1(1), Rm1(2), Rm1(3));
2893  //Pm = GlobalVector(CLHEP::dot(Pm0,ex), CLHEP::dot(Pm0, ey), CLHEP::dot(Pm0,ez));// =T*Pm0
2894  Pm = GlobalVector(CLHEP_dot(Pm0, ex), CLHEP_dot(Pm0, ey), CLHEP_dot(Pm0, ez)); // =T*Pm0
2895 
2896  // Local muon parameters after misalignment
2897  const Surface& refSurface = tsosMuon.surface();
2898  CLHEP::HepMatrix Tl(3, 3, 0);
2899 
2900  Tl(1, 1) = refSurface.rotation().xx();
2901  Tl(1, 2) = refSurface.rotation().xy();
2902  Tl(1, 3) = refSurface.rotation().xz();
2903 
2904  Tl(2, 1) = refSurface.rotation().yx();
2905  Tl(2, 2) = refSurface.rotation().yy();
2906  Tl(2, 3) = refSurface.rotation().yz();
2907 
2908  Tl(3, 1) = refSurface.rotation().zx();
2909  Tl(3, 2) = refSurface.rotation().zy();
2910  Tl(3, 3) = refSurface.rotation().zz();
2911 
2912  CLHEP::HepVector R0(3, 0), newR0(3, 0), newRl(3, 0), newPl(3, 0);
2913  R0(1) = refSurface.position().x();
2914  R0(2) = refSurface.position().y();
2915  R0(3) = refSurface.position().z();
2916 
2917  newRl = Tl * (Rm1 - R0);
2918  newPl = Tl * Pm1;
2919 
2920  Vm(0) = newPl(1) / newPl(3);
2921  Vm(1) = newPl(2) / newPl(3);
2922  Vm(2) = newRl(1);
2923  Vm(3) = newRl(2);
2924 
2925 #ifdef CHECK_DERIVATIVES_LOCAL_VS_ANGLES
2926  float del = 0.0001;
2927  //int ii = 1; T(3,2) +=del; T(2,3) -=del;
2928  int ii = 2;
2929  T(3, 1) -= del;
2930  T(1, 3) += del;
2931  //int ii = 3; T(1,2) -=del; T(2,1) +=del;
2932  AlgebraicVector4 Vm0 = Vm;
2933  CLHEP::HepVector Rm10 = Rm1, Pm10 = Pm1;
2934 
2935  CLHEP::HepMatrix newTl = Tl * T;
2936  Ti = T.T();
2937  di = -Ti * d;
2938  ex = Ti * ex0;
2939  ey = Ti * ey0;
2940  ez = Ti * ez0;
2941  TiN = Ti * Norm;
2942  si = CLHEP_dot((Ti * Rm0 - Rm0 + di), TiN) / CLHEP_dot(Pm0, TiN);
2943  Rm1(1) = CLHEP_dot(ex, Rm0 + si * Pm0 - di);
2944  Rm1(2) = CLHEP_dot(ey, Rm0 + si * Pm0 - di);
2945  Rm1(3) = CLHEP_dot(ez, Rm0 + si * Pm0 - di);
2946  Pm1 = T * Pm0;
2947 
2948  newRl = Tl * (Rm1 - R0);
2949  newPl = Tl * Pm1;
2950 
2951  Vm(0) = newPl(1) / newPl(3);
2952  Vm(1) = newPl(2) / newPl(3);
2953  Vm(2) = newRl(1);
2954  Vm(3) = newRl(2);
2955  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ========= dVm/da" << ii << " " << (Vm - Vm0) * (1. / del);
2956  if (debug_)
2957  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug")
2958  << " dRm/da3 " << ((Rm1 - Rm10) * (1. / del)).T() << " " << ((Pm1 - Pm10) * (1. / del)).T();
2959 #endif
2960 
2961  if (debug_) {
2962  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " Norm " << Norm.T();
2963  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ex " << (Tl.T() * ex0).T();
2964  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ey " << (Tl.T() * ey0).T();
2965  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ez " << (Tl.T() * ez0).T();
2966 
2967  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " pxyz/p gl 0 " << (Pm0 / Pm0.norm()).T();
2968  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " pxyz/p loc0 " << (Tl * Pm0 / (Tl * Pm0).norm()).T();
2969  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " pxyz/p glob " << (Pm1 / Pm1.norm()).T();
2970  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " pxyz/p loc " << (newPl / newPl.norm()).T();
2971 
2972  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " Rot " << refSurface.rotation();
2973  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " Tl " << Tl;
2974  edm::LogVerbatim("GlobalTrackerMuonAlignment")
2975  << " ---- phi g0 g1 l0 l1 " << atan2(Pm0(2), Pm0(1)) << " " << atan2(Pm1(2), Pm1(1)) << " "
2976  << atan2((Tl * Pm0)(2), (Tl * Pm0)(1)) << " " << atan2(newPl(2), newPl(1));
2977  edm::LogVerbatim("GlobalTrackerMuonAlignment")
2978  << " ---- angl Gl Loc " << atan2(Pm1(2), Pm1(1)) - atan2(Pm0(2), Pm0(1)) << " "
2979  << atan2(newPl(2), newPl(1)) - atan2((Tl * Pm0)(2), (Tl * Pm0)(1)) << " "
2980  << atan2(newPl(3), newPl(2)) - atan2((Tl * Pm0)(3), (Tl * Pm0)(2)) << " "
2981  << atan2(newPl(1), newPl(3)) - atan2((Tl * Pm0)(1), (Tl * Pm0)(3)) << " ";
2982  }
2983 
2984  if (debug_) {
2985  CLHEP::HepMatrix newTl(3, 3, 0);
2986  CLHEP::HepVector R0(3, 0), newRl(3, 0), newPl(3, 0);
2987  newTl = Tl * Ti.T();
2988  newR0 = Ti * R0 + di;
2989 
2990  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " N " << Norm.T();
2991  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " dtxl yl " << Vm(0) - tsosMuon.localParameters().vector()(1)
2992  << " " << Vm(1) - tsosMuon.localParameters().vector()(2);
2993  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " dXl dYl " << Vm(2) - tsosMuon.localParameters().vector()(3)
2994  << " " << Vm(3) - tsosMuon.localParameters().vector()(4);
2995  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " localM " << tsosMuon.localParameters().vector();
2996  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Vm " << Vm;
2997 
2998  CLHEP::HepVector Rvc(3, 0), Pvc(3, 0), Rvg(3, 0), Pvg(3, 0);
2999  Rvc(1) = Vm(2);
3000  Rvc(2) = Vm(3);
3001  Pvc(3) = tsosMuon.localParameters().pzSign() * Pm0.norm() / std::sqrt(1 + Vm(0) * Vm(0) + Vm(1) * Vm(1));
3002  Pvc(1) = Pvc(3) * Vm(0);
3003  Pvc(2) = Pvc(3) * Vm(1);
3004 
3005  Rvg = newTl.T() * Rvc + newR0;
3006  Pvg = newTl.T() * Pvc;
3007 
3008  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " newPl " << newPl.T();
3009  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Pvc " << Pvc.T();
3010  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Rvg " << Rvg.T();
3011  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Rm1 " << Rm1.T();
3012  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Pvg " << Pvg.T();
3013  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Pm1 " << Pm1.T();
3014  }
3015 
3016  if (debug_) { // ---------- how to transform cartesian from shifted to correct
3017 
3018  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ----- Pm " << Pm;
3019  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " T*Pm0 " << Pm1.T();
3020  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " Ti*T*Pm0 " << (Ti * (T * Pm0)).T();
3021 
3022  //CLHEP::HepVector Rml = (Rm0 + si*Pm0 - di) + di;
3023  CLHEP::HepVector Rml = Rm1(1) * ex + Rm1(2) * ey + Rm1(3) * ez + di;
3024 
3025  CLHEP::HepVector Rt0(3);
3026  Rt0(1) = Rt.x();
3027  Rt0(2) = Rt.y();
3028  Rt0(3) = Rt.z();
3029 
3030  //double sl = CLHEP::dot(T*(Rt0 - Rml), T*Norm) / CLHEP::dot(Ti * Pm1, Norm);
3031  double sl = CLHEP_dot(T * (Rt0 - Rml), T * Norm) / CLHEP_dot(Ti * Pm1, Norm);
3032  CLHEP::HepVector Rl = Rml + sl * (Ti * Pm1);
3033 
3034  //double A = CLHEP::dot(Ti*Pm1, Norm);
3035  //double B = CLHEP::dot(Rt0, Norm) - CLHEP::dot(Ti*Rm1, Norm)
3036  //+ CLHEP::dot(Ti*d, Norm);
3037  double A = CLHEP_dot(Ti * Pm1, Norm);
3038  double B = CLHEP_dot(Rt0, Norm) - CLHEP_dot(Ti * Rm1, Norm) + CLHEP_dot(Ti * d, Norm);
3039  double s = B / A;
3040  CLHEP::HepVector Dr = Ti * Rm1 - Ti * d + s * (Ti * Pm1) - Rm0;
3041  CLHEP::HepVector Dp = Ti * Pm1 - Pm0;
3042 
3043  CLHEP::HepVector TiR = Ti * Rm0 + di;
3044  CLHEP::HepVector Ri = T * TiR + d;
3045 
3046  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " --TTi-0 " << (Ri - Rm0).T();
3047  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -- Dr " << Dr.T();
3048  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -- Dp " << Dp.T();
3049  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ex " << ex.T();
3050  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ey " << ey.T();
3051  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ez " << ez.T();
3052  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ---- T ---- " << T;
3053  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ---- Ti ---- " << Ti;
3054  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ---- d ---- " << d.T();
3055  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ---- di ---- " << di.T();
3056  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -- si sl s " << si << " " << sl << " " << s;
3057  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- si sl " << si << " " << sl;
3058  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- si s " << si << " " << s;
3059  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -- Rml-(Rm0+si*Pm0) " << (Rml - Rm0 - si * Pm0).T();
3060  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- Rm0 " << Rm0.T();
3061  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- Rm1 " << Rm1.T();
3062  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- Rmi " << Rmi.T();
3063  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --siPm " << (si * Pm0).T();
3064  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --s2Pm " << (s2 * (T * Pm1)).T();
3065  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --R1-0 " << (Rm1 - Rm0).T();
3066  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --Ri-0 " << (Rmi - Rm0).T();
3067  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " --Rl-0 " << (Rl - Rm0).T();
3068  edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --Pi-0 " << (Pmi - Pm0).T();
3069  } // -------- end of debug
3070 
3071  return;
3072 
3073 } // end misalignMuonL
3074 
3075 // ---- refit any direction of transient track ------
3076 void GlobalTrackerMuonAlignment::trackFitter(reco::TrackRef alongTr,
3077  reco::TransientTrack& alongTTr,
3078  PropagationDirection direction,
3079  TrajectoryStateOnSurface& trackFittedTSOS) {
3080  bool info = false;
3081  bool smooth = false;
3082 
3083  if (info)
3084  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ......... start of trackFitter ......... ";
3085  if (false && info) {
3086  this->debugTrackHit(" Tracker track hits ", alongTr);
3087  this->debugTrackHit(" Tracker TransientTrack hits ", alongTTr);
3088  }
3089 
3090  edm::OwnVector<TrackingRecHit> recHit;
3091  DetId trackDetId = DetId(alongTr->innerDetId());
3092  for (auto const& hit : alongTr->recHits())
3093  recHit.push_back(hit->clone());
3094  if (direction != alongMomentum) {
3095  edm::OwnVector<TrackingRecHit> recHitAlong = recHit;
3096  recHit.clear();
3097  for (unsigned int ihit = recHitAlong.size() - 1; ihit + 1 > 0; ihit--) {
3098  recHit.push_back(recHitAlong[ihit]);
3099  }
3100  recHitAlong.clear();
3101  }
3102  if (info)
3103  edm::LogVerbatim("GlobalTrackerMuonAlignment")
3104  << " ... Own recHit.size() " << recHit.size() << " " << trackDetId.rawId();
3105 
3106  //MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitTrack;
3107  TransientTrackingRecHit::RecHitContainer recHitTrack;
3108  for (auto const& hit : alongTr->recHits())
3109  recHitTrack.push_back(TTRHBuilder->build(hit));
3110 
3111  if (info)
3112  edm::LogVerbatim("GlobalTrackerMuonAlignment")
3113  << " ... recHitTrack.size() " << recHit.size() << " " << trackDetId.rawId() << " ======> recHitMu ";
3114 
3115  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu = recHitTrack;
3116  /*
3117  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu =
3118  MuRHBuilder->build(alongTTr.recHits().begin(), alongTTr.recHits().end());
3119  */
3120  if (info)
3121  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ...along.... recHitMu.size() " << recHitMu.size();
3122  if (direction != alongMomentum) {
3123  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMuAlong = recHitMu;
3124  recHitMu.clear();
3125  for (unsigned int ihit = recHitMuAlong.size() - 1; ihit + 1 > 0; ihit--) {
3126  recHitMu.push_back(recHitMuAlong[ihit]);
3127  }
3128  recHitMuAlong.clear();
3129  }
3130  if (info)
3131  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ...opposite... recHitMu.size() " << recHitMu.size();
3132 
3133  TrajectoryStateOnSurface firstTSOS;
3134  if (direction == alongMomentum)
3135  firstTSOS = alongTTr.innermostMeasurementState();
3136  else
3137  firstTSOS = alongTTr.outermostMeasurementState();
3138 
3139  AlgebraicSymMatrix55 CovLoc;
3140  CovLoc(0, 0) = 1. * firstTSOS.localError().matrix()(0, 0);
3141  CovLoc(1, 1) = 10. * firstTSOS.localError().matrix()(1, 1);
3142  CovLoc(2, 2) = 10. * firstTSOS.localError().matrix()(2, 2);
3143  CovLoc(3, 3) = 100. * firstTSOS.localError().matrix()(3, 3);
3144  CovLoc(4, 4) = 100. * firstTSOS.localError().matrix()(4, 4);
3145  TrajectoryStateOnSurface initialTSOS(
3146  firstTSOS.localParameters(), LocalTrajectoryError(CovLoc), firstTSOS.surface(), &*magneticField_);
3147 
3148  PTrajectoryStateOnDet PTraj = trajectoryStateTransform::persistentState(initialTSOS, trackDetId.rawId());
3149  //const TrajectorySeed seedT(*PTraj, recHit, alongMomentum);
3150  const TrajectorySeed seedT(PTraj, recHit, direction);
3151 
3152  std::vector<Trajectory> trajVec;
3153  if (direction == alongMomentum)
3154  trajVec = theFitter->fit(seedT, recHitMu, initialTSOS);
3155  else
3156  trajVec = theFitterOp->fit(seedT, recHitMu, initialTSOS);
3157 
3158  if (info) {
3159  this->debugTrajectorySOSv(" innermostTSOS of TransientTrack", alongTTr.innermostMeasurementState());
3160  this->debugTrajectorySOSv(" outermostTSOS of TransientTrack", alongTTr.outermostMeasurementState());
3161  this->debugTrajectorySOS(" initialTSOS for theFitter ", initialTSOS);
3162  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " . . . . . trajVec.size() " << trajVec.size();
3163  if (!trajVec.empty())
3164  this->debugTrajectory(" theFitter trajectory ", trajVec[0]);
3165  }
3166 
3167  if (!smooth) {
3168  if (!trajVec.empty())
3169  trackFittedTSOS = trajVec[0].lastMeasurement().updatedState();
3170  } else {
3171  std::vector<Trajectory> trajSVec;
3172  if (!trajVec.empty()) {
3173  if (direction == alongMomentum)
3174  trajSVec = theSmoother->trajectories(*(trajVec.begin()));
3175  else
3176  trajSVec = theSmootherOp->trajectories(*(trajVec.begin()));
3177  }
3178  if (info)
3179  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " . . . . trajSVec.size() " << trajSVec.size();
3180  if (!trajSVec.empty())
3181  this->debugTrajectorySOSv("smoothed geom InnermostState", trajSVec[0].geometricalInnermostState());
3182  if (!trajSVec.empty())
3183  trackFittedTSOS = trajSVec[0].geometricalInnermostState();
3184  }
3185 
3186  if (info)
3187  this->debugTrajectorySOSv(" trackFittedTSOS ", trackFittedTSOS);
3188 
3189 } // end of trackFitter
3190 
3191 // ---- refit any direction of muon transient track ------
3192 void GlobalTrackerMuonAlignment::muonFitter(reco::TrackRef alongTr,
3193  reco::TransientTrack& alongTTr,
3194  PropagationDirection direction,
3195  TrajectoryStateOnSurface& trackFittedTSOS) {
3196  bool info = false;
3197  bool smooth = false;
3198 
3199  if (info)
3200  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ......... start of muonFitter ........";
3201  if (false && info) {
3202  this->debugTrackHit(" Muon track hits ", alongTr);
3203  this->debugTrackHit(" Muon TransientTrack hits ", alongTTr);
3204  }
3205 
3206  edm::OwnVector<TrackingRecHit> recHit;
3207  DetId trackDetId = DetId(alongTr->innerDetId());
3208  for (auto const& hit : alongTr->recHits())
3209  recHit.push_back(hit->clone());
3210  if (direction != alongMomentum) {
3211  edm::OwnVector<TrackingRecHit> recHitAlong = recHit;
3212  recHit.clear();
3213  for (unsigned int ihit = recHitAlong.size() - 1; ihit + 1 > 0; ihit--) {
3214  recHit.push_back(recHitAlong[ihit]);
3215  }
3216  recHitAlong.clear();
3217  }
3218  if (info)
3219  edm::LogVerbatim("GlobalTrackerMuonAlignment")
3220  << " ... Own recHit.size() DetId==Muon " << recHit.size() << " " << trackDetId.rawId();
3221 
3222  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu =
3223  MuRHBuilder->build(alongTTr.recHitsBegin(), alongTTr.recHitsEnd());
3224  if (info)
3225  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ...along.... recHitMu.size() " << recHitMu.size();
3226  if (direction != alongMomentum) {
3227  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMuAlong = recHitMu;
3228  recHitMu.clear();
3229  for (unsigned int ihit = recHitMuAlong.size() - 1; ihit + 1 > 0; ihit--) {
3230  recHitMu.push_back(recHitMuAlong[ihit]);
3231  }
3232  recHitMuAlong.clear();
3233  }
3234  if (info)
3235  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ...opposite... recHitMu.size() " << recHitMu.size();
3236 
3237  TrajectoryStateOnSurface firstTSOS;
3238  if (direction == alongMomentum)
3239  firstTSOS = alongTTr.innermostMeasurementState();
3240  else
3241  firstTSOS = alongTTr.outermostMeasurementState();
3242 
3243  AlgebraicSymMatrix55 CovLoc;
3244  CovLoc(0, 0) = 1. * firstTSOS.localError().matrix()(0, 0);
3245  CovLoc(1, 1) = 10. * firstTSOS.localError().matrix()(1, 1);
3246  CovLoc(2, 2) = 10. * firstTSOS.localError().matrix()(2, 2);
3247  CovLoc(3, 3) = 100. * firstTSOS.localError().matrix()(3, 3);
3248  CovLoc(4, 4) = 100. * firstTSOS.localError().matrix()(4, 4);
3249  TrajectoryStateOnSurface initialTSOS(
3250  firstTSOS.localParameters(), LocalTrajectoryError(CovLoc), firstTSOS.surface(), &*magneticField_);
3251 
3252  PTrajectoryStateOnDet PTraj = trajectoryStateTransform::persistentState(initialTSOS, trackDetId.rawId());
3253  //const TrajectorySeed seedT(*PTraj, recHit, alongMomentum);
3254  const TrajectorySeed seedT(PTraj, recHit, direction);
3255 
3256  std::vector<Trajectory> trajVec;
3257  if (direction == alongMomentum)
3258  trajVec = theFitter->fit(seedT, recHitMu, initialTSOS);
3259  else
3260  trajVec = theFitterOp->fit(seedT, recHitMu, initialTSOS);
3261 
3262  if (info) {
3263  this->debugTrajectorySOSv(" innermostTSOS of TransientTrack", alongTTr.innermostMeasurementState());
3264  this->debugTrajectorySOSv(" outermostTSOS of TransientTrack", alongTTr.outermostMeasurementState());
3265  this->debugTrajectorySOS(" initialTSOS for theFitter ", initialTSOS);
3266  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " . . . . . trajVec.size() " << trajVec.size();
3267  if (!trajVec.empty())
3268  this->debugTrajectory(" theFitter trajectory ", trajVec[0]);
3269  }
3270 
3271  if (!smooth) {
3272  if (!trajVec.empty())
3273  trackFittedTSOS = trajVec[0].lastMeasurement().updatedState();
3274  } else {
3275  std::vector<Trajectory> trajSVec;
3276  if (!trajVec.empty()) {
3277  if (direction == alongMomentum)
3278  trajSVec = theSmoother->trajectories(*(trajVec.begin()));
3279  else
3280  trajSVec = theSmootherOp->trajectories(*(trajVec.begin()));
3281  }
3282  if (info)
3283  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " . . . . trajSVec.size() " << trajSVec.size();
3284  if (!trajSVec.empty())
3285  this->debugTrajectorySOSv("smoothed geom InnermostState", trajSVec[0].geometricalInnermostState());
3286  if (!trajSVec.empty())
3287  trackFittedTSOS = trajSVec[0].geometricalInnermostState();
3288  }
3289 
3290 } // end of muonFitter
3291 
3292 // ---- debug printout of hits from TransientTrack ------
3293 void GlobalTrackerMuonAlignment::debugTrackHit(const std::string title, TransientTrack& alongTr) {
3294  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ------- " << title << " --------";
3295  int nHit = 1;
3296  for (trackingRecHit_iterator i = alongTr.recHitsBegin(); i != alongTr.recHitsEnd(); i++) {
3297  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Hit " << nHit++ << " DetId " << (*i)->geographicalId().det();
3298  if ((*i)->geographicalId().det() == DetId::Tracker)
3299  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Tracker ";
3300  else if ((*i)->geographicalId().det() == DetId::Muon)
3301  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Muon ";
3302  else
3303  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Unknown ";
3304  if (!(*i)->isValid())
3305  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " not valid ";
3306  else
3307  edm::LogVerbatim("GlobalTrackerMuonAlignment");
3308  }
3309 }
3310 
3311 // ---- debug printout of hits from TrackRef ------
3312 void GlobalTrackerMuonAlignment::debugTrackHit(const std::string title, reco::TrackRef alongTr) {
3313  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ------- " << title << " --------";
3314  int nHit = 1;
3315  for (auto const& hit : alongTr->recHits()) {
3316  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Hit " << nHit++ << " DetId " << hit->geographicalId().det();
3317  if (hit->geographicalId().det() == DetId::Tracker)
3318  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Tracker ";
3319  else if (hit->geographicalId().det() == DetId::Muon)
3320  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Muon ";
3321  else
3322  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Unknown ";
3323  if (!hit->isValid())
3324  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " not valid ";
3325  else
3326  edm::LogVerbatim("GlobalTrackerMuonAlignment");
3327  }
3328 }
3329 
3330 // ---- debug printout TrajectoryStateOnSurface ------
3331 void GlobalTrackerMuonAlignment::debugTrajectorySOS(const std::string title, TrajectoryStateOnSurface& trajSOS) {
3332  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " --- " << title << " --- ";
3333  if (!trajSOS.isValid()) {
3334  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Not valid !!!! ";
3335  return;
3336  }
3337  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " R |p| " << trajSOS.globalPosition().perp() << " "
3338  << trajSOS.globalMomentum().mag() << " charge " << trajSOS.charge();
3339  edm::LogVerbatim("GlobalTrackerMuonAlignment")
3340  << " x p " << trajSOS.globalParameters().vector()(0) << " " << trajSOS.globalParameters().vector()(1) << " "
3341  << trajSOS.globalParameters().vector()(2) << " " << trajSOS.globalParameters().vector()(3) << " "
3342  << trajSOS.globalParameters().vector()(4) << " " << trajSOS.globalParameters().vector()(5);
3343  edm::LogVerbatim("GlobalTrackerMuonAlignment")
3344  << " +/- " << std::sqrt(trajSOS.cartesianError().matrix()(0, 0)) << " "
3345  << std::sqrt(trajSOS.cartesianError().matrix()(1, 1)) << " " << std::sqrt(trajSOS.cartesianError().matrix()(2, 2))
3346  << " " << std::sqrt(trajSOS.cartesianError().matrix()(3, 3)) << " "
3347  << std::sqrt(trajSOS.cartesianError().matrix()(4, 4)) << " "
3348  << std::sqrt(trajSOS.cartesianError().matrix()(5, 5));
3349  edm::LogVerbatim("GlobalTrackerMuonAlignment")
3350  << " q/p dxy/dz xy " << trajSOS.localParameters().vector()(0) << " " << trajSOS.localParameters().vector()(1)
3351  << " " << trajSOS.localParameters().vector()(2) << " " << trajSOS.localParameters().vector()(3) << " "
3352  << trajSOS.localParameters().vector()(4);
3353  edm::LogVerbatim("GlobalTrackerMuonAlignment")
3354  << " +/- error " << std::sqrt(trajSOS.localError().matrix()(0, 0)) << " "
3355  << std::sqrt(trajSOS.localError().matrix()(1, 1)) << " " << std::sqrt(trajSOS.localError().matrix()(2, 2)) << " "
3356  << std::sqrt(trajSOS.localError().matrix()(3, 3)) << " " << std::sqrt(trajSOS.localError().matrix()(4, 4)) << " ";
3357 }
3358 
3359 // ---- debug printout TrajectoryStateOnSurface ------
3360 void GlobalTrackerMuonAlignment::debugTrajectorySOSv(const std::string title, TrajectoryStateOnSurface trajSOS) {
3361  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " --- " << title << " --- ";
3362  if (!trajSOS.isValid()) {
3363  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Not valid !!!! ";
3364  return;
3365  }
3366  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " R |p| " << trajSOS.globalPosition().perp() << " "
3367  << trajSOS.globalMomentum().mag() << " charge " << trajSOS.charge();
3368  edm::LogVerbatim("GlobalTrackerMuonAlignment")
3369  << " x p " << trajSOS.globalParameters().vector()(0) << " " << trajSOS.globalParameters().vector()(1) << " "
3370  << trajSOS.globalParameters().vector()(2) << " " << trajSOS.globalParameters().vector()(3) << " "
3371  << trajSOS.globalParameters().vector()(4) << " " << trajSOS.globalParameters().vector()(5);
3372  edm::LogVerbatim("GlobalTrackerMuonAlignment")
3373  << " +/- " << std::sqrt(trajSOS.cartesianError().matrix()(0, 0)) << " "
3374  << std::sqrt(trajSOS.cartesianError().matrix()(1, 1)) << " " << std::sqrt(trajSOS.cartesianError().matrix()(2, 2))
3375  << " " << std::sqrt(trajSOS.cartesianError().matrix()(3, 3)) << " "
3376  << std::sqrt(trajSOS.cartesianError().matrix()(4, 4)) << " "
3377  << std::sqrt(trajSOS.cartesianError().matrix()(5, 5));
3378  edm::LogVerbatim("GlobalTrackerMuonAlignment")
3379  << " q/p dxy/dz xy " << trajSOS.localParameters().vector()(0) << " " << trajSOS.localParameters().vector()(1)
3380  << " " << trajSOS.localParameters().vector()(2) << " " << trajSOS.localParameters().vector()(3) << " "
3381  << trajSOS.localParameters().vector()(4);
3382  edm::LogVerbatim("GlobalTrackerMuonAlignment")
3383  << " +/- error " << std::sqrt(trajSOS.localError().matrix()(0, 0)) << " "
3384  << std::sqrt(trajSOS.localError().matrix()(1, 1)) << " " << std::sqrt(trajSOS.localError().matrix()(2, 2)) << " "
3385  << std::sqrt(trajSOS.localError().matrix()(3, 3)) << " " << std::sqrt(trajSOS.localError().matrix()(4, 4)) << " ";
3386 }
3387 
3388 // ---- debug printout Trajectory ------
3389 void GlobalTrackerMuonAlignment::debugTrajectory(const std::string title, Trajectory& traj) {
3390  edm::LogVerbatim("GlobalTrackerMuonAlignment") << "\n"
3391  << " ...... " << title << " ...... ";
3392  if (!traj.isValid()) {
3393  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Not valid !!!!!!!! ";
3394  return;
3395  }
3396  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " chi2/Nhit " << traj.chiSquared() << " / " << traj.foundHits();
3397  if (traj.direction() == alongMomentum)
3398  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " alongMomentum >>>>";
3399  else
3400  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " oppositeToMomentum <<<<";
3401  this->debugTrajectorySOSv(" firstMeasurementTSOS ", traj.firstMeasurement().updatedState());
3402  //this->debugTrajectorySOSv(" firstMeasurementTSOS ",traj.firstMeasurement().predictedState());
3403  this->debugTrajectorySOSv(" lastMeasurementTSOS ", traj.lastMeasurement().updatedState());
3404  //this->debugTrajectorySOSv(" geom InnermostState", traj.geometricalInnermostState());
3405  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " . . . . . . . . . . . . . . . . . . . . . . . . . . . . \n";
3406 }
3407 
3408 // ---- write GlobalPositionRcd ------
3409 void GlobalTrackerMuonAlignment::writeGlPosRcd(CLHEP::HepVector& paramVec) {
3410  //paramVec(1) = 0.1; paramVec(2) = 0.2; paramVec(3) = 0.3; //0123
3411  //paramVec(4) = 0.0001; paramVec(5) = 0.0002; paramVec(6) = 0.0003;
3412  //paramVec(1) = 0.3; paramVec(2) = 0.4; paramVec(3) = 0.5; //03450123
3413  //paramVec(4) = 0.0001; paramVec(5) = 0.0002; paramVec(6) = 0.0003;
3414  //paramVec(1) = 2.; paramVec(2) = 2.; paramVec(3) = 2.; //222
3415  //paramVec(4) = 0.02; paramVec(5) = 0.02; paramVec(6) = 0.02;
3416  //paramVec(1) = -10.; paramVec(2) = -20.; paramVec(3) = -30.; //102030
3417  //paramVec(4) = -0.1; paramVec(5) = -0.2; paramVec(6) = -0.3;
3418  //paramVec(1) = 0.; paramVec(2) = 0.; paramVec(3) = 1.; //1z
3419  //paramVec(4) = 0.0000; paramVec(5) = 0.0000; paramVec(6) = 0.0000;
3420 
3421  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " paramVector " << paramVec.T();
3422 
3423  CLHEP::Hep3Vector colX, colY, colZ;
3424 
3425 #ifdef NOT_EXACT_ROTATION_MATRIX
3426  colX = CLHEP::Hep3Vector(1., -paramVec(6), paramVec(5));
3427  colY = CLHEP::Hep3Vector(paramVec(6), 1., -paramVec(4));
3428  colZ = CLHEP::Hep3Vector(-paramVec(5), paramVec(4), 1.);
3429 #else
3430  double s1 = std::sin(paramVec(4)), c1 = std::cos(paramVec(4));
3431  double s2 = std::sin(paramVec(5)), c2 = std::cos(paramVec(5));
3432  double s3 = std::sin(paramVec(6)), c3 = std::cos(paramVec(6));
3433  colX = CLHEP::Hep3Vector(c2 * c3, -c2 * s3, s2);
3434  colY = CLHEP::Hep3Vector(c1 * s3 + s1 * s2 * c3, c1 * c3 - s1 * s2 * s3, -s1 * c2);
3435  colZ = CLHEP::Hep3Vector(s1 * s3 - c1 * s2 * c3, s1 * c3 + c1 * s2 * s3, c1 * c2);
3436 #endif
3437 
3438  CLHEP::HepVector param0(6, 0);
3439 
3440  Alignments globalPositions{};
3441 
3442  //Tracker
3443  AlignTransform tracker(
3444  iteratorTrackerRcd->translation(), iteratorTrackerRcd->rotation(), DetId(DetId::Tracker).rawId());
3445  // Muon
3446  CLHEP::Hep3Vector posMuGlRcd = iteratorMuonRcd->translation();
3447  CLHEP::HepRotation rotMuGlRcd = iteratorMuonRcd->rotation();
3448  CLHEP::HepEulerAngles angMuGlRcd = iteratorMuonRcd->rotation().eulerAngles();
3449  if (debug_)
3450  edm::LogVerbatim("GlobalTrackerMuonAlignment")
3451  << " Old muon Rcd Euler angles " << angMuGlRcd.phi() << " " << angMuGlRcd.theta() << " " << angMuGlRcd.psi();
3452  AlignTransform muon;
3453  if ((angMuGlRcd.phi() == 0.) && (angMuGlRcd.theta() == 0.) && (angMuGlRcd.psi() == 0.) && (posMuGlRcd.x() == 0.) &&
3454  (posMuGlRcd.y() == 0.) && (posMuGlRcd.z() == 0.)) {
3455  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " New muon parameters are stored in Rcd";
3456 
3457  AlignTransform muonNew(AlignTransform::Translation(paramVec(1), paramVec(2), paramVec(3)),
3458  AlignTransform::Rotation(colX, colY, colZ),
3459  DetId(DetId::Muon).rawId());
3460  muon = muonNew;
3461  } else if ((paramVec(1) == 0.) && (paramVec(2) == 0.) && (paramVec(3) == 0.) && (paramVec(4) == 0.) &&
3462  (paramVec(5) == 0.) && (paramVec(6) == 0.)) {
3463  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Old muon parameters are stored in Rcd";
3464 
3465  AlignTransform muonNew(iteratorMuonRcd->translation(), iteratorMuonRcd->rotation(), DetId(DetId::Muon).rawId());
3466  muon = muonNew;
3467  } else {
3468  edm::LogVerbatim("GlobalTrackerMuonAlignment") << " New + Old muon parameters are stored in Rcd";
3469  CLHEP::Hep3Vector posMuGlRcdThis = CLHEP::Hep3Vector(paramVec(1), paramVec(2), paramVec(3));
3470  CLHEP::HepRotation rotMuGlRcdThis = CLHEP::HepRotation(colX, colY, colZ);
3471  CLHEP::Hep3Vector posMuGlRcdNew = rotMuGlRcdThis * posMuGlRcd + posMuGlRcdThis;
3472  CLHEP::HepRotation rotMuGlRcdNew = rotMuGlRcdThis * rotMuGlRcd;
3473 
3474  AlignTransform muonNew(posMuGlRcdNew, rotMuGlRcdNew, DetId(DetId::Muon).rawId());
3475  muon = muonNew;
3476  }
3477 
3478  // Ecal
3479  AlignTransform ecal(iteratorEcalRcd->translation(), iteratorEcalRcd->rotation(), DetId(DetId::Ecal).rawId());
3480  // Hcal
3481  AlignTransform hcal(iteratorHcalRcd->translation(), iteratorHcalRcd->rotation(), DetId(DetId::Hcal).rawId());
3482  // Calo
3483  AlignTransform calo(AlignTransform::Translation(param0(1), param0(2), param0(3)),
3484  AlignTransform::EulerAngles(param0(4), param0(5), param0(6)),
3485  DetId(DetId::Calo).rawId());
3486 
3487  edm::LogVerbatim("GlobalTrackerMuonAlignment")
3488  << "Tracker (" << tracker.rawId() << ") at " << tracker.translation() << " " << tracker.rotation().eulerAngles();
3489  edm::LogVerbatim("GlobalTrackerMuonAlignment") << tracker.rotation();
3490 
3491  edm::LogVerbatim("GlobalTrackerMuonAlignment")
3492  << "Muon (" << muon.rawId() << ") at " << muon.translation() << " " << muon.rotation().eulerAngles();
3493  edm::LogVerbatim("GlobalTrackerMuonAlignment")
3494  << " rotations angles around x,y,z "
3495  << " ( " << -muon.rotation().zy() << " " << muon.rotation().zx() << " " << -muon.rotation().yx() << " )";
3496  edm::LogVerbatim("GlobalTrackerMuonAlignment") << muon.rotation();
3497 
3498  edm::LogVerbatim("GlobalTrackerMuonAlignment")
3499  << "Ecal (" << ecal.rawId() << ") at " << ecal.translation() << " " << ecal.rotation().eulerAngles();
3500  edm::LogVerbatim("GlobalTrackerMuonAlignment") << ecal.rotation();
3501 
3502  edm::LogVerbatim("GlobalTrackerMuonAlignment")
3503  << "Hcal (" << hcal.rawId() << ") at " << hcal.translation() << " " << hcal.rotation().eulerAngles();
3504  edm::LogVerbatim("GlobalTrackerMuonAlignment") << hcal.rotation();
3505 
3506  edm::LogVerbatim("GlobalTrackerMuonAlignment")
3507  << "Calo (" << calo.rawId() << ") at " << calo.translation() << " " << calo.rotation().eulerAngles();
3508  edm::LogVerbatim("GlobalTrackerMuonAlignment") << calo.rotation();
3509 
3510  globalPositions.m_align.push_back(tracker);
3511  globalPositions.m_align.push_back(muon);
3512  globalPositions.m_align.push_back(ecal);
3513  globalPositions.m_align.push_back(hcal);
3514  globalPositions.m_align.push_back(calo);
3515 
3516  edm::LogVerbatim("GlobalTrackerMuonAlignment") << "Uploading to the database...";
3517 
3518  edm::Service<cond::service::PoolDBOutputService> poolDbService;
3519 
3520  if (!poolDbService.isAvailable())
3521  throw cms::Exception("NotAvailable") << "PoolDBOutputService not available";
3522 
3523  // if (poolDbService->isNewTagRequest("GlobalPositionRcd")) {
3524  // poolDbService->createOneIOV<Alignments>(globalPositions, poolDbService->endOfTime(), "GlobalPositionRcd");
3525  // } else {
3526  // poolDbService->appendOneIOV<Alignments>(globalPositions, poolDbService->currentTime(), "GlobalPositionRcd");
3527  // }
3528  poolDbService->writeOneIOV<Alignments>(globalPositions, poolDbService->currentTime(), "GlobalPositionRcd");
3529  edm::LogVerbatim("GlobalTrackerMuonAlignment") << "done!";
3530 
3531  return;
3532 }
3533 
3534 //define this as a plug-in
3535 DEFINE_FWK_MODULE(GlobalTrackerMuonAlignment);
LocalTrajectoryParameters::pzSign
float pzSign() const
Sign of the z-component of the momentum in the local frame.
Definition: LocalTrajectoryParameters.h:146
rpcPointValidation_cfi.recHit
recHit
Definition: rpcPointValidation_cfi.py:7
GlobalTrackerMuonAlignment::trackingGeometry_
const GlobalTrackingGeometry * trackingGeometry_
Definition: GlobalTrackerMuonAlignment.cc:195
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:134
GlobalTrackerMuonAlignment::CLHEP_dot
double CLHEP_dot(const CLHEP::HepVector &a, const CLHEP::HepVector &b)
Definition: GlobalTrackerMuonAlignment.cc:148
GlobalTrackerMuonAlignment::watchGlobalPositionRcd_
edm::ESWatcher< GlobalPositionRcd > watchGlobalPositionRcd_
Definition: GlobalTrackerMuonAlignment.cc:202
GlobalTrackerMuonAlignment::theTrackingGeometry
edm::ESHandle< GlobalTrackingGeometry > theTrackingGeometry
Definition: GlobalTrackerMuonAlignment.cc:194
GlobalTrackerMuonAlignment::debugTrajectorySOSv
void debugTrajectorySOSv(const std::string, TrajectoryStateOnSurface)
Definition: GlobalTrackerMuonAlignment.cc:3360
ProducerED_cfi.InputTag
InputTag
Definition: ProducerED_cfi.py:5
deDxTools::esConsumes
ESGetTokenH3DDVariant esConsumes(std::string const &Record, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
GlobalTrackerMuonAlignment::smuonIsolateToken_
const edm::EDGetTokenT< reco::MuonCollection > smuonIsolateToken_
Definition: GlobalTrackerMuonAlignment.cc:183
info
static const TGPicture * info(bool iBackgroundIsBlack)
Definition: FWCollectionSummaryWidget.cc:153
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:69
mps_fire.i
i
Definition: mps_fire.py:429
SmartPropagator::propagate
TrajectoryStateOnSurface propagate(STA const &state, SUR const &surface) const
Definition: Propagator.h:50
AlgebraicMatrix33
ROOT::Math::SMatrix< double, 3, 3, ROOT::Math::MatRepStd< double, 3, 3 > > AlgebraicMatrix33
Definition: AlgebraicROOTObjects.h:41
GlobalTrackerMuonAlignment::gradientGlobalAlg
void gradientGlobalAlg(GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, AlgebraicSymMatrix66 &)
Definition: GlobalTrackerMuonAlignment.cc:2095
Trajectory::isValid
bool isValid() const
Definition: Trajectory.h:257
GlobalTrackerMuonAlignment::gradientLocal
void gradientLocal(GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, CLHEP::HepSymMatrix &, CLHEP::HepMatrix &, CLHEP::HepVector &, AlgebraicVector4 &)
Definition: GlobalTrackerMuonAlignment.cc:2379
TrajectorySmoother::trajectories
virtual TrajectoryContainer trajectories(const Trajectory &traj) const
Definition: TrajectorySmoother.h:18
B
Definition: APVGainStruct.h:7
GlobalTrackerMuonAlignment::m_ttrhBuilderToken
const edm::ESGetToken< TransientTrackingRecHitBuilder, TransientRecHitRecord > m_ttrhBuilderToken
Definition: GlobalTrackerMuonAlignment.cc:162
TkRotation::xx
T xx() const
Definition: extTkRotation.h:193
nano_mu_digi_cff.rawId
rawId
Definition: nano_mu_digi_cff.py:58
TrajectoryStateOnSurface::localError
const LocalTrajectoryError & localError() const
Definition: TrajectoryStateOnSurface.h:77
Surface::tangentPlane
virtual ConstReferenceCountingPointer< TangentPlane > tangentPlane(const GlobalPoint &) const =0
w
T w() const
Definition: extBasic3DVector.h:229
PV3DBase::z
T z() const
Definition: PV3DBase.h:61
reco::TransientTrack::recHitsEnd
trackingRecHit_iterator recHitsEnd() const
last iterator to RecHits
Definition: TransientTrack.h:126
TrackCandidateProducer_cfi.propagator
propagator
Definition: TrackCandidateProducer_cfi.py:17
MuonTransientTrackingRecHitBuilder::ConstRecHitContainer
TransientTrackingRecHit::ConstRecHitContainer ConstRecHitContainer
Definition: MuonTransientTrackingRecHitBuilder.h:15
Vector3DBase::dot
PreciseFloatType< T, U >::Type dot(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:99
TkRotation::xy
T xy() const
Definition: extTkRotation.h:194
GlobalTrackerMuonAlignment::smuonToken_
const edm::EDGetTokenT< reco::MuonCollection > smuonToken_
Definition: GlobalTrackerMuonAlignment.cc:191
SteppingHelixPropagator_cfi.SteppingHelixPropagator
SteppingHelixPropagator
Definition: SteppingHelixPropagator_cfi.py:3
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
GlobalTrackerMuonAlignment::gmuonCosmicToken_
const edm::EDGetTokenT< reco::TrackCollection > gmuonCosmicToken_
Definition: GlobalTrackerMuonAlignment.cc:186
TrajectoryStateOnSurface::globalParameters
const GlobalTrajectoryParameters & globalParameters() const
Definition: TrajectoryStateOnSurface.h:64
AlignTransform::Translation
CLHEP::Hep3Vector Translation
Definition: AlignTransform.h:18
TkRotation::zz
T zz() const
Definition: extTkRotation.h:201
reco::TrackCollection
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:14
GlobalTrackerMuonAlignment::m_globalPosToken
const edm::ESGetToken< Alignments, GlobalPositionRcd > m_globalPosToken
Definition: GlobalTrackerMuonAlignment.cc:160
trajectoryStateTransform::persistentState
PTrajectoryStateOnDet persistentState(const TrajectoryStateOnSurface &ts, unsigned int detid)
Definition: TrajectoryStateTransform.cc:14
GlobalTrackerMuonAlignment::misalignMuon
void misalignMuon(GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &)
Definition: GlobalTrackerMuonAlignment.cc:2642
GlobalTrackerMuonAlignment::iteratorHcalRcd
std::vector< AlignTransform >::const_iterator iteratorHcalRcd
Definition: GlobalTrackerMuonAlignment.cc:230
TrackingRecHit::RecHitContainer
std::vector< ConstRecHitPointer > RecHitContainer
Definition: TrackingRecHit.h:31
Trajectory::chiSquared
float chiSquared() const
Definition: Trajectory.h:241
TkRotation::yy
T yy() const
Definition: extTkRotation.h:197
GlobalTrackerMuonAlignment::gmuonIsolateToken_
const edm::EDGetTokenT< reco::TrackCollection > gmuonIsolateToken_
Definition: GlobalTrackerMuonAlignment.cc:182
HLT_2024v14_cff.muon
muon
Definition: HLT_2024v14_cff.py:9566
TrajectoryStateOnSurface::localParameters
const LocalTrajectoryParameters & localParameters() const
Definition: TrajectoryStateOnSurface.h:73
AlignTransform::EulerAngles
CLHEP::HepEulerAngles EulerAngles
Definition: AlignTransform.h:17
TkRotation::yz
T yz() const
Definition: extTkRotation.h:198
GlobalTrackerMuonAlignment::gmuonToken_
const edm::EDGetTokenT< reco::TrackCollection > gmuonToken_
Definition: GlobalTrackerMuonAlignment.cc:190
PropagationDirection
PropagationDirection
Definition: PropagationDirection.h:4
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
reco::TransientTrack::recHitsBegin
trackingRecHit_iterator recHitsBegin() const
first iterator to RecHits
Definition: TransientTrack.h:124
TrajectoryStateOnSurface::surface
const SurfaceType & surface() const
Definition: TrajectoryStateOnSurface.h:78
GlobalTrackerMuonAlignment::gradientGlobal
void gradientGlobal(GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, AlgebraicSymMatrix66 &)
Definition: GlobalTrackerMuonAlignment.cc:2187
Trajectory::foundHits
int foundHits() const
Definition: Trajectory.h:206
DiMuonV_cfg.muons
muons
the two sets of parameters below are mutually exclusive, depending if RECO or ALCARECO is used the us...
Definition: DiMuonV_cfg.py:214
Alignments::m_align
std::vector< AlignTransform > m_align
Definition: Alignments.h:19
Plane
Definition: Plane.h:16
analyze
example_stream void analyze(const edm::Event &, const edm::EventSetup &) override
Trajectory::lastMeasurement
TrajectoryMeasurement const & lastMeasurement() const
Definition: Trajectory.h:150
GlobalTrackerMuonAlignment::m_propToken
const edm::ESGetToken< Propagator, TrackingComponentsRecord > m_propToken
Definition: GlobalTrackerMuonAlignment.cc:161
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
TrajectoryStateOnSurface::cartesianError
const CartesianTrajectoryError cartesianError() const
Definition: TrajectoryStateOnSurface.h:71
reco::MuonCollection
std::vector< Muon > MuonCollection
collection of Muon objects
Definition: MuonFwd.h:9
bk::beginJob
void beginJob()
Definition: Breakpoints.cc:14
PbPb_ZMuSkimMuonDPG_cff.tracker
tracker
Definition: PbPb_ZMuSkimMuonDPG_cff.py:60
dumpMFGeometry_cfg.delta
delta
Definition: dumpMFGeometry_cfg.py:25
DiMuonV_cfg.tracks
tracks
Definition: DiMuonV_cfg.py:215
edm::eventsetup::Itr
RecordProviders::iterator Itr
Definition: EventSetupProvider.cc:359
TkRotation::zx
T zx() const
Definition: extTkRotation.h:199
MuonTransientTrackingRecHitBuilder::build
RecHitPointer build(const TrackingRecHit *p, edm::ESHandle< GlobalTrackingGeometry > trackingGeometry) const
Call the MuonTransientTrackingRecHit::specificBuild.
Definition: MuonTransientTrackingRecHitBuilder.cc:20
PV3DBase::x
T x() const
Definition: PV3DBase.h:59
PV3DBase::y
T y() const
Definition: PV3DBase.h:60
Chi2MeasurementEstimator_cfi.Chi2MeasurementEstimator
Chi2MeasurementEstimator
Definition: Chi2MeasurementEstimator_cfi.py:5
iEvent
int iEvent
Definition: GenABIO.cc:224
cms::cuda::V
uint32_t T const *__restrict__ uint32_t const *__restrict__ int32_t int Histo::index_type cudaStream_t V
Definition: HistoContainer.h:51
TrajectoryStateOnSurface::globalPosition
GlobalPoint globalPosition() const
Definition: TrajectoryStateOnSurface.h:65
GlobalTrackerMuonAlignment::debugTrackHit
void debugTrackHit(const std::string, reco::TrackRef)
Definition: GlobalTrackerMuonAlignment.cc:3312
reco::TransientTrack::outermostMeasurementState
TrajectoryStateOnSurface outermostMeasurementState() const
Definition: TransientTrack.h:84
GlobalTrackerMuonAlignment::watchTrackingGeometry_
edm::ESWatcher< GlobalTrackingGeometryRecord > watchTrackingGeometry_
Definition: GlobalTrackerMuonAlignment.cc:193
bsc_activity_cfg.ecal
ecal
Definition: bsc_activity_cfg.py:25
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:23
TransientTrackingRecHitBuilder::build
virtual RecHitPointer build(const TrackingRecHit *p) const =0
build a tracking rechit from an existing rechit
Trajectory::direction
PropagationDirection const & direction() const
Definition: Trajectory.cc:133
LocalTrajectoryParameters::vector
AlgebraicVector5 vector() const
Definition: LocalTrajectoryParameters.h:120
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
edm::OwnVector::size
size_type size() const
Definition: OwnVector.h:300
visualization-live-secondInstance_cfg.m
m
Definition: visualization-live-secondInstance_cfg.py:84
edm::OwnVector::clear
void clear()
Definition: OwnVector.h:481
TkRotation::zy
T zy() const
Definition: extTkRotation.h:200
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:64
GlobalTrackerMuonAlignment::iteratorMuonRcd
std::vector< AlignTransform >::const_iterator iteratorMuonRcd
Definition: GlobalTrackerMuonAlignment.cc:228
GlobalTrackerMuonAlignment::muonToken_
const edm::EDGetTokenT< reco::TrackCollection > muonToken_
Definition: GlobalTrackerMuonAlignment.cc:189
cond::service::PoolDBOutputService::writeOneIOV
Hash writeOneIOV(const T &payload, Time_t time, const std::string &recordName)
Definition: PoolDBOutputService.h:64
AlgebraicVector5
ROOT::Math::SVector< double, 5 > AlgebraicVector5
Definition: AlgebraicROOTObjects.h:14
GlobalTrackerMuonAlignment::analyzeTrackTrajectory
void analyzeTrackTrajectory(const edm::Event &, const edm::EventSetup &)
Definition: GlobalTrackerMuonAlignment.cc:1309
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
reco::TransientTrack::innermostMeasurementState
TrajectoryStateOnSurface innermostMeasurementState() const
Definition: TransientTrack.h:86
AlgebraicVector4
ROOT::Math::SVector< double, 4 > AlgebraicVector4
Definition: AlgebraicROOTObjects.h:13
edm::EventSetup::getHandle
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
Definition: EventSetup.h:130
GlobalTrackerMuonAlignment::trackToken_
const edm::EDGetTokenT< reco::TrackCollection > trackToken_
Definition: GlobalTrackerMuonAlignment.cc:188
ztail.d
d
Definition: ztail.py:151
GlobalTrackerMuonAlignment::trackCosmicToken_
const edm::EDGetTokenT< reco::TrackCollection > trackCosmicToken_
Definition: GlobalTrackerMuonAlignment.cc:184
cuy.ii
ii
Definition: cuy.py:589
edm::ESHandleBase::isValid
bool isValid() const
Definition: ESHandle.h:44
GlobalTrackerMuonAlignment::misalignMuonL
void misalignMuonL(GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, AlgebraicVector4 &, TrajectoryStateOnSurface &, TrajectoryStateOnSurface &)
Definition: GlobalTrackerMuonAlignment.cc:2791
GlobalTrackerMuonAlignment::analyze
void analyze(const edm::Event &, const edm::EventSetup &) override
Definition: GlobalTrackerMuonAlignment.cc:359
GlobalTrackerMuonAlignment::GlobalTrackerMuonAlignment
GlobalTrackerMuonAlignment(const edm::ParameterSet &)
Definition: GlobalTrackerMuonAlignment.cc:290
AlgebraicSymMatrix66
ROOT::Math::SMatrix< double, 6, 6, ROOT::Math::MatRepSym< double, 6 > > AlgebraicSymMatrix66
Definition: AlgebraicROOTObjects.h:24
GlobalTrackerMuonAlignment::writeGlPosRcd
void writeGlPosRcd(CLHEP::HepVector &d3)
Definition: GlobalTrackerMuonAlignment.cc:3409
GlobalTrackerMuonAlignment::iteratorEcalRcd
std::vector< AlignTransform >::const_iterator iteratorEcalRcd
Definition: GlobalTrackerMuonAlignment.cc:229
DetId
Definition: DetId.h:17
GlobalTrackerMuonAlignment::muonIsolateToken_
const edm::EDGetTokenT< reco::TrackCollection > muonIsolateToken_
Definition: GlobalTrackerMuonAlignment.cc:181
MillePedeFileConverter_cfg.out
out
Definition: MillePedeFileConverter_cfg.py:31
GlobalTrackerMuonAlignment::trackFitter
void trackFitter(reco::TrackRef, reco::TransientTrack &, PropagationDirection, TrajectoryStateOnSurface &)
Definition: GlobalTrackerMuonAlignment.cc:3076
mag
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
Definition: Basic3DVectorLD.h:130
fillDescriptions
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
GloballyPositioned::position
const PositionType & position() const
Definition: GloballyPositioned.h:36
DetId::rawId
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
GlobalTrackerMuonAlignment::TTRHBuilder
const TransientTrackingRecHitBuilder * TTRHBuilder
Definition: GlobalTrackerMuonAlignment.cc:211
SmartPropagator_cfi.SmartPropagator
SmartPropagator
Definition: SmartPropagator_cfi.py:3
HGC3DClusterGenMatchSelector_cfi.dR
dR
Definition: HGC3DClusterGenMatchSelector_cfi.py:7
AlgebraicSymMatrix55
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
Definition: AlgebraicROOTObjects.h:23
GlobalTrackerMuonAlignment::fillDescriptions
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
Definition: GlobalTrackerMuonAlignment.cc:337
b
double b
Definition: hdecay.h:120
GlobalTrackerMuonAlignment::watchMagneticFieldRecord_
edm::ESWatcher< IdealMagneticFieldRecord > watchMagneticFieldRecord_
Definition: GlobalTrackerMuonAlignment.cc:197
edm::ConfigurationDescriptions::add
void add(std::string const &label, ParameterSetDescription const &psetDescription)
Definition: ConfigurationDescriptions.cc:57
nano_mu_local_reco_cff.bool
bool
Definition: nano_mu_local_reco_cff.py:14
TrajectoryStateOnSurface::globalMomentum
GlobalVector globalMomentum() const
Definition: TrajectoryStateOnSurface.h:66
HLT_2024v14_cff.magneticField
magneticField
Definition: HLT_2024v14_cff.py:291
edm::ESWatcher::check
bool check(const edm::EventSetup &iSetup)
Definition: ESWatcher.h:57
TrajectoryMeasurement::updatedState
TrajectoryStateOnSurface const & updatedState() const
Definition: TrajectoryMeasurement.h:184
GlobalTrackerMuonAlignment::smuonCosmicToken_
const edm::EDGetTokenT< reco::MuonCollection > smuonCosmicToken_
Definition: GlobalTrackerMuonAlignment.cc:187
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:46
edm
HLT enums.
Definition: AlignableModifier.h:19
Trajectory::firstMeasurement
TrajectoryMeasurement const & firstMeasurement() const
Definition: Trajectory.h:166
a
double a
Definition: hdecay.h:121
GlobalTrackerMuonAlignment::debugTrajectory
void debugTrajectory(const std::string, Trajectory &)
Definition: GlobalTrackerMuonAlignment.cc:3389
edm::OwnVector::const_iterator
Definition: OwnVector.h:41
LocalTrajectoryError::matrix
const AlgebraicSymMatrix55 & matrix() const
Definition: LocalTrajectoryError.h:60
GlobalTrackerMuonAlignment::m_MagFieldToken
const edm::ESGetToken< MagneticField, IdealMagneticFieldRecord > m_MagFieldToken
Definition: GlobalTrackerMuonAlignment.cc:159
AlgebraicVector3
ROOT::Math::SVector< double, 3 > AlgebraicVector3
Definition: AlgebraicROOTObjects.h:12
CartesianTrajectoryError::matrix
const AlgebraicSymMatrix66 & matrix() const
Definition: CartesianTrajectoryError.h:28
TrajectoryStateOnSurface::freeState
FreeTrajectoryState const * freeState(bool withErrors=true) const
Definition: TrajectoryStateOnSurface.h:58
AlgebraicVector6
ROOT::Math::SVector< double, 6 > AlgebraicVector6
Definition: AlgebraicROOTObjects.h:15
AlgebraicSymMatrix33
ROOT::Math::SMatrix< double, 3, 3, ROOT::Math::MatRepSym< double, 3 > > AlgebraicSymMatrix33
Definition: AlgebraicROOTObjects.h:21
Plane::normalVector
GlobalVector normalVector() const
Definition: Plane.h:41
GlobalTrackerMuonAlignment::iteratorTrackerRcd
std::vector< AlignTransform >::const_iterator iteratorTrackerRcd
Definition: GlobalTrackerMuonAlignment.cc:227
GlobalTrackerMuonAlignment::trackIsolateToken_
const edm::EDGetTokenT< reco::TrackCollection > trackIsolateToken_
Definition: GlobalTrackerMuonAlignment.cc:180
HLT_2024v14_cff.R0
R0
Definition: HLT_2024v14_cff.py:7225
GloballyPositioned::rotation
const RotationType & rotation() const
Definition: GloballyPositioned.h:38
edm::Service::isAvailable
bool isAvailable() const
Definition: Service.h:40
A
Definition: APVGainStruct.h:7
l1tPhase2CaloJetEmulator_cfi.calo
calo
Definition: l1tPhase2CaloJetEmulator_cfi.py:23
GlobalTrackerMuonAlignment::m_TkGeometryToken
const edm::ESGetToken< GlobalTrackingGeometry, GlobalTrackingGeometryRecord > m_TkGeometryToken
Definition: GlobalTrackerMuonAlignment.cc:158
GlobalTrackerMuonAlignment::watchTrackingComponentsRecord_
edm::ESWatcher< TrackingComponentsRecord > watchTrackingComponentsRecord_
Definition: GlobalTrackerMuonAlignment.cc:200
T
long double T
Definition: Basic3DVectorLD.h:48
hltEgammaHLTExtra_cfi.hcal
hcal
Definition: hltEgammaHLTExtra_cfi.py:42
GlobalTrackerMuonAlignment::analyzeTrackTrack
void analyzeTrackTrack(const edm::Event &, const edm::EventSetup &)
Definition: GlobalTrackerMuonAlignment.cc:625
TkRotation::xz
T xz() const
Definition: extTkRotation.h:195
GlobalTrackerMuonAlignment::muonCosmicToken_
const edm::EDGetTokenT< reco::TrackCollection > muonCosmicToken_
Definition: GlobalTrackerMuonAlignment.cc:185
GlobalTrackerMuonAlignment::muonFitter
void muonFitter(reco::TrackRef, reco::TransientTrack &, PropagationDirection, TrajectoryStateOnSurface &)
Definition: GlobalTrackerMuonAlignment.cc:3192
TkRotation::yx
T yx() const
Definition: extTkRotation.h:196
calo
Definition: Common.h:9
AlCaHLTBitMon_ParallelJobs.p
def p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
GlobalTrackerMuonAlignment::MuRHBuilder
MuonTransientTrackingRecHitBuilder * MuRHBuilder
Definition: GlobalTrackerMuonAlignment.cc:210
GlobalVector
Global3DVector GlobalVector
Definition: GlobalVector.h:10
GlobalTrackerMuonAlignment::debugTrajectorySOS
void debugTrajectorySOS(const std::string, TrajectoryStateOnSurface &)
Definition: GlobalTrackerMuonAlignment.cc:3331
AlignTransform::Rotation
CLHEP::HepRotation Rotation
Definition: AlignTransform.h:19
Generated for CMSSW Reference Manual by doxygen 1.8.14