Skip to content

Commit e3e2477

Browse files
jikim1290jikim1290
committed
add process with event plane
1 parent 2dd7aca commit e3e2477

1 file changed

Lines changed: 229 additions & 2 deletions

File tree

‎PWGLF/Tasks/Resonances/higherMassResonances.cxx‎

Lines changed: 229 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -20,13 +20,15 @@
2020

2121
#include "Common/CCDB/EventSelectionParams.h"
2222
#include "Common/CCDB/RCTSelectionFlags.h"
23+
#include "Common/Core/EventPlaneHelper.h"
2324
#include "Common/Core/RecoDecay.h"
2425
#include "Common/DataModel/Centrality.h"
2526
#include "Common/DataModel/EventSelection.h" //
2627
#include "Common/DataModel/Multiplicity.h"
2728
#include "Common/DataModel/PIDResponseTOF.h" //
2829
#include "Common/DataModel/PIDResponseTPC.h" //
2930
#include "Common/DataModel/TrackSelectionTables.h"
31+
#include "Common/DataModel/Qvectors.h"
3032

3133
#include <CommonConstants/MathConstants.h>
3234
#include <CommonConstants/PhysicsConstants.h>
@@ -69,7 +71,8 @@ using namespace o2::framework;
6971
using namespace o2::framework::expressions;
7072
using namespace o2::soa;
7173
using namespace o2::aod::rctsel;
72-
// using namespace o2::constants::physics;
74+
using namespace o2::constants::physics;
75+
7376
using std::array;
7477
// FixME: Add INEL>0 selection for the derived data
7578

@@ -181,10 +184,16 @@ struct HigherMassResonances {
181184
// Configurable<float> tpcCrossedrowsOverfcls{"tpcCrossedrowsOverfcls", 0.8, "TPC crossed rows over findable clusters"};
182185
Configurable<int> rotationalCut{"rotationalCut", 10, "Cut value (Rotation angle pi - pi/cut and pi + pi/cut)"};
183186

187+
// event plane configurables
188+
Configurable<int> cfgnTotalSystem{"cfgnTotalSystem", 7, "Total number of detectors in qVectorsTable"};
189+
Configurable<std::string> cfgDetName{"cfgDetName", "FT0C", "The name of detector to be analyzed"};
190+
Configurable<std::string> cfgRefAName{"cfgRefAName", "TPCPos", "The name of detector for reference A"};
191+
Configurable<std::string> cfgRefBName{"cfgRefBName", "TPCNeg", "The name of detector for reference B"};
192+
184193
// // Mass and pT axis as configurables
185194
ConfigurableAxis binsCent{"binsCent", {VARIABLE_WIDTH, 0., 5., 10., 30., 50., 70., 100., 110., 150.}, "Binning of the centrality axis"};
186195
ConfigurableAxis configThnAxisPOL{"configThnAxisPOL", {20, -1.0, 1.0}, "Costheta axis"};
187-
ConfigurableAxis configThnAxisPhi{"configThnAxisPhi", {70, 0.0f, 7.0f}, "Phi axis"}; // 0 to 2pi
196+
ConfigurableAxis configThnAxisPhi{"configThnAxisPhi", {12, 0.0f, o2::constants::math::TwoPI}, "Phi axis"}; // 0 to 2pi
188197
ConfigurableAxis ksMassBins{"ksMassBins", {200, 0.45f, 0.55f}, "K0s invariant mass axis"};
189198
ConfigurableAxis cGlueMassBins{"cGlueMassBins", {200, 0.9f, 3.0f}, "Glueball invariant mass axis"};
190199
ConfigurableAxis cPtBins{"cPtBins", {200, 0.0f, 20.0f}, "Glueball pT axis"};
@@ -218,6 +227,12 @@ struct HigherMassResonances {
218227
// const double massK0s = o2::constants::physics::MassK0Short;
219228
bool isMix = false;
220229

230+
EventPlaneHelper helperEP;
231+
int detId;
232+
int refAId;
233+
int refBId;
234+
float minQvecAmp = 1e-5;
235+
221236
void init(InitContext const&)
222237
{
223238
rctChecker.init(rctCut.cfgEvtRCTFlagCheckerLabel, rctCut.cfgEvtRCTFlagCheckerZDCCheck, rctCut.cfgEvtRCTFlagCheckerLimitAcceptAsBad);
@@ -233,6 +248,9 @@ struct HigherMassResonances {
233248
AxisSpec deltaMAxis = {config.configAxisDeltaM, "#Delta M (GeV/#it{c}^{2})"};
234249
AxisSpec angleSepAxis = {config.configAxisAngleSep, "Angular separation between V0s"};
235250
AxisSpec ptCorrAxis = {config.configAxisPtCorr, "Pt correlation between two K0s"};
251+
AxisSpec axisCentQA = {100, 0, 100, ""};
252+
AxisSpec axisEvtPlQA = {100, -o2::constants::math::PI, o2::constants::math::PI, ""};
253+
AxisSpec axisEvtResPlQA = {102, -1.02, 1.02, ""};
236254

237255
// THnSparses
238256
std::array<int, 5> sparses = {config.activateHelicityFrame, config.activateCollinsSoperFrame, config.activateProductionFrame, config.activateBeamAxisFrame, config.activateRandomFrame};
@@ -341,6 +359,10 @@ struct HigherMassResonances {
341359
hglue.add("h3glueInvMassDS", "h3glueInvMassDS", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, deltaMAxis, angleSepAxis, ptCorrAxis}, true);
342360
hglue.add("h3glueInvMassME", "h3glueInvMassME", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, deltaMAxis, angleSepAxis, ptCorrAxis}, true);
343361
hglue.add("h3glueInvMassRot", "h3glueInvMassRot", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, deltaMAxis, angleSepAxis, ptCorrAxis}, true);
362+
if (doprocessSEEP) {
363+
hglue.add("h3glueInvMassEPDS", "h3glueInvMassEPDS", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPhi}, true);
364+
hglue.add("h3glueInvMassEPRot", "h3glueInvMassEPRot", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPhi}, true);
365+
}
344366
}
345367

346368
// K0s topological/PID cuts
@@ -417,6 +439,42 @@ struct HigherMassResonances {
417439
hMChists.add("MCcorrections/hSignalLossNumerator3", "Kstar generated after event selection", kTH2F, {{ptAxis}, {multiplicityAxis}});
418440
hMChists.add("MCcorrections/hMultvsCent", "Kstar generated after event selection", kTH2F, {{multiplicityAxis}, {multiplicityAxis}});
419441
}
442+
443+
if (doprocessSEEP) {
444+
detId = getdetId(config.cfgDetName);
445+
refAId = getdetId(config.cfgRefAName);
446+
refBId = getdetId(config.cfgRefBName);
447+
448+
hglue.add("EpDet", "", {HistType::kTH2F, {axisCentQA, axisEvtPlQA}});
449+
hglue.add("EpRefA", "", {HistType::kTH2F, {axisCentQA, axisEvtPlQA}});
450+
hglue.add("EpRefB", "", {HistType::kTH2F, {axisCentQA, axisEvtPlQA}});
451+
452+
hglue.add("EpResDetRefA", "", {HistType::kTH2F, {axisCentQA, axisEvtResPlQA}});
453+
hglue.add("EpResDetRefB", "", {HistType::kTH2F, {axisCentQA, axisEvtResPlQA}});
454+
hglue.add("EpResRefARefB", "", {HistType::kTH2F, {axisCentQA, axisEvtResPlQA}});
455+
}
456+
}
457+
458+
template <typename T>
459+
int getdetId(const T& name)
460+
{
461+
if (name.value == "FT0C") {
462+
return 0;
463+
} else if (name.value == "FT0A") {
464+
return 1;
465+
} else if (name.value == "FT0M") {
466+
return 2;
467+
} else if (name.value == "FV0A") {
468+
return 3;
469+
} else if (name.value == "TPCPos") {
470+
return 4;
471+
} else if (name.value == "TPCNeg") {
472+
return 5;
473+
} else if (name.value == "TPCTot") {
474+
return 6;
475+
} else {
476+
return 0;
477+
}
420478
}
421479

422480
template <typename Coll>
@@ -799,6 +857,7 @@ struct HigherMassResonances {
799857
using EventCandidates = soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::FT0Mults, aod::FV0Mults, aod::MultZeqs, aod::CentFT0Ms, aod::CentFT0As, aod::CentFT0Cs, aod::CentFV0As, aod::Mults, aod::PVMults>>;
800858
using TrackCandidates = soa::Filtered<soa::Join<aod::Tracks, aod::TrackSelection, aod::TracksExtra, aod::TracksDCA, aod::pidTPCFullPi, aod::pidTOFFullPi>>;
801859
using V0TrackCandidate = soa::Join<aod::V0Datas, aod::V0TOFPIDs, aod::V0TOFNSigmas>;
860+
using EventCandidateEPs = soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::FT0Mults, aod::FV0Mults, aod::MultZeqs, aod::CentFT0Ms, aod::CentFT0As, aod::CentFT0Cs, aod::CentFV0As, aod::Mults, aod::PVMults, aod::Qvectors>>;
802861
// For Monte Carlo
803862
using EventCandidatesMC = soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels, aod::CentFT0Ms, aod::CentFT0As, aod::CentFT0Cs, aod::MultZeqs, aod::FT0Mults, aod::PVMults, aod::CentFV0As>;
804863
using TrackCandidatesMC = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtra, aod::TracksDCA, aod::TrackSelection, aod::pidTPCFullKa, aod::pidTOFFullKa, aod::pidTPCFullPi, aod::pidTOFFullPi, aod::McTrackLabels>>;
@@ -1157,6 +1216,174 @@ struct HigherMassResonances {
11571216
}
11581217
PROCESS_SWITCH(HigherMassResonances, processSE, "same event process", true);
11591218

1219+
void processSEEP(EventCandidateEPs::iterator const& collision, TrackCandidates const& /*tracks*/, V0TrackCandidate const& V0s)
1220+
{
1221+
multiplicity = 0.0;
1222+
1223+
if (config.cSelectMultEstimator == kFT0M) {
1224+
multiplicity = collision.centFT0M();
1225+
} else if (config.cSelectMultEstimator == kFT0A) {
1226+
multiplicity = collision.centFT0A();
1227+
} else if (config.cSelectMultEstimator == kFT0C) {
1228+
multiplicity = collision.centFT0C();
1229+
} else if (config.cSelectMultEstimator == kFV0A) {
1230+
multiplicity = collision.centFV0A();
1231+
} else {
1232+
multiplicity = collision.centFT0C(); // default
1233+
}
1234+
1235+
if (!selectionEvent(collision, true)) {
1236+
return;
1237+
}
1238+
1239+
if (collision.qvecAmp()[detId] < minQvecAmp || collision.qvecAmp()[refAId] < minQvecAmp || collision.qvecAmp()[refBId] < minQvecAmp) {
1240+
return;
1241+
}
1242+
1243+
if (config.qAevents) {
1244+
rEventSelection.fill(HIST("hVertexZRec"), collision.posZ());
1245+
rEventSelection.fill(HIST("hmultiplicity"), multiplicity);
1246+
}
1247+
1248+
float eps[3] = {0.};
1249+
eps[0] = helperEP.GetEventPlane(collision.qvecRe()[4 * detId + 3], collision.qvecIm()[4 * detId + 3], 2);
1250+
eps[1] = helperEP.GetEventPlane(collision.qvecRe()[4 * refAId + 3], collision.qvecIm()[4 * refAId + 3], 2);
1251+
eps[2] = helperEP.GetEventPlane(collision.qvecRe()[4 * refBId + 3], collision.qvecIm()[4 * refBId + 3], 2);
1252+
1253+
float resNumA = helperEP.GetResolution(eps[0], eps[1], 2);
1254+
float resNumB = helperEP.GetResolution(eps[0], eps[2], 2);
1255+
float resDenom = helperEP.GetResolution(eps[1], eps[2], 2);
1256+
1257+
hglue.fill(HIST("EpDet"), multiplicity, eps[0]);
1258+
hglue.fill(HIST("EpRefA"), multiplicity, eps[1]);
1259+
hglue.fill(HIST("EpRefB"), multiplicity, eps[2]);
1260+
1261+
hglue.fill(HIST("EpResDetRefA"), multiplicity, resNumA);
1262+
hglue.fill(HIST("EpResDetRefB"), multiplicity, resNumB);
1263+
hglue.fill(HIST("EpResRefARefB"), multiplicity, resDenom);
1264+
1265+
std::vector<int> v0indexes;
1266+
bool allConditionsMet = 0;
1267+
1268+
for (const auto& [v1, v2] : combinations(CombinationsFullIndexPolicy(V0s, V0s))) {
1269+
1270+
if (v1.size() == 0 || v2.size() == 0) {
1271+
continue;
1272+
}
1273+
1274+
if (!selectionV0(collision, v1, multiplicity)) {
1275+
continue;
1276+
}
1277+
if (!selectionV0(collision, v2, multiplicity)) {
1278+
continue;
1279+
}
1280+
1281+
auto postrack1 = v1.template posTrack_as<TrackCandidates>();
1282+
auto negtrack1 = v1.template negTrack_as<TrackCandidates>();
1283+
auto postrack2 = v2.template posTrack_as<TrackCandidates>();
1284+
auto negtrack2 = v2.template negTrack_as<TrackCandidates>();
1285+
1286+
double nTPCSigmaPos1{postrack1.tpcNSigmaPi()};
1287+
double nTPCSigmaNeg1{negtrack1.tpcNSigmaPi()};
1288+
double nTPCSigmaPos2{postrack2.tpcNSigmaPi()};
1289+
double nTPCSigmaNeg2{negtrack2.tpcNSigmaPi()};
1290+
1291+
if (!(isSelectedV0Daughter(negtrack1, -1, nTPCSigmaNeg1, v1) && isSelectedV0Daughter(postrack1, 1, nTPCSigmaPos1, v1))) {
1292+
continue;
1293+
}
1294+
if (!(isSelectedV0Daughter(postrack2, 1, nTPCSigmaPos2, v2) && isSelectedV0Daughter(negtrack2, -1, nTPCSigmaNeg2, v2))) {
1295+
continue;
1296+
}
1297+
1298+
if (std::find(v0indexes.begin(), v0indexes.end(), v1.globalIndex()) == v0indexes.end()) {
1299+
v0indexes.push_back(v1.globalIndex());
1300+
}
1301+
1302+
if (v2.globalIndex() <= v1.globalIndex()) {
1303+
continue;
1304+
}
1305+
1306+
if (postrack1.globalIndex() == postrack2.globalIndex()) {
1307+
continue;
1308+
}
1309+
if (negtrack1.globalIndex() == negtrack2.globalIndex()) {
1310+
continue;
1311+
}
1312+
1313+
double deltaRDaugherPos = std::sqrt(TVector2::Phi_mpi_pi(postrack1.phi() - negtrack1.phi()) * TVector2::Phi_mpi_pi(postrack1.phi() - negtrack1.phi()) + (postrack1.eta() - negtrack1.eta()) * (postrack1.eta() - negtrack1.eta()));
1314+
double deltaRDaugherNeg = std::sqrt(TVector2::Phi_mpi_pi(postrack2.phi() - negtrack2.phi()) * TVector2::Phi_mpi_pi(postrack2.phi() - negtrack2.phi()) + (postrack2.eta() - negtrack2.eta()) * (postrack2.eta() - negtrack2.eta()));
1315+
1316+
if (config.qAv0) {
1317+
rKzeroShort.fill(HIST("hDauDeltaR"), deltaRDaugherPos, deltaRDaugherNeg);
1318+
}
1319+
1320+
if (deltaRDaugherPos < config.deltaRDaugherCut || deltaRDaugherNeg < config.deltaRDaugherCut) {
1321+
continue;
1322+
}
1323+
1324+
if (config.isApplyEtaCutK0s && (v1.eta() < config.confDaughEta || v2.eta() < config.confDaughEta)) {
1325+
continue;
1326+
}
1327+
1328+
allConditionsMet = 1;
1329+
daughter1 = ROOT::Math::PxPyPzMVector(v1.px(), v1.py(), v1.pz(), o2::constants::physics::MassK0Short); // Kshort
1330+
daughter2 = ROOT::Math::PxPyPzMVector(v2.px(), v2.py(), v2.pz(), o2::constants::physics::MassK0Short); // Kshort
1331+
1332+
mother = daughter1 + daughter2; // invariant mass of Kshort pair
1333+
isMix = false;
1334+
1335+
const double deltaMass = deltaM(v1.mK0Short(), v2.mK0Short());
1336+
if (config.qAv0) {
1337+
rKzeroShort.fill(HIST("hK0ShortMassCorr"), v1.mK0Short(), v2.mK0Short(), deltaMass);
1338+
}
1339+
1340+
if (!config.qAOptimisation) {
1341+
if (deltaMass > config.cMaxDeltaM) {
1342+
continue;
1343+
}
1344+
}
1345+
1346+
const double ptCorr = (mother.Pt() - daughter1.Pt() != 0.) ? daughter1.Pt() / (mother.Pt() - daughter1.Pt()) : 0.;
1347+
if (config.qAv0) {
1348+
rKzeroShort.fill(HIST("hK0sPtCorrelation"), ptCorr);
1349+
}
1350+
1351+
double deltaRvalue = std::sqrt(TVector2::Phi_mpi_pi(v1.phi() - v2.phi()) * TVector2::Phi_mpi_pi(v1.phi() - v2.phi()) + (v1.eta() - v2.eta()) * (v1.eta() - v2.eta()));
1352+
1353+
if (!config.qAOptimisation) {
1354+
if (deltaRvalue < config.deltaRK0sCut) {
1355+
continue;
1356+
}
1357+
}
1358+
1359+
if (!config.isselectTWOKsOnly && config.qAOptimisation) {
1360+
1361+
if (std::abs(mother.Rapidity()) < config.rapidityMotherData) {
1362+
hglue.fill(HIST("h3glueInvMassEPDS"), multiplicity, mother.Pt(), mother.M(), RecoDecay::constrainAngle(2.0 * mother.Phi() - 2.0 * eps[0]));
1363+
}
1364+
1365+
for (int i = 0; i < config.cRotations; i++) {
1366+
double theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut);
1367+
1368+
daughterRot = ROOT::Math::PxPyPzMVector(daughter1.Px() * std::cos(theta2) - daughter1.Py() * std::sin(theta2), daughter1.Px() * std::sin(theta2) + daughter1.Py() * std::cos(theta2), daughter1.Pz(), daughter1.M());
1369+
1370+
motherRot = daughterRot + daughter2;
1371+
1372+
if (motherRot.Rapidity() < config.rapidityMotherData) {
1373+
hglue.fill(HIST("h3glueInvMassEPRot"), multiplicity, motherRot.Pt(), motherRot.M(), RecoDecay::constrainAngle(2.0 * motherRot.Phi() - 2.0 * eps[0]));
1374+
}
1375+
}
1376+
}
1377+
}
1378+
int sizeofv0indexes = v0indexes.size();
1379+
rKzeroShort.fill(HIST("NksProduced"), sizeofv0indexes);
1380+
if (config.isselectTWOKsOnly && sizeofv0indexes == config.noOfDaughters && allConditionsMet) {
1381+
fillInvMass(mother, multiplicity, daughter1, daughter2, false);
1382+
}
1383+
v0indexes.clear();
1384+
}
1385+
PROCESS_SWITCH(HigherMassResonances, processSEEP, "same event process with EP table", true);
1386+
11601387
ConfigurableAxis axisVertex{"axisVertex", {20, -10, 10}, "vertex axis for ME mixing"};
11611388
// ConfigurableAxis axisMultiplicityClass{"axisMultiplicityClass", {10, 0, 100}, "multiplicity percentile for ME mixing"};
11621389
ConfigurableAxis axisMultiplicity{"axisMultiplicity", {2000, 0, 10000}, "TPC multiplicity axis for ME mixing"};

0 commit comments

Comments
 (0)