C++ TVectorD类代码示例

Parabola2D::Parabola2D(TVectorD x, TVectorD y, TVectorD z) : Parabola() {

  // Prepare the members
  int n = x.GetNoElements();
  double* vx = x.GetMatrixArray();
  double* vy = y.GetMatrixArray();
  double* vz = z.GetMatrixArray();
  vector<double> entry_x;
  vector<double> entry_y;
  vector<double> entry_z;
  for(int i=0; i<n; i++){
    entry_x.push_back(vx[i]);
    entry_y.push_back(vy[i]);
    entry_z.push_back(vz[i]);
  }
  entries.push_back(entry_z); // _Has_ to be the first.
  entries.push_back(entry_x);
  entries.push_back(entry_y);

  // And do the magic
  FillMatrix();    // Do the fitting...
  FindMinimum();   // Find the Minimum...
  
  // Write out the information. The order that the coefficients appear has to be in increasing key,
  // as explained above. THIS IS IMPORTANT.
  FitFunction = new TF2("FitFnc","[5]*y*y+[4]*y*x+[3]*y+[2]*x*x+[1]*x+[0]");
  for (int i = 0; i<6; i++) FitFunction->SetParameter(i,fit_coef_vec[i]);
  FitGraph = new TGraph2D(n,vx,vy,vz);

}

void plotr0vstheta(){
	TFile *f;
        int isSum=0;
	if(isSum==0){
	    f = TFile::Open("mergedV_Prod.root");
	}
	else{
	    f = TFile::Open("mergedV_Sum.root");
	}
	int xbin=0;

//	TVectorD* vecVmean = (TVectorD*)f->Get(Form("D_%d/Vmean",xbin));
        TVectorD* vecV = (TVectorD*)f->Get(Form("D_%d/D_0/r01",xbin));
        
        double *V = vecV->GetMatrixArray();
        double Vmean = getmean(V,ntheta);
        double theta[ntheta];
        for(int itheta=0;itheta<ntheta;itheta++){
            theta[itheta]=itheta*TMath::Pi()/ntheta/nn;

        }
        int maxper10 = findmaxper(V,ntheta,Vmean);
        double maxper = (double)(maxper10+2)/10;
        TGraph *gV2theta = new TGraph(ntheta,theta,V);
        gV2theta->SetTitle("");
        gV2theta->GetXaxis()->SetTitle("#theta");
        gV2theta->GetYaxis()->SetTitle("r_{0}");
        gV2theta->GetYaxis()->SetTitleOffset(1.1);
        gV2theta->GetXaxis()->SetTitleSize(0.04);
        gV2theta->GetYaxis()->SetTitleSize(0.04);
        gV2theta->SetMarkerStyle(20);
        gV2theta->SetMarkerSize(1.3);
        gV2theta->SetMarkerColor(1);
        gV2theta->SetLineColor(1);
        gV2theta->SetMinimum(Vmean*(1-maxper*1.5));
        gV2theta->SetMaximum(Vmean*(1+maxper*1.5));
	gV2theta->Draw("AP");
        TLine *lup = new TLine(gV2theta->GetXaxis()->GetXmin(),Vmean*(1+maxper), gV2theta->GetXaxis()->GetXmax(),Vmean*(1+maxper));
        TLine *ldown = new TLine(gV2theta->GetXaxis()->GetXmin(),Vmean*(1-maxper), gV2theta->GetXaxis()->GetXmax(),Vmean*(1-maxper));
        TLine *l = new TLine(gV2theta->GetXaxis()->GetXmin(),Vmean, gV2theta->GetXaxis()->GetXmax(),Vmean);
        l->SetLineStyle(2);
        lup->SetLineStyle(2);
        ldown->SetLineStyle(2);
        l->SetLineWidth(1.2);
        lup->SetLineWidth(1.2);
        ldown->SetLineWidth(1.2);
        TLatex *tl = new TLatex();
 //     tl->SetNDC();
        tl->SetTextFont(42);
        tl->SetTextSize(0.04);
 //     tl->SetBorderStyle(0);
        tl->DrawLatex(0,Vmean*(1+maxper),Form("mean up %.f%%",maxper*100));
        tl->DrawLatex(0,Vmean*(1-maxper),Form("mean down %.f%%",maxper*100));
        tl->DrawLatex(1,0.2,Form("Multiplicity %d to %d",120,150));
        l->Draw("same");
        lup->Draw("same");
        ldown->Draw("same");
	if(isSum==0)c1->SaveAs("hr0theta_Prod.png");
	else c1->SaveAs("hr0theta_Sum.png");
}

TGraph *plotGF(int isSum, int xtheta, double *r0_theta, double *G2_theta, int marker, int color){

if(isSum)
TFile *f = TFile::Open("mergedV_Sum.root");
else
TFile *f = TFile::Open("mergedV_Prod.root");

TVectorD *vecDr = f->Get(Form("D_%d/r",xbin));
TVectorD *vecDg2 = f->Get(Form("D_%d/D_0/D_%d/G2",xbin,xtheta));
TVectorD *vecDr0 = f->Get(Form("D_%d/D_0/r0",xbin,xtheta));
TVectorD *vecDr01 = f->Get(Form("D_%d/D_0/r01",xbin,xtheta));

double *r = vecDr->GetMatrixArray();
double *g2 = vecDg2->GetMatrixArray();
double *r0 = vecDr0->GetMatrixArray();
double *r01 = vecDr01->GetMatrixArray();
(*r0_theta) = r0[xtheta];
for(int ir=0;ir<nstepr;ir++)
    if(r[ir] == r01[xtheta]) break;
(*G2_theta) = g2[ir];
TGraph *gr = new TGraph(nstepr,r,g2);
gr->SetMarkerSize(0.5);
gr->SetMarkerColor(color);
gr->SetMarkerStyle(marker);
f->Close();
return gr;
}

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

double* getEPR(TString indir){
	int ibin=0;
	TFile *f = TFile::Open(Form("%s/mergedVobs.root",indir.Data()));
	TVectorD *vecDEPR = (TVectorD*)f->Get(Form("D_%d/EPR",ibin));
	double *EPR = vecDEPR->GetMatrixArray();
	return EPR;
}

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

void loadVectorsFromFile(const char *filename, 
			 const char *scanfmt,
			 TVectorD&   v1,
			 TVectorD&   v2,
			 TVectorD&   v3)
{
  char linein[LINELEN];
  vector<double> v;

  if (TString(scanfmt).Contains("%f")) {
    cerr << "Must use %lf format for doubles, sorry " << endl;
    exit(-1);
  }

  FILE *fp = fopen(filename, "r");

  if (!fp) {
    cerr << "File failed to open, " << filename << endl;
    exit(-1);
  }

  if (gl_verbose)
    cout << "Loading vectors from file " << filename
	 << " with scan format \"" << scanfmt << "\"" << endl;

  while (!feof(fp) && fgets(linein,LINELEN,fp)) {
    double x1, x2, x3;
    if( linein[0]=='#' ) continue;                // comments are welcome
    if( sscanf(linein, scanfmt, &x1, &x2, &x3) != 3 ) {
      cerr << "scan failed, file " << filename << ", line = " << linein;
      cerr << ", scanfmt = " << scanfmt;
      cerr << ", skipping" << endl;
      //return;
      continue;
    }
    else {
      v.push_back(x1); v.push_back(x2); v.push_back(x3);
    }
  }

  int vecsize = v.size()/3;
  v1.ResizeTo(vecsize);
  v2.ResizeTo(vecsize);
  v3.ResizeTo(vecsize);

  if (gl_verbose) cout << "; read " << vecsize << " lines" << endl;

  for (int i=0; i<vecsize; i++) {
    v1[i] = v[3*i];
    v2[i] = v[3*i+1];
    v3[i] = v[3*i+2];
  }

  fclose(fp);
}                                                 // loadVectorsFromFile

void loadVectorsFromFile(const char *filename, 
			 const char *scanfmt,
			 TVectorD&   vx,
			 TVectorD&   vy,
			 char        xheader[],
			 char        yheader[])
{
  char linein[LINELEN];
  vector<double> v;

  FILE *fp = fopen(filename, "r");

  if (!fp) {
    cerr << "File failed to open, " << filename << endl;
    return;
  }

  if (gl_verbose)
    cout << "Loading vectors from file " << filename;

  while (!feof(fp) && fgets(linein,LINELEN,fp)) {
    double x, y;
    if( linein[0]=='#' ) {
      if (!v.size()) {           // first line, try to read headers
	TString hscanfmt(scanfmt);
	hscanfmt.ReplaceAll("lf","s");
	if( sscanf(&linein[1],hscanfmt.Data(),xheader,yheader) != 2 ) {
	  cerr << "failed to read in column headers" << endl;
	}
      }
      continue;                // comments are welcome
    }

    if( sscanf(linein, scanfmt, &x, &y) != 2 ) {
      cerr << "scan failed, file " << filename << ", line = " << linein << endl;
      return;
    }
    else {
      v.push_back(x); v.push_back(y);
    }
  }

  int vecsize = v.size()/2;
  vx.ResizeTo(vecsize);
  vy.ResizeTo(vecsize);

  if (gl_verbose) cout << "; read " << vecsize << " lines" << endl;

  for (int i=0; i<vecsize; i++) {
    vx[i] = v[2*i];
    vy[i] = v[2*i+1];
  }
}                                                 // loadVectorsFromFile

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;

}

void writeDataBackgroundHistosForModel(const std::map<TString,TGraph*>& m_bkgds,
				       const std::vector<TH1D *>& vchans,
				       TFile  *allHistFile)
{
  for (std::map<TString,TGraph*>::const_iterator it = m_bkgds.begin();
       it != m_bkgds.end();
       it++) {
    const TString& name = it->first;

    // determine binning from the signal histogram for this channel
    // - (sigh) have to find it first...
    //
    TString channame = name.Copy().Remove(0,strlen("Bckgrndtot_"));
    TH1D *sigh=(TH1D*)0;
    for (int ichan=0; ichan<NUMCHAN; ichan++) {
      sigh = vchans.at(ichan);
      if (strstr(sigh->GetName(),channame.Data()))
	break;
    }
    assert (sigh);


    // for variable binning - all histos must have the same binning per channel
    TAxis *xax = sigh->GetXaxis();

    TVectorD xbins   = TVectorD(sigh->GetNbinsX(),xax->GetXbins()->GetArray());

    int lobin = xax->FindFixBin(sumwinmin);
    int hibin = xax->FindFixBin(sumwinmax)-1;
    int nbins = hibin-lobin+1;

    TVectorD xwindow = xbins.GetSub(lobin-1,hibin);
    
    printf("Booking TH1D(%s,%s,%d,xwindowarray)\n",
	   name.Data(),name.Data(),nbins);
    TH1D *h = new TH1D(name.Data(),name.Data(),nbins,xwindow.GetMatrixArray());

    for (int ibin=1; ibin <= nbins; ibin++)
      h->SetBinContent(ibin,
		       it->second->Eval(h->GetBinCenter(ibin))
		       * h->GetBinWidth(ibin)
		       );

    allHistFile->WriteTObject(h);
  }

}                                            //  writeDataBackgroundHistosForModel

void dumpElements(TVectorD& a)
{
  cout << endl << endl;
  a.Print();
  cout << endl << endl;

  return;
}

void plotV2vstheta(){
	TFile *f;
	if(isSum==0){
	    f = TFile::Open("mergedV_Prod.root");
	}
	else{
	    f = TFile::Open("mergedV_Sum.root");
	}
	int xbin=0;

	TVectorD* vecVmean = (TVectorD*)f->Get(Form("D_%d/Vmean",xbin));
        TVectorD* vecV = (TVectorD*)f->Get(Form("D_%d/D_0/V",xbin));

        double Vmean = (*vecVmean)[0];
        cout<<Vmean<<endl;
        double *V = vecV->GetMatrixArray();
        double theta[ntheta];
        for(int itheta=0;itheta<ntheta;itheta++){
            theta[itheta]=itheta*TMath::Pi()/ntheta/nn;
            cout<<V[itheta]<<endl;
        }
        TGraph *gV2theta = new TGraph(ntheta,theta,V);
        gV2theta->SetMarkerStyle(20);
        gV2theta->SetMarkerSize(1.3);
        gV2theta->SetMarkerColor(1);
        gV2theta->SetLineColor(1);
        gV2theta->SetMinimum(Vmean*0.98);
        gV2theta->SetMaximum(Vmean*1.02);
	gV2theta->Draw("AP");
        TLine *lup = new TLine(gV2theta->GetXaxis()->GetXmin(),Vmean*1.01, gV2theta->GetXaxis()->GetXmax(),Vmean*1.01);
        TLine *ldown = new TLine(gV2theta->GetXaxis()->GetXmin(),Vmean*0.99, gV2theta->GetXaxis()->GetXmax(),Vmean*0.99);
        TLine *l = new TLine(gV2theta->GetXaxis()->GetXmin(),Vmean, gV2theta->GetXaxis()->GetXmax(),Vmean);
        l->SetLineStyle(2);
        lup->SetLineStyle(2);
        ldown->SetLineStyle(2);
        l->SetLineWidth(1.2);
        lup->SetLineWidth(1.2);
        ldown->SetLineWidth(1.2);
        l->Draw("same");
        lup->Draw("same");
        ldown->Draw("same");
	if(isSum==0)c1->SaveAs("hV2theta_Prod.png");
	else c1->SaveAs("hV2theta_Sum.png");
}

TGraph *plotGF(int isSum, int xtheta, double *r0_theta, double *G2_theta, int marker, int color){

TFile *f = TFile::Open("mergedV.root");

TVectorD *vecDr = f->Get(Form("r"));
TVectorD *vecDGRe = f->Get(Form("GRe"));
TVectorD *vecDGIm = f->Get(Form("GIm"));
TVectorD *vecDG2 = f->Get(Form("G2"));

double *r = vecDr->GetMatrixArray();
double *GRe = vecDGRe->GetMatrixArray();
double *GIm = vecDGIm->GetMatrixArray();
double *G2 = vecDG2->GetMatrixArray();
TGraph *gr = new TGraph(nstepr,r,G2);
if(isSum){
gcl->SetNpx(10000);
//gcl->SetParameters(sigma2[xtheta]-inV2*inV2*avgmult[xbin]*avgmult[xbin],inV2*avgmult[xbin]);
//linv->SetX1(j01/avgmult[xbin]/inV2);
//linv->SetX2(j01/avgmult[xbin]/inV2);
//linv->SetY2(5e-9);
}
gr->SetMarkerSize(0.5);
gr->SetMarkerColor(color);
gr->SetMarkerStyle(marker);
f->Close();
return gr;
}

void plotV2vstheta(){
    gStyle->SetOptFit(1);
    gStyle->SetOptStat(0);
    gStyle->SetOptTitle(0);
    gStyle->SetErrorX(0);
    int xbin = 0;
    TFile *f = new TFile("mergedV_Prod.root");
    TVectorD *vecV = (TVectorD*)f->Get(Form("D_%d/V"));
    TVectorD *vecdeltaV = (TVectorD*)f->Get(Form("D_%d/deltaV"));
    TVectorD *vecVmean = (TVectorD*)f->Get(Form("Vmean"));
    double *V = vecV->GetMatrixArray();
    double *deltaV = vecdeltaV->GetMatrixArray();
    double *Vmean = vecVmean->GetMatrixArray();
    double theta[ntheta];
    for(int itheta=0;itheta<ntheta;itheta++){
        theta[itheta]=itheta*TMath::Pi()/ntheta/nn;
    }
    TH1D *hFrame = new TH1D("","",300,-1,2);
    hFrame->GetXaxis()->SetTitle("#theta");
    hFrame->GetYaxis()->SetTitle("referenceV_{2}");	
    hFrame->GetXaxis()->SetTitleSize(0.04);
    hFrame->GetYaxis()->SetTitleSize(0.04);
    hFrame->GetXaxis()->SetRangeUser(-0.1,1.5);
    hFrame->SetMaximum(0.055);
    hFrame->SetMinimum(0.045);
    hFrame->Draw();
    TGraphErrors *gr = new TGraphErrors(ntheta,theta,V,0,deltaV);
    gr->SetMarkerSize(1.2);
    gr->SetMarkerStyle(20);
    gr->Draw("Psame");
    TLine *l = new TLine(0,inV2,1.4,inV2);
    l->SetLineStyle(2);
    l->Draw("same");
    c1->Print("V2vstheta.png");
}

	NBD(TString datastr, TString Glauberstr,TString histonamev){
		dataname.SetString(datastr.Data());
		Glaubername.SetString(Glauberstr.Data());
		histoname.SetString(histonamev.Data());
	        mubest.ResizeTo(1);  kbest.ResizeTo(1);
	        chis.ResizeTo(1);  Ndf.ResizeTo(1);
	};