KV2Body.h

/***************************************************************************
                          KV2Body.h  -  description
                             -------------------
    begin                : 28/11/2003
    copyright            : (C) 2002 by J.D. Frankland
    email                : frankland@ganil.fr
 
$Id: KV2Body.h,v 1.5 2009/02/02 13:52:29 ebonnet Exp $
 ***************************************************************************/
 
/***************************************************************************
 *                                                                         *
 *   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.                                   *
 *                                                                         *
 ***************************************************************************/
 
#ifndef KV2BODY_H
#define KV2BODY_H
 
#include "KVBase.h"
#include <vector>
#include <optional>
#include "TVector3.h"
#include "KVNucleus.h"
#include "TF1.h"
#include "TGraph.h"
#include "KVDrawable.h"
 
 
 
class KV2Body: public KVBase {
 
    std::vector<KVNucleus> fNuclei;
    Double_t fEDiss;
 
    TVector3 VCM;
    Double_t BCM;
    Double_t WLT;
    Double_t WCT;
    Double_t WC[5];
 
    Double_t VC[5];
    Double_t EC[5];
    Double_t K[5];
    Double_t TETAMAX[5];
    Double_t TETAMIN[5];
 
    mutable TF1* fKoxReactionXSec;
    mutable TF1* fEqbmChargeState;
    mutable TF1* fEqbmChargeStateShSol;
    mutable TF1* fEqbmChargeStateShGas;
 
   void Set4thNucleus();
   Double_t ThetaLabVsThetaCM(Double_t*, Double_t*);
   Double_t ELabVsThetaCM(Double_t*, Double_t*);
   Double_t ELabVsThetaLab(Double_t*, Double_t*);
   Double_t XSecRuthLab(Double_t*, Double_t*);
   Double_t XSecRuthLabInt(Double_t*, Double_t*);
   Double_t XSecRuthCM(Double_t*, Double_t*);
   Double_t XSecRuthCMVsThetaCM(Double_t*, Double_t*);
 
   mutable TF1* fThetaLabVsThetaCM[5];
   mutable TF1* fELabVsThetaCM[5];
   mutable TF1* fELabVsThetaLab[5];
 
   mutable TF1* fXSecRuthLabIntegral[5];
   mutable TF1* fXSecRuthLab[5];
 
    Bool_t fSetOutgoing;
 
    Double_t fIntPrec;
 
   struct KinematicSolutions {
      std::array<std::optional<double>, 2> solutions;
      void iterate_solutions(std::function<void(double)> f)
      {
         for (auto& s : solutions)
            if (s)
               f(*s);
      }
      double val(int i)
      {
         return solutions[i].value_or(0.);
      }
      std::optional<double> first() const
      {
         return solutions[0];
      }
      std::optional<double> second() const
      {
         return solutions[1];
      }
   };
   using KinematicSolutions_v = std::optional<KinematicSolutions>;
   KinematicSolutions_v FindRoots(TF1*, Double_t) const;
 
public:
   enum nucleus_of_interest {
      projectile = 1,
      target,
      projectile_like,
      target_like
   };
   enum kinematic_solution {
      high_E_branch,
      low_E_branch,
      both_branches
   };
 
   void init();
   KV2Body();
   KV2Body(const Char_t* systemname);
   KV2Body(const KVNucleus& compound, double Exx = 0.0);
   KV2Body(const KVNucleus& proj, const KVNucleus& targ, double Ediss = 0.0);
   KV2Body(const KVNucleus& proj, const KVNucleus& targ, const KVNucleus& outgoing, double Ediss = 0.0);
   KV2Body(KVNucleus*, KVNucleus* = nullptr, KVNucleus* = nullptr, Double_t = 0.0);
   virtual ~ KV2Body();
 
   void CalculateKinematics();
   TString GetReactionEquation() const;
 
   static Double_t GetVelocity(Double_t mass, Double_t E);
 
   void SetProjectile(const KVNucleus&);
   void SetProjectile(Int_t z, Int_t a = 0);
   void SetTarget(const KVNucleus&);
   void SetTarget(Int_t z, Int_t a = 0);
   void SetOutgoing(const KVNucleus& proj_out);
 
   void SetExcitEnergy(Double_t ex)
   {
      fEDiss = ex;
   }
   Double_t GetExcitEnergy() const
   {
      return fEDiss;
   }
   void SetEDiss(Double_t ex)
   {
      SetExcitEnergy(ex);
   }
   Double_t GetEDiss() const
   {
      return GetExcitEnergy();
   }
   KVNucleus* GetNucleus(nucleus_of_interest i) const;
 
   Double_t GetQReaction() const;
   Double_t GetQGroundStates() const;
   Double_t GetLabGrazingAngle(nucleus_of_interest i = projectile) const;
 
   Double_t GetCMEnergy() const;
   Double_t GetCMEnergy(nucleus_of_interest i) const;
 
   Double_t GetMaxAngleLab(nucleus_of_interest i) const;
   Double_t GetMinAngleLab(nucleus_of_interest i) const;
 
   TVector3 GetCMVelocity() const;
   Double_t GetCMVelocity(nucleus_of_interest i) const;
   Double_t GetCMGamma() const
   {
      Double_t gamma =
         1.0 - TMath::Power(BCM, 2.);
      if (gamma > 0.)
         gamma = 1. / TMath::Sqrt(gamma);
      else {
         KVError::Warning(this, "GetCMGamma", "1 - (beta)**2 = %f ... strange!", gamma);
         gamma = 0.0;
      }
      return gamma;
   }
 
   KVDrawable<TF1> GetThetaLabVsThetaCMFunc(nucleus_of_interest OfNucleus) const;
   KVDrawable<TF1> GetELabVsThetaCMFunc(nucleus_of_interest OfNucleus) const;
   KVDrawable<TF1> GetELabVsThetaLabFunc(nucleus_of_interest OfNucleus, kinematic_solution kine_sol = kinematic_solution::high_E_branch) const;
 
   Double_t GetThetaLab(Double_t ThetaCM, nucleus_of_interest OfNucleus) const
   {
      return const_cast<KV2Body*>(this)->GetThetaLabVsThetaCMFunc(OfNucleus)->Eval(ThetaCM);
   }
   Double_t GetELab(Double_t ThetaCM, nucleus_of_interest OfNucleus) const
   {
      return const_cast<KV2Body*>(this)->GetELabVsThetaCMFunc(OfNucleus)->Eval(ThetaCM);
   }
   KinematicSolutions_v GetThetaCM(Double_t ThetaLab, nucleus_of_interest OfNucleus) const;
   std::optional<double> GetThetaCM(Double_t ThetaLab, nucleus_of_interest OfNucleus, kinematic_solution kine_sol) const;
   Double_t GetThetaCM(nucleus_of_interest OfNucleus, Double_t theta, nucleus_of_interest OtherNucleus) const
   {
      if (TMath::Abs(OfNucleus - OtherNucleus) % 2) return 180. - theta;
      return theta;
   }
 
   KVDrawable<TGraph> GetGraphELabVsThetaLab(nucleus_of_interest OfNucleus = projectile_like, int npoints = 50) const;
 
   KinematicSolutions_v GetELab(nucleus_of_interest OfNucleus, Double_t ThetaLab, std::optional<nucleus_of_interest> AngleNucleus = {}) const;
   KinematicSolutions_v GetVLab(nucleus_of_interest OfNucleus, Double_t ThetaLab, std::optional<nucleus_of_interest> AngleNucleus = {}) const;
   KinematicSolutions_v GetThetaLab(nucleus_of_interest OfNucleus, Double_t ThetaLab, nucleus_of_interest AngleNucleus) const;
 
   Double_t GetXSecRuthLab(Double_t ThetaLab, nucleus_of_interest OfNucleus = projectile_like, kinematic_solution kineSol = high_E_branch) const;
   Double_t GetXSecRuthCM(Double_t ThetaLab, nucleus_of_interest OfNucleus = projectile_like, kinematic_solution kineSol = high_E_branch) const;
 
   Double_t GetIntegratedXSecRuthLab(Float_t th1, Float_t th2, Float_t phi1 = -1, Float_t phi2 = -1, nucleus_of_interest OfNucleus = projectile_like);
 
   KVDrawable<TF1> GetXSecRuthCMFunc(nucleus_of_interest OfNucleus = projectile_like, Double_t theta_cm_min = 1., Double_t theta_cm_max = 179.) const;
   KVDrawable<TF1> GetXSecRuthLabFunc(nucleus_of_interest OfNucleus = projectile_like, kinematic_solution kineSol = high_E_branch, Double_t theta_min = 1., Double_t theta_max = 179.) const;
   KVDrawable<TF1> GetXSecRuthLabIntegralFunc(nucleus_of_interest OfNucleus = projectile_like, Double_t theta_min = 1., Double_t theta_max = 179., kinematic_solution kineSol = high_E_branch) const;
 
   void Print(Option_t* opt = "") const override;
 
   Double_t BassIntBarrier();
   Double_t KoxReactionXSec(Double_t*, Double_t*);
   KVDrawable<TF1> GetKoxReactionXSecFunc() const;
 
   Double_t GetSphereDureReactionXSec(Double_t r0 = 1.05);
   Double_t GetBmaxFromReactionXSec(Double_t ReacXsec);
   Double_t GetIntegratedXsec(Double_t b1, Double_t b2);
 
   Double_t EqbmChargeState(Double_t* t, Double_t*);
   KVDrawable<TF1> GetEqbmChargeStateFunc() const;
   Double_t eqbm_charge_state_shiwietz_solid(Double_t* t, Double_t*);
   KVDrawable<TF1> GetShiwietzEqbmChargeStateFuncForSolidTargets() const;
   Double_t eqbm_charge_state_shiwietz_gas(Double_t* t, Double_t*);
   KVDrawable<TF1> GetShiwietzEqbmChargeStateFuncForGasTargets() const;
 
   Double_t GetIntegralPrecision()
   {
      return fIntPrec;
   }
   void SetIntegralPrecision(Double_t precision)
   {
      fIntPrec = precision;
   }
 
   Bool_t IsInverseKinematics() const
   {
      return GetNucleus(projectile)->GetA() > GetNucleus(target)->GetA();
   }
 
   ClassDefOverride(KV2Body, 0)         //Relativistic binary kinematical calculation
};
 
#endif