125 lines
4.1 KiB
C++
125 lines
4.1 KiB
C++
/***************************************************************************
|
|
|
|
PGbGLF.cpp
|
|
|
|
Author: Andreas Suter
|
|
e-mail: andreas.suter@psi.ch
|
|
|
|
***************************************************************************/
|
|
|
|
/***************************************************************************
|
|
* Copyright (C) 2007-2023 by Andreas Suter *
|
|
* andreas.suter@psi.ch *
|
|
* *
|
|
* This program is free software; you can redistribute it and/or modify *
|
|
* it under the terms of the GNU General Public License as published by *
|
|
* the Free Software Foundation; either version 2 of the License, or *
|
|
* (at your option) any later version. *
|
|
* *
|
|
* This program is distributed in the hope that it will be useful, *
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
|
* GNU General Public License for more details. *
|
|
* *
|
|
* You should have received a copy of the GNU General Public License *
|
|
* along with this program; if not, write to the *
|
|
* Free Software Foundation, Inc., *
|
|
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
|
|
***************************************************************************/
|
|
|
|
#ifdef HAVE_CONFIG_H
|
|
#include "config.h"
|
|
#endif
|
|
|
|
#include <cassert>
|
|
#include <iostream>
|
|
#include <cmath>
|
|
|
|
#include "PMusr.h"
|
|
#include "PGbGLF.h"
|
|
|
|
#define TWO_PI 6.28318530717958647692528676656
|
|
|
|
ClassImp(PGbGLF)
|
|
|
|
//--------------------------------------------------------------------------
|
|
// operator()
|
|
//--------------------------------------------------------------------------
|
|
/**
|
|
* <p>Method returning the function value at a given time for a given set of parameters.
|
|
*
|
|
* \param t time
|
|
* \param par
|
|
*/
|
|
Double_t PGbGLF::operator()(Double_t t, const std::vector<Double_t> &par) const
|
|
{
|
|
// param: [0] Bext (G), [1] Delta0 (1/us), [2] Rb (1)
|
|
assert(par.size()==3);
|
|
|
|
if (t<0.0)
|
|
return 1.0;
|
|
|
|
Double_t dval = 0.0;
|
|
Double_t wExt = TWO_PI * GAMMA_BAR_MUON * par[0];
|
|
Double_t s0 = par[1];
|
|
Double_t s1 = s0*par[2]; // sigma0 * Rb
|
|
|
|
Double_t wExt2 = wExt*wExt;
|
|
Double_t s02 = s0*s0;
|
|
Double_t s12 = s1*s1;
|
|
Double_t t2 = t*t;
|
|
|
|
if (par[0] < 3.0) { // take ZF solution for Bext <= 3G. This is good enough
|
|
Double_t aa = 1.0 + s12*t2;
|
|
dval = 0.333333333333333333 + 0.666666666666667 / pow(aa, 1.5)*(1.0-s02*t2/aa)*exp(-0.5*s02*t2/aa);
|
|
} else {
|
|
// P_z^LF (GbG, 1st part)
|
|
Double_t aa = 1.0+t2*s12;
|
|
dval = 1.0 - 2.0*(s02+s12)/wExt2 + 2.0*(s02+s12*aa)/(wExt2*pow(aa,2.5))*cos(wExt*t)*exp(-0.5*s02*t2/aa);
|
|
|
|
// P_z^LF (GbG, 2nd part)
|
|
Double_t dt = t;
|
|
Int_t n=1;
|
|
Double_t sumT = dt * pz_GbG_2(t, par) * 0.5;
|
|
Double_t sumM = 0.0;
|
|
Double_t tt = 0.0;
|
|
do {
|
|
sumM = 0.0;
|
|
for (Int_t i=0; i<n-1; i++) {
|
|
tt = (static_cast<Double_t>(i) + 0.5) * dt;
|
|
sumM += pz_GbG_2(tt, par);
|
|
}
|
|
sumM *= dt;
|
|
sumT = (sumT + sumM)*0.5;
|
|
dt /= 2.0;
|
|
n *= 2;
|
|
} while ((fabs(sumT-sumM) > 1.0e-5) && (n < 8192));
|
|
|
|
dval += sumT;
|
|
}
|
|
|
|
return dval;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------
|
|
// pz_GbG_2 (private)
|
|
//--------------------------------------------------------------------------
|
|
/**
|
|
* <p>Integrand of the non-analytic part
|
|
*
|
|
* \param t time
|
|
*/
|
|
Double_t PGbGLF::pz_GbG_2(Double_t t, const std::vector<Double_t> &par) const
|
|
{
|
|
Double_t wExt = TWO_PI * GAMMA_BAR_MUON * par[0];
|
|
Double_t s0 = par[1];
|
|
Double_t s1 = s0*par[2]; // sigma0 * Rb
|
|
|
|
Double_t s02 = s0*s0;
|
|
Double_t s12 = s1*s1;
|
|
Double_t t2 = t*t;
|
|
Double_t aa = 1.0+t2*s12;
|
|
|
|
return 2.0*(s02*s02+3.0*s12*s12*aa*aa+6.0*s02*s12*aa)/(pow(wExt,3.0)*pow(aa,4.5))*exp(-0.5*s02*t2/aa)*sin(wExt*t);
|
|
}
|