CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 ()
 
- Public Member Functions inherited from edm::EDAnalyzer
void callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)
 
 EDAnalyzer ()
 
ModuleDescription const & moduleDescription () const
 
std::string workerType () const
 
virtual ~EDAnalyzer ()
 
- Public Member Functions inherited from edm::EDConsumerBase
 EDConsumerBase ()
 
ProductHolderIndexAndSkipBit indexFrom (EDGetToken, BranchType, TypeID const &) const
 
void itemsMayGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const
 
void itemsToGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const
 
std::vector
< ProductHolderIndexAndSkipBit >
const & 		itemsToGetFromEvent () const
 
void labelsForToken (EDGetToken iToken, Labels &oLabels) const
 
void modulesDependentUpon (const std::string &iProcessName, std::vector< const char * > &oModuleLabels) const
 
bool registeredToConsume (ProductHolderIndex, bool, BranchType) const
 
bool registeredToConsumeMany (TypeID const &, BranchType) const
 
void updateLookup (BranchType iBranchType, ProductHolderIndexHelper const &)
 
virtual ~EDConsumerBase ()
 

Private Member Functions

virtual void analyze (const edm::Event &, const edm::EventSetup &) override
 
virtual void beginJob () override
 
virtual 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
 
- Static Public Member Functions inherited from edm::EDAnalyzer
static const std::string & baseType ()
 
static void fillDescriptions (ConfigurationDescriptions &descriptions)
 
static void prevalidate (ConfigurationDescriptions &)
 
- 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 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 271 of file GlobalTrackerMuonAlignment.cc.

References cosmicMuonMode_, edm::hlt::Exception, and isolatedMuonMode_.

272  :trackTags_(iConfig.getParameter<edm::InputTag>("tracks")),
273  muonTags_(iConfig.getParameter<edm::InputTag>("muons")),
274  gmuonTags_(iConfig.getParameter<edm::InputTag>("gmuons")),
275  smuonTags_(iConfig.getParameter<edm::InputTag>("smuons")),
276  propagator_(iConfig.getParameter<string>("Propagator")),
277 
278  cosmicMuonMode_(iConfig.getParameter<bool>("cosmics")),
279  isolatedMuonMode_(iConfig.getParameter<bool>("isolated")),
280 
281  refitMuon_(iConfig.getParameter<bool>("refitmuon")),
282  refitTrack_(iConfig.getParameter<bool>("refittrack")),
283 
284  rootOutFile_(iConfig.getUntrackedParameter<string>("rootOutFile")),
285  txtOutFile_(iConfig.getUntrackedParameter<string>("txtOutFile")),
286  writeDB_(iConfig.getUntrackedParameter<bool>("writeDB")),
287  debug_(iConfig.getUntrackedParameter<bool>("debug")),
288  defineFitter(true)
289 {
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 ( )

Definition at line 302 of file GlobalTrackerMuonAlignment.cc.

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

Member Function Documentation

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

Implements edm::EDAnalyzer.

Definition at line 316 of file GlobalTrackerMuonAlignment.cc.

References analyzeTrackTrajectory().

317 {
318  //GlobalTrackerMuonAlignment::analyzeTrackTrack(iEvent, iSetup);
319  GlobalTrackerMuonAlignment::analyzeTrackTrajectory(iEvent, iSetup);
320 }
GlobalTrackerMuonAlignment::analyzeTrackTrajectory
void analyzeTrackTrajectory(const edm::Event &, const edm::EventSetup &)
Definition: GlobalTrackerMuonAlignment.cc:1242
void GlobalTrackerMuonAlignment::analyzeTrackTrack ( const edm::Event &  iEvent,
const edm::EventSetup &  iSetup 
)

Definition at line 565 of file GlobalTrackerMuonAlignment.cc.

References alongMomentum, edm::EventBase::bunchCrossing(), TrajectoryStateOnSurface::cartesianError(), gather_cfg::cout, Vector3DBase< T, FrameTag >::dot(), alignCSCRings::e, DetId::Ecal, epsilon, TrajectoryStateOnSurface::freeState(), edm::EventSetup::get(), edm::Event::getByLabel(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), gradientGlobal(), gradientGlobalAlg(), gradientLocal(), DetId::Hcal, timingPdfMaker::histo, i, info(), TrajectoryStateOnSurface::isValid(), j, TrajectoryStateOnSurface::localError(), TrajectoryStateOnSurface::localParameters(), m, PV3DBase< T, PVType, FrameType >::mag(), mag(), CartesianTrajectoryError::matrix(), LocalTrajectoryError::matrix(), misalignMuonL(), DetId::Muon, patZpeak::muons, Plane::normalVector(), oppositeToMomentum, PV3DBase< T, PVType, FrameType >::perp(), PI, GloballyPositioned< T >::position(), defaultRKPropagator::Product::propagator, reco::tau::disc::Pt(), DetId::rawId(), dt_dqm_sourceclient_common_cff::reco, GloballyPositioned< T >::rotation(), edm::Event::run(), SmartPropagator_cfi::SmartPropagator, mathSSE::sqrt(), SteppingHelixPropagator_cfi::SteppingHelixPropagator, TrajectoryStateOnSurface::surface(), Surface::tangentPlane(), DetId::Tracker, testEve_cfg::tracks, LocalTrajectoryParameters::vector(), x, PV3DBase< T, PVType, FrameType >::x(), TkRotation< T >::xx(), TkRotation< T >::xy(), TkRotation< T >::xz(), detailsBasic3DVector::y, PV3DBase< T, PVType, FrameType >::y(), TkRotation< T >::yx(), TkRotation< T >::yy(), TkRotation< T >::yz(), detailsBasic3DVector::z, PV3DBase< T, PVType, FrameType >::z(), TkRotation< T >::zx(), TkRotation< T >::zy(), and TkRotation< T >::zz().

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

References alongMomentum, edm::EventBase::bunchCrossing(), TrajectoryStateOnSurface::cartesianError(), Chi2MeasurementEstimatorESProducer_cfi::Chi2MeasurementEstimator, gather_cfg::cout, Vector3DBase< T, FrameTag >::dot(), alignCSCRings::e, DetId::Ecal, epsilon, TrajectoryStateOnSurface::freeState(), edm::EventSetup::get(), edm::Event::getByLabel(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), gradientGlobal(), gradientGlobalAlg(), gradientLocal(), DetId::Hcal, timingPdfMaker::histo, i, info(), TrajectoryStateOnSurface::isValid(), j, KFTrajectoryFitterESProducer_cfi::KFTrajectoryFitter, KFTrajectorySmootherESProducer_cfi::KFTrajectorySmoother, TrajectoryStateOnSurface::localError(), TrajectoryStateOnSurface::localParameters(), m, PV3DBase< T, PVType, FrameType >::mag(), mag(), CartesianTrajectoryError::matrix(), LocalTrajectoryError::matrix(), misalignMuonL(), DetId::Muon, muonFitter(), patZpeak::muons, Plane::normalVector(), oppositeToMomentum, PV3DBase< T, PVType, FrameType >::perp(), PI, GloballyPositioned< T >::position(), defaultRKPropagator::Product::propagator, reco::tau::disc::Pt(), DetId::rawId(), dt_dqm_sourceclient_common_cff::reco, GloballyPositioned< T >::rotation(), edm::Event::run(), SmartPropagator_cfi::SmartPropagator, mathSSE::sqrt(), SteppingHelixPropagator_cfi::SteppingHelixPropagator, TrajectoryStateOnSurface::surface(), Surface::tangentPlane(), DetId::Tracker, trackFitter(), testEve_cfg::tracks, MixedLayerPairs_cfi::TTRHBuilder, LocalTrajectoryParameters::vector(), x, PV3DBase< T, PVType, FrameType >::x(), TkRotation< T >::xx(), TkRotation< T >::xy(), TkRotation< T >::xz(), detailsBasic3DVector::y, PV3DBase< T, PVType, FrameType >::y(), TkRotation< T >::yx(), TkRotation< T >::yy(), TkRotation< T >::yz(), detailsBasic3DVector::z, PV3DBase< T, PVType, FrameType >::z(), TkRotation< T >::zx(), TkRotation< T >::zy(), and TkRotation< T >::zz().

Referenced by analyze().

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

Reimplemented from edm::EDAnalyzer.

Definition at line 323 of file GlobalTrackerMuonAlignment.cc.

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

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

Definition at line 460 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().

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

Definition at line 141 of file GlobalTrackerMuonAlignment.cc.

References a, and b.

142  {
143  return a(1)*b(1)+a(2)*b(2)+a(3)*b(3);
144  }
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 3131 of file GlobalTrackerMuonAlignment.cc.

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

3132 {
3133  std::cout<<" ------- "<<title<<" --------"<<std::endl;
3134  int nHit = 1;
3135  TransientTrackingRecHit::RecHitContainer recHit;
3136  for (trackingRecHit_iterator i=alongTr->recHitsBegin();i!=alongTr->recHitsEnd(); i++) {
3137  std::cout<<" Hit "<<nHit++<<" DetId "<<(*i)->geographicalId().det();
3138  if( (*i)->geographicalId().det() == DetId::Tracker ) std::cout<<" Tracker ";
3139  else if ( (*i)->geographicalId().det() == DetId::Muon ) std::cout<<" Muon ";
3140  else std::cout<<" Unknown ";
3141  if(!(*i)->isValid()) std::cout<<" not valid "<<std::endl;
3142  else std::cout<<std::endl;
3143  }
3144 }
i
int i
Definition: DBlmapReader.cc:9
TrackingRecHit::RecHitContainer
std::vector< ConstRecHitPointer > RecHitContainer
Definition: TrackingRecHit.h:39
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
void GlobalTrackerMuonAlignment::debugTrackHit ( const std::string  title,
reco::TransientTrack &  alongTr 
)

Definition at line 3114 of file GlobalTrackerMuonAlignment.cc.

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

3115 {
3116  std::cout<<" ------- "<<title<<" --------"<<std::endl;
3117  int nHit = 1;
3118  TransientTrackingRecHit::RecHitContainer recHit;
3119  for (trackingRecHit_iterator i=alongTr.recHitsBegin();i!=alongTr.recHitsEnd(); i++) {
3120  std::cout<<" Hit "<<nHit++<<" DetId "<<(*i)->geographicalId().det();
3121  if( (*i)->geographicalId().det() == DetId::Tracker ) std::cout<<" Tracker ";
3122  else if ( (*i)->geographicalId().det() == DetId::Muon ) std::cout<<" Muon ";
3123  else std::cout<<" Unknown ";
3124  if(!(*i)->isValid()) std::cout<<" not valid "<<std::endl;
3125  else std::cout<<std::endl;
3126  }
3127 }
i
int i
Definition: DBlmapReader.cc:9
TrackingRecHit::RecHitContainer
std::vector< ConstRecHitPointer > RecHitContainer
Definition: TrackingRecHit.h:39
reco::TransientTrack::recHitsEnd
trackingRecHit_iterator recHitsEnd() const
last iterator to RecHits
Definition: TransientTrack.h:126
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
reco::TransientTrack::recHitsBegin
trackingRecHit_iterator recHitsBegin() const
first iterator to RecHits
Definition: TransientTrack.h:124
void GlobalTrackerMuonAlignment::debugTrajectory ( const std::string  title,
Trajectory &  traj 
)

Definition at line 3211 of file GlobalTrackerMuonAlignment.cc.

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

3212 {
3213  std::cout<<"\n"<<" ...... "<<title<<" ...... "<<std::endl;
3214  if(!traj.isValid()) {std::cout<<" Not valid !!!!!!!! "<<std::endl; return;}
3215  std::cout<<" chi2/Nhit "<<traj.chiSquared()<<" / "<<traj.foundHits();
3216  if(traj.direction() == alongMomentum) std::cout<<" alongMomentum >>>>"<<std::endl;
3217  else std::cout<<" oppositeToMomentum <<<<"<<std::endl;
3218  this->debugTrajectorySOSv(" firstMeasurementTSOS ",traj.firstMeasurement().updatedState());
3219  //this->debugTrajectorySOSv(" firstMeasurementTSOS ",traj.firstMeasurement().predictedState());
3220  this->debugTrajectorySOSv(" lastMeasurementTSOS ",traj.lastMeasurement().updatedState());
3221  //this->debugTrajectorySOSv(" geom InnermostState", traj.geometricalInnermostState());
3222  std::cout<<" . . . . . . . . . . . . . . . . . . . . . . . . . . . . \n"<<std::endl;
3223 }
Trajectory::foundHits
int foundHits() const
Definition: Trajectory.h:234
GlobalTrackerMuonAlignment::debugTrajectorySOSv
void debugTrajectorySOSv(const std::string, TrajectoryStateOnSurface)
Definition: GlobalTrackerMuonAlignment.cc:3179
Trajectory::direction
PropagationDirection const & direction() const
Definition: Trajectory.cc:118
Trajectory::lastMeasurement
TrajectoryMeasurement const & lastMeasurement() const
Definition: Trajectory.h:181
Trajectory::isValid
bool isValid() const
Definition: Trajectory.h:269
Trajectory::firstMeasurement
TrajectoryMeasurement const & firstMeasurement() const
Definition: Trajectory.h:194
Trajectory::chiSquared
float chiSquared() const
Definition: Trajectory.h:252
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
TrajectoryMeasurement::updatedState
TrajectoryStateOnSurface const & updatedState() const
Definition: TrajectoryMeasurement.h:190
void GlobalTrackerMuonAlignment::debugTrajectorySOS ( const std::string  title,
TrajectoryStateOnSurface &  trajSOS 
)

Definition at line 3147 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().

3149 {
3150  std::cout<<" --- "<<title<<" --- "<<std::endl;
3151  if(!trajSOS.isValid()) {std::cout<<" Not valid !!!! "<<std::endl; return;}
3152  std::cout<<" R |p| "<<trajSOS.globalPosition().perp()<<" "
3153  <<trajSOS.globalMomentum().mag()<<" charge "<<trajSOS.charge()<<std::endl;
3154  std::cout<<" x p "<<trajSOS.globalParameters().vector()(0)<<" "
3155  <<trajSOS.globalParameters().vector()(1)<<" "
3156  <<trajSOS.globalParameters().vector()(2)<<" "
3157  <<trajSOS.globalParameters().vector()(3)<<" "
3158  <<trajSOS.globalParameters().vector()(4)<<" "
3159  <<trajSOS.globalParameters().vector()(5)<<std::endl;
3160  std::cout<<" +/- "<<sqrt(trajSOS.cartesianError().matrix()(0,0))<<" "
3161  <<sqrt(trajSOS.cartesianError().matrix()(1,1))<<" "
3162  <<sqrt(trajSOS.cartesianError().matrix()(2,2))<<" "
3163  <<sqrt(trajSOS.cartesianError().matrix()(3,3))<<" "
3164  <<sqrt(trajSOS.cartesianError().matrix()(4,4))<<" "
3165  <<sqrt(trajSOS.cartesianError().matrix()(5,5))<<std::endl;
3166  std::cout<<" q/p dxy/dz xy "<<trajSOS.localParameters().vector()(0)<<" "
3167  <<trajSOS.localParameters().vector()(1)<<" "
3168  <<trajSOS.localParameters().vector()(2)<<" "
3169  <<trajSOS.localParameters().vector()(3)<<" "
3170  <<trajSOS.localParameters().vector()(4)<<std::endl;
3171  std::cout<<" +/- error "<<sqrt(trajSOS.localError().matrix()(0,0))<<" "
3172  <<sqrt(trajSOS.localError().matrix()(1,1))<<" "
3173  <<sqrt(trajSOS.localError().matrix()(2,2))<<" "
3174  <<sqrt(trajSOS.localError().matrix()(3,3))<<" "
3175  <<sqrt(trajSOS.localError().matrix()(4,4))<<" "<<std::endl;
3176 }
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:72
TrajectoryStateOnSurface::localParameters
const LocalTrajectoryParameters & localParameters() const
Definition: TrajectoryStateOnSurface.h:113
TrajectoryStateOnSurface::cartesianError
const CartesianTrajectoryError cartesianError() const
Definition: TrajectoryStateOnSurface.h:107
TrajectoryStateOnSurface::globalPosition
GlobalPoint globalPosition() const
Definition: TrajectoryStateOnSurface.h:89
LocalTrajectoryParameters::vector
AlgebraicVector5 vector() const
Definition: LocalTrajectoryParameters.h:120
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:67
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
LocalTrajectoryError::matrix
const AlgebraicSymMatrix55 & matrix() const
Definition: LocalTrajectoryError.h:61
TrajectoryStateOnSurface::localError
const LocalTrajectoryError & localError() const
Definition: TrajectoryStateOnSurface.h:125
CartesianTrajectoryError::matrix
const AlgebraicSymMatrix66 & matrix() const
Definition: CartesianTrajectoryError.h:29
TrajectoryStateOnSurface::globalParameters
const GlobalTrajectoryParameters & globalParameters() const
Definition: TrajectoryStateOnSurface.h:86
TrajectoryStateOnSurface::globalMomentum
GlobalVector globalMomentum() const
Definition: TrajectoryStateOnSurface.h:92
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
GlobalTrajectoryParameters::vector
AlgebraicVector6 vector() const
Definition: GlobalTrajectoryParameters.h:109
void GlobalTrackerMuonAlignment::debugTrajectorySOSv ( const std::string  title,
TrajectoryStateOnSurface  trajSOS 
)

Definition at line 3179 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().

3181 {
3182  std::cout<<" --- "<<title<<" --- "<<std::endl;
3183  if(!trajSOS.isValid()) {std::cout<<" Not valid !!!! "<<std::endl; return;}
3184  std::cout<<" R |p| "<<trajSOS.globalPosition().perp()<<" "
3185  <<trajSOS.globalMomentum().mag()<<" charge "<<trajSOS.charge()<<std::endl;
3186  std::cout<<" x p "<<trajSOS.globalParameters().vector()(0)<<" "
3187  <<trajSOS.globalParameters().vector()(1)<<" "
3188  <<trajSOS.globalParameters().vector()(2)<<" "
3189  <<trajSOS.globalParameters().vector()(3)<<" "
3190  <<trajSOS.globalParameters().vector()(4)<<" "
3191  <<trajSOS.globalParameters().vector()(5)<<std::endl;
3192  std::cout<<" +/- "<<sqrt(trajSOS.cartesianError().matrix()(0,0))<<" "
3193  <<sqrt(trajSOS.cartesianError().matrix()(1,1))<<" "
3194  <<sqrt(trajSOS.cartesianError().matrix()(2,2))<<" "
3195  <<sqrt(trajSOS.cartesianError().matrix()(3,3))<<" "
3196  <<sqrt(trajSOS.cartesianError().matrix()(4,4))<<" "
3197  <<sqrt(trajSOS.cartesianError().matrix()(5,5))<<std::endl;
3198  std::cout<<" q/p dxy/dz xy "<<trajSOS.localParameters().vector()(0)<<" "
3199  <<trajSOS.localParameters().vector()(1)<<" "
3200  <<trajSOS.localParameters().vector()(2)<<" "
3201  <<trajSOS.localParameters().vector()(3)<<" "
3202  <<trajSOS.localParameters().vector()(4)<<std::endl;
3203  std::cout<<" +/- error "<<sqrt(trajSOS.localError().matrix()(0,0))<<" "
3204  <<sqrt(trajSOS.localError().matrix()(1,1))<<" "
3205  <<sqrt(trajSOS.localError().matrix()(2,2))<<" "
3206  <<sqrt(trajSOS.localError().matrix()(3,3))<<" "
3207  <<sqrt(trajSOS.localError().matrix()(4,4))<<" "<<std::endl;
3208 }
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:72
TrajectoryStateOnSurface::localParameters
const LocalTrajectoryParameters & localParameters() const
Definition: TrajectoryStateOnSurface.h:113
TrajectoryStateOnSurface::cartesianError
const CartesianTrajectoryError cartesianError() const
Definition: TrajectoryStateOnSurface.h:107
TrajectoryStateOnSurface::globalPosition
GlobalPoint globalPosition() const
Definition: TrajectoryStateOnSurface.h:89
LocalTrajectoryParameters::vector
AlgebraicVector5 vector() const
Definition: LocalTrajectoryParameters.h:120
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:67
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
LocalTrajectoryError::matrix
const AlgebraicSymMatrix55 & matrix() const
Definition: LocalTrajectoryError.h:61
TrajectoryStateOnSurface::localError
const LocalTrajectoryError & localError() const
Definition: TrajectoryStateOnSurface.h:125
CartesianTrajectoryError::matrix
const AlgebraicSymMatrix66 & matrix() const
Definition: CartesianTrajectoryError.h:29
TrajectoryStateOnSurface::globalParameters
const GlobalTrajectoryParameters & globalParameters() const
Definition: TrajectoryStateOnSurface.h:86
TrajectoryStateOnSurface::globalMomentum
GlobalVector globalMomentum() const
Definition: TrajectoryStateOnSurface.h:92
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
GlobalTrajectoryParameters::vector
AlgebraicVector6 vector() const
Definition: GlobalTrajectoryParameters.h:109
void GlobalTrackerMuonAlignment::endJob ( void  )
overrideprivatevirtual

Reimplemented from edm::EDAnalyzer.

Definition at line 362 of file GlobalTrackerMuonAlignment.cc.

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

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

Definition at line 539 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().

539  {
540 
541  histo7->Fit("gaus","Q");
542 
543  histo12->Fit("gaus","Q");
544  histo13->Fit("gaus","Q");
545  histo14->Fit("gaus","Q");
546  histo15->Fit("gaus","Q");
547  histo16->Fit("gaus","Q");
548  histo17->Fit("gaus","Q");
549 
550  histo21->Fit("gaus","Q");
551  histo22->Fit("gaus","Q");
552  histo23->Fit("gaus","Q");
553  histo24->Fit("gaus","Q");
554  histo25->Fit("gaus","Q");
555  histo26->Fit("gaus","Q");
556 
557  histo32->Fit("gaus","Q");
558  histo33->Fit("gaus","Q");
559  histo34->Fit("gaus","Q");
560  histo35->Fit("gaus","Q");
561 }
void GlobalTrackerMuonAlignment::gradientGlobal ( GlobalVector &  GRt,
GlobalVector &  GPt,
GlobalVector &  GRm,
GlobalVector &  GPm,
GlobalVector &  GNorm,
AlgebraicSymMatrix66 &  GCov 
)

Definition at line 2091 of file GlobalTrackerMuonAlignment.cc.

References funct::A, a, gather_cfg::cout, delta, alignCSCRings::e, i, cuy::ii, info(), j, 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().

2094 {
2095  // we search for 6D global correction vector (d, a), where
2096  // 3D vector of shihts d
2097  // 3D vector of rotation angles a
2098 
2099  //bool alarm = true;
2100  bool alarm = false;
2101  bool info = false;
2102 
2103  int Nd = 6; // dimension of vector of alignment pararmeters, d
2104 
2105  //double PtMom = GPt.mag();
2106  CLHEP::HepSymMatrix w(Nd,0);
2107  for (int i=1; i <= Nd; i++)
2108  for (int j=1; j <= Nd; j++){
2109  if(j <= i ) w(i,j) = GCov(i-1, j-1);
2110  //if(i >= 3) w(i,j) /= PtMom;
2111  //if(j >= 3) w(i,j) /= PtMom;
2112  if((i == j) && (i<=3) && (GCov(i-1, j-1) < 1.e-20)) w(i,j) = 1.e20; // w=0
2113  if(i != j) w(i,j) = 0.; // use diaginal elements
2114  }
2115 
2116  //GPt /= GPt.mag();
2117  //GPm /= GPm.mag(); // end of transform
2118 
2119  CLHEP::HepVector V(Nd), Rt(3), Pt(3), Rm(3), Pm(3), Norm(3);
2120  Rt(1) = GRt.x(); Rt(2) = GRt.y(); Rt(3) = GRt.z();
2121  Pt(1) = GPt.x(); Pt(2) = GPt.y(); Pt(3) = GPt.z();
2122  Rm(1) = GRm.x(); Rm(2) = GRm.y(); Rm(3) = GRm.z();
2123  Pm(1) = GPm.x(); Pm(2) = GPm.y(); Pm(3) = GPm.z();
2124  Norm(1) = GNorm.x(); Norm(2) = GNorm.y(); Norm(3) = GNorm.z();
2125 
2126  V = dsum(Rm - Rt, Pm - Pt) ; //std::cout<<" V "<<V.T()<<std::endl;
2127 
2128  //double PmN = CLHEP::dot(Pm, Norm);
2129  double PmN = CLHEP_dot(Pm, Norm);
2130 
2131  CLHEP::HepMatrix Jac(Nd,Nd,0);
2132  for (int i=1; i <= 3; i++)
2133  for (int j=1; j <= 3; j++){
2134  Jac(i,j) = Pm(i)*Norm(j) / PmN;
2135  if(i == j ) Jac(i,j) -= 1.;
2136  }
2137 
2138  // dp/da
2139  Jac(4,4) = 0.; // dpx/dax
2140  Jac(5,4) = -Pm(3); // dpy/dax
2141  Jac(6,4) = Pm(2); // dpz/dax
2142  Jac(4,5) = Pm(3); // dpx/day
2143  Jac(5,5) = 0.; // dpy/day
2144  Jac(6,5) = -Pm(1); // dpz/day
2145  Jac(4,6) = -Pm(2); // dpx/daz
2146  Jac(5,6) = Pm(1); // dpy/daz
2147  Jac(6,6) = 0.; // dpz/daz
2148 
2149  CLHEP::HepVector dsda(3);
2150  dsda(1) = (Norm(2)*Rm(3) - Norm(3)*Rm(2)) / PmN;
2151  dsda(2) = (Norm(3)*Rm(1) - Norm(1)*Rm(3)) / PmN;
2152  dsda(3) = (Norm(1)*Rm(2) - Norm(2)*Rm(1)) / PmN;
2153 
2154  // dr/da
2155  Jac(1,4) = Pm(1)*dsda(1); // drx/dax
2156  Jac(2,4) = -Rm(3) + Pm(2)*dsda(1); // dry/dax
2157  Jac(3,4) = Rm(2) + Pm(3)*dsda(1); // drz/dax
2158 
2159  Jac(1,5) = Rm(3) + Pm(1)*dsda(2); // drx/day
2160  Jac(2,5) = Pm(2)*dsda(2); // dry/day
2161  Jac(3,5) = -Rm(1) + Pm(3)*dsda(2); // drz/day
2162 
2163  Jac(1,6) = -Rm(2) + Pm(1)*dsda(3); // drx/daz
2164  Jac(2,6) = Rm(1) + Pm(2)*dsda(3); // dry/daz
2165  Jac(3,6) = Pm(3)*dsda(3); // drz/daz
2166 
2167  CLHEP::HepSymMatrix W(Nd,0);
2168  int ierr;
2169  W = w.inverse(ierr);
2170  if(ierr != 0) { if(alarm)
2171  std::cout<<" gradientGlobal: inversion of matrix w fail "<<std::endl;
2172  return;
2173  }
2174 
2175  CLHEP::HepMatrix W_Jac(Nd,Nd,0);
2176  W_Jac = Jac.T() * W;
2177 
2178  CLHEP::HepVector grad3(Nd);
2179  grad3 = W_Jac * V;
2180 
2181  CLHEP::HepMatrix hess3(Nd,Nd);
2182  hess3 = Jac.T() * W * Jac;
2183  //hess3(4,4) = 1.e-10; hess3(5,5) = 1.e-10; hess3(6,6) = 1.e-10; //????????????????
2184 
2185  Grad += grad3;
2186  Hess += hess3;
2187 
2188  CLHEP::HepVector d3I = CLHEP::solve(Hess, -Grad);
2189  int iEr3I;
2190  CLHEP::HepMatrix Errd3I = Hess.inverse(iEr3I);
2191  if( iEr3I != 0) { if(alarm)
2192  std::cout<<" gradientGlobal error inverse Hess matrix !!!!!!!!!!!"<<std::endl;
2193  }
2194 
2195  if(info || debug_){
2196  std::cout<<" dG "<<d3I(1)<<" "<<d3I(2)<<" "<<d3I(3)<<" "<<d3I(4)<<" "
2197  <<d3I(5)<<" "<<d3I(6)<<std::endl;;
2198  std::cout<<" +- "<<sqrt(Errd3I(1,1))<<" "<<sqrt(Errd3I(2,2))<<" "<<sqrt(Errd3I(3,3))
2199  <<" "<<sqrt(Errd3I(4,4))<<" "<<sqrt(Errd3I(5,5))<<" "<<sqrt(Errd3I(6,6))
2200  <<std::endl;
2201  }
2202 
2203 #ifdef CHECK_OF_DERIVATIVES
2204  // -------------------- check of derivatives
2205 
2206  CLHEP::HepVector d(3,0), a(3,0);
2207  CLHEP::HepMatrix T(3,3,1);
2208 
2209  CLHEP::HepMatrix Ti = T.T();
2210  //double A = CLHEP::dot(Ti*Pm, Norm);
2211  //double B = CLHEP::dot((Rt -Ti*Rm + Ti*d), Norm);
2212  double A = CLHEP_dot(Ti*Pm, Norm);
2213  double B = CLHEP_dot((Rt -Ti*Rm + Ti*d), Norm);
2214  double s0 = B / A;
2215 
2216  CLHEP::HepVector r0(3,0), p0(3,0);
2217  r0 = Ti*Rm - Ti*d + s0*(Ti*Pm) - Rt;
2218  p0 = Ti*Pm - Pt;
2219 
2220  double delta = 0.0001;
2221 
2222  int ii = 3; d(ii) += delta; // d
2223  //T(2,3) += delta; T(3,2) -= delta; int ii = 1; // a1
2224  //T(3,1) += delta; T(1,3) -= delta; int ii = 2; // a2
2225  //T(1,2) += delta; T(2,1) -= delta; int ii = 3; // a2
2226  Ti = T.T();
2227  //A = CLHEP::dot(Ti*Pm, Norm);
2228  //B = CLHEP::dot((Rt -Ti*Rm + Ti*d), Norm);
2229  A = CLHEP_dot(Ti*Pm, Norm);
2230  B = CLHEP_dot((Rt -Ti*Rm + Ti*d), Norm);
2231  double s = B / A;
2232 
2233  CLHEP::HepVector r(3,0), p(3,0);
2234  r = Ti*Rm - Ti*d + s*(Ti*Pm) - Rt;
2235  p = Ti*Pm - Pt;
2236 
2237  std::cout<<" s0 s num dsda("<<ii<<") "<<s0<<" "<<s<<" "
2238  <<(s-s0)/delta<<" "<<dsda(ii)<<endl;
2239  // d(r,p) / d shift
2240  std::cout<<" -- An d(r,p)/d("<<ii<<") "<<Jac(1,ii)<<" "<<Jac(2,ii)<<" "
2241  <<Jac(3,ii)<<" "<<Jac(4,ii)<<" "<<Jac(5,ii)<<" "
2242  <<Jac(6,ii)<<std::endl;
2243  std::cout<<" Nu d(r,p)/d("<<ii<<") "<<(r(1)-r0(1))/delta<<" "
2244  <<(r(2)-r0(2))/delta<<" "<<(r(3)-r0(3))/delta<<" "
2245  <<(p(1)-p0(1))/delta<<" "<<(p(2)-p0(2))/delta<<" "
2246  <<(p(3)-p0(3))/delta<<std::endl;
2247  // d(r,p) / d angle
2248  std::cout<<" -- An d(r,p)/a("<<ii<<") "<<Jac(1,ii+3)<<" "<<Jac(2,ii+3)<<" "
2249  <<Jac(3,ii+3)<<" "<<Jac(4,ii+3)<<" "<<Jac(5,ii+3)<<" "
2250  <<Jac(6,ii+3)<<std::endl;
2251  std::cout<<" Nu d(r,p)/a("<<ii<<") "<<(r(1)-r0(1))/delta<<" "
2252  <<(r(2)-r0(2))/delta<<" "<<(r(3)-r0(3))/delta<<" "
2253  <<(p(1)-p0(1))/delta<<" "<<(p(2)-p0(2))/delta<<" "
2254  <<(p(3)-p0(3))/delta<<std::endl;
2255  // ----------------------------- end of check
2256 #endif
2257 
2258  return;
2259 } // 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:141
i
int i
Definition: DBlmapReader.cc:9
info
static const TGPicture * info(bool iBackgroundIsBlack)
Definition: FWCollectionSummaryWidget.cc:170
PV3DBase::y
T y() const
Definition: PV3DBase.h:63
cuy.ii
int ii
Definition: cuy.py:588
double_binary
Definition: DDStreamer.cc:227
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
PV3DBase::z
T z() const
Definition: PV3DBase.h:64
j
int j
Definition: DBlmapReader.cc:9
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
a
double a
Definition: hdecay.h:121
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
w
T w() const
Definition: extBasic3DVector.h:245
T
long double T
Definition: Basic3DVectorLD.h:57
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 2012 of file GlobalTrackerMuonAlignment.cc.

References gather_cfg::cout, PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi::dR, alignCSCRings::e, i, j, 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().

2015 {
2016  // ---------------------------- Calculate Information matrix and Gfree vector
2017  // Information == Hessian , Gfree == gradient of the objective function
2018 
2019  AlgebraicMatrix33 Jac;
2020  AlgebraicVector3 Wi, R_m, R_t, P_t, Norm, dR;
2021 
2022  R_m(0) = Rm.x(); R_m(1) = Rm.y(); R_m(2) = Rm.z();
2023  R_t(0) = Rt.x(); R_t(1) = Rt.y(); R_t(2) = Rt.z();
2024  P_t(0) = Pt.x(); P_t(1) = Pt.y(); P_t(2) = Pt.z();
2025  Norm(0) = Nl.x(); Norm(1) = Nl.y(); Norm(2) = Nl.z();
2026 
2027  for(int i=0; i<=2; i++){
2028  if(Cm(i,i) > 1.e-20)
2029  Wi(i) = 1./Cm(i,i);
2030  else Wi(i) = 1.e-10;
2031  dR(i) = R_m(i) - R_t(i);
2032  }
2033 
2034  float PtN = P_t(0)*Norm(0) + P_t(1)*Norm(1) + P_t(2)*Norm(2);
2035 
2036  Jac(0,0) = 1. - P_t(0)*Norm(0)/PtN;
2037  Jac(0,1) = - P_t(0)*Norm(1)/PtN;
2038  Jac(0,2) = - P_t(0)*Norm(2)/PtN;
2039 
2040  Jac(1,0) = - P_t(1)*Norm(0)/PtN;
2041  Jac(1,1) = 1. - P_t(1)*Norm(1)/PtN;
2042  Jac(1,2) = - P_t(1)*Norm(2)/PtN;
2043 
2044  Jac(2,0) = - P_t(2)*Norm(0)/PtN;
2045  Jac(2,1) = - P_t(2)*Norm(1)/PtN;
2046  Jac(2,2) = 1. - P_t(2)*Norm(2)/PtN;
2047 
2048  AlgebraicSymMatrix33 Itr;
2049 
2050  for(int i=0; i<=2; i++)
2051  for(int j=0; j<=2; j++){
2052  if(j < i) continue;
2053  Itr(i,j) = 0.;
2054  //std::cout<<" ij "<<i<<" "<<j<<std::endl;
2055  for(int k=0; k<=2; k++){
2056  Itr(i,j) += Jac(k,i)*Wi(k)*Jac(k,j);}}
2057 
2058  for(int i=0; i<=2; i++)
2059  for(int j=0; j<=2; j++){
2060  if(j < i) continue;
2061  Inf(i,j) += Itr(i,j);}
2062 
2063  AlgebraicVector3 Gtr(0., 0., 0.);
2064  for(int i=0; i<=2; i++)
2065  for(int k=0; k<=2; k++) Gtr(i) += dR(k)*Wi(k)*Jac(k,i);
2066  for(int i=0; i<=2; i++) Gfr(i) += Gtr(i);
2067 
2068  if(debug_){
2069  std::cout<<" Wi "<<Wi<<std::endl;
2070  std::cout<<" N "<<Norm<<std::endl;
2071  std::cout<<" P_t "<<P_t<<std::endl;
2072  std::cout<<" (Pt*N) "<<PtN<<std::endl;
2073  std::cout<<" dR "<<dR<<std::endl;
2074  std::cout<<" +/- "<<1./sqrt(Wi(0))<<" "<<1./sqrt(Wi(1))<<" "<<1./sqrt(Wi(2))
2075  <<" "<<std::endl;
2076  std::cout<<" Jacobian dr/ddx "<<std::endl;
2077  std::cout<<Jac<<std::endl;
2078  std::cout<<" G-- "<<Gtr<<std::endl;
2079  std::cout<<" Itrack "<<std::endl;
2080  std::cout<<Itr<<std::endl;
2081  std::cout<<" Gfr "<<Gfr<<std::endl;
2082  std::cout<<" -- Inf --"<<std::endl;
2083  std::cout<<Inf<<std::endl;
2084  }
2085 
2086  return;
2087 }
i
int i
Definition: DBlmapReader.cc:9
PV3DBase::y
T y() const
Definition: PV3DBase.h:63
AlgebraicSymMatrix33
ROOT::Math::SMatrix< double, 3, 3, ROOT::Math::MatRepSym< double, 3 > > AlgebraicSymMatrix33
Definition: AlgebraicROOTObjects.h:29
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
PV3DBase::z
T z() const
Definition: PV3DBase.h:64
AlgebraicVector3
ROOT::Math::SVector< double, 3 > AlgebraicVector3
Definition: AlgebraicROOTObjects.h:20
j
int j
Definition: DBlmapReader.cc:9
PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi.dR
tuple dR
Definition: PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi.py:19
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:79
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
PV3DBase::x
T x() const
Definition: PV3DBase.h:62
edm::eventsetup::Itr
Providers::iterator Itr
Definition: EventSetupProvider.cc:383
AlgebraicMatrix33
ROOT::Math::SMatrix< double, 3, 3, ROOT::Math::MatRepStd< double, 3, 3 > > AlgebraicMatrix33
Definition: AlgebraicROOTObjects.h:49
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 2264 of file GlobalTrackerMuonAlignment.cc.

References gather_cfg::cout, delta, i, cuy::ii, info(), j, reco::tau::disc::Pt(), mathSSE::sqrt(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by analyzeTrackTrack(), and analyzeTrackTrajectory().

2271 {
2272  // we search for 6D global correction vector (d, a), where
2273  // 3D vector of shihts d
2274  // 3D vector of rotation angles a
2275 
2276  bool alarm = true;
2277  //bool alarm = false;
2278  bool info = false;
2279 
2280  if(debug_)
2281  std::cout<<" gradientLocal "<<std::endl;
2282 
2283  /*
2284  const Surface& refSurface = tsosMuon.surface();
2285 
2286  CLHEP::HepMatrix rotLoc (3,3,0);
2287  rotLoc(1,1) = refSurface.rotation().xx();
2288  rotLoc(1,2) = refSurface.rotation().xy();
2289  rotLoc(1,3) = refSurface.rotation().xz();
2290 
2291  rotLoc(2,1) = refSurface.rotation().yx();
2292  rotLoc(2,2) = refSurface.rotation().yy();
2293  rotLoc(2,3) = refSurface.rotation().yz();
2294 
2295  rotLoc(3,1) = refSurface.rotation().zx();
2296  rotLoc(3,2) = refSurface.rotation().zy();
2297  rotLoc(3,3) = refSurface.rotation().zz();
2298  */
2299 
2300  CLHEP::HepVector Rt(3), Pt(3), Rm(3), Pm(3), Norm(3);
2301  Rt(1) = GRt.x(); Rt(2) = GRt.y(); Rt(3) = GRt.z();
2302  Pt(1) = GPt.x(); Pt(2) = GPt.y(); Pt(3) = GPt.z();
2303  Rm(1) = GRm.x(); Rm(2) = GRm.y(); Rm(3) = GRm.z();
2304  Pm(1) = GPm.x(); Pm(2) = GPm.y(); Pm(3) = GPm.z();
2305  Norm(1) = GNorm.x(); Norm(2) = GNorm.y(); Norm(3) = GNorm.z();
2306 
2307  CLHEP::HepVector V(4), Rml(3), Pml(3), Rtl(3), Ptl(3);
2308 
2309  /*
2310  R0(1) = refSurface.position().x();
2311  R0(2) = refSurface.position().y();
2312  R0(3) = refSurface.position().z();
2313  */
2314 
2315  Rml = rotLoc * (Rm - R0);
2316  Rtl = rotLoc * (Rt - R0);
2317  Pml = rotLoc * Pm;
2318  Ptl = rotLoc * Pt;
2319 
2320  V(1) = LPRm(0) - Ptl(1)/Ptl(3);
2321  V(2) = LPRm(1) - Ptl(2)/Ptl(3);
2322  V(3) = LPRm(2) - Rtl(1);
2323  V(4) = LPRm(3) - Rtl(2);
2324 
2325  /*
2326  CLHEP::HepSymMatrix Cov(4,0), W(4,0);
2327  for(int i=1; i<=4; i++)
2328  for(int j=1; j<=i; j++){
2329  Cov(i,j) = (tsosMuon.localError().matrix()
2330  + tsosTrack.localError().matrix())(i,j);
2331  //if(i != j) Cov(i,j) = 0.;
2332  //if((i == j) && ((i==1) || (i==2))) Cov(i,j) = 100.;
2333  //if((i == j) && ((i==3) || (i==4))) Cov(i,j) = 10000.;
2334  }
2335  W = Cov;
2336  */
2337 
2338  CLHEP::HepSymMatrix W = covLoc;
2339 
2340  int ierr;
2341  W.invert(ierr);
2342  if(ierr != 0) { if(alarm)
2343  std::cout<<" gradientLocal: inversion of matrix W fail "<<std::endl;
2344  return;
2345  }
2346 
2347  // JacobianCartesianToLocal
2348 
2349  //AlgebraicMatrix56 jacobian // differ from calculation above
2350  //= JacobianCartesianToLocal::JacobianCartesianToLocal
2351  //(refSurface, tsosTrack.localParameters()).jacobian();
2352  //for(int i=1; i<=4; i++) for(int j=1; j<=6; j++){
2353  //int j1 = j - 1; JacToLoc(i,j) = jacobian(i, j1);}
2354 
2355  CLHEP::HepMatrix JacToLoc(4,6,0);
2356  for(int i=1; i<=2; i++)
2357  for(int j=1; j<=3; j++){
2358  JacToLoc(i,j+3) = (rotLoc(i,j) - rotLoc(3,j)*Pml(i)/Pml(3))/Pml(3);
2359  JacToLoc(i+2,j) = rotLoc(i,j);
2360  }
2361 
2362  // JacobianCorrectionsToCartesian
2363  //double PmN = CLHEP::dot(Pm, Norm);
2364  double PmN = CLHEP_dot(Pm, Norm);
2365 
2366  CLHEP::HepMatrix Jac(6,6,0);
2367  for (int i=1; i <= 3; i++)
2368  for (int j=1; j <= 3; j++){
2369  Jac(i,j) = Pm(i)*Norm(j) / PmN;
2370  if(i == j ) Jac(i,j) -= 1.;
2371  }
2372 
2373  // dp/da
2374  Jac(4,4) = 0.; // dpx/dax
2375  Jac(5,4) = -Pm(3); // dpy/dax
2376  Jac(6,4) = Pm(2); // dpz/dax
2377  Jac(4,5) = Pm(3); // dpx/day
2378  Jac(5,5) = 0.; // dpy/day
2379  Jac(6,5) = -Pm(1); // dpz/day
2380  Jac(4,6) = -Pm(2); // dpx/daz
2381  Jac(5,6) = Pm(1); // dpy/daz
2382  Jac(6,6) = 0.; // dpz/daz
2383 
2384  CLHEP::HepVector dsda(3);
2385  dsda(1) = (Norm(2)*Rm(3) - Norm(3)*Rm(2)) / PmN;
2386  dsda(2) = (Norm(3)*Rm(1) - Norm(1)*Rm(3)) / PmN;
2387  dsda(3) = (Norm(1)*Rm(2) - Norm(2)*Rm(1)) / PmN;
2388 
2389  // dr/da
2390  Jac(1,4) = Pm(1)*dsda(1); // drx/dax
2391  Jac(2,4) = -Rm(3) + Pm(2)*dsda(1); // dry/dax
2392  Jac(3,4) = Rm(2) + Pm(3)*dsda(1); // drz/dax
2393 
2394  Jac(1,5) = Rm(3) + Pm(1)*dsda(2); // drx/day
2395  Jac(2,5) = Pm(2)*dsda(2); // dry/day
2396  Jac(3,5) = -Rm(1) + Pm(3)*dsda(2); // drz/day
2397 
2398  Jac(1,6) = -Rm(2) + Pm(1)*dsda(3); // drx/daz
2399  Jac(2,6) = Rm(1) + Pm(2)*dsda(3); // dry/daz
2400  Jac(3,6) = Pm(3)*dsda(3); // drz/daz
2401 
2402  // JacobianCorrectionToLocal
2403  CLHEP::HepMatrix JacCorLoc(4,6,0);
2404  JacCorLoc = JacToLoc * Jac;
2405 
2406  // gradient and Hessian
2407  CLHEP::HepMatrix W_Jac(6,4,0);
2408  W_Jac = JacCorLoc.T() * W;
2409 
2410  CLHEP::HepVector gradL(6);
2411  gradL = W_Jac * V;
2412 
2413  CLHEP::HepMatrix hessL(6,6);
2414  hessL = JacCorLoc.T() * W * JacCorLoc;
2415 
2416  GradL += gradL;
2417  HessL += hessL;
2418 
2419  CLHEP::HepVector dLI = CLHEP::solve(HessL, -GradL);
2420  int iErI;
2421  CLHEP::HepMatrix ErrdLI = HessL.inverse(iErI);
2422  if( iErI != 0) { if(alarm)
2423  std::cout<<" gradLocal Error inverse Hess matrix !!!!!!!!!!!"<<std::endl;
2424  }
2425 
2426  if(info || debug_){
2427  std::cout<<" dL "<<dLI(1)<<" "<<dLI(2)<<" "<<dLI(3)<<" "<<dLI(4)<<" "
2428  <<dLI(5)<<" "<<dLI(6)<<std::endl;;
2429  std::cout<<" +- "<<sqrt(ErrdLI(1,1))<<" "<<sqrt(ErrdLI(2,2))<<" "<<sqrt(ErrdLI(3,3))
2430  <<" "<<sqrt(ErrdLI(4,4))<<" "<<sqrt(ErrdLI(5,5))<<" "<<sqrt(ErrdLI(6,6))
2431  <<std::endl;
2432  }
2433 
2434  if(debug_){
2435  //std::cout<<" dV(da3) {"<<-JacCorLoc(1,6)*0.002<<" "<<-JacCorLoc(2,6)*0.002<<" "
2436  // <<-JacCorLoc(3,6)*0.002<<" "<<-JacCorLoc(4,6)*0.002<<"}"<<std::endl;
2437  //std::cout<<" JacCLo {"<<JacCorLoc(1,6)<<" "<<JacCorLoc(2,6)<<" "
2438  // <<JacCorLoc(3,6)<<" "<<JacCorLoc(4,6)<<"}"<<std::endl;
2439  //std::cout<<"Jpx/yx {"<<Jac(4,6)/Pm(3)<<" "<<Jac(5,6)/Pm(3)<<" "
2440  // <<Jac(1,6)<<" "<<Jac(2,6)<<"}"<<std::endl;
2441  //std::cout<<"Jac(,a3){"<<Jac(1,6)<<" "<<Jac(2,6)<<" "<<Jac(3,6)<<" "<<Jac(4,6)<<" "
2442  // <<Jac(5,6)<<" "<<Jac(6,6)<<std::endl;
2443  int i=5;
2444  if(GNorm.z() > 0.95)
2445  std::cout<<" Ecap1 N "<<GNorm<<std::endl;
2446  else if(GNorm.z() < -0.95)
2447  std::cout<<" Ecap2 N "<<GNorm<<std::endl;
2448  else
2449  std::cout<<" Barrel N "<<GNorm<<std::endl;
2450  std::cout<<" JacCLo(i,"<<i<<") = {"<<JacCorLoc(1,i)<<" "<<JacCorLoc(2,i)<<" "
2451  <<JacCorLoc(3,i)<<" "<<JacCorLoc(4,i)<<"}"<<std::endl;
2452  std::cout<<" rotLoc "<<rotLoc<<std::endl;
2453  std::cout<<" position "<<R0<<std::endl;
2454  std::cout<<" Pm,l "<<Pm.T()<<" "<<Pml(1)/Pml(3)<<" "<<Pml(2)/Pml(3)<<std::endl;
2455  std::cout<<" Pt,l "<<Pt.T()<<" "<<Ptl(1)/Ptl(3)<<" "<<Ptl(2)/Ptl(3)<<std::endl;
2456  std::cout<<" V "<<V.T()<<std::endl;
2457  std::cout<<" Cov \n"<<covLoc<<std::endl;
2458  std::cout<<" W*Cov "<<W*covLoc<<std::endl;
2459  //std::cout<<" JacCarToLoc = drldrc \n"<<
2460  //JacobianCartesianToLocal::JacobianCartesianToLocal
2461  //(refSurface, tsosTrack.localParameters()).jacobian()<<std::endl;
2462  std::cout<<" JacToLoc "<<JacToLoc<<std::endl;
2463 
2464  }
2465 
2466 #ifdef CHECK_OF_JACOBIAN_CARTESIAN_TO_LOCAL
2467  //---------------------- check of derivatives
2468  CLHEP::HepVector V0(4,0);
2469  V0(1) = Pml(1)/Pml(3) - Ptl(1)/Ptl(3);
2470  V0(2) = Pml(2)/Pml(3) - Ptl(2)/Ptl(3);
2471  V0(3) = Rml(1) - Rtl(1);
2472  V0(4) = Rml(2) - Rtl(2);
2473  int ii = 3; float delta = 0.01;
2474  CLHEP::HepVector V1(4,0);
2475  if(ii <= 3) {Rm(ii) += delta; Rml = rotLoc * (Rm - R0);}
2476  else {Pm(ii-3) += delta; Pml = rotLoc * Pm;}
2477  //if(ii <= 3) {Rt(ii) += delta; Rtl = rotLoc * (Rt - R0);}
2478  //else {Pt(ii-3) += delta; Ptl = rotLoc * Pt;}
2479  V1(1) = Pml(1)/Pml(3) - Ptl(1)/Ptl(3);
2480  V1(2) = Pml(2)/Pml(3) - Ptl(2)/Ptl(3);
2481  V1(3) = Rml(1) - Rtl(1);
2482  V1(4) = Rml(2) - Rtl(2);
2483 
2484  if(GNorm.z() > 0.95)
2485  std::cout<<" Ecap1 N "<<GNorm<<std::endl;
2486  else if(GNorm.z() < -0.95)
2487  std::cout<<" Ecap2 N "<<GNorm<<std::endl;
2488  else
2489  std::cout<<" Barrel N "<<GNorm<<std::endl;
2490  std::cout<<" dldc Num (i,"<<ii<<") "<<(V1(1)-V0(1))/delta<<" "<<(V1(2)-V0(2))/delta<<" "
2491  <<(V1(3)-V0(3))/delta<<" "<<(V1(4)-V0(4))/delta<<std::endl;
2492  std::cout<<" dldc Ana (i,"<<ii<<") "<<JacToLoc(1,ii)<<" "<<JacToLoc(2,ii)<<" "
2493  <<JacToLoc(3,ii)<<" "<<JacToLoc(4,ii)<<std::endl;
2494  //float dtxdpx = (rotLoc(1,1) - rotLoc(3,1)*Pml(1)/Pml(3))/Pml(3);
2495  //float dtydpx = (rotLoc(2,1) - rotLoc(3,2)*Pml(2)/Pml(3))/Pml(3);
2496  //std::cout<<" dtx/dpx dty/ "<<dtxdpx<<" "<<dtydpx<<std::endl;
2497 #endif
2498 
2499  return;
2500 } // 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:141
i
int i
Definition: DBlmapReader.cc:9
info
static const TGPicture * info(bool iBackgroundIsBlack)
Definition: FWCollectionSummaryWidget.cc:170
PV3DBase::y
T y() const
Definition: PV3DBase.h:63
cuy.ii
int ii
Definition: cuy.py:588
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
PV3DBase::z
T z() const
Definition: PV3DBase.h:64
j
int j
Definition: DBlmapReader.cc:9
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
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 2505 of file GlobalTrackerMuonAlignment.cc.

References funct::A, a, alignmentValidation::c1, funct::cos(), gather_cfg::cout, alignCSCRings::s, indexGen::s2, funct::sin(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

2507 {
2508  CLHEP::HepVector d(3,0), a(3,0), Norm(3), Rm0(3), Pm0(3), Rm1(3), Pm1(3),
2509  Rmi(3), Pmi(3) ;
2510 
2511  d(1)=0.0; d(2)=0.0; d(3)=0.0; a(1)=0.0000; a(2)=0.0000; a(3)=0.0000; // zero
2512  //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
2513  //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
2514  //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
2515  //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
2516  //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
2517  //d(1)=10.; d(2)=20.; d(3)=30.; a(1)=0.0100; a(2)=0.0200; a(3)=0.0300; // huge
2518  //d(1)=10.; d(2)=20.; d(3)=30.; a(1)=0.2000; a(2)=0.2500; a(3)=0.3000; // huge angles
2519  //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
2520  //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
2521  //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
2522 
2523  MuGlShift = d; MuGlAngle = a;
2524 
2525  if((d(1) == 0.) & (d(2) == 0.) & (d(3) == 0.) &
2526  (a(1) == 0.) & (a(2) == 0.) & (a(3) == 0.)){
2527  Rm = GRm;
2528  Pm = GPm;
2529  if(debug_)
2530  std::cout<<" ...... without misalignment "<<std::endl;
2531  return;
2532  }
2533 
2534  Rm0(1) = GRm.x(); Rm0(2) = GRm.y(); Rm0(3) = GRm.z();
2535  Pm0(1) = GPm.x(); Pm0(2) = GPm.y(); Pm0(3) = GPm.z();
2536  Norm(1) = Nl.x(); Norm(2) = Nl.y(); Norm(3) = Nl.z();
2537 
2538  CLHEP::HepMatrix T(3,3,1);
2539 
2540  //T(1,2) = a(3); T(1,3) = -a(2);
2541  //T(2,1) = -a(3); T(2,3) = a(1);
2542  //T(3,1) = a(2); T(3,2) = -a(1);
2543 
2544  double s1 = std::sin(a(1)), c1 = std::cos(a(1));
2545  double s2 = std::sin(a(2)), c2 = std::cos(a(2));
2546  double s3 = std::sin(a(3)), c3 = std::cos(a(3));
2547  T(1,1) = c2 * c3;
2548  T(1,2) = c1 * s3 + s1 * s2 * c3;
2549  T(1,3) = s1 * s3 - c1 * s2 * c3;
2550  T(2,1) = -c2 * s3;
2551  T(2,2) = c1 * c3 - s1 * s2 * s3;
2552  T(2,3) = s1 * c3 + c1 * s2 * s3;
2553  T(3,1) = s2;
2554  T(3,2) = -s1 * c2;
2555  T(3,3) = c1 * c2;
2556 
2557  //int terr;
2558  //CLHEP::HepMatrix Ti = T.inverse(terr);
2559  CLHEP::HepMatrix Ti = T.T();
2560  CLHEP::HepVector di = -Ti * d;
2561 
2562  CLHEP::HepVector ex0(3,0), ey0(3,0), ez0(3,0), ex(3,0), ey(3,0), ez(3,0);
2563  ex0(1) = 1.; ey0(2) = 1.; ez0(3) = 1.;
2564  ex = Ti*ex0; ey = Ti*ey0; ez = Ti*ez0;
2565 
2566  CLHEP::HepVector TiN = Ti * Norm;
2567  //double si = CLHEP::dot((Ti*Rm0 - Rm0 + di), TiN) / CLHEP::dot(Pm0, TiN);
2568  //Rm1(1) = CLHEP::dot(ex, Rm0 + si*Pm0 - di);
2569  //Rm1(2) = CLHEP::dot(ey, Rm0 + si*Pm0 - di);
2570  //Rm1(3) = CLHEP::dot(ez, Rm0 + si*Pm0 - di);
2571  double si = CLHEP_dot((Ti*Rm0 - Rm0 + di), TiN) / CLHEP_dot(Pm0, TiN);
2572  Rm1(1) = CLHEP_dot(ex, Rm0 + si*Pm0 - di);
2573  Rm1(2) = CLHEP_dot(ey, Rm0 + si*Pm0 - di);
2574  Rm1(3) = CLHEP_dot(ez, Rm0 + si*Pm0 - di);
2575  Pm1 = T * Pm0;
2576 
2577  Rm = GlobalVector(Rm1(1), Rm1(2), Rm1(3));
2578  //Pm = GlobalVector(CLHEP::dot(Pm0,ex), CLHEP::dot(Pm0, ey), CLHEP::dot(Pm0,ez));// =T*Pm0
2579  Pm = GlobalVector(CLHEP_dot(Pm0,ex), CLHEP_dot(Pm0, ey), CLHEP_dot(Pm0,ez));// =T*Pm0
2580 
2581  if(debug_){ // ------------- debug tranformation
2582 
2583  std::cout<<" ----- Pm "<<Pm<<std::endl;
2584  std::cout<<" T*Pm0 "<<Pm1.T()<<std::endl;
2585  //std::cout<<" Ti*T*Pm0 "<<(Ti*(T*Pm0)).T()<<std::endl;
2586 
2587  //CLHEP::HepVector Rml = (Rm0 + si*Pm0 - di) + di;
2588  CLHEP::HepVector Rml = Rm1(1)*ex + Rm1(2)*ey + Rm1(3)*ez + di;
2589 
2590  CLHEP::HepVector Rt0(3);
2591  Rt0(1)=Rt.x(); Rt0(2)=Rt.y(); Rt0(3)=Rt.z();
2592 
2593  //double sl = CLHEP::dot(T*(Rt0 - Rml), T*Norm) / CLHEP::dot(Ti * Pm1, Norm);
2594  double sl = CLHEP_dot(T*(Rt0 - Rml), T*Norm) / CLHEP_dot(Ti * Pm1, Norm);
2595  CLHEP::HepVector Rl = Rml + sl*(Ti*Pm1);
2596 
2597  //double A = CLHEP::dot(Ti*Pm1, Norm);
2598  //double B = CLHEP::dot(Rt0, Norm) - CLHEP::dot(Ti*Rm1, Norm)
2599  //+ CLHEP::dot(Ti*d, Norm);
2600  double A = CLHEP_dot(Ti*Pm1, Norm);
2601  double B = CLHEP_dot(Rt0, Norm) - CLHEP_dot(Ti*Rm1, Norm)
2602  + CLHEP_dot(Ti*d, Norm);
2603  double s = B/A;
2604  CLHEP::HepVector Dr = Ti*Rm1 - Ti*d + s*(Ti*Pm1) - Rm0;
2605  CLHEP::HepVector Dp = Ti*Pm1 - Pm0;
2606 
2607  CLHEP::HepVector TiR = Ti * Rm0 + di;
2608  CLHEP::HepVector Ri = T * TiR + d;
2609 
2610  std::cout<<" --TTi-0 "<<(Ri-Rm0).T()<<std::endl;
2611  std::cout<<" -- Dr "<<Dr.T()<<endl;
2612  std::cout<<" -- Dp "<<Dp.T()<<endl;
2613  //std::cout<<" ex "<<ex.T()<<endl;
2614  //std::cout<<" ey "<<ey.T()<<endl;
2615  //std::cout<<" ez "<<ez.T()<<endl;
2616  //std::cout<<" ---- T ---- "<<T<<std::endl;
2617  //std::cout<<" ---- Ti ---- "<<Ti<<std::endl;
2618  //std::cout<<" ---- d ---- "<<d.T()<<std::endl;
2619  //std::cout<<" ---- di ---- "<<di.T()<<std::endl;
2620  std::cout<<" -- si sl s "<<si<<" "<<sl<<" "<<s<<std::endl;
2621  //std::cout<<" -- si sl "<<si<<" "<<sl<<endl;
2622  //std::cout<<" -- si s "<<si<<" "<<s<<endl;
2623  std::cout<<" -- Rml-(Rm0+si*Pm0) "<<(Rml - Rm0 - si*Pm0).T()<<std::endl;
2624  //std::cout<<" -- Rm0 "<<Rm0.T()<<std::endl;
2625  //std::cout<<" -- Rm1 "<<Rm1.T()<<std::endl;
2626  //std::cout<<" -- Rmi "<<Rmi.T()<<std::endl;
2627  //std::cout<<" --siPm "<<(si*Pm0).T()<<std::endl;
2628  //std::cout<<" --s2Pm "<<(s2*(T * Pm1)).T()<<std::endl;
2629  //std::cout<<" --R1-0 "<<(Rm1-Rm0).T()<<std::endl;
2630  //std::cout<<" --Ri-0 "<<(Rmi-Rm0).T()<<std::endl;
2631  std::cout<<" --Rl-0 "<<(Rl-Rm0).T()<<std::endl;
2632  //std::cout<<" --Pi-0 "<<(Pmi-Pm0).T()<<std::endl;
2633  } // -------- end of debug
2634 
2635  return;
2636 
2637 } // end of misalignMuon
GlobalTrackerMuonAlignment::CLHEP_dot
double CLHEP_dot(const CLHEP::HepVector &a, const CLHEP::HepVector &b)
Definition: GlobalTrackerMuonAlignment.cc:141
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
PV3DBase::y
T y() const
Definition: PV3DBase.h:63
indexGen.s2
tuple s2
Definition: indexGen.py:106
double_binary
Definition: DDStreamer.cc:227
PV3DBase::z
T z() const
Definition: PV3DBase.h:64
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
a
double a
Definition: hdecay.h:121
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
T
long double T
Definition: Basic3DVectorLD.h:57
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 2641 of file GlobalTrackerMuonAlignment.cc.

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

Referenced by analyzeTrackTrack(), and analyzeTrackTrajectory().

2646 {
2647  CLHEP::HepVector d(3,0), a(3,0), Norm(3), Rm0(3), Pm0(3), Rm1(3), Pm1(3),
2648  Rmi(3), Pmi(3);
2649 
2650  d(1)=0.0; d(2)=0.0; d(3)=0.0; a(1)=0.0000; a(2)=0.0000; a(3)=0.0000; // zero
2651  //d(1)=0.0000001; d(2)=0.0; d(3)=0.0; a(1)=0.0000; a(2)=0.0000; a(3)=0.0000; // nonzero
2652  //d(1)=0.0; d(2)=0.0; d(3)=0.0; a(1)=0.0020; a(2)=0.0000; a(3)=0.0000; // a1=.002
2653  //d(1)=0.0; d(2)=0.0; d(3)=0.0; a(1)=0.0000; a(2)=0.0020; a(3)=0.0000; // a2=.002
2654  //d(1)=0.0; d(2)=0.0; d(3)=0.0; a(1)=0.0000; a(2)=0.0000; a(3)=0.0020; // a3=.002
2655  //d(1)=0.0; d(2)=0.0; d(3)=0.0; a(1)=0.0000; a(2)=0.0000; a(3)=0.0200; // a3=.020
2656  //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
2657  //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
2658  //d(1)=10.; d(2)=20.; d(3)=30.; a(1)=0.0100; a(2)=0.0200; a(3)=0.0300; // huge
2659  //d(1)=10.; d(2)=20.; d(3)=30.; a(1)=0.2000; a(2)=0.2500; a(3)=0.3000; // huge angles
2660  //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
2661  //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
2662  //d(1)=0.1; d(2)=0.2; d(3)=0.3; a(1)=0.0003; a(2)=0.0004; a(3)=0.0005; // 0123 0345
2663 
2664  MuGlShift = d; MuGlAngle = a;
2665 
2666  if((d(1) == 0.) & (d(2) == 0.) & (d(3) == 0.) &
2667  (a(1) == 0.) & (a(2) == 0.) & (a(3) == 0.)){
2668  Rm = GRm;
2669  Pm = GPm;
2670  AlgebraicVector4 Vm0;
2671  Vm = AlgebraicVector4(tsosMuon.localParameters().vector()(1),
2672  tsosMuon.localParameters().vector()(2),
2673  tsosMuon.localParameters().vector()(3),
2674  tsosMuon.localParameters().vector()(4));
2675  if(debug_)
2676  std::cout<<" ...... without misalignment "<<std::endl;
2677  return;
2678  }
2679 
2680  Rm0(1) = GRm.x(); Rm0(2) = GRm.y(); Rm0(3) = GRm.z();
2681  Pm0(1) = GPm.x(); Pm0(2) = GPm.y(); Pm0(3) = GPm.z();
2682  Norm(1) = Nl.x(); Norm(2) = Nl.y(); Norm(3) = Nl.z();
2683 
2684  CLHEP::HepMatrix T(3,3,1);
2685 
2686  //T(1,2) = a(3); T(1,3) = -a(2);
2687  //T(2,1) = -a(3); T(2,3) = a(1);
2688  //T(3,1) = a(2); T(3,2) = -a(1);
2689 
2690  double s1 = std::sin(a(1)), c1 = std::cos(a(1));
2691  double s2 = std::sin(a(2)), c2 = std::cos(a(2));
2692  double s3 = std::sin(a(3)), c3 = std::cos(a(3));
2693  T(1,1) = c2 * c3;
2694  T(1,2) = c1 * s3 + s1 * s2 * c3;
2695  T(1,3) = s1 * s3 - c1 * s2 * c3;
2696  T(2,1) = -c2 * s3;
2697  T(2,2) = c1 * c3 - s1 * s2 * s3;
2698  T(2,3) = s1 * c3 + c1 * s2 * s3;
2699  T(3,1) = s2;
2700  T(3,2) = -s1 * c2;
2701  T(3,3) = c1 * c2;
2702 
2703  // Ti, di what we have to apply for misalignment
2704  //int terr;
2705  //CLHEP::HepMatrix Ti = T.inverse(terr);
2706  CLHEP::HepMatrix Ti = T.T();
2707  CLHEP::HepVector di = -Ti * d;
2708 
2709  CLHEP::HepVector ex0(3,0), ey0(3,0), ez0(3,0), ex(3,0), ey(3,0), ez(3,0);
2710  ex0(1) = 1.; ey0(2) = 1.; ez0(3) = 1.;
2711  ex = Ti*ex0; ey = Ti*ey0; ez = Ti*ez0;
2712 
2713  CLHEP::HepVector TiN = Ti * Norm;
2714  //double si = CLHEP::dot((Ti*Rm0 - Rm0 + di), TiN) / CLHEP::dot(Pm0, TiN);
2715  //Rm1(1) = CLHEP::dot(ex, Rm0 + si*Pm0 - di);
2716  //Rm1(2) = CLHEP::dot(ey, Rm0 + si*Pm0 - di);
2717  //Rm1(3) = CLHEP::dot(ez, Rm0 + si*Pm0 - di);
2718  double si = CLHEP_dot((Ti*Rm0 - Rm0 + di), TiN) / CLHEP_dot(Pm0, TiN);
2719  Rm1(1) = CLHEP_dot(ex, Rm0 + si*Pm0 - di);
2720  Rm1(2) = CLHEP_dot(ey, Rm0 + si*Pm0 - di);
2721  Rm1(3) = CLHEP_dot(ez, Rm0 + si*Pm0 - di);
2722  Pm1 = T * Pm0;
2723 
2724  // Global muon parameters after misalignment
2725  Rm = GlobalVector(Rm1(1), Rm1(2), Rm1(3));
2726  //Pm = GlobalVector(CLHEP::dot(Pm0,ex), CLHEP::dot(Pm0, ey), CLHEP::dot(Pm0,ez));// =T*Pm0
2727  Pm = GlobalVector(CLHEP_dot(Pm0,ex), CLHEP_dot(Pm0, ey), CLHEP_dot(Pm0,ez));// =T*Pm0
2728 
2729  // Local muon parameters after misalignment
2730  const Surface& refSurface = tsosMuon.surface();
2731  CLHEP::HepMatrix Tl(3,3,0);
2732 
2733  Tl(1,1) = refSurface.rotation().xx();
2734  Tl(1,2) = refSurface.rotation().xy();
2735  Tl(1,3) = refSurface.rotation().xz();
2736 
2737  Tl(2,1) = refSurface.rotation().yx();
2738  Tl(2,2) = refSurface.rotation().yy();
2739  Tl(2,3) = refSurface.rotation().yz();
2740 
2741  Tl(3,1) = refSurface.rotation().zx();
2742  Tl(3,2) = refSurface.rotation().zy();
2743  Tl(3,3) = refSurface.rotation().zz();
2744 
2745  CLHEP::HepVector R0(3,0), newR0(3,0), newRl(3,0), newPl(3,0);
2746  R0(1) = refSurface.position().x();
2747  R0(2) = refSurface.position().y();
2748  R0(3) = refSurface.position().z();
2749 
2750  newRl = Tl * (Rm1 - R0);
2751  newPl = Tl * Pm1;
2752 
2753  Vm(0) = newPl(1)/newPl(3);
2754  Vm(1) = newPl(2)/newPl(3);
2755  Vm(2) = newRl(1);
2756  Vm(3) = newRl(2);
2757 
2758 #ifdef CHECK_DERIVATIVES_LOCAL_VS_ANGLES
2759  float del = 0.0001;
2760  //int ii = 1; T(3,2) +=del; T(2,3) -=del;
2761  int ii = 2; T(3,1) -=del; T(1,3) +=del;
2762  //int ii = 3; T(1,2) -=del; T(2,1) +=del;
2763  AlgebraicVector4 Vm0 = Vm;
2764  CLHEP::HepVector Rm10 = Rm1, Pm10 = Pm1;;
2765  CLHEP::HepMatrix newTl = Tl*T;
2766  Ti = T.T();
2767  di = -Ti * d;
2768  ex = Ti*ex0; ey = Ti*ey0; ez = Ti*ez0;
2769  TiN = Ti * Norm;
2770  si = CLHEP_dot((Ti*Rm0 - Rm0 + di), TiN) / CLHEP_dot(Pm0, TiN);
2771  Rm1(1) = CLHEP_dot(ex, Rm0 + si*Pm0 - di);
2772  Rm1(2) = CLHEP_dot(ey, Rm0 + si*Pm0 - di);
2773  Rm1(3) = CLHEP_dot(ez, Rm0 + si*Pm0 - di);
2774  Pm1 = T * Pm0;
2775 
2776  newRl = Tl * (Rm1 - R0);
2777  newPl = Tl * Pm1;
2778 
2779  Vm(0) = newPl(1)/newPl(3);
2780  Vm(1) = newPl(2)/newPl(3);
2781  Vm(2) = newRl(1);
2782  Vm(3) = newRl(2);
2783  std::cout<<" ========= dVm/da"<<ii<<" "<<(Vm-Vm0)*(1./del)<<std::endl;
2784 #endif
2785 
2786  if(debug_){
2787  //std::cout<<" dRm/da3 "<<((Rm1-Rm10)*(1./del)).T()<<" "<<((Pm1-Pm10)*(1./del)).T()<<std::endl;
2788  //std::cout<<" Norm "<<Norm.T()<<std::endl;
2789  std::cout<<" ex "<<(Tl.T()*ex0).T()<<std::endl;
2790  std::cout<<" ey "<<(Tl.T()*ey0).T()<<std::endl;
2791  std::cout<<" ez "<<(Tl.T()*ez0).T()<<std::endl;
2792  //std::cpot<<"
2793 
2794  std::cout<<" pxyz/p gl 0 "<<(Pm0/Pm0.norm()).T()<<std::endl;
2795  std::cout<<" pxyz/p loc0 "<<(Tl*Pm0/(Tl*Pm0).norm()).T()<<std::endl;
2796  std::cout<<" pxyz/p glob "<<(Pm1/Pm1.norm()).T()<<std::endl;
2797  std::cout<<" pxyz/p loc "<<(newPl/newPl.norm()).T()<<std::endl;
2798 
2799  //std::cout<<" Rot "<<refSurface.rotation()<<endl;
2800  //std::cout<<" Tl "<<Tl<<endl;
2801  std::cout<<" ---- phi g0 g1 l0 l1 "
2802  <<atan2(Pm0(2),Pm0(1))<<" "<<atan2(Pm1(2),Pm1(1))<<" "
2803  <<atan2((Tl*Pm0)(2),(Tl*Pm0)(1))<<" "
2804  <<atan2(newPl(2),newPl(1))<<std::endl;
2805  std::cout<<" ---- angl Gl Loc "<<atan2(Pm1(2),Pm1(1))-atan2(Pm0(2),Pm0(1))<<" "
2806  <<atan2(newPl(2),newPl(1))-atan2((Tl*Pm0)(2),(Tl*Pm0)(1))<<" "
2807  <<atan2(newPl(3),newPl(2))-atan2((Tl*Pm0)(3),(Tl*Pm0)(2))<<" "
2808  <<atan2(newPl(1),newPl(3))-atan2((Tl*Pm0)(1),(Tl*Pm0)(3))<<" "<<std::endl;
2809  }
2810 
2811  if(debug_){
2812  CLHEP::HepMatrix newTl(3,3,0);
2813  CLHEP::HepVector R0(3,0), newRl(3,0), newPl(3,0);
2814  newTl = Tl * Ti.T();
2815  newR0 = Ti * R0 + di;
2816 
2817  std::cout<<" N "<<Norm.T()<<std::endl;
2818  std::cout<<" dtxl yl "<<Vm(0)-tsosMuon.localParameters().vector()(1)<<" "
2819  <<Vm(1)-tsosMuon.localParameters().vector()(2)<<std::endl;
2820  std::cout<<" dXl dYl "<<Vm(2)-tsosMuon.localParameters().vector()(3)<<" "
2821  <<Vm(3)-tsosMuon.localParameters().vector()(4)<<std::endl;
2822  std::cout<<" localM "<<tsosMuon.localParameters().vector()<<std::endl;
2823  std::cout<<" Vm "<<Vm<<std::endl;
2824 
2825 
2826  CLHEP::HepVector Rvc(3,0), Pvc(3,0), Rvg(3,0), Pvg(3,0);
2827  Rvc(1) = Vm(2); Rvc(2) = Vm(3);
2828  Pvc(3) = tsosMuon.localParameters().pzSign() * Pm0.norm() /
2829  sqrt(1 + Vm(0)*Vm(0) + Vm(1)*Vm(1));
2830  Pvc(1) = Pvc(3) * Vm(0);
2831  Pvc(2) = Pvc(3) * Vm(1);
2832 
2833  Rvg = newTl.T() * Rvc + newR0;
2834  Pvg = newTl.T() * Pvc;
2835 
2836  std::cout<<" newPl "<<newPl.T()<<std::endl;
2837  std::cout<<" Pvc "<<Pvc.T()<<std::endl;
2838  std::cout<<" Rvg "<<Rvg.T()<<std::endl;
2839  std::cout<<" Rm1 "<<Rm1.T()<<std::endl;
2840  std::cout<<" Pvg "<<Pvg.T()<<std::endl;
2841  std::cout<<" Pm1 "<<Pm1.T()<<std::endl;
2842  }
2843 
2844  if(debug_){ // ---------- how to transform cartesian from shifted to correct
2845 
2846  std::cout<<" ----- Pm "<<Pm<<std::endl;
2847  std::cout<<" T*Pm0 "<<Pm1.T()<<std::endl;
2848  //std::cout<<" Ti*T*Pm0 "<<(Ti*(T*Pm0)).T()<<std::endl;
2849 
2850  //CLHEP::HepVector Rml = (Rm0 + si*Pm0 - di) + di;
2851  CLHEP::HepVector Rml = Rm1(1)*ex + Rm1(2)*ey + Rm1(3)*ez + di;
2852 
2853  CLHEP::HepVector Rt0(3);
2854  Rt0(1)=Rt.x(); Rt0(2)=Rt.y(); Rt0(3)=Rt.z();
2855 
2856  //double sl = CLHEP::dot(T*(Rt0 - Rml), T*Norm) / CLHEP::dot(Ti * Pm1, Norm);
2857  double sl = CLHEP_dot(T*(Rt0 - Rml), T*Norm) / CLHEP_dot(Ti * Pm1, Norm);
2858  CLHEP::HepVector Rl = Rml + sl*(Ti*Pm1);
2859 
2860  //double A = CLHEP::dot(Ti*Pm1, Norm);
2861  //double B = CLHEP::dot(Rt0, Norm) - CLHEP::dot(Ti*Rm1, Norm)
2862  //+ CLHEP::dot(Ti*d, Norm);
2863  double A = CLHEP_dot(Ti*Pm1, Norm);
2864  double B = CLHEP_dot(Rt0, Norm) - CLHEP_dot(Ti*Rm1, Norm)
2865  + CLHEP_dot(Ti*d, Norm);
2866  double s = B/A;
2867  CLHEP::HepVector Dr = Ti*Rm1 - Ti*d + s*(Ti*Pm1) - Rm0;
2868  CLHEP::HepVector Dp = Ti*Pm1 - Pm0;
2869 
2870  CLHEP::HepVector TiR = Ti * Rm0 + di;
2871  CLHEP::HepVector Ri = T * TiR + d;
2872 
2873  std::cout<<" --TTi-0 "<<(Ri-Rm0).T()<<std::endl;
2874  std::cout<<" -- Dr "<<Dr.T()<<endl;
2875  std::cout<<" -- Dp "<<Dp.T()<<endl;
2876  //std::cout<<" ex "<<ex.T()<<endl;
2877  //std::cout<<" ey "<<ey.T()<<endl;
2878  //std::cout<<" ez "<<ez.T()<<endl;
2879  //std::cout<<" ---- T ---- "<<T<<std::endl;
2880  //std::cout<<" ---- Ti ---- "<<Ti<<std::endl;
2881  //std::cout<<" ---- d ---- "<<d.T()<<std::endl;
2882  //std::cout<<" ---- di ---- "<<di.T()<<std::endl;
2883  std::cout<<" -- si sl s "<<si<<" "<<sl<<" "<<s<<std::endl;
2884  //std::cout<<" -- si sl "<<si<<" "<<sl<<endl;
2885  //std::cout<<" -- si s "<<si<<" "<<s<<endl;
2886  std::cout<<" -- Rml-(Rm0+si*Pm0) "<<(Rml - Rm0 - si*Pm0).T()<<std::endl;
2887  //std::cout<<" -- Rm0 "<<Rm0.T()<<std::endl;
2888  //std::cout<<" -- Rm1 "<<Rm1.T()<<std::endl;
2889  //std::cout<<" -- Rmi "<<Rmi.T()<<std::endl;
2890  //std::cout<<" --siPm "<<(si*Pm0).T()<<std::endl;
2891  //std::cout<<" --s2Pm "<<(s2*(T * Pm1)).T()<<std::endl;
2892  //std::cout<<" --R1-0 "<<(Rm1-Rm0).T()<<std::endl;
2893  //std::cout<<" --Ri-0 "<<(Rmi-Rm0).T()<<std::endl;
2894  std::cout<<" --Rl-0 "<<(Rl-Rm0).T()<<std::endl;
2895  //std::cout<<" --Pi-0 "<<(Pmi-Pm0).T()<<std::endl;
2896  } // -------- end of debug
2897 
2898  return;
2899 
2900 } // end misalignMuonL
TkRotation::xx
T xx() const
Definition: extTkRotation.h:251
GlobalTrackerMuonAlignment::CLHEP_dot
double CLHEP_dot(const CLHEP::HepVector &a, const CLHEP::HepVector &b)
Definition: GlobalTrackerMuonAlignment.cc:141
LocalTrajectoryParameters::pzSign
double pzSign() const
Sign of the z-component of the momentum in the local frame.
Definition: LocalTrajectoryParameters.h:146
TrajectoryStateOnSurface::localParameters
const LocalTrajectoryParameters & localParameters() const
Definition: TrajectoryStateOnSurface.h:113
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
PV3DBase::y
T y() const
Definition: PV3DBase.h:63
TkRotation::yx
T yx() const
Definition: extTkRotation.h:254
cuy.ii
int ii
Definition: cuy.py:588
indexGen.s2
tuple s2
Definition: indexGen.py:106
double_binary
Definition: DDStreamer.cc:227
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:120
TkRotation::zz
T zz() const
Definition: extTkRotation.h:259
TrajectoryStateOnSurface::surface
const SurfaceType & surface() const
Definition: TrajectoryStateOnSurface.h:128
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
PV3DBase::z
T z() const
Definition: PV3DBase.h:64
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
TkRotation::zy
T zy() const
Definition: extTkRotation.h:258
TkRotation::yy
T yy() const
Definition: extTkRotation.h:255
a
double a
Definition: hdecay.h:121
TkRotation::xz
T xz() const
Definition: extTkRotation.h:253
GloballyPositioned::rotation
const RotationType & rotation() const
Definition: GloballyPositioned.h:44
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
AlgebraicVector4
ROOT::Math::SVector< double, 4 > AlgebraicVector4
Definition: AlgebraicROOTObjects.h:21
T
long double T
Definition: Basic3DVectorLD.h:57
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::muonFitter ( reco::TrackRef  alongTr,
reco::TransientTrack &  alongTTr,
PropagationDirection  direction,
TrajectoryStateOnSurface &  trackFittedTSOS 
)

Definition at line 3019 of file GlobalTrackerMuonAlignment.cc.

References alongMomentum, edm::OwnVector< T, P >::clear(), gather_cfg::cout, i, info(), reco::TransientTrack::innermostMeasurementState(), TrajectoryStateOnSurface::localError(), TrajectoryStateOnSurface::localParameters(), LocalTrajectoryError::matrix(), reco::TransientTrack::outermostMeasurementState(), trajectoryStateTransform::persistentState(), edm::OwnVector< T, P >::push_back(), DetId::rawId(), reco::TransientTrack::recHitsBegin(), reco::TransientTrack::recHitsEnd(), edm::OwnVector< T, P >::size(), and TrajectoryStateOnSurface::surface().

Referenced by analyzeTrackTrajectory().

3022 {
3023  bool info = false;
3024  bool smooth = false;
3025 
3026  if(info) std::cout<<" ......... start of muonFitter ........"<<std::endl;
3027  if(false && info){
3028  this->debugTrackHit(" Muon track hits ", alongTr);
3029  this->debugTrackHit(" Muon TransientTrack hits ", alongTTr);
3030  }
3031 
3032  edm::OwnVector<TrackingRecHit> recHit;
3033  DetId trackDetId = DetId(alongTr->innerDetId());
3034  for (trackingRecHit_iterator i=alongTr->recHitsBegin();i!=alongTr->recHitsEnd(); i++) {
3035  recHit.push_back((*i)->clone());
3036  }
3037  if(direction != alongMomentum)
3038  {
3039  edm::OwnVector<TrackingRecHit> recHitAlong = recHit;
3040  recHit.clear();
3041  for(unsigned int ihit = recHitAlong.size()-1; ihit+1>0; ihit--){
3042  recHit.push_back(recHitAlong[ihit]);
3043  }
3044  recHitAlong.clear();
3045  }
3046  if(info)
3047  std::cout<<" ... Own recHit.size() DetId==Muon "<<recHit.size()<<" "<<trackDetId.rawId()<<std::endl;
3048 
3049  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu =
3050  MuRHBuilder->build(alongTTr.recHitsBegin(), alongTTr.recHitsEnd());
3051  if(info)
3052  std::cout<<" ...along.... recHitMu.size() "<<recHitMu.size()<<std::endl;
3053  if(direction != alongMomentum){
3054  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMuAlong = recHitMu;
3055  recHitMu.clear();
3056  for(unsigned int ihit = recHitMuAlong.size()-1; ihit+1>0; ihit--){
3057  recHitMu.push_back(recHitMuAlong[ihit]);
3058  }
3059  recHitMuAlong.clear();
3060  }
3061  if(info)
3062  std::cout<<" ...opposite... recHitMu.size() "<<recHitMu.size()<<std::endl;
3063 
3064  TrajectoryStateOnSurface firstTSOS;
3065  if(direction == alongMomentum) firstTSOS = alongTTr.innermostMeasurementState();
3066  else firstTSOS = alongTTr.outermostMeasurementState();
3067 
3068  AlgebraicSymMatrix55 CovLoc;
3069  CovLoc(0,0) = 1. * firstTSOS.localError().matrix()(0,0);
3070  CovLoc(1,1) = 10. * firstTSOS.localError().matrix()(1,1);
3071  CovLoc(2,2) = 10. * firstTSOS.localError().matrix()(2,2);
3072  CovLoc(3,3) = 100. * firstTSOS.localError().matrix()(3,3);
3073  CovLoc(4,4) = 100. * firstTSOS.localError().matrix()(4,4);
3074  TrajectoryStateOnSurface initialTSOS(firstTSOS.localParameters(), LocalTrajectoryError(CovLoc),
3075  firstTSOS.surface(), &*magneticField_);
3076 
3077  PTrajectoryStateOnDet PTraj =
3078  trajectoryStateTransform::persistentState(initialTSOS, trackDetId.rawId());
3079  //const TrajectorySeed seedT(*PTraj, recHit, alongMomentum);
3080  const TrajectorySeed seedT(PTraj, recHit, direction);
3081 
3082  std::vector<Trajectory> trajVec;
3083  if(direction == alongMomentum) trajVec = theFitter->fit(seedT, recHitMu, initialTSOS);
3084  else trajVec = theFitterOp->fit(seedT, recHitMu, initialTSOS);
3085 
3086  if(info){
3087  this->debugTrajectorySOSv(" innermostTSOS of TransientTrack",
3088  alongTTr.innermostMeasurementState());
3089  this->debugTrajectorySOSv(" outermostTSOS of TransientTrack",
3090  alongTTr.outermostMeasurementState());
3091  this->debugTrajectorySOS(" initialTSOS for theFitter ", initialTSOS);
3092  std::cout<<" . . . . . trajVec.size() "<<trajVec.size()<<std::endl;
3093  if(trajVec.size()) this->debugTrajectory(" theFitter trajectory ", trajVec[0]);
3094  }
3095 
3096  if(!smooth){
3097  if(trajVec.size()) trackFittedTSOS = trajVec[0].lastMeasurement().updatedState();}
3098  else{
3099  std::vector<Trajectory> trajSVec;
3100  if (trajVec.size()){
3101  if(direction == alongMomentum) trajSVec = theSmoother->trajectories(*(trajVec.begin()));
3102  else trajSVec = theSmootherOp->trajectories(*(trajVec.begin()));
3103  }
3104  if(info) std::cout<<" . . . . trajSVec.size() "<<trajSVec.size()<<std::endl;
3105  if(trajSVec.size()) this->debugTrajectorySOSv("smoothed geom InnermostState",
3106  trajSVec[0].geometricalInnermostState());
3107  if(trajSVec.size()) trackFittedTSOS = trajSVec[0].geometricalInnermostState();
3108  }
3109 
3110 } // end of muonFitter
i
int i
Definition: DBlmapReader.cc:9
MuonTransientTrackingRecHitBuilder::build
RecHitPointer build(const TrackingRecHit *p, edm::ESHandle< GlobalTrackingGeometry > trackingGeometry) const
Call the MuonTransientTrackingRecHit::specificBuild.
Definition: MuonTransientTrackingRecHitBuilder.cc:22
GlobalTrackerMuonAlignment::debugTrajectorySOSv
void debugTrajectorySOSv(const std::string, TrajectoryStateOnSurface)
Definition: GlobalTrackerMuonAlignment.cc:3179
info
static const TGPicture * info(bool iBackgroundIsBlack)
Definition: FWCollectionSummaryWidget.cc:170
TrajectoryStateOnSurface::localParameters
const LocalTrajectoryParameters & localParameters() const
Definition: TrajectoryStateOnSurface.h:113
MuonTransientTrackingRecHitBuilder::ConstRecHitContainer
TransientTrackingRecHit::ConstRecHitContainer ConstRecHitContainer
Definition: MuonTransientTrackingRecHitBuilder.h:17
edm::OwnVector::size
size_type size() const
Definition: OwnVector.h:248
trajectoryStateTransform::persistentState
PTrajectoryStateOnDet persistentState(const TrajectoryStateOnSurface &ts, unsigned int detid)
Definition: TrajectoryStateTransform.cc:16
AlgebraicSymMatrix55
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
Definition: AlgebraicROOTObjects.h:31
reco::TransientTrack::innermostMeasurementState
TrajectoryStateOnSurface innermostMeasurementState() const
Definition: TransientTrack.h:81
DetId::rawId
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
edm::OwnVector::push_back
void push_back(D *&d)
Definition: OwnVector.h:274
TrajectoryStateOnSurface::surface
const SurfaceType & surface() const
Definition: TrajectoryStateOnSurface.h:128
GlobalTrackerMuonAlignment::debugTrackHit
void debugTrackHit(const std::string, reco::TrackRef)
Definition: GlobalTrackerMuonAlignment.cc:3131
edm::OwnVector::clear
void clear()
Definition: OwnVector.h:371
LocalTrajectoryError::matrix
const AlgebraicSymMatrix55 & matrix() const
Definition: LocalTrajectoryError.h:61
TrajectoryStateOnSurface::localError
const LocalTrajectoryError & localError() const
Definition: TrajectoryStateOnSurface.h:125
reco::TransientTrack::outermostMeasurementState
TrajectoryStateOnSurface outermostMeasurementState() const
Definition: TransientTrack.h:78
reco::TransientTrack::recHitsEnd
trackingRecHit_iterator recHitsEnd() const
last iterator to RecHits
Definition: TransientTrack.h:126
DetId
Definition: DetId.h:18
GlobalTrackerMuonAlignment::debugTrajectory
void debugTrajectory(const std::string, Trajectory &)
Definition: GlobalTrackerMuonAlignment.cc:3211
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
reco::TransientTrack::recHitsBegin
trackingRecHit_iterator recHitsBegin() const
first iterator to RecHits
Definition: TransientTrack.h:124
GlobalTrackerMuonAlignment::MuRHBuilder
MuonTransientTrackingRecHitBuilder * MuRHBuilder
Definition: GlobalTrackerMuonAlignment.cc:187
GlobalTrackerMuonAlignment::debugTrajectorySOS
void debugTrajectorySOS(const std::string, TrajectoryStateOnSurface &)
Definition: GlobalTrackerMuonAlignment.cc:3147
void GlobalTrackerMuonAlignment::trackFitter ( reco::TrackRef  alongTr,
reco::TransientTrack &  alongTTr,
PropagationDirection  direction,
TrajectoryStateOnSurface &  trackFittedTSOS 
)

Definition at line 2905 of file GlobalTrackerMuonAlignment.cc.

References alongMomentum, edm::OwnVector< T, P >::clear(), gather_cfg::cout, i, info(), reco::TransientTrack::innermostMeasurementState(), TrajectoryStateOnSurface::localError(), TrajectoryStateOnSurface::localParameters(), LocalTrajectoryError::matrix(), reco::TransientTrack::outermostMeasurementState(), trajectoryStateTransform::persistentState(), edm::OwnVector< T, P >::push_back(), DetId::rawId(), edm::OwnVector< T, P >::size(), TrajectoryStateOnSurface::surface(), and MixedLayerPairs_cfi::TTRHBuilder.

Referenced by analyzeTrackTrajectory().

2908 {
2909  bool info = false;
2910  bool smooth = false;
2911 
2912  if(info)
2913  std::cout<<" ......... start of trackFitter ......... "<<std::endl;
2914  if(false && info){
2915  this->debugTrackHit(" Tracker track hits ", alongTr);
2916  this->debugTrackHit(" Tracker TransientTrack hits ", alongTTr);
2917  }
2918 
2919  edm::OwnVector<TrackingRecHit> recHit;
2920  DetId trackDetId = DetId(alongTr->innerDetId());
2921  for (trackingRecHit_iterator i=alongTr->recHitsBegin();i!=alongTr->recHitsEnd(); i++) {
2922  recHit.push_back((*i)->clone());
2923  }
2924  if(direction != alongMomentum)
2925  {
2926  edm::OwnVector<TrackingRecHit> recHitAlong = recHit;
2927  recHit.clear();
2928  for(unsigned int ihit = recHitAlong.size()-1; ihit+1>0; ihit--){
2929  recHit.push_back(recHitAlong[ihit]);
2930  }
2931  recHitAlong.clear();
2932  }
2933  if(info)
2934  std::cout<<" ... Own recHit.size() "<<recHit.size()<<" "<<trackDetId.rawId()<<std::endl;
2935 
2936  //MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitTrack;
2937  TransientTrackingRecHit::RecHitContainer recHitTrack;
2938  for (trackingRecHit_iterator i=alongTr->recHitsBegin();i!=alongTr->recHitsEnd(); i++) {
2939  if((*i)->isValid()){
2940  recHitTrack.push_back(TTRHBuilder->build(&**i ));}
2941  else{
2942  recHitTrack.push_back(TTRHBuilder->build(&**i ));}
2943  }
2944 
2945  if(info)
2946  std::cout<<" ... recHitTrack.size() "<<recHit.size()<<" "<<trackDetId.rawId()
2947  <<" ======> recHitMu "<<std::endl;
2948 
2949  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu = recHitTrack;
2950  /*
2951  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu =
2952  MuRHBuilder->build(alongTTr.recHitsBegin(), alongTTr.recHitsEnd());
2953  */
2954  if(info)
2955  std::cout<<" ...along.... recHitMu.size() "<<recHitMu.size()<<std::endl;
2956  if(direction != alongMomentum){
2957  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMuAlong = recHitMu;
2958  recHitMu.clear();
2959  for(unsigned int ihit = recHitMuAlong.size()-1; ihit+1>0; ihit--){
2960  recHitMu.push_back(recHitMuAlong[ihit]);
2961  }
2962  recHitMuAlong.clear();
2963  }
2964  if(info)
2965  std::cout<<" ...opposite... recHitMu.size() "<<recHitMu.size()<<std::endl;
2966 
2967  TrajectoryStateOnSurface firstTSOS;
2968  if(direction == alongMomentum) firstTSOS = alongTTr.innermostMeasurementState();
2969  else firstTSOS = alongTTr.outermostMeasurementState();
2970 
2971  AlgebraicSymMatrix55 CovLoc;
2972  CovLoc(0,0) = 1. * firstTSOS.localError().matrix()(0,0);
2973  CovLoc(1,1) = 10. * firstTSOS.localError().matrix()(1,1);
2974  CovLoc(2,2) = 10. * firstTSOS.localError().matrix()(2,2);
2975  CovLoc(3,3) = 100. * firstTSOS.localError().matrix()(3,3);
2976  CovLoc(4,4) = 100. * firstTSOS.localError().matrix()(4,4);
2977  TrajectoryStateOnSurface initialTSOS(firstTSOS.localParameters(), LocalTrajectoryError(CovLoc),
2978  firstTSOS.surface(), &*magneticField_);
2979 
2980  PTrajectoryStateOnDet PTraj =
2981  trajectoryStateTransform::persistentState(initialTSOS, trackDetId.rawId());
2982  //const TrajectorySeed seedT(*PTraj, recHit, alongMomentum);
2983  const TrajectorySeed seedT(PTraj, recHit, direction);
2984 
2985  std::vector<Trajectory> trajVec;
2986  if(direction == alongMomentum) trajVec = theFitter->fit(seedT, recHitMu, initialTSOS);
2987  else trajVec = theFitterOp->fit(seedT, recHitMu, initialTSOS);
2988 
2989  if(info){
2990  this->debugTrajectorySOSv(" innermostTSOS of TransientTrack",
2991  alongTTr.innermostMeasurementState());
2992  this->debugTrajectorySOSv(" outermostTSOS of TransientTrack",
2993  alongTTr.outermostMeasurementState());
2994  this->debugTrajectorySOS(" initialTSOS for theFitter ", initialTSOS);
2995  std::cout<<" . . . . . trajVec.size() "<<trajVec.size()<<std::endl;
2996  if(trajVec.size()) this->debugTrajectory(" theFitter trajectory ", trajVec[0]);
2997  }
2998 
2999  if(!smooth){
3000  if(trajVec.size()) trackFittedTSOS = trajVec[0].lastMeasurement().updatedState();}
3001  else{
3002  std::vector<Trajectory> trajSVec;
3003  if (trajVec.size()){
3004  if(direction == alongMomentum) trajSVec = theSmoother->trajectories(*(trajVec.begin()));
3005  else trajSVec = theSmootherOp->trajectories(*(trajVec.begin()));
3006  }
3007  if(info) std::cout<<" . . . . trajSVec.size() "<<trajSVec.size()<<std::endl;
3008  if(trajSVec.size()) this->debugTrajectorySOSv("smoothed geom InnermostState",
3009  trajSVec[0].geometricalInnermostState());
3010  if(trajSVec.size()) trackFittedTSOS = trajSVec[0].geometricalInnermostState();
3011  }
3012 
3013  if(info) this->debugTrajectorySOSv(" trackFittedTSOS ", trackFittedTSOS);
3014 
3015 } // end of trackFitter
i
int i
Definition: DBlmapReader.cc:9
GlobalTrackerMuonAlignment::debugTrajectorySOSv
void debugTrajectorySOSv(const std::string, TrajectoryStateOnSurface)
Definition: GlobalTrackerMuonAlignment.cc:3179
info
static const TGPicture * info(bool iBackgroundIsBlack)
Definition: FWCollectionSummaryWidget.cc:170
TrajectoryStateOnSurface::localParameters
const LocalTrajectoryParameters & localParameters() const
Definition: TrajectoryStateOnSurface.h:113
MuonTransientTrackingRecHitBuilder::ConstRecHitContainer
TransientTrackingRecHit::ConstRecHitContainer ConstRecHitContainer
Definition: MuonTransientTrackingRecHitBuilder.h:17
edm::OwnVector::size
size_type size() const
Definition: OwnVector.h:248
trajectoryStateTransform::persistentState
PTrajectoryStateOnDet persistentState(const TrajectoryStateOnSurface &ts, unsigned int detid)
Definition: TrajectoryStateTransform.cc:16
TrackingRecHit::RecHitContainer
std::vector< ConstRecHitPointer > RecHitContainer
Definition: TrackingRecHit.h:39
AlgebraicSymMatrix55
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
Definition: AlgebraicROOTObjects.h:31
reco::TransientTrack::innermostMeasurementState
TrajectoryStateOnSurface innermostMeasurementState() const
Definition: TransientTrack.h:81
DetId::rawId
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
edm::OwnVector::push_back
void push_back(D *&d)
Definition: OwnVector.h:274
TrajectoryStateOnSurface::surface
const SurfaceType & surface() const
Definition: TrajectoryStateOnSurface.h:128
GlobalTrackerMuonAlignment::debugTrackHit
void debugTrackHit(const std::string, reco::TrackRef)
Definition: GlobalTrackerMuonAlignment.cc:3131
TransientTrackingRecHitBuilder::build
virtual RecHitPointer build(const TrackingRecHit *p) const =0
build a tracking rechit from an existing rechit
edm::OwnVector::clear
void clear()
Definition: OwnVector.h:371
LocalTrajectoryError::matrix
const AlgebraicSymMatrix55 & matrix() const
Definition: LocalTrajectoryError.h:61
TrajectoryStateOnSurface::localError
const LocalTrajectoryError & localError() const
Definition: TrajectoryStateOnSurface.h:125
reco::TransientTrack::outermostMeasurementState
TrajectoryStateOnSurface outermostMeasurementState() const
Definition: TransientTrack.h:78
DetId
Definition: DetId.h:18
GlobalTrackerMuonAlignment::TTRHBuilder
const TransientTrackingRecHitBuilder * TTRHBuilder
Definition: GlobalTrackerMuonAlignment.cc:188
GlobalTrackerMuonAlignment::debugTrajectory
void debugTrajectory(const std::string, Trajectory &)
Definition: GlobalTrackerMuonAlignment.cc:3211
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
GlobalTrackerMuonAlignment::debugTrajectorySOS
void debugTrajectorySOS(const std::string, TrajectoryStateOnSurface &)
Definition: GlobalTrackerMuonAlignment.cc:3147
void GlobalTrackerMuonAlignment::writeGlPosRcd ( CLHEP::HepVector &  d3)

Definition at line 3227 of file GlobalTrackerMuonAlignment.cc.

References alignmentValidation::c1, DetId::Calo, funct::cos(), gather_cfg::cout, cond::service::PoolDBOutputService::currentTime(), DetId::Ecal, patCandidatesForDimuonsSequences_cff::ecal, DetId::Hcal, patCandidatesForDimuonsSequences_cff::hcal, edm::Service< T >::isAvailable(), Alignments::m_align, DetId::Muon, metsig::muon, DetId::rawId(), AlignTransform::rawId(), AlignTransform::rotation(), indexGen::s2, funct::sin(), DetId::Tracker, patCandidatesForDimuonsSequences_cff::tracker, AlignTransform::translation(), and cond::service::PoolDBOutputService::writeOne().

Referenced by endJob().

3228 {
3229  //paramVec(1) = 0.1; paramVec(2) = 0.2; paramVec(3) = 0.3; //0123
3230  //paramVec(4) = 0.0001; paramVec(5) = 0.0002; paramVec(6) = 0.0003;
3231  //paramVec(1) = 0.3; paramVec(2) = 0.4; paramVec(3) = 0.5; //03450123
3232  //paramVec(4) = 0.0001; paramVec(5) = 0.0002; paramVec(6) = 0.0003;
3233  //paramVec(1) = 2.; paramVec(2) = 2.; paramVec(3) = 2.; //222
3234  //paramVec(4) = 0.02; paramVec(5) = 0.02; paramVec(6) = 0.02;
3235  //paramVec(1) = -10.; paramVec(2) = -20.; paramVec(3) = -30.; //102030
3236  //paramVec(4) = -0.1; paramVec(5) = -0.2; paramVec(6) = -0.3;
3237  //paramVec(1) = 0.; paramVec(2) = 0.; paramVec(3) = 1.; //1z
3238  //paramVec(4) = 0.0000; paramVec(5) = 0.0000; paramVec(6) = 0.0000;
3239 
3240  std::cout<<" paramVector "<<paramVec.T()<<std::endl;
3241 
3242  CLHEP::Hep3Vector colX, colY, colZ;
3243 
3244 #ifdef NOT_EXACT_ROTATION_MATRIX
3245  colX = CLHEP::Hep3Vector( 1., -paramVec(6), paramVec(5));
3246  colY = CLHEP::Hep3Vector( paramVec(6), 1., -paramVec(4));
3247  colZ = CLHEP::Hep3Vector(-paramVec(5), paramVec(4), 1.);
3248 #else
3249  double s1 = std::sin(paramVec(4)), c1 = std::cos(paramVec(4));
3250  double s2 = std::sin(paramVec(5)), c2 = std::cos(paramVec(5));
3251  double s3 = std::sin(paramVec(6)), c3 = std::cos(paramVec(6));
3252  colX = CLHEP::Hep3Vector( c2 * c3, -c2 * s3, s2);
3253  colY = CLHEP::Hep3Vector(c1 * s3 + s1 * s2 * c3, c1 * c3 - s1 * s2 * s3, -s1 * c2);
3254  colZ = CLHEP::Hep3Vector(s1 * s3 - c1 * s2 * c3, s1 * c3 + c1 * s2 * s3, c1 * c2);
3255 #endif
3256 
3257  CLHEP::HepVector param0(6,0);
3258 
3259  Alignments* globalPositions = new Alignments();
3260 
3261  //Tracker
3262  AlignTransform tracker(iteratorTrackerRcd->translation(),
3263  iteratorTrackerRcd->rotation(),
3264  DetId(DetId::Tracker).rawId());
3265  // Muon
3266  CLHEP::Hep3Vector posMuGlRcd = iteratorMuonRcd->translation();
3267  CLHEP::HepRotation rotMuGlRcd = iteratorMuonRcd->rotation();
3268  CLHEP::HepEulerAngles angMuGlRcd = iteratorMuonRcd->rotation().eulerAngles();
3269  if(debug_)
3270  std::cout<<" Old muon Rcd Euler angles "<<angMuGlRcd.phi()<<" "<<angMuGlRcd.theta()
3271  <<" "<<angMuGlRcd.psi()<<std::endl;
3272  AlignTransform muon;
3273  if((angMuGlRcd.phi() == 0.) && (angMuGlRcd.theta() == 0.) && (angMuGlRcd.psi() == 0.) &&
3274  (posMuGlRcd.x() == 0.) && (posMuGlRcd.y() == 0.) && (posMuGlRcd.z() == 0.)){
3275  std::cout<<" New muon parameters are stored in Rcd"<<std::endl;
3276 
3277  AlignTransform muonNew(AlignTransform::Translation(paramVec(1),
3278  paramVec(2),
3279  paramVec(3)),
3280  AlignTransform::Rotation (colX,
3281  colY,
3282  colZ),
3283  DetId(DetId::Muon).rawId());
3284  muon = muonNew;
3285  }
3286  else if((paramVec(1) == 0.) && (paramVec(2) == 0.) && (paramVec(3) == 0.) &&
3287  (paramVec(4) == 0.) && (paramVec(5) == 0.) && (paramVec(6) == 0.)){
3288  std::cout<<" Old muon parameters are stored in Rcd"<<std::endl;
3289 
3290  AlignTransform muonNew(iteratorMuonRcd->translation(),
3291  iteratorMuonRcd->rotation(),
3292  DetId(DetId::Muon).rawId());
3293  muon = muonNew;
3294  }
3295  else{
3296  std::cout<<" New + Old muon parameters are stored in Rcd"<<std::endl;
3297  CLHEP::Hep3Vector posMuGlRcdThis = CLHEP::Hep3Vector(paramVec(1), paramVec(2), paramVec(3));
3298  CLHEP::HepRotation rotMuGlRcdThis = CLHEP::HepRotation(colX, colY, colZ);
3299  CLHEP::Hep3Vector posMuGlRcdNew = rotMuGlRcdThis * posMuGlRcd + posMuGlRcdThis;
3300  CLHEP::HepRotation rotMuGlRcdNew = rotMuGlRcdThis * rotMuGlRcd;
3301 
3302  AlignTransform muonNew(posMuGlRcdNew,
3303  rotMuGlRcdNew,
3304  DetId(DetId::Muon).rawId());
3305  muon = muonNew;
3306  }
3307 
3308  // Ecal
3309  AlignTransform ecal(iteratorEcalRcd->translation(),
3310  iteratorEcalRcd->rotation(),
3311  DetId(DetId::Ecal).rawId());
3312  // Hcal
3313  AlignTransform hcal(iteratorHcalRcd->translation(),
3314  iteratorHcalRcd->rotation(),
3315  DetId(DetId::Hcal).rawId());
3316  // Calo
3317  AlignTransform calo(AlignTransform::Translation(param0(1),
3318  param0(2),
3319  param0(3)),
3320  AlignTransform::EulerAngles(param0(4),
3321  param0(5),
3322  param0(6)),
3323  DetId(DetId::Calo).rawId());
3324 
3325  std::cout << "Tracker (" << tracker.rawId() << ") at " << tracker.translation()
3326  << " " << tracker.rotation().eulerAngles() << std::endl;
3327  std::cout << tracker.rotation() << std::endl;
3328 
3329  std::cout << "Muon (" << muon.rawId() << ") at " << muon.translation()
3330  << " " << muon.rotation().eulerAngles() << std::endl;
3331  std::cout << " rotations angles around x,y,z "
3332  << " ( " << -muon.rotation().zy() << " " << muon.rotation().zx()
3333  << " " << -muon.rotation().yx() << " )" << std::endl;
3334  std::cout << muon.rotation() << std::endl;
3335 
3336  std::cout << "Ecal (" << ecal.rawId() << ") at " << ecal.translation()
3337  << " " << ecal.rotation().eulerAngles() << std::endl;
3338  std::cout << ecal.rotation() << std::endl;
3339 
3340  std::cout << "Hcal (" << hcal.rawId() << ") at " << hcal.translation()
3341  << " " << hcal.rotation().eulerAngles() << std::endl;
3342  std::cout << hcal.rotation() << std::endl;
3343 
3344  std::cout << "Calo (" << calo.rawId() << ") at " << calo.translation()
3345  << " " << calo.rotation().eulerAngles() << std::endl;
3346  std::cout << calo.rotation() << std::endl;
3347 
3348  globalPositions->m_align.push_back(tracker);
3349  globalPositions->m_align.push_back(muon);
3350  globalPositions->m_align.push_back(ecal);
3351  globalPositions->m_align.push_back(hcal);
3352  globalPositions->m_align.push_back(calo);
3353 
3354  std::cout << "Uploading to the database..." << std::endl;
3355 
3356  edm::Service<cond::service::PoolDBOutputService> poolDbService;
3357 
3358  if (!poolDbService.isAvailable())
3359  throw cms::Exception("NotAvailable") << "PoolDBOutputService not available";
3360 
3361  // if (poolDbService->isNewTagRequest("GlobalPositionRcd")) {
3362 // poolDbService->createNewIOV<Alignments>(&(*globalPositions), poolDbService->endOfTime(), "GlobalPositionRcd");
3363 // } else {
3364 // poolDbService->appendSinceTime<Alignments>(&(*globalPositions), poolDbService->currentTime(), "GlobalPositionRcd");
3365 // }
3366  poolDbService->writeOne<Alignments>(&(*globalPositions),
3367  poolDbService->currentTime(),
3368  //poolDbService->beginOfTime(),
3369  "GlobalPositionRcd");
3370  std::cout << "done!" << std::endl;
3371 
3372  return;
3373 }
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
AlignTransform::Translation
CLHEP::Hep3Vector Translation
Definition: AlignTransform.h:19
GlobalTrackerMuonAlignment::iteratorHcalRcd
std::vector< AlignTransform >::const_iterator iteratorHcalRcd
Definition: GlobalTrackerMuonAlignment.cc:210
AlignTransform::EulerAngles
CLHEP::HepEulerAngles EulerAngles
Definition: AlignTransform.h:18
Alignments::m_align
std::vector< AlignTransform > m_align
Definition: Alignments.h:19
indexGen.s2
tuple s2
Definition: indexGen.py:106
DetId::rawId
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
AlignTransform::translation
const Translation & translation() const
Definition: AlignTransform.h:42
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
edm::Service::isAvailable
bool isAvailable() const
Definition: Service.h:46
GlobalTrackerMuonAlignment::iteratorMuonRcd
std::vector< AlignTransform >::const_iterator iteratorMuonRcd
Definition: GlobalTrackerMuonAlignment.cc:206
patCandidatesForDimuonsSequences_cff.tracker
tuple tracker
Definition: patCandidatesForDimuonsSequences_cff.py:69
cond::service::PoolDBOutputService::writeOne
void writeOne(T *payload, Time_t time, const std::string &recordName, bool withlogging=false)
Definition: PoolDBOutputService.h:58
GlobalTrackerMuonAlignment::iteratorEcalRcd
std::vector< AlignTransform >::const_iterator iteratorEcalRcd
Definition: GlobalTrackerMuonAlignment.cc:208
DetId
Definition: DetId.h:18
patCandidatesForDimuonsSequences_cff.ecal
tuple ecal
Definition: patCandidatesForDimuonsSequences_cff.py:70
AlignTransform::rotation
Rotation rotation() const
Definition: AlignTransform.h:47
GlobalTrackerMuonAlignment::iteratorTrackerRcd
std::vector< AlignTransform >::const_iterator iteratorTrackerRcd
Definition: GlobalTrackerMuonAlignment.cc:204
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
AlignTransform::rawId
align::ID rawId() const
Do not expose Euler angles since we may change its type later.
Definition: AlignTransform.h:45
patCandidatesForDimuonsSequences_cff.hcal
tuple hcal
Definition: patCandidatesForDimuonsSequences_cff.py:71
AlignTransform::Rotation
CLHEP::HepRotation Rotation
Definition: AlignTransform.h:20

Member Data Documentation

bool GlobalTrackerMuonAlignment::collectionCosmic
private

Definition at line 161 of file GlobalTrackerMuonAlignment.cc.

Referenced by beginJob(), and endJob().

bool GlobalTrackerMuonAlignment::collectionIsolated
private

Definition at line 162 of file GlobalTrackerMuonAlignment.cc.

Referenced by beginJob(), and endJob().

bool GlobalTrackerMuonAlignment::cosmicMuonMode_
private

Definition at line 159 of file GlobalTrackerMuonAlignment.cc.

Referenced by beginJob(), and GlobalTrackerMuonAlignment().

bool GlobalTrackerMuonAlignment::debug_
private

Definition at line 168 of file GlobalTrackerMuonAlignment.cc.

Referenced by endJob().

bool GlobalTrackerMuonAlignment::defineFitter
private

Definition at line 186 of file GlobalTrackerMuonAlignment.cc.

TFile* GlobalTrackerMuonAlignment::file
private

Definition at line 219 of file GlobalTrackerMuonAlignment.cc.

Referenced by beginJob(), and endJob().

AlgebraicVector3 GlobalTrackerMuonAlignment::Gfr
private

Definition at line 191 of file GlobalTrackerMuonAlignment.cc.

Referenced by beginJob(), and endJob().

const Alignments* GlobalTrackerMuonAlignment::globalPositionRcd_
private

Definition at line 180 of file GlobalTrackerMuonAlignment.cc.

edm::InputTag GlobalTrackerMuonAlignment::gmuonTags_
private

Definition at line 156 of file GlobalTrackerMuonAlignment.cc.

CLHEP::HepVector GlobalTrackerMuonAlignment::Grad
private

Definition at line 194 of file GlobalTrackerMuonAlignment.cc.

Referenced by beginJob(), and endJob().

CLHEP::HepVector GlobalTrackerMuonAlignment::GradL
private

Definition at line 197 of file GlobalTrackerMuonAlignment.cc.

Referenced by beginJob(), and endJob().

CLHEP::HepMatrix GlobalTrackerMuonAlignment::Hess
private

Definition at line 195 of file GlobalTrackerMuonAlignment.cc.

Referenced by beginJob(), and endJob().

CLHEP::HepMatrix GlobalTrackerMuonAlignment::HessL
private

Definition at line 198 of file GlobalTrackerMuonAlignment.cc.

Referenced by beginJob(), and endJob().

TH1F* GlobalTrackerMuonAlignment::histo
private

Definition at line 220 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist().

TH2F* GlobalTrackerMuonAlignment::histo101
private

Definition at line 255 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist().

TH2F* GlobalTrackerMuonAlignment::histo102
private

Definition at line 256 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist().

TH1F* GlobalTrackerMuonAlignment::histo11
private

Definition at line 229 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist().

TH1F* GlobalTrackerMuonAlignment::histo12
private

Definition at line 230 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist(), and fitHist().

TH1F* GlobalTrackerMuonAlignment::histo13
private

Definition at line 231 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist(), and fitHist().

TH1F* GlobalTrackerMuonAlignment::histo14
private

Definition at line 232 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist(), and fitHist().

TH1F* GlobalTrackerMuonAlignment::histo15
private

Definition at line 233 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist(), and fitHist().

TH1F* GlobalTrackerMuonAlignment::histo16
private

Definition at line 234 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist(), and fitHist().

TH1F* GlobalTrackerMuonAlignment::histo17
private

Definition at line 235 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist(), and fitHist().

TH1F* GlobalTrackerMuonAlignment::histo18
private

Definition at line 236 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist().

TH1F* GlobalTrackerMuonAlignment::histo19
private

Definition at line 237 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist().

TH1F* GlobalTrackerMuonAlignment::histo2
private

Definition at line 221 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist().

TH1F* GlobalTrackerMuonAlignment::histo20
private

Definition at line 238 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist().

TH1F* GlobalTrackerMuonAlignment::histo21
private

Definition at line 239 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist(), and fitHist().

TH1F* GlobalTrackerMuonAlignment::histo22
private

Definition at line 240 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist(), and fitHist().

TH1F* GlobalTrackerMuonAlignment::histo23
private

Definition at line 241 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist(), and fitHist().

TH1F* GlobalTrackerMuonAlignment::histo24
private

Definition at line 242 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist(), and fitHist().

TH1F* GlobalTrackerMuonAlignment::histo25
private

Definition at line 243 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist(), and fitHist().

TH1F* GlobalTrackerMuonAlignment::histo26
private

Definition at line 244 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist(), and fitHist().

TH1F* GlobalTrackerMuonAlignment::histo27
private

Definition at line 245 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist().

TH1F* GlobalTrackerMuonAlignment::histo28
private

Definition at line 246 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist().

TH1F* GlobalTrackerMuonAlignment::histo29
private

Definition at line 247 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist().

TH1F* GlobalTrackerMuonAlignment::histo3
private

Definition at line 222 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist().

TH1F* GlobalTrackerMuonAlignment::histo30
private

Definition at line 248 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist().

TH1F* GlobalTrackerMuonAlignment::histo31
private

Definition at line 249 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist().

TH1F* GlobalTrackerMuonAlignment::histo32
private

Definition at line 250 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist(), and fitHist().

TH1F* GlobalTrackerMuonAlignment::histo33
private

Definition at line 251 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist(), and fitHist().

TH1F* GlobalTrackerMuonAlignment::histo34
private

Definition at line 252 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist(), and fitHist().

TH1F* GlobalTrackerMuonAlignment::histo35
private

Definition at line 253 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist(), and fitHist().

TH1F* GlobalTrackerMuonAlignment::histo4
private

Definition at line 223 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist().

TH1F* GlobalTrackerMuonAlignment::histo5
private

Definition at line 224 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist().

TH1F* GlobalTrackerMuonAlignment::histo6
private

Definition at line 225 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist().

TH1F* GlobalTrackerMuonAlignment::histo7
private

Definition at line 226 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist(), and fitHist().

TH1F* GlobalTrackerMuonAlignment::histo8
private

Definition at line 227 of file GlobalTrackerMuonAlignment.cc.

Referenced by bookHist().

AlgebraicSymMatrix33 GlobalTrackerMuonAlignment::Inf
private

Definition at line 192 of file GlobalTrackerMuonAlignment.cc.

Referenced by beginJob(), and endJob().

bool GlobalTrackerMuonAlignment::isolatedMuonMode_
private

Definition at line 160 of file GlobalTrackerMuonAlignment.cc.

Referenced by beginJob(), and GlobalTrackerMuonAlignment().

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

Definition at line 208 of file GlobalTrackerMuonAlignment.cc.

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

Definition at line 210 of file GlobalTrackerMuonAlignment.cc.

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

Definition at line 206 of file GlobalTrackerMuonAlignment.cc.

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

Definition at line 204 of file GlobalTrackerMuonAlignment.cc.

const MagneticField* GlobalTrackerMuonAlignment::magneticField_
private

Definition at line 175 of file GlobalTrackerMuonAlignment.cc.

CLHEP::HepVector GlobalTrackerMuonAlignment::MuGlAngle
private

Definition at line 213 of file GlobalTrackerMuonAlignment.cc.

Referenced by endJob().

CLHEP::HepVector GlobalTrackerMuonAlignment::MuGlShift
private

Definition at line 212 of file GlobalTrackerMuonAlignment.cc.

Referenced by endJob().

edm::InputTag GlobalTrackerMuonAlignment::muonTags_
private

Definition at line 155 of file GlobalTrackerMuonAlignment.cc.

MuonTransientTrackingRecHitBuilder* GlobalTrackerMuonAlignment::MuRHBuilder
private

Definition at line 187 of file GlobalTrackerMuonAlignment.cc.

std::string GlobalTrackerMuonAlignment::MuSelect
private

Definition at line 215 of file GlobalTrackerMuonAlignment.cc.

Referenced by endJob().

int GlobalTrackerMuonAlignment::N_event
private

Definition at line 200 of file GlobalTrackerMuonAlignment.cc.

Referenced by beginJob(), and endJob().

int GlobalTrackerMuonAlignment::N_track
private

Definition at line 201 of file GlobalTrackerMuonAlignment.cc.

Referenced by beginJob(), and endJob().

ofstream GlobalTrackerMuonAlignment::OutGlobalTxt
private

Definition at line 217 of file GlobalTrackerMuonAlignment.cc.

Referenced by endJob().

string GlobalTrackerMuonAlignment::propagator_
private

Definition at line 158 of file GlobalTrackerMuonAlignment.cc.

bool GlobalTrackerMuonAlignment::refitMuon_
private

Definition at line 163 of file GlobalTrackerMuonAlignment.cc.

bool GlobalTrackerMuonAlignment::refitTrack_
private

Definition at line 164 of file GlobalTrackerMuonAlignment.cc.

string GlobalTrackerMuonAlignment::rootOutFile_
private

Definition at line 165 of file GlobalTrackerMuonAlignment.cc.

Referenced by beginJob().

edm::InputTag GlobalTrackerMuonAlignment::smuonTags_
private

Definition at line 157 of file GlobalTrackerMuonAlignment.cc.

Referenced by endJob().

KFTrajectoryFitter* GlobalTrackerMuonAlignment::theFitter
private

Definition at line 182 of file GlobalTrackerMuonAlignment.cc.

KFTrajectoryFitter* GlobalTrackerMuonAlignment::theFitterOp
private

Definition at line 184 of file GlobalTrackerMuonAlignment.cc.

KFTrajectorySmoother* GlobalTrackerMuonAlignment::theSmoother
private

Definition at line 183 of file GlobalTrackerMuonAlignment.cc.

KFTrajectorySmoother* GlobalTrackerMuonAlignment::theSmootherOp
private

Definition at line 185 of file GlobalTrackerMuonAlignment.cc.

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

Definition at line 171 of file GlobalTrackerMuonAlignment.cc.

const GlobalTrackingGeometry* GlobalTrackerMuonAlignment::trackingGeometry_
private

Definition at line 172 of file GlobalTrackerMuonAlignment.cc.

edm::InputTag GlobalTrackerMuonAlignment::trackTags_
private

Definition at line 154 of file GlobalTrackerMuonAlignment.cc.

const TransientTrackingRecHitBuilder* GlobalTrackerMuonAlignment::TTRHBuilder
private

Definition at line 188 of file GlobalTrackerMuonAlignment.cc.

string GlobalTrackerMuonAlignment::txtOutFile_
private

Definition at line 166 of file GlobalTrackerMuonAlignment.cc.

Referenced by endJob().

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

Definition at line 179 of file GlobalTrackerMuonAlignment.cc.

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

Definition at line 174 of file GlobalTrackerMuonAlignment.cc.

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

Definition at line 177 of file GlobalTrackerMuonAlignment.cc.

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

Definition at line 170 of file GlobalTrackerMuonAlignment.cc.

bool GlobalTrackerMuonAlignment::writeDB_
private

Definition at line 167 of file GlobalTrackerMuonAlignment.cc.

Referenced by endJob().

Generated for CMSSW Reference Manual by doxygen 1.8.5