diff --git a/PWGLF/TableProducer/Strangeness/strangenessbuilder.cxx b/PWGLF/TableProducer/Strangeness/strangenessbuilder.cxx index 31a015e28e2..151f3fe3daa 100644 --- a/PWGLF/TableProducer/Strangeness/strangenessbuilder.cxx +++ b/PWGLF/TableProducer/Strangeness/strangenessbuilder.cxx @@ -9,14 +9,14 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. // -// Strangeness builder task +// Strangeness builder task // ======================== -// -// This task produces all tables that may be necessary for -// strangeness analyses. A single device is provided to -// ensure better computing resource (memory) management. // -// process functions: +// This task produces all tables that may be necessary for +// strangeness analyses. A single device is provided to +// ensure better computing resource (memory) management. +// +// process functions: // -- processPreselectTPCPID ..: pre-selects TPC dE/dx-compatible candidates. // -- processRealData .........: use this OR processSimulation but NOT both // -- processSimulation .......: use this OR processRealData but NOT both @@ -38,94 +38,160 @@ using namespace o2; using namespace o2::framework; static constexpr int nParameters = 1; -static const std::vector tableNames{"V0Indices", //.0 (standard analysis: V0Data) - "V0CoresBase", //.1 (standard analyses: V0Data) - "V0Covs", //.2 - "CascIndices", //.3 (standard analyses: CascData) - "KFCascIndices", //.4 (standard analyses: KFCascData) - "TraCascIndices", //.5 (standard analyses: TraCascData) - "StoredCascCores", //.6 (standard analyses: CascData) - "StoredKFCascCores", //.7 (standard analyses: KFCascData) - "StoredTraCascCores", //.8 (standard analyses: TraCascData) - "CascCovs", //.9 - "KFCascCovs", //.10 - "TraCascCovs", //.11 - "V0TrackXs", //.12 - "CascTrackXs", //.13 - "CascBBs", //.14 - "V0DauCovs", //.15 (requested: tracking studies) - "V0DauCovIUs", //.16 (requested: tracking studies) - "V0TraPosAtDCAs", //.17 (requested: tracking studies) - "V0TraPosAtIUs", //.18 (requested: tracking studies) - "V0Ivanovs", //.19 (requested: tracking studies) - "McV0Labels", //.20 (MC/standard analysis) - "V0MCCores", //.21 (MC) - "V0CoreMCLabels", //.22 (MC) - "V0MCCollRefs", //.23 (MC) - "McCascLabels", // 24 (MC/standard analysis) - "McKFCascLabels", // 25 (MC) - "McTraCascLabels", // 26 (MC) - "McCascBBTags", // 27 (MC) - "CascMCCores", // 28 (MC) - "CascCoreMCLabels", // 29 (MC) - "CascMCCollRefs", // 30 (MC) - "StraCollision", // 31 (derived) - "StraCollLabels", // 32 (derived) - "StraMCCollisions", // 33 (MC/derived) - "StraMCCollMults", // 34 (MC/derived) - "StraCents", // 35 (derived) - "StraEvSels", // 36 (derived) - "StraStamps", // 37 (derived) - "V0CollRefs", // 38 (derived) - "CascCollRefs", // 39 (derived) - "KFCascCollRefs", // 40 (derived) - "TraCascCollRefs", // 41 (derived) - "DauTrackExtras", // 42 (derived) - "DauTrackMCIds", // 43 (MC/derived) - "DauTrackTPCPIDs", // 44 (derived) - "DauTrackTOFPIDs", // 45 (derived) - "V0Extras", // 46 (derived) - "CascExtras", // 47 (derived) - "StraTrackExtras", // 48 (derived) - "CascToTraRefs", //.49 (interlink) - "CascToKFRefs", //.50 (interlink) - "TraToCascRefs", //.51 (interlink) - "KFToCascRefs", //.52 (interlink) - "V0MCMothers", // 53 (MC/derived) - "CascMCMothers", // 54 (MC/derived) - "MotherMCParts", // 55 (MC/derived) - "StraFT0AQVs", // 56 (derived) - "StraFT0CQVs", // 57 (derived) - "StraFT0MQVs", // 58 (derived) - "StraFV0AQVs", // 59 (derived) - "StraTPCQVs", // 60 (derived) - "StraFT0CQVsEv", // 61 (derived) - "StraZDCSP", // 62 (derived) - "GeK0Short", // 63 (MC/derived) - "GeLambda", // 64 (MC/derived) - "GeAntiLambda", // 65 (MC/derived) - "GeXiMinus", // 66 (MC/derived) - "GeXiPlus", // 67 (MC/derived) - "GeOmegaMinus", // 68 (MC/derived) - "GeOmegaPlus", // 69 (MC/derived) - "V0FoundTags", // 70 (MC/derived) - "CascFoundTags", // 71 (MC/derived) - "StraOrigins" // 72 (derived) - }; +static const std::vector tableNames{ + "V0Indices", //.0 (standard analysis: V0Data) + "V0CoresBase", //.1 (standard analyses: V0Data) + "V0Covs", //.2 + "CascIndices", //.3 (standard analyses: CascData) + "KFCascIndices", //.4 (standard analyses: KFCascData) + "TraCascIndices", //.5 (standard analyses: TraCascData) + "StoredCascCores", //.6 (standard analyses: CascData) + "StoredKFCascCores", //.7 (standard analyses: KFCascData) + "StoredTraCascCores", //.8 (standard analyses: TraCascData) + "CascCovs", //.9 + "KFCascCovs", //.10 + "TraCascCovs", //.11 + "V0TrackXs", //.12 + "CascTrackXs", //.13 + "CascBBs", //.14 + "V0DauCovs", //.15 (requested: tracking studies) + "V0DauCovIUs", //.16 (requested: tracking studies) + "V0TraPosAtDCAs", //.17 (requested: tracking studies) + "V0TraPosAtIUs", //.18 (requested: tracking studies) + "V0Ivanovs", //.19 (requested: tracking studies) + "McV0Labels", //.20 (MC/standard analysis) + "V0MCCores", //.21 (MC) + "V0CoreMCLabels", //.22 (MC) + "V0MCCollRefs", //.23 (MC) + "McCascLabels", // 24 (MC/standard analysis) + "McKFCascLabels", // 25 (MC) + "McTraCascLabels", // 26 (MC) + "McCascBBTags", // 27 (MC) + "CascMCCores", // 28 (MC) + "CascCoreMCLabels", // 29 (MC) + "CascMCCollRefs", // 30 (MC) + "StraCollision", // 31 (derived) + "StraCollLabels", // 32 (derived) + "StraMCCollisions", // 33 (MC/derived) + "StraMCCollMults", // 34 (MC/derived) + "StraCents", // 35 (derived) + "StraEvSels", // 36 (derived) + "StraStamps", // 37 (derived) + "V0CollRefs", // 38 (derived) + "CascCollRefs", // 39 (derived) + "KFCascCollRefs", // 40 (derived) + "TraCascCollRefs", // 41 (derived) + "DauTrackExtras", // 42 (derived) + "DauTrackMCIds", // 43 (MC/derived) + "DauTrackTPCPIDs", // 44 (derived) + "DauTrackTOFPIDs", // 45 (derived) + "V0Extras", // 46 (derived) + "CascExtras", // 47 (derived) + "StraTrackExtras", // 48 (derived) + "CascToTraRefs", //.49 (interlink) + "CascToKFRefs", //.50 (interlink) + "TraToCascRefs", //.51 (interlink) + "KFToCascRefs", //.52 (interlink) + "V0MCMothers", // 53 (MC/derived) + "CascMCMothers", // 54 (MC/derived) + "MotherMCParts", // 55 (MC/derived) + "StraFT0AQVs", // 56 (derived) + "StraFT0CQVs", // 57 (derived) + "StraFT0MQVs", // 58 (derived) + "StraFV0AQVs", // 59 (derived) + "StraTPCQVs", // 60 (derived) + "StraFT0CQVsEv", // 61 (derived) + "StraZDCSP", // 62 (derived) + "GeK0Short", // 63 (MC/derived) + "GeLambda", // 64 (MC/derived) + "GeAntiLambda", // 65 (MC/derived) + "GeXiMinus", // 66 (MC/derived) + "GeXiPlus", // 67 (MC/derived) + "GeOmegaMinus", // 68 (MC/derived) + "GeOmegaPlus", // 69 (MC/derived) + "V0FoundTags", // 70 (MC/derived) + "CascFoundTags", // 71 (MC/derived) + "StraOrigins" // 72 (derived) +}; static constexpr int nTablesConst = 73; static const std::vector parameterNames{"enable"}; static const int defaultParameters[nTablesConst][nParameters]{ - {-1}, {1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, //0-9 - {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, //10-19 - {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, //20-29 - {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, //30-39 - {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, //40-49 - {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, //50-59 - {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, {-1}, //60-69 - {-1}, {-1}, {-1} //70-72 - }; + {-1}, + {1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, // 0-9 + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, // 10-19 + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, // 20-29 + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, // 30-39 + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, // 40-49 + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, // 50-59 + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, // 60-69 + {-1}, + {-1}, + {-1} // 70-72 +}; // use parameters + cov mat non-propagated, aux info + (extension propagated) using FullTracksExt = soa::Join; @@ -139,7 +205,7 @@ using TracksExtraWithPID = soa::Join v0indices; // standard part of V0Datas - Produces v0cores; // standard part of V0Datas - Produces v0covs; // for decay chain reco + Produces v0indices; // standard part of V0Datas + Produces v0cores; // standard part of V0Datas + Produces v0covs; // for decay chain reco //__________________________________________________ // cascade tables - Produces cascidx; // standard part of CascDatas - Produces kfcascidx; // standard part of KFCascDatas - Produces tracascidx; // standard part of TraCascDatas - Produces cascdata; // standard part of CascDatas - Produces kfcascdata; // standard part of KFCascDatas - Produces tracascdata; // standard part of TraCascDatas - Produces casccovs; // for decay chain reco - Produces kfcasccovs; // for decay chain reco - Produces tracasccovs; // for decay chain reco + Produces cascidx; // standard part of CascDatas + Produces kfcascidx; // standard part of KFCascDatas + Produces tracascidx; // standard part of TraCascDatas + Produces cascdata; // standard part of CascDatas + Produces kfcascdata; // standard part of KFCascDatas + Produces tracascdata; // standard part of TraCascDatas + Produces casccovs; // for decay chain reco + Produces kfcasccovs; // for decay chain reco + Produces tracasccovs; // for decay chain reco //__________________________________________________ // interlink tables @@ -244,11 +310,11 @@ struct StrangenessBuilder { //__________________________________________________ // secondary auxiliary tables - Produces v0trackXs; // for decay chain reco - Produces cascTrackXs; // for decay chain reco + Produces v0trackXs; // for decay chain reco + Produces cascTrackXs; // for decay chain reco // further auxiliary / optional if desired - Produces cascbb; + Produces cascbb; Produces v0daucovs; // covariances of daughter tracks Produces v0daucovIUs; // covariances of daughter tracks Produces v0dauPositions; // auxiliary debug information @@ -261,17 +327,17 @@ struct StrangenessBuilder { Produces v0mccores; // mc info storage Produces v0CoreMCLabels; // interlink V0Cores -> V0MCCores Produces v0mccollref; // references collisions from V0MCCores - + // MC information: Cascades - Produces casclabels; // MC labels for cascades - Produces kfcasclabels; // MC labels for KF cascades - Produces tracasclabels; // MC labels for tracked cascades - Produces bbtags; // bb tags (inv structure tagging in mc) - Produces cascmccores; // mc info storage - Produces cascCoreMClabels; // interlink CascCores -> CascMCCores - Produces cascmccollrefs; // references MC collisions from MC cascades - -//__________________________________________________ + Produces casclabels; // MC labels for cascades + Produces kfcasclabels; // MC labels for KF cascades + Produces tracasclabels; // MC labels for tracked cascades + Produces bbtags; // bb tags (inv structure tagging in mc) + Produces cascmccores; // mc info storage + Produces cascCoreMClabels; // interlink CascCores -> CascMCCores + Produces cascmccollrefs; // references MC collisions from MC cascades + + //__________________________________________________ // cascade interlinks Produces cascToTraRefs; // cascades -> tracked Produces cascToKFRefs; // cascades -> KF @@ -358,7 +424,7 @@ struct StrangenessBuilder { std::string prefix = "v0BuilderOpts"; Configurable generatePhotonCandidates{"generatePhotonCandidates", false, "generate gamma conversion candidates (V0s using TPC-only tracks)"}; - // MC builder options + // MC builder options Configurable mc_populateV0MCCoresSymmetric{"mc_populateV0MCCoresSymmetric", false, "populate V0MCCores table for derived data analysis, keep V0MCCores joinable with V0Cores"}; Configurable mc_populateV0MCCoresAsymmetric{"mc_populateV0MCCoresAsymmetric", true, "populate V0MCCores table for derived data analysis, create V0Cores -> V0MCCores interlink. Saves only labeled V0s."}; Configurable mc_treatPiToMuDecays{"mc_treatPiToMuDecays", true, "if true, will correctly capture pi -> mu and V0 label will still point to originating V0 decay in those cases. Nota bene: prong info will still be for the muon!"}; @@ -382,7 +448,7 @@ struct StrangenessBuilder { 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"}; - // MC builder options + // MC builder options Configurable mc_populateCascMCCoresSymmetric{"mc_populateCascMCCoresSymmetric", false, "populate CascMCCores table for derived data analysis, keep CascMCCores joinable with CascCores"}; Configurable mc_populateCascMCCoresAsymmetric{"mc_populateCascMCCoresAsymmetric", true, "populate CascMCCores table for derived data analysis, create CascCores -> CascMCCores interlink. Saves only labeled Cascades."}; Configurable mc_addGeneratedXiMinus{"mc_addGeneratedXiMinus", false, "add CascMCCore entry for generated, not-recoed XiMinus"}; @@ -403,16 +469,16 @@ struct StrangenessBuilder { std::vector v0sFromCascades; // Vector of v0 candidates used in cascades std::vector v0Map; // index to relate V0s -> v0sFromCascades - // for establishing CascData/KFData/TraCascData interlinks + // for establishing CascData/KFData/TraCascData interlinks - struct{ + struct { std::vector cascCoreToCascades; std::vector kfCascCoreToCascades; std::vector traCascCoreToCascades; std::vector cascadeToCascCores; std::vector cascadeToKFCascCores; std::vector cascadeToTraCascCores; - } interlinks; + } interlinks; //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* // Helper struct to contain V0MCCore information prior to filling @@ -467,7 +533,7 @@ struct StrangenessBuilder { void init(InitContext& context) { mRunNumber = 0; - + mEnabledTables.resize(nTables, 0); LOGF(info, "Configuring tables to generate"); @@ -478,7 +544,7 @@ struct StrangenessBuilder { if (f == 1) { mEnabledTables[i] = 1; } - if( f == -1){ + if (f == -1) { // autodetect this table in other devices for (DeviceSpec const& device : workflows.devices) { // Step 1: check if this device subscribed to the V0data table @@ -500,7 +566,7 @@ struct StrangenessBuilder { LOGF(info, "*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*"); for (int i = 0; i < nTables; i++) { // printout to be improved in the future - if(mEnabledTables[i]){ + if (mEnabledTables[i]) { LOGF(info, "Table enabled: %s", tableNames[i]); } } @@ -554,14 +620,15 @@ struct StrangenessBuilder { straHelper.lut = lut; LOG(info) << "Fully configured for run: " << bc.runNumber(); - // mmark this run as configured + // mmark this run as configured mRunNumber = bc.runNumber(); return true; } //__________________________________________________ - void resetInterlinks(){ + void resetInterlinks() + { interlinks.cascCoreToCascades.clear(); interlinks.kfCascCoreToCascades.clear(); interlinks.traCascCoreToCascades.clear(); @@ -571,23 +638,24 @@ struct StrangenessBuilder { } //__________________________________________________ - void populateCascadeInterlinks(){ - if(mEnabledTables[kCascToKFRefs]){ + void populateCascadeInterlinks() + { + if (mEnabledTables[kCascToKFRefs]) { for (auto& cascCore : interlinks.cascCoreToCascades) { cascToKFRefs(interlinks.cascadeToKFCascCores[cascCore]); } } - if(mEnabledTables[kCascToTraRefs]){ + if (mEnabledTables[kCascToTraRefs]) { for (auto& cascCore : interlinks.cascCoreToCascades) { cascToTraRefs(interlinks.cascadeToTraCascCores[cascCore]); } } - if(mEnabledTables[kKFToCascRefs]){ + if (mEnabledTables[kKFToCascRefs]) { for (auto& kfCascCore : interlinks.kfCascCoreToCascades) { kfToCascRefs(interlinks.cascadeToCascCores[kfCascCore]); } } - if(mEnabledTables[kTraToCascRefs]){ + if (mEnabledTables[kTraToCascRefs]) { for (auto& traCascCore : interlinks.traCascCoreToCascades) { traToCascRefs(interlinks.cascadeToCascCores[traCascCore]); } @@ -621,37 +689,37 @@ struct StrangenessBuilder { int nV0s = 0; // Loops over all V0s in the time frame for (auto& v0 : v0s) { - if(!mEnabledTables[kV0CoresBase] && v0Map[v0.globalIndex()] == -2){ + if (!mEnabledTables[kV0CoresBase] && v0Map[v0.globalIndex()] == -2) { // this v0 hasn't been used by cascades and we're not generating V0s, so skip it v0dataLink(-1, -1); - continue; + continue; } // Get tracks and generate candidate auto const& collision = v0.collision(); auto const& posTrack = v0.template posTrack_as(); auto const& negTrack = v0.template negTrack_as(); - if(!straHelper.buildV0Candidate(collision, posTrack, negTrack, v0.isCollinearV0(), mEnabledTables[kV0Covs])){ + if (!straHelper.buildV0Candidate(collision, posTrack, negTrack, v0.isCollinearV0(), mEnabledTables[kV0Covs])) { v0dataLink(-1, -1); continue; } nV0s++; - if(v0Map[v0.globalIndex()]==-1){ + if (v0Map[v0.globalIndex()] == -1) { v0Map[v0.globalIndex()] = v0sFromCascades.size(); // provide actual valid index in buffer v0sFromCascades.push_back(straHelper.v0); } // fill requested cursors only if type is not 0 if (v0.v0Type() == 1 || (v0.v0Type() == 1 && v0BuilderOpts.generatePhotonCandidates)) { - if(mEnabledTables[kV0Indices]){ + if (mEnabledTables[kV0Indices]) { // for referencing (especially - but not only - when using derived data) v0indices(v0.posTrackId(), v0.negTrackId(), v0.collisionId(), v0.globalIndex()); } - if(mEnabledTables[kV0TrackXs]){ + if (mEnabledTables[kV0TrackXs]) { // further decay chains may need this v0trackXs(straHelper.v0.positiveTrackX, straHelper.v0.negativeTrackX); } - if(mEnabledTables[kV0CoresBase]){ + if (mEnabledTables[kV0CoresBase]) { // standard analysis v0cores(straHelper.v0.position[0], straHelper.v0.position[1], straHelper.v0.position[2], straHelper.v0.positiveMomentum[0], straHelper.v0.positiveMomentum[1], straHelper.v0.positiveMomentum[2], @@ -664,13 +732,12 @@ struct StrangenessBuilder { v0.v0Type()); v0dataLink(v0cores.lastIndex(), -1); } - if(mEnabledTables[kV0TraPosAtDCAs]){ + if (mEnabledTables[kV0TraPosAtDCAs]) { // for tracking studies v0dauPositions(straHelper.v0.positivePosition[0], straHelper.v0.positivePosition[1], straHelper.v0.positivePosition[2], straHelper.v0.negativePosition[0], straHelper.v0.negativePosition[1], straHelper.v0.negativePosition[2]); - } - if(mEnabledTables[kV0TraPosAtIUs]){ + if (mEnabledTables[kV0TraPosAtIUs]) { // for tracking studies std::array positivePositionIU; std::array negativePositionIU; @@ -681,7 +748,7 @@ struct StrangenessBuilder { v0dauPositionsIU(positivePositionIU[0], positivePositionIU[1], positivePositionIU[2], negativePositionIU[0], negativePositionIU[1], negativePositionIU[2]); } - if(mEnabledTables[kV0Covs]){ + if (mEnabledTables[kV0Covs]) { v0covs(straHelper.v0.positionCovariance, straHelper.v0.momentumCovariance); } @@ -689,7 +756,7 @@ struct StrangenessBuilder { // MC handling part if constexpr (requires { TMCParticles::iterator; }) { // only worry about this if someone else worried about this - if((mEnabledTables[kV0MCCores] || mEnabledTables[kMcV0Labels] || mEnabledTables[kV0MCCollRefs])){ + if ((mEnabledTables[kV0MCCores] || mEnabledTables[kMcV0Labels] || mEnabledTables[kV0MCCollRefs])) { thisInfo.label = -1; thisInfo.motherLabel = -1; thisInfo.pdgCode = 0; @@ -756,7 +823,7 @@ struct StrangenessBuilder { } // end association check // Construct label table (note: this will be joinable with V0Datas!) - if(mEnabledTables[kMcV0Labels]){ + if (mEnabledTables[kMcV0Labels]) { v0labels(thisInfo.label, thisInfo.motherLabel); } @@ -770,7 +837,7 @@ struct StrangenessBuilder { // this is the most pedagogical approach, but it is also more limited // and it might use more disk space unnecessarily. if (v0BuilderOpts.mc_populateV0MCCoresSymmetric) { - if(mEnabledTables[kV0MCCores]){ + if (mEnabledTables[kV0MCCores]) { v0mccores( thisInfo.label, thisInfo.pdgCode, thisInfo.pdgCodeMother, thisInfo.pdgCodePositive, thisInfo.pdgCodeNegative, @@ -779,14 +846,14 @@ struct StrangenessBuilder { thisInfo.negP[0], thisInfo.negP[1], thisInfo.negP[2], thisInfo.momentum[0], thisInfo.momentum[1], thisInfo.momentum[2]); } - if(mEnabledTables[kV0MCCollRefs]){ + if (mEnabledTables[kV0MCCollRefs]) { v0mccollref(thisInfo.mcCollision); } // n.b. placing the interlink index here allows for the writing of // code that is agnostic with respect to the joinability of // V0Cores and V0MCCores (always dereference -> safe) - if(mEnabledTables[kV0CoreMCLabels]){ + if (mEnabledTables[kV0CoreMCLabels]) { v0CoreMCLabels(v0.globalIndex()); // interlink index } } @@ -801,7 +868,7 @@ struct StrangenessBuilder { for (uint32_t ii = 0; ii < mcV0infos.size(); ii++) { if (thisInfo.label == mcV0infos[ii].label && mcV0infos[ii].label > -1) { thisV0MCCoreIndex = ii; - break; // this exists already in list + break; // this exists already in list } } if (thisV0MCCoreIndex < 0 && thisInfo.label > -1) { @@ -809,11 +876,11 @@ struct StrangenessBuilder { thisV0MCCoreIndex = mcV0infos.size(); mcV0infos.push_back(thisInfo); } - if(mEnabledTables[kV0CoreMCLabels]){ + if (mEnabledTables[kV0CoreMCLabels]) { v0CoreMCLabels(thisV0MCCoreIndex); // interlink index } } - } // enabled tables check + } // enabled tables check } // constexpr requires check } } @@ -897,7 +964,7 @@ struct StrangenessBuilder { } for (auto info : mcV0infos) { - if(mEnabledTables[kV0MCCores]){ + if (mEnabledTables[kV0MCCores]) { v0mccores( info.label, info.pdgCode, info.pdgCodeMother, info.pdgCodePositive, info.pdgCodeNegative, @@ -906,7 +973,7 @@ struct StrangenessBuilder { info.negP[0], info.negP[1], info.negP[2], info.momentum[0], info.momentum[1], info.momentum[2]); } - if(mEnabledTables[kV0MCCollRefs]){ + if (mEnabledTables[kV0MCCollRefs]) { v0mccollref(info.mcCollision); } } @@ -1008,7 +1075,7 @@ struct StrangenessBuilder { } // end v0 mother loop } // end has_mothers check for V0 } // end conditional of pos/neg originating being the same - } // end association check + } // end association check } //__________________________________________________ @@ -1024,7 +1091,7 @@ struct StrangenessBuilder { mcParticleIsReco.resize(mcParticles.size(), false); } - if(!mEnabledTables[kStoredCascCores]){ + if (!mEnabledTables[kStoredCascCores]) { return; // don't do if no request for cascades in place } int nCascades = 0; @@ -1036,14 +1103,14 @@ struct StrangenessBuilder { auto const& posTrack = v0.template posTrack_as(); auto const& negTrack = v0.template negTrack_as(); auto const& bachTrack = cascade.template bachelor_as(); - if(!straHelper.buildCascadeCandidate(collision, - v0sFromCascades[v0Map[v0.globalIndex()]], - posTrack, - negTrack, - bachTrack, - mEnabledTables[kCascBBs], - cascadeBuilderOpts.useCascadeMomentumAtPrimVtx, - mEnabledTables[kCascCovs])){ + if (!straHelper.buildCascadeCandidate(collision, + v0sFromCascades[v0Map[v0.globalIndex()]], + posTrack, + negTrack, + bachTrack, + mEnabledTables[kCascBBs], + cascadeBuilderOpts.useCascadeMomentumAtPrimVtx, + mEnabledTables[kCascCovs])) { cascdataLink(-1); interlinks.cascadeToCascCores.push_back(-1); continue; // didn't work out, skip @@ -1051,12 +1118,12 @@ struct StrangenessBuilder { nCascades++; // generate analysis tables as required - if(mEnabledTables[kCascIndices]){ + if (mEnabledTables[kCascIndices]) { cascidx(cascade.globalIndex(), straHelper.cascade.positiveTrack, straHelper.cascade.negativeTrack, straHelper.cascade.bachelorTrack, straHelper.cascade.collisionId); } - if (mEnabledTables[kStoredCascCores]){ + if (mEnabledTables[kStoredCascCores]) { cascdata(straHelper.cascade.charge, straHelper.cascade.massXi, straHelper.cascade.massOmega, straHelper.cascade.cascadePosition[0], straHelper.cascade.cascadePosition[1], straHelper.cascade.cascadePosition[2], straHelper.cascade.v0Position[0], straHelper.cascade.v0Position[1], straHelper.cascade.v0Position[2], @@ -1076,7 +1143,7 @@ struct StrangenessBuilder { if (mEnabledTables[kCascTrackXs]) { cascTrackXs(straHelper.cascade.positiveTrackX, straHelper.cascade.negativeTrackX, straHelper.cascade.bachelorTrackX); } - if( mEnabledTables[kCascBBs]){ + if (mEnabledTables[kCascBBs]) { cascbb(straHelper.cascade.bachBaryonCosPA, straHelper.cascade.bachBaryonDCAxyToPV); } if (mEnabledTables[kCascCovs]) { @@ -1087,11 +1154,11 @@ struct StrangenessBuilder { // MC handling part if constexpr (requires { TMCParticles::iterator; }) { // only worry about this if someone else worried about this - if((mEnabledTables[kCascMCCores] || mEnabledTables[kMcCascLabels] || mEnabledTables[kCascMCCollRefs])){ + if ((mEnabledTables[kCascMCCores] || mEnabledTables[kMcCascLabels] || mEnabledTables[kCascMCCollRefs])) { extractMonteCarloProperties(posTrack, negTrack, bachTrack, mcParticles); // Construct label table (note: this will be joinable with CascDatas) - if(mEnabledTables[kMcCascLabels]){ + if (mEnabledTables[kMcCascLabels]) { casclabels( thisCascInfo.label, thisCascInfo.motherLabel); } @@ -1102,18 +1169,18 @@ struct StrangenessBuilder { } if (cascadeBuilderOpts.mc_populateCascMCCoresSymmetric) { - if(mEnabledTables[kCascMCCores]){ - cascmccores( - thisCascInfo.pdgCode, thisCascInfo.pdgCodeMother, thisCascInfo.pdgCodeV0, thisCascInfo.isPhysicalPrimary, - thisCascInfo.pdgCodePositive, thisCascInfo.pdgCodeNegative, thisCascInfo.pdgCodeBachelor, - thisCascInfo.xyz[0], thisCascInfo.xyz[1], thisCascInfo.xyz[2], - thisCascInfo.lxyz[0], thisCascInfo.lxyz[1], thisCascInfo.lxyz[2], - thisCascInfo.posP[0], thisCascInfo.posP[1], thisCascInfo.posP[2], - thisCascInfo.negP[0], thisCascInfo.negP[1], thisCascInfo.negP[2], - thisCascInfo.bachP[0], thisCascInfo.bachP[1], thisCascInfo.bachP[2], - thisCascInfo.momentum[0], thisCascInfo.momentum[1], thisCascInfo.momentum[2]); + if (mEnabledTables[kCascMCCores]) { + cascmccores( + thisCascInfo.pdgCode, thisCascInfo.pdgCodeMother, thisCascInfo.pdgCodeV0, thisCascInfo.isPhysicalPrimary, + thisCascInfo.pdgCodePositive, thisCascInfo.pdgCodeNegative, thisCascInfo.pdgCodeBachelor, + thisCascInfo.xyz[0], thisCascInfo.xyz[1], thisCascInfo.xyz[2], + thisCascInfo.lxyz[0], thisCascInfo.lxyz[1], thisCascInfo.lxyz[2], + thisCascInfo.posP[0], thisCascInfo.posP[1], thisCascInfo.posP[2], + thisCascInfo.negP[0], thisCascInfo.negP[1], thisCascInfo.negP[2], + thisCascInfo.bachP[0], thisCascInfo.bachP[1], thisCascInfo.bachP[2], + thisCascInfo.momentum[0], thisCascInfo.momentum[1], thisCascInfo.momentum[2]); } - if(mEnabledTables[kCascMCCollRefs]){ + if (mEnabledTables[kCascMCCollRefs]) { cascmccollrefs(thisCascInfo.mcCollision); } } @@ -1133,15 +1200,15 @@ struct StrangenessBuilder { thisCascMCCoreIndex = mcCascinfos.size(); mcCascinfos.push_back(thisCascInfo); } - if(mEnabledTables[kCascCoreMCLabels]){ + if (mEnabledTables[kCascCoreMCLabels]) { cascCoreMClabels(thisCascMCCoreIndex); // interlink: reconstructed -> MC index } } - } // enabled tables check + } // enabled tables check // if BB tags requested, generate them now - if(mEnabledTables[kMcCascBBTags]){ + if (mEnabledTables[kMcCascBBTags]) { bool bbTag = false; if (bachTrack.has_mcParticle()) { auto bachelorParticle = bachTrack.template mcParticle_as(); @@ -1158,7 +1225,7 @@ struct StrangenessBuilder { } } } - } // end if-pion + } // end if-pion if (bachelorParticle.pdgCode() == -211) { // pi-, look for proton in positive prong if (posTrack.has_mcParticle()) { auto baryonParticle = posTrack.template mcParticle_as(); @@ -1173,18 +1240,18 @@ struct StrangenessBuilder { } } } // end if-pion - } // end bachelor has mcparticle + } // end bachelor has mcparticle // Construct label table (note: this will be joinable with CascDatas) bbtags(bbTag); } // end BB tag table enabled check - } // constexpr requires mcParticles check + } // constexpr requires mcParticles check } // cascades loop //_________________________________________________________ // MC handling part if constexpr (requires { TMCParticles::iterator; }) { - if((mEnabledTables[kCascMCCores] || mEnabledTables[kMcCascLabels] || mEnabledTables[kCascMCCollRefs])){ + if ((mEnabledTables[kCascMCCores] || mEnabledTables[kMcCascLabels] || mEnabledTables[kCascMCCollRefs])) { // now populate V0MCCores if in asymmetric mode if (cascadeBuilderOpts.mc_populateCascMCCoresAsymmetric) { // first step: add any un-recoed v0mmcores that were requested @@ -1279,25 +1346,25 @@ struct StrangenessBuilder { } for (auto thisInfoToFill : mcCascinfos) { - if(mEnabledTables[kCascMCCores]){ - cascmccores( // a lot of the info below will be compressed in case of not-recoed MC (good!) - thisInfoToFill.pdgCode, thisInfoToFill.pdgCodeMother, thisInfoToFill.pdgCodeV0, thisInfoToFill.isPhysicalPrimary, - thisInfoToFill.pdgCodePositive, thisInfoToFill.pdgCodeNegative, thisInfoToFill.pdgCodeBachelor, - thisInfoToFill.xyz[0], thisInfoToFill.xyz[1], thisInfoToFill.xyz[2], - thisInfoToFill.lxyz[0], thisInfoToFill.lxyz[1], thisInfoToFill.lxyz[2], - thisInfoToFill.posP[0], thisInfoToFill.posP[1], thisInfoToFill.posP[2], - thisInfoToFill.negP[0], thisInfoToFill.negP[1], thisInfoToFill.negP[2], - thisInfoToFill.bachP[0], thisInfoToFill.bachP[1], thisInfoToFill.bachP[2], - thisInfoToFill.momentum[0], thisInfoToFill.momentum[1], thisInfoToFill.momentum[2]); + if (mEnabledTables[kCascMCCores]) { + cascmccores( // a lot of the info below will be compressed in case of not-recoed MC (good!) + thisInfoToFill.pdgCode, thisInfoToFill.pdgCodeMother, thisInfoToFill.pdgCodeV0, thisInfoToFill.isPhysicalPrimary, + thisInfoToFill.pdgCodePositive, thisInfoToFill.pdgCodeNegative, thisInfoToFill.pdgCodeBachelor, + thisInfoToFill.xyz[0], thisInfoToFill.xyz[1], thisInfoToFill.xyz[2], + thisInfoToFill.lxyz[0], thisInfoToFill.lxyz[1], thisInfoToFill.lxyz[2], + thisInfoToFill.posP[0], thisInfoToFill.posP[1], thisInfoToFill.posP[2], + thisInfoToFill.negP[0], thisInfoToFill.negP[1], thisInfoToFill.negP[2], + thisInfoToFill.bachP[0], thisInfoToFill.bachP[1], thisInfoToFill.bachP[2], + thisInfoToFill.momentum[0], thisInfoToFill.momentum[1], thisInfoToFill.momentum[2]); } - if(mEnabledTables[kCascMCCollRefs]){ + if (mEnabledTables[kCascMCCollRefs]) { cascmccollrefs(thisInfoToFill.mcCollision); } } } } // enabled tables check } // constexpr requires mcParticles check - + LOGF(info, "Cascades in DF: %i, cascades built: %i", cascades.size(), nCascades); } @@ -1305,7 +1372,7 @@ struct StrangenessBuilder { template void buildKFCascades(TCascades const& cascades, TMCParticles const& mcParticles) { - if(!mEnabledTables[kStoredKFCascCores]){ + if (!mEnabledTables[kStoredKFCascCores]) { return; // don't do if no request for cascades in place } int nCascades = 0; @@ -1317,17 +1384,17 @@ struct StrangenessBuilder { auto const& posTrack = v0.template posTrack_as(); auto const& negTrack = v0.template negTrack_as(); auto const& bachTrack = cascade.template bachelor_as(); - if(!straHelper.buildCascadeCandidateWithKF(collision, - posTrack, - negTrack, - bachTrack, - mEnabledTables[kCascBBs], - cascadeBuilderOpts.kfConstructMethod, - cascadeBuilderOpts.kfTuneForOmega, - cascadeBuilderOpts.kfUseV0MassConstraint, - cascadeBuilderOpts.kfUseCascadeMassConstraint, - cascadeBuilderOpts.kfDoDCAFitterPreMinimV0, - cascadeBuilderOpts.kfDoDCAFitterPreMinimCasc)){ + if (!straHelper.buildCascadeCandidateWithKF(collision, + posTrack, + negTrack, + bachTrack, + mEnabledTables[kCascBBs], + cascadeBuilderOpts.kfConstructMethod, + cascadeBuilderOpts.kfTuneForOmega, + cascadeBuilderOpts.kfUseV0MassConstraint, + cascadeBuilderOpts.kfUseCascadeMassConstraint, + cascadeBuilderOpts.kfDoDCAFitterPreMinimV0, + cascadeBuilderOpts.kfDoDCAFitterPreMinimCasc)) { kfcascdataLink(-1); interlinks.cascadeToKFCascCores.push_back(-1); continue; // didn't work out, skip @@ -1335,12 +1402,12 @@ struct StrangenessBuilder { nCascades++; // generate analysis tables as required - if(mEnabledTables[kKFCascIndices]){ + if (mEnabledTables[kKFCascIndices]) { kfcascidx(cascade.globalIndex(), straHelper.cascade.positiveTrack, straHelper.cascade.negativeTrack, straHelper.cascade.bachelorTrack, straHelper.cascade.collisionId); } - if (mEnabledTables[kStoredKFCascCores]){ + if (mEnabledTables[kStoredKFCascCores]) { kfcascdata(straHelper.cascade.charge, straHelper.cascade.massXi, straHelper.cascade.massOmega, straHelper.cascade.cascadePosition[0], straHelper.cascade.cascadePosition[1], straHelper.cascade.cascadePosition[2], straHelper.cascade.v0Position[0], straHelper.cascade.v0Position[1], straHelper.cascade.v0Position[2], @@ -1368,15 +1435,15 @@ struct StrangenessBuilder { // MC handling part (labels only) if constexpr (requires { TMCParticles::iterator; }) { // only worry about this if someone else worried about this - if((mEnabledTables[kMcKFCascLabels])){ + if ((mEnabledTables[kMcKFCascLabels])) { extractMonteCarloProperties(posTrack, negTrack, bachTrack, mcParticles); // Construct label table (note: this will be joinable with KFCascDatas) kfcasclabels(thisCascInfo.label, thisCascInfo.motherLabel); - } // enabled tables check + } // enabled tables check } // constexpr requires mcParticles check } // end loop over cascades - + LOGF(info, "KF Cascades in DF: %i, KF cascades built: %i", cascades.size(), nCascades); } @@ -1384,7 +1451,7 @@ struct StrangenessBuilder { template void buildTrackedCascades(TStrangeTracks const& cascadeTracks, TMCParticles const& mcParticles) { - if(!mEnabledTables[kStoredTraCascCores]){ + if (!mEnabledTables[kStoredTraCascCores]) { return; // don't do if no request for cascades in place } int nCascades = 0; @@ -1401,13 +1468,13 @@ struct StrangenessBuilder { auto const& posTrack = v0.template posTrack_as(); auto const& negTrack = v0.template negTrack_as(); auto const& bachTrack = cascade.template bachelor_as(); - if(!straHelper.buildCascadeCandidate(collision, - posTrack, - negTrack, - bachTrack, - mEnabledTables[kCascBBs], - cascadeBuilderOpts.useCascadeMomentumAtPrimVtx, - mEnabledTables[kCascCovs])){ + if (!straHelper.buildCascadeCandidate(collision, + posTrack, + negTrack, + bachTrack, + mEnabledTables[kCascBBs], + cascadeBuilderOpts.useCascadeMomentumAtPrimVtx, + mEnabledTables[kCascCovs])) { tracascdataLink(-1); interlinks.cascadeToTraCascCores.push_back(-1); continue; // didn't work out, skip @@ -1418,8 +1485,8 @@ struct StrangenessBuilder { gpu::gpustd::array dcaInfo; strangeTrackParCov.setPID(o2::track::PID::XiMinus); // FIXME: not OK for omegas o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, strangeTrackParCov, 2.f, straHelper.fitter.getMatCorrType(), &dcaInfo); - straHelper.cascade.cascadeDCAxy = dcaInfo[0]; - straHelper.cascade.cascadeDCAz = dcaInfo[1]; + straHelper.cascade.cascadeDCAxy = dcaInfo[0]; + straHelper.cascade.cascadeDCAz = dcaInfo[1]; // get momentum from strange track (should not be very different) strangeTrackParCov.getPxPyPzGlo(straHelper.cascade.cascadeMomentum); @@ -1428,23 +1495,23 @@ struct StrangenessBuilder { nCascades++; // generate analysis tables as required - if(mEnabledTables[kTraCascIndices]){ + if (mEnabledTables[kTraCascIndices]) { tracascidx(cascade.globalIndex(), straHelper.cascade.positiveTrack, straHelper.cascade.negativeTrack, straHelper.cascade.bachelorTrack, straHelper.cascade.collisionId); } - if (mEnabledTables[kStoredTraCascCores]){ + if (mEnabledTables[kStoredTraCascCores]) { tracascdata(straHelper.cascade.charge, cascadeTrack.xiMass(), cascadeTrack.omegaMass(), - cascadeTrack.decayX(), cascadeTrack.decayY(), cascadeTrack.decayZ(), - straHelper.cascade.v0Position[0], straHelper.cascade.v0Position[1], straHelper.cascade.v0Position[2], - straHelper.cascade.positiveMomentum[0], straHelper.cascade.positiveMomentum[1], straHelper.cascade.positiveMomentum[2], - straHelper.cascade.negativeMomentum[0], straHelper.cascade.negativeMomentum[1], straHelper.cascade.negativeMomentum[2], - straHelper.cascade.bachelorMomentum[0], straHelper.cascade.bachelorMomentum[1], straHelper.cascade.bachelorMomentum[2], - straHelper.cascade.cascadeMomentum[0], straHelper.cascade.cascadeMomentum[1], straHelper.cascade.cascadeMomentum[2], - straHelper.cascade.v0DaughterDCA, straHelper.cascade.cascadeDaughterDCA, - straHelper.cascade.positiveDCAxy, straHelper.cascade.negativeDCAxy, - straHelper.cascade.bachelorDCAxy, straHelper.cascade.cascadeDCAxy, straHelper.cascade.cascadeDCAz, - cascadeTrack.matchingChi2(), cascadeTrack.topologyChi2(), cascadeTrack.itsClsSize()); + cascadeTrack.decayX(), cascadeTrack.decayY(), cascadeTrack.decayZ(), + straHelper.cascade.v0Position[0], straHelper.cascade.v0Position[1], straHelper.cascade.v0Position[2], + straHelper.cascade.positiveMomentum[0], straHelper.cascade.positiveMomentum[1], straHelper.cascade.positiveMomentum[2], + straHelper.cascade.negativeMomentum[0], straHelper.cascade.negativeMomentum[1], straHelper.cascade.negativeMomentum[2], + straHelper.cascade.bachelorMomentum[0], straHelper.cascade.bachelorMomentum[1], straHelper.cascade.bachelorMomentum[2], + straHelper.cascade.cascadeMomentum[0], straHelper.cascade.cascadeMomentum[1], straHelper.cascade.cascadeMomentum[2], + straHelper.cascade.v0DaughterDCA, straHelper.cascade.cascadeDaughterDCA, + straHelper.cascade.positiveDCAxy, straHelper.cascade.negativeDCAxy, + straHelper.cascade.bachelorDCAxy, straHelper.cascade.cascadeDCAxy, straHelper.cascade.cascadeDCAz, + cascadeTrack.matchingChi2(), cascadeTrack.topologyChi2(), cascadeTrack.itsClsSize()); // interlink always produced if base core table generated tracascdataLink(tracascdata.lastIndex()); interlinks.traCascCoreToCascades.push_back(cascade.globalIndex()); @@ -1464,12 +1531,12 @@ struct StrangenessBuilder { // MC handling part (labels only) if constexpr (requires { TMCParticles::iterator; }) { // only worry about this if someone else worried about this - if((mEnabledTables[kMcTraCascLabels])){ + if ((mEnabledTables[kMcTraCascLabels])) { extractMonteCarloProperties(posTrack, negTrack, bachTrack, mcParticles); // Construct label table (note: this will be joinable with KFCascDatas) tracasclabels(thisCascInfo.label, thisCascInfo.motherLabel); - } // enabled tables check + } // enabled tables check } // constexpr requires mcParticles check } // end loop over cascades LOGF(info, "Tracked cascades in DF: %i, tracked cascades built: %i", cascadeTracks.size(), nCascades); @@ -1510,7 +1577,8 @@ struct StrangenessBuilder { template void dataProcess(TCollisions const& collisions, TV0s const& v0s, TCascades const& cascades, TTrackedCascades const& trackedCascades, TTracks const&, TBCs const& bcs, TMCParticles const& mcParticles) { - if(!initCCDB(bcs, collisions)) return; + if (!initCCDB(bcs, collisions)) + return; // reset vectors for cascade interlinks resetInterlinks(); @@ -1518,9 +1586,9 @@ struct StrangenessBuilder { // mark V0s that will be buffered for the cascade building markV0sUsedInCascades(v0s, cascades); - // build V0s + // build V0s buildV0s(v0s, mcParticles); - + // build cascades buildCascades(cascades, mcParticles); buildKFCascades(cascades, mcParticles); @@ -1535,12 +1603,12 @@ struct StrangenessBuilder { void processRealData(aod::Collisions const& collisions, aod::V0s const& v0s, aod::Cascades const& cascades, aod::TrackedCascades const& trackedCascades, FullTracksExtIU const& tracks, aod::BCsWithTimestamps const& bcs) { - dataProcess(collisions, v0s, cascades, trackedCascades, tracks, bcs, (TObject*) nullptr); + dataProcess(collisions, v0s, cascades, trackedCascades, tracks, bcs, (TObject*)nullptr); } void processRealDataRun2(aod::Collisions const& collisions, aod::V0s const& v0s, aod::Cascades const& cascades, FullTracksExt const& tracks, aod::BCsWithTimestamps const& bcs) { - dataProcess(collisions, v0s, cascades, (TObject*) nullptr, tracks, bcs, (TObject*) nullptr); + dataProcess(collisions, v0s, cascades, (TObject*)nullptr, tracks, bcs, (TObject*)nullptr); } void processMonteCarlo(aod::Collisions const& collisions, aod::V0s const& v0s, aod::Cascades const& cascades, aod::TrackedCascades const& trackedCascades, FullTracksExtLabeledIU const& tracks, aod::BCsWithTimestamps const& bcs, aod::McParticles const& mcParticles) @@ -1550,12 +1618,12 @@ struct StrangenessBuilder { void processMonteCarloRun2(aod::Collisions const& collisions, aod::V0s const& v0s, aod::Cascades const& cascades, FullTracksExtLabeled const& tracks, aod::BCsWithTimestamps const& bcs, aod::McParticles const& mcParticles) { - dataProcess(collisions, v0s, cascades, (TObject*) nullptr, tracks, bcs, mcParticles); + dataProcess(collisions, v0s, cascades, (TObject*)nullptr, tracks, bcs, mcParticles); } void processSimulationFindable(aod::Collisions const& collisions, aod::FindableV0s const& v0s, aod::Cascades const& cascades, FullTracksExtIU const& tracks, aod::BCsWithTimestamps const& bcs) { - dataProcess(collisions, v0s, cascades, (TObject*) nullptr, tracks, bcs, (TObject*) nullptr); + dataProcess(collisions, v0s, cascades, (TObject*)nullptr, tracks, bcs, (TObject*)nullptr); } PROCESS_SWITCH(StrangenessBuilder, processRealData, "process real data", true); diff --git a/PWGLF/Utils/strangenessBuilderHelper.h b/PWGLF/Utils/strangenessBuilderHelper.h index 37db56280ab..eda4ec6a556 100644 --- a/PWGLF/Utils/strangenessBuilderHelper.h +++ b/PWGLF/Utils/strangenessBuilderHelper.h @@ -39,7 +39,7 @@ namespace pwglf { //__________________________________________ // V0 information storage -struct v0candidate{ +struct v0candidate { // indexing int collisionId = -1; int negativeTrack = -1; @@ -50,17 +50,17 @@ struct v0candidate{ std::array negativeMomentum = {0.0f, 0.0f, 0.0f}; std::array positivePosition = {0.0f, 0.0f, 0.0f}; std::array negativePosition = {0.0f, 0.0f, 0.0f}; - float positiveTrackX = 0.0f; + float positiveTrackX = 0.0f; float negativeTrackX = 0.0f; float positiveDCAxy = 0.0f; float negativeDCAxy = 0.0f; - + // V0 properties std::array position = {0.0f, 0.0f, 0.0f}; float daughterDCA = 1000.0f; float pointingAngle = 0.0f; float dcaXY = 0.0f; - + // calculated masses for convenience float massGamma; float massK0Short; @@ -92,7 +92,7 @@ struct cascadeCandidate { float positiveDCAxy = 0.0f; float negativeDCAxy = 0.0f; float bachelorDCAxy = 0.0f; - float positiveTrackX = 0.0f; + float positiveTrackX = 0.0f; float negativeTrackX = 0.0f; float bachelorTrackX = 0.0f; @@ -100,7 +100,7 @@ struct cascadeCandidate { int charge = -1; // default: []Minus float cascadeDaughterDCA = 1000.0f; float v0DaughterDCA = 1000.0f; - + float pointingAngle = 0.0f; float cascadeDCAxy = 0.0f; float cascadeDCAz = 0.0f; @@ -111,8 +111,8 @@ struct cascadeCandidate { float bachBaryonCosPA = 0.0f; float bachBaryonDCAxyToPV = 1e+3; - float massXi = 0.0f; - float massOmega = 0.0f; + float massXi = 0.0f; + float massOmega = 0.0f; // KF-specific variables float kfV0Chi2 = 0.0f; @@ -130,716 +130,718 @@ struct cascadeCandidate { // builder helper class class strangenessBuilderHelper { - public: - strangenessBuilderHelper(){ - // standards hardcoded in builder ... - // ...but can be changed easily since fitter is public - fitter.setPropagateToPCA(true); - fitter.setMaxR(200.); - fitter.setMinParamChange(1e-3); - fitter.setMinRelChi2Change(0.9); - fitter.setMaxDZIni(1e9); - fitter.setMaxDXYIni(4.0f); - fitter.setMaxChi2(1e9); - fitter.setUseAbsDCA(true); - fitter.setWeightedFinalPCA(false); - - // LUT has to be loaded later - lut = nullptr; - fitter.setMatCorrType(o2::base::Propagator::MatCorrType::USEMatCorrLUT); - - // mag field has to be set later - fitter.setBz(-999.9f); // will NOT make sense if not changed - }; - - template - bool buildV0Candidate(o2::aod::Collision const& collision, - TTrack const& positiveTrack, - TTrack const& negativeTrack, - bool useCollinearFit = false, - bool calculateCovariance = false){ - // Calculate DCA with respect to the collision associated to the V0, not individual tracks - gpu::gpustd::array dcaInfo; - - auto posTrackPar = getTrackPar(positiveTrack); - o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, posTrackPar, 2.f, fitter.getMatCorrType(), &dcaInfo); - v0.positiveDCAxy = dcaInfo[0]; - - auto negTrackPar = getTrackPar(negativeTrack); - o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, negTrackPar, 2.f, fitter.getMatCorrType(), &dcaInfo); - v0.negativeDCAxy = dcaInfo[0]; - - o2::track::TrackParCov positiveTrackParam = getTrackParCov(positiveTrack); - o2::track::TrackParCov negativeTrackParam = getTrackParCov(negativeTrack); - - // Perform DCA fit - int nCand = 0; - fitter.setCollinear(useCollinearFit); - try { - nCand = fitter.process(positiveTrackParam, negativeTrackParam); - } catch (...) { - return false; - } - if (nCand == 0) { - return false; - } + public: + strangenessBuilderHelper() + { + // standards hardcoded in builder ... + // ...but can be changed easily since fitter is public + fitter.setPropagateToPCA(true); + fitter.setMaxR(200.); + fitter.setMinParamChange(1e-3); + fitter.setMinRelChi2Change(0.9); + fitter.setMaxDZIni(1e9); + fitter.setMaxDXYIni(4.0f); + fitter.setMaxChi2(1e9); + fitter.setUseAbsDCA(true); + fitter.setWeightedFinalPCA(false); + + // LUT has to be loaded later + lut = nullptr; + fitter.setMatCorrType(o2::base::Propagator::MatCorrType::USEMatCorrLUT); + + // mag field has to be set later + fitter.setBz(-999.9f); // will NOT make sense if not changed + }; + + template + bool buildV0Candidate(o2::aod::Collision const& collision, + TTrack const& positiveTrack, + TTrack const& negativeTrack, + bool useCollinearFit = false, + bool calculateCovariance = false) + { + // Calculate DCA with respect to the collision associated to the V0, not individual tracks + gpu::gpustd::array dcaInfo; + + auto posTrackPar = getTrackPar(positiveTrack); + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, posTrackPar, 2.f, fitter.getMatCorrType(), &dcaInfo); + v0.positiveDCAxy = dcaInfo[0]; + + auto negTrackPar = getTrackPar(negativeTrack); + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, negTrackPar, 2.f, fitter.getMatCorrType(), &dcaInfo); + v0.negativeDCAxy = dcaInfo[0]; + + o2::track::TrackParCov positiveTrackParam = getTrackParCov(positiveTrack); + o2::track::TrackParCov negativeTrackParam = getTrackParCov(negativeTrack); + + // Perform DCA fit + int nCand = 0; + fitter.setCollinear(useCollinearFit); + try { + nCand = fitter.process(positiveTrackParam, negativeTrackParam); + } catch (...) { + return false; + } + if (nCand == 0) { + return false; + } - v0.positiveTrackX = fitter.getTrack(0).getX(); - v0.negativeTrackX = fitter.getTrack(1).getX(); - positiveTrackParam = fitter.getTrack(0); - negativeTrackParam = fitter.getTrack(1); - positiveTrackParam.getPxPyPzGlo(v0.positiveMomentum); - negativeTrackParam.getPxPyPzGlo(v0.negativeMomentum); - positiveTrackParam.getXYZGlo(v0.positivePosition); - negativeTrackParam.getXYZGlo(v0.negativePosition); - - // get decay vertex coordinates - const auto& vtx = fitter.getPCACandidate(); - for (int i = 0; i < 3; i++) { - v0.position[i] = vtx[i]; - } + v0.positiveTrackX = fitter.getTrack(0).getX(); + v0.negativeTrackX = fitter.getTrack(1).getX(); + positiveTrackParam = fitter.getTrack(0); + negativeTrackParam = fitter.getTrack(1); + positiveTrackParam.getPxPyPzGlo(v0.positiveMomentum); + negativeTrackParam.getPxPyPzGlo(v0.negativeMomentum); + positiveTrackParam.getXYZGlo(v0.positivePosition); + negativeTrackParam.getXYZGlo(v0.negativePosition); + + // get decay vertex coordinates + const auto& vtx = fitter.getPCACandidate(); + for (int i = 0; i < 3; i++) { + v0.position[i] = vtx[i]; + } - v0.daughterDCA = TMath::Sqrt(fitter.getChi2AtPCACandidate()); - v0.pointingAngle = TMath::ACos(RecoDecay::cpa( - std::array{collision.posX(), collision.posY(), collision.posZ()}, - std::array{v0.position[0], v0.position[1], v0.position[2]}, - std::array{v0.positiveMomentum[0] + v0.negativeMomentum[0], v0.positiveMomentum[1] + v0.negativeMomentum[1], v0.positiveMomentum[2] + v0.negativeMomentum[2]} - )); - - v0.dcaXY = CalculateDCAStraightToPV( - v0.position[0], v0.position[1], v0.position[2], - v0.positiveMomentum[0] + v0.negativeMomentum[0], - v0.positiveMomentum[1] + v0.negativeMomentum[1], - v0.positiveMomentum[2] + v0.negativeMomentum[2], - collision.posX(), collision.posY(), collision.posZ()); - - // Calculate masses - v0.massGamma = RecoDecay::m(std::array{ - std::array{v0.positiveMomentum[0], v0.positiveMomentum[1], v0.positiveMomentum[2]}, - std::array{v0.negativeMomentum[0], v0.negativeMomentum[1], v0.negativeMomentum[2]}}, - std::array{o2::constants::physics::MassElectron, o2::constants::physics::MassElectron}); - v0.massK0Short = RecoDecay::m(std::array{ - std::array{v0.positiveMomentum[0], v0.positiveMomentum[1], v0.positiveMomentum[2]}, - std::array{v0.negativeMomentum[0], v0.negativeMomentum[1], v0.negativeMomentum[2]}}, - std::array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassPionCharged}); - v0.massLambda = RecoDecay::m(std::array{ - std::array{v0.positiveMomentum[0], v0.positiveMomentum[1], v0.positiveMomentum[2]}, - std::array{v0.negativeMomentum[0], v0.negativeMomentum[1], v0.negativeMomentum[2]}}, - std::array{o2::constants::physics::MassProton, o2::constants::physics::MassPionCharged}); - v0.massAntiLambda = RecoDecay::m(std::array{ - std::array{v0.positiveMomentum[0], v0.positiveMomentum[1], v0.positiveMomentum[2]}, - std::array{v0.negativeMomentum[0], v0.negativeMomentum[1], v0.negativeMomentum[2]}}, - std::array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassProton}); - - // calculate covariance if requested - if (calculateCovariance) { - // Calculate position covariance matrix - auto covVtxV = fitter.calcPCACovMatrix(0); - // std::array positionCovariance; - v0.positionCovariance[0] = covVtxV(0, 0); - v0.positionCovariance[1] = covVtxV(1, 0); - v0.positionCovariance[2] = covVtxV(1, 1); - v0.positionCovariance[3] = covVtxV(2, 0); - v0.positionCovariance[4] = covVtxV(2, 1); - v0.positionCovariance[5] = covVtxV(2, 2); - std::array covTpositive = {0.}; - std::array covTnegative = {0.}; - positiveTrackParam.getCovXYZPxPyPzGlo(covTpositive); - negativeTrackParam.getCovXYZPxPyPzGlo(covTnegative); - constexpr int MomInd[6] = {9, 13, 14, 18, 19, 20}; // cov matrix elements for momentum component - for (int i = 0; i < 6; i++) { - v0.momentumCovariance[i] = covTpositive[MomInd[i]] + covTnegative[MomInd[i]]; - } + v0.daughterDCA = TMath::Sqrt(fitter.getChi2AtPCACandidate()); + v0.pointingAngle = TMath::ACos(RecoDecay::cpa( + std::array{collision.posX(), collision.posY(), collision.posZ()}, + std::array{v0.position[0], v0.position[1], v0.position[2]}, + std::array{v0.positiveMomentum[0] + v0.negativeMomentum[0], v0.positiveMomentum[1] + v0.negativeMomentum[1], v0.positiveMomentum[2] + v0.negativeMomentum[2]})); + + v0.dcaXY = CalculateDCAStraightToPV( + v0.position[0], v0.position[1], v0.position[2], + v0.positiveMomentum[0] + v0.negativeMomentum[0], + v0.positiveMomentum[1] + v0.negativeMomentum[1], + v0.positiveMomentum[2] + v0.negativeMomentum[2], + collision.posX(), collision.posY(), collision.posZ()); + + // Calculate masses + v0.massGamma = RecoDecay::m(std::array{ + std::array{v0.positiveMomentum[0], v0.positiveMomentum[1], v0.positiveMomentum[2]}, + std::array{v0.negativeMomentum[0], v0.negativeMomentum[1], v0.negativeMomentum[2]}}, + std::array{o2::constants::physics::MassElectron, o2::constants::physics::MassElectron}); + v0.massK0Short = RecoDecay::m(std::array{ + std::array{v0.positiveMomentum[0], v0.positiveMomentum[1], v0.positiveMomentum[2]}, + std::array{v0.negativeMomentum[0], v0.negativeMomentum[1], v0.negativeMomentum[2]}}, + std::array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassPionCharged}); + v0.massLambda = RecoDecay::m(std::array{ + std::array{v0.positiveMomentum[0], v0.positiveMomentum[1], v0.positiveMomentum[2]}, + std::array{v0.negativeMomentum[0], v0.negativeMomentum[1], v0.negativeMomentum[2]}}, + std::array{o2::constants::physics::MassProton, o2::constants::physics::MassPionCharged}); + v0.massAntiLambda = RecoDecay::m(std::array{ + std::array{v0.positiveMomentum[0], v0.positiveMomentum[1], v0.positiveMomentum[2]}, + std::array{v0.negativeMomentum[0], v0.negativeMomentum[1], v0.negativeMomentum[2]}}, + std::array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassProton}); + + // calculate covariance if requested + if (calculateCovariance) { + // Calculate position covariance matrix + auto covVtxV = fitter.calcPCACovMatrix(0); + // std::array positionCovariance; + v0.positionCovariance[0] = covVtxV(0, 0); + v0.positionCovariance[1] = covVtxV(1, 0); + v0.positionCovariance[2] = covVtxV(1, 1); + v0.positionCovariance[3] = covVtxV(2, 0); + v0.positionCovariance[4] = covVtxV(2, 1); + v0.positionCovariance[5] = covVtxV(2, 2); + std::array covTpositive = {0.}; + std::array covTnegative = {0.}; + positiveTrackParam.getCovXYZPxPyPzGlo(covTpositive); + negativeTrackParam.getCovXYZPxPyPzGlo(covTnegative); + constexpr int MomInd[6] = {9, 13, 14, 18, 19, 20}; // cov matrix elements for momentum component + for (int i = 0; i < 6; i++) { + v0.momentumCovariance[i] = covTpositive[MomInd[i]] + covTnegative[MomInd[i]]; } - - // information validated, V0 built successfully. Signal OK - return true; } - // cascade builder creating a cascade from plain tracks - template - bool buildCascadeCandidate(o2::aod::Collision const& collision, - TTrack const& positiveTrack, - TTrack const& negativeTrack, - TTrack const& bachelorTrack, - bool calculateBachelorBaryonVariables = false, - bool useCascadeMomentumAtPV = false, - bool processCovariances = false) - { - if(!buildV0Candidate(collision, positiveTrack, negativeTrack, false, processCovariances)){ - return false; - } - if(!buildCascadeCandidate(collision, v0, positiveTrack, negativeTrack, bachelorTrack, calculateBachelorBaryonVariables, useCascadeMomentumAtPV, processCovariances)){ - return false; + // information validated, V0 built successfully. Signal OK + return true; + } + + // cascade builder creating a cascade from plain tracks + template + bool buildCascadeCandidate(o2::aod::Collision const& collision, + TTrack const& positiveTrack, + TTrack const& negativeTrack, + TTrack const& bachelorTrack, + bool calculateBachelorBaryonVariables = false, + bool useCascadeMomentumAtPV = false, + bool processCovariances = false) + { + if (!buildV0Candidate(collision, positiveTrack, negativeTrack, false, processCovariances)) { + return false; + } + if (!buildCascadeCandidate(collision, v0, positiveTrack, negativeTrack, bachelorTrack, calculateBachelorBaryonVariables, useCascadeMomentumAtPV, processCovariances)) { + return false; + } + return true; + } + + // cascade builder using pre-fabricated information, thus not calling + // the DCAfitter again for the V0 contained in the cascade + // if generating from scratch, prefer the other variant + template + bool buildCascadeCandidate(o2::aod::Collision const& collision, + v0candidate const& v0input, + TTrack const& positiveTrack, + TTrack const& negativeTrack, + TTrack const& bachelorTrack, + bool calculateBachelorBaryonVariables = false, + bool useCascadeMomentumAtPV = false, + bool processCovariances = false) + { + if (calculateBachelorBaryonVariables) { + // Calculates properties of the V0 comprised of bachelor and baryon in the cascade + // baryon: distinguished via bachelor charge + if (bachelorTrack.sign() < 0) { + processBachBaryonVariables(collision, bachelorTrack, positiveTrack); + } else { + processBachBaryonVariables(collision, bachelorTrack, negativeTrack); } - return true; } - // cascade builder using pre-fabricated information, thus not calling - // the DCAfitter again for the V0 contained in the cascade - // if generating from scratch, prefer the other variant - template - bool buildCascadeCandidate(o2::aod::Collision const& collision, - v0candidate const& v0input, - TTrack const& positiveTrack, - TTrack const& negativeTrack, - TTrack const& bachelorTrack, - bool calculateBachelorBaryonVariables = false, - bool useCascadeMomentumAtPV = false, - bool processCovariances = false) - { - if (calculateBachelorBaryonVariables) { - // Calculates properties of the V0 comprised of bachelor and baryon in the cascade - // baryon: distinguished via bachelor charge - if (bachelorTrack.sign() < 0) { - processBachBaryonVariables(collision, bachelorTrack, positiveTrack); - } else { - processBachBaryonVariables(collision, bachelorTrack, negativeTrack); - } - } + // Overall cascade charge + cascade.charge = bachelorTrack.signed1Pt() > 0 ? +1 : -1; - // Overall cascade charge - cascade.charge = bachelorTrack.signed1Pt() > 0 ? +1 : -1; + // bachelor DCA track to PV + // Calculate DCA with respect to the collision associated to the V0, not individual tracks + gpu::gpustd::array dcaInfo; - // bachelor DCA track to PV - // Calculate DCA with respect to the collision associated to the V0, not individual tracks - gpu::gpustd::array dcaInfo; + auto bachTrackPar = getTrackPar(bachelorTrack); + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, bachTrackPar, 2.f, fitter.getMatCorrType(), &dcaInfo); + cascade.bachelorDCAxy = dcaInfo[0]; - auto bachTrackPar = getTrackPar(bachelorTrack); - o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, bachTrackPar, 2.f, fitter.getMatCorrType(), &dcaInfo); - cascade.bachelorDCAxy = dcaInfo[0]; + // Do actual minimization + auto lBachelorTrack = getTrackParCov(bachelorTrack); - // Do actual minimization - auto lBachelorTrack = getTrackParCov(bachelorTrack); + // Set up covariance matrices (should in fact be optional) + std::array covV = {0.}; + constexpr int MomInd[6] = {9, 13, 14, 18, 19, 20}; // cov matrix elements for momentum component + for (int i = 0; i < 6; i++) { + covV[MomInd[i]] = v0input.momentumCovariance[i]; + covV[i] = v0input.positionCovariance[i]; + } + auto lV0Track = o2::track::TrackParCov( + {v0input.position[0], v0input.position[1], v0input.position[2]}, + {v0input.positiveMomentum[0] + v0input.negativeMomentum[0], v0input.positiveMomentum[1] + v0input.negativeMomentum[1], v0input.positiveMomentum[2] + v0input.negativeMomentum[2]}, + covV, 0, true); + lV0Track.setAbsCharge(0); + lV0Track.setPID(o2::track::PID::Lambda); + + //---/---/---/ + // Move close to minima + int nCand = 0; + try { + nCand = fitter.process(lV0Track, lBachelorTrack); + } catch (...) { + return false; + } + if (nCand == 0) + return false; - // Set up covariance matrices (should in fact be optional) - std::array covV = {0.}; - constexpr int MomInd[6] = {9, 13, 14, 18, 19, 20}; // cov matrix elements for momentum component - for (int i = 0; i < 6; i++) { - covV[MomInd[i]] = v0input.momentumCovariance[i]; - covV[i] = v0input.positionCovariance[i]; - } - auto lV0Track = o2::track::TrackParCov( - {v0input.position[0], v0input.position[1], v0input.position[2]}, - {v0input.positiveMomentum[0] + v0input.negativeMomentum[0], v0input.positiveMomentum[1] + v0input.negativeMomentum[1], v0input.positiveMomentum[2] + v0input.negativeMomentum[2]}, - covV, 0, true); - lV0Track.setAbsCharge(0); - lV0Track.setPID(o2::track::PID::Lambda); - - //---/---/---/ - // Move close to minima - int nCand = 0; - try { - nCand = fitter.process(lV0Track, lBachelorTrack); - } catch (...) { - return false; - } - if (nCand == 0) - return false; + lV0Track = fitter.getTrack(0); + lBachelorTrack = fitter.getTrack(1); - lV0Track = fitter.getTrack(0); - lBachelorTrack = fitter.getTrack(1); + // DCA between cascade daughters + cascade.cascadeDaughterDCA = TMath::Sqrt(fitter.getChi2AtPCACandidate()); - // DCA between cascade daughters - cascade.cascadeDaughterDCA = TMath::Sqrt(fitter.getChi2AtPCACandidate()); + lBachelorTrack.getPxPyPzGlo(cascade.bachelorMomentum); + // get decay vertex coordinates + const auto& vtx = fitter.getPCACandidate(); + for (int i = 0; i < 3; i++) { + cascade.cascadePosition[i] = vtx[i]; + } - lBachelorTrack.getPxPyPzGlo(cascade.bachelorMomentum); - // get decay vertex coordinates - const auto& vtx = fitter.getPCACandidate(); - for (int i = 0; i < 3; i++) { - cascade.cascadePosition[i] = vtx[i]; - } + cascade.pointingAngle = TMath::ACos(RecoDecay::cpa( + std::array{collision.posX(), collision.posY(), collision.posZ()}, + std::array{cascade.cascadePosition[0], cascade.cascadePosition[1], cascade.cascadePosition[2]}, + std::array{v0input.positiveMomentum[0] + v0input.negativeMomentum[0] + cascade.bachelorMomentum[0], + v0input.positiveMomentum[0] + v0input.negativeMomentum[1] + cascade.bachelorMomentum[1], + v0input.positiveMomentum[2] + v0input.negativeMomentum[2] + cascade.bachelorMomentum[2]})); + + // Calculate DCAxy of the cascade (with bending) + auto lCascadeTrack = fitter.createParentTrackParCov(); + lCascadeTrack.setAbsCharge(cascade.charge); // to be sure + lCascadeTrack.setPID(o2::track::PID::XiMinus); // FIXME: not OK for omegas + dcaInfo[0] = 999; + dcaInfo[1] = 999; + + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, lCascadeTrack, 2.f, fitter.getMatCorrType(), &dcaInfo); + cascade.cascadeDCAxy = dcaInfo[0]; + cascade.cascadeDCAz = dcaInfo[1]; + + // Calculate masses a priori + cascade.massXi = RecoDecay::m(std::array{std::array{cascade.bachelorMomentum[0], cascade.bachelorMomentum[1], cascade.bachelorMomentum[2]}, + std::array{v0input.positiveMomentum[0] + v0input.negativeMomentum[0], v0input.positiveMomentum[1] + v0input.negativeMomentum[1], v0input.positiveMomentum[2] + v0input.negativeMomentum[2]}}, + std::array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassLambda}); + cascade.massOmega = RecoDecay::m(std::array{std::array{cascade.bachelorMomentum[0], cascade.bachelorMomentum[1], cascade.bachelorMomentum[2]}, + std::array{v0input.positiveMomentum[0] + v0input.negativeMomentum[0], v0input.positiveMomentum[1] + v0input.negativeMomentum[1], v0input.positiveMomentum[2] + v0input.negativeMomentum[2]}}, + std::array{o2::constants::physics::MassKaonCharged, o2::constants::physics::MassLambda}); + + // Populate information + // cascadecandidate.v0Id = v0index.globalIndex(); + cascade.positiveTrack = positiveTrack.globalIndex(); + cascade.negativeTrack = negativeTrack.globalIndex(); + cascade.bachelorTrack = bachelorTrack.globalIndex(); + cascade.positiveTrackX = v0input.positiveTrackX; + cascade.negativeTrackX = v0input.negativeTrackX; + cascade.bachelorTrackX = lBachelorTrack.getX(); // from this minimization + cascade.v0Position[0] = v0input.position[0]; + cascade.v0Position[1] = v0input.position[1]; + cascade.v0Position[2] = v0input.position[2]; + cascade.v0Momentum[0] = v0input.positiveMomentum[0] + v0input.negativeMomentum[0]; + cascade.v0Momentum[1] = v0input.positiveMomentum[1] + v0input.negativeMomentum[1]; + cascade.v0Momentum[2] = v0input.positiveMomentum[2] + v0input.negativeMomentum[2]; + cascade.positiveMomentum[0] = v0input.positiveMomentum[0]; + cascade.positiveMomentum[1] = v0input.positiveMomentum[1]; + cascade.positiveMomentum[2] = v0input.positiveMomentum[2]; + cascade.negativeMomentum[0] = v0input.negativeMomentum[0]; + cascade.negativeMomentum[1] = v0input.negativeMomentum[1]; + cascade.negativeMomentum[2] = v0input.negativeMomentum[2]; + cascade.v0DaughterDCA = v0input.daughterDCA; + cascade.positiveDCAxy = v0input.positiveDCAxy; + cascade.negativeDCAxy = v0input.negativeDCAxy; + + if (useCascadeMomentumAtPV) { + lCascadeTrack.getPxPyPzGlo(cascade.cascadeMomentum); + } else { + cascade.cascadeMomentum[0] = cascade.bachelorMomentum[0] + v0input.positiveMomentum[0] + v0input.negativeMomentum[0]; + cascade.cascadeMomentum[1] = cascade.bachelorMomentum[1] + v0input.positiveMomentum[1] + v0input.negativeMomentum[1]; + cascade.cascadeMomentum[2] = cascade.bachelorMomentum[2] + v0input.positiveMomentum[2] + v0input.negativeMomentum[2]; + } - cascade.pointingAngle = TMath::ACos(RecoDecay::cpa( - std::array{collision.posX(), collision.posY(), collision.posZ()}, - std::array{cascade.cascadePosition[0], cascade.cascadePosition[1], cascade.cascadePosition[2]}, - std::array{v0input.positiveMomentum[0] + v0input.negativeMomentum[0] + cascade.bachelorMomentum[0], - v0input.positiveMomentum[0] + v0input.negativeMomentum[1] + cascade.bachelorMomentum[1], - v0input.positiveMomentum[2] + v0input.negativeMomentum[2] + cascade.bachelorMomentum[2]})); - - // Calculate DCAxy of the cascade (with bending) - auto lCascadeTrack = fitter.createParentTrackParCov(); - lCascadeTrack.setAbsCharge(cascade.charge); // to be sure - lCascadeTrack.setPID(o2::track::PID::XiMinus); // FIXME: not OK for omegas - dcaInfo[0] = 999; - dcaInfo[1] = 999; - - o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, lCascadeTrack, 2.f, fitter.getMatCorrType(), &dcaInfo); - cascade.cascadeDCAxy = dcaInfo[0]; - cascade.cascadeDCAz = dcaInfo[1]; - - // Calculate masses a priori - cascade.massXi = RecoDecay::m(std::array{std::array{cascade.bachelorMomentum[0], cascade.bachelorMomentum[1], cascade.bachelorMomentum[2]}, - std::array{v0input.positiveMomentum[0] + v0input.negativeMomentum[0], v0input.positiveMomentum[1] + v0input.negativeMomentum[1], v0input.positiveMomentum[2] + v0input.negativeMomentum[2]}}, - std::array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassLambda}); - cascade.massOmega = RecoDecay::m(std::array{std::array{cascade.bachelorMomentum[0], cascade.bachelorMomentum[1], cascade.bachelorMomentum[2]}, - std::array{v0input.positiveMomentum[0] + v0input.negativeMomentum[0], v0input.positiveMomentum[1] + v0input.negativeMomentum[1], v0input.positiveMomentum[2] + v0input.negativeMomentum[2]}}, - std::array{o2::constants::physics::MassKaonCharged, o2::constants::physics::MassLambda}); - - // Populate information - // cascadecandidate.v0Id = v0index.globalIndex(); - cascade.positiveTrack = positiveTrack.globalIndex(); - cascade.negativeTrack = negativeTrack.globalIndex(); - cascade.bachelorTrack = bachelorTrack.globalIndex(); - cascade.positiveTrackX = v0input.positiveTrackX; - cascade.negativeTrackX = v0input.negativeTrackX; - cascade.bachelorTrackX = lBachelorTrack.getX(); // from this minimization - cascade.v0Position[0] = v0input.position[0]; - cascade.v0Position[1] = v0input.position[1]; - cascade.v0Position[2] = v0input.position[2]; - cascade.v0Momentum[0] = v0input.positiveMomentum[0] + v0input.negativeMomentum[0]; - cascade.v0Momentum[1] = v0input.positiveMomentum[1] + v0input.negativeMomentum[1]; - cascade.v0Momentum[2] = v0input.positiveMomentum[2] + v0input.negativeMomentum[2]; - cascade.positiveMomentum[0] = v0input.positiveMomentum[0]; - cascade.positiveMomentum[1] = v0input.positiveMomentum[1]; - cascade.positiveMomentum[2] = v0input.positiveMomentum[2]; - cascade.negativeMomentum[0] = v0input.negativeMomentum[0]; - cascade.negativeMomentum[1] = v0input.negativeMomentum[1]; - cascade.negativeMomentum[2] = v0input.negativeMomentum[2]; - cascade.v0DaughterDCA = v0input.daughterDCA; - cascade.positiveDCAxy = v0input.positiveDCAxy; - cascade.negativeDCAxy = v0input.negativeDCAxy; - - if (useCascadeMomentumAtPV) { - lCascadeTrack.getPxPyPzGlo(cascade.cascadeMomentum); - } else { - cascade.cascadeMomentum[0] = cascade.bachelorMomentum[0] + v0input.positiveMomentum[0] + v0input.negativeMomentum[0]; - cascade.cascadeMomentum[1] = cascade.bachelorMomentum[1] + v0input.positiveMomentum[1] + v0input.negativeMomentum[1]; - cascade.cascadeMomentum[2] = cascade.bachelorMomentum[2] + v0input.positiveMomentum[2] + v0input.negativeMomentum[2]; + if (processCovariances) { + // Calculate position covariance matrix + auto covVtxV = fitter.calcPCACovMatrix(0); + // std::array positionCovariance; + float positionCovariance[6]; + positionCovariance[0] = covVtxV(0, 0); + positionCovariance[1] = covVtxV(1, 0); + positionCovariance[2] = covVtxV(1, 1); + positionCovariance[3] = covVtxV(2, 0); + positionCovariance[4] = covVtxV(2, 1); + positionCovariance[5] = covVtxV(2, 2); + // store momentum covariance matrix + std::array covTv0 = {0.}; + std::array covTbachelor = {0.}; + // std::array momentumCovariance; + lV0Track.getCovXYZPxPyPzGlo(covTv0); + lBachelorTrack.getCovXYZPxPyPzGlo(covTbachelor); + constexpr int MomInd[6] = {9, 13, 14, 18, 19, 20}; // cov matrix elements for momentum component + for (int i = 0; i < 21; i++) { + cascade.covariance[i] = 0.0f; } - - if(processCovariances){ - // Calculate position covariance matrix - auto covVtxV = fitter.calcPCACovMatrix(0); - // std::array positionCovariance; - float positionCovariance[6]; - positionCovariance[0] = covVtxV(0, 0); - positionCovariance[1] = covVtxV(1, 0); - positionCovariance[2] = covVtxV(1, 1); - positionCovariance[3] = covVtxV(2, 0); - positionCovariance[4] = covVtxV(2, 1); - positionCovariance[5] = covVtxV(2, 2); - // store momentum covariance matrix - std::array covTv0 = {0.}; - std::array covTbachelor = {0.}; - // std::array momentumCovariance; - lV0Track.getCovXYZPxPyPzGlo(covTv0); - lBachelorTrack.getCovXYZPxPyPzGlo(covTbachelor); - constexpr int MomInd[6] = {9, 13, 14, 18, 19, 20}; // cov matrix elements for momentum component - for (int i = 0; i < 21; i++) { - cascade.covariance[i] = 0.0f; - } - for (int i = 0; i < 6; i++) { - cascade.covariance[i] = positionCovariance[i]; - cascade.covariance[MomInd[i]] = covTv0[MomInd[i]] + covTbachelor[MomInd[i]]; - } + for (int i = 0; i < 6; i++) { + cascade.covariance[i] = positionCovariance[i]; + cascade.covariance[MomInd[i]] = covTv0[MomInd[i]] + covTbachelor[MomInd[i]]; } - - // Final outcome is YES if I got here! - return true; } - template - bool buildCascadeCandidateWithKF(o2::aod::Collision const& collision, - TTrack const& positiveTrack, - TTrack const& negativeTrack, - TTrack const& bachelorTrack, - bool calculateBachelorBaryonVariables = false, - int kfConstructMethod = 2, - bool kfTuneForOmega = false, - bool kfUseV0MassConstraint = true, - bool kfUseCascadeMassConstraint = false, - bool kfDoDCAFitterPreMinimV0 = false, - bool kfDoDCAFitterPreMinimCasc = false) - { - //*>~<*>~<*>~<*>~<*>~<*>~<*>~<*>~<*>~<* - // KF particle based rebuilding - // dispenses prior V0 generation, uses constrained (re-)fit based on bachelor charge - //*>~<*>~<*>~<*>~<*>~<*>~<*>~<*>~<*>~<* - - if (calculateBachelorBaryonVariables) { - // Calculates properties of the V0 comprised of bachelor and baryon in the cascade - // baryon: distinguished via bachelor charge - if (bachelorTrack.sign() < 0) { - processBachBaryonVariables(collision, bachelorTrack, positiveTrack); - } else { - processBachBaryonVariables(collision, bachelorTrack, negativeTrack); - } - } - - // Overall cascade charge - cascade.charge = bachelorTrack.signed1Pt() > 0 ? +1 : -1; - - // bachelor DCA track to PV - // Calculate DCA with respect to the collision associated to the V0, not individual tracks - gpu::gpustd::array dcaInfo; - - auto bachTrackPar = getTrackPar(bachelorTrack); - o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, bachTrackPar, 2.f, fitter.getMatCorrType(), &dcaInfo); - cascade.bachelorDCAxy = dcaInfo[0]; - o2::track::TrackParCov posTrackParCovForDCA = getTrackParCov(positiveTrack); - o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, posTrackParCovForDCA, 2.f, fitter.getMatCorrType(), &dcaInfo); - cascade.positiveDCAxy = dcaInfo[0]; - o2::track::TrackParCov negTrackParCovForDCA = getTrackParCov(negativeTrack); - o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, negTrackParCovForDCA, 2.f, fitter.getMatCorrType(), &dcaInfo); - cascade.negativeDCAxy = dcaInfo[0]; - - o2::track::TrackParCov lBachelorTrack = getTrackParCov(bachelorTrack); - o2::track::TrackParCov posTrackParCov = getTrackParCov(positiveTrack); - o2::track::TrackParCov negTrackParCov = getTrackParCov(negativeTrack); - - float massPosTrack, massNegTrack; - if (cascade.charge < 0) { - massPosTrack = o2::constants::physics::MassProton; - massNegTrack = o2::constants::physics::MassPionCharged; + // Final outcome is YES if I got here! + return true; + } + + template + bool buildCascadeCandidateWithKF(o2::aod::Collision const& collision, + TTrack const& positiveTrack, + TTrack const& negativeTrack, + TTrack const& bachelorTrack, + bool calculateBachelorBaryonVariables = false, + int kfConstructMethod = 2, + bool kfTuneForOmega = false, + bool kfUseV0MassConstraint = true, + bool kfUseCascadeMassConstraint = false, + bool kfDoDCAFitterPreMinimV0 = false, + bool kfDoDCAFitterPreMinimCasc = false) + { + //*>~<*>~<*>~<*>~<*>~<*>~<*>~<*>~<*>~<* + // KF particle based rebuilding + // dispenses prior V0 generation, uses constrained (re-)fit based on bachelor charge + //*>~<*>~<*>~<*>~<*>~<*>~<*>~<*>~<*>~<* + + if (calculateBachelorBaryonVariables) { + // Calculates properties of the V0 comprised of bachelor and baryon in the cascade + // baryon: distinguished via bachelor charge + if (bachelorTrack.sign() < 0) { + processBachBaryonVariables(collision, bachelorTrack, positiveTrack); } else { - massPosTrack = o2::constants::physics::MassPionCharged; - massNegTrack = o2::constants::physics::MassProton; - } - - //__________________________________________ - //*>~<* step 1 : V0 with dca fitter, uses material corrections implicitly - // This is optional - move close to minima and therefore take material - if (kfDoDCAFitterPreMinimV0) { - int nCand = 0; - try { - nCand = fitter.process(posTrackParCov, negTrackParCov); - } catch (...) { - LOG(error) << "Exception caught in DCA fitter process call!"; - return false; - } - if (nCand == 0) { - return false; - } - // save classical DCA daughters - cascade.v0DaughterDCA = TMath::Sqrt(fitter.getChi2AtPCACandidate()); - - // re-acquire from DCA fitter - posTrackParCov = fitter.getTrack(0); - negTrackParCov = fitter.getTrack(1); - } - - //__________________________________________ - //*>~<* step 2 : V0 with KF - // create KFParticle objects from trackParCovs - KFParticle kfpPos = createKFParticleFromTrackParCov(posTrackParCov, posTrackParCov.getCharge(), massPosTrack); - KFParticle kfpNeg = createKFParticleFromTrackParCov(negTrackParCov, negTrackParCov.getCharge(), massNegTrack); - const KFParticle* V0Daughters[2] = {&kfpPos, &kfpNeg}; - - // construct V0 - KFParticle KFV0; - KFV0.SetConstructMethod(kfConstructMethod); - try { - KFV0.Construct(V0Daughters, 2); - } catch (std::runtime_error& e) { - LOG(debug) << "Failed to construct cascade V0 from daughter tracks: " << e.what(); - return false; - } - - if (kfUseV0MassConstraint) { - KFV0.SetNonlinearMassConstraint(o2::constants::physics::MassLambda); + processBachBaryonVariables(collision, bachelorTrack, negativeTrack); } + } - // V0 constructed, now recovering TrackParCov for dca fitter minimization (with material correction) - KFV0.TransportToDecayVertex(); - o2::track::TrackParCov v0TrackParCov = getTrackParCovFromKFP(KFV0, o2::track::PID::Lambda, 0); - v0TrackParCov.setAbsCharge(0); // to be sure - - //__________________________________________ - //*>~<* step 3 : Cascade with dca fitter (with material corrections) - if (kfDoDCAFitterPreMinimCasc) { - int nCandCascade = 0; - try { - nCandCascade = fitter.process(v0TrackParCov, lBachelorTrack); - } catch (...) { - LOG(error) << "Exception caught in DCA fitter process call!"; - return false; - } - if (nCandCascade == 0) - return false; - - // save classical DCA daughters - // cascadecandidate.dcacascdau = TMath::Sqrt(fitter.getChi2AtPCACandidate()); - // if (cascadecandidate.dcacascdau > dcacascdau) - // return false; - - v0TrackParCov = fitter.getTrack(0); - lBachelorTrack = fitter.getTrack(1); - } + // Overall cascade charge + cascade.charge = bachelorTrack.signed1Pt() > 0 ? +1 : -1; + + // bachelor DCA track to PV + // Calculate DCA with respect to the collision associated to the V0, not individual tracks + gpu::gpustd::array dcaInfo; + + auto bachTrackPar = getTrackPar(bachelorTrack); + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, bachTrackPar, 2.f, fitter.getMatCorrType(), &dcaInfo); + cascade.bachelorDCAxy = dcaInfo[0]; + o2::track::TrackParCov posTrackParCovForDCA = getTrackParCov(positiveTrack); + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, posTrackParCovForDCA, 2.f, fitter.getMatCorrType(), &dcaInfo); + cascade.positiveDCAxy = dcaInfo[0]; + o2::track::TrackParCov negTrackParCovForDCA = getTrackParCov(negativeTrack); + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, negTrackParCovForDCA, 2.f, fitter.getMatCorrType(), &dcaInfo); + cascade.negativeDCAxy = dcaInfo[0]; + + o2::track::TrackParCov lBachelorTrack = getTrackParCov(bachelorTrack); + o2::track::TrackParCov posTrackParCov = getTrackParCov(positiveTrack); + o2::track::TrackParCov negTrackParCov = getTrackParCov(negativeTrack); + + float massPosTrack, massNegTrack; + if (cascade.charge < 0) { + massPosTrack = o2::constants::physics::MassProton; + massNegTrack = o2::constants::physics::MassPionCharged; + } else { + massPosTrack = o2::constants::physics::MassPionCharged; + massNegTrack = o2::constants::physics::MassProton; + } - //__________________________________________ - //*>~<* step 4 : Cascade with KF particle (potentially mass-constrained if asked) - float massBachelorPion = o2::constants::physics::MassPionCharged; - float massBachelorKaon = o2::constants::physics::MassKaonCharged; - - KFParticle kfpV0 = createKFParticleFromTrackParCov(v0TrackParCov, 0, o2::constants::physics::MassLambda); - KFParticle kfpBachPion = createKFParticleFromTrackParCov(lBachelorTrack, cascade.charge, massBachelorPion); - KFParticle kfpBachKaon = createKFParticleFromTrackParCov(lBachelorTrack, cascade.charge, massBachelorKaon); - const KFParticle* XiDaugthers[2] = {&kfpBachPion, &kfpV0}; - const KFParticle* OmegaDaugthers[2] = {&kfpBachKaon, &kfpV0}; - - // construct mother - KFParticle KFXi, KFOmega; - KFXi.SetConstructMethod(kfConstructMethod); - KFOmega.SetConstructMethod(kfConstructMethod); + //__________________________________________ + //*>~<* step 1 : V0 with dca fitter, uses material corrections implicitly + // This is optional - move close to minima and therefore take material + if (kfDoDCAFitterPreMinimV0) { + int nCand = 0; try { - KFXi.Construct(XiDaugthers, 2); - } catch (std::runtime_error& e) { - LOG(debug) << "Failed to construct xi from V0 and bachelor track: " << e.what(); + nCand = fitter.process(posTrackParCov, negTrackParCov); + } catch (...) { + LOG(error) << "Exception caught in DCA fitter process call!"; return false; } - try { - KFOmega.Construct(OmegaDaugthers, 2); - } catch (std::runtime_error& e) { - LOG(debug) << "Failed to construct omega from V0 and bachelor track: " << e.what(); + if (nCand == 0) { return false; } - if (kfUseCascadeMassConstraint) { - // set mass constraint if requested - // WARNING: this is only adequate for decay chains, i.e. XiC -> Xi or OmegaC -> Omega - KFXi.SetNonlinearMassConstraint(o2::constants::physics::MassXiMinus); - KFOmega.SetNonlinearMassConstraint(o2::constants::physics::MassOmegaMinus); - } - KFXi.TransportToDecayVertex(); - KFOmega.TransportToDecayVertex(); + // save classical DCA daughters + cascade.v0DaughterDCA = TMath::Sqrt(fitter.getChi2AtPCACandidate()); - // get DCA of daughters at vertex - cascade.cascadeDaughterDCA = kfpBachPion.GetDistanceFromParticle(kfpV0); + // re-acquire from DCA fitter + posTrackParCov = fitter.getTrack(0); + negTrackParCov = fitter.getTrack(1); + } - //__________________________________________ - //*>~<* step 5 : propagate cascade to primary vertex with material corrections if asked + //__________________________________________ + //*>~<* step 2 : V0 with KF + // create KFParticle objects from trackParCovs + KFParticle kfpPos = createKFParticleFromTrackParCov(posTrackParCov, posTrackParCov.getCharge(), massPosTrack); + KFParticle kfpNeg = createKFParticleFromTrackParCov(negTrackParCov, negTrackParCov.getCharge(), massNegTrack); + const KFParticle* V0Daughters[2] = {&kfpPos, &kfpNeg}; + + // construct V0 + KFParticle KFV0; + KFV0.SetConstructMethod(kfConstructMethod); + try { + KFV0.Construct(V0Daughters, 2); + } catch (std::runtime_error& e) { + LOG(debug) << "Failed to construct cascade V0 from daughter tracks: " << e.what(); + return false; + } - o2::track::TrackParCov lCascadeTrack; - if (!kfTuneForOmega) { - lCascadeTrack = getTrackParCovFromKFP(KFXi, o2::track::PID::XiMinus, cascade.charge); - } else { - lCascadeTrack = getTrackParCovFromKFP(KFOmega, o2::track::PID::OmegaMinus, cascade.charge); - } - dcaInfo[0] = 999; - dcaInfo[1] = 999; - o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, lCascadeTrack, 2.f, fitter.getMatCorrType(), &dcaInfo); - cascade.cascadeDCAxy = dcaInfo[0]; - cascade.cascadeDCAz = dcaInfo[1]; - - //__________________________________________ - //*>~<* step 6 : acquire all parameters for analysis - - // basic indices - cascade.v0Id = -1; - cascade.positiveTrack = positiveTrack.globalIndex(); - cascade.negativeTrack = negativeTrack.globalIndex(); - cascade.bachelorTrack = bachelorTrack.globalIndex(); - - // KF chi2 - cascade.kfV0Chi2 = KFV0.GetChi2(); - cascade.kfCascadeChi2 = KFXi.GetChi2(); - if (kfTuneForOmega) - cascade.kfCascadeChi2 = KFOmega.GetChi2(); - - // Daughter momentum not KF-updated FIXME --> but DCA fitter updated if pre-minimisation is used - lBachelorTrack.getPxPyPzGlo(cascade.bachelorMomentum); - posTrackParCov.getPxPyPzGlo(cascade.positiveMomentum); - negTrackParCov.getPxPyPzGlo(cascade.negativeMomentum); - - // Daughter track position at vertex not KF-updated FIXME --> but DCA fitter updated if pre-minimisation is used - posTrackParCov.getXYZGlo(cascade.positivePosition); - negTrackParCov.getXYZGlo(cascade.negativePosition); - - // mother position information from KF - cascade.v0Position[0] = KFV0.GetX(); - cascade.v0Position[1] = KFV0.GetY(); - cascade.v0Position[2] = KFV0.GetZ(); - - // mother momentumm information from KF - cascade.v0Momentum[0] = KFV0.GetPx(); - cascade.v0Momentum[1] = KFV0.GetPy(); - cascade.v0Momentum[2] = KFV0.GetPz(); - - // Mother position + momentum is KF updated - if (!kfTuneForOmega) { - cascade.cascadePosition[0] = KFXi.GetX(); - cascade.cascadePosition[1] = KFXi.GetY(); - cascade.cascadePosition[2] = KFXi.GetZ(); - cascade.cascadeMomentum[0] = KFXi.GetPx(); - cascade.cascadeMomentum[1] = KFXi.GetPy(); - cascade.cascadeMomentum[2] = KFXi.GetPz(); - } else { - cascade.cascadePosition[0] = KFOmega.GetX(); - cascade.cascadePosition[1] = KFOmega.GetY(); - cascade.cascadePosition[2] = KFOmega.GetZ(); - cascade.cascadeMomentum[0] = KFOmega.GetPx(); - cascade.cascadeMomentum[1] = KFOmega.GetPy(); - cascade.cascadeMomentum[2] = KFOmega.GetPz(); - } + if (kfUseV0MassConstraint) { + KFV0.SetNonlinearMassConstraint(o2::constants::physics::MassLambda); + } - // KF-aware cosPA - cascade.pointingAngle = TMath::ACos(RecoDecay::cpa( - std::array{collision.posX(), collision.posY(), collision.posZ()}, - std::array{cascade.cascadePosition[0], cascade.cascadePosition[1], cascade.cascadePosition[2]}, - std::array{cascade.cascadeMomentum[0], cascade.cascadeMomentum[1], cascade.cascadeMomentum[2]})); - - // Calculate masses a priori - float MLambda, SigmaLambda, MXi, SigmaXi, MOmega, SigmaOmega; - KFV0.GetMass(MLambda, SigmaLambda); - KFXi.GetMass(MXi, SigmaXi); - KFOmega.GetMass(MOmega, SigmaOmega); - cascade.kfMLambda = MLambda; - cascade.massXi = MXi; - cascade.massOmega = MOmega; - - // KF Cascade covariance matrix - o2::gpu::gpustd::array covCascKF; - for (int i = 0; i < 21; i++) { // get covariance matrix elements (lower triangle) - covCascKF[i] = KFXi.GetCovariance(i); - cascade.covariance[i] = covCascKF[i]; - } + // V0 constructed, now recovering TrackParCov for dca fitter minimization (with material correction) + KFV0.TransportToDecayVertex(); + o2::track::TrackParCov v0TrackParCov = getTrackParCovFromKFP(KFV0, o2::track::PID::Lambda, 0); + v0TrackParCov.setAbsCharge(0); // to be sure - // KF V0 covariance matrix - o2::gpu::gpustd::array covV0KF; - for (int i = 0; i < 21; i++) { // get covariance matrix elements (lower triangle) - covV0KF[i] = KFV0.GetCovariance(i); - cascade.kfTrackCovarianceV0[i] = covV0KF[i]; + //__________________________________________ + //*>~<* step 3 : Cascade with dca fitter (with material corrections) + if (kfDoDCAFitterPreMinimCasc) { + int nCandCascade = 0; + try { + nCandCascade = fitter.process(v0TrackParCov, lBachelorTrack); + } catch (...) { + LOG(error) << "Exception caught in DCA fitter process call!"; + return false; } + if (nCandCascade == 0) + return false; - // V0 daughter covariance matrices - std::array cvPosKF, cvNegKF; - posTrackParCov.getCovXYZPxPyPzGlo(cvPosKF); - negTrackParCov.getCovXYZPxPyPzGlo(cvNegKF); - for (int i = 0; i < 21; i++) { - cascade.kfTrackCovariancePos[i] = cvPosKF[i]; - cascade.kfTrackCovarianceNeg[i] = cvNegKF[i]; - } + // save classical DCA daughters + // cascadecandidate.dcacascdau = TMath::Sqrt(fitter.getChi2AtPCACandidate()); + // if (cascadecandidate.dcacascdau > dcacascdau) + // return false; - return true; + v0TrackParCov = fitter.getTrack(0); + lBachelorTrack = fitter.getTrack(1); } - o2::base::MatLayerCylSet* lut; // material LUT for DCA fitter - o2::vertexing::DCAFitterN<2> fitter; // 2-prong o2 dca fitter - - v0candidate v0; // storage for V0 candidate properties - cascadeCandidate cascade; // storage for cascade candidate properties - - private: - // internal helper to calculate DCAxy of a straight line to a given PV analytically - float CalculateDCAStraightToPV(float X, float Y, float Z, float Px, float Py, float Pz, float pvX, float pvY, float pvZ){ - return std::sqrt((std::pow((pvY - Y) * Pz - (pvZ - Z) * Py, 2) + std::pow((pvX - X) * Pz - (pvZ - Z) * Px, 2) + std::pow((pvX - X) * Py - (pvY - Y) * Px, 2)) / (Px * Px + Py * Py + Pz * Pz)); + //__________________________________________ + //*>~<* step 4 : Cascade with KF particle (potentially mass-constrained if asked) + float massBachelorPion = o2::constants::physics::MassPionCharged; + float massBachelorKaon = o2::constants::physics::MassKaonCharged; + + KFParticle kfpV0 = createKFParticleFromTrackParCov(v0TrackParCov, 0, o2::constants::physics::MassLambda); + KFParticle kfpBachPion = createKFParticleFromTrackParCov(lBachelorTrack, cascade.charge, massBachelorPion); + KFParticle kfpBachKaon = createKFParticleFromTrackParCov(lBachelorTrack, cascade.charge, massBachelorKaon); + const KFParticle* XiDaugthers[2] = {&kfpBachPion, &kfpV0}; + const KFParticle* OmegaDaugthers[2] = {&kfpBachKaon, &kfpV0}; + + // construct mother + KFParticle KFXi, KFOmega; + KFXi.SetConstructMethod(kfConstructMethod); + KFOmega.SetConstructMethod(kfConstructMethod); + try { + KFXi.Construct(XiDaugthers, 2); + } catch (std::runtime_error& e) { + LOG(debug) << "Failed to construct xi from V0 and bachelor track: " << e.what(); + return false; } + try { + KFOmega.Construct(OmegaDaugthers, 2); + } catch (std::runtime_error& e) { + LOG(debug) << "Failed to construct omega from V0 and bachelor track: " << e.what(); + return false; + } + if (kfUseCascadeMassConstraint) { + // set mass constraint if requested + // WARNING: this is only adequate for decay chains, i.e. XiC -> Xi or OmegaC -> Omega + KFXi.SetNonlinearMassConstraint(o2::constants::physics::MassXiMinus); + KFOmega.SetNonlinearMassConstraint(o2::constants::physics::MassOmegaMinus); + } + KFXi.TransportToDecayVertex(); + KFOmega.TransportToDecayVertex(); - template - void processBachBaryonVariables(TCollision const& collision, TTrack const& track1, TTrack const& track2) - { - cascade.bachBaryonCosPA = 0; // would ordinarily accept all - cascade.bachBaryonDCAxyToPV = 1e+3; // would ordinarily accept all + // get DCA of daughters at vertex + cascade.cascadeDaughterDCA = kfpBachPion.GetDistanceFromParticle(kfpV0); - // create tracks from table rows - o2::track::TrackParCov tr1 = getTrackParCov(track1); - o2::track::TrackParCov tr2 = getTrackParCov(track2); + //__________________________________________ + //*>~<* step 5 : propagate cascade to primary vertex with material corrections if asked - //---/---/---/ - // Move close to minima - int nCand = 0; - try { - nCand = fitter.process(tr1, tr2); - } catch (...) { - return; - } - if (nCand == 0) - return; // variables are such that candidate is accepted (not obvious...) - - // Calculate DCAxy of the cascade (with bending) - o2::track::TrackPar wrongV0 = fitter.createParentTrackPar(); - wrongV0.setAbsCharge(0); // charge zero - gpu::gpustd::array dcaInfo; - dcaInfo[0] = 999; - dcaInfo[1] = 999; - - // bachelor-baryon DCAxy to PV - o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, wrongV0, 2.f, fitter.getMatCorrType(), &dcaInfo); - cascade.bachBaryonDCAxyToPV = dcaInfo[0]; - - const auto& vtx = fitter.getPCACandidate(); - if (!fitter.isPropagateTracksToVertexDone()) - return; - - std::array tr1p; - std::array tr2p; - - fitter.getTrack(1).getPxPyPzGlo(tr1p); - fitter.getTrack(2).getPxPyPzGlo(tr2p); - - // bachelor-baryon CosPA - cascade.bachBaryonCosPA = RecoDecay::cpa( - std::array{collision.posX(), collision.posY(), collision.posZ()}, - std::array{vtx[0], vtx[1], vtx[2]}, - std::array{tr1p[0] + tr2p[0], tr1p[1] + tr2p[1], tr1p[2] + tr2p[2]}); - - // Potentially also to be considered: bachelor-baryon DCA (between the two tracks) - // to be added here as complementary information in the future + o2::track::TrackParCov lCascadeTrack; + if (!kfTuneForOmega) { + lCascadeTrack = getTrackParCovFromKFP(KFXi, o2::track::PID::XiMinus, cascade.charge); + } else { + lCascadeTrack = getTrackParCovFromKFP(KFOmega, o2::track::PID::OmegaMinus, cascade.charge); + } + dcaInfo[0] = 999; + dcaInfo[1] = 999; + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, lCascadeTrack, 2.f, fitter.getMatCorrType(), &dcaInfo); + cascade.cascadeDCAxy = dcaInfo[0]; + cascade.cascadeDCAz = dcaInfo[1]; + + //__________________________________________ + //*>~<* step 6 : acquire all parameters for analysis + + // basic indices + cascade.v0Id = -1; + cascade.positiveTrack = positiveTrack.globalIndex(); + cascade.negativeTrack = negativeTrack.globalIndex(); + cascade.bachelorTrack = bachelorTrack.globalIndex(); + + // KF chi2 + cascade.kfV0Chi2 = KFV0.GetChi2(); + cascade.kfCascadeChi2 = KFXi.GetChi2(); + if (kfTuneForOmega) + cascade.kfCascadeChi2 = KFOmega.GetChi2(); + + // Daughter momentum not KF-updated FIXME --> but DCA fitter updated if pre-minimisation is used + lBachelorTrack.getPxPyPzGlo(cascade.bachelorMomentum); + posTrackParCov.getPxPyPzGlo(cascade.positiveMomentum); + negTrackParCov.getPxPyPzGlo(cascade.negativeMomentum); + + // Daughter track position at vertex not KF-updated FIXME --> but DCA fitter updated if pre-minimisation is used + posTrackParCov.getXYZGlo(cascade.positivePosition); + negTrackParCov.getXYZGlo(cascade.negativePosition); + + // mother position information from KF + cascade.v0Position[0] = KFV0.GetX(); + cascade.v0Position[1] = KFV0.GetY(); + cascade.v0Position[2] = KFV0.GetZ(); + + // mother momentumm information from KF + cascade.v0Momentum[0] = KFV0.GetPx(); + cascade.v0Momentum[1] = KFV0.GetPy(); + cascade.v0Momentum[2] = KFV0.GetPz(); + + // Mother position + momentum is KF updated + if (!kfTuneForOmega) { + cascade.cascadePosition[0] = KFXi.GetX(); + cascade.cascadePosition[1] = KFXi.GetY(); + cascade.cascadePosition[2] = KFXi.GetZ(); + cascade.cascadeMomentum[0] = KFXi.GetPx(); + cascade.cascadeMomentum[1] = KFXi.GetPy(); + cascade.cascadeMomentum[2] = KFXi.GetPz(); + } else { + cascade.cascadePosition[0] = KFOmega.GetX(); + cascade.cascadePosition[1] = KFOmega.GetY(); + cascade.cascadePosition[2] = KFOmega.GetZ(); + cascade.cascadeMomentum[0] = KFOmega.GetPx(); + cascade.cascadeMomentum[1] = KFOmega.GetPy(); + cascade.cascadeMomentum[2] = KFOmega.GetPz(); } - // TrackParCov to KF converter - // FIXME: could be an utility somewhere else - // from Carolina Reetz (thank you!) - template - KFParticle createKFParticleFromTrackParCov(const o2::track::TrackParametrizationWithError& trackparCov, int charge, float mass) - { - std::array xyz, pxpypz; - float xyzpxpypz[6]; - trackparCov.getPxPyPzGlo(pxpypz); - trackparCov.getXYZGlo(xyz); - for (int i{0}; i < 3; ++i) { - xyzpxpypz[i] = xyz[i]; - xyzpxpypz[i + 3] = pxpypz[i]; - } + // KF-aware cosPA + cascade.pointingAngle = TMath::ACos(RecoDecay::cpa( + std::array{collision.posX(), collision.posY(), collision.posZ()}, + std::array{cascade.cascadePosition[0], cascade.cascadePosition[1], cascade.cascadePosition[2]}, + std::array{cascade.cascadeMomentum[0], cascade.cascadeMomentum[1], cascade.cascadeMomentum[2]})); + + // Calculate masses a priori + float MLambda, SigmaLambda, MXi, SigmaXi, MOmega, SigmaOmega; + KFV0.GetMass(MLambda, SigmaLambda); + KFXi.GetMass(MXi, SigmaXi); + KFOmega.GetMass(MOmega, SigmaOmega); + cascade.kfMLambda = MLambda; + cascade.massXi = MXi; + cascade.massOmega = MOmega; + + // KF Cascade covariance matrix + o2::gpu::gpustd::array covCascKF; + for (int i = 0; i < 21; i++) { // get covariance matrix elements (lower triangle) + covCascKF[i] = KFXi.GetCovariance(i); + cascade.covariance[i] = covCascKF[i]; + } - std::array cv; - try { - trackparCov.getCovXYZPxPyPzGlo(cv); - } catch (std::runtime_error& e) { - LOG(debug) << "Failed to get cov matrix from TrackParCov" << e.what(); - } + // KF V0 covariance matrix + o2::gpu::gpustd::array covV0KF; + for (int i = 0; i < 21; i++) { // get covariance matrix elements (lower triangle) + covV0KF[i] = KFV0.GetCovariance(i); + cascade.kfTrackCovarianceV0[i] = covV0KF[i]; + } - KFParticle kfPart; - float Mini, SigmaMini, M, SigmaM; - kfPart.GetMass(Mini, SigmaMini); - LOG(debug) << "Daughter KFParticle mass before creation: " << Mini << " +- " << SigmaMini; + // V0 daughter covariance matrices + std::array cvPosKF, cvNegKF; + posTrackParCov.getCovXYZPxPyPzGlo(cvPosKF); + negTrackParCov.getCovXYZPxPyPzGlo(cvNegKF); + for (int i = 0; i < 21; i++) { + cascade.kfTrackCovariancePos[i] = cvPosKF[i]; + cascade.kfTrackCovarianceNeg[i] = cvNegKF[i]; + } - try { - kfPart.Create(xyzpxpypz, cv.data(), charge, mass); - } catch (std::runtime_error& e) { - LOG(debug) << "Failed to create KFParticle from daughter TrackParCov" << e.what(); - } + return true; + } + + o2::base::MatLayerCylSet* lut; // material LUT for DCA fitter + o2::vertexing::DCAFitterN<2> fitter; // 2-prong o2 dca fitter + + v0candidate v0; // storage for V0 candidate properties + cascadeCandidate cascade; // storage for cascade candidate properties + + private: + // internal helper to calculate DCAxy of a straight line to a given PV analytically + float CalculateDCAStraightToPV(float X, float Y, float Z, float Px, float Py, float Pz, float pvX, float pvY, float pvZ) + { + return std::sqrt((std::pow((pvY - Y) * Pz - (pvZ - Z) * Py, 2) + std::pow((pvX - X) * Pz - (pvZ - Z) * Px, 2) + std::pow((pvX - X) * Py - (pvY - Y) * Px, 2)) / (Px * Px + Py * Py + Pz * Pz)); + } + + template + void processBachBaryonVariables(TCollision const& collision, TTrack const& track1, TTrack const& track2) + { + cascade.bachBaryonCosPA = 0; // would ordinarily accept all + cascade.bachBaryonDCAxyToPV = 1e+3; // would ordinarily accept all + + // create tracks from table rows + o2::track::TrackParCov tr1 = getTrackParCov(track1); + o2::track::TrackParCov tr2 = getTrackParCov(track2); + + //---/---/---/ + // Move close to minima + int nCand = 0; + try { + nCand = fitter.process(tr1, tr2); + } catch (...) { + return; + } + if (nCand == 0) + return; // variables are such that candidate is accepted (not obvious...) + + // Calculate DCAxy of the cascade (with bending) + o2::track::TrackPar wrongV0 = fitter.createParentTrackPar(); + wrongV0.setAbsCharge(0); // charge zero + gpu::gpustd::array dcaInfo; + dcaInfo[0] = 999; + dcaInfo[1] = 999; + + // bachelor-baryon DCAxy to PV + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, wrongV0, 2.f, fitter.getMatCorrType(), &dcaInfo); + cascade.bachBaryonDCAxyToPV = dcaInfo[0]; + + const auto& vtx = fitter.getPCACandidate(); + if (!fitter.isPropagateTracksToVertexDone()) + return; + + std::array tr1p; + std::array tr2p; + + fitter.getTrack(1).getPxPyPzGlo(tr1p); + fitter.getTrack(2).getPxPyPzGlo(tr2p); + + // bachelor-baryon CosPA + cascade.bachBaryonCosPA = RecoDecay::cpa( + std::array{collision.posX(), collision.posY(), collision.posZ()}, + std::array{vtx[0], vtx[1], vtx[2]}, + std::array{tr1p[0] + tr2p[0], tr1p[1] + tr2p[1], tr1p[2] + tr2p[2]}); + + // Potentially also to be considered: bachelor-baryon DCA (between the two tracks) + // to be added here as complementary information in the future + } + + // TrackParCov to KF converter + // FIXME: could be an utility somewhere else + // from Carolina Reetz (thank you!) + template + KFParticle createKFParticleFromTrackParCov(const o2::track::TrackParametrizationWithError& trackparCov, int charge, float mass) + { + std::array xyz, pxpypz; + float xyzpxpypz[6]; + trackparCov.getPxPyPzGlo(pxpypz); + trackparCov.getXYZGlo(xyz); + for (int i{0}; i < 3; ++i) { + xyzpxpypz[i] = xyz[i]; + xyzpxpypz[i + 3] = pxpypz[i]; + } - kfPart.GetMass(M, SigmaM); - LOG(debug) << "Daughter KFParticle mass after creation: " << M << " +- " << SigmaM; - return kfPart; + std::array cv; + try { + trackparCov.getCovXYZPxPyPzGlo(cv); + } catch (std::runtime_error& e) { + LOG(debug) << "Failed to get cov matrix from TrackParCov" << e.what(); } - // KF to TrackParCov converter - // FIXME: could be an utility somewhere else - // from Carolina Reetz (thank you!) - o2::track::TrackParCov getTrackParCovFromKFP(const KFParticle& kfParticle, const o2::track::PID pid, const int sign) - { - o2::gpu::gpustd::array xyz, pxpypz; - o2::gpu::gpustd::array cv; - - // get parameters from kfParticle - xyz[0] = kfParticle.GetX(); - xyz[1] = kfParticle.GetY(); - xyz[2] = kfParticle.GetZ(); - pxpypz[0] = kfParticle.GetPx(); - pxpypz[1] = kfParticle.GetPy(); - pxpypz[2] = kfParticle.GetPz(); - - // set covariance matrix elements (lower triangle) - for (int i = 0; i < 21; i++) { - cv[i] = kfParticle.GetCovariance(i); - } + KFParticle kfPart; + float Mini, SigmaMini, M, SigmaM; + kfPart.GetMass(Mini, SigmaMini); + LOG(debug) << "Daughter KFParticle mass before creation: " << Mini << " +- " << SigmaMini; - // create TrackParCov track - o2::track::TrackParCov track = o2::track::TrackParCov(xyz, pxpypz, cv, sign, true, pid); - return track; + try { + kfPart.Create(xyzpxpypz, cv.data(), charge, mass); + } catch (std::runtime_error& e) { + LOG(debug) << "Failed to create KFParticle from daughter TrackParCov" << e.what(); } + + kfPart.GetMass(M, SigmaM); + LOG(debug) << "Daughter KFParticle mass after creation: " << M << " +- " << SigmaM; + return kfPart; + } + + // KF to TrackParCov converter + // FIXME: could be an utility somewhere else + // from Carolina Reetz (thank you!) + o2::track::TrackParCov getTrackParCovFromKFP(const KFParticle& kfParticle, const o2::track::PID pid, const int sign) + { + o2::gpu::gpustd::array xyz, pxpypz; + o2::gpu::gpustd::array cv; + + // get parameters from kfParticle + xyz[0] = kfParticle.GetX(); + xyz[1] = kfParticle.GetY(); + xyz[2] = kfParticle.GetZ(); + pxpypz[0] = kfParticle.GetPx(); + pxpypz[1] = kfParticle.GetPy(); + pxpypz[2] = kfParticle.GetPz(); + + // set covariance matrix elements (lower triangle) + for (int i = 0; i < 21; i++) { + cv[i] = kfParticle.GetCovariance(i); + } + + // create TrackParCov track + o2::track::TrackParCov track = o2::track::TrackParCov(xyz, pxpypz, cv, sign, true, pid); + return track; + } }; } // namespace pwglf