diff --git a/PWGLF/DataModel/LFStrangenessTables.h b/PWGLF/DataModel/LFStrangenessTables.h index 3ab94c27276..1e6f3b990b9 100644 --- a/PWGLF/DataModel/LFStrangenessTables.h +++ b/PWGLF/DataModel/LFStrangenessTables.h @@ -20,6 +20,7 @@ #include "Common/DataModel/Multiplicity.h" #include "Common/DataModel/Centrality.h" #include "Common/DataModel/Qvectors.h" +#include "Common/DataModel/McCollisionExtra.h" #include "PWGLF/DataModel/EPCalibrationTables.h" namespace o2::aod @@ -294,6 +295,9 @@ DECLARE_SOA_COLUMN(PzPosMC, pzPosMC, float); //! V0 positive DECLARE_SOA_COLUMN(PxNegMC, pxNegMC, float); //! V0 positive daughter px (GeV/c) DECLARE_SOA_COLUMN(PyNegMC, pyNegMC, float); //! V0 positive daughter py (GeV/c) DECLARE_SOA_COLUMN(PzNegMC, pzNegMC, float); //! V0 positive daughter pz (GeV/c) +DECLARE_SOA_COLUMN(PxMC, pxMC, float); //! V0 px (GeV/c) +DECLARE_SOA_COLUMN(PyMC, pyMC, float); //! V0 py (GeV/c) +DECLARE_SOA_COLUMN(PzMC, pzMC, float); //! V0 pz (GeV/c) //______________________________________________________ // Binned content for generated particles: derived data @@ -464,6 +468,22 @@ DECLARE_SOA_DYNAMIC_COLUMN(IsPhotonTPConly, isPhotonTPConly, //! is tpc-only pho [](uint8_t V0Type) -> bool { return V0Type & (1 << 1); }); DECLARE_SOA_DYNAMIC_COLUMN(IsCollinear, isCollinear, //! is collinear V0 [](uint8_t V0Type) -> bool { return V0Type & (1 << 2); }); + +DECLARE_SOA_DYNAMIC_COLUMN(RapidityMC, rapidityMC, //! rapidity (0:K0, 1:L, 2:Lbar) + [](float PxMC, float PyMC, float PzMC, int value) -> float { + if (value == 0) + return RecoDecay::y(std::array{PxMC, PyMC, PzMC}, o2::constants::physics::MassKaonNeutral); + if (value == 1 || value == 2) + return RecoDecay::y(std::array{PxMC, PyMC, PzMC}, o2::constants::physics::MassLambda); + return 0.0f; + }); + +DECLARE_SOA_DYNAMIC_COLUMN(NegativePtMC, negativeptMC, //! negative daughter pT + [](float pxnegMC, float pynegMC) -> float { return RecoDecay::sqrtSumOfSquares(pxnegMC, pynegMC); }); +DECLARE_SOA_DYNAMIC_COLUMN(PositivePtMC, positiveptMC, //! positive daughter pT + [](float pxposMC, float pyposMC) -> float { return RecoDecay::sqrtSumOfSquares(pxposMC, pyposMC); }); +DECLARE_SOA_DYNAMIC_COLUMN(PtMC, ptMC, //! V0 pT + [](float pxMC, float pyMC) -> float { return RecoDecay::sqrtSumOfSquares(pxMC, pyMC); }); } // namespace v0data DECLARE_SOA_TABLE(V0Indices, "AOD", "V0INDEX", //! index table when using AO2Ds @@ -666,6 +686,19 @@ DECLARE_SOA_TABLE_VERSIONED(V0MCCores_001, "AOD", "V0MCCORE", 1, //! debug infor v0data::PxPosMC, v0data::PyPosMC, v0data::PzPosMC, v0data::PxNegMC, v0data::PyNegMC, v0data::PzNegMC); +DECLARE_SOA_TABLE_VERSIONED(V0MCCores_002, "AOD", "V0MCCORE", 2, //! debug information + v0data::ParticleIdMC, //! MC properties of the V0 for posterior analysis + v0data::PDGCode, v0data::PDGCodeMother, + v0data::PDGCodePositive, v0data::PDGCodeNegative, + v0data::IsPhysicalPrimary, v0data::XMC, v0data::YMC, v0data::ZMC, + v0data::PxPosMC, v0data::PyPosMC, v0data::PzPosMC, + v0data::PxNegMC, v0data::PyNegMC, v0data::PzNegMC, + v0data::PxMC, v0data::PyMC, v0data::PzMC, + v0data::RapidityMC, + v0data::NegativePtMC, + v0data::PositivePtMC, + v0data::PtMC); + DECLARE_SOA_TABLE(StoredV0MCCores_000, "AOD", "V0MCCORE", //! MC properties of the V0 for posterior analysis v0data::PDGCode, v0data::PDGCodeMother, v0data::PDGCodePositive, v0data::PDGCodeNegative, @@ -683,6 +716,16 @@ DECLARE_SOA_TABLE_VERSIONED(StoredV0MCCores_001, "AOD", "V0MCCORE", 1, //! debug v0data::PxNegMC, v0data::PyNegMC, v0data::PzNegMC, o2::soa::Marker<1>); +DECLARE_SOA_TABLE_VERSIONED(StoredV0MCCores_002, "AOD", "V0MCCORE", 2, //! debug information + v0data::ParticleIdMC, //! MC properties of the V0 for posterior analysis + v0data::PDGCode, v0data::PDGCodeMother, + v0data::PDGCodePositive, v0data::PDGCodeNegative, + v0data::IsPhysicalPrimary, v0data::XMC, v0data::YMC, v0data::ZMC, + v0data::PxPosMC, v0data::PyPosMC, v0data::PzPosMC, + v0data::PxNegMC, v0data::PyNegMC, v0data::PzNegMC, + v0data::PxMC, v0data::PyMC, v0data::PzMC, + o2::soa::Marker<1>); + DECLARE_SOA_TABLE(V0MCCollRefs, "AOD", "V0MCCOLLREF", //! refers MC candidate back to proper MC Collision o2::soa::Index<>, v0data::StraMCCollisionId, o2::soa::Marker<2>); @@ -696,8 +739,8 @@ DECLARE_SOA_TABLE(V0MCMothers, "AOD", "V0MCMOTHER", //! optional table for MC mo DECLARE_SOA_TABLE(StoredV0MCMothers, "AOD1", "V0MCMOTHER", //! optional table for MC mothers o2::soa::Index<>, v0data::MotherMCPartId, o2::soa::Marker<1>); -using V0MCCores = V0MCCores_001; -using StoredV0MCCores = StoredV0MCCores_001; +using V0MCCores = V0MCCores_002; +using StoredV0MCCores = StoredV0MCCores_002; using V0Index = V0Indices::iterator; using V0Core = V0Cores::iterator; @@ -708,6 +751,7 @@ using V0fCDatas = soa::Join; using V0fCData = V0fCDatas::iterator; using V0MCDatas = soa::Join; using V0MCData = V0MCDatas::iterator; +using V0MCCore = V0MCCores::iterator; // definitions of indices for interlink tables namespace v0data @@ -997,6 +1041,24 @@ DECLARE_SOA_DYNAMIC_COLUMN(BachelorEta, bacheloreta, //! bachelor daughter eta [](float PxPos, float PyPos, float PzPos) -> float { return RecoDecay::eta(std::array{PxPos, PyPos, PzPos}); }); DECLARE_SOA_DYNAMIC_COLUMN(BachelorPhi, bachelorphi, //! bachelor daughter phi [](float PxPos, float PyPos) -> float { return RecoDecay::phi(PxPos, PyPos); }); + +DECLARE_SOA_DYNAMIC_COLUMN(RapidityMC, rapidityMC, //! rapidity (0, 1: Xi; 2, 3: Omega) + [](float PxMC, float PyMC, float PzMC, int value) -> float { + if (value == 0 || value == 1) + return RecoDecay::y(std::array{PxMC, PyMC, PzMC}, o2::constants::physics::MassXiMinus); + if (value == 2 || value == 3) + return RecoDecay::y(std::array{PxMC, PyMC, PzMC}, o2::constants::physics::MassOmegaMinus); + return 0.0f; + }); + +DECLARE_SOA_DYNAMIC_COLUMN(NegativePtMC, negativeptMC, //! negative daughter pT + [](float pxNegMC, float pyNegMC) -> float { return RecoDecay::sqrtSumOfSquares(pxNegMC, pyNegMC); }); +DECLARE_SOA_DYNAMIC_COLUMN(PositivePtMC, positiveptMC, //! positive daughter pT + [](float pxPosMC, float pyPosMC) -> float { return RecoDecay::sqrtSumOfSquares(pxPosMC, pyPosMC); }); +DECLARE_SOA_DYNAMIC_COLUMN(BachelorPtMC, bachelorptMC, //! bachelor daughter pT + [](float pxBachMC, float pyBachMC) -> float { return RecoDecay::sqrtSumOfSquares(pxBachMC, pyBachMC); }); +DECLARE_SOA_DYNAMIC_COLUMN(PtMC, ptMC, //! cascade pT + [](float pxMC, float pyMC) -> float { return RecoDecay::sqrtSumOfSquares(pxMC, pyMC); }); } // namespace cascdata //______________________________________________________ @@ -1190,7 +1252,12 @@ DECLARE_SOA_TABLE(CascMCCores, "AOD", "CASCMCCORE", //! bachelor-baryon correlat cascdata::PxPosMC, cascdata::PyPosMC, cascdata::PzPosMC, cascdata::PxNegMC, cascdata::PyNegMC, cascdata::PzNegMC, cascdata::PxBachMC, cascdata::PyBachMC, cascdata::PzBachMC, - cascdata::PxMC, cascdata::PyMC, cascdata::PzMC); + cascdata::PxMC, cascdata::PyMC, cascdata::PzMC, + cascdata::RapidityMC, + cascdata::NegativePtMC, + cascdata::PositivePtMC, + cascdata::BachelorPtMC, + cascdata::PtMC); namespace cascdata { diff --git a/PWGLF/TableProducer/Strangeness/Converters/CMakeLists.txt b/PWGLF/TableProducer/Strangeness/Converters/CMakeLists.txt index 660365da5ac..defa3b1abf9 100644 --- a/PWGLF/TableProducer/Strangeness/Converters/CMakeLists.txt +++ b/PWGLF/TableProducer/Strangeness/Converters/CMakeLists.txt @@ -27,4 +27,9 @@ o2physics_add_dpl_workflow(strarawcentsconverter2v4 o2physics_add_dpl_workflow(v0coresconverter SOURCES v0coresconverter.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(v0coresconverter2 + SOURCES v0coresconverter2.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) \ No newline at end of file diff --git a/PWGLF/TableProducer/Strangeness/Converters/v0coresconverter2.cxx b/PWGLF/TableProducer/Strangeness/Converters/v0coresconverter2.cxx new file mode 100644 index 00000000000..c918f227525 --- /dev/null +++ b/PWGLF/TableProducer/Strangeness/Converters/v0coresconverter2.cxx @@ -0,0 +1,52 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" + +using namespace o2; +using namespace o2::framework; + +// Converts V0 version 001 to 002 +struct v0coresconverter2 { + Produces v0MCCores_002; + + void process(aod::V0MCCores_001 const& v0MCCores_001) + { + for (auto& values : v0MCCores_001) { + v0MCCores_002(0, + values.pdgCode(), + values.pdgCodeMother(), + values.pdgCodePositive(), + values.pdgCodeNegative(), + values.isPhysicalPrimary(), + values.xMC(), + values.yMC(), + values.zMC(), + values.pxPosMC(), + values.pyPosMC(), + values.pzPosMC(), + values.pxNegMC(), + values.pyNegMC(), + values.pzNegMC(), + values.pxPosMC() + values.pxNegMC(), + values.pyPosMC() + values.pyNegMC(), + values.pzPosMC() + values.pzNegMC()); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc)}; +} diff --git a/PWGLF/TableProducer/Strangeness/cascadebuilder.cxx b/PWGLF/TableProducer/Strangeness/cascadebuilder.cxx index e1c0cfed9e3..82af77421f0 100644 --- a/PWGLF/TableProducer/Strangeness/cascadebuilder.cxx +++ b/PWGLF/TableProducer/Strangeness/cascadebuilder.cxx @@ -194,6 +194,8 @@ struct cascadeBuilder { Configurable kfDoDCAFitterPreMinimV0{"kfDoDCAFitterPreMinimV0", true, "KF: do DCAFitter pre-optimization before KF fit to include material corrections for V0"}; Configurable kfDoDCAFitterPreMinimCasc{"kfDoDCAFitterPreMinimCasc", true, "KF: do DCAFitter pre-optimization before KF fit to include material corrections for Xi"}; + // for using cascade momentum at prim. vtx + Configurable useCascadeMomentumAtPrimVtx{"useCascadeMomentumAtPrimVtx", false, "if enabled, store cascade momentum at prim. vtx instead of decay point (= default)"}; // for topo var QA struct : ConfigurableGroup { ConfigurableAxis axisTopoVarPointingAngle{"axisConfigurations.axisTopoVarPointingAngle", {50, 0.0, 1.0}, "pointing angle"}; @@ -1090,6 +1092,13 @@ struct cascadeBuilder { o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, lCascadeTrack, 2.f, matCorrCascade, &dcaInfo); cascadecandidate.cascDCAxy = dcaInfo[0]; cascadecandidate.cascDCAz = dcaInfo[1]; + if (useCascadeMomentumAtPrimVtx) { + lCascadeTrack.getPxPyPzGlo(cascadecandidate.cascademom); + } else { + cascadecandidate.cascademom[0] = cascadecandidate.bachP[0] + cascadecandidate.v0mompos[0] + cascadecandidate.v0momneg[0]; + cascadecandidate.cascademom[1] = cascadecandidate.bachP[1] + cascadecandidate.v0mompos[1] + cascadecandidate.v0momneg[1]; + cascadecandidate.cascademom[2] = cascadecandidate.bachP[2] + cascadecandidate.v0mompos[2] + cascadecandidate.v0momneg[2]; + } // Calculate masses a priori cascadecandidate.mXi = RecoDecay::m(array{array{cascadecandidate.bachP[0], cascadecandidate.bachP[1], cascadecandidate.bachP[2]}, array{v0.pxpos() + v0.pxneg(), v0.pypos() + v0.pyneg(), v0.pzpos() + v0.pzneg()}}, array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassLambda}); @@ -1555,9 +1564,7 @@ struct cascadeBuilder { cascadecandidate.v0mompos[0], cascadecandidate.v0mompos[1], cascadecandidate.v0mompos[2], cascadecandidate.v0momneg[0], cascadecandidate.v0momneg[1], cascadecandidate.v0momneg[2], cascadecandidate.bachP[0], cascadecandidate.bachP[1], cascadecandidate.bachP[2], - cascadecandidate.bachP[0] + cascadecandidate.v0mompos[0] + cascadecandidate.v0momneg[0], // <--- redundant but ok - cascadecandidate.bachP[1] + cascadecandidate.v0mompos[1] + cascadecandidate.v0momneg[1], // <--- redundant but ok - cascadecandidate.bachP[2] + cascadecandidate.v0mompos[2] + cascadecandidate.v0momneg[2], // <--- redundant but ok + cascadecandidate.cascademom[0], cascadecandidate.cascademom[1], cascadecandidate.cascademom[2], cascadecandidate.v0dcadau, cascadecandidate.dcacascdau, cascadecandidate.v0dcapostopv, cascadecandidate.v0dcanegtopv, cascadecandidate.bachDCAxy, cascadecandidate.cascDCAxy, cascadecandidate.cascDCAz); // <--- no corresponding stratrack information available diff --git a/PWGLF/TableProducer/Strangeness/cascademcbuilder.cxx b/PWGLF/TableProducer/Strangeness/cascademcbuilder.cxx index 4084afdcae4..729ac3e27da 100644 --- a/PWGLF/TableProducer/Strangeness/cascademcbuilder.cxx +++ b/PWGLF/TableProducer/Strangeness/cascademcbuilder.cxx @@ -76,6 +76,9 @@ struct cascademcbuilder { int pdgCodeNegative; int pdgCodeBachelor; bool isPhysicalPrimary; + int processPositive = -1; + int processNegative = -1; + int processBachelor = -1; std::array xyz; std::array lxyz; std::array posP; @@ -89,8 +92,6 @@ struct cascademcbuilder { mcCascinfo thisInfo; //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* - void init(InitContext const&) {} - template void generateCascadeMCinfo(TCascadeTable cascTable, TMCParticleTable mcParticles) { @@ -142,6 +143,9 @@ struct cascademcbuilder { thisInfo.bachP[0] = lMCBachTrack.px(); thisInfo.bachP[1] = lMCBachTrack.py(); thisInfo.bachP[2] = lMCBachTrack.pz(); + thisInfo.processPositive = lMCPosTrack.getProcess(); + thisInfo.processNegative = lMCNegTrack.getProcess(); + thisInfo.processBachelor = lMCBachTrack.getProcess(); // Step 1: check if the mother is the same, go up a level if (lMCNegTrack.has_mothers() && lMCPosTrack.has_mothers()) { @@ -194,8 +198,8 @@ struct cascademcbuilder { thisInfo.label, thisInfo.motherLabel); // Mark mcParticle as recoed (no searching necessary afterwards) - if (thisInfo.motherLabel > -1) { - mcParticleIsReco[thisInfo.motherLabel] = true; + if (thisInfo.label > -1) { + mcParticleIsReco[thisInfo.label] = true; } if (populateCascMCCoresSymmetric) { @@ -216,10 +220,7 @@ struct cascademcbuilder { // step 1: check if this element is already provided in the table // using the packedIndices variable calculated above for (uint32_t ii = 0; ii < mcCascinfos.size(); ii++) { - if ( - thisInfo.mcParticlePositive == mcCascinfos[ii].mcParticlePositive && mcCascinfos[ii].mcParticlePositive > 0 && - thisInfo.mcParticleNegative == mcCascinfos[ii].mcParticleNegative && mcCascinfos[ii].mcParticleNegative > 0 && - thisInfo.mcParticleBachelor == mcCascinfos[ii].mcParticleBachelor && mcCascinfos[ii].mcParticleBachelor > 0) { + if (thisInfo.label == mcCascinfos[ii].label && mcCascinfos[ii].label > -1) { thisCascMCCoreIndex = ii; break; // this exists already in list } @@ -272,6 +273,54 @@ struct cascademcbuilder { thisInfo.momentum[0] = mcParticle.px(); thisInfo.momentum[1] = mcParticle.py(); thisInfo.momentum[2] = mcParticle.pz(); + thisInfo.label = mcParticle.globalIndex(); + + if (mcParticle.has_daughters()) { + auto const& daughters = mcParticle.template daughters_as(); + for (auto& dau : daughters) { + if (dau.getProcess() != 4) // check whether the daughter comes from a decay + continue; + + if (TMath::Abs(dau.pdgCode()) == 211 || TMath::Abs(dau.pdgCode()) == 321) { + thisInfo.pdgCodeBachelor = dau.pdgCode(); + thisInfo.bachP[0] = dau.px(); + thisInfo.bachP[1] = dau.py(); + thisInfo.bachP[2] = dau.pz(); + thisInfo.xyz[0] = dau.vx(); + thisInfo.xyz[1] = dau.vy(); + thisInfo.xyz[2] = dau.vz(); + thisInfo.mcParticleBachelor = dau.globalIndex(); + } + if (TMath::Abs(dau.pdgCode()) == 2212) { + thisInfo.pdgCodeV0 = dau.pdgCode(); + + for (auto& v0Dau : dau.template daughters_as()) { + if (v0Dau.getProcess() != 4) + continue; + + if (v0Dau.pdgCode() > 0) { + thisInfo.pdgCodePositive = v0Dau.pdgCode(); + thisInfo.processPositive = v0Dau.getProcess(); + thisInfo.posP[0] = v0Dau.px(); + thisInfo.posP[1] = v0Dau.py(); + thisInfo.posP[2] = v0Dau.pz(); + thisInfo.lxyz[0] = v0Dau.vx(); + thisInfo.lxyz[1] = v0Dau.vy(); + thisInfo.lxyz[2] = v0Dau.vz(); + thisInfo.mcParticlePositive = v0Dau.globalIndex(); + } + if (v0Dau.pdgCode() < 0) { + thisInfo.pdgCodeNegative = v0Dau.pdgCode(); + thisInfo.processNegative = v0Dau.getProcess(); + thisInfo.negP[0] = v0Dau.px(); + thisInfo.negP[1] = v0Dau.py(); + thisInfo.negP[2] = v0Dau.pz(); + thisInfo.mcParticleNegative = v0Dau.globalIndex(); + } + } + } + } + } // if I got here, it means this MC particle was not recoed and is of interest. Add it please mcCascinfos.push_back(thisInfo); diff --git a/PWGLF/TableProducer/Strangeness/lambdakzeromcbuilder.cxx b/PWGLF/TableProducer/Strangeness/lambdakzeromcbuilder.cxx index 8fe57bc144f..a2324681040 100644 --- a/PWGLF/TableProducer/Strangeness/lambdakzeromcbuilder.cxx +++ b/PWGLF/TableProducer/Strangeness/lambdakzeromcbuilder.cxx @@ -77,22 +77,40 @@ struct lambdakzeromcbuilder { h->GetXaxis()->SetBinLabel(2, "V0MCCores population"); h->GetXaxis()->SetBinLabel(3, "x check: duplicates"); h->GetXaxis()->SetBinLabel(4, "x check: unique"); + + auto hK0s = histos.add("hStatisticsK0s", "hBuildingStatisticsK0s", kTH1F, {{3, -0.5, 2.5f}}); + hK0s->GetXaxis()->SetBinLabel(1, "MC associated to reco."); + hK0s->GetXaxis()->SetBinLabel(2, "Not originating from decay"); + hK0s->GetXaxis()->SetBinLabel(3, "MC with un-recoed V0"); + + auto hLambda = histos.add("hStatisticsLambda", "hBuildingStatisticsLambda", kTH1F, {{3, -0.5, 2.5f}}); + hLambda->GetXaxis()->SetBinLabel(1, "MC associated to reco."); + hLambda->GetXaxis()->SetBinLabel(2, "Not originating from decay"); + hLambda->GetXaxis()->SetBinLabel(3, "MC with un-recoed V0"); + + auto hAlambda = histos.add("hStatisticsAlambda", "hBuildingStatisticsAlambda", kTH1F, {{3, -0.5, 2.5f}}); + hAlambda->GetXaxis()->SetBinLabel(1, "MC associated to reco."); + hAlambda->GetXaxis()->SetBinLabel(2, "Not originating from decay"); + hAlambda->GetXaxis()->SetBinLabel(3, "MC with un-recoed V0"); } //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* // Helper struct to contain V0MCCore information prior to filling struct mcV0info { - int label; - int motherLabel; - int pdgCode; - int pdgCodeMother; - int pdgCodePositive; - int pdgCodeNegative; - int mcCollision; - bool isPhysicalPrimary; + int label = -1; + int motherLabel = -1; + int pdgCode = 0; + int pdgCodeMother = 0; + int pdgCodePositive = 0; + int pdgCodeNegative = 0; + int mcCollision = -1; + bool isPhysicalPrimary = false; + int processPositive = -1; + int processNegative = -1; std::array xyz; std::array posP; std::array negP; + std::array momentum; uint64_t packedMcParticleIndices; }; mcV0info thisInfo; @@ -124,17 +142,21 @@ struct lambdakzeromcbuilder { thisInfo.xyz[0] = thisInfo.xyz[1] = thisInfo.xyz[2] = 0.0f; thisInfo.posP[0] = thisInfo.posP[1] = thisInfo.posP[2] = 0.0f; thisInfo.negP[0] = thisInfo.negP[1] = thisInfo.negP[2] = 0.0f; + thisInfo.momentum[0] = thisInfo.momentum[1] = thisInfo.momentum[2] = 0.0f; auto lNegTrack = v0.negTrack_as(); auto lPosTrack = v0.posTrack_as(); // Association check // There might be smarter ways of doing this in the future if (lNegTrack.has_mcParticle() && lPosTrack.has_mcParticle()) { - thisInfo.packedMcParticleIndices = combineProngIndices(lPosTrack.mcParticleId(), lNegTrack.mcParticleId()); auto lMCNegTrack = lNegTrack.mcParticle_as(); auto lMCPosTrack = lPosTrack.mcParticle_as(); + + thisInfo.packedMcParticleIndices = combineProngIndices(lPosTrack.mcParticleId(), lNegTrack.mcParticleId()); thisInfo.pdgCodePositive = lMCPosTrack.pdgCode(); thisInfo.pdgCodeNegative = lMCNegTrack.pdgCode(); + thisInfo.processPositive = lMCPosTrack.getProcess(); + thisInfo.processNegative = lMCNegTrack.getProcess(); thisInfo.posP[0] = lMCPosTrack.px(); thisInfo.posP[1] = lMCPosTrack.py(); thisInfo.posP[2] = lMCPosTrack.pz(); @@ -146,17 +168,50 @@ struct lambdakzeromcbuilder { for (auto& lPosMother : lMCPosTrack.mothers_as()) { if (lNegMother.globalIndex() == lPosMother.globalIndex()) { thisInfo.label = lNegMother.globalIndex(); + thisInfo.xyz[0] = lMCPosTrack.vx(); + thisInfo.xyz[1] = lMCPosTrack.vy(); + thisInfo.xyz[2] = lMCPosTrack.vz(); + + // MC pos. and neg. daughters are the same! Looking for replacement... + if (lMCPosTrack.globalIndex() == lMCNegTrack.globalIndex()) { + auto const& daughters = lNegMother.daughters_as(); + for (auto& ldau : daughters) { + // check if the candidate originate from a decay + // if not, this is not a suitable candidate for one of the decay daughters + if (ldau.getProcess() != 4) // see TMCProcess.h + continue; + + if (lMCPosTrack.pdgCode() < 0 && ldau.pdgCode() > 0) { // the positive track needs to be changed + thisInfo.pdgCodePositive = ldau.pdgCode(); + thisInfo.processPositive = ldau.getProcess(); + thisInfo.posP[0] = ldau.px(); + thisInfo.posP[1] = ldau.py(); + thisInfo.posP[2] = ldau.pz(); + thisInfo.xyz[0] = ldau.vx(); + thisInfo.xyz[1] = ldau.vy(); + thisInfo.xyz[2] = ldau.vz(); + } + if (lMCNegTrack.pdgCode() > 0 && ldau.pdgCode() < 0) { // the negative track needs to be changed + thisInfo.pdgCodeNegative = ldau.pdgCode(); + thisInfo.processNegative = ldau.getProcess(); + thisInfo.negP[0] = ldau.px(); + thisInfo.negP[1] = ldau.py(); + thisInfo.negP[2] = ldau.pz(); + } + } + } if (lNegMother.has_mcCollision()) { thisInfo.mcCollision = lNegMother.mcCollisionId(); // save this reference, please } // acquire information - thisInfo.xyz[0] = lMCPosTrack.vx(); - thisInfo.xyz[1] = lMCPosTrack.vy(); - thisInfo.xyz[2] = lMCPosTrack.vz(); thisInfo.pdgCode = lNegMother.pdgCode(); thisInfo.isPhysicalPrimary = lNegMother.isPhysicalPrimary(); + thisInfo.momentum[0] = lNegMother.px(); + thisInfo.momentum[1] = lNegMother.py(); + thisInfo.momentum[2] = lNegMother.pz(); + if (lNegMother.has_mothers()) { for (auto& lNegGrandMother : lNegMother.mothers_as()) { thisInfo.pdgCodeMother = lNegGrandMother.pdgCode(); @@ -173,8 +228,8 @@ struct lambdakzeromcbuilder { thisInfo.label, thisInfo.motherLabel); // Mark mcParticle as recoed (no searching necessary afterwards) - if (thisInfo.motherLabel > -1) { - mcParticleIsReco[thisInfo.motherLabel] = true; + if (thisInfo.label > -1) { + mcParticleIsReco[thisInfo.label] = true; } // ---] Symmetric populate [--- @@ -187,7 +242,8 @@ struct lambdakzeromcbuilder { thisInfo.pdgCodeMother, thisInfo.pdgCodePositive, thisInfo.pdgCodeNegative, thisInfo.isPhysicalPrimary, thisInfo.xyz[0], thisInfo.xyz[1], thisInfo.xyz[2], thisInfo.posP[0], thisInfo.posP[1], thisInfo.posP[2], - thisInfo.negP[0], thisInfo.negP[1], thisInfo.negP[2]); + thisInfo.negP[0], thisInfo.negP[1], thisInfo.negP[2], + thisInfo.momentum[0], thisInfo.momentum[1], thisInfo.momentum[2]); v0mccollref(thisInfo.mcCollision); // n.b. placing the interlink index here allows for the writing of @@ -204,17 +260,45 @@ struct lambdakzeromcbuilder { // step 1: check if this element is already provided in the table // using the packedIndices variable calculated above for (uint32_t ii = 0; ii < mcV0infos.size(); ii++) { - if (thisInfo.packedMcParticleIndices == mcV0infos[ii].packedMcParticleIndices && mcV0infos[ii].packedMcParticleIndices > 0) { + if (thisInfo.label == mcV0infos[ii].label && mcV0infos[ii].label > -1) { thisV0MCCoreIndex = ii; histos.fill(HIST("hBuildingStatistics"), 2.0f); // found break; // this exists already in list } } - if (thisV0MCCoreIndex < 0) { + if (thisV0MCCoreIndex < 0 && thisInfo.label > -1) { // this V0MCCore does not exist yet. Create it and reference it histos.fill(HIST("hBuildingStatistics"), 3.0f); // new thisV0MCCoreIndex = mcV0infos.size(); mcV0infos.push_back(thisInfo); + + // For bookkeeping + if (thisInfo.label > -1 && thisInfo.isPhysicalPrimary) { + float ymc = 1e3; + if (thisInfo.pdgCode == 310) + ymc = RecoDecay::y(std::array{thisInfo.posP[0] + thisInfo.negP[0], thisInfo.posP[1] + thisInfo.negP[1], thisInfo.posP[2] + thisInfo.negP[2]}, o2::constants::physics::MassKaonNeutral); + else if (TMath::Abs(thisInfo.pdgCode) == 3122) + ymc = RecoDecay::y(std::array{thisInfo.posP[0] + thisInfo.negP[0], thisInfo.posP[1] + thisInfo.negP[1], thisInfo.posP[2] + thisInfo.negP[2]}, o2::constants::physics::MassLambda); + + if (thisInfo.pdgCode == 310 && TMath::Abs(ymc) < rapidityWindow) { + histos.fill(HIST("hStatisticsK0s"), 0.0f); // found + if (thisInfo.processPositive != 4 || thisInfo.processNegative != 4) { + histos.fill(HIST("hStatisticsK0s"), 1.0f); // Not originating from decay + } + } + if (thisInfo.pdgCode == 3122 && TMath::Abs(ymc) < rapidityWindow) { + histos.fill(HIST("hStatisticsLambda"), 0.0f); // found + if (thisInfo.processPositive != 4 || thisInfo.processNegative != 4) { + histos.fill(HIST("hStatisticsLambda"), 1.0f); // Not originating from decay + } + } + if (thisInfo.pdgCode == -3122 && TMath::Abs(ymc) < rapidityWindow) { + histos.fill(HIST("hStatisticsAlambda"), 0.0f); // found + if (thisInfo.processPositive != 4 || thisInfo.processNegative != 4) { + histos.fill(HIST("hStatisticsAlambda"), 1.0f); // Not originating from decay + } + } + } } v0CoreMCLabels(thisV0MCCoreIndex); // interlink index } @@ -228,13 +312,14 @@ struct lambdakzeromcbuilder { thisInfo.label = -1; thisInfo.motherLabel = -1; thisInfo.pdgCode = 0; - thisInfo.pdgCodeMother = 0; - thisInfo.pdgCodePositive = 0; - thisInfo.pdgCodeNegative = 0; + thisInfo.pdgCodeMother = -1; + thisInfo.pdgCodePositive = -1; + thisInfo.pdgCodeNegative = -1; thisInfo.mcCollision = -1; thisInfo.xyz[0] = thisInfo.xyz[1] = thisInfo.xyz[2] = 0.0f; thisInfo.posP[0] = thisInfo.posP[1] = thisInfo.posP[2] = 0.0f; thisInfo.negP[0] = thisInfo.negP[1] = thisInfo.negP[2] = 0.0f; + thisInfo.momentum[0] = thisInfo.momentum[1] = thisInfo.momentum[2] = 0.0f; if (mcParticleIsReco[mcParticle.globalIndex()] == true) continue; // skip if already created in list @@ -249,16 +334,65 @@ struct lambdakzeromcbuilder { (addGeneratedGamma && mcParticle.pdgCode() == 22)) { thisInfo.pdgCode = mcParticle.pdgCode(); thisInfo.isPhysicalPrimary = mcParticle.isPhysicalPrimary(); + thisInfo.label = mcParticle.globalIndex(); if (mcParticle.has_mcCollision()) { thisInfo.mcCollision = mcParticle.mcCollisionId(); // save this reference, please } - // guarantee compressibility: keep momentum entirely in the positive prong - // WARNING: THIS IS AS ARBITRARY AS IT GETS but should be ok - thisInfo.posP[0] = mcParticle.px(); - thisInfo.posP[1] = mcParticle.py(); - thisInfo.posP[2] = mcParticle.pz(); + // + thisInfo.momentum[0] = mcParticle.px(); + thisInfo.momentum[1] = mcParticle.py(); + thisInfo.momentum[2] = mcParticle.pz(); + + if (mcParticle.has_mothers()) { + auto const& mother = mcParticle.mothers_first_as(); + thisInfo.pdgCodeMother = mother.pdgCode(); + thisInfo.motherLabel = mother.globalIndex(); + } + if (mcParticle.has_daughters()) { + auto const& daughters = mcParticle.daughters_as(); + + for (auto& dau : daughters) { + if (dau.getProcess() != 4) + continue; + + if (dau.pdgCode() > 0) { + thisInfo.pdgCodePositive = dau.pdgCode(); + thisInfo.processPositive = dau.getProcess(); + thisInfo.posP[0] = dau.px(); + thisInfo.posP[1] = dau.py(); + thisInfo.posP[2] = dau.pz(); + thisInfo.xyz[0] = dau.vx(); + thisInfo.xyz[1] = dau.vy(); + thisInfo.xyz[2] = dau.vz(); + } + if (dau.pdgCode() < 0) { + thisInfo.pdgCodeNegative = dau.pdgCode(); + thisInfo.processNegative = dau.getProcess(); + thisInfo.negP[0] = dau.px(); + thisInfo.negP[1] = dau.py(); + thisInfo.negP[2] = dau.pz(); + } + } + } + + // For bookkeeping + float ymc = 1e3; + if (mcParticle.pdgCode() == 310) + ymc = RecoDecay::y(std::array{thisInfo.posP[0] + thisInfo.negP[0], thisInfo.posP[1] + thisInfo.negP[0], thisInfo.posP[2] + thisInfo.negP[2]}, o2::constants::physics::MassKaonNeutral); + else if (TMath::Abs(mcParticle.pdgCode()) == 3122) + ymc = RecoDecay::y(std::array{thisInfo.posP[0] + thisInfo.negP[0], thisInfo.posP[1] + thisInfo.negP[0], thisInfo.posP[2] + thisInfo.negP[2]}, o2::constants::physics::MassLambda); + + if (mcParticle.pdgCode() == 310 && mcParticle.isPhysicalPrimary() && TMath::Abs(ymc) < rapidityWindow) { + histos.fill(HIST("hStatisticsK0s"), 2.0f); // found + } + if (mcParticle.pdgCode() == 3122 && mcParticle.isPhysicalPrimary() && TMath::Abs(ymc) < rapidityWindow) { + histos.fill(HIST("hStatisticsLambda"), 2.0f); // found + } + if (mcParticle.pdgCode() == -3122 && mcParticle.isPhysicalPrimary() && TMath::Abs(ymc) < rapidityWindow) { + histos.fill(HIST("hStatisticsAlambda"), 2.0f); // found + } // if I got here, it means this MC particle was not recoed and is of interest. Add it please mcV0infos.push_back(thisInfo); @@ -271,8 +405,9 @@ struct lambdakzeromcbuilder { info.pdgCodeMother, info.pdgCodePositive, info.pdgCodeNegative, info.isPhysicalPrimary, info.xyz[0], info.xyz[1], info.xyz[2], info.posP[0], info.posP[1], info.posP[2], - info.negP[0], info.negP[1], info.negP[2]); - v0mccollref(thisInfo.mcCollision); + info.negP[0], info.negP[1], info.negP[2], + info.momentum[0], info.momentum[1], info.momentum[2]); + v0mccollref(info.mcCollision); } } diff --git a/PWGLF/TableProducer/Strangeness/strangederivedbuilder.cxx b/PWGLF/TableProducer/Strangeness/strangederivedbuilder.cxx index e6bfeee0ced..08526c22d5d 100644 --- a/PWGLF/TableProducer/Strangeness/strangederivedbuilder.cxx +++ b/PWGLF/TableProducer/Strangeness/strangederivedbuilder.cxx @@ -211,8 +211,8 @@ struct strangederivedbuilder { // Creation of histograms: MC generated for (Int_t i = 0; i < nSpecies; i++) { - histos.add(Form("hGen%s", particleNames[i].data()), Form("hGen%s", particleNames[i].data()), kTH1D, {axisPt}); - histos.add(Form("h2dGen%s", particleNames[i].data()), Form("h2dGen%s", particleNames[i].data()), kTH2D, {axisCentrality, axisPt}); + histos.add(Form("hGenerated%s", particleNames[i].data()), Form("hGenerated%s", particleNames[i].data()), kTH1D, {axisPt}); + histos.add(Form("h2dGenerated%s", particleNames[i].data()), Form("h2dGenerated%s", particleNames[i].data()), kTH2D, {axisCentrality, axisPt}); } histos.add("h2dNVerticesVsCentrality", "h2dNVerticesVsCentrality", kTH2D, {axisCentrality, axisNVertices}); @@ -237,7 +237,7 @@ struct strangederivedbuilder { if (doprocessBinnedGenerated) { // reserve space for generated vectors if that process enabled - auto hBinFinder = histos.get(HIST("h2dGenK0Short")); + auto hBinFinder = histos.get(HIST("h2dGeneratedK0Short")); LOGF(info, "Binned generated processing enabled. Initialising with %i elements...", hBinFinder->GetNcells()); genK0Short.resize(hBinFinder->GetNcells(), 0); genLambda.resize(hBinFinder->GetNcells(), 0); @@ -364,7 +364,7 @@ struct strangederivedbuilder { tracasccollref(TraCascadeCollIndices[casc.globalIndex()]); } - void processCollisionsMC(soa::Join const& collisions, soa::Join const& V0s, soa::Join const& /*V0MCCores*/, soa::Join const& Cascades, aod::KFCascDatas const& KFCascades, aod::TraCascDatas const& TraCascades, aod::BCsWithTimestamps const&, soa::Join const& mcCollisions, aod::McParticles const&) + void processCollisionsMC(soa::Join const& collisions, soa::Join const& V0s, soa::Join const& /*V0MCCores*/, soa::Join const& Cascades, aod::KFCascDatas const& KFCascades, aod::TraCascDatas const& TraCascades, aod::BCsWithTimestamps const&, soa::Join const& mcCollisions, aod::McParticles const&) { // create collision indices beforehand std::vector V0CollIndices(V0s.size(), -1); // index -1: no collision @@ -470,7 +470,7 @@ struct strangederivedbuilder { // Figure out the numbering of the new tracks table // assume filling per order int nTracks = 0; - for (int i = 0; i < trackMap.size(); i++) { + for (int i = 0; i < int(trackMap.size()); i++) { if (trackMap[i] >= 0) { trackMap[i] = nTracks++; } @@ -545,7 +545,7 @@ struct strangederivedbuilder { // Figure out the numbering of the new tracks table // assume filling per order int nTracks = 0; - for (int i = 0; i < trackMap.size(); i++) { + for (int i = 0; i < int(trackMap.size()); i++) { if (trackMap[i] >= 0) { trackMap[i] = nTracks++; } @@ -637,7 +637,7 @@ struct strangederivedbuilder { // Figure out the numbering of the new mcMother table // assume filling per order int nParticles = 0; - for (int i = 0; i < motherReference.size(); i++) { + for (int i = 0; i < int(motherReference.size()); i++) { if (motherReference[i] >= 0) { motherReference[i] = nParticles++; // count particles of interest } @@ -710,7 +710,7 @@ struct strangederivedbuilder { static_for<0, nSpecies - 1>([&](auto i) { constexpr int index = i.value; if (mcp.pdgCode() == particlePDGCodes[index] && bitcheck(enabledBits, index)) { - histos.fill(HIST("hGen") + HIST(particleNamesConstExpr[index]), mcp.pt()); + histos.fill(HIST("hGenerated") + HIST(particleNamesConstExpr[index]), mcp.pt()); } }); } @@ -740,7 +740,7 @@ struct strangederivedbuilder { static_for<0, nSpecies - 1>([&](auto i) { constexpr int index = i.value; if (mcp.pdgCode() == particlePDGCodes[index] && bitcheck(enabledBits, index)) { - histos.fill(HIST("h2dGen") + HIST(particleNamesConstExpr[index]), bestCentrality, mcp.pt()); + histos.fill(HIST("h2dGenerated") + HIST(particleNamesConstExpr[index]), bestCentrality, mcp.pt()); } }); } @@ -763,7 +763,7 @@ struct strangederivedbuilder { const uint64_t mcCollIndex = mcCollision.globalIndex(); // use one of the generated histograms as the bin finder - auto hBinFinder = histos.get(HIST("h2dGenK0Short")); + auto hBinFinder = histos.get(HIST("h2dGeneratedK0Short")); auto mcParticles = mcParticlesEntireTable.sliceBy(mcParticlePerMcCollision, mcCollIndex); for (auto& mcp : mcParticles) { diff --git a/PWGLF/Tasks/Strangeness/derivedlambdakzeroanalysis.cxx b/PWGLF/Tasks/Strangeness/derivedlambdakzeroanalysis.cxx index ac1835fa9ea..0c163e41a9f 100644 --- a/PWGLF/Tasks/Strangeness/derivedlambdakzeroanalysis.cxx +++ b/PWGLF/Tasks/Strangeness/derivedlambdakzeroanalysis.cxx @@ -58,7 +58,8 @@ using std::array; using dauTracks = soa::Join; using dauMCTracks = soa::Join; using v0Candidates = soa::Join; -using v0MCCandidates = soa::Join; +// using v0MCCandidates = soa::Join; +using v0MCCandidates = soa::Join; // simple checkers, but ensure 64 bit integers #define bitset(var, nbit) ((var) |= (static_cast(1) << static_cast(nbit))) @@ -80,11 +81,6 @@ struct derivedlambdakzeroanalysis { Configurable requireIsVertexTOFmatched{"requireIsVertexTOFmatched", false, "require events with at least one of vertex contributors matched to TOF"}; Configurable requireIsVertexTRDmatched{"requireIsVertexTRDmatched", false, "require events with at least one of vertex contributors matched to TRD"}; Configurable rejectSameBunchPileup{"rejectSameBunchPileup", true, "reject collisions in case of pileup with another collision in the same foundBC"}; - Configurable requireNoHighOccupancyAgressive{"requireNoHighOccupancyAgressive", false, "reject collisions with high occupancies according to the aggressive cuts"}; - Configurable requireNoHighOccupancyStrict{"requireNoHighOccupancyStrict", false, "reject collisions with high occupancies according to the strict cuts"}; - Configurable requireNoHighOccupancyMedium{"requireNoHighOccupancyMedium", false, "reject collisions with high occupancies according to the medium cuts"}; - Configurable requireNoHighOccupancyRelaxed{"requireNoHighOccupancyRelaxed", false, "reject collisions with high occupancies according to the relaxed cuts"}; - Configurable requireNoHighOccupancyGentle{"requireNoHighOccupancyGentle", false, "reject collisions with high occupancies according to the gentle cuts"}; Configurable requireNoCollInTimeRangeStd{"requireNoCollInTimeRangeStd", true, "reject collisions corrupted by the cannibalism, with other collisions within +/- 10 microseconds"}; Configurable requireNoCollInTimeRangeNarrow{"requireNoCollInTimeRangeNarrow", false, "reject collisions corrupted by the cannibalism, with other collisions within +/- 10 microseconds"}; @@ -180,6 +176,10 @@ struct derivedlambdakzeroanalysis { // MC coll assoc QA axis ConfigurableAxis axisMonteCarloNch{"axisMonteCarloNch", {300, 0.0f, 3000.0f}, "N_{ch} MC"}; + // For manual sliceBy + // Preslice> perMcCollision = aod::v0data::straMCCollisionId; + PresliceUnsorted> perMcCollision = aod::v0data::straMCCollisionId; + enum selection : uint64_t { selCosPA = 0, selRadius, selRadiusMax, @@ -322,15 +322,10 @@ struct derivedlambdakzeroanalysis { histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(8, "kIsVertexTOFmatched"); histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(9, "kIsVertexTRDmatched"); histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(10, "kNoSameBunchPileup"); - histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(11, "kNoHighOccupancyAgressive"); - histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(12, "kNoHighOccupancyStrict"); - histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(13, "kNoHighOccupancyMedium"); - histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(14, "kNoHighOccupancyRelaxed"); - histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(15, "kNoHighOccupancyGentle"); - histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(16, "kNoCollInTimeRangeStd"); - histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(17, "kNoCollInTimeRangeNarrow"); - histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(18, "Below min occup."); - histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(19, "Above max occup."); + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(11, "kNoCollInTimeRangeStd"); + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(12, "kNoCollInTimeRangeNarrow"); + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(13, "Below min occup."); + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(14, "Above max occup."); histos.add("hEventCentrality", "hEventCentrality", kTH1F, {{100, 0.0f, +100.0f}}); histos.add("hCentralityVsNch", "hCentralityVsNch", kTH2F, {axisCentrality, axisNch}); @@ -542,7 +537,17 @@ struct derivedlambdakzeroanalysis { histos.add("GeneralQA/h2dArmenterosSelected", "h2dArmenterosSelected", kTH2F, {axisAPAlpha, axisAPQt}); // Creation of histograms: MC generated - if (doprocessBinnedGenerated) { + if (doprocessGenerated) { + histos.add("hGenEvents", "hGenEvents", kTH2F, {{axisNch}, {2, -0.5f, +1.5f}}); + histos.get(HIST("hGenEvents"))->GetYaxis()->SetBinLabel(1, "All gen. events"); + histos.get(HIST("hGenEvents"))->GetYaxis()->SetBinLabel(2, "Gen. with at least 1 rec. events"); + histos.add("hGenEventCentrality", "hGenEventCentrality", kTH1F, {{100, 0.0f, +100.0f}}); + + histos.add("hCentralityVsNcoll_beforeEvSel", "hCentralityVsNcoll_beforeEvSel", kTH2F, {axisCentrality, {50, -0.5f, 49.5f}}); + histos.add("hCentralityVsNcoll_afterEvSel", "hCentralityVsNcoll_afterEvSel", kTH2F, {axisCentrality, {50, -0.5f, 49.5f}}); + + histos.add("hCentralityVsMultMC", "hCentralityVsMultMC", kTH2F, {{100, 0.0f, 100.0f}, axisNch}); + histos.add("h2dGenK0Short", "h2dGenK0Short", kTH2D, {axisCentrality, axisPt}); histos.add("h2dGenLambda", "h2dGenLambda", kTH2D, {axisCentrality, axisPt}); histos.add("h2dGenAntiLambda", "h2dGenAntiLambda", kTH2D, {axisCentrality, axisPt}); @@ -550,6 +555,23 @@ struct derivedlambdakzeroanalysis { histos.add("h2dGenXiPlus", "h2dGenXiPlus", kTH2D, {axisCentrality, axisPt}); histos.add("h2dGenOmegaMinus", "h2dGenOmegaMinus", kTH2D, {axisCentrality, axisPt}); histos.add("h2dGenOmegaPlus", "h2dGenOmegaPlus", kTH2D, {axisCentrality, axisPt}); + + histos.add("h2dGenK0ShortVsMultMC", "h2dGenK0ShortVsMultMC", kTH2D, {axisNch, axisPt}); + histos.add("h2dGenLambdaVsMultMC", "h2dGenLambdaVsMultMC", kTH2D, {axisNch, axisPt}); + histos.add("h2dGenAntiLambdaVsMultMC", "h2dGenAntiLambdaVsMultMC", kTH2D, {axisNch, axisPt}); + histos.add("h2dGenXiMinusVsMultMC", "h2dGenXiMinusVsMultMC", kTH2D, {axisNch, axisPt}); + histos.add("h2dGenXiPlusVsMultMC", "h2dGenXiPlusVsMultMC", kTH2D, {axisNch, axisPt}); + histos.add("h2dGenOmegaMinusVsMultMC", "h2dGenOmegaMinusVsMultMC", kTH2D, {axisNch, axisPt}); + histos.add("h2dGenOmegaPlusVsMultMC", "h2dGenOmegaPlusVsMultMC", kTH2D, {axisNch, axisPt}); + } + if (doprocessBinnedGenerated) { + histos.add("h2dGeneratedK0Short", "h2dGeneratedK0Short", kTH2D, {axisCentrality, axisPt}); + histos.add("h2dGeneratedLambda", "h2dGeneratedLambda", kTH2D, {axisCentrality, axisPt}); + histos.add("h2dGeneratedAntiLambda", "h2dGeneratedAntiLambda", kTH2D, {axisCentrality, axisPt}); + histos.add("h2dGeneratedXiMinus", "h2dGeneratedXiMinus", kTH2D, {axisCentrality, axisPt}); + histos.add("h2dGeneratedXiPlus", "h2dGeneratedXiPlus", kTH2D, {axisCentrality, axisPt}); + histos.add("h2dGeneratedOmegaMinus", "h2dGeneratedOmegaMinus", kTH2D, {axisCentrality, axisPt}); + histos.add("h2dGeneratedOmegaPlus", "h2dGeneratedOmegaPlus", kTH2D, {axisCentrality, axisPt}); } // inspect histogram sizes, please @@ -1103,49 +1125,24 @@ struct derivedlambdakzeroanalysis { } histos.fill(HIST("hEventSelection"), 9 /* Not at same bunch pile-up */); - if (requireNoHighOccupancyAgressive && !collision.selection_bit(o2::aod::evsel::kNoHighOccupancyAgressive)) { - return; - } - histos.fill(HIST("hEventSelection"), 10 /* No occupancy according to the aggressive cuts */); - - if (requireNoHighOccupancyStrict && !collision.selection_bit(o2::aod::evsel::kNoHighOccupancyStrict)) { - return; - } - histos.fill(HIST("hEventSelection"), 11 /* No occupancy according to the strict cuts */); - - if (requireNoHighOccupancyMedium && !collision.selection_bit(o2::aod::evsel::kNoHighOccupancyMedium)) { - return; - } - histos.fill(HIST("hEventSelection"), 12 /* No occupancy according to the medium cuts */); - - if (requireNoHighOccupancyRelaxed && !collision.selection_bit(o2::aod::evsel::kNoHighOccupancyRelaxed)) { - return; - } - histos.fill(HIST("hEventSelection"), 13 /* No occupancy according to the relaxed cuts */); - - if (requireNoHighOccupancyGentle && !collision.selection_bit(o2::aod::evsel::kNoHighOccupancyGentle)) { - return; - } - histos.fill(HIST("hEventSelection"), 14 /* No occupancy according to the gentle cuts */); - if (requireNoCollInTimeRangeStd && !collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard)) { return; } - histos.fill(HIST("hEventSelection"), 15 /* No other collision within +/- 10 microseconds */); + histos.fill(HIST("hEventSelection"), 10 /* No other collision within +/- 10 microseconds */); if (requireNoCollInTimeRangeNarrow && !collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeNarrow)) { return; } - histos.fill(HIST("hEventSelection"), 16 /* No other collision within +/- 4 microseconds */); + histos.fill(HIST("hEventSelection"), 11 /* No other collision within +/- 4 microseconds */); if (minOccupancy > 0 && collision.trackOccupancyInTimeRange() < minOccupancy) { return; } - histos.fill(HIST("hEventSelection"), 17 /* Below min occupancy */); + histos.fill(HIST("hEventSelection"), 12 /* Below min occupancy */); if (maxOccupancy > 0 && collision.trackOccupancyInTimeRange() > maxOccupancy) { return; } - histos.fill(HIST("hEventSelection"), 18 /* Above max occupancy */); + histos.fill(HIST("hEventSelection"), 13 /* Above max occupancy */); float centrality = collision.centFT0C(); if (qaCentrality) { @@ -1184,7 +1181,7 @@ struct derivedlambdakzeroanalysis { // ______________________________________________________ // Simulated processing (subscribes to MC information too) - void processMonteCarlo(soa::Join::iterator const& collision, v0MCCandidates const& fullV0s, dauTracks const&, aod::MotherMCParts const&, soa::Join const& /*mccollisions*/) + void processMonteCarlo(soa::Join::iterator const& collision, v0MCCandidates const& fullV0s, dauTracks const&, aod::MotherMCParts const&, soa::Join const& /*mccollisions*/, soa::Join const&) { histos.fill(HIST("hEventSelection"), 0. /* all collisions */); if (!collision.sel8()) { @@ -1232,49 +1229,24 @@ struct derivedlambdakzeroanalysis { } histos.fill(HIST("hEventSelection"), 9 /* Not at same bunch pile-up */); - if (requireNoHighOccupancyAgressive && !collision.selection_bit(o2::aod::evsel::kNoHighOccupancyAgressive)) { - return; - } - histos.fill(HIST("hEventSelection"), 10 /* No occupancy according to the aggressive cuts */); - - if (requireNoHighOccupancyStrict && !collision.selection_bit(o2::aod::evsel::kNoHighOccupancyStrict)) { - return; - } - histos.fill(HIST("hEventSelection"), 11 /* No occupancy according to the strict cuts */); - - if (requireNoHighOccupancyMedium && !collision.selection_bit(o2::aod::evsel::kNoHighOccupancyMedium)) { - return; - } - histos.fill(HIST("hEventSelection"), 12 /* No occupancy according to the medium cuts */); - - if (requireNoHighOccupancyRelaxed && !collision.selection_bit(o2::aod::evsel::kNoHighOccupancyRelaxed)) { - return; - } - histos.fill(HIST("hEventSelection"), 13 /* No occupancy according to the relaxed cuts */); - - if (requireNoHighOccupancyGentle && !collision.selection_bit(o2::aod::evsel::kNoHighOccupancyGentle)) { - return; - } - histos.fill(HIST("hEventSelection"), 14 /* No occupancy according to the gentle cuts */); - if (requireNoCollInTimeRangeStd && !collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard)) { return; } - histos.fill(HIST("hEventSelection"), 15 /* No other collision within +/- 10 microseconds */); + histos.fill(HIST("hEventSelection"), 10 /* No other collision within +/- 10 microseconds */); if (requireNoCollInTimeRangeNarrow && !collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeNarrow)) { return; } - histos.fill(HIST("hEventSelection"), 16 /* No other collision within +/- 4 microseconds */); + histos.fill(HIST("hEventSelection"), 11 /* No other collision within +/- 4 microseconds */); if (minOccupancy > 0 && collision.trackOccupancyInTimeRange() < minOccupancy) { return; } - histos.fill(HIST("hEventSelection"), 17 /* Below min occupancy */); + histos.fill(HIST("hEventSelection"), 12 /* Below min occupancy */); if (maxOccupancy > 0 && collision.trackOccupancyInTimeRange() > maxOccupancy) { return; } - histos.fill(HIST("hEventSelection"), 18 /* Above max occupancy */); + histos.fill(HIST("hEventSelection"), 13 /* Above max occupancy */); float centrality = collision.centFT0C(); if (qaCentrality) { @@ -1295,18 +1267,23 @@ struct derivedlambdakzeroanalysis { if (std::abs(v0.negativeeta()) > daughterEtaCut || std::abs(v0.positiveeta()) > daughterEtaCut) continue; // remove acceptance that's badly reproduced by MC / superfluous in future + if (!v0.has_v0MCCore()) + continue; + + auto v0MC = v0.v0MCCore_as>(); + // fill AP plot for all V0s histos.fill(HIST("GeneralQA/h2dArmenterosAll"), v0.alpha(), v0.qtarm()); - float ptmc = RecoDecay::sqrtSumOfSquares(v0.pxPosMC() + v0.pxNegMC(), v0.pyPosMC() + v0.pyNegMC()); + float ptmc = RecoDecay::sqrtSumOfSquares(v0MC.pxPosMC() + v0MC.pxNegMC(), v0MC.pyPosMC() + v0MC.pyNegMC()); float ymc = 1e-3; - if (v0.pdgCode() == 310) - ymc = RecoDecay::y(std::array{v0.pxPosMC() + v0.pxNegMC(), v0.pyPosMC() + v0.pyNegMC(), v0.pzPosMC() + v0.pzNegMC()}, o2::constants::physics::MassKaonNeutral); - else if (TMath::Abs(v0.pdgCode()) == 3122) - ymc = RecoDecay::y(std::array{v0.pxPosMC() + v0.pxNegMC(), v0.pyPosMC() + v0.pyNegMC(), v0.pzPosMC() + v0.pzNegMC()}, o2::constants::physics::MassLambda); + if (v0MC.pdgCode() == 310) + ymc = RecoDecay::y(std::array{v0MC.pxPosMC() + v0MC.pxNegMC(), v0MC.pyPosMC() + v0MC.pyNegMC(), v0MC.pzPosMC() + v0MC.pzNegMC()}, o2::constants::physics::MassKaonNeutral); + else if (TMath::Abs(v0MC.pdgCode()) == 3122) + ymc = RecoDecay::y(std::array{v0MC.pxPosMC() + v0MC.pxNegMC(), v0MC.pyPosMC() + v0MC.pyNegMC(), v0MC.pzPosMC() + v0MC.pzNegMC()}, o2::constants::physics::MassLambda); uint64_t selMap = computeReconstructionBitmap(v0, collision, ymc, ymc, ptmc); - selMap = selMap | computeMCAssociation(v0); + selMap = selMap | computeMCAssociation(v0MC); // feeddown matrix always with association if (calculateFeeddownMatrix) @@ -1327,7 +1304,7 @@ struct derivedlambdakzeroanalysis { if (collision.has_straMCCollision()) { auto mcCollision = collision.straMCCollision_as>(); mcNch = mcCollision.multMCNParticlesEta05(); - correctCollision = (v0.straMCCollisionId() == mcCollision.globalIndex()); + correctCollision = (v0MC.straMCCollisionId() == mcCollision.globalIndex()); } analyseCollisionAssociation(v0, ptmc, mcNch, correctCollision, selMap); } @@ -1335,6 +1312,176 @@ struct derivedlambdakzeroanalysis { } // end v0 loop } + // ______________________________________________________ + // Simulated processing (subscribes to MC information too) + void processGenerated(soa::Join const& mcCollisions, soa::Join const& V0MCCores, soa::Join const& CascMCCores, soa::Join const& collisions) + { + std::vector listBestCollisionIdx = fillGenEventHist(mcCollisions, collisions); + for (auto const& v0MC : V0MCCores) { + if (!v0MC.has_straMCCollision()) + continue; + + if (!v0MC.isPhysicalPrimary()) + continue; + + float ptmc = RecoDecay::sqrtSumOfSquares(v0MC.pxPosMC() + v0MC.pxNegMC(), v0MC.pyPosMC() + v0MC.pyNegMC()); + float ymc = 1e3; + if (v0MC.pdgCode() == 310) + ymc = RecoDecay::y(std::array{v0MC.pxPosMC() + v0MC.pxNegMC(), v0MC.pyPosMC() + v0MC.pyNegMC(), v0MC.pzPosMC() + v0MC.pzNegMC()}, o2::constants::physics::MassKaonNeutral); + else if (TMath::Abs(v0MC.pdgCode()) == 3122) + ymc = RecoDecay::y(std::array{v0MC.pxPosMC() + v0MC.pxNegMC(), v0MC.pyPosMC() + v0MC.pyNegMC(), v0MC.pzPosMC() + v0MC.pzNegMC()}, o2::constants::physics::MassLambda); + + if (TMath::Abs(ymc) > rapidityCut) + continue; + + auto mcCollision = v0MC.straMCCollision_as>(); + float centrality = 100.5f; + if (listBestCollisionIdx[mcCollision.globalIndex()] > -1) { + auto collision = collisions.iteratorAt(listBestCollisionIdx[mcCollision.globalIndex()]); + centrality = collision.centFT0C(); + } + + if (v0MC.pdgCode() == 310) { + histos.fill(HIST("h2dGenK0Short"), centrality, ptmc); + histos.fill(HIST("h2dGenK0ShortVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc); + } + if (v0MC.pdgCode() == 3122) { + histos.fill(HIST("h2dGenLambda"), centrality, ptmc); + histos.fill(HIST("h2dGenLambdaVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc); + } + if (v0MC.pdgCode() == -3122) { + histos.fill(HIST("h2dGenAntiLambda"), centrality, ptmc); + histos.fill(HIST("h2dGenAntiLambdaVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc); + } + } + + for (auto const& cascMC : CascMCCores) { + if (!cascMC.has_straMCCollision()) + continue; + + if (!cascMC.isPhysicalPrimary()) + continue; + + float ptmc = RecoDecay::sqrtSumOfSquares(cascMC.pxMC(), cascMC.pyMC()); + float ymc = 1e3; + if (TMath::Abs(cascMC.pdgCode()) == 3312) + ymc = RecoDecay::y(std::array{cascMC.pxMC(), cascMC.pyMC(), cascMC.pzMC()}, o2::constants::physics::MassXiMinus); + else if (TMath::Abs(cascMC.pdgCode()) == 3334) + ymc = RecoDecay::y(std::array{cascMC.pxMC(), cascMC.pyMC(), cascMC.pzMC()}, o2::constants::physics::MassOmegaMinus); + + if (TMath::Abs(ymc) > rapidityCut) + continue; + + auto mcCollision = cascMC.straMCCollision_as>(); + float centrality = 100.5f; + if (listBestCollisionIdx[mcCollision.globalIndex()] > -1) { + auto collision = collisions.iteratorAt(listBestCollisionIdx[mcCollision.globalIndex()]); + centrality = collision.centFT0C(); + } + + if (cascMC.pdgCode() == 3312) { + histos.fill(HIST("h2dGenXiMinus"), centrality, ptmc); + histos.fill(HIST("h2dGenXiMinusVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc); + } + if (cascMC.pdgCode() == -3312) { + histos.fill(HIST("h2dGenXiPlus"), centrality, ptmc); + histos.fill(HIST("h2dGenXiPlusVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc); + } + if (cascMC.pdgCode() == 3334) { + histos.fill(HIST("h2dGenOmegaMinus"), centrality, ptmc); + histos.fill(HIST("h2dGenOmegaMinusVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc); + } + if (cascMC.pdgCode() == -3334) { + histos.fill(HIST("h2dGenOmegaPlus"), centrality, ptmc); + histos.fill(HIST("h2dGenOmegaPlusVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc); + } + } + } + + // ______________________________________________________ + // Simulated processing + // Fill event information (for event loss estimation) and return the index to the recoed collision associated to a given MC collision. + std::vector fillGenEventHist(soa::Join const& mcCollisions, soa::Join const& collisions) + { + std::vector listBestCollisionIdx(mcCollisions.size()); + for (auto const& mcCollision : mcCollisions) { + histos.fill(HIST("hGenEvents"), mcCollision.multMCNParticlesEta05(), 0 /* all gen. events*/); + + auto groupedCollisions = collisions.sliceBy(perMcCollision, mcCollision.globalIndex()); + // Check if there is at least one of the reconstructed collisions associated to this MC collision + // If so, we consider it + bool atLeastOne = false; + int biggestNContribs = -1; + int bestCollisionIndex = -1; + float centrality = 100.5f; + int nCollisions = 0; + for (auto const& collision : groupedCollisions) { + if (!collision.sel8()) { + continue; + } + if (std::abs(collision.posZ()) > 10.f) { + continue; + } + if (rejectITSROFBorder && !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) { + continue; + } + if (rejectTFBorder && !collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) { + continue; + } + if (requireIsVertexITSTPC && !collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) { + continue; + } + if (requireIsGoodZvtxFT0VsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { + continue; + } + if (requireIsVertexTOFmatched && !collision.selection_bit(o2::aod::evsel::kIsVertexTOFmatched)) { + continue; + } + if (requireIsVertexTRDmatched && !collision.selection_bit(o2::aod::evsel::kIsVertexTRDmatched)) { + continue; + } + if (rejectSameBunchPileup && !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { + continue; + } + if (requireNoCollInTimeRangeStd && !collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard)) { + continue; + } + if (requireNoCollInTimeRangeNarrow && !collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeNarrow)) { + continue; + } + + if (minOccupancy > 0 && collision.trackOccupancyInTimeRange() < minOccupancy) { + continue; + } + if (maxOccupancy > 0 && collision.trackOccupancyInTimeRange() > maxOccupancy) { + continue; + } + + if (biggestNContribs < collision.multPVTotalContributors()) { + biggestNContribs = collision.multPVTotalContributors(); + bestCollisionIndex = collision.globalIndex(); + centrality = collision.centFT0C(); + } + nCollisions++; + + atLeastOne = true; + } + listBestCollisionIdx[mcCollision.globalIndex()] = bestCollisionIndex; + + histos.fill(HIST("hCentralityVsNcoll_beforeEvSel"), centrality, groupedCollisions.size()); + histos.fill(HIST("hCentralityVsNcoll_afterEvSel"), centrality, nCollisions); + + histos.fill(HIST("hCentralityVsMultMC"), centrality, mcCollision.multMCNParticlesEta05()); + + if (atLeastOne) { + histos.fill(HIST("hGenEvents"), mcCollision.multMCNParticlesEta05(), 1 /* at least 1 rec. event*/); + + histos.fill(HIST("hGenEventCentrality"), centrality); + } + } + return listBestCollisionIdx; + } + // ______________________________________________________ // Simulated processing (subscribes to MC information too) void processBinnedGenerated( @@ -1342,59 +1489,59 @@ struct derivedlambdakzeroanalysis { aod::GeXiMinus const& geXiMinus, aod::GeXiPlus const& geXiPlus, aod::GeOmegaMinus const& geOmegaMinus, aod::GeOmegaPlus const& geOmegaPlus) { - auto hK0Short = histos.get(HIST("h2dGenK0Short")); - auto hLambda = histos.get(HIST("h2dGenLambda")); - auto hAntiLambda = histos.get(HIST("h2dGenAntiLambda")); - auto hXiMinus = histos.get(HIST("h2dGenXiMinus")); - auto hXiPlus = histos.get(HIST("h2dGenXiPlus")); - auto hOmegaMinus = histos.get(HIST("h2dGenOmegaMinus")); - auto hOmegaPlus = histos.get(HIST("h2dGenOmegaPlus")); + auto hK0Short = histos.get(HIST("h2dGeneratedK0Short")); + auto hLambda = histos.get(HIST("h2dGeneratedLambda")); + auto hAntiLambda = histos.get(HIST("h2dGeneratedAntiLambda")); + auto hXiMinus = histos.get(HIST("h2dGeneratedXiMinus")); + auto hXiPlus = histos.get(HIST("h2dGeneratedXiPlus")); + auto hOmegaMinus = histos.get(HIST("h2dGeneratedOmegaMinus")); + auto hOmegaPlus = histos.get(HIST("h2dGeneratedOmegaPlus")); for (auto& gVec : geK0Short) { - if (gVec.generatedK0Short().size() != hK0Short->GetNcells()) + if (int(gVec.generatedK0Short().size()) != hK0Short->GetNcells()) LOGF(fatal, "K0Short: Number of elements in generated array and number of cells in receiving histogram differ: %i vs %i!", gVec.generatedK0Short().size(), hK0Short->GetNcells()); - for (uint64_t iv = 0; iv < hK0Short->GetNcells(); iv++) { + for (int iv = 0; iv < hK0Short->GetNcells(); iv++) { hK0Short->SetBinContent(iv, hK0Short->GetBinContent(iv) + gVec.generatedK0Short()[iv]); } } for (auto& gVec : geLambda) { - if (gVec.generatedLambda().size() != hLambda->GetNcells()) + if (int(gVec.generatedLambda().size()) != hLambda->GetNcells()) LOGF(fatal, "Lambda: Number of elements in generated array and number of cells in receiving histogram differ: %i vs %i!", gVec.generatedLambda().size(), hLambda->GetNcells()); - for (uint64_t iv = 0; iv < hLambda->GetNcells(); iv++) { + for (int iv = 0; iv < hLambda->GetNcells(); iv++) { hLambda->SetBinContent(iv, hLambda->GetBinContent(iv) + gVec.generatedLambda()[iv]); } } for (auto& gVec : geAntiLambda) { - if (gVec.generatedAntiLambda().size() != hAntiLambda->GetNcells()) + if (int(gVec.generatedAntiLambda().size()) != hAntiLambda->GetNcells()) LOGF(fatal, "AntiLambda: Number of elements in generated array and number of cells in receiving histogram differ: %i vs %i!", gVec.generatedAntiLambda().size(), hAntiLambda->GetNcells()); - for (uint64_t iv = 0; iv < hAntiLambda->GetNcells(); iv++) { + for (int iv = 0; iv < hAntiLambda->GetNcells(); iv++) { hAntiLambda->SetBinContent(iv, hAntiLambda->GetBinContent(iv) + gVec.generatedAntiLambda()[iv]); } } for (auto& gVec : geXiMinus) { - if (gVec.generatedXiMinus().size() != hXiMinus->GetNcells()) + if (int(gVec.generatedXiMinus().size()) != hXiMinus->GetNcells()) LOGF(fatal, "XiMinus: Number of elements in generated array and number of cells in receiving histogram differ: %i vs %i!", gVec.generatedXiMinus().size(), hXiMinus->GetNcells()); - for (uint64_t iv = 0; iv < hXiMinus->GetNcells(); iv++) { + for (int iv = 0; iv < hXiMinus->GetNcells(); iv++) { hXiMinus->SetBinContent(iv, hXiMinus->GetBinContent(iv) + gVec.generatedXiMinus()[iv]); } } for (auto& gVec : geXiPlus) { - if (gVec.generatedXiPlus().size() != hXiPlus->GetNcells()) + if (int(gVec.generatedXiPlus().size()) != hXiPlus->GetNcells()) LOGF(fatal, "XiPlus: Number of elements in generated array and number of cells in receiving histogram differ: %i vs %i!", gVec.generatedXiPlus().size(), hXiPlus->GetNcells()); - for (uint64_t iv = 0; iv < hXiPlus->GetNcells(); iv++) { + for (int iv = 0; iv < hXiPlus->GetNcells(); iv++) { hXiPlus->SetBinContent(iv, hXiPlus->GetBinContent(iv) + gVec.generatedXiPlus()[iv]); } } for (auto& gVec : geOmegaMinus) { - if (gVec.generatedOmegaMinus().size() != hOmegaMinus->GetNcells()) + if (int(gVec.generatedOmegaMinus().size()) != hOmegaMinus->GetNcells()) LOGF(fatal, "OmegaMinus: Number of elements in generated array and number of cells in receiving histogram differ: %i vs %i!", gVec.generatedOmegaMinus().size(), hOmegaMinus->GetNcells()); - for (uint64_t iv = 0; iv < hOmegaMinus->GetNcells(); iv++) { + for (int iv = 0; iv < hOmegaMinus->GetNcells(); iv++) { hOmegaMinus->SetBinContent(iv, hOmegaMinus->GetBinContent(iv) + gVec.generatedOmegaMinus()[iv]); } } for (auto& gVec : geOmegaPlus) { - if (gVec.generatedOmegaPlus().size() != hOmegaPlus->GetNcells()) + if (int(gVec.generatedOmegaPlus().size()) != hOmegaPlus->GetNcells()) LOGF(fatal, "OmegaPlus: Number of elements in generated array and number of cells in receiving histogram differ: %i vs %i!", gVec.generatedOmegaPlus().size(), hOmegaPlus->GetNcells()); - for (uint64_t iv = 0; iv < hOmegaPlus->GetNcells(); iv++) { + for (int iv = 0; iv < hOmegaPlus->GetNcells(); iv++) { hOmegaPlus->SetBinContent(iv, hOmegaPlus->GetBinContent(iv) + gVec.generatedOmegaPlus()[iv]); } } @@ -1403,6 +1550,7 @@ struct derivedlambdakzeroanalysis { PROCESS_SWITCH(derivedlambdakzeroanalysis, processRealData, "process as if real data", true); PROCESS_SWITCH(derivedlambdakzeroanalysis, processMonteCarlo, "process as if MC", false); PROCESS_SWITCH(derivedlambdakzeroanalysis, processBinnedGenerated, "process MC generated", false); + PROCESS_SWITCH(derivedlambdakzeroanalysis, processGenerated, "process MC generated", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)