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