Skip to content

Commit 2cf4029

Browse files
ChiaraDeMartin95Chiara De Martin
ChiaraDeMartin95
and
Chiara De Martin
authored
PWGLF: Fix bug in MCGen process function of cascadeflow task (AliceO2Group#7949)
* remove iterator to avoid temporary issue with index sorting in derived data * add NEvents histogram for MC collisions --------- Co-authored-by: Chiara De Martin <[email protected]>
1 parent 875e677 commit 2cf4029

File tree

1 file changed

+89
-58
lines changed

1 file changed

+89
-58
lines changed

‎PWGLF/TableProducer/Strangeness/cascadeflow.cxx‎

Lines changed: 89 additions & 58 deletions
Original file line numberDiff line numberDiff line change
@@ -37,8 +37,9 @@ using std::array;
3737
using DauTracks = soa::Join<aod::DauTrackExtras, aod::DauTrackTPCPIDs>;
3838
using CollEventPlane = soa::Join<aod::StraCollisions, aod::StraCents, aod::StraEvSels, aod::StraFT0CQVs, aod::StraRawCents, aod::StraFT0CQVsEv, aod::StraTPCQVs>::iterator;
3939
using CollEventPlaneCentralFW = soa::Join<aod::StraCollisions, aod::StraCents, aod::StraEvSels, aod::StraFT0CQVs, aod::StraRawCents, aod::StraTPCQVs>::iterator;
40+
using MCCollisionsStra = soa::Join<aod::StraMCCollisions, aod::StraMCCollMults>;
4041
using CascCandidates = soa::Join<aod::CascCollRefs, aod::CascCores, aod::CascExtras, aod::CascBBs>;
41-
using CascMCCandidates = soa::Join<aod::CascCollRefs, aod::CascCores, aod::CascExtras, aod::CascBBs, aod::CascMCMothers, aod::CascCoreMCLabels>;
42+
using CascMCCandidates = soa::Join<aod::CascCollRefs, aod::CascCores, aod::CascExtras, aod::CascBBs, aod::CascCoreMCLabels>;
4243

4344
namespace cascadev2
4445
{
@@ -147,6 +148,7 @@ struct cascadeFlow {
147148
Configurable<float> nsigmatpcPi{"nsigmatpcPi", 5, "nsigmatpcPi"};
148149
Configurable<float> mintpccrrows{"mintpccrrows", 70, "mintpccrrows"};
149150
Configurable<float> etaCascMCGen{"etaCascMCGen", 0.8, "etaCascMCGen"};
151+
Configurable<float> yCascMCGen{"yCascMCGen", 0.5, "yCascMCGen"};
150152

151153
Configurable<std::string> ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
152154
Configurable<std::vector<std::string>> modelPathsCCDBXi{"modelPathsCCDBXi", std::vector<std::string>{"Users/c/chdemart/CascadesFlow"}, "Paths of models on CCDB"};
@@ -257,9 +259,9 @@ struct cascadeFlow {
257259
return phi;
258260
}
259261

260-
HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};
261-
HistogramRegistry histosMCGen{"histosMCGen", {}, OutputObjHandlingPolicy::AnalysisObject};
262-
HistogramRegistry resolution{"resolution", {}, OutputObjHandlingPolicy::AnalysisObject};
262+
HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject, false, true};
263+
HistogramRegistry histosMCGen{"histosMCGen", {}, OutputObjHandlingPolicy::AnalysisObject, false, true};
264+
HistogramRegistry resolution{"resolution", {}, OutputObjHandlingPolicy::AnalysisObject, false, true};
263265

264266
// Tables to produce
265267
Produces<aod::CascTraining> trainingSample;
@@ -352,6 +354,7 @@ struct cascadeFlow {
352354
const AxisSpec v2Axis{200, -1., 1., "#it{v}_{2}"};
353355
const AxisSpec CentAxis{18, 0., 90., "FT0C centrality percentile"};
354356
TString hNEventsLabels[8] = {"All", "sel8", "z vrtx", "kNoSameBunchPileup", "kIsGoodZvtxFT0vsPV", "trackOccupancyInTimeRange", "kNoCollInTimeRange", "kIsGoodEventEP"};
357+
TString hNEventsLabelsMC[5] = {"All", "z vtx", ">=1RecoColl", "1Reco", "2Reco"};
355358

356359
resolution.add("QVectorsT0CTPCA", "QVectorsT0CTPCA", HistType::kTH2F, {axisQVs, CentAxis});
357360
resolution.add("QVectorsT0CTPCC", "QVectorsT0CTPCC", HistType::kTH2F, {axisQVs, CentAxis});
@@ -386,9 +389,16 @@ struct cascadeFlow {
386389
histos.add("hv2CEPvsFT0C", "hv2CEPvsFT0C", HistType::kTH2F, {CentAxis, {100, -1, 1}});
387390
histos.add("hv2CEPvsv2CSP", "hv2CEPvsV2CSP", HistType::kTH2F, {{100, -1, 1}, {100, -1, 1}});
388391

389-
histosMCGen.add("h2DGenXi", "h2DGenXi", HistType::kTH2F, {{100, 0, 00}, {200, 0, 20}});
390-
histosMCGen.add("h2DGenOmega", "h2DGenOmega", HistType::kTH2F, {{100, 0, 100}, {200, 0, 20}});
391-
histosMCGen.add("hGenEta", "hGenEta", HistType::kTH1F, {{100, -1, 1}});
392+
histosMCGen.add("h2DGenXiEta08", "h2DGenXiEta08", HistType::kTH2F, {{100, 0, 00}, {200, 0, 20}});
393+
histosMCGen.add("h2DGenOmegaEta08", "h2DGenOmegaEta08", HistType::kTH2F, {{100, 0, 100}, {200, 0, 20}});
394+
histosMCGen.add("h2DGenXiY05", "h2DGenXiY05", HistType::kTH2F, {{100, 0, 00}, {200, 0, 20}});
395+
histosMCGen.add("h2DGenOmegaY05", "h2DGenOmegaY05", HistType::kTH2F, {{100, 0, 100}, {200, 0, 20}});
396+
histosMCGen.add("hGenXiY", "hGenXiY", HistType::kTH1F, {{100, -1, 1}});
397+
histosMCGen.add("hGenOmegaY", "hGenOmegaY", HistType::kTH1F, {{100, -1, 1}});
398+
histosMCGen.add("hNEventsMC", "hNEventsMC", {HistType::kTH1F, {{5, 0.f, 5.f}}});
399+
for (Int_t n = 1; n <= histosMCGen.get<TH1>(HIST("hNEventsMC"))->GetNbinsX(); n++) {
400+
histosMCGen.get<TH1>(HIST("hNEventsMC"))->GetXaxis()->SetBinLabel(n, hNEventsLabelsMC[n - 1]);
401+
}
392402

393403
for (int iS{0}; iS < 2; ++iS) {
394404
cascadev2::hMassBeforeSelVsPt[iS] = histos.add<TH2>(Form("hMassBeforeSelVsPt%s", cascadev2::speciesNames[iS].data()), "hMassBeforeSelVsPt", HistType::kTH2F, {massCascAxis[iS], ptAxis});
@@ -773,6 +783,7 @@ struct cascadeFlow {
773783
continue;
774784

775785
auto cascMC = casc.cascMCCore_as<soa::Join<aod::CascMCCores, aod::CascMCCollRefs>>();
786+
776787
int pdgCode{cascMC.pdgCode()};
777788
if (!(std::abs(pdgCode) == 3312 && std::abs(cascMC.pdgCodeV0()) == 3122 && std::abs(cascMC.pdgCodeBachelor()) == 211) // Xi
778789
&& !(std::abs(pdgCode) == 3334 && std::abs(cascMC.pdgCodeV0()) == 3122 && std::abs(cascMC.pdgCodeBachelor()) == 321)) // Omega
@@ -852,63 +863,83 @@ struct cascadeFlow {
852863
fillAnalysedTable(coll, casc, v2CSP, v2CEP, PsiT0C, BDTresponse[0], BDTresponse[1], pdgCode);
853864
}
854865
}
855-
void processMCGen(soa::Join<aod::StraMCCollisions, aod::StraMCCollMults>::iterator const& mcCollision, soa::Join<aod::StraCollisions, aod::StraCents, aod::StraEvSels, aod::StraCollLabels> const& collisions, soa::Join<aod::CascMCCores, aod::CascMCCollRefs> const& CascMCCores)
856-
{
857-
// Generated with accepted z vertex
858-
if (TMath::Abs(mcCollision.posZ()) > cutzvertex) {
859-
return;
860-
}
861866

862-
// Check if there is at least one of the reconstructed collisions associated to this MC collision
863-
auto groupedCollisions = collisions.sliceBy(perMcCollision, mcCollision.globalIndex());
864-
int biggestNContribs = -1;
865-
int bestCollisionIndex = -1;
866-
float centrality = 100.5f;
867-
int nCollisions = 0;
868-
for (auto const& collision : groupedCollisions) {
867+
void processMCGen(MCCollisionsStra const& mcCollisions, soa::Join<aod::StraCollisions, aod::StraCents, aod::StraEvSels, aod::StraCollLabels> const& collisions, soa::Join<aod::CascMCCores, aod::CascMCCollRefs> const& CascMCCores)
868+
{
869869

870-
if (!AcceptEvent(collision)) {
871-
continue;
872-
}
873-
if (biggestNContribs < collision.multPVTotalContributors()) {
874-
biggestNContribs = collision.multPVTotalContributors();
875-
bestCollisionIndex = collision.globalIndex();
876-
centrality = collision.centFT0C();
870+
for (auto const& mcCollision : mcCollisions) {
871+
histosMCGen.fill(HIST("hNEventsMC"), 0.5);
872+
// Generated with accepted z vertex
873+
if (TMath::Abs(mcCollision.posZ()) > cutzvertex) {
874+
return;
877875
}
878-
nCollisions++;
879-
}
880-
if (nCollisions < 1) {
881-
return;
882-
}
883-
for (auto const& cascMC : CascMCCores) {
884-
if (!cascMC.has_straMCCollision())
885-
continue;
876+
histosMCGen.fill(HIST("hNEventsMC"), 1.5);
877+
// Check if there is at least one of the reconstructed collisions associated to this MC collision
886878

887-
if (!cascMC.isPhysicalPrimary())
888-
continue;
889-
890-
float ptmc = RecoDecay::sqrtSumOfSquares(cascMC.pxMC(), cascMC.pyMC());
891-
892-
float theta = std::atan(ptmc / cascMC.pzMC()); //-pi/2 < theta < pi/2
879+
auto groupedCollisions = collisions.sliceBy(perMcCollision, mcCollision.globalIndex());
893880

894-
float theta1 = 0;
895-
896-
// if pz is positive (i.e. positive rapidity): 0 < theta < pi/2
897-
if (theta > 0)
898-
theta1 = theta; // 0 < theta1/2 < pi/4 --> 0 < tan (theta1/2) < 1 --> positive eta
899-
// if pz is negative (i.e. negative rapidity): -pi/2 < theta < 0 --> we need 0 < theta1/2 < pi/2 for the ln to be defined
900-
else
901-
theta1 = TMath::Pi() + theta; // pi/2 < theta1 < pi --> pi/4 < theta1/2 < pi/2 --> 1 < tan (theta1/2) --> negative eta
902-
903-
float cascMCeta = -log(std::tan(theta1 / 2));
904-
if (TMath::Abs(cascMCeta) > etaCascMCGen)
905-
continue;
906-
histosMCGen.fill(HIST("hGenEta"), cascMCeta);
881+
int biggestNContribs = -1;
882+
int bestCollisionIndex = -1;
883+
float centrality = 100.5f;
884+
int nCollisions = 0;
885+
for (auto const& collision : groupedCollisions) {
907886

908-
if (TMath::Abs(cascMC.pdgCode()) == 3312) {
909-
histosMCGen.fill(HIST("h2dGenXi"), centrality, ptmc);
910-
} else if (TMath::Abs(cascMC.pdgCode() == 3334)) {
911-
histosMCGen.fill(HIST("h2dGenOmega"), centrality, ptmc);
887+
if (!AcceptEvent(collision)) {
888+
continue;
889+
}
890+
if (biggestNContribs < collision.multPVTotalContributors()) {
891+
biggestNContribs = collision.multPVTotalContributors();
892+
bestCollisionIndex = collision.globalIndex();
893+
centrality = collision.centFT0C();
894+
}
895+
nCollisions++;
896+
}
897+
if (nCollisions < 1) {
898+
return;
899+
}
900+
histosMCGen.fill(HIST("hNEventsMC"), 2.5);
901+
if (nCollisions == 1)
902+
histosMCGen.fill(HIST("hNEventsMC"), 3.5);
903+
else if (nCollisions == 2)
904+
histosMCGen.fill(HIST("hNEventsMC"), 4.5);
905+
for (auto const& cascMC : CascMCCores) {
906+
if (!cascMC.has_straMCCollision())
907+
continue;
908+
909+
if (!cascMC.isPhysicalPrimary())
910+
continue;
911+
912+
float ptmc = RecoDecay::sqrtSumOfSquares(cascMC.pxMC(), cascMC.pyMC());
913+
914+
float theta = std::atan(ptmc / cascMC.pzMC()); //-pi/2 < theta < pi/2
915+
916+
float theta1 = 0;
917+
918+
// if pz is positive (i.e. positive rapidity): 0 < theta < pi/2
919+
if (theta > 0)
920+
theta1 = theta; // 0 < theta1/2 < pi/4 --> 0 < tan (theta1/2) < 1 --> positive eta
921+
// if pz is negative (i.e. negative rapidity): -pi/2 < theta < 0 --> we need 0 < theta1/2 < pi/2 for the ln to be defined
922+
else
923+
theta1 = TMath::Pi() + theta; // pi/2 < theta1 < pi --> pi/4 < theta1/2 < pi/2 --> 1 < tan (theta1/2) --> negative eta
924+
925+
float cascMCeta = -log(std::tan(theta1 / 2));
926+
float cascMCy = 0;
927+
928+
if (TMath::Abs(cascMC.pdgCode()) == 3312) {
929+
cascMCy = RecoDecay::y(std::array{cascMC.pxMC(), cascMC.pyMC(), cascMC.pzMC()}, constants::physics::MassXiMinus);
930+
if (TMath::Abs(cascMCeta) < etaCascMCGen)
931+
histosMCGen.fill(HIST("h2DGenXiEta08"), centrality, ptmc);
932+
if (TMath::Abs(cascMCy) < yCascMCGen)
933+
histosMCGen.fill(HIST("h2DGenXiY05"), centrality, ptmc);
934+
histosMCGen.fill(HIST("hGenXiY"), cascMCy);
935+
} else if (TMath::Abs(cascMC.pdgCode() == 3334)) {
936+
cascMCy = RecoDecay::y(std::array{cascMC.pxMC(), cascMC.pyMC(), cascMC.pzMC()}, constants::physics::MassOmegaMinus);
937+
if (TMath::Abs(cascMCeta) < etaCascMCGen)
938+
histosMCGen.fill(HIST("h2DGenOmegaEta08"), centrality, ptmc);
939+
if (TMath::Abs(cascMCy) < yCascMCGen)
940+
histosMCGen.fill(HIST("h2DGenOmegaY05"), centrality, ptmc);
941+
histosMCGen.fill(HIST("hGenOmegaY"), cascMCy);
942+
}
912943
}
913944
}
914945
}

0 commit comments

Comments
 (0)