Skip to content

Commit 7c33906

Browse files
dsekihatdsekihat
committed
PWGEM/PhotonMeson: implement searching window for PCM
1 parent 385aec0 commit 7c33906

3 files changed

Lines changed: 218 additions & 110 deletions

File tree

‎PWGEM/PhotonMeson/TableProducer/createEMReducedEvent.cxx‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -46,7 +46,7 @@ struct createEMReducedEvent {
4646

4747
HistogramRegistry registry{"registry"};
4848

49-
Preslice<aod::V0Photons> perCollision_pcm = aod::v0photon::collisionId;
49+
Preslice<aod::V0PhotonsKF> perCollision_pcm = aod::v0photon::collisionId;
5050
Preslice<aod::PHOSClusters> perCollision_phos = aod::skimmedcluster::collisionId;
5151
Preslice<aod::SkimEMCClusters> perCollision_emc = aod::skimmedcluster::collisionId;
5252

@@ -130,7 +130,7 @@ struct createEMReducedEvent {
130130

131131
} // end of process
132132

133-
void process_PCM(MyCollisions const& collisions, aod::BCs const& bcs, aod::V0Photons const& v0photons)
133+
void process_PCM(MyCollisions const& collisions, aod::BCs const& bcs, aod::V0PhotonsKF const& v0photons)
134134
{
135135
process<kPCM>(collisions, bcs, v0photons, nullptr, nullptr);
136136
}
@@ -142,12 +142,12 @@ struct createEMReducedEvent {
142142
{
143143
process<kEMC>(collisions, bcs, nullptr, nullptr, emcclusters);
144144
}
145-
void process_PCM_PHOS(MyCollisions const& collisions, aod::BCs const& bcs, aod::V0Photons const& v0photons, aod::PHOSClusters const& phosclusters)
145+
void process_PCM_PHOS(MyCollisions const& collisions, aod::BCs const& bcs, aod::V0PhotonsKF const& v0photons, aod::PHOSClusters const& phosclusters)
146146
{
147147
const uint8_t sysflag = kPCM | kPHOS;
148148
process<sysflag>(collisions, bcs, v0photons, phosclusters, nullptr);
149149
}
150-
void process_PCM_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::V0Photons const& v0photons, aod::SkimEMCClusters const& emcclusters)
150+
void process_PCM_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::V0PhotonsKF const& v0photons, aod::SkimEMCClusters const& emcclusters)
151151
{
152152
const uint8_t sysflag = kPCM | kEMC;
153153
process<sysflag>(collisions, bcs, v0photons, nullptr, emcclusters);
@@ -157,7 +157,7 @@ struct createEMReducedEvent {
157157
const uint8_t sysflag = kPHOS | kEMC;
158158
process<sysflag>(collisions, bcs, nullptr, phosclusters, emcclusters);
159159
}
160-
void process_PCM_PHOS_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::V0Photons const& v0photons, aod::PHOSClusters const& phosclusters, aod::SkimEMCClusters const& emcclusters)
160+
void process_PCM_PHOS_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::V0PhotonsKF const& v0photons, aod::PHOSClusters const& phosclusters, aod::SkimEMCClusters const& emcclusters)
161161
{
162162
const uint8_t sysflag = kPCM | kPHOS | kEMC;
163163
process<sysflag>(collisions, bcs, v0photons, phosclusters, emcclusters);

‎PWGEM/PhotonMeson/TableProducer/createPCM.cxx‎

Lines changed: 212 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,8 @@
1515
// Please write to: daiki.sekihata@cern.ch
1616

1717
#include <array>
18+
#include <vector>
19+
#include <algorithm>
1820
#include "Framework/runDataProcessing.h"
1921
#include "Framework/AnalysisTask.h"
2022
#include "Framework/AnalysisDataModel.h"
@@ -31,6 +33,7 @@
3133
#include "CCDB/BasicCCDBManager.h"
3234
#include "PWGEM/PhotonMeson/DataModel/gammaTables.h"
3335
#include "PWGLF/DataModel/LFStrangenessTables.h"
36+
#include "PWGEM/PhotonMeson/Utils/PCMUtilities.h"
3437

3538
using namespace o2;
3639
using namespace o2::framework;
@@ -52,6 +55,7 @@ struct createPCM {
5255
"createPCM",
5356
{
5457
{"hEventCounter", "hEventCounter", {HistType::kTH1F, {{5, 0.5f, 5.5f}}}},
58+
// {"hAP", "AP plot", {HistType::kTH2F, {{200, -1, +1}, {250, 0, 0.25}}}},
5559
},
5660
};
5761

@@ -73,7 +77,9 @@ struct createPCM {
7377
Configurable<float> v0Rmax{"v0Rmax", 180.0, "v0Rmax"};
7478
Configurable<float> dcamin{"dcamin", 0.1, "dcamin"};
7579
Configurable<float> dcamax{"dcamax", 1e+10, "dcamax"};
76-
Configurable<float> maxX{"maxX", 90.0, "maximum X (starting point of track X)"}; // maxX is equal to or smaller than minPropagationDistance in trackPropagation.cxx for DCA
80+
Configurable<int> nsw{"nsw", 1, "number of searching window in collisions"};
81+
Configurable<float> maxX{"maxX", 83.1, "maximum X (starting point X of track iu)"};
82+
Configurable<float> maxY{"maxY", 20.0, "maximum Y (starting point Y of track iu)"};
7783
Configurable<float> minpt{"minpt", 0.01, "min pT for single track in GeV/c"};
7884
Configurable<float> maxeta{"maxeta", 0.9, "eta acceptance for single track"};
7985
Configurable<int> mincrossedrows{"mincrossedrows", 10, "min crossed rows"};
@@ -174,8 +180,56 @@ struct createPCM {
174180
}
175181
}
176182

183+
float v0_alpha(float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg)
184+
{
185+
float momTot = RecoDecay::p(pxpos + pxneg, pypos + pyneg, pzpos + pzneg);
186+
float lQlNeg = RecoDecay::dotProd(array{pxneg, pyneg, pzneg}, array{pxpos + pxneg, pypos + pyneg, pzpos + pzneg}) / momTot;
187+
float lQlPos = RecoDecay::dotProd(array{pxpos, pypos, pzpos}, array{pxpos + pxneg, pypos + pyneg, pzpos + pzneg}) / momTot;
188+
return (lQlPos - lQlNeg) / (lQlPos + lQlNeg);
189+
}
190+
float v0_qt(float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg)
191+
{
192+
float momTot = RecoDecay::p2(pxpos + pxneg, pypos + pyneg, pzpos + pzneg);
193+
float dp = RecoDecay::dotProd(array{pxneg, pyneg, pzneg}, array{pxpos + pxneg, pypos + pyneg, pzpos + pzneg});
194+
return std::sqrt(RecoDecay::p2(pxneg, pyneg, pzneg) - dp * dp / momTot);
195+
}
196+
197+
template <typename TTrack>
198+
bool reconstructV0(TTrack const& ele, TTrack const& pos)
199+
{
200+
// fitter is memeber variable.
201+
auto pTrack = getTrackParCov(pos); // positive
202+
auto nTrack = getTrackParCov(ele); // negative
203+
array<float, 3> svpos = {0.}; // secondary vertex position
204+
array<float, 3> pvec0 = {0.};
205+
array<float, 3> pvec1 = {0.};
206+
207+
int nCand = fitter.process(pTrack, nTrack);
208+
if (nCand != 0) {
209+
fitter.propagateTracksToVertex();
210+
const auto& vtx = fitter.getPCACandidate();
211+
for (int i = 0; i < 3; i++) {
212+
svpos[i] = vtx[i];
213+
}
214+
fitter.getTrack(0).getPxPyPzGlo(pvec0); // positive
215+
fitter.getTrack(1).getPxPyPzGlo(pvec1); // negative
216+
} else {
217+
return false;
218+
}
219+
220+
float v0dca = fitter.getChi2AtPCACandidate(); // distance between 2 legs.
221+
if (v0dca > maxdcav0dau) {
222+
return false;
223+
}
224+
225+
if (!checkAP(v0_alpha(pvec0[0], pvec0[1], pvec0[2], pvec1[0], pvec1[1], pvec1[2]), v0_qt(pvec0[0], pvec0[1], pvec0[2], pvec1[0], pvec1[1], pvec1[2]))) { // store only photon conversions
226+
return false;
227+
}
228+
return true;
229+
}
230+
177231
template <typename TCollision, typename TTrack>
178-
void fillV0Table(TCollision const& collision, TTrack const& ele, TTrack const& pos)
232+
void fillV0Table(TCollision const& collision, TTrack const& ele, TTrack const& pos, const bool filltable)
179233
{
180234
array<float, 3> pVtx = {collision.posX(), collision.posY(), collision.posZ()};
181235
array<float, 3> svpos = {0.}; // secondary vertex position
@@ -209,19 +263,31 @@ struct createPCM {
209263
if (v0dca > maxdcav0dau) {
210264
return;
211265
}
212-
if (v0radius < v0Rmin || v0Rmax < v0radius) {
266+
if (v0CosinePA < minv0cospa) {
213267
return;
214268
}
215-
if (v0CosinePA < minv0cospa) {
269+
// registry.fill(HIST("hAP"), v0_alpha(pvec0[0], pvec0[1], pvec0[2], pvec1[0], pvec1[1], pvec1[2]), v0_qt(pvec0[0], pvec0[1], pvec0[2], pvec1[0], pvec1[1], pvec1[2]));
270+
271+
if (!checkAP(v0_alpha(pvec0[0], pvec0[1], pvec0[2], pvec1[0], pvec1[1], pvec1[2]), v0_qt(pvec0[0], pvec0[1], pvec0[2], pvec1[0], pvec1[1], pvec1[2]))) { // store only photon conversions
216272
return;
217273
}
218274

219-
v0data(pos.globalIndex(), ele.globalIndex(), collision.globalIndex(), -1,
220-
fitter.getTrack(0).getX(), fitter.getTrack(1).getX(),
221-
svpos[0], svpos[1], svpos[2],
222-
pvec0[0], pvec0[1], pvec0[2],
223-
pvec1[0], pvec1[1], pvec1[2],
224-
v0dca, pos.dcaXY(), ele.dcaXY());
275+
if (filltable) {
276+
if (v0radius < v0Rmin || v0Rmax < v0radius) {
277+
return;
278+
}
279+
v0data(pos.globalIndex(), ele.globalIndex(), collision.globalIndex(), -1,
280+
fitter.getTrack(0).getX(), fitter.getTrack(1).getX(),
281+
svpos[0], svpos[1], svpos[2],
282+
pvec0[0], pvec0[1], pvec0[2],
283+
pvec1[0], pvec1[1], pvec1[2],
284+
v0dca, pos.dcaXY(), ele.dcaXY());
285+
286+
} else {
287+
// LOGF(info, "storing: collision.globalIndex() = %d , pos.globalIndex() = %d , ele.globalIndex() = %d, cospa = %f", collision.globalIndex(), pos.globalIndex(), ele.globalIndex(), v0CosinePA);
288+
pca_map[std::make_tuple(pos.globalIndex(), ele.globalIndex(), collision.globalIndex())] = v0dca;
289+
cospa_map[std::make_tuple(pos.globalIndex(), ele.globalIndex(), collision.globalIndex())] = v0CosinePA;
290+
} // store indices
225291
}
226292

227293
template <typename TTrack>
@@ -248,32 +314,152 @@ struct createPCM {
248314
return true;
249315
}
250316

251-
Filter trackFilter = o2::aod::track::x < maxX && o2::aod::track::pt > minpt&& nabs(o2::aod::track::eta) < maxeta&& dcamin < nabs(o2::aod::track::dcaXY) && nabs(o2::aod::track::dcaXY) < dcamax;
317+
Filter trackFilter = o2::aod::track::x < maxX && nabs(o2::aod::track::y) < maxY && o2::aod::track::pt > minpt&& nabs(o2::aod::track::eta) < maxeta&& dcamin < nabs(o2::aod::track::dcaXY) && nabs(o2::aod::track::dcaXY) < dcamax && ((min_tpcdEdx < o2::aod::track::tpcSignal && o2::aod::track::tpcSignal < max_tpcdEdx) || o2::aod::track::tpcSignal < -10.f);
252318
using MyFilteredTracks = soa::Filtered<FullTracksExtIU>;
319+
320+
std::map<std::tuple<int32_t, int32_t, int32_t>, float> pca_map;
321+
std::map<std::tuple<int32_t, int32_t, int32_t>, float> cospa_map;
322+
323+
// Partition<MyFilteredTracks> orphan_posTracks = o2::aod::track::signed1Pt > 0.f && o2::aod::track::collisionId < int32_t(0);
324+
// Partition<MyFilteredTracks> orphan_negTracks = o2::aod::track::signed1Pt < 0.f && o2::aod::track::collisionId < int32_t(0);
253325
Partition<MyFilteredTracks> posTracks = o2::aod::track::signed1Pt > 0.f;
254326
Partition<MyFilteredTracks> negTracks = o2::aod::track::signed1Pt < 0.f;
327+
vector<decltype(negTracks->sliceByCached(o2::aod::track::collisionId, 0, cache))> negTracks_sw;
328+
vector<decltype(posTracks->sliceByCached(o2::aod::track::collisionId, 0, cache))> posTracks_sw;
255329

256330
void processSA(MyFilteredTracks const& tracks, aod::Collisions const& collisions, aod::BCsWithTimestamps const&)
257331
{
258-
for (auto& collision : collisions) {
259-
registry.fill(HIST("hEventCounter"), 1);
332+
// LOGF(info, "collisions.size() = %d, tracks.size() = %d", collisions.size(), tracks.size());
333+
for (int64_t icoll = 0; icoll < collisions.size(); icoll += nsw) { // don't repeat the same collision
334+
auto collision = collisions.rawIteratorAt(icoll);
335+
// LOGF(info, "collision.globalIndex() = %d", collision.globalIndex());
260336

261337
auto bc = collision.bc_as<aod::BCsWithTimestamps>();
262338
initCCDB(bc);
339+
// registry.fill(HIST("hEventCounter"), 1);
340+
341+
int32_t min_sw = std::max(int64_t(0), collision.globalIndex());
342+
int32_t max_sw = std::min(int64_t(min_sw + nsw), int64_t(collisions.size()));
343+
344+
// LOGF(info, "orphan_posTracks.size() = %d, orphan_negTracks.size() = %d", orphan_posTracks.size(), orphan_negTracks.size());
345+
negTracks_sw.reserve(max_sw - min_sw);
346+
posTracks_sw.reserve(max_sw - min_sw);
347+
348+
int npos = 0, nneg = 0;
349+
for (int32_t isw = min_sw; isw < max_sw; isw++) {
350+
negTracks_sw.emplace_back(negTracks->sliceByCached(o2::aod::track::collisionId, isw, cache));
351+
posTracks_sw.emplace_back(posTracks->sliceByCached(o2::aod::track::collisionId, isw, cache));
352+
npos += posTracks_sw.back().size();
353+
nneg += negTracks_sw.back().size();
354+
// LOGF(info, "collision.globalIndex() = %d , posTracks_sw.back().size() = %d , negTracks_sw.back().size() = %d", collision.globalIndex(), posTracks_sw.back().size(), negTracks_sw.back().size());
355+
}
356+
// LOGF(info, "min_sw = %d , max_sw = %d , collision.globalIndex() = %d , n posTracks_sw = %d , n negTracks_sw = %d", min_sw, max_sw, collision.globalIndex(), npos, nneg);
357+
358+
for (auto& negTracks_coll : negTracks_sw) {
359+
for (auto& posTracks_coll : posTracks_sw) {
360+
for (auto& [ele, pos] : combinations(CombinationsFullIndexPolicy(negTracks_coll, posTracks_coll))) {
361+
if (!isSelected(ele) || !isSelected(pos)) {
362+
continue;
363+
}
364+
if (!reconstructV0(ele, pos)) { // this is needed for speed-up.
365+
continue;
366+
}
367+
368+
for (int32_t isw = min_sw; isw < max_sw; isw++) {
369+
auto collision_in_sw = collisions.rawIteratorAt(isw);
370+
371+
if (ele.isPVContributor() && isw != ele.collisionId()) {
372+
continue;
373+
}
374+
if (pos.isPVContributor() && isw != pos.collisionId()) {
375+
continue;
376+
}
377+
378+
// if ((ele.isPVContributor() || ele.itsNClsInnerBarrel() > 0) && isw != ele.collisionId()) {
379+
// continue;
380+
// }
381+
// if ((pos.isPVContributor() || pos.itsNClsInnerBarrel() > 0) && isw != pos.collisionId()) {
382+
// continue;
383+
// }
384+
385+
// LOGF(info, "pairing: collision_in_sw.globalIndex() = %d , ele.collisionId() = %d , pos.collisionId() = %d ele.globalIndex() = %d , pos.globalIndex() = %d",
386+
// collision_in_sw.globalIndex(), ele.collisionId(), pos.collisionId(), ele.globalIndex(), pos.globalIndex());
387+
fillV0Table(collision_in_sw, ele, pos, false);
388+
} // end of searching window loop
389+
} // end of pairing loop
390+
} // end of pos track loop in sw
391+
} // end of pos track loop in sw
392+
393+
// LOGF(info, "possible number of V0 = %d", cospa_map.size());
394+
std::map<std::pair<uint32_t, uint32_t>, bool> used_pair_map;
395+
396+
for (const auto& [key, value] : cospa_map) {
397+
auto pos = tracks.rawIteratorAt(std::get<0>(key));
398+
auto ele = tracks.rawIteratorAt(std::get<1>(key));
399+
400+
// LOGF(info, "candidate : pos.globalIndex() = %d , ele.globalIndex() = %d , collision.globalIndex() = %d , cospa = %f , pca = %f", std::get<0>(key), std::get<1>(key), std::get<2>(key), value, pca_map[key]);
401+
402+
std::vector<float> vec_cospa; // vector for each searching window
403+
vec_cospa.reserve(max_sw - min_sw);
404+
for (int32_t isw = min_sw; isw < max_sw; isw++) {
405+
auto collision_in_sw = collisions.rawIteratorAt(isw);
406+
if (cospa_map.find(std::make_tuple(pos.globalIndex(), ele.globalIndex(), collision_in_sw.globalIndex())) != cospa_map.end()) {
407+
vec_cospa.emplace_back(cospa_map[std::make_tuple(pos.globalIndex(), ele.globalIndex(), collision_in_sw.globalIndex())]);
408+
} else {
409+
vec_cospa.emplace_back(-999.f);
410+
}
411+
} // end of searching window loop
412+
413+
// search for the most probable collision where V0 belongs by maximal cospa.
414+
int32_t collision_id_most_prob = std::distance(vec_cospa.begin(), std::max_element(vec_cospa.begin(), vec_cospa.end())) + min_sw;
415+
auto collision_most_prob = collisions.rawIteratorAt(collision_id_most_prob);
416+
// float max_cospa = *std::max_element(vec_cospa.begin(), vec_cospa.end());
417+
// LOGF(info, "max cospa is found! collision_most_prob.globalIndex() = %d , pos.collisionId() = %d , ele.collisionId() = %d, max_cospa = %f", collision_most_prob.globalIndex(), pos.collisionId(), ele.collisionId(), max_cospa);
418+
vec_cospa.clear();
419+
vec_cospa.shrink_to_fit();
420+
421+
// next, check pca between 2 legs in this searching window and select V0s that have the smallest pca to avoid double counting of legs.
422+
float v0pca = pca_map[std::make_tuple(pos.globalIndex(), ele.globalIndex(), collision_most_prob.globalIndex())];
423+
bool is_closest_v0 = true;
424+
for (const auto& [key_tmp, value_tmp] : pca_map) {
425+
auto pos_tmp = tracks.rawIteratorAt(std::get<0>(key_tmp));
426+
auto ele_tmp = tracks.rawIteratorAt(std::get<1>(key_tmp));
427+
428+
float v0pca_tmp = value_tmp;
429+
// float v0pca_tmp = 999.f;
430+
// if(pca_map.find(std::make_tuple(pos_tmp.globalIndex(), ele_tmp.globalIndex(), collision_most_prob.globalIndex())) != pca_map.end()){
431+
// v0pca_tmp = pca_map[std::make_tuple(pos_tmp.globalIndex(), ele_tmp.globalIndex(), collision_most_prob.globalIndex())];
432+
// }
433+
434+
if (ele.globalIndex() == ele_tmp.globalIndex() && pos.globalIndex() == pos_tmp.globalIndex()) { // skip exactly the same V0
435+
continue;
436+
}
437+
if ((ele.globalIndex() == ele_tmp.globalIndex() || pos.globalIndex() == pos_tmp.globalIndex()) && v0pca > v0pca_tmp) {
438+
// LOGF(info, "!reject! | collision id = %d | posid1 = %d , eleid1 = %d , posid2 = %d , eleid2 = %d , pca1 = %f , pca2 = %f",
439+
// collision.globalIndex(), pos.globalIndex(), ele.globalIndex(), pos_tmp.globalIndex(), ele_tmp.globalIndex(), v0pca, v0pca_tmp);
440+
is_closest_v0 = false;
441+
break;
442+
}
443+
} // end of pca_map loop
444+
445+
if (is_closest_v0 && used_pair_map.find(std::make_pair(pos.globalIndex(), ele.globalIndex())) == used_pair_map.end()) {
446+
// LOGF(info, "store : pos.globalIndex() = %d , ele.globalIndex() = %d , collision.globalIndex() = %d , cospa = %f , pca = %f", std::get<0>(key), std::get<1>(key), std::get<2>(key), value, pca_map[key]);
447+
fillV0Table(collision_most_prob, ele, pos, true);
448+
used_pair_map[std::make_pair(pos.globalIndex(), ele.globalIndex())] = true;
449+
}
450+
} // end of pca_map loop
451+
used_pair_map.clear();
263452

264-
auto negTracks_coll = negTracks->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache);
265-
auto posTracks_coll = posTracks->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache);
266-
267-
// LOGF(info, "collision.globalIndex() = %d , negTracks_coll.size() = %d , posTracks_coll.size() = %d", collision.globalIndex(), negTracks_coll.size(), posTracks_coll.size());
453+
pca_map.clear();
454+
cospa_map.clear();
268455

269-
for (auto& [ele, pos] : combinations(CombinationsFullIndexPolicy(negTracks_coll, posTracks_coll))) {
270-
if (!isSelected(ele) || !isSelected(pos)) {
271-
continue;
272-
}
273-
fillV0Table(collision, ele, pos);
274-
}
456+
negTracks_sw.clear();
457+
posTracks_sw.clear();
458+
negTracks_sw.shrink_to_fit();
459+
posTracks_sw.shrink_to_fit();
275460
} // end of collision loop
276-
} // end of process
461+
462+
} // end of process
277463
PROCESS_SWITCH(createPCM, processSA, "create V0s with stand-alone way", true);
278464

279465
Preslice<aod::TrackAssoc> trackIndicesPerCollision = aod::track_association::collisionId;
@@ -301,9 +487,9 @@ struct createPCM {
301487
}
302488

303489
if (ele.sign() < 0) {
304-
fillV0Table(collision, ele, pos);
490+
fillV0Table(collision, ele, pos, true);
305491
} else {
306-
fillV0Table(collision, pos, ele);
492+
fillV0Table(collision, pos, ele, true);
307493
}
308494
}
309495
} // end of collision loop

0 commit comments

Comments
 (0)