[PWGDQ] Adding scripts to run DQ simulations in Pb-Pb collisions

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

@@ -492,6 +492,152 @@ class O2_GeneratorParamChiC2 : public GeneratorTGenerator
   GeneratorParam* paramChiC2 = nullptr;
 };
 
+class O2_GeneratorParamJpsiPbPb5TeV : public GeneratorTGenerator
+{
+
+ public:
+  O2_GeneratorParamJpsiPbPb5TeV() : GeneratorTGenerator("ParamJpsi")
+  {
+    paramJpsi = new GeneratorParam(1, -1, PtJPsiPbPb5TeV, YJPsiPbPb5TeV, V2JPsiPbPb5TeV, IpJPsiPbPb5TeV);
+    paramJpsi->SetMomentumRange(0., 1.e6);
+    paramJpsi->SetPtRange(0, 999.);
+    paramJpsi->SetYRange(-4.2, -2.3);
+    paramJpsi->SetPhiRange(0., 360.);
+    paramJpsi->SetDecayer(new TPythia6Decayer());
+    paramJpsi->SetForceDecay(kNoDecay); // particle left undecayed
+    // - - paramJpsi->SetTrackingFlag(1);  // (from AliGenParam) -> check this
+    setTGenerator(paramJpsi);
+  };
+
+  ~O2_GeneratorParamJpsiPbPb5TeV()
+  {
+    delete paramJpsi;
+  };
+
+  Bool_t Init() override
+  {
+    GeneratorTGenerator::Init();
+    paramJpsi->Init();
+    return true;
+  }
+
+  void SetNSignalPerEvent(Int_t nsig) { paramJpsi->SetNumberParticles(nsig); }
+
+  //-------------------------------------------------------------------------//
+  static Double_t PtJPsiPbPb5TeV(const Double_t* px, const Double_t* /*dummy*/)
+  {
+    // jpsi pT in PbPb, tuned on data (2015) -> Castillo embedding https://alice.its.cern.ch/jira/browse/ALIROOT-8174?jql=text%20~%20%22LHC19a2%22
+    Double_t x = *px;
+    Float_t p0, p1, p2, p3;
+    p0 = 1.00715e6;
+    p1 = 3.50274;
+    p2 = 1.93403;
+    p3 = 3.96363;
+    return p0 * x / TMath::Power(1. + TMath::Power(x / p1, p2), p3);
+  }
+
+  //-------------------------------------------------------------------------//
+  static Double_t YJPsiPbPb5TeV(const Double_t* py, const Double_t* /*dummy*/)
+  {
+    // jpsi y in PbPb, tuned on data (2015) -> Castillo embedding https://alice.its.cern.ch/jira/browse/ALIROOT-8174?jql=text%20~%20%22LHC19a2%22
+    Double_t y = *py;
+    Float_t p0, p1, p2;
+    p0 = 1.09886e6;
+    p1 = 0;
+    p2 = 2.12568;
+    return p0 * TMath::Exp(-(1. / 2.) * TMath::Power(((y - p1) / p2), 2));
+  }
+
+  //-------------------------------------------------------------------------//
+  static Double_t V2JPsiPbPb5TeV(const Double_t* /*dummy*/, const Double_t* /*dummy*/)
+  {
+    // jpsi v2
+    return 0.;
+  }
+
+  //-------------------------------------------------------------------------//
+  static Int_t IpJPsiPbPb5TeV(TRandom*)
+  {
+    return 443;
+  }
+
+ private:
+  GeneratorParam* paramJpsi = nullptr;
+};
+
+class O2_GeneratorParamPsiPbPb5TeV : public GeneratorTGenerator
+{
+
+ public:
+  O2_GeneratorParamPsiPbPb5TeV() : GeneratorTGenerator("ParamPsi")
+  {
+    paramPsi = new GeneratorParam(1, -1, PtPsiPbPb5TeV, YPsiPbPb5TeV, V2PsiPbPb5TeV, IpPsiPbPb5TeV);
+    paramPsi->SetMomentumRange(0., 1.e6);
+    paramPsi->SetPtRange(0, 999.);
+    paramPsi->SetYRange(-4.2, -2.3);
+    paramPsi->SetPhiRange(0., 360.);
+    paramPsi->SetDecayer(new TPythia6Decayer());
+    paramPsi->SetForceDecay(kNoDecay); // particle left undecayed
+    // - - paramJpsi->SetTrackingFlag(1);  // check this
+    setTGenerator(paramPsi);
+  };
+
+  ~O2_GeneratorParamPsiPbPb5TeV()
+  {
+    delete paramPsi;
+  };
+
+  Bool_t Init() override
+  {
+    GeneratorTGenerator::Init();
+    paramPsi->Init();
+    return true;
+  }
+
+  void SetNSignalPerEvent(Int_t nsig) { paramPsi->SetNumberParticles(nsig); }
+
+  //-------------------------------------------------------------------------//
+  static Double_t PtPsiPbPb5TeV(const Double_t* px, const Double_t* /*dummy*/)
+  {
+    // jpsi pT in PbPb, tuned on data (2015) -> Castillo embedding https://alice.its.cern.ch/jira/browse/ALIROOT-8174?jql=text%20~%20%22LHC19a2%22
+    Double_t x = *px;
+    Float_t p0, p1, p2, p3;
+    p0 = 1.00715e6;
+    p1 = 3.50274;
+    p2 = 1.93403;
+    p3 = 3.96363;
+    return p0 * x / TMath::Power(1. + TMath::Power(x / p1, p2), p3);
+  }
+
+  //-------------------------------------------------------------------------//
+  static Double_t YPsiPbPb5TeV(const Double_t* py, const Double_t* /*dummy*/)
+  {
+    // jpsi y in PbPb, tuned on data (2015) -> Castillo embedding https://alice.its.cern.ch/jira/browse/ALIROOT-8174?jql=text%20~%20%22LHC19a2%22
+    Double_t y = *py;
+    Float_t p0, p1, p2;
+    p0 = 1.09886e6;
+    p1 = 0;
+    p2 = 2.12568;
+    return p0 * TMath::Exp(-(1. / 2.) * TMath::Power(((y - p1) / p2), 2));
+  }
+
+  //-------------------------------------------------------------------------//
+  static Double_t V2PsiPbPb5TeV(const Double_t* /*dummy*/, const Double_t* /*dummy*/)
+  {
+    // jpsi v2
+    return 0.;
+  }
+
+  //-------------------------------------------------------------------------//
+  static Int_t IpPsiPbPb5TeV(TRandom*)
+  {
+    return 100443;
+  }
+
+ private:
+  GeneratorParam* paramPsi = nullptr;
+};
+
 } // namespace eventgen
 } // namespace o2
 
@@ -671,6 +817,32 @@ FairGenerator*
   return genCocktailEvtGen;
 }
 
+FairGenerator* 
+  GeneratorCocktailPromptCharmoniaToMuonEvtGen_PbPb5TeV()
+{
+  auto genCocktailEvtGen = new o2::eventgen::GeneratorEvtGen<GeneratorCocktail>();
+
+  auto genJpsi = new o2::eventgen::O2_GeneratorParamJpsiPbPb5TeV;
+  genJpsi->SetNSignalPerEvent(4); // 4 J/psi generated per event by GeneratorParam
+  auto genPsi = new o2::eventgen::O2_GeneratorParamPsiPbPb5TeV;
+  genPsi->SetNSignalPerEvent(2);               // 2 Psi(2S) generated per event by GeneratorParam
+  genCocktailEvtGen->AddGenerator(genJpsi, 1); // 2/3 J/psi
+  genCocktailEvtGen->AddGenerator(genPsi, 1);  // 1/3 Psi(2S)
+
+  TString pdgs = "443;100443";
+  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_withInjectedPromptSignals_gaptriggered_dq.C

@@ -42,6 +42,9 @@ public:
       case 6: // generate prompt ChiC1 + ChiC2 cocktail at midrapidity
         mGeneratorParam = (Generator*)GeneratorCocktailChiCToElectronEvtGen_pp13TeV(); 
         break;
+      case 7: // generate prompt charmonia cocktail at forward rapidity at 5TeV
+        mGeneratorParam = (Generator*)GeneratorCocktailPromptCharmoniaToMuonEvtGen_PbPb5TeV(); 
+        break;
 
       }
     mGeneratorParam->Init();  

MC/config/PWGDQ/ini/Generator_InjectedPromptCharmoniaFwdy_TriggerGap_PbPb5TeV.ini

@@ -0,0 +1,7 @@
+### The external generator derives from GeneratorPythia8.
+[GeneratorExternal]
+fileName=${O2DPG_ROOT}/MC/config/PWGDQ/external/generator/generator_pythia8_withInjectedPromptSignals_gaptriggered_dq.C
+funcName=GeneratorPythia8InjectedPromptCharmoniaGapTriggered(5,7)
+
+[GeneratorPythia8]
+config=${O2DPG_ROOT}/MC/config/PWGDQ/pythia8/generator/pythia8_inel_triggerGap_5TeV.cfg

MC/config/PWGDQ/ini/tests/Generator_InjectedPromptCharmoniaFwdy_TriggerGap_PbPb5TeV.C

@@ -0,0 +1,85 @@
+int External()
+{
+    int checkPdgSignal[] = {443,100443};
+    int checkPdgDecay = 13;
+    double rapiditymin = -4.3; double rapiditymax = -2.3;
+    std::string path{"o2sim_Kine.root"};
+    std::cout << "Check for\nsignal PDG " << checkPdgSignal << "\ndecay PDG " << checkPdgDecay << "\n";
+    TFile file(path.c_str(), "READ");
+    if (file.IsZombie()) {
+        std::cerr << "Cannot open ROOT file " << path << "\n";
+        return 1;
+    }
+
+    auto tree = (TTree*)file.Get("o2sim");
+    std::vector<o2::MCTrack>* tracks{};
+    tree->SetBranchAddress("MCTrack", &tracks);
+
+    int nLeptons{};
+    int nAntileptons{};
+    int nLeptonPairs{};
+    int nLeptonPairsToBeDone{};
+    int nSignalJpsi{};
+    int nSignalPsi2S{};
+    int nSignalJpsiWithinAcc{};
+    int nSignalPsi2SWithinAcc{};
+    auto nEvents = tree->GetEntries();
+    o2::steer::MCKinematicsReader mcreader("o2sim", o2::steer::MCKinematicsReader::Mode::kMCKine);
+    Bool_t isInjected = kFALSE;
+
+    for (int i = 0; i < nEvents; i++) {
+        tree->GetEntry(i);
+        for (auto& track : *tracks) {
+            auto pdg = track.GetPdgCode();
+	    auto rapidity =  track.GetRapidity();
+	    auto idMoth = track.getMotherTrackId();
+            if (pdg == checkPdgDecay) {
+                // count leptons
+                nLeptons++;
+            } else if(pdg == -checkPdgDecay) {
+                // count anti-leptons
+                nAntileptons++;
+            } else if (pdg == checkPdgSignal[0] || pdg == checkPdgSignal[1]) {
+                if(idMoth < 0){
+		// count signal PDG 
+		pdg == checkPdgSignal[0] ? nSignalJpsi++ : nSignalPsi2S++;
+                // count signal PDG within acceptance 
+		if(rapidity > rapiditymin && rapidity < rapiditymax) { pdg == checkPdgSignal[0] ? nSignalJpsiWithinAcc++ : nSignalPsi2SWithinAcc++;}
+		}
+		auto child0 = o2::mcutils::MCTrackNavigator::getDaughter0(track, *tracks);
+                auto child1 = o2::mcutils::MCTrackNavigator::getDaughter1(track, *tracks);
+                if (child0 != nullptr && child1 != nullptr) {
+                    // check for parent-child relations
+                    auto pdg0 = child0->GetPdgCode();
+                    auto pdg1 = child1->GetPdgCode();
+                    std::cout << "First and last children of parent " << checkPdgSignal << " are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n";
+                    if (std::abs(pdg0) == checkPdgDecay && std::abs(pdg1) == checkPdgDecay && pdg0 == -pdg1) {
+                        nLeptonPairs++;
+                        if (child0->getToBeDone() && child1->getToBeDone()) {
+                            nLeptonPairsToBeDone++;
+                        }
+                    }
+                }
+            }
+        }
+    }
+    std::cout << "#events: " << nEvents << "\n"
+              << "#leptons: " << nLeptons << "\n"
+              << "#antileptons: " << nAntileptons << "\n"
+              << "#signal (prompt Jpsi): " << nSignalJpsi << "; within acceptance " << rapiditymin << " < y < " << rapiditymax << " : " << nSignalJpsiWithinAcc << "\n"
+              << "#signal (prompt Psi(2S)): " << nSignalPsi2S << "; within acceptance " << rapiditymin << " < y < " << rapiditymax << " : " << nSignalPsi2SWithinAcc << "\n"
+              << "#lepton pairs: " << nLeptonPairs << "\n"
+              << "#lepton pairs to be done: " << nLeptonPairs << "\n";
+
+
+    if (nLeptonPairs == 0 || nLeptons == 0 || nAntileptons == 0) {
+        std::cerr << "Number of leptons, number of anti-leptons as well as number of lepton pairs should all be greater than 1.\n";
+        return 1;
+    }
+    if (nLeptonPairs != nLeptonPairsToBeDone) {
+        std::cerr << "The number of lepton pairs should be the same as the number of lepton pairs which should be transported.\n";
+        return 1;
+    }
+
+    return 0;
+}

MC/config/PWGDQ/pythia8/generator/pythia8_inel_triggerGap_5TeV.cfg

@@ -0,0 +1,11 @@
+### beams
+Beams:idA 2212			# proton
+Beams:idB 2212 			# proton
+Beams:eCM 5360. 		# GeV
+
+### processes
+SoftQCD:inelastic on		# all inelastic processes
+
+### decays
+ParticleDecays:limitTau0 on	
+ParticleDecays:tau0Max 10.