CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/Alignment/CommonAlignmentProducer/plugins/GlobalTrackerMuonAlignment.cc

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    GlobalTrackerMuonAlignment
00004 // Class:      GlobalTrackerMuonAlignment
00005 // 
00015 //
00016 // Original Author:  Alexandre Spiridonov
00017 //         Created:  Fri Oct 16 15:59:05 CEST 2009
00018 //
00019 // $Id: GlobalTrackerMuonAlignment.cc,v 1.9 2011/12/22 20:02:51 innocent Exp $
00020 //
00021 
00022 // system include files
00023 #include <memory>
00024 #include <string>
00025 #include <map>
00026 #include <vector>
00027 #include <fstream>
00028 #include <typeinfo>
00029 
00030 // user include files
00031 
00032 // Framework
00033 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00034 #include "FWCore/Framework/interface/EDAnalyzer.h"
00035 #include "FWCore/Framework/interface/ESWatcher.h"
00036 #include "FWCore/Framework/interface/Frameworkfwd.h"
00037 //#include "FWCore/Framework/interface/EDProducer.h"
00038 #include "FWCore/Framework/interface/Event.h"
00039 #include "FWCore/Framework/interface/EventSetup.h"
00040 #include "FWCore/Framework/interface/MakerMacros.h"
00041 #include "FWCore/Utilities/interface/InputTag.h"
00042 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00043 #include "FWCore/ServiceRegistry/interface/Service.h"
00044 //#include "FWCore/Framework/interface/EventSetupRecordImplementation.h"
00045 
00046 // references
00047 #include "DataFormats/DetId/interface/DetId.h"
00048 #include "DataFormats/MuonReco/interface/Muon.h"
00049 #include "DataFormats/MuonReco/interface/MuonFwd.h"
00050 #include "DataFormats/TrackReco/interface/Track.h"
00051 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00052 #include "DataFormats/GeometrySurface/interface/TangentPlane.h"
00053 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00054 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
00055 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
00056 #include "Geometry/TrackingGeometryAligner/interface/GeometryAligner.h"
00057 #include "MagneticField/Engine/interface/MagneticField.h"
00058 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00059 
00060 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBuilder.h"
00061 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00062 #include "TrackingTools/TrackRefitter/interface/TrackTransformer.h"
00063 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00064 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
00065 #include "TrackingTools/GeomPropagators/interface/SmartPropagator.h"
00066 #include "TrackingTools/GeomPropagators/interface/PropagationDirectionChooser.h" 
00067 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
00068 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
00069 #include "TrackPropagation/RungeKutta/interface/RKTestPropagator.h" 
00070 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00071 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
00072 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCartesianToLocal.h"
00073 
00074 // Database
00075 #include "FWCore/Framework/interface/ESHandle.h"
00076 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
00077 
00078 // Alignment
00079 #include "CondFormats/Alignment/interface/Alignments.h"
00080 #include "CondFormats/Alignment/interface/AlignTransform.h"
00081 #include "CondFormats/AlignmentRecord/interface/GlobalPositionRcd.h"
00082 
00083 // Refit
00084 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
00085 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
00086 #include "TrackingTools/TrackFitters/interface/KFTrajectoryFitter.h"
00087 #include "TrackingTools/TrackFitters/interface/KFTrajectorySmoother.h"
00088 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
00089 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
00090 
00091 #include <TH1.h>
00092 #include <TH2.h>
00093 #include <TFile.h>
00094 
00095 #include<iostream>
00096 #include<cmath>
00097 #include<assert.h>
00098 #include<fstream>
00099 #include<iomanip>
00100 
00101 using namespace edm; 
00102 using namespace std;
00103 using namespace reco;
00104 
00105 //
00106 // class declaration
00107 //
00108 
00109 class GlobalTrackerMuonAlignment : public edm::EDAnalyzer {
00110  public:
00111   explicit GlobalTrackerMuonAlignment(const edm::ParameterSet&);
00112   ~GlobalTrackerMuonAlignment();
00113 
00114   void analyzeTrackTrack(const edm::Event&, const edm::EventSetup&);
00115   void analyzeTrackTrajectory(const edm::Event&, const edm::EventSetup&);
00116   void bookHist();
00117   void fitHist();
00118   void trackFitter(reco::TrackRef, reco::TransientTrack&, PropagationDirection, 
00119                   TrajectoryStateOnSurface&);
00120   void muonFitter(reco::TrackRef, reco::TransientTrack&, PropagationDirection, 
00121                   TrajectoryStateOnSurface&);
00122   void debugTrackHit(const std::string, reco::TrackRef);
00123   void debugTrackHit(const std::string, reco::TransientTrack&);
00124   void debugTrajectorySOS(const std::string, TrajectoryStateOnSurface&);
00125   void debugTrajectorySOSv(const std::string, TrajectoryStateOnSurface);
00126   void debugTrajectory(const std::string, Trajectory&);
00127 
00128   void gradientGlobal(GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&, 
00129                       GlobalVector&, AlgebraicSymMatrix66&);
00130   void gradientLocal(GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&, 
00131                      GlobalVector&, CLHEP::HepSymMatrix&, CLHEP::HepMatrix&, 
00132                      CLHEP::HepVector&, AlgebraicVector4&);
00133   void gradientGlobalAlg(GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&,  
00134                          AlgebraicSymMatrix66&);
00135   void misalignMuon(GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&, 
00136                     GlobalVector&, GlobalVector&);
00137   void misalignMuonL(GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&, 
00138                      GlobalVector&, GlobalVector&, AlgebraicVector4&, 
00139                      TrajectoryStateOnSurface&, TrajectoryStateOnSurface&);
00140   void writeGlPosRcd(CLHEP::HepVector& d3); 
00141   inline double CLHEP_dot(CLHEP::HepVector a, CLHEP::HepVector b)
00142   {
00143     return a(1)*b(1)+a(2)*b(2)+a(3)*b(3);
00144   }
00145 
00146  private:
00147 
00148   virtual void beginJob() ;
00149   virtual void analyze(const edm::Event&, const edm::EventSetup&);
00150   virtual void endJob() ;
00151   
00152   // ----------member data ---------------------------
00153 
00154   edm::InputTag trackTags_;     // used to select what tracks to read from configuration file
00155   edm::InputTag muonTags_;      // used to select what standalone muons
00156   edm::InputTag gmuonTags_;     // used to select what global muons
00157   edm::InputTag smuonTags_;     // used to select what selected muons
00158   string propagator_;           // name of the propagator
00159   bool cosmicMuonMode_;   
00160   bool isolatedMuonMode_;
00161   bool collectionCosmic;        // if true, read muon collection expected for CosmicMu
00162   bool collectionIsolated;      //                                            IsolateMu
00163   bool refitMuon_;              // if true, refit stand alone muon track
00164   bool refitTrack_;             //                tracker track
00165   string rootOutFile_;
00166   string txtOutFile_;
00167   bool writeDB_;                // write results to DB
00168   bool debug_;                  // run in debugging mode
00169 
00170   edm::ESWatcher<GlobalTrackingGeometryRecord> watchTrackingGeometry_;
00171   edm::ESHandle<GlobalTrackingGeometry> theTrackingGeometry;
00172   const GlobalTrackingGeometry* trackingGeometry_;
00173 
00174   edm::ESWatcher<IdealMagneticFieldRecord> watchMagneticFieldRecord_;
00175   const MagneticField* magneticField_;
00176 
00177   edm::ESWatcher<TrackingComponentsRecord> watchTrackingComponentsRecord_;
00178 
00179   edm::ESWatcher<GlobalPositionRcd> watchGlobalPositionRcd_;
00180   const Alignments* globalPositionRcd_;
00181 
00182   KFTrajectoryFitter* theFitter;
00183   KFTrajectorySmoother* theSmoother;
00184   KFTrajectoryFitter* theFitterOp;
00185   KFTrajectorySmoother* theSmootherOp;
00186   bool defineFitter;
00187   MuonTransientTrackingRecHitBuilder* MuRHBuilder;
00188   const TransientTrackingRecHitBuilder* TTRHBuilder;
00189 
00190   //                            // LSF for global d(3):    Inf * d = Gfr 
00191   AlgebraicVector3 Gfr;         // free terms 
00192   AlgebraicSymMatrix33 Inf;     // information matrix 
00193   //                            // LSF for global d(6):    Hess * d = Grad 
00194   CLHEP::HepVector Grad;        // gradient of the objective function in global parameters
00195   CLHEP::HepMatrix Hess;        // Hessian                  -- // ---     
00196 
00197   CLHEP::HepVector GradL;       // gradient of the objective function in local parameters
00198   CLHEP::HepMatrix HessL;       // Hessian                  -- // ---     
00199 
00200   int N_event;                  // processed events
00201   int N_track;                  // selected tracks
00202 
00203   std::vector<AlignTransform>::const_iterator 
00204   iteratorTrackerRcd;           // location of Tracker in container globalPositionRcd_->m_aligm
00205   std::vector<AlignTransform>::const_iterator 
00206   iteratorMuonRcd;              //              Muon
00207   std::vector<AlignTransform>::const_iterator
00208   iteratorEcalRcd;              //              Ecal
00209   std::vector<AlignTransform>::const_iterator
00210   iteratorHcalRcd;              //              Hcal
00211 
00212   CLHEP::HepVector MuGlShift;   // evaluated global muon shifts  
00213   CLHEP::HepVector MuGlAngle;   // evaluated global muon angles  
00214 
00215   std::string MuSelect;         // what part of muon system is selected for 1st hit 
00216 
00217   ofstream OutGlobalTxt;        // output the vector of global alignment as text     
00218 
00219   TFile * file;
00220   TH1F * histo;
00221   TH1F * histo2; // outerP
00222   TH1F * histo3; // outerLambda = PI/2-outerTheta
00223   TH1F * histo4; // phi
00224   TH1F * histo5; // outerR
00225   TH1F * histo6; // outerZ
00226   TH1F * histo7; // s
00227   TH1F * histo8; // chi^2 of muon-track
00228 
00229   TH1F * histo11; // |Rm-Rt| 
00230   TH1F * histo12; // Xm-Xt 
00231   TH1F * histo13; // Ym-Yt 
00232   TH1F * histo14; // Zm-Zt 
00233   TH1F * histo15; // Nxm-Nxt 
00234   TH1F * histo16; // Nym-Nyt 
00235   TH1F * histo17; // Nzm-Nzt 
00236   TH1F * histo18; // Error X of inner track  
00237   TH1F * histo19; // Error X of muon 
00238   TH1F * histo20; // Error of Xm-Xt 
00239   TH1F * histo21; // pull(Xm-Xt) 
00240   TH1F * histo22; // pull(Ym-Yt) 
00241   TH1F * histo23; // pull(Zm-Zt) 
00242   TH1F * histo24; // pull(PXm-PXt) 
00243   TH1F * histo25; // pull(PYm-Pyt) 
00244   TH1F * histo26; // pull(PZm-PZt) 
00245   TH1F * histo27; // Nx of tangent plane  
00246   TH1F * histo28; // Ny of tangent plane  
00247   TH1F * histo29; // lenght of inner track  
00248   TH1F * histo30; // lenght of muon track  
00249   TH1F * histo31; // chi2 local  
00250   TH1F * histo32; // pull(Pxm/Pzm - Pxt/Pzt) local  
00251   TH1F * histo33; // pull(Pym/Pzm - Pyt/Pzt) local
00252   TH1F * histo34; // pull(Xm - Xt) local
00253   TH1F * histo35; // pull(Ym - Yt) local
00254 
00255   TH2F * histo101; // Rtrack/muon vs Ztrack/muon
00256   TH2F * histo102; // Ytrack/muon vs Xtrack/muon
00257 };
00258 
00259 //
00260 // constants, enums and typedefs
00261 //
00262 
00263 
00264 //
00265 // static data member definitions
00266 //
00267 
00268 //
00269 // constructors and destructor
00270 //
00271 GlobalTrackerMuonAlignment::GlobalTrackerMuonAlignment(const edm::ParameterSet& iConfig)
00272   :trackTags_(iConfig.getParameter<edm::InputTag>("tracks")),
00273    muonTags_(iConfig.getParameter<edm::InputTag>("muons")),
00274    gmuonTags_(iConfig.getParameter<edm::InputTag>("gmuons")),
00275    smuonTags_(iConfig.getParameter<edm::InputTag>("smuons")),
00276    propagator_(iConfig.getParameter<string>("Propagator")),
00277 
00278    cosmicMuonMode_(iConfig.getParameter<bool>("cosmics")),
00279    isolatedMuonMode_(iConfig.getParameter<bool>("isolated")),
00280 
00281    refitMuon_(iConfig.getParameter<bool>("refitmuon")),
00282    refitTrack_(iConfig.getParameter<bool>("refittrack")),
00283 
00284    rootOutFile_(iConfig.getUntrackedParameter<string>("rootOutFile")),
00285    txtOutFile_(iConfig.getUntrackedParameter<string>("txtOutFile")),
00286    writeDB_(iConfig.getUntrackedParameter<bool>("writeDB")),
00287    debug_(iConfig.getUntrackedParameter<bool>("debug")),
00288    defineFitter(true)
00289 {
00290 #ifdef NO_FALSE_FALSE_MODE 
00291   if (cosmicMuonMode_==false && isolatedMuonMode_==false) {
00292     throw cms::Exception("BadConfig") << "Muon collection not set! "
00293                                       << "Use from GlobalTrackerMuonAlignment_test_cfg.py.";
00294   }
00295 #endif  
00296   if (cosmicMuonMode_==true && isolatedMuonMode_==true) {
00297     throw cms::Exception("BadConfig") << "Muon collection can be either cosmic or isolated! "
00298                                       << "Please set only one to true.";
00299   }
00300 }
00301 
00302 GlobalTrackerMuonAlignment::~GlobalTrackerMuonAlignment()
00303 {
00304  
00305    // do anything here that needs to be done at desctruction time
00306    // (e.g. close files, deallocate resources etc.)
00307 
00308 }
00309 
00310 //
00311 // member functions
00312 //
00313 
00314 // ------------ method called to for each event  ------------
00315 void
00316 GlobalTrackerMuonAlignment::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00317 {  
00318   //GlobalTrackerMuonAlignment::analyzeTrackTrack(iEvent, iSetup);  
00319   GlobalTrackerMuonAlignment::analyzeTrackTrajectory(iEvent, iSetup);  
00320 }
00321 
00322 // ------------ method called once for job just before starting event loop  ------------
00323 void GlobalTrackerMuonAlignment::beginJob()
00324 {
00325   N_event = 0;
00326   N_track = 0;
00327   
00328   if (cosmicMuonMode_==true && isolatedMuonMode_==false) {
00329     collectionCosmic = true; collectionIsolated = false;}
00330   else if (cosmicMuonMode_==false && isolatedMuonMode_==true) {
00331     collectionCosmic = false; collectionIsolated = true;}
00332   else {
00333     collectionCosmic = false; collectionIsolated = false;}
00334 
00335   for(int i=0; i<=2; i++){
00336     Gfr(i) = 0.;
00337     for(int j=0; j<=2; j++){
00338       Inf(i,j) = 0.;
00339     }}
00340 
00341   Grad =  CLHEP::HepVector(6,0);
00342   Hess =  CLHEP::HepSymMatrix(6,0);
00343 
00344   GradL =  CLHEP::HepVector(6,0);
00345   HessL =  CLHEP::HepSymMatrix(6,0);
00346 
00347   // histograms
00348   TDirectory * dirsave = gDirectory;
00349 
00350   file = new TFile(rootOutFile_.c_str(), "recreate"); 
00351   const bool oldAddDir = TH1::AddDirectoryStatus(); 
00352 
00353   TH1::AddDirectory(true);
00354 
00355   this->bookHist();
00356 
00357   TH1::AddDirectory(oldAddDir);
00358   dirsave->cd();
00359 }
00360 
00361 // ------------ method called once each job just after ending the event loop  ------------
00362 void GlobalTrackerMuonAlignment::endJob() {
00363   bool alarm = false;
00364 
00365   this->fitHist();
00366 
00367   AlgebraicVector3 d(0., 0., 0.);   // ------------  alignmnet  Global Algebraic
00368 
00369   AlgebraicSymMatrix33 InfI;   // inverse it   
00370   for(int i=0; i<=2; i++) 
00371     for(int j=0; j<=2; j++){
00372       if(j < i) continue;
00373       InfI(i,j) += Inf(i,j);}
00374   bool ierr = !InfI.Invert();
00375   if( ierr ) { if(alarm)
00376       std::cout<<" Error inverse  Inf matrix !!!!!!!!!!!"<<std::endl;}
00377   
00378   for(int i=0; i<=2; i++) 
00379     for(int k=0; k<=2; k++) d(i) -= InfI(i,k)*Gfr(k);
00380   // end of Global Algebraic
00381         
00382   //                                ---------------  alignment Global CLHEP
00383   CLHEP::HepVector d3 = CLHEP::solve(Hess, -Grad);
00384   int iEr3;
00385   CLHEP::HepMatrix Errd3 = Hess.inverse(iEr3);
00386   if( iEr3 != 0) { if(alarm)
00387     std::cout<<" endJob Error inverse Hess matrix !!!!!!!!!!!"<<std::endl;
00388   }
00389   // end of Global CLHEP
00390   
00391   //                                ----------------- alignment Local CLHEP
00392   CLHEP::HepVector dLI = CLHEP::solve(HessL, -GradL);
00393   int iErI;
00394   CLHEP::HepMatrix ErrdLI = HessL.inverse(iErI);
00395   if( iErI != 0) { if(alarm)
00396       std::cout<<" endJob Error inverse HessL matrix !!!!!!!!!!!"<<std::endl;
00397   }
00398   // end of Local CLHEP
00399 
00400   // printout of final parameters
00401   std::cout<<" ---- "<<N_event<<" event "<<N_track<<" tracks "
00402            <<MuSelect<<" ---- "<<std::endl;
00403   if (collectionIsolated) std::cout<<" ALCARECOMuAlCalIsolatedMu "<<std::endl;
00404   else if (collectionCosmic) std::cout<<"  ALCARECOMuAlGlobalCosmics "<<std::endl; 
00405   else std::cout<<smuonTags_<<std::endl;
00406   std::cout<<" Similated shifts[cm] "
00407            <<MuGlShift(1)<<" "<<MuGlShift(2)<<" "<<MuGlShift(3)<<" "
00408            <<" angles[rad] "<<MuGlAngle(1)<<" "<<MuGlAngle(2)<<" "<<MuGlAngle(3)<<" "
00409            <<std::endl;
00410   std::cout<<" d    "<<-d<<std::endl;;
00411   std::cout<<" +-   "<<sqrt(InfI(0,0))<<" "<<sqrt(InfI(1,1))<<" "<<sqrt(InfI(2,2))
00412            <<std::endl;
00413   std::cout<<" dG   "<<d3(1)<<" "<<d3(2)<<" "<<d3(3)<<" "<<d3(4)<<" "
00414            <<d3(5)<<" "<<d3(6)<<std::endl;;
00415   std::cout<<" +-   "<<sqrt(Errd3(1,1))<<" "<<sqrt(Errd3(2,2))<<" "<<sqrt(Errd3(3,3))
00416            <<" "<<sqrt(Errd3(4,4))<<" "<<sqrt(Errd3(5,5))<<" "<<sqrt(Errd3(6,6))
00417            <<std::endl;
00418   std::cout<<" dL   "<<dLI(1)<<" "<<dLI(2)<<" "<<dLI(3)<<" "<<dLI(4)<<" "
00419            <<dLI(5)<<" "<<dLI(6)<<std::endl;;
00420   std::cout<<" +-   "<<sqrt(ErrdLI(1,1))<<" "<<sqrt(ErrdLI(2,2))<<" "<<sqrt(ErrdLI(3,3))
00421            <<" "<<sqrt(ErrdLI(4,4))<<" "<<sqrt(ErrdLI(5,5))<<" "<<sqrt(ErrdLI(6,6))
00422            <<std::endl;
00423 
00424   // what do we write to DB
00425   CLHEP::HepVector vectorToDb(6,0), vectorErrToDb(6,0);
00426   //vectorToDb = d3;
00427   //for(unsigned int i=1; i<=6; i++) vectorErrToDb(i) = sqrt(Errd3(i,i));
00428   vectorToDb = -dLI;
00429   for(unsigned int i=1; i<=6; i++) vectorErrToDb(i) = sqrt(ErrdLI(i,i));
00430 
00431   // write histograms to root file
00432   file->Write(); 
00433   file->Close(); 
00434 
00435   // write global parameters to text file 
00436   OutGlobalTxt.open(txtOutFile_.c_str(), ios::out);
00437   if(!OutGlobalTxt.is_open()) std::cout<<" outglobal.txt is not open !!!!!"<<std::endl;
00438   else{
00439     OutGlobalTxt<<vectorToDb(1)<<" "<<vectorToDb(2)<<" "<<vectorToDb(3)<<" "<<vectorToDb(4)<<" "
00440                 <<vectorToDb(5)<<" "<<vectorToDb(6)<<" muon Global.\n";
00441     OutGlobalTxt<<vectorErrToDb(1)<<" "<<vectorErrToDb(1)<<" "<<vectorErrToDb(1)<<" "
00442                 <<vectorErrToDb(1)<<" "<<vectorErrToDb(1)<<" "<<vectorErrToDb(1)<<" errors.\n";
00443     OutGlobalTxt<<N_event<<" events are processed.\n";
00444 
00445     if (collectionIsolated) OutGlobalTxt<<"ALCARECOMuAlCalIsolatedMu.\n";
00446     else if (collectionCosmic) OutGlobalTxt<<"  ALCARECOMuAlGlobalCosmics.\n"; 
00447     else OutGlobalTxt<<smuonTags_<<".\n";
00448     OutGlobalTxt.close();
00449     std::cout<<" Write to the file outglobal.txt done  "<<std::endl;
00450   }
00451 
00452   // write new GlobalPositionRcd to DB
00453   if(debug_)
00454     std::cout<<" writeBD_ "<<writeDB_<<std::endl;
00455   if (writeDB_) 
00456     GlobalTrackerMuonAlignment::writeGlPosRcd(vectorToDb);
00457 }
00458 
00459 // ------------ book histogram  ------------
00460 void GlobalTrackerMuonAlignment::bookHist()
00461 {
00462   double PI = 3.1415927;
00463   histo = new TH1F("Pt","pt",1000,0,100); 
00464   histo2 = new TH1F("P","P [GeV/c]",400,0.,400.);       
00465   histo2->GetXaxis()->SetTitle("momentum [GeV/c]");
00466   histo3 = new TH1F("outerLambda","#lambda outer",100,-PI/2.,PI/2.);    
00467   histo3->GetXaxis()->SetTitle("#lambda outer");
00468   histo4 = new TH1F("phi","#phi [rad]",100,-PI,PI);       
00469   histo4->GetXaxis()->SetTitle("#phi [rad]");
00470   histo5 = new TH1F("Rmuon","inner muon hit R [cm]",100,0.,800.);       
00471   histo5->GetXaxis()->SetTitle("R of muon [cm]");
00472   histo6 = new TH1F("Zmuon","inner muon hit Z[cm]",100,-1000.,1000.);       
00473   histo6->GetXaxis()->SetTitle("Z of muon [cm]");
00474   histo7 = new TH1F("(Pm-Pt)/Pt"," (Pmuon-Ptrack)/Ptrack", 100,-2.,2.);       
00475   histo7->GetXaxis()->SetTitle("(Pmuon-Ptrack)/Ptrack");
00476   histo8 = new TH1F("chi muon-track","#chi^{2}(muon-track)", 1000,0.,1000.);       
00477   histo8->GetXaxis()->SetTitle("#chi^{2} of muon w.r.t. propagated track");
00478   histo11 = new TH1F("distance muon-track","distance muon w.r.t track [cm]",100,0.,30.);
00479   histo11->GetXaxis()->SetTitle("distance of muon w.r.t. track [cm]");
00480   histo12 = new TH1F("Xmuon-Xtrack","Xmuon-Xtrack [cm]",200,-20.,20.);
00481   histo12->GetXaxis()->SetTitle("Xmuon - Xtrack [cm]"); 
00482   histo13 = new TH1F("Ymuon-Ytrack","Ymuon-Ytrack [cm]",200,-20.,20.);
00483   histo13->GetXaxis()->SetTitle("Ymuon - Ytrack [cm]");
00484   histo14 = new TH1F("Zmuon-Ztrack","Zmuon-Ztrack [cm]",200,-20.,20.);       
00485   histo14->GetXaxis()->SetTitle("Zmuon-Ztrack [cm]");
00486   histo15 = new TH1F("NXmuon-NXtrack","NXmuon-NXtrack [rad]",200,-.1,.1);       
00487   histo15->GetXaxis()->SetTitle("N_{X}(muon)-N_{X}(track) [rad]");
00488   histo16 = new TH1F("NYmuon-NYtrack","NYmuon-NYtrack [rad]",200,-.1,.1);       
00489   histo16->GetXaxis()->SetTitle("N_{Y}(muon)-N_{Y}(track) [rad]");
00490   histo17 = new TH1F("NZmuon-NZtrack","NZmuon-NZtrack [rad]",200,-.1,.1);       
00491   histo17->GetXaxis()->SetTitle("N_{Z}(muon)-N_{Z}(track) [rad]");
00492   histo18 = new TH1F("expected error of Xinner","outer hit of inner tracker",100,0,.01);
00493   histo18->GetXaxis()->SetTitle("expected error of Xinner [cm]");
00494   histo19 = new TH1F("expected error of Xmuon","inner hit of muon",100,0,.1);       
00495   histo19->GetXaxis()->SetTitle("expected error of Xmuon [cm]");
00496   histo20 = new TH1F("expected error of Xmuon-Xtrack","muon w.r.t. propagated track",100,0.,10.);
00497   histo20->GetXaxis()->SetTitle("expected error of Xmuon-Xtrack [cm]");
00498   histo21 = new TH1F("pull of Xmuon-Xtrack","pull of Xmuon-Xtrack",100,-10.,10.);
00499   histo21->GetXaxis()->SetTitle("(Xmuon-Xtrack)/expected error");
00500   histo22 = new TH1F("pull of Ymuon-Ytrack","pull of Ymuon-Ytrack",100,-10.,10.);       
00501   histo22->GetXaxis()->SetTitle("(Ymuon-Ytrack)/expected error");
00502   histo23 = new TH1F("pull of Zmuon-Ztrack","pull of Zmuon-Ztrack",100,-10.,10.);       
00503   histo23->GetXaxis()->SetTitle("(Zmuon-Ztrack)/expected error");
00504   histo24 = new TH1F("pull of PXmuon-PXtrack","pull of PXmuon-PXtrack",100,-10.,10.);       
00505   histo24->GetXaxis()->SetTitle("(P_{X}(muon)-P_{X}(track))/expected error");
00506   histo25 = new TH1F("pull of PYmuon-PYtrack","pull of PYmuon-PYtrack",100,-10.,10.);       
00507   histo25->GetXaxis()->SetTitle("(P_{Y}(muon)-P_{Y}(track))/expected error");
00508   histo26 = new TH1F("pull of PZmuon-PZtrack","pull of PZmuon-PZtrack",100,-10.,10.);       
00509   histo26->GetXaxis()->SetTitle("(P_{Z}(muon)-P_{Z}(track))/expected error");
00510   histo27 = new TH1F("N_x","Nx of tangent plane",120,-1.1,1.1); 
00511   histo27->GetXaxis()->SetTitle("normal vector projection N_{X}");
00512   histo28 = new TH1F("N_y","Ny of tangent plane",120,-1.1,1.1); 
00513   histo28->GetXaxis()->SetTitle("normal vector projection N_{Y}");
00514   histo29 = new TH1F("lenght of track","lenght of track",200,0.,400); 
00515   histo29->GetXaxis()->SetTitle("lenght of track [cm]");
00516   histo30 = new TH1F("lenght of muon","lenght of muon",200,0.,800); 
00517   histo30->GetXaxis()->SetTitle("lenght of muon [cm]");
00518 
00519   histo31 = new TH1F("local chi muon-track","#local chi^{2}(muon-track)", 1000,0.,1000.);       
00520   histo31->GetXaxis()->SetTitle("#local chi^{2} of muon w.r.t. propagated track");
00521   histo32 = new TH1F("pull of Px/Pz local","pull of Px/Pz local",100,-10.,10.);
00522   histo32->GetXaxis()->SetTitle("local (Px/Pz(muon) - Px/Pz(track))/expected error");
00523   histo33 = new TH1F("pull of Py/Pz local","pull of Py/Pz local",100,-10.,10.);
00524   histo33->GetXaxis()->SetTitle("local (Py/Pz(muon) - Py/Pz(track))/expected error");
00525   histo34 = new TH1F("pull of X local","pull of X local",100,-10.,10.);
00526   histo34->GetXaxis()->SetTitle("local (Xmuon - Xtrack)/expected error");
00527   histo35 = new TH1F("pull of Y local","pull of Y local",100,-10.,10.);
00528   histo35->GetXaxis()->SetTitle("local (Ymuon - Ytrack)/expected error");
00529 
00530   histo101 = new TH2F("Rtr/mu vs Ztr/mu","hit of track/muon",100,-800.,800.,100,0.,600.);
00531   histo101->GetXaxis()->SetTitle("Z of track/muon [cm]");
00532   histo101->GetYaxis()->SetTitle("R of track/muon [cm]");
00533   histo102 = new TH2F("Ytr/mu vs Xtr/mu","hit of track/muon",100,-600.,600.,100,-600.,600.);
00534   histo102->GetXaxis()->SetTitle("X of track/muon [cm]");
00535   histo102->GetYaxis()->SetTitle("Y of track/muon [cm]");
00536 }
00537 
00538 // ------------ fit histogram  ------------
00539 void GlobalTrackerMuonAlignment::fitHist() {
00540 
00541   histo7->Fit("gaus","Q");
00542 
00543   histo12->Fit("gaus","Q");
00544   histo13->Fit("gaus","Q");
00545   histo14->Fit("gaus","Q");
00546   histo15->Fit("gaus","Q");
00547   histo16->Fit("gaus","Q");
00548   histo17->Fit("gaus","Q");
00549 
00550   histo21->Fit("gaus","Q");
00551   histo22->Fit("gaus","Q");
00552   histo23->Fit("gaus","Q");
00553   histo24->Fit("gaus","Q");
00554   histo25->Fit("gaus","Q");
00555   histo26->Fit("gaus","Q");
00556 
00557   histo32->Fit("gaus","Q");
00558   histo33->Fit("gaus","Q");
00559   histo34->Fit("gaus","Q");
00560   histo35->Fit("gaus","Q");
00561 }
00562 
00563 // ------------ method to analyze recoTrack & recoMuon of Global Muon  ------------
00564 void GlobalTrackerMuonAlignment::analyzeTrackTrack
00565 (const edm::Event& iEvent, const edm::EventSetup& iSetup)
00566 {
00567   using namespace edm;
00568   using namespace reco;  
00569   //debug_ = true;
00570   bool info = false;
00571   bool alarm = false;
00572   //bool alarm = true;
00573   double PI = 3.1415927;
00574  
00575   cosmicMuonMode_ = true; // true: both Cosmic and IsolatedMuon are treated with 1,2 tracks
00576   isolatedMuonMode_ = !cosmicMuonMode_; //true: only IsolatedMuon are treated with 1 track
00577  
00578   //if(N_event == 1)
00579   //GlobalPositionRcdRead1::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup);
00580 
00581   N_event++;
00582   if (info || debug_) {
00583     std::cout << "----- event " << N_event << " -- tracks "<<N_track<<" ---";
00584     if (cosmicMuonMode_) std::cout << " treated as CosmicMu ";
00585     if (isolatedMuonMode_) std::cout <<" treated as IsolatedMu ";
00586     std::cout << std::endl;
00587   }
00588 
00589   Handle<reco::TrackCollection> tracks;
00590   Handle<reco::TrackCollection> muons;
00591   Handle<reco::TrackCollection> gmuons;
00592   Handle<reco::MuonCollection> smuons;
00593 
00594   if (collectionIsolated){
00595     iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "TrackerOnly",  tracks);
00596     iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "StandAlone",  muons);
00597     iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "GlobalMuon",  gmuons);
00598     iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "SelectedMuons",  smuons);
00599   }
00600   else if (collectionCosmic){
00601     iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "TrackerOnly",  tracks);
00602     iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "StandAlone",  muons);
00603     iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "GlobalMuon",  gmuons);
00604     iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "SelectedMuons",  smuons);
00605   }
00606   else{
00607     iEvent.getByLabel(trackTags_,tracks);
00608     iEvent.getByLabel(muonTags_,muons);
00609     iEvent.getByLabel(gmuonTags_,gmuons);
00610     iEvent.getByLabel(smuonTags_,smuons);
00611   }  
00612 
00613   if (debug_) {
00614     std::cout << " ievBunch " << iEvent.bunchCrossing()
00615               << " runN " << (int)iEvent.run() <<std::endl;
00616     std::cout << " N tracks s/amu gmu selmu "<<tracks->size() << " " <<muons->size()
00617               << " "<<gmuons->size() << " " << smuons->size() << std::endl;
00618     for (MuonCollection::const_iterator itMuon = smuons->begin();
00619         itMuon != smuons->end();
00620          ++itMuon) {
00621       std::cout << " is isolatValid Matches " << itMuon->isIsolationValid()
00622                 << " " <<itMuon->isMatchesValid() << std::endl;
00623     }
00624   }
00625 
00626   if (isolatedMuonMode_) {                // ---- Only 1 Isolated Muon --------
00627     if (tracks->size() != 1) return;  
00628     if (muons->size()  != 1) return;  
00629     if (gmuons->size() != 1) return; 
00630     if (smuons->size() != 1) return; 
00631   }
00632   
00633   if (cosmicMuonMode_){                // ---- 1,2 Cosmic Muon --------
00634     if (smuons->size() > 2) return; 
00635     if (tracks->size() != smuons->size()) return;  
00636     if (muons->size() != smuons->size()) return;  
00637     if (gmuons->size() != smuons->size()) return;  
00638   }
00639   
00640   // ok mc_isolated_mu
00641   //edm::ESHandle<TrackerGeometry> trackerGeometry;
00642   //iSetup.get<TrackerDigiGeometryRecord>().get(trackerGeometry);
00643   // ok mc_isolated_mu
00644   //edm::ESHandle<DTGeometry> dtGeometry;
00645   //iSetup.get<MuonGeometryRecord>().get(dtGeometry);
00646   // don't work
00647   //edm::ESHandle<CSCGeometry> cscGeometry;
00648   //iSetup.get<MuonGeometryRecord>().get(cscGeometry);
00649   
00650   if (watchTrackingGeometry_.check(iSetup) || !trackingGeometry_) {
00651     edm::ESHandle<GlobalTrackingGeometry> trackingGeometry;
00652     iSetup.get<GlobalTrackingGeometryRecord>().get(trackingGeometry);
00653     trackingGeometry_ = &*trackingGeometry;
00654   }
00655   
00656   if (watchMagneticFieldRecord_.check(iSetup) || !magneticField_) {
00657     edm::ESHandle<MagneticField> magneticField;
00658     iSetup.get<IdealMagneticFieldRecord>().get(magneticField);
00659     magneticField_ = &*magneticField;
00660   }
00661 
00662   if (watchGlobalPositionRcd_.check(iSetup) || !globalPositionRcd_) {
00663     edm::ESHandle<Alignments> globalPositionRcd;
00664     iSetup.get<GlobalPositionRcd>().get(globalPositionRcd);
00665     globalPositionRcd_ = &*globalPositionRcd;
00666     for (std::vector<AlignTransform>::const_iterator i 
00667            = globalPositionRcd_->m_align.begin();
00668          i != globalPositionRcd_->m_align.end();  ++i){
00669       if(DetId(DetId::Tracker).rawId() == i->rawId()) iteratorTrackerRcd = i;
00670       if(DetId(DetId::Muon).rawId() == i->rawId()) iteratorMuonRcd = i;
00671       if(DetId(DetId::Ecal).rawId() == i->rawId()) iteratorEcalRcd = i;
00672       if(DetId(DetId::Hcal).rawId() == i->rawId()) iteratorHcalRcd = i;
00673     }
00674     if(true || debug_){
00675       std::cout << "=== iteratorMuonRcd " << iteratorMuonRcd->rawId()<<" ====\n" 
00676                 << " translation " << iteratorMuonRcd->translation()<<"\n" 
00677                 << " angles " << iteratorMuonRcd->rotation().eulerAngles() << std::endl;
00678       std::cout << iteratorMuonRcd->rotation() << std::endl;
00679     }   
00680   } // end of GlobalPositionRcd
00681   
00682   ESHandle<Propagator> propagator;
00683   iSetup.get<TrackingComponentsRecord>().get(propagator_, propagator);
00684   
00685   SteppingHelixPropagator alongStHePr = 
00686     SteppingHelixPropagator(magneticField_, alongMomentum);  
00687   SteppingHelixPropagator oppositeStHePr = 
00688     SteppingHelixPropagator(magneticField_, oppositeToMomentum);  
00689 
00690   //double tolerance = 5.e-5;
00691   RKTestPropagator alongRKPr = 
00692     RKTestPropagator(magneticField_, alongMomentum);
00693   RKTestPropagator oppositeRKPr = 
00694     RKTestPropagator(magneticField_, oppositeToMomentum);
00695 
00696   float epsilon = 5.;
00697   SmartPropagator alongSmPr = 
00698     SmartPropagator(alongRKPr, alongStHePr, magneticField_, alongMomentum, epsilon);
00699   SmartPropagator oppositeSmPr = 
00700     SmartPropagator(oppositeRKPr, oppositeStHePr, magneticField_, oppositeToMomentum, epsilon);
00701   
00702   // ................................................ selected/global muon
00703   //itMuon -->  Jim's globalMuon  
00704   for(MuonCollection::const_iterator itMuon = smuons->begin();
00705       itMuon != smuons->end(); ++itMuon) {
00706     
00707     if(debug_){
00708       std::cout<<" mu gM is GM Mu SaM tM "<<itMuon->isGlobalMuon()<<" "
00709                <<itMuon->isMuon()<<" "<<itMuon->isStandAloneMuon()
00710                <<" "<<itMuon->isTrackerMuon()<<" "
00711                <<std::endl;
00712     }
00713     
00714     // get information about the innermost muon hit -------------------------
00715     TransientTrack muTT(itMuon->outerTrack(), magneticField_, trackingGeometry_);
00716     TrajectoryStateOnSurface innerMuTSOS = muTT.innermostMeasurementState();
00717     TrajectoryStateOnSurface outerMuTSOS = muTT.outermostMeasurementState();
00718     
00719     // get information about the outermost tracker hit -----------------------
00720     TransientTrack trackTT(itMuon->track(), magneticField_, trackingGeometry_);
00721     TrajectoryStateOnSurface outerTrackTSOS = trackTT.outermostMeasurementState();
00722     TrajectoryStateOnSurface innerTrackTSOS = trackTT.innermostMeasurementState();
00723     
00724     GlobalPoint pointTrackIn  = innerTrackTSOS.globalPosition();
00725     GlobalPoint pointTrackOut = outerTrackTSOS.globalPosition();
00726     float lenghtTrack = (pointTrackOut-pointTrackIn).mag();
00727     GlobalPoint pointMuonIn  = innerMuTSOS.globalPosition();
00728     GlobalPoint pointMuonOut = outerMuTSOS.globalPosition();
00729     float lenghtMuon = (pointMuonOut - pointMuonIn).mag();
00730     GlobalVector momentumTrackOut = outerTrackTSOS.globalMomentum();
00731     GlobalVector momentumTrackIn = innerTrackTSOS.globalMomentum();
00732     float distanceInIn   = (pointTrackIn  - pointMuonIn).mag();
00733     float distanceInOut  = (pointTrackIn  - pointMuonOut).mag();
00734     float distanceOutIn  = (pointTrackOut - pointMuonIn).mag();
00735     float distanceOutOut = (pointTrackOut - pointMuonOut).mag();
00736 
00737     if(debug_){
00738       std::cout<<" pointTrackIn "<<pointTrackIn<<std::endl;
00739       std::cout<<"          Out "<<pointTrackOut<<" lenght "<<lenghtTrack<<std::endl;
00740       std::cout<<" point MuonIn "<<pointMuonIn<<std::endl;
00741       std::cout<<"          Out "<<pointMuonOut<<" lenght "<<lenghtMuon<<std::endl;
00742       std::cout<<" momeTrackIn Out "<<momentumTrackIn<<" "<<momentumTrackOut<<std::endl;
00743       std::cout<<" doi io ii oo "<<distanceOutIn<<" "<<distanceInOut<<" "
00744                <<distanceInIn<<" "<<distanceOutOut<<std::endl; 
00745   }   
00746 
00747     if(lenghtTrack < 90.) continue;
00748     if(lenghtMuon < 300.) continue;
00749     if(momentumTrackIn.mag() < 15. || momentumTrackOut.mag() < 15.) continue;
00750     if(trackTT.charge() != muTT.charge()) continue;
00751 
00752     if(debug_)
00753       std::cout<<" passed lenght momentum cuts"<<std::endl;
00754 
00755     GlobalVector GRm, GPm, Nl, Rm, Pm, Rt, Pt, Rt0;
00756     AlgebraicSymMatrix66 Cm, C0, Ce, C1;
00757 
00758     TrajectoryStateOnSurface extrapolationT;
00759     int tsosMuonIf = 0;
00760 
00761     if( isolatedMuonMode_ ){      //------------------------------- Isolated Muon -----
00762       const Surface& refSurface = innerMuTSOS.surface(); 
00763       ReferenceCountingPointer<TangentPlane> 
00764         tpMuLocal(refSurface.tangentPlane(innerMuTSOS.localPosition()));
00765       Nl = tpMuLocal->normalVector();
00766       ReferenceCountingPointer<TangentPlane> 
00767       tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
00768       GlobalVector Ng = tpMuGlobal->normalVector();
00769       Surface* surf = (Surface*)&refSurface;
00770       const Plane* refPlane = dynamic_cast<Plane*>(surf); 
00771       GlobalVector Nlp = refPlane->normalVector();
00772       
00773       if(debug_){
00774         std::cout<<" Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
00775          <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
00776         std::cout<<"  global "<<Ng.x()<<" "<<Ng.y()<<" "<<Ng.z()<<std::endl;
00777         std::cout<<"      lp "<<Nlp.x()<<" "<<Nlp.y()<<" "<<Nlp.z()<<std::endl;
00778         //std::cout<<" Nlocal Nglobal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
00779         // <<Ng.x()<<" "<<Ng.y()<<" "<<Ng.z()
00780         //<<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
00781       }
00782 
00783       //                          extrapolation to innermost muon hit     
00784       extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
00785       //extrapolationT = propagator->propagate(outerTrackTSOS, refSurface);
00786       
00787       if(!extrapolationT.isValid()){
00788         if(false & alarm)
00789           std::cout<<" !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
00790             //<<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
00791                    <<std::endl;
00792         continue;
00793       }
00794       tsosMuonIf = 1;
00795       Rt = GlobalVector((extrapolationT.globalPosition()).x(),
00796                         (extrapolationT.globalPosition()).y(),
00797                         (extrapolationT.globalPosition()).z());
00798       
00799       Pt = extrapolationT.globalMomentum();
00800       //                          global parameters of muon
00801       GRm = GlobalVector((innerMuTSOS.globalPosition()).x(),
00802                          (innerMuTSOS.globalPosition()).y(),
00803                          (innerMuTSOS.globalPosition()).z());
00804       GPm = innerMuTSOS.globalMomentum();
00805       Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
00806                          (outerTrackTSOS.globalPosition()).y(),
00807                          (outerTrackTSOS.globalPosition()).z());
00808       Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() + 
00809                                  innerMuTSOS.cartesianError().matrix());
00810       C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());                 
00811       Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix()); 
00812       C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
00813      
00814     }  //                    -------------------------------   end Isolated Muon -----
00815     
00816     
00817     if( cosmicMuonMode_ ){      //------------------------------- Cosmic Muon -----
00818 
00819       if((distanceOutIn <= distanceInOut) & (distanceOutIn <= distanceInIn) & 
00820          (distanceOutIn <= distanceOutOut)){             // ----- Out - In ------    
00821         if(debug_) std::cout<<" -----    Out - In -----"<<std::endl;
00822 
00823         const Surface& refSurface = innerMuTSOS.surface(); 
00824         ReferenceCountingPointer<TangentPlane> 
00825           tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
00826         Nl = tpMuGlobal->normalVector();
00827 
00828         //                          extrapolation to innermost muon hit     
00829         extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
00830         //extrapolationT = propagator->propagate(outerTrackTSOS, refSurface);
00831 
00832         if(!extrapolationT.isValid()){
00833           if(false & alarm)
00834             std::cout<<" !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
00835               //<<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
00836                      <<std::endl;
00837           continue;
00838         }
00839         if(debug_)
00840           std::cout<<" extrapolationT.isValid "<<extrapolationT.isValid()<<std::endl; 
00841 
00842         tsosMuonIf = 1;
00843         Rt = GlobalVector((extrapolationT.globalPosition()).x(),
00844                           (extrapolationT.globalPosition()).y(),
00845                           (extrapolationT.globalPosition()).z());
00846         
00847         Pt = extrapolationT.globalMomentum();
00848         //                          global parameters of muon
00849         GRm = GlobalVector((innerMuTSOS.globalPosition()).x(),
00850                            (innerMuTSOS.globalPosition()).y(),
00851                            (innerMuTSOS.globalPosition()).z());
00852         GPm = innerMuTSOS.globalMomentum();
00853         Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
00854                            (outerTrackTSOS.globalPosition()).y(),
00855                            (outerTrackTSOS.globalPosition()).z());
00856         Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() + 
00857                                    innerMuTSOS.cartesianError().matrix());
00858         C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());               
00859         Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix()); 
00860         C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
00861         
00862         if(false & debug_){
00863           //std::cout<<" ->propDir "<<propagator->propagationDirection()<<std::endl;
00864           PropagationDirectionChooser Chooser;
00865           std::cout<<" propDirCh "
00866                    <<Chooser.operator()(*outerTrackTSOS.freeState(), refSurface)
00867                    <<" Ch == along "<<(alongMomentum == 
00868                                        Chooser.operator()(*outerTrackTSOS.freeState(), refSurface))
00869                    <<std::endl;
00870           std::cout<<" --- Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
00871                    <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
00872         }
00873       }           //                                       enf of ---- Out - In -----
00874       
00875 
00876       else if((distanceInOut <= distanceInIn) & (distanceInOut <= distanceOutIn) & 
00877               (distanceInOut <= distanceOutOut)){            // ----- In - Out ------    
00878         if(debug_) std::cout<<" -----    In - Out -----"<<std::endl;
00879 
00880         const Surface& refSurface = outerMuTSOS.surface(); 
00881         ReferenceCountingPointer<TangentPlane> 
00882           tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
00883         Nl = tpMuGlobal->normalVector();
00884                 
00885         //                          extrapolation to outermost muon hit     
00886         extrapolationT = oppositeSmPr.propagate(innerTrackTSOS, refSurface);
00887 
00888         if(!extrapolationT.isValid()){
00889           if(false & alarm)
00890             std::cout<<" !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
00891               <<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
00892                      <<std::endl;
00893           continue;
00894         }
00895         if(debug_)
00896           std::cout<<" extrapolationT.isValid "<<extrapolationT.isValid()<<std::endl; 
00897         
00898         tsosMuonIf = 2;
00899         Rt = GlobalVector((extrapolationT.globalPosition()).x(),
00900                           (extrapolationT.globalPosition()).y(),
00901                           (extrapolationT.globalPosition()).z());
00902         
00903         Pt = extrapolationT.globalMomentum();
00904         //                          global parameters of muon
00905         GRm = GlobalVector((outerMuTSOS.globalPosition()).x(),
00906                            (outerMuTSOS.globalPosition()).y(),
00907                            (outerMuTSOS.globalPosition()).z());
00908         GPm = outerMuTSOS.globalMomentum();
00909         Rt0 = GlobalVector((innerTrackTSOS.globalPosition()).x(),
00910                            (innerTrackTSOS.globalPosition()).y(),
00911                            (innerTrackTSOS.globalPosition()).z());
00912         Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() + 
00913                                    outerMuTSOS.cartesianError().matrix());
00914         C0 = AlgebraicSymMatrix66(innerTrackTSOS.cartesianError().matrix());               
00915         Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix()); 
00916         C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
00917         
00918         if(false & debug_){
00919           //std::cout<<" ->propDir "<<propagator->propagationDirection()<<std::endl;
00920           PropagationDirectionChooser Chooser;
00921           std::cout<<" propDirCh "
00922                    <<Chooser.operator()(*innerTrackTSOS.freeState(), refSurface)
00923                    <<" Ch == oppisite "<<(oppositeToMomentum == 
00924                                           Chooser.operator()(*innerTrackTSOS.freeState(), refSurface))
00925                    <<std::endl;
00926           std::cout<<" --- Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
00927                    <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
00928         }
00929       }           //                                        enf of ---- In - Out -----
00930 
00931       else if((distanceOutOut <= distanceInOut) & (distanceOutOut <= distanceInIn) & 
00932               (distanceOutOut <= distanceOutIn)){          // ----- Out - Out ------    
00933         if(debug_) std::cout<<" -----    Out - Out -----"<<std::endl;
00934 
00935         // reject: momentum of track has opposite direction to muon track
00936         continue;
00937 
00938         const Surface& refSurface = outerMuTSOS.surface(); 
00939         ReferenceCountingPointer<TangentPlane> 
00940           tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
00941         Nl = tpMuGlobal->normalVector();
00942                 
00943         //                          extrapolation to outermost muon hit     
00944         extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
00945 
00946         if(!extrapolationT.isValid()){
00947           if(alarm)
00948             std::cout<<" !!!!!!!!! Catastrophe Out-Out extrapolationT.isValid "
00949                      <<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
00950                      <<std::endl;
00951           continue;
00952         }
00953         if(debug_)
00954           std::cout<<" extrapolationT.isValid "<<extrapolationT.isValid()<<std::endl; 
00955 
00956         tsosMuonIf = 2;
00957         Rt = GlobalVector((extrapolationT.globalPosition()).x(),
00958                           (extrapolationT.globalPosition()).y(),
00959                           (extrapolationT.globalPosition()).z());
00960         
00961         Pt = extrapolationT.globalMomentum();
00962         //                          global parameters of muon
00963         GRm = GlobalVector((outerMuTSOS.globalPosition()).x(),
00964                            (outerMuTSOS.globalPosition()).y(),
00965                            (outerMuTSOS.globalPosition()).z());
00966         GPm = outerMuTSOS.globalMomentum();
00967         Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
00968                            (outerTrackTSOS.globalPosition()).y(),
00969                            (outerTrackTSOS.globalPosition()).z());
00970         Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() + 
00971                                    outerMuTSOS.cartesianError().matrix());
00972         C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());               
00973         Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix()); 
00974         C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
00975         
00976         if(debug_){
00977           PropagationDirectionChooser Chooser;
00978           std::cout<<" propDirCh "
00979                    <<Chooser.operator()(*outerTrackTSOS.freeState(), refSurface)
00980                    <<" Ch == along "<<(alongMomentum == 
00981                                        Chooser.operator()(*outerTrackTSOS.freeState(), refSurface))
00982                    <<std::endl;
00983           std::cout<<" --- Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
00984                    <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
00985           std::cout<<"     Nornal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
00986                    <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
00987         }
00988       }           //                                       enf of ---- Out - Out -----
00989       else{
00990         if(alarm)
00991           std::cout<<" ------- !!!!!!!!!!!!!!! No proper Out - In \a\a\a\a\a\a\a"<<std::endl; 
00992         continue;
00993       }
00994       
00995     }  //                     -------------------------------   end Cosmic Muon -----
00996     
00997     if(tsosMuonIf == 0) {if(info) {std::cout<<"No tsosMuon !!!!!!"<<std::endl; continue;}}
00998     TrajectoryStateOnSurface tsosMuon;
00999     if(tsosMuonIf == 1) tsosMuon = muTT.innermostMeasurementState();
01000     else tsosMuon = muTT.outermostMeasurementState();
01001 
01002     //GlobalTrackerMuonAlignment::misalignMuon(GRm, GPm, Nl, Rt, Rm, Pm);
01003     AlgebraicVector4 LPRm;  // muon local (dx/dz, dy/dz, x, y)
01004     GlobalTrackerMuonAlignment::misalignMuonL
01005     (GRm, GPm, Nl, Rt, Rm, Pm, LPRm, extrapolationT, tsosMuon);
01006 
01007     GlobalVector resR = Rm - Rt;
01008     GlobalVector resP0 = Pm - Pt;
01009     GlobalVector resP = Pm / Pm.mag() - Pt / Pt.mag();
01010     float RelMomResidual = (Pm.mag() - Pt.mag()) / (Pt.mag() + 1.e-6);;
01011 
01012     AlgebraicVector6 Vm;
01013     Vm(0) = resR.x();  Vm(1) = resR.y();  Vm(2) = resR.z(); 
01014     Vm(3) = resP0.x();  Vm(4) = resP0.y();  Vm(5) = resP0.z(); 
01015     float Rmuon = Rm.perp();
01016     float Zmuon = Rm.z();
01017     float alfa_x = atan2(Nl.x(),Nl.y())*57.29578;
01018     
01019     if(debug_){
01020       std::cout<<" Nx Ny Nz alfa_x "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "<<alfa_x<<std::endl;
01021       //std::cout<<" Rm "<<Rm<<std::endl<<" Rt "<<Rt<<std::endl<<" resR "<<resR<<std::endl
01022       //       <<" resP "<<resP<<" dp/p "<<RelMomResidual<<std::endl;
01023     }
01024 
01025     double chi_d = 0;
01026     for(int i=0; i<=5; i++) chi_d += Vm(i)*Vm(i)/Cm(i,i);
01027         
01028     AlgebraicVector5 Vml(tsosMuon.localParameters().vector() 
01029                        - extrapolationT.localParameters().vector());
01030     AlgebraicSymMatrix55 m(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
01031     AlgebraicSymMatrix55 Cml(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
01032     bool ierrLoc = !m.Invert();
01033     if (ierrLoc && debug_ && info) { 
01034       std::cout<< " ==== Error inverting Local covariance matrix ==== "<<std::cout;
01035       continue;}
01036     double chi_Loc = ROOT::Math::Similarity(Vml,m);
01037     if(debug_)
01038       std::cout<<" chi_Loc px/pz/err "<<chi_Loc<<" "<<Vml(1)/sqrt(Cml(1,1))<<std::endl;
01039 
01040     if(Pt.mag() < 15.) continue;
01041     if(Pm.mag() <  5.) continue;
01042  
01043     //if(Pt.mag() < 30.) continue;          // momenum cut  < 30GeV
01044     //if(Pt.mag() < 60.) continue;          // momenum cut  < 30GeV
01045     //if(Pt.mag() > 50.) continue;          //  momenum cut > 50GeV
01046     //if(Pt.mag() > 100.) continue;          //  momenum cut > 100GeV
01047     //if(trackTT.charge() < 0) continue;    // select positive charge
01048     //if(trackTT.charge() > 0) continue;    // select negative charge
01049     
01050     //if(fabs(resR.x()) > 5.) continue;     // strong cut  X 
01051     //if(fabs(resR.y()) > 5.) continue;     //             Y 
01052     //if(fabs(resR.z()) > 5.) continue;     //             Z
01053     //if(fabs(resR.mag()) > 7.5) continue;   //            dR
01054 
01055     /*
01056     //if(fabs(RelMomResidual) > 0.5) continue;
01057     if(fabs(resR.x()) > 20.) continue;
01058     if(fabs(resR.y()) > 20.) continue;
01059     if(fabs(resR.z()) > 20.) continue;
01060     if(fabs(resR.mag()) > 30.) continue;
01061     if(fabs(resP.x()) > 0.06) continue;
01062     if(fabs(resP.y()) > 0.06) continue;
01063     if(fabs(resP.z()) > 0.06) continue;
01064     if(chi_d > 40.) continue;
01065     */
01066 
01067     //                                            select Barrel 
01068     //if(Rmuon < 400. || Rmuon > 450.) continue; 
01069     //if(Zmuon < -600. || Zmuon > 600.) continue;
01070     //if(fabs(Nl.z()) > 0.95) continue;
01071     //MuSelect = " Barrel";
01072     //                                                  EndCap1
01073     //if(Rmuon < 120. || Rmuon > 450.) continue;
01074     //if(Zmuon < -720.) continue;
01075     //if(Zmuon > -580.) continue;
01076     //if(fabs(Nl.z()) < 0.95) continue;  
01077     //MuSelect = " EndCap1";
01078     //                                                  EndCap2
01079     //if(Rmuon < 120. || Rmuon > 450.) continue;
01080     //if(Zmuon >  720.) continue;
01081     //if(Zmuon <  580.) continue;
01082     //if(fabs(Nl.z()) < 0.95) continue;  
01083     //MuSelect = " EndCap2";
01084     //                                                 select All
01085     if(Rmuon < 120. || Rmuon > 450.) continue;  
01086     if(Zmuon < -720. || Zmuon > 720.) continue;
01087     MuSelect = " Barrel+EndCaps";
01088 
01089     
01090     if(debug_)
01091       std::cout<<" .............. passed all cuts"<<std::endl;
01092         
01093     N_track++;
01094     //                     gradient and Hessian for each track
01095     
01096     GlobalTrackerMuonAlignment::gradientGlobalAlg(Rt, Pt, Rm, Nl, Cm);
01097     GlobalTrackerMuonAlignment::gradientGlobal(Rt, Pt, Rm, Pm, Nl, Cm); 
01098 
01099     CLHEP::HepSymMatrix covLoc(4,0);
01100     for(int i=1; i<=4; i++)
01101       for(int j=1; j<=i; j++){
01102         covLoc(i,j) = (tsosMuon.localError().matrix()
01103                        + extrapolationT.localError().matrix())(i,j);
01104         //if(i != j) Cov(i,j) = 0.;
01105       }
01106 
01107     const Surface& refSurface = tsosMuon.surface(); 
01108     CLHEP::HepMatrix rotLoc (3,3,0);
01109     rotLoc(1,1) = refSurface.rotation().xx();
01110     rotLoc(1,2) = refSurface.rotation().xy();
01111     rotLoc(1,3) = refSurface.rotation().xz();
01112     
01113     rotLoc(2,1) = refSurface.rotation().yx();
01114     rotLoc(2,2) = refSurface.rotation().yy();
01115     rotLoc(2,3) = refSurface.rotation().yz();
01116     
01117     rotLoc(3,1) = refSurface.rotation().zx();
01118     rotLoc(3,2) = refSurface.rotation().zy();
01119     rotLoc(3,3) = refSurface.rotation().zz();
01120     
01121     CLHEP::HepVector posLoc(3);
01122     posLoc(1) = refSurface.position().x();
01123     posLoc(2) = refSurface.position().y();
01124     posLoc(3) = refSurface.position().z();
01125 
01126     GlobalTrackerMuonAlignment::gradientLocal(Rt, Pt, Rm, Pm, Nl, covLoc, rotLoc, posLoc, LPRm);
01127 
01128     if(debug_){
01129       std::cout<<" Norm   "<<Nl<<std::endl;
01130       std::cout<<" posLoc "<<posLoc.T()<<std::endl;
01131       std::cout<<" rotLoc "<<rotLoc<<std::endl;
01132     }
01133 
01134     // -----------------------------------------------------  fill histogram 
01135     histo->Fill(itMuon->track()->pt()); 
01136     
01137     //histo2->Fill(itMuon->track()->outerP()); 
01138     histo2->Fill(Pt.mag()); 
01139     histo3->Fill((PI/2.-itMuon->track()->outerTheta())); 
01140     histo4->Fill(itMuon->track()->phi()); 
01141     histo5->Fill(Rmuon); 
01142     histo6->Fill(Zmuon); 
01143     histo7->Fill(RelMomResidual); 
01144     //histo8->Fill(chi); 
01145     histo8->Fill(chi_d); 
01146     
01147     histo101->Fill(Zmuon, Rmuon); 
01148     histo101->Fill(Rt0.z(), Rt0.perp()); 
01149     histo102->Fill(Rt0.x(), Rt0.y()); 
01150     histo102->Fill(Rm.x(), Rm.y()); 
01151     
01152     histo11->Fill(resR.mag()); 
01153     if(fabs(Nl.x()) < 0.98) histo12->Fill(resR.x()); 
01154     if(fabs(Nl.y()) < 0.98) histo13->Fill(resR.y()); 
01155     if(fabs(Nl.z()) < 0.98) histo14->Fill(resR.z()); 
01156     histo15->Fill(resP.x()); 
01157     histo16->Fill(resP.y()); 
01158     histo17->Fill(resP.z()); 
01159     
01160     if((fabs(Nl.x()) < 0.98) && (fabs(Nl.y()) < 0.98) &&(fabs(Nl.z()) < 0.98))
01161       { 
01162         histo18->Fill(sqrt(C0(0,0))); 
01163         histo19->Fill(sqrt(C1(0,0))); 
01164         histo20->Fill(sqrt(C1(0,0)+Ce(0,0)));              
01165       }
01166     if(fabs(Nl.x()) < 0.98) histo21->Fill(Vm(0)/sqrt(Cm(0,0))); 
01167     if(fabs(Nl.y()) < 0.98) histo22->Fill(Vm(1)/sqrt(Cm(1,1))); 
01168     if(fabs(Nl.z()) < 0.98) histo23->Fill(Vm(2)/sqrt(Cm(2,2))); 
01169     histo24->Fill(Vm(3)/sqrt(C1(3,3)+Ce(3,3))); 
01170     histo25->Fill(Vm(4)/sqrt(C1(4,4)+Ce(4,4))); 
01171     histo26->Fill(Vm(5)/sqrt(C1(5,5)+Ce(5,5))); 
01172     histo27->Fill(Nl.x()); 
01173     histo28->Fill(Nl.y()); 
01174     histo29->Fill(lenghtTrack); 
01175     histo30->Fill(lenghtMuon);     
01176     histo31->Fill(chi_Loc);
01177     histo32->Fill(Vml(1)/sqrt(Cml(1,1))); 
01178     histo33->Fill(Vml(2)/sqrt(Cml(2,2))); 
01179     histo34->Fill(Vml(3)/sqrt(Cml(3,3))); 
01180     histo35->Fill(Vml(4)/sqrt(Cml(4,4))); 
01181 
01182     if (debug_) {   //--------------------------------- debug print ----------
01183       
01184       std::cout<<" diag 0[ "<<C0(0,0)<<" "<<C0(1,1)<<" "<<C0(2,2)<<" "<<C0(3,3)<<" "
01185                <<C0(4,4)<<" "<<C0(5,5)<<" ]"<<std::endl;
01186       std::cout<<" diag e[ "<<Ce(0,0)<<" "<<Ce(1,1)<<" "<<Ce(2,2)<<" "<<Ce(3,3)<<" "
01187                <<Ce(4,4)<<" "<<Ce(5,5)<<" ]"<<std::endl;
01188       std::cout<<" diag 1[ "<<C1(0,0)<<" "<<C1(1,1)<<" "<<C1(2,2)<<" "<<C1(3,3)<<" "
01189                <<C1(4,4)<<" "<<C1(5,5)<<" ]"<<std::endl;
01190       std::cout<<" Rm   "<<Rm.x()<<" "<<Rm.y()<<" "<<Rm.z()
01191                <<" Pm   "<<Pm.x()<<" "<<Pm.y()<<" "<<Pm.z()<<std::endl;
01192       std::cout<<" Rt   "<<Rt.x()<<" "<<Rt.y()<<" "<<Rt.z()
01193                <<" Pt   "<<Pt.x()<<" "<<Pt.y()<<" "<<Pt.z()<<std::endl;
01194       std::cout<<" Nl*(Rm-Rt) "<<Nl.dot(Rm-Rt)<<std::endl;
01195       std::cout<<" resR "<<resR.x()<<" "<<resR.y()<<" "<<resR.z()
01196                <<" resP "<<resP.x()<<" "<<resP.y()<<" "<<resP.z()<<std::endl;       
01197       std::cout<<" Rm-t "<<(Rm-Rt).x()<<" "<<(Rm-Rt).y()<<" "<<(Rm-Rt).z()
01198                <<" Pm-t "<<(Pm-Pt).x()<<" "<<(Pm-Pt).y()<<" "<<(Pm-Pt).z()<<std::endl;       
01199       std::cout<<" Vm   "<<Vm<<std::endl;
01200       std::cout<<" +-   "<<sqrt(Cm(0,0))<<" "<<sqrt(Cm(1,1))<<" "<<sqrt(Cm(2,2))<<" "
01201                <<sqrt(Cm(3,3))<<" "<<sqrt(Cm(4,4))<<" "<<sqrt(Cm(5,5))<<std::endl;
01202       std::cout<<" Pmuon Ptrack dP/Ptrack "<<itMuon->outerTrack()->p()<<" "
01203                <<itMuon->track()->outerP()<<" "
01204                <<(itMuon->outerTrack()->p() - 
01205                   itMuon->track()->outerP())/itMuon->track()->outerP()<<std::endl;
01206       std::cout<<" cov matrix "<<std::endl;
01207       std::cout<<Cm<<std::endl;
01208       std::cout<<" diag [ "<<Cm(0,0)<<" "<<Cm(1,1)<<" "<<Cm(2,2)<<" "<<Cm(3,3)<<" "
01209                <<Cm(4,4)<<" "<<Cm(5,5)<<" ]"<<std::endl;
01210       
01211       static AlgebraicSymMatrix66 Ro;
01212       double Diag[6];
01213       for(int i=0; i<=5; i++) Diag[i] = sqrt(Cm(i,i));
01214       for(int i=0; i<=5; i++)
01215         for(int j=0; j<=5; j++)
01216           Ro(i,j) = Cm(i,j)/Diag[i]/Diag[j];
01217       std::cout<<" correlation matrix "<<std::endl;
01218       std::cout<<Ro<<std::endl;
01219       
01220       AlgebraicSymMatrix66 CmI;
01221       for(int i=0; i<=5; i++)
01222         for(int j=0; j<=5; j++)
01223           CmI(i,j) = Cm(i,j);
01224       
01225       bool ierr = !CmI.Invert();
01226       if( ierr ) { if(alarm || debug_)
01227           std::cout<<" Error inverse covariance matrix !!!!!!!!!!!"<<std::endl;
01228         continue;
01229       }       
01230       std::cout<<" inverse cov matrix "<<std::endl;
01231       std::cout<<Cm<<std::endl;
01232       
01233       double chi = ROOT::Math::Similarity(Vm, CmI);
01234       std::cout<<" chi chi_d "<<chi<<" "<<chi_d<<std::endl;
01235     }  // end of debug_ printout  --------------------------------------------
01236     
01237   } // end loop on selected muons, i.e. Jim's globalMuon
01238 
01239 } //end of analyzeTrackTrack
01240 
01241 
01242 // ------- method to analyze recoTrack & trajectoryMuon of Global Muon ------
01243 void GlobalTrackerMuonAlignment::analyzeTrackTrajectory
01244 (const edm::Event& iEvent, const edm::EventSetup& iSetup)
01245 {
01246   using namespace edm;
01247   using namespace reco;  
01248   //debug_ = true;
01249   bool info = false;
01250   bool alarm = false;
01251   //bool alarm = true;
01252   double PI = 3.1415927;
01253  
01254   //-*-*-*-*-*-*-
01255   cosmicMuonMode_ = true; // true: both Cosmic and IsolatedMuon are treated with 1,2 tracks
01256   //cosmicMuonMode_ = false; // true: both Cosmic and IsolatedMuon are treated with 1,2 tracks
01257   //-*-*-*-*-*-*-
01258   isolatedMuonMode_ = !cosmicMuonMode_; //true: only IsolatedMuon are treated with 1 track
01259  
01260   N_event++;
01261   if (info || debug_) {
01262     std::cout << "----- event " << N_event << " -- tracks "<<N_track<<" ---";
01263     if (cosmicMuonMode_) std::cout << " treated as CosmicMu ";
01264     if (isolatedMuonMode_) std::cout <<" treated as IsolatedMu ";
01265     std::cout << std::endl;
01266   }
01267 
01268   Handle<reco::TrackCollection> tracks;
01269   Handle<reco::TrackCollection> muons;
01270   Handle<reco::TrackCollection> gmuons;
01271   Handle<reco::MuonCollection> smuons;
01272 
01273   if (collectionIsolated){
01274     iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "TrackerOnly",  tracks);
01275     iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "StandAlone",  muons);
01276     iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "GlobalMuon",  gmuons);
01277     iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "SelectedMuons",  smuons);
01278   }
01279   else if (collectionCosmic){
01280     iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "TrackerOnly",  tracks);
01281     iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "StandAlone",  muons);
01282     iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "GlobalMuon",  gmuons);
01283     iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "SelectedMuons",  smuons);
01284   }
01285   else{
01286     iEvent.getByLabel(trackTags_,tracks);
01287     iEvent.getByLabel(muonTags_,muons);
01288     iEvent.getByLabel(gmuonTags_,gmuons);
01289     iEvent.getByLabel(smuonTags_,smuons);
01290   }  
01291 
01292   if (debug_) {
01293     std::cout << " ievBunch " << iEvent.bunchCrossing()
01294               << " runN " << (int)iEvent.run() <<std::endl;
01295     std::cout << " N tracks s/amu gmu selmu "<<tracks->size() << " " <<muons->size()
01296               << " "<<gmuons->size() << " " << smuons->size() << std::endl;
01297     for (MuonCollection::const_iterator itMuon = smuons->begin();
01298         itMuon != smuons->end();
01299          ++itMuon) {
01300       std::cout << " is isolatValid Matches " << itMuon->isIsolationValid()
01301                 << " " <<itMuon->isMatchesValid() << std::endl;
01302     }
01303   }
01304 
01305   if (isolatedMuonMode_) {                // ---- Only 1 Isolated Muon --------
01306     if (tracks->size() != 1) return;  
01307     if (muons->size()  != 1) return;  
01308     if (gmuons->size() != 1) return; 
01309     if (smuons->size() != 1) return; 
01310   }
01311   
01312   if (cosmicMuonMode_){                // ---- 1,2 Cosmic Muon --------
01313     if (smuons->size() > 2) return; 
01314     if (tracks->size() != smuons->size()) return;  
01315     if (muons->size() != smuons->size()) return;  
01316     if (gmuons->size() != smuons->size()) return;  
01317   }
01318   
01319   // ok mc_isolated_mu
01320   //edm::ESHandle<TrackerGeometry> trackerGeometry;
01321   //iSetup.get<TrackerDigiGeometryRecord>().get(trackerGeometry);
01322   // ok mc_isolated_mu
01323   //edm::ESHandle<DTGeometry> dtGeometry;
01324   //iSetup.get<MuonGeometryRecord>().get(dtGeometry);
01325   // don't work
01326   //edm::ESHandle<CSCGeometry> cscGeometry;
01327   //iSetup.get<MuonGeometryRecord>().get(cscGeometry);
01328   
01329   if (watchTrackingGeometry_.check(iSetup) || !trackingGeometry_) {
01330     edm::ESHandle<GlobalTrackingGeometry> trackingGeometry;
01331     iSetup.get<GlobalTrackingGeometryRecord>().get(trackingGeometry);
01332     trackingGeometry_ = &*trackingGeometry;
01333     theTrackingGeometry = trackingGeometry;
01334   }
01335   
01336   if (watchMagneticFieldRecord_.check(iSetup) || !magneticField_) {
01337     edm::ESHandle<MagneticField> magneticField;
01338     iSetup.get<IdealMagneticFieldRecord>().get(magneticField);
01339     magneticField_ = &*magneticField;
01340   }
01341 
01342   if (watchGlobalPositionRcd_.check(iSetup) || !globalPositionRcd_) {
01343     edm::ESHandle<Alignments> globalPositionRcd;
01344     iSetup.get<GlobalPositionRcd>().get(globalPositionRcd);
01345     globalPositionRcd_ = &*globalPositionRcd;
01346     for (std::vector<AlignTransform>::const_iterator i 
01347            = globalPositionRcd_->m_align.begin();
01348          i != globalPositionRcd_->m_align.end();  ++i){
01349       if(DetId(DetId::Tracker).rawId() == i->rawId()) iteratorTrackerRcd = i;
01350       if(DetId(DetId::Muon).rawId() == i->rawId()) iteratorMuonRcd = i;
01351       if(DetId(DetId::Ecal).rawId() == i->rawId()) iteratorEcalRcd = i;
01352       if(DetId(DetId::Hcal).rawId() == i->rawId()) iteratorHcalRcd = i;
01353     }
01354     if(debug_){
01355       std::cout << "=== iteratorTrackerRcd " << iteratorTrackerRcd->rawId()<<" ====\n" 
01356                 << " translation " << iteratorTrackerRcd->translation()<<"\n" 
01357                 << " angles " << iteratorTrackerRcd->rotation().eulerAngles() << std::endl;
01358       std::cout << iteratorTrackerRcd->rotation() << std::endl;
01359       std::cout << "=== iteratorMuonRcd " << iteratorMuonRcd->rawId()<<" ====\n" 
01360                 << " translation " << iteratorMuonRcd->translation()<<"\n" 
01361                 << " angles " << iteratorMuonRcd->rotation().eulerAngles() << std::endl;
01362       std::cout << iteratorMuonRcd->rotation() << std::endl;
01363     }   
01364   } // end of GlobalPositionRcd
01365   
01366   ESHandle<Propagator> propagator;
01367   iSetup.get<TrackingComponentsRecord>().get(propagator_, propagator);
01368 
01369   SteppingHelixPropagator alongStHePr = 
01370     SteppingHelixPropagator(magneticField_, alongMomentum);  
01371   SteppingHelixPropagator oppositeStHePr = 
01372     SteppingHelixPropagator(magneticField_, oppositeToMomentum);  
01373 
01374   //double tolerance = 5.e-5;
01375   RKTestPropagator alongRKPr = 
01376     RKTestPropagator(magneticField_, alongMomentum);
01377   RKTestPropagator oppositeRKPr = 
01378     RKTestPropagator(magneticField_, oppositeToMomentum);
01379 
01380   float epsilon = 5.;
01381   SmartPropagator alongSmPr = 
01382     SmartPropagator(alongRKPr, alongStHePr, magneticField_, alongMomentum, epsilon);
01383   SmartPropagator oppositeSmPr = 
01384     SmartPropagator(oppositeRKPr, oppositeStHePr, magneticField_, oppositeToMomentum, epsilon);
01385 
01386   if(defineFitter){
01387     if(debug_)
01388       std::cout<<" ............... DEFINE FITTER ..................."<<std::endl;
01389     KFUpdator* theUpdator = new KFUpdator();
01390     //Chi2MeasurementEstimator* theEstimator = new Chi2MeasurementEstimator(30);
01391     Chi2MeasurementEstimator* theEstimator = new Chi2MeasurementEstimator(100000,100000);
01392     theFitter = new KFTrajectoryFitter(alongSmPr, 
01393                                        *theUpdator, 
01394                                        *theEstimator);
01395     theSmoother = new KFTrajectorySmoother(alongSmPr,
01396                                            *theUpdator,     
01397                                            *theEstimator);
01398     theFitterOp = new KFTrajectoryFitter(oppositeSmPr, 
01399                                          *theUpdator, 
01400                                          *theEstimator);
01401     theSmootherOp = new KFTrajectorySmoother(oppositeSmPr,
01402                                              *theUpdator,     
01403                                              *theEstimator);
01404     
01405     edm::ESHandle<TransientTrackingRecHitBuilder> builder;
01406     iSetup.get<TransientRecHitRecord>().get("WithTrackAngle",builder);
01407     this->TTRHBuilder = &(*builder);
01408     this->MuRHBuilder = new MuonTransientTrackingRecHitBuilder(theTrackingGeometry);
01409     if(debug_){
01410       std::cout<<" theTrackingGeometry.isValid() "<<theTrackingGeometry.isValid()<<std::endl;
01411       std::cout<< "get also the MuonTransientTrackingRecHitBuilder" << "\n";
01412       std::cout<< "get also the TransientTrackingRecHitBuilder" << "\n";
01413     }
01414     defineFitter = false;
01415   }
01416 
01417   // ................................................ selected/global muon
01418   //itMuon -->  Jim's globalMuon  
01419   for(MuonCollection::const_iterator itMuon = smuons->begin();
01420       itMuon != smuons->end(); ++itMuon) {
01421     
01422     if(debug_){
01423       std::cout<<" mu gM is GM Mu SaM tM "<<itMuon->isGlobalMuon()<<" "
01424                <<itMuon->isMuon()<<" "<<itMuon->isStandAloneMuon()
01425                <<" "<<itMuon->isTrackerMuon()<<" "
01426                <<std::endl;
01427     }
01428     
01429     // get information about the innermost muon hit -------------------------
01430     TransientTrack muTT(itMuon->outerTrack(), magneticField_, trackingGeometry_);
01431     TrajectoryStateOnSurface innerMuTSOS = muTT.innermostMeasurementState();
01432     TrajectoryStateOnSurface outerMuTSOS = muTT.outermostMeasurementState();
01433     
01434     // get information about the outermost tracker hit -----------------------
01435     TransientTrack trackTT(itMuon->track(), magneticField_, trackingGeometry_);
01436     TrajectoryStateOnSurface outerTrackTSOS = trackTT.outermostMeasurementState();
01437     TrajectoryStateOnSurface innerTrackTSOS = trackTT.innermostMeasurementState();
01438     
01439     GlobalPoint pointTrackIn  = innerTrackTSOS.globalPosition();
01440     GlobalPoint pointTrackOut = outerTrackTSOS.globalPosition();
01441     float lenghtTrack = (pointTrackOut-pointTrackIn).mag();
01442     GlobalPoint pointMuonIn  = innerMuTSOS.globalPosition();
01443     GlobalPoint pointMuonOut = outerMuTSOS.globalPosition();
01444     float lenghtMuon = (pointMuonOut - pointMuonIn).mag();
01445     GlobalVector momentumTrackOut = outerTrackTSOS.globalMomentum();
01446     GlobalVector momentumTrackIn = innerTrackTSOS.globalMomentum();
01447     float distanceInIn   = (pointTrackIn  - pointMuonIn).mag();
01448     float distanceInOut  = (pointTrackIn  - pointMuonOut).mag();
01449     float distanceOutIn  = (pointTrackOut - pointMuonIn).mag();
01450     float distanceOutOut = (pointTrackOut - pointMuonOut).mag();
01451 
01452     if(debug_){
01453       std::cout<<" pointTrackIn "<<pointTrackIn<<std::endl;
01454       std::cout<<"          Out "<<pointTrackOut<<" lenght "<<lenghtTrack<<std::endl;
01455       std::cout<<" point MuonIn "<<pointMuonIn<<std::endl;
01456       std::cout<<"          Out "<<pointMuonOut<<" lenght "<<lenghtMuon<<std::endl;
01457       std::cout<<" momeTrackIn Out "<<momentumTrackIn<<" "<<momentumTrackOut<<std::endl;
01458       std::cout<<" doi io ii oo "<<distanceOutIn<<" "<<distanceInOut<<" "
01459                <<distanceInIn<<" "<<distanceOutOut<<std::endl; 
01460   }   
01461 
01462     if(lenghtTrack < 90.) continue;
01463     if(lenghtMuon < 300.) continue;
01464     if(innerMuTSOS.globalMomentum().mag() < 5. || outerMuTSOS.globalMomentum().mag() < 5.) continue;
01465     if(momentumTrackIn.mag() < 15. || momentumTrackOut.mag() < 15.) continue;
01466     if(trackTT.charge() != muTT.charge()) continue;
01467 
01468     if(debug_)
01469       std::cout<<" passed lenght momentum cuts"<<std::endl;
01470 
01471     GlobalVector GRm, GPm, Nl, Rm, Pm, Rt, Pt, Rt0;
01472     AlgebraicSymMatrix66 Cm, C0, Ce, C1;
01473 
01474     //extrapolationT = propagator->propagate(outerTrackTSOS, refSurface);
01475     TrajectoryStateOnSurface extrapolationT;
01476     int tsosMuonIf = 0;
01477 
01478     TrajectoryStateOnSurface muonFittedTSOS;
01479     TrajectoryStateOnSurface trackFittedTSOS;
01480 
01481     if( isolatedMuonMode_ ){      //------------------------------- Isolated Muon --- Out-In --
01482       if(debug_) std::cout<<" ------ Isolated (out-in) ----- "<<std::endl;
01483       const Surface& refSurface = innerMuTSOS.surface(); 
01484       ReferenceCountingPointer<TangentPlane> 
01485         tpMuLocal(refSurface.tangentPlane(innerMuTSOS.localPosition()));
01486       Nl = tpMuLocal->normalVector();
01487       ReferenceCountingPointer<TangentPlane> 
01488       tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
01489       GlobalVector Ng = tpMuGlobal->normalVector();
01490       Surface* surf = (Surface*)&refSurface;
01491       const Plane* refPlane = dynamic_cast<Plane*>(surf); 
01492       GlobalVector Nlp = refPlane->normalVector();
01493       
01494       if(debug_){
01495         std::cout<<" Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
01496          <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
01497         std::cout<<"  global "<<Ng.x()<<" "<<Ng.y()<<" "<<Ng.z()<<std::endl;
01498         std::cout<<"      lp "<<Nlp.x()<<" "<<Nlp.y()<<" "<<Nlp.z()<<std::endl;
01499         //std::cout<<" Nlocal Nglobal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
01500         // <<Ng.x()<<" "<<Ng.y()<<" "<<Ng.z()
01501         //<<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
01502       }
01503 
01504       //                          extrapolation to innermost muon hit     
01505 
01506       //extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
01507       if(!refitTrack_)
01508         extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
01509       else{
01510         GlobalTrackerMuonAlignment::trackFitter(itMuon->track(),
01511                                                 trackTT,alongMomentum,trackFittedTSOS);
01512         if(trackFittedTSOS.isValid())
01513           extrapolationT = alongSmPr.propagate(trackFittedTSOS, refSurface);       
01514       }   
01515       
01516       if(!extrapolationT.isValid()){
01517         if(false & alarm)
01518           std::cout<<" !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
01519             //<<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
01520                    <<std::endl;
01521         continue;
01522       }
01523       tsosMuonIf = 1;
01524       Rt = GlobalVector((extrapolationT.globalPosition()).x(),
01525                         (extrapolationT.globalPosition()).y(),
01526                         (extrapolationT.globalPosition()).z());
01527       
01528       Pt = extrapolationT.globalMomentum();
01529       //                          global parameters of muon
01530       GRm = GlobalVector((innerMuTSOS.globalPosition()).x(),
01531                          (innerMuTSOS.globalPosition()).y(),
01532                          (innerMuTSOS.globalPosition()).z());
01533       GPm = innerMuTSOS.globalMomentum();
01534       
01535       Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
01536                          (outerTrackTSOS.globalPosition()).y(),
01537                          (outerTrackTSOS.globalPosition()).z());
01538       Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() + 
01539                                  innerMuTSOS.cartesianError().matrix());
01540       C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());                 
01541       Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix()); 
01542       C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
01543 
01544       if(refitMuon_)
01545         GlobalTrackerMuonAlignment::muonFitter(itMuon->outerTrack(),muTT,oppositeToMomentum,muonFittedTSOS);
01546 
01547     }  //                    -------------------------------   end Isolated Muon -- Out - In  ---
01548     
01549     
01550     if( cosmicMuonMode_ ){      //------------------------------- Cosmic Muon -----
01551       
01552       if((distanceOutIn <= distanceInOut) & (distanceOutIn <= distanceInIn) & 
01553          (distanceOutIn <= distanceOutOut)){             // ----- Out - In ------    
01554         if(debug_) std::cout<<" -----    Out - In -----"<<std::endl;
01555 
01556         const Surface& refSurface = innerMuTSOS.surface(); 
01557         ReferenceCountingPointer<TangentPlane> 
01558           tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
01559         Nl = tpMuGlobal->normalVector();
01560 
01561         //                          extrapolation to innermost muon hit     
01562         //extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
01563         if(!refitTrack_)
01564           extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
01565         else{
01566           GlobalTrackerMuonAlignment::trackFitter(itMuon->track(),
01567                                                   trackTT,alongMomentum,trackFittedTSOS);
01568           if(trackFittedTSOS.isValid())
01569             extrapolationT = alongSmPr.propagate(trackFittedTSOS, refSurface);     
01570         }  
01571 
01572         if(!extrapolationT.isValid()){
01573           if(false & alarm)
01574             std::cout<<" !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
01575               //<<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
01576                      <<std::endl;
01577           continue;
01578         }
01579         if(debug_)
01580           std::cout<<" extrapolationT.isValid "<<extrapolationT.isValid()<<std::endl; 
01581 
01582         tsosMuonIf = 1;
01583         Rt = GlobalVector((extrapolationT.globalPosition()).x(),
01584                           (extrapolationT.globalPosition()).y(),
01585                           (extrapolationT.globalPosition()).z());
01586         
01587         Pt = extrapolationT.globalMomentum();
01588         //                          global parameters of muon
01589         GRm = GlobalVector((innerMuTSOS.globalPosition()).x(),
01590                            (innerMuTSOS.globalPosition()).y(),
01591                            (innerMuTSOS.globalPosition()).z());
01592         GPm = innerMuTSOS.globalMomentum();
01593         Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
01594                            (outerTrackTSOS.globalPosition()).y(),
01595                            (outerTrackTSOS.globalPosition()).z());
01596         Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() + 
01597                                    innerMuTSOS.cartesianError().matrix());
01598         C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());               
01599         Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix()); 
01600         C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
01601         
01602         if(refitMuon_)
01603           GlobalTrackerMuonAlignment::muonFitter(itMuon->outerTrack(),muTT,oppositeToMomentum,muonFittedTSOS);
01604 
01605         if(false & debug_){
01606           //std::cout<<" ->propDir "<<propagator->propagationDirection()<<std::endl;
01607           PropagationDirectionChooser Chooser;
01608           std::cout<<" propDirCh "
01609                    <<Chooser.operator()(*outerTrackTSOS.freeState(), refSurface)
01610                    <<" Ch == along "<<(alongMomentum == 
01611                                        Chooser.operator()(*outerTrackTSOS.freeState(), refSurface))
01612                    <<std::endl;
01613           std::cout<<" --- Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
01614                    <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
01615         }
01616       }           //                                       enf of ---- Out - In -----
01617       
01618 
01619       else if((distanceInOut <= distanceInIn) & (distanceInOut <= distanceOutIn) & 
01620               (distanceInOut <= distanceOutOut)){            // ----- In - Out ------    
01621         if(debug_) std::cout<<" -----    In - Out -----"<<std::endl;
01622 
01623         const Surface& refSurface = outerMuTSOS.surface(); 
01624         ReferenceCountingPointer<TangentPlane> 
01625           tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
01626         Nl = tpMuGlobal->normalVector();
01627                 
01628         //                          extrapolation to outermost muon hit     
01629         //extrapolationT = oppositeSmPr.propagate(innerTrackTSOS, refSurface);
01630         if(!refitTrack_)
01631           extrapolationT = oppositeSmPr.propagate(innerTrackTSOS, refSurface);
01632         else{
01633           GlobalTrackerMuonAlignment::trackFitter(itMuon->track(),
01634                                                 trackTT,oppositeToMomentum,trackFittedTSOS);
01635         if(trackFittedTSOS.isValid())
01636           extrapolationT = oppositeSmPr.propagate(trackFittedTSOS, refSurface);    
01637       }   
01638 
01639         if(!extrapolationT.isValid()){
01640           if(false & alarm)
01641             std::cout<<" !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
01642               <<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
01643                      <<std::endl;
01644           continue;
01645         }
01646         if(debug_)
01647           std::cout<<" extrapolationT.isValid "<<extrapolationT.isValid()<<std::endl; 
01648         
01649         tsosMuonIf = 2;
01650         Rt = GlobalVector((extrapolationT.globalPosition()).x(),
01651                           (extrapolationT.globalPosition()).y(),
01652                           (extrapolationT.globalPosition()).z());
01653         
01654         Pt = extrapolationT.globalMomentum();
01655         //                          global parameters of muon
01656         GRm = GlobalVector((outerMuTSOS.globalPosition()).x(),
01657                            (outerMuTSOS.globalPosition()).y(),
01658                            (outerMuTSOS.globalPosition()).z());
01659         GPm = outerMuTSOS.globalMomentum();
01660         Rt0 = GlobalVector((innerTrackTSOS.globalPosition()).x(),
01661                            (innerTrackTSOS.globalPosition()).y(),
01662                            (innerTrackTSOS.globalPosition()).z());
01663         Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() + 
01664                                    outerMuTSOS.cartesianError().matrix());
01665         C0 = AlgebraicSymMatrix66(innerTrackTSOS.cartesianError().matrix());               
01666         Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix()); 
01667         C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
01668         
01669         if(refitMuon_)
01670           GlobalTrackerMuonAlignment::muonFitter(itMuon->outerTrack(),muTT,alongMomentum,muonFittedTSOS);
01671 
01672         if(false & debug_){
01673           //std::cout<<" ->propDir "<<propagator->propagationDirection()<<std::endl;
01674           PropagationDirectionChooser Chooser;
01675           std::cout<<" propDirCh "
01676                    <<Chooser.operator()(*innerTrackTSOS.freeState(), refSurface)
01677                    <<" Ch == oppisite "<<(oppositeToMomentum == 
01678                                           Chooser.operator()(*innerTrackTSOS.freeState(), refSurface))
01679                    <<std::endl;
01680           std::cout<<" --- Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
01681                    <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
01682         }
01683       }           //                                        enf of ---- In - Out -----
01684 
01685       else if((distanceOutOut <= distanceInOut) & (distanceOutOut <= distanceInIn) & 
01686               (distanceOutOut <= distanceOutIn)){          // ----- Out - Out ------    
01687         if(debug_) std::cout<<" -----    Out - Out -----"<<std::endl;
01688 
01689         // reject: momentum of track has opposite direction to muon track
01690         continue;
01691 
01692         const Surface& refSurface = outerMuTSOS.surface(); 
01693         ReferenceCountingPointer<TangentPlane> 
01694           tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
01695         Nl = tpMuGlobal->normalVector();
01696                 
01697         //                          extrapolation to outermost muon hit     
01698         extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
01699 
01700         if(!extrapolationT.isValid()){
01701           if(alarm)
01702             std::cout<<" !!!!!!!!! Catastrophe Out-Out extrapolationT.isValid "
01703                      <<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
01704                      <<std::endl;
01705           continue;
01706         }
01707         if(debug_)
01708           std::cout<<" extrapolationT.isValid "<<extrapolationT.isValid()<<std::endl; 
01709 
01710         tsosMuonIf = 2;
01711         Rt = GlobalVector((extrapolationT.globalPosition()).x(),
01712                           (extrapolationT.globalPosition()).y(),
01713                           (extrapolationT.globalPosition()).z());
01714         
01715         Pt = extrapolationT.globalMomentum();
01716         //                          global parameters of muon
01717         GRm = GlobalVector((outerMuTSOS.globalPosition()).x(),
01718                            (outerMuTSOS.globalPosition()).y(),
01719                            (outerMuTSOS.globalPosition()).z());
01720         GPm = outerMuTSOS.globalMomentum();
01721         Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
01722                            (outerTrackTSOS.globalPosition()).y(),
01723                            (outerTrackTSOS.globalPosition()).z());
01724         Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() + 
01725                                    outerMuTSOS.cartesianError().matrix());
01726         C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());               
01727         Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix()); 
01728         C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
01729         
01730         if(debug_){
01731           PropagationDirectionChooser Chooser;
01732           std::cout<<" propDirCh "
01733                    <<Chooser.operator()(*outerTrackTSOS.freeState(), refSurface)
01734                    <<" Ch == along "<<(alongMomentum == 
01735                                        Chooser.operator()(*outerTrackTSOS.freeState(), refSurface))
01736                    <<std::endl;
01737           std::cout<<" --- Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
01738                    <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
01739           std::cout<<"     Nornal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
01740                    <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
01741         }
01742       }           //                                       enf of ---- Out - Out -----
01743       else{
01744         if(alarm)
01745           std::cout<<" ------- !!!!!!!!!!!!!!! No proper Out - In \a\a\a\a\a\a\a"<<std::endl; 
01746         continue;
01747       }
01748       
01749     }  //                     -------------------------------   end Cosmic Muon -----
01750     
01751     if(tsosMuonIf == 0) {if(info) {std::cout<<"No tsosMuon !!!!!!"<<std::endl; continue;}}
01752     TrajectoryStateOnSurface tsosMuon;
01753     if(tsosMuonIf == 1) tsosMuon = muTT.innermostMeasurementState();
01754     else tsosMuon = muTT.outermostMeasurementState();
01755 
01756     //GlobalTrackerMuonAlignment::misalignMuon(GRm, GPm, Nl, Rt, Rm, Pm);
01757     AlgebraicVector4 LPRm;  // muon local (dx/dz, dy/dz, x, y)
01758     GlobalTrackerMuonAlignment::misalignMuonL
01759     (GRm, GPm, Nl, Rt, Rm, Pm, LPRm, extrapolationT, tsosMuon);
01760 
01761     if(refitTrack_){
01762       if(!trackFittedTSOS.isValid()){ 
01763         if(info) std::cout<<" =================  trackFittedTSOS notValid !!!!!!!! "<<std::endl;
01764         continue;}
01765       if(debug_) this->debugTrajectorySOS(" trackFittedTSOS ", trackFittedTSOS);
01766     }
01767 
01768     if(refitMuon_){
01769       if(!muonFittedTSOS.isValid()){ 
01770         if(info) std::cout<<" =================  muonFittedTSOS notValid !!!!!!!! "<<std::endl;
01771         continue;}
01772       if(debug_) this->debugTrajectorySOS(" muonFittedTSOS ", muonFittedTSOS);
01773       Rm = GlobalVector((muonFittedTSOS.globalPosition()).x(),
01774                         (muonFittedTSOS.globalPosition()).y(),
01775                         (muonFittedTSOS.globalPosition()).z());
01776       Pm = muonFittedTSOS.globalMomentum();
01777       LPRm = AlgebraicVector4(muonFittedTSOS.localParameters().vector()(1),
01778                               muonFittedTSOS.localParameters().vector()(2),
01779                               muonFittedTSOS.localParameters().vector()(3),
01780                               muonFittedTSOS.localParameters().vector()(4)); 
01781     }
01782     GlobalVector resR = Rm - Rt;
01783     GlobalVector resP0 = Pm - Pt;
01784     GlobalVector resP = Pm / Pm.mag() - Pt / Pt.mag();
01785     float RelMomResidual = (Pm.mag() - Pt.mag()) / (Pt.mag() + 1.e-6);;
01786 
01787     AlgebraicVector6 Vm;
01788     Vm(0) = resR.x();  Vm(1) = resR.y();  Vm(2) = resR.z(); 
01789     Vm(3) = resP0.x();  Vm(4) = resP0.y();  Vm(5) = resP0.z(); 
01790     float Rmuon = Rm.perp();
01791     float Zmuon = Rm.z();
01792     float alfa_x = atan2(Nl.x(),Nl.y())*57.29578;
01793     
01794     if(debug_){
01795       std::cout<<" Nx Ny Nz alfa_x "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "<<alfa_x<<std::endl;
01796       //std::cout<<" Rm "<<Rm<<std::endl<<" Rt "<<Rt<<std::endl<<" resR "<<resR<<std::endl
01797       //       <<" resP "<<resP<<" dp/p "<<RelMomResidual<<std::endl;
01798     }
01799 
01800     double chi_d = 0;
01801     for(int i=0; i<=5; i++) chi_d += Vm(i)*Vm(i)/Cm(i,i);
01802         
01803     AlgebraicVector5 Vml(tsosMuon.localParameters().vector() 
01804                        - extrapolationT.localParameters().vector());
01805     AlgebraicSymMatrix55 m(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
01806     AlgebraicSymMatrix55 Cml(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
01807     bool ierrLoc = !m.Invert();
01808     if (ierrLoc && debug_ && info) { 
01809       std::cout<< " ==== Error inverting Local covariance matrix ==== "<<std::cout;
01810       continue;}
01811     double chi_Loc = ROOT::Math::Similarity(Vml,m);
01812     if(debug_)
01813       std::cout<<" chi_Loc px/pz/err "<<chi_Loc<<" "<<Vml(1)/sqrt(Cml(1,1))<<std::endl;
01814 
01815     if(Pt.mag() < 15.) continue;
01816     if(Pm.mag() <  5.) continue;
01817  
01818     //if(Pt.mag() < 30.) continue;          // momenum cut  < 30GeV
01819     //if(Pt.mag() < 60.) continue;          // momenum cut  < 30GeV
01820     //if(Pt.mag() > 50.) continue;          //  momenum cut > 50GeV
01821     //if(Pt.mag() > 100.) continue;          //  momenum cut > 100GeV
01822     //if(trackTT.charge() < 0) continue;    // select positive charge
01823     //if(trackTT.charge() > 0) continue;    // select negative charge
01824     
01825     //if(fabs(resR.x()) > 5.) continue;     // strong cut  X 
01826     //if(fabs(resR.y()) > 5.) continue;     //             Y 
01827     //if(fabs(resR.z()) > 5.) continue;     //             Z
01828     //if(fabs(resR.mag()) > 7.5) continue;   //            dR
01829 
01830     //if(fabs(RelMomResidual) > 0.5) continue;
01831     if(fabs(resR.x()) > 20.) continue;
01832     if(fabs(resR.y()) > 20.) continue;
01833     if(fabs(resR.z()) > 20.) continue;
01834     if(fabs(resR.mag()) > 30.) continue;
01835     if(fabs(resP.x()) > 0.06) continue;
01836     if(fabs(resP.y()) > 0.06) continue;
01837     if(fabs(resP.z()) > 0.06) continue;
01838     if(chi_d > 40.) continue;
01839 
01840     //                                            select Barrel 
01841     //if(Rmuon < 400. || Rmuon > 450.) continue; 
01842     //if(Zmuon < -600. || Zmuon > 600.) continue;
01843     //if(fabs(Nl.z()) > 0.95) continue;  
01844     //MuSelect = " Barrel";
01845     //                                                  EndCap1
01846     //if(Rmuon < 120. || Rmuon > 450.) continue;
01847     //if(Zmuon < -720.) continue;
01848     //if(Zmuon > -580.) continue;
01849     //if(fabs(Nl.z()) < 0.95) continue;  
01850     //MuSelect = " EndCap1";
01851     //                                                  EndCap2
01852     //if(Rmuon < 120. || Rmuon > 450.) continue;
01853     //if(Zmuon >  720.) continue;
01854     //if(Zmuon <  580.) continue;
01855     //if(fabs(Nl.z()) < 0.95) continue;  
01856     //MuSelect = " EndCap2";
01857     //                                                 select All
01858     if(Rmuon < 120. || Rmuon > 450.) continue;  
01859     if(Zmuon < -720. || Zmuon > 720.) continue;
01860     MuSelect = " Barrel+EndCaps";
01861 
01862     if(debug_)
01863       std::cout<<" .............. passed all cuts"<<std::endl;
01864         
01865     N_track++;
01866     //                     gradient and Hessian for each track
01867     
01868     GlobalTrackerMuonAlignment::gradientGlobalAlg(Rt, Pt, Rm, Nl, Cm);
01869     GlobalTrackerMuonAlignment::gradientGlobal(Rt, Pt, Rm, Pm, Nl, Cm); 
01870 
01871     CLHEP::HepSymMatrix covLoc(4,0);
01872     for(int i=1; i<=4; i++)
01873       for(int j=1; j<=i; j++){
01874         covLoc(i,j) = (tsosMuon.localError().matrix()
01875                        + extrapolationT.localError().matrix())(i,j);
01876         //if(i != j) Cov(i,j) = 0.;
01877       }
01878 
01879     const Surface& refSurface = tsosMuon.surface(); 
01880     CLHEP::HepMatrix rotLoc (3,3,0);
01881     rotLoc(1,1) = refSurface.rotation().xx();
01882     rotLoc(1,2) = refSurface.rotation().xy();
01883     rotLoc(1,3) = refSurface.rotation().xz();
01884     
01885     rotLoc(2,1) = refSurface.rotation().yx();
01886     rotLoc(2,2) = refSurface.rotation().yy();
01887     rotLoc(2,3) = refSurface.rotation().yz();
01888     
01889     rotLoc(3,1) = refSurface.rotation().zx();
01890     rotLoc(3,2) = refSurface.rotation().zy();
01891     rotLoc(3,3) = refSurface.rotation().zz();
01892     
01893     CLHEP::HepVector posLoc(3);
01894     posLoc(1) = refSurface.position().x();
01895     posLoc(2) = refSurface.position().y();
01896     posLoc(3) = refSurface.position().z();
01897 
01898     GlobalTrackerMuonAlignment::gradientLocal(Rt, Pt, Rm, Pm, Nl, covLoc, rotLoc, posLoc, LPRm);
01899 
01900     if(debug_){
01901       std::cout<<" Norm   "<<Nl<<std::endl;
01902       std::cout<<" posLoc "<<posLoc.T()<<std::endl;
01903       std::cout<<" rotLoc "<<rotLoc<<std::endl;
01904     }
01905 
01906     // -----------------------------------------------------  fill histogram 
01907     histo->Fill(itMuon->track()->pt()); 
01908     
01909     //histo2->Fill(itMuon->track()->outerP()); 
01910     histo2->Fill(Pt.mag()); 
01911     histo3->Fill((PI/2.-itMuon->track()->outerTheta())); 
01912     histo4->Fill(itMuon->track()->phi()); 
01913     histo5->Fill(Rmuon); 
01914     histo6->Fill(Zmuon); 
01915     histo7->Fill(RelMomResidual); 
01916     //histo8->Fill(chi); 
01917     histo8->Fill(chi_d); 
01918     
01919     histo101->Fill(Zmuon, Rmuon); 
01920     histo101->Fill(Rt0.z(), Rt0.perp()); 
01921     histo102->Fill(Rt0.x(), Rt0.y()); 
01922     histo102->Fill(Rm.x(), Rm.y()); 
01923     
01924     histo11->Fill(resR.mag()); 
01925     if(fabs(Nl.x()) < 0.98) histo12->Fill(resR.x()); 
01926     if(fabs(Nl.y()) < 0.98) histo13->Fill(resR.y()); 
01927     if(fabs(Nl.z()) < 0.98) histo14->Fill(resR.z()); 
01928     histo15->Fill(resP.x()); 
01929     histo16->Fill(resP.y()); 
01930     histo17->Fill(resP.z()); 
01931     
01932     if((fabs(Nl.x()) < 0.98) && (fabs(Nl.y()) < 0.98) &&(fabs(Nl.z()) < 0.98))
01933       { 
01934         histo18->Fill(sqrt(C0(0,0))); 
01935         histo19->Fill(sqrt(C1(0,0))); 
01936         histo20->Fill(sqrt(C1(0,0)+Ce(0,0)));              
01937       }
01938     if(fabs(Nl.x()) < 0.98) histo21->Fill(Vm(0)/sqrt(Cm(0,0))); 
01939     if(fabs(Nl.y()) < 0.98) histo22->Fill(Vm(1)/sqrt(Cm(1,1))); 
01940     if(fabs(Nl.z()) < 0.98) histo23->Fill(Vm(2)/sqrt(Cm(2,2))); 
01941     histo24->Fill(Vm(3)/sqrt(C1(3,3)+Ce(3,3))); 
01942     histo25->Fill(Vm(4)/sqrt(C1(4,4)+Ce(4,4))); 
01943     histo26->Fill(Vm(5)/sqrt(C1(5,5)+Ce(5,5))); 
01944     histo27->Fill(Nl.x()); 
01945     histo28->Fill(Nl.y()); 
01946     histo29->Fill(lenghtTrack); 
01947     histo30->Fill(lenghtMuon);     
01948     histo31->Fill(chi_Loc);
01949     histo32->Fill(Vml(1)/sqrt(Cml(1,1))); 
01950     histo33->Fill(Vml(2)/sqrt(Cml(2,2))); 
01951     histo34->Fill(Vml(3)/sqrt(Cml(3,3))); 
01952     histo35->Fill(Vml(4)/sqrt(Cml(4,4))); 
01953 
01954     if (debug_) {   //--------------------------------- debug print ----------
01955       
01956       std::cout<<" diag 0[ "<<C0(0,0)<<" "<<C0(1,1)<<" "<<C0(2,2)<<" "<<C0(3,3)<<" "
01957                <<C0(4,4)<<" "<<C0(5,5)<<" ]"<<std::endl;
01958       std::cout<<" diag e[ "<<Ce(0,0)<<" "<<Ce(1,1)<<" "<<Ce(2,2)<<" "<<Ce(3,3)<<" "
01959                <<Ce(4,4)<<" "<<Ce(5,5)<<" ]"<<std::endl;
01960       std::cout<<" diag 1[ "<<C1(0,0)<<" "<<C1(1,1)<<" "<<C1(2,2)<<" "<<C1(3,3)<<" "
01961                <<C1(4,4)<<" "<<C1(5,5)<<" ]"<<std::endl;
01962       std::cout<<" Rm   "<<Rm.x()<<" "<<Rm.y()<<" "<<Rm.z()
01963                <<" Pm   "<<Pm.x()<<" "<<Pm.y()<<" "<<Pm.z()<<std::endl;
01964       std::cout<<" Rt   "<<Rt.x()<<" "<<Rt.y()<<" "<<Rt.z()
01965                <<" Pt   "<<Pt.x()<<" "<<Pt.y()<<" "<<Pt.z()<<std::endl;
01966       std::cout<<" Nl*(Rm-Rt) "<<Nl.dot(Rm-Rt)<<std::endl;
01967       std::cout<<" resR "<<resR.x()<<" "<<resR.y()<<" "<<resR.z()
01968                <<" resP "<<resP.x()<<" "<<resP.y()<<" "<<resP.z()<<std::endl;       
01969       std::cout<<" Rm-t "<<(Rm-Rt).x()<<" "<<(Rm-Rt).y()<<" "<<(Rm-Rt).z()
01970                <<" Pm-t "<<(Pm-Pt).x()<<" "<<(Pm-Pt).y()<<" "<<(Pm-Pt).z()<<std::endl;       
01971       std::cout<<" Vm   "<<Vm<<std::endl;
01972       std::cout<<" +-   "<<sqrt(Cm(0,0))<<" "<<sqrt(Cm(1,1))<<" "<<sqrt(Cm(2,2))<<" "
01973                <<sqrt(Cm(3,3))<<" "<<sqrt(Cm(4,4))<<" "<<sqrt(Cm(5,5))<<std::endl;
01974       std::cout<<" Pmuon Ptrack dP/Ptrack "<<itMuon->outerTrack()->p()<<" "
01975                <<itMuon->track()->outerP()<<" "
01976                <<(itMuon->outerTrack()->p() - 
01977                   itMuon->track()->outerP())/itMuon->track()->outerP()<<std::endl;
01978       std::cout<<" cov matrix "<<std::endl;
01979       std::cout<<Cm<<std::endl;
01980       std::cout<<" diag [ "<<Cm(0,0)<<" "<<Cm(1,1)<<" "<<Cm(2,2)<<" "<<Cm(3,3)<<" "
01981                <<Cm(4,4)<<" "<<Cm(5,5)<<" ]"<<std::endl;
01982       
01983       static AlgebraicSymMatrix66 Ro;
01984       double Diag[6];
01985       for(int i=0; i<=5; i++) Diag[i] = sqrt(Cm(i,i));
01986       for(int i=0; i<=5; i++)
01987         for(int j=0; j<=5; j++)
01988           Ro(i,j) = Cm(i,j)/Diag[i]/Diag[j];
01989       std::cout<<" correlation matrix "<<std::endl;
01990       std::cout<<Ro<<std::endl;
01991       
01992       AlgebraicSymMatrix66 CmI;
01993       for(int i=0; i<=5; i++)
01994         for(int j=0; j<=5; j++)
01995           CmI(i,j) = Cm(i,j);
01996       
01997       bool ierr = !CmI.Invert();
01998       if( ierr ) { if(alarm || debug_)
01999           std::cout<<" Error inverse covariance matrix !!!!!!!!!!!"<<std::endl;
02000         continue;
02001       }       
02002       std::cout<<" inverse cov matrix "<<std::endl;
02003       std::cout<<Cm<<std::endl;
02004       
02005       double chi = ROOT::Math::Similarity(Vm, CmI);
02006       std::cout<<" chi chi_d "<<chi<<" "<<chi_d<<std::endl;
02007     }  // end of debug_ printout  --------------------------------------------
02008     
02009   } // end loop on selected muons, i.e. Jim's globalMuon
02010 
02011 } //end of analyzeTrackTrajectory
02012 
02013 
02014 // ----  calculate gradient and Hessian matrix (algebraic) to search global shifts ------
02015 void 
02016 GlobalTrackerMuonAlignment::gradientGlobalAlg(GlobalVector& Rt, GlobalVector& Pt, 
02017                                      GlobalVector& Rm, GlobalVector& Nl, 
02018                                      AlgebraicSymMatrix66& Cm)
02019 {
02020   // ---------------------------- Calculate Information matrix and Gfree vector
02021   // Information == Hessian , Gfree == gradient of the objective function
02022   
02023   AlgebraicMatrix33 Jac; 
02024   AlgebraicVector3 Wi, R_m, R_t, P_t, Norm, dR;   
02025   
02026   R_m(0)  = Rm.x();  R_m(1)  = Rm.y();  R_m(2)  = Rm.z(); 
02027   R_t(0)  = Rt.x();  R_t(1)  = Rt.y();  R_t(2)  = Rt.z(); 
02028   P_t(0)  = Pt.x();  P_t(1)  = Pt.y();  P_t(2)  = Pt.z(); 
02029   Norm(0) = Nl.x();  Norm(1) = Nl.y();  Norm(2) = Nl.z(); 
02030   
02031   for(int i=0; i<=2; i++){
02032     if(Cm(i,i) > 1.e-20)
02033       Wi(i) = 1./Cm(i,i);
02034     else Wi(i) = 1.e-10;
02035     dR(i) = R_m(i) - R_t(i);
02036   }
02037   
02038   float PtN = P_t(0)*Norm(0) + P_t(1)*Norm(1) + P_t(2)*Norm(2);
02039   
02040   Jac(0,0) = 1. - P_t(0)*Norm(0)/PtN;
02041   Jac(0,1) =    - P_t(0)*Norm(1)/PtN;
02042   Jac(0,2) =    - P_t(0)*Norm(2)/PtN;
02043   
02044   Jac(1,0) =    - P_t(1)*Norm(0)/PtN;
02045   Jac(1,1) = 1. - P_t(1)*Norm(1)/PtN;
02046   Jac(1,2) =    - P_t(1)*Norm(2)/PtN;
02047   
02048   Jac(2,0) =    - P_t(2)*Norm(0)/PtN;
02049   Jac(2,1) =    - P_t(2)*Norm(1)/PtN;
02050   Jac(2,2) = 1. - P_t(2)*Norm(2)/PtN;
02051   
02052   AlgebraicSymMatrix33 Itr;
02053   
02054   for(int i=0; i<=2; i++) 
02055     for(int j=0; j<=2; j++){
02056       if(j < i) continue;
02057       Itr(i,j) = 0.;
02058       //std::cout<<" ij "<<i<<" "<<j<<std::endl;         
02059       for(int k=0; k<=2; k++){
02060         Itr(i,j) += Jac(k,i)*Wi(k)*Jac(k,j);}}
02061   
02062   for(int i=0; i<=2; i++) 
02063     for(int j=0; j<=2; j++){
02064       if(j < i) continue;
02065       Inf(i,j) += Itr(i,j);}
02066   
02067   AlgebraicVector3 Gtr(0., 0., 0.);
02068   for(int i=0; i<=2; i++)
02069     for(int k=0; k<=2; k++) Gtr(i) += dR(k)*Wi(k)*Jac(k,i);
02070   for(int i=0; i<=2; i++) Gfr(i) += Gtr(i); 
02071 
02072   if(debug_){
02073    std::cout<<" Wi  "<<Wi<<std::endl;     
02074    std::cout<<" N   "<<Norm<<std::endl;
02075    std::cout<<" P_t "<<P_t<<std::endl;
02076    std::cout<<" (Pt*N) "<<PtN<<std::endl;  
02077    std::cout<<" dR  "<<dR<<std::endl;
02078    std::cout<<" +/- "<<1./sqrt(Wi(0))<<" "<<1./sqrt(Wi(1))<<" "<<1./sqrt(Wi(2))
02079          <<" "<<std::endl;
02080    std::cout<<" Jacobian dr/ddx "<<std::endl;
02081    std::cout<<Jac<<std::endl;
02082    std::cout<<" G-- "<<Gtr<<std::endl;
02083    std::cout<<" Itrack "<<std::endl;
02084    std::cout<<Itr<<std::endl;
02085    std::cout<<" Gfr "<<Gfr<<std::endl;
02086    std::cout<<" -- Inf --"<<std::endl;
02087    std::cout<<Inf<<std::endl;
02088   }
02089 
02090   return;
02091 }
02092 
02093 // ----  calculate gradient and Hessian matrix in global parameters ------
02094 void 
02095 GlobalTrackerMuonAlignment::gradientGlobal(GlobalVector& GRt, GlobalVector& GPt, 
02096                                   GlobalVector& GRm, GlobalVector& GPm, 
02097                                   GlobalVector& GNorm, AlgebraicSymMatrix66 & GCov)
02098 {
02099   // we search for 6D global correction vector (d, a), where
02100   //                        3D vector of shihts d 
02101   //               3D vector of rotation angles    a   
02102 
02103   //bool alarm = true;
02104   bool alarm = false;
02105   bool info = false;
02106 
02107   int Nd = 6;     // dimension of vector of alignment pararmeters, d 
02108   
02109   //double PtMom = GPt.mag();
02110   CLHEP::HepSymMatrix w(Nd,0);
02111   for (int i=1; i <= Nd; i++)
02112     for (int j=1; j <= Nd; j++){
02113       if(j <= i ) w(i,j) = GCov(i-1, j-1);
02114       //if(i >= 3) w(i,j) /= PtMom;
02115       //if(j >= 3) w(i,j) /= PtMom;
02116       if((i == j) && (i<=3) && (GCov(i-1, j-1) < 1.e-20))  w(i,j) = 1.e20; // w=0  
02117       if(i != j) w(i,j) = 0.;                // use diaginal elements    
02118     }
02119 
02120   //GPt /= GPt.mag();
02121   //GPm /= GPm.mag();   // end of transform
02122 
02123   CLHEP::HepVector V(Nd), Rt(3), Pt(3), Rm(3), Pm(3), Norm(3);
02124   Rt(1) = GRt.x(); Rt(2) = GRt.y();  Rt(3) = GRt.z();
02125   Pt(1) = GPt.x(); Pt(2) = GPt.y();  Pt(3) = GPt.z();
02126   Rm(1) = GRm.x(); Rm(2) = GRm.y();  Rm(3) = GRm.z();
02127   Pm(1) = GPm.x(); Pm(2) = GPm.y();  Pm(3) = GPm.z();
02128    Norm(1) = GNorm.x(); Norm(2) = GNorm.y(); Norm(3) = GNorm.z(); 
02129 
02130   V = dsum(Rm - Rt, Pm - Pt) ;   //std::cout<<" V "<<V.T()<<std::endl;
02131   
02132   //double PmN = CLHEP::dot(Pm, Norm);
02133   double PmN = CLHEP_dot(Pm, Norm);
02134 
02135   CLHEP::HepMatrix Jac(Nd,Nd,0);
02136   for (int i=1; i <= 3; i++)
02137     for (int j=1; j <= 3; j++){
02138       Jac(i,j) =  Pm(i)*Norm(j) / PmN;
02139       if(i == j ) Jac(i,j) -= 1.;
02140     } 
02141 
02142   //                                            dp/da                   
02143   Jac(4,4) =     0.;   // dpx/dax
02144   Jac(5,4) = -Pm(3);   // dpy/dax
02145   Jac(6,4) =  Pm(2);   // dpz/dax
02146   Jac(4,5) =  Pm(3);   // dpx/day
02147   Jac(5,5) =     0.;   // dpy/day
02148   Jac(6,5) = -Pm(1);   // dpz/day
02149   Jac(4,6) = -Pm(2);   // dpx/daz
02150   Jac(5,6) =  Pm(1);   // dpy/daz
02151   Jac(6,6) =     0.;   // dpz/daz
02152 
02153   CLHEP::HepVector dsda(3);
02154   dsda(1) = (Norm(2)*Rm(3) - Norm(3)*Rm(2)) / PmN;
02155   dsda(2) = (Norm(3)*Rm(1) - Norm(1)*Rm(3)) / PmN;
02156   dsda(3) = (Norm(1)*Rm(2) - Norm(2)*Rm(1)) / PmN;
02157 
02158   //                                             dr/da  
02159   Jac(1,4) =          Pm(1)*dsda(1);  // drx/dax
02160   Jac(2,4) = -Rm(3) + Pm(2)*dsda(1);  // dry/dax
02161   Jac(3,4) =  Rm(2) + Pm(3)*dsda(1);  // drz/dax
02162 
02163   Jac(1,5) =  Rm(3) + Pm(1)*dsda(2);  // drx/day
02164   Jac(2,5) =          Pm(2)*dsda(2);  // dry/day
02165   Jac(3,5) = -Rm(1) + Pm(3)*dsda(2);  // drz/day
02166 
02167   Jac(1,6) = -Rm(2) + Pm(1)*dsda(3);  // drx/daz
02168   Jac(2,6) =  Rm(1) + Pm(2)*dsda(3);  // dry/daz
02169   Jac(3,6) =          Pm(3)*dsda(3);  // drz/daz
02170 
02171   CLHEP::HepSymMatrix W(Nd,0);
02172   int ierr;
02173   W = w.inverse(ierr);
02174   if(ierr != 0) { if(alarm)
02175       std::cout<<" gradientGlobal: inversion of matrix w fail "<<std::endl;
02176     return;
02177   }
02178 
02179   CLHEP::HepMatrix W_Jac(Nd,Nd,0);
02180   W_Jac = Jac.T() * W;
02181   
02182   CLHEP::HepVector grad3(Nd);
02183   grad3 = W_Jac * V;
02184 
02185   CLHEP::HepMatrix hess3(Nd,Nd);
02186   hess3 = Jac.T() * W * Jac;
02187   //hess3(4,4) = 1.e-10;  hess3(5,5) = 1.e-10;  hess3(6,6) = 1.e-10; //????????????????
02188 
02189   Grad += grad3;
02190   Hess += hess3;
02191 
02192   CLHEP::HepVector d3I = CLHEP::solve(Hess, -Grad);
02193   int iEr3I;
02194   CLHEP::HepMatrix Errd3I = Hess.inverse(iEr3I);
02195   if( iEr3I != 0) { if(alarm)
02196       std::cout<<" gradientGlobal error inverse  Hess matrix !!!!!!!!!!!"<<std::endl;
02197   }
02198 
02199   if(info || debug_){
02200     std::cout<<" dG   "<<d3I(1)<<" "<<d3I(2)<<" "<<d3I(3)<<" "<<d3I(4)<<" "
02201              <<d3I(5)<<" "<<d3I(6)<<std::endl;;
02202     std::cout<<" +-   "<<sqrt(Errd3I(1,1))<<" "<<sqrt(Errd3I(2,2))<<" "<<sqrt(Errd3I(3,3))
02203              <<" "<<sqrt(Errd3I(4,4))<<" "<<sqrt(Errd3I(5,5))<<" "<<sqrt(Errd3I(6,6))
02204              <<std::endl;
02205   }
02206 
02207 #ifdef CHECK_OF_DERIVATIVES
02208   //              --------------------    check of derivatives
02209   
02210   CLHEP::HepVector d(3,0), a(3,0); 
02211   CLHEP::HepMatrix T(3,3,1);
02212 
02213   CLHEP::HepMatrix Ti = T.T();
02214   //double A = CLHEP::dot(Ti*Pm, Norm);
02215   //double B = CLHEP::dot((Rt -Ti*Rm + Ti*d), Norm);
02216   double A = CLHEP_dot(Ti*Pm, Norm);
02217   double B = CLHEP_dot((Rt -Ti*Rm + Ti*d), Norm);
02218   double s0 = B / A;
02219 
02220   CLHEP::HepVector r0(3,0), p0(3,0);
02221   r0 = Ti*Rm - Ti*d + s0*(Ti*Pm) - Rt;
02222   p0 = Ti*Pm - Pt;
02223 
02224   double delta = 0.0001;
02225 
02226   int ii = 3; d(ii) += delta; // d
02227   //T(2,3) += delta;   T(3,2) -= delta;  int ii = 1; // a1
02228   //T(3,1) += delta;   T(1,3) -= delta;   int ii = 2; // a2
02229   //T(1,2) += delta;   T(2,1) -= delta;   int ii = 3; // a2
02230   Ti = T.T();
02231   //A = CLHEP::dot(Ti*Pm, Norm);
02232   //B = CLHEP::dot((Rt -Ti*Rm + Ti*d), Norm);
02233   A = CLHEP_dot(Ti*Pm, Norm);
02234   B = CLHEP_dot((Rt -Ti*Rm + Ti*d), Norm);
02235   double s = B / A;
02236 
02237   CLHEP::HepVector r(3,0), p(3,0);
02238   r = Ti*Rm - Ti*d + s*(Ti*Pm) - Rt;
02239   p = Ti*Pm - Pt;
02240 
02241   std::cout<<" s0 s num dsda("<<ii<<") "<<s0<<" "<<s<<" "
02242            <<(s-s0)/delta<<" "<<dsda(ii)<<endl;
02243   // d(r,p) / d shift
02244   std::cout<<" -- An d(r,p)/d("<<ii<<") "<<Jac(1,ii)<<" "<<Jac(2,ii)<<" "
02245      <<Jac(3,ii)<<" "<<Jac(4,ii)<<" "<<Jac(5,ii)<<" "
02246      <<Jac(6,ii)<<std::endl;
02247   std::cout<<"    Nu d(r,p)/d("<<ii<<") "<<(r(1)-r0(1))/delta<<" "
02248      <<(r(2)-r0(2))/delta<<" "<<(r(3)-r0(3))/delta<<" "
02249      <<(p(1)-p0(1))/delta<<" "<<(p(2)-p0(2))/delta<<" "
02250      <<(p(3)-p0(3))/delta<<std::endl;
02251   // d(r,p) / d angle
02252   std::cout<<" -- An d(r,p)/a("<<ii<<") "<<Jac(1,ii+3)<<" "<<Jac(2,ii+3)<<" "
02253            <<Jac(3,ii+3)<<" "<<Jac(4,ii+3)<<" "<<Jac(5,ii+3)<<" "
02254            <<Jac(6,ii+3)<<std::endl;
02255   std::cout<<"    Nu d(r,p)/a("<<ii<<") "<<(r(1)-r0(1))/delta<<" "
02256            <<(r(2)-r0(2))/delta<<" "<<(r(3)-r0(3))/delta<<" "
02257            <<(p(1)-p0(1))/delta<<" "<<(p(2)-p0(2))/delta<<" "
02258            <<(p(3)-p0(3))/delta<<std::endl;
02259   //               -----------------------------  end of check
02260 #endif
02261 
02262   return;
02263 } // end gradientGlobal
02264 
02265 
02266 // ----  calculate gradient and Hessian matrix in local parameters ------
02267 void 
02268 GlobalTrackerMuonAlignment::gradientLocal(GlobalVector& GRt, GlobalVector& GPt, 
02269                                           GlobalVector& GRm, GlobalVector& GPm, 
02270                                           GlobalVector& GNorm, 
02271                                           CLHEP::HepSymMatrix& covLoc,
02272                                           CLHEP::HepMatrix& rotLoc,
02273                                           CLHEP::HepVector& R0,
02274                                           AlgebraicVector4& LPRm) 
02275 {
02276   // we search for 6D global correction vector (d, a), where
02277   //                        3D vector of shihts d 
02278   //               3D vector of rotation angles    a   
02279 
02280   bool alarm = true;
02281   //bool alarm = false;
02282   bool info = false;
02283 
02284   if(debug_)
02285     std::cout<<" gradientLocal "<<std::endl; 
02286 
02287   /*
02288   const Surface& refSurface = tsosMuon.surface();
02289 
02290   CLHEP::HepMatrix rotLoc (3,3,0);
02291   rotLoc(1,1) = refSurface.rotation().xx();
02292   rotLoc(1,2) = refSurface.rotation().xy();
02293   rotLoc(1,3) = refSurface.rotation().xz();
02294   
02295   rotLoc(2,1) = refSurface.rotation().yx();
02296   rotLoc(2,2) = refSurface.rotation().yy();
02297   rotLoc(2,3) = refSurface.rotation().yz();
02298   
02299   rotLoc(3,1) = refSurface.rotation().zx();
02300   rotLoc(3,2) = refSurface.rotation().zy();
02301   rotLoc(3,3) = refSurface.rotation().zz();
02302   */  
02303 
02304   CLHEP::HepVector Rt(3), Pt(3), Rm(3), Pm(3), Norm(3);
02305   Rt(1) = GRt.x(); Rt(2) = GRt.y();  Rt(3) = GRt.z();
02306   Pt(1) = GPt.x(); Pt(2) = GPt.y();  Pt(3) = GPt.z();
02307   Rm(1) = GRm.x(); Rm(2) = GRm.y();  Rm(3) = GRm.z();
02308   Pm(1) = GPm.x(); Pm(2) = GPm.y();  Pm(3) = GPm.z();
02309   Norm(1) = GNorm.x(); Norm(2) = GNorm.y(); Norm(3) = GNorm.z(); 
02310 
02311   CLHEP::HepVector V(4), Rml(3), Pml(3), Rtl(3), Ptl(3);
02312 
02313   /*
02314   R0(1) = refSurface.position().x();
02315   R0(2) = refSurface.position().y();
02316   R0(3) = refSurface.position().z();
02317   */ 
02318 
02319   Rml = rotLoc * (Rm - R0);
02320   Rtl = rotLoc * (Rt - R0);
02321   Pml = rotLoc * Pm;
02322   Ptl = rotLoc * Pt;
02323   
02324   V(1) = LPRm(0) - Ptl(1)/Ptl(3); 
02325   V(2) = LPRm(1) - Ptl(2)/Ptl(3);   
02326   V(3) = LPRm(2) - Rtl(1);   
02327   V(4) = LPRm(3) - Rtl(2); 
02328 
02329   /*
02330   CLHEP::HepSymMatrix Cov(4,0), W(4,0);
02331   for(int i=1; i<=4; i++)
02332     for(int j=1; j<=i; j++){
02333       Cov(i,j) = (tsosMuon.localError().matrix()
02334                   + tsosTrack.localError().matrix())(i,j);
02335       //if(i != j) Cov(i,j) = 0.;
02336       //if((i == j) && ((i==1) || (i==2))) Cov(i,j) = 100.;
02337       //if((i == j) && ((i==3) || (i==4))) Cov(i,j) = 10000.;
02338     }
02339   W = Cov;
02340   */
02341 
02342   CLHEP::HepSymMatrix W = covLoc;
02343 
02344   int ierr;
02345   W.invert(ierr);
02346   if(ierr != 0) { if(alarm)
02347       std::cout<<" gradientLocal: inversion of matrix W fail "<<std::endl;
02348     return;
02349   }
02350   
02351   //                               JacobianCartesianToLocal
02352  
02353   //AlgebraicMatrix56 jacobian   // differ from calculation above 
02354   //= JacobianCartesianToLocal::JacobianCartesianToLocal 
02355   //(refSurface, tsosTrack.localParameters()).jacobian();  
02356   //for(int i=1; i<=4; i++) for(int j=1; j<=6; j++){
02357   //int j1 = j - 1; JacToLoc(i,j) = jacobian(i, j1);}
02358   
02359   CLHEP::HepMatrix JacToLoc(4,6,0);
02360   for(int i=1; i<=2; i++)
02361     for(int j=1; j<=3; j++){
02362       JacToLoc(i,j+3) = (rotLoc(i,j) - rotLoc(3,j)*Pml(i)/Pml(3))/Pml(3); 
02363       JacToLoc(i+2,j) = rotLoc(i,j);
02364     }
02365 
02366   //                                    JacobianCorrectionsToCartesian
02367   //double PmN = CLHEP::dot(Pm, Norm);
02368   double PmN = CLHEP_dot(Pm, Norm);
02369 
02370   CLHEP::HepMatrix Jac(6,6,0);
02371   for (int i=1; i <= 3; i++)
02372     for (int j=1; j <= 3; j++){
02373       Jac(i,j) =  Pm(i)*Norm(j) / PmN;
02374       if(i == j ) Jac(i,j) -= 1.;
02375     }  
02376 
02377   //                                            dp/da                   
02378   Jac(4,4) =     0.;   // dpx/dax
02379   Jac(5,4) = -Pm(3);   // dpy/dax
02380   Jac(6,4) =  Pm(2);   // dpz/dax
02381   Jac(4,5) =  Pm(3);   // dpx/day
02382   Jac(5,5) =     0.;   // dpy/day
02383   Jac(6,5) = -Pm(1);   // dpz/day
02384   Jac(4,6) = -Pm(2);   // dpx/daz
02385   Jac(5,6) =  Pm(1);   // dpy/daz
02386   Jac(6,6) =     0.;   // dpz/daz
02387 
02388   CLHEP::HepVector dsda(3);
02389   dsda(1) = (Norm(2)*Rm(3) - Norm(3)*Rm(2)) / PmN;
02390   dsda(2) = (Norm(3)*Rm(1) - Norm(1)*Rm(3)) / PmN;
02391   dsda(3) = (Norm(1)*Rm(2) - Norm(2)*Rm(1)) / PmN;
02392 
02393   //                                             dr/da  
02394   Jac(1,4) =          Pm(1)*dsda(1);  // drx/dax
02395   Jac(2,4) = -Rm(3) + Pm(2)*dsda(1);  // dry/dax
02396   Jac(3,4) =  Rm(2) + Pm(3)*dsda(1);  // drz/dax
02397 
02398   Jac(1,5) =  Rm(3) + Pm(1)*dsda(2);  // drx/day
02399   Jac(2,5) =          Pm(2)*dsda(2);  // dry/day
02400   Jac(3,5) = -Rm(1) + Pm(3)*dsda(2);  // drz/day
02401 
02402   Jac(1,6) = -Rm(2) + Pm(1)*dsda(3);  // drx/daz
02403   Jac(2,6) =  Rm(1) + Pm(2)*dsda(3);  // dry/daz
02404   Jac(3,6) =          Pm(3)*dsda(3);  // drz/daz
02405 
02406   //                                   JacobianCorrectionToLocal
02407   CLHEP::HepMatrix JacCorLoc(4,6,0);
02408   JacCorLoc = JacToLoc * Jac;
02409 
02410   //                                   gradient and Hessian 
02411   CLHEP::HepMatrix W_Jac(6,4,0);
02412   W_Jac = JacCorLoc.T() * W;
02413   
02414   CLHEP::HepVector gradL(6);
02415   gradL = W_Jac * V;
02416 
02417   CLHEP::HepMatrix hessL(6,6);
02418   hessL = JacCorLoc.T() * W * JacCorLoc;
02419 
02420   GradL += gradL;
02421   HessL += hessL;
02422 
02423   CLHEP::HepVector dLI = CLHEP::solve(HessL, -GradL);
02424   int iErI;
02425   CLHEP::HepMatrix ErrdLI = HessL.inverse(iErI);
02426   if( iErI != 0) { if(alarm)
02427       std::cout<<" gradLocal Error inverse  Hess matrix !!!!!!!!!!!"<<std::endl;
02428   }
02429 
02430   if(info || debug_){
02431     std::cout<<" dL   "<<dLI(1)<<" "<<dLI(2)<<" "<<dLI(3)<<" "<<dLI(4)<<" "
02432              <<dLI(5)<<" "<<dLI(6)<<std::endl;;
02433     std::cout<<" +-   "<<sqrt(ErrdLI(1,1))<<" "<<sqrt(ErrdLI(2,2))<<" "<<sqrt(ErrdLI(3,3))
02434              <<" "<<sqrt(ErrdLI(4,4))<<" "<<sqrt(ErrdLI(5,5))<<" "<<sqrt(ErrdLI(6,6))
02435              <<std::endl;
02436   }
02437  
02438   if(debug_){
02439     //std::cout<<" dV(da3) {"<<-JacCorLoc(1,6)*0.002<<" "<<-JacCorLoc(2,6)*0.002<<" "
02440     //   <<-JacCorLoc(3,6)*0.002<<" "<<-JacCorLoc(4,6)*0.002<<"}"<<std::endl;
02441     //std::cout<<" JacCLo {"<<JacCorLoc(1,6)<<" "<<JacCorLoc(2,6)<<" "
02442     //   <<JacCorLoc(3,6)<<" "<<JacCorLoc(4,6)<<"}"<<std::endl;
02443     //std::cout<<"Jpx/yx  {"<<Jac(4,6)/Pm(3)<<" "<<Jac(5,6)/Pm(3)<<" "
02444     //   <<Jac(1,6)<<" "<<Jac(2,6)<<"}"<<std::endl;
02445     //std::cout<<"Jac(,a3){"<<Jac(1,6)<<" "<<Jac(2,6)<<" "<<Jac(3,6)<<" "<<Jac(4,6)<<" "
02446     //   <<Jac(5,6)<<" "<<Jac(6,6)<<std::endl;
02447     int i=5;  
02448     if(GNorm.z() > 0.95)
02449       std::cout<<" Ecap1   N "<<GNorm<<std::endl; 
02450     else if(GNorm.z() < -0.95)
02451       std::cout<<" Ecap2   N "<<GNorm<<std::endl; 
02452     else
02453       std::cout<<" Barrel  N "<<GNorm<<std::endl; 
02454     std::cout<<" JacCLo(i,"<<i<<") = {"<<JacCorLoc(1,i)<<" "<<JacCorLoc(2,i)<<" "
02455      <<JacCorLoc(3,i)<<" "<<JacCorLoc(4,i)<<"}"<<std::endl;
02456     std::cout<<"   rotLoc "<<rotLoc<<std::endl;
02457     std::cout<<" position "<<R0<<std::endl;
02458     std::cout<<" Pm,l "<<Pm.T()<<" "<<Pml(1)/Pml(3)<<" "<<Pml(2)/Pml(3)<<std::endl;
02459     std::cout<<" Pt,l "<<Pt.T()<<" "<<Ptl(1)/Ptl(3)<<" "<<Ptl(2)/Ptl(3)<<std::endl;
02460     std::cout<<" V  "<<V.T()<<std::endl;  
02461     std::cout<<" Cov \n"<<covLoc<<std::endl;
02462     std::cout<<" W*Cov "<<W*covLoc<<std::endl;
02463     //std::cout<<" JacCarToLoc = drldrc \n"<<
02464     //JacobianCartesianToLocal::JacobianCartesianToLocal
02465     //(refSurface, tsosTrack.localParameters()).jacobian()<<std::endl;
02466     std::cout<<" JacToLoc "<<JacToLoc<<std::endl;
02467     
02468   }
02469 
02470 #ifdef CHECK_OF_JACOBIAN_CARTESIAN_TO_LOCAL
02471   //---------------------- check of derivatives
02472   CLHEP::HepVector V0(4,0);
02473   V0(1) = Pml(1)/Pml(3) - Ptl(1)/Ptl(3); 
02474   V0(2) = Pml(2)/Pml(3) - Ptl(2)/Ptl(3);   
02475   V0(3) = Rml(1) - Rtl(1);   
02476   V0(4) = Rml(2) - Rtl(2); 
02477   int ii = 3; float delta = 0.01;
02478   CLHEP::HepVector V1(4,0);
02479   if(ii <= 3) {Rm(ii) += delta; Rml = rotLoc * (Rm - R0);}
02480   else {Pm(ii-3) += delta; Pml = rotLoc * Pm;}
02481   //if(ii <= 3) {Rt(ii) += delta; Rtl = rotLoc * (Rt - R0);}
02482   //else {Pt(ii-3) += delta; Ptl = rotLoc * Pt;}
02483   V1(1) = Pml(1)/Pml(3) - Ptl(1)/Ptl(3); 
02484   V1(2) = Pml(2)/Pml(3) - Ptl(2)/Ptl(3);   
02485   V1(3) = Rml(1) - Rtl(1);   
02486   V1(4) = Rml(2) - Rtl(2); 
02487 
02488   if(GNorm.z() > 0.95)
02489     std::cout<<" Ecap1   N "<<GNorm<<std::endl; 
02490   else if(GNorm.z() < -0.95)
02491     std::cout<<" Ecap2   N "<<GNorm<<std::endl; 
02492   else
02493     std::cout<<" Barrel  N "<<GNorm<<std::endl; 
02494   std::cout<<" dldc Num (i,"<<ii<<") "<<(V1(1)-V0(1))/delta<<" "<<(V1(2)-V0(2))/delta<<" "
02495            <<(V1(3)-V0(3))/delta<<" "<<(V1(4)-V0(4))/delta<<std::endl;
02496   std::cout<<" dldc Ana (i,"<<ii<<") "<<JacToLoc(1,ii)<<" "<<JacToLoc(2,ii)<<" "
02497            <<JacToLoc(3,ii)<<" "<<JacToLoc(4,ii)<<std::endl;
02498   //float dtxdpx = (rotLoc(1,1) - rotLoc(3,1)*Pml(1)/Pml(3))/Pml(3);
02499   //float dtydpx = (rotLoc(2,1) - rotLoc(3,2)*Pml(2)/Pml(3))/Pml(3);
02500   //std::cout<<" dtx/dpx dty/ "<<dtxdpx<<" "<<dtydpx<<std::endl;
02501 #endif                                            
02502 
02503   return;
02504 } // end gradientLocal
02505 
02506 
02507 // ----  simulate a misalignment of muon system   ------
02508 void 
02509 GlobalTrackerMuonAlignment::misalignMuon(GlobalVector& GRm, GlobalVector& GPm, GlobalVector& Nl, 
02510                                 GlobalVector& Rt, GlobalVector& Rm, GlobalVector& Pm)
02511 {
02512   CLHEP::HepVector d(3,0), a(3,0), Norm(3), Rm0(3), Pm0(3), Rm1(3), Pm1(3), 
02513     Rmi(3), Pmi(3) ; 
02514   
02515   d(1)=0.0; d(2)=0.0; d(3)=0.0; a(1)=0.0000; a(2)=0.0000; a(3)=0.0000; // zero     
02516   //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
02517   //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
02518   //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    
02519   //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     
02520   //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    
02521   //d(1)=10.; d(2)=20.; d(3)=30.; a(1)=0.0100; a(2)=0.0200; a(3)=0.0300; // huge      
02522   //d(1)=10.; d(2)=20.; d(3)=30.; a(1)=0.2000; a(2)=0.2500; a(3)=0.3000; // huge angles      
02523   //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
02524   //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
02525   //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
02526 
02527   MuGlShift = d; MuGlAngle = a;
02528 
02529   if((d(1) == 0.) & (d(2) == 0.) & (d(3) == 0.) & 
02530      (a(1) == 0.) & (a(2) == 0.) & (a(3) == 0.)){
02531     Rm = GRm;
02532     Pm = GPm;
02533     if(debug_)
02534       std::cout<<" ......   without misalignment "<<std::endl;
02535     return;
02536   }
02537   
02538   Rm0(1) = GRm.x(); Rm0(2) = GRm.y(); Rm0(3) = GRm.z();
02539   Pm0(1) = GPm.x(); Pm0(2) = GPm.y(); Pm0(3) = GPm.z();  
02540   Norm(1) = Nl.x(); Norm(2) = Nl.y(); Norm(3) = Nl.z();     
02541   
02542   CLHEP::HepMatrix T(3,3,1);
02543  
02544   //T(1,2) =  a(3); T(1,3) = -a(2);
02545   //T(2,1) = -a(3); T(2,3) =  a(1);
02546   //T(3,1) =  a(2); T(3,2) = -a(1);
02547   
02548   double s1 = std::sin(a(1)), c1 = std::cos(a(1));
02549   double s2 = std::sin(a(2)), c2 = std::cos(a(2));
02550   double s3 = std::sin(a(3)), c3 = std::cos(a(3));     
02551   T(1,1) =  c2 * c3;  
02552   T(1,2) =  c1 * s3 + s1 * s2 * c3;
02553   T(1,3) =  s1 * s3 - c1 * s2 * c3;
02554   T(2,1) = -c2 * s3; 
02555   T(2,2) =  c1 * c3 - s1 * s2 * s3;
02556   T(2,3) =  s1 * c3 + c1 * s2 * s3;
02557   T(3,1) =  s2;
02558   T(3,2) = -s1 * c2;
02559   T(3,3) =  c1 * c2;
02560   
02561   //int terr;
02562   //CLHEP::HepMatrix Ti = T.inverse(terr);
02563   CLHEP::HepMatrix Ti = T.T();
02564   CLHEP::HepVector di = -Ti * d;
02565   
02566   CLHEP::HepVector ex0(3,0), ey0(3,0), ez0(3,0), ex(3,0), ey(3,0), ez(3,0);
02567   ex0(1) = 1.; ey0(2) = 1.;  ez0(3) = 1.;
02568   ex = Ti*ex0; ey = Ti*ey0; ez = Ti*ez0;
02569   
02570   CLHEP::HepVector TiN = Ti * Norm; 
02571   //double si = CLHEP::dot((Ti*Rm0 - Rm0 + di), TiN) / CLHEP::dot(Pm0, TiN); 
02572   //Rm1(1) = CLHEP::dot(ex, Rm0 + si*Pm0 - di);
02573   //Rm1(2) = CLHEP::dot(ey, Rm0 + si*Pm0 - di);
02574   //Rm1(3) = CLHEP::dot(ez, Rm0 + si*Pm0 - di);
02575   double si = CLHEP_dot((Ti*Rm0 - Rm0 + di), TiN) / CLHEP_dot(Pm0, TiN); 
02576   Rm1(1) = CLHEP_dot(ex, Rm0 + si*Pm0 - di);
02577   Rm1(2) = CLHEP_dot(ey, Rm0 + si*Pm0 - di);
02578   Rm1(3) = CLHEP_dot(ez, Rm0 + si*Pm0 - di);
02579   Pm1 = T * Pm0;
02580   
02581   Rm = GlobalVector(Rm1(1), Rm1(2), Rm1(3));  
02582   //Pm = GlobalVector(CLHEP::dot(Pm0,ex), CLHEP::dot(Pm0, ey), CLHEP::dot(Pm0,ez));// =T*Pm0 
02583   Pm = GlobalVector(CLHEP_dot(Pm0,ex), CLHEP_dot(Pm0, ey), CLHEP_dot(Pm0,ez));// =T*Pm0 
02584   
02585   if(debug_){    //    -------------  debug tranformation 
02586     
02587     std::cout<<" ----- Pm "<<Pm<<std::endl;
02588     std::cout<<"    T*Pm0 "<<Pm1.T()<<std::endl;
02589     //std::cout<<" Ti*T*Pm0 "<<(Ti*(T*Pm0)).T()<<std::endl;
02590     
02591     //CLHEP::HepVector Rml = (Rm0 + si*Pm0 - di) + di;
02592     CLHEP::HepVector Rml = Rm1(1)*ex + Rm1(2)*ey + Rm1(3)*ez + di;
02593     
02594     CLHEP::HepVector Rt0(3); 
02595     Rt0(1)=Rt.x();  Rt0(2)=Rt.y();  Rt0(3)=Rt.z();      
02596     
02597     //double sl = CLHEP::dot(T*(Rt0 - Rml), T*Norm) / CLHEP::dot(Ti * Pm1, Norm);     
02598     double sl = CLHEP_dot(T*(Rt0 - Rml), T*Norm) / CLHEP_dot(Ti * Pm1, Norm);     
02599     CLHEP::HepVector Rl = Rml + sl*(Ti*Pm1);
02600     
02601     //double A = CLHEP::dot(Ti*Pm1, Norm);
02602     //double B = CLHEP::dot(Rt0, Norm) - CLHEP::dot(Ti*Rm1, Norm) 
02603     //+ CLHEP::dot(Ti*d, Norm); 
02604     double A = CLHEP_dot(Ti*Pm1, Norm);
02605     double B = CLHEP_dot(Rt0, Norm) - CLHEP_dot(Ti*Rm1, Norm) 
02606       + CLHEP_dot(Ti*d, Norm); 
02607     double s = B/A;
02608     CLHEP::HepVector Dr = Ti*Rm1 - Ti*d  + s*(Ti*Pm1) - Rm0;
02609     CLHEP::HepVector Dp = Ti*Pm1 - Pm0;
02610 
02611     CLHEP::HepVector TiR = Ti * Rm0 + di;
02612     CLHEP::HepVector Ri = T * TiR + d;
02613     
02614     std::cout<<" --TTi-0 "<<(Ri-Rm0).T()<<std::endl;
02615     std::cout<<" --  Dr  "<<Dr.T()<<endl;
02616     std::cout<<" --  Dp  "<<Dp.T()<<endl;
02617     //std::cout<<" ex "<<ex.T()<<endl;
02618     //std::cout<<" ey "<<ey.T()<<endl;
02619     //std::cout<<" ez "<<ez.T()<<endl;
02620     //std::cout<<" ---- T  ---- "<<T<<std::endl;
02621     //std::cout<<" ---- Ti ---- "<<Ti<<std::endl;
02622     //std::cout<<" ---- d  ---- "<<d.T()<<std::endl;
02623     //std::cout<<" ---- di ---- "<<di.T()<<std::endl;
02624     std::cout<<" -- si sl s "<<si<<" "<<sl<<" "<<s<<std::endl;
02625     //std::cout<<" -- si sl "<<si<<" "<<sl<<endl;
02626     //std::cout<<" -- si s "<<si<<" "<<s<<endl;
02627     std::cout<<" -- Rml-(Rm0+si*Pm0) "<<(Rml - Rm0 - si*Pm0).T()<<std::endl;
02628     //std::cout<<" -- Rm0 "<<Rm0.T()<<std::endl;
02629     //std::cout<<" -- Rm1 "<<Rm1.T()<<std::endl;
02630     //std::cout<<" -- Rmi "<<Rmi.T()<<std::endl;
02631     //std::cout<<" --siPm "<<(si*Pm0).T()<<std::endl;
02632     //std::cout<<" --s2Pm "<<(s2*(T * Pm1)).T()<<std::endl;
02633     //std::cout<<" --R1-0 "<<(Rm1-Rm0).T()<<std::endl;
02634     //std::cout<<" --Ri-0 "<<(Rmi-Rm0).T()<<std::endl;
02635     std::cout<<" --Rl-0 "<<(Rl-Rm0).T()<<std::endl;
02636     //std::cout<<" --Pi-0 "<<(Pmi-Pm0).T()<<std::endl;
02637   }      //                                --------   end of debug
02638   
02639   return;
02640 
02641 } // end of misalignMuon
02642 
02643 // ----  simulate a misalignment of muon system with Local  ------
02644 void 
02645 GlobalTrackerMuonAlignment::misalignMuonL(GlobalVector& GRm, GlobalVector& GPm, GlobalVector& Nl, 
02646                                           GlobalVector& Rt, GlobalVector& Rm, GlobalVector& Pm,
02647                                           AlgebraicVector4& Vm, 
02648                                           TrajectoryStateOnSurface& tsosTrack,
02649                                           TrajectoryStateOnSurface& tsosMuon)
02650 {
02651   CLHEP::HepVector d(3,0), a(3,0), Norm(3), Rm0(3), Pm0(3), Rm1(3), Pm1(3), 
02652     Rmi(3), Pmi(3); 
02653   
02654   d(1)=0.0; d(2)=0.0; d(3)=0.0; a(1)=0.0000; a(2)=0.0000; a(3)=0.0000; // zero     
02655   //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
02656   //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
02657   //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
02658   //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    
02659   //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    
02660   //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     
02661   //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    
02662   //d(1)=10.; d(2)=20.; d(3)=30.; a(1)=0.0100; a(2)=0.0200; a(3)=0.0300; // huge      
02663   //d(1)=10.; d(2)=20.; d(3)=30.; a(1)=0.2000; a(2)=0.2500; a(3)=0.3000; // huge angles      
02664   //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
02665   //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
02666   //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
02667 
02668   MuGlShift = d; MuGlAngle = a;
02669 
02670   if((d(1) == 0.) & (d(2) == 0.) & (d(3) == 0.) & 
02671      (a(1) == 0.) & (a(2) == 0.) & (a(3) == 0.)){
02672     Rm = GRm;
02673     Pm = GPm;
02674     AlgebraicVector4 Vm0;
02675     Vm = AlgebraicVector4(tsosMuon.localParameters().vector()(1),
02676                           tsosMuon.localParameters().vector()(2),
02677                           tsosMuon.localParameters().vector()(3),
02678                           tsosMuon.localParameters().vector()(4)); 
02679     if(debug_)
02680       std::cout<<" ......   without misalignment "<<std::endl;
02681     return;
02682   }
02683   
02684   Rm0(1) = GRm.x(); Rm0(2) = GRm.y(); Rm0(3) = GRm.z();
02685   Pm0(1) = GPm.x(); Pm0(2) = GPm.y(); Pm0(3) = GPm.z();  
02686   Norm(1) = Nl.x(); Norm(2) = Nl.y(); Norm(3) = Nl.z();     
02687   
02688   CLHEP::HepMatrix T(3,3,1);
02689  
02690   //T(1,2) =  a(3); T(1,3) = -a(2);
02691   //T(2,1) = -a(3); T(2,3) =  a(1);
02692   //T(3,1) =  a(2); T(3,2) = -a(1);
02693   
02694   double s1 = std::sin(a(1)), c1 = std::cos(a(1));
02695   double s2 = std::sin(a(2)), c2 = std::cos(a(2));
02696   double s3 = std::sin(a(3)), c3 = std::cos(a(3));     
02697   T(1,1) =  c2 * c3;  
02698   T(1,2) =  c1 * s3 + s1 * s2 * c3;
02699   T(1,3) =  s1 * s3 - c1 * s2 * c3;
02700   T(2,1) = -c2 * s3; 
02701   T(2,2) =  c1 * c3 - s1 * s2 * s3;
02702   T(2,3) =  s1 * c3 + c1 * s2 * s3;
02703   T(3,1) =  s2;
02704   T(3,2) = -s1 * c2;
02705   T(3,3) =  c1 * c2;
02706 
02707   //                    Ti, di what we have to apply for misalignment  
02708   //int terr;
02709   //CLHEP::HepMatrix Ti = T.inverse(terr);
02710   CLHEP::HepMatrix Ti = T.T();
02711   CLHEP::HepVector di = -Ti * d;
02712   
02713   CLHEP::HepVector ex0(3,0), ey0(3,0), ez0(3,0), ex(3,0), ey(3,0), ez(3,0);
02714   ex0(1) = 1.; ey0(2) = 1.;  ez0(3) = 1.;
02715   ex = Ti*ex0; ey = Ti*ey0; ez = Ti*ez0;
02716 
02717   CLHEP::HepVector TiN = Ti * Norm; 
02718   //double si = CLHEP::dot((Ti*Rm0 - Rm0 + di), TiN) / CLHEP::dot(Pm0, TiN); 
02719   //Rm1(1) = CLHEP::dot(ex, Rm0 + si*Pm0 - di);
02720   //Rm1(2) = CLHEP::dot(ey, Rm0 + si*Pm0 - di);
02721   //Rm1(3) = CLHEP::dot(ez, Rm0 + si*Pm0 - di);
02722   double si = CLHEP_dot((Ti*Rm0 - Rm0 + di), TiN) / CLHEP_dot(Pm0, TiN); 
02723   Rm1(1) = CLHEP_dot(ex, Rm0 + si*Pm0 - di);
02724   Rm1(2) = CLHEP_dot(ey, Rm0 + si*Pm0 - di);
02725   Rm1(3) = CLHEP_dot(ez, Rm0 + si*Pm0 - di);
02726   Pm1 = T * Pm0;
02727     
02728   // Global muon parameters after misalignment
02729   Rm = GlobalVector(Rm1(1), Rm1(2), Rm1(3));  
02730   //Pm = GlobalVector(CLHEP::dot(Pm0,ex), CLHEP::dot(Pm0, ey), CLHEP::dot(Pm0,ez));// =T*Pm0 
02731   Pm = GlobalVector(CLHEP_dot(Pm0,ex), CLHEP_dot(Pm0, ey), CLHEP_dot(Pm0,ez));// =T*Pm0 
02732 
02733   // Local muon parameters after misalignment
02734   const Surface& refSurface = tsosMuon.surface();
02735   CLHEP::HepMatrix Tl(3,3,0);
02736 
02737   Tl(1,1) = refSurface.rotation().xx();
02738   Tl(1,2) = refSurface.rotation().xy();
02739   Tl(1,3) = refSurface.rotation().xz();
02740   
02741   Tl(2,1) = refSurface.rotation().yx();
02742   Tl(2,2) = refSurface.rotation().yy();
02743   Tl(2,3) = refSurface.rotation().yz();
02744   
02745   Tl(3,1) = refSurface.rotation().zx();
02746   Tl(3,2) = refSurface.rotation().zy();
02747   Tl(3,3) = refSurface.rotation().zz();
02748 
02749   CLHEP::HepVector R0(3,0), newR0(3,0), newRl(3,0), newPl(3,0);
02750   R0(1) = refSurface.position().x();
02751   R0(2) = refSurface.position().y();
02752   R0(3) = refSurface.position().z();
02753 
02754   newRl = Tl * (Rm1 - R0);
02755   newPl = Tl * Pm1;
02756 
02757   Vm(0) = newPl(1)/newPl(3);
02758   Vm(1) = newPl(2)/newPl(3);
02759   Vm(2) = newRl(1);
02760   Vm(3) = newRl(2);
02761   
02762 #ifdef CHECK_DERIVATIVES_LOCAL_VS_ANGLES
02763   float del = 0.0001;
02764   //int ii = 1; T(3,2) +=del; T(2,3) -=del;
02765   int ii = 2; T(3,1) -=del; T(1,3) +=del;
02766   //int ii = 3; T(1,2) -=del; T(2,1) +=del;
02767   AlgebraicVector4 Vm0 = Vm;
02768   CLHEP::HepVector Rm10 = Rm1, Pm10 = Pm1;;
02769   CLHEP::HepMatrix newTl = Tl*T;
02770   Ti = T.T();
02771   di = -Ti * d;  
02772   ex = Ti*ex0; ey = Ti*ey0; ez = Ti*ez0;
02773   TiN = Ti * Norm; 
02774   si = CLHEP_dot((Ti*Rm0 - Rm0 + di), TiN) / CLHEP_dot(Pm0, TiN); 
02775   Rm1(1) = CLHEP_dot(ex, Rm0 + si*Pm0 - di);
02776   Rm1(2) = CLHEP_dot(ey, Rm0 + si*Pm0 - di);
02777   Rm1(3) = CLHEP_dot(ez, Rm0 + si*Pm0 - di);
02778   Pm1 = T * Pm0;
02779 
02780   newRl = Tl * (Rm1 - R0);
02781   newPl = Tl * Pm1;
02782 
02783   Vm(0) = newPl(1)/newPl(3);
02784   Vm(1) = newPl(2)/newPl(3);
02785   Vm(2) = newRl(1);
02786   Vm(3) = newRl(2);
02787   std::cout<<" ========= dVm/da"<<ii<<" "<<(Vm-Vm0)*(1./del)<<std::endl;
02788 #endif 
02789 
02790   if(debug_){
02791   //std::cout<<" dRm/da3 "<<((Rm1-Rm10)*(1./del)).T()<<" "<<((Pm1-Pm10)*(1./del)).T()<<std::endl; 
02792   //std::cout<<"             Norm "<<Norm.T()<<std::endl;
02793   std::cout<<"             ex  "<<(Tl.T()*ex0).T()<<std::endl;
02794   std::cout<<"             ey  "<<(Tl.T()*ey0).T()<<std::endl;
02795   std::cout<<"             ez  "<<(Tl.T()*ez0).T()<<std::endl;
02796   //std::cpot<<" 
02797 
02798   std::cout<<"    pxyz/p gl 0 "<<(Pm0/Pm0.norm()).T()<<std::endl;
02799   std::cout<<"    pxyz/p loc0 "<<(Tl*Pm0/(Tl*Pm0).norm()).T()<<std::endl;
02800   std::cout<<"    pxyz/p glob  "<<(Pm1/Pm1.norm()).T()<<std::endl;
02801   std::cout<<"    pxyz/p  loc  "<<(newPl/newPl.norm()).T()<<std::endl;
02802 
02803   //std::cout<<"             Rot   "<<refSurface.rotation()<<endl;
02804   //std::cout<<"             Tl   "<<Tl<<endl;
02805   std::cout<<" ---- phi g0 g1 l0 l1  "
02806            <<atan2(Pm0(2),Pm0(1))<<" "<<atan2(Pm1(2),Pm1(1))<<" "
02807            <<atan2((Tl*Pm0)(2),(Tl*Pm0)(1))<<" "
02808            <<atan2(newPl(2),newPl(1))<<std::endl;
02809   std::cout<<" ---- angl Gl Loc "<<atan2(Pm1(2),Pm1(1))-atan2(Pm0(2),Pm0(1))<<" "
02810            <<atan2(newPl(2),newPl(1))-atan2((Tl*Pm0)(2),(Tl*Pm0)(1))<<" "
02811            <<atan2(newPl(3),newPl(2))-atan2((Tl*Pm0)(3),(Tl*Pm0)(2))<<" "
02812            <<atan2(newPl(1),newPl(3))-atan2((Tl*Pm0)(1),(Tl*Pm0)(3))<<" "<<std::endl;
02813   }
02814 
02815   if(debug_){
02816     CLHEP::HepMatrix newTl(3,3,0);
02817     CLHEP::HepVector R0(3,0), newRl(3,0), newPl(3,0);
02818     newTl = Tl * Ti.T();
02819     newR0 = Ti * R0 + di;
02820 
02821     std::cout<<" N "<<Norm.T()<<std::endl;
02822     std::cout<<" dtxl yl "<<Vm(0)-tsosMuon.localParameters().vector()(1)<<" "
02823              <<Vm(1)-tsosMuon.localParameters().vector()(2)<<std::endl;
02824     std::cout<<" dXl dYl "<<Vm(2)-tsosMuon.localParameters().vector()(3)<<" "
02825              <<Vm(3)-tsosMuon.localParameters().vector()(4)<<std::endl;
02826     std::cout<<" localM    "<<tsosMuon.localParameters().vector()<<std::endl;  
02827     std::cout<<"       Vm  "<<Vm<<std::endl;
02828     
02829     
02830     CLHEP::HepVector Rvc(3,0), Pvc(3,0), Rvg(3,0), Pvg(3,0);
02831     Rvc(1) = Vm(2); Rvc(2) = Vm(3);
02832     Pvc(3) = tsosMuon.localParameters().pzSign() * Pm0.norm() /
02833       sqrt(1 + Vm(0)*Vm(0) + Vm(1)*Vm(1));
02834     Pvc(1) = Pvc(3) * Vm(0);
02835     Pvc(2) = Pvc(3) * Vm(1);
02836     
02837     Rvg = newTl.T() * Rvc + newR0;
02838     Pvg = newTl.T() * Pvc;
02839     
02840     std::cout<<" newPl     "<<newPl.T()<<std::endl;
02841     std::cout<<"    Pvc    "<<Pvc.T()<<std::endl;
02842     std::cout<<"    Rvg    "<<Rvg.T()<<std::endl;
02843     std::cout<<"    Rm1    "<<Rm1.T()<<std::endl;
02844     std::cout<<"    Pvg    "<<Pvg.T()<<std::endl;
02845     std::cout<<"    Pm1    "<<Pm1.T()<<std::endl;
02846   }
02847   
02848   if(debug_){    //    ----------  how to transform cartesian from shifted to correct 
02849     
02850     std::cout<<" ----- Pm "<<Pm<<std::endl;
02851     std::cout<<"    T*Pm0 "<<Pm1.T()<<std::endl;
02852     //std::cout<<" Ti*T*Pm0 "<<(Ti*(T*Pm0)).T()<<std::endl;
02853     
02854     //CLHEP::HepVector Rml = (Rm0 + si*Pm0 - di) + di;
02855     CLHEP::HepVector Rml = Rm1(1)*ex + Rm1(2)*ey + Rm1(3)*ez + di;
02856     
02857     CLHEP::HepVector Rt0(3); 
02858     Rt0(1)=Rt.x();  Rt0(2)=Rt.y();  Rt0(3)=Rt.z();      
02859     
02860     //double sl = CLHEP::dot(T*(Rt0 - Rml), T*Norm) / CLHEP::dot(Ti * Pm1, Norm);     
02861     double sl = CLHEP_dot(T*(Rt0 - Rml), T*Norm) / CLHEP_dot(Ti * Pm1, Norm);     
02862     CLHEP::HepVector Rl = Rml + sl*(Ti*Pm1);
02863     
02864     //double A = CLHEP::dot(Ti*Pm1, Norm);
02865     //double B = CLHEP::dot(Rt0, Norm) - CLHEP::dot(Ti*Rm1, Norm) 
02866     //+ CLHEP::dot(Ti*d, Norm); 
02867     double A = CLHEP_dot(Ti*Pm1, Norm);
02868     double B = CLHEP_dot(Rt0, Norm) - CLHEP_dot(Ti*Rm1, Norm) 
02869       + CLHEP_dot(Ti*d, Norm); 
02870     double s = B/A;
02871     CLHEP::HepVector Dr = Ti*Rm1 - Ti*d  + s*(Ti*Pm1) - Rm0;
02872     CLHEP::HepVector Dp = Ti*Pm1 - Pm0;
02873 
02874     CLHEP::HepVector TiR = Ti * Rm0 + di;
02875     CLHEP::HepVector Ri = T * TiR + d;
02876     
02877     std::cout<<" --TTi-0 "<<(Ri-Rm0).T()<<std::endl;
02878     std::cout<<" --  Dr  "<<Dr.T()<<endl;
02879     std::cout<<" --  Dp  "<<Dp.T()<<endl;
02880     //std::cout<<" ex "<<ex.T()<<endl;
02881     //std::cout<<" ey "<<ey.T()<<endl;
02882     //std::cout<<" ez "<<ez.T()<<endl;
02883     //std::cout<<" ---- T  ---- "<<T<<std::endl;
02884     //std::cout<<" ---- Ti ---- "<<Ti<<std::endl;
02885     //std::cout<<" ---- d  ---- "<<d.T()<<std::endl;
02886     //std::cout<<" ---- di ---- "<<di.T()<<std::endl;
02887     std::cout<<" -- si sl s "<<si<<" "<<sl<<" "<<s<<std::endl;
02888     //std::cout<<" -- si sl "<<si<<" "<<sl<<endl;
02889     //std::cout<<" -- si s "<<si<<" "<<s<<endl;
02890     std::cout<<" -- Rml-(Rm0+si*Pm0) "<<(Rml - Rm0 - si*Pm0).T()<<std::endl;
02891     //std::cout<<" -- Rm0 "<<Rm0.T()<<std::endl;
02892     //std::cout<<" -- Rm1 "<<Rm1.T()<<std::endl;
02893     //std::cout<<" -- Rmi "<<Rmi.T()<<std::endl;
02894     //std::cout<<" --siPm "<<(si*Pm0).T()<<std::endl;
02895     //std::cout<<" --s2Pm "<<(s2*(T * Pm1)).T()<<std::endl;
02896     //std::cout<<" --R1-0 "<<(Rm1-Rm0).T()<<std::endl;
02897     //std::cout<<" --Ri-0 "<<(Rmi-Rm0).T()<<std::endl;
02898     std::cout<<" --Rl-0 "<<(Rl-Rm0).T()<<std::endl;
02899     //std::cout<<" --Pi-0 "<<(Pmi-Pm0).T()<<std::endl;
02900   }      //                                --------   end of debug
02901   
02902   return;
02903 
02904 }  // end misalignMuonL
02905 
02906 
02907 // ----  refit any direction of transient track ------
02908 void 
02909 GlobalTrackerMuonAlignment::trackFitter(reco::TrackRef alongTr, reco::TransientTrack& alongTTr,
02910                                        PropagationDirection direction,
02911                                        TrajectoryStateOnSurface& trackFittedTSOS)
02912 {
02913   bool info = false;
02914   bool smooth = false;
02915   
02916   if(info)
02917     std::cout<<" ......... start of trackFitter ......... "<<std::endl;  
02918   if(false && info){ 
02919     this->debugTrackHit(" Tracker track hits ", alongTr); 
02920     this->debugTrackHit(" Tracker TransientTrack hits ", alongTTr); 
02921   }
02922 
02923   edm::OwnVector<TrackingRecHit> recHit;
02924   DetId trackDetId = DetId(alongTr->innerDetId());
02925   for (trackingRecHit_iterator i=alongTr->recHitsBegin();i!=alongTr->recHitsEnd(); i++) {
02926     recHit.push_back((*i)->clone());
02927   }
02928   if(direction != alongMomentum)
02929     {
02930       edm::OwnVector<TrackingRecHit> recHitAlong = recHit;
02931       recHit.clear();
02932       for(unsigned int ihit = recHitAlong.size()-1; ihit+1>0; ihit--){
02933         recHit.push_back(recHitAlong[ihit]);
02934       }
02935       recHitAlong.clear();
02936     }
02937   if(info)
02938     std::cout<<" ... Own recHit.size() "<<recHit.size()<<" "<<trackDetId.rawId()<<std::endl;
02939 
02940   //MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitTrack;
02941   TransientTrackingRecHit::RecHitContainer recHitTrack;
02942   for (trackingRecHit_iterator i=alongTr->recHitsBegin();i!=alongTr->recHitsEnd(); i++) {
02943     if((*i)->isValid()){
02944       recHitTrack.push_back(TTRHBuilder->build(&**i ));} 
02945     else{
02946       recHitTrack.push_back(TTRHBuilder->build(&**i ));}  
02947   }
02948   
02949   if(info)
02950     std::cout<<" ... recHitTrack.size() "<<recHit.size()<<" "<<trackDetId.rawId()
02951              <<" ======> recHitMu "<<std::endl;
02952 
02953   MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu = recHitTrack;
02954   /*
02955     MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu =
02956     MuRHBuilder->build(alongTTr.recHitsBegin(), alongTTr.recHitsEnd());
02957   */
02958   if(info)
02959     std::cout<<" ...along.... recHitMu.size() "<<recHitMu.size()<<std::endl;
02960   if(direction != alongMomentum){
02961     MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMuAlong = recHitMu;
02962     recHitMu.clear();
02963     for(unsigned int ihit = recHitMuAlong.size()-1; ihit+1>0; ihit--){
02964       recHitMu.push_back(recHitMuAlong[ihit]);
02965     }
02966     recHitMuAlong.clear();
02967   }
02968   if(info)
02969     std::cout<<" ...opposite... recHitMu.size() "<<recHitMu.size()<<std::endl;
02970 
02971   TrajectoryStateOnSurface firstTSOS;
02972   if(direction == alongMomentum) firstTSOS = alongTTr.innermostMeasurementState();
02973   else                           firstTSOS = alongTTr.outermostMeasurementState();
02974 
02975   AlgebraicSymMatrix55 CovLoc;
02976   CovLoc(0,0) = 1.   * firstTSOS.localError().matrix()(0,0);
02977   CovLoc(1,1) = 10.  * firstTSOS.localError().matrix()(1,1);
02978   CovLoc(2,2) = 10.  * firstTSOS.localError().matrix()(2,2);
02979   CovLoc(3,3) = 100. * firstTSOS.localError().matrix()(3,3);
02980   CovLoc(4,4) = 100. * firstTSOS.localError().matrix()(4,4);
02981   TrajectoryStateOnSurface initialTSOS(firstTSOS.localParameters(), LocalTrajectoryError(CovLoc),
02982                                        firstTSOS.surface(), &*magneticField_);
02983   
02984   PTrajectoryStateOnDet  PTraj = 
02985     trajectoryStateTransform::persistentState(initialTSOS, trackDetId.rawId());
02986   //const TrajectorySeed seedT(*PTraj, recHit, alongMomentum);  
02987   const TrajectorySeed seedT(PTraj, recHit, direction);  
02988 
02989   std::vector<Trajectory> trajVec;
02990   if(direction == alongMomentum) trajVec = theFitter->fit(seedT, recHitMu, initialTSOS);
02991   else                           trajVec = theFitterOp->fit(seedT, recHitMu, initialTSOS);
02992   
02993   if(info){
02994     this->debugTrajectorySOSv(" innermostTSOS of TransientTrack", 
02995                              alongTTr.innermostMeasurementState());
02996     this->debugTrajectorySOSv(" outermostTSOS of TransientTrack", 
02997                              alongTTr.outermostMeasurementState());
02998     this->debugTrajectorySOS(" initialTSOS for theFitter ", initialTSOS);
02999     std::cout<<" . . . . . trajVec.size() "<<trajVec.size()<<std::endl;
03000     if(trajVec.size()) this->debugTrajectory(" theFitter trajectory ", trajVec[0]);
03001   }
03002 
03003   if(!smooth){
03004     if(trajVec.size()) trackFittedTSOS = trajVec[0].lastMeasurement().updatedState();}
03005   else{
03006     std::vector<Trajectory> trajSVec;
03007     if (trajVec.size()){
03008       if(direction == alongMomentum) trajSVec = theSmoother->trajectories(*(trajVec.begin()));
03009       else                           trajSVec = theSmootherOp->trajectories(*(trajVec.begin()));
03010     }
03011     if(info) std::cout<<" . . . . trajSVec.size() "<<trajSVec.size()<<std::endl;
03012     if(trajSVec.size()) this->debugTrajectorySOSv("smoothed geom InnermostState",
03013                                                   trajSVec[0].geometricalInnermostState());
03014     if(trajSVec.size()) trackFittedTSOS = trajSVec[0].geometricalInnermostState();
03015   }
03016 
03017   if(info) this->debugTrajectorySOSv(" trackFittedTSOS ", trackFittedTSOS);
03018 
03019 } // end of trackFitter
03020 
03021 // ----  refit any direction of muon transient track ------
03022 void 
03023 GlobalTrackerMuonAlignment::muonFitter(reco::TrackRef alongTr, reco::TransientTrack& alongTTr,
03024                                        PropagationDirection direction,
03025                                        TrajectoryStateOnSurface& trackFittedTSOS)
03026 {  
03027   bool info = false;
03028   bool smooth = false;
03029   
03030   if(info)  std::cout<<" ......... start of muonFitter ........"<<std::endl;  
03031   if(false && info){  
03032     this->debugTrackHit(" Muon track hits ", alongTr); 
03033     this->debugTrackHit(" Muon TransientTrack hits ", alongTTr); 
03034   }
03035 
03036   edm::OwnVector<TrackingRecHit> recHit;
03037   DetId trackDetId = DetId(alongTr->innerDetId());
03038   for (trackingRecHit_iterator i=alongTr->recHitsBegin();i!=alongTr->recHitsEnd(); i++) {
03039     recHit.push_back((*i)->clone());
03040   }
03041   if(direction != alongMomentum)
03042     {
03043       edm::OwnVector<TrackingRecHit> recHitAlong = recHit;
03044       recHit.clear();
03045       for(unsigned int ihit = recHitAlong.size()-1; ihit+1>0; ihit--){
03046         recHit.push_back(recHitAlong[ihit]);
03047       }
03048       recHitAlong.clear();
03049     }
03050   if(info)
03051     std::cout<<" ... Own recHit.size() DetId==Muon "<<recHit.size()<<" "<<trackDetId.rawId()<<std::endl;
03052 
03053   MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu =
03054     MuRHBuilder->build(alongTTr.recHitsBegin(), alongTTr.recHitsEnd());
03055   if(info)
03056     std::cout<<" ...along.... recHitMu.size() "<<recHitMu.size()<<std::endl;
03057   if(direction != alongMomentum){
03058     MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMuAlong = recHitMu;
03059     recHitMu.clear();
03060     for(unsigned int ihit = recHitMuAlong.size()-1; ihit+1>0; ihit--){
03061       recHitMu.push_back(recHitMuAlong[ihit]);
03062     }
03063     recHitMuAlong.clear();
03064   }
03065   if(info)
03066     std::cout<<" ...opposite... recHitMu.size() "<<recHitMu.size()<<std::endl;
03067 
03068   TrajectoryStateOnSurface firstTSOS;
03069   if(direction == alongMomentum) firstTSOS = alongTTr.innermostMeasurementState();
03070   else                           firstTSOS = alongTTr.outermostMeasurementState();
03071 
03072   AlgebraicSymMatrix55 CovLoc;
03073   CovLoc(0,0) = 1.   * firstTSOS.localError().matrix()(0,0);
03074   CovLoc(1,1) = 10.  * firstTSOS.localError().matrix()(1,1);
03075   CovLoc(2,2) = 10.  * firstTSOS.localError().matrix()(2,2);
03076   CovLoc(3,3) = 100. * firstTSOS.localError().matrix()(3,3);
03077   CovLoc(4,4) = 100. * firstTSOS.localError().matrix()(4,4);
03078   TrajectoryStateOnSurface initialTSOS(firstTSOS.localParameters(), LocalTrajectoryError(CovLoc),
03079                                        firstTSOS.surface(), &*magneticField_);
03080   
03081   PTrajectoryStateOnDet PTraj = 
03082     trajectoryStateTransform::persistentState(initialTSOS, trackDetId.rawId());
03083   //const TrajectorySeed seedT(*PTraj, recHit, alongMomentum);  
03084   const TrajectorySeed seedT(PTraj, recHit, direction);  
03085 
03086   std::vector<Trajectory> trajVec;
03087   if(direction == alongMomentum) trajVec = theFitter->fit(seedT, recHitMu, initialTSOS);
03088   else                           trajVec = theFitterOp->fit(seedT, recHitMu, initialTSOS);
03089   
03090   if(info){
03091     this->debugTrajectorySOSv(" innermostTSOS of TransientTrack", 
03092                              alongTTr.innermostMeasurementState());
03093     this->debugTrajectorySOSv(" outermostTSOS of TransientTrack", 
03094                              alongTTr.outermostMeasurementState());
03095     this->debugTrajectorySOS(" initialTSOS for theFitter ", initialTSOS);
03096     std::cout<<" . . . . . trajVec.size() "<<trajVec.size()<<std::endl;
03097     if(trajVec.size()) this->debugTrajectory(" theFitter trajectory ", trajVec[0]);
03098   }
03099 
03100   if(!smooth){
03101     if(trajVec.size()) trackFittedTSOS = trajVec[0].lastMeasurement().updatedState();}
03102   else{
03103     std::vector<Trajectory> trajSVec;
03104     if (trajVec.size()){
03105       if(direction == alongMomentum) trajSVec = theSmoother->trajectories(*(trajVec.begin()));
03106       else                           trajSVec = theSmootherOp->trajectories(*(trajVec.begin()));
03107     }
03108     if(info) std::cout<<" . . . . trajSVec.size() "<<trajSVec.size()<<std::endl;
03109     if(trajSVec.size()) this->debugTrajectorySOSv("smoothed geom InnermostState",
03110                                                   trajSVec[0].geometricalInnermostState());
03111     if(trajSVec.size()) trackFittedTSOS = trajSVec[0].geometricalInnermostState();
03112   }
03113 
03114 } // end of muonFitter
03115 
03116 
03117 // ----  debug printout of hits from TransientTrack  ------
03118 void GlobalTrackerMuonAlignment::debugTrackHit(const std::string title, TransientTrack& alongTr) 
03119 {
03120   std::cout<<" ------- "<<title<<" --------"<<std::endl;
03121   int nHit = 1;
03122   TransientTrackingRecHit::RecHitContainer recHit;
03123   for (trackingRecHit_iterator i=alongTr.recHitsBegin();i!=alongTr.recHitsEnd(); i++) {
03124     std::cout<<" Hit "<<nHit++<<" DetId "<<(*i)->geographicalId().det();
03125     if( (*i)->geographicalId().det() == DetId::Tracker ) std::cout<<" Tracker ";
03126     else if ( (*i)->geographicalId().det() == DetId::Muon ) std::cout<<" Muon "; 
03127     else std::cout<<" Unknown ";
03128     if(!(*i)->isValid()) std::cout<<" not valid "<<std::endl;  
03129     else std::cout<<std::endl;
03130   }
03131 }
03132 
03133 
03134 // ----  debug printout of hits from TrackRef   ------
03135 void GlobalTrackerMuonAlignment::debugTrackHit(const std::string title, reco::TrackRef alongTr) 
03136 {
03137   std::cout<<" ------- "<<title<<" --------"<<std::endl;
03138   int nHit = 1;
03139   TransientTrackingRecHit::RecHitContainer recHit;
03140   for (trackingRecHit_iterator i=alongTr->recHitsBegin();i!=alongTr->recHitsEnd(); i++) {
03141     std::cout<<" Hit "<<nHit++<<" DetId "<<(*i)->geographicalId().det();
03142     if( (*i)->geographicalId().det() == DetId::Tracker ) std::cout<<" Tracker ";
03143     else if ( (*i)->geographicalId().det() == DetId::Muon ) std::cout<<" Muon "; 
03144     else std::cout<<" Unknown ";
03145     if(!(*i)->isValid()) std::cout<<" not valid "<<std::endl;  
03146     else std::cout<<std::endl;
03147   }
03148 }
03149 
03150 // ----  debug printout TrajectoryStateOnSurface  ------
03151 void GlobalTrackerMuonAlignment::debugTrajectorySOS(const std::string title, 
03152                                                      TrajectoryStateOnSurface& trajSOS) 
03153 {
03154   std::cout<<"    --- "<<title<<" --- "<<std::endl;
03155   if(!trajSOS.isValid()) {std::cout<<"      Not valid !!!! "<<std::endl; return;}
03156   std::cout<<" R |p| "<<trajSOS.globalPosition().perp()<<" "
03157            <<trajSOS.globalMomentum().mag()<<" charge "<<trajSOS.charge()<<std::endl; 
03158   std::cout<<" x p  "<<trajSOS.globalParameters().vector()(0)<<" "
03159            <<trajSOS.globalParameters().vector()(1)<<" "
03160            <<trajSOS.globalParameters().vector()(2)<<" "
03161            <<trajSOS.globalParameters().vector()(3)<<" "
03162            <<trajSOS.globalParameters().vector()(4)<<" "
03163            <<trajSOS.globalParameters().vector()(5)<<std::endl;
03164   std::cout<<" +/- "<<sqrt(trajSOS.cartesianError().matrix()(0,0))<<" "
03165            <<sqrt(trajSOS.cartesianError().matrix()(1,1))<<" "
03166            <<sqrt(trajSOS.cartesianError().matrix()(2,2))<<" "
03167            <<sqrt(trajSOS.cartesianError().matrix()(3,3))<<" "
03168            <<sqrt(trajSOS.cartesianError().matrix()(4,4))<<" "
03169            <<sqrt(trajSOS.cartesianError().matrix()(5,5))<<std::endl;
03170   std::cout<<" q/p dxy/dz xy "<<trajSOS.localParameters().vector()(0)<<" "
03171            <<trajSOS.localParameters().vector()(1)<<" "
03172            <<trajSOS.localParameters().vector()(2)<<" "
03173            <<trajSOS.localParameters().vector()(3)<<" "
03174            <<trajSOS.localParameters().vector()(4)<<std::endl;
03175   std::cout<<"   +/- error  "<<sqrt(trajSOS.localError().matrix()(0,0))<<" "
03176            <<sqrt(trajSOS.localError().matrix()(1,1))<<" "
03177            <<sqrt(trajSOS.localError().matrix()(2,2))<<" "
03178            <<sqrt(trajSOS.localError().matrix()(3,3))<<" "
03179            <<sqrt(trajSOS.localError().matrix()(4,4))<<" "<<std::endl;
03180 }
03181 
03182 // ----  debug printout TrajectoryStateOnSurface  ------
03183 void GlobalTrackerMuonAlignment::debugTrajectorySOSv(const std::string title, 
03184                                                      TrajectoryStateOnSurface trajSOS) 
03185 {
03186   std::cout<<"    --- "<<title<<" --- "<<std::endl;
03187   if(!trajSOS.isValid()) {std::cout<<"      Not valid !!!! "<<std::endl; return;}
03188   std::cout<<" R |p| "<<trajSOS.globalPosition().perp()<<" "
03189            <<trajSOS.globalMomentum().mag()<<" charge "<<trajSOS.charge()<<std::endl; 
03190   std::cout<<" x p  "<<trajSOS.globalParameters().vector()(0)<<" "
03191            <<trajSOS.globalParameters().vector()(1)<<" "
03192            <<trajSOS.globalParameters().vector()(2)<<" "
03193            <<trajSOS.globalParameters().vector()(3)<<" "
03194            <<trajSOS.globalParameters().vector()(4)<<" "
03195            <<trajSOS.globalParameters().vector()(5)<<std::endl;
03196   std::cout<<" +/- "<<sqrt(trajSOS.cartesianError().matrix()(0,0))<<" "
03197            <<sqrt(trajSOS.cartesianError().matrix()(1,1))<<" "
03198            <<sqrt(trajSOS.cartesianError().matrix()(2,2))<<" "
03199            <<sqrt(trajSOS.cartesianError().matrix()(3,3))<<" "
03200            <<sqrt(trajSOS.cartesianError().matrix()(4,4))<<" "
03201            <<sqrt(trajSOS.cartesianError().matrix()(5,5))<<std::endl;
03202   std::cout<<" q/p dxy/dz xy "<<trajSOS.localParameters().vector()(0)<<" "
03203            <<trajSOS.localParameters().vector()(1)<<" "
03204            <<trajSOS.localParameters().vector()(2)<<" "
03205            <<trajSOS.localParameters().vector()(3)<<" "
03206            <<trajSOS.localParameters().vector()(4)<<std::endl;
03207   std::cout<<"   +/- error  "<<sqrt(trajSOS.localError().matrix()(0,0))<<" "
03208            <<sqrt(trajSOS.localError().matrix()(1,1))<<" "
03209            <<sqrt(trajSOS.localError().matrix()(2,2))<<" "
03210            <<sqrt(trajSOS.localError().matrix()(3,3))<<" "
03211            <<sqrt(trajSOS.localError().matrix()(4,4))<<" "<<std::endl;
03212 }
03213 
03214 // ----  debug printout Trajectory   ------
03215 void GlobalTrackerMuonAlignment::debugTrajectory(const std::string title, Trajectory& traj) 
03216 {
03217   std::cout<<"\n"<<"    ...... "<<title<<" ...... "<<std::endl;
03218   if(!traj.isValid()) {std::cout<<"          Not valid !!!!!!!!  "<<std::endl; return;}
03219   std::cout<<" chi2/Nhit "<<traj.chiSquared()<<" / "<<traj.foundHits();
03220   if(traj.direction() == alongMomentum) std::cout<<" alongMomentum  >>>>"<<std::endl;
03221   else                                  std::cout<<" oppositeToMomentum  <<<<"<<std::endl;
03222   this->debugTrajectorySOSv(" firstMeasurementTSOS ",traj.firstMeasurement().updatedState());
03223   //this->debugTrajectorySOSv(" firstMeasurementTSOS ",traj.firstMeasurement().predictedState());
03224   this->debugTrajectorySOSv(" lastMeasurementTSOS ",traj.lastMeasurement().updatedState());
03225   //this->debugTrajectorySOSv(" geom InnermostState", traj.geometricalInnermostState());
03226   std::cout<<"    . . . . . . . . . . . . . . . . . . . . . . . . . . . . \n"<<std::endl;
03227 }
03228 
03229 
03230 // ----  write GlobalPositionRcd   ------
03231 void GlobalTrackerMuonAlignment::writeGlPosRcd(CLHEP::HepVector& paramVec) 
03232 {
03233   //paramVec(1) = 0.1; paramVec(2) = 0.2; paramVec(3) = 0.3;         //0123
03234   //paramVec(4) = 0.0001; paramVec(5) = 0.0002; paramVec(6) = 0.0003;
03235   //paramVec(1) = 0.3; paramVec(2) = 0.4; paramVec(3) = 0.5;         //03450123
03236   //paramVec(4) = 0.0001; paramVec(5) = 0.0002; paramVec(6) = 0.0003;
03237   //paramVec(1) = 2.; paramVec(2) = 2.; paramVec(3) = 2.;               //222
03238   //paramVec(4) = 0.02; paramVec(5) = 0.02; paramVec(6) = 0.02;
03239   //paramVec(1) = -10.; paramVec(2) = -20.; paramVec(3) = -30.;         //102030
03240   //paramVec(4) = -0.1; paramVec(5) = -0.2; paramVec(6) = -0.3;
03241   //paramVec(1) = 0.; paramVec(2) = 0.; paramVec(3) = 1.;               //1z
03242   //paramVec(4) = 0.0000; paramVec(5) = 0.0000; paramVec(6) = 0.0000;
03243 
03244   std::cout<<" paramVector "<<paramVec.T()<<std::endl;
03245   
03246   CLHEP::Hep3Vector colX, colY, colZ;
03247 
03248 #ifdef NOT_EXACT_ROTATION_MATRIX     
03249   colX = CLHEP::Hep3Vector(          1.,  -paramVec(6),  paramVec(5));
03250   colY = CLHEP::Hep3Vector( paramVec(6),            1., -paramVec(4));
03251   colZ = CLHEP::Hep3Vector(-paramVec(5),   paramVec(4),           1.);
03252 #else 
03253   double s1 = std::sin(paramVec(4)), c1 = std::cos(paramVec(4));
03254   double s2 = std::sin(paramVec(5)), c2 = std::cos(paramVec(5));
03255   double s3 = std::sin(paramVec(6)), c3 = std::cos(paramVec(6));
03256   colX = CLHEP::Hep3Vector(               c2 * c3,               -c2 * s3,       s2);
03257   colY = CLHEP::Hep3Vector(c1 * s3 + s1 * s2 * c3, c1 * c3 - s1 * s2 * s3, -s1 * c2);
03258   colZ = CLHEP::Hep3Vector(s1 * s3 - c1 * s2 * c3, s1 * c3 + c1 * s2 * s3,  c1 * c2);
03259 #endif
03260 
03261   CLHEP::HepVector  param0(6,0);
03262 
03263   Alignments* globalPositions = new Alignments();  
03264 
03265   //Tracker
03266   AlignTransform tracker(iteratorTrackerRcd->translation(),
03267                          iteratorTrackerRcd->rotation(),
03268                          DetId(DetId::Tracker).rawId());
03269   // Muon
03270   CLHEP::Hep3Vector posMuGlRcd = iteratorMuonRcd->translation();
03271   CLHEP::HepRotation rotMuGlRcd = iteratorMuonRcd->rotation();
03272   CLHEP::HepEulerAngles angMuGlRcd = iteratorMuonRcd->rotation().eulerAngles();
03273   if(debug_)
03274     std::cout<<" Old muon Rcd Euler angles "<<angMuGlRcd.phi()<<" "<<angMuGlRcd.theta()
03275              <<" "<<angMuGlRcd.psi()<<std::endl;
03276   AlignTransform muon;
03277   if((angMuGlRcd.phi() == 0.) && (angMuGlRcd.theta() == 0.) && (angMuGlRcd.psi() == 0.) &&
03278      (posMuGlRcd.x() == 0.) && (posMuGlRcd.y() == 0.) && (posMuGlRcd.z() == 0.)){
03279     std::cout<<" New muon parameters are stored in Rcd"<<std::endl; 
03280 
03281     AlignTransform muonNew(AlignTransform::Translation(paramVec(1),
03282                                                        paramVec(2),
03283                                                        paramVec(3)),
03284                            AlignTransform::Rotation   (colX,
03285                                                        colY,
03286                                                        colZ),
03287                            DetId(DetId::Muon).rawId());
03288     muon = muonNew;
03289   }
03290   else if((paramVec(1) == 0.) && (paramVec(2) == 0.) && (paramVec(3) == 0.) && 
03291           (paramVec(4) == 0.) && (paramVec(5) == 0.) && (paramVec(6) == 0.)){
03292     std::cout<<" Old muon parameters are stored in Rcd"<<std::endl; 
03293     
03294     AlignTransform muonNew(iteratorMuonRcd->translation(),
03295                            iteratorMuonRcd->rotation(),
03296                            DetId(DetId::Muon).rawId());
03297     muon = muonNew;  
03298   }
03299   else{
03300     std::cout<<" New + Old muon parameters are stored in Rcd"<<std::endl; 
03301     CLHEP::Hep3Vector posMuGlRcdThis = CLHEP::Hep3Vector(paramVec(1), paramVec(2), paramVec(3));   
03302     CLHEP::HepRotation rotMuGlRcdThis = CLHEP::HepRotation(colX, colY, colZ);   
03303     CLHEP::Hep3Vector posMuGlRcdNew = rotMuGlRcdThis * posMuGlRcd + posMuGlRcdThis;
03304     CLHEP::HepRotation rotMuGlRcdNew = rotMuGlRcdThis * rotMuGlRcd;
03305 
03306     AlignTransform muonNew(posMuGlRcdNew,
03307                            rotMuGlRcdNew, 
03308                            DetId(DetId::Muon).rawId());
03309     muon = muonNew;
03310   }
03311 
03312   // Ecal
03313   AlignTransform ecal(iteratorEcalRcd->translation(),
03314                          iteratorEcalRcd->rotation(),
03315                          DetId(DetId::Ecal).rawId());
03316   // Hcal
03317   AlignTransform hcal(iteratorHcalRcd->translation(),
03318                          iteratorHcalRcd->rotation(),
03319                          DetId(DetId::Hcal).rawId());
03320   // Calo
03321   AlignTransform calo(AlignTransform::Translation(param0(1),
03322                                                   param0(2),
03323                                                   param0(3)),
03324                       AlignTransform::EulerAngles(param0(4),
03325                                                   param0(5),
03326                                                   param0(6)),
03327                       DetId(DetId::Calo).rawId());
03328 
03329   std::cout << "Tracker (" << tracker.rawId() << ") at " << tracker.translation() 
03330             << " " << tracker.rotation().eulerAngles() << std::endl;
03331   std::cout << tracker.rotation() << std::endl;
03332   
03333   std::cout << "Muon (" << muon.rawId() << ") at " << muon.translation() 
03334             << " " << muon.rotation().eulerAngles() << std::endl;
03335   std::cout << "          rotations angles around x,y,z "
03336             << " ( " << -muon.rotation().zy() << " " << muon.rotation().zx() 
03337             << " " << -muon.rotation().yx() << " )" << std::endl;
03338    std::cout << muon.rotation() << std::endl;
03339    
03340    std::cout << "Ecal (" << ecal.rawId() << ") at " << ecal.translation() 
03341              << " " << ecal.rotation().eulerAngles() << std::endl;
03342    std::cout << ecal.rotation() << std::endl;
03343    
03344    std::cout << "Hcal (" << hcal.rawId() << ") at " << hcal.translation()
03345              << " " << hcal.rotation().eulerAngles() << std::endl;
03346    std::cout << hcal.rotation() << std::endl;
03347    
03348    std::cout << "Calo (" << calo.rawId() << ") at " << calo.translation()
03349              << " " << calo.rotation().eulerAngles() << std::endl;
03350    std::cout << calo.rotation() << std::endl;
03351 
03352    globalPositions->m_align.push_back(tracker);
03353    globalPositions->m_align.push_back(muon);
03354    globalPositions->m_align.push_back(ecal);
03355    globalPositions->m_align.push_back(hcal);
03356    globalPositions->m_align.push_back(calo);
03357 
03358    std::cout << "Uploading to the database..." << std::endl;
03359 
03360    edm::Service<cond::service::PoolDBOutputService> poolDbService;
03361    
03362    if (!poolDbService.isAvailable())
03363      throw cms::Exception("NotAvailable") << "PoolDBOutputService not available";
03364                   
03365    //    if (poolDbService->isNewTagRequest("GlobalPositionRcd")) {
03366 //       poolDbService->createNewIOV<Alignments>(&(*globalPositions), poolDbService->endOfTime(), "GlobalPositionRcd");
03367 //    } else {
03368 //       poolDbService->appendSinceTime<Alignments>(&(*globalPositions), poolDbService->currentTime(), "GlobalPositionRcd");
03369 //    }
03370    poolDbService->writeOne<Alignments>(&(*globalPositions), 
03371                                        poolDbService->currentTime(),
03372                                        //poolDbService->beginOfTime(),
03373                                        "GlobalPositionRcd");
03374    std::cout << "done!" << std::endl;
03375    
03376    return;
03377 }
03378 
03379 
03380 //define this as a plug-in
03381 DEFINE_FWK_MODULE(GlobalTrackerMuonAlignment);