File indexing completed on 2021-02-14 12:45:17
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 #include "DataFormats/GeometrySurface/interface/Surface.h"
0012 #include "Alignment/MillePedeAlignmentAlgorithm/interface/MillePedeMonitor.h"
0013
0014 #include "DataFormats/TrackReco/interface/Track.h"
0015
0016 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0017
0018 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0019 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
0020 #include <Geometry/CommonDetUnit/interface/GeomDetType.h>
0021
0022 #include "Alignment/CommonAlignment/interface/Alignable.h"
0023 #include "Alignment/CommonAlignment/interface/AlignableBeamSpot.h"
0024 #include "Alignment/CommonAlignment/interface/AlignableDetOrUnitPtr.h"
0025 #include "Alignment/CommonAlignmentParametrization/interface/FrameToFrameDerivative.h"
0026
0027 #include "DataFormats/DetId/interface/DetId.h"
0028 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0029 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0030 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0031 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0032 const int kBPIX = PixelSubdetector::PixelBarrel;
0033 const int kFPIX = PixelSubdetector::PixelEndcap;
0034
0035 #include <TProfile2D.h>
0036 #include <TFile.h>
0037 #include <TDirectory.h>
0038 #include <TMath.h>
0039
0040 #include "Alignment/MillePedeAlignmentAlgorithm/interface/PedeLabelerBase.h"
0041
0042 typedef TransientTrackingRecHit::ConstRecHitPointer ConstRecHitPointer;
0043
0044
0045 MillePedeMonitor::MillePedeMonitor(const TrackerTopology *tTopo, const char *rootFileName)
0046 : myRootDir(nullptr), myDeleteDir(false), trackerTopology(tTopo) {
0047 myRootDir = TFile::Open(rootFileName, "recreate");
0048 myDeleteDir = true;
0049
0050 this->init(myRootDir);
0051 }
0052
0053
0054 MillePedeMonitor::MillePedeMonitor(TDirectory *rootDir, const TrackerTopology *tTopo)
0055 : myRootDir(nullptr), myDeleteDir(false), trackerTopology(tTopo) {
0056
0057
0058 myRootDir = rootDir;
0059 myDeleteDir = false;
0060
0061 this->init(myRootDir);
0062 }
0063
0064
0065 MillePedeMonitor::~MillePedeMonitor() {
0066 myRootDir->Write();
0067 if (myDeleteDir)
0068 delete myRootDir;
0069 }
0070
0071
0072 bool MillePedeMonitor::init(TDirectory *directory) {
0073 if (!directory)
0074 return false;
0075 TDirectory *oldDir = gDirectory;
0076
0077 const int kNumBins = 20;
0078 double binsPt[kNumBins + 1] = {0.};
0079 if (!this->equidistLogBins(binsPt, kNumBins, 0.8, 100.)) {
0080
0081 }
0082 const int nHits = 35;
0083
0084 myTrackHists1D.push_back(new TH1F("ptTrackLogBins", "p_{t}(track);p_{t} [GeV]", kNumBins, binsPt));
0085
0086 myTrackHists1D.push_back(new TH1F("ptTrack", "p_{t}(track);p_{t} [GeV]", kNumBins, binsPt[0], binsPt[kNumBins]));
0087 myTrackHists1D.push_back(new TH1F("pTrack", "p(track);p [GeV]", kNumBins, binsPt[0], 1.3 * binsPt[kNumBins]));
0088 myTrackHists1D.push_back(new TH1F("etaTrack", "#eta(track);#eta", 26, -2.6, 2.6));
0089 myTrackHists1D.push_back(new TH1F("thetaTrack", "#theta(track);#theta", 100, 0., TMath::Pi()));
0090 myTrackHists1D.push_back(new TH1F("phiTrack", "#phi(track);#phi", 15, -TMath::Pi(), TMath::Pi()));
0091
0092 myTrackHists1D.push_back(new TH1F("nHitTrack", "N_{hit}(track);N_{hit}", 40, 5., 45.));
0093 myTrackHists1D.push_back(new TH1F("nHitBPIXTrack", "N_{hit, BPIX}(track);N_{hit}", 45, 0., 45.));
0094 myTrackHists1D.push_back(new TH1F("nHitFPIXTrack", "N_{hit, FPIX}(track);N_{hit}", 45, 0., 45.));
0095 myTrackHists1D.push_back(new TH1F("nHitFPIXplusTrack", "N_{hit, BPIXplus}(track);N_{hit}", 45, 0., 45.));
0096 myTrackHists1D.push_back(new TH1F("nHitFPIXminusTrack", "N_{hit, BPIXminus}(track);N_{hit}", 45, 0., 45.));
0097 myTrackHists1D.push_back(new TH1F("nHitPIXELTrack", "N_{hit, PIXEL}(track);N_{hit}", 45, 0., 45.));
0098 myTrackHists1D.push_back(new TH1F("nHitTIBTrack", "N_{hit, TIB}(track);N_{hit}", 45, 0., 45.));
0099 myTrackHists1D.push_back(new TH1F("nHitTOBTrack", "N_{hit, TOB}(track);N_{hit}", 45, 0., 45.));
0100 myTrackHists1D.push_back(new TH1F("nHitTIDplusTrack", "N_{hit, TIDplus}(track);N_{hit}", 45, 0., 45.));
0101 myTrackHists1D.push_back(new TH1F("nHitTIDminusTrack", "N_{hit, TIDminus}(track);N_{hit}", 45, 0., 45.));
0102 myTrackHists1D.push_back(new TH1F("nHitTIDTrack", "N_{hit, TID}(track);N_{hit}", 45, 0., 45.));
0103 myTrackHists1D.push_back(new TH1F("nHitTECplusTrack", "N_{hit, TECplus}(track);N_{hit}", 45, 0., 45.));
0104 myTrackHists1D.push_back(new TH1F("nHitTECminusTrack", "N_{hit, TECminus}(track);N_{hit}", 45, 0., 45.));
0105 myTrackHists1D.push_back(new TH1F("nHitTECTrack", "N_{hit, TEC}(track);N_{hit}", 45, 0., 45.));
0106 myTrackHists1D.push_back(new TH1F("nHitENDCAPplusTrack", "N_{hit, ENDCAPplus}(track);N_{hit}", 45, 0., 45.));
0107 myTrackHists1D.push_back(new TH1F("nHitENDCAPminusTrack", "N_{hit, ENDCAPminus}(track);N_{hit}", 45, 0., 45.));
0108 myTrackHists1D.push_back(new TH1F("nHitENDCAPTrack", "N_{hit, ENDCAP}(track);N_{hit}", 45, 0., 45.));
0109 myTrackHists2D.push_back(new TH2F(
0110 "nHitENDCAPTrackMinusVsPlus", "N_{hit, ENDCAP}(track);N_{hit, plus};N_{hit, minus}", 45, 0., 45., 45, 0., 45.));
0111 myTrackHists1D.push_back(new TH1F("nHitInvalidTrack", "N_{hit, invalid}(track)", 5, 0., 5.));
0112 myTrackHists1D.push_back(new TH1F("r1Track", "r(1st hit);r [cm]", 20, 0., 20.));
0113 myTrackHists1D.push_back(new TH1F("phi1Track", "#phi(1st hit);#phi", 30, -TMath::Pi(), TMath::Pi()));
0114 myTrackHists1D.push_back(new TH1F("y1Track", "y(1st hit);y [cm]", 40, -120., +120.));
0115 myTrackHists1D.push_back(new TH1F("z1Track", "z(1st hit);z [cm]", 20, -50, +50));
0116 myTrackHists1D.push_back(new TH1F("r1TrackSigned", "r(1st hit);r_{#pm} [cm]", 40, -120., 120.));
0117 myTrackHists1D.push_back(new TH1F("z1TrackFull", "z(1st hit);z [cm]", 20, -300., +300.));
0118 myTrackHists1D.push_back(new TH1F("rLastTrack", "r(last hit);r [cm]", 20, 20., 120.));
0119 myTrackHists1D.push_back(new TH1F("phiLastTrack", "#phi(last hit);#phi", 30, -TMath::Pi(), TMath::Pi()));
0120 myTrackHists1D.push_back(new TH1F("yLastTrack", "y(last hit);y [cm]", 40, -120., +120.));
0121 myTrackHists1D.push_back(new TH1F("zLastTrack", "z(last hit);z [cm]", 30, -300., +300.));
0122 myTrackHists1D.push_back(new TH1F("chi2PerNdf", "#chi^{2}/ndf;#chi^{2}/ndf", 500, 0., 50.));
0123 myTrackHists1D.push_back(new TH1F("impParZ", "impact parameter in z", 20, -20., 20.));
0124 myTrackHists1D.push_back(new TH1F("impParErrZ", "error of impact parameter in z", 40, 0., 0.06));
0125 myTrackHists1D.push_back(new TH1F("impParRphi", "impact parameter in r#phi", 51, -0.05, .05));
0126 myTrackHists1D.push_back(new TH1F("impParErrRphi", "error of impact parameter in r#phi", 50, 0., 0.01));
0127
0128 myTrackHists2D.push_back(
0129 new TH2F("rz1TrackFull", "rz(1st hit);z [cm]; r_{#pm} [cm]", 40, -282., +282., 40, -115., +115.));
0130 myTrackHists2D.push_back(new TH2F("xy1Track", "xy(1st hit);x [cm]; y [cm]", 40, -115., +115., 40, -115., +115.));
0131
0132 TDirectory *dirTracks = directory->mkdir("trackHists", "input tracks");
0133 this->addToDirectory(myTrackHists1D, dirTracks);
0134 this->addToDirectory(myTrackHists2D, dirTracks);
0135
0136
0137 myUsedTrackHists1D = this->cloneHists(myTrackHists1D, "used", " (used tracks)");
0138 myUsedTrackHists2D = this->cloneHists(myTrackHists2D, "used", " (used tracks)");
0139
0140 myUsedTrackHists1D.push_back(new TH1F("usedHitsX", "n(x-hits) transferred to pede;n(x-hits)", nHits, 0., nHits));
0141 myUsedTrackHists1D.push_back(new TH1F("usedHitsY", "n(y-hits) transferred to pede;n(y-hits)", 10, 0., 10));
0142
0143 TDirectory *dirUsedTracks = directory->mkdir("usedTrackHists", "used tracks");
0144 this->addToDirectory(myUsedTrackHists1D, dirUsedTracks);
0145 this->addToDirectory(myUsedTrackHists2D, dirUsedTracks);
0146
0147
0148 myTrajectoryHists1D.push_back(new TH1F("validRefTraj", "validity of ReferenceTrajectory", 2, 0., 2.));
0149
0150 myTrajectoryHists2D.push_back(new TProfile2D(
0151 "profCorr", "mean of |#rho|, #rho#neq0;hit x/y;hit x/y;", 2 * nHits, 0., nHits, 2 * nHits, 0., nHits));
0152 myTrajectoryHists2D.push_back(new TProfile2D("profCorrOffXy",
0153 "mean of |#rho|, #rho#neq0, no xy_{hit};hit x/y;hit x/y;",
0154 2 * nHits,
0155 0.,
0156 nHits,
0157 2 * nHits,
0158 0.,
0159 nHits));
0160
0161 myTrajectoryHists2D.push_back(new TH2F(
0162 "hitCorrOffBlock", "hit correlations: off-block-diagonals;N(hit);#rho", 2 * nHits, 0., nHits, 81, -.06, .06));
0163 TArrayD logBins(102);
0164 this->equidistLogBins(logBins.GetArray(), logBins.GetSize() - 1, 1.E-11, .1);
0165 myTrajectoryHists2D.push_back(new TH2F("hitCorrOffBlockLog",
0166 "hit correlations: off-block-diagonals;N(hit);|#rho|",
0167 2 * nHits,
0168 0.,
0169 nHits,
0170 logBins.GetSize() - 1,
0171 logBins.GetArray()));
0172
0173 myTrajectoryHists2D.push_back(
0174 new TH2F("hitCorrXy", "hit correlations: xy;N(hit);#rho", nHits, 0., nHits, 81, -.5, .5));
0175 myTrajectoryHists2D.push_back(
0176 new TH2F("hitCorrXyValid", "hit correlations: xy, 2D-det.;N(hit);#rho", nHits, 0., nHits, 81, -.02, .02));
0177 this->equidistLogBins(logBins.GetArray(), logBins.GetSize() - 1, 1.E-10, .5);
0178 myTrajectoryHists2D.push_back(new TH2F(
0179 "hitCorrXyLog", "hit correlations: xy;N(hit);|#rho|", nHits, 0., nHits, logBins.GetSize() - 1, logBins.GetArray()));
0180 myTrajectoryHists2D.push_back(new TH2F("hitCorrXyLogValid",
0181 "hit correlations: xy, 2D-det.;N(hit);|#rho|",
0182 nHits,
0183 0.,
0184 nHits,
0185 logBins.GetSize() - 1,
0186 logBins.GetArray()));
0187
0188 myTrajectoryHists1D.push_back(new TH1F("measLocX", "local x measurements;x", 101, -6., 6.));
0189 myTrajectoryHists1D.push_back(new TH1F("measLocY", "local y measurements, 2D-det.;y", 101, -10., 10.));
0190 myTrajectoryHists1D.push_back(new TH1F("trajLocX", "local x trajectory;x", 101, -6., 6.));
0191 myTrajectoryHists1D.push_back(new TH1F("trajLocY", "local y trajectory, 2D-det.;y", 101, -10., 10.));
0192
0193 myTrajectoryHists1D.push_back(new TH1F("residLocX", "local x residual;#Deltax", 101, -.75, .75));
0194 myTrajectoryHists1D.push_back(new TH1F("residLocY", "local y residual, 2D-det.;#Deltay", 101, -2., 2.));
0195 myTrajectoryHists1D.push_back(new TH1F("reduResidLocX", "local x reduced residual;#Deltax/#sigma", 101, -20., 20.));
0196 myTrajectoryHists1D.push_back(
0197 new TH1F("reduResidLocY", "local y reduced residual, 2D-det.;#Deltay/#sigma", 101, -20., 20.));
0198
0199
0200 myTrajectoryHists2D.push_back(
0201 new TH2F("measLocXvsHit", "local x measurements;hit;x", nHits, 0., nHits, 101, -6., 6.));
0202 myTrajectoryHists2D.push_back(
0203 new TH2F("measLocYvsHit", "local y measurements, 2D-det.;hit;y", nHits, 0., nHits, 101, -10., 10.));
0204 myTrajectoryHists2D.push_back(new TH2F("trajLocXvsHit", "local x trajectory;hit;x", nHits, 0., nHits, 101, -6., 6.));
0205 myTrajectoryHists2D.push_back(
0206 new TH2F("trajLocYvsHit", "local y trajectory, 2D-det.;hit;y", nHits, 0., nHits, 101, -10., 10.));
0207
0208 myTrajectoryHists2D.push_back(
0209 new TH2F("residLocXvsHit", "local x residual;hit;#Deltax", nHits, 0., nHits, 101, -.75, .75));
0210 myTrajectoryHists2D.push_back(
0211 new TH2F("residLocYvsHit", "local y residual, 2D-det.;hit;#Deltay", nHits, 0., nHits, 101, -1., 1.));
0212 myTrajectoryHists2D.push_back(
0213 new TH2F("reduResidLocXvsHit", "local x reduced residual;hit;#Deltax/#sigma", nHits, 0., nHits, 101, -20., 20.));
0214 myTrajectoryHists2D.push_back(new TH2F(
0215 "reduResidLocYvsHit", "local y reduced residual, 2D-det.;hit;#Deltay/#sigma", nHits, 0., nHits, 101, -20., 20.));
0216
0217 myTrajectoryHists2D.push_back(
0218 new TProfile2D("profDerivatives", "mean derivatives;hit x/y;parameter;", 2 * nHits, 0., nHits, 10, 0., 10.));
0219
0220 myTrajectoryHists2D.push_back(new TH2F(
0221 "derivatives", "derivative;parameter;#partial(x/y)_{local}/#partial(param)", 10, 0., 10., 101, -20., 20.));
0222 this->equidistLogBins(logBins.GetArray(), logBins.GetSize() - 1, 1.E-12, 100.);
0223 myTrajectoryHists2D.push_back(new TH2F("derivativesLog",
0224 "|derivative|;parameter;|#partial(x/y)_{local}/#partial(param)|",
0225 10,
0226 0.,
0227 10.,
0228 logBins.GetSize() - 1,
0229 logBins.GetArray()));
0230 myTrajectoryHists2D.push_back(new TH2F("derivativesVsPhi",
0231 "derivatives vs. #phi;#phi(geomDet);#partial(x/y)_{local}/#partial(param)",
0232 50,
0233 -TMath::Pi(),
0234 TMath::Pi(),
0235 101,
0236 -300.,
0237 300.));
0238
0239
0240 TDirectory *dirTraject = directory->mkdir("refTrajectoryHists", "ReferenceTrajectory's");
0241 this->addToDirectory(myTrajectoryHists2D, dirTraject);
0242 this->addToDirectory(myTrajectoryHists1D, dirTraject);
0243
0244
0245 myDerivHists2D.push_back(new TH2F("localDerivsPar",
0246 "local derivatives vs. paramater;parameter;#partial/#partial(param)",
0247 6,
0248 0.,
0249 6.,
0250 101,
0251 -200.,
0252 200.));
0253 myDerivHists2D.push_back(new TH2F("localDerivsPhi",
0254 "local derivatives vs. #phi(det);#phi(det);#partial/#partial(param)",
0255 51,
0256 -TMath::Pi(),
0257 TMath::Pi(),
0258 101,
0259 -150.,
0260 150.));
0261 this->equidistLogBins(logBins.GetArray(), logBins.GetSize() - 1, 1.E-13, 150.);
0262 myDerivHists2D.push_back(new TH2F("localDerivsParLog",
0263 "local derivatives (#neq 0) vs. parameter;parameter;|#partial/#partial(param)|",
0264 6,
0265 0.,
0266 6.,
0267 logBins.GetSize() - 1,
0268 logBins.GetArray()));
0269 myDerivHists2D.push_back(new TH2F("localDerivsPhiLog",
0270 "local derivatives (#neq 0) vs. #phi(det);#phi(det);|#partial/#partial(param)|",
0271 51,
0272 -TMath::Pi(),
0273 TMath::Pi(),
0274 logBins.GetSize() - 1,
0275 logBins.GetArray()));
0276 const unsigned int maxParNum = PedeLabelerBase::theMaxNumParam;
0277 myDerivHists2D.push_back(new TH2F("globalDerivsPar",
0278 "global derivatives vs. paramater;parameter;#partial/#partial(param)",
0279 maxParNum,
0280 0.,
0281 maxParNum,
0282 100,
0283 -200,
0284 200));
0285 myDerivHists2D.push_back(new TH2F("globalDerivsPhi",
0286 "global derivatives vs. #phi(det);#phi(det);#partial/#partial(param)",
0287 51,
0288 -TMath::Pi(),
0289 TMath::Pi(),
0290 102,
0291 -1.02,
0292 1.02));
0293 this->equidistLogBins(logBins.GetArray(), logBins.GetSize() - 1, 1.E-7, 5.e2);
0294 myDerivHists2D.push_back(new TH2F("globalDerivsParLog",
0295 "global derivatives (#neq 0) vs. parameter;parameter;|#partial/#partial(param)|",
0296 maxParNum,
0297 0.,
0298 maxParNum,
0299 logBins.GetSize() - 1,
0300 logBins.GetArray()));
0301 myDerivHists2D.push_back(new TH2F("globalDerivsPhiLog",
0302 "global derivatives (#neq 0) vs. #phi(det);#phi(det);|#partial/#partial(param)|",
0303 51,
0304 -TMath::Pi(),
0305 TMath::Pi(),
0306 logBins.GetSize() - 1,
0307 logBins.GetArray()));
0308
0309
0310
0311
0312
0313
0314 TDirectory *dirDerivs = directory->mkdir("derivatives", "derivatives etc.");
0315 this->addToDirectory(myDerivHists2D, dirDerivs);
0316
0317
0318 myResidHists2D.push_back(new TH2F(
0319 "residPhi", "residuum vs. #phi(det);#phi(det);residuum[cm]", 51, -TMath::Pi(), TMath::Pi(), 101, -.5, .5));
0320 myResidHists2D.push_back(
0321 new TH2F("sigmaPhi", "#sigma vs. #phi(det);#phi(det);#sigma[cm]", 51, -TMath::Pi(), TMath::Pi(), 101, .0, .1));
0322 myResidHists2D.push_back(new TH2F("reduResidPhi",
0323 "residuum/#sigma vs. #phi(det);#phi(det);residuum/#sigma",
0324 51,
0325 -TMath::Pi(),
0326 TMath::Pi(),
0327 101,
0328 -14.,
0329 14.));
0330
0331
0332
0333
0334
0335
0336
0337
0338
0339
0340
0341
0342
0343
0344
0345
0346
0347
0348
0349
0350
0351
0352
0353
0354
0355
0356
0357
0358
0359
0360
0361
0362
0363
0364
0365
0366
0367
0368 TDirectory *dirResid = directory->mkdir("residuals", "hit residuals, sigma,...");
0369 this->addToDirectory(myResidHists2D, dirResid);
0370
0371
0372
0373
0374 std::vector<TH1 *> allResidHistsX;
0375 allResidHistsX.push_back(new TH1F("resid", "hit residuals;residuum [cm]", 101, -.5, .5));
0376
0377 allResidHistsX.push_back(new TH1F("sigma", "hit uncertainties;#sigma [cm]", 100, 0., 1.));
0378
0379 allResidHistsX.push_back(
0380 new TH1F("reduResid", "reduced hit residuals;res./#sigma", 101, -10., 10.));
0381
0382 allResidHistsX.push_back(
0383 new TH1F("angle", "#phi_{tr} wrt normal (sens. plane);#phi_{n}^{sens}", 50, 0., TMath::PiOver2()));
0384 allResidHistsX.push_back(new TH2F("residVsAngle",
0385 "residuum vs. #phi_{tr} wrt normal (sens. plane);#phi_{n}^{sens};residuum [cm]",
0386 50,
0387 0.,
0388 TMath::PiOver2(),
0389 51,
0390 -1.,
0391 1.));
0392 this->equidistLogBins(logBins.GetArray(), logBins.GetSize() - 1, 1.E-6, 100.);
0393 allResidHistsX.push_back(new TH2F("sigmaVsAngle",
0394 "#sigma vs. #phi_{tr} wrt normal (sens. plane);#phi_{n}^{sens};#sigma [cm]",
0395 50,
0396 0.,
0397 TMath::PiOver2(),
0398 logBins.GetSize() - 1,
0399 logBins.GetArray()));
0400 allResidHistsX.push_back(
0401 new TH2F("reduResidVsAngle",
0402 "reduced residuum vs. #phi_{tr} wrt normal (sens. plane);#phi_{n}^{sens};res./#sigma",
0403 50,
0404 0.,
0405 TMath::PiOver2(),
0406 51,
0407 -15.,
0408 15.));
0409
0410 allResidHistsX.push_back(
0411 new TH1F("residGt45", "hit residuals (#phi_{n}^{sens}>45#circ);residuum [cm]", 101, -.5, .5));
0412
0413 allResidHistsX.push_back(
0414 new TH1F("sigmaGt45", "hit uncertainties(#phi_{n}^{sens}>45#circ);#sigma [cm]", 100, 0., 1.));
0415
0416 allResidHistsX.push_back(new TH1F(
0417 "reduResidGt45", "reduced hit residuals(#phi_{n}^{sens}>45#circ);res./#sigma", 101, -10., 10.));
0418
0419 allResidHistsX.push_back(
0420 new TH1F("residLt45", "hit residuals (#phi_{n}^{sens}<45#circ);residuum [cm]", 101, -.5, .5));
0421
0422 allResidHistsX.push_back(
0423 new TH1F("sigmaLt45", "hit uncertainties(#phi_{n}^{sens}<45#circ);#sigma [cm]", 100, 0., 1.));
0424
0425 allResidHistsX.push_back(new TH1F(
0426 "reduResidLt45", "reduced hit residuals(#phi_{n}^{sens}<45#circ);res./#sigma", 101, -10., 10.));
0427 myResidHistsVec1DX.push_back(allResidHistsX);
0428
0429 myResidHistsVec1DX.push_back(this->cloneHists(allResidHistsX, "TPB", " x-coord. in pixel barrel"));
0430 myResidHistsVec1DX.push_back(this->cloneHists(allResidHistsX, "TPE", " x-coord. in pixel discs"));
0431 myResidHistsVec1DX.push_back(this->cloneHists(allResidHistsX, "TIB", " x-coord. in TIB"));
0432 myResidHistsVec1DX.push_back(this->cloneHists(allResidHistsX, "TID", " x-coord. in TID"));
0433 myResidHistsVec1DX.push_back(this->cloneHists(allResidHistsX, "TOB", " x-coord. in TOB"));
0434 myResidHistsVec1DX.push_back(this->cloneHists(allResidHistsX, "TEC", " x-coord. in TEC"));
0435
0436 for (unsigned int iHit = 0; iHit < 30; ++iHit) {
0437 for (unsigned int iHist = 0; iHist < 4 && iHist < allResidHistsX.size(); ++iHist) {
0438 TH1 *h = allResidHistsX[iHist];
0439 myResidHitHists1DX.push_back(static_cast<TH1 *>(h->Clone(Form("%s_%d", h->GetName(), iHit))));
0440 myResidHitHists1DX.back()->SetTitle(Form("%s, hit %d", h->GetTitle(), iHit));
0441 }
0442 }
0443
0444 TDirectory *dirResidX = (dirResid ? dirResid : directory)->mkdir("X", "hit residuals etc. for x measurements");
0445 this->addToDirectory(myResidHitHists1DX, dirResidX);
0446 for (std::vector<std::vector<TH1 *> >::iterator vecIter = myResidHistsVec1DX.begin(),
0447 vecIterEnd = myResidHistsVec1DX.end();
0448 vecIter != vecIterEnd;
0449 ++vecIter) {
0450 this->addToDirectory(*vecIter, dirResidX);
0451 }
0452
0453
0454 myResidHistsVec1DY.push_back(this->cloneHists(allResidHistsX, "", " y-coord."));
0455 std::vector<TH1 *> &yHists1D = allResidHistsX;
0456 myResidHistsVec1DY.push_back(this->cloneHists(yHists1D, "TPB", " y-coord. in pixel barrel"));
0457 myResidHistsVec1DY.push_back(this->cloneHists(yHists1D, "TPE", " y-coord. in pixel discs"));
0458 myResidHistsVec1DY.push_back(this->cloneHists(yHists1D, "TIB", " y-coord. in TIB"));
0459 myResidHistsVec1DY.push_back(this->cloneHists(yHists1D, "TID", " y-coord. in TID"));
0460 myResidHistsVec1DY.push_back(this->cloneHists(yHists1D, "TOB", " y-coord. in TOB"));
0461 myResidHistsVec1DY.push_back(this->cloneHists(yHists1D, "TEC", " y-coord. in TEC"));
0462 myResidHitHists1DY = this->cloneHists(myResidHitHists1DX, "", " y-coord.");
0463
0464 TDirectory *dirResidY = (dirResid ? dirResid : directory)->mkdir("Y", "hit residuals etc. for y measurements");
0465 this->addToDirectory(myResidHitHists1DY, dirResidY);
0466 for (std::vector<std::vector<TH1 *> >::iterator vecIter = myResidHistsVec1DY.begin(),
0467 vecIterEnd = myResidHistsVec1DY.end();
0468 vecIter != vecIterEnd;
0469 ++vecIter) {
0470 this->addToDirectory(*vecIter, dirResidY);
0471 }
0472
0473
0474 myFrame2FrameHists2D.push_back(
0475 new TProfile2D("frame2frame", "mean frame to frame derivatives;col;row", 6, 0., 6., 6, 0., 6.));
0476 myFrame2FrameHists2D.push_back(
0477 new TProfile2D("frame2frameAbs", "mean |frame to frame derivatives|, #neq0;col;row", 6, 0., 6., 6, 0., 6.));
0478
0479 this->equidistLogBins(logBins.GetArray(), logBins.GetSize() - 1, 1.E-36, 100.);
0480 for (unsigned int i = 0; i < 6; ++i) {
0481 for (unsigned int j = 0; j < 6; ++j) {
0482 myFrame2FrameHists2D.push_back(new TH2F(Form("frame2framePhi%d%d", i, j),
0483 Form("frame to frame derivatives, %d%d;#phi(aliDet);deriv", i, j),
0484 51,
0485 -TMath::Pi(),
0486 TMath::Pi(),
0487 10,
0488 0.,
0489 1.));
0490 myFrame2FrameHists2D.back()->SetCanExtend(TH1::kAllAxes);
0491 myFrame2FrameHists2D.push_back(new TH2F(Form("frame2frameR%d%d", i, j),
0492 Form("frame to frame derivatives, %d%d;r(aliDet);deriv", i, j),
0493 51,
0494 0.,
0495 110.,
0496 10,
0497 0.,
0498 1.));
0499 myFrame2FrameHists2D.back()->SetCanExtend(TH1::kAllAxes);
0500
0501 myFrame2FrameHists2D.push_back(
0502 new TH2F(Form("frame2framePhiLog%d%d", i, j),
0503 Form("frame to frame |derivatives|, %d%d, #neq0;#phi(aliDet);deriv", i, j),
0504 51,
0505 -TMath::Pi(),
0506 TMath::Pi(),
0507 logBins.GetSize() - 1,
0508 logBins.GetArray()));
0509 myFrame2FrameHists2D.push_back(new TH2F(Form("frame2frameRLog%d%d", i, j),
0510 Form("frame to frame |derivatives|, %d%d, #neq0;r(aliDet);deriv", i, j),
0511 51,
0512 0.,
0513 110.,
0514 logBins.GetSize() - 1,
0515 logBins.GetArray()));
0516 }
0517 }
0518
0519 TDirectory *dirF2f = directory->mkdir("frame2FrameHists", "derivatives etc.");
0520 this->addToDirectory(myFrame2FrameHists2D, dirF2f);
0521
0522 myCorrHists.push_back(new TH1F("xyCorrTPB", "#rho_{xy} in pixel barrel", 50, -.5, .5));
0523 myCorrHists.push_back(new TH1F("xyCorrTPE", "#rho_{xy} in forward pixel", 50, -.5, .5));
0524 myCorrHists.push_back(new TH1F("xyCorrTIB", "#rho_{xy} in TIB", 50, -.5, .5));
0525 myCorrHists.push_back(new TH1F("xyCorrTID", "#rho_{xy} in TID", 50, -1., 1.));
0526 myCorrHists.push_back(new TH1F("xyCorrTOB", "#rho_{xy} in TOB", 50, -.5, .5));
0527 myCorrHists.push_back(new TH1F("xyCorrTEC", "#rho_{xy} in TEC", 50, -1., 1.));
0528 TDirectory *dirCorr = directory->mkdir("hitCorrelationHists", "correlations");
0529 this->addToDirectory(myCorrHists, dirCorr);
0530
0531
0532 myPxbSurveyHists.push_back(new TH1F("PxbSurvChi2_lo", "#chi^{2} from PXB survey", 25, 0, 25));
0533 myPxbSurveyHists.push_back(new TH1F("PxbSurvChi2_md", "#chi^{2} from PXB survey", 25, 0, 500));
0534 myPxbSurveyHists.push_back(new TH1F("PxbSurvChi2_hi", "#chi^{2} from PXB survey", 25, 0, 50000));
0535 myPxbSurveyHists.push_back(new TH1F("PxbSurvChi2prob", "Math::Prob(#chi^{2},4) from PXB survey", 25, 0, 1));
0536 myPxbSurveyHists.push_back(new TH1F("PxbSurv_a0", "a_{0} from PXB survey", 100, -3000, 3000));
0537 myPxbSurveyHists.push_back(new TH1F("PxbSurv_a0_abs", "fabs(a_{0}) from PXB survey", 100, 0, 3000));
0538 myPxbSurveyHists.push_back(new TH1F("PxbSurv_a1", "a_{1} from PXB survey", 100, -6000, 6000));
0539 myPxbSurveyHists.push_back(new TH1F("PxbSurv_a1_abs", "fabs(a_{1}) from PXB survey", 100, 0, 6000));
0540 myPxbSurveyHists.push_back(
0541 new TH1F("PxbSurv_scale", "scale (#sqrt{a_{2}^{2}+a_{3}^{2}}) from PXB survey", 100, 0, 1500));
0542 myPxbSurveyHists.push_back(new TH1F("PxbSurv_phi", "angle(#atan{a_{3}/a_{4}}) from PXB survey", 100, -.05, .05));
0543 TDirectory *dirPxbSurvey = directory->mkdir("PxbSurveyHists", "PxbSurvey");
0544 this->addToDirectory(myPxbSurveyHists, dirPxbSurvey);
0545
0546 oldDir->cd();
0547 return true;
0548 }
0549
0550
0551 bool MillePedeMonitor::equidistLogBins(double *bins, int nBins, double first, double last) const {
0552
0553
0554
0555
0556
0557 if (nBins < 1 || first <= 0. || last <= 0.)
0558 return false;
0559
0560 bins[0] = first;
0561 bins[nBins] = last;
0562 const double firstLog = TMath::Log10(bins[0]);
0563 const double lastLog = TMath::Log10(bins[nBins]);
0564 for (int i = 1; i < nBins; ++i) {
0565 bins[i] = TMath::Power(10., firstLog + i * (lastLog - firstLog) / (nBins));
0566 }
0567
0568 return true;
0569 }
0570
0571
0572 void MillePedeMonitor::fillTrack(const reco::Track *track) { this->fillTrack(track, myTrackHists1D, myTrackHists2D); }
0573
0574
0575 void MillePedeMonitor::fillUsedTrack(const reco::Track *track, unsigned int nHitX, unsigned int nHitY) {
0576
0577 static const int iUsedX = this->GetIndex(myUsedTrackHists1D, "usedHitsX");
0578 myUsedTrackHists1D[iUsedX]->Fill(nHitX);
0579 static const int iUsedY = this->GetIndex(myUsedTrackHists1D, "usedHitsY");
0580 myUsedTrackHists1D[iUsedY]->Fill(nHitY);
0581
0582 if (!track)
0583 return;
0584 this->fillTrack(track, myUsedTrackHists1D, myUsedTrackHists2D);
0585 }
0586
0587
0588 void MillePedeMonitor::fillTrack(const reco::Track *track,
0589 std::vector<TH1 *> &trackHists1D,
0590 std::vector<TH2 *> &trackHists2D) {
0591 if (!track)
0592 return;
0593
0594 const reco::TrackBase::Vector p(track->momentum());
0595
0596 static const int iPtLog = this->GetIndex(trackHists1D, "ptTrackLogBins");
0597 trackHists1D[iPtLog]->Fill(track->pt());
0598 static const int iPt = this->GetIndex(trackHists1D, "ptTrack");
0599 trackHists1D[iPt]->Fill(track->pt());
0600 static const int iP = this->GetIndex(trackHists1D, "pTrack");
0601 trackHists1D[iP]->Fill(p.R());
0602 static const int iEta = this->GetIndex(trackHists1D, "etaTrack");
0603 trackHists1D[iEta]->Fill(p.Eta());
0604 static const int iTheta = this->GetIndex(trackHists1D, "thetaTrack");
0605 trackHists1D[iTheta]->Fill(p.Theta());
0606 static const int iPhi = this->GetIndex(trackHists1D, "phiTrack");
0607 trackHists1D[iPhi]->Fill(p.Phi());
0608
0609 static const int iNhit = this->GetIndex(trackHists1D, "nHitTrack");
0610 trackHists1D[iNhit]->Fill(track->numberOfValidHits());
0611 static const int iNhitInvalid = this->GetIndex(trackHists1D, "nHitInvalidTrack");
0612 trackHists1D[iNhitInvalid]->Fill(track->numberOfLostHits());
0613
0614 int nhitinTIB = 0, nhitinTOB = 0, nhitinTID = 0;
0615 int nhitinTEC = 0, nhitinBPIX = 0, nhitinFPIX = 0, nhitinPIXEL = 0;
0616 int nhitinENDCAP = 0, nhitinENDCAPplus = 0, nhitinENDCAPminus = 0;
0617 int nhitinTIDplus = 0, nhitinTIDminus = 0;
0618 int nhitinFPIXplus = 0, nhitinFPIXminus = 0;
0619 int nhitinTECplus = 0, nhitinTECminus = 0;
0620 unsigned int thishit = 0;
0621
0622 for (auto const &hit : track->recHits()) {
0623 thishit++;
0624 const DetId detId(hit->geographicalId());
0625 const int subdetId = detId.subdetId();
0626
0627 if (!hit->isValid())
0628 continue;
0629 if (detId.det() != DetId::Tracker) {
0630 edm::LogError("DetectorMismatch") << "@SUB=MillePedeMonitor::fillTrack"
0631 << "DetId.det() != DetId::Tracker (=" << DetId::Tracker << "), but "
0632 << detId.det() << ".";
0633 }
0634
0635 if (SiStripDetId::TIB == subdetId)
0636 ++nhitinTIB;
0637 else if (SiStripDetId::TOB == subdetId)
0638 ++nhitinTOB;
0639 else if (SiStripDetId::TID == subdetId) {
0640 ++nhitinTID;
0641 ++nhitinENDCAP;
0642
0643 if (trackerTopology->tidIsZMinusSide(detId)) {
0644 ++nhitinTIDminus;
0645 ++nhitinENDCAPminus;
0646 } else if (trackerTopology->tidIsZPlusSide(detId)) {
0647 ++nhitinTIDplus;
0648 ++nhitinENDCAPplus;
0649 }
0650 } else if (SiStripDetId::TEC == subdetId) {
0651 ++nhitinTEC;
0652 ++nhitinENDCAP;
0653
0654 if (trackerTopology->tecIsZMinusSide(detId)) {
0655 ++nhitinTECminus;
0656 ++nhitinENDCAPminus;
0657 } else if (trackerTopology->tecIsZPlusSide(detId)) {
0658 ++nhitinTECplus;
0659 ++nhitinENDCAPplus;
0660 }
0661 } else if (kBPIX == subdetId) {
0662 ++nhitinBPIX;
0663 ++nhitinPIXEL;
0664 } else if (kFPIX == subdetId) {
0665 ++nhitinFPIX;
0666 ++nhitinPIXEL;
0667
0668 if (trackerTopology->pxfSide(detId) == 1)
0669 ++nhitinFPIXminus;
0670 else if (trackerTopology->pxfSide(detId) == 2)
0671 ++nhitinFPIXplus;
0672 }
0673
0674 }
0675
0676 static const int iNhit01 = this->GetIndex(trackHists1D, "nHitBPIXTrack");
0677 trackHists1D[iNhit01]->Fill(nhitinBPIX);
0678 static const int iNhit02 = this->GetIndex(trackHists1D, "nHitFPIXplusTrack");
0679 trackHists1D[iNhit02]->Fill(nhitinFPIXplus);
0680 static const int iNhit03 = this->GetIndex(trackHists1D, "nHitFPIXminusTrack");
0681 trackHists1D[iNhit03]->Fill(nhitinFPIXminus);
0682 static const int iNhit04 = this->GetIndex(trackHists1D, "nHitFPIXTrack");
0683 trackHists1D[iNhit04]->Fill(nhitinFPIX);
0684 static const int iNhit05 = this->GetIndex(trackHists1D, "nHitPIXELTrack");
0685 trackHists1D[iNhit05]->Fill(nhitinPIXEL);
0686 static const int iNhit06 = this->GetIndex(trackHists1D, "nHitTIBTrack");
0687 trackHists1D[iNhit06]->Fill(nhitinTIB);
0688 static const int iNhit07 = this->GetIndex(trackHists1D, "nHitTOBTrack");
0689 trackHists1D[iNhit07]->Fill(nhitinTOB);
0690 static const int iNhit08 = this->GetIndex(trackHists1D, "nHitTIDplusTrack");
0691 trackHists1D[iNhit08]->Fill(nhitinTIDplus);
0692 static const int iNhit09 = this->GetIndex(trackHists1D, "nHitTIDminusTrack");
0693 trackHists1D[iNhit09]->Fill(nhitinTIDminus);
0694 static const int iNhit10 = this->GetIndex(trackHists1D, "nHitTIDTrack");
0695 trackHists1D[iNhit10]->Fill(nhitinTID);
0696 static const int iNhit11 = this->GetIndex(trackHists1D, "nHitTECplusTrack");
0697 trackHists1D[iNhit11]->Fill(nhitinTECplus);
0698 static const int iNhit12 = this->GetIndex(trackHists1D, "nHitTECminusTrack");
0699 trackHists1D[iNhit12]->Fill(nhitinTECminus);
0700 static const int iNhit13 = this->GetIndex(trackHists1D, "nHitTECTrack");
0701 trackHists1D[iNhit13]->Fill(nhitinTEC);
0702 static const int iNhit14 = this->GetIndex(trackHists1D, "nHitENDCAPplusTrack");
0703 trackHists1D[iNhit14]->Fill(nhitinENDCAPplus);
0704 static const int iNhit15 = this->GetIndex(trackHists1D, "nHitENDCAPminusTrack");
0705 trackHists1D[iNhit15]->Fill(nhitinENDCAPminus);
0706 static const int iNhit16 = this->GetIndex(trackHists1D, "nHitENDCAPTrack");
0707 trackHists1D[iNhit16]->Fill(nhitinENDCAP);
0708 static const int iNhit17 = this->GetIndex(trackHists2D, "nHitENDCAPTrackMinusVsPlus");
0709 trackHists2D[iNhit17]->Fill(nhitinENDCAPplus, nhitinENDCAPminus);
0710
0711 if (track->innerOk()) {
0712 const reco::TrackBase::Point &firstPoint(track->innerPosition());
0713 static const int iR1 = this->GetIndex(trackHists1D, "r1Track");
0714 trackHists1D[iR1]->Fill(firstPoint.Rho());
0715 const double rSigned1 = (firstPoint.y() > 0 ? firstPoint.Rho() : -firstPoint.Rho());
0716 static const int iR1Signed = this->GetIndex(trackHists1D, "r1TrackSigned");
0717 trackHists1D[iR1Signed]->Fill(rSigned1);
0718 static const int iZ1 = this->GetIndex(trackHists1D, "z1Track");
0719 trackHists1D[iZ1]->Fill(firstPoint.Z());
0720 static const int iZ1Full = this->GetIndex(trackHists1D, "z1TrackFull");
0721 trackHists1D[iZ1Full]->Fill(firstPoint.Z());
0722 static const int iY1 = this->GetIndex(trackHists1D, "y1Track");
0723 trackHists1D[iY1]->Fill(firstPoint.Y());
0724 static const int iPhi1 = this->GetIndex(trackHists1D, "phi1Track");
0725 trackHists1D[iPhi1]->Fill(firstPoint.phi());
0726 static const int iRz1Full = this->GetIndex(trackHists2D, "rz1TrackFull");
0727 trackHists2D[iRz1Full]->Fill(firstPoint.Z(), rSigned1);
0728 static const int iXy1 = this->GetIndex(trackHists2D, "xy1Track");
0729 trackHists2D[iXy1]->Fill(firstPoint.X(), firstPoint.Y());
0730 }
0731
0732 if (track->outerOk()) {
0733 const reco::TrackBase::Point &lastPoint(track->outerPosition());
0734 static const int iRlast = this->GetIndex(trackHists1D, "rLastTrack");
0735 trackHists1D[iRlast]->Fill(lastPoint.Rho());
0736 static const int iZlast = this->GetIndex(trackHists1D, "zLastTrack");
0737 trackHists1D[iZlast]->Fill(lastPoint.Z());
0738 static const int iYlast = this->GetIndex(trackHists1D, "yLastTrack");
0739 trackHists1D[iYlast]->Fill(lastPoint.Y());
0740 static const int iPhiLast = this->GetIndex(trackHists1D, "phiLastTrack");
0741 trackHists1D[iPhiLast]->Fill(lastPoint.phi());
0742 }
0743
0744 static const int iChi2Ndf = this->GetIndex(trackHists1D, "chi2PerNdf");
0745 trackHists1D[iChi2Ndf]->Fill(track->normalizedChi2());
0746
0747 static const int iImpZ = this->GetIndex(trackHists1D, "impParZ");
0748 trackHists1D[iImpZ]->Fill(track->dz());
0749 static const int iImpZerr = this->GetIndex(trackHists1D, "impParErrZ");
0750 trackHists1D[iImpZerr]->Fill(track->dzError());
0751
0752 static const int iImpRphi = this->GetIndex(trackHists1D, "impParRphi");
0753 trackHists1D[iImpRphi]->Fill(track->d0());
0754 static const int iImpRphiErr = this->GetIndex(trackHists1D, "impParErrRphi");
0755 trackHists1D[iImpRphiErr]->Fill(track->d0Error());
0756 }
0757
0758
0759 void MillePedeMonitor::fillRefTrajectory(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr) {
0760 static const int iValid = this->GetIndex(myTrajectoryHists1D, "validRefTraj");
0761 if (refTrajPtr->isValid()) {
0762 myTrajectoryHists1D[iValid]->Fill(1.);
0763 } else {
0764 myTrajectoryHists1D[iValid]->Fill(0.);
0765 return;
0766 }
0767
0768 const AlgebraicSymMatrix &covMeasLoc = refTrajPtr->measurementErrors();
0769 const AlgebraicVector &measurementsLoc = refTrajPtr->measurements();
0770 const AlgebraicVector &trajectoryLoc = refTrajPtr->trajectoryPositions();
0771 const TransientTrackingRecHit::ConstRecHitContainer &recHits = refTrajPtr->recHits();
0772 const AlgebraicMatrix &derivatives = refTrajPtr->derivatives();
0773
0774
0775 const int nRow = refTrajPtr->numberOfHitMeas();
0776
0777 for (int iRow = 0; iRow < nRow; ++iRow) {
0778 const double residuum = measurementsLoc[iRow] - trajectoryLoc[iRow];
0779 const double covMeasLocRow = covMeasLoc[iRow][iRow];
0780 const bool is2DhitRow = (!recHits[iRow / 2]->detUnit()
0781 || recHits[iRow / 2]->detUnit()->type().isTrackerPixel());
0782
0783 if (TMath::Even(iRow)) {
0784 static const int iMeasLocX = this->GetIndex(myTrajectoryHists1D, "measLocX");
0785 myTrajectoryHists1D[iMeasLocX]->Fill(measurementsLoc[iRow]);
0786 static const int iTrajLocX = this->GetIndex(myTrajectoryHists1D, "trajLocX");
0787 myTrajectoryHists1D[iTrajLocX]->Fill(trajectoryLoc[iRow]);
0788 static const int iResidLocX = this->GetIndex(myTrajectoryHists1D, "residLocX");
0789 myTrajectoryHists1D[iResidLocX]->Fill(residuum);
0790
0791 static const int iMeasLocXvsHit = this->GetIndex(myTrajectoryHists2D, "measLocXvsHit");
0792 myTrajectoryHists2D[iMeasLocXvsHit]->Fill(iRow / 2, measurementsLoc[iRow]);
0793 static const int iTrajLocXvsHit = this->GetIndex(myTrajectoryHists2D, "trajLocXvsHit");
0794 myTrajectoryHists2D[iTrajLocXvsHit]->Fill(iRow / 2, trajectoryLoc[iRow]);
0795 static const int iResidLocXvsHit = this->GetIndex(myTrajectoryHists2D, "residLocXvsHit");
0796 myTrajectoryHists2D[iResidLocXvsHit]->Fill(iRow / 2, residuum);
0797
0798 if (covMeasLocRow > 0.) {
0799 static const int iReduResidLocX = this->GetIndex(myTrajectoryHists1D, "reduResidLocX");
0800 myTrajectoryHists1D[iReduResidLocX]->Fill(residuum / TMath::Sqrt(covMeasLocRow));
0801 static const int iReduResidLocXvsHit = this->GetIndex(myTrajectoryHists2D, "reduResidLocXvsHit");
0802 myTrajectoryHists2D[iReduResidLocXvsHit]->Fill(iRow / 2, residuum / TMath::Sqrt(covMeasLocRow));
0803 }
0804 } else if (is2DhitRow) {
0805
0806 static const int iMeasLocY = this->GetIndex(myTrajectoryHists1D, "measLocY");
0807 myTrajectoryHists1D[iMeasLocY]->Fill(measurementsLoc[iRow]);
0808 static const int iTrajLocY = this->GetIndex(myTrajectoryHists1D, "trajLocY");
0809 myTrajectoryHists1D[iTrajLocY]->Fill(trajectoryLoc[iRow]);
0810 static const int iResidLocY = this->GetIndex(myTrajectoryHists1D, "residLocY");
0811 myTrajectoryHists1D[iResidLocY]->Fill(residuum);
0812
0813 static const int iMeasLocYvsHit = this->GetIndex(myTrajectoryHists2D, "measLocYvsHit");
0814 myTrajectoryHists2D[iMeasLocYvsHit]->Fill(iRow / 2, measurementsLoc[iRow]);
0815 static const int iTrajLocYvsHit = this->GetIndex(myTrajectoryHists2D, "trajLocYvsHit");
0816 myTrajectoryHists2D[iTrajLocYvsHit]->Fill(iRow / 2, trajectoryLoc[iRow]);
0817 static const int iResidLocYvsHit = this->GetIndex(myTrajectoryHists2D, "residLocYvsHit");
0818 myTrajectoryHists2D[iResidLocYvsHit]->Fill(iRow / 2, residuum);
0819
0820 if (covMeasLocRow > 0.) {
0821 static const int iReduResidLocY = this->GetIndex(myTrajectoryHists1D, "reduResidLocY");
0822 myTrajectoryHists1D[iReduResidLocY]->Fill(residuum / TMath::Sqrt(covMeasLocRow));
0823 static const int iReduResidLocYvsHit = this->GetIndex(myTrajectoryHists2D, "reduResidLocYvsHit");
0824 myTrajectoryHists2D[iReduResidLocYvsHit]->Fill(iRow / 2, residuum / TMath::Sqrt(covMeasLocRow));
0825 }
0826 }
0827
0828 float nHitRow = iRow / 2;
0829 if (TMath::Odd(iRow))
0830 nHitRow += 0.5;
0831
0832 for (int iCol = iRow + 1; iCol < nRow; ++iCol) {
0833 double rho = TMath::Sqrt(covMeasLocRow + covMeasLoc[iCol][iCol]);
0834 rho = (0. == rho ? -2 : covMeasLoc[iRow][iCol] / rho);
0835 float nHitCol = iCol / 2;
0836 if (TMath::Odd(iCol))
0837 nHitCol += 0.5;
0838
0839 if (0. == rho)
0840 continue;
0841 static const int iProfileCorr = this->GetIndex(myTrajectoryHists2D, "profCorr");
0842 myTrajectoryHists2D[iProfileCorr]->Fill(nHitRow, nHitCol, TMath::Abs(rho));
0843 if (iRow + 1 == iCol && TMath::Even(iRow)) {
0844
0845
0846 static const int iHitCorrXy = this->GetIndex(myTrajectoryHists2D, "hitCorrXy");
0847 myTrajectoryHists2D[iHitCorrXy]->Fill(iRow / 2, rho);
0848 static const int iHitCorrXyLog = this->GetIndex(myTrajectoryHists2D, "hitCorrXyLog");
0849 myTrajectoryHists2D[iHitCorrXyLog]->Fill(iRow / 2, TMath::Abs(rho));
0850 if (is2DhitRow) {
0851 static const int iHitCorrXyValid = this->GetIndex(myTrajectoryHists2D, "hitCorrXyValid");
0852 myTrajectoryHists2D[iHitCorrXyValid]->Fill(iRow / 2, rho);
0853 static const int iHitCorrXyLogValid = this->GetIndex(myTrajectoryHists2D, "hitCorrXyLogValid");
0854 myTrajectoryHists2D[iHitCorrXyLogValid]->Fill(iRow / 2, TMath::Abs(rho));
0855 }
0856 } else {
0857 static const int iProfCorrOffXy = this->GetIndex(myTrajectoryHists2D, "profCorrOffXy");
0858 myTrajectoryHists2D[iProfCorrOffXy]->Fill(nHitRow, nHitCol, TMath::Abs(rho));
0859 static const int iHitCorrOffBlock = this->GetIndex(myTrajectoryHists2D, "hitCorrOffBlock");
0860 myTrajectoryHists2D[iHitCorrOffBlock]->Fill(nHitRow, rho);
0861 static const int iHitCorOffBlkLg = this->GetIndex(myTrajectoryHists2D, "hitCorrOffBlockLog");
0862 myTrajectoryHists2D[iHitCorOffBlkLg]->Fill(nHitRow, TMath::Abs(rho));
0863 }
0864 }
0865
0866
0867 for (int iCol = 0; iCol < derivatives.num_col(); ++iCol) {
0868 static const int iProfDerivatives = this->GetIndex(myTrajectoryHists2D, "profDerivatives");
0869 myTrajectoryHists2D[iProfDerivatives]->Fill(nHitRow, iCol, derivatives[iRow][iCol]);
0870 static const int iDerivatives = this->GetIndex(myTrajectoryHists2D, "derivatives");
0871 myTrajectoryHists2D[iDerivatives]->Fill(iCol, derivatives[iRow][iCol]);
0872 static const int iDerivativesLog = this->GetIndex(myTrajectoryHists2D, "derivativesLog");
0873 myTrajectoryHists2D[iDerivativesLog]->Fill(iCol, TMath::Abs(derivatives[iRow][iCol]));
0874 static const int iDerivativesPhi = this->GetIndex(myTrajectoryHists2D, "derivativesVsPhi");
0875 myTrajectoryHists2D[iDerivativesPhi]->Fill(recHits[iRow / 2]->det()->position().phi(), derivatives[iRow][iCol]);
0876 }
0877 }
0878 }
0879
0880
0881 void MillePedeMonitor::fillDerivatives(const ConstRecHitPointer &recHit,
0882 const float *localDerivs,
0883 unsigned int nLocal,
0884 const float *globalDerivs,
0885 unsigned int nGlobal,
0886 const int *labels) {
0887 const double phi = recHit->det()->position().phi();
0888
0889 static const int iLocPar = this->GetIndex(myDerivHists2D, "localDerivsPar");
0890 static const int iLocPhi = this->GetIndex(myDerivHists2D, "localDerivsPhi");
0891 static const int iLocParLog = this->GetIndex(myDerivHists2D, "localDerivsParLog");
0892 static const int iLocPhiLog = this->GetIndex(myDerivHists2D, "localDerivsPhiLog");
0893 for (unsigned int i = 0; i < nLocal; ++i) {
0894 myDerivHists2D[iLocPar]->Fill(i, localDerivs[i]);
0895 myDerivHists2D[iLocPhi]->Fill(phi, localDerivs[i]);
0896 if (localDerivs[i]) {
0897 myDerivHists2D[iLocParLog]->Fill(i, TMath::Abs(localDerivs[i]));
0898 myDerivHists2D[iLocPhiLog]->Fill(phi, TMath::Abs(localDerivs[i]));
0899 }
0900 }
0901
0902 static const int iGlobPar = this->GetIndex(myDerivHists2D, "globalDerivsPar");
0903 static const int iGlobPhi = this->GetIndex(myDerivHists2D, "globalDerivsPhi");
0904 static const int iGlobParLog = this->GetIndex(myDerivHists2D, "globalDerivsParLog");
0905 static const int iGlobPhiLog = this->GetIndex(myDerivHists2D, "globalDerivsPhiLog");
0906
0907 for (unsigned int i = 0; i < nGlobal; ++i) {
0908 const unsigned int parNum = (labels ? (labels[i] % PedeLabelerBase::theMaxNumParam) - 1 : i);
0909 myDerivHists2D[iGlobPar]->Fill(parNum, globalDerivs[i]);
0910 myDerivHists2D[iGlobPhi]->Fill(phi, globalDerivs[i]);
0911 if (globalDerivs[i]) {
0912 myDerivHists2D[iGlobParLog]->Fill(parNum, TMath::Abs(globalDerivs[i]));
0913 myDerivHists2D[iGlobPhiLog]->Fill(phi, TMath::Abs(globalDerivs[i]));
0914
0915 }
0916 }
0917 }
0918
0919
0920 void MillePedeMonitor::fillResiduals(const ConstRecHitPointer &recHit,
0921 const TrajectoryStateOnSurface &tsos,
0922 unsigned int nHit,
0923 float residuum,
0924 float sigma,
0925 bool isY) {
0926
0927
0928 const GlobalPoint detPos(recHit->det()->position());
0929 const double phi = detPos.phi();
0930
0931
0932 static const int iResPhi = this->GetIndex(myResidHists2D, "residPhi");
0933 myResidHists2D[iResPhi]->Fill(phi, residuum);
0934 static const int iSigPhi = this->GetIndex(myResidHists2D, "sigmaPhi");
0935 myResidHists2D[iSigPhi]->Fill(phi, sigma);
0936 if (sigma) {
0937 static const int iReduResPhi = this->GetIndex(myResidHists2D, "reduResidPhi");
0938 myResidHists2D[iReduResPhi]->Fill(phi, residuum / sigma);
0939 }
0940
0941
0942
0943
0944
0945
0946
0947
0948
0949
0950
0951
0952
0953
0954
0955
0956
0957
0958
0959
0960
0961
0962
0963
0964
0965
0966
0967
0968
0969
0970
0971
0972
0973 const LocalVector mom(tsos.localDirection());
0974 const float phiSensToNorm = TMath::ATan(TMath::Abs((isY ? mom.y() : mom.x()) / mom.z()));
0975
0976 std::vector<std::vector<TH1 *> > &histArrayVec = (isY ? myResidHistsVec1DY : myResidHistsVec1DX);
0977 std::vector<TH1 *> &hitHists = (isY ? myResidHitHists1DY : myResidHitHists1DX);
0978
0979
0980 this->fillResidualHists(histArrayVec[0], phiSensToNorm, residuum, sigma);
0981 this->fillResidualHitHists(hitHists, phiSensToNorm, residuum, sigma, nHit);
0982
0983 const DetId detId(recHit->det()->geographicalId());
0984 if (detId.det() == DetId::Tracker) {
0985
0986 unsigned int subDetNum = detId.subdetId();
0987 if (subDetNum < histArrayVec.size() && subDetNum > 0) {
0988 this->fillResidualHists(histArrayVec[subDetNum], phiSensToNorm, residuum, sigma);
0989 } else {
0990 if (detId != AlignableBeamSpot::detId())
0991 edm::LogWarning("Alignment") << "@SUB=MillePedeMonitor::fillResiduals"
0992 << "Expect subDetNum from 1 to 6, got " << subDetNum;
0993 }
0994 }
0995 }
0996
0997
0998 void MillePedeMonitor::fillResidualHists(const std::vector<TH1 *> &hists,
0999 float phiSensToNorm,
1000 float residuum,
1001 float sigma) {
1002
1003
1004
1005
1006 static const int iRes = this->GetIndex(hists, "resid");
1007 hists[iRes]->Fill(residuum);
1008 static const int iSigma = this->GetIndex(hists, "sigma");
1009 hists[iSigma]->Fill(sigma);
1010 static const int iSigmaVsAngle = this->GetIndex(hists, "sigmaVsAngle");
1011 hists[iSigmaVsAngle]->Fill(phiSensToNorm, sigma);
1012 static const int iResidVsAngle = this->GetIndex(hists, "residVsAngle");
1013 hists[iResidVsAngle]->Fill(phiSensToNorm, residuum);
1014 static const int iReduRes = this->GetIndex(hists, "reduResid");
1015 static const int iReduResidVsAngle = this->GetIndex(hists, "reduResidVsAngle");
1016 if (sigma) {
1017 hists[iReduRes]->Fill(residuum / sigma);
1018 hists[iReduResidVsAngle]->Fill(phiSensToNorm, residuum / sigma);
1019 }
1020 static const int iAngle = this->GetIndex(hists, "angle");
1021 hists[iAngle]->Fill(phiSensToNorm);
1022
1023 if (phiSensToNorm > TMath::DegToRad() * 45.) {
1024 static const int iResGt45 = this->GetIndex(hists, "residGt45");
1025 hists[iResGt45]->Fill(residuum);
1026 static const int iSigmaGt45 = this->GetIndex(hists, "sigmaGt45");
1027 hists[iSigmaGt45]->Fill(sigma);
1028 static const int iReduResGt45 = this->GetIndex(hists, "reduResidGt45");
1029 if (sigma)
1030 hists[iReduResGt45]->Fill(residuum / sigma);
1031 } else {
1032 static const int iResLt45 = this->GetIndex(hists, "residLt45");
1033 hists[iResLt45]->Fill(residuum);
1034 static const int iSigmaLt45 = this->GetIndex(hists, "sigmaLt45");
1035 hists[iSigmaLt45]->Fill(sigma);
1036 static const int iReduResLt45 = this->GetIndex(hists, "reduResidLt45");
1037 if (sigma)
1038 hists[iReduResLt45]->Fill(residuum / sigma);
1039 }
1040 }
1041
1042
1043 void MillePedeMonitor::fillResidualHitHists(
1044 const std::vector<TH1 *> &hists, float angle, float residuum, float sigma, unsigned int nHit) {
1045
1046
1047
1048
1049 static constexpr unsigned int maxNhit = 29;
1050 static const auto iResHit = indexArray1D<TH1, maxNhit + 1>(hists, "resid_%d");
1051 static const auto iSigmaHit = indexArray1D<TH1, maxNhit + 1>(hists, "sigma_%d");
1052 static const auto iReduResHit = indexArray1D<TH1, maxNhit + 1>(hists, "reduResid_%d");
1053 static const auto iAngleHit = indexArray1D<TH1, maxNhit + 1>(hists, "angle_%d");
1054 if (nHit > maxNhit)
1055 nHit = maxNhit;
1056
1057 hists[iResHit[nHit]]->Fill(residuum);
1058 hists[iSigmaHit[nHit]]->Fill(sigma);
1059 if (sigma) {
1060 hists[iReduResHit[nHit]]->Fill(residuum / sigma);
1061 }
1062 hists[iAngleHit[nHit]]->Fill(angle);
1063 }
1064
1065
1066 void MillePedeMonitor::fillFrameToFrame(const AlignableDetOrUnitPtr &aliDet, const Alignable *ali) {
1067
1068 FrameToFrameDerivative ftfd;
1069 const AlgebraicMatrix frameDeriv = ftfd.frameToFrameDerivative(aliDet, ali);
1070
1071
1072 static const int iF2f = this->GetIndex(myFrame2FrameHists2D, "frame2frame");
1073 static const int iF2fAbs = this->GetIndex(myFrame2FrameHists2D, "frame2frameAbs");
1074 static const auto iF2fIjPhi = indexArray2D<TH2, 6>(myFrame2FrameHists2D, "frame2framePhi%d%d");
1075 static const auto iF2fIjPhiLog = indexArray2D<TH2, 6>(myFrame2FrameHists2D, "frame2framePhiLog%d%d");
1076 static const auto iF2fIjR = indexArray2D<TH2, 6>(myFrame2FrameHists2D, "frame2frameR%d%d");
1077 static const auto iF2fIjRLog = indexArray2D<TH2, 6>(myFrame2FrameHists2D, "frame2frameRLog%d%d");
1078
1079 const double phi = aliDet->globalPosition().phi();
1080 const double r = aliDet->globalPosition().perp();
1081 for (unsigned int i = 0; i < 6; ++i) {
1082 for (unsigned int j = 0; j < 6; ++j) {
1083 myFrame2FrameHists2D[iF2f]->Fill(i, j, frameDeriv[i][j]);
1084 myFrame2FrameHists2D[iF2fIjPhi[i][j]]->Fill(phi, frameDeriv[i][j]);
1085 myFrame2FrameHists2D[iF2fIjR[i][j]]->Fill(r, frameDeriv[i][j]);
1086 if (frameDeriv[i][j]) {
1087 myFrame2FrameHists2D[iF2fAbs]->Fill(i, j, TMath::Abs(frameDeriv[i][j]));
1088 myFrame2FrameHists2D[iF2fIjPhiLog[i][j]]->Fill(phi, TMath::Abs(frameDeriv[i][j]));
1089 myFrame2FrameHists2D[iF2fIjRLog[i][j]]->Fill(r, TMath::Abs(frameDeriv[i][j]));
1090 }
1091 }
1092 }
1093 }
1094
1095
1096 void MillePedeMonitor::fillCorrelations2D(float corr, const ConstRecHitPointer &recHit) {
1097 const DetId detId(recHit->det()->geographicalId());
1098 if (detId.det() != DetId::Tracker)
1099 return;
1100
1101 if ((detId.subdetId() < 1 || detId.subdetId() > 6) && detId != AlignableBeamSpot::detId()) {
1102 edm::LogWarning("Alignment") << "@SUB=MillePedeMonitor::fillCorrelations2D"
1103 << "Expect subdetId from 1 to 6, got " << detId.subdetId();
1104 return;
1105 }
1106
1107 myCorrHists[detId.subdetId() - 1]->Fill(corr);
1108 }
1109
1110
1111 void MillePedeMonitor::fillPxbSurveyHistsChi2(const float &chi2) {
1112 static const int iPxbSurvChi2_lo = this->GetIndex(myPxbSurveyHists, "PxbSurvChi2_lo");
1113 myPxbSurveyHists[iPxbSurvChi2_lo]->Fill(chi2);
1114 static const int iPxbSurvChi2_md = this->GetIndex(myPxbSurveyHists, "PxbSurvChi2_md");
1115 myPxbSurveyHists[iPxbSurvChi2_md]->Fill(chi2);
1116 static const int iPxbSurvChi2_hi = this->GetIndex(myPxbSurveyHists, "PxbSurvChi2_hi");
1117 myPxbSurveyHists[iPxbSurvChi2_hi]->Fill(chi2);
1118 static const int iPxbSurvChi2prob = this->GetIndex(myPxbSurveyHists, "PxbSurvChi2prob");
1119 myPxbSurveyHists[iPxbSurvChi2prob]->Fill(TMath::Prob(chi2, 4));
1120 }
1121
1122
1123 void MillePedeMonitor::fillPxbSurveyHistsLocalPars(const float &a0, const float &a1, const float &S, const float &phi) {
1124 static const int iPxbSurv_a0 = this->GetIndex(myPxbSurveyHists, "PxbSurv_a0");
1125 myPxbSurveyHists[iPxbSurv_a0]->Fill(a0);
1126 static const int iPxbSurv_a0_abs = this->GetIndex(myPxbSurveyHists, "PxbSurv_a0_abs");
1127 myPxbSurveyHists[iPxbSurv_a0_abs]->Fill(fabs(a0));
1128 static const int iPxbSurv_a1 = this->GetIndex(myPxbSurveyHists, "PxbSurv_a1");
1129 myPxbSurveyHists[iPxbSurv_a1]->Fill(a1);
1130 static const int iPxbSurv_a1_abs = this->GetIndex(myPxbSurveyHists, "PxbSurv_a1_abs");
1131 myPxbSurveyHists[iPxbSurv_a1_abs]->Fill(fabs(a1));
1132 static const int iPxbSurv_scale = this->GetIndex(myPxbSurveyHists, "PxbSurv_scale");
1133 myPxbSurveyHists[iPxbSurv_scale]->Fill(S);
1134 static const int iPxbSurv_phi = this->GetIndex(myPxbSurveyHists, "PxbSurv_phi");
1135 myPxbSurveyHists[iPxbSurv_phi]->Fill(phi);
1136 }