diff --git a/.gitignore b/.gitignore index 7ecc949..9a108d9 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ ad++ testJSymbol testJSymbol.cpp +test.C diff --git a/ad++.cpp b/ad++.cpp index d3101e7..1933ad8 100644 --- a/ad++.cpp +++ b/ad++.cpp @@ -18,6 +18,8 @@ int main(int argc, char **argv){ const int Ji = 5; const int Jf = 4; + + double sigma = 1; //====================================== @@ -26,7 +28,7 @@ int main(int argc, char **argv){ //=============================== 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"); diff --git a/ad++.h b/ad++.h index 3b99ba7..bff3944 100644 --- a/ad++.h +++ b/ad++.h @@ -57,40 +57,48 @@ void PrintRk(int k){ } } -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; - } - +double w(int M, double 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 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, double J, double sigma){ -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; 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; @@ -100,9 +108,11 @@ double Bk(int k, int J){ double TheoryAD(double * x, double *par){ /// par[0] = a; /// par[1] = delta; + /// par[2] = sigma; + /// par[3] = J; 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] ); + 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]; @@ -113,7 +123,8 @@ TGraphErrors * fitAD( std::vector> data_deg_count_error, double detRadius_cm, double targetDistance_cm , double detThickness_cm , - int Ji, int Jf){ + int Ji, int Jf, + double sigma = 0){ /// Calculate coefficient double * Qk = QK(energy_keV, detRadius_cm, targetDistance_cm, detThickness_cm); @@ -127,15 +138,14 @@ TGraphErrors * fitAD( std::vector> data_deg_count_error, 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]); - + for( int k = 0; k <=4; k += 2){ + 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 const int dataSize = (int) data_deg_count_error.size(); @@ -158,16 +168,26 @@ TGraphErrors * fitAD( std::vector> data_deg_count_error, printf("======================\n"); 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->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->SetLineWidth(2); fit->SetNpx(1000); fit->SetParameter(0, 3000); 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->GetYaxis()->SetTitle("Data"); @@ -179,6 +199,7 @@ TGraphErrors * fitAD( std::vector> data_deg_count_error, 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); + printf("Best fit sigma = %f(%f) \n", paraA2[3], paraE2[3]); TCanvas * canvas = new TCanvas("canvas", "Fitting Angular Distribution", 600, 400); canvas->cd(1); @@ -190,13 +211,16 @@ TGraphErrors * fitAD( std::vector> data_deg_count_error, 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.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.75, Form(" #sigma: %5.1f(%5.1f)", paraA2[2], paraE2[2])); text.SetTextColor(2); - 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; } @@ -207,7 +231,8 @@ void fitOldAD(std::vector> data_deg_count_error, double detRadius_cm, double targetDistance_cm , double detThickness_cm , - int Ji, int Jf){ + int Ji, int Jf, + double sigma){ /// Calculate coefficient double * Qk = QK(energy_keV, detRadius_cm, targetDistance_cm, detThickness_cm); @@ -221,14 +246,15 @@ void fitOldAD(std::vector> data_deg_count_error, if( L == 0 ) L = 1; 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); 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]); + + for( int k = 0; k <=4; k += 2){ + 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 const int dataSize = (int) data_deg_count_error.size();