File indexing completed on 2024-04-06 11:56:37
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
0621 for (auto const &hit : track->recHits()) {
0622 const DetId detId(hit->geographicalId());
0623 const int subdetId = detId.subdetId();
0624
0625 if (!hit->isValid())
0626 continue;
0627 if (detId.det() != DetId::Tracker) {
0628 edm::LogError("DetectorMismatch") << "@SUB=MillePedeMonitor::fillTrack"
0629 << "DetId.det() != DetId::Tracker (=" << DetId::Tracker << "), but "
0630 << detId.det() << ".";
0631 }
0632
0633 if (SiStripDetId::TIB == subdetId)
0634 ++nhitinTIB;
0635 else if (SiStripDetId::TOB == subdetId)
0636 ++nhitinTOB;
0637 else if (SiStripDetId::TID == subdetId) {
0638 ++nhitinTID;
0639 ++nhitinENDCAP;
0640
0641 if (trackerTopology->tidIsZMinusSide(detId)) {
0642 ++nhitinTIDminus;
0643 ++nhitinENDCAPminus;
0644 } else if (trackerTopology->tidIsZPlusSide(detId)) {
0645 ++nhitinTIDplus;
0646 ++nhitinENDCAPplus;
0647 }
0648 } else if (SiStripDetId::TEC == subdetId) {
0649 ++nhitinTEC;
0650 ++nhitinENDCAP;
0651
0652 if (trackerTopology->tecIsZMinusSide(detId)) {
0653 ++nhitinTECminus;
0654 ++nhitinENDCAPminus;
0655 } else if (trackerTopology->tecIsZPlusSide(detId)) {
0656 ++nhitinTECplus;
0657 ++nhitinENDCAPplus;
0658 }
0659 } else if (kBPIX == subdetId) {
0660 ++nhitinBPIX;
0661 ++nhitinPIXEL;
0662 } else if (kFPIX == subdetId) {
0663 ++nhitinFPIX;
0664 ++nhitinPIXEL;
0665
0666 if (trackerTopology->pxfSide(detId) == 1)
0667 ++nhitinFPIXminus;
0668 else if (trackerTopology->pxfSide(detId) == 2)
0669 ++nhitinFPIXplus;
0670 }
0671
0672 }
0673
0674 static const int iNhit01 = this->GetIndex(trackHists1D, "nHitBPIXTrack");
0675 trackHists1D[iNhit01]->Fill(nhitinBPIX);
0676 static const int iNhit02 = this->GetIndex(trackHists1D, "nHitFPIXplusTrack");
0677 trackHists1D[iNhit02]->Fill(nhitinFPIXplus);
0678 static const int iNhit03 = this->GetIndex(trackHists1D, "nHitFPIXminusTrack");
0679 trackHists1D[iNhit03]->Fill(nhitinFPIXminus);
0680 static const int iNhit04 = this->GetIndex(trackHists1D, "nHitFPIXTrack");
0681 trackHists1D[iNhit04]->Fill(nhitinFPIX);
0682 static const int iNhit05 = this->GetIndex(trackHists1D, "nHitPIXELTrack");
0683 trackHists1D[iNhit05]->Fill(nhitinPIXEL);
0684 static const int iNhit06 = this->GetIndex(trackHists1D, "nHitTIBTrack");
0685 trackHists1D[iNhit06]->Fill(nhitinTIB);
0686 static const int iNhit07 = this->GetIndex(trackHists1D, "nHitTOBTrack");
0687 trackHists1D[iNhit07]->Fill(nhitinTOB);
0688 static const int iNhit08 = this->GetIndex(trackHists1D, "nHitTIDplusTrack");
0689 trackHists1D[iNhit08]->Fill(nhitinTIDplus);
0690 static const int iNhit09 = this->GetIndex(trackHists1D, "nHitTIDminusTrack");
0691 trackHists1D[iNhit09]->Fill(nhitinTIDminus);
0692 static const int iNhit10 = this->GetIndex(trackHists1D, "nHitTIDTrack");
0693 trackHists1D[iNhit10]->Fill(nhitinTID);
0694 static const int iNhit11 = this->GetIndex(trackHists1D, "nHitTECplusTrack");
0695 trackHists1D[iNhit11]->Fill(nhitinTECplus);
0696 static const int iNhit12 = this->GetIndex(trackHists1D, "nHitTECminusTrack");
0697 trackHists1D[iNhit12]->Fill(nhitinTECminus);
0698 static const int iNhit13 = this->GetIndex(trackHists1D, "nHitTECTrack");
0699 trackHists1D[iNhit13]->Fill(nhitinTEC);
0700 static const int iNhit14 = this->GetIndex(trackHists1D, "nHitENDCAPplusTrack");
0701 trackHists1D[iNhit14]->Fill(nhitinENDCAPplus);
0702 static const int iNhit15 = this->GetIndex(trackHists1D, "nHitENDCAPminusTrack");
0703 trackHists1D[iNhit15]->Fill(nhitinENDCAPminus);
0704 static const int iNhit16 = this->GetIndex(trackHists1D, "nHitENDCAPTrack");
0705 trackHists1D[iNhit16]->Fill(nhitinENDCAP);
0706 static const int iNhit17 = this->GetIndex(trackHists2D, "nHitENDCAPTrackMinusVsPlus");
0707 trackHists2D[iNhit17]->Fill(nhitinENDCAPplus, nhitinENDCAPminus);
0708
0709 if (track->innerOk()) {
0710 const reco::TrackBase::Point &firstPoint(track->innerPosition());
0711 static const int iR1 = this->GetIndex(trackHists1D, "r1Track");
0712 trackHists1D[iR1]->Fill(firstPoint.Rho());
0713 const double rSigned1 = (firstPoint.y() > 0 ? firstPoint.Rho() : -firstPoint.Rho());
0714 static const int iR1Signed = this->GetIndex(trackHists1D, "r1TrackSigned");
0715 trackHists1D[iR1Signed]->Fill(rSigned1);
0716 static const int iZ1 = this->GetIndex(trackHists1D, "z1Track");
0717 trackHists1D[iZ1]->Fill(firstPoint.Z());
0718 static const int iZ1Full = this->GetIndex(trackHists1D, "z1TrackFull");
0719 trackHists1D[iZ1Full]->Fill(firstPoint.Z());
0720 static const int iY1 = this->GetIndex(trackHists1D, "y1Track");
0721 trackHists1D[iY1]->Fill(firstPoint.Y());
0722 static const int iPhi1 = this->GetIndex(trackHists1D, "phi1Track");
0723 trackHists1D[iPhi1]->Fill(firstPoint.phi());
0724 static const int iRz1Full = this->GetIndex(trackHists2D, "rz1TrackFull");
0725 trackHists2D[iRz1Full]->Fill(firstPoint.Z(), rSigned1);
0726 static const int iXy1 = this->GetIndex(trackHists2D, "xy1Track");
0727 trackHists2D[iXy1]->Fill(firstPoint.X(), firstPoint.Y());
0728 }
0729
0730 if (track->outerOk()) {
0731 const reco::TrackBase::Point &lastPoint(track->outerPosition());
0732 static const int iRlast = this->GetIndex(trackHists1D, "rLastTrack");
0733 trackHists1D[iRlast]->Fill(lastPoint.Rho());
0734 static const int iZlast = this->GetIndex(trackHists1D, "zLastTrack");
0735 trackHists1D[iZlast]->Fill(lastPoint.Z());
0736 static const int iYlast = this->GetIndex(trackHists1D, "yLastTrack");
0737 trackHists1D[iYlast]->Fill(lastPoint.Y());
0738 static const int iPhiLast = this->GetIndex(trackHists1D, "phiLastTrack");
0739 trackHists1D[iPhiLast]->Fill(lastPoint.phi());
0740 }
0741
0742 static const int iChi2Ndf = this->GetIndex(trackHists1D, "chi2PerNdf");
0743 trackHists1D[iChi2Ndf]->Fill(track->normalizedChi2());
0744
0745 static const int iImpZ = this->GetIndex(trackHists1D, "impParZ");
0746 trackHists1D[iImpZ]->Fill(track->dz());
0747 static const int iImpZerr = this->GetIndex(trackHists1D, "impParErrZ");
0748 trackHists1D[iImpZerr]->Fill(track->dzError());
0749
0750 static const int iImpRphi = this->GetIndex(trackHists1D, "impParRphi");
0751 trackHists1D[iImpRphi]->Fill(track->d0());
0752 static const int iImpRphiErr = this->GetIndex(trackHists1D, "impParErrRphi");
0753 trackHists1D[iImpRphiErr]->Fill(track->d0Error());
0754 }
0755
0756
0757 void MillePedeMonitor::fillRefTrajectory(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr) {
0758 static const int iValid = this->GetIndex(myTrajectoryHists1D, "validRefTraj");
0759 if (refTrajPtr->isValid()) {
0760 myTrajectoryHists1D[iValid]->Fill(1.);
0761 } else {
0762 myTrajectoryHists1D[iValid]->Fill(0.);
0763 return;
0764 }
0765
0766 const AlgebraicSymMatrix &covMeasLoc = refTrajPtr->measurementErrors();
0767 const AlgebraicVector &measurementsLoc = refTrajPtr->measurements();
0768 const AlgebraicVector &trajectoryLoc = refTrajPtr->trajectoryPositions();
0769 const TransientTrackingRecHit::ConstRecHitContainer &recHits = refTrajPtr->recHits();
0770 const AlgebraicMatrix &derivatives = refTrajPtr->derivatives();
0771
0772
0773 const int nRow = refTrajPtr->numberOfHitMeas();
0774
0775 for (int iRow = 0; iRow < nRow; ++iRow) {
0776 const double residuum = measurementsLoc[iRow] - trajectoryLoc[iRow];
0777 const double covMeasLocRow = covMeasLoc[iRow][iRow];
0778 const bool is2DhitRow = (!recHits[iRow / 2]->detUnit()
0779 || recHits[iRow / 2]->detUnit()->type().isTrackerPixel());
0780
0781 if (TMath::Even(iRow)) {
0782 static const int iMeasLocX = this->GetIndex(myTrajectoryHists1D, "measLocX");
0783 myTrajectoryHists1D[iMeasLocX]->Fill(measurementsLoc[iRow]);
0784 static const int iTrajLocX = this->GetIndex(myTrajectoryHists1D, "trajLocX");
0785 myTrajectoryHists1D[iTrajLocX]->Fill(trajectoryLoc[iRow]);
0786 static const int iResidLocX = this->GetIndex(myTrajectoryHists1D, "residLocX");
0787 myTrajectoryHists1D[iResidLocX]->Fill(residuum);
0788
0789 static const int iMeasLocXvsHit = this->GetIndex(myTrajectoryHists2D, "measLocXvsHit");
0790 myTrajectoryHists2D[iMeasLocXvsHit]->Fill(iRow / 2, measurementsLoc[iRow]);
0791 static const int iTrajLocXvsHit = this->GetIndex(myTrajectoryHists2D, "trajLocXvsHit");
0792 myTrajectoryHists2D[iTrajLocXvsHit]->Fill(iRow / 2, trajectoryLoc[iRow]);
0793 static const int iResidLocXvsHit = this->GetIndex(myTrajectoryHists2D, "residLocXvsHit");
0794 myTrajectoryHists2D[iResidLocXvsHit]->Fill(iRow / 2, residuum);
0795
0796 if (covMeasLocRow > 0.) {
0797 static const int iReduResidLocX = this->GetIndex(myTrajectoryHists1D, "reduResidLocX");
0798 myTrajectoryHists1D[iReduResidLocX]->Fill(residuum / TMath::Sqrt(covMeasLocRow));
0799 static const int iReduResidLocXvsHit = this->GetIndex(myTrajectoryHists2D, "reduResidLocXvsHit");
0800 myTrajectoryHists2D[iReduResidLocXvsHit]->Fill(iRow / 2, residuum / TMath::Sqrt(covMeasLocRow));
0801 }
0802 } else if (is2DhitRow) {
0803
0804 static const int iMeasLocY = this->GetIndex(myTrajectoryHists1D, "measLocY");
0805 myTrajectoryHists1D[iMeasLocY]->Fill(measurementsLoc[iRow]);
0806 static const int iTrajLocY = this->GetIndex(myTrajectoryHists1D, "trajLocY");
0807 myTrajectoryHists1D[iTrajLocY]->Fill(trajectoryLoc[iRow]);
0808 static const int iResidLocY = this->GetIndex(myTrajectoryHists1D, "residLocY");
0809 myTrajectoryHists1D[iResidLocY]->Fill(residuum);
0810
0811 static const int iMeasLocYvsHit = this->GetIndex(myTrajectoryHists2D, "measLocYvsHit");
0812 myTrajectoryHists2D[iMeasLocYvsHit]->Fill(iRow / 2, measurementsLoc[iRow]);
0813 static const int iTrajLocYvsHit = this->GetIndex(myTrajectoryHists2D, "trajLocYvsHit");
0814 myTrajectoryHists2D[iTrajLocYvsHit]->Fill(iRow / 2, trajectoryLoc[iRow]);
0815 static const int iResidLocYvsHit = this->GetIndex(myTrajectoryHists2D, "residLocYvsHit");
0816 myTrajectoryHists2D[iResidLocYvsHit]->Fill(iRow / 2, residuum);
0817
0818 if (covMeasLocRow > 0.) {
0819 static const int iReduResidLocY = this->GetIndex(myTrajectoryHists1D, "reduResidLocY");
0820 myTrajectoryHists1D[iReduResidLocY]->Fill(residuum / TMath::Sqrt(covMeasLocRow));
0821 static const int iReduResidLocYvsHit = this->GetIndex(myTrajectoryHists2D, "reduResidLocYvsHit");
0822 myTrajectoryHists2D[iReduResidLocYvsHit]->Fill(iRow / 2, residuum / TMath::Sqrt(covMeasLocRow));
0823 }
0824 }
0825
0826 float nHitRow = iRow / 2;
0827 if (TMath::Odd(iRow))
0828 nHitRow += 0.5;
0829
0830 for (int iCol = iRow + 1; iCol < nRow; ++iCol) {
0831 double rho = TMath::Sqrt(covMeasLocRow + covMeasLoc[iCol][iCol]);
0832 rho = (0. == rho ? -2 : covMeasLoc[iRow][iCol] / rho);
0833 float nHitCol = iCol / 2;
0834 if (TMath::Odd(iCol))
0835 nHitCol += 0.5;
0836
0837 if (0. == rho)
0838 continue;
0839 static const int iProfileCorr = this->GetIndex(myTrajectoryHists2D, "profCorr");
0840 myTrajectoryHists2D[iProfileCorr]->Fill(nHitRow, nHitCol, TMath::Abs(rho));
0841 if (iRow + 1 == iCol && TMath::Even(iRow)) {
0842
0843
0844 static const int iHitCorrXy = this->GetIndex(myTrajectoryHists2D, "hitCorrXy");
0845 myTrajectoryHists2D[iHitCorrXy]->Fill(iRow / 2, rho);
0846 static const int iHitCorrXyLog = this->GetIndex(myTrajectoryHists2D, "hitCorrXyLog");
0847 myTrajectoryHists2D[iHitCorrXyLog]->Fill(iRow / 2, TMath::Abs(rho));
0848 if (is2DhitRow) {
0849 static const int iHitCorrXyValid = this->GetIndex(myTrajectoryHists2D, "hitCorrXyValid");
0850 myTrajectoryHists2D[iHitCorrXyValid]->Fill(iRow / 2, rho);
0851 static const int iHitCorrXyLogValid = this->GetIndex(myTrajectoryHists2D, "hitCorrXyLogValid");
0852 myTrajectoryHists2D[iHitCorrXyLogValid]->Fill(iRow / 2, TMath::Abs(rho));
0853 }
0854 } else {
0855 static const int iProfCorrOffXy = this->GetIndex(myTrajectoryHists2D, "profCorrOffXy");
0856 myTrajectoryHists2D[iProfCorrOffXy]->Fill(nHitRow, nHitCol, TMath::Abs(rho));
0857 static const int iHitCorrOffBlock = this->GetIndex(myTrajectoryHists2D, "hitCorrOffBlock");
0858 myTrajectoryHists2D[iHitCorrOffBlock]->Fill(nHitRow, rho);
0859 static const int iHitCorOffBlkLg = this->GetIndex(myTrajectoryHists2D, "hitCorrOffBlockLog");
0860 myTrajectoryHists2D[iHitCorOffBlkLg]->Fill(nHitRow, TMath::Abs(rho));
0861 }
0862 }
0863
0864
0865 for (int iCol = 0; iCol < derivatives.num_col(); ++iCol) {
0866 static const int iProfDerivatives = this->GetIndex(myTrajectoryHists2D, "profDerivatives");
0867 myTrajectoryHists2D[iProfDerivatives]->Fill(nHitRow, iCol, derivatives[iRow][iCol]);
0868 static const int iDerivatives = this->GetIndex(myTrajectoryHists2D, "derivatives");
0869 myTrajectoryHists2D[iDerivatives]->Fill(iCol, derivatives[iRow][iCol]);
0870 static const int iDerivativesLog = this->GetIndex(myTrajectoryHists2D, "derivativesLog");
0871 myTrajectoryHists2D[iDerivativesLog]->Fill(iCol, TMath::Abs(derivatives[iRow][iCol]));
0872 static const int iDerivativesPhi = this->GetIndex(myTrajectoryHists2D, "derivativesVsPhi");
0873 myTrajectoryHists2D[iDerivativesPhi]->Fill(recHits[iRow / 2]->det()->position().phi(), derivatives[iRow][iCol]);
0874 }
0875 }
0876 }
0877
0878
0879 void MillePedeMonitor::fillDerivatives(const ConstRecHitPointer &recHit,
0880 const float *localDerivs,
0881 unsigned int nLocal,
0882 const float *globalDerivs,
0883 unsigned int nGlobal,
0884 const int *labels) {
0885 const double phi = recHit->det()->position().phi();
0886
0887 static const int iLocPar = this->GetIndex(myDerivHists2D, "localDerivsPar");
0888 static const int iLocPhi = this->GetIndex(myDerivHists2D, "localDerivsPhi");
0889 static const int iLocParLog = this->GetIndex(myDerivHists2D, "localDerivsParLog");
0890 static const int iLocPhiLog = this->GetIndex(myDerivHists2D, "localDerivsPhiLog");
0891 for (unsigned int i = 0; i < nLocal; ++i) {
0892 myDerivHists2D[iLocPar]->Fill(i, localDerivs[i]);
0893 myDerivHists2D[iLocPhi]->Fill(phi, localDerivs[i]);
0894 if (localDerivs[i]) {
0895 myDerivHists2D[iLocParLog]->Fill(i, TMath::Abs(localDerivs[i]));
0896 myDerivHists2D[iLocPhiLog]->Fill(phi, TMath::Abs(localDerivs[i]));
0897 }
0898 }
0899
0900 static const int iGlobPar = this->GetIndex(myDerivHists2D, "globalDerivsPar");
0901 static const int iGlobPhi = this->GetIndex(myDerivHists2D, "globalDerivsPhi");
0902 static const int iGlobParLog = this->GetIndex(myDerivHists2D, "globalDerivsParLog");
0903 static const int iGlobPhiLog = this->GetIndex(myDerivHists2D, "globalDerivsPhiLog");
0904
0905 for (unsigned int i = 0; i < nGlobal; ++i) {
0906 const unsigned int parNum = (labels ? (labels[i] % PedeLabelerBase::theMaxNumParam) - 1 : i);
0907 myDerivHists2D[iGlobPar]->Fill(parNum, globalDerivs[i]);
0908 myDerivHists2D[iGlobPhi]->Fill(phi, globalDerivs[i]);
0909 if (globalDerivs[i]) {
0910 myDerivHists2D[iGlobParLog]->Fill(parNum, TMath::Abs(globalDerivs[i]));
0911 myDerivHists2D[iGlobPhiLog]->Fill(phi, TMath::Abs(globalDerivs[i]));
0912
0913 }
0914 }
0915 }
0916
0917
0918 void MillePedeMonitor::fillResiduals(const ConstRecHitPointer &recHit,
0919 const TrajectoryStateOnSurface &tsos,
0920 unsigned int nHit,
0921 float residuum,
0922 float sigma,
0923 bool isY) {
0924
0925
0926 const GlobalPoint detPos(recHit->det()->position());
0927 const double phi = detPos.phi();
0928
0929
0930 static const int iResPhi = this->GetIndex(myResidHists2D, "residPhi");
0931 myResidHists2D[iResPhi]->Fill(phi, residuum);
0932 static const int iSigPhi = this->GetIndex(myResidHists2D, "sigmaPhi");
0933 myResidHists2D[iSigPhi]->Fill(phi, sigma);
0934 if (sigma) {
0935 static const int iReduResPhi = this->GetIndex(myResidHists2D, "reduResidPhi");
0936 myResidHists2D[iReduResPhi]->Fill(phi, residuum / sigma);
0937 }
0938
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 const LocalVector mom(tsos.localDirection());
0972 const float phiSensToNorm = TMath::ATan(TMath::Abs((isY ? mom.y() : mom.x()) / mom.z()));
0973
0974 std::vector<std::vector<TH1 *> > &histArrayVec = (isY ? myResidHistsVec1DY : myResidHistsVec1DX);
0975 std::vector<TH1 *> &hitHists = (isY ? myResidHitHists1DY : myResidHitHists1DX);
0976
0977
0978 this->fillResidualHists(histArrayVec[0], phiSensToNorm, residuum, sigma);
0979 this->fillResidualHitHists(hitHists, phiSensToNorm, residuum, sigma, nHit);
0980
0981 const DetId detId(recHit->det()->geographicalId());
0982 if (detId.det() == DetId::Tracker) {
0983
0984 unsigned int subDetNum = detId.subdetId();
0985 if (subDetNum < histArrayVec.size() && subDetNum > 0) {
0986 this->fillResidualHists(histArrayVec[subDetNum], phiSensToNorm, residuum, sigma);
0987 } else {
0988 if (detId != AlignableBeamSpot::detId())
0989 edm::LogWarning("Alignment") << "@SUB=MillePedeMonitor::fillResiduals"
0990 << "Expect subDetNum from 1 to 6, got " << subDetNum;
0991 }
0992 }
0993 }
0994
0995
0996 void MillePedeMonitor::fillResidualHists(const std::vector<TH1 *> &hists,
0997 float phiSensToNorm,
0998 float residuum,
0999 float sigma) {
1000
1001
1002
1003
1004 static const int iRes = this->GetIndex(hists, "resid");
1005 hists[iRes]->Fill(residuum);
1006 static const int iSigma = this->GetIndex(hists, "sigma");
1007 hists[iSigma]->Fill(sigma);
1008 static const int iSigmaVsAngle = this->GetIndex(hists, "sigmaVsAngle");
1009 hists[iSigmaVsAngle]->Fill(phiSensToNorm, sigma);
1010 static const int iResidVsAngle = this->GetIndex(hists, "residVsAngle");
1011 hists[iResidVsAngle]->Fill(phiSensToNorm, residuum);
1012 static const int iReduRes = this->GetIndex(hists, "reduResid");
1013 static const int iReduResidVsAngle = this->GetIndex(hists, "reduResidVsAngle");
1014 if (sigma) {
1015 hists[iReduRes]->Fill(residuum / sigma);
1016 hists[iReduResidVsAngle]->Fill(phiSensToNorm, residuum / sigma);
1017 }
1018 static const int iAngle = this->GetIndex(hists, "angle");
1019 hists[iAngle]->Fill(phiSensToNorm);
1020
1021 if (phiSensToNorm > TMath::DegToRad() * 45.) {
1022 static const int iResGt45 = this->GetIndex(hists, "residGt45");
1023 hists[iResGt45]->Fill(residuum);
1024 static const int iSigmaGt45 = this->GetIndex(hists, "sigmaGt45");
1025 hists[iSigmaGt45]->Fill(sigma);
1026 static const int iReduResGt45 = this->GetIndex(hists, "reduResidGt45");
1027 if (sigma)
1028 hists[iReduResGt45]->Fill(residuum / sigma);
1029 } else {
1030 static const int iResLt45 = this->GetIndex(hists, "residLt45");
1031 hists[iResLt45]->Fill(residuum);
1032 static const int iSigmaLt45 = this->GetIndex(hists, "sigmaLt45");
1033 hists[iSigmaLt45]->Fill(sigma);
1034 static const int iReduResLt45 = this->GetIndex(hists, "reduResidLt45");
1035 if (sigma)
1036 hists[iReduResLt45]->Fill(residuum / sigma);
1037 }
1038 }
1039
1040
1041 void MillePedeMonitor::fillResidualHitHists(
1042 const std::vector<TH1 *> &hists, float angle, float residuum, float sigma, unsigned int nHit) {
1043
1044
1045
1046
1047 static constexpr unsigned int maxNhit = 29;
1048 static const auto iResHit = indexArray1D<TH1, maxNhit + 1>(hists, "resid_%d");
1049 static const auto iSigmaHit = indexArray1D<TH1, maxNhit + 1>(hists, "sigma_%d");
1050 static const auto iReduResHit = indexArray1D<TH1, maxNhit + 1>(hists, "reduResid_%d");
1051 static const auto iAngleHit = indexArray1D<TH1, maxNhit + 1>(hists, "angle_%d");
1052 if (nHit > maxNhit)
1053 nHit = maxNhit;
1054
1055 hists[iResHit[nHit]]->Fill(residuum);
1056 hists[iSigmaHit[nHit]]->Fill(sigma);
1057 if (sigma) {
1058 hists[iReduResHit[nHit]]->Fill(residuum / sigma);
1059 }
1060 hists[iAngleHit[nHit]]->Fill(angle);
1061 }
1062
1063
1064 void MillePedeMonitor::fillFrameToFrame(const AlignableDetOrUnitPtr &aliDet, const Alignable *ali) {
1065
1066 FrameToFrameDerivative ftfd;
1067 const AlgebraicMatrix frameDeriv = ftfd.frameToFrameDerivative(aliDet, ali);
1068
1069
1070 static const int iF2f = this->GetIndex(myFrame2FrameHists2D, "frame2frame");
1071 static const int iF2fAbs = this->GetIndex(myFrame2FrameHists2D, "frame2frameAbs");
1072 static const auto iF2fIjPhi = indexArray2D<TH2, 6>(myFrame2FrameHists2D, "frame2framePhi%d%d");
1073 static const auto iF2fIjPhiLog = indexArray2D<TH2, 6>(myFrame2FrameHists2D, "frame2framePhiLog%d%d");
1074 static const auto iF2fIjR = indexArray2D<TH2, 6>(myFrame2FrameHists2D, "frame2frameR%d%d");
1075 static const auto iF2fIjRLog = indexArray2D<TH2, 6>(myFrame2FrameHists2D, "frame2frameRLog%d%d");
1076
1077 const double phi = aliDet->globalPosition().phi();
1078 const double r = aliDet->globalPosition().perp();
1079 for (unsigned int i = 0; i < 6; ++i) {
1080 for (unsigned int j = 0; j < 6; ++j) {
1081 myFrame2FrameHists2D[iF2f]->Fill(i, j, frameDeriv[i][j]);
1082 myFrame2FrameHists2D[iF2fIjPhi[i][j]]->Fill(phi, frameDeriv[i][j]);
1083 myFrame2FrameHists2D[iF2fIjR[i][j]]->Fill(r, frameDeriv[i][j]);
1084 if (frameDeriv[i][j]) {
1085 myFrame2FrameHists2D[iF2fAbs]->Fill(i, j, TMath::Abs(frameDeriv[i][j]));
1086 myFrame2FrameHists2D[iF2fIjPhiLog[i][j]]->Fill(phi, TMath::Abs(frameDeriv[i][j]));
1087 myFrame2FrameHists2D[iF2fIjRLog[i][j]]->Fill(r, TMath::Abs(frameDeriv[i][j]));
1088 }
1089 }
1090 }
1091 }
1092
1093
1094 void MillePedeMonitor::fillCorrelations2D(float corr, const ConstRecHitPointer &recHit) {
1095 const DetId detId(recHit->det()->geographicalId());
1096 if (detId.det() != DetId::Tracker)
1097 return;
1098
1099 if ((detId.subdetId() < 1 || detId.subdetId() > 6) && detId != AlignableBeamSpot::detId()) {
1100 edm::LogWarning("Alignment") << "@SUB=MillePedeMonitor::fillCorrelations2D"
1101 << "Expect subdetId from 1 to 6, got " << detId.subdetId();
1102 return;
1103 }
1104
1105 myCorrHists[detId.subdetId() - 1]->Fill(corr);
1106 }
1107
1108
1109 void MillePedeMonitor::fillPxbSurveyHistsChi2(const float &chi2) {
1110 static const int iPxbSurvChi2_lo = this->GetIndex(myPxbSurveyHists, "PxbSurvChi2_lo");
1111 myPxbSurveyHists[iPxbSurvChi2_lo]->Fill(chi2);
1112 static const int iPxbSurvChi2_md = this->GetIndex(myPxbSurveyHists, "PxbSurvChi2_md");
1113 myPxbSurveyHists[iPxbSurvChi2_md]->Fill(chi2);
1114 static const int iPxbSurvChi2_hi = this->GetIndex(myPxbSurveyHists, "PxbSurvChi2_hi");
1115 myPxbSurveyHists[iPxbSurvChi2_hi]->Fill(chi2);
1116 static const int iPxbSurvChi2prob = this->GetIndex(myPxbSurveyHists, "PxbSurvChi2prob");
1117 myPxbSurveyHists[iPxbSurvChi2prob]->Fill(TMath::Prob(chi2, 4));
1118 }
1119
1120
1121 void MillePedeMonitor::fillPxbSurveyHistsLocalPars(const float &a0, const float &a1, const float &S, const float &phi) {
1122 static const int iPxbSurv_a0 = this->GetIndex(myPxbSurveyHists, "PxbSurv_a0");
1123 myPxbSurveyHists[iPxbSurv_a0]->Fill(a0);
1124 static const int iPxbSurv_a0_abs = this->GetIndex(myPxbSurveyHists, "PxbSurv_a0_abs");
1125 myPxbSurveyHists[iPxbSurv_a0_abs]->Fill(fabs(a0));
1126 static const int iPxbSurv_a1 = this->GetIndex(myPxbSurveyHists, "PxbSurv_a1");
1127 myPxbSurveyHists[iPxbSurv_a1]->Fill(a1);
1128 static const int iPxbSurv_a1_abs = this->GetIndex(myPxbSurveyHists, "PxbSurv_a1_abs");
1129 myPxbSurveyHists[iPxbSurv_a1_abs]->Fill(fabs(a1));
1130 static const int iPxbSurv_scale = this->GetIndex(myPxbSurveyHists, "PxbSurv_scale");
1131 myPxbSurveyHists[iPxbSurv_scale]->Fill(S);
1132 static const int iPxbSurv_phi = this->GetIndex(myPxbSurveyHists, "PxbSurv_phi");
1133 myPxbSurveyHists[iPxbSurv_phi]->Fill(phi);
1134 }