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