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