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