KVWilckeReactionParameters.h

 
#ifndef __KVWILCKEREACTIONPARAMETERS_H
#define __KVWILCKEREACTIONPARAMETERS_H
 
#define THIRD (1./3.)
#include "KVNucleus.h"
#include "TF1.h"
#include "TH1.h"
#include "TPad.h"
#include <TF2.h>
 
class KVWilckeReactionParameters {
   /* The following variables are exactly copied from pp. 395-400 of Wilcke et al. */
    Int_t ZP, ZT;
    Int_t NP, NT;
    Int_t AP, AT;
   Int_t ZC;
   Int_t NC;
   Int_t AC;
 
   Double_t MU;
   Double_t RINT;
   Double_t R0;
    Double_t RP, RT;
    Double_t CP, CT;
   Double_t CBAR;
    Double_t RCP, RCT, RCTOTAL;
    static Double_t e2_Wilcke;
    static Double_t mu_Wilcke;
    Double_t V0;
   Double_t BSS_K, BSS_N;
    Double_t VC_RINT;
    Double_t FISSIONTKE;
    Double_t ASYMMFISSIONTKE;
    Double_t GAMMA;
    Double_t PROXFACTOR;
    Double_t RBARRIER;
    Double_t LCRIT;
    Double_t VRB;
 
   /* Numerical constants from W.D. Myers, Phys. Lett. B 30, 451 (1969) */
    static Double_t a1_Myers;
    static Double_t a2_Myers;
    static Double_t J_Myers;
    static Double_t K_Myers;
    static Double_t L_Myers;
    static Double_t M_Myers;
    static Double_t Q_Myers;
    static Double_t c1_Myers;
    static Double_t r0_Myers;
 
    TF1* fBSS;
    TF1* fProx;
    TF1* fPotential;
    TF1* fCMThetaQuart;
    TF1* fLmax;
    TF1* fSigmaR;
    TF1* fSigmaFus;
   TF1* fCentPot;
 
   void init();
 
public:
   KVWilckeReactionParameters();
   KVWilckeReactionParameters(const KVNucleus& proj, const KVNucleus& targ);
   virtual ~KVWilckeReactionParameters();
 
   void SetEntranceChannel(const KVNucleus& proj, const KVNucleus& targ);
 
   static Double_t InteractionRadius(Int_t aproj, Int_t atarg)
   {
      return MatterHalfDensityRadius(atarg) + MatterHalfDensityRadius(aproj) + 4.49 -
             (MatterHalfDensityRadius(atarg) + MatterHalfDensityRadius(aproj)) / 6.35;
   }
   static Double_t r0_Wilcke(Int_t aproj, Int_t atarg)
   {
      return InteractionRadius(aproj, atarg) / (pow(aproj, THIRD) + pow(atarg, THIRD));
   }
   static Double_t SharpRadius(Int_t A)
   {
      return (1.28 * pow(A, THIRD) - 0.76 + 0.8 * pow(A, -THIRD));
   }
   static Double_t MatterHalfDensityRadius(Int_t A)
   {
      return (SharpRadius(A) * (1. - pow(SharpRadius(A), -2.)));
   }
   static Double_t epsilon_bar_Myers(Int_t Z, Int_t A);
   static Double_t delta_bar_Myers(Int_t Z, Int_t A);
   static Double_t ChargeRadius_Myers(Int_t Z, Int_t A)
   {
      return r0_Myers * pow(2.*Z / (1. - 3.*epsilon_bar_Myers(Z, A)) / (1. - delta_bar_Myers(Z, A)), THIRD);
   }
   static Double_t BSS_V0(Int_t zp, Int_t ap, Int_t zt, Int_t at)
   {
      Double_t v0  = pow((zp + zt), 2.) / pow(pow(ChargeRadius_Myers(zt, at), 3.) + pow(ChargeRadius_Myers(zp, ap), 3.), THIRD);
      v0 -= (pow(zt, 2.) / ChargeRadius_Myers(zt, at) + pow(zp, 2.) / ChargeRadius_Myers(zp, ap));
      v0 *= 3.*e2_Wilcke / 5.;
      return v0;
   }
   Double_t VC(Double_t* r, Double_t*)
   {
      if (*r >= RCTOTAL) {
         return (ZP * ZT * e2_Wilcke / (*r));
      }
      if (*r == 0) return V0;
      return (V0 - BSS_K * pow(*r, BSS_N));
   }
   Double_t ProxPot(Double_t* r, Double_t*)
   {
      static Double_t zeta1 = 1.2511;
      static Double_t zeta0 = 2.54;
      static Double_t Kprox = 0.0852;
 
      Double_t Phi;
    Double_t zeta = *r - CP - CT;
      Double_t dzeta = zeta - zeta0;
      if (zeta <= zeta1) {
         Phi = -0.5 * pow(dzeta, 2) - Kprox * pow(dzeta, 3);
      }
      else {
         Phi = -3.437 * exp(-zeta / 0.75);
      }
      return PROXFACTOR * Phi;
   }
   Double_t Potential(Double_t* r, Double_t*)
   {
      return ProxPot(r, 0) + VC(r, 0);
   }
   Double_t CentrifugalPotential(Double_t* x, Double_t* l)
   {
 
      Double_t R = TMath::Max(0.1, x[0]);
      Double_t Vcent = l[0] * (l[0] + 1.) * pow(KVNucleus::hbar, 2) / (2.*MU * mu_Wilcke * R * R);
      return ProxPot(&x[0], 0) + VC(&x[0], 0) + Vcent;
   }
   Double_t GetAngularMomentumFromImpactParameter(Double_t e_sur_a, Double_t b) const
   {
      return b * k(e_sur_a);
   }
   Double_t GetImpactParameterFromAngularMomentum(Double_t e_sur_a, Double_t l) const
   {
      return l / k(e_sur_a);
   }
   Double_t GetCrossSectionFromMaxAngularMomentum(Double_t e_sur_a, Double_t lmax) const
   {
      return 10.*(TMath::Pi() / pow(k(e_sur_a), 2)) * pow(lmax + 0.5, 2);
   }
   Double_t GetMaxAngularMomentumFromCrossSection(Double_t e_sur_a, Double_t sigma) const
   {
      return sqrt(sigma / (10.*(TMath::Pi() / pow(k(e_sur_a), 2)))) - 0.5;
   }
   Double_t GetCrossSectionFromMaxImpactParameter(Double_t bmax) const
   {
      return 10.*TMath::Pi() * pow(bmax, 2);
   }
   Double_t GetMaxImpactParameterFromCrossSection(Double_t sigma) const
   {
      return sqrt(sigma / (10.*TMath::Pi()));
   }
   TF1* GetCentrifugalPotential(Double_t e_sur_a, Double_t b) const
   {
      Double_t l = GetAngularMomentumFromImpactParameter(e_sur_a, b);
      fCentPot->SetParameter(0, l);
      return fCentPot;
   }
 
   TF1* GetCentrifugalPotential(Double_t l) const
   {
      fCentPot->SetParameter(0, l);
      return fCentPot;
   }
   TF1* GetBSSCoulombPotential() const
   {
      return fBSS;
   }
   TF1* GetNuclearProximityPotential() const
   {
      return fProx;
   }
   TF1* GetTotalPotential() const
   {
      return fPotential;
   }
   static Double_t TKESymFiss(Int_t Z, Int_t A)
   {
      return (0.1071 * Z * Z / pow(A, THIRD) + 22.3);
   }
   Double_t GetFullDampingTKE() const
   {
      return ASYMMFISSIONTKE;
   }
   static Double_t NLDSurfaceTensionCoefficient(Int_t Z, Int_t A)
   {
      Double_t I = (A - 2.*Z) / (1.*A);
      return 0.9517 * (1. - 1.7826 * I * I);
   }
   static Double_t SWaveFusionBarrierRadius(Int_t zp, Int_t ap, Int_t zt, Int_t at)
   {
      Double_t zpzt = zp * zt;
      Double_t D;
      if (zpzt < 500) {
         D = 0.3117 * pow(zpzt, 0.2122);
      }
      else {
         D = 1.096 + 1.391e-04 * zpzt;
      }
      return InteractionRadius(ap, at) - D;
   }
   static Double_t RLDCriticalAngularMomentum(Int_t z, Int_t a);
   Double_t Eta(Double_t e_sur_a) const
   {
      return 0.15746 * ZP * ZT / pow(e_sur_a, 0.5);
   }
   Double_t k(Double_t e_sur_a) const
   {
      return 0.2187 * AT * AP * pow(e_sur_a, 0.5) / (1.*AC);
   }
   Double_t QuarterPointAngle(Double_t* x, Double_t*) const
   {
      Double_t krint = k(*x) * RINT;
      Double_t eta = Eta(*x);
      if (krint < 2.*eta) return TMath::Pi();
      return 2.*asin(eta / (krint - eta));
   }
   Double_t ProjectileLabQP(Double_t e) const
   {
      Double_t QP = QuarterPointAngle(&e, 0);
      return atan(sin(QP) / (cos(QP) + AP / (1.*AT)));
   }
   Double_t ProjectileLabEQP(Double_t e) const
   {
      return AP * e * (AP * AP + 2.*AP * AT * cos(QuarterPointAngle(&e, 0)) + AT * AT) / (1.*AC * AC);
   }
 
   Double_t Lmax(Double_t* x, Double_t*) const
   {
      return Eta(*x) / tan(QuarterPointAngle(x, 0) / 2.);
   }
   TF1* GetCMQuarterPointAngle() const
   {
      return fCMThetaQuart;
   }
   TF1* GetLmax() const
   {
      return fLmax;
   }
   Double_t SigmaR(Double_t* x, Double_t*) const
   {
      return GetCrossSectionFromMaxAngularMomentum(*x, Lmax(x, 0));
   }
   Double_t ECM(Double_t e_sur_a) const
   {
      return MU * e_sur_a;
   }
   TF1* GetReactionCrossSection() const
   {
      return fSigmaR;
   }
   Double_t SigmaFus(Double_t* e_sur_a, Double_t*) const
   {
      if (*e_sur_a <= 0) return 0.;
      Double_t e = VRB / ECM(*e_sur_a);
      Double_t S1 = (e > 1. ? 0. : GetCrossSectionFromMaxImpactParameter(RBARRIER) * (1. - e));
      Double_t S2 = GetCrossSectionFromMaxAngularMomentum(*e_sur_a, LCRIT);
      Double_t S = TMath::Min(S1, S2);
      return S;
   }
   TF1* GetFusionCrossSection() const
   {
      return fSigmaFus;
   }
   void DrawAllPotentials(Double_t l = 0) const
   {
      TF1* totpot;
      if (l > 0) {
         totpot = GetCentrifugalPotential(l);
         totpot->SetTitle(Form("HIPOT l=%3.0f", l));
      }
      else totpot = GetTotalPotential();
      totpot->SetNpx(1000);
      GetBSSCoulombPotential()->SetNpx(1000);
      GetNuclearProximityPotential()->SetNpx(1000);
      totpot->SetLineColor(kBlack);
      totpot->Draw();
      GetBSSCoulombPotential()->SetLineColor(kRed);
      GetBSSCoulombPotential()->Draw("same");
      Double_t max = TMath::Max(GetBSSCoulombPotential()->GetMaximum() * 2., totpot->GetMinimum() * 10);
      GetNuclearProximityPotential()->SetLineColor(kBlue);
      GetNuclearProximityPotential()->Draw("same");
      Double_t min = GetNuclearProximityPotential()->GetMinimum() * 1.2;
      totpot->GetHistogram()->GetYaxis()->SetRangeUser(min, max);
      ((TPad*)gPad)->BuildLegend();
   }
 
   void Print() const;
   Double_t PotentialPocketRadius(Double_t l);
   Double_t GetMaximumAngularMomentumWithPocket();
   Double_t PotentialMaximumRadius(Double_t l);
 
   Double_t GetMu() const
   {
      return MU;
   }
   Int_t GetAP() const
   {
      return AP;
   }
   Int_t GetAT() const
   {
      return AT;
   }
   Int_t GetAC() const
   {
      return AC;
   }
   Double_t GetCP() const
   {
      return CP;
   }
   Double_t GetCT() const
   {
      return CT;
   }
   Double_t GetRint() const
   {
      return RINT;
   }
   Double_t GetBassReactionCrossSection(Double_t e_sur_a);
 
   ClassDef(KVWilckeReactionParameters, 1) //Reaction parameters for heavy-ion collisions (Wilcke et al)
};
 
#endif