rtang/ADPlusPlus
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

change w(m), fitAD can fit sigma

Browse Source
This commit is contained in:
Ryan@iMac 2022-10-12 18:24:52 -04:00
parent ddfbf67404
commit 5ae128573d
3 changed files with 73 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored
View File

@ -1,3 +1,4 @@
ad++ ad++
testJSymbol testJSymbol
testJSymbol.cpp testJSymbol.cpp
test.C

4
ad++.cpp
View File

@ -18,6 +18,8 @@ int main(int argc, char **argv){
const int Ji = 5; const int Ji = 5;
const int Jf = 4; const int Jf = 4;
double sigma = 1;
//====================================== //======================================
@ -26,7 +28,7 @@ int main(int argc, char **argv){
//=============================== Old AD code way //=============================== Old AD code way
fitOldAD(data, energy_keV, detRadius_cm, targetDistance_cm, detThickness_cm, Ji, Jf); fitOldAD(data, energy_keV, detRadius_cm, targetDistance_cm, detThickness_cm, Ji, Jf, sigma);
printf("===================== Crtl+C to end.\n"); printf("===================== Crtl+C to end.\n");

112
ad++.h
View File

@ -57,40 +57,48 @@ void PrintRk(int k){
} }
} }
double w(int M, int J){ double w(int M, double sigma){
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); 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. //double wFix(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 w(M, 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);
// }
//}
for( int J = 0; J < 11; J++){ double Bk(int k, double J, double sigma){
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 norm_w = 0;
for( int m = -J; m < J+1; m++) norm_w += w(m, sigma);
double sum = 0; double sum = 0;
for(int m = -J; m < J+1 ; m++){ 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); sum += ( w(m, sigma) / norm_w ) * pow(-1, J-m) * pow(2*J+1,0.5) * CGcoeff(k, 0, J, m, J, -m);
} }
return sum; return sum;
@ -100,9 +108,11 @@ double Bk(int k, int J){
double TheoryAD(double * x, double *par){ double TheoryAD(double * x, double *par){
/// par[0] = a; /// par[0] = a;
/// par[1] = delta; /// par[1] = delta;
/// par[2] = sigma;
/// par[3] = J;
double result = 0; double result = 0;
for( int k = 0; k <= 4; k += 2){ 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] ); result += Q[k]*Bk(k, par[3], par[2]) * 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]; return result * par[0];
@ -113,7 +123,8 @@ TGraphErrors * fitAD( std::vector<std::vector<double>> data_deg_count_error,
double detRadius_cm, double detRadius_cm,
double targetDistance_cm , double targetDistance_cm ,
double detThickness_cm , double detThickness_cm ,
int Ji, int Jf){ int Ji, int Jf,
double sigma = 0){
/// Calculate coefficient /// Calculate coefficient
double * Qk = QK(energy_keV, detRadius_cm, targetDistance_cm, detThickness_cm); double * Qk = QK(energy_keV, detRadius_cm, targetDistance_cm, detThickness_cm);
@ -127,15 +138,14 @@ TGraphErrors * fitAD( std::vector<std::vector<double>> data_deg_count_error,
if( L == 0 ) L = 1; if( L == 0 ) L = 1;
for( int k = 0; k <= 4; k += 2){ for( int k = 0; k <= 4; k += 2){
B[k] = Bk(k, Ji);
R1[k] = Rk(k, L , L , Ji, Jf); R1[k] = Rk(k, L , L , Ji, Jf);
R2[k] = Rk(k, L , L+1, Ji, Jf); R2[k] = Rk(k, L , L+1, Ji, Jf);
R3[k] = Rk(k, L+1, L+1, Ji, Jf); R3[k] = Rk(k, L+1, L+1, Ji, Jf);
} }
printf("Qk2 : %f \n", Q[2]); for( int k = 0; k <=4; k += 2){
printf("Qk4 : %f \n", Q[4]); printf("Qk(%d) : %10.8f, R1(%d) : %10.8f, R2(%d) : %10.8f, R3(%d) : %10.8f\n", k, Q[k], k, R1[k], k, R2[k], k, R3[k]);
}
/// Create TGraphEorrs /// Create TGraphEorrs
const int dataSize = (int) data_deg_count_error.size(); const int dataSize = (int) data_deg_count_error.size();
@ -158,16 +168,26 @@ TGraphErrors * fitAD( std::vector<std::vector<double>> data_deg_count_error,
printf("======================\n"); printf("======================\n");
TGraphErrors * gExp = new TGraphErrors( dataSize, x, y, ex, ey); TGraphErrors * gExp = new TGraphErrors( dataSize, x, y, ex, ey);
gExp->SetTitle(""); gExp->SetTitle(Form("%.0f keV, detRadius = %.1f cm, detThickness = %.1f cm, distance = %.1f cm", energy_keV, detRadius_cm, detThickness_cm, targetDistance_cm));
gExp->GetXaxis()->SetTitle("Angle [rad]"); gExp->GetXaxis()->SetTitle("Angle [rad]");
gExp->GetYaxis()->SetTitle("Data"); gExp->GetYaxis()->SetTitle("Data");
TF1 * fit = new TF1("fit", TheoryAD, 0, PI, 2); TF1 * fit = new TF1("fit", TheoryAD, 0, PI, 4);
fit->SetLineColor(2); fit->SetLineColor(2);
fit->SetLineWidth(2); fit->SetLineWidth(2);
fit->SetNpx(1000); fit->SetNpx(1000);
fit->SetParameter(0, 3000); fit->SetParameter(0, 3000);
fit->SetParameter(1, 1); fit->SetParameter(1, 1);
fit->SetParameter(2, 1);
fit->FixParameter(3, Ji);
fit->SetParLimits(1, -30, 30);
fit->SetParLimits(2, 0.2, Ji);
if( sigma > 0 ) fit->FixParameter(2, sigma);
fit->GetXaxis()->SetTitle("Angle [rad]"); fit->GetXaxis()->SetTitle("Angle [rad]");
fit->GetYaxis()->SetTitle("Data"); fit->GetYaxis()->SetTitle("Data");
@ -179,6 +199,7 @@ TGraphErrors * fitAD( std::vector<std::vector<double>> data_deg_count_error,
printf("===================== \n"); printf("===================== \n");
printf("Best fit Amp = %f(%f)\n", paraA2[0], paraE2[0]); 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); printf("Best fit delta = %f(%f) = %f(%f) deg\n", paraA2[1], paraE2[1], atan(paraA2[1]) * 180/PI, atan(paraE2[1])*180/PI);
printf("Best fit sigma = %f(%f) \n", paraA2[3], paraE2[3]);
TCanvas * canvas = new TCanvas("canvas", "Fitting Angular Distribution", 600, 400); TCanvas * canvas = new TCanvas("canvas", "Fitting Angular Distribution", 600, 400);
canvas->cd(1); canvas->cd(1);
@ -190,13 +211,16 @@ TGraphErrors * fitAD( std::vector<std::vector<double>> data_deg_count_error,
text.SetTextFont(82); text.SetTextFont(82);
text.SetTextSize(0.04); 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.85, Form(" #delta: %5.2f(%5.2f) = %5.1f(%5.1f) deg", paraA2[1], paraE2[1], 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.DrawLatex(0.12, 0.80, Form(" Amp: %5.1f(%5.1f)", paraA2[0], paraE2[0]));
text.DrawLatex(0.12, 0.75, Form(" #sigma: %5.1f(%5.1f)", paraA2[2], paraE2[2]));
text.SetTextColor(2); text.SetTextColor(2);
text.DrawLatex(0.8, 0.8, Form("%d->%d", Ji, Jf)); text.DrawLatex(0.8, 0.8, Form("%d->%d", Ji, Jf));
for( int k = 0; k <= 4; k += 2) printf("Bk(%d) : %10.7f\n", k, Bk(k, Ji, paraA2[2]));
return gExp; return gExp;
} }
@ -207,7 +231,8 @@ void fitOldAD(std::vector<std::vector<double>> data_deg_count_error,
double detRadius_cm, double detRadius_cm,
double targetDistance_cm , double targetDistance_cm ,
double detThickness_cm , double detThickness_cm ,
int Ji, int Jf){ int Ji, int Jf,
double sigma){
/// Calculate coefficient /// Calculate coefficient
double * Qk = QK(energy_keV, detRadius_cm, targetDistance_cm, detThickness_cm); double * Qk = QK(energy_keV, detRadius_cm, targetDistance_cm, detThickness_cm);
@ -221,14 +246,15 @@ void fitOldAD(std::vector<std::vector<double>> data_deg_count_error,
if( L == 0 ) L = 1; if( L == 0 ) L = 1;
for( int k = 0; k <= 4; k += 2){ for( int k = 0; k <= 4; k += 2){
B[k] = Bk(k, Ji); B[k] = Bk(k, Ji, sigma);
R1[k] = Rk(k, L , L , Ji, Jf); R1[k] = Rk(k, L , L , Ji, Jf);
R2[k] = Rk(k, L , L+1, Ji, Jf); R2[k] = Rk(k, L , L+1, Ji, Jf);
R3[k] = Rk(k, L+1, L+1, Ji, Jf); R3[k] = Rk(k, L+1, L+1, Ji, Jf);
} }
printf("Qk2 : %f \n", Q[2]); for( int k = 0; k <=4; k += 2){
printf("Qk4 : %f \n", Q[4]); printf("Qk(%d) : %f, Bk(%d) : %f, R1(%d) : %f, R2(%d) : %f, R3(%d) : %f\n", k, Q[k], k, B[k], k, R1[k], k, R2[k], k, R3[k]);
}
/// Create TGraphEorrs /// Create TGraphEorrs
const int dataSize = (int) data_deg_count_error.size(); const int dataSize = (int) data_deg_count_error.size();