Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:58:34 

0001 /****************************************************************************
0002 * Authors: 
0003 *  Jan Kašpar (jan.kaspar@gmail.com) 
0004 ****************************************************************************/
0005 
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008 
0009 #include "CalibPPS/AlignmentRelative/interface/JanAlignmentAlgorithm.h"
0010 #include "CalibPPS/AlignmentRelative/interface/AlignmentTask.h"
0011 #include "CalibPPS/AlignmentRelative/interface/Utilities.h"
0012 
0013 #include "TMatrixDSymEigen.h"
0014 #include "TDecompSVD.h"
0015 #include "TFile.h"
0016 #include "TCanvas.h"
0017 #include "TH2D.h"
0018 
0019 #include <cmath>
0020 
0021 //#define DEBUG 1
0022 
0023 using namespace std;
0024 using namespace edm;
0025 
0026 //----------------------------------------------------------------------------------------------------
0027 
0028 JanAlignmentAlgorithm::JanAlignmentAlgorithm(const ParameterSet &ps, AlignmentTask *_t)
0029     : AlignmentAlgorithm(ps, _t), Sc(nullptr), Mc(nullptr) {
0030   const ParameterSet &lps = ps.getParameterSet("JanAlignmentAlgorithm");
0031   weakLimit = lps.getParameter<double>("weakLimit");
0032   stopOnSingularModes = lps.getParameter<bool>("stopOnSingularModes");
0033   buildDiagnosticPlots = lps.getParameter<bool>("buildDiagnosticPlots");
0034 }
0035 
0036 //----------------------------------------------------------------------------------------------------
0037 
0038 JanAlignmentAlgorithm::~JanAlignmentAlgorithm() {}
0039 
0040 //----------------------------------------------------------------------------------------------------
0041 
0042 void JanAlignmentAlgorithm::begin(const CTPPSGeometry *geometryReal, const CTPPSGeometry *geometryMisaligned) {
0043   // initialize M and S components
0044   Mc = new TVectorD[task->quantityClasses.size()];
0045   Sc = new TMatrixD *[task->quantityClasses.size()];
0046   for (unsigned int i = 0; i < task->quantityClasses.size(); i++) {
0047     unsigned int rows = task->quantitiesOfClass(task->quantityClasses[i]);
0048 
0049     Mc[i].ResizeTo(rows);
0050     Mc[i].Zero();
0051 
0052     Sc[i] = new TMatrixD[task->quantityClasses.size()];
0053     for (unsigned int j = 0; j < task->quantityClasses.size(); j++) {
0054       unsigned int cols = task->quantitiesOfClass(task->quantityClasses[j]);
0055       Sc[i][j].ResizeTo(rows, cols);
0056       Sc[i][j].Zero();
0057     }
0058   }
0059 
0060   // prepare statistics plots
0061   if (buildDiagnosticPlots) {
0062     for (const auto &it : task->geometry.getSensorMap()) {
0063       unsigned int id = it.first;
0064       char buf[50];
0065       DetStat s;
0066 
0067       sprintf(buf, "%u: m distribution", id);
0068       s.m_dist = new TH1D(buf, ";u or v   (mm)", 100, -25., 25.);
0069 
0070       sprintf(buf, "%u: R distribution", id);
0071       s.R_dist = new TH1D(buf, ";R   (mm)", 500, -0.5, 0.5);
0072 
0073       for (unsigned int c = 0; c < task->quantityClasses.size(); c++) {
0074         sprintf(buf, "%u: coef, %s", id, task->quantityClassTag(task->quantityClasses[c]).c_str());
0075         s.coefHist.push_back(new TH1D(buf, ";coefficient", 100, -2., +2.));
0076 
0077         sprintf(buf, "%u: R vs. coef, %s", id, task->quantityClassTag(task->quantityClasses[c]).c_str());
0078         TGraph *g = new TGraph();
0079         g->SetName(buf);
0080         g->SetTitle(";coefficient;residual   (mm)");
0081         s.resVsCoef.push_back(g);
0082       }
0083 
0084       statistics[id] = s;
0085     }
0086   }
0087 
0088   events = 0;
0089 }
0090 
0091 //----------------------------------------------------------------------------------------------------
0092 
0093 void JanAlignmentAlgorithm::feed(const HitCollection &selection, const LocalTrackFit &trackFit) {
0094   if (verbosity > 9)
0095     printf("\n>> JanAlignmentAlgorithm::Feed\n");
0096 
0097   events++;
0098 
0099   // prepare fit - make z0 compatible
0100   double hax = trackFit.ax;
0101   double hay = trackFit.ay;
0102   double hbx = trackFit.bx + trackFit.ax * (task->geometry.z0 - trackFit.z0);
0103   double hby = trackFit.by + trackFit.ay * (task->geometry.z0 - trackFit.z0);
0104 
0105   // prepare Gamma matrices (full of zeros)
0106   TMatrixD *Ga = new TMatrixD[task->quantityClasses.size()];
0107   for (unsigned int i = 0; i < task->quantityClasses.size(); i++) {
0108     Ga[i].ResizeTo(selection.size(), Mc[i].GetNrows());
0109     Ga[i].Zero();
0110   }
0111 
0112   TMatrixD A(selection.size(), 4);
0113   TMatrixD Vi(selection.size(), selection.size());
0114   TVectorD m(selection.size());
0115 
0116   set<unsigned int> rpSet;
0117   if (buildDiagnosticPlots) {
0118     for (const auto &hit : selection) {
0119       CTPPSDetId detId(hit.id);
0120       const unsigned int rpDecId = 100 * detId.arm() + 10 * detId.station() + detId.rp();
0121       rpSet.insert(rpDecId);
0122     }
0123   }
0124 
0125   // fill fit matrix and Gamma matrices
0126   unsigned int j = 0;
0127 
0128   for (HitCollection::const_iterator hit = selection.begin(); hit != selection.end(); ++hit, ++j) {
0129     unsigned int id = hit->id;
0130 
0131     const DetGeometry &d = task->geometry.get(id);
0132     const auto &dirData = d.getDirectionData(hit->dirIdx);
0133 
0134     A(j, 0) = hit->z * dirData.dx;
0135     A(j, 1) = dirData.dx;
0136     A(j, 2) = hit->z * dirData.dy;
0137     A(j, 3) = dirData.dy;
0138 
0139     m(j) = hit->position + dirData.s - (hit->z - d.z) * dirData.dz;  // in mm
0140 
0141     Vi(j, j) = 1. / hit->sigma / hit->sigma;
0142 
0143     double C = dirData.dx, S = dirData.dy;
0144 
0145     double hx = hax * hit->z + hbx;  // in mm
0146     double hy = hay * hit->z + hby;
0147     double R = m(j) - (hx * C + hy * S);  // (standard) residual
0148 
0149     if (buildDiagnosticPlots) {
0150       statistics[id].m_dist->Fill(m(j));
0151       statistics[id].R_dist->Fill(R);
0152     }
0153 
0154     for (unsigned int i = 0; i < task->quantityClasses.size(); i++) {
0155       // check compatibility
0156       signed int matrixIndex = task->getMeasurementIndex(task->quantityClasses[i], hit->id, hit->dirIdx);
0157       if (matrixIndex < 0)
0158         continue;
0159 
0160       matrixIndex = task->getQuantityIndex(task->quantityClasses[i], hit->id);
0161 
0162       switch (task->quantityClasses[i]) {
0163         case AlignmentTask::qcShR1:
0164           Ga[i][j][matrixIndex] = -1.;
0165           break;
0166 
0167         case AlignmentTask::qcShR2:
0168           Ga[i][j][matrixIndex] = -1.;
0169           break;
0170 
0171         case AlignmentTask::qcShZ:
0172           Ga[i][j][matrixIndex] = hax * C + hay * S;
0173           break;
0174 
0175         case AlignmentTask::qcRotZ:
0176           Ga[i][j][matrixIndex] = (hax * hit->z + hbx - d.sx) * (-S) + (hay * hit->z + hby - d.sy) * C;
0177           break;
0178       }
0179 
0180       if (buildDiagnosticPlots) {
0181         double c = Ga[i][j][matrixIndex];
0182         DetStat &s = statistics[id];
0183         s.coefHist[i]->Fill(c);
0184         s.resVsCoef[i]->SetPoint(s.resVsCoef[i]->GetN(), c, R);
0185 
0186         if (task->quantityClasses[i] == AlignmentTask::qcRotZ) {
0187           map<set<unsigned int>, ScatterPlot>::iterator hit = s.resVsCoefRot_perRPSet.find(rpSet);
0188           if (hit == s.resVsCoefRot_perRPSet.end()) {
0189             ScatterPlot sp;
0190             sp.g = new TGraph();
0191             sp.h = new TH2D("", "", 40, -20., +20., 60, -0.15, +0.15);
0192             hit = s.resVsCoefRot_perRPSet.insert(pair<set<unsigned int>, ScatterPlot>(rpSet, sp)).first;
0193           }
0194           hit->second.g->SetPoint(hit->second.g->GetN(), c, R);
0195           hit->second.h->Fill(c, R);
0196         }
0197       }
0198     }
0199   }
0200 
0201   // sigma matrix
0202   TMatrixD AT(TMatrixD::kTransposed, A);
0203   TMatrixD ATViA(4, 4);
0204   ATViA = AT * Vi * A;
0205   TMatrixD ATViAI(ATViA);
0206   ATViAI = ATViA.Invert();
0207   TMatrixD sigma(Vi);
0208   sigma -= Vi * A * ATViAI * AT * Vi;
0209 
0210   // traspose Gamma matrices
0211   TMatrixD *GaT = new TMatrixD[task->quantityClasses.size()];
0212   for (unsigned int i = 0; i < task->quantityClasses.size(); i++) {
0213     GaT[i].ResizeTo(Mc[i].GetNrows(), selection.size());
0214     GaT[i].Transpose(Ga[i]);
0215   }
0216 
0217   // normalized residuals
0218   TVectorD r(selection.size());
0219   r = sigma * m;
0220 
0221   // increment M
0222   for (unsigned int i = 0; i < task->quantityClasses.size(); i++) {
0223     if (Mc[i].GetNrows() < 1)
0224       continue;
0225 
0226     Mc[i] += GaT[i] * r;
0227   }
0228 
0229   // increment S
0230   for (unsigned int i = 0; i < task->quantityClasses.size(); i++) {
0231     for (unsigned int j = 0; j < task->quantityClasses.size(); j++) {
0232       if (Sc[i][j].GetNrows() < 1 || Sc[i][j].GetNcols() < 1)
0233         continue;
0234 
0235       Sc[i][j] += GaT[i] * sigma * Ga[j];
0236     }
0237   }
0238 
0239 #ifdef DEBUG
0240   printf("* checking normalized residuals, selection.size = %u\n", selection.size());
0241   r.Print();
0242 
0243   for (unsigned int i = 0; i < task->quantityClasses.size(); i++) {
0244     printf("- class %u\n", i);
0245     TVectorD t(Mc[i].GetNrows());
0246     for (int j = 0; j < t.GetNrows(); j++)
0247       t[j] = 1.;
0248     t.Print();
0249 
0250     Ga[i].Print();
0251 
0252     TVectorD tt(selection.size());
0253     tt = sigma * Ga[i] * t;
0254 
0255     double ttn = sqrt(tt.Norm2Sqr());
0256     printf("|tt| = %E\n", ttn);
0257     if (ttn > 1E-8)
0258       tt.Print();
0259   }
0260 #endif
0261 
0262   delete[] Ga;
0263   delete[] GaT;
0264 }
0265 
0266 //----------------------------------------------------------------------------------------------------
0267 
0268 void JanAlignmentAlgorithm::analyze() {
0269   if (verbosity > 2)
0270     printf("\n>> JanAlignmentAlgorithm::Analyze\n");
0271 
0272   // calculate full dimension
0273   unsigned int dim = 0;
0274   for (unsigned int i = 0; i < task->quantityClasses.size(); i++)
0275     dim += Mc[i].GetNrows();
0276 
0277   if (verbosity > 2) {
0278     printf("\tsensors: %u\n", task->geometry.getNumberOfDetectors());
0279     printf("\tfull dimension: %u\n", dim);
0280     printf("\tquantity classes: %lu\n", task->quantityClasses.size());
0281   }
0282 
0283   // build full M
0284   M.ResizeTo(dim);
0285   unsigned int offset = 0;
0286   for (unsigned int i = 0; i < task->quantityClasses.size(); i++) {
0287     M.SetSub(offset, Mc[i]);
0288     offset += Mc[i].GetNrows();
0289   }
0290 
0291   // build full S
0292   S.ResizeTo(dim, dim);
0293   unsigned int r_offset = 0, c_offset = 0;
0294   for (unsigned int i = 0; i < task->quantityClasses.size(); i++) {
0295     c_offset = 0;
0296     unsigned int r_size = 0, c_size = 0;
0297     for (unsigned int j = 0; j < task->quantityClasses.size(); j++) {
0298       r_size = Sc[i][j].GetNrows();
0299       c_size = Sc[i][j].GetNcols();
0300 
0301       if (r_size < 1 || c_size < 1)
0302         continue;
0303 
0304       TMatrixDSub(S, r_offset, r_offset + r_size - 1, c_offset, c_offset + c_size - 1) = Sc[i][j];
0305       c_offset += c_size;
0306     }
0307     r_offset += r_size;
0308   }
0309 
0310   // analyze symmetricity
0311   if (verbosity >= 3) {
0312     double maxDiff = 0., maxElem = 0.;
0313     for (unsigned int i = 0; i < dim; i++) {
0314       for (unsigned int j = 0; j < dim; j++) {
0315         double diff = S[i][j] - S[j][i];
0316         if (fabs(diff) > maxDiff)
0317           maxDiff = diff;
0318         if (S[i][j] > maxElem)
0319           maxElem = S[i][j];
0320       }
0321     }
0322 
0323     printf("\n* S matrix:\n\tdimension = %i\n\tmaximum asymmetry: %E\t(ratio to maximum element %E)\n",
0324            dim,
0325            maxDiff,
0326            maxDiff / maxElem);
0327   }
0328 
0329   // make a symmetric copy
0330   TMatrixDSym S_sym(dim);
0331   for (unsigned int j = 0; j < dim; j++) {
0332     for (unsigned int i = 0; i < dim; i++) {
0333       S_sym[i][j] = S[i][j];
0334     }
0335   }
0336 
0337   // eigen analysis of S
0338   TMatrixDSymEigen S_eig(S_sym);
0339   const TVectorD &S_eigVal_temp = S_eig.GetEigenValues();
0340   S_eigVal.ResizeTo(S_eigVal_temp.GetNrows());
0341   S_eigVal = S_eigVal_temp;
0342   const TMatrixD &S_eigVec_temp = S_eig.GetEigenVectors();
0343   S_eigVec.ResizeTo(S_eigVec_temp);
0344   S_eigVec = S_eigVec_temp;
0345 
0346   // identify singular modes
0347   for (int i = 0; i < S_eigVal.GetNrows(); i++) {
0348     double nev = S_eigVal[i] / events;
0349     if (fabs(nev) < singularLimit) {
0350       SingularMode sM{S_eigVal[i], TMatrixDColumn(S_eigVec, i), i};
0351       singularModes.push_back(sM);
0352     }
0353   }
0354 
0355 #if DEBUG
0356   // print singular vectors
0357   if (singularModes.size() > 0) {
0358     printf("\n* S singular modes\n   | ");
0359     for (unsigned int i = 0; i < singularModes.size(); i++)
0360       printf("%+10.3E   ", singularModes[i].val);
0361     printf("\n-- | ");
0362 
0363     for (unsigned int i = 0; i < singularModes.size(); i++)
0364       printf("----------   ");
0365     printf("\n");
0366 
0367     for (unsigned int j = 0; j < dim; j++) {
0368       printf("%2u | ", j);
0369       for (unsigned int i = 0; i < singularModes.size(); i++) {
0370         printf("%+10.3E   ", singularModes[i].vec[j]);
0371       }
0372       printf("\n");
0373     }
0374   } else
0375     printf("\n* S has no singular modes\n");
0376 #endif
0377 }
0378 
0379 //----------------------------------------------------------------------------------------------------
0380 
0381 unsigned int JanAlignmentAlgorithm::solve(const std::vector<AlignmentConstraint> &constraints,
0382                                           map<unsigned int, AlignmentResult> &results,
0383                                           TDirectory *dir) {
0384   if (verbosity)
0385     printf(">> JanAlignmentAlgorithm::Solve\n");
0386 
0387   results.clear();
0388 
0389   // build C matrix
0390   unsigned int dim = S.GetNrows();
0391   TMatrixD C(dim, constraints.size());
0392   TMatrixD C2(dim, constraints.size());
0393   for (unsigned int i = 0; i < constraints.size(); i++) {
0394     unsigned int offset = 0;
0395     for (auto &quantityClass : task->quantityClasses) {
0396       const TVectorD &cv = constraints[i].coef.find(quantityClass)->second;
0397       for (int k = 0; k < cv.GetNrows(); k++) {
0398         C[offset][i] = events * cv[k];
0399         C2[offset][i] = events * cv[k] * 1E3;
0400         offset++;
0401       }
0402     }
0403   }
0404 
0405 #ifdef DEBUG
0406   printf("\n* constraint matrix\n");
0407   Print(C);
0408 #endif
0409 
0410   // build E matrix (singular vectors of S as its columns)
0411   TMatrixD E(S.GetNrows(), singularModes.size());
0412   for (unsigned int i = 0; i < singularModes.size(); i++)
0413     for (int j = 0; j < S.GetNrows(); j++)
0414       E(j, i) = singularModes[i].vec[j];
0415 
0416   // build CS matrix
0417   TMatrixDSym CS(dim + constraints.size());
0418   TMatrixDSym CS2(dim + constraints.size());
0419   CS.Zero();
0420   CS2.Zero();
0421   for (unsigned int j = 0; j < dim; j++) {
0422     for (unsigned int i = 0; i < dim; i++) {
0423       CS[i][j] = S[i][j];
0424       CS2[i][j] = S[i][j];
0425     }
0426   }
0427 
0428   for (unsigned int i = 0; i < constraints.size(); i++) {
0429     for (unsigned int j = 0; j < dim; j++) {
0430       CS[j][dim + i] = CS[dim + i][j] = C(j, i);
0431       CS2[j][dim + i] = CS2[dim + i][j] = C2(j, i);
0432     }
0433   }
0434 
0435   // eigen analysis of CS matrix
0436   TMatrixDSymEigen CS_eig(CS);
0437   TVectorD CS_eigVal = CS_eig.GetEigenValues();
0438   TMatrixD CS_eigVec = CS_eig.GetEigenVectors();
0439 
0440   // check regularity of CS matrix
0441   if (verbosity >= 2) {
0442     printf("\n* eigen values of CS and S matrices (events = %u)\n", events);
0443     printf("   #          CS    norm. CS               S     norm. S\n");
0444   }
0445 
0446   unsigned int singularModeCount = 0;
0447   vector<unsigned int> weakModeIdx;
0448   for (int i = 0; i < CS_eigVal.GetNrows(); i++) {
0449     const double CS_nev = CS_eigVal[i] / events;
0450 
0451     if (fabs(CS_nev) < singularLimit)
0452       singularModeCount++;
0453 
0454     if (verbosity >= 2) {
0455       printf("%4i%+12.2E%+12.2E", i, CS_eigVal[i], CS_nev);
0456       if (fabs(CS_nev) < singularLimit) {
0457         singularModeCount++;
0458         printf(" (S)");
0459       } else {
0460         if (fabs(CS_nev) < weakLimit) {
0461           weakModeIdx.push_back(i);
0462           printf(" (W)");
0463         } else {
0464           printf("    ");
0465         }
0466       }
0467 
0468       if (i < S_eigVal.GetNrows()) {
0469         double S_nev = S_eigVal[i] / events;
0470         printf("%+12.2E%+12.2E", S_eigVal[i], S_nev);
0471         if (fabs(S_nev) < singularLimit)
0472           printf(" (S)");
0473         else if (fabs(S_nev) < weakLimit)
0474           printf(" (W)");
0475       }
0476 
0477       printf("\n");
0478     }
0479   }
0480 
0481   if (verbosity >= 2) {
0482     // print weak vectors
0483     if (!weakModeIdx.empty()) {
0484       unsigned int columns = 10;
0485       unsigned int first = 0;
0486 
0487       while (first < weakModeIdx.size()) {
0488         unsigned int last = first + columns;
0489         if (last >= weakModeIdx.size())
0490           last = weakModeIdx.size();
0491 
0492         printf("\n* CS weak modes\n    | ");
0493         for (unsigned int i = first; i < last; i++)
0494           printf("%+10.3E   ", CS_eigVal[weakModeIdx[i]]);
0495         printf("\n--- | ");
0496 
0497         for (unsigned int i = first; i < last; i++)
0498           printf("----------   ");
0499         printf("\n");
0500 
0501         // determine maximum elements
0502         vector<double> maxs;
0503         for (unsigned int i = first; i < last; i++) {
0504           double max = 0;
0505           for (unsigned int j = 0; j < dim + constraints.size(); j++) {
0506             double v = fabs(CS_eigVec(weakModeIdx[i], j));
0507             if (v > max)
0508               max = v;
0509           }
0510           maxs.push_back(max);
0511         }
0512 
0513         for (unsigned int j = 0; j < dim + constraints.size(); j++) {
0514           printf("%3u | ", j);
0515           for (unsigned int i = first; i < last; i++) {
0516             double v = CS_eigVec(weakModeIdx[i], j);
0517             if (fabs(v) / maxs[i - first] > 1E-3)
0518               printf("%+10.3E   ", v);
0519             else
0520               printf("         .   ");
0521           }
0522           printf("\n");
0523         }
0524 
0525         first = last;
0526       }
0527     } else
0528       printf("\n* CS has no weak modes\n");
0529   }
0530 
0531   // check the regularity of C^T E
0532   if (verbosity >= 2) {
0533     if (E.GetNcols() == C.GetNcols()) {
0534       TMatrixD CTE(C, TMatrixD::kTransposeMult, E);
0535       print(CTE, "* CTE matrix:");
0536       const double &det = CTE.Determinant();
0537       printf(
0538           "\n* det(CTE) = %E, max(CTE) = %E, det(CTE)/max(CTE) = %E\n\tmax(C) = %E, max(E) = %E, "
0539           "det(CTE)/max(C)/max(E) = %E\n",
0540           det,
0541           CTE.Max(),
0542           det / CTE.Max(),
0543           C.Max(),
0544           E.Max(),
0545           det / C.Max() / E.Max());
0546     } else {
0547       printf(">> JanAlignmentAlgorithm::Solve > WARNING: C matrix has %u, while E matrix %u columns.\n",
0548              C.GetNcols(),
0549              E.GetNcols());
0550     }
0551   }
0552 
0553   // stop if CS is singular
0554   if (singularModeCount > 0 && stopOnSingularModes) {
0555     LogError("PPS") << "\n>> JanAlignmentAlgorithm::Solve > ERROR: There are " << singularModeCount
0556                     << " singular modes in CS matrix. Stopping.";
0557     return 1;
0558   }
0559 
0560   // build MV vector
0561   TVectorD MV(dim + constraints.size());
0562   for (unsigned int i = 0; i < dim; i++)
0563     MV[i] = M[i];
0564   for (unsigned int i = 0; i < constraints.size(); i++)
0565     MV[dim + i] = events * constraints[i].val;
0566 
0567   // perform inversion and solution
0568   TMatrixD CSI(TMatrixD::kInverted, CS);
0569   TMatrixD CS2I(TMatrixD::kInverted, CS2);
0570   TVectorD AL(MV);
0571   AL = CSI * MV;
0572 
0573   // evaluate error matrix
0574   TMatrixD S0(S);  // new parts full of zeros
0575   S0.ResizeTo(dim + constraints.size(), dim + constraints.size());
0576   TMatrixD EM(CS);
0577   EM = CSI * S0 * CSI;
0578 
0579   TMatrixD EM2(CS2);
0580   EM2 = CS2I * S0 * CS2I;
0581 
0582   TMatrixD EMdiff(EM2 - EM);
0583 
0584   if (verbosity >= 3) {
0585     double max1 = -1., max2 = -1., maxDiff = -1.;
0586     for (int i = 0; i < EMdiff.GetNrows(); i++) {
0587       for (int j = 0; j < EMdiff.GetNcols(); j++) {
0588         if (maxDiff < fabs(EMdiff(i, j)))
0589           maxDiff = fabs(EMdiff(i, j));
0590 
0591         if (max1 < fabs(EM(i, j)))
0592           max1 = fabs(EM(i, j));
0593 
0594         if (max2 < fabs(EM2(i, j)))
0595           max2 = fabs(EM2(i, j));
0596       }
0597     }
0598 
0599     printf("EM max = %E, EM2 max = %E, EM2 - EM max = %E\n", max1, max2, maxDiff);
0600   }
0601 
0602   // tests
0603   TMatrixD &U = CS_eigVec;
0604   TMatrixD UT(TMatrixD::kTransposed, U);
0605   //TMatrixD CSEi(CS);
0606   //CSEi = UT * CS * U;
0607   //Print(CSEi, "CSEi");
0608 
0609   TMatrixD EMEi(EM);
0610   EMEi = UT * EM * U;
0611   //Print(EMEi, "*EMEi");
0612 
0613   if (verbosity >= 3) {
0614     double max = -1.;
0615     for (int i = 0; i < EMEi.GetNrows(); i++) {
0616       for (int j = 0; j < EMEi.GetNcols(); j++) {
0617         if (max < EMEi(i, j))
0618           max = EMEi(i, j);
0619       }
0620     }
0621 
0622     printf("max = %E\n", max);
0623   }
0624 
0625   // print lambda values
0626   if (verbosity >= 3) {
0627     printf("\n* Lambda (from the contribution of singular modes to MV)\n");
0628     for (unsigned int i = 0; i < constraints.size(); i++) {
0629       printf("\t%u (%25s)\t%+10.1E +- %10.1E\n",
0630              i,
0631              constraints[i].name.c_str(),
0632              AL[dim + i] * 1E3,
0633              sqrt(EM[i + dim][i + dim]) * 1E3);
0634     }
0635   }
0636 
0637   // fill results
0638   unsigned int offset = 0;
0639   vector<unsigned int> offsets;
0640   for (unsigned int i = 0; i < task->quantityClasses.size(); i++) {
0641     offsets.push_back(offset);
0642     offset += Mc[i].GetNrows();
0643   }
0644 
0645   for (const auto &dit : task->geometry.getSensorMap()) {
0646     AlignmentResult r;
0647 
0648     for (unsigned int i = 0; i < task->quantityClasses.size(); i++) {
0649       signed idx = task->getQuantityIndex(task->quantityClasses[i], dit.first);
0650       if (idx < 0)
0651         continue;
0652 
0653       unsigned int fi = offsets[i] + idx;
0654       double v = AL[fi];
0655       double e = sqrt(EM[fi][fi]);
0656       switch (task->quantityClasses[i]) {
0657         case AlignmentTask::qcShR1:
0658           r.setShR1(v);
0659           r.setShR1Unc(e);
0660           break;
0661         case AlignmentTask::qcShR2:
0662           r.setShR2(v);
0663           r.setShR2Unc(e);
0664           break;
0665         case AlignmentTask::qcShZ:
0666           r.setShZ(v);
0667           r.setShZUnc(e);
0668           break;
0669         case AlignmentTask::qcRotZ:
0670           r.setRotZ(v);
0671           r.setRotZUnc(e);
0672           break;
0673       }
0674     }
0675 
0676     results[dit.first] = r;
0677   }
0678 
0679   // save matrices, eigen data, ...
0680   if (dir) {
0681     dir->cd();
0682 
0683     S.Write("S");
0684     S_eigVal.Write("S_eigen_values");
0685     S_eigVec.Write("S_eigen_vectors");
0686 
0687     E.Write("E");
0688     C.Write("C");
0689 
0690     CS.Write("CS");
0691     CS_eigVal.Write("CS_eigen_values");
0692     CS_eigVec.Write("CS_eigen_vectors");
0693 
0694     MV.Write("MV");
0695     AL.Write("AL");
0696 
0697     S0.Write("S0");
0698     EM.Write("EM");
0699   }
0700 
0701   // success
0702   return 0;
0703 }
0704 
0705 //----------------------------------------------------------------------------------------------------
0706 
0707 void JanAlignmentAlgorithm::end() {
0708   delete[] Mc;
0709 
0710   for (unsigned int i = 0; i < task->quantityClasses.size(); i++) {
0711     delete[] Sc[i];
0712   }
0713   delete[] Sc;
0714 }
0715 
0716 //----------------------------------------------------------------------------------------------------
0717 
0718 void JanAlignmentAlgorithm::saveDiagnostics(TDirectory *dir) {
0719   if (!buildDiagnosticPlots)
0720     return;
0721 
0722   for (map<unsigned int, DetStat>::iterator it = statistics.begin(); it != statistics.end(); ++it) {
0723     char buf[50];
0724     sprintf(buf, "%u", it->first);
0725     gDirectory = dir->mkdir(buf);
0726 
0727     it->second.m_dist->Write();
0728     it->second.R_dist->Write();
0729 
0730     for (unsigned int c = 0; c < task->quantityClasses.size(); c++) {
0731       it->second.coefHist[c]->Write();
0732       it->second.resVsCoef[c]->Write();
0733     }
0734 
0735     gDirectory = gDirectory->mkdir("R vs. rot. coef, per RP set");
0736     TCanvas *c = new TCanvas;
0737     c->SetName("R vs. rot. coef, overlapped");
0738     TH2D *frame = new TH2D("frame", "frame", 100, -20., +20., 100, -0.15, +0.15);
0739     frame->Draw();
0740     unsigned int idx = 0;
0741     for (map<set<unsigned int>, ScatterPlot>::iterator iit = it->second.resVsCoefRot_perRPSet.begin();
0742          iit != it->second.resVsCoefRot_perRPSet.end();
0743          ++iit, ++idx) {
0744       string label;
0745       bool first = true;
0746       for (set<unsigned int>::iterator sit = iit->first.begin(); sit != iit->first.end(); ++sit) {
0747         char buf[50];
0748         sprintf(buf, "%u", *sit);
0749 
0750         if (first) {
0751           label = buf;
0752           first = false;
0753         } else {
0754           label = label + ", " + buf;
0755         }
0756       }
0757 
0758       iit->second.g->SetTitle(";rotation coefficient   (mm);residual   (mm)");
0759       iit->second.g->SetMarkerColor(idx + 1);
0760       iit->second.g->SetName(label.c_str());
0761       iit->second.g->Draw((idx == 0) ? "p" : "p");
0762       iit->second.g->Write();
0763 
0764       iit->second.h->SetName((label + " (hist)").c_str());
0765       iit->second.h->SetTitle(";rotation coefficient   (mm);residual   (mm)");
0766       iit->second.h->Write();
0767     }
0768 
0769     gDirectory->cd("..");
0770     c->Write();
0771   }
0772 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!