From 78576a70d2f84b58c6d76100a49a17bdc8be2008 Mon Sep 17 00:00:00 2001 From: "Ryan@iMac" Date: Wed, 12 Oct 2022 13:48:17 -0400 Subject: [PATCH] pull out the functions to ad++.h, so that it can be loaded as root marco --- ad++.cpp | 280 ++---------------------------------------------- ad++.h | 319 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ makefile | 2 +- 3 files changed, 329 insertions(+), 272 deletions(-) create mode 100644 ad++.h diff --git a/ad++.cpp b/ad++.cpp index db41150..d3101e7 100644 --- a/ad++.cpp +++ b/ad++.cpp @@ -1,34 +1,4 @@ -#include "jSymbol.h" -#include "Qk.h" - -#include -#include -#include -#include -#include -#include -#include -#include - - -#define PI 3.14159265358979323846 - -double Racah(int j1, int j2, int J, int j3, int j12, int j23); -double Rk(int k, int L1, int L2, int J1, int J2); -void PrintRk(int k); -double w(int M, int J); -void Print_w_sum(); -double Bk(int k, int J); -double YE(double * x , double *par);/// for root Fit - - -/// This is for general fit a, delta -double Q[5] = {0}; -double B[5] = {0}; -double R1[5] = {0}; -double R2[5] = {0}; -double R3[5] = {0}; -double Fit(double * x, double *par); +#include "ad++.h" int main(int argc, char **argv){ @@ -40,162 +10,24 @@ int main(int argc, char **argv){ double targetDistance_cm = 20; double detThickness_cm = 10; - double data[][3] = { { 150.00, 2441.44, 122}, - { 131.75, 2580.11, 129}, - { 90.00, 4652.08, 232}, - { 48.75, 3023.17, 151} }; + std::vector> data = { { 150.00, 2441.44, 122}, + { 131.75, 2580.11, 129}, + { 90.00, 4652.08, 232}, + { 48.75, 3023.17, 151} }; const int Ji = 5; const int Jf = 4; //====================================== - /// detector acceptance - double * Qk = QK(energy_keV, detRadius_cm, targetDistance_cm, detThickness_cm); - Q[0] = 1; - Q[2] = Qk[0]; - Q[4] = Qk[1]; - - printf("Qk2 : %f \n", Q[2]); - printf("Qk4 : %f \n", Q[4]); - const int dataSize = sizeof(data)/sizeof(double)/3; - - /// for TGraphErrors - double x[dataSize]; - double y[dataSize]; - double ex[dataSize]; - double ey[dataSize]; - - printf("============= Data :\n"); - for( int i = 0; i < dataSize; i++){ - printf("%2d | %8.2f, %8.2f(%4.0f) \n", i, data[i][0], data[i][1], data[i][2]); - - x[i] = data[i][0] * PI/180; - y[i] = data[i][1]; - ey[i] = data[i][2]; - ex[i] = 0.; - } - printf("======================\n"); - - TGraphErrors * gExp = new TGraphErrors( dataSize, x, y, ex, ey); - gExp->SetTitle(""); - gExp->GetXaxis()->SetTitle("Angle [rad]"); - gExp->GetYaxis()->SetTitle("Data"); - - TCanvas * c1 = new TCanvas("c1", "c1", 1000, 500); - c1->Divide(2, 1); - c1->cd(1); - gExp->Draw("AP*"); - - c1->Modified(); - gSystem->ProcessEvents(); + //=============================== Fit a, delta at once + TGraphErrors * gExp = fitAD(data, energy_keV, detRadius_cm, targetDistance_cm, detThickness_cm, Ji, Jf); - ///======== Fitting the experimental distribution with a0+ a2*P(2,cos(theta)) + a4 * P(4, cos(theta)) - TF1 * f1 = new TF1("f1", YE, 0, PI, 3); + //=============================== Old AD code way + fitOldAD(data, energy_keV, detRadius_cm, targetDistance_cm, detThickness_cm, Ji, Jf); - f1->SetLineColor(4); - f1->SetLineWidth(2); - f1->SetNpx(1000); - f1->SetParameter(0, 2000); - f1->SetParameter(1, -2000); - f1->SetParameter(2, 2000); - - gExp->Fit("f1", ""); - - const Double_t* paraE = f1->GetParErrors(); - const Double_t* paraA = f1->GetParameters(); - - double A0 = paraA[0]; - double A2 = paraA[1]; - double A4 = paraA[2]; - - ///=================================== Fit with Theritical - - - int L = abs(Ji - Jf); - if( L == 0 ) L = 1; - - - for( int k = 0; k <= 4; k += 2){ - B[k] = Bk(k, Ji); - R1[k] = Rk(k, L , L , Ji, Jf); - R2[k] = Rk(k, L , L+1, Ji, Jf); - R3[k] = Rk(k, L+1, L+1, Ji, Jf); - } - - std::vector deltaDeg; - std::vector LogChiSq; - - for( float deltaAngle = -90; deltaAngle <= 90 ; deltaAngle += 2. ){ - - double delta = tan(deltaAngle * PI/ 180.); - - double chiSq = 0; - for( int i = 0; i < dataSize; i++){ - - double YT = 0; - - for( int k = 0; k <= 4; k += 2){ - YT += Q[k] * B[k] * LegendreP(k, x[i]) * ( R1[k] + 2 * delta * R2[k] + delta*delta* R3[k] ) / (1 + delta * delta) ; - } - - //printf(" YT : %f , YE : %f \n", A0 * YT, y[i]); - - chiSq += pow( A0 * YT - y[i], 2)/ dataSize / ey[i] / ey[i]; - - } - - deltaDeg.push_back(deltaAngle); - LogChiSq.push_back(log(chiSq)); - //printf(" %6.2f deg, %8.3f \n", deltaAngle, log(chiSq)); - } - - c1->cd(2); - TGraph * gDelta = new TGraph((int) deltaDeg.size(), &deltaDeg[0], &LogChiSq[0]); - gDelta->SetTitle(""); - gDelta->GetXaxis()->SetTitle("aTan(Mixing Ratio) [Deg]"); - gDelta->GetYaxis()->SetTitle("Log(chi-sq)"); - - gDelta->Draw("APL*"); - - c1->Modified(); - gSystem->ProcessEvents(); - - ///=============================== Fit a, delta at once - - TF1 * fit = new TF1("fit", Fit, 0, PI, 2); - fit->SetLineColor(2); - fit->SetLineWidth(2); - fit->SetNpx(1000); - fit->SetParameter(0, 3000); - fit->SetParameter(1, 1); - - gExp->Fit("fit", "q"); - - const Double_t * paraE2 = fit->GetParErrors(); - const Double_t * paraA2 = fit->GetParameters(); - - printf("===================== \n"); - printf("Best fit Amp = %f(%f)\n", paraA2[0], paraE2[0]); - printf("Best fit delta = %f(%f) = %f(%f) deg\n", paraA2[1], paraE2[1], atan(paraA2[1]) * 180/PI, atan(paraE2[1])*180/PI); - - c1->cd(1); - fit->Draw("same"); - f1->Draw("same"); - - TLatex text; - text.SetNDC(); - text.SetTextFont(82); - text.SetTextSize(0.04); - - text.DrawLatex(0.12, 0.85, Form("#delta : %5.1f(%5.1f) deg", atan(paraA2[1]) * 180/PI, atan(paraE2[1])*180/PI)); - text.DrawLatex(0.12, 0.80, Form("Amp: %5.1f(%5.1f)", paraA2[0], paraE2[0])); - - text.SetTextColor(2); - - text.DrawLatex(0.8, 0.8, Form("%d->%d", Ji, Jf)); printf("===================== Crtl+C to end.\n"); @@ -204,97 +36,3 @@ int main(int argc, char **argv){ return 0; } - -//######################### -///use for root fit -double YE(double * x , double *par){ - /// x[0] = angle in radian; - /// par[0] = a0; - /// par[1] = a2; - /// par[2] = a4; - - return par[0] + par[1] * LegendreP(2, x[0]) + par[2] * LegendreP(4, x[0]); -} - -///use for fit a, delta -double Fit(double * x, double *par){ - /// par[0] = a; - /// par[1] = delta; - - double result = 0; - for( int k = 0; k <= 4; k += 2){ - - result += Q[k]*B[k] * LegendreP(k, x[0]) * ( R1[k] + 2*par[1]*R2[k] + par[1]*par[1]*R3[k] ) / (1 + par[1]*par[1] ); - } - - return result * par[0]; -} - - -double Racah(int j1, int j2, int J, int j3, int j12, int j23){ - - return pow(-1, j1+j2+j3+J) * SixJSymbol(j1, j2, j12, j3, J, j23); - -} - -double Rk(int k, int L1, int L2, int J1, int J2){ - - return pow(-1, 1+J1-J2+L2-L1-k) * pow((2*J1+1)*(2*L1+1)*(2*L2+1), 0.5) * CGcoeff(k, 0, L1, 1, L2, -1) * Racah(J1, J1, L1, L2, k, J2); - -} - -void PrintRk(int k){ - - for( int J1 = 1; J1 < 8; J1 ++){ - printf("==============================\n"); - for( int J2 = 0; J2 < 8; J2++){ - - int L = abs(J1 - J2); - if( L == 0 ) L = 1; - printf("%d %d | %10.6f, %10.6f, %10.6f \n", J1, J2, Rk(k, L, L, J1, J2), Rk(k, L, L+1, J1, J2), Rk(k, L+1, L+1, J1, J2)); - - } - } -} - -double w(int M, int J){ - double sigma = 1; - switch (J) { - case 0 : sigma = 0.3989422804014327 ; break; - case 1 : sigma = 0.5723377817486753 ; break; - case 2 : sigma = 0.7013915463848625 ; break; - case 3 : sigma = 0.8091713162791643 ; break; - case 4 : sigma = 0.9037290722944527 ; break; - case 5 : sigma = 0.9890249035482789 ; break; - case 6 : sigma = 1.0673592868302038 ; break; - case 7 : sigma = 1.140210831444403 ; break; - case 8 : sigma = 1.208599155456379 ; break; - case 9 : sigma = 1.273313297516925 ; break; - case 10 : sigma = 1.335665551821612 ; break; - } - - return 1./sqrt(2 * PI) / sigma * exp( - M*M / 2. / sigma/sigma); - -} - -void Print_w_sum(){ /// Check the normalization of w(m), sum(w(m)) = 1 for all J. - - for( int J = 0; J < 11; J++){ - double w_sum = 0; - for( int m = -J ; m < J+1; m++){ - w_sum += w(m, J); - } - printf("%2d %.5f\n", J, w_sum); - } -} - -double Bk(int k, int J){ - - double sum = 0; - for(int m = -J; m < J+1 ; m++){ - sum += w(m, J) * pow(-1, J-m) * pow(2*J+1,0.5) * CGcoeff(k, 0, J, m, J, -m); - } - - return sum; - -} diff --git a/ad++.h b/ad++.h new file mode 100644 index 0000000..3b99ba7 --- /dev/null +++ b/ad++.h @@ -0,0 +1,319 @@ +#ifndef AD_LIBARAY_H +#define AD_LIBARAY_H + +#include "jSymbol.h" +#include "Qk.h" + +#include +#include +#include +#include +#include +#include +#include +#include + + +#define PI 3.14159265358979323846 + + +/// These are the coefficients for theritical distribution, index 1, 3 are not used. +double Q[5] = {0}; +double B[5] = {0}; +double R1[5] = {0}; +double R2[5] = {0}; +double R3[5] = {0}; + +///use for root fit +double YE(double * x , double *par){ + /// x[0] = angle in radian; + /// par[0] = a0; + /// par[1] = a2; + /// par[2] = a4; + + return par[0] + par[1] * LegendreP(2, x[0]) + par[2] * LegendreP(4, x[0]); +} + + + +double Racah(int j1, int j2, int J, int j3, int j12, int j23){ + return pow(-1, j1+j2+j3+J) * SixJSymbol(j1, j2, j12, j3, J, j23); +} + +double Rk(int k, int L1, int L2, int J1, int J2){ + return pow(-1, 1+J1-J2+L2-L1-k) * pow((2*J1+1)*(2*L1+1)*(2*L2+1), 0.5) * CGcoeff(k, 0, L1, 1, L2, -1) * Racah(J1, J1, L1, L2, k, J2); +} + +void PrintRk(int k){ + for( int J1 = 1; J1 < 8; J1 ++){ + printf("==============================\n"); + for( int J2 = 0; J2 < 8; J2++){ + + int L = abs(J1 - J2); + if( L == 0 ) L = 1; + printf("%d %d | %10.6f, %10.6f, %10.6f \n", J1, J2, Rk(k, L, L, J1, J2), Rk(k, L, L+1, J1, J2), Rk(k, L+1, L+1, J1, J2)); + + } + } +} + +double w(int M, int J){ + double sigma = 1; + switch (J) { + case 0 : sigma = 0.3989422804014327 ; break; + case 1 : sigma = 0.5723377817486753 ; break; + case 2 : sigma = 0.7013915463848625 ; break; + case 3 : sigma = 0.8091713162791643 ; break; + case 4 : sigma = 0.9037290722944527 ; break; + case 5 : sigma = 0.9890249035482789 ; break; + case 6 : sigma = 1.0673592868302038 ; break; + case 7 : sigma = 1.140210831444403 ; break; + case 8 : sigma = 1.208599155456379 ; break; + case 9 : sigma = 1.273313297516925 ; break; + case 10 : sigma = 1.335665551821612 ; break; + } + + return 1./sqrt(2 * PI) / sigma * exp( - M*M / 2. / sigma/sigma); +} + +void Print_w_sum(){ /// Check the normalization of w(m), sum(w(m)) = 1 for all J. + + for( int J = 0; J < 11; J++){ + double w_sum = 0; + for( int m = -J ; m < J+1; m++){ + w_sum += w(m, J); + } + printf("%2d %.5f\n", J, w_sum); + } +} + +double Bk(int k, int J){ + double sum = 0; + for(int m = -J; m < J+1 ; m++){ + sum += w(m, J) * pow(-1, J-m) * pow(2*J+1,0.5) * CGcoeff(k, 0, J, m, J, -m); + } + + return sum; +} + +///use for fit a, delta +double TheoryAD(double * x, double *par){ + /// par[0] = a; + /// par[1] = delta; + double result = 0; + for( int k = 0; k <= 4; k += 2){ + result += Q[k]*B[k] * LegendreP(k, x[0]) * ( R1[k] + 2*par[1]*R2[k] + par[1]*par[1]*R3[k] ) / (1 + par[1]*par[1] ); + } + + return result * par[0]; +} + +TGraphErrors * fitAD( std::vector> data_deg_count_error, + double energy_keV, + double detRadius_cm, + double targetDistance_cm , + double detThickness_cm , + int Ji, int Jf){ + + /// Calculate coefficient + double * Qk = QK(energy_keV, detRadius_cm, targetDistance_cm, detThickness_cm); + Q[0] = 1; + Q[2] = Qk[0]; + Q[4] = Qk[1]; + + delete Qk; + + int L = abs(Ji - Jf); + if( L == 0 ) L = 1; + + for( int k = 0; k <= 4; k += 2){ + B[k] = Bk(k, Ji); + R1[k] = Rk(k, L , L , Ji, Jf); + R2[k] = Rk(k, L , L+1, Ji, Jf); + R3[k] = Rk(k, L+1, L+1, Ji, Jf); + } + + printf("Qk2 : %f \n", Q[2]); + printf("Qk4 : %f \n", Q[4]); + + + /// Create TGraphEorrs + const int dataSize = (int) data_deg_count_error.size(); + + /// for TGraphErrors + double x[dataSize]; + double y[dataSize]; + double ex[dataSize]; + double ey[dataSize]; + + printf("============= Data :\n"); + for( int i = 0; i < dataSize; i++){ + printf("%2d | %8.2f, %8.2f(%4.0f) \n", i, data_deg_count_error[i][0], data_deg_count_error[i][1], data_deg_count_error[i][2]); + + x[i] = data_deg_count_error[i][0] * PI/180; + y[i] = data_deg_count_error[i][1]; + ey[i] = data_deg_count_error[i][2]; + ex[i] = 0.; + } + printf("======================\n"); + + TGraphErrors * gExp = new TGraphErrors( dataSize, x, y, ex, ey); + gExp->SetTitle(""); + gExp->GetXaxis()->SetTitle("Angle [rad]"); + gExp->GetYaxis()->SetTitle("Data"); + + TF1 * fit = new TF1("fit", TheoryAD, 0, PI, 2); + fit->SetLineColor(2); + fit->SetLineWidth(2); + fit->SetNpx(1000); + fit->SetParameter(0, 3000); + fit->SetParameter(1, 1); + fit->GetXaxis()->SetTitle("Angle [rad]"); + fit->GetYaxis()->SetTitle("Data"); + + gExp->Fit("fit"); + + const Double_t * paraE2 = fit->GetParErrors(); + const Double_t * paraA2 = fit->GetParameters(); + + printf("===================== \n"); + printf("Best fit Amp = %f(%f)\n", paraA2[0], paraE2[0]); + printf("Best fit delta = %f(%f) = %f(%f) deg\n", paraA2[1], paraE2[1], atan(paraA2[1]) * 180/PI, atan(paraE2[1])*180/PI); + + TCanvas * canvas = new TCanvas("canvas", "Fitting Angular Distribution", 600, 400); + canvas->cd(1); + gExp->Draw("AP*"); + fit->Draw("same"); + + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.04); + + text.DrawLatex(0.12, 0.85, Form("#delta : %5.1f(%5.1f) deg", atan(paraA2[1]) * 180/PI, atan(paraE2[1])*180/PI)); + text.DrawLatex(0.12, 0.80, Form("Amp: %5.1f(%5.1f)", paraA2[0], paraE2[0])); + + text.SetTextColor(2); + + text.DrawLatex(0.8, 0.8, Form("%d->%d", Ji, Jf)); + + return gExp; + +} + + +void fitOldAD(std::vector> data_deg_count_error, + double energy_keV, + double detRadius_cm, + double targetDistance_cm , + double detThickness_cm , + int Ji, int Jf){ + + /// Calculate coefficient + double * Qk = QK(energy_keV, detRadius_cm, targetDistance_cm, detThickness_cm); + Q[0] = 1; + Q[2] = Qk[0]; + Q[4] = Qk[1]; + + delete Qk; + + int L = abs(Ji - Jf); + if( L == 0 ) L = 1; + + for( int k = 0; k <= 4; k += 2){ + B[k] = Bk(k, Ji); + R1[k] = Rk(k, L , L , Ji, Jf); + R2[k] = Rk(k, L , L+1, Ji, Jf); + R3[k] = Rk(k, L+1, L+1, Ji, Jf); + } + + printf("Qk2 : %f \n", Q[2]); + printf("Qk4 : %f \n", Q[4]); + + /// Create TGraphEorrs + const int dataSize = (int) data_deg_count_error.size(); + + /// for TGraphErrors + double x[dataSize]; + double y[dataSize]; + double ex[dataSize]; + double ey[dataSize]; + + printf("============= Data :\n"); + for( int i = 0; i < dataSize; i++){ + printf("%2d | %8.2f, %8.2f(%4.0f) \n", i, data_deg_count_error[i][0], data_deg_count_error[i][1], data_deg_count_error[i][2]); + + x[i] = data_deg_count_error[i][0] * PI/180; + y[i] = data_deg_count_error[i][1]; + ey[i] = data_deg_count_error[i][2]; + ex[i] = 0.; + } + printf("======================\n"); + + TGraphErrors * gExp = new TGraphErrors( dataSize, x, y, ex, ey); + gExp->SetTitle(""); + gExp->GetXaxis()->SetTitle("Angle [rad]"); + gExp->GetYaxis()->SetTitle("Data"); + + + ///======== Fitting the experimental distribution with a0+ a2*P(2,cos(theta)) + a4 * P(4, cos(theta)) + TF1 * f1 = new TF1("f1", YE, 0, PI, 3); + + f1->SetLineColor(4); + f1->SetLineWidth(2); + f1->SetNpx(1000); + f1->SetParameter(0, 2000); + f1->SetParameter(1, -2000); + f1->SetParameter(2, 2000); + + gExp->Fit("f1", ""); + + const Double_t* paraE = f1->GetParErrors(); + const Double_t* paraA = f1->GetParameters(); + + double A0 = paraA[0]; + + ///=================================== Find Chi-sq + std::vector deltaDeg; + std::vector LogChiSq; + + for( float deltaAngle = -90; deltaAngle <= 90 ; deltaAngle += 2. ){ + + double delta = tan(deltaAngle * PI/ 180.); + + double chiSq = 0; + for( int i = 0; i < dataSize; i++){ + + double YT = 0; + + for( int k = 0; k <= 4; k += 2){ + YT += Q[k] * B[k] * LegendreP(k, x[i] * PI/180.) * ( R1[k] + 2 * delta * R2[k] + delta*delta* R3[k] ) / (1 + delta * delta) ; + } + + chiSq += pow( A0 * YT - y[i], 2)/ dataSize / ey[i] / ey[i]; + + } + + deltaDeg.push_back(deltaAngle); + LogChiSq.push_back(log(chiSq)); + //printf(" %6.2f deg, %8.3f \n", deltaAngle, log(chiSq)); + } + + TGraph * gDelta = new TGraph((int) deltaDeg.size(), &deltaDeg[0], &LogChiSq[0]); + gDelta->SetTitle(""); + gDelta->GetXaxis()->SetTitle("aTan(Mixing Ratio) [Deg]"); + gDelta->GetYaxis()->SetTitle("Log(chi-sq)"); + + ///======== Plot + + TCanvas * c1 = new TCanvas("c1", "Chi-sq of mixing ratio", 1000, 500); + c1->Divide(2, 1); + c1->cd(1); + gExp->Draw("AP*"); + + c1->cd(2); + gDelta->Draw("APL*"); + +} + +#endif diff --git a/makefile b/makefile index 7fccf09..915f837 100644 --- a/makefile +++ b/makefile @@ -1,2 +1,2 @@ -ad++ : ad++.cpp Qk.h jSymbol.h +ad++ : ad++.cpp Qk.h jSymbol.h ad++.h g++ ad++.cpp -o ad++ `root-config --cflags --glibs`