Skip to content
Snippets Groups Projects

WIP: Dis test

Closed Whitney Armstrong requested to merge dis_test into master
1 unresolved thread
4 files
+ 421
45
Compare changes
  • Side-by-side
  • Inline
Files
4
+ 131
31
@@ -44,7 +44,9 @@ inline double get_pdg_mass(std::string_view part) {
return 3.0969;
} else if (part == "upsilon") {
return 9.49630;
} else {
} else if (part == "proton"){
return 0.938272;
}else {
throw unknown_particle_error{part};
}
}
@@ -117,44 +119,142 @@ find_decay_pair(const std::vector<ROOT::Math::PxPyPzMVector>& parts,
return {parts[first], parts[second]};
}
// Calculate the invariant mass of a given pair of particles
inline double
get_im(const std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>&
particle_pair) {
return (particle_pair.first + particle_pair.second).mass();
}
// Calculate the transverse momentum of a given pair of particles
inline double
get_pt(const std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>&
particle_pair) {
struct pair_kinematics{
double mass, Pt, phi, y;
};
inline pair_kinematics calc_pair_kinematics_rec(const std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>&
particle_pair){
double px_pair = (particle_pair.first + particle_pair.second).px();
double py_pair = (particle_pair.first + particle_pair.second).py();
return sqrt(px_pair*px_pair + py_pair*py_pair);
double py_pair = (particle_pair.first + particle_pair.second).py();
double pz_pair = (particle_pair.first + particle_pair.second).pz();
double mass_pair = (particle_pair.first + particle_pair.second).mass();
double energy_pair = sqrt(mass_pair*mass_pair + px_pair*px_pair + py_pair*py_pair + pz_pair*pz_pair);
pair_kinematics kinematics = {mass_pair, sqrt(px_pair*px_pair + py_pair*py_pair), std::atan2(py_pair,px_pair), 0.5*log((energy_pair + pz_pair)/(energy_pair - pz_pair))};
return kinematics;
}
inline pair_kinematics calc_pair_kinematics_sim(const std::vector<ROOT::Math::PxPyPzMVector>& parts){
ROOT::Math::PxPyPzMVector pair_4p(parts[7] + parts[8]);
double px_pair = pair_4p.px();
double py_pair = pair_4p.py();
double pz_pair = pair_4p.pz();
double mass_pair = pair_4p.mass();
double energy_pair = sqrt(mass_pair*mass_pair + px_pair*px_pair + py_pair*py_pair + pz_pair*pz_pair);
pair_kinematics kinematics = {mass_pair, sqrt(px_pair*px_pair + py_pair*py_pair), std::atan2(py_pair,px_pair), 0.5*log((energy_pair + pz_pair)/(energy_pair - pz_pair))};
return kinematics;
}
// Calculate the azimuthal angle of a given pair of particles
inline double
get_phi(const std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>&
particle_pair) {
double px_pair = (particle_pair.first + particle_pair.second).px();
double py_pair = (particle_pair.first + particle_pair.second).py();
double phi_pair = std::atan2(py_pair,px_pair);
//if(py_pair <= 0.) phi_pair = - phi_pair;
return phi_pair;
get_im(pair_kinematics kinematics) {
return kinematics.mass;
}
// Calculate the rapidity of a given pair of particles
inline double
get_y(const std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>&
particle_pair) {
double px_pair = (particle_pair.first + particle_pair.second).px();
double py_pair = (particle_pair.first + particle_pair.second).py();
double pz_pair = (particle_pair.first + particle_pair.second).pz();
double mass_pair = (particle_pair.first + particle_pair.second).mass();
double energy_pair = sqrt(mass_pair*mass_pair + px_pair*px_pair + py_pair*py_pair + pz_pair*pz_pair);
return 0.5*log((energy_pair + pz_pair)/(energy_pair - pz_pair));
get_pt(pair_kinematics kinematics) {
return kinematics.Pt;
}
inline double
get_phi(pair_kinematics kinematics) {
return kinematics.phi;
}
inline double
get_y(pair_kinematics kinematics) {
return kinematics.y;
}
//========================================================================================================
//for structure functions
struct inv_quant{ //add more when needed
double nu, Q2, x, t;
};
//for simu
inline inv_quant calc_inv_quant_simu(const std::vector<ROOT::Math::PxPyPzMVector>& parts){
ROOT::Math::PxPyPzMVector q(parts[0] - parts[2]);
ROOT::Math::PxPyPzMVector P(parts[3]);
ROOT::Math::PxPyPzMVector Delta(parts[6] - parts[3]);
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};
return quantities;
}
//for tracking
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;
}
inline double get_nu(inv_quant quantities) {
return quantities.nu/1000.;
}
inline double get_Q2(inv_quant quantities) {
return quantities.Q2;
}
inline double get_x(inv_quant quantities) {
return quantities.x;
}
inline double get_t(inv_quant quantities) {
return quantities.t;
}
//for tracking, add later
//=========================================================================================================
} // namespace util
Loading