TPC/groot/HoughBinSize.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 <iostream>
0009 #include <cmath>
0010
0011 TH2D* HoughHistMC = 0;
0012
0013 using namespace std;
0014
0015 void FillHoughHist()
0016 {
0017   if (!HoughHistMC)
0018     {
0019       HoughHistMC = new TH2D("HoughHistMC", "HoughHistMC" , 15 , -2.5 , 2.5, 12, 200, 2500);
0020     }
0021   //HoughHistMC->Reset();
0022   groot* Tree=groot::instance();
0023   int Npoints=0, k1=0; //number of points/blobs in an event. considering all tracks.
0024   for (int i=0; i<Nr;i++)
0025     {
0026       Npoints+=Tree->theBlobs[i].size();
0027       if (Tree->theBlobs[i].size()==1)
0028     k1++;
0029     }
0030
0031   cout<<"checking............";
0032   double points[2][k1];
0033   for (int i=0; i<k1+1; i++)
0034     {
0035       points[0][i]=0;
0036       points[1][i]=0;
0037     }
0038   cout<<"\n\n\n\t\t\t"<<Npoints<<"\n\n";
0039   int k=0;
0040   for (int i=0; i<Nr; i++)
0041     {
0042       for (int j=0; j< Tree->theBlobs[i].size(); j++)
0043     { if (Tree->theBlobs[i].size()==1)
0044         {
0045           points[0][k] = (Tree->theBlobs[i])[j]->CentroidX(); //X value
0046           points[1][k]= (Tree->theBlobs[i])[j]->CentroidY(); // Y value
0047           //cout<<"\n"<<points[0][k]<<"\t\t"<<points[1][k];
0048           k++;
0049         }
0050
0051     }
0052     }
0053
0054
0055   cout <<"\n N" << Npoints <<"\t k1"<<k1<<"\t K"<<k;
0056   double m, c;
0057
0058   static std::vector <double> inverse_slope;
0059   static std::vector <double> offset;
0060   double y1,y2,x1,x2;
0061   for (int i=0; i<k; i++)
0062     {
0063       for (int j=i+1; j<k; j++)
0064     {
0065       x1=points[0][i];
0066       x2=points[0][j];
0067       y1=points[1][i];
0068       y2=points[1][j];
0069       m= (y2-y1)/(x2-x1);
0070       c= (y2*x1 - y1*x2)/(x1-x2);
0071       inverse_slope.push_back(1/m);
0072       offset.push_back(c);
0073       //HoughHistMC->Fill(1/m,c);
0074       //cout<<"\n"<<m<<"\t"<<c;
0075     }
0076     }
0077   int size = static_cast<int>(inverse_slope.size());
0078   cout <<"\n Size "<< size;
0079
0080   //Calculating Median of Slope and Offset for Bin Size determination.
0081   double Median_InverseSlope, Median_Offset;
0082   sort(inverse_slope.begin(), inverse_slope.end());
0083   if (size  % 2 != 0)
0084     Median_InverseSlope = inverse_slope[(size+1)/2];
0085   else
0086     Median_InverseSlope =  (inverse_slope[(size)/2] + inverse_slope[(size+2)/2])/2.0 ;
0087   cout <<"\n Median_ Inverse Slope"<<Median_InverseSlope;
0088
0089   sort(offset.begin(), offset.end());
0090   if (size  % 2 != 0)
0091     Median_Offset = offset[(size+1)/2];
0092   else
0093     Median_Offset =  (offset[(size)/2] + offset[(size+2)/2])/2.0 ;
0094   cout <<"\nMedian_Offset"<<Median_Offset <<"\n";
0095   /*
0096   for(int i=0; i<size; ++i)
0097   std::cout <<inverse_slope[i] << ' ';*/
0098
0099
0100   //Calculating Std.devaition of Slope and Offset for Bin Size determination.
0101
0102
0103
0104
0105   std::vector <double> accum_IS;
0106   std::for_each (std::begin(inverse_slope), std::end(inverse_slope), [&](const double d) {
0107       accum_IS.push_back(abs(d - Median_InverseSlope));
0108     });
0109   cout<<"\n Size of accum vector without static" << accum_IS.size();
0110   /*
0111   cout<<"\n the median deviation list\n";
0112   for(int i=0; i<size; ++i)
0113   std::cout <<accum[i] << ' ';*/
0114
0115   double MAD_InverseSlope = 0;
0116   sort(accum_IS.begin(), accum_IS.end());
0117   if (size  % 2 != 0)
0118     MAD_InverseSlope = accum_IS[(size+1)/2];
0119   else
0120     MAD_InverseSlope = accum_IS[size/2];
0121     //Med_dev =  (accum[(size)/2] + accum[(size+2)/2])/2.0 ;
0122
0123   cout <<"\n MAD_InverseSlope"<< MAD_InverseSlope;
0124
0125
0126   std::vector <double> accum_O;
0127   std::for_each (std::begin(inverse_slope), std::end(inverse_slope), [&](const double d) {
0128       accum_O.push_back(abs(d - Median_Offset));
0129     });
0130
0131   cout<<"\n Size of accum vector with static" << accum_O.size();
0132   /*
0133   cout<<"\n the median deviation list\n";
0134   for(int i=0; i<size; ++i)
0135   std::cout <<accum[i] << ' ';*/
0136
0137   double MAD_Offset = 0;
0138   sort(accum_O.begin(), accum_O.end());
0139   if (size  % 2 != 0)
0140     MAD_Offset = accum_O[(size+1)/2];
0141   else
0142     MAD_Offset = accum_O[size/2];
0143     //Med_dev =  (accum[(size)/2] + accum[(size+2)/2])/2.0 ;
0144
0145   cout <<"\n MAD_Offset"<< MAD_Offset;
0146
0147
0148
0149
0150
0151
0152   //ATrack *theTracks;
0153   //int binmax = HoughHistMC->GetMaximumBin();
0154   //double Var[2];
0155   //Var[0] = 1/(HoughHistMC->GetXaxis()->GetBinCenter(binmax));
0156   //Var[1] = HoughHistMC->GetYaxis()->GetBinCenter(binmax);
0157   //Tree->theTracks->SetSlope(Var[0]);
0158   //Tree->theTracks->SetOffset(Var[1]);
0159   // cout<<"\nFinal\t"<<Tree->theTracks.Slope<<"\t"<<Tree->theTracks.Offset;
0160   //cout <<"\ncheck........"<<(Tree->theTracks)->check();
0161
0162
0163 }