Adding Upsilon Generator

MC/config/PWGDQ/external/generator/GeneratorBottomonia.C

@@ -0,0 +1,280 @@
+//
+// generators for bottomonia considering at midrapidity and forward rapidity
+//
+
+R__ADD_INCLUDE_PATH($O2DPG_ROOT/MC/config/PWGDQ/EvtGen)
+R__LOAD_LIBRARY(libpythia6)
+#include "GeneratorCocktail.C"
+#include "GeneratorEvtGen.C"
+
+namespace o2
+{
+namespace eventgen
+{
+
+/////////////////////////////////////////////////////////////////////////////
+class O2_GeneratorParamUpsilon1SFwdY : public GeneratorTGenerator
+{
+
+ public:
+  O2_GeneratorParamUpsilon1SFwdY() : GeneratorTGenerator("ParamUpsilon1S")
+  {
+    paramUpsilon1S = new GeneratorParam(1, -1, PtUpsilon1Spp13TeV, YUpsilon1Spp13TeV, V2Upsilon1Spp13TeV, IpUpsilon1Spp13TeV);
+    paramUpsilon1S->SetMomentumRange(0., 1.e6);
+    paramUpsilon1S->SetPtRange(0, 999.);
+    paramUpsilon1S->SetYRange(-4.2, -2.3);
+    paramUpsilon1S->SetPhiRange(0., 360.);
+    paramUpsilon1S->SetDecayer(new TPythia6Decayer());
+    paramUpsilon1S->SetForceDecay(kNoDecay); // particle left undecayed
+    // - - paramUpsilon1S->SetTrackingFlag(1);  // (from AliGenParam) -> check this
+    setTGenerator(paramUpsilon1S);
+  };
+
+  ~O2_GeneratorParamUpsilon1SFwdY()
+  {
+    delete paramUpsilon1S;
+  };
+
+  Bool_t Init() override
+  {
+    GeneratorTGenerator::Init();
+    paramUpsilon1S->Init();
+    return true;
+  }
+
+  void SetNSignalPerEvent(Int_t nsig) { paramUpsilon1S->SetNumberParticles(nsig); }
+
+  //-------------------------------------------------------------------------//
+  static Double_t PtUpsilon1Spp13TeV(const Double_t* px, const Double_t* /*dummy*/)
+  {
+    // Upsilon1S pt shape from LHCb pp@13TeV arXiv:1804.09214
+    Double_t x = *px;
+    Float_t p0, p1, p2, p3;
+
+    p0 = 4.11195e+02;
+    p1 = 1.03097e+01;
+    p2 = 1.62309e+00;
+    p3 = 4.84709e+00;
+
+    return (p0 * x / TMath::Power(1. + TMath::Power(x / p1, p2), p3));
+  }
+
+  //-------------------------------------------------------------------------//
+  static Double_t YUpsilon1Spp13TeV(const Double_t* py, const Double_t* /*dummy*/)
+  {
+    // Upsilon1S y shape from LHCb pp@13TeV arXiv:1804.09214
+    Double_t x = *py;
+    Float_t p0, p1;
+
+    p0 =  3.07931e+03;
+    p1 = -3.53102e-02;
+
+    return (p0 * (1. + p1 * x * x));
+  }
+
+  //-------------------------------------------------------------------------//
+  static Double_t V2Upsilon1Spp13TeV(const Double_t* /*dummy*/, const Double_t* /*dummy*/)
+  {
+    // Upsilon(1S) v2
+    return 0.;
+  }
+
+  //-------------------------------------------------------------------------//
+  static Int_t IpUpsilon1Spp13TeV(TRandom*)
+  {
+    return 553;
+  }
+
+ private:
+  GeneratorParam* paramUpsilon1S = nullptr;
+};
+
+/////////////////////////////////////////////////////////////////////////////
+class O2_GeneratorParamUpsilon2SFwdY : public GeneratorTGenerator
+{
+
+ public:
+  O2_GeneratorParamUpsilon2SFwdY() : GeneratorTGenerator("ParamUpsilon2S")
+  {
+    paramUpsilon2S = new GeneratorParam(1, -1, PtUpsilon2Spp13TeV, YUpsilon2Spp13TeV, V2Upsilon2Spp13TeV, IpUpsilon2Spp13TeV);
+    paramUpsilon2S->SetMomentumRange(0., 1.e6);
+    paramUpsilon2S->SetPtRange(0, 999.);
+    paramUpsilon2S->SetYRange(-4.2, -2.3);
+    paramUpsilon2S->SetPhiRange(0., 360.);
+    paramUpsilon2S->SetDecayer(new TPythia6Decayer());
+    paramUpsilon2S->SetForceDecay(kNoDecay); // particle left undecayed
+    // - - paramUpsilon2S->SetTrackingFlag(1);  // (from AliGenParam) -> check this
+    setTGenerator(paramUpsilon2S);
+  };
+
+  ~O2_GeneratorParamUpsilon2SFwdY()
+  {
+    delete paramUpsilon2S;
+  };
+
+  Bool_t Init() override
+  {
+    GeneratorTGenerator::Init();
+    paramUpsilon2S->Init();
+    return true;
+  }
+
+  void SetNSignalPerEvent(Int_t nsig) { paramUpsilon2S->SetNumberParticles(nsig); }
+
+  //-------------------------------------------------------------------------//
+  static Double_t PtUpsilon2Spp13TeV(const Double_t* px, const Double_t* /*dummy*/)
+  {
+    // Upsilon2S pt shape from LHCb pp@13TeV arXiv:1804.09214
+    Double_t x = *px;
+    Float_t p0, p1, p2, p3;
+
+    p0 = 8.15699e+01;
+    p1 = 1.48060e+01;
+    p2 = 1.50018e+00;
+    p3 = 6.34208e+00;
+
+    return (p0 * x / TMath::Power(1. + TMath::Power(x / p1, p2), p3));
+  }
+
+  //-------------------------------------------------------------------------//
+  static Double_t YUpsilon2Spp13TeV(const Double_t* py, const Double_t* /*dummy*/)
+  {
+    // Upsilon2s y shape from LHCb pp@13TeV arXiv:1804.09214
+    Double_t x = *py;
+    Float_t p0, p1;
+
+    p0 =  7.50409e+02;
+    p1 = -3.57039e-02;
+
+    return (p0 * (1. + p1 * x * x));
+  }
+
+  //-------------------------------------------------------------------------//
+  static Double_t V2Upsilon2Spp13TeV(const Double_t* /*dummy*/, const Double_t* /*dummy*/)
+  {
+    // Upsilon(2S) v2
+    return 0.;
+  }
+
+  //-------------------------------------------------------------------------//
+  static Int_t IpUpsilon2Spp13TeV(TRandom*)
+  {
+    return 100553;
+  }
+
+ private:
+  GeneratorParam* paramUpsilon2S = nullptr;
+};
+
+/////////////////////////////////////////////////////////////////////////////
+class O2_GeneratorParamUpsilon3SFwdY : public GeneratorTGenerator
+{
+
+ public:
+  O2_GeneratorParamUpsilon3SFwdY() : GeneratorTGenerator("ParamUpsilon3S")
+  {
+    paramUpsilon3S = new GeneratorParam(1, -1, PtUpsilon3Spp13TeV, YUpsilon3Spp13TeV, V2Upsilon3Spp13TeV, IpUpsilon3Spp13TeV);
+    paramUpsilon3S->SetMomentumRange(0., 1.e6);
+    paramUpsilon3S->SetPtRange(0, 999.);
+    paramUpsilon3S->SetYRange(-4.2, -2.3);
+    paramUpsilon3S->SetPhiRange(0., 360.);
+    paramUpsilon3S->SetDecayer(new TPythia6Decayer());
+    paramUpsilon3S->SetForceDecay(kNoDecay); // particle left undecayed
+    // - - paramUpsilon3S->SetTrackingFlag(1);  // (from AliGenParam) -> check this
+    setTGenerator(paramUpsilon3S);
+  };
+
+  ~O2_GeneratorParamUpsilon3SFwdY()
+  {
+    delete paramUpsilon3S;
+  };
+
+  Bool_t Init() override
+  {
+    GeneratorTGenerator::Init();
+    paramUpsilon3S->Init();
+    return true;
+  }
+
+  void SetNSignalPerEvent(Int_t nsig) { paramUpsilon3S->SetNumberParticles(nsig); }
+
+  //-------------------------------------------------------------------------//
+  static Double_t PtUpsilon3Spp13TeV(const Double_t* px, const Double_t* /*dummy*/)
+  {
+    // Upsilon3S pt shape from LHCb pp@13TeV arXiv:1804.09214
+    Double_t x = *px;
+    Float_t p0, p1, p2, p3;
+
+    p0 = 3.51590e+01;
+    p1 = 2.30813e+01;
+    p2 = 1.40822e+00;
+    p3 = 9.38026e+00;
+
+    return (p0 * x / TMath::Power(1. + TMath::Power(x / p1, p2), p3));
+  }
+
+  //-------------------------------------------------------------------------//
+  static Double_t YUpsilon3Spp13TeV(const Double_t* py, const Double_t* /*dummy*/)
+  {
+    // Upsilon3s y shape from LHCb pp@13TeV arXiv:1804.09214
+    Double_t x = *py;
+    Float_t p0, p1;
+
+    p0 =  3.69961e+02;
+    p1 = -3.54650e-02;
+
+    return (p0 * (1. + p1 * x * x));
+  }
+
+  //-------------------------------------------------------------------------//
+  static Double_t V2Upsilon3Spp13TeV(const Double_t* /*dummy*/, const Double_t* /*dummy*/)
+  {
+    // Upsilon(3S) v2
+    return 0.;
+  }
+
+  //-------------------------------------------------------------------------//
+  static Int_t IpUpsilon3Spp13TeV(TRandom*)
+  {
+    return 200553;
+  }
+
+ private:
+  GeneratorParam* paramUpsilon3S = nullptr;
+};
+
+
+} // namespace eventgen
+} // namespace o2
+
+FairGenerator* GeneratorCocktailBottomoniaToMuonEvtGen_pp13TeV()
+{
+
+  auto genCocktailEvtGen = new o2::eventgen::GeneratorEvtGen<GeneratorCocktail>();
+
+  auto genUpsilon1S = new o2::eventgen::O2_GeneratorParamUpsilon1SFwdY;
+  genUpsilon1S->SetNSignalPerEvent(1); // 1 Upsilon(1S) generated per event by GeneratorParam
+
+  auto genUpsilon2S = new o2::eventgen::O2_GeneratorParamUpsilon2SFwdY;
+  genUpsilon2S->SetNSignalPerEvent(1); // 1 Upsilon(2S) generated per event by GeneratorParam
+
+  auto genUpsilon3S = new o2::eventgen::O2_GeneratorParamUpsilon3SFwdY;
+  genUpsilon3S->SetNSignalPerEvent(1); // 1 Upsilon(3S) generated per event by GeneratorParam
+
+  genCocktailEvtGen->AddGenerator(genUpsilon1S, 1); // add Upsilon(1S) generator
+  genCocktailEvtGen->AddGenerator(genUpsilon2S, 1); // add Upsilon(2S) generator
+  genCocktailEvtGen->AddGenerator(genUpsilon3S, 1); // add Upsilon(3S) generator
+
+  TString pdgs = "553;100553;200553";
+  std::string spdg;
+  TObjArray* obj = pdgs.Tokenize(";");
+  genCocktailEvtGen->SetSizePdg(obj->GetEntriesFast());
+  for (int i = 0; i < obj->GetEntriesFast(); i++) {
+    spdg = obj->At(i)->GetName();
+    genCocktailEvtGen->AddPdg(std::stoi(spdg), i);
+    printf("PDG %d \n", std::stoi(spdg));
+  }
+  genCocktailEvtGen->SetForceDecay(kEvtDiMuon);
+
+  return genCocktailEvtGen;
+}

MC/config/PWGDQ/external/generator/generator_pythia8_withInjectedBottomoniaSignals_gaptriggered_dq.C

@@ -0,0 +1,76 @@
+#include "FairGenerator.h"
+#include "Generators/GeneratorPythia8.h"
+#include "Pythia8/Pythia.h"
+#include "TRandom.h"
+#include "GeneratorBottomonia.C"
+#include <string>
+
+using namespace o2::eventgen;
+using namespace Pythia8;
+
+class GeneratorPythia8BottomoniaInjectedGapTriggeredDQ : public o2::eventgen::GeneratorPythia8 {
+public:
+
+  /// default constructor
+  GeneratorPythia8BottomoniaInjectedGapTriggeredDQ() = default;
+
+  /// constructor
+  GeneratorPythia8BottomoniaInjectedGapTriggeredDQ(int inputTriggerRatio = 5, int gentype = 0) {
+
+    mGeneratedEvents = 0;
+    mGeneratorParam = 0x0;
+    mInverseTriggerRatio = inputTriggerRatio; 
+    switch (gentype) {
+      case 0: // generate bottomonia cocktail at forward rapidity
+        mGeneratorParam = (Generator*)GeneratorCocktailBottomoniaToMuonEvtGen_pp13TeV(); 
+        break;
+    }
+    mGeneratorParam->Init();  
+  }
+
+  ///  Destructor
+  ~GeneratorPythia8BottomoniaInjectedGapTriggeredDQ() = default;
+
+protected:
+
+Bool_t importParticles() override
+  {
+    GeneratorPythia8::importParticles();
+    bool genOk = false;
+    if (mGeneratedEvents % mInverseTriggerRatio == 0) { // add injected prompt signals to the stack
+      bool genOk = false;
+      while (!genOk) {    
+        genOk =  (mGeneratorParam->generateEvent() && mGeneratorParam->importParticles()) ? true : false ;
+      }
+      int originalSize = mParticles.size();
+      for (int ipart=0; ipart < mGeneratorParam->getParticles().size(); ipart++) {
+        TParticle part = TParticle(mGeneratorParam->getParticles().at(ipart));
+        if(part.GetFirstMother() >= 0) part.SetFirstMother(part.GetFirstMother() + originalSize);
+        if(part.GetFirstDaughter() >= 0) part.SetFirstDaughter(part.GetFirstDaughter() + originalSize);
+        if(part.GetLastDaughter() >= 0) part.SetLastDaughter(part.GetLastDaughter() + originalSize);
+        mParticles.push_back(part); 
+        // encodeParticleStatusAndTracking method already called in GeneratorEvtGen.C 
+      }	   
+      mGeneratorParam->clearParticles(); 
+    }
+
+    mGeneratedEvents++;
+    return true;
+  }
+
+
+private:
+  Generator* mGeneratorParam; 
+  // Control gap-triggering
+  unsigned long long mGeneratedEvents;
+  int mInverseTriggerRatio;
+};
+
+// Predefined generators:
+FairGenerator *GeneratorPythia8InjectedBottomoniaGapTriggered(int inputTriggerRatio, int gentype) {
+  auto myGen = new GeneratorPythia8BottomoniaInjectedGapTriggeredDQ(inputTriggerRatio,gentype);
+  auto seed = (gRandom->TRandom::GetSeed() % 900000000);
+  myGen->readString("Random:setSeed on");
+  myGen->readString("Random:seed " + std::to_string(seed));
+  return myGen;
+}

MC/config/PWGDQ/ini/Generator_InjectedBottomoniaFwdy_TriggerGap.ini

@@ -0,0 +1,7 @@
+### The external generator derives from GeneratorPythia8.
+[GeneratorExternal]
+fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/external/generator/generator_pythia8_withInjectedBottomoniaSignals_gaptriggered_dq.C
+funcName=GeneratorPythia8InjectedBottomoniaGapTriggered(0,3)
+
+[GeneratorPythia8]
+config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/pythia8/generator/pythia8_inel_triggerGap.cfg