tunable ratio of ccbar and bbbar for embedding

[wuctlby]

Hi @xinyepeng , @gluparel , @stefanopolitano , @zhangbiao-phy , and @mfaggin ,

As discussed in the previous D2H meeting, I implemented a new generator with a tunable ratio of ccbar and bbbar for embedding the HF event into the signal event and the corresponding ini file, as well as the modification of 0.36000 0 2212 311 ### Λc+ -> p K0S 1.59% ==> 1 0.36000 0 2212 310 ### Λc+ -> p K0S 1.59% in the Pythia config file.

Specifically,

• generator_pythia8_embed_hf.C
  -- new hf_generators enum GapHFRatio to trigger GeneratorPythia8GapTriggeredHF with a given quarkPdgList
  -- BuildQuarkListFromBOverC, which builds the quarkPdgList from bOverCRatio
  -- GeneratorPythia8EmbedHFRatio, which accepts the new parameter bOverCRatio for generating the quarkPdgList. For example, if bOverCRatio=3, the BuildQuarkListFromBOverC will generate a list [5,5,5,4], which will embed the ccbar and bbbar in the order of this list

Please let me know if you have any comments or suggestions. Thanks a lot!

[github-actions]

REQUEST FOR PRODUCTION RELEASES:
To request your PR to be included in production software, please add the corresponding labels called "async-" to your PR. Add the labels directly (if you have the permissions) or add a comment of the form (note that labels are separated by a ",")

+async-label <label1>, <label2>, !<label3> ...

This will add <label1> and <label2> and removes <label3>.

The following labels are available
async-2023-pbpb-apass4
async-2023-pp-apass4
async-2024-pp-apass1
async-2022-pp-apass7
async-2024-pp-cpass0
async-2024-PbPb-apass1
async-2024-ppRef-apass1
async-2024-PbPb-apass2
async-2023-PbPb-apass5

[mfaggin]

Suggested change

	if (bOverCRatio <= 0.f) {
	LOG(fatal) << "bOverCRatio (b/c) must be > 0";
	}
	if (bOverCRatio > 19.f) {
	bOverCRatio = 19.f;
	LOG(warn) << "bOverCRatio (b/c) too large, using 19:1";
	}else if (bOverCRatio < 1.f/19.f) {
	bOverCRatio = 1.f/19.f;
	LOG(warn) << "bOverCRatio (b/c) too small, using 1:19";
	}
	int nC = 1, nB = 1;
	float bestErr = 1e9f;
	for (int c = 1; c <= 19; ++c) {
	int b = (int)std::lround(bOverCRatio * c);
	if (b < 1 || c + b > 20) {
	continue;
	}
	float err = std::fabs((float)b / c - bOverCRatio);
	if (err < bestErr) {
	bestErr = err;
	nC = c;
	nB = b;
	}
	}
	const int iterNMax = 19;
	if (bOverCRatio <= 0.f) {
	LOG(fatal) << "bOverCRatio (b/c) must be > 0";
	}
	if (bOverCRatio > iterNMax*1.f) {
	bOverCRatio = iterNMax*1.f;
	LOG(warn) << "bOverCRatio (b/c) too large, using 19:1";
	}else if (bOverCRatio < 1.f/iterNMax) {
	bOverCRatio = 1.f/iterNMax;
	LOG(warn) << "bOverCRatio (b/c) too small, using 1:19";
	}
	int nC = 1, nB = 1;
	float bestErr = 1e9f;
	for (int c = 1; c <= iterNMax; ++c) {
	int b = static_cast<int>(std::lround(bOverCRatio * c));
	if (b < 1 || c + b > 20) {
	continue;
	}
	float err = std::fabs(static_cast<float>(b) / c - bOverCRatio);
	if (err < bestErr) {
	/// This check here is needed in case bOverCRatio*c is not integer (it can happen with e.g. bOverCRatio=4./9.)
	/// In this case, b is its truncation and the desired ratio is not obtained
	/// It means that one needs to continue iterating until bOverCRatio*c is integer, i.e. b not truncated
	///
	/// Possible cases:
	/// 1. we want nB = nC*R, with R integer
	/// -> we enter here in the first loop iteration
	/// 2. we want more charm than beauty by an integer amount, i.e. nB = nC*R with R=1./N, with N integer
	/// -> we enter here after N iterations, when c=N
	/// 3. we want either more charm or beauty, but with a factor R that is not integer, as well as its inverse (e.g. bOverCRatio=4./9.)
	/// -> In this case, bOverCRatio*c is not integer, namely b is its truncation and the desired ratio is not obtained
	/// The code iterates at most until c becomes equal to the denominator of the fraction
	/// 4. we want one of the previous cases, but we assign to bOverCRatio a value that is not rational
	/// or such as we do not enter here within 19 iterations
	/// -> nC and nB, are not touched, therefore we do not have the desired fraction. We'll need to throw a fatal (*)
	///
	bestErr = err;
	nC = c;
	nB = b;
	/// If we are at this point, we reached already the desired ratio between b and c.
	/// Let's break the loop
	break;
	}
	}
	// (*) check if we have the desired fraction
	bool isRatioUnity = std::fabs(bOverCRatio-1) < 1e-05;
	if (!isRatioUnity && nC==1 && nB==1) {
	LOG(fatal) << "nC=" << nC << ", nB=" << nB << " but you ask bOverCRatio to be " << bOverCRatio <<", which is different from nB/nC. It means either that bOverCRatio is not rational, or that you need more than " << iterNMax << " iterations. Change it!";
	}

[mfaggin]

this approach here is nice, but in my opinion as written here it can be improved.

1. The loop can be stopped the first time we have c and b such that their ratio corresponds to the desired one, namely we need to break; the loop to avoid useless iterations.
   In fact, with the current implementation always 19 iterations are done, even if only 1 is needed.

Your code gives what follows, assuming bOverCRatio=3.f:

iter
1          c=1, b=3   OK -> nC=1, nB=3    <-- we can stop here already
2          c=2, b=6   OK -> nC=2, nB=6
3          c=3, b=9   OK -> nC=3, nB=9
4          c=4, b=12  OK -> nC=4, nB=12
5          c=5, b=15  OK -> nC=5, nB=15   <-- this is the last time nC and nB are updated
6          c=6, b=18  NOT OK (c + b > 20)
7          c=7, b=21  NOT OK (c + b > 20)
.
.
.
19         c=19, b=57  NOT OK (c + b > 20)

After this loop you have then nC=5 and nB=15, therefore you create a vector of 20 elements equal to {5,5,5,4,5,5,5,4,5,5,5,4,5,5,5,4}.
I propose to add a break; to solve this, avoiding the useless iterations and reducing in this case to have the vector {5,5,5,4}, which is enough for our purpose.

2. The code does not crash in case the the desired ratio is not guaranteed. This can happen in case as input one puts a value that is not rational, and for (very) particular values of this ratio that do not allow to converge within 19 iterations (I do not dislike the idea to have a max. number of iterations, to keep the running time limited). In this case we need the code to crash and ask the user to change the value, since it's very likely that the input value is not correct. I propose to add this check after the loop.

3. I propose to add some comments to make the code more readable

I propose these code changes, it would be good if you can give a few trials and let me know what you think

[mfaggin]

In principle this should work in pythia, but did you check that this works?
The physics wants that only 4122:addChannel = 1 0.36000 0 2212 311 is correct, since K0, K0bar are the only states with a well-defined strangeness content, that materialize as K0s or K0L only in decay.
Did you check that the correct and expected decays are present, i.e. the fact that we have Lc decaying into p and K0s with K0s decaying in 2 pions?

[mfaggin]

thanks @wuctlby , I left a couple of comments with a few suggestions and a question. Thank you!