diff --git a/PWGCF/Femto/Core/closePairRejection.h b/PWGCF/Femto/Core/closePairRejection.h index 9460fad49dc..51470519c2d 100644 --- a/PWGCF/Femto/Core/closePairRejection.h +++ b/PWGCF/Femto/Core/closePairRejection.h @@ -26,6 +26,8 @@ #include #include +#include + #include #include #include @@ -69,6 +71,8 @@ struct ConfCpr : o2::framework::ConfigurableGroup { o2::framework::Configurable plotAngularCorrelation{"plotAngularCorrelation", false, "Plot angular correlation of particles (eta1 vs eta2 & phi1 vs phi2"}; o2::framework::Configurable detaMax{"detaMax", 0.01f, "Maximium deta"}; o2::framework::Configurable dphistarMax{"dphistarMax", 0.01f, "Maximum dphistar"}; + o2::framework::Configurable ktDependent{"ktDependent", false, "If true, use dphistarMaxFormula instead of dphistarMax"}; + o2::framework::Configurable dphistarMaxFormula{"dphistarMaxFormula", "0.01", "Maximum dphistar as TFormula, x=kT"}; o2::framework::Configurable detaCenter{"detaCenter", 0.f, "Center of deta cut"}; o2::framework::Configurable dphistarCenter{"dphistarCenter", 0.f, "Center of dphistar cut"}; o2::framework::Configurable kinematicMin{"kinematicMin", -1.f, "Minimum kstar/Q3 of pair/triplet for plotting (Set to negative value to turn off the cut)"}; @@ -191,9 +195,14 @@ class CloseTrackRejection { mDetaMax = confCpr.detaMax.value; mDphistarMax = confCpr.dphistarMax.value; + mKtDependent = confCpr.ktDependent.value; + + if (mKtDependent) { + mDphistarMaxFunc = TF1((std::string(prefix) + "DphistarMax").c_str(), confCpr.dphistarMaxFormula.value.c_str(), 0.f, 10.f); + } // check the limits - if (mDetaMax <= 0 || mDphistarMax <= 0) { + if (mDetaMax <= 0 || (!mKtDependent && mDphistarMax <= 0)) { LOG(fatal) << "Limits for Close Pair Rejection are invalid (0 or negative). Breaking..."; } @@ -292,6 +301,12 @@ class CloseTrackRejection mAverageDphistar = 0.f; // if computation at all radii fail, set it 0 } + if (mKtDependent) { + const float ktx = t1.pt() * std::cos(t1.phi()) + t2.pt() * std::cos(t2.phi()); + const float kty = t1.pt() * std::sin(t1.phi()) + t2.pt() * std::sin(t2.phi()); + mDphistarMax = mDphistarMaxFunc.Eval(0.5f * std::hypot(ktx, kty)); + } + if (mPlotAngularCorrelation) { mPhi1 = t1.phi(); mPhi2 = t2.phi(); @@ -414,6 +429,8 @@ class CloseTrackRejection float mDphistarMax = 0.f; float mDetaCenter = 0.f; float mDphistarCenter = 0.f; + bool mKtDependent = false; + TF1 mDphistarMaxFunc; float mAverageDphistar = 0.f; float mDeta = 0.f;