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/RKTestPropagator.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 RKTestPropagator alongRKPr =
00692 RKTestPropagator(magneticField_, alongMomentum);
00693 RKTestPropagator oppositeRKPr =
00694 RKTestPropagator(magneticField_, oppositeToMomentum);
00695
00696 float epsilon = 5.;
00697 SmartPropagator alongSmPr =
00698 SmartPropagator(alongRKPr, alongStHePr, magneticField_, alongMomentum, epsilon);
00699 SmartPropagator oppositeSmPr =
00700 SmartPropagator(oppositeRKPr, oppositeStHePr, magneticField_, oppositeToMomentum, epsilon);
00701
00702
00703
00704 for(MuonCollection::const_iterator itMuon = smuons->begin();
00705 itMuon != smuons->end(); ++itMuon) {
00706
00707 if(debug_){
00708 std::cout<<" mu gM is GM Mu SaM tM "<<itMuon->isGlobalMuon()<<" "
00709 <<itMuon->isMuon()<<" "<<itMuon->isStandAloneMuon()
00710 <<" "<<itMuon->isTrackerMuon()<<" "
00711 <<std::endl;
00712 }
00713
00714
00715 TransientTrack muTT(itMuon->outerTrack(), magneticField_, trackingGeometry_);
00716 TrajectoryStateOnSurface innerMuTSOS = muTT.innermostMeasurementState();
00717 TrajectoryStateOnSurface outerMuTSOS = muTT.outermostMeasurementState();
00718
00719
00720 TransientTrack trackTT(itMuon->track(), magneticField_, trackingGeometry_);
00721 TrajectoryStateOnSurface outerTrackTSOS = trackTT.outermostMeasurementState();
00722 TrajectoryStateOnSurface innerTrackTSOS = trackTT.innermostMeasurementState();
00723
00724 GlobalPoint pointTrackIn = innerTrackTSOS.globalPosition();
00725 GlobalPoint pointTrackOut = outerTrackTSOS.globalPosition();
00726 float lenghtTrack = (pointTrackOut-pointTrackIn).mag();
00727 GlobalPoint pointMuonIn = innerMuTSOS.globalPosition();
00728 GlobalPoint pointMuonOut = outerMuTSOS.globalPosition();
00729 float lenghtMuon = (pointMuonOut - pointMuonIn).mag();
00730 GlobalVector momentumTrackOut = outerTrackTSOS.globalMomentum();
00731 GlobalVector momentumTrackIn = innerTrackTSOS.globalMomentum();
00732 float distanceInIn = (pointTrackIn - pointMuonIn).mag();
00733 float distanceInOut = (pointTrackIn - pointMuonOut).mag();
00734 float distanceOutIn = (pointTrackOut - pointMuonIn).mag();
00735 float distanceOutOut = (pointTrackOut - pointMuonOut).mag();
00736
00737 if(debug_){
00738 std::cout<<" pointTrackIn "<<pointTrackIn<<std::endl;
00739 std::cout<<" Out "<<pointTrackOut<<" lenght "<<lenghtTrack<<std::endl;
00740 std::cout<<" point MuonIn "<<pointMuonIn<<std::endl;
00741 std::cout<<" Out "<<pointMuonOut<<" lenght "<<lenghtMuon<<std::endl;
00742 std::cout<<" momeTrackIn Out "<<momentumTrackIn<<" "<<momentumTrackOut<<std::endl;
00743 std::cout<<" doi io ii oo "<<distanceOutIn<<" "<<distanceInOut<<" "
00744 <<distanceInIn<<" "<<distanceOutOut<<std::endl;
00745 }
00746
00747 if(lenghtTrack < 90.) continue;
00748 if(lenghtMuon < 300.) continue;
00749 if(momentumTrackIn.mag() < 15. || momentumTrackOut.mag() < 15.) continue;
00750 if(trackTT.charge() != muTT.charge()) continue;
00751
00752 if(debug_)
00753 std::cout<<" passed lenght momentum cuts"<<std::endl;
00754
00755 GlobalVector GRm, GPm, Nl, Rm, Pm, Rt, Pt, Rt0;
00756 AlgebraicSymMatrix66 Cm, C0, Ce, C1;
00757
00758 TrajectoryStateOnSurface extrapolationT;
00759 int tsosMuonIf = 0;
00760
00761 if( isolatedMuonMode_ ){
00762 const Surface& refSurface = innerMuTSOS.surface();
00763 ReferenceCountingPointer<TangentPlane>
00764 tpMuLocal(refSurface.tangentPlane(innerMuTSOS.localPosition()));
00765 Nl = tpMuLocal->normalVector();
00766 ReferenceCountingPointer<TangentPlane>
00767 tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
00768 GlobalVector Ng = tpMuGlobal->normalVector();
00769 Surface* surf = (Surface*)&refSurface;
00770 const Plane* refPlane = dynamic_cast<Plane*>(surf);
00771 GlobalVector Nlp = refPlane->normalVector();
00772
00773 if(debug_){
00774 std::cout<<" Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
00775 <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
00776 std::cout<<" global "<<Ng.x()<<" "<<Ng.y()<<" "<<Ng.z()<<std::endl;
00777 std::cout<<" lp "<<Nlp.x()<<" "<<Nlp.y()<<" "<<Nlp.z()<<std::endl;
00778
00779
00780
00781 }
00782
00783
00784 extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
00785
00786
00787 if(!extrapolationT.isValid()){
00788 if(false & alarm)
00789 std::cout<<" !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
00790
00791 <<std::endl;
00792 continue;
00793 }
00794 tsosMuonIf = 1;
00795 Rt = GlobalVector((extrapolationT.globalPosition()).x(),
00796 (extrapolationT.globalPosition()).y(),
00797 (extrapolationT.globalPosition()).z());
00798
00799 Pt = extrapolationT.globalMomentum();
00800
00801 GRm = GlobalVector((innerMuTSOS.globalPosition()).x(),
00802 (innerMuTSOS.globalPosition()).y(),
00803 (innerMuTSOS.globalPosition()).z());
00804 GPm = innerMuTSOS.globalMomentum();
00805 Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
00806 (outerTrackTSOS.globalPosition()).y(),
00807 (outerTrackTSOS.globalPosition()).z());
00808 Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() +
00809 innerMuTSOS.cartesianError().matrix());
00810 C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
00811 Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
00812 C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
00813
00814 }
00815
00816
00817 if( cosmicMuonMode_ ){
00818
00819 if((distanceOutIn <= distanceInOut) & (distanceOutIn <= distanceInIn) &
00820 (distanceOutIn <= distanceOutOut)){
00821 if(debug_) std::cout<<" ----- Out - In -----"<<std::endl;
00822
00823 const Surface& refSurface = innerMuTSOS.surface();
00824 ReferenceCountingPointer<TangentPlane>
00825 tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
00826 Nl = tpMuGlobal->normalVector();
00827
00828
00829 extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
00830
00831
00832 if(!extrapolationT.isValid()){
00833 if(false & alarm)
00834 std::cout<<" !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
00835
00836 <<std::endl;
00837 continue;
00838 }
00839 if(debug_)
00840 std::cout<<" extrapolationT.isValid "<<extrapolationT.isValid()<<std::endl;
00841
00842 tsosMuonIf = 1;
00843 Rt = GlobalVector((extrapolationT.globalPosition()).x(),
00844 (extrapolationT.globalPosition()).y(),
00845 (extrapolationT.globalPosition()).z());
00846
00847 Pt = extrapolationT.globalMomentum();
00848
00849 GRm = GlobalVector((innerMuTSOS.globalPosition()).x(),
00850 (innerMuTSOS.globalPosition()).y(),
00851 (innerMuTSOS.globalPosition()).z());
00852 GPm = innerMuTSOS.globalMomentum();
00853 Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
00854 (outerTrackTSOS.globalPosition()).y(),
00855 (outerTrackTSOS.globalPosition()).z());
00856 Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() +
00857 innerMuTSOS.cartesianError().matrix());
00858 C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
00859 Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
00860 C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
00861
00862 if(false & debug_){
00863
00864 PropagationDirectionChooser Chooser;
00865 std::cout<<" propDirCh "
00866 <<Chooser.operator()(*outerTrackTSOS.freeState(), refSurface)
00867 <<" Ch == along "<<(alongMomentum ==
00868 Chooser.operator()(*outerTrackTSOS.freeState(), refSurface))
00869 <<std::endl;
00870 std::cout<<" --- Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
00871 <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
00872 }
00873 }
00874
00875
00876 else if((distanceInOut <= distanceInIn) & (distanceInOut <= distanceOutIn) &
00877 (distanceInOut <= distanceOutOut)){
00878 if(debug_) std::cout<<" ----- In - Out -----"<<std::endl;
00879
00880 const Surface& refSurface = outerMuTSOS.surface();
00881 ReferenceCountingPointer<TangentPlane>
00882 tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
00883 Nl = tpMuGlobal->normalVector();
00884
00885
00886 extrapolationT = oppositeSmPr.propagate(innerTrackTSOS, refSurface);
00887
00888 if(!extrapolationT.isValid()){
00889 if(false & alarm)
00890 std::cout<<" !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
00891 <<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
00892 <<std::endl;
00893 continue;
00894 }
00895 if(debug_)
00896 std::cout<<" extrapolationT.isValid "<<extrapolationT.isValid()<<std::endl;
00897
00898 tsosMuonIf = 2;
00899 Rt = GlobalVector((extrapolationT.globalPosition()).x(),
00900 (extrapolationT.globalPosition()).y(),
00901 (extrapolationT.globalPosition()).z());
00902
00903 Pt = extrapolationT.globalMomentum();
00904
00905 GRm = GlobalVector((outerMuTSOS.globalPosition()).x(),
00906 (outerMuTSOS.globalPosition()).y(),
00907 (outerMuTSOS.globalPosition()).z());
00908 GPm = outerMuTSOS.globalMomentum();
00909 Rt0 = GlobalVector((innerTrackTSOS.globalPosition()).x(),
00910 (innerTrackTSOS.globalPosition()).y(),
00911 (innerTrackTSOS.globalPosition()).z());
00912 Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() +
00913 outerMuTSOS.cartesianError().matrix());
00914 C0 = AlgebraicSymMatrix66(innerTrackTSOS.cartesianError().matrix());
00915 Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
00916 C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
00917
00918 if(false & debug_){
00919
00920 PropagationDirectionChooser Chooser;
00921 std::cout<<" propDirCh "
00922 <<Chooser.operator()(*innerTrackTSOS.freeState(), refSurface)
00923 <<" Ch == oppisite "<<(oppositeToMomentum ==
00924 Chooser.operator()(*innerTrackTSOS.freeState(), refSurface))
00925 <<std::endl;
00926 std::cout<<" --- Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
00927 <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
00928 }
00929 }
00930
00931 else if((distanceOutOut <= distanceInOut) & (distanceOutOut <= distanceInIn) &
00932 (distanceOutOut <= distanceOutIn)){
00933 if(debug_) std::cout<<" ----- Out - Out -----"<<std::endl;
00934
00935
00936 continue;
00937
00938 const Surface& refSurface = outerMuTSOS.surface();
00939 ReferenceCountingPointer<TangentPlane>
00940 tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
00941 Nl = tpMuGlobal->normalVector();
00942
00943
00944 extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
00945
00946 if(!extrapolationT.isValid()){
00947 if(alarm)
00948 std::cout<<" !!!!!!!!! Catastrophe Out-Out extrapolationT.isValid "
00949 <<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
00950 <<std::endl;
00951 continue;
00952 }
00953 if(debug_)
00954 std::cout<<" extrapolationT.isValid "<<extrapolationT.isValid()<<std::endl;
00955
00956 tsosMuonIf = 2;
00957 Rt = GlobalVector((extrapolationT.globalPosition()).x(),
00958 (extrapolationT.globalPosition()).y(),
00959 (extrapolationT.globalPosition()).z());
00960
00961 Pt = extrapolationT.globalMomentum();
00962
00963 GRm = GlobalVector((outerMuTSOS.globalPosition()).x(),
00964 (outerMuTSOS.globalPosition()).y(),
00965 (outerMuTSOS.globalPosition()).z());
00966 GPm = outerMuTSOS.globalMomentum();
00967 Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
00968 (outerTrackTSOS.globalPosition()).y(),
00969 (outerTrackTSOS.globalPosition()).z());
00970 Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() +
00971 outerMuTSOS.cartesianError().matrix());
00972 C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
00973 Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
00974 C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
00975
00976 if(debug_){
00977 PropagationDirectionChooser Chooser;
00978 std::cout<<" propDirCh "
00979 <<Chooser.operator()(*outerTrackTSOS.freeState(), refSurface)
00980 <<" Ch == along "<<(alongMomentum ==
00981 Chooser.operator()(*outerTrackTSOS.freeState(), refSurface))
00982 <<std::endl;
00983 std::cout<<" --- Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
00984 <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
00985 std::cout<<" Nornal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
00986 <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
00987 }
00988 }
00989 else{
00990 if(alarm)
00991 std::cout<<" ------- !!!!!!!!!!!!!!! No proper Out - In \a\a\a\a\a\a\a"<<std::endl;
00992 continue;
00993 }
00994
00995 }
00996
00997 if(tsosMuonIf == 0) {if(info) {std::cout<<"No tsosMuon !!!!!!"<<std::endl; continue;}}
00998 TrajectoryStateOnSurface tsosMuon;
00999 if(tsosMuonIf == 1) tsosMuon = muTT.innermostMeasurementState();
01000 else tsosMuon = muTT.outermostMeasurementState();
01001
01002
01003 AlgebraicVector4 LPRm;
01004 GlobalTrackerMuonAlignment::misalignMuonL
01005 (GRm, GPm, Nl, Rt, Rm, Pm, LPRm, extrapolationT, tsosMuon);
01006
01007 GlobalVector resR = Rm - Rt;
01008 GlobalVector resP0 = Pm - Pt;
01009 GlobalVector resP = Pm / Pm.mag() - Pt / Pt.mag();
01010 float RelMomResidual = (Pm.mag() - Pt.mag()) / (Pt.mag() + 1.e-6);;
01011
01012 AlgebraicVector6 Vm;
01013 Vm(0) = resR.x(); Vm(1) = resR.y(); Vm(2) = resR.z();
01014 Vm(3) = resP0.x(); Vm(4) = resP0.y(); Vm(5) = resP0.z();
01015 float Rmuon = Rm.perp();
01016 float Zmuon = Rm.z();
01017 float alfa_x = atan2(Nl.x(),Nl.y())*57.29578;
01018
01019 if(debug_){
01020 std::cout<<" Nx Ny Nz alfa_x "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "<<alfa_x<<std::endl;
01021
01022
01023 }
01024
01025 double chi_d = 0;
01026 for(int i=0; i<=5; i++) chi_d += Vm(i)*Vm(i)/Cm(i,i);
01027
01028 AlgebraicVector5 Vml(tsosMuon.localParameters().vector()
01029 - extrapolationT.localParameters().vector());
01030 AlgebraicSymMatrix55 m(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
01031 AlgebraicSymMatrix55 Cml(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
01032 bool ierrLoc = !m.Invert();
01033 if (ierrLoc && debug_ && info) {
01034 std::cout<< " ==== Error inverting Local covariance matrix ==== "<<std::cout;
01035 continue;}
01036 double chi_Loc = ROOT::Math::Similarity(Vml,m);
01037 if(debug_)
01038 std::cout<<" chi_Loc px/pz/err "<<chi_Loc<<" "<<Vml(1)/sqrt(Cml(1,1))<<std::endl;
01039
01040 if(Pt.mag() < 15.) continue;
01041 if(Pm.mag() < 5.) continue;
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
01084
01085 if(Rmuon < 120. || Rmuon > 450.) continue;
01086 if(Zmuon < -720. || Zmuon > 720.) continue;
01087 MuSelect = " Barrel+EndCaps";
01088
01089
01090 if(debug_)
01091 std::cout<<" .............. passed all cuts"<<std::endl;
01092
01093 N_track++;
01094
01095
01096 GlobalTrackerMuonAlignment::gradientGlobalAlg(Rt, Pt, Rm, Nl, Cm);
01097 GlobalTrackerMuonAlignment::gradientGlobal(Rt, Pt, Rm, Pm, Nl, Cm);
01098
01099 CLHEP::HepSymMatrix covLoc(4,0);
01100 for(int i=1; i<=4; i++)
01101 for(int j=1; j<=i; j++){
01102 covLoc(i,j) = (tsosMuon.localError().matrix()
01103 + extrapolationT.localError().matrix())(i,j);
01104
01105 }
01106
01107 const Surface& refSurface = tsosMuon.surface();
01108 CLHEP::HepMatrix rotLoc (3,3,0);
01109 rotLoc(1,1) = refSurface.rotation().xx();
01110 rotLoc(1,2) = refSurface.rotation().xy();
01111 rotLoc(1,3) = refSurface.rotation().xz();
01112
01113 rotLoc(2,1) = refSurface.rotation().yx();
01114 rotLoc(2,2) = refSurface.rotation().yy();
01115 rotLoc(2,3) = refSurface.rotation().yz();
01116
01117 rotLoc(3,1) = refSurface.rotation().zx();
01118 rotLoc(3,2) = refSurface.rotation().zy();
01119 rotLoc(3,3) = refSurface.rotation().zz();
01120
01121 CLHEP::HepVector posLoc(3);
01122 posLoc(1) = refSurface.position().x();
01123 posLoc(2) = refSurface.position().y();
01124 posLoc(3) = refSurface.position().z();
01125
01126 GlobalTrackerMuonAlignment::gradientLocal(Rt, Pt, Rm, Pm, Nl, covLoc, rotLoc, posLoc, LPRm);
01127
01128 if(debug_){
01129 std::cout<<" Norm "<<Nl<<std::endl;
01130 std::cout<<" posLoc "<<posLoc.T()<<std::endl;
01131 std::cout<<" rotLoc "<<rotLoc<<std::endl;
01132 }
01133
01134
01135 histo->Fill(itMuon->track()->pt());
01136
01137
01138 histo2->Fill(Pt.mag());
01139 histo3->Fill((PI/2.-itMuon->track()->outerTheta()));
01140 histo4->Fill(itMuon->track()->phi());
01141 histo5->Fill(Rmuon);
01142 histo6->Fill(Zmuon);
01143 histo7->Fill(RelMomResidual);
01144
01145 histo8->Fill(chi_d);
01146
01147 histo101->Fill(Zmuon, Rmuon);
01148 histo101->Fill(Rt0.z(), Rt0.perp());
01149 histo102->Fill(Rt0.x(), Rt0.y());
01150 histo102->Fill(Rm.x(), Rm.y());
01151
01152 histo11->Fill(resR.mag());
01153 if(fabs(Nl.x()) < 0.98) histo12->Fill(resR.x());
01154 if(fabs(Nl.y()) < 0.98) histo13->Fill(resR.y());
01155 if(fabs(Nl.z()) < 0.98) histo14->Fill(resR.z());
01156 histo15->Fill(resP.x());
01157 histo16->Fill(resP.y());
01158 histo17->Fill(resP.z());
01159
01160 if((fabs(Nl.x()) < 0.98) && (fabs(Nl.y()) < 0.98) &&(fabs(Nl.z()) < 0.98))
01161 {
01162 histo18->Fill(sqrt(C0(0,0)));
01163 histo19->Fill(sqrt(C1(0,0)));
01164 histo20->Fill(sqrt(C1(0,0)+Ce(0,0)));
01165 }
01166 if(fabs(Nl.x()) < 0.98) histo21->Fill(Vm(0)/sqrt(Cm(0,0)));
01167 if(fabs(Nl.y()) < 0.98) histo22->Fill(Vm(1)/sqrt(Cm(1,1)));
01168 if(fabs(Nl.z()) < 0.98) histo23->Fill(Vm(2)/sqrt(Cm(2,2)));
01169 histo24->Fill(Vm(3)/sqrt(C1(3,3)+Ce(3,3)));
01170 histo25->Fill(Vm(4)/sqrt(C1(4,4)+Ce(4,4)));
01171 histo26->Fill(Vm(5)/sqrt(C1(5,5)+Ce(5,5)));
01172 histo27->Fill(Nl.x());
01173 histo28->Fill(Nl.y());
01174 histo29->Fill(lenghtTrack);
01175 histo30->Fill(lenghtMuon);
01176 histo31->Fill(chi_Loc);
01177 histo32->Fill(Vml(1)/sqrt(Cml(1,1)));
01178 histo33->Fill(Vml(2)/sqrt(Cml(2,2)));
01179 histo34->Fill(Vml(3)/sqrt(Cml(3,3)));
01180 histo35->Fill(Vml(4)/sqrt(Cml(4,4)));
01181
01182 if (debug_) {
01183
01184 std::cout<<" diag 0[ "<<C0(0,0)<<" "<<C0(1,1)<<" "<<C0(2,2)<<" "<<C0(3,3)<<" "
01185 <<C0(4,4)<<" "<<C0(5,5)<<" ]"<<std::endl;
01186 std::cout<<" diag e[ "<<Ce(0,0)<<" "<<Ce(1,1)<<" "<<Ce(2,2)<<" "<<Ce(3,3)<<" "
01187 <<Ce(4,4)<<" "<<Ce(5,5)<<" ]"<<std::endl;
01188 std::cout<<" diag 1[ "<<C1(0,0)<<" "<<C1(1,1)<<" "<<C1(2,2)<<" "<<C1(3,3)<<" "
01189 <<C1(4,4)<<" "<<C1(5,5)<<" ]"<<std::endl;
01190 std::cout<<" Rm "<<Rm.x()<<" "<<Rm.y()<<" "<<Rm.z()
01191 <<" Pm "<<Pm.x()<<" "<<Pm.y()<<" "<<Pm.z()<<std::endl;
01192 std::cout<<" Rt "<<Rt.x()<<" "<<Rt.y()<<" "<<Rt.z()
01193 <<" Pt "<<Pt.x()<<" "<<Pt.y()<<" "<<Pt.z()<<std::endl;
01194 std::cout<<" Nl*(Rm-Rt) "<<Nl.dot(Rm-Rt)<<std::endl;
01195 std::cout<<" resR "<<resR.x()<<" "<<resR.y()<<" "<<resR.z()
01196 <<" resP "<<resP.x()<<" "<<resP.y()<<" "<<resP.z()<<std::endl;
01197 std::cout<<" Rm-t "<<(Rm-Rt).x()<<" "<<(Rm-Rt).y()<<" "<<(Rm-Rt).z()
01198 <<" Pm-t "<<(Pm-Pt).x()<<" "<<(Pm-Pt).y()<<" "<<(Pm-Pt).z()<<std::endl;
01199 std::cout<<" Vm "<<Vm<<std::endl;
01200 std::cout<<" +- "<<sqrt(Cm(0,0))<<" "<<sqrt(Cm(1,1))<<" "<<sqrt(Cm(2,2))<<" "
01201 <<sqrt(Cm(3,3))<<" "<<sqrt(Cm(4,4))<<" "<<sqrt(Cm(5,5))<<std::endl;
01202 std::cout<<" Pmuon Ptrack dP/Ptrack "<<itMuon->outerTrack()->p()<<" "
01203 <<itMuon->track()->outerP()<<" "
01204 <<(itMuon->outerTrack()->p() -
01205 itMuon->track()->outerP())/itMuon->track()->outerP()<<std::endl;
01206 std::cout<<" cov matrix "<<std::endl;
01207 std::cout<<Cm<<std::endl;
01208 std::cout<<" diag [ "<<Cm(0,0)<<" "<<Cm(1,1)<<" "<<Cm(2,2)<<" "<<Cm(3,3)<<" "
01209 <<Cm(4,4)<<" "<<Cm(5,5)<<" ]"<<std::endl;
01210
01211 static AlgebraicSymMatrix66 Ro;
01212 double Diag[6];
01213 for(int i=0; i<=5; i++) Diag[i] = sqrt(Cm(i,i));
01214 for(int i=0; i<=5; i++)
01215 for(int j=0; j<=5; j++)
01216 Ro(i,j) = Cm(i,j)/Diag[i]/Diag[j];
01217 std::cout<<" correlation matrix "<<std::endl;
01218 std::cout<<Ro<<std::endl;
01219
01220 AlgebraicSymMatrix66 CmI;
01221 for(int i=0; i<=5; i++)
01222 for(int j=0; j<=5; j++)
01223 CmI(i,j) = Cm(i,j);
01224
01225 bool ierr = !CmI.Invert();
01226 if( ierr ) { if(alarm || debug_)
01227 std::cout<<" Error inverse covariance matrix !!!!!!!!!!!"<<std::endl;
01228 continue;
01229 }
01230 std::cout<<" inverse cov matrix "<<std::endl;
01231 std::cout<<Cm<<std::endl;
01232
01233 double chi = ROOT::Math::Similarity(Vm, CmI);
01234 std::cout<<" chi chi_d "<<chi<<" "<<chi_d<<std::endl;
01235 }
01236
01237 }
01238
01239 }
01240
01241
01242
01243 void GlobalTrackerMuonAlignment::analyzeTrackTrajectory
01244 (const edm::Event& iEvent, const edm::EventSetup& iSetup)
01245 {
01246 using namespace edm;
01247 using namespace reco;
01248
01249 bool info = false;
01250 bool alarm = false;
01251
01252 double PI = 3.1415927;
01253
01254
01255 cosmicMuonMode_ = true;
01256
01257
01258 isolatedMuonMode_ = !cosmicMuonMode_;
01259
01260 N_event++;
01261 if (info || debug_) {
01262 std::cout << "----- event " << N_event << " -- tracks "<<N_track<<" ---";
01263 if (cosmicMuonMode_) std::cout << " treated as CosmicMu ";
01264 if (isolatedMuonMode_) std::cout <<" treated as IsolatedMu ";
01265 std::cout << std::endl;
01266 }
01267
01268 Handle<reco::TrackCollection> tracks;
01269 Handle<reco::TrackCollection> muons;
01270 Handle<reco::TrackCollection> gmuons;
01271 Handle<reco::MuonCollection> smuons;
01272
01273 if (collectionIsolated){
01274 iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "TrackerOnly", tracks);
01275 iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "StandAlone", muons);
01276 iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "GlobalMuon", gmuons);
01277 iEvent.getByLabel("ALCARECOMuAlCalIsolatedMu", "SelectedMuons", smuons);
01278 }
01279 else if (collectionCosmic){
01280 iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "TrackerOnly", tracks);
01281 iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "StandAlone", muons);
01282 iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "GlobalMuon", gmuons);
01283 iEvent.getByLabel("ALCARECOMuAlGlobalCosmics", "SelectedMuons", smuons);
01284 }
01285 else{
01286 iEvent.getByLabel(trackTags_,tracks);
01287 iEvent.getByLabel(muonTags_,muons);
01288 iEvent.getByLabel(gmuonTags_,gmuons);
01289 iEvent.getByLabel(smuonTags_,smuons);
01290 }
01291
01292 if (debug_) {
01293 std::cout << " ievBunch " << iEvent.bunchCrossing()
01294 << " runN " << (int)iEvent.run() <<std::endl;
01295 std::cout << " N tracks s/amu gmu selmu "<<tracks->size() << " " <<muons->size()
01296 << " "<<gmuons->size() << " " << smuons->size() << std::endl;
01297 for (MuonCollection::const_iterator itMuon = smuons->begin();
01298 itMuon != smuons->end();
01299 ++itMuon) {
01300 std::cout << " is isolatValid Matches " << itMuon->isIsolationValid()
01301 << " " <<itMuon->isMatchesValid() << std::endl;
01302 }
01303 }
01304
01305 if (isolatedMuonMode_) {
01306 if (tracks->size() != 1) return;
01307 if (muons->size() != 1) return;
01308 if (gmuons->size() != 1) return;
01309 if (smuons->size() != 1) return;
01310 }
01311
01312 if (cosmicMuonMode_){
01313 if (smuons->size() > 2) return;
01314 if (tracks->size() != smuons->size()) return;
01315 if (muons->size() != smuons->size()) return;
01316 if (gmuons->size() != smuons->size()) return;
01317 }
01318
01319
01320
01321
01322
01323
01324
01325
01326
01327
01328
01329 if (watchTrackingGeometry_.check(iSetup) || !trackingGeometry_) {
01330 edm::ESHandle<GlobalTrackingGeometry> trackingGeometry;
01331 iSetup.get<GlobalTrackingGeometryRecord>().get(trackingGeometry);
01332 trackingGeometry_ = &*trackingGeometry;
01333 theTrackingGeometry = trackingGeometry;
01334 }
01335
01336 if (watchMagneticFieldRecord_.check(iSetup) || !magneticField_) {
01337 edm::ESHandle<MagneticField> magneticField;
01338 iSetup.get<IdealMagneticFieldRecord>().get(magneticField);
01339 magneticField_ = &*magneticField;
01340 }
01341
01342 if (watchGlobalPositionRcd_.check(iSetup) || !globalPositionRcd_) {
01343 edm::ESHandle<Alignments> globalPositionRcd;
01344 iSetup.get<GlobalPositionRcd>().get(globalPositionRcd);
01345 globalPositionRcd_ = &*globalPositionRcd;
01346 for (std::vector<AlignTransform>::const_iterator i
01347 = globalPositionRcd_->m_align.begin();
01348 i != globalPositionRcd_->m_align.end(); ++i){
01349 if(DetId(DetId::Tracker).rawId() == i->rawId()) iteratorTrackerRcd = i;
01350 if(DetId(DetId::Muon).rawId() == i->rawId()) iteratorMuonRcd = i;
01351 if(DetId(DetId::Ecal).rawId() == i->rawId()) iteratorEcalRcd = i;
01352 if(DetId(DetId::Hcal).rawId() == i->rawId()) iteratorHcalRcd = i;
01353 }
01354 if(debug_){
01355 std::cout << "=== iteratorTrackerRcd " << iteratorTrackerRcd->rawId()<<" ====\n"
01356 << " translation " << iteratorTrackerRcd->translation()<<"\n"
01357 << " angles " << iteratorTrackerRcd->rotation().eulerAngles() << std::endl;
01358 std::cout << iteratorTrackerRcd->rotation() << std::endl;
01359 std::cout << "=== iteratorMuonRcd " << iteratorMuonRcd->rawId()<<" ====\n"
01360 << " translation " << iteratorMuonRcd->translation()<<"\n"
01361 << " angles " << iteratorMuonRcd->rotation().eulerAngles() << std::endl;
01362 std::cout << iteratorMuonRcd->rotation() << std::endl;
01363 }
01364 }
01365
01366 ESHandle<Propagator> propagator;
01367 iSetup.get<TrackingComponentsRecord>().get(propagator_, propagator);
01368
01369 SteppingHelixPropagator alongStHePr =
01370 SteppingHelixPropagator(magneticField_, alongMomentum);
01371 SteppingHelixPropagator oppositeStHePr =
01372 SteppingHelixPropagator(magneticField_, oppositeToMomentum);
01373
01374
01375 RKTestPropagator alongRKPr =
01376 RKTestPropagator(magneticField_, alongMomentum);
01377 RKTestPropagator oppositeRKPr =
01378 RKTestPropagator(magneticField_, oppositeToMomentum);
01379
01380 float epsilon = 5.;
01381 SmartPropagator alongSmPr =
01382 SmartPropagator(alongRKPr, alongStHePr, magneticField_, alongMomentum, epsilon);
01383 SmartPropagator oppositeSmPr =
01384 SmartPropagator(oppositeRKPr, oppositeStHePr, magneticField_, oppositeToMomentum, epsilon);
01385
01386 if(defineFitter){
01387 if(debug_)
01388 std::cout<<" ............... DEFINE FITTER ..................."<<std::endl;
01389 KFUpdator* theUpdator = new KFUpdator();
01390
01391 Chi2MeasurementEstimator* theEstimator = new Chi2MeasurementEstimator(100000,100000);
01392 theFitter = new KFTrajectoryFitter(alongSmPr,
01393 *theUpdator,
01394 *theEstimator);
01395 theSmoother = new KFTrajectorySmoother(alongSmPr,
01396 *theUpdator,
01397 *theEstimator);
01398 theFitterOp = new KFTrajectoryFitter(oppositeSmPr,
01399 *theUpdator,
01400 *theEstimator);
01401 theSmootherOp = new KFTrajectorySmoother(oppositeSmPr,
01402 *theUpdator,
01403 *theEstimator);
01404
01405 edm::ESHandle<TransientTrackingRecHitBuilder> builder;
01406 iSetup.get<TransientRecHitRecord>().get("WithTrackAngle",builder);
01407 this->TTRHBuilder = &(*builder);
01408 this->MuRHBuilder = new MuonTransientTrackingRecHitBuilder(theTrackingGeometry);
01409 if(debug_){
01410 std::cout<<" theTrackingGeometry.isValid() "<<theTrackingGeometry.isValid()<<std::endl;
01411 std::cout<< "get also the MuonTransientTrackingRecHitBuilder" << "\n";
01412 std::cout<< "get also the TransientTrackingRecHitBuilder" << "\n";
01413 }
01414 defineFitter = false;
01415 }
01416
01417
01418
01419 for(MuonCollection::const_iterator itMuon = smuons->begin();
01420 itMuon != smuons->end(); ++itMuon) {
01421
01422 if(debug_){
01423 std::cout<<" mu gM is GM Mu SaM tM "<<itMuon->isGlobalMuon()<<" "
01424 <<itMuon->isMuon()<<" "<<itMuon->isStandAloneMuon()
01425 <<" "<<itMuon->isTrackerMuon()<<" "
01426 <<std::endl;
01427 }
01428
01429
01430 TransientTrack muTT(itMuon->outerTrack(), magneticField_, trackingGeometry_);
01431 TrajectoryStateOnSurface innerMuTSOS = muTT.innermostMeasurementState();
01432 TrajectoryStateOnSurface outerMuTSOS = muTT.outermostMeasurementState();
01433
01434
01435 TransientTrack trackTT(itMuon->track(), magneticField_, trackingGeometry_);
01436 TrajectoryStateOnSurface outerTrackTSOS = trackTT.outermostMeasurementState();
01437 TrajectoryStateOnSurface innerTrackTSOS = trackTT.innermostMeasurementState();
01438
01439 GlobalPoint pointTrackIn = innerTrackTSOS.globalPosition();
01440 GlobalPoint pointTrackOut = outerTrackTSOS.globalPosition();
01441 float lenghtTrack = (pointTrackOut-pointTrackIn).mag();
01442 GlobalPoint pointMuonIn = innerMuTSOS.globalPosition();
01443 GlobalPoint pointMuonOut = outerMuTSOS.globalPosition();
01444 float lenghtMuon = (pointMuonOut - pointMuonIn).mag();
01445 GlobalVector momentumTrackOut = outerTrackTSOS.globalMomentum();
01446 GlobalVector momentumTrackIn = innerTrackTSOS.globalMomentum();
01447 float distanceInIn = (pointTrackIn - pointMuonIn).mag();
01448 float distanceInOut = (pointTrackIn - pointMuonOut).mag();
01449 float distanceOutIn = (pointTrackOut - pointMuonIn).mag();
01450 float distanceOutOut = (pointTrackOut - pointMuonOut).mag();
01451
01452 if(debug_){
01453 std::cout<<" pointTrackIn "<<pointTrackIn<<std::endl;
01454 std::cout<<" Out "<<pointTrackOut<<" lenght "<<lenghtTrack<<std::endl;
01455 std::cout<<" point MuonIn "<<pointMuonIn<<std::endl;
01456 std::cout<<" Out "<<pointMuonOut<<" lenght "<<lenghtMuon<<std::endl;
01457 std::cout<<" momeTrackIn Out "<<momentumTrackIn<<" "<<momentumTrackOut<<std::endl;
01458 std::cout<<" doi io ii oo "<<distanceOutIn<<" "<<distanceInOut<<" "
01459 <<distanceInIn<<" "<<distanceOutOut<<std::endl;
01460 }
01461
01462 if(lenghtTrack < 90.) continue;
01463 if(lenghtMuon < 300.) continue;
01464 if(innerMuTSOS.globalMomentum().mag() < 5. || outerMuTSOS.globalMomentum().mag() < 5.) continue;
01465 if(momentumTrackIn.mag() < 15. || momentumTrackOut.mag() < 15.) continue;
01466 if(trackTT.charge() != muTT.charge()) continue;
01467
01468 if(debug_)
01469 std::cout<<" passed lenght momentum cuts"<<std::endl;
01470
01471 GlobalVector GRm, GPm, Nl, Rm, Pm, Rt, Pt, Rt0;
01472 AlgebraicSymMatrix66 Cm, C0, Ce, C1;
01473
01474
01475 TrajectoryStateOnSurface extrapolationT;
01476 int tsosMuonIf = 0;
01477
01478 TrajectoryStateOnSurface muonFittedTSOS;
01479 TrajectoryStateOnSurface trackFittedTSOS;
01480
01481 if( isolatedMuonMode_ ){
01482 if(debug_) std::cout<<" ------ Isolated (out-in) ----- "<<std::endl;
01483 const Surface& refSurface = innerMuTSOS.surface();
01484 ReferenceCountingPointer<TangentPlane>
01485 tpMuLocal(refSurface.tangentPlane(innerMuTSOS.localPosition()));
01486 Nl = tpMuLocal->normalVector();
01487 ReferenceCountingPointer<TangentPlane>
01488 tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
01489 GlobalVector Ng = tpMuGlobal->normalVector();
01490 Surface* surf = (Surface*)&refSurface;
01491 const Plane* refPlane = dynamic_cast<Plane*>(surf);
01492 GlobalVector Nlp = refPlane->normalVector();
01493
01494 if(debug_){
01495 std::cout<<" Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
01496 <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
01497 std::cout<<" global "<<Ng.x()<<" "<<Ng.y()<<" "<<Ng.z()<<std::endl;
01498 std::cout<<" lp "<<Nlp.x()<<" "<<Nlp.y()<<" "<<Nlp.z()<<std::endl;
01499
01500
01501
01502 }
01503
01504
01505
01506
01507 if(!refitTrack_)
01508 extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
01509 else{
01510 GlobalTrackerMuonAlignment::trackFitter(itMuon->track(),
01511 trackTT,alongMomentum,trackFittedTSOS);
01512 if(trackFittedTSOS.isValid())
01513 extrapolationT = alongSmPr.propagate(trackFittedTSOS, refSurface);
01514 }
01515
01516 if(!extrapolationT.isValid()){
01517 if(false & alarm)
01518 std::cout<<" !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
01519
01520 <<std::endl;
01521 continue;
01522 }
01523 tsosMuonIf = 1;
01524 Rt = GlobalVector((extrapolationT.globalPosition()).x(),
01525 (extrapolationT.globalPosition()).y(),
01526 (extrapolationT.globalPosition()).z());
01527
01528 Pt = extrapolationT.globalMomentum();
01529
01530 GRm = GlobalVector((innerMuTSOS.globalPosition()).x(),
01531 (innerMuTSOS.globalPosition()).y(),
01532 (innerMuTSOS.globalPosition()).z());
01533 GPm = innerMuTSOS.globalMomentum();
01534
01535 Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
01536 (outerTrackTSOS.globalPosition()).y(),
01537 (outerTrackTSOS.globalPosition()).z());
01538 Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() +
01539 innerMuTSOS.cartesianError().matrix());
01540 C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
01541 Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
01542 C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
01543
01544 if(refitMuon_)
01545 GlobalTrackerMuonAlignment::muonFitter(itMuon->outerTrack(),muTT,oppositeToMomentum,muonFittedTSOS);
01546
01547 }
01548
01549
01550 if( cosmicMuonMode_ ){
01551
01552 if((distanceOutIn <= distanceInOut) & (distanceOutIn <= distanceInIn) &
01553 (distanceOutIn <= distanceOutOut)){
01554 if(debug_) std::cout<<" ----- Out - In -----"<<std::endl;
01555
01556 const Surface& refSurface = innerMuTSOS.surface();
01557 ReferenceCountingPointer<TangentPlane>
01558 tpMuGlobal(refSurface.tangentPlane(innerMuTSOS.globalPosition()));
01559 Nl = tpMuGlobal->normalVector();
01560
01561
01562
01563 if(!refitTrack_)
01564 extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
01565 else{
01566 GlobalTrackerMuonAlignment::trackFitter(itMuon->track(),
01567 trackTT,alongMomentum,trackFittedTSOS);
01568 if(trackFittedTSOS.isValid())
01569 extrapolationT = alongSmPr.propagate(trackFittedTSOS, refSurface);
01570 }
01571
01572 if(!extrapolationT.isValid()){
01573 if(false & alarm)
01574 std::cout<<" !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
01575
01576 <<std::endl;
01577 continue;
01578 }
01579 if(debug_)
01580 std::cout<<" extrapolationT.isValid "<<extrapolationT.isValid()<<std::endl;
01581
01582 tsosMuonIf = 1;
01583 Rt = GlobalVector((extrapolationT.globalPosition()).x(),
01584 (extrapolationT.globalPosition()).y(),
01585 (extrapolationT.globalPosition()).z());
01586
01587 Pt = extrapolationT.globalMomentum();
01588
01589 GRm = GlobalVector((innerMuTSOS.globalPosition()).x(),
01590 (innerMuTSOS.globalPosition()).y(),
01591 (innerMuTSOS.globalPosition()).z());
01592 GPm = innerMuTSOS.globalMomentum();
01593 Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
01594 (outerTrackTSOS.globalPosition()).y(),
01595 (outerTrackTSOS.globalPosition()).z());
01596 Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() +
01597 innerMuTSOS.cartesianError().matrix());
01598 C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
01599 Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
01600 C1 = AlgebraicSymMatrix66(innerMuTSOS.cartesianError().matrix());
01601
01602 if(refitMuon_)
01603 GlobalTrackerMuonAlignment::muonFitter(itMuon->outerTrack(),muTT,oppositeToMomentum,muonFittedTSOS);
01604
01605 if(false & debug_){
01606
01607 PropagationDirectionChooser Chooser;
01608 std::cout<<" propDirCh "
01609 <<Chooser.operator()(*outerTrackTSOS.freeState(), refSurface)
01610 <<" Ch == along "<<(alongMomentum ==
01611 Chooser.operator()(*outerTrackTSOS.freeState(), refSurface))
01612 <<std::endl;
01613 std::cout<<" --- Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
01614 <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
01615 }
01616 }
01617
01618
01619 else if((distanceInOut <= distanceInIn) & (distanceInOut <= distanceOutIn) &
01620 (distanceInOut <= distanceOutOut)){
01621 if(debug_) std::cout<<" ----- In - Out -----"<<std::endl;
01622
01623 const Surface& refSurface = outerMuTSOS.surface();
01624 ReferenceCountingPointer<TangentPlane>
01625 tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
01626 Nl = tpMuGlobal->normalVector();
01627
01628
01629
01630 if(!refitTrack_)
01631 extrapolationT = oppositeSmPr.propagate(innerTrackTSOS, refSurface);
01632 else{
01633 GlobalTrackerMuonAlignment::trackFitter(itMuon->track(),
01634 trackTT,oppositeToMomentum,trackFittedTSOS);
01635 if(trackFittedTSOS.isValid())
01636 extrapolationT = oppositeSmPr.propagate(trackFittedTSOS, refSurface);
01637 }
01638
01639 if(!extrapolationT.isValid()){
01640 if(false & alarm)
01641 std::cout<<" !!!!!!!!! Catastrophe Out-In extrapolationT.isValid "
01642 <<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
01643 <<std::endl;
01644 continue;
01645 }
01646 if(debug_)
01647 std::cout<<" extrapolationT.isValid "<<extrapolationT.isValid()<<std::endl;
01648
01649 tsosMuonIf = 2;
01650 Rt = GlobalVector((extrapolationT.globalPosition()).x(),
01651 (extrapolationT.globalPosition()).y(),
01652 (extrapolationT.globalPosition()).z());
01653
01654 Pt = extrapolationT.globalMomentum();
01655
01656 GRm = GlobalVector((outerMuTSOS.globalPosition()).x(),
01657 (outerMuTSOS.globalPosition()).y(),
01658 (outerMuTSOS.globalPosition()).z());
01659 GPm = outerMuTSOS.globalMomentum();
01660 Rt0 = GlobalVector((innerTrackTSOS.globalPosition()).x(),
01661 (innerTrackTSOS.globalPosition()).y(),
01662 (innerTrackTSOS.globalPosition()).z());
01663 Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() +
01664 outerMuTSOS.cartesianError().matrix());
01665 C0 = AlgebraicSymMatrix66(innerTrackTSOS.cartesianError().matrix());
01666 Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
01667 C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
01668
01669 if(refitMuon_)
01670 GlobalTrackerMuonAlignment::muonFitter(itMuon->outerTrack(),muTT,alongMomentum,muonFittedTSOS);
01671
01672 if(false & debug_){
01673
01674 PropagationDirectionChooser Chooser;
01675 std::cout<<" propDirCh "
01676 <<Chooser.operator()(*innerTrackTSOS.freeState(), refSurface)
01677 <<" Ch == oppisite "<<(oppositeToMomentum ==
01678 Chooser.operator()(*innerTrackTSOS.freeState(), refSurface))
01679 <<std::endl;
01680 std::cout<<" --- Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
01681 <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
01682 }
01683 }
01684
01685 else if((distanceOutOut <= distanceInOut) & (distanceOutOut <= distanceInIn) &
01686 (distanceOutOut <= distanceOutIn)){
01687 if(debug_) std::cout<<" ----- Out - Out -----"<<std::endl;
01688
01689
01690 continue;
01691
01692 const Surface& refSurface = outerMuTSOS.surface();
01693 ReferenceCountingPointer<TangentPlane>
01694 tpMuGlobal(refSurface.tangentPlane(outerMuTSOS.globalPosition()));
01695 Nl = tpMuGlobal->normalVector();
01696
01697
01698 extrapolationT = alongSmPr.propagate(outerTrackTSOS, refSurface);
01699
01700 if(!extrapolationT.isValid()){
01701 if(alarm)
01702 std::cout<<" !!!!!!!!! Catastrophe Out-Out extrapolationT.isValid "
01703 <<"\a\a\a\a\a\a\a\a"<<extrapolationT.isValid()
01704 <<std::endl;
01705 continue;
01706 }
01707 if(debug_)
01708 std::cout<<" extrapolationT.isValid "<<extrapolationT.isValid()<<std::endl;
01709
01710 tsosMuonIf = 2;
01711 Rt = GlobalVector((extrapolationT.globalPosition()).x(),
01712 (extrapolationT.globalPosition()).y(),
01713 (extrapolationT.globalPosition()).z());
01714
01715 Pt = extrapolationT.globalMomentum();
01716
01717 GRm = GlobalVector((outerMuTSOS.globalPosition()).x(),
01718 (outerMuTSOS.globalPosition()).y(),
01719 (outerMuTSOS.globalPosition()).z());
01720 GPm = outerMuTSOS.globalMomentum();
01721 Rt0 = GlobalVector((outerTrackTSOS.globalPosition()).x(),
01722 (outerTrackTSOS.globalPosition()).y(),
01723 (outerTrackTSOS.globalPosition()).z());
01724 Cm = AlgebraicSymMatrix66 (extrapolationT.cartesianError().matrix() +
01725 outerMuTSOS.cartesianError().matrix());
01726 C0 = AlgebraicSymMatrix66(outerTrackTSOS.cartesianError().matrix());
01727 Ce = AlgebraicSymMatrix66(extrapolationT.cartesianError().matrix());
01728 C1 = AlgebraicSymMatrix66(outerMuTSOS.cartesianError().matrix());
01729
01730 if(debug_){
01731 PropagationDirectionChooser Chooser;
01732 std::cout<<" propDirCh "
01733 <<Chooser.operator()(*outerTrackTSOS.freeState(), refSurface)
01734 <<" Ch == along "<<(alongMomentum ==
01735 Chooser.operator()(*outerTrackTSOS.freeState(), refSurface))
01736 <<std::endl;
01737 std::cout<<" --- Nlocal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
01738 <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
01739 std::cout<<" Nornal "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "
01740 <<" alfa "<<int(asin(Nl.x())*57.29578)<<std::endl;
01741 }
01742 }
01743 else{
01744 if(alarm)
01745 std::cout<<" ------- !!!!!!!!!!!!!!! No proper Out - In \a\a\a\a\a\a\a"<<std::endl;
01746 continue;
01747 }
01748
01749 }
01750
01751 if(tsosMuonIf == 0) {if(info) {std::cout<<"No tsosMuon !!!!!!"<<std::endl; continue;}}
01752 TrajectoryStateOnSurface tsosMuon;
01753 if(tsosMuonIf == 1) tsosMuon = muTT.innermostMeasurementState();
01754 else tsosMuon = muTT.outermostMeasurementState();
01755
01756
01757 AlgebraicVector4 LPRm;
01758 GlobalTrackerMuonAlignment::misalignMuonL
01759 (GRm, GPm, Nl, Rt, Rm, Pm, LPRm, extrapolationT, tsosMuon);
01760
01761 if(refitTrack_){
01762 if(!trackFittedTSOS.isValid()){
01763 if(info) std::cout<<" ================= trackFittedTSOS notValid !!!!!!!! "<<std::endl;
01764 continue;}
01765 if(debug_) this->debugTrajectorySOS(" trackFittedTSOS ", trackFittedTSOS);
01766 }
01767
01768 if(refitMuon_){
01769 if(!muonFittedTSOS.isValid()){
01770 if(info) std::cout<<" ================= muonFittedTSOS notValid !!!!!!!! "<<std::endl;
01771 continue;}
01772 if(debug_) this->debugTrajectorySOS(" muonFittedTSOS ", muonFittedTSOS);
01773 Rm = GlobalVector((muonFittedTSOS.globalPosition()).x(),
01774 (muonFittedTSOS.globalPosition()).y(),
01775 (muonFittedTSOS.globalPosition()).z());
01776 Pm = muonFittedTSOS.globalMomentum();
01777 LPRm = AlgebraicVector4(muonFittedTSOS.localParameters().vector()(1),
01778 muonFittedTSOS.localParameters().vector()(2),
01779 muonFittedTSOS.localParameters().vector()(3),
01780 muonFittedTSOS.localParameters().vector()(4));
01781 }
01782 GlobalVector resR = Rm - Rt;
01783 GlobalVector resP0 = Pm - Pt;
01784 GlobalVector resP = Pm / Pm.mag() - Pt / Pt.mag();
01785 float RelMomResidual = (Pm.mag() - Pt.mag()) / (Pt.mag() + 1.e-6);;
01786
01787 AlgebraicVector6 Vm;
01788 Vm(0) = resR.x(); Vm(1) = resR.y(); Vm(2) = resR.z();
01789 Vm(3) = resP0.x(); Vm(4) = resP0.y(); Vm(5) = resP0.z();
01790 float Rmuon = Rm.perp();
01791 float Zmuon = Rm.z();
01792 float alfa_x = atan2(Nl.x(),Nl.y())*57.29578;
01793
01794 if(debug_){
01795 std::cout<<" Nx Ny Nz alfa_x "<<Nl.x()<<" "<<Nl.y()<<" "<<Nl.z()<<" "<<alfa_x<<std::endl;
01796
01797
01798 }
01799
01800 double chi_d = 0;
01801 for(int i=0; i<=5; i++) chi_d += Vm(i)*Vm(i)/Cm(i,i);
01802
01803 AlgebraicVector5 Vml(tsosMuon.localParameters().vector()
01804 - extrapolationT.localParameters().vector());
01805 AlgebraicSymMatrix55 m(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
01806 AlgebraicSymMatrix55 Cml(tsosMuon.localError().matrix() + extrapolationT.localError().matrix());
01807 bool ierrLoc = !m.Invert();
01808 if (ierrLoc && debug_ && info) {
01809 std::cout<< " ==== Error inverting Local covariance matrix ==== "<<std::cout;
01810 continue;}
01811 double chi_Loc = ROOT::Math::Similarity(Vml,m);
01812 if(debug_)
01813 std::cout<<" chi_Loc px/pz/err "<<chi_Loc<<" "<<Vml(1)/sqrt(Cml(1,1))<<std::endl;
01814
01815 if(Pt.mag() < 15.) continue;
01816 if(Pm.mag() < 5.) continue;
01817
01818
01819
01820
01821
01822
01823
01824
01825
01826
01827
01828
01829
01830
01831 if(fabs(resR.x()) > 20.) continue;
01832 if(fabs(resR.y()) > 20.) continue;
01833 if(fabs(resR.z()) > 20.) continue;
01834 if(fabs(resR.mag()) > 30.) continue;
01835 if(fabs(resP.x()) > 0.06) continue;
01836 if(fabs(resP.y()) > 0.06) continue;
01837 if(fabs(resP.z()) > 0.06) continue;
01838 if(chi_d > 40.) continue;
01839
01840
01841
01842
01843
01844
01845
01846
01847
01848
01849
01850
01851
01852
01853
01854
01855
01856
01857
01858 if(Rmuon < 120. || Rmuon > 450.) continue;
01859 if(Zmuon < -720. || Zmuon > 720.) continue;
01860 MuSelect = " Barrel+EndCaps";
01861
01862 if(debug_)
01863 std::cout<<" .............. passed all cuts"<<std::endl;
01864
01865 N_track++;
01866
01867
01868 GlobalTrackerMuonAlignment::gradientGlobalAlg(Rt, Pt, Rm, Nl, Cm);
01869 GlobalTrackerMuonAlignment::gradientGlobal(Rt, Pt, Rm, Pm, Nl, Cm);
01870
01871 CLHEP::HepSymMatrix covLoc(4,0);
01872 for(int i=1; i<=4; i++)
01873 for(int j=1; j<=i; j++){
01874 covLoc(i,j) = (tsosMuon.localError().matrix()
01875 + extrapolationT.localError().matrix())(i,j);
01876
01877 }
01878
01879 const Surface& refSurface = tsosMuon.surface();
01880 CLHEP::HepMatrix rotLoc (3,3,0);
01881 rotLoc(1,1) = refSurface.rotation().xx();
01882 rotLoc(1,2) = refSurface.rotation().xy();
01883 rotLoc(1,3) = refSurface.rotation().xz();
01884
01885 rotLoc(2,1) = refSurface.rotation().yx();
01886 rotLoc(2,2) = refSurface.rotation().yy();
01887 rotLoc(2,3) = refSurface.rotation().yz();
01888
01889 rotLoc(3,1) = refSurface.rotation().zx();
01890 rotLoc(3,2) = refSurface.rotation().zy();
01891 rotLoc(3,3) = refSurface.rotation().zz();
01892
01893 CLHEP::HepVector posLoc(3);
01894 posLoc(1) = refSurface.position().x();
01895 posLoc(2) = refSurface.position().y();
01896 posLoc(3) = refSurface.position().z();
01897
01898 GlobalTrackerMuonAlignment::gradientLocal(Rt, Pt, Rm, Pm, Nl, covLoc, rotLoc, posLoc, LPRm);
01899
01900 if(debug_){
01901 std::cout<<" Norm "<<Nl<<std::endl;
01902 std::cout<<" posLoc "<<posLoc.T()<<std::endl;
01903 std::cout<<" rotLoc "<<rotLoc<<std::endl;
01904 }
01905
01906
01907 histo->Fill(itMuon->track()->pt());
01908
01909
01910 histo2->Fill(Pt.mag());
01911 histo3->Fill((PI/2.-itMuon->track()->outerTheta()));
01912 histo4->Fill(itMuon->track()->phi());
01913 histo5->Fill(Rmuon);
01914 histo6->Fill(Zmuon);
01915 histo7->Fill(RelMomResidual);
01916
01917 histo8->Fill(chi_d);
01918
01919 histo101->Fill(Zmuon, Rmuon);
01920 histo101->Fill(Rt0.z(), Rt0.perp());
01921 histo102->Fill(Rt0.x(), Rt0.y());
01922 histo102->Fill(Rm.x(), Rm.y());
01923
01924 histo11->Fill(resR.mag());
01925 if(fabs(Nl.x()) < 0.98) histo12->Fill(resR.x());
01926 if(fabs(Nl.y()) < 0.98) histo13->Fill(resR.y());
01927 if(fabs(Nl.z()) < 0.98) histo14->Fill(resR.z());
01928 histo15->Fill(resP.x());
01929 histo16->Fill(resP.y());
01930 histo17->Fill(resP.z());
01931
01932 if((fabs(Nl.x()) < 0.98) && (fabs(Nl.y()) < 0.98) &&(fabs(Nl.z()) < 0.98))
01933 {
01934 histo18->Fill(sqrt(C0(0,0)));
01935 histo19->Fill(sqrt(C1(0,0)));
01936 histo20->Fill(sqrt(C1(0,0)+Ce(0,0)));
01937 }
01938 if(fabs(Nl.x()) < 0.98) histo21->Fill(Vm(0)/sqrt(Cm(0,0)));
01939 if(fabs(Nl.y()) < 0.98) histo22->Fill(Vm(1)/sqrt(Cm(1,1)));
01940 if(fabs(Nl.z()) < 0.98) histo23->Fill(Vm(2)/sqrt(Cm(2,2)));
01941 histo24->Fill(Vm(3)/sqrt(C1(3,3)+Ce(3,3)));
01942 histo25->Fill(Vm(4)/sqrt(C1(4,4)+Ce(4,4)));
01943 histo26->Fill(Vm(5)/sqrt(C1(5,5)+Ce(5,5)));
01944 histo27->Fill(Nl.x());
01945 histo28->Fill(Nl.y());
01946 histo29->Fill(lenghtTrack);
01947 histo30->Fill(lenghtMuon);
01948 histo31->Fill(chi_Loc);
01949 histo32->Fill(Vml(1)/sqrt(Cml(1,1)));
01950 histo33->Fill(Vml(2)/sqrt(Cml(2,2)));
01951 histo34->Fill(Vml(3)/sqrt(Cml(3,3)));
01952 histo35->Fill(Vml(4)/sqrt(Cml(4,4)));
01953
01954 if (debug_) {
01955
01956 std::cout<<" diag 0[ "<<C0(0,0)<<" "<<C0(1,1)<<" "<<C0(2,2)<<" "<<C0(3,3)<<" "
01957 <<C0(4,4)<<" "<<C0(5,5)<<" ]"<<std::endl;
01958 std::cout<<" diag e[ "<<Ce(0,0)<<" "<<Ce(1,1)<<" "<<Ce(2,2)<<" "<<Ce(3,3)<<" "
01959 <<Ce(4,4)<<" "<<Ce(5,5)<<" ]"<<std::endl;
01960 std::cout<<" diag 1[ "<<C1(0,0)<<" "<<C1(1,1)<<" "<<C1(2,2)<<" "<<C1(3,3)<<" "
01961 <<C1(4,4)<<" "<<C1(5,5)<<" ]"<<std::endl;
01962 std::cout<<" Rm "<<Rm.x()<<" "<<Rm.y()<<" "<<Rm.z()
01963 <<" Pm "<<Pm.x()<<" "<<Pm.y()<<" "<<Pm.z()<<std::endl;
01964 std::cout<<" Rt "<<Rt.x()<<" "<<Rt.y()<<" "<<Rt.z()
01965 <<" Pt "<<Pt.x()<<" "<<Pt.y()<<" "<<Pt.z()<<std::endl;
01966 std::cout<<" Nl*(Rm-Rt) "<<Nl.dot(Rm-Rt)<<std::endl;
01967 std::cout<<" resR "<<resR.x()<<" "<<resR.y()<<" "<<resR.z()
01968 <<" resP "<<resP.x()<<" "<<resP.y()<<" "<<resP.z()<<std::endl;
01969 std::cout<<" Rm-t "<<(Rm-Rt).x()<<" "<<(Rm-Rt).y()<<" "<<(Rm-Rt).z()
01970 <<" Pm-t "<<(Pm-Pt).x()<<" "<<(Pm-Pt).y()<<" "<<(Pm-Pt).z()<<std::endl;
01971 std::cout<<" Vm "<<Vm<<std::endl;
01972 std::cout<<" +- "<<sqrt(Cm(0,0))<<" "<<sqrt(Cm(1,1))<<" "<<sqrt(Cm(2,2))<<" "
01973 <<sqrt(Cm(3,3))<<" "<<sqrt(Cm(4,4))<<" "<<sqrt(Cm(5,5))<<std::endl;
01974 std::cout<<" Pmuon Ptrack dP/Ptrack "<<itMuon->outerTrack()->p()<<" "
01975 <<itMuon->track()->outerP()<<" "
01976 <<(itMuon->outerTrack()->p() -
01977 itMuon->track()->outerP())/itMuon->track()->outerP()<<std::endl;
01978 std::cout<<" cov matrix "<<std::endl;
01979 std::cout<<Cm<<std::endl;
01980 std::cout<<" diag [ "<<Cm(0,0)<<" "<<Cm(1,1)<<" "<<Cm(2,2)<<" "<<Cm(3,3)<<" "
01981 <<Cm(4,4)<<" "<<Cm(5,5)<<" ]"<<std::endl;
01982
01983 static AlgebraicSymMatrix66 Ro;
01984 double Diag[6];
01985 for(int i=0; i<=5; i++) Diag[i] = sqrt(Cm(i,i));
01986 for(int i=0; i<=5; i++)
01987 for(int j=0; j<=5; j++)
01988 Ro(i,j) = Cm(i,j)/Diag[i]/Diag[j];
01989 std::cout<<" correlation matrix "<<std::endl;
01990 std::cout<<Ro<<std::endl;
01991
01992 AlgebraicSymMatrix66 CmI;
01993 for(int i=0; i<=5; i++)
01994 for(int j=0; j<=5; j++)
01995 CmI(i,j) = Cm(i,j);
01996
01997 bool ierr = !CmI.Invert();
01998 if( ierr ) { if(alarm || debug_)
01999 std::cout<<" Error inverse covariance matrix !!!!!!!!!!!"<<std::endl;
02000 continue;
02001 }
02002 std::cout<<" inverse cov matrix "<<std::endl;
02003 std::cout<<Cm<<std::endl;
02004
02005 double chi = ROOT::Math::Similarity(Vm, CmI);
02006 std::cout<<" chi chi_d "<<chi<<" "<<chi_d<<std::endl;
02007 }
02008
02009 }
02010
02011 }
02012
02013
02014
02015 void
02016 GlobalTrackerMuonAlignment::gradientGlobalAlg(GlobalVector& Rt, GlobalVector& Pt,
02017 GlobalVector& Rm, GlobalVector& Nl,
02018 AlgebraicSymMatrix66& Cm)
02019 {
02020
02021
02022
02023 AlgebraicMatrix33 Jac;
02024 AlgebraicVector3 Wi, R_m, R_t, P_t, Norm, dR;
02025
02026 R_m(0) = Rm.x(); R_m(1) = Rm.y(); R_m(2) = Rm.z();
02027 R_t(0) = Rt.x(); R_t(1) = Rt.y(); R_t(2) = Rt.z();
02028 P_t(0) = Pt.x(); P_t(1) = Pt.y(); P_t(2) = Pt.z();
02029 Norm(0) = Nl.x(); Norm(1) = Nl.y(); Norm(2) = Nl.z();
02030
02031 for(int i=0; i<=2; i++){
02032 if(Cm(i,i) > 1.e-20)
02033 Wi(i) = 1./Cm(i,i);
02034 else Wi(i) = 1.e-10;
02035 dR(i) = R_m(i) - R_t(i);
02036 }
02037
02038 float PtN = P_t(0)*Norm(0) + P_t(1)*Norm(1) + P_t(2)*Norm(2);
02039
02040 Jac(0,0) = 1. - P_t(0)*Norm(0)/PtN;
02041 Jac(0,1) = - P_t(0)*Norm(1)/PtN;
02042 Jac(0,2) = - P_t(0)*Norm(2)/PtN;
02043
02044 Jac(1,0) = - P_t(1)*Norm(0)/PtN;
02045 Jac(1,1) = 1. - P_t(1)*Norm(1)/PtN;
02046 Jac(1,2) = - P_t(1)*Norm(2)/PtN;
02047
02048 Jac(2,0) = - P_t(2)*Norm(0)/PtN;
02049 Jac(2,1) = - P_t(2)*Norm(1)/PtN;
02050 Jac(2,2) = 1. - P_t(2)*Norm(2)/PtN;
02051
02052 AlgebraicSymMatrix33 Itr;
02053
02054 for(int i=0; i<=2; i++)
02055 for(int j=0; j<=2; j++){
02056 if(j < i) continue;
02057 Itr(i,j) = 0.;
02058
02059 for(int k=0; k<=2; k++){
02060 Itr(i,j) += Jac(k,i)*Wi(k)*Jac(k,j);}}
02061
02062 for(int i=0; i<=2; i++)
02063 for(int j=0; j<=2; j++){
02064 if(j < i) continue;
02065 Inf(i,j) += Itr(i,j);}
02066
02067 AlgebraicVector3 Gtr(0., 0., 0.);
02068 for(int i=0; i<=2; i++)
02069 for(int k=0; k<=2; k++) Gtr(i) += dR(k)*Wi(k)*Jac(k,i);
02070 for(int i=0; i<=2; i++) Gfr(i) += Gtr(i);
02071
02072 if(debug_){
02073 std::cout<<" Wi "<<Wi<<std::endl;
02074 std::cout<<" N "<<Norm<<std::endl;
02075 std::cout<<" P_t "<<P_t<<std::endl;
02076 std::cout<<" (Pt*N) "<<PtN<<std::endl;
02077 std::cout<<" dR "<<dR<<std::endl;
02078 std::cout<<" +/- "<<1./sqrt(Wi(0))<<" "<<1./sqrt(Wi(1))<<" "<<1./sqrt(Wi(2))
02079 <<" "<<std::endl;
02080 std::cout<<" Jacobian dr/ddx "<<std::endl;
02081 std::cout<<Jac<<std::endl;
02082 std::cout<<" G-- "<<Gtr<<std::endl;
02083 std::cout<<" Itrack "<<std::endl;
02084 std::cout<<Itr<<std::endl;
02085 std::cout<<" Gfr "<<Gfr<<std::endl;
02086 std::cout<<" -- Inf --"<<std::endl;
02087 std::cout<<Inf<<std::endl;
02088 }
02089
02090 return;
02091 }
02092
02093
02094 void
02095 GlobalTrackerMuonAlignment::gradientGlobal(GlobalVector& GRt, GlobalVector& GPt,
02096 GlobalVector& GRm, GlobalVector& GPm,
02097 GlobalVector& GNorm, AlgebraicSymMatrix66 & GCov)
02098 {
02099
02100
02101
02102
02103
02104 bool alarm = false;
02105 bool info = false;
02106
02107 int Nd = 6;
02108
02109
02110 CLHEP::HepSymMatrix w(Nd,0);
02111 for (int i=1; i <= Nd; i++)
02112 for (int j=1; j <= Nd; j++){
02113 if(j <= i ) w(i,j) = GCov(i-1, j-1);
02114
02115
02116 if((i == j) && (i<=3) && (GCov(i-1, j-1) < 1.e-20)) w(i,j) = 1.e20;
02117 if(i != j) w(i,j) = 0.;
02118 }
02119
02120
02121
02122
02123 CLHEP::HepVector V(Nd), Rt(3), Pt(3), Rm(3), Pm(3), Norm(3);
02124 Rt(1) = GRt.x(); Rt(2) = GRt.y(); Rt(3) = GRt.z();
02125 Pt(1) = GPt.x(); Pt(2) = GPt.y(); Pt(3) = GPt.z();
02126 Rm(1) = GRm.x(); Rm(2) = GRm.y(); Rm(3) = GRm.z();
02127 Pm(1) = GPm.x(); Pm(2) = GPm.y(); Pm(3) = GPm.z();
02128 Norm(1) = GNorm.x(); Norm(2) = GNorm.y(); Norm(3) = GNorm.z();
02129
02130 V = dsum(Rm - Rt, Pm - Pt) ;
02131
02132
02133 double PmN = CLHEP_dot(Pm, Norm);
02134
02135 CLHEP::HepMatrix Jac(Nd,Nd,0);
02136 for (int i=1; i <= 3; i++)
02137 for (int j=1; j <= 3; j++){
02138 Jac(i,j) = Pm(i)*Norm(j) / PmN;
02139 if(i == j ) Jac(i,j) -= 1.;
02140 }
02141
02142
02143 Jac(4,4) = 0.;
02144 Jac(5,4) = -Pm(3);
02145 Jac(6,4) = Pm(2);
02146 Jac(4,5) = Pm(3);
02147 Jac(5,5) = 0.;
02148 Jac(6,5) = -Pm(1);
02149 Jac(4,6) = -Pm(2);
02150 Jac(5,6) = Pm(1);
02151 Jac(6,6) = 0.;
02152
02153 CLHEP::HepVector dsda(3);
02154 dsda(1) = (Norm(2)*Rm(3) - Norm(3)*Rm(2)) / PmN;
02155 dsda(2) = (Norm(3)*Rm(1) - Norm(1)*Rm(3)) / PmN;
02156 dsda(3) = (Norm(1)*Rm(2) - Norm(2)*Rm(1)) / PmN;
02157
02158
02159 Jac(1,4) = Pm(1)*dsda(1);
02160 Jac(2,4) = -Rm(3) + Pm(2)*dsda(1);
02161 Jac(3,4) = Rm(2) + Pm(3)*dsda(1);
02162
02163 Jac(1,5) = Rm(3) + Pm(1)*dsda(2);
02164 Jac(2,5) = Pm(2)*dsda(2);
02165 Jac(3,5) = -Rm(1) + Pm(3)*dsda(2);
02166
02167 Jac(1,6) = -Rm(2) + Pm(1)*dsda(3);
02168 Jac(2,6) = Rm(1) + Pm(2)*dsda(3);
02169 Jac(3,6) = Pm(3)*dsda(3);
02170
02171 CLHEP::HepSymMatrix W(Nd,0);
02172 int ierr;
02173 W = w.inverse(ierr);
02174 if(ierr != 0) { if(alarm)
02175 std::cout<<" gradientGlobal: inversion of matrix w fail "<<std::endl;
02176 return;
02177 }
02178
02179 CLHEP::HepMatrix W_Jac(Nd,Nd,0);
02180 W_Jac = Jac.T() * W;
02181
02182 CLHEP::HepVector grad3(Nd);
02183 grad3 = W_Jac * V;
02184
02185 CLHEP::HepMatrix hess3(Nd,Nd);
02186 hess3 = Jac.T() * W * Jac;
02187
02188
02189 Grad += grad3;
02190 Hess += hess3;
02191
02192 CLHEP::HepVector d3I = CLHEP::solve(Hess, -Grad);
02193 int iEr3I;
02194 CLHEP::HepMatrix Errd3I = Hess.inverse(iEr3I);
02195 if( iEr3I != 0) { if(alarm)
02196 std::cout<<" gradientGlobal error inverse Hess matrix !!!!!!!!!!!"<<std::endl;
02197 }
02198
02199 if(info || debug_){
02200 std::cout<<" dG "<<d3I(1)<<" "<<d3I(2)<<" "<<d3I(3)<<" "<<d3I(4)<<" "
02201 <<d3I(5)<<" "<<d3I(6)<<std::endl;;
02202 std::cout<<" +- "<<sqrt(Errd3I(1,1))<<" "<<sqrt(Errd3I(2,2))<<" "<<sqrt(Errd3I(3,3))
02203 <<" "<<sqrt(Errd3I(4,4))<<" "<<sqrt(Errd3I(5,5))<<" "<<sqrt(Errd3I(6,6))
02204 <<std::endl;
02205 }
02206
02207 #ifdef CHECK_OF_DERIVATIVES
02208
02209
02210 CLHEP::HepVector d(3,0), a(3,0);
02211 CLHEP::HepMatrix T(3,3,1);
02212
02213 CLHEP::HepMatrix Ti = T.T();
02214
02215
02216 double A = CLHEP_dot(Ti*Pm, Norm);
02217 double B = CLHEP_dot((Rt -Ti*Rm + Ti*d), Norm);
02218 double s0 = B / A;
02219
02220 CLHEP::HepVector r0(3,0), p0(3,0);
02221 r0 = Ti*Rm - Ti*d + s0*(Ti*Pm) - Rt;
02222 p0 = Ti*Pm - Pt;
02223
02224 double delta = 0.0001;
02225
02226 int ii = 3; d(ii) += delta;
02227
02228
02229
02230 Ti = T.T();
02231
02232
02233 A = CLHEP_dot(Ti*Pm, Norm);
02234 B = CLHEP_dot((Rt -Ti*Rm + Ti*d), Norm);
02235 double s = B / A;
02236
02237 CLHEP::HepVector r(3,0), p(3,0);
02238 r = Ti*Rm - Ti*d + s*(Ti*Pm) - Rt;
02239 p = Ti*Pm - Pt;
02240
02241 std::cout<<" s0 s num dsda("<<ii<<") "<<s0<<" "<<s<<" "
02242 <<(s-s0)/delta<<" "<<dsda(ii)<<endl;
02243
02244 std::cout<<" -- An d(r,p)/d("<<ii<<") "<<Jac(1,ii)<<" "<<Jac(2,ii)<<" "
02245 <<Jac(3,ii)<<" "<<Jac(4,ii)<<" "<<Jac(5,ii)<<" "
02246 <<Jac(6,ii)<<std::endl;
02247 std::cout<<" Nu d(r,p)/d("<<ii<<") "<<(r(1)-r0(1))/delta<<" "
02248 <<(r(2)-r0(2))/delta<<" "<<(r(3)-r0(3))/delta<<" "
02249 <<(p(1)-p0(1))/delta<<" "<<(p(2)-p0(2))/delta<<" "
02250 <<(p(3)-p0(3))/delta<<std::endl;
02251
02252 std::cout<<" -- An d(r,p)/a("<<ii<<") "<<Jac(1,ii+3)<<" "<<Jac(2,ii+3)<<" "
02253 <<Jac(3,ii+3)<<" "<<Jac(4,ii+3)<<" "<<Jac(5,ii+3)<<" "
02254 <<Jac(6,ii+3)<<std::endl;
02255 std::cout<<" Nu d(r,p)/a("<<ii<<") "<<(r(1)-r0(1))/delta<<" "
02256 <<(r(2)-r0(2))/delta<<" "<<(r(3)-r0(3))/delta<<" "
02257 <<(p(1)-p0(1))/delta<<" "<<(p(2)-p0(2))/delta<<" "
02258 <<(p(3)-p0(3))/delta<<std::endl;
02259
02260 #endif
02261
02262 return;
02263 }
02264
02265
02266
02267 void
02268 GlobalTrackerMuonAlignment::gradientLocal(GlobalVector& GRt, GlobalVector& GPt,
02269 GlobalVector& GRm, GlobalVector& GPm,
02270 GlobalVector& GNorm,
02271 CLHEP::HepSymMatrix& covLoc,
02272 CLHEP::HepMatrix& rotLoc,
02273 CLHEP::HepVector& R0,
02274 AlgebraicVector4& LPRm)
02275 {
02276
02277
02278
02279
02280 bool alarm = true;
02281
02282 bool info = false;
02283
02284 if(debug_)
02285 std::cout<<" gradientLocal "<<std::endl;
02286
02287
02288
02289
02290
02291
02292
02293
02294
02295
02296
02297
02298
02299
02300
02301
02302
02303
02304 CLHEP::HepVector Rt(3), Pt(3), Rm(3), Pm(3), Norm(3);
02305 Rt(1) = GRt.x(); Rt(2) = GRt.y(); Rt(3) = GRt.z();
02306 Pt(1) = GPt.x(); Pt(2) = GPt.y(); Pt(3) = GPt.z();
02307 Rm(1) = GRm.x(); Rm(2) = GRm.y(); Rm(3) = GRm.z();
02308 Pm(1) = GPm.x(); Pm(2) = GPm.y(); Pm(3) = GPm.z();
02309 Norm(1) = GNorm.x(); Norm(2) = GNorm.y(); Norm(3) = GNorm.z();
02310
02311 CLHEP::HepVector V(4), Rml(3), Pml(3), Rtl(3), Ptl(3);
02312
02313
02314
02315
02316
02317
02318
02319 Rml = rotLoc * (Rm - R0);
02320 Rtl = rotLoc * (Rt - R0);
02321 Pml = rotLoc * Pm;
02322 Ptl = rotLoc * Pt;
02323
02324 V(1) = LPRm(0) - Ptl(1)/Ptl(3);
02325 V(2) = LPRm(1) - Ptl(2)/Ptl(3);
02326 V(3) = LPRm(2) - Rtl(1);
02327 V(4) = LPRm(3) - Rtl(2);
02328
02329
02330
02331
02332
02333
02334
02335
02336
02337
02338
02339
02340
02341
02342 CLHEP::HepSymMatrix W = covLoc;
02343
02344 int ierr;
02345 W.invert(ierr);
02346 if(ierr != 0) { if(alarm)
02347 std::cout<<" gradientLocal: inversion of matrix W fail "<<std::endl;
02348 return;
02349 }
02350
02351
02352
02353
02354
02355
02356
02357
02358
02359 CLHEP::HepMatrix JacToLoc(4,6,0);
02360 for(int i=1; i<=2; i++)
02361 for(int j=1; j<=3; j++){
02362 JacToLoc(i,j+3) = (rotLoc(i,j) - rotLoc(3,j)*Pml(i)/Pml(3))/Pml(3);
02363 JacToLoc(i+2,j) = rotLoc(i,j);
02364 }
02365
02366
02367
02368 double PmN = CLHEP_dot(Pm, Norm);
02369
02370 CLHEP::HepMatrix Jac(6,6,0);
02371 for (int i=1; i <= 3; i++)
02372 for (int j=1; j <= 3; j++){
02373 Jac(i,j) = Pm(i)*Norm(j) / PmN;
02374 if(i == j ) Jac(i,j) -= 1.;
02375 }
02376
02377
02378 Jac(4,4) = 0.;
02379 Jac(5,4) = -Pm(3);
02380 Jac(6,4) = Pm(2);
02381 Jac(4,5) = Pm(3);
02382 Jac(5,5) = 0.;
02383 Jac(6,5) = -Pm(1);
02384 Jac(4,6) = -Pm(2);
02385 Jac(5,6) = Pm(1);
02386 Jac(6,6) = 0.;
02387
02388 CLHEP::HepVector dsda(3);
02389 dsda(1) = (Norm(2)*Rm(3) - Norm(3)*Rm(2)) / PmN;
02390 dsda(2) = (Norm(3)*Rm(1) - Norm(1)*Rm(3)) / PmN;
02391 dsda(3) = (Norm(1)*Rm(2) - Norm(2)*Rm(1)) / PmN;
02392
02393
02394 Jac(1,4) = Pm(1)*dsda(1);
02395 Jac(2,4) = -Rm(3) + Pm(2)*dsda(1);
02396 Jac(3,4) = Rm(2) + Pm(3)*dsda(1);
02397
02398 Jac(1,5) = Rm(3) + Pm(1)*dsda(2);
02399 Jac(2,5) = Pm(2)*dsda(2);
02400 Jac(3,5) = -Rm(1) + Pm(3)*dsda(2);
02401
02402 Jac(1,6) = -Rm(2) + Pm(1)*dsda(3);
02403 Jac(2,6) = Rm(1) + Pm(2)*dsda(3);
02404 Jac(3,6) = Pm(3)*dsda(3);
02405
02406
02407 CLHEP::HepMatrix JacCorLoc(4,6,0);
02408 JacCorLoc = JacToLoc * Jac;
02409
02410
02411 CLHEP::HepMatrix W_Jac(6,4,0);
02412 W_Jac = JacCorLoc.T() * W;
02413
02414 CLHEP::HepVector gradL(6);
02415 gradL = W_Jac * V;
02416
02417 CLHEP::HepMatrix hessL(6,6);
02418 hessL = JacCorLoc.T() * W * JacCorLoc;
02419
02420 GradL += gradL;
02421 HessL += hessL;
02422
02423 CLHEP::HepVector dLI = CLHEP::solve(HessL, -GradL);
02424 int iErI;
02425 CLHEP::HepMatrix ErrdLI = HessL.inverse(iErI);
02426 if( iErI != 0) { if(alarm)
02427 std::cout<<" gradLocal Error inverse Hess matrix !!!!!!!!!!!"<<std::endl;
02428 }
02429
02430 if(info || debug_){
02431 std::cout<<" dL "<<dLI(1)<<" "<<dLI(2)<<" "<<dLI(3)<<" "<<dLI(4)<<" "
02432 <<dLI(5)<<" "<<dLI(6)<<std::endl;;
02433 std::cout<<" +- "<<sqrt(ErrdLI(1,1))<<" "<<sqrt(ErrdLI(2,2))<<" "<<sqrt(ErrdLI(3,3))
02434 <<" "<<sqrt(ErrdLI(4,4))<<" "<<sqrt(ErrdLI(5,5))<<" "<<sqrt(ErrdLI(6,6))
02435 <<std::endl;
02436 }
02437
02438 if(debug_){
02439
02440
02441
02442
02443
02444
02445
02446
02447 int i=5;
02448 if(GNorm.z() > 0.95)
02449 std::cout<<" Ecap1 N "<<GNorm<<std::endl;
02450 else if(GNorm.z() < -0.95)
02451 std::cout<<" Ecap2 N "<<GNorm<<std::endl;
02452 else
02453 std::cout<<" Barrel N "<<GNorm<<std::endl;
02454 std::cout<<" JacCLo(i,"<<i<<") = {"<<JacCorLoc(1,i)<<" "<<JacCorLoc(2,i)<<" "
02455 <<JacCorLoc(3,i)<<" "<<JacCorLoc(4,i)<<"}"<<std::endl;
02456 std::cout<<" rotLoc "<<rotLoc<<std::endl;
02457 std::cout<<" position "<<R0<<std::endl;
02458 std::cout<<" Pm,l "<<Pm.T()<<" "<<Pml(1)/Pml(3)<<" "<<Pml(2)/Pml(3)<<std::endl;
02459 std::cout<<" Pt,l "<<Pt.T()<<" "<<Ptl(1)/Ptl(3)<<" "<<Ptl(2)/Ptl(3)<<std::endl;
02460 std::cout<<" V "<<V.T()<<std::endl;
02461 std::cout<<" Cov \n"<<covLoc<<std::endl;
02462 std::cout<<" W*Cov "<<W*covLoc<<std::endl;
02463
02464
02465
02466 std::cout<<" JacToLoc "<<JacToLoc<<std::endl;
02467
02468 }
02469
02470 #ifdef CHECK_OF_JACOBIAN_CARTESIAN_TO_LOCAL
02471
02472 CLHEP::HepVector V0(4,0);
02473 V0(1) = Pml(1)/Pml(3) - Ptl(1)/Ptl(3);
02474 V0(2) = Pml(2)/Pml(3) - Ptl(2)/Ptl(3);
02475 V0(3) = Rml(1) - Rtl(1);
02476 V0(4) = Rml(2) - Rtl(2);
02477 int ii = 3; float delta = 0.01;
02478 CLHEP::HepVector V1(4,0);
02479 if(ii <= 3) {Rm(ii) += delta; Rml = rotLoc * (Rm - R0);}
02480 else {Pm(ii-3) += delta; Pml = rotLoc * Pm;}
02481
02482
02483 V1(1) = Pml(1)/Pml(3) - Ptl(1)/Ptl(3);
02484 V1(2) = Pml(2)/Pml(3) - Ptl(2)/Ptl(3);
02485 V1(3) = Rml(1) - Rtl(1);
02486 V1(4) = Rml(2) - Rtl(2);
02487
02488 if(GNorm.z() > 0.95)
02489 std::cout<<" Ecap1 N "<<GNorm<<std::endl;
02490 else if(GNorm.z() < -0.95)
02491 std::cout<<" Ecap2 N "<<GNorm<<std::endl;
02492 else
02493 std::cout<<" Barrel N "<<GNorm<<std::endl;
02494 std::cout<<" dldc Num (i,"<<ii<<") "<<(V1(1)-V0(1))/delta<<" "<<(V1(2)-V0(2))/delta<<" "
02495 <<(V1(3)-V0(3))/delta<<" "<<(V1(4)-V0(4))/delta<<std::endl;
02496 std::cout<<" dldc Ana (i,"<<ii<<") "<<JacToLoc(1,ii)<<" "<<JacToLoc(2,ii)<<" "
02497 <<JacToLoc(3,ii)<<" "<<JacToLoc(4,ii)<<std::endl;
02498
02499
02500
02501 #endif
02502
02503 return;
02504 }
02505
02506
02507
02508 void
02509 GlobalTrackerMuonAlignment::misalignMuon(GlobalVector& GRm, GlobalVector& GPm, GlobalVector& Nl,
02510 GlobalVector& Rt, GlobalVector& Rm, GlobalVector& Pm)
02511 {
02512 CLHEP::HepVector d(3,0), a(3,0), Norm(3), Rm0(3), Pm0(3), Rm1(3), Pm1(3),
02513 Rmi(3), Pmi(3) ;
02514
02515 d(1)=0.0; d(2)=0.0; d(3)=0.0; a(1)=0.0000; a(2)=0.0000; a(3)=0.0000;
02516
02517
02518
02519
02520
02521
02522
02523
02524
02525
02526
02527 MuGlShift = d; MuGlAngle = a;
02528
02529 if((d(1) == 0.) & (d(2) == 0.) & (d(3) == 0.) &
02530 (a(1) == 0.) & (a(2) == 0.) & (a(3) == 0.)){
02531 Rm = GRm;
02532 Pm = GPm;
02533 if(debug_)
02534 std::cout<<" ...... without misalignment "<<std::endl;
02535 return;
02536 }
02537
02538 Rm0(1) = GRm.x(); Rm0(2) = GRm.y(); Rm0(3) = GRm.z();
02539 Pm0(1) = GPm.x(); Pm0(2) = GPm.y(); Pm0(3) = GPm.z();
02540 Norm(1) = Nl.x(); Norm(2) = Nl.y(); Norm(3) = Nl.z();
02541
02542 CLHEP::HepMatrix T(3,3,1);
02543
02544
02545
02546
02547
02548 double s1 = std::sin(a(1)), c1 = std::cos(a(1));
02549 double s2 = std::sin(a(2)), c2 = std::cos(a(2));
02550 double s3 = std::sin(a(3)), c3 = std::cos(a(3));
02551 T(1,1) = c2 * c3;
02552 T(1,2) = c1 * s3 + s1 * s2 * c3;
02553 T(1,3) = s1 * s3 - c1 * s2 * c3;
02554 T(2,1) = -c2 * s3;
02555 T(2,2) = c1 * c3 - s1 * s2 * s3;
02556 T(2,3) = s1 * c3 + c1 * s2 * s3;
02557 T(3,1) = s2;
02558 T(3,2) = -s1 * c2;
02559 T(3,3) = c1 * c2;
02560
02561
02562
02563 CLHEP::HepMatrix Ti = T.T();
02564 CLHEP::HepVector di = -Ti * d;
02565
02566 CLHEP::HepVector ex0(3,0), ey0(3,0), ez0(3,0), ex(3,0), ey(3,0), ez(3,0);
02567 ex0(1) = 1.; ey0(2) = 1.; ez0(3) = 1.;
02568 ex = Ti*ex0; ey = Ti*ey0; ez = Ti*ez0;
02569
02570 CLHEP::HepVector TiN = Ti * Norm;
02571
02572
02573
02574
02575 double si = CLHEP_dot((Ti*Rm0 - Rm0 + di), TiN) / CLHEP_dot(Pm0, TiN);
02576 Rm1(1) = CLHEP_dot(ex, Rm0 + si*Pm0 - di);
02577 Rm1(2) = CLHEP_dot(ey, Rm0 + si*Pm0 - di);
02578 Rm1(3) = CLHEP_dot(ez, Rm0 + si*Pm0 - di);
02579 Pm1 = T * Pm0;
02580
02581 Rm = GlobalVector(Rm1(1), Rm1(2), Rm1(3));
02582
02583 Pm = GlobalVector(CLHEP_dot(Pm0,ex), CLHEP_dot(Pm0, ey), CLHEP_dot(Pm0,ez));
02584
02585 if(debug_){
02586
02587 std::cout<<" ----- Pm "<<Pm<<std::endl;
02588 std::cout<<" T*Pm0 "<<Pm1.T()<<std::endl;
02589
02590
02591
02592 CLHEP::HepVector Rml = Rm1(1)*ex + Rm1(2)*ey + Rm1(3)*ez + di;
02593
02594 CLHEP::HepVector Rt0(3);
02595 Rt0(1)=Rt.x(); Rt0(2)=Rt.y(); Rt0(3)=Rt.z();
02596
02597
02598 double sl = CLHEP_dot(T*(Rt0 - Rml), T*Norm) / CLHEP_dot(Ti * Pm1, Norm);
02599 CLHEP::HepVector Rl = Rml + sl*(Ti*Pm1);
02600
02601
02602
02603
02604 double A = CLHEP_dot(Ti*Pm1, Norm);
02605 double B = CLHEP_dot(Rt0, Norm) - CLHEP_dot(Ti*Rm1, Norm)
02606 + CLHEP_dot(Ti*d, Norm);
02607 double s = B/A;
02608 CLHEP::HepVector Dr = Ti*Rm1 - Ti*d + s*(Ti*Pm1) - Rm0;
02609 CLHEP::HepVector Dp = Ti*Pm1 - Pm0;
02610
02611 CLHEP::HepVector TiR = Ti * Rm0 + di;
02612 CLHEP::HepVector Ri = T * TiR + d;
02613
02614 std::cout<<" --TTi-0 "<<(Ri-Rm0).T()<<std::endl;
02615 std::cout<<" -- Dr "<<Dr.T()<<endl;
02616 std::cout<<" -- Dp "<<Dp.T()<<endl;
02617
02618
02619
02620
02621
02622
02623
02624 std::cout<<" -- si sl s "<<si<<" "<<sl<<" "<<s<<std::endl;
02625
02626
02627 std::cout<<" -- Rml-(Rm0+si*Pm0) "<<(Rml - Rm0 - si*Pm0).T()<<std::endl;
02628
02629
02630
02631
02632
02633
02634
02635 std::cout<<" --Rl-0 "<<(Rl-Rm0).T()<<std::endl;
02636
02637 }
02638
02639 return;
02640
02641 }
02642
02643
02644 void
02645 GlobalTrackerMuonAlignment::misalignMuonL(GlobalVector& GRm, GlobalVector& GPm, GlobalVector& Nl,
02646 GlobalVector& Rt, GlobalVector& Rm, GlobalVector& Pm,
02647 AlgebraicVector4& Vm,
02648 TrajectoryStateOnSurface& tsosTrack,
02649 TrajectoryStateOnSurface& tsosMuon)
02650 {
02651 CLHEP::HepVector d(3,0), a(3,0), Norm(3), Rm0(3), Pm0(3), Rm1(3), Pm1(3),
02652 Rmi(3), Pmi(3);
02653
02654 d(1)=0.0; d(2)=0.0; d(3)=0.0; a(1)=0.0000; a(2)=0.0000; a(3)=0.0000;
02655
02656
02657
02658
02659
02660
02661
02662
02663
02664
02665
02666
02667
02668 MuGlShift = d; MuGlAngle = a;
02669
02670 if((d(1) == 0.) & (d(2) == 0.) & (d(3) == 0.) &
02671 (a(1) == 0.) & (a(2) == 0.) & (a(3) == 0.)){
02672 Rm = GRm;
02673 Pm = GPm;
02674 AlgebraicVector4 Vm0;
02675 Vm = AlgebraicVector4(tsosMuon.localParameters().vector()(1),
02676 tsosMuon.localParameters().vector()(2),
02677 tsosMuon.localParameters().vector()(3),
02678 tsosMuon.localParameters().vector()(4));
02679 if(debug_)
02680 std::cout<<" ...... without misalignment "<<std::endl;
02681 return;
02682 }
02683
02684 Rm0(1) = GRm.x(); Rm0(2) = GRm.y(); Rm0(3) = GRm.z();
02685 Pm0(1) = GPm.x(); Pm0(2) = GPm.y(); Pm0(3) = GPm.z();
02686 Norm(1) = Nl.x(); Norm(2) = Nl.y(); Norm(3) = Nl.z();
02687
02688 CLHEP::HepMatrix T(3,3,1);
02689
02690
02691
02692
02693
02694 double s1 = std::sin(a(1)), c1 = std::cos(a(1));
02695 double s2 = std::sin(a(2)), c2 = std::cos(a(2));
02696 double s3 = std::sin(a(3)), c3 = std::cos(a(3));
02697 T(1,1) = c2 * c3;
02698 T(1,2) = c1 * s3 + s1 * s2 * c3;
02699 T(1,3) = s1 * s3 - c1 * s2 * c3;
02700 T(2,1) = -c2 * s3;
02701 T(2,2) = c1 * c3 - s1 * s2 * s3;
02702 T(2,3) = s1 * c3 + c1 * s2 * s3;
02703 T(3,1) = s2;
02704 T(3,2) = -s1 * c2;
02705 T(3,3) = c1 * c2;
02706
02707
02708
02709
02710 CLHEP::HepMatrix Ti = T.T();
02711 CLHEP::HepVector di = -Ti * d;
02712
02713 CLHEP::HepVector ex0(3,0), ey0(3,0), ez0(3,0), ex(3,0), ey(3,0), ez(3,0);
02714 ex0(1) = 1.; ey0(2) = 1.; ez0(3) = 1.;
02715 ex = Ti*ex0; ey = Ti*ey0; ez = Ti*ez0;
02716
02717 CLHEP::HepVector TiN = Ti * Norm;
02718
02719
02720
02721
02722 double si = CLHEP_dot((Ti*Rm0 - Rm0 + di), TiN) / CLHEP_dot(Pm0, TiN);
02723 Rm1(1) = CLHEP_dot(ex, Rm0 + si*Pm0 - di);
02724 Rm1(2) = CLHEP_dot(ey, Rm0 + si*Pm0 - di);
02725 Rm1(3) = CLHEP_dot(ez, Rm0 + si*Pm0 - di);
02726 Pm1 = T * Pm0;
02727
02728
02729 Rm = GlobalVector(Rm1(1), Rm1(2), Rm1(3));
02730
02731 Pm = GlobalVector(CLHEP_dot(Pm0,ex), CLHEP_dot(Pm0, ey), CLHEP_dot(Pm0,ez));
02732
02733
02734 const Surface& refSurface = tsosMuon.surface();
02735 CLHEP::HepMatrix Tl(3,3,0);
02736
02737 Tl(1,1) = refSurface.rotation().xx();
02738 Tl(1,2) = refSurface.rotation().xy();
02739 Tl(1,3) = refSurface.rotation().xz();
02740
02741 Tl(2,1) = refSurface.rotation().yx();
02742 Tl(2,2) = refSurface.rotation().yy();
02743 Tl(2,3) = refSurface.rotation().yz();
02744
02745 Tl(3,1) = refSurface.rotation().zx();
02746 Tl(3,2) = refSurface.rotation().zy();
02747 Tl(3,3) = refSurface.rotation().zz();
02748
02749 CLHEP::HepVector R0(3,0), newR0(3,0), newRl(3,0), newPl(3,0);
02750 R0(1) = refSurface.position().x();
02751 R0(2) = refSurface.position().y();
02752 R0(3) = refSurface.position().z();
02753
02754 newRl = Tl * (Rm1 - R0);
02755 newPl = Tl * Pm1;
02756
02757 Vm(0) = newPl(1)/newPl(3);
02758 Vm(1) = newPl(2)/newPl(3);
02759 Vm(2) = newRl(1);
02760 Vm(3) = newRl(2);
02761
02762 #ifdef CHECK_DERIVATIVES_LOCAL_VS_ANGLES
02763 float del = 0.0001;
02764
02765 int ii = 2; T(3,1) -=del; T(1,3) +=del;
02766
02767 AlgebraicVector4 Vm0 = Vm;
02768 CLHEP::HepVector Rm10 = Rm1, Pm10 = Pm1;;
02769 CLHEP::HepMatrix newTl = Tl*T;
02770 Ti = T.T();
02771 di = -Ti * d;
02772 ex = Ti*ex0; ey = Ti*ey0; ez = Ti*ez0;
02773 TiN = Ti * Norm;
02774 si = CLHEP_dot((Ti*Rm0 - Rm0 + di), TiN) / CLHEP_dot(Pm0, TiN);
02775 Rm1(1) = CLHEP_dot(ex, Rm0 + si*Pm0 - di);
02776 Rm1(2) = CLHEP_dot(ey, Rm0 + si*Pm0 - di);
02777 Rm1(3) = CLHEP_dot(ez, Rm0 + si*Pm0 - di);
02778 Pm1 = T * Pm0;
02779
02780 newRl = Tl * (Rm1 - R0);
02781 newPl = Tl * Pm1;
02782
02783 Vm(0) = newPl(1)/newPl(3);
02784 Vm(1) = newPl(2)/newPl(3);
02785 Vm(2) = newRl(1);
02786 Vm(3) = newRl(2);
02787 std::cout<<" ========= dVm/da"<<ii<<" "<<(Vm-Vm0)*(1./del)<<std::endl;
02788 #endif
02789
02790 if(debug_){
02791
02792
02793 std::cout<<" ex "<<(Tl.T()*ex0).T()<<std::endl;
02794 std::cout<<" ey "<<(Tl.T()*ey0).T()<<std::endl;
02795 std::cout<<" ez "<<(Tl.T()*ez0).T()<<std::endl;
02796
02797
02798 std::cout<<" pxyz/p gl 0 "<<(Pm0/Pm0.norm()).T()<<std::endl;
02799 std::cout<<" pxyz/p loc0 "<<(Tl*Pm0/(Tl*Pm0).norm()).T()<<std::endl;
02800 std::cout<<" pxyz/p glob "<<(Pm1/Pm1.norm()).T()<<std::endl;
02801 std::cout<<" pxyz/p loc "<<(newPl/newPl.norm()).T()<<std::endl;
02802
02803
02804
02805 std::cout<<" ---- phi g0 g1 l0 l1 "
02806 <<atan2(Pm0(2),Pm0(1))<<" "<<atan2(Pm1(2),Pm1(1))<<" "
02807 <<atan2((Tl*Pm0)(2),(Tl*Pm0)(1))<<" "
02808 <<atan2(newPl(2),newPl(1))<<std::endl;
02809 std::cout<<" ---- angl Gl Loc "<<atan2(Pm1(2),Pm1(1))-atan2(Pm0(2),Pm0(1))<<" "
02810 <<atan2(newPl(2),newPl(1))-atan2((Tl*Pm0)(2),(Tl*Pm0)(1))<<" "
02811 <<atan2(newPl(3),newPl(2))-atan2((Tl*Pm0)(3),(Tl*Pm0)(2))<<" "
02812 <<atan2(newPl(1),newPl(3))-atan2((Tl*Pm0)(1),(Tl*Pm0)(3))<<" "<<std::endl;
02813 }
02814
02815 if(debug_){
02816 CLHEP::HepMatrix newTl(3,3,0);
02817 CLHEP::HepVector R0(3,0), newRl(3,0), newPl(3,0);
02818 newTl = Tl * Ti.T();
02819 newR0 = Ti * R0 + di;
02820
02821 std::cout<<" N "<<Norm.T()<<std::endl;
02822 std::cout<<" dtxl yl "<<Vm(0)-tsosMuon.localParameters().vector()(1)<<" "
02823 <<Vm(1)-tsosMuon.localParameters().vector()(2)<<std::endl;
02824 std::cout<<" dXl dYl "<<Vm(2)-tsosMuon.localParameters().vector()(3)<<" "
02825 <<Vm(3)-tsosMuon.localParameters().vector()(4)<<std::endl;
02826 std::cout<<" localM "<<tsosMuon.localParameters().vector()<<std::endl;
02827 std::cout<<" Vm "<<Vm<<std::endl;
02828
02829
02830 CLHEP::HepVector Rvc(3,0), Pvc(3,0), Rvg(3,0), Pvg(3,0);
02831 Rvc(1) = Vm(2); Rvc(2) = Vm(3);
02832 Pvc(3) = tsosMuon.localParameters().pzSign() * Pm0.norm() /
02833 sqrt(1 + Vm(0)*Vm(0) + Vm(1)*Vm(1));
02834 Pvc(1) = Pvc(3) * Vm(0);
02835 Pvc(2) = Pvc(3) * Vm(1);
02836
02837 Rvg = newTl.T() * Rvc + newR0;
02838 Pvg = newTl.T() * Pvc;
02839
02840 std::cout<<" newPl "<<newPl.T()<<std::endl;
02841 std::cout<<" Pvc "<<Pvc.T()<<std::endl;
02842 std::cout<<" Rvg "<<Rvg.T()<<std::endl;
02843 std::cout<<" Rm1 "<<Rm1.T()<<std::endl;
02844 std::cout<<" Pvg "<<Pvg.T()<<std::endl;
02845 std::cout<<" Pm1 "<<Pm1.T()<<std::endl;
02846 }
02847
02848 if(debug_){
02849
02850 std::cout<<" ----- Pm "<<Pm<<std::endl;
02851 std::cout<<" T*Pm0 "<<Pm1.T()<<std::endl;
02852
02853
02854
02855 CLHEP::HepVector Rml = Rm1(1)*ex + Rm1(2)*ey + Rm1(3)*ez + di;
02856
02857 CLHEP::HepVector Rt0(3);
02858 Rt0(1)=Rt.x(); Rt0(2)=Rt.y(); Rt0(3)=Rt.z();
02859
02860
02861 double sl = CLHEP_dot(T*(Rt0 - Rml), T*Norm) / CLHEP_dot(Ti * Pm1, Norm);
02862 CLHEP::HepVector Rl = Rml + sl*(Ti*Pm1);
02863
02864
02865
02866
02867 double A = CLHEP_dot(Ti*Pm1, Norm);
02868 double B = CLHEP_dot(Rt0, Norm) - CLHEP_dot(Ti*Rm1, Norm)
02869 + CLHEP_dot(Ti*d, Norm);
02870 double s = B/A;
02871 CLHEP::HepVector Dr = Ti*Rm1 - Ti*d + s*(Ti*Pm1) - Rm0;
02872 CLHEP::HepVector Dp = Ti*Pm1 - Pm0;
02873
02874 CLHEP::HepVector TiR = Ti * Rm0 + di;
02875 CLHEP::HepVector Ri = T * TiR + d;
02876
02877 std::cout<<" --TTi-0 "<<(Ri-Rm0).T()<<std::endl;
02878 std::cout<<" -- Dr "<<Dr.T()<<endl;
02879 std::cout<<" -- Dp "<<Dp.T()<<endl;
02880
02881
02882
02883
02884
02885
02886
02887 std::cout<<" -- si sl s "<<si<<" "<<sl<<" "<<s<<std::endl;
02888
02889
02890 std::cout<<" -- Rml-(Rm0+si*Pm0) "<<(Rml - Rm0 - si*Pm0).T()<<std::endl;
02891
02892
02893
02894
02895
02896
02897
02898 std::cout<<" --Rl-0 "<<(Rl-Rm0).T()<<std::endl;
02899
02900 }
02901
02902 return;
02903
02904 }
02905
02906
02907
02908 void
02909 GlobalTrackerMuonAlignment::trackFitter(reco::TrackRef alongTr, reco::TransientTrack& alongTTr,
02910 PropagationDirection direction,
02911 TrajectoryStateOnSurface& trackFittedTSOS)
02912 {
02913 bool info = false;
02914 bool smooth = false;
02915
02916 if(info)
02917 std::cout<<" ......... start of trackFitter ......... "<<std::endl;
02918 if(false && info){
02919 this->debugTrackHit(" Tracker track hits ", alongTr);
02920 this->debugTrackHit(" Tracker TransientTrack hits ", alongTTr);
02921 }
02922
02923 edm::OwnVector<TrackingRecHit> recHit;
02924 DetId trackDetId = DetId(alongTr->innerDetId());
02925 for (trackingRecHit_iterator i=alongTr->recHitsBegin();i!=alongTr->recHitsEnd(); i++) {
02926 recHit.push_back((*i)->clone());
02927 }
02928 if(direction != alongMomentum)
02929 {
02930 edm::OwnVector<TrackingRecHit> recHitAlong = recHit;
02931 recHit.clear();
02932 for(unsigned int ihit = recHitAlong.size()-1; ihit+1>0; ihit--){
02933 recHit.push_back(recHitAlong[ihit]);
02934 }
02935 recHitAlong.clear();
02936 }
02937 if(info)
02938 std::cout<<" ... Own recHit.size() "<<recHit.size()<<" "<<trackDetId.rawId()<<std::endl;
02939
02940
02941 TransientTrackingRecHit::RecHitContainer recHitTrack;
02942 for (trackingRecHit_iterator i=alongTr->recHitsBegin();i!=alongTr->recHitsEnd(); i++) {
02943 if((*i)->isValid()){
02944 recHitTrack.push_back(TTRHBuilder->build(&**i ));}
02945 else{
02946 recHitTrack.push_back(TTRHBuilder->build(&**i ));}
02947 }
02948
02949 if(info)
02950 std::cout<<" ... recHitTrack.size() "<<recHit.size()<<" "<<trackDetId.rawId()
02951 <<" ======> recHitMu "<<std::endl;
02952
02953 MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu = recHitTrack;
02954
02955
02956
02957
02958 if(info)
02959 std::cout<<" ...along.... recHitMu.size() "<<recHitMu.size()<<std::endl;
02960 if(direction != alongMomentum){
02961 MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMuAlong = recHitMu;
02962 recHitMu.clear();
02963 for(unsigned int ihit = recHitMuAlong.size()-1; ihit+1>0; ihit--){
02964 recHitMu.push_back(recHitMuAlong[ihit]);
02965 }
02966 recHitMuAlong.clear();
02967 }
02968 if(info)
02969 std::cout<<" ...opposite... recHitMu.size() "<<recHitMu.size()<<std::endl;
02970
02971 TrajectoryStateOnSurface firstTSOS;
02972 if(direction == alongMomentum) firstTSOS = alongTTr.innermostMeasurementState();
02973 else firstTSOS = alongTTr.outermostMeasurementState();
02974
02975 AlgebraicSymMatrix55 CovLoc;
02976 CovLoc(0,0) = 1. * firstTSOS.localError().matrix()(0,0);
02977 CovLoc(1,1) = 10. * firstTSOS.localError().matrix()(1,1);
02978 CovLoc(2,2) = 10. * firstTSOS.localError().matrix()(2,2);
02979 CovLoc(3,3) = 100. * firstTSOS.localError().matrix()(3,3);
02980 CovLoc(4,4) = 100. * firstTSOS.localError().matrix()(4,4);
02981 TrajectoryStateOnSurface initialTSOS(firstTSOS.localParameters(), LocalTrajectoryError(CovLoc),
02982 firstTSOS.surface(), &*magneticField_);
02983
02984 PTrajectoryStateOnDet PTraj =
02985 trajectoryStateTransform::persistentState(initialTSOS, trackDetId.rawId());
02986
02987 const TrajectorySeed seedT(PTraj, recHit, direction);
02988
02989 std::vector<Trajectory> trajVec;
02990 if(direction == alongMomentum) trajVec = theFitter->fit(seedT, recHitMu, initialTSOS);
02991 else trajVec = theFitterOp->fit(seedT, recHitMu, initialTSOS);
02992
02993 if(info){
02994 this->debugTrajectorySOSv(" innermostTSOS of TransientTrack",
02995 alongTTr.innermostMeasurementState());
02996 this->debugTrajectorySOSv(" outermostTSOS of TransientTrack",
02997 alongTTr.outermostMeasurementState());
02998 this->debugTrajectorySOS(" initialTSOS for theFitter ", initialTSOS);
02999 std::cout<<" . . . . . trajVec.size() "<<trajVec.size()<<std::endl;
03000 if(trajVec.size()) this->debugTrajectory(" theFitter trajectory ", trajVec[0]);
03001 }
03002
03003 if(!smooth){
03004 if(trajVec.size()) trackFittedTSOS = trajVec[0].lastMeasurement().updatedState();}
03005 else{
03006 std::vector<Trajectory> trajSVec;
03007 if (trajVec.size()){
03008 if(direction == alongMomentum) trajSVec = theSmoother->trajectories(*(trajVec.begin()));
03009 else trajSVec = theSmootherOp->trajectories(*(trajVec.begin()));
03010 }
03011 if(info) std::cout<<" . . . . trajSVec.size() "<<trajSVec.size()<<std::endl;
03012 if(trajSVec.size()) this->debugTrajectorySOSv("smoothed geom InnermostState",
03013 trajSVec[0].geometricalInnermostState());
03014 if(trajSVec.size()) trackFittedTSOS = trajSVec[0].geometricalInnermostState();
03015 }
03016
03017 if(info) this->debugTrajectorySOSv(" trackFittedTSOS ", trackFittedTSOS);
03018
03019 }
03020
03021
03022 void
03023 GlobalTrackerMuonAlignment::muonFitter(reco::TrackRef alongTr, reco::TransientTrack& alongTTr,
03024 PropagationDirection direction,
03025 TrajectoryStateOnSurface& trackFittedTSOS)
03026 {
03027 bool info = false;
03028 bool smooth = false;
03029
03030 if(info) std::cout<<" ......... start of muonFitter ........"<<std::endl;
03031 if(false && info){
03032 this->debugTrackHit(" Muon track hits ", alongTr);
03033 this->debugTrackHit(" Muon TransientTrack hits ", alongTTr);
03034 }
03035
03036 edm::OwnVector<TrackingRecHit> recHit;
03037 DetId trackDetId = DetId(alongTr->innerDetId());
03038 for (trackingRecHit_iterator i=alongTr->recHitsBegin();i!=alongTr->recHitsEnd(); i++) {
03039 recHit.push_back((*i)->clone());
03040 }
03041 if(direction != alongMomentum)
03042 {
03043 edm::OwnVector<TrackingRecHit> recHitAlong = recHit;
03044 recHit.clear();
03045 for(unsigned int ihit = recHitAlong.size()-1; ihit+1>0; ihit--){
03046 recHit.push_back(recHitAlong[ihit]);
03047 }
03048 recHitAlong.clear();
03049 }
03050 if(info)
03051 std::cout<<" ... Own recHit.size() DetId==Muon "<<recHit.size()<<" "<<trackDetId.rawId()<<std::endl;
03052
03053 MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMu =
03054 MuRHBuilder->build(alongTTr.recHitsBegin(), alongTTr.recHitsEnd());
03055 if(info)
03056 std::cout<<" ...along.... recHitMu.size() "<<recHitMu.size()<<std::endl;
03057 if(direction != alongMomentum){
03058 MuonTransientTrackingRecHitBuilder::ConstRecHitContainer recHitMuAlong = recHitMu;
03059 recHitMu.clear();
03060 for(unsigned int ihit = recHitMuAlong.size()-1; ihit+1>0; ihit--){
03061 recHitMu.push_back(recHitMuAlong[ihit]);
03062 }
03063 recHitMuAlong.clear();
03064 }
03065 if(info)
03066 std::cout<<" ...opposite... recHitMu.size() "<<recHitMu.size()<<std::endl;
03067
03068 TrajectoryStateOnSurface firstTSOS;
03069 if(direction == alongMomentum) firstTSOS = alongTTr.innermostMeasurementState();
03070 else firstTSOS = alongTTr.outermostMeasurementState();
03071
03072 AlgebraicSymMatrix55 CovLoc;
03073 CovLoc(0,0) = 1. * firstTSOS.localError().matrix()(0,0);
03074 CovLoc(1,1) = 10. * firstTSOS.localError().matrix()(1,1);
03075 CovLoc(2,2) = 10. * firstTSOS.localError().matrix()(2,2);
03076 CovLoc(3,3) = 100. * firstTSOS.localError().matrix()(3,3);
03077 CovLoc(4,4) = 100. * firstTSOS.localError().matrix()(4,4);
03078 TrajectoryStateOnSurface initialTSOS(firstTSOS.localParameters(), LocalTrajectoryError(CovLoc),
03079 firstTSOS.surface(), &*magneticField_);
03080
03081 PTrajectoryStateOnDet PTraj =
03082 trajectoryStateTransform::persistentState(initialTSOS, trackDetId.rawId());
03083
03084 const TrajectorySeed seedT(PTraj, recHit, direction);
03085
03086 std::vector<Trajectory> trajVec;
03087 if(direction == alongMomentum) trajVec = theFitter->fit(seedT, recHitMu, initialTSOS);
03088 else trajVec = theFitterOp->fit(seedT, recHitMu, initialTSOS);
03089
03090 if(info){
03091 this->debugTrajectorySOSv(" innermostTSOS of TransientTrack",
03092 alongTTr.innermostMeasurementState());
03093 this->debugTrajectorySOSv(" outermostTSOS of TransientTrack",
03094 alongTTr.outermostMeasurementState());
03095 this->debugTrajectorySOS(" initialTSOS for theFitter ", initialTSOS);
03096 std::cout<<" . . . . . trajVec.size() "<<trajVec.size()<<std::endl;
03097 if(trajVec.size()) this->debugTrajectory(" theFitter trajectory ", trajVec[0]);
03098 }
03099
03100 if(!smooth){
03101 if(trajVec.size()) trackFittedTSOS = trajVec[0].lastMeasurement().updatedState();}
03102 else{
03103 std::vector<Trajectory> trajSVec;
03104 if (trajVec.size()){
03105 if(direction == alongMomentum) trajSVec = theSmoother->trajectories(*(trajVec.begin()));
03106 else trajSVec = theSmootherOp->trajectories(*(trajVec.begin()));
03107 }
03108 if(info) std::cout<<" . . . . trajSVec.size() "<<trajSVec.size()<<std::endl;
03109 if(trajSVec.size()) this->debugTrajectorySOSv("smoothed geom InnermostState",
03110 trajSVec[0].geometricalInnermostState());
03111 if(trajSVec.size()) trackFittedTSOS = trajSVec[0].geometricalInnermostState();
03112 }
03113
03114 }
03115
03116
03117
03118 void GlobalTrackerMuonAlignment::debugTrackHit(const std::string title, TransientTrack& alongTr)
03119 {
03120 std::cout<<" ------- "<<title<<" --------"<<std::endl;
03121 int nHit = 1;
03122 TransientTrackingRecHit::RecHitContainer recHit;
03123 for (trackingRecHit_iterator i=alongTr.recHitsBegin();i!=alongTr.recHitsEnd(); i++) {
03124 std::cout<<" Hit "<<nHit++<<" DetId "<<(*i)->geographicalId().det();
03125 if( (*i)->geographicalId().det() == DetId::Tracker ) std::cout<<" Tracker ";
03126 else if ( (*i)->geographicalId().det() == DetId::Muon ) std::cout<<" Muon ";
03127 else std::cout<<" Unknown ";
03128 if(!(*i)->isValid()) std::cout<<" not valid "<<std::endl;
03129 else std::cout<<std::endl;
03130 }
03131 }
03132
03133
03134
03135 void GlobalTrackerMuonAlignment::debugTrackHit(const std::string title, reco::TrackRef alongTr)
03136 {
03137 std::cout<<" ------- "<<title<<" --------"<<std::endl;
03138 int nHit = 1;
03139 TransientTrackingRecHit::RecHitContainer recHit;
03140 for (trackingRecHit_iterator i=alongTr->recHitsBegin();i!=alongTr->recHitsEnd(); i++) {
03141 std::cout<<" Hit "<<nHit++<<" DetId "<<(*i)->geographicalId().det();
03142 if( (*i)->geographicalId().det() == DetId::Tracker ) std::cout<<" Tracker ";
03143 else if ( (*i)->geographicalId().det() == DetId::Muon ) std::cout<<" Muon ";
03144 else std::cout<<" Unknown ";
03145 if(!(*i)->isValid()) std::cout<<" not valid "<<std::endl;
03146 else std::cout<<std::endl;
03147 }
03148 }
03149
03150
03151 void GlobalTrackerMuonAlignment::debugTrajectorySOS(const std::string title,
03152 TrajectoryStateOnSurface& trajSOS)
03153 {
03154 std::cout<<" --- "<<title<<" --- "<<std::endl;
03155 if(!trajSOS.isValid()) {std::cout<<" Not valid !!!! "<<std::endl; return;}
03156 std::cout<<" R |p| "<<trajSOS.globalPosition().perp()<<" "
03157 <<trajSOS.globalMomentum().mag()<<" charge "<<trajSOS.charge()<<std::endl;
03158 std::cout<<" x p "<<trajSOS.globalParameters().vector()(0)<<" "
03159 <<trajSOS.globalParameters().vector()(1)<<" "
03160 <<trajSOS.globalParameters().vector()(2)<<" "
03161 <<trajSOS.globalParameters().vector()(3)<<" "
03162 <<trajSOS.globalParameters().vector()(4)<<" "
03163 <<trajSOS.globalParameters().vector()(5)<<std::endl;
03164 std::cout<<" +/- "<<sqrt(trajSOS.cartesianError().matrix()(0,0))<<" "
03165 <<sqrt(trajSOS.cartesianError().matrix()(1,1))<<" "
03166 <<sqrt(trajSOS.cartesianError().matrix()(2,2))<<" "
03167 <<sqrt(trajSOS.cartesianError().matrix()(3,3))<<" "
03168 <<sqrt(trajSOS.cartesianError().matrix()(4,4))<<" "
03169 <<sqrt(trajSOS.cartesianError().matrix()(5,5))<<std::endl;
03170 std::cout<<" q/p dxy/dz xy "<<trajSOS.localParameters().vector()(0)<<" "
03171 <<trajSOS.localParameters().vector()(1)<<" "
03172 <<trajSOS.localParameters().vector()(2)<<" "
03173 <<trajSOS.localParameters().vector()(3)<<" "
03174 <<trajSOS.localParameters().vector()(4)<<std::endl;
03175 std::cout<<" +/- error "<<sqrt(trajSOS.localError().matrix()(0,0))<<" "
03176 <<sqrt(trajSOS.localError().matrix()(1,1))<<" "
03177 <<sqrt(trajSOS.localError().matrix()(2,2))<<" "
03178 <<sqrt(trajSOS.localError().matrix()(3,3))<<" "
03179 <<sqrt(trajSOS.localError().matrix()(4,4))<<" "<<std::endl;
03180 }
03181
03182
03183 void GlobalTrackerMuonAlignment::debugTrajectorySOSv(const std::string title,
03184 TrajectoryStateOnSurface trajSOS)
03185 {
03186 std::cout<<" --- "<<title<<" --- "<<std::endl;
03187 if(!trajSOS.isValid()) {std::cout<<" Not valid !!!! "<<std::endl; return;}
03188 std::cout<<" R |p| "<<trajSOS.globalPosition().perp()<<" "
03189 <<trajSOS.globalMomentum().mag()<<" charge "<<trajSOS.charge()<<std::endl;
03190 std::cout<<" x p "<<trajSOS.globalParameters().vector()(0)<<" "
03191 <<trajSOS.globalParameters().vector()(1)<<" "
03192 <<trajSOS.globalParameters().vector()(2)<<" "
03193 <<trajSOS.globalParameters().vector()(3)<<" "
03194 <<trajSOS.globalParameters().vector()(4)<<" "
03195 <<trajSOS.globalParameters().vector()(5)<<std::endl;
03196 std::cout<<" +/- "<<sqrt(trajSOS.cartesianError().matrix()(0,0))<<" "
03197 <<sqrt(trajSOS.cartesianError().matrix()(1,1))<<" "
03198 <<sqrt(trajSOS.cartesianError().matrix()(2,2))<<" "
03199 <<sqrt(trajSOS.cartesianError().matrix()(3,3))<<" "
03200 <<sqrt(trajSOS.cartesianError().matrix()(4,4))<<" "
03201 <<sqrt(trajSOS.cartesianError().matrix()(5,5))<<std::endl;
03202 std::cout<<" q/p dxy/dz xy "<<trajSOS.localParameters().vector()(0)<<" "
03203 <<trajSOS.localParameters().vector()(1)<<" "
03204 <<trajSOS.localParameters().vector()(2)<<" "
03205 <<trajSOS.localParameters().vector()(3)<<" "
03206 <<trajSOS.localParameters().vector()(4)<<std::endl;
03207 std::cout<<" +/- error "<<sqrt(trajSOS.localError().matrix()(0,0))<<" "
03208 <<sqrt(trajSOS.localError().matrix()(1,1))<<" "
03209 <<sqrt(trajSOS.localError().matrix()(2,2))<<" "
03210 <<sqrt(trajSOS.localError().matrix()(3,3))<<" "
03211 <<sqrt(trajSOS.localError().matrix()(4,4))<<" "<<std::endl;
03212 }
03213
03214
03215 void GlobalTrackerMuonAlignment::debugTrajectory(const std::string title, Trajectory& traj)
03216 {
03217 std::cout<<"\n"<<" ...... "<<title<<" ...... "<<std::endl;
03218 if(!traj.isValid()) {std::cout<<" Not valid !!!!!!!! "<<std::endl; return;}
03219 std::cout<<" chi2/Nhit "<<traj.chiSquared()<<" / "<<traj.foundHits();
03220 if(traj.direction() == alongMomentum) std::cout<<" alongMomentum >>>>"<<std::endl;
03221 else std::cout<<" oppositeToMomentum <<<<"<<std::endl;
03222 this->debugTrajectorySOSv(" firstMeasurementTSOS ",traj.firstMeasurement().updatedState());
03223
03224 this->debugTrajectorySOSv(" lastMeasurementTSOS ",traj.lastMeasurement().updatedState());
03225
03226 std::cout<<" . . . . . . . . . . . . . . . . . . . . . . . . . . . . \n"<<std::endl;
03227 }
03228
03229
03230
03231 void GlobalTrackerMuonAlignment::writeGlPosRcd(CLHEP::HepVector& paramVec)
03232 {
03233
03234
03235
03236
03237
03238
03239
03240
03241
03242
03243
03244 std::cout<<" paramVector "<<paramVec.T()<<std::endl;
03245
03246 CLHEP::Hep3Vector colX, colY, colZ;
03247
03248 #ifdef NOT_EXACT_ROTATION_MATRIX
03249 colX = CLHEP::Hep3Vector( 1., -paramVec(6), paramVec(5));
03250 colY = CLHEP::Hep3Vector( paramVec(6), 1., -paramVec(4));
03251 colZ = CLHEP::Hep3Vector(-paramVec(5), paramVec(4), 1.);
03252 #else
03253 double s1 = std::sin(paramVec(4)), c1 = std::cos(paramVec(4));
03254 double s2 = std::sin(paramVec(5)), c2 = std::cos(paramVec(5));
03255 double s3 = std::sin(paramVec(6)), c3 = std::cos(paramVec(6));
03256 colX = CLHEP::Hep3Vector( c2 * c3, -c2 * s3, s2);
03257 colY = CLHEP::Hep3Vector(c1 * s3 + s1 * s2 * c3, c1 * c3 - s1 * s2 * s3, -s1 * c2);
03258 colZ = CLHEP::Hep3Vector(s1 * s3 - c1 * s2 * c3, s1 * c3 + c1 * s2 * s3, c1 * c2);
03259 #endif
03260
03261 CLHEP::HepVector param0(6,0);
03262
03263 Alignments* globalPositions = new Alignments();
03264
03265
03266 AlignTransform tracker(iteratorTrackerRcd->translation(),
03267 iteratorTrackerRcd->rotation(),
03268 DetId(DetId::Tracker).rawId());
03269
03270 CLHEP::Hep3Vector posMuGlRcd = iteratorMuonRcd->translation();
03271 CLHEP::HepRotation rotMuGlRcd = iteratorMuonRcd->rotation();
03272 CLHEP::HepEulerAngles angMuGlRcd = iteratorMuonRcd->rotation().eulerAngles();
03273 if(debug_)
03274 std::cout<<" Old muon Rcd Euler angles "<<angMuGlRcd.phi()<<" "<<angMuGlRcd.theta()
03275 <<" "<<angMuGlRcd.psi()<<std::endl;
03276 AlignTransform muon;
03277 if((angMuGlRcd.phi() == 0.) && (angMuGlRcd.theta() == 0.) && (angMuGlRcd.psi() == 0.) &&
03278 (posMuGlRcd.x() == 0.) && (posMuGlRcd.y() == 0.) && (posMuGlRcd.z() == 0.)){
03279 std::cout<<" New muon parameters are stored in Rcd"<<std::endl;
03280
03281 AlignTransform muonNew(AlignTransform::Translation(paramVec(1),
03282 paramVec(2),
03283 paramVec(3)),
03284 AlignTransform::Rotation (colX,
03285 colY,
03286 colZ),
03287 DetId(DetId::Muon).rawId());
03288 muon = muonNew;
03289 }
03290 else if((paramVec(1) == 0.) && (paramVec(2) == 0.) && (paramVec(3) == 0.) &&
03291 (paramVec(4) == 0.) && (paramVec(5) == 0.) && (paramVec(6) == 0.)){
03292 std::cout<<" Old muon parameters are stored in Rcd"<<std::endl;
03293
03294 AlignTransform muonNew(iteratorMuonRcd->translation(),
03295 iteratorMuonRcd->rotation(),
03296 DetId(DetId::Muon).rawId());
03297 muon = muonNew;
03298 }
03299 else{
03300 std::cout<<" New + Old muon parameters are stored in Rcd"<<std::endl;
03301 CLHEP::Hep3Vector posMuGlRcdThis = CLHEP::Hep3Vector(paramVec(1), paramVec(2), paramVec(3));
03302 CLHEP::HepRotation rotMuGlRcdThis = CLHEP::HepRotation(colX, colY, colZ);
03303 CLHEP::Hep3Vector posMuGlRcdNew = rotMuGlRcdThis * posMuGlRcd + posMuGlRcdThis;
03304 CLHEP::HepRotation rotMuGlRcdNew = rotMuGlRcdThis * rotMuGlRcd;
03305
03306 AlignTransform muonNew(posMuGlRcdNew,
03307 rotMuGlRcdNew,
03308 DetId(DetId::Muon).rawId());
03309 muon = muonNew;
03310 }
03311
03312
03313 AlignTransform ecal(iteratorEcalRcd->translation(),
03314 iteratorEcalRcd->rotation(),
03315 DetId(DetId::Ecal).rawId());
03316
03317 AlignTransform hcal(iteratorHcalRcd->translation(),
03318 iteratorHcalRcd->rotation(),
03319 DetId(DetId::Hcal).rawId());
03320
03321 AlignTransform calo(AlignTransform::Translation(param0(1),
03322 param0(2),
03323 param0(3)),
03324 AlignTransform::EulerAngles(param0(4),
03325 param0(5),
03326 param0(6)),
03327 DetId(DetId::Calo).rawId());
03328
03329 std::cout << "Tracker (" << tracker.rawId() << ") at " << tracker.translation()
03330 << " " << tracker.rotation().eulerAngles() << std::endl;
03331 std::cout << tracker.rotation() << std::endl;
03332
03333 std::cout << "Muon (" << muon.rawId() << ") at " << muon.translation()
03334 << " " << muon.rotation().eulerAngles() << std::endl;
03335 std::cout << " rotations angles around x,y,z "
03336 << " ( " << -muon.rotation().zy() << " " << muon.rotation().zx()
03337 << " " << -muon.rotation().yx() << " )" << std::endl;
03338 std::cout << muon.rotation() << std::endl;
03339
03340 std::cout << "Ecal (" << ecal.rawId() << ") at " << ecal.translation()
03341 << " " << ecal.rotation().eulerAngles() << std::endl;
03342 std::cout << ecal.rotation() << std::endl;
03343
03344 std::cout << "Hcal (" << hcal.rawId() << ") at " << hcal.translation()
03345 << " " << hcal.rotation().eulerAngles() << std::endl;
03346 std::cout << hcal.rotation() << std::endl;
03347
03348 std::cout << "Calo (" << calo.rawId() << ") at " << calo.translation()
03349 << " " << calo.rotation().eulerAngles() << std::endl;
03350 std::cout << calo.rotation() << std::endl;
03351
03352 globalPositions->m_align.push_back(tracker);
03353 globalPositions->m_align.push_back(muon);
03354 globalPositions->m_align.push_back(ecal);
03355 globalPositions->m_align.push_back(hcal);
03356 globalPositions->m_align.push_back(calo);
03357
03358 std::cout << "Uploading to the database..." << std::endl;
03359
03360 edm::Service<cond::service::PoolDBOutputService> poolDbService;
03361
03362 if (!poolDbService.isAvailable())
03363 throw cms::Exception("NotAvailable") << "PoolDBOutputService not available";
03364
03365
03366
03367
03368
03369
03370 poolDbService->writeOne<Alignments>(&(*globalPositions),
03371 poolDbService->currentTime(),
03372
03373 "GlobalPositionRcd");
03374 std::cout << "done!" << std::endl;
03375
03376 return;
03377 }
03378
03379
03380
03381 DEFINE_FWK_MODULE(GlobalTrackerMuonAlignment);
1.7.3