See SN0089 for the full write-up of 3d method.

Write-up of 2d method is at:

Code is at:

/star/u/rjreed/PPV2009/v1/CenterFinderv2.C
 

I still need to fix it so that the code accepts the proper data structure, but it appears to be working on a test case of 1.

 for (int i = 0;i<N;++i){
  
    float x;float y;float px;float py;float sig;
//read in values here

    float pt = sqrt(px*px+py*py);
    cout<<"pt = "<<pt<<endl;
    float px_u = px/pt; //formula uses unit vector in p
    float py_u = py/pt;
    cout<<"px_u = "<<px_u<<" py_u = "<<py_u<<endl;
 
    M[0][0] = M[0][0] + py_u*py_u/(sig*sig);
    M[0][1] = M[0][1] - px_u*py_u/(sig*sig);
    M[1][0] = M[1][0] - px_u*py_u/(sig*sig);
    M[1][1] = M[1][1] + px_u*px_u/(sig*sig);

    b[0] = b[0] + (x*py_u*py_u - y*px_u*py_u)/(sig*sig);
    b[1] = b[1] + (-x*px_u*py_u + y*px_u*px_u)/(sig*sig);}
  for(int i =0;i<2;i++){
    for(int j = 0;j<2;j++){
      cout<<"M["<<i<<"]["<<j<<"]= "<<M[i][j]<<endl;}}
  cout<<"b = "<<b[0]<<" "<<b[1]<<endl;
  float Minv[2][2];
  float detM = M[0][0]*M[1][1]-M[1][0]*M[0][1];
  Minv[0][0] = M[1][1]/detM;
  Minv[1][1] = M[0][0]/detM;
  Minv[1][0] = -M[1][0]/detM;
  Minv[0][1] = -M[0][1]/detM;

  for(int i =0;i<2;i++){
    for(int j = 0;j<2;j++){
      cout<<"Minv["<<i<<"]["<<j<<"]= "<<Minv[i][j]<<endl;}}

  float Vx = Minv[0][0]*b[0]+Minv[0][1]*b[1];
  float Vy = Minv[1][0]*b[0]+Minv[1][1]*b[1];

Results from running one test case:

x = -1.05 y = 0.48 px = 6.27 py = 13.51
pt = 14.8941
px_u = 0.420973 py_u = 0.907073
x = -0.02 y = -0.3 px = -14.01 py = 1.17
pt = 14.0588
px_u = -0.996531 py_u = 0.0832221
x = -1 y = -0.69 px = 5.39 py = -7.81
pt = 9.48937
px_u = 0.568004 py_u = -0.823026
x = -0.76 y = -0.52 px = -5.65 py = 8.27
pt = 10.0158
px_u = -0.564111 py_u = 0.825699
x = 0.03 y = 0.09 px = 25.33 py = -8.03
pt = 26.5723
px_u = 0.953247 py_u = -0.302194
M[0][0]= 181.984
M[0][1]= 14.5591
M[1][0]= 14.5591
M[1][1]= 132.123
b = -244.106 -38.4834
Minv[0][0]= 0.00554387
Minv[0][1]= -0.000610898
Minv[1][0]= -0.000610898
Minv[1][1]= 0.007636

Vx = -1.32979 Vy = -0.144735

Geant vertex is at:  GVER z=89.7 x=-1.22 y=-0.18

To Do:

Allow proper input

Add switches for cuts on p and Rxy

Itterate to remove outliers