diff --git a/Scripts/Mathematica/CMs/InputParamsDoubleOmegaPhi.m b/Scripts/Mathematica/CMs/InputParamsDoubleOmegaPhi.m
new file mode 100644
index 0000000000000000000000000000000000000000..47f7bf36d04c7ce0b3cc066f7721ca801fd0cf97
--- /dev/null
+++ b/Scripts/Mathematica/CMs/InputParamsDoubleOmegaPhi.m
@@ -0,0 +1,69 @@
+(* Copyright 2025
+ Bertram Kopf (bertram@ep1.rub.de)
+ Ruhr-Universit<C3><A4>t Bochum
+
+ This file is part of Pawian.
+
+ Pawian 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 3 of the License, or
+ (at your option) any later version.
+
+ Pawian 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 Pawian. If not, see <http://www.gnu.org/licenses/>
+
+***********************************************************************************************
+ *)
+
+(* Precision[12/10]; *)
+
+(* masses of the decay products *)
+mpi:=0.13957;
+mpi0:=0.1349768;
+
+mKp:=0.493677;
+mKm:=0.493677;
+
+valImag := 0.00;
+
+(* Masses of the unstable resonance*)
+MBW1:=0.78266;
+MBW2:=1.019461;
+
+(*Widths of unstable resonance*)
+GammaBW1:=0.00868;
+GammaBW2:=0.004249;
+
+(* Masses of the decay particles *)
+mDec11:=2 mpi;
+mDec12:=mpi0;
+
+mDec21:=mKp;
+mDec22:=mKm;
+
+
+orbMomBW1 = 1;
+orbMomBW2 = 1;
+
+f2PoleBB1 := 0.02275031992477571786636110849745549994;
+m2PoleBB1 := 0.61238097137475222594286512618182349074;
+
+f2PoleBB2 := 0.67110752879828435237795796451561483007;
+m2PoleBB2 := 1.05187282169299207701471933416642866246;
+
+f2PoleDudek1 := 0.0225865133682695854973481933446531801;
+m2PoleDudek1 := 0.60741971755426057356854570645121124506;
+
+f2PoleDudek2 := 0.5596430172817635954719483569269226501;
+m2PoleDudek2 := 1.03873920225401985970974681359698101954;
+
+
+scalingFacBBSmearing:= 1. / 3.287318346901323;
+scalingFacDudekSmearing:=1. / 3.8866039897112796;
+scalingFacBWSmearing:= 1. / 0.7454720957505188;
+
diff --git a/Scripts/Mathematica/CMs/InputParamsDoubleSigmaSigma.m b/Scripts/Mathematica/CMs/InputParamsDoubleSigmaSigma.m
new file mode 100644
index 0000000000000000000000000000000000000000..d5d9b3224194a33d0308e20abac50e8cbdea23c1
--- /dev/null
+++ b/Scripts/Mathematica/CMs/InputParamsDoubleSigmaSigma.m
@@ -0,0 +1,66 @@
+(* Copyright 2025
+ Bertram Kopf (bertram@ep1.rub.de)
+ Ruhr-Universit<C3><A4>t Bochum
+
+ This file is part of Pawian.
+
+ Pawian 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 3 of the License, or
+ (at your option) any later version.
+
+ Pawian 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 Pawian. If not, see <http://www.gnu.org/licenses/>
+
+***********************************************************************************************
+ *)
+
+(* Precision[12/10]; *)
+
+(* masses of the decay products *)
+mpi:=0.13957;
+mpi0:=0.1349768;
+
+valImag := 0.00;
+
+(* Masses of the unstable resonance*)
+MBW1:=0.475;
+MBW2:=0.475;
+
+(*Widths of unstable resonance*)
+GammaBW1:=0.550;
+GammaBW2:=0.550;
+
+(* Masses of the decay particles *)
+mDec11:=mpi;
+mDec12:=mpi;
+
+mDec21:=mpi0;
+mDec22:=mpi0;
+
+
+orbMomBW1 = 0;
+orbMomBW2 = 0;
+
+f2PoleBB1 := 0.53525549369732492365844555006462551715;
+m2PoleBB1 := 0.02072166534172146370023379334682632026;
+
+f2PoleBB2 := 0.52537598909653228136099257201590318318;
+m2PoleBB2 := 0.01303711145733995864252490339790527583;
+
+f2PoleDudek1 := 0.53525549368142545098154277118417779544;
+m2PoleDudek1 := -0.32003256520041713522826607205769323751;
+
+f2PoleDudek2 := 0.52537598910691176848144924060406008584;
+m2PoleDudek2 := -0.32142763113897706864276934163696308259;
+
+
+scalingFacBBSmearing:= 1. / 1.1927781712395298;
+scalingFacDudekSmearing:=1. / 1.1927781669814008;
+scalingFacBWSmearing:= 1./ 1.1245037531410065;
+
diff --git a/Scripts/Mathematica/CMs/InputParamsKstar892Pi.m b/Scripts/Mathematica/CMs/InputParamsKstar892Pi.m
new file mode 100644
index 0000000000000000000000000000000000000000..c0dae7396f67bca60c876d005d5880343a129ae3
--- /dev/null
+++ b/Scripts/Mathematica/CMs/InputParamsKstar892Pi.m
@@ -0,0 +1,58 @@
+(* Copyright 2025
+ Bertram Kopf (bertram@ep1.rub.de)
+ Ruhr-Universit<C3><A4>t Bochum
+
+ This file is part of Pawian.
+
+ Pawian 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 3 of the License, or
+ (at your option) any later version.
+
+ Pawian 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 Pawian. If not, see <http://www.gnu.org/licenses/>
+
+***********************************************************************************************
+ *)
+
+(* Precision[12/10]; *)
+
+(* masses of the decay products *)
+mK:=0.49761
+mpi:=0.13957;
+mpi0:=0.1349768;
+
+valImag := 0.00;
+
+(**** parameters K*(892)0 pi0 -> (K- pi+) pi0 based on BB CM ****)
+(* Mass of the unstable resonance*)
+MBW:=0.8955;
+
+(*Width of unstable resonance*)
+GammaBW:=0.0514;
+
+(* Masses of the decay particles *)
+mDec1:=mK;
+mDec2:=mpi;
+
+(* Mass of the decay recoil particle *)
+mRec=mpi0;
+
+orbMomBW = 1;
+
+
+f2PoleBB := 0.31133809685087025051448922330820490963;
+m2PoleBB := 0.8238329330807920326010290727540961481;
+
+f2PoleDudek := 0.28759507372364078839913006596975963494;
+m2PoleDudek := 0.78197101092932291878039859817243872515;
+
+scalingFacBBSmearing:= 1. / 1.105361310555928;
+scalingFacDudekSmearing:=1. / 1.105361310555928;
+scalingFacBWSmearing:= 1.;
+(* scalingFacBWSmearing:= 1. / 0.08802808902167991; *)
diff --git a/Scripts/Mathematica/CMs/InputParamsPi1300Pi.m b/Scripts/Mathematica/CMs/InputParamsPi1300Pi.m
new file mode 100644
index 0000000000000000000000000000000000000000..7cb33edd0434cdd7938fcefc8cbbc935ce07f104
--- /dev/null
+++ b/Scripts/Mathematica/CMs/InputParamsPi1300Pi.m
@@ -0,0 +1,60 @@
+(* Copyright 2023
+ Bertram Kopf (bertram@ep1.rub.de)
+ Ruhr-Universit<C3><A4>t Bochum
+
+ This file is part of Pawian.
+
+ Pawian 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 3 of the License, or
+ (at your option) any later version.
+
+ Pawian 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 Pawian. If not, see <http://www.gnu.org/licenses/>
+
+***********************************************************************************************
+ *)
+
+(* Precision[12/10]; *)
+
+(* masses of the decay products *)
+mpi:=0.13957;
+mpi0:=0.1349768;
+
+valImag := 0.00;
+
+
+
+(**** CM parameters for pi(1300)0 pi0 w/ pi(1300)0 -> (pi+ pi0) pi- ****)
+
+(* Mass of the unstable resonance*)
+MBW:=1.300;
+
+(*Width of unstable resonance*)
+GammaBW:=0.400;
+
+(* Masses of the decay particles *)
+mDec1:=mpi+mpi0;
+mDec2:=mpi;
+
+(* Mass of the decay recoil particle *)
+mRec=mpi0;
+
+(* Momentum of the decay pi(1300) -> (pi+ pi0) pi- *)
+orbMomBW = 0;
+
+f2PoleBB := 0.62562280402090093702928976311097271904;
+m2PoleBB := 1.35581006495341548877557531562277255115;
+
+f2PoleDudek := 0.62562280398394301015550122301071064861;
+m2PoleDudek := 0.97194445334541372054996097989249217186;
+
+scalingFacBBSmearing:= 1.;
+scalingFacDudekSmearing:= 1.;
+scalingFacBWSmearing:= 1. / 0.8336427551137701;
+(* scalingFacBWSmearing:= 1. / 0.43630464471242414; *)
diff --git a/Scripts/Mathematica/CMs/InputParamsRho770Pi.m b/Scripts/Mathematica/CMs/InputParamsRho770Pi.m
new file mode 100644
index 0000000000000000000000000000000000000000..11fbfc0f1843c28c6f7083c8be39beb2b34c982f
--- /dev/null
+++ b/Scripts/Mathematica/CMs/InputParamsRho770Pi.m
@@ -0,0 +1,57 @@
+(* Copyright 2025
+ Bertram Kopf (bertram@ep1.rub.de)
+ Ruhr-Universit<C3><A4>t Bochum
+
+ This file is part of Pawian.
+
+ Pawian 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 3 of the License, or
+ (at your option) any later version.
+
+ Pawian 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 Pawian. If not, see <http://www.gnu.org/licenses/>
+
+***********************************************************************************************
+ *)
+
+(* Precision[12/10]; *)
+
+(* masses of the decay products *)
+mpi:=0.13957;
+mpi0:=0.1349768;
+
+valImag := 0.00;
+
+(**** parameters rho(700)+ pi- -> (pi+ pi0) pi- based on BB CM ****)
+(* Mass of the unstable resonance*)
+MBW:=0.77511;
+
+(*Width of unstable resonance*)
+GammaBW:=0.1491;
+
+(* Masses of the decay particles *)
+mDec1:=mpi+mpi0;
+mDec2:=mpi;
+
+(* Mass of the decay recoil particle *)
+mRec=mpi0;
+
+orbMomBW = 1;
+
+
+f2PoleBB := 0.4207051462764408283327900815775571806;
+m2PoleBB := 0.60605788493333204795276935269551359375;
+
+f2PoleDudek := 0.36966524147228358938226916212290303466;
+m2PoleDudek := 0.52303101506895045228388878328917860629;
+
+scalingFacBBSmearing:= 1./ 1.0655618168689487;
+scalingFacDudekSmearing:=1. / 1.1786373003089863;
+scalingFacBWSmearing:= 1. / 0.7697543609311899;
+(* scalingFacBWSmearing:= 1. / 0.08802808902167991; *)
diff --git a/Scripts/Mathematica/CMs/Integrals.m b/Scripts/Mathematica/CMs/Integrals.m
new file mode 100644
index 0000000000000000000000000000000000000000..36d1ecc6823ba2e3901cde8e98d9b26cf2271d7a
--- /dev/null
+++ b/Scripts/Mathematica/CMs/Integrals.m
@@ -0,0 +1,157 @@
+(* Copyright 2025
+ Bertram Kopf (bertram@ep1.rub.de)
+ Ruhr-Universit<C3><A4>t Bochum
+
+ This file is part of Pawian.
+
+ Pawian 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 3 of the License, or
+ (at your option) any later version.
+
+ Pawian 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 Pawian. If not, see <http://www.gnu.org/licenses/>
+
+***********************************************************************************************
+ *)
+
+Get[FileNameJoin[{Directory[], "../KinUtils.m"}]]
+
+(* Precision[1.2]; *)
+(* Precision[12/10]; *)
+
+dsDudek[sprime_,m1_,m2_,f2R_,m2R_, s0_, orbL_] := m2R - sprime + f2R BarrierSimple[sprime, (m1+m2)^2, s0, orbL] cmDudek[sprime,m1,m2];
+
+ImagSigmads2Dudek[sprime_,m1_,m2_,f2R_,m2R_, s0_, orbL_] := If[Abs[dsDudek[sprime,m1,m2,f2R,m2R, s0, orbL]]^2 < 10^30, - BarrierSimple[sprime, (m1+m2)^2, s0, orbL] Im[cmDudek[sprime, m1, m2]] / Abs[dsDudek[sprime,m1,m2,f2R,m2R, s0, orbL]]^2, 0.];
+
+RNIntDudek[m1_, m2_, f2R_, m2R_, s0_, orbL_, scalingFac_, smax_] :=
+ scalingFac /\[Pi] NIntegrate[ f2R ImagSigmads2Dudek[sb,m1,m2,f2R,m2R, s0, orbL]
+ , {sb, (m1 + m2)^2, smax},
+ Method -> {"GaussKronrodRule", "SymbolicProcessing" -> 0},
+ MaxRecursion -> 20,
+ AccuracyGoal -> 5];
+
+
+dsBB[sprime_,m1_,m2_,f2R_,m2R_, s0_, orbL_] := m2R - sprime - f2R BarrierSimple[sprime, (m1+m2)^2, s0, orbL] cmReid[sprime,m1,m2];
+
+ImagSigmads2BB[sprime_,m1_,m2_,f2R_,m2R_, s0_, orbL_] := If[Abs[dsBB[sprime,m1,m2,f2R,m2R, s0, orbL]]^2 < 10^30, BarrierSimple[sprime, (m1+m2)^2, s0, orbL] Im[cmReid[sprime, m1, m2]] / Abs[dsBB[sprime,m1,m2,f2R,m2R, s0, orbL]]^2, 0.];
+
+RNIntBB[m1_, m2_, f2R_, m2R_, s0_, orbL_, scalingFac_, smax_] :=
+ scalingFac /\[Pi] NIntegrate[ f2R ImagSigmads2BB[sb,m1,m2,f2R, m2R, s0, orbL]
+ , {sb, (m1 + m2)^2, smax},
+ Method -> {"GaussKronrodRule", "SymbolicProcessing" -> 0},
+ MaxRecursion -> 20,
+ AccuracyGoal -> 5];
+
+RNIntBW[M0_, Width_, m1_, m2_, orbMom_, scalingFac_, smax_] :=
+ scalingFac /\[Pi] NIntegrate[ rho[sb, m1, m2] Abs[BreitWignerBlattWRelNorm[sb, orbMom, M0, Width, m1, m2]]^2
+ , {sb, (m1 + m2)^2, smax},
+ Method -> {"GaussKronrodRule", "SymbolicProcessing" -> 0},
+ MaxRecursion -> 20,
+ AccuracyGoal -> 5];
+
+RNIntDudekviaBB[s_][m1_, m2_, mRecoil_, f2R_, m2R_, s0_, orbL_, scalingFac_] :=
+ scalingFac /\[Pi] NIntegrate[ f2R ImagSigmads2BB[sb,m1,m2,f2R, m2R, s0, orbL] cmDudek[s,Sqrt[sb], mRecoil]
+ , {sb, (m1 + m2)^2, Infinity},
+ Method -> {"GaussKronrodRule", "SymbolicProcessing" -> 0},
+ MaxRecursion -> 20,
+ AccuracyGoal -> 5];
+
+RNIntDudekviaDudek[s_][m1_, m2_, mRecoil_, f2R_, m2R_, s0_, orbL_, scalingFac_] :=
+ scalingFac /\[Pi] NIntegrate[ f2R ImagSigmads2Dudek[sb,m1,m2,f2R, m2R, s0, orbL] cmDudek[s,Sqrt[sb], mRecoil]
+ , {sb, (m1 + m2)^2, Infinity},
+ Method -> {"GaussKronrodRule", "SymbolicProcessing" -> 0},
+ MaxRecursion -> 20,
+ AccuracyGoal -> 5];
+
+RNIntDudekviaBW[s_][m1_, m2_, mrec_, orbMom_, M0_, Width_, scale_] :=
+ scale /\[Pi] NIntegrate[ rho[sb, m1, m2] Abs[BreitWignerBlattWRelNorm[sb, orbMom, M0, Width, m1, m2]]^2 cmDudek[s, Sqrt[sb], mrec]
+ , {sb, (m1 + m2)^2, Infinity},
+ Method -> {"GaussKronrodRule", "SymbolicProcessing" -> 0},
+ MaxRecursion -> 20,
+ AccuracyGoal -> 5];
+
+RNIntBBviaBB[s_][m1_, m2_, mRecoil_, f2R_, m2R_, s0_, orbL_, scalingFac_] :=
+ scalingFac /\[Pi] NIntegrate[ f2R ImagSigmads2BB[sb,m1,m2,f2R, m2R, s0, orbL] cmReid[s,Sqrt[sb], mRecoil]
+ , {sb, (m1 + m2)^2, Infinity},
+ Method -> {"GaussKronrodRule", "SymbolicProcessing" -> 0},
+ MaxRecursion -> 20,
+ AccuracyGoal -> 5];
+
+RNIntBBviaDudek[s_][m1_, m2_, mRecoil_, f2R_, m2R_, s0_, orbL_, scalingFac_] :=
+ scalingFac /\[Pi] NIntegrate[ f2R ImagSigmads2Dudek[sb,m1,m2,f2R, m2R, s0, orbL] cmReid[s,Sqrt[sb], mRecoil]
+ , {sb, (m1 + m2)^2, Infinity},
+ Method -> {"GaussKronrodRule", "SymbolicProcessing" -> 0},
+ MaxRecursion -> 20,
+ AccuracyGoal -> 5];
+
+RNIntBBviaBW[s_][m1_, m2_, mrec_, orbMom_, M0_, Width_, scale_] :=
+ scale /\[Pi] NIntegrate[ rho[sb, m1, m2] Abs[BreitWignerBlattWRelNorm[sb, orbMom, M0, Width, m1, m2]]^2 cmReid[s, Sqrt[sb], mrec]
+ , {sb, (m1 + m2)^2, Infinity},
+ Method -> {"GaussKronrodRule", "SymbolicProcessing" -> 0},
+ MaxRecursion -> 20,
+ AccuracyGoal -> 5];
+
+
+
+
+RNIntDudekDouble[m11_, m12_, f2R1_, m2R1_, orbL1_, s01_, m21_, m22_, f2R2_, m2R2_, orbL2_, s02_, scalingFac_, smax_] :=
+ scalingFac /\[Pi] NIntegrate[ f2R1 ImagSigmads2Dudek[sb1,m11,m12,f2R1,m2R1, s01, orbL1]
+ f2R2 ImagSigmads2Dudek[sb2,m21,m22,f2R2,m2R2, s02, orbL2]
+ , {sb1, (m11 + m12)^2, smax}, {sb2, (m21 + m22)^2, smax},
+ Method -> {"GaussKronrodRule", "SymbolicProcessing" -> 0},
+ MaxRecursion -> 20,
+ AccuracyGoal -> 5];
+
+RNIntDudekDoubleViaDudek[s_][m11_, m12_, f2R1_, m2R1_, orbL1_, s01_, m21_, m22_, f2R2_, m2R2_, orbL2_, s02_, scalingFac_] :=
+ scalingFac /\[Pi] NIntegrate[ f2R1 ImagSigmads2Dudek[sb1,m11,m12,f2R1,m2R1, s01, orbL1]
+ f2R2 ImagSigmads2Dudek[sb2,m21,m22,f2R2,m2R2, s02, orbL2]
+ cmDudek[s, Sqrt[sb1], Sqrt[sb2]]
+ , {sb1, (m11 + m12)^2, Infinity}, {sb2, (m21 + m22)^2, Infinity},
+ Method -> {"GaussKronrodRule", "SymbolicProcessing" -> 0},
+ MaxRecursion -> 20,
+ AccuracyGoal -> 5];
+
+
+
+
+RNIntBBDouble[m11_, m12_, f2R1_, m2R1_, orbL1_, s01_, m21_, m22_, f2R2_, m2R2_, orbL2_, s02_, scalingFac_, smax_] :=
+ scalingFac /\[Pi] NIntegrate[ f2R1 ImagSigmads2BB[sb1,m11,m12,f2R1,m2R1, s01, orbL1]
+ f2R2 ImagSigmads2BB[sb2,m21,m22,f2R2,m2R2, s02, orbL2]
+ , {sb1, (m11 + m12)^2, smax}, {sb2, (m21 + m22)^2, smax},
+ Method -> {"GaussKronrodRule", "SymbolicProcessing" -> 0},
+ MaxRecursion -> 20,
+ AccuracyGoal -> 5];
+
+
+RNIntDudekDoubleViaBB[s_][m11_, m12_, f2R1_, m2R1_, orbL1_, s01_, m21_, m22_, f2R2_, m2R2_, orbL2_, s02_, scalingFac_] :=
+ scalingFac /\[Pi] NIntegrate[ f2R1 ImagSigmads2BB[sb1,m11,m12,f2R1,m2R1, s01, orbL1]
+ f2R2 ImagSigmads2BB[sb2,m21,m22,f2R2,m2R2, s02, orbL2]
+ cmDudek[s, Sqrt[sb1], Sqrt[sb2]]
+ , {sb1, (m11 + m12)^2, Infinity}, {sb2, (m21 + m22)^2, Infinity},
+ Method -> {"GaussKronrodRule", "SymbolicProcessing" -> 0},
+ MaxRecursion -> 20,
+ AccuracyGoal -> 5];
+
+RNIntBWDouble[M01_, Width1_, m11_, m12_, orbMom1_, M02_, Width2_, m21_, m22_, orbMom2_,scalingFac_, smax_] :=
+ scalingFac /\[Pi] NIntegrate[ rho[sb1, m11, m12] Abs[BreitWignerBlattWRelNorm[sb1, orbMom1, M01, Width1, m11, m12]]^2
+ rho[sb2, m21, m22] Abs[BreitWignerBlattWRelNorm[sb2, orbMom2, M02, Width2, m21, m22]]^2
+ , {sb1, (m11 + m12)^2, smax}, {sb2, (m21 + m22)^2, smax},
+ Method -> {"GaussKronrodRule", "SymbolicProcessing" -> 0},
+ MaxRecursion -> 20,
+ AccuracyGoal -> 5];
+
+
+RNIntDudekDoubleViaBW[s_][M01_, Width1_, m11_, m12_, orbMom1_, M02_, Width2_, m21_, m22_, orbMom2_,scalingFac_] :=
+ scalingFac /\[Pi] NIntegrate[ rho[sb1, m11, m12] Abs[BreitWignerBlattWRelNorm[sb1, orbMom1, M01, Width1, m11, m12]]^2
+ rho[sb2, m21, m22] Abs[BreitWignerBlattWRelNorm[sb2, orbMom2, M02, Width2, m21, m22]]^2
+ cmDudek[s, Sqrt[sb1], Sqrt[sb2]]
+ , {sb1, (m11 + m12)^2, 10000.}, {sb2, (m21 + m22)^2, 10000.},
+ Method -> {"GaussKronrodRule", "SymbolicProcessing" -> 0},
+ MaxRecursion -> 20,
+ AccuracyGoal -> 5];
diff --git a/Scripts/Mathematica/CMs/KinUtils.m b/Scripts/Mathematica/CMs/KinUtils.m
new file mode 100644
index 0000000000000000000000000000000000000000..8ef6e028230c837a5742d024aa789140055033e8
--- /dev/null
+++ b/Scripts/Mathematica/CMs/KinUtils.m
@@ -0,0 +1,114 @@
+(* Copyright 2025
+ Bertram Kopf (bertram@ep1.rub.de)
+ Ruhr-Universität Bochum
+
+ This file is part of Pawian.
+
+ Pawian 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 3 of the License, or
+ (at your option) any later version.
+
+ Pawian 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 Pawian. If not, see <http://www.gnu.org/licenses/>
+
+***********************************************************************************************
+ This Mathematica-script is a colletion of functions which are needed to perform calculations with the K-matrix formalism.
+ *)
+
+(* phase space factor *)
+rho[s_,m1_,m2_]:= Sqrt[1-(m1+m2)^2/s] Sqrt[1-(m1-m2)^2/s];
+
+(* break up momentum *)
+BreakupMomQDefaultFromS[s_, s1_, s2_]:=
+ Module[{result},
+ result=rho[s, Sqrt[s1], Sqrt[s2]] Sqrt[s]/2.;
+ (* result = Sqrt[s*s + s1*s1 + s2*s2 - 2. * ( s*s1 + s*s2 + s1*s2)]/Sqrt[4.*s]; *)
+ result
+ ];
+
+(* CM Function based on Dudek et al Phys. Rev. D91, 054008 (2015), Appendix B; subtraction constant at s=s_thresh CM(s_thresh)=0 *)
+xi[s_,m1_,m2_]:=1-(m1+m2)^2/s;
+cmDudekPre[s_,m_,u_] := rho[s,m,u]/\[Pi] Log[(xi[s,m,u]+rho[s,m,u])/(xi[s,m,u]-rho[s,m,u])]-xi[s,m,u]/\[Pi] (u-m)/(u+m) Log[u/m];
+
+(* making sure to be in the 1 quadrant for s=real *)
+cmDudek[s_,m_,u_] := cmDudekPre[s + I * If[Im[s] == 0, 10.^-8, 0.], m, u];
+
+
+(* CM Function based on Reid with subtraction constant at s=0 CM(0)=0 *)
+(* cmReid[s_,m_,u_]:= - ( cmDudek[s,m,u] - 1/\[Pi] (1. + 2 m u / (u u - m m) Log[u/m] ) ); *)
+
+(* proteted case for u==m with L'Hospital *)
+cmReid[s_,m_,u_]:= - ( cmDudek[s,m,u] - If[m != u, 1/\[Pi] (1. + 2 m u / (u u - m m) Log[u/m]), 2./\[Pi]] );
+
+(* cmReid[s_,m_,u_]:= - ( cmDudek[s,m,u] - If[m != u, 1/\[Pi] (1. + 2 m u / (u u - m m) Log[u/m] ), 2./\[Pi] *)
+
+
+(* Blatt-Weisskopf factors *)
+BlattWeisskopf[l_, z_]:=
+ Module[{result},
+ result = Which[l==0, 1.0 + 0. I, l==1, Sqrt[2. z / (z+1.)] + 0. I, l==2, Sqrt[13. z z / (z z + 3. z + 9.)]+ 0. I];
+ result
+ ];
+
+BWBarrier[l_, q_, q0_]:=
+ Module[{result},
+ qR:=0.1973;
+ z:=(Norm[q] Norm[q])/(Norm[qR] Norm[qR]);
+ (* Print["z: ", z]; *)
+ z0:=(Norm[q0] Norm[q0])/(Norm[qR] Norm[qR]);
+ (* Print["z0: ", z0]; *)
+ result:= BlattWeisskopf[l,z]/BlattWeisskopf[l,z0];
+ result
+ ];
+
+(* elasticity *)
+Sii[s_, m1_, m2_, Tii_] := 1. + 2. I Sqrt[Re[rho[s,m1,m2]]] Tii Sqrt[Re[rho[s,m1,m2]]];
+
+
+(* extract delta phase of an Argand plot *)
+deltaArgand[mass_, m1_, m2_, Tii_]:=
+ Module[{delta},
+ Tiirho :=rho[mass*mass,m1,m2] Tii;
+ ReE := Re[Tiirho];
+ ImE := Im[Tiirho] - 0.5;
+ tmpdelta := ArcTan[Abs[ReE], ImE];
+ delta := 0.5 tmpdelta 180./Pi + 45.0;
+ Which[ReE < 0.0, delta = 180.0 - delta];
+ delta
+ ];
+
+(* Breit-Wigner function w/ Blatt-Weisskopf factors *)
+BreitWignerBlattWRel[sCurrent_, orbMom_, M0_, Width_, Ma_, Mb_]:=
+ Module[{BW},
+ s0=M0 M0;
+ rho0=rho[s0, Ma, Mb];
+ rhoCurrent=rho[sCurrent, Ma, Mb];
+ momQ0=BreakupMomQDefaultFromS[s0, Ma Ma, Mb Mb];
+ momQCurrent=BreakupMomQDefaultFromS[sCurrent, Ma Ma, Mb Mb];
+ bwBarrierFac=If[orbMom==0, 1., BWBarrier[orbMom, momQCurrent, momQ0]];
+ BW = M0 Width bwBarrierFac / (s0 - sCurrent - I (rhoCurrent/rho0) M0 Width bwBarrierFac bwBarrierFac);
+ BW
+ ];
+
+BreitWignerBlattWRelNorm[sCurrent_, orbMom_, M0_, Width_, Ma_, Mb_]:=
+ Module[{BW},
+ s0=M0 M0;
+ rho0=rho[s0, Ma, Mb];
+ rhoCurrent=rho[sCurrent, Ma, Mb];
+ momQ0=BreakupMomQDefaultFromS[s0, Ma Ma, Mb Mb];
+ momQCurrent=BreakupMomQDefaultFromS[sCurrent, Ma Ma, Mb Mb];
+ bwBarrierFac=If[orbMom==0, 1., BWBarrier[orbMom, momQCurrent, momQ0]];
+ GammaM= Width (rhoCurrent/rho0) bwBarrierFac bwBarrierFac;
+ BW = Sqrt[M0 GammaM] / (s0 - sCurrent - I M0 GammaM);
+ BW
+ ];
+
+(* sinplified barrier function based on s *)
+BarrierSimple[s_, sthresh_, s0_, orbMom_]:= If[orbMom==0, 1., ( (s-sthresh) / (s+s0) )^orbMom ];
+
diff --git a/Scripts/Mathematica/CMs/LUTsDouble/CMDudekDoubleViaBBLUT.m b/Scripts/Mathematica/CMs/LUTsDouble/CMDudekDoubleViaBBLUT.m
new file mode 100644
index 0000000000000000000000000000000000000000..0ec76c8ce728fd5216a4f2fedcbef311119df15a
--- /dev/null
+++ b/Scripts/Mathematica/CMs/LUTsDouble/CMDudekDoubleViaBBLUT.m
@@ -0,0 +1,88 @@
+(* Copyright 2025
+ Bertram Kopf (bertram@ep1.rub.de)
+ Ruhr-Universit<C3><A4>t Bochum
+
+ This file is part of Pawian.
+
+ Pawian 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 3 of the License, or
+ (at your option) any later version.
+
+ Pawian 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 Pawian. If not, see <http://www.gnu.org/licenses/>
+
+***********************************************************************************************
+ *)
+
+Get[FileNameJoin[{Directory[], "../Integrals.m"}]]
+Get[FileNameJoin[{Directory[], "../InputParamsDoubleSigmaSigma.m"}]]
+
+(* Precision[1.2]; *)
+(* Precision[12/10]; *)
+
+(*
+mrealstepsize:= 0.0005;
+mrealmin:= 0.005;
+mrealmax:= 3.5;
+
+mimagstepsize := 0.002;
+numOfCores := 1;
+mimagstepsizeNoCtimes := numOfCores mimagstepsize;
+mimagmin := -0.6;
+mimagmax := 0.002;
+ *)
+
+(*
+mrealstepsize:= 0.001;
+mrealmin:= 0.026;
+mrealmax:= 3.6;
+
+mimagstepsize := 0.01;
+numOfCores := 8;
+mimagstepsizeNoCtimes := numOfCores mimagstepsize;
+mimagmin := -0.6;
+mimagmax := 0.01;
+ *)
+
+mrealstepsize:= 0.01;
+mrealmin:= 0.026;
+mrealmax:= 3.6;
+
+mimagstepsize := 0.01;
+numOfCores := 8;
+mimagstepsizeNoCtimes := numOfCores mimagstepsize;
+mimagmin := -0.1;
+mimagmax := 0.01;
+
+(*
+{t1,b1} =
+ AbsoluteTiming[
+ *)
+For[ mreal=mrealmin, mreal<=mrealmax, mreal+=mrealstepsize,
+ Print["mreal: ", mreal];
+ ofile=OpenAppend["./CMDudekDoubleViaBBLUT.dat"];
+ For[ mimag = mimagmax, mimag >= mimagmin, mimag -= mimagstepsizeNoCtimes,{
+ SetSharedVariable[list];
+ list = {};
+ ParallelDo[{
+ currentCM = RNIntDudekDoubleViaBB[(mreal + I mimag) (mreal + I mimag)][mDec11, mDec12, f2PoleBB1, m2PoleBB1, orbMomBW1, (mDec11+mDec12)^2, mDec21, mDec22, f2PoleBB2, m2PoleBB2, orbMomBW2,(mDec21+mDec22)^2 , scalingFacBBSmearing];
+ AppendTo[list, {mimag, currentCM}]},
+ {mimag, mimag, mimag-mimagstepsizeNoCtimes+mimagstepsize, -mimagstepsize}
+ ],
+ list = SortBy[list, First, Greater];
+ WriteString[ofile, NumberForm[mreal,5], "\t", NumberForm[#[[1]], 5], "\t", NumberForm[Re[#[[2]]], 16], "\t", NumberForm[Im[#[[2]]], 16], "\n"] & /@ list;
+ };
+ ]
+ Close[ofile];
+ ];
+ (*
+];
+
+Print["t1: ", t1,"\tb1: ", b1];
+ *)
diff --git a/Scripts/Mathematica/CMs/LUTsDouble/CMDudekDoubleViaBWLUT.m b/Scripts/Mathematica/CMs/LUTsDouble/CMDudekDoubleViaBWLUT.m
new file mode 100644
index 0000000000000000000000000000000000000000..642bab463df4a821813b8eb2256fd9fd88dfe16d
--- /dev/null
+++ b/Scripts/Mathematica/CMs/LUTsDouble/CMDudekDoubleViaBWLUT.m
@@ -0,0 +1,89 @@
+(* Copyright 2025
+ Bertram Kopf (bertram@ep1.rub.de)
+ Ruhr-Universit<C3><A4>t Bochum
+
+ This file is part of Pawian.
+
+ Pawian 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 3 of the License, or
+ (at your option) any later version.
+
+ Pawian 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 Pawian. If not, see <http://www.gnu.org/licenses/>
+
+***********************************************************************************************
+ *)
+
+Get[FileNameJoin[{Directory[], "../Integrals.m"}]]
+Get[FileNameJoin[{Directory[], "../InputParamsDoubleSigmaSigma.m"}]]
+
+(* Precision[1.2]; *)
+(* Precision[12/10]; *)
+
+(*
+mrealstepsize:= 0.0005;
+mrealmin:= 0.005;
+mrealmax:= 3.5;
+
+mimagstepsize := 0.002;
+numOfCores := 1;
+mimagstepsizeNoCtimes := numOfCores mimagstepsize;
+mimagmin := -0.6;
+mimagmax := 0.002;
+ *)
+
+
+mrealstepsize:= 0.001;
+mrealmin:= 0.026;
+mrealmax:= 3.6;
+
+mimagstepsize := 0.01;
+numOfCores := 8;
+mimagstepsizeNoCtimes := numOfCores mimagstepsize;
+mimagmin := -0.6;
+mimagmax := 0.01;
+
+(*
+mrealstepsize:= 0.01;
+mrealmin:= 0.026;
+mrealmax:= 3.6;
+
+mimagstepsize := 0.01;
+numOfCores := 8;
+mimagstepsizeNoCtimes := numOfCores mimagstepsize;
+mimagmin := -0.1;
+mimagmax := 0.01;
+ *)
+
+(*
+{t1,b1} =
+ AbsoluteTiming[
+ *)
+For[ mreal=mrealmin, mreal<=mrealmax, mreal+=mrealstepsize,
+ Print["mreal: ", mreal];
+ ofile=OpenAppend["./CMDudekDoubleViaBWLUT.dat"];
+ For[ mimag = mimagmax, mimag >= mimagmin, mimag -= mimagstepsizeNoCtimes,{
+ SetSharedVariable[list];
+ list = {};
+ ParallelDo[{
+ currentCM = RNIntDudekDoubleViaBW[(mreal + I mimag) (mreal + I mimag)][MBW1, GammaBW1, mDec11, mDec12, orbMomBW1, MBW2, GammaBW2, mDec21, mDec22, orbMomBW2, scalingFacBWSmearing];
+ AppendTo[list, {mimag, currentCM}]},
+ {mimag, mimag, mimag-mimagstepsizeNoCtimes+mimagstepsize, -mimagstepsize}
+ ],
+ list = SortBy[list, First, Greater];
+ WriteString[ofile, NumberForm[mreal,5], "\t", NumberForm[#[[1]], 5], "\t", NumberForm[Re[#[[2]]], 16], "\t", NumberForm[Im[#[[2]]], 16], "\n"] & /@ list;
+ };
+ ]
+ Close[ofile];
+ ];
+ (*
+];
+
+Print["t1: ", t1,"\tb1: ", b1];
+ *)
diff --git a/Scripts/Mathematica/CMs/LUTsDouble/CMDudekDoubleViaDudekLUT.m b/Scripts/Mathematica/CMs/LUTsDouble/CMDudekDoubleViaDudekLUT.m
new file mode 100644
index 0000000000000000000000000000000000000000..2c745b6c924f633d8508d4c9aaf9e676340c476b
--- /dev/null
+++ b/Scripts/Mathematica/CMs/LUTsDouble/CMDudekDoubleViaDudekLUT.m
@@ -0,0 +1,67 @@
+(* Copyright 2025
+ Bertram Kopf (bertram@ep1.rub.de)
+ Ruhr-Universit<C3><A4>t Bochum
+
+ This file is part of Pawian.
+
+ Pawian 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 3 of the License, or
+ (at your option) any later version.
+
+ Pawian 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 Pawian. If not, see <http://www.gnu.org/licenses/>
+
+***********************************************************************************************
+ *)
+
+Get[FileNameJoin[{Directory[], "../Integrals.m"}]]
+(* Get[FileNameJoin[{Directory[], "../InputParamsDoubleSigmaSigma.m"}]] *)
+Get[FileNameJoin[{Directory[], "../InputParamsDoubleOmegaPhi.m"}]]
+
+(* Precision[1.2]; *)
+(* Precision[12/10]; *)
+
+
+
+mrealstepsize:= 0.0005;
+mrealmin:= 0.026;
+mrealmax:= 2.8;
+
+mimagstepsize := 0.01;
+numOfCores := 8;
+mimagstepsizeNoCtimes := numOfCores mimagstepsize;
+mimagmin := -0.5;
+mimagmax := 0.01;
+
+(*
+{t1,b1} =
+ AbsoluteTiming[
+ *)
+For[ mreal=mrealmin, mreal<=mrealmax, mreal+=mrealstepsize,
+ Print["mreal: ", mreal];
+ ofile=OpenAppend["./CMDudekDoubleViaDudekLUT.dat"];
+ For[ mimag = mimagmax, mimag >= mimagmin, mimag -= mimagstepsizeNoCtimes,{
+ SetSharedVariable[list];
+ list = {};
+ ParallelDo[{
+ currentCM = RNIntDudekDoubleViaDudek[(mreal + I mimag) (mreal + I mimag)][mDec11, mDec12, f2PoleDudek1, m2PoleDudek1, orbMomBW1, (mDec11+mDec12)^2, mDec21, mDec22, f2PoleDudek2, m2PoleDudek2, orbMomBW2,(mDec21+mDec22)^2 , scalingFacDudekSmearing];
+ AppendTo[list, {mimag, currentCM}]},
+ {mimag, mimag, mimag-mimagstepsizeNoCtimes+mimagstepsize, -mimagstepsize}
+ ],
+ list = SortBy[list, First, Greater];
+ WriteString[ofile, NumberForm[mreal,5], "\t", NumberForm[#[[1]], 5], "\t", NumberForm[Re[#[[2]]], 16], "\t", NumberForm[Im[#[[2]]], 16], "\n"] & /@ list;
+ };
+ ]
+ Close[ofile];
+ ];
+ (*
+];
+
+Print["t1: ", t1,"\tb1: ", b1];
+ *)
diff --git a/Scripts/Mathematica/CMs/LUTsSingle/CMDudekViaBBLUT.m b/Scripts/Mathematica/CMs/LUTsSingle/CMDudekViaBBLUT.m
new file mode 100644
index 0000000000000000000000000000000000000000..bd195056968aa6928ff3d09ffe6e56225eec7606
--- /dev/null
+++ b/Scripts/Mathematica/CMs/LUTsSingle/CMDudekViaBBLUT.m
@@ -0,0 +1,72 @@
+(* Copyright 2025
+ Bertram Kopf (bertram@ep1.rub.de)
+ Ruhr-Universit<C3><A4>t Bochum
+
+ This file is part of Pawian.
+
+ Pawian 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 3 of the License, or
+ (at your option) any later version.
+
+ Pawian 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 Pawian. If not, see <http://www.gnu.org/licenses/>
+
+***********************************************************************************************
+ *)
+
+Get[FileNameJoin[{Directory[], "../Integrals.m"}]]
+Get[FileNameJoin[{Directory[], "../InputParamsRho770Pi.m"}]]
+
+(* Precision[1.2]; *)
+(* Precision[12/10]; *)
+
+(*
+mrealstepsize:= 0.0005;
+mrealmin:= 0.005;
+mrealmax:= 3.5;
+
+mimagstepsize := 0.001;
+numOfCores := 1;
+mimagstepsizeNoCtimes := numOfCores mimagstepsize;
+mimagmin := -0.6;
+mimagmax := 0.005;
+ *)
+
+mrealstepsize:= 0.1;
+mrealmin:= 0.005;
+mrealmax:= 3.5;
+
+mimagstepsize := 0.1;
+numOfCores := 1;
+mimagstepsizeNoCtimes := numOfCores mimagstepsize;
+mimagmin := -0.6;
+mimagmax := 0.005;
+
+{t1,b1} =
+ AbsoluteTiming[
+For[ mreal=mrealmin, mreal<=mrealmax, mreal+=mrealstepsize,
+ Print["mreal: ", mreal];
+ ofile=OpenAppend["./CMDudekViaBBLUT.dat"];
+ For[ mimag = mimagmax, mimag >= mimagmin, mimag -= mimagstepsizeNoCtimes,{
+ SetSharedVariable[list];
+ list = {};
+ ParallelDo[{
+ currentCM = RNIntDudekviaBB[(mreal + I mimag) (mreal + I mimag)][mDec1, mDec2, mRec, f2PoleBB, m2PoleBB, (MBW MBW), orbMomBW, scalingFacBBSmearing];
+ AppendTo[list, {mimag, currentCM}]},
+ {mimag, mimag, mimag-mimagstepsizeNoCtimes+mimagstepsize, -mimagstepsize}
+ ],
+ list = SortBy[list, First, Greater];
+ WriteString[ofile, NumberForm[mreal,5], "\t", NumberForm[#[[1]], 5], "\t", NumberForm[Re[#[[2]]], 16], "\t", NumberForm[Im[#[[2]]], 16], "\n"] & /@ list;
+ };
+ ]
+ Close[ofile];
+ ];
+];
+
+Print["t1: ", t1,"\tb1: ", b1];
diff --git a/Scripts/Mathematica/CMs/LUTsSingle/CMDudekViaBWLUT.m b/Scripts/Mathematica/CMs/LUTsSingle/CMDudekViaBWLUT.m
new file mode 100644
index 0000000000000000000000000000000000000000..5e39192dca43d0ed47cf783bc0d89f79c22c2e1b
--- /dev/null
+++ b/Scripts/Mathematica/CMs/LUTsSingle/CMDudekViaBWLUT.m
@@ -0,0 +1,72 @@
+(* Copyright 2025
+ Bertram Kopf (bertram@ep1.rub.de)
+ Ruhr-Universit<C3><A4>t Bochum
+
+ This file is part of Pawian.
+
+ Pawian 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 3 of the License, or
+ (at your option) any later version.
+
+ Pawian 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 Pawian. If not, see <http://www.gnu.org/licenses/>
+
+***********************************************************************************************
+ *)
+
+Get[FileNameJoin[{Directory[], "../Integrals.m"}]]
+Get[FileNameJoin[{Directory[], "../InputParamsRho770Pi.m"}]]
+
+(* Precision[1.2]; *)
+(* Precision[12/10]; *)
+
+(*
+mrealstepsize:= 0.0005;
+mrealmin:= 0.005;
+mrealmax:= 3.5;
+
+mimagstepsize := 0.001;
+numOfCores := 1;
+mimagstepsizeNoCtimes := numOfCores mimagstepsize;
+mimagmin := -0.6;
+mimagmax := 0.005;
+ *)
+
+mrealstepsize:= 0.1;
+mrealmin:= 0.005;
+mrealmax:= 3.5;
+
+mimagstepsize := 0.1;
+numOfCores := 1;
+mimagstepsizeNoCtimes := numOfCores mimagstepsize;
+mimagmin := -0.6;
+mimagmax := 0.005;
+
+{t1,b1} =
+ AbsoluteTiming[
+For[ mreal=mrealmin, mreal<=mrealmax, mreal+=mrealstepsize,
+ Print["mreal: ", mreal];
+ ofile=OpenAppend["./CMDudekViaBWLUT.dat"];
+ For[ mimag = mimagmax, mimag >= mimagmin, mimag -= mimagstepsizeNoCtimes,{
+ SetSharedVariable[list];
+ list = {};
+ ParallelDo[{
+ currentCM = RNIntDudekviaBW[(mreal + I mimag) (mreal + I mimag)][mDec1, mDec2, mRec, orbMomBW, MBW, GammaBW, scalingFacBWSmearing];
+ AppendTo[list, {mimag, currentCM}]},
+ {mimag, mimag, mimag-mimagstepsizeNoCtimes+mimagstepsize, -mimagstepsize}
+ ],
+ list = SortBy[list, First, Greater];
+ WriteString[ofile, NumberForm[mreal,5], "\t", NumberForm[#[[1]], 5], "\t", NumberForm[Re[#[[2]]], 16], "\t", NumberForm[Im[#[[2]]], 16], "\n"] & /@ list;
+ };
+ ]
+ Close[ofile];
+ ];
+];
+
+Print["t1: ", t1,"\tb1: ", b1];
diff --git a/Scripts/Mathematica/CMs/LUTsSingle/CMDudekViaDudekLUT.m b/Scripts/Mathematica/CMs/LUTsSingle/CMDudekViaDudekLUT.m
new file mode 100644
index 0000000000000000000000000000000000000000..c61b3bf0e3b8ebf27168a6405190ab9ff438b900
--- /dev/null
+++ b/Scripts/Mathematica/CMs/LUTsSingle/CMDudekViaDudekLUT.m
@@ -0,0 +1,74 @@
+(* Copyright 2025
+ Bertram Kopf (bertram@ep1.rub.de)
+ Ruhr-Universit<C3><A4>t Bochum
+
+ This file is part of Pawian.
+
+ Pawian 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 3 of the License, or
+ (at your option) any later version.
+
+ Pawian 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 Pawian. If not, see <http://www.gnu.org/licenses/>
+
+***********************************************************************************************
+ *)
+
+Get[FileNameJoin[{Directory[], "../Integrals.m"}]]
+Get[FileNameJoin[{Directory[], "../InputParamsRho770Pi.m"}]]
+
+(* Precision[1.2]; *)
+(* Precision[12/10]; *)
+
+mrealstepsize:= 0.0005;
+mrealmin:= 0.005;
+mrealmax:= 3.5;
+
+mimagstepsize := 0.002;
+numOfCores := 1;
+mimagstepsizeNoCtimes := numOfCores mimagstepsize;
+mimagmin := -0.6;
+mimagmax := 0.002;
+
+(*
+mrealstepsize:= 0.1;
+mrealmin:= 0.005;
+mrealmax:= 3.5;
+
+mimagstepsize := 0.1;
+numOfCores := 1;
+mimagstepsizeNoCtimes := numOfCores mimagstepsize;
+mimagmin := -0.6;
+mimagmax := 0.005;
+
+{t1,b1} =
+ AbsoluteTiming[
+ *)
+For[ mreal=mrealmin, mreal<=mrealmax, mreal+=mrealstepsize,
+ Print["mreal: ", mreal];
+ ofile=OpenAppend["./CMDudekViaDudekLUT.dat"];
+ For[ mimag = mimagmax, mimag >= mimagmin, mimag -= mimagstepsizeNoCtimes,{
+ SetSharedVariable[list];
+ list = {};
+ ParallelDo[{
+ currentCM = RNIntDudekviaDudek[(mreal + I mimag) (mreal + I mimag)][mDec1, mDec2, mRec, f2PoleDudek, m2PoleDudek, (MBW MBW), orbMomBW, scalingFacDudekSmearing];
+ AppendTo[list, {mimag, currentCM}]},
+ {mimag, mimag, mimag-mimagstepsizeNoCtimes+mimagstepsize, -mimagstepsize}
+ ],
+ list = SortBy[list, First, Greater];
+ WriteString[ofile, NumberForm[mreal,5], "\t", NumberForm[#[[1]], 5], "\t", NumberForm[Re[#[[2]]], 16], "\t", NumberForm[Im[#[[2]]], 16], "\n"] & /@ list;
+ };
+ ]
+ Close[ofile];
+ ];
+ (*
+];
+
+Print["t1: ", t1,"\tb1: ", b1];
+ *)
diff --git a/Scripts/Mathematica/CMs/ParamExtraction/PoleExtractDouble_CM.m b/Scripts/Mathematica/CMs/ParamExtraction/PoleExtractDouble_CM.m
new file mode 100644
index 0000000000000000000000000000000000000000..80924d55705eb3e353498dbc54b64c8ac3a3f445
--- /dev/null
+++ b/Scripts/Mathematica/CMs/ParamExtraction/PoleExtractDouble_CM.m
@@ -0,0 +1,90 @@
+(* Copyright 2020
+ Meike Kuessner
+ mail: mkuessner@ep1.rub.de
+ Ruhr-Universität Bochum
+
+ modified by Bertram Kopf (bertram@ep1.rub.de)
+
+ This file is part of Pawian.
+
+ Pawian 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 3 of the License, or
+ (at your option) any later version.
+
+ Pawian 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 Pawian. If not, see <http://www.gnu.org/licenses/>
+
+***********************************************************************
+This Script computes the bare masses and bare coupling strengths of a resonance
+by its physical mass, width and masses of the decay products in case of a
+two-body decay.
+The calculation is based on K-matrix descriptions with Chew-Mandelstam functions
+ as defined in Phys.Rev.D 19 (1979) 239 ("BB") and in
+ Dudek et al Phys. Rev. D91, 054008 (2015), Appendix B ("Dudek")
+***********************************************************************
+*)
+
+Get[FileNameJoin[{Directory[], "../KinUtils.m"}]]
+(* Get[FileNameJoin[{Directory[], "../InputParamsDoubleSigmaSigma.m"}]] *)
+Get[FileNameJoin[{Directory[], "../InputParamsDoubleOmegaPhi.m"}]]
+
+Print["MBW1: ", MBW1,"\tGammaBW1: ", GammaBW1,"\tmDec11: ", mDec11,"\tmDec12: ", mDec12,"\torbMomBW1: ", orbMomBW1];
+Print["MBW2: ", MBW2,"\tGammaBW2: ", GammaBW2,"\tmDec21: ", mDec21,"\tmDec22: ", mDec22,"\torbMomBW2: ", orbMomBW2];
+
+
+ds2BB[s_,m1_,m2_,fR2_,mR2_, s0_, orbL_]:= Abs[mR2 - s - fR2 BarrierSimple[s, (m1+m2)^2, s0, orbL] Conjugate[cmReid[s,m1,m2]]]^2;
+
+
+
+(* Determine Minimum and save *)
+resultBB1=FindMinimum[SetPrecision[ds2BB[(MBW1 - I GammaBW1 /2.)^2,mDec11,mDec12,a,b, (MBW1 MBW1), orbMomBW1],30],{a,b},WorkingPrecision->20];
+{fsqBB1,msqBB1}={a,b}/.resultBB1[[2]];
+Print["m^2(BB1) = ", msqBB1, " ", "f^2(BB1) = ", fsqBB1];
+Print["m(BB1) = ", Sqrt[msqBB1], " ", "f(BB1) = ", Sqrt[fsqBB1]];
+
+(* Plot Contour around Minimum and Save Plot, change plot style here if needed or comment to not save the plot*)
+plotBB1=ContourPlot[Log[ds2BB[(MBW1-I 2 GammaBW1)^2,mDec11,mDec12,f,m, (MBW1 MBW1), orbMomBW1]],{f,fsqBB1-0.1,fsqBB1+0.1},{m,msqBB1-0.1,msqBB1+0.1},Contours->100,PlotLegends->Automatic,WorkingPrecision->30,ContourStyle->None, ColorFunction->"Rainbow", Epilog->({Black,PointSize[.03],Point[{fsqBB,msqBB}],Text["minimum BB 1",{fsqBB1,msqBB1*0.975}]}),FrameLabel->{f^2PoleBB1,m^2PoleBB1}];
+
+Export["PoleBB1.pdf",plotBB1];
+
+resultBB2=FindMinimum[SetPrecision[ds2BB[(MBW2 - I GammaBW2 /2.)^2,mDec21,mDec22,a,b, (MBW2 MBW2), orbMomBW2],30],{a,b},WorkingPrecision->20];
+{fsqBB2,msqBB2}={a,b}/.resultBB2[[2]];
+Print["m^2(BB2) = ", msqBB2, " ", "f^2(BB2) = ", fsqBB2];
+Print["m(BB2) = ", Sqrt[msqBB2], " ", "f(BB2) = ", Sqrt[fsqBB2]];
+
+(* Plot Contour around Minimum and Save Plot, change plot style here if needed or comment to not save the plot*)
+plotBB2=ContourPlot[Log[ds2BB[(MBW2-I 2 GammaBW2)^2,mDec21,mDec22,f,m, (MBW2 MBW2), orbMomBW2]],{f,fsqBB2-0.1,fsqBB2+0.1},{m,msqBB2-0.1,msqBB2+0.1},Contours->100,PlotLegends->Automatic,WorkingPrecision->30,ContourStyle->None, ColorFunction->"Rainbow", Epilog->({Black,PointSize[.03],Point[{fsqBB,msqBB}],Text["minimum BB 2",{fsqBB2,msqBB2*0.975}]}),FrameLabel->{f^2PoleBB2,m^2PoleBB2}];
+
+Export["PoleBB2.pdf",plotBB2];
+
+
+(* Extraction for Dudek CM *)
+ds2Dudek[s_,m1_,m2_,fR2_,mR2_, s0_, orbL_]:= Abs[mR2 - s + fR2 BarrierSimple[s, (m1+m2)^2, s0, orbL] Conjugate[cmDudek[s,m1,m2]]]^2;
+
+(* Determine Minimum and save *)
+resultDudek1=FindMinimum[SetPrecision[ds2Dudek[(MBW1 - I GammaBW1 /2.)^2,mDec11,mDec12,a,b, (MBW1 MBW1), orbMomBW1],30],{a,b},WorkingPrecision->20];
+{fsqDudek1,msqDudek1}={a,b}/.resultDudek1[[2]];
+Print["m^2(Dudek1)= ", msqDudek1, " ", "f^2(Dudek1)= ", fsqDudek1];
+Print["m(Dudek1)= ", Sqrt[msqDudek1], " ", "f(Dudek1)= ", Sqrt[fsqDudek1]];
+
+(* Plot Contour around Minimum and Save Plot, change plot style here if needed or comment to not save the plot*)
+plotDudek1=ContourPlot[Log[ds2[(MBW1-I 2 GammaBW1)^2,mDec11,mDec12,f,m, (MBW1 MBW1), orbMomBW1]],{f,fsqDudek1-0.1,fsqDudek1+0.1},{m,msqDudek1-0.1,msqDudek1+0.1},Contours->100,PlotLegends->Automatic,WorkingPrecision->30,ContourStyle->None, ColorFunction->"Rainbow", Epilog->({Black,PointSize[.03],Point[{fsqDudek1,msqDudek1}],Text["minimum Dudek 1",{fsqDudek1,msqDudek1*0.975}]}),FrameLabel->{f^2,m^2}];
+
+Export["PoleDudek1.pdf",plotDudek1];
+
+
+resultDudek2=FindMinimum[SetPrecision[ds2Dudek[(MBW2 - I GammaBW2 /2.)^2,mDec21,mDec22,a,b, (MBW2 MBW2), orbMomBW2],30],{a,b},WorkingPrecision->20];
+{fsqDudek2,msqDudek2}={a,b}/.resultDudek2[[2]];
+Print["m^2(Dudek2)= ", msqDudek2, " ", "f^2(Dudek2)= ", fsqDudek2];
+Print["m(Dudek2)= ", Sqrt[msqDudek2], " ", "f(Dudek2)= ", Sqrt[fsqDudek2]];
+
+(* Plot Contour around Minimum and Save Plot, change plot style here if needed or comment to not save the plot*)
+plotDudek2=ContourPlot[Log[ds2[(MBW2-I 2 GammaBW2)^2,mDec21,mDec22,f,m, (MBW2 MBW2), orbMomBW2]],{f,fsqDudek2-0.1,fsqDudek2+0.1},{m,msqDudek2-0.1,msqDudek2+0.1},Contours->100,PlotLegends->Automatic,WorkingPrecision->30,ContourStyle->None, ColorFunction->"Rainbow", Epilog->({Black,PointSize[.03],Point[{fsqDudek2,msqDudek2}],Text["minimum Dudek 2",{fsqDudek2,msqDudek2*0.975}]}),FrameLabel->{f^2,m^2}];
+
+Export["PoleDudek2.pdf",plotDudek2];
diff --git a/Scripts/Mathematica/CMs/ParamExtraction/PoleExtract_CM.m b/Scripts/Mathematica/CMs/ParamExtraction/PoleExtract_CM.m
new file mode 100644
index 0000000000000000000000000000000000000000..936ab68e524f013b7430f807fe2a7c9bedb64f21
--- /dev/null
+++ b/Scripts/Mathematica/CMs/ParamExtraction/PoleExtract_CM.m
@@ -0,0 +1,71 @@
+(* Copyright 2020
+ Meike Kuessner
+ mail: mkuessner@ep1.rub.de
+ Ruhr-Universität Bochum
+
+ modified by Bertram Kopf (bertram@ep1.rub.de)
+
+ This file is part of Pawian.
+
+ Pawian 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 3 of the License, or
+ (at your option) any later version.
+
+ Pawian 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 Pawian. If not, see <http://www.gnu.org/licenses/>
+
+***********************************************************************
+This Script computes the bare masses and bare coupling strengths of a resonance
+by its physical mass, width and masses of the decay products in case of a
+two-body decay.
+The calculation is based on K-matrix descriptions with Chew-Mandelstam functions
+ as defined in Phys.Rev.D 19 (1979) 239 ("BB") and in
+ Dudek et al Phys. Rev. D91, 054008 (2015), Appendix B ("Dudek")
+***********************************************************************
+*)
+
+Get[FileNameJoin[{Directory[], "../KinUtils.m"}]]
+(* Get[FileNameJoin[{Directory[], "../InputParamsPi1300Pi.m"}]] *)
+(* Get[FileNameJoin[{Directory[], "../InputParamsRho770Pi.m"}]] *)
+
+Get[FileNameJoin[{Directory[], "../InputParamsKstar892Pi.m"}]]
+
+Print["MBW: ", MBW,"\tGammaBW: ", GammaBW,"\tmDec1: ", mDec1,"\tmDec2: ", mDec2,"\torbL: ", orbMomBW];
+
+(* Singularity of denominator |d(s)|^2 of dispersion relation defines pole position, see Basdevant & Berger, Phys.Rev.D 19 (1979) 239 eq(12) and eq(14) *)
+
+ds2BB[s_,m1_,m2_,fR2_,mR2_, s0_, orbL_]:= Abs[mR2 - s - fR2 BarrierSimple[s, (m1+m2)^2, s0, orbL] Conjugate[cmReid[s,m1,m2]]]^2;
+
+
+
+(* Determine Minimum and save *)
+resultBB=FindMinimum[SetPrecision[ds2BB[(MBW - I GammaBW /2.)^2,mDec1,mDec2,a,b, (MBW MBW), orbMomBW],30],{a,b},WorkingPrecision->20];
+{fsqBB,msqBB}={a,b}/.resultBB[[2]];
+Print["m^2(BB) = ", msqBB, " ", "f^2(BB) = ", fsqBB];
+Print["m(BB) = ", Sqrt[msqBB], " ", "f(BB) = ", Sqrt[fsqBB]];
+
+(* Plot Contour around Minimum and Save Plot, change plot style here if needed or comment to not save the plot*)
+plotBB=ContourPlot[Log[ds2BB[(MBW-I 2 GammaBW)^2,mDec1,mDec2,f,m, (MBW MBW), orbMomBW]],{f,fsqBB-0.1,fsqBB+0.1},{m,msqBB-0.1,msqBB+0.1},Contours->100,PlotLegends->Automatic,WorkingPrecision->30,ContourStyle->None, ColorFunction->"Rainbow", Epilog->({Black,PointSize[.03],Point[{fsqBB,msqBB}],Text["minimum BB",{fsqBB,msqBB*0.975}]}),FrameLabel->{f^2PoleBB,m^2PoleBB}];
+
+Export["PoleBB.pdf",plotBB];
+
+
+(* Extraction for Dudek CM *)
+ds2Dudek[s_,m1_,m2_,fR2_,mR2_, s0_, orbL_]:= Abs[mR2 - s + fR2 BarrierSimple[s, (m1+m2)^2, s0, orbL] Conjugate[cmDudek[s,m1,m2]]]^2;
+
+(* Determine Minimum and save *)
+resultDudek=FindMinimum[SetPrecision[ds2Dudek[(MBW - I GammaBW /2.)^2,mDec1,mDec2,a,b, (MBW MBW), orbMomBW],30],{a,b},WorkingPrecision->20];
+{fsqDudek,msqDudek}={a,b}/.resultDudek[[2]];
+Print["m^2(Dudek)= ", msqDudek, " ", "f^2(Dudek)= ", fsqDudek];
+Print["m(Dudek)= ", Sqrt[msqDudek], " ", "f(Dudek)= ", Sqrt[fsqDudek]];
+
+(* Plot Contour around Minimum and Save Plot, change plot style here if needed or comment to not save the plot*)
+plotDudek=ContourPlot[Log[ds2[(MBW-I 2 GammaBW)^2,mDec1,mDec2,f,m, (MBW MBW), orbMomBW]],{f,fsqDudek-0.1,fsqDudek+0.1},{m,msqDudek-0.1,msqDudek+0.1},Contours->100,PlotLegends->Automatic,WorkingPrecision->30,ContourStyle->None, ColorFunction->"Rainbow", Epilog->({Black,PointSize[.03],Point[{fsqDudek,msqDudek}],Text["minimum Dudek",{fsqDudek,msqDudek*0.975}]}),FrameLabel->{f^2,m^2}];
+
+Export["PoleDudek.pdf",plotDudek];
diff --git a/Scripts/Mathematica/CMs/SmearDouble/CMDudekDoubleUnstable.m b/Scripts/Mathematica/CMs/SmearDouble/CMDudekDoubleUnstable.m
new file mode 100644
index 0000000000000000000000000000000000000000..c1428ddf269508c4f5843e088a5f74f3b750eb0e
--- /dev/null
+++ b/Scripts/Mathematica/CMs/SmearDouble/CMDudekDoubleUnstable.m
@@ -0,0 +1,34 @@
+(* Copyright 2025
+ Bertram Kopf (bertram@ep1.rub.de)
+ Ruhr-Universit<C3><A4>t Bochum
+
+ This file is part of Pawian.
+
+ Pawian 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 3 of the License, or
+ (at your option) any later version.
+
+ Pawian 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 Pawian. If not, see <http://www.gnu.org/licenses/>
+
+***********************************************************************************************
+ *)
+Get[FileNameJoin[{Directory[], "../KinUtils.m"}]]
+Get[FileNameJoin[{Directory[], "../Integrals.m"}]]
+Get[FileNameJoin[{Directory[], "../InputParamsDoubleSigmaSigma.m"}]]
+
+(* Precision[1.2]; *)
+(* Precision[12/10]; *)
+
+Print["RNIntDudekDoubleViaDudek[(1.0 - I 0.1) (1.0 - I 0.1)][mDec11, mDec12, f2PoleDudek1, m2PoleDudek1, orbMomBW1, (mDec11+mDec12)^2, mDec21, mDec22, f2PoleDudek2, m2PoleDudek2, orbMomBW2,(mDec21+mDec22)^2 , scalingFacDudekSmearing]: " , RNIntDudekDoubleViaDudek[(1.0 - I 0.1) (1.0 - I 0.1)][mDec11, mDec12, f2PoleDudek1, m2PoleDudek1, orbMomBW1, (mDec11+mDec12)^2, mDec21, mDec22, f2PoleDudek2, m2PoleDudek2, orbMomBW2,(mDec21+mDec22)^2 , scalingFacDudekSmearing]];
+
+plotCMDudekDoubleviaDudekUnstable = ReImPlot[RNIntDudekDoubleViaDudek[(b + I valImag) (b + I valImag)][mDec11, mDec12, f2PoleDudek1, m2PoleDudek1, orbMomBW1, (mDec11+mDec12)^2, mDec21, mDec22, f2PoleDudek2, m2PoleDudek2, orbMomBW2,(mDec21+mDec22)^2 , scalingFacDudekSmearing], {b, 0., 2.5}, PlotStyle -> {Purple,{Purple,Dashed}}, PlotLegends->Automatic, AxesLabel->{Sqrt[s] [GeV/(c c)], CM Dudek double via Dudek unstable}];
+
+Export["CMDudekDoubleViaDudek.pdf", plotCMDudekDoubleviaDudekUnstable];
+
diff --git a/Scripts/Mathematica/CMs/SmearDouble/PlotSmearingDoubleBB.m b/Scripts/Mathematica/CMs/SmearDouble/PlotSmearingDoubleBB.m
new file mode 100644
index 0000000000000000000000000000000000000000..9b373e43efaebf9d16ae821f195d43a1b6796eda
--- /dev/null
+++ b/Scripts/Mathematica/CMs/SmearDouble/PlotSmearingDoubleBB.m
@@ -0,0 +1,41 @@
+(* Copyright 2025
+ Bertram Kopf (bertram@ep1.rub.de)
+ Ruhr-Universit<C3><A4>t Bochum
+
+ This file is part of Pawian.
+
+ Pawian 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 3 of the License, or
+ (at your option) any later version.
+
+ Pawian 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 Pawian. If not, see <http://www.gnu.org/licenses/>
+
+***********************************************************************************************
+ *)
+Get[FileNameJoin[{Directory[], "../KinUtils.m"}]]
+Get[FileNameJoin[{Directory[], "../Integrals.m"}]]
+(* Get[FileNameJoin[{Directory[], "../InputParamsDoubleSigmaSigma.m"}]] *)
+Get[FileNameJoin[{Directory[], "../InputParamsDoubleOmegaPhi.m"}]]
+(* Precision[1.2]; *)
+(* Precision[12/10]; *)
+
+plotIntegrFuncBB3Mass1 = Plot[ {f2PoleBB1 ImagSigmads2BB[b b, mDec11, mDec12, f2PoleBB1, m2PoleBB1, (MBW1 MBW1), orbMomBW1]},{b, (mDec11+mDec12), 3.}, AxesLabel->{M [GeV/(c c)], integr} ];
+Export["IntegrFuncBB3Mass1.pdf", plotIntegrFuncBB3Mass1];
+
+plotIntegrFuncBB3Mass2 = Plot[ {f2PoleBB2 ImagSigmads2BB[b b, mDec21, mDec22, f2PoleBB2, m2PoleBB2, (MBW2 MBW2), orbMomBW2]},{b, (mDec21+mDec22), 3.}, AxesLabel->{M [GeV/(c c)], integr} ];
+Export["IntegrFuncBB3Mass2.pdf", plotIntegrFuncBB3Mass2];
+
+Print["RNIntBBDouble[mDec11, mDec12, f2PoleBB1, m2PoleBB1, orbMomBW1, (mDec11+mDec12)^2, mDec21, mDec22, f2PoleBB2, m2PoleBB2, orbMomBW2,(mDec21+mDec22)^2 , scalingFacBBSmearing, 10.]: ", RNIntBBDouble[mDec11, mDec12, f2PoleBB1, m2PoleBB1, orbMomBW1, (mDec11+mDec12)^2, mDec21, mDec22, f2PoleBB2, m2PoleBB2, orbMomBW2,(mDec21+mDec22)^2 , scalingFacBBSmearing, 10.]];
+
+Print["RNIntBBDouble[mDec11, mDec12, f2PoleBB1, m2PoleBB1, orbMomBW1, (mDec11+mDec12)^2, mDec21, mDec22, f2PoleBB2, m2PoleBB2, orbMomBW2,(mDec21+mDec22)^2 , scalingFacBBSmearing, 1000.]: ", RNIntBBDouble[mDec11, mDec12, f2PoleBB1, m2PoleBB1, orbMomBW1, (mDec11+mDec12)^2, mDec21, mDec22, f2PoleBB2, m2PoleBB2, orbMomBW2,(mDec21+mDec22)^2 , scalingFacBBSmearing, 1000.]];
+
+Print["RNIntBBDouble[mDec11, mDec12, f2PoleBB1, m2PoleBB1, orbMomBW1, (mDec11+mDec12)^2, mDec21, mDec22, f2PoleBB2, m2PoleBB2, orbMomBW2,(mDec21+mDec22)^2 , scalingFacBBSmearing, 100000.]: ", RNIntBBDouble[mDec11, mDec12, f2PoleBB1, m2PoleBB1, orbMomBW1, (mDec11+mDec12)^2, mDec21, mDec22, f2PoleBB2, m2PoleBB2, orbMomBW2,(mDec21+mDec22)^2 , scalingFacBBSmearing, 100000.]];
+
+
diff --git a/Scripts/Mathematica/CMs/SmearDouble/PlotSmearingDoubleBW.m b/Scripts/Mathematica/CMs/SmearDouble/PlotSmearingDoubleBW.m
new file mode 100644
index 0000000000000000000000000000000000000000..acc50f03875fa7207b50a45de7dec212487271c8
--- /dev/null
+++ b/Scripts/Mathematica/CMs/SmearDouble/PlotSmearingDoubleBW.m
@@ -0,0 +1,45 @@
+(* Copyright 2025
+ Bertram Kopf (bertram@ep1.rub.de)
+ Ruhr-Universit<C3><A4>t Bochum
+
+ This file is part of Pawian.
+
+ Pawian 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 3 of the License, or
+ (at your option) any later version.
+
+ Pawian 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 Pawian. If not, see <http://www.gnu.org/licenses/>
+
+***********************************************************************************************
+ *)
+Get[FileNameJoin[{Directory[], "../KinUtils.m"}]]
+Get[FileNameJoin[{Directory[], "../Integrals.m"}]]
+(* Get[FileNameJoin[{Directory[], "../InputParamsDoubleSigmaSigma.m"}]] *)
+Get[FileNameJoin[{Directory[], "../InputParamsDoubleOmegaPhi.m"}]]
+
+(* Precision[1.2]; *)
+(* Precision[12/10]; *)
+
+plotIntegrFuncBW3Mass1 = Plot[ {rho[b b, mDec11, mDec12] Abs[BreitWignerBlattWRelNorm[b b, orbMomBW1, MBW1, GammaBW1, mDec11, mDec12]]^2},{b, (mDec11+mDec12), 3.}, AxesLabel->{M [GeV/(c c)], integr} ];
+Export["IntegrFuncBW3Mass1.pdf", plotIntegrFuncBW3Mass1];
+
+plotIntegrFuncBW3Mass2 = Plot[ {rho[b b, mDec21, mDec22] Abs[BreitWignerBlattWRelNorm[b b, orbMomBW2, MBW2, GammaBW2, mDec21, mDec22]]^2},{b, (mDec21+mDec22), 3.}, AxesLabel->{M [GeV/(c c)], integr} ];
+Export["IntegrFuncBW3Mass2.pdf", plotIntegrFuncBW3Mass2];
+
+
+Print["RNIntBWDouble[MBW1, GammaBW1, mDec11, mDec12, orbMomBW1, MBW2, GammaBW2, mDec21, mDec22, orbMomBW2, scalingFacBWSmearing, 10.]: ", RNIntBWDouble[MBW1, GammaBW1, mDec11, mDec12, orbMomBW1, MBW2, GammaBW2, mDec21, mDec22, orbMomBW2, scalingFacBWSmearing, 10.]];
+
+
+Print["RNIntBWDouble[MBW1, GammaBW1, mDec11, mDec12, orbMomBW1, MBW2, GammaBW2, mDec21, mDec22, orbMomBW2, scalingFacBWSmearing, 1000.]: ", RNIntBWDouble[MBW1, GammaBW1, mDec11, mDec12, orbMomBW1, MBW2, GammaBW2, mDec21, mDec22, orbMomBW2, scalingFacBWSmearing, 1000.]];
+
+Print["RNIntBWDouble[MBW1, GammaBW1, mDec11, mDec12, orbMomBW1, MBW2, GammaBW2, mDec21, mDec22, orbMomBW2, scalingFacBWSmearing, 100000.]: ", RNIntBWDouble[MBW1, GammaBW1, mDec11, mDec12, orbMomBW1, MBW2, GammaBW2, mDec21, mDec22, orbMomBW2, scalingFacBWSmearing, 100000.]];
+
+
+
diff --git a/Scripts/Mathematica/CMs/SmearDouble/PlotSmearingDoubleDudek.m b/Scripts/Mathematica/CMs/SmearDouble/PlotSmearingDoubleDudek.m
new file mode 100644
index 0000000000000000000000000000000000000000..7a7df434f0c59736a443bdb1d593dbb939c7487e
--- /dev/null
+++ b/Scripts/Mathematica/CMs/SmearDouble/PlotSmearingDoubleDudek.m
@@ -0,0 +1,42 @@
+(* Copyright 2025
+ Bertram Kopf (bertram@ep1.rub.de)
+ Ruhr-Universit<C3><A4>t Bochum
+
+ This file is part of Pawian.
+
+ Pawian 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 3 of the License, or
+ (at your option) any later version.
+
+ Pawian 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 Pawian. If not, see <http://www.gnu.org/licenses/>
+
+***********************************************************************************************
+ *)
+Get[FileNameJoin[{Directory[], "../KinUtils.m"}]]
+Get[FileNameJoin[{Directory[], "../Integrals.m"}]]
+(* Get[FileNameJoin[{Directory[], "../InputParamsDoubleSigmaSigma.m"}]] *)
+Get[FileNameJoin[{Directory[], "../InputParamsDoubleOmegaPhi.m"}]]
+
+(* Precision[1.2]; *)
+(* Precision[12/10]; *)
+
+plotIntegrFuncDudek3Mass1 = Plot[ {f2PoleDudek1 ImagSigmads2Dudek[b b, mDec11, mDec12, f2PoleDudek1, m2PoleDudek1, (MBW1 MBW1), orbMomBW1]},{b, (mDec11+mDec12), 3.}, AxesLabel->{M [GeV/(c c)], integr} ];
+Export["IntegrFuncDudek3Mass1.pdf", plotIntegrFuncDudek3Mass1];
+
+plotIntegrFuncDudek3Mass2 = Plot[ {f2PoleDudek2 ImagSigmads2Dudek[b b, mDec21, mDec22, f2PoleDudek2, m2PoleDudek2, (MBW2 MBW2), orbMomBW2]},{b, (mDec21+mDec22), 3.}, AxesLabel->{M [GeV/(c c)], integr} ];
+Export["IntegrFuncDudek3Mass2.pdf", plotIntegrFuncDudek3Mass2];
+
+Print["RNIntDudekDouble[mDec11, mDec12, f2PoleDudek1, m2PoleDudek1, orbMomBW1, (mDec11+mDec12)^2, mDec21, mDec22, f2PoleDudek2, m2PoleDudek2, orbMomBW2,(mDec21+mDec22)^2 , scalingFacDudekSmearing, 10.]: ", RNIntDudekDouble[mDec11, mDec12, f2PoleDudek1, m2PoleDudek1, orbMomBW1, (mDec11+mDec12)^2, mDec21, mDec22, f2PoleDudek2, m2PoleDudek2, orbMomBW2,(mDec21+mDec22)^2 , scalingFacDudekSmearing, 10.]];
+
+Print["RNIntDudekDouble[mDec11, mDec12, f2PoleDudek1, m2PoleDudek1, orbMomBW1, (mDec11+mDec12)^2, mDec21, mDec22, f2PoleDudek2, m2PoleDudek2, orbMomBW2,(mDec21+mDec22)^2 , scalingFacDudekSmearing, 1000.]: ", RNIntDudekDouble[mDec11, mDec12, f2PoleDudek1, m2PoleDudek1, orbMomBW1, (mDec11+mDec12)^2, mDec21, mDec22, f2PoleDudek2, m2PoleDudek2, orbMomBW2,(mDec21+mDec22)^2 , scalingFacDudekSmearing, 1000.]];
+
+Print["RNIntDudekDouble[mDec11, mDec12, f2PoleDudek1, m2PoleDudek1, orbMomBW1, (mDec11+mDec12)^2, mDec21, mDec22, f2PoleDudek2, m2PoleDudek2, orbMomBW2,(mDec21+mDec22)^2 , scalingFacDudekSmearing, 100000.]: ", RNIntDudekDouble[mDec11, mDec12, f2PoleDudek1, m2PoleDudek1, orbMomBW1, (mDec11+mDec12)^2, mDec21, mDec22, f2PoleDudek2, m2PoleDudek2, orbMomBW2,(mDec21+mDec22)^2 , scalingFacDudekSmearing, 100000.]];
+
+
diff --git a/Scripts/Mathematica/CMs/SmearSingle/CMBBUnstable.m b/Scripts/Mathematica/CMs/SmearSingle/CMBBUnstable.m
new file mode 100644
index 0000000000000000000000000000000000000000..8dbca1cdb590cec953c9d839a19c42de01d3026b
--- /dev/null
+++ b/Scripts/Mathematica/CMs/SmearSingle/CMBBUnstable.m
@@ -0,0 +1,56 @@
+(* Copyright 2025
+ Bertram Kopf (bertram@ep1.rub.de)
+ Ruhr-Universit<C3><A4>t Bochum
+
+ This file is part of Pawian.
+
+ Pawian 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 3 of the License, or
+ (at your option) any later version.
+
+ Pawian 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 Pawian. If not, see <http://www.gnu.org/licenses/>
+
+***********************************************************************************************
+ *)
+Get[FileNameJoin[{Directory[], "../KinUtils.m"}]]
+Get[FileNameJoin[{Directory[], "../Integrals.m"}]]
+(* Get[FileNameJoin[{Directory[], "../InputParamsPi1300Pi.m"}]] *)
+Get[FileNameJoin[{Directory[], "../InputParamsRho770Pi.m"}]]
+
+(* Precision[1.2]; *)
+(* Precision[12/10]; *)
+
+plotCMBBStable = ReImPlot[cmReid[(b + I valImag) (b + I valImag), mRec, MBW], {b, 0., 4.0}, PlotStyle -> {Red,{Red,Dashed}}, PlotLegends->Automatic, AxesLabel->{Sqrt[s] [GeV/(c c)], CM BB stable}];
+Export["CMBBStable.pdf", plotCMBBStable];
+
+plotCMBBviaBBUnstable = ReImPlot[RNIntBBviaBB[(b + I valImag) (b + I valImag)][mDec1, mDec2, mRec, f2PoleBB, m2PoleBB, (MBW MBW), orbMomBW, scalingFacBBSmearing], {b, 0., 4.0}, PlotStyle -> {Blue,{Blue,Dashed}}, PlotLegends->Automatic, AxesLabel->{Sqrt[s] [GeV/(c c)], CM BB via BB unstable}];
+Export["CMBBViaBBUnstable.pdf", plotCMBBviaBBUnstable];
+
+plCMBBviaBBComp=Show[plotCMBBStable, plotCMBBviaBBUnstable];
+Export["CMBBViaBB_compare.pdf", plCMBBviaBBComp];
+
+ plotCMBBviaDudekUnstable = ReImPlot[RNIntBBviaDudek[(b + I valImag) (b + I valImag)][mDec1, mDec2, mRec, f2PoleDudek, m2PoleDudek, (MBW MBW), orbMomBW, scalingFacDudekSmearing], {b, 0., 4.0}, PlotStyle -> {Purple,{Purple,Dashed}}, PlotLegends->Automatic, AxesLabel->{Sqrt[s] [GeV/(c c)], CM BB via BB unstable}];
+Export["CMBBViaDudekUnstable.pdf", plotCMBBviaDudekUnstable];
+
+plCMBBviaDudekComp=Show[plotCMBBStable, plotCMBBviaDudekUnstable];
+Export["CMBBViaDudek_compare.pdf", plCMBBviaDudekComp];
+
+plCMBBviaDudekvsBBComp=Show[plotCMBBviaBBUnstable, plotCMBBviaDudekUnstable];
+Export["CMBBviaDudekvsBB_compare.pdf", plCMBBviaDudekvsBBComp];
+
+
+plotCMBBviaBWUnstable = ReImPlot[RNIntBBviaBW[(b + I valImag) (b + I valImag)][mDec1, mDec2, mRec, orbMomBW, MBW, GammaBW, scalingFacBWSmearing], {b, 0., 4.0}, PlotStyle -> {Purple,{Purple,Dashed}}, PlotLegends->Automatic, AxesLabel->{Sqrt[s] [GeV/(c c)], CM BB via BB unstable}];
+Export["CMBBViaBWUnstable.pdf", plotCMBBviaBWUnstable];
+
+plCMBBviaBWComp=Show[plotCMBBStable, plotCMBBviaBWUnstable];
+Export["CMBBViaBW_compare.pdf", plCMBBviaBWComp];
+
+plCMBBviaDudekvsBWComp=Show[plotCMBBviaBWUnstable, plotCMBBviaDudekUnstable];
+Export["CMBBviaDudekvsBW_compare.pdf", plCMBBviaDudekvsBWComp];
diff --git a/Scripts/Mathematica/CMs/SmearSingle/CMDudekUnstable.m b/Scripts/Mathematica/CMs/SmearSingle/CMDudekUnstable.m
new file mode 100644
index 0000000000000000000000000000000000000000..49581635e54b4ad8df2198f9e4104698f911e21f
--- /dev/null
+++ b/Scripts/Mathematica/CMs/SmearSingle/CMDudekUnstable.m
@@ -0,0 +1,57 @@
+(* Copyright 2025
+ Bertram Kopf (bertram@ep1.rub.de)
+ Ruhr-Universit<C3><A4>t Bochum
+
+ This file is part of Pawian.
+
+ Pawian 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 3 of the License, or
+ (at your option) any later version.
+
+ Pawian 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 Pawian. If not, see <http://www.gnu.org/licenses/>
+
+***********************************************************************************************
+ *)
+Get[FileNameJoin[{Directory[], "../KinUtils.m"}]]
+Get[FileNameJoin[{Directory[], "../Integrals.m"}]]
+(* Get[FileNameJoin[{Directory[], "../InputParamsPi1300Pi.m"}]] *)
+Get[FileNameJoin[{Directory[], "../InputParamsRho770Pi.m"}]]
+
+(* Precision[1.2]; *)
+(* Precision[12/10]; *)
+
+plotCMDudekStable = ReImPlot[cmDudek[(b + I valImag) (b + I valImag), mRec, MBW], {b, 0., 4.0}, PlotStyle -> {Red,{Red,Dashed}}, PlotLegends->Automatic, AxesLabel->{Sqrt[s] [GeV/(c c)], CM Dudek stable}];
+Export["CMDudekStable.pdf", plotCMDudekStable];
+
+plotCMDudekviaBBUnstable = ReImPlot[RNIntDudekviaBB[(b + I valImag) (b + I valImag)][mDec1, mDec2, mRec, f2PoleBB, m2PoleBB, (MBW MBW), orbMomBW, scalingFacBBSmearing], {b, 0., 4.0}, PlotStyle -> {Blue,{Blue,Dashed}}, PlotLegends->Automatic, AxesLabel->{Sqrt[s] [GeV/(c c)], CM Dudek via BB unstable}];
+Export["CMDudekViaBBUnstable.pdf", plotCMDudekviaBBUnstable];
+
+plCMDudekviaBBComp=Show[plotCMDudekStable, plotCMDudekviaBBUnstable];
+Export["CMDudekViaBB_compare.pdf", plCMDudekviaBBComp];
+
+ plotCMDudekviaDudekUnstable = ReImPlot[RNIntDudekviaDudek[(b + I valImag) (b + I valImag)][mDec1, mDec2, mRec, f2PoleDudek, m2PoleDudek, (MBW MBW), orbMomBW, scalingFacDudekSmearing], {b, 0., 4.0}, PlotStyle -> {Purple,{Purple,Dashed}}, PlotLegends->Automatic, AxesLabel->{Sqrt[s] [GeV/(c c)], CM Dudek via Dudek unstable}];
+Export["CMDudekViaDudekUnstable.pdf", plotCMDudekviaDudekUnstable];
+
+plCMDudekviaDudekComp=Show[plotCMDudekStable, plotCMDudekviaDudekUnstable];
+Export["CMDudekViaDudek_compare.pdf", plCMDudekviaDudekComp];
+
+
+plCMDudekviaDudekvsBBComp=Show[plotCMDudekviaBBUnstable, plotCMDudekviaDudekUnstable];
+Export["CMDudekviaDudekvsBB_compare.pdf", plCMDudekviaDudekvsBBComp];
+
+
+plotCMDudekviaBWUnstable = ReImPlot[RNIntDudekviaBW[(b + I valImag) (b + I valImag)][mDec1, mDec2, mRec, orbMomBW, MBW, GammaBW, scalingFacBWSmearing], {b, 0., 4.0}, PlotStyle -> {Purple,{Purple,Dashed}}, PlotLegends->Automatic, AxesLabel->{Sqrt[s] [GeV/(c c)], CM Dudek via BW unstable}];
+Export["CMDudekViaBWUnstable.pdf", plotCMDudekviaBWUnstable];
+
+plCMDudekviaBWComp=Show[plotCMDudekStable, plotCMDudekviaBWUnstable];
+Export["CMDudekViaBW_compare.pdf", plCMDudekviaBWComp];
+
+plCMDudekviaDudekvsBWComp=Show[plotCMDudekviaBWUnstable, plotCMDudekviaDudekUnstable];
+Export["CMDudekviaDudekvsBW_compare.pdf", plCMDudekviaDudekvsBWComp];
diff --git a/Scripts/Mathematica/CMs/SmearSingle/PlotSmearingBB.m b/Scripts/Mathematica/CMs/SmearSingle/PlotSmearingBB.m
new file mode 100644
index 0000000000000000000000000000000000000000..4d44f9b1c1b63f91cb9d6c384ce37e1cf4c29f0a
--- /dev/null
+++ b/Scripts/Mathematica/CMs/SmearSingle/PlotSmearingBB.m
@@ -0,0 +1,62 @@
+(* Copyright 2025
+ Bertram Kopf (bertram@ep1.rub.de)
+ Ruhr-Universit<C3><A4>t Bochum
+
+ This file is part of Pawian.
+
+ Pawian 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 3 of the License, or
+ (at your option) any later version.
+
+ Pawian 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 Pawian. If not, see <http://www.gnu.org/licenses/>
+
+***********************************************************************************************
+ *)
+Get[FileNameJoin[{Directory[], "../KinUtils.m"}]]
+Get[FileNameJoin[{Directory[], "../Integrals.m"}]]
+(* Get[FileNameJoin[{Directory[], "../InputParamsPi1300Pi.m"}]] *)
+Get[FileNameJoin[{Directory[], "../InputParamsRho770Pi.m"}]]
+(* Precision[1.2]; *)
+(* Precision[12/10]; *)
+
+
+plotIntegrFuncBB3Mass = Plot[ {f2PoleBB ImagSigmads2BB[b b, mDec1, mDec2, f2PoleBB, m2PoleBB, (MBW MBW), orbMomBW]},{b, (mDec1+mDec2), 3.}, AxesLabel->{M [GeV/(c c)], integr} ];
+Export["IntegrFuncBB3Mass.pdf", plotIntegrFuncBB3Mass];
+
+plotIntegrFuncBB3 = Plot[ {f2PoleBB ImagSigmads2BB[b, mDec1, mDec2, f2PoleBB, m2PoleBB, (MBW MBW), orbMomBW]},{b, (mDec1+mDec2)^2, 3.}, AxesLabel->{sprime [GeV/(c c)], integr} ];
+Export["IntegrFuncBB3.pdf", plotIntegrFuncBB3];
+
+plotIntegrFuncBB10 = Plot[ {f2PoleBB ImagSigmads2BB[b, mDec1, mDec2, f2PoleBB, m2PoleBB, (MBW MBW), orbMomBW]},{b, (mDec1+mDec2)^2, 10.}, AxesLabel->{sprime [GeV/(c c c c)], integr} ];
+Export["IntegrFuncBB10.pdf", plotIntegrFuncBB10];
+
+plotIntegrFuncBB100 = Plot[ {f2PoleBB ImagSigmads2BB[b, mDec1, mDec2, f2PoleBB, m2PoleBB, (MBW MBW), orbMomBW]},{b, (mDec1+mDec2)^2, 100.}, AxesLabel->{sprime [GeV GeV/(c c c c)], integr} ];
+Export["IntegrFuncBB100.pdf", plotIntegrFuncBB100];
+
+plotIntegrFuncBB1000 = Plot[ {f2PoleBB ImagSigmads2BB[b, mDec1, mDec2, f2PoleBB, m2PoleBB, (MBW MBW), orbMomBW]},{b, (mDec1+mDec2)^2, 1000.}, AxesLabel->{sprime [GeV GeV/(c c c c)], integr} ];
+Export["IntegrFuncBB1000.pdf", plotIntegrFuncBB1000];
+
+plotIntegrFuncBB10000 = Plot[ {f2PoleBB ImagSigmads2BB[b, mDec1, mDec2, f2PoleBB, m2PoleBB, (MBW MBW), orbMomBW]},{b, (mDec1+mDec2)^2, 10000.}, AxesLabel->{sprime [GeV GeV/(c c c c)], integr} ];
+Export["IntegrFuncBB10000.pdf", plotIntegrFuncBB10000];
+
+Print["intergal s<3: ", RNIntBB[mDec1, mDec2, f2PoleBB, m2PoleBB, (MBW MBW), orbMomBW, scalingFacBBSmearing, 3.]];
+Print["intergal s<100000: ", RNIntBB[mDec1, mDec2, f2PoleBB, m2PoleBB, (MBW MBW), orbMomBW, scalingFacBBSmearing, 100000.]];
+
+plotIntegralBBsmax10 = Plot[{RNIntBB[mDec1, mDec2, f2PoleBB, m2PoleBB, (MBW MBW), orbMomBW, scalingFacBBSmearing, b]}, {b, 3., 10.}, AxesLabel->{smax [GeV GeV/(c c c c)], integral} ];
+Export["IntegralBBsmax10.pdf", plotIntegralBBsmax10];
+
+plotIntegralBBsmax100 = Plot[{RNIntBB[mDec1, mDec2, f2PoleBB, m2PoleBB, (MBW MBW), orbMomBW, scalingFacBBSmearing, b]}, {b, 3., 100.}, AxesLabel->{smax [GeV GeV/(c c c c)], integral} ];
+Export["IntegralBBsmax100.pdf", plotIntegralBBsmax100];
+
+plotIntegralBBsmax1000 = Plot[{RNIntBB[mDec1, mDec2, f2PoleBB, m2PoleBB, (MBW MBW), orbMomBW, scalingFacBBSmearing, b]}, {b, 3., 1000.}, AxesLabel->{smax [GeV GeV/(c c c c)], integral} ];
+Export["IntegralBBsmax1000.pdf", plotIntegralBBsmax1000];
+
+plotIntegralBBsmax10000 = Plot[{RNIntBB[mDec1, mDec2, f2PoleBB, m2PoleBB, (MBW MBW), orbMomBW, scalingFacBBSmearing, b]}, {b, 3., 10000.}, AxesLabel->{smax [GeV GeV/(c c c c)], integral} ];
+Export["IntegralBBsmax10000.pdf", plotIntegralBBsmax10000];
+
diff --git a/Scripts/Mathematica/CMs/SmearSingle/PlotSmearingBW.m b/Scripts/Mathematica/CMs/SmearSingle/PlotSmearingBW.m
new file mode 100644
index 0000000000000000000000000000000000000000..e33cc12d6e753308471742d6c32a8a0c3326fc18
--- /dev/null
+++ b/Scripts/Mathematica/CMs/SmearSingle/PlotSmearingBW.m
@@ -0,0 +1,50 @@
+(* Copyright 2025
+ Bertram Kopf (bertram@ep1.rub.de)
+ Ruhr-Universit<C3><A4>t Bochum
+
+ This file is part of Pawian.
+
+ Pawian 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 3 of the License, or
+ (at your option) any later version.
+
+ Pawian 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 Pawian. If not, see <http://www.gnu.org/licenses/>
+
+***********************************************************************************************
+ *)
+Get[FileNameJoin[{Directory[], "../KinUtils.m"}]]
+Get[FileNameJoin[{Directory[], "../Integrals.m"}]]
+(* Get[FileNameJoin[{Directory[], "../InputParamsPi1300Pi.m"}]] *)
+Get[FileNameJoin[{Directory[], "../InputParamsRho770Pi.m"}]]
+
+(* Precision[1.2]; *)
+(* Precision[12/10]; *)
+
+plotIntegrFuncBW3Mass = Plot[ {scalingFacBWSmearing rho[b b, mDec1, mDec2] Abs[BreitWignerBlattWRelNorm[b b, orbMomBW, MBW, GammaBW, mDec1, mDec2]]^2},{b, (mDec1+mDec2), 3.}, AxesLabel->{M [GeV/(c c)], integr} ];
+Export["IntegrFuncBW3Mass.pdf", plotIntegrFuncBW3Mass];
+
+
+Print["intergal s<3: ", RNIntBW[MBW, GammaBW, mDec1, mDec2, orbMomBW, scalingFacBWSmearing, 3.]];
+Print["intergal s<100000: ", RNIntBW[MBW, GammaBW, mDec1, mDec2, orbMomBW, scalingFacBWSmearing, 100000.]];
+
+
+plotIntegralBWsmax10 = Plot[{RNIntBW[MBW, GammaBW, mDec1, mDec2, orbMomBW, scalingFacBWSmearing, b]}, {b, 3., 10.}, AxesLabel->{smax [GeV GeV/(c c c c)], integral} ];
+Export["IntegralBWsmax10.pdf", plotIntegralBWsmax10];
+
+plotIntegralBWsmax100 = Plot[{RNIntBW[MBW, GammaBW, mDec1, mDec2, orbMomBW, scalingFacBWSmearing,b]}, {b, 3., 100.}, AxesLabel->{smax [GeV GeV/(c c c c)], integral} ];
+Export["IntegralBWsmax100.pdf", plotIntegralBWsmax100];
+
+plotIntegralBWsmax1000 = Plot[{RNIntBW[MBW, GammaBW, mDec1, mDec2, orbMomBW, scalingFacBWSmearing, b]}, {b, 3., 1000.}, AxesLabel->{smax [GeV GeV/(c c c c)], integral} ];
+Export["IntegralBWsmax1000.pdf", plotIntegralBWsmax1000];
+
+plotIntegralBWsmax10000 = Plot[{RNIntBW[MBW, GammaBW, mDec1, mDec2, orbMomBW, scalingFacBWSmearing, b]}, {b, 3., 10000.}, AxesLabel->{smax [GeV GeV/(c c c c)], integral} ];
+Export["IntegralBWsmax10000.pdf", plotIntegralBWsmax10000];
+
+
diff --git a/Scripts/Mathematica/CMs/SmearSingle/PlotSmearingDudek.m b/Scripts/Mathematica/CMs/SmearSingle/PlotSmearingDudek.m
new file mode 100644
index 0000000000000000000000000000000000000000..8dfcb65d6a34891041397fe811ffee417da5a33f
--- /dev/null
+++ b/Scripts/Mathematica/CMs/SmearSingle/PlotSmearingDudek.m
@@ -0,0 +1,64 @@
+(* Copyright 2025
+ Bertram Kopf (bertram@ep1.rub.de)
+ Ruhr-Universit<C3><A4>t Bochum
+
+ This file is part of Pawian.
+
+ Pawian 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 3 of the License, or
+ (at your option) any later version.
+
+ Pawian 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 Pawian. If not, see <http://www.gnu.org/licenses/>
+
+***********************************************************************************************
+ *)
+Get[FileNameJoin[{Directory[], "../KinUtils.m"}]]
+Get[FileNameJoin[{Directory[], "../Integrals.m"}]]
+(* Get[FileNameJoin[{Directory[], "../InputParamsPi1300Pi.m"}]] *)
+(* Get[FileNameJoin[{Directory[], "../InputParamsRho770Pi.m"}]] *)
+
+Get[FileNameJoin[{Directory[], "../InputParamsKstar892Pi.m"}]]
+(* Precision[1.2]; *)
+(* Precision[12/10]; *)
+
+plotIntegrFuncDudek3Mass = Plot[ {f2PoleDudek ImagSigmads2Dudek[b b, mDec1, mDec2, f2PoleDudek, m2PoleDudek, (MBW MBW), orbMomBW]},{b, (mDec1+mDec2), 3.}, AxesLabel->{M [GeV/(c c)], integr} ];
+Export["IntegrFuncDudek3Mass.pdf", plotIntegrFuncDudek3Mass];
+
+
+plotIntegrFuncDudek3 = Plot[ {f2PoleDudek ImagSigmads2Dudek[b, mDec1, mDec2, f2PoleDudek, m2PoleDudek, (MBW MBW), orbMomBW]},{b, (mDec1+mDec2)^2, 3.}, AxesLabel->{sprime [GeV/(c c)], integr} ];
+Export["IntegrFuncDudek3.pdf", plotIntegrFuncDudek3];
+
+plotIntegrFuncDudek10 = Plot[ {f2PoleDudek ImagSigmads2Dudek[b, mDec1, mDec2, f2PoleDudek, m2PoleDudek, (MBW MBW), orbMomBW]},{b, (mDec1+mDec2)^2, 10.}, AxesLabel->{sprime [GeV/(c c c c)], integr} ];
+Export["IntegrFuncDudek10.pdf", plotIntegrFuncDudek10];
+
+plotIntegrFuncDudek100 = Plot[ {f2PoleDudek ImagSigmads2Dudek[b, mDec1, mDec2, f2PoleDudek, m2PoleDudek, (MBW MBW), orbMomBW]},{b, (mDec1+mDec2)^2, 100.}, AxesLabel->{sprime [GeV GeV/(c c c c)], integr} ];
+Export["IntegrFuncDudek100.pdf", plotIntegrFuncDudek100];
+
+plotIntegrFuncDudek1000 = Plot[ {f2PoleDudek ImagSigmads2Dudek[b, mDec1, mDec2, f2PoleDudek, m2PoleDudek, (MBW MBW), orbMomBW]},{b, (mDec1+mDec2)^2, 1000.}, AxesLabel->{sprime [GeV GeV/(c c c c)], integr} ];
+Export["IntegrFuncDudek1000.pdf", plotIntegrFuncDudek1000];
+
+plotIntegrFuncDudek10000 = Plot[ {f2PoleDudek ImagSigmads2Dudek[b, mDec1, mDec2, f2PoleDudek, m2PoleDudek, (MBW MBW), orbMomBW]},{b, (mDec1+mDec2)^2, 10000.}, AxesLabel->{sprime [GeV GeV/(c c c c)], integr} ];
+Export["IntegrFuncDudek10000.pdf", plotIntegrFuncDudek10000];
+
+Print["intergal s<3: ", RNIntDudek[mDec1, mDec2, f2PoleDudek, m2PoleDudek, (MBW MBW), orbMomBW, scalingFacDudekSmearing, 3.]];
+Print["intergal s<100000: ", RNIntDudek[mDec1, mDec2, f2PoleDudek, m2PoleDudek, (MBW MBW), orbMomBW, scalingFacDudekSmearing, 100000.]];
+
+ plotIntegralDudeksmax10 = Plot[{RNIntDudek[mDec1, mDec2, f2PoleDudek, m2PoleDudek, (MBW MBW), orbMomBW, scalingFacDudekSmearing, b]}, {b, 3., 10.}, AxesLabel->{smax [GeV GeV/(c c c c)], integral} ];
+Export["IntegralDudeksmax10.pdf", plotIntegralDudeksmax10];
+
+plotIntegralDudeksmax100 = Plot[{RNIntDudek[mDec1, mDec2, f2PoleDudek, m2PoleDudek, (MBW MBW), orbMomBW, scalingFacDudekSmearing, b]}, {b, 3., 100.}, AxesLabel->{smax [GeV GeV/(c c c c)], integral} ];
+Export["IntegralDudeksmax100.pdf", plotIntegralDudeksmax100];
+
+plotIntegralDudeksmax1000 = Plot[{RNIntDudek[mDec1, mDec2, f2PoleDudek, m2PoleDudek, (MBW MBW), orbMomBW, scalingFacDudekSmearing, b]}, {b, 3., 1000.}, AxesLabel->{smax [GeV GeV/(c c c c)], integral} ];
+Export["IntegralDudeksmax1000.pdf", plotIntegralDudeksmax1000];
+
+plotIntegralDudeksmax10000 = Plot[{RNIntDudek[mDec1, mDec2, f2PoleDudek, m2PoleDudek, (MBW MBW), orbMomBW, scalingFacDudekSmearing, b]}, {b, 3., 10000.}, AxesLabel->{smax [GeV GeV/(c c c c)], integral} ];
+Export["IntegralDudeksmax10000.pdf", plotIntegralDudeksmax10000];
+