TrackDirectionModule.cxx

Go to the documentation of this file.

//
//
// TrackDirectionModule
//
// Package: DetectorAlignment
//
// Contact: rustem@fnal.gov
//
// Created on: Fri Jul 23 15:26:05 2005
//

//LOCAL package headers
#include "TrackDirectionModule.h"
#include "NtpAlignmentRecord.h"
#include "NtpAlignTrackStrip.h"

//MINOS headers
#include "JobControl/JobCModuleRegistry.h"
#include "MessageService/MsgService.h"
#include "MessageService/MsgFormat.h"
#include "MinosObjectMap/MomNavigator.h"
#include "Validity/VldContext.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "AstroUtil/AstCoordinate.h"

//ROOT headers
#include "TFile.h"
#include "TTree.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TStopwatch.h"

//C++ headers
#include <cassert>
#include <sstream>
#include <cmath>

CVSID("$Id: TrackDirectionModule.cxx,v 1.2 2005/07/30 21:25:31 rustem Exp $");

JOBMODULE(TrackDirectionModule,"TrackDirectionModule","TrackDirectionModule");

using namespace std;

TrackDirectionModule::TrackDirectionModule()
   :fTrackDirectionFile(0),
    fRun(0),
    fSubRun(0),
    fhCosX(0),
    fhCosY(0),
    fhCosZ(0),
    fhTrackLength(0),
    f2hCosXvsCosZ(0),
    f2hCosYvsCosZ(0),
    f2hCosYvsCosX(0),
    fMaxTrackChargeCut(60000.0),       //cut on max 2d track charge used for alignment
    fMinTrackChargeCut(4000.0),        //cut on min 2d track charge used for alignment
    fMaxSigmaOfTPosCut(0.03),          //cut on uncertainty in tpos from 2d track fit
    fTrackChargeRatioCut(0.5),         //min 3d track's charge ratio to total snarl charge
    fMinCosZCut(0.25),                 //min fabs(cosz) of 3d track angle with Z axis
    fMinTrackStripChargeCut(100.0),    //min charge in a alternative track strip 
    fMaxTrackStripChargeCut(5000.0),   //max charge in a hit alignment strip
    fApplyCuts(false),
    fNRecords(0),
    fNFailedRead(0),
    fFailedCut(0),
    fNFailedCutMaxVStripCharge(0),
    fNFailedCutMaxUStripCharge(0),
    fNFailedTrackChargeCut(0),
    fNFailedCosZCut(0),
    fNFailedChargeRatioCut(0),
    fNFailedSigmaTPosCut(0)
{
}

void TrackDirectionModule::BeginJob()
{ 
   MSG("TrackDir", Msg::kDebug) << "TrackDirectionModule::BeginJob()"<< endl;
   MSG("TrackDir", Msg::kInfo) << "Initializing variables and histograms..." << endl;

   string OutfilePath;
   const char* strENV_OUTFILE_DIR = getenv("ENV_OUTFILE_DIR");
   if ( strENV_OUTFILE_DIR == 0 || strlen(strENV_OUTFILE_DIR) == 0 )
   {
      MSG("TrackDir", Msg::kFatal) << "ENV_OUTFILE_DIR is not set!"
                                << " Program is terminating now." <<endl;
      abort();
   } else 
      OutfilePath = string(strENV_OUTFILE_DIR);
   
   const char* strENV_APPLY_TRACK_CUTS = getenv("ENV_APPLY_TRACK_CUTS");
   if ( strENV_APPLY_TRACK_CUTS == 0 || strlen(strENV_APPLY_TRACK_CUTS) == 0 )
   {
      MSG("TrackDir", Msg::kWarning) << "ENV_APPLY_TRACK_CUTS is not set!"
                                  << " Using default value: fApplyCuts=" <<fApplyCuts<<endl;
   } else 
      if(atoi(strENV_APPLY_TRACK_CUTS) == 1)
         fApplyCuts = true;
      else
         fApplyCuts = false;
   
   
   MSG("TrackDir", Msg::kInfo) << "OutfilePath = " << OutfilePath << endl
                            << "fApplyCuts = " << fApplyCuts << endl << endl;
   
   stringstream ss;
   ss << OutfilePath << "/track_direction.root";

   fTrackDirectionFile = new TFile(ss.str().c_str(), "RECREATE");
   fTrackDirectionTree = new TTree("TrackDirection", "TrackDirection");

   fTrackDirectionTree -> Branch("Run", &fRun, "fRun/D"); 
   fTrackDirectionTree -> Branch("SubRun", &fSubRun, "fSubRun/D"); 
   fTrackDirectionTree -> Branch("CosX", &fCosX, "fCosX/D"); 
   fTrackDirectionTree -> Branch("CosY", &fCosY, "fCosY/D"); 
   fTrackDirectionTree -> Branch("CosZ", &fCosZ, "fCosZ/D"); 
   fTrackDirectionTree -> Branch("Length", &fLength, "fLength/D"); 
   fTrackDirectionTree -> Branch("Sigma", &fSigma, "fSigma/D"); 
   fTrackDirectionTree -> Branch("Time", &fTime, "fTime/D"); 
   fTrackDirectionTree -> Branch("Altitude", &fAltitude, "fAltitude/D"); 
   fTrackDirectionTree -> Branch("Azimuth", &fAzimuth, "fAzimuth/D"); 
   fTrackDirectionTree -> Branch("Charge", &fCharge, "fCharge/D"); 
   fTrackDirectionTree -> Branch("NStrip", &fNStrip, "fNStrip/D"); 
   fTrackDirectionTree -> Branch("Year", &fYear, "fYear/i"); 
   fTrackDirectionTree -> Branch("Month", &fMonth, "fMonth/i"); 
   fTrackDirectionTree -> Branch("Day", &fDay, "fDay/i"); 
   fTrackDirectionTree -> Branch("Hour", &fHour, "fHour/i"); 
   fTrackDirectionTree -> Branch("Minute", &fMinute, "fMinute/i"); 
   fTrackDirectionTree -> Branch("Sec", &fSec, "fSec/i"); 

   fhCosX = new TH1D("CosX", "Directional cosine with x axis", 100, -1.0, 1.0);
   fhCosY = new TH1D("CosY", "Directional cosine with y axis", 100,  0.0, 1.0);
   fhCosZ = new TH1D("CosZ", "Directional cosine with z axis", 100, -1.0, 1.0);
   fhTrackLength = new TH1D("TrackLength", "Track length in meter", 500, 0.0, 16.0);
   
   f2hCosXvsCosZ = new TH2D("CosXvsCosZ", "CosX vs CosZ", 100, -1.0, 1.0, 100, -1.0, 1.0);
   f2hCosYvsCosZ = new TH2D("CosYvsCosZ", "CosY vs CosZ", 100, -1.0, 1.0, 100,  0.0, 1.0);
   f2hCosYvsCosX = new TH2D("CosYvsCosX", "CosY vs CosX", 100, -1.0, 1.0, 100,  0.0, 1.0);   

   fTimer.Start();
   fTimerInterval.Start();
}


JobCResult TrackDirectionModule::Ana(const MomNavigator* mom)
{

   NtpAlignmentRecord* ntprec =
      dynamic_cast<NtpAlignmentRecord*>(mom->GetFragment("NtpAlignmentRecord"));
   
   if(!ntprec)
      return JobCResult::kAOK;   

   const RecHeader hdr = ntprec->GetHeader();
   const VldContext vldc = hdr.GetVldContext();
   fRun = (ntprec -> GetHeader()).GetRun();
   fSubRun = (ntprec -> GetHeader()).GetSubRun();
   if(fNRecords % 10000 == 0){
      MsgFormat f("%10d");      
      MSG("TrackDir", Msg::kInfo) << "Found NtpAlignmentRecord #" << f(fNRecords)
                                  << " from run #" << fRun << ", subrun #" << fSubRun
                                  << ", cpu time = " << fTimerInterval.RealTime() <<endl;
      fTimerInterval.Reset();
      fTimerInterval.Start();
   }   
   
   fNRecords++;      
   
   if(!GetTrackRecoQuality(ntprec)){
      fFailedCut++;
      return JobCResult::kAOK;   
   }      
      
   //Read residuals and other strips' info to vectors of AlignmentStrips
   //This is the main function where lots of work is done
   if(!ProcessRecord(ntprec, vldc)){
      fFailedCut++;
      return JobCResult::kAOK;      
   }   
   
   if(GetDirectionalCosines(fCosX, fCosY, fCosZ, fLength, fSigma, vldc)){
      fTime = vldc.GetTimeStamp().GetSec();
      vldc.GetTimeStamp().GetDate(kTRUE, 0, &fYear, &fMonth, &fDay);
      vldc.GetTimeStamp().GetTime(kTRUE, 0, &fHour, &fMinute, &fSec);
      fCharge = ntprec -> ucharge + ntprec -> vcharge;
      fNStrip = ntprec -> ntrackustrip + ntprec -> ntrackvstrip;
      AstUtil::LocalToHorizon(fCosX, fCosY, fCosZ, vldc.GetDetector(), fAltitude, fAzimuth);
      fTrackDirectionTree -> Fill();
      
      fhCosX -> Fill(fCosX);
      fhCosY -> Fill(fCosY);
      fhCosZ -> Fill(fCosZ);
      fhTrackLength -> Fill(fLength);
      
      f2hCosXvsCosZ -> Fill(fCosZ, fCosX);
      f2hCosYvsCosZ -> Fill(fCosZ, fCosY);
      f2hCosYvsCosX -> Fill(fCosX, fCosY);
   }
   
   return JobCResult::kAOK;   
}

void TrackDirectionModule::EndJob()
{ 
   PrintJobStatistics();   

   fTrackDirectionFile -> Write();
   fTrackDirectionFile -> Close();

   MSG("TrackDir", Msg::kInfo) << "TrackDirectionModule is done." <<endl;
}

bool TrackDirectionModule::GetTrackRecoQuality(const NtpAlignmentRecord* ntprec)
{
   //Check that: 1) Ratio of track charge to total charge in the current 
   //            record is greater than fTrackChargeRatioCut
   //            2) Uncertainty in tpos of a fit is less that fSigmaOfTPosCut
   //            3) Charge of 2d tracks is less that fTrackChargeCut

   if(!fApplyCuts)
      return true;

   //Cut on 2d track charge
   if(ntprec->vcharge > fMaxTrackChargeCut || ntprec->vcharge < fMinTrackChargeCut ||
      ntprec->ucharge > fMaxTrackChargeCut || ntprec->ucharge < fMinTrackChargeCut)
   {
      fNFailedTrackChargeCut++;
      return false;
   }

   //Cut in tpos uncertainty from linear fit to a track
   if(ntprec->vsigmaoftpos > fMaxSigmaOfTPosCut || ntprec->usigmaoftpos > fMaxSigmaOfTPosCut)
   {
      fNFailedSigmaTPosCut++;
      return false;
   }

   //Cut on 3d track Hough cos with horizontal z axis
   if(fabs(ntprec->hcosz) < fMinCosZCut)
   {
      fNFailedCosZCut++;
      return false;
   }

   const double allcharge = ntprec -> vcharge + ntprec -> ucharge
      + ntprec -> vcandcharge + ntprec -> ucandcharge;

   if(allcharge < 1.0)
      return false;
   
   Double_t ratio = (ntprec -> vcharge + ntprec -> ucharge)/allcharge;
   
   if(ratio < fTrackChargeRatioCut)
   {
      fNFailedChargeRatioCut++;
      return false;
   }   

   return true;
}


bool TrackDirectionModule::ProcessRecord(const NtpAlignmentRecord* ntprec, 
                                        const VldContext& vldc)
{
   fTrackVStrip.clear();
   fTrackUStrip.clear();
   
   UgliGeomHandle ugh(vldc);   

   for(unsigned int i = 0; i < ntprec -> ntrackvstrip; ++i){
      NtpAlignTrackStrip* strip = dynamic_cast<NtpAlignTrackStrip*>(ntprec->trackvstrip->At(i));
      assert(strip);
      
      AlignmentStrip astrip(strip);
      if(fApplyCuts && astrip.charge > fMaxTrackStripChargeCut){
         fNFailedCutMaxVStripCharge++;
         return false;
      }
      
      PlexStripEndId plexid(astrip.plexseid);
      UgliStripHandle ush = ugh.GetStripHandle(plexid);

      astrip.plane = plexid.GetPlane();
      astrip.strip = plexid.GetStrip();
      astrip.tpos = ush.GetTPos();
      astrip.zpos = ush.GetScintPlnHandle().GetZ0();
      astrip.tposrelmdl = ush.GetTPosRelMdl();
      astrip.lposrelmdl = ush.GetLPosRelMdl();

      //use fit results in opposite view to get longitudinal position
      astrip.ghitpos = (ntprec -> ufita)*astrip.zpos + ntprec -> ufitb;      

      ConvertToLocal(astrip, plexid, ugh, ush);
      fTrackVStrip.push_back(astrip);
   }


   for(unsigned int i = 0; i < ntprec -> ntrackustrip; ++i){
      NtpAlignTrackStrip* strip = dynamic_cast<NtpAlignTrackStrip*>(ntprec->trackustrip->At(i));
      assert(strip);
      
      AlignmentStrip astrip(strip);      
      if(fApplyCuts && astrip.charge > fMaxTrackStripChargeCut){
         fNFailedCutMaxUStripCharge++;
         return false;
      }

      PlexStripEndId plexid(astrip.plexseid);
      UgliStripHandle ush = ugh.GetStripHandle(plexid);

      astrip.plane = plexid.GetPlane();
      astrip.strip = plexid.GetStrip(); 
      astrip.tpos = ush.GetTPos();
      astrip.zpos = ush.GetScintPlnHandle().GetZ0();
      astrip.tposrelmdl = ush.GetTPosRelMdl();
      astrip.lposrelmdl = ush.GetLPosRelMdl();

      //use fit results in opposite view to get longitudinal position
      astrip.ghitpos = (ntprec -> vfita)*astrip.zpos + ntprec -> vfitb;
      
      ConvertToLocal(astrip, plexid, ugh, ush);
      fTrackUStrip.push_back(astrip);
   }

   return true;
}


bool TrackDirectionModule::GetDirectionalCosines(Double_t &cosx, Double_t &cosy, Double_t &cosz, 
                                                 Double_t &length, Double_t &sigma, const VldContext &vld)
{
   UgliGeomHandle ugh(vld);
   map<int, TVector3> track;
   for(unsigned int i = 0; i < fTrackVStrip.size(); ++i){
      AlignmentStrip &astrip = fTrackVStrip[i];
      TVector3 uvz(astrip.ghitpos, astrip.tpos, astrip.zpos);
      TVector3 xyz = ugh.uvz2xyz(uvz);
      pair <int, TVector3> p(astrip.plane, xyz);
      track.insert(p);
   }
   
   for(unsigned int i = 0; i < fTrackUStrip.size(); ++i){
      AlignmentStrip &astrip = fTrackUStrip[i];
      TVector3 uvz(astrip.tpos, astrip.ghitpos, astrip.zpos);
      TVector3 xyz = ugh.uvz2xyz(uvz);
      pair <int, TVector3> p(astrip.plane, xyz);
      track.insert(p);
   }
   
   Double_t A=0.0, B=0.0, D=0.0;
   Double_t Cx=0.0, Ex=0.0, Fx=0.0;
   Double_t Cy=0.0, Ey=0.0, Fy=0.0;

   int firstplane=1000, lastplane=0;
   for(map<int, TVector3>::const_iterator it = track.begin(); it != track.end(); ++it){
      const TVector3 &v = it->second;
      const Double_t &x = v.x();
      const Double_t &y = v.y();
      const Double_t &z = v.z();
      A += z, B += 1.0, D += z*z;
      Cx += x, Ex += z*x, Fx += x*x;
      Cy += y, Ey += z*y, Fy += y*y;
      const int &plane = it -> first;
      if(firstplane > plane)
         firstplane = plane;
      if(lastplane < plane)
         lastplane = plane;
   }

   Double_t det = B*D - A*A;

   if(fabs(det) < 0.000001){
      MSG("TrackDir", Msg::kWarning) << "Linear fit failed!" << endl;
      return false;
   }
   
   Double_t FitAx = (Ex*B-Cx*A)/det;
   Double_t FitAy = (Ey*B-Cy*A)/det;
   Double_t FitBx = (D*Cx-Ex*A)/det;
   Double_t FitBy = (D*Cy-Ey*A)/det;

   if(FitAy < 0.0)
      cosz = - pow(1.0/(FitAx*FitAx + FitAy*FitAy + 1.0), 0.5);
   else 
      cosz = pow(1.0/(FitAx*FitAx + FitAy*FitAy + 1.0), 0.5);

   cosy = fabs(FitAy*cosz);
   cosx = FitAx*fabs(cosz);
   
   const TVector3 &fs = track[firstplane];
   const TVector3 &ls = track[lastplane];
   TVector3 tl(fs.x()-ls.x(), fs.y()-ls.y(), fs.z()-ls.z());
   length = tl.Mag();

   Double_t sumx=0.0, sumy=0.0, npoints=0.0;
   for(map<int, TVector3>::const_iterator it = track.begin(); it != track.end(); ++it){
      const TVector3 &v = it->second;
      const Double_t &x = v.x();
      const Double_t &y = v.y();
      const Double_t &z = v.z();
      sumx += (x-FitAx*z-FitBx)*(x-FitAx*z-FitBx);
      sumy += (y-FitAy*z-FitBy)*(y-FitAy*z-FitBy);
      npoints += 1.0;
   }
   
   npoints = npoints - 2.0;
   if(npoints < 0.9){
      MSG("TrackDir", Msg::kWarning) << "Linear fit failed!" << endl;
      return false;
   }
   sigma = pow(sumx/npoints, 0.5) + pow(sumy/npoints, 0.5);
   return true;
}


void TrackDirectionModule::ConvertToLocal(AlignmentStrip &astrip,
                                         PlexStripEndId &plexid, 
                                         UgliGeomHandle &ugh,
                                         UgliStripHandle &ush)
{
   // calculate the position of strip ends in (strip) local co-ordinates
   // x = distance along the strip, y=distance in +tpos coordinate
   // the "logical" strip is the one as-if there were no coil/bypass
   Double_t halflen = ush.GetHalfLength();
   //   TVector3 localCenter(0,0,0);  // center of the "logical" strip
   //   TVector3 localEastEnd(-halflen,0,0);
   TVector3 localWestEnd(+halflen, 0.0, 0.0);
   //   TVector3 localHitLPos(astrip.lhitpos,0,0);
   
   // calculate the positions in global XYZ space
   //TVector3 xyz0 = ush.LocalToGlobal(localCenter);
   //TVector3 xyzE = ush.LocalToGlobal(localEastEnd);
   //TVector3 xyzW = ush.LocalToGlobal(localWestEnd);

   Double_t u = 0.0, v = 0.0;

   if(plexid.GetPlaneView() == PlaneView::kU){
      v = astrip.ghitpos;
      u = astrip.residual + astrip.tpos;
   } else {
      u = astrip.ghitpos;
      v = astrip.residual + astrip.tpos;
   }
   

   TVector3 ghituvz(u, v, astrip.zpos);
   TVector3 ghitxyz = ugh.uvz2xyz(ghituvz);
   TVector3 lhitxyz = ush.GlobalToLocal(ghitxyz);
      
   astrip.length = 2.0*halflen;
   astrip.lhitpos = lhitxyz.x();
   astrip.pigtail = ush.WlsPigtail(StripEnd::kWest);
   astrip.wlsbypass = ush.WlsBypass();

   if(fabs(lhitxyz.x()) > halflen)
      astrip.goodhit = false;

   return;
/*   
     double y = lhitxyz.y();
     double z = lhitxyz.z();
     fHistogram -> FillRoundingErrors(plexid, y, z);
     if(fabs(y + z) > 1.0e-4)
     MSG("TrackDir", Msg::kWarning) <<"In local coordinates y + z = " << fabs(y + z) << endl;
*/
}


void TrackDirectionModule::PrintJobStatistics()
{
   MSG("TrackDir", Msg::kInfo) << "TrackDirectionModule::EndJob()"<< endl << endl;    

   MSG("TrackDir", Msg::kInfo) << "Cpu time = " << fTimer.CpuTime() <<endl;
   MSG("TrackDir", Msg::kInfo) << "Real time = " << fTimer.RealTime() <<endl;
   MSG("TrackDir", Msg::kInfo) << "fNRecords  = " << fNRecords << endl;
   MSG("TrackDir", Msg::kInfo) << "fFailedCut = " << fFailedCut << endl;
   MSG("TrackDir", Msg::kInfo) << "fNFailedRead = " << fNFailedRead << endl;
   MSG("TrackDir", Msg::kInfo) << "fNFailedCutMaxVStripCharge = " << fNFailedCutMaxVStripCharge <<endl;
   MSG("TrackDir", Msg::kInfo) << "fNFailedCutMaxUStripCharge = " << fNFailedCutMaxUStripCharge <<endl;
   MSG("TrackDir", Msg::kInfo) << "fNFailedTrackChargeCut  = " << fNFailedTrackChargeCut << endl;
   MSG("TrackDir", Msg::kInfo) << "fNFailedCosZCut  = "        << fNFailedCosZCut << endl;
   MSG("TrackDir", Msg::kInfo) << "fNFailedSigmaTPosCut  = "   << fNFailedSigmaTPosCut << endl;
   MSG("TrackDir", Msg::kInfo) << "fNFailedChargeRatioCut  = " << fNFailedChargeRatioCut << endl;
}

---