C++ TMatrixD::Print方法代码示例

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

//==========================================================================
//
//  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 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);
}

//***############## main fitting Fxn ################ *****//
void FitPlotAndSave( char *Ifile ){
 
/** Plot Options***/	
//gROOT->Reset();
// gROOT->Clear();
gROOT->SetStyle("Plain") ;
gROOT->SetBatch(kFALSE);
gStyle->SetOptTitle(1);
gStyle->SetOptStat(0);
gStyle->SetOptFit(1);
gStyle->SetStatX(.89);
gStyle->SetStatY(.89) ;
gStyle->SetStatBorderSize(0);
//gStyle->SetOptStat(1111111)
gStyle->SetCanvasColor(kWhite);   // background is no longer mouse-dropping white
gStyle->SetPalette(1);        // blue to red false color palette. Use 9 for b/w
gStyle->SetCanvasBorderMode(0);     // turn off canvas borders
gStyle->SetPadBorderMode(0);
gStyle->SetPaintTextFormat("5.2f");  // What precision to put numbers if plotted with "TEXT"

// For publishing:
gStyle->SetLineWidth(2);
gStyle->SetTextSize(1.1);
gStyle->SetLabelSize(0.06,"xy");
gStyle->SetTitleSize(0.08,"xy");
gStyle->SetTitleOffset(1.2,"x");
gStyle->SetTitleOffset(1.0,"y");
gStyle->SetPadTopMargin(0.1);
gStyle->SetPadRightMargin(0.1);
gStyle->SetPadBottomMargin(0.16);
gStyle->SetPadLeftMargin(0.12);
TGaxis::SetMaxDigits(1); // Set Axis to be of the form 0.11 10^N


       TFile *ifile  = new TFile(Ifile);
        
	TF1 *fitFcn  = new TF1("fitFcn", mygaus, FitLowRange, FitHighRange, 3 );
	fitFcn->SetNpx(500);
	fitFcn->SetLineWidth(4);
	fitFcn->SetLineStyle(5);
	fitFcn->SetLineColor(kBlue);
        cout <<" Calling Fitting Fxntion" << endl;
	TH1F*h_Seed_TimeEBEB = (TH1F*)ifile->Get("EBEB/seed time");
	if(h_Seed_TimeEBEB == 0){ std::cout  <<"!! Histogram Does not exist!!" << std::endl; throw 1;}
        
	h_Seed_TimeEBEB->SetTitle("Seed Time[ns]");   
        h_Seed_TimeEBEB->SetMarkerStyle(20);
        h_Seed_TimeEBEB->SetMarkerSize(0.8);
        h_Seed_TimeEBEB->SetStats(1);
        h_Seed_TimeEBEB->SetTitleSize(0.08, "x");   
        h_Seed_TimeEBEB->SetTitleOffset(1.0, "x");    
        h_Seed_TimeEBEB->SetTitleSize(0.06, "y"); 
        h_Seed_TimeEBEB->SetTitleOffset(0.95, "y");    
        h_Seed_TimeEBEB->SetYTitle("Number of Seeds/0.05ns"); 
        h_Seed_TimeEBEB->SetXTitle("t_{seed}[ns]"); 
        h_Seed_TimeEBEB->GetXaxis()->SetRangeUser(FitLowRange, FitHighRange);   

       /** Set parms as parms of Fit Fxn **/
	fitFcn->SetParameters(500, h_Seed_TimeEBEB->GetMean(), h_Seed_TimeEBEB->GetRMS() );
	fitFcn->SetParNames("CONST", "#mu(ns)", "#sigma(ns)");
	h_Seed_TimeEBEB->Fit("fitFcn", "LL"); /**Fit with improved LL**/
	std::cout << "Printing Fit Parameters for EBEB ......   " << std::endl;
        printf("Integral of function in EBEB = %g\n", fitFcn->Integral( FitLowRange, FitHighRange));

        //*** retrive fit results***//
        int npar = fitFcn->GetNpar();
        TVirtualFitter *fit = TVirtualFitter::GetFitter();
        fit->PrintResults(2,0.);
        TMatrixD *CovMatrix = new TMatrixD ( npar, npar, fit->GetCovarianceMatrix() );
	CovMatrix->Print();
        TCanvas *c1 = new TCanvas("c1","EB-EB",200,10,800,900);
     	c1->SetGridx();
     	c1->SetGridy();
	c1->GetFrame()->SetFillColor(21);
	c1->GetFrame()->SetBorderMode(-1);
	c1->GetFrame()->SetBorderSize(5);
        /* c1->Divide(2,1);  */
	c1->cd();
	h_Seed_TimeEBEB->Draw();
	fitFcn->Draw("sames");
        c1->SetLogy(0);
	//  draw the legend
    	TLegend *leg = new TLegend(0.15,0.72,0.3,0.85);
       	leg->SetTextFont(72);
        leg->SetTextSize(0.04);
        leg->AddEntry(h_Seed_TimeEBEB,"EB","lpe");
        leg->AddEntry(fitFcn,"GAUS","l");
        leg->Draw();
	c1->SaveAs("Seed_Time_DoubleElectron_Run2012A-EB-EB.png");
}       

//.........这里部分代码省略.........
				CovPion[0][0] = cov[0];
				CovPion[0][1] = cov[1];
				CovPion[0][2] = cov[2];
				CovPion[0][3] = cov[3];
				CovPion[0][4] = cov[4];

				CovPion[1][1] = cov[5];
				CovPion[1][2] = cov[6];
				CovPion[1][3] = cov[7];
				CovPion[1][4] = cov[8];

				CovPion[2][2] = cov[9];
				CovPion[2][3] = cov[10];
				CovPion[2][4] = cov[11];

				CovPion[3][3] = cov[12];
				CovPion[3][4] = cov[13];

				CovPion[4][4] = cov[14];

				for (int i=0; i<5; i++){
					for(int j=0; j<5; j++){
						CovPion[j][i]=CovPion[i][j];
					}
				}

				cout << "Covariance matrix for pion: " << endl;

				for(int i=0; i<7; i++){
					for(int j=0; j<7; j++){
						if(TMath::Abs(piCov[i][j])<1e-6) piCov[i][j]=0;
					}
				}
				piCov.Print();

				cout << "First Par Covariance matrix for pion: " << endl;
				CovPion.Print();

			}

			//Get proton

			for (int j=0; j<proton.GetLength(); ++j){

				bool truth = theAnalysis->McTruthMatch(proton[j]);

				if(truth){
					TMatrixD protCov = proton[j]->Cov7();

					FairTrackParP paramFirst = theAnalysis->GetFirstPar(proton[j]);

					double cov[15];
					paramFirst.GetCov(cov);

					for (int i=0; i<15; i++){
						if(TMath::Abs(cov[i]<1e-6)) cov[i]=0;
					}

					CovProton[0][0] = cov[0];
					CovProton[0][1] = cov[1];
					CovProton[0][2] = cov[2];
					CovProton[0][3] = cov[3];
					CovProton[0][4] = cov[4];

					CovProton[1][1] = cov[5];
					CovProton[1][2] = cov[6];

//.........这里部分代码省略.........
   
   c1->cd();

   TPad *pad2 = new TPad("pad2","pad2",0,0,1,0.5);
   pad2->SetTopMargin(0);
   pad2->SetBottomMargin(0.4);
   pad2->Draw();
   pad2->cd();
   pad2->SetGrid();
   h2->SetStats(0);
   h2->Divide(h1);
   //h2->SetMarkerStyle(21);
   h2->Draw("ep");
   h2->GetYaxis()->SetLabelSize(0.1);
   h2->GetYaxis()->SetRangeUser(-0.5, 2.5);// ,yTopLimit);
   h2->GetYaxis()->SetTitle("[email protected]+PY8 / MG5+PY6");
   h2->GetYaxis()->SetTitleSize(0.06);
   h2->GetYaxis()->SetTitleOffset(0.7);
   h2->GetXaxis()->SetLabelSize(0.1);
   h2->GetXaxis()->SetTitle(x_title);
   h2->GetXaxis()->SetTitleSize(0.16);
   h2->GetXaxis()->SetTitleOffset(0.9);
 //  Double_t matrix[4][4];
   h2->Fit("pol3","","",50.0,1200.0);
   TF1 *ratio = h2->GetFunction("pol3");
   TVirtualFitter *fitter = TVirtualFitter::GetFitter();
   TMatrixD matrix(4,4,fitter->GetCovarianceMatrix());
   Double_t errorPar00 = fitter->GetCovarianceMatrixElement(0,0);
   Double_t errorPar11 = fitter->GetCovarianceMatrixElement(1,1);
   Double_t errorPar22 = fitter->GetCovarianceMatrixElement(2,2);
   Double_t errorPar33 = fitter->GetCovarianceMatrixElement(3,3);
//   c1->cd();

   matrix.Print();

   //const TMatrixDSym m = matrix;
   const TMatrixDEigen eigen(matrix);
   const TMatrixD eigenVal = eigen.GetEigenValues();
   const TMatrixD V = eigen.GetEigenVectors(); 

   cout << "V" << endl;

   V.Print();

   cout << "eigenVal" << endl;

   eigenVal.Print();



   cout << "Recomputed diag" << endl;

   //const TMatrixD Vt(TMatrixD::kTransposed,V);
   //const TMatrixD Vinv = V.Invert();
   const TMatrixD Vt(TMatrixD::kTransposed,V);
   //cout << "V-1" << endl;
   //Vinv.Print();
   cout << "Vt" << endl;
   Vt.Print();

   const TMatrixD VAVt = Vt*matrix*V;
   VAVt.Print();


   const TVectorD FittedParam(4);
   FittedParam(0) = fitter->GetParameter(0);

// Open a FITS file and retrieve the first plane of the image array 
// as a TASImage object
void FITS_tutorial1()
{
   printf("\n\n--------------------------------\n");
   printf("WELCOME TO FITS tutorial #1 !!!!\n");
   printf("--------------------------------\n");
   printf("We're gonna open a FITS file that contains only the\n");
   printf("primary HDU, consisting on an image.\n");
   printf("The object you will see is a snapshot of the NGC7662 nebula,\n");
   printf("which was taken by the author on November 2009 in Barcelona (CATALONIA).\n\n");
      
   if (!gROOT->IsBatch()) {
      //printf("Press ENTER to start..."); getchar();
   }
   
   // Open primary HDU from file
   TFITSHDU *hdu = new TFITSHDU("sample1.fits");
   if (hdu == 0) {
      printf("ERROR: could not access the HDU\n"); return;
   }
   printf("File successfully open!\n");
   
   // Dump the HDUs within the FITS file
   // and also their metadata
   //printf("Press ENTER to see summary of all data stored in the file:"); getchar();
   
   hdu->Print("F+");
   
   printf("....................................\n");
   // Here we get the exposure time.
   //printf("Press ENTER to retrieve the exposure time from the HDU metadata..."); getchar();
   printf("Exposure time = %s\n", hdu->GetKeywordValue("EXPTIME").Data());

   
   // Read the primary array as a matrix,
   // selecting only layer 0.
   // This function may be useful to
   // do image processing.
   printf("....................................\n");
   printf("We can read the image as a matrix of values.\n");
   printf("This feature is useful to do image processing, e.g:\n");
   printf("histogram equalization, custom filtering, ...\n");
   //printf("Press ENTER to continue..."); getchar();
   
   TMatrixD *mat = hdu->ReadAsMatrix(0);
   mat->Print();
   delete mat;
   
   // Read the primary array as an image,
   // selecting only layer 0.
   printf("....................................\n");
   printf("Now the primary array will be read both as an image and as a histogram,\n");
   printf("and they will be shown in a canvas.\n");
   //printf("Press ENTER to continue..."); getchar();
   
   TASImage *im = hdu->ReadAsImage(0);
   
   // Read the primary array as a histogram.
   // Depending on array dimensions, returned
   // histogram will be 1D, 2D or 3D
   TH1 *hist = hdu->ReadAsHistogram();
   
   
   TCanvas *c = new TCanvas("c1", "FITS tutorial #1", 800, 300);
   c->Divide(2,1);
   c->cd(1);
   im->Draw();
   c->cd(2);
   hist->Draw("COL");

   // Clean up
   delete hdu;
}