ROOT: Visualization of conformal mapping of points with isochrone radii - including a comparison between a simplified approach and a correct one

confmaptest.cpp

#include <iostream>
#include <algorithm>
#include <vector>
#include <cfloat>

void readPoints(string filename, std::vector<double> &x, std::vector<double> &y, std::vector<double> &r, int upToLineNumber = 2) {
	ifstream file(filename.c_str());
	float tempX, tempY, tempZ, tempR, tempI;
	char tempChar[10];
	int i = 1;
	while(i <= upToLineNumber && file >> tempX >> tempY >> tempZ >> tempR >> tempI >> tempChar) {
		x.push_back(tempX);
		y.push_back(tempY);
		r.push_back(tempR);
		i++;
	}
	file.close();
}

bool epsilon_equal_to(double value, double valueToCompareTo) {
	return (value >= (valueToCompareTo - DBL_EPSILON) && value <= (valueToCompareTo + DBL_EPSILON) );
}

double function_confMap(double value, double x, double y) {
	double x2y2 = x*x + y*y;
	return value/x2y2;
}
double confMap_x(double x, double y) {
	return function_confMap(x, x, y);
}
double confMap_y(double x, double y) {
	return function_confMap(-y, x, y);
}
double confMap_r(double r, double x, double y) {
	return function_confMap(r, x, y);
}

TGraph * pointsToGraph(const std::vector<double> x, const std::vector<double> y) {
	TGraph * thisGraph = new TGraph();
	for (int i = 0; i < x.size(); ++i) {
		thisGraph->SetPoint(i, x[i], y[i]);
	}
	return thisGraph;
}

std::vector<double> vPointsOnCircle_x(const double x, const double r, const int everyXDegrees) {
	std::vector<double> vX;
	for (int i = 0; i < 360/everyXDegrees; ++i) {
		double circleX = x + r * cos(i*everyXDegrees*TMath::DegToRad());
		vX.push_back(circleX);
	}
	return vX;
}
std::vector<double> vPointsOnCircle_y(const double y, const double r, const int everyXDegrees) {
	std::vector<double> vY;
	for (int i = 0; i < 360/everyXDegrees; ++i) {
		double circleY = y + r * sin(i*everyXDegrees*TMath::DegToRad());
		vY.push_back(circleY);
	}
	return vY;
}
TGraph * circleAroundPoint(const double x, const double y, const double r, const int everyXDegrees) {
	return pointsToGraph(vPointsOnCircle_x(x, r, everyXDegrees), vPointsOnCircle_y(y, r, everyXDegrees));
}

TGraph * cfCircleAroundPoint(const double x, const double y, const double r, const int everyXDegrees) {
	std::vector<double> vX, vY;

	/**
	 * take every points on the original circle and conformal transform it according to its central (x,y) point
	 */
	vX = vPointsOnCircle_x(x, r, everyXDegrees);
	vY = vPointsOnCircle_y(y, r, everyXDegrees);

	std::vector<double> cfVX, cfVY;

	for (int i = 0; i < vX.size(); ++i) {
		// double cfX = confMap_r(vX[i], x, y);
		double cfX = confMap_r(vX[i], vX[i], vY[i]);
		// double cfY = confMap_r(-vY[i], x, y);
		double cfY = confMap_r(-vY[i], vX[i], vY[i]);
		cfVX.push_back(cfX);
		cfVY.push_back(cfY);
	}
	return pointsToGraph(cfVX, cfVY);
}

int confmaptest2() {

	int verbose = 1;
	std::vector<double> x, y, r;

	// readPoints("data_confmap.dat", x, y, r, 18);
	readPoints("tubePos_1GeV_ideal_radius.dat", x, y, r, 18);
	 
	std::vector<double> cfX, cfY, cfR;
	for (int i = 0; i < x.size(); ++i) {
		cfX.push_back(confMap_x(x[i], y[i]));
		cfY.push_back(confMap_y(x[i], y[i]));
		cfR.push_back(confMap_r(r[i], x[i], y[i])); // Yuties method
		if (verbose >= 1) std::cout << "cfX[" << i << "] = " << cfX[i] << ", cfY = " << cfY[i] << ", cfR = " << cfR[i] << std::endl;
	}

	TGraph * originalPoints = pointsToGraph(x, y);
	originalPoints->SetMarkerStyle(kFullDotMedium);
	originalPoints->SetMarkerColor(kBlue);

	TGraph * cfPoints = pointsToGraph(cfX, cfY);
	cfPoints->SetMarkerStyle(kFullDotMedium);
	cfPoints->SetMarkerColor(kRed);

	TMultiGraph * mgOriginalSpace = new TMultiGraph();
	mgOriginalSpace->Add(originalPoints, "P");

	TMultiGraph * mgCfSpace = new TMultiGraph();
	mgCfSpace->Add(cfPoints, "P");

	int everyXDegrees = 10;
	for (int i = 0; i < x.size(); ++i) {
		TGraph * circleAroundOriginalPointI = circleAroundPoint(x[i], y[i], r[i], everyXDegrees);
		circleAroundOriginalPointI->SetLineColor(originalPoints->GetMarkerColor()-3);
		mgOriginalSpace->Add(circleAroundOriginalPointI, "L");
	}
	// Yuties method of transforming r
	std::vector<TGraph*> theCfCircles;
	for (int i = 0; i < x.size(); ++i) {
		TGraph * circleAroundCfPointI = circleAroundPoint(cfX[i], cfY[i], cfR[i], everyXDegrees);
		circleAroundCfPointI->SetLineColor(cfPoints->GetMarkerColor()-3);
		theCfCircles.push_back(circleAroundCfPointI);
		mgCfSpace->Add(circleAroundCfPointI, "L");
	}
	// My method of transforming r
	std::vector<TGraph*> theMyCfCircles;
	for (int i = 0; i < x.size(); ++i) {
		TGraph * myCircleAroundCFPointI = cfCircleAroundPoint(x[i], y[i], r[i], everyXDegrees);
		myCircleAroundCFPointI->SetLineColor(kOrange-3);
		theMyCfCircles.push_back(myCircleAroundCFPointI);
		mgCfSpace->Add(myCircleAroundCFPointI, "L");
	}

	// Calculate differences
	int printPointNumber = 2; // if verbose is > 0, which point should be printed as debug info
	std::vector<std::vector<double> > theTheDeltaRadii;
	std::vector<std::vector<double> > theTheNormalizedDeltaRadii;
	for (int i = 0; i < x.size(); ++i) {
		std::vector<double> theDeltaRadii;
		std::vector<double> theNormalizedDeltaRadii;
		double yRadius = cfR[i];
		for (int j = 0; j < theMyCfCircles[i]->GetN(); ++j) {
			double myX, myY;
			theMyCfCircles[i]->GetPoint(j, myX, myY);
			double myRadius = sqrt(pow(myX - cfX[i], 2) + pow(myY - cfY[i], 2));
			if (verbose > 0 && printPointNumber == i) std::cout << "("<< myRadius << " - " << yRadius << ")/" << yRadius << " = " << fabs(myRadius - yRadius)/yRadius << std::endl;
			theDeltaRadii.push_back(myRadius - yRadius);
			theNormalizedDeltaRadii.push_back((myRadius - yRadius)/yRadius); // the deviation of the simplified value to the calculated value wrt the calculated value
		}
		theTheDeltaRadii.push_back(theDeltaRadii);
		theTheNormalizedDeltaRadii.push_back(theNormalizedDeltaRadii);
	}
	// Calculate Means
	std::vector<double> theMeanDeltaRadii;
	std::vector<double> theMeanNormalizedDeltaRadii;
	std::vector<double> theExtremeNormalizedDeltaRadii;
	for (int i = 0; i < theTheDeltaRadii.size(); ++i) {
		double tempRadii[100];
		double tempNormalizedRadii[100];
		for (int j = 0; j < theTheDeltaRadii[i].size(); ++j) {
			tempRadii[j] = theTheDeltaRadii[i][j];
			tempNormalizedRadii[j] = theTheNormalizedDeltaRadii[i][j];
		}
		theMeanDeltaRadii.push_back(TMath::Mean(theTheDeltaRadii[i].size(), tempRadii));
		theMeanNormalizedDeltaRadii.push_back(TMath::Mean(theTheDeltaRadii[i].size(), tempNormalizedRadii));
		theExtremeNormalizedDeltaRadii.push_back(*(std::max_element(theTheNormalizedDeltaRadii[i].begin(), theTheNormalizedDeltaRadii[i].end())));
		if (verbose > 0 && printPointNumber == i) std::cout << theMeanNormalizedDeltaRadii[i] << std::endl;
	}

	// Display stuff vs. R = sqrt(x^2+y^2) of central (original) point
	TGraph * gMeansVsR = new TGraph();
	TGraph * gNormMeansVsR = new TGraph();
	TGraph * gNormExtremeVsR = new TGraph();
	for (int i = 0; i < theMeanDeltaRadii.size(); ++i) {
		double R = sqrt(pow(x[i], 2) + pow(y[i], 2));
		gMeansVsR->SetPoint(i, R, theMeanDeltaRadii[i]);
		gNormMeansVsR->SetPoint(i, R, theMeanNormalizedDeltaRadii[i]*100);
		gNormExtremeVsR->SetPoint(i, R, theExtremeNormalizedDeltaRadii[i]*100);
	}
	// gMeansVsR->SetLineColor(kCyan+2);
	gNormMeansVsR->SetLineColor(kMagenta+2);
	gNormExtremeVsR->SetLineColor(kCyan+2);
	gNormMeansVsR->SetMarkerStyle(kFullDotMedium);
	gNormExtremeVsR->SetMarkerStyle(kFullDotMedium);
	gNormMeansVsR->GetXaxis()->SetTitle("Distance to center");
	gNormMeansVsR->GetYaxis()->SetTitle("Mean, norm, #Delta r_{i} / %");
	gNormMeansVsR->GetYaxis()->SetTitleOffset(1.45);
	gNormExtremeVsR->GetXaxis()->SetTitle("Distance to center");
	gNormExtremeVsR->GetYaxis()->SetTitle("Max norm #Delta r_{i} / %");

	TCanvas * c1 = new TCanvas("c1","c1",0,0,500,500);
	mgOriginalSpace->Draw("ALP");

	TCanvas * c2 = new TCanvas("c2","c2",505,0,500,500);
	mgCfSpace->Draw("ALP");

	TCanvas * c3 = new TCanvas("c3","c3",0,200,500,500);
	gNormMeansVsR->Draw("ALP");

	TCanvas * c4 = new TCanvas("c4","c4",505,200,500,500);
	gNormExtremeVsR->Draw("ALP");

	return 1;
}