draw_matrix()
{

	TH2D *par_histo[44];
	TGraph2D *par_graph[44];

	for(int i=0; i<44; i++)
	{
		TString phname = "par_histo_"; phname+=i;
//		par_histo[i] = new TH2D(phname, phname, 10, 1200, 1300, 10, 0, TMath::TwoPi());
//		TString gname = "par";
//		gname+=i;
//		par_graph[i]=new TGraph2D(gname, gname);
		par_graph[i]=new TGraph2D;
		par_histo[i] = new TH2D(phname, phname, 200, 0, TMath::TwoPi(), 200, 0, 1400);
	}

	// Open file with parameter values
	Int_t pnum, ccdx, ccdy;
	// a_shift, r_shift - slight shift of points, to ensure that the star of that coordinates is inside the interpolation range
	Double_t pval, a_shift=0, r_shift=0;

	FILE *fpm = fopen("params_matrix", "r");
	while(fscanf(fpm, "%d\t%d\t%d\t%lf", &pnum, &ccdx, &ccdy, &pval)!=EOF)
	{
		Double_t axis = TMath::ATan2((ccdy-1031),(ccdx-1024))-TMath::Pi()/2;
		if(axis<0) axis+=TMath::TwoPi();
		Double_t r = sqrt((ccdx-1024)*(ccdx-1024)+(ccdy-1031)*(ccdy-1031));
//		par_histo[pnum]->SetBinContent(par_histo[pnum]->FindBin(r, axis), pval);

		par_graph[pnum]->SetPoint(par_graph[pnum]->GetN(), r+r_shift, axis+a_shift, pval);
		par_graph[pnum]->SetPoint(par_graph[pnum]->GetN(), r+r_shift, (axis+TMath::TwoPi())+a_shift, pval);
		par_graph[pnum]->SetPoint(par_graph[pnum]->GetN(), r+r_shift, (axis-TMath::TwoPi())+a_shift, pval);
		
		// Extend points for smalles and largest r
		// Correct extrapolation is only visible on unzoomed range of the graph!
		// choose points from largest r border
		
		// just three points with easy cutoff in current parameters matrix
		if(r>1290)
		{
			par_graph[pnum]->SetPoint(par_graph[pnum]->GetN(), 1500, axis+a_shift, pval);
			par_graph[pnum]->SetPoint(par_graph[pnum]->GetN(), 1500, (axis+TMath::TwoPi())+a_shift, pval);
			par_graph[pnum]->SetPoint(par_graph[pnum]->GetN(), 1500, (axis-TMath::TwoPi())+a_shift, pval);
		}
		
		// chose points from lowest r border
		if(r<60)
		{
			par_graph[pnum]->SetPoint(par_graph[pnum]->GetN(), -10, axis+a_shift, pval);
			par_graph[pnum]->SetPoint(par_graph[pnum]->GetN(), -10, (axis+TMath::TwoPi())+a_shift, pval);
			par_graph[pnum]->SetPoint(par_graph[pnum]->GetN(), -10, (axis-TMath::TwoPi())+a_shift, pval);
		}
		
//		par_graph[pnum]->SetPoint(par_graph[pnum]->GetN(), ccdx, ccdy, pval);

		// Manually set shifts to increase parameter space, not decrease it (taken from points graph)
/*		if(r<1215 && axis<1) { r_shift=-2; a_shift=-0.001;}
		else if(r<1190 && axis>5) { r_shift=-2; a_shift=0.001;}
		else if(r>1240) { r_shift=2; a_shift=0;}
		else if(r>1210 && axis>3.5) { r_shift=2; a_shift=0.001;}		
		par_graph[pnum]->SetPoint(par_graph[pnum]->GetN(), r+r_shift, axis+a_shift, pval); */
//		par_graph[pnum]->SetPoint(par_graph[pnum]->GetN(), ccdx, ccdy, pval);
//		cout << ccdx << " " << ccdy << " " << r << " " << axis << " " << pval << endl;
	}
		
	gStyle->SetPalette(1,0);

	TCanvas *c1[44];

//	for(int i=0; i<44; i++)
	for(int i=8; i<9; i++)
	{	
		if(par_graph[i]->GetN()>0)
		{
		
			// Fill histograms with interpolated values
			
			for(int x=0; x<par_histo[i]->GetNbinsX(); x++) for(int y=0; y<par_histo[i]->GetNbinsY(); y++)
			{
				Double_t grx=1400./par_histo[i]->GetNbinsX()*x;
				Double_t gry = TMath::TwoPi()/par_histo[i]->GetNbinsY()*y;
				Double_t grval = par_graph[i]->Interpolate(grx, gry);
				if(grval==0) cout << grx << " " << gry << endl;
				par_histo[i]->Fill(gry, grx, grval);
//				cout << grx << " " << gry << " " << grval << endl;
			}
		
			TString cname = "par";
			cname+=i;
			c1[i] = new TCanvas(cname, cname, 600, 500);
			c1[i]->SetTheta(90);
			c1[i]->SetPhi(-89.999);			
//			par_graph[i]->SetMarkerStyle(20);			
			par_graph[i]->Draw("tri2 p0");
			par_graph[i]->Draw("pcolz same");
			par_graph[i]->SetMarkerStyle(20);
			par_graph[i]->GetYaxis()->SetTitle("Angle [rad]");
			par_graph[i]->GetYaxis()->SetTitleOffset(1.2);
			par_graph[i]->GetXaxis()->SetTitleOffset(1.2);
			par_graph[i]->GetYaxis()->SetRangeUser(-1, (TMath::TwoPi()+1));
			par_graph[i]->GetXaxis()->SetRangeUser(-400, 1600);
//			gc->SetMarkerStyle(20);
			par_graph[i]->SetTitle(cname);
//			TGraph2D *gc = (TGraph2D*)par_graph[i]->DrawClone("pcolz same");
			
	
			c1[i]->Modified(); c1[i]->Update();
	
			new TCanvas();
			par_histo[i]->Draw("colz");
			par_histo[i]->SetName("par_histo");
//			par_histo[i]->SaveAs(TString("par_histo/")+cname+TString("_rad.root"));
			c1[i]->Modified(); c1[i]->Update();
cout << par_graph[i]->Interpolate(425, 1.91) << endl;	
			cname+="_rad.png";
//			c1[i]->SaveAs(TString("par_img/")+cname);
		}
	}
	
	cout << par_histo[pnum]->Interpolate(1350, 2);	

}