diff --git a/PWGLF/Tasks/Resonances/higherMassResonances.cxx b/PWGLF/Tasks/Resonances/higherMassResonances.cxx index 3a34c13fad2..7b13c8e0596 100644 --- a/PWGLF/Tasks/Resonances/higherMassResonances.cxx +++ b/PWGLF/Tasks/Resonances/higherMassResonances.cxx @@ -51,6 +51,7 @@ #include #include #include +#include #include #include @@ -69,50 +70,57 @@ struct HigherMassResonances { HistogramRegistry hglue{"hglueball", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; HistogramRegistry hMChists{"hMChists", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; - struct RCTCut : ConfigurableGroup { + struct : ConfigurableGroup { Configurable requireRCTFlagChecker{"requireRCTFlagChecker", true, "Check event quality in run condition table"}; Configurable cfgEvtRCTFlagCheckerLabel{"cfgEvtRCTFlagCheckerLabel", "CBT_hadronPID", "Evt sel: RCT flag checker label"}; Configurable cfgEvtRCTFlagCheckerZDCCheck{"cfgEvtRCTFlagCheckerZDCCheck", false, "Evt sel: RCT flag checker ZDC check"}; Configurable cfgEvtRCTFlagCheckerLimitAcceptAsBad{"cfgEvtRCTFlagCheckerLimitAcceptAsBad", true, "Evt sel: RCT flag checker treat Limited Acceptance As Bad"}; - - RCTFlagsChecker rctChecker; + } rctCut; + RCTFlagsChecker rctChecker; + + enum MultEstimator { + kFT0M, + kFT0A, + kFT0C, + kFV0A, + kFV0C, + kFV0M, + kNEstimators // useful if you want to iterate or size things }; - RCTCut rctCut; struct : ConfigurableGroup { // PID and QA Configurable qAv0{"qAv0", false, "qAv0"}; Configurable qAPID{"qAPID", true, "qAPID"}; - // Configurable qAv0Daughters{"qAv0Daughters", false, "QA of v0 daughters"}; Configurable qAevents{"qAevents", false, "QA of events"}; - // Configurable invMass1D{"invMass1D", false, "1D invariant mass histograms"}; Configurable correlation2Dhist{"correlation2Dhist", true, "Lamda K0 mass correlation"}; - Configurable cDCAv0topv{"cDCAv0topv", false, "DCA V0 to PV"}; - // Configurable armcut{"armcut", false, "arm cut"}; - Configurable globalTracks{"globalTracks", false, "Global tracks"}; + Configurable isApplyDCAv0topv{"isApplyDCAv0topv", false, "DCA V0 to PV"}; Configurable hasTPC{"hasTPC", false, "TPC"}; - Configurable selectTWOKsOnly{"selectTWOKsOnly", true, "Select only events with two K0s"}; - Configurable applyPairRapidityRec{"applyPairRapidityRec", false, "Apply pair rapidity cut on reconstructed mother (after already applying rapidity cut on generated mother)"}; - Configurable applyPairRapidityGen{"applyPairRapidityGen", false, "Apply pair rapidity cut on generated mother (before applying rapidity cut on reconstructed mother)"}; + Configurable isselectTWOKsOnly{"isselectTWOKsOnly", true, "Select only events with two K0s"}; + Configurable isapplyPairRapidityRec{"isapplyPairRapidityRec", false, "Apply pair rapidity cut on reconstructed mother (after already applying rapidity cut on generated mother)"}; + Configurable isapplyPairRapidityGen{"isapplyPairRapidityGen", false, "Apply pair rapidity cut on generated mother (before applying rapidity cut on reconstructed mother)"}; + Configurable cSelectMultEstimator{"cSelectMultEstimator", 0, "Select multiplicity estimator: 0 - FT0M, 1 - FT0A, 2 - FT0C"}; + Configurable configOccCut{"configOccCut", 1000, "Occupancy cut"}; + Configurable isVertexTOFMatched{"isVertexTOFMatched", false, "Vertex TOF Matched"}; + Configurable isNoCollInTimeRangeStandard{"isNoCollInTimeRangeStandard", false, "No collision in time range standard"}; // Configurables for event selection + // Configurable isINELgt0{"isINELgt0", true, "INEL>0 selection"}; + Configurable isTriggerTVX{"isTriggerTVX", false, "TriggerTVX"}; + Configurable isGoodZvtxFT0vsPV{"isGoodZvtxFT0vsPV", false, "IsGoodZvtxFT0vsPV"}; + Configurable isApplyOccCut{"isApplyOccCut", true, "Apply occupancy cut"}; Configurable cutzvertex{"cutzvertex", 10.0f, "Accepted z-vertex range (cm)"}; - Configurable cfgETAcut{"cfgETAcut", 0.8f, "Track ETA cut"}; Configurable timFrameEvsel{"timFrameEvsel", true, "TPC Time frame boundary cut"}; - // Configurable piluprejection{"piluprejection", false, "Pileup rejection"}; - // Configurable goodzvertex{"goodzvertex", false, "removes collisions with large differences between z of PV by tracks and z of PV from FT0 A-C time difference."}; - // Configurable itstpctracks{"itstpctracks", false, "selects collisions with at least one ITS-TPC track,"}; - // Configurable additionalEvsel{"additionalEvsel", false, "Additional event selcection"}; - // Configurable applyOccupancyCut{"applyOccupancyCut", false, "Apply occupancy cut"}; - // Configurable occupancyCut{"occupancyCut", 1000, "Mimimum Occupancy cut"}; + Configurable isNoSameBunchPileup{"isNoSameBunchPileup", true, "kNoSameBunchPileup"}; + Configurable isAllLayersGoodITS{"isAllLayersGoodITS", true, "Require all ITS layers to be good"}; + Configurable isNoTimeFrameBorder{"isNoTimeFrameBorder", true, "kNoTimeFrameBorder"}; + Configurable isNoITSROFrameBorder{"isNoITSROFrameBorder", true, "kNoITSROFrameBorder"}; // Configurable parameters for V0 selection Configurable confV0DCADaughMax{"confV0DCADaughMax", 1.0f, "DCA b/w V0 daughters"}; Configurable v0settingDcapostopv{"v0settingDcapostopv", 0.06, "DCA Pos To PV"}; Configurable v0settingDcanegtopv{"v0settingDcanegtopv", 0.06, "DCA Neg To PV"}; - Configurable cMaxV0DCA{"cMaxV0DCA", 1, "DCA V0 to PV"}; - // Configurable isStandarv0{"isStandarv0", false, "Standard V0"}; - // Configurable ConfDaughDCAMin{"ConfDaughDCAMin", 0.06f, "V0 Daugh sel: Max. DCA Daugh to PV (cm)"}; // same as DCA pos to pv and neg to pv + Configurable cMaxV0DCA{"cMaxV0DCA", 0.3, "DCA V0 to PV"}; Configurable confV0PtMin{"confV0PtMin", 0.f, "Minimum transverse momentum of V0"}; Configurable confV0CPAMin{"confV0CPAMin", 0.97f, "Minimum CPA of V0"}; Configurable confV0TranRadV0Min{"confV0TranRadV0Min", 0.5f, "Minimum transverse radius"}; @@ -122,27 +130,27 @@ struct HigherMassResonances { Configurable cWidthKs0{"cWidthKs0", 0.005, "Width of KS0"}; Configurable confDaughEta{"confDaughEta", 0.8f, "V0 Daugh sel: max eta"}; Configurable confDaughTPCnclsMin{"confDaughTPCnclsMin", 70.f, "V0 Daugh sel: Min. nCls TPC"}; - Configurable confDaughPIDCuts{"confDaughPIDCuts", 5, "PID selections for KS0 daughters"}; - // Configurable confarmcut{"confarmcut", 0.2f, "Armenteros cut"}; + Configurable confDaughPIDCutTPC{"confDaughPIDCutTPC", 5, "PID selections for KS0 daughters"}; + Configurable confDaughPIDCutTOF{"confDaughPIDCutTOF", 5, "PID selections for KS0 daughters in TOF"}; Configurable confKsrapidity{"confKsrapidity", 0.5f, "Rapidity cut on K0s"}; - // Configurable lowmasscutks0{"lowmasscutks0", 0.497 - 4 * 0.005, "Low mass cut on K0s"}; - // Configurable highmasscutks0{"highmasscutks0", 0.497 + 4 * 0.005, "High mass cut on K0s"}; - Configurable applyAngSepCut{"applyAngSepCut", false, "Apply angular separation cut"}; Configurable angSepCut{"angSepCut", 0.01f, "Angular separation cut"}; + Configurable isapplyAngSepCut{"isapplyAngSepCut", false, "Apply angular separation cut"}; + Configurable isStandardV0{"isStandardV0", false, "Standard V0 selection"}; + Configurable isApplyEtaCutK0s{"isApplyEtaCutK0s", false, "Apply eta cut on K0s daughters"}; + Configurable cfgETAcut{"cfgETAcut", 0.8f, "Track ETA cut"}; // Configurable for track selection and multiplicity Configurable cfgPTcut{"cfgPTcut", 0.2f, "Track PT cut"}; Configurable cfgNmixedEvents{"cfgNmixedEvents", 5, "Number of mixed events"}; - Configurable cfgMultFOTM{"cfgMultFOTM", true, "Use FOTM multiplicity if pp else use 0 here for PbPb (FT0C)"}; ConfigurableAxis binsCent{"binsCent", {VARIABLE_WIDTH, 0., 5., 10., 30., 50., 70., 100., 110., 150.}, "Binning of the centrality axis"}; // Configurable for MC Configurable isMC{"isMC", false, "Is MC"}; - Configurable allGenCollisions{"allGenCollisions", true, "To fill all generated collisions for the signal loss calculations"}; - Configurable cTVXEvsel{"cTVXEvsel", true, "Triggger selection"}; - Configurable avoidsplitrackMC{"avoidsplitrackMC", false, "avoid split track in MC"}; + Configurable isallGenCollisions{"isallGenCollisions", true, "To fill all generated collisions for the signal loss calculations"}; + Configurable iscTVXEvsel{"iscTVXEvsel", true, "Triggger selection"}; + Configurable isavoidsplitrackMC{"isavoidsplitrackMC", false, "avoid split track in MC"}; + Configurable isapplyRapidityMC{"isapplyRapidityMC", true, "Apply rapidity cut on generated and reconstructed particles"}; Configurable selectMCparticles{"selectMCparticles", 1, "0: f0(1710), 1: f2(1525), 2: a2(1320), 3: f0(1370), 4: f0(1500)"}; - Configurable applyRapidityMC{"applyRapidityMC", true, "Apply rapidity cut on generated and reconstructed particles"}; std::vector pdgCodes = {10331, 335, 115, 10221, 9030221}; // output THnSparses @@ -153,10 +161,8 @@ struct HigherMassResonances { Configurable cRotations{"cRotations", 3, "Number of random rotations in the rotational background"}; // Other cuts on Ks and glueball - // Configurable rapidityks{"rapidityks", true, "rapidity cut on K0s"}; - Configurable applyCompetingcut{"applyCompetingcut", false, "Competing cascade rejection cut"}; + Configurable isapplyCompetingcut{"isapplyCompetingcut", false, "Competing cascade rejection cut"}; Configurable competingcascrejlambda{"competingcascrejlambda", 0.005, "rejecting competing cascade lambda"}; - // Configurable competingcascrejlambdaanti{"competingcascrejlambdaanti", 0.005, "rejecting competing cascade anti-lambda"}; // If one of the pions is misidentified as a proton, then instead of Ks we reconstruct lambda, therefore the competing cascade rejection cut is applied in which if the reconstrcted mass of a pion and proton (which we are assuming to be misidentified as proton) is close to lambda or anti-lambda, then the track is rejected Configurable tpcCrossedrows{"tpcCrossedrows", 70, "TPC crossed rows"}; Configurable tpcCrossedrowsOverfcls{"tpcCrossedrowsOverfcls", 0.8, "TPC crossed rows over findable clusters"}; @@ -167,9 +173,6 @@ struct HigherMassResonances { ConfigurableAxis ksMassBins{"ksMassBins", {200, 0.45f, 0.55f}, "K0s invariant mass axis"}; ConfigurableAxis cGlueMassBins{"cGlueMassBins", {200, 0.9f, 3.0f}, "Glueball invariant mass axis"}; ConfigurableAxis cPtBins{"cPtBins", {200, 0.0f, 20.0f}, "Glueball pT axis"}; - // ConfigurableAxis axisdEdx{"axisdEdx", {20000, 0.0f, 200.0f}, "dE/dx (a.u.)"}; - // ConfigurableAxis axisPtfordEbydx{"axisPtfordEbydx", {2000, 0, 20}, "pT (GeV/c)"}; - // ConfigurableAxis axisMultdist{"axisMultdist", {3500, 0, 70000}, "Multiplicity distribution"}; // fixed variables float rapidityMotherData = 0.5; @@ -199,10 +202,7 @@ struct HigherMassResonances { void init(InitContext const&) { - rctCut.rctChecker.init( - rctCut.cfgEvtRCTFlagCheckerLabel, - rctCut.cfgEvtRCTFlagCheckerZDCCheck, - rctCut.cfgEvtRCTFlagCheckerLimitAcceptAsBad); + rctChecker.init(rctCut.cfgEvtRCTFlagCheckerLabel, rctCut.cfgEvtRCTFlagCheckerZDCCheck, rctCut.cfgEvtRCTFlagCheckerLimitAcceptAsBad); // Axes AxisSpec k0ShortMassAxis = {config.ksMassBins, "#it{M}_{inv} [GeV/#it{c}^{2}]"}; @@ -216,8 +216,6 @@ struct HigherMassResonances { // THnSparses std::array sparses = {config.activateTHnSparseCosThStarHelicity, config.activateTHnSparseCosThStarProduction, config.activateTHnSparseCosThStarBeam, config.activateTHnSparseCosThStarRandom}; - // std::array sparses = {config.activateTHnSparseCosThStarHelicity}; - if (std::accumulate(sparses.begin(), sparses.end(), 0) == 0) { LOGP(fatal, "No output THnSparses enabled"); } else { @@ -239,31 +237,71 @@ struct HigherMassResonances { if (config.qAevents) { rEventSelection.add("hVertexZRec", "hVertexZRec", {HistType::kTH1F, {vertexZAxis}}); rEventSelection.add("hmultiplicity", "multiplicity percentile distribution", {HistType::kTH1F, {{150, 0.0f, 150.0f}}}); - // rEventSelection.add("multdist_FT0M", "FT0M Multiplicity distribution", kTH1F, {config.axisMultdist}); - // rEventSelection.add("multdist_FT0A", "FT0A Multiplicity distribution", kTH1F, {config.axisMultdist}); - // rEventSelection.add("multdist_FT0C", "FT0C Multiplicity distribution", kTH1F, {config.axisMultdist}); - // rEventSelection.add("hNcontributor", "Number of primary vertex contributor", kTH1F, {{2000, 0.0f, 10000.0f}}); + hglue.add("htrackscheck_v0", "htrackscheck_v0", kTH1I, {{15, 0, 15}}); + hglue.add("htrackscheck_v0_daughters", "htrackscheck_v0_daughters", kTH1I, {{15, 0, 15}}); + hMChists.add("events_check", "No. of events in the generated MC", kTH1I, {{20, 0, 20}}); + hMChists.add("events_checkrec", "No. of events in the reconstructed MC", kTH1I, {{25, 0, 25}}); + + rEventSelection.add("hEventCut", "No. of event after cuts", kTH1I, {{20, 0, 20}}); + std::shared_ptr hCutFlow = rEventSelection.get(HIST("hEventCut")); + hCutFlow->GetXaxis()->SetBinLabel(1, "All Events"); + hCutFlow->GetXaxis()->SetBinLabel(2, "|Vz| < cut"); + hCutFlow->GetXaxis()->SetBinLabel(3, "sel8"); + hCutFlow->GetXaxis()->SetBinLabel(4, "kNoTimeFrameBorder"); + hCutFlow->GetXaxis()->SetBinLabel(5, "kNoITSROFrameBorder"); + hCutFlow->GetXaxis()->SetBinLabel(6, "kNoSameBunchPileup"); + hCutFlow->GetXaxis()->SetBinLabel(7, "kIsGoodITSLayersAll"); + hCutFlow->GetXaxis()->SetBinLabel(8, "Occupancy Cut"); + hCutFlow->GetXaxis()->SetBinLabel(9, "rctChecker"); + hCutFlow->GetXaxis()->SetBinLabel(10, "kIsTriggerTVX"); + hCutFlow->GetXaxis()->SetBinLabel(11, "kIsGoodZvtxFT0vsPV"); + hCutFlow->GetXaxis()->SetBinLabel(12, "IsINELgt0"); + hCutFlow->GetXaxis()->SetBinLabel(13, "isVertexITSTPC"); + hCutFlow->GetXaxis()->SetBinLabel(14, "isVertexTOFMatched"); + + std::shared_ptr hv0label = rEventSelection.get(HIST("htrackscheck_v0")); + hv0label->GetXaxis()->SetBinLabel(1, "All Tracks"); + hv0label->GetXaxis()->SetBinLabel(2, "DCA V0 to PV"); + hv0label->GetXaxis()->SetBinLabel(3, "y K0s"); + hv0label->GetXaxis()->SetBinLabel(4, "Min V0 pT"); + hv0label->GetXaxis()->SetBinLabel(5, "Daughter DCA"); + hv0label->GetXaxis()->SetBinLabel(6, "CosPA"); + hv0label->GetXaxis()->SetBinLabel(7, "Decay Radius"); + hv0label->GetXaxis()->SetBinLabel(8, "Lifetime"); + hv0label->GetXaxis()->SetBinLabel(9, "CompetingCascade"); + hv0label->GetXaxis()->SetBinLabel(10, "Standard V0"); + hv0label->GetXaxis()->SetBinLabel(11, "Mass Tolerance"); + + std::shared_ptr hv0DauLabel = rEventSelection.get(HIST("htrackscheck_v0_daughters")); + hv0DauLabel->GetXaxis()->SetBinLabel(1, "AllDau Tracks"); + hv0DauLabel->GetXaxis()->SetBinLabel(2, "has TPC"); + hv0DauLabel->GetXaxis()->SetBinLabel(3, "TPC CrossedRows"); + hv0DauLabel->GetXaxis()->SetBinLabel(4, "TPC CRFC"); + hv0DauLabel->GetXaxis()->SetBinLabel(5, "TPC Chi2NCL"); + hv0DauLabel->GetXaxis()->SetBinLabel(6, "Charge"); + hv0DauLabel->GetXaxis()->SetBinLabel(7, "Eta"); + hv0DauLabel->GetXaxis()->SetBinLabel(8, "PID TPC"); + + std::shared_ptr hv0labelmcrec = rEventSelection.get(HIST("events_checkrec")); + hv0labelmcrec->GetXaxis()->SetBinLabel(1, "All Tracks"); + hv0labelmcrec->GetXaxis()->SetBinLabel(2, "V0Daughter Sel."); + hv0labelmcrec->GetXaxis()->SetBinLabel(3, "V0 Sel."); + hv0labelmcrec->GetXaxis()->SetBinLabel(4, "V0 PDG"); + hv0labelmcrec->GetXaxis()->SetBinLabel(5, "Mother PDG"); + hv0labelmcrec->GetXaxis()->SetBinLabel(6, "Same Mother"); + hv0labelmcrec->GetXaxis()->SetBinLabel(7, "Split Track"); + hv0labelmcrec->GetXaxis()->SetBinLabel(8, "Global Index"); + hv0labelmcrec->GetXaxis()->SetBinLabel(9, "Generator"); + hv0labelmcrec->GetXaxis()->SetBinLabel(10, "Rapidity"); } hglue.add("h3glueInvMassDS", "h3glueInvMassDS", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL, thnAxisPhi}, true); hglue.add("h3glueInvMassME", "h3glueInvMassME", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL, thnAxisPhi}, true); hglue.add("h3glueInvMassRot", "h3glueInvMassRot", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL, thnAxisPhi}, true); - hglue.add("heventscheck", "heventscheck", kTH1I, {{10, 0, 10}}); - hglue.add("htrackscheck_v0", "htrackscheck_v0", kTH1I, {{15, 0, 15}}); - hglue.add("htrackscheck_v0_daughters", "htrackscheck_v0_daughters", kTH1I, {{15, 0, 15}}); // K0s topological/PID cuts if (config.correlation2Dhist) { rKzeroShort.add("mass_lambda_kshort_before", "mass under lambda hypotheses and Kshort mass", kTH2F, {{100, 0.2, 0.8}, {100, 0.9, 1.5}}); - // rKzeroShort.add("mass_lambda_kshort_after1", "mass under lambda hypotheses and Kshort mass", kTH2F, {{100, 0.2, 0.8}, {100, 0.9, 1.5}}); - // rKzeroShort.add("mass_lambda_kshort_after2", "mass under lambda hypotheses and Kshort mass", kTH2F, {{100, 0.2, 0.8}, {100, 0.9, 1.5}}); - // rKzeroShort.add("mass_lambda_kshort_after3", "mass under lambda hypotheses and Kshort mass", kTH2F, {{100, 0.2, 0.8}, {100, 0.9, 1.5}}); - // rKzeroShort.add("mass_lambda_kshort_after4", "mass under lambda hypotheses and Kshort mass", kTH2F, {{100, 0.2, 0.8}, {100, 0.9, 1.5}}); - // rKzeroShort.add("mass_lambda_kshort_after5", "mass under lambda hypotheses and Kshort mass", kTH2F, {{100, 0.2, 0.8}, {100, 0.9, 1.5}}); - // rKzeroShort.add("mass_lambda_kshort_after6", "mass under lambda hypotheses and Kshort mass", kTH2F, {{100, 0.2, 0.8}, {100, 0.9, 1.5}}); - // rKzeroShort.add("mass_lambda_kshort_after7", "mass under lambda hypotheses and Kshort mass", kTH2F, {{100, 0.2, 0.8}, {100, 0.9, 1.5}}); - // rKzeroShort.add("mass_lambda_kshort_after8", "mass under lambda hypotheses and Kshort mass", kTH2F, {{100, 0.2, 0.8}, {100, 0.9, 1.5}}); - // rKzeroShort.add("mass_lambda_kshort_after9", "mass under lambda hypotheses and Kshort mass", kTH2F, {{100, 0.2, 0.8}, {100, 0.9, 1.5}}); rKzeroShort.add("mass_lambda_kshort_after10", "mass under lambda hypotheses and Kshort mass", kTH2F, {{100, 0.2, 0.8}, {100, 0.9, 1.5}}); } if (config.qAv0) { @@ -276,25 +314,6 @@ struct HigherMassResonances { rKzeroShort.add("hV0CosPA", "hV0CosPA", {HistType::kTH1F, {{100, 0.96f, 1.1f}}}); rKzeroShort.add("hLT", "hLT", {HistType::kTH1F, {{100, 0.0f, 50.0f}}}); rKzeroShort.add("angularSeparation", "Angular distribution between two K0s vs pT", {HistType::kTH1F, {{200, 0.0f, 4.0f}}}); - // rKzeroShort.add("Mass_lambda", "Mass under lambda hypothesis", kTH1F, {glueballMassAxis}); - // rKzeroShort.add("mass_AntiLambda", "Mass under anti-lambda hypothesis", kTH1F, {glueballMassAxis}); - // rKzeroShort.add("mass_Gamma", "Mass under Gamma hypothesis", kTH1F, {glueballMassAxis}); - - // rKzeroShort.add("mass_Hypertriton", "Mass under hypertriton hypothesis", kTH1F, {glueballMassAxis}); - // rKzeroShort.add("mass_AnitHypertriton", "Mass under anti-hypertriton hypothesis", kTH1F, {glueballMassAxis}); - // rKzeroShort.add("rapidity", "Rapidity distribution", kTH1F, {{100, -1.0f, 1.0f}}); - // rKzeroShort.add("hv0radius", "hv0radius", kTH1F, {{100, 0.0f, 200.0f}}); - // rKzeroShort.add("hDCApostopv", "DCA positive daughter to PV", kTH1F, {{1000, -10.0f, 10.0f}}); - // rKzeroShort.add("hDCAnegtopv", "DCA negative daughter to PV", kTH1F, {{1000, -10.0f, 10.0f}}); - // rKzeroShort.add("hcDCAv0topv", "DCA V0 to PV", kTH1F, {{60, -3.0f, 3.0f}}); - // rKzeroShort.add("halpha", "Armenteros alpha", kTH1F, {{100, -5.0f, 5.0f}}); - // rKzeroShort.add("hqtarmbyalpha", "qtarm/alpha", kTH1F, {{100, 0.0f, 1.0f}}); - // rKzeroShort.add("hpsipair", "psi pair angle", kTH1F, {{100, -5.0f, 5.0f}}); - - // // Topological histograms (before the selection) - // rKzeroShort.add("hDCAV0Daughters_before", "DCA between v0 daughters before the selection", {HistType::kTH1F, {{60, -3.0f, 3.0f}}}); - // rKzeroShort.add("hV0CosPA_before", "hV0CosPA_before", {HistType::kTH1F, {{200, 0.91f, 1.1f}}}); - // rKzeroShort.add("hLT_before", "hLT_before", {HistType::kTH1F, {{100, 0.0f, 50.0f}}}); } rKzeroShort.add("NksProduced", "Number of K0s produced", kTH1I, {{15, -0.5, 14.5}}); @@ -305,29 +324,16 @@ struct HigherMassResonances { rKzeroShort.add("hNSigmaNegPionK0s_after", "hNSigmaNegPionK0s_after", {HistType::kTH2F, {{ptAxis}, {100, -5.f, 5.f}}}); // rKzeroShort.add("dE_by_dx_TPC", "dE/dx signal in the TPC as a function of pT", kTH2F, {config.axisPtfordEbydx, config.axisdEdx}); } - // if (config.qAv0Daughters) { - // rKzeroShort.add("negative_pt", "Negative daughter pT", kTH1F, {ptAxis}); - // rKzeroShort.add("positive_pt", "Positive daughter pT", kTH1F, {ptAxis}); - // rKzeroShort.add("negative_eta", "Negative daughter eta", kTH1F, {{100, -1.0f, 1.0f}}); - // rKzeroShort.add("positive_eta", "Positive daughter eta", kTH1F, {{100, -1.0f, 1.0f}}); - // rKzeroShort.add("negative_phi", "Negative daughter phi", kTH1F, {{70, 0.0f, 7.0f}}); - // rKzeroShort.add("positive_phi", "Positive daughter phi", kTH1F, {{70, 0.0f, 7.0f}}); - // } // For MC if (config.isMC) { - hMChists.add("events_check", "No. of events in the generated MC", kTH1I, {{20, 0, 20}}); - hMChists.add("events_checkrec", "No. of events in the reconstructed MC", kTH1I, {{25, 0, 25}}); hMChists.add("Genf1710", "Gen f_{0}(1710)", kTHnSparseF, {multiplicityAxis, ptAxis, thnAxisPOL}); hMChists.add("Genf17102", "Gen f_{0}(1710)", kTHnSparseF, {multiplicityAxis, ptAxis, thnAxisPOL}); hMChists.add("Recf1710_pt1", "Rec f_{0}(1710) p_{T}", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL}); hMChists.add("Recf1710_pt2", "Rec f_{0}(1710) p_{T}", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL}); - // hMChists.add("Recf1710_p", "Rec f_{0}(1710) p", kTH1F, {ptAxis}); hMChists.add("h1Recsplit", "Rec p_{T}2", kTH1F, {ptAxis}); - // hMChists.add("Recf1710_mass", "Rec f_{0}(1710) mass", kTH1F, {glueballMassAxis}); hMChists.add("Genf1710_mass", "Gen f_{0}(1710) mass", kTH1F, {glueballMassAxis}); hMChists.add("Genf1710_mass2", "Gen f_{0}(1710) mass", kTH1F, {glueballMassAxis}); - // hMChists.add("Genf1710_pt2", "Gen f_{0}(1710) p_{T}", kTH1F, {ptAxis}); hMChists.add("GenPhi", "Gen Phi", kTH1F, {{70, 0.0, 7.0f}}); hMChists.add("GenPhi2", "Gen Phi", kTH1F, {{70, 0.0, 7.0f}}); hMChists.add("GenEta", "Gen Eta", kTHnSparseF, {{150, -1.5f, 1.5f}}); @@ -342,41 +348,85 @@ struct HigherMassResonances { hMChists.add("RecRapidity2", "Rec Rapidity", kTH1F, {{100, -1.0f, 1.0f}}); hMChists.add("Rec_Multiplicity", "Multiplicity in MC", kTH1F, {multiplicityAxis}); hMChists.add("MC_mult_after_event_sel", "Multiplicity in MC", kTH1F, {multiplicityAxis}); - // hMChists.add("GenPx", "Gen Px", kTH1F, {{100, -10.0f, 10.0f}}); - // hMChists.add("GenPy", "Gen Py", kTH1F, {{100, -10.0f, 10.0f}}); - // hMChists.add("GenPz", "Gen Pz", kTH1F, {{100, -10.0f, 10.0f}}); } } - template - bool eventselection(Collision const& collision) + template + bool selectionEvent(const Coll& collision, bool fillHist = true) { - hglue.fill(HIST("heventscheck"), 1.5); + if (fillHist) + rEventSelection.fill(HIST("hEventCut"), 0); - if (config.timFrameEvsel && (!collision.selection_bit(aod::evsel::kNoTimeFrameBorder) || !collision.selection_bit(aod::evsel::kNoITSROFrameBorder))) { + if (std::abs(collision.posZ()) > config.cutzvertex) return false; - } - hglue.fill(HIST("heventscheck"), 2.5); + if (fillHist) + rEventSelection.fill(HIST("hEventCut"), 1); - if (!collision.sel8()) { + if (!collision.sel8()) return false; - } - hglue.fill(HIST("heventscheck"), 3.5); + if (fillHist) + rEventSelection.fill(HIST("hEventCut"), 2); + + if (config.isNoTimeFrameBorder && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) + return false; + if (fillHist) + rEventSelection.fill(HIST("hEventCut"), 3); + + if (config.isNoITSROFrameBorder && !collision.selection_bit(aod::evsel::kNoITSROFrameBorder)) + return false; + if (fillHist) + rEventSelection.fill(HIST("hEventCut"), 4); + + if (config.isNoSameBunchPileup && (!collision.selection_bit(aod::evsel::kNoSameBunchPileup))) + return false; + if (fillHist) + rEventSelection.fill(HIST("hEventCut"), 5); + + if (config.isAllLayersGoodITS && !collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) + return false; + if (fillHist) + rEventSelection.fill(HIST("hEventCut"), 6); - // if (config.piluprejection && !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { + // if (config.isNoCollInTimeRangeStandard && (!collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard))) // return false; - // } - // hglue.fill(HIST("heventscheck"), 4.5); - // if (config.goodzvertex && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { + if (config.isApplyOccCut && (std::abs(collision.trackOccupancyInTimeRange()) > config.configOccCut)) + return false; + if (fillHist) + rEventSelection.fill(HIST("hEventCut"), 7); + + if (rctCut.requireRCTFlagChecker && !rctChecker(collision)) + return false; + if (fillHist) + rEventSelection.fill(HIST("hEventCut"), 8); + + if (config.isTriggerTVX && !collision.selection_bit(aod::evsel::kIsTriggerTVX)) + return false; + if (fillHist) + rEventSelection.fill(HIST("hEventCut"), 9); + + if (config.isGoodZvtxFT0vsPV && !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) + return false; + if (fillHist) + rEventSelection.fill(HIST("hEventCut"), 10); + + // if (config.isINELgt0 && !collision.isInelGt0()) { // return false; // } - // hglue.fill(HIST("heventscheck"), 5.5); + // if (fillHist) + rEventSelection.fill(HIST("hEventCut"), 11); - // if (config.itstpctracks && !collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) { + // if (config.isVertexITSTPC && !collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) { // return false; // } - // hglue.fill(HIST("heventscheck"), 6.5); + if (fillHist) + rEventSelection.fill(HIST("hEventCut"), 12); + + if (config.isVertexTOFMatched && !collision.selection_bit(aod::evsel::kIsVertexTOFmatched)) { + return false; + } + if (fillHist) + rEventSelection.fill(HIST("hEventCut"), 13); return true; } @@ -395,124 +445,89 @@ struct HigherMassResonances { float ctauK0s = candidate.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * o2::constants::physics::MassK0Short; float lowmasscutks0 = o2::constants::physics::MassKPlus - config.cWidthKs0 * config.cSigmaMassKs0; float highmasscutks0 = o2::constants::physics::MassKPlus + config.cWidthKs0 * config.cSigmaMassKs0; - // float decayLength = candidate.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * RecoDecay::sqrtSumOfSquares(candidate.px(), candidate.py(), candidate.pz()); if (config.qAv0) { rKzeroShort.fill(HIST("hMassK0Shortbefore"), candidate.mK0Short(), candidate.pt()); rKzeroShort.fill(HIST("hLT"), ctauK0s); rKzeroShort.fill(HIST("hDCAV0Daughters"), candidate.dcaV0daughters()); rKzeroShort.fill(HIST("hV0CosPA"), candidate.v0cosPA()); - // rKzeroShort.fill(HIST("Mass_lambda"), candidate.mLambda()); - // rKzeroShort.fill(HIST("mass_AntiLambda"), candidate.mAntiLambda()); - // rKzeroShort.fill(HIST("mass_Gamma"), candidate.mGamma()); - // rKzeroShort.fill(HIST("mass_Hypertriton"), candidate.mHypertriton()); - // rKzeroShort.fill(HIST("mass_AnitHypertriton"), candidate.mAntiHypertriton()); - // rKzeroShort.fill(HIST("rapidity"), candidate.yK0Short()); - // rKzeroShort.fill(HIST("hv0radius"), candidate.v0radius()); - // rKzeroShort.fill(HIST("hDCApostopv"), candidate.dcapostopv()); - // rKzeroShort.fill(HIST("hDCAnegtopv"), candidate.dcanegtopv()); - // rKzeroShort.fill(HIST("hcDCAv0topv"), candidate.dcav0topv()); - // rKzeroShort.fill(HIST("halpha"), candidate.alpha()); - // rKzeroShort.fill(HIST("hqtarmbyalpha"), arm); - // rKzeroShort.fill(HIST("hpsipair"), candidate.psipair()); } if (config.correlation2Dhist) rKzeroShort.fill(HIST("mass_lambda_kshort_before"), candidate.mK0Short(), candidate.mLambda()); hglue.fill(HIST("htrackscheck_v0"), 0.5); - if (config.cDCAv0topv && std::fabs(candidate.dcav0topv()) > config.cMaxV0DCA) { + if (config.isApplyDCAv0topv && std::fabs(candidate.dcav0topv()) > config.cMaxV0DCA) { return false; } hglue.fill(HIST("htrackscheck_v0"), 1.5); - // if (config.correlation2Dhist) - // rKzeroShort.fill(HIST("mass_lambda_kshort_after1"), candidate.mK0Short(), candidate.mLambda()); - // if (config.rapidityks && std::abs(candidate.yK0Short()) >= config.confKsrapidity) { - // return false; - // } - if (std::abs(candidate.yK0Short()) >= config.confKsrapidity) { + if (std::abs(candidate.rapidity(0)) >= config.confKsrapidity) { return false; } hglue.fill(HIST("htrackscheck_v0"), 2.5); - // if (config.correlation2Dhist) - // rKzeroShort.fill(HIST("mass_lambda_kshort_after2"), candidate.mK0Short(), candidate.mLambda()); if (pT < config.confV0PtMin) { return false; } hglue.fill(HIST("htrackscheck_v0"), 3.5); - // if (config.correlation2Dhist) - // rKzeroShort.fill(HIST("mass_lambda_kshort_after3"), candidate.mK0Short(), candidate.mLambda()); if (dcaDaughv0 > config.confV0DCADaughMax) { return false; } hglue.fill(HIST("htrackscheck_v0"), 4.5); - // if (config.correlation2Dhist) - // rKzeroShort.fill(HIST("mass_lambda_kshort_after4"), candidate.mK0Short(), candidate.mLambda()); if (cpav0 < config.confV0CPAMin) { return false; } hglue.fill(HIST("htrackscheck_v0"), 5.5); - // if (config.correlation2Dhist) - // rKzeroShort.fill(HIST("mass_lambda_kshort_after5"), candidate.mK0Short(), candidate.mLambda()); if (tranRad < config.confV0TranRadV0Min) { return false; } hglue.fill(HIST("htrackscheck_v0"), 6.5); - // if (config.correlation2Dhist) - // rKzeroShort.fill(HIST("mass_lambda_kshort_after6"), candidate.mK0Short(), candidate.mLambda()); if (tranRad > config.confV0TranRadV0Max) { return false; } hglue.fill(HIST("htrackscheck_v0"), 7.5); - // if (config.correlation2Dhist) - // rKzeroShort.fill(HIST("mass_lambda_kshort_after7"), candidate.mK0Short(), candidate.mLambda()); if (std::fabs(ctauK0s) > config.cMaxV0LifeTime) { return false; } hglue.fill(HIST("htrackscheck_v0"), 8.5); - // if (config.correlation2Dhist) - // rKzeroShort.fill(HIST("mass_lambda_kshort_after8"), candidate.mK0Short(), candidate.mLambda()); - - // if (config.armcut && arm < config.confarmcut) { - // return false; - // } - hglue.fill(HIST("htrackscheck_v0"), 9.5); - // if (config.correlation2Dhist) - // rKzeroShort.fill(HIST("mass_lambda_kshort_after9"), candidate.mK0Short(), candidate.mLambda()); - // if (config.applyCompetingcut && (std::abs(candidate.mLambda() - PDGdatabase->Mass(3122)) <= config.competingcascrejlambda || std::abs(candidate.mAntiLambda() - PDGdatabase->Mass(-3122)) <= config.competingcascrejlambdaanti)) - if (config.applyCompetingcut && (std::abs(candidate.mLambda() - o2::constants::physics::MassLambda0) <= config.competingcascrejlambda || std::abs(candidate.mAntiLambda() - o2::constants::physics::MassLambda0) <= config.competingcascrejlambda)) { + if (config.isapplyCompetingcut && (std::abs(candidate.mLambda() - o2::constants::physics::MassLambda0) <= config.competingcascrejlambda || std::abs(candidate.mAntiLambda() - o2::constants::physics::MassLambda0) <= config.competingcascrejlambda)) { return false; } - hglue.fill(HIST("htrackscheck_v0"), 10.5); + hglue.fill(HIST("htrackscheck_v0"), 9.5); + if (config.correlation2Dhist) rKzeroShort.fill(HIST("mass_lambda_kshort_after10"), candidate.mK0Short(), candidate.mLambda()); if (config.qAv0) { rKzeroShort.fill(HIST("hMassK0ShortSelected"), candidate.mK0Short(), candidate.pt()); - // rKzeroShort.fill(HIST("mass_lambda_kshort_after"), candidate.mK0Short(), candidate.mLambda()); } + if (config.isStandardV0 && candidate.v0Type() != 1) { + return false; // Only standard V0s are selected + } + hglue.fill(HIST("htrackscheck_v0"), 10.5); + if (candidate.mK0Short() < lowmasscutks0 || candidate.mK0Short() > highmasscutks0) { return false; } + hglue.fill(HIST("htrackscheck_v0"), 11.5); + return true; } template - bool isSelectedV0Daughter(T const& track, float charge, double nsigmaV0Daughter, V0s const& /*candidate*/) + bool isSelectedV0Daughter(T const& track, float charge, double nsigmaV0DaughterTPC, V0s const& /*candidate*/) { if (config.qAPID) { // Filling the PID of the V0 daughters in the region of the K0 peak. (charge == 1) ? rKzeroShort.fill(HIST("hNSigmaPosPionK0s_before"), track.tpcInnerParam(), track.tpcNSigmaPi()) : rKzeroShort.fill(HIST("hNSigmaNegPionK0s_before"), track.tpcInnerParam(), track.tpcNSigmaPi()); - // rKzeroShort.fill(HIST("dE_by_dx_TPC"), track.p(), track.tpcSignal()); } const auto eta = track.eta(); const auto tpcNClsF = track.tpcNClsFound(); @@ -524,24 +539,18 @@ struct HigherMassResonances { return false; hglue.fill(HIST("htrackscheck_v0_daughters"), 1.5); - if (!config.globalTracks) { - if (track.tpcNClsCrossedRows() < config.tpcCrossedrows) - return false; - hglue.fill(HIST("htrackscheck_v0_daughters"), 2.5); + if (track.tpcNClsCrossedRows() < config.tpcCrossedrows) + return false; + hglue.fill(HIST("htrackscheck_v0_daughters"), 2.5); - if (track.tpcCrossedRowsOverFindableCls() < config.tpcCrossedrowsOverfcls) - return false; - hglue.fill(HIST("htrackscheck_v0_daughters"), 3.5); + if (track.tpcCrossedRowsOverFindableCls() < config.tpcCrossedrowsOverfcls) + return false; + hglue.fill(HIST("htrackscheck_v0_daughters"), 3.5); - if (tpcNClsF < config.confDaughTPCnclsMin) { - return false; - } - hglue.fill(HIST("htrackscheck_v0_daughters"), 4.5); - } else { - if (!track.isGlobalTrack()) - return false; - hglue.fill(HIST("htrackscheck_v0_daughters"), 4.5); + if (tpcNClsF < config.confDaughTPCnclsMin) { + return false; } + hglue.fill(HIST("htrackscheck_v0_daughters"), 4.5); if (charge < 0 && sign > 0) { return false; @@ -558,11 +567,13 @@ struct HigherMassResonances { } hglue.fill(HIST("htrackscheck_v0_daughters"), 7.5); - if (std::abs(nsigmaV0Daughter) > config.confDaughPIDCuts) { + if (std::abs(nsigmaV0DaughterTPC) > config.confDaughPIDCutTPC) { return false; } hglue.fill(HIST("htrackscheck_v0_daughters"), 8.5); + // if (std::abs()) + if (config.qAPID) { (charge == 1) ? rKzeroShort.fill(HIST("hNSigmaPosPionK0s_after"), track.tpcInnerParam(), track.tpcNSigmaPi()) : rKzeroShort.fill(HIST("hNSigmaNegPionK0s_after"), track.tpcInnerParam(), track.tpcNSigmaPi()); } @@ -584,7 +595,7 @@ struct HigherMassResonances { if (config.qAv0) { rKzeroShort.fill(HIST("angularSeparation"), angle); } - if (config.applyAngSepCut && angle > config.angSepCut) { + if (config.isapplyAngSepCut && angle > config.angSepCut) { return false; } return true; @@ -599,11 +610,11 @@ struct HigherMassResonances { Filter preFilterV0 = (nabs(aod::v0data::dcapostopv) > config.v0settingDcapostopv && nabs(aod::v0data::dcanegtopv) > config.v0settingDcanegtopv); // Defining the type of the daughter tracks - using EventCandidates = soa::Filtered>; + using EventCandidates = soa::Filtered>; using TrackCandidates = soa::Filtered>; using V0TrackCandidate = aod::V0Datas; // For Monte Carlo - using EventCandidatesMC = soa::Join; + using EventCandidatesMC = soa::Join; using TrackCandidatesMC = soa::Filtered>; using V0TrackCandidatesMC = soa::Filtered>; // zBeam direction in lab frame @@ -752,20 +763,26 @@ struct HigherMassResonances { void processSE(EventCandidates::iterator const& collision, TrackCandidates const& /*tracks*/, aod::V0Datas const& V0s) { - hglue.fill(HIST("heventscheck"), 0.5); multiplicity = 0.0; - if (config.cfgMultFOTM) { + + if (config.cSelectMultEstimator == kFT0M) { multiplicity = collision.centFT0M(); - } else { + } else if (config.cSelectMultEstimator == kFT0A) { + multiplicity = collision.centFT0A(); + } else if (config.cSelectMultEstimator == kFT0C) { multiplicity = collision.centFT0C(); - } - if (!eventselection(collision)) { - return; + } else if (config.cSelectMultEstimator == kFV0A) { + multiplicity = collision.centFV0A(); + } else { + multiplicity = collision.centFT0M(); // default } - if (rctCut.requireRCTFlagChecker && !rctCut.rctChecker(collision)) { + if (!selectionEvent(collision, true)) { return; } + // if (rctCut.requireRCTFlagChecker && !rctCut.rctChecker(collision)) { + // return; + // } // auto occupancyNumber = collision.trackOccupancyInTimeRange(); // if (applyOccupancyCut && occupancyNumber < occupancyCut) { @@ -814,26 +831,25 @@ struct HigherMassResonances { continue; } + if (postrack1.hasTOF() && negtrack1.hasTOF() && postrack2.hasTOF() && negtrack2.hasTOF()) { + double nTOFSigmaPos1{postrack1.tofNSigmaPi()}; + double nTOFSigmaNeg1{negtrack1.tofNSigmaPi()}; + double nTOFSigmaPos2{postrack2.tofNSigmaPi()}; + double nTOFSigmaNeg2{negtrack2.tofNSigmaPi()}; + if ((std::abs(nTOFSigmaPos1) > config.confDaughPIDCutTPC) || (std::abs(nTOFSigmaNeg1) > config.confDaughPIDCutTPC) || + (std::abs(nTOFSigmaPos2) > config.confDaughPIDCutTPC) || (std::abs(nTOFSigmaNeg2) > config.confDaughPIDCutTPC)) { + continue; + } + } + if (std::find(v0indexes.begin(), v0indexes.end(), v1.globalIndex()) == v0indexes.end()) { v0indexes.push_back(v1.globalIndex()); } - // if (!(std::find(v0indexes.begin(), v0indexes.end(), v2.globalIndex()) != v0indexes.end())) { - // v0indexes.push_back(v2.globalIndex()); - // } if (v2.globalIndex() <= v1.globalIndex()) { continue; } - // if (config.qAv0Daughters) { - // rKzeroShort.fill(HIST("negative_pt"), negtrack1.pt()); - // rKzeroShort.fill(HIST("positive_pt"), postrack1.pt()); - // rKzeroShort.fill(HIST("negative_eta"), negtrack1.eta()); - // rKzeroShort.fill(HIST("positive_eta"), postrack1.eta()); - // rKzeroShort.fill(HIST("negative_phi"), negtrack1.phi()); - // rKzeroShort.fill(HIST("positive_phi"), postrack1.phi()); - // } - if (postrack1.globalIndex() == postrack2.globalIndex()) { continue; } @@ -845,6 +861,10 @@ struct HigherMassResonances { continue; } + if (config.isApplyEtaCutK0s && (v1.eta() < config.confDaughEta || v2.eta() < config.confDaughEta)) { + continue; + } + if (config.qAv0) { rKzeroShort.fill(HIST("hMasscorrelationbefore"), v1.mK0Short(), v2.mK0Short()); } @@ -855,12 +875,12 @@ struct HigherMassResonances { mother = daughter1 + daughter2; // invariant mass of Kshort pair isMix = false; - if (!config.selectTWOKsOnly) + if (!config.isselectTWOKsOnly) fillInvMass(mother, multiplicity, daughter1, daughter2, isMix); } int sizeofv0indexes = v0indexes.size(); rKzeroShort.fill(HIST("NksProduced"), sizeofv0indexes); - if (config.selectTWOKsOnly && sizeofv0indexes == config.noOfDaughters && allConditionsMet) { + if (config.isselectTWOKsOnly && sizeofv0indexes == config.noOfDaughters && allConditionsMet) { fillInvMass(mother, multiplicity, daughter1, daughter2, false); } v0indexes.clear(); @@ -871,315 +891,41 @@ struct HigherMassResonances { using V0CandidatesDerivedData = soa::Join; using DauTracks = soa::Join; - void processSEderived(EventCandidatesDerivedData::iterator const& collision, TrackCandidates const& /*tracks*/, aod::V0Datas const& V0s) - { - hglue.fill(HIST("heventscheck"), 0.5); - multiplicity = 0.0; - if (config.cfgMultFOTM) { - multiplicity = collision.centFT0M(); - } else { - multiplicity = collision.centFT0C(); - } - if (!eventselection(collision)) { - return; - } + ConfigurableAxis axisVertex{"axisVertex", {20, -10, 10}, "vertex axis for ME mixing"}; + // ConfigurableAxis axisMultiplicityClass{"axisMultiplicityClass", {10, 0, 100}, "multiplicity percentile for ME mixing"}; + ConfigurableAxis axisMultiplicity{"axisMultiplicity", {2000, 0, 10000}, "TPC multiplicity axis for ME mixing"}; - if (rctCut.requireRCTFlagChecker && !rctCut.rctChecker(collision)) { - return; - } + // using BinningTypeTPCMultiplicity = ColumnBinningPolicy; + using BinningTypeFT0M = ColumnBinningPolicy; + using BinningTypeFT0A = ColumnBinningPolicy; + using BinningTypeFT0C = ColumnBinningPolicy; + using BinningTypeFV0A = ColumnBinningPolicy; - // auto occupancyNumber = collision.trackOccupancyInTimeRange(); - // if (applyOccupancyCut && occupancyNumber < occupancyCut) { - // return; - // } + BinningTypeFT0M binningOnFT0M{{axisVertex, axisMultiplicity}, true}; + BinningTypeFT0A binningOnFT0A{{axisVertex, axisMultiplicity}, true}; + BinningTypeFT0C binningOnFT0C{{axisVertex, axisMultiplicity}, true}; + BinningTypeFV0A binningOnFV0A{{axisVertex, axisMultiplicity}, true}; - if (config.qAevents) { - rEventSelection.fill(HIST("hVertexZRec"), collision.posZ()); - rEventSelection.fill(HIST("hmultiplicity"), multiplicity); - // rEventSelection.fill(HIST("multdist_FT0M"), collision.multFT0M()); - // rEventSelection.fill(HIST("multdist_FT0A"), collision.multFT0A()); - // rEventSelection.fill(HIST("multdist_FT0C"), collision.multFT0C()); - // rEventSelection.fill(HIST("hNcontributor"), collision.numContrib()); - } - - std::vector v0indexes; - bool allConditionsMet = 0; - - for (const auto& [v1, v2] : combinations(CombinationsFullIndexPolicy(V0s, V0s))) { - - if (v1.size() == 0 || v2.size() == 0) { - continue; - } - - if (!selectionV0(collision, v1, multiplicity)) { - continue; - } - if (!selectionV0(collision, v2, multiplicity)) { - continue; - } - - auto postrack1 = v1.template posTrack_as(); - auto negtrack1 = v1.template negTrack_as(); - auto postrack2 = v2.template posTrack_as(); - auto negtrack2 = v2.template negTrack_as(); - - double nTPCSigmaPos1{postrack1.tpcNSigmaPi()}; - double nTPCSigmaNeg1{negtrack1.tpcNSigmaPi()}; - double nTPCSigmaPos2{postrack2.tpcNSigmaPi()}; - double nTPCSigmaNeg2{negtrack2.tpcNSigmaPi()}; - - if (!(isSelectedV0Daughter(negtrack1, -1, nTPCSigmaNeg1, v1) && isSelectedV0Daughter(postrack1, 1, nTPCSigmaPos1, v1))) { - continue; - } - if (!(isSelectedV0Daughter(postrack2, 1, nTPCSigmaPos2, v2) && isSelectedV0Daughter(negtrack2, -1, nTPCSigmaNeg2, v2))) { - continue; - } - - if (std::find(v0indexes.begin(), v0indexes.end(), v1.globalIndex()) == v0indexes.end()) { - v0indexes.push_back(v1.globalIndex()); - } - // if (!(std::find(v0indexes.begin(), v0indexes.end(), v2.globalIndex()) != v0indexes.end())) { - // v0indexes.push_back(v2.globalIndex()); - // } - - if (v2.globalIndex() <= v1.globalIndex()) { - continue; - } - - // if (config.qAv0Daughters) { - // rKzeroShort.fill(HIST("negative_pt"), negtrack1.pt()); - // rKzeroShort.fill(HIST("positive_pt"), postrack1.pt()); - // rKzeroShort.fill(HIST("negative_eta"), negtrack1.eta()); - // rKzeroShort.fill(HIST("positive_eta"), postrack1.eta()); - // rKzeroShort.fill(HIST("negative_phi"), negtrack1.phi()); - // rKzeroShort.fill(HIST("positive_phi"), postrack1.phi()); - // } - - if (postrack1.globalIndex() == postrack2.globalIndex()) { - continue; - } - if (negtrack1.globalIndex() == negtrack2.globalIndex()) { - continue; - } - - if (!applyAngSep(v1, v2)) { - continue; - } - - if (config.qAv0) { - rKzeroShort.fill(HIST("hMasscorrelationbefore"), v1.mK0Short(), v2.mK0Short()); - } - allConditionsMet = 1; - daughter1 = ROOT::Math::PxPyPzMVector(v1.px(), v1.py(), v1.pz(), o2::constants::physics::MassK0Short); // Kshort - daughter2 = ROOT::Math::PxPyPzMVector(v2.px(), v2.py(), v2.pz(), o2::constants::physics::MassK0Short); // Kshort - - mother = daughter1 + daughter2; // invariant mass of Kshort pair - isMix = false; - - if (!config.selectTWOKsOnly) - fillInvMass(mother, multiplicity, daughter1, daughter2, isMix); - } - int sizeofv0indexes = v0indexes.size(); - rKzeroShort.fill(HIST("NksProduced"), sizeofv0indexes); - if (config.selectTWOKsOnly && sizeofv0indexes == config.noOfDaughters && allConditionsMet) { - fillInvMass(mother, multiplicity, daughter1, daughter2, false); - } - v0indexes.clear(); - } - PROCESS_SWITCH(HigherMassResonances, processSEderived, "same event process in strangeness derived data", true); - - ConfigurableAxis mevz = {"mevz", {10, -10., 10.}, "mixed event vertex z binning"}; - ConfigurableAxis memult = {"memult", {20, 0, 100}, "mixed event multiplicity binning"}; - - // Processing Event Mixing using BinningType = ColumnBinningPolicy; - BinningType colBinning{{mevz, memult}, true}; + BinningType colBinning{{axisVertex, axisMultiplicity}, true}; // for derived data only Preslice tracksPerCollisionV0Mixed = o2::aod::v0data::straCollisionId; // for derived data only - void processMEderived(EventCandidatesDerivedData const& collisions, TrackCandidates const& /*tracks*/, V0CandidatesDerivedData const& v0s) - { - // auto tracksTuple = std::make_tuple(v0s); - // BinningTypeVertexContributor binningOnPositions1{{mevz, memult}, true}; - // BinningTypeCentralityM binningOnPositions2{{mevz, memult}, true}; - - // SameKindPair pair1{binningOnPositions1, config.cfgNmixedEvents, -1, collisions, tracksTuple, &cache}; // for PbPb - // SameKindPair pair2{binningOnPositions2, config.cfgNmixedEvents, -1, collisions, tracksTuple, &cache}; // for pp - - // if (config.cfgMultFOTM) { - for (const auto& [c1, c2] : selfCombinations(colBinning, config.cfgNmixedEvents, -1, collisions, collisions)) // two different centrality c1 and c2 and tracks corresponding to them - { - - multiplicity = 0.0; - multiplicity = c1.centFT0M(); - - if (!eventselection(c1) || !eventselection(c2)) { - continue; - } - // auto occupancyNumber = c1.trackOccupancyInTimeRange(); - // auto occupancyNumber2 = c2.trackOccupancyInTimeRange(); - // if (applyOccupancyCut && (occupancyNumber < occupancyCut || occupancyNumber2 < occupancyCut)) { - // return; - // } - - if (rctCut.requireRCTFlagChecker && !rctCut.rctChecker(c1)) { - return; - } - if (rctCut.requireRCTFlagChecker && !rctCut.rctChecker(c2)) { - return; - } - auto groupV01 = v0s.sliceBy(tracksPerCollisionV0Mixed, c1.index()); - auto groupV02 = v0s.sliceBy(tracksPerCollisionV0Mixed, c2.index()); - for (const auto& [t1, t2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(groupV01, groupV02))) { - - if (t1.size() == 0 || t2.size() == 0) { - continue; - } - - if (!selectionV0(c1, t1, multiplicity)) - continue; - if (!selectionV0(c2, t2, multiplicity)) - continue; - - auto postrack1 = t1.template posTrackExtra_as(); - auto negtrack1 = t1.template negTrackExtra_as(); - auto postrack2 = t2.template posTrackExtra_as(); - auto negtrack2 = t2.template negTrackExtra_as(); - - if (postrack1.globalIndex() == postrack2.globalIndex()) { - continue; - } - if (negtrack1.globalIndex() == negtrack2.globalIndex()) { - continue; - } - double nTPCSigmaPos1{postrack1.tpcNSigmaPi()}; - double nTPCSigmaNeg1{negtrack1.tpcNSigmaPi()}; - double nTPCSigmaPos2{postrack2.tpcNSigmaPi()}; - double nTPCSigmaNeg2{negtrack2.tpcNSigmaPi()}; - - if (!isSelectedV0Daughter(postrack1, 1, nTPCSigmaPos1, t1)) { - continue; - } - if (!isSelectedV0Daughter(postrack2, 1, nTPCSigmaPos2, t2)) { - continue; - } - if (!isSelectedV0Daughter(negtrack1, -1, nTPCSigmaNeg1, t1)) { - continue; - } - if (!isSelectedV0Daughter(negtrack2, -1, nTPCSigmaNeg2, t2)) { - continue; - } - - daughter1 = ROOT::Math::PxPyPzMVector(t1.px(), t1.py(), t1.pz(), o2::constants::physics::MassK0Short); // Kshort - daughter2 = ROOT::Math::PxPyPzMVector(t2.px(), t2.py(), t2.pz(), o2::constants::physics::MassK0Short); // Kshort - - mother = daughter1 + daughter2; // invariant mass of Kshort pair - isMix = true; - fillInvMass(mother, multiplicity, daughter1, daughter2, isMix); - } - } - // } - // else { - // for (const auto& [c1, tracks1, c2, tracks2] : pair1) // two different centrality c1 and c2 and tracks corresponding to them - // { - // multiplicity = 0.0f; - // multiplicity = c1.centFT0C(); - - // if (!eventselection(c1) || !eventselection(c2)) { - // continue; - // } - // // auto occupancyNumber = c1.trackOccupancyInTimeRange(); - // // auto occupancyNumber2 = c2.trackOccupancyInTimeRange(); - // // if (applyOccupancyCut && (occupancyNumber < occupancyCut || occupancyNumber2 < occupancyCut)) { - // // return; - // // } - - // for (const auto& [t1, t2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { - // if (t1.size() == 0 || t2.size() == 0) { - // continue; - // } - - // if (!selectionV0(c1, t1, multiplicity)) - // continue; - // if (!selectionV0(c2, t2, multiplicity)) - // continue; - - // auto postrack1 = t1.template posTrack_as(); - // auto negtrack1 = t1.template negTrack_as(); - // auto postrack2 = t2.template posTrack_as(); - // auto negtrack2 = t2.template negTrack_as(); - // if (postrack1.globalIndex() == postrack2.globalIndex()) { - // continue; - // } - // if (negtrack1.globalIndex() == negtrack2.globalIndex()) { - // continue; - // } - // double nTPCSigmaPos1{postrack1.tpcNSigmaPi()}; - // double nTPCSigmaNeg1{negtrack1.tpcNSigmaPi()}; - // double nTPCSigmaPos2{postrack2.tpcNSigmaPi()}; - // double nTPCSigmaNeg2{negtrack2.tpcNSigmaPi()}; - - // if (!isSelectedV0Daughter(postrack1, 1, nTPCSigmaPos1, t1)) { - // continue; - // } - // if (!isSelectedV0Daughter(postrack2, 1, nTPCSigmaPos2, t2)) { - // continue; - // } - // if (!isSelectedV0Daughter(negtrack1, -1, nTPCSigmaNeg1, t1)) { - // continue; - // } - // if (!isSelectedV0Daughter(negtrack2, -1, nTPCSigmaNeg2, t2)) { - // continue; - // } - // daughter1 = ROOT::Math::PxPyPzMVector(t1.px(), t1.py(), t1.pz(), o2::constants::physics::MassK0Short); // Kshort - // daughter2 = ROOT::Math::PxPyPzMVector(t2.px(), t2.py(), t2.pz(), o2::constants::physics::MassK0Short); // Kshort - - // mother = daughter1 + daughter2; // invariant mass of Kshort pair - // isMix = true; - // fillInvMass(mother, multiplicity, daughter1, daughter2, isMix); - // } - // } - // } - } - PROCESS_SWITCH(HigherMassResonances, processMEderived, "mixed event process in derived data", true); - - array pvec0; - array pvec1; - // use any one of 3 alias depending on the dataset. If pp then FT0M and if pbpb then FTOC - using BinningTypeTPCMultiplicity = ColumnBinningPolicy; - using BinningTypeCentralityM = ColumnBinningPolicy; - using BinningTypeVertexContributor = ColumnBinningPolicy; - void processME(EventCandidates const& collisions, TrackCandidates const& /*tracks*/, V0TrackCandidate const& v0s) { auto tracksTuple = std::make_tuple(v0s); - BinningTypeVertexContributor binningOnPositions1{{mevz, memult}, true}; - BinningTypeCentralityM binningOnPositions2{{mevz, memult}, true}; + SameKindPair pair1{binningOnFT0M, config.cfgNmixedEvents, -1, collisions, tracksTuple, &cache}; + SameKindPair pair2{binningOnFT0A, config.cfgNmixedEvents, -1, collisions, tracksTuple, &cache}; + SameKindPair pair3{binningOnFT0C, config.cfgNmixedEvents, -1, collisions, tracksTuple, &cache}; + SameKindPair pair4{binningOnFV0A, config.cfgNmixedEvents, -1, collisions, tracksTuple, &cache}; - SameKindPair pair1{binningOnPositions1, config.cfgNmixedEvents, -1, collisions, tracksTuple, &cache}; // for PbPb - SameKindPair pair2{binningOnPositions2, config.cfgNmixedEvents, -1, collisions, tracksTuple, &cache}; // for pp - - if (config.cfgMultFOTM) { - for (const auto& [c1, tracks1, c2, tracks2] : pair2) // two different centrality c1 and c2 and tracks corresponding to them - { + auto runMixing = [&](auto& pair, auto multiplicityGetter) { + for (const auto& [c1, tracks1, c2, tracks2] : pair) { - multiplicity = 0.0; - multiplicity = c1.centFT0M(); + multiplicity = multiplicityGetter(c1); - if (!eventselection(c1) || !eventselection(c2)) { + if (!selectionEvent(c1, false) || !selectionEvent(c2, false)) { continue; } - // auto occupancyNumber = c1.trackOccupancyInTimeRange(); - // auto occupancyNumber2 = c2.trackOccupancyInTimeRange(); - // if (applyOccupancyCut && (occupancyNumber < occupancyCut || occupancyNumber2 < occupancyCut)) { - // return; - // } - - if (rctCut.requireRCTFlagChecker && !rctCut.rctChecker(c1)) { - return; - } - if (rctCut.requireRCTFlagChecker && !rctCut.rctChecker(c2)) { - return; - } for (const auto& [t1, t2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { @@ -1228,66 +974,16 @@ struct HigherMassResonances { fillInvMass(mother, multiplicity, daughter1, daughter2, isMix); } } - } else { - for (const auto& [c1, tracks1, c2, tracks2] : pair1) // two different centrality c1 and c2 and tracks corresponding to them - { - multiplicity = 0.0f; - multiplicity = c1.centFT0C(); - - if (!eventselection(c1) || !eventselection(c2)) { - continue; - } - // auto occupancyNumber = c1.trackOccupancyInTimeRange(); - // auto occupancyNumber2 = c2.trackOccupancyInTimeRange(); - // if (applyOccupancyCut && (occupancyNumber < occupancyCut || occupancyNumber2 < occupancyCut)) { - // return; - // } - - for (const auto& [t1, t2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { - if (t1.size() == 0 || t2.size() == 0) { - continue; - } - - if (!selectionV0(c1, t1, multiplicity)) - continue; - if (!selectionV0(c2, t2, multiplicity)) - continue; - - auto postrack1 = t1.template posTrack_as(); - auto negtrack1 = t1.template negTrack_as(); - auto postrack2 = t2.template posTrack_as(); - auto negtrack2 = t2.template negTrack_as(); - if (postrack1.globalIndex() == postrack2.globalIndex()) { - continue; - } - if (negtrack1.globalIndex() == negtrack2.globalIndex()) { - continue; - } - double nTPCSigmaPos1{postrack1.tpcNSigmaPi()}; - double nTPCSigmaNeg1{negtrack1.tpcNSigmaPi()}; - double nTPCSigmaPos2{postrack2.tpcNSigmaPi()}; - double nTPCSigmaNeg2{negtrack2.tpcNSigmaPi()}; - - if (!isSelectedV0Daughter(postrack1, 1, nTPCSigmaPos1, t1)) { - continue; - } - if (!isSelectedV0Daughter(postrack2, 1, nTPCSigmaPos2, t2)) { - continue; - } - if (!isSelectedV0Daughter(negtrack1, -1, nTPCSigmaNeg1, t1)) { - continue; - } - if (!isSelectedV0Daughter(negtrack2, -1, nTPCSigmaNeg2, t2)) { - continue; - } - daughter1 = ROOT::Math::PxPyPzMVector(t1.px(), t1.py(), t1.pz(), o2::constants::physics::MassK0Short); // Kshort - daughter2 = ROOT::Math::PxPyPzMVector(t2.px(), t2.py(), t2.pz(), o2::constants::physics::MassK0Short); // Kshort - - mother = daughter1 + daughter2; // invariant mass of Kshort pair - isMix = true; - fillInvMass(mother, multiplicity, daughter1, daughter2, isMix); - } - } + }; + // Call mixing based on selected estimator + if (config.cSelectMultEstimator == kFT0M) { + runMixing(pair1, [](const auto& c) { return c.centFT0M(); }); + } else if (config.cSelectMultEstimator == kFT0A) { + runMixing(pair2, [](const auto& c) { return c.centFT0A(); }); + } else if (config.cSelectMultEstimator == kFT0C) { + runMixing(pair3, [](const auto& c) { return c.centFT0C(); }); + } else if (config.cSelectMultEstimator == kFV0A) { + runMixing(pair4, [](const auto& c) { return c.centFV0A(); }); } } PROCESS_SWITCH(HigherMassResonances, processME, "mixed event process", true); @@ -1304,45 +1000,36 @@ struct HigherMassResonances { return; } hMChists.fill(HIST("events_check"), 0.5); - if (std::abs(mcCollision.posZ()) < config.cutzvertex) { - hMChists.fill(HIST("events_check"), 1.5); - } - // int Nchinel = 0; - // for (const auto& mcParticle : mcParticles) { - // auto pdgcode = std::abs(mcParticle.pdgCode()); - // if (mcParticle.isPhysicalPrimary() && (pdgcode == 211 || pdgcode == 321 || pdgcode == 2212 || pdgcode == 11 || pdgcode == 13)) { - // if (std::abs(mcParticle.eta()) < 1.0) { - // Nchinel = Nchinel + 1; - // } - // } - // } - // if (Nchinel > 0 && std::abs(mcCollision.posZ()) < config.cutzvertex) - hMChists.fill(HIST("events_check"), 2.5); std::vector selectedEvents(collisions.size()); int nevts = 0; - multiplicityGen = 0.0; + multiplicityGen = -999.0; for (const auto& collision : collisions) { - if (std::abs(collision.mcCollision().posZ()) > config.cutzvertex) { - continue; - } - if (config.timFrameEvsel && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) { - continue; + // multiplicityGen = collision.centFT0M(); + if (config.cSelectMultEstimator == kFT0M) { + multiplicityGen = collision.centFT0M(); + } else if (config.cSelectMultEstimator == kFT0A) { + multiplicityGen = collision.centFT0A(); + } else if (config.cSelectMultEstimator == kFT0C) { + multiplicityGen = collision.centFT0C(); + } else if (config.cSelectMultEstimator == kFV0A) { + multiplicityGen = collision.centFV0A(); + } else { + multiplicityGen = collision.centFT0M(); // default } - if (config.cTVXEvsel && (!collision.selection_bit(aod::evsel::kIsTriggerTVX))) { + + if (!selectionEvent(collision, true)) { continue; } - multiplicityGen = collision.centFT0M(); - selectedEvents[nevts++] = collision.mcCollision_as().globalIndex(); } selectedEvents.resize(nevts); hMChists.fill(HIST("events_check"), 3.5); const auto evtReconstructedAndSelected = std::find(selectedEvents.begin(), selectedEvents.end(), mcCollision.globalIndex()) != selectedEvents.end(); - if (!config.allGenCollisions && !evtReconstructedAndSelected) { // Check that the event is reconstructed and that the reconstructed events pass the selection + if (!config.isallGenCollisions && !evtReconstructedAndSelected) { // Check that the event is reconstructed and that the reconstructed events pass the selection return; } hMChists.fill(HIST("events_check"), 4.5); @@ -1354,15 +1041,11 @@ struct HigherMassResonances { } hMChists.fill(HIST("events_check"), 5.5); - if (config.applyRapidityMC && std::abs(mcParticle.y()) >= config.rapidityMotherData) { + if (config.isapplyRapidityMC && std::abs(mcParticle.y()) >= config.rapidityMotherData) { continue; } hMChists.fill(HIST("events_check"), 6.5); - // if (counter < 1e3) - // std::cout << "px " << mcParticle.px() << " py " << mcParticle.py() << " pz " << mcParticle.pz() << " y " << mcParticle.y() << std::endl; - // counter++; - auto kDaughters = mcParticle.daughters_as(); if (kDaughters.size() != config.noOfDaughters) { continue; @@ -1405,7 +1088,7 @@ struct HigherMassResonances { hMChists.fill(HIST("GenEta"), mcParticle.eta()); hMChists.fill(HIST("GenPhi"), mcParticle.phi()); - if (config.applyPairRapidityGen && std::abs(lResonanceGen1.Rapidity()) >= config.rapidityMotherData) { + if (config.isapplyPairRapidityGen && std::abs(lResonanceGen1.Rapidity()) >= config.rapidityMotherData) { continue; } @@ -1428,42 +1111,35 @@ struct HigherMassResonances { return; } - auto multiplicity = collision.centFT0M(); - hMChists.fill(HIST("Rec_Multiplicity"), multiplicity); - - hMChists.fill(HIST("events_checkrec"), 0.5); - if (!collision.has_mcCollision()) { - return; + auto multiplicity = -999.0; + if (config.cSelectMultEstimator == kFT0M) { + multiplicity = collision.centFT0M(); + } else if (config.cSelectMultEstimator == kFT0A) { + multiplicity = collision.centFT0A(); + } else if (config.cSelectMultEstimator == kFT0C) { + multiplicity = collision.centFT0C(); + } else if (config.cSelectMultEstimator == kFV0A) { + multiplicity = collision.centFV0A(); + } else { + multiplicity = collision.centFT0M(); // default } - hMChists.fill(HIST("events_checkrec"), 1.5); - // // if (std::abs(collision.mcCollision().posZ()) > config.cutzvertex || !collision.sel8()) { - if (std::abs(collision.mcCollision().posZ()) > config.cutzvertex) { + + if (!selectionEvent(collision, false)) { return; } - hMChists.fill(HIST("events_checkrec"), 2.5); - - // if (config.timFrameEvsel && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) { - // return; - // } - // hMChists.fill(HIST("events_checkrec"), 3.5); - // if (config.cTVXEvsel && (!collision.selection_bit(aod::evsel::kIsTriggerTVX))) { - // return; - // } + hMChists.fill(HIST("Rec_Multiplicity"), multiplicity); - if (!collision.sel8()) { + if (!collision.has_mcCollision()) { return; } - hMChists.fill(HIST("events_checkrec"), 4.5); + hMChists.fill(HIST("MC_mult_after_event_sel"), multiplicity); eventCounter++; - // auto oldindex = -999; for (const auto& v01 : V0s) { for (const auto& v02 : V0s) { - hMChists.fill(HIST("events_checkrec"), 5.5); - if (v02.index() <= v01.index()) { continue; } @@ -1471,7 +1147,7 @@ struct HigherMassResonances { if (!v01.has_mcParticle() || !v02.has_mcParticle()) { continue; } - hMChists.fill(HIST("events_checkrec"), 6.5); + hMChists.fill(HIST("events_checkrec"), 0.5); auto postrack1 = v01.template posTrack_as(); auto negtrack1 = v01.template negTrack_as(); @@ -1481,11 +1157,9 @@ struct HigherMassResonances { if (!postrack1.has_mcParticle() || !postrack2.has_mcParticle()) continue; // Checking that the daughter tracks come from particles and are not fake - hMChists.fill(HIST("events_checkrec"), 7.5); if (!negtrack1.has_mcParticle() || !negtrack2.has_mcParticle()) continue; - hMChists.fill(HIST("events_checkrec"), 8.5); double nTPCSigmaPos1[1]{postrack1.tpcNSigmaPi()}; double nTPCSigmaNeg1[1]{negtrack1.tpcNSigmaPi()}; @@ -1495,17 +1169,16 @@ struct HigherMassResonances { if (!isSelectedV0Daughter(postrack1, 1, nTPCSigmaPos1[0], v01) || !isSelectedV0Daughter(postrack2, 1, nTPCSigmaPos2[0], v02)) { continue; } - hMChists.fill(HIST("events_checkrec"), 9.5); if (!isSelectedV0Daughter(negtrack1, -1, nTPCSigmaNeg1[0], v01) || !isSelectedV0Daughter(negtrack2, -1, nTPCSigmaNeg2[0], v02)) { continue; } - hMChists.fill(HIST("events_checkrec"), 10.5); + hMChists.fill(HIST("events_checkrec"), 1.5); if (!selectionV0(collision, v01, multiplicity) || !selectionV0(collision, v02, multiplicity)) { continue; } - hMChists.fill(HIST("events_checkrec"), 11.5); + hMChists.fill(HIST("events_checkrec"), 2.5); auto mctrackv01 = v01.mcParticle(); auto mctrackv02 = v02.mcParticle(); @@ -1516,7 +1189,7 @@ struct HigherMassResonances { if (std::abs(trackv0PDG1) != PDG_t::kK0Short || std::abs(trackv0PDG2) != PDG_t::kK0Short) { continue; } - hMChists.fill(HIST("events_checkrec"), 12.5); + hMChists.fill(HIST("events_checkrec"), 3.5); for (const auto& mothertrack1 : mctrackv01.mothers_as()) { @@ -1530,17 +1203,15 @@ struct HigherMassResonances { for (const auto& mothertrack2 : mctrackv02.mothers_as()) { - hMChists.fill(HIST("events_checkrec"), 13.5); - if (mothertrack1.pdgCode() != config.pdgCodes[config.selectMCparticles]) { continue; } - hMChists.fill(HIST("events_checkrec"), 14.5); + hMChists.fill(HIST("events_checkrec"), 4.5); if (mothertrack1.pdgCode() != mothertrack2.pdgCode()) { continue; } - hMChists.fill(HIST("events_checkrec"), 15.5); + hMChists.fill(HIST("events_checkrec"), 5.5); gindex2.push_back(mothertrack2.globalIndex()); if (gindex2.size() > 1) { @@ -1548,24 +1219,24 @@ struct HigherMassResonances { continue; } } - hMChists.fill(HIST("events_checkrec"), 16.5); + hMChists.fill(HIST("events_checkrec"), 6.5); if (mothertrack1.globalIndex() != mothertrack2.globalIndex()) { continue; } - hMChists.fill(HIST("events_checkrec"), 17.5); + hMChists.fill(HIST("events_checkrec"), 7.5); if (!mothertrack1.producedByGenerator()) { continue; } - hMChists.fill(HIST("events_checkrec"), 18.5); + hMChists.fill(HIST("events_checkrec"), 8.5); - if (config.applyRapidityMC && std::abs(mothertrack1.y()) >= config.rapidityMotherData) { + if (config.isapplyRapidityMC && std::abs(mothertrack1.y()) >= config.rapidityMotherData) { continue; } - hMChists.fill(HIST("events_checkrec"), 19.5); + hMChists.fill(HIST("events_checkrec"), 9.5); - // if (config.avoidsplitrackMC && oldindex == mothertrack1.globalIndex()) { + // if (config.isavoidsplitrackMC && oldindex == mothertrack1.globalIndex()) { // hMChists.fill(HIST("h1Recsplit"), mothertrack1.pt()); // continue; // } @@ -1592,7 +1263,7 @@ struct HigherMassResonances { hMChists.fill(HIST("RecPhi"), mothertrack1.phi()); hMChists.fill(HIST("RecEta"), mothertrack1.eta()); - if (config.applyPairRapidityRec && std::abs(mother.Rapidity()) >= config.rapidityMotherData) { + if (config.isapplyPairRapidityRec && std::abs(mother.Rapidity()) >= config.rapidityMotherData) { continue; } @@ -1608,6 +1279,181 @@ struct HigherMassResonances { } } PROCESS_SWITCH(HigherMassResonances, processRec, "Process Reconstructed", false); + + void processSEderived(EventCandidatesDerivedData::iterator const& collision, TrackCandidates const& /*tracks*/, aod::V0Datas const& V0s) + { + if (config.cSelectMultEstimator == kFT0M) { + multiplicity = collision.centFT0M(); + } else if (config.cSelectMultEstimator == kFT0A) { + multiplicity = collision.centFT0A(); + } else if (config.cSelectMultEstimator == kFT0C) { + multiplicity = collision.centFT0C(); + } else if (config.cSelectMultEstimator == kFV0A) { + multiplicity = collision.centFV0A(); + } else { + multiplicity = collision.centFT0M(); // default + } + + if (!selectionEvent(collision, true)) { + return; + } + + if (config.qAevents) { + rEventSelection.fill(HIST("hVertexZRec"), collision.posZ()); + rEventSelection.fill(HIST("hmultiplicity"), multiplicity); + } + + std::vector v0indexes; + bool allConditionsMet = 0; + + for (const auto& [v1, v2] : combinations(CombinationsFullIndexPolicy(V0s, V0s))) { + + if (v1.size() == 0 || v2.size() == 0) { + continue; + } + + if (!selectionV0(collision, v1, multiplicity)) { + continue; + } + if (!selectionV0(collision, v2, multiplicity)) { + continue; + } + + auto postrack1 = v1.template posTrack_as(); + auto negtrack1 = v1.template negTrack_as(); + auto postrack2 = v2.template posTrack_as(); + auto negtrack2 = v2.template negTrack_as(); + + double nTPCSigmaPos1{postrack1.tpcNSigmaPi()}; + double nTPCSigmaNeg1{negtrack1.tpcNSigmaPi()}; + double nTPCSigmaPos2{postrack2.tpcNSigmaPi()}; + double nTPCSigmaNeg2{negtrack2.tpcNSigmaPi()}; + + if (!(isSelectedV0Daughter(negtrack1, -1, nTPCSigmaNeg1, v1) && isSelectedV0Daughter(postrack1, 1, nTPCSigmaPos1, v1))) { + continue; + } + if (!(isSelectedV0Daughter(postrack2, 1, nTPCSigmaPos2, v2) && isSelectedV0Daughter(negtrack2, -1, nTPCSigmaNeg2, v2))) { + continue; + } + + if (std::find(v0indexes.begin(), v0indexes.end(), v1.globalIndex()) == v0indexes.end()) { + v0indexes.push_back(v1.globalIndex()); + } + // if (!(std::find(v0indexes.begin(), v0indexes.end(), v2.globalIndex()) != v0indexes.end())) { + // v0indexes.push_back(v2.globalIndex()); + // } + + if (v2.globalIndex() <= v1.globalIndex()) { + continue; + } + + // if (config.qAv0Daughters) { + // rKzeroShort.fill(HIST("negative_pt"), negtrack1.pt()); + // rKzeroShort.fill(HIST("positive_pt"), postrack1.pt()); + // rKzeroShort.fill(HIST("negative_eta"), negtrack1.eta()); + // rKzeroShort.fill(HIST("positive_eta"), postrack1.eta()); + // rKzeroShort.fill(HIST("negative_phi"), negtrack1.phi()); + // rKzeroShort.fill(HIST("positive_phi"), postrack1.phi()); + // } + + if (postrack1.globalIndex() == postrack2.globalIndex()) { + continue; + } + if (negtrack1.globalIndex() == negtrack2.globalIndex()) { + continue; + } + + if (!applyAngSep(v1, v2)) { + continue; + } + + if (config.qAv0) { + rKzeroShort.fill(HIST("hMasscorrelationbefore"), v1.mK0Short(), v2.mK0Short()); + } + allConditionsMet = 1; + daughter1 = ROOT::Math::PxPyPzMVector(v1.px(), v1.py(), v1.pz(), o2::constants::physics::MassK0Short); // Kshort + daughter2 = ROOT::Math::PxPyPzMVector(v2.px(), v2.py(), v2.pz(), o2::constants::physics::MassK0Short); // Kshort + + mother = daughter1 + daughter2; // invariant mass of Kshort pair + isMix = false; + + if (!config.isselectTWOKsOnly) + fillInvMass(mother, multiplicity, daughter1, daughter2, isMix); + } + int sizeofv0indexes = v0indexes.size(); + rKzeroShort.fill(HIST("NksProduced"), sizeofv0indexes); + if (config.isselectTWOKsOnly && sizeofv0indexes == config.noOfDaughters && allConditionsMet) { + fillInvMass(mother, multiplicity, daughter1, daughter2, false); + } + v0indexes.clear(); + } + PROCESS_SWITCH(HigherMassResonances, processSEderived, "same event process in strangeness derived data", false); + + void processMEderived(EventCandidatesDerivedData const& collisions, TrackCandidates const& /*tracks*/, V0CandidatesDerivedData const& v0s) + { + + for (const auto& [c1, c2] : selfCombinations(colBinning, config.cfgNmixedEvents, -1, collisions, collisions)) // two different centrality c1 and c2 and tracks corresponding to them + { + + multiplicity = 0.0; + multiplicity = c1.centFT0M(); + + if (!selectionEvent(c1, false) || !selectionEvent(c2, false)) { + continue; + } + + auto groupV01 = v0s.sliceBy(tracksPerCollisionV0Mixed, c1.index()); + auto groupV02 = v0s.sliceBy(tracksPerCollisionV0Mixed, c2.index()); + for (const auto& [t1, t2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(groupV01, groupV02))) { + + if (t1.size() == 0 || t2.size() == 0) { + continue; + } + + if (!selectionV0(c1, t1, multiplicity)) + continue; + if (!selectionV0(c2, t2, multiplicity)) + continue; + + auto postrack1 = t1.template posTrackExtra_as(); + auto negtrack1 = t1.template negTrackExtra_as(); + auto postrack2 = t2.template posTrackExtra_as(); + auto negtrack2 = t2.template negTrackExtra_as(); + + if (postrack1.globalIndex() == postrack2.globalIndex()) { + continue; + } + if (negtrack1.globalIndex() == negtrack2.globalIndex()) { + continue; + } + double nTPCSigmaPos1{postrack1.tpcNSigmaPi()}; + double nTPCSigmaNeg1{negtrack1.tpcNSigmaPi()}; + double nTPCSigmaPos2{postrack2.tpcNSigmaPi()}; + double nTPCSigmaNeg2{negtrack2.tpcNSigmaPi()}; + + if (!isSelectedV0Daughter(postrack1, 1, nTPCSigmaPos1, t1)) { + continue; + } + if (!isSelectedV0Daughter(postrack2, 1, nTPCSigmaPos2, t2)) { + continue; + } + if (!isSelectedV0Daughter(negtrack1, -1, nTPCSigmaNeg1, t1)) { + continue; + } + if (!isSelectedV0Daughter(negtrack2, -1, nTPCSigmaNeg2, t2)) { + continue; + } + + daughter1 = ROOT::Math::PxPyPzMVector(t1.px(), t1.py(), t1.pz(), o2::constants::physics::MassK0Short); // Kshort + daughter2 = ROOT::Math::PxPyPzMVector(t2.px(), t2.py(), t2.pz(), o2::constants::physics::MassK0Short); // Kshort + + mother = daughter1 + daughter2; // invariant mass of Kshort pair + isMix = true; + fillInvMass(mother, multiplicity, daughter1, daughter2, isMix); + } + } + } + PROCESS_SWITCH(HigherMassResonances, processMEderived, "mixed event process in derived data", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)