@@ -3202,6 +3202,125 @@ void VarManager::FillPairME(T1 const& t1, T2 const& t2, float* values)
32023202 values[kDeltaEtaPair2 ] = v1.Eta () - v2.Eta ();
32033203 }
32043204
3205+ // polarization parameters
3206+ bool useHE = fgUsedVars[kCosThetaHE ] || fgUsedVars[kPhiHE ]; // helicity frame
3207+ bool useCS = fgUsedVars[kCosThetaCS ] || fgUsedVars[kPhiCS ]; // Collins-Soper frame
3208+ bool usePP = fgUsedVars[kCosThetaPP ]; // production plane frame
3209+ bool useRM = fgUsedVars[kCosThetaRM ]; // Random frame
3210+
3211+ if (useHE || useCS || usePP || useRM) {
3212+ ROOT::Math::Boost boostv12{v12.BoostToCM ()};
3213+ ROOT::Math::XYZVectorF v1_CM{(boostv12 (v1).Vect ()).Unit ()};
3214+ ROOT::Math::XYZVectorF v2_CM{(boostv12 (v2).Vect ()).Unit ()};
3215+ ROOT::Math::XYZVectorF Beam1_CM{(boostv12 (fgBeamA).Vect ()).Unit ()};
3216+ ROOT::Math::XYZVectorF Beam2_CM{(boostv12 (fgBeamC).Vect ()).Unit ()};
3217+
3218+ // using positive sign convention for the first track
3219+ ROOT::Math::XYZVectorF v_CM = (t1.sign () > 0 ? v1_CM : v2_CM);
3220+
3221+ if (useHE) {
3222+ ROOT::Math::XYZVectorF zaxis_HE{(v12.Vect ()).Unit ()};
3223+ ROOT::Math::XYZVectorF yaxis_HE{(Beam1_CM.Cross (Beam2_CM)).Unit ()};
3224+ ROOT::Math::XYZVectorF xaxis_HE{(yaxis_HE.Cross (zaxis_HE)).Unit ()};
3225+ if (fgUsedVars[kCosThetaHE ])
3226+ values[kCosThetaHE ] = zaxis_HE.Dot (v_CM);
3227+ if (fgUsedVars[kPhiHE ]) {
3228+ values[kPhiHE ] = TMath::ATan2 (yaxis_HE.Dot (v_CM), xaxis_HE.Dot (v_CM));
3229+ if (values[kPhiHE ] < 0 ) {
3230+ values[kPhiHE ] += 2 * TMath::Pi (); // ensure phi is in [0, 2pi]
3231+ }
3232+ }
3233+ if (fgUsedVars[kPhiTildeHE ]) {
3234+ if (fgUsedVars[kCosThetaHE ] && fgUsedVars[kPhiHE ]) {
3235+ if (values[kCosThetaHE ] > 0 ) {
3236+ values[kPhiTildeHE ] = values[kPhiHE ] - 0.25 * TMath::Pi (); // phi_tilde = phi - pi/4
3237+ if (values[kPhiTildeHE ] < 0 ) {
3238+ values[kPhiTildeHE ] += 2 * TMath::Pi (); // ensure phi_tilde is in [0, 2pi]
3239+ }
3240+ } else {
3241+ values[kPhiTildeHE ] = values[kPhiHE ] - 0.75 * TMath::Pi (); // phi_tilde = phi - 3pi/4
3242+ if (values[kPhiTildeHE ] < 0 ) {
3243+ values[kPhiTildeHE ] += 2 * TMath::Pi (); // ensure phi_tilde is in [0, 2pi]
3244+ }
3245+ }
3246+ } else {
3247+ values[kPhiTildeHE ] = -999 ; // not computable
3248+ }
3249+ }
3250+ }
3251+
3252+ if (useCS) {
3253+ ROOT::Math::XYZVectorF zaxis_CS{(Beam1_CM - Beam2_CM).Unit ()};
3254+ ROOT::Math::XYZVectorF yaxis_CS{(Beam1_CM.Cross (Beam2_CM)).Unit ()};
3255+ ROOT::Math::XYZVectorF xaxis_CS{(yaxis_CS.Cross (zaxis_CS)).Unit ()};
3256+ if (fgUsedVars[kCosThetaCS ])
3257+ values[kCosThetaCS ] = zaxis_CS.Dot (v_CM);
3258+ if (fgUsedVars[kPhiCS ]) {
3259+ values[kPhiCS ] = TMath::ATan2 (yaxis_CS.Dot (v_CM), xaxis_CS.Dot (v_CM));
3260+ if (values[kPhiCS ] < 0 ) {
3261+ values[kPhiCS ] += 2 * TMath::Pi (); // ensure phi is in [0, 2pi]
3262+ }
3263+ }
3264+ if (fgUsedVars[kPhiTildeCS ]) {
3265+ if (fgUsedVars[kCosThetaCS ] && fgUsedVars[kPhiCS ]) {
3266+ if (values[kCosThetaCS ] > 0 ) {
3267+ values[kPhiTildeCS ] = values[kPhiCS ] - 0.25 * TMath::Pi (); // phi_tilde = phi - pi/4
3268+ if (values[kPhiTildeCS ] < 0 ) {
3269+ values[kPhiTildeCS ] += 2 * TMath::Pi (); // ensure phi_tilde is in [0, 2pi]
3270+ }
3271+ } else {
3272+ values[kPhiTildeCS ] = values[kPhiCS ] - 0.75 * TMath::Pi (); // phi_tilde = phi - 3pi/4
3273+ if (values[kPhiTildeCS ] < 0 ) {
3274+ values[kPhiTildeCS ] += 2 * TMath::Pi (); // ensure phi_tilde is in [0, 2pi]
3275+ }
3276+ }
3277+ } else {
3278+ values[kPhiTildeCS ] = -999 ; // not computable
3279+ }
3280+ }
3281+ }
3282+
3283+ if (usePP) {
3284+ ROOT::Math::XYZVector zaxis_PP = ROOT::Math::XYZVector (v12.Py (), -v12.Px (), 0 .f );
3285+ ROOT::Math::XYZVector yaxis_PP{(v12.Vect ()).Unit ()};
3286+ ROOT::Math::XYZVector xaxis_PP{(yaxis_PP.Cross (zaxis_PP)).Unit ()};
3287+ if (fgUsedVars[kCosThetaPP ]) {
3288+ values[kCosThetaPP ] = zaxis_PP.Dot (v_CM) / std::sqrt (zaxis_PP.Mag2 ());
3289+ }
3290+ if (fgUsedVars[kPhiPP ]) {
3291+ values[kPhiPP ] = TMath::ATan2 (yaxis_PP.Dot (v_CM), xaxis_PP.Dot (v_CM));
3292+ if (values[kPhiPP ] < 0 ) {
3293+ values[kPhiPP ] += 2 * TMath::Pi (); // ensure phi is in [0, 2pi]
3294+ }
3295+ }
3296+ if (fgUsedVars[kPhiTildePP ]) {
3297+ if (fgUsedVars[kCosThetaPP ] && fgUsedVars[kPhiPP ]) {
3298+ if (values[kCosThetaPP ] > 0 ) {
3299+ values[kPhiTildePP ] = values[kPhiPP ] - 0.25 * TMath::Pi (); // phi_tilde = phi - pi/4
3300+ if (values[kPhiTildePP ] < 0 ) {
3301+ values[kPhiTildePP ] += 2 * TMath::Pi (); // ensure phi_tilde is in [0, 2pi]
3302+ }
3303+ } else {
3304+ values[kPhiTildePP ] = values[kPhiPP ] - 0.75 * TMath::Pi (); // phi_tilde = phi - 3pi/4
3305+ if (values[kPhiTildePP ] < 0 ) {
3306+ values[kPhiTildePP ] += 2 * TMath::Pi (); // ensure phi_tilde is in [0, 2pi]
3307+ }
3308+ }
3309+ } else {
3310+ values[kPhiTildePP ] = -999 ; // not computable
3311+ }
3312+ }
3313+ }
3314+
3315+ if (useRM) {
3316+ double randomCostheta = gRandom ->Uniform (-1 ., 1 .);
3317+ double randomPhi = gRandom ->Uniform (0 ., 2 . * TMath::Pi ());
3318+ ROOT::Math::XYZVectorF zaxis_RM (randomCostheta, std::sqrt (1 - randomCostheta * randomCostheta) * std::cos (randomPhi), std::sqrt (1 - randomCostheta * randomCostheta) * std::sin (randomPhi));
3319+ if (fgUsedVars[kCosThetaRM ])
3320+ values[kCosThetaRM ] = zaxis_RM.Dot (v_CM);
3321+ }
3322+ }
3323+
32053324 if constexpr ((fillMap & ReducedEventQvector) > 0 || (fillMap & CollisionQvect) > 0 ) {
32063325 // TODO: provide different computations for vn
32073326 // Compute the scalar product UQ for two muon from different event using Q-vector from A, for second and third harmonic
0 commit comments