本文整理汇总了C++中TMatrixD::Invert方法的典型用法代码示例。如果您正苦于以下问题：C++ TMatrixD::Invert方法的具体用法？C++ TMatrixD::Invert怎么用？C++ TMatrixD::Invert使用的例子？那么, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类TMatrixD的用法示例。
在下文中一共展示了TMatrixD::Invert方法的7个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: calculateInvertedMatrixErrors
void calculateInvertedMatrixErrors(TMatrixD &T, TMatrixD &TErrPos, TMatrixD &TErrNeg,
TMatrixD &TinvErr){
// Calculate errors on the inverted matrix by the Monte Carlo
// method
Double_t det;
int nRow = T.GetNrows();
int nCol = T.GetNcols();
TMatrixD TinvSum(nRow,nCol);
TMatrixD TinvSumSquares(nRow,nCol);
// Reset Matrix where we will be accumulating RMS/sigma:
TinvSum = 0;
TinvSumSquares = 0;
TinvErr = 0;
// Do many tries, accumulate RMS
int N = 100000;
for(int iTry = 0; iTry<N; iTry++){
// Find the smeared matrix
TMatrixD Tsmeared = T;
for(int i = 0; i<nRow; i++){
for(int j = 0; j<nCol; j++){
double central = T(i,j);
double sigPos = TErrPos(i,j);
double sigNeg = TErrNeg(i,j);
// Switch to symmetric errors: approximation, but much simpler
double sig = (sigPos+sigNeg)/2.0;
Tsmeared(i,j) = gRandom->Gaus(central,sig);
}
}
// Find the inverted to smeared matrix
TMatrixD TinvSmeared = Tsmeared;
TinvSmeared.Invert(&det);
// Accumulate sum and sum of squares for each element
for(int i2 = 0; i2<nRow; i2++){
for(int j2 = 0; j2<nCol; j2++){
TinvSum (i2,j2) += TinvSmeared(i2,j2);
TinvSumSquares(i2,j2) += TinvSmeared(i2,j2)*TinvSmeared(i2,j2);
}
}
}
// Calculate the error matrix
TMatrixD TinvAverage = TinvSum;
for(int i = 0; i<nRow; i++){
for(int j = 0; j<nCol; j++){
TinvErr(i,j) = sqrt( TinvSumSquares(i,j)/double(N)
- (TinvSum(i,j)/double(N))*(TinvSum(i,j)/double(N)) );
TinvAverage(i,j) = TinvSum(i,j)/double(N);
}
}
return;
}示例2: MultiGaus
void MultiGaus(const TVectorD& parMeans, const TMatrixDSym& covMatrix, TVectorD& genPars)
{
TRandom3 rnd(0);
int nPars = parMeans.GetNrows();
if(nPars <= 0) {
Error("MultiGaus", "Must have >0 pars");
return;
}
if(covMatrix.GetNrows() != nPars) {
Error("MultiGaus", "parMeans.GetNrows() != covMatrix.GetNrows()");
return;
}
// Check that covMatrix is symmetric
for(int iRow = 0; iRow < nPars; iRow++) {
for(int iCol = iRow; iCol < nPars; iCol++) {
if(covMatrix(iRow, iCol) != covMatrix(iCol, iRow)) {
Error("MultiGaus", "malformed cov matrix at row %d, col %d", iRow, iCol);
return;
}
}
}
genPars.ResizeTo(nPars);
TMatrixDSymEigen eigenvariances(covMatrix);
TMatrixD V = eigenvariances.GetEigenVectors();
TVectorD rotParMeans = V * parMeans;
for(int iPar = 0; iPar < nPars; iPar++) {
double variance = eigenvariances.GetEigenValues()[iPar];
// check for positive-definiteness of covMatrix
if(variance < 0) {
Error("MultiGaus", "Got a negative eigenvariance (%f) on iPar = %d", variance, iPar);
}
genPars[iPar] = rnd.Gaus(rotParMeans[iPar], sqrt(variance));
}
V.Invert();
genPars = V * genPars;
}示例3: bigmatrix_corr
//.........这里部分代码省略.........
dumpElements(jesSys_ele);
dumpElements(jesSys_muo);
// print out measurement value
// dumpElements(measurement);
// // print out error matrix component
dumpElements(errorM);
// assign the correlation between different bins due to JES
// 100% correlation
for(int ibin=0; ibin < nbins; ibin++){
for(int jbin=0; jbin<nbins; jbin++)
{
if(jbin!=ibin){
// electron channel
errorM(ibin,jbin)= correlation*jesSys_ele(ibin)*jesSys_ele(jbin);
// muon channel
errorM(nbins+ibin,nbins+jbin)= correlation*jesSys_muo(ibin)*
jesSys_muo(jbin);
// electron-muon channel
errorM(ibin,nbins+jbin)= correlation*jesSys_ele(ibin)*jesSys_muo(jbin);
// muon-electron channel
errorM(nbins+ibin,jbin)= correlation*jesSys_muo(ibin)*jesSys_ele(jbin);
}
}
}
dumpElements(errorM);
TMatrixD errorInverse = errorM;
double* det;
errorInverse.Invert(det);
dumpElements(errorInverse);
// TMatrixD Unit = errorInverse*errorM;
// dumpElements(Unit);
TMatrixD matrixRight(nbins,NELE);
matrixRight = transposeU*errorInverse;
// dumpElements(matrixRight);
TMatrixD matrixLeft(nbins,nbins);
matrixLeft = transposeU*(errorInverse*U);
// dumpElements(matrixLeft);
TMatrixD matrixLeftInverse = matrixLeft;
double* det2;
matrixLeftInverse.Invert(det2);
// dumpElements(matrixLeftInverse);
TMatrixD lambda(nbins,NELE);
lambda= matrixLeftInverse*matrixRight;
// dumpElements(lambda);
TMatrixD transposeLambda(NELE,nbins);
transposeLambda.Transpose(lambda);
// dumpElements(transposeLambda);
// dumpElements(lambda);
TVectorD combined_value(nbins);
combined_value = lambda*measurement;
// dumpElements(combined_value);示例4: plotDYUnfoldingMatrix
//.........这里部分代码省略.........
} else {
DetCorrFactor(i) = 0;
DetCorrFactorErrPos(i) = 0;
DetCorrFactorErrNeg(i) = 0;
}
}
// Find response matrix, which is simply the normalized migration matrix
for(int ifsr = 0; ifsr < DetMigration.GetNrows(); ifsr++){
double nEventsInFsrMassBin = 0;
for(int ireco = 0; ireco < DetMigration.GetNcols(); ireco++)
nEventsInFsrMassBin += DetMigration(ifsr,ireco);
// Now normalize each element and find errors
for(int ireco = 0; ireco < DetMigration.GetNcols(); ireco++){
tCentral = 0;
tErrPos = 0;
tErrNeg = 0;
// Note: the event counts supposedly are dominated with events with weight "1"
// coming from the primary MC sample, so the error is assumed Poissonian in
// the call for efficiency-calculating function below.
if( nEventsInFsrMassBin != 0 ){
EfficiencyDivide(DetMigration(ifsr,ireco), nEventsInFsrMassBin, tCentral, tErrNeg, tErrPos);
// BayesEfficiency does not seem to be working in newer ROOT versions,
// so it is replaced by simler method
// BayesEfficiency( DetMigration(ifsr,ireco), nEventsInFsrMassBin, tCentral, tErrNeg, tErrPos);
}
DetResponse (ifsr,ireco) = tCentral;
DetResponseErrPos(ifsr,ireco) = tErrPos;
DetResponseErrNeg(ifsr,ireco) = tErrNeg;
}
}
// Find inverted response matrix
TMatrixD DetInvertedResponse = DetResponse;
Double_t det;
DetInvertedResponse.Invert(&det);
TMatrixD DetInvertedResponseErr(DetInvertedResponse.GetNrows(), DetInvertedResponse.GetNcols());
calculateInvertedMatrixErrors(DetResponse, DetResponseErrPos, DetResponseErrNeg, DetInvertedResponseErr);
// Package bin limits into TVectorD for storing in a file
TVectorD BinLimitsArray(DYTools::nMassBins+1);
for(int i=0; i<DYTools::nMassBins+1; i++)
BinLimitsArray(i) = DYTools::massBinLimits[i];
// Store constants in the file
TString outputDir(TString("../root_files/constants/")+dirTag);
if((systematicsMode==DYTools::RESOLUTION_STUDY) || (systematicsMode==DYTools::FSR_STUDY))
outputDir = TString("../root_files/systematics/")+dirTag;
gSystem->mkdir(outputDir,kTRUE);
TString unfoldingConstFileName(outputDir+TString("/unfolding_constants.root"));
if(systematicsMode==DYTools::RESOLUTION_STUDY){
unfoldingConstFileName = outputDir+TString("/unfolding_constants_seed_");
unfoldingConstFileName += seed;
unfoldingConstFileName += ".root";
}
if(systematicsMode==DYTools::FSR_STUDY){
unfoldingConstFileName = outputDir+TString("/unfolding_constants_reweight_");
unfoldingConstFileName += int(100*reweightFsr);
unfoldingConstFileName += ".root";
}
TFile fConst(unfoldingConstFileName, "recreate" );
DetResponse .Write("DetResponse");
DetInvertedResponse .Write("DetInvertedResponse");
DetInvertedResponseErr .Write("DetInvertedResponseErr");
BinLimitsArray .Write("BinLimitsArray");示例5: PredictionStraightLine
//#############################################################################
double* TrTrackA::PredictionStraightLine(int index){
static double pred[7];
#pragma omp threadprivate (pred)
int _Nhit = _Hit.size();
// Consistency check on the number of hits
if (_Nhit<3) return NULL;
// Get hit positions and uncertainties
// Scale errors (we will use sigmas in microns)
double hits[_Nhit][3];
double sigma[_Nhit][3];
for (int i=0;i<_Nhit;i++){
for (int j=0;j<3;j++){
hits[i][j] = _Hit[i].Coo[j];
sigma[i][j] = 1.e4*_Hit[i].ECoo[j];
if (i==index) sigma[i][j] *= 1.e2;
}
}
// Lenghts
double lenz[_Nhit];
for (int i=0;i<_Nhit;i++) lenz[i] = hits[i][2] - _Hit[0].Coo[2];
// F and G matrices
const int idim = 4;
double d[2*_Nhit][idim];
for (int i=0;i<_Nhit;i++) {
int ix = i;
int iy = i+_Nhit;
for (int j=0;j<idim;j++) { d[ix][j] = 0; d[iy][j] = 0;}
d[ix][0] = 1.;
d[iy][1] = 1.;
d[ix][2] = lenz[i];
d[iy][3] = lenz[i];
}
// F*S_x*x + G*S_y*y
double dx[idim];
for (int j=0;j<idim;j++) {
dx[j] = 0.;
for (int l=0;l<_Nhit;l++) {
dx[j] += d[l][j]/sigma[l][0]/sigma[l][0]*hits[l][0];
dx[j] += d[l+_Nhit][j]/sigma[l][1]/sigma[l][1]*hits[l][1];
}
}
// (F*S_x*F + G*S_y*G)
double Param[idim];
double InvCov[idim][idim];
for (int j=0;j<idim;j++) {
for (int k=0;k<idim;k++) {
InvCov[j][k] = 0.;
for (int l=0;l<_Nhit;l++) {
InvCov[j][k] += d[l][j]/sigma[l][0]/sigma[l][0]*d[l][k];
InvCov[j][k] += d[l+_Nhit][j]/sigma[l][1]/sigma[l][1]*d[l+_Nhit][k];
}
}
}
// (F*S_x*F + G*S_y*G)**{-1}
double determ = 0.0;
TMatrixD ParaCovariance = TMatrixD(idim,idim,(Double_t*)InvCov," ");
ParaCovariance = ParaCovariance.Invert(&determ);
if (determ<=0) return NULL;
// Solution
for (int k=0;k<idim;k++) {
Param[k] = 0.;
for (int i=0;i<idim;i++) {
Param[k] += ParaCovariance(k,i)*dx[i];
}
}
// Chi2 (xl and yl in microns, since sigmas are in microns too)
pred[0] = 0;
pred[1] = 0;
pred[2] = _Hit[index].Coo[2];
pred[3] = acos(-sqrt(1-Param[2]*Param[2]-Param[3]*Param[3]));
pred[4] = atan2(Param[3],Param[2]);
pred[5] = 0.;
pred[6] = 0.;
for (int k=0;k<idim;k++) {
pred[0] += d[index][k]*Param[k];
pred[1] += d[index+_Nhit][k]*Param[k];
for (int l=0;l<idim;l++) {
pred[5] += d[index][k]*d[index][l]*ParaCovariance(k,l);
pred[6] += d[index+_Nhit][k]*d[index+_Nhit][l]*ParaCovariance(k,l);
}
}
if (pred[5]>0.) pred[5] = 1.e-4*sqrt(pred[5]); else pred[5]=0.;
if (pred[6]>0.) pred[6] = 1.e-4*sqrt(pred[6]); else pred[6]=0.;
return pred;
}示例6: Prediction
//#############################################################################
double* TrTrackA::Prediction(int index){
static double pred[7];
#pragma omp threadprivate (pred)
int _Nhit = _Hit.size();
// Consistency check on the number of hits
if (_Nhit<4) return NULL;
// Scale errors (we will use sigmas in microns)
double hits[_Nhit][3];
double sigma[_Nhit][3];
for (int i=0;i<_Nhit;i++){
for (int j=0;j<3;j++){
hits[i][j] = _Hit[i].Coo[j];
sigma[i][j] = 1.e4*_Hit[i].ECoo[j];
if (i==index) sigma[i][j] *= 1.e2;
}
}
if (!_PathIntExist) SetPathInt();
// F and G matrices
const int idim = 5;
double d[2*_Nhit][idim];
for (int i=0;i<_Nhit;i++) {
int ix = i;
int iy = i+_Nhit;
for (int j=0;j<idim;j++) { d[ix][j] = 0; d[iy][j] = 0;}
d[ix][0] = 1.;
d[iy][1] = 1.;
for (int k=0;k<=i;k++) {
d[ix][2] += _PathLength[k];
d[iy][3] += _PathLength[k];
d[ix][4] += _PathLength[k]*_PathLength[k]*_PathIntegral_x[0][k];
d[iy][4] += _PathLength[k]*_PathLength[k]*_PathIntegral_x[1][k];
for (int l=k+1;l<=i;l++) {
d[ix][4] += _PathLength[k]*_PathLength[l]*_PathIntegral_u[0][k];
d[iy][4] += _PathLength[k]*_PathLength[l]*_PathIntegral_u[1][k];
}
}
}
// F*S_x*x + G*S_y*y
double dx[idim];
for (int j=0;j<idim;j++) {
dx[j] = 0.;
for (int l=0;l<_Nhit;l++) {
dx[j] += d[l][j]/sigma[l][0]/sigma[l][0]*hits[l][0];
dx[j] += d[l+_Nhit][j]/sigma[l][1]/sigma[l][1]*hits[l][1];
}
}
// (F*S_x*F + G*S_y*G)
double Param[idim];
double InvCov[idim][idim];
for (int j=0;j<idim;j++) {
for (int k=0;k<idim;k++) {
InvCov[j][k] = 0.;
for (int l=0;l<_Nhit;l++) {
InvCov[j][k] += d[l][j]/sigma[l][0]/sigma[l][0]*d[l][k];
InvCov[j][k] += d[l+_Nhit][j]/sigma[l][1]/sigma[l][1]*d[l+_Nhit][k];
}
}
}
// (F*S_x*F + G*S_y*G)**{-1}
double determ = 0.0;
TMatrixD ParaCovariance = TMatrixD(idim,idim,(Double_t*)InvCov," ");
ParaCovariance = ParaCovariance.Invert(&determ);
if (determ<=0) return NULL;
// Solution
//printf(">>>>>>>>>>>>> Cov AFTER:\n");
for (int k=0;k<idim;k++) {
Param[k] = 0.;
for (int i=0;i<idim;i++) {
Param[k] += ParaCovariance[k][i]*dx[i];
//printf(" %f", ParaCovariance[k][i]);
}
//printf("\n");
}
// Chi2 (xl and yl in microns, since sigmas are in microns too)
pred[0] = 0;
pred[1] = 0;
pred[2] = hits[index][2];
pred[3] = acos(-sqrt(1-Param[2]*Param[2]-Param[3]*Param[3]));
pred[4] = atan2(Param[3],Param[2]);
pred[5] = 0.;
pred[6] = 0.;
for (int k=0;k<idim;k++) {
pred[0] += d[index][k]*Param[k];
pred[1] += d[index+_Nhit][k]*Param[k];
for (int l=0;l<idim;l++) {
pred[5] += d[index][k]*d[index][l]*ParaCovariance(k,l);
pred[6] += d[index+_Nhit][k]*d[index+_Nhit][l]*ParaCovariance(k,l);
//.........这里部分代码省略.........
示例7: StraightLineFit
//#############################################################################
int TrTrackA::StraightLineFit(){
int _Nhit = _Hit.size();
// Reset fit values
Chi2 = FLT_MAX;
Ndof = 2*_Nhit-4;
Theta = 0.0;
Phi = 0.0;
U[0] = 0.0;
U[1] = 0.0;
U[2] = 1.0;
P0[0] = 0.0;
P0[1] = 0.0;
P0[2] = 0.0;
Rigidity = FLT_MAX;
// Consistency check on the number of hits
if (_Nhit<3) return -2;
P0[2] = _Hit[0].Coo[2];
// Get hit positions and uncertainties
// Scale errors (we will use sigmas in microns)
double hits[_Nhit][3];
double sigma[_Nhit][3];
for (int i=0;i<_Nhit;i++){
for (int j=0;j<3;j++){
hits[i][j] = _Hit[i].Coo[j];
sigma[i][j] = 1.e4*_Hit[i].ECoo[j];
}
}
// Lenghts
double lenz[_Nhit];
for (int i=0;i<_Nhit;i++) lenz[i] = hits[i][2] - P0[2];
// F and G matrices
const int idim = 4;
double d[2*_Nhit][idim];
for (int i=0;i<_Nhit;i++) {
int ix = i;
int iy = i+_Nhit;
for (int j=0;j<idim;j++) { d[ix][j] = 0; d[iy][j] = 0;}
d[ix][0] = 1.;
d[iy][1] = 1.;
d[ix][2] = lenz[i];
d[iy][3] = lenz[i];
}
// F*S_x*x + G*S_y*y
double dx[idim];
for (int j=0;j<idim;j++) {
dx[j] = 0.;
for (int l=0;l<_Nhit;l++) {
dx[j] += d[l][j]/sigma[l][0]/sigma[l][0]*hits[l][0];
dx[j] += d[l+_Nhit][j]/sigma[l][1]/sigma[l][1]*hits[l][1];
}
}
// (F*S_x*F + G*S_y*G)
double Param[idim];
double InvCov[idim][idim];
for (int j=0;j<idim;j++) {
for (int k=0;k<idim;k++) {
InvCov[j][k] = 0.;
for (int l=0;l<_Nhit;l++) {
InvCov[j][k] += d[l][j]/sigma[l][0]/sigma[l][0]*d[l][k];
InvCov[j][k] += d[l+_Nhit][j]/sigma[l][1]/sigma[l][1]*d[l+_Nhit][k];
}
}
}
// (F*S_x*F + G*S_y*G)**{-1}
double determ = 0.0;
TMatrixD ParaCovariance = TMatrixD(idim,idim,(Double_t*)InvCov," ");
ParaCovariance = ParaCovariance.Invert(&determ);
if (determ<=0) return -1;
// Solution
for (int k=0;k<idim;k++) {
Param[k] = 0.;
for (int i=0;i<idim;i++) {
Param[k] += ParaCovariance(k,i)*dx[i];
}
}
// Chi2 (xl and yl in microns, since sigmas are in microns too)
Chi2 = 0.;
for (int l=0;l<_Nhit;l++) {
double xl = hits[l][0]*1.e4;
double yl = hits[l][1]*1.e4;
for (int k=0;k<idim;k++) {
xl -= d[l][k]*Param[k]*1.e4;
yl -= d[l+_Nhit][k]*Param[k]*1.e4;
}
Chi2 += xl/sigma[l][0]/sigma[l][0]*xl + yl/sigma[l][1]/sigma[l][1]*yl;
}
// Final result
//.........这里部分代码省略.........