Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:56:37 

0001 /**
0002  * \file MillePedeMonitor.cc
0003  *
0004  *  \author    : Gero Flucke
0005  *  date       : October 2006
0006  *  $Revision: 1.24 $
0007  *  $Date: 2011/02/16 13:11:57 $
0008  *  (last update by $Author: mussgill $)
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   //  cout << "MillePedeMonitor using input TDirectory" << endl;
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;  //hists are deleted with their directory
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.};  // fully initialised with 0.
0079   if (!this->equidistLogBins(binsPt, kNumBins, 0.8, 100.)) {
0080     //     cerr << "MillePedeMonitor::init: problem with log bins" << endl;
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   // used track
0137   myUsedTrackHists1D = this->cloneHists(myTrackHists1D, "used", " (used tracks)");
0138   myUsedTrackHists2D = this->cloneHists(myTrackHists2D, "used", " (used tracks)");
0139   // must be after clone: index in vector!
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   // ReferenceTrajectory
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   // 2D vs. hit
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   //  myTrajectoryHists2D.back()->SetCanExtend(TH1::kAllAxes);
0239 
0240   TDirectory *dirTraject = directory->mkdir("refTrajectoryHists", "ReferenceTrajectory's");
0241   this->addToDirectory(myTrajectoryHists2D, dirTraject);
0242   this->addToDirectory(myTrajectoryHists1D, dirTraject);
0243 
0244   // derivatives hists
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   //   this->equidistLogBins(logBins.GetArray(), logBins.GetSize()-1, 1.e-40, 1.E-35);
0309   //   myDerivHists2D.push_back
0310   //     (new TH2F("globalDerivsPhiLog2",
0311   //               "global derivatives (#neq 0) vs. #phi(det);#phi(det);|#partial/#partial(param)|",
0312   //               51, -TMath::Pi(), TMath::Pi(), logBins.GetSize()-1, logBins.GetArray()));
0313 
0314   TDirectory *dirDerivs = directory->mkdir("derivatives", "derivatives etc.");
0315   this->addToDirectory(myDerivHists2D, dirDerivs);
0316 
0317   // residual hists
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   //   myResidHists2D.push_back(new TProfile2D("residXProfXy",
0332   //                                        "mean |residuum| (u);x [cm];y [cm];#LT|residuum|#GT",
0333   //                                        25, -110., 110., 25, -110., 110.));
0334   //   myResidHists2D.push_back(new TProfile2D("residXProfZr",
0335   //                                        "mean |residuum| (u);z [cm];r_{#pm} [cm];#LT|residuum|#GT",
0336   //                                        25, -275., 275., 25, -110., 110.));
0337   //   myResidHists2D.push_back(new TProfile2D("residYProfXy",
0338   //                                        "mean |residuum| (v);x [cm];y [cm];#LT|residuum|#GT",
0339   //                                        25, -110., 110., 25, -110., 110.));
0340   //   myResidHists2D.push_back(new TProfile2D("residYProfZr",
0341   //                                        "mean |residuum| (v);z [cm];r_{#pm} [cm];#LT|residuum|#GT",
0342   //                                        25, -275., 275., 25, -110., 110.));
0343   //   myResidHists2D.push_back(new TProfile2D("reduResidXProfXy",
0344   //                                        "mean |residuum/#sigma| (u);x [cm];y [cm];#LT|res./#sigma|#GT",
0345   //                                        25, -110., 110., 25, -110., 110.));
0346   //   myResidHists2D.push_back(new TProfile2D("reduResidXProfZr",
0347   //                                        "mean |residuum/#sigma| (u);z [cm];r_{#pm} [cm];#LT|res./#sigma|#GT",
0348   //                                        25, -275., 275., 25, -110., 110.));
0349   //   myResidHists2D.push_back(new TProfile2D("reduResidYProfXy",
0350   //                                        "mean |residuum/#sigma| (v);x [cm];y [cm];#LT|res./#sigma|#GT",
0351   //                                        25, -110., 110., 25, -110., 110.));
0352   //   myResidHists2D.push_back(new TProfile2D("reduResidYProfZr",
0353   //                                        "mean |residuum/#sigma| (v);z [cm];r_{#pm} [cm];#LT|res./#sigma|#GT",
0354   //                                        25, -275., 275., 25, -110., 110.));
0355   //   myResidHists2D.push_back(new TProfile2D("sigmaXProfXy",
0356   //                                        "mean sigma (u);x [cm];y [cm];#LT#sigma#GT",
0357   //                                        25, -110., 110., 25, -110., 110.));
0358   //   myResidHists2D.push_back(new TProfile2D("sigmaXProfZr",
0359   //                                        "mean sigma (u);z [cm];r_{#pm} [cm];#LT#sigma#GT",
0360   //                                        25, -275., 275., 25, -110., 110.));
0361   //   myResidHists2D.push_back(new TProfile2D("sigmaYProfXy",
0362   //                                        "mean sigma (v);x [cm];y [cm];#LT#sigma#GT",
0363   //                                        25, -110., 110., 25, -110., 110.));
0364   //   myResidHists2D.push_back(new TProfile2D("sigmaYProfZr",
0365   //                                        "mean sigma (v);z [cm];r_{#pm} [cm];#LT#sigma#GT",
0366   //                                        25, -275., 275., 25, -110., 110.));
0367 
0368   TDirectory *dirResid = directory->mkdir("residuals", "hit residuals, sigma,...");
0369   this->addToDirectory(myResidHists2D, dirResid);
0370 
0371   // Residuum, hit sigma and res./sigma for all sensor/track angles and separated for large and
0372   // small angles with respect to the sensor normal in sensitive direction.
0373   // Here for x-measurements:
0374   std::vector<TH1 *> allResidHistsX;
0375   allResidHistsX.push_back(new TH1F("resid", "hit residuals;residuum [cm]", 101, -.5, .5));  //51,-.05, .05));
0376   //allResidHistsX.back()->SetCanExtend(TH1::kAllAxes);
0377   allResidHistsX.push_back(new TH1F("sigma", "hit uncertainties;#sigma [cm]", 100, 0., 1.));  //50, 0., .02));
0378   //allResidHistsX.back()->SetCanExtend(TH1::kAllAxes);
0379   allResidHistsX.push_back(
0380       new TH1F("reduResid", "reduced hit residuals;res./#sigma", 101, -10., 10.));  //51, -3., 3.));
0381   //  allResidHistsX.back()->SetCanExtend(TH1::kAllAxes);
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));  //51, -.05, .05));
0412   // allResidHistsX.back()->SetCanExtend(TH1::kAllAxes);
0413   allResidHistsX.push_back(
0414       new TH1F("sigmaGt45", "hit uncertainties(#phi_{n}^{sens}>45#circ);#sigma [cm]", 100, 0., 1.));  //50, 0., .02));
0415   // allResidHistsX.back()->SetCanExtend(TH1::kAllAxes);
0416   allResidHistsX.push_back(new TH1F(
0417       "reduResidGt45", "reduced hit residuals(#phi_{n}^{sens}>45#circ);res./#sigma", 101, -10., 10.));  //51,-3.,3.));
0418   // allResidHistsX.back()->SetCanExtend(TH1::kAllAxes);
0419   allResidHistsX.push_back(
0420       new TH1F("residLt45", "hit residuals (#phi_{n}^{sens}<45#circ);residuum [cm]", 101, -.5, .5));  //51, -.15, .15));
0421   // allResidHistsX.back()->SetCanExtend(TH1::kAllAxes);
0422   allResidHistsX.push_back(
0423       new TH1F("sigmaLt45", "hit uncertainties(#phi_{n}^{sens}<45#circ);#sigma [cm]", 100, 0., 1.));  //50, 0., .01));
0424   // allResidHistsX.back()->SetCanExtend(TH1::kAllAxes);
0425   allResidHistsX.push_back(new TH1F(
0426       "reduResidLt45", "reduced hit residuals(#phi_{n}^{sens}<45#circ);res./#sigma", 101, -10., 10.));  //51,-3.,3.));
0427   myResidHistsVec1DX.push_back(allResidHistsX);  // at [0] for all subdets together...
0428   // ... then separately - indices/order like DetId.subdetId() in tracker:
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   // finally, differential in hit number (but subdet independent)
0436   for (unsigned int iHit = 0; iHit < 30; ++iHit) {  // 4: for each hit fill only angle independent plots
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   // Now clone the same as above for y-ccordinate:
0454   myResidHistsVec1DY.push_back(this->cloneHists(allResidHistsX, "", " y-coord."));  // all subdets
0455   std::vector<TH1 *> &yHists1D = allResidHistsX;  //myResidHistsVec1DY.back(); crashes? ROOT?
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.");  // diff. in nHit
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   // farme-to-frame derivatives
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   // PXB Survey
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   // Filling 'bins' with borders of 'nBins' bins between 'first' and 'last'
0553   // that are equidistant when viewed in log scale,
0554   // so 'bins' must have length nBins+1;
0555   // If 'first', 'last' or 'nBins' are not positive, failure is reported.
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   // these hist exist only for 'used track' hists:
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;  // only real hits count as in track->numberOfValidHits()
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   }  // end loop on hits
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   // CHK
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()  // FIXME: as in MillePedeAlignmentAlgorithm::is2D()
0779                              || recHits[iRow / 2]->detUnit()->type().isTrackerPixel());
0780     //the GeomDetUnit* is zero for composite hits (matched hits in the tracker,...).
0781     if (TMath::Even(iRow)) {  // local x
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) {  // local y, 2D detectors only
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;  // '/2' not '/2.'!
0827     if (TMath::Odd(iRow))
0828       nHitRow += 0.5;  // y-hit gets 0.5
0829     // correlations
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;  //cf. comment nHitRow
0834       if (TMath::Odd(iCol))
0835         nHitCol += 0.5;  // dito
0836       //      myProfileCorr->Fill(nHitRow, nHitCol, TMath::Abs(rho));
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)) {  // i.e. if iRow is x and iCol the same hit's y
0842         //      static const int iProfileCorrXy = this->GetIndex(myTrajectoryHists2D,"profCorrOffXy");
0843         //  myTrajectoryHists2D[iProfileCorrOffXy]->Fill(iRow/2, rho);
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);  // nhitRow??
0851           static const int iHitCorrXyLogValid = this->GetIndex(myTrajectoryHists2D, "hitCorrXyLogValid");
0852           myTrajectoryHists2D[iHitCorrXyLogValid]->Fill(iRow / 2, TMath::Abs(rho));  // nhitRow??
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     }  // end loop on columns of covariance
0863 
0864     // derivatives
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   }  // end loop on rows of covarianvce
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   //   static const int iGlobPhiLog2 = this->GetIndex(myDerivHists2D, "globalDerivsPhiLog2");
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       //       myDerivHists2D[iGlobPhiLog2]->Fill(phi, TMath::Abs(globalDerivs[i]));
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   // isY == false: x-measurements
0925   // isY == true:  y-measurements
0926   const GlobalPoint detPos(recHit->det()->position());
0927   const double phi = detPos.phi();
0928   //  const double rSigned = (detPos.y() > 0. ? detPos.perp() : -detPos.perp());
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   //   if (isY) {
0940   //     static const int iResYXy = this->GetIndex(myResidHists2D, "residYProfXy");
0941   //     myResidHists2D[iResYXy]->Fill(detPos.x(), detPos.y(), TMath::Abs(residuum));
0942   //     static const int iResYZr = this->GetIndex(myResidHists2D, "residYProfZr");
0943   //     myResidHists2D[iResYZr]->Fill(detPos.z(), rSigned, TMath::Abs(residuum));
0944   //     static const int iSigmaYXy = this->GetIndex(myResidHists2D, "sigmaYProfXy");
0945   //     myResidHists2D[iSigmaYXy]->Fill(detPos.x(), detPos.y(), sigma);
0946   //     static const int iSigmaYZr = this->GetIndex(myResidHists2D, "sigmaYProfZr");
0947   //     myResidHists2D[iSigmaYZr]->Fill(detPos.z(), rSigned, sigma);
0948   //     if (sigma) {
0949   //       static const int iReduResYXy = this->GetIndex(myResidHists2D, "reduResidYProfXy");
0950   //       myResidHists2D[iReduResYXy]->Fill(detPos.x(), detPos.y(), TMath::Abs(residuum/sigma));
0951   //       static const int iReduResYZr = this->GetIndex(myResidHists2D, "reduResidYProfZr");
0952   //       myResidHists2D[iReduResYZr]->Fill(detPos.z(), rSigned, TMath::Abs(residuum/sigma));
0953   //     }
0954   //   } else {
0955   //     static const int iResXXy = this->GetIndex(myResidHists2D, "residXProfXy");
0956   //     myResidHists2D[iResXXy]->Fill(detPos.x(), detPos.y(), TMath::Abs(residuum));
0957   //     static const int iResXZr = this->GetIndex(myResidHists2D, "residXProfZr");
0958   //     myResidHists2D[iResXZr]->Fill(detPos.z(), rSigned, TMath::Abs(residuum));
0959   //     static const int iSigmaXXy = this->GetIndex(myResidHists2D, "sigmaXProfXy");
0960   //     myResidHists2D[iSigmaXXy]->Fill(detPos.x(), detPos.y(), sigma);
0961   //     static const int iSigmaXZr = this->GetIndex(myResidHists2D, "sigmaXProfZr");
0962   //     myResidHists2D[iSigmaXZr]->Fill(detPos.z(), rSigned, sigma);
0963   //     if (sigma) {
0964   //       static const int iReduResXXy = this->GetIndex(myResidHists2D, "reduResidXProfXy");
0965   //       myResidHists2D[iReduResXXy]->Fill(detPos.x(), detPos.y(), TMath::Abs(residuum/sigma));
0966   //       static const int iReduResXZr = this->GetIndex(myResidHists2D, "reduResidXProfZr");
0967   //       myResidHists2D[iReduResXZr]->Fill(detPos.z(), rSigned, TMath::Abs(residuum/sigma));
0968   //     }
0969   //   }
0970 
0971   const LocalVector mom(tsos.localDirection());  // mom.z()==0. impossible for TSOS:
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   // call with histArrayVec[0] first: 'static' inside is referring to "subdet-less" names (X/Y irrelevant)
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     //   const GeomDet::SubDetector subDetNum = recHit->det()->subDetector();
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   // Histogram indices are calculated at first call, so make sure the vector of hists at the first
1001   // call has the correct names (i.e. without subdet extension) and all later calls have the
1002   // same order of hists...
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   // Histogram indices are calculated at first call, so make sure the vector of hists at the first
1044   // call has the correct names (i.e. without subdet extension) and all later calls have the
1045   // same order of hists...
1046 
1047   static constexpr unsigned int maxNhit = 29;  // 0...29 foreseen in initialisation...
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;  // limit of hists
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   // get derivative of higher level structure w.r.t. det
1066   FrameToFrameDerivative ftfd;
1067   const AlgebraicMatrix frameDeriv = ftfd.frameToFrameDerivative(aliDet, ali);
1068   //const AlgebraicMatrix frameDeriv = ftfd.frameToFrameDerivativeAtOrgRot(aliDet, ali);
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();  // after misalignment...
1078   const double r = aliDet->globalPosition().perp();   // after misalignment...
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!