Skip to content

Commit e6a853e

Browse files
miedema-11jokonig
miedema-11
authored and
jokonig
committed
[PWGUD] fix ptdiff (AliceO2Group#17274)
resolve merge conflicts
1 parent 4d0c041 commit e6a853e

3 files changed

Lines changed: 213 additions & 31 deletions

File tree

‎PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx‎

Lines changed: 111 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -114,6 +114,7 @@ struct PhotonConversionBuilder {
114114
// Operation and minimisation criteria
115115
Configurable<double> d_bz_input{"d_bz", -999, "bz field, -999 is automatic"};
116116
Configurable<int> useMatCorrType{"useMatCorrType", 0, "0: none, 1: TGeo, 2: LUT"};
117+
Configurable<int> modeTrackPropagation{"modeTrackPropagation", 0, "0: use real track propagation, including material, 1: use fast approximation using only geometry, 2: Use real track propagation and make comparison to fast propagation (only for debugging and testing)"};
117118

118119
// single track cuts
119120
Configurable<int> min_ncluster_tpc{"min_ncluster_tpc", 0, "min ncluster tpc"};
@@ -325,6 +326,14 @@ struct PhotonConversionBuilder {
325326
registry.add("V0/hBDTScoreAfterCutVsPt", "BDT score after cut vs pT; pT (GeV/c); BDT score", {HistType::kTH2F, {{1000, 0.0f, 20.0f}, {1000, 0.0f, 1.0f}}});
326327
}
327328
}
329+
330+
// Compare proper propagation and fast geometrical propagation
331+
if (modeTrackPropagation == 2) {
332+
registry.add("V0Leg/hDCAxyPropagationCompare", "Comparison of DCA_{xy} propagation;DCA_{xy} (cm) proper propagation; DCA_{xy} (cm) geom. propagation", {HistType::kTH2F, {{200, -10., 10.}, {200, -10., 10.}}});
333+
registry.add("V0Leg/hDCAzPropagationCompare", "Comparison of DCA_{z} propagation;DCA_{z} (cm) proper propagation; DCA_{z} (cm) geom. propagation", {HistType::kTH2F, {{200, -10., 10.}, {200, -10., 10.}}});
334+
registry.add("V0/hPhivPropagationCompare", "Comparison of #phi_{v};#phi_{v} proper propagation; #phi_{v} geom. propagation", {HistType::kTH2F, {{100, 0., 3.2}, {100, 0., 3.2}}});
335+
registry.add("V0/hPsiPairPropagationCompare", "Comparison of #Psi_{pair};#Psi_{pair} proper propagation; #Psi_{pair} geom. propagation", {HistType::kTH2F, {{100, 0., 1.6}, {100, 0., 1.6}}});
336+
}
328337
}
329338

330339
void initCCDB(aod::BCsWithTimestamps::iterator const& bc)
@@ -555,8 +564,23 @@ struct PhotonConversionBuilder {
555564
return;
556565
}
557566
auto pTrackC = pTrack;
567+
o2::math_utils::Point3D<float> vtxPrim{
568+
collision.posX(),
569+
collision.posY(),
570+
collision.posZ()};
571+
if (modeTrackPropagation > 0) {
572+
dcaInfo = CalculateDCAFast(pTrackC, vtxPrim, d_bz);
573+
}
558574
pTrackC.setPID(o2::track::PID::Electron);
559-
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, pTrackC, 2.f, matCorr, &dcaInfo);
575+
576+
std::array<float, 2> dcaInfoFast = dcaInfo;
577+
if (modeTrackPropagation != 1) {
578+
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, pTrackC, 2.f, matCorr, &dcaInfo);
579+
if (modeTrackPropagation == 2) {
580+
registry.fill(HIST("V0Leg/hDCAxyPropagationCompare"), dcaInfo[0], dcaInfoFast[0]);
581+
registry.fill(HIST("V0Leg/hDCAzPropagationCompare"), dcaInfo[1], dcaInfoFast[1]);
582+
}
583+
}
560584
auto posdcaXY = dcaInfo[0];
561585
auto posdcaZ = dcaInfo[1];
562586

@@ -566,8 +590,19 @@ struct PhotonConversionBuilder {
566590
return;
567591
}
568592
auto nTrackC = nTrack;
593+
if (modeTrackPropagation > 0) {
594+
dcaInfo = CalculateDCAFast(nTrackC, vtxPrim, d_bz);
595+
}
569596
nTrackC.setPID(o2::track::PID::Electron);
570-
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, nTrackC, 2.f, matCorr, &dcaInfo);
597+
598+
dcaInfoFast = dcaInfo;
599+
if (modeTrackPropagation != 1) {
600+
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, nTrackC, 2.f, matCorr, &dcaInfo);
601+
if (modeTrackPropagation == 2) {
602+
registry.fill(HIST("V0Leg/hDCAxyPropagationCompare"), dcaInfo[0], dcaInfoFast[0]);
603+
registry.fill(HIST("V0Leg/hDCAzPropagationCompare"), dcaInfo[1], dcaInfoFast[1]);
604+
}
605+
}
571606
auto eledcaXY = dcaInfo[0];
572607
auto eledcaZ = dcaInfo[1];
573608

@@ -587,30 +622,76 @@ struct PhotonConversionBuilder {
587622

588623
float phiv = 999.f;
589624
float psipair = 999.f;
625+
float phivFast = 999.f;
626+
float psipairFast = 999.f;
590627
float baseR = std::hypot(xyz[0], xyz[1]);
591-
std::array<float, 3> offsetsR = {propV0LegsRadius, 30.f, 10.f};
592-
bool pPropagatedSuccess = false;
593-
bool nPropagatedSuccess = false;
594-
auto pTrackProp = pTrack;
595-
auto nTrackProp = nTrack;
596-
for (const float& offsetR : offsetsR) {
597-
pTrackProp = pTrack;
598-
pTrackProp.setPID(o2::track::PID::Electron);
599-
nTrackProp = nTrack;
600-
nTrackProp.setPID(o2::track::PID::Electron);
601-
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, pTrackProp, 2.f, matCorr, &dcaInfo);
602-
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, nTrackProp, 2.f, matCorr, &dcaInfo);
603-
pPropagatedSuccess = o2::base::Propagator::Instance()->propagateToR(pTrackProp, baseR + offsetR);
604-
nPropagatedSuccess = o2::base::Propagator::Instance()->propagateToR(nTrackProp, baseR + offsetR);
605-
if (pPropagatedSuccess && nPropagatedSuccess) {
606-
KFPTrack kfp_track_posProp = createKFPTrackFromTrackParCov(pTrackProp, pos.sign(), pos.tpcNClsFound(), pos.tpcChi2NCl());
607-
KFPTrack kfp_track_eleProp = createKFPTrackFromTrackParCov(nTrackProp, ele.sign(), ele.tpcNClsFound(), ele.tpcChi2NCl());
608-
phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(kfp_track_posProp.GetPx(), kfp_track_posProp.GetPy(), kfp_track_posProp.GetPz(), kfp_track_eleProp.GetPx(), kfp_track_eleProp.GetPy(), kfp_track_eleProp.GetPz(), pos.sign(), ele.sign(), d_bz);
609-
psipair = o2::aod::pwgem::dilepton::utils::pairutil::getPsiPair(kfp_track_posProp.GetPx(), kfp_track_posProp.GetPy(), kfp_track_posProp.GetPz(), kfp_track_eleProp.GetPx(), kfp_track_eleProp.GetPy(), kfp_track_eleProp.GetPz());
610-
break;
628+
// This method uses the track Helix instead of the full propagation.
629+
// Hence, it is only an approximation but much faster
630+
if (modeTrackPropagation > 0) {
631+
632+
o2::track::TrackAuxPar helixPosEle(nTrack, d_bz);
633+
o2::track::TrackAuxPar helixPosPos(pTrack, d_bz);
634+
635+
float diffX = helixPosEle.xC - helixPosPos.xC;
636+
float diffY = helixPosEle.yC - helixPosPos.yC;
637+
float phiHelix = std::atan2(diffY, diffX);
638+
phiHelix -= o2::constants::math::PI / 2.;
639+
if (phiHelix > 2 * o2::constants::math::PI) {
640+
phiHelix -= o2::constants::math::PI * 2;
641+
}
642+
if (phiHelix < 0) {
643+
phiHelix += o2::constants::math::PI * 2;
644+
}
645+
646+
// Electron
647+
float arcLenghtEle = helixPosEle.rC * 0.9 > propV0LegsRadius ? std::asin(propV0LegsRadius / helixPosEle.rC) * helixPosEle.rC : o2::constants::math::PI / 2.2 * helixPosEle.rC; // This assumes that the photon momentum vector is a tangent of the circle
648+
auto propTrackEle = getPropMomentumFromTrackHelix(arcLenghtEle, ele.pt(), d_bz / 10., ele.phi(), ele.tgl(), phiHelix - ele.phi());
649+
// Positron
650+
float arcLenghtPos = helixPosPos.rC * 0.9 > propV0LegsRadius ? std::asin(propV0LegsRadius / helixPosPos.rC) * helixPosPos.rC : o2::constants::math::PI / 2.2 * helixPosPos.rC; // This assumes that the photon momentum vector is a tangent of the circle
651+
auto propTrackPos = getPropMomentumFromTrackHelix(arcLenghtPos, pos.pt(), d_bz / 10., pos.phi(), pos.tgl(), phiHelix - pos.phi());
652+
653+
phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(propTrackPos[0], propTrackPos[1], propTrackPos[2], propTrackEle[0], propTrackEle[1], propTrackEle[2], pos.sign(), ele.sign(), d_bz);
654+
psipair = o2::aod::pwgem::dilepton::utils::pairutil::getPsiPair(propTrackPos[0], propTrackPos[1], propTrackPos[2], propTrackEle[0], propTrackEle[1], propTrackEle[2]);
655+
// Store values for later comparison
656+
phivFast = phiv;
657+
psipairFast = psipair;
658+
}
659+
// This uses the full propagation including material effects
660+
if (modeTrackPropagation != 1) {
661+
float offsetsR[3] = {propV0LegsRadius, 30.f, 10.f};
662+
bool pPropagatedSuccess = false;
663+
bool nPropagatedSuccess = false;
664+
auto pTrackProp = pTrack;
665+
auto nTrackProp = nTrack;
666+
for (float offsetR : offsetsR) {
667+
pTrackProp = pTrack;
668+
pTrackProp.setPID(o2::track::PID::Electron);
669+
nTrackProp = nTrack;
670+
nTrackProp.setPID(o2::track::PID::Electron);
671+
672+
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, pTrackProp, 2.f, matCorr, &dcaInfo);
673+
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, nTrackProp, 2.f, matCorr, &dcaInfo);
674+
675+
pPropagatedSuccess = o2::base::Propagator::Instance()->propagateToR(pTrackProp, baseR + offsetR);
676+
nPropagatedSuccess = o2::base::Propagator::Instance()->propagateToR(nTrackProp, baseR + offsetR);
677+
678+
if (pPropagatedSuccess && nPropagatedSuccess) {
679+
KFPTrack kfp_track_posProp = createKFPTrackFromTrackParCov(pTrackProp, pos.sign(), pos.tpcNClsFound(), pos.tpcChi2NCl());
680+
KFPTrack kfp_track_eleProp = createKFPTrackFromTrackParCov(nTrackProp, ele.sign(), ele.tpcNClsFound(), ele.tpcChi2NCl());
681+
phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(kfp_track_posProp.GetPx(), kfp_track_posProp.GetPy(), kfp_track_posProp.GetPz(), kfp_track_eleProp.GetPx(), kfp_track_eleProp.GetPy(), kfp_track_eleProp.GetPz(), pos.sign(), ele.sign(), d_bz);
682+
psipair = o2::aod::pwgem::dilepton::utils::pairutil::getPsiPair(kfp_track_posProp.GetPx(), kfp_track_posProp.GetPy(), kfp_track_posProp.GetPz(), kfp_track_eleProp.GetPx(), kfp_track_eleProp.GetPy(), kfp_track_eleProp.GetPz());
683+
break;
684+
} else {
685+
LOG(debug) << "Propagation to offset" << offsetR << " cm failed for " << (pPropagatedSuccess ? "negative" : "positive") << " track. Trying smaller offset.";
686+
}
687+
}
688+
if (modeTrackPropagation == 2) {
689+
registry.fill(HIST("V0/hPhivPropagationCompare"), phiv, phivFast);
690+
registry.fill(HIST("V0/hPsiPairPropagationCompare"), psipair, psipairFast);
611691
}
612692
LOG(debug) << "Propagation to offset" << offsetR << " cm failed for " << (pPropagatedSuccess ? "negative" : "positive") << " track. Trying smaller offset.";
613693
}
694+
614695
if (phiv == 999.f || psipair == 999.f) {
615696
LOG(debug) << "Propagation failed for all radii (" << propV0LegsRadius << ", 30, 10 cm). Using default values for phiv and psipair (999.f).";
616697
}
@@ -985,14 +1066,14 @@ struct PhotonConversionBuilder {
9851066
fillV0Table<isMC, TBCs, TCollisions, TTracks>(v0, true);
9861067
} // end of fullv0Id loop
9871068

988-
for (const auto& collision : collisions) {
989-
if constexpr (isMC) {
990-
if (!collision.has_mcCollision()) {
991-
continue;
992-
}
993-
}
994-
// events_ngpcm(nv0_map[collision.globalIndex()]);
995-
} // end of collision loop
1069+
// for (const auto& collision : collisions) {
1070+
// if constexpr (isMC) {
1071+
// if (!collision.has_mcCollision()) {
1072+
// continue;
1073+
// }
1074+
// }
1075+
// // events_ngpcm(nv0_map[collision.globalIndex()]);
1076+
// } // end of collision loop
9961077

9971078
pca_map.clear();
9981079
cospa_map.clear();

‎PWGEM/PhotonMeson/Utils/PCMUtilities.h‎

Lines changed: 90 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -104,10 +104,100 @@ inline void Vtx_recalculationParCov(o2::base::Propagator* prop, const o2::track:
104104

105105
xyz[2] = (trackPosInformationCopy.getZ() * helixNeg.rC + trackNegInformationCopy.getZ() * helixPos.rC) / (helixPos.rC + helixNeg.rC);
106106
}
107+
108+
//_______________________________________________________________________
109+
/// \brief Calculate DCA for tracks using the track helix and the inclination angl tan(lambda)
110+
/// \param trk track parameterization to obtain helix parameters
111+
/// \param vtx primary vertex position
112+
/// \param magField magnetic field strenght of L3
113+
/// \return DCAxy, DCAz
114+
template <typename TrackPrecision = float>
115+
std::array<float, 2> CalculateDCAFast(const o2::track::TrackParametrizationWithError<TrackPrecision>& trk, const o2::math_utils::Point3D<float>& vtx, const float magField)
116+
{
117+
118+
std::array<float, 2> dca;
119+
120+
// obtain circle from track in x-y plane
121+
const o2::track::TrackAuxPar helixPos(trk, magField);
122+
123+
// obtain position in global coordinates and tan(lambda)
124+
const auto posTrack = trk.getXYZGlo();
125+
const float trX = posTrack.X();
126+
const float trY = posTrack.Y();
127+
const float trZ = posTrack.Z();
128+
const float tangentLambda = trk.getTgl();
129+
130+
// Calculate DCAxy
131+
// Use distance in x and y between circle center and vtx, afterwards subtract radius of circle
132+
const float distX = helixPos.xC - vtx.X();
133+
const float distY = helixPos.yC - vtx.Y();
134+
const float trackCircCenter = sqrt(distX * distX + distY * distY);
135+
dca[0] = helixPos.rC - trackCircCenter; // this has to be changed
136+
137+
// Calculate DCAz
138+
// First step: Calculate arc lenght of circle between current position and primary vertex in x-y
139+
const float theta0 = std::atan2(trY - helixPos.yC, trX - helixPos.xC);
140+
const float thetav = std::atan2(vtx.Y() - helixPos.yC, vtx.X() - helixPos.xC);
141+
142+
// Make sure angle is between -pi and pi
143+
float dtheta = thetav - theta0;
144+
while (dtheta > M_PI) {
145+
dtheta -= 2.0 * M_PI;
146+
}
147+
while (dtheta <= -M_PI) {
148+
dtheta += 2.0 * M_PI;
149+
}
150+
151+
// arc-lenght along helix
152+
const float arcLenght = std::fabs(helixPos.rC * dtheta);
153+
154+
// get global z-position at DCA
155+
const float ZPosGlo = trZ - arcLenght * tangentLambda;
156+
157+
// DCA calculated from
158+
dca[1] = ZPosGlo - vtx.Z();
159+
160+
return dca;
161+
}
162+
163+
//_______________________________________________________________________
164+
/// \brief Calculate the track momentum at a different place of the track Helix.
165+
/// \param s arc-lenght where track should be propagated to
166+
/// \param pt pt of the track
167+
/// \param bz magnetic field strenght of L3
168+
/// \param phi0 initial phi angle from track
169+
/// \param tgl tan(lambda) of the track
170+
/// \param addPhi optional additional rotation of the track in phi-direction
171+
/// \return track momentum vector at new position
172+
inline std::array<float, 3> getPropMomentumFromTrackHelix(float s, float pt, float bz, float phi0, float tgl, float addPhi = 0.f)
173+
{
174+
175+
// Calculate the change in phi considering the track radius and the track arc lenght s
176+
float dphi = 0.3f * bz * s / 100. / pt; // s in cm
177+
178+
// Calculate the phi at the secondary vertex
179+
float phi = phi0 + dphi + addPhi;
180+
if (phi > 2 * o2::constants::math::PI) {
181+
phi -= o2::constants::math::PI * 2;
182+
}
183+
if (phi < 0) {
184+
phi += o2::constants::math::PI * 2;
185+
}
186+
187+
// Calculate px,y,z at the new propagated vertex
188+
std::array<float, 3> trackP;
189+
trackP[0] = std::cos(phi) * pt;
190+
trackP[1] = std::sin(phi) * pt;
191+
trackP[2] = tgl * pt;
192+
193+
return trackP;
194+
}
195+
107196
//_______________________________________________________________________
108197
template <typename TrackPrecision = float, typename T1, typename T2>
109198
inline void Vtx_recalculation(o2::base::Propagator* prop, T1 lTrackPos, T2 lTrackNeg, float xyz[3], o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE)
110199
{
200+
111201
// o2::track::TrackParametrizationWithError<TrackPrecision> = TrackParCov, I use the full version to have control over the data type
112202
o2::track::TrackParametrizationWithError<TrackPrecision> trackPosInformation = getTrackParCov(lTrackPos); // first get an object that stores Track information (positive)
113203
o2::track::TrackParametrizationWithError<TrackPrecision> trackNegInformation = getTrackParCov(lTrackNeg); // first get an object that stores Track information (negative)

‎PWGUD/Tasks/flowCumulantsUpc.cxx‎

Lines changed: 12 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -304,9 +304,20 @@ struct FlowCumulantsUpc {
304304
oba->Add(new TNamed(Form("Ch10Gap24_pt_%i", i + 1), "Ch10Gap24_pTDiff"));
305305
std::vector<std::string> userDefineGFWCorr = cfgUserDefineGFWCorr;
306306
std::vector<std::string> userDefineGFWName = cfgUserDefineGFWName;
307+
if (userDefineGFWName.size() != userDefineGFWCorr.size()) {
308+
LOGF(fatal, "The GFWConfig names you provided are NOT matching with configurations. userDefineGFWName.size(): %d, userDefineGFWCorr.size(): %d", userDefineGFWName.size(), userDefineGFWCorr.size());
309+
}
310+
LOGF(info, "User adding FlowContainer Array:");
307311
if (!userDefineGFWCorr.empty() && !userDefineGFWName.empty()) {
308312
for (uint i = 0; i < userDefineGFWName.size(); i++) {
309-
oba->Add(new TNamed(userDefineGFWName.at(i).c_str(), userDefineGFWName.at(i).c_str()));
313+
if (userDefineGFWCorr.at(i).find("poi") != std::string::npos) {
314+
LOGF(info, "%d: pT-diff array %s", i, userDefineGFWName.at(i).c_str());
315+
for (auto iPt = 0; iPt < fPtAxis->GetNbins(); iPt++)
316+
oba->Add(new TNamed(Form("%s_pt_%i", userDefineGFWName.at(i).c_str(), iPt + 1), Form("%s_pTDiff", userDefineGFWName.at(i).c_str())));
317+
} else {
318+
LOGF(info, "%d: %s", i, userDefineGFWName.at(i).c_str());
319+
oba->Add(new TNamed(userDefineGFWName.at(i).c_str(), userDefineGFWName.at(i).c_str()));
320+
}
310321
}
311322
}
312323
fFC->SetName("FlowContainer");

0 commit comments

Comments
 (0)