13 :squareMode(squareMode_), countEvents(0),mineta(mineta_), maxeta(maxeta_), minphi(minphi_), maxphi(maxphi_)
33 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)
40 std::vector<float> totalSolution(
Nchannels,1.);
41 std::vector<float> iterSolution(
Nchannels,1.);
42 std::vector<std::vector<float> > myEventMatrix(eventMatrix);
43 std::vector<float> myEnergyVector(energyVector);
46 for (
int iter=1;iter<=nIter;iter++)
49 for (
unsigned int ireg=0; ireg<
regMinEta.size(); ireg++)
51 std::vector<float> regIterSolution, regEnergyVector;
52 std::vector<int> regVmaxCeta, regVmaxCphi;
53 std::vector<std::vector<float> > regEventMatrix;
59 for (
unsigned int ia=0; ia<VmaxCeta.size(); ia++)
65 regVmaxCeta.push_back(VmaxCeta[ia]);
66 regVmaxCphi.push_back(VmaxCphi[ia]);
68 std::vector<float> regEvent = myEventMatrix[ia];
69 float regEnergy = energyVector[ia];
70 for (
int i2=0; i2<
Nxtals; i2++)
72 int iFullReg = regionalHH.
indexSqr2Reg(i2, VmaxCeta[ia], VmaxCphi[ia]);
75 regEnergy -= regEvent[i2];
79 regEventMatrix.push_back(regEvent);
80 regEnergyVector.push_back(regEnergy);
87 regIterSolution = regionalHH.
iterate(regEventMatrix, regVmaxCeta, regVmaxCphi, regEnergyVector);
92 for (
unsigned int i1=0; i1<regIterSolution.size(); i1++)
94 int regFrame = regLength/2;
96 int currRegEta = i1 / currRegPhiRange +
regMinEta[ireg];
97 int currRegPhi = i1 % currRegPhiRange + regMinPhi[ireg];
100 if ( (currRegEta >= (regMinEta[ireg]+regFrame*(!(regMinEta[ireg]==
mineta))) ) &&
102 (currRegPhi >= (regMinPhi[ireg]+regFrame*(!(regMinPhi[ireg]==
minphi))) ) &&
106 iterSolution[newindex] = regIterSolution[i1];
111 if (iterSolution.empty())
return iterSolution;
114 for (
int ievent = 0; ievent<
Nevents; ievent++)
116 myEventMatrix[ievent] =
recalibrateEvent(myEventMatrix[ievent], VmaxCeta[ievent], VmaxCphi[ievent], iterSolution);
122 totalSolution[
i] *= iterSolution[
i];
127 return totalSolution;
132 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)
136 std::vector<float> totalSolution(
Nchannels,1.);
137 std::vector<float> iterSolution;
138 std::vector<std::vector<float> > myEventMatrix(eventMatrix);
139 std::vector<float> myEnergyVector(energyVector);
144 for (
int iter=1;iter<=nIter;iter++)
156 for (j=0; j<
Nxtals; j++) {sumOverEnergy += myEventMatrix[
i][j];}
157 sumOverEnergy /= myEnergyVector[
i];
158 scale += sumOverEnergy;
168 iterSolution =
iterate(myEventMatrix,VmaxCeta,VmaxCphi,myEnergyVector);
170 if (iterSolution.empty())
return iterSolution;
173 for (
int ievent = 0; ievent<
Nevents; ievent++)
175 myEventMatrix[ievent] =
recalibrateEvent(myEventMatrix[ievent], VmaxCeta[ievent], VmaxCphi[ievent], iterSolution);
181 totalSolution[
i] *= iterSolution[
i];
186 return totalSolution;
191 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)
193 std::vector<float> solution;
199 std::cout <<
"Householder::runIter(): more channels to calibrate than events available. " << std::endl;
202 std::cout <<
" ****************** ERROR *********************" << std::endl;
211 std::cout <<
"Householder::runIter(): matrix dimensions non-conformant. " << std::endl;
212 std::cout <<
" energyVectorOrig.size()=" << energyVectorOrig.size() << std::endl;
214 std::cout <<
" ****************** ERROR *********************" << std::endl;
230 if( !decomposeSuccess )
232 std::cout <<
"Householder::runIter(): Failed: Singular condition in decomposition."<< std::endl;
233 std::cout <<
"***************** PROBLEM in DECOMPOSITION *************************"<<std::endl;
238 float mydbleps = 2.22045e-16;
239 float etasqr = mydbleps*mydbleps;
273 normy0 += solution[
i] * solution[
i];
274 norme1 += e[
i] * e[
i];
283 if (norme1>(0.0625*normy0))
298 while (norme1>(etasqr*normy0))
324 if (norme1>(0.0625*norme0))
break;
348 float beta,sigma,alphak,eventMatrixkk;
426 if (eventMatrixkk < 0.)
433 beta = 1 / (sigma - eventMatrixkk * alphak);
492 for (i=Nchannels-2; i>=0; i--)
514 std::vector<float> newEventSquare(eventSquare);
521 newEventSquare[
i] *= recalibrateVector[iFull];
523 return newEventSquare;
533 if (curr_eta * maxCeta <= 0) {
if (maxCeta > 0) curr_eta--;
else curr_eta++; }
536 if (curr_phi < 1) curr_phi += 360;
537 if (curr_phi > 360) curr_phi -= 360;
545 if ( (!negPhiDirection && (curr_phi >=
minphi && curr_phi <=
maxphi)) ||
546 (negPhiDirection && !(curr_phi >=
minphi && curr_phi <=
maxphi)) )
548 iFullphi = curr_phi -
minphi;
549 if (iFullphi < 0) iFullphi += 360;
550 regionIndex = (curr_eta -
mineta) * (
maxphi - minphi + 1 + 360*negPhiDirection) + iFullphi;
558 std::vector< std::vector<float> > fullMatrix;
569 foo[iFull] = eventMatrix[
i][
k];
571 fullMatrix.push_back(foo);
582 int remEta =
Neta % regLength;
583 if (remEta > regLength/2) remEta -= regLength;
584 int numSubRegEta =
Neta / regLength + (remEta < 0)*1;
586 int addtoEta = remEta / numSubRegEta;
587 int addtomoreEta = remEta % numSubRegEta;
590 int remPhi =
Nphi % regLength;
591 if (remPhi > regLength/2) remPhi -= regLength;
592 int numSubRegPhi =
Nphi / regLength + (remPhi < 0)*1;
594 int addtoPhi = remPhi / numSubRegPhi;
595 int addtomorePhi = remPhi % numSubRegPhi;
598 int addedLengthEta = 0;
599 int addedLengthPhi = 0;
600 int startIndexEta =
mineta;
604 for (
int i=0;
i < numSubRegEta;
i++)
606 addedLengthEta = regLength + addtoEta + addtomoreEta/
abs(addtomoreEta)*(
i<
abs(addtomoreEta));
607 endIndexEta = startIndexEta + addedLengthEta - 1;
611 for (
int j=0; j < numSubRegPhi; j++)
613 addedLengthPhi = regLength + addtoPhi + addtomorePhi/
abs(addtomorePhi)*(j<
abs(addtomorePhi));
614 endIndexPhi = startIndexPhi + addedLengthPhi - 1;
616 regMinPhi.push_back(startIndexPhi - regFrame*(j!=0) );
617 regMaxPhi.push_back(endIndexPhi + regFrame*(j!=(numSubRegPhi-1)) );
618 regMinEta.push_back(startIndexEta - regFrame*(
i!=0) );
619 regMaxEta.push_back(endIndexEta + regFrame*(
i!=(numSubRegEta-1)) );
621 startIndexPhi = endIndexPhi + 1;
623 startIndexEta = endIndexEta + 1;
void solve(std::vector< float > &y)
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.
void makeRegions(const int ®Length)
std::vector< int > regMinPhi
std::vector< int > regMinEta
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< int > regMaxEta
int indexSqr2Reg(const int &sqrIndex, const int &maxCeta, const int &maxCphi)
Method to translate from square indices to region indices.
std::vector< std::vector< float > > eventMatrixProc
Abs< T >::type abs(const T &t)
std::vector< std::vector< float > > eventMatrixOrig
std::vector< float > energyVectorProc
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 ®Length=5)
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.
~HouseholderDecomposition()
Destructor.
HouseholderDecomposition(int squareMode_=5, int mineta_=1, int maxeta_=85, int minphi_=1, int maxphi_=20)
Default constructor.
std::vector< float > alpha
std::vector< int > regMaxPhi