209#include "TRestAxionField.h"
213#include <gsl/gsl_integration.h>
251 Double_t lengthInMeters = Lcoh *
units(
"m");
255 Double_t sol = lengthInMeters * Bmag * tm;
269 Double_t lengthInMeters = Lcoh *
units(
"m");
273 Double_t sol = lengthInMeters * Bmag * tm / 2;
274 sol = sol * sol * 1.0e-20;
296 Double_t photonMass = mg;
300 RESTDebug <<
"+--------------------------------------------------------------------------+" << RESTendl;
301 RESTDebug <<
" TRestAxionField::GammaTransmissionProbability. Parameter summary" << RESTendl;
302 RESTDebug <<
" Photon mass : " << photonMass <<
" eV" << RESTendl;
303 RESTDebug <<
" Axion mass : " << ma <<
" eV" << RESTendl;
304 RESTDebug <<
" Axion energy : " <<
fEa <<
" keV" << RESTendl;
305 RESTDebug <<
" Lcoh : " <<
fLcoh <<
" mm" << RESTendl;
306 RESTDebug <<
" Bmag : " <<
fBmag <<
" T" << RESTendl;
307 RESTDebug <<
"+--------------------------------------------------------------------------+" << RESTendl;
311 Double_t q = (ma * ma - photonMass * photonMass) / 2. / (
fEa *
units(
"eV"));
313 Double_t phi = q * l;
315 Double_t Gamma = absLength;
317 Double_t GammaL = Gamma * cohLength *
units(
"cm") /
units(
"m");
320 std::cout <<
"+------------------------+" << std::endl;
321 std::cout <<
" Intermediate calculations" << std::endl;
322 std::cout <<
" q : " << q <<
" eV" << std::endl;
323 std::cout <<
" l : " << l <<
" eV-1" << std::endl;
324 std::cout <<
" phi : " << phi << std::endl;
325 std::cout <<
"Gamma : " << Gamma << std::endl;
326 std::cout <<
"GammaL : " << GammaL << std::endl;
327 std::cout <<
"+------------------------+" << std::endl;
330 Double_t MFactor = phi * phi + GammaL * GammaL / 4.0;
331 MFactor = 1.0 / MFactor;
334 std::cout <<
"Mfactor : " << MFactor << std::endl;
336 std::cout <<
"cos(phi) : " << cos(phi) << std::endl;
337 std::cout <<
"Exp(-GammaL) : " << exp(-GammaL) << std::endl;
341 (1 + exp(-GammaL) - 2 * exp(-GammaL / 2) * cos(phi)));
343 if (fDebug) std::cout <<
"Axion-photon transmission probability : " << sol << std::endl;
353 Double_t mg, Double_t absLength) {
381 Double_t Ea, Double_t ma, Double_t mg,
382 Double_t absLength) {
384 Double_t Lcoh = (Bmag.size() - 1) * deltaL;
385 Double_t cohLength = Lcoh *
units(
"m");
387 Double_t photonMass = mg;
391 Double_t fieldAverage = 0;
392 if (Bmag.size() > 0) fieldAverage = std::accumulate(Bmag.begin(), Bmag.end(), 0.0) / Bmag.size();
395 std::cout <<
"+--------------------------------------------------------------------------+"
397 std::cout <<
" TRestAxionField::GammaTransmissionProbability. Parameter summary" << std::endl;
398 std::cout <<
" Photon mass : " << photonMass <<
" eV" << std::endl;
399 std::cout <<
" Axion mass : " << ma <<
" eV" << std::endl;
400 std::cout <<
" Axion energy : " << Ea <<
" keV" << std::endl;
401 std::cout <<
" Lcoh : " << cohLength <<
" m" << std::endl;
402 std::cout <<
" Bmag average : " << fieldAverage <<
" T" << std::endl;
403 std::cout <<
"+--------------------------------------------------------------------------+"
408 if (ma == 0.0 && photonMass == 0.0)
return BLHalfSquared(fieldAverage, Lcoh);
410 Double_t q = (ma * ma - photonMass * photonMass) / 2. / (Ea *
units(
"eV"));
412 Double_t phi = q * l;
414 Double_t Gamma = absLength;
416 Double_t GammaL = Gamma * cohLength *
units(
"cm") /
units(
"m");
419 std::cout <<
"+------------------------+" << std::endl;
420 std::cout <<
" Intermediate calculations" << std::endl;
421 std::cout <<
" q : " << q <<
" eV" << std::endl;
422 std::cout <<
" l : " << l <<
" eV-1" << std::endl;
423 std::cout <<
" phi : " << phi << std::endl;
424 std::cout <<
"Gamma : " << Gamma << std::endl;
425 std::cout <<
"GammaL : " << GammaL << std::endl;
426 std::cout <<
"+------------------------+" << std::endl;
429 Double_t MFactor = phi * phi + GammaL * GammaL / 4.0;
430 MFactor = 1.0 / MFactor;
433 std::cout <<
"Mfactor : " << MFactor << std::endl;
434 std::cout <<
"(BL/2)^2 : " <<
BLHalfSquared(fieldAverage, Lcoh) << std::endl;
435 std::cout <<
"cos(phi) : " << cos(phi) << std::endl;
436 std::cout <<
"Exp(-GammaL) : " << exp(-GammaL) << std::endl;
442 for (
unsigned int n = 0; n < Bmag.size() - 1; n++) {
443 Double_t Bmiddle = 0.5 * (Bmag[n] + Bmag[n + 1]);
445 Double_t lStepIneV = ((double)n + 0.5) * deltaIneV;
446 Double_t lStepInCm = ((double)n + 0.5) * deltaL *
units(
"cm");
448 TComplex qCgC(0.5 * Gamma * lStepInCm, -q * lStepIneV);
449 qCgC = TComplex::Exp(qCgC);
451 TComplex integrand = Bmiddle * deltaL * qCgC;
456 Double_t sol = exp(-GammaL) * sum.Rho2() *
BLHalfSquared(1, 1);
459 if (fDebug) std::cout <<
"Axion-photon transmission probability : " << sol << std::endl;
461 return (Double_t)sol;
486 RESTError <<
"TRestAxionField::GammaTransmissionFieldMapProbability requires a magnetic field map!"
488 RESTError <<
"Use TRestAxionField::AssignMagneticField method to assign one" << RESTendl;
492 double photonMass = 0;
496 std::cout <<
"+--------------------------------------------------------------------------+"
498 std::cout <<
" TRestAxionField::GammaTransmissionProbability. Parameter summary" << std::endl;
499 std::cout <<
" Photon mass : " << photonMass <<
" eV" << std::endl;
500 std::cout <<
" Axion mass : " << ma <<
" eV" << std::endl;
501 std::cout <<
" Axion energy : " << Ea <<
" keV" << std::endl;
502 std::cout <<
"+--------------------------------------------------------------------------+"
506 double q = (ma * ma - photonMass * photonMass) / 2. / (Ea *
units(
"eV"));
513 std::cout <<
"+------------------------+" << std::endl;
514 std::cout <<
" Intermediate calculations" << std::endl;
515 std::cout <<
" q : " << q <<
" eV" << std::endl;
516 std::cout <<
"Gamma : " << Gamma << std::endl;
517 std::cout <<
"+------------------------+" << std::endl;
545 Int_t num_intervals) {
548 std::pair<TRestAxionMagneticField*, double> params = {
fMagneticField, Gamma};
550 gsl_integration_workspace* workspace = gsl_integration_workspace_alloc(num_intervals);
556 auto start = std::chrono::system_clock::now();
558 gsl_integration_qag(&F, 0,
fMagneticField->GetTrackLength(), accuracy, accuracy, num_intervals,
559 GSL_INTEG_GAUSS61, workspace, &reprob, &rerr);
561 auto end = std::chrono::system_clock::now();
562 auto seconds = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
567 double prob = C * reprob * reprob;
568 double proberr = 2 * C * reprob * rerr;
571 std::cout <<
" ---- TRestAxionField::ComputeResonanceIntegral (QAG) ----" << std::endl;
572 std::cout <<
"Gamma: " << Gamma <<
" mm-1" << std::endl;
573 std::cout <<
"accuracy: " << accuracy << std::endl;
574 std::cout <<
"num_intervals: " << num_intervals << std::endl;
575 std::cout <<
" -------" << std::endl;
576 std::cout <<
"Probability: " << prob << std::endl;
577 std::cout <<
"Probability error: " << proberr << std::endl;
578 std::cout <<
"Computing time: " << seconds.count() << std::endl;
579 std::cout <<
" -------" << std::endl;
582 std::pair<Double_t, Double_t> sol = {prob, proberr};
607 double improb, imerr;
609 std::pair<TRestAxionMagneticField*, double> params = {
fMagneticField, Gamma};
611 gsl_integration_workspace* workspace = gsl_integration_workspace_alloc(num_intervals);
617 auto start = std::chrono::system_clock::now();
619 gsl_integration_qawo_table* wf =
620 gsl_integration_qawo_table_alloc(q,
fMagneticField->GetTrackLength(), GSL_INTEG_COSINE, qawo_levels);
621 gsl_integration_qawo(&F, 0, accuracy, accuracy, num_intervals, workspace, wf, &reprob, &rerr);
623 gsl_integration_qawo_table_set(wf, q,
fMagneticField->GetTrackLength(), GSL_INTEG_SINE);
624 gsl_integration_qawo(&F, 0, accuracy, accuracy, num_intervals, workspace, wf, &improb, &imerr);
626 gsl_integration_qawo_table_free(wf);
628 auto end = std::chrono::system_clock::now();
629 auto seconds = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
634 double prob = C * (reprob * reprob + improb * improb);
635 double proberr = 2 * C * TMath::Sqrt(reprob * reprob * rerr * rerr + improb * improb * imerr * imerr);
638 std::cout <<
" ---- TRestAxionField::ComputeOffResonanceIntegral (QAWO) ----" << std::endl;
639 std::cout <<
"Gamma: " << Gamma <<
" mm-1" << std::endl;
640 std::cout <<
"q: " << q <<
"mm-1" << std::endl;
641 std::cout <<
"accuracy: " << accuracy << std::endl;
642 std::cout <<
"num_intervals: " << num_intervals << std::endl;
643 std::cout <<
"qawo_levels: " << qawo_levels << std::endl;
644 std::cout <<
" -------" << std::endl;
645 std::cout <<
"Probability: " << prob << std::endl;
646 std::cout <<
"Probability error: " << proberr << std::endl;
647 std::cout <<
"Computing time: " << seconds.count() << std::endl;
648 std::cout <<
" -------" << std::endl;
651 std::pair<Double_t, Double_t> sol = {prob, proberr};
671 Double_t photonMass = mg;
675 RESTDebug <<
"+--------------------------------------------------------------------------+"
677 RESTDebug <<
" TRestAxionField::AxionAbsorptionProbability. Parameter summary" << RESTendl;
678 RESTDebug <<
" Photon mass : " << photonMass <<
" eV" << RESTendl;
679 RESTDebug <<
" Axion mass : " << ma <<
" eV" << RESTendl;
680 RESTDebug <<
" Axion energy : " <<
fEa <<
" keV" << RESTendl;
681 RESTDebug <<
" Lcoh : " <<
fLcoh <<
" mm" << RESTendl;
682 RESTDebug <<
" Bmag : " <<
fBmag <<
" T" << RESTendl;
683 RESTDebug <<
"+--------------------------------------------------------------------------+"
689 Double_t q = (ma * ma - photonMass * photonMass) / 2. / (
fEa *
units(
"eV"));
691 Double_t phi = q * l;
693 Double_t Gamma = absLength;
695 Double_t GammaL = Gamma * cohLength *
units(
"cm") /
units(
"m");
698 RESTDebug <<
"+------------------------+" << RESTendl;
699 RESTDebug <<
" Intermediate calculations" << RESTendl;
700 RESTDebug <<
" q : " << q <<
" eV" << RESTendl;
701 RESTDebug <<
" l : " << l <<
" eV-1" << RESTendl;
702 RESTDebug <<
" phi : " << phi << RESTendl;
703 RESTDebug <<
"Gamma : " << Gamma << RESTendl;
704 RESTDebug <<
"GammaL : " << GammaL << RESTendl;
705 RESTDebug <<
"+------------------------+" << RESTendl;
708 Double_t MFactor = phi * phi + GammaL * GammaL / 4.0;
709 MFactor = 1.0 / MFactor;
712 RESTDebug <<
"Mfactor : " << MFactor << RESTendl;
714 RESTDebug <<
"cos(phi) : " << cos(phi) << RESTendl;
715 RESTDebug <<
"Exp(-GammaL) : " << exp(-GammaL) << RESTendl;
720 if (fDebug) RESTDebug <<
"Axion-photon absorption probability : " << sol << RESTendl;
730 Double_t mg, Double_t absLength) {
751 Double_t maxMass = 10;
753 Double_t resonanceMass = 0;
757 Double_t scanMass = resonanceMass;
761 if (scanMass > maxMass) {
762 RESTError <<
"TRestAxionField::GammaTransmissionProbability. Something went wrong when "
769 Double_t fwhm = scanMass - resonanceMass;
771 RESTError <<
"TRestAxionField::GammaTransmissionProbability. Problem calculating FWHM!" << RESTendl;
803 std::vector<std::pair<Double_t, Double_t>> massDensityPairs;
819 Double_t ma = firstMass;
823 massDensityPairs.push_back(std::make_pair(ma, density));
826 Double_t factor = TMath::Exp(-ma * rampDown) + 1;
832 massDensityPairs.push_back(std::make_pair(ma, density));
838 return massDensityPairs;
A metadata class to define the gas properties used in axion search calculations.
Double_t GetPhotonMass(double en)
It returns the equivalent photon mass (in eV) for the gas mixture at the given input energy expressed...
Double_t GetDensityForMass(double m_gamma, double en=4.2)
It returns the equivalent gas density for a given photon mass expressed in eV and a given axion energ...
void SetGasDensity(TString gasName, Double_t density)
It adds a new gas component to the mixture. If it already exists it will update its density.
Double_t GetPhotonAbsorptionLength(Double_t energy)
It returns the inverse of the absorption lenght, for the gas mixture, in cm-1, for the given energy i...
A basic class to define analytical axion-photon conversion calculations for axion helioscopes.
std::pair< Double_t, Double_t > GammaTransmissionFieldMapProbability(Double_t Ea, Double_t ma, Double_t accuracy=1.e-1, Int_t num_intervals=100, Int_t qawo_levels=20)
Performs the calculation of axion-photon conversion probability using directly equation (28) from J....
Double_t BL(Double_t Bmag, Double_t Lcoh)
Performs the calculation of (BL) factor in natural units.
Double_t AxionAbsorptionProbability(Double_t ma, Double_t mg=0, Double_t absLength=0)
Performs the calculation of axion-photon absorption probability using directly equation (18) from van...
void Initialize()
Initialization of TRestAxionField class.
static double Integrand(double x, void *params)
Integrand used for axion-photon probability integration.
Double_t GammaTransmissionProbability(Double_t ma, Double_t mg=0, Double_t absLength=0)
Performs the calculation of axion-photon conversion probability using directly equation (11) from van...
void AssignBufferGas(TRestAxionBufferGas *buffGas)
It assigns a gas buffer medium to the calculation.
~TRestAxionField()
Default destructor.
Double_t BLHalfSquared(Double_t Bmag, Double_t Lcoh)
Performs the calculation of (BL/2)^2 factor in natural units.
Double_t fEa
The energy of the axion in keV.
std::vector< std::pair< Double_t, Double_t > > GetMassDensityScanning(std::string gasName="He", Double_t maMax=0.15, Double_t rampDown=5.0)
This method determines the proper densities to be used in an axion helioscope experiment in order to ...
std::pair< Double_t, Double_t > ComputeOffResonanceIntegral(Double_t q, Double_t Gamma, Double_t accuracy, Int_t num_intervals, Int_t qawo_levels)
Performs the calculation of axion-photon conversion probability using directly equation (28) from J....
TRestAxionMagneticField * fMagneticField
A pointer to the magnetic field definition.
TRestAxionBufferGas * fBufferGas
A pointer to the buffer gas definition.
Double_t fBmag
The magnetic field in Teslas (used for constant field formulas)
Double_t fLcoh
The coherence lenght (in mm) where the magnetic field is defined (for constant field)
TRestAxionField()
Default constructor.
Double_t GammaTransmissionFWHM(Double_t step=0.00001)
Performs the calculation of the FWHM for the axion-photon conversion probability computed in TRestAxi...
std::pair< Double_t, Double_t > ComputeResonanceIntegral(Double_t Gamma, Double_t accuracy, Int_t num_intervals)
Performs the calculation of axion-photon conversion probability using directly equation (28) from J....
constexpr double lightSpeed
Speed of light in m/s.
constexpr double PhMeterIneV
A meter in eV.
constexpr double naturalElectron
Electron charge in natural units.