//experiments with fit. //also tried to draw the covariance ellipse: plot c1. //not finished but almost. still good enough to present(i think). #include #include #include #include "draw_ellipse.h" #include "TMath.h" #include "TH2F.h" using namespace std; TEllipse* draw_ellipse1(TMatrixD *means,TMatrixD *var, int x_dim,int y_dim){ //Temporary variance matrix double x_l,y_l; TMatrixD tempVar(2,2); tempVar(0,0)=var[0](x_dim,x_dim); tempVar(1,1)=var[0](y_dim,y_dim); tempVar(1,0)=var[0](x_dim,y_dim); tempVar(0,1)=var[0](x_dim,y_dim); TMatrixDEigen eigen(tempVar); TMatrixD eigen_vec=eigen.GetEigenVectors(); TMatrixD eigen_val=eigen.GetEigenValues(); double lamda_x,lamda_y; lamda_x=abs(eigen_val(0,0)); lamda_y=abs(eigen_val(1,1)); std::cout<SetStyle("Plain"); gStyle->SetOptFit(111111); gStyle->SetFrameBorderMode(0); gStyle->SetFillColor(0); gStyle->SetCanvasColor(0); gStyle->SetFrameFillColor(0); gStyle->SetCanvasBorderMode(0); // create a random number generator gRandom = new TRandom3(); // create a histogram TH1D * hist = new TH1D("data", ";x;N", 200, -10, 10); // fill in the histogram for (int i = 0; i < 1000; ++i) hist->Fill(gRandom->Gaus(0,1)); // define a fit function = gauss TF1 * f1 = new TF1("gauss", "[0] / sqrt(2.0 * TMath::Pi()) / [2] * exp(-(x-[1])*(x-[1])/2./[2]/[2])", 0, 100); //set parameter start values (mandatory). f1->SetParNames("Constant","Mean","Sigma"); f1->SetParameters(5.,hist->GetMean(),hist->GetRMS()); f1->SetParLimits(0, 10, 50); f1->SetParLimits(1, 0, 1); f1->SetParLimits(2, 0.1, 10.0); f1->SetLineWidth(2); f1->SetLineColor(2); // create a canvas to draw the histogram TCanvas * c1= new TCanvas("c1", "fitted data",5,5,800,600); // perform fit hist->Fit("gauss"); hist->Draw(); hist->SaveAs("fit.eps"); } void draw_ellipse_example(){ TMatrixD *sd=new TMatrixD(2,2); sd[0](0,0)=1; sd[0](1,1)=10; TMatrixD *u=new TMatrixD(2,1); TMatrixD *point=new TMatrixD(2,1); u[0](0,0)=1; u[0](1,0)=1; point[0](0,0)=2; point[0](1,0)=5; int dims=2; double denominator=pow(2*M_PI,dims)*sd[0].Determinant(); cout<Divide(3,3); // divides the canvas into three rows and three columns c1->cd(1); TGraph *h1 = new TGraph(); TGraph *h2 = new TGraph(); TGraph *h3 = new TGraph(); TGraph *h4 = new TGraph(); h1->SetMarkerStyle(23); h2->SetMarkerStyle(23); h1->SetMarkerColor(kRed); h2->SetMarkerColor(kBlue); double x1[2*30000]; double y1[2*30000]; double x2[30]; double y2[30]; int label[30*2]; TRandom randgen; for(int i=0;i<60000;i++){ x1[i]=randgen.Gaus(1,1); y1[i]=randgen.Gaus(1,10); } for(int i=0;i<60000;i++){ h1->SetPoint(h1->GetN(), x1[i],y1[i]); } h1->Draw("AP"); h1->SetMarkerStyle(2); h1->SetMarkerColor(kRed); //TGraph *g = new TGraph(2*30,x1,y1); //where x , y are 2 arrays of n points and draw this graph with option āpā //g->SetMarkerStyle(24); //g->Draw("ap"); TEllipse *el; cout<<"noe"<SetFillColorAlpha (1, 0); el->Draw(); double inx1=x1[0]; double inx2=x1[2]; double iny1=y1[0]; double iny2=y1[2]; int label1[2*30]; double max_x=TMath::MaxElement(60000,x1); double max_y=TMath::MaxElement(60000,y1); double min_x=TMath::MinElement(60000,x1); double min_y=TMath::MinElement(60000,y1); TCanvas *c2=new TCanvas(); TH2D *hist2=new TH2D("hist2", "noe", 50, min_x, max_x, 50, min_y, max_y); //TH2D *hist2=new TH2D("hist2"); hist2->FillN(60000,x1,y1,NULL,1); // hist2->Draw("colz"); //hist2->SetContour(5,contours); //hist2->Draw(); TCanvas *c7=new TCanvas(); TH1D *hist4=new TH1D("hist4", "noe1", 50, min_y, max_y); hist4->FillN(60000,y1,0); //hist4->Draw("colz"); //TF1 *f1 = new TF1("f1","gaus",1,1,1); //f1->SetParameters(hist4->GetMaximum(), hist4->GetMean(), hist4->GetRMS() ); //TF1 * f1 = new TF1("ga", "[0] / sqrt(2.0 * TMath::Pi()) / [2] * exp(-(x-[1])*(x-[1])/2./[2]/[2])", -10, 10); TF1 * f1 = new TF1("f1","gaus",-10,10); //set parameter start values (mandatory). f1->SetParNames("Constant","Mean","Sigma"); f1->SetParameters(4000,hist4->GetMean(),hist4->GetRMS()); f1->SetParLimits(0, 1, 1000000); f1->SetParLimits(1, -10, 10); f1->SetParLimits(2, 0.1, 10.0); f1->SetLineWidth(2); f1->SetLineColor(2); //f1->Print(); int fitStatus=hist4->Fit("f1"); hist4->Draw(); cout<GetParameters()<Draw(); //hist4->Draw("colz"); //TF1 *f1 = new TF1("f1","gaus",1,1,1); //f1->SetParameters(hist4->GetMaximum(), hist4->GetMean(), hist4->GetRMS() ); gRandom = new TRandom3(); TCanvas *c8=new TCanvas(); // create a histogram TH1D * hist5 = new TH1D("data", ";x;N", 200, -30, 30); // fill in the histogram for (int i = 0; i < 1000000; ++i){ if(i%3==0){ hist5->Fill(gRandom->Gaus(5,1));} else if(i%3==1){ hist5->Fill(gRandom->Gaus(-5,5));} else{//hist5->Fill(gRandom->Gaus(10,2)); } } //TString gauss= TF1 * f2 = new TF1("f2", "[0] / sqrt(2.0 * TMath::Pi()) / [2] * exp(-(x-[1])*(x-[1])/2./[2]/[2])+[3] / sqrt(2.0 * TMath::Pi()) / [5] * exp(-(x-[4])*(x-[4])/2./[5]/[5])",-30,30); //+[6] / sqrt(2.0 * TMath::Pi()) / [8] * exp(-(x-[7])*(x-[7])/2./[8]/[8]))", -30, 30); //TF1 * f2 = new TF1("f2","gaus(1)+gaus(1)",-30,30); //set parameter start values (mandatory). //f2->SetParNames("Constant","Mean","Sigma"); f2->SetParameters(4000,hist5->GetMean(),hist5->GetRMS(),4000,hist5->GetMean(),hist5->GetRMS());//,4000,hist5->GetMean(),hist5->GetRMS()); f2->SetParLimits(0, 1, 1000000); f2->SetParLimits(1, -20, 20); f2->SetParLimits(2, 0.1, 0.0); f2->SetParLimits(3, 1, 1000000); f2->SetParLimits(4, -20, 20); f2->SetParLimits(5, 0.1, 10.0); //f2->SetParLimits(6, 1, 1000000); //f2->SetParLimits(7, -20, 20); //f2->SetParLimits(8, 0.1, 10.0); f2->SetLineWidth(2); f2->SetLineColor(2); //f1->Print(); int fitStatus2=hist5->Fit("f2"); hist5->Draw(); //cout<GetParameters()<Draw(); //gauss-2D: //TF1 * f2 = new TF1("gauss2", "[0] / sqrt(2.0 * TMath::Pi()) / [2] * exp(-(x-[1])*(x-[1])/2./[2]/[2])+[3] / sqrt(2.0 * TMath::Pi()) / [5] * exp(-(x-[4])*(x-[4])/2./[5]/[5])",-30,30); //+[6] / sqrt(2.0 * TMath::Pi()) / [8] * exp(-(x-[7])*(x-[7])/2./[8]/[8]))", -30, 30); }