CSCChamberFitter Class Reference

#include <CSCChamberFitter.h>

Public Member Functions
	CSCChamberFitter (const edm::ParameterSet &iConfig, std::vector< CSCPairResidualsConstraint * > &residualsConstraints)
bool	fit (std::vector< CSCAlignmentCorrections * > &corrections) const
void	radiusCorrection (AlignableNavigator *alignableNavigator, AlignmentParameterStore *alignmentParameterStore, bool combineME11) const
virtual	~CSCChamberFitter ()

Protected Member Functions
long	alignableId (std::string alignable) const
double	chi2 (const AlgebraicVector &A, double lambda) const
double	hessian (int k, int l, double lambda) const
int	index (std::string alignable) const
bool	isFrame (int i) const
double	lhsVector (int k) const
void	walk (std::map< int, bool > &touched, int alignable) const

Protected Attributes
std::vector< std::string >	m_alignables
std::vector< CSCPairConstraint * >	m_constraints
int	m_fixed
std::vector< int >	m_frames
std::string	m_name

Detailed Description

Date:

2010/05/06 05:57:33

Revision:

1.2

Author

J. Pivarski - Texas A&M University pivar.nosp@m.ski@.nosp@m.physi.nosp@m.cs.t.nosp@m.amu.e.nosp@m.du

Definition at line 24 of file CSCChamberFitter.h.

Constructor & Destructor Documentation

CSCChamberFitter::CSCChamberFitter	(	const edm::ParameterSet &	iConfig,
		std::vector< CSCPairResidualsConstraint * > &	residualsConstraints
	)

Definition at line 13 of file CSCChamberFitter.cc.

References alignableId(), CSCDetId::chamber(), HLT_25ns14e33_v1_cff::constraint, createBeamHaloJobs::constraints, CSCDetId::endcap(), relativeConstraints::error, edm::hlt::Exception, edm::ParameterSet::getParameter(), i, index(), j, m_alignables, m_constraints, m_fixed, m_frames, m_name, CSCDetId::ring(), CSCDetId::station(), AlCaHLTBitMon_QueryRunRegistry::string, relativeConstraints::value, and walk().

                                                                                                                               {
  m_name = iConfig.getParameter<std::string>("name");
  m_alignables = iConfig.getParameter<std::vector<std::string> >("alignables");
  if (m_alignables.size() == 0) {
    throw cms::Exception("BadConfig") << "Fitter " << m_name << " has no alignables!" << std::endl;
  }

  int i = 0;
  for (std::vector<std::string>::const_iterator alignable = m_alignables.begin();  alignable != m_alignables.end();  ++alignable) {
    if (alignableId(*alignable) == -1) m_frames.push_back(i);
    i++;
  }

  m_fixed = -1;
  std::string fixed = iConfig.getParameter<std::string>("fixed");
  if (fixed != "") {
    int i = 0;
    for (std::vector<std::string>::const_iterator alignable = m_alignables.begin();  alignable != m_alignables.end();  ++alignable) {
      if (fixed == *alignable) {
    m_fixed = i;
      }
      i++;
    }
    if (m_fixed == -1) throw cms::Exception("BadConfig") << "Cannot fix unrecognized alignable " << fixed << std::endl;
  }

  int numConstraints = 0;
  std::vector<edm::ParameterSet> constraints = iConfig.getParameter<std::vector<edm::ParameterSet> >("constraints");
  for (std::vector<edm::ParameterSet>::const_iterator constraint = constraints.begin();  constraint != constraints.end();  ++constraint) {
    int i = index(constraint->getParameter<std::string>("i"));
    int j = index(constraint->getParameter<std::string>("j"));
    double value = constraint->getParameter<double>("value");
    double error = constraint->getParameter<double>("error");
    
    if (i < 0) throw cms::Exception("BadConfig") << "Unrecognized alignable " << constraint->getParameter<std::string>("i") << " in constraint " << numConstraints << " of fitter " << m_name << std::endl;
    if (j < 0) throw cms::Exception("BadConfig") << "Unrecognized alignable " << constraint->getParameter<std::string>("j") << " in constraint " << numConstraints << " of fitter " << m_name << std::endl;
    if (error <= 0.) throw cms::Exception("BadConfig") << "Non-positive uncertainty in constraint " << numConstraints << " of fitter " << m_name << std::endl;
    if (i == j) throw cms::Exception("BadConfig") << "Self-connection from " << constraint->getParameter<std::string>("i") << " to " << constraint->getParameter<std::string>("j") << " is not allowed in constraint " << numConstraints << " of fitter " << m_name << std::endl;

    m_constraints.push_back(new CSCPairConstraint(i, j, value, error));
    numConstraints++;
  }

  // insert CSCPairResidualsConstraints
  for (unsigned int i = 0;  i < m_alignables.size();  i++) {
    std::string alignable_i = m_alignables[i];
    long id_i = alignableId(alignable_i);
    if (id_i != -1) {
      CSCDetId cscid_i(id_i);

      for (unsigned int j = 0;  j < m_alignables.size();  j++) {
    std::string alignable_j = m_alignables[j];
    long id_j = alignableId(alignable_j);
    if (i != j  &&  id_j != -1) {
      CSCDetId cscid_j(id_j);

      if (!(cscid_i.station() == 1  &&  cscid_i.ring() == 3  &&  cscid_j.station() == 1  &&  cscid_j.ring() == 3)) {

         int next_chamber = cscid_i.chamber() + 1;
         if (cscid_i.station() > 1  &&  cscid_i.ring() == 1  &&  next_chamber == 19) next_chamber = 1;
         else if (!(cscid_i.station() > 1  &&  cscid_i.ring() == 1)  &&  next_chamber == 37) next_chamber = 1;
         if (cscid_i.endcap() == cscid_j.endcap()  &&  cscid_i.station() == cscid_j.station()  &&  cscid_i.ring() == cscid_j.ring()  &&  next_chamber == cscid_j.chamber()) {
        
        CSCPairResidualsConstraint *residualsConstraint = new CSCPairResidualsConstraint(residualsConstraints.size(), i, j, cscid_i, cscid_j);
        m_constraints.push_back(residualsConstraint);
        residualsConstraints.push_back(residualsConstraint);
        numConstraints++;       
         }
      }
    }
      }
    }
  }

  std::map<int,bool> touched;
  for (unsigned int i = 0;  i < m_alignables.size();  i++) touched[i] = false;
  walk(touched, 0);
  for (unsigned int i = 0;  i < m_alignables.size();  i++) {
    if (!touched[i]) throw cms::Exception("BadConfig") << "Fitter " << m_name << " is not a connected graph (no way to get to " << m_alignables[i] << " from " << m_alignables[0] << ", for instance)" << std::endl;
  }
}

virtual CSCChamberFitter::~CSCChamberFitter ( )

inlinevirtual

Definition at line 27 of file CSCChamberFitter.h.

{};

Member Function Documentation

long CSCChamberFitter::alignableId ( std::string  alignable ) const

protected

Definition at line 115 of file CSCChamberFitter.cc.

References CSCDetId, Reference_intrackfit_cff::endcap, relativeConstraints::ring, and relativeConstraints::station.

Referenced by CSCChamberFitter(), and fit().

                                                            {
  if (alignable.size() != 9) return -1;

  if (alignable[0] == 'M'  &&  alignable[1] == 'E') {
    int endcap = -1;
    if (alignable[2] == '+') endcap = 1;
    else if (alignable[2] == '-') endcap = 2;

    if (endcap != -1) {
      int station = -1;
      if (alignable[3] == '1') station = 1;
      else if (alignable[3] == '2') station = 2;
      else if (alignable[3] == '3') station = 3;
      else if (alignable[3] == '4') station = 4;

      if (alignable[4] == '/'  &&  station != -1) {
    int ring = -1;
    if (alignable[5] == '1') ring = 1;
    else if (alignable[5] == '2') ring = 2;
    else if (alignable[5] == '3') ring = 3;
    else if (alignable[5] == '4') ring = 4;
    if (station > 1  &&  ring > 2) return -1;

    if (alignable[6] == '/'  &&  ring != -1) {
      int chamber = -1;
      if      (alignable[7] == '0'  &&  alignable[8] == '1') chamber = 1;
      else if (alignable[7] == '0'  &&  alignable[8] == '2') chamber = 2;
      else if (alignable[7] == '0'  &&  alignable[8] == '3') chamber = 3;
      else if (alignable[7] == '0'  &&  alignable[8] == '4') chamber = 4;
      else if (alignable[7] == '0'  &&  alignable[8] == '5') chamber = 5;
      else if (alignable[7] == '0'  &&  alignable[8] == '6') chamber = 6;
      else if (alignable[7] == '0'  &&  alignable[8] == '7') chamber = 7;
      else if (alignable[7] == '0'  &&  alignable[8] == '8') chamber = 8;
      else if (alignable[7] == '0'  &&  alignable[8] == '9') chamber = 9;
      else if (alignable[7] == '1'  &&  alignable[8] == '0') chamber = 10;
      else if (alignable[7] == '1'  &&  alignable[8] == '1') chamber = 11;
      else if (alignable[7] == '1'  &&  alignable[8] == '2') chamber = 12;
      else if (alignable[7] == '1'  &&  alignable[8] == '3') chamber = 13;
      else if (alignable[7] == '1'  &&  alignable[8] == '4') chamber = 14;
      else if (alignable[7] == '1'  &&  alignable[8] == '5') chamber = 15;
      else if (alignable[7] == '1'  &&  alignable[8] == '6') chamber = 16;
      else if (alignable[7] == '1'  &&  alignable[8] == '7') chamber = 17;
      else if (alignable[7] == '1'  &&  alignable[8] == '8') chamber = 18;
      else if (alignable[7] == '1'  &&  alignable[8] == '9') chamber = 19;
      else if (alignable[7] == '2'  &&  alignable[8] == '0') chamber = 20;
      else if (alignable[7] == '2'  &&  alignable[8] == '1') chamber = 21;
      else if (alignable[7] == '2'  &&  alignable[8] == '2') chamber = 22;
      else if (alignable[7] == '2'  &&  alignable[8] == '3') chamber = 23;
      else if (alignable[7] == '2'  &&  alignable[8] == '4') chamber = 24;
      else if (alignable[7] == '2'  &&  alignable[8] == '5') chamber = 25;
      else if (alignable[7] == '2'  &&  alignable[8] == '6') chamber = 26;
      else if (alignable[7] == '2'  &&  alignable[8] == '7') chamber = 27;
      else if (alignable[7] == '2'  &&  alignable[8] == '8') chamber = 28;
      else if (alignable[7] == '2'  &&  alignable[8] == '9') chamber = 29;
      else if (alignable[7] == '3'  &&  alignable[8] == '0') chamber = 30;
      else if (alignable[7] == '3'  &&  alignable[8] == '1') chamber = 31;
      else if (alignable[7] == '3'  &&  alignable[8] == '2') chamber = 32;
      else if (alignable[7] == '3'  &&  alignable[8] == '3') chamber = 33;
      else if (alignable[7] == '3'  &&  alignable[8] == '4') chamber = 34;
      else if (alignable[7] == '3'  &&  alignable[8] == '5') chamber = 35;
      else if (alignable[7] == '3'  &&  alignable[8] == '6') chamber = 36;

      if (station > 1  &&  ring == 1  &&  chamber > 18) return -1;

      if (chamber != -1) {
        return CSCDetId(endcap, station, ring, chamber, 0).rawId();
      }
    }
      }
    }
  }

  return -1;
}

double CSCChamberFitter::chi2 ( const AlgebraicVector &  A, double  lambda  ) const

protected

Definition at line 197 of file CSCChamberFitter.cc.

References HLT_25ns14e33_v1_cff::constraint, i, isFrame(), m_alignables, m_constraints, m_fixed, funct::pow(), and alignCSCRings::s.

Referenced by fit().

                                                                           {
  double sumFixed = 0.;

  if (m_fixed == -1) {
    for (unsigned int i = 0;  i < m_alignables.size();  i++) {
      if (!isFrame(i)) {
    sumFixed += A[i];
      }
    }
  }
  else {
    sumFixed = A[m_fixed];
  }

  double s = lambda * sumFixed*sumFixed;
  for (std::vector<CSCPairConstraint*>::const_iterator constraint = m_constraints.begin();  constraint != m_constraints.end();  ++constraint) {
     if ((*constraint)->valid()) {
    s += pow((*constraint)->value() - A[(*constraint)->i()] + A[(*constraint)->j()], 2)/(*constraint)->error()/(*constraint)->error();
     }
  }
  return s;
}

bool CSCChamberFitter::fit

(

std::vector< CSCAlignmentCorrections * > &

corrections

)

const

Definition at line 254 of file CSCChamberFitter.cc.

References HLT_25ns14e33_v1_cff::A, alignableId(), chi2(), HLT_25ns14e33_v1_cff::constraint, CSCDetId, hessian(), i, infinity, CSCAlignmentCorrections::insertCorrection(), CSCAlignmentCorrections::insertMode(), CSCAlignmentCorrections::insertResidual(), isFrame(), j, relval_steps::k, prof2calltree::l, lhsVector(), m_alignables, m_constraints, m_name, and mathSSE::sqrt().

                                                                                 {
  double lambda = 1./infinity/infinity;

  AlgebraicVector A(m_alignables.size());
  AlgebraicVector V(m_alignables.size());
  AlgebraicMatrix M(m_alignables.size(), m_alignables.size());

  for (unsigned int k = 0;  k < m_alignables.size();  k++) {
    A[k] = 0.;
    V[k] = lhsVector(k);
    for (unsigned int l = 0;  l < m_alignables.size();  l++) {
      M[k][l] = hessian(k, l, lambda);
    }
  }

  double oldchi2 = chi2(A, lambda);

  int ierr;
  M.invert(ierr);
  if (ierr != 0) {
    edm::LogError("CSCOverlapsAlignmentAlgorithm") << "Matrix inversion failed for fitter " << m_name << " matrix is " << M << std::endl;
    return false;
  }

  A = M * V;  // that's the alignment step

  CSCAlignmentCorrections *correction = new CSCAlignmentCorrections(m_name, oldchi2, chi2(A, lambda));

  for (unsigned int i = 0;  i < m_alignables.size();  i++) {
    if (!isFrame(i)) {
      correction->insertCorrection(m_alignables[i], CSCDetId(alignableId(m_alignables[i])), A[i]);
    }
  }

  // we have to switch to a completely different linear algebrea
  // package because CLHEP doesn't compute
  // eigenvectors/diagonalization (?!?)
  TMatrixD tmatrix(m_alignables.size(), m_alignables.size());
  for (unsigned int i = 0;  i < m_alignables.size();  i++) {
    for (unsigned int j = 0;  j < m_alignables.size();  j++) {
      tmatrix[i][j] = M[i][j];
    }
  }
  TMatrixDEigen tmatrixdeigen(tmatrix); 
  TMatrixD basis = tmatrixdeigen.GetEigenVectors();
  TMatrixD invbasis = tmatrixdeigen.GetEigenVectors();
  invbasis.Invert();
  TMatrixD diagonalized = invbasis * (tmatrix * basis);

  for (unsigned int i = 0;  i < m_alignables.size();  i++) {
    std::vector<double> coefficient;
    std::vector<std::string> modename;
    std::vector<long> modeid;
    for (unsigned int j = 0;  j < m_alignables.size();  j++) {
      coefficient.push_back(invbasis[i][j]);
      modename.push_back(m_alignables[j]);
      modeid.push_back(alignableId(m_alignables[j]));
    }

    correction->insertMode(coefficient, modename, modeid, sqrt(2. * fabs(diagonalized[i][i])) * (diagonalized[i][i] >= 0. ? 1. : -1.));
  }

  for (std::vector<CSCPairConstraint*>::const_iterator constraint = m_constraints.begin();  constraint != m_constraints.end();  ++constraint) {
     if ((*constraint)->valid()) {
    double residual = (*constraint)->value() - A[(*constraint)->i()] + A[(*constraint)->j()];
    correction->insertResidual(m_alignables[(*constraint)->i()], m_alignables[(*constraint)->j()], (*constraint)->value(), (*constraint)->error(), residual, residual/(*constraint)->error());
     }
  }

  corrections.push_back(correction);
  return true;
}

double CSCChamberFitter::hessian ( int  k, int  l, double  lambda  ) const

protected

Definition at line 232 of file CSCChamberFitter.cc.

References HLT_25ns14e33_v1_cff::constraint, ztail::d, infinity, isFrame(), relval_steps::k, prof2calltree::l, m_constraints, m_fixed, and alignCSCRings::s.

Referenced by fit().

                                                                  {
  double s = 0.;

  if (m_fixed == -1) {
    if (!isFrame(k)  &&  !isFrame(l)) s += 2.*lambda;
  }
  else {
    if (k == l  &&  l == m_fixed) s += 2.*lambda;
  }

  for (std::vector<CSCPairConstraint*>::const_iterator constraint = m_constraints.begin();  constraint != m_constraints.end();  ++constraint) {
     double d = 2./infinity/infinity;
     if ((*constraint)->valid()) {
    d = 2./(*constraint)->error()/(*constraint)->error();
     }

     if (k == l  &&  ((*constraint)->i() == k  ||  (*constraint)->j() == k)) s += d;
     if (((*constraint)->i() == k  &&  (*constraint)->j() == l)  ||  ((*constraint)->j() == k  &&  (*constraint)->i() == l)) s -= d;
  }
  return s;
}

int CSCChamberFitter::index ( std::string  alignable ) const

protected

Definition at line 95 of file CSCChamberFitter.cc.

References a, i, and m_alignables.

Referenced by BeautifulSoup.PageElement::_invert(), and CSCChamberFitter().

                                                     {
  int i = 0;
  for (std::vector<std::string>::const_iterator a = m_alignables.begin();  a != m_alignables.end();  ++a) {
    if (*a == alignable) return i;
    i++;
  }
  return -1;
}

bool CSCChamberFitter::isFrame ( int  i ) const

protected

Definition at line 190 of file CSCChamberFitter.cc.

References m_frames.

Referenced by chi2(), fit(), and hessian().

                                          {
  for (std::vector<int>::const_iterator frame = m_frames.begin();  frame != m_frames.end();  ++frame) {
    if (i == *frame) return true;
  }
  return false;
}

double CSCChamberFitter::lhsVector ( int  k ) const

protected

Definition at line 220 of file CSCChamberFitter.cc.

References HLT_25ns14e33_v1_cff::constraint, ztail::d, relval_steps::k, m_constraints, and alignCSCRings::s.

Referenced by fit().

                                              {
  double s = 0.;
  for (std::vector<CSCPairConstraint*>::const_iterator constraint = m_constraints.begin();  constraint != m_constraints.end();  ++constraint) {
     if ((*constraint)->valid()) {
    double d = 2.*(*constraint)->value()/(*constraint)->error()/(*constraint)->error();
    if ((*constraint)->i() == k) s += d;
    if ((*constraint)->j() == k) s -= d;
     }
  }
  return s;
}

void CSCChamberFitter::radiusCorrection	(	AlignableNavigator *	alignableNavigator,
		AlignmentParameterStore *	alignmentParameterStore,
		bool	combineME11
	)		const

Definition at line 328 of file CSCChamberFitter.cc.

References AlignableDetOrUnitPtr::alignable(), AlignableNavigator::alignableFromDetId(), Alignable::alignmentParameters(), AlignmentParameterStore::applyParameters(), CSCDetId::chamber(), AlignmentParameters::cloneFromSelected(), HLT_25ns14e33_v1_cff::constraint, alignCSCRings::e, CSCDetId::endcap(), CSCPairResidualsConstraint::id_i(), m_constraints, M_PI, NULL, applyRadialCorrections::radial_correction, CSCPairResidualsConstraint::radius(), CSCDetId::ring(), Alignable::setAlignmentParameters(), AlignmentParameters::setValid(), CSCDetId::station(), CSCPairResidualsConstraint::valid(), and CSCPairResidualsConstraint::value().

                                                                                                                                                        {
   double sum_phipos_residuals = 0.;
   double num_valid = 0.;
   double sum_radius = 0.;
   double num_total = 0.;
   for (std::vector<CSCPairConstraint*>::const_iterator constraint = m_constraints.begin();  constraint != m_constraints.end();  ++constraint) {
      CSCPairResidualsConstraint *residualsConstraint = dynamic_cast<CSCPairResidualsConstraint*>(*constraint);
      if (residualsConstraint != NULL) {

     if (residualsConstraint->valid()) {
        sum_phipos_residuals += residualsConstraint->value();
        num_valid += 1.;
     }

     sum_radius += residualsConstraint->radius(true);
     num_total += 1.;
      }
   }
   if (num_valid == 0.  ||  num_total == 0.) return;
   double average_phi_residual = sum_phipos_residuals / num_valid;
   double average_radius = sum_radius / num_total;

   double radial_correction = average_phi_residual * average_radius * num_total / (2.*M_PI);

   for (std::vector<CSCPairConstraint*>::const_iterator constraint = m_constraints.begin();  constraint != m_constraints.end();  ++constraint) {
      CSCPairResidualsConstraint *residualsConstraint = dynamic_cast<CSCPairResidualsConstraint*>(*constraint);
      if (residualsConstraint != NULL) {

     const DetId id(residualsConstraint->id_i());
     Alignable *alignable = alignableNavigator->alignableFromDetId(id).alignable();
     Alignable *also = NULL;
     if (combineME11  &&  residualsConstraint->id_i().station() == 1  &&  residualsConstraint->id_i().ring() == 1) {
        CSCDetId alsoid(residualsConstraint->id_i().endcap(), 1, 4, residualsConstraint->id_i().chamber(), 0);
        const DetId alsoid2(alsoid);
        also = alignableNavigator->alignableFromDetId(alsoid2).alignable();
     }
      
     AlgebraicVector params(6);
     AlgebraicSymMatrix cov(6);

     params[1] = radial_correction;
     cov[1][1] = 1e-6;

     AlignmentParameters *parnew = alignable->alignmentParameters()->cloneFromSelected(params, cov);
     alignable->setAlignmentParameters(parnew);
     alignmentParameterStore->applyParameters(alignable);
     alignable->alignmentParameters()->setValid(true);
     if (also != NULL) {
        AlignmentParameters *parnew2 = also->alignmentParameters()->cloneFromSelected(params, cov);
        also->setAlignmentParameters(parnew2);
        alignmentParameterStore->applyParameters(also);
        also->alignmentParameters()->setValid(true);
     }
      }
   }
}

void CSCChamberFitter::walk ( std::map< int, bool > &  touched, int  alignable  ) const

protected

Definition at line 104 of file CSCChamberFitter.cc.

References HLT_25ns14e33_v1_cff::constraint, and m_constraints.

Referenced by CSCChamberFitter().

                                                                          {
  touched[alignable] = true;

  for (std::vector<CSCPairConstraint*>::const_iterator constraint = m_constraints.begin();  constraint != m_constraints.end();  ++constraint) {
    if (alignable == (*constraint)->i()  ||  alignable == (*constraint)->j()) {
      if (!touched[(*constraint)->i()]) walk(touched, (*constraint)->i());
      if (!touched[(*constraint)->j()]) walk(touched, (*constraint)->j());
    }
  }
}

Member Data Documentation

std::vector<std::string> CSCChamberFitter::m_alignables

protected

Definition at line 42 of file CSCChamberFitter.h.

Referenced by chi2(), CSCChamberFitter(), fit(), and index().

std::vector<CSCPairConstraint*> CSCChamberFitter::m_constraints

protected

Definition at line 45 of file CSCChamberFitter.h.

Referenced by chi2(), CSCChamberFitter(), fit(), hessian(), lhsVector(), radiusCorrection(), and walk().

int CSCChamberFitter::m_fixed

protected

Definition at line 44 of file CSCChamberFitter.h.

Referenced by chi2(), CSCChamberFitter(), and hessian().

std::vector<int> CSCChamberFitter::m_frames

protected

Definition at line 43 of file CSCChamberFitter.h.

Referenced by CSCChamberFitter(), and isFrame().

std::string CSCChamberFitter::m_name

protected

Definition at line 41 of file CSCChamberFitter.h.

Referenced by CSCChamberFitter(), and fit().