Skip to content
Snippets Groups Projects

Add J/psi resolution; re-write RC with Dummy instead of rc tracks

Merged Add J/psi resolution; re-write RC with Dummy instead of rc tracks
ziyue_work_branch into master
Merged Ziyue Zhang requested to merge ziyue_work_branch into master
Compare
and
4 files
+ 146
− 35
Compare changes
  • Side-by-side
  • Inline
Files
4
benchmarks/dvmp/analysis/dvmp.h
+ 73
− 11
@@ -46,7 +46,63 @@ namespace util {
}
//for tracking
inline inv_quant calc_inv_quant_rec(const std::vector<ROOT::Math::PxPyPzMVector>& parts, const double pdg_mass){
// inline inv_quant calc_inv_quant_rec(const std::vector<ROOT::Math::PxPyPzMVector>& parts, const double pdg_mass){
// int first = -1;
// int second = -1;
// double best_mass = -1;
// // go through all particle combinatorics, calculate the invariant mass
// // for each combination, and remember which combination is the closest
// // to the desired pdg_mass
// for (int i = 0; i < parts.size(); ++i) {
// for (int j = i + 1; j < parts.size(); ++j) {
// const double new_mass{(parts[i] + parts[j]).mass()};
// if (fabs(new_mass - pdg_mass) < fabs(best_mass - pdg_mass)) {
// first = i;
// second = j;
// best_mass = new_mass;
// }
// }
// }
// if (first < 0 || parts.size() < 3 ){
// inv_quant quantities = {-999., -999., -999., -999.};
// return quantities;
// }
// ROOT::Math::PxPyPzMVector pair_4p(parts[first] + parts[second]);
// ROOT::Math::PxPyPzMVector e1, P;
// double e1_Energy = sqrt(10.*10. + get_pdg_mass("electron")*get_pdg_mass("electron"));
// double P_Energy = sqrt(100.*100. + get_pdg_mass("proton")*get_pdg_mass("proton"));
// e1.SetPxPyPzE(0., 0., -10., e1_Energy);
// P.SetPxPyPzE(0., 0., 100., P_Energy);
// int scatteredIdx = -1;
// float dp = 10.;
// for(int i = 0 ; i < parts.size(); i++){
// if(i==first || i==second) continue;
// ROOT::Math::PxPyPzMVector k_prime(parts[i]);
// float ptmp = sqrt(parts[i].px()*parts[i].px() + parts[i].py()*parts[i].py() + parts[i].pz()*parts[i].pz());
// if( (k_prime.px()) * (pair_4p.px()) + (k_prime.py()) * (pair_4p.py()) + (k_prime.pz()) * (pair_4p.pz()) > 0. || ptmp >= 10.) continue; //angle between jpsi and scattered electron < pi/2, 3-momentum mag < 10.
// if(dp > 10.- ptmp){ //if there are more than one candidate of scattered electron, choose the one with highest 3-momentum mag
// scatteredIdx = i;
// dp = 10. - ptmp;
// }
// }
// if(scatteredIdx ==-1){
// inv_quant quantities = {-999., -999., -999., -999.};
// return quantities;
// }
// ROOT::Math::PxPyPzMVector q(e1 - parts[scatteredIdx]);
// ROOT::Math::PxPyPzMVector Delta(q - pair_4p);
// double nu = q.Dot(P) / P.mass();
// double Q2 = - q.Dot(q);
// double t = Delta.Dot(Delta);
// inv_quant quantities = {nu, Q2, Q2/2./P.mass()/nu, t};
// //inv_quant quantities = {-999., -999., -999., -999.};
// return quantities;
// }
//for tracking
inline inv_quant calc_inv_quant_rec(const std::vector<ROOT::Math::PxPyPzMVector>& parts, const double pdg_mass, const double daughter_mass){
int first = -1;
int second = -1;
double best_mass = -1;
@@ -55,7 +111,9 @@ namespace util {
// for each combination, and remember which combination is the closest
// to the desired pdg_mass
for (int i = 0; i < parts.size(); ++i) {
if( fabs(parts[i].mass() - daughter_mass)/daughter_mass > 0.001) continue;
for (int j = i + 1; j < parts.size(); ++j) {
if( fabs(parts[j].mass() - daughter_mass)/daughter_mass > 0.001) continue;
const double new_mass{(parts[i] + parts[j]).mass()};
if (fabs(new_mass - pdg_mass) < fabs(best_mass - pdg_mass)) {
first = i;
@@ -64,10 +122,14 @@ namespace util {
}
}
}
if (first < 0 || parts.size() < 3 ){
inv_quant quantities = {-999., -999., -999., -999.};
//find l-pair for jpsi
if (first < 0 || parts.size() < 3 ){ //fewer than 3 candidates
inv_quant quantities = {-99999., -99999., -99999., -99999.};
return quantities;
}
//construct the beam kinematics
ROOT::Math::PxPyPzMVector pair_4p(parts[first] + parts[second]);
ROOT::Math::PxPyPzMVector e1, P;
double e1_Energy = sqrt(10.*10. + get_pdg_mass("electron")*get_pdg_mass("electron"));
@@ -75,19 +137,20 @@ namespace util {
e1.SetPxPyPzE(0., 0., -10., e1_Energy);
P.SetPxPyPzE(0., 0., 100., P_Energy);
int scatteredIdx = -1;
float dp = 10.;
for(int i = 0 ; i < parts.size(); i++){
if(i==first || i==second) continue;
ROOT::Math::PxPyPzMVector k_prime(parts[i]);
if(i==first || i==second) continue; //skip the paired leptons
ROOT::Math::PxPyPzMVector k_prime(parts[i]); //scattered
float ptmp = sqrt(parts[i].px()*parts[i].px() + parts[i].py()*parts[i].py() + parts[i].pz()*parts[i].pz());
if( (k_prime.px()) * (pair_4p.px()) + (k_prime.py()) * (pair_4p.py()) + (k_prime.pz()) * (pair_4p.pz()) > 0. || ptmp >= 10.) continue; //angle between jpsi and scattered electron < pi/2, 3-momentum mag < 10.
if(dp > 10.- ptmp){ //if there are more than one candidate of scattered electron, choose the one with highest 3-momentum mag
//if( (k_prime.px()) * (pair_4p.px()) + (k_prime.py()) * (pair_4p.py()) + (k_prime.pz()) * (pair_4p.pz()) > 0. || ptmp >= 10.) continue; //angle between jpsi and scattered electron < pi/2, 3-momentum mag < 10.
//if(dp > 10.- ptmp){ //if there are more than one candidate of scattered electron, choose the one with highest 3-momentum mag
scatteredIdx = i;
dp = 10. - ptmp;
}
// dp = 10. - ptmp;
//}
}
if(scatteredIdx ==-1){
inv_quant quantities = {-999., -999., -999., -999.};
inv_quant quantities = {-99999., -99999., -99999., -99999.};
return quantities;
}
ROOT::Math::PxPyPzMVector q(e1 - parts[scatteredIdx]);
@@ -100,7 +163,6 @@ namespace util {
//inv_quant quantities = {-999., -999., -999., -999.};
return quantities;
}