TPC/groot/ReconstructTracks.C

0001 #include "ReconstructTracks.h"
0002 #include "groot.h"
0003 #include "ABlob.h"
0004 #include "ABlob.h"
0005
0006 #include "ATrack.h"
0007
0008 #include "TH2D.h"
0009 #include "TMath.h"
0010
0011 #include <iostream>
0012 #include <cmath>
0013
0014 using namespace std;
0015
0016
0017 static TH2D *trkDeltaOri; // accumalation of uncorrected data per run

0018 static TH2D *trkDeltaRes; // accumalation of corrected and normalized data per run

0019
0020 void ReconstructTracks()
0021 {
0022      groot* Tree = groot::instance();
0023      Tree->theTracks.clear();
0024
0025      //

0026      //  OK Track Fans.  Now that we have some studies

0027      //  done of the beam at FNAL, we find that the

0028      //  TRUE beam is pretty well columnated.  For this reason,

0029      //  we can a priori reject blob combinations that

0030      //  do not follow the standard pointing.

0031      //

0032      //  We shall apply hard-coded fits to the rough

0033      //  correlation between the track directions and

0034      //  thereby be able to make cuts that will result in

0035      //  only pretty good tracks.

0036      //

0037      //  Wish us luck!

0038      //

0039      //                           TKH 10-12-2013

0040      //

0041      //

0042
0043      static bool trkCreateOutput = true;
0044      double xmu, ymu, xsg, ysg, xysg;
0045
0046      if(trkCreateOutput) { // create new containers. Once per running program

0047        trkDeltaOri = new TH2D("trkDeltaOri","trkDeltaOri;dX (cm);dY (cm)",768,-10,+10,768,-10,+10);
0048        trkDeltaRes = new TH2D("trkDeltaRes","trkDeltaRes;Normalized dX; Normalized dY",120,-6,+6,120,-6,+6);
0049        trkCreateOutput = false;
0050      }
0051
0052      // gathering blob info

0053      const int mxcl1 = 40;
0054      const int mxcl2 = 90; //mxcl2>mxcl1!

0055      int cl1[mxcl1], cl2[mxcl2];
0056      int ncl1=0, ncl2=0;
0057      for( int i=0; i< Tree->theBlob.size(); i++)
0058        {
0059      if ( Tree->theBlob.at( i )->ID() == 0 )
0060        cl1[ncl1++] = i;
0061      if ( Tree->theBlob.at( i )->ID() == 1 )
0062        cl2[ncl2++] = i;
0063      if((ncl1+1)>mxcl1)
0064        {
0065          printf("We found more than %d blobs!!\n",mxcl1);
0066          break;
0067        }
0068      if((ncl2+1)>mxcl2)
0069        {
0070          printf("We found more than %d blobs!!\n",mxcl2);
0071          break;
0072        }
0073        }
0074      // making all combinations

0075      xysg = ReconstructTracks_Ref(xmu,xsg,ymu,ysg);
0076      double var[mxcl1][mxcl2];
0077      double vc3[mxcl1][mxcl2];
0078      for(int i=0; i!=ncl1; ++i)
0079        for(int j=0; j<ncl2; j++) {
0080      double x1 = Tree->theBlob[cl1[i]]->BestX();
0081      double y1 = Tree->theBlob[cl1[i]]->BestY();
0082      double x2 = Tree->theBlob[cl2[j]]->BestX();
0083      double y2 = Tree->theBlob[cl2[j]]->BestY();
0084      double c1 = Tree->theBlob[cl1[i]]->GetNSigma();
0085      double c2 = Tree->theBlob[cl2[j]]->GetNSigma();
0086      double dx = x2-x1;
0087      double dy = y2-y1;
0088      trkDeltaOri->Fill(dx/10000,dy/10000); //in cm

0089      double nSigmaX = (dx-xmu)/xsg;
0090      double nSigmaY = (dy-ymu)/ysg;
0091      trkDeltaRes->Fill(nSigmaX,nSigmaY);
0092      double c3 = nSigmaX*nSigmaX + nSigmaY*nSigmaY;
0093      var[i][j] = c1*c1+c2*c2+c3;
0094      vc3[i][j] = TMath::Sqrt(c3);
0095        }
0096      // sorting

0097      int nMAXtrk=TMath::Min(ncl1,ncl2);
0098      int trk[mxcl1][2];
0099      int ntrk=0;
0100      for(;ntrk!=nMAXtrk;++ntrk) {
0101        double min = 1e12; // to save time and hopping no nsigma is bigger than this =)

0102        int pi, pj;
0103        double c3;
0104        for(int i=0; i!=ncl1; ++i) {
0105      bool skip = false;
0106      for(int n=0; n!=ntrk; ++n) if(trk[n][0]==cl1[i]) skip=true; // skipping combinations with used blobs from TR1

0107      if(skip) continue;
0108      //cout << "loop " << cl1[i] << endl;

0109      for(int j=0; j!=ncl2; ++j) {
0110        skip=false;
0111        for(int n=0; n!=ntrk; ++n) if(trk[n][1]==cl2[j]) skip=true; // skipping combinations with used blobs from TR2

0112        if(skip) continue;
0113        //cout << "  loop " << cl2[j] << " quality^2 " << var[i][j] << endl;

0114        if(var[i][j]<min) {
0115          min = var[i][j];
0116          c3 = vc3[i][j];
0117          pi = cl1[i];
0118          pj = cl2[j];
0119        }
0120      }
0121        }
0122        trk[ntrk][0] = pi;
0123        trk[ntrk][1] = pj;
0124        vector<ABlob*> tempvec;
0125        tempvec.push_back( Tree->theBlob.at( pi ) ); //tr1

0126        tempvec.push_back( Tree->theBlob.at( pj ) ); //tr2

0127        Tree->theTracks.push_back ( new ATrack( tempvec,
0128                           Tree->theBlob.at(pi)->GetNSigma(),
0129                           Tree->theBlob.at(pj)->GetNSigma(),
0130                           c3) );
0131        //cout << "Found track " << ntrk << " made out of blobs " << pi << " and " << pj << " with quality " << Tree->theTracks.at(ntrk)->GetOverallQuality() << endl;

0132      }
0133
0134 }