Skip to content

Commit bfde7d0

Browse files
alcalivaalcaliva
committed
[PWGLF] add histograms for syst uncertainty on primary fraction
1 parent 3703f7a commit bfde7d0

File tree

1 file changed

+148
-0
lines changed

1 file changed

+148
-0
lines changed

‎PWGLF/Tasks/Nuspex/antinucleiInJets.cxx‎

Lines changed: 148 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -142,6 +142,22 @@ struct AntinucleiInJets {
142142
Configurable<bool> cfgCompensatePIDinTracking{"cfgCompensatePIDinTracking", false, "If true, divide tpcInnerParam by the electric charge"};
143143
Configurable<std::array<double, 5>> cfgBetheBlochParams{"cfgBetheBlochParams", {0.6539, 1.591, 0.8225, 2.363, 0.09}, "TPC Bethe-Bloch parameterisation for He3"};
144144

145+
// Configuration parameters for CCDB access and reweighting input files
146+
Configurable<std::string> urlToCcdb{"urlToCcdb", "http://alice-ccdb.cern.ch", "url of the personal ccdb"};
147+
Configurable<std::string> pathToFile{"pathToFile", "Users/a/alcaliva/PrimaryFractionAntip", "path to file with reweighting"};
148+
Configurable<std::string> nameReweightingAntiprotons{"weights_proton", "", "nameReweightingAntiprotons"};
149+
Configurable<std::string> nameReweightingAntiLambda{"weights_lambda", "", "nameReweightingAntiLambda"};
150+
Configurable<std::string> nameReweightingAntiSigma{"weights_sigma", "", "nameReweightingAntiSigma"};
151+
Configurable<std::string> nameReweightingAntiXi{"weights_xi", "", "nameReweightingAntiXi"};
152+
Configurable<std::string> nameReweightingAntiOmega{"weights_omega", "", "nameReweightingAntiOmega"};
153+
154+
// Reweighting histograms
155+
TH1F* primaryAntiprotons;
156+
TH1F* primaryAntiLambda;
157+
TH1F* primaryAntiSigma;
158+
TH1F* primaryAntiXi;
159+
TH1F* primaryAntiOmega;
160+
145161
// CCDB manager service for accessing condition data
146162
Service<o2::ccdb::BasicCCDBManager> ccdb;
147163

@@ -179,6 +195,16 @@ struct AntinucleiInJets {
179195
itsResponse.setMCDefaultParameters();
180196
}
181197

198+
// Load reweighting histograms from CCDB if antinuclei efficiency processing is enabled
199+
if (doprocessAntinucleiEfficiency) {
200+
ccdb->setURL(urlToCcdb.value);
201+
ccdb->setCaching(true);
202+
ccdb->setLocalObjectValidityChecking();
203+
ccdb->setCreatedNotAfter(std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count());
204+
ccdb->setFatalWhenNull(false);
205+
getReweightingHistograms(ccdb, TString(pathToFile), TString(nameReweightingAntiprotons), TString(nameReweightingAntiLambda), TString(nameReweightingAntiSigma), TString(nameReweightingAntiXi), TString(nameReweightingAntiOmega));
206+
}
207+
182208
// Binning
183209
double min = 0.0;
184210
double max = 6.0;
@@ -361,6 +387,17 @@ struct AntinucleiInJets {
361387

362388
// nsigmaITS for antiproton candidates
363389
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}"}});
390+
391+
// Systematics on the fraction of primary antiprotons
392+
registryMC.add("antip_prim_pythia", "antip_prim_pythia", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
393+
registryMC.add("antip_prim_std", "antip_prim_std", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
394+
registryMC.add("antip_prim_up", "antip_prim_up", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
395+
registryMC.add("antip_prim_low", "antip_prim_low", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
396+
397+
registryMC.add("antip_sec_pythia", "antip_sec_pythia", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
398+
registryMC.add("antip_sec_std", "antip_sec_std", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
399+
registryMC.add("antip_sec_up", "antip_sec_up", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
400+
registryMC.add("antip_sec_low", "antip_sec_low", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
364401
}
365402

366403
// Systematic uncertainties (Data)
@@ -442,6 +479,29 @@ struct AntinucleiInJets {
442479
}
443480
}
444481

482+
void getReweightingHistograms(o2::framework::Service<o2::ccdb::BasicCCDBManager> const& ccdbObj,
483+
TString filepath, TString antip, TString antilambda, TString antisigma,
484+
TString antixi, TString antiomega)
485+
{
486+
TList* list = ccdbObj->get<TList>(filepath.Data());
487+
if (!list) {
488+
LOGP(error, "Could not retrieve the list from CCDB");
489+
return;
490+
}
491+
492+
primaryAntiprotons = static_cast<TH1F*>(list->FindObject(antip));
493+
primaryAntiLambda = static_cast<TH1F*>(list->FindObject(antilambda));
494+
primaryAntiSigma = static_cast<TH1F*>(list->FindObject(antisigma));
495+
primaryAntiXi = static_cast<TH1F*>(list->FindObject(antixi));
496+
primaryAntiOmega = static_cast<TH1F*>(list->FindObject(antiomega));
497+
498+
if (!primaryAntiprotons || !primaryAntiSigma || !primaryAntiLambda || !primaryAntiXi || !primaryAntiOmega) {
499+
LOGP(error, "Missing one or more reweighting histograms in CCDB list");
500+
}
501+
502+
LOGP(info, "Successfully loaded reweighting histograms from CCDB path");
503+
}
504+
445505
// Compute two unit vectors perpendicular to p
446506
void getPerpendicularAxis(const TVector3& p, TVector3& u, double sign)
447507
{
@@ -1344,6 +1404,94 @@ struct AntinucleiInJets {
13441404
continue;
13451405
const auto particle = track.mcParticle();
13461406

1407+
// ****************************************************************************************************************
1408+
1409+
// Systematic uncertainty on primary fraction
1410+
if (track.sign() < 0 && particle.pdgCode() == PDG_t::kProtonBar) {
1411+
1412+
// Primary antiprotons
1413+
if (particle.isPhysicalPrimary()) {
1414+
1415+
// Initialize weights
1416+
double wPrimStd(1.0), wPrimUp(1.0), wPrimLow(1.0);
1417+
1418+
// Weight assignment
1419+
if (particle.pt() < primaryAntiprotons -> GetXaxis() -> GetXmax()) {
1420+
int ipt = primaryAntiprotons -> FindBin(particle.pt());
1421+
wPrimStd = primaryAntiprotons -> GetBinContent(ipt);
1422+
wPrimUp = wPrimStd + primaryAntiprotons -> GetBinError(ipt);
1423+
wPrimLow = wPrimStd - primaryAntiprotons -> GetBinError(ipt);
1424+
}
1425+
1426+
// Fill histograms
1427+
registryMC.fill(HIST("antip_prim_pythia"), track.pt());
1428+
registryMC.fill(HIST("antip_prim_std"), track.pt(), wPrimStd);
1429+
registryMC.fill(HIST("antip_prim_up"), track.pt(), wPrimUp);
1430+
registryMC.fill(HIST("antip_prim_low"), track.pt(), wPrimLow);
1431+
}
1432+
1433+
// Secondary antiprotons from week decays
1434+
if (!particle.isPhysicalPrimary() && particle.has_mothers()) {
1435+
1436+
// Initialize weights
1437+
double wSecStd(1.0), wSecUp(1.0), wSecLow(1.0);
1438+
auto mother = mcParticles.iteratorAt(particle.mothersIds()[0]);
1439+
1440+
// Antiprotons from sigma
1441+
if (std::abs(mother.pdgCode()) == PDG_t::kSigmaBarMinus) {
1442+
if (mother.pt() < primaryAntiSigma -> GetXaxis() -> GetXmax()) {
1443+
int ipt = primaryAntiSigma -> FindBin(mother.pt());
1444+
wSecStd = primaryAntiSigma -> GetBinContent(ipt);
1445+
wSecUp = wSecStd + primaryAntiSigma -> GetBinError(ipt);
1446+
wSecLow = wSecStd - primaryAntiSigma -> GetBinError(ipt);
1447+
}
1448+
}
1449+
1450+
// Antiprotons from primary Lambda0
1451+
if (std::abs(mother.pdgCode()) == kLambda0Bar) {
1452+
if (mother.isPhysicalPrimary()) {
1453+
if (mother.pt() < primaryAntiLambda -> GetXaxis() -> GetXmax()) {
1454+
int ipt = primaryAntiLambda -> FindBin(mother.pt());
1455+
wSecStd = primaryAntiLambda -> GetBinContent(ipt);
1456+
wSecUp = wSecStd + primaryAntiLambda -> GetBinError(ipt);
1457+
wSecLow = wSecStd - primaryAntiLambda -> GetBinError(ipt);
1458+
}
1459+
}
1460+
1461+
// Antiprotons from secondary Lambda0 (Xi -> Lambda0)
1462+
if (!mother.isPhysicalPrimary()) {
1463+
auto grandmother = mcParticles.iteratorAt(mother.mothersIds()[0]);
1464+
if (std::abs(grandmother.pdgCode()) == kXiMinus) {
1465+
if (grandmother.pt() < primaryAntiXi -> GetXaxis() -> GetXmax()) {
1466+
int ipt = primaryAntiXi -> FindBin(grandmother.pt());
1467+
wSecStd = primaryAntiXi -> GetBinContent(ipt);
1468+
wSecUp = wSecStd + primaryAntiXi -> GetBinError(ipt);
1469+
wSecLow = wSecStd - primaryAntiXi -> GetBinError(ipt);
1470+
}
1471+
}
1472+
1473+
// Antiprotons from secondary Lambda0 (Omega -> Lambda0)
1474+
if (std::abs(grandmother.pdgCode()) == kOmegaMinus) {
1475+
if (grandmother.pt() < primaryAntiOmega -> GetXaxis() -> GetXmax()) {
1476+
int ipt = primaryAntiOmega -> FindBin(grandmother.pt());
1477+
wSecStd = primaryAntiOmega -> GetBinContent(ipt);
1478+
wSecUp = wSecStd + primaryAntiOmega -> GetBinError(ipt);
1479+
wSecLow = wSecStd - primaryAntiOmega -> GetBinError(ipt);
1480+
}
1481+
}
1482+
}
1483+
}
1484+
1485+
// Fill histograms
1486+
registryMC.fill(HIST("antip_sec_pythia"), track.pt());
1487+
registryMC.fill(HIST("antip_sec_std"), track.pt(), wSecStd);
1488+
registryMC.fill(HIST("antip_sec_up"), track.pt(), wSecUp);
1489+
registryMC.fill(HIST("antip_sec_low"), track.pt(), wSecLow);
1490+
}
1491+
}
1492+
1493+
// ****************************************************************************************************************
1494+
13471495
// Select only physical primary particles
13481496
if (!particle.isPhysicalPrimary())
13491497
continue;

0 commit comments

Comments
 (0)