[PWGLF] add histograms for syst uncertainty on primary fraction (#13272)

Author: alcaliva

PWGLF/Tasks/Nuspex/antinucleiInJets.cxx

@@ -142,6 +142,22 @@ struct AntinucleiInJets {
   Configurable<bool> cfgCompensatePIDinTracking{"cfgCompensatePIDinTracking", false, "If true, divide tpcInnerParam by the electric charge"};
   Configurable<std::array<double, 5>> cfgBetheBlochParams{"cfgBetheBlochParams", {0.6539, 1.591, 0.8225, 2.363, 0.09}, "TPC Bethe-Bloch parameterisation for He3"};
 
+  // Configuration parameters for CCDB access and reweighting input files
+  Configurable<std::string> urlToCcdb{"urlToCcdb", "http://alice-ccdb.cern.ch", "url of the personal ccdb"};
+  Configurable<std::string> pathToFile{"pathToFile", "Users/a/alcaliva/PrimaryFractionAntip", "path to file with reweighting"};
+  Configurable<std::string> weightsProton{"weightsProton", "", "weightsProton"};
+  Configurable<std::string> weightsLambda{"weightsLambda", "", "weightsLambda"};
+  Configurable<std::string> weightsSigma{"weightsSigma", "", "weightsSigma"};
+  Configurable<std::string> weightsXi{"weightsXi", "", "weightsXi"};
+  Configurable<std::string> weightsOmega{"weightsOmega", "", "weightsOmega"};
+
+  // Reweighting histograms
+  TH1F* primaryAntiprotons;
+  TH1F* primaryAntiLambda;
+  TH1F* primaryAntiSigma;
+  TH1F* primaryAntiXi;
+  TH1F* primaryAntiOmega;
+
   // CCDB manager service for accessing condition data
   Service<o2::ccdb::BasicCCDBManager> ccdb;
 
@@ -179,6 +195,16 @@ struct AntinucleiInJets {
       itsResponse.setMCDefaultParameters();
     }
 
+    // Load reweighting histograms from CCDB if antinuclei efficiency processing is enabled
+    if (doprocessAntinucleiEfficiency) {
+      ccdb->setURL(urlToCcdb.value);
+      ccdb->setCaching(true);
+      ccdb->setLocalObjectValidityChecking();
+      ccdb->setCreatedNotAfter(std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count());
+      ccdb->setFatalWhenNull(false);
+      getReweightingHistograms(ccdb, TString(pathToFile), TString(weightsProton), TString(weightsLambda), TString(weightsSigma), TString(weightsXi), TString(weightsOmega));
+    }
+
     // Binning
     double min = 0.0;
     double max = 6.0;
@@ -361,6 +387,17 @@ struct AntinucleiInJets {
 
       // nsigmaITS for antiproton candidates
       registryMC.add("antiproton_nsigma_its_mc", "antiproton_nsigma_its_mc", HistType::kTH2F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {400, -20.0, 20.0, "n#sigma_{ITS}"}});
+
+      // Systematics on the fraction of primary antiprotons
+      registryMC.add("antip_prim_pythia", "antip_prim_pythia", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
+      registryMC.add("antip_prim_std", "antip_prim_std", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
+      registryMC.add("antip_prim_up", "antip_prim_up", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
+      registryMC.add("antip_prim_low", "antip_prim_low", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
+
+      registryMC.add("antip_sec_pythia", "antip_sec_pythia", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
+      registryMC.add("antip_sec_std", "antip_sec_std", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
+      registryMC.add("antip_sec_up", "antip_sec_up", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
+      registryMC.add("antip_sec_low", "antip_sec_low", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
     }
 
     // Systematic uncertainties (Data)
@@ -442,6 +479,27 @@ struct AntinucleiInJets {
     }
   }
 
+  void getReweightingHistograms(o2::framework::Service<o2::ccdb::BasicCCDBManager> const& ccdbObj, TString filepath, TString antip, TString antilambda, TString antisigma, TString antixi, TString antiomega)
+  {
+    TList* list = ccdbObj->get<TList>(filepath.Data());
+    if (!list) {
+      LOGP(error, "Could not retrieve the list from CCDB");
+      return;
+    }
+
+    primaryAntiprotons = static_cast<TH1F*>(list->FindObject(antip));
+    primaryAntiLambda = static_cast<TH1F*>(list->FindObject(antilambda));
+    primaryAntiSigma = static_cast<TH1F*>(list->FindObject(antisigma));
+    primaryAntiXi = static_cast<TH1F*>(list->FindObject(antixi));
+    primaryAntiOmega = static_cast<TH1F*>(list->FindObject(antiomega));
+
+    if (!primaryAntiprotons || !primaryAntiSigma || !primaryAntiLambda || !primaryAntiXi || !primaryAntiOmega) {
+      LOGP(error, "Missing one or more reweighting histograms in CCDB list");
+    }
+
+    LOGP(info, "Successfully loaded reweighting histograms from CCDB path");
+  }
+
   // Compute two unit vectors perpendicular to p
   void getPerpendicularAxis(const TVector3& p, TVector3& u, double sign)
   {
@@ -1344,6 +1402,94 @@ struct AntinucleiInJets {
           continue;
         const auto particle = track.mcParticle();
 
+        // ****************************************************************************************************************
+
+        // Systematic uncertainty on primary fraction
+        if (track.sign() < 0 && particle.pdgCode() == PDG_t::kProtonBar) {
+
+          // Primary antiprotons
+          if (particle.isPhysicalPrimary()) {
+
+            // Initialize weights
+            double wPrimStd(1.0), wPrimUp(1.0), wPrimLow(1.0);
+
+            // Weight assignment
+            if (particle.pt() < primaryAntiprotons->GetXaxis()->GetXmax()) {
+              int ipt = primaryAntiprotons->FindBin(particle.pt());
+              wPrimStd = primaryAntiprotons->GetBinContent(ipt);
+              wPrimUp = wPrimStd + primaryAntiprotons->GetBinError(ipt);
+              wPrimLow = wPrimStd - primaryAntiprotons->GetBinError(ipt);
+            }
+
+            // Fill histograms
+            registryMC.fill(HIST("antip_prim_pythia"), track.pt());
+            registryMC.fill(HIST("antip_prim_std"), track.pt(), wPrimStd);
+            registryMC.fill(HIST("antip_prim_up"), track.pt(), wPrimUp);
+            registryMC.fill(HIST("antip_prim_low"), track.pt(), wPrimLow);
+          }
+
+          // Secondary antiprotons from week decays
+          if (!particle.isPhysicalPrimary() && particle.has_mothers()) {
+
+            // Initialize weights
+            double wSecStd(1.0), wSecUp(1.0), wSecLow(1.0);
+            auto mother = mcParticles.iteratorAt(particle.mothersIds()[0]);
+
+            // Antiprotons from sigma
+            if (std::abs(mother.pdgCode()) == PDG_t::kSigmaBarMinus) {
+              if (mother.pt() < primaryAntiSigma->GetXaxis()->GetXmax()) {
+                int ipt = primaryAntiSigma->FindBin(mother.pt());
+                wSecStd = primaryAntiSigma->GetBinContent(ipt);
+                wSecUp = wSecStd + primaryAntiSigma->GetBinError(ipt);
+                wSecLow = wSecStd - primaryAntiSigma->GetBinError(ipt);
+              }
+            }
+
+            // Antiprotons from primary Lambda0
+            if (std::abs(mother.pdgCode()) == kLambda0Bar) {
+              if (mother.isPhysicalPrimary()) {
+                if (mother.pt() < primaryAntiLambda->GetXaxis()->GetXmax()) {
+                  int ipt = primaryAntiLambda->FindBin(mother.pt());
+                  wSecStd = primaryAntiLambda->GetBinContent(ipt);
+                  wSecUp = wSecStd + primaryAntiLambda->GetBinError(ipt);
+                  wSecLow = wSecStd - primaryAntiLambda->GetBinError(ipt);
+                }
+              }
+
+              // Antiprotons from secondary Lambda0 (Xi -> Lambda0)
+              if (!mother.isPhysicalPrimary()) {
+                auto grandmother = mcParticles.iteratorAt(mother.mothersIds()[0]);
+                if (std::abs(grandmother.pdgCode()) == kXiMinus) {
+                  if (grandmother.pt() < primaryAntiXi->GetXaxis()->GetXmax()) {
+                    int ipt = primaryAntiXi->FindBin(grandmother.pt());
+                    wSecStd = primaryAntiXi->GetBinContent(ipt);
+                    wSecUp = wSecStd + primaryAntiXi->GetBinError(ipt);
+                    wSecLow = wSecStd - primaryAntiXi->GetBinError(ipt);
+                  }
+                }
+
+                // Antiprotons from secondary Lambda0 (Omega -> Lambda0)
+                if (std::abs(grandmother.pdgCode()) == kOmegaMinus) {
+                  if (grandmother.pt() < primaryAntiOmega->GetXaxis()->GetXmax()) {
+                    int ipt = primaryAntiOmega->FindBin(grandmother.pt());
+                    wSecStd = primaryAntiOmega->GetBinContent(ipt);
+                    wSecUp = wSecStd + primaryAntiOmega->GetBinError(ipt);
+                    wSecLow = wSecStd - primaryAntiOmega->GetBinError(ipt);
+                  }
+                }
+              }
+            }
+
+            // Fill histograms
+            registryMC.fill(HIST("antip_sec_pythia"), track.pt());
+            registryMC.fill(HIST("antip_sec_std"), track.pt(), wSecStd);
+            registryMC.fill(HIST("antip_sec_up"), track.pt(), wSecUp);
+            registryMC.fill(HIST("antip_sec_low"), track.pt(), wSecLow);
+          }
+        }
+
+        // ****************************************************************************************************************
+
         // Select only physical primary particles
         if (!particle.isPhysicalPrimary())
           continue;