CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
MillePedeMonitor.cc
Go to the documentation of this file.
1 
11 #include "DataFormats/GeometrySurface/interface/Surface.h"
12 #include "Alignment/MillePedeAlignmentAlgorithm/interface/MillePedeMonitor.h"
13 
14 #include "DataFormats/TrackReco/interface/Track.h"
15 
16 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
17 
18 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
19 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
20 #include <Geometry/CommonDetUnit/interface/GeomDetType.h>
21 
22 #include "Alignment/CommonAlignment/interface/Alignable.h"
23 #include "Alignment/CommonAlignment/interface/AlignableBeamSpot.h"
24 #include "Alignment/CommonAlignment/interface/AlignableDetOrUnitPtr.h"
25 #include "Alignment/CommonAlignmentParametrization/interface/FrameToFrameDerivative.h"
26 
27 #include "DataFormats/DetId/interface/DetId.h"
28 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
29 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
30 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
31 #include "Geometry/Records/interface/IdealGeometryRecord.h"
32 const int kBPIX = PixelSubdetector::PixelBarrel;
33 const int kFPIX = PixelSubdetector::PixelEndcap;
34 
35 #include <TProfile2D.h>
36 #include <TFile.h>
37 #include <TDirectory.h>
38 #include <TMath.h>
39 
40 #include "Alignment/MillePedeAlignmentAlgorithm/interface/PedeLabelerBase.h"
41 
42 typedef TransientTrackingRecHit::ConstRecHitPointer ConstRecHitPointer;
43 
44 //__________________________________________________________________
45 MillePedeMonitor::MillePedeMonitor(const TrackerTopology *tTopo, const char *rootFileName)
46  : myRootDir(nullptr), myDeleteDir(false), trackerTopology(tTopo) {
47  myRootDir = TFile::Open(rootFileName, "recreate");
48  myDeleteDir = true;
49 
50  this->init(myRootDir);
51 }
52 
53 //__________________________________________________________________
54 MillePedeMonitor::MillePedeMonitor(TDirectory *rootDir, const TrackerTopology *tTopo)
55  : myRootDir(nullptr), myDeleteDir(false), trackerTopology(tTopo) {
56  // cout << "MillePedeMonitor using input TDirectory" << endl;
57 
58  myRootDir = rootDir;
59  myDeleteDir = false;
60 
61  this->init(myRootDir);
62 }
63 
64 //__________________________________________________________________
65 MillePedeMonitor::~MillePedeMonitor() {
66  myRootDir->Write();
67  if (myDeleteDir)
68  delete myRootDir; //hists are deleted with their directory
69 }
70 
71 //__________________________________________________________________
72 bool MillePedeMonitor::init(TDirectory *directory) {
73  if (!directory)
74  return false;
75  TDirectory *oldDir = gDirectory;
76 
77  const int kNumBins = 20;
78  double binsPt[kNumBins + 1] = {0.}; // fully initialised with 0.
79  if (!this->equidistLogBins(binsPt, kNumBins, 0.8, 100.)) {
80  // cerr << "MillePedeMonitor::init: problem with log bins" << endl;
81  }
82  const int nHits = 35;
83 
84  myTrackHists1D.push_back(new TH1F("ptTrackLogBins", "p_{t}(track);p_{t} [GeV]", kNumBins, binsPt));
85 
86  myTrackHists1D.push_back(new TH1F("ptTrack", "p_{t}(track);p_{t} [GeV]", kNumBins, binsPt[0], binsPt[kNumBins]));
87  myTrackHists1D.push_back(new TH1F("pTrack", "p(track);p [GeV]", kNumBins, binsPt[0], 1.3 * binsPt[kNumBins]));
88  myTrackHists1D.push_back(new TH1F("etaTrack", "#eta(track);#eta", 26, -2.6, 2.6));
89  myTrackHists1D.push_back(new TH1F("thetaTrack", "#theta(track);#theta", 100, 0., TMath::Pi()));
90  myTrackHists1D.push_back(new TH1F("phiTrack", "#phi(track);#phi", 15, -TMath::Pi(), TMath::Pi()));
91 
92  myTrackHists1D.push_back(new TH1F("nHitTrack", "N_{hit}(track);N_{hit}", 40, 5., 45.));
93  myTrackHists1D.push_back(new TH1F("nHitBPIXTrack", "N_{hit, BPIX}(track);N_{hit}", 45, 0., 45.));
94  myTrackHists1D.push_back(new TH1F("nHitFPIXTrack", "N_{hit, FPIX}(track);N_{hit}", 45, 0., 45.));
95  myTrackHists1D.push_back(new TH1F("nHitFPIXplusTrack", "N_{hit, BPIXplus}(track);N_{hit}", 45, 0., 45.));
96  myTrackHists1D.push_back(new TH1F("nHitFPIXminusTrack", "N_{hit, BPIXminus}(track);N_{hit}", 45, 0., 45.));
97  myTrackHists1D.push_back(new TH1F("nHitPIXELTrack", "N_{hit, PIXEL}(track);N_{hit}", 45, 0., 45.));
98  myTrackHists1D.push_back(new TH1F("nHitTIBTrack", "N_{hit, TIB}(track);N_{hit}", 45, 0., 45.));
99  myTrackHists1D.push_back(new TH1F("nHitTOBTrack", "N_{hit, TOB}(track);N_{hit}", 45, 0., 45.));
100  myTrackHists1D.push_back(new TH1F("nHitTIDplusTrack", "N_{hit, TIDplus}(track);N_{hit}", 45, 0., 45.));
101  myTrackHists1D.push_back(new TH1F("nHitTIDminusTrack", "N_{hit, TIDminus}(track);N_{hit}", 45, 0., 45.));
102  myTrackHists1D.push_back(new TH1F("nHitTIDTrack", "N_{hit, TID}(track);N_{hit}", 45, 0., 45.));
103  myTrackHists1D.push_back(new TH1F("nHitTECplusTrack", "N_{hit, TECplus}(track);N_{hit}", 45, 0., 45.));
104  myTrackHists1D.push_back(new TH1F("nHitTECminusTrack", "N_{hit, TECminus}(track);N_{hit}", 45, 0., 45.));
105  myTrackHists1D.push_back(new TH1F("nHitTECTrack", "N_{hit, TEC}(track);N_{hit}", 45, 0., 45.));
106  myTrackHists1D.push_back(new TH1F("nHitENDCAPplusTrack", "N_{hit, ENDCAPplus}(track);N_{hit}", 45, 0., 45.));
107  myTrackHists1D.push_back(new TH1F("nHitENDCAPminusTrack", "N_{hit, ENDCAPminus}(track);N_{hit}", 45, 0., 45.));
108  myTrackHists1D.push_back(new TH1F("nHitENDCAPTrack", "N_{hit, ENDCAP}(track);N_{hit}", 45, 0., 45.));
109  myTrackHists2D.push_back(new TH2F(
110  "nHitENDCAPTrackMinusVsPlus", "N_{hit, ENDCAP}(track);N_{hit, plus};N_{hit, minus}", 45, 0., 45., 45, 0., 45.));
111  myTrackHists1D.push_back(new TH1F("nHitInvalidTrack", "N_{hit, invalid}(track)", 5, 0., 5.));
112  myTrackHists1D.push_back(new TH1F("r1Track", "r(1st hit);r [cm]", 20, 0., 20.));
113  myTrackHists1D.push_back(new TH1F("phi1Track", "#phi(1st hit);#phi", 30, -TMath::Pi(), TMath::Pi()));
114  myTrackHists1D.push_back(new TH1F("y1Track", "y(1st hit);y [cm]", 40, -120., +120.));
115  myTrackHists1D.push_back(new TH1F("z1Track", "z(1st hit);z [cm]", 20, -50, +50));
116  myTrackHists1D.push_back(new TH1F("r1TrackSigned", "r(1st hit);r_{#pm} [cm]", 40, -120., 120.));
117  myTrackHists1D.push_back(new TH1F("z1TrackFull", "z(1st hit);z [cm]", 20, -300., +300.));
118  myTrackHists1D.push_back(new TH1F("rLastTrack", "r(last hit);r [cm]", 20, 20., 120.));
119  myTrackHists1D.push_back(new TH1F("phiLastTrack", "#phi(last hit);#phi", 30, -TMath::Pi(), TMath::Pi()));
120  myTrackHists1D.push_back(new TH1F("yLastTrack", "y(last hit);y [cm]", 40, -120., +120.));
121  myTrackHists1D.push_back(new TH1F("zLastTrack", "z(last hit);z [cm]", 30, -300., +300.));
122  myTrackHists1D.push_back(new TH1F("chi2PerNdf", "#chi^{2}/ndf;#chi^{2}/ndf", 500, 0., 50.));
123  myTrackHists1D.push_back(new TH1F("impParZ", "impact parameter in z", 20, -20., 20.));
124  myTrackHists1D.push_back(new TH1F("impParErrZ", "error of impact parameter in z", 40, 0., 0.06));
125  myTrackHists1D.push_back(new TH1F("impParRphi", "impact parameter in r#phi", 51, -0.05, .05));
126  myTrackHists1D.push_back(new TH1F("impParErrRphi", "error of impact parameter in r#phi", 50, 0., 0.01));
127 
128  myTrackHists2D.push_back(
129  new TH2F("rz1TrackFull", "rz(1st hit);z [cm]; r_{#pm} [cm]", 40, -282., +282., 40, -115., +115.));
130  myTrackHists2D.push_back(new TH2F("xy1Track", "xy(1st hit);x [cm]; y [cm]", 40, -115., +115., 40, -115., +115.));
131 
132  TDirectory *dirTracks = directory->mkdir("trackHists", "input tracks");
133  this->addToDirectory(myTrackHists1D, dirTracks);
134  this->addToDirectory(myTrackHists2D, dirTracks);
135 
136  // used track
137  myUsedTrackHists1D = this->cloneHists(myTrackHists1D, "used", " (used tracks)");
138  myUsedTrackHists2D = this->cloneHists(myTrackHists2D, "used", " (used tracks)");
139  // must be after clone: index in vector!
140  myUsedTrackHists1D.push_back(new TH1F("usedHitsX", "n(x-hits) transferred to pede;n(x-hits)", nHits, 0., nHits));
141  myUsedTrackHists1D.push_back(new TH1F("usedHitsY", "n(y-hits) transferred to pede;n(y-hits)", 10, 0., 10));
142 
143  TDirectory *dirUsedTracks = directory->mkdir("usedTrackHists", "used tracks");
144  this->addToDirectory(myUsedTrackHists1D, dirUsedTracks);
145  this->addToDirectory(myUsedTrackHists2D, dirUsedTracks);
146 
147  // ReferenceTrajectory
148  myTrajectoryHists1D.push_back(new TH1F("validRefTraj", "validity of ReferenceTrajectory", 2, 0., 2.));
149 
150  myTrajectoryHists2D.push_back(new TProfile2D(
151  "profCorr", "mean of |#rho|, #rho#neq0;hit x/y;hit x/y;", 2 * nHits, 0., nHits, 2 * nHits, 0., nHits));
152  myTrajectoryHists2D.push_back(new TProfile2D("profCorrOffXy",
153  "mean of |#rho|, #rho#neq0, no xy_{hit};hit x/y;hit x/y;",
154  2 * nHits,
155  0.,
156  nHits,
157  2 * nHits,
158  0.,
159  nHits));
160 
161  myTrajectoryHists2D.push_back(new TH2F(
162  "hitCorrOffBlock", "hit correlations: off-block-diagonals;N(hit);#rho", 2 * nHits, 0., nHits, 81, -.06, .06));
163  TArrayD logBins(102);
164  this->equidistLogBins(logBins.GetArray(), logBins.GetSize() - 1, 1.E-11, .1);
165  myTrajectoryHists2D.push_back(new TH2F("hitCorrOffBlockLog",
166  "hit correlations: off-block-diagonals;N(hit);|#rho|",
167  2 * nHits,
168  0.,
169  nHits,
170  logBins.GetSize() - 1,
171  logBins.GetArray()));
172 
173  myTrajectoryHists2D.push_back(
174  new TH2F("hitCorrXy", "hit correlations: xy;N(hit);#rho", nHits, 0., nHits, 81, -.5, .5));
175  myTrajectoryHists2D.push_back(
176  new TH2F("hitCorrXyValid", "hit correlations: xy, 2D-det.;N(hit);#rho", nHits, 0., nHits, 81, -.02, .02));
177  this->equidistLogBins(logBins.GetArray(), logBins.GetSize() - 1, 1.E-10, .5);
178  myTrajectoryHists2D.push_back(new TH2F(
179  "hitCorrXyLog", "hit correlations: xy;N(hit);|#rho|", nHits, 0., nHits, logBins.GetSize() - 1, logBins.GetArray()));
180  myTrajectoryHists2D.push_back(new TH2F("hitCorrXyLogValid",
181  "hit correlations: xy, 2D-det.;N(hit);|#rho|",
182  nHits,
183  0.,
184  nHits,
185  logBins.GetSize() - 1,
186  logBins.GetArray()));
187 
188  myTrajectoryHists1D.push_back(new TH1F("measLocX", "local x measurements;x", 101, -6., 6.));
189  myTrajectoryHists1D.push_back(new TH1F("measLocY", "local y measurements, 2D-det.;y", 101, -10., 10.));
190  myTrajectoryHists1D.push_back(new TH1F("trajLocX", "local x trajectory;x", 101, -6., 6.));
191  myTrajectoryHists1D.push_back(new TH1F("trajLocY", "local y trajectory, 2D-det.;y", 101, -10., 10.));
192 
193  myTrajectoryHists1D.push_back(new TH1F("residLocX", "local x residual;#Deltax", 101, -.75, .75));
194  myTrajectoryHists1D.push_back(new TH1F("residLocY", "local y residual, 2D-det.;#Deltay", 101, -2., 2.));
195  myTrajectoryHists1D.push_back(new TH1F("reduResidLocX", "local x reduced residual;#Deltax/#sigma", 101, -20., 20.));
196  myTrajectoryHists1D.push_back(
197  new TH1F("reduResidLocY", "local y reduced residual, 2D-det.;#Deltay/#sigma", 101, -20., 20.));
198 
199  // 2D vs. hit
200  myTrajectoryHists2D.push_back(
201  new TH2F("measLocXvsHit", "local x measurements;hit;x", nHits, 0., nHits, 101, -6., 6.));
202  myTrajectoryHists2D.push_back(
203  new TH2F("measLocYvsHit", "local y measurements, 2D-det.;hit;y", nHits, 0., nHits, 101, -10., 10.));
204  myTrajectoryHists2D.push_back(new TH2F("trajLocXvsHit", "local x trajectory;hit;x", nHits, 0., nHits, 101, -6., 6.));
205  myTrajectoryHists2D.push_back(
206  new TH2F("trajLocYvsHit", "local y trajectory, 2D-det.;hit;y", nHits, 0., nHits, 101, -10., 10.));
207 
208  myTrajectoryHists2D.push_back(
209  new TH2F("residLocXvsHit", "local x residual;hit;#Deltax", nHits, 0., nHits, 101, -.75, .75));
210  myTrajectoryHists2D.push_back(
211  new TH2F("residLocYvsHit", "local y residual, 2D-det.;hit;#Deltay", nHits, 0., nHits, 101, -1., 1.));
212  myTrajectoryHists2D.push_back(
213  new TH2F("reduResidLocXvsHit", "local x reduced residual;hit;#Deltax/#sigma", nHits, 0., nHits, 101, -20., 20.));
214  myTrajectoryHists2D.push_back(new TH2F(
215  "reduResidLocYvsHit", "local y reduced residual, 2D-det.;hit;#Deltay/#sigma", nHits, 0., nHits, 101, -20., 20.));
216 
217  myTrajectoryHists2D.push_back(
218  new TProfile2D("profDerivatives", "mean derivatives;hit x/y;parameter;", 2 * nHits, 0., nHits, 10, 0., 10.));
219 
220  myTrajectoryHists2D.push_back(new TH2F(
221  "derivatives", "derivative;parameter;#partial(x/y)_{local}/#partial(param)", 10, 0., 10., 101, -20., 20.));
222  this->equidistLogBins(logBins.GetArray(), logBins.GetSize() - 1, 1.E-12, 100.);
223  myTrajectoryHists2D.push_back(new TH2F("derivativesLog",
224  "|derivative|;parameter;|#partial(x/y)_{local}/#partial(param)|",
225  10,
226  0.,
227  10.,
228  logBins.GetSize() - 1,
229  logBins.GetArray()));
230  myTrajectoryHists2D.push_back(new TH2F("derivativesVsPhi",
231  "derivatives vs. #phi;#phi(geomDet);#partial(x/y)_{local}/#partial(param)",
232  50,
233  -TMath::Pi(),
234  TMath::Pi(),
235  101,
236  -300.,
237  300.));
238  // myTrajectoryHists2D.back()->SetCanExtend(TH1::kAllAxes);
239 
240  TDirectory *dirTraject = directory->mkdir("refTrajectoryHists", "ReferenceTrajectory's");
241  this->addToDirectory(myTrajectoryHists2D, dirTraject);
242  this->addToDirectory(myTrajectoryHists1D, dirTraject);
243 
244  // derivatives hists
245  myDerivHists2D.push_back(new TH2F("localDerivsPar",
246  "local derivatives vs. paramater;parameter;#partial/#partial(param)",
247  6,
248  0.,
249  6.,
250  101,
251  -200.,
252  200.));
253  myDerivHists2D.push_back(new TH2F("localDerivsPhi",
254  "local derivatives vs. #phi(det);#phi(det);#partial/#partial(param)",
255  51,
256  -TMath::Pi(),
257  TMath::Pi(),
258  101,
259  -150.,
260  150.));
261  this->equidistLogBins(logBins.GetArray(), logBins.GetSize() - 1, 1.E-13, 150.);
262  myDerivHists2D.push_back(new TH2F("localDerivsParLog",
263  "local derivatives (#neq 0) vs. parameter;parameter;|#partial/#partial(param)|",
264  6,
265  0.,
266  6.,
267  logBins.GetSize() - 1,
268  logBins.GetArray()));
269  myDerivHists2D.push_back(new TH2F("localDerivsPhiLog",
270  "local derivatives (#neq 0) vs. #phi(det);#phi(det);|#partial/#partial(param)|",
271  51,
272  -TMath::Pi(),
273  TMath::Pi(),
274  logBins.GetSize() - 1,
275  logBins.GetArray()));
276  const unsigned int maxParNum = PedeLabelerBase::theMaxNumParam;
277  myDerivHists2D.push_back(new TH2F("globalDerivsPar",
278  "global derivatives vs. paramater;parameter;#partial/#partial(param)",
279  maxParNum,
280  0.,
281  maxParNum,
282  100,
283  -200,
284  200));
285  myDerivHists2D.push_back(new TH2F("globalDerivsPhi",
286  "global derivatives vs. #phi(det);#phi(det);#partial/#partial(param)",
287  51,
288  -TMath::Pi(),
289  TMath::Pi(),
290  102,
291  -1.02,
292  1.02));
293  this->equidistLogBins(logBins.GetArray(), logBins.GetSize() - 1, 1.E-7, 5.e2);
294  myDerivHists2D.push_back(new TH2F("globalDerivsParLog",
295  "global derivatives (#neq 0) vs. parameter;parameter;|#partial/#partial(param)|",
296  maxParNum,
297  0.,
298  maxParNum,
299  logBins.GetSize() - 1,
300  logBins.GetArray()));
301  myDerivHists2D.push_back(new TH2F("globalDerivsPhiLog",
302  "global derivatives (#neq 0) vs. #phi(det);#phi(det);|#partial/#partial(param)|",
303  51,
304  -TMath::Pi(),
305  TMath::Pi(),
306  logBins.GetSize() - 1,
307  logBins.GetArray()));
308  // this->equidistLogBins(logBins.GetArray(), logBins.GetSize()-1, 1.e-40, 1.E-35);
309  // myDerivHists2D.push_back
310  // (new TH2F("globalDerivsPhiLog2",
311  // "global derivatives (#neq 0) vs. #phi(det);#phi(det);|#partial/#partial(param)|",
312  // 51, -TMath::Pi(), TMath::Pi(), logBins.GetSize()-1, logBins.GetArray()));
313 
314  TDirectory *dirDerivs = directory->mkdir("derivatives", "derivatives etc.");
315  this->addToDirectory(myDerivHists2D, dirDerivs);
316 
317  // residual hists
318  myResidHists2D.push_back(new TH2F(
319  "residPhi", "residuum vs. #phi(det);#phi(det);residuum[cm]", 51, -TMath::Pi(), TMath::Pi(), 101, -.5, .5));
320  myResidHists2D.push_back(
321  new TH2F("sigmaPhi", "#sigma vs. #phi(det);#phi(det);#sigma[cm]", 51, -TMath::Pi(), TMath::Pi(), 101, .0, .1));
322  myResidHists2D.push_back(new TH2F("reduResidPhi",
323  "residuum/#sigma vs. #phi(det);#phi(det);residuum/#sigma",
324  51,
325  -TMath::Pi(),
326  TMath::Pi(),
327  101,
328  -14.,
329  14.));
330 
331  // myResidHists2D.push_back(new TProfile2D("residXProfXy",
332  // "mean |residuum| (u);x [cm];y [cm];#LT|residuum|#GT",
333  // 25, -110., 110., 25, -110., 110.));
334  // myResidHists2D.push_back(new TProfile2D("residXProfZr",
335  // "mean |residuum| (u);z [cm];r_{#pm} [cm];#LT|residuum|#GT",
336  // 25, -275., 275., 25, -110., 110.));
337  // myResidHists2D.push_back(new TProfile2D("residYProfXy",
338  // "mean |residuum| (v);x [cm];y [cm];#LT|residuum|#GT",
339  // 25, -110., 110., 25, -110., 110.));
340  // myResidHists2D.push_back(new TProfile2D("residYProfZr",
341  // "mean |residuum| (v);z [cm];r_{#pm} [cm];#LT|residuum|#GT",
342  // 25, -275., 275., 25, -110., 110.));
343  // myResidHists2D.push_back(new TProfile2D("reduResidXProfXy",
344  // "mean |residuum/#sigma| (u);x [cm];y [cm];#LT|res./#sigma|#GT",
345  // 25, -110., 110., 25, -110., 110.));
346  // myResidHists2D.push_back(new TProfile2D("reduResidXProfZr",
347  // "mean |residuum/#sigma| (u);z [cm];r_{#pm} [cm];#LT|res./#sigma|#GT",
348  // 25, -275., 275., 25, -110., 110.));
349  // myResidHists2D.push_back(new TProfile2D("reduResidYProfXy",
350  // "mean |residuum/#sigma| (v);x [cm];y [cm];#LT|res./#sigma|#GT",
351  // 25, -110., 110., 25, -110., 110.));
352  // myResidHists2D.push_back(new TProfile2D("reduResidYProfZr",
353  // "mean |residuum/#sigma| (v);z [cm];r_{#pm} [cm];#LT|res./#sigma|#GT",
354  // 25, -275., 275., 25, -110., 110.));
355  // myResidHists2D.push_back(new TProfile2D("sigmaXProfXy",
356  // "mean sigma (u);x [cm];y [cm];#LT#sigma#GT",
357  // 25, -110., 110., 25, -110., 110.));
358  // myResidHists2D.push_back(new TProfile2D("sigmaXProfZr",
359  // "mean sigma (u);z [cm];r_{#pm} [cm];#LT#sigma#GT",
360  // 25, -275., 275., 25, -110., 110.));
361  // myResidHists2D.push_back(new TProfile2D("sigmaYProfXy",
362  // "mean sigma (v);x [cm];y [cm];#LT#sigma#GT",
363  // 25, -110., 110., 25, -110., 110.));
364  // myResidHists2D.push_back(new TProfile2D("sigmaYProfZr",
365  // "mean sigma (v);z [cm];r_{#pm} [cm];#LT#sigma#GT",
366  // 25, -275., 275., 25, -110., 110.));
367 
368  TDirectory *dirResid = directory->mkdir("residuals", "hit residuals, sigma,...");
369  this->addToDirectory(myResidHists2D, dirResid);
370 
371  // Residuum, hit sigma and res./sigma for all sensor/track angles and separated for large and
372  // small angles with respect to the sensor normal in sensitive direction.
373  // Here for x-measurements:
374  std::vector<TH1 *> allResidHistsX;
375  allResidHistsX.push_back(new TH1F("resid", "hit residuals;residuum [cm]", 101, -.5, .5)); //51,-.05, .05));
376  //allResidHistsX.back()->SetCanExtend(TH1::kAllAxes);
377  allResidHistsX.push_back(new TH1F("sigma", "hit uncertainties;#sigma [cm]", 100, 0., 1.)); //50, 0., .02));
378  //allResidHistsX.back()->SetCanExtend(TH1::kAllAxes);
379  allResidHistsX.push_back(
380  new TH1F("reduResid", "reduced hit residuals;res./#sigma", 101, -10., 10.)); //51, -3., 3.));
381  // allResidHistsX.back()->SetCanExtend(TH1::kAllAxes);
382  allResidHistsX.push_back(
383  new TH1F("angle", "#phi_{tr} wrt normal (sens. plane);#phi_{n}^{sens}", 50, 0., TMath::PiOver2()));
384  allResidHistsX.push_back(new TH2F("residVsAngle",
385  "residuum vs. #phi_{tr} wrt normal (sens. plane);#phi_{n}^{sens};residuum [cm]",
386  50,
387  0.,
388  TMath::PiOver2(),
389  51,
390  -1.,
391  1.));
392  this->equidistLogBins(logBins.GetArray(), logBins.GetSize() - 1, 1.E-6, 100.);
393  allResidHistsX.push_back(new TH2F("sigmaVsAngle",
394  "#sigma vs. #phi_{tr} wrt normal (sens. plane);#phi_{n}^{sens};#sigma [cm]",
395  50,
396  0.,
397  TMath::PiOver2(),
398  logBins.GetSize() - 1,
399  logBins.GetArray()));
400  allResidHistsX.push_back(
401  new TH2F("reduResidVsAngle",
402  "reduced residuum vs. #phi_{tr} wrt normal (sens. plane);#phi_{n}^{sens};res./#sigma",
403  50,
404  0.,
405  TMath::PiOver2(),
406  51,
407  -15.,
408  15.));
409 
410  allResidHistsX.push_back(
411  new TH1F("residGt45", "hit residuals (#phi_{n}^{sens}>45#circ);residuum [cm]", 101, -.5, .5)); //51, -.05, .05));
412  // allResidHistsX.back()->SetCanExtend(TH1::kAllAxes);
413  allResidHistsX.push_back(
414  new TH1F("sigmaGt45", "hit uncertainties(#phi_{n}^{sens}>45#circ);#sigma [cm]", 100, 0., 1.)); //50, 0., .02));
415  // allResidHistsX.back()->SetCanExtend(TH1::kAllAxes);
416  allResidHistsX.push_back(new TH1F(
417  "reduResidGt45", "reduced hit residuals(#phi_{n}^{sens}>45#circ);res./#sigma", 101, -10., 10.)); //51,-3.,3.));
418  // allResidHistsX.back()->SetCanExtend(TH1::kAllAxes);
419  allResidHistsX.push_back(
420  new TH1F("residLt45", "hit residuals (#phi_{n}^{sens}<45#circ);residuum [cm]", 101, -.5, .5)); //51, -.15, .15));
421  // allResidHistsX.back()->SetCanExtend(TH1::kAllAxes);
422  allResidHistsX.push_back(
423  new TH1F("sigmaLt45", "hit uncertainties(#phi_{n}^{sens}<45#circ);#sigma [cm]", 100, 0., 1.)); //50, 0., .01));
424  // allResidHistsX.back()->SetCanExtend(TH1::kAllAxes);
425  allResidHistsX.push_back(new TH1F(
426  "reduResidLt45", "reduced hit residuals(#phi_{n}^{sens}<45#circ);res./#sigma", 101, -10., 10.)); //51,-3.,3.));
427  myResidHistsVec1DX.push_back(allResidHistsX); // at [0] for all subdets together...
428  // ... then separately - indices/order like DetId.subdetId() in tracker:
429  myResidHistsVec1DX.push_back(this->cloneHists(allResidHistsX, "TPB", " x-coord. in pixel barrel"));
430  myResidHistsVec1DX.push_back(this->cloneHists(allResidHistsX, "TPE", " x-coord. in pixel discs"));
431  myResidHistsVec1DX.push_back(this->cloneHists(allResidHistsX, "TIB", " x-coord. in TIB"));
432  myResidHistsVec1DX.push_back(this->cloneHists(allResidHistsX, "TID", " x-coord. in TID"));
433  myResidHistsVec1DX.push_back(this->cloneHists(allResidHistsX, "TOB", " x-coord. in TOB"));
434  myResidHistsVec1DX.push_back(this->cloneHists(allResidHistsX, "TEC", " x-coord. in TEC"));
435  // finally, differential in hit number (but subdet independent)
436  for (unsigned int iHit = 0; iHit < 30; ++iHit) { // 4: for each hit fill only angle independent plots
437  for (unsigned int iHist = 0; iHist < 4 && iHist < allResidHistsX.size(); ++iHist) {
438  TH1 *h = allResidHistsX[iHist];
439  myResidHitHists1DX.push_back(static_cast<TH1 *>(h->Clone(Form("%s_%d", h->GetName(), iHit))));
440  myResidHitHists1DX.back()->SetTitle(Form("%s, hit %d", h->GetTitle(), iHit));
441  }
442  }
443 
444  TDirectory *dirResidX = (dirResid ? dirResid : directory)->mkdir("X", "hit residuals etc. for x measurements");
445  this->addToDirectory(myResidHitHists1DX, dirResidX);
446  for (std::vector<std::vector<TH1 *> >::iterator vecIter = myResidHistsVec1DX.begin(),
447  vecIterEnd = myResidHistsVec1DX.end();
448  vecIter != vecIterEnd;
449  ++vecIter) {
450  this->addToDirectory(*vecIter, dirResidX);
451  }
452 
453  // Now clone the same as above for y-ccordinate:
454  myResidHistsVec1DY.push_back(this->cloneHists(allResidHistsX, "", " y-coord.")); // all subdets
455  std::vector<TH1 *> &yHists1D = allResidHistsX; //myResidHistsVec1DY.back(); crashes? ROOT?
456  myResidHistsVec1DY.push_back(this->cloneHists(yHists1D, "TPB", " y-coord. in pixel barrel"));
457  myResidHistsVec1DY.push_back(this->cloneHists(yHists1D, "TPE", " y-coord. in pixel discs"));
458  myResidHistsVec1DY.push_back(this->cloneHists(yHists1D, "TIB", " y-coord. in TIB"));
459  myResidHistsVec1DY.push_back(this->cloneHists(yHists1D, "TID", " y-coord. in TID"));
460  myResidHistsVec1DY.push_back(this->cloneHists(yHists1D, "TOB", " y-coord. in TOB"));
461  myResidHistsVec1DY.push_back(this->cloneHists(yHists1D, "TEC", " y-coord. in TEC"));
462  myResidHitHists1DY = this->cloneHists(myResidHitHists1DX, "", " y-coord."); // diff. in nHit
463 
464  TDirectory *dirResidY = (dirResid ? dirResid : directory)->mkdir("Y", "hit residuals etc. for y measurements");
465  this->addToDirectory(myResidHitHists1DY, dirResidY);
466  for (std::vector<std::vector<TH1 *> >::iterator vecIter = myResidHistsVec1DY.begin(),
467  vecIterEnd = myResidHistsVec1DY.end();
468  vecIter != vecIterEnd;
469  ++vecIter) {
470  this->addToDirectory(*vecIter, dirResidY);
471  }
472 
473  // farme-to-frame derivatives
474  myFrame2FrameHists2D.push_back(
475  new TProfile2D("frame2frame", "mean frame to frame derivatives;col;row", 6, 0., 6., 6, 0., 6.));
476  myFrame2FrameHists2D.push_back(
477  new TProfile2D("frame2frameAbs", "mean |frame to frame derivatives|, #neq0;col;row", 6, 0., 6., 6, 0., 6.));
478 
479  this->equidistLogBins(logBins.GetArray(), logBins.GetSize() - 1, 1.E-36, 100.);
480  for (unsigned int i = 0; i < 6; ++i) {
481  for (unsigned int j = 0; j < 6; ++j) {
482  myFrame2FrameHists2D.push_back(new TH2F(Form("frame2framePhi%d%d", i, j),
483  Form("frame to frame derivatives, %d%d;#phi(aliDet);deriv", i, j),
484  51,
485  -TMath::Pi(),
486  TMath::Pi(),
487  10,
488  0.,
489  1.));
490  myFrame2FrameHists2D.back()->SetCanExtend(TH1::kAllAxes);
491  myFrame2FrameHists2D.push_back(new TH2F(Form("frame2frameR%d%d", i, j),
492  Form("frame to frame derivatives, %d%d;r(aliDet);deriv", i, j),
493  51,
494  0.,
495  110.,
496  10,
497  0.,
498  1.));
499  myFrame2FrameHists2D.back()->SetCanExtend(TH1::kAllAxes);
500 
501  myFrame2FrameHists2D.push_back(
502  new TH2F(Form("frame2framePhiLog%d%d", i, j),
503  Form("frame to frame |derivatives|, %d%d, #neq0;#phi(aliDet);deriv", i, j),
504  51,
505  -TMath::Pi(),
506  TMath::Pi(),
507  logBins.GetSize() - 1,
508  logBins.GetArray()));
509  myFrame2FrameHists2D.push_back(new TH2F(Form("frame2frameRLog%d%d", i, j),
510  Form("frame to frame |derivatives|, %d%d, #neq0;r(aliDet);deriv", i, j),
511  51,
512  0.,
513  110.,
514  logBins.GetSize() - 1,
515  logBins.GetArray()));
516  }
517  }
518 
519  TDirectory *dirF2f = directory->mkdir("frame2FrameHists", "derivatives etc.");
520  this->addToDirectory(myFrame2FrameHists2D, dirF2f);
521 
522  myCorrHists.push_back(new TH1F("xyCorrTPB", "#rho_{xy} in pixel barrel", 50, -.5, .5));
523  myCorrHists.push_back(new TH1F("xyCorrTPE", "#rho_{xy} in forward pixel", 50, -.5, .5));
524  myCorrHists.push_back(new TH1F("xyCorrTIB", "#rho_{xy} in TIB", 50, -.5, .5));
525  myCorrHists.push_back(new TH1F("xyCorrTID", "#rho_{xy} in TID", 50, -1., 1.));
526  myCorrHists.push_back(new TH1F("xyCorrTOB", "#rho_{xy} in TOB", 50, -.5, .5));
527  myCorrHists.push_back(new TH1F("xyCorrTEC", "#rho_{xy} in TEC", 50, -1., 1.));
528  TDirectory *dirCorr = directory->mkdir("hitCorrelationHists", "correlations");
529  this->addToDirectory(myCorrHists, dirCorr);
530 
531  // PXB Survey
532  myPxbSurveyHists.push_back(new TH1F("PxbSurvChi2_lo", "#chi^{2} from PXB survey", 25, 0, 25));
533  myPxbSurveyHists.push_back(new TH1F("PxbSurvChi2_md", "#chi^{2} from PXB survey", 25, 0, 500));
534  myPxbSurveyHists.push_back(new TH1F("PxbSurvChi2_hi", "#chi^{2} from PXB survey", 25, 0, 50000));
535  myPxbSurveyHists.push_back(new TH1F("PxbSurvChi2prob", "Math::Prob(#chi^{2},4) from PXB survey", 25, 0, 1));
536  myPxbSurveyHists.push_back(new TH1F("PxbSurv_a0", "a_{0} from PXB survey", 100, -3000, 3000));
537  myPxbSurveyHists.push_back(new TH1F("PxbSurv_a0_abs", "fabs(a_{0}) from PXB survey", 100, 0, 3000));
538  myPxbSurveyHists.push_back(new TH1F("PxbSurv_a1", "a_{1} from PXB survey", 100, -6000, 6000));
539  myPxbSurveyHists.push_back(new TH1F("PxbSurv_a1_abs", "fabs(a_{1}) from PXB survey", 100, 0, 6000));
540  myPxbSurveyHists.push_back(
541  new TH1F("PxbSurv_scale", "scale (#sqrt{a_{2}^{2}+a_{3}^{2}}) from PXB survey", 100, 0, 1500));
542  myPxbSurveyHists.push_back(new TH1F("PxbSurv_phi", "angle(#atan{a_{3}/a_{4}}) from PXB survey", 100, -.05, .05));
543  TDirectory *dirPxbSurvey = directory->mkdir("PxbSurveyHists", "PxbSurvey");
544  this->addToDirectory(myPxbSurveyHists, dirPxbSurvey);
545 
546  oldDir->cd();
547  return true;
548 }
549 
550 //__________________________________________________________________
551 bool MillePedeMonitor::equidistLogBins(double *bins, int nBins, double first, double last) const {
552  // Filling 'bins' with borders of 'nBins' bins between 'first' and 'last'
553  // that are equidistant when viewed in log scale,
554  // so 'bins' must have length nBins+1;
555  // If 'first', 'last' or 'nBins' are not positive, failure is reported.
556 
557  if (nBins < 1 || first <= 0. || last <= 0.)
558  return false;
559 
560  bins[0] = first;
561  bins[nBins] = last;
562  const double firstLog = TMath::Log10(bins[0]);
563  const double lastLog = TMath::Log10(bins[nBins]);
564  for (int i = 1; i < nBins; ++i) {
565  bins[i] = TMath::Power(10., firstLog + i * (lastLog - firstLog) / (nBins));
566  }
567 
568  return true;
569 }
570 
571 //__________________________________________________________________
572 void MillePedeMonitor::fillTrack(const reco::Track *track) { this->fillTrack(track, myTrackHists1D, myTrackHists2D); }
573 
574 //__________________________________________________________________
575 void MillePedeMonitor::fillUsedTrack(const reco::Track *track, unsigned int nHitX, unsigned int nHitY) {
576  // these hist exist only for 'used track' hists:
577  static const int iUsedX = this->GetIndex(myUsedTrackHists1D, "usedHitsX");
578  myUsedTrackHists1D[iUsedX]->Fill(nHitX);
579  static const int iUsedY = this->GetIndex(myUsedTrackHists1D, "usedHitsY");
580  myUsedTrackHists1D[iUsedY]->Fill(nHitY);
581 
582  if (!track)
583  return;
584  this->fillTrack(track, myUsedTrackHists1D, myUsedTrackHists2D);
585 }
586 
587 //__________________________________________________________________
588 void MillePedeMonitor::fillTrack(const reco::Track *track,
589  std::vector<TH1 *> &trackHists1D,
590  std::vector<TH2 *> &trackHists2D) {
591  if (!track)
592  return;
593 
594  const reco::TrackBase::Vector p(track->momentum());
595 
596  static const int iPtLog = this->GetIndex(trackHists1D, "ptTrackLogBins");
597  trackHists1D[iPtLog]->Fill(track->pt());
598  static const int iPt = this->GetIndex(trackHists1D, "ptTrack");
599  trackHists1D[iPt]->Fill(track->pt());
600  static const int iP = this->GetIndex(trackHists1D, "pTrack");
601  trackHists1D[iP]->Fill(p.R());
602  static const int iEta = this->GetIndex(trackHists1D, "etaTrack");
603  trackHists1D[iEta]->Fill(p.Eta());
604  static const int iTheta = this->GetIndex(trackHists1D, "thetaTrack");
605  trackHists1D[iTheta]->Fill(p.Theta());
606  static const int iPhi = this->GetIndex(trackHists1D, "phiTrack");
607  trackHists1D[iPhi]->Fill(p.Phi());
608 
609  static const int iNhit = this->GetIndex(trackHists1D, "nHitTrack");
610  trackHists1D[iNhit]->Fill(track->numberOfValidHits());
611  static const int iNhitInvalid = this->GetIndex(trackHists1D, "nHitInvalidTrack");
612  trackHists1D[iNhitInvalid]->Fill(track->numberOfLostHits());
613 
614  int nhitinTIB = 0, nhitinTOB = 0, nhitinTID = 0;
615  int nhitinTEC = 0, nhitinBPIX = 0, nhitinFPIX = 0, nhitinPIXEL = 0;
616  int nhitinENDCAP = 0, nhitinENDCAPplus = 0, nhitinENDCAPminus = 0;
617  int nhitinTIDplus = 0, nhitinTIDminus = 0;
618  int nhitinFPIXplus = 0, nhitinFPIXminus = 0;
619  int nhitinTECplus = 0, nhitinTECminus = 0;
620  unsigned int thishit = 0;
621 
622  for (auto const &hit : track->recHits()) {
623  thishit++;
624  const DetId detId(hit->geographicalId());
625  const int subdetId = detId.subdetId();
626 
627  if (!hit->isValid())
628  continue; // only real hits count as in track->numberOfValidHits()
629  if (detId.det() != DetId::Tracker) {
630  edm::LogError("DetectorMismatch") << "@SUB=MillePedeMonitor::fillTrack"
631  << "DetId.det() != DetId::Tracker (=" << DetId::Tracker << "), but "
632  << detId.det() << ".";
633  }
634 
635  if (SiStripDetId::TIB == subdetId)
636  ++nhitinTIB;
637  else if (SiStripDetId::TOB == subdetId)
638  ++nhitinTOB;
639  else if (SiStripDetId::TID == subdetId) {
640  ++nhitinTID;
641  ++nhitinENDCAP;
642 
643  if (trackerTopology->tidIsZMinusSide(detId)) {
644  ++nhitinTIDminus;
645  ++nhitinENDCAPminus;
646  } else if (trackerTopology->tidIsZPlusSide(detId)) {
647  ++nhitinTIDplus;
648  ++nhitinENDCAPplus;
649  }
650  } else if (SiStripDetId::TEC == subdetId) {
651  ++nhitinTEC;
652  ++nhitinENDCAP;
653 
654  if (trackerTopology->tecIsZMinusSide(detId)) {
655  ++nhitinTECminus;
656  ++nhitinENDCAPminus;
657  } else if (trackerTopology->tecIsZPlusSide(detId)) {
658  ++nhitinTECplus;
659  ++nhitinENDCAPplus;
660  }
661  } else if (kBPIX == subdetId) {
662  ++nhitinBPIX;
663  ++nhitinPIXEL;
664  } else if (kFPIX == subdetId) {
665  ++nhitinFPIX;
666  ++nhitinPIXEL;
667 
668  if (trackerTopology->pxfSide(detId) == 1)
669  ++nhitinFPIXminus;
670  else if (trackerTopology->pxfSide(detId) == 2)
671  ++nhitinFPIXplus;
672  }
673 
674  } // end loop on hits
675 
676  static const int iNhit01 = this->GetIndex(trackHists1D, "nHitBPIXTrack");
677  trackHists1D[iNhit01]->Fill(nhitinBPIX);
678  static const int iNhit02 = this->GetIndex(trackHists1D, "nHitFPIXplusTrack");
679  trackHists1D[iNhit02]->Fill(nhitinFPIXplus);
680  static const int iNhit03 = this->GetIndex(trackHists1D, "nHitFPIXminusTrack");
681  trackHists1D[iNhit03]->Fill(nhitinFPIXminus);
682  static const int iNhit04 = this->GetIndex(trackHists1D, "nHitFPIXTrack");
683  trackHists1D[iNhit04]->Fill(nhitinFPIX);
684  static const int iNhit05 = this->GetIndex(trackHists1D, "nHitPIXELTrack");
685  trackHists1D[iNhit05]->Fill(nhitinPIXEL);
686  static const int iNhit06 = this->GetIndex(trackHists1D, "nHitTIBTrack");
687  trackHists1D[iNhit06]->Fill(nhitinTIB);
688  static const int iNhit07 = this->GetIndex(trackHists1D, "nHitTOBTrack");
689  trackHists1D[iNhit07]->Fill(nhitinTOB);
690  static const int iNhit08 = this->GetIndex(trackHists1D, "nHitTIDplusTrack");
691  trackHists1D[iNhit08]->Fill(nhitinTIDplus);
692  static const int iNhit09 = this->GetIndex(trackHists1D, "nHitTIDminusTrack");
693  trackHists1D[iNhit09]->Fill(nhitinTIDminus);
694  static const int iNhit10 = this->GetIndex(trackHists1D, "nHitTIDTrack");
695  trackHists1D[iNhit10]->Fill(nhitinTID);
696  static const int iNhit11 = this->GetIndex(trackHists1D, "nHitTECplusTrack");
697  trackHists1D[iNhit11]->Fill(nhitinTECplus);
698  static const int iNhit12 = this->GetIndex(trackHists1D, "nHitTECminusTrack");
699  trackHists1D[iNhit12]->Fill(nhitinTECminus);
700  static const int iNhit13 = this->GetIndex(trackHists1D, "nHitTECTrack");
701  trackHists1D[iNhit13]->Fill(nhitinTEC);
702  static const int iNhit14 = this->GetIndex(trackHists1D, "nHitENDCAPplusTrack");
703  trackHists1D[iNhit14]->Fill(nhitinENDCAPplus);
704  static const int iNhit15 = this->GetIndex(trackHists1D, "nHitENDCAPminusTrack");
705  trackHists1D[iNhit15]->Fill(nhitinENDCAPminus);
706  static const int iNhit16 = this->GetIndex(trackHists1D, "nHitENDCAPTrack");
707  trackHists1D[iNhit16]->Fill(nhitinENDCAP);
708  static const int iNhit17 = this->GetIndex(trackHists2D, "nHitENDCAPTrackMinusVsPlus");
709  trackHists2D[iNhit17]->Fill(nhitinENDCAPplus, nhitinENDCAPminus);
710 
711  if (track->innerOk()) {
712  const reco::TrackBase::Point &firstPoint(track->innerPosition());
713  static const int iR1 = this->GetIndex(trackHists1D, "r1Track");
714  trackHists1D[iR1]->Fill(firstPoint.Rho());
715  const double rSigned1 = (firstPoint.y() > 0 ? firstPoint.Rho() : -firstPoint.Rho());
716  static const int iR1Signed = this->GetIndex(trackHists1D, "r1TrackSigned");
717  trackHists1D[iR1Signed]->Fill(rSigned1);
718  static const int iZ1 = this->GetIndex(trackHists1D, "z1Track");
719  trackHists1D[iZ1]->Fill(firstPoint.Z());
720  static const int iZ1Full = this->GetIndex(trackHists1D, "z1TrackFull");
721  trackHists1D[iZ1Full]->Fill(firstPoint.Z());
722  static const int iY1 = this->GetIndex(trackHists1D, "y1Track");
723  trackHists1D[iY1]->Fill(firstPoint.Y());
724  static const int iPhi1 = this->GetIndex(trackHists1D, "phi1Track");
725  trackHists1D[iPhi1]->Fill(firstPoint.phi());
726  static const int iRz1Full = this->GetIndex(trackHists2D, "rz1TrackFull");
727  trackHists2D[iRz1Full]->Fill(firstPoint.Z(), rSigned1);
728  static const int iXy1 = this->GetIndex(trackHists2D, "xy1Track");
729  trackHists2D[iXy1]->Fill(firstPoint.X(), firstPoint.Y());
730  }
731 
732  if (track->outerOk()) {
733  const reco::TrackBase::Point &lastPoint(track->outerPosition());
734  static const int iRlast = this->GetIndex(trackHists1D, "rLastTrack");
735  trackHists1D[iRlast]->Fill(lastPoint.Rho());
736  static const int iZlast = this->GetIndex(trackHists1D, "zLastTrack");
737  trackHists1D[iZlast]->Fill(lastPoint.Z());
738  static const int iYlast = this->GetIndex(trackHists1D, "yLastTrack");
739  trackHists1D[iYlast]->Fill(lastPoint.Y());
740  static const int iPhiLast = this->GetIndex(trackHists1D, "phiLastTrack");
741  trackHists1D[iPhiLast]->Fill(lastPoint.phi());
742  }
743 
744  static const int iChi2Ndf = this->GetIndex(trackHists1D, "chi2PerNdf");
745  trackHists1D[iChi2Ndf]->Fill(track->normalizedChi2());
746 
747  static const int iImpZ = this->GetIndex(trackHists1D, "impParZ");
748  trackHists1D[iImpZ]->Fill(track->dz());
749  static const int iImpZerr = this->GetIndex(trackHists1D, "impParErrZ");
750  trackHists1D[iImpZerr]->Fill(track->dzError());
751 
752  static const int iImpRphi = this->GetIndex(trackHists1D, "impParRphi");
753  trackHists1D[iImpRphi]->Fill(track->d0());
754  static const int iImpRphiErr = this->GetIndex(trackHists1D, "impParErrRphi");
755  trackHists1D[iImpRphiErr]->Fill(track->d0Error());
756 }
757 
758 //__________________________________________________________________
759 void MillePedeMonitor::fillRefTrajectory(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr) {
760  static const int iValid = this->GetIndex(myTrajectoryHists1D, "validRefTraj");
761  if (refTrajPtr->isValid()) {
762  myTrajectoryHists1D[iValid]->Fill(1.);
763  } else {
764  myTrajectoryHists1D[iValid]->Fill(0.);
765  return;
766  }
767 
768  const AlgebraicSymMatrix &covMeasLoc = refTrajPtr->measurementErrors();
769  const AlgebraicVector &measurementsLoc = refTrajPtr->measurements();
770  const AlgebraicVector &trajectoryLoc = refTrajPtr->trajectoryPositions();
771  const TransientTrackingRecHit::ConstRecHitContainer &recHits = refTrajPtr->recHits();
772  const AlgebraicMatrix &derivatives = refTrajPtr->derivatives();
773 
774  // CHK
775  const int nRow = refTrajPtr->numberOfHitMeas();
776 
777  for (int iRow = 0; iRow < nRow; ++iRow) {
778  const double residuum = measurementsLoc[iRow] - trajectoryLoc[iRow];
779  const double covMeasLocRow = covMeasLoc[iRow][iRow];
780  const bool is2DhitRow = (!recHits[iRow / 2]->detUnit() // FIXME: as in MillePedeAlignmentAlgorithm::is2D()
781  || recHits[iRow / 2]->detUnit()->type().isTrackerPixel());
782  //the GeomDetUnit* is zero for composite hits (matched hits in the tracker,...).
783  if (TMath::Even(iRow)) { // local x
784  static const int iMeasLocX = this->GetIndex(myTrajectoryHists1D, "measLocX");
785  myTrajectoryHists1D[iMeasLocX]->Fill(measurementsLoc[iRow]);
786  static const int iTrajLocX = this->GetIndex(myTrajectoryHists1D, "trajLocX");
787  myTrajectoryHists1D[iTrajLocX]->Fill(trajectoryLoc[iRow]);
788  static const int iResidLocX = this->GetIndex(myTrajectoryHists1D, "residLocX");
789  myTrajectoryHists1D[iResidLocX]->Fill(residuum);
790 
791  static const int iMeasLocXvsHit = this->GetIndex(myTrajectoryHists2D, "measLocXvsHit");
792  myTrajectoryHists2D[iMeasLocXvsHit]->Fill(iRow / 2, measurementsLoc[iRow]);
793  static const int iTrajLocXvsHit = this->GetIndex(myTrajectoryHists2D, "trajLocXvsHit");
794  myTrajectoryHists2D[iTrajLocXvsHit]->Fill(iRow / 2, trajectoryLoc[iRow]);
795  static const int iResidLocXvsHit = this->GetIndex(myTrajectoryHists2D, "residLocXvsHit");
796  myTrajectoryHists2D[iResidLocXvsHit]->Fill(iRow / 2, residuum);
797 
798  if (covMeasLocRow > 0.) {
799  static const int iReduResidLocX = this->GetIndex(myTrajectoryHists1D, "reduResidLocX");
800  myTrajectoryHists1D[iReduResidLocX]->Fill(residuum / TMath::Sqrt(covMeasLocRow));
801  static const int iReduResidLocXvsHit = this->GetIndex(myTrajectoryHists2D, "reduResidLocXvsHit");
802  myTrajectoryHists2D[iReduResidLocXvsHit]->Fill(iRow / 2, residuum / TMath::Sqrt(covMeasLocRow));
803  }
804  } else if (is2DhitRow) { // local y, 2D detectors only
805 
806  static const int iMeasLocY = this->GetIndex(myTrajectoryHists1D, "measLocY");
807  myTrajectoryHists1D[iMeasLocY]->Fill(measurementsLoc[iRow]);
808  static const int iTrajLocY = this->GetIndex(myTrajectoryHists1D, "trajLocY");
809  myTrajectoryHists1D[iTrajLocY]->Fill(trajectoryLoc[iRow]);
810  static const int iResidLocY = this->GetIndex(myTrajectoryHists1D, "residLocY");
811  myTrajectoryHists1D[iResidLocY]->Fill(residuum);
812 
813  static const int iMeasLocYvsHit = this->GetIndex(myTrajectoryHists2D, "measLocYvsHit");
814  myTrajectoryHists2D[iMeasLocYvsHit]->Fill(iRow / 2, measurementsLoc[iRow]);
815  static const int iTrajLocYvsHit = this->GetIndex(myTrajectoryHists2D, "trajLocYvsHit");
816  myTrajectoryHists2D[iTrajLocYvsHit]->Fill(iRow / 2, trajectoryLoc[iRow]);
817  static const int iResidLocYvsHit = this->GetIndex(myTrajectoryHists2D, "residLocYvsHit");
818  myTrajectoryHists2D[iResidLocYvsHit]->Fill(iRow / 2, residuum);
819 
820  if (covMeasLocRow > 0.) {
821  static const int iReduResidLocY = this->GetIndex(myTrajectoryHists1D, "reduResidLocY");
822  myTrajectoryHists1D[iReduResidLocY]->Fill(residuum / TMath::Sqrt(covMeasLocRow));
823  static const int iReduResidLocYvsHit = this->GetIndex(myTrajectoryHists2D, "reduResidLocYvsHit");
824  myTrajectoryHists2D[iReduResidLocYvsHit]->Fill(iRow / 2, residuum / TMath::Sqrt(covMeasLocRow));
825  }
826  }
827 
828  float nHitRow = iRow / 2; // '/2' not '/2.'!
829  if (TMath::Odd(iRow))
830  nHitRow += 0.5; // y-hit gets 0.5
831  // correlations
832  for (int iCol = iRow + 1; iCol < nRow; ++iCol) {
833  double rho = TMath::Sqrt(covMeasLocRow + covMeasLoc[iCol][iCol]);
834  rho = (0. == rho ? -2 : covMeasLoc[iRow][iCol] / rho);
835  float nHitCol = iCol / 2; //cf. comment nHitRow
836  if (TMath::Odd(iCol))
837  nHitCol += 0.5; // dito
838  // myProfileCorr->Fill(nHitRow, nHitCol, TMath::Abs(rho));
839  if (0. == rho)
840  continue;
841  static const int iProfileCorr = this->GetIndex(myTrajectoryHists2D, "profCorr");
842  myTrajectoryHists2D[iProfileCorr]->Fill(nHitRow, nHitCol, TMath::Abs(rho));
843  if (iRow + 1 == iCol && TMath::Even(iRow)) { // i.e. if iRow is x and iCol the same hit's y
844  // static const int iProfileCorrXy = this->GetIndex(myTrajectoryHists2D,"profCorrOffXy");
845  // myTrajectoryHists2D[iProfileCorrOffXy]->Fill(iRow/2, rho);
846  static const int iHitCorrXy = this->GetIndex(myTrajectoryHists2D, "hitCorrXy");
847  myTrajectoryHists2D[iHitCorrXy]->Fill(iRow / 2, rho);
848  static const int iHitCorrXyLog = this->GetIndex(myTrajectoryHists2D, "hitCorrXyLog");
849  myTrajectoryHists2D[iHitCorrXyLog]->Fill(iRow / 2, TMath::Abs(rho));
850  if (is2DhitRow) {
851  static const int iHitCorrXyValid = this->GetIndex(myTrajectoryHists2D, "hitCorrXyValid");
852  myTrajectoryHists2D[iHitCorrXyValid]->Fill(iRow / 2, rho); // nhitRow??
853  static const int iHitCorrXyLogValid = this->GetIndex(myTrajectoryHists2D, "hitCorrXyLogValid");
854  myTrajectoryHists2D[iHitCorrXyLogValid]->Fill(iRow / 2, TMath::Abs(rho)); // nhitRow??
855  }
856  } else {
857  static const int iProfCorrOffXy = this->GetIndex(myTrajectoryHists2D, "profCorrOffXy");
858  myTrajectoryHists2D[iProfCorrOffXy]->Fill(nHitRow, nHitCol, TMath::Abs(rho));
859  static const int iHitCorrOffBlock = this->GetIndex(myTrajectoryHists2D, "hitCorrOffBlock");
860  myTrajectoryHists2D[iHitCorrOffBlock]->Fill(nHitRow, rho);
861  static const int iHitCorOffBlkLg = this->GetIndex(myTrajectoryHists2D, "hitCorrOffBlockLog");
862  myTrajectoryHists2D[iHitCorOffBlkLg]->Fill(nHitRow, TMath::Abs(rho));
863  }
864  } // end loop on columns of covariance
865 
866  // derivatives
867  for (int iCol = 0; iCol < derivatives.num_col(); ++iCol) {
868  static const int iProfDerivatives = this->GetIndex(myTrajectoryHists2D, "profDerivatives");
869  myTrajectoryHists2D[iProfDerivatives]->Fill(nHitRow, iCol, derivatives[iRow][iCol]);
870  static const int iDerivatives = this->GetIndex(myTrajectoryHists2D, "derivatives");
871  myTrajectoryHists2D[iDerivatives]->Fill(iCol, derivatives[iRow][iCol]);
872  static const int iDerivativesLog = this->GetIndex(myTrajectoryHists2D, "derivativesLog");
873  myTrajectoryHists2D[iDerivativesLog]->Fill(iCol, TMath::Abs(derivatives[iRow][iCol]));
874  static const int iDerivativesPhi = this->GetIndex(myTrajectoryHists2D, "derivativesVsPhi");
875  myTrajectoryHists2D[iDerivativesPhi]->Fill(recHits[iRow / 2]->det()->position().phi(), derivatives[iRow][iCol]);
876  }
877  } // end loop on rows of covarianvce
878 }
879 
880 //____________________________________________________________________
881 void MillePedeMonitor::fillDerivatives(const ConstRecHitPointer &recHit,
882  const float *localDerivs,
883  unsigned int nLocal,
884  const float *globalDerivs,
885  unsigned int nGlobal,
886  const int *labels) {
887  const double phi = recHit->det()->position().phi();
888 
889  static const int iLocPar = this->GetIndex(myDerivHists2D, "localDerivsPar");
890  static const int iLocPhi = this->GetIndex(myDerivHists2D, "localDerivsPhi");
891  static const int iLocParLog = this->GetIndex(myDerivHists2D, "localDerivsParLog");
892  static const int iLocPhiLog = this->GetIndex(myDerivHists2D, "localDerivsPhiLog");
893  for (unsigned int i = 0; i < nLocal; ++i) {
894  myDerivHists2D[iLocPar]->Fill(i, localDerivs[i]);
895  myDerivHists2D[iLocPhi]->Fill(phi, localDerivs[i]);
896  if (localDerivs[i]) {
897  myDerivHists2D[iLocParLog]->Fill(i, TMath::Abs(localDerivs[i]));
898  myDerivHists2D[iLocPhiLog]->Fill(phi, TMath::Abs(localDerivs[i]));
899  }
900  }
901 
902  static const int iGlobPar = this->GetIndex(myDerivHists2D, "globalDerivsPar");
903  static const int iGlobPhi = this->GetIndex(myDerivHists2D, "globalDerivsPhi");
904  static const int iGlobParLog = this->GetIndex(myDerivHists2D, "globalDerivsParLog");
905  static const int iGlobPhiLog = this->GetIndex(myDerivHists2D, "globalDerivsPhiLog");
906  // static const int iGlobPhiLog2 = this->GetIndex(myDerivHists2D, "globalDerivsPhiLog2");
907  for (unsigned int i = 0; i < nGlobal; ++i) {
908  const unsigned int parNum = (labels ? (labels[i] % PedeLabelerBase::theMaxNumParam) - 1 : i);
909  myDerivHists2D[iGlobPar]->Fill(parNum, globalDerivs[i]);
910  myDerivHists2D[iGlobPhi]->Fill(phi, globalDerivs[i]);
911  if (globalDerivs[i]) {
912  myDerivHists2D[iGlobParLog]->Fill(parNum, TMath::Abs(globalDerivs[i]));
913  myDerivHists2D[iGlobPhiLog]->Fill(phi, TMath::Abs(globalDerivs[i]));
914  // myDerivHists2D[iGlobPhiLog2]->Fill(phi, TMath::Abs(globalDerivs[i]));
915  }
916  }
917 }
918 
919 //____________________________________________________________________
920 void MillePedeMonitor::fillResiduals(const ConstRecHitPointer &recHit,
921  const TrajectoryStateOnSurface &tsos,
922  unsigned int nHit,
923  float residuum,
924  float sigma,
925  bool isY) {
926  // isY == false: x-measurements
927  // isY == true: y-measurements
928  const GlobalPoint detPos(recHit->det()->position());
929  const double phi = detPos.phi();
930  // const double rSigned = (detPos.y() > 0. ? detPos.perp() : -detPos.perp());
931 
932  static const int iResPhi = this->GetIndex(myResidHists2D, "residPhi");
933  myResidHists2D[iResPhi]->Fill(phi, residuum);
934  static const int iSigPhi = this->GetIndex(myResidHists2D, "sigmaPhi");
935  myResidHists2D[iSigPhi]->Fill(phi, sigma);
936  if (sigma) {
937  static const int iReduResPhi = this->GetIndex(myResidHists2D, "reduResidPhi");
938  myResidHists2D[iReduResPhi]->Fill(phi, residuum / sigma);
939  }
940 
941  // if (isY) {
942  // static const int iResYXy = this->GetIndex(myResidHists2D, "residYProfXy");
943  // myResidHists2D[iResYXy]->Fill(detPos.x(), detPos.y(), TMath::Abs(residuum));
944  // static const int iResYZr = this->GetIndex(myResidHists2D, "residYProfZr");
945  // myResidHists2D[iResYZr]->Fill(detPos.z(), rSigned, TMath::Abs(residuum));
946  // static const int iSigmaYXy = this->GetIndex(myResidHists2D, "sigmaYProfXy");
947  // myResidHists2D[iSigmaYXy]->Fill(detPos.x(), detPos.y(), sigma);
948  // static const int iSigmaYZr = this->GetIndex(myResidHists2D, "sigmaYProfZr");
949  // myResidHists2D[iSigmaYZr]->Fill(detPos.z(), rSigned, sigma);
950  // if (sigma) {
951  // static const int iReduResYXy = this->GetIndex(myResidHists2D, "reduResidYProfXy");
952  // myResidHists2D[iReduResYXy]->Fill(detPos.x(), detPos.y(), TMath::Abs(residuum/sigma));
953  // static const int iReduResYZr = this->GetIndex(myResidHists2D, "reduResidYProfZr");
954  // myResidHists2D[iReduResYZr]->Fill(detPos.z(), rSigned, TMath::Abs(residuum/sigma));
955  // }
956  // } else {
957  // static const int iResXXy = this->GetIndex(myResidHists2D, "residXProfXy");
958  // myResidHists2D[iResXXy]->Fill(detPos.x(), detPos.y(), TMath::Abs(residuum));
959  // static const int iResXZr = this->GetIndex(myResidHists2D, "residXProfZr");
960  // myResidHists2D[iResXZr]->Fill(detPos.z(), rSigned, TMath::Abs(residuum));
961  // static const int iSigmaXXy = this->GetIndex(myResidHists2D, "sigmaXProfXy");
962  // myResidHists2D[iSigmaXXy]->Fill(detPos.x(), detPos.y(), sigma);
963  // static const int iSigmaXZr = this->GetIndex(myResidHists2D, "sigmaXProfZr");
964  // myResidHists2D[iSigmaXZr]->Fill(detPos.z(), rSigned, sigma);
965  // if (sigma) {
966  // static const int iReduResXXy = this->GetIndex(myResidHists2D, "reduResidXProfXy");
967  // myResidHists2D[iReduResXXy]->Fill(detPos.x(), detPos.y(), TMath::Abs(residuum/sigma));
968  // static const int iReduResXZr = this->GetIndex(myResidHists2D, "reduResidXProfZr");
969  // myResidHists2D[iReduResXZr]->Fill(detPos.z(), rSigned, TMath::Abs(residuum/sigma));
970  // }
971  // }
972 
973  const LocalVector mom(tsos.localDirection()); // mom.z()==0. impossible for TSOS:
974  const float phiSensToNorm = TMath::ATan(TMath::Abs((isY ? mom.y() : mom.x()) / mom.z()));
975 
976  std::vector<std::vector<TH1 *> > &histArrayVec = (isY ? myResidHistsVec1DY : myResidHistsVec1DX);
977  std::vector<TH1 *> &hitHists = (isY ? myResidHitHists1DY : myResidHitHists1DX);
978 
979  // call with histArrayVec[0] first: 'static' inside is referring to "subdet-less" names (X/Y irrelevant)
980  this->fillResidualHists(histArrayVec[0], phiSensToNorm, residuum, sigma);
981  this->fillResidualHitHists(hitHists, phiSensToNorm, residuum, sigma, nHit);
982 
983  const DetId detId(recHit->det()->geographicalId());
984  if (detId.det() == DetId::Tracker) {
985  // const GeomDet::SubDetector subDetNum = recHit->det()->subDetector();
986  unsigned int subDetNum = detId.subdetId();
987  if (subDetNum < histArrayVec.size() && subDetNum > 0) {
988  this->fillResidualHists(histArrayVec[subDetNum], phiSensToNorm, residuum, sigma);
989  } else {
990  if (detId != AlignableBeamSpot::detId())
991  edm::LogWarning("Alignment") << "@SUB=MillePedeMonitor::fillResiduals"
992  << "Expect subDetNum from 1 to 6, got " << subDetNum;
993  }
994  }
995 }
996 
997 //____________________________________________________________________
998 void MillePedeMonitor::fillResidualHists(const std::vector<TH1 *> &hists,
999  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 }
MillePedeMonitor::myResidHistsVec1DY
std::vector< std::vector< TH1 * > > myResidHistsVec1DY
[0]=all [1]=TPB [2]=TPE [3]=TIB [4]=TID [5]=TOB [6]=TEC
Definition: MillePedeMonitor.h:109
Pi
const double Pi
Definition: CosmicMuonParameters.h:18
MillePedeMonitor::myResidHists2D
std::vector< TH2 * > myResidHists2D
Definition: MillePedeMonitor.h:107
MillePedeMonitor::GetIndex
int GetIndex(const std::vector< OBJECT_TYPE * > &vec, const TString &name)
Definition: MillePedeMonitor.h:120
reco::TrackBase::d0Error
double d0Error() const
error on d0
Definition: TrackBase.h:772
MillePedeMonitor::fillCorrelations2D
void fillCorrelations2D(float corr, const TransientTrackingRecHit::ConstRecHitPointer &hit)
Definition: MillePedeMonitor.cc:1096
mps_fire.i
i
Definition: mps_fire.py:428
MillePedeMonitor::myUsedTrackHists1D
std::vector< TH1 * > myUsedTrackHists1D
Definition: MillePedeMonitor.h:102
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:69
reco::TrackBase::d0
double d0() const
dxy parameter in perigee convention (d0 = -dxy)
Definition: TrackBase.h:611
MillePedeMonitor::myCorrHists
std::vector< TH1 * > myCorrHists
Definition: MillePedeMonitor.h:113
MillePedeMonitor::fillFrameToFrame
void fillFrameToFrame(const AlignableDetOrUnitPtr &aliDet, const Alignable *ali)
Definition: MillePedeMonitor.cc:1066
MillePedeMonitor::myPxbSurveyHists
std::vector< TH1 * > myPxbSurveyHists
Definition: MillePedeMonitor.h:114
reco::TrackBase::normalizedChi2
double normalizedChi2() const
chi-squared divided by n.d.o.f. (or chi-squared * 1e6 if n.d.o.f. is zero)
Definition: TrackBase.h:593
SiStripDetId::TID
static constexpr auto TID
Definition: SiStripDetId.h:38
TrajectoryStateOnSurface::localDirection
LocalVector localDirection() const
Definition: TrajectoryStateOnSurface.h:76
MillePedeMonitor::fillPxbSurveyHistsLocalPars
void fillPxbSurveyHistsLocalPars(const float &a0, const float &a1, const float &S, const float &phi)
Definition: MillePedeMonitor.cc:1123
PV3DBase::phi
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
reco::Track::innerOk
bool innerOk() const
return true if the innermost hit is valid
Definition: Track.h:53
reco::TrackBase::numberOfLostHits
unsigned short numberOfLostHits() const
number of cases where track crossed a layer without getting a hit.
Definition: TrackBase.h:801
HLT_FULL_cff.track
tuple track
Definition: HLT_FULL_cff.py:11953
MillePedeMonitor::trackerTopology
const TrackerTopology * trackerTopology
Definition: MillePedeMonitor.h:116
reco::TrackBase::momentum
const Vector & momentum() const
track momentum vector
Definition: TrackBase.h:664
edm::LogError
Log< level::Error, false > LogError
Definition: MessageLogger.h:123
MillePedeMonitor::cloneHists
std::vector< OBJECT_TYPE * > cloneHists(const std::vector< OBJECT_TYPE * > &orgs, const TString &namAd, const TString &titAd) const
Definition: MillePedeMonitor.h:133
reco::Track::outerPosition
const math::XYZPoint & outerPosition() const
position of the outermost hit
Definition: Track.h:62
MillePedeMonitor::myResidHitHists1DX
std::vector< TH1 * > myResidHitHists1DX
[0]=all [1]=TPB [2]=TPE [3]=TIB [4]=TID [5]=TOB [6]=TEC
Definition: MillePedeMonitor.h:110
MillePedeMonitor::MillePedeMonitor
MillePedeMonitor(const TrackerTopology *tTopo, const char *rootFile="trackMonitor.root")
Definition: MillePedeMonitor.cc:45
MillePedeMonitor::fillResiduals
void fillResiduals(const TransientTrackingRecHit::ConstRecHitPointer &recHit, const TrajectoryStateOnSurface &tsos, unsigned int nHit, float residuum, float sigma, bool isY)
Definition: MillePedeMonitor.cc:920
TrackerTopology::tecIsZMinusSide
bool tecIsZMinusSide(const DetId &id) const
Definition: TrackerTopology.h:262
MillePedeMonitor::fillPxbSurveyHistsChi2
void fillPxbSurveyHistsChi2(const float &chi2)
Definition: MillePedeMonitor.cc:1111
MillePedeMonitor::myUsedTrackHists2D
std::vector< TH2 * > myUsedTrackHists2D
Definition: MillePedeMonitor.h:103
PedeLabelerBase::theMaxNumParam
static const unsigned int theMaxNumParam
Definition: PedeLabelerBase.h:109
TrackerTopology::tidIsZMinusSide
bool tidIsZMinusSide(const DetId &id) const
Definition: TrackerTopology.h:259
reco::Track::innerPosition
const math::XYZPoint & innerPosition() const
position of the innermost hit
Definition: Track.h:56
kFPIX
const int kFPIX
Definition: AlignmentTrackSelector.cc:26
ConstRecHitPointer
TransientTrackingRecHit::ConstRecHitPointer ConstRecHitPointer
Definition: MillePedeAlignmentAlgorithm.cc:70
MillePedeMonitor::addToDirectory
void addToDirectory(const std::vector< OBJECT_TYPE * > &objs, TDirectory *dir) const
Definition: MillePedeMonitor.h:154
dtT0Analyzer_cfg.rootFileName
tuple rootFileName
Definition: dtT0Analyzer_cfg.py:26
MillePedeMonitor::myTrackHists2D
std::vector< TH2 * > myTrackHists2D
Definition: MillePedeMonitor.h:101
MillePedeMonitor::myDerivHists2D
std::vector< TH2 * > myDerivHists2D
Definition: MillePedeMonitor.h:106
HLT_FULL_cff.chi2
tuple chi2
Definition: HLT_FULL_cff.py:10024
AlgebraicMatrix
CLHEP::HepMatrix AlgebraicMatrix
Definition: AlgebraicObjects.h:14
MillePedeMonitor::fillUsedTrack
void fillUsedTrack(const reco::Track *track, unsigned int nHitX, unsigned int nHitY)
Definition: MillePedeMonitor.cc:575
TrackingRecHit::ConstRecHitPointer
std::shared_ptr< TrackingRecHit const > ConstRecHitPointer
Definition: TrackingRecHit.h:25
MillePedeMonitor::myTrajectoryHists1D
std::vector< TH1 * > myTrajectoryHists1D
Definition: MillePedeMonitor.h:104
reco::TrackBase::pt
double pt() const
track transverse momentum
Definition: TrackBase.h:637
MillePedeMonitor::myRootDir
TDirectory * myRootDir
Definition: MillePedeMonitor.h:97
Abs
T Abs(T a)
Definition: MathUtil.h:49
DetId::subdetId
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:48
reco::Track::recHits
auto recHits() const
Access to reconstructed hits on the track.
Definition: Track.h:85
reco::TrackBase::numberOfValidHits
unsigned short numberOfValidHits() const
number of valid hits found
Definition: TrackBase.h:798
MillePedeMonitor::myTrajectoryHists2D
std::vector< TH2 * > myTrajectoryHists2D
Definition: MillePedeMonitor.h:105
MillePedeMonitor::init
bool init(TDirectory *directory)
Definition: MillePedeMonitor.cc:72
reco::TrackBase::Point
math::XYZPoint Point
point in the space
Definition: TrackBase.h:80
FrameToFrameDerivative::frameToFrameDerivative
AlgebraicMatrix frameToFrameDerivative(const Alignable *object, const Alignable *composedObject) const
Definition: FrameToFrameDerivative.cc:16
SiStripDetId::TOB
static constexpr auto TOB
Definition: SiStripDetId.h:39
TrackingRecHit::ConstRecHitContainer
std::vector< ConstRecHitPointer > ConstRecHitContainer
Definition: TrackingRecHit.h:32
reco::TrackBase::dz
double dz() const
dz parameter (= dsz/cos(lambda)). This is the track z0 w.r.t (0,0,0) only if the refPoint is close to...
Definition: TrackBase.h:622
reco::TrackBase::dzError
double dzError() const
error on dz
Definition: TrackBase.h:778
DetId
Definition: DetId.h:17
AlgebraicVector
CLHEP::HepVector AlgebraicVector
Definition: AlgebraicObjects.h:13
MillePedeMonitor::myFrame2FrameHists2D
std::vector< TH2 * > myFrame2FrameHists2D
Definition: MillePedeMonitor.h:112
MillePedeMonitor::equidistLogBins
bool equidistLogBins(double *bins, int nBins, double first, double last) const
Definition: MillePedeMonitor.cc:551
AlignableBeamSpot::detId
static const DetId detId()
Definition: AlignableBeamSpot.h:90
TrackerTopology::tidIsZPlusSide
bool tidIsZPlusSide(const DetId &id) const
Definition: TrackerTopology.h:258
reco::Track::outerOk
bool outerOk() const
return true if the outermost hit is valid
Definition: Track.h:50
a0
static constexpr float a0
Definition: L1EGammaCrystalsEmulatorProducer.cc:82
HLT_FULL_cff.recHits
tuple recHits
Definition: HLT_FULL_cff.py:85944
MillePedeMonitor::fillResidualHists
void fillResidualHists(const std::vector< TH1 * > &hists, float phiSensToNorm, float residuum, float sigma)
Definition: MillePedeMonitor.cc:998
MillePedeMonitor::fillTrack
void fillTrack(const reco::Track *track)
Definition: MillePedeMonitor.cc:572
nHits
caConstants::TupleMultiplicity const CAHitNtupletGeneratorKernelsGPU::HitToTuple const cms::cuda::AtomicPairCounter GPUCACell const *__restrict__ uint32_t const *__restrict__ gpuPixelDoublets::CellNeighborsVector const gpuPixelDoublets::CellTracksVector const GPUCACell::OuterHitOfCell const int32_t nHits
Definition: CAHitNtupletGeneratorKernelsImpl.h:43
SiStripDetId::TIB
static constexpr auto TIB
Definition: SiStripDetId.h:37
TrackerTopology::tecIsZPlusSide
bool tecIsZPlusSide(const DetId &id) const
Definition: TrackerTopology.h:261
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
compare.hists
list hists
Definition: compare.py:318
MillePedeMonitor::myTrackHists1D
std::vector< TH1 * > myTrackHists1D
Definition: MillePedeMonitor.h:100
eostools.mkdir
def mkdir
Definition: eostools.py:251
MillePedeMonitor::myResidHitHists1DY
std::vector< TH1 * > myResidHitHists1DY
Definition: MillePedeMonitor.h:111
position
static int position[264][3]
Definition: ReadPGInfo.cc:289
MillePedeMonitor::fillRefTrajectory
void fillRefTrajectory(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr)
Definition: MillePedeMonitor.cc:759
AlgebraicSymMatrix
CLHEP::HepSymMatrix AlgebraicSymMatrix
Definition: AlgebraicObjects.h:15
TrackerTopology::pxfSide
unsigned int pxfSide(const DetId &id) const
Definition: TrackerTopology.h:192
MillePedeMonitor::fillDerivatives
void fillDerivatives(const TransientTrackingRecHit::ConstRecHitPointer &recHit, const float *localDerivs, unsigned int nLocal, const float *globalDerivs, unsigned int nGlobal, const int *labels)
Definition: MillePedeMonitor.cc:881
Alignable::globalPosition
const PositionType & globalPosition() const
Return the global position of the object.
Definition: Alignable.h:135
dqmdumpme.last
tuple last
Definition: dqmdumpme.py:56
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:122
MillePedeMonitor::myResidHistsVec1DX
std::vector< std::vector< TH1 * > > myResidHistsVec1DX
Definition: MillePedeMonitor.h:108
h
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
reco::TrackBase::Vector
math::XYZVector Vector
spatial vector
Definition: TrackBase.h:77
HLT_FULL_cff.labels
tuple labels
Definition: HLT_FULL_cff.py:9999
SiStripDetId::TEC
static constexpr auto TEC
Definition: SiStripDetId.h:40
MillePedeMonitor::fillResidualHitHists
void fillResidualHitHists(const std::vector< TH1 * > &hists, float angle, float residuum, float sigma, unsigned int nHit)
Definition: MillePedeMonitor.cc:1043
kBPIX
const int kBPIX
Definition: AlignmentTrackSelector.cc:25
angle
T angle(T x1, T y1, T z1, T x2, T y2, T z2)
Definition: angle.h:11
Generated for CMSSW Reference Manual by doxygen 1.8.5