CMS 3D CMS Logo
List of all members | Public Member Functions | Private Member Functions | Private Attributes
GlobalTrackerMuonAlignment Class Reference

#include <Alignment/GlobalTrackerMuonAlignment/src/GlobalTrackerMuonAlignment.cc>

Inheritance diagram for GlobalTrackerMuonAlignment:
edm::EDAnalyzer edm::EDConsumerBase

Public Member Functions

void analyzeTrackTrack (const edm::Event &, const edm::EventSetup &)
 
void analyzeTrackTrajectory (const edm::Event &, const edm::EventSetup &)
 
void bookHist ()
 
double CLHEP_dot (const CLHEP::HepVector &a, const CLHEP::HepVector &b)
 
void debugTrackHit (const std::string, reco::TrackRef)
 
void debugTrackHit (const std::string, reco::TransientTrack &)
 
void debugTrajectory (const std::string, Trajectory &)
 
void debugTrajectorySOS (const std::string, TrajectoryStateOnSurface &)
 
void debugTrajectorySOSv (const std::string, TrajectoryStateOnSurface)
 
void fitHist ()
 
 GlobalTrackerMuonAlignment (const edm::ParameterSet &)
 
void gradientGlobal (GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, AlgebraicSymMatrix66 &)
 
void gradientGlobalAlg (GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, AlgebraicSymMatrix66 &)
 
void gradientLocal (GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, CLHEP::HepSymMatrix &, CLHEP::HepMatrix &, CLHEP::HepVector &, AlgebraicVector4 &)
 
void misalignMuon (GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &)
 
void misalignMuonL (GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, AlgebraicVector4 &, TrajectoryStateOnSurface &, TrajectoryStateOnSurface &)
 
void muonFitter (reco::TrackRef, reco::TransientTrack &, PropagationDirection, TrajectoryStateOnSurface &)
 
void trackFitter (reco::TrackRef, reco::TransientTrack &, PropagationDirection, TrajectoryStateOnSurface &)
 
void writeGlPosRcd (CLHEP::HepVector &d3)
 
 ~GlobalTrackerMuonAlignment () override
 
- Public Member Functions inherited from edm::EDAnalyzer
void callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)
 
 EDAnalyzer ()
 
SerialTaskQueue * globalLuminosityBlocksQueue ()
 
SerialTaskQueue * globalRunsQueue ()
 
ModuleDescription const & moduleDescription () const
 
std::string workerType () const
 
 ~EDAnalyzer () override
 
- Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfo > consumesInfo () const
 
void convertCurrentProcessAlias (std::string const &processName)
 Convert "@currentProcess" in InputTag process names to the actual current process name. More...
 
 EDConsumerBase ()
 
 EDConsumerBase (EDConsumerBase &&)=default
 
 EDConsumerBase (EDConsumerBase const &)=delete
 
ESProxyIndex const * esGetTokenIndices (edm::Transition iTrans) const
 
std::vector< ESProxyIndex > const & esGetTokenIndicesVector (edm::Transition iTrans) const
 
std::vector< ESRecordIndex > const & esGetTokenRecordIndicesVector (edm::Transition iTrans) const
 
ProductResolverIndexAndSkipBit indexFrom (EDGetToken, BranchType, TypeID const &) const
 
void itemsMayGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const
 
void itemsToGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const
 
std::vector< ProductResolverIndexAndSkipBit > const & itemsToGetFrom (BranchType iType) const
 
void labelsForToken (EDGetToken iToken, Labels &oLabels) const
 
void modulesWhoseProductsAreConsumed (std::vector< ModuleDescription const * > &modules, ProductRegistry const &preg, std::map< std::string, ModuleDescription const * > const &labelsToDesc, std::string const &processName) const
 
EDConsumerBase & operator= (EDConsumerBase &&)=default
 
EDConsumerBase const & operator= (EDConsumerBase const &)=delete
 
bool registeredToConsume (ProductResolverIndex, bool, BranchType) const
 
bool registeredToConsumeMany (TypeID const &, BranchType) const
 
ProductResolverIndexAndSkipBit uncheckedIndexFrom (EDGetToken) const
 
void updateLookup (BranchType iBranchType, ProductResolverIndexHelper const &, bool iPrefetchMayGet)
 
void updateLookup (eventsetup::ESRecordsToProxyIndices const &)
 
virtual ~EDConsumerBase () noexcept(false)
 

Private Member Functions

void analyze (const edm::Event &, const edm::EventSetup &) override
 
void beginJob () override
 
void endJob () override
 

Private Attributes

bool collectionCosmic
 
bool collectionIsolated
 
bool cosmicMuonMode_
 
bool debug_
 
bool defineFitter
 
TFile * file
 
AlgebraicVector3 Gfr
 
const Alignments * globalPositionRcd_
 
edm::InputTag gmuonTags_
 
CLHEP::HepVector Grad
 
CLHEP::HepVector GradL
 
CLHEP::HepMatrix Hess
 
CLHEP::HepMatrix HessL
 
TH1F * histo
 
TH2F * histo101
 
TH2F * histo102
 
TH1F * histo11
 
TH1F * histo12
 
TH1F * histo13
 
TH1F * histo14
 
TH1F * histo15
 
TH1F * histo16
 
TH1F * histo17
 
TH1F * histo18
 
TH1F * histo19
 
TH1F * histo2
 
TH1F * histo20
 
TH1F * histo21
 
TH1F * histo22
 
TH1F * histo23
 
TH1F * histo24
 
TH1F * histo25
 
TH1F * histo26
 
TH1F * histo27
 
TH1F * histo28
 
TH1F * histo29
 
TH1F * histo3
 
TH1F * histo30
 
TH1F * histo31
 
TH1F * histo32
 
TH1F * histo33
 
TH1F * histo34
 
TH1F * histo35
 
TH1F * histo4
 
TH1F * histo5
 
TH1F * histo6
 
TH1F * histo7
 
TH1F * histo8
 
AlgebraicSymMatrix33 Inf
 
bool isolatedMuonMode_
 
std::vector< AlignTransform >::const_iterator iteratorEcalRcd
 
std::vector< AlignTransform >::const_iterator iteratorHcalRcd
 
std::vector< AlignTransform >::const_iterator iteratorMuonRcd
 
std::vector< AlignTransform >::const_iterator iteratorTrackerRcd
 
const MagneticField * magneticField_
 
CLHEP::HepVector MuGlAngle
 
CLHEP::HepVector MuGlShift
 
edm::InputTag muonTags_
 
MuonTransientTrackingRecHitBuilder * MuRHBuilder
 
std::string MuSelect
 
int N_event
 
int N_track
 
ofstream OutGlobalTxt
 
string propagator_
 
bool refitMuon_
 
bool refitTrack_
 
string rootOutFile_
 
edm::InputTag smuonTags_
 
KFTrajectoryFitter * theFitter
 
KFTrajectoryFitter * theFitterOp
 
KFTrajectorySmoother * theSmoother
 
KFTrajectorySmoother * theSmootherOp
 
edm::ESHandle< GlobalTrackingGeometry > theTrackingGeometry
 
const GlobalTrackingGeometry * trackingGeometry_
 
edm::InputTag trackTags_
 
const TransientTrackingRecHitBuilder * TTRHBuilder
 
string txtOutFile_
 
edm::ESWatcher< GlobalPositionRcd > watchGlobalPositionRcd_
 
edm::ESWatcher< IdealMagneticFieldRecord > watchMagneticFieldRecord_
 
edm::ESWatcher< TrackingComponentsRecord > watchTrackingComponentsRecord_
 
edm::ESWatcher< GlobalTrackingGeometryRecord > watchTrackingGeometry_
 
bool writeDB_
 

Additional Inherited Members

- Public Types inherited from edm::EDAnalyzer
typedef EDAnalyzer ModuleType
 
- Public Types inherited from edm::EDConsumerBase
typedef ProductLabels Labels
 
- Static Public Member Functions inherited from edm::EDAnalyzer
static const std::string & baseType ()
 
static void fillDescriptions (ConfigurationDescriptions &descriptions)
 
static void prevalidate (ConfigurationDescriptions &)
 
static bool wantsGlobalLuminosityBlocks ()
 
static bool wantsGlobalRuns ()
 
static bool wantsInputProcessBlocks ()
 
static bool wantsProcessBlocks ()
 
static bool wantsStreamLuminosityBlocks ()
 
static bool wantsStreamRuns ()
 
- Protected Member Functions inherited from edm::EDConsumerBase
EDGetToken consumes (const TypeToGet &id, edm::InputTag const &tag)
 
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType > consumes (edm::InputTag const &tag)
 
template<BranchType B = InEvent>
EDConsumerBaseAdaptor< B > consumes (edm::InputTag tag) noexcept
 
template<BranchType B>
EDGetToken consumes (TypeToGet const &id, edm::InputTag const &tag)
 
ConsumesCollector consumesCollector ()
 Use a ConsumesCollector to gather consumes information from helper functions. More...
 
template<typename ProductType , BranchType B = InEvent>
void consumesMany ()
 
void consumesMany (const TypeToGet &id)
 
template<BranchType B>
void consumesMany (const TypeToGet &id)
 
template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto esConsumes ()
 
template<Transition Tr = Transition::Event>
constexpr auto esConsumes () noexcept
 
template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto esConsumes (ESInputTag const &tag)
 
template<Transition Tr = Transition::Event>
auto esConsumes (ESInputTag tag) noexcept
 
template<Transition Tr = Transition::Event>
ESGetTokenGeneric esConsumes (eventsetup::EventSetupRecordKey const &iRecord, eventsetup::DataKey const &iKey)
 Used with EventSetupRecord::doGet. More...
 
EDGetToken mayConsume (const TypeToGet &id, edm::InputTag const &tag)
 
template<BranchType B>
EDGetToken mayConsume (const TypeToGet &id, edm::InputTag const &tag)
 
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType > mayConsume (edm::InputTag const &tag)
 

Detailed Description

Description: Producer of relative tracker and muon system alignment

Implementation: A sample of global muons is used for the aligning tracker and muon system relatively as "rigid bodies", i.e. determining offset and rotation (6 numbers)

Definition at line 109 of file GlobalTrackerMuonAlignment.cc.

Constructor & Destructor Documentation

◆ GlobalTrackerMuonAlignment()

GlobalTrackerMuonAlignment::GlobalTrackerMuonAlignment ( const edm::ParameterSet &  iConfig)
explicit

Definition at line 272 of file GlobalTrackerMuonAlignment.cc.

273  : trackTags_(iConfig.getParameter<edm::InputTag>("tracks")),
274  muonTags_(iConfig.getParameter<edm::InputTag>("muons")),
275  gmuonTags_(iConfig.getParameter<edm::InputTag>("gmuons")),
276  smuonTags_(iConfig.getParameter<edm::InputTag>("smuons")),
277  propagator_(iConfig.getParameter<string>("Propagator")),
278 
279  cosmicMuonMode_(iConfig.getParameter<bool>("cosmics")),
280  isolatedMuonMode_(iConfig.getParameter<bool>("isolated")),
281 
282  refitMuon_(iConfig.getParameter<bool>("refitmuon")),
283  refitTrack_(iConfig.getParameter<bool>("refittrack")),
284 
285  rootOutFile_(iConfig.getUntrackedParameter<string>("rootOutFile")),
286  txtOutFile_(iConfig.getUntrackedParameter<string>("txtOutFile")),
287  writeDB_(iConfig.getUntrackedParameter<bool>("writeDB")),
288  debug_(iConfig.getUntrackedParameter<bool>("debug")),
289  defineFitter(true) {
290 #ifdef NO_FALSE_FALSE_MODE
291  if (cosmicMuonMode_ == false && isolatedMuonMode_ == false) {
292  throw cms::Exception("BadConfig") << "Muon collection not set! "
293  << "Use from GlobalTrackerMuonAlignment_test_cfg.py.";
294  }
295 #endif
296  if (cosmicMuonMode_ == true && isolatedMuonMode_ == true) {
297  throw cms::Exception("BadConfig") << "Muon collection can be either cosmic or isolated! "
298  << "Please set only one to true.";
299  }
300 }

References cosmicMuonMode_, Exception, and isolatedMuonMode_.

◆ ~GlobalTrackerMuonAlignment()

GlobalTrackerMuonAlignment::~GlobalTrackerMuonAlignment ( )
override

Definition at line 302 of file GlobalTrackerMuonAlignment.cc.

302  {
303  // do anything here that needs to be done at desctruction time
304  // (e.g. close files, deallocate resources etc.)
305 }

Member Function Documentation

◆ analyze()

void GlobalTrackerMuonAlignment::analyze ( const edm::Event &  iEvent,
const edm::EventSetup &  iSetup 
)
overrideprivatevirtual

Implements edm::EDAnalyzer.

Definition at line 312 of file GlobalTrackerMuonAlignment.cc.

312  {
313  //GlobalTrackerMuonAlignment::analyzeTrackTrack(iEvent, iSetup);
314  GlobalTrackerMuonAlignment::analyzeTrackTrajectory(iEvent, iSetup);
315 }

References analyzeTrackTrajectory(), and iEvent.

◆ analyzeTrackTrack()

void GlobalTrackerMuonAlignment::analyzeTrackTrack ( const edm::Event &  iEvent,
const edm::EventSetup &  iSetup 
)

Definition at line 573 of file GlobalTrackerMuonAlignment.cc.

573  {
574  using namespace edm;
575  using namespace reco;
576  //debug_ = true;
577  bool info = false;
578  bool alarm = false;
579  //bool alarm = true;
580  double PI = 3.1415927;
581 
582  cosmicMuonMode_ = true; // true: both Cosmic and IsolatedMuon are treated with 1,2 tracks
583  isolatedMuonMode_ = !cosmicMuonMode_; //true: only IsolatedMuon are treated with 1 track
584 
585  //if(N_event == 1)
586  //GlobalPositionRcdRead1::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup);
587 
588  N_event++;
589  if (info || debug_) {
590  std::cout << "----- event " << N_event << " -- tracks " << N_track << " ---";
591  if (cosmicMuonMode_)
592  std::cout << " treated as CosmicMu ";
593  if (isolatedMuonMode_)
594  std::cout << " treated as IsolatedMu ";
595  std::cout << std::endl;
596  }
597 
598  Handle<reco::TrackCollection> tracks;
599  Handle<reco::TrackCollection> muons;
600  Handle<reco::TrackCollection> gmuons;
601  Handle<reco::MuonCollection> smuons;
602 
603  if (collectionIsolated) {
604  iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "TrackerOnly", tracks);
605  iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "StandAlone", muons);
606  iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "GlobalMuon", gmuons);
607  iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "SelectedMuons", smuons);
608  } else if (collectionCosmic) {
609  iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "TrackerOnly", tracks);
610  iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "StandAlone", muons);
611  iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "GlobalMuon", gmuons);
612  iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "SelectedMuons", smuons);
613  } else {
614  iEvent.getByLabel(trackTags_, tracks);
615  iEvent.getByLabel(muonTags_, muons);
616  iEvent.getByLabel(gmuonTags_, gmuons);
617  iEvent.getByLabel(smuonTags_, smuons);
618  }
619 
620  if (debug_) {
621  std::cout << " ievBunch " << iEvent.bunchCrossing() << " runN " << (int)iEvent.run() << std::endl;
622  std::cout << " N tracks s/amu gmu selmu " << tracks->size() << " " << muons->size() << " " << gmuons->size() << " "
623  << smuons->size() << std::endl;
624  for (MuonCollection::const_iterator itMuon = smuons->begin(); itMuon != smuons->end(); ++itMuon) {
625  std::cout << " is isolatValid Matches " << itMuon->isIsolationValid() << " " << itMuon->isMatchesValid()
626  << std::endl;
627  }
628  }
629 
630  if (isolatedMuonMode_) { // ---- Only 1 Isolated Muon --------
631  if (tracks->size() != 1)
632  return;
633  if (muons->size() != 1)
634  return;
635  if (gmuons->size() != 1)
636  return;
637  if (smuons->size() != 1)
638  return;
639  }
640 
641  if (cosmicMuonMode_) { // ---- 1,2 Cosmic Muon --------
642  if (smuons->size() > 2)
643  return;
644  if (tracks->size() != smuons->size())
645  return;
646  if (muons->size() != smuons->size())
647  return;
648  if (gmuons->size() != smuons->size())
649  return;
650  }
651 
652  // ok mc_isolated_mu
653  //edm::ESHandle<TrackerGeometry> trackerGeometry;
654  //iSetup.get<TrackerDigiGeometryRecord>().get(trackerGeometry);
655  // ok mc_isolated_mu
656  //edm::ESHandle<DTGeometry> dtGeometry;
657  //iSetup.get<MuonGeometryRecord>().get(dtGeometry);
658  // don't work
659  //edm::ESHandle<CSCGeometry> cscGeometry;
660  //iSetup.get<MuonGeometryRecord>().get(cscGeometry);
661 
662  if (watchTrackingGeometry_.check(iSetup) || !trackingGeometry_) {
663  edm::ESHandle<GlobalTrackingGeometry> trackingGeometry;
664  iSetup.get<GlobalTrackingGeometryRecord>().get(trackingGeometry);
665  trackingGeometry_ = &*trackingGeometry;
666  }
667 
668  if (watchMagneticFieldRecord_.check(iSetup) || !magneticField_) {
669  edm::ESHandle<MagneticField> magneticField;
670  iSetup.get<IdealMagneticFieldRecord>().get(magneticField);
671  magneticField_ = &*magneticField;
672  }
673 
674  if (watchGlobalPositionRcd_.check(iSetup) || !globalPositionRcd_) {
675  edm::ESHandle<Alignments> globalPositionRcd;
676  iSetup.get<GlobalPositionRcd>().get(globalPositionRcd);
677  globalPositionRcd_ = &*globalPositionRcd;
678  for (std::vector<AlignTransform>::const_iterator i = globalPositionRcd_->m_align.begin();
679  i != globalPositionRcd_->m_align.end();
680  ++i) {
681  if (DetId(DetId::Tracker).rawId() == i->rawId())
682  iteratorTrackerRcd = i;
683  if (DetId(DetId::Muon).rawId() == i->rawId())
684  iteratorMuonRcd = i;
685  if (DetId(DetId::Ecal).rawId() == i->rawId())
686  iteratorEcalRcd = i;
687  if (DetId(DetId::Hcal).rawId() == i->rawId())
688  iteratorHcalRcd = i;
689  }
690  if (true || debug_) {
691  std::cout << "=== iteratorMuonRcd " << iteratorMuonRcd->rawId() << " ====\n"
692  << " translation " << iteratorMuonRcd->translation() << "\n"
693  << " angles " << iteratorMuonRcd->rotation().eulerAngles() << std::endl;
694  std::cout << iteratorMuonRcd->rotation() << std::endl;
695  }
696  } // end of GlobalPositionRcd
697 
698  ESHandle<Propagator> propagator;
699  iSetup.get<TrackingComponentsRecord>().get(propagator_, propagator);
700 
701  SteppingHelixPropagator alongStHePr = SteppingHelixPropagator(magneticField_, alongMomentum);
702  SteppingHelixPropagator oppositeStHePr = SteppingHelixPropagator(magneticField_, oppositeToMomentum);
703 
704  //double tolerance = 5.e-5;
705  defaultRKPropagator::Product aprod(magneticField_, alongMomentum, 5.e-5);
706  auto& alongRKPr = aprod.propagator;
707  defaultRKPropagator::Product oprod(magneticField_, oppositeToMomentum, 5.e-5);
708  auto& oppositeRKPr = oprod.propagator;
709 
710  float epsilon = 5.;
711  SmartPropagator alongSmPr = SmartPropagator(alongRKPr, alongStHePr, magneticField_, alongMomentum, epsilon);
712  SmartPropagator oppositeSmPr =
713  SmartPropagator(oppositeRKPr, oppositeStHePr, magneticField_, oppositeToMomentum, epsilon);
714 
715  // ................................................ selected/global muon
716  //itMuon --> Jim's globalMuon
717  for (MuonCollection::const_iterator itMuon = smuons->begin(); itMuon != smuons->end(); ++itMuon) {
718  if (debug_) {
719  std::cout << " mu gM is GM Mu SaM tM " << itMuon->isGlobalMuon() << " " << itMuon->isMuon() << " "
720  << itMuon->isStandAloneMuon() << " " << itMuon->isTrackerMuon() << " " << std::endl;
721  }
722 
723  // get information about the innermost muon hit -------------------------
724  TransientTrack muTT(itMuon->outerTrack(), magneticField_, trackingGeometry_);
725  TrajectoryStateOnSurface innerMuTSOS = muTT.innermostMeasurementState();
726  TrajectoryStateOnSurface outerMuTSOS = muTT.outermostMeasurementState();
727 
728  // get information about the outermost tracker hit -----------------------
729  TransientTrack trackTT(itMuon->track(), magneticField_, trackingGeometry_);
730  TrajectoryStateOnSurface outerTrackTSOS = trackTT.outermostMeasurementState();
731  TrajectoryStateOnSurface innerTrackTSOS = trackTT.innermostMeasurementState();
732 
733  GlobalPoint pointTrackIn = innerTrackTSOS.globalPosition();
734  GlobalPoint pointTrackOut = outerTrackTSOS.globalPosition();
735  float lenghtTrack = (pointTrackOut - pointTrackIn).mag();
736  GlobalPoint pointMuonIn = innerMuTSOS.globalPosition();
737  GlobalPoint pointMuonOut = outerMuTSOS.globalPosition();
738  float lenghtMuon = (pointMuonOut - pointMuonIn).mag();
739  GlobalVector momentumTrackOut = outerTrackTSOS.globalMomentum();
740  GlobalVector momentumTrackIn = innerTrackTSOS.globalMomentum();
741  float distanceInIn = (pointTrackIn - pointMuonIn).mag();
742  float distanceInOut = (pointTrackIn - pointMuonOut).mag();
743  float distanceOutIn = (pointTrackOut - pointMuonIn).mag();
744  float distanceOutOut = (pointTrackOut - pointMuonOut).mag();
745 
746  if (debug_) {
747  std::cout << " pointTrackIn " << pointTrackIn << std::endl;
748  std::cout << " Out " << pointTrackOut << " lenght " << lenghtTrack << std::endl;
749  std::cout << " point MuonIn " << pointMuonIn << std::endl;
750  std::cout << " Out " << pointMuonOut << " lenght " << lenghtMuon << std::endl;
751  std::cout << " momeTrackIn Out " << momentumTrackIn << " " << momentumTrackOut << std::endl;
752  std::cout << " doi io ii oo " << distanceOutIn << " " << distanceInOut << " " << distanceInIn << " "
753  << distanceOutOut << std::endl;
754  }
755 
756  if (lenghtTrack < 90.)
757  continue;
758  if (lenghtMuon < 300.)
759  continue;
760  if (momentumTrackIn.mag() < 15. || momentumTrackOut.mag() < 15.)
761  continue;
762  if (trackTT.charge() != muTT.charge())
763  continue;
764 
765  if (debug_)
766  std::cout << " passed lenght momentum cuts" << std::endl;
767 
768  GlobalVector GRm, GPm, Nl, Rm, Pm, Rt, Pt, Rt0;
769  AlgebraicSymMatrix66 Cm, C0, Ce, C1;
770 
771  TrajectoryStateOnSurface extrapolationT;
772  int tsosMuonIf = 0;
773 
774  if (isolatedMuonMode_) { //------------------------------- Isolated Muon -----
775  const Surface& refSurface = innerMuTSOS.surface();
776  ConstReferenceCountingPointer<TangentPlane> tpMuLocal(refSurface.tangentPlane(innerMuTSOS.localPosition()));
777  Nl = tpMuLocal->normalVector();
778  ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
779  GlobalVector Ng = tpMuGlobal->normalVector();
780  Surface* surf = (Surface*)&refSurface;
781  const Plane* refPlane = dynamic_cast<Plane*>(surf);
782  GlobalVector Nlp = refPlane->normalVector();
783 
784  if (debug_) {
785  std::cout << " Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
786  << " alfa " << int(asin(Nl.x()) * 57.29578) << std::endl;
787  std::cout << " global " << Ng.x() << " " << Ng.y() << " " << Ng.z() << std::endl;
788  std::cout << " lp " << Nlp.x() << " " << Nlp.y() << " " << Nlp.z() << std::endl;
789  //std::cout<<" Nlocal Nglobal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
790  // <<Ng.x()<<" "<<Ng.y()<<" "<<Ng.z()
791  //<<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
792  }
793 
794  // extrapolation to innermost muon hit
795  extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
796  //extrapolationT = propagator->propagate(outerTrackTSOS, refSurface);
797 
798  if (!extrapolationT.isValid()) {
799  if (false & alarm)
800  std::cout << " !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
801  //<<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
802  << std::endl;
803  continue;
804  }
805  tsosMuonIf = 1;
806  Rt = GlobalVector((extrapolationT.globalPosition()).x(),
807  (extrapolationT.globalPosition()).y(),
808  (extrapolationT.globalPosition()).z());
809 
810  Pt = extrapolationT.globalMomentum();
811  // global parameters of muon
812  GRm = GlobalVector(
813  (innerMuTSOS.globalPosition()).x(), (innerMuTSOS.globalPosition()).y(), (innerMuTSOS.globalPosition()).z());
814  GPm = innerMuTSOS.globalMomentum();
815  Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
816  (outerTrackTSOS.globalPosition()).y(),
817  (outerTrackTSOS.globalPosition()).z());
818  Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + innerMuTSOS.cartesianError().matrix());
819  C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
820  Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
821  C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
822 
823  } // ------------------------------- end Isolated Muon -----
824 
825  if (cosmicMuonMode_) { //------------------------------- Cosmic Muon -----
826 
827  if ((distanceOutIn <= distanceInOut) & (distanceOutIn <= distanceInIn) &
828  (distanceOutIn <= distanceOutOut)) { // ----- Out - In ------
829  if (debug_)
830  std::cout << " ----- Out - In -----" << std::endl;
831 
832  const Surface& refSurface = innerMuTSOS.surface();
833  ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
834  Nl = tpMuGlobal->normalVector();
835 
836  // extrapolation to innermost muon hit
837  extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
838  //extrapolationT = propagator->propagate(outerTrackTSOS, refSurface);
839 
840  if (!extrapolationT.isValid()) {
841  if (false & alarm)
842  std::cout << " !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
843  //<<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
844  << std::endl;
845  continue;
846  }
847  if (debug_)
848  std::cout << " extrapolationT.isValid " << extrapolationT.isValid() << std::endl;
849 
850  tsosMuonIf = 1;
851  Rt = GlobalVector((extrapolationT.globalPosition()).x(),
852  (extrapolationT.globalPosition()).y(),
853  (extrapolationT.globalPosition()).z());
854 
855  Pt = extrapolationT.globalMomentum();
856  // global parameters of muon
857  GRm = GlobalVector(
858  (innerMuTSOS.globalPosition()).x(), (innerMuTSOS.globalPosition()).y(), (innerMuTSOS.globalPosition()).z());
859  GPm = innerMuTSOS.globalMomentum();
860  Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
861  (outerTrackTSOS.globalPosition()).y(),
862  (outerTrackTSOS.globalPosition()).z());
863  Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + innerMuTSOS.cartesianError().matrix());
864  C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
865  Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
866  C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
867 
868  if (false & debug_) {
869  //std::cout<<" ->propDir "<<propagator->propagationDirection()<<std::endl;
870  PropagationDirectionChooser Chooser;
871  std::cout << " propDirCh " << Chooser.operator()(*outerTrackTSOS.freeState(), refSurface) << " Ch == along "
872  << (alongMomentum == Chooser.operator()(*outerTrackTSOS.freeState(), refSurface)) << std::endl;
873  std::cout << " --- Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
874  << " alfa " << int(asin(Nl.x()) * 57.29578) << std::endl;
875  }
876  } // enf of ---- Out - In -----
877 
878  else if ((distanceInOut <= distanceInIn) & (distanceInOut <= distanceOutIn) &
879  (distanceInOut <= distanceOutOut)) { // ----- In - Out ------
880  if (debug_)
881  std::cout << " ----- In - Out -----" << std::endl;
882 
883  const Surface& refSurface = outerMuTSOS.surface();
884  ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
885  Nl = tpMuGlobal->normalVector();
886 
887  // extrapolation to outermost muon hit
888  extrapolationT = oppositeSmPr.propagate(innerTrackTSOS, refSurface);
889 
890  if (!extrapolationT.isValid()) {
891  if (false & alarm)
892  std::cout << " !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
893  << "\a\a\a\a\a\a\a\a" << extrapolationT.isValid() << std::endl;
894  continue;
895  }
896  if (debug_)
897  std::cout << " extrapolationT.isValid " << extrapolationT.isValid() << std::endl;
898 
899  tsosMuonIf = 2;
900  Rt = GlobalVector((extrapolationT.globalPosition()).x(),
901  (extrapolationT.globalPosition()).y(),
902  (extrapolationT.globalPosition()).z());
903 
904  Pt = extrapolationT.globalMomentum();
905  // global parameters of muon
906  GRm = GlobalVector(
907  (outerMuTSOS.globalPosition()).x(), (outerMuTSOS.globalPosition()).y(), (outerMuTSOS.globalPosition()).z());
908  GPm = outerMuTSOS.globalMomentum();
909  Rt0 = GlobalVector((innerTrackTSOS.globalPosition()).x(),
910  (innerTrackTSOS.globalPosition()).y(),
911  (innerTrackTSOS.globalPosition()).z());
912  Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + outerMuTSOS.cartesianError().matrix());
913  C0 = AlgebraicSymMatrix66(innerTrackTSOS.cartesianError().matrix());
914  Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
915  C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
916 
917  if (false & debug_) {
918  //std::cout<<" ->propDir "<<propagator->propagationDirection()<<std::endl;
919  PropagationDirectionChooser Chooser;
920  std::cout << " propDirCh " << Chooser.operator()(*innerTrackTSOS.freeState(), refSurface)
921  << " Ch == oppisite "
922  << (oppositeToMomentum == Chooser.operator()(*innerTrackTSOS.freeState(), refSurface)) << std::endl;
923  std::cout << " --- Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
924  << " alfa " << int(asin(Nl.x()) * 57.29578) << std::endl;
925  }
926  } // enf of ---- In - Out -----
927 
928  else if ((distanceOutOut <= distanceInOut) & (distanceOutOut <= distanceInIn) &
929  (distanceOutOut <= distanceOutIn)) { // ----- Out - Out ------
930  if (debug_)
931  std::cout << " ----- Out - Out -----" << std::endl;
932 
933  // reject: momentum of track has opposite direction to muon track
934  continue;
935 
936  const Surface& refSurface = outerMuTSOS.surface();
937  ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
938  Nl = tpMuGlobal->normalVector();
939 
940  // extrapolation to outermost muon hit
941  extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
942 
943  if (!extrapolationT.isValid()) {
944  if (alarm)
945  std::cout << " !!!!!!!!! Catastrophe Out-Out extrapolationT.isValid "
946  << "\a\a\a\a\a\a\a\a" << extrapolationT.isValid() << std::endl;
947  continue;
948  }
949  if (debug_)
950  std::cout << " extrapolationT.isValid " << extrapolationT.isValid() << std::endl;
951 
952  tsosMuonIf = 2;
953  Rt = GlobalVector((extrapolationT.globalPosition()).x(),
954  (extrapolationT.globalPosition()).y(),
955  (extrapolationT.globalPosition()).z());
956 
957  Pt = extrapolationT.globalMomentum();
958  // global parameters of muon
959  GRm = GlobalVector(
960  (outerMuTSOS.globalPosition()).x(), (outerMuTSOS.globalPosition()).y(), (outerMuTSOS.globalPosition()).z());
961  GPm = outerMuTSOS.globalMomentum();
962  Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
963  (outerTrackTSOS.globalPosition()).y(),
964  (outerTrackTSOS.globalPosition()).z());
965  Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + outerMuTSOS.cartesianError().matrix());
966  C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
967  Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
968  C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
969 
970  if (debug_) {
971  PropagationDirectionChooser Chooser;
972  std::cout << " propDirCh " << Chooser.operator()(*outerTrackTSOS.freeState(), refSurface) << " Ch == along "
973  << (alongMomentum == Chooser.operator()(*outerTrackTSOS.freeState(), refSurface)) << std::endl;
974  std::cout << " --- Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
975  << " alfa " << int(asin(Nl.x()) * 57.29578) << std::endl;
976  std::cout << " Nornal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
977  << " alfa " << int(asin(Nl.x()) * 57.29578) << std::endl;
978  }
979  } // enf of ---- Out - Out -----
980  else {
981  if (alarm)
982  std::cout << " ------- !!!!!!!!!!!!!!! No proper Out - In \a\a\a\a\a\a\a" << std::endl;
983  continue;
984  }
985 
986  } // ------------------------------- end Cosmic Muon -----
987 
988  if (tsosMuonIf == 0) {
989  if (info) {
990  std::cout << "No tsosMuon !!!!!!" << std::endl;
991  continue;
992  }
993  }
994  TrajectoryStateOnSurface tsosMuon;
995  if (tsosMuonIf == 1)
996  tsosMuon = muTT.innermostMeasurementState();
997  else
998  tsosMuon = muTT.outermostMeasurementState();
999 
1000  //GlobalTrackerMuonAlignment::misalignMuon(GRm, GPm, Nl, Rt, Rm, Pm);
1001  AlgebraicVector4 LPRm; // muon local (dx/dz, dy/dz, x, y)
1002  GlobalTrackerMuonAlignment::misalignMuonL(GRm, GPm, Nl, Rt, Rm, Pm, LPRm, extrapolationT, tsosMuon);
1003 
1004  GlobalVector resR = Rm - Rt;
1005  GlobalVector resP0 = Pm - Pt;
1006  GlobalVector resP = Pm / Pm.mag() - Pt / Pt.mag();
1007  float RelMomResidual = (Pm.mag() - Pt.mag()) / (Pt.mag() + 1.e-6);
1008  ;
1009 
1010  AlgebraicVector6 Vm;
1011  Vm(0) = resR.x();
1012  Vm(1) = resR.y();
1013  Vm(2) = resR.z();
1014  Vm(3) = resP0.x();
1015  Vm(4) = resP0.y();
1016  Vm(5) = resP0.z();
1017  float Rmuon = Rm.perp();
1018  float Zmuon = Rm.z();
1019  float alfa_x = atan2(Nl.x(), Nl.y()) * 57.29578;
1020 
1021  if (debug_) {
1022  std::cout << " Nx Ny Nz alfa_x " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " " << alfa_x << std::endl;
1023  //std::cout<<" Rm "<<Rm<<std::endl<<" Rt "<<Rt<<std::endl<<" resR "<<resR<<std::endl
1024  // <<" resP "<<resP<<" dp/p "<<RelMomResidual<<std::endl;
1025  }
1026 
1027  double chi_d = 0;
1028  for (int i = 0; i <= 5; i++)
1029  chi_d += Vm(i) * Vm(i) / Cm(i, i);
1030 
1031  AlgebraicVector5 Vml(tsosMuon.localParameters().vector() - extrapolationT.localParameters().vector());
1032  AlgebraicSymMatrix55 m(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
1033  AlgebraicSymMatrix55 Cml(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
1034  bool ierrLoc = !m.Invert();
1035  if (ierrLoc && debug_ && info) {
1036  std::cout << " ==== Error inverting Local covariance matrix ==== " << std::endl;
1037  continue;
1038  }
1039  double chi_Loc = ROOT::Math::Similarity(Vml, m);
1040  if (debug_)
1041  std::cout << " chi_Loc px/pz/err " << chi_Loc << " " << Vml(1) / sqrt(Cml(1, 1)) << std::endl;
1042 
1043  if (Pt.mag() < 15.)
1044  continue;
1045  if (Pm.mag() < 5.)
1046  continue;
1047 
1048  //if(Pt.mag() < 30.) continue; // momenum cut < 30GeV
1049  //if(Pt.mag() < 60.) continue; // momenum cut < 30GeV
1050  //if(Pt.mag() > 50.) continue; // momenum cut > 50GeV
1051  //if(Pt.mag() > 100.) continue; // momenum cut > 100GeV
1052  //if(trackTT.charge() < 0) continue; // select positive charge
1053  //if(trackTT.charge() > 0) continue; // select negative charge
1054 
1055  //if(fabs(resR.x()) > 5.) continue; // strong cut X
1056  //if(fabs(resR.y()) > 5.) continue; // Y
1057  //if(fabs(resR.z()) > 5.) continue; // Z
1058  //if(fabs(resR.mag()) > 7.5) continue; // dR
1059 
1060  /*
1061  //if(fabs(RelMomResidual) > 0.5) continue;
1062  if(fabs(resR.x()) > 20.) continue;
1063  if(fabs(resR.y()) > 20.) continue;
1064  if(fabs(resR.z()) > 20.) continue;
1065  if(fabs(resR.mag()) > 30.) continue;
1066  if(fabs(resP.x()) > 0.06) continue;
1067  if(fabs(resP.y()) > 0.06) continue;
1068  if(fabs(resP.z()) > 0.06) continue;
1069  if(chi_d > 40.) continue;
1070  */
1071 
1072  // select Barrel
1073  //if(Rmuon < 400. || Rmuon > 450.) continue;
1074  //if(Zmuon < -600. || Zmuon > 600.) continue;
1075  //if(fabs(Nl.z()) > 0.95) continue;
1076  //MuSelect = " Barrel";
1077  // EndCap1
1078  //if(Rmuon < 120. || Rmuon > 450.) continue;
1079  //if(Zmuon < -720.) continue;
1080  //if(Zmuon > -580.) continue;
1081  //if(fabs(Nl.z()) < 0.95) continue;
1082  //MuSelect = " EndCap1";
1083  // EndCap2
1084  //if(Rmuon < 120. || Rmuon > 450.) continue;
1085  //if(Zmuon > 720.) continue;
1086  //if(Zmuon < 580.) continue;
1087  //if(fabs(Nl.z()) < 0.95) continue;
1088  //MuSelect = " EndCap2";
1089  // select All
1090  if (Rmuon < 120. || Rmuon > 450.)
1091  continue;
1092  if (Zmuon < -720. || Zmuon > 720.)
1093  continue;
1094  MuSelect = " Barrel+EndCaps";
1095 
1096  if (debug_)
1097  std::cout << " .............. passed all cuts" << std::endl;
1098 
1099  N_track++;
1100  // gradient and Hessian for each track
1101 
1102  GlobalTrackerMuonAlignment::gradientGlobalAlg(Rt, Pt, Rm, Nl, Cm);
1103  GlobalTrackerMuonAlignment::gradientGlobal(Rt, Pt, Rm, Pm, Nl, Cm);
1104 
1105  CLHEP::HepSymMatrix covLoc(4, 0);
1106  for (int i = 1; i <= 4; i++)
1107  for (int j = 1; j <= i; j++) {
1108  covLoc(i, j) = (tsosMuon.localError().matrix() + extrapolationT.localError().matrix())(i, j);
1109  //if(i != j) Cov(i,j) = 0.;
1110  }
1111 
1112  const Surface& refSurface = tsosMuon.surface();
1113  CLHEP::HepMatrix rotLoc(3, 3, 0);
1114  rotLoc(1, 1) = refSurface.rotation().xx();
1115  rotLoc(1, 2) = refSurface.rotation().xy();
1116  rotLoc(1, 3) = refSurface.rotation().xz();
1117 
1118  rotLoc(2, 1) = refSurface.rotation().yx();
1119  rotLoc(2, 2) = refSurface.rotation().yy();
1120  rotLoc(2, 3) = refSurface.rotation().yz();
1121 
1122  rotLoc(3, 1) = refSurface.rotation().zx();
1123  rotLoc(3, 2) = refSurface.rotation().zy();
1124  rotLoc(3, 3) = refSurface.rotation().zz();
1125 
1126  CLHEP::HepVector posLoc(3);
1127  posLoc(1) = refSurface.position().x();
1128  posLoc(2) = refSurface.position().y();
1129  posLoc(3) = refSurface.position().z();
1130 
1131  GlobalTrackerMuonAlignment::gradientLocal(Rt, Pt, Rm, Pm, Nl, covLoc, rotLoc, posLoc, LPRm);
1132 
1133  if (debug_) {
1134  std::cout << " Norm " << Nl << std::endl;
1135  std::cout << " posLoc " << posLoc.T() << std::endl;
1136  std::cout << " rotLoc " << rotLoc << std::endl;
1137  }
1138 
1139  // ----------------------------------------------------- fill histogram
1140  histo->Fill(itMuon->track()->pt());
1141 
1142  //histo2->Fill(itMuon->track()->outerP());
1143  histo2->Fill(Pt.mag());
1144  histo3->Fill((PI / 2. - itMuon->track()->outerTheta()));
1145  histo4->Fill(itMuon->track()->phi());
1146  histo5->Fill(Rmuon);
1147  histo6->Fill(Zmuon);
1148  histo7->Fill(RelMomResidual);
1149  //histo8->Fill(chi);
1150  histo8->Fill(chi_d);
1151 
1152  histo101->Fill(Zmuon, Rmuon);
1153  histo101->Fill(Rt0.z(), Rt0.perp());
1154  histo102->Fill(Rt0.x(), Rt0.y());
1155  histo102->Fill(Rm.x(), Rm.y());
1156 
1157  histo11->Fill(resR.mag());
1158  if (fabs(Nl.x()) < 0.98)
1159  histo12->Fill(resR.x());
1160  if (fabs(Nl.y()) < 0.98)
1161  histo13->Fill(resR.y());
1162  if (fabs(Nl.z()) < 0.98)
1163  histo14->Fill(resR.z());
1164  histo15->Fill(resP.x());
1165  histo16->Fill(resP.y());
1166  histo17->Fill(resP.z());
1167 
1168  if ((fabs(Nl.x()) < 0.98) && (fabs(Nl.y()) < 0.98) && (fabs(Nl.z()) < 0.98)) {
1169  histo18->Fill(sqrt(C0(0, 0)));
1170  histo19->Fill(sqrt(C1(0, 0)));
1171  histo20->Fill(sqrt(C1(0, 0) + Ce(0, 0)));
1172  }
1173  if (fabs(Nl.x()) < 0.98)
1174  histo21->Fill(Vm(0) / sqrt(Cm(0, 0)));
1175  if (fabs(Nl.y()) < 0.98)
1176  histo22->Fill(Vm(1) / sqrt(Cm(1, 1)));
1177  if (fabs(Nl.z()) < 0.98)
1178  histo23->Fill(Vm(2) / sqrt(Cm(2, 2)));
1179  histo24->Fill(Vm(3) / sqrt(C1(3, 3) + Ce(3, 3)));
1180  histo25->Fill(Vm(4) / sqrt(C1(4, 4) + Ce(4, 4)));
1181  histo26->Fill(Vm(5) / sqrt(C1(5, 5) + Ce(5, 5)));
1182  histo27->Fill(Nl.x());
1183  histo28->Fill(Nl.y());
1184  histo29->Fill(lenghtTrack);
1185  histo30->Fill(lenghtMuon);
1186  histo31->Fill(chi_Loc);
1187  histo32->Fill(Vml(1) / sqrt(Cml(1, 1)));
1188  histo33->Fill(Vml(2) / sqrt(Cml(2, 2)));
1189  histo34->Fill(Vml(3) / sqrt(Cml(3, 3)));
1190  histo35->Fill(Vml(4) / sqrt(Cml(4, 4)));
1191 
1192  if (debug_) { //--------------------------------- debug print ----------
1193 
1194  std::cout << " diag 0[ " << C0(0, 0) << " " << C0(1, 1) << " " << C0(2, 2) << " " << C0(3, 3) << " " << C0(4, 4)
1195  << " " << C0(5, 5) << " ]" << std::endl;
1196  std::cout << " diag e[ " << Ce(0, 0) << " " << Ce(1, 1) << " " << Ce(2, 2) << " " << Ce(3, 3) << " " << Ce(4, 4)
1197  << " " << Ce(5, 5) << " ]" << std::endl;
1198  std::cout << " diag 1[ " << C1(0, 0) << " " << C1(1, 1) << " " << C1(2, 2) << " " << C1(3, 3) << " " << C1(4, 4)
1199  << " " << C1(5, 5) << " ]" << std::endl;
1200  std::cout << " Rm " << Rm.x() << " " << Rm.y() << " " << Rm.z() << " Pm " << Pm.x() << " " << Pm.y() << " "
1201  << Pm.z() << std::endl;
1202  std::cout << " Rt " << Rt.x() << " " << Rt.y() << " " << Rt.z() << " Pt " << Pt.x() << " " << Pt.y() << " "
1203  << Pt.z() << std::endl;
1204  std::cout << " Nl*(Rm-Rt) " << Nl.dot(Rm - Rt) << std::endl;
1205  std::cout << " resR " << resR.x() << " " << resR.y() << " " << resR.z() << " resP " << resP.x() << " " << resP.y()
1206  << " " << resP.z() << std::endl;
1207  std::cout << " Rm-t " << (Rm - Rt).x() << " " << (Rm - Rt).y() << " " << (Rm - Rt).z() << " Pm-t "
1208  << (Pm - Pt).x() << " " << (Pm - Pt).y() << " " << (Pm - Pt).z() << std::endl;
1209  std::cout << " Vm " << Vm << std::endl;
1210  std::cout << " +- " << sqrt(Cm(0, 0)) << " " << sqrt(Cm(1, 1)) << " " << sqrt(Cm(2, 2)) << " " << sqrt(Cm(3, 3))
1211  << " " << sqrt(Cm(4, 4)) << " " << sqrt(Cm(5, 5)) << std::endl;
1212  std::cout << " Pmuon Ptrack dP/Ptrack " << itMuon->outerTrack()->p() << " " << itMuon->track()->outerP() << " "
1213  << (itMuon->outerTrack()->p() - itMuon->track()->outerP()) / itMuon->track()->outerP() << std::endl;
1214  std::cout << " cov matrix " << std::endl;
1215  std::cout << Cm << std::endl;
1216  std::cout << " diag [ " << Cm(0, 0) << " " << Cm(1, 1) << " " << Cm(2, 2) << " " << Cm(3, 3) << " " << Cm(4, 4)
1217  << " " << Cm(5, 5) << " ]" << std::endl;
1218 
1219  AlgebraicSymMatrix66 Ro;
1220  double Diag[6];
1221  for (int i = 0; i <= 5; i++)
1222  Diag[i] = sqrt(Cm(i, i));
1223  for (int i = 0; i <= 5; i++)
1224  for (int j = 0; j <= 5; j++)
1225  Ro(i, j) = Cm(i, j) / Diag[i] / Diag[j];
1226  std::cout << " correlation matrix " << std::endl;
1227  std::cout << Ro << std::endl;
1228 
1229  AlgebraicSymMatrix66 CmI;
1230  for (int i = 0; i <= 5; i++)
1231  for (int j = 0; j <= 5; j++)
1232  CmI(i, j) = Cm(i, j);
1233 
1234  bool ierr = !CmI.Invert();
1235  if (ierr) {
1236  if (alarm || debug_)
1237  std::cout << " Error inverse covariance matrix !!!!!!!!!!!" << std::endl;
1238  continue;
1239  }
1240  std::cout << " inverse cov matrix " << std::endl;
1241  std::cout << Cm << std::endl;
1242 
1243  double chi = ROOT::Math::Similarity(Vm, CmI);
1244  std::cout << " chi chi_d " << chi << " " << chi_d << std::endl;
1245  } // end of debug_ printout --------------------------------------------
1246 
1247  } // end loop on selected muons, i.e. Jim's globalMuon
1248 
1249 } //end of analyzeTrackTrack

References alongMomentum, TrajectoryStateOnSurface::cartesianError(), edm::ESWatcher< T >::check(), collectionCosmic, collectionIsolated, cosmicMuonMode_, gather_cfg::cout, debug_, Vector3DBase< T, FrameTag >::dot(), MillePedeFileConverter_cfg::e, DetId::Ecal, geometryDiff::epsilon, TrajectoryStateOnSurface::freeState(), edm::EventSetup::get(), get, TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), globalPositionRcd_, gmuonTags_, gradientGlobal(), gradientGlobalAlg(), gradientLocal(), DetId::Hcal, histo, histo101, histo102, histo11, histo12, histo13, histo14, histo15, histo16, histo17, histo18, histo19, histo2, histo20, histo21, histo22, histo23, histo24, histo25, histo26, histo27, histo28, histo29, histo3, histo30, histo31, histo32, histo33, histo34, histo35, histo4, histo5, histo6, histo7, histo8, mps_fire::i, iEvent, info(), createfilelist::int, isolatedMuonMode_, TrajectoryStateOnSurface::isValid(), iteratorEcalRcd, iteratorHcalRcd, iteratorMuonRcd, iteratorTrackerRcd, dqmiolumiharvest::j, TrajectoryStateOnSurface::localError(), TrajectoryStateOnSurface::localParameters(), visualization-live-secondInstance_cfg::m, Alignments::m_align, PV3DBase< T, PVType, FrameType >::mag(), mag(), HLT_FULL_cff::magneticField, magneticField_, CartesianTrajectoryError::matrix(), LocalTrajectoryError::matrix(), misalignMuonL(), DetId::Muon, PDWG_BPHSkim_cff::muons, muonTags_, MuSelect, N_event, N_track, Plane::normalVector(), oppositeToMomentum, PV3DBase< T, PVType, FrameType >::perp(), GloballyPositioned< T >::position(), SmartPropagator::propagate(), TrackCandidateProducer_cfi::propagator, defaultRKPropagator::Product::propagator, propagator_, DetId::rawId(), GloballyPositioned< T >::rotation(), SmartPropagator_cfi::SmartPropagator, smuonTags_, mathSSE::sqrt(), SteppingHelixPropagator_cfi::SteppingHelixPropagator, TrajectoryStateOnSurface::surface(), Surface::tangentPlane(), DetId::Tracker, trackingGeometry_, PDWG_EXOHSCP_cff::tracks, trackTags_, LocalTrajectoryParameters::vector(), watchGlobalPositionRcd_, watchMagneticFieldRecord_, watchTrackingGeometry_, x, PV3DBase< T, PVType, FrameType >::x(), TkRotation< T >::xx(), TkRotation< T >::xy(), TkRotation< T >::xz(), y, PV3DBase< T, PVType, FrameType >::y(), TkRotation< T >::yx(), TkRotation< T >::yy(), TkRotation< T >::yz(), z, PV3DBase< T, PVType, FrameType >::z(), TkRotation< T >::zx(), TkRotation< T >::zy(), and TkRotation< T >::zz().

◆ analyzeTrackTrajectory()

void GlobalTrackerMuonAlignment::analyzeTrackTrajectory ( const edm::Event &  iEvent,
const edm::EventSetup &  iSetup 
)

Definition at line 1252 of file GlobalTrackerMuonAlignment.cc.

1252  {
1253  using namespace edm;
1254  using namespace reco;
1255  //debug_ = true;
1256  bool info = false;
1257  bool alarm = false;
1258  //bool alarm = true;
1259  double PI = 3.1415927;
1260 
1261  //-*-*-*-*-*-*-
1262  cosmicMuonMode_ = true; // true: both Cosmic and IsolatedMuon are treated with 1,2 tracks
1263  //cosmicMuonMode_ = false; // true: both Cosmic and IsolatedMuon are treated with 1,2 tracks
1264  //-*-*-*-*-*-*-
1265  isolatedMuonMode_ = !cosmicMuonMode_; //true: only IsolatedMuon are treated with 1 track
1266 
1267  N_event++;
1268  if (info || debug_) {
1269  std::cout << "----- event " << N_event << " -- tracks " << N_track << " ---";
1270  if (cosmicMuonMode_)
1271  std::cout << " treated as CosmicMu ";
1272  if (isolatedMuonMode_)
1273  std::cout << " treated as IsolatedMu ";
1274  std::cout << std::endl;
1275  }
1276 
1277  Handle<reco::TrackCollection> tracks;
1278  Handle<reco::TrackCollection> muons;
1279  Handle<reco::TrackCollection> gmuons;
1280  Handle<reco::MuonCollection> smuons;
1281 
1282  if (collectionIsolated) {
1283  iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "TrackerOnly", tracks);
1284  iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "StandAlone", muons);
1285  iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "GlobalMuon", gmuons);
1286  iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "SelectedMuons", smuons);
1287  } else if (collectionCosmic) {
1288  iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "TrackerOnly", tracks);
1289  iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "StandAlone", muons);
1290  iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "GlobalMuon", gmuons);
1291  iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "SelectedMuons", smuons);
1292  } else {
1293  iEvent.getByLabel(trackTags_, tracks);
1294  iEvent.getByLabel(muonTags_, muons);
1295  iEvent.getByLabel(gmuonTags_, gmuons);
1296  iEvent.getByLabel(smuonTags_, smuons);
1297  }
1298 
1299  if (debug_) {
1300  std::cout << " ievBunch " << iEvent.bunchCrossing() << " runN " << (int)iEvent.run() << std::endl;
1301  std::cout << " N tracks s/amu gmu selmu " << tracks->size() << " " << muons->size() << " " << gmuons->size() << " "
1302  << smuons->size() << std::endl;
1303  for (MuonCollection::const_iterator itMuon = smuons->begin(); itMuon != smuons->end(); ++itMuon) {
1304  std::cout << " is isolatValid Matches " << itMuon->isIsolationValid() << " " << itMuon->isMatchesValid()
1305  << std::endl;
1306  }
1307  }
1308 
1309  if (isolatedMuonMode_) { // ---- Only 1 Isolated Muon --------
1310  if (tracks->size() != 1)
1311  return;
1312  if (muons->size() != 1)
1313  return;
1314  if (gmuons->size() != 1)
1315  return;
1316  if (smuons->size() != 1)
1317  return;
1318  }
1319 
1320  if (cosmicMuonMode_) { // ---- 1,2 Cosmic Muon --------
1321  if (smuons->size() > 2)
1322  return;
1323  if (tracks->size() != smuons->size())
1324  return;
1325  if (muons->size() != smuons->size())
1326  return;
1327  if (gmuons->size() != smuons->size())
1328  return;
1329  }
1330 
1331  // ok mc_isolated_mu
1332  //edm::ESHandle<TrackerGeometry> trackerGeometry;
1333  //iSetup.get<TrackerDigiGeometryRecord>().get(trackerGeometry);
1334  // ok mc_isolated_mu
1335  //edm::ESHandle<DTGeometry> dtGeometry;
1336  //iSetup.get<MuonGeometryRecord>().get(dtGeometry);
1337  // don't work
1338  //edm::ESHandle<CSCGeometry> cscGeometry;
1339  //iSetup.get<MuonGeometryRecord>().get(cscGeometry);
1340 
1341  if (watchTrackingGeometry_.check(iSetup) || !trackingGeometry_) {
1342  edm::ESHandle<GlobalTrackingGeometry> trackingGeometry;
1343  iSetup.get<GlobalTrackingGeometryRecord>().get(trackingGeometry);
1344  trackingGeometry_ = &*trackingGeometry;
1345  theTrackingGeometry = trackingGeometry;
1346  }
1347 
1348  if (watchMagneticFieldRecord_.check(iSetup) || !magneticField_) {
1349  edm::ESHandle<MagneticField> magneticField;
1350  iSetup.get<IdealMagneticFieldRecord>().get(magneticField);
1351  magneticField_ = &*magneticField;
1352  }
1353 
1354  if (watchGlobalPositionRcd_.check(iSetup) || !globalPositionRcd_) {
1355  edm::ESHandle<Alignments> globalPositionRcd;
1356  iSetup.get<GlobalPositionRcd>().get(globalPositionRcd);
1357  globalPositionRcd_ = &*globalPositionRcd;
1358  for (std::vector<AlignTransform>::const_iterator i = globalPositionRcd_->m_align.begin();
1359  i != globalPositionRcd_->m_align.end();
1360  ++i) {
1361  if (DetId(DetId::Tracker).rawId() == i->rawId())
1362  iteratorTrackerRcd = i;
1363  if (DetId(DetId::Muon).rawId() == i->rawId())
1364  iteratorMuonRcd = i;
1365  if (DetId(DetId::Ecal).rawId() == i->rawId())
1366  iteratorEcalRcd = i;
1367  if (DetId(DetId::Hcal).rawId() == i->rawId())
1368  iteratorHcalRcd = i;
1369  }
1370  if (debug_) {
1371  std::cout << "=== iteratorTrackerRcd " << iteratorTrackerRcd->rawId() << " ====\n"
1372  << " translation " << iteratorTrackerRcd->translation() << "\n"
1373  << " angles " << iteratorTrackerRcd->rotation().eulerAngles() << std::endl;
1374  std::cout << iteratorTrackerRcd->rotation() << std::endl;
1375  std::cout << "=== iteratorMuonRcd " << iteratorMuonRcd->rawId() << " ====\n"
1376  << " translation " << iteratorMuonRcd->translation() << "\n"
1377  << " angles " << iteratorMuonRcd->rotation().eulerAngles() << std::endl;
1378  std::cout << iteratorMuonRcd->rotation() << std::endl;
1379  }
1380  } // end of GlobalPositionRcd
1381 
1382  ESHandle<Propagator> propagator;
1383  iSetup.get<TrackingComponentsRecord>().get(propagator_, propagator);
1384 
1385  SteppingHelixPropagator alongStHePr = SteppingHelixPropagator(magneticField_, alongMomentum);
1386  SteppingHelixPropagator oppositeStHePr = SteppingHelixPropagator(magneticField_, oppositeToMomentum);
1387 
1388  //double tolerance = 5.e-5;
1389  defaultRKPropagator::Product aprod(magneticField_, alongMomentum, 5.e-5);
1390  auto& alongRKPr = aprod.propagator;
1391  defaultRKPropagator::Product oprod(magneticField_, oppositeToMomentum, 5.e-5);
1392  auto& oppositeRKPr = oprod.propagator;
1393 
1394  float epsilon = 5.;
1395  SmartPropagator alongSmPr = SmartPropagator(alongRKPr, alongStHePr, magneticField_, alongMomentum, epsilon);
1396  SmartPropagator oppositeSmPr =
1397  SmartPropagator(oppositeRKPr, oppositeStHePr, magneticField_, oppositeToMomentum, epsilon);
1398 
1399  if (defineFitter) {
1400  if (debug_)
1401  std::cout << " ............... DEFINE FITTER ..................." << std::endl;
1402  KFUpdator* theUpdator = new KFUpdator();
1403  //Chi2MeasurementEstimator* theEstimator = new Chi2MeasurementEstimator(30);
1404  Chi2MeasurementEstimator* theEstimator = new Chi2MeasurementEstimator(100000, 100000);
1405  theFitter = new KFTrajectoryFitter(alongSmPr, *theUpdator, *theEstimator);
1406  theSmoother = new KFTrajectorySmoother(alongSmPr, *theUpdator, *theEstimator);
1407  theFitterOp = new KFTrajectoryFitter(oppositeSmPr, *theUpdator, *theEstimator);
1408  theSmootherOp = new KFTrajectorySmoother(oppositeSmPr, *theUpdator, *theEstimator);
1409 
1410  edm::ESHandle<TransientTrackingRecHitBuilder> builder;
1411  iSetup.get<TransientRecHitRecord>().get("WithTrackAngle", builder);
1412  this->TTRHBuilder = &(*builder);
1413  this->MuRHBuilder = new MuonTransientTrackingRecHitBuilder(theTrackingGeometry);
1414  if (debug_) {
1415  std::cout << " theTrackingGeometry.isValid() " << theTrackingGeometry.isValid() << std::endl;
1416  std::cout << "get also the MuonTransientTrackingRecHitBuilder"
1417  << "\n";
1418  std::cout << "get also the TransientTrackingRecHitBuilder"
1419  << "\n";
1420  }
1421  defineFitter = false;
1422  }
1423 
1424  // ................................................ selected/global muon
1425  //itMuon --> Jim's globalMuon
1426  for (MuonCollection::const_iterator itMuon = smuons->begin(); itMuon != smuons->end(); ++itMuon) {
1427  if (debug_) {
1428  std::cout << " mu gM is GM Mu SaM tM " << itMuon->isGlobalMuon() << " " << itMuon->isMuon() << " "
1429  << itMuon->isStandAloneMuon() << " " << itMuon->isTrackerMuon() << " " << std::endl;
1430  }
1431 
1432  // get information about the innermost muon hit -------------------------
1433  TransientTrack muTT(itMuon->outerTrack(), magneticField_, trackingGeometry_);
1434  TrajectoryStateOnSurface innerMuTSOS = muTT.innermostMeasurementState();
1435  TrajectoryStateOnSurface outerMuTSOS = muTT.outermostMeasurementState();
1436 
1437  // get information about the outermost tracker hit -----------------------
1438  TransientTrack trackTT(itMuon->track(), magneticField_, trackingGeometry_);
1439  TrajectoryStateOnSurface outerTrackTSOS = trackTT.outermostMeasurementState();
1440  TrajectoryStateOnSurface innerTrackTSOS = trackTT.innermostMeasurementState();
1441 
1442  GlobalPoint pointTrackIn = innerTrackTSOS.globalPosition();
1443  GlobalPoint pointTrackOut = outerTrackTSOS.globalPosition();
1444  float lenghtTrack = (pointTrackOut - pointTrackIn).mag();
1445  GlobalPoint pointMuonIn = innerMuTSOS.globalPosition();
1446  GlobalPoint pointMuonOut = outerMuTSOS.globalPosition();
1447  float lenghtMuon = (pointMuonOut - pointMuonIn).mag();
1448  GlobalVector momentumTrackOut = outerTrackTSOS.globalMomentum();
1449  GlobalVector momentumTrackIn = innerTrackTSOS.globalMomentum();
1450  float distanceInIn = (pointTrackIn - pointMuonIn).mag();
1451  float distanceInOut = (pointTrackIn - pointMuonOut).mag();
1452  float distanceOutIn = (pointTrackOut - pointMuonIn).mag();
1453  float distanceOutOut = (pointTrackOut - pointMuonOut).mag();
1454 
1455  if (debug_) {
1456  std::cout << " pointTrackIn " << pointTrackIn << std::endl;
1457  std::cout << " Out " << pointTrackOut << " lenght " << lenghtTrack << std::endl;
1458  std::cout << " point MuonIn " << pointMuonIn << std::endl;
1459  std::cout << " Out " << pointMuonOut << " lenght " << lenghtMuon << std::endl;
1460  std::cout << " momeTrackIn Out " << momentumTrackIn << " " << momentumTrackOut << std::endl;
1461  std::cout << " doi io ii oo " << distanceOutIn << " " << distanceInOut << " " << distanceInIn << " "
1462  << distanceOutOut << std::endl;
1463  }
1464 
1465  if (lenghtTrack < 90.)
1466  continue;
1467  if (lenghtMuon < 300.)
1468  continue;
1469  if (innerMuTSOS.globalMomentum().mag() < 5. || outerMuTSOS.globalMomentum().mag() < 5.)
1470  continue;
1471  if (momentumTrackIn.mag() < 15. || momentumTrackOut.mag() < 15.)
1472  continue;
1473  if (trackTT.charge() != muTT.charge())
1474  continue;
1475 
1476  if (debug_)
1477  std::cout << " passed lenght momentum cuts" << std::endl;
1478 
1479  GlobalVector GRm, GPm, Nl, Rm, Pm, Rt, Pt, Rt0;
1480  AlgebraicSymMatrix66 Cm, C0, Ce, C1;
1481 
1482  //extrapolationT = propagator->propagate(outerTrackTSOS, refSurface);
1483  TrajectoryStateOnSurface extrapolationT;
1484  int tsosMuonIf = 0;
1485 
1486  TrajectoryStateOnSurface muonFittedTSOS;
1487  TrajectoryStateOnSurface trackFittedTSOS;
1488 
1489  if (isolatedMuonMode_) { //------------------------------- Isolated Muon --- Out-In --
1490  if (debug_)
1491  std::cout << " ------ Isolated (out-in) ----- " << std::endl;
1492  const Surface& refSurface = innerMuTSOS.surface();
1493  ConstReferenceCountingPointer<TangentPlane> tpMuLocal(refSurface.tangentPlane(innerMuTSOS.localPosition()));
1494  Nl = tpMuLocal->normalVector();
1495  ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
1496  GlobalVector Ng = tpMuGlobal->normalVector();
1497  Surface* surf = (Surface*)&refSurface;
1498  const Plane* refPlane = dynamic_cast<Plane*>(surf);
1499  GlobalVector Nlp = refPlane->normalVector();
1500 
1501  if (debug_) {
1502  std::cout << " Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1503  << " alfa " << int(asin(Nl.x()) * 57.29578) << std::endl;
1504  std::cout << " global " << Ng.x() << " " << Ng.y() << " " << Ng.z() << std::endl;
1505  std::cout << " lp " << Nlp.x() << " " << Nlp.y() << " " << Nlp.z() << std::endl;
1506  //std::cout<<" Nlocal Nglobal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
1507  // <<Ng.x()<<" "<<Ng.y()<<" "<<Ng.z()
1508  //<<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
1509  }
1510 
1511  // extrapolation to innermost muon hit
1512 
1513  //extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
1514  if (!refitTrack_)
1515  extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
1516  else {
1517  GlobalTrackerMuonAlignment::trackFitter(itMuon->track(), trackTT, alongMomentum, trackFittedTSOS);
1518  if (trackFittedTSOS.isValid())
1519  extrapolationT = alongSmPr.propagate(trackFittedTSOS, refSurface);
1520  }
1521 
1522  if (!extrapolationT.isValid()) {
1523  if (false & alarm)
1524  std::cout << " !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
1525  //<<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
1526  << std::endl;
1527  continue;
1528  }
1529  tsosMuonIf = 1;
1530  Rt = GlobalVector((extrapolationT.globalPosition()).x(),
1531  (extrapolationT.globalPosition()).y(),
1532  (extrapolationT.globalPosition()).z());
1533 
1534  Pt = extrapolationT.globalMomentum();
1535  // global parameters of muon
1536  GRm = GlobalVector(
1537  (innerMuTSOS.globalPosition()).x(), (innerMuTSOS.globalPosition()).y(), (innerMuTSOS.globalPosition()).z());
1538  GPm = innerMuTSOS.globalMomentum();
1539 
1540  Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
1541  (outerTrackTSOS.globalPosition()).y(),
1542  (outerTrackTSOS.globalPosition()).z());
1543  Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + innerMuTSOS.cartesianError().matrix());
1544  C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
1545  Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
1546  C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
1547 
1548  if (refitMuon_)
1549  GlobalTrackerMuonAlignment::muonFitter(itMuon->outerTrack(), muTT, oppositeToMomentum, muonFittedTSOS);
1550 
1551  } // ------------------------------- end Isolated Muon -- Out - In ---
1552 
1553  if (cosmicMuonMode_) { //------------------------------- Cosmic Muon -----
1554 
1555  if ((distanceOutIn <= distanceInOut) & (distanceOutIn <= distanceInIn) &
1556  (distanceOutIn <= distanceOutOut)) { // ----- Out - In ------
1557  if (debug_)
1558  std::cout << " ----- Out - In -----" << std::endl;
1559 
1560  const Surface& refSurface = innerMuTSOS.surface();
1561  ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
1562  Nl = tpMuGlobal->normalVector();
1563 
1564  // extrapolation to innermost muon hit
1565  //extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
1566  if (!refitTrack_)
1567  extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
1568  else {
1569  GlobalTrackerMuonAlignment::trackFitter(itMuon->track(), trackTT, alongMomentum, trackFittedTSOS);
1570  if (trackFittedTSOS.isValid())
1571  extrapolationT = alongSmPr.propagate(trackFittedTSOS, refSurface);
1572  }
1573 
1574  if (!extrapolationT.isValid()) {
1575  if (false & alarm)
1576  std::cout << " !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
1577  //<<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
1578  << std::endl;
1579  continue;
1580  }
1581  if (debug_)
1582  std::cout << " extrapolationT.isValid " << extrapolationT.isValid() << std::endl;
1583 
1584  tsosMuonIf = 1;
1585  Rt = GlobalVector((extrapolationT.globalPosition()).x(),
1586  (extrapolationT.globalPosition()).y(),
1587  (extrapolationT.globalPosition()).z());
1588 
1589  Pt = extrapolationT.globalMomentum();
1590  // global parameters of muon
1591  GRm = GlobalVector(
1592  (innerMuTSOS.globalPosition()).x(), (innerMuTSOS.globalPosition()).y(), (innerMuTSOS.globalPosition()).z());
1593  GPm = innerMuTSOS.globalMomentum();
1594  Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
1595  (outerTrackTSOS.globalPosition()).y(),
1596  (outerTrackTSOS.globalPosition()).z());
1597  Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + innerMuTSOS.cartesianError().matrix());
1598  C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
1599  Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
1600  C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
1601 
1602  if (refitMuon_)
1603  GlobalTrackerMuonAlignment::muonFitter(itMuon->outerTrack(), muTT, oppositeToMomentum, muonFittedTSOS);
1604 
1605  if (false & debug_) {
1606  //std::cout<<" ->propDir "<<propagator->propagationDirection()<<std::endl;
1607  PropagationDirectionChooser Chooser;
1608  std::cout << " propDirCh " << Chooser.operator()(*outerTrackTSOS.freeState(), refSurface) << " Ch == along "
1609  << (alongMomentum == Chooser.operator()(*outerTrackTSOS.freeState(), refSurface)) << std::endl;
1610  std::cout << " --- Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1611  << " alfa " << int(asin(Nl.x()) * 57.29578) << std::endl;
1612  }
1613  } // enf of ---- Out - In -----
1614 
1615  else if ((distanceInOut <= distanceInIn) & (distanceInOut <= distanceOutIn) &
1616  (distanceInOut <= distanceOutOut)) { // ----- In - Out ------
1617  if (debug_)
1618  std::cout << " ----- In - Out -----" << std::endl;
1619 
1620  const Surface& refSurface = outerMuTSOS.surface();
1621  ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
1622  Nl = tpMuGlobal->normalVector();
1623 
1624  // extrapolation to outermost muon hit
1625  //extrapolationT = oppositeSmPr.propagate(innerTrackTSOS, refSurface);
1626  if (!refitTrack_)
1627  extrapolationT = oppositeSmPr.propagate(innerTrackTSOS, refSurface);
1628  else {
1629  GlobalTrackerMuonAlignment::trackFitter(itMuon->track(), trackTT, oppositeToMomentum, trackFittedTSOS);
1630  if (trackFittedTSOS.isValid())
1631  extrapolationT = oppositeSmPr.propagate(trackFittedTSOS, refSurface);
1632  }
1633 
1634  if (!extrapolationT.isValid()) {
1635  if (false & alarm)
1636  std::cout << " !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
1637  << "\a\a\a\a\a\a\a\a" << extrapolationT.isValid() << std::endl;
1638  continue;
1639  }
1640  if (debug_)
1641  std::cout << " extrapolationT.isValid " << extrapolationT.isValid() << std::endl;
1642 
1643  tsosMuonIf = 2;
1644  Rt = GlobalVector((extrapolationT.globalPosition()).x(),
1645  (extrapolationT.globalPosition()).y(),
1646  (extrapolationT.globalPosition()).z());
1647 
1648  Pt = extrapolationT.globalMomentum();
1649  // global parameters of muon
1650  GRm = GlobalVector(
1651  (outerMuTSOS.globalPosition()).x(), (outerMuTSOS.globalPosition()).y(), (outerMuTSOS.globalPosition()).z());
1652  GPm = outerMuTSOS.globalMomentum();
1653  Rt0 = GlobalVector((innerTrackTSOS.globalPosition()).x(),
1654  (innerTrackTSOS.globalPosition()).y(),
1655  (innerTrackTSOS.globalPosition()).z());
1656  Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + outerMuTSOS.cartesianError().matrix());
1657  C0 = AlgebraicSymMatrix66(innerTrackTSOS.cartesianError().matrix());
1658  Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
1659  C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
1660 
1661  if (refitMuon_)
1662  GlobalTrackerMuonAlignment::muonFitter(itMuon->outerTrack(), muTT, alongMomentum, muonFittedTSOS);
1663 
1664  if (false & debug_) {
1665  //std::cout<<" ->propDir "<<propagator->propagationDirection()<<std::endl;
1666  PropagationDirectionChooser Chooser;
1667  std::cout << " propDirCh " << Chooser.operator()(*innerTrackTSOS.freeState(), refSurface)
1668  << " Ch == oppisite "
1669  << (oppositeToMomentum == Chooser.operator()(*innerTrackTSOS.freeState(), refSurface)) << std::endl;
1670  std::cout << " --- Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1671  << " alfa " << int(asin(Nl.x()) * 57.29578) << std::endl;
1672  }
1673  } // enf of ---- In - Out -----
1674 
1675  else if ((distanceOutOut <= distanceInOut) & (distanceOutOut <= distanceInIn) &
1676  (distanceOutOut <= distanceOutIn)) { // ----- Out - Out ------
1677  if (debug_)
1678  std::cout << " ----- Out - Out -----" << std::endl;
1679 
1680  // reject: momentum of track has opposite direction to muon track
1681  continue;
1682 
1683  const Surface& refSurface = outerMuTSOS.surface();
1684  ConstReferenceCountingPointer<TangentPlane> tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
1685  Nl = tpMuGlobal->normalVector();
1686 
1687  // extrapolation to outermost muon hit
1688  extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
1689 
1690  if (!extrapolationT.isValid()) {
1691  if (alarm)
1692  std::cout << " !!!!!!!!! Catastrophe Out-Out extrapolationT.isValid "
1693  << "\a\a\a\a\a\a\a\a" << extrapolationT.isValid() << std::endl;
1694  continue;
1695  }
1696  if (debug_)
1697  std::cout << " extrapolationT.isValid " << extrapolationT.isValid() << std::endl;
1698 
1699  tsosMuonIf = 2;
1700  Rt = GlobalVector((extrapolationT.globalPosition()).x(),
1701  (extrapolationT.globalPosition()).y(),
1702  (extrapolationT.globalPosition()).z());
1703 
1704  Pt = extrapolationT.globalMomentum();
1705  // global parameters of muon
1706  GRm = GlobalVector(
1707  (outerMuTSOS.globalPosition()).x(), (outerMuTSOS.globalPosition()).y(), (outerMuTSOS.globalPosition()).z());
1708  GPm = outerMuTSOS.globalMomentum();
1709  Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
1710  (outerTrackTSOS.globalPosition()).y(),
1711  (outerTrackTSOS.globalPosition()).z());
1712  Cm = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix() + outerMuTSOS.cartesianError().matrix());
1713  C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
1714  Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
1715  C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
1716 
1717  if (debug_) {
1718  PropagationDirectionChooser Chooser;
1719  std::cout << " propDirCh " << Chooser.operator()(*outerTrackTSOS.freeState(), refSurface) << " Ch == along "
1720  << (alongMomentum == Chooser.operator()(*outerTrackTSOS.freeState(), refSurface)) << std::endl;
1721  std::cout << " --- Nlocal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1722  << " alfa " << int(asin(Nl.x()) * 57.29578) << std::endl;
1723  std::cout << " Nornal " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " "
1724  << " alfa " << int(asin(Nl.x()) * 57.29578) << std::endl;
1725  }
1726  } // enf of ---- Out - Out -----
1727  else {
1728  if (alarm)
1729  std::cout << " ------- !!!!!!!!!!!!!!! No proper Out - In \a\a\a\a\a\a\a" << std::endl;
1730  continue;
1731  }
1732 
1733  } // ------------------------------- end Cosmic Muon -----
1734 
1735  if (tsosMuonIf == 0) {
1736  if (info) {
1737  std::cout << "No tsosMuon !!!!!!" << std::endl;
1738  continue;
1739  }
1740  }
1741  TrajectoryStateOnSurface tsosMuon;
1742  if (tsosMuonIf == 1)
1743  tsosMuon = muTT.innermostMeasurementState();
1744  else
1745  tsosMuon = muTT.outermostMeasurementState();
1746 
1747  //GlobalTrackerMuonAlignment::misalignMuon(GRm, GPm, Nl, Rt, Rm, Pm);
1748  AlgebraicVector4 LPRm; // muon local (dx/dz, dy/dz, x, y)
1749  GlobalTrackerMuonAlignment::misalignMuonL(GRm, GPm, Nl, Rt, Rm, Pm, LPRm, extrapolationT, tsosMuon);
1750 
1751  if (refitTrack_) {
1752  if (!trackFittedTSOS.isValid()) {
1753  if (info)
1754  std::cout << " ================= trackFittedTSOS notValid !!!!!!!! " << std::endl;
1755  continue;
1756  }
1757  if (debug_)
1758  this->debugTrajectorySOS(" trackFittedTSOS ", trackFittedTSOS);
1759  }
1760 
1761  if (refitMuon_) {
1762  if (!muonFittedTSOS.isValid()) {
1763  if (info)
1764  std::cout << " ================= muonFittedTSOS notValid !!!!!!!! " << std::endl;
1765  continue;
1766  }
1767  if (debug_)
1768  this->debugTrajectorySOS(" muonFittedTSOS ", muonFittedTSOS);
1769  Rm = GlobalVector((muonFittedTSOS.globalPosition()).x(),
1770  (muonFittedTSOS.globalPosition()).y(),
1771  (muonFittedTSOS.globalPosition()).z());
1772  Pm = muonFittedTSOS.globalMomentum();
1773  LPRm = AlgebraicVector4(muonFittedTSOS.localParameters().vector()(1),
1774  muonFittedTSOS.localParameters().vector()(2),
1775  muonFittedTSOS.localParameters().vector()(3),
1776  muonFittedTSOS.localParameters().vector()(4));
1777  }
1778  GlobalVector resR = Rm - Rt;
1779  GlobalVector resP0 = Pm - Pt;
1780  GlobalVector resP = Pm / Pm.mag() - Pt / Pt.mag();
1781  float RelMomResidual = (Pm.mag() - Pt.mag()) / (Pt.mag() + 1.e-6);
1782  ;
1783 
1784  AlgebraicVector6 Vm;
1785  Vm(0) = resR.x();
1786  Vm(1) = resR.y();
1787  Vm(2) = resR.z();
1788  Vm(3) = resP0.x();
1789  Vm(4) = resP0.y();
1790  Vm(5) = resP0.z();
1791  float Rmuon = Rm.perp();
1792  float Zmuon = Rm.z();
1793  float alfa_x = atan2(Nl.x(), Nl.y()) * 57.29578;
1794 
1795  if (debug_) {
1796  std::cout << " Nx Ny Nz alfa_x " << Nl.x() << " " << Nl.y() << " " << Nl.z() << " " << alfa_x << std::endl;
1797  //std::cout<<" Rm "<<Rm<<std::endl<<" Rt "<<Rt<<std::endl<<" resR "<<resR<<std::endl
1798  // <<" resP "<<resP<<" dp/p "<<RelMomResidual<<std::endl;
1799  }
1800 
1801  double chi_d = 0;
1802  for (int i = 0; i <= 5; i++)
1803  chi_d += Vm(i) * Vm(i) / Cm(i, i);
1804 
1805  AlgebraicVector5 Vml(tsosMuon.localParameters().vector() - extrapolationT.localParameters().vector());
1806  AlgebraicSymMatrix55 m(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
1807  AlgebraicSymMatrix55 Cml(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
1808  bool ierrLoc = !m.Invert();
1809  if (ierrLoc && debug_ && info) {
1810  std::cout << " ==== Error inverting Local covariance matrix ==== " << std::endl;
1811  continue;
1812  }
1813  double chi_Loc = ROOT::Math::Similarity(Vml, m);
1814  if (debug_)
1815  std::cout << " chi_Loc px/pz/err " << chi_Loc << " " << Vml(1) / sqrt(Cml(1, 1)) << std::endl;
1816 
1817  if (Pt.mag() < 15.)
1818  continue;
1819  if (Pm.mag() < 5.)
1820  continue;
1821 
1822  //if(Pt.mag() < 30.) continue; // momenum cut < 30GeV
1823  //if(Pt.mag() < 60.) continue; // momenum cut < 30GeV
1824  //if(Pt.mag() > 50.) continue; // momenum cut > 50GeV
1825  //if(Pt.mag() > 100.) continue; // momenum cut > 100GeV
1826  //if(trackTT.charge() < 0) continue; // select positive charge
1827  //if(trackTT.charge() > 0) continue; // select negative charge
1828 
1829  //if(fabs(resR.x()) > 5.) continue; // strong cut X
1830  //if(fabs(resR.y()) > 5.) continue; // Y
1831  //if(fabs(resR.z()) > 5.) continue; // Z
1832  //if(fabs(resR.mag()) > 7.5) continue; // dR
1833 
1834  //if(fabs(RelMomResidual) > 0.5) continue;
1835  if (fabs(resR.x()) > 20.)
1836  continue;
1837  if (fabs(resR.y()) > 20.)
1838  continue;
1839  if (fabs(resR.z()) > 20.)
1840  continue;
1841  if (fabs(resR.mag()) > 30.)
1842  continue;
1843  if (fabs(resP.x()) > 0.06)
1844  continue;
1845  if (fabs(resP.y()) > 0.06)
1846  continue;
1847  if (fabs(resP.z()) > 0.06)
1848  continue;
1849  if (chi_d > 40.)
1850  continue;
1851 
1852  // select Barrel
1853  //if(Rmuon < 400. || Rmuon > 450.) continue;
1854  //if(Zmuon < -600. || Zmuon > 600.) continue;
1855  //if(fabs(Nl.z()) > 0.95) continue;
1856  //MuSelect = " Barrel";
1857  // EndCap1
1858  //if(Rmuon < 120. || Rmuon > 450.) continue;
1859  //if(Zmuon < -720.) continue;
1860  //if(Zmuon > -580.) continue;
1861  //if(fabs(Nl.z()) < 0.95) continue;
1862  //MuSelect = " EndCap1";
1863  // EndCap2
1864  //if(Rmuon < 120. || Rmuon > 450.) continue;
1865  //if(Zmuon > 720.) continue;
1866  //if(Zmuon < 580.) continue;
1867  //if(fabs(Nl.z()) < 0.95) continue;
1868  //MuSelect = " EndCap2";
1869  // select All
1870  if (Rmuon < 120. || Rmuon > 450.)
1871  continue;
1872  if (Zmuon < -720. || Zmuon > 720.)
1873  continue;
1874  MuSelect = " Barrel+EndCaps";
1875 
1876  if (debug_)
1877  std::cout << " .............. passed all cuts" << std::endl;
1878 
1879  N_track++;
1880  // gradient and Hessian for each track
1881 
1882  GlobalTrackerMuonAlignment::gradientGlobalAlg(Rt, Pt, Rm, Nl, Cm);
1883  GlobalTrackerMuonAlignment::gradientGlobal(Rt, Pt, Rm, Pm, Nl, Cm);
1884 
1885  CLHEP::HepSymMatrix covLoc(4, 0);
1886  for (int i = 1; i <= 4; i++)
1887  for (int j = 1; j <= i; j++) {
1888  covLoc(i, j) = (tsosMuon.localError().matrix() + extrapolationT.localError().matrix())(i, j);
1889  //if(i != j) Cov(i,j) = 0.;
1890  }
1891 
1892  const Surface& refSurface = tsosMuon.surface();
1893  CLHEP::HepMatrix rotLoc(3, 3, 0);
1894  rotLoc(1, 1) = refSurface.rotation().xx();
1895  rotLoc(1, 2) = refSurface.rotation().xy();
1896  rotLoc(1, 3) = refSurface.rotation().xz();
1897 
1898  rotLoc(2, 1) = refSurface.rotation().yx();
1899  rotLoc(2, 2) = refSurface.rotation().yy();
1900  rotLoc(2, 3) = refSurface.rotation().yz();
1901 
1902  rotLoc(3, 1) = refSurface.rotation().zx();
1903  rotLoc(3, 2) = refSurface.rotation().zy();
1904  rotLoc(3, 3) = refSurface.rotation().zz();
1905 
1906  CLHEP::HepVector posLoc(3);
1907  posLoc(1) = refSurface.position().x();
1908  posLoc(2) = refSurface.position().y();
1909  posLoc(3) = refSurface.position().z();
1910 
1911  GlobalTrackerMuonAlignment::gradientLocal(Rt, Pt, Rm, Pm, Nl, covLoc, rotLoc, posLoc, LPRm);
1912 
1913  if (debug_) {
1914  std::cout << " Norm " << Nl << std::endl;
1915  std::cout << " posLoc " << posLoc.T() << std::endl;
1916  std::cout << " rotLoc " << rotLoc << std::endl;
1917  }
1918 
1919  // ----------------------------------------------------- fill histogram
1920  histo->Fill(itMuon->track()->pt());
1921 
1922  //histo2->Fill(itMuon->track()->outerP());
1923  histo2->Fill(Pt.mag());
1924  histo3->Fill((PI / 2. - itMuon->track()->outerTheta()));
1925  histo4->Fill(itMuon->track()->phi());
1926  histo5->Fill(Rmuon);
1927  histo6->Fill(Zmuon);
1928  histo7->Fill(RelMomResidual);
1929  //histo8->Fill(chi);
1930  histo8->Fill(chi_d);
1931 
1932  histo101->Fill(Zmuon, Rmuon);
1933  histo101->Fill(Rt0.z(), Rt0.perp());
1934  histo102->Fill(Rt0.x(), Rt0.y());
1935  histo102->Fill(Rm.x(), Rm.y());
1936 
1937  histo11->Fill(resR.mag());
1938  if (fabs(Nl.x()) < 0.98)
1939  histo12->Fill(resR.x());
1940  if (fabs(Nl.y()) < 0.98)
1941  histo13->Fill(resR.y());
1942  if (fabs(Nl.z()) < 0.98)
1943  histo14->Fill(resR.z());
1944  histo15->Fill(resP.x());
1945  histo16->Fill(resP.y());
1946  histo17->Fill(resP.z());
1947 
1948  if ((fabs(Nl.x()) < 0.98) && (fabs(Nl.y()) < 0.98) && (fabs(Nl.z()) < 0.98)) {
1949  histo18->Fill(sqrt(C0(0, 0)));
1950  histo19->Fill(sqrt(C1(0, 0)));
1951  histo20->Fill(sqrt(C1(0, 0) + Ce(0, 0)));
1952  }
1953  if (fabs(Nl.x()) < 0.98)
1954  histo21->Fill(Vm(0) / sqrt(Cm(0, 0)));
1955  if (fabs(Nl.y()) < 0.98)
1956  histo22->Fill(Vm(1) / sqrt(Cm(1, 1)));
1957  if (fabs(Nl.z()) < 0.98)
1958  histo23->Fill(Vm(2) / sqrt(Cm(2, 2)));
1959  histo24->Fill(Vm(3) / sqrt(C1(3, 3) + Ce(3, 3)));
1960  histo25->Fill(Vm(4) / sqrt(C1(4, 4) + Ce(4, 4)));
1961  histo26->Fill(Vm(5) / sqrt(C1(5, 5) + Ce(5, 5)));
1962  histo27->Fill(Nl.x());
1963  histo28->Fill(Nl.y());
1964  histo29->Fill(lenghtTrack);
1965  histo30->Fill(lenghtMuon);
1966  histo31->Fill(chi_Loc);
1967  histo32->Fill(Vml(1) / sqrt(Cml(1, 1)));
1968  histo33->Fill(Vml(2) / sqrt(Cml(2, 2)));
1969  histo34->Fill(Vml(3) / sqrt(Cml(3, 3)));
1970  histo35->Fill(Vml(4) / sqrt(Cml(4, 4)));
1971 
1972  if (debug_) { //--------------------------------- debug print ----------
1973 
1974  std::cout << " diag 0[ " << C0(0, 0) << " " << C0(1, 1) << " " << C0(2, 2) << " " << C0(3, 3) << " " << C0(4, 4)
1975  << " " << C0(5, 5) << " ]" << std::endl;
1976  std::cout << " diag e[ " << Ce(0, 0) << " " << Ce(1, 1) << " " << Ce(2, 2) << " " << Ce(3, 3) << " " << Ce(4, 4)
1977  << " " << Ce(5, 5) << " ]" << std::endl;
1978  std::cout << " diag 1[ " << C1(0, 0) << " " << C1(1, 1) << " " << C1(2, 2) << " " << C1(3, 3) << " " << C1(4, 4)
1979  << " " << C1(5, 5) << " ]" << std::endl;
1980  std::cout << " Rm " << Rm.x() << " " << Rm.y() << " " << Rm.z() << " Pm " << Pm.x() << " " << Pm.y() << " "
1981  << Pm.z() << std::endl;
1982  std::cout << " Rt " << Rt.x() << " " << Rt.y() << " " << Rt.z() << " Pt " << Pt.x() << " " << Pt.y() << " "
1983  << Pt.z() << std::endl;
1984  std::cout << " Nl*(Rm-Rt) " << Nl.dot(Rm - Rt) << std::endl;
1985  std::cout << " resR " << resR.x() << " " << resR.y() << " " << resR.z() << " resP " << resP.x() << " " << resP.y()
1986  << " " << resP.z() << std::endl;
1987  std::cout << " Rm-t " << (Rm - Rt).x() << " " << (Rm - Rt).y() << " " << (Rm - Rt).z() << " Pm-t "
1988  << (Pm - Pt).x() << " " << (Pm - Pt).y() << " " << (Pm - Pt).z() << std::endl;
1989  std::cout << " Vm " << Vm << std::endl;
1990  std::cout << " +- " << sqrt(Cm(0, 0)) << " " << sqrt(Cm(1, 1)) << " " << sqrt(Cm(2, 2)) << " " << sqrt(Cm(3, 3))
1991  << " " << sqrt(Cm(4, 4)) << " " << sqrt(Cm(5, 5)) << std::endl;
1992  std::cout << " Pmuon Ptrack dP/Ptrack " << itMuon->outerTrack()->p() << " " << itMuon->track()->outerP() << " "
1993  << (itMuon->outerTrack()->p() - itMuon->track()->outerP()) / itMuon->track()->outerP() << std::endl;
1994  std::cout << " cov matrix " << std::endl;
1995  std::cout << Cm << std::endl;
1996  std::cout << " diag [ " << Cm(0, 0) << " " << Cm(1, 1) << " " << Cm(2, 2) << " " << Cm(3, 3) << " " << Cm(4, 4)
1997  << " " << Cm(5, 5) << " ]" << std::endl;
1998 
1999  AlgebraicSymMatrix66 Ro;
2000  double Diag[6];
2001  for (int i = 0; i <= 5; i++)
2002  Diag[i] = sqrt(Cm(i, i));
2003  for (int i = 0; i <= 5; i++)
2004  for (int j = 0; j <= 5; j++)
2005  Ro(i, j) = Cm(i, j) / Diag[i] / Diag[j];
2006  std::cout << " correlation matrix " << std::endl;
2007  std::cout << Ro << std::endl;
2008 
2009  AlgebraicSymMatrix66 CmI;
2010  for (int i = 0; i <= 5; i++)
2011  for (int j = 0; j <= 5; j++)
2012  CmI(i, j) = Cm(i, j);
2013 
2014  bool ierr = !CmI.Invert();
2015  if (ierr) {
2016  if (alarm || debug_)
2017  std::cout << " Error inverse covariance matrix !!!!!!!!!!!" << std::endl;
2018  continue;
2019  }
2020  std::cout << " inverse cov matrix " << std::endl;
2021  std::cout << Cm << std::endl;
2022 
2023  double chi = ROOT::Math::Similarity(Vm, CmI);
2024  std::cout << " chi chi_d " << chi << " " << chi_d << std::endl;
2025  } // end of debug_ printout --------------------------------------------
2026 
2027  } // end loop on selected muons, i.e. Jim's globalMuon
2028 
2029 } //end of analyzeTrackTrajectory

References alongMomentum, TrajectoryStateOnSurface::cartesianError(), edm::ESWatcher< T >::check(), Chi2MeasurementEstimator_cfi::Chi2MeasurementEstimator, collectionCosmic, collectionIsolated, cosmicMuonMode_, gather_cfg::cout, debug_, debugTrajectorySOS(), defineFitter, Vector3DBase< T, FrameTag >::dot(), MillePedeFileConverter_cfg::e, DetId::Ecal, geometryDiff::epsilon, TrajectoryStateOnSurface::freeState(), edm::EventSetup::get(), get, TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), globalPositionRcd_, gmuonTags_, gradientGlobal(), gradientGlobalAlg(), gradientLocal(), DetId::Hcal, histo, histo101, histo102, histo11, histo12, histo13, histo14, histo15, histo16, histo17, histo18, histo19, histo2, histo20, histo21, histo22, histo23, histo24, histo25, histo26, histo27, histo28, histo29, histo3, histo30, histo31, histo32, histo33, histo34, histo35, histo4, histo5, histo6, histo7, histo8, mps_fire::i, iEvent, info(), createfilelist::int, isolatedMuonMode_, edm::ESHandleBase::isValid(), TrajectoryStateOnSurface::isValid(), iteratorEcalRcd, iteratorHcalRcd, iteratorMuonRcd, iteratorTrackerRcd, dqmiolumiharvest::j, TrajectoryStateOnSurface::localError(), TrajectoryStateOnSurface::localParameters(), visualization-live-secondInstance_cfg::m, Alignments::m_align, PV3DBase< T, PVType, FrameType >::mag(), mag(), HLT_FULL_cff::magneticField, magneticField_, CartesianTrajectoryError::matrix(), LocalTrajectoryError::matrix(), misalignMuonL(), DetId::Muon, muonFitter(), PDWG_BPHSkim_cff::muons, muonTags_, MuRHBuilder, MuSelect, N_event, N_track, Plane::normalVector(), oppositeToMomentum, PV3DBase< T, PVType, FrameType >::perp(), GloballyPositioned< T >::position(), SmartPropagator::propagate(), TrackCandidateProducer_cfi::propagator, defaultRKPropagator::Product::propagator, propagator_, DetId::rawId(), refitMuon_, refitTrack_, GloballyPositioned< T >::rotation(), SmartPropagator_cfi::SmartPropagator, smuonTags_, mathSSE::sqrt(), SteppingHelixPropagator_cfi::SteppingHelixPropagator, TrajectoryStateOnSurface::surface(), Surface::tangentPlane(), theFitter, theFitterOp, theSmoother, theSmootherOp, theTrackingGeometry, DetId::Tracker, trackFitter(), trackingGeometry_, PDWG_EXOHSCP_cff::tracks, trackTags_, TTRHBuilder, LocalTrajectoryParameters::vector(), watchGlobalPositionRcd_, watchMagneticFieldRecord_, watchTrackingGeometry_, x, PV3DBase< T, PVType, FrameType >::x(), TkRotation< T >::xx(), TkRotation< T >::xy(), TkRotation< T >::xz(), y, PV3DBase< T, PVType, FrameType >::y(), TkRotation< T >::yx(), TkRotation< T >::yy(), TkRotation< T >::yz(), z, PV3DBase< T, PVType, FrameType >::z(), TkRotation< T >::zx(), TkRotation< T >::zy(), and TkRotation< T >::zz().

Referenced by analyze().

◆ beginJob()

void GlobalTrackerMuonAlignment::beginJob ( void  )
overrideprivatevirtual

Reimplemented from edm::EDAnalyzer.

Definition at line 318 of file GlobalTrackerMuonAlignment.cc.

318  {
319  N_event = 0;
320  N_track = 0;
321 
322  if (cosmicMuonMode_ == true && isolatedMuonMode_ == false) {
323  collectionCosmic = true;
324  collectionIsolated = false;
325  } else if (cosmicMuonMode_ == false && isolatedMuonMode_ == true) {
326  collectionCosmic = false;
327  collectionIsolated = true;
328  } else {
329  collectionCosmic = false;
330  collectionIsolated = false;
331  }
332 
333  for (int i = 0; i <= 2; i++) {
334  Gfr(i) = 0.;
335  for (int j = 0; j <= 2; j++) {
336  Inf(i, j) = 0.;
337  }
338  }
339 
340  Grad = CLHEP::HepVector(6, 0);
341  Hess = CLHEP::HepSymMatrix(6, 0);
342 
343  GradL = CLHEP::HepVector(6, 0);
344  HessL = CLHEP::HepSymMatrix(6, 0);
345 
346  // histograms
347  TDirectory* dirsave = gDirectory;
348 
349  file = new TFile(rootOutFile_.c_str(), "recreate");
350  const bool oldAddDir = TH1::AddDirectoryStatus();
351 
352  TH1::AddDirectory(true);
353 
354  this->bookHist();
355 
356  TH1::AddDirectory(oldAddDir);
357  dirsave->cd();
358 }

References bookHist(), collectionCosmic, collectionIsolated, cosmicMuonMode_, file, Gfr, Grad, GradL, Hess, HessL, mps_fire::i, Inf, isolatedMuonMode_, dqmiolumiharvest::j, N_event, N_track, and rootOutFile_.

◆ bookHist()

void GlobalTrackerMuonAlignment::bookHist ( )

Definition at line 471 of file GlobalTrackerMuonAlignment.cc.

471  {
472  double PI = 3.1415927;
473  histo = new TH1F("Pt", "pt", 1000, 0, 100);
474  histo2 = new TH1F("P", "P [GeV/c]", 400, 0., 400.);
475  histo2->GetXaxis()->SetTitle("momentum [GeV/c]");
476  histo3 = new TH1F("outerLambda", "#lambda outer", 100, -PI / 2., PI / 2.);
477  histo3->GetXaxis()->SetTitle("#lambda outer");
478  histo4 = new TH1F("phi", "#phi [rad]", 100, -PI, PI);
479  histo4->GetXaxis()->SetTitle("#phi [rad]");
480  histo5 = new TH1F("Rmuon", "inner muon hit R [cm]", 100, 0., 800.);
481  histo5->GetXaxis()->SetTitle("R of muon [cm]");
482  histo6 = new TH1F("Zmuon", "inner muon hit Z[cm]", 100, -1000., 1000.);
483  histo6->GetXaxis()->SetTitle("Z of muon [cm]");
484  histo7 = new TH1F("(Pm-Pt)/Pt", " (Pmuon-Ptrack)/Ptrack", 100, -2., 2.);
485  histo7->GetXaxis()->SetTitle("(Pmuon-Ptrack)/Ptrack");
486  histo8 = new TH1F("chi muon-track", "#chi^{2}(muon-track)", 1000, 0., 1000.);
487  histo8->GetXaxis()->SetTitle("#chi^{2} of muon w.r.t. propagated track");
488  histo11 = new TH1F("distance muon-track", "distance muon w.r.t track [cm]", 100, 0., 30.);
489  histo11->GetXaxis()->SetTitle("distance of muon w.r.t. track [cm]");
490  histo12 = new TH1F("Xmuon-Xtrack", "Xmuon-Xtrack [cm]", 200, -20., 20.);
491  histo12->GetXaxis()->SetTitle("Xmuon - Xtrack [cm]");
492  histo13 = new TH1F("Ymuon-Ytrack", "Ymuon-Ytrack [cm]", 200, -20., 20.);
493  histo13->GetXaxis()->SetTitle("Ymuon - Ytrack [cm]");
494  histo14 = new TH1F("Zmuon-Ztrack", "Zmuon-Ztrack [cm]", 200, -20., 20.);
495  histo14->GetXaxis()->SetTitle("Zmuon-Ztrack [cm]");
496  histo15 = new TH1F("NXmuon-NXtrack", "NXmuon-NXtrack [rad]", 200, -.1, .1);
497  histo15->GetXaxis()->SetTitle("N_{X}(muon)-N_{X}(track) [rad]");
498  histo16 = new TH1F("NYmuon-NYtrack", "NYmuon-NYtrack [rad]", 200, -.1, .1);
499  histo16->GetXaxis()->SetTitle("N_{Y}(muon)-N_{Y}(track) [rad]");
500  histo17 = new TH1F("NZmuon-NZtrack", "NZmuon-NZtrack [rad]", 200, -.1, .1);
501  histo17->GetXaxis()->SetTitle("N_{Z}(muon)-N_{Z}(track) [rad]");
502  histo18 = new TH1F("expected error of Xinner", "outer hit of inner tracker", 100, 0, .01);
503  histo18->GetXaxis()->SetTitle("expected error of Xinner [cm]");
504  histo19 = new TH1F("expected error of Xmuon", "inner hit of muon", 100, 0, .1);
505  histo19->GetXaxis()->SetTitle("expected error of Xmuon [cm]");
506  histo20 = new TH1F("expected error of Xmuon-Xtrack", "muon w.r.t. propagated track", 100, 0., 10.);
507  histo20->GetXaxis()->SetTitle("expected error of Xmuon-Xtrack [cm]");
508  histo21 = new TH1F("pull of Xmuon-Xtrack", "pull of Xmuon-Xtrack", 100, -10., 10.);
509  histo21->GetXaxis()->SetTitle("(Xmuon-Xtrack)/expected error");
510  histo22 = new TH1F("pull of Ymuon-Ytrack", "pull of Ymuon-Ytrack", 100, -10., 10.);
511  histo22->GetXaxis()->SetTitle("(Ymuon-Ytrack)/expected error");
512  histo23 = new TH1F("pull of Zmuon-Ztrack", "pull of Zmuon-Ztrack", 100, -10., 10.);
513  histo23->GetXaxis()->SetTitle("(Zmuon-Ztrack)/expected error");
514  histo24 = new TH1F("pull of PXmuon-PXtrack", "pull of PXmuon-PXtrack", 100, -10., 10.);
515  histo24->GetXaxis()->SetTitle("(P_{X}(muon)-P_{X}(track))/expected error");
516  histo25 = new TH1F("pull of PYmuon-PYtrack", "pull of PYmuon-PYtrack", 100, -10., 10.);
517  histo25->GetXaxis()->SetTitle("(P_{Y}(muon)-P_{Y}(track))/expected error");
518  histo26 = new TH1F("pull of PZmuon-PZtrack", "pull of PZmuon-PZtrack", 100, -10., 10.);
519  histo26->GetXaxis()->SetTitle("(P_{Z}(muon)-P_{Z}(track))/expected error");
520  histo27 = new TH1F("N_x", "Nx of tangent plane", 120, -1.1, 1.1);
521  histo27->GetXaxis()->SetTitle("normal vector projection N_{X}");
522  histo28 = new TH1F("N_y", "Ny of tangent plane", 120, -1.1, 1.1);
523  histo28->GetXaxis()->SetTitle("normal vector projection N_{Y}");
524  histo29 = new TH1F("lenght of track", "lenght of track", 200, 0., 400);
525  histo29->GetXaxis()->SetTitle("lenght of track [cm]");
526  histo30 = new TH1F("lenght of muon", "lenght of muon", 200, 0., 800);
527  histo30->GetXaxis()->SetTitle("lenght of muon [cm]");
528 
529  histo31 = new TH1F("local chi muon-track", "#local chi^{2}(muon-track)", 1000, 0., 1000.);
530  histo31->GetXaxis()->SetTitle("#local chi^{2} of muon w.r.t. propagated track");
531  histo32 = new TH1F("pull of Px/Pz local", "pull of Px/Pz local", 100, -10., 10.);
532  histo32->GetXaxis()->SetTitle("local (Px/Pz(muon) - Px/Pz(track))/expected error");
533  histo33 = new TH1F("pull of Py/Pz local", "pull of Py/Pz local", 100, -10., 10.);
534  histo33->GetXaxis()->SetTitle("local (Py/Pz(muon) - Py/Pz(track))/expected error");
535  histo34 = new TH1F("pull of X local", "pull of X local", 100, -10., 10.);
536  histo34->GetXaxis()->SetTitle("local (Xmuon - Xtrack)/expected error");
537  histo35 = new TH1F("pull of Y local", "pull of Y local", 100, -10., 10.);
538  histo35->GetXaxis()->SetTitle("local (Ymuon - Ytrack)/expected error");
539 
540  histo101 = new TH2F("Rtr/mu vs Ztr/mu", "hit of track/muon", 100, -800., 800., 100, 0., 600.);
541  histo101->GetXaxis()->SetTitle("Z of track/muon [cm]");
542  histo101->GetYaxis()->SetTitle("R of track/muon [cm]");
543  histo102 = new TH2F("Ytr/mu vs Xtr/mu", "hit of track/muon", 100, -600., 600., 100, -600., 600.);
544  histo102->GetXaxis()->SetTitle("X of track/muon [cm]");
545  histo102->GetYaxis()->SetTitle("Y of track/muon [cm]");
546 }

References histo, histo101, histo102, histo11, histo12, histo13, histo14, histo15, histo16, histo17, histo18, histo19, histo2, histo20, histo21, histo22, histo23, histo24, histo25, histo26, histo27, histo28, histo29, histo3, histo30, histo31, histo32, histo33, histo34, histo35, histo4, histo5, histo6, histo7, and histo8.

Referenced by beginJob().

◆ CLHEP_dot()

double GlobalTrackerMuonAlignment::CLHEP_dot ( const CLHEP::HepVector &  a,
const CLHEP::HepVector &  b 
)
inline

Definition at line 148 of file GlobalTrackerMuonAlignment.cc.

148  {
149  return a(1) * b(1) + a(2) * b(2) + a(3) * b(3);
150  }

References a, and b.

Referenced by gradientGlobal(), gradientLocal(), misalignMuon(), and misalignMuonL().

◆ debugTrackHit() [1/2]

void GlobalTrackerMuonAlignment::debugTrackHit ( const std::string  title,
reco::TrackRef  alongTr 
)

Definition at line 3231 of file GlobalTrackerMuonAlignment.cc.

3231  {
3232  std::cout << " ------- " << title << " --------" << std::endl;
3233  int nHit = 1;
3234  for (auto const& hit : alongTr->recHits()) {
3235  std::cout << " Hit " << nHit++ << " DetId " << hit->geographicalId().det();
3236  if (hit->geographicalId().det() == DetId::Tracker)
3237  std::cout << " Tracker ";
3238  else if (hit->geographicalId().det() == DetId::Muon)
3239  std::cout << " Muon ";
3240  else
3241  std::cout << " Unknown ";
3242  if (!hit->isValid())
3243  std::cout << " not valid " << std::endl;
3244  else
3245  std::cout << std::endl;
3246  }
3247 }

References gather_cfg::cout, DetId::Muon, runGCPTkAlMap::title, and DetId::Tracker.

Referenced by muonFitter(), and trackFitter().

◆ debugTrackHit() [2/2]

void GlobalTrackerMuonAlignment::debugTrackHit ( const std::string  title,
reco::TransientTrack &  alongTr 
)

Definition at line 3212 of file GlobalTrackerMuonAlignment.cc.

3212  {
3213  std::cout << " ------- " << title << " --------" << std::endl;
3214  int nHit = 1;
3215  for (trackingRecHit_iterator i = alongTr.recHitsBegin(); i != alongTr.recHitsEnd(); i++) {
3216  std::cout << " Hit " << nHit++ << " DetId " << (*i)->geographicalId().det();
3217  if ((*i)->geographicalId().det() == DetId::Tracker)
3218  std::cout << " Tracker ";
3219  else if ((*i)->geographicalId().det() == DetId::Muon)
3220  std::cout << " Muon ";
3221  else
3222  std::cout << " Unknown ";
3223  if (!(*i)->isValid())
3224  std::cout << " not valid " << std::endl;
3225  else
3226  std::cout << std::endl;
3227  }
3228 }

References gather_cfg::cout, mps_fire::i, DetId::Muon, reco::TransientTrack::recHitsBegin(), reco::TransientTrack::recHitsEnd(), runGCPTkAlMap::title, and DetId::Tracker.

◆ debugTrajectory()

void GlobalTrackerMuonAlignment::debugTrajectory ( const std::string  title,
Trajectory &  traj 
)

Definition at line 3302 of file GlobalTrackerMuonAlignment.cc.

3302  {
3303  std::cout << "\n"
3304  << " ...... " << title << " ...... " << std::endl;
3305  if (!traj.isValid()) {
3306  std::cout << " Not valid !!!!!!!! " << std::endl;
3307  return;
3308  }
3309  std::cout << " chi2/Nhit " << traj.chiSquared() << " / " << traj.foundHits();
3310  if (traj.direction() == alongMomentum)
3311  std::cout << " alongMomentum >>>>" << std::endl;
3312  else
3313  std::cout << " oppositeToMomentum <<<<" << std::endl;
3314  this->debugTrajectorySOSv(" firstMeasurementTSOS ", traj.firstMeasurement().updatedState());
3315  //this->debugTrajectorySOSv(" firstMeasurementTSOS ",traj.firstMeasurement().predictedState());
3316  this->debugTrajectorySOSv(" lastMeasurementTSOS ", traj.lastMeasurement().updatedState());
3317  //this->debugTrajectorySOSv(" geom InnermostState", traj.geometricalInnermostState());
3318  std::cout << " . . . . . . . . . . . . . . . . . . . . . . . . . . . . \n" << std::endl;
3319 }

References alongMomentum, Trajectory::chiSquared(), gather_cfg::cout, debugTrajectorySOSv(), Trajectory::direction(), Trajectory::firstMeasurement(), Trajectory::foundHits(), Trajectory::isValid(), Trajectory::lastMeasurement(), runGCPTkAlMap::title, and TrajectoryMeasurement::updatedState().

Referenced by muonFitter(), and trackFitter().

◆ debugTrajectorySOS()

void GlobalTrackerMuonAlignment::debugTrajectorySOS ( const std::string  title,
TrajectoryStateOnSurface &  trajSOS 
)

Definition at line 3250 of file GlobalTrackerMuonAlignment.cc.

3250  {
3251  std::cout << " --- " << title << " --- " << std::endl;
3252  if (!trajSOS.isValid()) {
3253  std::cout << " Not valid !!!! " << std::endl;
3254  return;
3255  }
3256  std::cout << " R |p| " << trajSOS.globalPosition().perp() << " " << trajSOS.globalMomentum().mag() << " charge "
3257  << trajSOS.charge() << std::endl;
3258  std::cout << " x p " << trajSOS.globalParameters().vector()(0) << " " << trajSOS.globalParameters().vector()(1)
3259  << " " << trajSOS.globalParameters().vector()(2) << " " << trajSOS.globalParameters().vector()(3) << " "
3260  << trajSOS.globalParameters().vector()(4) << " " << trajSOS.globalParameters().vector()(5) << std::endl;
3261  std::cout << " +/- " << sqrt(trajSOS.cartesianError().matrix()(0, 0)) << " "
3262  << sqrt(trajSOS.cartesianError().matrix()(1, 1)) << " " << sqrt(trajSOS.cartesianError().matrix()(2, 2))
3263  << " " << sqrt(trajSOS.cartesianError().matrix()(3, 3)) << " "
3264  << sqrt(trajSOS.cartesianError().matrix()(4, 4)) << " " << sqrt(trajSOS.cartesianError().matrix()(5, 5))
3265  << std::endl;
3266  std::cout << " q/p dxy/dz xy " << trajSOS.localParameters().vector()(0) << " "
3267  << trajSOS.localParameters().vector()(1) << " " << trajSOS.localParameters().vector()(2) << " "
3268  << trajSOS.localParameters().vector()(3) << " " << trajSOS.localParameters().vector()(4) << std::endl;
3269  std::cout << " +/- error " << sqrt(trajSOS.localError().matrix()(0, 0)) << " "
3270  << sqrt(trajSOS.localError().matrix()(1, 1)) << " " << sqrt(trajSOS.localError().matrix()(2, 2)) << " "
3271  << sqrt(trajSOS.localError().matrix()(3, 3)) << " " << sqrt(trajSOS.localError().matrix()(4, 4)) << " "
3272  << std::endl;
3273 }

References TrajectoryStateOnSurface::cartesianError(), TrajectoryStateOnSurface::charge(), gather_cfg::cout, TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalParameters(), TrajectoryStateOnSurface::globalPosition(), TrajectoryStateOnSurface::isValid(), TrajectoryStateOnSurface::localError(), TrajectoryStateOnSurface::localParameters(), PV3DBase< T, PVType, FrameType >::mag(), CartesianTrajectoryError::matrix(), LocalTrajectoryError::matrix(), PV3DBase< T, PVType, FrameType >::perp(), mathSSE::sqrt(), runGCPTkAlMap::title, GlobalTrajectoryParameters::vector(), and LocalTrajectoryParameters::vector().

Referenced by analyzeTrackTrajectory(), muonFitter(), and trackFitter().

◆ debugTrajectorySOSv()

void GlobalTrackerMuonAlignment::debugTrajectorySOSv ( const std::string  title,
TrajectoryStateOnSurface  trajSOS 
)

Definition at line 3276 of file GlobalTrackerMuonAlignment.cc.

3276  {
3277  std::cout << " --- " << title << " --- " << std::endl;
3278  if (!trajSOS.isValid()) {
3279  std::cout << " Not valid !!!! " << std::endl;
3280  return;
3281  }
3282  std::cout << " R |p| " << trajSOS.globalPosition().perp() << " " << trajSOS.globalMomentum().mag() << " charge "
3283  << trajSOS.charge() << std::endl;
3284  std::cout << " x p " << trajSOS.globalParameters().vector()(0) << " " << trajSOS.globalParameters().vector()(1)
3285  << " " << trajSOS.globalParameters().vector()(2) << " " << trajSOS.globalParameters().vector()(3) << " "
3286  << trajSOS.globalParameters().vector()(4) << " " << trajSOS.globalParameters().vector()(5) << std::endl;
3287  std::cout << " +/- " << sqrt(trajSOS.cartesianError().matrix()(0, 0)) << " "
3288  << sqrt(trajSOS.cartesianError().matrix()(1, 1)) << " " << sqrt(trajSOS.cartesianError().matrix()(2, 2))
3289  << " " << sqrt(trajSOS.cartesianError().matrix()(3, 3)) << " "
3290  << sqrt(trajSOS.cartesianError().matrix()(4, 4)) << " " << sqrt(trajSOS.cartesianError().matrix()(5, 5))
3291  << std::endl;
3292  std::cout << " q/p dxy/dz xy " << trajSOS.localParameters().vector()(0) << " "
3293  << trajSOS.localParameters().vector()(1) << " " << trajSOS.localParameters().vector()(2) << " "
3294  << trajSOS.localParameters().vector()(3) << " " << trajSOS.localParameters().vector()(4) << std::endl;
3295  std::cout << " +/- error " << sqrt(trajSOS.localError().matrix()(0, 0)) << " "
3296  << sqrt(trajSOS.localError().matrix()(1, 1)) << " " << sqrt(trajSOS.localError().matrix()(2, 2)) << " "
3297  << sqrt(trajSOS.localError().matrix()(3, 3)) << " " << sqrt(trajSOS.localError().matrix()(4, 4)) << " "
3298  << std::endl;
3299 }

References TrajectoryStateOnSurface::cartesianError(), TrajectoryStateOnSurface::charge(), gather_cfg::cout, TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalParameters(), TrajectoryStateOnSurface::globalPosition(), TrajectoryStateOnSurface::isValid(), TrajectoryStateOnSurface::localError(), TrajectoryStateOnSurface::localParameters(), PV3DBase< T, PVType, FrameType >::mag(), CartesianTrajectoryError::matrix(), LocalTrajectoryError::matrix(), PV3DBase< T, PVType, FrameType >::perp(), mathSSE::sqrt(), runGCPTkAlMap::title, GlobalTrajectoryParameters::vector(), and LocalTrajectoryParameters::vector().

Referenced by debugTrajectory(), muonFitter(), and trackFitter().

◆ endJob()

void GlobalTrackerMuonAlignment::endJob ( void  )
overrideprivatevirtual

Reimplemented from edm::EDAnalyzer.

Definition at line 361 of file GlobalTrackerMuonAlignment.cc.

361  {
362  bool alarm = false;
363 
364  this->fitHist();
365 
366  AlgebraicVector3 d(0., 0., 0.); // ------------ alignmnet Global Algebraic
367 
368  AlgebraicSymMatrix33 InfI; // inverse it
369  for (int i = 0; i <= 2; i++)
370  for (int j = 0; j <= 2; j++) {
371  if (j < i)
372  continue;
373  InfI(i, j) += Inf(i, j);
374  }
375  bool ierr = !InfI.Invert();
376  if (ierr) {
377  if (alarm)
378  std::cout << " Error inverse Inf matrix !!!!!!!!!!!" << std::endl;
379  }
380 
381  for (int i = 0; i <= 2; i++)
382  for (int k = 0; k <= 2; k++)
383  d(i) -= InfI(i, k) * Gfr(k);
384  // end of Global Algebraic
385 
386  // --------------- alignment Global CLHEP
387  CLHEP::HepVector d3 = CLHEP::solve(Hess, -Grad);
388  int iEr3;
389  CLHEP::HepMatrix Errd3 = Hess.inverse(iEr3);
390  if (iEr3 != 0) {
391  if (alarm)
392  std::cout << " endJob Error inverse Hess matrix !!!!!!!!!!!" << std::endl;
393  }
394  // end of Global CLHEP
395 
396  // ----------------- alignment Local CLHEP
397  CLHEP::HepVector dLI = CLHEP::solve(HessL, -GradL);
398  int iErI;
399  CLHEP::HepMatrix ErrdLI = HessL.inverse(iErI);
400  if (iErI != 0) {
401  if (alarm)
402  std::cout << " endJob Error inverse HessL matrix !!!!!!!!!!!" << std::endl;
403  }
404  // end of Local CLHEP
405 
406  // printout of final parameters
407  std::cout << " ---- " << N_event << " event " << N_track << " tracks " << MuSelect << " ---- " << std::endl;
408  if (collectionIsolated)
409  std::cout << " ALCARECOMuAlCalIsolatedMu " << std::endl;
410  else if (collectionCosmic)
411  std::cout << " ALCARECOMuAlGlobalCosmics " << std::endl;
412  else
413  std::cout << smuonTags_ << std::endl;
414  std::cout << " Similated shifts[cm] " << MuGlShift(1) << " " << MuGlShift(2) << " " << MuGlShift(3) << " "
415  << " angles[rad] " << MuGlAngle(1) << " " << MuGlAngle(2) << " " << MuGlAngle(3) << " " << std::endl;
416  std::cout << " d " << -d << std::endl;
417  ;
418  std::cout << " +- " << sqrt(InfI(0, 0)) << " " << sqrt(InfI(1, 1)) << " " << sqrt(InfI(2, 2)) << std::endl;
419  std::cout << " dG " << d3(1) << " " << d3(2) << " " << d3(3) << " " << d3(4) << " " << d3(5) << " " << d3(6)
420  << std::endl;
421  ;
422  std::cout << " +- " << sqrt(Errd3(1, 1)) << " " << sqrt(Errd3(2, 2)) << " " << sqrt(Errd3(3, 3)) << " "
423  << sqrt(Errd3(4, 4)) << " " << sqrt(Errd3(5, 5)) << " " << sqrt(Errd3(6, 6)) << std::endl;
424  std::cout << " dL " << dLI(1) << " " << dLI(2) << " " << dLI(3) << " " << dLI(4) << " " << dLI(5) << " " << dLI(6)
425  << std::endl;
426  ;
427  std::cout << " +- " << sqrt(ErrdLI(1, 1)) << " " << sqrt(ErrdLI(2, 2)) << " " << sqrt(ErrdLI(3, 3)) << " "
428  << sqrt(ErrdLI(4, 4)) << " " << sqrt(ErrdLI(5, 5)) << " " << sqrt(ErrdLI(6, 6)) << std::endl;
429 
430  // what do we write to DB
431  CLHEP::HepVector vectorToDb(6, 0), vectorErrToDb(6, 0);
432  //vectorToDb = d3;
433  //for(unsigned int i=1; i<=6; i++) vectorErrToDb(i) = sqrt(Errd3(i,i));
434  vectorToDb = -dLI;
435  for (unsigned int i = 1; i <= 6; i++)
436  vectorErrToDb(i) = sqrt(ErrdLI(i, i));
437 
438  // write histograms to root file
439  file->Write();
440  file->Close();
441 
442  // write global parameters to text file
443  OutGlobalTxt.open(txtOutFile_.c_str(), ios::out);
444  if (!OutGlobalTxt.is_open())
445  std::cout << " outglobal.txt is not open !!!!!" << std::endl;
446  else {
447  OutGlobalTxt << vectorToDb(1) << " " << vectorToDb(2) << " " << vectorToDb(3) << " " << vectorToDb(4) << " "
448  << vectorToDb(5) << " " << vectorToDb(6) << " muon Global.\n";
449  OutGlobalTxt << vectorErrToDb(1) << " " << vectorErrToDb(1) << " " << vectorErrToDb(1) << " " << vectorErrToDb(1)
450  << " " << vectorErrToDb(1) << " " << vectorErrToDb(1) << " errors.\n";
451  OutGlobalTxt << N_event << " events are processed.\n";
452 
453  if (collectionIsolated)
454  OutGlobalTxt << "ALCARECOMuAlCalIsolatedMu.\n";
455  else if (collectionCosmic)
456  OutGlobalTxt << " ALCARECOMuAlGlobalCosmics.\n";
457  else
458  OutGlobalTxt << smuonTags_ << ".\n";
459  OutGlobalTxt.close();
460  std::cout << " Write to the file outglobal.txt done " << std::endl;
461  }
462 
463  // write new GlobalPositionRcd to DB
464  if (debug_)
465  std::cout << " writeBD_ " << writeDB_ << std::endl;
466  if (writeDB_)
467  GlobalTrackerMuonAlignment::writeGlPosRcd(vectorToDb);
468 }

References collectionCosmic, collectionIsolated, gather_cfg::cout, ztail::d, debug_, file, fitHist(), Gfr, Grad, GradL, Hess, HessL, mps_fire::i, Inf, dqmiolumiharvest::j, dqmdumpme::k, MuGlAngle, MuGlShift, MuSelect, N_event, N_track, MillePedeFileConverter_cfg::out, OutGlobalTxt, smuonTags_, mathSSE::sqrt(), txtOutFile_, writeDB_, and writeGlPosRcd().

Referenced by o2olib.O2ORunMgr::executeJob().

◆ fitHist()

void GlobalTrackerMuonAlignment::fitHist ( )

Definition at line 549 of file GlobalTrackerMuonAlignment.cc.

549  {
550  histo7->Fit("gaus", "Q");
551 
552  histo12->Fit("gaus", "Q");
553  histo13->Fit("gaus", "Q");
554  histo14->Fit("gaus", "Q");
555  histo15->Fit("gaus", "Q");
556  histo16->Fit("gaus", "Q");
557  histo17->Fit("gaus", "Q");
558 
559  histo21->Fit("gaus", "Q");
560  histo22->Fit("gaus", "Q");
561  histo23->Fit("gaus", "Q");
562  histo24->Fit("gaus", "Q");
563  histo25->Fit("gaus", "Q");
564  histo26->Fit("gaus", "Q");
565 
566  histo32->Fit("gaus", "Q");
567  histo33->Fit("gaus", "Q");
568  histo34->Fit("gaus", "Q");
569  histo35->Fit("gaus", "Q");
570 }

References histo12, histo13, histo14, histo15, histo16, histo17, histo21, histo22, histo23, histo24, histo25, histo26, histo32, histo33, histo34, histo35, and histo7.

Referenced by endJob().

◆ gradientGlobal()

void GlobalTrackerMuonAlignment::gradientGlobal ( GlobalVector &  GRt,
GlobalVector &  GPt,
GlobalVector &  GRm,
GlobalVector &  GPm,
GlobalVector &  GNorm,
AlgebraicSymMatrix66 &  GCov 
)

Definition at line 2123 of file GlobalTrackerMuonAlignment.cc.

2128  {
2129  // we search for 6D global correction vector (d, a), where
2130  // 3D vector of shihts d
2131  // 3D vector of rotation angles a
2132 
2133  //bool alarm = true;
2134  bool alarm = false;
2135  bool info = false;
2136 
2137  int Nd = 6; // dimension of vector of alignment pararmeters, d
2138 
2139  //double PtMom = GPt.mag();
2140  CLHEP::HepSymMatrix w(Nd, 0);
2141  for (int i = 1; i <= Nd; i++)
2142  for (int j = 1; j <= Nd; j++) {
2143  if (j <= i)
2144  w(i, j) = GCov(i - 1, j - 1);
2145  //if(i >= 3) w(i,j) /= PtMom;
2146  //if(j >= 3) w(i,j) /= PtMom;
2147  if ((i == j) && (i <= 3) && (GCov(i - 1, j - 1) < 1.e-20))
2148  w(i, j) = 1.e20; // w=0
2149  if (i != j)
2150  w(i, j) = 0.; // use diaginal elements
2151  }
2152 
2153  //GPt /= GPt.mag();
2154  //GPm /= GPm.mag(); // end of transform
2155 
2156  CLHEP::HepVector V(Nd), Rt(3), Pt(3), Rm(3), Pm(3), Norm(3);
2157  Rt(1) = GRt.x();
2158  Rt(2) = GRt.y();
2159  Rt(3) = GRt.z();
2160  Pt(1) = GPt.x();
2161  Pt(2) = GPt.y();
2162  Pt(3) = GPt.z();
2163  Rm(1) = GRm.x();
2164  Rm(2) = GRm.y();
2165  Rm(3) = GRm.z();
2166  Pm(1) = GPm.x();
2167  Pm(2) = GPm.y();
2168  Pm(3) = GPm.z();
2169  Norm(1) = GNorm.x();
2170  Norm(2) = GNorm.y();
2171  Norm(3) = GNorm.z();
2172 
2173  V = dsum(Rm - Rt, Pm - Pt); //std::cout<<" V "<<V.T()<<std::endl;
2174 
2175  //double PmN = CLHEP::dot(Pm, Norm);
2176  double PmN = CLHEP_dot(Pm, Norm);
2177 
2178  CLHEP::HepMatrix Jac(Nd, Nd, 0);
2179  for (int i = 1; i <= 3; i++)
2180  for (int j = 1; j <= 3; j++) {
2181  Jac(i, j) = Pm(i) * Norm(j) / PmN;
2182  if (i == j)
2183  Jac(i, j) -= 1.;
2184  }
2185 
2186  // dp/da
2187  Jac(4, 4) = 0.; // dpx/dax
2188  Jac(5, 4) = -Pm(3); // dpy/dax
2189  Jac(6, 4) = Pm(2); // dpz/dax
2190  Jac(4, 5) = Pm(3); // dpx/day
2191  Jac(5, 5) = 0.; // dpy/day
2192  Jac(6, 5) = -Pm(1); // dpz/day
2193  Jac(4, 6) = -Pm(2); // dpx/daz
2194  Jac(5, 6) = Pm(1); // dpy/daz
2195  Jac(6, 6) = 0.; // dpz/daz
2196 
2197  CLHEP::HepVector dsda(3);
2198  dsda(1) = (Norm(2) * Rm(3) - Norm(3) * Rm(2)) / PmN;
2199  dsda(2) = (Norm(3) * Rm(1) - Norm(1) * Rm(3)) / PmN;
2200  dsda(3) = (Norm(1) * Rm(2) - Norm(2) * Rm(1)) / PmN;
2201 
2202  // dr/da
2203  Jac(1, 4) = Pm(1) * dsda(1); // drx/dax
2204  Jac(2, 4) = -Rm(3) + Pm(2) * dsda(1); // dry/dax
2205  Jac(3, 4) = Rm(2) + Pm(3) * dsda(1); // drz/dax
2206 
2207  Jac(1, 5) = Rm(3) + Pm(1) * dsda(2); // drx/day
2208  Jac(2, 5) = Pm(2) * dsda(2); // dry/day
2209  Jac(3, 5) = -Rm(1) + Pm(3) * dsda(2); // drz/day
2210 
2211  Jac(1, 6) = -Rm(2) + Pm(1) * dsda(3); // drx/daz
2212  Jac(2, 6) = Rm(1) + Pm(2) * dsda(3); // dry/daz
2213  Jac(3, 6) = Pm(3) * dsda(3); // drz/daz
2214 
2215  CLHEP::HepSymMatrix W(Nd, 0);
2216  int ierr;
2217  W = w.inverse(ierr);
2218  if (ierr != 0) {
2219  if (alarm)
2220  std::cout << " gradientGlobal: inversion of matrix w fail " << std::endl;
2221  return;
2222  }
2223 
2224  CLHEP::HepMatrix W_Jac(Nd, Nd, 0);
2225  W_Jac = Jac.T() * W;
2226 
2227  CLHEP::HepVector grad3(Nd);
2228  grad3 = W_Jac * V;
2229 
2230  CLHEP::HepMatrix hess3(Nd, Nd);
2231  hess3 = Jac.T() * W * Jac;
2232  //hess3(4,4) = 1.e-10; hess3(5,5) = 1.e-10; hess3(6,6) = 1.e-10; //????????????????
2233 
2234  Grad += grad3;
2235  Hess += hess3;
2236 
2237  CLHEP::HepVector d3I = CLHEP::solve(Hess, -Grad);
2238  int iEr3I;
2239  CLHEP::HepMatrix Errd3I = Hess.inverse(iEr3I);
2240  if (iEr3I != 0) {
2241  if (alarm)
2242  std::cout << " gradientGlobal error inverse Hess matrix !!!!!!!!!!!" << std::endl;
2243  }
2244 
2245  if (info || debug_) {
2246  std::cout << " dG " << d3I(1) << " " << d3I(2) << " " << d3I(3) << " " << d3I(4) << " " << d3I(5) << " " << d3I(6)
2247  << std::endl;
2248  ;
2249  std::cout << " +- " << sqrt(Errd3I(1, 1)) << " " << sqrt(Errd3I(2, 2)) << " " << sqrt(Errd3I(3, 3)) << " "
2250  << sqrt(Errd3I(4, 4)) << " " << sqrt(Errd3I(5, 5)) << " " << sqrt(Errd3I(6, 6)) << std::endl;
2251  }
2252 
2253 #ifdef CHECK_OF_DERIVATIVES
2254  // -------------------- check of derivatives
2255 
2256  CLHEP::HepVector d(3, 0), a(3, 0);
2257  CLHEP::HepMatrix T(3, 3, 1);
2258 
2259  CLHEP::HepMatrix Ti = T.T();
2260  //double A = CLHEP::dot(Ti*Pm, Norm);
2261  //double B = CLHEP::dot((Rt -Ti*Rm + Ti*d), Norm);
2262  double A = CLHEP_dot(Ti * Pm, Norm);
2263  double B = CLHEP_dot((Rt - Ti * Rm + Ti * d), Norm);
2264  double s0 = B / A;
2265 
2266  CLHEP::HepVector r0(3, 0), p0(3, 0);
2267  r0 = Ti * Rm - Ti * d + s0 * (Ti * Pm) - Rt;
2268  p0 = Ti * Pm - Pt;
2269 
2270  double delta = 0.0001;
2271 
2272  int ii = 3;
2273  d(ii) += delta; // d
2274  //T(2,3) += delta; T(3,2) -= delta; int ii = 1; // a1
2275  //T(3,1) += delta; T(1,3) -= delta; int ii = 2; // a2
2276  //T(1,2) += delta; T(2,1) -= delta; int ii = 3; // a2
2277  Ti = T.T();
2278  //A = CLHEP::dot(Ti*Pm, Norm);
2279  //B = CLHEP::dot((Rt -Ti*Rm + Ti*d), Norm);
2280  A = CLHEP_dot(Ti * Pm, Norm);
2281  B = CLHEP_dot((Rt - Ti * Rm + Ti * d), Norm);
2282  double s = B / A;
2283 
2284  CLHEP::HepVector r(3, 0), p(3, 0);
2285  r = Ti * Rm - Ti * d + s * (Ti * Pm) - Rt;
2286  p = Ti * Pm - Pt;
2287 
2288  std::cout << " s0 s num dsda(" << ii << ") " << s0 << " " << s << " " << (s - s0) / delta << " " << dsda(ii) << endl;
2289  // d(r,p) / d shift
2290  std::cout << " -- An d(r,p)/d(" << ii << ") " << Jac(1, ii) << " " << Jac(2, ii) << " " << Jac(3, ii) << " "
2291  << Jac(4, ii) << " " << Jac(5, ii) << " " << Jac(6, ii) << std::endl;
2292  std::cout << " Nu d(r,p)/d(" << ii << ") " << (r(1) - r0(1)) / delta << " " << (r(2) - r0(2)) / delta << " "
2293  << (r(3) - r0(3)) / delta << " " << (p(1) - p0(1)) / delta << " " << (p(2) - p0(2)) / delta << " "
2294  << (p(3) - p0(3)) / delta << std::endl;
2295  // d(r,p) / d angle
2296  std::cout << " -- An d(r,p)/a(" << ii << ") " << Jac(1, ii + 3) << " " << Jac(2, ii + 3) << " " << Jac(3, ii + 3)
2297  << " " << Jac(4, ii + 3) << " " << Jac(5, ii + 3) << " " << Jac(6, ii + 3) << std::endl;
2298  std::cout << " Nu d(r,p)/a(" << ii << ") " << (r(1) - r0(1)) / delta << " " << (r(2) - r0(2)) / delta << " "
2299  << (r(3) - r0(3)) / delta << " " << (p(1) - p0(1)) / delta << " " << (p(2) - p0(2)) / delta << " "
2300  << (p(3) - p0(3)) / delta << std::endl;
2301  // ----------------------------- end of check
2302 #endif
2303 
2304  return;
2305 } // end gradientGlobal

References MaterialEffects_cfi::A, a, TtFullHadDaughter::B, CLHEP_dot(), gather_cfg::cout, ztail::d, debug_, dumpMFGeometry_cfg::delta, MillePedeFileConverter_cfg::e, Grad, Hess, mps_fire::i, cuy::ii, info(), dqmiolumiharvest::j, AlCaHLTBitMon_ParallelJobs::p, alignCSCRings::r, alignCSCRings::s, mathSSE::sqrt(), cms::cuda::V, w, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by analyzeTrackTrack(), and analyzeTrackTrajectory().

◆ gradientGlobalAlg()

void GlobalTrackerMuonAlignment::gradientGlobalAlg ( GlobalVector &  Rt,
GlobalVector &  Pt,
GlobalVector &  Rm,
GlobalVector &  Nl,
AlgebraicSymMatrix66 &  Cm 
)

Definition at line 2032 of file GlobalTrackerMuonAlignment.cc.

2033  {
2034  // ---------------------------- Calculate Information matrix and Gfree vector
2035  // Information == Hessian , Gfree == gradient of the objective function
2036 
2037  AlgebraicMatrix33 Jac;
2038  AlgebraicVector3 Wi, R_m, R_t, P_t, Norm, dR;
2039 
2040  R_m(0) = Rm.x();
2041  R_m(1) = Rm.y();
2042  R_m(2) = Rm.z();
2043  R_t(0) = Rt.x();
2044  R_t(1) = Rt.y();
2045  R_t(2) = Rt.z();
2046  P_t(0) = Pt.x();
2047  P_t(1) = Pt.y();
2048  P_t(2) = Pt.z();
2049  Norm(0) = Nl.x();
2050  Norm(1) = Nl.y();
2051  Norm(2) = Nl.z();
2052 
2053  for (int i = 0; i <= 2; i++) {
2054  if (Cm(i, i) > 1.e-20)
2055  Wi(i) = 1. / Cm(i, i);
2056  else
2057  Wi(i) = 1.e-10;
2058  dR(i) = R_m(i) - R_t(i);
2059  }
2060 
2061  float PtN = P_t(0) * Norm(0) + P_t(1) * Norm(1) + P_t(2) * Norm(2);
2062 
2063  Jac(0, 0) = 1. - P_t(0) * Norm(0) / PtN;
2064  Jac(0, 1) = -P_t(0) * Norm(1) / PtN;
2065  Jac(0, 2) = -P_t(0) * Norm(2) / PtN;
2066 
2067  Jac(1, 0) = -P_t(1) * Norm(0) / PtN;
2068  Jac(1, 1) = 1. - P_t(1) * Norm(1) / PtN;
2069  Jac(1, 2) = -P_t(1) * Norm(2) / PtN;
2070 
2071  Jac(2, 0) = -P_t(2) * Norm(0) / PtN;
2072  Jac(2, 1) = -P_t(2) * Norm(1) / PtN;
2073  Jac(2, 2) = 1. - P_t(2) * Norm(2) / PtN;
2074 
2075  AlgebraicSymMatrix33 Itr;
2076 
2077  for (int i = 0; i <= 2; i++)
2078  for (int j = 0; j <= 2; j++) {
2079  if (j < i)
2080  continue;
2081  Itr(i, j) = 0.;
2082  //std::cout<<" ij "<<i<<" "<<j<<std::endl;
2083  for (int k = 0; k <= 2; k++) {
2084  Itr(i, j) += Jac(k, i) * Wi(k) * Jac(k, j);
2085  }
2086  }
2087 
2088  for (int i = 0; i <= 2; i++)
2089  for (int j = 0; j <= 2; j++) {
2090  if (j < i)
2091  continue;
2092  Inf(i, j) += Itr(i, j);
2093  }
2094 
2095  AlgebraicVector3 Gtr(0., 0., 0.);
2096  for (int i = 0; i <= 2; i++)
2097  for (int k = 0; k <= 2; k++)
2098  Gtr(i) += dR(k) * Wi(k) * Jac(k, i);
2099  for (int i = 0; i <= 2; i++)
2100  Gfr(i) += Gtr(i);
2101 
2102  if (debug_) {
2103  std::cout << " Wi " << Wi << std::endl;
2104  std::cout << " N " << Norm << std::endl;
2105  std::cout << " P_t " << P_t << std::endl;
2106  std::cout << " (Pt*N) " << PtN << std::endl;
2107  std::cout << " dR " << dR << std::endl;
2108  std::cout << " +/- " << 1. / sqrt(Wi(0)) << " " << 1. / sqrt(Wi(1)) << " " << 1. / sqrt(Wi(2)) << " " << std::endl;
2109  std::cout << " Jacobian dr/ddx " << std::endl;
2110  std::cout << Jac << std::endl;
2111  std::cout << " G-- " << Gtr << std::endl;
2112  std::cout << " Itrack " << std::endl;
2113  std::cout << Itr << std::endl;
2114  std::cout << " Gfr " << Gfr << std::endl;
2115  std::cout << " -- Inf --" << std::endl;
2116  std::cout << Inf << std::endl;
2117  }
2118 
2119  return;
2120 }

References gather_cfg::cout, debug_, HGC3DClusterGenMatchSelector_cfi::dR, MillePedeFileConverter_cfg::e, Gfr, mps_fire::i, Inf, dqmiolumiharvest::j, dqmdumpme::k, mathSSE::sqrt(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by analyzeTrackTrack(), and analyzeTrackTrajectory().

◆ gradientLocal()

void GlobalTrackerMuonAlignment::gradientLocal ( GlobalVector &  GRt,
GlobalVector &  GPt,
GlobalVector &  GRm,
GlobalVector &  GPm,
GlobalVector &  GNorm,
CLHEP::HepSymMatrix &  covLoc,
CLHEP::HepMatrix &  rotLoc,
CLHEP::HepVector &  R0,
AlgebraicVector4 &  LPRm 
)

Definition at line 2308 of file GlobalTrackerMuonAlignment.cc.

2316  {
2317  // we search for 6D global correction vector (d, a), where
2318  // 3D vector of shihts d
2319  // 3D vector of rotation angles a
2320 
2321  bool alarm = true;
2322  //bool alarm = false;
2323  bool info = false;
2324 
2325  if (debug_)
2326  std::cout << " gradientLocal " << std::endl;
2327 
2328  /*
2329  const Surface& refSurface = tsosMuon.surface();
2330 
2331  CLHEP::HepMatrix rotLoc (3,3,0);
2332  rotLoc(1,1) = refSurface.rotation().xx();
2333  rotLoc(1,2) = refSurface.rotation().xy();
2334  rotLoc(1,3) = refSurface.rotation().xz();
2335 
2336  rotLoc(2,1) = refSurface.rotation().yx();
2337  rotLoc(2,2) = refSurface.rotation().yy();
2338  rotLoc(2,3) = refSurface.rotation().yz();
2339 
2340  rotLoc(3,1) = refSurface.rotation().zx();
2341  rotLoc(3,2) = refSurface.rotation().zy();
2342  rotLoc(3,3) = refSurface.rotation().zz();
2343  */
2344 
2345  CLHEP::HepVector Rt(3), Pt(3), Rm(3), Pm(3), Norm(3);
2346  Rt(1) = GRt.x();
2347  Rt(2) = GRt.y();
2348  Rt(3) = GRt.z();
2349  Pt(1) = GPt.x();
2350  Pt(2) = GPt.y();
2351  Pt(3) = GPt.z();
2352  Rm(1) = GRm.x();
2353  Rm(2) = GRm.y();
2354  Rm(3) = GRm.z();
2355  Pm(1) = GPm.x();
2356  Pm(2) = GPm.y();
2357  Pm(3) = GPm.z();
2358  Norm(1) = GNorm.x();
2359  Norm(2) = GNorm.y();
2360  Norm(3) = GNorm.z();
2361 
2362  CLHEP::HepVector V(4), Rml(3), Pml(3), Rtl(3), Ptl(3);
2363 
2364  /*
2365  R0(1) = refSurface.position().x();
2366  R0(2) = refSurface.position().y();
2367  R0(3) = refSurface.position().z();
2368  */
2369 
2370  Rml = rotLoc * (Rm - R0);
2371  Rtl = rotLoc * (Rt - R0);
2372  Pml = rotLoc * Pm;
2373  Ptl = rotLoc * Pt;
2374 
2375  V(1) = LPRm(0) - Ptl(1) / Ptl(3);
2376  V(2) = LPRm(1) - Ptl(2) / Ptl(3);
2377  V(3) = LPRm(2) - Rtl(1);
2378  V(4) = LPRm(3) - Rtl(2);
2379 
2380  /*
2381  CLHEP::HepSymMatrix Cov(4,0), W(4,0);
2382  for(int i=1; i<=4; i++)
2383  for(int j=1; j<=i; j++){
2384  Cov(i,j) = (tsosMuon.localError().matrix()
2385  + tsosTrack.localError().matrix())(i,j);
2386  //if(i != j) Cov(i,j) = 0.;
2387  //if((i == j) && ((i==1) || (i==2))) Cov(i,j) = 100.;
2388  //if((i == j) && ((i==3) || (i==4))) Cov(i,j) = 10000.;
2389  }
2390  W = Cov;
2391  */
2392 
2393  CLHEP::HepSymMatrix W = covLoc;
2394 
2395  int ierr;
2396  W.invert(ierr);
2397  if (ierr != 0) {
2398  if (alarm)
2399  std::cout << " gradientLocal: inversion of matrix W fail " << std::endl;
2400  return;
2401  }
2402 
2403  // JacobianCartesianToLocal
2404 
2405  //AlgebraicMatrix56 jacobian // differ from calculation above
2406  //= JacobianCartesianToLocal::JacobianCartesianToLocal
2407  //(refSurface, tsosTrack.localParameters()).jacobian();
2408  //for(int i=1; i<=4; i++) for(int j=1; j<=6; j++){
2409  //int j1 = j - 1; JacToLoc(i,j) = jacobian(i, j1);}
2410 
2411  CLHEP::HepMatrix JacToLoc(4, 6, 0);
2412  for (int i = 1; i <= 2; i++)
2413  for (int j = 1; j <= 3; j++) {
2414  JacToLoc(i, j + 3) = (rotLoc(i, j) - rotLoc(3, j) * Pml(i) / Pml(3)) / Pml(3);
2415  JacToLoc(i + 2, j) = rotLoc(i, j);
2416  }
2417 
2418  // JacobianCorrectionsToCartesian
2419  //double PmN = CLHEP::dot(Pm, Norm);
2420  double PmN = CLHEP_dot(Pm, Norm);
2421 
2422  CLHEP::HepMatrix Jac(6, 6, 0);
2423  for (int i = 1; i <= 3; i++)
2424  for (int j = 1; j <= 3; j++) {
2425  Jac(i, j) = Pm(i) * Norm(j) / PmN;
2426  if (i == j)
2427  Jac(i, j) -= 1.;
2428  }
2429 
2430  // dp/da
2431  Jac(4, 4) = 0.; // dpx/dax
2432  Jac(5, 4) = -Pm(3); // dpy/dax
2433  Jac(6, 4) = Pm(2); // dpz/dax
2434  Jac(4, 5) = Pm(3); // dpx/day
2435  Jac(5, 5) = 0.; // dpy/day
2436  Jac(6, 5) = -Pm(1); // dpz/day
2437  Jac(4, 6) = -Pm(2); // dpx/daz
2438  Jac(5, 6) = Pm(1); // dpy/daz
2439  Jac(6, 6) = 0.; // dpz/daz
2440 
2441  CLHEP::HepVector dsda(3);
2442  dsda(1) = (Norm(2) * Rm(3) - Norm(3) * Rm(2)) / PmN;
2443  dsda(2) = (Norm(3) * Rm(1) - Norm(1) * Rm(3)) / PmN;
2444  dsda(3) = (Norm(1) * Rm(2) - Norm(2) * Rm(1)) / PmN;
2445 
2446  // dr/da
2447  Jac(1, 4) = Pm(1) * dsda(1); // drx/dax
2448  Jac(2, 4) = -Rm(3) + Pm(2) * dsda(1); // dry/dax
2449  Jac(3, 4) = Rm(2) + Pm(3) * dsda(1); // drz/dax
2450 
2451  Jac(1, 5) = Rm(3) + Pm(1) * dsda(2); // drx/day
2452  Jac(2, 5) = Pm(2) * dsda(2); // dry/day
2453  Jac(3, 5) = -Rm(1) + Pm(3) * dsda(2); // drz/day
2454 
2455  Jac(1, 6) = -Rm(2) + Pm(1) * dsda(3); // drx/daz
2456  Jac(2, 6) = Rm(1) + Pm(2) * dsda(3); // dry/daz
2457  Jac(3, 6) = Pm(3) * dsda(3); // drz/daz
2458 
2459  // JacobianCorrectionToLocal
2460  CLHEP::HepMatrix JacCorLoc(4, 6, 0);
2461  JacCorLoc = JacToLoc * Jac;
2462 
2463  // gradient and Hessian
2464  CLHEP::HepMatrix W_Jac(6, 4, 0);
2465  W_Jac = JacCorLoc.T() * W;
2466 
2467  CLHEP::HepVector gradL(6);
2468  gradL = W_Jac * V;
2469 
2470  CLHEP::HepMatrix hessL(6, 6);
2471  hessL = JacCorLoc.T() * W * JacCorLoc;
2472 
2473  GradL += gradL;
2474  HessL += hessL;
2475 
2476  CLHEP::HepVector dLI = CLHEP::solve(HessL, -GradL);
2477  int iErI;
2478  CLHEP::HepMatrix ErrdLI = HessL.inverse(iErI);
2479  if (iErI != 0) {
2480  if (alarm)
2481  std::cout << " gradLocal Error inverse Hess matrix !!!!!!!!!!!" << std::endl;
2482  }
2483 
2484  if (info || debug_) {
2485  std::cout << " dL " << dLI(1) << " " << dLI(2) << " " << dLI(3) << " " << dLI(4) << " " << dLI(5) << " " << dLI(6)
2486  << std::endl;
2487  ;
2488  std::cout << " +- " << sqrt(ErrdLI(1, 1)) << " " << sqrt(ErrdLI(2, 2)) << " " << sqrt(ErrdLI(3, 3)) << " "
2489  << sqrt(ErrdLI(4, 4)) << " " << sqrt(ErrdLI(5, 5)) << " " << sqrt(ErrdLI(6, 6)) << std::endl;
2490  }
2491 
2492  if (debug_) {
2493  //std::cout<<" dV(da3) {"<<-JacCorLoc(1,6)*0.002<<" "<<-JacCorLoc(2,6)*0.002<<" "
2494  // <<-JacCorLoc(3,6)*0.002<<" "<<-JacCorLoc(4,6)*0.002<<"}"<<std::endl;
2495  //std::cout<<" JacCLo {"<<JacCorLoc(1,6)<<" "<<JacCorLoc(2,6)<<" "
2496  // <<JacCorLoc(3,6)<<" "<<JacCorLoc(4,6)<<"}"<<std::endl;
2497  //std::cout<<"Jpx/yx {"<<Jac(4,6)/Pm(3)<<" "<<Jac(5,6)/Pm(3)<<" "
2498  // <<Jac(1,6)<<" "<<Jac(2,6)<<"}"<<std::endl;
2499  //std::cout<<"Jac(,a3){"<<Jac(1,6)<<" "<<Jac(2,6)<<" "<<Jac(3,6)<<" "<<Jac(4,6)<<" "
2500  // <<Jac(5,6)<<" "<<Jac(6,6)<<std::endl;
2501  int i = 5;
2502  if (GNorm.z() > 0.95)
2503  std::cout << " Ecap1 N " << GNorm << std::endl;
2504  else if (GNorm.z() < -0.95)
2505  std::cout << " Ecap2 N " << GNorm << std::endl;
2506  else
2507  std::cout << " Barrel N " << GNorm << std::endl;
2508  std::cout << " JacCLo(i," << i << ") = {" << JacCorLoc(1, i) << " " << JacCorLoc(2, i) << " " << JacCorLoc(3, i)
2509  << " " << JacCorLoc(4, i) << "}" << std::endl;
2510  std::cout << " rotLoc " << rotLoc << std::endl;
2511  std::cout << " position " << R0 << std::endl;
2512  std::cout << " Pm,l " << Pm.T() << " " << Pml(1) / Pml(3) << " " << Pml(2) / Pml(3) << std::endl;
2513  std::cout << " Pt,l " << Pt.T() << " " << Ptl(1) / Ptl(3) << " " << Ptl(2) / Ptl(3) << std::endl;
2514  std::cout << " V " << V.T() << std::endl;
2515  std::cout << " Cov \n" << covLoc << std::endl;
2516  std::cout << " W*Cov " << W * covLoc << std::endl;
2517  //std::cout<<" JacCarToLoc = drldrc \n"<<
2518  //JacobianCartesianToLocal::JacobianCartesianToLocal
2519  //(refSurface, tsosTrack.localParameters()).jacobian()<<std::endl;
2520  std::cout << " JacToLoc " << JacToLoc << std::endl;
2521  }
2522 
2523 #ifdef CHECK_OF_JACOBIAN_CARTESIAN_TO_LOCAL
2524  //---------------------- check of derivatives
2525  CLHEP::HepVector V0(4, 0);
2526  V0(1) = Pml(1) / Pml(3) - Ptl(1) / Ptl(3);
2527  V0(2) = Pml(2) / Pml(3) - Ptl(2) / Ptl(3);
2528  V0(3) = Rml(1) - Rtl(1);
2529  V0(4) = Rml(2) - Rtl(2);
2530  int ii = 3;
2531  float delta = 0.01;
2532  CLHEP::HepVector V1(4, 0);
2533  if (ii <= 3) {
2534  Rm(ii) += delta;
2535  Rml = rotLoc * (Rm - R0);
2536  } else {
2537  Pm(ii - 3) += delta;
2538  Pml = rotLoc * Pm;
2539  }
2540  //if(ii <= 3) {Rt(ii) += delta; Rtl = rotLoc * (Rt - R0);}
2541  //else {Pt(ii-3) += delta; Ptl = rotLoc * Pt;}
2542  V1(1) = Pml(1) / Pml(3) - Ptl(1) / Ptl(3);
2543  V1(2) = Pml(2) / Pml(3) - Ptl(2) / Ptl(3);
2544  V1(3) = Rml(1) - Rtl(1);
2545  V1(4) = Rml(2) - Rtl(2);
2546 
2547  if (GNorm.z() > 0.95)
2548  std::cout << " Ecap1 N " << GNorm << std::endl;
2549  else if (GNorm.z() < -0.95)
2550  std::cout << " Ecap2 N " << GNorm << std::endl;
2551  else
2552  std::cout << " Barrel N " << GNorm << std::endl;
2553  std::cout << " dldc Num (i," << ii << ") " << (V1(1) - V0(1)) / delta << " " << (V1(2) - V0(2)) / delta << " "
2554  << (V1(3) - V0(3)) / delta << " " << (V1(4) - V0(4)) / delta << std::endl;
2555  std::cout << " dldc Ana (i," << ii << ") " << JacToLoc(1, ii) << " " << JacToLoc(2, ii) << " " << JacToLoc(3, ii)
2556  << " " << JacToLoc(4, ii) << std::endl;
2557  //float dtxdpx = (rotLoc(1,1) - rotLoc(3,1)*Pml(1)/Pml(3))/Pml(3);
2558  //float dtydpx = (rotLoc(2,1) - rotLoc(3,2)*Pml(2)/Pml(3))/Pml(3);
2559  //std::cout<<" dtx/dpx dty/ "<<dtxdpx<<" "<<dtydpx<<std::endl;
2560 #endif
2561 
2562  return;
2563 } // end gradientLocal

References CLHEP_dot(), gather_cfg::cout, debug_, dumpMFGeometry_cfg::delta, GradL, HessL, mps_fire::i, cuy::ii, info(), dqmiolumiharvest::j, HLT_FULL_cff::R0, mathSSE::sqrt(), cms::cuda::V, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by analyzeTrackTrack(), and analyzeTrackTrajectory().

◆ misalignMuon()

void GlobalTrackerMuonAlignment::misalignMuon ( GlobalVector &  GRm,
GlobalVector &  GPm,
GlobalVector &  Nl,
GlobalVector &  Rt,
GlobalVector &  Rm,
GlobalVector &  Pm 
)

Definition at line 2566 of file GlobalTrackerMuonAlignment.cc.

2567  {
2568  CLHEP::HepVector d(3, 0), a(3, 0), Norm(3), Rm0(3), Pm0(3), Rm1(3), Pm1(3), Rmi(3), Pmi(3);
2569 
2570  d(1) = 0.0;
2571  d(2) = 0.0;
2572  d(3) = 0.0;
2573  a(1) = 0.0000;
2574  a(2) = 0.0000;
2575  a(3) = 0.0000; // zero
2576  //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
2577  //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
2578  //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
2579  //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
2580  //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
2581  //d(1)=10.; d(2)=20.; d(3)=30.; a(1)=0.0100; a(2)=0.0200; a(3)=0.0300; // huge
2582  //d(1)=10.; d(2)=20.; d(3)=30.; a(1)=0.2000; a(2)=0.2500; a(3)=0.3000; // huge angles
2583  //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
2584  //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
2585  //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
2586 
2587  MuGlShift = d;
2588  MuGlAngle = a;
2589 
2590  if ((d(1) == 0.) & (d(2) == 0.) & (d(3) == 0.) & (a(1) == 0.) & (a(2) == 0.) & (a(3) == 0.)) {
2591  Rm = GRm;
2592  Pm = GPm;
2593  if (debug_)
2594  std::cout << " ...... without misalignment " << std::endl;
2595  return;
2596  }
2597 
2598  Rm0(1) = GRm.x();
2599  Rm0(2) = GRm.y();
2600  Rm0(3) = GRm.z();
2601  Pm0(1) = GPm.x();
2602  Pm0(2) = GPm.y();
2603  Pm0(3) = GPm.z();
2604  Norm(1) = Nl.x();
2605  Norm(2) = Nl.y();
2606  Norm(3) = Nl.z();
2607 
2608  CLHEP::HepMatrix T(3, 3, 1);
2609 
2610  //T(1,2) = a(3); T(1,3) = -a(2);
2611  //T(2,1) = -a(3); T(2,3) = a(1);
2612  //T(3,1) = a(2); T(3,2) = -a(1);
2613 
2614  double s1 = std::sin(a(1)), c1 = std::cos(a(1));
2615  double s2 = std::sin(a(2)), c2 = std::cos(a(2));
2616  double s3 = std::sin(a(3)), c3 = std::cos(a(3));
2617  T(1, 1) = c2 * c3;
2618  T(1, 2) = c1 * s3 + s1 * s2 * c3;
2619  T(1, 3) = s1 * s3 - c1 * s2 * c3;
2620  T(2, 1) = -c2 * s3;
2621  T(2, 2) = c1 * c3 - s1 * s2 * s3;
2622  T(2, 3) = s1 * c3 + c1 * s2 * s3;
2623  T(3, 1) = s2;
2624  T(3, 2) = -s1 * c2;
2625  T(3, 3) = c1 * c2;
2626 
2627  //int terr;
2628  //CLHEP::HepMatrix Ti = T.inverse(terr);
2629  CLHEP::HepMatrix Ti = T.T();
2630  CLHEP::HepVector di = -Ti * d;
2631 
2632  CLHEP::HepVector ex0(3, 0), ey0(3, 0), ez0(3, 0), ex(3, 0), ey(3, 0), ez(3, 0);
2633  ex0(1) = 1.;
2634  ey0(2) = 1.;
2635  ez0(3) = 1.;
2636  ex = Ti * ex0;
2637  ey = Ti * ey0;
2638  ez = Ti * ez0;
2639 
2640  CLHEP::HepVector TiN = Ti * Norm;
2641  //double si = CLHEP::dot((Ti*Rm0 - Rm0 + di), TiN) / CLHEP::dot(Pm0, TiN);
2642  //Rm1(1) = CLHEP::dot(ex, Rm0 + si*Pm0 - di);
2643  //Rm1(2) = CLHEP::dot(ey, Rm0 + si*Pm0 - di);
2644  //Rm1(3) = CLHEP::dot(ez, Rm0 + si*Pm0 - di);
2645  double si = CLHEP_dot((Ti * Rm0 - Rm0 + di), TiN) / CLHEP_dot(Pm0, TiN);
2646  Rm1(1) = CLHEP_dot(ex, Rm0 + si * Pm0 - di);
2647  Rm1(2) = CLHEP_dot(ey, Rm0 + si * Pm0 - di);
2648  Rm1(3) = CLHEP_dot(ez, Rm0 + si * Pm0 - di);
2649  Pm1 = T * Pm0;
2650 
2651  Rm = GlobalVector(Rm1(1), Rm1(2), Rm1(3));
2652  //Pm = GlobalVector(CLHEP::dot(Pm0,ex), CLHEP::dot(Pm0, ey), CLHEP::dot(Pm0,ez));// =T*Pm0
2653  Pm = GlobalVector(CLHEP_dot(Pm0, ex), CLHEP_dot(Pm0, ey), CLHEP_dot(Pm0, ez)); // =T*Pm0
2654 
2655  if (debug_) { // ------------- debug tranformation
2656 
2657  std::cout << " ----- Pm " << Pm << std::endl;
2658  std::cout << " T*Pm0 " << Pm1.T() << std::endl;
2659  //std::cout<<" Ti*T*Pm0 "<<(Ti*(T*Pm0)).T()<<std::endl;
2660 
2661  //CLHEP::HepVector Rml = (Rm0 + si*Pm0 - di) + di;
2662  CLHEP::HepVector Rml = Rm1(1) * ex + Rm1(2) * ey + Rm1(3) * ez + di;
2663 
2664  CLHEP::HepVector Rt0(3);
2665  Rt0(1) = Rt.x();
2666  Rt0(2) = Rt.y();
2667  Rt0(3) = Rt.z();
2668 
2669  //double sl = CLHEP::dot(T*(Rt0 - Rml), T*Norm) / CLHEP::dot(Ti * Pm1, Norm);
2670  double sl = CLHEP_dot(T * (Rt0 - Rml), T * Norm) / CLHEP_dot(Ti * Pm1, Norm);
2671  CLHEP::HepVector Rl = Rml + sl * (Ti * Pm1);
2672 
2673  //double A = CLHEP::dot(Ti*Pm1, Norm);
2674  //double B = CLHEP::dot(Rt0, Norm) - CLHEP::dot(Ti*Rm1, Norm)
2675  //+ CLHEP::dot(Ti*d, Norm);
2676  double A = CLHEP_dot(Ti * Pm1, Norm);
2677  double B = CLHEP_dot(Rt0, Norm) - CLHEP_dot(Ti * Rm1, Norm) + CLHEP_dot(Ti * d, Norm);
2678  double s = B / A;
2679  CLHEP::HepVector Dr = Ti * Rm1 - Ti * d + s * (Ti * Pm1) - Rm0;
2680  CLHEP::HepVector Dp = Ti * Pm1 - Pm0;
2681 
2682  CLHEP::HepVector TiR = Ti * Rm0 + di;
2683  CLHEP::HepVector Ri = T * TiR + d;
2684 
2685  std::cout << " --TTi-0 " << (Ri - Rm0).T() << std::endl;
2686  std::cout << " -- Dr " << Dr.T() << endl;
2687  std::cout << " -- Dp " << Dp.T() << endl;
2688  //std::cout<<" ex "<<ex.T()<<endl;
2689  //std::cout<<" ey "<<ey.T()<<endl;
2690  //std::cout<<" ez "<<ez.T()<<endl;
2691  //std::cout<<" ---- T ---- "<<T<<std::endl;
2692  //std::cout<<" ---- Ti ---- "<<Ti<<std::endl;
2693  //std::cout<<" ---- d ---- "<<d.T()<<std::endl;
2694  //std::cout<<" ---- di ---- "<<di.T()<<std::endl;
2695  std::cout << " -- si sl s " << si << " " << sl << " " << s << std::endl;
2696  //std::cout<<" -- si sl "<<si<<" "<<sl<<endl;
2697  //std::cout<<" -- si s "<<si<<" "<<s<<endl;
2698  std::cout << " -- Rml-(Rm0+si*Pm0) " << (Rml - Rm0 - si * Pm0).T() << std::endl;
2699  //std::cout<<" -- Rm0 "<<Rm0.T()<<std::endl;
2700  //std::cout<<" -- Rm1 "<<Rm1.T()<<std::endl;
2701  //std::cout<<" -- Rmi "<<Rmi.T()<<std::endl;
2702  //std::cout<<" --siPm "<<(si*Pm0).T()<<std::endl;
2703  //std::cout<<" --s2Pm "<<(s2*(T * Pm1)).T()<<std::endl;
2704  //std::cout<<" --R1-0 "<<(Rm1-Rm0).T()<<std::endl;
2705  //std::cout<<" --Ri-0 "<<(Rmi-Rm0).T()<<std::endl;
2706  std::cout << " --Rl-0 " << (Rl - Rm0).T() << std::endl;
2707  //std::cout<<" --Pi-0 "<<(Pmi-Pm0).T()<<std::endl;
2708  } // -------- end of debug
2709 
2710  return;
2711 
2712 } // end of misalignMuon

References MaterialEffects_cfi::A, a, TtFullHadDaughter::B, alignmentValidation::c1, CLHEP_dot(), funct::cos(), gather_cfg::cout, ztail::d, debug_, MuGlAngle, MuGlShift, alignCSCRings::s, indexGen::s2, funct::sin(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

◆ misalignMuonL()

void GlobalTrackerMuonAlignment::misalignMuonL ( GlobalVector &  GRm,
GlobalVector &  GPm,
GlobalVector &  Nl,
GlobalVector &  Rt,
GlobalVector &  Rm,
GlobalVector &  Pm,
AlgebraicVector4 &  Vm,
TrajectoryStateOnSurface &  tsosTrack,
TrajectoryStateOnSurface &  tsosMuon 
)

Definition at line 2715 of file GlobalTrackerMuonAlignment.cc.

2723  {
2724  CLHEP::HepVector d(3, 0), a(3, 0), Norm(3), Rm0(3), Pm0(3), Rm1(3), Pm1(3), Rmi(3), Pmi(3);
2725 
2726  d(1) = 0.0;
2727  d(2) = 0.0;
2728  d(3) = 0.0;
2729  a(1) = 0.0000;
2730  a(2) = 0.0000;
2731  a(3) = 0.0000; // zero
2732  //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
2733  //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
2734  //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
2735  //d(1)=0.0; d(2)=0.0; d(3)=0.0; a(1)=0.0000; a(2)=0.0000; a(3)=0.0020; // a3=.002
2736  //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
2737  //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
2738  //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
2739  //d(1)=10.; d(2)=20.; d(3)=30.; a(1)=0.0100; a(2)=0.0200; a(3)=0.0300; // huge
2740  //d(1)=10.; d(2)=20.; d(3)=30.; a(1)=0.2000; a(2)=0.2500; a(3)=0.3000; // huge angles
2741  //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
2742  //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
2743  //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
2744 
2745  MuGlShift = d;
2746  MuGlAngle = a;
2747 
2748  if ((d(1) == 0.) & (d(2) == 0.) & (d(3) == 0.) & (a(1) == 0.) & (a(2) == 0.) & (a(3) == 0.)) {
2749  Rm = GRm;
2750  Pm = GPm;
2751  AlgebraicVector4 Vm0;
2752  Vm = AlgebraicVector4(tsosMuon.localParameters().vector()(1),
2753  tsosMuon.localParameters().vector()(2),
2754  tsosMuon.localParameters().vector()(3),
2755  tsosMuon.localParameters().vector()(4));
2756  if (debug_)
2757  std::cout << " ...... without misalignment " << std::endl;
2758  return;
2759  }
2760 
2761  Rm0(1) = GRm.x();
2762  Rm0(2) = GRm.y();
2763  Rm0(3) = GRm.z();
2764  Pm0(1) = GPm.x();
2765  Pm0(2) = GPm.y();
2766  Pm0(3) = GPm.z();
2767  Norm(1) = Nl.x();
2768  Norm(2) = Nl.y();
2769  Norm(3) = Nl.z();
2770 
2771  CLHEP::HepMatrix T(3, 3, 1);
2772 
2773  //T(1,2) = a(3); T(1,3) = -a(2);
2774  //T(2,1) = -a(3); T(2,3) = a(1);
2775  //T(3,1) = a(2); T(3,2) = -a(1);
2776 
2777  double s1 = std::sin(a(1)), c1 = std::cos(a(1));
2778  double s2 = std::sin(a(2)), c2 = std::cos(a(2));
2779  double s3 = std::sin(a(3)), c3 = std::cos(a(3));
2780  T(1, 1) = c2 * c3;
2781  T(1, 2) = c1 * s3 + s1 * s2 * c3;
2782  T(1, 3) = s1 * s3 - c1 * s2 * c3;
2783  T(2, 1) = -c2 * s3;
2784  T(2, 2) = c1 * c3 - s1 * s2 * s3;
2785  T(2, 3) = s1 * c3 + c1 * s2 * s3;
2786  T(3, 1) = s2;
2787  T(3, 2) = -s1 * c2;
2788  T(3, 3) = c1 * c2;
2789 
2790  // Ti, di what we have to apply for misalignment
2791  //int terr;
2792  //CLHEP::HepMatrix Ti = T.inverse(terr);
2793  CLHEP::HepMatrix Ti = T.T();
2794  CLHEP::HepVector di = -Ti * d;
2795 
2796  CLHEP::HepVector ex0(3, 0), ey0(3, 0), ez0(3, 0), ex(3, 0), ey(3, 0), ez(3, 0);
2797  ex0(1) = 1.;
2798  ey0(2) = 1.;
2799  ez0(3) = 1.;
2800  ex = Ti * ex0;
2801  ey = Ti * ey0;
2802  ez = Ti * ez0;
2803 
2804  CLHEP::HepVector TiN = Ti * Norm;
2805  //double si = CLHEP::dot((Ti*Rm0 - Rm0 + di), TiN) / CLHEP::dot(Pm0, TiN);
2806  //Rm1(1) = CLHEP::dot(ex, Rm0 + si*Pm0 - di);
2807  //Rm1(2) = CLHEP::dot(ey, Rm0 + si*Pm0 - di);
2808  //Rm1(3) = CLHEP::dot(ez, Rm0 + si*Pm0 - di);
2809  double si = CLHEP_dot((Ti * Rm0 - Rm0 + di), TiN) / CLHEP_dot(Pm0, TiN);
2810  Rm1(1) = CLHEP_dot(ex, Rm0 + si * Pm0 - di);
2811  Rm1(2) = CLHEP_dot(ey, Rm0 + si * Pm0 - di);
2812  Rm1(3) = CLHEP_dot(ez, Rm0 + si * Pm0 - di);
2813  Pm1 = T * Pm0;
2814 
2815  // Global muon parameters after misalignment
2816  Rm = GlobalVector(Rm1(1), Rm1(2), Rm1(3));
2817  //Pm = GlobalVector(CLHEP::dot(Pm0,ex), CLHEP::dot(Pm0, ey), CLHEP::dot(Pm0,ez));// =T*Pm0
2818  Pm = GlobalVector(CLHEP_dot(Pm0, ex), CLHEP_dot(Pm0, ey), CLHEP_dot(Pm0, ez)); // =T*Pm0
2819 
2820  // Local muon parameters after misalignment
2821  const Surface& refSurface = tsosMuon.surface();
2822  CLHEP::HepMatrix Tl(3, 3, 0);
2823 
2824  Tl(1, 1) = refSurface.rotation().xx();
2825  Tl(1, 2) = refSurface.rotation().xy();
2826  Tl(1, 3) = refSurface.rotation().xz();
2827 
2828  Tl(2, 1) = refSurface.rotation().yx();
2829  Tl(2, 2) = refSurface.rotation().yy();
2830  Tl(2, 3) = refSurface.rotation().yz();
2831 
2832  Tl(3, 1) = refSurface.rotation().zx();
2833  Tl(3, 2) = refSurface.rotation().zy();
2834  Tl(3, 3) = refSurface.rotation().zz();
2835 
2836  CLHEP::HepVector R0(3, 0), newR0(3, 0), newRl(3, 0), newPl(3, 0);
2837  R0(1) = refSurface.position().x();
2838  R0(2) = refSurface.position().y();
2839  R0(3) = refSurface.position().z();
2840 
2841  newRl = Tl * (Rm1 - R0);
2842  newPl = Tl * Pm1;
2843 
2844  Vm(0) = newPl(1) / newPl(3);
2845  Vm(1) = newPl(2) / newPl(3);
2846  Vm(2) = newRl(1);
2847  Vm(3) = newRl(2);
2848 
2849 #ifdef CHECK_DERIVATIVES_LOCAL_VS_ANGLES
2850  float del = 0.0001;
2851  //int ii = 1; T(3,2) +=del; T(2,3) -=del;
2852  int ii = 2;
2853  T(3, 1) -= del;
2854  T(1, 3) += del;
2855  //int ii = 3; T(1,2) -=del; T(2,1) +=del;
2856  AlgebraicVector4 Vm0 = Vm;
2857  CLHEP::HepVector Rm10 = Rm1, Pm10 = Pm1;
2858  ;
2859  CLHEP::HepMatrix newTl = Tl * T;
2860  Ti = T.T();
2861  di = -Ti * d;
2862  ex = Ti * ex0;
2863  ey = Ti * ey0;
2864  ez = Ti * ez0;
2865  TiN = Ti * Norm;
2866  si = CLHEP_dot((Ti * Rm0 - Rm0 + di), TiN) / CLHEP_dot(Pm0, TiN);
2867  Rm1(1) = CLHEP_dot(ex, Rm0 + si * Pm0 - di);
2868  Rm1(2) = CLHEP_dot(ey, Rm0 + si * Pm0 - di);
2869  Rm1(3) = CLHEP_dot(ez, Rm0 + si * Pm0 - di);
2870  Pm1 = T * Pm0;
2871 
2872  newRl = Tl * (Rm1 - R0);
2873  newPl = Tl * Pm1;
2874 
2875  Vm(0) = newPl(1) / newPl(3);
2876  Vm(1) = newPl(2) / newPl(3);
2877  Vm(2) = newRl(1);
2878  Vm(3) = newRl(2);
2879  std::cout << " ========= dVm/da" << ii << " " << (Vm - Vm0) * (1. / del) << std::endl;
2880 #endif
2881 
2882  if (debug_) {
2883  //std::cout<<" dRm/da3 "<<((Rm1-Rm10)*(1./del)).T()<<" "<<((Pm1-Pm10)*(1./del)).T()<<std::endl;
2884  //std::cout<<" Norm "<<Norm.T()<<std::endl;
2885  std::cout << " ex " << (Tl.T() * ex0).T() << std::endl;
2886  std::cout << " ey " << (Tl.T() * ey0).T() << std::endl;
2887  std::cout << " ez " << (Tl.T() * ez0).T() << std::endl;
2888  //std::cpot<<"
2889 
2890  std::cout << " pxyz/p gl 0 " << (Pm0 / Pm0.norm()).T() << std::endl;
2891  std::cout << " pxyz/p loc0 " << (Tl * Pm0 / (Tl * Pm0).norm()).T() << std::endl;
2892  std::cout << " pxyz/p glob " << (Pm1 / Pm1.norm()).T() << std::endl;
2893  std::cout << " pxyz/p loc " << (newPl / newPl.norm()).T() << std::endl;
2894 
2895  //std::cout<<" Rot "<<refSurface.rotation()<<endl;
2896  //std::cout<<" Tl "<<Tl<<endl;
2897  std::cout << " ---- phi g0 g1 l0 l1 " << atan2(Pm0(2), Pm0(1)) << " " << atan2(Pm1(2), Pm1(1)) << " "
2898  << atan2((Tl * Pm0)(2), (Tl * Pm0)(1)) << " " << atan2(newPl(2), newPl(1)) << std::endl;
2899  std::cout << " ---- angl Gl Loc " << atan2(Pm1(2), Pm1(1)) - atan2(Pm0(2), Pm0(1)) << " "
2900  << atan2(newPl(2), newPl(1)) - atan2((Tl * Pm0)(2), (Tl * Pm0)(1)) << " "
2901  << atan2(newPl(3), newPl(2)) - atan2((Tl * Pm0)(3), (Tl * Pm0)(2)) << " "
2902  << atan2(newPl(1), newPl(3)) - atan2((Tl * Pm0)(1), (Tl * Pm0)(3)) << " " << std::endl;
2903  }
2904 
2905  if (debug_) {
2906  CLHEP::HepMatrix newTl(3, 3, 0);
2907  CLHEP::HepVector R0(3, 0), newRl(3, 0), newPl(3, 0);
2908  newTl = Tl * Ti.T();
2909  newR0 = Ti * R0 + di;
2910 
2911  std::cout << " N " << Norm.T() << std::endl;
2912  std::cout << " dtxl yl " << Vm(0) - tsosMuon.localParameters().vector()(1) << " "
2913  << Vm(1) - tsosMuon.localParameters().vector()(2) << std::endl;
2914  std::cout << " dXl dYl " << Vm(2) - tsosMuon.localParameters().vector()(3) << " "
2915  << Vm(3) - tsosMuon.localParameters().vector()(4) << std::endl;
2916  std::cout << " localM " << tsosMuon.localParameters().vector() << std::endl;
2917  std::cout << " Vm " << Vm << std::endl;
2918 
2919  CLHEP::HepVector Rvc(3, 0), Pvc(3, 0), Rvg(3, 0), Pvg(3, 0);
2920  Rvc(1) = Vm(2);
2921  Rvc(2) = Vm(3);
2922  Pvc(3) = tsosMuon.localParameters().pzSign() * Pm0.norm() / sqrt(1 + Vm(0) * Vm(0) + Vm(1) * Vm(1));
2923  Pvc(1) = Pvc(3) * Vm(0);
2924  Pvc(2) = Pvc(3) * Vm(1);
2925 
2926  Rvg = newTl.T() * Rvc + newR0;
2927  Pvg = newTl.T() * Pvc;
2928 
2929  std::cout << " newPl " << newPl.T() << std::endl;
2930  std::cout << " Pvc " << Pvc.T() << std::endl;
2931  std::cout << " Rvg " << Rvg.T() << std::endl;
2932  std::cout << " Rm1 " << Rm1.T() << std::endl;
2933  std::cout << " Pvg " << Pvg.T() << std::endl;
2934  std::cout << " Pm1 " << Pm1.T() << std::endl;
2935  }
2936 
2937  if (debug_) { // ---------- how to transform cartesian from shifted to correct
2938 
2939  std::cout << " ----- Pm " << Pm << std::endl;
2940  std::cout << " T*Pm0 " << Pm1.T() << std::endl;
2941  //std::cout<<" Ti*T*Pm0 "<<(Ti*(T*Pm0)).T()<<std::endl;
2942 
2943  //CLHEP::HepVector Rml = (Rm0 + si*Pm0 - di) + di;
2944  CLHEP::HepVector Rml = Rm1(1) * ex + Rm1(2) * ey + Rm1(3) * ez + di;
2945 
2946  CLHEP::HepVector Rt0(3);
2947  Rt0(1) = Rt.x();
2948  Rt0(2) = Rt.y();
2949  Rt0(3) = Rt.z();
2950 
2951  //double sl = CLHEP::dot(T*(Rt0 - Rml), T*Norm) / CLHEP::dot(Ti * Pm1, Norm);
2952  double sl = CLHEP_dot(T * (Rt0 - Rml), T * Norm) / CLHEP_dot(Ti * Pm1, Norm);
2953  CLHEP::HepVector Rl = Rml + sl * (Ti * Pm1);
2954 
2955  //double A = CLHEP::dot(Ti*Pm1, Norm);
2956  //double B = CLHEP::dot(Rt0, Norm) - CLHEP::dot(Ti*Rm1, Norm)
2957  //+ CLHEP::dot(Ti*d, Norm);
2958  double A = CLHEP_dot(Ti * Pm1, Norm);
2959  double B = CLHEP_dot(Rt0, Norm) - CLHEP_dot(Ti * Rm1, Norm) + CLHEP_dot(Ti * d, Norm);
2960  double s = B / A;
2961  CLHEP::HepVector Dr = Ti * Rm1 - Ti * d + s * (Ti * Pm1) - Rm0;
2962  CLHEP::HepVector Dp = Ti * Pm1 - Pm0;
2963 
2964  CLHEP::HepVector TiR = Ti * Rm0 + di;
2965  CLHEP::HepVector Ri = T * TiR + d;
2966 
2967  std::cout << " --TTi-0 " << (Ri - Rm0).T() << std::endl;
2968  std::cout << " -- Dr " << Dr.T() << endl;
2969  std::cout << " -- Dp " << Dp.T() << endl;
2970  //std::cout<<" ex "<<ex.T()<<endl;
2971  //std::cout<<" ey "<<ey.T()<<endl;
2972  //std::cout<<" ez "<<ez.T()<<endl;
2973  //std::cout<<" ---- T ---- "<<T<<std::endl;
2974  //std::cout<<" ---- Ti ---- "<<Ti<<std::endl;
2975  //std::cout<<" ---- d ---- "<<d.T()<<std::endl;
2976  //std::cout<<" ---- di ---- "<<di.T()<<std::endl;
2977  std::cout << " -- si sl s " << si << " " << sl << " " << s << std::endl;
2978  //std::cout<<" -- si sl "<<si<<" "<<sl<<endl;
2979  //std::cout<<" -- si s "<<si<<" "<<s<<endl;
2980  std::cout << " -- Rml-(Rm0+si*Pm0) " << (Rml - Rm0 - si * Pm0).T() << std::endl;
2981  //std::cout<<" -- Rm0 "<<Rm0.T()<<std::endl;
2982  //std::cout<<" -- Rm1 "<<Rm1.T()<<std::endl;
2983  //std::cout<<" -- Rmi "<<Rmi.T()<<std::endl;
2984  //std::cout<<" --siPm "<<(si*Pm0).T()<<std::endl;
2985  //std::cout<<" --s2Pm "<<(s2*(T * Pm1)).T()<<std::endl;
2986  //std::cout<<" --R1-0 "<<(Rm1-Rm0).T()<<std::endl;
2987  //std::cout<<" --Ri-0 "<<(Rmi-Rm0).T()<<std::endl;
2988  std::cout << " --Rl-0 " << (Rl - Rm0).T() << std::endl;
2989  //std::cout<<" --Pi-0 "<<(Pmi-Pm0).T()<<std::endl;
2990  } // -------- end of debug
2991 
2992  return;
2993 
2994 } // end misalignMuonL

References MaterialEffects_cfi::A, a, TtFullHadDaughter::B, alignmentValidation::c1, CLHEP_dot(), funct::cos(), gather_cfg::cout, ztail::d, debug_, cuy::ii, TrajectoryStateOnSurface::localParameters(), MuGlAngle, MuGlShift, GloballyPositioned< T >::position(), LocalTrajectoryParameters::pzSign(), HLT_FULL_cff::R0, GloballyPositioned< T >::rotation(), alignCSCRings::s, indexGen::s2, funct::sin(), mathSSE::sqrt(), TrajectoryStateOnSurface::surface(), LocalTrajectoryParameters::vector(), PV3DBase< T, PVType, FrameType >::x(), TkRotation< T >::xx(), TkRotation< T >::xy(), TkRotation< T >::xz(), PV3DBase< T, PVType, FrameType >::y(), TkRotation< T >::yx(), TkRotation< T >::yy(), TkRotation< T >::yz(), PV3DBase< T, PVType, FrameType >::z(), TkRotation< T >::zx(), TkRotation< T >::zy(), and TkRotation< T >::zz().

Referenced by analyzeTrackTrack(), and analyzeTrackTrajectory().

◆ muonFitter()

void GlobalTrackerMuonAlignment::muonFitter ( reco::TrackRef  alongTr,
reco::TransientTrack &  alongTTr,
PropagationDirection  direction,
TrajectoryStateOnSurface &  trackFittedTSOS 
)

Definition at line 3112 of file GlobalTrackerMuonAlignment.cc.

3115  {
3116  bool info = false;
3117  bool smooth = false;
3118 
3119  if (info)
3120  std::cout << " ......... start of muonFitter ........" << std::endl;
3121  if (false && info) {
3122  this->debugTrackHit(" Muon track hits ", alongTr);
3123  this->debugTrackHit(" Muon TransientTrack hits ", alongTTr);
3124  }
3125 
3126  edm::OwnVector<TrackingRecHit> recHit;
3127  DetId trackDetId = DetId(alongTr->innerDetId());
3128  for (auto const& hit : alongTr->recHits())
3129  recHit.push_back(hit->clone());
3130  if (direction != alongMomentum) {
3131  edm::OwnVector<TrackingRecHit> recHitAlong = recHit;
3132  recHit.clear();
3133  for (unsigned int ihit = recHitAlong.size() - 1; ihit + 1 > 0; ihit--) {
3134  recHit.push_back(recHitAlong[ihit]);
3135  }
3136  recHitAlong.clear();
3137  }
3138  if (info)
3139  std::cout << " ... Own recHit.size() DetId==Muon " << recHit.size() << " " << trackDetId.rawId() << std::endl;
3140 
3141  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu =
3142  MuRHBuilder->build(alongTTr.recHitsBegin(), alongTTr.recHitsEnd());
3143  if (info)
3144  std::cout << " ...along.... recHitMu.size() " << recHitMu.size() << std::endl;
3145  if (direction != alongMomentum) {
3146  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMuAlong = recHitMu;
3147  recHitMu.clear();
3148  for (unsigned int ihit = recHitMuAlong.size() - 1; ihit + 1 > 0; ihit--) {
3149  recHitMu.push_back(recHitMuAlong[ihit]);
3150  }
3151  recHitMuAlong.clear();
3152  }
3153  if (info)
3154  std::cout << " ...opposite... recHitMu.size() " << recHitMu.size() << std::endl;
3155 
3156  TrajectoryStateOnSurface firstTSOS;
3157  if (direction == alongMomentum)
3158  firstTSOS = alongTTr.innermostMeasurementState();
3159  else
3160  firstTSOS = alongTTr.outermostMeasurementState();
3161 
3162  AlgebraicSymMatrix55 CovLoc;
3163  CovLoc(0, 0) = 1. * firstTSOS.localError().matrix()(0, 0);
3164  CovLoc(1, 1) = 10. * firstTSOS.localError().matrix()(1, 1);
3165  CovLoc(2, 2) = 10. * firstTSOS.localError().matrix()(2, 2);
3166  CovLoc(3, 3) = 100. * firstTSOS.localError().matrix()(3, 3);
3167  CovLoc(4, 4) = 100. * firstTSOS.localError().matrix()(4, 4);
3168  TrajectoryStateOnSurface initialTSOS(
3169  firstTSOS.localParameters(), LocalTrajectoryError(CovLoc), firstTSOS.surface(), &*magneticField_);
3170 
3171  PTrajectoryStateOnDet PTraj = trajectoryStateTransform::persistentState(initialTSOS, trackDetId.rawId());
3172  //const TrajectorySeed seedT(*PTraj, recHit, alongMomentum);
3173  const TrajectorySeed seedT(PTraj, recHit, direction);
3174 
3175  std::vector<Trajectory> trajVec;
3176  if (direction == alongMomentum)
3177  trajVec = theFitter->fit(seedT, recHitMu, initialTSOS);
3178  else
3179  trajVec = theFitterOp->fit(seedT, recHitMu, initialTSOS);
3180 
3181  if (info) {
3182  this->debugTrajectorySOSv(" innermostTSOS of TransientTrack", alongTTr.innermostMeasurementState());
3183  this->debugTrajectorySOSv(" outermostTSOS of TransientTrack", alongTTr.outermostMeasurementState());
3184  this->debugTrajectorySOS(" initialTSOS for theFitter ", initialTSOS);
3185  std::cout << " . . . . . trajVec.size() " << trajVec.size() << std::endl;
3186  if (!trajVec.empty())
3187  this->debugTrajectory(" theFitter trajectory ", trajVec[0]);
3188  }
3189 
3190  if (!smooth) {
3191  if (!trajVec.empty())
3192  trackFittedTSOS = trajVec[0].lastMeasurement().updatedState();
3193  } else {
3194  std::vector<Trajectory> trajSVec;
3195  if (!trajVec.empty()) {
3196  if (direction == alongMomentum)
3197  trajSVec = theSmoother->trajectories(*(trajVec.begin()));
3198  else
3199  trajSVec = theSmootherOp->trajectories(*(trajVec.begin()));
3200  }
3201  if (info)
3202  std::cout << " . . . . trajSVec.size() " << trajSVec.size() << std::endl;
3203  if (!trajSVec.empty())
3204  this->debugTrajectorySOSv("smoothed geom InnermostState", trajSVec[0].geometricalInnermostState());
3205  if (!trajSVec.empty())
3206  trackFittedTSOS = trajSVec[0].geometricalInnermostState();
3207  }
3208 
3209 } // end of muonFitter

References alongMomentum, MuonTransientTrackingRecHitBuilder::build(), edm::OwnVector< T, P >::clear(), gather_cfg::cout, debugTrackHit(), debugTrajectory(), debugTrajectorySOS(), debugTrajectorySOSv(), info(), reco::TransientTrack::innermostMeasurementState(), TrajectoryStateOnSurface::localError(), TrajectoryStateOnSurface::localParameters(), magneticField_, LocalTrajectoryError::matrix(), MuRHBuilder, reco::TransientTrack::outermostMeasurementState(), trajectoryStateTransform::persistentState(), DetId::rawId(), rpcPointValidation_cfi::recHit, reco::TransientTrack::recHitsBegin(), reco::TransientTrack::recHitsEnd(), edm::OwnVector< T, P >::size(), TrajectoryStateOnSurface::surface(), theFitter, theFitterOp, theSmoother, theSmootherOp, and TrajectorySmoother::trajectories().

Referenced by analyzeTrackTrajectory().

◆ trackFitter()

void GlobalTrackerMuonAlignment::trackFitter ( reco::TrackRef  alongTr,
reco::TransientTrack &  alongTTr,
PropagationDirection  direction,
TrajectoryStateOnSurface &  trackFittedTSOS 
)

Definition at line 2997 of file GlobalTrackerMuonAlignment.cc.

3000  {
3001  bool info = false;
3002  bool smooth = false;
3003 
3004  if (info)
3005  std::cout << " ......... start of trackFitter ......... " << std::endl;
3006  if (false && info) {
3007  this->debugTrackHit(" Tracker track hits ", alongTr);
3008  this->debugTrackHit(" Tracker TransientTrack hits ", alongTTr);
3009  }
3010 
3011  edm::OwnVector<TrackingRecHit> recHit;
3012  DetId trackDetId = DetId(alongTr->innerDetId());
3013  for (auto const& hit : alongTr->recHits())
3014  recHit.push_back(hit->clone());
3015  if (direction != alongMomentum) {
3016  edm::OwnVector<TrackingRecHit> recHitAlong = recHit;
3017  recHit.clear();
3018  for (unsigned int ihit = recHitAlong.size() - 1; ihit + 1 > 0; ihit--) {
3019  recHit.push_back(recHitAlong[ihit]);
3020  }
3021  recHitAlong.clear();
3022  }
3023  if (info)
3024  std::cout << " ... Own recHit.size() " << recHit.size() << " " << trackDetId.rawId() << std::endl;
3025 
3026  //MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitTrack;
3027  TransientTrackingRecHit::RecHitContainer recHitTrack;
3028  for (auto const& hit : alongTr->recHits())
3029  recHitTrack.push_back(TTRHBuilder->build(hit));
3030 
3031  if (info)
3032  std::cout << " ... recHitTrack.size() " << recHit.size() << " " << trackDetId.rawId() << " ======> recHitMu "
3033  << std::endl;
3034 
3035  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu = recHitTrack;
3036  /*
3037  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu =
3038  MuRHBuilder->build(alongTTr.recHits().begin(), alongTTr.recHits().end());
3039  */
3040  if (info)
3041  std::cout << " ...along.... recHitMu.size() " << recHitMu.size() << std::endl;
3042  if (direction != alongMomentum) {
3043  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMuAlong = recHitMu;
3044  recHitMu.clear();
3045  for (unsigned int ihit = recHitMuAlong.size() - 1; ihit + 1 > 0; ihit--) {
3046  recHitMu.push_back(recHitMuAlong[ihit]);
3047  }
3048  recHitMuAlong.clear();
3049  }
3050  if (info)
3051  std::cout << " ...opposite... recHitMu.size() " << recHitMu.size() << std::endl;
3052 
3053  TrajectoryStateOnSurface firstTSOS;
3054  if (direction == alongMomentum)
3055  firstTSOS = alongTTr.innermostMeasurementState();
3056  else
3057  firstTSOS = alongTTr.outermostMeasurementState();
3058 
3059  AlgebraicSymMatrix55 CovLoc;
3060  CovLoc(0, 0) = 1. * firstTSOS.localError().matrix()(0, 0);
3061  CovLoc(1, 1) = 10. * firstTSOS.localError().matrix()(1, 1);
3062  CovLoc(2, 2) = 10. * firstTSOS.localError().matrix()(2, 2);
3063  CovLoc(3, 3) = 100. * firstTSOS.localError().matrix()(3, 3);
3064  CovLoc(4, 4) = 100. * firstTSOS.localError().matrix()(4, 4);
3065  TrajectoryStateOnSurface initialTSOS(
3066  firstTSOS.localParameters(), LocalTrajectoryError(CovLoc), firstTSOS.surface(), &*magneticField_);
3067 
3068  PTrajectoryStateOnDet PTraj = trajectoryStateTransform::persistentState(initialTSOS, trackDetId.rawId());
3069  //const TrajectorySeed seedT(*PTraj, recHit, alongMomentum);
3070  const TrajectorySeed seedT(PTraj, recHit, direction);
3071 
3072  std::vector<Trajectory> trajVec;
3073  if (direction == alongMomentum)
3074  trajVec = theFitter->fit(seedT, recHitMu, initialTSOS);
3075  else
3076  trajVec = theFitterOp->fit(seedT, recHitMu, initialTSOS);
3077 
3078  if (info) {
3079  this->debugTrajectorySOSv(" innermostTSOS of TransientTrack", alongTTr.innermostMeasurementState());
3080  this->debugTrajectorySOSv(" outermostTSOS of TransientTrack", alongTTr.outermostMeasurementState());
3081  this->debugTrajectorySOS(" initialTSOS for theFitter ", initialTSOS);
3082  std::cout << " . . . . . trajVec.size() " << trajVec.size() << std::endl;
3083  if (!trajVec.empty())
3084  this->debugTrajectory(" theFitter trajectory ", trajVec[0]);
3085  }
3086 
3087  if (!smooth) {
3088  if (!trajVec.empty())
3089  trackFittedTSOS = trajVec[0].lastMeasurement().updatedState();
3090  } else {
3091  std::vector<Trajectory> trajSVec;
3092  if (!trajVec.empty()) {
3093  if (direction == alongMomentum)
3094  trajSVec = theSmoother->trajectories(*(trajVec.begin()));
3095  else
3096  trajSVec = theSmootherOp->trajectories(*(trajVec.begin()));
3097  }
3098  if (info)
3099  std::cout << " . . . . trajSVec.size() " << trajSVec.size() << std::endl;
3100  if (!trajSVec.empty())
3101  this->debugTrajectorySOSv("smoothed geom InnermostState", trajSVec[0].geometricalInnermostState());
3102  if (!trajSVec.empty())
3103  trackFittedTSOS = trajSVec[0].geometricalInnermostState();
3104  }
3105 
3106  if (info)
3107  this->debugTrajectorySOSv(" trackFittedTSOS ", trackFittedTSOS);
3108 
3109 } // end of trackFitter

References alongMomentum, TransientTrackingRecHitBuilder::build(), edm::OwnVector< T, P >::clear(), gather_cfg::cout, debugTrackHit(), debugTrajectory(), debugTrajectorySOS(), debugTrajectorySOSv(), info(), reco::TransientTrack::innermostMeasurementState(), TrajectoryStateOnSurface::localError(), TrajectoryStateOnSurface::localParameters(), magneticField_, LocalTrajectoryError::matrix(), reco::TransientTrack::outermostMeasurementState(), trajectoryStateTransform::persistentState(), DetId::rawId(), rpcPointValidation_cfi::recHit, edm::OwnVector< T, P >::size(), TrajectoryStateOnSurface::surface(), theFitter, theFitterOp, theSmoother, theSmootherOp, TrajectorySmoother::trajectories(), and TTRHBuilder.

Referenced by analyzeTrackTrajectory().

◆ writeGlPosRcd()

void GlobalTrackerMuonAlignment::writeGlPosRcd ( CLHEP::HepVector &  d3)

Definition at line 3322 of file GlobalTrackerMuonAlignment.cc.

3322  {
3323  //paramVec(1) = 0.1; paramVec(2) = 0.2; paramVec(3) = 0.3; //0123
3324  //paramVec(4) = 0.0001; paramVec(5) = 0.0002; paramVec(6) = 0.0003;
3325  //paramVec(1) = 0.3; paramVec(2) = 0.4; paramVec(3) = 0.5; //03450123
3326  //paramVec(4) = 0.0001; paramVec(5) = 0.0002; paramVec(6) = 0.0003;
3327  //paramVec(1) = 2.; paramVec(2) = 2.; paramVec(3) = 2.; //222
3328  //paramVec(4) = 0.02; paramVec(5) = 0.02; paramVec(6) = 0.02;
3329  //paramVec(1) = -10.; paramVec(2) = -20.; paramVec(3) = -30.; //102030
3330  //paramVec(4) = -0.1; paramVec(5) = -0.2; paramVec(6) = -0.3;
3331  //paramVec(1) = 0.; paramVec(2) = 0.; paramVec(3) = 1.; //1z
3332  //paramVec(4) = 0.0000; paramVec(5) = 0.0000; paramVec(6) = 0.0000;
3333 
3334  std::cout << " paramVector " << paramVec.T() << std::endl;
3335 
3336  CLHEP::Hep3Vector colX, colY, colZ;
3337 
3338 #ifdef NOT_EXACT_ROTATION_MATRIX
3339  colX = CLHEP::Hep3Vector(1., -paramVec(6), paramVec(5));
3340  colY = CLHEP::Hep3Vector(paramVec(6), 1., -paramVec(4));
3341  colZ = CLHEP::Hep3Vector(-paramVec(5), paramVec(4), 1.);
3342 #else
3343  double s1 = std::sin(paramVec(4)), c1 = std::cos(paramVec(4));
3344  double s2 = std::sin(paramVec(5)), c2 = std::cos(paramVec(5));
3345  double s3 = std::sin(paramVec(6)), c3 = std::cos(paramVec(6));
3346  colX = CLHEP::Hep3Vector(c2 * c3, -c2 * s3, s2);
3347  colY = CLHEP::Hep3Vector(c1 * s3 + s1 * s2 * c3, c1 * c3 - s1 * s2 * s3, -s1 * c2);
3348  colZ = CLHEP::Hep3Vector(s1 * s3 - c1 * s2 * c3, s1 * c3 + c1 * s2 * s3, c1 * c2);
3349 #endif
3350 
3351  CLHEP::HepVector param0(6, 0);
3352 
3353  Alignments* globalPositions = new Alignments();
3354 
3355  //Tracker
3356  AlignTransform tracker(
3357  iteratorTrackerRcd->translation(), iteratorTrackerRcd->rotation(), DetId(DetId::Tracker).rawId());
3358  // Muon
3359  CLHEP::Hep3Vector posMuGlRcd = iteratorMuonRcd->translation();
3360  CLHEP::HepRotation rotMuGlRcd = iteratorMuonRcd->rotation();
3361  CLHEP::HepEulerAngles angMuGlRcd = iteratorMuonRcd->rotation().eulerAngles();
3362  if (debug_)
3363  std::cout << " Old muon Rcd Euler angles " << angMuGlRcd.phi() << " " << angMuGlRcd.theta() << " "
3364  << angMuGlRcd.psi() << std::endl;
3365  AlignTransform muon;
3366  if ((angMuGlRcd.phi() == 0.) && (angMuGlRcd.theta() == 0.) && (angMuGlRcd.psi() == 0.) && (posMuGlRcd.x() == 0.) &&
3367  (posMuGlRcd.y() == 0.) && (posMuGlRcd.z() == 0.)) {
3368  std::cout << " New muon parameters are stored in Rcd" << std::endl;
3369 
3370  AlignTransform muonNew(AlignTransform::Translation(paramVec(1), paramVec(2), paramVec(3)),
3371  AlignTransform::Rotation(colX, colY, colZ),
3372  DetId(DetId::Muon).rawId());
3373  muon = muonNew;
3374  } else if ((paramVec(1) == 0.) && (paramVec(2) == 0.) && (paramVec(3) == 0.) && (paramVec(4) == 0.) &&
3375  (paramVec(5) == 0.) && (paramVec(6) == 0.)) {
3376  std::cout << " Old muon parameters are stored in Rcd" << std::endl;
3377 
3378  AlignTransform muonNew(iteratorMuonRcd->translation(), iteratorMuonRcd->rotation(), DetId(DetId::Muon).rawId());
3379  muon = muonNew;
3380  } else {
3381  std::cout << " New + Old muon parameters are stored in Rcd" << std::endl;
3382  CLHEP::Hep3Vector posMuGlRcdThis = CLHEP::Hep3Vector(paramVec(1), paramVec(2), paramVec(3));
3383  CLHEP::HepRotation rotMuGlRcdThis = CLHEP::HepRotation(colX, colY, colZ);
3384  CLHEP::Hep3Vector posMuGlRcdNew = rotMuGlRcdThis * posMuGlRcd + posMuGlRcdThis;
3385  CLHEP::HepRotation rotMuGlRcdNew = rotMuGlRcdThis * rotMuGlRcd;
3386 
3387  AlignTransform muonNew(posMuGlRcdNew, rotMuGlRcdNew, DetId(DetId::Muon).rawId());
3388  muon = muonNew;
3389  }
3390 
3391  // Ecal
3392  AlignTransform ecal(iteratorEcalRcd->translation(), iteratorEcalRcd->rotation(), DetId(DetId::Ecal).rawId());
3393  // Hcal
3394  AlignTransform hcal(iteratorHcalRcd->translation(), iteratorHcalRcd->rotation(), DetId(DetId::Hcal).rawId());
3395  // Calo
3396  AlignTransform calo(AlignTransform::Translation(param0(1), param0(2), param0(3)),
3397  AlignTransform::EulerAngles(param0(4), param0(5), param0(6)),
3398  DetId(DetId::Calo).rawId());
3399 
3400  std::cout << "Tracker (" << tracker.rawId() << ") at " << tracker.translation() << " "
3401  << tracker.rotation().eulerAngles() << std::endl;
3402  std::cout << tracker.rotation() << std::endl;
3403 
3404  std::cout << "Muon (" << muon.rawId() << ") at " << muon.translation() << " " << muon.rotation().eulerAngles()
3405  << std::endl;
3406  std::cout << " rotations angles around x,y,z "
3407  << " ( " << -muon.rotation().zy() << " " << muon.rotation().zx() << " " << -muon.rotation().yx() << " )"
3408  << std::endl;
3409  std::cout << muon.rotation() << std::endl;
3410 
3411  std::cout << "Ecal (" << ecal.rawId() << ") at " << ecal.translation() << " " << ecal.rotation().eulerAngles()
3412  << std::endl;
3413  std::cout << ecal.rotation() << std::endl;
3414 
3415  std::cout << "Hcal (" << hcal.rawId() << ") at " << hcal.translation() << " " << hcal.rotation().eulerAngles()
3416  << std::endl;
3417  std::cout << hcal.rotation() << std::endl;
3418 
3419  std::cout << "Calo (" << calo.rawId() << ") at " << calo.translation() << " " << calo.rotation().eulerAngles()
3420  << std::endl;
3421  std::cout << calo.rotation() << std::endl;
3422 
3423  globalPositions->m_align.push_back(tracker);
3424  globalPositions->m_align.push_back(muon);
3425  globalPositions->m_align.push_back(ecal);
3426  globalPositions->m_align.push_back(hcal);
3427  globalPositions->m_align.push_back(calo);
3428 
3429  std::cout << "Uploading to the database..." << std::endl;
3430 
3431  edm::Service<cond::service::PoolDBOutputService> poolDbService;
3432 
3433  if (!poolDbService.isAvailable())
3434  throw cms::Exception("NotAvailable") << "PoolDBOutputService not available";
3435 
3436  // if (poolDbService->isNewTagRequest("GlobalPositionRcd")) {
3437  // poolDbService->createNewIOV<Alignments>(&(*globalPositions), poolDbService->endOfTime(), "GlobalPositionRcd");
3438  // } else {
3439  // poolDbService->appendSinceTime<Alignments>(&(*globalPositions), poolDbService->currentTime(), "GlobalPositionRcd");
3440  // }
3441  poolDbService->writeOne<Alignments>(&(*globalPositions),
3442  poolDbService->currentTime(),
3443  //poolDbService->beginOfTime(),
3444  "GlobalPositionRcd");
3445  std::cout << "done!" << std::endl;
3446 
3447  return;
3448 }

References alignmentValidation::c1, DetId::Calo, L1TowerCalibrationProducer_cfi::calo, funct::cos(), gather_cfg::cout, cond::service::PoolDBOutputService::currentTime(), debug_, bsc_activity_cfg::ecal, DetId::Ecal, DetId::Hcal, patCandidatesForDimuonsSequences_cff::hcal, edm::Service< T >::isAvailable(), iteratorEcalRcd, iteratorHcalRcd, iteratorMuonRcd, iteratorTrackerRcd, Alignments::m_align, DetId::Muon, HLT_FULL_cff::muon, DetId::rawId(), indexGen::s2, funct::sin(), DetId::Tracker, PbPb_ZMuSkimMuonDPG_cff::tracker, and cond::service::PoolDBOutputService::writeOne().

Referenced by endJob().

Member Data Documentation

◆ collectionCosmic

bool GlobalTrackerMuonAlignment::collectionCosmic
private

Definition at line 166 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), beginJob(), and endJob().

◆ collectionIsolated

bool GlobalTrackerMuonAlignment::collectionIsolated
private

Definition at line 167 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), beginJob(), and endJob().

◆ cosmicMuonMode_

bool GlobalTrackerMuonAlignment::cosmicMuonMode_
private

Definition at line 164 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), beginJob(), and GlobalTrackerMuonAlignment().

◆ debug_

bool GlobalTrackerMuonAlignment::debug_
private

Definition at line 173 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), endJob(), gradientGlobal(), gradientGlobalAlg(), gradientLocal(), misalignMuon(), misalignMuonL(), and writeGlPosRcd().

◆ defineFitter

bool GlobalTrackerMuonAlignment::defineFitter
private

Definition at line 191 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrajectory().

◆ file

TFile* GlobalTrackerMuonAlignment::file
private

Definition at line 221 of file GlobalTrackerMuonAlignment.cc.

Referenced by beginJob(), python.rootplot.utilities.RootFile::cd(), endJob(), ztee.GZipLog::finish(), ztee.GZipLog::flush_block(), python.rootplot.utilities.RootFile::get(), tfile.TFileService::stop(), core.VertexHistograms.VertexHistograms::write(), core.TreeAnalyzerNumpy.TreeAnalyzerNumpy::write(), and ztee.GZipLog::write_block().

◆ Gfr

AlgebraicVector3 GlobalTrackerMuonAlignment::Gfr
private

Definition at line 196 of file GlobalTrackerMuonAlignment.cc.

Referenced by beginJob(), endJob(), and gradientGlobalAlg().

◆ globalPositionRcd_

const Alignments* GlobalTrackerMuonAlignment::globalPositionRcd_
private

Definition at line 185 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), and analyzeTrackTrajectory().

◆ gmuonTags_

edm::InputTag GlobalTrackerMuonAlignment::gmuonTags_
private

Definition at line 161 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), and analyzeTrackTrajectory().

◆ Grad

CLHEP::HepVector GlobalTrackerMuonAlignment::Grad
private

Definition at line 199 of file GlobalTrackerMuonAlignment.cc.

Referenced by beginJob(), endJob(), and gradientGlobal().

◆ GradL

CLHEP::HepVector GlobalTrackerMuonAlignment::GradL
private

Definition at line 202 of file GlobalTrackerMuonAlignment.cc.

Referenced by beginJob(), endJob(), and gradientLocal().

◆ Hess

CLHEP::HepMatrix GlobalTrackerMuonAlignment::Hess
private

Definition at line 200 of file GlobalTrackerMuonAlignment.cc.

Referenced by beginJob(), endJob(), and gradientGlobal().

◆ HessL

CLHEP::HepMatrix GlobalTrackerMuonAlignment::HessL
private

Definition at line 203 of file GlobalTrackerMuonAlignment.cc.

Referenced by beginJob(), endJob(), and gradientLocal().

◆ histo

TH1F* GlobalTrackerMuonAlignment::histo
private

Definition at line 222 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), and bookHist().

◆ histo101

TH2F* GlobalTrackerMuonAlignment::histo101
private

Definition at line 257 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), and bookHist().

◆ histo102

TH2F* GlobalTrackerMuonAlignment::histo102
private

Definition at line 258 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), and bookHist().

◆ histo11

TH1F* GlobalTrackerMuonAlignment::histo11
private

Definition at line 231 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), and bookHist().

◆ histo12

TH1F* GlobalTrackerMuonAlignment::histo12
private

Definition at line 232 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), bookHist(), and fitHist().

◆ histo13

TH1F* GlobalTrackerMuonAlignment::histo13
private

Definition at line 233 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), bookHist(), and fitHist().

◆ histo14

TH1F* GlobalTrackerMuonAlignment::histo14
private

Definition at line 234 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), bookHist(), and fitHist().

◆ histo15

TH1F* GlobalTrackerMuonAlignment::histo15
private

Definition at line 235 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), bookHist(), and fitHist().

◆ histo16

TH1F* GlobalTrackerMuonAlignment::histo16
private

Definition at line 236 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), bookHist(), and fitHist().

◆ histo17

TH1F* GlobalTrackerMuonAlignment::histo17
private

Definition at line 237 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), bookHist(), and fitHist().

◆ histo18

TH1F* GlobalTrackerMuonAlignment::histo18
private

Definition at line 238 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), and bookHist().

◆ histo19

TH1F* GlobalTrackerMuonAlignment::histo19
private

Definition at line 239 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), and bookHist().

◆ histo2

TH1F* GlobalTrackerMuonAlignment::histo2
private

Definition at line 223 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), and bookHist().

◆ histo20

TH1F* GlobalTrackerMuonAlignment::histo20
private

Definition at line 240 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), and bookHist().

◆ histo21

TH1F* GlobalTrackerMuonAlignment::histo21
private

Definition at line 241 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), bookHist(), and fitHist().

◆ histo22

TH1F* GlobalTrackerMuonAlignment::histo22
private

Definition at line 242 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), bookHist(), and fitHist().

◆ histo23

TH1F* GlobalTrackerMuonAlignment::histo23
private

Definition at line 243 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), bookHist(), and fitHist().

◆ histo24

TH1F* GlobalTrackerMuonAlignment::histo24
private

Definition at line 244 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), bookHist(), and fitHist().

◆ histo25

TH1F* GlobalTrackerMuonAlignment::histo25
private

Definition at line 245 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), bookHist(), and fitHist().

◆ histo26

TH1F* GlobalTrackerMuonAlignment::histo26
private

Definition at line 246 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), bookHist(), and fitHist().

◆ histo27

TH1F* GlobalTrackerMuonAlignment::histo27
private

Definition at line 247 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), and bookHist().

◆ histo28

TH1F* GlobalTrackerMuonAlignment::histo28
private

Definition at line 248 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), and bookHist().

◆ histo29

TH1F* GlobalTrackerMuonAlignment::histo29
private

Definition at line 249 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), and bookHist().

◆ histo3

TH1F* GlobalTrackerMuonAlignment::histo3
private

Definition at line 224 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), and bookHist().

◆ histo30

TH1F* GlobalTrackerMuonAlignment::histo30
private

Definition at line 250 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), and bookHist().

◆ histo31

TH1F* GlobalTrackerMuonAlignment::histo31
private

Definition at line 251 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), and bookHist().

◆ histo32

TH1F* GlobalTrackerMuonAlignment::histo32
private

Definition at line 252 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), bookHist(), and fitHist().

◆ histo33

TH1F* GlobalTrackerMuonAlignment::histo33
private

Definition at line 253 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), bookHist(), and fitHist().

◆ histo34

TH1F* GlobalTrackerMuonAlignment::histo34
private

Definition at line 254 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), bookHist(), and fitHist().

◆ histo35

TH1F* GlobalTrackerMuonAlignment::histo35
private

Definition at line 255 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), bookHist(), and fitHist().

◆ histo4

TH1F* GlobalTrackerMuonAlignment::histo4
private

Definition at line 225 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), and bookHist().

◆ histo5

TH1F* GlobalTrackerMuonAlignment::histo5
private

Definition at line 226 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), and bookHist().

◆ histo6

TH1F* GlobalTrackerMuonAlignment::histo6
private

Definition at line 227 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), and bookHist().

◆ histo7

TH1F* GlobalTrackerMuonAlignment::histo7
private

Definition at line 228 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), bookHist(), and fitHist().

◆ histo8

TH1F* GlobalTrackerMuonAlignment::histo8
private

Definition at line 229 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), and bookHist().

◆ Inf

AlgebraicSymMatrix33 GlobalTrackerMuonAlignment::Inf
private

Definition at line 197 of file GlobalTrackerMuonAlignment.cc.

Referenced by beginJob(), endJob(), and gradientGlobalAlg().

◆ isolatedMuonMode_

bool GlobalTrackerMuonAlignment::isolatedMuonMode_
private

Definition at line 165 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), beginJob(), and GlobalTrackerMuonAlignment().

◆ iteratorEcalRcd

std::vector<AlignTransform>::const_iterator GlobalTrackerMuonAlignment::iteratorEcalRcd
private

Definition at line 211 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), and writeGlPosRcd().

◆ iteratorHcalRcd

std::vector<AlignTransform>::const_iterator GlobalTrackerMuonAlignment::iteratorHcalRcd
private

Definition at line 212 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), and writeGlPosRcd().

◆ iteratorMuonRcd

std::vector<AlignTransform>::const_iterator GlobalTrackerMuonAlignment::iteratorMuonRcd
private

Definition at line 210 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), and writeGlPosRcd().

◆ iteratorTrackerRcd

std::vector<AlignTransform>::const_iterator GlobalTrackerMuonAlignment::iteratorTrackerRcd
private

Definition at line 209 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), and writeGlPosRcd().

◆ magneticField_

const MagneticField* GlobalTrackerMuonAlignment::magneticField_
private

Definition at line 180 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), muonFitter(), and trackFitter().

◆ MuGlAngle

CLHEP::HepVector GlobalTrackerMuonAlignment::MuGlAngle
private

Definition at line 215 of file GlobalTrackerMuonAlignment.cc.

Referenced by endJob(), misalignMuon(), and misalignMuonL().

◆ MuGlShift

CLHEP::HepVector GlobalTrackerMuonAlignment::MuGlShift
private

Definition at line 214 of file GlobalTrackerMuonAlignment.cc.

Referenced by endJob(), misalignMuon(), and misalignMuonL().

◆ muonTags_

edm::InputTag GlobalTrackerMuonAlignment::muonTags_
private

Definition at line 160 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), and analyzeTrackTrajectory().

◆ MuRHBuilder

MuonTransientTrackingRecHitBuilder* GlobalTrackerMuonAlignment::MuRHBuilder
private

Definition at line 192 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrajectory(), and muonFitter().

◆ MuSelect

std::string GlobalTrackerMuonAlignment::MuSelect
private

Definition at line 217 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), and endJob().

◆ N_event

int GlobalTrackerMuonAlignment::N_event
private

Definition at line 205 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), beginJob(), and endJob().

◆ N_track

int GlobalTrackerMuonAlignment::N_track
private

Definition at line 206 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), beginJob(), and endJob().

◆ OutGlobalTxt

ofstream GlobalTrackerMuonAlignment::OutGlobalTxt
private

Definition at line 219 of file GlobalTrackerMuonAlignment.cc.

Referenced by endJob().

◆ propagator_

string GlobalTrackerMuonAlignment::propagator_
private

Definition at line 163 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), and analyzeTrackTrajectory().

◆ refitMuon_

bool GlobalTrackerMuonAlignment::refitMuon_
private

Definition at line 168 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrajectory().

◆ refitTrack_

bool GlobalTrackerMuonAlignment::refitTrack_
private

Definition at line 169 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrajectory().

◆ rootOutFile_

string GlobalTrackerMuonAlignment::rootOutFile_
private

Definition at line 170 of file GlobalTrackerMuonAlignment.cc.

Referenced by beginJob().

◆ smuonTags_

edm::InputTag GlobalTrackerMuonAlignment::smuonTags_
private

Definition at line 162 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), analyzeTrackTrajectory(), and endJob().

◆ theFitter

KFTrajectoryFitter* GlobalTrackerMuonAlignment::theFitter
private

Definition at line 187 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrajectory(), muonFitter(), and trackFitter().

◆ theFitterOp

KFTrajectoryFitter* GlobalTrackerMuonAlignment::theFitterOp
private

Definition at line 189 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrajectory(), muonFitter(), and trackFitter().

◆ theSmoother

KFTrajectorySmoother* GlobalTrackerMuonAlignment::theSmoother
private

Definition at line 188 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrajectory(), muonFitter(), and trackFitter().

◆ theSmootherOp

KFTrajectorySmoother* GlobalTrackerMuonAlignment::theSmootherOp
private

Definition at line 190 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrajectory(), muonFitter(), and trackFitter().

◆ theTrackingGeometry

edm::ESHandle<GlobalTrackingGeometry> GlobalTrackerMuonAlignment::theTrackingGeometry
private

Definition at line 176 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrajectory().

◆ trackingGeometry_

const GlobalTrackingGeometry* GlobalTrackerMuonAlignment::trackingGeometry_
private

Definition at line 177 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), and analyzeTrackTrajectory().

◆ trackTags_

edm::InputTag GlobalTrackerMuonAlignment::trackTags_
private

Definition at line 159 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), and analyzeTrackTrajectory().

◆ TTRHBuilder

const TransientTrackingRecHitBuilder* GlobalTrackerMuonAlignment::TTRHBuilder
private

Definition at line 193 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrajectory(), and trackFitter().

◆ txtOutFile_

string GlobalTrackerMuonAlignment::txtOutFile_
private

Definition at line 171 of file GlobalTrackerMuonAlignment.cc.

Referenced by endJob().

◆ watchGlobalPositionRcd_

edm::ESWatcher<GlobalPositionRcd> GlobalTrackerMuonAlignment::watchGlobalPositionRcd_
private

Definition at line 184 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), and analyzeTrackTrajectory().

◆ watchMagneticFieldRecord_

edm::ESWatcher<IdealMagneticFieldRecord> GlobalTrackerMuonAlignment::watchMagneticFieldRecord_
private

Definition at line 179 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), and analyzeTrackTrajectory().

◆ watchTrackingComponentsRecord_

edm::ESWatcher<TrackingComponentsRecord> GlobalTrackerMuonAlignment::watchTrackingComponentsRecord_
private

Definition at line 182 of file GlobalTrackerMuonAlignment.cc.

◆ watchTrackingGeometry_

edm::ESWatcher<GlobalTrackingGeometryRecord> GlobalTrackerMuonAlignment::watchTrackingGeometry_
private

Definition at line 175 of file GlobalTrackerMuonAlignment.cc.

Referenced by analyzeTrackTrack(), and analyzeTrackTrajectory().

◆ writeDB_

bool GlobalTrackerMuonAlignment::writeDB_
private

Definition at line 172 of file GlobalTrackerMuonAlignment.cc.

Referenced by endJob().

Vector3DBase
Definition: Vector3DBase.h:8
GlobalTrackerMuonAlignment::magneticField_
const MagneticField * magneticField_
Definition: GlobalTrackerMuonAlignment.cc:180
edm::ESWatcher::check
bool check(const edm::EventSetup &iSetup)
Definition: ESWatcher.h:52
PI
Definition: PayloadInspector.h:20
cms::cuda::V
cudaStream_t T uint32_t const T *__restrict__ const uint32_t *__restrict__ uint32_t int cudaStream_t V
Definition: HistoContainer.h:99
PDWG_BPHSkim_cff.muons
muons
Definition: PDWG_BPHSkim_cff.py:47
DDAxes::y
AlgebraicVector3
ROOT::Math::SVector< double, 3 > AlgebraicVector3
Definition: AlgebraicROOTObjects.h:12
PDWG_EXOHSCP_cff.tracks
tracks
Definition: PDWG_EXOHSCP_cff.py:28
GlobalTrackerMuonAlignment::analyzeTrackTrajectory
void analyzeTrackTrajectory(const edm::Event &, const edm::EventSetup &)
Definition: GlobalTrackerMuonAlignment.cc:1252
AlgebraicSymMatrix33
ROOT::Math::SMatrix< double, 3, 3, ROOT::Math::MatRepSym< double, 3 > > AlgebraicSymMatrix33
Definition: AlgebraicROOTObjects.h:21
GlobalTrackerMuonAlignment::histo33
TH1F * histo33
Definition: GlobalTrackerMuonAlignment.cc:253
runGCPTkAlMap.title
string title
Definition: runGCPTkAlMap.py:94
mps_fire.i
i
Definition: mps_fire.py:428
GlobalTrackerMuonAlignment::histo18
TH1F * histo18
Definition: GlobalTrackerMuonAlignment.cc:238
TkRotation::zy
T zy() const
Definition: extTkRotation.h:258
GlobalTrackerMuonAlignment::histo7
TH1F * histo7
Definition: GlobalTrackerMuonAlignment.cc:228
defaultRKPropagator::Product
Definition: defaultRKPropagator.h:41
GlobalTrackerMuonAlignment::theTrackingGeometry
edm::ESHandle< GlobalTrackingGeometry > theTrackingGeometry
Definition: GlobalTrackerMuonAlignment.cc:176
L1TowerCalibrationProducer_cfi.calo
calo
Definition: L1TowerCalibrationProducer_cfi.py:59
GlobalTrackerMuonAlignment::trackTags_
edm::InputTag trackTags_
Definition: GlobalTrackerMuonAlignment.cc:159
GlobalTrackerMuonAlignment::collectionIsolated
bool collectionIsolated
Definition: GlobalTrackerMuonAlignment.cc:167
TrajectoryStateOnSurface::cartesianError
const CartesianTrajectoryError cartesianError() const
Definition: TrajectoryStateOnSurface.h:71
GlobalTrackerMuonAlignment::iteratorTrackerRcd
std::vector< AlignTransform >::const_iterator iteratorTrackerRcd
Definition: GlobalTrackerMuonAlignment.cc:209
muon
Definition: MuonCocktails.h:17
Trajectory::chiSquared
float chiSquared() const
Definition: Trajectory.h:241
PV3DBase::x
T x() const
Definition: PV3DBase.h:59
Trajectory::direction
PropagationDirection const & direction() const
Definition: Trajectory.cc:133
GlobalTrackerMuonAlignment::histo32
TH1F * histo32
Definition: GlobalTrackerMuonAlignment.cc:252
GlobalTrackerMuonAlignment::Inf
AlgebraicSymMatrix33 Inf
Definition: GlobalTrackerMuonAlignment.cc:197
hcal
Definition: ConfigurationDatabase.cc:13
edm
HLT enums.
Definition: AlignableModifier.h:19
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
trajectoryStateTransform::persistentState
PTrajectoryStateOnDet persistentState(const TrajectoryStateOnSurface &ts, unsigned int detid)
Definition: TrajectoryStateTransform.cc:14
GlobalTrackerMuonAlignment::histo34
TH1F * histo34
Definition: GlobalTrackerMuonAlignment.cc:254
GlobalTrackerMuonAlignment::histo15
TH1F * histo15
Definition: GlobalTrackerMuonAlignment.cc:235
GlobalTrackerMuonAlignment::Gfr
AlgebraicVector3 Gfr
Definition: GlobalTrackerMuonAlignment.cc:196
AlignTransform::Rotation
CLHEP::HepRotation Rotation
Definition: AlignTransform.h:19
TrajectoryStateOnSurface::globalPosition
GlobalPoint globalPosition() const
Definition: TrajectoryStateOnSurface.h:65
MuonTransientTrackingRecHitBuilder::ConstRecHitContainer
TransientTrackingRecHit::ConstRecHitContainer ConstRecHitContainer
Definition: MuonTransientTrackingRecHitBuilder.h:15
gather_cfg.cout
cout
Definition: gather_cfg.py:144
GlobalTrackerMuonAlignment::theSmoother
KFTrajectorySmoother * theSmoother
Definition: GlobalTrackerMuonAlignment.cc:188
DetId::Hcal
Definition: DetId.h:28
GlobalTrackerMuonAlignment::MuGlAngle
CLHEP::HepVector MuGlAngle
Definition: GlobalTrackerMuonAlignment.cc:215
GlobalTrackerMuonAlignment::misalignMuonL
void misalignMuonL(GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, AlgebraicVector4 &, TrajectoryStateOnSurface &, TrajectoryStateOnSurface &)
Definition: GlobalTrackerMuonAlignment.cc:2715
GlobalTrackerMuonAlignment::histo26
TH1F * histo26
Definition: GlobalTrackerMuonAlignment.cc:246
TransientRecHitRecord
Definition: TransientRecHitRecord.h:14
TrajectoryStateOnSurface::charge
TrackCharge charge() const
Definition: TrajectoryStateOnSurface.h:68
GlobalTrackerMuonAlignment::histo28
TH1F * histo28
Definition: GlobalTrackerMuonAlignment.cc:248
GlobalTrackerMuonAlignment::histo24
TH1F * histo24
Definition: GlobalTrackerMuonAlignment.cc:244
GlobalTrackerMuonAlignment::Grad
CLHEP::HepVector Grad
Definition: GlobalTrackerMuonAlignment.cc:199
GlobalTrackerMuonAlignment::histo11
TH1F * histo11
Definition: GlobalTrackerMuonAlignment.cc:231
GlobalTrackerMuonAlignment::txtOutFile_
string txtOutFile_
Definition: GlobalTrackerMuonAlignment.cc:171
oppositeToMomentum
Definition: PropagationDirection.h:4
AlignTransform::EulerAngles
CLHEP::HepEulerAngles EulerAngles
Definition: AlignTransform.h:17
ConstReferenceCountingPointer
Definition: ReferenceCounted.h:67
GlobalTrackerMuonAlignment::histo13
TH1F * histo13
Definition: GlobalTrackerMuonAlignment.cc:233
AlignTransform::Translation
CLHEP::Hep3Vector Translation
Definition: AlignTransform.h:18
GlobalTrackerMuonAlignment::watchTrackingGeometry_
edm::ESWatcher< GlobalTrackingGeometryRecord > watchTrackingGeometry_
Definition: GlobalTrackerMuonAlignment.cc:175
GlobalTrackingGeometryRecord
Definition: GlobalTrackingGeometryRecord.h:17
GlobalTrackerMuonAlignment::trackingGeometry_
const GlobalTrackingGeometry * trackingGeometry_
Definition: GlobalTrackerMuonAlignment.cc:177
Surface
Definition: Surface.h:36
info
static const TGPicture * info(bool iBackgroundIsBlack)
Definition: FWCollectionSummaryWidget.cc:153
SteppingHelixPropagator_cfi.SteppingHelixPropagator
SteppingHelixPropagator
Definition: SteppingHelixPropagator_cfi.py:3
AlgebraicMatrix33
ROOT::Math::SMatrix< double, 3, 3, ROOT::Math::MatRepStd< double, 3, 3 > > AlgebraicMatrix33
Definition: AlgebraicROOTObjects.h:41
TrajectoryMeasurement::updatedState
TrajectoryStateOnSurface const & updatedState() const
Definition: TrajectoryMeasurement.h:184
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
LocalTrajectoryError::matrix
const AlgebraicSymMatrix55 & matrix() const
Definition: LocalTrajectoryError.h:60
DDAxes::x
GlobalTrackerMuonAlignment::muonFitter
void muonFitter(reco::TrackRef, reco::TransientTrack &, PropagationDirection, TrajectoryStateOnSurface &)
Definition: GlobalTrackerMuonAlignment.cc:3112
indexGen.s2
s2
Definition: indexGen.py:107
GlobalTrackerMuonAlignment::writeDB_
bool writeDB_
Definition: GlobalTrackerMuonAlignment.cc:172
GlobalTrackerMuonAlignment::histo35
TH1F * histo35
Definition: GlobalTrackerMuonAlignment.cc:255
GlobalTrackerMuonAlignment::histo16
TH1F * histo16
Definition: GlobalTrackerMuonAlignment.cc:236
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:45
Trajectory::foundHits
int foundHits() const
Definition: Trajectory.h:206
GlobalVector
Global3DVector GlobalVector
Definition: GlobalVector.h:10
HLT_FULL_cff.magneticField
magneticField
Definition: HLT_FULL_cff.py:348
reco::TransientTrack::recHitsBegin
trackingRecHit_iterator recHitsBegin() const
first iterator to RecHits
Definition: TransientTrack.h:126
TransientTrackingRecHitBuilder::build
virtual RecHitPointer build(const TrackingRecHit *p) const =0
build a tracking rechit from an existing rechit
GlobalTrackerMuonAlignment::TTRHBuilder
const TransientTrackingRecHitBuilder * TTRHBuilder
Definition: GlobalTrackerMuonAlignment.cc:193
edm::Handle< reco::TrackCollection >
GlobalTrackerMuonAlignment::histo2
TH1F * histo2
Definition: GlobalTrackerMuonAlignment.cc:223
DetId::Calo
Definition: DetId.h:29
edm::Service::isAvailable
bool isAvailable() const
Definition: Service.h:40
GlobalTrackerMuonAlignment::histo30
TH1F * histo30
Definition: GlobalTrackerMuonAlignment.cc:250
GlobalPositionRcd
Definition: GlobalPositionRcd.h:6
rpcPointValidation_cfi.recHit
recHit
Definition: rpcPointValidation_cfi.py:7
GlobalTrackerMuonAlignment::N_event
int N_event
Definition: GlobalTrackerMuonAlignment.cc:205
GlobalTrackerMuonAlignment::histo6
TH1F * histo6
Definition: GlobalTrackerMuonAlignment.cc:227
geometryDiff.epsilon
int epsilon
Definition: geometryDiff.py:26
TrajectorySmoother::trajectories
virtual TrajectoryContainer trajectories(const Trajectory &traj) const
Definition: TrajectorySmoother.h:18
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
GlobalTrackerMuonAlignment::theSmootherOp
KFTrajectorySmoother * theSmootherOp
Definition: GlobalTrackerMuonAlignment.cc:190
PV3DBase::z
T z() const
Definition: PV3DBase.h:61
patCandidatesForDimuonsSequences_cff.hcal
hcal
Definition: patCandidatesForDimuonsSequences_cff.py:37
IdealMagneticFieldRecord
Definition: IdealMagneticFieldRecord.h:11
DetId
Definition: DetId.h:17
GlobalTrackerMuonAlignment::cosmicMuonMode_
bool cosmicMuonMode_
Definition: GlobalTrackerMuonAlignment.cc:164
TrajectoryStateOnSurface
Definition: TrajectoryStateOnSurface.h:16
alignCSCRings.s
s
Definition: alignCSCRings.py:92
GlobalTrackerMuonAlignment::histo101
TH2F * histo101
Definition: GlobalTrackerMuonAlignment.cc:257
SteppingHelixPropagator
Definition: SteppingHelixPropagator.h:36
GlobalTrackerMuonAlignment::refitTrack_
bool refitTrack_
Definition: GlobalTrackerMuonAlignment.cc:169
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
TrajectoryStateOnSurface::freeState
FreeTrajectoryState const * freeState(bool withErrors=true) const
Definition: TrajectoryStateOnSurface.h:58
edm::EventSetup::get
T get() const
Definition: EventSetup.h:80
TrackCandidateProducer_cfi.propagator
propagator
Definition: TrackCandidateProducer_cfi.py:17
GlobalTrackerMuonAlignment::collectionCosmic
bool collectionCosmic
Definition: GlobalTrackerMuonAlignment.cc:166
GlobalTrackerMuonAlignment::N_track
int N_track
Definition: GlobalTrackerMuonAlignment.cc:206
TrackingRecHit::RecHitContainer
std::vector< ConstRecHitPointer > RecHitContainer
Definition: TrackingRecHit.h:31
GlobalTrackerMuonAlignment::isolatedMuonMode_
bool isolatedMuonMode_
Definition: GlobalTrackerMuonAlignment.cc:165
GlobalTrackerMuonAlignment::histo14
TH1F * histo14
Definition: GlobalTrackerMuonAlignment.cc:234
HLT_FULL_cff.muon
muon
Definition: HLT_FULL_cff.py:11710
GlobalTrackerMuonAlignment::histo21
TH1F * histo21
Definition: GlobalTrackerMuonAlignment.cc:241
TkRotation::xz
T xz() const
Definition: extTkRotation.h:253
reco::TransientTrack::recHitsEnd
trackingRecHit_iterator recHitsEnd() const
last iterator to RecHits
Definition: TransientTrack.h:128
w
const double w
Definition: UKUtility.cc:23
visualization-live-secondInstance_cfg.m
m
Definition: visualization-live-secondInstance_cfg.py:72
GlobalTrackerMuonAlignment::histo19
TH1F * histo19
Definition: GlobalTrackerMuonAlignment.cc:239
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
GlobalTrackerMuonAlignment::smuonTags_
edm::InputTag smuonTags_
Definition: GlobalTrackerMuonAlignment.cc:162
AlgebraicSymMatrix66
ROOT::Math::SMatrix< double, 6, 6, ROOT::Math::MatRepSym< double, 6 > > AlgebraicSymMatrix66
Definition: AlgebraicROOTObjects.h:24
DDAxes::z
Chi2MeasurementEstimator_cfi.Chi2MeasurementEstimator
Chi2MeasurementEstimator
Definition: Chi2MeasurementEstimator_cfi.py:5
GlobalTrackerMuonAlignment::histo8
TH1F * histo8
Definition: GlobalTrackerMuonAlignment.cc:229
GlobalTrackerMuonAlignment::bookHist
void bookHist()
Definition: GlobalTrackerMuonAlignment.cc:471
edm::ESHandle< GlobalTrackingGeometry >
GlobalTrackerMuonAlignment::gradientGlobalAlg
void gradientGlobalAlg(GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, AlgebraicSymMatrix66 &)
Definition: GlobalTrackerMuonAlignment.cc:2032
GlobalTrajectoryParameters::vector
AlgebraicVector6 vector() const
Definition: GlobalTrajectoryParameters.h:97
GlobalTrackerMuonAlignment::rootOutFile_
string rootOutFile_
Definition: GlobalTrackerMuonAlignment.cc:170
LocalTrajectoryParameters::vector
AlgebraicVector5 vector() const
Definition: LocalTrajectoryParameters.h:120
GlobalTrackerMuonAlignment::OutGlobalTxt
ofstream OutGlobalTxt
Definition: GlobalTrackerMuonAlignment.cc:219
dqmdumpme.k
k
Definition: dqmdumpme.py:60
Point3DBase< float, GlobalTag >
GlobalTrackerMuonAlignment::CLHEP_dot
double CLHEP_dot(const CLHEP::HepVector &a, const CLHEP::HepVector &b)
Definition: GlobalTrackerMuonAlignment.cc:148
b
double b
Definition: hdecay.h:118
GlobalTrackerMuonAlignment::iteratorHcalRcd
std::vector< AlignTransform >::const_iterator iteratorHcalRcd
Definition: GlobalTrackerMuonAlignment.cc:212
edm::OwnVector::const_iterator
Definition: OwnVector.h:41
AlignTransform
Definition: AlignTransform.h:15
TkRotation::yx
T yx() const
Definition: extTkRotation.h:254
TrajectoryStateOnSurface::localParameters
const LocalTrajectoryParameters & localParameters() const
Definition: TrajectoryStateOnSurface.h:73
GlobalTrackerMuonAlignment::debug_
bool debug_
Definition: GlobalTrackerMuonAlignment.cc:173
PbPb_ZMuSkimMuonDPG_cff.tracker
tracker
Definition: PbPb_ZMuSkimMuonDPG_cff.py:60
GlobalTrackerMuonAlignment::Hess
CLHEP::HepMatrix Hess
Definition: GlobalTrackerMuonAlignment.cc:200
LocalTrajectoryError
Definition: LocalTrajectoryError.h:20
LocalTrajectoryParameters::pzSign
float pzSign() const
Sign of the z-component of the momentum in the local frame.
Definition: LocalTrajectoryParameters.h:146
GlobalTrackerMuonAlignment::MuRHBuilder
MuonTransientTrackingRecHitBuilder * MuRHBuilder
Definition: GlobalTrackerMuonAlignment.cc:192
Surface::tangentPlane
virtual ConstReferenceCountingPointer< TangentPlane > tangentPlane(const GlobalPoint &) const =0
GlobalTrackerMuonAlignment::histo12
TH1F * histo12
Definition: GlobalTrackerMuonAlignment.cc:232
GlobalTrackerMuonAlignment::MuSelect
std::string MuSelect
Definition: GlobalTrackerMuonAlignment.cc:217
a
double a
Definition: hdecay.h:119
DetId::Tracker
Definition: DetId.h:25
TkRotation::xx
T xx() const
Definition: extTkRotation.h:251
GlobalTrackerMuonAlignment::iteratorMuonRcd
std::vector< AlignTransform >::const_iterator iteratorMuonRcd
Definition: GlobalTrackerMuonAlignment.cc:210
SmartPropagator::propagate
TrajectoryStateOnSurface propagate(STA const &state, SUR const &surface) const
Definition: Propagator.h:50
GlobalTrackerMuonAlignment::histo4
TH1F * histo4
Definition: GlobalTrackerMuonAlignment.cc:225
Trajectory::lastMeasurement
TrajectoryMeasurement const & lastMeasurement() const
Definition: Trajectory.h:150
GlobalTrackerMuonAlignment::muonTags_
edm::InputTag muonTags_
Definition: GlobalTrackerMuonAlignment.cc:160
GlobalTrackerMuonAlignment::gmuonTags_
edm::InputTag gmuonTags_
Definition: GlobalTrackerMuonAlignment.cc:161
dumpMFGeometry_cfg.delta
delta
Definition: dumpMFGeometry_cfg.py:25
GlobalTrackerMuonAlignment::histo5
TH1F * histo5
Definition: GlobalTrackerMuonAlignment.cc:226
GlobalTrackerMuonAlignment::MuGlShift
CLHEP::HepVector MuGlShift
Definition: GlobalTrackerMuonAlignment.cc:214
SmartPropagator_cfi.SmartPropagator
SmartPropagator
Definition: SmartPropagator_cfi.py:3
PV3DBase::y
T y() const
Definition: PV3DBase.h:60
edm::Service< cond::service::PoolDBOutputService >
createfilelist.int
int
Definition: createfilelist.py:10
iEvent
int iEvent
Definition: GenABIO.cc:224
alignmentValidation.c1
c1
do drawing
Definition: alignmentValidation.py:1025
GlobalTrackerMuonAlignment::propagator_
string propagator_
Definition: GlobalTrackerMuonAlignment.cc:163
GlobalTrackerMuonAlignment::refitMuon_
bool refitMuon_
Definition: GlobalTrackerMuonAlignment.cc:168
MuonTransientTrackingRecHitBuilder::build
RecHitPointer build(const TrackingRecHit *p, edm::ESHandle< GlobalTrackingGeometry > trackingGeometry) const
Call the MuonTransientTrackingRecHit::specificBuild.
Definition: MuonTransientTrackingRecHitBuilder.cc:20
TkRotation::xy
T xy() const
Definition: extTkRotation.h:252
TkRotation::zz
T zz() const
Definition: extTkRotation.h:259
reco::TransientTrack::outermostMeasurementState
TrajectoryStateOnSurface outermostMeasurementState() const
Definition: TransientTrack.h:86
PropagationDirectionChooser
Definition: PropagationDirectionChooser.h:15
DetId::Ecal
Definition: DetId.h:27
AlgebraicVector5
ROOT::Math::SVector< double, 5 > AlgebraicVector5
Definition: AlgebraicROOTObjects.h:14
get
#define get
AlgebraicVector6
ROOT::Math::SVector< double, 6 > AlgebraicVector6
Definition: AlgebraicROOTObjects.h:15
GlobalTrackerMuonAlignment::writeGlPosRcd
void writeGlPosRcd(CLHEP::HepVector &d3)
Definition: GlobalTrackerMuonAlignment.cc:3322
Vector3DBase::dot
PreciseFloatType< T, U >::Type dot(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:99
edm::ESHandleBase::isValid
bool isValid() const
Definition: ESHandle.h:44
CartesianTrajectoryError::matrix
const AlgebraicSymMatrix66 & matrix() const
Definition: CartesianTrajectoryError.h:28
GlobalTrackerMuonAlignment::GradL
CLHEP::HepVector GradL
Definition: GlobalTrackerMuonAlignment.cc:202
GlobalTrackerMuonAlignment::iteratorEcalRcd
std::vector< AlignTransform >::const_iterator iteratorEcalRcd
Definition: GlobalTrackerMuonAlignment.cc:211
SmartPropagator
Definition: SmartPropagator.h:34
alignCSCRings.r
r
Definition: alignCSCRings.py:93
GlobalTrackerMuonAlignment::debugTrajectory
void debugTrajectory(const std::string, Trajectory &)
Definition: GlobalTrackerMuonAlignment.cc:3302
TkRotation::yz
T yz() const
Definition: extTkRotation.h:256
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:64
GloballyPositioned::position
const PositionType & position() const
Definition: GloballyPositioned.h:36
TrajectoryStateOnSurface::globalMomentum
GlobalVector globalMomentum() const
Definition: TrajectoryStateOnSurface.h:66
mag
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
Definition: Basic3DVectorLD.h:127
TtFullHadDaughter::B
static const std::string B
Definition: TtFullHadronicEvent.h:9
AlgebraicVector4
ROOT::Math::SVector< double, 4 > AlgebraicVector4
Definition: AlgebraicROOTObjects.h:13
MaterialEffects_cfi.A
A
Definition: MaterialEffects_cfi.py:11
GlobalTrackerMuonAlignment::histo31
TH1F * histo31
Definition: GlobalTrackerMuonAlignment.cc:251
GlobalTrackerMuonAlignment::theFitter
KFTrajectoryFitter * theFitter
Definition: GlobalTrackerMuonAlignment.cc:187
DetId::rawId
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
Trajectory::firstMeasurement
TrajectoryMeasurement const & firstMeasurement() const
Definition: Trajectory.h:166
reco::TransientTrack
Definition: TransientTrack.h:19
GlobalTrackerMuonAlignment::histo17
TH1F * histo17
Definition: GlobalTrackerMuonAlignment.cc:237
GlobalTrackerMuonAlignment::debugTrajectorySOSv
void debugTrajectorySOSv(const std::string, TrajectoryStateOnSurface)
Definition: GlobalTrackerMuonAlignment.cc:3276
T
long double T
Definition: Basic3DVectorLD.h:48
KFTrajectorySmoother
Definition: KFTrajectorySmoother.h:20
TkRotation::yy
T yy() const
Definition: extTkRotation.h:255
Exception
Definition: hltDiff.cc:246
cond::service::PoolDBOutputService::writeOne
Hash writeOne(const T *payload, Time_t time, const std::string &recordName)
Definition: PoolDBOutputService.h:63
GlobalTrackerMuonAlignment::histo3
TH1F * histo3
Definition: GlobalTrackerMuonAlignment.cc:224
GlobalTrackerMuonAlignment::trackFitter
void trackFitter(reco::TrackRef, reco::TransientTrack &, PropagationDirection, TrajectoryStateOnSurface &)
Definition: GlobalTrackerMuonAlignment.cc:2997
Plane
Definition: Plane.h:16
GlobalTrackerMuonAlignment::histo22
TH1F * histo22
Definition: GlobalTrackerMuonAlignment.cc:242
GlobalTrackerMuonAlignment::histo102
TH2F * histo102
Definition: GlobalTrackerMuonAlignment.cc:258
GlobalTrackerMuonAlignment::debugTrackHit
void debugTrackHit(const std::string, reco::TrackRef)
Definition: GlobalTrackerMuonAlignment.cc:3231
TrajectorySeed
Definition: TrajectorySeed.h:18
GlobalTrackerMuonAlignment::HessL
CLHEP::HepMatrix HessL
Definition: GlobalTrackerMuonAlignment.cc:203
GlobalTrackerMuonAlignment::histo23
TH1F * histo23
Definition: GlobalTrackerMuonAlignment.cc:243
GlobalTrackerMuonAlignment::gradientLocal
void gradientLocal(GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, CLHEP::HepSymMatrix &, CLHEP::HepMatrix &, CLHEP::HepVector &, AlgebraicVector4 &)
Definition: GlobalTrackerMuonAlignment.cc:2308
GlobalTrackerMuonAlignment::gradientGlobal
void gradientGlobal(GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, AlgebraicSymMatrix66 &)
Definition: GlobalTrackerMuonAlignment.cc:2123
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
GlobalTrackerMuonAlignment::histo20
TH1F * histo20
Definition: GlobalTrackerMuonAlignment.cc:240
TrajectoryStateOnSurface::surface
const SurfaceType & surface() const
Definition: TrajectoryStateOnSurface.h:78
bsc_activity_cfg.ecal
ecal
Definition: bsc_activity_cfg.py:25
GlobalTrackerMuonAlignment::debugTrajectorySOS
void debugTrajectorySOS(const std::string, TrajectoryStateOnSurface &)
Definition: GlobalTrackerMuonAlignment.cc:3250
MillePedeFileConverter_cfg.out
out
Definition: MillePedeFileConverter_cfg.py:31
Alignments::m_align
std::vector< AlignTransform > m_align
Definition: Alignments.h:19
GlobalTrackerMuonAlignment::watchMagneticFieldRecord_
edm::ESWatcher< IdealMagneticFieldRecord > watchMagneticFieldRecord_
Definition: GlobalTrackerMuonAlignment.cc:179
ztail.d
d
Definition: ztail.py:151
GlobalTrackerMuonAlignment::file
TFile * file
Definition: GlobalTrackerMuonAlignment.cc:221
GlobalTrackerMuonAlignment::defineFitter
bool defineFitter
Definition: GlobalTrackerMuonAlignment.cc:191
cms::Exception
Definition: Exception.h:70
DetId::Muon
Definition: DetId.h:26
TrajectoryStateOnSurface::globalParameters
const GlobalTrajectoryParameters & globalParameters() const
Definition: TrajectoryStateOnSurface.h:64
PTrajectoryStateOnDet
Definition: PTrajectoryStateOnDet.h:10
HGC3DClusterGenMatchSelector_cfi.dR
dR
Definition: HGC3DClusterGenMatchSelector_cfi.py:7
calo
Definition: Common.h:9
GlobalTrackerMuonAlignment::watchGlobalPositionRcd_
edm::ESWatcher< GlobalPositionRcd > watchGlobalPositionRcd_
Definition: GlobalTrackerMuonAlignment.cc:184
GloballyPositioned::rotation
const RotationType & rotation() const
Definition: GloballyPositioned.h:38
GlobalTrackerMuonAlignment::globalPositionRcd_
const Alignments * globalPositionRcd_
Definition: GlobalTrackerMuonAlignment.cc:185
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
GlobalTrackerMuonAlignment::histo29
TH1F * histo29
Definition: GlobalTrackerMuonAlignment.cc:249
HLT_FULL_cff.R0
R0
Definition: HLT_FULL_cff.py:8629
GlobalTrackerMuonAlignment::histo25
TH1F * histo25
Definition: GlobalTrackerMuonAlignment.cc:245
TrajectoryStateOnSurface::localError
const LocalTrajectoryError & localError() const
Definition: TrajectoryStateOnSurface.h:77
MuonTransientTrackingRecHitBuilder
Definition: MuonTransientTrackingRecHitBuilder.h:12
AlgebraicSymMatrix55
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
Definition: AlgebraicROOTObjects.h:23
TkRotation::zx
T zx() const
Definition: extTkRotation.h:257
GlobalTrackerMuonAlignment::histo
TH1F * histo
Definition: GlobalTrackerMuonAlignment.cc:222
edm::OwnVector::size
size_type size() const
Definition: OwnVector.h:300
A
Alignments
Definition: Alignments.h:10
edm::OwnVector::clear
void clear()
Definition: OwnVector.h:481
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:69
cuy.ii
ii
Definition: cuy.py:590
edm::InputTag
Definition: InputTag.h:15
alongMomentum
Definition: PropagationDirection.h:4
KFUpdator
Definition: KFUpdator.h:32
cond::service::PoolDBOutputService::currentTime
cond::Time_t currentTime() const
Definition: PoolDBOutputService.cc:217
GlobalTrackerMuonAlignment::fitHist
void fitHist()
Definition: GlobalTrackerMuonAlignment.cc:549
TrajectoryStateOnSurface::isValid
bool isValid() const
Definition: TrajectoryStateOnSurface.h:54
Trajectory::isValid
bool isValid() const
Definition: Trajectory.h:257
hit
Definition: SiStripHitEffFromCalibTree.cc:88
edm::OwnVector< TrackingRecHit >
edm::eventsetup::Itr
RecordProviders::iterator Itr
Definition: EventSetupProvider.cc:351
GlobalTrackerMuonAlignment::theFitterOp
KFTrajectoryFitter * theFitterOp
Definition: GlobalTrackerMuonAlignment.cc:189
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
reco::TransientTrack::innermostMeasurementState
TrajectoryStateOnSurface innermostMeasurementState() const
Definition: TransientTrack.h:88
GlobalTrackerMuonAlignment::histo27
TH1F * histo27
Definition: GlobalTrackerMuonAlignment.cc:247
TrackingComponentsRecord
Definition: TrackingComponentsRecord.h:12
Plane::normalVector
GlobalVector normalVector() const
Definition: Plane.h:41
Generated for CMSSW Reference Manual by doxygen 1.8.16