math/photosensorSphere/FindSphereCenter.C

0001 Float_t sphx,sphy,sphz;
0002 TTree *tr;
0003 TH1D *hist;
0004
0005 // given starting coordinates `start{x,y,z}`, scan through a
0006 // lattice of points, with size `size`; for each lattice point,
0007 // compute the RMS of the distance from that point to all the
0008 // photosensors; the closer this lattice point is to the spherical
0009 // center, the smaller the RMS will be; coordinates `end{x,y,z}` will
0010 // be the best lattice point, `enddist` will be the mean distance
0011 // too the photosensor points, and `endrms` will be the rms
0012 void FindCenter(Float_t startx, Float_t starty, Float_t startz,
0013   Float_t size,
0014   Float_t &endx, Float_t &endy, Float_t &endz,
0015   Float_t &enddist, Float_t &endrms) {
0016
0017   endrms = 1e6;
0018   enddist = 1e6;
0019
0020   Float_t currx,curry,currz;
0021   Float_t currrms;
0022
0023   int ndiv = 5; // number of lattice points in 1 dimension; keep it odd
0024   Float_t div = size/ndiv; // lattice step size
0025   int rc = (ndiv-1)/2;
0026
0027   // lattice loop
0028   for(int x=0; x<ndiv; x++) {
0029   for(int y=0; y<ndiv; y++) {
0030   for(int z=0; z<ndiv; z++) {
0031
0032     // current lattice point
0033     currx = startx + (x-rc)*div;
0034     curry = starty + (y-rc)*div;
0035     currz = startz + (z-rc)*div;
0036     //printf("%f %f %f\n",currx,curry,currz);
0037
0038     // calculate RMS of distances between this
0039     // lattice point and the photosensors
0040     hist->Reset();
0041     for(int e=0; e<tr->GetEntries(); e++) {
0042       tr->GetEntry(e);
0043       hist->Fill(
0044           TMath::Sqrt(
0045             TMath::Power(sphx-currx,2)+
0046             TMath::Power(sphy-curry,2)+
0047             TMath::Power(sphz-currz,2)
0048             )
0049           );
0050     };
0051     currrms = hist->GetStdDev();
0052     //printf("   rms = %f\n",currrms);
0053
0054     // if this is the best RMS, set result values
0055     if(currrms<endrms) {
0056       endx = currx;
0057       endy = curry;
0058       endz = currz;
0059       endrms = currrms;
0060       enddist = hist->GetMean();
0061     };
0062
0063   };};}; // end lattice loop
0064
0065   // print result
0066   printf("result: (%.2f, %.2f, %.2f)  rms=%.4f  dist=%.4f\n",
0067     endx,endy,endz,endrms,enddist);
0068 };
0069
0070 //////////////////////////////////////////////////
0071
0072 void FindSphereCenter(TString treeFile="photosensors.dat") {
0073
0074   // read photosensor cooridnates
0075   tr = new TTree();
0076   tr->ReadFile(treeFile,"x/F:y/F:z/F");
0077   tr->SetBranchAddress("x",&sphx);
0078   tr->SetBranchAddress("y",&sphy);
0079   tr->SetBranchAddress("z",&sphz);
0080
0081   // histogram for calculating RMS
0082   hist = new TH1D("hist","hist",2000,0,2000);
0083
0084   Float_t guessx,guessy,guessz;
0085   Float_t bestx,besty,bestz;
0086   Float_t size,dist,rms;
0087
0088   // hyperparameters
0089   const int NIter = 1000; // number of iterations
0090   Float_t initsize = 300; // initial lattice size
0091   guessx = 0; // initial guess
0092   guessy = 0;
0093   guessz = 0;
0094
0095   // iteratively search for the center: start with
0096   // initial guess and lattice size, and get the best lattice point;
0097   // then start with this best lattice point, and search again with
0098   // a lattice half the size; since each subsequent lattice is smaller
0099   // and the number of lattice points is constant, each lattice will
0100   // be more dense, and hopefully we will converge on the true sphere
0101   // center and radius
0102   for(int d=0; d<NIter; d++) {
0103     size = d>0 ? initsize/(2*d) : initsize;
0104     FindCenter(
0105       guessx,guessy,guessz,
0106       size,
0107       bestx,besty,bestz,
0108       dist,rms
0109       );
0110     guessx = bestx;
0111     guessy = besty;
0112     guessz = bestz;
0113   };
0114 };