Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-08-13 05:00:11 

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