StraightTrackAlignment.cxx

Go to the documentation of this file.

#include "Alignment/StraightTrackAlignment.h"
#include "Alignment/ScintModule.h"
#include "Alignment/AlignHists.h"

#include "RecoBase/CandStripHandle.h"
#include "RecoBase/CandTrackHandle.h"
#include "UgliGeometry/UgliScintPlnHandle.h"
#include "Navigation/NavSet.h"

#include <fstream>
#include <vector>
#include <cassert>
using std::vector;
using std::list;

using std::cerr;
using std::cout;
using std::endl;

StraightTrackAlignment::StraightTrackAlignment()
    : fHister(0)
      ,fMinPlanes(14)            // hard code for now.
      ,fHistName(0)
{
}

StraightTrackAlignment::~StraightTrackAlignment()
{
    if (fHister) delete fHister;
}

// Create the histogrammer
void StraightTrackAlignment::Init(UgliGeomHandle& ugh, const char* histname)
{
    fHistName = histname;
    vector<UgliScintPlnHandle> planevec = ugh.GetScintPlnHandleVector();
    fNplanes = planevec.size();
    fNmdlperplane = planevec[0].GetScintMdlHandleVector().size();
}

// return the ScintModule containing the given USH.
ScintModule* StraightTrackAlignment::LookupScintModule(UgliStripHandle ush)
{
    UgliScintMdlHandle usmh = ush.GetScintMdlHandle();
    PlexScintMdlId smid = usmh.GetPlexScintMdlId();

    ScintModuleLookup::iterator it = fSMLookup.find(smid);
    if (it == fSMLookup.end()) {  // D.N.E.
        ScintModule* smod = new ScintModule(usmh);
        assert(smod);
        fSMLookup[smid] = smod;
        return smod;
    }

    assert(it->second);
    return it->second;
}

// Save any tracks
void StraightTrackAlignment::AddTrack (const CandTrackHandle& track_handle)
{
    int nhits = track_handle.GetNDaughters();
    if (nhits < fMinPlanes) return;

//    cerr << nhits << " hits\n";

    UgliGeomHandle ugh(*(track_handle.GetVldContext()));
    CandStripHandleItr sitr(track_handle.GetDaughterIterator());
    
    StraightTrack st[2];
    CandStripHandle* csh = 0;
    while ((csh = sitr())) {
        cerr << ".";
        PlexStripEndId seid = csh->GetStripEndId();
        int view = 0;
        if (seid.GetPlaneView() != PlaneView::kV) view = 1;
        st[view].push_back(ugh.GetStripHandle(seid));
    }
    cerr << endl;

    if (st[0].size()) {
        fTracks[0].push_back(st[0]);
    }
    if (st[1].size()) {
        fTracks[1].push_back(st[1]);
    }
}


// Fits a linear track to the StraighTrack's hits between "lit" and
// "end" while skipping the hit "skip".  The predicted transverse
// position of "skip" is fed to skip's ScintModule.
double
StraightTrackAlignment::FitTrackLessOne(StraightTrack::iterator lit,
                                        StraightTrack::iterator end,
                                        StraightTrack::iterator skip)
{
    // Linear regression parameters.  See for eg. Num Rec Ch 15.
    // x is Z, y is TPos (U or V).
    double S, Sx, Sy, Sxx, Sxy;
    S=Sx=Sy=Sxx=Sxy=0.0;

    int count = 0;

    for (;lit != end; ++lit) {
        if (lit == skip) continue;

        UgliStripHandle ush = (*lit);
        ScintModule* smod = this->LookupScintModule(ush);

        float x = ush.GlobalPos(0)(2);
        float y = ush.GetTPos() + smod->GetOffset();
        
        Sx += x;
        Sy += y;
        Sxx += x*x;
        Sxy += x*y;        

        ++count;
    }
    S = count;

    // Note, since I assume a constant uncertainty in the transverse
    // position measurement, it drops out of hte final result.

    // The linear fit: y = a + bx
    double delta = S*Sxx - Sx*Sx;
    double a = (Sxx*Sy - Sx*Sxy)/delta;
    double b = (S*Sxy - Sx*Sy)/delta;

    double predict = a + b*(*skip).GlobalPos(0)(2);
    predict -= (*skip).GetTPos();
    ScintModule* smod = this->LookupScintModule(*skip);
    if (smod) smod->AccumResid(predict);
    else {
        cerr << "Failed to lookup scint module!\n";
    }
    return predict;
}

// Fits track once for each digit.
double StraightTrackAlignment::ApplyTrack(StraightTrack st, 
                                          PlaneView::PlaneView_t view)
{
    double tot = 0;
    StraightTrack::iterator lit, first=st.begin(), end = st.end();

    vector<int> planev, mdlv, stripv;
    vector<double> residv;

    int count = 0;
    for (lit = first; lit != end; ++lit) {
//        cerr << ".";
        double resid = this->FitTrackLessOne(first,end,lit);
        if (fHister) {
            UgliScintMdlHandle usmh = (*lit).GetScintMdlHandle();
            PlexPlaneId ppid = usmh.GetPlexPlaneId();
            residv.push_back(resid);
            planev.push_back(ppid.GetPlane());
            stripv.push_back((*lit).GetSEId().GetStrip());
            mdlv.push_back(usmh.GetModuleNum());
        }
        tot += resid;
        ++count;
    }
    if (fHister) fHister->ApplyTrack(view,stripv,planev,mdlv,residv);
    
    return count ? tot/count : 0;
}


double StraightTrackAlignment::ApplyAllTracks(PlaneView::PlaneView_t the_view)
{
    int view = 0;
    if (the_view == PlaneView::kV) view = 1;
    list<StraightTrack>::iterator lit, end = fTracks[view].end();

    cerr << "Applying " << fTracks[view].size() << " tracks in view " 
         << view << endl;

    double tot = 0;
    int count = 0;
    for (lit=fTracks[view].begin(); lit != end; ++lit) {
        tot += this->ApplyTrack(*lit,the_view);
        ++ count;
    }
    cerr << "Done applying tracks\n";
    return count ? tot/count : 0;
}
double StraightTrackAlignment::ApplyAllTracks()
{
    static bool been_here = false;
    if (! been_here) {
        fHister = new AlignHists(fNplanes,fNmdlperplane,fHistName);
        been_here = true;
    }
                                 

    double ret = 0;
    ret += this->ApplyAllTracks(PlaneView::kU);
    ret += this->ApplyAllTracks(PlaneView::kV);
    return ret/2.0;
}

void StraightTrackAlignment::ApplyAllOffsets(void)
{
    ScintModuleLookup::iterator mit, end = fSMLookup.end();
    for (mit = fSMLookup.begin(); mit != end; ++mit) {
        ScintModule* sm = mit->second;
        if (!sm) {
            cerr << "Bad ScintModule in lookup\n";
            continue;
        }
        sm->UpdateOffsets();
        sm->ClearStatistics();
    }
    if (fHister) fHister->IncrementIteration();
}

---