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