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