Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-02-14 12:45:17 

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