diff --git a/PWGEM/PhotonMeson/Core/HistogramsLibrary.cxx b/PWGEM/PhotonMeson/Core/HistogramsLibrary.cxx index 01339f98502..0ad7d2e5139 100644 --- a/PWGEM/PhotonMeson/Core/HistogramsLibrary.cxx +++ b/PWGEM/PhotonMeson/Core/HistogramsLibrary.cxx @@ -60,7 +60,7 @@ void o2::aod::emphotonhistograms::DefineHistograms(THashList* list, const char* list->Add(new TH1F("hNclsITS", "number of ITS clusters", 8, -0.5, 7.5)); list->Add(new TH1F("hChi2ITS", "chi2/number of ITS clusters", 360, 0, 36)); list->Add(new TH2F("hXY", "X vs. Y;X;Y", 100, 0, 100, 100, -50, 50)); - list->Add(new TH2F("hZX", "Z vs. X;Z;X", 1000, -100, 100, 100, 0, 100)); + list->Add(new TH2F("hZX", "Z vs. X;Z;X", 400, -100, 100, 100, 0, 100)); list->Add(new TH2F("hZY", "Z vs. Y;Z;Y", 200, -100, 100, 100, -50, 50)); list->Add(new TH2F("hDCAxyEta", "DCAxy vs. #eta;#eta;DCA_{xy} (cm)", 400, -2, +2, 100, -50, 50)); list->Add(new TH2F("hDCAxyZ", "DCAxy vs. Z;Z (cm);DCA_{xy} (cm)", 200, -100, +100, 100, -50, 50)); @@ -85,6 +85,7 @@ void o2::aod::emphotonhistograms::DefineHistograms(THashList* list, const char* list->Add(new TH2F("hKFChi2vsY", "KF chi2 vs. recalc. conversion point in Y;Y (cm);KF chi2/NDF", 400, -200.0f, 200.0f, 100, 0.f, 100.0f)); list->Add(new TH2F("hKFChi2vsZ", "KF chi2 vs. recalc. conversion point in Z;Z (cm);KF chi2/NDF", 500, -250.0f, 250.0f, 100, 0.f, 100.0f)); list->Add(new TH1F("hNgamma", "Number of #gamma candidates per collision", 101, -0.5f, 100.5f)); + list->Add(new TH2F("hV0R_minTrackX", "V0R vs. trackiu x;R_{xy} of V0 (cm);min TrackIU X (cm)", 100, 0.0f, 100.0f, 200, 0.f, 100.0f)); const int nrxy = 102; double rxy[nrxy] = {0.f}; diff --git a/PWGEM/PhotonMeson/Core/V0PhotonCut.h b/PWGEM/PhotonMeson/Core/V0PhotonCut.h index 8c76b69c2e4..505fbea4e96 100644 --- a/PWGEM/PhotonMeson/Core/V0PhotonCut.h +++ b/PWGEM/PhotonMeson/Core/V0PhotonCut.h @@ -116,9 +116,13 @@ class V0PhotonCut : public TNamed auto pos = v0.template posTrack_as(); auto ele = v0.template negTrack_as(); - // float pos_rxy = sqrt(pow(pos.x(),2) + pow(pos.y(),2)); - // float ele_rxy = sqrt(pow(ele.x(),2) + pow(ele.y(),2)); - // if (v0.recalculatedVtxR() > std::min(pos_rxy, ele_rxy)) { + // if(pos.hasITS() && ele.hasITS()){ + // bool ret = sqrt(pow(v0.alpha() / 0.95, 2) + pow(v0.qtarm() / 0.02, 2) ) < 1.0; + // if (!ret){ + // return false; + // } + // } + // if (v0.recalculatedVtxR() > std::min(pos.x(), ele.x()) + 7.f && (pos.x() > 1.f && ele.x() > 1.f)) { // return false; // } diff --git a/PWGEM/PhotonMeson/TableProducer/CMakeLists.txt b/PWGEM/PhotonMeson/TableProducer/CMakeLists.txt index d6f8e937231..84737f85461 100644 --- a/PWGEM/PhotonMeson/TableProducer/CMakeLists.txt +++ b/PWGEM/PhotonMeson/TableProducer/CMakeLists.txt @@ -53,3 +53,8 @@ o2physics_add_dpl_workflow(produce-meson-calo SOURCES produceMesonCalo.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(check-v0-mc + SOURCES checkV0MC.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::DCAFitter O2Physics::AnalysisCore + COMPONENT_NAME Analysis) diff --git a/PWGEM/PhotonMeson/TableProducer/checkV0MC.cxx b/PWGEM/PhotonMeson/TableProducer/checkV0MC.cxx new file mode 100644 index 00000000000..dd77bb8580d --- /dev/null +++ b/PWGEM/PhotonMeson/TableProducer/checkV0MC.cxx @@ -0,0 +1,155 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \brief check MC true V0s +/// \author daiki.sekihata@cern.ch + +#include +#include +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "Common/Core/trackUtilities.h" +#include "PWGEM/PhotonMeson/DataModel/gammaTables.h" +#include "PWGEM/PhotonMeson/Utils/PCMUtilities.h" +#include "PWGEM/PhotonMeson/Utils/MCUtilities.h" + +using namespace o2; +using namespace o2::soa; +using namespace o2::framework; +using namespace o2::framework::expressions; + +using MyTracks = soa::Join; +using MyTracksMC = soa::Join; + +struct checkV0MC { + + Configurable maxX{"maxX", 83.1, "maximum X (starting point of track X)"}; + Configurable maxY{"maxY", 20.0, "maximum X (starting point of track X)"}; + Configurable minpt{"minpt", 0.01, "min pt for track"}; + Configurable maxeta{"maxeta", 0.9, "eta acceptance"}; + Configurable dcamin{"dcamin", 0.1, "dcamin"}; + Configurable dcamax{"dcamax", 1e+10, "dcamax"}; + + HistogramRegistry fRegistry{ + "fRegistry", + { + {"hDiffCollId", "difference in collision id between ele/pos", {HistType::kTH1F, {{101, -50.5f, +50.5f}}}}, + {"hOrphanTrack", "Number of orphan track from true photon conversions", {HistType::kTH1F, {{1, 0.5f, +1.5f}}}}, + {"hV0PhotonCorrX", "V0 photon track iu X correlation;track iu X for ele;track iu X for pos", {HistType::kTH2F, {{100, 0, +100}, {100, 0, +100}}}}, + {"hV0PhotonCorrY", "V0 photon track iu Y correlation;track iu Y for ele;track iu Y diff", {HistType::kTH2F, {{200, -100, +100}, {400, -20, +20}}}}, + {"hV0PhotonCorrZ", "V0 photon track iu Z correlation;track iu Z for ele;track iu Z diff", {HistType::kTH2F, {{200, -100, +100}, {400, -20, +20}}}}, + {"hV0PhotonCorrYvsX", "V0 photon track iu Y correlation;track iu Y for ele;track iu Y diff", {HistType::kTH2F, {{100, 0, +100}, {400, -20, +20}}}}, + {"hV0PhotonCorrZvsX", "V0 photon track iu Z correlation;track iu Z for ele;track iu Z diff", {HistType::kTH2F, {{100, 0, +100}, {400, -20, +20}}}}, + }, + }; + + template + bool isSelected(TTrack const& track) + { + if (track.pt() < minpt || abs(track.eta()) > maxeta) { + return false; + } + if (abs(track.dcaXY()) < dcamin || dcamax < abs(track.dcaXY())) { + return false; + } + if (!track.hasITS() && !track.hasTPC()) { + return false; + } + + if (track.hasTPC()) { + if (track.tpcSignal() < 40.f || 110.f < track.tpcSignal()) { + return false; + } + } + return true; + } + + 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; + using MyFilteredTracksMC = soa::Filtered; + Partition posTracks = o2::aod::track::signed1Pt > 0.f; + Partition negTracks = o2::aod::track::signed1Pt < 0.f; + Partition orphan_posTracks = o2::aod::track::signed1Pt > 0.f && o2::aod::track::collisionId < 0; + Partition orphan_negTracks = o2::aod::track::signed1Pt < 0.f && o2::aod::track::collisionId < 0; + + int ndf = 0; + + using MyCollisions = soa::Join; + void processMC(MyCollisions const& collisions, aod::McCollisions const&, aod::BCsWithTimestamps const& bcs, MyFilteredTracksMC const& tracks, aod::McParticles const& mcTracks) + { + LOGF(info, "n pos tracks = %d , n neg tracks = %d", posTracks.size(), negTracks.size()); + ndf++; + if (ndf > 3) { // 3 DFs are enough. + return; + } + + for (auto& [ele, pos] : combinations(CombinationsFullIndexPolicy(negTracks, posTracks))) { + if (!ele.has_mcParticle() || !pos.has_mcParticle()) { + continue; + } + + if (!isSelected(ele) || !isSelected(pos)) { + continue; + } + + auto elemc = ele.template mcParticle_as(); + auto posmc = pos.template mcParticle_as(); + + int photonid = FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 22, mcTracks); + if (photonid < 0) { // check swap, true electron is reconstructed as positron and vice versa. + photonid = FindCommonMotherFrom2Prongs(posmc, elemc, 11, -11, 22, mcTracks); + } + if (photonid < 0) { + continue; + } + // At this point, there are only true photon conversions. + auto photonmc = mcTracks.iteratorAt(photonid); + if (!IsPhysicalPrimary(photonmc.mcCollision(), photonmc, mcTracks)) { + continue; + } + + bool has_coll_ele = ele.has_collision(); + bool has_coll_pos = pos.has_collision(); + + if (!has_coll_ele) { + fRegistry.fill(HIST("hOrphanTrack"), 1); + } + if (!has_coll_pos) { + fRegistry.fill(HIST("hOrphanTrack"), 1); + } + + if (!has_coll_ele || !has_coll_pos) { + continue; + } + + auto collision_ele = ele.collision(); + auto collision_pos = pos.collision(); + int diff = collision_ele.globalIndex() - collision_pos.globalIndex(); + fRegistry.fill(HIST("hDiffCollId"), diff); + + fRegistry.fill(HIST("hV0PhotonCorrX"), ele.x(), pos.x()); + fRegistry.fill(HIST("hV0PhotonCorrY"), ele.y(), pos.y() - ele.y()); + fRegistry.fill(HIST("hV0PhotonCorrZ"), ele.z(), pos.z() - ele.z()); + fRegistry.fill(HIST("hV0PhotonCorrYvsX"), ele.x(), pos.y() - ele.y()); + fRegistry.fill(HIST("hV0PhotonCorrZvsX"), ele.x(), pos.z() - ele.z()); + + } // end of pairing loop + } + PROCESS_SWITCH(checkV0MC, processMC, "process mc truth info", true); + + void processDummy(soa::Join const& collisions) {} + PROCESS_SWITCH(checkV0MC, processDummy, "process dummy", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc, TaskName{"check-v0-mc"})}; +} diff --git a/PWGEM/PhotonMeson/TableProducer/createEMReducedEvent.cxx b/PWGEM/PhotonMeson/TableProducer/createEMReducedEvent.cxx index 52624c541ff..4b1f6574ada 100644 --- a/PWGEM/PhotonMeson/TableProducer/createEMReducedEvent.cxx +++ b/PWGEM/PhotonMeson/TableProducer/createEMReducedEvent.cxx @@ -46,7 +46,7 @@ struct createEMReducedEvent { HistogramRegistry registry{"registry"}; - Preslice perCollision_pcm = aod::v0photon::collisionId; + Preslice perCollision_pcm = aod::v0photon::collisionId; Preslice perCollision_phos = aod::skimmedcluster::collisionId; Preslice perCollision_emc = aod::skimmedcluster::collisionId; @@ -130,7 +130,7 @@ struct createEMReducedEvent { } // end of process - void process_PCM(MyCollisions const& collisions, aod::BCs const& bcs, aod::V0Photons const& v0photons) + void process_PCM(MyCollisions const& collisions, aod::BCs const& bcs, aod::V0PhotonsKF const& v0photons) { process(collisions, bcs, v0photons, nullptr, nullptr); } @@ -142,12 +142,12 @@ struct createEMReducedEvent { { process(collisions, bcs, nullptr, nullptr, emcclusters); } - void process_PCM_PHOS(MyCollisions const& collisions, aod::BCs const& bcs, aod::V0Photons const& v0photons, aod::PHOSClusters const& phosclusters) + void process_PCM_PHOS(MyCollisions const& collisions, aod::BCs const& bcs, aod::V0PhotonsKF const& v0photons, aod::PHOSClusters const& phosclusters) { const uint8_t sysflag = kPCM | kPHOS; process(collisions, bcs, v0photons, phosclusters, nullptr); } - void process_PCM_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::V0Photons const& v0photons, aod::SkimEMCClusters const& emcclusters) + void process_PCM_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::V0PhotonsKF const& v0photons, aod::SkimEMCClusters const& emcclusters) { const uint8_t sysflag = kPCM | kEMC; process(collisions, bcs, v0photons, nullptr, emcclusters); @@ -157,7 +157,7 @@ struct createEMReducedEvent { const uint8_t sysflag = kPHOS | kEMC; process(collisions, bcs, nullptr, phosclusters, emcclusters); } - void process_PCM_PHOS_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::V0Photons const& v0photons, aod::PHOSClusters const& phosclusters, aod::SkimEMCClusters const& emcclusters) + void process_PCM_PHOS_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::V0PhotonsKF const& v0photons, aod::PHOSClusters const& phosclusters, aod::SkimEMCClusters const& emcclusters) { const uint8_t sysflag = kPCM | kPHOS | kEMC; process(collisions, bcs, v0photons, phosclusters, emcclusters); diff --git a/PWGEM/PhotonMeson/TableProducer/createPCM.cxx b/PWGEM/PhotonMeson/TableProducer/createPCM.cxx index fced55e42ee..60a6298659e 100644 --- a/PWGEM/PhotonMeson/TableProducer/createPCM.cxx +++ b/PWGEM/PhotonMeson/TableProducer/createPCM.cxx @@ -15,6 +15,8 @@ // Please write to: daiki.sekihata@cern.ch #include +#include +#include #include "Framework/runDataProcessing.h" #include "Framework/AnalysisTask.h" #include "Framework/AnalysisDataModel.h" @@ -31,6 +33,7 @@ #include "CCDB/BasicCCDBManager.h" #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" #include "PWGLF/DataModel/LFStrangenessTables.h" +#include "PWGEM/PhotonMeson/Utils/PCMUtilities.h" using namespace o2; using namespace o2::framework; @@ -52,6 +55,7 @@ struct createPCM { "createPCM", { {"hEventCounter", "hEventCounter", {HistType::kTH1F, {{5, 0.5f, 5.5f}}}}, + // {"hAP", "AP plot", {HistType::kTH2F, {{200, -1, +1}, {250, 0, 0.25}}}}, }, }; @@ -73,13 +77,16 @@ struct createPCM { Configurable v0Rmax{"v0Rmax", 180.0, "v0Rmax"}; Configurable dcamin{"dcamin", 0.1, "dcamin"}; Configurable dcamax{"dcamax", 1e+10, "dcamax"}; - Configurable maxX{"maxX", 90.0, "maximum X (starting point of track X)"}; // maxX is equal to or smaller than minPropagationDistance in trackPropagation.cxx for DCA + Configurable nsw{"nsw", 1, "number of searching window in collisions"}; + Configurable maxX{"maxX", 83.1, "maximum X (starting point X of track iu)"}; + Configurable maxY{"maxY", 20.0, "maximum Y (starting point Y of track iu)"}; Configurable minpt{"minpt", 0.01, "min pT for single track in GeV/c"}; Configurable maxeta{"maxeta", 0.9, "eta acceptance for single track"}; Configurable mincrossedrows{"mincrossedrows", 10, "min crossed rows"}; Configurable maxchi2tpc{"maxchi2tpc", 4.0, "max chi2/NclsTPC"}; Configurable min_tpcdEdx{"min_tpcdEdx", 30.0, "min TPC dE/dx"}; Configurable max_tpcdEdx{"max_tpcdEdx", 110.0, "max TPC dE/dx"}; + Configurable margin_r{"margin_r", 7.0, "margin for r cut"}; int mRunNumber; float d_bz; @@ -174,18 +181,68 @@ struct createPCM { } } + float v0_alpha(float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg) + { + float momTot = RecoDecay::p(pxpos + pxneg, pypos + pyneg, pzpos + pzneg); + float lQlNeg = RecoDecay::dotProd(array{pxneg, pyneg, pzneg}, array{pxpos + pxneg, pypos + pyneg, pzpos + pzneg}) / momTot; + float lQlPos = RecoDecay::dotProd(array{pxpos, pypos, pzpos}, array{pxpos + pxneg, pypos + pyneg, pzpos + pzneg}) / momTot; + return (lQlPos - lQlNeg) / (lQlPos + lQlNeg); + } + float v0_qt(float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg) + { + float momTot = RecoDecay::p2(pxpos + pxneg, pypos + pyneg, pzpos + pzneg); + float dp = RecoDecay::dotProd(array{pxneg, pyneg, pzneg}, array{pxpos + pxneg, pypos + pyneg, pzpos + pzneg}); + return std::sqrt(RecoDecay::p2(pxneg, pyneg, pzneg) - dp * dp / momTot); + } + template - bool IsTPConlyTrack(TTrack const& track) + bool reconstructV0(TTrack const& ele, TTrack const& pos) { - if (track.hasTPC() && (!track.hasITS() && !track.hasTOF() && !track.hasTRD())) { - return true; + // fitter is memeber variable. + auto pTrack = getTrackParCov(pos); // positive + auto nTrack = getTrackParCov(ele); // negative + array svpos = {0.}; // secondary vertex position + array pvec0 = {0.}; + array pvec1 = {0.}; + + int nCand = fitter.process(pTrack, nTrack); + if (nCand != 0) { + fitter.propagateTracksToVertex(); + const auto& vtx = fitter.getPCACandidate(); + for (int i = 0; i < 3; i++) { + svpos[i] = vtx[i]; + } + fitter.getTrack(0).getPxPyPzGlo(pvec0); // positive + fitter.getTrack(1).getPxPyPzGlo(pvec1); // negative } else { return false; } + + float v0dca = fitter.getChi2AtPCACandidate(); // distance between 2 legs. + if (v0dca > maxdcav0dau) { + return false; + } + + 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]), 0.95, 0.05)) { // store only photon conversions + return false; + } + if (ele.hasITS() && pos.hasITS() && !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]), 0.95, 0.02)) { // store only photon conversions + return false; + } + + float xyz[3] = {0.f, 0.f, 0.f}; + Vtx_recalculation(o2::base::Propagator::Instance(), pos, ele, xyz); + float recalculatedVtxR = sqrt(pow(xyz[0], 2) + pow(xyz[1], 2)); + // LOGF(info, "recalculated vtx : x = %f , y = %f , z = %f", xyz[0], xyz[1], xyz[2]); + if (recalculatedVtxR > std::min(pos.x(), ele.x()) + margin_r && (pos.x() > 1.f && ele.x() > 1.f)) { + return false; + } + + return true; } template - void fillV0Table(TCollision const& collision, TTrack const& ele, TTrack const& pos) + void fillV0Table(TCollision const& collision, TTrack const& ele, TTrack const& pos, const bool filltable) { array pVtx = {collision.posX(), collision.posY(), collision.posZ()}; array svpos = {0.}; // secondary vertex position @@ -219,19 +276,31 @@ struct createPCM { if (v0dca > maxdcav0dau) { return; } - if (v0radius < v0Rmin || v0Rmax < v0radius) { + if (v0CosinePA < minv0cospa) { return; } - if (v0CosinePA < minv0cospa) { + // 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])); + + 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 return; } - v0data(pos.globalIndex(), ele.globalIndex(), collision.globalIndex(), -1, - fitter.getTrack(0).getX(), fitter.getTrack(1).getX(), - svpos[0], svpos[1], svpos[2], - pvec0[0], pvec0[1], pvec0[2], - pvec1[0], pvec1[1], pvec1[2], - v0dca, pos.dcaXY(), ele.dcaXY()); + if (filltable) { + if (v0radius < v0Rmin || v0Rmax < v0radius) { + return; + } + v0data(pos.globalIndex(), ele.globalIndex(), collision.globalIndex(), -1, + fitter.getTrack(0).getX(), fitter.getTrack(1).getX(), + svpos[0], svpos[1], svpos[2], + pvec0[0], pvec0[1], pvec0[2], + pvec1[0], pvec1[1], pvec1[2], + v0dca, pos.dcaXY(), ele.dcaXY()); + + } else { + // LOGF(info, "storing: collision.globalIndex() = %d , pos.globalIndex() = %d , ele.globalIndex() = %d, cospa = %f", collision.globalIndex(), pos.globalIndex(), ele.globalIndex(), v0CosinePA); + pca_map[std::make_tuple(pos.globalIndex(), ele.globalIndex(), collision.globalIndex())] = v0dca; + cospa_map[std::make_tuple(pos.globalIndex(), ele.globalIndex(), collision.globalIndex())] = v0CosinePA; + } // store indices } template @@ -258,38 +327,145 @@ struct createPCM { return true; } - 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; + 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); using MyFilteredTracks = soa::Filtered; + + std::map, float> pca_map; + std::map, float> cospa_map; + + // Partition orphan_posTracks = o2::aod::track::signed1Pt > 0.f && o2::aod::track::collisionId < int32_t(0); + // Partition orphan_negTracks = o2::aod::track::signed1Pt < 0.f && o2::aod::track::collisionId < int32_t(0); Partition posTracks = o2::aod::track::signed1Pt > 0.f; Partition negTracks = o2::aod::track::signed1Pt < 0.f; + vectorsliceByCached(o2::aod::track::collisionId, 0, cache))> negTracks_sw; + vectorsliceByCached(o2::aod::track::collisionId, 0, cache))> posTracks_sw; void processSA(MyFilteredTracks const& tracks, aod::Collisions const& collisions, aod::BCsWithTimestamps const&) { - for (auto& collision : collisions) { - registry.fill(HIST("hEventCounter"), 1); + // LOGF(info, "collisions.size() = %d, tracks.size() = %d", collisions.size(), tracks.size()); + for (int64_t icoll = 0; icoll < collisions.size(); icoll += nsw) { // don't repeat the same collision + auto collision = collisions.rawIteratorAt(icoll); + // LOGF(info, "collision.globalIndex() = %d", collision.globalIndex()); auto bc = collision.bc_as(); initCCDB(bc); + // registry.fill(HIST("hEventCounter"), 1); + + int32_t min_sw = std::max(int64_t(0), collision.globalIndex()); + int32_t max_sw = std::min(int64_t(min_sw + nsw), int64_t(collisions.size())); + + // LOGF(info, "orphan_posTracks.size() = %d, orphan_negTracks.size() = %d", orphan_posTracks.size(), orphan_negTracks.size()); + negTracks_sw.reserve(max_sw - min_sw); + posTracks_sw.reserve(max_sw - min_sw); + + int npos = 0, nneg = 0; + for (int32_t isw = min_sw; isw < max_sw; isw++) { + negTracks_sw.emplace_back(negTracks->sliceByCached(o2::aod::track::collisionId, isw, cache)); + posTracks_sw.emplace_back(posTracks->sliceByCached(o2::aod::track::collisionId, isw, cache)); + npos += posTracks_sw.back().size(); + nneg += negTracks_sw.back().size(); + // 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()); + } + // 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); + + for (auto& negTracks_coll : negTracks_sw) { + for (auto& posTracks_coll : posTracks_sw) { + for (auto& [ele, pos] : combinations(CombinationsFullIndexPolicy(negTracks_coll, posTracks_coll))) { + if (!isSelected(ele) || !isSelected(pos)) { + continue; + } + if (!reconstructV0(ele, pos)) { // this is needed for speed-up. + continue; + } + + for (int32_t isw = min_sw; isw < max_sw; isw++) { + auto collision_in_sw = collisions.rawIteratorAt(isw); + + if (ele.isPVContributor() && isw != ele.collisionId()) { + continue; + } + if (pos.isPVContributor() && isw != pos.collisionId()) { + continue; + } + + // LOGF(info, "pairing: collision_in_sw.globalIndex() = %d , ele.collisionId() = %d , pos.collisionId() = %d ele.globalIndex() = %d , pos.globalIndex() = %d", + // collision_in_sw.globalIndex(), ele.collisionId(), pos.collisionId(), ele.globalIndex(), pos.globalIndex()); + fillV0Table(collision_in_sw, ele, pos, false); + } // end of searching window loop + } // end of pairing loop + } // end of pos track loop in sw + } // end of pos track loop in sw + + // LOGF(info, "possible number of V0 = %d", cospa_map.size()); + std::map, bool> used_pair_map; + + for (const auto& [key, value] : cospa_map) { + auto pos = tracks.rawIteratorAt(std::get<0>(key)); + auto ele = tracks.rawIteratorAt(std::get<1>(key)); + + // 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]); + + std::vector vec_cospa; // vector for each searching window + vec_cospa.reserve(max_sw - min_sw); + for (int32_t isw = min_sw; isw < max_sw; isw++) { + auto collision_in_sw = collisions.rawIteratorAt(isw); + if (cospa_map.find(std::make_tuple(pos.globalIndex(), ele.globalIndex(), collision_in_sw.globalIndex())) != cospa_map.end()) { + vec_cospa.emplace_back(cospa_map[std::make_tuple(pos.globalIndex(), ele.globalIndex(), collision_in_sw.globalIndex())]); + } else { + vec_cospa.emplace_back(-999.f); + } + } // end of searching window loop + + // search for the most probable collision where V0 belongs by maximal cospa. + int32_t collision_id_most_prob = std::distance(vec_cospa.begin(), std::max_element(vec_cospa.begin(), vec_cospa.end())) + min_sw; + auto collision_most_prob = collisions.rawIteratorAt(collision_id_most_prob); + // float max_cospa = *std::max_element(vec_cospa.begin(), vec_cospa.end()); + // 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); + vec_cospa.clear(); + vec_cospa.shrink_to_fit(); + + // next, check pca between 2 legs in this searching window and select V0s that have the smallest pca to avoid double counting of legs. + float v0pca = pca_map[std::make_tuple(pos.globalIndex(), ele.globalIndex(), collision_most_prob.globalIndex())]; + bool is_closest_v0 = true; + for (const auto& [key_tmp, value_tmp] : pca_map) { + auto pos_tmp = tracks.rawIteratorAt(std::get<0>(key_tmp)); + auto ele_tmp = tracks.rawIteratorAt(std::get<1>(key_tmp)); + + float v0pca_tmp = value_tmp; + // float v0pca_tmp = 999.f; + // if(pca_map.find(std::make_tuple(pos_tmp.globalIndex(), ele_tmp.globalIndex(), collision_most_prob.globalIndex())) != pca_map.end()){ + // v0pca_tmp = pca_map[std::make_tuple(pos_tmp.globalIndex(), ele_tmp.globalIndex(), collision_most_prob.globalIndex())]; + // } + + if (ele.globalIndex() == ele_tmp.globalIndex() && pos.globalIndex() == pos_tmp.globalIndex()) { // skip exactly the same V0 + continue; + } + if ((ele.globalIndex() == ele_tmp.globalIndex() || pos.globalIndex() == pos_tmp.globalIndex()) && v0pca > v0pca_tmp) { + // LOGF(info, "!reject! | collision id = %d | posid1 = %d , eleid1 = %d , posid2 = %d , eleid2 = %d , pca1 = %f , pca2 = %f", + // collision.globalIndex(), pos.globalIndex(), ele.globalIndex(), pos_tmp.globalIndex(), ele_tmp.globalIndex(), v0pca, v0pca_tmp); + is_closest_v0 = false; + break; + } + } // end of pca_map loop + + if (is_closest_v0 && used_pair_map.find(std::make_pair(pos.globalIndex(), ele.globalIndex())) == used_pair_map.end()) { + // 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]); + fillV0Table(collision_most_prob, ele, pos, true); + used_pair_map[std::make_pair(pos.globalIndex(), ele.globalIndex())] = true; + } + } // end of pca_map loop + used_pair_map.clear(); - auto negTracks_coll = negTracks->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache); - auto posTracks_coll = posTracks->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache); - - // LOGF(info, "collision.globalIndex() = %d , negTracks_coll.size() = %d , posTracks_coll.size() = %d", collision.globalIndex(), negTracks_coll.size(), posTracks_coll.size()); + pca_map.clear(); + cospa_map.clear(); - for (auto& [ele, pos] : combinations(CombinationsFullIndexPolicy(negTracks_coll, posTracks_coll))) { - if (!isSelected(ele) || !isSelected(pos)) { - continue; - } - // if(ele.hasITS() && (!ele.hasTPC() && !ele.hasTOF() && !ele.hasTRD())){ - // LOGF(info, "ele : ITSonly track %d", ele.globalIndex()); - // } - // if(pos.hasITS() && (!pos.hasTPC() && !pos.hasTOF() && !pos.hasTRD())){ - // LOGF(info, "pos : ITSonly track %d", pos.globalIndex()); - // } - fillV0Table(collision, ele, pos); - } + negTracks_sw.clear(); + posTracks_sw.clear(); + negTracks_sw.shrink_to_fit(); + posTracks_sw.shrink_to_fit(); } // end of collision loop - } // end of process + + } // end of process PROCESS_SWITCH(createPCM, processSA, "create V0s with stand-alone way", true); Preslice trackIndicesPerCollision = aod::track_association::collisionId; @@ -317,9 +493,9 @@ struct createPCM { } if (ele.sign() < 0) { - fillV0Table(collision, ele, pos); + fillV0Table(collision, ele, pos, true); } else { - fillV0Table(collision, pos, ele); + fillV0Table(collision, pos, ele, true); } } } // end of collision loop diff --git a/PWGEM/PhotonMeson/TableProducer/skimmerGammaConversions.cxx b/PWGEM/PhotonMeson/TableProducer/skimmerGammaConversions.cxx index a1fa08a60db..3ae1850283d 100644 --- a/PWGEM/PhotonMeson/TableProducer/skimmerGammaConversions.cxx +++ b/PWGEM/PhotonMeson/TableProducer/skimmerGammaConversions.cxx @@ -235,30 +235,6 @@ struct skimmerGammaConversions { return true; } - // template - // bool checkV0leg(TLeg const& leg) - //{ - // if (abs(leg.eta()) > maxeta) { - // return false; - // } - // if (abs(leg.tpcNSigmaEl()) > maxTPCNsigmaEl) { - // return false; - // } - // if (leg.tpcChi2NCl() > maxchi2tpc) { - // return false; - // } - // if (leg.tpcNClsCrossedRows() < mincrossedrows) { - // return false; - // } - // if (abs(leg.dcaXY()) < dcamin) { - // return false; - // } - // if (dcamax < abs(leg.dcaXY())) { - // return false; - // } - // return true; - // } - template void fillV0KF(TCollision const& collision, TV0 const& v0, recalculatedVertexParameters recalculatedVertex) { @@ -315,7 +291,7 @@ struct skimmerGammaConversions { // ============================ FUNCTION DEFINITIONS ==================================================== - Preslice perCollision = aod::v0data::collisionId; + PresliceUnsorted perCollision = aod::v0data::collisionId; void processRec(aod::Collisions const& collisions, aod::BCsWithTimestamps const& bcs, @@ -338,35 +314,13 @@ struct skimmerGammaConversions { continue; } - bool flag_closer = true; - for (auto& v0tmp : groupedV0s) { - if (!checkAP(v0tmp.alpha(), v0tmp.qtarm())) { // store only photon conversions - continue; - } - auto pos_tmp = v0tmp.template posTrack_as(); // positive daughter - auto ele_tmp = v0tmp.template negTrack_as(); // negative daughter - if (!checkV0leg(pos_tmp) || !checkV0leg(ele_tmp)) { - continue; - } - - if (v0.index() == v0tmp.index()) { // don't check onviously, exactly the same v0. - // LOGF(info, "don't check the exactly the same 2 V0s"); - continue; - } - if ((ele.globalIndex() == ele_tmp.globalIndex() || pos.globalIndex() == pos_tmp.globalIndex()) && v0.dcaV0daughters() > v0tmp.dcaV0daughters()) { - // LOGF(info, "!reject! | collision id = %d | g1 id = %d , g2 id = %d , posid1 = %d , eleid1 = %d , posid2 = %d , eleid2 = %d , pca1 = %f , pca2 = %f", - // collision.globalIndex(), v0.index(), v0tmp.index(), pos.globalIndex(), ele.globalIndex(), pos_tmp.globalIndex(), ele_tmp.globalIndex(), v0.dcaV0daughters(), v0tmp.dcaV0daughters()); - flag_closer = false; - break; - } - } // end of v0tmp loop - - if (!flag_closer) { - continue; - } - + float xyz[3] = {0.f, 0.f, 0.f}; + Vtx_recalculation(o2::base::Propagator::Instance(), pos, ele, xyz); recalculatedVertexParameters recalculatedVertex; - Vtx_recalculation(pos, ele, &recalculatedVertex); + recalculatedVertex.recalculatedConversionPoint[0] = xyz[0]; + recalculatedVertex.recalculatedConversionPoint[1] = xyz[1]; + recalculatedVertex.recalculatedConversionPoint[2] = xyz[2]; + v0photons(collision.globalIndex(), v0legs.lastIndex() + 1, v0legs.lastIndex() + 2, v0.x(), v0.y(), v0.z(), v0.pxpos(), v0.pypos(), v0.pzpos(), @@ -418,43 +372,19 @@ struct skimmerGammaConversions { continue; } - bool flag_closer = true; - for (auto& v0tmp : lGroupedV0s) { - if (!checkAP(v0tmp.alpha(), v0tmp.qtarm())) { // store only photon conversions - continue; - } - auto pos_tmp = v0tmp.template posTrack_as(); // positive daughter - auto ele_tmp = v0tmp.template negTrack_as(); // negative daughter - if (!checkV0leg(pos_tmp) || !checkV0leg(ele_tmp)) { - continue; - } - - if (v0.index() == v0tmp.index()) { // don't check onviously, exactly the same v0. - // LOGF(info, "don't check the exactly the same 2 V0s"); - continue; - } - if ((ele.globalIndex() == ele_tmp.globalIndex() || pos.globalIndex() == pos_tmp.globalIndex()) && v0.dcaV0daughters() > v0tmp.dcaV0daughters()) { - // LOGF(info, "!reject! | collision id = %d | g1 id = %d , g2 id = %d , posid1 = %d , eleid1 = %d , posid2 = %d , eleid2 = %d , pca1 = %f , pca2 = %f", - // collision.globalIndex(), v0.index(), v0tmp.index(), pos.globalIndex(), ele.globalIndex(), pos_tmp.globalIndex(), ele_tmp.globalIndex(), v0.dcaV0daughters(), v0tmp.dcaV0daughters()); - flag_closer = false; - break; - } - } // end of v0tmp loop - - if (!flag_closer) { - continue; - } - if (!ele.has_mcParticle() || !pos.has_mcParticle()) { continue; // If no MC particle is found, skip the v0 } eV0Confirmation v0Status = isTrueV0(v0, pos, ele); - fRegistry.get(HIST("hV0Confirmation"))->Fill(v0Status); + float xyz[3] = {0.f, 0.f, 0.f}; + Vtx_recalculation(o2::base::Propagator::Instance(), pos, ele, xyz); recalculatedVertexParameters recalculatedVertex; - Vtx_recalculation(pos, ele, &recalculatedVertex); + recalculatedVertex.recalculatedConversionPoint[0] = xyz[0]; + recalculatedVertex.recalculatedConversionPoint[1] = xyz[1]; + recalculatedVertex.recalculatedConversionPoint[2] = xyz[2]; v0photons(collision.globalIndex(), v0legs.lastIndex() + 1, v0legs.lastIndex() + 2, v0.x(), v0.y(), v0.z(), @@ -579,84 +509,6 @@ struct skimmerGammaConversions { } return kGoodMcMother; } - - template - void Vtx_recalculation(T lTrackPos, T lTrackNeg, recalculatedVertexParameters* recalculatedVertex) - { - o2::base::Propagator* prop = o2::base::Propagator::Instance(); // This singleton propagator requires some initialisation of the CCDB object. - float bz = prop->getNominalBz(); - - //******************************************************* - - // o2::track::TrackParametrizationWithError = TrackParCov, I use the full version to have control over the data type - o2::track::TrackParametrizationWithError trackPosInformation = getTrackParCov(lTrackPos); // first get an object that stores Track information (positive) - o2::track::TrackParametrizationWithError trackNegInformation = getTrackParCov(lTrackNeg); // first get an object that stores Track information (negative) - - o2::track::TrackAuxPar helixPos(trackPosInformation, bz); // This object is a descendant of a CircleXY and stores cirlce information with respect to the magnetic field. This object uses functions and information of the o2::track::TrackParametrizationWithError object (positive) - o2::track::TrackAuxPar helixNeg(trackNegInformation, bz); // This object is a descendant of a CircleXY and stores cirlce information with respect to the magnetic field. This object uses functions and information of the o2::track::TrackParametrizationWithError object (negative) - - recalculatedVertex->recalculatedConversionPoint[0] = (helixPos.xC * helixNeg.rC + helixNeg.xC * helixPos.rC) / (helixPos.rC + helixNeg.rC); // This calculates the coordinates of the conversion point as an weighted average of the two helix centers. xC and yC should be the global coordinates for the helix center as far as I understand. But you can double check the code of trackPosInformation.getCircleParamsLoc - recalculatedVertex->recalculatedConversionPoint[1] = (helixPos.yC * helixNeg.rC + helixNeg.yC * helixPos.rC) / (helixPos.rC + helixNeg.rC); // If this calculation doesn't work check if the rotateZ function, because the "documentation" says I get global coordinates but maybe i don't. - - // I am unsure about the Z calculation but this is how it is done in AliPhysics as far as I understand - o2::track::TrackParametrizationWithError trackPosInformationCopy = o2::track::TrackParametrizationWithError(trackPosInformation); - o2::track::TrackParametrizationWithError trackNegInformationCopy = o2::track::TrackParametrizationWithError(trackNegInformation); - - // I think this calculation gets the closest point on the track to the conversion point - // This alpha is a different alpha than the usual alpha and I think it is the angle between X axis and conversion point - Double_t alphaPos = TMath::Pi() + TMath::ATan2(-(recalculatedVertex->recalculatedConversionPoint[1] - helixPos.yC), (recalculatedVertex->recalculatedConversionPoint[0] - helixPos.xC)); - Double_t alphaNeg = TMath::Pi() + TMath::ATan2(-(recalculatedVertex->recalculatedConversionPoint[1] - helixNeg.yC), (recalculatedVertex->recalculatedConversionPoint[0] - helixNeg.xC)); - - Double_t vertexXPos = helixPos.xC + helixPos.rC * TMath::Cos(alphaPos); - Double_t vertexYPos = helixPos.yC + helixPos.rC * TMath::Sin(alphaPos); - Double_t vertexXNeg = helixNeg.xC + helixNeg.rC * TMath::Cos(alphaNeg); - Double_t vertexYNeg = helixNeg.yC + helixNeg.rC * TMath::Sin(alphaNeg); - - TVector2 vertexPos(vertexXPos, vertexYPos); - TVector2 vertexNeg(vertexXNeg, vertexYNeg); - - // Convert to local coordinate system - TVector2 vertexPosRot = vertexPos.Rotate(-trackPosInformationCopy.getAlpha()); - TVector2 vertexNegRot = vertexNeg.Rotate(-trackNegInformationCopy.getAlpha()); - - prop->propagateToX(trackPosInformationCopy, - vertexPosRot.X(), - bz, - o2::base::PropagatorImpl::MAX_SIN_PHI, - o2::base::PropagatorImpl::MAX_STEP, - o2::base::PropagatorImpl::MatCorrType::USEMatCorrNONE); - // o2::base::PropagatorImpl::MatCorrType::USEMatCorrLUT); - prop->propagateToX(trackNegInformationCopy, - vertexNegRot.X(), - bz, - o2::base::PropagatorImpl::MAX_SIN_PHI, - o2::base::PropagatorImpl::MAX_STEP, - o2::base::PropagatorImpl::MatCorrType::USEMatCorrNONE); - // o2::base::PropagatorImpl::MatCorrType::USEMatCorrLUT); - - // TODO: This is still off and needs to be checked... - recalculatedVertex->recalculatedConversionPoint[2] = (trackPosInformationCopy.getZ() * helixNeg.rC + trackNegInformationCopy.getZ() * helixPos.rC) / (helixPos.rC + helixNeg.rC); - - KFPTrack kFTrackPos = createKFPTrackFromTrackParCov(trackPosInformationCopy, lTrackPos.sign(), lTrackPos.tpcNClsFound(), lTrackPos.tpcChi2NCl()); - int pdg_ePlus = -11; // e+ - KFParticle kFParticleEPlus(kFTrackPos, pdg_ePlus); - - KFPTrack kFTrackNeg = createKFPTrackFromTrackParCov(trackNegInformationCopy, lTrackNeg.sign(), lTrackNeg.tpcNClsFound(), lTrackNeg.tpcChi2NCl()); - int pdg_eMinus = 11; // e- - KFParticle kFParticleEMinus(kFTrackNeg, pdg_eMinus); - - KFParticle gammaKF; - gammaKF.SetConstructMethod(2); - gammaKF.AddDaughter(kFParticleEPlus); - gammaKF.AddDaughter(kFParticleEMinus); - gammaKF.SetNonlinearMassConstraint(kfMassConstrain); - - if (gammaKF.GetNDF() == 0) { - recalculatedVertex->KFParticleChi2DividedByNDF = -1.f; - } else { - recalculatedVertex->KFParticleChi2DividedByNDF = gammaKF.GetChi2() / gammaKF.GetNDF(); - } - } }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGEM/PhotonMeson/Tasks/pcmQC.cxx b/PWGEM/PhotonMeson/Tasks/pcmQC.cxx index 6a8c37ae866..cc0730b7094 100644 --- a/PWGEM/PhotonMeson/Tasks/pcmQC.cxx +++ b/PWGEM/PhotonMeson/Tasks/pcmQC.cxx @@ -163,6 +163,7 @@ struct PCMQC { if (cut.IsSelected(v0)) { o2::aod::emphotonhistograms::FillHistClass(list_v0_cut, "", v0); nv0++; + reinterpret_cast(fMainList->FindObject("V0")->FindObject(cut.GetName())->FindObject("hV0R_minTrackX"))->Fill(v0.recalculatedVtxR(), std::min(pos.x(), ele.x())); for (auto& leg : {pos, ele}) { o2::aod::emphotonhistograms::FillHistClass(list_v0leg_cut, "", leg); } diff --git a/PWGEM/PhotonMeson/Tasks/pcmQCMC.cxx b/PWGEM/PhotonMeson/Tasks/pcmQCMC.cxx index 612375efdb8..60602dc3432 100644 --- a/PWGEM/PhotonMeson/Tasks/pcmQCMC.cxx +++ b/PWGEM/PhotonMeson/Tasks/pcmQCMC.cxx @@ -216,6 +216,7 @@ struct PCMQCMC { reinterpret_cast(fMainList->FindObject("V0")->FindObject(cut.GetName())->FindObject("hConvPoint_diffZ_recalc"))->Fill(elemc.vz(), v0.recalculatedVtxZ() - elemc.vz()); nv0++; + reinterpret_cast(fMainList->FindObject("V0")->FindObject(cut.GetName())->FindObject("hV0R_minTrackX"))->Fill(v0.recalculatedVtxR(), std::min(pos.x(), ele.x())); for (auto& leg : {pos, ele}) { o2::aod::emphotonhistograms::FillHistClass(list_v0leg_cut, "", leg); } diff --git a/PWGEM/PhotonMeson/Utils/PCMUtilities.h b/PWGEM/PhotonMeson/Utils/PCMUtilities.h index 3c027f821e4..db2314c34fc 100644 --- a/PWGEM/PhotonMeson/Utils/PCMUtilities.h +++ b/PWGEM/PhotonMeson/Utils/PCMUtilities.h @@ -15,13 +15,15 @@ #ifndef PWGEM_PHOTONMESON_UTILS_PCMUTILITIES_H_ #define PWGEM_PHOTONMESON_UTILS_PCMUTILITIES_H_ +#include +#include "DCAFitter/HelixHelper.h" +#include "DetectorsBase/Propagator.h" +#include "Common/Core/trackUtilities.h" #include "Framework/AnalysisTask.h" //_______________________________________________________________________ -bool checkAP(float alpha, float qt) +bool checkAP(const float alpha, const float qt, const float alpha_max = 0.95, const float qt_max = 0.05) { - const float alpha_max = 0.95; - const float qt_max = 0.05; float ellipse = pow(alpha / alpha_max, 2) + pow(qt / qt_max, 2); if (ellipse < 1.0) { return true; @@ -30,6 +32,62 @@ bool checkAP(float alpha, float qt) } } //_______________________________________________________________________ +template +void Vtx_recalculation(o2::base::Propagator* prop, T lTrackPos, T lTrackNeg, float xyz[3]) +{ + float bz = prop->getNominalBz(); + + //******************************************************* + + // o2::track::TrackParametrizationWithError = TrackParCov, I use the full version to have control over the data type + o2::track::TrackParametrizationWithError trackPosInformation = getTrackParCov(lTrackPos); // first get an object that stores Track information (positive) + o2::track::TrackParametrizationWithError trackNegInformation = getTrackParCov(lTrackNeg); // first get an object that stores Track information (negative) + + o2::track::TrackAuxPar helixPos(trackPosInformation, bz); // This object is a descendant of a CircleXY and stores cirlce information with respect to the magnetic field. This object uses functions and information of the o2::track::TrackParametrizationWithError object (positive) + o2::track::TrackAuxPar helixNeg(trackNegInformation, bz); // This object is a descendant of a CircleXY and stores cirlce information with respect to the magnetic field. This object uses functions and information of the o2::track::TrackParametrizationWithError object (negative) + + xyz[0] = (helixPos.xC * helixNeg.rC + helixNeg.xC * helixPos.rC) / (helixPos.rC + helixNeg.rC); // This calculates the coordinates of the conversion point as an weighted average of the two helix centers. xC and yC should be the global coordinates for the helix center as far as I understand. But you can double check the code of trackPosInformation.getCircleParamsLoc + xyz[1] = (helixPos.yC * helixNeg.rC + helixNeg.yC * helixPos.rC) / (helixPos.rC + helixNeg.rC); // If this calculation doesn't work check if the rotateZ function, because the "documentation" says I get global coordinates but maybe i don't. + + // I am unsure about the Z calculation but this is how it is done in AliPhysics as far as I understand + o2::track::TrackParametrizationWithError trackPosInformationCopy = o2::track::TrackParametrizationWithError(trackPosInformation); + o2::track::TrackParametrizationWithError trackNegInformationCopy = o2::track::TrackParametrizationWithError(trackNegInformation); + + // I think this calculation gets the closest point on the track to the conversion point + // This alpha is a different alpha than the usual alpha and I think it is the angle between X axis and conversion point + Double_t alphaPos = TMath::Pi() + TMath::ATan2(-(xyz[1] - helixPos.yC), (xyz[0] - helixPos.xC)); + Double_t alphaNeg = TMath::Pi() + TMath::ATan2(-(xyz[1] - helixNeg.yC), (xyz[0] - helixNeg.xC)); + + Double_t vertexXPos = helixPos.xC + helixPos.rC * TMath::Cos(alphaPos); + Double_t vertexYPos = helixPos.yC + helixPos.rC * TMath::Sin(alphaPos); + Double_t vertexXNeg = helixNeg.xC + helixNeg.rC * TMath::Cos(alphaNeg); + Double_t vertexYNeg = helixNeg.yC + helixNeg.rC * TMath::Sin(alphaNeg); + + TVector2 vertexPos(vertexXPos, vertexYPos); + TVector2 vertexNeg(vertexXNeg, vertexYNeg); + + // Convert to local coordinate system + TVector2 vertexPosRot = vertexPos.Rotate(-trackPosInformationCopy.getAlpha()); + TVector2 vertexNegRot = vertexNeg.Rotate(-trackNegInformationCopy.getAlpha()); + + prop->propagateToX(trackPosInformationCopy, + vertexPosRot.X(), + bz, + o2::base::PropagatorImpl::MAX_SIN_PHI, + o2::base::PropagatorImpl::MAX_STEP, + o2::base::PropagatorImpl::MatCorrType::USEMatCorrNONE); + // o2::base::PropagatorImpl::MatCorrType::USEMatCorrLUT); + prop->propagateToX(trackNegInformationCopy, + vertexNegRot.X(), + bz, + o2::base::PropagatorImpl::MAX_SIN_PHI, + o2::base::PropagatorImpl::MAX_STEP, + o2::base::PropagatorImpl::MatCorrType::USEMatCorrNONE); + // o2::base::PropagatorImpl::MatCorrType::USEMatCorrLUT); + + // TODO: This is still off and needs to be checked... + xyz[2] = (trackPosInformationCopy.getZ() * helixNeg.rC + trackNegInformationCopy.getZ() * helixPos.rC) / (helixPos.rC + helixNeg.rC); +} //_______________________________________________________________________ //_______________________________________________________________________ #endif // PWGEM_PHOTONMESON_UTILS_PCMUTILITIES_H_