DetectorAlignment.cxx

Go to the documentation of this file.

//
// DetectorAlignment
//
// Package: DetectorAlignment
//
// An alignment engine 
//
// Contact: rustem@fnal.gov
//
// Created on: Fri Sep 24 15:26:05 2004
//


//DetectorAlignment classes
#include "DetectorAlignment.h"
#include "NtpAlignmentRecord.h"
#include "NtpAlignTrackStrip.h"
#include "NtpAlignCandStrip.h"

//MINOSSOFT
#include "MessageService/MsgService.h"
#include "MessageService/MsgFormat.h"
#include "RecoBase/CandStripHandle.h"
#include "RecoBase/CandStripListHandle.h"
#include "CandTrackSR/CandTrackSRHandle.h"
#include "CandTrackSR/CandTrackSRListHandle.h"

//ROOT
#include "TClonesArray.h"

//C++
#include <algorithm>
#include <iostream>
#include <cmath>

CVSID("$Id: DetectorAlignment.cxx,v 1.5 2009/02/28 21:46:10 gmieg Exp $");

using std::endl;

DetectorAlignment::DetectorAlignment()
   :f2dTrackCharge(0.0),
    fCandVCharge(0.0),
    fCandUCharge(0.0),
    f2dTrackResidualRMS(0.0),
    fUncertaintyInTPosofStrips(0.0),
    fFita(0.0),
    fFitb(0.0),
    fSigmaOfa(0.0),
    fSigmaOfb(0.0),
    fCovab(0.0),
    fDeltaT(0.0)
{
   MSG("Align", Msg::kDebug) << "Constructor DetectorAlignment()" << endl;
}

bool DetectorAlignment::RunAlignment(const CandTrackSRHandle* tkh, 
                                     const CandStripListHandle* cslh,
                                     NtpAlignmentRecord* ntpalignrec)
{

   if(!tkh || !cslh || !ntpalignrec) return false;

   //clear vectors that store AlignmentStrips and track's time inteval
   fTrackTimes.clear();

   //Save strip data from CandTrackSR to vector<AlignmentStrip>, one for each view
   GetTrackStrips(tkh);

   //Save data from candidate strips to vector<AlignmentStrip>, one for each view
   //Ignore strips that belong to a track
   GetCandStrips(cslh);

   //Get directional cosines for 3d track
   ntpalignrec -> hcosu = tkh -> GetHoughDirCosU();
   ntpalignrec -> hcosv = tkh -> GetHoughDirCosV();
   ntpalignrec -> hcosz = tkh -> GetHoughDirCosZ();

   //Fit V view track with a straight line and save fit information
   vector<AlignmentStrip>::iterator begitr = fTrackVStrip.begin();
   vector<AlignmentStrip>::iterator enditr = fTrackVStrip.end();

   //V view: calculate residuals, track charge and residual rms
   MakeAlignmentTrack(begitr, enditr);

   FitTrackLessOne(begitr, enditr, enditr);

   //Fill ntuples record for V view track
   ntpalignrec -> vcharge      = f2dTrackCharge;
   ntpalignrec -> vcandcharge  = fCandVCharge;
   ntpalignrec -> vtrackrms    = f2dTrackResidualRMS;
   ntpalignrec -> vfita        = fFita;
   ntpalignrec -> vfitb        = fFitb;
   ntpalignrec -> vsigmaofa    = fSigmaOfa;
   ntpalignrec -> vsigmaofb    = fSigmaOfb;
   ntpalignrec -> vsigmaoftpos = fUncertaintyInTPosofStrips;
   ntpalignrec -> vcovab       = fCovab;
   ntpalignrec -> vdt          = fDeltaT;
   

   //Fit U view track with a straight line and save fit information
   begitr = fTrackUStrip.begin(); 
   enditr = fTrackUStrip.end();
   
   //U view: calculate residuals, track charge and residual rms
   MakeAlignmentTrack(begitr, enditr);

   FitTrackLessOne(begitr, enditr, enditr);

   //Fill ntuples record for U view track
   ntpalignrec -> ucharge      = f2dTrackCharge;
   ntpalignrec -> ucandcharge  = fCandUCharge;
   ntpalignrec -> utrackrms    = f2dTrackResidualRMS;
   ntpalignrec -> ufita        = fFita;
   ntpalignrec -> ufitb        = fFitb;
   ntpalignrec -> usigmaofa    = fSigmaOfa;
   ntpalignrec -> usigmaofb    = fSigmaOfb;
   ntpalignrec -> usigmaoftpos = fUncertaintyInTPosofStrips;
   ntpalignrec -> ucovab       = fCovab;
   ntpalignrec -> udt          = fDeltaT;

   //Fill TClonesArrays with candidate strips and track strips
   TClonesArray& trackvstrips = *(ntpalignrec -> trackvstrip);
   for(unsigned int i = 0; i < fTrackVStrip.size(); ++i){
      NtpAlignTrackStrip* ntpalign = new(trackvstrips[i])NtpAlignTrackStrip(fTrackVStrip[i]);
      if(!ntpalign)
         MSG("Align", Msg::kError) <<"Failed to fill TClonesArray!" << endl;
   }
   
   TClonesArray& trackustrips = *(ntpalignrec -> trackustrip);
   for(unsigned int i = 0; i < fTrackUStrip.size(); ++i){
      NtpAlignTrackStrip* ntpalign = new(trackustrips[i])NtpAlignTrackStrip(fTrackUStrip[i]);
      if(!ntpalign)
         MSG("Align", Msg::kError) <<"Failed to fill TClonesArray!" << endl;
   }
   
   TClonesArray& candvstrips = *(ntpalignrec -> candvstrip);
   for(unsigned int i = 0; i < fCandVStrip.size(); ++i){
      NtpAlignCandStrip* ntpalign = new(candvstrips[i])NtpAlignCandStrip(fCandVStrip[i]);
      if(!ntpalign)
         MSG("Align", Msg::kError) <<"Failed to fill TClonesArray!" << endl;
   }

   TClonesArray& candustrips = *(ntpalignrec -> candustrip);
   for(unsigned int i = 0; i <  fCandUStrip.size(); ++i){
      NtpAlignCandStrip* ntpalign = new(candustrips[i])NtpAlignCandStrip(fCandUStrip[i]);
      if(!ntpalign)
         MSG("Align", Msg::kError) <<"Failed to fill TClonesArray!" << endl;
   }
   
   ntpalignrec -> nstrip       = cslh -> GetNDaughters();
   ntpalignrec -> ntrackvstrip = fTrackVStrip.size();
   ntpalignrec -> ntrackustrip = fTrackUStrip.size();
   ntpalignrec -> ncandvstrip  = fCandVStrip.size();
   ntpalignrec -> ncandustrip  = fCandUStrip.size();

   sort(fTrackTimes.begin(), fTrackTimes.end());

   if(fTrackTimes.size() > 0)
      ntpalignrec ->  dt = fTrackTimes[fTrackTimes.size()-1]-fTrackTimes[0];
   else
      ntpalignrec ->  dt = 0.0;

   const int sum = fTrackVStrip.size() + fTrackUStrip.size() + fCandVStrip.size() + fCandUStrip.size();
   const int diff = cslh->GetNDaughters() - sum;
   if(diff != 0){
      MSG("Align", Msg::kDebug)<< "# of cand and track strips do not add up! "
                               << "diff = " << diff << endl;
      return false;
   }   
   return true;
}


//Get charges of hit strips in CandTrack and fill corresponding histograms
void DetectorAlignment::GetTrackStrips(const CandTrackSRHandle* tkh)
{
   MSG("Align", Msg::kVerbose) << "Getting CandStripHandles from CandTrack: " 
                               << "tkh->GetNDaughters() = "
                               << tkh->GetNDaughters()
                               << ", tkh->GetTitle() = " << tkh->GetTitle()<<endl;

   fTrackVStrip.clear();
   fTrackUStrip.clear();

   int nstrip = 0;
   TIter sitr(tkh->GetDaughterIterator());
   while( CandStripHandle* csh = dynamic_cast<CandStripHandle*>(sitr()) ){
      AlignmentStrip astrip(csh);
      nstrip++;
      MSG("Align", Msg::kVerbose) << "Adding " <<PlexStripEndId(astrip.plexseid) 
                                  <<" to track record" <<endl;
      if( csh -> GetPlaneView() == PlaneView::kV)
         fTrackVStrip.push_back(astrip);
      else
         fTrackUStrip.push_back(astrip);      
   }

   MSG("Align", Msg::kVerbose) <<"Added "<<nstrip<<" track strips." << endl;
   return;
}


void DetectorAlignment::GetCandStrips(const CandStripListHandle* cslh)
{
   MSG("Align", Msg::kVerbose) << "Separating CandStripHandles from track strips" <<endl
                               << "cslh->GetNDaughters() = "<< cslh -> GetNDaughters() 
                               << ", cslh->GetTitle() = " << cslh -> GetTitle() <<endl;
   fCandVStrip.clear();
   fCandUStrip.clear();   
   
   fCandVCharge = 0.0;
   fCandUCharge = 0.0;

   MsgFormat ffmt("%10.9f");
   //store ony candidate handles that are not part of CandTrackSRHandle
   TIter sitr(cslh->GetDaughterIterator());
   while( CandStripHandle* csh = dynamic_cast<CandStripHandle*>(sitr()) ){
      MSG("Align", Msg::kVerbose) << "CandStrip: " << csh->GetStripEndId() 
                                  << " time: " << ffmt(csh->GetBegTime())
                                  << " charge = " << csh->GetCharge(CalDigitType::kNone) <<endl;  
      AlignmentStrip astrip(csh);
      if(csh->GetPlaneView() == PlaneView::kV){
         bool result = StripBelongsToVTrack(astrip); 
         if(!result){
            fCandVStrip.push_back(astrip);
            fCandVCharge += astrip.charge;
            MSG("Align", Msg::kVerbose) <<"Adding candstrip " <<PlexStripEndId(astrip.plexseid)<<endl;
         }
         continue;
      }
      
      if(csh->GetPlaneView() == PlaneView::kU){
         bool result = StripBelongsToUTrack(astrip); 
         if(!result){
            fCandUStrip.push_back(astrip);
            fCandUCharge += astrip.charge;
            MSG("Align", Msg::kVerbose) <<"Adding candstrip " <<PlexStripEndId(astrip.plexseid)<<endl;
         }
         continue;
      } 
      MSG("Align", Msg::kError) << "Strip is not in V or U view!!!"<<endl;
   }
}


void DetectorAlignment::MakeAlignmentTrack(vector<AlignmentStrip>::iterator begin,
                                           vector<AlignmentStrip>::iterator end)
{
   MSG("Align", Msg::kVerbose) << "MakeAlignmentTrack()..."<<endl;

   f2dTrackCharge      = 0.0;
   f2dTrackResidualRMS = 0.0;
   fDeltaT             = 0.0;   

   vector<AlignmentStrip>::iterator itr;   
   double size = 0.0;

   vector<Double_t> begtime;
   vector<Double_t> endtime;
   
   //Filing vectors with data from 2d track
   for (itr = begin; itr != end; ++itr) {
      size += 1.0;
      const double residual = FitTrackLessOne(begin, end, itr);
      const double &charge   = itr -> charge;

      //store residual in AlignStrip
      itr -> residual = residual;

      //calculating "rms" of a current track
      f2dTrackResidualRMS += residual*residual;
      
      begtime.push_back(itr -> begtime);
      endtime.push_back(itr -> endtime);
      
      //add charge as we iterate along 2d track:
      if(charge > 0.0)
         f2dTrackCharge += itr->charge;
      else
         MSG("Align", Msg::kVerbose) <<PlexStripEndId(itr->plexseid)<<" charge = "<<charge<<endl;      
      
      MSG("Align", Msg::kVerbose) <<PlexStripEndId(itr->plexseid)<<" residual = "<<residual<<endl;
   }

   //sorting begtime and endtime to get time time interval of 2d track
   if(begtime.size() > 0 && endtime.size() > 0){
      sort(begtime.begin(), begtime.end());
      sort(endtime.begin(), endtime.end());
      fDeltaT = endtime[endtime.size()-1] - begtime[0];
      fTrackTimes.push_back(endtime[endtime.size()-1]);
      fTrackTimes.push_back(begtime[0]);
   }
   
   if(size > 0.0) f2dTrackResidualRMS = pow(f2dTrackResidualRMS/size, 0.5);

   MSG("Align", Msg::kVerbose) << "MakeAlignmentTrack() DONE!" << endl;
}


// Fit 2d track with a straight line using least squares method
double DetectorAlignment::FitTrackLessOne(vector<AlignmentStrip>::const_iterator begin,
                                          vector<AlignmentStrip>::const_iterator end,
                                          vector<AlignmentStrip>::const_iterator skip)
{
   MSG("Align", Msg::kVerbose) << "FitTrack..." << endl;

   fFita = 0.0;
   fFitb = 0.0;
   fSigmaOfa = 0.0;
   fSigmaOfb = 0.0;
   fCovab = 0.0;
   fUncertaintyInTPosofStrips = 0.0;

   // Linear regression parameters least squares method. See for eg. Num Rec Ch 15
   // Fit y = ax + b, as in W. R. Leo, Techniques for ..., Section 4.7.2
   // x is strip's Z position, y is strip's transverse position (U or V).
   double A = 0.0, B=0.0, C=0.0, D=0.0, E=0.0, F=0.0, S=0.0;
   
   for(vector<AlignmentStrip>::const_iterator run = begin; run != end; ++run)
   {
      if(run != skip)
      {
         const double &x = run -> zpos;
         const double &y = run -> tpos;
         A += x;
         B += 1.0;
         C += y;
         D += x*x;
         E += x*y;
         F += y*y;
      }
   }
   
   const double det = B*D - A*A; //determinant of error matrix
   if(det < 0.000001)
   {
      MSG("Align", Msg::kError) << "Linear fit failed!" << endl;
      return -100.0;
   }
   
   // The linear fitparameters in y = ax + b
   const double a = (E*B-C*A)/det;
   const double b = (D*C-E*A)/det;

   if(skip != end)
   {
      const double skipped_x = skip -> zpos;
      const double skipped_y = skip -> tpos;
      const double predict   = a*skipped_x + b;       
      return (predict - skipped_y);
   }
   
   fFita     = a;
   fFitb     = b;
   fSigmaOfa = pow(B/det, 0.5);
   fSigmaOfb = pow(D/det, 0.5);
   fCovab    = -A/det;   

   double npoints = 0.0;
   for(vector<AlignmentStrip>::const_iterator run = begin; run != end; ++run)
   {
      const double &x = run -> zpos;
      const double &y = run -> tpos;
      S += (y-a*x-b)*(y-a*x-b);
      npoints += 1.0;
   }
   
   if(npoints - 2.0 < 0.9){
      MSG("Align", Msg::kError) << "Linear fit failed!" << endl;
      return -100.0;
   }
   
   //square root of least squares sum
   fUncertaintyInTPosofStrips = pow(S/(npoints-2.0), 0.5);

   return fUncertaintyInTPosofStrips;
}


bool  DetectorAlignment::StripBelongsToVTrack(const AlignmentStrip& astrip)
{
   //if there is a track strip such that: PlexStripEndIds are the same
   //                                     beggining times are within 10ns of each other
   //then return true
   //
   //otherwise return false
   
   for(unsigned int i = 0; i < fTrackVStrip.size(); ++i){
      if(fTrackVStrip[i].plexseid == astrip.plexseid){
         double diff = fabs(fTrackVStrip[i].begtime - astrip.begtime);
         if(diff < 0.000000002){
            MsgFormat ffmt("%10.9f");
            MsgFormat fi("%3d");
            MSG("Align", Msg::kVerbose) << "Removing "
                                        << PlexStripEndId(astrip.plexseid) << " V track strip #"<<fi(i) <<endl
                                        << "    Charge = " << astrip.charge 
                                        << " lhs t = " << ffmt(fTrackVStrip[i].begtime)
                                        << " rhs t = " << ffmt(astrip.begtime)
                                        << " diff = "<< ffmt(diff) <<endl;
            
            return true;  
         }
      }
   }
   return false;
}

bool  DetectorAlignment::StripBelongsToUTrack(const AlignmentStrip& astrip)
{
   //if there is a track strip such that: PlexStripEndIds are the same
   //                                     beggining times are within 1ns of each other
   //then return true
   //
   //otherwise return false
   
   for(unsigned int i = 0; i < fTrackUStrip.size(); ++i){
      if(fTrackUStrip[i].plexseid == astrip.plexseid){
         double diff = fabs(fTrackUStrip[i].begtime - astrip.begtime);
         if(diff < 0.000000002){
            MsgFormat ffmt("%10.9f");
            MsgFormat fi("%3d");
            MSG("Align", Msg::kVerbose) << "Removing "
                                        << PlexStripEndId(astrip.plexseid) << " U track strip #"<<fi(i)<<endl
                                        << "    Charge = " << astrip.charge 
                                        << " lhs t = " << ffmt(fTrackUStrip[i].begtime)
                                        << " rhs t = " << ffmt(astrip.begtime) 
                                        << " diff = "<< ffmt(diff) <<endl;
            return true;  
         }   
      }
   }
   return false;
}

---