Reference Manual

Go to the documentation of this file.
00001 
00011 #include "DataFormats/GeometrySurface/interface/Surface.h" 
00012 #include "Alignment/MillePedeAlignmentAlgorithm/interface/MillePedeMonitor.h"
00013 
00014 #include "DataFormats/TrackReco/interface/Track.h"
00015 
00016 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00017 
00018 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00019 #include <Geometry/CommonDetUnit/interface/GeomDetUnit.h>
00020 #include <Geometry/CommonDetUnit/interface/GeomDetType.h>
00021 
00022 #include "Alignment/CommonAlignment/interface/Alignable.h"
00023 #include "Alignment/CommonAlignment/interface/AlignableBeamSpot.h"
00024 #include "Alignment/CommonAlignment/interface/AlignableDetOrUnitPtr.h"
00025 #include "Alignment/CommonAlignmentParametrization/interface/FrameToFrameDerivative.h"
00026 
00027 #include "DataFormats/DetId/interface/DetId.h"
00028 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
00029 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
00030 #include "DataFormats/SiStripDetId/interface/TIDDetId.h"
00031 #include "DataFormats/SiStripDetId/interface/TECDetId.h"
00032 #include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
00033 const int kBPIX = PixelSubdetector::PixelBarrel;
00034 const int kFPIX = PixelSubdetector::PixelEndcap;
00035 
00036 #include <TProfile2D.h>
00037 #include <TFile.h>
00038 #include <TDirectory.h>
00039 #include <TMath.h>
00040 
00041 #include "Alignment/MillePedeAlignmentAlgorithm/interface/PedeLabelerBase.h"
00042 
00043 typedef TransientTrackingRecHit::ConstRecHitPointer   ConstRecHitPointer;
00044 
00045 //__________________________________________________________________
00046 MillePedeMonitor::MillePedeMonitor(const char *rootFileName)
00047   : myRootDir(0), myDeleteDir(false)
00048 {
00049   myRootDir = TFile::Open(rootFileName, "recreate");
00050   myDeleteDir = true;
00051 
00052   this->init(myRootDir);
00053 }
00054 
00055 //__________________________________________________________________
00056 MillePedeMonitor::MillePedeMonitor(TDirectory *rootDir) 
00057   : myRootDir(0), myDeleteDir(false)
00058 {
00059   //  cout << "MillePedeMonitor using input TDirectory" << endl;
00060 
00061   myRootDir = rootDir;
00062   myDeleteDir = false;
00063 
00064   this->init(myRootDir);
00065 }
00066 
00067 //__________________________________________________________________
00068 MillePedeMonitor::~MillePedeMonitor()
00069 {
00070 
00071   myRootDir->Write();
00072   if (myDeleteDir) delete myRootDir; //hists are deleted with their directory
00073 }
00074 
00075 //__________________________________________________________________
00076 bool MillePedeMonitor::init(TDirectory *directory)
00077 {
00078   if (!directory) return false;
00079   TDirectory *oldDir = gDirectory;
00080 
00081   const int kNumBins = 20;
00082   double binsPt[kNumBins+1] = {0.}; // fully initialised with 0.
00083   if (!this->equidistLogBins(binsPt, kNumBins, 0.8, 100.)) {
00084 //     cerr << "MillePedeMonitor::init: problem with log bins" << endl;
00085   }
00086   const int nHits = 35;
00087 
00088   myTrackHists1D.push_back(new TH1F("ptTrackLogBins",  "p_{t}(track);p_{t} [GeV]",
00089                                     kNumBins, binsPt));
00090 
00091   myTrackHists1D.push_back(new TH1F("ptTrack",  "p_{t}(track);p_{t} [GeV]",
00092                                     kNumBins, binsPt[0], binsPt[kNumBins]));
00093   myTrackHists1D.push_back(new TH1F("pTrack",  "p(track);p [GeV]",
00094                                     kNumBins, binsPt[0], 1.3*binsPt[kNumBins]));
00095   myTrackHists1D.push_back(new TH1F("etaTrack", "#eta(track);#eta", 26, -2.6, 2.6));
00096   myTrackHists1D.push_back(new TH1F("thetaTrack", "#theta(track);#theta", 100, 0., TMath::Pi()));
00097   myTrackHists1D.push_back(new TH1F("phiTrack", "#phi(track);#phi", 15, -TMath::Pi(), TMath::Pi()));
00098 
00099   myTrackHists1D.push_back(new TH1F("nHitTrack", "N_{hit}(track);N_{hit}", 40, 5., 45.));
00100   myTrackHists1D.push_back(new TH1F("nHitBPIXTrack", "N_{hit, BPIX}(track);N_{hit}", 45, 0., 45.));
00101   myTrackHists1D.push_back(new TH1F("nHitFPIXTrack", "N_{hit, FPIX}(track);N_{hit}", 45, 0., 45.));
00102   myTrackHists1D.push_back(new TH1F("nHitFPIXplusTrack", "N_{hit, BPIXplus}(track);N_{hit}", 45, 0., 45.));
00103   myTrackHists1D.push_back(new TH1F("nHitFPIXminusTrack", "N_{hit, BPIXminus}(track);N_{hit}", 45, 0., 45.));
00104   myTrackHists1D.push_back(new TH1F("nHitPIXELTrack", "N_{hit, PIXEL}(track);N_{hit}", 45, 0., 45.));
00105   myTrackHists1D.push_back(new TH1F("nHitTIBTrack", "N_{hit, TIB}(track);N_{hit}", 45, 0., 45.));
00106   myTrackHists1D.push_back(new TH1F("nHitTOBTrack", "N_{hit, TOB}(track);N_{hit}", 45, 0., 45.));
00107   myTrackHists1D.push_back(new TH1F("nHitTIDplusTrack", "N_{hit, TIDplus}(track);N_{hit}", 45, 0., 45.));
00108   myTrackHists1D.push_back(new TH1F("nHitTIDminusTrack", "N_{hit, TIDminus}(track);N_{hit}", 45, 0., 45.));
00109   myTrackHists1D.push_back(new TH1F("nHitTIDTrack", "N_{hit, TID}(track);N_{hit}", 45, 0., 45.));
00110   myTrackHists1D.push_back(new TH1F("nHitTECplusTrack", "N_{hit, TECplus}(track);N_{hit}", 45, 0., 45.));
00111   myTrackHists1D.push_back(new TH1F("nHitTECminusTrack", "N_{hit, TECminus}(track);N_{hit}", 45, 0., 45.));
00112   myTrackHists1D.push_back(new TH1F("nHitTECTrack", "N_{hit, TEC}(track);N_{hit}", 45, 0., 45.));
00113   myTrackHists1D.push_back(new TH1F("nHitENDCAPplusTrack", "N_{hit, ENDCAPplus}(track);N_{hit}", 45, 0., 45.));
00114   myTrackHists1D.push_back(new TH1F("nHitENDCAPminusTrack", "N_{hit, ENDCAPminus}(track);N_{hit}", 45, 0., 45.));
00115   myTrackHists1D.push_back(new TH1F("nHitENDCAPTrack", "N_{hit, ENDCAP}(track);N_{hit}", 45, 0., 45.));
00116   myTrackHists2D.push_back(new TH2F("nHitENDCAPTrackMinusVsPlus", "N_{hit, ENDCAP}(track);N_{hit, plus};N_{hit, minus}",
00117                                     45, 0., 45.,45, 0., 45.));
00118   myTrackHists1D.push_back(new TH1F("nHitInvalidTrack", "N_{hit, invalid}(track)", 5, 0., 5.));
00119   myTrackHists1D.push_back(new TH1F("r1Track", "r(1st hit);r [cm]", 20, 0., 20.));
00120   myTrackHists1D.push_back(new TH1F("phi1Track", "#phi(1st hit);#phi",
00121                                     30, -TMath::Pi(), TMath::Pi()));
00122   myTrackHists1D.push_back(new TH1F("y1Track", "y(1st hit);y [cm]", 40, -120., +120.));
00123   myTrackHists1D.push_back(new TH1F("z1Track", "z(1st hit);z [cm]", 20, -50, +50));
00124   myTrackHists1D.push_back(new TH1F("r1TrackSigned", "r(1st hit);r_{#pm} [cm]",
00125                                     40, -120., 120.));
00126   myTrackHists1D.push_back(new TH1F("z1TrackFull", "z(1st hit);z [cm]", 20, -300., +300.));
00127   myTrackHists1D.push_back(new TH1F("rLastTrack", "r(last hit);r [cm]", 20, 20., 120.));
00128   myTrackHists1D.push_back(new TH1F("phiLastTrack", "#phi(last hit);#phi",
00129                                     30, -TMath::Pi(), TMath::Pi()));
00130   myTrackHists1D.push_back(new TH1F("yLastTrack", "y(last hit);y [cm]", 40, -120., +120.));
00131   myTrackHists1D.push_back(new TH1F("zLastTrack", "z(last hit);z [cm]", 30, -300., +300.));
00132   myTrackHists1D.push_back(new TH1F("chi2PerNdf", "#chi^{2}/ndf;#chi^{2}/ndf", 500, 0., 50.));
00133   myTrackHists1D.push_back(new TH1F("impParZ", "impact parameter in z", 20, -20., 20.));
00134   myTrackHists1D.push_back(new TH1F("impParErrZ", "error of impact parameter in z",
00135                                     40, 0., 0.06));  
00136   myTrackHists1D.push_back(new TH1F("impParRphi", "impact parameter in r#phi", 51, -0.05, .05));
00137   myTrackHists1D.push_back(new TH1F("impParErrRphi", "error of impact parameter in r#phi",
00138                                     50, 0., 0.01));
00139 
00140   myTrackHists2D.push_back(new TH2F("rz1TrackFull", "rz(1st hit);z [cm]; r_{#pm} [cm]",
00141                                     40, -282., +282., 40, -115., +115.));
00142   myTrackHists2D.push_back(new TH2F("xy1Track", "xy(1st hit);x [cm]; y [cm]",
00143                                     40, -115., +115., 40, -115., +115.));
00144   
00145   TDirectory *dirTracks = directory->mkdir("trackHists", "input tracks");
00146   this->addToDirectory(myTrackHists1D, dirTracks);
00147   this->addToDirectory(myTrackHists2D, dirTracks);
00148 
00149 // used track 
00150   myUsedTrackHists1D = this->cloneHists(myTrackHists1D, "used", " (used tracks)");
00151   myUsedTrackHists2D = this->cloneHists(myTrackHists2D, "used", " (used tracks)");
00152   // must be after clone: index in vector!
00153   myUsedTrackHists1D.push_back(new TH1F("usedHitsX", "n(x-hits) transferred to pede;n(x-hits)", nHits, 0., nHits));
00154   myUsedTrackHists1D.push_back(new TH1F("usedHitsY", "n(y-hits) transferred to pede;n(y-hits)", 10, 0., 10));
00155 
00156   TDirectory *dirUsedTracks = directory->mkdir("usedTrackHists", "used tracks");
00157   this->addToDirectory(myUsedTrackHists1D, dirUsedTracks);
00158   this->addToDirectory(myUsedTrackHists2D, dirUsedTracks);
00159 
00160 // ReferenceTrajectory
00161   myTrajectoryHists1D.push_back(new TH1F("validRefTraj", "validity of ReferenceTrajectory",
00162                                          2, 0., 2.));
00163 
00164   myTrajectoryHists2D.push_back(new TProfile2D("profCorr",
00165                                                "mean of |#rho|, #rho#neq0;hit x/y;hit x/y;",
00166                                                2*nHits, 0., nHits, 2*nHits, 0., nHits));
00167   myTrajectoryHists2D.push_back
00168     (new TProfile2D("profCorrOffXy", "mean of |#rho|, #rho#neq0, no xy_{hit};hit x/y;hit x/y;",
00169                     2*nHits, 0., nHits, 2*nHits, 0., nHits));
00170 
00171   myTrajectoryHists2D.push_back(new TH2F("hitCorrOffBlock",
00172                                          "hit correlations: off-block-diagonals;N(hit);#rho",
00173                                          2*nHits, 0., nHits, 81, -.06, .06));
00174   TArrayD logBins(102);
00175   this->equidistLogBins(logBins.GetArray(), logBins.GetSize()-1, 1.E-11, .1);
00176   myTrajectoryHists2D.push_back(new TH2F("hitCorrOffBlockLog",
00177                                          "hit correlations: off-block-diagonals;N(hit);|#rho|",
00178                                          2*nHits, 0., nHits, logBins.GetSize()-1, logBins.GetArray()));
00179 
00180   myTrajectoryHists2D.push_back(new TH2F("hitCorrXy", "hit correlations: xy;N(hit);#rho",
00181                                          nHits, 0., nHits, 81, -.5, .5));
00182   myTrajectoryHists2D.push_back
00183     (new TH2F("hitCorrXyValid", "hit correlations: xy, 2D-det.;N(hit);#rho",
00184               nHits, 0., nHits, 81, -.02, .02));
00185   this->equidistLogBins(logBins.GetArray(), logBins.GetSize()-1, 1.E-10, .5);
00186   myTrajectoryHists2D.push_back(new TH2F("hitCorrXyLog", "hit correlations: xy;N(hit);|#rho|",
00187                                          nHits, 0., nHits, logBins.GetSize()-1, logBins.GetArray()));
00188   myTrajectoryHists2D.push_back
00189     (new TH2F("hitCorrXyLogValid", "hit correlations: xy, 2D-det.;N(hit);|#rho|",
00190               nHits, 0., nHits, logBins.GetSize()-1, logBins.GetArray()));
00191 
00192 
00193   myTrajectoryHists1D.push_back(new TH1F("measLocX", "local x measurements;x", 101, -6., 6.));
00194   myTrajectoryHists1D.push_back(new TH1F("measLocY", "local y measurements, 2D-det.;y",
00195                                          101, -10., 10.));
00196   myTrajectoryHists1D.push_back(new TH1F("trajLocX", "local x trajectory;x", 101, -6., 6.));
00197   myTrajectoryHists1D.push_back(new TH1F("trajLocY", "local y trajectory, 2D-det.;y",
00198                                          101, -10., 10.));
00199 
00200   myTrajectoryHists1D.push_back(new TH1F("residLocX", "local x residual;#Deltax", 101, -.75, .75));
00201   myTrajectoryHists1D.push_back(new TH1F("residLocY", "local y residual, 2D-det.;#Deltay",
00202                                          101, -2., 2.));
00203   myTrajectoryHists1D.push_back(new TH1F("reduResidLocX", "local x reduced residual;#Deltax/#sigma",
00204                                          101, -20., 20.));
00205   myTrajectoryHists1D.push_back
00206     (new TH1F("reduResidLocY", "local y reduced residual, 2D-det.;#Deltay/#sigma", 101, -20., 20.));
00207 
00208   // 2D vs. hit
00209   myTrajectoryHists2D.push_back(new TH2F("measLocXvsHit", "local x measurements;hit;x", 
00210                                          nHits, 0., nHits, 101, -6., 6.));
00211   myTrajectoryHists2D.push_back(new TH2F("measLocYvsHit", "local y measurements, 2D-det.;hit;y",
00212                                          nHits, 0., nHits, 101, -10., 10.));
00213   myTrajectoryHists2D.push_back(new TH2F("trajLocXvsHit", "local x trajectory;hit;x",
00214                                          nHits, 0., nHits, 101, -6., 6.));
00215   myTrajectoryHists2D.push_back(new TH2F("trajLocYvsHit", "local y trajectory, 2D-det.;hit;y",
00216                                          nHits, 0., nHits,  101, -10., 10.));
00217 
00218   myTrajectoryHists2D.push_back(new TH2F("residLocXvsHit", "local x residual;hit;#Deltax",
00219                                          nHits, 0., nHits, 101, -.75, .75));
00220   myTrajectoryHists2D.push_back(new TH2F("residLocYvsHit", "local y residual, 2D-det.;hit;#Deltay",
00221                                          nHits, 0., nHits, 101, -1., 1.));
00222   myTrajectoryHists2D.push_back
00223     (new TH2F("reduResidLocXvsHit", "local x reduced residual;hit;#Deltax/#sigma",
00224               nHits, 0., nHits, 101, -20., 20.));
00225   myTrajectoryHists2D.push_back
00226     (new TH2F("reduResidLocYvsHit", "local y reduced residual, 2D-det.;hit;#Deltay/#sigma",
00227               nHits, 0., nHits, 101, -20., 20.));
00228 
00229 
00230   myTrajectoryHists2D.push_back(new TProfile2D("profDerivatives",
00231                                                "mean derivatives;hit x/y;parameter;",
00232                                                2*nHits, 0., nHits, 10, 0., 10.));
00233 
00234   myTrajectoryHists2D.push_back
00235     (new TH2F("derivatives", "derivative;parameter;#partial(x/y)_{local}/#partial(param)",
00236               10, 0., 10., 101, -20., 20.));
00237   this->equidistLogBins(logBins.GetArray(), logBins.GetSize()-1, 1.E-12, 100.);
00238   myTrajectoryHists2D.push_back
00239     (new TH2F("derivativesLog", "|derivative|;parameter;|#partial(x/y)_{local}/#partial(param)|",
00240               10, 0., 10., logBins.GetSize()-1, logBins.GetArray()));
00241   myTrajectoryHists2D.push_back
00242     (new TH2F("derivativesVsPhi", 
00243               "derivatives vs. #phi;#phi(geomDet);#partial(x/y)_{local}/#partial(param)",
00244               50, -TMath::Pi(), TMath::Pi(), 101, -300., 300.));
00245   //  myTrajectoryHists2D.back()->SetBit(TH1::kCanRebin);
00246 
00247   TDirectory *dirTraject = directory->mkdir("refTrajectoryHists", "ReferenceTrajectory's");
00248   this->addToDirectory(myTrajectoryHists2D, dirTraject);
00249   this->addToDirectory(myTrajectoryHists1D, dirTraject);
00250 
00251 // derivatives hists
00252   myDerivHists2D.push_back
00253     (new TH2F("localDerivsPar","local derivatives vs. paramater;parameter;#partial/#partial(param)",
00254               6, 0., 6., 101, -200., 200.));
00255   myDerivHists2D.push_back
00256     (new TH2F("localDerivsPhi","local derivatives vs. #phi(det);#phi(det);#partial/#partial(param)",
00257               51, -TMath::Pi(), TMath::Pi(), 101, -150., 150.));
00258   this->equidistLogBins(logBins.GetArray(), logBins.GetSize()-1, 1.E-13, 150.);
00259   myDerivHists2D.push_back
00260     (new TH2F("localDerivsParLog",
00261               "local derivatives (#neq 0) vs. parameter;parameter;|#partial/#partial(param)|",
00262               6, 0., 6., logBins.GetSize()-1, logBins.GetArray()));
00263   myDerivHists2D.push_back
00264     (new TH2F("localDerivsPhiLog",
00265               "local derivatives (#neq 0) vs. #phi(det);#phi(det);|#partial/#partial(param)|",
00266               51, -TMath::Pi(), TMath::Pi(), logBins.GetSize()-1, logBins.GetArray()));
00267   const unsigned int maxParNum = PedeLabelerBase::theMaxNumParam;
00268   myDerivHists2D.push_back
00269     (new TH2F("globalDerivsPar",
00270               "global derivatives vs. paramater;parameter;#partial/#partial(param)",
00271               maxParNum, 0., maxParNum, 100, -200, 200));
00272   myDerivHists2D.push_back
00273     (new TH2F("globalDerivsPhi",
00274               "global derivatives vs. #phi(det);#phi(det);#partial/#partial(param)",
00275               51, -TMath::Pi(), TMath::Pi(), 102, -1.02, 1.02));
00276   this->equidistLogBins(logBins.GetArray(), logBins.GetSize()-1, 1.E-7, 5.e2);
00277   myDerivHists2D.push_back
00278     (new TH2F("globalDerivsParLog",
00279               "global derivatives (#neq 0) vs. parameter;parameter;|#partial/#partial(param)|",
00280               maxParNum, 0., maxParNum, logBins.GetSize()-1, logBins.GetArray()));
00281   myDerivHists2D.push_back
00282     (new TH2F("globalDerivsPhiLog",
00283               "global derivatives (#neq 0) vs. #phi(det);#phi(det);|#partial/#partial(param)|",
00284               51, -TMath::Pi(), TMath::Pi(), logBins.GetSize()-1, logBins.GetArray()));
00285 //   this->equidistLogBins(logBins.GetArray(), logBins.GetSize()-1, 1.e-40, 1.E-35);
00286 //   myDerivHists2D.push_back
00287 //     (new TH2F("globalDerivsPhiLog2",
00288 //               "global derivatives (#neq 0) vs. #phi(det);#phi(det);|#partial/#partial(param)|",
00289 //               51, -TMath::Pi(), TMath::Pi(), logBins.GetSize()-1, logBins.GetArray()));
00290 
00291   TDirectory *dirDerivs = directory->mkdir("derivatives", "derivatives etc.");
00292   this->addToDirectory(myDerivHists2D, dirDerivs);
00293 
00294 // residual hists
00295   myResidHists2D.push_back(new TH2F("residPhi","residuum vs. #phi(det);#phi(det);residuum[cm]",
00296                                     51, -TMath::Pi(), TMath::Pi(), 101, -.5, .5));
00297   myResidHists2D.push_back(new TH2F("sigmaPhi","#sigma vs. #phi(det);#phi(det);#sigma[cm]",
00298                                     51, -TMath::Pi(), TMath::Pi(), 101, .0, .1));
00299   myResidHists2D.push_back(new TH2F("reduResidPhi",
00300                                     "residuum/#sigma vs. #phi(det);#phi(det);residuum/#sigma",
00301                                     51, -TMath::Pi(), TMath::Pi(), 101, -14., 14.));
00302   
00303 //   myResidHists2D.push_back(new TProfile2D("residXProfXy",
00304 //                                        "mean |residuum| (u);x [cm];y [cm];#LT|residuum|#GT",
00305 //                                        25, -110., 110., 25, -110., 110.));
00306 //   myResidHists2D.push_back(new TProfile2D("residXProfZr",
00307 //                                        "mean |residuum| (u);z [cm];r_{#pm} [cm];#LT|residuum|#GT",
00308 //                                        25, -275., 275., 25, -110., 110.));
00309 //   myResidHists2D.push_back(new TProfile2D("residYProfXy",
00310 //                                        "mean |residuum| (v);x [cm];y [cm];#LT|residuum|#GT",
00311 //                                        25, -110., 110., 25, -110., 110.));
00312 //   myResidHists2D.push_back(new TProfile2D("residYProfZr",
00313 //                                        "mean |residuum| (v);z [cm];r_{#pm} [cm];#LT|residuum|#GT",
00314 //                                        25, -275., 275., 25, -110., 110.));
00315 //   myResidHists2D.push_back(new TProfile2D("reduResidXProfXy",
00316 //                                        "mean |residuum/#sigma| (u);x [cm];y [cm];#LT|res./#sigma|#GT",
00317 //                                        25, -110., 110., 25, -110., 110.));
00318 //   myResidHists2D.push_back(new TProfile2D("reduResidXProfZr",
00319 //                                        "mean |residuum/#sigma| (u);z [cm];r_{#pm} [cm];#LT|res./#sigma|#GT",
00320 //                                        25, -275., 275., 25, -110., 110.));
00321 //   myResidHists2D.push_back(new TProfile2D("reduResidYProfXy",
00322 //                                        "mean |residuum/#sigma| (v);x [cm];y [cm];#LT|res./#sigma|#GT",
00323 //                                        25, -110., 110., 25, -110., 110.));
00324 //   myResidHists2D.push_back(new TProfile2D("reduResidYProfZr",
00325 //                                        "mean |residuum/#sigma| (v);z [cm];r_{#pm} [cm];#LT|res./#sigma|#GT",
00326 //                                        25, -275., 275., 25, -110., 110.));
00327 //   myResidHists2D.push_back(new TProfile2D("sigmaXProfXy",
00328 //                                        "mean sigma (u);x [cm];y [cm];#LT#sigma#GT",
00329 //                                        25, -110., 110., 25, -110., 110.));
00330 //   myResidHists2D.push_back(new TProfile2D("sigmaXProfZr",
00331 //                                        "mean sigma (u);z [cm];r_{#pm} [cm];#LT#sigma#GT",
00332 //                                        25, -275., 275., 25, -110., 110.));
00333 //   myResidHists2D.push_back(new TProfile2D("sigmaYProfXy",
00334 //                                        "mean sigma (v);x [cm];y [cm];#LT#sigma#GT",
00335 //                                        25, -110., 110., 25, -110., 110.));
00336 //   myResidHists2D.push_back(new TProfile2D("sigmaYProfZr",
00337 //                                        "mean sigma (v);z [cm];r_{#pm} [cm];#LT#sigma#GT",
00338 //                                        25, -275., 275., 25, -110., 110.));
00339   
00340   TDirectory *dirResid = directory->mkdir("residuals", "hit residuals, sigma,...");
00341   this->addToDirectory(myResidHists2D, dirResid);
00342 
00343 // Residuum, hit sigma and res./sigma for all sensor/track angles and separated for large and
00344 // small angles with respect to the sensor normal in sensitive direction. 
00345   // Here for x-measurements:
00346   std::vector<TH1*> allResidHistsX;
00347   allResidHistsX.push_back(new TH1F("resid", "hit residuals;residuum [cm]", 101,-.5,.5));//51,-.05, .05));
00348   //allResidHistsX.back()->SetBit(TH1::kCanRebin);
00349   allResidHistsX.push_back(new TH1F("sigma", "hit uncertainties;#sigma [cm]", 100,0.,1.));//50, 0., .02));
00350   //allResidHistsX.back()->SetBit(TH1::kCanRebin);
00351   allResidHistsX.push_back(new TH1F("reduResid", "reduced hit residuals;res./#sigma",
00352                                     101, -10., 10.));//51, -3., 3.));
00353   //  allResidHistsX.back()->SetBit(TH1::kCanRebin);
00354   allResidHistsX.push_back(new TH1F("angle", "#phi_{tr} wrt normal (sens. plane);#phi_{n}^{sens}",
00355                                     50, 0., TMath::PiOver2()));
00356   allResidHistsX.push_back(new TH2F("residVsAngle",
00357                                     "residuum vs. #phi_{tr} wrt normal (sens. plane);#phi_{n}^{sens};residuum [cm]",
00358                                     50, 0., TMath::PiOver2(), 51, -1., 1.));
00359   this->equidistLogBins(logBins.GetArray(), logBins.GetSize()-1, 1.E-6, 100.);
00360   allResidHistsX.push_back(new TH2F("sigmaVsAngle",
00361                                     "#sigma vs. #phi_{tr} wrt normal (sens. plane);#phi_{n}^{sens};#sigma [cm]",
00362                                     50, 0., TMath::PiOver2(), logBins.GetSize()-1, logBins.GetArray()));
00363   allResidHistsX.push_back(new TH2F("reduResidVsAngle",
00364                                     "reduced residuum vs. #phi_{tr} wrt normal (sens. plane);#phi_{n}^{sens};res./#sigma",
00365                                     50, 0., TMath::PiOver2(), 51, -15., 15.));
00366 
00367   allResidHistsX.push_back(new TH1F("residGt45",
00368                                     "hit residuals (#phi_{n}^{sens}>45#circ);residuum [cm]",
00369                                     101, -.5, .5));//51, -.05, .05));
00370   // allResidHistsX.back()->SetBit(TH1::kCanRebin);
00371   allResidHistsX.push_back(new TH1F("sigmaGt45",
00372                                     "hit uncertainties(#phi_{n}^{sens}>45#circ);#sigma [cm]",
00373                                      100, 0., 1.));//50, 0., .02));
00374   // allResidHistsX.back()->SetBit(TH1::kCanRebin);
00375   allResidHistsX.push_back(new TH1F("reduResidGt45",
00376                                     "reduced hit residuals(#phi_{n}^{sens}>45#circ);res./#sigma",
00377                                     101, -10., 10.));//51,-3.,3.));
00378   // allResidHistsX.back()->SetBit(TH1::kCanRebin);
00379   allResidHistsX.push_back(new TH1F("residLt45",
00380                                     "hit residuals (#phi_{n}^{sens}<45#circ);residuum [cm]",
00381                                     101, -.5, .5));//51, -.15, .15));
00382   // allResidHistsX.back()->SetBit(TH1::kCanRebin);
00383   allResidHistsX.push_back(new TH1F("sigmaLt45",
00384                                     "hit uncertainties(#phi_{n}^{sens}<45#circ);#sigma [cm]",
00385                                     100, 0., 1.));//50, 0., .01));
00386   // allResidHistsX.back()->SetBit(TH1::kCanRebin);
00387   allResidHistsX.push_back(new TH1F("reduResidLt45",
00388                                     "reduced hit residuals(#phi_{n}^{sens}<45#circ);res./#sigma",
00389                                     101, -10., 10.));//51,-3.,3.));
00390   myResidHistsVec1DX.push_back(allResidHistsX); // at [0] for all subdets together...
00391   // ... then separately - indices/order like DetId.subdetId() in tracker:
00392   myResidHistsVec1DX.push_back(this->cloneHists(allResidHistsX, "TPB", " x-coord. in pixel barrel"));
00393   myResidHistsVec1DX.push_back(this->cloneHists(allResidHistsX, "TPE", " x-coord. in pixel discs"));
00394   myResidHistsVec1DX.push_back(this->cloneHists(allResidHistsX, "TIB", " x-coord. in TIB"));
00395   myResidHistsVec1DX.push_back(this->cloneHists(allResidHistsX, "TID", " x-coord. in TID"));
00396   myResidHistsVec1DX.push_back(this->cloneHists(allResidHistsX, "TOB", " x-coord. in TOB"));
00397   myResidHistsVec1DX.push_back(this->cloneHists(allResidHistsX, "TEC", " x-coord. in TEC"));
00398   // finally, differential in hit number (but subdet independent)
00399   for (unsigned int iHit = 0; iHit < 30; ++iHit) { // 4: for each hit fill only angle independent plots
00400     for (unsigned int iHist = 0; iHist < 4 && iHist < allResidHistsX.size(); ++iHist) {
00401       TH1 *h = allResidHistsX[iHist];
00402       myResidHitHists1DX.push_back(static_cast<TH1*>(h->Clone(Form("%s_%d", h->GetName(), iHit))));
00403       myResidHitHists1DX.back()->SetTitle(Form("%s, hit %d", h->GetTitle(), iHit));
00404     }
00405   }
00406 
00407   TDirectory *dirResidX =
00408     (dirResid ? dirResid : directory)->mkdir("X", "hit residuals etc. for x measurements");
00409   this->addToDirectory(myResidHitHists1DX, dirResidX);
00410   for (std::vector<std::vector<TH1*> >::iterator vecIter = myResidHistsVec1DX.begin(),
00411          vecIterEnd = myResidHistsVec1DX.end(); vecIter != vecIterEnd; ++vecIter) {
00412     this->addToDirectory(*vecIter, dirResidX);
00413   }
00414 
00415   // Now clone the same as above for y-ccordinate:
00416   myResidHistsVec1DY.push_back(this->cloneHists(allResidHistsX, "", " y-coord."));// all subdets
00417   std::vector<TH1*> &yHists1D = allResidHistsX;//myResidHistsVec1DY.back(); crashes? ROOT?
00418   myResidHistsVec1DY.push_back(this->cloneHists(yHists1D, "TPB", " y-coord. in pixel barrel"));
00419   myResidHistsVec1DY.push_back(this->cloneHists(yHists1D, "TPE", " y-coord. in pixel discs"));
00420   myResidHistsVec1DY.push_back(this->cloneHists(yHists1D, "TIB", " y-coord. in TIB"));
00421   myResidHistsVec1DY.push_back(this->cloneHists(yHists1D, "TID", " y-coord. in TID"));
00422   myResidHistsVec1DY.push_back(this->cloneHists(yHists1D, "TOB", " y-coord. in TOB"));
00423   myResidHistsVec1DY.push_back(this->cloneHists(yHists1D, "TEC", " y-coord. in TEC"));
00424   myResidHitHists1DY = this->cloneHists(myResidHitHists1DX, "", " y-coord.");// diff. in nHit
00425 
00426   TDirectory *dirResidY = 
00427     (dirResid ? dirResid : directory)->mkdir("Y", "hit residuals etc. for y measurements");
00428   this->addToDirectory(myResidHitHists1DY, dirResidY);
00429   for (std::vector<std::vector<TH1*> >::iterator vecIter = myResidHistsVec1DY.begin(),
00430          vecIterEnd = myResidHistsVec1DY.end(); vecIter != vecIterEnd; ++vecIter) {
00431     this->addToDirectory(*vecIter, dirResidY);
00432   }
00433 
00434   // farme-to-frame derivatives
00435   myFrame2FrameHists2D.push_back(new TProfile2D("frame2frame",
00436                                                 "mean frame to frame derivatives;col;row",
00437                                                 6, 0., 6., 6, 0., 6.));
00438   myFrame2FrameHists2D.push_back(new TProfile2D("frame2frameAbs",
00439                                                 "mean |frame to frame derivatives|, #neq0;col;row",
00440                                                 6, 0., 6., 6, 0., 6.));
00441 
00442   this->equidistLogBins(logBins.GetArray(), logBins.GetSize()-1, 1.E-36, 100.);
00443   for (unsigned int i = 0; i < 6; ++i) {
00444     for (unsigned int j = 0; j < 6; ++j) {
00445       myFrame2FrameHists2D.push_back
00446         (new TH2F(Form("frame2framePhi%d%d", i, j),
00447                   Form("frame to frame derivatives, %d%d;#phi(aliDet);deriv",i,j),
00448                   51, -TMath::Pi(), TMath::Pi(), 10, 0., 1.));
00449       myFrame2FrameHists2D.back()->SetBit(TH1::kCanRebin);
00450       myFrame2FrameHists2D.push_back
00451         (new TH2F(Form("frame2frameR%d%d", i, j),
00452                   Form("frame to frame derivatives, %d%d;r(aliDet);deriv",i,j),
00453                   51, 0., 110., 10, 0., 1.));
00454       myFrame2FrameHists2D.back()->SetBit(TH1::kCanRebin);
00455 
00456       myFrame2FrameHists2D.push_back
00457         (new TH2F(Form("frame2framePhiLog%d%d", i, j),
00458                   Form("frame to frame |derivatives|, %d%d, #neq0;#phi(aliDet);deriv",i,j),
00459                   51, -TMath::Pi(), TMath::Pi(), logBins.GetSize()-1, logBins.GetArray()));
00460       myFrame2FrameHists2D.push_back
00461         (new TH2F(Form("frame2frameRLog%d%d", i, j),
00462                   Form("frame to frame |derivatives|, %d%d, #neq0;r(aliDet);deriv",i,j),
00463                   51, 0., 110., logBins.GetSize()-1, logBins.GetArray()));
00464     }
00465   }
00466 
00467   TDirectory *dirF2f = directory->mkdir("frame2FrameHists", "derivatives etc.");
00468   this->addToDirectory(myFrame2FrameHists2D, dirF2f);
00469 
00470 
00471   myCorrHists.push_back(new TH1F("xyCorrTPB", "#rho_{xy} in pixel barrel", 50, -.5, .5));
00472   myCorrHists.push_back(new TH1F("xyCorrTPE", "#rho_{xy} in forward pixel", 50, -.5, .5));
00473   myCorrHists.push_back(new TH1F("xyCorrTIB", "#rho_{xy} in TIB", 50, -.5, .5));
00474   myCorrHists.push_back(new TH1F("xyCorrTID", "#rho_{xy} in TID", 50, -1., 1.));
00475   myCorrHists.push_back(new TH1F("xyCorrTOB", "#rho_{xy} in TOB", 50, -.5, .5));
00476   myCorrHists.push_back(new TH1F("xyCorrTEC", "#rho_{xy} in TEC", 50, -1., 1.));
00477   TDirectory *dirCorr = directory->mkdir("hitCorrelationHists", "correlations");
00478   this->addToDirectory( myCorrHists, dirCorr);
00479 
00480   // PXB Survey
00481   myPxbSurveyHists.push_back(new TH1F("PxbSurvChi2_lo", "#chi^{2} from PXB survey", 25, 0, 25));
00482   myPxbSurveyHists.push_back(new TH1F("PxbSurvChi2_md", "#chi^{2} from PXB survey", 25, 0, 500));
00483   myPxbSurveyHists.push_back(new TH1F("PxbSurvChi2_hi", "#chi^{2} from PXB survey", 25, 0, 50000));
00484   myPxbSurveyHists.push_back(new TH1F("PxbSurvChi2prob", "Math::Prob(#chi^{2},4) from PXB survey", 25, 0, 1));
00485   myPxbSurveyHists.push_back(new TH1F("PxbSurv_a0", "a_{0} from PXB survey", 100, -3000, 3000));
00486   myPxbSurveyHists.push_back(new TH1F("PxbSurv_a0_abs", "fabs(a_{0}) from PXB survey", 100, 0, 3000));
00487   myPxbSurveyHists.push_back(new TH1F("PxbSurv_a1", "a_{1} from PXB survey", 100, -6000, 6000));
00488   myPxbSurveyHists.push_back(new TH1F("PxbSurv_a1_abs", "fabs(a_{1}) from PXB survey", 100, 0, 6000));
00489   myPxbSurveyHists.push_back(new TH1F("PxbSurv_scale", "scale (#sqrt{a_{2}^{2}+a_{3}^{2}}) from PXB survey", 100, 0, 1500));
00490   myPxbSurveyHists.push_back(new TH1F("PxbSurv_phi", "angle(#atan{a_{3}/a_{4}}) from PXB survey", 100, -.05, .05));
00491   TDirectory *dirPxbSurvey = directory->mkdir("PxbSurveyHists", "PxbSurvey");
00492   this->addToDirectory( myPxbSurveyHists, dirPxbSurvey);
00493 
00494   oldDir->cd();
00495   return true;
00496 }
00497 
00498 
00499 //__________________________________________________________________
00500 bool MillePedeMonitor::equidistLogBins(double* bins, int nBins, 
00501                                         double first, double last) const
00502 {
00503   // Filling 'bins' with borders of 'nBins' bins between 'first' and 'last'
00504   // that are equidistant when viewed in log scale,
00505   // so 'bins' must have length nBins+1;
00506   // If 'first', 'last' or 'nBins' are not positive, failure is reported.
00507 
00508   if (nBins < 1 || first <= 0. || last <= 0.) return false;
00509 
00510   bins[0] = first;
00511   bins[nBins] = last;
00512   const double firstLog = TMath::Log10(bins[0]);
00513   const double lastLog  = TMath::Log10(bins[nBins]);
00514   for (int i = 1; i < nBins; ++i) {
00515     bins[i] = TMath::Power(10., firstLog + i*(lastLog-firstLog)/(nBins));
00516   }
00517 
00518   return true;
00519 }
00520 
00521 //__________________________________________________________________
00522 void MillePedeMonitor::fillTrack(const reco::Track *track)
00523 {
00524   this->fillTrack(track, myTrackHists1D, myTrackHists2D);
00525 }
00526 
00527 //__________________________________________________________________
00528 void MillePedeMonitor::fillUsedTrack(const reco::Track *track, unsigned int nHitX,
00529                                      unsigned int nHitY)
00530 {
00531   // these hist exist only for 'used track' hists:
00532   static const int iUsedX = this->GetIndex(myUsedTrackHists1D, "usedHitsX");
00533   myUsedTrackHists1D[iUsedX]->Fill(nHitX);
00534   static const int iUsedY = this->GetIndex(myUsedTrackHists1D, "usedHitsY");
00535   myUsedTrackHists1D[iUsedY]->Fill(nHitY);
00536 
00537   if (!track) return;
00538   this->fillTrack(track, myUsedTrackHists1D, myUsedTrackHists2D);
00539 }
00540 
00541 //__________________________________________________________________
00542 void MillePedeMonitor::fillTrack(const reco::Track *track, std::vector<TH1*> &trackHists1D,
00543                                  std::vector<TH2*> &trackHists2D)
00544 {
00545   if (!track) return;
00546 
00547   const reco::TrackBase::Vector p(track->momentum());
00548 
00549   static const int iPtLog = this->GetIndex(trackHists1D, "ptTrackLogBins");
00550   trackHists1D[iPtLog]->Fill(track->pt());
00551   static const int iPt = this->GetIndex(trackHists1D, "ptTrack");
00552   trackHists1D[iPt]->Fill(track->pt());
00553   static const int iP = this->GetIndex(trackHists1D, "pTrack");
00554   trackHists1D[iP]->Fill(p.R());
00555   static const int iEta = this->GetIndex(trackHists1D, "etaTrack");
00556   trackHists1D[iEta]->Fill(p.Eta());
00557   static const int iTheta = this->GetIndex(trackHists1D, "thetaTrack");
00558   trackHists1D[iTheta]->Fill(p.Theta());
00559   static const int iPhi = this->GetIndex(trackHists1D, "phiTrack");
00560   trackHists1D[iPhi]->Fill(p.Phi());
00561 
00562   static const int iNhit = this->GetIndex(trackHists1D, "nHitTrack");
00563   trackHists1D[iNhit]->Fill(track->numberOfValidHits());
00564   static const int iNhitInvalid = this->GetIndex(trackHists1D, "nHitInvalidTrack");
00565   trackHists1D[iNhitInvalid]->Fill(track->numberOfLostHits());
00566 
00567   int nhitinTIB = 0, nhitinTOB = 0, nhitinTID = 0;
00568   int nhitinTEC = 0, nhitinBPIX = 0, nhitinFPIX = 0, nhitinPIXEL=0;
00569   int nhitinENDCAP = 0, nhitinENDCAPplus = 0, nhitinENDCAPminus = 0;
00570   int nhitinTIDplus = 0, nhitinTIDminus = 0;
00571   int nhitinFPIXplus = 0, nhitinFPIXminus = 0;
00572   int nhitinTECplus = 0, nhitinTECminus = 0;
00573   unsigned int thishit = 0;
00574 
00575   for (trackingRecHit_iterator iHit = track->recHitsBegin(); iHit != track->recHitsEnd(); ++iHit) {
00576     thishit++;
00577     const DetId detId((*iHit)->geographicalId());
00578     const int subdetId = detId.subdetId(); 
00579 
00580     if (!(*iHit)->isValid()) continue; // only real hits count as in track->numberOfValidHits()
00581     if (detId.det() != DetId::Tracker) {
00582       edm::LogError("DetectorMismatch") << "@SUB=MillePedeMonitor::fillTrack"
00583                                         << "DetId.det() != DetId::Tracker (=" << DetId::Tracker
00584                                         << "), but " << detId.det() << ".";
00585     }
00586 
00587     if      (SiStripDetId::TIB == subdetId) ++nhitinTIB;
00588     else if (SiStripDetId::TOB == subdetId) ++nhitinTOB;
00589     else if (SiStripDetId::TID == subdetId) {
00590       ++nhitinTID;
00591       ++nhitinENDCAP;
00592       TIDDetId tidId(detId);
00593       if (tidId.isZMinusSide()) {
00594         ++nhitinTIDminus;
00595         ++nhitinENDCAPminus;
00596       }
00597       else if (tidId.isZPlusSide()) {
00598         ++nhitinTIDplus;
00599         ++nhitinENDCAPplus;
00600       }
00601     }
00602     else if (SiStripDetId::TEC == subdetId) {
00603       ++nhitinTEC;
00604       ++nhitinENDCAP;
00605       TECDetId tecId(detId);
00606       if (tecId.isZMinusSide()) {
00607         ++nhitinTECminus;
00608         ++nhitinENDCAPminus;
00609       }
00610       else if (tecId.isZPlusSide()) {
00611         ++nhitinTECplus;
00612         ++nhitinENDCAPplus;
00613       }
00614     }
00615     else if (            kBPIX == subdetId) {++nhitinBPIX;++nhitinPIXEL;}
00616     else if (            kFPIX == subdetId) {
00617       ++nhitinFPIX;
00618       ++nhitinPIXEL;
00619       PXFDetId fpixId(detId);
00620       if (fpixId.side()==1) ++nhitinFPIXminus;
00621       else if (fpixId.side()==2) ++nhitinFPIXplus;
00622     }
00623 
00624   } // end loop on hits
00625 
00626   static const int iNhit01 = this->GetIndex(trackHists1D, "nHitBPIXTrack");
00627   trackHists1D[iNhit01]->Fill(nhitinBPIX);
00628   static const int iNhit02 = this->GetIndex(trackHists1D, "nHitFPIXplusTrack");
00629   trackHists1D[iNhit02]->Fill(nhitinFPIXplus);
00630   static const int iNhit03 = this->GetIndex(trackHists1D, "nHitFPIXminusTrack");
00631   trackHists1D[iNhit03]->Fill(nhitinFPIXminus);
00632   static const int iNhit04 = this->GetIndex(trackHists1D, "nHitFPIXTrack");
00633   trackHists1D[iNhit04]->Fill(nhitinFPIX);
00634   static const int iNhit05 = this->GetIndex(trackHists1D, "nHitPIXELTrack");
00635   trackHists1D[iNhit05]->Fill(nhitinPIXEL);
00636   static const int iNhit06 = this->GetIndex(trackHists1D, "nHitTIBTrack");
00637   trackHists1D[iNhit06]->Fill(nhitinTIB);
00638   static const int iNhit07 = this->GetIndex(trackHists1D, "nHitTOBTrack");
00639   trackHists1D[iNhit07]->Fill(nhitinTOB);
00640   static const int iNhit08 = this->GetIndex(trackHists1D, "nHitTIDplusTrack");
00641   trackHists1D[iNhit08]->Fill(nhitinTIDplus);
00642   static const int iNhit09 = this->GetIndex(trackHists1D, "nHitTIDminusTrack");
00643   trackHists1D[iNhit09]->Fill(nhitinTIDminus);
00644   static const int iNhit10 = this->GetIndex(trackHists1D, "nHitTIDTrack");
00645   trackHists1D[iNhit10]->Fill(nhitinTID);
00646   static const int iNhit11 = this->GetIndex(trackHists1D, "nHitTECplusTrack");
00647   trackHists1D[iNhit11]->Fill(nhitinTECplus);
00648   static const int iNhit12 = this->GetIndex(trackHists1D, "nHitTECminusTrack");
00649   trackHists1D[iNhit12]->Fill(nhitinTECminus);
00650   static const int iNhit13 = this->GetIndex(trackHists1D, "nHitTECTrack");
00651   trackHists1D[iNhit13]->Fill(nhitinTEC);
00652   static const int iNhit14 = this->GetIndex(trackHists1D, "nHitENDCAPplusTrack");
00653   trackHists1D[iNhit14]->Fill(nhitinENDCAPplus);
00654   static const int iNhit15 = this->GetIndex(trackHists1D, "nHitENDCAPminusTrack");
00655   trackHists1D[iNhit15]->Fill(nhitinENDCAPminus);
00656   static const int iNhit16 = this->GetIndex(trackHists1D, "nHitENDCAPTrack");
00657   trackHists1D[iNhit16]->Fill(nhitinENDCAP);
00658   static const int iNhit17 = this->GetIndex(trackHists2D, "nHitENDCAPTrackMinusVsPlus");
00659   trackHists2D[iNhit17]->Fill(nhitinENDCAPplus,nhitinENDCAPminus);
00660 
00661   if (track->innerOk()) {
00662     const reco::TrackBase::Point firstPoint(track->innerPosition());
00663     static const int iR1 = this->GetIndex(trackHists1D, "r1Track");
00664     trackHists1D[iR1]->Fill(firstPoint.Rho());
00665     const double rSigned1 = (firstPoint.y() > 0 ? firstPoint.Rho() : -firstPoint.Rho());
00666     static const int iR1Signed = this->GetIndex(trackHists1D, "r1TrackSigned");
00667     trackHists1D[iR1Signed]->Fill(rSigned1);
00668     static const int iZ1 = this->GetIndex(trackHists1D, "z1Track");
00669     trackHists1D[iZ1]->Fill(firstPoint.Z());
00670     static const int iZ1Full = this->GetIndex(trackHists1D, "z1TrackFull");
00671     trackHists1D[iZ1Full]->Fill(firstPoint.Z());
00672     static const int iY1 = this->GetIndex(trackHists1D, "y1Track");
00673     trackHists1D[iY1]->Fill(firstPoint.Y());
00674     static const int iPhi1 = this->GetIndex(trackHists1D, "phi1Track");
00675     trackHists1D[iPhi1]->Fill(firstPoint.phi());
00676     static const int iRz1Full = this->GetIndex(trackHists2D, "rz1TrackFull");
00677     trackHists2D[iRz1Full]->Fill(firstPoint.Z(), rSigned1);
00678     static const int iXy1 = this->GetIndex(trackHists2D, "xy1Track");
00679     trackHists2D[iXy1]->Fill(firstPoint.X(), firstPoint.Y());
00680   }
00681 
00682   if (track->outerOk()) {
00683     const reco::TrackBase::Point lastPoint(track->outerPosition());
00684     static const int iRlast = this->GetIndex(trackHists1D, "rLastTrack");
00685     trackHists1D[iRlast]->Fill(lastPoint.Rho());
00686     static const int iZlast = this->GetIndex(trackHists1D, "zLastTrack");
00687     trackHists1D[iZlast]->Fill(lastPoint.Z());
00688     static const int iYlast = this->GetIndex(trackHists1D, "yLastTrack");
00689     trackHists1D[iYlast]->Fill(lastPoint.Y());
00690     static const int iPhiLast = this->GetIndex(trackHists1D, "phiLastTrack");
00691     trackHists1D[iPhiLast]->Fill(lastPoint.phi());
00692   }
00693 
00694   static const int iChi2Ndf = this->GetIndex(trackHists1D, "chi2PerNdf");
00695   trackHists1D[iChi2Ndf]->Fill(track->normalizedChi2());
00696 
00697   static const int iImpZ = this->GetIndex(trackHists1D, "impParZ");
00698   trackHists1D[iImpZ]->Fill(track->dz());
00699   static const int iImpZerr = this->GetIndex(trackHists1D, "impParErrZ");
00700   trackHists1D[iImpZerr]->Fill(track->dzError());
00701 
00702 
00703   static const int iImpRphi = this->GetIndex(trackHists1D, "impParRphi");
00704   trackHists1D[iImpRphi]->Fill(track->d0());
00705   static const int iImpRphiErr = this->GetIndex(trackHists1D, "impParErrRphi");
00706   trackHists1D[iImpRphiErr]->Fill(track->d0Error());
00707 
00708 }
00709 
00710 //__________________________________________________________________
00711 void MillePedeMonitor::fillRefTrajectory(const ReferenceTrajectoryBase::ReferenceTrajectoryPtr &refTrajPtr)
00712 {
00713 
00714   static const int iValid = this->GetIndex(myTrajectoryHists1D, "validRefTraj");
00715   if (refTrajPtr->isValid()) {
00716     myTrajectoryHists1D[iValid]->Fill(1.);
00717   } else {
00718     myTrajectoryHists1D[iValid]->Fill(0.);
00719     return;
00720   }
00721 
00722   const AlgebraicSymMatrix &covMeasLoc = refTrajPtr->measurementErrors();
00723   const AlgebraicVector &measurementsLoc = refTrajPtr->measurements();
00724   const AlgebraicVector &trajectoryLoc = refTrajPtr->trajectoryPositions();
00725   const TransientTrackingRecHit::ConstRecHitContainer &recHits = refTrajPtr->recHits();
00726   const AlgebraicMatrix &derivatives = refTrajPtr->derivatives();
00727   
00728 // CHK
00729   const int nRow = refTrajPtr->numberOfHitMeas(); 
00730 
00731   for (int iRow = 0; iRow < nRow; ++iRow) {
00732     const double residuum = measurementsLoc[iRow] - trajectoryLoc[iRow];
00733     const double covMeasLocRow = covMeasLoc[iRow][iRow];
00734     const bool is2DhitRow = (!recHits[iRow/2]->detUnit() // FIXME: as in MillePedeAlignmentAlgorithm::is2D()
00735                              || recHits[iRow/2]->detUnit()->type().isTrackerPixel());
00736     //the GeomDetUnit* is zero for composite hits (matched hits in the tracker,...). 
00737     if (TMath::Even(iRow)) { // local x
00738       static const int iMeasLocX = this->GetIndex(myTrajectoryHists1D, "measLocX");
00739       myTrajectoryHists1D[iMeasLocX]->Fill(measurementsLoc[iRow]);
00740       static const int iTrajLocX = this->GetIndex(myTrajectoryHists1D, "trajLocX");
00741       myTrajectoryHists1D[iTrajLocX]->Fill(trajectoryLoc[iRow]);
00742       static const int iResidLocX = this->GetIndex(myTrajectoryHists1D, "residLocX");
00743       myTrajectoryHists1D[iResidLocX]->Fill(residuum);
00744 
00745       static const int iMeasLocXvsHit = this->GetIndex(myTrajectoryHists2D, "measLocXvsHit");
00746       myTrajectoryHists2D[iMeasLocXvsHit]->Fill(iRow/2, measurementsLoc[iRow]);
00747       static const int iTrajLocXvsHit = this->GetIndex(myTrajectoryHists2D, "trajLocXvsHit");
00748       myTrajectoryHists2D[iTrajLocXvsHit]->Fill(iRow/2, trajectoryLoc[iRow]);
00749       static const int iResidLocXvsHit = this->GetIndex(myTrajectoryHists2D, "residLocXvsHit");
00750       myTrajectoryHists2D[iResidLocXvsHit]->Fill(iRow/2, residuum);
00751 
00752       if (covMeasLocRow > 0.) { 
00753         static const int iReduResidLocX = this->GetIndex(myTrajectoryHists1D, "reduResidLocX");
00754         myTrajectoryHists1D[iReduResidLocX]->Fill(residuum / TMath::Sqrt(covMeasLocRow));
00755         static const int iReduResidLocXvsHit = 
00756           this->GetIndex(myTrajectoryHists2D, "reduResidLocXvsHit");
00757         myTrajectoryHists2D[iReduResidLocXvsHit]->Fill(iRow/2, residuum/TMath::Sqrt(covMeasLocRow));
00758       }
00759     } else if (is2DhitRow) { // local y, 2D detectors only
00760 
00761       static const int iMeasLocY = this->GetIndex(myTrajectoryHists1D, "measLocY");
00762       myTrajectoryHists1D[iMeasLocY]->Fill(measurementsLoc[iRow]);
00763       static const int iTrajLocY = this->GetIndex(myTrajectoryHists1D, "trajLocY");
00764       myTrajectoryHists1D[iTrajLocY]->Fill(trajectoryLoc[iRow]);
00765       static const int iResidLocY = this->GetIndex(myTrajectoryHists1D, "residLocY");
00766       myTrajectoryHists1D[iResidLocY]->Fill(residuum);
00767 
00768       static const int iMeasLocYvsHit = this->GetIndex(myTrajectoryHists2D, "measLocYvsHit");
00769       myTrajectoryHists2D[iMeasLocYvsHit]->Fill(iRow/2, measurementsLoc[iRow]);
00770       static const int iTrajLocYvsHit = this->GetIndex(myTrajectoryHists2D, "trajLocYvsHit");
00771       myTrajectoryHists2D[iTrajLocYvsHit]->Fill(iRow/2, trajectoryLoc[iRow]);
00772       static const int iResidLocYvsHit = this->GetIndex(myTrajectoryHists2D, "residLocYvsHit");
00773       myTrajectoryHists2D[iResidLocYvsHit]->Fill(iRow/2, residuum);
00774 
00775       if (covMeasLocRow > 0.) {
00776         static const int iReduResidLocY = this->GetIndex(myTrajectoryHists1D, "reduResidLocY");
00777         myTrajectoryHists1D[iReduResidLocY]->Fill(residuum / TMath::Sqrt(covMeasLocRow));
00778         static const int iReduResidLocYvsHit = this->GetIndex(myTrajectoryHists2D,
00779                                                               "reduResidLocYvsHit");
00780         myTrajectoryHists2D[iReduResidLocYvsHit]->Fill(iRow/2, residuum/TMath::Sqrt(covMeasLocRow));
00781       }
00782     }
00783 
00784     float nHitRow = iRow/2; // '/2' not '/2.'!
00785     if (TMath::Odd(iRow)) nHitRow += 0.5; // y-hit gets 0.5
00786     // correlations
00787     for (int iCol = iRow+1; iCol < nRow; ++iCol) {
00788       double rho = TMath::Sqrt(covMeasLocRow + covMeasLoc[iCol][iCol]);
00789       rho = (0. == rho ? -2 : covMeasLoc[iRow][iCol] / rho);
00790       float nHitCol = iCol/2; //cf. comment nHitRow
00791       if (TMath::Odd(iCol)) nHitCol += 0.5; // dito
00792       //      myProfileCorr->Fill(nHitRow, nHitCol, TMath::Abs(rho));
00793       if (0. == rho) continue;
00794       static const int iProfileCorr = this->GetIndex(myTrajectoryHists2D, "profCorr");
00795       myTrajectoryHists2D[iProfileCorr]->Fill(nHitRow, nHitCol, TMath::Abs(rho));
00796       if (iRow+1 == iCol && TMath::Even(iRow)) { // i.e. if iRow is x and iCol the same hit's y
00797 //      static const int iProfileCorrXy = this->GetIndex(myTrajectoryHists2D,"profCorrOffXy");
00798 //      myTrajectoryHists2D[iProfileCorrOffXy]->Fill(iRow/2, rho);
00799         static const int iHitCorrXy = this->GetIndex(myTrajectoryHists2D, "hitCorrXy");
00800         myTrajectoryHists2D[iHitCorrXy]->Fill(iRow/2, rho);
00801         static const int iHitCorrXyLog = this->GetIndex(myTrajectoryHists2D, "hitCorrXyLog");
00802         myTrajectoryHists2D[iHitCorrXyLog]->Fill(iRow/2, TMath::Abs(rho));
00803         if (is2DhitRow) {
00804           static const int iHitCorrXyValid = this->GetIndex(myTrajectoryHists2D, "hitCorrXyValid");
00805           myTrajectoryHists2D[iHitCorrXyValid]->Fill(iRow/2, rho); // nhitRow??
00806           static const int iHitCorrXyLogValid = this->GetIndex(myTrajectoryHists2D,
00807                                                                "hitCorrXyLogValid");
00808           myTrajectoryHists2D[iHitCorrXyLogValid]->Fill(iRow/2, TMath::Abs(rho)); // nhitRow??
00809         }
00810       } else {
00811         static const int iProfCorrOffXy = this->GetIndex(myTrajectoryHists2D, "profCorrOffXy");
00812         myTrajectoryHists2D[iProfCorrOffXy]->Fill(nHitRow, nHitCol, TMath::Abs(rho));
00813         static const int iHitCorrOffBlock = this->GetIndex(myTrajectoryHists2D, "hitCorrOffBlock");
00814         myTrajectoryHists2D[iHitCorrOffBlock]->Fill(nHitRow, rho);
00815         static const int iHitCorOffBlkLg = this->GetIndex(myTrajectoryHists2D,"hitCorrOffBlockLog");
00816         myTrajectoryHists2D[iHitCorOffBlkLg]->Fill(nHitRow, TMath::Abs(rho));
00817       }
00818     } // end loop on columns of covariance
00819     
00820     // derivatives
00821     for (int iCol = 0; iCol < derivatives.num_col(); ++iCol) {
00822       static const int iProfDerivatives = this->GetIndex(myTrajectoryHists2D, "profDerivatives");
00823       myTrajectoryHists2D[iProfDerivatives]->Fill(nHitRow, iCol, derivatives[iRow][iCol]);
00824       static const int iDerivatives = this->GetIndex(myTrajectoryHists2D, "derivatives");
00825       myTrajectoryHists2D[iDerivatives]->Fill(iCol, derivatives[iRow][iCol]);
00826       static const int iDerivativesLog = this->GetIndex(myTrajectoryHists2D, "derivativesLog");
00827       myTrajectoryHists2D[iDerivativesLog]->Fill(iCol, TMath::Abs(derivatives[iRow][iCol]));
00828       static const int iDerivativesPhi = this->GetIndex(myTrajectoryHists2D, "derivativesVsPhi");
00829       myTrajectoryHists2D[iDerivativesPhi]->Fill(recHits[iRow/2]->det()->position().phi(),
00830                                                  derivatives[iRow][iCol]);
00831     }
00832   } // end loop on rows of covarianvce
00833 
00834 }
00835 
00836 //____________________________________________________________________
00837 void MillePedeMonitor::fillDerivatives(const ConstRecHitPointer &recHit,
00838                                        const float *localDerivs, unsigned int nLocal,
00839                                        const float *globalDerivs, unsigned int nGlobal,
00840                                        const int *labels)
00841 {
00842   const double phi = recHit->det()->position().phi();
00843 
00844   static const int iLocPar = this->GetIndex(myDerivHists2D, "localDerivsPar");
00845   static const int iLocPhi = this->GetIndex(myDerivHists2D, "localDerivsPhi");
00846   static const int iLocParLog = this->GetIndex(myDerivHists2D, "localDerivsParLog");
00847   static const int iLocPhiLog = this->GetIndex(myDerivHists2D, "localDerivsPhiLog");
00848   for (unsigned int i = 0; i < nLocal; ++i) {
00849     myDerivHists2D[iLocPar]->Fill(i, localDerivs[i]);
00850     myDerivHists2D[iLocPhi]->Fill(phi, localDerivs[i]);
00851     if (localDerivs[i]) {
00852       myDerivHists2D[iLocParLog]->Fill(i, TMath::Abs(localDerivs[i]));
00853       myDerivHists2D[iLocPhiLog]->Fill(phi, TMath::Abs(localDerivs[i]));
00854     }
00855   }
00856 
00857   static const int iGlobPar = this->GetIndex(myDerivHists2D, "globalDerivsPar");
00858   static const int iGlobPhi = this->GetIndex(myDerivHists2D, "globalDerivsPhi");
00859   static const int iGlobParLog = this->GetIndex(myDerivHists2D, "globalDerivsParLog");
00860   static const int iGlobPhiLog = this->GetIndex(myDerivHists2D, "globalDerivsPhiLog");
00861 //   static const int iGlobPhiLog2 = this->GetIndex(myDerivHists2D, "globalDerivsPhiLog2");
00862   for (unsigned int i = 0; i < nGlobal; ++i) {
00863     const unsigned int parNum = (labels ? (labels[i]%PedeLabelerBase::theMaxNumParam)-1 : i);
00864     myDerivHists2D[iGlobPar]->Fill(parNum, globalDerivs[i]);
00865     myDerivHists2D[iGlobPhi]->Fill(phi, globalDerivs[i]);
00866     if (globalDerivs[i]) {
00867       myDerivHists2D[iGlobParLog]->Fill(parNum, TMath::Abs(globalDerivs[i]));
00868       myDerivHists2D[iGlobPhiLog]->Fill(phi, TMath::Abs(globalDerivs[i]));
00869 //       myDerivHists2D[iGlobPhiLog2]->Fill(phi, TMath::Abs(globalDerivs[i]));
00870     }
00871   }
00872 
00873 }
00874 
00875 
00876 //____________________________________________________________________
00877 void MillePedeMonitor::fillResiduals(const ConstRecHitPointer &recHit,
00878                                      const TrajectoryStateOnSurface &tsos,
00879                                      unsigned int nHit, float residuum, float sigma, bool isY)
00880 {
00881   // isY == false: x-measurements
00882   // isY == true:  y-measurements
00883   const GlobalPoint detPos(recHit->det()->position());
00884   const double phi = detPos.phi();
00885 //  const double rSigned = (detPos.y() > 0. ? detPos.perp() : -detPos.perp());
00886 
00887   static const int iResPhi = this->GetIndex(myResidHists2D, "residPhi");
00888   myResidHists2D[iResPhi]->Fill(phi, residuum);
00889   static const int iSigPhi = this->GetIndex(myResidHists2D, "sigmaPhi");
00890   myResidHists2D[iSigPhi]->Fill(phi, sigma);
00891   if (sigma) {
00892     static const int iReduResPhi = this->GetIndex(myResidHists2D, "reduResidPhi");
00893     myResidHists2D[iReduResPhi]->Fill(phi, residuum/sigma);
00894   }
00895 
00896 //   if (isY) {
00897 //     static const int iResYXy = this->GetIndex(myResidHists2D, "residYProfXy");
00898 //     myResidHists2D[iResYXy]->Fill(detPos.x(), detPos.y(), TMath::Abs(residuum));
00899 //     static const int iResYZr = this->GetIndex(myResidHists2D, "residYProfZr");
00900 //     myResidHists2D[iResYZr]->Fill(detPos.z(), rSigned, TMath::Abs(residuum));
00901 //     static const int iSigmaYXy = this->GetIndex(myResidHists2D, "sigmaYProfXy");
00902 //     myResidHists2D[iSigmaYXy]->Fill(detPos.x(), detPos.y(), sigma);
00903 //     static const int iSigmaYZr = this->GetIndex(myResidHists2D, "sigmaYProfZr");
00904 //     myResidHists2D[iSigmaYZr]->Fill(detPos.z(), rSigned, sigma);
00905 //     if (sigma) {
00906 //       static const int iReduResYXy = this->GetIndex(myResidHists2D, "reduResidYProfXy");
00907 //       myResidHists2D[iReduResYXy]->Fill(detPos.x(), detPos.y(), TMath::Abs(residuum/sigma));
00908 //       static const int iReduResYZr = this->GetIndex(myResidHists2D, "reduResidYProfZr");
00909 //       myResidHists2D[iReduResYZr]->Fill(detPos.z(), rSigned, TMath::Abs(residuum/sigma));
00910 //     }
00911 //   } else {
00912 //     static const int iResXXy = this->GetIndex(myResidHists2D, "residXProfXy");
00913 //     myResidHists2D[iResXXy]->Fill(detPos.x(), detPos.y(), TMath::Abs(residuum));
00914 //     static const int iResXZr = this->GetIndex(myResidHists2D, "residXProfZr");
00915 //     myResidHists2D[iResXZr]->Fill(detPos.z(), rSigned, TMath::Abs(residuum));
00916 //     static const int iSigmaXXy = this->GetIndex(myResidHists2D, "sigmaXProfXy");
00917 //     myResidHists2D[iSigmaXXy]->Fill(detPos.x(), detPos.y(), sigma);
00918 //     static const int iSigmaXZr = this->GetIndex(myResidHists2D, "sigmaXProfZr");
00919 //     myResidHists2D[iSigmaXZr]->Fill(detPos.z(), rSigned, sigma);
00920 //     if (sigma) {
00921 //       static const int iReduResXXy = this->GetIndex(myResidHists2D, "reduResidXProfXy");
00922 //       myResidHists2D[iReduResXXy]->Fill(detPos.x(), detPos.y(), TMath::Abs(residuum/sigma));
00923 //       static const int iReduResXZr = this->GetIndex(myResidHists2D, "reduResidXProfZr");
00924 //       myResidHists2D[iReduResXZr]->Fill(detPos.z(), rSigned, TMath::Abs(residuum/sigma));
00925 //     }
00926 //   }
00927 
00928   const LocalVector mom(tsos.localDirection()); // mom.z()==0. impossible for TSOS:
00929   const float phiSensToNorm = TMath::ATan(TMath::Abs((isY ? mom.y() : mom.x())/mom.z()));
00930 
00931   std::vector<std::vector<TH1*> > &histArrayVec = (isY ? myResidHistsVec1DY : myResidHistsVec1DX);
00932   std::vector<TH1*> &hitHists =                   (isY ? myResidHitHists1DY : myResidHitHists1DX);
00933 
00934   // call with histArrayVec[0] first: 'static' inside is referring to "subdet-less" names (X/Y irrelevant)
00935   this->fillResidualHists(histArrayVec[0], phiSensToNorm, residuum, sigma);
00936   this->fillResidualHitHists(hitHists, phiSensToNorm, residuum, sigma, nHit);
00937 
00938   const DetId detId(recHit->det()->geographicalId());
00939   if (detId.det() == DetId::Tracker) {
00940     //   const GeomDet::SubDetector subDetNum = recHit->det()->subDetector();
00941     unsigned int subDetNum = detId.subdetId();
00942     if (subDetNum < histArrayVec.size() && subDetNum > 0) {
00943       this->fillResidualHists(histArrayVec[subDetNum], phiSensToNorm, residuum, sigma);
00944     } else {
00945       if (detId!=AlignableBeamSpot::detId())
00946         edm::LogWarning("Alignment") << "@SUB=MillePedeMonitor::fillResiduals"
00947                                      << "Expect subDetNum from 1 to 6, got " << subDetNum;
00948     }
00949   }
00950 }
00951 
00952 //____________________________________________________________________
00953 void MillePedeMonitor::fillResidualHists(const std::vector<TH1*> &hists,
00954                                          float phiSensToNorm, float residuum, float sigma)
00955 {
00956   // Histogram indices are calculated at first call, so make sure the vector of hists at the first
00957   // call has the correct names (i.e. without subdet extension) and all later calls have the
00958   // same order of hists...
00959   
00960   static const int iRes = this->GetIndex(hists, "resid");
00961   hists[iRes]->Fill(residuum);
00962   static const int iSigma = this->GetIndex(hists, "sigma");
00963   hists[iSigma]->Fill(sigma);
00964   static const int iSigmaVsAngle = this->GetIndex(hists, "sigmaVsAngle");
00965   hists[iSigmaVsAngle]->Fill(phiSensToNorm, sigma);
00966   static const int iResidVsAngle = this->GetIndex(hists, "residVsAngle");
00967   hists[iResidVsAngle]->Fill(phiSensToNorm, residuum);
00968   static const int iReduRes = this->GetIndex(hists, "reduResid");
00969   static const int iReduResidVsAngle = this->GetIndex(hists, "reduResidVsAngle");
00970   if (sigma) {
00971     hists[iReduRes]->Fill(residuum/sigma);
00972     hists[iReduResidVsAngle]->Fill(phiSensToNorm, residuum/sigma);
00973   }
00974   static const int iAngle = this->GetIndex(hists, "angle");
00975   hists[iAngle]->Fill(phiSensToNorm);
00976   
00977   if (phiSensToNorm > TMath::DegToRad()*45.) {
00978     static const int iResGt45 = this->GetIndex(hists, "residGt45");
00979     hists[iResGt45]->Fill(residuum);
00980     static const int iSigmaGt45 = this->GetIndex(hists, "sigmaGt45");
00981     hists[iSigmaGt45]->Fill(sigma);
00982     static const int iReduResGt45 = this->GetIndex(hists, "reduResidGt45");
00983     if (sigma) hists[iReduResGt45]->Fill(residuum/sigma);
00984   } else {
00985     static const int iResLt45 = this->GetIndex(hists, "residLt45");
00986     hists[iResLt45]->Fill(residuum);
00987     static const int iSigmaLt45 = this->GetIndex(hists, "sigmaLt45");
00988     hists[iSigmaLt45]->Fill(sigma);
00989     static const int iReduResLt45 = this->GetIndex(hists, "reduResidLt45");
00990     if (sigma) hists[iReduResLt45]->Fill(residuum/sigma);
00991   }
00992 }
00993 
00994 //____________________________________________________________________
00995 void MillePedeMonitor::fillResidualHitHists(const std::vector<TH1*> &hists, float angle,
00996                                             float residuum, float sigma, unsigned int nHit)
00997 {
00998   // Histogram indices are calculated at first call, so make sure the vector of hists at the first
00999   // call has the correct names (i.e. without subdet extension) and all later calls have the
01000   // same order of hists...
01001   
01002   static const unsigned int maxNhit = 29;  // 0...29 foreseen in initialisation...
01003   static int iResHit[maxNhit+1] = {-1};
01004   static int iSigmaHit[maxNhit+1] = {-1};
01005   static int iReduResHit[maxNhit+1] = {-1};
01006   static int iAngleHit[maxNhit+1] = {-1};
01007   if (iResHit[0] == -1) { // first event only...
01008     for (unsigned int i = 0; i <= maxNhit; ++i) {
01009       iResHit[i] = this->GetIndex(hists, Form("resid_%d", i));
01010       iSigmaHit[i] = this->GetIndex(hists, Form("sigma_%d", i));
01011       iReduResHit[i] = this->GetIndex(hists, Form("reduResid_%d", i));
01012       iAngleHit[i] = this->GetIndex(hists, Form("angle_%d", i));
01013     }
01014   }
01015   if (nHit > maxNhit) nHit = maxNhit; // limit of hists
01016 
01017   hists[iResHit[nHit]]->Fill(residuum);
01018   hists[iSigmaHit[nHit]]->Fill(sigma);
01019   if (sigma) {
01020     hists[iReduResHit[nHit]]->Fill(residuum/sigma);
01021   }
01022   hists[iAngleHit[nHit]]->Fill(angle);
01023 }
01024 
01025 //____________________________________________________________________
01026 void MillePedeMonitor::fillFrameToFrame(const AlignableDetOrUnitPtr &aliDet, const Alignable *ali)
01027 {
01028   // get derivative of higher level structure w.r.t. det
01029   FrameToFrameDerivative ftfd;
01030   const AlgebraicMatrix frameDeriv = ftfd.frameToFrameDerivative(aliDet, ali);
01031   //const AlgebraicMatrix frameDeriv = ftfd.frameToFrameDerivativeAtOrgRot(aliDet, ali);
01032 
01033   static const int iF2f = this->GetIndex(myFrame2FrameHists2D, "frame2frame");
01034   static const int iF2fAbs = this->GetIndex(myFrame2FrameHists2D, "frame2frameAbs");
01035   static int iF2fIjPhi[6][6], iF2fIjPhiLog[6][6], iF2fIjR[6][6], iF2fIjRLog[6][6];
01036   static bool first = true;
01037   if (first) {
01038     for (unsigned int i = 0; i < 6; ++i) {
01039       for (unsigned int j = 0; j < 6; ++j) {
01040         iF2fIjPhi[i][j] = this->GetIndex(myFrame2FrameHists2D, Form("frame2framePhi%d%d", i, j));
01041         iF2fIjPhiLog[i][j]=this->GetIndex(myFrame2FrameHists2D, Form("frame2framePhiLog%d%d",i,j));
01042         iF2fIjR[i][j] = this->GetIndex(myFrame2FrameHists2D, Form("frame2frameR%d%d", i, j));
01043         iF2fIjRLog[i][j]=this->GetIndex(myFrame2FrameHists2D, Form("frame2frameRLog%d%d",i,j));
01044       }
01045     }
01046     first = false;
01047   }
01048   
01049   const double phi = aliDet->globalPosition().phi(); // after misalignment...
01050   const double r = aliDet->globalPosition().perp(); // after misalignment...
01051   for (unsigned int i = 0; i < 6; ++i) {
01052     for (unsigned int j = 0; j < 6; ++j) {
01053       myFrame2FrameHists2D[iF2f]->Fill(i, j, frameDeriv[i][j]);
01054       myFrame2FrameHists2D[iF2fIjPhi[i][j]]->Fill(phi, frameDeriv[i][j]);
01055       myFrame2FrameHists2D[iF2fIjR[i][j]]->Fill(r, frameDeriv[i][j]);
01056       if (frameDeriv[i][j]) {
01057         myFrame2FrameHists2D[iF2fAbs]->Fill(i, j, TMath::Abs(frameDeriv[i][j]));
01058         myFrame2FrameHists2D[iF2fIjPhiLog[i][j]]->Fill(phi, TMath::Abs(frameDeriv[i][j]));
01059         myFrame2FrameHists2D[iF2fIjRLog[i][j]]->Fill(r, TMath::Abs(frameDeriv[i][j]));
01060       }
01061     }
01062   }
01063 }
01064 
01065 //____________________________________________________________________
01066 void MillePedeMonitor::fillCorrelations2D(float corr, const ConstRecHitPointer &recHit)
01067 {
01068   const DetId detId(recHit->det()->geographicalId());
01069   if (detId.det() != DetId::Tracker) return;
01070 
01071   if ((detId.subdetId() < 1 || detId.subdetId() > 6) &&
01072       detId!=AlignableBeamSpot::detId()) {
01073     edm::LogWarning("Alignment") << "@SUB=MillePedeMonitor::fillCorrelations2D"
01074                                  << "Expect subdetId from 1 to 6, got " << detId.subdetId();
01075     return;
01076   }
01077 
01078   myCorrHists[detId.subdetId()-1]->Fill(corr);
01079 }
01080 
01081 //____________________________________________________________________
01082 void MillePedeMonitor::fillPxbSurveyHistsChi2(const float &chi2)
01083 {
01084   static const int iPxbSurvChi2_lo = this->GetIndex(myPxbSurveyHists,"PxbSurvChi2_lo");
01085   myPxbSurveyHists[iPxbSurvChi2_lo]->Fill(chi2);
01086   static const int iPxbSurvChi2_md = this->GetIndex(myPxbSurveyHists,"PxbSurvChi2_md");
01087   myPxbSurveyHists[iPxbSurvChi2_md]->Fill(chi2);
01088   static const int iPxbSurvChi2_hi = this->GetIndex(myPxbSurveyHists,"PxbSurvChi2_hi");
01089   myPxbSurveyHists[iPxbSurvChi2_hi]->Fill(chi2);
01090   static const int iPxbSurvChi2prob = this->GetIndex(myPxbSurveyHists,"PxbSurvChi2prob");
01091   myPxbSurveyHists[iPxbSurvChi2prob]->Fill(TMath::Prob(chi2,4));
01092 }
01093 
01094 //____________________________________________________________________
01095 void MillePedeMonitor::fillPxbSurveyHistsLocalPars(const float &a0, const float &a1, const float &S, const float &phi)
01096 {
01097   static const int iPxbSurv_a0 = this->GetIndex(myPxbSurveyHists,"PxbSurv_a0");
01098   myPxbSurveyHists[iPxbSurv_a0]->Fill(a0);
01099   static const int iPxbSurv_a0_abs = this->GetIndex(myPxbSurveyHists,"PxbSurv_a0_abs");
01100   myPxbSurveyHists[iPxbSurv_a0_abs]->Fill(fabs(a0));
01101   static const int iPxbSurv_a1 = this->GetIndex(myPxbSurveyHists,"PxbSurv_a1");
01102   myPxbSurveyHists[iPxbSurv_a1]->Fill(a1);
01103   static const int iPxbSurv_a1_abs = this->GetIndex(myPxbSurveyHists,"PxbSurv_a1_abs");
01104   myPxbSurveyHists[iPxbSurv_a1_abs]->Fill(fabs(a1));
01105   static const int iPxbSurv_scale = this->GetIndex(myPxbSurveyHists,"PxbSurv_scale");
01106   myPxbSurveyHists[iPxbSurv_scale]->Fill(S);
01107   static const int iPxbSurv_phi = this->GetIndex(myPxbSurveyHists,"PxbSurv_phi");
01108   myPxbSurveyHists[iPxbSurv_phi]->Fill(phi);
01109 }
01110 
01111
/data/refman/pasoursint/CMSSW_5_3_0/src/Alignment/MillePedeAlignmentAlgorithm/src/MillePedeMonitor.cc
