TPC/groot/FillHoughHist.C

0001 #include "FillHoughHist.h"
0002 #include "TH1D.h"
0003 #include "TH2D.h"
0004 #include "ABlob.h"
0005 #include "ATrack.h"
0006 #include "groot.h"
0007
0008 #include "AZigzag.h"
0009
0010 #include <iostream>
0011 #include <cmath>
0012
0013 TH2D* HoughHistMC  = 0; //  This is in relative coordinates to "tune" the Hough Space.
0014 TH2D* HoughHistABS = 0; //  This is in absolute coordinates to "solve" the Hough Space.
0015
0016 using namespace std;
0017
0018 void FillHoughHist()
0019 {
0020   groot* Tree=groot::instance();
0021
0022   if (!HoughHistMC)
0023     {
0024       HoughHistMC = new TH2D("HoughHistMC", "HoughHistMC" , 10, (Med_inverseSlope - 10*MAD_inverseSlope),(Med_inverseSlope + 10*MAD_inverseSlope), 10 , (Med_Offset - 10*MAD_Offset), (Med_Offset + 10*MAD_Offset));
0025
0026
0027       double x00 = Tree->ZigzagMap2[0][0]->XCenter();
0028       double y00 = Tree->ZigzagMap2[0][0]->YCenter();
0029       double x11 = Tree->ZigzagMap2[Nr][Nphi]->XCenter();
0030       double y11 = Tree->ZigzagMap2[Nr][Nphi]->YCenter();
0031
0032       double mi = abs(inverseSlope(x00, y00, x11, y11));
0033       double c  = abs(intercept(x00, y00, x11, y11));
0034
0035       int NHOUGHX = 2*mi/(HFACTOR*MAD_inverseSlope);
0036       int NHOUGHY = 2*c/(HFACTOR*MAD_Offset);
0037
0038       HoughHistABS = new TH2D("HoughHistABS", "HoughHistABS" , NHOUGHX, -mi, mi, NHOUGHY, -c, c);
0039     }
0040   HoughHistMC->Reset();
0041   HoughHistABS->Reset();
0042   int Npoints=0;
0043   for (int i=0; i<Nr;i++)
0044     {
0045       Npoints+=Tree->theBlobs[i].size();
0046     }
0047
0048   //cout<<"checking............";
0049   double points[2][Npoints];
0050   for (int i=0; i<Npoints; i++)
0051     {
0052       points[0][i]=0;
0053       points[1][i]=0;
0054     }
0055   //cout<<"\n\n\n\t\t\t"<<Npoints<<"\n\n";
0056   int k=0;
0057   for (int i=0; i<Nr; i++)
0058     {
0059       for (int j=0; j< Tree->theBlobs[i].size(); j++)
0060     {
0061       points[0][k] = (Tree->theBlobs[i])[j]->CentroidX(); //X value
0062       points[1][k]= (Tree->theBlobs[i])[j]->CentroidY(); // Y value
0063       //cout<<"\n"<<points[0][k]<<"\t\t"<<points[1][k];
0064       k++;
0065     }
0066
0067     }
0068
0069
0070   vector<double> MI;
0071   vector<double> C;
0072   double y1,y2,x1,x2;
0073   for (int i=0; i<k; i++)
0074     {
0075       for (int j=i+1; j<k; j++)
0076     {
0077       x1=points[0][i];
0078       x2=points[0][j];
0079       y1=points[1][i];
0080       y2=points[1][j];
0081       double mi = inverseSlope(x1, y1, x2, y2);
0082       double c  = intercept   (x1, y1, x2, y2);
0083       HoughHistABS->Fill(mi, c);
0084       MI.push_back(mi);
0085       C.push_back(c);
0086     }
0087     }
0088
0089
0090   for (int i=0; i<k; i++)
0091     {
0092       for (int j=i+1; j<k; j++)
0093     {
0094       x1=points[0][i];
0095       x2=points[0][j];
0096       y1=points[1][i];
0097       y2=points[1][j];
0098       double mi = inverseSlope(x1, y1, x2, y2);
0099       double c  = intercept   (x1, y1, x2, y2);
0100       HoughHistMC->Fill(mi, c);
0101     }
0102     }
0103
0104   //ATrack *theTracks;
0105   //int binmax = HoughHistMC->GetMaximumBin();
0106   //double Var[2];
0107   //Var[0] = 1/(HoughHistMC->GetXaxis()->GetBinCenter(binmax));
0108   //Var[1] = HoughHistMC->GetYaxis()->GetBinCenter(binmax);
0109   //Tree->theTracks->SetSlope(Var[0]);
0110   //Tree->theTracks->SetOffset(Var[1]);
0111   // cout<<"\nFinal\t"<<Tree->theTracks.Slope<<"\t"<<Tree->theTracks.Offset;
0112   //cout <<"\ncheck........"<<(Tree->theTracks)->check();
0113
0114
0115 }
0116
0117 double inverseSlope(double x1, double y1, double x2, double y2)
0118 {
0119   double mi = (x2-x1)/(y2-y1);
0120   return mi;
0121 }
0122
0123 double intercept   (double x1, double y1, double x2, double y2)
0124 {
0125   double c= (y2*x1 - y1*x2)/(x1-x2);
0126   return c;
0127 }