C++ TMatrixD类代码示例

void plot_covariance(TMatrixD cov, string path, string prefix, string name, string legend_position = "top_left", bool detail = false)
{
//declaring the canvas
    if (detail) { cout << "Ploting " << name << endl; }
    TCanvas *c1 = new TCanvas("c1","Canvas",0,29,1450,870);
    gStyle->SetOptStat(0);
    gStyle->SetOptTitle(kFALSE);
    gStyle->SetPalette(1);
    gStyle->SetPaintTextFormat("4.2g");
    gPad->SetFillColor(0);
    gPad->SetBorderMode(0);
    gPad->SetBorderSize(2);
    gPad->SetLeftMargin(0.10);
    gPad->SetRightMargin(0.10);
    gPad->SetTopMargin(0.01);
    gPad->SetFrameBorderMode(0);
    gPad->SetLogz();

    cov.Draw("colz");
    cov.Draw("text same");


//setting the output files
   string fileout = prefix + name;
   string out_png = path + "png/" + fileout + ".png";
   string out_c = path + "c/" + fileout + ".C";
   string out_eps = path + "eps/" + fileout + ".eps";
    
//save the file and close the canvas
    c1->Print( out_png.c_str() );
    c1->Print( out_c.c_str() );
    c1->Print( out_eps.c_str() );
    c1->Close();

}

   void CheckPoints(){
      std::vector<int> RemovedPoints;
      TMatrixD y_cov_new;
      int j = 0;
    
      for(int i = 0; i < y_val.size(); i++) {
         //         std::cout << "i: " << i << "   j: " << j << std::endl;
         if( y_val.at(i) == 0) {
            x_val.erase(x_val.begin()+i);
            y_val.erase(y_val.begin()+i);
            RemovedPoints.push_back(j);
            i = i-1;
         }
         j++;
      }    
    //   for(int i = 0; i < x_val.size(); i++) {
//           std::cout << "x: " << x_val.at(i) << std::endl;

//       }

      //  std::cout<< "Removed Points: " << RemovedPoints.size() << std::endl;
      //  std::cout<< "Remaining Points: " << x_val.size() << std::endl;

      y_cov_new.ResizeTo(x_val.size(),x_val.size());
      for( int i=0; i < x_val.size(); i++) {
         for(int k= 0; k < x_val.size(); k++) {
            y_cov_new(i,k) = y_cov(i+RemovedPoints.size(),k+RemovedPoints.size());
         }
      }
      y_cov.ResizeTo(0,0);
      y_cov.ResizeTo(x_val.size(),x_val.size());
      y_cov = y_cov_new;  
   }

TLorentzVector leptonic_fitter_algebraic::lv_with_mass( const TMatrixD& mat, double mass )
{
  assert( mat.GetNcols() == 1 || mat.GetNrows() == 1 );
  TLorentzVector lv;
  const double *array = mat.GetMatrixArray();
  lv.SetXYZM( array[0], array[1], array[2], mass );
  return lv;
}

TMatrix EUTelState::computePropagationJacobianFromLocalStateToNextLocalState(TVector3 positionEnd, TVector3 momentumEnd, float arcLength,float nextPlaneID) {
	streamlog_out(DEBUG2) << "-------------------------------EUTelState::computePropagationJacobianFromStateToThisZLocation()-------------------------BEGIN" << std::endl;
	if(arcLength == 0 or arcLength < 0 ){ 
		throw(lcio::Exception( Utility::outputColourString("The arc length is less than or equal to zero. ","RED"))); 
	}
	TMatrixD curvilinearJacobian = geo::gGeometry().getPropagationJacobianCurvilinear(arcLength,getOmega(), computeCartesianMomentum().Unit(),momentumEnd.Unit());
	streamlog_out(DEBUG0)<<"This is the curvilinear jacobian at sensor:" << std::endl; 
	streamlog_message( DEBUG0, curvilinearJacobian.Print();, std::endl; );

int addInQuadrature(const TMatrixD &sf, const TMatrixD &sfErr, const TMatrixD &sfSystRelErr, TMatrixD &totErr) {
  for (int ir=0; ir<sf.GetNrows(); ++ir) {
    for (int ic=0; ic<sf.GetNcols(); ++ic) {
      double err2=pow(sfErr(ir,ic),2) + pow(sf(ir,ic)*sfSystRelErr(ir,ic),2);
      totErr(ir,ic)=sqrt(err2);
    }
  }
  return 1;
}

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;
}

//==========================================================================
//
//  Dump element content of TMatrixD and TVectorD
//
//==========================================================================
void dumpElements(TMatrixD& a)
{
  cout << endl << endl;
  const int nrows = a.GetNrows();
  const int ncols = a.GetNcols();
  if(nrows==ncols)
    cout << "determinent = " << a.Determinant() << endl;
  a.Print();
  cout << endl << endl;

  return;
}

void leptonic_fitter_algebraic::real_eigenvalues( const TMatrixD& mat )
{
  Int_t ncol = mat.GetNcols();
  assert( ncol < 10 && mat.GetNrows() == ncol );

  _real_eigs.clear();

  TMatrixDEigen eigen( mat );
  const TVectorD& e_res = eigen.GetEigenValuesRe();
  const TVectorD& e_ims = eigen.GetEigenValuesIm();
  for( int ir=0; ir < ncol; ++ir ) {
    if( e_ims[ ir ] == 0 ) _real_eigs.push_back( e_res[ ir ] );
  }
}

double leptonic_fitter_algebraic::cofactor( const TMatrixD& mat, int i, int j )
{
  assert( mat.GetNcols() == 3 && mat.GetNrows() == 3 );
  assert( i >= 0 && i < 3 && j >= 0 && j < 3 );

  int i0 = ( i      ) ? 0 : 1;
  int i1 = ( i == 2 ) ? 1 : 2;
  int j0 = ( j      ) ? 0 : 1;
  int j1 = ( j == 2 ) ? 1 : 2;

  int sign = ( ( i+j ) % 2 ) ? -1 : 1;

  return sign * ( mat[i0][j0] * mat[i1][j1] - mat[i0][j1] * mat[i1][j0] );
}

void example() {
   ROOT::R::TRInterface &r = ROOT::R::TRInterface::Instance();
   // print R version
   r.Execute("print(version$version.string)");

   // compute standard deviation of 1000 vector normal numbers

   double std_dev_r  = r.Eval("sd(rnorm(10000))");
   std::vector<double> v = r.Eval("rnorm(10000)");
   double std_dev_root = TMath::StdDev(v.begin(),v.end());
   std::cout << "standard deviation from R    = " << std_dev_r << std::endl;
   std::cout << "standard deviation from ROOT = " <<  std_dev_root << std::endl;
   if (!TMath::AreEqualAbs(std_dev_r,std_dev_root,0.1))
      Error("ROOT-R-Example","Different std-dev found");

   // use << to execute the R command instead of Execute
   r << "mat<-matrix(c(1,2,3,4,5,6),2,3,byrow=TRUE)";
   TMatrixD m = r["mat"];
   std::array<double,6> a = r.Eval("seq(1:6)");
   TMatrixD m2(2,3,a.data());

   if  (!(m==m2)) {
      Error("ROOT-R-Example","Different matrix  found");
      m.Print();
      m2.Print();
   }

   // example on how to pass ROOT objects to R
   std::vector<double> v_root{1,2,3,4,5,6,7,8};
   r["v"] = v_root;
   r << "v2<-seq(1:8)";
   bool isEqual = r.Eval("prod(v==v2)");
   if (!isEqual) {
      Error("ROOT-R-Example","Different vector created");
      r << "print(v)";
      r << "print(v2)";
   }

   // example on how to pass functions to R

   r["gaus_pdf"] = ROOT::Math::normal_pdf;

   r << "y<-gaus_pdf(0,1,1)";
   double value_r = r["y"];
   double value_root = ROOT::Math::normal_pdf(0,1,1);
   std::cout << "Function gaussian(0,1,1) evaluated in  R    = " << value_r << std::endl;
   std::cout << "Function gaussian(0,1,1) evaluated in ROOT  = " <<  value_root << std::endl;
   if (value_r != value_root)
      Error("ROOT-R-Example","Different function value forund in r = %f and ROOT = %f", value_r, value_root);
}

void leptonic_fitter_algebraic::swap_x_y( const TMatrixD& in, TMatrixD& out )
{
  if( in.GetNoElements() != 9 ) 
    throw std::runtime_error( "ERROR! leptonic_fitter_algebraic::swap_x_y handles only 3 by 3 matrices");

  if( &in == &out )
    throw std::runtime_error( "ERROR! leptonic_fitter_algebraic::swap_x_y needs different in and out objects");

  static int remap[ 9 ]={ 4, 3, 5,   1, 0, 2,  7, 6, 8 };
  static double vec_out[9];

  const double *vec_in = in.GetMatrixArray();
  for( int ic=0;ic<9;++ic ) vec_out[ ic ] = vec_in[ remap[ ic ] ];
  out.Use( 3, 3, vec_out );
}

int GetMassProfile(const TMatrixD &m, int yIdx, TVectorD &v) {
  v=0;
  for (int i=0; i<m.GetNrows(); ++i) {
    v[i] = m(i,yIdx);
  }
  return 1;
}

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;

}

//==========================================================================
//
//  Dump element content of TMatrixD and TVectorD
//
//==========================================================================
void dumpElements(TMatrixD& a)
{
  cout << endl << endl;
  const int nrows = a.GetNrows();
  const int ncols = a.GetNcols();
  if(nrows==ncols)
    cout << "determinent = " << a.Determinant() << endl;
  cout << "------------------------------------------------------------------"
       << endl;
  for(int i=0; i< nrows; i++)
    {
      for(int j=0; j< ncols; j++)
	cout << a(i,j) << ", \t";    
	
      cout << endl;
    }
  cout << "------------------------------------------------------------------"
       << endl;
  cout << endl << endl;

  return;
}

int LoadMatrix(const TString &fname, TMatrixD &M, TMatrixD &Merr, const TString &field, const TString &fieldErr) {
  TFile f(fname,"read");
  TMatrixD *Mptr=(TMatrixD*)f.Get(field);
  TMatrixD *MerrPtr=(TMatrixD*)f.Get(fieldErr);
  f.Close();
  if (!Mptr || !MerrPtr) {
    std::cout << "error in loading from <" << fname << ">\n";
    if (!Mptr) std::cout << " failed to load <" << field << ">\n";
    if (!MerrPtr) std::cout << " failed to load <" << fieldErr << ">\n";
    return 0;
  }
  if ((Mptr->GetNrows() != M.GetNrows()) ||
      (Mptr->GetNcols() != M.GetNcols()) ) {
    std::cout << "dim mismatch in <" << fname << ">\n";
    return 0;
  }
  M = *Mptr;
  Merr= *MerrPtr;
  delete Mptr;
  delete MerrPtr;
  return 1;
}