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