Skip to content

Instantly share code, notes, and snippets.

@AndiH
Created November 27, 2012 16:15
Show Gist options
  • Save AndiH/4155157 to your computer and use it in GitHub Desktop.
Save AndiH/4155157 to your computer and use it in GitHub Desktop.
ROOT: Visualization of conformal mapping of points with isochrone radii - including a comparison between a simplified approach and a correct one
#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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment