A full PWGHF simulation example for embedding

Author: preghenella

MC/config/PWGHF/external/generator/GeneratorHF.C

@@ -0,0 +1,118 @@
+/// \author R+Preghenella - July 2020
+
+// Example of an implementation of an event generator
+// that provides HF signals for embedding in background
+
+#include "Pythia8/Pythia.h"
+
+namespace o2 {
+namespace eventgen {
+
+class GeneratorHF : public GeneratorPythia8
+{
+  
+public:
+
+  GeneratorHF() : GeneratorPythia8() { };
+  ~GeneratorHF() = default;
+
+  // We initialise the local Pythia8 event where we store the particles
+  // of the signal event that is the sum of multiple Pythia8 events
+  // generated according to the generateEvent() function below.
+  Bool_t Init() override {
+    mOutputEvent.init("(GeneratorHF output event)", &mPythia.particleData);
+    return GeneratorPythia8::Init();
+  }
+  
+  // This function is called by the primary generator
+  // for each event in case we are in embedding mode.
+  // We use it to setup the number of signal events
+  // to be generated and to be embedded on the background.
+  void notifyEmbedding(const FairMCEventHeader *bkgHeader) override {
+    mEvents = mFormula.Eval(bkgHeader->GetB());
+    std::cout << " --- notify embedding: impact parameter is " << bkgHeader->GetB() << ", generating " << mEvents << " signal events " << std::endl;
+  };
+
+  // We override this function to be able to generate multiple
+  // events and build an output event that is the sum of them
+  // where we have stripped out only the sub-event starting from
+  // the c-cbar ancestor particle
+  Bool_t generateEvent() override {
+
+    // reset counter and event
+    mOutputEvent.reset();
+    
+    // loop over number of events to be generated 
+    int nEvents = 0;
+    while (nEvents < mEvents) {
+
+      // generate event
+      if (!GeneratorPythia8::generateEvent()) return false;
+
+      // find the c-cbar ancestor
+      auto ancestor = findAncestor(mPythia.event);
+      if (ancestor < 0) continue;
+
+      // append ancestor and its daughters to the output event
+      selectFromAncestor(ancestor, mPythia.event, mOutputEvent);
+      nEvents++;
+    }
+
+    if (mVerbose) mOutputEvent.list();
+    
+    return true;
+  };
+
+  // We override this event to import the particles from the
+  // output event that we have constructed as the sum of multiple
+  // Pythia8 sub-events as generated above
+  Bool_t importParticles() override {
+    return GeneratorPythia8::importParticles(mOutputEvent);
+  }
+
+  // search for c-cbar mother with at least one c at midrapidity
+  int findAncestor(Pythia8::Event& event) {
+    for (int ipa = 0; ipa < event.size(); ++ipa) {
+      auto daughterList = event[ipa].daughterList();
+      bool hasc = false, hascbar = false, atmidy = false;
+      for (auto ida : daughterList) {
+	if (event[ida].id() == 4) hasc = true;
+	if (event[ida].id() == -4) hascbar = true;
+	if (fabs(event[ida].y()) < mRapidity) atmidy = true;
+      }
+      if (hasc && hascbar && atmidy)
+	return ipa;
+    }
+    return -1;
+  };
+  
+  void setRapidity(double val) { mRapidity = val; };
+  void setVerbose(bool val) { mVerbose = val; };
+  void setFormula(std::string val) { mFormula.Compile(val.c_str()); };
+  
+private:
+
+  TFormula mFormula;
+  int mEvents = 1;
+  Pythia8::Event mOutputEvent;
+  double mRapidity = 1.5;
+  bool mVerbose = false;
+  
+};
+
+}}
+  
+/** generator instance and settings **/
+
+FairGenerator*
+GeneratorHF(double rapidity = 1.5, bool verbose = false)
+{
+  auto gen = new o2::eventgen::GeneratorHF();
+  gen->setRapidity(rapidity);
+  gen->setVerbose(verbose);
+  gen->setFormula("max(1.,120.*(x<5.)+80.*(1.-x/20.)*(x>5.)*(x<11.)+240.*(1.-x/13.)*(x>11.))");
+  
+  return gen;
+}
+
+

MC/config/PWGHF/external/trigger/.gitkeep

@@ -0,0 +1,21 @@
+### The external generator derives from GeneratorPythia8.
+### This part configures the bits of the interface: configuration and user hooks
+
+[Pythia8]
+config = ${O2DPG_ROOT}/MC/config/common/pythia8/generator/pythia8.cfg
+hooksFileName = ${O2DPG_ROOT}/MC/config/PWGHF/pythia8/hooks/pythia8_userhooks_ccbar.C
+hooksFuncName = pythia8_userhooks_ccbar(1.5)
+
+### The setup uses the base configuration of the decayer which is loaded from the file specified by config[0].
+### On top of the base configuration, two more sets of settings are loaded sequentially from config[1] and [2].
+
+[DecayerPythia8]
+config[0] = ${O2DPG_ROOT}/MC/config/common/pythia8/decayer/base.cfg
+config[1] = ${O2DPG_ROOT}/MC/config/PWGHF/pythia8/decayer/force_hadronic_D.cfg
+config[2] = ${O2DPG_ROOT}/MC/config/PWGHF/pythia8/decayer/force_hadronic_D_use4bodies.cfg
+
+### The setup forces some particles to be decayed by the external decayer instead of Geant.
+### The PDG list of the particles is specified below.
+
+[SimUserDecay]
+pdglist = 411 421 431 4112 4122 4232 4132

MC/config/PWGHF/ini/GeneratorHF.ini

@@ -0,0 +1,96 @@
+### author: Roberto Preghenella (preghenella@bo.infn.it)
+### since:  July 2020
+
+### Force golden D decay modes
+### the origin is from AliRoot AliDecayerPythia8::ForceHadronicD (Latest commit c509466 on May 17)
+###
+### This file reproduces the configuration achieved with ForceHadronicD(0,0,0)
+
+### add D+ decays absent in PYTHIA8 decay table and set BRs from PDG for other
+411:oneChannel = 1 0.0752 0 -321 211 211
+411:addChannel = 1 0.0104 0 -313 211
+411:addChannel = 1 0.0156 0 311 211
+411:addChannel = 1 0.00276 0 333 211
+## add Lc decays absent in PYTHIA8 decay table and set BRs from PDG for other
+4122:oneChannel = 1 0.0196 100 2212 -313
+4122:addChannel = 1 0.0108 100 2224 -321
+4122:addChannel = 1 0.022 100 3124 211
+4122:addChannel = 1 0.035 0 2212 -321 211
+4122:addChannel = 1 0.0159 0 2212 311
+4122:addChannel = 1 0.0130 0 3122 211
+### add Xic+ decays absent in PYTHIA8 decay table
+4232:addChannel = 1 0.2 0 2212 -313
+4232:addChannel = 1 0.2 0 2212 -321 211
+4232:addChannel = 1 0.2 0 3324 211
+4232:addChannel = 1 0.2 0 3312 211 211
+### add Xic0 decays absent in PYTHIA8 decay table
+4132:addChannel = 1 0.2 0 3312 211
+
+### K* -> K pi
+313:onMode = off
+313:onIfAll = 321 211
+### for Ds -> Phi pi+
+333:onMode = off
+333:onIfAll = 321 321
+### for D0 -> rho0 pi+ k-
+113:onMode = off
+113:onIfAll = 211 211
+### for Lambda_c -> Delta++ K-
+2224:onMode = off
+2224:onIfAll = 2212 211
+### for Lambda_c -> Lambda(1520) K-
+3124:onMode = off
+3124:onIfAll = 2212 321
+
+### Omega_c -> Omega pi
+4332:onMode = off
+4332:onIfMatch = 3334 211
+
+### switch off all decay channels
+411:onMode = off
+421:onMode = off
+431:onMode = off
+4112:onMode = off
+4122:onMode = off
+4232:onMode = off
+4132:onMode = off
+
+### D+/- -> K pi pi 
+411:onIfMatch = 321 211 211
+### D+/- -> K* pi
+411:onIfMatch = 313 211
+### D+/- -> phi pi
+411:onIfMatch = 333 211
+
+### D0 -> K pi
+421:onIfMatch = 321 211
+
+### D_s -> K K*
+431:onIfMatch = 321 313
+### D_s -> Phi pi
+431:onIfMatch = 333 211
+
+### Lambda_c -> p K*
+4122:onIfMatch = 2212 313
+### Lambda_c -> Delta K
+4122:onIfMatch = 2224 321
+### Lambda_c -> Lambda(1520) pi
+4122:onIfMatch = 3124 211
+### Lambda_c -> p K pi
+4122:onIfMatch = 2212 321 211
+### Lambda_c -> Lambda pi 
+4122:onIfMatch = 3122 211
+### Lambda_c -> p K0
+4122:onIfMatch = 2212 311
+
+### Xic+ -> pK*0
+4232:onIfMatch = 2212 313
+### Xic+ -> p K- pi+
+4232:onIfMatch = 2212 321 211
+### Xic+ -> Xi*0 pi+, Xi*->Xi- pi+
+4232:onIfMatch = 3324 211
+### Xic+ -> Xi- pi+ pi+
+4232:onIfMatch = 3312 211 211
+    
+### Xic0 -> Xi- pi+
+4132:onIfMatch = 3312 211

MC/config/PWGHF/pythia8/decayer/force_hadronic_D.cfg

@@ -0,0 +1,12 @@
+### author: Roberto Preghenella (preghenella@bo.infn.it)
+### since:  July 2020
+
+### Force golden D decay modes (force Lc channel 1)
+### the origin is from AliRoot AliDecayerPythia8::ForceHadronicD (Latest commit c509466 on May 17)
+###
+### This file has to be used in conjunction with force_hadronic_D.cfg and loaded
+### afterwards to reproduce the configuration achieved with ForceHadronicD(0,0,1)
+
+### force Lc -> p K pi
+4122:onMode = off
+4122:onIfMatch = 2212 321 211

MC/config/PWGHF/pythia8/decayer/force_hadronic_D_forceLcChannel1.cfg

@@ -0,0 +1,12 @@
+### author: Roberto Preghenella (preghenella@bo.infn.it)
+### since:  July 2020
+
+### Force golden D decay modes (force Lc channel 2)
+### the origin is from AliRoot AliDecayerPythia8::ForceHadronicD (Latest commit c509466 on May 17)
+###
+### This file has to be used in conjunction with force_hadronic_D.cfg and loaded
+### afterwards to reproduce the configuration achieved with ForceHadronicD(0,0,2)
+
+### force Lc -> p K0s
+4122:onMode = off
+4122:onIfMatch = 2212 311

MC/config/PWGHF/pythia8/decayer/force_hadronic_D_forceLcChannel2.cfg

@@ -0,0 +1,15 @@
+### author: Roberto Preghenella (preghenella@bo.infn.it)
+### since:  July 2020
+
+### Force golden D decay modes (use 4 bodies option)
+### the origin is from AliRoot AliDecayerPythia8::ForceHadronicD (Latest commit c509466 on May 17)
+###
+### This file has to be used in conjunction with force_hadronic_D.cfg and loaded
+### afterwards to reproduce the configuration achieved with ForceHadronicD(1,0,0)
+
+### D0 -> K pi pi pi
+421:onIfMatch = 321 211 211 211
+### D0 -> K pi rho
+421:onIfMatch = 321 211 113
+### D0 -> K*0 pi pi
+421:onIfMatch = 313 211 211

MC/config/PWGHF/pythia8/decayer/force_hadronic_D_use4bodies.cfg

@@ -0,0 +1,13 @@
+### author: Roberto Preghenella (preghenella@bo.infn.it)
+### since:  July 2020
+
+### Force golden D decay modes (use D to V0 option)
+### the origin is from AliRoot AliDecayerPythia8::ForceHadronicD (Latest commit c509466 on May 17)
+###
+### This file has to be used in conjunction with force_hadronic_D.cfg and loaded
+### afterwards to reproduce the configuration achieved with ForceHadronicD(0,1,0)
+
+### Ds -> K0K
+431:onIfMatch = 311 321
+### Ds -> K0pi
+411:onIfMatch = 311 211

MC/config/PWGHF/pythia8/decayer/force_hadronic_D_useDtoV0.cfg

@@ -0,0 +1,47 @@
+/// \author R+Preghenella - July 2020
+
+/// This Pythia8 UserHooks can veto the processing at parton level.
+/// The partonic event is scanned searching for a c-cbar mother
+/// with at least one of the c quarks produced withing a fiducial
+/// window around midrapidity that can be specified by the user.
+
+#include "Pythia8/Pythia.h"
+
+class UserHooks_ccbar : public Pythia8::UserHooks
+{
+  
+ public:
+  UserHooks_ccbar() = default;
+  ~UserHooks_ccbar() = default;
+  bool canVetoPartonLevel() override { return true; };
+  bool doVetoPartonLevel(const Pythia8::Event& event) override {
+    // search for c-cbar mother with at least one c at midrapidity
+    for (int ipa = 0; ipa < event.size(); ++ipa) {
+      auto daughterList = event[ipa].daughterList();
+      bool hasc = false, hascbar = false, atmidy = false;
+      for (auto ida : daughterList) {
+	if (event[ida].id() == 4) hasc = true;
+	if (event[ida].id() == -4) hascbar = true;
+	if (fabs(event[ida].y()) < mRapidity) atmidy = true;
+      }
+      if (hasc && hascbar && atmidy)
+	return false; // found it, do not veto event
+    }
+    return true; // did not find it, veto event
+  };
+
+  void setRapidity(double val) { mRapidity = val; };
+  
+private:
+
+  double mRapidity = 1.5;
+  
+};
+
+Pythia8::UserHooks*
+  pythia8_userhooks_ccbar(double rapidity = 1.5)
+{
+  auto hooks = new UserHooks_ccbar();
+  hooks->setRapidity(rapidity);
+  return hooks;
+}