Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-12-17 23:35:47 

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