File indexing completed on 2023-03-17 10:39:17
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016 #include "FWCore/Framework/interface/Frameworkfwd.h"
0017 #include "FWCore/Framework/interface/one/EDProducer.h"
0018 #include "FWCore/Framework/interface/ESHandle.h"
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/Framework/interface/EventSetup.h"
0021 #include "FWCore/Framework/interface/MakerMacros.h"
0022 #include "FWCore/Framework/interface/Run.h"
0023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0024 #include "FWCore/Utilities/interface/InputTag.h"
0025 #include "FWCore/ServiceRegistry/interface/Service.h"
0026 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0027
0028 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
0029 #include <Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h>
0030 #include <Geometry/Records/interface/TrackerDigiGeometryRecord.h>
0031 #include <MagneticField/Engine/interface/MagneticField.h>
0032 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0033
0034
0035
0036 #include "DataFormats/Provenance/interface/BranchType.h"
0037
0038
0039 #include "DataFormats/Alignment/interface/TkLasBeam.h"
0040 #include "DataFormats/Alignment/interface/TkFittedLasBeam.h"
0041
0042
0043 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0044 #include "Alignment/LaserAlignment/interface/TsosVectorCollection.h"
0045 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0046 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0047
0048 #include <iostream>
0049 #include "TMinuit.h"
0050 #include "TGraphErrors.h"
0051 #include "TF1.h"
0052 #include "TH1.h"
0053 #include "TH2.h"
0054
0055 using namespace edm;
0056 using namespace std;
0057
0058
0059
0060
0061
0062 class TkLasBeamFitter : public edm::one::EDProducer<edm::EndRunProducer> {
0063 public:
0064 explicit TkLasBeamFitter(const edm::ParameterSet &config);
0065 ~TkLasBeamFitter() override;
0066
0067
0068 void produce(edm::Event &event, const edm::EventSetup &setup) override;
0069
0070 void endRunProduce(edm::Run &run, const edm::EventSetup &setup) override;
0071
0072
0073 private:
0074
0075
0076 void getLasBeams(TkFittedLasBeam &beam, vector<TrajectoryStateOnSurface> &tsosLas);
0077 void getLasHits(TkFittedLasBeam &beam,
0078 const SiStripLaserRecHit2D &hit,
0079 vector<const GeomDetUnit *> &gd,
0080 vector<GlobalPoint> &globHit,
0081 unsigned int &hitsAtTecPlus);
0082
0083
0084
0085
0086 static double tecPlusFunction(double *x, double *par);
0087 static double tecMinusFunction(double *x, double *par);
0088 static double atFunction(double *x, double *par);
0089
0090 void fitter(TkFittedLasBeam &beam,
0091 AlgebraicSymMatrix &covMatrix,
0092 unsigned int &hitsAtTecPlus,
0093 unsigned int &nFitParams,
0094 std::vector<double> &hitPhi,
0095 std::vector<double> &hitPhiError,
0096 std::vector<double> &hitZprimeError,
0097 double &zMin,
0098 double &zMax,
0099 double &bsAngleParam,
0100 double &offset,
0101 double &offsetError,
0102 double &slope,
0103 double &slopeError,
0104 double &phiAtMinusParam,
0105 double &phiAtPlusParam,
0106 double &atThetaSplitParam);
0107
0108 void trackPhi(TkFittedLasBeam &beam,
0109 unsigned int &hit,
0110 double &trackPhi,
0111 double &trackPhiRef,
0112 double &offset,
0113 double &slope,
0114 double &bsAngleParam,
0115 double &phiAtMinusParam,
0116 double &phiAtPlusParam,
0117 double &atThetaSplitParam,
0118 std::vector<GlobalPoint> &globHit);
0119
0120 void globalTrackPoint(TkFittedLasBeam &beam,
0121 unsigned int &hit,
0122 unsigned int &hitsAtTecPlus,
0123 double &trackPhi,
0124 double &trackPhiRef,
0125 std::vector<GlobalPoint> &globHit,
0126 std::vector<GlobalPoint> &globPtrack,
0127 GlobalPoint &globPref,
0128 std::vector<double> &hitPhiError);
0129
0130 void buildTrajectory(TkFittedLasBeam &beam,
0131 unsigned int &hit,
0132 vector<const GeomDetUnit *> &gd,
0133 std::vector<GlobalPoint> &globPtrack,
0134 vector<TrajectoryStateOnSurface> &tsosLas,
0135 GlobalPoint &globPref);
0136
0137 bool fitBeam(TkFittedLasBeam &beam,
0138 AlgebraicSymMatrix &covMatrix,
0139 unsigned int &hitsAtTecPlus,
0140 unsigned int &nFitParams,
0141 double &offset,
0142 double &slope,
0143 vector<GlobalPoint> &globPtrack,
0144 double &bsAngleParam,
0145 double &chi2);
0146
0147
0148
0149
0150 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken_;
0151 const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomToken_;
0152
0153
0154 edm::Handle<TkLasBeamCollection> laserBeams;
0155 edm::ESHandle<MagneticField> fieldHandle;
0156 edm::ESHandle<TrackerGeometry> geometry;
0157
0158 const edm::InputTag src_;
0159 bool fitBeamSplitters_;
0160 unsigned int nAtParameters_;
0161
0162 edm::Service<TFileService> fs;
0163
0164
0165 static vector<double> gHitZprime;
0166 static vector<double> gBarrelModuleRadius;
0167 static vector<double> gBarrelModuleOffset;
0168 static float gTIBparam;
0169 static float gTOBparam;
0170 static double gBeamR;
0171 static double gBeamZ0;
0172 static double gBeamSplitterZprime;
0173 static unsigned int gHitsAtTecMinus;
0174 static double gBSparam;
0175 static bool gFitBeamSplitters;
0176 static bool gIsInnerBarrel;
0177
0178
0179 TH1F *h_bsAngle, *h_hitX, *h_hitXTecPlus, *h_hitXTecMinus, *h_hitXAt, *h_chi2, *h_chi2ndof, *h_pull, *h_res,
0180 *h_resTecPlus, *h_resTecMinus, *h_resAt;
0181 TH2F *h_bsAngleVsBeam, *h_hitXvsZTecPlus, *h_hitXvsZTecMinus, *h_hitXvsZAt, *h_resVsZTecPlus, *h_resVsZTecMinus,
0182 *h_resVsZAt, *h_resVsHitTecPlus, *h_resVsHitTecMinus, *h_resVsHitAt;
0183 };
0184
0185
0186
0187
0188
0189
0190
0191
0192
0193
0194 vector<double> TkLasBeamFitter::gHitZprime;
0195 vector<double> TkLasBeamFitter::gBarrelModuleRadius;
0196 vector<double> TkLasBeamFitter::gBarrelModuleOffset;
0197 float TkLasBeamFitter::gTIBparam = 0.097614;
0198 float TkLasBeamFitter::gTOBparam = 0.034949;
0199 double TkLasBeamFitter::gBeamR = 0.0;
0200 double TkLasBeamFitter::gBeamZ0 = 0.0;
0201 double TkLasBeamFitter::gBeamSplitterZprime = 0.0;
0202 unsigned int TkLasBeamFitter::gHitsAtTecMinus = 0;
0203 double TkLasBeamFitter::gBSparam = 0.0;
0204 bool TkLasBeamFitter::gFitBeamSplitters = false;
0205 bool TkLasBeamFitter::gIsInnerBarrel = false;
0206
0207
0208
0209
0210 TkLasBeamFitter::TkLasBeamFitter(const edm::ParameterSet &iConfig)
0211 : magFieldToken_(esConsumes<edm::Transition::EndRun>()),
0212 geomToken_(esConsumes<edm::Transition::EndRun>()),
0213 src_(iConfig.getParameter<edm::InputTag>("src")),
0214 fitBeamSplitters_(iConfig.getParameter<bool>("fitBeamSplitters")),
0215 nAtParameters_(iConfig.getParameter<unsigned int>("numberOfFittedAtParameters")),
0216 h_bsAngle(nullptr),
0217 h_hitX(nullptr),
0218 h_hitXTecPlus(nullptr),
0219 h_hitXTecMinus(nullptr),
0220 h_hitXAt(nullptr),
0221 h_chi2(nullptr),
0222 h_chi2ndof(nullptr),
0223 h_pull(nullptr),
0224 h_res(nullptr),
0225 h_resTecPlus(nullptr),
0226 h_resTecMinus(nullptr),
0227 h_resAt(nullptr),
0228 h_bsAngleVsBeam(nullptr),
0229 h_hitXvsZTecPlus(nullptr),
0230 h_hitXvsZTecMinus(nullptr),
0231 h_hitXvsZAt(nullptr),
0232 h_resVsZTecPlus(nullptr),
0233 h_resVsZTecMinus(nullptr),
0234 h_resVsZAt(nullptr),
0235 h_resVsHitTecPlus(nullptr),
0236 h_resVsHitTecMinus(nullptr),
0237 h_resVsHitAt(nullptr) {
0238
0239 this->produces<TkFittedLasBeamCollection, edm::Transition::EndRun>();
0240 this->produces<TsosVectorCollection, edm::Transition::EndRun>();
0241
0242
0243 }
0244
0245
0246 TkLasBeamFitter::~TkLasBeamFitter() {
0247
0248
0249 }
0250
0251
0252
0253
0254
0255
0256
0257 void TkLasBeamFitter::produce(edm::Event &iEvent, const edm::EventSetup &setup) {
0258
0259 }
0260
0261
0262
0263 void TkLasBeamFitter::endRunProduce(edm::Run &run, const edm::EventSetup &setup) {
0264
0265
0266
0267
0268
0269
0270
0271
0272
0273
0274 h_hitX = fs->make<TH1F>("hitX", "local x of LAS hits;local x [cm];N", 100, -0.5, 0.5);
0275 h_hitXTecPlus = fs->make<TH1F>("hitXTecPlus", "local x of LAS hits in TECplus;local x [cm];N", 100, -0.5, 0.5);
0276 h_hitXTecMinus = fs->make<TH1F>("hitXTecMinus", "local x of LAS hits in TECminus;local x [cm];N", 100, -0.5, 0.5);
0277 h_hitXAt = fs->make<TH1F>("hitXAt", "local x of LAS hits in ATs;local x [cm];N", 100, -2.5, 2.5);
0278 h_hitXvsZTecPlus =
0279 fs->make<TH2F>("hitXvsZTecPlus", "local x vs z in TECplus;z [cm];local x [cm]", 80, 120, 280, 100, -0.5, 0.5);
0280 h_hitXvsZTecMinus =
0281 fs->make<TH2F>("hitXvsZTecMinus", "local x vs z in TECMinus;z [cm];local x [cm]", 80, -280, -120, 100, -0.5, 0.5);
0282 h_hitXvsZAt = fs->make<TH2F>("hitXvsZAt", "local x vs z in ATs;z [cm];local x [cm]", 200, -200, 200, 100, -0.5, 0.5);
0283 h_chi2 = fs->make<TH1F>("chi2", "#chi^{2};#chi^{2};N", 100, 0, 2000);
0284 h_chi2ndof = fs->make<TH1F>("chi2ndof", "#chi^{2} per degree of freedom;#chi^{2}/N_{dof};N", 100, 0, 300);
0285 h_pull = fs->make<TH1F>("pull", "pulls of #phi residuals;pull;N", 50, -10, 10);
0286 h_res = fs->make<TH1F>("res", "#phi residuals;#phi_{track} - #phi_{hit} [rad];N", 60, -0.0015, 0.0015);
0287 h_resTecPlus =
0288 fs->make<TH1F>("resTecPlus", "#phi residuals in TECplus;#phi_{track} - #phi_{hit} [rad];N", 30, -0.0015, 0.0015);
0289 h_resTecMinus = fs->make<TH1F>(
0290 "resTecMinus", "#phi residuals in TECminus;#phi_{track} - #phi_{hit} [rad];N", 60, -0.0015, 0.0015);
0291 h_resAt = fs->make<TH1F>("resAt", "#phi residuals in ATs;#phi_{track} - #phi_{hit} [rad];N", 30, -0.0015, 0.0015);
0292 h_resVsZTecPlus = fs->make<TH2F>("resVsZTecPlus",
0293 "phi residuals vs. z in TECplus;z [cm];#phi_{track} - #phi_{hit} [rad]",
0294 80,
0295 120,
0296 280,
0297 100,
0298 -0.0015,
0299 0.0015);
0300 h_resVsZTecMinus = fs->make<TH2F>("resVsZTecMinus",
0301 "phi residuals vs. z in TECminus;z [cm];#phi_{track} - #phi_{hit} [rad]",
0302 80,
0303 -280,
0304 -120,
0305 100,
0306 -0.0015,
0307 0.0015);
0308 h_resVsZAt = fs->make<TH2F>(
0309 "resVsZAt", "#phi residuals vs. z in ATs;N;#phi_{track} - #phi_{hit} [rad]", 200, -200, 200, 100, -0.0015, 0.0015);
0310 h_resVsHitTecPlus = fs->make<TH2F>("resVsHitTecPlus",
0311 "#phi residuals vs. hits in TECplus;hit no.;#phi_{track} - #phi_{hit} [rad]",
0312 144,
0313 0,
0314 144,
0315 100,
0316 -0.0015,
0317 0.0015);
0318 h_resVsHitTecMinus = fs->make<TH2F>("resVsHitTecMinus",
0319 "#phi residuals vs. hits in TECminus;hit no.;#phi_{track} - #phi_{hit} [rad]",
0320 144,
0321 0,
0322 144,
0323 100,
0324 -0.0015,
0325 0.0015);
0326 h_resVsHitAt = fs->make<TH2F>("resVsHitAt",
0327 "#phi residuals vs. hits in ATs;hit no.;#phi_{track} - #phi_{hit} [rad]",
0328 176,
0329 0,
0330 176,
0331 100,
0332 -0.0015,
0333 0.0015);
0334 h_bsAngle = fs->make<TH1F>("bsAngle", "fitted beam splitter angle;BS angle [rad];N", 40, -0.004, 0.004);
0335 h_bsAngleVsBeam = fs->make<TH2F>(
0336 "bsAngleVsBeam", "fitted beam splitter angle per beam;Beam no.;BS angle [rad]", 40, 0, 300, 100, -0.004, 0.004);
0337
0338
0339
0340 auto fittedBeams = std::make_unique<TkFittedLasBeamCollection>();
0341
0342 auto tsosesVec = std::make_unique<TsosVectorCollection>();
0343
0344
0345 run.getByLabel("LaserAlignment", "tkLaserBeams", laserBeams);
0346 geometry = setup.getHandle(geomToken_);
0347 fieldHandle = setup.getHandle(magFieldToken_);
0348
0349
0350
0351
0352
0353
0354
0355
0356
0357
0358
0359 double bsParams[40] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0360 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
0361
0362
0363 unsigned int beamNo(0);
0364
0365 gFitBeamSplitters = fitBeamSplitters_;
0366 if (fitBeamSplitters_)
0367 cout << "Fitting BS!" << endl;
0368 else
0369 cout << "BS fixed, not fitted!" << endl;
0370
0371
0372 for (TkLasBeamCollection::const_iterator iBeam = laserBeams->begin(), iEnd = laserBeams->end(); iBeam != iEnd;
0373 ++iBeam) {
0374 TkFittedLasBeam beam(*iBeam);
0375 vector<TrajectoryStateOnSurface> tsosLas;
0376
0377
0378 gBSparam = bsParams[beamNo];
0379
0380
0381 this->getLasBeams(beam, tsosLas);
0382
0383
0384 fittedBeams->push_back(beam);
0385 tsosesVec->push_back(tsosLas);
0386
0387
0388
0389
0390
0391
0392
0393
0394 beamNo++;
0395 }
0396
0397
0398 run.put(std::move(fittedBeams));
0399 run.put(std::move(tsosesVec));
0400 }
0401
0402
0403
0404
0405 void TkLasBeamFitter::getLasBeams(TkFittedLasBeam &beam, vector<TrajectoryStateOnSurface> &tsosLas) {
0406 cout << "---------------------------------------" << endl;
0407 cout << "beam id: " << beam.getBeamId()
0408 << " isTec+: " << (beam.isTecInternal(1) ? "Y" : "N") << " isTec-: " << (beam.isTecInternal(-1) ? "Y" : "N")
0409 << " isAt: " << (beam.isAlignmentTube() ? "Y" : "N") << " isR6: " << (beam.isRing6() ? "Y" : "N") << endl;
0410
0411
0412 gHitsAtTecMinus = 0;
0413 gHitZprime.clear();
0414 gBarrelModuleRadius.clear();
0415 gBarrelModuleOffset.clear();
0416
0417
0418 gBeamR = beam.isRing6() ? 84.0 : 56.4;
0419
0420 vector<const GeomDetUnit *> gd;
0421 vector<GlobalPoint> globHit;
0422 unsigned int hitsAtTecPlus(0);
0423 double sumZ(0.);
0424
0425
0426 for (TkLasBeam::const_iterator iHit = beam.begin(); iHit < beam.end(); ++iHit) {
0427
0428
0429 const SiStripLaserRecHit2D &hit(*iHit);
0430
0431 this->getLasHits(beam, hit, gd, globHit, hitsAtTecPlus);
0432 sumZ += globHit.back().z();
0433
0434
0435 h_hitX->Fill(hit.localPosition().x());
0436
0437 if (beam.isTecInternal(1)) {
0438 h_hitXTecPlus->Fill(hit.localPosition().x());
0439 h_hitXvsZTecPlus->Fill(globHit.back().z(), hit.localPosition().x());
0440 }
0441
0442 else if (beam.isTecInternal(-1)) {
0443 h_hitXTecMinus->Fill(hit.localPosition().x());
0444 h_hitXvsZTecMinus->Fill(globHit.back().z(), hit.localPosition().x());
0445 }
0446
0447 else {
0448 h_hitXAt->Fill(hit.localPosition().x());
0449 h_hitXvsZAt->Fill(globHit.back().z(), hit.localPosition().x());
0450 }
0451 }
0452
0453 gBeamZ0 = sumZ / globHit.size();
0454 double zMin(0.), zMax(0.);
0455
0456 if (beam.isTecInternal(1)) {
0457 gBeamSplitterZprime = 205.75 - gBeamZ0;
0458 zMin = 120.0 - gBeamZ0;
0459 zMax = 280.0 - gBeamZ0;
0460 }
0461
0462 else if (beam.isTecInternal(-1)) {
0463 gBeamSplitterZprime = -205.75 - gBeamZ0;
0464 zMin = -280.0 - gBeamZ0;
0465 zMax = -120.0 - gBeamZ0;
0466 }
0467
0468 else {
0469 gBeamSplitterZprime = 112.3 - gBeamZ0;
0470 zMin = -200.0 - gBeamZ0;
0471 zMax = 200.0 - gBeamZ0;
0472 }
0473
0474
0475 vector<double> hitPhi, hitPhiError, hitZprimeError;
0476
0477 for (unsigned int hit = 0; hit < globHit.size(); ++hit) {
0478 hitPhi.push_back(static_cast<double>(globHit[hit].phi()));
0479
0480 hitPhiError.push_back(0.003 / globHit[hit].perp());
0481
0482 hitZprimeError.push_back(0.0);
0483
0484 if (beam.isAlignmentTube() && abs(globHit[hit].z()) < 112.3) {
0485 gBarrelModuleRadius.push_back(globHit[hit].perp());
0486 gBarrelModuleOffset.push_back(gBarrelModuleRadius.back() - gBeamR);
0487
0488 if (gBarrelModuleOffset.back() < 0.0) {
0489 gIsInnerBarrel = true;
0490 } else {
0491 gIsInnerBarrel = false;
0492 }
0493 gHitZprime.push_back(globHit[hit].z() - gBeamZ0 - abs(gBarrelModuleOffset.back()));
0494 }
0495
0496 else {
0497 gHitZprime.push_back(globHit[hit].z() - gBeamZ0);
0498 }
0499 }
0500
0501
0502 unsigned int tecParams(3), atParams(0);
0503 if (nAtParameters_ == 3)
0504 atParams = 3;
0505 else if (nAtParameters_ == 5)
0506 atParams = 5;
0507 else
0508 atParams = 6;
0509 unsigned int nFitParams(0);
0510 if (!fitBeamSplitters_ || (hitsAtTecPlus == 0 && beam.isAlignmentTube())) {
0511 tecParams = tecParams - 1;
0512 atParams = atParams - 1;
0513 }
0514 if (beam.isTecInternal()) {
0515 nFitParams = tecParams;
0516 } else {
0517 nFitParams = atParams;
0518 }
0519
0520
0521 double offset(0.), offsetError(0.), slope(0.), slopeError(0.), bsAngleParam(0.), phiAtMinusParam(0.),
0522 phiAtPlusParam(0.), atThetaSplitParam(0.);
0523 AlgebraicSymMatrix covMatrix;
0524 if (!fitBeamSplitters_ || (beam.isAlignmentTube() && hitsAtTecPlus == 0)) {
0525 covMatrix = AlgebraicSymMatrix(nFitParams, 1);
0526 } else {
0527 covMatrix = AlgebraicSymMatrix(nFitParams - 1, 1);
0528 }
0529
0530 this->fitter(beam,
0531 covMatrix,
0532 hitsAtTecPlus,
0533 nFitParams,
0534 hitPhi,
0535 hitPhiError,
0536 hitZprimeError,
0537 zMin,
0538 zMax,
0539 bsAngleParam,
0540 offset,
0541 offsetError,
0542 slope,
0543 slopeError,
0544 phiAtMinusParam,
0545 phiAtPlusParam,
0546 atThetaSplitParam);
0547
0548 vector<GlobalPoint> globPtrack;
0549 GlobalPoint globPref;
0550 double chi2(0.);
0551
0552 for (unsigned int hit = 0; hit < gHitZprime.size(); ++hit) {
0553
0554 double trackPhi(0.), trackPhiRef(0.);
0555
0556 this->trackPhi(beam,
0557 hit,
0558 trackPhi,
0559 trackPhiRef,
0560 offset,
0561 slope,
0562 bsAngleParam,
0563 phiAtMinusParam,
0564 phiAtPlusParam,
0565 atThetaSplitParam,
0566 globHit);
0567
0568 cout << "track phi = " << trackPhi << ", hit phi = " << hitPhi[hit] << ", zPrime = " << gHitZprime[hit]
0569 << " r = " << globHit[hit].perp() << endl;
0570
0571 this->globalTrackPoint(beam, hit, hitsAtTecPlus, trackPhi, trackPhiRef, globHit, globPtrack, globPref, hitPhiError);
0572
0573
0574 const double phiResidual = globPtrack[hit].phi() - globHit[hit].phi();
0575
0576 const double phiResidualPull = phiResidual / hitPhiError[hit];
0577
0578
0579
0580
0581 chi2 += phiResidual * phiResidual / (hitPhiError[hit] * hitPhiError[hit]);
0582
0583
0584 h_res->Fill(phiResidual);
0585
0586 if (beam.isTecInternal(1)) {
0587 h_pull->Fill(phiResidualPull);
0588 h_resTecPlus->Fill(phiResidual);
0589 h_resVsZTecPlus->Fill(globPtrack[hit].z(), phiResidual);
0590
0591 if (beam.isRing6()) {
0592 h_resVsHitTecPlus->Fill(hit + (beam.getBeamId() - 1) / 10 * 9 + 72, phiResidual);
0593 }
0594
0595 else {
0596 h_resVsHitTecPlus->Fill(hit + beam.getBeamId() / 10 * 9, phiResidual);
0597 }
0598 }
0599
0600 else if (beam.isTecInternal(-1)) {
0601 h_pull->Fill(phiResidualPull);
0602 h_resTecMinus->Fill(phiResidual);
0603 h_resVsZTecMinus->Fill(globPtrack[hit].z(), phiResidual);
0604
0605 if (beam.isRing6()) {
0606 h_resVsHitTecMinus->Fill(hit + (beam.getBeamId() - 101) / 10 * 9 + 72, phiResidual);
0607 }
0608
0609 else {
0610 h_resVsHitTecMinus->Fill(hit + (beam.getBeamId() - 100) / 10 * 9, phiResidual);
0611 }
0612 }
0613
0614 else {
0615 h_pull->Fill(phiResidualPull);
0616 h_resAt->Fill(phiResidual);
0617 h_resVsZAt->Fill(globPtrack[hit].z(), phiResidual);
0618 h_resVsHitAt->Fill(hit + (beam.getBeamId() - 200) / 10 * 22, phiResidual);
0619 }
0620
0621 this->buildTrajectory(beam, hit, gd, globPtrack, tsosLas, globPref);
0622 }
0623
0624 cout << "chi^2 = " << chi2 << ", chi^2/ndof = " << chi2 / (gHitZprime.size() - nFitParams) << endl;
0625 this->fitBeam(beam, covMatrix, hitsAtTecPlus, nFitParams, offset, slope, globPtrack, bsAngleParam, chi2);
0626
0627 cout << "bsAngleParam = " << bsAngleParam << endl;
0628
0629
0630
0631 h_chi2->Fill(chi2);
0632 h_chi2ndof->Fill(chi2 / (gHitZprime.size() - nFitParams));
0633 if (bsAngleParam != 0.0) {
0634 h_bsAngle->Fill(2.0 * atan(0.5 * bsAngleParam));
0635 h_bsAngleVsBeam->Fill(beam.getBeamId(), 2.0 * atan(0.5 * bsAngleParam));
0636 }
0637 }
0638
0639
0640 void TkLasBeamFitter::getLasHits(TkFittedLasBeam &beam,
0641 const SiStripLaserRecHit2D &hit,
0642 vector<const GeomDetUnit *> &gd,
0643 vector<GlobalPoint> &globHit,
0644 unsigned int &hitsAtTecPlus) {
0645
0646 gd.push_back(geometry->idToDetUnit(hit.getDetId()));
0647 GlobalPoint globPtemp(gd.back()->toGlobal(hit.localPosition()));
0648
0649
0650 globHit.push_back(globPtemp);
0651
0652 if (beam.isAlignmentTube()) {
0653 if (abs(globPtemp.z()) > 112.3) {
0654 if (globPtemp.z() < 112.3)
0655 gHitsAtTecMinus++;
0656 else
0657 hitsAtTecPlus++;
0658 }
0659 }
0660 }
0661
0662
0663 double TkLasBeamFitter::tecPlusFunction(double *x, double *par) {
0664 double z = x[0];
0665
0666 if (z < gBeamSplitterZprime) {
0667 return par[0] + par[1] * z;
0668 } else {
0669 if (gFitBeamSplitters) {
0670
0671 return par[0] + par[1] * z - par[2] * (z - gBeamSplitterZprime) / gBeamR;
0672 } else {
0673 return par[0] + par[1] * z - gBSparam * (z - gBeamSplitterZprime) / gBeamR;
0674 }
0675 }
0676 }
0677
0678 double TkLasBeamFitter::tecMinusFunction(double *x, double *par) {
0679 double z = x[0];
0680
0681 if (z > gBeamSplitterZprime) {
0682 return par[0] + par[1] * z;
0683 } else {
0684 if (gFitBeamSplitters) {
0685
0686 return par[0] + par[1] * z + par[2] * (z - gBeamSplitterZprime) / gBeamR;
0687 } else {
0688 return par[0] + par[1] * z + gBSparam * (z - gBeamSplitterZprime) / gBeamR;
0689 }
0690 }
0691 }
0692
0693 double TkLasBeamFitter::atFunction(double *x, double *par) {
0694 double z = x[0];
0695
0696 if (z < -gBeamSplitterZprime - 2.0 * gBeamZ0) {
0697 return par[0] + par[1] * z;
0698 }
0699
0700 else if (-gBeamSplitterZprime - 2.0 * gBeamZ0 < z && z < -gBeamZ0) {
0701
0702
0703 if (!gIsInnerBarrel) {
0704 return par[0] + par[1] * z + gTOBparam * (par[2] + par[4]);
0705 }
0706
0707 else {
0708 return par[0] + par[1] * z - gTIBparam * (par[2] - par[4]);
0709 }
0710 }
0711
0712 else if (-gBeamZ0 < z && z < gBeamSplitterZprime) {
0713
0714
0715 if (!gIsInnerBarrel) {
0716 return par[0] + par[1] * z + gTOBparam * (par[3] - par[4]);
0717 }
0718
0719 else {
0720 return par[0] + par[1] * z - gTIBparam * (par[3] + par[4]);
0721 }
0722 }
0723
0724 else {
0725 if (gFitBeamSplitters) {
0726
0727 return par[0] + par[1] * z - par[5] * (z - gBeamSplitterZprime) / gBeamR;
0728 } else {
0729 return par[0] + par[1] * z - gBSparam * (z - gBeamSplitterZprime) / gBeamR;
0730 }
0731 }
0732 }
0733
0734
0735 void TkLasBeamFitter::fitter(TkFittedLasBeam &beam,
0736 AlgebraicSymMatrix &covMatrix,
0737 unsigned int &hitsAtTecPlus,
0738 unsigned int &nFitParams,
0739 vector<double> &hitPhi,
0740 vector<double> &hitPhiError,
0741 vector<double> &hitZprimeError,
0742 double &zMin,
0743 double &zMax,
0744 double &bsAngleParam,
0745 double &offset,
0746 double &offsetError,
0747 double &slope,
0748 double &slopeError,
0749 double &phiAtMinusParam,
0750 double &phiAtPlusParam,
0751 double &atThetaSplitParam) {
0752 TGraphErrors *lasData =
0753 new TGraphErrors(gHitZprime.size(), &(gHitZprime[0]), &(hitPhi[0]), &(hitZprimeError[0]), &(hitPhiError[0]));
0754
0755
0756 if (beam.isTecInternal(1)) {
0757 TF1 tecPlus("tecPlus", tecPlusFunction, zMin, zMax, nFitParams);
0758 tecPlus.SetParameter(1, 0);
0759 tecPlus.SetParameter(nFitParams - 1, 0);
0760 lasData->Fit(&tecPlus, "R");
0761 } else if (beam.isTecInternal(-1)) {
0762 TF1 tecMinus("tecMinus", tecMinusFunction, zMin, zMax, nFitParams);
0763 tecMinus.SetParameter(1, 0);
0764 tecMinus.SetParameter(nFitParams - 1, 0);
0765 lasData->Fit(&tecMinus, "R");
0766 } else {
0767 TF1 at("at", atFunction, zMin, zMax, nFitParams);
0768 at.SetParameter(1, 0);
0769 at.SetParameter(nFitParams - 1, 0);
0770 lasData->Fit(&at, "R");
0771 }
0772
0773
0774 gMinuit->GetParameter(0, offset, offsetError);
0775 gMinuit->GetParameter(1, slope, slopeError);
0776
0777
0778
0779 double bsAngleParamError(0.), phiAtMinusParamError(0.), phiAtPlusParamError(0.), atThetaSplitParamError(0.);
0780
0781 if (beam.isAlignmentTube()) {
0782 gMinuit->GetParameter(2, phiAtMinusParam, phiAtMinusParamError);
0783 gMinuit->GetParameter(3, phiAtPlusParam, phiAtPlusParamError);
0784 gMinuit->GetParameter(4, atThetaSplitParam, atThetaSplitParamError);
0785 }
0786
0787 if (fitBeamSplitters_) {
0788 if (beam.isAlignmentTube() && hitsAtTecPlus == 0) {
0789 bsAngleParam = gBSparam;
0790 } else {
0791 gMinuit->GetParameter(nFitParams - 1, bsAngleParam, bsAngleParamError);
0792 }
0793 } else {
0794 bsAngleParam = gBSparam;
0795 }
0796
0797
0798 vector<double> vec(covMatrix.num_col() * covMatrix.num_col());
0799 gMinuit->mnemat(&vec[0], covMatrix.num_col());
0800 for (int col = 0; col < covMatrix.num_col(); col++) {
0801 for (int row = 0; row < covMatrix.num_col(); row++) {
0802 covMatrix[col][row] = vec[row + covMatrix.num_col() * col];
0803 }
0804 }
0805
0806
0807
0808 delete lasData;
0809 }
0810
0811
0812 void TkLasBeamFitter::trackPhi(TkFittedLasBeam &beam,
0813 unsigned int &hit,
0814 double &trackPhi,
0815 double &trackPhiRef,
0816 double &offset,
0817 double &slope,
0818 double &bsAngleParam,
0819 double &phiAtMinusParam,
0820 double &phiAtPlusParam,
0821 double &atThetaSplitParam,
0822 vector<GlobalPoint> &globHit) {
0823
0824 if (beam.isTecInternal(1)) {
0825 if (gHitZprime[hit] < gBeamSplitterZprime) {
0826 trackPhi = offset + slope * gHitZprime[hit];
0827 trackPhiRef = offset + slope * (gHitZprime[hit] + 1.0);
0828 } else {
0829 trackPhi = offset + slope * gHitZprime[hit] - bsAngleParam * (gHitZprime[hit] - gBeamSplitterZprime) / gBeamR;
0830 trackPhiRef = offset + slope * (gHitZprime[hit] + 1.0) -
0831 bsAngleParam * ((gHitZprime[hit] + 1.0) - gBeamSplitterZprime) / gBeamR;
0832 }
0833 }
0834
0835 else if (beam.isTecInternal(-1)) {
0836 if (gHitZprime[hit] > gBeamSplitterZprime) {
0837 trackPhi = offset + slope * gHitZprime[hit];
0838 trackPhiRef = offset + slope * (gHitZprime[hit] + 1.0);
0839 } else {
0840 trackPhi = offset + slope * gHitZprime[hit] + bsAngleParam * (gHitZprime[hit] - gBeamSplitterZprime) / gBeamR;
0841 trackPhiRef = offset + slope * (gHitZprime[hit] + 1.0) +
0842 bsAngleParam * ((gHitZprime[hit] + 1.0) - gBeamSplitterZprime) / gBeamR;
0843 }
0844 }
0845
0846 else {
0847
0848 if (gHitZprime[hit] < -gBeamSplitterZprime - 2.0 * gBeamZ0) {
0849 trackPhi = offset + slope * gHitZprime[hit];
0850 trackPhiRef = offset + slope * (gHitZprime[hit] + 1.0);
0851 }
0852
0853 else if (-gBeamSplitterZprime - 2.0 * gBeamZ0 < gHitZprime[hit] && gHitZprime[hit] < -gBeamZ0) {
0854 if (!gIsInnerBarrel) {
0855 trackPhi = offset + slope * gHitZprime[hit] + gTOBparam * (phiAtMinusParam + atThetaSplitParam);
0856 } else {
0857 trackPhi = offset + slope * gHitZprime[hit] - gTIBparam * (phiAtMinusParam - atThetaSplitParam);
0858 }
0859 trackPhiRef = offset + slope * (gHitZprime[hit] + abs(gBarrelModuleOffset[hit - gHitsAtTecMinus]));
0860 }
0861
0862 else if (-gBeamZ0 < gHitZprime[hit] && gHitZprime[hit] < gBeamSplitterZprime) {
0863 if (!gIsInnerBarrel) {
0864 trackPhi = offset + slope * gHitZprime[hit] + gTOBparam * (phiAtPlusParam - atThetaSplitParam);
0865 } else {
0866 trackPhi = offset + slope * gHitZprime[hit] - gTIBparam * (phiAtPlusParam + atThetaSplitParam);
0867 }
0868 trackPhiRef = offset + slope * (gHitZprime[hit] + abs(gBarrelModuleOffset[hit - gHitsAtTecMinus]));
0869 }
0870
0871 else {
0872 trackPhi = offset + slope * gHitZprime[hit] - bsAngleParam * (gHitZprime[hit] - gBeamSplitterZprime) / gBeamR;
0873 trackPhiRef = offset + slope * (gHitZprime[hit] + 1.0) -
0874 bsAngleParam * ((gHitZprime[hit] + 1.0) - gBeamSplitterZprime) / gBeamR;
0875 }
0876 }
0877 }
0878
0879
0880 void TkLasBeamFitter::globalTrackPoint(TkFittedLasBeam &beam,
0881 unsigned int &hit,
0882 unsigned int &hitsAtTecPlus,
0883 double &trackPhi,
0884 double &trackPhiRef,
0885 vector<GlobalPoint> &globHit,
0886 vector<GlobalPoint> &globPtrack,
0887 GlobalPoint &globPref,
0888 vector<double> &hitPhiError) {
0889
0890 if (beam.isTecInternal(0)) {
0891 globPtrack.push_back(GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhi, globHit[hit].z())));
0892 globPref = GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhiRef, globHit[hit].z() + 1.0));
0893 }
0894
0895 else {
0896
0897 if (hit < gHitsAtTecMinus) {
0898 globPtrack.push_back(GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhi, globHit[hit].z())));
0899 globPref = GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhiRef, globHit[hit].z() + 1.0));
0900 }
0901
0902 else if (hit > gHitZprime.size() - hitsAtTecPlus - 1) {
0903 globPtrack.push_back(GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhi, globHit[hit].z())));
0904 globPref = GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhiRef, globHit[hit].z() + 1.0));
0905 }
0906
0907 else {
0908 globPtrack.push_back(GlobalPoint(GlobalPoint::Cylindrical(globHit[hit].perp(), trackPhi, globHit[hit].z())));
0909 globPref = GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhiRef, globHit[hit].z()));
0910 }
0911 }
0912 }
0913
0914
0915 void TkLasBeamFitter::buildTrajectory(TkFittedLasBeam &beam,
0916 unsigned int &hit,
0917 vector<const GeomDetUnit *> &gd,
0918 vector<GlobalPoint> &globPtrack,
0919 vector<TrajectoryStateOnSurface> &tsosLas,
0920 GlobalPoint &globPref) {
0921 const MagneticField *magneticField = fieldHandle.product();
0922 GlobalVector trajectoryState;
0923
0924
0925 if (beam.isTecInternal(1)) {
0926 trajectoryState = GlobalVector(globPref - globPtrack[hit]);
0927 }
0928
0929 else if (beam.isTecInternal(-1)) {
0930 trajectoryState = GlobalVector(globPtrack[hit] - globPref);
0931 }
0932
0933 else {
0934
0935 if (gHitZprime[hit] < -gBeamSplitterZprime - 2.0 * gBeamZ0) {
0936 trajectoryState = GlobalVector(globPtrack[hit] - globPref);
0937 }
0938
0939 else if (gHitZprime[hit] > gBeamSplitterZprime) {
0940 trajectoryState = GlobalVector(globPref - globPtrack[hit]);
0941 }
0942
0943 else {
0944 trajectoryState = GlobalVector(globPtrack[hit] - globPref);
0945 }
0946 }
0947
0948 const FreeTrajectoryState ftsLas = FreeTrajectoryState(globPtrack[hit], trajectoryState, 0, magneticField);
0949 tsosLas.push_back(TrajectoryStateOnSurface(ftsLas, gd[hit]->surface(), SurfaceSideDefinition::beforeSurface));
0950 }
0951
0952
0953 bool TkLasBeamFitter::fitBeam(TkFittedLasBeam &beam,
0954 AlgebraicSymMatrix &covMatrix,
0955 unsigned int &hitsAtTecPlus,
0956 unsigned int &nFitParams,
0957 double &offset,
0958 double &slope,
0959 vector<GlobalPoint> &globPtrack,
0960 double &bsAngleParam,
0961 double &chi2) {
0962
0963 unsigned int paramType(0);
0964 if (!fitBeamSplitters_)
0965 paramType = 1;
0966 if (beam.isAlignmentTube() && hitsAtTecPlus == 0)
0967 paramType = 0;
0968
0969
0970
0971 const unsigned int nPedeParams(nFitParams);
0972
0973
0974 std::vector<TkFittedLasBeam::Scalar> params(nPedeParams);
0975 params[0] = offset;
0976 params[1] = slope;
0977
0978
0979
0980 AlgebraicMatrix derivatives(gHitZprime.size(), nPedeParams);
0981
0982 for (unsigned int hit = 0; hit < gHitZprime.size(); ++hit) {
0983
0984 derivatives[hit][0] = 1.0;
0985
0986
0987
0988 if (beam.isTecInternal(1)) {
0989 derivatives[hit][1] = globPtrack[hit].z();
0990
0991
0992
0993
0994
0995
0996 }
0997
0998 else if (beam.isTecInternal(-1)) {
0999 derivatives[hit][1] = globPtrack[hit].z();
1000
1001
1002
1003
1004
1005
1006 }
1007
1008 else {
1009
1010 if (gHitZprime[hit] < -gBeamSplitterZprime - 2.0 * gBeamZ0) {
1011 derivatives[hit][1] = globPtrack[hit].z();
1012
1013
1014
1015 }
1016
1017 else if (gHitZprime[hit] > gBeamSplitterZprime) {
1018 derivatives[hit][1] = globPtrack[hit].z();
1019
1020
1021
1022 }
1023
1024 else {
1025 derivatives[hit][1] = globPtrack[hit].z() - gBarrelModuleOffset[hit - gHitsAtTecMinus];
1026
1027
1028
1029 }
1030 }
1031 }
1032
1033 unsigned int firstFixedParam(covMatrix.num_col());
1034
1035
1036
1037
1038 beam.setParameters(paramType, params, covMatrix, derivatives, firstFixedParam, chi2);
1039
1040 return true;
1041 }
1042
1043
1044
1045 DEFINE_FWK_MODULE(TkLasBeamFitter);