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

Public Member Functions

void analyzeTrackTrack (const edm::Event &, const edm::EventSetup &)
 
void analyzeTrackTrajectory (const edm::Event &, const edm::EventSetup &)
 
void bookHist ()
 
double CLHEP_dot (CLHEP::HepVector a, 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
 EDAnalyzer ()
 
std::string workerType () const
 
virtual ~EDAnalyzer ()
 

Private Member Functions

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

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
 
typedef WorkerT< EDAnalyzer > WorkerType
 
- 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::EDAnalyzer
CurrentProcessingContext const * currentContext () const
 

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 
)
privatevirtual

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:1244
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(), DetId::Ecal, epsilon, TrajectoryStateOnSurface::freeState(), edm::EventSetup::get(), edm::Event::getByLabel(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), gradientGlobal(), gradientGlobalAlg(), gradientLocal(), DetId::Hcal, interpolateCardsSimple::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(), LargeD0_PixelPairStep_cff::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  RKTestPropagator alongRKPr =
692  RKTestPropagator(magneticField_, alongMomentum);
693  RKTestPropagator oppositeRKPr =
694  RKTestPropagator(magneticField_, oppositeToMomentum);
695 
696  float epsilon = 5.;
697  SmartPropagator alongSmPr =
698  SmartPropagator(alongRKPr, alongStHePr, magneticField_, alongMomentum, epsilon);
699  SmartPropagator oppositeSmPr =
700  SmartPropagator(oppositeRKPr, oppositeStHePr, magneticField_, oppositeToMomentum, epsilon);
701 
702  // ................................................ selected/global muon
703  //itMuon --> Jim's globalMuon
704  for(MuonCollection::const_iterator itMuon = smuons->begin();
705  itMuon != smuons->end(); ++itMuon) {
706 
707  if(debug_){
708  std::cout<<" mu gM is GM Mu SaM tM "<<itMuon->isGlobalMuon()<<" "
709  <<itMuon->isMuon()<<" "<<itMuon->isStandAloneMuon()
710  <<" "<<itMuon->isTrackerMuon()<<" "
711  <<std::endl;
712  }
713 
714  // get information about the innermost muon hit -------------------------
715  TransientTrack muTT(itMuon->outerTrack(), magneticField_, trackingGeometry_);
716  TrajectoryStateOnSurface innerMuTSOS = muTT.innermostMeasurementState();
717  TrajectoryStateOnSurface outerMuTSOS = muTT.outermostMeasurementState();
718 
719  // get information about the outermost tracker hit -----------------------
720  TransientTrack trackTT(itMuon->track(), magneticField_, trackingGeometry_);
721  TrajectoryStateOnSurface outerTrackTSOS = trackTT.outermostMeasurementState();
722  TrajectoryStateOnSurface innerTrackTSOS = trackTT.innermostMeasurementState();
723 
724  GlobalPoint pointTrackIn = innerTrackTSOS.globalPosition();
725  GlobalPoint pointTrackOut = outerTrackTSOS.globalPosition();
726  float lenghtTrack = (pointTrackOut-pointTrackIn).mag();
727  GlobalPoint pointMuonIn = innerMuTSOS.globalPosition();
728  GlobalPoint pointMuonOut = outerMuTSOS.globalPosition();
729  float lenghtMuon = (pointMuonOut - pointMuonIn).mag();
730  GlobalVector momentumTrackOut = outerTrackTSOS.globalMomentum();
731  GlobalVector momentumTrackIn = innerTrackTSOS.globalMomentum();
732  float distanceInIn = (pointTrackIn - pointMuonIn).mag();
733  float distanceInOut = (pointTrackIn - pointMuonOut).mag();
734  float distanceOutIn = (pointTrackOut - pointMuonIn).mag();
735  float distanceOutOut = (pointTrackOut - pointMuonOut).mag();
736 
737  if(debug_){
738  std::cout<<" pointTrackIn "<<pointTrackIn<<std::endl;
739  std::cout<<" Out "<<pointTrackOut<<" lenght "<<lenghtTrack<<std::endl;
740  std::cout<<" point MuonIn "<<pointMuonIn<<std::endl;
741  std::cout<<" Out "<<pointMuonOut<<" lenght "<<lenghtMuon<<std::endl;
742  std::cout<<" momeTrackIn Out "<<momentumTrackIn<<" "<<momentumTrackOut<<std::endl;
743  std::cout<<" doi io ii oo "<<distanceOutIn<<" "<<distanceInOut<<" "
744  <<distanceInIn<<" "<<distanceOutOut<<std::endl;
745  }
746 
747  if(lenghtTrack < 90.) continue;
748  if(lenghtMuon < 300.) continue;
749  if(momentumTrackIn.mag() < 15. || momentumTrackOut.mag() < 15.) continue;
750  if(trackTT.charge() != muTT.charge()) continue;
751 
752  if(debug_)
753  std::cout<<" passed lenght momentum cuts"<<std::endl;
754 
755  GlobalVector GRm, GPm, Nl, Rm, Pm, Rt, Pt, Rt0;
756  AlgebraicSymMatrix66 Cm, C0, Ce, C1;
757 
758  TrajectoryStateOnSurface extrapolationT;
759  int tsosMuonIf = 0;
760 
761  if( isolatedMuonMode_ ){ //------------------------------- Isolated Muon -----
762  const Surface& refSurface = innerMuTSOS.surface();
763  ReferenceCountingPointer<TangentPlane>
764  tpMuLocal(refSurface.tangentPlane(innerMuTSOS.localPosition()));
765  Nl = tpMuLocal->normalVector();
766  ReferenceCountingPointer<TangentPlane>
767  tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
768  GlobalVector Ng = tpMuGlobal->normalVector();
769  Surface* surf = (Surface*)&refSurface;
770  const Plane* refPlane = dynamic_cast<Plane*>(surf);
771  GlobalVector Nlp = refPlane->normalVector();
772 
773  if(debug_){
774  std::cout<<" Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
775  <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
776  std::cout<<" global "<<Ng.x()<<" "<<Ng.y()<<" "<<Ng.z()<<std::endl;
777  std::cout<<" lp "<<Nlp.x()<<" "<<Nlp.y()<<" "<<Nlp.z()<<std::endl;
778  //std::cout<<" Nlocal Nglobal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
779  // <<Ng.x()<<" "<<Ng.y()<<" "<<Ng.z()
780  //<<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
781  }
782 
783  // extrapolation to innermost muon hit
784  extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
785  //extrapolationT = propagator->propagate(outerTrackTSOS, refSurface);
786 
787  if(!extrapolationT.isValid()){
788  if(false & alarm)
789  std::cout<<" !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
790  //<<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
791  <<std::endl;
792  continue;
793  }
794  tsosMuonIf = 1;
795  Rt = GlobalVector((extrapolationT.globalPosition()).x(),
796  (extrapolationT.globalPosition()).y(),
797  (extrapolationT.globalPosition()).z());
798 
799  Pt = extrapolationT.globalMomentum();
800  // global parameters of muon
801  GRm = GlobalVector((innerMuTSOS.globalPosition()).x(),
802  (innerMuTSOS.globalPosition()).y(),
803  (innerMuTSOS.globalPosition()).z());
804  GPm = innerMuTSOS.globalMomentum();
805  Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
806  (outerTrackTSOS.globalPosition()).y(),
807  (outerTrackTSOS.globalPosition()).z());
808  Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() +
809  innerMuTSOS.cartesianError().matrix());
810  C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
811  Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
812  C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
813 
814  } // ------------------------------- end Isolated Muon -----
815 
816 
817  if( cosmicMuonMode_ ){ //------------------------------- Cosmic Muon -----
818 
819  if((distanceOutIn <= distanceInOut) & (distanceOutIn <= distanceInIn) &
820  (distanceOutIn <= distanceOutOut)){ // ----- Out - In ------
821  if(debug_) std::cout<<" ----- Out - In -----"<<std::endl;
822 
823  const Surface& refSurface = innerMuTSOS.surface();
824  ReferenceCountingPointer<TangentPlane>
825  tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
826  Nl = tpMuGlobal->normalVector();
827 
828  // extrapolation to innermost muon hit
829  extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
830  //extrapolationT = propagator->propagate(outerTrackTSOS, refSurface);
831 
832  if(!extrapolationT.isValid()){
833  if(false & alarm)
834  std::cout<<" !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
835  //<<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
836  <<std::endl;
837  continue;
838  }
839  if(debug_)
840  std::cout<<" extrapolationT.isValid "<<extrapolationT.isValid()<<std::endl;
841 
842  tsosMuonIf = 1;
843  Rt = GlobalVector((extrapolationT.globalPosition()).x(),
844  (extrapolationT.globalPosition()).y(),
845  (extrapolationT.globalPosition()).z());
846 
847  Pt = extrapolationT.globalMomentum();
848  // global parameters of muon
849  GRm = GlobalVector((innerMuTSOS.globalPosition()).x(),
850  (innerMuTSOS.globalPosition()).y(),
851  (innerMuTSOS.globalPosition()).z());
852  GPm = innerMuTSOS.globalMomentum();
853  Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
854  (outerTrackTSOS.globalPosition()).y(),
855  (outerTrackTSOS.globalPosition()).z());
856  Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() +
857  innerMuTSOS.cartesianError().matrix());
858  C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
859  Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
860  C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
861 
862  if(false & debug_){
863  //std::cout<<" ->propDir "<<propagator->propagationDirection()<<std::endl;
864  PropagationDirectionChooser Chooser;
865  std::cout<<" propDirCh "
866  <<Chooser.operator()(*outerTrackTSOS.freeState(), refSurface)
867  <<" Ch == along "<<(alongMomentum ==
868  Chooser.operator()(*outerTrackTSOS.freeState(), refSurface))
869  <<std::endl;
870  std::cout<<" --- Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
871  <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
872  }
873  } // enf of ---- Out - In -----
874 
875 
876  else if((distanceInOut <= distanceInIn) & (distanceInOut <= distanceOutIn) &
877  (distanceInOut <= distanceOutOut)){ // ----- In - Out ------
878  if(debug_) std::cout<<" ----- In - Out -----"<<std::endl;
879 
880  const Surface& refSurface = outerMuTSOS.surface();
881  ReferenceCountingPointer<TangentPlane>
882  tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
883  Nl = tpMuGlobal->normalVector();
884 
885  // extrapolation to outermost muon hit
886  extrapolationT = oppositeSmPr.propagate(innerTrackTSOS, refSurface);
887 
888  if(!extrapolationT.isValid()){
889  if(false & alarm)
890  std::cout<<" !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
891  <<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
892  <<std::endl;
893  continue;
894  }
895  if(debug_)
896  std::cout<<" extrapolationT.isValid "<<extrapolationT.isValid()<<std::endl;
897 
898  tsosMuonIf = 2;
899  Rt = GlobalVector((extrapolationT.globalPosition()).x(),
900  (extrapolationT.globalPosition()).y(),
901  (extrapolationT.globalPosition()).z());
902 
903  Pt = extrapolationT.globalMomentum();
904  // global parameters of muon
905  GRm = GlobalVector((outerMuTSOS.globalPosition()).x(),
906  (outerMuTSOS.globalPosition()).y(),
907  (outerMuTSOS.globalPosition()).z());
908  GPm = outerMuTSOS.globalMomentum();
909  Rt0 = GlobalVector((innerTrackTSOS.globalPosition()).x(),
910  (innerTrackTSOS.globalPosition()).y(),
911  (innerTrackTSOS.globalPosition()).z());
912  Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() +
913  outerMuTSOS.cartesianError().matrix());
914  C0 = AlgebraicSymMatrix66(innerTrackTSOS.cartesianError().matrix());
915  Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
916  C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
917 
918  if(false & debug_){
919  //std::cout<<" ->propDir "<<propagator->propagationDirection()<<std::endl;
920  PropagationDirectionChooser Chooser;
921  std::cout<<" propDirCh "
922  <<Chooser.operator()(*innerTrackTSOS.freeState(), refSurface)
923  <<" Ch == oppisite "<<(oppositeToMomentum ==
924  Chooser.operator()(*innerTrackTSOS.freeState(), refSurface))
925  <<std::endl;
926  std::cout<<" --- Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
927  <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
928  }
929  } // enf of ---- In - Out -----
930 
931  else if((distanceOutOut <= distanceInOut) & (distanceOutOut <= distanceInIn) &
932  (distanceOutOut <= distanceOutIn)){ // ----- Out - Out ------
933  if(debug_) std::cout<<" ----- Out - Out -----"<<std::endl;
934 
935  // reject: momentum of track has opposite direction to muon track
936  continue;
937 
938  const Surface& refSurface = outerMuTSOS.surface();
939  ReferenceCountingPointer<TangentPlane>
940  tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
941  Nl = tpMuGlobal->normalVector();
942 
943  // extrapolation to outermost muon hit
944  extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
945 
946  if(!extrapolationT.isValid()){
947  if(alarm)
948  std::cout<<" !!!!!!!!! Catastrophe Out-Out extrapolationT.isValid "
949  <<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
950  <<std::endl;
951  continue;
952  }
953  if(debug_)
954  std::cout<<" extrapolationT.isValid "<<extrapolationT.isValid()<<std::endl;
955 
956  tsosMuonIf = 2;
957  Rt = GlobalVector((extrapolationT.globalPosition()).x(),
958  (extrapolationT.globalPosition()).y(),
959  (extrapolationT.globalPosition()).z());
960 
961  Pt = extrapolationT.globalMomentum();
962  // global parameters of muon
963  GRm = GlobalVector((outerMuTSOS.globalPosition()).x(),
964  (outerMuTSOS.globalPosition()).y(),
965  (outerMuTSOS.globalPosition()).z());
966  GPm = outerMuTSOS.globalMomentum();
967  Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
968  (outerTrackTSOS.globalPosition()).y(),
969  (outerTrackTSOS.globalPosition()).z());
970  Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() +
971  outerMuTSOS.cartesianError().matrix());
972  C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
973  Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
974  C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
975 
976  if(debug_){
977  PropagationDirectionChooser Chooser;
978  std::cout<<" propDirCh "
979  <<Chooser.operator()(*outerTrackTSOS.freeState(), refSurface)
980  <<" Ch == along "<<(alongMomentum ==
981  Chooser.operator()(*outerTrackTSOS.freeState(), refSurface))
982  <<std::endl;
983  std::cout<<" --- Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
984  <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
985  std::cout<<" Nornal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
986  <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
987  }
988  } // enf of ---- Out - Out -----
989  else{
990  if(alarm)
991  std::cout<<" ------- !!!!!!!!!!!!!!! No proper Out - In \a\a\a\a\a\a\a"<<std::endl;
992  continue;
993  }
994 
995  } // ------------------------------- end Cosmic Muon -----
996 
997  if(tsosMuonIf == 0) {if(info) {std::cout<<"No tsosMuon !!!!!!"<<std::endl; continue;}}
998  TrajectoryStateOnSurface tsosMuon;
999  if(tsosMuonIf == 1) tsosMuon = muTT.innermostMeasurementState();
1000  else tsosMuon = muTT.outermostMeasurementState();
1001 
1002  //GlobalTrackerMuonAlignment::misalignMuon(GRm, GPm, Nl, Rt, Rm, Pm);
1003  AlgebraicVector4 LPRm; // muon local (dx/dz, dy/dz, x, y)
1004  GlobalTrackerMuonAlignment::misalignMuonL
1005  (GRm, GPm, Nl, Rt, Rm, Pm, LPRm, extrapolationT, tsosMuon);
1006 
1007  GlobalVector resR = Rm - Rt;
1008  GlobalVector resP0 = Pm - Pt;
1009  GlobalVector resP = Pm / Pm.mag() - Pt / Pt.mag();
1010  float RelMomResidual = (Pm.mag() - Pt.mag()) / (Pt.mag() + 1.e-6);;
1011 
1012  AlgebraicVector6 Vm;
1013  Vm(0) = resR.x(); Vm(1) = resR.y(); Vm(2) = resR.z();
1014  Vm(3) = resP0.x(); Vm(4) = resP0.y(); Vm(5) = resP0.z();
1015  float Rmuon = Rm.perp();
1016  float Zmuon = Rm.z();
1017  float alfa_x = atan2(Nl.x(),Nl.y())*57.29578;
1018 
1019  if(debug_){
1020  std::cout<<" Nx Ny Nz alfa_x "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "<<alfa_x<<std::endl;
1021  //std::cout<<" Rm "<<Rm<<std::endl<<" Rt "<<Rt<<std::endl<<" resR "<<resR<<std::endl
1022  // <<" resP "<<resP<<" dp/p "<<RelMomResidual<<std::endl;
1023  }
1024 
1025  double chi_d = 0;
1026  for(int i=0; i<=5; i++) chi_d += Vm(i)*Vm(i)/Cm(i,i);
1027 
1028  AlgebraicVector5 Vml(tsosMuon.localParameters().vector()
1029  - extrapolationT.localParameters().vector());
1030  AlgebraicSymMatrix55 m(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
1031  AlgebraicSymMatrix55 Cml(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
1032  bool ierrLoc = !m.Invert();
1033  if (ierrLoc && debug_ && info) {
1034  std::cout<< " ==== Error inverting Local covariance matrix ==== "<<std::cout;
1035  continue;}
1036  double chi_Loc = ROOT::Math::Similarity(Vml,m);
1037  if(debug_)
1038  std::cout<<" chi_Loc px/pz/err "<<chi_Loc<<" "<<Vml(1)/sqrt(Cml(1,1))<<std::endl;
1039 
1040  if(Pt.mag() < 15.) continue;
1041  if(Pm.mag() < 5.) continue;
1042 
1043  //if(Pt.mag() < 30.) continue; // momenum cut < 30GeV
1044  //if(Pt.mag() < 60.) continue; // momenum cut < 30GeV
1045  //if(Pt.mag() > 50.) continue; // momenum cut > 50GeV
1046  //if(Pt.mag() > 100.) continue; // momenum cut > 100GeV
1047  //if(trackTT.charge() < 0) continue; // select positive charge
1048  //if(trackTT.charge() > 0) continue; // select negative charge
1049 
1050  //if(fabs(resR.x()) > 5.) continue; // strong cut X
1051  //if(fabs(resR.y()) > 5.) continue; // Y
1052  //if(fabs(resR.z()) > 5.) continue; // Z
1053  //if(fabs(resR.mag()) > 7.5) continue; // dR
1054 
1055  /*
1056  //if(fabs(RelMomResidual) > 0.5) continue;
1057  if(fabs(resR.x()) > 20.) continue;
1058  if(fabs(resR.y()) > 20.) continue;
1059  if(fabs(resR.z()) > 20.) continue;
1060  if(fabs(resR.mag()) > 30.) continue;
1061  if(fabs(resP.x()) > 0.06) continue;
1062  if(fabs(resP.y()) > 0.06) continue;
1063  if(fabs(resP.z()) > 0.06) continue;
1064  if(chi_d > 40.) continue;
1065  */
1066 
1067  // select Barrel
1068  //if(Rmuon < 400. || Rmuon > 450.) continue;
1069  //if(Zmuon < -600. || Zmuon > 600.) continue;
1070  //if(fabs(Nl.z()) > 0.95) continue;
1071  //MuSelect = " Barrel";
1072  // EndCap1
1073  //if(Rmuon < 120. || Rmuon > 450.) continue;
1074  //if(Zmuon < -720.) continue;
1075  //if(Zmuon > -580.) continue;
1076  //if(fabs(Nl.z()) < 0.95) continue;
1077  //MuSelect = " EndCap1";
1078  // EndCap2
1079  //if(Rmuon < 120. || Rmuon > 450.) continue;
1080  //if(Zmuon > 720.) continue;
1081  //if(Zmuon < 580.) continue;
1082  //if(fabs(Nl.z()) < 0.95) continue;
1083  //MuSelect = " EndCap2";
1084  // select All
1085  if(Rmuon < 120. || Rmuon > 450.) continue;
1086  if(Zmuon < -720. || Zmuon > 720.) continue;
1087  MuSelect = " Barrel+EndCaps";
1088 
1089 
1090  if(debug_)
1091  std::cout<<" .............. passed all cuts"<<std::endl;
1092 
1093  N_track++;
1094  // gradient and Hessian for each track
1095 
1096  GlobalTrackerMuonAlignment::gradientGlobalAlg(Rt, Pt, Rm, Nl, Cm);
1097  GlobalTrackerMuonAlignment::gradientGlobal(Rt, Pt, Rm, Pm, Nl, Cm);
1098 
1099  CLHEP::HepSymMatrix covLoc(4,0);
1100  for(int i=1; i<=4; i++)
1101  for(int j=1; j<=i; j++){
1102  covLoc(i,j) = (tsosMuon.localError().matrix()
1103  + extrapolationT.localError().matrix())(i,j);
1104  //if(i != j) Cov(i,j) = 0.;
1105  }
1106 
1107  const Surface& refSurface = tsosMuon.surface();
1108  CLHEP::HepMatrix rotLoc (3,3,0);
1109  rotLoc(1,1) = refSurface.rotation().xx();
1110  rotLoc(1,2) = refSurface.rotation().xy();
1111  rotLoc(1,3) = refSurface.rotation().xz();
1112 
1113  rotLoc(2,1) = refSurface.rotation().yx();
1114  rotLoc(2,2) = refSurface.rotation().yy();
1115  rotLoc(2,3) = refSurface.rotation().yz();
1116 
1117  rotLoc(3,1) = refSurface.rotation().zx();
1118  rotLoc(3,2) = refSurface.rotation().zy();
1119  rotLoc(3,3) = refSurface.rotation().zz();
1120 
1121  CLHEP::HepVector posLoc(3);
1122  posLoc(1) = refSurface.position().x();
1123  posLoc(2) = refSurface.position().y();
1124  posLoc(3) = refSurface.position().z();
1125 
1126  GlobalTrackerMuonAlignment::gradientLocal(Rt, Pt, Rm, Pm, Nl, covLoc, rotLoc, posLoc, LPRm);
1127 
1128  if(debug_){
1129  std::cout<<" Norm "<<Nl<<std::endl;
1130  std::cout<<" posLoc "<<posLoc.T()<<std::endl;
1131  std::cout<<" rotLoc "<<rotLoc<<std::endl;
1132  }
1133 
1134  // ----------------------------------------------------- fill histogram
1135  histo->Fill(itMuon->track()->pt());
1136 
1137  //histo2->Fill(itMuon->track()->outerP());
1138  histo2->Fill(Pt.mag());
1139  histo3->Fill((PI/2.-itMuon->track()->outerTheta()));
1140  histo4->Fill(itMuon->track()->phi());
1141  histo5->Fill(Rmuon);
1142  histo6->Fill(Zmuon);
1143  histo7->Fill(RelMomResidual);
1144  //histo8->Fill(chi);
1145  histo8->Fill(chi_d);
1146 
1147  histo101->Fill(Zmuon, Rmuon);
1148  histo101->Fill(Rt0.z(), Rt0.perp());
1149  histo102->Fill(Rt0.x(), Rt0.y());
1150  histo102->Fill(Rm.x(), Rm.y());
1151 
1152  histo11->Fill(resR.mag());
1153  if(fabs(Nl.x()) < 0.98) histo12->Fill(resR.x());
1154  if(fabs(Nl.y()) < 0.98) histo13->Fill(resR.y());
1155  if(fabs(Nl.z()) < 0.98) histo14->Fill(resR.z());
1156  histo15->Fill(resP.x());
1157  histo16->Fill(resP.y());
1158  histo17->Fill(resP.z());
1159 
1160  if((fabs(Nl.x()) < 0.98) && (fabs(Nl.y()) < 0.98) &&(fabs(Nl.z()) < 0.98))
1161  {
1162  histo18->Fill(sqrt(C0(0,0)));
1163  histo19->Fill(sqrt(C1(0,0)));
1164  histo20->Fill(sqrt(C1(0,0)+Ce(0,0)));
1165  }
1166  if(fabs(Nl.x()) < 0.98) histo21->Fill(Vm(0)/sqrt(Cm(0,0)));
1167  if(fabs(Nl.y()) < 0.98) histo22->Fill(Vm(1)/sqrt(Cm(1,1)));
1168  if(fabs(Nl.z()) < 0.98) histo23->Fill(Vm(2)/sqrt(Cm(2,2)));
1169  histo24->Fill(Vm(3)/sqrt(C1(3,3)+Ce(3,3)));
1170  histo25->Fill(Vm(4)/sqrt(C1(4,4)+Ce(4,4)));
1171  histo26->Fill(Vm(5)/sqrt(C1(5,5)+Ce(5,5)));
1172  histo27->Fill(Nl.x());
1173  histo28->Fill(Nl.y());
1174  histo29->Fill(lenghtTrack);
1175  histo30->Fill(lenghtMuon);
1176  histo31->Fill(chi_Loc);
1177  histo32->Fill(Vml(1)/sqrt(Cml(1,1)));
1178  histo33->Fill(Vml(2)/sqrt(Cml(2,2)));
1179  histo34->Fill(Vml(3)/sqrt(Cml(3,3)));
1180  histo35->Fill(Vml(4)/sqrt(Cml(4,4)));
1181 
1182  if (debug_) { //--------------------------------- debug print ----------
1183 
1184  std::cout<<" diag 0[ "<<C0(0,0)<<" "<<C0(1,1)<<" "<<C0(2,2)<<" "<<C0(3,3)<<" "
1185  <<C0(4,4)<<" "<<C0(5,5)<<" ]"<<std::endl;
1186  std::cout<<" diag e[ "<<Ce(0,0)<<" "<<Ce(1,1)<<" "<<Ce(2,2)<<" "<<Ce(3,3)<<" "
1187  <<Ce(4,4)<<" "<<Ce(5,5)<<" ]"<<std::endl;
1188  std::cout<<" diag 1[ "<<C1(0,0)<<" "<<C1(1,1)<<" "<<C1(2,2)<<" "<<C1(3,3)<<" "
1189  <<C1(4,4)<<" "<<C1(5,5)<<" ]"<<std::endl;
1190  std::cout<<" Rm "<<Rm.x()<<" "<<Rm.y()<<" "<<Rm.z()
1191  <<" Pm "<<Pm.x()<<" "<<Pm.y()<<" "<<Pm.z()<<std::endl;
1192  std::cout<<" Rt "<<Rt.x()<<" "<<Rt.y()<<" "<<Rt.z()
1193  <<" Pt "<<Pt.x()<<" "<<Pt.y()<<" "<<Pt.z()<<std::endl;
1194  std::cout<<" Nl*(Rm-Rt) "<<Nl.dot(Rm-Rt)<<std::endl;
1195  std::cout<<" resR "<<resR.x()<<" "<<resR.y()<<" "<<resR.z()
1196  <<" resP "<<resP.x()<<" "<<resP.y()<<" "<<resP.z()<<std::endl;
1197  std::cout<<" Rm-t "<<(Rm-Rt).x()<<" "<<(Rm-Rt).y()<<" "<<(Rm-Rt).z()
1198  <<" Pm-t "<<(Pm-Pt).x()<<" "<<(Pm-Pt).y()<<" "<<(Pm-Pt).z()<<std::endl;
1199  std::cout<<" Vm "<<Vm<<std::endl;
1200  std::cout<<" +- "<<sqrt(Cm(0,0))<<" "<<sqrt(Cm(1,1))<<" "<<sqrt(Cm(2,2))<<" "
1201  <<sqrt(Cm(3,3))<<" "<<sqrt(Cm(4,4))<<" "<<sqrt(Cm(5,5))<<std::endl;
1202  std::cout<<" Pmuon Ptrack dP/Ptrack "<<itMuon->outerTrack()->p()<<" "
1203  <<itMuon->track()->outerP()<<" "
1204  <<(itMuon->outerTrack()->p() -
1205  itMuon->track()->outerP())/itMuon->track()->outerP()<<std::endl;
1206  std::cout<<" cov matrix "<<std::endl;
1207  std::cout<<Cm<<std::endl;
1208  std::cout<<" diag [ "<<Cm(0,0)<<" "<<Cm(1,1)<<" "<<Cm(2,2)<<" "<<Cm(3,3)<<" "
1209  <<Cm(4,4)<<" "<<Cm(5,5)<<" ]"<<std::endl;
1210 
1211  static AlgebraicSymMatrix66 Ro;
1212  double Diag[6];
1213  for(int i=0; i<=5; i++) Diag[i] = sqrt(Cm(i,i));
1214  for(int i=0; i<=5; i++)
1215  for(int j=0; j<=5; j++)
1216  Ro(i,j) = Cm(i,j)/Diag[i]/Diag[j];
1217  std::cout<<" correlation matrix "<<std::endl;
1218  std::cout<<Ro<<std::endl;
1219 
1220  AlgebraicSymMatrix66 CmI;
1221  for(int i=0; i<=5; i++)
1222  for(int j=0; j<=5; j++)
1223  CmI(i,j) = Cm(i,j);
1224 
1225  bool ierr = !CmI.Invert();
1226  if( ierr ) { if(alarm || debug_)
1227  std::cout<<" Error inverse covariance matrix !!!!!!!!!!!"<<std::endl;
1228  continue;
1229  }
1230  std::cout<<" inverse cov matrix "<<std::endl;
1231  std::cout<<Cm<<std::endl;
1232 
1233  double chi = ROOT::Math::Similarity(Vm, CmI);
1234  std::cout<<" chi chi_d "<<chi<<" "<<chi_d<<std::endl;
1235  } // end of debug_ printout --------------------------------------------
1236 
1237  } // end loop on selected muons, i.e. Jim's globalMuon
1238 
1239 } //end of analyzeTrackTrack
TkRotation::xx
T xx() const
Definition: newTkRotation.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
GlobalTrackerMuonAlignment::gradientGlobalAlg
void gradientGlobalAlg(GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, AlgebraicSymMatrix66 &)
Definition: GlobalTrackerMuonAlignment.cc:2016
GlobalTrackerMuonAlignment::gradientLocal
void gradientLocal(GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, CLHEP::HepSymMatrix &, CLHEP::HepMatrix &, CLHEP::HepVector &, AlgebraicVector4 &)
Definition: GlobalTrackerMuonAlignment.cc:2268
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:71
TrajectoryStateOnSurface::localParameters
const LocalTrajectoryParameters & localParameters() const
Definition: TrajectoryStateOnSurface.h:167
PI
#define PI
Definition: HcalBeamMonitor.cc:14
mag
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
Definition: Basic3DVectorLD.h:142
Plane::normalVector
GlobalVector normalVector() const
Definition: Plane.h:47
AlgebraicSymMatrix66
ROOT::Math::SMatrix< double, 6, 6, ROOT::Math::MatRepSym< double, 6 > > AlgebraicSymMatrix66
Definition: AlgebraicROOTObjects.h:22
TrajectoryStateOnSurface::cartesianError
const CartesianTrajectoryError cartesianError() const
Definition: TrajectoryStateOnSurface.h:161
PV3DBase::y
T y() const
Definition: PV3DBase.h:62
TkRotation::yx
T yx() const
Definition: newTkRotation.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:143
Vector3DBase::dot
PreciseFloatType< T, U >::Type dot(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:107
LargeD0_PixelPairStep_cff.propagator
tuple propagator
Definition: LargeD0_PixelPairStep_cff.py:84
AlgebraicSymMatrix55
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
Definition: AlgebraicROOTObjects.h:21
GlobalTrackerMuonAlignment::gradientGlobal
void gradientGlobal(GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, AlgebraicSymMatrix66 &)
Definition: GlobalTrackerMuonAlignment.cc:2095
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:105
Alignments::m_align
std::vector< AlignTransform > m_align
Definition: Alignments.h:14
AlgebraicVector6
ROOT::Math::SVector< double, 6 > AlgebraicVector6
Definition: AlgebraicROOTObjects.h:14
Plane
Definition: Plane.h:17
detailsBasic3DVector::z
double double double z
Definition: newBasic3DVector.h:17
DetId::rawId
uint32_t rawId() const
get the raw id
Definition: DetId.h:45
TkRotation::zx
T zx() const
Definition: newTkRotation.h:257
TkRotation::xy
T xy() const
Definition: newTkRotation.h:252
LocalTrajectoryParameters::vector
AlgebraicVector5 vector() const
Definition: LocalTrajectoryParameters.h:120
TkRotation::zz
T zz() const
Definition: newTkRotation.h:259
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:66
SteppingHelixPropagator_cfi.SteppingHelixPropagator
tuple SteppingHelixPropagator
Definition: SteppingHelixPropagator_cfi.py:3
TrajectoryStateOnSurface::freeState
FreeTrajectoryState * freeState(bool withErrors=true) const
Definition: TrajectoryStateOnSurface.h:130
GlobalTrackerMuonAlignment::watchTrackingGeometry_
edm::ESWatcher< GlobalTrackingGeometryRecord > watchTrackingGeometry_
Definition: GlobalTrackerMuonAlignment.cc:170
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
PV3DBase::z
T z() const
Definition: PV3DBase.h:63
edm::Event::run
RunNumber_t run() const
Definition: Event.h:67
GlobalTrackerMuonAlignment::iteratorMuonRcd
std::vector< AlignTransform >::const_iterator iteratorMuonRcd
Definition: GlobalTrackerMuonAlignment.cc:206
TkRotation::zy
T zy() const
Definition: newTkRotation.h:258
j
int j
Definition: DBlmapReader.cc:9
LocalTrajectoryError::matrix
const AlgebraicSymMatrix55 & matrix() const
Definition: LocalTrajectoryError.h:60
TrajectoryStateOnSurface::localError
const LocalTrajectoryError & localError() const
Definition: TrajectoryStateOnSurface.h:179
TkRotation::yy
T yy() const
Definition: newTkRotation.h:255
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:356
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:2645
GlobalTrackerMuonAlignment::iteratorEcalRcd
std::vector< AlignTransform >::const_iterator iteratorEcalRcd
Definition: GlobalTrackerMuonAlignment.cc:208
DetId
Definition: DetId.h:20
Surface::tangentPlane
virtual ReferenceCountingPointer< TangentPlane > tangentPlane(const GlobalPoint &) const =0
AlgebraicVector5
ROOT::Math::SVector< double, 5 > AlgebraicVector5
Definition: AlgebraicROOTObjects.h:13
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:59
TrajectoryStateOnSurface::globalMomentum
GlobalVector globalMomentum() const
Definition: TrajectoryStateOnSurface.h:146
patZpeak.muons
tuple muons
Definition: patZpeak.py:38
TkRotation::xz
T xz() const
Definition: newTkRotation.h:253
TrajectoryStateOnSurface::surface
const Surface & surface() const
Definition: TrajectoryStateOnSurface.h:182
GloballyPositioned::rotation
const RotationType & rotation() const
Definition: GloballyPositioned.h:43
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:12
PV3DBase::x
T x() const
Definition: PV3DBase.h:61
GloballyPositioned::position
const PositionType & position() const
Definition: GloballyPositioned.h:41
SmartPropagator_cfi.SmartPropagator
tuple SmartPropagator
Definition: SmartPropagator_cfi.py:3
TkRotation::yz
T yz() const
Definition: newTkRotation.h:256
GlobalVector
Global3DVector GlobalVector
Definition: GlobalVector.h:10
void GlobalTrackerMuonAlignment::analyzeTrackTrajectory ( const edm::Event &  iEvent,
const edm::EventSetup &  iSetup 
)

Definition at line 1244 of file GlobalTrackerMuonAlignment.cc.

References alongMomentum, edm::EventBase::bunchCrossing(), TrajectoryStateOnSurface::cartesianError(), Chi2MeasurementEstimatorESProducer_cfi::Chi2MeasurementEstimator, gather_cfg::cout, Vector3DBase< T, FrameTag >::dot(), DetId::Ecal, epsilon, TrajectoryStateOnSurface::freeState(), edm::EventSetup::get(), edm::Event::getByLabel(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), gradientGlobal(), gradientGlobalAlg(), gradientLocal(), DetId::Hcal, interpolateCardsSimple::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(), LargeD0_PixelPairStep_cff::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, RSFinalFitAnalytical_cff::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().

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

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 ( CLHEP::HepVector  a,
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 3135 of file GlobalTrackerMuonAlignment.cc.

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

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

Definition at line 3118 of file GlobalTrackerMuonAlignment.cc.

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

3119 {
3120  std::cout<<" ------- "<<title<<" --------"<<std::endl;
3121  int nHit = 1;
3122  TransientTrackingRecHit::RecHitContainer recHit;
3123  for (trackingRecHit_iterator i=alongTr.recHitsBegin();i!=alongTr.recHitsEnd(); i++) {
3124  std::cout<<" Hit "<<nHit++<<" DetId "<<(*i)->geographicalId().det();
3125  if( (*i)->geographicalId().det() == DetId::Tracker ) std::cout<<" Tracker ";
3126  else if ( (*i)->geographicalId().det() == DetId::Muon ) std::cout<<" Muon ";
3127  else std::cout<<" Unknown ";
3128  if(!(*i)->isValid()) std::cout<<" not valid "<<std::endl;
3129  else std::cout<<std::endl;
3130  }
3131 }
i
int i
Definition: DBlmapReader.cc:9
TransientTrackingRecHit::RecHitContainer
std::vector< ConstRecHitPointer > RecHitContainer
Definition: TransientTrackingRecHit.h:29
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 3215 of file GlobalTrackerMuonAlignment.cc.

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

3216 {
3217  std::cout<<"\n"<<" ...... "<<title<<" ...... "<<std::endl;
3218  if(!traj.isValid()) {std::cout<<" Not valid !!!!!!!! "<<std::endl; return;}
3219  std::cout<<" chi2/Nhit "<<traj.chiSquared()<<" / "<<traj.foundHits();
3220  if(traj.direction() == alongMomentum) std::cout<<" alongMomentum >>>>"<<std::endl;
3221  else std::cout<<" oppositeToMomentum <<<<"<<std::endl;
3222  this->debugTrajectorySOSv(" firstMeasurementTSOS ",traj.firstMeasurement().updatedState());
3223  //this->debugTrajectorySOSv(" firstMeasurementTSOS ",traj.firstMeasurement().predictedState());
3224  this->debugTrajectorySOSv(" lastMeasurementTSOS ",traj.lastMeasurement().updatedState());
3225  //this->debugTrajectorySOSv(" geom InnermostState", traj.geometricalInnermostState());
3226  std::cout<<" . . . . . . . . . . . . . . . . . . . . . . . . . . . . \n"<<std::endl;
3227 }
Trajectory::foundHits
int foundHits() const
Definition: Trajectory.h:224
GlobalTrackerMuonAlignment::debugTrajectorySOSv
void debugTrajectorySOSv(const std::string, TrajectoryStateOnSurface)
Definition: GlobalTrackerMuonAlignment.cc:3183
Trajectory::direction
PropagationDirection const & direction() const
Definition: Trajectory.cc:195
Trajectory::lastMeasurement
TrajectoryMeasurement const & lastMeasurement() const
Definition: Trajectory.h:181
TrajectoryMeasurement::updatedState
TrajectoryStateOnSurface updatedState() const
Definition: TrajectoryMeasurement.h:187
Trajectory::isValid
bool isValid() const
Definition: Trajectory.h:259
Trajectory::firstMeasurement
TrajectoryMeasurement const & firstMeasurement() const
Definition: Trajectory.h:194
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
Trajectory::chiSquared
double chiSquared() const
Definition: Trajectory.h:242
void GlobalTrackerMuonAlignment::debugTrajectorySOS ( const std::string  title,
TrajectoryStateOnSurface &  trajSOS 
)

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

3153 {
3154  std::cout<<" --- "<<title<<" --- "<<std::endl;
3155  if(!trajSOS.isValid()) {std::cout<<" Not valid !!!! "<<std::endl; return;}
3156  std::cout<<" R |p| "<<trajSOS.globalPosition().perp()<<" "
3157  <<trajSOS.globalMomentum().mag()<<" charge "<<trajSOS.charge()<<std::endl;
3158  std::cout<<" x p "<<trajSOS.globalParameters().vector()(0)<<" "
3159  <<trajSOS.globalParameters().vector()(1)<<" "
3160  <<trajSOS.globalParameters().vector()(2)<<" "
3161  <<trajSOS.globalParameters().vector()(3)<<" "
3162  <<trajSOS.globalParameters().vector()(4)<<" "
3163  <<trajSOS.globalParameters().vector()(5)<<std::endl;
3164  std::cout<<" +/- "<<sqrt(trajSOS.cartesianError().matrix()(0,0))<<" "
3165  <<sqrt(trajSOS.cartesianError().matrix()(1,1))<<" "
3166  <<sqrt(trajSOS.cartesianError().matrix()(2,2))<<" "
3167  <<sqrt(trajSOS.cartesianError().matrix()(3,3))<<" "
3168  <<sqrt(trajSOS.cartesianError().matrix()(4,4))<<" "
3169  <<sqrt(trajSOS.cartesianError().matrix()(5,5))<<std::endl;
3170  std::cout<<" q/p dxy/dz xy "<<trajSOS.localParameters().vector()(0)<<" "
3171  <<trajSOS.localParameters().vector()(1)<<" "
3172  <<trajSOS.localParameters().vector()(2)<<" "
3173  <<trajSOS.localParameters().vector()(3)<<" "
3174  <<trajSOS.localParameters().vector()(4)<<std::endl;
3175  std::cout<<" +/- error "<<sqrt(trajSOS.localError().matrix()(0,0))<<" "
3176  <<sqrt(trajSOS.localError().matrix()(1,1))<<" "
3177  <<sqrt(trajSOS.localError().matrix()(2,2))<<" "
3178  <<sqrt(trajSOS.localError().matrix()(3,3))<<" "
3179  <<sqrt(trajSOS.localError().matrix()(4,4))<<" "<<std::endl;
3180 }
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:71
TrajectoryStateOnSurface::localParameters
const LocalTrajectoryParameters & localParameters() const
Definition: TrajectoryStateOnSurface.h:167
TrajectoryStateOnSurface::cartesianError
const CartesianTrajectoryError cartesianError() const
Definition: TrajectoryStateOnSurface.h:161
TrajectoryStateOnSurface::globalPosition
GlobalPoint globalPosition() const
Definition: TrajectoryStateOnSurface.h:143
LocalTrajectoryParameters::vector
AlgebraicVector5 vector() const
Definition: LocalTrajectoryParameters.h:120
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:66
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
LocalTrajectoryError::matrix
const AlgebraicSymMatrix55 & matrix() const
Definition: LocalTrajectoryError.h:60
TrajectoryStateOnSurface::localError
const LocalTrajectoryError & localError() const
Definition: TrajectoryStateOnSurface.h:179
CartesianTrajectoryError::matrix
const AlgebraicSymMatrix66 & matrix() const
Definition: CartesianTrajectoryError.h:29
TrajectoryStateOnSurface::globalParameters
const GlobalTrajectoryParameters & globalParameters() const
Definition: TrajectoryStateOnSurface.h:140
TrajectoryStateOnSurface::globalMomentum
GlobalVector globalMomentum() const
Definition: TrajectoryStateOnSurface.h:146
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
GlobalTrajectoryParameters::vector
AlgebraicVector6 vector() const
Definition: GlobalTrajectoryParameters.h:100
void GlobalTrackerMuonAlignment::debugTrajectorySOSv ( const std::string  title,
TrajectoryStateOnSurface  trajSOS 
)

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

3185 {
3186  std::cout<<" --- "<<title<<" --- "<<std::endl;
3187  if(!trajSOS.isValid()) {std::cout<<" Not valid !!!! "<<std::endl; return;}
3188  std::cout<<" R |p| "<<trajSOS.globalPosition().perp()<<" "
3189  <<trajSOS.globalMomentum().mag()<<" charge "<<trajSOS.charge()<<std::endl;
3190  std::cout<<" x p "<<trajSOS.globalParameters().vector()(0)<<" "
3191  <<trajSOS.globalParameters().vector()(1)<<" "
3192  <<trajSOS.globalParameters().vector()(2)<<" "
3193  <<trajSOS.globalParameters().vector()(3)<<" "
3194  <<trajSOS.globalParameters().vector()(4)<<" "
3195  <<trajSOS.globalParameters().vector()(5)<<std::endl;
3196  std::cout<<" +/- "<<sqrt(trajSOS.cartesianError().matrix()(0,0))<<" "
3197  <<sqrt(trajSOS.cartesianError().matrix()(1,1))<<" "
3198  <<sqrt(trajSOS.cartesianError().matrix()(2,2))<<" "
3199  <<sqrt(trajSOS.cartesianError().matrix()(3,3))<<" "
3200  <<sqrt(trajSOS.cartesianError().matrix()(4,4))<<" "
3201  <<sqrt(trajSOS.cartesianError().matrix()(5,5))<<std::endl;
3202  std::cout<<" q/p dxy/dz xy "<<trajSOS.localParameters().vector()(0)<<" "
3203  <<trajSOS.localParameters().vector()(1)<<" "
3204  <<trajSOS.localParameters().vector()(2)<<" "
3205  <<trajSOS.localParameters().vector()(3)<<" "
3206  <<trajSOS.localParameters().vector()(4)<<std::endl;
3207  std::cout<<" +/- error "<<sqrt(trajSOS.localError().matrix()(0,0))<<" "
3208  <<sqrt(trajSOS.localError().matrix()(1,1))<<" "
3209  <<sqrt(trajSOS.localError().matrix()(2,2))<<" "
3210  <<sqrt(trajSOS.localError().matrix()(3,3))<<" "
3211  <<sqrt(trajSOS.localError().matrix()(4,4))<<" "<<std::endl;
3212 }
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:71
TrajectoryStateOnSurface::localParameters
const LocalTrajectoryParameters & localParameters() const
Definition: TrajectoryStateOnSurface.h:167
TrajectoryStateOnSurface::cartesianError
const CartesianTrajectoryError cartesianError() const
Definition: TrajectoryStateOnSurface.h:161
TrajectoryStateOnSurface::globalPosition
GlobalPoint globalPosition() const
Definition: TrajectoryStateOnSurface.h:143
LocalTrajectoryParameters::vector
AlgebraicVector5 vector() const
Definition: LocalTrajectoryParameters.h:120
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:66
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
LocalTrajectoryError::matrix
const AlgebraicSymMatrix55 & matrix() const
Definition: LocalTrajectoryError.h:60
TrajectoryStateOnSurface::localError
const LocalTrajectoryError & localError() const
Definition: TrajectoryStateOnSurface.h:179
CartesianTrajectoryError::matrix
const AlgebraicSymMatrix66 & matrix() const
Definition: CartesianTrajectoryError.h:29
TrajectoryStateOnSurface::globalParameters
const GlobalTrajectoryParameters & globalParameters() const
Definition: TrajectoryStateOnSurface.h:140
TrajectoryStateOnSurface::globalMomentum
GlobalVector globalMomentum() const
Definition: TrajectoryStateOnSurface.h:146
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
GlobalTrajectoryParameters::vector
AlgebraicVector6 vector() const
Definition: GlobalTrajectoryParameters.h:100
void GlobalTrackerMuonAlignment::endJob ( void  )
privatevirtual

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:19
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
AlgebraicVector3
ROOT::Math::SVector< double, 3 > AlgebraicVector3
Definition: AlgebraicROOTObjects.h:11
j
int j
Definition: DBlmapReader.cc:9
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:81
dbtoconf.out
tuple out
Definition: dbtoconf.py:99
GlobalTrackerMuonAlignment::writeGlPosRcd
void writeGlPosRcd(CLHEP::HepVector &d3)
Definition: GlobalTrackerMuonAlignment.cc:3231
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 2095 of file GlobalTrackerMuonAlignment.cc.

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

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

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

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

References gather_cfg::cout, delta, i, 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().

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

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

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

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

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

3026 {
3027  bool info = false;
3028  bool smooth = false;
3029 
3030  if(info) std::cout<<" ......... start of muonFitter ........"<<std::endl;
3031  if(false && info){
3032  this->debugTrackHit(" Muon track hits ", alongTr);
3033  this->debugTrackHit(" Muon TransientTrack hits ", alongTTr);
3034  }
3035 
3036  edm::OwnVector<TrackingRecHit> recHit;
3037  DetId trackDetId = DetId(alongTr->innerDetId());
3038  for (trackingRecHit_iterator i=alongTr->recHitsBegin();i!=alongTr->recHitsEnd(); i++) {
3039  recHit.push_back((*i)->clone());
3040  }
3041  if(direction != alongMomentum)
3042  {
3043  edm::OwnVector<TrackingRecHit> recHitAlong = recHit;
3044  recHit.clear();
3045  for(unsigned int ihit = recHitAlong.size()-1; ihit+1>0; ihit--){
3046  recHit.push_back(recHitAlong[ihit]);
3047  }
3048  recHitAlong.clear();
3049  }
3050  if(info)
3051  std::cout<<" ... Own recHit.size() DetId==Muon "<<recHit.size()<<" "<<trackDetId.rawId()<<std::endl;
3052 
3053  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu =
3054  MuRHBuilder->build(alongTTr.recHitsBegin(), alongTTr.recHitsEnd());
3055  if(info)
3056  std::cout<<" ...along.... recHitMu.size() "<<recHitMu.size()<<std::endl;
3057  if(direction != alongMomentum){
3058  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMuAlong = recHitMu;
3059  recHitMu.clear();
3060  for(unsigned int ihit = recHitMuAlong.size()-1; ihit+1>0; ihit--){
3061  recHitMu.push_back(recHitMuAlong[ihit]);
3062  }
3063  recHitMuAlong.clear();
3064  }
3065  if(info)
3066  std::cout<<" ...opposite... recHitMu.size() "<<recHitMu.size()<<std::endl;
3067 
3068  TrajectoryStateOnSurface firstTSOS;
3069  if(direction == alongMomentum) firstTSOS = alongTTr.innermostMeasurementState();
3070  else firstTSOS = alongTTr.outermostMeasurementState();
3071 
3072  AlgebraicSymMatrix55 CovLoc;
3073  CovLoc(0,0) = 1. * firstTSOS.localError().matrix()(0,0);
3074  CovLoc(1,1) = 10. * firstTSOS.localError().matrix()(1,1);
3075  CovLoc(2,2) = 10. * firstTSOS.localError().matrix()(2,2);
3076  CovLoc(3,3) = 100. * firstTSOS.localError().matrix()(3,3);
3077  CovLoc(4,4) = 100. * firstTSOS.localError().matrix()(4,4);
3078  TrajectoryStateOnSurface initialTSOS(firstTSOS.localParameters(), LocalTrajectoryError(CovLoc),
3079  firstTSOS.surface(), &*magneticField_);
3080 
3081  PTrajectoryStateOnDet PTraj =
3082  trajectoryStateTransform::persistentState(initialTSOS, trackDetId.rawId());
3083  //const TrajectorySeed seedT(*PTraj, recHit, alongMomentum);
3084  const TrajectorySeed seedT(PTraj, recHit, direction);
3085 
3086  std::vector<Trajectory> trajVec;
3087  if(direction == alongMomentum) trajVec = theFitter->fit(seedT, recHitMu, initialTSOS);
3088  else trajVec = theFitterOp->fit(seedT, recHitMu, initialTSOS);
3089 
3090  if(info){
3091  this->debugTrajectorySOSv(" innermostTSOS of TransientTrack",
3092  alongTTr.innermostMeasurementState());
3093  this->debugTrajectorySOSv(" outermostTSOS of TransientTrack",
3094  alongTTr.outermostMeasurementState());
3095  this->debugTrajectorySOS(" initialTSOS for theFitter ", initialTSOS);
3096  std::cout<<" . . . . . trajVec.size() "<<trajVec.size()<<std::endl;
3097  if(trajVec.size()) this->debugTrajectory(" theFitter trajectory ", trajVec[0]);
3098  }
3099 
3100  if(!smooth){
3101  if(trajVec.size()) trackFittedTSOS = trajVec[0].lastMeasurement().updatedState();}
3102  else{
3103  std::vector<Trajectory> trajSVec;
3104  if (trajVec.size()){
3105  if(direction == alongMomentum) trajSVec = theSmoother->trajectories(*(trajVec.begin()));
3106  else trajSVec = theSmootherOp->trajectories(*(trajVec.begin()));
3107  }
3108  if(info) std::cout<<" . . . . trajSVec.size() "<<trajSVec.size()<<std::endl;
3109  if(trajSVec.size()) this->debugTrajectorySOSv("smoothed geom InnermostState",
3110  trajSVec[0].geometricalInnermostState());
3111  if(trajSVec.size()) trackFittedTSOS = trajSVec[0].geometricalInnermostState();
3112  }
3113 
3114 } // 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:24
GlobalTrackerMuonAlignment::debugTrajectorySOSv
void debugTrajectorySOSv(const std::string, TrajectoryStateOnSurface)
Definition: GlobalTrackerMuonAlignment.cc:3183
TrajectoryStateOnSurface::localParameters
const LocalTrajectoryParameters & localParameters() const
Definition: TrajectoryStateOnSurface.h:167
MuonTransientTrackingRecHitBuilder::ConstRecHitContainer
TransientTrackingRecHit::ConstRecHitContainer ConstRecHitContainer
Definition: MuonTransientTrackingRecHitBuilder.h:17
edm::OwnVector::size
size_type size() const
Definition: OwnVector.h:247
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:21
reco::TransientTrack::innermostMeasurementState
TrajectoryStateOnSurface innermostMeasurementState() const
Definition: TransientTrack.h:81
DetId::rawId
uint32_t rawId() const
get the raw id
Definition: DetId.h:45
KFTrajectorySmoother::trajectories
virtual std::vector< Trajectory > trajectories(const Trajectory &aTraj) const
Definition: KFTrajectorySmoother.cc:26
edm::OwnVector::push_back
void push_back(D *&d)
Definition: OwnVector.h:273
GlobalTrackerMuonAlignment::debugTrackHit
void debugTrackHit(const std::string, reco::TrackRef)
Definition: GlobalTrackerMuonAlignment.cc:3135
KFTrajectoryFitter::fit
virtual std::vector< Trajectory > fit(const Trajectory &aTraj) const
Definition: KFTrajectoryFitter.cc:20
edm::OwnVector::clear
void clear()
Definition: OwnVector.h:370
LocalTrajectoryError::matrix
const AlgebraicSymMatrix55 & matrix() const
Definition: LocalTrajectoryError.h:60
TrajectoryStateOnSurface::localError
const LocalTrajectoryError & localError() const
Definition: TrajectoryStateOnSurface.h:179
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:20
GlobalTrackerMuonAlignment::debugTrajectory
void debugTrajectory(const std::string, Trajectory &)
Definition: GlobalTrackerMuonAlignment.cc:3215
TrajectoryStateOnSurface::surface
const Surface & surface() const
Definition: TrajectoryStateOnSurface.h:182
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:3151
void GlobalTrackerMuonAlignment::trackFitter ( reco::TrackRef  alongTr,
reco::TransientTrack &  alongTTr,
PropagationDirection  direction,
TrajectoryStateOnSurface &  trackFittedTSOS 
)

Definition at line 2909 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 RSFinalFitAnalytical_cff::TTRHBuilder.

Referenced by analyzeTrackTrajectory().

2912 {
2913  bool info = false;
2914  bool smooth = false;
2915 
2916  if(info)
2917  std::cout<<" ......... start of trackFitter ......... "<<std::endl;
2918  if(false && info){
2919  this->debugTrackHit(" Tracker track hits ", alongTr);
2920  this->debugTrackHit(" Tracker TransientTrack hits ", alongTTr);
2921  }
2922 
2923  edm::OwnVector<TrackingRecHit> recHit;
2924  DetId trackDetId = DetId(alongTr->innerDetId());
2925  for (trackingRecHit_iterator i=alongTr->recHitsBegin();i!=alongTr->recHitsEnd(); i++) {
2926  recHit.push_back((*i)->clone());
2927  }
2928  if(direction != alongMomentum)
2929  {
2930  edm::OwnVector<TrackingRecHit> recHitAlong = recHit;
2931  recHit.clear();
2932  for(unsigned int ihit = recHitAlong.size()-1; ihit+1>0; ihit--){
2933  recHit.push_back(recHitAlong[ihit]);
2934  }
2935  recHitAlong.clear();
2936  }
2937  if(info)
2938  std::cout<<" ... Own recHit.size() "<<recHit.size()<<" "<<trackDetId.rawId()<<std::endl;
2939 
2940  //MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitTrack;
2941  TransientTrackingRecHit::RecHitContainer recHitTrack;
2942  for (trackingRecHit_iterator i=alongTr->recHitsBegin();i!=alongTr->recHitsEnd(); i++) {
2943  if((*i)->isValid()){
2944  recHitTrack.push_back(TTRHBuilder->build(&**i ));}
2945  else{
2946  recHitTrack.push_back(TTRHBuilder->build(&**i ));}
2947  }
2948 
2949  if(info)
2950  std::cout<<" ... recHitTrack.size() "<<recHit.size()<<" "<<trackDetId.rawId()
2951  <<" ======> recHitMu "<<std::endl;
2952 
2953  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu = recHitTrack;
2954  /*
2955  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu =
2956  MuRHBuilder->build(alongTTr.recHitsBegin(), alongTTr.recHitsEnd());
2957  */
2958  if(info)
2959  std::cout<<" ...along.... recHitMu.size() "<<recHitMu.size()<<std::endl;
2960  if(direction != alongMomentum){
2961  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMuAlong = recHitMu;
2962  recHitMu.clear();
2963  for(unsigned int ihit = recHitMuAlong.size()-1; ihit+1>0; ihit--){
2964  recHitMu.push_back(recHitMuAlong[ihit]);
2965  }
2966  recHitMuAlong.clear();
2967  }
2968  if(info)
2969  std::cout<<" ...opposite... recHitMu.size() "<<recHitMu.size()<<std::endl;
2970 
2971  TrajectoryStateOnSurface firstTSOS;
2972  if(direction == alongMomentum) firstTSOS = alongTTr.innermostMeasurementState();
2973  else firstTSOS = alongTTr.outermostMeasurementState();
2974 
2975  AlgebraicSymMatrix55 CovLoc;
2976  CovLoc(0,0) = 1. * firstTSOS.localError().matrix()(0,0);
2977  CovLoc(1,1) = 10. * firstTSOS.localError().matrix()(1,1);
2978  CovLoc(2,2) = 10. * firstTSOS.localError().matrix()(2,2);
2979  CovLoc(3,3) = 100. * firstTSOS.localError().matrix()(3,3);
2980  CovLoc(4,4) = 100. * firstTSOS.localError().matrix()(4,4);
2981  TrajectoryStateOnSurface initialTSOS(firstTSOS.localParameters(), LocalTrajectoryError(CovLoc),
2982  firstTSOS.surface(), &*magneticField_);
2983 
2984  PTrajectoryStateOnDet PTraj =
2985  trajectoryStateTransform::persistentState(initialTSOS, trackDetId.rawId());
2986  //const TrajectorySeed seedT(*PTraj, recHit, alongMomentum);
2987  const TrajectorySeed seedT(PTraj, recHit, direction);
2988 
2989  std::vector<Trajectory> trajVec;
2990  if(direction == alongMomentum) trajVec = theFitter->fit(seedT, recHitMu, initialTSOS);
2991  else trajVec = theFitterOp->fit(seedT, recHitMu, initialTSOS);
2992 
2993  if(info){
2994  this->debugTrajectorySOSv(" innermostTSOS of TransientTrack",
2995  alongTTr.innermostMeasurementState());
2996  this->debugTrajectorySOSv(" outermostTSOS of TransientTrack",
2997  alongTTr.outermostMeasurementState());
2998  this->debugTrajectorySOS(" initialTSOS for theFitter ", initialTSOS);
2999  std::cout<<" . . . . . trajVec.size() "<<trajVec.size()<<std::endl;
3000  if(trajVec.size()) this->debugTrajectory(" theFitter trajectory ", trajVec[0]);
3001  }
3002 
3003  if(!smooth){
3004  if(trajVec.size()) trackFittedTSOS = trajVec[0].lastMeasurement().updatedState();}
3005  else{
3006  std::vector<Trajectory> trajSVec;
3007  if (trajVec.size()){
3008  if(direction == alongMomentum) trajSVec = theSmoother->trajectories(*(trajVec.begin()));
3009  else trajSVec = theSmootherOp->trajectories(*(trajVec.begin()));
3010  }
3011  if(info) std::cout<<" . . . . trajSVec.size() "<<trajSVec.size()<<std::endl;
3012  if(trajSVec.size()) this->debugTrajectorySOSv("smoothed geom InnermostState",
3013  trajSVec[0].geometricalInnermostState());
3014  if(trajSVec.size()) trackFittedTSOS = trajSVec[0].geometricalInnermostState();
3015  }
3016 
3017  if(info) this->debugTrajectorySOSv(" trackFittedTSOS ", trackFittedTSOS);
3018 
3019 } // end of trackFitter
i
int i
Definition: DBlmapReader.cc:9
GlobalTrackerMuonAlignment::debugTrajectorySOSv
void debugTrajectorySOSv(const std::string, TrajectoryStateOnSurface)
Definition: GlobalTrackerMuonAlignment.cc:3183
TrajectoryStateOnSurface::localParameters
const LocalTrajectoryParameters & localParameters() const
Definition: TrajectoryStateOnSurface.h:167
MuonTransientTrackingRecHitBuilder::ConstRecHitContainer
TransientTrackingRecHit::ConstRecHitContainer ConstRecHitContainer
Definition: MuonTransientTrackingRecHitBuilder.h:17
edm::OwnVector::size
size_type size() const
Definition: OwnVector.h:247
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:21
TransientTrackingRecHit::RecHitContainer
std::vector< ConstRecHitPointer > RecHitContainer
Definition: TransientTrackingRecHit.h:29
reco::TransientTrack::innermostMeasurementState
TrajectoryStateOnSurface innermostMeasurementState() const
Definition: TransientTrack.h:81
DetId::rawId
uint32_t rawId() const
get the raw id
Definition: DetId.h:45
KFTrajectorySmoother::trajectories
virtual std::vector< Trajectory > trajectories(const Trajectory &aTraj) const
Definition: KFTrajectorySmoother.cc:26
edm::OwnVector::push_back
void push_back(D *&d)
Definition: OwnVector.h:273
GlobalTrackerMuonAlignment::debugTrackHit
void debugTrackHit(const std::string, reco::TrackRef)
Definition: GlobalTrackerMuonAlignment.cc:3135
KFTrajectoryFitter::fit
virtual std::vector< Trajectory > fit(const Trajectory &aTraj) const
Definition: KFTrajectoryFitter.cc:20
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:370
LocalTrajectoryError::matrix
const AlgebraicSymMatrix55 & matrix() const
Definition: LocalTrajectoryError.h:60
TrajectoryStateOnSurface::localError
const LocalTrajectoryError & localError() const
Definition: TrajectoryStateOnSurface.h:179
reco::TransientTrack::outermostMeasurementState
TrajectoryStateOnSurface outermostMeasurementState() const
Definition: TransientTrack.h:78
DetId
Definition: DetId.h:20
GlobalTrackerMuonAlignment::TTRHBuilder
const TransientTrackingRecHitBuilder * TTRHBuilder
Definition: GlobalTrackerMuonAlignment.cc:188
GlobalTrackerMuonAlignment::debugTrajectory
void debugTrajectory(const std::string, Trajectory &)
Definition: GlobalTrackerMuonAlignment.cc:3215
TrajectoryStateOnSurface::surface
const Surface & surface() const
Definition: TrajectoryStateOnSurface.h:182
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
GlobalTrackerMuonAlignment::debugTrajectorySOS
void debugTrajectorySOS(const std::string, TrajectoryStateOnSurface &)
Definition: GlobalTrackerMuonAlignment.cc:3151
void GlobalTrackerMuonAlignment::writeGlPosRcd ( CLHEP::HepVector &  d3)

Definition at line 3231 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, AlignTransform::rawId(), DetId::rawId(), AlignTransform::rotation(), indexGen::s2, funct::sin(), DetId::Tracker, patCandidatesForDimuonsSequences_cff::tracker, AlignTransform::translation(), and cond::service::PoolDBOutputService::writeOne().

Referenced by endJob().

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

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