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
 
char MuSelect [100]
 
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)
 
- 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:1249
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, trackerHits::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  Surface* surf = (Surface*)&refSurface;
829  const Plane* refPlane = dynamic_cast<Plane*>(surf);
830  GlobalVector Nlp = refPlane->normalVector();
831 
832  // extrapolation to innermost muon hit
833  extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
834  //extrapolationT = propagator->propagate(outerTrackTSOS, refSurface);
835 
836  if(!extrapolationT.isValid()){
837  if(false & alarm)
838  std::cout<<" !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
839  //<<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
840  <<std::endl;
841  continue;
842  }
843  if(debug_)
844  std::cout<<" extrapolationT.isValid "<<extrapolationT.isValid()<<std::endl;
845 
846  tsosMuonIf = 1;
847  Rt = GlobalVector((extrapolationT.globalPosition()).x(),
848  (extrapolationT.globalPosition()).y(),
849  (extrapolationT.globalPosition()).z());
850 
851  Pt = extrapolationT.globalMomentum();
852  // global parameters of muon
853  GRm = GlobalVector((innerMuTSOS.globalPosition()).x(),
854  (innerMuTSOS.globalPosition()).y(),
855  (innerMuTSOS.globalPosition()).z());
856  GPm = innerMuTSOS.globalMomentum();
857  Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
858  (outerTrackTSOS.globalPosition()).y(),
859  (outerTrackTSOS.globalPosition()).z());
860  Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() +
861  innerMuTSOS.cartesianError().matrix());
862  C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
863  Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
864  C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
865 
866  if(false & debug_){
867  //std::cout<<" ->propDir "<<propagator->propagationDirection()<<std::endl;
868  PropagationDirectionChooser Chooser;
869  std::cout<<" propDirCh "
870  <<Chooser.operator()(*outerTrackTSOS.freeState(), refSurface)
871  <<" Ch == along "<<(alongMomentum ==
872  Chooser.operator()(*outerTrackTSOS.freeState(), refSurface))
873  <<std::endl;
874  std::cout<<" --- Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
875  <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
876  }
877  } // enf of ---- Out - In -----
878 
879 
880  else if((distanceInOut <= distanceInIn) & (distanceInOut <= distanceOutIn) &
881  (distanceInOut <= distanceOutOut)){ // ----- In - Out ------
882  if(debug_) std::cout<<" ----- In - Out -----"<<std::endl;
883 
884  const Surface& refSurface = outerMuTSOS.surface();
885  ReferenceCountingPointer<TangentPlane>
886  tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
887  Nl = tpMuGlobal->normalVector();
888 
889  // extrapolation to outermost muon hit
890  extrapolationT = oppositeSmPr.propagate(innerTrackTSOS, refSurface);
891 
892  if(!extrapolationT.isValid()){
893  if(false & alarm)
894  std::cout<<" !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
895  <<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
896  <<std::endl;
897  continue;
898  }
899  if(debug_)
900  std::cout<<" extrapolationT.isValid "<<extrapolationT.isValid()<<std::endl;
901 
902  tsosMuonIf = 2;
903  Rt = GlobalVector((extrapolationT.globalPosition()).x(),
904  (extrapolationT.globalPosition()).y(),
905  (extrapolationT.globalPosition()).z());
906 
907  Pt = extrapolationT.globalMomentum();
908  // global parameters of muon
909  GRm = GlobalVector((outerMuTSOS.globalPosition()).x(),
910  (outerMuTSOS.globalPosition()).y(),
911  (outerMuTSOS.globalPosition()).z());
912  GPm = outerMuTSOS.globalMomentum();
913  Rt0 = GlobalVector((innerTrackTSOS.globalPosition()).x(),
914  (innerTrackTSOS.globalPosition()).y(),
915  (innerTrackTSOS.globalPosition()).z());
916  Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() +
917  outerMuTSOS.cartesianError().matrix());
918  C0 = AlgebraicSymMatrix66(innerTrackTSOS.cartesianError().matrix());
919  Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
920  C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
921 
922  if(false & debug_){
923  //std::cout<<" ->propDir "<<propagator->propagationDirection()<<std::endl;
924  PropagationDirectionChooser Chooser;
925  std::cout<<" propDirCh "
926  <<Chooser.operator()(*innerTrackTSOS.freeState(), refSurface)
927  <<" Ch == oppisite "<<(oppositeToMomentum ==
928  Chooser.operator()(*innerTrackTSOS.freeState(), refSurface))
929  <<std::endl;
930  std::cout<<" --- Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
931  <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
932  }
933  } // enf of ---- In - Out -----
934 
935  else if((distanceOutOut <= distanceInOut) & (distanceOutOut <= distanceInIn) &
936  (distanceOutOut <= distanceOutIn)){ // ----- Out - Out ------
937  if(debug_) std::cout<<" ----- Out - Out -----"<<std::endl;
938 
939  // reject: momentum of track has opposite direction to muon track
940  continue;
941 
942  const Surface& refSurface = outerMuTSOS.surface();
943  ReferenceCountingPointer<TangentPlane>
944  tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
945  Nl = tpMuGlobal->normalVector();
946 
947  // extrapolation to outermost muon hit
948  extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
949 
950  if(!extrapolationT.isValid()){
951  if(alarm)
952  std::cout<<" !!!!!!!!! Catastrophe Out-Out extrapolationT.isValid "
953  <<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
954  <<std::endl;
955  continue;
956  }
957  if(debug_)
958  std::cout<<" extrapolationT.isValid "<<extrapolationT.isValid()<<std::endl;
959 
960  tsosMuonIf = 2;
961  Rt = GlobalVector((extrapolationT.globalPosition()).x(),
962  (extrapolationT.globalPosition()).y(),
963  (extrapolationT.globalPosition()).z());
964 
965  Pt = extrapolationT.globalMomentum();
966  // global parameters of muon
967  GRm = GlobalVector((outerMuTSOS.globalPosition()).x(),
968  (outerMuTSOS.globalPosition()).y(),
969  (outerMuTSOS.globalPosition()).z());
970  GPm = outerMuTSOS.globalMomentum();
971  Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
972  (outerTrackTSOS.globalPosition()).y(),
973  (outerTrackTSOS.globalPosition()).z());
974  Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() +
975  outerMuTSOS.cartesianError().matrix());
976  C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
977  Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
978  C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
979 
980  if(debug_){
981  PropagationDirectionChooser Chooser;
982  std::cout<<" propDirCh "
983  <<Chooser.operator()(*outerTrackTSOS.freeState(), refSurface)
984  <<" Ch == along "<<(alongMomentum ==
985  Chooser.operator()(*outerTrackTSOS.freeState(), refSurface))
986  <<std::endl;
987  std::cout<<" --- Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
988  <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
989  std::cout<<" Nornal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
990  <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
991  }
992  } // enf of ---- Out - Out -----
993  else{
994  if(alarm)
995  std::cout<<" ------- !!!!!!!!!!!!!!! No proper Out - In \a\a\a\a\a\a\a"<<std::endl;
996  continue;
997  }
998 
999  } // ------------------------------- end Cosmic Muon -----
1000 
1001  if(tsosMuonIf == 0) {if(info) {std::cout<<"No tsosMuon !!!!!!"<<std::endl; continue;}}
1002  TrajectoryStateOnSurface tsosMuon;
1003  if(tsosMuonIf == 1) tsosMuon = muTT.innermostMeasurementState();
1004  else tsosMuon = muTT.outermostMeasurementState();
1005 
1006  //GlobalTrackerMuonAlignment::misalignMuon(GRm, GPm, Nl, Rt, Rm, Pm);
1007  AlgebraicVector4 LPRm; // muon local (dx/dz, dy/dz, x, y)
1008  GlobalTrackerMuonAlignment::misalignMuonL
1009  (GRm, GPm, Nl, Rt, Rm, Pm, LPRm, extrapolationT, tsosMuon);
1010 
1011  GlobalVector resR = Rm - Rt;
1012  GlobalVector resP0 = Pm - Pt;
1013  GlobalVector resP = Pm / Pm.mag() - Pt / Pt.mag();
1014  float RelMomResidual = (Pm.mag() - Pt.mag()) / (Pt.mag() + 1.e-6);;
1015 
1016  AlgebraicVector6 Vm;
1017  Vm(0) = resR.x(); Vm(1) = resR.y(); Vm(2) = resR.z();
1018  Vm(3) = resP0.x(); Vm(4) = resP0.y(); Vm(5) = resP0.z();
1019  float Rmuon = Rm.perp();
1020  float Zmuon = Rm.z();
1021  float alfa_x = atan2(Nl.x(),Nl.y())*57.29578;
1022 
1023  if(debug_){
1024  std::cout<<" Nx Ny Nz alfa_x "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "<<alfa_x<<std::endl;
1025  //std::cout<<" Rm "<<Rm<<std::endl<<" Rt "<<Rt<<std::endl<<" resR "<<resR<<std::endl
1026  // <<" resP "<<resP<<" dp/p "<<RelMomResidual<<std::endl;
1027  }
1028 
1029  double chi_d = 0;
1030  for(int i=0; i<=5; i++) chi_d += Vm(i)*Vm(i)/Cm(i,i);
1031 
1032  AlgebraicVector5 Vml(tsosMuon.localParameters().vector()
1033  - extrapolationT.localParameters().vector());
1034  AlgebraicSymMatrix55 m(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
1035  AlgebraicSymMatrix55 Cml(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
1036  bool ierrLoc = !m.Invert();
1037  if (ierrLoc && debug_ && info) {
1038  std::cout<< " ==== Error inverting Local covariance matrix ==== "<<std::cout;
1039  continue;}
1040  double chi_Loc = ROOT::Math::Similarity(Vml,m);
1041  if(debug_)
1042  std::cout<<" chi_Loc px/pz/err "<<chi_Loc<<" "<<Vml(1)/sqrt(Cml(1,1))<<std::endl;
1043 
1044  if(Pt.mag() < 15.) continue;
1045  if(Pm.mag() < 5.) continue;
1046 
1047  //if(Pt.mag() < 30.) continue; // momenum cut < 30GeV
1048  //if(Pt.mag() < 60.) continue; // momenum cut < 30GeV
1049  //if(Pt.mag() > 50.) continue; // momenum cut > 50GeV
1050  //if(Pt.mag() > 100.) continue; // momenum cut > 100GeV
1051  //if(trackTT.charge() < 0) continue; // select positive charge
1052  //if(trackTT.charge() > 0) continue; // select negative charge
1053 
1054  //if(fabs(resR.x()) > 5.) continue; // strong cut X
1055  //if(fabs(resR.y()) > 5.) continue; // Y
1056  //if(fabs(resR.z()) > 5.) continue; // Z
1057  //if(fabs(resR.mag()) > 7.5) continue; // dR
1058 
1059  /*
1060  //if(fabs(RelMomResidual) > 0.5) continue;
1061  if(fabs(resR.x()) > 20.) continue;
1062  if(fabs(resR.y()) > 20.) continue;
1063  if(fabs(resR.z()) > 20.) continue;
1064  if(fabs(resR.mag()) > 30.) continue;
1065  if(fabs(resP.x()) > 0.06) continue;
1066  if(fabs(resP.y()) > 0.06) continue;
1067  if(fabs(resP.z()) > 0.06) continue;
1068  if(chi_d > 40.) continue;
1069  */
1070 
1071  // select Barrel
1072  //if(Rmuon < 400. || Rmuon > 450.) continue;
1073  //if(Zmuon < -600. || Zmuon > 600.) continue;
1074  //if(fabs(Nl.z()) > 0.95) continue;
1075  //std::sprintf(MuSelect, " Barrel");
1076  // EndCap1
1077  //if(Rmuon < 120. || Rmuon > 450.) continue;
1078  //if(Zmuon < -720.) continue;
1079  //if(Zmuon > -580.) continue;
1080  //if(fabs(Nl.z()) < 0.95) continue;
1081  //std::sprintf(MuSelect, " EndCap1");
1082  // EndCap2
1083  //if(Rmuon < 120. || Rmuon > 450.) continue;
1084  //if(Zmuon > 720.) continue;
1085  //if(Zmuon < 580.) continue;
1086  //if(fabs(Nl.z()) < 0.95) continue;
1087  //std::sprintf(MuSelect, " EndCap2");
1088  // select All
1089  if(Rmuon < 120. || Rmuon > 450.) continue;
1090  if(Zmuon < -720. || Zmuon > 720.) continue;
1091  std::sprintf(MuSelect, " Barrel+EndCaps");
1092 
1093 
1094 
1095  if(debug_)
1096  std::cout<<" .............. passed all cuts"<<std::endl;
1097 
1098  N_track++;
1099  // gradient and Hessian for each track
1100 
1101  GlobalTrackerMuonAlignment::gradientGlobalAlg(Rt, Pt, Rm, Nl, Cm);
1102  GlobalTrackerMuonAlignment::gradientGlobal(Rt, Pt, Rm, Pm, Nl, Cm);
1103 
1104  CLHEP::HepSymMatrix covLoc(4,0);
1105  for(int i=1; i<=4; i++)
1106  for(int j=1; j<=i; j++){
1107  covLoc(i,j) = (tsosMuon.localError().matrix()
1108  + extrapolationT.localError().matrix())(i,j);
1109  //if(i != j) Cov(i,j) = 0.;
1110  }
1111 
1112  const Surface& refSurface = tsosMuon.surface();
1113  CLHEP::HepMatrix rotLoc (3,3,0);
1114  rotLoc(1,1) = refSurface.rotation().xx();
1115  rotLoc(1,2) = refSurface.rotation().xy();
1116  rotLoc(1,3) = refSurface.rotation().xz();
1117 
1118  rotLoc(2,1) = refSurface.rotation().yx();
1119  rotLoc(2,2) = refSurface.rotation().yy();
1120  rotLoc(2,3) = refSurface.rotation().yz();
1121 
1122  rotLoc(3,1) = refSurface.rotation().zx();
1123  rotLoc(3,2) = refSurface.rotation().zy();
1124  rotLoc(3,3) = refSurface.rotation().zz();
1125 
1126  CLHEP::HepVector posLoc(3);
1127  posLoc(1) = refSurface.position().x();
1128  posLoc(2) = refSurface.position().y();
1129  posLoc(3) = refSurface.position().z();
1130 
1131  GlobalTrackerMuonAlignment::gradientLocal(Rt, Pt, Rm, Pm, Nl, covLoc, rotLoc, posLoc, LPRm);
1132 
1133  if(debug_){
1134  std::cout<<" Norm "<<Nl<<std::endl;
1135  std::cout<<" posLoc "<<posLoc.T()<<std::endl;
1136  std::cout<<" rotLoc "<<rotLoc<<std::endl;
1137  }
1138 
1139  // ----------------------------------------------------- fill histogram
1140  histo->Fill(itMuon->track()->pt());
1141 
1142  //histo2->Fill(itMuon->track()->outerP());
1143  histo2->Fill(Pt.mag());
1144  histo3->Fill((PI/2.-itMuon->track()->outerTheta()));
1145  histo4->Fill(itMuon->track()->phi());
1146  histo5->Fill(Rmuon);
1147  histo6->Fill(Zmuon);
1148  histo7->Fill(RelMomResidual);
1149  //histo8->Fill(chi);
1150  histo8->Fill(chi_d);
1151 
1152  histo101->Fill(Zmuon, Rmuon);
1153  histo101->Fill(Rt0.z(), Rt0.perp());
1154  histo102->Fill(Rt0.x(), Rt0.y());
1155  histo102->Fill(Rm.x(), Rm.y());
1156 
1157  histo11->Fill(resR.mag());
1158  if(fabs(Nl.x()) < 0.98) histo12->Fill(resR.x());
1159  if(fabs(Nl.y()) < 0.98) histo13->Fill(resR.y());
1160  if(fabs(Nl.z()) < 0.98) histo14->Fill(resR.z());
1161  histo15->Fill(resP.x());
1162  histo16->Fill(resP.y());
1163  histo17->Fill(resP.z());
1164 
1165  if((fabs(Nl.x()) < 0.98) && (fabs(Nl.y()) < 0.98) &&(fabs(Nl.z()) < 0.98))
1166  {
1167  histo18->Fill(sqrt(C0(0,0)));
1168  histo19->Fill(sqrt(C1(0,0)));
1169  histo20->Fill(sqrt(C1(0,0)+Ce(0,0)));
1170  }
1171  if(fabs(Nl.x()) < 0.98) histo21->Fill(Vm(0)/sqrt(Cm(0,0)));
1172  if(fabs(Nl.y()) < 0.98) histo22->Fill(Vm(1)/sqrt(Cm(1,1)));
1173  if(fabs(Nl.z()) < 0.98) histo23->Fill(Vm(2)/sqrt(Cm(2,2)));
1174  histo24->Fill(Vm(3)/sqrt(C1(3,3)+Ce(3,3)));
1175  histo25->Fill(Vm(4)/sqrt(C1(4,4)+Ce(4,4)));
1176  histo26->Fill(Vm(5)/sqrt(C1(5,5)+Ce(5,5)));
1177  histo27->Fill(Nl.x());
1178  histo28->Fill(Nl.y());
1179  histo29->Fill(lenghtTrack);
1180  histo30->Fill(lenghtMuon);
1181  histo31->Fill(chi_Loc);
1182  histo32->Fill(Vml(1)/sqrt(Cml(1,1)));
1183  histo33->Fill(Vml(2)/sqrt(Cml(2,2)));
1184  histo34->Fill(Vml(3)/sqrt(Cml(3,3)));
1185  histo35->Fill(Vml(4)/sqrt(Cml(4,4)));
1186 
1187  if (debug_) { //--------------------------------- debug print ----------
1188 
1189  std::cout<<" diag 0[ "<<C0(0,0)<<" "<<C0(1,1)<<" "<<C0(2,2)<<" "<<C0(3,3)<<" "
1190  <<C0(4,4)<<" "<<C0(5,5)<<" ]"<<std::endl;
1191  std::cout<<" diag e[ "<<Ce(0,0)<<" "<<Ce(1,1)<<" "<<Ce(2,2)<<" "<<Ce(3,3)<<" "
1192  <<Ce(4,4)<<" "<<Ce(5,5)<<" ]"<<std::endl;
1193  std::cout<<" diag 1[ "<<C1(0,0)<<" "<<C1(1,1)<<" "<<C1(2,2)<<" "<<C1(3,3)<<" "
1194  <<C1(4,4)<<" "<<C1(5,5)<<" ]"<<std::endl;
1195  std::cout<<" Rm "<<Rm.x()<<" "<<Rm.y()<<" "<<Rm.z()
1196  <<" Pm "<<Pm.x()<<" "<<Pm.y()<<" "<<Pm.z()<<std::endl;
1197  std::cout<<" Rt "<<Rt.x()<<" "<<Rt.y()<<" "<<Rt.z()
1198  <<" Pt "<<Pt.x()<<" "<<Pt.y()<<" "<<Pt.z()<<std::endl;
1199  std::cout<<" Nl*(Rm-Rt) "<<Nl.dot(Rm-Rt)<<std::endl;
1200  std::cout<<" resR "<<resR.x()<<" "<<resR.y()<<" "<<resR.z()
1201  <<" resP "<<resP.x()<<" "<<resP.y()<<" "<<resP.z()<<std::endl;
1202  std::cout<<" Rm-t "<<(Rm-Rt).x()<<" "<<(Rm-Rt).y()<<" "<<(Rm-Rt).z()
1203  <<" Pm-t "<<(Pm-Pt).x()<<" "<<(Pm-Pt).y()<<" "<<(Pm-Pt).z()<<std::endl;
1204  std::cout<<" Vm "<<Vm<<std::endl;
1205  std::cout<<" +- "<<sqrt(Cm(0,0))<<" "<<sqrt(Cm(1,1))<<" "<<sqrt(Cm(2,2))<<" "
1206  <<sqrt(Cm(3,3))<<" "<<sqrt(Cm(4,4))<<" "<<sqrt(Cm(5,5))<<std::endl;
1207  std::cout<<" Pmuon Ptrack dP/Ptrack "<<itMuon->outerTrack()->p()<<" "
1208  <<itMuon->track()->outerP()<<" "
1209  <<(itMuon->outerTrack()->p() -
1210  itMuon->track()->outerP())/itMuon->track()->outerP()<<std::endl;
1211  std::cout<<" cov matrix "<<std::endl;
1212  std::cout<<Cm<<std::endl;
1213  std::cout<<" diag [ "<<Cm(0,0)<<" "<<Cm(1,1)<<" "<<Cm(2,2)<<" "<<Cm(3,3)<<" "
1214  <<Cm(4,4)<<" "<<Cm(5,5)<<" ]"<<std::endl;
1215 
1216  static AlgebraicSymMatrix66 Ro;
1217  double Diag[6];
1218  for(int i=0; i<=5; i++) Diag[i] = sqrt(Cm(i,i));
1219  for(int i=0; i<=5; i++)
1220  for(int j=0; j<=5; j++)
1221  Ro(i,j) = Cm(i,j)/Diag[i]/Diag[j];
1222  std::cout<<" correlation matrix "<<std::endl;
1223  std::cout<<Ro<<std::endl;
1224 
1225  AlgebraicSymMatrix66 CmI;
1226  for(int i=0; i<=5; i++)
1227  for(int j=0; j<=5; j++)
1228  CmI(i,j) = Cm(i,j);
1229 
1230  bool ierr = !CmI.Invert();
1231  if( ierr ) { if(alarm || debug_)
1232  std::cout<<" Error inverse covariance matrix !!!!!!!!!!!"<<std::endl;
1233  continue;
1234  }
1235  std::cout<<" inverse cov matrix "<<std::endl;
1236  std::cout<<Cm<<std::endl;
1237 
1238  double chi = ROOT::Math::Similarity(Vm, CmI);
1239  std::cout<<" chi chi_d "<<chi<<" "<<chi_d<<std::endl;
1240  } // end of debug_ printout --------------------------------------------
1241 
1242  } // end loop on selected muons, i.e. Jim's globalMuon
1243 
1244 } //end of analyzeTrackTrack
TkRotation::xx
T xx() const
Definition: newTkRotation.h:217
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:2025
GlobalTrackerMuonAlignment::gradientLocal
void gradientLocal(GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, CLHEP::HepSymMatrix &, CLHEP::HepMatrix &, CLHEP::HepVector &, AlgebraicVector4 &)
Definition: GlobalTrackerMuonAlignment.cc:2277
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:66
TrajectoryStateOnSurface::localParameters
const LocalTrajectoryParameters & localParameters() const
Definition: TrajectoryStateOnSurface.h:152
PI
#define PI
Definition: HcalBeamMonitor.cc:14
mag
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
Definition: Basic3DVectorLD.h:138
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
PV3DBase::y
T y() const
Definition: PV3DBase.h:57
TkRotation::yx
T yx() const
Definition: newTkRotation.h:220
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:128
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:2104
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:101
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:223
TkRotation::xy
T xy() const
Definition: newTkRotation.h:218
LocalTrajectoryParameters::vector
AlgebraicVector5 vector() const
Definition: LocalTrajectoryParameters.h:147
TkRotation::zz
T zz() const
Definition: newTkRotation.h:225
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:61
SteppingHelixPropagator_cfi.SteppingHelixPropagator
tuple SteppingHelixPropagator
Definition: SteppingHelixPropagator_cfi.py:3
TrajectoryStateOnSurface::freeState
FreeTrajectoryState * freeState(bool withErrors=true) const
Definition: TrajectoryStateOnSurface.h:115
TrajectoryStateOnSurface::cartesianError
const CartesianTrajectoryError & cartesianError() const
Definition: TrajectoryStateOnSurface.h:146
GlobalTrackerMuonAlignment::watchTrackingGeometry_
edm::ESWatcher< GlobalTrackingGeometryRecord > watchTrackingGeometry_
Definition: GlobalTrackerMuonAlignment.cc:170
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:28
PV3DBase::z
T z() const
Definition: PV3DBase.h:58
edm::Event::run
RunNumber_t run() const
Definition: Event.h:66
GlobalTrackerMuonAlignment::iteratorMuonRcd
std::vector< AlignTransform >::const_iterator iteratorMuonRcd
Definition: GlobalTrackerMuonAlignment.cc:206
TkRotation::zy
T zy() const
Definition: newTkRotation.h:224
j
int j
Definition: DBlmapReader.cc:9
LocalTrajectoryError::matrix
const AlgebraicSymMatrix55 & matrix() const
Definition: LocalTrajectoryError.h:50
TrajectoryStateOnSurface::localError
const LocalTrajectoryError & localError() const
Definition: TrajectoryStateOnSurface.h:164
TkRotation::yy
T yy() const
Definition: newTkRotation.h:221
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:355
CartesianTrajectoryError::matrix
const AlgebraicSymMatrix66 & matrix() const
Definition: CartesianTrajectoryError.h:40
GlobalTrackerMuonAlignment::misalignMuonL
void misalignMuonL(GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, GlobalVector &, AlgebraicVector4 &, TrajectoryStateOnSurface &, TrajectoryStateOnSurface &)
Definition: GlobalTrackerMuonAlignment.cc:2654
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:131
patZpeak.muons
tuple muons
Definition: patZpeak.py:38
TkRotation::xz
T xz() const
Definition: newTkRotation.h:219
TrajectoryStateOnSurface::surface
const Surface & surface() const
Definition: TrajectoryStateOnSurface.h:167
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:41
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:56
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:222
GlobalVector
Global3DVector GlobalVector
Definition: GlobalVector.h:10
void GlobalTrackerMuonAlignment::analyzeTrackTrajectory ( const edm::Event &  iEvent,
const edm::EventSetup &  iSetup 
)

Definition at line 1249 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, trackerHits::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().

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

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

3145 {
3146  std::cout<<" ------- "<<title<<" --------"<<std::endl;
3147  int nHit = 1;
3148  TransientTrackingRecHit::RecHitContainer recHit;
3149  for (trackingRecHit_iterator i=alongTr->recHitsBegin();i!=alongTr->recHitsEnd(); i++) {
3150  std::cout<<" Hit "<<nHit++<<" DetId "<<(*i)->geographicalId().det();
3151  if( (*i)->geographicalId().det() == DetId::Tracker ) std::cout<<" Tracker ";
3152  else if ( (*i)->geographicalId().det() == DetId::Muon ) std::cout<<" Muon ";
3153  else std::cout<<" Unknown ";
3154  if(!(*i)->isValid()) std::cout<<" not valid "<<std::endl;
3155  else std::cout<<std::endl;
3156  }
3157 }
i
int i
Definition: DBlmapReader.cc:9
TransientTrackingRecHit::RecHitContainer
std::vector< ConstRecHitPointer > RecHitContainer
Definition: TransientTrackingRecHit.h:22
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:41
void GlobalTrackerMuonAlignment::debugTrackHit ( const std::string  title,
reco::TransientTrack &  alongTr 
)

Definition at line 3127 of file GlobalTrackerMuonAlignment.cc.

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

3128 {
3129  std::cout<<" ------- "<<title<<" --------"<<std::endl;
3130  int nHit = 1;
3131  TransientTrackingRecHit::RecHitContainer recHit;
3132  for (trackingRecHit_iterator i=alongTr.recHitsBegin();i!=alongTr.recHitsEnd(); i++) {
3133  std::cout<<" Hit "<<nHit++<<" DetId "<<(*i)->geographicalId().det();
3134  if( (*i)->geographicalId().det() == DetId::Tracker ) std::cout<<" Tracker ";
3135  else if ( (*i)->geographicalId().det() == DetId::Muon ) std::cout<<" Muon ";
3136  else std::cout<<" Unknown ";
3137  if(!(*i)->isValid()) std::cout<<" not valid "<<std::endl;
3138  else std::cout<<std::endl;
3139  }
3140 }
i
int i
Definition: DBlmapReader.cc:9
TransientTrackingRecHit::RecHitContainer
std::vector< ConstRecHitPointer > RecHitContainer
Definition: TransientTrackingRecHit.h:22
reco::TransientTrack::recHitsEnd
trackingRecHit_iterator recHitsEnd() const
last iterator to RecHits
Definition: TransientTrack.h:96
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:41
reco::TransientTrack::recHitsBegin
trackingRecHit_iterator recHitsBegin() const
first iterator to RecHits
Definition: TransientTrack.h:94
void GlobalTrackerMuonAlignment::debugTrajectory ( const std::string  title,
Trajectory &  traj 
)

Definition at line 3224 of file GlobalTrackerMuonAlignment.cc.

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

3225 {
3226  std::cout<<"\n"<<" ...... "<<title<<" ...... "<<std::endl;
3227  if(!traj.isValid()) {std::cout<<" Not valid !!!!!!!! "<<std::endl; return;}
3228  std::cout<<" chi2/Nhit "<<traj.chiSquared()<<" / "<<traj.foundHits();
3229  if(traj.direction() == alongMomentum) std::cout<<" alongMomentum >>>>"<<std::endl;
3230  else std::cout<<" oppositeToMomentum <<<<"<<std::endl;
3231  this->debugTrajectorySOSv(" firstMeasurementTSOS ",traj.firstMeasurement().updatedState());
3232  //this->debugTrajectorySOSv(" firstMeasurementTSOS ",traj.firstMeasurement().predictedState());
3233  this->debugTrajectorySOSv(" lastMeasurementTSOS ",traj.lastMeasurement().updatedState());
3234  //this->debugTrajectorySOSv(" geom InnermostState", traj.geometricalInnermostState());
3235  std::cout<<" . . . . . . . . . . . . . . . . . . . . . . . . . . . . \n"<<std::endl;
3236 }
Trajectory::foundHits
int foundHits() const
Definition: Trajectory.h:190
GlobalTrackerMuonAlignment::debugTrajectorySOSv
void debugTrajectorySOSv(const std::string, TrajectoryStateOnSurface)
Definition: GlobalTrackerMuonAlignment.cc:3192
Trajectory::direction
PropagationDirection const & direction() const
Definition: Trajectory.cc:195
Trajectory::lastMeasurement
TrajectoryMeasurement const & lastMeasurement() const
Definition: Trajectory.h:147
TrajectoryMeasurement::updatedState
TrajectoryStateOnSurface updatedState() const
Definition: TrajectoryMeasurement.h:162
Trajectory::isValid
bool isValid() const
Definition: Trajectory.h:225
Trajectory::firstMeasurement
TrajectoryMeasurement const & firstMeasurement() const
Definition: Trajectory.h:160
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:41
Trajectory::chiSquared
double chiSquared() const
Definition: Trajectory.h:208
void GlobalTrackerMuonAlignment::debugTrajectorySOS ( const std::string  title,
TrajectoryStateOnSurface &  trajSOS 
)

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

3162 {
3163  std::cout<<" --- "<<title<<" --- "<<std::endl;
3164  if(!trajSOS.isValid()) {std::cout<<" Not valid !!!! "<<std::endl; return;}
3165  std::cout<<" R |p| "<<trajSOS.globalPosition().perp()<<" "
3166  <<trajSOS.globalMomentum().mag()<<" charge "<<trajSOS.charge()<<std::endl;
3167  std::cout<<" x p "<<trajSOS.globalParameters().vector()(0)<<" "
3168  <<trajSOS.globalParameters().vector()(1)<<" "
3169  <<trajSOS.globalParameters().vector()(2)<<" "
3170  <<trajSOS.globalParameters().vector()(3)<<" "
3171  <<trajSOS.globalParameters().vector()(4)<<" "
3172  <<trajSOS.globalParameters().vector()(5)<<std::endl;
3173  std::cout<<" +/- "<<sqrt(trajSOS.cartesianError().matrix()(0,0))<<" "
3174  <<sqrt(trajSOS.cartesianError().matrix()(1,1))<<" "
3175  <<sqrt(trajSOS.cartesianError().matrix()(2,2))<<" "
3176  <<sqrt(trajSOS.cartesianError().matrix()(3,3))<<" "
3177  <<sqrt(trajSOS.cartesianError().matrix()(4,4))<<" "
3178  <<sqrt(trajSOS.cartesianError().matrix()(5,5))<<std::endl;
3179  std::cout<<" q/p dxy/dz xy "<<trajSOS.localParameters().vector()(0)<<" "
3180  <<trajSOS.localParameters().vector()(1)<<" "
3181  <<trajSOS.localParameters().vector()(2)<<" "
3182  <<trajSOS.localParameters().vector()(3)<<" "
3183  <<trajSOS.localParameters().vector()(4)<<std::endl;
3184  std::cout<<" +/- error "<<sqrt(trajSOS.localError().matrix()(0,0))<<" "
3185  <<sqrt(trajSOS.localError().matrix()(1,1))<<" "
3186  <<sqrt(trajSOS.localError().matrix()(2,2))<<" "
3187  <<sqrt(trajSOS.localError().matrix()(3,3))<<" "
3188  <<sqrt(trajSOS.localError().matrix()(4,4))<<" "<<std::endl;
3189 }
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:66
TrajectoryStateOnSurface::localParameters
const LocalTrajectoryParameters & localParameters() const
Definition: TrajectoryStateOnSurface.h:152
TrajectoryStateOnSurface::globalPosition
GlobalPoint globalPosition() const
Definition: TrajectoryStateOnSurface.h:128
LocalTrajectoryParameters::vector
AlgebraicVector5 vector() const
Definition: LocalTrajectoryParameters.h:147
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:61
TrajectoryStateOnSurface::cartesianError
const CartesianTrajectoryError & cartesianError() const
Definition: TrajectoryStateOnSurface.h:146
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:28
LocalTrajectoryError::matrix
const AlgebraicSymMatrix55 & matrix() const
Definition: LocalTrajectoryError.h:50
TrajectoryStateOnSurface::localError
const LocalTrajectoryError & localError() const
Definition: TrajectoryStateOnSurface.h:164
CartesianTrajectoryError::matrix
const AlgebraicSymMatrix66 & matrix() const
Definition: CartesianTrajectoryError.h:40
TrajectoryStateOnSurface::globalParameters
const GlobalTrajectoryParameters & globalParameters() const
Definition: TrajectoryStateOnSurface.h:125
TrajectoryStateOnSurface::globalMomentum
GlobalVector globalMomentum() const
Definition: TrajectoryStateOnSurface.h:131
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:41
GlobalTrajectoryParameters::vector
AlgebraicVector6 vector() const
Definition: GlobalTrajectoryParameters.h:90
void GlobalTrackerMuonAlignment::debugTrajectorySOSv ( const std::string  title,
TrajectoryStateOnSurface  trajSOS 
)

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

3194 {
3195  std::cout<<" --- "<<title<<" --- "<<std::endl;
3196  if(!trajSOS.isValid()) {std::cout<<" Not valid !!!! "<<std::endl; return;}
3197  std::cout<<" R |p| "<<trajSOS.globalPosition().perp()<<" "
3198  <<trajSOS.globalMomentum().mag()<<" charge "<<trajSOS.charge()<<std::endl;
3199  std::cout<<" x p "<<trajSOS.globalParameters().vector()(0)<<" "
3200  <<trajSOS.globalParameters().vector()(1)<<" "
3201  <<trajSOS.globalParameters().vector()(2)<<" "
3202  <<trajSOS.globalParameters().vector()(3)<<" "
3203  <<trajSOS.globalParameters().vector()(4)<<" "
3204  <<trajSOS.globalParameters().vector()(5)<<std::endl;
3205  std::cout<<" +/- "<<sqrt(trajSOS.cartesianError().matrix()(0,0))<<" "
3206  <<sqrt(trajSOS.cartesianError().matrix()(1,1))<<" "
3207  <<sqrt(trajSOS.cartesianError().matrix()(2,2))<<" "
3208  <<sqrt(trajSOS.cartesianError().matrix()(3,3))<<" "
3209  <<sqrt(trajSOS.cartesianError().matrix()(4,4))<<" "
3210  <<sqrt(trajSOS.cartesianError().matrix()(5,5))<<std::endl;
3211  std::cout<<" q/p dxy/dz xy "<<trajSOS.localParameters().vector()(0)<<" "
3212  <<trajSOS.localParameters().vector()(1)<<" "
3213  <<trajSOS.localParameters().vector()(2)<<" "
3214  <<trajSOS.localParameters().vector()(3)<<" "
3215  <<trajSOS.localParameters().vector()(4)<<std::endl;
3216  std::cout<<" +/- error "<<sqrt(trajSOS.localError().matrix()(0,0))<<" "
3217  <<sqrt(trajSOS.localError().matrix()(1,1))<<" "
3218  <<sqrt(trajSOS.localError().matrix()(2,2))<<" "
3219  <<sqrt(trajSOS.localError().matrix()(3,3))<<" "
3220  <<sqrt(trajSOS.localError().matrix()(4,4))<<" "<<std::endl;
3221 }
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:66
TrajectoryStateOnSurface::localParameters
const LocalTrajectoryParameters & localParameters() const
Definition: TrajectoryStateOnSurface.h:152
TrajectoryStateOnSurface::globalPosition
GlobalPoint globalPosition() const
Definition: TrajectoryStateOnSurface.h:128
LocalTrajectoryParameters::vector
AlgebraicVector5 vector() const
Definition: LocalTrajectoryParameters.h:147
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:61
TrajectoryStateOnSurface::cartesianError
const CartesianTrajectoryError & cartesianError() const
Definition: TrajectoryStateOnSurface.h:146
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:28
LocalTrajectoryError::matrix
const AlgebraicSymMatrix55 & matrix() const
Definition: LocalTrajectoryError.h:50
TrajectoryStateOnSurface::localError
const LocalTrajectoryError & localError() const
Definition: TrajectoryStateOnSurface.h:164
CartesianTrajectoryError::matrix
const AlgebraicSymMatrix66 & matrix() const
Definition: CartesianTrajectoryError.h:40
TrajectoryStateOnSurface::globalParameters
const GlobalTrajectoryParameters & globalParameters() const
Definition: TrajectoryStateOnSurface.h:125
TrajectoryStateOnSurface::globalMomentum
GlobalVector globalMomentum() const
Definition: TrajectoryStateOnSurface.h:131
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:41
GlobalTrajectoryParameters::vector
AlgebraicVector6 vector() const
Definition: GlobalTrajectoryParameters.h:90
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:28
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: Pythia6Hadronizer.cc:86
dbtoconf.out
tuple out
Definition: dbtoconf.py:99
GlobalTrackerMuonAlignment::writeGlPosRcd
void writeGlPosRcd(CLHEP::HepVector &d3)
Definition: GlobalTrackerMuonAlignment.cc:3240
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:41
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 2104 of file GlobalTrackerMuonAlignment.cc.

References funct::A, a, gather_cfg::cout, delta, i, info, j, L1TEmulatorMonitor_cff::p, reco::tau::disc::Pt(), csvReporter::r, asciidump::s, mathSSE::sqrt(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by analyzeTrackTrack(), and analyzeTrackTrajectory().

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

Definition at line 2025 of file GlobalTrackerMuonAlignment.cc.

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

2028 {
2029  // ---------------------------- Calculate Information matrix and Gfree vector
2030  // Information == Hessian , Gfree == gradient of the objective function
2031 
2032  AlgebraicMatrix33 Jac;
2033  AlgebraicVector3 Wi, R_m, R_t, P_t, Norm, dR;
2034 
2035  R_m(0) = Rm.x(); R_m(1) = Rm.y(); R_m(2) = Rm.z();
2036  R_t(0) = Rt.x(); R_t(1) = Rt.y(); R_t(2) = Rt.z();
2037  P_t(0) = Pt.x(); P_t(1) = Pt.y(); P_t(2) = Pt.z();
2038  Norm(0) = Nl.x(); Norm(1) = Nl.y(); Norm(2) = Nl.z();
2039 
2040  for(int i=0; i<=2; i++){
2041  if(Cm(i,i) > 1.e-20)
2042  Wi(i) = 1./Cm(i,i);
2043  else Wi(i) = 1.e-10;
2044  dR(i) = R_m(i) - R_t(i);
2045  }
2046 
2047  float PtN = P_t(0)*Norm(0) + P_t(1)*Norm(1) + P_t(2)*Norm(2);
2048 
2049  Jac(0,0) = 1. - P_t(0)*Norm(0)/PtN;
2050  Jac(0,1) = - P_t(0)*Norm(1)/PtN;
2051  Jac(0,2) = - P_t(0)*Norm(2)/PtN;
2052 
2053  Jac(1,0) = - P_t(1)*Norm(0)/PtN;
2054  Jac(1,1) = 1. - P_t(1)*Norm(1)/PtN;
2055  Jac(1,2) = - P_t(1)*Norm(2)/PtN;
2056 
2057  Jac(2,0) = - P_t(2)*Norm(0)/PtN;
2058  Jac(2,1) = - P_t(2)*Norm(1)/PtN;
2059  Jac(2,2) = 1. - P_t(2)*Norm(2)/PtN;
2060 
2061  AlgebraicSymMatrix33 Itr;
2062 
2063  for(int i=0; i<=2; i++)
2064  for(int j=0; j<=2; j++){
2065  if(j < i) continue;
2066  Itr(i,j) = 0.;
2067  //std::cout<<" ij "<<i<<" "<<j<<std::endl;
2068  for(int k=0; k<=2; k++){
2069  Itr(i,j) += Jac(k,i)*Wi(k)*Jac(k,j);}}
2070 
2071  for(int i=0; i<=2; i++)
2072  for(int j=0; j<=2; j++){
2073  if(j < i) continue;
2074  Inf(i,j) += Itr(i,j);}
2075 
2076  AlgebraicVector3 Gtr(0., 0., 0.);
2077  for(int i=0; i<=2; i++)
2078  for(int k=0; k<=2; k++) Gtr(i) += dR(k)*Wi(k)*Jac(k,i);
2079  for(int i=0; i<=2; i++) Gfr(i) += Gtr(i);
2080 
2081  if(debug_){
2082  std::cout<<" Wi "<<Wi<<std::endl;
2083  std::cout<<" N "<<Norm<<std::endl;
2084  std::cout<<" P_t "<<P_t<<std::endl;
2085  std::cout<<" (Pt*N) "<<PtN<<std::endl;
2086  std::cout<<" dR "<<dR<<std::endl;
2087  std::cout<<" +/- "<<1./sqrt(Wi(0))<<" "<<1./sqrt(Wi(1))<<" "<<1./sqrt(Wi(2))
2088  <<" "<<std::endl;
2089  std::cout<<" Jacobian dr/ddx "<<std::endl;
2090  std::cout<<Jac<<std::endl;
2091  std::cout<<" G-- "<<Gtr<<std::endl;
2092  std::cout<<" Itrack "<<std::endl;
2093  std::cout<<Itr<<std::endl;
2094  std::cout<<" Gfr "<<Gfr<<std::endl;
2095  std::cout<<" -- Inf --"<<std::endl;
2096  std::cout<<Inf<<std::endl;
2097  }
2098 
2099  return;
2100 }
i
int i
Definition: DBlmapReader.cc:9
PV3DBase::y
T y() const
Definition: PV3DBase.h:57
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:28
PV3DBase::z
T z() const
Definition: PV3DBase.h:58
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: Pythia6Hadronizer.cc:86
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:41
PV3DBase::x
T x() const
Definition: PV3DBase.h:56
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 2277 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().

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

Definition at line 2518 of file GlobalTrackerMuonAlignment.cc.

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

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

References funct::A, a, alignmentValidation::c1, funct::cos(), gather_cfg::cout, TrajectoryStateOnSurface::localParameters(), lumiNorm::norm, GloballyPositioned< T >::position(), LocalTrajectoryParameters::pzSign(), GloballyPositioned< T >::rotation(), asciidump::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().

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

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

3035 {
3036  bool info = false;
3037  bool smooth = false;
3038 
3039  if(info) std::cout<<" ......... start of muonFitter ........"<<std::endl;
3040  if(false && info){
3041  this->debugTrackHit(" Muon track hits ", alongTr);
3042  this->debugTrackHit(" Muon TransientTrack hits ", alongTTr);
3043  }
3044 
3045  edm::OwnVector<TrackingRecHit> recHit;
3046  DetId trackDetId = DetId(alongTr->innerDetId());
3047  for (trackingRecHit_iterator i=alongTr->recHitsBegin();i!=alongTr->recHitsEnd(); i++) {
3048  recHit.push_back((*i)->clone());
3049  }
3050  if(direction != alongMomentum)
3051  {
3052  edm::OwnVector<TrackingRecHit> recHitAlong = recHit;
3053  recHit.clear();
3054  for(unsigned int ihit = recHitAlong.size()-1; ihit+1>0; ihit--){
3055  recHit.push_back(recHitAlong[ihit]);
3056  }
3057  recHitAlong.clear();
3058  }
3059  if(info)
3060  std::cout<<" ... Own recHit.size() DetId==Muon "<<recHit.size()<<" "<<trackDetId.rawId()<<std::endl;
3061 
3062  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu =
3063  MuRHBuilder->build(alongTTr.recHitsBegin(), alongTTr.recHitsEnd());
3064  if(info)
3065  std::cout<<" ...along.... recHitMu.size() "<<recHitMu.size()<<std::endl;
3066  if(direction != alongMomentum){
3067  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMuAlong = recHitMu;
3068  recHitMu.clear();
3069  for(unsigned int ihit = recHitMuAlong.size()-1; ihit+1>0; ihit--){
3070  recHitMu.push_back(recHitMuAlong[ihit]);
3071  }
3072  recHitMuAlong.clear();
3073  }
3074  if(info)
3075  std::cout<<" ...opposite... recHitMu.size() "<<recHitMu.size()<<std::endl;
3076 
3077  TrajectoryStateOnSurface firstTSOS;
3078  if(direction == alongMomentum) firstTSOS = alongTTr.innermostMeasurementState();
3079  else firstTSOS = alongTTr.outermostMeasurementState();
3080 
3081  AlgebraicSymMatrix CovLoc(5,1);
3082  CovLoc(1,1) = 1. * firstTSOS.localError().matrix()(0,0);
3083  CovLoc(2,2) = 10. * firstTSOS.localError().matrix()(1,1);
3084  CovLoc(3,3) = 10. * firstTSOS.localError().matrix()(2,2);
3085  CovLoc(4,4) = 100. * firstTSOS.localError().matrix()(3,3);
3086  CovLoc(5,5) = 100. * firstTSOS.localError().matrix()(4,4);
3087  TrajectoryStateOnSurface initialTSOS(firstTSOS.localParameters(), LocalTrajectoryError(CovLoc),
3088  firstTSOS.surface(), &*magneticField_);
3089  TrajectoryStateTransform transformer;
3090  PTrajectoryStateOnDet* PTraj =
3091  transformer.persistentState(initialTSOS, trackDetId.rawId());
3092  //const TrajectorySeed seedT(*PTraj, recHit, alongMomentum);
3093  const TrajectorySeed seedT(*PTraj, recHit, direction);
3094 
3095  std::vector<Trajectory> trajVec;
3096  if(direction == alongMomentum) trajVec = theFitter->fit(seedT, recHitMu, initialTSOS);
3097  else trajVec = theFitterOp->fit(seedT, recHitMu, initialTSOS);
3098 
3099  if(info){
3100  this->debugTrajectorySOSv(" innermostTSOS of TransientTrack",
3101  alongTTr.innermostMeasurementState());
3102  this->debugTrajectorySOSv(" outermostTSOS of TransientTrack",
3103  alongTTr.outermostMeasurementState());
3104  this->debugTrajectorySOS(" initialTSOS for theFitter ", initialTSOS);
3105  std::cout<<" . . . . . trajVec.size() "<<trajVec.size()<<std::endl;
3106  if(trajVec.size()) this->debugTrajectory(" theFitter trajectory ", trajVec[0]);
3107  }
3108 
3109  if(!smooth){
3110  if(trajVec.size()) trackFittedTSOS = trajVec[0].lastMeasurement().updatedState();}
3111  else{
3112  std::vector<Trajectory> trajSVec;
3113  if (trajVec.size()){
3114  if(direction == alongMomentum) trajSVec = theSmoother->trajectories(*(trajVec.begin()));
3115  else trajSVec = theSmootherOp->trajectories(*(trajVec.begin()));
3116  }
3117  if(info) std::cout<<" . . . . trajSVec.size() "<<trajSVec.size()<<std::endl;
3118  if(trajSVec.size()) this->debugTrajectorySOSv("smoothed geom InnermostState",
3119  trajSVec[0].geometricalInnermostState());
3120  if(trajSVec.size()) trackFittedTSOS = trajSVec[0].geometricalInnermostState();
3121  }
3122 
3123 } // 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:3192
TrajectoryStateOnSurface::localParameters
const LocalTrajectoryParameters & localParameters() const
Definition: TrajectoryStateOnSurface.h:152
MuonTransientTrackingRecHitBuilder::ConstRecHitContainer
TransientTrackingRecHit::ConstRecHitContainer ConstRecHitContainer
Definition: MuonTransientTrackingRecHitBuilder.h:17
edm::OwnVector::size
size_type size() const
Definition: OwnVector.h:262
reco::TransientTrack::innermostMeasurementState
TrajectoryStateOnSurface innermostMeasurementState() const
Definition: TransientTrack.h:51
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:290
GlobalTrackerMuonAlignment::debugTrackHit
void debugTrackHit(const std::string, reco::TrackRef)
Definition: GlobalTrackerMuonAlignment.cc:3144
KFTrajectoryFitter::fit
virtual std::vector< Trajectory > fit(const Trajectory &aTraj) const
Definition: KFTrajectoryFitter.cc:20
edm::OwnVector::clear
void clear()
Definition: OwnVector.h:399
LocalTrajectoryError::matrix
const AlgebraicSymMatrix55 & matrix() const
Definition: LocalTrajectoryError.h:50
TrajectoryStateOnSurface::localError
const LocalTrajectoryError & localError() const
Definition: TrajectoryStateOnSurface.h:164
reco::TransientTrack::outermostMeasurementState
TrajectoryStateOnSurface outermostMeasurementState() const
Definition: TransientTrack.h:48
reco::TransientTrack::recHitsEnd
trackingRecHit_iterator recHitsEnd() const
last iterator to RecHits
Definition: TransientTrack.h:96
DetId
Definition: DetId.h:20
GlobalTrackerMuonAlignment::debugTrajectory
void debugTrajectory(const std::string, Trajectory &)
Definition: GlobalTrackerMuonAlignment.cc:3224
TrajectoryStateOnSurface::surface
const Surface & surface() const
Definition: TrajectoryStateOnSurface.h:167
AlgebraicSymMatrix
CLHEP::HepSymMatrix AlgebraicSymMatrix
Definition: AlgebraicObjects.h:14
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:41
reco::TransientTrack::recHitsBegin
trackingRecHit_iterator recHitsBegin() const
first iterator to RecHits
Definition: TransientTrack.h:94
GlobalTrackerMuonAlignment::MuRHBuilder
MuonTransientTrackingRecHitBuilder * MuRHBuilder
Definition: GlobalTrackerMuonAlignment.cc:187
GlobalTrackerMuonAlignment::debugTrajectorySOS
void debugTrajectorySOS(const std::string, TrajectoryStateOnSurface &)
Definition: GlobalTrackerMuonAlignment.cc:3160
void GlobalTrackerMuonAlignment::trackFitter ( reco::TrackRef  alongTr,
reco::TransientTrack &  alongTTr,
PropagationDirection  direction,
TrajectoryStateOnSurface &  trackFittedTSOS 
)

Definition at line 2918 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(), edm::OwnVector< T, P >::push_back(), DetId::rawId(), edm::OwnVector< T, P >::size(), TrajectoryStateOnSurface::surface(), and RSFinalFitAnalytical_cff::TTRHBuilder.

Referenced by analyzeTrackTrajectory().

2921 {
2922  bool info = false;
2923  bool smooth = false;
2924 
2925  if(info)
2926  std::cout<<" ......... start of trackFitter ......... "<<std::endl;
2927  if(false && info){
2928  this->debugTrackHit(" Tracker track hits ", alongTr);
2929  this->debugTrackHit(" Tracker TransientTrack hits ", alongTTr);
2930  }
2931 
2932  edm::OwnVector<TrackingRecHit> recHit;
2933  DetId trackDetId = DetId(alongTr->innerDetId());
2934  for (trackingRecHit_iterator i=alongTr->recHitsBegin();i!=alongTr->recHitsEnd(); i++) {
2935  recHit.push_back((*i)->clone());
2936  }
2937  if(direction != alongMomentum)
2938  {
2939  edm::OwnVector<TrackingRecHit> recHitAlong = recHit;
2940  recHit.clear();
2941  for(unsigned int ihit = recHitAlong.size()-1; ihit+1>0; ihit--){
2942  recHit.push_back(recHitAlong[ihit]);
2943  }
2944  recHitAlong.clear();
2945  }
2946  if(info)
2947  std::cout<<" ... Own recHit.size() "<<recHit.size()<<" "<<trackDetId.rawId()<<std::endl;
2948 
2949  //MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitTrack;
2950  TransientTrackingRecHit::RecHitContainer recHitTrack;
2951  for (trackingRecHit_iterator i=alongTr->recHitsBegin();i!=alongTr->recHitsEnd(); i++) {
2952  if((*i)->isValid()){
2953  recHitTrack.push_back(TTRHBuilder->build(&**i ));}
2954  else{
2955  recHitTrack.push_back(TTRHBuilder->build(&**i ));}
2956  }
2957 
2958  if(info)
2959  std::cout<<" ... recHitTrack.size() "<<recHit.size()<<" "<<trackDetId.rawId()
2960  <<" ======> recHitMu "<<std::endl;
2961 
2962  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu = recHitTrack;
2963  /*
2964  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu =
2965  MuRHBuilder->build(alongTTr.recHitsBegin(), alongTTr.recHitsEnd());
2966  */
2967  if(info)
2968  std::cout<<" ...along.... recHitMu.size() "<<recHitMu.size()<<std::endl;
2969  if(direction != alongMomentum){
2970  MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMuAlong = recHitMu;
2971  recHitMu.clear();
2972  for(unsigned int ihit = recHitMuAlong.size()-1; ihit+1>0; ihit--){
2973  recHitMu.push_back(recHitMuAlong[ihit]);
2974  }
2975  recHitMuAlong.clear();
2976  }
2977  if(info)
2978  std::cout<<" ...opposite... recHitMu.size() "<<recHitMu.size()<<std::endl;
2979 
2980  TrajectoryStateOnSurface firstTSOS;
2981  if(direction == alongMomentum) firstTSOS = alongTTr.innermostMeasurementState();
2982  else firstTSOS = alongTTr.outermostMeasurementState();
2983 
2984  AlgebraicSymMatrix CovLoc(5,1);
2985  CovLoc(1,1) = 1. * firstTSOS.localError().matrix()(0,0);
2986  CovLoc(2,2) = 10. * firstTSOS.localError().matrix()(1,1);
2987  CovLoc(3,3) = 10. * firstTSOS.localError().matrix()(2,2);
2988  CovLoc(4,4) = 100. * firstTSOS.localError().matrix()(3,3);
2989  CovLoc(5,5) = 100. * firstTSOS.localError().matrix()(4,4);
2990  TrajectoryStateOnSurface initialTSOS(firstTSOS.localParameters(), LocalTrajectoryError(CovLoc),
2991  firstTSOS.surface(), &*magneticField_);
2992  TrajectoryStateTransform transformer;
2993  PTrajectoryStateOnDet* PTraj =
2994  transformer.persistentState(initialTSOS, trackDetId.rawId());
2995  //const TrajectorySeed seedT(*PTraj, recHit, alongMomentum);
2996  const TrajectorySeed seedT(*PTraj, recHit, direction);
2997 
2998  std::vector<Trajectory> trajVec;
2999  if(direction == alongMomentum) trajVec = theFitter->fit(seedT, recHitMu, initialTSOS);
3000  else trajVec = theFitterOp->fit(seedT, recHitMu, initialTSOS);
3001 
3002  if(info){
3003  this->debugTrajectorySOSv(" innermostTSOS of TransientTrack",
3004  alongTTr.innermostMeasurementState());
3005  this->debugTrajectorySOSv(" outermostTSOS of TransientTrack",
3006  alongTTr.outermostMeasurementState());
3007  this->debugTrajectorySOS(" initialTSOS for theFitter ", initialTSOS);
3008  std::cout<<" . . . . . trajVec.size() "<<trajVec.size()<<std::endl;
3009  if(trajVec.size()) this->debugTrajectory(" theFitter trajectory ", trajVec[0]);
3010  }
3011 
3012  if(!smooth){
3013  if(trajVec.size()) trackFittedTSOS = trajVec[0].lastMeasurement().updatedState();}
3014  else{
3015  std::vector<Trajectory> trajSVec;
3016  if (trajVec.size()){
3017  if(direction == alongMomentum) trajSVec = theSmoother->trajectories(*(trajVec.begin()));
3018  else trajSVec = theSmootherOp->trajectories(*(trajVec.begin()));
3019  }
3020  if(info) std::cout<<" . . . . trajSVec.size() "<<trajSVec.size()<<std::endl;
3021  if(trajSVec.size()) this->debugTrajectorySOSv("smoothed geom InnermostState",
3022  trajSVec[0].geometricalInnermostState());
3023  if(trajSVec.size()) trackFittedTSOS = trajSVec[0].geometricalInnermostState();
3024  }
3025 
3026  if(info) this->debugTrajectorySOSv(" trackFittedTSOS ", trackFittedTSOS);
3027 
3028 } // end of trackFitter
i
int i
Definition: DBlmapReader.cc:9
GlobalTrackerMuonAlignment::debugTrajectorySOSv
void debugTrajectorySOSv(const std::string, TrajectoryStateOnSurface)
Definition: GlobalTrackerMuonAlignment.cc:3192
TrajectoryStateOnSurface::localParameters
const LocalTrajectoryParameters & localParameters() const
Definition: TrajectoryStateOnSurface.h:152
MuonTransientTrackingRecHitBuilder::ConstRecHitContainer
TransientTrackingRecHit::ConstRecHitContainer ConstRecHitContainer
Definition: MuonTransientTrackingRecHitBuilder.h:17
edm::OwnVector::size
size_type size() const
Definition: OwnVector.h:262
TransientTrackingRecHit::RecHitContainer
std::vector< ConstRecHitPointer > RecHitContainer
Definition: TransientTrackingRecHit.h:22
reco::TransientTrack::innermostMeasurementState
TrajectoryStateOnSurface innermostMeasurementState() const
Definition: TransientTrack.h:51
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:290
GlobalTrackerMuonAlignment::debugTrackHit
void debugTrackHit(const std::string, reco::TrackRef)
Definition: GlobalTrackerMuonAlignment.cc:3144
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:399
LocalTrajectoryError::matrix
const AlgebraicSymMatrix55 & matrix() const
Definition: LocalTrajectoryError.h:50
TrajectoryStateOnSurface::localError
const LocalTrajectoryError & localError() const
Definition: TrajectoryStateOnSurface.h:164
reco::TransientTrack::outermostMeasurementState
TrajectoryStateOnSurface outermostMeasurementState() const
Definition: TransientTrack.h:48
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:3224
TrajectoryStateOnSurface::surface
const Surface & surface() const
Definition: TrajectoryStateOnSurface.h:167
AlgebraicSymMatrix
CLHEP::HepSymMatrix AlgebraicSymMatrix
Definition: AlgebraicObjects.h:14
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:41
GlobalTrackerMuonAlignment::debugTrajectorySOS
void debugTrajectorySOS(const std::string, TrajectoryStateOnSurface &)
Definition: GlobalTrackerMuonAlignment.cc:3160
void GlobalTrackerMuonAlignment::writeGlPosRcd ( CLHEP::HepVector &  d3)

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

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

char GlobalTrackerMuonAlignment::MuSelect[100]
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