pull out the functions to ad++.h, so that it can be loaded as root marco
This commit is contained in:
Ryan@iMac
parent
3b64ba8edd
commit
78576a70d2
280
ad++.cpp
280
ad++.cpp
|
|
@ -1,34 +1,4 @@
|
|||
#include "jSymbol.h"
|
||||
#include "Qk.h"
|
||||
|
||||
#include <TROOT.h>
|
||||
#include <TSystem.h>
|
||||
#include <TAxis.h>
|
||||
#include <TF1.h>
|
||||
#include <TLatex.h>
|
||||
#include <TGraphErrors.h>
|
||||
#include <TApplication.h>
|
||||
#include <TCanvas.h>
|
||||
|
||||
|
||||
#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<std::vector<double>> 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<double> deltaDeg;
|
||||
std::vector<double> 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;
|
||||
|
||||
}
|
||||
|
|
|
|||
319
ad++.h
Normal file
319
ad++.h
Normal file
|
|
@ -0,0 +1,319 @@
|
|||
#ifndef AD_LIBARAY_H
|
||||
#define AD_LIBARAY_H
|
||||
|
||||
#include "jSymbol.h"
|
||||
#include "Qk.h"
|
||||
|
||||
#include <TROOT.h>
|
||||
#include <TSystem.h>
|
||||
#include <TAxis.h>
|
||||
#include <TF1.h>
|
||||
#include <TLatex.h>
|
||||
#include <TGraphErrors.h>
|
||||
#include <TApplication.h>
|
||||
#include <TCanvas.h>
|
||||
|
||||
|
||||
#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<std::vector<double>> 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<std::vector<double>> 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<double> deltaDeg;
|
||||
std::vector<double> 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
|
||||
Loading…
Reference in New Issue
Block a user