New cumac class added.

DanielWielanek · Nov 13, 2024 · 4496791 · 4496791
1 parent 9dddb61
commit 4496791
Show file tree

Hide file tree

Showing 5 changed files with 72 additions and 6 deletions.
diff --git a/analysis/femto/HalFemtoLinkDef.h b/analysis/femto/HalFemtoLinkDef.h
@@ -84,6 +84,7 @@
 #pragma link C++ class Hal::CorrFit1DCFCumacK0K0 + ;
 #pragma link C++ class Hal::CorrFit1DCFCumacK0Kch + ;
 #pragma link C++ class Hal::CorrFit1DCFCumacDLam + ;
+#pragma link C++ class Hal::CorrFit1DCFCumacStrong + ;
 
 
 #pragma link C++ class Hal::CorrFitParamsSetup + ;

diff --git a/analysis/femto/HalFemtoLinkDef2.h b/analysis/femto/HalFemtoLinkDef2.h
@@ -81,6 +81,7 @@
 #pragma link C++ class Hal::CorrFit1DCFCumacK0K0 + ;
 #pragma link C++ class Hal::CorrFit1DCFCumacK0Kch + ;
 #pragma link C++ class Hal::CorrFit1DCFCumacDLam + ;
+#pragma link C++ class Hal::CorrFit1DCFCumacStrong + ;
 
 #pragma link C++ class Hal::CorrFitParamsSetup + ;
 #pragma link C++ class Hal::CorrFit3DCFMultiDim + ;

diff --git a/analysis/femto/corrfit/CorrFit1DCF.h b/analysis/femto/corrfit/CorrFit1DCF.h
@@ -30,11 +30,10 @@ namespace Hal {
 
   private:
     Double_t EvalDenominator(Double_t x) const;
-    Int_t fRinvParIndex   = {1};
-    Int_t fLambdaParIndex = {2};
-
 
   protected:
+    Int_t fRinvParIndex   = {1};
+    Int_t fLambdaParIndex = {2};
     /**
      * pointer of acutally calculated bin during chi-square computation
      */

diff --git a/analysis/femto/corrfit/fittingfunctions/CorrFit1DCFCumac.cxx b/analysis/femto/corrfit/fittingfunctions/CorrFit1DCFCumac.cxx
@@ -25,6 +25,9 @@ namespace Hal {
     fZ = z;
     return Integr();
   }
+
+  Double_t CorrFit1DCFCumac::F(Double_t d, Double_t r0) const { return 1.0 - d / (2.0 * fPis * r0); }
+
   Double_t CorrFit1DCFCumac::F2(Double_t z) const { return (1.0 - TMath::Exp(-z * z)) / z; }
 
   Double_t CorrFit1DCFCumac::Ff1(Double_t x) const { return TMath::Exp(x * x - fZ * fZ) / fZ; }
@@ -495,8 +498,6 @@ namespace Hal {
     FixParameter(QuartetEffectiveRadiusID(), 3.22);
   }
 
-  Double_t CorrFit1DCFCumacDLam::F(Double_t d, Double_t r0) const { return 1.0 - d / (2.0 * fPis * r0); }
-
 
   Double_t CorrFit1DCFCumac::Get(Double_t q, Double_t r) {
     Double_t val[1];
@@ -508,4 +509,50 @@ namespace Hal {
     return CalculateCF(val, params);
   }
 
+  Double_t CorrFit1DCFCumacStrong::CalculateCF(const Double_t* x, const Double_t* params) const {
+    Double_t ak = 0;
+    switch (fKinematics) {
+      case Femto::EKinematics::kPRF: ak = x[0]; break;
+      case Femto::EKinematics::kLCMS: ak = x[0] * 0.5; break;
+      default: ak = 0.0; break;
+    }
+    Double_t r0    = Femto::FmToGeV(params[RadiusID()]);
+    Double_t q     = 2.0 * ak;
+    auto component = [&](double f0, double d0) {
+      f0 = Femto::FmToGeV(f0);
+      d0 = Femto::FmToGeV(d0);
+      // TComplex tmp = TComplex(d0 * ak * ak * 0.5, -ak) - 1.0 / f0;  // https://arxiv.org/pdf/2005.05012
+      TComplex tmp = TComplex(d0 * ak * ak * 0.5, -ak) + 1.0 / f0;  // Yu Hu
+      auto fk      = TComplex(1.0, 0) / tmp;
+      return (fk * fk).Rho() * F(d0, r0) / (2.0 * r0 * r0) + 2.0 * fk.Re() / (fPis * r0) * F1(q * r0) - fk.Im() / r0 * F2(q * r0);
+    };
+    double cf = 0;
+    for (int iSpin = 0; iSpin < fSpinStates; iSpin++) {
+      cf += params[SpinWeightID(iSpin)] * component(params[F0ID(iSpin)], params[D0ID(iSpin)]);
+    }
+    return params[NormID()] * (1.0 + params[LambdaID()] * cf);
+  }
+
+  CorrFit1DCFCumacStrong::CorrFit1DCFCumacStrong(Int_t spinStates) :
+    CorrFit1DCFCumac(3 + spinStates * 3), fSpinStates(spinStates) {
+    fNormParIndex   = 0;
+    fLambdaParIndex = 1;
+    fRinvParIndex   = 2;
+    SetParameterName(LambdaID(), "#lambda");
+    SetParameterName(NormID(), "N");
+    SetParameterName(RadiusID(), "R");
+    int start_par = 3;
+    for (int i = 0; i < fSpinStates; i++) {
+      SetParameterName(start_par++, Form("f0_%i", i));
+      SetParameterName(start_par++, Form("d0_%i", i));
+      SetParameterName(start_par++, Form("frac_%i", i));
+    }
+  }
+
+  Int_t CorrFit1DCFCumacStrong::SpinWeightID(Int_t spin) const { return 3 + spin * 3 + 2; }
+
+  Int_t CorrFit1DCFCumacStrong::D0ID(Int_t spin) const { return 3 + spin * 3 + 1; }
+
+  Int_t CorrFit1DCFCumacStrong::F0ID(Int_t spin) const { return 3 + spin * 3; }
+
 }  // namespace Hal
diff --git a/analysis/femto/corrfit/fittingfunctions/CorrFit1DCFCumac.h b/analysis/femto/corrfit/fittingfunctions/CorrFit1DCFCumac.h
@@ -37,6 +37,7 @@ namespace Hal {
     Simps(Double_t a1, Double_t b1, Double_t h1, Double_t reps1, Double_t aeps1, Double_t x, Double_t aih, Double_t aiabs) const;
     Double_t F1(Double_t Z) const;
     Double_t F2(Double_t Z) const;
+    Double_t F(Double_t d, Double_t r0) const;  // for d-lam
     CorrFit1DCFCumac(Int_t params);
 
   public:
@@ -94,7 +95,6 @@ namespace Hal {
    */
   class CorrFit1DCFCumacDLam : public CorrFit1DCFCumac {
   protected:
-    Double_t F(Double_t d, Double_t r0) const;
     virtual Double_t CalculateCF(const Double_t* x, const Double_t* params) const;
 
   public:
@@ -134,6 +134,7 @@ namespace Hal {
     virtual ~CorrFit1DCFCumacK0K0() {};
     ClassDef(CorrFit1DCFCumacK0K0, 1)
   };
+
   class CorrFit1DCFCumacK0Kch : public CorrFit1DCFCumac {
   protected:
     Hal::Femto::CorrFitGammaCalc fGammaCalc;
@@ -154,5 +155,22 @@ namespace Hal {
     virtual ~CorrFit1DCFCumacK0Kch() {};
     ClassDef(CorrFit1DCFCumacK0Kch, 1)
   };
+
+  class CorrFit1DCFCumacStrong : public CorrFit1DCFCumac {
+    const Int_t fSpinStates;
+
+  protected:
+    virtual Double_t CalculateCF(const Double_t* x, const Double_t* params) const;
+
+  public:
+    CorrFit1DCFCumacStrong(Int_t spinStates);
+    Int_t SpinWeightID(Int_t spin) const;
+    Int_t D0ID(Int_t spin) const;
+    Int_t F0ID(Int_t spin) const;
+    virtual ~CorrFit1DCFCumacStrong() {}
+    ClassDef(CorrFit1DCFCumacStrong, 1)
+  };
+
+
 }  // namespace Hal
 #endif /* HALCORRFIT1DCFCUMAC_H_ */