MakeAlignmentModule.cxx

Go to the documentation of this file.

//
//
// MakeAlignmentModule
//
// Package: DetectorAlignment
//
// Alignment engine that processes ntuples files
//
// Contact: rustem@fnal.gov
//
// Created on: Fri Sep 24 15:26:05 2004
//

//Local package headers
#include "MakeAlignmentModule.h"
#include "NtpAlignmentRecord.h"
#include "NtpAlignTrackStrip.h"
#include "NtpAlignCandStrip.h"

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

//ROOT headers
#include "TFile.h"

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


CVSID("$Id: MakeAlignmentModule.cxx,v 1.6 2005/08/10 21:07:39 rustem Exp $");

JOBMODULE(MakeAlignmentModule,"MakeAlignmentModule","MakeAlignmentModule");

using namespace std;

MakeAlignmentModule::MakeAlignmentModule()
   :fAlgorithm(0),
    fHistogram(0),
    fNearCoilOn(2004, 12, 9, 0, 0, 0), //Time when ND coil was first turned on
    fIteration(0),                      //Alignment iteration number
    fCurrentRun(0),
    fCurrentSubRun(0),
    fNProcessedSubruns(0),
    fNSubRunsPerSubset(41),
    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
    fSpecialNDRunNumber(5845),          //ND run # used in 1st ND alignment with magnetic field
    fRecalculateResidual(false),
    fApplyCuts(false),
    fStatisticalSubsetSize(1000000),    //Number of records in a statistical subset
    fNFailedFit(0),
    fNRecords(0),
    fNPassedRecords(0),
    fNSubsetRecords(0),
    fNFailedRead(0),
    fNMissedVPlane(0),
    fNMissedUPlane(0),
    fTotalVPlaneHits(0),
    fTotalUPlaneHits(0),
    fFailedCut(0),
    fNHitsOutsideStrip(0),
    fNFailedCutMaxVStripCharge(0),
    fNFailedCutMaxUStripCharge(0),
    fMissedVPlaneMeanCharge(0.0),
    fMissedUPlaneMeanCharge(0.0),
    fMissedVPlaneChargeSigma(0.0),
    fMissedUPlaneChargeSigma(0.0),
    fNFailedTrackChargeCut(0),
    fNFailedStripChargeCut(0),
    fNFailedCosZCut(0),
    fNFailedChargeRatioCut(0),
    fNFailedSigmaTPosCut(0)
{
}

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

   const char* strENV_ALIGN_ITER = getenv("ENV_ALIGN_ITER");
   if ( strENV_ALIGN_ITER == 0 || strlen(strENV_ALIGN_ITER) == 0 )
   {
      MSG("Align", Msg::kFatal) << "ENV_ALIGN_ITER is not set!" <<endl
                                << " Program is terminating now." <<endl;
      abort();
   } else 
      fIteration = atoi(strENV_ALIGN_ITER);
   
   string alignmentpath;
   const char* strENV_ALIGN_DIR = getenv("ENV_ALIGN_DIR");
   if ( strENV_ALIGN_DIR == 0 || strlen(strENV_ALIGN_DIR) == 0 )
   {
      MSG("Align", Msg::kFatal) << "ENV_ALIGN_DIR is not set!"
                                << " Program is terminating now." <<endl;
      abort();
   } else
      alignmentpath = string(strENV_ALIGN_DIR);
   
   const char* strENV_ALIGN_APPLYCUTS = getenv("ENV_ALIGN_APPLYCUTS");
   if ( strENV_ALIGN_APPLYCUTS == 0 || strlen(strENV_ALIGN_APPLYCUTS) == 0 )
   {
      MSG("Align", Msg::kWarning) << "ENV_ALIGN_APPLYCUTS is not set!"
                                  << " Using default value: fApplyCuts=false" <<endl;
   } else 
      if(atoi(strENV_ALIGN_APPLYCUTS) == 1)
         fApplyCuts = true;
      
   const char* strENV_ALIGN_RECALCULATE = getenv("ENV_ALIGN_RECALCULATE");
   if ( strENV_ALIGN_RECALCULATE == 0 || strlen(strENV_ALIGN_RECALCULATE) == 0 )
   {
      MSG("Align", Msg::kWarning) << "ENV_ALIGN_RECALCULATE is not set!"
                                  << " Using default value: fRecalculateResidual=false" <<endl;
   } else 
      if(atoi(strENV_ALIGN_RECALCULATE)==1)
         fRecalculateResidual = true;

   MSG("Align", Msg::kInfo) << endl << "fIteration = " << fIteration << endl
                            << "Alignment path: " << alignmentpath << endl
                            << "fApplyCuts = " << fApplyCuts <<endl
                            << "fRecalculateResidual = " << fRecalculateResidual <<endl << endl;
   
   stringstream ss;
   if(fApplyCuts)
      ss << alignmentpath << "/new_alignment_" << fIteration << ".root";
   else
      ss << alignmentpath << "/new_alignment_no_cuts_" << fIteration << ".root";         
   
   fRootFile = new TFile(ss.str().c_str(),"RECREATE");
   
   fTimer.Start();
   fTimerInterval.Start();
}


JobCResult MakeAlignmentModule::Ana(const MomNavigator* mom)
{
   fNRecords++;        

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

   const RecCandHeader candheader = ntprec->GetHeader();
   const int run    = candheader.GetRun();
   const int subrun = candheader.GetSubRun();
   VldContext vldc  = candheader.GetVldContext();

   if(fNRecords % 10000 == 0)
   {
      MsgFormat f("%10d");      
      MSG("Align", Msg::kInfo) << "Found NtpAlignmentRecord #" << f(fNRecords)
                               << " run #" << run << ", subrun #" << subrun
                               << ", cpu time = " << fTimerInterval.RealTime() <<endl;
      fTimerInterval.Reset();
      fTimerInterval.Start();
   }    

   if(!fAlgorithm && !fHistogram)
   { 
      fValidityContext = vldc;
      fHistogram = new AlignmentHistograms(fRootFile, fValidityContext);
      fAlgorithm = new AlignmentAlgorithm(fRootFile, fValidityContext);
   }
   
   if(fValidityContext.GetDetector() != vldc.GetDetector())
   {
      MSG("Align", Msg::kError) << "The data in current record is for the different detector!" << endl;
      return JobCResult::kAOK;
   } else
      fValidityContext = vldc;   

   //Fill run summary object
   fRunSummary.Fill(ntprec); 

   if(fApplyCuts && !GetTrackQuality(ntprec))
   {
      fFailedCut++;
      return JobCResult::kAOK;   
   }      
      
   //Read residuals and other strips' info to vectors of AlignmentStrips
   if(!ReadRecord(ntprec))
   {
      fFailedCut++;
      fNFailedRead++;
      return JobCResult::kAOK;      
   }   

   if(fRecalculateResidual && !RecalculateResiduals()) 
   {
      return JobCResult::kAOK;
   }
   
   ProcessRecord(ntprec);

   fHistogram -> FillPlanePlex(fTrackVStrip, fCandVStrip);
   fHistogram -> FillPlanePlex(fTrackUStrip, fCandUStrip);

   fNPassedRecords++;
   fNSubsetRecords++;
   
   if(fCurrentRun != run || fCurrentSubRun != subrun)
   {
      fCurrentRun    = run;
      fCurrentSubRun = subrun;
      fNProcessedSubruns++;      
      if(fNProcessedSubruns % fNSubRunsPerSubset == 0)
      {
         if(fCurrentRun != fSpecialNDRunNumber ||
            fValidityContext.GetDetector() != DetectorType::kNear)
         {
            fAlgorithm -> ResetSubset(fValidityContext, fNSubsetRecords);
            fNSubsetRecords = 0;
         }      
      }
      if(fCurrentRun == fSpecialNDRunNumber &&
         fValidityContext.GetDetector() == DetectorType::kNear)
         if(fCurrentSubRun == 0 || fCurrentSubRun == 24)
         {
            fAlgorithm -> ResetSubset(fValidityContext, fNSubsetRecords);
            fNSubsetRecords = 0;
         }
   }
   
   return JobCResult::kAOK;   
}

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

   if(fHistogram) delete fHistogram;

   if(fAlgorithm)
   {
      fAlgorithm -> ResetSubset(fValidityContext, fNSubsetRecords);
      fAlgorithm -> RunAlgorithm();
      delete fAlgorithm;
   }

   fRunSummary.FillSummaryTree(fRootFile);

   MSG("Align", Msg::kInfo) << "Saving root file " << fRootFile -> GetName() << endl;
   fRootFile -> Write();
   fRootFile -> Close();

   MSG("Align", Msg::kInfo) << "MakeAlignmentModule is done." << endl;
}


//__________________________________________________________________________________________________
bool MakeAlignmentModule::GetTrackQuality(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

   //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 MakeAlignmentModule::ReadRecord(const NtpAlignmentRecord* ntprec)
{
   fTrackVStrip.clear();
   fTrackUStrip.clear();
   fCandVStrip.clear();
   fCandUStrip.clear();

   UgliGeomHandle ugh(fValidityContext);

   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++;
         fNFailedStripChargeCut++;
         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();
      astrip.charge_pc  = astrip.charge * ntprec -> hcosz;

      //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++;
         fNFailedStripChargeCut++;
         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();
      astrip.charge_pc  = astrip.charge * ntprec -> hcosz;

      //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);
   }


   for(unsigned int i = 0; i < ntprec -> ncandvstrip; ++i)
   {
      NtpAlignCandStrip* strip = dynamic_cast<NtpAlignCandStrip*>(ntprec->candvstrip->At(i));
      assert(strip);
      
      AlignmentStrip astrip(strip);
      PlexStripEndId plexid(astrip.plexseid);
      UgliStripHandle ush = ugh.GetStripHandle(plexid);

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

      fAlgorithm -> FillCandStrip(astrip);
   }

   for(unsigned int i = 0; i < ntprec -> ncandustrip; ++i)
   {
      NtpAlignCandStrip* strip = dynamic_cast<NtpAlignCandStrip*>(ntprec->candustrip->At(i));
      assert(strip);

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

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

      fAlgorithm -> FillCandStrip(astrip);
   }

   fTotalVPlaneHits += fTrackVStrip.size();
   fTotalUPlaneHits += fTrackUStrip.size();
   return true;
}


//__________________________________________________________________________________________________
void MakeAlignmentModule::ProcessRecord(const NtpAlignmentRecord* ntprec)
{
   for(vector<AlignmentStrip>::const_iterator vit = fTrackVStrip.begin(); 
       vit != fTrackVStrip.end(); ++vit)
   {
      const AlignmentStrip &astrip = *vit;
      fAlgorithm -> Fill(astrip);
      fHistogram -> Fill(astrip, ntprec);
   }

   for(vector<AlignmentStrip>::const_iterator uit = fTrackUStrip.begin(); 
       uit != fTrackUStrip.end(); ++uit)
   {
      const AlignmentStrip &astrip = *uit;
      fAlgorithm -> Fill(astrip);
      fHistogram -> Fill(astrip, ntprec);
   }

}

//__________________________________________________________________________________________________
bool MakeAlignmentModule::RecalculateResiduals()
{
   
   UgliGeomHandle ugh(fValidityContext);

   for(unsigned int i = 0; i < fTrackVStrip.size(); ++i){
      AlignmentStrip &astrip = fTrackVStrip[i];
      PlexStripEndId plexid(astrip.plexseid);
      UgliStripHandle ush = ugh.GetStripHandle(plexid);    
      TVector3 local(astrip.lhitpos, 0, 0);
      TVector3 globalxyz = ush.GlobalPos(astrip.lhitpos);
      TVector3 globaluvz = ugh.xyz2uvz(globalxyz);
      double tpos = globaluvz.y();
      double zpos = globaluvz.z();
      if(fabs(astrip.tpos-tpos) > 0.001){
         MsgFormat f("%10.9f");
         MSG("Align",Msg::kError) << "global tpos - local tpos = " << f(astrip.tpos-tpos) <<endl;
      }
      if(fabs(astrip.zpos-zpos) > 0.001){
         MsgFormat f("%10.9f");
         MSG("Align",Msg::kError) << "zpos before - zpos after = " << f(astrip.zpos-zpos) <<endl;
      }
      astrip.tpos = globaluvz.y();
   }
   
   for(unsigned int i = 0; i < fTrackUStrip.size(); ++i){
      AlignmentStrip &astrip = fTrackUStrip[i];
      PlexStripEndId plexid(astrip.plexseid);
      UgliStripHandle ush = ugh.GetStripHandle(plexid);      
      TVector3 local(astrip.lhitpos, 0, 0);
      TVector3 globalxyz = ush.GlobalPos(astrip.lhitpos);
      TVector3 globaluvz = ugh.xyz2uvz(globalxyz);
      double tpos = globaluvz.x();
      double zpos = globaluvz.z();
      if(fabs(astrip.tpos-tpos) > 0.001){
         MsgFormat f("%10.9f");
         MSG("Align",Msg::kError) << "global tpos - local tpos = " << f(astrip.tpos-tpos) <<endl;
      }
      if(fabs(astrip.zpos-zpos) > 0.001){
         MsgFormat f("%10.9f");
         MSG("Align",Msg::kError) << "zpos before - zpos after = " << f(astrip.zpos-zpos) <<endl;
      }

      astrip.tpos = globaluvz.x();
   }

   for(vector<AlignmentStrip>::iterator vit = fTrackVStrip.begin(); 
       vit != fTrackVStrip.end(); ++vit){
      double residual = FitTrackLessOne(fTrackVStrip.begin(),   fTrackVStrip.end(), vit);
      if(fabs(residual) > 1.0)
         return false;
      AlignmentStrip &astrip = *vit;
      if(fabs(astrip.residual-residual) > 0.00001){
         MsgFormat ffrt("%7.6f");
         MSG("Align", Msg::kError) << "V residual diff after tpos correction for " 
                                   << PlexStripEndId(astrip.plexseid) << " = "
                                   << ffrt(residual - astrip.residual) << endl;
      }
      astrip.residual = residual;
   }
   
   for(vector<AlignmentStrip>::iterator uit = fTrackUStrip.begin(); 
       uit != fTrackUStrip.end(); ++uit){
      double residual = FitTrackLessOne(fTrackUStrip.begin(),   fTrackUStrip.end(), uit);
      if(fabs(residual) > 1.0)
         return false;
      AlignmentStrip &astrip = *uit;
      if(fabs(astrip.residual-residual) > 0.00001){
         MsgFormat ffrt("%7.6f");
         MSG("Align", Msg::kError) << "U residual diff after tpos correction for " 
                                   << PlexStripEndId(astrip.plexseid) << " = "
                                   << ffrt(residual - astrip.residual) << endl;
      }

      astrip.residual = residual;
   }
   return true;
}

// Fits a linear track to the straigh track's hits between "run" and
// "end" while skipping the hit "skip".  
//__________________________________________________________________________________________________
double MakeAlignmentModule::FitTrackLessOne(vector<AlignmentStrip>::const_iterator begin,
                                            vector<AlignmentStrip>::const_iterator end, 
                                            vector<AlignmentStrip>::const_iterator skip)
{

   
   MSG("Align", Msg::kVerbose) << "FitTrackLessOne..." << endl;
   // 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;

   for (vector<AlignmentStrip>::const_iterator run = begin; run != end; ++run) 
      if (run != skip){
         double x = run -> zpos;
         double y = run -> tpos;
         A += x;
         B += 1.0;
         C += y;
         D += x*x;
         E += x*y;
         F += y*y;
      }
   
   double det = B*D - A*A; //determinant of error matrix
   if(fabs(det) < 0.000001){
      MSG("Align", Msg::kWarning) << "Linear fit failed!" << endl;
      return -100.0;
   }
   
   // The linear fit: y = ax + b
   double a  = (E*B-C*A)/det;
   double b  = (D*C-E*A)/det;

   double skipped_x = skip -> zpos;
   double skipped_y = skip -> tpos;
   double predict   = a*skipped_x + b;       
    
   MSG("Align", Msg::kVerbose) << "FitTrackLessOne Done!" << endl;
   return (predict - skipped_y);
    
}

//__________________________________________________________________________________________________
void MakeAlignmentModule::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;
      fNHitsOutsideStrip++;
   }
}


//__________________________________________________________________________________________________
void MakeAlignmentModule::PrintJobStatistics()
{
   MSG("Align", Msg::kInfo) << "MakeAlignmentModule::EndJob()"<< endl << endl;    

   MSG("Align", Msg::kInfo) << "Cpu time = " << fTimer.CpuTime() <<endl;
   MSG("Align", Msg::kInfo) << "Real time = " << fTimer.RealTime() <<endl;
   MSG("Align", Msg::kInfo) << "fNPassedRecords  = " << fNPassedRecords << endl;
   MSG("Align", Msg::kInfo) << "fNRecords  = " << fNRecords << endl;
   MSG("Align", Msg::kInfo) << "fFailedCut = " << fFailedCut << endl;
   MSG("Align", Msg::kInfo) << "fNFailedFit = " << fNFailedFit << endl;
   MSG("Align", Msg::kInfo) << "fNFailedRead = " << fNFailedRead << endl;
   MSG("Align", Msg::kInfo) << "fNMissedVPlane = " << fNMissedVPlane << endl;
   MSG("Align", Msg::kInfo) << "fNMissedUPlane = " << fNMissedUPlane << endl;
   MSG("Align", Msg::kInfo) << "fTotalVPlaneHits = " << fTotalVPlaneHits << endl;
   MSG("Align", Msg::kInfo) << "fTotalUPlaneHits = " << fTotalUPlaneHits << endl;
   MSG("Align", Msg::kInfo) << "fNHitsOutsideStrip = "      << fNHitsOutsideStrip << endl;
   if(fApplyCuts)
   {
      MSG("Align", Msg::kInfo) << "fNFailedCutMaxVStripCharge = " << fNFailedCutMaxVStripCharge <<endl;
      MSG("Align", Msg::kInfo) << "fNFailedCutMaxUStripCharge = " << fNFailedCutMaxUStripCharge <<endl;
      MSG("Align", Msg::kInfo) << "fNFailedTrackChargeCut  = " << fNFailedTrackChargeCut << endl;
      MSG("Align", Msg::kInfo) << "fNFailedCosZCut  = "        << fNFailedCosZCut << endl;
      MSG("Align", Msg::kInfo) << "fNFailedSigmaTPosCut  = "   << fNFailedSigmaTPosCut << endl;
      MSG("Align", Msg::kInfo) << "fNFailedChargeRatioCut  = " << fNFailedChargeRatioCut << endl;
      MSG("Align", Msg::kInfo) << "fNFailedStripChargeCut  = " << fNFailedStripChargeCut << endl;
   }

   MSG("Align", Msg::kInfo) << endl;
   
   if(fNMissedVPlane > 1){
      Double_t size = fNMissedVPlane;
      fMissedVPlaneMeanCharge = fMissedVPlaneMeanCharge/size;
      fMissedVPlaneChargeSigma = pow(fMissedVPlaneChargeSigma/size -
                                     fMissedVPlaneMeanCharge*fMissedVPlaneMeanCharge, 0.5);
      MSG("Align", Msg::kInfo) << "fMissedVPlaneMeanCharge = " << fMissedVPlaneMeanCharge
                               << " +/- " << fMissedVPlaneChargeSigma <<endl;
   }
   
   if(fNMissedUPlane > 1){
      Double_t size = fNMissedUPlane;
      fMissedUPlaneMeanCharge = fMissedUPlaneMeanCharge/size;
      fMissedUPlaneChargeSigma = pow(fMissedUPlaneChargeSigma/size -
                                     fMissedUPlaneMeanCharge*fMissedUPlaneMeanCharge, 0.5);
      MSG("Align", Msg::kInfo) << "fMissedUPlaneMeanCharge = " << fMissedUPlaneMeanCharge
                               << " +/- " << fMissedUPlaneChargeSigma <<endl;
   }     

    MSG("Align", Msg::kInfo) << endl;
}


//__________________________________________________________________________________________________
void  MakeAlignmentModule::PrintAlignmentStrip(const AlignmentStrip& astrip)
{   
   //Format msg service so that time and charge can be properly deisplayed
   MsgFormat ffmt("%10.9f");
   MsgFormat ffrt("%5.4f");
   MsgFormat fimt("%5d");
   
   MSG("Align", Msg::kVerbose) << PlexStripEndId(astrip.plexseid)
                               << " charge = " << fimt(astrip.charge)
                               << " begtime = " << ffmt(astrip.begtime)
                               << " ndigit = " << astrip.ndigit
                               << " residual = " << ffrt(astrip.residual) << endl;
}

---