HouseholderDecomposition Class Reference

#include <HouseholderDecomposition.h>

Public Member Functions
	HouseholderDecomposition (int squareMode_=5, int mineta_=1, int maxeta_=85, int minphi_=1, int maxphi_=20)
	Default constructor. More...
int	indexSqr2Reg (const int &sqrIndex, const int &maxCeta, const int &maxCphi)
	Method to translate from square indices to region indices. More...
std::vector< float >	iterate (const std::vector< std::vector< float > > &eventMatrix, const std::vector< int > &VmaxCeta, const std::vector< int > &VmaxCphi, const std::vector< float > &energyVector, const int &nIter, const bool &normalizeFlag=false)
std::vector< float >	iterate (const std::vector< std::vector< float > > &eventMatrix, const std::vector< int > &VmaxCeta, const std::vector< int > &VmaxCphi, const std::vector< float > &energyVectorOrig)
	Run the Householder Algorithm. Returns the vector of calibration coefficients. More...
std::vector< float >	recalibrateEvent (const std::vector< float > &eventSquare, const int &maxCeta, const int &maxCphi, const std::vector< float > &recalibrateVector)
	Recalibrate before next iteration: give previous solution vector as argument. More...
std::vector< float >	runRegional (const std::vector< std::vector< float > > &eventMatrix, const std::vector< int > &VmaxCeta, const std::vector< int > &VmaxCphi, const std::vector< float > &energyVector, const int &nIter, const int &regLength=5)
	~HouseholderDecomposition ()
	Destructor. More...

Private Member Functions
bool	decompose ()
void	makeRegions (const int &regLength)
void	solve (std::vector< float > &y)
std::vector< std::vector< float > >	unzipMatrix (const std::vector< std::vector< float > > &eventMatrix, const std::vector< int > &VmaxCeta, const std::vector< int > &VmaxCphi)
	Unzips the skimmed matrix into a full matrix. More...

Private Attributes
std::vector< float >	alpha
int	countEvents
std::vector< float >	energyVectorProc
std::vector< std::vector< float > >	eventMatrixOrig
std::vector< std::vector< float > >	eventMatrixProc
int	maxeta
int	maxphi
int	mineta
int	minphi
int	Nchannels
int	Neta
int	Nevents
int	Nphi
int	Nxtals
std::vector< int >	pivot
std::vector< int >	regMaxEta
std::vector< int >	regMaxPhi
std::vector< int >	regMinEta
std::vector< int >	regMinPhi
float	sigmaReplacement
int	squareMode

Detailed Description

Implementation of the QR decomposition of a matrix using Householder transformation

13.03.2007: R.Ofierzynski

• using a reduced matrix
• implements Regional Householder Algorithm

Author

Lorenzo Agostino, R.Ofierzynski, CERN

Definition at line 17 of file HouseholderDecomposition.h.

Constructor & Destructor Documentation

HouseholderDecomposition::HouseholderDecomposition	(	int	squareMode_ = 5,
		int	mineta_ = 1,
		int	maxeta_ = 85,
		int	minphi_ = 1,
		int	maxphi_ = 20
	)

Default constructor.

Definition at line 12 of file HouseholderDecomposition.cc.

References maxeta, maxphi, mineta, minphi, Nchannels, Neta, Nphi, Nxtals, sigmaReplacement, and squareMode.

  :squareMode(squareMode_), countEvents(0),mineta(mineta_), maxeta(maxeta_), minphi(minphi_), maxphi(maxphi_)
{
  Neta = maxeta - mineta + 1;
  if (mineta * maxeta < 0) Neta--; // there's no eta index = 0
  Nphi = maxphi - minphi + 1;
  if (Nphi <0) Nphi += 360;
  
  Nchannels = Neta * Nphi; // no. of channels, get it from edges of the region

  Nxtals = squareMode*squareMode; // no. of xtals in one event

  sigmaReplacement = 0.00001; // the sum of columns is replaced by this value in case it is zero (e.g. dead crystal)

}

HouseholderDecomposition::~HouseholderDecomposition

(

)

Destructor.

Definition at line 28 of file HouseholderDecomposition.cc.

{
}

Member Function Documentation

bool HouseholderDecomposition::decompose ( )

private

Make decomposition input: m=number of events, n=number of channels, qr=event matrix output: qr = transformed event matrix, alpha, pivot returns a boolean value, true if decomposition worked, false if it didn't

Definition at line 344 of file HouseholderDecomposition.cc.

References alpha, beta, eventMatrixProc, mps_fire::i, gen::k, Nchannels, Nevents, pivot, sigmaReplacement, mathSSE::sqrt(), and y.

Referenced by iterate().

{
  int i,j,jbar,k;
  float beta,sigma,alphak,eventMatrixkk;
  std::vector<float> y(Nchannels);
  std::vector<float> sum(Nchannels);

  //  std::cout << "Householder::decompose() started" << std::endl;
  
  for (j=0;j<Nchannels;j++) 
    {
      // jth column sum: squared sum for each crystal
      sum[j]=0.;
      for (i=0;i<Nevents;i++)
    sum[j]+=eventMatrixProc[i][j]*eventMatrixProc[i][j];
      
      // bookkeeping vector
      pivot[j] = j;
    }
  
  for (k=0;k<Nchannels;k++) 
    {
      // kth Householder transformation
      sigma = sum[k];
      jbar = k;
      
      // go through all following columns
      // find the largest sumSquared in the following columns
      for (j=k+1;j<Nchannels;j++) 
    {
      if (sum[j] > sigma) 
        {
          sigma = sum[j];
          jbar=j;
        }
    }

      if (jbar != k) 
    {
      // column interchange:
      //     interchange within: bookkeeping vector, squaredSum, eventMatrixProc 

      i = pivot[k];
      pivot[k]=pivot[jbar];
      pivot[jbar]=i;

      sum[jbar]=sum[k];
      sum[k]=sigma;
      
      for (i=0;i<Nevents;i++) 
        {
          sigma=eventMatrixProc[i][k];
          eventMatrixProc[i][k]=eventMatrixProc[i][jbar];
          eventMatrixProc[i][jbar]=sigma;
        }
    } // end column interchange

      // now store in sigma the squared sum of the readoutEnergies for this column(crystal)
      sigma=0.;
      for (i=k;i<Nevents;i++)
    {
      sigma+=eventMatrixProc[i][k]*eventMatrixProc[i][k];
    }

      // found a zero-column, bail out
      if (sigma == 0.) 
    {
//    std::cout << "Householder::decompose() failed" << std::endl;
//    return false;
      // rof 14.12.2006: workaround to avoid failure of algorithm because of dead crystals:
      sigma = sigmaReplacement;
      //      std::cout << "Householder::decompose - found zero column " << jbar << ", replacing sum of column elements by " << sigma << std::endl;
    }


      // the following paragraph acts only on the diagonal element:
      // if element=0, then calculate alpha & beta

      // take the diagonal element
      eventMatrixkk = eventMatrixProc[k][k];

      if (eventMatrixkk < 0.) 
    alpha[k] = sqrt(sigma);
      else
    alpha[k] = sqrt(sigma) * (-1.);

      alphak = alpha[k];

      beta = 1 / (sigma - eventMatrixkk * alphak);
      // replace it
      eventMatrixProc[k][k] = eventMatrixkk - alphak;

      for (j=k+1; j<Nchannels; j++) 
    {
      y[j] = 0.;

      for (i=k; i<Nevents; i++)
        {
          y[j] += eventMatrixProc[i][k] * eventMatrixProc[i][j];
        }

      y[j] *= beta;
    }

      for (j=k+1; j<Nchannels; j++) 
    {
      for (i=k; i<Nevents; i++) 
        {
          eventMatrixProc[i][j] -= eventMatrixProc[i][k] * y[j];
          sum[j] -= eventMatrixProc[k][j] * eventMatrixProc[k][j];
        }
    }
    } // end of kth householder transformation
  
  //  std::cout << "Householder::decompose() finished" << std::endl;
  
  return true;
}

int HouseholderDecomposition::indexSqr2Reg	(	const int &	sqrIndex,
		const int &	maxCeta,
		const int &	maxCphi
	)

Method to translate from square indices to region indices.

Definition at line 526 of file HouseholderDecomposition.cc.

References maxeta, maxphi, mineta, minphi, and squareMode.

Referenced by recalibrateEvent(), runRegional(), and unzipMatrix().

{
  int regionIndex;

  // get the current eta, phi indices
  int curr_eta = maxCeta - squareMode/2 + sqrIndex%squareMode;
  if (curr_eta * maxCeta <= 0) {if (maxCeta > 0) curr_eta--; else curr_eta++; }  // JUMP over 0

  int curr_phi = maxCphi - squareMode/2 + sqrIndex/squareMode;
  if (curr_phi < 1) curr_phi += 360;
  if (curr_phi > 360) curr_phi -= 360;

  bool negPhiDirection = (maxphi < minphi);
  int iFullphi;

  regionIndex = -1;

  if (curr_eta >= mineta && curr_eta <= maxeta)
    if ( (!negPhiDirection && (curr_phi >= minphi && curr_phi <= maxphi)) ||
     (negPhiDirection && !(curr_phi >= minphi && curr_phi <= maxphi))      ) 
      {
    iFullphi = curr_phi - minphi;
    if (iFullphi < 0) iFullphi += 360;
    regionIndex = (curr_eta - mineta) * (maxphi - minphi + 1 + 360*negPhiDirection) + iFullphi;
      }

  return regionIndex;
}

std::vector< float > HouseholderDecomposition::iterate	(	const std::vector< std::vector< float > > &	eventMatrix,
		const std::vector< int > &	VmaxCeta,
		const std::vector< int > &	VmaxCphi,
		const std::vector< float > &	energyVector,
		const int &	nIter,
		const bool &	normalizeFlag = false
	)

Run the Householder Algorithm several times (nIter). Returns the final vector of calibration coefficients. Comment from author: unless you do a new selection in between the iterations you don't need to run more than once; a second iteration on the same events does not improve the result in this case

Definition at line 131 of file HouseholderDecomposition.cc.

References mps_fire::i, Nchannels, Nevents, Nxtals, recalibrateEvent(), and Scenarios_cff::scale.

Referenced by ElectronCalibration::endJob(), ElectronCalibrationUniv::endJob(), and runRegional().

{
  Nevents = eventMatrix.size(); // Number of events to calibrate with

  std::vector<float> totalSolution(Nchannels,1.);
  std::vector<float> iterSolution;
  std::vector<std::vector<float> > myEventMatrix(eventMatrix);
  std::vector<float> myEnergyVector(energyVector);

  int i, j;

  // Iterate the correction
  for (int iter=1;iter<=nIter;iter++) 
    {

      // if normalization flag is set, normalize energies
      float sumOverEnergy;
      if (normalizeFlag)
    {
      float scale = 0.;
      
      for (i=0; i<Nevents; i++)
        {
          sumOverEnergy = 0.;
          for (j=0; j<Nxtals; j++) {sumOverEnergy += myEventMatrix[i][j];}
          sumOverEnergy /= myEnergyVector[i];
          scale += sumOverEnergy;
        }
      scale /= Nevents;
      
      for (i=0; i<Nevents; i++) {myEnergyVector[i] *= scale;}     
    } // end normalize energies



      // now the real work starts:
      iterSolution = iterate(myEventMatrix,VmaxCeta,VmaxCphi,myEnergyVector);

      if (iterSolution.empty()) return iterSolution;

      // re-calibrate eventMatrix with solution
      for (int ievent = 0; ievent<Nevents; ievent++)
    {
      myEventMatrix[ievent] = recalibrateEvent(myEventMatrix[ievent], VmaxCeta[ievent], VmaxCphi[ievent], iterSolution);
    }

      for (int i=0; i<Nchannels; i++) 
    {
      // save solution into theCalibVector
      totalSolution[i] *= iterSolution[i];
    }

    } // end iterate correction

  return totalSolution;
}

std::vector< float > HouseholderDecomposition::iterate	(	const std::vector< std::vector< float > > &	eventMatrix,
		const std::vector< int > &	VmaxCeta,
		const std::vector< int > &	VmaxCphi,
		const std::vector< float > &	energyVectorOrig
	)

Run the Householder Algorithm. Returns the vector of calibration coefficients.

Definition at line 190 of file HouseholderDecomposition.cc.

References alpha, gather_cfg::cout, decompose(), MillePedeFileConverter_cfg::e, energyVectorProc, eventMatrixOrig, eventMatrixProc, mps_fire::i, Nchannels, Nevents, pivot, solve(), and unzipMatrix().

{
  std::vector<float> solution; 

  Nevents=eventMatrix.size();      // Number of events to calibrate with

  if (Nchannels > Nevents)
    {
      std::cout << "Householder::runIter(): more channels to calibrate than events available. " << std::endl;
      std::cout << "  Nchannels=" << Nchannels << std::endl;
      std::cout << "  Nevents=" << Nevents << std::endl;
      std::cout << " ******************    ERROR   *********************" << std::endl;
      return solution; // empty vector
    }

  // input: eventMatrixOrig - the unzipped matrix
  eventMatrixOrig = unzipMatrix(eventMatrix,VmaxCeta,VmaxCphi);

 if (eventMatrixOrig.size() != energyVectorOrig.size())
    {
      std::cout << "Householder::runIter(): matrix dimensions non-conformant. " << std::endl;
      std::cout << "  energyVectorOrig.size()=" << energyVectorOrig.size() << std::endl;
      std::cout << "  eventMatrixOrig.size()=" << eventMatrixOrig.size() << std::endl;
      std::cout << " ******************    ERROR   *********************" << std::endl;
      return solution; // empty vector
    }


  int i,j;
  eventMatrixProc = eventMatrixOrig;
  energyVectorProc = energyVectorOrig; // copy energyVectorOrig vector
  std::vector<float> e(Nchannels);
  alpha.assign(Nchannels,0.);
  pivot.assign(Nchannels,0);
  

  //--------------------
  bool decomposeSuccess = decompose();

  if( !decomposeSuccess ) 
    {
      std::cout << "Householder::runIter(): Failed: Singular condition in decomposition."<< std::endl;
      std::cout << "***************** PROBLEM in DECOMPOSITION *************************"<<std::endl;
      return solution; // empty vector
    }

  /* DBL_EPSILON: Difference between 1.0 and the minimum float greater than 1.0 */
  float mydbleps = 2.22045e-16;  //DBL_EPSILON;
  float etasqr = mydbleps*mydbleps; 
  //  std::cout << "LOOK at DBL_EPSILON :" << mydbleps <<std::endl;


  //--------------------
  // apply transformations to rhs - find solution vector
  solution.assign(Nchannels,0.);
  solve(solution);


  // compute residual vector energyVectorProc
  for (i=0;i<Nevents;i++) 
    {
      energyVectorProc[i] = energyVectorOrig[i];
      for (j=0; j<Nchannels; j++)
    {
      energyVectorProc[i]-=eventMatrixOrig[i][j]*solution[j];
    }
    }


  //--------------------
  // compute first correction vector e
  solve(e);

  float normy0=0.;
  float norme1=0.;
  float norme0;
    



  for (i=0;i<Nchannels;i++) 
    {
      normy0 += solution[i] * solution[i];
      norme1 += e[i] * e[i];
    }
  
//  std::cout << "Householder::runIter(): applying first correction";
//  std::cout << " normy0 = " << normy0;
//  std::cout << " norme1 = " << norme1 << std::endl;

  // not attempt at obtaining the solution is made unless the norm of the first
  // correction  is significantly smaller than the norm of the initial solution
  if (norme1>(0.0625*normy0)) 
    {
      //      std::cout << "Householder::runIter(): first correction is too large. Failed." << std::endl;
    }

  // improve the solution
  for (i=0; i<Nchannels; i++)
    {
      solution[i]+=e[i];
    }

  //  std::cout << "Householder::runIter(): improving solution...." << std::endl;

  //--------------------
  // only continue iteration if the correction was significant
  while (norme1>(etasqr*normy0)) 
    {
      //      std::cout << "Householder::runIter(): norme1 = " << norme1 << std::endl;
    
      for (i=0; i<Nevents; i++) 
    {
      energyVectorProc[i] = energyVectorOrig[i];
      for (j=0; j<Nchannels; j++)
        {
          energyVectorProc[i] -= eventMatrixOrig[i][j] * solution[j];
        }
    }

    // compute next correction vector
    solve(e);

    norme0 = norme1;
    norme1 = 0.;

    for (i=0;i<Nchannels;i++)
      {
    norme1+=e[i]*e[i];
      }

    // terminate iteration if the norm of the new correction failed to decrease
    // significantly compared to the norm of the previous correction
    if (norme1>(0.0625*norme0)) break;

    // apply correction vector
    for (i=0;i<Nchannels;i++)
      {
    solution[i]+=e[i];
      }
  }

  //clean up
  eventMatrixOrig.clear();
  eventMatrixProc.clear();
  energyVectorProc.clear();
  alpha.clear();
  pivot.clear();

  
  return solution;
}

void HouseholderDecomposition::makeRegions ( const int &  regLength )

private

Determines the regions used for splitting of the full matrix and calibrating separately used by the public runRegional method

Definition at line 576 of file HouseholderDecomposition.cc.

References funct::abs(), mps_fire::i, mineta, minphi, Neta, Nphi, regMaxEta, regMaxPhi, regMinEta, regMinPhi, and squareMode.

Referenced by runRegional().

{
  //  int regFrame = regLength/2;
  int regFrame = squareMode/2;
  // first eta:
  int remEta = Neta % regLength;
  if (remEta > regLength/2) remEta -= regLength;
  int numSubRegEta = Neta / regLength + (remEta < 0)*1;

  int addtoEta = remEta / numSubRegEta;
  int addtomoreEta = remEta % numSubRegEta; // add "addtomore" number of times (addto+1), remaining times add just (addto)

  // then phi:
  int remPhi = Nphi % regLength;
  if (remPhi > regLength/2) remPhi -= regLength;
  int numSubRegPhi = Nphi / regLength + (remPhi < 0)*1;

  int addtoPhi = remPhi / numSubRegPhi;
  int addtomorePhi = remPhi % numSubRegPhi; // add "addtomore" number of times (addto+-1), remaining times add just (addto)

  // now put it all together
  int addedLengthEta = 0;
  int addedLengthPhi = 0;
  int startIndexEta = mineta;
  int startIndexPhi;
  int endIndexEta = 0;
  int endIndexPhi;
  for (int i=0; i < numSubRegEta; i++)
    {
      addedLengthEta = regLength + addtoEta + addtomoreEta/abs(addtomoreEta)*(i<abs(addtomoreEta));
      endIndexEta = startIndexEta + addedLengthEta - 1;

      startIndexPhi = minphi;
      endIndexPhi = 0;
      for (int j=0; j < numSubRegPhi; j++)
    {
      addedLengthPhi = regLength + addtoPhi + addtomorePhi/abs(addtomorePhi)*(j<abs(addtomorePhi));
      endIndexPhi = startIndexPhi + addedLengthPhi - 1;

      regMinPhi.push_back(startIndexPhi - regFrame*(j!=0) );
      regMaxPhi.push_back(endIndexPhi + regFrame*(j!=(numSubRegPhi-1)) );
      regMinEta.push_back(startIndexEta - regFrame*(i!=0) );
      regMaxEta.push_back(endIndexEta + regFrame*(i!=(numSubRegEta-1)) );

      startIndexPhi = endIndexPhi + 1;
    }
      startIndexEta = endIndexEta + 1;
    }

//  // print it all
//  std::cout << "Householder::makeRegions created the following regions for calibration:" << std::endl;
//  for (int i=0; i<regMinEta.size(); i++)
//    std::cout << "Region " << i << ": eta = " << regMinEta[i] << " to " << regMaxEta[i] << ", phi = " << regMinPhi[i] << " to " << regMaxPhi[i] << std::endl;

}

std::vector< float > HouseholderDecomposition::recalibrateEvent	(	const std::vector< float > &	eventSquare,
		const int &	maxCeta,
		const int &	maxCphi,
		const std::vector< float > &	recalibrateVector
	)

Recalibrate before next iteration: give previous solution vector as argument.

Definition at line 511 of file HouseholderDecomposition.cc.

References mps_fire::i, indexSqr2Reg(), and Nxtals.

Referenced by iterate(), and runRegional().

{
  std::vector<float> newEventSquare(eventSquare);
  int iFull;

  for (int i=0; i<Nxtals; i++) 
    {
      iFull = indexSqr2Reg(i, maxCeta, maxCphi);
      if (iFull >=0)
    newEventSquare[i] *= recalibrateVector[iFull];
    }
  return newEventSquare;
}

std::vector< float > HouseholderDecomposition::runRegional	(	const std::vector< std::vector< float > > &	eventMatrix,
		const std::vector< int > &	VmaxCeta,
		const std::vector< int > &	VmaxCphi,
		const std::vector< float > &	energyVector,
		const int &	nIter,
		const int &	regLength = 5
	)

Run Regional HouseholderAlgorithm (fast version), that splits matrix into regional matrices and inverts them separately. Returns the final vector of calibration coefficients. input: eventMatrix - the skimmed event matrix, VmaxCeta, VmaxCphi - vectors containing eta and phi indices of the maximum containment crystal for each event, energyVector - the energy vector, nIter - number of iterations to be performed, regLength - default length of the region (in eta- and phi-indices), regLength=5 recommended Comment from author: if you use the same events, 2 iterations are recommended; the second iteration gives corrections of the order of 0.0001

Definition at line 33 of file HouseholderDecomposition.cc.

References mps_fire::i, indexSqr2Reg(), iterate(), makeRegions(), maxeta, maxphi, mineta, minphi, Nchannels, Nevents, Nphi, Nxtals, recalibrateEvent(), regMaxEta, regMaxPhi, regMinEta, regMinPhi, and squareMode.

Referenced by ElectronCalibration::endJob(), and ElectronCalibrationUniv::endJob().

{
  // make regions
  makeRegions(regLength);

  Nevents = eventMatrix.size(); // Number of events to calibrate with

  std::vector<float> totalSolution(Nchannels,1.);
  std::vector<float> iterSolution(Nchannels,1.);
  std::vector<std::vector<float> > myEventMatrix(eventMatrix);

  // loop over nIter
  for (int iter=1;iter<=nIter;iter++) 
    {
      // loop over regions
      for (unsigned int ireg=0; ireg<regMinEta.size(); ireg++)
    {
      std::vector<float> regIterSolution, regEnergyVector;
      std::vector<int> regVmaxCeta, regVmaxCphi;
      std::vector<std::vector<float> > regEventMatrix;

      // initialize new instance with regional min,max indices
      HouseholderDecomposition regionalHH(squareMode,regMinEta[ireg],regMaxEta[ireg],regMinPhi[ireg],regMaxPhi[ireg]);

      // copy all events in region into new eventmatrix, energyvector, VmaxCeta, VmaxCphi
      for (unsigned int ia=0; ia<VmaxCeta.size(); ia++)
        {
          if ((VmaxCeta[ia] >= regMinEta[ireg]) && (VmaxCeta[ia] <= regMaxEta[ireg])
          && (VmaxCphi[ia] >= regMinPhi[ireg]) && (VmaxCphi[ia] <= regMaxPhi[ireg]))
        {
          // save event, calculate new eventmatrix(truncated) and energy
          regVmaxCeta.push_back(VmaxCeta[ia]);
          regVmaxCphi.push_back(VmaxCphi[ia]);

          std::vector<float> regEvent = myEventMatrix[ia];
          float regEnergy = energyVector[ia];
          for (int i2=0; i2<Nxtals; i2++)
            {
              int iFullReg = regionalHH.indexSqr2Reg(i2, VmaxCeta[ia], VmaxCphi[ia]);
              if (iFullReg <0) // crystal outside
            {
              regEnergy -= regEvent[i2];
              regEvent[i2] = 0.;
            }
            }
          regEventMatrix.push_back(regEvent);
          regEnergyVector.push_back(regEnergy);
        }
        }

      // calibrate
      //      std::cout << "HouseholderDecomposition::runRegional - Starting calibration of region " << ireg << ": eta " 
      //           << regMinEta[ireg] << " to " << regMaxEta[ireg] << ", phi " << regMinPhi[ireg] << " to " << regMaxPhi[ireg] << std::endl;
      regIterSolution = regionalHH.iterate(regEventMatrix, regVmaxCeta, regVmaxCphi, regEnergyVector);
      //      std::cout << "HouseholderDecomposition::runRegional - calibration of region finished. " << std::endl;

      // save solution into global iterSolution
      // don't forget to delete the ones that are on the border !
      for (unsigned int i1=0; i1<regIterSolution.size(); i1++)
        {
          int regFrame = regLength/2;
          int currRegPhiRange = regMaxPhi[ireg] - regMinPhi[ireg] + 1;
          int currRegEta = i1 / currRegPhiRange + regMinEta[ireg];
          int currRegPhi = i1 % currRegPhiRange + regMinPhi[ireg];
          int newindex = -100;
          // if crystal well inside:
          if ( (currRegEta >= (regMinEta[ireg]+regFrame*(!(regMinEta[ireg]==mineta))) ) &&
           (currRegEta <= (regMaxEta[ireg]-regFrame*(!(regMaxEta[ireg]==maxeta))) ) &&
           (currRegPhi >= (regMinPhi[ireg]+regFrame*(!(regMinPhi[ireg]==minphi))) ) &&
           (currRegPhi <= (regMaxPhi[ireg]-regFrame*(!(regMaxPhi[ireg]==maxphi))) ) )
        {
          newindex = (currRegEta-mineta)*Nphi + currRegPhi-minphi;
          iterSolution[newindex] = regIterSolution[i1];
        }
        }
    } // end loop over regions

      if (iterSolution.empty()) return iterSolution;
      
      // re-calibrate eventMatrix with solution
      for (int ievent = 0; ievent<Nevents; ievent++)
    {
      myEventMatrix[ievent] = recalibrateEvent(myEventMatrix[ievent], VmaxCeta[ievent], VmaxCphi[ievent], iterSolution);
    }
      
      // save solution into theCalibVector
      for (int i=0; i<Nchannels; i++) 
    {
      totalSolution[i] *= iterSolution[i];
    }
      
    } // end loop over nIter

  return totalSolution;

}

void HouseholderDecomposition::solve ( std::vector< float > &  y )

private

Apply transformations to rhs output: r = residual vector (energy vector), y = solution

Definition at line 464 of file HouseholderDecomposition.cc.

References alpha, energyVectorProc, eventMatrixProc, CustomPhysics_cfi::gamma, mps_fire::i, Nchannels, Nevents, pivot, and z.

Referenced by iterate().

{
  std::vector<float> z(Nchannels,0.);

  float gamma;
  int i,j;

  //  std::cout << "Householder::solve() begin" << std::endl;

  for (j=0; j<Nchannels; j++) 
    {
      // apply jth transformation to the right hand side
      gamma = 0.;
      for (i=j; i<Nevents; i++)
    {
      gamma += eventMatrixProc[i][j] * energyVectorProc[i];
    }
      gamma /= (alpha[j] * eventMatrixProc[j][j]);

      for (i=j; i<Nevents; i++)
    {
      energyVectorProc[i] += gamma * eventMatrixProc[i][j];
    }
    }

  z[Nchannels-1] = energyVectorProc[Nchannels-1] / alpha[Nchannels-1];

  for (i=Nchannels-2; i>=0; i--) 
    {
      z[i] = energyVectorProc[i];
      for (j=i+1; j<Nchannels; j++)
    {
      z[i] -= eventMatrixProc[i][j]*z[j];
    }
      z[i] /= alpha[i];
    }
  
  for (i=0; i<Nchannels; i++)
    {
      y[pivot[i]] = z[i];
    }
  
  //  std::cout << "Householder::solve() finished." << std::endl;
  
}

std::vector< std::vector< float > > HouseholderDecomposition::unzipMatrix ( const std::vector< std::vector< float > > &  eventMatrix, const std::vector< int > &  VmaxCeta, const std::vector< int > &  VmaxCphi  )

private

Unzips the skimmed matrix into a full matrix.

Definition at line 555 of file HouseholderDecomposition.cc.

References mps_fire::i, indexSqr2Reg(), gen::k, Nchannels, Nevents, and Nxtals.

Referenced by iterate().

{
  std::vector< std::vector<float> > fullMatrix;

  int iFull;

  for (int i=0; i<Nevents; i++)
    {
      std::vector<float> foo(Nchannels,0.);
      for (int k=0; k<Nxtals; k++)
    {
      iFull = indexSqr2Reg(k, VmaxCeta[i], VmaxCphi[i]);
      if (iFull >=0)
        foo[iFull] = eventMatrix[i][k];
    }
      fullMatrix.push_back(foo);
    }

  return fullMatrix;
}

Member Data Documentation

std::vector<float> HouseholderDecomposition::alpha

private

Definition at line 75 of file HouseholderDecomposition.h.

Referenced by decompose(), iterate(), and solve().

int HouseholderDecomposition::countEvents

private

Definition at line 69 of file HouseholderDecomposition.h.

std::vector<float> HouseholderDecomposition::energyVectorProc

private

Definition at line 74 of file HouseholderDecomposition.h.

Referenced by iterate(), and solve().

std::vector< std::vector<float> > HouseholderDecomposition::eventMatrixOrig

private

Definition at line 72 of file HouseholderDecomposition.h.

Referenced by iterate().

std::vector< std::vector<float> > HouseholderDecomposition::eventMatrixProc

private

Definition at line 73 of file HouseholderDecomposition.h.

Referenced by decompose(), iterate(), and solve().

int HouseholderDecomposition::maxeta

private

Definition at line 70 of file HouseholderDecomposition.h.

Referenced by HouseholderDecomposition(), indexSqr2Reg(), and runRegional().

int HouseholderDecomposition::maxphi

private

Definition at line 70 of file HouseholderDecomposition.h.

Referenced by HouseholderDecomposition(), indexSqr2Reg(), and runRegional().

int HouseholderDecomposition::mineta

private

Definition at line 70 of file HouseholderDecomposition.h.

Referenced by HouseholderDecomposition(), indexSqr2Reg(), makeRegions(), and runRegional().

int HouseholderDecomposition::minphi

private

Definition at line 70 of file HouseholderDecomposition.h.

Referenced by HouseholderDecomposition(), indexSqr2Reg(), makeRegions(), and runRegional().

int HouseholderDecomposition::Nchannels

private

Definition at line 71 of file HouseholderDecomposition.h.

Referenced by decompose(), HouseholderDecomposition(), iterate(), runRegional(), solve(), and unzipMatrix().

int HouseholderDecomposition::Neta

private

Definition at line 70 of file HouseholderDecomposition.h.

Referenced by HouseholderDecomposition(), and makeRegions().

int HouseholderDecomposition::Nevents

private

Definition at line 71 of file HouseholderDecomposition.h.

Referenced by decompose(), iterate(), runRegional(), solve(), and unzipMatrix().

int HouseholderDecomposition::Nphi

private

Definition at line 70 of file HouseholderDecomposition.h.

Referenced by HouseholderDecomposition(), makeRegions(), and runRegional().

int HouseholderDecomposition::Nxtals

private

Definition at line 71 of file HouseholderDecomposition.h.

Referenced by HouseholderDecomposition(), iterate(), recalibrateEvent(), runRegional(), and unzipMatrix().

std::vector<int> HouseholderDecomposition::pivot

private

Definition at line 76 of file HouseholderDecomposition.h.

Referenced by decompose(), iterate(), and solve().

std::vector<int> HouseholderDecomposition::regMaxEta

private

Definition at line 78 of file HouseholderDecomposition.h.

Referenced by makeRegions(), and runRegional().

std::vector<int> HouseholderDecomposition::regMaxPhi

private

Definition at line 78 of file HouseholderDecomposition.h.

Referenced by makeRegions(), and runRegional().

std::vector<int> HouseholderDecomposition::regMinEta

private

Definition at line 78 of file HouseholderDecomposition.h.

Referenced by makeRegions(), and runRegional().

std::vector<int> HouseholderDecomposition::regMinPhi

private

Definition at line 78 of file HouseholderDecomposition.h.

Referenced by makeRegions(), and runRegional().

float HouseholderDecomposition::sigmaReplacement

private

Definition at line 79 of file HouseholderDecomposition.h.

Referenced by decompose(), and HouseholderDecomposition().

int HouseholderDecomposition::squareMode

private

Definition at line 69 of file HouseholderDecomposition.h.

Referenced by HouseholderDecomposition(), indexSqr2Reg(), makeRegions(), and runRegional().