Problem saving canvas with multiple TGraphs in pdf- or eps-format

Root_for_Optics · February 8, 2022, 11:26am

Hi,
I have a plot consisting of 2 TGraphs each on 2 distinct TPads (one of them transparent), so that I can use 2 different axes (TGaxis). While I can save the canvas as png- or jpg-files (with correspondingly low quality) just fine, the eps-format is missing one of the four graphs, and the pdf-format is devoid of any graphs (and has only the standard, left-side axis). Has anyone experienced such a problem before? Solving the eps-format problem would be sufficient for me.

Part of the code:

	TGraph *g_gdd1 = new TGraph(x.size(),&(x[0]),&(y_gdd1[0]));
	TGraph *g_tod1 = new TGraph(x_tod1.size(),&(x_tod1[0]),&(y_tod1[0]));
	TGraph *g_gdd2 = new TGraph(x.size(),&(x[0]),&(y_gdd2[0]));
	TGraph *g_tod2 = new TGraph(x_tod2.size(),&(x_tod2[0]),&(y_tod2[0]));

	// plot:
	TCanvas *canvas = getCanvas();	
	int gdd_color = kBlue;
	int tod_color = kRed;
	double dx = (aoi2-aoi1)/0.8; // 10 per cent margins left and right
	double dy_gdd = (gdd_max-gdd_min)/0.8; // 10 per cent margins left and right
	double dy_tod = (tod_max-tod_min)/0.8; // 10 per cent margins left and right
	TPad *pad1 = new TPad("pad1","",0,0,1,1);
	pad1->Range(aoi1,gdd_min-0.1*dy_gdd,aoi2,gdd_max+0.1*dy_gdd);
	TPad *pad2 = new TPad("pad2","",0,0,1,1);
	pad2->Range(aoi1-0.1*dx,tod_min-0.1*dy_tod,aoi2+0.1*dx,tod_max+0.1*dy_tod);
	pad2->SetFillStyle(4000); //will be transparent
	
	pad1->Draw();
	pad1->cd();
	g_gdd1->Draw();
	g_gdd1->SetLineColor(gdd_color);
	g_gdd1->GetXaxis()->SetRangeUser(aoi1,aoi2);
	g_gdd1->GetXaxis()->SetTitle("AOI [#circ]");
	g_gdd1->GetXaxis()->CenterTitle();
	g_gdd1->GetYaxis()->SetRangeUser(gdd_min,gdd_max);
	g_gdd1->GetYaxis()->SetAxisColor(gdd_color);
	g_gdd1->GetYaxis()->SetLabelColor(gdd_color);
	g_gdd1->GetYaxis()->SetTitle("|GDD| [fs^{2}]");
	g_gdd1->GetYaxis()->SetTitleSize(0.04);
	g_gdd1->GetYaxis()->SetTitleColor(gdd_color);
	g_gdd1->GetYaxis()->CenterTitle();
	string title = "Simultaneous matching of GDD and TOD for #lambda = " + roundToNDigits(lambda,3) + "#mum";
	g_gdd1->SetTitle(title.c_str());
	g_gdd1->Draw();
	canvas->Update();
	g_gdd2->Draw("same");
	g_gdd2->SetLineColor(gdd_color);
	g_gdd1->SetLineWidth(2);
	g_gdd2->SetLineWidth(2);
	g_gdd2->Draw("same");
	g_gdd1->Draw("same");
	canvas->Update();
	
	TLine *l_gdd = new TLine(aoi1,g_gdd2->Eval(aoi_match2),max(aoi_match1,aoi_match2),g_gdd1->Eval(aoi_match1));
	l_gdd->SetLineStyle(2);
	l_gdd->SetLineColor(gdd_color);
	l_gdd->Draw();
	pad1->Update();
	
	pad2->Draw();
	pad2->cd();	
	g_tod2->Draw("same");	
	g_tod1->Draw("same");	
	g_tod1->GetYaxis()->SetRangeUser(tod_min,tod_max);
	g_tod2->GetYaxis()->SetRangeUser(tod_min,tod_max);
	g_tod1->SetLineColor(tod_color);
	g_tod2->SetLineColor(tod_color);
	g_tod1->SetLineWidth(2);
	g_tod2->SetLineWidth(2);
	//g_tod1->Draw("same");
	//g_tod2->Draw("same");
	pad2->Update();
	
	//double hack_factor = 1.00035;
    //TGaxis *axis = new TGaxis(aoi2*hack_factor,tod_min,aoi2*hack_factor,tod_max,tod_min,tod_max,50510,"+L");
    TGaxis *axis = new TGaxis(aoi2,tod_min,aoi2,tod_max,tod_min,tod_max,50510,"+L");
	axis->SetLabelColor(tod_color);
	axis->SetLabelSize(0.035);
	axis->SetLabelFont(42);
	axis->SetLineColor(tod_color);
	axis->SetTitle("TOD [fs^{3}]");
	axis->SetTitleSize(0.04);
	axis->SetTitleFont(42);
	axis->SetTitleColor(tod_color);
	axis->CenterTitle();
	axis->Draw();
	//pad2->Update();

	TLine *l_match1 = new TLine(aoi_match1,tod_min,aoi_match1,tod_min+0.6*(tod_max-tod_min));
	l_match1->SetLineStyle(2);
	l_match1->Draw();
	TLine *l_match2 = new TLine(aoi_match2,tod_min,aoi_match2,tod_min+0.6*(tod_max-tod_min));
	l_match2->SetLineStyle(2);
	l_match2->Draw();
	TLine *l_tod = new TLine(min(aoi_match1,aoi_match2),g_tod1->Eval(aoi_match1),aoi2,g_tod1->Eval(aoi_match1));
	l_tod->SetLineStyle(2);
	l_tod->SetLineColor(tod_color);
	l_tod->Draw();
	//pad2->Update();

	TMarker *tm1 = new TMarker();
	tm1->SetNDC(true);
	tm1->SetMarkerStyle(2);
	tm1->SetMarkerSize(3);
	tm1->DrawMarker(littrow1,tod1_littrow);
	TMarker *tm2 = new TMarker();
	tm2->SetNDC(true);
	tm2->SetMarkerStyle(2);
	tm2->SetMarkerSize(3);
	tm2->DrawMarker(littrow2,tod2_littrow);

	bool p1Right = (period1<period2);
	TPaveText *pt1 = new TPaveText((p1Right)?0.775:0.125,0.6,(p1Right)?0.875:0.225,0.75,"NDC");
	pt1->SetFillColor(0);
	pt1->SetFillStyle(4000);
	pt1->SetTextSize(0.025);
	pt1->AddText(getLineDensity(period1).c_str());
	string aoi_match1_str = "AOI = " + roundToNDigits(aoi_match1,1) + " #circ";
	pt1->AddText(aoi_match1_str.c_str());
	string d1_str = "d = " + roundToNDigits(d1,1) + " mm";
	pt1->AddText(d1_str.c_str());
	pt1->Draw("same");

	TPaveText *pt2 = new TPaveText((p1Right)?0.125:0.775,0.6,(p1Right)?0.225:0.875,0.75,"NDC");
	pt2->SetFillColor(0);
	pt2->SetFillStyle(4000);
	pt2->SetTextSize(0.025);
	pt2->AddText(getLineDensity(period2).c_str());
	string aoi_match2_str = "AOI = " + roundToNDigits(aoi_match2,1) + " #circ";
	pt2->AddText(aoi_match2_str.c_str());
	string d2_str = "d = " + roundToNDigits(d2,1) + " mm";
	pt2->AddText(d2_str.c_str());
	pt2->Draw("same");
	pad2->Update();
	g_tod1->Draw("same");	
	pad2->Update();

ROOT Version: 6.22/06
Platform: win32
Compiler: Not Provided

couet · February 8, 2022, 1:03pm

By change do you have some complete code we can run reproducing the problem ?

Root_for_Optics · February 8, 2022, 1:14pm

Not really. I mean I have 3 files that I could attach or upload; they run standalone.

Root_for_Optics · February 8, 2022, 1:25pm

Should I upload those 3 files?

couet · February 8, 2022, 1:27pm

I need a reproducer. to see the problem.

Root_for_Optics · February 8, 2022, 1:50pm

It says as a new user I can’t put links in the post …

couet · February 8, 2022, 1:55pm

The macro you posted earlier cannot be run because there is no data to build the graphs. If you modify just a bit to make it “runable” and if the problem can be reproducer with this little example, that’s fine.

Root_for_Optics · February 8, 2022, 2:38pm

Ok, I significantly shortened my macro and included some library functions, for me it works standalone now:

const double AOI_MIN = 5.0; // initialize minimum AOI to 5° for graph
const double AOI_MAX = 85.0; // initialize maximum AOI to 85° for graph and for calculation
const int LIGHTSPEED = 299792458; // c in vacuum in m/s
const int NUM_X1 = 500; // number of sampling points for graphs (first independent variable)
const int NUM_X2 = 100; // number of sampling points for graphs (second independent variable)
const double DUMMY_LAMBDA_MIN = -1;
const double DUMMY_LAMBDA_MAX = 123456789;
const double PRECISION = 0.00001; // 0.01 nm precision required for iterative calculation of lambda
const double LITTROW_MAX = 75.0; // set maximum Littrow to 75° for plotting (to avoid divergence)
const int DO_MAX = 10; // maximum diffraction order

enum MEDIUM{AIR,SUBSTRATE};

int canvas_counter = 1;

enum DISPERSION_ORDER{GDD,TOD,TODoverGDD};

const double approx_refractive_index = 1.45; // n = 1.45 for Suprasil 3002 between 752.5nm...1400nm

TGraph *Suprasil3002; // refractive index for Suprasil 3002

void initializeSuprasil3002()
{
	// data from Heraeus datasheet
	cout << "initializing refractive index of Suprasil 3002 ..." << endl;
	double x[] = {0.36548,0.380,0.400,0.40465,0.43583,0.4416,0.4471,0.48613,0.488,0.5145,0.532,0.54607,0.58756,0.6328,0.65627,0.6943,0.7525,0.800,0.850,0.900,0.905,1.000,1.064,1.153,1.200,1.319,1.400,1.600};
	double y[] = {1.4748,1.4728,1.4704,1.4699,1.4670,1.4665,1.4661,1.4634,1.4633,1.4619,1.4610,1.4604,1.4588,1.4573,1.4567,1.4557,1.4545,1.4536,1.4528,1.4520,1.4520,1.4507,1.4499,1.4489,1.4483,1.4470,1.4461,1.4437};
	Suprasil3002 = new TGraph(28,x,y);
	Suprasil3002->Draw();		
}

TCanvas* getCanvas()
{
	
	string canvasName = "c";
	canvasName += to_string(canvas_counter);
	//TCanvas *canvas = (TCanvas*)gROOT->GetListOfCanvases()->FindObject(canvasName);
	//TCanvas *canvas = new TCanvas(canvasName.c_str(),canvasName.c_str(),600,400);
	TCanvas *canvas = new TCanvas(canvasName.c_str(),canvasName.c_str(),1200,800);
	canvas_counter++;
	return canvas;
	
}

double getRefractiveIndex(double lambda)
{
	// lambda in mum
	if (Suprasil3002 == NULL) {
		initializeSuprasil3002();
	}
	//if ((lambda<0.365) || (lambda>1.6))
		//cout << "Achtung: value of refractive index of Suprasil 3002 at lambda = " << roundToNDigits(lambda,3) << " mum is not reliable!" << endl;
	
	//cout << Suprasil3002->Eval(lambda) << endl;
	//return 1.45;
	return Suprasil3002->Eval(lambda);	
}

string roundToNDigits(double x, int n, bool include_sign)
{
	double rounded = floor(x*TMath::Power(10,n)+0.5)/TMath::Power(10,n);
	string longString = to_string(rounded);
	string shortString = longString.substr(0,longString.find(".")+n+1);
	if (include_sign && (x>=0.0))
		shortString.insert(0,"+");
	return shortString;	
}

string roundToNDigits(double x, int n)
{
	return roundToNDigits(x,n,false);
}

string roundToNSignedDigits(double x, int n)
{
	return roundToNDigits(x,n,true);
}

string getLineDensity(int period) 
{
	// period in nm
	int num_digits = 1;
	double line_density = (double)TMath::Power(10,6)/period; // line density in 1/mm
	string text = roundToNDigits(line_density,num_digits) + " l/mm";
	return text;
}

double getLambdaMinWithoutHigherOrderDiffraction(int period)
{
	// minimum lambda before onset of higher order diffraction
	// lambda_min in mum, period in nm
	double DO_minusTwo_in_air = 2*(double)period/1000/3; // onset of diffraction order m = -2 (or +1) in air
	double initial_index = approx_refractive_index;
	double lambda = initial_index*DO_minusTwo_in_air;
	while (TMath::Abs(getRefractiveIndex(lambda)*DO_minusTwo_in_air - lambda) > PRECISION) {
		lambda = getRefractiveIndex(lambda)*DO_minusTwo_in_air;
		//cout << lambda << endl;
	}
	return lambda;
}

double getLambdaMaxAtLittrowMax(int period)
{
	// lambda_max in mum, period in nm
	double littrow_max = LITTROW_MAX/180.0*TMath::Pi();
	double max = 2*(double)period/1000*sin(littrow_max);
	return max;
}

double getLittrowAngle(int period, double lambda)
{
	// period in nm, lambda in mum
	const int m = -1; // diffraction order -1
	double lambda_nm = lambda*1000;
	double littrow = TMath::ASin(-m*lambda_nm/2/period);
	return littrow/TMath::Pi()*180.0; // littrow angle in degrees	
}

double aod_vs_aoi(double *var, double *par)
{
	// var[0] = external angle of incidence in degrees
	// par[0] = lambda in mum
	// par[1] = period in mum
	// par[2] = transmission order
	// par[3] = MEDIUM (AIR or SUBSTRATE)

	double input_angle = var[0]*TMath::Pi()/180.0; // convert angle from degrees to radians
	double lambda = par[0];
	int m = (int)par[2];
	MEDIUM medium = (MEDIUM)par[3];
	double index = (medium==SUBSTRATE) ? getRefractiveIndex(lambda) : 1;
	double output_angle = TMath::ASin((TMath::Sin(input_angle)+m*lambda/par[1])/index);
	return output_angle/TMath::Pi()*180.0; // aod in degress
}

double getDiffractionAngle(double aoi, double lambda, int period, int order = -1, MEDIUM medium = AIR)
{
	// aoi in degrees, lambda in mum, period in nm
	TF1 *f1 = new TF1("f_aod_vs_aoi",aod_vs_aoi,0,90,4);
	f1->SetParameter(0,lambda);
	f1->SetParameter(1,(double)period/1000);
	f1->SetParameter(2,(double)order);
	f1->SetParameter(3,(double)medium); 
	return f1->Eval(aoi);
}

double getIncidentAngle(double aod, double lambda, int period)
{
	// aod in degrees, lambda in mum, period in nm, diffraction order set to -1
	const int m = (aod<=0.0) ? -1 : 1; // sign convention for transission order
	return getDiffractionAngle(aod,lambda,period,-m,AIR); // AOI(AOD) = AOD(AOI) with m replaced by -m
}

double gdd_littrow(double *var, double *par)
{
	// var[0] = distance d in mm
	// var[1] = lambda in micrometer
	// par[0] = period in mum

	// Littrow cinfiguration
	const int m = -1; // transmission order
	const int passes = 2;
	double d = var[0]*1000; // convert d from mm to mum
	double lambda = var[1];
	const double v_c = LIGHTSPEED*TMath::Power(10,-9); // c in micrometer/femtosecond
	double f = -passes*2*d*lambda/TMath::Pi()/v_c/v_c*m*m*lambda*lambda/4/par[0]/par[0]*TMath::Power(1-m*m*lambda*lambda/4/par[0]/par[0],-1.5);
	return f;
	
}

double gdd(double *var, double *par)
{
	// var[0] = lambda in micrometer
	// var[1] = aoi in degrees
	// par[0] = period in mum
	// par[1] = distance d in mm

	const int m = -1; // transmission order
	const int passes = 2;
	double d = par[1]*1000; // convert d from mm to mum
	double lambda = var[0];
	double sin_aoi = TMath::Sin(var[1]*TMath::Pi()/180); // sin(aoi)
	const double v_c = LIGHTSPEED*TMath::Power(10,-9); // c in micrometer/femtosecond
	double f = -passes*d*lambda*lambda*lambda*m*m/2/TMath::Pi()/v_c/v_c/par[0]/par[0]*TMath::Power(1-(sin_aoi+m*lambda/par[0])*(sin_aoi+m*lambda/par[0]),-1.5);
	return f;
	
}

double tod_over_gdd_littrow(double *var, double *par)
{
	// var[0] = lambda in micrometer
	// par[0] = period in mum

	// Littrow cinfiguration
	const int m = -1; // transmission order
	double lambda = var[0];
	const double v_c = LIGHTSPEED*TMath::Power(10,-9); // c in micrometer/femtosecond
	double f = -3*lambda/2/TMath::Pi()/v_c*(1+m*m*lambda*lambda/4/par[0]/par[0])/(1-m*m*lambda*lambda/4/par[0]/par[0]);
	return f;
	
}

double tod_over_gdd(double *var, double *par)
{
	// var[0] = lambda in micrometer
	// var[1] = aoi in degrees
	// par[0] = period in mum

	const int m = -1; // transmission order
	double lambda = var[0];
	double sin_aoi = TMath::Sin(var[1]*TMath::Pi()/180); // sin(aoi)
	const double v_c = LIGHTSPEED*TMath::Power(10,-9); // c in micrometer/femtosecond
	double f = -3*lambda/2/TMath::Pi()/v_c*(1+m*lambda/par[0]*(sin_aoi+m*lambda/par[0])/(1-(sin_aoi+m*lambda/par[0])*(sin_aoi+m*lambda/par[0])));
	return f;
	
}

double calculateAOIforTODoverGDDmatching(double lambda, int period1, int period2, double aoi1)
{
	const int m = -1; // diffraction order = -1
	double p1 = (double)period1/1000; // period_1 in mum
	double p2 = (double)period2/1000; // period_2 in mum
	double sin_aoi1 = TMath::Sin(aoi1*TMath::Pi()/180); // sin(aoi_1)
	double sin_aod1 = sin_aoi1 + m*lambda/p1; // sin(aod_1)
	double b = 1-sin_aod1*sin_aod1; // expression in brackets
	double sin_aod2 = (TMath::Sqrt(p1*p1*b*b + 4*p2*p2*sin_aod1*sin_aod1) - p1*b)/2/p2/sin_aod1; // sin(aod_2) for negative aod
	double sin_aoi2 = sin_aod2 - m*lambda/p2; // sin(aoi_2)
	double aoi2 = TMath::ASin(sin_aoi2)*180/TMath::Pi(); // aoi_2 in degrees
	return aoi2;
}

double calculateAOIforGivenTODoverGDD(int period, double lambda, double gdd, double tod)
{
	const int m = -1; // diffraction order = -1
	double p = (double)period/1000; // period in mum
	const double v_c = LIGHTSPEED*TMath::Power(10,-9); // c in micrometer/femtosecond
	double a = -p/m/lambda*(1 + 2*TMath::Pi()*v_c/3/lambda*tod/gdd);
	double sin_aod = (TMath::Sqrt(1+4*a*a)-1)/2/a;
	double sin_aoi = sin_aod - m*lambda/p;
	double aoi = TMath::ASin(sin_aoi)*180/TMath::Pi();
	return aoi;
}

double calculateDistanceForGivenAOIandGDD(int period, double lambda, double gdd, double aoi_deg)
{
	const int m = -1; // diffraction order = -1
	double p = (double)period/1000; // period in mum
	const double v_c = LIGHTSPEED*TMath::Power(10,-9); // c in micrometer/femtosecond
	const int passes = 2;
	double aoi = aoi_deg/180*TMath::Pi();
	double sin_aod = TMath::Sin(aoi) + m*lambda/p;
	double d = -gdd*2*TMath::Pi()*v_c*v_c*p*p/passes/m/m/lambda/lambda/lambda*TMath::Power(1-sin_aod*sin_aod,1.5);
	return d/1000; // convert to mm
}

/////////////////////////////////////  plot macros  /////////////////////////////////////

void matchGDDandTOD(double lambda, int period1, int period2, double targetGDD, double targetTOD = DUMMY_LAMBDA_MIN)
{

	double lambda_min = getLambdaMinWithoutHigherOrderDiffraction(period1);
	double lambda_max = getLambdaMaxAtLittrowMax(period1);
	double littrow1 = getLittrowAngle(period1, lambda);
	double littrow2 = getLittrowAngle(period2, lambda);
	// find x-axis range:
	double x_margin = 0.1;
	// if targetTOD is not specified, match AODs (use arithmetic average as AOD for period1):
	double aoi_match1 = (targetTOD==DUMMY_LAMBDA_MIN) ? getIncidentAngle(-(littrow1+littrow2)/2,lambda,period1) : calculateAOIforGivenTODoverGDD(period1,lambda,targetGDD,targetTOD);
	double aoi_match2 = calculateAOIforTODoverGDDmatching(lambda,period1,period2,aoi_match1);
	double delta_aoi = TMath::Abs(aoi_match2-aoi_match1);
	double aoi1 = min(aoi_match1,aoi_match2) - x_margin*delta_aoi;
	double aoi2 = max(aoi_match1,aoi_match2) + x_margin*delta_aoi;
	cout << "AOI range: " << aoi1 << "..." << aoi2 << endl;
	double d1 = calculateDistanceForGivenAOIandGDD(period1,lambda,targetGDD,aoi_match1);
	double d2 = calculateDistanceForGivenAOIandGDD(period2,lambda,targetGDD,aoi_match2);
	cout << "GDD(TG" << period1 << "," << lambda << "mum)==GDD(TG" << period2 << "," << lambda << "mum) AND TOD(TG" << period1 << "," << lambda << "mum)==TOD(TG" << period2 << "," << lambda << "mum) for:" << endl;
	cout << "TG" << period1 << ": AOI = " << aoi_match1 << " deg, d = " << d1 << " mm" << endl;
	cout << "TG" << period2 << ": AOI = " << aoi_match2 << " deg, d = " << d2 << " mm" << endl;
	// find y-axes range:
	double y_margin = 0.15;
	TF2 *f2_gdd = new TF2("f_gdd",gdd,lambda_min,lambda_max,aoi1,aoi2,2);
	TF2 *f2_tog = new TF2("f_todovergdd",tod_over_gdd,lambda_min,lambda_max,aoi1,aoi2,1);
	f2_gdd->SetParameter(0,(double)period1/1000);
	f2_gdd->SetParameter(1,d1);
	f2_tog->SetParameter(0,(double)period1/1000);
	double gdd1_littrow = TMath::Abs(f2_gdd->Eval(lambda,littrow1));
	double tod1_littrow = f2_gdd->Eval(lambda,littrow1)*f2_tog->Eval(lambda,littrow1);
	f2_gdd->SetParameter(0,(double)period2/1000);
	f2_gdd->SetParameter(1,d2);
	f2_tog->SetParameter(0,(double)period2/1000);
	double gdd2_littrow = TMath::Abs(f2_gdd->Eval(lambda,littrow2));
	double tod2_littrow = f2_gdd->Eval(lambda,littrow2)*f2_tog->Eval(lambda,littrow2);
	double delta_gdd = TMath::Abs(gdd2_littrow-gdd1_littrow);
	double gdd_min = min(gdd1_littrow,gdd2_littrow) - y_margin*delta_gdd;
	double gdd_max = max(gdd1_littrow,gdd2_littrow) + y_margin*delta_gdd;
	cout << "GDD range: " << gdd_min << "..." << gdd_max << endl;
	double delta_tod = TMath::Abs(tod2_littrow-tod1_littrow);
	double tod_min = min(tod1_littrow,tod2_littrow) - y_margin*delta_tod;
	double tod_max = max(tod1_littrow,tod2_littrow) + y_margin*delta_tod;
	cout << "TOD range: " << tod_min << "..." << tod_max << endl;
	// make graphs:
	vector<double> x, x_tod1, x_tod2, y_gdd1, y_gdd2, y_tod1, y_tod2;
	for (int j=0; j<NUM_X1; j++) {
		double aoi = aoi1 + (aoi2-aoi1)/(NUM_X1-1)*j;
		f2_gdd->SetParameter(0,(double)period1/1000);
		f2_gdd->SetParameter(1,d1);
		double gdd1 = TMath::Abs(f2_gdd->Eval(lambda,aoi));
		f2_gdd->SetParameter(0,(double)period2/1000);
		f2_gdd->SetParameter(1,d2);
		double gdd2 = TMath::Abs(f2_gdd->Eval(lambda,aoi));
		x.push_back(aoi);
		y_gdd1.push_back(gdd1);	
		y_gdd2.push_back(gdd2);		
	}
	for (int j=0; j<NUM_X1; j++) {
		double aoi = aoi1 + (aoi2-aoi1)/(NUM_X1-1)*j;
		f2_gdd->SetParameter(0,(double)period1/1000);
		f2_gdd->SetParameter(1,d1);
		f2_tog->SetParameter(0,(double)period1/1000);
		double tod = f2_gdd->Eval(lambda,aoi)*f2_tog->Eval(lambda,aoi);
		if ((tod<tod_min) || (tod>tod_max))
			continue;
		x_tod1.push_back(aoi);
		y_tod1.push_back(tod);	
	}
	for (int j=0; j<NUM_X1; j++) {
		double aoi = aoi1 + (aoi2-aoi1)/(NUM_X1-1)*j;
		f2_gdd->SetParameter(0,(double)period2/1000);
		f2_gdd->SetParameter(1,d2);
		f2_tog->SetParameter(0,(double)period2/1000);
		double tod = f2_gdd->Eval(lambda,aoi)*f2_tog->Eval(lambda,aoi);
		if ((tod<tod_min) || (tod>tod_max))
			continue;
		x_tod2.push_back(aoi);
		y_tod2.push_back(tod);	
	}

	TGraph *g_gdd1 = new TGraph(x.size(),&(x[0]),&(y_gdd1[0]));
	TGraph *g_tod1 = new TGraph(x_tod1.size(),&(x_tod1[0]),&(y_tod1[0]));
	TGraph *g_gdd2 = new TGraph(x.size(),&(x[0]),&(y_gdd2[0]));
	TGraph *g_tod2 = new TGraph(x_tod2.size(),&(x_tod2[0]),&(y_tod2[0]));

	// plot:
	TCanvas *canvas = getCanvas();	
	int gdd_color = kBlue;
	int tod_color = kRed;
	double dx = (aoi2-aoi1)/0.8; // 10 per cent margins left and right
	double dy_gdd = (gdd_max-gdd_min)/0.8; // 10 per cent margins left and right
	double dy_tod = (tod_max-tod_min)/0.8; // 10 per cent margins left and right
	TPad *pad1 = new TPad("pad1","",0,0,1,1);
	pad1->Range(aoi1,gdd_min-0.1*dy_gdd,aoi2,gdd_max+0.1*dy_gdd);
	TPad *pad2 = new TPad("pad2","",0,0,1,1);
	pad2->Range(aoi1-0.1*dx,tod_min-0.1*dy_tod,aoi2+0.1*dx,tod_max+0.1*dy_tod);
	pad2->SetFillStyle(4000); //will be transparent
	
	pad1->Draw();
	pad1->cd();
	g_gdd1->Draw();
	g_gdd1->SetLineColor(gdd_color);
	g_gdd1->GetXaxis()->SetRangeUser(aoi1,aoi2);
	g_gdd1->GetXaxis()->SetTitle("AOI [#circ]");
	g_gdd1->GetXaxis()->CenterTitle();
	g_gdd1->GetYaxis()->SetRangeUser(gdd_min,gdd_max);
	g_gdd1->GetYaxis()->SetAxisColor(gdd_color);
	g_gdd1->GetYaxis()->SetLabelColor(gdd_color);
	g_gdd1->GetYaxis()->SetTitle("|GDD| [fs^{2}]");
	g_gdd1->GetYaxis()->SetTitleSize(0.04);
	g_gdd1->GetYaxis()->SetTitleColor(gdd_color);
	g_gdd1->GetYaxis()->CenterTitle();
	string title = "Simultaneous matching of GDD and TOD for #lambda = " + roundToNDigits(lambda,3) + "#mum";
	g_gdd1->SetTitle(title.c_str());
	g_gdd1->Draw();
	canvas->Update();
	g_gdd2->Draw("same");
	g_gdd2->SetLineColor(gdd_color);
	g_gdd1->SetLineWidth(2);
	g_gdd2->SetLineWidth(2);
	g_gdd2->Draw("same");
	g_gdd1->Draw("same");
	canvas->Update();
	
	TLine *l_gdd = new TLine(aoi1,g_gdd2->Eval(aoi_match2),max(aoi_match1,aoi_match2),g_gdd1->Eval(aoi_match1));
	l_gdd->SetLineStyle(2);
	l_gdd->SetLineColor(gdd_color);
	l_gdd->Draw();
	pad1->Update();
	
	pad2->Draw();
	pad2->cd();	
	g_tod2->Draw("same");	
	g_tod1->Draw("same");	
	g_tod1->GetYaxis()->SetRangeUser(tod_min,tod_max);
	g_tod2->GetYaxis()->SetRangeUser(tod_min,tod_max);
	g_tod1->SetLineColor(tod_color);
	g_tod2->SetLineColor(tod_color);
	g_tod1->SetLineWidth(2);
	g_tod2->SetLineWidth(2);
	//g_tod1->Draw("same");
	//g_tod2->Draw("same");
	pad2->Update();
	
	//double hack_factor = 1.00035;
    //TGaxis *axis = new TGaxis(aoi2*hack_factor,tod_min,aoi2*hack_factor,tod_max,tod_min,tod_max,50510,"+L");
    TGaxis *axis = new TGaxis(aoi2,tod_min,aoi2,tod_max,tod_min,tod_max,50510,"+L");
	axis->SetLabelColor(tod_color);
	axis->SetLabelSize(0.035);
	axis->SetLabelFont(42);
	axis->SetLineColor(tod_color);
	axis->SetTitle("TOD [fs^{3}]");
	axis->SetTitleSize(0.04);
	axis->SetTitleFont(42);
	axis->SetTitleColor(tod_color);
	axis->CenterTitle();
	axis->Draw();
	//pad2->Update();

	TLine *l_match1 = new TLine(aoi_match1,tod_min,aoi_match1,tod_min+0.6*(tod_max-tod_min));
	l_match1->SetLineStyle(2);
	l_match1->Draw();
	TLine *l_match2 = new TLine(aoi_match2,tod_min,aoi_match2,tod_min+0.6*(tod_max-tod_min));
	l_match2->SetLineStyle(2);
	l_match2->Draw();
	TLine *l_tod = new TLine(min(aoi_match1,aoi_match2),g_tod1->Eval(aoi_match1),aoi2,g_tod1->Eval(aoi_match1));
	l_tod->SetLineStyle(2);
	l_tod->SetLineColor(tod_color);
	l_tod->Draw();
	//pad2->Update();

	TMarker *tm1 = new TMarker();
	tm1->SetNDC(true);
	tm1->SetMarkerStyle(2);
	tm1->SetMarkerSize(3);
	tm1->DrawMarker(littrow1,tod1_littrow);
	TMarker *tm2 = new TMarker();
	tm2->SetNDC(true);
	tm2->SetMarkerStyle(2);
	tm2->SetMarkerSize(3);
	tm2->DrawMarker(littrow2,tod2_littrow);

	bool p1Right = (period1<period2);
	TPaveText *pt1 = new TPaveText((p1Right)?0.775:0.125,0.6,(p1Right)?0.875:0.225,0.75,"NDC");
	pt1->SetFillColor(0);
	pt1->SetFillStyle(4000);
	pt1->SetTextSize(0.025);
	pt1->AddText(getLineDensity(period1).c_str());
	string aoi_match1_str = "AOI = " + roundToNDigits(aoi_match1,1) + " #circ";
	pt1->AddText(aoi_match1_str.c_str());
	string d1_str = "d = " + roundToNDigits(d1,1) + " mm";
	pt1->AddText(d1_str.c_str());
	pt1->Draw("same");

	TPaveText *pt2 = new TPaveText((p1Right)?0.125:0.775,0.6,(p1Right)?0.225:0.875,0.75,"NDC");
	pt2->SetFillColor(0);
	pt2->SetFillStyle(4000);
	pt2->SetTextSize(0.025);
	pt2->AddText(getLineDensity(period2).c_str());
	string aoi_match2_str = "AOI = " + roundToNDigits(aoi_match2,1) + " #circ";
	pt2->AddText(aoi_match2_str.c_str());
	string d2_str = "d = " + roundToNDigits(d2,1) + " mm";
	pt2->AddText(d2_str.c_str());
	pt2->Draw("same");
	pad2->Update();
	g_tod1->Draw("same");	
	pad2->Update();
	
}

Root_for_Optics · February 8, 2022, 2:39pm

Ah, load it as a macro and then do e.g.:
matchGDDandTOD(1.03,543,575,-100000)

couet · February 8, 2022, 3:46pm

I get this plot:

It sounds a bit suspicious as some plots are outside the axis frame…
When I save as pdf I get a better plot.

couet · February 8, 2022, 5:06pm

It seems your first graph is corrupted. I did:

      printf("%g\n",gdd1); // print the y value of the 1st graph
      y_gdd1.push_back(gdd1);

I get:

TOD range: 191365...6.10093e+06
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
nan
6.10331e+08
1.60608e+08
8.02346e+07
5.00493e+07
3.49918e+07
2.62283e+07
2.06049e+07
1.67442e+07
1.3959e+07
1.18718e+07
1.026e+07
8.98453e+06
7.95472e+06
7.109e+06
6.40435e+06
5.80985e+06
5.3028e+06
4.86617e+06
4.48697e+06
4.15514e+06
3.86278e+06
3.60361e+06
3.37257e+06
3.16556e+06
2.97921e+06
2.81074e+06
2.65784e+06
2.51856e+06
2.39125e+06
2.27453e+06
2.16718e+06
2.0682e+06
[...]

There is many “Nan” in your graph.
With Nan the result is unpredictable.

couet · February 8, 2022, 6:29pm

Here is a working version of your macro.

const double AOI_MIN = 5.0; // initialize minimum AOI to 5° for graph
const double AOI_MAX = 85.0; // initialize maximum AOI to 85° for graph and for calculation
const int LIGHTSPEED = 299792458; // c in vacuum in m/s
const int NUM_X1 = 500; // number of sampling points for graphs (first independent variable)
const int NUM_X2 = 100; // number of sampling points for graphs (second independent variable)
const double DUMMY_LAMBDA_MIN = -1;
const double DUMMY_LAMBDA_MAX = 123456789;
const double PRECISION = 0.00001; // 0.01 nm precision required for iterative calculation of lambda
const double LITTROW_MAX = 75.0; // set maximum Littrow to 75° for plotting (to avoid divergence)
const int DO_MAX = 10; // maximum diffraction order

enum MEDIUM{AIR,SUBSTRATE};

int canvas_counter = 1;

enum DISPERSION_ORDER{GDD,TOD,TODoverGDD};

const double approx_refractive_index = 1.45; // n = 1.45 for Suprasil 3002 between 752.5nm...1400nm

TGraph *Suprasil3002; // refractive index for Suprasil 3002

void initializeSuprasil3002()
{
   // data from Heraeus datasheet
   cout << "initializing refractive index of Suprasil 3002 ..." << endl;
   double x[] = {0.36548,0.380,0.400,0.40465,0.43583,0.4416,0.4471,0.48613,0.488,0.5145,0.532,0.54607,0.58756,0.6328,0.65627,0.6943,0.7525,0.800,0.850,0.900,0.905,1.000,1.064,1.153,1.200,1.319,1.400,1.600};
   double y[] = {1.4748,1.4728,1.4704,1.4699,1.4670,1.4665,1.4661,1.4634,1.4633,1.4619,1.4610,1.4604,1.4588,1.4573,1.4567,1.4557,1.4545,1.4536,1.4528,1.4520,1.4520,1.4507,1.4499,1.4489,1.4483,1.4470,1.4461,1.4437};
   Suprasil3002 = new TGraph(28,x,y);
   Suprasil3002->Draw();
}

TCanvas* getCanvas()
{

   string canvasName = "c";
   canvasName += to_string(canvas_counter);
   TCanvas *canvas = new TCanvas(canvasName.c_str(),canvasName.c_str(),1200,800);
   canvas_counter++;
   return canvas;

}

double getRefractiveIndex(double lambda)
{
   // lambda in mum
   if (Suprasil3002 == NULL) {
      initializeSuprasil3002();
   }
   //if ((lambda<0.365) || (lambda>1.6))
      //cout << "Achtung: value of refractive index of Suprasil 3002 at lambda = " << roundToNDigits(lambda,3) << " mum is not reliable!" << endl;

   //cout << Suprasil3002->Eval(lambda) << endl;
   //return 1.45;
   return Suprasil3002->Eval(lambda);
}

string roundToNDigits(double x, int n, bool include_sign)
{
   double rounded = floor(x*TMath::Power(10,n)+0.5)/TMath::Power(10,n);
   string longString = to_string(rounded);
   string shortString = longString.substr(0,longString.find(".")+n+1);
   if (include_sign && (x>=0.0))
      shortString.insert(0,"+");
   return shortString;
}

string roundToNDigits(double x, int n)
{
   return roundToNDigits(x,n,false);
}

string roundToNSignedDigits(double x, int n)
{
   return roundToNDigits(x,n,true);
}

string getLineDensity(int period)
{
   // period in nm
   int num_digits = 1;
   double line_density = (double)TMath::Power(10,6)/period; // line density in 1/mm
   string text = roundToNDigits(line_density,num_digits) + " l/mm";
   return text;
}

double getLambdaMinWithoutHigherOrderDiffraction(int period)
{
   // minimum lambda before onset of higher order diffraction
   // lambda_min in mum, period in nm
   double DO_minusTwo_in_air = 2*(double)period/1000/3; // onset of diffraction order m = -2 (or +1) in air
   double initial_index = approx_refractive_index;
   double lambda = initial_index*DO_minusTwo_in_air;
   while (TMath::Abs(getRefractiveIndex(lambda)*DO_minusTwo_in_air - lambda) > PRECISION) {
      lambda = getRefractiveIndex(lambda)*DO_minusTwo_in_air;
      //cout << lambda << endl;
   }
   return lambda;
}

double getLambdaMaxAtLittrowMax(int period)
{
   // lambda_max in mum, period in nm
   double littrow_max = LITTROW_MAX/180.0*TMath::Pi();
   double max = 2*(double)period/1000*sin(littrow_max);
   return max;
}

double getLittrowAngle(int period, double lambda)
{
   // period in nm, lambda in mum
   const int m = -1; // diffraction order -1
   double lambda_nm = lambda*1000;
   double littrow = TMath::ASin(-m*lambda_nm/2/period);
   return littrow/TMath::Pi()*180.0; // littrow angle in degrees
}

double aod_vs_aoi(double *var, double *par)
{
   // var[0] = external angle of incidence in degrees
   // par[0] = lambda in mum
   // par[1] = period in mum
   // par[2] = transmission order
   // par[3] = MEDIUM (AIR or SUBSTRATE)

   double input_angle = var[0]*TMath::Pi()/180.0; // convert angle from degrees to radians
   double lambda = par[0];
   int m = (int)par[2];
   MEDIUM medium = (MEDIUM)par[3];
   double index = (medium==SUBSTRATE) ? getRefractiveIndex(lambda) : 1;
   double output_angle = TMath::ASin((TMath::Sin(input_angle)+m*lambda/par[1])/index);
   return output_angle/TMath::Pi()*180.0; // aod in degress
}

double getDiffractionAngle(double aoi, double lambda, int period, int order = -1, MEDIUM medium = AIR)
{
   // aoi in degrees, lambda in mum, period in nm
   TF1 *f1 = new TF1("f_aod_vs_aoi",aod_vs_aoi,0,90,4);
   f1->SetParameter(0,lambda);
   f1->SetParameter(1,(double)period/1000);
   f1->SetParameter(2,(double)order);
   f1->SetParameter(3,(double)medium);
   return f1->Eval(aoi);
}

double getIncidentAngle(double aod, double lambda, int period)
{
   // aod in degrees, lambda in mum, period in nm, diffraction order set to -1
   const int m = (aod<=0.0) ? -1 : 1; // sign convention for transission order
   return getDiffractionAngle(aod,lambda,period,-m,AIR); // AOI(AOD) = AOD(AOI) with m replaced by -m
}

double gdd_littrow(double *var, double *par)
{
   // var[0] = distance d in mm
   // var[1] = lambda in micrometer
   // par[0] = period in mum

   // Littrow cinfiguration
   const int m = -1; // transmission order
   const int passes = 2;
   double d = var[0]*1000; // convert d from mm to mum
   double lambda = var[1];
   const double v_c = LIGHTSPEED*TMath::Power(10,-9); // c in micrometer/femtosecond
   double f = -passes*2*d*lambda/TMath::Pi()/v_c/v_c*m*m*lambda*lambda/4/par[0]/par[0]*TMath::Power(1-m*m*lambda*lambda/4/par[0]/par[0],-1.5);
   return f;

}

double gdd(double *var, double *par)
{
   // var[0] = lambda in micrometer
   // var[1] = aoi in degrees
   // par[0] = period in mum
   // par[1] = distance d in mm

   const int m = -1; // transmission order
   const int passes = 2;
   double d = par[1]*1000; // convert d from mm to mum
   double lambda = var[0];
   double sin_aoi = TMath::Sin(var[1]*TMath::Pi()/180); // sin(aoi)
   const double v_c = LIGHTSPEED*TMath::Power(10,-9); // c in micrometer/femtosecond
   double f = -passes*d*lambda*lambda*lambda*m*m/2/TMath::Pi()/v_c/v_c/par[0]/par[0]*TMath::Power(1-(sin_aoi+m*lambda/par[0])*(sin_aoi+m*lambda/par[0]),-1.5);
   if(TMath::IsNaN(f)) {
      return 0;
   } else {
      return f;
   }
}

double tod_over_gdd_littrow(double *var, double *par)
{
   // var[0] = lambda in micrometer
   // par[0] = period in mum

   // Littrow cinfiguration
   const int m = -1; // transmission order
   double lambda = var[0];
   const double v_c = LIGHTSPEED*TMath::Power(10,-9); // c in micrometer/femtosecond
   double f = -3*lambda/2/TMath::Pi()/v_c*(1+m*m*lambda*lambda/4/par[0]/par[0])/(1-m*m*lambda*lambda/4/par[0]/par[0]);
   return f;

}

double tod_over_gdd(double *var, double *par)
{
   // var[0] = lambda in micrometer
   // var[1] = aoi in degrees
   // par[0] = period in mum

   const int m = -1; // transmission order
   double lambda = var[0];
   double sin_aoi = TMath::Sin(var[1]*TMath::Pi()/180); // sin(aoi)
   const double v_c = LIGHTSPEED*TMath::Power(10,-9); // c in micrometer/femtosecond
   double f = -3*lambda/2/TMath::Pi()/v_c*(1+m*lambda/par[0]*(sin_aoi+m*lambda/par[0])/(1-(sin_aoi+m*lambda/par[0])*(sin_aoi+m*lambda/par[0])));
   return f;

}

double calculateAOIforTODoverGDDmatching(double lambda, int period1, int period2, double aoi1)
{
   const int m = -1; // diffraction order = -1
   double p1 = (double)period1/1000; // period_1 in mum
   double p2 = (double)period2/1000; // period_2 in mum
   double sin_aoi1 = TMath::Sin(aoi1*TMath::Pi()/180); // sin(aoi_1)
   double sin_aod1 = sin_aoi1 + m*lambda/p1; // sin(aod_1)
   double b = 1-sin_aod1*sin_aod1; // expression in brackets
   double sin_aod2 = (TMath::Sqrt(p1*p1*b*b + 4*p2*p2*sin_aod1*sin_aod1) - p1*b)/2/p2/sin_aod1; // sin(aod_2) for negative aod
   double sin_aoi2 = sin_aod2 - m*lambda/p2; // sin(aoi_2)
   double aoi2 = TMath::ASin(sin_aoi2)*180/TMath::Pi(); // aoi_2 in degrees
   return aoi2;
}

double calculateAOIforGivenTODoverGDD(int period, double lambda, double gdd, double tod)
{
   const int m = -1; // diffraction order = -1
   double p = (double)period/1000; // period in mum
   const double v_c = LIGHTSPEED*TMath::Power(10,-9); // c in micrometer/femtosecond
   double a = -p/m/lambda*(1 + 2*TMath::Pi()*v_c/3/lambda*tod/gdd);
   double sin_aod = (TMath::Sqrt(1+4*a*a)-1)/2/a;
   double sin_aoi = sin_aod - m*lambda/p;
   double aoi = TMath::ASin(sin_aoi)*180/TMath::Pi();
   return aoi;
}

double calculateDistanceForGivenAOIandGDD(int period, double lambda, double gdd, double aoi_deg)
{
   const int m = -1; // diffraction order = -1
   double p = (double)period/1000; // period in mum
   const double v_c = LIGHTSPEED*TMath::Power(10,-9); // c in micrometer/femtosecond
   const int passes = 2;
   double aoi = aoi_deg/180*TMath::Pi();
   double sin_aod = TMath::Sin(aoi) + m*lambda/p;
   double d = -gdd*2*TMath::Pi()*v_c*v_c*p*p/passes/m/m/lambda/lambda/lambda*TMath::Power(1-sin_aod*sin_aod,1.5);
   return d/1000; // convert to mm
}

/////////////////////////////////////  plot macros  /////////////////////////////////////

void matchGDDandTOD(double lambda, int period1, int period2, double targetGDD, double targetTOD = DUMMY_LAMBDA_MIN)
{

   double lambda_min = getLambdaMinWithoutHigherOrderDiffraction(period1);
   double lambda_max = getLambdaMaxAtLittrowMax(period1);
   double littrow1 = getLittrowAngle(period1, lambda);
   double littrow2 = getLittrowAngle(period2, lambda);
   // find x-axis range:
   double x_margin = 0.1;
   // if targetTOD is not specified, match AODs (use arithmetic average as AOD for period1):
   double aoi_match1 = (targetTOD==DUMMY_LAMBDA_MIN) ? getIncidentAngle(-(littrow1+littrow2)/2,lambda,period1) : calculateAOIforGivenTODoverGDD(period1,lambda,targetGDD,targetTOD);
   double aoi_match2 = calculateAOIforTODoverGDDmatching(lambda,period1,period2,aoi_match1);
   double delta_aoi = TMath::Abs(aoi_match2-aoi_match1);
   double aoi1 = min(aoi_match1,aoi_match2) - x_margin*delta_aoi;
   double aoi2 = max(aoi_match1,aoi_match2) + x_margin*delta_aoi;
   cout << "AOI range: " << aoi1 << "..." << aoi2 << endl;
   double d1 = calculateDistanceForGivenAOIandGDD(period1,lambda,targetGDD,aoi_match1);
   double d2 = calculateDistanceForGivenAOIandGDD(period2,lambda,targetGDD,aoi_match2);
   cout << "GDD(TG" << period1 << "," << lambda << "mum)==GDD(TG" << period2 << "," << lambda << "mum) AND TOD(TG" << period1 << "," << lambda << "mum)==TOD(TG" << period2 << "," << lambda << "mum) for:" << endl;
   cout << "TG" << period1 << ": AOI = " << aoi_match1 << " deg, d = " << d1 << " mm" << endl;
   cout << "TG" << period2 << ": AOI = " << aoi_match2 << " deg, d = " << d2 << " mm" << endl;
   // find y-axes range:
   double y_margin = 0.15;
   TF2 *f2_gdd = new TF2("f_gdd",gdd,lambda_min,lambda_max,aoi1,aoi2,2);
   TF2 *f2_tog = new TF2("f_todovergdd",tod_over_gdd,lambda_min,lambda_max,aoi1,aoi2,1);
   f2_gdd->SetParameter(0,(double)period1/1000);
   f2_gdd->SetParameter(1,d1);
   f2_tog->SetParameter(0,(double)period1/1000);
   double gdd1_littrow = TMath::Abs(f2_gdd->Eval(lambda,littrow1));
   double tod1_littrow = f2_gdd->Eval(lambda,littrow1)*f2_tog->Eval(lambda,littrow1);
   f2_gdd->SetParameter(0,(double)period2/1000);
   f2_gdd->SetParameter(1,d2);
   f2_tog->SetParameter(0,(double)period2/1000);
   double gdd2_littrow = TMath::Abs(f2_gdd->Eval(lambda,littrow2));
   double tod2_littrow = f2_gdd->Eval(lambda,littrow2)*f2_tog->Eval(lambda,littrow2);
   double delta_gdd = TMath::Abs(gdd2_littrow-gdd1_littrow);
   double gdd_min = min(gdd1_littrow,gdd2_littrow) - y_margin*delta_gdd;
   double gdd_max = max(gdd1_littrow,gdd2_littrow) + y_margin*delta_gdd;
   cout << "GDD range: " << gdd_min << "..." << gdd_max << endl;
   double delta_tod = TMath::Abs(tod2_littrow-tod1_littrow);
   double tod_min = min(tod1_littrow,tod2_littrow) - y_margin*delta_tod;
   double tod_max = max(tod1_littrow,tod2_littrow) + y_margin*delta_tod;
   cout << "TOD range: " << tod_min << "..." << tod_max << endl;
   // make graphs:
   vector<double> x, x_tod1, x_tod2, y_gdd1, y_gdd2, y_tod1, y_tod2;
   for (int j=0; j<NUM_X1; j++) {
      double aoi = aoi1 + (aoi2-aoi1)/(NUM_X1-1)*j;
      f2_gdd->SetParameter(0,(double)period1/1000);
      f2_gdd->SetParameter(1,d1);
      double gdd1 = TMath::Abs(f2_gdd->Eval(lambda,aoi));
      f2_gdd->SetParameter(0,(double)period2/1000);
      f2_gdd->SetParameter(1,d2);
      double gdd2 = TMath::Abs(f2_gdd->Eval(lambda,aoi));
      x.push_back(aoi);
      y_gdd1.push_back(gdd1);
      y_gdd2.push_back(gdd2);
   }
   for (int j=0; j<NUM_X1; j++) {
      double aoi = aoi1 + (aoi2-aoi1)/(NUM_X1-1)*j;
      f2_gdd->SetParameter(0,(double)period1/1000);
      f2_gdd->SetParameter(1,d1);
      f2_tog->SetParameter(0,(double)period1/1000);
      double tod = f2_gdd->Eval(lambda,aoi)*f2_tog->Eval(lambda,aoi);
      if ((tod<tod_min) || (tod>tod_max))
         continue;
      x_tod1.push_back(aoi);
      y_tod1.push_back(tod);
   }
   for (int j=0; j<NUM_X1; j++) {
      double aoi = aoi1 + (aoi2-aoi1)/(NUM_X1-1)*j;
      f2_gdd->SetParameter(0,(double)period2/1000);
      f2_gdd->SetParameter(1,d2);
      f2_tog->SetParameter(0,(double)period2/1000);
      double tod = f2_gdd->Eval(lambda,aoi)*f2_tog->Eval(lambda,aoi);
      if ((tod<tod_min) || (tod>tod_max))
         continue;
      x_tod2.push_back(aoi);
      y_tod2.push_back(tod);
   }

   TGraph *g_gdd1 = new TGraph(x.size(),&(x[0]),&(y_gdd1[0]));
   TGraph *g_tod1 = new TGraph(x_tod1.size(),&(x_tod1[0]),&(y_tod1[0]));
   TGraph *g_gdd2 = new TGraph(x.size(),&(x[0]),&(y_gdd2[0]));
   TGraph *g_tod2 = new TGraph(x_tod2.size(),&(x_tod2[0]),&(y_tod2[0]));

   // plot:
   TCanvas *canvas = getCanvas();
   int gdd_color = kBlue;
   int tod_color = kRed;
   double dx = (aoi2-aoi1)/0.8; // 10 per cent margins left and right
   double dy_gdd = (gdd_max-gdd_min)/0.8; // 10 per cent margins left and right
   double dy_tod = (tod_max-tod_min)/0.8; // 10 per cent margins left and right

   g_gdd1->SetLineColor(gdd_color);
   g_gdd1->SetLineWidth(2);
   g_gdd1->GetXaxis()->SetRangeUser(aoi1,aoi2);
   g_gdd1->GetXaxis()->SetTitle("AOI [#circ]");
   g_gdd1->GetXaxis()->CenterTitle();
   g_gdd1->GetYaxis()->SetRangeUser(gdd_min,gdd_max);
   g_gdd1->GetYaxis()->SetAxisColor(gdd_color);
   g_gdd1->GetYaxis()->SetLabelColor(gdd_color);
   g_gdd1->GetYaxis()->SetTitle("|GDD| [fs^{2}]");
   g_gdd1->GetYaxis()->SetTitleSize(0.04);
   g_gdd1->GetYaxis()->SetTitleColor(gdd_color);
   g_gdd1->GetYaxis()->CenterTitle();
   string title = "Simultaneous matching of GDD and TOD for #lambda = " + roundToNDigits(lambda,3) + "#mum";
   g_gdd1->SetTitle(title.c_str());
   g_gdd1->Draw("AL");
   g_gdd2->SetLineColor(gdd_color);
   g_gdd2->SetLineWidth(2);
   g_gdd2->Draw("L");

   TLine *l_gdd = new TLine(aoi1,g_gdd2->Eval(aoi_match2),max(aoi_match1,aoi_match2),g_gdd1->Eval(aoi_match1));
   l_gdd->SetLineStyle(2);
   l_gdd->SetLineColor(gdd_color);
   l_gdd->Draw();

   TPad *pad2 = new TPad("pad2","",0,0,1,1);
   pad2->SetFillColor(0);
   pad2->SetFillStyle(4000);
   pad2->SetFrameFillStyle(0);
   pad2->Draw();
   pad2->cd();

   g_tod2->SetTitle("");
   g_tod2->Draw("AL");
   g_tod2->GetYaxis()->SetRangeUser(tod_min,tod_max);
   g_tod2->GetXaxis()->SetRangeUser(aoi1,aoi2);
   g_tod2->GetYaxis()->SetLabelSize(0);
   g_tod2->GetYaxis()->SetTickSize(0);

   pad2->Update();
   g_tod1->Draw("L");
   g_tod1->SetLineColor(tod_color);
   g_tod2->SetLineColor(tod_color);
   g_tod1->SetLineWidth(2);
   g_tod2->SetLineWidth(2);
   pad2->Update();

   TGaxis *axis = new TGaxis(aoi2,tod_min,aoi2,tod_max,tod_min,tod_max,50510,"+L");
   axis->SetLabelColor(tod_color);
   axis->SetLabelSize(0.035);
   axis->SetLabelFont(42);
   axis->SetLineColor(tod_color);
   axis->SetTitle("TOD [fs^{3}]");
   axis->SetTitleSize(0.04);
   axis->SetTitleFont(42);
   axis->SetTitleColor(tod_color);
   axis->CenterTitle();
   axis->Draw();

   TLine *l_match1 = new TLine(aoi_match1,tod_min,aoi_match1,tod_min+0.6*(tod_max-tod_min));
   l_match1->SetLineStyle(2);
   l_match1->Draw();
   TLine *l_match2 = new TLine(aoi_match2,tod_min,aoi_match2,tod_min+0.6*(tod_max-tod_min));
   l_match2->SetLineStyle(2);
   l_match2->Draw();
   TLine *l_tod = new TLine(min(aoi_match1,aoi_match2),g_tod1->Eval(aoi_match1),aoi2,g_tod1->Eval(aoi_match1));
   l_tod->SetLineStyle(2);
   l_tod->SetLineColor(tod_color);
   l_tod->Draw();

   TMarker *tm1 = new TMarker();
   tm1->SetNDC(true);
   tm1->SetMarkerStyle(2);
   tm1->SetMarkerSize(3);
   tm1->DrawMarker(littrow1,tod1_littrow);
   TMarker *tm2 = new TMarker();
   tm2->SetNDC(true);
   tm2->SetMarkerStyle(2);
   tm2->SetMarkerSize(3);
   tm2->DrawMarker(littrow2,tod2_littrow);

   bool p1Right = (period1<period2);
   TPaveText *pt1 = new TPaveText((p1Right)?0.775:0.125,0.6,(p1Right)?0.875:0.225,0.75,"NDC");
   pt1->SetFillColor(0);
   pt1->SetFillStyle(0);
   pt1->SetTextSize(0.025);
   pt1->AddText(getLineDensity(period1).c_str());
   string aoi_match1_str = "AOI = " + roundToNDigits(aoi_match1,1) + " #circ";
   pt1->AddText(aoi_match1_str.c_str());
   string d1_str = "d = " + roundToNDigits(d1,1) + " mm";
   pt1->AddText(d1_str.c_str());
   pt1->Draw();

   TPaveText *pt2 = new TPaveText((p1Right)?0.125:0.775,0.6,(p1Right)?0.225:0.875,0.75,"NDC");
   pt2->SetFillColor(0);
   pt2->SetFillStyle(0);
   pt2->SetTextSize(0.025);
   pt2->AddText(getLineDensity(period2).c_str());
   string aoi_match2_str = "AOI = " + roundToNDigits(aoi_match2,1) + " #circ";
   pt2->AddText(aoi_match2_str.c_str());
   string d2_str = "d = " + roundToNDigits(d2,1) + " mm";
   pt2->AddText(d2_str.c_str());
   pt2->Draw();
   pad2->Update();
}

I restructured a lot the plotting part (it is simpler). And also I added a protection against NaN : now when a Nan is found, 0 is returned. You surely have a better fix for that. I did not know what would be the best fix for that. Returning 0 is just to avoid NaN in TGraph. The vertical line in the middle of the plot is due to that. The y values suddenly change from 0 to something big outside the vertical range.

Note on the option “SAME”: in your macro you used the option “SAME” in many places.This option is only valid for histograms. For instance it is not valid for TGraph or TPaveStats.

Root_for_Optics · February 9, 2022, 8:23am

Dear Olivier,

thanks a lot for your detailed investigation. The plot you fixed looks almost exactly like mine (except for the vertical blue line). I have to look into why the graph is corrupted for you but not for me. I also need to check your code of course, thanks again for putting so much effort into it.
But did you try to save it as eps or pdf. What did it look like?

couet · February 9, 2022, 8:45am

Yes, see the following screen dump. left to right: screen, ps, pdf

Root_for_Optics · February 9, 2022, 8:54am

Ok, so couldn’t reproduce the problem. I have to thoroughly investigate the differences between your and my code (I can’t do it right now); hopefully there’s a hint or even the solution to the problem.
I will update you once I have more to say.
Thanks again.

couet · February 9, 2022, 8:58am

What do you get with my code ? are ps and pdf ok with my code ?

Root_for_Optics · February 9, 2022, 9:44am

Yes, both look good, halleluja!
Thanks so much.
As I said, I still need to investigate what’s the differnece exactly…

couet · February 9, 2022, 1:13pm

Mainly two:

The protection against NaN → up to you to decide what to do in case in NaN appears
The plotting part has been re-tructured completly:

No need of the “SAME” option with graphs. The first graph needs option “A” (see the doc for details.
No need for the pad “pad1”
Remove some useless “Draw”
…

system · February 23, 2022, 1:14pm

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.