AlignChamZNoBF.cxx

Go to the documentation of this file.

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

#include "Alignment/AlignChamZNoBF.h"

#include "Config/CfgConfig.h"
#include "Config/CfgParam.h"
#include "EventDataModel/EDMEventHandle.h"
#include "TrkRBase/TrackSeg.h"
#include "TrkRBase/TrkRUtil.h"

#include <TMinuit.h>
#include <TTree.h>

#include <iostream>
#include <cmath>

MODULE_DECL(AlignChamZNoBF);
using namespace std;

static AlignChamZNoBF* gModule = 0;

//......................................................................

AlignChamZNoBF::AlignChamZNoBF(const char* version) :
  JobCModule("AlignChamZNoBF")
{
  // Register module to be updated with a given config version
  SetCfgVersion(version);
  // Make sure that configuration has been loaded
  CheckInit();
}

//......................................................................

AlignChamZNoBF::~AlignChamZNoBF()
{
}

//......................................................................

JobCResult AlignChamZNoBF::Ana(const EDMEventHandle& evt)
{
  evt.Raw();
  evt.Cal();
  std::vector<const TrackSeg*> seg(0);
  try {
    evt.Reco().Get(fSegFolder.c_str(), seg); 
    if (seg.empty()) return kFailed;
    if (fDebug) cout << "Got " << seg.size() << " tracks" << endl;
  }
  catch (...) { return kFailed; }

  TrkInfo_t ti;
  double resid[6][4];

  for (unsigned int i = 0; i < seg.size(); ++i) {
    const TrackSeg* trk = seg[i];

    if (trk->NClust() < fMinNClust) continue;

    int nGoodClu = 0;

    memset(&ti, 0, sizeof(ti));

    for (int c = kDC1; c < kNCham; ++c) {
      for (int p = 0; p < kNPlane; ++p) {
    const WireClust* wc = trk->Clust(c, p);
    if (!wc || wc->NWire() > 1) continue;
    nGoodClu++;
    ti.fWire[c - kDC1][p] = wc->AvgWire();
    ti.fWeight[c - kDC1][p] = 1.0 / (wc->USig() * wc->USig());
      }
    }

    if (nGoodClu < fMinNClust) continue;

    // Make sure that track residuals are reasonable
    Fit(ti, resid);
    for (int c = 0; c < 6 && nGoodClu > 0; ++c) {
      for (int p = 0; p < kNPlane; ++p) {
    if (ti.fWeight[c][p] > 1 && 
        fabs(resid[c][p]) > 1.5 * GChamGeo::Pitch(c + kDC1)) {
      nGoodClu = 0;
      break;
    }
      }
    }

    if (nGoodClu < fMinNClust) continue;

    fTrackList.push_back(ti);
  }

  return kPassed;
}

//......................................................................
  
void AlignChamZNoBF::NewRun(int run, int /*subrun*/)
{
  GMIPPGeo::Instance(run);
}

//......................................................................

double AlignChamZNoBF::ComputeF()
{
  if (fUseCorr) return ComputeCorrF();
  else return ComputeChi2F();
}

//......................................................................

double AlignChamZNoBF::ComputeChi2F()
{
  double f = 0;

  double u[24];
  double w[24];
  double z[24];
  double cosA[24];
  double sinA[24];
  double z1 = GMIPPGeo::Instance().Cham(kDC1).ZGas() - 50;
  double z2 = GMIPPGeo::Instance().Cham(kPWC6).ZGas() + 50;
  double par[4];
  double parErr[4];

  for (unsigned int i = 0; i < fTrackList.size(); ++i) {
    TrkInfo_t& ti = fTrackList[i];

    int nClu = 0;
    
    for (int c = 0; c < 6; ++c) {
      GChamGeo& geo = GMIPPGeo::Instance().Cham(c + kDC1);
      for (int p = 0; p < kNPlane; ++p) {
    if (ti.fWeight[c][p] < 1e-3) continue;
    
    geo.WireToU(p + 1, ti.fWire[c][p], u[nClu]);
    w[nClu] = ti.fWeight[c][p];
    z[nClu] = geo.Z(p + 1);
    cosA[nClu] = geo.GAngleCos(p + 1);
    sinA[nClu] = geo.GAngleSin(p + 1);
    
    nClu++;
      }
    }
    
    f += TrkRUtil::FitLineFromWires(nClu, u, w, z, cosA, sinA, z1, z2, 
                    par, parErr);
  }
  return f;
}

//......................................................................

double AlignChamZNoBF::ComputeCorrF() {
  double resid[6][kNPlane];
  int n[6][kNPlane];
  double x[6][kNPlane];
  double y[6][kNPlane];
  double xx[6][kNPlane];
  double yy[6][kNPlane];
  double xy[6][kNPlane];
  memset(x, 0, sizeof(x));
  memset(n, 0, sizeof(n));
  memset(y, 0, sizeof(y));
  memset(xx, 0, sizeof(xx));
  memset(yy, 0, sizeof(yy));
  memset(xy, 0, sizeof(xy));

  for (unsigned int i = 0; i < fTrackList.size(); ++i) {
    TrkInfo_t& ti = fTrackList[i];
    double dxdz = Fit(ti, resid);

    for (int c = 0; c < 6; ++c) {
      for (int p = 0; p < kNPlane; ++p) {
    if (ti.fWeight[c][p] < 1) continue;
        n[c][p]++;
        x[c][p] += dxdz;
        xx[c][p] += dxdz * dxdz;
        y[c][p] += resid[c][p];
        yy[c][p] += resid[c][p] * resid[c][p];
        xy[c][p] += dxdz * resid[c][p];
      }
    }
  }
  double f = 0;
  for (int c = 0; c < 6; ++c) {
    for (int p = 0; p < kNPlane; ++p) {
      if (n[c][p] < 10) continue;

      // Add square of correlation coefficient for this plane
      double corr = (n[c][p] * xy[c][p] - x[c][p] * y[c][p]);
      corr *= corr;
      corr /= (n[c][p] * xx[c][p] - x[c][p] * x[c][p]) *
        (n[c][p] * yy[c][p] - y[c][p] * y[c][p]);
      assert(fabs(corr) < 1);

      f += corr;
    }
  }
  return f;
}

//......................................................................

double AlignChamZNoBF::Fit(TrkInfo_t& ti, double resid[6][4], float* dydz,
             float xc[6][4], float yc[6][4])
{
  double u[24];
  double w[24];
  double z[24];
  double cosA[24];
  double sinA[24];
  double z1 = GMIPPGeo::Instance().Cham(kDC1).ZGas() - 50;
  double z2 = GMIPPGeo::Instance().Cham(kPWC6).ZGas() + 50;
  double par[4];
  double parErr[4];
  int plane[24];

  int nClu = 0;

  for (int c = 0; c < 6; ++c) {
    GChamGeo& geo = GMIPPGeo::Instance().Cham(c + kDC1);
    for (int p = 0; p < kNPlane; ++p) {
      if (ti.fWeight[c][p] < 1e-3) {
    resid[c][p] = 1e6;
    continue;
      }

      geo.WireToU(p + 1, ti.fWire[c][p], u[nClu]);
      w[nClu] = ti.fWeight[c][p];
      z[nClu] = geo.Z(p + 1);
      cosA[nClu] = geo.GAngleCos(p + 1);
      sinA[nClu] = geo.GAngleSin(p + 1);
      plane[nClu] = c * 4 + p;

      nClu++;
    }
  }

  TrkRUtil::FitLineFromWires(nClu, u, w, z, cosA, sinA, z1, z2, 
                 par, parErr);

  double x1[3] = {par[0], par[1], z1};
  double x2[3] = {par[2], par[3], z2};
  double x[3];

  double dxdz = (par[2] - par[0]) / (z2 - z1);
  if (dydz) *dydz = (par[3] - par[1]) / (z2 - z1);
  if (xc && yc) {
    for (int c = 0; c < 6; ++c) {
      GChamGeo& geo = GMIPPGeo::Instance().Cham(c + kDC1);
      for (int p = 0; p < 4; ++p) {
    x[2] = geo.Z(p + 1);
    TrkRUtil::PredictXY(x1, 0, x2, 0, x, 0);
    xc[c][p] = x[0];
    yc[c][p] = x[1];
      }
    }
  }

  for (int i = 0; i < nClu; ++i) {
    double saveW = w[i];
    w[i] = 0;
    TrkRUtil::FitLineFromWires(nClu, u, w, z, cosA, sinA, z1, z2, 
                   par, parErr);
    w[i] = saveW;

    x1[0] = par[0];
    x1[1] = par[1];
    x2[0] = par[2];
    x2[1] = par[3];
    x[2] = z[i];
    TrkRUtil::PredictXY(x1, 0, x2, 0, x, 0);
    
    int c = plane[i] / 4;
    int p = plane[i] % 4;
    GMIPPGeo::Instance().Cham(c + kDC1).XYToU(p + 1, x[0], x[1], 
                          resid[c][p]);
    resid[c][p] -= u[i];
  }

  return dxdz;
}

//......................................................................

static void gsFCN(int& /*npar*/, double* /*gin*/, double& f,
                  double* par, int /*iflag*/)
{
  // Setup geometry based on the new constants
  for (int i = 0; i < 6; ++i) {
    GMIPPGeo::Instance().Cham(i + kDC1).SetZGas(par[i]);
    GMIPPGeo::Instance().Cham(i + kDC1).SetGAngle(par[i + 6] * 1e-3);
  }

  f = gModule->ComputeF();
}

//......................................................................
  
void AlignChamZNoBF::EndRun(int /*run*/, int /*subrun*/)
{
  if (fTrackList.size() < 100) {
    cout << "Can't do alignment with " << fTrackList.size() << " tracks"
     << endl;
    return;
  }

  cout << "Doing no B-field alignment with " << fTrackList.size()
       << " tracks" << endl;

  MakeTree("pre");

  gModule = this;
  TMinuit minu(12);
  minu.SetFCN(gsFCN);
  
  char name[256];
  int ierr;
  double savePos[6], saveAng[6];
  for (int i = 0; i < 6; ++i) {
    GChamGeo& geo = GMIPPGeo::Instance().Cham(i + kDC1);
    sprintf(name, "%s_z", geo.Name());
    minu.mnparm(i, name, geo.ZGas(), 0.1, 
        geo.ZGas() - 10, geo.ZGas() + 10, ierr);
    savePos[i] = geo.ZGas();

    sprintf(name, "%s_angle", geo.Name());
    saveAng[i] = geo.GAngle(1) - geo.LAngle(1);
    minu.mnparm(i + 6, name, saveAng[i] * 1e3, 1, -5e2, 5e2, ierr);
  }

  for (unsigned int i = 0; i < fFixCham.size(); ++i) 
    minu.FixParameter(fFixCham[i]);

  minu.Migrad();

  MakeTree("post");

  double par, pe;
  for (int i = 0; i < 6; ++i) {
    minu.GetParameter(i, par, pe);
    cout << GChamGeo::Name(i + kDC1) << " z shift: " << par - savePos[i]
     << " cm, err: " << pe << endl;

    minu.GetParameter(i + 6, par, pe);
    cout << GChamGeo::Name(i + kDC1) << " angle shift: " 
     << par * 1e-3 - saveAng[i] << " deg, err: " << pe << endl;
  }
}

//......................................................................

void AlignChamZNoBF::MakeTree(const char* name)
{
  char str[256];
  sprintf(str, "%s-alignment Tree", name);
  TTree* tree = new TTree(name, str);

  double u[6][4];
  float x[6][4];
  float y[6][4];
  double uperp[6][4];
  double resid[6][4];
  float dxdz;
  float dydz;
  tree->Branch("u", &u, "u[6][4]/D");
  tree->Branch("x", &x, "x[6][4]/D");
  tree->Branch("y", &y, "y[6][4]/D");
  tree->Branch("uperp", &uperp, "uperp[6][4]/D");
  tree->Branch("resid", &resid, "resid[6][4]/D");
  tree->Branch("dxdz", &dxdz, "dxdz/F");
  tree->Branch("dydz", &dydz, "dydz/F");
  
  for (unsigned int i = 0; i < fTrackList.size(); ++i) {
    TrkInfo_t& ti = fTrackList[i];

    dxdz = Fit(ti, resid, &dydz, x, y);

    for (int c = 0; c < 6; ++c) {
      GChamGeo& geo = GMIPPGeo::Instance().Cham(c + kDC1);
      for (int p = 0; p < 4; ++p) {
    if (ti.fWeight[c][p] < 1e-3) u[c][p] = 1e6;
    else geo.WireToU(p + 1, ti.fWire[c][p], u[c][p]);
    geo.XYToUperp(p + 1, x[c][p], y[c][p], uperp[c][p]);
      }
    }
    tree->Fill();
  }

  tree->Write();
  delete tree;
}

//......................................................................

void AlignChamZNoBF::Update(const CfgConfig& c)
{
  fFixCham.clear();
  c("Debug").Get(fDebug);
  c("SegFolder").Get(fSegFolder);
  c("MinNClust").Get(fMinNClust);
  try { c("FixCham").Get(fFixCham); }
  catch (...) {
    int i;
    c("FixCham").Get(i);
    fFixCham.push_back(i);
  }
  c("UseCorr").Get(fUseCorr);

  // Change the flag to signify that our configuration has been read
  fIsInit = true;
}

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

---