15 {
16
17
18 double rmol = 10.00;
19 double dzchar = 100.0;
20
21
22 double aE = a.centroid().E();
23 double aX = a.centroid().Px();
24 double aY = a.centroid().Py();
25 double aZ = a.centroid().Pz();
26
27 double bE = b.centroid().E();
28 double bX = b.centroid().Px();
29 double bY = b.centroid().Py();
30 double bZ = b.centroid().Pz();
31
32 double dijz;
33 if (aE >= bE) {
34 dijz = bZ - aZ;
35 } else {
36 dijz = aZ - bZ;
37 }
38
39 double dijT = pow(pow(aX - bX, 2) + pow(aY - bY, 2), 0.5);
40
41 double weightT = exp(pow(dijT / rmol, 2)) - 1;
42 double weightZ = (exp(abs(dijz) / dzchar) - 1);
43
44
45 if (weightT <= weightZ) {
46 return weightZ;
47 } else {
48 return weightT;
49 }
50 }