CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
List of all members | Classes | Public Member Functions | Private Attributes
JanAlignmentAlgorithm Class Reference

Jan's alignment algorithm. More...

#include <JanAlignmentAlgorithm.h>

Inheritance diagram for JanAlignmentAlgorithm:
AlignmentAlgorithm

Classes

struct  DetStat
 structure holding statistical information for one detector More...
 
struct  ScatterPlot
 a scatter plot, with graph and histogram representations More...
 

Public Member Functions

void analyze () override
 analyzes the data collected More...
 
void begin (const CTPPSGeometry *geometryReal, const CTPPSGeometry *geometryMisaligned) override
 prepare for processing More...
 
void end () override
 cleans up after processing More...
 
void feed (const HitCollection &, const LocalTrackFit &) override
 process one track More...
 
std::string getName () override
 
bool hasErrorEstimate () override
 returns whether this algorithm is capable of estimating result uncertainties More...
 
 JanAlignmentAlgorithm ()
 dummy constructor (not to be used) More...
 
 JanAlignmentAlgorithm (const edm::ParameterSet &ps, AlignmentTask *_t)
 normal constructor More...
 
void saveDiagnostics (TDirectory *) override
 saves diagnostic histograms/plots More...
 
unsigned int solve (const std::vector< AlignmentConstraint > &, std::map< unsigned int, AlignmentResult > &results, TDirectory *dir) override
 
 ~JanAlignmentAlgorithm () override
 
- Public Member Functions inherited from AlignmentAlgorithm
 AlignmentAlgorithm ()
 dummy constructor (not to be used) More...
 
 AlignmentAlgorithm (const edm::ParameterSet &ps, AlignmentTask *_t)
 normal constructor More...
 
virtual ~AlignmentAlgorithm ()
 

Private Attributes

bool buildDiagnosticPlots
 flag whether to build statistical plots More...
 
unsigned int events
 event count More...
 
TVectorD M
 final M vector More...
 
TVectorD * Mc
 
TMatrixD S
 final S matrix More...
 
TVectorD S_eigVal
 eigen values of the S matrix More...
 
TMatrixD S_eigVec
 matrix of S eigenvectors More...
 
TMatrixD ** Sc
 
std::vector< SingularModesingularModes
 a list of the singular modes of the S matrix More...
 
std::map< unsigned int, DetStatstatistics
 statistical data collection More...
 
bool stopOnSingularModes
 whether to stop when singular modes are identified More...
 
double weakLimit
 normalized eigen value below which the (CS) eigen vectors are considered as weak More...
 

Additional Inherited Members

- Protected Attributes inherited from AlignmentAlgorithm
double singularLimit
 eigenvalues in (-singularLimit, singularLimit) are treated as singular More...
 
AlignmentTasktask
 the tasked to be completed More...
 
unsigned int verbosity
 

Detailed Description

Jan's alignment algorithm.

Definition at line 25 of file JanAlignmentAlgorithm.h.

Constructor & Destructor Documentation

◆ JanAlignmentAlgorithm() [1/2]

JanAlignmentAlgorithm::JanAlignmentAlgorithm ( )
inline

dummy constructor (not to be used)

Definition at line 83 of file JanAlignmentAlgorithm.h.

83 {}

◆ JanAlignmentAlgorithm() [2/2]

JanAlignmentAlgorithm::JanAlignmentAlgorithm ( const edm::ParameterSet ps,
AlignmentTask _t 
)

normal constructor

Definition at line 28 of file JanAlignmentAlgorithm.cc.

References buildDiagnosticPlots, edm::ParameterSet::getParameter(), edm::ParameterSet::getParameterSet(), stopOnSingularModes, and weakLimit.

29  : AlignmentAlgorithm(ps, _t), Sc(nullptr), Mc(nullptr) {
30  const ParameterSet &lps = ps.getParameterSet("JanAlignmentAlgorithm");
31  weakLimit = lps.getParameter<double>("weakLimit");
32  stopOnSingularModes = lps.getParameter<bool>("stopOnSingularModes");
33  buildDiagnosticPlots = lps.getParameter<bool>("buildDiagnosticPlots");
34 }
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
ParameterSet const & getParameterSet(std::string const &) const
bool stopOnSingularModes
whether to stop when singular modes are identified
bool buildDiagnosticPlots
flag whether to build statistical plots
AlignmentAlgorithm()
dummy constructor (not to be used)
double weakLimit
normalized eigen value below which the (CS) eigen vectors are considered as weak

◆ ~JanAlignmentAlgorithm()

JanAlignmentAlgorithm::~JanAlignmentAlgorithm ( )
override

Definition at line 38 of file JanAlignmentAlgorithm.cc.

38 {}

Member Function Documentation

◆ analyze()

void JanAlignmentAlgorithm::analyze ( )
overridevirtual

analyzes the data collected

Implements AlignmentAlgorithm.

Definition at line 268 of file JanAlignmentAlgorithm.cc.

References change_name::diff, events, AlignmentTask::geometry, AlignmentGeometry::getNumberOfDetectors(), mps_fire::i, dqmiolumiharvest::j, M, reco::castor::maxDiff(), Mc, makeMEIFBenchmarkPlots::nev, hltrates_dqm_sourceclient-live_cfg::offset, AlignmentTask::quantityClasses, S_eigVal, S_eigVec, Sc, AlignmentAlgorithm::singularLimit, singularModes, AlignmentAlgorithm::task, heppy_batch::val, and AlignmentAlgorithm::verbosity.

268  {
269  if (verbosity > 2)
270  printf("\n>> JanAlignmentAlgorithm::Analyze\n");
271 
272  // calculate full dimension
273  unsigned int dim = 0;
274  for (unsigned int i = 0; i < task->quantityClasses.size(); i++)
275  dim += Mc[i].GetNrows();
276 
277  if (verbosity > 2) {
278  printf("\tsensors: %u\n", task->geometry.getNumberOfDetectors());
279  printf("\tfull dimension: %u\n", dim);
280  printf("\tquantity classes: %lu\n", task->quantityClasses.size());
281  }
282 
283  // build full M
284  M.ResizeTo(dim);
285  unsigned int offset = 0;
286  for (unsigned int i = 0; i < task->quantityClasses.size(); i++) {
287  M.SetSub(offset, Mc[i]);
288  offset += Mc[i].GetNrows();
289  }
290 
291  // build full S
292  S.ResizeTo(dim, dim);
293  unsigned int r_offset = 0, c_offset = 0;
294  for (unsigned int i = 0; i < task->quantityClasses.size(); i++) {
295  c_offset = 0;
296  unsigned int r_size = 0, c_size = 0;
297  for (unsigned int j = 0; j < task->quantityClasses.size(); j++) {
298  r_size = Sc[i][j].GetNrows();
299  c_size = Sc[i][j].GetNcols();
300 
301  if (r_size < 1 || c_size < 1)
302  continue;
303 
304  TMatrixDSub(S, r_offset, r_offset + r_size - 1, c_offset, c_offset + c_size - 1) = Sc[i][j];
305  c_offset += c_size;
306  }
307  r_offset += r_size;
308  }
309 
310  // analyze symmetricity
311  if (verbosity >= 3) {
312  double maxDiff = 0., maxElem = 0.;
313  for (unsigned int i = 0; i < dim; i++) {
314  for (unsigned int j = 0; j < dim; j++) {
315  double diff = S[i][j] - S[j][i];
316  if (fabs(diff) > maxDiff)
317  maxDiff = diff;
318  if (S[i][j] > maxElem)
319  maxElem = S[i][j];
320  }
321  }
322 
323  printf("\n* S matrix:\n\tdimension = %i\n\tmaximum asymmetry: %E\t(ratio to maximum element %E)\n",
324  dim,
325  maxDiff,
326  maxDiff / maxElem);
327  }
328 
329  // make a symmetric copy
330  TMatrixDSym S_sym(dim);
331  for (unsigned int j = 0; j < dim; j++) {
332  for (unsigned int i = 0; i < dim; i++) {
333  S_sym[i][j] = S[i][j];
334  }
335  }
336 
337  // eigen analysis of S
338  TMatrixDSymEigen S_eig(S_sym);
339  const TVectorD &S_eigVal_temp = S_eig.GetEigenValues();
340  S_eigVal.ResizeTo(S_eigVal_temp.GetNrows());
341  S_eigVal = S_eigVal_temp;
342  const TMatrixD &S_eigVec_temp = S_eig.GetEigenVectors();
343  S_eigVec.ResizeTo(S_eigVec_temp);
344  S_eigVec = S_eigVec_temp;
345 
346  // identify singular modes
347  for (int i = 0; i < S_eigVal.GetNrows(); i++) {
348  double nev = S_eigVal[i] / events;
349  if (fabs(nev) < singularLimit) {
350  SingularMode sM{S_eigVal[i], TMatrixDColumn(S_eigVec, i), i};
351  singularModes.push_back(sM);
352  }
353  }
354 
355 #if DEBUG
356  // print singular vectors
357  if (singularModes.size() > 0) {
358  printf("\n* S singular modes\n | ");
359  for (unsigned int i = 0; i < singularModes.size(); i++)
360  printf("%+10.3E ", singularModes[i].val);
361  printf("\n-- | ");
362 
363  for (unsigned int i = 0; i < singularModes.size(); i++)
364  printf("---------- ");
365  printf("\n");
366 
367  for (unsigned int j = 0; j < dim; j++) {
368  printf("%2u | ", j);
369  for (unsigned int i = 0; i < singularModes.size(); i++) {
370  printf("%+10.3E ", singularModes[i].vec[j]);
371  }
372  printf("\n");
373  }
374  } else
375  printf("\n* S has no singular modes\n");
376 #endif
377 }
AlignmentTask * task
the tasked to be completed
TMatrixD S_eigVec
matrix of S eigenvectors
TVectorD S_eigVal
eigen values of the S matrix
std::vector< QuantityClass > quantityClasses
list of quantity classes to be optimized
Definition: AlignmentTask.h:68
AlignmentGeometry geometry
the geometry for this task
Definition: AlignmentTask.h:48
double singularLimit
eigenvalues in (-singularLimit, singularLimit) are treated as singular
float maxDiff(float one, float two, float three, float four)
unsigned int getNumberOfDetectors() const
returns the number of detectors in the collection
TVectorD M
final M vector
std::vector< SingularMode > singularModes
a list of the singular modes of the S matrix
unsigned int events
event count

◆ begin()

void JanAlignmentAlgorithm::begin ( const CTPPSGeometry geometryReal,
const CTPPSGeometry geometryMisaligned 
)
overridevirtual

prepare for processing

Implements AlignmentAlgorithm.

Definition at line 42 of file JanAlignmentAlgorithm.cc.

References visDQMUpload::buf, buildDiagnosticPlots, HltBtagPostValidation_cff::c, events, g, AlignmentTask::geometry, AlignmentGeometry::getSensorMap(), mps_fire::i, l1ctLayer2EG_cff::id, ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, dqmiolumiharvest::j, Mc, AlignmentTask::quantitiesOfClass(), AlignmentTask::quantityClasses, AlignmentTask::quantityClassTag(), mysort::rows, alignCSCRings::s, Sc, statistics, and AlignmentAlgorithm::task.

42  {
43  // initialize M and S components
44  Mc = new TVectorD[task->quantityClasses.size()];
45  Sc = new TMatrixD *[task->quantityClasses.size()];
46  for (unsigned int i = 0; i < task->quantityClasses.size(); i++) {
47  unsigned int rows = task->quantitiesOfClass(task->quantityClasses[i]);
48 
49  Mc[i].ResizeTo(rows);
50  Mc[i].Zero();
51 
52  Sc[i] = new TMatrixD[task->quantityClasses.size()];
53  for (unsigned int j = 0; j < task->quantityClasses.size(); j++) {
54  unsigned int cols = task->quantitiesOfClass(task->quantityClasses[j]);
55  Sc[i][j].ResizeTo(rows, cols);
56  Sc[i][j].Zero();
57  }
58  }
59 
60  // prepare statistics plots
62  for (const auto &it : task->geometry.getSensorMap()) {
63  unsigned int id = it.first;
64  char buf[50];
65  DetStat s;
66 
67  sprintf(buf, "%u: m distribution", id);
68  s.m_dist = new TH1D(buf, ";u or v (mm)", 100, -25., 25.);
69 
70  sprintf(buf, "%u: R distribution", id);
71  s.R_dist = new TH1D(buf, ";R (mm)", 500, -0.5, 0.5);
72 
73  for (unsigned int c = 0; c < task->quantityClasses.size(); c++) {
74  sprintf(buf, "%u: coef, %s", id, task->quantityClassTag(task->quantityClasses[c]).c_str());
75  s.coefHist.push_back(new TH1D(buf, ";coefficient", 100, -2., +2.));
76 
77  sprintf(buf, "%u: R vs. coef, %s", id, task->quantityClassTag(task->quantityClasses[c]).c_str());
78  TGraph *g = new TGraph();
79  g->SetName(buf);
80  g->SetTitle(";coefficient;residual (mm)");
81  s.resVsCoef.push_back(g);
82  }
83 
84  statistics[id] = s;
85  }
86  }
87 
88  events = 0;
89 }
AlignmentTask * task
the tasked to be completed
std::string quantityClassTag(QuantityClass) const
returns a string tag for the given quantity class
const std::map< unsigned int, DetGeometry > & getSensorMap() const
std::vector< QuantityClass > quantityClasses
list of quantity classes to be optimized
Definition: AlignmentTask.h:68
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
Definition: Activities.doc:4
bool buildDiagnosticPlots
flag whether to build statistical plots
AlignmentGeometry geometry
the geometry for this task
Definition: AlignmentTask.h:48
unsigned int quantitiesOfClass(QuantityClass) const
returns the number of quantities of the given class
std::map< unsigned int, DetStat > statistics
statistical data collection
unsigned int events
event count
rows
Definition: mysort.py:12

◆ end()

void JanAlignmentAlgorithm::end ( )
overridevirtual

cleans up after processing

Implements AlignmentAlgorithm.

Definition at line 707 of file JanAlignmentAlgorithm.cc.

References mps_fire::i, Mc, AlignmentTask::quantityClasses, Sc, and AlignmentAlgorithm::task.

Referenced by Types.LuminosityBlockRange::cppID(), and Types.EventRange::cppID().

707  {
708  delete[] Mc;
709 
710  for (unsigned int i = 0; i < task->quantityClasses.size(); i++) {
711  delete[] Sc[i];
712  }
713  delete[] Sc;
714 }
AlignmentTask * task
the tasked to be completed
std::vector< QuantityClass > quantityClasses
list of quantity classes to be optimized
Definition: AlignmentTask.h:68

◆ feed()

void JanAlignmentAlgorithm::feed ( const HitCollection ,
const LocalTrackFit  
)
overridevirtual

process one track

Implements AlignmentAlgorithm.

Definition at line 93 of file JanAlignmentAlgorithm.cc.

References A, LocalTrackFit::ax, LocalTrackFit::ay, buildDiagnosticPlots, LocalTrackFit::bx, LocalTrackFit::by, correctionTermsCaloMet_cff::C, HltBtagPostValidation_cff::c, ztail::d, hcalRecHitTable_cff::detId, events, JanAlignmentAlgorithm::ScatterPlot::g, AlignmentTask::geometry, AlignmentGeometry::get(), AlignmentTask::getMeasurementIndex(), AlignmentTask::getQuantityIndex(), JanAlignmentAlgorithm::ScatterPlot::h, mps_fire::i, l1ctLayer2EG_cff::id, hit::id, dqmiolumiharvest::j, visualization-live-secondInstance_cfg::m, Mc, AlignmentTask::qcRotZ, AlignmentTask::qcShR1, AlignmentTask::qcShR2, AlignmentTask::qcShZ, AlignmentTask::quantityClasses, dttmaxenums::R, alignCSCRings::r, S, alignCSCRings::s, Sc, corrVsCorr::selection, mathSSE::sqrt(), statistics, submitPVValidationJobs::t, AlignmentAlgorithm::task, groupFilesInBlocks::tt, AlignmentAlgorithm::verbosity, hit::z, LocalTrackFit::z0, and AlignmentGeometry::z0.

Referenced by Page1Parser.Page1Parser::_Parse().

93  {
94  if (verbosity > 9)
95  printf("\n>> JanAlignmentAlgorithm::Feed\n");
96 
97  events++;
98 
99  // prepare fit - make z0 compatible
100  double hax = trackFit.ax;
101  double hay = trackFit.ay;
102  double hbx = trackFit.bx + trackFit.ax * (task->geometry.z0 - trackFit.z0);
103  double hby = trackFit.by + trackFit.ay * (task->geometry.z0 - trackFit.z0);
104 
105  // prepare Gamma matrices (full of zeros)
106  TMatrixD *Ga = new TMatrixD[task->quantityClasses.size()];
107  for (unsigned int i = 0; i < task->quantityClasses.size(); i++) {
108  Ga[i].ResizeTo(selection.size(), Mc[i].GetNrows());
109  Ga[i].Zero();
110  }
111 
112  TMatrixD A(selection.size(), 4);
113  TMatrixD Vi(selection.size(), selection.size());
114  TVectorD m(selection.size());
115 
116  set<unsigned int> rpSet;
117  if (buildDiagnosticPlots) {
118  for (const auto &hit : selection) {
120  const unsigned int rpDecId = 100 * detId.arm() + 10 * detId.station() + detId.rp();
121  rpSet.insert(rpDecId);
122  }
123  }
124 
125  // fill fit matrix and Gamma matrices
126  unsigned int j = 0;
127 
128  for (HitCollection::const_iterator hit = selection.begin(); hit != selection.end(); ++hit, ++j) {
129  unsigned int id = hit->id;
130 
131  const DetGeometry &d = task->geometry.get(id);
132  const auto &dirData = d.getDirectionData(hit->dirIdx);
133 
134  A(j, 0) = hit->z * dirData.dx;
135  A(j, 1) = dirData.dx;
136  A(j, 2) = hit->z * dirData.dy;
137  A(j, 3) = dirData.dy;
138 
139  m(j) = hit->position + dirData.s - (hit->z - d.z) * dirData.dz; // in mm
140 
141  Vi(j, j) = 1. / hit->sigma / hit->sigma;
142 
143  double C = dirData.dx, S = dirData.dy;
144 
145  double hx = hax * hit->z + hbx; // in mm
146  double hy = hay * hit->z + hby;
147  double R = m(j) - (hx * C + hy * S); // (standard) residual
148 
149  if (buildDiagnosticPlots) {
150  statistics[id].m_dist->Fill(m(j));
151  statistics[id].R_dist->Fill(R);
152  }
153 
154  for (unsigned int i = 0; i < task->quantityClasses.size(); i++) {
155  // check compatibility
156  signed int matrixIndex = task->getMeasurementIndex(task->quantityClasses[i], hit->id, hit->dirIdx);
157  if (matrixIndex < 0)
158  continue;
159 
160  matrixIndex = task->getQuantityIndex(task->quantityClasses[i], hit->id);
161 
162  switch (task->quantityClasses[i]) {
164  Ga[i][j][matrixIndex] = -1.;
165  break;
166 
168  Ga[i][j][matrixIndex] = -1.;
169  break;
170 
172  Ga[i][j][matrixIndex] = hax * C + hay * S;
173  break;
174 
176  Ga[i][j][matrixIndex] = (hax * hit->z + hbx - d.sx) * (-S) + (hay * hit->z + hby - d.sy) * C;
177  break;
178  }
179 
180  if (buildDiagnosticPlots) {
181  double c = Ga[i][j][matrixIndex];
182  DetStat &s = statistics[id];
183  s.coefHist[i]->Fill(c);
184  s.resVsCoef[i]->SetPoint(s.resVsCoef[i]->GetN(), c, R);
185 
187  map<set<unsigned int>, ScatterPlot>::iterator hit = s.resVsCoefRot_perRPSet.find(rpSet);
188  if (hit == s.resVsCoefRot_perRPSet.end()) {
189  ScatterPlot sp;
190  sp.g = new TGraph();
191  sp.h = new TH2D("", "", 40, -20., +20., 60, -0.15, +0.15);
192  hit = s.resVsCoefRot_perRPSet.insert(pair<set<unsigned int>, ScatterPlot>(rpSet, sp)).first;
193  }
194  hit->second.g->SetPoint(hit->second.g->GetN(), c, R);
195  hit->second.h->Fill(c, R);
196  }
197  }
198  }
199  }
200 
201  // sigma matrix
202  TMatrixD AT(TMatrixD::kTransposed, A);
203  TMatrixD ATViA(4, 4);
204  ATViA = AT * Vi * A;
205  TMatrixD ATViAI(ATViA);
206  ATViAI = ATViA.Invert();
207  TMatrixD sigma(Vi);
208  sigma -= Vi * A * ATViAI * AT * Vi;
209 
210  // traspose Gamma matrices
211  TMatrixD *GaT = new TMatrixD[task->quantityClasses.size()];
212  for (unsigned int i = 0; i < task->quantityClasses.size(); i++) {
213  GaT[i].ResizeTo(Mc[i].GetNrows(), selection.size());
214  GaT[i].Transpose(Ga[i]);
215  }
216 
217  // normalized residuals
218  TVectorD r(selection.size());
219  r = sigma * m;
220 
221  // increment M
222  for (unsigned int i = 0; i < task->quantityClasses.size(); i++) {
223  if (Mc[i].GetNrows() < 1)
224  continue;
225 
226  Mc[i] += GaT[i] * r;
227  }
228 
229  // increment S
230  for (unsigned int i = 0; i < task->quantityClasses.size(); i++) {
231  for (unsigned int j = 0; j < task->quantityClasses.size(); j++) {
232  if (Sc[i][j].GetNrows() < 1 || Sc[i][j].GetNcols() < 1)
233  continue;
234 
235  Sc[i][j] += GaT[i] * sigma * Ga[j];
236  }
237  }
238 
239 #ifdef DEBUG
240  printf("* checking normalized residuals, selection.size = %u\n", selection.size());
241  r.Print();
242 
243  for (unsigned int i = 0; i < task->quantityClasses.size(); i++) {
244  printf("- class %u\n", i);
245  TVectorD t(Mc[i].GetNrows());
246  for (int j = 0; j < t.GetNrows(); j++)
247  t[j] = 1.;
248  t.Print();
249 
250  Ga[i].Print();
251 
252  TVectorD tt(selection.size());
253  tt = sigma * Ga[i] * t;
254 
255  double ttn = sqrt(tt.Norm2Sqr());
256  printf("|tt| = %E\n", ttn);
257  if (ttn > 1E-8)
258  tt.Print();
259  }
260 #endif
261 
262  delete[] Ga;
263  delete[] GaT;
264 }
AlignmentTask * task
the tasked to be completed
A structure to hold relevant geometrical information about one detector/sensor.
selection
main part
Definition: corrVsCorr.py:100
TMatrixD S
final S matrix
signed int getMeasurementIndex(QuantityClass cl, unsigned int detId, unsigned int dirIdx) const
returns measurement index (if non-existent, returns -1)
std::vector< QuantityClass > quantityClasses
list of quantity classes to be optimized
Definition: AlignmentTask.h:68
Definition: TTTypes.h:54
detector shifts in first readout direction
Definition: AlignmentTask.h:61
bool buildDiagnosticPlots
flag whether to build statistical plots
AlignmentGeometry geometry
the geometry for this task
Definition: AlignmentTask.h:48
T sqrt(T t)
Definition: SSEVec.h:19
d
Definition: ztail.py:151
detector shifts in second readout direction
Definition: AlignmentTask.h:62
unsigned int id
signed int getQuantityIndex(QuantityClass cl, unsigned int detId) const
returns measurement index (if non-existent, returns -1)
detector rotations around z
Definition: AlignmentTask.h:64
Base class for CTPPS detector IDs.
Definition: CTPPSDetId.h:32
const DetGeometry & get(unsigned int id) const
retrieves sensor geometry
std::map< unsigned int, DetStat > statistics
statistical data collection
Definition: APVGainStruct.h:7
double z0
a characteristic z in mm
unsigned int events
event count
detector shifts in z
Definition: AlignmentTask.h:63

◆ getName()

std::string JanAlignmentAlgorithm::getName ( )
inlineoverridevirtual

Reimplemented from AlignmentAlgorithm.

Definition at line 90 of file JanAlignmentAlgorithm.h.

Referenced by plotting.Plot::draw().

90 { return "Jan"; }

◆ hasErrorEstimate()

bool JanAlignmentAlgorithm::hasErrorEstimate ( )
inlineoverridevirtual

returns whether this algorithm is capable of estimating result uncertainties

Implements AlignmentAlgorithm.

Definition at line 92 of file JanAlignmentAlgorithm.h.

92 { return true; }

◆ saveDiagnostics()

void JanAlignmentAlgorithm::saveDiagnostics ( TDirectory *  )
overridevirtual

saves diagnostic histograms/plots

Implements AlignmentAlgorithm.

Definition at line 718 of file JanAlignmentAlgorithm.cc.

References visDQMUpload::buf, buildDiagnosticPlots, HltBtagPostValidation_cff::c, DeadROC_duringRun::dir, dqmdumpme::first, heavyIonCSV_trainingSettings::idx, ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, label, genParticles_cff::map, AlignmentTask::quantityClasses, statistics, and AlignmentAlgorithm::task.

718  {
720  return;
721 
722  for (map<unsigned int, DetStat>::iterator it = statistics.begin(); it != statistics.end(); ++it) {
723  char buf[50];
724  sprintf(buf, "%u", it->first);
725  gDirectory = dir->mkdir(buf);
726 
727  it->second.m_dist->Write();
728  it->second.R_dist->Write();
729 
730  for (unsigned int c = 0; c < task->quantityClasses.size(); c++) {
731  it->second.coefHist[c]->Write();
732  it->second.resVsCoef[c]->Write();
733  }
734 
735  gDirectory = gDirectory->mkdir("R vs. rot. coef, per RP set");
736  TCanvas *c = new TCanvas;
737  c->SetName("R vs. rot. coef, overlapped");
738  TH2D *frame = new TH2D("frame", "frame", 100, -20., +20., 100, -0.15, +0.15);
739  frame->Draw();
740  unsigned int idx = 0;
741  for (map<set<unsigned int>, ScatterPlot>::iterator iit = it->second.resVsCoefRot_perRPSet.begin();
742  iit != it->second.resVsCoefRot_perRPSet.end();
743  ++iit, ++idx) {
744  string label;
745  bool first = true;
746  for (set<unsigned int>::iterator sit = iit->first.begin(); sit != iit->first.end(); ++sit) {
747  char buf[50];
748  sprintf(buf, "%u", *sit);
749 
750  if (first) {
751  label = buf;
752  first = false;
753  } else {
754  label = label + ", " + buf;
755  }
756  }
757 
758  iit->second.g->SetTitle(";rotation coefficient (mm);residual (mm)");
759  iit->second.g->SetMarkerColor(idx + 1);
760  iit->second.g->SetName(label.c_str());
761  iit->second.g->Draw((idx == 0) ? "p" : "p");
762  iit->second.g->Write();
763 
764  iit->second.h->SetName((label + " (hist)").c_str());
765  iit->second.h->SetTitle(";rotation coefficient (mm);residual (mm)");
766  iit->second.h->Write();
767  }
768 
769  gDirectory->cd("..");
770  c->Write();
771  }
772 }
AlignmentTask * task
the tasked to be completed
std::vector< QuantityClass > quantityClasses
list of quantity classes to be optimized
Definition: AlignmentTask.h:68
char const * label
bool buildDiagnosticPlots
flag whether to build statistical plots
std::map< unsigned int, DetStat > statistics
statistical data collection

◆ solve()

unsigned int JanAlignmentAlgorithm::solve ( const std::vector< AlignmentConstraint > &  ,
std::map< unsigned int, AlignmentResult > &  results,
TDirectory *  dir 
)
overridevirtual

solves the alignment problem with the given constraints

Parameters
dira directory (in StraightTrackAlignment::taskDataFileName) where intermediate results can be stored

Implements AlignmentAlgorithm.

Definition at line 381 of file JanAlignmentAlgorithm.cc.

References correctionTermsCaloMet_cff::C, mps_check::columns, createBeamHaloJobs::constraints, cuy::cv, DeadROC_duringRun::dir, MillePedeFileConverter_cfg::e, events, dqmdumpme::first, AlignmentTask::geometry, AlignmentTask::getQuantityIndex(), AlignmentGeometry::getSensorMap(), mps_fire::i, heavyIonCSV_trainingSettings::idx, dqmiolumiharvest::j, dqmdumpme::k, dqmdumpme::last, M, SiStripPI::max, reco::castor::maxDiff(), Mc, Skims_PA_cff::name, hltrates_dqm_sourceclient-live_cfg::offset, unpackBuffers-CaloStage1::offsets, print(), AlignmentTask::qcRotZ, AlignmentTask::qcShR1, AlignmentTask::qcShR2, AlignmentTask::qcShZ, AlignmentTask::quantityClasses, alignCSCRings::r, mysort::results, S_eigVal, S_eigVec, AlignmentAlgorithm::singularLimit, singularModes, mathSSE::sqrt(), stopOnSingularModes, AlignmentAlgorithm::task, mitigatedMETSequence_cff::U, findQualityFiles::v, heppy_batch::val, AlignmentAlgorithm::verbosity, and weakLimit.

383  {
384  if (verbosity)
385  printf(">> JanAlignmentAlgorithm::Solve\n");
386 
387  results.clear();
388 
389  // build C matrix
390  unsigned int dim = S.GetNrows();
391  TMatrixD C(dim, constraints.size());
392  TMatrixD C2(dim, constraints.size());
393  for (unsigned int i = 0; i < constraints.size(); i++) {
394  unsigned int offset = 0;
395  for (auto &quantityClass : task->quantityClasses) {
396  const TVectorD &cv = constraints[i].coef.find(quantityClass)->second;
397  for (int k = 0; k < cv.GetNrows(); k++) {
398  C[offset][i] = events * cv[k];
399  C2[offset][i] = events * cv[k] * 1E3;
400  offset++;
401  }
402  }
403  }
404 
405 #ifdef DEBUG
406  printf("\n* constraint matrix\n");
407  Print(C);
408 #endif
409 
410  // build E matrix (singular vectors of S as its columns)
411  TMatrixD E(S.GetNrows(), singularModes.size());
412  for (unsigned int i = 0; i < singularModes.size(); i++)
413  for (int j = 0; j < S.GetNrows(); j++)
414  E(j, i) = singularModes[i].vec[j];
415 
416  // build CS matrix
417  TMatrixDSym CS(dim + constraints.size());
418  TMatrixDSym CS2(dim + constraints.size());
419  CS.Zero();
420  CS2.Zero();
421  for (unsigned int j = 0; j < dim; j++) {
422  for (unsigned int i = 0; i < dim; i++) {
423  CS[i][j] = S[i][j];
424  CS2[i][j] = S[i][j];
425  }
426  }
427 
428  for (unsigned int i = 0; i < constraints.size(); i++) {
429  for (unsigned int j = 0; j < dim; j++) {
430  CS[j][dim + i] = CS[dim + i][j] = C(j, i);
431  CS2[j][dim + i] = CS2[dim + i][j] = C2(j, i);
432  }
433  }
434 
435  // eigen analysis of CS matrix
436  TMatrixDSymEigen CS_eig(CS);
437  TVectorD CS_eigVal = CS_eig.GetEigenValues();
438  TMatrixD CS_eigVec = CS_eig.GetEigenVectors();
439 
440  // check regularity of CS matrix
441  if (verbosity >= 2) {
442  printf("\n* eigen values of CS and S matrices (events = %u)\n", events);
443  printf(" # CS norm. CS S norm. S\n");
444  }
445 
446  unsigned int singularModeCount = 0;
447  vector<unsigned int> weakModeIdx;
448  for (int i = 0; i < CS_eigVal.GetNrows(); i++) {
449  const double CS_nev = CS_eigVal[i] / events;
450 
451  if (fabs(CS_nev) < singularLimit)
452  singularModeCount++;
453 
454  if (verbosity >= 2) {
455  printf("%4i%+12.2E%+12.2E", i, CS_eigVal[i], CS_nev);
456  if (fabs(CS_nev) < singularLimit) {
457  singularModeCount++;
458  printf(" (S)");
459  } else {
460  if (fabs(CS_nev) < weakLimit) {
461  weakModeIdx.push_back(i);
462  printf(" (W)");
463  } else {
464  printf(" ");
465  }
466  }
467 
468  if (i < S_eigVal.GetNrows()) {
469  double S_nev = S_eigVal[i] / events;
470  printf("%+12.2E%+12.2E", S_eigVal[i], S_nev);
471  if (fabs(S_nev) < singularLimit)
472  printf(" (S)");
473  else if (fabs(S_nev) < weakLimit)
474  printf(" (W)");
475  }
476 
477  printf("\n");
478  }
479  }
480 
481  if (verbosity >= 2) {
482  // print weak vectors
483  if (!weakModeIdx.empty()) {
484  unsigned int columns = 10;
485  unsigned int first = 0;
486 
487  while (first < weakModeIdx.size()) {
488  unsigned int last = first + columns;
489  if (last >= weakModeIdx.size())
490  last = weakModeIdx.size();
491 
492  printf("\n* CS weak modes\n | ");
493  for (unsigned int i = first; i < last; i++)
494  printf("%+10.3E ", CS_eigVal[weakModeIdx[i]]);
495  printf("\n--- | ");
496 
497  for (unsigned int i = first; i < last; i++)
498  printf("---------- ");
499  printf("\n");
500 
501  // determine maximum elements
502  vector<double> maxs;
503  for (unsigned int i = first; i < last; i++) {
504  double max = 0;
505  for (unsigned int j = 0; j < dim + constraints.size(); j++) {
506  double v = fabs(CS_eigVec(weakModeIdx[i], j));
507  if (v > max)
508  max = v;
509  }
510  maxs.push_back(max);
511  }
512 
513  for (unsigned int j = 0; j < dim + constraints.size(); j++) {
514  printf("%3u | ", j);
515  for (unsigned int i = first; i < last; i++) {
516  double v = CS_eigVec(weakModeIdx[i], j);
517  if (fabs(v) / maxs[i - first] > 1E-3)
518  printf("%+10.3E ", v);
519  else
520  printf(" . ");
521  }
522  printf("\n");
523  }
524 
525  first = last;
526  }
527  } else
528  printf("\n* CS has no weak modes\n");
529  }
530 
531  // check the regularity of C^T E
532  if (verbosity >= 2) {
533  if (E.GetNcols() == C.GetNcols()) {
534  TMatrixD CTE(C, TMatrixD::kTransposeMult, E);
535  print(CTE, "* CTE matrix:");
536  const double &det = CTE.Determinant();
537  printf(
538  "\n* det(CTE) = %E, max(CTE) = %E, det(CTE)/max(CTE) = %E\n\tmax(C) = %E, max(E) = %E, "
539  "det(CTE)/max(C)/max(E) = %E\n",
540  det,
541  CTE.Max(),
542  det / CTE.Max(),
543  C.Max(),
544  E.Max(),
545  det / C.Max() / E.Max());
546  } else {
547  printf(">> JanAlignmentAlgorithm::Solve > WARNING: C matrix has %u, while E matrix %u columns.\n",
548  C.GetNcols(),
549  E.GetNcols());
550  }
551  }
552 
553  // stop if CS is singular
554  if (singularModeCount > 0 && stopOnSingularModes) {
555  LogError("PPS") << "\n>> JanAlignmentAlgorithm::Solve > ERROR: There are " << singularModeCount
556  << " singular modes in CS matrix. Stopping.";
557  return 1;
558  }
559 
560  // build MV vector
561  TVectorD MV(dim + constraints.size());
562  for (unsigned int i = 0; i < dim; i++)
563  MV[i] = M[i];
564  for (unsigned int i = 0; i < constraints.size(); i++)
565  MV[dim + i] = events * constraints[i].val;
566 
567  // perform inversion and solution
568  TMatrixD CSI(TMatrixD::kInverted, CS);
569  TMatrixD CS2I(TMatrixD::kInverted, CS2);
570  TVectorD AL(MV);
571  AL = CSI * MV;
572 
573  // evaluate error matrix
574  TMatrixD S0(S); // new parts full of zeros
575  S0.ResizeTo(dim + constraints.size(), dim + constraints.size());
576  TMatrixD EM(CS);
577  EM = CSI * S0 * CSI;
578 
579  TMatrixD EM2(CS2);
580  EM2 = CS2I * S0 * CS2I;
581 
582  TMatrixD EMdiff(EM2 - EM);
583 
584  if (verbosity >= 3) {
585  double max1 = -1., max2 = -1., maxDiff = -1.;
586  for (int i = 0; i < EMdiff.GetNrows(); i++) {
587  for (int j = 0; j < EMdiff.GetNcols(); j++) {
588  if (maxDiff < fabs(EMdiff(i, j)))
589  maxDiff = fabs(EMdiff(i, j));
590 
591  if (max1 < fabs(EM(i, j)))
592  max1 = fabs(EM(i, j));
593 
594  if (max2 < fabs(EM2(i, j)))
595  max2 = fabs(EM2(i, j));
596  }
597  }
598 
599  printf("EM max = %E, EM2 max = %E, EM2 - EM max = %E\n", max1, max2, maxDiff);
600  }
601 
602  // tests
603  TMatrixD &U = CS_eigVec;
604  TMatrixD UT(TMatrixD::kTransposed, U);
605  //TMatrixD CSEi(CS);
606  //CSEi = UT * CS * U;
607  //Print(CSEi, "CSEi");
608 
609  TMatrixD EMEi(EM);
610  EMEi = UT * EM * U;
611  //Print(EMEi, "*EMEi");
612 
613  if (verbosity >= 3) {
614  double max = -1.;
615  for (int i = 0; i < EMEi.GetNrows(); i++) {
616  for (int j = 0; j < EMEi.GetNcols(); j++) {
617  if (max < EMEi(i, j))
618  max = EMEi(i, j);
619  }
620  }
621 
622  printf("max = %E\n", max);
623  }
624 
625  // print lambda values
626  if (verbosity >= 3) {
627  printf("\n* Lambda (from the contribution of singular modes to MV)\n");
628  for (unsigned int i = 0; i < constraints.size(); i++) {
629  printf("\t%u (%25s)\t%+10.1E +- %10.1E\n",
630  i,
631  constraints[i].name.c_str(),
632  AL[dim + i] * 1E3,
633  sqrt(EM[i + dim][i + dim]) * 1E3);
634  }
635  }
636 
637  // fill results
638  unsigned int offset = 0;
639  vector<unsigned int> offsets;
640  for (unsigned int i = 0; i < task->quantityClasses.size(); i++) {
641  offsets.push_back(offset);
642  offset += Mc[i].GetNrows();
643  }
644 
645  for (const auto &dit : task->geometry.getSensorMap()) {
647 
648  for (unsigned int i = 0; i < task->quantityClasses.size(); i++) {
649  signed idx = task->getQuantityIndex(task->quantityClasses[i], dit.first);
650  if (idx < 0)
651  continue;
652 
653  unsigned int fi = offsets[i] + idx;
654  double v = AL[fi];
655  double e = sqrt(EM[fi][fi]);
656  switch (task->quantityClasses[i]) {
658  r.setShR1(v);
659  r.setShR1Unc(e);
660  break;
662  r.setShR2(v);
663  r.setShR2Unc(e);
664  break;
666  r.setShZ(v);
667  r.setShZUnc(e);
668  break;
670  r.setRotZ(v);
671  r.setRotZUnc(e);
672  break;
673  }
674  }
675 
676  results[dit.first] = r;
677  }
678 
679  // save matrices, eigen data, ...
680  if (dir) {
681  dir->cd();
682 
683  S.Write("S");
684  S_eigVal.Write("S_eigen_values");
685  S_eigVec.Write("S_eigen_vectors");
686 
687  E.Write("E");
688  C.Write("C");
689 
690  CS.Write("CS");
691  CS_eigVal.Write("CS_eigen_values");
692  CS_eigVec.Write("CS_eigen_vectors");
693 
694  MV.Write("MV");
695  AL.Write("AL");
696 
697  S0.Write("S0");
698  EM.Write("EM");
699  }
700 
701  // success
702  return 0;
703 }
AlignmentTask * task
the tasked to be completed
TMatrixD S_eigVec
matrix of S eigenvectors
TVectorD S_eigVal
eigen values of the S matrix
const std::map< unsigned int, DetGeometry > & getSensorMap() const
cv
Definition: cuy.py:363
Log< level::Error, false > LogError
std::vector< QuantityClass > quantityClasses
list of quantity classes to be optimized
Definition: AlignmentTask.h:68
bool stopOnSingularModes
whether to stop when singular modes are identified
detector shifts in first readout direction
Definition: AlignmentTask.h:61
AlignmentGeometry geometry
the geometry for this task
Definition: AlignmentTask.h:48
T sqrt(T t)
Definition: SSEVec.h:19
void print(TMatrixD &m, const char *label=nullptr, bool mathematicaFormat=false)
Definition: Utilities.cc:47
detector shifts in second readout direction
Definition: AlignmentTask.h:62
double singularLimit
eigenvalues in (-singularLimit, singularLimit) are treated as singular
signed int getQuantityIndex(QuantityClass cl, unsigned int detId) const
returns measurement index (if non-existent, returns -1)
detector rotations around z
Definition: AlignmentTask.h:64
float maxDiff(float one, float two, float three, float four)
Result of CTPPS track-based alignment.
results
Definition: mysort.py:8
TVectorD M
final M vector
std::vector< SingularMode > singularModes
a list of the singular modes of the S matrix
unsigned int events
event count
detector shifts in z
Definition: AlignmentTask.h:63
double weakLimit
normalized eigen value below which the (CS) eigen vectors are considered as weak

Member Data Documentation

◆ buildDiagnosticPlots

bool JanAlignmentAlgorithm::buildDiagnosticPlots
private

flag whether to build statistical plots

Definition at line 79 of file JanAlignmentAlgorithm.h.

Referenced by begin(), feed(), JanAlignmentAlgorithm(), and saveDiagnostics().

◆ events

unsigned int JanAlignmentAlgorithm::events
private

◆ M

TVectorD JanAlignmentAlgorithm::M
private

final M vector

Definition at line 55 of file JanAlignmentAlgorithm.h.

Referenced by analyze(), and solve().

◆ Mc

TVectorD* JanAlignmentAlgorithm::Mc
private

M vector components indeces correspond to the qClToOpt list

Definition at line 49 of file JanAlignmentAlgorithm.h.

Referenced by analyze(), begin(), end(), feed(), and solve().

◆ S

TMatrixD JanAlignmentAlgorithm::S
private

final S matrix

Definition at line 52 of file JanAlignmentAlgorithm.h.

Referenced by feed().

◆ S_eigVal

TVectorD JanAlignmentAlgorithm::S_eigVal
private

eigen values of the S matrix

Definition at line 58 of file JanAlignmentAlgorithm.h.

Referenced by analyze(), and solve().

◆ S_eigVec

TMatrixD JanAlignmentAlgorithm::S_eigVec
private

matrix of S eigenvectors

Definition at line 61 of file JanAlignmentAlgorithm.h.

Referenced by analyze(), and solve().

◆ Sc

TMatrixD** JanAlignmentAlgorithm::Sc
private

S matrix components indeces correspond to the qClToOpt list

Definition at line 45 of file JanAlignmentAlgorithm.h.

Referenced by analyze(), begin(), end(), and feed().

◆ singularModes

std::vector<SingularMode> JanAlignmentAlgorithm::singularModes
private

a list of the singular modes of the S matrix

Definition at line 64 of file JanAlignmentAlgorithm.h.

Referenced by analyze(), and solve().

◆ statistics

std::map<unsigned int, DetStat> JanAlignmentAlgorithm::statistics
private

statistical data collection

Definition at line 76 of file JanAlignmentAlgorithm.h.

Referenced by begin(), feed(), and saveDiagnostics().

◆ stopOnSingularModes

bool JanAlignmentAlgorithm::stopOnSingularModes
private

whether to stop when singular modes are identified

Definition at line 67 of file JanAlignmentAlgorithm.h.

Referenced by JanAlignmentAlgorithm(), and solve().

◆ weakLimit

double JanAlignmentAlgorithm::weakLimit
private

normalized eigen value below which the (CS) eigen vectors are considered as weak

Definition at line 70 of file JanAlignmentAlgorithm.h.

Referenced by JanAlignmentAlgorithm(), and solve().