AlignBC.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// $Id: AlignBC.cxx,v 1.3 2007/07/05 21:55:40 jpaley Exp $
// 
// lebedev@physics.harvard.edu
////////////////////////////////////////////////////////////////////////

#include "Alignment/AlignBC.h"

#include "ChamCalib/ChamDrift.h"
#include "Config/CfgConfig.h"
#include "Config/CfgParam.h"
#include "EventDataModel/EDMEventHandle.h"
#include "TrkRBase/TrackSeg.h"
#include "TrkRBase/WireCross.h"
#include "TrkRBase/WireClust.h"
#include "TrkRBase/TrkRUtil.h"
#include "TrkRBase/TrkChamGeo.h"

#include <TDirectory.h>
#include <TTree.h>
#include <TH1F.h>
#include <TF1.h>

#include <iostream>
#include <cassert>
#include <cmath>

MODULE_DECL(AlignBC);
ClassImp(BCAlignNt);

using namespace std;
using namespace TrkRUtil;

AlignBC* AlignBC::fsModule = 0;
double AlignBC::fZ1 = 0;
double AlignBC::fZ2 = 0;

//----------------------------------------------------------------------

AlignBC::AlignBC(const char* version) :
  JobCModule("AlignBC")
{
  SetCfgVersion(version);
  CheckInit();
}

//----------------------------------------------------------------------

JobCResult AlignBC::Ana(const EDMEventHandle& evt)
{
  vector<const TrackSeg*> vec(0);
  // Make sure raw and calibrated data buckets are loaded
  evt.Raw(); 
  evt.Cal();
  try {evt.Reco().Get(fTrkFolder.c_str(), vec); }
  catch (...) {return kFailed;}

  int nTrk = vec.size();
  if (nTrk == 0 || nTrk > fMaxNTrackEvt) return kFailed;

  TrkInfo_t *trkInfo = 0;
  
  // Copy vector info into our private struct
  for (int i = 0; i < nTrk; ++i) {
    const TrackSeg* seg = vec[i];

    if (!trkInfo) trkInfo = new TrkInfo_t;

    int nGoodClu = seg->NClust();
    // Skip tracks with too few clusters
    if (nGoodClu < fMinNClust) continue;

    for (int cham = 0; cham < kNBC; ++cham) {
      for (int p = 0; p < kNPlane; ++p) {
    const WireClust* wc = seg->Clust(cham, p);
    if (!wc || wc->NWire() > 1) {
      trkInfo->fWeight[cham][p] = -1;
      nGoodClu--;
      if (nGoodClu < fMinNClust) {
        cham = kNBC;
        break;
      }
    }
    else {
      assert(wc->Plane() == p + 1);
      trkInfo->fWire[cham][p] = wc->AvgWire();
      trkInfo->fWeight[cham][p] = 1 / (wc->USig() * wc->USig());
      trkInfo->fT[cham][p] = wc->T();
    }
      }
    }
    if (nGoodClu < fMinNClust) continue;
    fTrkInfo.push_back(trkInfo);
    trkInfo = 0;
  }

  return kPassed;
}

//----------------------------------------------------------------------

void AlignBC::NewRun(int /*run*/, int /*subrun*/)
{
}

//----------------------------------------------------------------------

void AlignBC::EndRun(int run, int /*subrun*/)
{
  if ((int) fTrkInfo.size() < fMinNTrack) {
    cout << "AlignBC: " << fTrkInfo.size() << " is not enough tracks"
     << endl;
    return;
  }

  fZ1 = GMIPPGeo::Instance().Cham(kBC1).ZGas();
  fZ2 = GMIPPGeo::Instance().Cham(kBC3).ZGas();
  
  double startPar[kNBC][kNPlane];
  for (int c = 0; c < kNBC; ++c) {
    GChamGeo& geo = GMIPPGeo::Instance().Cham(c);
    for (int p = 0; p < kNPlane; ++p) startPar[c][p] = geo.WireZero(p + 1);
  }

  if (!RoughAlign()) return;

  for (int c = 0; c < kNBC; ++c) {
    GChamGeo& geo = GMIPPGeo::Instance().Cham(c);
    cout << geo.Name() << " wire 0 changes: ";
    for (int p = 0; p < kNPlane; ++p) {
      cout << geo.WireZero(p + 1) - startPar[c][p] << " ";
    }
    cout << endl;
  }

  if (fFillNt) FillNt("nt1", false);

  if (fDoFineAlign) {
    FineAlign();

    for (int c = 0; c < kNBC; ++c) {
      GChamGeo& geo = GMIPPGeo::Instance().Cham(c);
      cout << geo.Name() << " wire 0 changes: ";
      for (int p = 0; p < kNPlane; ++p) {
    cout << geo.WireZero(p + 1) - startPar[c][p] << " ";
      }
      cout << endl;
    }

    if (fFillNt) FillNt("nt2", true);
  }

  WriteToDB(run);
}

//----------------------------------------------------------------------

void AlignBC::FillNt(const char* name, bool resolve)
{
  BCAlignNt *nt = new BCAlignNt;
  TTree* tree = new TTree(name, "BC Alignment Ntuple");
  tree->Branch("nt", "BCAlignNt", &nt);

  double x1[3], x2[3], x3[3];
  x1[2] = fZ1;
  x2[2] = fZ2;

  for (unsigned int i = 0; i < fTrkInfo.size(); ++i) {
    if (!Fit(fTrkInfo[i], nt, resolve)) continue;

    x1[0] = nt->par[0];
    x1[1] = nt->par[1];
    x2[0] = nt->par[2];
    x2[1] = nt->par[3];

    for (int c = 0; c < kNBC; ++c) {
      GChamGeo& geo = GMIPPGeo::Instance().Cham(c);
      for (int p = 0; p < kNPlane; ++p) {
    x3[2] = geo.Z(p + 1);
    PredictXY(x1, 0, x2, 0, x3, 0);
    nt->uperp[c][p] = x3[0] * geo.GAngleSin(p + 1) + 
      x3[1] * geo.GAngleCos(p + 1);
    nt->ufit[c][p] = x3[0] * geo.GAngleCos(p + 1) -
      x3[1] * geo.GAngleSin(p + 1);
      }
    }
    tree->Fill();
  }
  tree->Write();
  delete tree;
  delete nt;
}

//----------------------------------------------------------------------

void AlignBC::Update(const CfgConfig& c)
{
  c("DoFineAlign").Get(fDoFineAlign); 
  c("MaxAlignAttempt").Get(fMaxAlignAttempt);
  c("MinNClust").Get(fMinNClust);
  c("MinNTrack").Get(fMinNTrack);
  c("MaxNTrackEvt").Get(fMaxNTrackEvt);
  c("TrkFolder").Get(fTrkFolder);
  c("OkRoughShift").Get(fOkRoughShift);
  c("OkFineShift").Get(fOkFineShift);
  c("MinTime").Get(fMinTime);
  c("MaxTime").Get(fMaxTime);
  c("FillNt").Get(fFillNt);

  fIsInit = true;
}

//----------------------------------------------------------------------

bool AlignBC::Fit(const TrkInfo_t* ti, BCAlignNt* nt, 
             bool resolve /* = false */)
{
  double w[12];
  double u[12];
  double du[12]; // Shift proportional to drift time
  double z[12];
  double cosA[12];
  double sinA[12];
  short  plane[12];

  ChamDrift& drift = ChamDrift::Instance();
  
  // Copy information from TrkInfo_t into arrays
  nt->nclust = 0;
  for (int cham = 0; cham < kNBC; ++cham) {
    GChamGeo& geo = GMIPPGeo::Instance().Cham(cham);
    for (int p = 0; p < kNPlane; ++p) {
      if (ti->fWeight[cham][p] <= 0) {
    nt->u[cham][p] = 1e6;
    continue;
      }

      if (resolve) {
    if(ti->fT[cham][p] < fMinTime || 
       ti->fT[cham][p] > fMaxTime) continue;
    du[nt->nclust] = 
      drift.TimeToDist(cham, p, ti->fT[cham][p], &w[nt->nclust]);
    w[nt->nclust] = 1 / (w[nt->nclust] * w[nt->nclust]);
      }
      else w[nt->nclust] = ti->fWeight[cham][p];
      
      geo.WireToU(p + 1, ti->fWire[cham][p], u[nt->nclust]);
      nt->u[cham][p] = u[nt->nclust];
      z[nt->nclust] = geo.Z(p + 1);
      cosA[nt->nclust] = geo.GAngleCos(p + 1);
      sinA[nt->nclust] = geo.GAngleSin(p + 1);
      plane[nt->nclust] = cham * kNPlane + p;
      nt->nclust++;
    }
  }

  if (nt->nclust < 5) return false;

  // Resolve ambiguity if we need to
  double par[4];
  if (resolve) {
    int bestMask = 0;
    nt->chi2 = 1e6;
    int nCombo = (1 << nt->nclust);

    double uCopy[12];
    for (int m = 0; m < nCombo; ++m) {
      bool okCombo = true; 

      for (int i = 0; i < nt->nclust; ++i) {
    if (m & (1 << i)) {
      // We don't have to double count combinations with du=0
      if (du[i] < 1e-6) {
        okCombo = false;
        break;
      }
      uCopy[i] = u[i] - du[i];
    }
    else uCopy[i] = u[i] + du[i];
      }
      if (!okCombo) continue;

      double c2 = FitLineFromWires(nt->nclust, uCopy, w, z, cosA, sinA, 
                   fZ1, fZ2, par, 0);
      if (c2 < nt->chi2) {
    nt->chi2 = c2;
    bestMask = m;
      }
    }

//     cout << "Resolved with best mask " << (void*) bestMask
//   << ", chi2 = " << nt->chi2 << "/" << nt->nclust - 4 << endl;
    for (int i = 0; i < nt->nclust; ++i) {
      if (bestMask & (1 << i)) u[i] -= du[i];
      else u[i] += du[i];
      nt->u[plane[i] / 4][plane[i] % 4] = u[i];
    }
  }

  // Save chi squared for the full fit
  nt->chi2 = FitLineFromWires(nt->nclust, u, w, z, cosA, sinA, 
                  fZ1, fZ2, par, 0);

  for (int i = 0; i < 4; ++i) nt->par[i] = par[i];
  nt->z1 = fZ1;
  nt->z2 = fZ2;

  double wCopy[12];

  double x1[3], x2[3], x3[3];
  x1[2] = fZ1;
  x2[2] = fZ2;

  // Now exclude one plane at a time, and save residuals
  for (int c = 0; c < kNBC; ++c) {
    GChamGeo& geo = GMIPPGeo::Instance().Cham(c);
    for (int p = 0; p < kNPlane; ++p) {
      if (ti->fWeight[c][p] <= 0) {
    nt->chi2rm[c][p] = 1e6;
    nt->resid[c][p] = 1e6;
    continue;
      }
      
      memcpy(wCopy, w, 12 * sizeof(double));
      
      int i = 0;
      for ( ; i < nt->nclust; ++i) {
    if (plane[i] == c * kNPlane + p) {
      wCopy[i] = 0;
      break;
    }
      }
      if (i >= nt->nclust) {
    nt->chi2rm[c][p] = 1e6;
    nt->resid[c][p] = 1e6;
    continue;
      }
      
      nt->chi2rm[c][p] = 
    FitLineFromWires(nt->nclust, u, wCopy, z, cosA, sinA, 
               fZ1, fZ2, par, 0);
      x1[0] = par[0];
      x1[1] = par[1];
      x2[0] = par[2];
      x2[1] = par[3];
      x3[2] = geo.Z(p + 1);

      PredictXY(x1, 0, x2, 0, x3, 0);
      nt->resid[c][p] = x3[0] * geo.GAngleCos(p + 1) - 
    x3[1] * geo.GAngleSin(p + 1) - u[i];
    }
  }

  return true;
}

//----------------------------------------------------------------------

bool AlignBC::RoughAlign()
{
  TH1F* resH[kNBC][kNPlane];
  
  char str[256];
  for (int c = 0; c < kNBC; ++c) {
    for (int p = 0; p < kNPlane; ++p) {
      sprintf(str, "resid%d%d", c, p);
      resH[c][p] = new TH1F(str, str, 400, -2, 2); // 1% wire spacing bins
      resH[c][p]->SetXTitle("Residual (wire spacing)");
    }
  }

  BCAlignNt nt;
  TDirectory* dir = gDirectory;

  for (int at = 0; ; ++at) {
    if (at >= fMaxAlignAttempt) return false;
    cout << "BC align attempt " << at << " with " << fTrkInfo.size()
     << " tracks" << endl;
    sprintf(str, "attempt%d", at);
    gDirectory->mkdir(str);
    gDirectory->cd(str);

    for (unsigned int i = 0; i < fTrkInfo.size(); ++i) {
      if (!Fit(fTrkInfo[i], &nt)) continue;
      for (int c = 0; c < kNBC; ++c) {
    for (int p = 0; p < kNPlane; ++p) { 
      if (nt.resid[c][p] > 100) continue;
      resH[c][p]->Fill(nt.resid[c][p] / GChamGeo::Pitch(c));
    }
      }
    }

    double res[3][4];
    double maxRes = 0;
    int maxResP = -1;
    for (int c = 0; c < kNBC; ++c) {
      GChamGeo& geo = GMIPPGeo::Instance().Cham(c);
      for (int p = 0; p < kNPlane; ++p) {
    if (resH[c][p]->GetEntries() > 50) {
      res[c][p] = resH[c][p]->GetMean();
      if (fabs(maxRes) < fabs(res[c][p])) {
        maxRes = res[c][p];
        maxResP = c * kNPlane + p;
      }
    }
    resH[c][p]->Write();
    resH[c][p]->Reset();

    // Shift wire plane by a fraction of residual
    geo.SetWireZero(p + 1, geo.WireZero(p + 1) -
            geo.WireDir(p + 1) * res[c][p] * 0.3);

      }
    }

    dir->cd();

    if (maxResP < 0) return false;

    cout << "Max residual in " << GChamGeo::Name(maxResP / 4) 
     << " P" << maxResP % 4 + 1 << ": " << maxRes << " w.s." << endl;
    if (fabs(maxRes) < 0.25 && PreenTrackList(res)) continue;
    if (fabs(maxRes) < fOkRoughShift) break;
  }

  for (int c = 0; c < kNBC; ++c) {
    for (int p = 0; p < kNPlane; ++p) delete resH[c][p];
  }
  return true;
}

//----------------------------------------------------------------------

bool AlignBC::PreenTrackList(double res[3][4])
{
  int removed = 0;
  BCAlignNt nt;
  int c, p;

  vector<TrkInfo_t*>::iterator itr(fTrkInfo.begin());
  int n = fTrkInfo.size();
  for ( ; itr != fTrkInfo.end(); ) {
    const TrkInfo_t* ti = *itr;
    n--;
    if (!Fit(ti, &nt)) {
      delete ti;
      itr = fTrkInfo.erase(itr);
      removed++;
      continue;
    }

    for (c = 0; c < 3; ++c) {
      for (p = 0; p < 4; ++p) {
    double r = fabs(nt.resid[c][p] / GChamGeo::Pitch(c) - res[c][p]) ;
    if (r < 1e3 && r > 1.5 + fabs(res[c][p]))
      break;
      }
      if (p != 4) break;
    }

    if (c != 3) {
      delete ti;
      itr = fTrkInfo.erase(itr);
      removed++;
    }
    else ++itr;
  }

  return (removed != 0);
}

//----------------------------------------------------------------------

bool AlignBC::FineAlign()
{

  TH1F* resH[kNBC][kNPlane];
  
  char str[256];
  for (int c = 0; c < kNBC; ++c) {
    for (int p = 0; p < kNPlane; ++p) {
      sprintf(str, "resid%d%d", c, p);
      resH[c][p] = new TH1F(str, str, 240, -0.3, 0.3); // 25 micron bins
      resH[c][p]->SetXTitle("Residual (cm)");
    }
  }

  TF1* fu = new TF1("fu", "gaus", -1, 1);

  BCAlignNt nt;
  TDirectory* dir = gDirectory;

  for (int at = 0; ; ++at) {
    cout << "BC Fine align attempt " << at << endl;
    sprintf(str, "fine%d", at);
    gDirectory->mkdir(str);
    gDirectory->cd(str);

    double maxShift = 0;

    int nTrk = 0;
    for (unsigned int i = 0; i < fTrkInfo.size(); ++i) {
      if (!Fit(fTrkInfo[i], &nt, true)) continue;
      nTrk++;
      for (int c = 0; c < kNBC; ++c) {
    for (int p = 0; p < kNPlane; ++p) { 
      if (fabs(nt.resid[c][p]) > 0.3) continue;
      resH[c][p]->Fill(nt.resid[c][p]);
    }
      }
    }

    if (nTrk < 100) return false;

    for (int c = 0; c < kNBC; ++c) {
      GChamGeo& geo = GMIPPGeo::Instance().Cham(c);
      for (int p = 0; p < kNPlane; ++p) {
    if (resH[c][p]->GetEntries() > 50) {
      fu->SetParameter(0, resH[c][p]->GetMaximum());
      fu->SetParLimits(0, resH[c][p]->GetMaximum() * 0.7,
               resH[c][p]->GetMaximum() * 1.5);
      fu->SetParameter(1, resH[c][p]->GetMean());
      fu->SetParLimits(1, resH[c][p]->GetMean() - 0.2, 
               resH[c][p]->GetMean() + 0.2);
      fu->SetParameter(2, resH[c][p]->GetRMS());
      resH[c][p]->Fit(fu, "Q0+");

      double res = fu->GetParameter(1);
      res *= 0.3 / geo.Pitch();
      if (maxShift < fabs(res)) maxShift = fabs(res);
      
      geo.SetWireZero(p + 1, geo.WireZero(p + 1) + res);
    }
    resH[c][p]->Write();
    resH[c][p]->Reset();
      }
    }

    dir->cd();
    cout << "Max shift: " << maxShift << endl;
    if (maxShift < fOkFineShift) break;
  }

  for (int c = 0; c < kNBC; ++c) {
    for (int p = 0; p < kNPlane; ++p) delete resH[c][p];
  }
  delete fu;

  return true;
}

//----------------------------------------------------------------------

void AlignBC::WriteToDB(int run)
{
  static int minRun = 10000000;
  static int maxRun = 0;

  if (strcmp(getenv("PGUSER"), "mippdbwrite") != 0) {
    cout << "User " << getenv("PGUSER")
     << " is not allowed to write to DB" << endl;
    return;
  }

  TrkChamGeo& db = TrkChamGeo::Instance();
  for (int c = 0; c < kNBC; ++c) db.SetCham(c);
  if (minRun > run) minRun = run;
  if (maxRun < run) maxRun = run;
  db.SetRunRange(minRun, maxRun);
  db.WriteToDB(TrkChamGeo::kAlignBC, 7); // Write BC alignment only

}

////////////////////////////////////////////////////////////////////////

---