Skip to content
Snippets Groups Projects
Commit 4ff1c087 authored by Ziyue Zhang's avatar Ziyue Zhang
Browse files

WIP: Replacing old functions

parent 1f5caf96
No related branches found
No related tags found
No related merge requests found
This commit is part of merge request !39. Comments created here will be created in the context of that merge request.
Hide whitespace changes
Inline Side-by-side
benchmarks/dvmp/analysis/dvmp.h
+ 8
− 93
Edit View file @ 4ff1c087
......@@ -19,28 +19,9 @@ namespace util {
//========================================================================================================
// for structure functions
struct inv_quant { // add more when needed
double nu, Q2, x, y, t;
};
// for simu
/*inline inv_quant calc_inv_quant_sim(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]);//exact P' - P
//ROOT::Math::PxPyPzMVector Delta(parts[0] - parts[2] - parts[5]);//e - e' - jpsi
//ROOT::Math::PxPyPzMVector Delta(parts[0] - parts[2] - parts[7] - parts[8]);//jpsi->l l' gamma, ignore gamma
//P' - P = e - e' - jpsi, jpsi != l + l' (there could be an additional photon)
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;
}*/
//import Geant4 and set the wanted particles in the intended order
//0:e0 1:p0 2:e1 3:p1 4:recoil system (without p1) 5:l1 from 4 6:l2 from 4
inline auto momenta_sort_sim(const std::vector<dd4pod::Geant4ParticleData>& parts, std::string_view mother, std::string_view daughter){//mother and daughter are not used yet; will be useful when generater is different and/or when the mcparticles doesn't follow the same order in all events
......@@ -54,10 +35,10 @@ namespace util {
double e = sqrt(px*px + py*py + pz*pz + mass*mass);
momenta[i].SetPxPyPzE(px, py, pz, e);
}
for(int i = 0 ; i < 7 ; i++){
cout<<Form("sim, idx = %d, P(px, py, pz, mass) = (%f, %f, %f, %f)", i, momenta[i].px(), momenta[i].py(), momenta[i].pz(), momenta[i].mass())<<endl;
}
cout<<"========================================="<<endl;
//for(int i = 0 ; i < 7 ; i++){
// cout<<Form("sim, idx = %d, P(px, py, pz, mass) = (%f, %f, %f, %f)", i, momenta[i].px(), momenta[i].py(), momenta[i].pz(), momenta[i].mass())<<endl;
//}
//cout<<"========================================="<<endl;
return momenta;
}
......@@ -132,79 +113,13 @@ namespace util {
// dp = 10. - ptmp;
//}
}
for(int i = 0 ; i < 7 ; i++){
cout<<Form("dum, idx = %d, P(px, py, pz, mass) = (%f, %f, %f, %f)", i, momenta[i].px(), momenta[i].py(), momenta[i].pz(), momenta[i].mass())<<endl;
}
cout<<"========================================="<<endl;
//for(int i = 0 ; i < 7 ; i++){
// cout<<Form("dum, idx = %d, P(px, py, pz, mass) = (%f, %f, %f, %f)", i, momenta[i].px(), momenta[i].py(), momenta[i].pz(), momenta[i].mass())<<endl;
//}
//cout<<"========================================="<<endl;
return momenta;
}
//for Dummy rc
/*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;
// 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) {
if( fabs(parts[i].mass() - daughter_mass)/daughter_mass > 0.01) continue;
for (int j = i + 1; j < parts.size(); ++j) {
if( fabs(parts[j].mass() - daughter_mass)/daughter_mass > 0.01) continue;
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 ){ //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"));
//double P_Energy = sqrt(100.*100. + get_pdg_mass("proton")*get_pdg_mass("proton"));
double e1_Energy = sqrt((1.305e-8 - 10.)*(1.305e-8 - 10.) + (0.0005109+9.888e-8)*(0.0005109+9.888e-8));
double P_Energy = sqrt((99.995598 + 1.313e-7)*(99.995598 + 1.313e-7) + (0.938272-1.23e-12)*(0.938272-1.23e-12));
e1.SetPxPyPzE(0., 0., 1.305e-8 - 10., e1_Energy);
P.SetPxPyPzE(0., 0., 99.995598 + 1.313e-7, P_Energy);
int scatteredIdx = -1;
//float dp = 10.;
for(int i = 0 ; i < parts.size(); i++){
if(i==first || i==second) continue; //skip the paired leptons
//if( fabs(parts[i].mass() - get_pdg_mass("electron"))/get_pdg_mass("electron") > 0.01) continue;
if( fabs(parts[i].mass() - 0.0005109 - 9.888e-8)/(0.0005109 + 9.888e-8) > 0.01) continue;
//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
scatteredIdx = i;
// dp = 10. - ptmp;
//}
}
if(scatteredIdx ==-1){
inv_quant quantities = {-99999., -99999., -99999., -99999.};
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};
return quantities;
}*/
inline inv_quant calc_inv_quant(const std::vector<ROOT::Math::PxPyPzMVector>& parts)
{
//0:e0 1:p0 2:e1 3:p1 4:recoil system (without p1) 5:l1 from 4 6:l2 from 4
......
benchmarks/dvmp/analysis/vm_invar.cxx
+ 4
− 9
Edit View file @ 4ff1c087
......@@ -76,25 +76,20 @@ int vm_invar(const std::string& config_name)
};*/
auto momenta_sort_sim = [vm_name, decay_name](const std::vector<dd4pod::Geant4ParticleData>& parts){
return util::momenta_sort_sim(parts, vm_name, decay_name);
};
auto momenta_sort_rec = [vm_name, decay_name](const std::vector<eic::ReconstructedParticleData>& parts){
};
return util::momenta_sort_rec(parts, vm_name, decay_name);
};
//====================================================================
// Define analysis flow
auto d_im = d.Define("p_rec", util::momenta_RC, {"DummyReconstructedParticles"})
.Define("p_sim", util::momenta_from_simulation, {"mcparticles2"})
auto d_im = d.//Define("p_rec", util::momenta_RC, {"DummyReconstructedParticles"})
//.Define("p_sim", util::momenta_from_simulation, {"mcparticles2"})
.Define("p_rec_sorted", momenta_sort_rec, {"DummyReconstructedParticles"})
.Define("p_sim_sorted", momenta_sort_sim, {"mcparticles2"})
.Define("N", "p_rec.size()")
.Define("N", "p_rec_sorted.size()")
.Define("invariant_quantities_rec", util::calc_inv_quant, {"p_rec_sorted"})
.Define("invariant_quantities_sim", util::calc_inv_quant, {"p_sim_sorted"})
//.Define("p_sim", util::momenta_from_simulation, {"mcparticles2"})
//================================================================
//.Define("invariant_quantities_rec", calc_inv_quant_rec, {"p_rec"})
......
options/tracker_reconstruction.py
+ 1
− 1
Edit View file @ 4ff1c087
......@@ -57,7 +57,7 @@ podioinput = PodioInput("PodioReader",
dummy = MC2DummyParticle("MC2Dummy",
inputCollection="mcparticles",
outputCollection="DummyReconstructedParticles",
smearing = 0.0)
smearing = 0.1)
## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
copier = MCCopier("MCCopier",
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment