HF cocktail with forced decays

MC/config/PWGEM/external/generator/Generator_pythia8_forcedDecays.C

@@ -0,0 +1,151 @@
+#include "Generators/DecayerPythia8.h"
+#include "Generators/GeneratorPythia8.h"
+#include "SimulationDataFormat/MCGenProperties.h"
+#include "SimulationDataFormat/ParticleStatus.h"
+#include "TLorentzVector.h"
+
+namespace o2 {
+namespace eventgen {
+
+class DecayerPythia8ForceDecays : public DecayerPythia8 {
+public:
+  DecayerPythia8ForceDecays() : DecayerPythia8(){};
+  ~DecayerPythia8ForceDecays() = default;
+
+  void calculateWeights(std::vector<int> &pdgs) {
+    TLorentzVector mom = TLorentzVector(0., 0., 0., 9999999.);
+    for (int pdg : pdgs) {
+      Decay(pdg, &mom); // do one fake decay to initalize everything correctly
+      auto particleData = mPythia.particleData.particleDataEntryPtr(pdg);
+      double weight = 0.;
+      for (int c = 0; c < particleData->sizeChannels(); c++) {
+        weight += particleData->channel(c).currentBR();
+      }
+      LOG(info) << "PDG = " << pdg
+                << ": sum of branching ratios of active decay channels = "
+                << weight;
+      mWeights[pdg] = weight;
+      mPythia.particleData.mayDecay(
+          pdg, false); // make sure it is not decayed in any other way (should
+                       // already be switched off in the decayer config file)
+    }
+  }
+
+  void forceDecays(std::vector<TParticle> &mParticles, int mother_pos) {
+    TParticle p = mParticles[mother_pos];
+    int pdg = p.GetPdgCode();
+    TLorentzVector mom = TLorentzVector(p.Px(), p.Py(), p.Pz(), p.Energy());
+    Decay(pdg, &mom);
+    TClonesArray daughters = TClonesArray("TParticle");
+    int nParticles = ImportParticles(&daughters);
+    p.SetStatusCode(o2::mcgenstatus::MCGenStatusEncoding(2, -11).fullEncoding);
+    p.SetBit(ParticleStatus::kToBeDone, false);
+    double mother_weight = p.GetWeight();
+    TParticle *mother = static_cast<TParticle *>(daughters[0]);
+    int mParticles_size = mParticles.size();
+    p.SetFirstDaughter(mother->GetFirstDaughter() + mParticles_size - 1);
+    p.SetLastDaughter(mother->GetLastDaughter() + mParticles_size - 1);
+    for (int j = 1; j < nParticles;
+         j++) { // start loop at 1 to not include mother
+      TParticle *d = static_cast<TParticle *>(daughters[j]);
+      double decay_weight = mWeights[abs(pdg)];
+      if (decay_weight == 0) {
+        LOG(error) << "Decaying particle (PDG = " << pdg
+                   << ") with decay weight = 0. Did you set the pdg codes for "
+                      "calculating weights correctly?";
+      }
+      d->SetWeight(decay_weight * mother_weight);
+      if (d->GetStatusCode() == 1) {
+        d->SetStatusCode(
+            o2::mcgenstatus::MCGenStatusEncoding(1, 91).fullEncoding);
+        d->SetBit(ParticleStatus::kToBeDone, true);
+      } else {
+        d->SetStatusCode(
+            o2::mcgenstatus::MCGenStatusEncoding(2, -91).fullEncoding);
+        d->SetBit(ParticleStatus::kToBeDone, false);
+      }
+      int firstmother = d->GetFirstMother();
+      int firstdaughter = d->GetFirstDaughter();
+      int lastdaughter = d->GetLastDaughter();
+      if (firstmother == 0) {
+        d->SetFirstMother(mother_pos);
+        d->SetLastMother(mother_pos);
+      } else {
+        d->SetFirstMother(firstmother + mParticles_size - 1);
+        d->SetLastMother(firstmother + mParticles_size - 1);
+      }
+      if (firstdaughter == 0) {
+        d->SetFirstDaughter(-1);
+      } else {
+        d->SetFirstDaughter(firstdaughter + mParticles_size - 1);
+      }
+      if (lastdaughter == 0) {
+        d->SetLastDaughter(-1);
+      } else {
+        d->SetLastDaughter(lastdaughter + mParticles_size - 1);
+      }
+      mParticles.push_back(*d);
+    }
+  }
+
+private:
+  std::map<int, double> mWeights;
+};
+
+class GeneratorPythia8ForcedDecays : public GeneratorPythia8 {
+
+public:
+  GeneratorPythia8ForcedDecays() : GeneratorPythia8(){};
+  ~GeneratorPythia8ForcedDecays() = default;
+
+  // overriden methods
+  bool Init() override { return GeneratorPythia8::Init() && InitDecayer(); };
+  bool importParticles() override {
+    return GeneratorPythia8::importParticles() && makeForcedDecays();
+  };
+
+  void setPDGs(TString pdgs) {
+    TObjArray *obj = pdgs.Tokenize(";");
+    for (int i = 0; i < obj->GetEntriesFast(); i++) {
+      std::string spdg = obj->At(i)->GetName();
+      int pdg = std::stoi(spdg);
+      mPdgCodes.push_back(pdg);
+      LOG(info) << "Force decay of PDG = " << pdg;
+    }
+  }
+
+protected:
+  bool InitDecayer() {
+    mDecayer = new DecayerPythia8ForceDecays();
+    mDecayer->Init();
+    mDecayer->calculateWeights(mPdgCodes);
+    return true;
+  }
+
+  bool makeForcedDecays() {
+    int mParticles_size = mParticles.size();
+    for (int i = 0; i < mParticles_size; i++) {
+      auto p = mParticles[i];
+      int pdg = p.GetPdgCode();
+      if (std::find(mPdgCodes.begin(), mPdgCodes.end(), abs(pdg)) !=
+          mPdgCodes.end()) {
+        mDecayer->forceDecays(mParticles, i);
+        mParticles_size = mParticles.size();
+      }
+    }
+    return true;
+  }
+
+private:
+  DecayerPythia8ForceDecays *mDecayer;
+  std::vector<int> mPdgCodes;
+};
+
+} // namespace eventgen
+} // namespace o2
+
+FairGenerator *GeneratePythia8ForcedDecays(TString pdgs) {
+  auto gen = new o2::eventgen::GeneratorPythia8ForcedDecays();
+  gen->setPDGs(pdgs);
+  return gen;
+}

MC/config/PWGEM/ini/Pythia8_Beauty_Cocktail.ini

@@ -0,0 +1,12 @@
+[GeneratorExternal]
+fileName = ${O2DPG_ROOT}/MC/config/PWGEM/external/generator/Generator_pythia8_forcedDecays.C
+funcName=GeneratePythia8ForcedDecays("411;421;431;4122;4232;4132;4332;511;521;531;5122;5132;5232;5332")
+
+[GeneratorPythia8]
+config = ${O2DPG_ROOT}/MC/config/PWGEM/pythia8/generator/pythia8_hf_cocktail.cfg
+hooksFileName = ${O2DPG_ROOT}/MC/config/PWGHF/pythia8/hooks/pythia8_userhooks_qqbar.C
+hooksFuncName = pythia8_userhooks_bbbar(-9999.,9999.)
+
+[DecayerPythia8]
+config[0] = ${O2DPG_ROOT}/MC/config/common/pythia8/decayer/base.cfg
+config[1] = ${O2DPG_ROOT}/MC/config/PWGEM/pythia8/decayer/force_semileptonic.cfg

MC/config/PWGEM/ini/Pythia8_Charm_Cocktail.ini

@@ -0,0 +1,12 @@
+[GeneratorExternal]
+fileName = ${O2DPG_ROOT}/MC/config/PWGEM/external/generator/Generator_pythia8_forcedDecays.C
+funcName=GeneratePythia8ForcedDecays("411;421;431;4122;4232;4132;4332")
+
+[GeneratorPythia8]
+config = ${O2DPG_ROOT}/MC/config/PWGEM/pythia8/generator/pythia8_hf_cocktail.cfg
+hooksFileName = ${O2DPG_ROOT}/MC/config/PWGHF/pythia8/hooks/pythia8_userhooks_qqbar.C
+hooksFuncName = pythia8_userhooks_ccbar(-9999.,9999.)
+
+[DecayerPythia8]
+config[0] = ${O2DPG_ROOT}/MC/config/common/pythia8/decayer/base.cfg
+config[1] = ${O2DPG_ROOT}/MC/config/PWGEM/pythia8/decayer/force_semileptonic.cfg

MC/config/PWGEM/pythia8/decayer/force_semileptonic.cfg

@@ -0,0 +1,80 @@
+### only semileptonic decays for charm
+### D+
+411:oneChannel = 1 0.087 0 -311 -11 12
+411:addChannel = 1 0.040 0 -321 211 -11 12
+411:addChannel = 1 0.037 0 -313 -11 12
+### D0
+421:oneChannel = 1 0.035 0 -321 -11 12
+421:addChannel = 1 0.022 0 -323 -11 12
+421:addChannel = 1 0.016 0 -321 111 -11 12
+421:addChannel = 1 0.014 0 -311 -211 -11 12
+### Ds
+431:oneChannel = 1 0.025 0 333 -11 12
+431:addChannel = 1 0.027 0 221 -11 12
+### Lambdac
+4122:oneChannel = 1 0.036 0 3122 -11 12
+### chi_{c}^{+}
+4232:oneChannel = 1 0.07 0 3322 -11 12
+### chi_{c}^{0}
+4132:oneChannel = 1 0.014 0 3312 -11 12
+### Omega_{c}
+4332:oneChannel = 1 0.01224 0 3334 -11 12
+
+### only semileptonic decays for beauty
+### B0
+511:oneChannel = 1 0.0207000 0 12 -11 -411
+511:addChannel = 1 0.0570000 0 12 -11 -413
+511:addChannel = 1 0.0023000 0 12 -11 -415
+511:addChannel = 1 0.0001330 0 12 -11 -211
+511:addChannel = 1 0.0002690 0 12 -11 -213
+511:addChannel = 1 0.0045000 0 12 -11 -10411
+511:addChannel = 1 0.0052000 0 12 -11 -10413
+511:addChannel = 1 0.0083000 0 12 -11 -20413
+
+### B+
+521:oneChannel = 1 0.0000720 0 12 -11 111
+521:addChannel = 1 0.0001450 0 12 -11 113
+521:addChannel = 1 0.0000840 0 12 -11 221
+521:addChannel = 1 0.0001450 0 12 -11 223
+521:addChannel = 1 0.0000840 0 12 -11 331
+521:addChannel = 1 0.0224000 0 12 -11 -421
+521:addChannel = 1 0.0617000 0 12 -11 -423
+521:addChannel = 1 0.0030000 0 12 -11 -425
+521:addChannel = 1 0.0049000 0 12 -11 -10421
+521:addChannel = 1 0.0056000 0 12 -11 -10423
+521:addChannel = 1 0.0090000 0 12 -11 -20423
+
+### Bs
+531:oneChannel = 1 0.0002000 0 12 -11 -321
+531:addChannel = 1 0.0003000 0 12 -11 -323
+531:addChannel = 1 0.0210000 0 12 -11 -431
+531:addChannel = 1 0.0490000 0 12 -11 -433
+531:addChannel = 1 0.0070000 0 12 -11 -435
+531:addChannel = 1 0.0003000 0 12 -11 -10323
+531:addChannel = 1 0.0040000 0 12 -11 -10431
+531:addChannel = 1 0.0070000 0 12 -11 -10433
+531:addChannel = 1 0.0002000 0 12 -11 -20323
+531:addChannel = 1 0.0040000 0 12 -11 -20433
+
+### Lambdab
+5122:oneChannel = 1 0.0546000 0 -12 11 4122
+5122:addChannel = 1 0.0096000 0 -12 11 4124
+5122:addChannel = 1 0.0128000 0 -12 11 14122
+
+### Chi_{b}^{-}
+5132:oneChannel = 1 0.1080010 0 -12 11 4 3101
+5132:addChannel = 1 0.0020000 0 -12 11 2 3101
+### Chi_{b}^{0}
+5232:oneChannel = 1 0.1080010 0 -12 11 4 3201
+5232:addChannel = 1 0.0020000 0 -12 11 2 3201
+### Omega_{b}^{-}
+5332:oneChannel = 1 0.1080010 1 -12 11 4 3303
+5332:oneChannel = 1 0.0020000 1 -12 11 2 3303
+
+
+# Correct OmegaC decay length (wrong in PYTHIA8 decay table) (mm/c)
+4332:tau0 = 0.08000000000
+# Correct Lb decay length (wrong in PYTHIA8 decay table)
+5122:tau0 = 4.41000e-01
+
+

MC/config/PWGEM/pythia8/generator/pythia8_hf_cocktail.cfg

@@ -0,0 +1,48 @@
+### beams
+Beams:idA 2212			# proton
+Beams:idB 2212                  # proton
+Beams:eCM 13600.                # GeV
+
+### processes
+SoftQCD:inelastic on
+
+### decays
+ParticleDecays:limitTau0 on
+ParticleDecays:tau0Max 10.
+
+### switch on color reconnection in mode 2 (https://arxiv.org/pdf/1505.01681.pdf)
+Tune:pp = 14
+ColourReconnection:mode = 1
+ColourReconnection:allowDoubleJunRem = off
+ColourReconnection:m0 = 0.3
+ColourReconnection:allowJunctions = on
+ColourReconnection:junctionCorrection = 1.20
+ColourReconnection:timeDilationMode = 2
+ColourReconnection:timeDilationPar = 0.18
+StringPT:sigma = 0.335
+StringZ:aLund = 0.36
+StringZ:bLund = 0.56
+StringFlav:probQQtoQ = 0.078
+StringFlav:ProbStoUD = 0.2
+StringFlav:probQQ1toQQ0join = 0.0275,0.0275,0.0275,0.0275
+MultiPartonInteractions:pT0Ref = 2.15
+BeamRemnants:remnantMode = 1
+BeamRemnants:saturation =5
+
+### switch off decays
+411:mayDecay  off   # D+
+421:mayDecay  off   # D0
+431:mayDecay  off   # D_s
+4122:mayDecay off   # Lambda_c
+4232:mayDecay off   # Xi_c
+4132:mayDecay off   # Xi_c_0
+4332:mayDecay off   # Omega_c
+511:mayDecay  off	# B0
+521:mayDecay  off	# B+
+531:mayDecay  off	# B_s0
+541:mayDecay  off	# B_c+
+5112:mayDecay off	# Sigma_b-
+5122:mayDecay off	# Lambda_b0
+5132:mayDecay off	# Xi_b-
+5232:mayDecay off	# Xi_b0
+5332:mayDecay off	# Omega_b-