AliceO2Group · preghenella · Aug 16, 2021 · Aug 4, 2021
diff --git a/examples/pythia8/pythia_nuclei.cfg b/examples/pythia8/pythia_nuclei.cfg
@@ -1,3 +1,7 @@
 #Config file to define the nuclei species that are not in vanilla pythia
 1000020030:all 3He 3Hebar 1 6 0 2.8094
 1000010030:all 3Tr 3Trbar 1 3 0 2.8089218
+2010010020:all 2CDeuteron 2CDeuteronbar 1 3 0 3.226
+# c-deuteron -> deuteron K- pi+
+2010010020:tau0=0.06000000000
+2010010020:addChannel = 1 .1 0 1000010020 -321 211
diff --git a/examples/scripts/common.py b/examples/scripts/common.py
@@ -2,6 +2,7 @@

 """
 Common header for AOD python scripts
+Author: Nicolo' Jacazio, nicolo.jacazio@cern.ch
 """

 import argparse

diff --git a/examples/scripts/createO2tables.py b/examples/scripts/createO2tables.py
@@ -1,7 +1,8 @@
 #! /usr/bin/env python3

 """
-Handler to run the DelphesO2 framework and to create O2 analysis tables
+Handler to run the DelphesO2 framework and to create O2 analysis tables.
+Author: Nicolo' Jacazio, nicolo.jacazio@cern.ch
 """

 import configparser

diff --git a/examples/scripts/default_configfile.ini b/examples/scripts/default_configfile.ini
@@ -74,23 +74,41 @@ generators = $DELPHESO2_ROOT/examples/pythia8/pythia8_XeXe.cfg
 custom_gen = INPUT_FILES /tmp/AnalysisResults_*.hepmc

 [GUN]
-custom_gen = rpythia8-gun --pdg 421 --px 1. --py 0. --pz 0. --xProd 0. --yProd 0. --zProd 0. --config $DELPHESO2_ROOT/examples/pythia8/decays/force_hadronic_D.cfg --decay
+custom_gen = rpythia8-gun --pdg 421 --px 1. --py 0. --pz 0. --xProd 0. --yProd 0. --zProd 0. --config $DELPHESO2_ROOT/examples/pythia8/decays/force_hadronic_D.cfg --decay --npart 100

 [BOX_pion]
 custom_gen = rpythia8-box --pdg 211 --etamin -2. --etamax 2. --phimin 0. --phimax 6.28 --pmin 0 --pmax 10 --xProd 0. --yProd 0. --zProd 0. --decay --npart 100

 [BOX_proton]
-custom_gen = rpythia8-box --pdg 2212 --etamin -2. --etamax 2. --phimin 0. --phimax 6.28 --pmin 0 --pmax 10 --xProd 0. --yProd 0. --zProd 0. --decay --npart 100
+custom_gen = rpythia8-box --pdg 2212 --etamin -2. --etamax 2. --phimin 0. --phimax 6.28 --pmin 0 --pmax 10 --xProd 0. --yProd 0. --zProd 0. --decay --npart 100 --config $DELPHESO2_ROOT/examples/pythia8/pythia_nuclei.cfg --background-config $O2_ROOT/share/Generators/egconfig/pythia8_hi.cfg
+
+[BOX_proton_nobkg]
+custom_gen = rpythia8-box --pdg 2212 --etamin -2. --etamax 2. --phimin 0. --phimax 6.28 --pmin 0 --pmax 10 --xProd 0. --yProd 0. --zProd 0. --decay --npart 100 --config $DELPHESO2_ROOT/examples/pythia8/pythia_nuclei.cfg

 [BOX_deuteron]
-custom_gen = rpythia8-box --pdg 1000010020 --etamin -2. --etamax 2. --phimin 0. --phimax 6.28 --pmin 0 --pmax 10 --xProd 0. --yProd 0. --zProd 0. --decay --npart 100
+custom_gen = rpythia8-box --pdg 1000010020 --etamin -2. --etamax 2. --phimin 0. --phimax 6.28 --pmin 0 --pmax 10 --xProd 0. --yProd 0. --zProd 0. --decay --npart 100 --config $DELPHESO2_ROOT/examples/pythia8/pythia_nuclei.cfg --background-config $O2_ROOT/share/Generators/egconfig/pythia8_hi.cfg
+
+[BOX_deuteron_nobkg]
+custom_gen = rpythia8-box --pdg 1000010020 --etamin -2. --etamax 2. --phimin 0. --phimax 6.28 --pmin 0 --pmax 10 --xProd 0. --yProd 0. --zProd 0. --decay --npart 100 --config $DELPHESO2_ROOT/examples/pythia8/pythia_nuclei.cfg

 [BOX_triton]
+custom_gen = rpythia8-box --pdg 1000010030 --etamin -2. --etamax 2. --phimin 0. --phimax 6.28 --pmin 0 --pmax 10 --xProd 0. --yProd 0. --zProd 0. --decay --npart 100 --config $DELPHESO2_ROOT/examples/pythia8/pythia_nuclei.cfg --background-config $O2_ROOT/share/Generators/egconfig/pythia8_hi.cfg
+
+[BOX_triton_nobkg]
 custom_gen = rpythia8-box --pdg 1000010030 --etamin -2. --etamax 2. --phimin 0. --phimax 6.28 --pmin 0 --pmax 10 --xProd 0. --yProd 0. --zProd 0. --decay --npart 100 --config $DELPHESO2_ROOT/examples/pythia8/pythia_nuclei.cfg

 [BOX_helium3]
+custom_gen = rpythia8-box --pdg 1000020030 --etamin -2. --etamax 2. --phimin 0. --phimax 6.28 --pmin 0 --pmax 10 --xProd 0. --yProd 0. --zProd 0. --decay --npart 100 --config $DELPHESO2_ROOT/examples/pythia8/pythia_nuclei.cfg --background-config $O2_ROOT/share/Generators/egconfig/pythia8_hi.cfg
+
+[BOX_helium3_nobkg]
 custom_gen = rpythia8-box --pdg 1000020030 --etamin -2. --etamax 2. --phimin 0. --phimax 6.28 --pmin 0 --pmax 10 --xProd 0. --yProd 0. --zProd 0. --decay --npart 100 --config $DELPHESO2_ROOT/examples/pythia8/pythia_nuclei.cfg

+[BOX_c_deuteron_nobkg]
+custom_gen = rpythia8-box --pdg 2010010020 --etamin -2. --etamax 2. --phimin 0. --phimax 6.28 --pmin 0 --pmax 10 --xProd 0. --yProd 0. --zProd 0. --decay --npart 100 --config /home/njacazio/alice/DelphesO2/examples/pythia8/pythia_nuclei.cfg
+
+[BOX_c_deuteron]
+custom_gen = rpythia8-box --pdg 2010010020 --etamin -2. --etamax 2. --phimin 0. --phimax 6.28 --pmin 0 --pmax 10 --xProd 0. --yProd 0. --zProd 0. --decay --npart 100 --config /home/njacazio/alice/DelphesO2/examples/pythia8/pythia_nuclei.cfg --background-config $O2_ROOT/share/Generators/egconfig/pythia8_hi.cfg
+
 [GUN_Lc_pKpi]
 custom_gen = rpythia8-gun --pdg 4122 --px 1. --py 0. --pz 0. --xProd 1. --yProd 0. --zProd 0. --config $O2DPG_ROOT/MC/config/PWGHF/pythia8/decayer/force_hadronic_D_forceLcChannel1.cfg --decay


diff --git a/src/TrackSmearer.cc b/src/TrackSmearer.cc
@@ -135,6 +135,9 @@ TrackSmearer::smearTrack(O2Track &o2track, int pid, float nch)
 {

  auto pt = o2track.getPt();
+  if (abs(pid) == 1000020030) {
+    pt *= 2.f;
+  }
  auto eta = o2track.getEta();
  auto lutEntry = getLUTEntry(pid, nch, 0., eta, pt);
  if (!lutEntry || !lutEntry->valid) return false;

diff --git a/src/lutWrite_aod.cc b/src/lutWrite_aod.cc
@@ -4,6 +4,8 @@
 /// @author: Nicolo' Jacazio
 /// @email: nicolo.jacazio@cern.ch

+/// \brief LUT writer with the AnalysisResults from AODs as input
+
 #include "TCanvas.h"
 #include "TDatabasePDG.h"
 #include "TFile.h"
@@ -12,32 +14,23 @@
 #include "TProfile2D.h"
 #include "TProfile3D.h"
 #include "TVectorD.h"
-#include "fwdRes/fwdRes.C"
+#include "DetectorK/DetectorK.h"
 #include "lutCovm.hh"
 #include <Riostream.h>
-
-void diagonalise(lutEntry_t& lutEntry);
-
-bool fwdSolve(float* covm, float pt = 0.1, float eta = 0.0,
-    float mass = 0.13957000)
-{
-  if (fwdRes(covm, pt, eta, mass) < 0)
-    return false;
-  return true;
-}
+#include "lutWrite.cc"

 void lutWrite_aod(const char* filename = "/tmp/lutCovm.pi.aod.dat",
-    int pdg = 211,
-    float field = 0.2, int layer = 0, int what = 0,
-    int efftype = 0,
-    const char* infilename = "/tmp/AnalysisResults_LUT.root",
-    float minPt = 0.f,
-    float maxPt = 80.f,
-    float minEta = -4.f,
-    float maxEta = 4.f)
+                  int pdg = 211,
+                  float field = 0.2, int layer = 0, int what = 0,
+                  int efftype = 0,
+                  const char* infilename = "/tmp/AnalysisResults_LUT.root",
+                  float minPt = 0.f,
+                  float maxPt = 80.f,
+                  float minEta = -4.f,
+                  float maxEta = 4.f)
 {

-  std::map<int, std::string> partname{ { 11, "electron" }, { 13, "muon" }, { 211, "pion" }, { 321, "kaon" }, { 2212, "proton" } };
+  std::map<int, std::string> partname{{11, "electron"}, {13, "muon"}, {211, "pion"}, {321, "kaon"}, {2212, "proton"}};
  const std::string dn = "alice3-lutmaker-" + partname[pdg];

  // Get the input from the analysis results
@@ -55,23 +48,23 @@ void lutWrite_aod(const char* filename = "/tmp/lutCovm.pi.aod.dat",
    return;
  }
  // d->ls();
-  std::map<std::string, TH1F*> h{ { "eta", nullptr }, { "pt", nullptr } };
-  std::map<std::string, TProfile2D*> m{ { "CovMat_cYY", nullptr },
-    { "CovMat_cZY", nullptr },
-    { "CovMat_cZZ", nullptr },
-    { "CovMat_cSnpY", nullptr },
-    { "CovMat_cSnpZ", nullptr },
-    { "CovMat_cSnpSnp", nullptr },
-    { "CovMat_cTglY", nullptr },
-    { "CovMat_cTglZ", nullptr },
-    { "CovMat_cTglSnp", nullptr },
-    { "CovMat_cTglTgl", nullptr },
-    { "CovMat_c1PtY", nullptr },
-    { "CovMat_c1PtZ", nullptr },
-    { "CovMat_c1PtSnp", nullptr },
-    { "CovMat_c1PtTgl", nullptr },
-    { "CovMat_c1Pt21Pt2", nullptr },
-    { "Efficiency", nullptr } };
+  std::map<std::string, TH1F*> h{{"eta", nullptr}, {"pt", nullptr}};
+  std::map<std::string, TProfile2D*> m{{"CovMat_cYY", nullptr},
+                                       {"CovMat_cZY", nullptr},
+                                       {"CovMat_cZZ", nullptr},
+                                       {"CovMat_cSnpY", nullptr},
+                                       {"CovMat_cSnpZ", nullptr},
+                                       {"CovMat_cSnpSnp", nullptr},
+                                       {"CovMat_cTglY", nullptr},
+                                       {"CovMat_cTglZ", nullptr},
+                                       {"CovMat_cTglSnp", nullptr},
+                                       {"CovMat_cTglTgl", nullptr},
+                                       {"CovMat_c1PtY", nullptr},
+                                       {"CovMat_c1PtZ", nullptr},
+                                       {"CovMat_c1PtSnp", nullptr},
+                                       {"CovMat_c1PtTgl", nullptr},
+                                       {"CovMat_c1Pt21Pt2", nullptr},
+                                       {"Efficiency", nullptr}};

  struct binning {
    int n = 0;
@@ -288,31 +281,3 @@ void lutWrite_aod(const char* filename = "/tmp/lutCovm.pi.aod.dat",
  can->cd(2);
  hetacalls->Draw("HIST");
 }
-
-void diagonalise(lutEntry_t& lutEntry)
-{
-  // Printf(" --- diagonalise: pt = %f, eta = %f", lutEntry.pt, lutEntry.eta);
-  TMatrixDSym m(5);
-  double fcovm[5][5];
-  for (int i = 0, k = 0; i < 5; ++i)
-    for (int j = 0; j < i + 1; ++j, ++k) {
-      fcovm[i][j] = lutEntry.covm[k];
-      fcovm[j][i] = lutEntry.covm[k];
-    }
-  m.SetMatrixArray((double*)fcovm);
-  TMatrixDSymEigen eigen(m);
-  // eigenvalues vector
-  TVectorD eigenVal = eigen.GetEigenValues();
-  for (int i = 0; i < 5; ++i)
-    lutEntry.eigval[i] = eigenVal[i];
-  // eigenvectors matrix
-  TMatrixD eigenVec = eigen.GetEigenVectors();
-  for (int i = 0; i < 5; ++i)
-    for (int j = 0; j < 5; ++j)
-      lutEntry.eigvec[i][j] = eigenVec[i][j];
-  // inverse eigenvectors matrix
-  eigenVec.Invert();
-  for (int i = 0; i < 5; ++i)
-    for (int j = 0; j < 5; ++j)
-      lutEntry.eiginv[i][j] = eigenVec[i][j];
-}