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 const &)=delete
 
 EDConsumerBase (EDConsumerBase &&)=default
 
ESProxyIndex const * esGetTokenIndices (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 const & operator= (EDConsumerBase const &)=delete
 
EDConsumerBase & operator= (EDConsumerBase &&)=default
 
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 wantsStreamLuminosityBlocks ()
 
static bool wantsStreamRuns ()
 
- Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType > consumes (edm::InputTag const &tag)
 
EDGetToken consumes (const TypeToGet &id, edm::InputTag const &tag)
 
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<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto esConsumes (ESInputTag const &tag)
 
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType > mayConsume (edm::InputTag const &tag)
 
EDGetToken mayConsume (const TypeToGet &id, edm::InputTag const &tag)
 
template<BranchType B>
EDGetToken mayConsume (const TypeToGet &id, 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 ( const edm::ParameterSet &  iConfig)
explicit

Definition at line 272 of file GlobalTrackerMuonAlignment.cc.

References cosmicMuonMode_, Exception, and isolatedMuonMode_.

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 }
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
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

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

Definition at line 312 of file GlobalTrackerMuonAlignment.cc.

References analyzeTrackTrajectory().

312  {
313  //GlobalTrackerMuonAlignment::analyzeTrackTrack(iEvent, iSetup);
314  GlobalTrackerMuonAlignment::analyzeTrackTrajectory(iEvent, iSetup);
315 }
GlobalTrackerMuonAlignment::analyzeTrackTrajectory
void analyzeTrackTrajectory(const edm::Event &, const edm::EventSetup &)
Definition: GlobalTrackerMuonAlignment.cc:1252
void GlobalTrackerMuonAlignment::analyzeTrackTrack ( const edm::Event &  iEvent,
const edm::EventSetup &  iSetup 
)

Definition at line 573 of file GlobalTrackerMuonAlignment.cc.

References alongMomentum, edm::EventBase::bunchCrossing(), 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(), edm::Event::getByLabel(), 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, info(), createfilelist::int, isolatedMuonMode_, TrajectoryStateOnSurface::isValid(), iteratorEcalRcd, iteratorHcalRcd, iteratorMuonRcd, iteratorTrackerRcd, TrajectoryStateOnSurface::localError(), TrajectoryStateOnSurface::localParameters(), funct::m, Alignments::m_align, PV3DBase< T, PVType, FrameType >::mag(), mag(), seedCreatorFromRegionConsecutiveHitsEDProducer_cff::magneticField, magneticField_, CartesianTrajectoryError::matrix(), LocalTrajectoryError::matrix(), misalignMuonL(), DetId::Muon, extraflags_cff::muons, muonTags_, MuSelect, N_event, N_track, Plane::normalVector(), oppositeToMomentum, PV3DBase< T, PVType, FrameType >::perp(), PI, GloballyPositioned< T >::position(), Propagator::propagate(), PhotonConversionTrajectorySeedProducerFromQuadruplets_cfi::propagator, defaultRKPropagator::Product::propagator, propagator_, reco::tau::disc::Pt(), DetId::rawId(), GloballyPositioned< T >::rotation(), edm::Event::run(), SmartPropagator_cfi::SmartPropagator, smuonTags_, mathSSE::sqrt(), SteppingHelixPropagator_cfi::SteppingHelixPropagator, TrajectoryStateOnSurface::surface(), Surface::tangentPlane(), DetId::Tracker, trackingGeometry_, l1t::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().

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
TkRotation::xx
T xx() const
Definition: extTkRotation.h:251
GlobalTrackerMuonAlignment::trackingGeometry_
const GlobalTrackingGeometry * trackingGeometry_
Definition: GlobalTrackerMuonAlignment.cc:177
GlobalTrackerMuonAlignment::watchGlobalPositionRcd_
edm::ESWatcher< GlobalPositionRcd > watchGlobalPositionRcd_
Definition: GlobalTrackerMuonAlignment.cc:184
info
static const TGPicture * info(bool iBackgroundIsBlack)
Definition: FWCollectionSummaryWidget.cc:170
Surface::tangentPlane
virtual ConstReferenceCountingPointer< TangentPlane > tangentPlane(const GlobalPoint &) const =0
mps_fire.i
i
Definition: mps_fire.py:338
GlobalTrackerMuonAlignment::gradientGlobalAlg
void gradientGlobalAlg(GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, AlgebraicSymMatrix66 &)
Definition: GlobalTrackerMuonAlignment.cc:2032
GlobalTrackerMuonAlignment::gradientLocal
void gradientLocal(GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, CLHEP::HepSymMatrix &, CLHEP::HepMatrix &, CLHEP::HepVector &, AlgebraicVector4 &)
Definition: GlobalTrackerMuonAlignment.cc:2308
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:72
TrajectoryStateOnSurface::localParameters
const LocalTrajectoryParameters & localParameters() const
Definition: TrajectoryStateOnSurface.h:107
mag
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
Definition: Basic3DVectorLD.h:156
Plane::normalVector
GlobalVector normalVector() const
Definition: Plane.h:41
SteppingHelixPropagator_cfi.SteppingHelixPropagator
SteppingHelixPropagator
Definition: SteppingHelixPropagator_cfi.py:3
AlgebraicSymMatrix66
ROOT::Math::SMatrix< double, 6, 6, ROOT::Math::MatRepSym< double, 6 > > AlgebraicSymMatrix66
Definition: AlgebraicROOTObjects.h:24
DetId::rawId
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:50
TrajectoryStateOnSurface::cartesianError
const CartesianTrajectoryError cartesianError() const
Definition: TrajectoryStateOnSurface.h:101
PV3DBase::y
T y() const
Definition: PV3DBase.h:63
TkRotation::yx
T yx() const
Definition: extTkRotation.h:254
edm::EventBase::bunchCrossing
int bunchCrossing() const
Definition: EventBase.h:64
GlobalTrackerMuonAlignment::iteratorHcalRcd
std::vector< AlignTransform >::const_iterator iteratorHcalRcd
Definition: GlobalTrackerMuonAlignment.cc:212
TrajectoryStateOnSurface::globalPosition
GlobalPoint globalPosition() const
Definition: TrajectoryStateOnSurface.h:83
Vector3DBase::dot
PreciseFloatType< T, U >::Type dot(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:107
AlgebraicSymMatrix55
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
Definition: AlgebraicROOTObjects.h:23
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
GlobalTrackerMuonAlignment::gradientGlobal
void gradientGlobal(GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, AlgebraicSymMatrix66 &)
Definition: GlobalTrackerMuonAlignment.cc:2123
Alignments::m_align
std::vector< AlignTransform > m_align
Definition: Alignments.h:19
AlgebraicVector6
ROOT::Math::SVector< double, 6 > AlgebraicVector6
Definition: AlgebraicROOTObjects.h:15
Plane
Definition: Plane.h:17
TkRotation::zx
T zx() const
Definition: extTkRotation.h:257
TkRotation::xy
T xy() const
Definition: extTkRotation.h:252
LocalTrajectoryParameters::vector
AlgebraicVector5 vector() const
Definition: LocalTrajectoryParameters.h:125
TkRotation::zz
T zz() const
Definition: extTkRotation.h:259
TrajectoryStateOnSurface::surface
const SurfaceType & surface() const
Definition: TrajectoryStateOnSurface.h:122
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:67
GlobalTrackerMuonAlignment::watchTrackingGeometry_
edm::ESWatcher< GlobalTrackingGeometryRecord > watchTrackingGeometry_
Definition: GlobalTrackerMuonAlignment.cc:175
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
TrajectoryStateOnSurface::freeState
FreeTrajectoryState const * freeState(bool withErrors=true) const
Definition: TrajectoryStateOnSurface.h:70
PV3DBase::z
T z() const
Definition: PV3DBase.h:64
edm::Event::run
RunNumber_t run() const
Definition: Event.h:101
seedCreatorFromRegionConsecutiveHitsEDProducer_cff.magneticField
magneticField
Definition: seedCreatorFromRegionConsecutiveHitsEDProducer_cff.py:4
GlobalTrackerMuonAlignment::iteratorMuonRcd
std::vector< AlignTransform >::const_iterator iteratorMuonRcd
Definition: GlobalTrackerMuonAlignment.cc:210
TkRotation::zy
T zy() const
Definition: extTkRotation.h:258
extraflags_cff.muons
muons
Definition: extraflags_cff.py:17
LocalTrajectoryError::matrix
const AlgebraicSymMatrix55 & matrix() const
Definition: LocalTrajectoryError.h:61
TrajectoryStateOnSurface::localError
const LocalTrajectoryError & localError() const
Definition: TrajectoryStateOnSurface.h:119
TkRotation::yy
T yy() const
Definition: extTkRotation.h:255
PI
#define PI
Definition: QcdUeDQM.h:36
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:480
CartesianTrajectoryError::matrix
const AlgebraicSymMatrix66 & matrix() const
Definition: CartesianTrajectoryError.h:29
GlobalTrackerMuonAlignment::misalignMuonL
void misalignMuonL(GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, AlgebraicVector4 &, TrajectoryStateOnSurface &, TrajectoryStateOnSurface &)
Definition: GlobalTrackerMuonAlignment.cc:2715
GlobalTrackerMuonAlignment::iteratorEcalRcd
std::vector< AlignTransform >::const_iterator iteratorEcalRcd
Definition: GlobalTrackerMuonAlignment.cc:211
DetId
Definition: DetId.h:18
AlgebraicVector5
ROOT::Math::SVector< double, 5 > AlgebraicVector5
Definition: AlgebraicROOTObjects.h:14
SmartPropagator_cfi.SmartPropagator
SmartPropagator
Definition: SmartPropagator_cfi.py:3
GlobalTrackerMuonAlignment::watchMagneticFieldRecord_
edm::ESWatcher< IdealMagneticFieldRecord > watchMagneticFieldRecord_
Definition: GlobalTrackerMuonAlignment.cc:179
Propagator::propagate
TrajectoryStateOnSurface propagate(STA const &state, SUR const &surface) const
Definition: Propagator.h:53
edm::ESWatcher::check
bool check(const edm::EventSetup &iSetup)
Definition: ESWatcher.h:52
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:44
edm
HLT enums.
Definition: AlignableModifier.h:17
TrajectoryStateOnSurface::globalMomentum
GlobalVector globalMomentum() const
Definition: TrajectoryStateOnSurface.h:86
TkRotation::xz
T xz() const
Definition: extTkRotation.h:253
edm::EventSetup::get
T get() const
Definition: EventSetup.h:71
GloballyPositioned::rotation
const RotationType & rotation() const
Definition: GloballyPositioned.h:44
GlobalTrackerMuonAlignment::iteratorTrackerRcd
std::vector< AlignTransform >::const_iterator iteratorTrackerRcd
Definition: GlobalTrackerMuonAlignment.cc:209
gather_cfg.cout
cout
Definition: gather_cfg.py:144
PhotonConversionTrajectorySeedProducerFromQuadruplets_cfi.propagator
propagator
Definition: PhotonConversionTrajectorySeedProducerFromQuadruplets_cfi.py:41
AlgebraicVector4
ROOT::Math::SVector< double, 4 > AlgebraicVector4
Definition: AlgebraicROOTObjects.h:13
PV3DBase::x
T x() const
Definition: PV3DBase.h:62
GloballyPositioned::position
const PositionType & position() const
Definition: GloballyPositioned.h:42
TkRotation::yz
T yz() const
Definition: extTkRotation.h:256
GlobalVector
Global3DVector GlobalVector
Definition: GlobalVector.h:10
void GlobalTrackerMuonAlignment::analyzeTrackTrajectory ( const edm::Event &  iEvent,
const edm::EventSetup &  iSetup 
)

Definition at line 1252 of file GlobalTrackerMuonAlignment.cc.

References alongMomentum, edm::EventBase::bunchCrossing(), 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(), edm::Event::getByLabel(), 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, info(), createfilelist::int, isolatedMuonMode_, edm::ESHandleBase::isValid(), TrajectoryStateOnSurface::isValid(), iteratorEcalRcd, iteratorHcalRcd, iteratorMuonRcd, iteratorTrackerRcd, TrajectoryStateOnSurface::localError(), TrajectoryStateOnSurface::localParameters(), funct::m, Alignments::m_align, PV3DBase< T, PVType, FrameType >::mag(), mag(), seedCreatorFromRegionConsecutiveHitsEDProducer_cff::magneticField, magneticField_, CartesianTrajectoryError::matrix(), LocalTrajectoryError::matrix(), misalignMuonL(), DetId::Muon, muonFitter(), extraflags_cff::muons, muonTags_, MuRHBuilder, MuSelect, N_event, N_track, Plane::normalVector(), oppositeToMomentum, PV3DBase< T, PVType, FrameType >::perp(), PI, GloballyPositioned< T >::position(), Propagator::propagate(), PhotonConversionTrajectorySeedProducerFromQuadruplets_cfi::propagator, defaultRKPropagator::Product::propagator, propagator_, reco::tau::disc::Pt(), DetId::rawId(), refitMuon_, refitTrack_, GloballyPositioned< T >::rotation(), edm::Event::run(), SmartPropagator_cfi::SmartPropagator, smuonTags_, mathSSE::sqrt(), SteppingHelixPropagator_cfi::SteppingHelixPropagator, TrajectoryStateOnSurface::surface(), Surface::tangentPlane(), theFitter, theFitterOp, theSmoother, theSmootherOp, theTrackingGeometry, DetId::Tracker, trackFitter(), trackingGeometry_, l1t::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().

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
TkRotation::xx
T xx() const
Definition: extTkRotation.h:251
GlobalTrackerMuonAlignment::trackingGeometry_
const GlobalTrackingGeometry * trackingGeometry_
Definition: GlobalTrackerMuonAlignment.cc:177
GlobalTrackerMuonAlignment::watchGlobalPositionRcd_
edm::ESWatcher< GlobalPositionRcd > watchGlobalPositionRcd_
Definition: GlobalTrackerMuonAlignment.cc:184
GlobalTrackerMuonAlignment::theTrackingGeometry
edm::ESHandle< GlobalTrackingGeometry > theTrackingGeometry
Definition: GlobalTrackerMuonAlignment.cc:176
info
static const TGPicture * info(bool iBackgroundIsBlack)
Definition: FWCollectionSummaryWidget.cc:170
Surface::tangentPlane
virtual ConstReferenceCountingPointer< TangentPlane > tangentPlane(const GlobalPoint &) const =0
mps_fire.i
i
Definition: mps_fire.py:338
GlobalTrackerMuonAlignment::gradientGlobalAlg
void gradientGlobalAlg(GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, AlgebraicSymMatrix66 &)
Definition: GlobalTrackerMuonAlignment.cc:2032
GlobalTrackerMuonAlignment::gradientLocal
void gradientLocal(GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, CLHEP::HepSymMatrix &, CLHEP::HepMatrix &, CLHEP::HepVector &, AlgebraicVector4 &)
Definition: GlobalTrackerMuonAlignment.cc:2308
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:72
TrajectoryStateOnSurface::localParameters
const LocalTrajectoryParameters & localParameters() const
Definition: TrajectoryStateOnSurface.h:107
mag
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
Definition: Basic3DVectorLD.h:156
Plane::normalVector
GlobalVector normalVector() const
Definition: Plane.h:41
SteppingHelixPropagator_cfi.SteppingHelixPropagator
SteppingHelixPropagator
Definition: SteppingHelixPropagator_cfi.py:3
AlgebraicSymMatrix66
ROOT::Math::SMatrix< double, 6, 6, ROOT::Math::MatRepSym< double, 6 > > AlgebraicSymMatrix66
Definition: AlgebraicROOTObjects.h:24
DetId::rawId
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:50
TrajectoryStateOnSurface::cartesianError
const CartesianTrajectoryError cartesianError() const
Definition: TrajectoryStateOnSurface.h:101
PV3DBase::y
T y() const
Definition: PV3DBase.h:63
TkRotation::yx
T yx() const
Definition: extTkRotation.h:254
edm::EventBase::bunchCrossing
int bunchCrossing() const
Definition: EventBase.h:64
GlobalTrackerMuonAlignment::iteratorHcalRcd
std::vector< AlignTransform >::const_iterator iteratorHcalRcd
Definition: GlobalTrackerMuonAlignment.cc:212
TrajectoryStateOnSurface::globalPosition
GlobalPoint globalPosition() const
Definition: TrajectoryStateOnSurface.h:83
Vector3DBase::dot
PreciseFloatType< T, U >::Type dot(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:107
AlgebraicSymMatrix55
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
Definition: AlgebraicROOTObjects.h:23
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
GlobalTrackerMuonAlignment::gradientGlobal
void gradientGlobal(GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, AlgebraicSymMatrix66 &)
Definition: GlobalTrackerMuonAlignment.cc:2123
Alignments::m_align
std::vector< AlignTransform > m_align
Definition: Alignments.h:19
AlgebraicVector6
ROOT::Math::SVector< double, 6 > AlgebraicVector6
Definition: AlgebraicROOTObjects.h:15
Plane
Definition: Plane.h:17
TkRotation::zx
T zx() const
Definition: extTkRotation.h:257
TkRotation::xy
T xy() const
Definition: extTkRotation.h:252
LocalTrajectoryParameters::vector
AlgebraicVector5 vector() const
Definition: LocalTrajectoryParameters.h:125
Chi2MeasurementEstimator_cfi.Chi2MeasurementEstimator
Chi2MeasurementEstimator
Definition: Chi2MeasurementEstimator_cfi.py:5
TkRotation::zz
T zz() const
Definition: extTkRotation.h:259
TrajectoryStateOnSurface::surface
const SurfaceType & surface() const
Definition: TrajectoryStateOnSurface.h:122
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:67
GlobalTrackerMuonAlignment::watchTrackingGeometry_
edm::ESWatcher< GlobalTrackingGeometryRecord > watchTrackingGeometry_
Definition: GlobalTrackerMuonAlignment.cc:175
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
TrajectoryStateOnSurface::freeState
FreeTrajectoryState const * freeState(bool withErrors=true) const
Definition: TrajectoryStateOnSurface.h:70
PV3DBase::z
T z() const
Definition: PV3DBase.h:64
edm::Event::run
RunNumber_t run() const
Definition: Event.h:101
seedCreatorFromRegionConsecutiveHitsEDProducer_cff.magneticField
magneticField
Definition: seedCreatorFromRegionConsecutiveHitsEDProducer_cff.py:4
GlobalTrackerMuonAlignment::iteratorMuonRcd
std::vector< AlignTransform >::const_iterator iteratorMuonRcd
Definition: GlobalTrackerMuonAlignment.cc:210
TkRotation::zy
T zy() const
Definition: extTkRotation.h:258
extraflags_cff.muons
muons
Definition: extraflags_cff.py:17
LocalTrajectoryError::matrix
const AlgebraicSymMatrix55 & matrix() const
Definition: LocalTrajectoryError.h:61
TrajectoryStateOnSurface::localError
const LocalTrajectoryError & localError() const
Definition: TrajectoryStateOnSurface.h:119
TkRotation::yy
T yy() const
Definition: extTkRotation.h:255
PI
#define PI
Definition: QcdUeDQM.h:36
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:480
CartesianTrajectoryError::matrix
const AlgebraicSymMatrix66 & matrix() const
Definition: CartesianTrajectoryError.h:29
GlobalTrackerMuonAlignment::misalignMuonL
void misalignMuonL(GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, AlgebraicVector4 &, TrajectoryStateOnSurface &, TrajectoryStateOnSurface &)
Definition: GlobalTrackerMuonAlignment.cc:2715
GlobalTrackerMuonAlignment::iteratorEcalRcd
std::vector< AlignTransform >::const_iterator iteratorEcalRcd
Definition: GlobalTrackerMuonAlignment.cc:211
DetId
Definition: DetId.h:18
GlobalTrackerMuonAlignment::trackFitter
void trackFitter(reco::TrackRef, reco::TransientTrack &, PropagationDirection, TrajectoryStateOnSurface &)
Definition: GlobalTrackerMuonAlignment.cc:2997
AlgebraicVector5
ROOT::Math::SVector< double, 5 > AlgebraicVector5
Definition: AlgebraicROOTObjects.h:14
GlobalTrackerMuonAlignment::TTRHBuilder
const TransientTrackingRecHitBuilder * TTRHBuilder
Definition: GlobalTrackerMuonAlignment.cc:193
SmartPropagator_cfi.SmartPropagator
SmartPropagator
Definition: SmartPropagator_cfi.py:3
GlobalTrackerMuonAlignment::watchMagneticFieldRecord_
edm::ESWatcher< IdealMagneticFieldRecord > watchMagneticFieldRecord_
Definition: GlobalTrackerMuonAlignment.cc:179
Propagator::propagate
TrajectoryStateOnSurface propagate(STA const &state, SUR const &surface) const
Definition: Propagator.h:53
edm::ESWatcher::check
bool check(const edm::EventSetup &iSetup)
Definition: ESWatcher.h:52
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:44
edm
HLT enums.
Definition: AlignableModifier.h:17
TrajectoryStateOnSurface::globalMomentum
GlobalVector globalMomentum() const
Definition: TrajectoryStateOnSurface.h:86
TkRotation::xz
T xz() const
Definition: extTkRotation.h:253
edm::EventSetup::get
T get() const
Definition: EventSetup.h:71
GloballyPositioned::rotation
const RotationType & rotation() const
Definition: GloballyPositioned.h:44
GlobalTrackerMuonAlignment::iteratorTrackerRcd
std::vector< AlignTransform >::const_iterator iteratorTrackerRcd
Definition: GlobalTrackerMuonAlignment.cc:209
gather_cfg.cout
cout
Definition: gather_cfg.py:144
PhotonConversionTrajectorySeedProducerFromQuadruplets_cfi.propagator
propagator
Definition: PhotonConversionTrajectorySeedProducerFromQuadruplets_cfi.py:41
edm::ESHandleBase::isValid
bool isValid() const
Definition: ESHandle.h:44
AlgebraicVector4
ROOT::Math::SVector< double, 4 > AlgebraicVector4
Definition: AlgebraicROOTObjects.h:13
PV3DBase::x
T x() const
Definition: PV3DBase.h:62
GloballyPositioned::position
const PositionType & position() const
Definition: GloballyPositioned.h:42
TkRotation::yz
T yz() const
Definition: extTkRotation.h:256
GlobalTrackerMuonAlignment::muonFitter
void muonFitter(reco::TrackRef, reco::TransientTrack &, PropagationDirection, TrajectoryStateOnSurface &)
Definition: GlobalTrackerMuonAlignment.cc:3112
GlobalTrackerMuonAlignment::MuRHBuilder
MuonTransientTrackingRecHitBuilder * MuRHBuilder
Definition: GlobalTrackerMuonAlignment.cc:192
GlobalVector
Global3DVector GlobalVector
Definition: GlobalVector.h:10
GlobalTrackerMuonAlignment::debugTrajectorySOS
void debugTrajectorySOS(const std::string, TrajectoryStateOnSurface &)
Definition: GlobalTrackerMuonAlignment.cc:3250
void GlobalTrackerMuonAlignment::beginJob ( void  )
overrideprivatevirtual

Reimplemented from edm::EDAnalyzer.

Definition at line 318 of file GlobalTrackerMuonAlignment.cc.

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

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 }
mps_fire.i
i
Definition: mps_fire.py:338
void GlobalTrackerMuonAlignment::bookHist ( )

Definition at line 471 of file GlobalTrackerMuonAlignment.cc.

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, histo8, and PI.

Referenced by beginJob().

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 }
PI
#define PI
Definition: QcdUeDQM.h:36
double GlobalTrackerMuonAlignment::CLHEP_dot ( const CLHEP::HepVector &  a,
const CLHEP::HepVector &  b 
)
inline

Definition at line 148 of file GlobalTrackerMuonAlignment.cc.

References a, analyze(), b, and bk::beginJob().

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

148  {
149  return a(1) * b(1) + a(2) * b(2) + a(3) * b(3);
150  }
b
double b
Definition: hdecay.h:120
a
double a
Definition: hdecay.h:121
void GlobalTrackerMuonAlignment::debugTrackHit ( const std::string  title,
reco::TrackRef  alongTr 
)

Definition at line 3231 of file GlobalTrackerMuonAlignment.cc.

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

Referenced by muonFitter(), and trackFitter().

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 }
gather_cfg.cout
cout
Definition: gather_cfg.py:144
void GlobalTrackerMuonAlignment::debugTrackHit ( const std::string  title,
reco::TransientTrack &  alongTr 
)

Definition at line 3212 of file GlobalTrackerMuonAlignment.cc.

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

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 }
mps_fire.i
i
Definition: mps_fire.py:338
reco::TransientTrack::recHitsEnd
trackingRecHit_iterator recHitsEnd() const
last iterator to RecHits
Definition: TransientTrack.h:132
edm::OwnVector::const_iterator
Definition: OwnVector.h:41
gather_cfg.cout
cout
Definition: gather_cfg.py:144
reco::TransientTrack::recHitsBegin
trackingRecHit_iterator recHitsBegin() const
first iterator to RecHits
Definition: TransientTrack.h:130
void GlobalTrackerMuonAlignment::debugTrajectory ( const std::string  title,
Trajectory &  traj 
)

Definition at line 3302 of file GlobalTrackerMuonAlignment.cc.

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

Referenced by muonFitter(), and trackFitter().

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 }
Trajectory::foundHits
int foundHits() const
Definition: Trajectory.h:225
GlobalTrackerMuonAlignment::debugTrajectorySOSv
void debugTrajectorySOSv(const std::string, TrajectoryStateOnSurface)
Definition: GlobalTrackerMuonAlignment.cc:3276
Trajectory::direction
PropagationDirection const & direction() const
Definition: Trajectory.cc:140
Trajectory::lastMeasurement
TrajectoryMeasurement const & lastMeasurement() const
Definition: Trajectory.h:174
Trajectory::isValid
bool isValid() const
Definition: Trajectory.h:279
Trajectory::firstMeasurement
TrajectoryMeasurement const & firstMeasurement() const
Definition: Trajectory.h:187
Trajectory::chiSquared
float chiSquared() const
Definition: Trajectory.h:262
TrajectoryMeasurement::updatedState
TrajectoryStateOnSurface const & updatedState() const
Definition: TrajectoryMeasurement.h:180
gather_cfg.cout
cout
Definition: gather_cfg.py:144
void GlobalTrackerMuonAlignment::debugTrajectorySOS ( const std::string  title,
TrajectoryStateOnSurface &  trajSOS 
)

Definition at line 3250 of file GlobalTrackerMuonAlignment.cc.

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(), GlobalTrajectoryParameters::vector(), and LocalTrajectoryParameters::vector().

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

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 }
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:72
TrajectoryStateOnSurface::localParameters
const LocalTrajectoryParameters & localParameters() const
Definition: TrajectoryStateOnSurface.h:107
TrajectoryStateOnSurface::cartesianError
const CartesianTrajectoryError cartesianError() const
Definition: TrajectoryStateOnSurface.h:101
TrajectoryStateOnSurface::globalPosition
GlobalPoint globalPosition() const
Definition: TrajectoryStateOnSurface.h:83
LocalTrajectoryParameters::vector
AlgebraicVector5 vector() const
Definition: LocalTrajectoryParameters.h:125
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:67
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
LocalTrajectoryError::matrix
const AlgebraicSymMatrix55 & matrix() const
Definition: LocalTrajectoryError.h:61
TrajectoryStateOnSurface::localError
const LocalTrajectoryError & localError() const
Definition: TrajectoryStateOnSurface.h:119
CartesianTrajectoryError::matrix
const AlgebraicSymMatrix66 & matrix() const
Definition: CartesianTrajectoryError.h:29
TrajectoryStateOnSurface::globalParameters
const GlobalTrajectoryParameters & globalParameters() const
Definition: TrajectoryStateOnSurface.h:80
TrajectoryStateOnSurface::globalMomentum
GlobalVector globalMomentum() const
Definition: TrajectoryStateOnSurface.h:86
gather_cfg.cout
cout
Definition: gather_cfg.py:144
GlobalTrajectoryParameters::vector
AlgebraicVector6 vector() const
Definition: GlobalTrajectoryParameters.h:120
void GlobalTrackerMuonAlignment::debugTrajectorySOSv ( const std::string  title,
TrajectoryStateOnSurface  trajSOS 
)

Definition at line 3276 of file GlobalTrackerMuonAlignment.cc.

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(), GlobalTrajectoryParameters::vector(), and LocalTrajectoryParameters::vector().

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

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 }
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:72
TrajectoryStateOnSurface::localParameters
const LocalTrajectoryParameters & localParameters() const
Definition: TrajectoryStateOnSurface.h:107
TrajectoryStateOnSurface::cartesianError
const CartesianTrajectoryError cartesianError() const
Definition: TrajectoryStateOnSurface.h:101
TrajectoryStateOnSurface::globalPosition
GlobalPoint globalPosition() const
Definition: TrajectoryStateOnSurface.h:83
LocalTrajectoryParameters::vector
AlgebraicVector5 vector() const
Definition: LocalTrajectoryParameters.h:125
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:67
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
LocalTrajectoryError::matrix
const AlgebraicSymMatrix55 & matrix() const
Definition: LocalTrajectoryError.h:61
TrajectoryStateOnSurface::localError
const LocalTrajectoryError & localError() const
Definition: TrajectoryStateOnSurface.h:119
CartesianTrajectoryError::matrix
const AlgebraicSymMatrix66 & matrix() const
Definition: CartesianTrajectoryError.h:29
TrajectoryStateOnSurface::globalParameters
const GlobalTrajectoryParameters & globalParameters() const
Definition: TrajectoryStateOnSurface.h:80
TrajectoryStateOnSurface::globalMomentum
GlobalVector globalMomentum() const
Definition: TrajectoryStateOnSurface.h:86
gather_cfg.cout
cout
Definition: gather_cfg.py:144
GlobalTrajectoryParameters::vector
AlgebraicVector6 vector() const
Definition: GlobalTrajectoryParameters.h:120
void GlobalTrackerMuonAlignment::endJob ( void  )
overrideprivatevirtual

Reimplemented from edm::EDAnalyzer.

Definition at line 361 of file GlobalTrackerMuonAlignment.cc.

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

Referenced by o2olib.O2ORunMgr::executeJob().

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 }
mps_fire.i
i
Definition: mps_fire.py:338
AlgebraicSymMatrix33
ROOT::Math::SMatrix< double, 3, 3, ROOT::Math::MatRepSym< double, 3 > > AlgebraicSymMatrix33
Definition: AlgebraicROOTObjects.h:21
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
AlgebraicVector3
ROOT::Math::SVector< double, 3 > AlgebraicVector3
Definition: AlgebraicROOTObjects.h:12
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:79
GlobalTrackerMuonAlignment::writeGlPosRcd
void writeGlPosRcd(CLHEP::HepVector &d3)
Definition: GlobalTrackerMuonAlignment.cc:3322
MillePedeFileConverter_cfg.out
out
Definition: MillePedeFileConverter_cfg.py:31
gather_cfg.cout
cout
Definition: gather_cfg.py:144
void GlobalTrackerMuonAlignment::fitHist ( )

Definition at line 549 of file GlobalTrackerMuonAlignment.cc.

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

Referenced by endJob().

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 }
void GlobalTrackerMuonAlignment::gradientGlobal ( GlobalVector &  GRt,
GlobalVector &  GPt,
GlobalVector &  GRm,
GlobalVector &  GPm,
GlobalVector &  GNorm,
AlgebraicSymMatrix66 &  GCov 
)

Definition at line 2123 of file GlobalTrackerMuonAlignment.cc.

References patCaloMETCorrections_cff::A, a, TtFullHadDaughter::B, CLHEP_dot(), gather_cfg::cout, edmIntegrityCheck::d, debug_, delta, MillePedeFileConverter_cfg::e, Grad, Hess, mps_fire::i, cuy::ii, info(), AlCaHLTBitMon_ParallelJobs::p, reco::tau::disc::Pt(), alignCSCRings::r, alignCSCRings::s, mathSSE::sqrt(), w, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by analyzeTrackTrack(), and analyzeTrackTrajectory().

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
delta
dbl * delta
Definition: mlp_gen.cc:36
GlobalTrackerMuonAlignment::CLHEP_dot
double CLHEP_dot(const CLHEP::HepVector &a, const CLHEP::HepVector &b)
Definition: GlobalTrackerMuonAlignment.cc:148
info
static const TGPicture * info(bool iBackgroundIsBlack)
Definition: FWCollectionSummaryWidget.cc:170
mps_fire.i
i
Definition: mps_fire.py:338
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
w
const double w
Definition: UKUtility.cc:23
PV3DBase::y
T y() const
Definition: PV3DBase.h:63
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
PV3DBase::z
T z() const
Definition: PV3DBase.h:64
TtFullHadDaughter::B
static const std::string B
Definition: TtFullHadronicEvent.h:9
cuy.ii
ii
Definition: cuy.py:590
a
double a
Definition: hdecay.h:121
patCaloMETCorrections_cff.A
A
Definition: patCaloMETCorrections_cff.py:43
gather_cfg.cout
cout
Definition: gather_cfg.py:144
T
long double T
Definition: Basic3DVectorLD.h:66
PV3DBase::x
T x() const
Definition: PV3DBase.h:62
void GlobalTrackerMuonAlignment::gradientGlobalAlg ( GlobalVector &  Rt,
GlobalVector &  Pt,
GlobalVector &  Rm,
GlobalVector &  Nl,
AlgebraicSymMatrix66 &  Cm 
)

Definition at line 2032 of file GlobalTrackerMuonAlignment.cc.

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

Referenced by analyzeTrackTrack(), and analyzeTrackTrajectory().

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 }
mps_fire.i
i
Definition: mps_fire.py:338
PV3DBase::y
T y() const
Definition: PV3DBase.h:63
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
edm::eventsetup::Itr
RecordProviders::iterator Itr
Definition: EventSetupProvider.cc:354
AlgebraicSymMatrix33
ROOT::Math::SMatrix< double, 3, 3, ROOT::Math::MatRepSym< double, 3 > > AlgebraicSymMatrix33
Definition: AlgebraicROOTObjects.h:21
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
PV3DBase::z
T z() const
Definition: PV3DBase.h:64
AlgebraicVector3
ROOT::Math::SVector< double, 3 > AlgebraicVector3
Definition: AlgebraicROOTObjects.h:12
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:79
gather_cfg.cout
cout
Definition: gather_cfg.py:144
PV3DBase::x
T x() const
Definition: PV3DBase.h:62
AlgebraicMatrix33
ROOT::Math::SMatrix< double, 3, 3, ROOT::Math::MatRepStd< double, 3, 3 > > AlgebraicMatrix33
Definition: AlgebraicROOTObjects.h:41
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.

References CLHEP_dot(), gather_cfg::cout, debug_, delta, GradL, HessL, mps_fire::i, cuy::ii, info(), reco::tau::disc::Pt(), pfBoostedDoubleSVAK8TagInfos_cfi::R0, mathSSE::sqrt(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by analyzeTrackTrack(), and analyzeTrackTrajectory().

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
delta
dbl * delta
Definition: mlp_gen.cc:36
GlobalTrackerMuonAlignment::CLHEP_dot
double CLHEP_dot(const CLHEP::HepVector &a, const CLHEP::HepVector &b)
Definition: GlobalTrackerMuonAlignment.cc:148
pfBoostedDoubleSVAK8TagInfos_cfi.R0
R0
Definition: pfBoostedDoubleSVAK8TagInfos_cfi.py:8
info
static const TGPicture * info(bool iBackgroundIsBlack)
Definition: FWCollectionSummaryWidget.cc:170
mps_fire.i
i
Definition: mps_fire.py:338
PV3DBase::y
T y() const
Definition: PV3DBase.h:63
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
PV3DBase::z
T z() const
Definition: PV3DBase.h:64
cuy.ii
ii
Definition: cuy.py:590
gather_cfg.cout
cout
Definition: gather_cfg.py:144
PV3DBase::x
T x() const
Definition: PV3DBase.h:62
void GlobalTrackerMuonAlignment::misalignMuon ( GlobalVector &  GRm,
GlobalVector &  GPm,
GlobalVector &  Nl,
GlobalVector &  Rt,
GlobalVector &  Rm,
GlobalVector &  Pm 
)

Definition at line 2566 of file GlobalTrackerMuonAlignment.cc.

References patCaloMETCorrections_cff::A, a, TtFullHadDaughter::B, alignmentValidation::c1, CLHEP_dot(), funct::cos(), gather_cfg::cout, edmIntegrityCheck::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().

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
GlobalTrackerMuonAlignment::CLHEP_dot
double CLHEP_dot(const CLHEP::HepVector &a, const CLHEP::HepVector &b)
Definition: GlobalTrackerMuonAlignment.cc:148
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
PV3DBase::y
T y() const
Definition: PV3DBase.h:63
PV3DBase::z
T z() const
Definition: PV3DBase.h:64
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
TtFullHadDaughter::B
static const std::string B
Definition: TtFullHadronicEvent.h:9
a
double a
Definition: hdecay.h:121
patCaloMETCorrections_cff.A
A
Definition: patCaloMETCorrections_cff.py:43
gather_cfg.cout
cout
Definition: gather_cfg.py:144
T
long double T
Definition: Basic3DVectorLD.h:66
PV3DBase::x
T x() const
Definition: PV3DBase.h:62
GlobalVector
Global3DVector GlobalVector
Definition: GlobalVector.h:10
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.

References patCaloMETCorrections_cff::A, a, TtFullHadDaughter::B, alignmentValidation::c1, CLHEP_dot(), funct::cos(), gather_cfg::cout, edmIntegrityCheck::d, debug_, cuy::ii, TrajectoryStateOnSurface::localParameters(), MuGlAngle, MuGlShift, GloballyPositioned< T >::position(), LocalTrajectoryParameters::pzSign(), pfBoostedDoubleSVAK8TagInfos_cfi::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