Skip to content

Commit 4ac4669

Browse files
yhambardyhambard
authored
[PWGEM] phosElId dispersion calculation fix (AliceO2Group#11372)
1 parent 6ab8317 commit 4ac4669

File tree

1 file changed

+75
-49
lines changed

1 file changed

+75
-49
lines changed

‎PWGEM/Tasks/phosElId.cxx‎

Lines changed: 75 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -81,6 +81,30 @@ enum CentEstimators { FV0A,
8181
FDDM,
8282
NTPV };
8383

84+
bool testLambda(float pt, float l1, float l2, float cutThreshold, bool useNegativeCrossTerm)
85+
{
86+
float l2Mean = 1.53126f + 9.50835e+06f / (1.f + 1.08728e+07f * pt + 1.73420e+06f * pt * pt);
87+
float l1Mean = 1.12365f + 0.123770f * std::exp(-pt * 0.246551f) + 5.30000e-03f * pt;
88+
float l2Sigma = 6.48260e-02f + 7.60261e+10f / (1.f + 1.53012e+11f * pt + 5.01265e+05f * pt * pt) + 9.00000e-03f * pt;
89+
float l1Sigma = 4.44719e-04f + 6.99839e-01f / (1.f + 1.22497e+00f * pt + 6.78604e-07f * pt * pt) + 9.00000e-03f * pt;
90+
float c = -0.35f - 0.550f * std::exp(-0.390730f * pt);
91+
if (l1Sigma == 0.f || l2Sigma == 0.f)
92+
return false;
93+
94+
float term1 = 0.5f * (l1 - l1Mean) * (l1 - l1Mean) / (l1Sigma * l1Sigma);
95+
float term2 = 0.5f * (l2 - l2Mean) * (l2 - l2Mean) / (l2Sigma * l2Sigma);
96+
float crossTerm = 0.5f * c * (l1 - l1Mean) * (l2 - l2Mean) / (l1Sigma * l2Sigma);
97+
98+
float rSquared;
99+
if (useNegativeCrossTerm) {
100+
rSquared = term1 + term2 - crossTerm;
101+
} else {
102+
rSquared = term1 + term2 + crossTerm;
103+
}
104+
105+
return rSquared < cutThreshold;
106+
}
107+
84108
struct PhosElId {
85109

86110
Produces<o2::aod::PHOSMatchindexTable> phosMatch;
@@ -89,7 +113,10 @@ struct PhosElId {
89113
using MyTracks = soa::Join<aod::Tracks, aod::TracksExtra, aod::TracksDCA,
90114
aod::TracksDCACov, aod::pidTOFFullEl, aod::pidTPCFullEl,
91115
aod::pidTPCFullPi, aod::pidTPCFullKa, aod::pidTPCFullPr>;
92-
Configurable<bool> isSel8{"isSel8", 1, "check if event is Single Event Latch-up 8"};
116+
Configurable<bool> isSel8{"isSel8", 1, "check if event is Single Event Latch-up 8"},
117+
mSwapM20M02ForTestLambda{"mSwapM20M02ForTestLambda", false, "Swap m20 and m02 arguments for testLambda (false for note's correct order, true for swapped/original incorrect order)"},
118+
mUseNegativeCrossTerm{"mUseNegativeCrossTerm", true, "Use negative sign for the cross-term in testLambda (true for analysis note version, false for old version)"};
119+
93120
Configurable<float> mColMaxZ{"mColMaxZ", 10.f, "maximum z accepted in analysis"},
94121
mMinCluE{"mMinCluE", 0.3, "Minimum cluster energy for analysis"},
95122
mMinCluTime{"minCluTime", -25.e-9, "Min. cluster time"},
@@ -125,7 +152,8 @@ struct PhosElId {
125152
TPCNSigmaKaMax{"TPCNSigmaKaMax", {4.f}, "max TPC nsigma kaon for exclusion"},
126153
TOFNSigmaElMin{"TOFNSigmaElMin", {-3.f}, "min TOF nsigma e for inclusion"},
127154
TOFNSigmaElMax{"TOFNSigmaElMax", {3.f}, "max TOF nsigma e for inclusion"},
128-
NsigmaTrackMatch{"NsigmaTrackMatch", {2.f}, "PHOS Track Matching Nsigma for inclusion"};
155+
NsigmaTrackMatch{"NsigmaTrackMatch", {2.f}, "PHOS Track Matching Nsigma for inclusion"},
156+
mShowerShapeCutValue{"mShowerShapeCutValue", 4.f, "Cut threshold for testLambda shower shape"};
129157

130158
Configurable<int> mEvSelTrig{"mEvSelTrig", kTVXinPHOS, "Select events with this trigger"},
131159
mAmountOfModules{"mAmountOfModules", 4, "amount of modules for PHOS"},
@@ -371,7 +399,11 @@ struct PhosElId {
371399
clu.time() > mMaxCluTime || clu.time() < mMinCluTime)
372400
continue;
373401

374-
bool isDispOK = testLambda(cluE, clu.m02(), clu.m20());
402+
bool isDispOK = false;
403+
if (mSwapM20M02ForTestLambda)
404+
isDispOK = testLambda(cluE, clu.m02(), clu.m20(), mShowerShapeCutValue, mUseNegativeCrossTerm);
405+
else
406+
isDispOK = testLambda(cluE, clu.m20(), clu.m02(), mShowerShapeCutValue, mUseNegativeCrossTerm);
375407
float posX = clu.x(), posZ = clu.z(), dX = trackX - posX, dZ = trackZ - posZ, Ep = cluE / trackMom;
376408

377409
mHistManager.fill(HIST("coordinateMatching/hdZpmod"), dZ, trackPT, module);
@@ -446,6 +478,11 @@ struct PhosElId {
446478
for (auto const& clu : clusters) {
447479
double cluE = clu.e(), cluTime = clu.time();
448480
int mod = clu.mod();
481+
bool isDispOK = false;
482+
if (mSwapM20M02ForTestLambda)
483+
isDispOK = testLambda(cluE, clu.m02(), clu.m20(), mShowerShapeCutValue, mUseNegativeCrossTerm);
484+
else
485+
isDispOK = testLambda(cluE, clu.m20(), clu.m02(), mShowerShapeCutValue, mUseNegativeCrossTerm);
449486
if (cluE > mMinCluE) {
450487
mHistManager.fill(HIST("clusterSpectra/hCluE_mod_energy_cut"), cluE, mod);
451488
mHistManager.fill(HIST("clusterSpectra/hCluE_v_mod_v_time"), cluE, cluTime * 1e9, mod);
@@ -455,7 +492,7 @@ struct PhosElId {
455492
mHistManager.fill(HIST("clusterSpectra/hCluE_mod_cell_cut"), cluE, mod);
456493
mHistManager.fill(HIST("coordinateMatching/hCluXZ_mod"), clu.x(), clu.z(), mod);
457494
mHistManager.fill(HIST("clusterSpectra/hCluE_ncells_mod"), cluE, clu.ncell(), mod);
458-
if (testLambda(cluE, clu.m02(), clu.m20()))
495+
if (isDispOK)
459496
mHistManager.fill(HIST("clusterSpectra/hCluE_mod_disp"), cluE, mod);
460497
}
461498
}
@@ -482,7 +519,7 @@ struct PhosElId {
482519
mHistManager.fill(HIST("singleLoop/trackdist/clusterSpectra/hCluE_v_pt_Nsigma"), cluE, matchedTrack.pt(), mod);
483520
mHistManager.fill(HIST("singleLoop/trackdist/energyMomentumRatio/hEp_v_pt_Nsigma"), cluE / matchedTrack.p(), matchedTrack.pt(), mod);
484521
mHistManager.fill(HIST("singleLoop/trackdist/energyMomentumRatio/hEp_v_E_Nsigma"), cluE / matchedTrack.p(), cluE, mod);
485-
if (testLambda(cluE, clu.m02(), clu.m20())) {
522+
if (isDispOK) {
486523
mHistManager.fill(HIST("singleLoop/trackdist/clusterSpectra/hCluE_v_pt_Nsigma_disp"), cluE, matchedTrack.pt(), mod);
487524
mHistManager.fill(HIST("singleLoop/trackdist/energyMomentumRatio/hEp_v_pt_Nsigma_disp"), cluE / matchedTrack.p(), matchedTrack.pt(), mod);
488525
mHistManager.fill(HIST("singleLoop/trackdist/energyMomentumRatio/hEp_v_E_Nsigma_disp"), cluE / matchedTrack.p(), cluE, mod);
@@ -505,13 +542,13 @@ struct PhosElId {
505542
isElectron = true;
506543
}
507544
if (isElectron) {
508-
mHistManager.fill(HIST("singleLoop/trackdist/clusterSpectra/hCluE_v_pt_Nsigma_TPC"), cluE, matchedTrack.pt(), mod);
509-
mHistManager.fill(HIST("singleLoop/trackdist/energyMomentumRatio/hEp_v_pt_Nsigma_TPC"), cluE / matchedTrack.p(), matchedTrack.pt(), mod);
510-
mHistManager.fill(HIST("singleLoop/trackdist/energyMomentumRatio/hEp_v_E_Nsigma_TPC"), cluE / matchedTrack.p(), cluE, mod);
511-
if (testLambda(cluE, clu.m02(), clu.m20())) {
512-
mHistManager.fill(HIST("singleLoop/trackdist/clusterSpectra/hCluE_v_pt_Nsigma_disp_TPC"), cluE, matchedTrack.pt(), mod);
513-
mHistManager.fill(HIST("singleLoop/trackdist/energyMomentumRatio/hEp_v_pt_Nsigma_disp_TPC"), cluE / matchedTrack.p(), matchedTrack.pt(), mod);
514-
mHistManager.fill(HIST("singleLoop/trackdist/energyMomentumRatio/hEp_v_E_Nsigma_disp_TPC"), cluE / matchedTrack.p(), cluE, mod);
545+
mHistManager.fill(HIST("singleLoop/trackdist/clusterSpectra/hCluE_v_pt_Nsigma_TPCel"), cluE, matchedTrack.pt(), mod);
546+
mHistManager.fill(HIST("singleLoop/trackdist/energyMomentumRatio/hEp_v_pt_Nsigma_TPCel"), cluE / matchedTrack.p(), matchedTrack.pt(), mod);
547+
mHistManager.fill(HIST("singleLoop/trackdist/energyMomentumRatio/hEp_v_E_Nsigma_TPCel"), cluE / matchedTrack.p(), cluE, mod);
548+
if (isDispOK) {
549+
mHistManager.fill(HIST("singleLoop/trackdist/clusterSpectra/hCluE_v_pt_Nsigma_disp_TPCel"), cluE, matchedTrack.pt(), mod);
550+
mHistManager.fill(HIST("singleLoop/trackdist/energyMomentumRatio/hEp_v_pt_Nsigma_disp_TPCel"), cluE / matchedTrack.p(), matchedTrack.pt(), mod);
551+
mHistManager.fill(HIST("singleLoop/trackdist/energyMomentumRatio/hEp_v_E_Nsigma_disp_TPCel"), cluE / matchedTrack.p(), cluE, mod);
515552
}
516553
}
517554
} // end of cluster loop
@@ -593,21 +630,6 @@ struct PhosElId {
593630
trackZ = posL[1];
594631
return true;
595632
}
596-
//_____________________________________________________________________________
597-
bool testLambda(float pt, float l1, float l2)
598-
{
599-
// Parameterization for full dispersion
600-
float l2Mean = 1.53126 + 9.50835e+06 / (1. + 1.08728e+07 * pt + 1.73420e+06 * pt * pt);
601-
float l1Mean = 1.12365 + 0.123770 * std::exp(-pt * 0.246551) + 5.30000e-03 * pt;
602-
float l2Sigma = 6.48260e-02 + 7.60261e+10 / (1. + 1.53012e+11 * pt + 5.01265e+05 * pt * pt) + 9.00000e-03 * pt;
603-
float l1Sigma = 4.44719e-04 + 6.99839e-01 / (1. + 1.22497e+00 * pt + 6.78604e-07 * pt * pt) + 9.00000e-03 * pt;
604-
float c = -0.35 - 0.550 * std::exp(-0.390730 * pt);
605-
606-
return 0.5 * (l1 - l1Mean) * (l1 - l1Mean) / l1Sigma / l1Sigma +
607-
0.5 * (l2 - l2Mean) * (l2 - l2Mean) / l2Sigma / l2Sigma +
608-
0.5 * c * (l1 - l1Mean) * (l2 - l2Mean) / l1Sigma / l2Sigma <
609-
4.;
610-
}
611633
};
612634

613635
struct MassSpectra {
@@ -863,7 +885,10 @@ struct TpcElIdMassSpectrum {
863885
using MyTracks = soa::Join<aod::Tracks, aod::TracksExtra, aod::TracksDCA,
864886
aod::TracksDCACov, aod::pidTOFFullEl, aod::pidTPCFullEl,
865887
aod::pidTPCFullPi, aod::pidTPCFullKa, aod::pidTPCFullPr>;
866-
Configurable<bool> isSel8{"isSel8", 1, "check if event is Single Event Latch-up 8"};
888+
Configurable<bool> isSel8{"isSel8", 1, "check if event is Single Event Latch-up 8"},
889+
mSwapM20M02ForTestLambda{"mSwapM20M02ForTestLambda", false, "Swap m20 and m02 arguments for testLambda (false for note's correct order, true for swapped/original incorrect order)"},
890+
mUseNegativeCrossTerm{"mUseNegativeCrossTerm", true, "Use negative sign for the cross-term in testLambda (true for analysis note version, false for old version)"};
891+
867892
Configurable<float> mColMaxZ{"mColMaxZ", 10.f, "maximum z accepted in analysis"},
868893
mMinCluE{"mMinCluE", 0.1, "Minimum cluster energy for photons in the analysis"},
869894
mCutMIPCluE{"mCutMIPCluE", 0.3, "Min cluster energy to reject MIPs in the analysis"},
@@ -903,7 +928,8 @@ struct TpcElIdMassSpectrum {
903928
eeMassMin{"eeMassMin", {2.9f}, "J/psi(e+e-) Mass corridor lower limit (for Chic selection)"},
904929
eeMassMax{"eeMassMax", {3.3f}, "J/psi(e+e-) Mass corridor upper limit (for Chic selection)"},
905930
JpsiMass{"JpsiMass", {3.097f}, "J/psi Mass constant"},
906-
mMassSpectrumLowerCutoff{"mMassSpectrumLowerCutoff", {0.01f}, "Used to exclude 0+0 masses"};
931+
mMassSpectrumLowerCutoff{"mMassSpectrumLowerCutoff", {0.01f}, "Used to exclude 0+0 masses"},
932+
mShowerShapeCutValue{"mShowerShapeCutValue", 4.f, "Cut threshold for testLambda shower shape"};
907933

908934
Configurable<int> mEvSelTrig{"mEvSelTrig", kTVXinPHOS, "Select events with this trigger"},
909935
CentBinning{"CentBinning", 10, "Binning for centrality"},
@@ -1223,7 +1249,11 @@ struct TpcElIdMassSpectrum {
12231249
continue;
12241250
bool matchFlag = false;
12251251
bool isJpsi = (pairMass > eeMassMin && pairMass < eeMassMax);
1226-
bool isDispOK = testLambda(cluE, gamma.m02(), gamma.m20());
1252+
bool isDispOK = false;
1253+
if (mSwapM20M02ForTestLambda)
1254+
isDispOK = testLambda(cluE, gamma.m02(), gamma.m20(), mShowerShapeCutValue, mUseNegativeCrossTerm);
1255+
else
1256+
isDispOK = testLambda(cluE, gamma.m20(), gamma.m02(), mShowerShapeCutValue, mUseNegativeCrossTerm);
12271257
for (auto const& match : matches) {
12281258
if (gamma.index() == match.caloClusterId()) {
12291259
matchFlag = true;
@@ -1333,10 +1363,16 @@ struct TpcElIdMassSpectrum {
13331363

13341364
for (auto const& gamma1 : clusters) {
13351365
float cluE1 = gamma1.e();
1336-
if (cluE1 < mMinCluE || cluE1 > mMaxCluE || gamma1.ncell() < mMinCluNcell || gamma1.time() > mMaxCluTime || gamma1.time() < mMinCluTime)
1366+
if (cluE1 < mMinCluE || gamma1.ncell() < mMinCluNcell || gamma1.time() > mMaxCluTime || gamma1.time() < mMinCluTime)
13371367
continue;
13381368
bool matchFlag1 = false;
1339-
if (!testLambda(cluE1, gamma1.m02(), gamma1.m20()))
1369+
1370+
bool isDispOKClu1 = false;
1371+
if (mSwapM20M02ForTestLambda)
1372+
isDispOKClu1 = testLambda(cluE1, gamma1.m02(), gamma1.m20(), mShowerShapeCutValue, mUseNegativeCrossTerm);
1373+
else
1374+
isDispOKClu1 = testLambda(cluE1, gamma1.m20(), gamma1.m02(), mShowerShapeCutValue, mUseNegativeCrossTerm);
1375+
if (!isDispOKClu1)
13401376
continue;
13411377
for (auto const& match : matches) {
13421378
if (gamma1.index() == match.caloClusterId()) {
@@ -1348,9 +1384,14 @@ struct TpcElIdMassSpectrum {
13481384
if (gamma1.index() >= gamma2.index())
13491385
continue;
13501386
float cluE2 = gamma2.e();
1351-
if (cluE2 < mMinCluE || cluE2 > mMaxCluE || gamma2.ncell() < mMinCluNcell || gamma2.time() > mMaxCluTime || gamma2.time() < mMinCluTime)
1387+
if (cluE2 < mMinCluE || gamma2.ncell() < mMinCluNcell || gamma2.time() > mMaxCluTime || gamma2.time() < mMinCluTime)
13521388
continue;
1353-
if (!testLambda(cluE2, gamma2.m02(), gamma2.m20()))
1389+
bool isDispOKClu2 = false;
1390+
if (mSwapM20M02ForTestLambda)
1391+
isDispOKClu2 = testLambda(cluE2, gamma2.m02(), gamma2.m20(), mShowerShapeCutValue, mUseNegativeCrossTerm);
1392+
else
1393+
isDispOKClu2 = testLambda(cluE2, gamma2.m20(), gamma2.m02(), mShowerShapeCutValue, mUseNegativeCrossTerm);
1394+
if (!isDispOKClu2)
13541395
continue;
13551396
bool matchFlag2 = false;
13561397
for (auto const& match : matches) {
@@ -1395,21 +1436,6 @@ struct TpcElIdMassSpectrum {
13951436
}
13961437
}
13971438
}
1398-
1399-
bool testLambda(float pt, float l1, float l2)
1400-
{
1401-
float l2Mean = 1.53126f + 9.50835e+06f / (1.f + 1.08728e+07f * pt + 1.73420e+06f * pt * pt);
1402-
float l1Mean = 1.12365f + 0.123770f * std::exp(-pt * 0.246551f) + 5.30000e-03f * pt;
1403-
float l2Sigma = 6.48260e-02f + 7.60261e+10f / (1.f + 1.53012e+11f * pt + 5.01265e+05f * pt * pt) + 9.00000e-03f * pt;
1404-
float l1Sigma = 4.44719e-04f + 6.99839e-01f / (1.f + 1.22497e+00f * pt + 6.78604e-07f * pt * pt) + 9.00000e-03f * pt;
1405-
float c = -0.35f - 0.550f * std::exp(-0.390730f * pt);
1406-
if (l1Sigma == 0.f || l2Sigma == 0.f)
1407-
return false;
1408-
return 0.5f * (l1 - l1Mean) * (l1 - l1Mean) / (l1Sigma * l1Sigma) +
1409-
0.5f * (l2 - l2Mean) * (l2 - l2Mean) / (l2Sigma * l2Sigma) +
1410-
0.5f * c * (l1 - l1Mean) * (l2 - l2Mean) / (l1Sigma * l2Sigma) <
1411-
4.f;
1412-
}
14131439
};
14141440

14151441
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)