/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

Public Member Functions | Private Member Functions | Static Private Member Functions | Private Attributes
CoarseScaleCalib Class Reference

#include <CoarseScaleCalib.h>

List of all members.

Public Member Functions

 CoarseScaleCalib (vector< double > adc, vector< double > scale, double adc1pedestal, double adc2pedestal, double adc1pedestalWidth, double adc2pedestalWidth)
void Fit ()
int GetEmpiricalDist (int scale, int adc)
double prob (double x, int scale)
void Get_quality (double &, double &)
double Get_mu ()
double Get_sigma ()
double Get_rho2 ()
double Get_M ()
double Get_N1 ()
double Get_N2 ()
double Get_ADC2_Smearing ()

Private Member Functions

double calc_mu (double sigma, double m, double rho2)
double adc1_likelihood (double sigma)
double adc2_likelihood (double m, double rho2)
double likelihood (double sigma, double m, double rho2)

Static Private Member Functions

static double q1_unnormalized (double x, double mu, double sigma, double eta1, double rho1)
static double p_q1_unnormalized (double x, void *parameters)
static double q1_integral (double x, double mu, double sigma, double eta1, double rho1)
static double q2_unnormalized (double x, double mu, double sigma, double m, double eta1, double eta2, double rho2)
static double p_q2_unnormalized (double x, void *parameters)
static double q2_integral (double x, double mu, double sigma, double m, double eta1, double eta2, double rho2)
static double constraint (double x, void *parameter)
static double p_adc1_likelihood (double x, void *parameters)
static double p_adc2_likelihood (const gsl_vector *v, void *parameters)
static double p_likelihood (const gsl_vector *v, void *parameters)

Private Attributes

vector< double > adc1
vector< double > adc2
double N
double N1
double N2
double PED1
double PED2
double PEDWID1
double PEDWID2
double M
double MU
double SIGMA
double RHO2
STATUS CalibStatus

Detailed Description

Definition at line 11 of file CoarseScaleCalib.h.


Constructor & Destructor Documentation

CoarseScaleCalib::CoarseScaleCalib ( vector< double >  adc,
vector< double >  scale,
double  adc1pedestal,
double  adc2pedestal,
double  adc1pedestalWidth,
double  adc2pedestalWidth 
)

Definition at line 29 of file CoarseScaleCalib.cc.

                                                                                       {

        if (adc.size() != scale.size()) throw("ADC and scale vectors must have the same size.");

        // Separate the adc vector into two spectra (by scale)
        for (unsigned int i=0; i<adc.size(); i++) {
                if (scale.at(i)==1 && adc.at(i)>=RANGE_LOW) adc1.push_back(adc.at(i));
                if (scale.at(i)==2 && adc.at(i)<=RANGE_HIGH && adc.at(i)>=0) adc2.push_back(adc.at(i));
        }

        

        // Counts in each spectrum
        N1 = adc1.size();
        N2 = adc2.size();
        N = N1+N2;

        // Pedestal-related quantities
        PED1 = adc1pedestal;
        PED2 = adc2pedestal;
        PEDWID1 = adc1pedestalWidth;
        PEDWID2 = adc2pedestalWidth;

        // Fit results initialized
        MU = 0;
        SIGMA = SIGMA_GUESS;
        M = M_GUESS;
        RHO2 = RHO2_GUESS;

        // Set status
        CalibStatus = NOTYETFIT;
}

Member Function Documentation

void CoarseScaleCalib::Fit ( )

Definition at line 334 of file CoarseScaleCalib.cc.

                           {

  if (N<MINCTS) {
    CalibStatus = INSUFFICIENTCOUNTS;
    throw("Not enough counts in ADC spectrum");
  }
  
  if (N1<MINFRAC*N || N2<MINFRAC*N) {
    CalibStatus = IMBALANCED;
    throw("Scales are not balanced (one scale has many more counts than the other)");
  }
  
  // ==================================================================
  //   Minimize the ADC1 likelihood first to get an estimate for sigma
  // ==================================================================
  
  CoarseScaleCalib* params[] = { this };
  gsl_function f;
  f.function = &CoarseScaleCalib::p_adc1_likelihood;
  f.params = params;
  
  const gsl_min_fminimizer_type* min_alg = gsl_min_fminimizer_goldensection;
  gsl_min_fminimizer* myMinimizer = gsl_min_fminimizer_alloc(min_alg);
  gsl_min_fminimizer_set(myMinimizer, &f, SIGMA_GUESS, 1, 1000);
  
  int status;
  int iter = 0, max_iter = 1000;
  double thismin, x_lo, x_hi;
  do {
    iter++;
    status = gsl_min_fminimizer_iterate(myMinimizer);
    
    thismin = gsl_min_fminimizer_x_minimum(myMinimizer);
    x_lo = gsl_min_fminimizer_x_lower(myMinimizer);
    x_hi = gsl_min_fminimizer_x_upper(myMinimizer);
    status = gsl_min_test_interval(x_lo, x_hi, 0.1, 0.0);
  } while(status == GSL_CONTINUE && iter < max_iter);
  
  gsl_min_fminimizer_free(myMinimizer);
  
  SIGMA = thismin;
  
  // =============================================================================================
  //   Use this value of SIGMA and minimize the ADC2 likelihood to get an estimate for M and RHO2
  // =============================================================================================
  
  gsl_multimin_function f2;
  CoarseScaleCalib* params2[] = { this };
  f2.n = 2;
  f2.f = &CoarseScaleCalib::p_adc2_likelihood;
  f2.params = params2;
  
  gsl_vector* first_guess;
  first_guess = gsl_vector_alloc(2);
  gsl_vector_set(first_guess, 0, M_GUESS);
  gsl_vector_set(first_guess, 1, RHO2_GUESS);
  
  gsl_vector* trial_step;
  trial_step = gsl_vector_alloc(2);
  gsl_vector_set(trial_step, 0, 2);
  gsl_vector_set(trial_step, 1, 20);
  
  const gsl_multimin_fminimizer_type* multimin_alg = gsl_multimin_fminimizer_nmsimplex;
  gsl_multimin_fminimizer* myMinimizer2 = gsl_multimin_fminimizer_alloc(multimin_alg, 2);
  gsl_multimin_fminimizer_set(myMinimizer2, &f2, first_guess, trial_step);
  
  iter = 0;
  do {
    iter++;
    status = gsl_multimin_fminimizer_iterate(myMinimizer2);
    if (status) break;
    status = gsl_multimin_test_size(gsl_multimin_fminimizer_size(myMinimizer2), 1);
    
  } while (status == GSL_CONTINUE && iter < max_iter);
  
  M = gsl_vector_get(myMinimizer2->x,0);
  RHO2 = gsl_vector_get(myMinimizer2->x,1);
  
  gsl_multimin_fminimizer_free(myMinimizer2);
  gsl_vector_free(first_guess);
  gsl_vector_free(trial_step);
  
  // ============================================================================
  //   Use the results so far as a first guess and minimize the total likelihood
  // ============================================================================
  
  gsl_multimin_function f3;
  CoarseScaleCalib* params3[] = { this };
  f3.n = 3;
  f3.f = &CoarseScaleCalib::p_likelihood;
  f3.params = params3;
  
  gsl_vector* guess;
  guess = gsl_vector_alloc(3);
  gsl_vector_set(guess, 0, SIGMA);
  gsl_vector_set(guess, 1, M);
  gsl_vector_set(guess, 2, RHO2);
  
  // Uses smaller initial step size since we have a pretty good initial guess
  gsl_vector* step;
  step = gsl_vector_alloc(3);
  gsl_vector_set(step, 0, 10);
  gsl_vector_set(step, 1, 0.2);
  gsl_vector_set(step, 2, 2);
  
  const gsl_multimin_fminimizer_type* multimin_alg3 = gsl_multimin_fminimizer_nmsimplex;
  gsl_multimin_fminimizer* myMinimizer3 = gsl_multimin_fminimizer_alloc(multimin_alg3, 3);
  gsl_multimin_fminimizer_set(myMinimizer3, &f3, guess, step);
  
  iter = 0;
  do {
    iter++;
    status = gsl_multimin_fminimizer_iterate(myMinimizer3);
    if (status) break;
    status = gsl_multimin_test_size(gsl_multimin_fminimizer_size(myMinimizer3), 0.1); // Lower tolerance
    
  } while (status == GSL_CONTINUE && iter < max_iter);
  
  SIGMA = gsl_vector_get(myMinimizer3->x,0);
  M = gsl_vector_get(myMinimizer3->x,1);
  RHO2 = gsl_vector_get(myMinimizer3->x,2);
  MU = calc_mu(SIGMA,M,RHO2);
  
  gsl_multimin_fminimizer_free(myMinimizer3);
  gsl_vector_free(guess);
  gsl_vector_free(step);
  
  MU = calc_mu(SIGMA,M,RHO2);
  CalibStatus = FITSUCCESS;
  
  /*
    cout << "************ Sigma=" << SIGMA << endl;
        cout << "************ M=" << M << endl;
        cout << "************ Rho2=" << RHO2 << endl;
        cout << "************ Mu=" << MU << endl;

        // Print the ADC1 spectrum
        for (int x=RANGE_LOW; x<=THRESHOLD; x+=1) {
                cout << x << '\t' << GetEmpiricalDist(1,x)/N1 << '\t' << prob(x,1) << endl;
        }

        // Print the ADC2 spectrum
        for (int x=0; x<=RANGE_HIGH; x+=1) {
                cout << x << '\t' << GetEmpiricalDist(2,x)/N2 << '\t' << prob(x,2) << endl;
        }
        */
}
int CoarseScaleCalib::GetEmpiricalDist ( int  scale,
int  adc 
)

Definition at line 482 of file CoarseScaleCalib.cc.

                                                         {
        if (scale==1) return (int) count(adc1.begin(), adc1.end(), adc);
        else if (scale==2) return (int) count(adc2.begin(), adc2.end(), adc);
        else throw("Invalid scale passed to CoarseScaleCalib::GetEmpiricalDist()");
}
double CoarseScaleCalib::prob ( double  x,
int  scale 
)

Definition at line 488 of file CoarseScaleCalib.cc.

                                                 {
        if (CalibStatus != FITSUCCESS) throw("No calibration fit");
        if (scale != 1 && scale != 2) throw("Not a valid scale number");

        if (scale==1) {
                if (x<RANGE_LOW || x>=THRESHOLD) return 0;
                return q1_unnormalized(x, MU, SIGMA, PED1, PEDWID1) / q1_integral(THRESHOLD, MU, SIGMA, PED1, PEDWID1);
        }
        if (scale==2) {
                if (x<0 || x>M*RANGE_HIGH) return 0;
                return q2_unnormalized(x, MU, SIGMA, M, PED1, PED2, RHO2) / q2_integral(RANGE_HIGH, MU, SIGMA, M, PED1, PED2, RHO2);
        }
        return -1; // Should never get here
}
void CoarseScaleCalib::Get_quality ( double &  quality1,
double &  quality2 
)

Definition at line 503 of file CoarseScaleCalib.cc.

                                                                     {
        if (CalibStatus != FITSUCCESS) throw("No calibration fit");

        double chisq=0;
        double dof = 0;
        double e,o;

        vector<double> hist1(int(THRESHOLD-RANGE_LOW),0);
        for (unsigned int i=0; i<adc1.size(); i++) {
                hist1.at(int(adc1.at(i)-RANGE_LOW)) += 1;
        }
        for (unsigned int j=0; j<hist1.size(); j++) {
                e = N1*prob(RANGE_LOW+j,1);
                o = hist1.at(j);
                if (e>1e-6) {
                        chisq += (e-o)*(e-o)/e;
                        dof += 1;
                }
        }
        quality1 = chisq/dof;

        chisq = 0; dof=0;
        vector<double> hist2(int(RANGE_HIGH+1),0);
        for (unsigned int i=0; i<adc2.size(); i++) {
                hist2.at(int(adc2.at(i))) += 1;
        }
        for (unsigned int j=0; j<hist2.size(); j++) {
                e = N2*prob(j,2);
                o = hist2.at(j);
                if (e>1e-6) {
                        chisq += (e-o)*(e-o)/e;
                        dof += 1;
                }
        }
        quality2 = chisq/dof;
}
double CoarseScaleCalib::Get_mu ( )

Definition at line 540 of file CoarseScaleCalib.cc.

                                {
        if (CalibStatus != FITSUCCESS) throw("No calibration fit");
        return MU;
}
double CoarseScaleCalib::Get_sigma ( )

Definition at line 545 of file CoarseScaleCalib.cc.

                                   {
        if (CalibStatus != FITSUCCESS) throw("No calibration fit");
        return SIGMA;
}
double CoarseScaleCalib::Get_rho2 ( )

Definition at line 550 of file CoarseScaleCalib.cc.

                                  {
        if (CalibStatus != FITSUCCESS) throw("No calibration fit");
        return RHO2;
}
double CoarseScaleCalib::Get_M ( )

Definition at line 555 of file CoarseScaleCalib.cc.

                               {
        if (CalibStatus != FITSUCCESS) throw("No calibration fit");
        return M;
}
double CoarseScaleCalib::Get_N1 ( )

Definition at line 560 of file CoarseScaleCalib.cc.

                                {
        return N1;
}
double CoarseScaleCalib::Get_N2 ( )

Definition at line 564 of file CoarseScaleCalib.cc.

                                {
        return N2;
}
double CoarseScaleCalib::Get_ADC2_Smearing ( )

Definition at line 568 of file CoarseScaleCalib.cc.

                                           {
        if (CalibStatus != FITSUCCESS) throw("No calibration fit");
        if (RHO2 <= PEDWID2) return 0;
        return sqrt(RHO2*RHO2-PEDWID2*PEDWID2);
}
double CoarseScaleCalib::q1_unnormalized ( double  x,
double  mu,
double  sigma,
double  eta1,
double  rho1 
) [static, private]

Definition at line 65 of file CoarseScaleCalib.cc.

                                                                                                    {

        double sigP1 = sqrt(sigma*sigma + rho1*rho1);
        double sigT1 = rho1*(sigP1/sigma);
        double muP1  = mu - eta1;
        double muT1  = eta1 + (rho1/sigma)*(rho1/sigma) - (sigP1/sigma)*(sigP1/sigma)*THRESHOLD;

        double f1 = gsl_ran_gaussian_pdf(x-muP1, sigP1);
        double f2 = gsl_cdf_gaussian_P(-x-muT1, sigT1);

        return f1*f2;
}
double CoarseScaleCalib::p_q1_unnormalized ( double  x,
void *  parameters 
) [static, private]

Definition at line 77 of file CoarseScaleCalib.cc.

                                                                     {
        double* p = (double*) parameters;
        return q1_unnormalized(x,p[0],p[1],p[2],p[3]);
}
double CoarseScaleCalib::q1_integral ( double  x,
double  mu,
double  sigma,
double  eta1,
double  rho1 
) [static, private]

Definition at line 83 of file CoarseScaleCalib.cc.

                                                                                                {

        // Speed up and stabilize the algorithm by shortcutting the calculation of the constraint for extreme values of x
        if (x<=RANGE_LOW) return 0;
        double n = (x-mu+eta1)/sqrt(sigma*sigma+rho1*rho1);
        if (n<-5) return 0;
        if (n>5) return q1_integral(4.95*sqrt(sigma*sigma+rho1*rho1)+mu-eta1, mu, sigma, eta1, rho1);

        double params[] = { mu, sigma, eta1, rho1 };

        gsl_function F;
        F.function = &p_q1_unnormalized;
        F.params = &params;

        double result;
        double abserr;

        gsl_integration_workspace* myIntegrator = gsl_integration_workspace_alloc(1000000);
        gsl_integration_qag(&F, RANGE_LOW, x, 1e-8, 0, 1000000, 6, myIntegrator, &result, &abserr);
        gsl_integration_workspace_free(myIntegrator);

        return result;
}
double CoarseScaleCalib::q2_unnormalized ( double  x,
double  mu,
double  sigma,
double  m,
double  eta1,
double  eta2,
double  rho2 
) [static, private]

Definition at line 108 of file CoarseScaleCalib.cc.

                                                                                                                           {

        double sigP2 = sqrt((sigma/m)*(sigma/m) + rho2*rho2);
        double sigT2 = m*m*sigma*sigP2*rho2;
        double muP2  = (mu - eta1)/m;
        double muT2  = (mu + m*eta2 - eta1)*m*m*rho2*rho2 - (THRESHOLD + m*eta2 - eta1)*m*m*sigP2*sigP2;

        double f1 = gsl_ran_gaussian_pdf(x-muP2, sigP2);
        double f2 = 1.0-gsl_cdf_gaussian_P(-m*sigma*sigma*(x+eta2)-muT2,sigT2);

        return f1*f2;
}
double CoarseScaleCalib::p_q2_unnormalized ( double  x,
void *  parameters 
) [static, private]

Definition at line 120 of file CoarseScaleCalib.cc.

                                                                     {
        double* p = (double*) parameters;
        return q2_unnormalized(x,p[0],p[1],p[2],p[3],p[4],p[5]);
}
double CoarseScaleCalib::q2_integral ( double  x,
double  mu,
double  sigma,
double  m,
double  eta1,
double  eta2,
double  rho2 
) [static, private]

Definition at line 126 of file CoarseScaleCalib.cc.

                                                                                                                       {

        // Speed up and stabilize the algorithm by shortcutting the calculation of the constraint for extreme values of x
        if (x<=0) return 0;
        double n = (x-(mu-eta1)/m)/sqrt(sigma*sigma+rho2*rho2);
        if (n<-5) return 0;
        if (n>5) return q2_integral(4.95*sqrt(sigma*sigma+rho2*rho2)+(mu-eta1)/m, mu, sigma, m, eta1, eta2, rho2);

        double params[] = { mu, sigma, m, eta1, eta2, rho2 };

        gsl_function F;
        F.function = &p_q2_unnormalized;
        F.params = &params;

        double result;
        double abserr;

        gsl_integration_workspace* myIntegrator = gsl_integration_workspace_alloc(1000000);
        gsl_integration_qag(&F, 0, x, 1e-8, 0, 1000000, 6, myIntegrator, &result, &abserr);
        gsl_integration_workspace_free(myIntegrator);

        return result;
}
double CoarseScaleCalib::constraint ( double  x,
void *  parameter 
) [static, private]

Definition at line 151 of file CoarseScaleCalib.cc.

                                                               {
        double* p = (double*) parameters;
        double sigma = p[0];
        double m     = p[1];
        double eta1  = p[2];
        double eta2  = p[3];
        double rho1  = p[4];
        double rho2  = p[5];
        double ratio = p[6];

        double num = q1_integral(THRESHOLD,    mu, sigma, eta1, rho1);
        double den = q2_integral(RANGE_HIGH, mu, sigma, m, eta1, eta2, rho2);

        if (num==0 || den==0) {
                if (mu < THRESHOLD) return 1e6;
                if (mu >= THRESHOLD) return -1e6;
        }

        return num/den - ratio;
}
double CoarseScaleCalib::calc_mu ( double  sigma,
double  m,
double  rho2 
) [private]

Definition at line 173 of file CoarseScaleCalib.cc.

                                                                    {

        double params[] = {sigma, m, PED1, PED2, PEDWID1, rho2, N1/N2};
        gsl_function f;
        f.function = &constraint;
        f.params = &params;

        const gsl_root_fsolver_type* root_finding_alg = gsl_root_fsolver_bisection;
        gsl_root_fsolver* mySolver = gsl_root_fsolver_alloc(root_finding_alg);
        gsl_root_fsolver_set(mySolver, &f, 1, 10000);

        int status;
        int iter = 0, max_iter = 1000;
        double root, x_lo, x_hi;
        do {
                iter++;
                status = gsl_root_fsolver_iterate(mySolver);

                root = gsl_root_fsolver_root(mySolver);
                x_lo = gsl_root_fsolver_x_lower(mySolver);
                x_hi = gsl_root_fsolver_x_upper(mySolver);
                status = gsl_root_test_interval (x_lo, x_hi, 0.001, 0);

        } while (status == GSL_CONTINUE && iter < max_iter);

        gsl_root_fsolver_free(mySolver);

        return root;
}
double CoarseScaleCalib::adc1_likelihood ( double  sigma) [private]

Definition at line 204 of file CoarseScaleCalib.cc.

                                                     {

        double mu    = calc_mu(sigma,M,RHO2); // Calculate mu from the constraint equation based on the current values of M and RHO2

        double sigP1 = sqrt(sigma*sigma + PEDWID1*PEDWID1);
        double sigT1 = PEDWID1*(sigP1/sigma);
        double muP1  = mu - PED1;
        double muT1  = PED1 + (PEDWID1/sigma)*(PEDWID1/sigma) - (sigP1/sigma)*(sigP1/sigma)*THRESHOLD;

        double term1 = -N1*log(sigP1);
        double term2 = 0;
        for (int i=0; i<N1; i++) {
                term2 -= (adc1.at(i)-muP1)*(adc1.at(i)-muP1);
        }
        term2 /= 2*sigP1*sigP1;
        double term3 = 0;
        for (int i=0; i<N1; i++) {
                term3 += log( gsl_cdf_gaussian_P(-adc1.at(i)-muT1, sigT1) );
        }
        double term4 = -q1_integral(THRESHOLD, mu, sigma, PED1, PEDWID1);

        double likelihood = term1+term2+term3+term4;
        if (likelihood != likelihood || likelihood < -numeric_limits<double>::max() || likelihood > numeric_limits<double>::max()) {
                return 1e6; // likelihood is nan, inf, or -inf
        }
        return -likelihood;
}
double CoarseScaleCalib::p_adc1_likelihood ( double  x,
void *  parameters 
) [static, private]

Definition at line 231 of file CoarseScaleCalib.cc.

                                                                     {
        CoarseScaleCalib** p = (CoarseScaleCalib**) parameters;
        return p[0]->adc1_likelihood(x);
}
double CoarseScaleCalib::adc2_likelihood ( double  m,
double  rho2 
) [private]

Definition at line 237 of file CoarseScaleCalib.cc.

                                                              {

        double mu    = calc_mu(SIGMA,m,rho2); // Calculate mu from the constraint equation based on the current values of SIGMA

        double sigP2 = sqrt((SIGMA/m)*(SIGMA/m) + rho2*rho2);
        double sigT2 = m*m*SIGMA*sigP2*rho2;
        double muP2  = (mu - PED1)/m;
        double muT2  = (mu + m*PED2 - PED1)*m*m*rho2*rho2 - (THRESHOLD + m*PED2 - PED1)*m*m*sigP2*sigP2;

        double term1 = -N2*log(sigP2);
        double term2 = 0;
        for (int i=0; i<N2; i++) {
                term2 -= (adc2.at(i)-muP2)*(adc2.at(i)-muP2);
        }
        term2 /= 2*sigP2*sigP2;
        double term3 = 0;
        for (int i=0; i<N2; i++) {
                term3 += log( 1-gsl_cdf_gaussian_P(-m*SIGMA*SIGMA*(adc2.at(i)+PED2)-muT2, sigT2) );
        }
        double term4 = -q2_integral(RANGE_HIGH, mu, SIGMA, m, PED1, PED2, rho2);

        double likelihood = term1+term2+term3+term4;
        if (likelihood != likelihood || likelihood < -numeric_limits<double>::max() || likelihood > numeric_limits<double>::max()) {
                return 1e6; // likelihood is nan, inf, or -inf
        }
        return -likelihood;
}
double CoarseScaleCalib::p_adc2_likelihood ( const gsl_vector *  v,
void *  parameters 
) [static, private]

Definition at line 264 of file CoarseScaleCalib.cc.

                                                                                {
        CoarseScaleCalib** p = (CoarseScaleCalib**) parameters;
        double m = gsl_vector_get(v,0);
        double rho2 = gsl_vector_get(v,1);
        return p[0]->adc2_likelihood(m,rho2);
}
double CoarseScaleCalib::likelihood ( double  sigma,
double  m,
double  rho2 
) [private]

Definition at line 272 of file CoarseScaleCalib.cc.

                                                                       {

        double mu = calc_mu(sigma,m,rho2);

        // ADC1 part
        double sigP1 = sqrt(sigma*sigma + PEDWID1*PEDWID1);
        double sigT1 = PEDWID1*sigma/sigP1;
        double muP1  = mu - PED1;
        double muT1;

        double term1 = -N1*log(sigP1);
        double term2 = 0;
        for (int i=0; i<N1; i++) {
                term2 -= (adc1.at(i)-muP1)*(adc1.at(i)-muP1);
        }
        term2 /= 2*sigP1*sigP1;
        double term3 = 0;
        for (int i=0; i<N1; i++) {
                muT1  = mu*(sigT1/sigma)*(sigT1/sigma) + (adc1.at(i)+PED1)*(sigma/sigP1)*(sigma/sigP1);
                term3 += log( gsl_cdf_gaussian_P(THRESHOLD-muT1, sigT1) );
        }
        double term4 = -q1_integral(THRESHOLD, mu, sigma, PED1, PEDWID1);

        double likelihood1 = term1+term2+term3+term4;

        // ADC2 part
        double sigP2 = sqrt((sigma/m)*(sigma/m) + rho2*rho2);
        double sigT2 = m*m*sigma*sigP2*rho2;
        double muP2  = (mu - PED1)/m;
        double muT2  = (mu + m*PED2 - PED1)*m*m*rho2*rho2 - (THRESHOLD + m*PED2 - PED1)*m*m*sigP2*sigP2;

        double term1a = -N2*log(sigP2);
        double term2a = 0;
        for (int i=0; i<N2; i++) {
                term2a -= (adc2.at(i)-muP2)*(adc2.at(i)-muP2);
        }
        term2a /= 2*sigP2*sigP2;
        double term3a = 0;
        for (int i=0; i<N2; i++) {
                term3a += log( 1-gsl_cdf_gaussian_P(-m*sigma*sigma*(adc2.at(i)+PED2)-muT2, sigT2) );
        }
        double term4a = -q2_integral(RANGE_HIGH, mu, sigma, m, PED1, PED2, rho2);

        double likelihood2 = term1a+term2a+term3a+term4a;

        // Total (negative) likelihood
        double likelihood = likelihood1+likelihood2;
        if (likelihood != likelihood || likelihood < -numeric_limits<double>::max() || likelihood > numeric_limits<double>::max()) {
                return 1e6; // likelihood is nan, inf, or -inf
        }
        return -likelihood;
}
double CoarseScaleCalib::p_likelihood ( const gsl_vector *  v,
void *  parameters 
) [static, private]

Definition at line 324 of file CoarseScaleCalib.cc.

                                                                           {
        CoarseScaleCalib** p = (CoarseScaleCalib**) parameters;
        double sigma = gsl_vector_get(v,0);
        double m = gsl_vector_get(v,1);
        double rho2 = gsl_vector_get(v,2);

        return p[0]->likelihood(sigma,m,rho2);
}

Member Data Documentation

vector<double> CoarseScaleCalib::adc1 [private]

Definition at line 50 of file CoarseScaleCalib.h.

vector<double> CoarseScaleCalib::adc2 [private]

Definition at line 50 of file CoarseScaleCalib.h.

double CoarseScaleCalib::N [private]

Definition at line 53 of file CoarseScaleCalib.h.

double CoarseScaleCalib::N1 [private]

Definition at line 53 of file CoarseScaleCalib.h.

double CoarseScaleCalib::N2 [private]

Definition at line 53 of file CoarseScaleCalib.h.

double CoarseScaleCalib::PED1 [private]

Definition at line 56 of file CoarseScaleCalib.h.

double CoarseScaleCalib::PED2 [private]

Definition at line 56 of file CoarseScaleCalib.h.

double CoarseScaleCalib::PEDWID1 [private]

Definition at line 56 of file CoarseScaleCalib.h.

double CoarseScaleCalib::PEDWID2 [private]

Definition at line 56 of file CoarseScaleCalib.h.

double CoarseScaleCalib::M [private]

Definition at line 59 of file CoarseScaleCalib.h.

double CoarseScaleCalib::MU [private]

Definition at line 59 of file CoarseScaleCalib.h.

double CoarseScaleCalib::SIGMA [private]

Definition at line 59 of file CoarseScaleCalib.h.

double CoarseScaleCalib::RHO2 [private]

Definition at line 59 of file CoarseScaleCalib.h.

Definition at line 60 of file CoarseScaleCalib.h.


The documentation for this class was generated from the following files:
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:14:01 for CoarseScaleCalibrationAlg by doxygen 1.7.4