add X3872 and Psi2S to Jpsi pi pi

MC/config/PWGDQ/EvtGen/DecayTablesEvtgen/PSITOJPSIPIPI.DEC

@@ -0,0 +1,8 @@
+Decay J/psi
+### from DECAY.DEC
+1.000 e+   e-                PHOTOS VLL;
+Enddecay
+Decay psi(2S)
+1.000 J/psi   pi+     pi-         VVPIPI;
+Enddecay
+End

MC/config/PWGDQ/EvtGen/DecayTablesEvtgen/X3872TOJPSIPIPI.DEC

@@ -0,0 +1,10 @@
+Decay X_1(3872)
+#### Breit-Wigner function for the pion distribution (S-Wave approximation)
+1.000 J/psi   pi+     pi-         XJPSIRHO0PIPI_S;
+Enddecay
+Decay J/psi
+### from DECAY.DEC
+1.000 	e+   	e-                PHOTOS VLL;
+Enddecay
+
+End

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

@@ -638,6 +638,75 @@ class O2_GeneratorParamPsiPbPb5TeV : public GeneratorTGenerator
   GeneratorParam* paramPsi = nullptr;
 };
 
+
+class O2_GeneratorParamX3872MidY : public GeneratorTGenerator
+{
+
+ public:
+  O2_GeneratorParamX3872MidY() : GeneratorTGenerator("ParamX3872MidY")
+  {
+    paramX3872 = new GeneratorParam(1, -1, PtX3872pp13TeV, YX3872pp13TeV, V2X3872pp13TeV, IpX3872pp13TeV);
+    paramX3872->SetMomentumRange(0., 1.e6);
+    paramX3872->SetPtRange(0., 1000.);
+    paramX3872->SetYRange(-1.0, 1.0);
+    paramX3872->SetPhiRange(0., 360.);
+    paramX3872->SetDecayer(new TPythia6Decayer()); // Pythia
+    paramX3872->SetForceDecay(kNoDecay);           // particle left undecayed
+    setTGenerator(paramX3872);
+  };
+
+  ~O2_GeneratorParamX3872MidY()
+  {
+    delete paramX3872;
+  };
+
+  Bool_t Init() override
+  {
+    GeneratorTGenerator::Init();
+    paramX3872->Init();
+    return true;
+  }
+
+  void SetNSignalPerEvent(Int_t nsig) { paramX3872->SetNumberParticles(nsig); }
+
+  //-------------------------------------------------------------------------//
+  static Double_t PtX3872pp13TeV(const Double_t* px, const Double_t* /*dummy*/)
+  {
+    // prompt X3872 pT
+    // pp, 13TeV (tuned LHCb pp 13 TeV)
+    //
+    const Double_t kC = 7.64519e+00 ;
+    const Double_t kpt0 = 5.30628e+00;
+    const Double_t kn = 3.30887e+00;
+    Double_t pt = px[0];
+    return kC * pt / TMath::Power((1. + (pt / kpt0) * (pt / kpt0)), kn);
+  }
+
+  //-------------------------------------------------------------------------//
+  static Double_t YX3872pp13TeV(const Double_t* py, const Double_t* /*dummy*/)
+  {
+    // flat rapidity distribution assumed at midrapidity
+    return 1.;
+  }
+
+  //-------------------------------------------------------------------------//
+  static Double_t V2X3872pp13TeV(const Double_t* /*dummy*/, const Double_t* /*dummy*/)
+  {
+    // X3872 v2
+    return 0.;
+  }
+
+  //-------------------------------------------------------------------------//
+  static Int_t IpX3872pp13TeV(TRandom*)
+  {
+    return 9920443;
+  }
+
+ private:
+  GeneratorParam* paramX3872 = nullptr;
+};
+
+
 } // namespace eventgen
 } // namespace o2
 
@@ -741,6 +810,31 @@ FairGenerator* GeneratorCocktailPromptCharmoniaToMuonEvtGen_pp13TeV()
   return genCocktailEvtGen;
 }
 
+
+FairGenerator*
+  GeneratorParamPromptPsiToJpsiPiPiEvtGen_pp13TeV(TString pdgs = "100443")
+{
+  auto gen = new o2::eventgen::GeneratorEvtGen<o2::eventgen::O2_GeneratorParamPsiMidY>();
+  gen->SetNSignalPerEvent(1); // number of jpsis per event
+
+  std::string spdg;
+  TObjArray* obj = pdgs.Tokenize(";");
+  gen->SetSizePdg(obj->GetEntriesFast());
+  for (int i = 0; i < obj->GetEntriesFast(); i++) {
+    spdg = obj->At(i)->GetName();
+    gen->AddPdg(std::stoi(spdg), i);
+    printf("PDG %d \n", std::stoi(spdg));
+  }
+  TString pathO2 = gSystem->ExpandPathName("$O2DPG_ROOT/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen");
+  gen->SetDecayTable(Form("%s/PSITOJPSIPIPI.DEC", pathO2.Data()));
+
+  // print debug
+  gen->PrintDebug();
+
+  return gen;
+}
+
+
 FairGenerator*
   GeneratorParamPromptJpsiToMuonEvtGen_pp13TeV(TString pdgs = "443")
 {
@@ -844,6 +938,29 @@ FairGenerator*
 }
 
 
+FairGenerator*
+  GeneratorParamX3872ToJpsiEvtGen_pp13TeV(TString pdgs = "9920443")
+{
+  auto gen = new o2::eventgen::GeneratorEvtGen<o2::eventgen::O2_GeneratorParamX3872MidY>();
+  gen->SetNSignalPerEvent(1); // number of jpsis per event
+
+  std::string spdg;
+  TObjArray* obj = pdgs.Tokenize(";");
+  gen->SetSizePdg(obj->GetEntriesFast());
+  for (int i = 0; i < obj->GetEntriesFast(); i++) {
+    spdg = obj->At(i)->GetName();
+    gen->AddPdg(std::stoi(spdg), i);
+    printf("PDG %d \n", std::stoi(spdg));
+  }
+  TString pathO2 = gSystem->ExpandPathName("$O2DPG_ROOT/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen");
+  gen->SetDecayTable(Form("%s/X3872TOJPSIPIPI.DEC", pathO2.Data()));
+
+  // print debug
+  gen->PrintDebug();
+
+  return gen;
+}
+
 
 
 

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

@@ -45,7 +45,12 @@ public:
       case 7: // generate prompt charmonia cocktail at forward rapidity at 5TeV
         mGeneratorParam = (Generator*)GeneratorCocktailPromptCharmoniaToMuonEvtGen_PbPb5TeV(); 
         break;
-
+      case 8: // generate prompt X_1(3872) to Jpsi pi pi at midrapidity
+        mGeneratorParam = (Generator*)GeneratorParamX3872ToJpsiEvtGen_pp13TeV("9920443");
+        break;
+      case 9: // generate prompt psi2S to Jpsi pi pi at midrapidity
+        mGeneratorParam = (Generator*)GeneratorParamPromptPsiToJpsiPiPiEvtGen_pp13TeV("100443");
+        break;
       }
     mGeneratorParam->Init();  
   }

MC/config/PWGDQ/ini/Generator_InjectedPromptPsi2SToJpsiPiPiMidy_TriggerGap.ini

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

MC/config/PWGDQ/ini/Generator_InjectedPromptX3782Midy_TriggerGap.ini

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

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

@@ -0,0 +1,99 @@
+int External()
+{
+    int checkPdgSignal[] = {100443};
+    int checkPdgDecay[] = {443, 211, -211}; 
+    int leptonPdg = 11;
+    Double_t rapidityWindow = 1.0;
+    std::string path{"o2sim_Kine.root"};
+    std::cout << "Check for\nsignal PDG " << checkPdgSignal[0] << "\n decay PDG " << checkPdgDecay[0] << ", " << checkPdgDecay[1] << ", " << checkPdgDecay[2] << "\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 nSignalPionsPos{};
+    int nSignalPionsNeg{};
+    int nSignalPsi2S{};
+    int nSignalJpsiWithinAcc{};
+    int nSignalPionsPosWithinAcc{};
+    int nSignalPionsNegWithinAcc{};
+    auto nEvents = tree->GetEntries();
+    o2::steer::MCKinematicsReader mcreader("o2sim", o2::steer::MCKinematicsReader::Mode::kMCKine);
+    Bool_t hasPsi2SMoth = 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();
+	    int idJpsi = -1; int IdChild0 = -1; int IdChild1 = -1;
+            if (pdg == leptonPdg) {
+                // count leptons
+                nLeptons++;
+            } else if(pdg == -leptonPdg) {
+                // count anti-leptons
+                nAntileptons++;
+            } else if (pdg == checkPdgSignal[0]) {
+		// check daughters
+                  std::cout << "Signal PDG: " << pdg << "\n";
+                  for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j) {
+                  auto pdgDau = tracks->at(j).GetPdgCode();
+                  std::cout << "Daughter " << j << " is: " << pdgDau << "\n";
+                    if(TMath::Abs(pdgDau) == checkPdgDecay[0] ) { nSignalJpsi++; if( std::abs(track.GetRapidity()) < rapidityWindow) nSignalJpsiWithinAcc++; idJpsi = j; }
+                    if(pdgDau == checkPdgDecay[1] ) { nSignalPionsPos++; if( std::abs(track.GetRapidity()) < rapidityWindow) nSignalPionsPosWithinAcc++; }
+                    if(pdgDau == checkPdgDecay[2] ) { nSignalPionsNeg++; if( std::abs(track.GetRapidity()) < rapidityWindow) nSignalPionsNegWithinAcc++; }
+		  }
+
+		  auto trackJpsi = tracks->at(idJpsi);
+		  for (int j{trackJpsi.getFirstDaughterTrackId()}; j <= trackJpsi.getLastDaughterTrackId(); ++j) {
+                  auto pdgDau = tracks->at(j).GetPdgCode();
+		   if(pdgDau == leptonPdg ) IdChild0 = j;
+		   if(pdgDau == -leptonPdg ) IdChild1 = j;
+		  }
+         	auto child0 = tracks->at(IdChild0);
+         	auto child1 = tracks->at(IdChild1);
+                    // check for parent-child relations
+                    auto pdg0 = child0.GetPdgCode();
+                    auto pdg1 = child1.GetPdgCode();
+                    std::cout << "Lepton daughter particles of mother  " << trackJpsi.GetPdgCode() << " are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n"; 
+                    if (std::abs(pdg0) == leptonPdg && std::abs(pdg1) == leptonPdg && pdg0 == -pdg1) {
+                        nLeptonPairs++;
+                        if (child0.getToBeDone() && child1.getToBeDone()) {
+                            nLeptonPairsToBeDone++;
+                        }
+                    }
+                //} 
+            }
+        }
+    }
+    std::cout << "#events: " << nEvents << "\n"
+              << "#leptons: " << nLeptons << "\n"
+              << "#antileptons: " << nAntileptons << "\n"
+              << "#signal (jpsi <- psi2S): " << nSignalJpsi << "; within acceptance (|y| < " << rapidityWindow << "): " << nSignalJpsiWithinAcc << "\n"
+              << "#signal (pi+ <- psi2S): " << nSignalPionsPos << "; within acceptance (|y| < " << rapidityWindow << "): " << nSignalPionsPosWithinAcc << "\n"
+              << "#signal (pi- <- psi2S): " << nSignalPionsNeg << "; within acceptance (|y| < " << rapidityWindow << "): " << nSignalPionsNegWithinAcc << "\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/ini/tests/Generator_InjectedPromptX3782Midy_TriggerGap.C

@@ -0,0 +1,99 @@
+int External()
+{
+    int checkPdgSignal[] = {9920443}; // pdg code X3872
+    int checkPdgDecay[] = {443, 211, -211}; 
+    int leptonPdg = 11;
+    Double_t rapidityWindow = 1.0;
+    std::string path{"o2sim_Kine.root"};
+    std::cout << "Check for\nsignal PDG " << checkPdgSignal[0] << "\n decay PDG " << checkPdgDecay[0] << ", " << checkPdgDecay[1] << ", " << checkPdgDecay[2] << "\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 nSignalX3872{};
+    int nSignalPionsPos{};
+    int nSignalPionsNeg{};
+    int nSignalPsi2S{};
+    int nSignalX3872WithinAcc{};
+    int nSignalPionsPosWithinAcc{};
+    int nSignalPionsNegWithinAcc{};
+    auto nEvents = tree->GetEntries();
+    o2::steer::MCKinematicsReader mcreader("o2sim", o2::steer::MCKinematicsReader::Mode::kMCKine);
+    Bool_t hasPsi2SMoth = 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();
+	    int idX3872 = -1; int IdChild0 = -1; int IdChild1 = -1;
+            if (pdg == leptonPdg) {
+                // count leptons
+                nLeptons++;
+            } else if(pdg == -leptonPdg) {
+                // count anti-leptons
+                nAntileptons++;
+            } else if (pdg == checkPdgSignal[0]) {
+		// check daughters
+                  std::cout << "Signal PDG: " << pdg << "\n";
+                  for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j) {
+                  auto pdgDau = tracks->at(j).GetPdgCode();
+                  std::cout << "Daughter " << j << " is: " << pdgDau << "\n";
+                    if(TMath::Abs(pdgDau) == checkPdgDecay[0] ) { nSignalX3872++; if( std::abs(track.GetRapidity()) < rapidityWindow) nSignalX3872WithinAcc++; idX3872 = j; }
+                    if(pdgDau == checkPdgDecay[1] ) { nSignalPionsPos++; if( std::abs(track.GetRapidity()) < rapidityWindow) nSignalPionsPosWithinAcc++; }
+                    if(pdgDau == checkPdgDecay[2] ) { nSignalPionsNeg++; if( std::abs(track.GetRapidity()) < rapidityWindow) nSignalPionsNegWithinAcc++; }
+		  }
+
+		  auto trackX3872 = tracks->at(idX3872);
+		  for (int j{trackX3872.getFirstDaughterTrackId()}; j <= trackX3872.getLastDaughterTrackId(); ++j) {
+                  auto pdgDau = tracks->at(j).GetPdgCode();
+		   if(pdgDau == leptonPdg ) IdChild0 = j;
+		   if(pdgDau == -leptonPdg ) IdChild1 = j;
+		  }
+         	auto child0 = tracks->at(IdChild0);
+         	auto child1 = tracks->at(IdChild1);
+                    // check for parent-child relations
+                    auto pdg0 = child0.GetPdgCode();
+                    auto pdg1 = child1.GetPdgCode();
+                    std::cout << "Lepton daughter particles of mother  " << trackX3872.GetPdgCode() << " are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n"; 
+                    if (std::abs(pdg0) == leptonPdg && std::abs(pdg1) == leptonPdg && pdg0 == -pdg1) {
+                        nLeptonPairs++;
+                        if (child0.getToBeDone() && child1.getToBeDone()) {
+                            nLeptonPairsToBeDone++;
+                        }
+                    }
+                //} 
+            }
+        }
+    }
+    std::cout << "#events: " << nEvents << "\n"
+              << "#leptons: " << nLeptons << "\n"
+              << "#antileptons: " << nAntileptons << "\n"
+              << "#signal (jpsi <- X3872): " << nSignalX3872 << "; within acceptance (|y| < " << rapidityWindow << "): " << nSignalX3872WithinAcc << "\n"
+              << "#signal (pi+ <- X3872): " << nSignalPionsPos << "; within acceptance (|y| < " << rapidityWindow << "): " << nSignalPionsPosWithinAcc << "\n"
+              << "#signal (pi- <- X3872): " << nSignalPionsNeg << "; within acceptance (|y| < " << rapidityWindow << "): " << nSignalPionsNegWithinAcc << "\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;
+}