Skip to content

Commit

Permalink
Improved verbosity control.
Browse files Browse the repository at this point in the history
  • Loading branch information
@jan-kaspar
jan-kaspar committed Aug 21, 2020
1 parent dbadb61 commit a33b681
Show file tree
Hide file tree
Showing 10 changed files with 201 additions and 195 deletions.
3 changes: 0 additions & 3 deletions Alignment/PPSTrackBased/interface/StraightTrackAlignment.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,9 +65,6 @@ class StraightTrackAlignment
/// verbosity level
unsigned int verbosity;

/// verbosity level while factorization
unsigned int factorizationVerbosity;

/// list of RPs for which the alignment parameters shall be optimized
std::vector<unsigned int> rpIds;

Expand Down
3 changes: 0 additions & 3 deletions Alignment/PPSTrackBased/plugins/PPSFastLocalSimulation.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -315,15 +315,12 @@ void PPSFastLocalSimulation::GenerateTrack(unsigned int idx, CLHEP::HepRandomEng
Ai = A.Invert();
TVectorD P(3);
P = Ai * B;
//printf(" de z = %f, p1 = %f, p2 = %f\n", P(0), P(1), P(2));

double de_z = P(0);
CTPPSGeometry::Vector h_glo(ax * de_z + bx, ay * de_z + by, de_z + z0);
//printf(" h_glo: x = %f, y = %f, z = %f\n", h_glo.x(), h_glo.y(), h_glo.z());

// hit in local coordinates
CTPPSGeometry::Vector h_loc = geometry->globalToLocal(detId, h_glo);
//printf(" h_loc: c1 = %f, c2 = %f, c3 = %f\n", h_loc.x(), h_loc.y(), h_loc.z());

// strips
if (detId.subdetId() == CTPPSDetId::sdTrackingStrip)
Expand Down
3 changes: 1 addition & 2 deletions Alignment/PPSTrackBased/python/ppsStraightTrackAligner_cfi.py
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,7 @@
import FWCore.ParameterSet.Config as cms

ppsStraightTrackAligner = cms.EDAnalyzer("PPSStraightTrackAligner",
verbosity = cms.untracked.uint32(0),
factorizationVerbosity = cms.untracked.uint32(0),
verbosity = cms.untracked.uint32(1),

# ---------- input and event selection ----------

Expand Down
1 change: 0 additions & 1 deletion Alignment/PPSTrackBased/src/IdealResult.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,6 @@ vector<SingularMode> IdealResult::analyze()
unsigned int IdealResult::solve(const std::vector<AlignmentConstraint> &constraints,
std::map<unsigned int, AlignmentResult> &results, TDirectory *dir)
{
printf(">> IdealResult::Solve\n\tvalues in mm and rad\n");
results.clear();

// determine dimension and offsets
Expand Down
252 changes: 132 additions & 120 deletions Alignment/PPSTrackBased/src/JanAlignmentAlgorithm.cc
View file
Original file line number Diff line number Diff line change
@@ -1,9 +1,6 @@
/****************************************************************************
*
* This is a part of TOTEM offline software.
* Authors:
* Jan Kašpar ([email protected])
*
****************************************************************************/

#include "FWCore/MessageLogger/interface/MessageLogger.h"
Expand Down Expand Up @@ -164,8 +161,6 @@ void JanAlignmentAlgorithm::feed(const HitCollection &selection, const LocalTrac
double hy = hay * hit->z + hby;
double R = m(j) - (hx*C + hy*S); // (standard) residual

//printf("%i\t%i\t%i\t%.3f\n", j, id, idx, measVec[j]);

if (buildDiagnosticPlots)
{
statistics[id].m_dist->Fill(m(j));
Expand Down Expand Up @@ -344,7 +339,7 @@ vector<SingularMode> JanAlignmentAlgorithm::analyze()
}

// analyze symmetricity
if (verbosity > 2)
if (verbosity >= 3)
{
double maxDiff = 0., maxElem = 0.;
for (unsigned int i = 0; i < dim; i++)
Expand Down Expand Up @@ -396,7 +391,7 @@ vector<SingularMode> JanAlignmentAlgorithm::analyze()
}
}

#if 0
#if DEBUG
// print singular vectors
if (singularModes.size() > 0)
{
Expand Down Expand Up @@ -430,7 +425,9 @@ vector<SingularMode> JanAlignmentAlgorithm::analyze()
unsigned int JanAlignmentAlgorithm::solve(const std::vector<AlignmentConstraint> &constraints,
map<unsigned int, AlignmentResult> &results, TDirectory *dir)
{
printf(">> JanAlignmentAlgorithm::Solve\n");
if (verbosity)
printf(">> JanAlignmentAlgorithm::Solve\n");

results.clear();

// build C matrix
Expand Down Expand Up @@ -492,114 +489,130 @@ unsigned int JanAlignmentAlgorithm::solve(const std::vector<AlignmentConstraint>
TMatrixD CS_eigVec = CS_eig.GetEigenVectors();

// check regularity of CS matrix
printf("\n* eigen values of CS and S matrices (events = %u)\n", events);
printf(" # CS norm. CS S norm. S\n");
if (verbosity >= 2)
{
printf("\n* eigen values of CS and S matrices (events = %u)\n", events);
printf(" # CS norm. CS S norm. S\n");
}

unsigned int singularModeCount = 0;
vector<unsigned int> weakModeIdx;
for (int i = 0; i < CS_eigVal.GetNrows(); i++)
{
double CS_nev = CS_eigVal[i]/events;
printf("%4i%+12.2E%+12.2E", i, CS_eigVal[i], CS_nev);
const double CS_nev = CS_eigVal[i]/events;

if (fabs(CS_nev) < singularLimit)
{
singularModeCount++;
printf(" (S)");
} else {
if (fabs(CS_nev) < weakLimit)
{
weakModeIdx.push_back(i);
printf(" (W)");
} else {
printf(" ");
}
}

if (i < S_eigVal.GetNrows())
if (verbosity >= 2)
{
double S_nev = S_eigVal[i]/events;
printf("%+12.2E%+12.2E", S_eigVal[i], S_nev);
if (fabs(S_nev) < singularLimit)
printf("%4i%+12.2E%+12.2E", i, CS_eigVal[i], CS_nev);
if (fabs(CS_nev) < singularLimit)
{
singularModeCount++;
printf(" (S)");
else
if (fabs(S_nev) < weakLimit)
} else {
if (fabs(CS_nev) < weakLimit)
{
weakModeIdx.push_back(i);
printf(" (W)");
}
} else {
printf(" ");
}
}

printf("\n");
if (i < S_eigVal.GetNrows())
{
double S_nev = S_eigVal[i]/events;
printf("%+12.2E%+12.2E", S_eigVal[i], S_nev);
if (fabs(S_nev) < singularLimit)
printf(" (S)");
else
if (fabs(S_nev) < weakLimit)
printf(" (W)");
}

printf("\n");
}
}

// print weak vectors
if (weakModeIdx.size() > 0)
if (verbosity >= 2)
{
unsigned int columns = 10;
unsigned int first = 0;

while (first < weakModeIdx.size())
// print weak vectors
if (weakModeIdx.size() > 0)
{
unsigned int last = first + columns;
if (last >= weakModeIdx.size())
last = weakModeIdx.size();

printf("\n* CS weak modes\n | ");
for (unsigned int i = first; i < last; i++)
printf("%+10.3E ", CS_eigVal[weakModeIdx[i]]);
printf("\n--- | ");

for (unsigned int i = first; i < last; i++)
printf("---------- ");
printf("\n");
unsigned int columns = 10;
unsigned int first = 0;

// determine maximum elements
vector<double> maxs;
for (unsigned int i = first; i < last; i++)
while (first < weakModeIdx.size())
{
double max = 0;
for (unsigned int j = 0; j < dim + constraints.size(); j++)
unsigned int last = first + columns;
if (last >= weakModeIdx.size())
last = weakModeIdx.size();

printf("\n* CS weak modes\n | ");
for (unsigned int i = first; i < last; i++)
printf("%+10.3E ", CS_eigVal[weakModeIdx[i]]);
printf("\n--- | ");

for (unsigned int i = first; i < last; i++)
printf("---------- ");
printf("\n");

// determine maximum elements
vector<double> maxs;
for (unsigned int i = first; i < last; i++)
{
double v = fabs(CS_eigVec(weakModeIdx[i], j));
if (v > max)
max = v;
double max = 0;
for (unsigned int j = 0; j < dim + constraints.size(); j++)
{
double v = fabs(CS_eigVec(weakModeIdx[i], j));
if (v > max)
max = v;
}
maxs.push_back(max);
}
maxs.push_back(max);
}

for (unsigned int j = 0; j < dim + constraints.size(); j++)
{
printf("%3u | ", j);
for (unsigned int i = first; i < last; i++)

for (unsigned int j = 0; j < dim + constraints.size(); j++)
{
double v = CS_eigVec(weakModeIdx[i], j);
if (fabs(v)/maxs[i-first] > 1E-3)
printf("%+10.3E ", v);
else
printf(" . ");
printf("%3u | ", j);
for (unsigned int i = first; i < last; i++)
{
double v = CS_eigVec(weakModeIdx[i], j);
if (fabs(v)/maxs[i-first] > 1E-3)
printf("%+10.3E ", v);
else
printf(" . ");
}
printf("\n");
}
printf("\n");
}

first = last;
}
} else
printf("\n* CS has no weak modes\n");
first = last;
}
} else
printf("\n* CS has no weak modes\n");
}

// check the regularity of C^T E
if (E.GetNcols() == C.GetNcols())
if (verbosity >= 2)
{
TMatrixD CTE(C, TMatrixD::kTransposeMult, E);
print(CTE, "* CTE matrix:");
const double &det = CTE.Determinant();
printf("\n* det(CTE) = %E, max(CTE) = %E, det(CTE)/max(CTE) = %E\n\tmax(C) = %E, max(E) = %E, det(CTE)/max(C)/max(E) = %E\n",
det, CTE.Max(), det/CTE.Max(), C.Max(), E.Max(), det/C.Max()/E.Max());
} else
printf(">> JanAlignmentAlgorithm::Solve > WARNING: C matrix has %u, while E matrix %u columns.\n", C.GetNcols(), E.GetNcols());
if (E.GetNcols() == C.GetNcols())
{
TMatrixD CTE(C, TMatrixD::kTransposeMult, E);
print(CTE, "* CTE matrix:");
const double &det = CTE.Determinant();
printf("\n* det(CTE) = %E, max(CTE) = %E, det(CTE)/max(CTE) = %E\n\tmax(C) = %E, max(E) = %E, det(CTE)/max(C)/max(E) = %E\n",
det, CTE.Max(), det/CTE.Max(), C.Max(), E.Max(), det/C.Max()/E.Max());
} else {
printf(">> JanAlignmentAlgorithm::Solve > WARNING: C matrix has %u, while E matrix %u columns.\n", C.GetNcols(), E.GetNcols());
}
}

// stop if CS is singular
if (singularModeCount > 0)
if (singularModeCount > 0 && stopOnSingularModes)
{
LogProblem("JanAlignmentAlgorithm") << "\n>> JanAlignmentAlgorithm::Solve > ERROR: There are "
<< singularModeCount << " singular modes in CS matrix.";
if (stopOnSingularModes)
return 1;
LogProblem("PPS") << "\n>> JanAlignmentAlgorithm::Solve > ERROR: There are " << singularModeCount << " singular modes in CS matrix. Stopping.";
return 1;
}

// build MV vector
Expand All @@ -625,25 +638,27 @@ unsigned int JanAlignmentAlgorithm::solve(const std::vector<AlignmentConstraint>
EM2 = CS2I * S0 * CS2I;

TMatrixD EMdiff(EM2 - EM);
//Print(EMdiff, "EMdiff");

double max1 = -1., max2 = -1., maxDiff = -1.;
for (int i = 0; i < EMdiff.GetNrows(); i++)
if (verbosity >= 3)
{
for (int j = 0; j < EMdiff.GetNcols(); j++)
double max1 = -1., max2 = -1., maxDiff = -1.;
for (int i = 0; i < EMdiff.GetNrows(); i++)
{
if (maxDiff < fabs(EMdiff(i, j)))
maxDiff = fabs(EMdiff(i, j));
for (int j = 0; j < EMdiff.GetNcols(); j++)
{
if (maxDiff < fabs(EMdiff(i, j)))
maxDiff = fabs(EMdiff(i, j));

if (max1 < fabs(EM(i, j)))
max1 = fabs(EM(i, j));

if (max2 < fabs(EM2(i, j)))
max2 = fabs(EM2(i, j));
if (max1 < fabs(EM(i, j)))
max1 = fabs(EM(i, j));

if (max2 < fabs(EM2(i, j)))
max2 = fabs(EM2(i, j));
}
}
}

printf("EM max = %E, EM2 max = %E, EM2 - EM max = %E\n", max1, max2, maxDiff);
printf("EM max = %E, EM2 max = %E, EM2 - EM max = %E\n", max1, max2, maxDiff);
}

// tests
TMatrixD &U = CS_eigVec;
Expand All @@ -652,38 +667,35 @@ unsigned int JanAlignmentAlgorithm::solve(const std::vector<AlignmentConstraint>
//CSEi = UT * CS * U;
//Print(CSEi, "CSEi");


TMatrixD EMEi(EM);
EMEi = UT * EM * U;
//Print(EMEi, "*EMEi");

double max = -1.;
for (int i = 0; i < EMEi.GetNrows(); i++)
if (verbosity >= 3)
{
for (int j = 0; j < EMEi.GetNcols(); j++)
double max = -1.;
for (int i = 0; i < EMEi.GetNrows(); i++)
{
if (max < EMEi(i, j))
max = EMEi(i, j);
for (int j = 0; j < EMEi.GetNcols(); j++)
{
if (max < EMEi(i, j))
max = EMEi(i, j);
}
}
}

printf("max = %E\n", max);
printf("max = %E\n", max);
}

/*
double th = max / 100;
for (int i = 0; i < EMEi.GetNrows(); i++)
for (int j = 0; j < EMEi.GetNcols(); j++)
if (th < EMEi(i, j))
printf("element %u, %u: %E\n", i, j, EMEi(i, j));
*/

// print lambda values
printf("\n* Lambda (from the contribution of singular modes to MV)\n");
for (unsigned int i = 0; i < constraints.size(); i++)
if (verbosity >= 3)
{
printf("\t%u (%25s)\t%+10.1E +- %10.1E\n", i, constraints[i].name.c_str(),
printf("\n* Lambda (from the contribution of singular modes to MV)\n");
for (unsigned int i = 0; i < constraints.size(); i++)
{
printf("\t%u (%25s)\t%+10.1E +- %10.1E\n", i, constraints[i].name.c_str(),
AL[dim+i]*1E3,
sqrt(EM[i+dim][i+dim])*1E3);
}
}

// fill results
Expand Down
Loading

0 comments on commit a33b681

Please sign in to comment.