Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 09:31:57 

0001 // -*- C++ -*-
0002 //
0003 // Package:    GlobalTrackerMuonAlignment
0004 // Class:      GlobalTrackerMuonAlignment
0005 //
0006 /**\class GlobalTrackerMuonAlignment GlobalTrackerMuonAlignment.cc Alignment/GlobalTrackerMuonAlignment/src/GlobalTrackerMuonAlignment.cc
0007 
0008  Description: Producer of relative tracker and muon system alignment
0009 
0010  Implementation:
0011  A sample of global muons is used for the aligning tracker and muon system 
0012  relatively as "rigid bodies", i.e. determining offset and rotation (6 numbers) 
0013     
0014 */
0015 //
0016 // Original Author:  Alexandre Spiridonov
0017 //         Created:  Fri Oct 16 15:59:05 CEST 2009
0018 //
0019 // $Id: GlobalTrackerMuonAlignment.cc,v 1.11 2012/07/16 12:17:53 eulisse Exp $
0020 //
0021 
0022 // system include files
0023 #include <memory>
0024 #include <string>
0025 #include <map>
0026 #include <vector>
0027 #include <fstream>
0028 #include <typeinfo>
0029 
0030 // user include files
0031 
0032 // Framework
0033 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0034 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0035 #include "FWCore/Framework/interface/ESWatcher.h"
0036 #include "FWCore/Framework/interface/Frameworkfwd.h"
0037 #include "FWCore/Framework/interface/Event.h"
0038 #include "FWCore/Framework/interface/EventSetup.h"
0039 #include "FWCore/Framework/interface/MakerMacros.h"
0040 #include "FWCore/Utilities/interface/InputTag.h"
0041 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0042 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0043 #include "FWCore/ServiceRegistry/interface/Service.h"
0044 //#include "FWCore/Framework/interface/EventSetupRecordImplementation.h"
0045 
0046 // references
0047 #include "DataFormats/DetId/interface/DetId.h"
0048 #include "DataFormats/MuonReco/interface/Muon.h"
0049 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0050 #include "DataFormats/TrackReco/interface/Track.h"
0051 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0052 #include "DataFormats/GeometrySurface/interface/TangentPlane.h"
0053 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0054 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0055 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0056 #include "Geometry/CommonTopologies/interface/GeometryAligner.h"
0057 #include "MagneticField/Engine/interface/MagneticField.h"
0058 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0059 
0060 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBuilder.h"
0061 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0062 #include "TrackingTools/TrackRefitter/interface/TrackTransformer.h"
0063 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0064 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0065 #include "TrackingTools/GeomPropagators/interface/SmartPropagator.h"
0066 #include "TrackingTools/GeomPropagators/interface/PropagationDirectionChooser.h"
0067 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0068 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
0069 #include "TrackPropagation/RungeKutta/interface/defaultRKPropagator.h"
0070 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0071 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0072 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCartesianToLocal.h"
0073 
0074 // Database
0075 #include "FWCore/Framework/interface/ESHandle.h"
0076 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0077 
0078 // Alignment
0079 #include "CondFormats/Alignment/interface/Alignments.h"
0080 #include "CondFormats/Alignment/interface/AlignTransform.h"
0081 #include "CondFormats/AlignmentRecord/interface/GlobalPositionRcd.h"
0082 
0083 // Refit
0084 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0085 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0086 #include "TrackingTools/TrackFitters/interface/KFTrajectoryFitter.h"
0087 #include "TrackingTools/TrackFitters/interface/KFTrajectorySmoother.h"
0088 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0089 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0090 
0091 #include <TH1.h>
0092 #include <TH2.h>
0093 #include <TFile.h>
0094 
0095 #include <iostream>
0096 #include <cmath>
0097 #include <cassert>
0098 #include <fstream>
0099 #include <iomanip>
0100 
0101 using namespace edm;
0102 using namespace std;
0103 using namespace reco;
0104 
0105 //
0106 // class declaration
0107 //
0108 
0109 class GlobalTrackerMuonAlignment : public edm::one::EDAnalyzer<> {
0110 public:
0111   explicit GlobalTrackerMuonAlignment(const edm::ParameterSet&);
0112   ~GlobalTrackerMuonAlignment() override = default;
0113   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0114 
0115   void analyzeTrackTrack(const edm::Event&, const edm::EventSetup&);
0116   void analyzeTrackTrajectory(const edm::Event&, const edm::EventSetup&);
0117   void bookHist();
0118   void fitHist();
0119   void trackFitter(reco::TrackRef, reco::TransientTrack&, PropagationDirection, TrajectoryStateOnSurface&);
0120   void muonFitter(reco::TrackRef, reco::TransientTrack&, PropagationDirection, TrajectoryStateOnSurface&);
0121   void debugTrackHit(const std::string, reco::TrackRef);
0122   void debugTrackHit(const std::string, reco::TransientTrack&);
0123   void debugTrajectorySOS(const std::string, TrajectoryStateOnSurface&);
0124   void debugTrajectorySOSv(const std::string, TrajectoryStateOnSurface);
0125   void debugTrajectory(const std::string, Trajectory&);
0126 
0127   void gradientGlobal(GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&, AlgebraicSymMatrix66&);
0128   void gradientLocal(GlobalVector&,
0129                      GlobalVector&,
0130                      GlobalVector&,
0131                      GlobalVector&,
0132                      GlobalVector&,
0133                      CLHEP::HepSymMatrix&,
0134                      CLHEP::HepMatrix&,
0135                      CLHEP::HepVector&,
0136                      AlgebraicVector4&);
0137   void gradientGlobalAlg(GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&, AlgebraicSymMatrix66&);
0138   void misalignMuon(GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&);
0139   void misalignMuonL(GlobalVector&,
0140                      GlobalVector&,
0141                      GlobalVector&,
0142                      GlobalVector&,
0143                      GlobalVector&,
0144                      GlobalVector&,
0145                      AlgebraicVector4&,
0146                      TrajectoryStateOnSurface&,
0147                      TrajectoryStateOnSurface&);
0148   void writeGlPosRcd(CLHEP::HepVector& d3);
0149   inline double CLHEP_dot(const CLHEP::HepVector& a, const CLHEP::HepVector& b) {
0150     return a(1) * b(1) + a(2) * b(2) + a(3) * b(3);
0151   }
0152 
0153 private:
0154   void beginJob() override;
0155   void analyze(const edm::Event&, const edm::EventSetup&) override;
0156   void endJob() override;
0157 
0158   // ----------member data ---------------------------
0159   const edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> m_TkGeometryToken;
0160   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> m_MagFieldToken;
0161   const edm::ESGetToken<Alignments, GlobalPositionRcd> m_globalPosToken;
0162   const edm::ESGetToken<Propagator, TrackingComponentsRecord> m_propToken;
0163   const edm::ESGetToken<TransientTrackingRecHitBuilder, TransientRecHitRecord> m_ttrhBuilderToken;
0164 
0165   const edm::InputTag trackTags_;  // used to select what tracks to read from configuration file
0166   const edm::InputTag muonTags_;   // used to select what standalone muons
0167   const edm::InputTag gmuonTags_;  // used to select what global muons
0168   const edm::InputTag smuonTags_;  // used to select what selected muons
0169   const std::string propagator_;   // name of the propagator
0170   bool cosmicMuonMode_;
0171   bool isolatedMuonMode_;
0172   bool collectionCosmic;    // if true, read muon collection expected for CosmicMu
0173   bool collectionIsolated;  //                                            IsolateMu
0174   bool refitMuon_;          // if true, refit stand alone muon track
0175   bool refitTrack_;         //                tracker track
0176   const std::string rootOutFile_;
0177   const std::string txtOutFile_;
0178   const bool writeDB_;  // write results to DB
0179   const bool debug_;    // run in debugging mode
0180 
0181   const edm::EDGetTokenT<reco::TrackCollection> trackIsolateToken_;
0182   const edm::EDGetTokenT<reco::TrackCollection> muonIsolateToken_;
0183   const edm::EDGetTokenT<reco::TrackCollection> gmuonIsolateToken_;
0184   const edm::EDGetTokenT<reco::MuonCollection> smuonIsolateToken_;
0185   const edm::EDGetTokenT<reco::TrackCollection> trackCosmicToken_;
0186   const edm::EDGetTokenT<reco::TrackCollection> muonCosmicToken_;
0187   const edm::EDGetTokenT<reco::TrackCollection> gmuonCosmicToken_;
0188   const edm::EDGetTokenT<reco::MuonCollection> smuonCosmicToken_;
0189   const edm::EDGetTokenT<reco::TrackCollection> trackToken_;
0190   const edm::EDGetTokenT<reco::TrackCollection> muonToken_;
0191   const edm::EDGetTokenT<reco::TrackCollection> gmuonToken_;
0192   const edm::EDGetTokenT<reco::MuonCollection> smuonToken_;
0193 
0194   edm::ESWatcher<GlobalTrackingGeometryRecord> watchTrackingGeometry_;
0195   edm::ESHandle<GlobalTrackingGeometry> theTrackingGeometry;
0196   const GlobalTrackingGeometry* trackingGeometry_;
0197 
0198   edm::ESWatcher<IdealMagneticFieldRecord> watchMagneticFieldRecord_;
0199   const MagneticField* magneticField_;
0200 
0201   edm::ESWatcher<TrackingComponentsRecord> watchTrackingComponentsRecord_;
0202 
0203   edm::ESWatcher<GlobalPositionRcd> watchGlobalPositionRcd_;
0204   const Alignments* globalPositionRcd_;
0205 
0206   KFTrajectoryFitter* theFitter;
0207   KFTrajectorySmoother* theSmoother;
0208   KFTrajectoryFitter* theFitterOp;
0209   KFTrajectorySmoother* theSmootherOp;
0210   bool defineFitter;
0211   MuonTransientTrackingRecHitBuilder* MuRHBuilder;
0212   const TransientTrackingRecHitBuilder* TTRHBuilder;
0213 
0214   //                            // LSF for global d(3):    Inf * d = Gfr
0215   AlgebraicVector3 Gfr;      // free terms
0216   AlgebraicSymMatrix33 Inf;  // information matrix
0217   //                            // LSF for global d(6):    Hess * d = Grad
0218   CLHEP::HepVector Grad;  // gradient of the objective function in global parameters
0219   CLHEP::HepMatrix Hess;  // Hessian                  -- // ---
0220 
0221   CLHEP::HepVector GradL;  // gradient of the objective function in local parameters
0222   CLHEP::HepMatrix HessL;  // Hessian                  -- // ---
0223 
0224   int N_event;  // processed events
0225   int N_track;  // selected tracks
0226 
0227   std::vector<AlignTransform>::const_iterator
0228       iteratorTrackerRcd;  // location of Tracker in container globalPositionRcd_->m_aligm
0229   std::vector<AlignTransform>::const_iterator iteratorMuonRcd;  //              Muon
0230   std::vector<AlignTransform>::const_iterator iteratorEcalRcd;  //              Ecal
0231   std::vector<AlignTransform>::const_iterator iteratorHcalRcd;  //              Hcal
0232 
0233   CLHEP::HepVector MuGlShift;  // evaluated global muon shifts
0234   CLHEP::HepVector MuGlAngle;  // evaluated global muon angles
0235 
0236   std::string MuSelect;  // what part of muon system is selected for 1st hit
0237 
0238   ofstream OutGlobalTxt;  // output the vector of global alignment as text
0239 
0240   TFile* file;
0241   TH1F* histo;
0242   TH1F* histo2;  // outerP
0243   TH1F* histo3;  // outerLambda = PI/2-outerTheta
0244   TH1F* histo4;  // phi
0245   TH1F* histo5;  // outerR
0246   TH1F* histo6;  // outerZ
0247   TH1F* histo7;  // s
0248   TH1F* histo8;  // chi^2 of muon-track
0249 
0250   TH1F* histo11;  // |Rm-Rt|
0251   TH1F* histo12;  // Xm-Xt
0252   TH1F* histo13;  // Ym-Yt
0253   TH1F* histo14;  // Zm-Zt
0254   TH1F* histo15;  // Nxm-Nxt
0255   TH1F* histo16;  // Nym-Nyt
0256   TH1F* histo17;  // Nzm-Nzt
0257   TH1F* histo18;  // Error X of inner track
0258   TH1F* histo19;  // Error X of muon
0259   TH1F* histo20;  // Error of Xm-Xt
0260   TH1F* histo21;  // pull(Xm-Xt)
0261   TH1F* histo22;  // pull(Ym-Yt)
0262   TH1F* histo23;  // pull(Zm-Zt)
0263   TH1F* histo24;  // pull(PXm-PXt)
0264   TH1F* histo25;  // pull(PYm-Pyt)
0265   TH1F* histo26;  // pull(PZm-PZt)
0266   TH1F* histo27;  // Nx of tangent plane
0267   TH1F* histo28;  // Ny of tangent plane
0268   TH1F* histo29;  // lenght of inner track
0269   TH1F* histo30;  // lenght of muon track
0270   TH1F* histo31;  // chi2 local
0271   TH1F* histo32;  // pull(Pxm/Pzm - Pxt/Pzt) local
0272   TH1F* histo33;  // pull(Pym/Pzm - Pyt/Pzt) local
0273   TH1F* histo34;  // pull(Xm - Xt) local
0274   TH1F* histo35;  // pull(Ym - Yt) local
0275 
0276   TH2F* histo101;  // Rtrack/muon vs Ztrack/muon
0277   TH2F* histo102;  // Ytrack/muon vs Xtrack/muon
0278 };
0279 
0280 //
0281 // constants, enums and typedefs
0282 //
0283 
0284 //
0285 // static data member definitions
0286 //
0287 
0288 //
0289 // constructors and destructor
0290 //
0291 GlobalTrackerMuonAlignment::GlobalTrackerMuonAlignment(const edm::ParameterSet& iConfig)
0292     : m_TkGeometryToken(esConsumes()),
0293       m_MagFieldToken(esConsumes()),
0294       m_globalPosToken(esConsumes()),
0295       m_propToken(esConsumes(edm::ESInputTag("", iConfig.getParameter<string>("Propagator")))),
0296       m_ttrhBuilderToken(esConsumes(edm::ESInputTag("", "witTrackAngle"))),
0297       trackTags_(iConfig.getParameter<edm::InputTag>("tracks")),
0298       muonTags_(iConfig.getParameter<edm::InputTag>("muons")),
0299       gmuonTags_(iConfig.getParameter<edm::InputTag>("gmuons")),
0300       smuonTags_(iConfig.getParameter<edm::InputTag>("smuons")),
0301       cosmicMuonMode_(iConfig.getParameter<bool>("cosmics")),
0302       isolatedMuonMode_(iConfig.getParameter<bool>("isolated")),
0303       refitMuon_(iConfig.getParameter<bool>("refitmuon")),
0304       refitTrack_(iConfig.getParameter<bool>("refittrack")),
0305       rootOutFile_(iConfig.getUntrackedParameter<string>("rootOutFile")),
0306       txtOutFile_(iConfig.getUntrackedParameter<string>("txtOutFile")),
0307       writeDB_(iConfig.getUntrackedParameter<bool>("writeDB")),
0308       debug_(iConfig.getUntrackedParameter<bool>("debug")),
0309       trackIsolateToken_(consumes<reco::TrackCollection>(edm::InputTag("ALCARECOMuAlCalIsolatedMu", "TrackerOnly"))),
0310       muonIsolateToken_(consumes<reco::TrackCollection>(edm::InputTag("ALCARECOMuAlCalIsolatedMu", "StandAlone"))),
0311       gmuonIsolateToken_(consumes<reco::TrackCollection>(edm::InputTag("ALCARECOMuAlCalIsolatedMu", "GlobalMuon"))),
0312       smuonIsolateToken_(consumes<reco::MuonCollection>(edm::InputTag("ALCARECOMuAlCalIsolatedMu", "SelectedMuons"))),
0313       trackCosmicToken_(
0314           consumes<reco::TrackCollection>(edm::InputTag("ALCARECOMuAlCOMMuGlobalCosmicsMu", "TrackerOnly"))),
0315       muonCosmicToken_(
0316           consumes<reco::TrackCollection>(edm::InputTag("ALCARECOMuAlCOMMuGlobalCosmicsMu", "StandAlone"))),
0317       gmuonCosmicToken_(
0318           consumes<reco::TrackCollection>(edm::InputTag("ALCARECOMuAlCOMMuGlobalCosmicsMu", "GlobalMuon"))),
0319       smuonCosmicToken_(
0320           consumes<reco::MuonCollection>(edm::InputTag("ALCARECOMuAlCOMMuGlobalCosmicsMu", "SelectedMuons"))),
0321       trackToken_(consumes<reco::TrackCollection>(trackTags_)),
0322       muonToken_(consumes<reco::TrackCollection>(muonTags_)),
0323       gmuonToken_(consumes<reco::TrackCollection>(gmuonTags_)),
0324       smuonToken_(consumes<reco::MuonCollection>(smuonTags_)),
0325       defineFitter(true) {
0326 #ifdef NO_FALSE_FALSE_MODE
0327   if (cosmicMuonMode_ == false && isolatedMuonMode_ == false) {
0328     throw cms::Exception("BadConfig") << "Muon collection not set! "
0329                                       << "Use from GlobalTrackerMuonAlignment_test_cfg.py.";
0330   }
0331 #endif
0332   if (cosmicMuonMode_ == true && isolatedMuonMode_ == true) {
0333     throw cms::Exception("BadConfig") << "Muon collection can be either cosmic or isolated! "
0334                                       << "Please set only one to true.";
0335   }
0336 }
0337 
0338 void GlobalTrackerMuonAlignment::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0339   edm::ParameterSetDescription desc;
0340   desc.add<std::string>("propagator", "SteppingHelixPropagator");
0341   desc.add<edm::InputTag>("tracks", edm::InputTag("ALCARECOMuAlGlobalCosmics:TrackerOnly"));
0342   desc.add<edm::InputTag>("muons", edm::InputTag("ALCARECOMuAlGlobalCosmics:StandAlone"));
0343   desc.add<edm::InputTag>("gmuons", edm::InputTag("ALCARECOMuAlGlobalCosmics:GlobalMuon"));
0344   desc.add<edm::InputTag>("smuons", edm::InputTag("ALCARECOMuAlGlobalCosmics:SelectedMuons"));
0345   desc.add<bool>("isolated", false);
0346   desc.add<bool>("cosmics", false);
0347   desc.add<bool>("refitmuon", false);
0348   desc.add<bool>("refittrack", false);
0349   desc.addUntracked<std::string>("rootOutFile", "outfile.root");
0350   desc.addUntracked<std::string>("txtOutFile", "outglobal.txt");
0351   desc.addUntracked<bool>("writeDB", false);
0352   desc.addUntracked<bool>("debug", false);
0353   descriptions.add("globalTrackerMuonAlignment", desc);
0354 }
0355 //
0356 // member functions
0357 //
0358 
0359 // ------------ method called to for each event  ------------
0360 void GlobalTrackerMuonAlignment::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0361   //GlobalTrackerMuonAlignment::analyzeTrackTrack(iEvent, iSetup);
0362   GlobalTrackerMuonAlignment::analyzeTrackTrajectory(iEvent, iSetup);
0363 }
0364 
0365 // ------------ method called once for job just before starting event loop  ------------
0366 void GlobalTrackerMuonAlignment::beginJob() {
0367   N_event = 0;
0368   N_track = 0;
0369 
0370   if (cosmicMuonMode_ == true && isolatedMuonMode_ == false) {
0371     collectionCosmic = true;
0372     collectionIsolated = false;
0373   } else if (cosmicMuonMode_ == false && isolatedMuonMode_ == true) {
0374     collectionCosmic = false;
0375     collectionIsolated = true;
0376   } else {
0377     collectionCosmic = false;
0378     collectionIsolated = false;
0379   }
0380 
0381   for (int i = 0; i <= 2; i++) {
0382     Gfr(i) = 0.;
0383     for (int j = 0; j <= 2; j++) {
0384       Inf(i, j) = 0.;
0385     }
0386   }
0387 
0388   Grad = CLHEP::HepVector(6, 0);
0389   Hess = CLHEP::HepSymMatrix(6, 0);
0390 
0391   GradL = CLHEP::HepVector(6, 0);
0392   HessL = CLHEP::HepSymMatrix(6, 0);
0393 
0394   // histograms
0395   TDirectory* dirsave = gDirectory;
0396 
0397   file = new TFile(rootOutFile_.c_str(), "recreate");
0398   const bool oldAddDir = TH1::AddDirectoryStatus();
0399 
0400   TH1::AddDirectory(true);
0401 
0402   this->bookHist();
0403 
0404   TH1::AddDirectory(oldAddDir);
0405   dirsave->cd();
0406 }
0407 
0408 // ------------ method called once each job just after ending the event loop  ------------
0409 void GlobalTrackerMuonAlignment::endJob() {
0410   bool alarm = false;
0411 
0412   this->fitHist();
0413 
0414   AlgebraicVector3 d(0., 0., 0.);  // ------------  alignmnet  Global Algebraic
0415 
0416   AlgebraicSymMatrix33 InfI;  // inverse it
0417   for (int i = 0; i <= 2; i++)
0418     for (int j = 0; j <= 2; j++) {
0419       if (j < i)
0420         continue;
0421       InfI(i, j) += Inf(i, j);
0422     }
0423   bool ierr = !InfI.Invert();
0424   if (ierr) {
0425     if (alarm)
0426       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Error inverse  Inf matrix !!!!!!!!!!!";
0427   }
0428 
0429   for (int i = 0; i <= 2; i++)
0430     for (int k = 0; k <= 2; k++)
0431       d(i) -= InfI(i, k) * Gfr(k);
0432   // end of Global Algebraic
0433 
0434   //                                ---------------  alignment Global CLHEP
0435   CLHEP::HepVector d3 = CLHEP::solve(Hess, -Grad);
0436   int iEr3;
0437   CLHEP::HepMatrix Errd3 = Hess.inverse(iEr3);
0438   if (iEr3 != 0) {
0439     if (alarm)
0440       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " endJob Error inverse Hess matrix !!!!!!!!!!!";
0441   }
0442   // end of Global CLHEP
0443 
0444   //                                ----------------- alignment Local CLHEP
0445   CLHEP::HepVector dLI = CLHEP::solve(HessL, -GradL);
0446   int iErI;
0447   CLHEP::HepMatrix ErrdLI = HessL.inverse(iErI);
0448   if (iErI != 0) {
0449     if (alarm)
0450       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " endJob Error inverse HessL matrix !!!!!!!!!!!";
0451   }
0452   // end of Local CLHEP
0453 
0454   // printout of final parameters
0455   edm::LogVerbatim("GlobalTrackerMuonAlignment")
0456       << " ---- " << N_event << " event " << N_track << " tracks " << MuSelect << " ---- ";
0457   if (collectionIsolated)
0458     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ALCARECOMuAlCalIsolatedMu";
0459   else if (collectionCosmic)
0460     edm::LogVerbatim("GlobalTrackerMuonAlignment") << "  ALCARECOMuAlGlobalCosmics";
0461   else
0462     edm::LogVerbatim("GlobalTrackerMuonAlignment") << smuonTags_;
0463   edm::LogVerbatim("GlobalTrackerMuonAlignment")
0464       << " Similated shifts[cm] " << MuGlShift(1) << " " << MuGlShift(2) << " " << MuGlShift(3) << " "
0465       << " angles[rad] " << MuGlAngle(1) << " " << MuGlAngle(2) << " " << MuGlAngle(3) << " ";
0466   edm::LogVerbatim("GlobalTrackerMuonAlignment") << " d    " << -d;
0467 
0468   edm::LogVerbatim("GlobalTrackerMuonAlignment")
0469       << " +-   " << std::sqrt(InfI(0, 0)) << " " << std::sqrt(InfI(1, 1)) << " " << std::sqrt(InfI(2, 2));
0470   edm::LogVerbatim("GlobalTrackerMuonAlignment")
0471       << " dG   " << d3(1) << " " << d3(2) << " " << d3(3) << " " << d3(4) << " " << d3(5) << " " << d3(6);
0472 
0473   edm::LogVerbatim("GlobalTrackerMuonAlignment")
0474       << " +-   " << std::sqrt(Errd3(1, 1)) << " " << std::sqrt(Errd3(2, 2)) << " " << std::sqrt(Errd3(3, 3)) << " "
0475       << std::sqrt(Errd3(4, 4)) << " " << std::sqrt(Errd3(5, 5)) << " " << std::sqrt(Errd3(6, 6));
0476   edm::LogVerbatim("GlobalTrackerMuonAlignment")
0477       << " dL   " << dLI(1) << " " << dLI(2) << " " << dLI(3) << " " << dLI(4) << " " << dLI(5) << " " << dLI(6);
0478 
0479   edm::LogVerbatim("GlobalTrackerMuonAlignment")
0480       << " +-   " << std::sqrt(ErrdLI(1, 1)) << " " << std::sqrt(ErrdLI(2, 2)) << " " << std::sqrt(ErrdLI(3, 3)) << " "
0481       << std::sqrt(ErrdLI(4, 4)) << " " << std::sqrt(ErrdLI(5, 5)) << " " << std::sqrt(ErrdLI(6, 6));
0482 
0483   // what do we write to DB
0484   CLHEP::HepVector vectorToDb(6, 0), vectorErrToDb(6, 0);
0485   //vectorToDb = d3;
0486   //for(unsigned int i=1; i<=6; i++) vectorErrToDb(i) = std::sqrt(Errd3(i,i));
0487   vectorToDb = -dLI;
0488   for (unsigned int i = 1; i <= 6; i++)
0489     vectorErrToDb(i) = std::sqrt(ErrdLI(i, i));
0490 
0491   // write histograms to root file
0492   file->Write();
0493   file->Close();
0494 
0495   // write global parameters to text file
0496   OutGlobalTxt.open(txtOutFile_.c_str(), ios::out);
0497   if (!OutGlobalTxt.is_open())
0498     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " outglobal.txt is not open !!!!!";
0499   else {
0500     OutGlobalTxt << vectorToDb(1) << " " << vectorToDb(2) << " " << vectorToDb(3) << " " << vectorToDb(4) << " "
0501                  << vectorToDb(5) << " " << vectorToDb(6) << " muon Global.\n";
0502     OutGlobalTxt << vectorErrToDb(1) << " " << vectorErrToDb(1) << " " << vectorErrToDb(1) << " " << vectorErrToDb(1)
0503                  << " " << vectorErrToDb(1) << " " << vectorErrToDb(1) << " errors.\n";
0504     OutGlobalTxt << N_event << " events are processed.\n";
0505 
0506     if (collectionIsolated)
0507       OutGlobalTxt << "ALCARECOMuAlCalIsolatedMu.\n";
0508     else if (collectionCosmic)
0509       OutGlobalTxt << "  ALCARECOMuAlGlobalCosmics.\n";
0510     else
0511       OutGlobalTxt << smuonTags_ << ".\n";
0512     OutGlobalTxt.close();
0513     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Write to the file outglobal.txt done  ";
0514   }
0515 
0516   // write new GlobalPositionRcd to DB
0517   if (debug_)
0518     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " writeBD_ " << writeDB_;
0519   if (writeDB_)
0520     GlobalTrackerMuonAlignment::writeGlPosRcd(vectorToDb);
0521 }
0522 
0523 // ------------ book histogram  ------------
0524 void GlobalTrackerMuonAlignment::bookHist() {
0525   double PI = 3.1415927;
0526   histo = new TH1F("Pt", "pt", 1000, 0, 100);
0527   histo2 = new TH1F("P", "P [GeV/c]", 400, 0., 400.);
0528   histo2->GetXaxis()->SetTitle("momentum [GeV/c]");
0529   histo3 = new TH1F("outerLambda", "#lambda outer", 100, -PI / 2., PI / 2.);
0530   histo3->GetXaxis()->SetTitle("#lambda outer");
0531   histo4 = new TH1F("phi", "#phi [rad]", 100, -PI, PI);
0532   histo4->GetXaxis()->SetTitle("#phi [rad]");
0533   histo5 = new TH1F("Rmuon", "inner muon hit R [cm]", 100, 0., 800.);
0534   histo5->GetXaxis()->SetTitle("R of muon [cm]");
0535   histo6 = new TH1F("Zmuon", "inner muon hit Z[cm]", 100, -1000., 1000.);
0536   histo6->GetXaxis()->SetTitle("Z of muon [cm]");
0537   histo7 = new TH1F("(Pm-Pt)/Pt", " (Pmuon-Ptrack)/Ptrack", 100, -2., 2.);
0538   histo7->GetXaxis()->SetTitle("(Pmuon-Ptrack)/Ptrack");
0539   histo8 = new TH1F("chi muon-track", "#chi^{2}(muon-track)", 1000, 0., 1000.);
0540   histo8->GetXaxis()->SetTitle("#chi^{2} of muon w.r.t. propagated track");
0541   histo11 = new TH1F("distance muon-track", "distance muon w.r.t track [cm]", 100, 0., 30.);
0542   histo11->GetXaxis()->SetTitle("distance of muon w.r.t. track [cm]");
0543   histo12 = new TH1F("Xmuon-Xtrack", "Xmuon-Xtrack [cm]", 200, -20., 20.);
0544   histo12->GetXaxis()->SetTitle("Xmuon - Xtrack [cm]");
0545   histo13 = new TH1F("Ymuon-Ytrack", "Ymuon-Ytrack [cm]", 200, -20., 20.);
0546   histo13->GetXaxis()->SetTitle("Ymuon - Ytrack [cm]");
0547   histo14 = new TH1F("Zmuon-Ztrack", "Zmuon-Ztrack [cm]", 200, -20., 20.);
0548   histo14->GetXaxis()->SetTitle("Zmuon-Ztrack [cm]");
0549   histo15 = new TH1F("NXmuon-NXtrack", "NXmuon-NXtrack [rad]", 200, -.1, .1);
0550   histo15->GetXaxis()->SetTitle("N_{X}(muon)-N_{X}(track) [rad]");
0551   histo16 = new TH1F("NYmuon-NYtrack", "NYmuon-NYtrack [rad]", 200, -.1, .1);
0552   histo16->GetXaxis()->SetTitle("N_{Y}(muon)-N_{Y}(track) [rad]");
0553   histo17 = new TH1F("NZmuon-NZtrack", "NZmuon-NZtrack [rad]", 200, -.1, .1);
0554   histo17->GetXaxis()->SetTitle("N_{Z}(muon)-N_{Z}(track) [rad]");
0555   histo18 = new TH1F("expected error of Xinner", "outer hit of inner tracker", 100, 0, .01);
0556   histo18->GetXaxis()->SetTitle("expected error of Xinner [cm]");
0557   histo19 = new TH1F("expected error of Xmuon", "inner hit of muon", 100, 0, .1);
0558   histo19->GetXaxis()->SetTitle("expected error of Xmuon [cm]");
0559   histo20 = new TH1F("expected error of Xmuon-Xtrack", "muon w.r.t. propagated track", 100, 0., 10.);
0560   histo20->GetXaxis()->SetTitle("expected error of Xmuon-Xtrack [cm]");
0561   histo21 = new TH1F("pull of Xmuon-Xtrack", "pull of Xmuon-Xtrack", 100, -10., 10.);
0562   histo21->GetXaxis()->SetTitle("(Xmuon-Xtrack)/expected error");
0563   histo22 = new TH1F("pull of Ymuon-Ytrack", "pull of Ymuon-Ytrack", 100, -10., 10.);
0564   histo22->GetXaxis()->SetTitle("(Ymuon-Ytrack)/expected error");
0565   histo23 = new TH1F("pull of Zmuon-Ztrack", "pull of Zmuon-Ztrack", 100, -10., 10.);
0566   histo23->GetXaxis()->SetTitle("(Zmuon-Ztrack)/expected error");
0567   histo24 = new TH1F("pull of PXmuon-PXtrack", "pull of PXmuon-PXtrack", 100, -10., 10.);
0568   histo24->GetXaxis()->SetTitle("(P_{X}(muon)-P_{X}(track))/expected error");
0569   histo25 = new TH1F("pull of PYmuon-PYtrack", "pull of PYmuon-PYtrack", 100, -10., 10.);
0570   histo25->GetXaxis()->SetTitle("(P_{Y}(muon)-P_{Y}(track))/expected error");
0571   histo26 = new TH1F("pull of PZmuon-PZtrack", "pull of PZmuon-PZtrack", 100, -10., 10.);
0572   histo26->GetXaxis()->SetTitle("(P_{Z}(muon)-P_{Z}(track))/expected error");
0573   histo27 = new TH1F("N_x", "Nx of tangent plane", 120, -1.1, 1.1);
0574   histo27->GetXaxis()->SetTitle("normal vector projection N_{X}");
0575   histo28 = new TH1F("N_y", "Ny of tangent plane", 120, -1.1, 1.1);
0576   histo28->GetXaxis()->SetTitle("normal vector projection N_{Y}");
0577   histo29 = new TH1F("lenght of track", "lenght of track", 200, 0., 400);
0578   histo29->GetXaxis()->SetTitle("lenght of track [cm]");
0579   histo30 = new TH1F("lenght of muon", "lenght of muon", 200, 0., 800);
0580   histo30->GetXaxis()->SetTitle("lenght of muon [cm]");
0581 
0582   histo31 = new TH1F("local chi muon-track", "#local chi^{2}(muon-track)", 1000, 0., 1000.);
0583   histo31->GetXaxis()->SetTitle("#local chi^{2} of muon w.r.t. propagated track");
0584   histo32 = new TH1F("pull of Px/Pz local", "pull of Px/Pz local", 100, -10., 10.);
0585   histo32->GetXaxis()->SetTitle("local (Px/Pz(muon) - Px/Pz(track))/expected error");
0586   histo33 = new TH1F("pull of Py/Pz local", "pull of Py/Pz local", 100, -10., 10.);
0587   histo33->GetXaxis()->SetTitle("local (Py/Pz(muon) - Py/Pz(track))/expected error");
0588   histo34 = new TH1F("pull of X local", "pull of X local", 100, -10., 10.);
0589   histo34->GetXaxis()->SetTitle("local (Xmuon - Xtrack)/expected error");
0590   histo35 = new TH1F("pull of Y local", "pull of Y local", 100, -10., 10.);
0591   histo35->GetXaxis()->SetTitle("local (Ymuon - Ytrack)/expected error");
0592 
0593   histo101 = new TH2F("Rtr/mu vs Ztr/mu", "hit of track/muon", 100, -800., 800., 100, 0., 600.);
0594   histo101->GetXaxis()->SetTitle("Z of track/muon [cm]");
0595   histo101->GetYaxis()->SetTitle("R of track/muon [cm]");
0596   histo102 = new TH2F("Ytr/mu vs Xtr/mu", "hit of track/muon", 100, -600., 600., 100, -600., 600.);
0597   histo102->GetXaxis()->SetTitle("X of track/muon [cm]");
0598   histo102->GetYaxis()->SetTitle("Y of track/muon [cm]");
0599 }
0600 
0601 // ------------ fit histogram  ------------
0602 void GlobalTrackerMuonAlignment::fitHist() {
0603   histo7->Fit("gaus", "Q");
0604 
0605   histo12->Fit("gaus", "Q");
0606   histo13->Fit("gaus", "Q");
0607   histo14->Fit("gaus", "Q");
0608   histo15->Fit("gaus", "Q");
0609   histo16->Fit("gaus", "Q");
0610   histo17->Fit("gaus", "Q");
0611 
0612   histo21->Fit("gaus", "Q");
0613   histo22->Fit("gaus", "Q");
0614   histo23->Fit("gaus", "Q");
0615   histo24->Fit("gaus", "Q");
0616   histo25->Fit("gaus", "Q");
0617   histo26->Fit("gaus", "Q");
0618 
0619   histo32->Fit("gaus", "Q");
0620   histo33->Fit("gaus", "Q");
0621   histo34->Fit("gaus", "Q");
0622   histo35->Fit("gaus", "Q");
0623 }
0624 
0625 // ------------ method to analyze recoTrack & recoMuon of Global Muon  ------------
0626 void GlobalTrackerMuonAlignment::analyzeTrackTrack(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0627   using namespace edm;
0628   using namespace reco;
0629   //debug_ = true;
0630   bool info = false;
0631   bool alarm = false;
0632   //bool alarm = true;
0633   double PI = 3.1415927;
0634 
0635   cosmicMuonMode_ = true;                // true: both Cosmic and IsolatedMuon are treated with 1,2 tracks
0636   isolatedMuonMode_ = !cosmicMuonMode_;  //true: only IsolatedMuon are treated with 1 track
0637 
0638   //if(N_event == 1)
0639   //GlobalPositionRcdRead1::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup);
0640 
0641   N_event++;
0642   if (info || debug_) {
0643     edm::LogVerbatim("GlobalTrackerMuonAlignment") << "----- event " << N_event << " -- tracks " << N_track << " ---";
0644     if (cosmicMuonMode_)
0645       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " treated as CosmicMu ";
0646     if (isolatedMuonMode_)
0647       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " treated as IsolatedMu ";
0648     edm::LogVerbatim("GlobalTrackerMuonAlignment");
0649   }
0650 
0651   const edm::Handle<reco::TrackCollection>& tracks =
0652       ((collectionIsolated)
0653            ? iEvent.getHandle(trackIsolateToken_)
0654            : ((collectionCosmic) ? iEvent.getHandle(trackCosmicToken_) : iEvent.getHandle(trackToken_)));
0655   const edm::Handle<reco::TrackCollection>& muons =
0656       ((collectionIsolated) ? iEvent.getHandle(muonIsolateToken_)
0657                             : ((collectionCosmic) ? iEvent.getHandle(muonCosmicToken_) : iEvent.getHandle(muonToken_)));
0658   const edm::Handle<reco::TrackCollection>& gmuons =
0659       ((collectionIsolated)
0660            ? iEvent.getHandle(gmuonIsolateToken_)
0661            : ((collectionCosmic) ? iEvent.getHandle(gmuonCosmicToken_) : iEvent.getHandle(gmuonToken_)));
0662   const edm::Handle<reco::MuonCollection>& smuons =
0663       ((collectionIsolated)
0664            ? iEvent.getHandle(smuonIsolateToken_)
0665            : ((collectionCosmic) ? iEvent.getHandle(smuonCosmicToken_) : iEvent.getHandle(smuonToken_)));
0666 
0667   if (debug_) {
0668     edm::LogVerbatim("GlobalTrackerMuonAlignment")
0669         << " ievBunch " << iEvent.bunchCrossing() << " runN " << (int)iEvent.run();
0670     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " N tracks s/amu gmu selmu " << tracks->size() << " "
0671                                                    << muons->size() << " " << gmuons->size() << " " << smuons->size();
0672     for (MuonCollection::const_iterator itMuon = smuons->begin(); itMuon != smuons->end(); ++itMuon) {
0673       edm::LogVerbatim("GlobalTrackerMuonAlignment")
0674           << " is isolatValid Matches " << itMuon->isIsolationValid() << " " << itMuon->isMatchesValid();
0675     }
0676   }
0677 
0678   if (isolatedMuonMode_) {  // ---- Only 1 Isolated Muon --------
0679     if (tracks->size() != 1)
0680       return;
0681     if (muons->size() != 1)
0682       return;
0683     if (gmuons->size() != 1)
0684       return;
0685     if (smuons->size() != 1)
0686       return;
0687   }
0688 
0689   if (cosmicMuonMode_) {  // ---- 1,2 Cosmic Muon --------
0690     if (smuons->size() > 2)
0691       return;
0692     if (tracks->size() != smuons->size())
0693       return;
0694     if (muons->size() != smuons->size())
0695       return;
0696     if (gmuons->size() != smuons->size())
0697       return;
0698   }
0699 
0700   // ok mc_isolated_mu
0701   //edm::ESHandle<TrackerGeometry> trackerGeometry;
0702   //iSetup.get<TrackerDigiGeometryRecord>().get(trackerGeometry);
0703   // ok mc_isolated_mu
0704   //edm::ESHandle<DTGeometry> dtGeometry;
0705   //iSetup.get<MuonGeometryRecord>().get(dtGeometry);
0706   // don't work
0707   //edm::ESHandle<CSCGeometry> cscGeometry;
0708   //iSetup.get<MuonGeometryRecord>().get(cscGeometry);
0709 
0710   if (watchTrackingGeometry_.check(iSetup) || !trackingGeometry_) {
0711     edm::ESHandle<GlobalTrackingGeometry> trackingGeometry = iSetup.getHandle(m_TkGeometryToken);
0712     trackingGeometry_ = &*trackingGeometry;
0713   }
0714 
0715   if (watchMagneticFieldRecord_.check(iSetup) || !magneticField_) {
0716     edm::ESHandle<MagneticField> magneticField = iSetup.getHandle(m_MagFieldToken);
0717     magneticField_ = &*magneticField;
0718   }
0719 
0720   if (watchGlobalPositionRcd_.check(iSetup) || !globalPositionRcd_) {
0721     edm::ESHandle<Alignments> globalPositionRcd = iSetup.getHandle(m_globalPosToken);
0722     globalPositionRcd_ = &*globalPositionRcd;
0723     for (std::vector<AlignTransform>::const_iterator i = globalPositionRcd_->m_align.begin();
0724          i != globalPositionRcd_->m_align.end();
0725          ++i) {
0726       if (DetId(DetId::Tracker).rawId() == i->rawId())
0727         iteratorTrackerRcd = i;
0728       if (DetId(DetId::Muon).rawId() == i->rawId())
0729         iteratorMuonRcd = i;
0730       if (DetId(DetId::Ecal).rawId() == i->rawId())
0731         iteratorEcalRcd = i;
0732       if (DetId(DetId::Hcal).rawId() == i->rawId())
0733         iteratorHcalRcd = i;
0734     }
0735     if (true || debug_) {
0736       edm::LogVerbatim("GlobalTrackerMuonAlignment") << "=== iteratorMuonRcd " << iteratorMuonRcd->rawId() << " ====\n"
0737                                                      << " translation " << iteratorMuonRcd->translation() << "\n"
0738                                                      << " angles " << iteratorMuonRcd->rotation().eulerAngles();
0739       edm::LogVerbatim("GlobalTrackerMuonAlignment") << iteratorMuonRcd->rotation();
0740     }
0741   }  // end of GlobalPositionRcd
0742 
0743   edm::ESHandle<Propagator> propagator = iSetup.getHandle(m_propToken);
0744 
0745   SteppingHelixPropagator alongStHePr = SteppingHelixPropagator(magneticField_, alongMomentum);
0746   SteppingHelixPropagator oppositeStHePr = SteppingHelixPropagator(magneticField_, oppositeToMomentum);
0747 
0748   //double tolerance = 5.e-5;
0749   defaultRKPropagator::Product aprod(magneticField_, alongMomentum, 5.e-5);
0750   auto& alongRKPr = aprod.propagator;
0751   defaultRKPropagator::Product oprod(magneticField_, oppositeToMomentum, 5.e-5);
0752   auto& oppositeRKPr = oprod.propagator;
0753 
0754   float epsilon = 5.;
0755   SmartPropagator alongSmPr = SmartPropagator(alongRKPr, alongStHePr, magneticField_, alongMomentum, epsilon);
0756   SmartPropagator oppositeSmPr =
0757       SmartPropagator(oppositeRKPr, oppositeStHePr, magneticField_, oppositeToMomentum, epsilon);
0758 
0759   // ................................................ selected/global muon
0760   //itMuon -->  Jim's globalMuon
0761   for (MuonCollection::const_iterator itMuon = smuons->begin(); itMuon != smuons->end(); ++itMuon) {
0762     if (debug_) {
0763       edm::LogVerbatim("GlobalTrackerMuonAlignment")
0764           << " mu gM is GM Mu SaM tM " << itMuon->isGlobalMuon() << " " << itMuon->isMuon() << " "
0765           << itMuon->isStandAloneMuon() << " " << itMuon->isTrackerMuon() << " ";
0766     }
0767 
0768     // get information about the innermost muon hit -------------------------
0769     TransientTrack muTT(itMuon->outerTrack(), magneticField_, trackingGeometry_);
0770     TrajectoryStateOnSurface innerMuTSOS = muTT.innermostMeasurementState();
0771     TrajectoryStateOnSurface outerMuTSOS = muTT.outermostMeasurementState();
0772 
0773     // get information about the outermost tracker hit -----------------------
0774     TransientTrack trackTT(itMuon->track(), magneticField_, trackingGeometry_);
0775     TrajectoryStateOnSurface outerTrackTSOS = trackTT.outermostMeasurementState();
0776     TrajectoryStateOnSurface innerTrackTSOS = trackTT.innermostMeasurementState();
0777 
0778     GlobalPoint pointTrackIn = innerTrackTSOS.globalPosition();
0779     GlobalPoint pointTrackOut = outerTrackTSOS.globalPosition();
0780     float lenghtTrack = (pointTrackOut - pointTrackIn).mag();
0781     GlobalPoint pointMuonIn = innerMuTSOS.globalPosition();
0782     GlobalPoint pointMuonOut = outerMuTSOS.globalPosition();
0783     float lenghtMuon = (pointMuonOut - pointMuonIn).mag();
0784     GlobalVector momentumTrackOut = outerTrackTSOS.globalMomentum();
0785     GlobalVector momentumTrackIn = innerTrackTSOS.globalMomentum();
0786     float distanceInIn = (pointTrackIn - pointMuonIn).mag();
0787     float distanceInOut = (pointTrackIn - pointMuonOut).mag();
0788     float distanceOutIn = (pointTrackOut - pointMuonIn).mag();
0789     float distanceOutOut = (pointTrackOut - pointMuonOut).mag();
0790 
0791     if (debug_) {
0792       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " pointTrackIn " << pointTrackIn;
0793       edm::LogVerbatim("GlobalTrackerMuonAlignment") << "          Out " << pointTrackOut << " lenght " << lenghtTrack;
0794       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " point MuonIn " << pointMuonIn;
0795       edm::LogVerbatim("GlobalTrackerMuonAlignment") << "          Out " << pointMuonOut << " lenght " << lenghtMuon;
0796       edm::LogVerbatim("GlobalTrackerMuonAlignment")
0797           << " momeTrackIn Out " << momentumTrackIn << " " << momentumTrackOut;
0798       edm::LogVerbatim("GlobalTrackerMuonAlignment")
0799           << " doi io ii oo " << distanceOutIn << " " << distanceInOut << " " << distanceInIn << " " << distanceOutOut;
0800     }
0801 
0802     if (lenghtTrack < 90.)
0803       continue;
0804     if (lenghtMuon < 300.)
0805       continue;
0806     if (momentumTrackIn.mag() < 15. || momentumTrackOut.mag() < 15.)
0807       continue;
0808     if (trackTT.charge() != muTT.charge())
0809       continue;
0810 
0811     if (debug_)
0812       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " passed lenght momentum cuts";
0813 
0814     GlobalVector GRm, GPm, Nl, Rm, Pm, Rt, Pt, Rt0;
0815     AlgebraicSymMatrix66 Cm, C0, Ce, C1;
0816 
0817     TrajectoryStateOnSurface extrapolationT;
0818     int tsosMuonIf = 0;
0819 
0820     if (isolatedMuonMode_) {  //------------------------------- Isolated Muon -----
0821       const Surface& refSurface = innerMuTSOS.surface();
0822       ConstReferenceCountingPointer<TangentPlane> tpMuLocal(refSurface.tangentPlane(innerMuTSOS.localPosition()));
0823       Nl = tpMuLocal->normalVector();
0824       ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
0825       GlobalVector Ng = tpMuGlobal->normalVector();
0826       Surface* surf = (Surface*)&refSurface;
0827       const Plane* refPlane = dynamic_cast<Plane*>(surf);
0828       GlobalVector Nlp = refPlane->normalVector();
0829 
0830       if (debug_) {
0831         edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
0832                                                        << " alfa " << int(asin(Nl.x()) * 57.29578);
0833         edm::LogVerbatim("GlobalTrackerMuonAlignment") << "  global " << Ng.x() << " " << Ng.y() << " " << Ng.z();
0834         edm::LogVerbatim("GlobalTrackerMuonAlignment") << "      lp " << Nlp.x() << " " << Nlp.y() << " " << Nlp.z();
0835         edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug")
0836             << " Nlocal Nglobal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " " << Ng.x() << " " << Ng.y() << " "
0837             << Ng.z() << " alfa " << static_cast<int>(asin(Nl.x()) * 57.29578);
0838       }
0839 
0840       //                          extrapolation to innermost muon hit
0841       extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
0842       //extrapolationT = propagator->propagate(outerTrackTSOS, refSurface);
0843 
0844       if (!extrapolationT.isValid()) {
0845         if (false & alarm)
0846           edm::LogVerbatim("GlobalTrackerMuonAlignment") << " !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
0847               //<<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
0848               ;
0849         continue;
0850       }
0851       tsosMuonIf = 1;
0852       Rt = GlobalVector((extrapolationT.globalPosition()).x(),
0853                         (extrapolationT.globalPosition()).y(),
0854                         (extrapolationT.globalPosition()).z());
0855 
0856       Pt = extrapolationT.globalMomentum();
0857       //                          global parameters of muon
0858       GRm = GlobalVector(
0859           (innerMuTSOS.globalPosition()).x(), (innerMuTSOS.globalPosition()).y(), (innerMuTSOS.globalPosition()).z());
0860       GPm = innerMuTSOS.globalMomentum();
0861       Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
0862                          (outerTrackTSOS.globalPosition()).y(),
0863                          (outerTrackTSOS.globalPosition()).z());
0864       Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + innerMuTSOS.cartesianError().matrix());
0865       C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
0866       Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
0867       C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
0868 
0869     }  //                    -------------------------------   end Isolated Muon -----
0870 
0871     if (cosmicMuonMode_) {  //------------------------------- Cosmic Muon -----
0872 
0873       if ((distanceOutIn <= distanceInOut) & (distanceOutIn <= distanceInIn) &
0874           (distanceOutIn <= distanceOutOut)) {  // ----- Out - In ------
0875         if (debug_)
0876           edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -----    Out - In -----";
0877 
0878         const Surface& refSurface = innerMuTSOS.surface();
0879         ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
0880         Nl = tpMuGlobal->normalVector();
0881 
0882         //                          extrapolation to innermost muon hit
0883         extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
0884         //extrapolationT = propagator->propagate(outerTrackTSOS, refSurface);
0885 
0886         if (!extrapolationT.isValid()) {
0887           if (false & alarm)
0888             edm::LogVerbatim("GlobalTrackerMuonAlignment") << " !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
0889                 //<<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
0890                 ;
0891           continue;
0892         }
0893         if (debug_)
0894           edm::LogVerbatim("GlobalTrackerMuonAlignment") << " extrapolationT.isValid " << extrapolationT.isValid();
0895 
0896         tsosMuonIf = 1;
0897         Rt = GlobalVector((extrapolationT.globalPosition()).x(),
0898                           (extrapolationT.globalPosition()).y(),
0899                           (extrapolationT.globalPosition()).z());
0900 
0901         Pt = extrapolationT.globalMomentum();
0902         //                          global parameters of muon
0903         GRm = GlobalVector(
0904             (innerMuTSOS.globalPosition()).x(), (innerMuTSOS.globalPosition()).y(), (innerMuTSOS.globalPosition()).z());
0905         GPm = innerMuTSOS.globalMomentum();
0906         Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
0907                            (outerTrackTSOS.globalPosition()).y(),
0908                            (outerTrackTSOS.globalPosition()).z());
0909         Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + innerMuTSOS.cartesianError().matrix());
0910         C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
0911         Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
0912         C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
0913 
0914         if (false & debug_) {
0915           edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ->propDir " << propagator->propagationDirection();
0916           PropagationDirectionChooser Chooser;
0917           edm::LogVerbatim("GlobalTrackerMuonAlignment")
0918               << " propDirCh " << Chooser.operator()(*outerTrackTSOS.freeState(), refSurface) << " Ch == along "
0919               << (alongMomentum == Chooser.operator()(*outerTrackTSOS.freeState(), refSurface));
0920           edm::LogVerbatim("GlobalTrackerMuonAlignment")
0921               << " --- Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
0922               << " alfa " << int(asin(Nl.x()) * 57.29578);
0923         }
0924       }  //                                       enf of ---- Out - In -----
0925 
0926       else if ((distanceInOut <= distanceInIn) & (distanceInOut <= distanceOutIn) &
0927                (distanceInOut <= distanceOutOut)) {  // ----- In - Out ------
0928         if (debug_)
0929           edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -----    In - Out -----";
0930 
0931         const Surface& refSurface = outerMuTSOS.surface();
0932         ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
0933         Nl = tpMuGlobal->normalVector();
0934 
0935         //                          extrapolation to outermost muon hit
0936         extrapolationT = oppositeSmPr.propagate(innerTrackTSOS, refSurface);
0937 
0938         if (!extrapolationT.isValid()) {
0939           if (false & alarm)
0940             edm::LogVerbatim("GlobalTrackerMuonAlignment") << " !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
0941                                                            << "\a\a\a\a\a\a\a\a" << extrapolationT.isValid();
0942           continue;
0943         }
0944         if (debug_)
0945           edm::LogVerbatim("GlobalTrackerMuonAlignment") << " extrapolationT.isValid " << extrapolationT.isValid();
0946 
0947         tsosMuonIf = 2;
0948         Rt = GlobalVector((extrapolationT.globalPosition()).x(),
0949                           (extrapolationT.globalPosition()).y(),
0950                           (extrapolationT.globalPosition()).z());
0951 
0952         Pt = extrapolationT.globalMomentum();
0953         //                          global parameters of muon
0954         GRm = GlobalVector(
0955             (outerMuTSOS.globalPosition()).x(), (outerMuTSOS.globalPosition()).y(), (outerMuTSOS.globalPosition()).z());
0956         GPm = outerMuTSOS.globalMomentum();
0957         Rt0 = GlobalVector((innerTrackTSOS.globalPosition()).x(),
0958                            (innerTrackTSOS.globalPosition()).y(),
0959                            (innerTrackTSOS.globalPosition()).z());
0960         Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + outerMuTSOS.cartesianError().matrix());
0961         C0 = AlgebraicSymMatrix66(innerTrackTSOS.cartesianError().matrix());
0962         Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
0963         C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
0964 
0965         if (false & debug_) {
0966           edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ->propDir " << propagator->propagationDirection();
0967           PropagationDirectionChooser Chooser;
0968           edm::LogVerbatim("GlobalTrackerMuonAlignment")
0969               << " propDirCh " << Chooser.operator()(*innerTrackTSOS.freeState(), refSurface) << " Ch == oppisite "
0970               << (oppositeToMomentum == Chooser.operator()(*innerTrackTSOS.freeState(), refSurface));
0971           edm::LogVerbatim("GlobalTrackerMuonAlignment")
0972               << " --- Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
0973               << " alfa " << int(asin(Nl.x()) * 57.29578);
0974         }
0975       }  //                                        enf of ---- In - Out -----
0976 
0977       else if ((distanceOutOut <= distanceInOut) & (distanceOutOut <= distanceInIn) &
0978                (distanceOutOut <= distanceOutIn)) {  // ----- Out - Out ------
0979         if (debug_)
0980           edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -----    Out - Out -----";
0981 
0982         // reject: momentum of track has opposite direction to muon track
0983         continue;
0984 
0985         const Surface& refSurface = outerMuTSOS.surface();
0986         ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
0987         Nl = tpMuGlobal->normalVector();
0988 
0989         //                          extrapolation to outermost muon hit
0990         extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
0991 
0992         if (!extrapolationT.isValid()) {
0993           if (alarm)
0994             edm::LogVerbatim("GlobalTrackerMuonAlignment") << " !!!!!!!!! Catastrophe Out-Out extrapolationT.isValid "
0995                                                            << "\a\a\a\a\a\a\a\a" << extrapolationT.isValid();
0996           continue;
0997         }
0998         if (debug_)
0999           edm::LogVerbatim("GlobalTrackerMuonAlignment") << " extrapolationT.isValid " << extrapolationT.isValid();
1000 
1001         tsosMuonIf = 2;
1002         Rt = GlobalVector((extrapolationT.globalPosition()).x(),
1003                           (extrapolationT.globalPosition()).y(),
1004                           (extrapolationT.globalPosition()).z());
1005 
1006         Pt = extrapolationT.globalMomentum();
1007         //                          global parameters of muon
1008         GRm = GlobalVector(
1009             (outerMuTSOS.globalPosition()).x(), (outerMuTSOS.globalPosition()).y(), (outerMuTSOS.globalPosition()).z());
1010         GPm = outerMuTSOS.globalMomentum();
1011         Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
1012                            (outerTrackTSOS.globalPosition()).y(),
1013                            (outerTrackTSOS.globalPosition()).z());
1014         Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + outerMuTSOS.cartesianError().matrix());
1015         C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
1016         Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
1017         C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
1018 
1019         if (debug_) {
1020           PropagationDirectionChooser Chooser;
1021           edm::LogVerbatim("GlobalTrackerMuonAlignment")
1022               << " propDirCh " << Chooser.operator()(*outerTrackTSOS.freeState(), refSurface) << " Ch == along "
1023               << (alongMomentum == Chooser.operator()(*outerTrackTSOS.freeState(), refSurface));
1024           edm::LogVerbatim("GlobalTrackerMuonAlignment")
1025               << " --- Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1026               << " alfa " << int(asin(Nl.x()) * 57.29578);
1027           edm::LogVerbatim("GlobalTrackerMuonAlignment")
1028               << "     Nornal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1029               << " alfa " << int(asin(Nl.x()) * 57.29578);
1030         }
1031       }  //                                       enf of ---- Out - Out -----
1032       else {
1033         if (alarm)
1034           edm::LogVerbatim("GlobalTrackerMuonAlignment")
1035               << " ------- !!!!!!!!!!!!!!! No proper Out - In \a\a\a\a\a\a\a";
1036         continue;
1037       }
1038 
1039     }  //                     -------------------------------   end Cosmic Muon -----
1040 
1041     if (tsosMuonIf == 0) {
1042       if (info) {
1043         edm::LogVerbatim("GlobalTrackerMuonAlignment") << "No tsosMuon !!!!!!";
1044         continue;
1045       }
1046     }
1047     TrajectoryStateOnSurface tsosMuon;
1048     if (tsosMuonIf == 1)
1049       tsosMuon = muTT.innermostMeasurementState();
1050     else
1051       tsosMuon = muTT.outermostMeasurementState();
1052 
1053     //GlobalTrackerMuonAlignment::misalignMuon(GRm, GPm, Nl, Rt, Rm, Pm);
1054     AlgebraicVector4 LPRm;  // muon local (dx/dz, dy/dz, x, y)
1055     GlobalTrackerMuonAlignment::misalignMuonL(GRm, GPm, Nl, Rt, Rm, Pm, LPRm, extrapolationT, tsosMuon);
1056 
1057     GlobalVector resR = Rm - Rt;
1058     GlobalVector resP0 = Pm - Pt;
1059     GlobalVector resP = Pm / Pm.mag() - Pt / Pt.mag();
1060     float RelMomResidual = (Pm.mag() - Pt.mag()) / (Pt.mag() + 1.e-6);
1061     ;
1062 
1063     AlgebraicVector6 Vm;
1064     Vm(0) = resR.x();
1065     Vm(1) = resR.y();
1066     Vm(2) = resR.z();
1067     Vm(3) = resP0.x();
1068     Vm(4) = resP0.y();
1069     Vm(5) = resP0.z();
1070     float Rmuon = Rm.perp();
1071     float Zmuon = Rm.z();
1072     float alfa_x = atan2(Nl.x(), Nl.y()) * 57.29578;
1073 
1074     if (debug_) {
1075       edm::LogVerbatim("GlobalTrackerMuonAlignment")
1076           << " Nx Ny Nz alfa_x " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " " << alfa_x;
1077       edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug")
1078           << " Rm " << Rm << "\n Rt " << Rt << "\n resR " << resR << "\n resP " << resP << " dp/p " << RelMomResidual;
1079     }
1080 
1081     double chi_d = 0;
1082     for (int i = 0; i <= 5; i++)
1083       chi_d += Vm(i) * Vm(i) / Cm(i, i);
1084 
1085     AlgebraicVector5 Vml(tsosMuon.localParameters().vector() - extrapolationT.localParameters().vector());
1086     AlgebraicSymMatrix55 m(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
1087     AlgebraicSymMatrix55 Cml(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
1088     bool ierrLoc = !m.Invert();
1089     if (ierrLoc && debug_ && info) {
1090       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ==== Error inverting Local covariance matrix ==== ";
1091       continue;
1092     }
1093     double chi_Loc = ROOT::Math::Similarity(Vml, m);
1094     if (debug_)
1095       edm::LogVerbatim("GlobalTrackerMuonAlignment")
1096           << " chi_Loc px/pz/err " << chi_Loc << " " << Vml(1) / std::sqrt(Cml(1, 1));
1097 
1098     if (Pt.mag() < 15.)
1099       continue;
1100     if (Pm.mag() < 5.)
1101       continue;
1102 
1103     //if(Pt.mag() < 30.) continue;          // momenum cut  < 30GeV
1104     //if(Pt.mag() < 60.) continue;          // momenum cut  < 30GeV
1105     //if(Pt.mag() > 50.) continue;          //  momenum cut > 50GeV
1106     //if(Pt.mag() > 100.) continue;          //  momenum cut > 100GeV
1107     //if(trackTT.charge() < 0) continue;    // select positive charge
1108     //if(trackTT.charge() > 0) continue;    // select negative charge
1109 
1110     //if(fabs(resR.x()) > 5.) continue;     // strong cut  X
1111     //if(fabs(resR.y()) > 5.) continue;     //             Y
1112     //if(fabs(resR.z()) > 5.) continue;     //             Z
1113     //if(fabs(resR.mag()) > 7.5) continue;   //            dR
1114 
1115     /*
1116     //if(fabs(RelMomResidual) > 0.5) continue;
1117     if(fabs(resR.x()) > 20.) continue;
1118     if(fabs(resR.y()) > 20.) continue;
1119     if(fabs(resR.z()) > 20.) continue;
1120     if(fabs(resR.mag()) > 30.) continue;
1121     if(fabs(resP.x()) > 0.06) continue;
1122     if(fabs(resP.y()) > 0.06) continue;
1123     if(fabs(resP.z()) > 0.06) continue;
1124     if(chi_d > 40.) continue;
1125     */
1126 
1127     //                                            select Barrel
1128     //if(Rmuon < 400. || Rmuon > 450.) continue;
1129     //if(Zmuon < -600. || Zmuon > 600.) continue;
1130     //if(fabs(Nl.z()) > 0.95) continue;
1131     //MuSelect = " Barrel";
1132     //                                                  EndCap1
1133     //if(Rmuon < 120. || Rmuon > 450.) continue;
1134     //if(Zmuon < -720.) continue;
1135     //if(Zmuon > -580.) continue;
1136     //if(fabs(Nl.z()) < 0.95) continue;
1137     //MuSelect = " EndCap1";
1138     //                                                  EndCap2
1139     //if(Rmuon < 120. || Rmuon > 450.) continue;
1140     //if(Zmuon >  720.) continue;
1141     //if(Zmuon <  580.) continue;
1142     //if(fabs(Nl.z()) < 0.95) continue;
1143     //MuSelect = " EndCap2";
1144     //                                                 select All
1145     if (Rmuon < 120. || Rmuon > 450.)
1146       continue;
1147     if (Zmuon < -720. || Zmuon > 720.)
1148       continue;
1149     MuSelect = " Barrel+EndCaps";
1150 
1151     if (debug_)
1152       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " .............. passed all cuts";
1153 
1154     N_track++;
1155     //                     gradient and Hessian for each track
1156 
1157     GlobalTrackerMuonAlignment::gradientGlobalAlg(Rt, Pt, Rm, Nl, Cm);
1158     GlobalTrackerMuonAlignment::gradientGlobal(Rt, Pt, Rm, Pm, Nl, Cm);
1159 
1160     CLHEP::HepSymMatrix covLoc(4, 0);
1161     for (int i = 1; i <= 4; i++)
1162       for (int j = 1; j <= i; j++) {
1163         covLoc(i, j) = (tsosMuon.localError().matrix() + extrapolationT.localError().matrix())(i, j);
1164         //if(i != j) Cov(i,j) = 0.;
1165       }
1166 
1167     const Surface& refSurface = tsosMuon.surface();
1168     CLHEP::HepMatrix rotLoc(3, 3, 0);
1169     rotLoc(1, 1) = refSurface.rotation().xx();
1170     rotLoc(1, 2) = refSurface.rotation().xy();
1171     rotLoc(1, 3) = refSurface.rotation().xz();
1172 
1173     rotLoc(2, 1) = refSurface.rotation().yx();
1174     rotLoc(2, 2) = refSurface.rotation().yy();
1175     rotLoc(2, 3) = refSurface.rotation().yz();
1176 
1177     rotLoc(3, 1) = refSurface.rotation().zx();
1178     rotLoc(3, 2) = refSurface.rotation().zy();
1179     rotLoc(3, 3) = refSurface.rotation().zz();
1180 
1181     CLHEP::HepVector posLoc(3);
1182     posLoc(1) = refSurface.position().x();
1183     posLoc(2) = refSurface.position().y();
1184     posLoc(3) = refSurface.position().z();
1185 
1186     GlobalTrackerMuonAlignment::gradientLocal(Rt, Pt, Rm, Pm, Nl, covLoc, rotLoc, posLoc, LPRm);
1187 
1188     if (debug_) {
1189       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Norm   " << Nl;
1190       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " posLoc " << posLoc.T();
1191       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " rotLoc " << rotLoc;
1192     }
1193 
1194     // -----------------------------------------------------  fill histogram
1195     histo->Fill(itMuon->track()->pt());
1196 
1197     //histo2->Fill(itMuon->track()->outerP());
1198     histo2->Fill(Pt.mag());
1199     histo3->Fill((PI / 2. - itMuon->track()->outerTheta()));
1200     histo4->Fill(itMuon->track()->phi());
1201     histo5->Fill(Rmuon);
1202     histo6->Fill(Zmuon);
1203     histo7->Fill(RelMomResidual);
1204     //histo8->Fill(chi);
1205     histo8->Fill(chi_d);
1206 
1207     histo101->Fill(Zmuon, Rmuon);
1208     histo101->Fill(Rt0.z(), Rt0.perp());
1209     histo102->Fill(Rt0.x(), Rt0.y());
1210     histo102->Fill(Rm.x(), Rm.y());
1211 
1212     histo11->Fill(resR.mag());
1213     if (fabs(Nl.x()) < 0.98)
1214       histo12->Fill(resR.x());
1215     if (fabs(Nl.y()) < 0.98)
1216       histo13->Fill(resR.y());
1217     if (fabs(Nl.z()) < 0.98)
1218       histo14->Fill(resR.z());
1219     histo15->Fill(resP.x());
1220     histo16->Fill(resP.y());
1221     histo17->Fill(resP.z());
1222 
1223     if ((fabs(Nl.x()) < 0.98) && (fabs(Nl.y()) < 0.98) && (fabs(Nl.z()) < 0.98)) {
1224       histo18->Fill(std::sqrt(C0(0, 0)));
1225       histo19->Fill(std::sqrt(C1(0, 0)));
1226       histo20->Fill(std::sqrt(C1(0, 0) + Ce(0, 0)));
1227     }
1228     if (fabs(Nl.x()) < 0.98)
1229       histo21->Fill(Vm(0) / std::sqrt(Cm(0, 0)));
1230     if (fabs(Nl.y()) < 0.98)
1231       histo22->Fill(Vm(1) / std::sqrt(Cm(1, 1)));
1232     if (fabs(Nl.z()) < 0.98)
1233       histo23->Fill(Vm(2) / std::sqrt(Cm(2, 2)));
1234     histo24->Fill(Vm(3) / std::sqrt(C1(3, 3) + Ce(3, 3)));
1235     histo25->Fill(Vm(4) / std::sqrt(C1(4, 4) + Ce(4, 4)));
1236     histo26->Fill(Vm(5) / std::sqrt(C1(5, 5) + Ce(5, 5)));
1237     histo27->Fill(Nl.x());
1238     histo28->Fill(Nl.y());
1239     histo29->Fill(lenghtTrack);
1240     histo30->Fill(lenghtMuon);
1241     histo31->Fill(chi_Loc);
1242     histo32->Fill(Vml(1) / std::sqrt(Cml(1, 1)));
1243     histo33->Fill(Vml(2) / std::sqrt(Cml(2, 2)));
1244     histo34->Fill(Vml(3) / std::sqrt(Cml(3, 3)));
1245     histo35->Fill(Vml(4) / std::sqrt(Cml(4, 4)));
1246 
1247     if (debug_) {  //--------------------------------- debug print ----------
1248 
1249       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " diag 0[ " << C0(0, 0) << " " << C0(1, 1) << " " << C0(2, 2)
1250                                                      << " " << C0(3, 3) << " " << C0(4, 4) << " " << C0(5, 5) << " ]";
1251       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " diag e[ " << Ce(0, 0) << " " << Ce(1, 1) << " " << Ce(2, 2)
1252                                                      << " " << Ce(3, 3) << " " << Ce(4, 4) << " " << Ce(5, 5) << " ]";
1253       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " diag 1[ " << C1(0, 0) << " " << C1(1, 1) << " " << C1(2, 2)
1254                                                      << " " << C1(3, 3) << " " << C1(4, 4) << " " << C1(5, 5) << " ]";
1255       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Rm   " << Rm.x() << " " << Rm.y() << " " << Rm.z() << " Pm   "
1256                                                      << Pm.x() << " " << Pm.y() << " " << Pm.z();
1257       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Rt   " << Rt.x() << " " << Rt.y() << " " << Rt.z() << " Pt   "
1258                                                      << Pt.x() << " " << Pt.y() << " " << Pt.z();
1259       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Nl*(Rm-Rt) " << Nl.dot(Rm - Rt);
1260       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " resR " << resR.x() << " " << resR.y() << " " << resR.z()
1261                                                      << " resP " << resP.x() << " " << resP.y() << " " << resP.z();
1262       edm::LogVerbatim("GlobalTrackerMuonAlignment")
1263           << " Rm-t " << (Rm - Rt).x() << " " << (Rm - Rt).y() << " " << (Rm - Rt).z() << " Pm-t " << (Pm - Pt).x()
1264           << " " << (Pm - Pt).y() << " " << (Pm - Pt).z();
1265       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Vm   " << Vm;
1266       edm::LogVerbatim("GlobalTrackerMuonAlignment")
1267           << " +-   " << std::sqrt(Cm(0, 0)) << " " << std::sqrt(Cm(1, 1)) << " " << std::sqrt(Cm(2, 2)) << " "
1268           << std::sqrt(Cm(3, 3)) << " " << std::sqrt(Cm(4, 4)) << " " << std::sqrt(Cm(5, 5));
1269       edm::LogVerbatim("GlobalTrackerMuonAlignment")
1270           << " Pmuon Ptrack dP/Ptrack " << itMuon->outerTrack()->p() << " " << itMuon->track()->outerP() << " "
1271           << (itMuon->outerTrack()->p() - itMuon->track()->outerP()) / itMuon->track()->outerP();
1272       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " cov matrix ";
1273       edm::LogVerbatim("GlobalTrackerMuonAlignment") << Cm;
1274       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " diag [ " << Cm(0, 0) << " " << Cm(1, 1) << " " << Cm(2, 2)
1275                                                      << " " << Cm(3, 3) << " " << Cm(4, 4) << " " << Cm(5, 5) << " ]";
1276 
1277       AlgebraicSymMatrix66 Ro;
1278       double Diag[6];
1279       for (int i = 0; i <= 5; i++)
1280         Diag[i] = std::sqrt(Cm(i, i));
1281       for (int i = 0; i <= 5; i++)
1282         for (int j = 0; j <= 5; j++)
1283           Ro(i, j) = Cm(i, j) / Diag[i] / Diag[j];
1284       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " correlation matrix ";
1285       edm::LogVerbatim("GlobalTrackerMuonAlignment") << Ro;
1286 
1287       AlgebraicSymMatrix66 CmI;
1288       for (int i = 0; i <= 5; i++)
1289         for (int j = 0; j <= 5; j++)
1290           CmI(i, j) = Cm(i, j);
1291 
1292       bool ierr = !CmI.Invert();
1293       if (ierr) {
1294         if (alarm || debug_)
1295           edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Error inverse covariance matrix !!!!!!!!!!!";
1296         continue;
1297       }
1298       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " inverse cov matrix ";
1299       edm::LogVerbatim("GlobalTrackerMuonAlignment") << Cm;
1300 
1301       double chi = ROOT::Math::Similarity(Vm, CmI);
1302       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " chi chi_d " << chi << " " << chi_d;
1303     }  // end of debug_ printout  --------------------------------------------
1304 
1305   }  // end loop on selected muons, i.e. Jim's globalMuon
1306 
1307 }  //end of analyzeTrackTrack
1308 
1309 // ------- method to analyze recoTrack & trajectoryMuon of Global Muon ------
1310 void GlobalTrackerMuonAlignment::analyzeTrackTrajectory(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
1311   using namespace edm;
1312   using namespace reco;
1313   //debug_ = true;
1314   bool info = false;
1315   bool alarm = false;
1316   //bool alarm = true;
1317   double PI = 3.1415927;
1318 
1319   //-*-*-*-*-*-*-
1320   cosmicMuonMode_ = true;  // true: both Cosmic and IsolatedMuon are treated with 1,2 tracks
1321   //cosmicMuonMode_ = false; // true: both Cosmic and IsolatedMuon are treated with 1,2 tracks
1322   //-*-*-*-*-*-*-
1323   isolatedMuonMode_ = !cosmicMuonMode_;  //true: only IsolatedMuon are treated with 1 track
1324 
1325   N_event++;
1326   if (info || debug_) {
1327     edm::LogVerbatim("GlobalTrackerMuonAlignment") << "----- event " << N_event << " -- tracks " << N_track << " ---";
1328     if (cosmicMuonMode_)
1329       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " treated as CosmicMu ";
1330     if (isolatedMuonMode_)
1331       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " treated as IsolatedMu ";
1332     edm::LogVerbatim("GlobalTrackerMuonAlignment");
1333   }
1334 
1335   const edm::Handle<reco::TrackCollection>& tracks =
1336       ((collectionIsolated)
1337            ? iEvent.getHandle(trackIsolateToken_)
1338            : ((collectionCosmic) ? iEvent.getHandle(trackCosmicToken_) : iEvent.getHandle(trackToken_)));
1339   const edm::Handle<reco::TrackCollection>& muons =
1340       ((collectionIsolated) ? iEvent.getHandle(muonIsolateToken_)
1341                             : ((collectionCosmic) ? iEvent.getHandle(muonCosmicToken_) : iEvent.getHandle(muonToken_)));
1342   const edm::Handle<reco::TrackCollection>& gmuons =
1343       ((collectionIsolated)
1344            ? iEvent.getHandle(gmuonIsolateToken_)
1345            : ((collectionCosmic) ? iEvent.getHandle(gmuonCosmicToken_) : iEvent.getHandle(gmuonToken_)));
1346   const edm::Handle<reco::MuonCollection>& smuons =
1347       ((collectionIsolated)
1348            ? iEvent.getHandle(smuonIsolateToken_)
1349            : ((collectionCosmic) ? iEvent.getHandle(smuonCosmicToken_) : iEvent.getHandle(smuonToken_)));
1350 
1351   if (debug_) {
1352     edm::LogVerbatim("GlobalTrackerMuonAlignment")
1353         << " ievBunch " << iEvent.bunchCrossing() << " runN " << (int)iEvent.run();
1354     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " N tracks s/amu gmu selmu " << tracks->size() << " "
1355                                                    << muons->size() << " " << gmuons->size() << " " << smuons->size();
1356     for (MuonCollection::const_iterator itMuon = smuons->begin(); itMuon != smuons->end(); ++itMuon) {
1357       edm::LogVerbatim("GlobalTrackerMuonAlignment")
1358           << " is isolatValid Matches " << itMuon->isIsolationValid() << " " << itMuon->isMatchesValid();
1359     }
1360   }
1361 
1362   if (isolatedMuonMode_) {  // ---- Only 1 Isolated Muon --------
1363     if (tracks->size() != 1)
1364       return;
1365     if (muons->size() != 1)
1366       return;
1367     if (gmuons->size() != 1)
1368       return;
1369     if (smuons->size() != 1)
1370       return;
1371   }
1372 
1373   if (cosmicMuonMode_) {  // ---- 1,2 Cosmic Muon --------
1374     if (smuons->size() > 2)
1375       return;
1376     if (tracks->size() != smuons->size())
1377       return;
1378     if (muons->size() != smuons->size())
1379       return;
1380     if (gmuons->size() != smuons->size())
1381       return;
1382   }
1383 
1384   // ok mc_isolated_mu
1385   //edm::ESHandle<TrackerGeometry> trackerGeometry;
1386   //iSetup.get<TrackerDigiGeometryRecord>().get(trackerGeometry);
1387   // ok mc_isolated_mu
1388   //edm::ESHandle<DTGeometry> dtGeometry;
1389   //iSetup.get<MuonGeometryRecord>().get(dtGeometry);
1390   // don't work
1391   //edm::ESHandle<CSCGeometry> cscGeometry;
1392   //iSetup.get<MuonGeometryRecord>().get(cscGeometry);
1393 
1394   if (watchTrackingGeometry_.check(iSetup) || !trackingGeometry_) {
1395     edm::ESHandle<GlobalTrackingGeometry> trackingGeometry = iSetup.getHandle(m_TkGeometryToken);
1396     trackingGeometry_ = &*trackingGeometry;
1397     theTrackingGeometry = trackingGeometry;
1398   }
1399 
1400   if (watchMagneticFieldRecord_.check(iSetup) || !magneticField_) {
1401     edm::ESHandle<MagneticField> magneticField = iSetup.getHandle(m_MagFieldToken);
1402     magneticField_ = &*magneticField;
1403   }
1404 
1405   if (watchGlobalPositionRcd_.check(iSetup) || !globalPositionRcd_) {
1406     edm::ESHandle<Alignments> globalPositionRcd = iSetup.getHandle(m_globalPosToken);
1407     globalPositionRcd_ = &*globalPositionRcd;
1408     for (std::vector<AlignTransform>::const_iterator i = globalPositionRcd_->m_align.begin();
1409          i != globalPositionRcd_->m_align.end();
1410          ++i) {
1411       if (DetId(DetId::Tracker).rawId() == i->rawId())
1412         iteratorTrackerRcd = i;
1413       if (DetId(DetId::Muon).rawId() == i->rawId())
1414         iteratorMuonRcd = i;
1415       if (DetId(DetId::Ecal).rawId() == i->rawId())
1416         iteratorEcalRcd = i;
1417       if (DetId(DetId::Hcal).rawId() == i->rawId())
1418         iteratorHcalRcd = i;
1419     }
1420     if (debug_) {
1421       edm::LogVerbatim("GlobalTrackerMuonAlignment")
1422           << "=== iteratorTrackerRcd " << iteratorTrackerRcd->rawId() << " ====\n"
1423           << " translation " << iteratorTrackerRcd->translation() << "\n"
1424           << " angles " << iteratorTrackerRcd->rotation().eulerAngles();
1425       edm::LogVerbatim("GlobalTrackerMuonAlignment") << iteratorTrackerRcd->rotation();
1426       edm::LogVerbatim("GlobalTrackerMuonAlignment") << "=== iteratorMuonRcd " << iteratorMuonRcd->rawId() << " ====\n"
1427                                                      << " translation " << iteratorMuonRcd->translation() << "\n"
1428                                                      << " angles " << iteratorMuonRcd->rotation().eulerAngles();
1429       edm::LogVerbatim("GlobalTrackerMuonAlignment") << iteratorMuonRcd->rotation();
1430     }
1431   }  // end of GlobalPositionRcd
1432 
1433   edm::ESHandle<Propagator> propagator = iSetup.getHandle(m_propToken);
1434 
1435   SteppingHelixPropagator alongStHePr = SteppingHelixPropagator(magneticField_, alongMomentum);
1436   SteppingHelixPropagator oppositeStHePr = SteppingHelixPropagator(magneticField_, oppositeToMomentum);
1437 
1438   //double tolerance = 5.e-5;
1439   defaultRKPropagator::Product aprod(magneticField_, alongMomentum, 5.e-5);
1440   auto& alongRKPr = aprod.propagator;
1441   defaultRKPropagator::Product oprod(magneticField_, oppositeToMomentum, 5.e-5);
1442   auto& oppositeRKPr = oprod.propagator;
1443 
1444   float epsilon = 5.;
1445   SmartPropagator alongSmPr = SmartPropagator(alongRKPr, alongStHePr, magneticField_, alongMomentum, epsilon);
1446   SmartPropagator oppositeSmPr =
1447       SmartPropagator(oppositeRKPr, oppositeStHePr, magneticField_, oppositeToMomentum, epsilon);
1448 
1449   if (defineFitter) {
1450     if (debug_)
1451       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ............... DEFINE FITTER ...................";
1452     KFUpdator* theUpdator = new KFUpdator();
1453     //Chi2MeasurementEstimator* theEstimator = new Chi2MeasurementEstimator(30);
1454     Chi2MeasurementEstimator* theEstimator = new Chi2MeasurementEstimator(100000, 100000);
1455     theFitter = new KFTrajectoryFitter(alongSmPr, *theUpdator, *theEstimator);
1456     theSmoother = new KFTrajectorySmoother(alongSmPr, *theUpdator, *theEstimator);
1457     theFitterOp = new KFTrajectoryFitter(oppositeSmPr, *theUpdator, *theEstimator);
1458     theSmootherOp = new KFTrajectorySmoother(oppositeSmPr, *theUpdator, *theEstimator);
1459 
1460     edm::ESHandle<TransientTrackingRecHitBuilder> builder = iSetup.getHandle(m_ttrhBuilderToken);
1461     this->TTRHBuilder = &(*builder);
1462     this->MuRHBuilder = new MuonTransientTrackingRecHitBuilder(theTrackingGeometry);
1463     if (debug_) {
1464       edm::LogVerbatim("GlobalTrackerMuonAlignment")
1465           << " theTrackingGeometry.isValid() " << theTrackingGeometry.isValid();
1466       edm::LogVerbatim("GlobalTrackerMuonAlignment") << "get also the MuonTransientTrackingRecHitBuilder"
1467                                                      << "\n";
1468       edm::LogVerbatim("GlobalTrackerMuonAlignment") << "get also the TransientTrackingRecHitBuilder"
1469                                                      << "\n";
1470     }
1471     defineFitter = false;
1472   }
1473 
1474   // ................................................ selected/global muon
1475   //itMuon -->  Jim's globalMuon
1476   for (MuonCollection::const_iterator itMuon = smuons->begin(); itMuon != smuons->end(); ++itMuon) {
1477     if (debug_) {
1478       edm::LogVerbatim("GlobalTrackerMuonAlignment")
1479           << " mu gM is GM Mu SaM tM " << itMuon->isGlobalMuon() << " " << itMuon->isMuon() << " "
1480           << itMuon->isStandAloneMuon() << " " << itMuon->isTrackerMuon() << " ";
1481     }
1482 
1483     // get information about the innermost muon hit -------------------------
1484     TransientTrack muTT(itMuon->outerTrack(), magneticField_, trackingGeometry_);
1485     TrajectoryStateOnSurface innerMuTSOS = muTT.innermostMeasurementState();
1486     TrajectoryStateOnSurface outerMuTSOS = muTT.outermostMeasurementState();
1487 
1488     // get information about the outermost tracker hit -----------------------
1489     TransientTrack trackTT(itMuon->track(), magneticField_, trackingGeometry_);
1490     TrajectoryStateOnSurface outerTrackTSOS = trackTT.outermostMeasurementState();
1491     TrajectoryStateOnSurface innerTrackTSOS = trackTT.innermostMeasurementState();
1492 
1493     GlobalPoint pointTrackIn = innerTrackTSOS.globalPosition();
1494     GlobalPoint pointTrackOut = outerTrackTSOS.globalPosition();
1495     float lenghtTrack = (pointTrackOut - pointTrackIn).mag();
1496     GlobalPoint pointMuonIn = innerMuTSOS.globalPosition();
1497     GlobalPoint pointMuonOut = outerMuTSOS.globalPosition();
1498     float lenghtMuon = (pointMuonOut - pointMuonIn).mag();
1499     GlobalVector momentumTrackOut = outerTrackTSOS.globalMomentum();
1500     GlobalVector momentumTrackIn = innerTrackTSOS.globalMomentum();
1501     float distanceInIn = (pointTrackIn - pointMuonIn).mag();
1502     float distanceInOut = (pointTrackIn - pointMuonOut).mag();
1503     float distanceOutIn = (pointTrackOut - pointMuonIn).mag();
1504     float distanceOutOut = (pointTrackOut - pointMuonOut).mag();
1505 
1506     if (debug_) {
1507       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " pointTrackIn " << pointTrackIn;
1508       edm::LogVerbatim("GlobalTrackerMuonAlignment") << "          Out " << pointTrackOut << " lenght " << lenghtTrack;
1509       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " point MuonIn " << pointMuonIn;
1510       edm::LogVerbatim("GlobalTrackerMuonAlignment") << "          Out " << pointMuonOut << " lenght " << lenghtMuon;
1511       edm::LogVerbatim("GlobalTrackerMuonAlignment")
1512           << " momeTrackIn Out " << momentumTrackIn << " " << momentumTrackOut;
1513       edm::LogVerbatim("GlobalTrackerMuonAlignment")
1514           << " doi io ii oo " << distanceOutIn << " " << distanceInOut << " " << distanceInIn << " " << distanceOutOut;
1515     }
1516 
1517     if (lenghtTrack < 90.)
1518       continue;
1519     if (lenghtMuon < 300.)
1520       continue;
1521     if (innerMuTSOS.globalMomentum().mag() < 5. || outerMuTSOS.globalMomentum().mag() < 5.)
1522       continue;
1523     if (momentumTrackIn.mag() < 15. || momentumTrackOut.mag() < 15.)
1524       continue;
1525     if (trackTT.charge() != muTT.charge())
1526       continue;
1527 
1528     if (debug_)
1529       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " passed lenght momentum cuts";
1530 
1531     GlobalVector GRm, GPm, Nl, Rm, Pm, Rt, Pt, Rt0;
1532     AlgebraicSymMatrix66 Cm, C0, Ce, C1;
1533 
1534     //extrapolationT = propagator->propagate(outerTrackTSOS, refSurface);
1535     TrajectoryStateOnSurface extrapolationT;
1536     int tsosMuonIf = 0;
1537 
1538     TrajectoryStateOnSurface muonFittedTSOS;
1539     TrajectoryStateOnSurface trackFittedTSOS;
1540 
1541     if (isolatedMuonMode_) {  //------------------------------- Isolated Muon --- Out-In --
1542       if (debug_)
1543         edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ------ Isolated (out-in) ----- ";
1544       const Surface& refSurface = innerMuTSOS.surface();
1545       ConstReferenceCountingPointer<TangentPlane> tpMuLocal(refSurface.tangentPlane(innerMuTSOS.localPosition()));
1546       Nl = tpMuLocal->normalVector();
1547       ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
1548       GlobalVector Ng = tpMuGlobal->normalVector();
1549       Surface* surf = (Surface*)&refSurface;
1550       const Plane* refPlane = dynamic_cast<Plane*>(surf);
1551       GlobalVector Nlp = refPlane->normalVector();
1552 
1553       if (debug_) {
1554         edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1555                                                        << " alfa " << int(asin(Nl.x()) * 57.29578);
1556         edm::LogVerbatim("GlobalTrackerMuonAlignment") << "  global " << Ng.x() << " " << Ng.y() << " " << Ng.z();
1557         edm::LogVerbatim("GlobalTrackerMuonAlignment") << "      lp " << Nlp.x() << " " << Nlp.y() << " " << Nlp.z();
1558         edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug")
1559             << " Nlocal Nglobal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " " << Ng.x() << " " << Ng.y() << " "
1560             << Ng.z() << " alfa " << static_cast<int>(asin(Nl.x()) * 57.29578);
1561       }
1562 
1563       //                          extrapolation to innermost muon hit
1564 
1565       //extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
1566       if (!refitTrack_)
1567         extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
1568       else {
1569         GlobalTrackerMuonAlignment::trackFitter(itMuon->track(), trackTT, alongMomentum, trackFittedTSOS);
1570         if (trackFittedTSOS.isValid())
1571           extrapolationT = alongSmPr.propagate(trackFittedTSOS, refSurface);
1572       }
1573 
1574       if (!extrapolationT.isValid()) {
1575         if (false & alarm)
1576           edm::LogVerbatim("GlobalTrackerMuonAlignment") << " !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
1577               //<<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
1578               ;
1579         continue;
1580       }
1581       tsosMuonIf = 1;
1582       Rt = GlobalVector((extrapolationT.globalPosition()).x(),
1583                         (extrapolationT.globalPosition()).y(),
1584                         (extrapolationT.globalPosition()).z());
1585 
1586       Pt = extrapolationT.globalMomentum();
1587       //                          global parameters of muon
1588       GRm = GlobalVector(
1589           (innerMuTSOS.globalPosition()).x(), (innerMuTSOS.globalPosition()).y(), (innerMuTSOS.globalPosition()).z());
1590       GPm = innerMuTSOS.globalMomentum();
1591 
1592       Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
1593                          (outerTrackTSOS.globalPosition()).y(),
1594                          (outerTrackTSOS.globalPosition()).z());
1595       Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + innerMuTSOS.cartesianError().matrix());
1596       C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
1597       Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
1598       C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
1599 
1600       if (refitMuon_)
1601         GlobalTrackerMuonAlignment::muonFitter(itMuon->outerTrack(), muTT, oppositeToMomentum, muonFittedTSOS);
1602 
1603     }  //                    -------------------------------   end Isolated Muon -- Out - In  ---
1604 
1605     if (cosmicMuonMode_) {  //------------------------------- Cosmic Muon -----
1606 
1607       if ((distanceOutIn <= distanceInOut) & (distanceOutIn <= distanceInIn) &
1608           (distanceOutIn <= distanceOutOut)) {  // ----- Out - In ------
1609         if (debug_)
1610           edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -----    Out - In -----";
1611 
1612         const Surface& refSurface = innerMuTSOS.surface();
1613         ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
1614         Nl = tpMuGlobal->normalVector();
1615 
1616         //                          extrapolation to innermost muon hit
1617         //extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
1618         if (!refitTrack_)
1619           extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
1620         else {
1621           GlobalTrackerMuonAlignment::trackFitter(itMuon->track(), trackTT, alongMomentum, trackFittedTSOS);
1622           if (trackFittedTSOS.isValid())
1623             extrapolationT = alongSmPr.propagate(trackFittedTSOS, refSurface);
1624         }
1625 
1626         if (!extrapolationT.isValid()) {
1627           if (false & alarm)
1628             edm::LogVerbatim("GlobalTrackerMuonAlignment") << " !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
1629                 //<<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
1630                 ;
1631           continue;
1632         }
1633         if (debug_)
1634           edm::LogVerbatim("GlobalTrackerMuonAlignment") << " extrapolationT.isValid " << extrapolationT.isValid();
1635 
1636         tsosMuonIf = 1;
1637         Rt = GlobalVector((extrapolationT.globalPosition()).x(),
1638                           (extrapolationT.globalPosition()).y(),
1639                           (extrapolationT.globalPosition()).z());
1640 
1641         Pt = extrapolationT.globalMomentum();
1642         //                          global parameters of muon
1643         GRm = GlobalVector(
1644             (innerMuTSOS.globalPosition()).x(), (innerMuTSOS.globalPosition()).y(), (innerMuTSOS.globalPosition()).z());
1645         GPm = innerMuTSOS.globalMomentum();
1646         Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
1647                            (outerTrackTSOS.globalPosition()).y(),
1648                            (outerTrackTSOS.globalPosition()).z());
1649         Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + innerMuTSOS.cartesianError().matrix());
1650         C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
1651         Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
1652         C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
1653 
1654         if (refitMuon_)
1655           GlobalTrackerMuonAlignment::muonFitter(itMuon->outerTrack(), muTT, oppositeToMomentum, muonFittedTSOS);
1656 
1657         if (false & debug_) {
1658           edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ->propDir " << propagator->propagationDirection();
1659           PropagationDirectionChooser Chooser;
1660           edm::LogVerbatim("GlobalTrackerMuonAlignment")
1661               << " propDirCh " << Chooser.operator()(*outerTrackTSOS.freeState(), refSurface) << " Ch == along "
1662               << (alongMomentum == Chooser.operator()(*outerTrackTSOS.freeState(), refSurface));
1663           edm::LogVerbatim("GlobalTrackerMuonAlignment")
1664               << " --- Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1665               << " alfa " << int(asin(Nl.x()) * 57.29578);
1666         }
1667       }  //                                       enf of ---- Out - In -----
1668 
1669       else if ((distanceInOut <= distanceInIn) & (distanceInOut <= distanceOutIn) &
1670                (distanceInOut <= distanceOutOut)) {  // ----- In - Out ------
1671         if (debug_)
1672           edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -----    In - Out -----";
1673 
1674         const Surface& refSurface = outerMuTSOS.surface();
1675         ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
1676         Nl = tpMuGlobal->normalVector();
1677 
1678         //                          extrapolation to outermost muon hit
1679         //extrapolationT = oppositeSmPr.propagate(innerTrackTSOS, refSurface);
1680         if (!refitTrack_)
1681           extrapolationT = oppositeSmPr.propagate(innerTrackTSOS, refSurface);
1682         else {
1683           GlobalTrackerMuonAlignment::trackFitter(itMuon->track(), trackTT, oppositeToMomentum, trackFittedTSOS);
1684           if (trackFittedTSOS.isValid())
1685             extrapolationT = oppositeSmPr.propagate(trackFittedTSOS, refSurface);
1686         }
1687 
1688         if (!extrapolationT.isValid()) {
1689           if (false & alarm)
1690             edm::LogVerbatim("GlobalTrackerMuonAlignment") << " !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
1691                                                            << "\a\a\a\a\a\a\a\a" << extrapolationT.isValid();
1692           continue;
1693         }
1694         if (debug_)
1695           edm::LogVerbatim("GlobalTrackerMuonAlignment") << " extrapolationT.isValid " << extrapolationT.isValid();
1696 
1697         tsosMuonIf = 2;
1698         Rt = GlobalVector((extrapolationT.globalPosition()).x(),
1699                           (extrapolationT.globalPosition()).y(),
1700                           (extrapolationT.globalPosition()).z());
1701 
1702         Pt = extrapolationT.globalMomentum();
1703         //                          global parameters of muon
1704         GRm = GlobalVector(
1705             (outerMuTSOS.globalPosition()).x(), (outerMuTSOS.globalPosition()).y(), (outerMuTSOS.globalPosition()).z());
1706         GPm = outerMuTSOS.globalMomentum();
1707         Rt0 = GlobalVector((innerTrackTSOS.globalPosition()).x(),
1708                            (innerTrackTSOS.globalPosition()).y(),
1709                            (innerTrackTSOS.globalPosition()).z());
1710         Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + outerMuTSOS.cartesianError().matrix());
1711         C0 = AlgebraicSymMatrix66(innerTrackTSOS.cartesianError().matrix());
1712         Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
1713         C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
1714 
1715         if (refitMuon_)
1716           GlobalTrackerMuonAlignment::muonFitter(itMuon->outerTrack(), muTT, alongMomentum, muonFittedTSOS);
1717 
1718         if (false & debug_) {
1719           edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ->propDir " << propagator->propagationDirection();
1720           PropagationDirectionChooser Chooser;
1721           edm::LogVerbatim("GlobalTrackerMuonAlignment")
1722               << " propDirCh " << Chooser.operator()(*innerTrackTSOS.freeState(), refSurface) << " Ch == oppisite "
1723               << (oppositeToMomentum == Chooser.operator()(*innerTrackTSOS.freeState(), refSurface));
1724           edm::LogVerbatim("GlobalTrackerMuonAlignment")
1725               << " --- Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1726               << " alfa " << int(asin(Nl.x()) * 57.29578);
1727         }
1728       }  //                                        enf of ---- In - Out -----
1729 
1730       else if ((distanceOutOut <= distanceInOut) & (distanceOutOut <= distanceInIn) &
1731                (distanceOutOut <= distanceOutIn)) {  // ----- Out - Out ------
1732         if (debug_)
1733           edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -----    Out - Out -----";
1734 
1735         // reject: momentum of track has opposite direction to muon track
1736         continue;
1737 
1738         const Surface& refSurface = outerMuTSOS.surface();
1739         ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
1740         Nl = tpMuGlobal->normalVector();
1741 
1742         //                          extrapolation to outermost muon hit
1743         extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
1744 
1745         if (!extrapolationT.isValid()) {
1746           if (alarm)
1747             edm::LogVerbatim("GlobalTrackerMuonAlignment") << " !!!!!!!!! Catastrophe Out-Out extrapolationT.isValid "
1748                                                            << "\a\a\a\a\a\a\a\a" << extrapolationT.isValid();
1749           continue;
1750         }
1751         if (debug_)
1752           edm::LogVerbatim("GlobalTrackerMuonAlignment") << " extrapolationT.isValid " << extrapolationT.isValid();
1753 
1754         tsosMuonIf = 2;
1755         Rt = GlobalVector((extrapolationT.globalPosition()).x(),
1756                           (extrapolationT.globalPosition()).y(),
1757                           (extrapolationT.globalPosition()).z());
1758 
1759         Pt = extrapolationT.globalMomentum();
1760         //                          global parameters of muon
1761         GRm = GlobalVector(
1762             (outerMuTSOS.globalPosition()).x(), (outerMuTSOS.globalPosition()).y(), (outerMuTSOS.globalPosition()).z());
1763         GPm = outerMuTSOS.globalMomentum();
1764         Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
1765                            (outerTrackTSOS.globalPosition()).y(),
1766                            (outerTrackTSOS.globalPosition()).z());
1767         Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + outerMuTSOS.cartesianError().matrix());
1768         C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
1769         Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
1770         C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
1771 
1772         if (debug_) {
1773           PropagationDirectionChooser Chooser;
1774           edm::LogVerbatim("GlobalTrackerMuonAlignment")
1775               << " propDirCh " << Chooser.operator()(*outerTrackTSOS.freeState(), refSurface) << " Ch == along "
1776               << (alongMomentum == Chooser.operator()(*outerTrackTSOS.freeState(), refSurface));
1777           edm::LogVerbatim("GlobalTrackerMuonAlignment")
1778               << " --- Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1779               << " alfa " << int(asin(Nl.x()) * 57.29578);
1780           edm::LogVerbatim("GlobalTrackerMuonAlignment")
1781               << "     Nornal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1782               << " alfa " << int(asin(Nl.x()) * 57.29578);
1783         }
1784       }  //                                       enf of ---- Out - Out -----
1785       else {
1786         if (alarm)
1787           edm::LogVerbatim("GlobalTrackerMuonAlignment")
1788               << " ------- !!!!!!!!!!!!!!! No proper Out - In \a\a\a\a\a\a\a";
1789         continue;
1790       }
1791 
1792     }  //                     -------------------------------   end Cosmic Muon -----
1793 
1794     if (tsosMuonIf == 0) {
1795       if (info) {
1796         edm::LogVerbatim("GlobalTrackerMuonAlignment") << "No tsosMuon !!!!!!";
1797         continue;
1798       }
1799     }
1800     TrajectoryStateOnSurface tsosMuon;
1801     if (tsosMuonIf == 1)
1802       tsosMuon = muTT.innermostMeasurementState();
1803     else
1804       tsosMuon = muTT.outermostMeasurementState();
1805 
1806     //GlobalTrackerMuonAlignment::misalignMuon(GRm, GPm, Nl, Rt, Rm, Pm);
1807     AlgebraicVector4 LPRm;  // muon local (dx/dz, dy/dz, x, y)
1808     GlobalTrackerMuonAlignment::misalignMuonL(GRm, GPm, Nl, Rt, Rm, Pm, LPRm, extrapolationT, tsosMuon);
1809 
1810     if (refitTrack_) {
1811       if (!trackFittedTSOS.isValid()) {
1812         if (info)
1813           edm::LogVerbatim("GlobalTrackerMuonAlignment") << " =================  trackFittedTSOS notValid !!!!!!!! ";
1814         continue;
1815       }
1816       if (debug_)
1817         this->debugTrajectorySOS(" trackFittedTSOS ", trackFittedTSOS);
1818     }
1819 
1820     if (refitMuon_) {
1821       if (!muonFittedTSOS.isValid()) {
1822         if (info)
1823           edm::LogVerbatim("GlobalTrackerMuonAlignment") << " =================  muonFittedTSOS notValid !!!!!!!! ";
1824         continue;
1825       }
1826       if (debug_)
1827         this->debugTrajectorySOS(" muonFittedTSOS ", muonFittedTSOS);
1828       Rm = GlobalVector((muonFittedTSOS.globalPosition()).x(),
1829                         (muonFittedTSOS.globalPosition()).y(),
1830                         (muonFittedTSOS.globalPosition()).z());
1831       Pm = muonFittedTSOS.globalMomentum();
1832       LPRm = AlgebraicVector4(muonFittedTSOS.localParameters().vector()(1),
1833                               muonFittedTSOS.localParameters().vector()(2),
1834                               muonFittedTSOS.localParameters().vector()(3),
1835                               muonFittedTSOS.localParameters().vector()(4));
1836     }
1837     GlobalVector resR = Rm - Rt;
1838     GlobalVector resP0 = Pm - Pt;
1839     GlobalVector resP = Pm / Pm.mag() - Pt / Pt.mag();
1840     float RelMomResidual = (Pm.mag() - Pt.mag()) / (Pt.mag() + 1.e-6);
1841     ;
1842 
1843     AlgebraicVector6 Vm;
1844     Vm(0) = resR.x();
1845     Vm(1) = resR.y();
1846     Vm(2) = resR.z();
1847     Vm(3) = resP0.x();
1848     Vm(4) = resP0.y();
1849     Vm(5) = resP0.z();
1850     float Rmuon = Rm.perp();
1851     float Zmuon = Rm.z();
1852     float alfa_x = atan2(Nl.x(), Nl.y()) * 57.29578;
1853 
1854     if (debug_) {
1855       edm::LogVerbatim("GlobalTrackerMuonAlignment")
1856           << " Nx Ny Nz alfa_x " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " " << alfa_x;
1857       edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug")
1858           << " Rm " << Rm << "\n Rt " << Rt << "\n resR " << resR << "\n resP " << resP << " dp/p " << RelMomResidual;
1859     }
1860 
1861     double chi_d = 0;
1862     for (int i = 0; i <= 5; i++)
1863       chi_d += Vm(i) * Vm(i) / Cm(i, i);
1864 
1865     AlgebraicVector5 Vml(tsosMuon.localParameters().vector() - extrapolationT.localParameters().vector());
1866     AlgebraicSymMatrix55 m(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
1867     AlgebraicSymMatrix55 Cml(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
1868     bool ierrLoc = !m.Invert();
1869     if (ierrLoc && debug_ && info) {
1870       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ==== Error inverting Local covariance matrix ==== ";
1871       continue;
1872     }
1873     double chi_Loc = ROOT::Math::Similarity(Vml, m);
1874     if (debug_)
1875       edm::LogVerbatim("GlobalTrackerMuonAlignment")
1876           << " chi_Loc px/pz/err " << chi_Loc << " " << Vml(1) / std::sqrt(Cml(1, 1));
1877 
1878     if (Pt.mag() < 15.)
1879       continue;
1880     if (Pm.mag() < 5.)
1881       continue;
1882 
1883     //if(Pt.mag() < 30.) continue;          // momenum cut  < 30GeV
1884     //if(Pt.mag() < 60.) continue;          // momenum cut  < 30GeV
1885     //if(Pt.mag() > 50.) continue;          //  momenum cut > 50GeV
1886     //if(Pt.mag() > 100.) continue;          //  momenum cut > 100GeV
1887     //if(trackTT.charge() < 0) continue;    // select positive charge
1888     //if(trackTT.charge() > 0) continue;    // select negative charge
1889 
1890     //if(fabs(resR.x()) > 5.) continue;     // strong cut  X
1891     //if(fabs(resR.y()) > 5.) continue;     //             Y
1892     //if(fabs(resR.z()) > 5.) continue;     //             Z
1893     //if(fabs(resR.mag()) > 7.5) continue;   //            dR
1894 
1895     //if(fabs(RelMomResidual) > 0.5) continue;
1896     if (fabs(resR.x()) > 20.)
1897       continue;
1898     if (fabs(resR.y()) > 20.)
1899       continue;
1900     if (fabs(resR.z()) > 20.)
1901       continue;
1902     if (fabs(resR.mag()) > 30.)
1903       continue;
1904     if (fabs(resP.x()) > 0.06)
1905       continue;
1906     if (fabs(resP.y()) > 0.06)
1907       continue;
1908     if (fabs(resP.z()) > 0.06)
1909       continue;
1910     if (chi_d > 40.)
1911       continue;
1912 
1913     //                                            select Barrel
1914     //if(Rmuon < 400. || Rmuon > 450.) continue;
1915     //if(Zmuon < -600. || Zmuon > 600.) continue;
1916     //if(fabs(Nl.z()) > 0.95) continue;
1917     //MuSelect = " Barrel";
1918     //                                                  EndCap1
1919     //if(Rmuon < 120. || Rmuon > 450.) continue;
1920     //if(Zmuon < -720.) continue;
1921     //if(Zmuon > -580.) continue;
1922     //if(fabs(Nl.z()) < 0.95) continue;
1923     //MuSelect = " EndCap1";
1924     //                                                  EndCap2
1925     //if(Rmuon < 120. || Rmuon > 450.) continue;
1926     //if(Zmuon >  720.) continue;
1927     //if(Zmuon <  580.) continue;
1928     //if(fabs(Nl.z()) < 0.95) continue;
1929     //MuSelect = " EndCap2";
1930     //                                                 select All
1931     if (Rmuon < 120. || Rmuon > 450.)
1932       continue;
1933     if (Zmuon < -720. || Zmuon > 720.)
1934       continue;
1935     MuSelect = " Barrel+EndCaps";
1936 
1937     if (debug_)
1938       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " .............. passed all cuts";
1939 
1940     N_track++;
1941     //                     gradient and Hessian for each track
1942 
1943     GlobalTrackerMuonAlignment::gradientGlobalAlg(Rt, Pt, Rm, Nl, Cm);
1944     GlobalTrackerMuonAlignment::gradientGlobal(Rt, Pt, Rm, Pm, Nl, Cm);
1945 
1946     CLHEP::HepSymMatrix covLoc(4, 0);
1947     for (int i = 1; i <= 4; i++)
1948       for (int j = 1; j <= i; j++) {
1949         covLoc(i, j) = (tsosMuon.localError().matrix() + extrapolationT.localError().matrix())(i, j);
1950         //if(i != j) Cov(i,j) = 0.;
1951       }
1952 
1953     const Surface& refSurface = tsosMuon.surface();
1954     CLHEP::HepMatrix rotLoc(3, 3, 0);
1955     rotLoc(1, 1) = refSurface.rotation().xx();
1956     rotLoc(1, 2) = refSurface.rotation().xy();
1957     rotLoc(1, 3) = refSurface.rotation().xz();
1958 
1959     rotLoc(2, 1) = refSurface.rotation().yx();
1960     rotLoc(2, 2) = refSurface.rotation().yy();
1961     rotLoc(2, 3) = refSurface.rotation().yz();
1962 
1963     rotLoc(3, 1) = refSurface.rotation().zx();
1964     rotLoc(3, 2) = refSurface.rotation().zy();
1965     rotLoc(3, 3) = refSurface.rotation().zz();
1966 
1967     CLHEP::HepVector posLoc(3);
1968     posLoc(1) = refSurface.position().x();
1969     posLoc(2) = refSurface.position().y();
1970     posLoc(3) = refSurface.position().z();
1971 
1972     GlobalTrackerMuonAlignment::gradientLocal(Rt, Pt, Rm, Pm, Nl, covLoc, rotLoc, posLoc, LPRm);
1973 
1974     if (debug_) {
1975       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Norm   " << Nl;
1976       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " posLoc " << posLoc.T();
1977       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " rotLoc " << rotLoc;
1978     }
1979 
1980     // -----------------------------------------------------  fill histogram
1981     histo->Fill(itMuon->track()->pt());
1982 
1983     //histo2->Fill(itMuon->track()->outerP());
1984     histo2->Fill(Pt.mag());
1985     histo3->Fill((PI / 2. - itMuon->track()->outerTheta()));
1986     histo4->Fill(itMuon->track()->phi());
1987     histo5->Fill(Rmuon);
1988     histo6->Fill(Zmuon);
1989     histo7->Fill(RelMomResidual);
1990     //histo8->Fill(chi);
1991     histo8->Fill(chi_d);
1992 
1993     histo101->Fill(Zmuon, Rmuon);
1994     histo101->Fill(Rt0.z(), Rt0.perp());
1995     histo102->Fill(Rt0.x(), Rt0.y());
1996     histo102->Fill(Rm.x(), Rm.y());
1997 
1998     histo11->Fill(resR.mag());
1999     if (fabs(Nl.x()) < 0.98)
2000       histo12->Fill(resR.x());
2001     if (fabs(Nl.y()) < 0.98)
2002       histo13->Fill(resR.y());
2003     if (fabs(Nl.z()) < 0.98)
2004       histo14->Fill(resR.z());
2005     histo15->Fill(resP.x());
2006     histo16->Fill(resP.y());
2007     histo17->Fill(resP.z());
2008 
2009     if ((fabs(Nl.x()) < 0.98) && (fabs(Nl.y()) < 0.98) && (fabs(Nl.z()) < 0.98)) {
2010       histo18->Fill(std::sqrt(C0(0, 0)));
2011       histo19->Fill(std::sqrt(C1(0, 0)));
2012       histo20->Fill(std::sqrt(C1(0, 0) + Ce(0, 0)));
2013     }
2014     if (fabs(Nl.x()) < 0.98)
2015       histo21->Fill(Vm(0) / std::sqrt(Cm(0, 0)));
2016     if (fabs(Nl.y()) < 0.98)
2017       histo22->Fill(Vm(1) / std::sqrt(Cm(1, 1)));
2018     if (fabs(Nl.z()) < 0.98)
2019       histo23->Fill(Vm(2) / std::sqrt(Cm(2, 2)));
2020     histo24->Fill(Vm(3) / std::sqrt(C1(3, 3) + Ce(3, 3)));
2021     histo25->Fill(Vm(4) / std::sqrt(C1(4, 4) + Ce(4, 4)));
2022     histo26->Fill(Vm(5) / std::sqrt(C1(5, 5) + Ce(5, 5)));
2023     histo27->Fill(Nl.x());
2024     histo28->Fill(Nl.y());
2025     histo29->Fill(lenghtTrack);
2026     histo30->Fill(lenghtMuon);
2027     histo31->Fill(chi_Loc);
2028     histo32->Fill(Vml(1) / std::sqrt(Cml(1, 1)));
2029     histo33->Fill(Vml(2) / std::sqrt(Cml(2, 2)));
2030     histo34->Fill(Vml(3) / std::sqrt(Cml(3, 3)));
2031     histo35->Fill(Vml(4) / std::sqrt(Cml(4, 4)));
2032 
2033     if (debug_) {  //--------------------------------- debug print ----------
2034 
2035       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " diag 0[ " << C0(0, 0) << " " << C0(1, 1) << " " << C0(2, 2)
2036                                                      << " " << C0(3, 3) << " " << C0(4, 4) << " " << C0(5, 5) << " ]";
2037       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " diag e[ " << Ce(0, 0) << " " << Ce(1, 1) << " " << Ce(2, 2)
2038                                                      << " " << Ce(3, 3) << " " << Ce(4, 4) << " " << Ce(5, 5) << " ]";
2039       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " diag 1[ " << C1(0, 0) << " " << C1(1, 1) << " " << C1(2, 2)
2040                                                      << " " << C1(3, 3) << " " << C1(4, 4) << " " << C1(5, 5) << " ]";
2041       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Rm   " << Rm.x() << " " << Rm.y() << " " << Rm.z() << " Pm   "
2042                                                      << Pm.x() << " " << Pm.y() << " " << Pm.z();
2043       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Rt   " << Rt.x() << " " << Rt.y() << " " << Rt.z() << " Pt   "
2044                                                      << Pt.x() << " " << Pt.y() << " " << Pt.z();
2045       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Nl*(Rm-Rt) " << Nl.dot(Rm - Rt);
2046       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " resR " << resR.x() << " " << resR.y() << " " << resR.z()
2047                                                      << " resP " << resP.x() << " " << resP.y() << " " << resP.z();
2048       edm::LogVerbatim("GlobalTrackerMuonAlignment")
2049           << " Rm-t " << (Rm - Rt).x() << " " << (Rm - Rt).y() << " " << (Rm - Rt).z() << " Pm-t " << (Pm - Pt).x()
2050           << " " << (Pm - Pt).y() << " " << (Pm - Pt).z();
2051       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Vm   " << Vm;
2052       edm::LogVerbatim("GlobalTrackerMuonAlignment")
2053           << " +-   " << std::sqrt(Cm(0, 0)) << " " << std::sqrt(Cm(1, 1)) << " " << std::sqrt(Cm(2, 2)) << " "
2054           << std::sqrt(Cm(3, 3)) << " " << std::sqrt(Cm(4, 4)) << " " << std::sqrt(Cm(5, 5));
2055       edm::LogVerbatim("GlobalTrackerMuonAlignment")
2056           << " Pmuon Ptrack dP/Ptrack " << itMuon->outerTrack()->p() << " " << itMuon->track()->outerP() << " "
2057           << (itMuon->outerTrack()->p() - itMuon->track()->outerP()) / itMuon->track()->outerP();
2058       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " cov matrix ";
2059       edm::LogVerbatim("GlobalTrackerMuonAlignment") << Cm;
2060       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " diag [ " << Cm(0, 0) << " " << Cm(1, 1) << " " << Cm(2, 2)
2061                                                      << " " << Cm(3, 3) << " " << Cm(4, 4) << " " << Cm(5, 5) << " ]";
2062 
2063       AlgebraicSymMatrix66 Ro;
2064       double Diag[6];
2065       for (int i = 0; i <= 5; i++)
2066         Diag[i] = std::sqrt(Cm(i, i));
2067       for (int i = 0; i <= 5; i++)
2068         for (int j = 0; j <= 5; j++)
2069           Ro(i, j) = Cm(i, j) / Diag[i] / Diag[j];
2070       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " correlation matrix ";
2071       edm::LogVerbatim("GlobalTrackerMuonAlignment") << Ro;
2072 
2073       AlgebraicSymMatrix66 CmI;
2074       for (int i = 0; i <= 5; i++)
2075         for (int j = 0; j <= 5; j++)
2076           CmI(i, j) = Cm(i, j);
2077 
2078       bool ierr = !CmI.Invert();
2079       if (ierr) {
2080         if (alarm || debug_)
2081           edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Error inverse covariance matrix !!!!!!!!!!!";
2082         continue;
2083       }
2084       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " inverse cov matrix ";
2085       edm::LogVerbatim("GlobalTrackerMuonAlignment") << Cm;
2086 
2087       double chi = ROOT::Math::Similarity(Vm, CmI);
2088       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " chi chi_d " << chi << " " << chi_d;
2089     }  // end of debug_ printout  --------------------------------------------
2090 
2091   }  // end loop on selected muons, i.e. Jim's globalMuon
2092 
2093 }  //end of analyzeTrackTrajectory
2094 
2095 // ----  calculate gradient and Hessian matrix (algebraic) to search global shifts ------
2096 void GlobalTrackerMuonAlignment::gradientGlobalAlg(
2097     GlobalVector& Rt, GlobalVector& Pt, GlobalVector& Rm, GlobalVector& Nl, AlgebraicSymMatrix66& Cm) {
2098   // ---------------------------- Calculate Information matrix and Gfree vector
2099   // Information == Hessian , Gfree == gradient of the objective function
2100 
2101   AlgebraicMatrix33 Jac;
2102   AlgebraicVector3 Wi, R_m, R_t, P_t, Norm, dR;
2103 
2104   R_m(0) = Rm.x();
2105   R_m(1) = Rm.y();
2106   R_m(2) = Rm.z();
2107   R_t(0) = Rt.x();
2108   R_t(1) = Rt.y();
2109   R_t(2) = Rt.z();
2110   P_t(0) = Pt.x();
2111   P_t(1) = Pt.y();
2112   P_t(2) = Pt.z();
2113   Norm(0) = Nl.x();
2114   Norm(1) = Nl.y();
2115   Norm(2) = Nl.z();
2116 
2117   for (int i = 0; i <= 2; i++) {
2118     if (Cm(i, i) > 1.e-20)
2119       Wi(i) = 1. / Cm(i, i);
2120     else
2121       Wi(i) = 1.e-10;
2122     dR(i) = R_m(i) - R_t(i);
2123   }
2124 
2125   float PtN = P_t(0) * Norm(0) + P_t(1) * Norm(1) + P_t(2) * Norm(2);
2126 
2127   Jac(0, 0) = 1. - P_t(0) * Norm(0) / PtN;
2128   Jac(0, 1) = -P_t(0) * Norm(1) / PtN;
2129   Jac(0, 2) = -P_t(0) * Norm(2) / PtN;
2130 
2131   Jac(1, 0) = -P_t(1) * Norm(0) / PtN;
2132   Jac(1, 1) = 1. - P_t(1) * Norm(1) / PtN;
2133   Jac(1, 2) = -P_t(1) * Norm(2) / PtN;
2134 
2135   Jac(2, 0) = -P_t(2) * Norm(0) / PtN;
2136   Jac(2, 1) = -P_t(2) * Norm(1) / PtN;
2137   Jac(2, 2) = 1. - P_t(2) * Norm(2) / PtN;
2138 
2139   AlgebraicSymMatrix33 Itr;
2140 
2141   for (int i = 0; i <= 2; i++)
2142     for (int j = 0; j <= 2; j++) {
2143       if (j < i)
2144         continue;
2145       Itr(i, j) = 0.;
2146       edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ij " << i << " " << j;
2147       for (int k = 0; k <= 2; k++) {
2148         Itr(i, j) += Jac(k, i) * Wi(k) * Jac(k, j);
2149       }
2150     }
2151 
2152   for (int i = 0; i <= 2; i++)
2153     for (int j = 0; j <= 2; j++) {
2154       if (j < i)
2155         continue;
2156       Inf(i, j) += Itr(i, j);
2157     }
2158 
2159   AlgebraicVector3 Gtr(0., 0., 0.);
2160   for (int i = 0; i <= 2; i++)
2161     for (int k = 0; k <= 2; k++)
2162       Gtr(i) += dR(k) * Wi(k) * Jac(k, i);
2163   for (int i = 0; i <= 2; i++)
2164     Gfr(i) += Gtr(i);
2165 
2166   if (debug_) {
2167     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Wi  " << Wi;
2168     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " N   " << Norm;
2169     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " P_t " << P_t;
2170     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " (Pt*N) " << PtN;
2171     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " dR  " << dR;
2172     edm::LogVerbatim("GlobalTrackerMuonAlignment")
2173         << " +/- " << 1. / std::sqrt(Wi(0)) << " " << 1. / std::sqrt(Wi(1)) << " " << 1. / std::sqrt(Wi(2)) << " ";
2174     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Jacobian dr/ddx ";
2175     edm::LogVerbatim("GlobalTrackerMuonAlignment") << Jac;
2176     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " G-- " << Gtr;
2177     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Itrack ";
2178     edm::LogVerbatim("GlobalTrackerMuonAlignment") << Itr;
2179     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Gfr " << Gfr;
2180     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -- Inf --";
2181     edm::LogVerbatim("GlobalTrackerMuonAlignment") << Inf;
2182   }
2183 
2184   return;
2185 }
2186 
2187 // ----  calculate gradient and Hessian matrix in global parameters ------
2188 void GlobalTrackerMuonAlignment::gradientGlobal(GlobalVector& GRt,
2189                                                 GlobalVector& GPt,
2190                                                 GlobalVector& GRm,
2191                                                 GlobalVector& GPm,
2192                                                 GlobalVector& GNorm,
2193                                                 AlgebraicSymMatrix66& GCov) {
2194   // we search for 6D global correction vector (d, a), where
2195   //                        3D vector of shihts d
2196   //               3D vector of rotation angles    a
2197 
2198   //bool alarm = true;
2199   bool alarm = false;
2200   bool info = false;
2201 
2202   int Nd = 6;  // dimension of vector of alignment pararmeters, d
2203 
2204   //double PtMom = GPt.mag();
2205   CLHEP::HepSymMatrix w(Nd, 0);
2206   for (int i = 1; i <= Nd; i++)
2207     for (int j = 1; j <= Nd; j++) {
2208       if (j <= i)
2209         w(i, j) = GCov(i - 1, j - 1);
2210       //if(i >= 3) w(i,j) /= PtMom;
2211       //if(j >= 3) w(i,j) /= PtMom;
2212       if ((i == j) && (i <= 3) && (GCov(i - 1, j - 1) < 1.e-20))
2213         w(i, j) = 1.e20;  // w=0
2214       if (i != j)
2215         w(i, j) = 0.;  // use diaginal elements
2216     }
2217 
2218   //GPt /= GPt.mag();
2219   //GPm /= GPm.mag();   // end of transform
2220 
2221   CLHEP::HepVector V(Nd), Rt(3), Pt(3), Rm(3), Pm(3), Norm(3);
2222   Rt(1) = GRt.x();
2223   Rt(2) = GRt.y();
2224   Rt(3) = GRt.z();
2225   Pt(1) = GPt.x();
2226   Pt(2) = GPt.y();
2227   Pt(3) = GPt.z();
2228   Rm(1) = GRm.x();
2229   Rm(2) = GRm.y();
2230   Rm(3) = GRm.z();
2231   Pm(1) = GPm.x();
2232   Pm(2) = GPm.y();
2233   Pm(3) = GPm.z();
2234   Norm(1) = GNorm.x();
2235   Norm(2) = GNorm.y();
2236   Norm(3) = GNorm.z();
2237 
2238   V = dsum(Rm - Rt, Pm - Pt);
2239   edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " V " << V.T();
2240 
2241   //double PmN = CLHEP::dot(Pm, Norm);
2242   double PmN = CLHEP_dot(Pm, Norm);
2243 
2244   CLHEP::HepMatrix Jac(Nd, Nd, 0);
2245   for (int i = 1; i <= 3; i++)
2246     for (int j = 1; j <= 3; j++) {
2247       Jac(i, j) = Pm(i) * Norm(j) / PmN;
2248       if (i == j)
2249         Jac(i, j) -= 1.;
2250     }
2251 
2252   //                                            dp/da
2253   Jac(4, 4) = 0.;      // dpx/dax
2254   Jac(5, 4) = -Pm(3);  // dpy/dax
2255   Jac(6, 4) = Pm(2);   // dpz/dax
2256   Jac(4, 5) = Pm(3);   // dpx/day
2257   Jac(5, 5) = 0.;      // dpy/day
2258   Jac(6, 5) = -Pm(1);  // dpz/day
2259   Jac(4, 6) = -Pm(2);  // dpx/daz
2260   Jac(5, 6) = Pm(1);   // dpy/daz
2261   Jac(6, 6) = 0.;      // dpz/daz
2262 
2263   CLHEP::HepVector dsda(3);
2264   dsda(1) = (Norm(2) * Rm(3) - Norm(3) * Rm(2)) / PmN;
2265   dsda(2) = (Norm(3) * Rm(1) - Norm(1) * Rm(3)) / PmN;
2266   dsda(3) = (Norm(1) * Rm(2) - Norm(2) * Rm(1)) / PmN;
2267 
2268   //                                             dr/da
2269   Jac(1, 4) = Pm(1) * dsda(1);           // drx/dax
2270   Jac(2, 4) = -Rm(3) + Pm(2) * dsda(1);  // dry/dax
2271   Jac(3, 4) = Rm(2) + Pm(3) * dsda(1);   // drz/dax
2272 
2273   Jac(1, 5) = Rm(3) + Pm(1) * dsda(2);   // drx/day
2274   Jac(2, 5) = Pm(2) * dsda(2);           // dry/day
2275   Jac(3, 5) = -Rm(1) + Pm(3) * dsda(2);  // drz/day
2276 
2277   Jac(1, 6) = -Rm(2) + Pm(1) * dsda(3);  // drx/daz
2278   Jac(2, 6) = Rm(1) + Pm(2) * dsda(3);   // dry/daz
2279   Jac(3, 6) = Pm(3) * dsda(3);           // drz/daz
2280 
2281   CLHEP::HepSymMatrix W(Nd, 0);
2282   int ierr;
2283   W = w.inverse(ierr);
2284   if (ierr != 0) {
2285     if (alarm)
2286       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " gradientGlobal: inversion of matrix w fail ";
2287     return;
2288   }
2289 
2290   CLHEP::HepMatrix W_Jac(Nd, Nd, 0);
2291   W_Jac = Jac.T() * W;
2292 
2293   CLHEP::HepVector grad3(Nd);
2294   grad3 = W_Jac * V;
2295 
2296   CLHEP::HepMatrix hess3(Nd, Nd);
2297   hess3 = Jac.T() * W * Jac;
2298   //hess3(4,4) = 1.e-10;  hess3(5,5) = 1.e-10;  hess3(6,6) = 1.e-10; //????????????????
2299 
2300   Grad += grad3;
2301   Hess += hess3;
2302 
2303   CLHEP::HepVector d3I = CLHEP::solve(Hess, -Grad);
2304   int iEr3I;
2305   CLHEP::HepMatrix Errd3I = Hess.inverse(iEr3I);
2306   if (iEr3I != 0) {
2307     if (alarm)
2308       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " gradientGlobal error inverse  Hess matrix !!!!!!!!!!!";
2309   }
2310 
2311   if (info || debug_) {
2312     edm::LogVerbatim("GlobalTrackerMuonAlignment")
2313         << " dG   " << d3I(1) << " " << d3I(2) << " " << d3I(3) << " " << d3I(4) << " " << d3I(5) << " " << d3I(6);
2314     ;
2315     edm::LogVerbatim("GlobalTrackerMuonAlignment")
2316         << " +-   " << std::sqrt(Errd3I(1, 1)) << " " << std::sqrt(Errd3I(2, 2)) << " " << std::sqrt(Errd3I(3, 3))
2317         << " " << std::sqrt(Errd3I(4, 4)) << " " << std::sqrt(Errd3I(5, 5)) << " " << std::sqrt(Errd3I(6, 6));
2318   }
2319 
2320 #ifdef CHECK_OF_DERIVATIVES
2321   //              --------------------    check of derivatives
2322 
2323   CLHEP::HepVector d(3, 0), a(3, 0);
2324   CLHEP::HepMatrix T(3, 3, 1);
2325 
2326   CLHEP::HepMatrix Ti = T.T();
2327   //double A = CLHEP::dot(Ti*Pm, Norm);
2328   //double B = CLHEP::dot((Rt -Ti*Rm + Ti*d), Norm);
2329   double A = CLHEP_dot(Ti * Pm, Norm);
2330   double B = CLHEP_dot((Rt - Ti * Rm + Ti * d), Norm);
2331   double s0 = B / A;
2332 
2333   CLHEP::HepVector r0(3, 0), p0(3, 0);
2334   r0 = Ti * Rm - Ti * d + s0 * (Ti * Pm) - Rt;
2335   p0 = Ti * Pm - Pt;
2336 
2337   double delta = 0.0001;
2338 
2339   int ii = 3;
2340   d(ii) += delta;  // d
2341   //T(2,3) += delta;   T(3,2) -= delta;  int ii = 1; // a1
2342   //T(3,1) += delta;   T(1,3) -= delta;   int ii = 2; // a2
2343   //T(1,2) += delta;   T(2,1) -= delta;   int ii = 3; // a2
2344   Ti = T.T();
2345   //A = CLHEP::dot(Ti*Pm, Norm);
2346   //B = CLHEP::dot((Rt -Ti*Rm + Ti*d), Norm);
2347   A = CLHEP_dot(Ti * Pm, Norm);
2348   B = CLHEP_dot((Rt - Ti * Rm + Ti * d), Norm);
2349   double s = B / A;
2350 
2351   CLHEP::HepVector r(3, 0), p(3, 0);
2352   r = Ti * Rm - Ti * d + s * (Ti * Pm) - Rt;
2353   p = Ti * Pm - Pt;
2354 
2355   edm::LogVerbatim("GlobalTrackerMuonAlignment")
2356       << " s0 s num dsda(" << ii << ") " << s0 << " " << s << " " << (s - s0) / delta << " " << dsda(ii);
2357   // d(r,p) / d shift
2358   edm::LogVerbatim("GlobalTrackerMuonAlignment")
2359       << " -- An d(r,p)/d(" << ii << ") " << Jac(1, ii) << " " << Jac(2, ii) << " " << Jac(3, ii) << " " << Jac(4, ii)
2360       << " " << Jac(5, ii) << " " << Jac(6, ii);
2361   edm::LogVerbatim("GlobalTrackerMuonAlignment")
2362       << "    Nu d(r,p)/d(" << ii << ") " << (r(1) - r0(1)) / delta << " " << (r(2) - r0(2)) / delta << " "
2363       << (r(3) - r0(3)) / delta << " " << (p(1) - p0(1)) / delta << " " << (p(2) - p0(2)) / delta << " "
2364       << (p(3) - p0(3)) / delta;
2365   // d(r,p) / d angle
2366   edm::LogVerbatim("GlobalTrackerMuonAlignment")
2367       << " -- An d(r,p)/a(" << ii << ") " << Jac(1, ii + 3) << " " << Jac(2, ii + 3) << " " << Jac(3, ii + 3) << " "
2368       << Jac(4, ii + 3) << " " << Jac(5, ii + 3) << " " << Jac(6, ii + 3);
2369   edm::LogVerbatim("GlobalTrackerMuonAlignment")
2370       << "    Nu d(r,p)/a(" << ii << ") " << (r(1) - r0(1)) / delta << " " << (r(2) - r0(2)) / delta << " "
2371       << (r(3) - r0(3)) / delta << " " << (p(1) - p0(1)) / delta << " " << (p(2) - p0(2)) / delta << " "
2372       << (p(3) - p0(3)) / delta;
2373   //               -----------------------------  end of check
2374 #endif
2375 
2376   return;
2377 }  // end gradientGlobal
2378 
2379 // ----  calculate gradient and Hessian matrix in local parameters ------
2380 void GlobalTrackerMuonAlignment::gradientLocal(GlobalVector& GRt,
2381                                                GlobalVector& GPt,
2382                                                GlobalVector& GRm,
2383                                                GlobalVector& GPm,
2384                                                GlobalVector& GNorm,
2385                                                CLHEP::HepSymMatrix& covLoc,
2386                                                CLHEP::HepMatrix& rotLoc,
2387                                                CLHEP::HepVector& R0,
2388                                                AlgebraicVector4& LPRm) {
2389   // we search for 6D global correction vector (d, a), where
2390   //                        3D vector of shihts d
2391   //               3D vector of rotation angles    a
2392 
2393   bool alarm = true;
2394   //bool alarm = false;
2395   bool info = false;
2396 
2397   if (debug_)
2398     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " gradientLocal ";
2399 
2400   /*
2401   const Surface& refSurface = tsosMuon.surface();
2402 
2403   CLHEP::HepMatrix rotLoc (3,3,0);
2404   rotLoc(1,1) = refSurface.rotation().xx();
2405   rotLoc(1,2) = refSurface.rotation().xy();
2406   rotLoc(1,3) = refSurface.rotation().xz();
2407   
2408   rotLoc(2,1) = refSurface.rotation().yx();
2409   rotLoc(2,2) = refSurface.rotation().yy();
2410   rotLoc(2,3) = refSurface.rotation().yz();
2411   
2412   rotLoc(3,1) = refSurface.rotation().zx();
2413   rotLoc(3,2) = refSurface.rotation().zy();
2414   rotLoc(3,3) = refSurface.rotation().zz();
2415   */
2416 
2417   CLHEP::HepVector Rt(3), Pt(3), Rm(3), Pm(3), Norm(3);
2418   Rt(1) = GRt.x();
2419   Rt(2) = GRt.y();
2420   Rt(3) = GRt.z();
2421   Pt(1) = GPt.x();
2422   Pt(2) = GPt.y();
2423   Pt(3) = GPt.z();
2424   Rm(1) = GRm.x();
2425   Rm(2) = GRm.y();
2426   Rm(3) = GRm.z();
2427   Pm(1) = GPm.x();
2428   Pm(2) = GPm.y();
2429   Pm(3) = GPm.z();
2430   Norm(1) = GNorm.x();
2431   Norm(2) = GNorm.y();
2432   Norm(3) = GNorm.z();
2433 
2434   CLHEP::HepVector V(4), Rml(3), Pml(3), Rtl(3), Ptl(3);
2435 
2436   /*
2437   R0(1) = refSurface.position().x();
2438   R0(2) = refSurface.position().y();
2439   R0(3) = refSurface.position().z();
2440   */
2441 
2442   Rml = rotLoc * (Rm - R0);
2443   Rtl = rotLoc * (Rt - R0);
2444   Pml = rotLoc * Pm;
2445   Ptl = rotLoc * Pt;
2446 
2447   V(1) = LPRm(0) - Ptl(1) / Ptl(3);
2448   V(2) = LPRm(1) - Ptl(2) / Ptl(3);
2449   V(3) = LPRm(2) - Rtl(1);
2450   V(4) = LPRm(3) - Rtl(2);
2451 
2452   /*
2453   CLHEP::HepSymMatrix Cov(4,0), W(4,0);
2454   for(int i=1; i<=4; i++)
2455     for(int j=1; j<=i; j++){
2456       Cov(i,j) = (tsosMuon.localError().matrix()
2457           + tsosTrack.localError().matrix())(i,j);
2458       //if(i != j) Cov(i,j) = 0.;
2459       //if((i == j) && ((i==1) || (i==2))) Cov(i,j) = 100.;
2460       //if((i == j) && ((i==3) || (i==4))) Cov(i,j) = 10000.;
2461     }
2462   W = Cov;
2463   */
2464 
2465   CLHEP::HepSymMatrix W = covLoc;
2466 
2467   int ierr;
2468   W.invert(ierr);
2469   if (ierr != 0) {
2470     if (alarm)
2471       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " gradientLocal: inversion of matrix W fail ";
2472     return;
2473   }
2474 
2475   //                               JacobianCartesianToLocal
2476 
2477   //AlgebraicMatrix56 jacobian   // differ from calculation above
2478   //= JacobianCartesianToLocal::JacobianCartesianToLocal
2479   //(refSurface, tsosTrack.localParameters()).jacobian();
2480   //for(int i=1; i<=4; i++) for(int j=1; j<=6; j++){
2481   //int j1 = j - 1; JacToLoc(i,j) = jacobian(i, j1);}
2482 
2483   CLHEP::HepMatrix JacToLoc(4, 6, 0);
2484   for (int i = 1; i <= 2; i++)
2485     for (int j = 1; j <= 3; j++) {
2486       JacToLoc(i, j + 3) = (rotLoc(i, j) - rotLoc(3, j) * Pml(i) / Pml(3)) / Pml(3);
2487       JacToLoc(i + 2, j) = rotLoc(i, j);
2488     }
2489 
2490   //                                    JacobianCorrectionsToCartesian
2491   //double PmN = CLHEP::dot(Pm, Norm);
2492   double PmN = CLHEP_dot(Pm, Norm);
2493 
2494   CLHEP::HepMatrix Jac(6, 6, 0);
2495   for (int i = 1; i <= 3; i++)
2496     for (int j = 1; j <= 3; j++) {
2497       Jac(i, j) = Pm(i) * Norm(j) / PmN;
2498       if (i == j)
2499         Jac(i, j) -= 1.;
2500     }
2501 
2502   //                                            dp/da
2503   Jac(4, 4) = 0.;      // dpx/dax
2504   Jac(5, 4) = -Pm(3);  // dpy/dax
2505   Jac(6, 4) = Pm(2);   // dpz/dax
2506   Jac(4, 5) = Pm(3);   // dpx/day
2507   Jac(5, 5) = 0.;      // dpy/day
2508   Jac(6, 5) = -Pm(1);  // dpz/day
2509   Jac(4, 6) = -Pm(2);  // dpx/daz
2510   Jac(5, 6) = Pm(1);   // dpy/daz
2511   Jac(6, 6) = 0.;      // dpz/daz
2512 
2513   CLHEP::HepVector dsda(3);
2514   dsda(1) = (Norm(2) * Rm(3) - Norm(3) * Rm(2)) / PmN;
2515   dsda(2) = (Norm(3) * Rm(1) - Norm(1) * Rm(3)) / PmN;
2516   dsda(3) = (Norm(1) * Rm(2) - Norm(2) * Rm(1)) / PmN;
2517 
2518   //                                             dr/da
2519   Jac(1, 4) = Pm(1) * dsda(1);           // drx/dax
2520   Jac(2, 4) = -Rm(3) + Pm(2) * dsda(1);  // dry/dax
2521   Jac(3, 4) = Rm(2) + Pm(3) * dsda(1);   // drz/dax
2522 
2523   Jac(1, 5) = Rm(3) + Pm(1) * dsda(2);   // drx/day
2524   Jac(2, 5) = Pm(2) * dsda(2);           // dry/day
2525   Jac(3, 5) = -Rm(1) + Pm(3) * dsda(2);  // drz/day
2526 
2527   Jac(1, 6) = -Rm(2) + Pm(1) * dsda(3);  // drx/daz
2528   Jac(2, 6) = Rm(1) + Pm(2) * dsda(3);   // dry/daz
2529   Jac(3, 6) = Pm(3) * dsda(3);           // drz/daz
2530 
2531   //                                   JacobianCorrectionToLocal
2532   CLHEP::HepMatrix JacCorLoc(4, 6, 0);
2533   JacCorLoc = JacToLoc * Jac;
2534 
2535   //                                   gradient and Hessian
2536   CLHEP::HepMatrix W_Jac(6, 4, 0);
2537   W_Jac = JacCorLoc.T() * W;
2538 
2539   CLHEP::HepVector gradL(6);
2540   gradL = W_Jac * V;
2541 
2542   CLHEP::HepMatrix hessL(6, 6);
2543   hessL = JacCorLoc.T() * W * JacCorLoc;
2544 
2545   GradL += gradL;
2546   HessL += hessL;
2547 
2548   CLHEP::HepVector dLI = CLHEP::solve(HessL, -GradL);
2549   int iErI;
2550   CLHEP::HepMatrix ErrdLI = HessL.inverse(iErI);
2551   if (iErI != 0) {
2552     if (alarm)
2553       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " gradLocal Error inverse  Hess matrix !!!!!!!!!!!";
2554   }
2555 
2556   if (info || debug_) {
2557     edm::LogVerbatim("GlobalTrackerMuonAlignment")
2558         << " dL   " << dLI(1) << " " << dLI(2) << " " << dLI(3) << " " << dLI(4) << " " << dLI(5) << " " << dLI(6);
2559     ;
2560     edm::LogVerbatim("GlobalTrackerMuonAlignment")
2561         << " +-   " << std::sqrt(ErrdLI(1, 1)) << " " << std::sqrt(ErrdLI(2, 2)) << " " << std::sqrt(ErrdLI(3, 3))
2562         << " " << std::sqrt(ErrdLI(4, 4)) << " " << std::sqrt(ErrdLI(5, 5)) << " " << std::sqrt(ErrdLI(6, 6));
2563   }
2564 
2565   if (debug_) {
2566     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug")
2567         << " dV(da3) {" << -JacCorLoc(1, 6) * 0.002 << " " << -JacCorLoc(2, 6) * 0.002 << " "
2568         << -JacCorLoc(3, 6) * 0.002 << " " << -JacCorLoc(4, 6) * 0.002 << "}";
2569     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " JacCLo {" << JacCorLoc(1, 6) << " " << JacCorLoc(2, 6)
2570                                                         << " " << JacCorLoc(3, 6) << " " << JacCorLoc(4, 6) << "}";
2571     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug")
2572         << "Jpx/yx  {" << Jac(4, 6) / Pm(3) << " " << Jac(5, 6) / Pm(3) << " " << Jac(1, 6) << " " << Jac(2, 6) << "}";
2573     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug")
2574         << "Jac(,a3){" << Jac(1, 6) << " " << Jac(2, 6) << " " << Jac(3, 6) << " " << Jac(4, 6) << " " << Jac(5, 6)
2575         << " " << Jac(6, 6);
2576     int i = 5;
2577     if (GNorm.z() > 0.95)
2578       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Ecap1   N " << GNorm;
2579     else if (GNorm.z() < -0.95)
2580       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Ecap2   N " << GNorm;
2581     else
2582       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Barrel  N " << GNorm;
2583     edm::LogVerbatim("GlobalTrackerMuonAlignment")
2584         << " JacCLo(i," << i << ") = {" << JacCorLoc(1, i) << " " << JacCorLoc(2, i) << " " << JacCorLoc(3, i) << " "
2585         << JacCorLoc(4, i) << "}";
2586     edm::LogVerbatim("GlobalTrackerMuonAlignment") << "   rotLoc " << rotLoc;
2587     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " position " << R0;
2588     edm::LogVerbatim("GlobalTrackerMuonAlignment")
2589         << " Pm,l " << Pm.T() << " " << Pml(1) / Pml(3) << " " << Pml(2) / Pml(3);
2590     edm::LogVerbatim("GlobalTrackerMuonAlignment")
2591         << " Pt,l " << Pt.T() << " " << Ptl(1) / Ptl(3) << " " << Ptl(2) / Ptl(3);
2592     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " V  " << V.T();
2593     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Cov \n" << covLoc;
2594     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " W*Cov " << W * covLoc;
2595     //edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " JacCarToLoc = drldrc \n" << JacobianCartesianToLocal::JacobianCartesianToLocal(refSurface, tsosTrack.localParameters()).jacobian();
2596     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " JacToLoc " << JacToLoc;
2597   }
2598 
2599 #ifdef CHECK_OF_JACOBIAN_CARTESIAN_TO_LOCAL
2600   //---------------------- check of derivatives
2601   CLHEP::HepVector V0(4, 0);
2602   V0(1) = Pml(1) / Pml(3) - Ptl(1) / Ptl(3);
2603   V0(2) = Pml(2) / Pml(3) - Ptl(2) / Ptl(3);
2604   V0(3) = Rml(1) - Rtl(1);
2605   V0(4) = Rml(2) - Rtl(2);
2606   int ii = 3;
2607   float delta = 0.01;
2608   CLHEP::HepVector V1(4, 0);
2609   if (ii <= 3) {
2610     Rm(ii) += delta;
2611     Rml = rotLoc * (Rm - R0);
2612   } else {
2613     Pm(ii - 3) += delta;
2614     Pml = rotLoc * Pm;
2615   }
2616   //if(ii <= 3) {Rt(ii) += delta; Rtl = rotLoc * (Rt - R0);}
2617   //else {Pt(ii-3) += delta; Ptl = rotLoc * Pt;}
2618   V1(1) = Pml(1) / Pml(3) - Ptl(1) / Ptl(3);
2619   V1(2) = Pml(2) / Pml(3) - Ptl(2) / Ptl(3);
2620   V1(3) = Rml(1) - Rtl(1);
2621   V1(4) = Rml(2) - Rtl(2);
2622 
2623   if (GNorm.z() > 0.95)
2624     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Ecap1   N " << GNorm;
2625   else if (GNorm.z() < -0.95)
2626     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Ecap2   N " << GNorm;
2627   else
2628     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Barrel  N " << GNorm;
2629   edm::LogVerbatim("GlobalTrackerMuonAlignment")
2630       << " dldc Num (i," << ii << ") " << (V1(1) - V0(1)) / delta << " " << (V1(2) - V0(2)) / delta << " "
2631       << (V1(3) - V0(3)) / delta << " " << (V1(4) - V0(4)) / delta;
2632   edm::LogVerbatim("GlobalTrackerMuonAlignment") << " dldc Ana (i," << ii << ") " << JacToLoc(1, ii) << " "
2633                                                  << JacToLoc(2, ii) << " " << JacToLoc(3, ii) << " " << JacToLoc(4, ii);
2634   float dtxdpx = (rotLoc(1, 1) - rotLoc(3, 1) * Pml(1) / Pml(3)) / Pml(3);
2635   float dtydpx = (rotLoc(2, 1) - rotLoc(3, 2) * Pml(2) / Pml(3)) / Pml(3);
2636   edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " dtx/dpx dty/ " << dtxdpx << " " << dtydpx;
2637 #endif
2638 
2639   return;
2640 }  // end gradientLocal
2641 
2642 // ----  simulate a misalignment of muon system   ------
2643 void GlobalTrackerMuonAlignment::misalignMuon(
2644     GlobalVector& GRm, GlobalVector& GPm, GlobalVector& Nl, GlobalVector& Rt, GlobalVector& Rm, GlobalVector& Pm) {
2645   CLHEP::HepVector d(3, 0), a(3, 0), Norm(3), Rm0(3), Pm0(3), Rm1(3), Pm1(3), Rmi(3), Pmi(3);
2646 
2647   d(1) = 0.0;
2648   d(2) = 0.0;
2649   d(3) = 0.0;
2650   a(1) = 0.0000;
2651   a(2) = 0.0000;
2652   a(3) = 0.0000;  // zero
2653   //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
2654   //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
2655   //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
2656   //d(1)=1.0; d(2)=2.0; d(3)=3.0; a(1)=0.0000; a(2)=0.0000; a(3)=0.0000; // 12300
2657   //d(1)=0.0; d(2)=0.0; d(3)=0.0; a(1)=0.0000; a(2)=0.0000; a(3)=0.0020; // a3=0.002
2658   //d(1)=10.; d(2)=20.; d(3)=30.; a(1)=0.0100; a(2)=0.0200; a(3)=0.0300; // huge
2659   //d(1)=10.; d(2)=20.; d(3)=30.; a(1)=0.2000; a(2)=0.2500; a(3)=0.3000; // huge angles
2660   //d(1)=0.3; d(2)=0.4; d(3)=0.5; a(1)=0.0005; a(2)=0.0006; a(3)=0.0007; // 0345 0567
2661   //d(1)=0.3; d(2)=0.4; d(3)=0.5; a(1)=0.0001; a(2)=0.0002; a(3)=0.0003; // 0345 0123
2662   //d(1)=0.2; d(2)=0.2; d(3)=0.2; a(1)=0.0002; a(2)=0.0002; a(3)=0.0002; // 0111 0111
2663 
2664   MuGlShift = d;
2665   MuGlAngle = a;
2666 
2667   if ((d(1) == 0.) && (d(2) == 0.) && (d(3) == 0.) && (a(1) == 0.) && (a(2) == 0.) && (a(3) == 0.)) {
2668     Rm = GRm;
2669     Pm = GPm;
2670     if (debug_)
2671       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ......   without misalignment ";
2672     return;
2673   }
2674 
2675   Rm0(1) = GRm.x();
2676   Rm0(2) = GRm.y();
2677   Rm0(3) = GRm.z();
2678   Pm0(1) = GPm.x();
2679   Pm0(2) = GPm.y();
2680   Pm0(3) = GPm.z();
2681   Norm(1) = Nl.x();
2682   Norm(2) = Nl.y();
2683   Norm(3) = Nl.z();
2684 
2685   CLHEP::HepMatrix T(3, 3, 1);
2686 
2687   //T(1,2) =  a(3); T(1,3) = -a(2);
2688   //T(2,1) = -a(3); T(2,3) =  a(1);
2689   //T(3,1) =  a(2); T(3,2) = -a(1);
2690 
2691   double s1 = std::sin(a(1)), c1 = std::cos(a(1));
2692   double s2 = std::sin(a(2)), c2 = std::cos(a(2));
2693   double s3 = std::sin(a(3)), c3 = std::cos(a(3));
2694   T(1, 1) = c2 * c3;
2695   T(1, 2) = c1 * s3 + s1 * s2 * c3;
2696   T(1, 3) = s1 * s3 - c1 * s2 * c3;
2697   T(2, 1) = -c2 * s3;
2698   T(2, 2) = c1 * c3 - s1 * s2 * s3;
2699   T(2, 3) = s1 * c3 + c1 * s2 * s3;
2700   T(3, 1) = s2;
2701   T(3, 2) = -s1 * c2;
2702   T(3, 3) = c1 * c2;
2703 
2704   //int terr;
2705   //CLHEP::HepMatrix Ti = T.inverse(terr);
2706   CLHEP::HepMatrix Ti = T.T();
2707   CLHEP::HepVector di = -Ti * d;
2708 
2709   CLHEP::HepVector ex0(3, 0), ey0(3, 0), ez0(3, 0), ex(3, 0), ey(3, 0), ez(3, 0);
2710   ex0(1) = 1.;
2711   ey0(2) = 1.;
2712   ez0(3) = 1.;
2713   ex = Ti * ex0;
2714   ey = Ti * ey0;
2715   ez = Ti * ez0;
2716 
2717   CLHEP::HepVector TiN = Ti * Norm;
2718   //double si = CLHEP::dot((Ti*Rm0 - Rm0 + di), TiN) / CLHEP::dot(Pm0, TiN);
2719   //Rm1(1) = CLHEP::dot(ex, Rm0 + si*Pm0 - di);
2720   //Rm1(2) = CLHEP::dot(ey, Rm0 + si*Pm0 - di);
2721   //Rm1(3) = CLHEP::dot(ez, Rm0 + si*Pm0 - di);
2722   double si = CLHEP_dot((Ti * Rm0 - Rm0 + di), TiN) / CLHEP_dot(Pm0, TiN);
2723   Rm1(1) = CLHEP_dot(ex, Rm0 + si * Pm0 - di);
2724   Rm1(2) = CLHEP_dot(ey, Rm0 + si * Pm0 - di);
2725   Rm1(3) = CLHEP_dot(ez, Rm0 + si * Pm0 - di);
2726   Pm1 = T * Pm0;
2727 
2728   Rm = GlobalVector(Rm1(1), Rm1(2), Rm1(3));
2729   //Pm = GlobalVector(CLHEP::dot(Pm0,ex), CLHEP::dot(Pm0, ey), CLHEP::dot(Pm0,ez));// =T*Pm0
2730   Pm = GlobalVector(CLHEP_dot(Pm0, ex), CLHEP_dot(Pm0, ey), CLHEP_dot(Pm0, ez));  // =T*Pm0
2731 
2732   if (debug_) {  //    -------------  debug tranformation
2733 
2734     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ----- Pm " << Pm;
2735     edm::LogVerbatim("GlobalTrackerMuonAlignment") << "    T*Pm0 " << Pm1.T();
2736     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " Ti*T*Pm0 " << (Ti * (T * Pm0)).T();
2737 
2738     //CLHEP::HepVector Rml = (Rm0 + si*Pm0 - di) + di;
2739     CLHEP::HepVector Rml = Rm1(1) * ex + Rm1(2) * ey + Rm1(3) * ez + di;
2740 
2741     CLHEP::HepVector Rt0(3);
2742     Rt0(1) = Rt.x();
2743     Rt0(2) = Rt.y();
2744     Rt0(3) = Rt.z();
2745 
2746     //double sl = CLHEP::dot(T*(Rt0 - Rml), T*Norm) / CLHEP::dot(Ti * Pm1, Norm);
2747     double sl = CLHEP_dot(T * (Rt0 - Rml), T * Norm) / CLHEP_dot(Ti * Pm1, Norm);
2748     CLHEP::HepVector Rl = Rml + sl * (Ti * Pm1);
2749 
2750     //double A = CLHEP::dot(Ti*Pm1, Norm);
2751     //double B = CLHEP::dot(Rt0, Norm) - CLHEP::dot(Ti*Rm1, Norm)
2752     //+ CLHEP::dot(Ti*d, Norm);
2753     double A = CLHEP_dot(Ti * Pm1, Norm);
2754     double B = CLHEP_dot(Rt0, Norm) - CLHEP_dot(Ti * Rm1, Norm) + CLHEP_dot(Ti * d, Norm);
2755     double s = B / A;
2756     CLHEP::HepVector Dr = Ti * Rm1 - Ti * d + s * (Ti * Pm1) - Rm0;
2757     CLHEP::HepVector Dp = Ti * Pm1 - Pm0;
2758 
2759     CLHEP::HepVector TiR = Ti * Rm0 + di;
2760     CLHEP::HepVector Ri = T * TiR + d;
2761 
2762     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " --TTi-0 " << (Ri - Rm0).T();
2763     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " --  Dr  " << Dr.T();
2764     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " --  Dp  " << Dp.T();
2765     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ex " << ex.T();
2766     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ey " << ey.T();
2767     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ez " << ez.T();
2768     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ---- T  ---- " << T;
2769     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ---- Ti ---- " << Ti;
2770     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ---- d  ---- " << d.T();
2771     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ---- di ---- " << di.T();
2772     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -- si sl s " << si << " " << sl << " " << s;
2773     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- si sl " << si << " " << sl;
2774     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- si s " << si << " " << s;
2775     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -- Rml-(Rm0+si*Pm0) " << (Rml - Rm0 - si * Pm0).T();
2776     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- Rm0 " << Rm0.T();
2777     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- Rm1 " << Rm1.T();
2778     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- Rmi " << Rmi.T();
2779     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --siPm " << (si * Pm0).T();
2780     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --s2Pm " << (s2 * (T * Pm1)).T();
2781     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --R1-0 " << (Rm1 - Rm0).T();
2782     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --Ri-0 " << (Rmi - Rm0).T();
2783     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " --Rl-0 " << (Rl - Rm0).T();
2784     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --Pi-0 " << (Pmi - Pm0).T();
2785   }  //                                --------   end of debug
2786 
2787   return;
2788 
2789 }  // end of misalignMuon
2790 
2791 // ----  simulate a misalignment of muon system with Local  ------
2792 void GlobalTrackerMuonAlignment::misalignMuonL(GlobalVector& GRm,
2793                                                GlobalVector& GPm,
2794                                                GlobalVector& Nl,
2795                                                GlobalVector& Rt,
2796                                                GlobalVector& Rm,
2797                                                GlobalVector& Pm,
2798                                                AlgebraicVector4& Vm,
2799                                                TrajectoryStateOnSurface& tsosTrack,
2800                                                TrajectoryStateOnSurface& tsosMuon) {
2801   CLHEP::HepVector d(3, 0), a(3, 0), Norm(3), Rm0(3), Pm0(3), Rm1(3), Pm1(3), Rmi(3), Pmi(3);
2802 
2803   d(1) = 0.0;
2804   d(2) = 0.0;
2805   d(3) = 0.0;
2806   a(1) = 0.0000;
2807   a(2) = 0.0000;
2808   a(3) = 0.0000;  // zero
2809   //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
2810   //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
2811   //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
2812   //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
2813   //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
2814   //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
2815   //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
2816   //d(1)=10.; d(2)=20.; d(3)=30.; a(1)=0.0100; a(2)=0.0200; a(3)=0.0300; // huge
2817   //d(1)=10.; d(2)=20.; d(3)=30.; a(1)=0.2000; a(2)=0.2500; a(3)=0.3000; // huge angles
2818   //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
2819   //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
2820   //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
2821 
2822   MuGlShift = d;
2823   MuGlAngle = a;
2824 
2825   if ((d(1) == 0.) && (d(2) == 0.) && (d(3) == 0.) && (a(1) == 0.) && (a(2) == 0.) && (a(3) == 0.)) {
2826     Rm = GRm;
2827     Pm = GPm;
2828     AlgebraicVector4 Vm0;
2829     Vm = AlgebraicVector4(tsosMuon.localParameters().vector()(1),
2830                           tsosMuon.localParameters().vector()(2),
2831                           tsosMuon.localParameters().vector()(3),
2832                           tsosMuon.localParameters().vector()(4));
2833     if (debug_)
2834       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ......   without misalignment ";
2835     return;
2836   }
2837 
2838   Rm0(1) = GRm.x();
2839   Rm0(2) = GRm.y();
2840   Rm0(3) = GRm.z();
2841   Pm0(1) = GPm.x();
2842   Pm0(2) = GPm.y();
2843   Pm0(3) = GPm.z();
2844   Norm(1) = Nl.x();
2845   Norm(2) = Nl.y();
2846   Norm(3) = Nl.z();
2847 
2848   CLHEP::HepMatrix T(3, 3, 1);
2849 
2850   //T(1,2) =  a(3); T(1,3) = -a(2);
2851   //T(2,1) = -a(3); T(2,3) =  a(1);
2852   //T(3,1) =  a(2); T(3,2) = -a(1);
2853 
2854   double s1 = std::sin(a(1)), c1 = std::cos(a(1));
2855   double s2 = std::sin(a(2)), c2 = std::cos(a(2));
2856   double s3 = std::sin(a(3)), c3 = std::cos(a(3));
2857   T(1, 1) = c2 * c3;
2858   T(1, 2) = c1 * s3 + s1 * s2 * c3;
2859   T(1, 3) = s1 * s3 - c1 * s2 * c3;
2860   T(2, 1) = -c2 * s3;
2861   T(2, 2) = c1 * c3 - s1 * s2 * s3;
2862   T(2, 3) = s1 * c3 + c1 * s2 * s3;
2863   T(3, 1) = s2;
2864   T(3, 2) = -s1 * c2;
2865   T(3, 3) = c1 * c2;
2866 
2867   //                    Ti, di what we have to apply for misalignment
2868   //int terr;
2869   //CLHEP::HepMatrix Ti = T.inverse(terr);
2870   CLHEP::HepMatrix Ti = T.T();
2871   CLHEP::HepVector di = -Ti * d;
2872 
2873   CLHEP::HepVector ex0(3, 0), ey0(3, 0), ez0(3, 0), ex(3, 0), ey(3, 0), ez(3, 0);
2874   ex0(1) = 1.;
2875   ey0(2) = 1.;
2876   ez0(3) = 1.;
2877   ex = Ti * ex0;
2878   ey = Ti * ey0;
2879   ez = Ti * ez0;
2880 
2881   CLHEP::HepVector TiN = Ti * Norm;
2882   //double si = CLHEP::dot((Ti*Rm0 - Rm0 + di), TiN) / CLHEP::dot(Pm0, TiN);
2883   //Rm1(1) = CLHEP::dot(ex, Rm0 + si*Pm0 - di);
2884   //Rm1(2) = CLHEP::dot(ey, Rm0 + si*Pm0 - di);
2885   //Rm1(3) = CLHEP::dot(ez, Rm0 + si*Pm0 - di);
2886   double si = CLHEP_dot((Ti * Rm0 - Rm0 + di), TiN) / CLHEP_dot(Pm0, TiN);
2887   Rm1(1) = CLHEP_dot(ex, Rm0 + si * Pm0 - di);
2888   Rm1(2) = CLHEP_dot(ey, Rm0 + si * Pm0 - di);
2889   Rm1(3) = CLHEP_dot(ez, Rm0 + si * Pm0 - di);
2890   Pm1 = T * Pm0;
2891 
2892   // Global muon parameters after misalignment
2893   Rm = GlobalVector(Rm1(1), Rm1(2), Rm1(3));
2894   //Pm = GlobalVector(CLHEP::dot(Pm0,ex), CLHEP::dot(Pm0, ey), CLHEP::dot(Pm0,ez));// =T*Pm0
2895   Pm = GlobalVector(CLHEP_dot(Pm0, ex), CLHEP_dot(Pm0, ey), CLHEP_dot(Pm0, ez));  // =T*Pm0
2896 
2897   // Local muon parameters after misalignment
2898   const Surface& refSurface = tsosMuon.surface();
2899   CLHEP::HepMatrix Tl(3, 3, 0);
2900 
2901   Tl(1, 1) = refSurface.rotation().xx();
2902   Tl(1, 2) = refSurface.rotation().xy();
2903   Tl(1, 3) = refSurface.rotation().xz();
2904 
2905   Tl(2, 1) = refSurface.rotation().yx();
2906   Tl(2, 2) = refSurface.rotation().yy();
2907   Tl(2, 3) = refSurface.rotation().yz();
2908 
2909   Tl(3, 1) = refSurface.rotation().zx();
2910   Tl(3, 2) = refSurface.rotation().zy();
2911   Tl(3, 3) = refSurface.rotation().zz();
2912 
2913   CLHEP::HepVector R0(3, 0), newR0(3, 0), newRl(3, 0), newPl(3, 0);
2914   R0(1) = refSurface.position().x();
2915   R0(2) = refSurface.position().y();
2916   R0(3) = refSurface.position().z();
2917 
2918   newRl = Tl * (Rm1 - R0);
2919   newPl = Tl * Pm1;
2920 
2921   Vm(0) = newPl(1) / newPl(3);
2922   Vm(1) = newPl(2) / newPl(3);
2923   Vm(2) = newRl(1);
2924   Vm(3) = newRl(2);
2925 
2926 #ifdef CHECK_DERIVATIVES_LOCAL_VS_ANGLES
2927   float del = 0.0001;
2928   //int ii = 1; T(3,2) +=del; T(2,3) -=del;
2929   int ii = 2;
2930   T(3, 1) -= del;
2931   T(1, 3) += del;
2932   //int ii = 3; T(1,2) -=del; T(2,1) +=del;
2933   AlgebraicVector4 Vm0 = Vm;
2934   CLHEP::HepVector Rm10 = Rm1, Pm10 = Pm1;
2935 
2936   CLHEP::HepMatrix newTl = Tl * T;
2937   Ti = T.T();
2938   di = -Ti * d;
2939   ex = Ti * ex0;
2940   ey = Ti * ey0;
2941   ez = Ti * ez0;
2942   TiN = Ti * Norm;
2943   si = CLHEP_dot((Ti * Rm0 - Rm0 + di), TiN) / CLHEP_dot(Pm0, TiN);
2944   Rm1(1) = CLHEP_dot(ex, Rm0 + si * Pm0 - di);
2945   Rm1(2) = CLHEP_dot(ey, Rm0 + si * Pm0 - di);
2946   Rm1(3) = CLHEP_dot(ez, Rm0 + si * Pm0 - di);
2947   Pm1 = T * Pm0;
2948 
2949   newRl = Tl * (Rm1 - R0);
2950   newPl = Tl * Pm1;
2951 
2952   Vm(0) = newPl(1) / newPl(3);
2953   Vm(1) = newPl(2) / newPl(3);
2954   Vm(2) = newRl(1);
2955   Vm(3) = newRl(2);
2956   edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ========= dVm/da" << ii << " " << (Vm - Vm0) * (1. / del);
2957   if (debug_)
2958     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug")
2959         << " dRm/da3 " << ((Rm1 - Rm10) * (1. / del)).T() << " " << ((Pm1 - Pm10) * (1. / del)).T();
2960 #endif
2961 
2962   if (debug_) {
2963     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << "             Norm " << Norm.T();
2964     edm::LogVerbatim("GlobalTrackerMuonAlignment") << "             ex  " << (Tl.T() * ex0).T();
2965     edm::LogVerbatim("GlobalTrackerMuonAlignment") << "             ey  " << (Tl.T() * ey0).T();
2966     edm::LogVerbatim("GlobalTrackerMuonAlignment") << "             ez  " << (Tl.T() * ez0).T();
2967 
2968     edm::LogVerbatim("GlobalTrackerMuonAlignment") << "    pxyz/p gl 0 " << (Pm0 / Pm0.norm()).T();
2969     edm::LogVerbatim("GlobalTrackerMuonAlignment") << "    pxyz/p loc0 " << (Tl * Pm0 / (Tl * Pm0).norm()).T();
2970     edm::LogVerbatim("GlobalTrackerMuonAlignment") << "    pxyz/p glob  " << (Pm1 / Pm1.norm()).T();
2971     edm::LogVerbatim("GlobalTrackerMuonAlignment") << "    pxyz/p  loc  " << (newPl / newPl.norm()).T();
2972 
2973     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << "             Rot   " << refSurface.rotation();
2974     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << "             Tl   " << Tl;
2975     edm::LogVerbatim("GlobalTrackerMuonAlignment")
2976         << " ---- phi g0 g1 l0 l1  " << atan2(Pm0(2), Pm0(1)) << " " << atan2(Pm1(2), Pm1(1)) << " "
2977         << atan2((Tl * Pm0)(2), (Tl * Pm0)(1)) << " " << atan2(newPl(2), newPl(1));
2978     edm::LogVerbatim("GlobalTrackerMuonAlignment")
2979         << " ---- angl Gl Loc " << atan2(Pm1(2), Pm1(1)) - atan2(Pm0(2), Pm0(1)) << " "
2980         << atan2(newPl(2), newPl(1)) - atan2((Tl * Pm0)(2), (Tl * Pm0)(1)) << " "
2981         << atan2(newPl(3), newPl(2)) - atan2((Tl * Pm0)(3), (Tl * Pm0)(2)) << " "
2982         << atan2(newPl(1), newPl(3)) - atan2((Tl * Pm0)(1), (Tl * Pm0)(3)) << " ";
2983   }
2984 
2985   if (debug_) {
2986     CLHEP::HepMatrix newTl(3, 3, 0);
2987     CLHEP::HepVector R0(3, 0), newRl(3, 0), newPl(3, 0);
2988     newTl = Tl * Ti.T();
2989     newR0 = Ti * R0 + di;
2990 
2991     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " N " << Norm.T();
2992     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " dtxl yl " << Vm(0) - tsosMuon.localParameters().vector()(1)
2993                                                    << " " << Vm(1) - tsosMuon.localParameters().vector()(2);
2994     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " dXl dYl " << Vm(2) - tsosMuon.localParameters().vector()(3)
2995                                                    << " " << Vm(3) - tsosMuon.localParameters().vector()(4);
2996     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " localM    " << tsosMuon.localParameters().vector();
2997     edm::LogVerbatim("GlobalTrackerMuonAlignment") << "       Vm  " << Vm;
2998 
2999     CLHEP::HepVector Rvc(3, 0), Pvc(3, 0), Rvg(3, 0), Pvg(3, 0);
3000     Rvc(1) = Vm(2);
3001     Rvc(2) = Vm(3);
3002     Pvc(3) = tsosMuon.localParameters().pzSign() * Pm0.norm() / std::sqrt(1 + Vm(0) * Vm(0) + Vm(1) * Vm(1));
3003     Pvc(1) = Pvc(3) * Vm(0);
3004     Pvc(2) = Pvc(3) * Vm(1);
3005 
3006     Rvg = newTl.T() * Rvc + newR0;
3007     Pvg = newTl.T() * Pvc;
3008 
3009     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " newPl     " << newPl.T();
3010     edm::LogVerbatim("GlobalTrackerMuonAlignment") << "    Pvc    " << Pvc.T();
3011     edm::LogVerbatim("GlobalTrackerMuonAlignment") << "    Rvg    " << Rvg.T();
3012     edm::LogVerbatim("GlobalTrackerMuonAlignment") << "    Rm1    " << Rm1.T();
3013     edm::LogVerbatim("GlobalTrackerMuonAlignment") << "    Pvg    " << Pvg.T();
3014     edm::LogVerbatim("GlobalTrackerMuonAlignment") << "    Pm1    " << Pm1.T();
3015   }
3016 
3017   if (debug_) {  //    ----------  how to transform cartesian from shifted to correct
3018 
3019     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ----- Pm " << Pm;
3020     edm::LogVerbatim("GlobalTrackerMuonAlignment") << "    T*Pm0 " << Pm1.T();
3021     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " Ti*T*Pm0 " << (Ti * (T * Pm0)).T();
3022 
3023     //CLHEP::HepVector Rml = (Rm0 + si*Pm0 - di) + di;
3024     CLHEP::HepVector Rml = Rm1(1) * ex + Rm1(2) * ey + Rm1(3) * ez + di;
3025 
3026     CLHEP::HepVector Rt0(3);
3027     Rt0(1) = Rt.x();
3028     Rt0(2) = Rt.y();
3029     Rt0(3) = Rt.z();
3030 
3031     //double sl = CLHEP::dot(T*(Rt0 - Rml), T*Norm) / CLHEP::dot(Ti * Pm1, Norm);
3032     double sl = CLHEP_dot(T * (Rt0 - Rml), T * Norm) / CLHEP_dot(Ti * Pm1, Norm);
3033     CLHEP::HepVector Rl = Rml + sl * (Ti * Pm1);
3034 
3035     //double A = CLHEP::dot(Ti*Pm1, Norm);
3036     //double B = CLHEP::dot(Rt0, Norm) - CLHEP::dot(Ti*Rm1, Norm)
3037     //+ CLHEP::dot(Ti*d, Norm);
3038     double A = CLHEP_dot(Ti * Pm1, Norm);
3039     double B = CLHEP_dot(Rt0, Norm) - CLHEP_dot(Ti * Rm1, Norm) + CLHEP_dot(Ti * d, Norm);
3040     double s = B / A;
3041     CLHEP::HepVector Dr = Ti * Rm1 - Ti * d + s * (Ti * Pm1) - Rm0;
3042     CLHEP::HepVector Dp = Ti * Pm1 - Pm0;
3043 
3044     CLHEP::HepVector TiR = Ti * Rm0 + di;
3045     CLHEP::HepVector Ri = T * TiR + d;
3046 
3047     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " --TTi-0 " << (Ri - Rm0).T();
3048     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " --  Dr  " << Dr.T();
3049     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " --  Dp  " << Dp.T();
3050     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ex " << ex.T();
3051     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ey " << ey.T();
3052     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ez " << ez.T();
3053     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ---- T  ---- " << T;
3054     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ---- Ti ---- " << Ti;
3055     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ---- d  ---- " << d.T();
3056     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " ---- di ---- " << di.T();
3057     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -- si sl s " << si << " " << sl << " " << s;
3058     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- si sl " << si << " " << sl;
3059     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- si s " << si << " " << s;
3060     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " -- Rml-(Rm0+si*Pm0) " << (Rml - Rm0 - si * Pm0).T();
3061     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- Rm0 " << Rm0.T();
3062     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- Rm1 " << Rm1.T();
3063     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " -- Rmi " << Rmi.T();
3064     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --siPm " << (si * Pm0).T();
3065     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --s2Pm " << (s2 * (T * Pm1)).T();
3066     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --R1-0 " << (Rm1 - Rm0).T();
3067     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --Ri-0 " << (Rmi - Rm0).T();
3068     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " --Rl-0 " << (Rl - Rm0).T();
3069     edm::LogVerbatim("GlobalTrackerMuonAlignmentDebug") << " --Pi-0 " << (Pmi - Pm0).T();
3070   }  //                                --------   end of debug
3071 
3072   return;
3073 
3074 }  // end misalignMuonL
3075 
3076 // ----  refit any direction of transient track ------
3077 void GlobalTrackerMuonAlignment::trackFitter(reco::TrackRef alongTr,
3078                                              reco::TransientTrack& alongTTr,
3079                                              PropagationDirection direction,
3080                                              TrajectoryStateOnSurface& trackFittedTSOS) {
3081   bool info = false;
3082   bool smooth = false;
3083 
3084   if (info)
3085     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ......... start of trackFitter ......... ";
3086   if (false && info) {
3087     this->debugTrackHit(" Tracker track hits ", alongTr);
3088     this->debugTrackHit(" Tracker TransientTrack hits ", alongTTr);
3089   }
3090 
3091   edm::OwnVector<TrackingRecHit> recHit;
3092   DetId trackDetId = DetId(alongTr->innerDetId());
3093   for (auto const& hit : alongTr->recHits())
3094     recHit.push_back(hit->clone());
3095   if (direction != alongMomentum) {
3096     edm::OwnVector<TrackingRecHit> recHitAlong = recHit;
3097     recHit.clear();
3098     for (unsigned int ihit = recHitAlong.size() - 1; ihit + 1 > 0; ihit--) {
3099       recHit.push_back(recHitAlong[ihit]);
3100     }
3101     recHitAlong.clear();
3102   }
3103   if (info)
3104     edm::LogVerbatim("GlobalTrackerMuonAlignment")
3105         << " ... Own recHit.size() " << recHit.size() << " " << trackDetId.rawId();
3106 
3107   //MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitTrack;
3108   TransientTrackingRecHit::RecHitContainer recHitTrack;
3109   for (auto const& hit : alongTr->recHits())
3110     recHitTrack.push_back(TTRHBuilder->build(hit));
3111 
3112   if (info)
3113     edm::LogVerbatim("GlobalTrackerMuonAlignment")
3114         << " ... recHitTrack.size() " << recHit.size() << " " << trackDetId.rawId() << " ======> recHitMu ";
3115 
3116   MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu = recHitTrack;
3117   /*
3118     MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu =
3119     MuRHBuilder->build(alongTTr.recHits().begin(), alongTTr.recHits().end());
3120   */
3121   if (info)
3122     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ...along.... recHitMu.size() " << recHitMu.size();
3123   if (direction != alongMomentum) {
3124     MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMuAlong = recHitMu;
3125     recHitMu.clear();
3126     for (unsigned int ihit = recHitMuAlong.size() - 1; ihit + 1 > 0; ihit--) {
3127       recHitMu.push_back(recHitMuAlong[ihit]);
3128     }
3129     recHitMuAlong.clear();
3130   }
3131   if (info)
3132     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ...opposite... recHitMu.size() " << recHitMu.size();
3133 
3134   TrajectoryStateOnSurface firstTSOS;
3135   if (direction == alongMomentum)
3136     firstTSOS = alongTTr.innermostMeasurementState();
3137   else
3138     firstTSOS = alongTTr.outermostMeasurementState();
3139 
3140   AlgebraicSymMatrix55 CovLoc;
3141   CovLoc(0, 0) = 1. * firstTSOS.localError().matrix()(0, 0);
3142   CovLoc(1, 1) = 10. * firstTSOS.localError().matrix()(1, 1);
3143   CovLoc(2, 2) = 10. * firstTSOS.localError().matrix()(2, 2);
3144   CovLoc(3, 3) = 100. * firstTSOS.localError().matrix()(3, 3);
3145   CovLoc(4, 4) = 100. * firstTSOS.localError().matrix()(4, 4);
3146   TrajectoryStateOnSurface initialTSOS(
3147       firstTSOS.localParameters(), LocalTrajectoryError(CovLoc), firstTSOS.surface(), &*magneticField_);
3148 
3149   PTrajectoryStateOnDet PTraj = trajectoryStateTransform::persistentState(initialTSOS, trackDetId.rawId());
3150   //const TrajectorySeed seedT(*PTraj, recHit, alongMomentum);
3151   const TrajectorySeed seedT(PTraj, recHit, direction);
3152 
3153   std::vector<Trajectory> trajVec;
3154   if (direction == alongMomentum)
3155     trajVec = theFitter->fit(seedT, recHitMu, initialTSOS);
3156   else
3157     trajVec = theFitterOp->fit(seedT, recHitMu, initialTSOS);
3158 
3159   if (info) {
3160     this->debugTrajectorySOSv(" innermostTSOS of TransientTrack", alongTTr.innermostMeasurementState());
3161     this->debugTrajectorySOSv(" outermostTSOS of TransientTrack", alongTTr.outermostMeasurementState());
3162     this->debugTrajectorySOS(" initialTSOS for theFitter ", initialTSOS);
3163     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " . . . . . trajVec.size() " << trajVec.size();
3164     if (!trajVec.empty())
3165       this->debugTrajectory(" theFitter trajectory ", trajVec[0]);
3166   }
3167 
3168   if (!smooth) {
3169     if (!trajVec.empty())
3170       trackFittedTSOS = trajVec[0].lastMeasurement().updatedState();
3171   } else {
3172     std::vector<Trajectory> trajSVec;
3173     if (!trajVec.empty()) {
3174       if (direction == alongMomentum)
3175         trajSVec = theSmoother->trajectories(*(trajVec.begin()));
3176       else
3177         trajSVec = theSmootherOp->trajectories(*(trajVec.begin()));
3178     }
3179     if (info)
3180       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " . . . . trajSVec.size() " << trajSVec.size();
3181     if (!trajSVec.empty())
3182       this->debugTrajectorySOSv("smoothed geom InnermostState", trajSVec[0].geometricalInnermostState());
3183     if (!trajSVec.empty())
3184       trackFittedTSOS = trajSVec[0].geometricalInnermostState();
3185   }
3186 
3187   if (info)
3188     this->debugTrajectorySOSv(" trackFittedTSOS ", trackFittedTSOS);
3189 
3190 }  // end of trackFitter
3191 
3192 // ----  refit any direction of muon transient track ------
3193 void GlobalTrackerMuonAlignment::muonFitter(reco::TrackRef alongTr,
3194                                             reco::TransientTrack& alongTTr,
3195                                             PropagationDirection direction,
3196                                             TrajectoryStateOnSurface& trackFittedTSOS) {
3197   bool info = false;
3198   bool smooth = false;
3199 
3200   if (info)
3201     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ......... start of muonFitter ........";
3202   if (false && info) {
3203     this->debugTrackHit(" Muon track hits ", alongTr);
3204     this->debugTrackHit(" Muon TransientTrack hits ", alongTTr);
3205   }
3206 
3207   edm::OwnVector<TrackingRecHit> recHit;
3208   DetId trackDetId = DetId(alongTr->innerDetId());
3209   for (auto const& hit : alongTr->recHits())
3210     recHit.push_back(hit->clone());
3211   if (direction != alongMomentum) {
3212     edm::OwnVector<TrackingRecHit> recHitAlong = recHit;
3213     recHit.clear();
3214     for (unsigned int ihit = recHitAlong.size() - 1; ihit + 1 > 0; ihit--) {
3215       recHit.push_back(recHitAlong[ihit]);
3216     }
3217     recHitAlong.clear();
3218   }
3219   if (info)
3220     edm::LogVerbatim("GlobalTrackerMuonAlignment")
3221         << " ... Own recHit.size() DetId==Muon " << recHit.size() << " " << trackDetId.rawId();
3222 
3223   MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu =
3224       MuRHBuilder->build(alongTTr.recHitsBegin(), alongTTr.recHitsEnd());
3225   if (info)
3226     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ...along.... recHitMu.size() " << recHitMu.size();
3227   if (direction != alongMomentum) {
3228     MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMuAlong = recHitMu;
3229     recHitMu.clear();
3230     for (unsigned int ihit = recHitMuAlong.size() - 1; ihit + 1 > 0; ihit--) {
3231       recHitMu.push_back(recHitMuAlong[ihit]);
3232     }
3233     recHitMuAlong.clear();
3234   }
3235   if (info)
3236     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ...opposite... recHitMu.size() " << recHitMu.size();
3237 
3238   TrajectoryStateOnSurface firstTSOS;
3239   if (direction == alongMomentum)
3240     firstTSOS = alongTTr.innermostMeasurementState();
3241   else
3242     firstTSOS = alongTTr.outermostMeasurementState();
3243 
3244   AlgebraicSymMatrix55 CovLoc;
3245   CovLoc(0, 0) = 1. * firstTSOS.localError().matrix()(0, 0);
3246   CovLoc(1, 1) = 10. * firstTSOS.localError().matrix()(1, 1);
3247   CovLoc(2, 2) = 10. * firstTSOS.localError().matrix()(2, 2);
3248   CovLoc(3, 3) = 100. * firstTSOS.localError().matrix()(3, 3);
3249   CovLoc(4, 4) = 100. * firstTSOS.localError().matrix()(4, 4);
3250   TrajectoryStateOnSurface initialTSOS(
3251       firstTSOS.localParameters(), LocalTrajectoryError(CovLoc), firstTSOS.surface(), &*magneticField_);
3252 
3253   PTrajectoryStateOnDet PTraj = trajectoryStateTransform::persistentState(initialTSOS, trackDetId.rawId());
3254   //const TrajectorySeed seedT(*PTraj, recHit, alongMomentum);
3255   const TrajectorySeed seedT(PTraj, recHit, direction);
3256 
3257   std::vector<Trajectory> trajVec;
3258   if (direction == alongMomentum)
3259     trajVec = theFitter->fit(seedT, recHitMu, initialTSOS);
3260   else
3261     trajVec = theFitterOp->fit(seedT, recHitMu, initialTSOS);
3262 
3263   if (info) {
3264     this->debugTrajectorySOSv(" innermostTSOS of TransientTrack", alongTTr.innermostMeasurementState());
3265     this->debugTrajectorySOSv(" outermostTSOS of TransientTrack", alongTTr.outermostMeasurementState());
3266     this->debugTrajectorySOS(" initialTSOS for theFitter ", initialTSOS);
3267     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " . . . . . trajVec.size() " << trajVec.size();
3268     if (!trajVec.empty())
3269       this->debugTrajectory(" theFitter trajectory ", trajVec[0]);
3270   }
3271 
3272   if (!smooth) {
3273     if (!trajVec.empty())
3274       trackFittedTSOS = trajVec[0].lastMeasurement().updatedState();
3275   } else {
3276     std::vector<Trajectory> trajSVec;
3277     if (!trajVec.empty()) {
3278       if (direction == alongMomentum)
3279         trajSVec = theSmoother->trajectories(*(trajVec.begin()));
3280       else
3281         trajSVec = theSmootherOp->trajectories(*(trajVec.begin()));
3282     }
3283     if (info)
3284       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " . . . . trajSVec.size() " << trajSVec.size();
3285     if (!trajSVec.empty())
3286       this->debugTrajectorySOSv("smoothed geom InnermostState", trajSVec[0].geometricalInnermostState());
3287     if (!trajSVec.empty())
3288       trackFittedTSOS = trajSVec[0].geometricalInnermostState();
3289   }
3290 
3291 }  // end of muonFitter
3292 
3293 // ----  debug printout of hits from TransientTrack  ------
3294 void GlobalTrackerMuonAlignment::debugTrackHit(const std::string title, TransientTrack& alongTr) {
3295   edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ------- " << title << " --------";
3296   int nHit = 1;
3297   for (trackingRecHit_iterator i = alongTr.recHitsBegin(); i != alongTr.recHitsEnd(); i++) {
3298     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Hit " << nHit++ << " DetId " << (*i)->geographicalId().det();
3299     if ((*i)->geographicalId().det() == DetId::Tracker)
3300       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Tracker ";
3301     else if ((*i)->geographicalId().det() == DetId::Muon)
3302       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Muon ";
3303     else
3304       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Unknown ";
3305     if (!(*i)->isValid())
3306       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " not valid ";
3307     else
3308       edm::LogVerbatim("GlobalTrackerMuonAlignment");
3309   }
3310 }
3311 
3312 // ----  debug printout of hits from TrackRef   ------
3313 void GlobalTrackerMuonAlignment::debugTrackHit(const std::string title, reco::TrackRef alongTr) {
3314   edm::LogVerbatim("GlobalTrackerMuonAlignment") << " ------- " << title << " --------";
3315   int nHit = 1;
3316   for (auto const& hit : alongTr->recHits()) {
3317     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Hit " << nHit++ << " DetId " << hit->geographicalId().det();
3318     if (hit->geographicalId().det() == DetId::Tracker)
3319       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Tracker ";
3320     else if (hit->geographicalId().det() == DetId::Muon)
3321       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Muon ";
3322     else
3323       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Unknown ";
3324     if (!hit->isValid())
3325       edm::LogVerbatim("GlobalTrackerMuonAlignment") << " not valid ";
3326     else
3327       edm::LogVerbatim("GlobalTrackerMuonAlignment");
3328   }
3329 }
3330 
3331 // ----  debug printout TrajectoryStateOnSurface  ------
3332 void GlobalTrackerMuonAlignment::debugTrajectorySOS(const std::string title, TrajectoryStateOnSurface& trajSOS) {
3333   edm::LogVerbatim("GlobalTrackerMuonAlignment") << "    --- " << title << " --- ";
3334   if (!trajSOS.isValid()) {
3335     edm::LogVerbatim("GlobalTrackerMuonAlignment") << "      Not valid !!!! ";
3336     return;
3337   }
3338   edm::LogVerbatim("GlobalTrackerMuonAlignment") << " R |p| " << trajSOS.globalPosition().perp() << " "
3339                                                  << trajSOS.globalMomentum().mag() << " charge " << trajSOS.charge();
3340   edm::LogVerbatim("GlobalTrackerMuonAlignment")
3341       << " x p  " << trajSOS.globalParameters().vector()(0) << " " << trajSOS.globalParameters().vector()(1) << " "
3342       << trajSOS.globalParameters().vector()(2) << " " << trajSOS.globalParameters().vector()(3) << " "
3343       << trajSOS.globalParameters().vector()(4) << " " << trajSOS.globalParameters().vector()(5);
3344   edm::LogVerbatim("GlobalTrackerMuonAlignment")
3345       << " +/- " << std::sqrt(trajSOS.cartesianError().matrix()(0, 0)) << " "
3346       << std::sqrt(trajSOS.cartesianError().matrix()(1, 1)) << " " << std::sqrt(trajSOS.cartesianError().matrix()(2, 2))
3347       << " " << std::sqrt(trajSOS.cartesianError().matrix()(3, 3)) << " "
3348       << std::sqrt(trajSOS.cartesianError().matrix()(4, 4)) << " "
3349       << std::sqrt(trajSOS.cartesianError().matrix()(5, 5));
3350   edm::LogVerbatim("GlobalTrackerMuonAlignment")
3351       << " q/p dxy/dz xy " << trajSOS.localParameters().vector()(0) << " " << trajSOS.localParameters().vector()(1)
3352       << " " << trajSOS.localParameters().vector()(2) << " " << trajSOS.localParameters().vector()(3) << " "
3353       << trajSOS.localParameters().vector()(4);
3354   edm::LogVerbatim("GlobalTrackerMuonAlignment")
3355       << "   +/- error  " << std::sqrt(trajSOS.localError().matrix()(0, 0)) << " "
3356       << std::sqrt(trajSOS.localError().matrix()(1, 1)) << " " << std::sqrt(trajSOS.localError().matrix()(2, 2)) << " "
3357       << std::sqrt(trajSOS.localError().matrix()(3, 3)) << " " << std::sqrt(trajSOS.localError().matrix()(4, 4)) << " ";
3358 }
3359 
3360 // ----  debug printout TrajectoryStateOnSurface  ------
3361 void GlobalTrackerMuonAlignment::debugTrajectorySOSv(const std::string title, TrajectoryStateOnSurface trajSOS) {
3362   edm::LogVerbatim("GlobalTrackerMuonAlignment") << "    --- " << title << " --- ";
3363   if (!trajSOS.isValid()) {
3364     edm::LogVerbatim("GlobalTrackerMuonAlignment") << "      Not valid !!!! ";
3365     return;
3366   }
3367   edm::LogVerbatim("GlobalTrackerMuonAlignment") << " R |p| " << trajSOS.globalPosition().perp() << " "
3368                                                  << trajSOS.globalMomentum().mag() << " charge " << trajSOS.charge();
3369   edm::LogVerbatim("GlobalTrackerMuonAlignment")
3370       << " x p  " << trajSOS.globalParameters().vector()(0) << " " << trajSOS.globalParameters().vector()(1) << " "
3371       << trajSOS.globalParameters().vector()(2) << " " << trajSOS.globalParameters().vector()(3) << " "
3372       << trajSOS.globalParameters().vector()(4) << " " << trajSOS.globalParameters().vector()(5);
3373   edm::LogVerbatim("GlobalTrackerMuonAlignment")
3374       << " +/- " << std::sqrt(trajSOS.cartesianError().matrix()(0, 0)) << " "
3375       << std::sqrt(trajSOS.cartesianError().matrix()(1, 1)) << " " << std::sqrt(trajSOS.cartesianError().matrix()(2, 2))
3376       << " " << std::sqrt(trajSOS.cartesianError().matrix()(3, 3)) << " "
3377       << std::sqrt(trajSOS.cartesianError().matrix()(4, 4)) << " "
3378       << std::sqrt(trajSOS.cartesianError().matrix()(5, 5));
3379   edm::LogVerbatim("GlobalTrackerMuonAlignment")
3380       << " q/p dxy/dz xy " << trajSOS.localParameters().vector()(0) << " " << trajSOS.localParameters().vector()(1)
3381       << " " << trajSOS.localParameters().vector()(2) << " " << trajSOS.localParameters().vector()(3) << " "
3382       << trajSOS.localParameters().vector()(4);
3383   edm::LogVerbatim("GlobalTrackerMuonAlignment")
3384       << "   +/- error  " << std::sqrt(trajSOS.localError().matrix()(0, 0)) << " "
3385       << std::sqrt(trajSOS.localError().matrix()(1, 1)) << " " << std::sqrt(trajSOS.localError().matrix()(2, 2)) << " "
3386       << std::sqrt(trajSOS.localError().matrix()(3, 3)) << " " << std::sqrt(trajSOS.localError().matrix()(4, 4)) << " ";
3387 }
3388 
3389 // ----  debug printout Trajectory   ------
3390 void GlobalTrackerMuonAlignment::debugTrajectory(const std::string title, Trajectory& traj) {
3391   edm::LogVerbatim("GlobalTrackerMuonAlignment") << "\n"
3392                                                  << "    ...... " << title << " ...... ";
3393   if (!traj.isValid()) {
3394     edm::LogVerbatim("GlobalTrackerMuonAlignment") << "          Not valid !!!!!!!!  ";
3395     return;
3396   }
3397   edm::LogVerbatim("GlobalTrackerMuonAlignment") << " chi2/Nhit " << traj.chiSquared() << " / " << traj.foundHits();
3398   if (traj.direction() == alongMomentum)
3399     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " alongMomentum  >>>>";
3400   else
3401     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " oppositeToMomentum  <<<<";
3402   this->debugTrajectorySOSv(" firstMeasurementTSOS ", traj.firstMeasurement().updatedState());
3403   //this->debugTrajectorySOSv(" firstMeasurementTSOS ",traj.firstMeasurement().predictedState());
3404   this->debugTrajectorySOSv(" lastMeasurementTSOS ", traj.lastMeasurement().updatedState());
3405   //this->debugTrajectorySOSv(" geom InnermostState", traj.geometricalInnermostState());
3406   edm::LogVerbatim("GlobalTrackerMuonAlignment") << "    . . . . . . . . . . . . . . . . . . . . . . . . . . . . \n";
3407 }
3408 
3409 // ----  write GlobalPositionRcd   ------
3410 void GlobalTrackerMuonAlignment::writeGlPosRcd(CLHEP::HepVector& paramVec) {
3411   //paramVec(1) = 0.1; paramVec(2) = 0.2; paramVec(3) = 0.3;         //0123
3412   //paramVec(4) = 0.0001; paramVec(5) = 0.0002; paramVec(6) = 0.0003;
3413   //paramVec(1) = 0.3; paramVec(2) = 0.4; paramVec(3) = 0.5;         //03450123
3414   //paramVec(4) = 0.0001; paramVec(5) = 0.0002; paramVec(6) = 0.0003;
3415   //paramVec(1) = 2.; paramVec(2) = 2.; paramVec(3) = 2.;               //222
3416   //paramVec(4) = 0.02; paramVec(5) = 0.02; paramVec(6) = 0.02;
3417   //paramVec(1) = -10.; paramVec(2) = -20.; paramVec(3) = -30.;         //102030
3418   //paramVec(4) = -0.1; paramVec(5) = -0.2; paramVec(6) = -0.3;
3419   //paramVec(1) = 0.; paramVec(2) = 0.; paramVec(3) = 1.;               //1z
3420   //paramVec(4) = 0.0000; paramVec(5) = 0.0000; paramVec(6) = 0.0000;
3421 
3422   edm::LogVerbatim("GlobalTrackerMuonAlignment") << " paramVector " << paramVec.T();
3423 
3424   CLHEP::Hep3Vector colX, colY, colZ;
3425 
3426 #ifdef NOT_EXACT_ROTATION_MATRIX
3427   colX = CLHEP::Hep3Vector(1., -paramVec(6), paramVec(5));
3428   colY = CLHEP::Hep3Vector(paramVec(6), 1., -paramVec(4));
3429   colZ = CLHEP::Hep3Vector(-paramVec(5), paramVec(4), 1.);
3430 #else
3431   double s1 = std::sin(paramVec(4)), c1 = std::cos(paramVec(4));
3432   double s2 = std::sin(paramVec(5)), c2 = std::cos(paramVec(5));
3433   double s3 = std::sin(paramVec(6)), c3 = std::cos(paramVec(6));
3434   colX = CLHEP::Hep3Vector(c2 * c3, -c2 * s3, s2);
3435   colY = CLHEP::Hep3Vector(c1 * s3 + s1 * s2 * c3, c1 * c3 - s1 * s2 * s3, -s1 * c2);
3436   colZ = CLHEP::Hep3Vector(s1 * s3 - c1 * s2 * c3, s1 * c3 + c1 * s2 * s3, c1 * c2);
3437 #endif
3438 
3439   CLHEP::HepVector param0(6, 0);
3440 
3441   Alignments globalPositions{};
3442 
3443   //Tracker
3444   AlignTransform tracker(
3445       iteratorTrackerRcd->translation(), iteratorTrackerRcd->rotation(), DetId(DetId::Tracker).rawId());
3446   // Muon
3447   CLHEP::Hep3Vector posMuGlRcd = iteratorMuonRcd->translation();
3448   CLHEP::HepRotation rotMuGlRcd = iteratorMuonRcd->rotation();
3449   CLHEP::HepEulerAngles angMuGlRcd = iteratorMuonRcd->rotation().eulerAngles();
3450   if (debug_)
3451     edm::LogVerbatim("GlobalTrackerMuonAlignment")
3452         << " Old muon Rcd Euler angles " << angMuGlRcd.phi() << " " << angMuGlRcd.theta() << " " << angMuGlRcd.psi();
3453   AlignTransform muon;
3454   if ((angMuGlRcd.phi() == 0.) && (angMuGlRcd.theta() == 0.) && (angMuGlRcd.psi() == 0.) && (posMuGlRcd.x() == 0.) &&
3455       (posMuGlRcd.y() == 0.) && (posMuGlRcd.z() == 0.)) {
3456     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " New muon parameters are stored in Rcd";
3457 
3458     AlignTransform muonNew(AlignTransform::Translation(paramVec(1), paramVec(2), paramVec(3)),
3459                            AlignTransform::Rotation(colX, colY, colZ),
3460                            DetId(DetId::Muon).rawId());
3461     muon = muonNew;
3462   } else if ((paramVec(1) == 0.) && (paramVec(2) == 0.) && (paramVec(3) == 0.) && (paramVec(4) == 0.) &&
3463              (paramVec(5) == 0.) && (paramVec(6) == 0.)) {
3464     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " Old muon parameters are stored in Rcd";
3465 
3466     AlignTransform muonNew(iteratorMuonRcd->translation(), iteratorMuonRcd->rotation(), DetId(DetId::Muon).rawId());
3467     muon = muonNew;
3468   } else {
3469     edm::LogVerbatim("GlobalTrackerMuonAlignment") << " New + Old muon parameters are stored in Rcd";
3470     CLHEP::Hep3Vector posMuGlRcdThis = CLHEP::Hep3Vector(paramVec(1), paramVec(2), paramVec(3));
3471     CLHEP::HepRotation rotMuGlRcdThis = CLHEP::HepRotation(colX, colY, colZ);
3472     CLHEP::Hep3Vector posMuGlRcdNew = rotMuGlRcdThis * posMuGlRcd + posMuGlRcdThis;
3473     CLHEP::HepRotation rotMuGlRcdNew = rotMuGlRcdThis * rotMuGlRcd;
3474 
3475     AlignTransform muonNew(posMuGlRcdNew, rotMuGlRcdNew, DetId(DetId::Muon).rawId());
3476     muon = muonNew;
3477   }
3478 
3479   // Ecal
3480   AlignTransform ecal(iteratorEcalRcd->translation(), iteratorEcalRcd->rotation(), DetId(DetId::Ecal).rawId());
3481   // Hcal
3482   AlignTransform hcal(iteratorHcalRcd->translation(), iteratorHcalRcd->rotation(), DetId(DetId::Hcal).rawId());
3483   // Calo
3484   AlignTransform calo(AlignTransform::Translation(param0(1), param0(2), param0(3)),
3485                       AlignTransform::EulerAngles(param0(4), param0(5), param0(6)),
3486                       DetId(DetId::Calo).rawId());
3487 
3488   edm::LogVerbatim("GlobalTrackerMuonAlignment")
3489       << "Tracker (" << tracker.rawId() << ") at " << tracker.translation() << " " << tracker.rotation().eulerAngles();
3490   edm::LogVerbatim("GlobalTrackerMuonAlignment") << tracker.rotation();
3491 
3492   edm::LogVerbatim("GlobalTrackerMuonAlignment")
3493       << "Muon (" << muon.rawId() << ") at " << muon.translation() << " " << muon.rotation().eulerAngles();
3494   edm::LogVerbatim("GlobalTrackerMuonAlignment")
3495       << "          rotations angles around x,y,z "
3496       << " ( " << -muon.rotation().zy() << " " << muon.rotation().zx() << " " << -muon.rotation().yx() << " )";
3497   edm::LogVerbatim("GlobalTrackerMuonAlignment") << muon.rotation();
3498 
3499   edm::LogVerbatim("GlobalTrackerMuonAlignment")
3500       << "Ecal (" << ecal.rawId() << ") at " << ecal.translation() << " " << ecal.rotation().eulerAngles();
3501   edm::LogVerbatim("GlobalTrackerMuonAlignment") << ecal.rotation();
3502 
3503   edm::LogVerbatim("GlobalTrackerMuonAlignment")
3504       << "Hcal (" << hcal.rawId() << ") at " << hcal.translation() << " " << hcal.rotation().eulerAngles();
3505   edm::LogVerbatim("GlobalTrackerMuonAlignment") << hcal.rotation();
3506 
3507   edm::LogVerbatim("GlobalTrackerMuonAlignment")
3508       << "Calo (" << calo.rawId() << ") at " << calo.translation() << " " << calo.rotation().eulerAngles();
3509   edm::LogVerbatim("GlobalTrackerMuonAlignment") << calo.rotation();
3510 
3511   globalPositions.m_align.push_back(tracker);
3512   globalPositions.m_align.push_back(muon);
3513   globalPositions.m_align.push_back(ecal);
3514   globalPositions.m_align.push_back(hcal);
3515   globalPositions.m_align.push_back(calo);
3516 
3517   edm::LogVerbatim("GlobalTrackerMuonAlignment") << "Uploading to the database...";
3518 
3519   edm::Service<cond::service::PoolDBOutputService> poolDbService;
3520 
3521   if (!poolDbService.isAvailable())
3522     throw cms::Exception("NotAvailable") << "PoolDBOutputService not available";
3523 
3524   //    if (poolDbService->isNewTagRequest("GlobalPositionRcd")) {
3525   //       poolDbService->createOneIOV<Alignments>(globalPositions, poolDbService->endOfTime(), "GlobalPositionRcd");
3526   //    } else {
3527   //       poolDbService->appendOneIOV<Alignments>(globalPositions, poolDbService->currentTime(), "GlobalPositionRcd");
3528   //    }
3529   poolDbService->writeOneIOV<Alignments>(globalPositions, poolDbService->currentTime(), "GlobalPositionRcd");
3530   edm::LogVerbatim("GlobalTrackerMuonAlignment") << "done!";
3531 
3532   return;
3533 }
3534 
3535 //define this as a plug-in
3536 DEFINE_FWK_MODULE(GlobalTrackerMuonAlignment);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!