00001
00002
00003
00004
00005
00015
00016
00017
00018
00019
00020
00021
00022
00023 #include <memory>
00024 #include <string>
00025 #include <map>
00026 #include <vector>
00027 #include <fstream>
00028 #include <typeinfo>
00029
00030
00031
00032
00033 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00034 #include "FWCore/Framework/interface/EDAnalyzer.h"
00035 #include "FWCore/Framework/interface/ESWatcher.h"
00036 #include "FWCore/Framework/interface/Frameworkfwd.h"
00037
00038 #include "FWCore/Framework/interface/Event.h"
00039 #include "FWCore/Framework/interface/EventSetup.h"
00040 #include "FWCore/Framework/interface/MakerMacros.h"
00041 #include "FWCore/Utilities/interface/InputTag.h"
00042 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00043 #include "FWCore/ServiceRegistry/interface/Service.h"
00044
00045
00046
00047 #include "DataFormats/DetId/interface/DetId.h"
00048 #include "DataFormats/MuonReco/interface/Muon.h"
00049 #include "DataFormats/MuonReco/interface/MuonFwd.h"
00050 #include "DataFormats/TrackReco/interface/Track.h"
00051 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00052 #include "DataFormats/GeometrySurface/interface/TangentPlane.h"
00053 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00054 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
00055 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
00056 #include "Geometry/TrackingGeometryAligner/interface/GeometryAligner.h"
00057 #include "MagneticField/Engine/interface/MagneticField.h"
00058 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00059
00060 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBuilder.h"
00061 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00062 #include "TrackingTools/TrackRefitter/interface/TrackTransformer.h"
00063 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00064 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
00065 #include "TrackingTools/GeomPropagators/interface/SmartPropagator.h"
00066 #include "TrackingTools/GeomPropagators/interface/PropagationDirectionChooser.h"
00067 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
00068 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
00069 #include "TrackPropagation/RungeKutta/interface/defaultRKPropagator.h"
00070 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00071 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
00072 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCartesianToLocal.h"
00073
00074
00075 #include "FWCore/Framework/interface/ESHandle.h"
00076 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
00077
00078
00079 #include "CondFormats/Alignment/interface/Alignments.h"
00080 #include "CondFormats/Alignment/interface/AlignTransform.h"
00081 #include "CondFormats/AlignmentRecord/interface/GlobalPositionRcd.h"
00082
00083
00084 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
00085 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
00086 #include "TrackingTools/TrackFitters/interface/KFTrajectoryFitter.h"
00087 #include "TrackingTools/TrackFitters/interface/KFTrajectorySmoother.h"
00088 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
00089 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
00090
00091 #include <TH1.h>
00092 #include <TH2.h>
00093 #include <TFile.h>
00094
00095 #include<iostream>
00096 #include<cmath>
00097 #include<assert.h>
00098 #include<fstream>
00099 #include<iomanip>
00100
00101 using namespace edm;
00102 using namespace std;
00103 using namespace reco;
00104
00105
00106
00107
00108
00109 class GlobalTrackerMuonAlignment : public edm::EDAnalyzer {
00110 public:
00111 explicit GlobalTrackerMuonAlignment(const edm::ParameterSet&);
00112 ~GlobalTrackerMuonAlignment();
00113
00114 void analyzeTrackTrack(const edm::Event&, const edm::EventSetup&);
00115 void analyzeTrackTrajectory(const edm::Event&, const edm::EventSetup&);
00116 void bookHist();
00117 void fitHist();
00118 void trackFitter(reco::TrackRef, reco::TransientTrack&, PropagationDirection,
00119 TrajectoryStateOnSurface&);
00120 void muonFitter(reco::TrackRef, reco::TransientTrack&, PropagationDirection,
00121 TrajectoryStateOnSurface&);
00122 void debugTrackHit(const std::string, reco::TrackRef);
00123 void debugTrackHit(const std::string, reco::TransientTrack&);
00124 void debugTrajectorySOS(const std::string, TrajectoryStateOnSurface&);
00125 void debugTrajectorySOSv(const std::string, TrajectoryStateOnSurface);
00126 void debugTrajectory(const std::string, Trajectory&);
00127
00128 void gradientGlobal(GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&,
00129 GlobalVector&, AlgebraicSymMatrix66&);
00130 void gradientLocal(GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&,
00131 GlobalVector&, CLHEP::HepSymMatrix&, CLHEP::HepMatrix&,
00132 CLHEP::HepVector&, AlgebraicVector4&);
00133 void gradientGlobalAlg(GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&,
00134 AlgebraicSymMatrix66&);
00135 void misalignMuon(GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&,
00136 GlobalVector&, GlobalVector&);
00137 void misalignMuonL(GlobalVector&, GlobalVector&, GlobalVector&, GlobalVector&,
00138 GlobalVector&, GlobalVector&, AlgebraicVector4&,
00139 TrajectoryStateOnSurface&, TrajectoryStateOnSurface&);
00140 void writeGlPosRcd(CLHEP::HepVector& d3);
00141 inline double CLHEP_dot(CLHEP::HepVector a, CLHEP::HepVector b)
00142 {
00143 return a(1)*b(1)+a(2)*b(2)+a(3)*b(3);
00144 }
00145
00146 private:
00147
00148 virtual void beginJob() ;
00149 virtual void analyze(const edm::Event&, const edm::EventSetup&);
00150 virtual void endJob() ;
00151
00152
00153
00154 edm::InputTag trackTags_;
00155 edm::InputTag muonTags_;
00156 edm::InputTag gmuonTags_;
00157 edm::InputTag smuonTags_;
00158 string propagator_;
00159 bool cosmicMuonMode_;
00160 bool isolatedMuonMode_;
00161 bool collectionCosmic;
00162 bool collectionIsolated;
00163 bool refitMuon_;
00164 bool refitTrack_;
00165 string rootOutFile_;
00166 string txtOutFile_;
00167 bool writeDB_;
00168 bool debug_;
00169
00170 edm::ESWatcher<GlobalTrackingGeometryRecord> watchTrackingGeometry_;
00171 edm::ESHandle<GlobalTrackingGeometry> theTrackingGeometry;
00172 const GlobalTrackingGeometry* trackingGeometry_;
00173
00174 edm::ESWatcher<IdealMagneticFieldRecord> watchMagneticFieldRecord_;
00175 const MagneticField* magneticField_;
00176
00177 edm::ESWatcher<TrackingComponentsRecord> watchTrackingComponentsRecord_;
00178
00179 edm::ESWatcher<GlobalPositionRcd> watchGlobalPositionRcd_;
00180 const Alignments* globalPositionRcd_;
00181
00182 KFTrajectoryFitter* theFitter;
00183 KFTrajectorySmoother* theSmoother;
00184 KFTrajectoryFitter* theFitterOp;
00185 KFTrajectorySmoother* theSmootherOp;
00186 bool defineFitter;
00187 MuonTransientTrackingRecHitBuilder* MuRHBuilder;
00188 const TransientTrackingRecHitBuilder* TTRHBuilder;
00189
00190
00191 AlgebraicVector3 Gfr;
00192 AlgebraicSymMatrix33 Inf;
00193
00194 CLHEP::HepVector Grad;
00195 CLHEP::HepMatrix Hess;
00196
00197 CLHEP::HepVector GradL;
00198 CLHEP::HepMatrix HessL;
00199
00200 int N_event;
00201 int N_track;
00202
00203 std::vector<AlignTransform>::const_iterator
00204 iteratorTrackerRcd;
00205 std::vector<AlignTransform>::const_iterator
00206 iteratorMuonRcd;
00207 std::vector<AlignTransform>::const_iterator
00208 iteratorEcalRcd;
00209 std::vector<AlignTransform>::const_iterator
00210 iteratorHcalRcd;
00211
00212 CLHEP::HepVector MuGlShift;
00213 CLHEP::HepVector MuGlAngle;
00214
00215 std::string MuSelect;
00216
00217 ofstream OutGlobalTxt;
00218
00219 TFile * file;
00220 TH1F * histo;
00221 TH1F * histo2;
00222 TH1F * histo3;
00223 TH1F * histo4;
00224 TH1F * histo5;
00225 TH1F * histo6;
00226 TH1F * histo7;
00227 TH1F * histo8;
00228
00229 TH1F * histo11;
00230 TH1F * histo12;
00231 TH1F * histo13;
00232 TH1F * histo14;
00233 TH1F * histo15;
00234 TH1F * histo16;
00235 TH1F * histo17;
00236 TH1F * histo18;
00237 TH1F * histo19;
00238 TH1F * histo20;
00239 TH1F * histo21;
00240 TH1F * histo22;
00241 TH1F * histo23;
00242 TH1F * histo24;
00243 TH1F * histo25;
00244 TH1F * histo26;
00245 TH1F * histo27;
00246 TH1F * histo28;
00247 TH1F * histo29;
00248 TH1F * histo30;
00249 TH1F * histo31;
00250 TH1F * histo32;
00251 TH1F * histo33;
00252 TH1F * histo34;
00253 TH1F * histo35;
00254
00255 TH2F * histo101;
00256 TH2F * histo102;
00257 };
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271 GlobalTrackerMuonAlignment::GlobalTrackerMuonAlignment(const edm::ParameterSet& iConfig)
00272 :trackTags_(iConfig.getParameter<edm::InputTag>("tracks")),
00273 muonTags_(iConfig.getParameter<edm::InputTag>("muons")),
00274 gmuonTags_(iConfig.getParameter<edm::InputTag>("gmuons")),
00275 smuonTags_(iConfig.getParameter<edm::InputTag>("smuons")),
00276 propagator_(iConfig.getParameter<string>("Propagator")),
00277
00278 cosmicMuonMode_(iConfig.getParameter<bool>("cosmics")),
00279 isolatedMuonMode_(iConfig.getParameter<bool>("isolated")),
00280
00281 refitMuon_(iConfig.getParameter<bool>("refitmuon")),
00282 refitTrack_(iConfig.getParameter<bool>("refittrack")),
00283
00284 rootOutFile_(iConfig.getUntrackedParameter<string>("rootOutFile")),
00285 txtOutFile_(iConfig.getUntrackedParameter<string>("txtOutFile")),
00286 writeDB_(iConfig.getUntrackedParameter<bool>("writeDB")),
00287 debug_(iConfig.getUntrackedParameter<bool>("debug")),
00288 defineFitter(true)
00289 {
00290 #ifdef NO_FALSE_FALSE_MODE
00291 if (cosmicMuonMode_==false && isolatedMuonMode_==false) {
00292 throw cms::Exception("BadConfig") << "Muon collection not set! "
00293 << "Use from GlobalTrackerMuonAlignment_test_cfg.py.";
00294 }
00295 #endif
00296 if (cosmicMuonMode_==true && isolatedMuonMode_==true) {
00297 throw cms::Exception("BadConfig") << "Muon collection can be either cosmic or isolated! "
00298 << "Please set only one to true.";
00299 }
00300 }
00301
00302 GlobalTrackerMuonAlignment::~GlobalTrackerMuonAlignment()
00303 {
00304
00305
00306
00307
00308 }
00309
00310
00311
00312
00313
00314
00315 void
00316 GlobalTrackerMuonAlignment::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00317 {
00318
00319 GlobalTrackerMuonAlignment::analyzeTrackTrajectory(iEvent, iSetup);
00320 }
00321
00322
00323 void GlobalTrackerMuonAlignment::beginJob()
00324 {
00325 N_event = 0;
00326 N_track = 0;
00327
00328 if (cosmicMuonMode_==true && isolatedMuonMode_==false) {
00329 collectionCosmic = true; collectionIsolated = false;}
00330 else if (cosmicMuonMode_==false && isolatedMuonMode_==true) {
00331 collectionCosmic = false; collectionIsolated = true;}
00332 else {
00333 collectionCosmic = false; collectionIsolated = false;}
00334
00335 for(int i=0; i<=2; i++){
00336 Gfr(i) = 0.;
00337 for(int j=0; j<=2; j++){
00338 Inf(i,j) = 0.;
00339 }}
00340
00341 Grad = CLHEP::HepVector(6,0);
00342 Hess = CLHEP::HepSymMatrix(6,0);
00343
00344 GradL = CLHEP::HepVector(6,0);
00345 HessL = CLHEP::HepSymMatrix(6,0);
00346
00347
00348 TDirectory * dirsave = gDirectory;
00349
00350 file = new TFile(rootOutFile_.c_str(), "recreate");
00351 const bool oldAddDir = TH1::AddDirectoryStatus();
00352
00353 TH1::AddDirectory(true);
00354
00355 this->bookHist();
00356
00357 TH1::AddDirectory(oldAddDir);
00358 dirsave->cd();
00359 }
00360
00361
00362 void GlobalTrackerMuonAlignment::endJob() {
00363 bool alarm = false;
00364
00365 this->fitHist();
00366
00367 AlgebraicVector3 d(0., 0., 0.);
00368
00369 AlgebraicSymMatrix33 InfI;
00370 for(int i=0; i<=2; i++)
00371 for(int j=0; j<=2; j++){
00372 if(j < i) continue;
00373 InfI(i,j) += Inf(i,j);}
00374 bool ierr = !InfI.Invert();
00375 if( ierr ) { if(alarm)
00376 std::cout<<" Error inverse Inf matrix !!!!!!!!!!!"<<std::endl;}
00377
00378 for(int i=0; i<=2; i++)
00379 for(int k=0; k<=2; k++) d(i) -= InfI(i,k)*Gfr(k);
00380
00381
00382
00383 CLHEP::HepVector d3 = CLHEP::solve(Hess, -Grad);
00384 int iEr3;
00385 CLHEP::HepMatrix Errd3 = Hess.inverse(iEr3);
00386 if( iEr3 != 0) { if(alarm)
00387 std::cout<<" endJob Error inverse Hess matrix !!!!!!!!!!!"<<std::endl;
00388 }
00389
00390
00391
00392 CLHEP::HepVector dLI = CLHEP::solve(HessL, -GradL);
00393 int iErI;
00394 CLHEP::HepMatrix ErrdLI = HessL.inverse(iErI);
00395 if( iErI != 0) { if(alarm)
00396 std::cout<<" endJob Error inverse HessL matrix !!!!!!!!!!!"<<std::endl;
00397 }
00398
00399
00400
00401 std::cout<<" ---- "<<N_event<<" event "<<N_track<<" tracks "
00402 <<MuSelect<<" ---- "<<std::endl;
00403 if (collectionIsolated) std::cout<<" ALCARECOMuAlCalIsolatedMu "<<std::endl;
00404 else if (collectionCosmic) std::cout<<" ALCARECOMuAlGlobalCosmics "<<std::endl;
00405 else std::cout<<smuonTags_<<std::endl;
00406 std::cout<<" Similated shifts[cm] "
00407 <<MuGlShift(1)<<" "<<MuGlShift(2)<<" "<<MuGlShift(3)<<" "
00408 <<" angles[rad] "<<MuGlAngle(1)<<" "<<MuGlAngle(2)<<" "<<MuGlAngle(3)<<" "
00409 <<std::endl;
00410 std::cout<<" d "<<-d<<std::endl;;
00411 std::cout<<" +- "<<sqrt(InfI(0,0))<<" "<<sqrt(InfI(1,1))<<" "<<sqrt(InfI(2,2))
00412 <<std::endl;
00413 std::cout<<" dG "<<d3(1)<<" "<<d3(2)<<" "<<d3(3)<<" "<<d3(4)<<" "
00414 <<d3(5)<<" "<<d3(6)<<std::endl;;
00415 std::cout<<" +- "<<sqrt(Errd3(1,1))<<" "<<sqrt(Errd3(2,2))<<" "<<sqrt(Errd3(3,3))
00416 <<" "<<sqrt(Errd3(4,4))<<" "<<sqrt(Errd3(5,5))<<" "<<sqrt(Errd3(6,6))
00417 <<std::endl;
00418 std::cout<<" dL "<<dLI(1)<<" "<<dLI(2)<<" "<<dLI(3)<<" "<<dLI(4)<<" "
00419 <<dLI(5)<<" "<<dLI(6)<<std::endl;;
00420 std::cout<<" +- "<<sqrt(ErrdLI(1,1))<<" "<<sqrt(ErrdLI(2,2))<<" "<<sqrt(ErrdLI(3,3))
00421 <<" "<<sqrt(ErrdLI(4,4))<<" "<<sqrt(ErrdLI(5,5))<<" "<<sqrt(ErrdLI(6,6))
00422 <<std::endl;
00423
00424
00425 CLHEP::HepVector vectorToDb(6,0), vectorErrToDb(6,0);
00426
00427
00428 vectorToDb = -dLI;
00429 for(unsigned int i=1; i<=6; i++) vectorErrToDb(i) = sqrt(ErrdLI(i,i));
00430
00431
00432 file->Write();
00433 file->Close();
00434
00435
00436 OutGlobalTxt.open(txtOutFile_.c_str(), ios::out);
00437 if(!OutGlobalTxt.is_open()) std::cout<<" outglobal.txt is not open !!!!!"<<std::endl;
00438 else{
00439 OutGlobalTxt<<vectorToDb(1)<<" "<<vectorToDb(2)<<" "<<vectorToDb(3)<<" "<<vectorToDb(4)<<" "
00440 <<vectorToDb(5)<<" "<<vectorToDb(6)<<" muon Global.\n";
00441 OutGlobalTxt<<vectorErrToDb(1)<<" "<<vectorErrToDb(1)<<" "<<vectorErrToDb(1)<<" "
00442 <<vectorErrToDb(1)<<" "<<vectorErrToDb(1)<<" "<<vectorErrToDb(1)<<" errors.\n";
00443 OutGlobalTxt<<N_event<<" events are processed.\n";
00444
00445 if (collectionIsolated) OutGlobalTxt<<"ALCARECOMuAlCalIsolatedMu.\n";
00446 else if (collectionCosmic) OutGlobalTxt<<" ALCARECOMuAlGlobalCosmics.\n";
00447 else OutGlobalTxt<<smuonTags_<<".\n";
00448 OutGlobalTxt.close();
00449 std::cout<<" Write to the file outglobal.txt done "<<std::endl;
00450 }
00451
00452
00453 if(debug_)
00454 std::cout<<" writeBD_ "<<writeDB_<<std::endl;
00455 if (writeDB_)
00456 GlobalTrackerMuonAlignment::writeGlPosRcd(vectorToDb);
00457 }
00458
00459
00460 void GlobalTrackerMuonAlignment::bookHist()
00461 {
00462 double PI = 3.1415927;
00463 histo = new TH1F("Pt","pt",1000,0,100);
00464 histo2 = new TH1F("P","P [GeV/c]",400,0.,400.);
00465 histo2->GetXaxis()->SetTitle("momentum [GeV/c]");
00466 histo3 = new TH1F("outerLambda","#lambda outer",100,-PI/2.,PI/2.);
00467 histo3->GetXaxis()->SetTitle("#lambda outer");
00468 histo4 = new TH1F("phi","#phi [rad]",100,-PI,PI);
00469 histo4->GetXaxis()->SetTitle("#phi [rad]");
00470 histo5 = new TH1F("Rmuon","inner muon hit R [cm]",100,0.,800.);
00471 histo5->GetXaxis()->SetTitle("R of muon [cm]");
00472 histo6 = new TH1F("Zmuon","inner muon hit Z[cm]",100,-1000.,1000.);
00473 histo6->GetXaxis()->SetTitle("Z of muon [cm]");
00474 histo7 = new TH1F("(Pm-Pt)/Pt"," (Pmuon-Ptrack)/Ptrack", 100,-2.,2.);
00475 histo7->GetXaxis()->SetTitle("(Pmuon-Ptrack)/Ptrack");
00476 histo8 = new TH1F("chi muon-track","#chi^{2}(muon-track)", 1000,0.,1000.);
00477 histo8->GetXaxis()->SetTitle("#chi^{2} of muon w.r.t. propagated track");
00478 histo11 = new TH1F("distance muon-track","distance muon w.r.t track [cm]",100,0.,30.);
00479 histo11->GetXaxis()->SetTitle("distance of muon w.r.t. track [cm]");
00480 histo12 = new TH1F("Xmuon-Xtrack","Xmuon-Xtrack [cm]",200,-20.,20.);
00481 histo12->GetXaxis()->SetTitle("Xmuon - Xtrack [cm]");
00482 histo13 = new TH1F("Ymuon-Ytrack","Ymuon-Ytrack [cm]",200,-20.,20.);
00483 histo13->GetXaxis()->SetTitle("Ymuon - Ytrack [cm]");
00484 histo14 = new TH1F("Zmuon-Ztrack","Zmuon-Ztrack [cm]",200,-20.,20.);
00485 histo14->GetXaxis()->SetTitle("Zmuon-Ztrack [cm]");
00486 histo15 = new TH1F("NXmuon-NXtrack","NXmuon-NXtrack [rad]",200,-.1,.1);
00487 histo15->GetXaxis()->SetTitle("N_{X}(muon)-N_{X}(track) [rad]");
00488 histo16 = new TH1F("NYmuon-NYtrack","NYmuon-NYtrack [rad]",200,-.1,.1);
00489 histo16->GetXaxis()->SetTitle("N_{Y}(muon)-N_{Y}(track) [rad]");
00490 histo17 = new TH1F("NZmuon-NZtrack","NZmuon-NZtrack [rad]",200,-.1,.1);
00491 histo17->GetXaxis()->SetTitle("N_{Z}(muon)-N_{Z}(track) [rad]");
00492 histo18 = new TH1F("expected error of Xinner","outer hit of inner tracker",100,0,.01);
00493 histo18->GetXaxis()->SetTitle("expected error of Xinner [cm]");
00494 histo19 = new TH1F("expected error of Xmuon","inner hit of muon",100,0,.1);
00495 histo19->GetXaxis()->SetTitle("expected error of Xmuon [cm]");
00496 histo20 = new TH1F("expected error of Xmuon-Xtrack","muon w.r.t. propagated track",100,0.,10.);
00497 histo20->GetXaxis()->SetTitle("expected error of Xmuon-Xtrack [cm]");
00498 histo21 = new TH1F("pull of Xmuon-Xtrack","pull of Xmuon-Xtrack",100,-10.,10.);
00499 histo21->GetXaxis()->SetTitle("(Xmuon-Xtrack)/expected error");
00500 histo22 = new TH1F("pull of Ymuon-Ytrack","pull of Ymuon-Ytrack",100,-10.,10.);
00501 histo22->GetXaxis()->SetTitle("(Ymuon-Ytrack)/expected error");
00502 histo23 = new TH1F("pull of Zmuon-Ztrack","pull of Zmuon-Ztrack",100,-10.,10.);
00503 histo23->GetXaxis()->SetTitle("(Zmuon-Ztrack)/expected error");
00504 histo24 = new TH1F("pull of PXmuon-PXtrack","pull of PXmuon-PXtrack",100,-10.,10.);
00505 histo24->GetXaxis()->SetTitle("(P_{X}(muon)-P_{X}(track))/expected error");
00506 histo25 = new TH1F("pull of PYmuon-PYtrack","pull of PYmuon-PYtrack",100,-10.,10.);
00507 histo25->GetXaxis()->SetTitle("(P_{Y}(muon)-P_{Y}(track))/expected error");
00508 histo26 = new TH1F("pull of PZmuon-PZtrack","pull of PZmuon-PZtrack",100,-10.,10.);
00509 histo26->GetXaxis()->SetTitle("(P_{Z}(muon)-P_{Z}(track))/expected error");
00510 histo27 = new TH1F("N_x","Nx of tangent plane",120,-1.1,1.1);
00511 histo27->GetXaxis()->SetTitle("normal vector projection N_{X}");
00512 histo28 = new TH1F("N_y","Ny of tangent plane",120,-1.1,1.1);
00513 histo28->GetXaxis()->SetTitle("normal vector projection N_{Y}");
00514 histo29 = new TH1F("lenght of track","lenght of track",200,0.,400);
00515 histo29->GetXaxis()->SetTitle("lenght of track [cm]");
00516 histo30 = new TH1F("lenght of muon","lenght of muon",200,0.,800);
00517 histo30->GetXaxis()->SetTitle("lenght of muon [cm]");
00518
00519 histo31 = new TH1F("local chi muon-track","#local chi^{2}(muon-track)", 1000,0.,1000.);
00520 histo31->GetXaxis()->SetTitle("#local chi^{2} of muon w.r.t. propagated track");
00521 histo32 = new TH1F("pull of Px/Pz local","pull of Px/Pz local",100,-10.,10.);
00522 histo32->GetXaxis()->SetTitle("local (Px/Pz(muon) - Px/Pz(track))/expected error");
00523 histo33 = new TH1F("pull of Py/Pz local","pull of Py/Pz local",100,-10.,10.);
00524 histo33->GetXaxis()->SetTitle("local (Py/Pz(muon) - Py/Pz(track))/expected error");
00525 histo34 = new TH1F("pull of X local","pull of X local",100,-10.,10.);
00526 histo34->GetXaxis()->SetTitle("local (Xmuon - Xtrack)/expected error");
00527 histo35 = new TH1F("pull of Y local","pull of Y local",100,-10.,10.);
00528 histo35->GetXaxis()->SetTitle("local (Ymuon - Ytrack)/expected error");
00529
00530 histo101 = new TH2F("Rtr/mu vs Ztr/mu","hit of track/muon",100,-800.,800.,100,0.,600.);
00531 histo101->GetXaxis()->SetTitle("Z of track/muon [cm]");
00532 histo101->GetYaxis()->SetTitle("R of track/muon [cm]");
00533 histo102 = new TH2F("Ytr/mu vs Xtr/mu","hit of track/muon",100,-600.,600.,100,-600.,600.);
00534 histo102->GetXaxis()->SetTitle("X of track/muon [cm]");
00535 histo102->GetYaxis()->SetTitle("Y of track/muon [cm]");
00536 }
00537
00538
00539 void GlobalTrackerMuonAlignment::fitHist() {
00540
00541 histo7->Fit("gaus","Q");
00542
00543 histo12->Fit("gaus","Q");
00544 histo13->Fit("gaus","Q");
00545 histo14->Fit("gaus","Q");
00546 histo15->Fit("gaus","Q");
00547 histo16->Fit("gaus","Q");
00548 histo17->Fit("gaus","Q");
00549
00550 histo21->Fit("gaus","Q");
00551 histo22->Fit("gaus","Q");
00552 histo23->Fit("gaus","Q");
00553 histo24->Fit("gaus","Q");
00554 histo25->Fit("gaus","Q");
00555 histo26->Fit("gaus","Q");
00556
00557 histo32->Fit("gaus","Q");
00558 histo33->Fit("gaus","Q");
00559 histo34->Fit("gaus","Q");
00560 histo35->Fit("gaus","Q");
00561 }
00562
00563
00564 void GlobalTrackerMuonAlignment::analyzeTrackTrack
00565 (const edm::Event& iEvent, const edm::EventSetup& iSetup)
00566 {
00567 using namespace edm;
00568 using namespace reco;
00569
00570 bool info = false;
00571 bool alarm = false;
00572
00573 double PI = 3.1415927;
00574
00575 cosmicMuonMode_ = true;
00576 isolatedMuonMode_ = !cosmicMuonMode_;
00577
00578
00579
00580
00581 N_event++;
00582 if (info || debug_) {
00583 std::cout << "----- event " << N_event << " -- tracks "<<N_track<<" ---";
00584 if (cosmicMuonMode_) std::cout << " treated as CosmicMu ";
00585 if (isolatedMuonMode_) std::cout <<" treated as IsolatedMu ";
00586 std::cout << std::endl;
00587 }
00588
00589 Handle<reco::TrackCollection> tracks;
00590 Handle<reco::TrackCollection> muons;
00591 Handle<reco::TrackCollection> gmuons;
00592 Handle<reco::MuonCollection> smuons;
00593
00594 if (collectionIsolated){
00595 iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "TrackerOnly", tracks);
00596 iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "StandAlone", muons);
00597 iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "GlobalMuon", gmuons);
00598 iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "SelectedMuons", smuons);
00599 }
00600 else if (collectionCosmic){
00601 iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "TrackerOnly", tracks);
00602 iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "StandAlone", muons);
00603 iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "GlobalMuon", gmuons);
00604 iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "SelectedMuons", smuons);
00605 }
00606 else{
00607 iEvent.getByLabel(trackTags_,tracks);
00608 iEvent.getByLabel(muonTags_,muons);
00609 iEvent.getByLabel(gmuonTags_,gmuons);
00610 iEvent.getByLabel(smuonTags_,smuons);
00611 }
00612
00613 if (debug_) {
00614 std::cout << " ievBunch " << iEvent.bunchCrossing()
00615 << " runN " << (int)iEvent.run() <<std::endl;
00616 std::cout << " N tracks s/amu gmu selmu "<<tracks->size() << " " <<muons->size()
00617 << " "<<gmuons->size() << " " << smuons->size() << std::endl;
00618 for (MuonCollection::const_iterator itMuon = smuons->begin();
00619 itMuon != smuons->end();
00620 ++itMuon) {
00621 std::cout << " is isolatValid Matches " << itMuon->isIsolationValid()
00622 << " " <<itMuon->isMatchesValid() << std::endl;
00623 }
00624 }
00625
00626 if (isolatedMuonMode_) {
00627 if (tracks->size() != 1) return;
00628 if (muons->size() != 1) return;
00629 if (gmuons->size() != 1) return;
00630 if (smuons->size() != 1) return;
00631 }
00632
00633 if (cosmicMuonMode_){
00634 if (smuons->size() > 2) return;
00635 if (tracks->size() != smuons->size()) return;
00636 if (muons->size() != smuons->size()) return;
00637 if (gmuons->size() != smuons->size()) return;
00638 }
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650 if (watchTrackingGeometry_.check(iSetup) || !trackingGeometry_) {
00651 edm::ESHandle<GlobalTrackingGeometry> trackingGeometry;
00652 iSetup.get<GlobalTrackingGeometryRecord>().get(trackingGeometry);
00653 trackingGeometry_ = &*trackingGeometry;
00654 }
00655
00656 if (watchMagneticFieldRecord_.check(iSetup) || !magneticField_) {
00657 edm::ESHandle<MagneticField> magneticField;
00658 iSetup.get<IdealMagneticFieldRecord>().get(magneticField);
00659 magneticField_ = &*magneticField;
00660 }
00661
00662 if (watchGlobalPositionRcd_.check(iSetup) || !globalPositionRcd_) {
00663 edm::ESHandle<Alignments> globalPositionRcd;
00664 iSetup.get<GlobalPositionRcd>().get(globalPositionRcd);
00665 globalPositionRcd_ = &*globalPositionRcd;
00666 for (std::vector<AlignTransform>::const_iterator i
00667 = globalPositionRcd_->m_align.begin();
00668 i != globalPositionRcd_->m_align.end(); ++i){
00669 if(DetId(DetId::Tracker).rawId() == i->rawId()) iteratorTrackerRcd = i;
00670 if(DetId(DetId::Muon).rawId() == i->rawId()) iteratorMuonRcd = i;
00671 if(DetId(DetId::Ecal).rawId() == i->rawId()) iteratorEcalRcd = i;
00672 if(DetId(DetId::Hcal).rawId() == i->rawId()) iteratorHcalRcd = i;
00673 }
00674 if(true || debug_){
00675 std::cout << "=== iteratorMuonRcd " << iteratorMuonRcd->rawId()<<" ====\n"
00676 << " translation " << iteratorMuonRcd->translation()<<"\n"
00677 << " angles " << iteratorMuonRcd->rotation().eulerAngles() << std::endl;
00678 std::cout << iteratorMuonRcd->rotation() << std::endl;
00679 }
00680 }
00681
00682 ESHandle<Propagator> propagator;
00683 iSetup.get<TrackingComponentsRecord>().get(propagator_, propagator);
00684
00685 SteppingHelixPropagator alongStHePr =
00686 SteppingHelixPropagator(magneticField_, alongMomentum);
00687 SteppingHelixPropagator oppositeStHePr =
00688 SteppingHelixPropagator(magneticField_, oppositeToMomentum);
00689
00690
00691 defaultRKPropagator::Product aprod( magneticField_, alongMomentum, 5.e-5); auto & alongRKPr = aprod.propagator;
00692 defaultRKPropagator::Product oprod( magneticField_, oppositeToMomentum, 5.e-5); auto & oppositeRKPr = oprod.propagator;
00693
00694 float epsilon = 5.;
00695 SmartPropagator alongSmPr =
00696 SmartPropagator(alongRKPr, alongStHePr, magneticField_, alongMomentum, epsilon);
00697 SmartPropagator oppositeSmPr =
00698 SmartPropagator(oppositeRKPr, oppositeStHePr, magneticField_, oppositeToMomentum, epsilon);
00699
00700
00701
00702 for(MuonCollection::const_iterator itMuon = smuons->begin();
00703 itMuon != smuons->end(); ++itMuon) {
00704
00705 if(debug_){
00706 std::cout<<" mu gM is GM Mu SaM tM "<<itMuon->isGlobalMuon()<<" "
00707 <<itMuon->isMuon()<<" "<<itMuon->isStandAloneMuon()
00708 <<" "<<itMuon->isTrackerMuon()<<" "
00709 <<std::endl;
00710 }
00711
00712
00713 TransientTrack muTT(itMuon->outerTrack(), magneticField_, trackingGeometry_);
00714 TrajectoryStateOnSurface innerMuTSOS = muTT.innermostMeasurementState();
00715 TrajectoryStateOnSurface outerMuTSOS = muTT.outermostMeasurementState();
00716
00717
00718 TransientTrack trackTT(itMuon->track(), magneticField_, trackingGeometry_);
00719 TrajectoryStateOnSurface outerTrackTSOS = trackTT.outermostMeasurementState();
00720 TrajectoryStateOnSurface innerTrackTSOS = trackTT.innermostMeasurementState();
00721
00722 GlobalPoint pointTrackIn = innerTrackTSOS.globalPosition();
00723 GlobalPoint pointTrackOut = outerTrackTSOS.globalPosition();
00724 float lenghtTrack = (pointTrackOut-pointTrackIn).mag();
00725 GlobalPoint pointMuonIn = innerMuTSOS.globalPosition();
00726 GlobalPoint pointMuonOut = outerMuTSOS.globalPosition();
00727 float lenghtMuon = (pointMuonOut - pointMuonIn).mag();
00728 GlobalVector momentumTrackOut = outerTrackTSOS.globalMomentum();
00729 GlobalVector momentumTrackIn = innerTrackTSOS.globalMomentum();
00730 float distanceInIn = (pointTrackIn - pointMuonIn).mag();
00731 float distanceInOut = (pointTrackIn - pointMuonOut).mag();
00732 float distanceOutIn = (pointTrackOut - pointMuonIn).mag();
00733 float distanceOutOut = (pointTrackOut - pointMuonOut).mag();
00734
00735 if(debug_){
00736 std::cout<<" pointTrackIn "<<pointTrackIn<<std::endl;
00737 std::cout<<" Out "<<pointTrackOut<<" lenght "<<lenghtTrack<<std::endl;
00738 std::cout<<" point MuonIn "<<pointMuonIn<<std::endl;
00739 std::cout<<" Out "<<pointMuonOut<<" lenght "<<lenghtMuon<<std::endl;
00740 std::cout<<" momeTrackIn Out "<<momentumTrackIn<<" "<<momentumTrackOut<<std::endl;
00741 std::cout<<" doi io ii oo "<<distanceOutIn<<" "<<distanceInOut<<" "
00742 <<distanceInIn<<" "<<distanceOutOut<<std::endl;
00743 }
00744
00745 if(lenghtTrack < 90.) continue;
00746 if(lenghtMuon < 300.) continue;
00747 if(momentumTrackIn.mag() < 15. || momentumTrackOut.mag() < 15.) continue;
00748 if(trackTT.charge() != muTT.charge()) continue;
00749
00750 if(debug_)
00751 std::cout<<" passed lenght momentum cuts"<<std::endl;
00752
00753 GlobalVector GRm, GPm, Nl, Rm, Pm, Rt, Pt, Rt0;
00754 AlgebraicSymMatrix66 Cm, C0, Ce, C1;
00755
00756 TrajectoryStateOnSurface extrapolationT;
00757 int tsosMuonIf = 0;
00758
00759 if( isolatedMuonMode_ ){
00760 const Surface& refSurface = innerMuTSOS.surface();
00761 ReferenceCountingPointer<TangentPlane>
00762 tpMuLocal(refSurface.tangentPlane(innerMuTSOS.localPosition()));
00763 Nl = tpMuLocal->normalVector();
00764 ReferenceCountingPointer<TangentPlane>
00765 tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
00766 GlobalVector Ng = tpMuGlobal->normalVector();
00767 Surface* surf = (Surface*)&refSurface;
00768 const Plane* refPlane = dynamic_cast<Plane*>(surf);
00769 GlobalVector Nlp = refPlane->normalVector();
00770
00771 if(debug_){
00772 std::cout<<" Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
00773 <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
00774 std::cout<<" global "<<Ng.x()<<" "<<Ng.y()<<" "<<Ng.z()<<std::endl;
00775 std::cout<<" lp "<<Nlp.x()<<" "<<Nlp.y()<<" "<<Nlp.z()<<std::endl;
00776
00777
00778
00779 }
00780
00781
00782 extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
00783
00784
00785 if(!extrapolationT.isValid()){
00786 if(false & alarm)
00787 std::cout<<" !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
00788
00789 <<std::endl;
00790 continue;
00791 }
00792 tsosMuonIf = 1;
00793 Rt = GlobalVector((extrapolationT.globalPosition()).x(),
00794 (extrapolationT.globalPosition()).y(),
00795 (extrapolationT.globalPosition()).z());
00796
00797 Pt = extrapolationT.globalMomentum();
00798
00799 GRm = GlobalVector((innerMuTSOS.globalPosition()).x(),
00800 (innerMuTSOS.globalPosition()).y(),
00801 (innerMuTSOS.globalPosition()).z());
00802 GPm = innerMuTSOS.globalMomentum();
00803 Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
00804 (outerTrackTSOS.globalPosition()).y(),
00805 (outerTrackTSOS.globalPosition()).z());
00806 Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() +
00807 innerMuTSOS.cartesianError().matrix());
00808 C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
00809 Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
00810 C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
00811
00812 }
00813
00814
00815 if( cosmicMuonMode_ ){
00816
00817 if((distanceOutIn <= distanceInOut) & (distanceOutIn <= distanceInIn) &
00818 (distanceOutIn <= distanceOutOut)){
00819 if(debug_) std::cout<<" ----- Out - In -----"<<std::endl;
00820
00821 const Surface& refSurface = innerMuTSOS.surface();
00822 ReferenceCountingPointer<TangentPlane>
00823 tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
00824 Nl = tpMuGlobal->normalVector();
00825
00826
00827 extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
00828
00829
00830 if(!extrapolationT.isValid()){
00831 if(false & alarm)
00832 std::cout<<" !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
00833
00834 <<std::endl;
00835 continue;
00836 }
00837 if(debug_)
00838 std::cout<<" extrapolationT.isValid "<<extrapolationT.isValid()<<std::endl;
00839
00840 tsosMuonIf = 1;
00841 Rt = GlobalVector((extrapolationT.globalPosition()).x(),
00842 (extrapolationT.globalPosition()).y(),
00843 (extrapolationT.globalPosition()).z());
00844
00845 Pt = extrapolationT.globalMomentum();
00846
00847 GRm = GlobalVector((innerMuTSOS.globalPosition()).x(),
00848 (innerMuTSOS.globalPosition()).y(),
00849 (innerMuTSOS.globalPosition()).z());
00850 GPm = innerMuTSOS.globalMomentum();
00851 Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
00852 (outerTrackTSOS.globalPosition()).y(),
00853 (outerTrackTSOS.globalPosition()).z());
00854 Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() +
00855 innerMuTSOS.cartesianError().matrix());
00856 C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
00857 Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
00858 C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
00859
00860 if(false & debug_){
00861
00862 PropagationDirectionChooser Chooser;
00863 std::cout<<" propDirCh "
00864 <<Chooser.operator()(*outerTrackTSOS.freeState(), refSurface)
00865 <<" Ch == along "<<(alongMomentum ==
00866 Chooser.operator()(*outerTrackTSOS.freeState(), refSurface))
00867 <<std::endl;
00868 std::cout<<" --- Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
00869 <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
00870 }
00871 }
00872
00873
00874 else if((distanceInOut <= distanceInIn) & (distanceInOut <= distanceOutIn) &
00875 (distanceInOut <= distanceOutOut)){
00876 if(debug_) std::cout<<" ----- In - Out -----"<<std::endl;
00877
00878 const Surface& refSurface = outerMuTSOS.surface();
00879 ReferenceCountingPointer<TangentPlane>
00880 tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
00881 Nl = tpMuGlobal->normalVector();
00882
00883
00884 extrapolationT = oppositeSmPr.propagate(innerTrackTSOS, refSurface);
00885
00886 if(!extrapolationT.isValid()){
00887 if(false & alarm)
00888 std::cout<<" !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
00889 <<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
00890 <<std::endl;
00891 continue;
00892 }
00893 if(debug_)
00894 std::cout<<" extrapolationT.isValid "<<extrapolationT.isValid()<<std::endl;
00895
00896 tsosMuonIf = 2;
00897 Rt = GlobalVector((extrapolationT.globalPosition()).x(),
00898 (extrapolationT.globalPosition()).y(),
00899 (extrapolationT.globalPosition()).z());
00900
00901 Pt = extrapolationT.globalMomentum();
00902
00903 GRm = GlobalVector((outerMuTSOS.globalPosition()).x(),
00904 (outerMuTSOS.globalPosition()).y(),
00905 (outerMuTSOS.globalPosition()).z());
00906 GPm = outerMuTSOS.globalMomentum();
00907 Rt0 = GlobalVector((innerTrackTSOS.globalPosition()).x(),
00908 (innerTrackTSOS.globalPosition()).y(),
00909 (innerTrackTSOS.globalPosition()).z());
00910 Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() +
00911 outerMuTSOS.cartesianError().matrix());
00912 C0 = AlgebraicSymMatrix66(innerTrackTSOS.cartesianError().matrix());
00913 Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
00914 C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
00915
00916 if(false & debug_){
00917
00918 PropagationDirectionChooser Chooser;
00919 std::cout<<" propDirCh "
00920 <<Chooser.operator()(*innerTrackTSOS.freeState(), refSurface)
00921 <<" Ch == oppisite "<<(oppositeToMomentum ==
00922 Chooser.operator()(*innerTrackTSOS.freeState(), refSurface))
00923 <<std::endl;
00924 std::cout<<" --- Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
00925 <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
00926 }
00927 }
00928
00929 else if((distanceOutOut <= distanceInOut) & (distanceOutOut <= distanceInIn) &
00930 (distanceOutOut <= distanceOutIn)){
00931 if(debug_) std::cout<<" ----- Out - Out -----"<<std::endl;
00932
00933
00934 continue;
00935
00936 const Surface& refSurface = outerMuTSOS.surface();
00937 ReferenceCountingPointer<TangentPlane>
00938 tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
00939 Nl = tpMuGlobal->normalVector();
00940
00941
00942 extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
00943
00944 if(!extrapolationT.isValid()){
00945 if(alarm)
00946 std::cout<<" !!!!!!!!! Catastrophe Out-Out extrapolationT.isValid "
00947 <<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
00948 <<std::endl;
00949 continue;
00950 }
00951 if(debug_)
00952 std::cout<<" extrapolationT.isValid "<<extrapolationT.isValid()<<std::endl;
00953
00954 tsosMuonIf = 2;
00955 Rt = GlobalVector((extrapolationT.globalPosition()).x(),
00956 (extrapolationT.globalPosition()).y(),
00957 (extrapolationT.globalPosition()).z());
00958
00959 Pt = extrapolationT.globalMomentum();
00960
00961 GRm = GlobalVector((outerMuTSOS.globalPosition()).x(),
00962 (outerMuTSOS.globalPosition()).y(),
00963 (outerMuTSOS.globalPosition()).z());
00964 GPm = outerMuTSOS.globalMomentum();
00965 Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
00966 (outerTrackTSOS.globalPosition()).y(),
00967 (outerTrackTSOS.globalPosition()).z());
00968 Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() +
00969 outerMuTSOS.cartesianError().matrix());
00970 C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
00971 Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
00972 C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
00973
00974 if(debug_){
00975 PropagationDirectionChooser Chooser;
00976 std::cout<<" propDirCh "
00977 <<Chooser.operator()(*outerTrackTSOS.freeState(), refSurface)
00978 <<" Ch == along "<<(alongMomentum ==
00979 Chooser.operator()(*outerTrackTSOS.freeState(), refSurface))
00980 <<std::endl;
00981 std::cout<<" --- Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
00982 <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
00983 std::cout<<" Nornal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
00984 <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
00985 }
00986 }
00987 else{
00988 if(alarm)
00989 std::cout<<" ------- !!!!!!!!!!!!!!! No proper Out - In \a\a\a\a\a\a\a"<<std::endl;
00990 continue;
00991 }
00992
00993 }
00994
00995 if(tsosMuonIf == 0) {if(info) {std::cout<<"No tsosMuon !!!!!!"<<std::endl; continue;}}
00996 TrajectoryStateOnSurface tsosMuon;
00997 if(tsosMuonIf == 1) tsosMuon = muTT.innermostMeasurementState();
00998 else tsosMuon = muTT.outermostMeasurementState();
00999
01000
01001 AlgebraicVector4 LPRm;
01002 GlobalTrackerMuonAlignment::misalignMuonL
01003 (GRm, GPm, Nl, Rt, Rm, Pm, LPRm, extrapolationT, tsosMuon);
01004
01005 GlobalVector resR = Rm - Rt;
01006 GlobalVector resP0 = Pm - Pt;
01007 GlobalVector resP = Pm / Pm.mag() - Pt / Pt.mag();
01008 float RelMomResidual = (Pm.mag() - Pt.mag()) / (Pt.mag() + 1.e-6);;
01009
01010 AlgebraicVector6 Vm;
01011 Vm(0) = resR.x(); Vm(1) = resR.y(); Vm(2) = resR.z();
01012 Vm(3) = resP0.x(); Vm(4) = resP0.y(); Vm(5) = resP0.z();
01013 float Rmuon = Rm.perp();
01014 float Zmuon = Rm.z();
01015 float alfa_x = atan2(Nl.x(),Nl.y())*57.29578;
01016
01017 if(debug_){
01018 std::cout<<" Nx Ny Nz alfa_x "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "<<alfa_x<<std::endl;
01019
01020
01021 }
01022
01023 double chi_d = 0;
01024 for(int i=0; i<=5; i++) chi_d += Vm(i)*Vm(i)/Cm(i,i);
01025
01026 AlgebraicVector5 Vml(tsosMuon.localParameters().vector()
01027 - extrapolationT.localParameters().vector());
01028 AlgebraicSymMatrix55 m(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
01029 AlgebraicSymMatrix55 Cml(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
01030 bool ierrLoc = !m.Invert();
01031 if (ierrLoc && debug_ && info) {
01032 std::cout<< " ==== Error inverting Local covariance matrix ==== "<<std::cout;
01033 continue;}
01034 double chi_Loc = ROOT::Math::Similarity(Vml,m);
01035 if(debug_)
01036 std::cout<<" chi_Loc px/pz/err "<<chi_Loc<<" "<<Vml(1)/sqrt(Cml(1,1))<<std::endl;
01037
01038 if(Pt.mag() < 15.) continue;
01039 if(Pm.mag() < 5.) continue;
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079
01080
01081
01082
01083 if(Rmuon < 120. || Rmuon > 450.) continue;
01084 if(Zmuon < -720. || Zmuon > 720.) continue;
01085 MuSelect = " Barrel+EndCaps";
01086
01087
01088 if(debug_)
01089 std::cout<<" .............. passed all cuts"<<std::endl;
01090
01091 N_track++;
01092
01093
01094 GlobalTrackerMuonAlignment::gradientGlobalAlg(Rt, Pt, Rm, Nl, Cm);
01095 GlobalTrackerMuonAlignment::gradientGlobal(Rt, Pt, Rm, Pm, Nl, Cm);
01096
01097 CLHEP::HepSymMatrix covLoc(4,0);
01098 for(int i=1; i<=4; i++)
01099 for(int j=1; j<=i; j++){
01100 covLoc(i,j) = (tsosMuon.localError().matrix()
01101 + extrapolationT.localError().matrix())(i,j);
01102
01103 }
01104
01105 const Surface& refSurface = tsosMuon.surface();
01106 CLHEP::HepMatrix rotLoc (3,3,0);
01107 rotLoc(1,1) = refSurface.rotation().xx();
01108 rotLoc(1,2) = refSurface.rotation().xy();
01109 rotLoc(1,3) = refSurface.rotation().xz();
01110
01111 rotLoc(2,1) = refSurface.rotation().yx();
01112 rotLoc(2,2) = refSurface.rotation().yy();
01113 rotLoc(2,3) = refSurface.rotation().yz();
01114
01115 rotLoc(3,1) = refSurface.rotation().zx();
01116 rotLoc(3,2) = refSurface.rotation().zy();
01117 rotLoc(3,3) = refSurface.rotation().zz();
01118
01119 CLHEP::HepVector posLoc(3);
01120 posLoc(1) = refSurface.position().x();
01121 posLoc(2) = refSurface.position().y();
01122 posLoc(3) = refSurface.position().z();
01123
01124 GlobalTrackerMuonAlignment::gradientLocal(Rt, Pt, Rm, Pm, Nl, covLoc, rotLoc, posLoc, LPRm);
01125
01126 if(debug_){
01127 std::cout<<" Norm "<<Nl<<std::endl;
01128 std::cout<<" posLoc "<<posLoc.T()<<std::endl;
01129 std::cout<<" rotLoc "<<rotLoc<<std::endl;
01130 }
01131
01132
01133 histo->Fill(itMuon->track()->pt());
01134
01135
01136 histo2->Fill(Pt.mag());
01137 histo3->Fill((PI/2.-itMuon->track()->outerTheta()));
01138 histo4->Fill(itMuon->track()->phi());
01139 histo5->Fill(Rmuon);
01140 histo6->Fill(Zmuon);
01141 histo7->Fill(RelMomResidual);
01142
01143 histo8->Fill(chi_d);
01144
01145 histo101->Fill(Zmuon, Rmuon);
01146 histo101->Fill(Rt0.z(), Rt0.perp());
01147 histo102->Fill(Rt0.x(), Rt0.y());
01148 histo102->Fill(Rm.x(), Rm.y());
01149
01150 histo11->Fill(resR.mag());
01151 if(fabs(Nl.x()) < 0.98) histo12->Fill(resR.x());
01152 if(fabs(Nl.y()) < 0.98) histo13->Fill(resR.y());
01153 if(fabs(Nl.z()) < 0.98) histo14->Fill(resR.z());
01154 histo15->Fill(resP.x());
01155 histo16->Fill(resP.y());
01156 histo17->Fill(resP.z());
01157
01158 if((fabs(Nl.x()) < 0.98) && (fabs(Nl.y()) < 0.98) &&(fabs(Nl.z()) < 0.98))
01159 {
01160 histo18->Fill(sqrt(C0(0,0)));
01161 histo19->Fill(sqrt(C1(0,0)));
01162 histo20->Fill(sqrt(C1(0,0)+Ce(0,0)));
01163 }
01164 if(fabs(Nl.x()) < 0.98) histo21->Fill(Vm(0)/sqrt(Cm(0,0)));
01165 if(fabs(Nl.y()) < 0.98) histo22->Fill(Vm(1)/sqrt(Cm(1,1)));
01166 if(fabs(Nl.z()) < 0.98) histo23->Fill(Vm(2)/sqrt(Cm(2,2)));
01167 histo24->Fill(Vm(3)/sqrt(C1(3,3)+Ce(3,3)));
01168 histo25->Fill(Vm(4)/sqrt(C1(4,4)+Ce(4,4)));
01169 histo26->Fill(Vm(5)/sqrt(C1(5,5)+Ce(5,5)));
01170 histo27->Fill(Nl.x());
01171 histo28->Fill(Nl.y());
01172 histo29->Fill(lenghtTrack);
01173 histo30->Fill(lenghtMuon);
01174 histo31->Fill(chi_Loc);
01175 histo32->Fill(Vml(1)/sqrt(Cml(1,1)));
01176 histo33->Fill(Vml(2)/sqrt(Cml(2,2)));
01177 histo34->Fill(Vml(3)/sqrt(Cml(3,3)));
01178 histo35->Fill(Vml(4)/sqrt(Cml(4,4)));
01179
01180 if (debug_) {
01181
01182 std::cout<<" diag 0[ "<<C0(0,0)<<" "<<C0(1,1)<<" "<<C0(2,2)<<" "<<C0(3,3)<<" "
01183 <<C0(4,4)<<" "<<C0(5,5)<<" ]"<<std::endl;
01184 std::cout<<" diag e[ "<<Ce(0,0)<<" "<<Ce(1,1)<<" "<<Ce(2,2)<<" "<<Ce(3,3)<<" "
01185 <<Ce(4,4)<<" "<<Ce(5,5)<<" ]"<<std::endl;
01186 std::cout<<" diag 1[ "<<C1(0,0)<<" "<<C1(1,1)<<" "<<C1(2,2)<<" "<<C1(3,3)<<" "
01187 <<C1(4,4)<<" "<<C1(5,5)<<" ]"<<std::endl;
01188 std::cout<<" Rm "<<Rm.x()<<" "<<Rm.y()<<" "<<Rm.z()
01189 <<" Pm "<<Pm.x()<<" "<<Pm.y()<<" "<<Pm.z()<<std::endl;
01190 std::cout<<" Rt "<<Rt.x()<<" "<<Rt.y()<<" "<<Rt.z()
01191 <<" Pt "<<Pt.x()<<" "<<Pt.y()<<" "<<Pt.z()<<std::endl;
01192 std::cout<<" Nl*(Rm-Rt) "<<Nl.dot(Rm-Rt)<<std::endl;
01193 std::cout<<" resR "<<resR.x()<<" "<<resR.y()<<" "<<resR.z()
01194 <<" resP "<<resP.x()<<" "<<resP.y()<<" "<<resP.z()<<std::endl;
01195 std::cout<<" Rm-t "<<(Rm-Rt).x()<<" "<<(Rm-Rt).y()<<" "<<(Rm-Rt).z()
01196 <<" Pm-t "<<(Pm-Pt).x()<<" "<<(Pm-Pt).y()<<" "<<(Pm-Pt).z()<<std::endl;
01197 std::cout<<" Vm "<<Vm<<std::endl;
01198 std::cout<<" +- "<<sqrt(Cm(0,0))<<" "<<sqrt(Cm(1,1))<<" "<<sqrt(Cm(2,2))<<" "
01199 <<sqrt(Cm(3,3))<<" "<<sqrt(Cm(4,4))<<" "<<sqrt(Cm(5,5))<<std::endl;
01200 std::cout<<" Pmuon Ptrack dP/Ptrack "<<itMuon->outerTrack()->p()<<" "
01201 <<itMuon->track()->outerP()<<" "
01202 <<(itMuon->outerTrack()->p() -
01203 itMuon->track()->outerP())/itMuon->track()->outerP()<<std::endl;
01204 std::cout<<" cov matrix "<<std::endl;
01205 std::cout<<Cm<<std::endl;
01206 std::cout<<" diag [ "<<Cm(0,0)<<" "<<Cm(1,1)<<" "<<Cm(2,2)<<" "<<Cm(3,3)<<" "
01207 <<Cm(4,4)<<" "<<Cm(5,5)<<" ]"<<std::endl;
01208
01209 AlgebraicSymMatrix66 Ro;
01210 double Diag[6];
01211 for(int i=0; i<=5; i++) Diag[i] = sqrt(Cm(i,i));
01212 for(int i=0; i<=5; i++)
01213 for(int j=0; j<=5; j++)
01214 Ro(i,j) = Cm(i,j)/Diag[i]/Diag[j];
01215 std::cout<<" correlation matrix "<<std::endl;
01216 std::cout<<Ro<<std::endl;
01217
01218 AlgebraicSymMatrix66 CmI;
01219 for(int i=0; i<=5; i++)
01220 for(int j=0; j<=5; j++)
01221 CmI(i,j) = Cm(i,j);
01222
01223 bool ierr = !CmI.Invert();
01224 if( ierr ) { if(alarm || debug_)
01225 std::cout<<" Error inverse covariance matrix !!!!!!!!!!!"<<std::endl;
01226 continue;
01227 }
01228 std::cout<<" inverse cov matrix "<<std::endl;
01229 std::cout<<Cm<<std::endl;
01230
01231 double chi = ROOT::Math::Similarity(Vm, CmI);
01232 std::cout<<" chi chi_d "<<chi<<" "<<chi_d<<std::endl;
01233 }
01234
01235 }
01236
01237 }
01238
01239
01240
01241 void GlobalTrackerMuonAlignment::analyzeTrackTrajectory
01242 (const edm::Event& iEvent, const edm::EventSetup& iSetup)
01243 {
01244 using namespace edm;
01245 using namespace reco;
01246
01247 bool info = false;
01248 bool alarm = false;
01249
01250 double PI = 3.1415927;
01251
01252
01253 cosmicMuonMode_ = true;
01254
01255
01256 isolatedMuonMode_ = !cosmicMuonMode_;
01257
01258 N_event++;
01259 if (info || debug_) {
01260 std::cout << "----- event " << N_event << " -- tracks "<<N_track<<" ---";
01261 if (cosmicMuonMode_) std::cout << " treated as CosmicMu ";
01262 if (isolatedMuonMode_) std::cout <<" treated as IsolatedMu ";
01263 std::cout << std::endl;
01264 }
01265
01266 Handle<reco::TrackCollection> tracks;
01267 Handle<reco::TrackCollection> muons;
01268 Handle<reco::TrackCollection> gmuons;
01269 Handle<reco::MuonCollection> smuons;
01270
01271 if (collectionIsolated){
01272 iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "TrackerOnly", tracks);
01273 iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "StandAlone", muons);
01274 iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "GlobalMuon", gmuons);
01275 iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "SelectedMuons", smuons);
01276 }
01277 else if (collectionCosmic){
01278 iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "TrackerOnly", tracks);
01279 iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "StandAlone", muons);
01280 iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "GlobalMuon", gmuons);
01281 iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "SelectedMuons", smuons);
01282 }
01283 else{
01284 iEvent.getByLabel(trackTags_,tracks);
01285 iEvent.getByLabel(muonTags_,muons);
01286 iEvent.getByLabel(gmuonTags_,gmuons);
01287 iEvent.getByLabel(smuonTags_,smuons);
01288 }
01289
01290 if (debug_) {
01291 std::cout << " ievBunch " << iEvent.bunchCrossing()
01292 << " runN " << (int)iEvent.run() <<std::endl;
01293 std::cout << " N tracks s/amu gmu selmu "<<tracks->size() << " " <<muons->size()
01294 << " "<<gmuons->size() << " " << smuons->size() << std::endl;
01295 for (MuonCollection::const_iterator itMuon = smuons->begin();
01296 itMuon != smuons->end();
01297 ++itMuon) {
01298 std::cout << " is isolatValid Matches " << itMuon->isIsolationValid()
01299 << " " <<itMuon->isMatchesValid() << std::endl;
01300 }
01301 }
01302
01303 if (isolatedMuonMode_) {
01304 if (tracks->size() != 1) return;
01305 if (muons->size() != 1) return;
01306 if (gmuons->size() != 1) return;
01307 if (smuons->size() != 1) return;
01308 }
01309
01310 if (cosmicMuonMode_){
01311 if (smuons->size() > 2) return;
01312 if (tracks->size() != smuons->size()) return;
01313 if (muons->size() != smuons->size()) return;
01314 if (gmuons->size() != smuons->size()) return;
01315 }
01316
01317
01318
01319
01320
01321
01322
01323
01324
01325
01326
01327 if (watchTrackingGeometry_.check(iSetup) || !trackingGeometry_) {
01328 edm::ESHandle<GlobalTrackingGeometry> trackingGeometry;
01329 iSetup.get<GlobalTrackingGeometryRecord>().get(trackingGeometry);
01330 trackingGeometry_ = &*trackingGeometry;
01331 theTrackingGeometry = trackingGeometry;
01332 }
01333
01334 if (watchMagneticFieldRecord_.check(iSetup) || !magneticField_) {
01335 edm::ESHandle<MagneticField> magneticField;
01336 iSetup.get<IdealMagneticFieldRecord>().get(magneticField);
01337 magneticField_ = &*magneticField;
01338 }
01339
01340 if (watchGlobalPositionRcd_.check(iSetup) || !globalPositionRcd_) {
01341 edm::ESHandle<Alignments> globalPositionRcd;
01342 iSetup.get<GlobalPositionRcd>().get(globalPositionRcd);
01343 globalPositionRcd_ = &*globalPositionRcd;
01344 for (std::vector<AlignTransform>::const_iterator i
01345 = globalPositionRcd_->m_align.begin();
01346 i != globalPositionRcd_->m_align.end(); ++i){
01347 if(DetId(DetId::Tracker).rawId() == i->rawId()) iteratorTrackerRcd = i;
01348 if(DetId(DetId::Muon).rawId() == i->rawId()) iteratorMuonRcd = i;
01349 if(DetId(DetId::Ecal).rawId() == i->rawId()) iteratorEcalRcd = i;
01350 if(DetId(DetId::Hcal).rawId() == i->rawId()) iteratorHcalRcd = i;
01351 }
01352 if(debug_){
01353 std::cout << "=== iteratorTrackerRcd " << iteratorTrackerRcd->rawId()<<" ====\n"
01354 << " translation " << iteratorTrackerRcd->translation()<<"\n"
01355 << " angles " << iteratorTrackerRcd->rotation().eulerAngles() << std::endl;
01356 std::cout << iteratorTrackerRcd->rotation() << std::endl;
01357 std::cout << "=== iteratorMuonRcd " << iteratorMuonRcd->rawId()<<" ====\n"
01358 << " translation " << iteratorMuonRcd->translation()<<"\n"
01359 << " angles " << iteratorMuonRcd->rotation().eulerAngles() << std::endl;
01360 std::cout << iteratorMuonRcd->rotation() << std::endl;
01361 }
01362 }
01363
01364 ESHandle<Propagator> propagator;
01365 iSetup.get<TrackingComponentsRecord>().get(propagator_, propagator);
01366
01367 SteppingHelixPropagator alongStHePr =
01368 SteppingHelixPropagator(magneticField_, alongMomentum);
01369 SteppingHelixPropagator oppositeStHePr =
01370 SteppingHelixPropagator(magneticField_, oppositeToMomentum);
01371
01372
01373 defaultRKPropagator::Product aprod( magneticField_, alongMomentum, 5.e-5); auto & alongRKPr = aprod.propagator;
01374 defaultRKPropagator::Product oprod( magneticField_, oppositeToMomentum, 5.e-5); auto & oppositeRKPr = oprod.propagator;
01375
01376 float epsilon = 5.;
01377 SmartPropagator alongSmPr =
01378 SmartPropagator(alongRKPr, alongStHePr, magneticField_, alongMomentum, epsilon);
01379 SmartPropagator oppositeSmPr =
01380 SmartPropagator(oppositeRKPr, oppositeStHePr, magneticField_, oppositeToMomentum, epsilon);
01381
01382 if(defineFitter){
01383 if(debug_)
01384 std::cout<<" ............... DEFINE FITTER ..................."<<std::endl;
01385 KFUpdator* theUpdator = new KFUpdator();
01386
01387 Chi2MeasurementEstimator* theEstimator = new Chi2MeasurementEstimator(100000,100000);
01388 theFitter = new KFTrajectoryFitter(alongSmPr,
01389 *theUpdator,
01390 *theEstimator);
01391 theSmoother = new KFTrajectorySmoother(alongSmPr,
01392 *theUpdator,
01393 *theEstimator);
01394 theFitterOp = new KFTrajectoryFitter(oppositeSmPr,
01395 *theUpdator,
01396 *theEstimator);
01397 theSmootherOp = new KFTrajectorySmoother(oppositeSmPr,
01398 *theUpdator,
01399 *theEstimator);
01400
01401 edm::ESHandle<TransientTrackingRecHitBuilder> builder;
01402 iSetup.get<TransientRecHitRecord>().get("WithTrackAngle",builder);
01403 this->TTRHBuilder = &(*builder);
01404 this->MuRHBuilder = new MuonTransientTrackingRecHitBuilder(theTrackingGeometry);
01405 if(debug_){
01406 std::cout<<" theTrackingGeometry.isValid() "<<theTrackingGeometry.isValid()<<std::endl;
01407 std::cout<< "get also the MuonTransientTrackingRecHitBuilder" << "\n";
01408 std::cout<< "get also the TransientTrackingRecHitBuilder" << "\n";
01409 }
01410 defineFitter = false;
01411 }
01412
01413
01414
01415 for(MuonCollection::const_iterator itMuon = smuons->begin();
01416 itMuon != smuons->end(); ++itMuon) {
01417
01418 if(debug_){
01419 std::cout<<" mu gM is GM Mu SaM tM "<<itMuon->isGlobalMuon()<<" "
01420 <<itMuon->isMuon()<<" "<<itMuon->isStandAloneMuon()
01421 <<" "<<itMuon->isTrackerMuon()<<" "
01422 <<std::endl;
01423 }
01424
01425
01426 TransientTrack muTT(itMuon->outerTrack(), magneticField_, trackingGeometry_);
01427 TrajectoryStateOnSurface innerMuTSOS = muTT.innermostMeasurementState();
01428 TrajectoryStateOnSurface outerMuTSOS = muTT.outermostMeasurementState();
01429
01430
01431 TransientTrack trackTT(itMuon->track(), magneticField_, trackingGeometry_);
01432 TrajectoryStateOnSurface outerTrackTSOS = trackTT.outermostMeasurementState();
01433 TrajectoryStateOnSurface innerTrackTSOS = trackTT.innermostMeasurementState();
01434
01435 GlobalPoint pointTrackIn = innerTrackTSOS.globalPosition();
01436 GlobalPoint pointTrackOut = outerTrackTSOS.globalPosition();
01437 float lenghtTrack = (pointTrackOut-pointTrackIn).mag();
01438 GlobalPoint pointMuonIn = innerMuTSOS.globalPosition();
01439 GlobalPoint pointMuonOut = outerMuTSOS.globalPosition();
01440 float lenghtMuon = (pointMuonOut - pointMuonIn).mag();
01441 GlobalVector momentumTrackOut = outerTrackTSOS.globalMomentum();
01442 GlobalVector momentumTrackIn = innerTrackTSOS.globalMomentum();
01443 float distanceInIn = (pointTrackIn - pointMuonIn).mag();
01444 float distanceInOut = (pointTrackIn - pointMuonOut).mag();
01445 float distanceOutIn = (pointTrackOut - pointMuonIn).mag();
01446 float distanceOutOut = (pointTrackOut - pointMuonOut).mag();
01447
01448 if(debug_){
01449 std::cout<<" pointTrackIn "<<pointTrackIn<<std::endl;
01450 std::cout<<" Out "<<pointTrackOut<<" lenght "<<lenghtTrack<<std::endl;
01451 std::cout<<" point MuonIn "<<pointMuonIn<<std::endl;
01452 std::cout<<" Out "<<pointMuonOut<<" lenght "<<lenghtMuon<<std::endl;
01453 std::cout<<" momeTrackIn Out "<<momentumTrackIn<<" "<<momentumTrackOut<<std::endl;
01454 std::cout<<" doi io ii oo "<<distanceOutIn<<" "<<distanceInOut<<" "
01455 <<distanceInIn<<" "<<distanceOutOut<<std::endl;
01456 }
01457
01458 if(lenghtTrack < 90.) continue;
01459 if(lenghtMuon < 300.) continue;
01460 if(innerMuTSOS.globalMomentum().mag() < 5. || outerMuTSOS.globalMomentum().mag() < 5.) continue;
01461 if(momentumTrackIn.mag() < 15. || momentumTrackOut.mag() < 15.) continue;
01462 if(trackTT.charge() != muTT.charge()) continue;
01463
01464 if(debug_)
01465 std::cout<<" passed lenght momentum cuts"<<std::endl;
01466
01467 GlobalVector GRm, GPm, Nl, Rm, Pm, Rt, Pt, Rt0;
01468 AlgebraicSymMatrix66 Cm, C0, Ce, C1;
01469
01470
01471 TrajectoryStateOnSurface extrapolationT;
01472 int tsosMuonIf = 0;
01473
01474 TrajectoryStateOnSurface muonFittedTSOS;
01475 TrajectoryStateOnSurface trackFittedTSOS;
01476
01477 if( isolatedMuonMode_ ){
01478 if(debug_) std::cout<<" ------ Isolated (out-in) ----- "<<std::endl;
01479 const Surface& refSurface = innerMuTSOS.surface();
01480 ReferenceCountingPointer<TangentPlane>
01481 tpMuLocal(refSurface.tangentPlane(innerMuTSOS.localPosition()));
01482 Nl = tpMuLocal->normalVector();
01483 ReferenceCountingPointer<TangentPlane>
01484 tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
01485 GlobalVector Ng = tpMuGlobal->normalVector();
01486 Surface* surf = (Surface*)&refSurface;
01487 const Plane* refPlane = dynamic_cast<Plane*>(surf);
01488 GlobalVector Nlp = refPlane->normalVector();
01489
01490 if(debug_){
01491 std::cout<<" Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
01492 <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
01493 std::cout<<" global "<<Ng.x()<<" "<<Ng.y()<<" "<<Ng.z()<<std::endl;
01494 std::cout<<" lp "<<Nlp.x()<<" "<<Nlp.y()<<" "<<Nlp.z()<<std::endl;
01495
01496
01497
01498 }
01499
01500
01501
01502
01503 if(!refitTrack_)
01504 extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
01505 else{
01506 GlobalTrackerMuonAlignment::trackFitter(itMuon->track(),
01507 trackTT,alongMomentum,trackFittedTSOS);
01508 if(trackFittedTSOS.isValid())
01509 extrapolationT = alongSmPr.propagate(trackFittedTSOS, refSurface);
01510 }
01511
01512 if(!extrapolationT.isValid()){
01513 if(false & alarm)
01514 std::cout<<" !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
01515
01516 <<std::endl;
01517 continue;
01518 }
01519 tsosMuonIf = 1;
01520 Rt = GlobalVector((extrapolationT.globalPosition()).x(),
01521 (extrapolationT.globalPosition()).y(),
01522 (extrapolationT.globalPosition()).z());
01523
01524 Pt = extrapolationT.globalMomentum();
01525
01526 GRm = GlobalVector((innerMuTSOS.globalPosition()).x(),
01527 (innerMuTSOS.globalPosition()).y(),
01528 (innerMuTSOS.globalPosition()).z());
01529 GPm = innerMuTSOS.globalMomentum();
01530
01531 Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
01532 (outerTrackTSOS.globalPosition()).y(),
01533 (outerTrackTSOS.globalPosition()).z());
01534 Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() +
01535 innerMuTSOS.cartesianError().matrix());
01536 C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
01537 Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
01538 C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
01539
01540 if(refitMuon_)
01541 GlobalTrackerMuonAlignment::muonFitter(itMuon->outerTrack(),muTT,oppositeToMomentum,muonFittedTSOS);
01542
01543 }
01544
01545
01546 if( cosmicMuonMode_ ){
01547
01548 if((distanceOutIn <= distanceInOut) & (distanceOutIn <= distanceInIn) &
01549 (distanceOutIn <= distanceOutOut)){
01550 if(debug_) std::cout<<" ----- Out - In -----"<<std::endl;
01551
01552 const Surface& refSurface = innerMuTSOS.surface();
01553 ReferenceCountingPointer<TangentPlane>
01554 tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
01555 Nl = tpMuGlobal->normalVector();
01556
01557
01558
01559 if(!refitTrack_)
01560 extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
01561 else{
01562 GlobalTrackerMuonAlignment::trackFitter(itMuon->track(),
01563 trackTT,alongMomentum,trackFittedTSOS);
01564 if(trackFittedTSOS.isValid())
01565 extrapolationT = alongSmPr.propagate(trackFittedTSOS, refSurface);
01566 }
01567
01568 if(!extrapolationT.isValid()){
01569 if(false & alarm)
01570 std::cout<<" !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
01571
01572 <<std::endl;
01573 continue;
01574 }
01575 if(debug_)
01576 std::cout<<" extrapolationT.isValid "<<extrapolationT.isValid()<<std::endl;
01577
01578 tsosMuonIf = 1;
01579 Rt = GlobalVector((extrapolationT.globalPosition()).x(),
01580 (extrapolationT.globalPosition()).y(),
01581 (extrapolationT.globalPosition()).z());
01582
01583 Pt = extrapolationT.globalMomentum();
01584
01585 GRm = GlobalVector((innerMuTSOS.globalPosition()).x(),
01586 (innerMuTSOS.globalPosition()).y(),
01587 (innerMuTSOS.globalPosition()).z());
01588 GPm = innerMuTSOS.globalMomentum();
01589 Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
01590 (outerTrackTSOS.globalPosition()).y(),
01591 (outerTrackTSOS.globalPosition()).z());
01592 Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() +
01593 innerMuTSOS.cartesianError().matrix());
01594 C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
01595 Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
01596 C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
01597
01598 if(refitMuon_)
01599 GlobalTrackerMuonAlignment::muonFitter(itMuon->outerTrack(),muTT,oppositeToMomentum,muonFittedTSOS);
01600
01601 if(false & debug_){
01602
01603 PropagationDirectionChooser Chooser;
01604 std::cout<<" propDirCh "
01605 <<Chooser.operator()(*outerTrackTSOS.freeState(), refSurface)
01606 <<" Ch == along "<<(alongMomentum ==
01607 Chooser.operator()(*outerTrackTSOS.freeState(), refSurface))
01608 <<std::endl;
01609 std::cout<<" --- Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
01610 <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
01611 }
01612 }
01613
01614
01615 else if((distanceInOut <= distanceInIn) & (distanceInOut <= distanceOutIn) &
01616 (distanceInOut <= distanceOutOut)){
01617 if(debug_) std::cout<<" ----- In - Out -----"<<std::endl;
01618
01619 const Surface& refSurface = outerMuTSOS.surface();
01620 ReferenceCountingPointer<TangentPlane>
01621 tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
01622 Nl = tpMuGlobal->normalVector();
01623
01624
01625
01626 if(!refitTrack_)
01627 extrapolationT = oppositeSmPr.propagate(innerTrackTSOS, refSurface);
01628 else{
01629 GlobalTrackerMuonAlignment::trackFitter(itMuon->track(),
01630 trackTT,oppositeToMomentum,trackFittedTSOS);
01631 if(trackFittedTSOS.isValid())
01632 extrapolationT = oppositeSmPr.propagate(trackFittedTSOS, refSurface);
01633 }
01634
01635 if(!extrapolationT.isValid()){
01636 if(false & alarm)
01637 std::cout<<" !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
01638 <<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
01639 <<std::endl;
01640 continue;
01641 }
01642 if(debug_)
01643 std::cout<<" extrapolationT.isValid "<<extrapolationT.isValid()<<std::endl;
01644
01645 tsosMuonIf = 2;
01646 Rt = GlobalVector((extrapolationT.globalPosition()).x(),
01647 (extrapolationT.globalPosition()).y(),
01648 (extrapolationT.globalPosition()).z());
01649
01650 Pt = extrapolationT.globalMomentum();
01651
01652 GRm = GlobalVector((outerMuTSOS.globalPosition()).x(),
01653 (outerMuTSOS.globalPosition()).y(),
01654 (outerMuTSOS.globalPosition()).z());
01655 GPm = outerMuTSOS.globalMomentum();
01656 Rt0 = GlobalVector((innerTrackTSOS.globalPosition()).x(),
01657 (innerTrackTSOS.globalPosition()).y(),
01658 (innerTrackTSOS.globalPosition()).z());
01659 Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() +
01660 outerMuTSOS.cartesianError().matrix());
01661 C0 = AlgebraicSymMatrix66(innerTrackTSOS.cartesianError().matrix());
01662 Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
01663 C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
01664
01665 if(refitMuon_)
01666 GlobalTrackerMuonAlignment::muonFitter(itMuon->outerTrack(),muTT,alongMomentum,muonFittedTSOS);
01667
01668 if(false & debug_){
01669
01670 PropagationDirectionChooser Chooser;
01671 std::cout<<" propDirCh "
01672 <<Chooser.operator()(*innerTrackTSOS.freeState(), refSurface)
01673 <<" Ch == oppisite "<<(oppositeToMomentum ==
01674 Chooser.operator()(*innerTrackTSOS.freeState(), refSurface))
01675 <<std::endl;
01676 std::cout<<" --- Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
01677 <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
01678 }
01679 }
01680
01681 else if((distanceOutOut <= distanceInOut) & (distanceOutOut <= distanceInIn) &
01682 (distanceOutOut <= distanceOutIn)){
01683 if(debug_) std::cout<<" ----- Out - Out -----"<<std::endl;
01684
01685
01686 continue;
01687
01688 const Surface& refSurface = outerMuTSOS.surface();
01689 ReferenceCountingPointer<TangentPlane>
01690 tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
01691 Nl = tpMuGlobal->normalVector();
01692
01693
01694 extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
01695
01696 if(!extrapolationT.isValid()){
01697 if(alarm)
01698 std::cout<<" !!!!!!!!! Catastrophe Out-Out extrapolationT.isValid "
01699 <<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
01700 <<std::endl;
01701 continue;
01702 }
01703 if(debug_)
01704 std::cout<<" extrapolationT.isValid "<<extrapolationT.isValid()<<std::endl;
01705
01706 tsosMuonIf = 2;
01707 Rt = GlobalVector((extrapolationT.globalPosition()).x(),
01708 (extrapolationT.globalPosition()).y(),
01709 (extrapolationT.globalPosition()).z());
01710
01711 Pt = extrapolationT.globalMomentum();
01712
01713 GRm = GlobalVector((outerMuTSOS.globalPosition()).x(),
01714 (outerMuTSOS.globalPosition()).y(),
01715 (outerMuTSOS.globalPosition()).z());
01716 GPm = outerMuTSOS.globalMomentum();
01717 Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
01718 (outerTrackTSOS.globalPosition()).y(),
01719 (outerTrackTSOS.globalPosition()).z());
01720 Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() +
01721 outerMuTSOS.cartesianError().matrix());
01722 C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
01723 Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
01724 C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
01725
01726 if(debug_){
01727 PropagationDirectionChooser Chooser;
01728 std::cout<<" propDirCh "
01729 <<Chooser.operator()(*outerTrackTSOS.freeState(), refSurface)
01730 <<" Ch == along "<<(alongMomentum ==
01731 Chooser.operator()(*outerTrackTSOS.freeState(), refSurface))
01732 <<std::endl;
01733 std::cout<<" --- Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
01734 <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
01735 std::cout<<" Nornal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
01736 <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
01737 }
01738 }
01739 else{
01740 if(alarm)
01741 std::cout<<" ------- !!!!!!!!!!!!!!! No proper Out - In \a\a\a\a\a\a\a"<<std::endl;
01742 continue;
01743 }
01744
01745 }
01746
01747 if(tsosMuonIf == 0) {if(info) {std::cout<<"No tsosMuon !!!!!!"<<std::endl; continue;}}
01748 TrajectoryStateOnSurface tsosMuon;
01749 if(tsosMuonIf == 1) tsosMuon = muTT.innermostMeasurementState();
01750 else tsosMuon = muTT.outermostMeasurementState();
01751
01752
01753 AlgebraicVector4 LPRm;
01754 GlobalTrackerMuonAlignment::misalignMuonL
01755 (GRm, GPm, Nl, Rt, Rm, Pm, LPRm, extrapolationT, tsosMuon);
01756
01757 if(refitTrack_){
01758 if(!trackFittedTSOS.isValid()){
01759 if(info) std::cout<<" ================= trackFittedTSOS notValid !!!!!!!! "<<std::endl;
01760 continue;}
01761 if(debug_) this->debugTrajectorySOS(" trackFittedTSOS ", trackFittedTSOS);
01762 }
01763
01764 if(refitMuon_){
01765 if(!muonFittedTSOS.isValid()){
01766 if(info) std::cout<<" ================= muonFittedTSOS notValid !!!!!!!! "<<std::endl;
01767 continue;}
01768 if(debug_) this->debugTrajectorySOS(" muonFittedTSOS ", muonFittedTSOS);
01769 Rm = GlobalVector((muonFittedTSOS.globalPosition()).x(),
01770 (muonFittedTSOS.globalPosition()).y(),
01771 (muonFittedTSOS.globalPosition()).z());
01772 Pm = muonFittedTSOS.globalMomentum();
01773 LPRm = AlgebraicVector4(muonFittedTSOS.localParameters().vector()(1),
01774 muonFittedTSOS.localParameters().vector()(2),
01775 muonFittedTSOS.localParameters().vector()(3),
01776 muonFittedTSOS.localParameters().vector()(4));
01777 }
01778 GlobalVector resR = Rm - Rt;
01779 GlobalVector resP0 = Pm - Pt;
01780 GlobalVector resP = Pm / Pm.mag() - Pt / Pt.mag();
01781 float RelMomResidual = (Pm.mag() - Pt.mag()) / (Pt.mag() + 1.e-6);;
01782
01783 AlgebraicVector6 Vm;
01784 Vm(0) = resR.x(); Vm(1) = resR.y(); Vm(2) = resR.z();
01785 Vm(3) = resP0.x(); Vm(4) = resP0.y(); Vm(5) = resP0.z();
01786 float Rmuon = Rm.perp();
01787 float Zmuon = Rm.z();
01788 float alfa_x = atan2(Nl.x(),Nl.y())*57.29578;
01789
01790 if(debug_){
01791 std::cout<<" Nx Ny Nz alfa_x "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "<<alfa_x<<std::endl;
01792
01793
01794 }
01795
01796 double chi_d = 0;
01797 for(int i=0; i<=5; i++) chi_d += Vm(i)*Vm(i)/Cm(i,i);
01798
01799 AlgebraicVector5 Vml(tsosMuon.localParameters().vector()
01800 - extrapolationT.localParameters().vector());
01801 AlgebraicSymMatrix55 m(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
01802 AlgebraicSymMatrix55 Cml(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
01803 bool ierrLoc = !m.Invert();
01804 if (ierrLoc && debug_ && info) {
01805 std::cout<< " ==== Error inverting Local covariance matrix ==== "<<std::cout;
01806 continue;}
01807 double chi_Loc = ROOT::Math::Similarity(Vml,m);
01808 if(debug_)
01809 std::cout<<" chi_Loc px/pz/err "<<chi_Loc<<" "<<Vml(1)/sqrt(Cml(1,1))<<std::endl;
01810
01811 if(Pt.mag() < 15.) continue;
01812 if(Pm.mag() < 5.) continue;
01813
01814
01815
01816
01817
01818
01819
01820
01821
01822
01823
01824
01825
01826
01827 if(fabs(resR.x()) > 20.) continue;
01828 if(fabs(resR.y()) > 20.) continue;
01829 if(fabs(resR.z()) > 20.) continue;
01830 if(fabs(resR.mag()) > 30.) continue;
01831 if(fabs(resP.x()) > 0.06) continue;
01832 if(fabs(resP.y()) > 0.06) continue;
01833 if(fabs(resP.z()) > 0.06) continue;
01834 if(chi_d > 40.) continue;
01835
01836
01837
01838
01839
01840
01841
01842
01843
01844
01845
01846
01847
01848
01849
01850
01851
01852
01853
01854 if(Rmuon < 120. || Rmuon > 450.) continue;
01855 if(Zmuon < -720. || Zmuon > 720.) continue;
01856 MuSelect = " Barrel+EndCaps";
01857
01858 if(debug_)
01859 std::cout<<" .............. passed all cuts"<<std::endl;
01860
01861 N_track++;
01862
01863
01864 GlobalTrackerMuonAlignment::gradientGlobalAlg(Rt, Pt, Rm, Nl, Cm);
01865 GlobalTrackerMuonAlignment::gradientGlobal(Rt, Pt, Rm, Pm, Nl, Cm);
01866
01867 CLHEP::HepSymMatrix covLoc(4,0);
01868 for(int i=1; i<=4; i++)
01869 for(int j=1; j<=i; j++){
01870 covLoc(i,j) = (tsosMuon.localError().matrix()
01871 + extrapolationT.localError().matrix())(i,j);
01872
01873 }
01874
01875 const Surface& refSurface = tsosMuon.surface();
01876 CLHEP::HepMatrix rotLoc (3,3,0);
01877 rotLoc(1,1) = refSurface.rotation().xx();
01878 rotLoc(1,2) = refSurface.rotation().xy();
01879 rotLoc(1,3) = refSurface.rotation().xz();
01880
01881 rotLoc(2,1) = refSurface.rotation().yx();
01882 rotLoc(2,2) = refSurface.rotation().yy();
01883 rotLoc(2,3) = refSurface.rotation().yz();
01884
01885 rotLoc(3,1) = refSurface.rotation().zx();
01886 rotLoc(3,2) = refSurface.rotation().zy();
01887 rotLoc(3,3) = refSurface.rotation().zz();
01888
01889 CLHEP::HepVector posLoc(3);
01890 posLoc(1) = refSurface.position().x();
01891 posLoc(2) = refSurface.position().y();
01892 posLoc(3) = refSurface.position().z();
01893
01894 GlobalTrackerMuonAlignment::gradientLocal(Rt, Pt, Rm, Pm, Nl, covLoc, rotLoc, posLoc, LPRm);
01895
01896 if(debug_){
01897 std::cout<<" Norm "<<Nl<<std::endl;
01898 std::cout<<" posLoc "<<posLoc.T()<<std::endl;
01899 std::cout<<" rotLoc "<<rotLoc<<std::endl;
01900 }
01901
01902
01903 histo->Fill(itMuon->track()->pt());
01904
01905
01906 histo2->Fill(Pt.mag());
01907 histo3->Fill((PI/2.-itMuon->track()->outerTheta()));
01908 histo4->Fill(itMuon->track()->phi());
01909 histo5->Fill(Rmuon);
01910 histo6->Fill(Zmuon);
01911 histo7->Fill(RelMomResidual);
01912
01913 histo8->Fill(chi_d);
01914
01915 histo101->Fill(Zmuon, Rmuon);
01916 histo101->Fill(Rt0.z(), Rt0.perp());
01917 histo102->Fill(Rt0.x(), Rt0.y());
01918 histo102->Fill(Rm.x(), Rm.y());
01919
01920 histo11->Fill(resR.mag());
01921 if(fabs(Nl.x()) < 0.98) histo12->Fill(resR.x());
01922 if(fabs(Nl.y()) < 0.98) histo13->Fill(resR.y());
01923 if(fabs(Nl.z()) < 0.98) histo14->Fill(resR.z());
01924 histo15->Fill(resP.x());
01925 histo16->Fill(resP.y());
01926 histo17->Fill(resP.z());
01927
01928 if((fabs(Nl.x()) < 0.98) && (fabs(Nl.y()) < 0.98) &&(fabs(Nl.z()) < 0.98))
01929 {
01930 histo18->Fill(sqrt(C0(0,0)));
01931 histo19->Fill(sqrt(C1(0,0)));
01932 histo20->Fill(sqrt(C1(0,0)+Ce(0,0)));
01933 }
01934 if(fabs(Nl.x()) < 0.98) histo21->Fill(Vm(0)/sqrt(Cm(0,0)));
01935 if(fabs(Nl.y()) < 0.98) histo22->Fill(Vm(1)/sqrt(Cm(1,1)));
01936 if(fabs(Nl.z()) < 0.98) histo23->Fill(Vm(2)/sqrt(Cm(2,2)));
01937 histo24->Fill(Vm(3)/sqrt(C1(3,3)+Ce(3,3)));
01938 histo25->Fill(Vm(4)/sqrt(C1(4,4)+Ce(4,4)));
01939 histo26->Fill(Vm(5)/sqrt(C1(5,5)+Ce(5,5)));
01940 histo27->Fill(Nl.x());
01941 histo28->Fill(Nl.y());
01942 histo29->Fill(lenghtTrack);
01943 histo30->Fill(lenghtMuon);
01944 histo31->Fill(chi_Loc);
01945 histo32->Fill(Vml(1)/sqrt(Cml(1,1)));
01946 histo33->Fill(Vml(2)/sqrt(Cml(2,2)));
01947 histo34->Fill(Vml(3)/sqrt(Cml(3,3)));
01948 histo35->Fill(Vml(4)/sqrt(Cml(4,4)));
01949
01950 if (debug_) {
01951
01952 std::cout<<" diag 0[ "<<C0(0,0)<<" "<<C0(1,1)<<" "<<C0(2,2)<<" "<<C0(3,3)<<" "
01953 <<C0(4,4)<<" "<<C0(5,5)<<" ]"<<std::endl;
01954 std::cout<<" diag e[ "<<Ce(0,0)<<" "<<Ce(1,1)<<" "<<Ce(2,2)<<" "<<Ce(3,3)<<" "
01955 <<Ce(4,4)<<" "<<Ce(5,5)<<" ]"<<std::endl;
01956 std::cout<<" diag 1[ "<<C1(0,0)<<" "<<C1(1,1)<<" "<<C1(2,2)<<" "<<C1(3,3)<<" "
01957 <<C1(4,4)<<" "<<C1(5,5)<<" ]"<<std::endl;
01958 std::cout<<" Rm "<<Rm.x()<<" "<<Rm.y()<<" "<<Rm.z()
01959 <<" Pm "<<Pm.x()<<" "<<Pm.y()<<" "<<Pm.z()<<std::endl;
01960 std::cout<<" Rt "<<Rt.x()<<" "<<Rt.y()<<" "<<Rt.z()
01961 <<" Pt "<<Pt.x()<<" "<<Pt.y()<<" "<<Pt.z()<<std::endl;
01962 std::cout<<" Nl*(Rm-Rt) "<<Nl.dot(Rm-Rt)<<std::endl;
01963 std::cout<<" resR "<<resR.x()<<" "<<resR.y()<<" "<<resR.z()
01964 <<" resP "<<resP.x()<<" "<<resP.y()<<" "<<resP.z()<<std::endl;
01965 std::cout<<" Rm-t "<<(Rm-Rt).x()<<" "<<(Rm-Rt).y()<<" "<<(Rm-Rt).z()
01966 <<" Pm-t "<<(Pm-Pt).x()<<" "<<(Pm-Pt).y()<<" "<<(Pm-Pt).z()<<std::endl;
01967 std::cout<<" Vm "<<Vm<<std::endl;
01968 std::cout<<" +- "<<sqrt(Cm(0,0))<<" "<<sqrt(Cm(1,1))<<" "<<sqrt(Cm(2,2))<<" "
01969 <<sqrt(Cm(3,3))<<" "<<sqrt(Cm(4,4))<<" "<<sqrt(Cm(5,5))<<std::endl;
01970 std::cout<<" Pmuon Ptrack dP/Ptrack "<<itMuon->outerTrack()->p()<<" "
01971 <<itMuon->track()->outerP()<<" "
01972 <<(itMuon->outerTrack()->p() -
01973 itMuon->track()->outerP())/itMuon->track()->outerP()<<std::endl;
01974 std::cout<<" cov matrix "<<std::endl;
01975 std::cout<<Cm<<std::endl;
01976 std::cout<<" diag [ "<<Cm(0,0)<<" "<<Cm(1,1)<<" "<<Cm(2,2)<<" "<<Cm(3,3)<<" "
01977 <<Cm(4,4)<<" "<<Cm(5,5)<<" ]"<<std::endl;
01978
01979 AlgebraicSymMatrix66 Ro;
01980 double Diag[6];
01981 for(int i=0; i<=5; i++) Diag[i] = sqrt(Cm(i,i));
01982 for(int i=0; i<=5; i++)
01983 for(int j=0; j<=5; j++)
01984 Ro(i,j) = Cm(i,j)/Diag[i]/Diag[j];
01985 std::cout<<" correlation matrix "<<std::endl;
01986 std::cout<<Ro<<std::endl;
01987
01988 AlgebraicSymMatrix66 CmI;
01989 for(int i=0; i<=5; i++)
01990 for(int j=0; j<=5; j++)
01991 CmI(i,j) = Cm(i,j);
01992
01993 bool ierr = !CmI.Invert();
01994 if( ierr ) { if(alarm || debug_)
01995 std::cout<<" Error inverse covariance matrix !!!!!!!!!!!"<<std::endl;
01996 continue;
01997 }
01998 std::cout<<" inverse cov matrix "<<std::endl;
01999 std::cout<<Cm<<std::endl;
02000
02001 double chi = ROOT::Math::Similarity(Vm, CmI);
02002 std::cout<<" chi chi_d "<<chi<<" "<<chi_d<<std::endl;
02003 }
02004
02005 }
02006
02007 }
02008
02009
02010
02011 void
02012 GlobalTrackerMuonAlignment::gradientGlobalAlg(GlobalVector& Rt, GlobalVector& Pt,
02013 GlobalVector& Rm, GlobalVector& Nl,
02014 AlgebraicSymMatrix66& Cm)
02015 {
02016
02017
02018
02019 AlgebraicMatrix33 Jac;
02020 AlgebraicVector3 Wi, R_m, R_t, P_t, Norm, dR;
02021
02022 R_m(0) = Rm.x(); R_m(1) = Rm.y(); R_m(2) = Rm.z();
02023 R_t(0) = Rt.x(); R_t(1) = Rt.y(); R_t(2) = Rt.z();
02024 P_t(0) = Pt.x(); P_t(1) = Pt.y(); P_t(2) = Pt.z();
02025 Norm(0) = Nl.x(); Norm(1) = Nl.y(); Norm(2) = Nl.z();
02026
02027 for(int i=0; i<=2; i++){
02028 if(Cm(i,i) > 1.e-20)
02029 Wi(i) = 1./Cm(i,i);
02030 else Wi(i) = 1.e-10;
02031 dR(i) = R_m(i) - R_t(i);
02032 }
02033
02034 float PtN = P_t(0)*Norm(0) + P_t(1)*Norm(1) + P_t(2)*Norm(2);
02035
02036 Jac(0,0) = 1. - P_t(0)*Norm(0)/PtN;
02037 Jac(0,1) = - P_t(0)*Norm(1)/PtN;
02038 Jac(0,2) = - P_t(0)*Norm(2)/PtN;
02039
02040 Jac(1,0) = - P_t(1)*Norm(0)/PtN;
02041 Jac(1,1) = 1. - P_t(1)*Norm(1)/PtN;
02042 Jac(1,2) = - P_t(1)*Norm(2)/PtN;
02043
02044 Jac(2,0) = - P_t(2)*Norm(0)/PtN;
02045 Jac(2,1) = - P_t(2)*Norm(1)/PtN;
02046 Jac(2,2) = 1. - P_t(2)*Norm(2)/PtN;
02047
02048 AlgebraicSymMatrix33 Itr;
02049
02050 for(int i=0; i<=2; i++)
02051 for(int j=0; j<=2; j++){
02052 if(j < i) continue;
02053 Itr(i,j) = 0.;
02054
02055 for(int k=0; k<=2; k++){
02056 Itr(i,j) += Jac(k,i)*Wi(k)*Jac(k,j);}}
02057
02058 for(int i=0; i<=2; i++)
02059 for(int j=0; j<=2; j++){
02060 if(j < i) continue;
02061 Inf(i,j) += Itr(i,j);}
02062
02063 AlgebraicVector3 Gtr(0., 0., 0.);
02064 for(int i=0; i<=2; i++)
02065 for(int k=0; k<=2; k++) Gtr(i) += dR(k)*Wi(k)*Jac(k,i);
02066 for(int i=0; i<=2; i++) Gfr(i) += Gtr(i);
02067
02068 if(debug_){
02069 std::cout<<" Wi "<<Wi<<std::endl;
02070 std::cout<<" N "<<Norm<<std::endl;
02071 std::cout<<" P_t "<<P_t<<std::endl;
02072 std::cout<<" (Pt*N) "<<PtN<<std::endl;
02073 std::cout<<" dR "<<dR<<std::endl;
02074 std::cout<<" +/- "<<1./sqrt(Wi(0))<<" "<<1./sqrt(Wi(1))<<" "<<1./sqrt(Wi(2))
02075 <<" "<<std::endl;
02076 std::cout<<" Jacobian dr/ddx "<<std::endl;
02077 std::cout<<Jac<<std::endl;
02078 std::cout<<" G-- "<<Gtr<<std::endl;
02079 std::cout<<" Itrack "<<std::endl;
02080 std::cout<<Itr<<std::endl;
02081 std::cout<<" Gfr "<<Gfr<<std::endl;
02082 std::cout<<" -- Inf --"<<std::endl;
02083 std::cout<<Inf<<std::endl;
02084 }
02085
02086 return;
02087 }
02088
02089
02090 void
02091 GlobalTrackerMuonAlignment::gradientGlobal(GlobalVector& GRt, GlobalVector& GPt,
02092 GlobalVector& GRm, GlobalVector& GPm,
02093 GlobalVector& GNorm, AlgebraicSymMatrix66 & GCov)
02094 {
02095
02096
02097
02098
02099
02100 bool alarm = false;
02101 bool info = false;
02102
02103 int Nd = 6;
02104
02105
02106 CLHEP::HepSymMatrix w(Nd,0);
02107 for (int i=1; i <= Nd; i++)
02108 for (int j=1; j <= Nd; j++){
02109 if(j <= i ) w(i,j) = GCov(i-1, j-1);
02110
02111
02112 if((i == j) && (i<=3) && (GCov(i-1, j-1) < 1.e-20)) w(i,j) = 1.e20;
02113 if(i != j) w(i,j) = 0.;
02114 }
02115
02116
02117
02118
02119 CLHEP::HepVector V(Nd), Rt(3), Pt(3), Rm(3), Pm(3), Norm(3);
02120 Rt(1) = GRt.x(); Rt(2) = GRt.y(); Rt(3) = GRt.z();
02121 Pt(1) = GPt.x(); Pt(2) = GPt.y(); Pt(3) = GPt.z();
02122 Rm(1) = GRm.x(); Rm(2) = GRm.y(); Rm(3) = GRm.z();
02123 Pm(1) = GPm.x(); Pm(2) = GPm.y(); Pm(3) = GPm.z();
02124 Norm(1) = GNorm.x(); Norm(2) = GNorm.y(); Norm(3) = GNorm.z();
02125
02126 V = dsum(Rm - Rt, Pm - Pt) ;
02127
02128
02129 double PmN = CLHEP_dot(Pm, Norm);
02130
02131 CLHEP::HepMatrix Jac(Nd,Nd,0);
02132 for (int i=1; i <= 3; i++)
02133 for (int j=1; j <= 3; j++){
02134 Jac(i,j) = Pm(i)*Norm(j) / PmN;
02135 if(i == j ) Jac(i,j) -= 1.;
02136 }
02137
02138
02139 Jac(4,4) = 0.;
02140 Jac(5,4) = -Pm(3);
02141 Jac(6,4) = Pm(2);
02142 Jac(4,5) = Pm(3);
02143 Jac(5,5) = 0.;
02144 Jac(6,5) = -Pm(1);
02145 Jac(4,6) = -Pm(2);
02146 Jac(5,6) = Pm(1);
02147 Jac(6,6) = 0.;
02148
02149 CLHEP::HepVector dsda(3);
02150 dsda(1) = (Norm(2)*Rm(3) - Norm(3)*Rm(2)) / PmN;
02151 dsda(2) = (Norm(3)*Rm(1) - Norm(1)*Rm(3)) / PmN;
02152 dsda(3) = (Norm(1)*Rm(2) - Norm(2)*Rm(1)) / PmN;
02153
02154
02155 Jac(1,4) = Pm(1)*dsda(1);
02156 Jac(2,4) = -Rm(3) + Pm(2)*dsda(1);
02157 Jac(3,4) = Rm(2) + Pm(3)*dsda(1);
02158
02159 Jac(1,5) = Rm(3) + Pm(1)*dsda(2);
02160 Jac(2,5) = Pm(2)*dsda(2);
02161 Jac(3,5) = -Rm(1) + Pm(3)*dsda(2);
02162
02163 Jac(1,6) = -Rm(2) + Pm(1)*dsda(3);
02164 Jac(2,6) = Rm(1) + Pm(2)*dsda(3);
02165 Jac(3,6) = Pm(3)*dsda(3);
02166
02167 CLHEP::HepSymMatrix W(Nd,0);
02168 int ierr;
02169 W = w.inverse(ierr);
02170 if(ierr != 0) { if(alarm)
02171 std::cout<<" gradientGlobal: inversion of matrix w fail "<<std::endl;
02172 return;
02173 }
02174
02175 CLHEP::HepMatrix W_Jac(Nd,Nd,0);
02176 W_Jac = Jac.T() * W;
02177
02178 CLHEP::HepVector grad3(Nd);
02179 grad3 = W_Jac * V;
02180
02181 CLHEP::HepMatrix hess3(Nd,Nd);
02182 hess3 = Jac.T() * W * Jac;
02183
02184
02185 Grad += grad3;
02186 Hess += hess3;
02187
02188 CLHEP::HepVector d3I = CLHEP::solve(Hess, -Grad);
02189 int iEr3I;
02190 CLHEP::HepMatrix Errd3I = Hess.inverse(iEr3I);
02191 if( iEr3I != 0) { if(alarm)
02192 std::cout<<" gradientGlobal error inverse Hess matrix !!!!!!!!!!!"<<std::endl;
02193 }
02194
02195 if(info || debug_){
02196 std::cout<<" dG "<<d3I(1)<<" "<<d3I(2)<<" "<<d3I(3)<<" "<<d3I(4)<<" "
02197 <<d3I(5)<<" "<<d3I(6)<<std::endl;;
02198 std::cout<<" +- "<<sqrt(Errd3I(1,1))<<" "<<sqrt(Errd3I(2,2))<<" "<<sqrt(Errd3I(3,3))
02199 <<" "<<sqrt(Errd3I(4,4))<<" "<<sqrt(Errd3I(5,5))<<" "<<sqrt(Errd3I(6,6))
02200 <<std::endl;
02201 }
02202
02203 #ifdef CHECK_OF_DERIVATIVES
02204
02205
02206 CLHEP::HepVector d(3,0), a(3,0);
02207 CLHEP::HepMatrix T(3,3,1);
02208
02209 CLHEP::HepMatrix Ti = T.T();
02210
02211
02212 double A = CLHEP_dot(Ti*Pm, Norm);
02213 double B = CLHEP_dot((Rt -Ti*Rm + Ti*d), Norm);
02214 double s0 = B / A;
02215
02216 CLHEP::HepVector r0(3,0), p0(3,0);
02217 r0 = Ti*Rm - Ti*d + s0*(Ti*Pm) - Rt;
02218 p0 = Ti*Pm - Pt;
02219
02220 double delta = 0.0001;
02221
02222 int ii = 3; d(ii) += delta;
02223
02224
02225
02226 Ti = T.T();
02227
02228
02229 A = CLHEP_dot(Ti*Pm, Norm);
02230 B = CLHEP_dot((Rt -Ti*Rm + Ti*d), Norm);
02231 double s = B / A;
02232
02233 CLHEP::HepVector r(3,0), p(3,0);
02234 r = Ti*Rm - Ti*d + s*(Ti*Pm) - Rt;
02235 p = Ti*Pm - Pt;
02236
02237 std::cout<<" s0 s num dsda("<<ii<<") "<<s0<<" "<<s<<" "
02238 <<(s-s0)/delta<<" "<<dsda(ii)<<endl;
02239
02240 std::cout<<" -- An d(r,p)/d("<<ii<<") "<<Jac(1,ii)<<" "<<Jac(2,ii)<<" "
02241 <<Jac(3,ii)<<" "<<Jac(4,ii)<<" "<<Jac(5,ii)<<" "
02242 <<Jac(6,ii)<<std::endl;
02243 std::cout<<" Nu d(r,p)/d("<<ii<<") "<<(r(1)-r0(1))/delta<<" "
02244 <<(r(2)-r0(2))/delta<<" "<<(r(3)-r0(3))/delta<<" "
02245 <<(p(1)-p0(1))/delta<<" "<<(p(2)-p0(2))/delta<<" "
02246 <<(p(3)-p0(3))/delta<<std::endl;
02247
02248 std::cout<<" -- An d(r,p)/a("<<ii<<") "<<Jac(1,ii+3)<<" "<<Jac(2,ii+3)<<" "
02249 <<Jac(3,ii+3)<<" "<<Jac(4,ii+3)<<" "<<Jac(5,ii+3)<<" "
02250 <<Jac(6,ii+3)<<std::endl;
02251 std::cout<<" Nu d(r,p)/a("<<ii<<") "<<(r(1)-r0(1))/delta<<" "
02252 <<(r(2)-r0(2))/delta<<" "<<(r(3)-r0(3))/delta<<" "
02253 <<(p(1)-p0(1))/delta<<" "<<(p(2)-p0(2))/delta<<" "
02254 <<(p(3)-p0(3))/delta<<std::endl;
02255
02256 #endif
02257
02258 return;
02259 }
02260
02261
02262
02263 void
02264 GlobalTrackerMuonAlignment::gradientLocal(GlobalVector& GRt, GlobalVector& GPt,
02265 GlobalVector& GRm, GlobalVector& GPm,
02266 GlobalVector& GNorm,
02267 CLHEP::HepSymMatrix& covLoc,
02268 CLHEP::HepMatrix& rotLoc,
02269 CLHEP::HepVector& R0,
02270 AlgebraicVector4& LPRm)
02271 {
02272
02273
02274
02275
02276 bool alarm = true;
02277
02278 bool info = false;
02279
02280 if(debug_)
02281 std::cout<<" gradientLocal "<<std::endl;
02282
02283
02284
02285
02286
02287
02288
02289
02290
02291
02292
02293
02294
02295
02296
02297
02298
02299
02300 CLHEP::HepVector Rt(3), Pt(3), Rm(3), Pm(3), Norm(3);
02301 Rt(1) = GRt.x(); Rt(2) = GRt.y(); Rt(3) = GRt.z();
02302 Pt(1) = GPt.x(); Pt(2) = GPt.y(); Pt(3) = GPt.z();
02303 Rm(1) = GRm.x(); Rm(2) = GRm.y(); Rm(3) = GRm.z();
02304 Pm(1) = GPm.x(); Pm(2) = GPm.y(); Pm(3) = GPm.z();
02305 Norm(1) = GNorm.x(); Norm(2) = GNorm.y(); Norm(3) = GNorm.z();
02306
02307 CLHEP::HepVector V(4), Rml(3), Pml(3), Rtl(3), Ptl(3);
02308
02309
02310
02311
02312
02313
02314
02315 Rml = rotLoc * (Rm - R0);
02316 Rtl = rotLoc * (Rt - R0);
02317 Pml = rotLoc * Pm;
02318 Ptl = rotLoc * Pt;
02319
02320 V(1) = LPRm(0) - Ptl(1)/Ptl(3);
02321 V(2) = LPRm(1) - Ptl(2)/Ptl(3);
02322 V(3) = LPRm(2) - Rtl(1);
02323 V(4) = LPRm(3) - Rtl(2);
02324
02325
02326
02327
02328
02329
02330
02331
02332
02333
02334
02335
02336
02337
02338 CLHEP::HepSymMatrix W = covLoc;
02339
02340 int ierr;
02341 W.invert(ierr);
02342 if(ierr != 0) { if(alarm)
02343 std::cout<<" gradientLocal: inversion of matrix W fail "<<std::endl;
02344 return;
02345 }
02346
02347
02348
02349
02350
02351
02352
02353
02354
02355 CLHEP::HepMatrix JacToLoc(4,6,0);
02356 for(int i=1; i<=2; i++)
02357 for(int j=1; j<=3; j++){
02358 JacToLoc(i,j+3) = (rotLoc(i,j) - rotLoc(3,j)*Pml(i)/Pml(3))/Pml(3);
02359 JacToLoc(i+2,j) = rotLoc(i,j);
02360 }
02361
02362
02363
02364 double PmN = CLHEP_dot(Pm, Norm);
02365
02366 CLHEP::HepMatrix Jac(6,6,0);
02367 for (int i=1; i <= 3; i++)
02368 for (int j=1; j <= 3; j++){
02369 Jac(i,j) = Pm(i)*Norm(j) / PmN;
02370 if(i == j ) Jac(i,j) -= 1.;
02371 }
02372
02373
02374 Jac(4,4) = 0.;
02375 Jac(5,4) = -Pm(3);
02376 Jac(6,4) = Pm(2);
02377 Jac(4,5) = Pm(3);
02378 Jac(5,5) = 0.;
02379 Jac(6,5) = -Pm(1);
02380 Jac(4,6) = -Pm(2);
02381 Jac(5,6) = Pm(1);
02382 Jac(6,6) = 0.;
02383
02384 CLHEP::HepVector dsda(3);
02385 dsda(1) = (Norm(2)*Rm(3) - Norm(3)*Rm(2)) / PmN;
02386 dsda(2) = (Norm(3)*Rm(1) - Norm(1)*Rm(3)) / PmN;
02387 dsda(3) = (Norm(1)*Rm(2) - Norm(2)*Rm(1)) / PmN;
02388
02389
02390 Jac(1,4) = Pm(1)*dsda(1);
02391 Jac(2,4) = -Rm(3) + Pm(2)*dsda(1);
02392 Jac(3,4) = Rm(2) + Pm(3)*dsda(1);
02393
02394 Jac(1,5) = Rm(3) + Pm(1)*dsda(2);
02395 Jac(2,5) = Pm(2)*dsda(2);
02396 Jac(3,5) = -Rm(1) + Pm(3)*dsda(2);
02397
02398 Jac(1,6) = -Rm(2) + Pm(1)*dsda(3);
02399 Jac(2,6) = Rm(1) + Pm(2)*dsda(3);
02400 Jac(3,6) = Pm(3)*dsda(3);
02401
02402
02403 CLHEP::HepMatrix JacCorLoc(4,6,0);
02404 JacCorLoc = JacToLoc * Jac;
02405
02406
02407 CLHEP::HepMatrix W_Jac(6,4,0);
02408 W_Jac = JacCorLoc.T() * W;
02409
02410 CLHEP::HepVector gradL(6);
02411 gradL = W_Jac * V;
02412
02413 CLHEP::HepMatrix hessL(6,6);
02414 hessL = JacCorLoc.T() * W * JacCorLoc;
02415
02416 GradL += gradL;
02417 HessL += hessL;
02418
02419 CLHEP::HepVector dLI = CLHEP::solve(HessL, -GradL);
02420 int iErI;
02421 CLHEP::HepMatrix ErrdLI = HessL.inverse(iErI);
02422 if( iErI != 0) { if(alarm)
02423 std::cout<<" gradLocal Error inverse Hess matrix !!!!!!!!!!!"<<std::endl;
02424 }
02425
02426 if(info || debug_){
02427 std::cout<<" dL "<<dLI(1)<<" "<<dLI(2)<<" "<<dLI(3)<<" "<<dLI(4)<<" "
02428 <<dLI(5)<<" "<<dLI(6)<<std::endl;;
02429 std::cout<<" +- "<<sqrt(ErrdLI(1,1))<<" "<<sqrt(ErrdLI(2,2))<<" "<<sqrt(ErrdLI(3,3))
02430 <<" "<<sqrt(ErrdLI(4,4))<<" "<<sqrt(ErrdLI(5,5))<<" "<<sqrt(ErrdLI(6,6))
02431 <<std::endl;
02432 }
02433
02434 if(debug_){
02435
02436
02437
02438
02439
02440
02441
02442
02443 int i=5;
02444 if(GNorm.z() > 0.95)
02445 std::cout<<" Ecap1 N "<<GNorm<<std::endl;
02446 else if(GNorm.z() < -0.95)
02447 std::cout<<" Ecap2 N "<<GNorm<<std::endl;
02448 else
02449 std::cout<<" Barrel N "<<GNorm<<std::endl;
02450 std::cout<<" JacCLo(i,"<<i<<") = {"<<JacCorLoc(1,i)<<" "<<JacCorLoc(2,i)<<" "
02451 <<JacCorLoc(3,i)<<" "<<JacCorLoc(4,i)<<"}"<<std::endl;
02452 std::cout<<" rotLoc "<<rotLoc<<std::endl;
02453 std::cout<<" position "<<R0<<std::endl;
02454 std::cout<<" Pm,l "<<Pm.T()<<" "<<Pml(1)/Pml(3)<<" "<<Pml(2)/Pml(3)<<std::endl;
02455 std::cout<<" Pt,l "<<Pt.T()<<" "<<Ptl(1)/Ptl(3)<<" "<<Ptl(2)/Ptl(3)<<std::endl;
02456 std::cout<<" V "<<V.T()<<std::endl;
02457 std::cout<<" Cov \n"<<covLoc<<std::endl;
02458 std::cout<<" W*Cov "<<W*covLoc<<std::endl;
02459
02460
02461
02462 std::cout<<" JacToLoc "<<JacToLoc<<std::endl;
02463
02464 }
02465
02466 #ifdef CHECK_OF_JACOBIAN_CARTESIAN_TO_LOCAL
02467
02468 CLHEP::HepVector V0(4,0);
02469 V0(1) = Pml(1)/Pml(3) - Ptl(1)/Ptl(3);
02470 V0(2) = Pml(2)/Pml(3) - Ptl(2)/Ptl(3);
02471 V0(3) = Rml(1) - Rtl(1);
02472 V0(4) = Rml(2) - Rtl(2);
02473 int ii = 3; float delta = 0.01;
02474 CLHEP::HepVector V1(4,0);
02475 if(ii <= 3) {Rm(ii) += delta; Rml = rotLoc * (Rm - R0);}
02476 else {Pm(ii-3) += delta; Pml = rotLoc * Pm;}
02477
02478
02479 V1(1) = Pml(1)/Pml(3) - Ptl(1)/Ptl(3);
02480 V1(2) = Pml(2)/Pml(3) - Ptl(2)/Ptl(3);
02481 V1(3) = Rml(1) - Rtl(1);
02482 V1(4) = Rml(2) - Rtl(2);
02483
02484 if(GNorm.z() > 0.95)
02485 std::cout<<" Ecap1 N "<<GNorm<<std::endl;
02486 else if(GNorm.z() < -0.95)
02487 std::cout<<" Ecap2 N "<<GNorm<<std::endl;
02488 else
02489 std::cout<<" Barrel N "<<GNorm<<std::endl;
02490 std::cout<<" dldc Num (i,"<<ii<<") "<<(V1(1)-V0(1))/delta<<" "<<(V1(2)-V0(2))/delta<<" "
02491 <<(V1(3)-V0(3))/delta<<" "<<(V1(4)-V0(4))/delta<<std::endl;
02492 std::cout<<" dldc Ana (i,"<<ii<<") "<<JacToLoc(1,ii)<<" "<<JacToLoc(2,ii)<<" "
02493 <<JacToLoc(3,ii)<<" "<<JacToLoc(4,ii)<<std::endl;
02494
02495
02496
02497 #endif
02498
02499 return;
02500 }
02501
02502
02503
02504 void
02505 GlobalTrackerMuonAlignment::misalignMuon(GlobalVector& GRm, GlobalVector& GPm, GlobalVector& Nl,
02506 GlobalVector& Rt, GlobalVector& Rm, GlobalVector& Pm)
02507 {
02508 CLHEP::HepVector d(3,0), a(3,0), Norm(3), Rm0(3), Pm0(3), Rm1(3), Pm1(3),
02509 Rmi(3), Pmi(3) ;
02510
02511 d(1)=0.0; d(2)=0.0; d(3)=0.0; a(1)=0.0000; a(2)=0.0000; a(3)=0.0000;
02512
02513
02514
02515
02516
02517
02518
02519
02520
02521
02522
02523 MuGlShift = d; MuGlAngle = a;
02524
02525 if((d(1) == 0.) & (d(2) == 0.) & (d(3) == 0.) &
02526 (a(1) == 0.) & (a(2) == 0.) & (a(3) == 0.)){
02527 Rm = GRm;
02528 Pm = GPm;
02529 if(debug_)
02530 std::cout<<" ...... without misalignment "<<std::endl;
02531 return;
02532 }
02533
02534 Rm0(1) = GRm.x(); Rm0(2) = GRm.y(); Rm0(3) = GRm.z();
02535 Pm0(1) = GPm.x(); Pm0(2) = GPm.y(); Pm0(3) = GPm.z();
02536 Norm(1) = Nl.x(); Norm(2) = Nl.y(); Norm(3) = Nl.z();
02537
02538 CLHEP::HepMatrix T(3,3,1);
02539
02540
02541
02542
02543
02544 double s1 = std::sin(a(1)), c1 = std::cos(a(1));
02545 double s2 = std::sin(a(2)), c2 = std::cos(a(2));
02546 double s3 = std::sin(a(3)), c3 = std::cos(a(3));
02547 T(1,1) = c2 * c3;
02548 T(1,2) = c1 * s3 + s1 * s2 * c3;
02549 T(1,3) = s1 * s3 - c1 * s2 * c3;
02550 T(2,1) = -c2 * s3;
02551 T(2,2) = c1 * c3 - s1 * s2 * s3;
02552 T(2,3) = s1 * c3 + c1 * s2 * s3;
02553 T(3,1) = s2;
02554 T(3,2) = -s1 * c2;
02555 T(3,3) = c1 * c2;
02556
02557
02558
02559 CLHEP::HepMatrix Ti = T.T();
02560 CLHEP::HepVector di = -Ti * d;
02561
02562 CLHEP::HepVector ex0(3,0), ey0(3,0), ez0(3,0), ex(3,0), ey(3,0), ez(3,0);
02563 ex0(1) = 1.; ey0(2) = 1.; ez0(3) = 1.;
02564 ex = Ti*ex0; ey = Ti*ey0; ez = Ti*ez0;
02565
02566 CLHEP::HepVector TiN = Ti * Norm;
02567
02568
02569
02570
02571 double si = CLHEP_dot((Ti*Rm0 - Rm0 + di), TiN) / CLHEP_dot(Pm0, TiN);
02572 Rm1(1) = CLHEP_dot(ex, Rm0 + si*Pm0 - di);
02573 Rm1(2) = CLHEP_dot(ey, Rm0 + si*Pm0 - di);
02574 Rm1(3) = CLHEP_dot(ez, Rm0 + si*Pm0 - di);
02575 Pm1 = T * Pm0;
02576
02577 Rm = GlobalVector(Rm1(1), Rm1(2), Rm1(3));
02578
02579 Pm = GlobalVector(CLHEP_dot(Pm0,ex), CLHEP_dot(Pm0, ey), CLHEP_dot(Pm0,ez));
02580
02581 if(debug_){
02582
02583 std::cout<<" ----- Pm "<<Pm<<std::endl;
02584 std::cout<<" T*Pm0 "<<Pm1.T()<<std::endl;
02585
02586
02587
02588 CLHEP::HepVector Rml = Rm1(1)*ex + Rm1(2)*ey + Rm1(3)*ez + di;
02589
02590 CLHEP::HepVector Rt0(3);
02591 Rt0(1)=Rt.x(); Rt0(2)=Rt.y(); Rt0(3)=Rt.z();
02592
02593
02594 double sl = CLHEP_dot(T*(Rt0 - Rml), T*Norm) / CLHEP_dot(Ti * Pm1, Norm);
02595 CLHEP::HepVector Rl = Rml + sl*(Ti*Pm1);
02596
02597
02598
02599
02600 double A = CLHEP_dot(Ti*Pm1, Norm);
02601 double B = CLHEP_dot(Rt0, Norm) - CLHEP_dot(Ti*Rm1, Norm)
02602 + CLHEP_dot(Ti*d, Norm);
02603 double s = B/A;
02604 CLHEP::HepVector Dr = Ti*Rm1 - Ti*d + s*(Ti*Pm1) - Rm0;
02605 CLHEP::HepVector Dp = Ti*Pm1 - Pm0;
02606
02607 CLHEP::HepVector TiR = Ti * Rm0 + di;
02608 CLHEP::HepVector Ri = T * TiR + d;
02609
02610 std::cout<<" --TTi-0 "<<(Ri-Rm0).T()<<std::endl;
02611 std::cout<<" -- Dr "<<Dr.T()<<endl;
02612 std::cout<<" -- Dp "<<Dp.T()<<endl;
02613
02614
02615
02616
02617
02618
02619
02620 std::cout<<" -- si sl s "<<si<<" "<<sl<<" "<<s<<std::endl;
02621
02622
02623 std::cout<<" -- Rml-(Rm0+si*Pm0) "<<(Rml - Rm0 - si*Pm0).T()<<std::endl;
02624
02625
02626
02627
02628
02629
02630
02631 std::cout<<" --Rl-0 "<<(Rl-Rm0).T()<<std::endl;
02632
02633 }
02634
02635 return;
02636
02637 }
02638
02639
02640 void
02641 GlobalTrackerMuonAlignment::misalignMuonL(GlobalVector& GRm, GlobalVector& GPm, GlobalVector& Nl,
02642 GlobalVector& Rt, GlobalVector& Rm, GlobalVector& Pm,
02643 AlgebraicVector4& Vm,
02644 TrajectoryStateOnSurface& tsosTrack,
02645 TrajectoryStateOnSurface& tsosMuon)
02646 {
02647 CLHEP::HepVector d(3,0), a(3,0), Norm(3), Rm0(3), Pm0(3), Rm1(3), Pm1(3),
02648 Rmi(3), Pmi(3);
02649
02650 d(1)=0.0; d(2)=0.0; d(3)=0.0; a(1)=0.0000; a(2)=0.0000; a(3)=0.0000;
02651
02652
02653
02654
02655
02656
02657
02658
02659
02660
02661
02662
02663
02664 MuGlShift = d; MuGlAngle = a;
02665
02666 if((d(1) == 0.) & (d(2) == 0.) & (d(3) == 0.) &
02667 (a(1) == 0.) & (a(2) == 0.) & (a(3) == 0.)){
02668 Rm = GRm;
02669 Pm = GPm;
02670 AlgebraicVector4 Vm0;
02671 Vm = AlgebraicVector4(tsosMuon.localParameters().vector()(1),
02672 tsosMuon.localParameters().vector()(2),
02673 tsosMuon.localParameters().vector()(3),
02674 tsosMuon.localParameters().vector()(4));
02675 if(debug_)
02676 std::cout<<" ...... without misalignment "<<std::endl;
02677 return;
02678 }
02679
02680 Rm0(1) = GRm.x(); Rm0(2) = GRm.y(); Rm0(3) = GRm.z();
02681 Pm0(1) = GPm.x(); Pm0(2) = GPm.y(); Pm0(3) = GPm.z();
02682 Norm(1) = Nl.x(); Norm(2) = Nl.y(); Norm(3) = Nl.z();
02683
02684 CLHEP::HepMatrix T(3,3,1);
02685
02686
02687
02688
02689
02690 double s1 = std::sin(a(1)), c1 = std::cos(a(1));
02691 double s2 = std::sin(a(2)), c2 = std::cos(a(2));
02692 double s3 = std::sin(a(3)), c3 = std::cos(a(3));
02693 T(1,1) = c2 * c3;
02694 T(1,2) = c1 * s3 + s1 * s2 * c3;
02695 T(1,3) = s1 * s3 - c1 * s2 * c3;
02696 T(2,1) = -c2 * s3;
02697 T(2,2) = c1 * c3 - s1 * s2 * s3;
02698 T(2,3) = s1 * c3 + c1 * s2 * s3;
02699 T(3,1) = s2;
02700 T(3,2) = -s1 * c2;
02701 T(3,3) = c1 * c2;
02702
02703
02704
02705
02706 CLHEP::HepMatrix Ti = T.T();
02707 CLHEP::HepVector di = -Ti * d;
02708
02709 CLHEP::HepVector ex0(3,0), ey0(3,0), ez0(3,0), ex(3,0), ey(3,0), ez(3,0);
02710 ex0(1) = 1.; ey0(2) = 1.; ez0(3) = 1.;
02711 ex = Ti*ex0; ey = Ti*ey0; ez = Ti*ez0;
02712
02713 CLHEP::HepVector TiN = Ti * Norm;
02714
02715
02716
02717
02718 double si = CLHEP_dot((Ti*Rm0 - Rm0 + di), TiN) / CLHEP_dot(Pm0, TiN);
02719 Rm1(1) = CLHEP_dot(ex, Rm0 + si*Pm0 - di);
02720 Rm1(2) = CLHEP_dot(ey, Rm0 + si*Pm0 - di);
02721 Rm1(3) = CLHEP_dot(ez, Rm0 + si*Pm0 - di);
02722 Pm1 = T * Pm0;
02723
02724
02725 Rm = GlobalVector(Rm1(1), Rm1(2), Rm1(3));
02726
02727 Pm = GlobalVector(CLHEP_dot(Pm0,ex), CLHEP_dot(Pm0, ey), CLHEP_dot(Pm0,ez));
02728
02729
02730 const Surface& refSurface = tsosMuon.surface();
02731 CLHEP::HepMatrix Tl(3,3,0);
02732
02733 Tl(1,1) = refSurface.rotation().xx();
02734 Tl(1,2) = refSurface.rotation().xy();
02735 Tl(1,3) = refSurface.rotation().xz();
02736
02737 Tl(2,1) = refSurface.rotation().yx();
02738 Tl(2,2) = refSurface.rotation().yy();
02739 Tl(2,3) = refSurface.rotation().yz();
02740
02741 Tl(3,1) = refSurface.rotation().zx();
02742 Tl(3,2) = refSurface.rotation().zy();
02743 Tl(3,3) = refSurface.rotation().zz();
02744
02745 CLHEP::HepVector R0(3,0), newR0(3,0), newRl(3,0), newPl(3,0);
02746 R0(1) = refSurface.position().x();
02747 R0(2) = refSurface.position().y();
02748 R0(3) = refSurface.position().z();
02749
02750 newRl = Tl * (Rm1 - R0);
02751 newPl = Tl * Pm1;
02752
02753 Vm(0) = newPl(1)/newPl(3);
02754 Vm(1) = newPl(2)/newPl(3);
02755 Vm(2) = newRl(1);
02756 Vm(3) = newRl(2);
02757
02758 #ifdef CHECK_DERIVATIVES_LOCAL_VS_ANGLES
02759 float del = 0.0001;
02760
02761 int ii = 2; T(3,1) -=del; T(1,3) +=del;
02762
02763 AlgebraicVector4 Vm0 = Vm;
02764 CLHEP::HepVector Rm10 = Rm1, Pm10 = Pm1;;
02765 CLHEP::HepMatrix newTl = Tl*T;
02766 Ti = T.T();
02767 di = -Ti * d;
02768 ex = Ti*ex0; ey = Ti*ey0; ez = Ti*ez0;
02769 TiN = Ti * Norm;
02770 si = CLHEP_dot((Ti*Rm0 - Rm0 + di), TiN) / CLHEP_dot(Pm0, TiN);
02771 Rm1(1) = CLHEP_dot(ex, Rm0 + si*Pm0 - di);
02772 Rm1(2) = CLHEP_dot(ey, Rm0 + si*Pm0 - di);
02773 Rm1(3) = CLHEP_dot(ez, Rm0 + si*Pm0 - di);
02774 Pm1 = T * Pm0;
02775
02776 newRl = Tl * (Rm1 - R0);
02777 newPl = Tl * Pm1;
02778
02779 Vm(0) = newPl(1)/newPl(3);
02780 Vm(1) = newPl(2)/newPl(3);
02781 Vm(2) = newRl(1);
02782 Vm(3) = newRl(2);
02783 std::cout<<" ========= dVm/da"<<ii<<" "<<(Vm-Vm0)*(1./del)<<std::endl;
02784 #endif
02785
02786 if(debug_){
02787
02788
02789 std::cout<<" ex "<<(Tl.T()*ex0).T()<<std::endl;
02790 std::cout<<" ey "<<(Tl.T()*ey0).T()<<std::endl;
02791 std::cout<<" ez "<<(Tl.T()*ez0).T()<<std::endl;
02792
02793
02794 std::cout<<" pxyz/p gl 0 "<<(Pm0/Pm0.norm()).T()<<std::endl;
02795 std::cout<<" pxyz/p loc0 "<<(Tl*Pm0/(Tl*Pm0).norm()).T()<<std::endl;
02796 std::cout<<" pxyz/p glob "<<(Pm1/Pm1.norm()).T()<<std::endl;
02797 std::cout<<" pxyz/p loc "<<(newPl/newPl.norm()).T()<<std::endl;
02798
02799
02800
02801 std::cout<<" ---- phi g0 g1 l0 l1 "
02802 <<atan2(Pm0(2),Pm0(1))<<" "<<atan2(Pm1(2),Pm1(1))<<" "
02803 <<atan2((Tl*Pm0)(2),(Tl*Pm0)(1))<<" "
02804 <<atan2(newPl(2),newPl(1))<<std::endl;
02805 std::cout<<" ---- angl Gl Loc "<<atan2(Pm1(2),Pm1(1))-atan2(Pm0(2),Pm0(1))<<" "
02806 <<atan2(newPl(2),newPl(1))-atan2((Tl*Pm0)(2),(Tl*Pm0)(1))<<" "
02807 <<atan2(newPl(3),newPl(2))-atan2((Tl*Pm0)(3),(Tl*Pm0)(2))<<" "
02808 <<atan2(newPl(1),newPl(3))-atan2((Tl*Pm0)(1),(Tl*Pm0)(3))<<" "<<std::endl;
02809 }
02810
02811 if(debug_){
02812 CLHEP::HepMatrix newTl(3,3,0);
02813 CLHEP::HepVector R0(3,0), newRl(3,0), newPl(3,0);
02814 newTl = Tl * Ti.T();
02815 newR0 = Ti * R0 + di;
02816
02817 std::cout<<" N "<<Norm.T()<<std::endl;
02818 std::cout<<" dtxl yl "<<Vm(0)-tsosMuon.localParameters().vector()(1)<<" "
02819 <<Vm(1)-tsosMuon.localParameters().vector()(2)<<std::endl;
02820 std::cout<<" dXl dYl "<<Vm(2)-tsosMuon.localParameters().vector()(3)<<" "
02821 <<Vm(3)-tsosMuon.localParameters().vector()(4)<<std::endl;
02822 std::cout<<" localM "<<tsosMuon.localParameters().vector()<<std::endl;
02823 std::cout<<" Vm "<<Vm<<std::endl;
02824
02825
02826 CLHEP::HepVector Rvc(3,0), Pvc(3,0), Rvg(3,0), Pvg(3,0);
02827 Rvc(1) = Vm(2); Rvc(2) = Vm(3);
02828 Pvc(3) = tsosMuon.localParameters().pzSign() * Pm0.norm() /
02829 sqrt(1 + Vm(0)*Vm(0) + Vm(1)*Vm(1));
02830 Pvc(1) = Pvc(3) * Vm(0);
02831 Pvc(2) = Pvc(3) * Vm(1);
02832
02833 Rvg = newTl.T() * Rvc + newR0;
02834 Pvg = newTl.T() * Pvc;
02835
02836 std::cout<<" newPl "<<newPl.T()<<std::endl;
02837 std::cout<<" Pvc "<<Pvc.T()<<std::endl;
02838 std::cout<<" Rvg "<<Rvg.T()<<std::endl;
02839 std::cout<<" Rm1 "<<Rm1.T()<<std::endl;
02840 std::cout<<" Pvg "<<Pvg.T()<<std::endl;
02841 std::cout<<" Pm1 "<<Pm1.T()<<std::endl;
02842 }
02843
02844 if(debug_){
02845
02846 std::cout<<" ----- Pm "<<Pm<<std::endl;
02847 std::cout<<" T*Pm0 "<<Pm1.T()<<std::endl;
02848
02849
02850
02851 CLHEP::HepVector Rml = Rm1(1)*ex + Rm1(2)*ey + Rm1(3)*ez + di;
02852
02853 CLHEP::HepVector Rt0(3);
02854 Rt0(1)=Rt.x(); Rt0(2)=Rt.y(); Rt0(3)=Rt.z();
02855
02856
02857 double sl = CLHEP_dot(T*(Rt0 - Rml), T*Norm) / CLHEP_dot(Ti * Pm1, Norm);
02858 CLHEP::HepVector Rl = Rml + sl*(Ti*Pm1);
02859
02860
02861
02862
02863 double A = CLHEP_dot(Ti*Pm1, Norm);
02864 double B = CLHEP_dot(Rt0, Norm) - CLHEP_dot(Ti*Rm1, Norm)
02865 + CLHEP_dot(Ti*d, Norm);
02866 double s = B/A;
02867 CLHEP::HepVector Dr = Ti*Rm1 - Ti*d + s*(Ti*Pm1) - Rm0;
02868 CLHEP::HepVector Dp = Ti*Pm1 - Pm0;
02869
02870 CLHEP::HepVector TiR = Ti * Rm0 + di;
02871 CLHEP::HepVector Ri = T * TiR + d;
02872
02873 std::cout<<" --TTi-0 "<<(Ri-Rm0).T()<<std::endl;
02874 std::cout<<" -- Dr "<<Dr.T()<<endl;
02875 std::cout<<" -- Dp "<<Dp.T()<<endl;
02876
02877
02878
02879
02880
02881
02882
02883 std::cout<<" -- si sl s "<<si<<" "<<sl<<" "<<s<<std::endl;
02884
02885
02886 std::cout<<" -- Rml-(Rm0+si*Pm0) "<<(Rml - Rm0 - si*Pm0).T()<<std::endl;
02887
02888
02889
02890
02891
02892
02893
02894 std::cout<<" --Rl-0 "<<(Rl-Rm0).T()<<std::endl;
02895
02896 }
02897
02898 return;
02899
02900 }
02901
02902
02903
02904 void
02905 GlobalTrackerMuonAlignment::trackFitter(reco::TrackRef alongTr, reco::TransientTrack& alongTTr,
02906 PropagationDirection direction,
02907 TrajectoryStateOnSurface& trackFittedTSOS)
02908 {
02909 bool info = false;
02910 bool smooth = false;
02911
02912 if(info)
02913 std::cout<<" ......... start of trackFitter ......... "<<std::endl;
02914 if(false && info){
02915 this->debugTrackHit(" Tracker track hits ", alongTr);
02916 this->debugTrackHit(" Tracker TransientTrack hits ", alongTTr);
02917 }
02918
02919 edm::OwnVector<TrackingRecHit> recHit;
02920 DetId trackDetId = DetId(alongTr->innerDetId());
02921 for (trackingRecHit_iterator i=alongTr->recHitsBegin();i!=alongTr->recHitsEnd(); i++) {
02922 recHit.push_back((*i)->clone());
02923 }
02924 if(direction != alongMomentum)
02925 {
02926 edm::OwnVector<TrackingRecHit> recHitAlong = recHit;
02927 recHit.clear();
02928 for(unsigned int ihit = recHitAlong.size()-1; ihit+1>0; ihit--){
02929 recHit.push_back(recHitAlong[ihit]);
02930 }
02931 recHitAlong.clear();
02932 }
02933 if(info)
02934 std::cout<<" ... Own recHit.size() "<<recHit.size()<<" "<<trackDetId.rawId()<<std::endl;
02935
02936
02937 TransientTrackingRecHit::RecHitContainer recHitTrack;
02938 for (trackingRecHit_iterator i=alongTr->recHitsBegin();i!=alongTr->recHitsEnd(); i++) {
02939 if((*i)->isValid()){
02940 recHitTrack.push_back(TTRHBuilder->build(&**i ));}
02941 else{
02942 recHitTrack.push_back(TTRHBuilder->build(&**i ));}
02943 }
02944
02945 if(info)
02946 std::cout<<" ... recHitTrack.size() "<<recHit.size()<<" "<<trackDetId.rawId()
02947 <<" ======> recHitMu "<<std::endl;
02948
02949 MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu = recHitTrack;
02950
02951
02952
02953
02954 if(info)
02955 std::cout<<" ...along.... recHitMu.size() "<<recHitMu.size()<<std::endl;
02956 if(direction != alongMomentum){
02957 MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMuAlong = recHitMu;
02958 recHitMu.clear();
02959 for(unsigned int ihit = recHitMuAlong.size()-1; ihit+1>0; ihit--){
02960 recHitMu.push_back(recHitMuAlong[ihit]);
02961 }
02962 recHitMuAlong.clear();
02963 }
02964 if(info)
02965 std::cout<<" ...opposite... recHitMu.size() "<<recHitMu.size()<<std::endl;
02966
02967 TrajectoryStateOnSurface firstTSOS;
02968 if(direction == alongMomentum) firstTSOS = alongTTr.innermostMeasurementState();
02969 else firstTSOS = alongTTr.outermostMeasurementState();
02970
02971 AlgebraicSymMatrix55 CovLoc;
02972 CovLoc(0,0) = 1. * firstTSOS.localError().matrix()(0,0);
02973 CovLoc(1,1) = 10. * firstTSOS.localError().matrix()(1,1);
02974 CovLoc(2,2) = 10. * firstTSOS.localError().matrix()(2,2);
02975 CovLoc(3,3) = 100. * firstTSOS.localError().matrix()(3,3);
02976 CovLoc(4,4) = 100. * firstTSOS.localError().matrix()(4,4);
02977 TrajectoryStateOnSurface initialTSOS(firstTSOS.localParameters(), LocalTrajectoryError(CovLoc),
02978 firstTSOS.surface(), &*magneticField_);
02979
02980 PTrajectoryStateOnDet PTraj =
02981 trajectoryStateTransform::persistentState(initialTSOS, trackDetId.rawId());
02982
02983 const TrajectorySeed seedT(PTraj, recHit, direction);
02984
02985 std::vector<Trajectory> trajVec;
02986 if(direction == alongMomentum) trajVec = theFitter->fit(seedT, recHitMu, initialTSOS);
02987 else trajVec = theFitterOp->fit(seedT, recHitMu, initialTSOS);
02988
02989 if(info){
02990 this->debugTrajectorySOSv(" innermostTSOS of TransientTrack",
02991 alongTTr.innermostMeasurementState());
02992 this->debugTrajectorySOSv(" outermostTSOS of TransientTrack",
02993 alongTTr.outermostMeasurementState());
02994 this->debugTrajectorySOS(" initialTSOS for theFitter ", initialTSOS);
02995 std::cout<<" . . . . . trajVec.size() "<<trajVec.size()<<std::endl;
02996 if(trajVec.size()) this->debugTrajectory(" theFitter trajectory ", trajVec[0]);
02997 }
02998
02999 if(!smooth){
03000 if(trajVec.size()) trackFittedTSOS = trajVec[0].lastMeasurement().updatedState();}
03001 else{
03002 std::vector<Trajectory> trajSVec;
03003 if (trajVec.size()){
03004 if(direction == alongMomentum) trajSVec = theSmoother->trajectories(*(trajVec.begin()));
03005 else trajSVec = theSmootherOp->trajectories(*(trajVec.begin()));
03006 }
03007 if(info) std::cout<<" . . . . trajSVec.size() "<<trajSVec.size()<<std::endl;
03008 if(trajSVec.size()) this->debugTrajectorySOSv("smoothed geom InnermostState",
03009 trajSVec[0].geometricalInnermostState());
03010 if(trajSVec.size()) trackFittedTSOS = trajSVec[0].geometricalInnermostState();
03011 }
03012
03013 if(info) this->debugTrajectorySOSv(" trackFittedTSOS ", trackFittedTSOS);
03014
03015 }
03016
03017
03018 void
03019 GlobalTrackerMuonAlignment::muonFitter(reco::TrackRef alongTr, reco::TransientTrack& alongTTr,
03020 PropagationDirection direction,
03021 TrajectoryStateOnSurface& trackFittedTSOS)
03022 {
03023 bool info = false;
03024 bool smooth = false;
03025
03026 if(info) std::cout<<" ......... start of muonFitter ........"<<std::endl;
03027 if(false && info){
03028 this->debugTrackHit(" Muon track hits ", alongTr);
03029 this->debugTrackHit(" Muon TransientTrack hits ", alongTTr);
03030 }
03031
03032 edm::OwnVector<TrackingRecHit> recHit;
03033 DetId trackDetId = DetId(alongTr->innerDetId());
03034 for (trackingRecHit_iterator i=alongTr->recHitsBegin();i!=alongTr->recHitsEnd(); i++) {
03035 recHit.push_back((*i)->clone());
03036 }
03037 if(direction != alongMomentum)
03038 {
03039 edm::OwnVector<TrackingRecHit> recHitAlong = recHit;
03040 recHit.clear();
03041 for(unsigned int ihit = recHitAlong.size()-1; ihit+1>0; ihit--){
03042 recHit.push_back(recHitAlong[ihit]);
03043 }
03044 recHitAlong.clear();
03045 }
03046 if(info)
03047 std::cout<<" ... Own recHit.size() DetId==Muon "<<recHit.size()<<" "<<trackDetId.rawId()<<std::endl;
03048
03049 MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu =
03050 MuRHBuilder->build(alongTTr.recHitsBegin(), alongTTr.recHitsEnd());
03051 if(info)
03052 std::cout<<" ...along.... recHitMu.size() "<<recHitMu.size()<<std::endl;
03053 if(direction != alongMomentum){
03054 MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMuAlong = recHitMu;
03055 recHitMu.clear();
03056 for(unsigned int ihit = recHitMuAlong.size()-1; ihit+1>0; ihit--){
03057 recHitMu.push_back(recHitMuAlong[ihit]);
03058 }
03059 recHitMuAlong.clear();
03060 }
03061 if(info)
03062 std::cout<<" ...opposite... recHitMu.size() "<<recHitMu.size()<<std::endl;
03063
03064 TrajectoryStateOnSurface firstTSOS;
03065 if(direction == alongMomentum) firstTSOS = alongTTr.innermostMeasurementState();
03066 else firstTSOS = alongTTr.outermostMeasurementState();
03067
03068 AlgebraicSymMatrix55 CovLoc;
03069 CovLoc(0,0) = 1. * firstTSOS.localError().matrix()(0,0);
03070 CovLoc(1,1) = 10. * firstTSOS.localError().matrix()(1,1);
03071 CovLoc(2,2) = 10. * firstTSOS.localError().matrix()(2,2);
03072 CovLoc(3,3) = 100. * firstTSOS.localError().matrix()(3,3);
03073 CovLoc(4,4) = 100. * firstTSOS.localError().matrix()(4,4);
03074 TrajectoryStateOnSurface initialTSOS(firstTSOS.localParameters(), LocalTrajectoryError(CovLoc),
03075 firstTSOS.surface(), &*magneticField_);
03076
03077 PTrajectoryStateOnDet PTraj =
03078 trajectoryStateTransform::persistentState(initialTSOS, trackDetId.rawId());
03079
03080 const TrajectorySeed seedT(PTraj, recHit, direction);
03081
03082 std::vector<Trajectory> trajVec;
03083 if(direction == alongMomentum) trajVec = theFitter->fit(seedT, recHitMu, initialTSOS);
03084 else trajVec = theFitterOp->fit(seedT, recHitMu, initialTSOS);
03085
03086 if(info){
03087 this->debugTrajectorySOSv(" innermostTSOS of TransientTrack",
03088 alongTTr.innermostMeasurementState());
03089 this->debugTrajectorySOSv(" outermostTSOS of TransientTrack",
03090 alongTTr.outermostMeasurementState());
03091 this->debugTrajectorySOS(" initialTSOS for theFitter ", initialTSOS);
03092 std::cout<<" . . . . . trajVec.size() "<<trajVec.size()<<std::endl;
03093 if(trajVec.size()) this->debugTrajectory(" theFitter trajectory ", trajVec[0]);
03094 }
03095
03096 if(!smooth){
03097 if(trajVec.size()) trackFittedTSOS = trajVec[0].lastMeasurement().updatedState();}
03098 else{
03099 std::vector<Trajectory> trajSVec;
03100 if (trajVec.size()){
03101 if(direction == alongMomentum) trajSVec = theSmoother->trajectories(*(trajVec.begin()));
03102 else trajSVec = theSmootherOp->trajectories(*(trajVec.begin()));
03103 }
03104 if(info) std::cout<<" . . . . trajSVec.size() "<<trajSVec.size()<<std::endl;
03105 if(trajSVec.size()) this->debugTrajectorySOSv("smoothed geom InnermostState",
03106 trajSVec[0].geometricalInnermostState());
03107 if(trajSVec.size()) trackFittedTSOS = trajSVec[0].geometricalInnermostState();
03108 }
03109
03110 }
03111
03112
03113
03114 void GlobalTrackerMuonAlignment::debugTrackHit(const std::string title, TransientTrack& alongTr)
03115 {
03116 std::cout<<" ------- "<<title<<" --------"<<std::endl;
03117 int nHit = 1;
03118 TransientTrackingRecHit::RecHitContainer recHit;
03119 for (trackingRecHit_iterator i=alongTr.recHitsBegin();i!=alongTr.recHitsEnd(); i++) {
03120 std::cout<<" Hit "<<nHit++<<" DetId "<<(*i)->geographicalId().det();
03121 if( (*i)->geographicalId().det() == DetId::Tracker ) std::cout<<" Tracker ";
03122 else if ( (*i)->geographicalId().det() == DetId::Muon ) std::cout<<" Muon ";
03123 else std::cout<<" Unknown ";
03124 if(!(*i)->isValid()) std::cout<<" not valid "<<std::endl;
03125 else std::cout<<std::endl;
03126 }
03127 }
03128
03129
03130
03131 void GlobalTrackerMuonAlignment::debugTrackHit(const std::string title, reco::TrackRef alongTr)
03132 {
03133 std::cout<<" ------- "<<title<<" --------"<<std::endl;
03134 int nHit = 1;
03135 TransientTrackingRecHit::RecHitContainer recHit;
03136 for (trackingRecHit_iterator i=alongTr->recHitsBegin();i!=alongTr->recHitsEnd(); i++) {
03137 std::cout<<" Hit "<<nHit++<<" DetId "<<(*i)->geographicalId().det();
03138 if( (*i)->geographicalId().det() == DetId::Tracker ) std::cout<<" Tracker ";
03139 else if ( (*i)->geographicalId().det() == DetId::Muon ) std::cout<<" Muon ";
03140 else std::cout<<" Unknown ";
03141 if(!(*i)->isValid()) std::cout<<" not valid "<<std::endl;
03142 else std::cout<<std::endl;
03143 }
03144 }
03145
03146
03147 void GlobalTrackerMuonAlignment::debugTrajectorySOS(const std::string title,
03148 TrajectoryStateOnSurface& trajSOS)
03149 {
03150 std::cout<<" --- "<<title<<" --- "<<std::endl;
03151 if(!trajSOS.isValid()) {std::cout<<" Not valid !!!! "<<std::endl; return;}
03152 std::cout<<" R |p| "<<trajSOS.globalPosition().perp()<<" "
03153 <<trajSOS.globalMomentum().mag()<<" charge "<<trajSOS.charge()<<std::endl;
03154 std::cout<<" x p "<<trajSOS.globalParameters().vector()(0)<<" "
03155 <<trajSOS.globalParameters().vector()(1)<<" "
03156 <<trajSOS.globalParameters().vector()(2)<<" "
03157 <<trajSOS.globalParameters().vector()(3)<<" "
03158 <<trajSOS.globalParameters().vector()(4)<<" "
03159 <<trajSOS.globalParameters().vector()(5)<<std::endl;
03160 std::cout<<" +/- "<<sqrt(trajSOS.cartesianError().matrix()(0,0))<<" "
03161 <<sqrt(trajSOS.cartesianError().matrix()(1,1))<<" "
03162 <<sqrt(trajSOS.cartesianError().matrix()(2,2))<<" "
03163 <<sqrt(trajSOS.cartesianError().matrix()(3,3))<<" "
03164 <<sqrt(trajSOS.cartesianError().matrix()(4,4))<<" "
03165 <<sqrt(trajSOS.cartesianError().matrix()(5,5))<<std::endl;
03166 std::cout<<" q/p dxy/dz xy "<<trajSOS.localParameters().vector()(0)<<" "
03167 <<trajSOS.localParameters().vector()(1)<<" "
03168 <<trajSOS.localParameters().vector()(2)<<" "
03169 <<trajSOS.localParameters().vector()(3)<<" "
03170 <<trajSOS.localParameters().vector()(4)<<std::endl;
03171 std::cout<<" +/- error "<<sqrt(trajSOS.localError().matrix()(0,0))<<" "
03172 <<sqrt(trajSOS.localError().matrix()(1,1))<<" "
03173 <<sqrt(trajSOS.localError().matrix()(2,2))<<" "
03174 <<sqrt(trajSOS.localError().matrix()(3,3))<<" "
03175 <<sqrt(trajSOS.localError().matrix()(4,4))<<" "<<std::endl;
03176 }
03177
03178
03179 void GlobalTrackerMuonAlignment::debugTrajectorySOSv(const std::string title,
03180 TrajectoryStateOnSurface trajSOS)
03181 {
03182 std::cout<<" --- "<<title<<" --- "<<std::endl;
03183 if(!trajSOS.isValid()) {std::cout<<" Not valid !!!! "<<std::endl; return;}
03184 std::cout<<" R |p| "<<trajSOS.globalPosition().perp()<<" "
03185 <<trajSOS.globalMomentum().mag()<<" charge "<<trajSOS.charge()<<std::endl;
03186 std::cout<<" x p "<<trajSOS.globalParameters().vector()(0)<<" "
03187 <<trajSOS.globalParameters().vector()(1)<<" "
03188 <<trajSOS.globalParameters().vector()(2)<<" "
03189 <<trajSOS.globalParameters().vector()(3)<<" "
03190 <<trajSOS.globalParameters().vector()(4)<<" "
03191 <<trajSOS.globalParameters().vector()(5)<<std::endl;
03192 std::cout<<" +/- "<<sqrt(trajSOS.cartesianError().matrix()(0,0))<<" "
03193 <<sqrt(trajSOS.cartesianError().matrix()(1,1))<<" "
03194 <<sqrt(trajSOS.cartesianError().matrix()(2,2))<<" "
03195 <<sqrt(trajSOS.cartesianError().matrix()(3,3))<<" "
03196 <<sqrt(trajSOS.cartesianError().matrix()(4,4))<<" "
03197 <<sqrt(trajSOS.cartesianError().matrix()(5,5))<<std::endl;
03198 std::cout<<" q/p dxy/dz xy "<<trajSOS.localParameters().vector()(0)<<" "
03199 <<trajSOS.localParameters().vector()(1)<<" "
03200 <<trajSOS.localParameters().vector()(2)<<" "
03201 <<trajSOS.localParameters().vector()(3)<<" "
03202 <<trajSOS.localParameters().vector()(4)<<std::endl;
03203 std::cout<<" +/- error "<<sqrt(trajSOS.localError().matrix()(0,0))<<" "
03204 <<sqrt(trajSOS.localError().matrix()(1,1))<<" "
03205 <<sqrt(trajSOS.localError().matrix()(2,2))<<" "
03206 <<sqrt(trajSOS.localError().matrix()(3,3))<<" "
03207 <<sqrt(trajSOS.localError().matrix()(4,4))<<" "<<std::endl;
03208 }
03209
03210
03211 void GlobalTrackerMuonAlignment::debugTrajectory(const std::string title, Trajectory& traj)
03212 {
03213 std::cout<<"\n"<<" ...... "<<title<<" ...... "<<std::endl;
03214 if(!traj.isValid()) {std::cout<<" Not valid !!!!!!!! "<<std::endl; return;}
03215 std::cout<<" chi2/Nhit "<<traj.chiSquared()<<" / "<<traj.foundHits();
03216 if(traj.direction() == alongMomentum) std::cout<<" alongMomentum >>>>"<<std::endl;
03217 else std::cout<<" oppositeToMomentum <<<<"<<std::endl;
03218 this->debugTrajectorySOSv(" firstMeasurementTSOS ",traj.firstMeasurement().updatedState());
03219
03220 this->debugTrajectorySOSv(" lastMeasurementTSOS ",traj.lastMeasurement().updatedState());
03221
03222 std::cout<<" . . . . . . . . . . . . . . . . . . . . . . . . . . . . \n"<<std::endl;
03223 }
03224
03225
03226
03227 void GlobalTrackerMuonAlignment::writeGlPosRcd(CLHEP::HepVector& paramVec)
03228 {
03229
03230
03231
03232
03233
03234
03235
03236
03237
03238
03239
03240 std::cout<<" paramVector "<<paramVec.T()<<std::endl;
03241
03242 CLHEP::Hep3Vector colX, colY, colZ;
03243
03244 #ifdef NOT_EXACT_ROTATION_MATRIX
03245 colX = CLHEP::Hep3Vector( 1., -paramVec(6), paramVec(5));
03246 colY = CLHEP::Hep3Vector( paramVec(6), 1., -paramVec(4));
03247 colZ = CLHEP::Hep3Vector(-paramVec(5), paramVec(4), 1.);
03248 #else
03249 double s1 = std::sin(paramVec(4)), c1 = std::cos(paramVec(4));
03250 double s2 = std::sin(paramVec(5)), c2 = std::cos(paramVec(5));
03251 double s3 = std::sin(paramVec(6)), c3 = std::cos(paramVec(6));
03252 colX = CLHEP::Hep3Vector( c2 * c3, -c2 * s3, s2);
03253 colY = CLHEP::Hep3Vector(c1 * s3 + s1 * s2 * c3, c1 * c3 - s1 * s2 * s3, -s1 * c2);
03254 colZ = CLHEP::Hep3Vector(s1 * s3 - c1 * s2 * c3, s1 * c3 + c1 * s2 * s3, c1 * c2);
03255 #endif
03256
03257 CLHEP::HepVector param0(6,0);
03258
03259 Alignments* globalPositions = new Alignments();
03260
03261
03262 AlignTransform tracker(iteratorTrackerRcd->translation(),
03263 iteratorTrackerRcd->rotation(),
03264 DetId(DetId::Tracker).rawId());
03265
03266 CLHEP::Hep3Vector posMuGlRcd = iteratorMuonRcd->translation();
03267 CLHEP::HepRotation rotMuGlRcd = iteratorMuonRcd->rotation();
03268 CLHEP::HepEulerAngles angMuGlRcd = iteratorMuonRcd->rotation().eulerAngles();
03269 if(debug_)
03270 std::cout<<" Old muon Rcd Euler angles "<<angMuGlRcd.phi()<<" "<<angMuGlRcd.theta()
03271 <<" "<<angMuGlRcd.psi()<<std::endl;
03272 AlignTransform muon;
03273 if((angMuGlRcd.phi() == 0.) && (angMuGlRcd.theta() == 0.) && (angMuGlRcd.psi() == 0.) &&
03274 (posMuGlRcd.x() == 0.) && (posMuGlRcd.y() == 0.) && (posMuGlRcd.z() == 0.)){
03275 std::cout<<" New muon parameters are stored in Rcd"<<std::endl;
03276
03277 AlignTransform muonNew(AlignTransform::Translation(paramVec(1),
03278 paramVec(2),
03279 paramVec(3)),
03280 AlignTransform::Rotation (colX,
03281 colY,
03282 colZ),
03283 DetId(DetId::Muon).rawId());
03284 muon = muonNew;
03285 }
03286 else if((paramVec(1) == 0.) && (paramVec(2) == 0.) && (paramVec(3) == 0.) &&
03287 (paramVec(4) == 0.) && (paramVec(5) == 0.) && (paramVec(6) == 0.)){
03288 std::cout<<" Old muon parameters are stored in Rcd"<<std::endl;
03289
03290 AlignTransform muonNew(iteratorMuonRcd->translation(),
03291 iteratorMuonRcd->rotation(),
03292 DetId(DetId::Muon).rawId());
03293 muon = muonNew;
03294 }
03295 else{
03296 std::cout<<" New + Old muon parameters are stored in Rcd"<<std::endl;
03297 CLHEP::Hep3Vector posMuGlRcdThis = CLHEP::Hep3Vector(paramVec(1), paramVec(2), paramVec(3));
03298 CLHEP::HepRotation rotMuGlRcdThis = CLHEP::HepRotation(colX, colY, colZ);
03299 CLHEP::Hep3Vector posMuGlRcdNew = rotMuGlRcdThis * posMuGlRcd + posMuGlRcdThis;
03300 CLHEP::HepRotation rotMuGlRcdNew = rotMuGlRcdThis * rotMuGlRcd;
03301
03302 AlignTransform muonNew(posMuGlRcdNew,
03303 rotMuGlRcdNew,
03304 DetId(DetId::Muon).rawId());
03305 muon = muonNew;
03306 }
03307
03308
03309 AlignTransform ecal(iteratorEcalRcd->translation(),
03310 iteratorEcalRcd->rotation(),
03311 DetId(DetId::Ecal).rawId());
03312
03313 AlignTransform hcal(iteratorHcalRcd->translation(),
03314 iteratorHcalRcd->rotation(),
03315 DetId(DetId::Hcal).rawId());
03316
03317 AlignTransform calo(AlignTransform::Translation(param0(1),
03318 param0(2),
03319 param0(3)),
03320 AlignTransform::EulerAngles(param0(4),
03321 param0(5),
03322 param0(6)),
03323 DetId(DetId::Calo).rawId());
03324
03325 std::cout << "Tracker (" << tracker.rawId() << ") at " << tracker.translation()
03326 << " " << tracker.rotation().eulerAngles() << std::endl;
03327 std::cout << tracker.rotation() << std::endl;
03328
03329 std::cout << "Muon (" << muon.rawId() << ") at " << muon.translation()
03330 << " " << muon.rotation().eulerAngles() << std::endl;
03331 std::cout << " rotations angles around x,y,z "
03332 << " ( " << -muon.rotation().zy() << " " << muon.rotation().zx()
03333 << " " << -muon.rotation().yx() << " )" << std::endl;
03334 std::cout << muon.rotation() << std::endl;
03335
03336 std::cout << "Ecal (" << ecal.rawId() << ") at " << ecal.translation()
03337 << " " << ecal.rotation().eulerAngles() << std::endl;
03338 std::cout << ecal.rotation() << std::endl;
03339
03340 std::cout << "Hcal (" << hcal.rawId() << ") at " << hcal.translation()
03341 << " " << hcal.rotation().eulerAngles() << std::endl;
03342 std::cout << hcal.rotation() << std::endl;
03343
03344 std::cout << "Calo (" << calo.rawId() << ") at " << calo.translation()
03345 << " " << calo.rotation().eulerAngles() << std::endl;
03346 std::cout << calo.rotation() << std::endl;
03347
03348 globalPositions->m_align.push_back(tracker);
03349 globalPositions->m_align.push_back(muon);
03350 globalPositions->m_align.push_back(ecal);
03351 globalPositions->m_align.push_back(hcal);
03352 globalPositions->m_align.push_back(calo);
03353
03354 std::cout << "Uploading to the database..." << std::endl;
03355
03356 edm::Service<cond::service::PoolDBOutputService> poolDbService;
03357
03358 if (!poolDbService.isAvailable())
03359 throw cms::Exception("NotAvailable") << "PoolDBOutputService not available";
03360
03361
03362
03363
03364
03365
03366 poolDbService->writeOne<Alignments>(&(*globalPositions),
03367 poolDbService->currentTime(),
03368
03369 "GlobalPositionRcd");
03370 std::cout << "done!" << std::endl;
03371
03372 return;
03373 }
03374
03375
03376
03377 DEFINE_FWK_MODULE(GlobalTrackerMuonAlignment);
1.7.3