C++ TVectorD::GetNrows方法代码示例

TGraphErrors* plotTG(int i,bool isSum, int ntheta, TString dirname, int marker, int color){
	if(ntheta==5){
		if(!isSum)
        		TFile *f = TFile::Open(Form("%s/M%d%d/mergedv_Prod2.root",dirname.Data(),trkpointmax[i],trkpointmin[i]));
		else
        		TFile *f = TFile::Open(Form("%s/M%d%d/mergedv_Prod.root",dirname.Data(),trkpointmax[i],trkpointmin[i]));
	}
	else{
	        if(!isSum)
                        TFile *f = TFile::Open(Form("theta%d/%s/M%d%d/mergedv_Prod2.root",ntheta,dirname.Data(),trkpointmax[i],trkpointmin[i]));
                else
                        TFile *f = TFile::Open(Form("theta%d/%s/M%d%d/mergedv_Prod.root",ntheta,dirname.Data(),trkpointmax[i],trkpointmin[i]));
	}

        TVectorD *vecDv2 = (TVectorD*)f->Get(Form("D_%d/vmean",ibin));
        TVectorD *vecDv2err = (TVectorD*)f->Get(Form("D_%d/deltavmean",ibin));
        TVectorD *vecDavgpt = (TVectorD*)f->Get(Form("D_%d/avgpt",ibin));

        double *avgpt = vecDavgpt->GetMatrixArray();
        double *v2 = vecDv2->GetMatrixArray();
        double *v2err = vecDv2err->GetMatrixArray();
        int npt = vecDavgpt->GetNrows();
	TGraphErrors *gr=new TGraphErrors(npt,avgpt,v2,0,v2err);
	gr->SetMarkerSize(1.3);
        gr->SetMarkerStyle(marker);
        gr->SetMarkerColor(color);
        gr->SetLineColor(color);
	f->Close();
	return gr;
}

void DrawCell( TMVA::PDEFoamCell *cell, TMVA::PDEFoam *foam,
	       Double_t x, Double_t y,
	       Double_t xscale,  Double_t yscale )
{
   // recursively draw cell and it's daughters

   Float_t xsize = xscale*1.5;
   Float_t ysize = yscale/3;
   if (xsize > 0.15) xsize=0.1; //xscale/2;
   if (cell->GetDau0() != NULL) {
      TLine *a1 = new TLine(x-xscale/4, y-ysize, x-xscale, y-ysize*2);
      a1->SetLineWidth(2);
      a1->Draw();
      DrawCell(cell->GetDau0(), foam, x-xscale, y-yscale, xscale/2, yscale);
   }
   if (cell->GetDau1() != NULL){
      TLine *a1 = new TLine(x+xscale/4, y-ysize, x+xscale, y-ysize*2);
      a1->SetLineWidth(2);
      a1->Draw();
      DrawCell(cell->GetDau1(), foam, x+xscale, y-yscale, xscale/2, yscale);
   }

   TPaveText *t = new TPaveText(x-xsize, y-ysize, x+xsize, y+ysize, "NDC");

   t->SetBorderSize(1);
   t->SetFillStyle(1);

   // draw all cell elements
   t->AddText( Form("Intg=%.5f", cell->GetIntg()) );
   t->AddText( Form("Var=%.5f", cell->GetDriv()) );
   TVectorD *vec = (TVectorD*) cell->GetElement();
   if (vec != NULL){
      for (Int_t i = 0; i < vec->GetNrows(); ++i) {
	 t->AddText( Form("E[%i]=%.5f", i, vec(i)) );
      }
   }

   if (cell->GetStat() != 1) {
      // cell is inactive --> draw split point
      t->SetFillColor( TColor::GetColor("#BBBBBB") );
      t->SetTextColor( TColor::GetColor("#000000") );

      // cell position and size
      TMVA::PDEFoamVect cellPosi(foam->GetTotDim()), cellSize(foam->GetTotDim());
      cell->GetHcub(cellPosi, cellSize);
      Int_t    kBest = cell->GetBest(); // best division variable
      Double_t xBest = cell->GetXdiv(); // best division point
      t->AddText( Form("dim=%i", kBest) );
      t->AddText( Form("cut=%.5g", foam->VarTransformInvers(kBest,cellPosi[kBest] + xBest*cellSize[kBest])) );
   } else {
      t->SetFillColor( TColor::GetColor("#DD0033") );
      t->SetTextColor( TColor::GetColor("#FFFFFF") );
   }

   t->Draw();

   return;
}

TMatrixD Chol (TVectorD& covV, TVectorD& newSig)
{
  int nCov = covV.GetNrows();
  int n = newSig.GetNrows();
  std::cout << nCov << " " << n << std::endl;
  if ( nCov != n*(n+1)/2. ) {
    std::cout << "vecTest: mismatch in inputs" << std::endl;
    return TMatrixD();
  }
  //
  // create modified vector (replacing std.dev.s)
  //
  TVectorD newCovV(covV);
  int ind(0);
  for ( int i=0; i<n; ++i ) {
    for ( int j=0; j<=i; ++j ) {
      if ( j==i )  newCovV[ind] = newSig(i);
      ++ind;
    }
  }
  return Chol(newCovV,newSig);
}

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;

}

   double operator() (double *x, double *p) {
      // 4 parameters
      int dim = X.GetNrows();
      int k = 0;
      for (int i = 0; i<dim; ++i) { X[i] = x[i] - p[k]; k++; }
      for (int i = 0; i<dim; ++i) {
         CovMat(i,i) = p[k]*p[k];
         k++;
      }
      for (int i = 0; i<dim; ++i) {
         for (int j = i+1; j<dim; ++j) {
            // p now are the correlations N(N-1)/2
               CovMat(i,j) = p[k]*sqrt(CovMat(i,i)*CovMat(j,j));
               CovMat(j,i) = CovMat(i,j);
               k++;
         }
      }
      if (debug) {
         X.Print();
         CovMat.Print();
      }

      double det = CovMat.Determinant();
      if (det <= 0) {
         Fatal("GausND","Determinant is <= 0 det = %f",det);
         CovMat.Print();
         return 0;
      }
      double norm = std::pow( 2. * TMath::Pi(), dim/2) * sqrt(det);
      // compute the gaussians
      CovMat.Invert();
      double fval  = std::exp( - 0.5 * CovMat.Similarity(X) )/ norm;

      if (debug) {
         std::cout << "det  " << det << std::endl;
         std::cout << "norm " << norm << std::endl;
         std::cout << "fval " << fval << std::endl;
      }

      return fval;
   }