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