AlignChamZ.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// $Id: AlignChamZ.cxx,v 1.4 2007/02/20 02:38:06 lebedev Exp $
//
// lebedev@physics.harvard.edu
////////////////////////////////////////////////////////////////////////

#include "Alignment/AlignChamZ.h"

#include "Bfield/bfield.h"
#include "Bfield/BFMagnet.h"
#include "Config/CfgConfig.h"
#include "Config/CfgParam.h"
#include "EventDataModel/EDMEventHandle.h"
#include "TrkRBase/TrkChamGeo.h"
#include "TrkRBase/TrkCand.h"
#include "TrkRBase/WireClust.h"
#include "TrkRBase/TrkRUtil.h"

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

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

MODULE_DECL(AlignChamZ);

using namespace std;
using namespace TrkRUtil;

static AlignChamZ* gModule;

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

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

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

AlignChamZ::~AlignChamZ()
{
}

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

JobCResult AlignChamZ::Ana(const EDMEventHandle& evt)
{
  vector<const TrkCand*> vec(0);

  try { evt.Reco().Get(fFolder.c_str(), vec); }
  catch (...) { return kFailed; }

  TrkInfo_t ti;

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

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

    int nClu = 0;
    for (int j = 0; j < trk->NClust(); ++j) {
      const WireClust* wc = trk->Clust(j);
      int cham = wc->Chamber();
      int p = wc->Plane() - 1;

      if (cham > kBC3) nClu++;

      ti.fWire[cham][p] = wc->AvgWire();
      ti.fWeight[cham][p] = 1.0 / (wc->USig() * wc->USig());
    }
    if (nClu < fMinNClust) continue;
    fTrackList.push_back(ti);
  }

  return kPassed;
}

//......................................................................
  
void AlignChamZ::NewRun(int run, int subrun)
{
  GMIPPGeo::Instance(run);
  config_bfield(run, subrun);
  if (fLoadAlignment) TrkChamGeo::Instance().Init(run);
  fZ0 = GMIPPGeo::Instance().Cham(kDC1).ZGas() - 20;
}

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

double AlignChamZ::ComputeF()
{
  CalcLambda(fZ0, fLambdaX, fLambdaY);
  
  if (fUseCorr) return ComputeCorrF();
  else return ComputeChi2F();
}

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

double AlignChamZ::ComputeChi2F()
{
  double u[kNCham * kNPlane];
  double w[kNCham * kNPlane];
  double z[kNCham * kNPlane];
  double cosA[kNCham * kNPlane];
  double sinA[kNCham * kNPlane];
  double lambdaX[kNCham * kNPlane];
  double lambdaY[kNCham * kNPlane];
  int chamber[kNCham * kNPlane];

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

    int nClust = 0;
    for (int cham = 0; cham < kNCham; ++cham) {
      GChamGeo& geo = GMIPPGeo::Instance().Cham(cham);
      for (int p = 0; p < kNPlane; ++p) {
    if (ti.fWeight[cham][p] < 1e-3) continue;
    geo.WireToU(p + 1, ti.fWire[cham][p], u[nClust]);
    w[nClust] = ti.fWeight[cham][p];
    z[nClust] = geo.Z(p + 1);
    cosA[nClust] = geo.GAngleCos(p + 1);
    sinA[nClust] = geo.GAngleSin(p + 1);
    lambdaX[nClust] = fLambdaX[cham][p];
    lambdaY[nClust] = fLambdaY[cham][p];
    chamber[nClust] = cham * kNPlane + p;
    nClust++;
      }
    }

    double par[5];
    FitCurvedTrack(nClust, u, w, z, cosA, sinA, fZ0, lambdaX, lambdaY,
           par, 0);
    
    // Reinitialize lambda's based on the first-try fit
    double lx[kNCham][kNPlane], ly[kNCham][kNPlane];
    CalcLambda(fZ0, lx, ly, par);
    for (int i = 0; i < nClust; ++i) {
      int cham =chamber[i] / kNPlane;
      int pl = chamber[i] % kNPlane;
      lambdaX[i] = lx[cham][pl];
      lambdaY[i] = ly[cham][pl];
    }
    double chi2 = FitCurvedTrack(nClust, u, w, z, cosA, sinA, fZ0, 
                 lambdaX, lambdaY, par, 0);
    if (chi2 < 1e5) f += chi2;
    else {
      cout << "Failed fit?: chi2=" << chi2 << " (" 
       <<  par[0] << ", " << par[1]
       << "; " << par[2] << ", " << par[3] << "; " << par[4]
       << endl;
    }
  }
  return f;
}

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

double AlignChamZ::ComputeCorrF()
{
  double resid[kNCham][kNPlane];
  int n[kNCham][kNPlane];
  double x[kNCham][kNPlane];
  double y[kNCham][kNPlane];
  double xx[kNCham][kNPlane];
  double yy[kNCham][kNPlane];
  double xy[kNCham][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));

  double dxdz[kNCham];

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

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

  double f = 0;
  for (int c = kDC1; c < kNCham; ++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);

      //cout << "Correlation coeff for " << c << p << ": " << corr << endl;
      
      f += corr;
    }
  }
  return f;
}

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

double AlignChamZ::Fit(const TrkInfo_t& ti, double resid[kNCham][kNPlane],
              double dxdz[kNCham])
{
  int nClust = 0;
  double u[kNCham * kNPlane];
  double w[kNCham * kNPlane];
  double z[kNCham * kNPlane];
  double cosA[kNCham * kNPlane];
  double sinA[kNCham * kNPlane];
  double lambdaX[kNCham * kNPlane];
  double lambdaY[kNCham * kNPlane];
  short chamber[kNCham * kNPlane];
  
  for (int cham = 0; cham < kNCham; ++cham) {
    GChamGeo& geo = GMIPPGeo::Instance().Cham(cham);
    for (int p = 0; p < kNPlane; ++p) {
      resid[cham][p] = 0;
      if (ti.fWeight[cham][p] < 1e-3) continue;
      geo.WireToU(p + 1, ti.fWire[cham][p], u[nClust]);
      w[nClust] = ti.fWeight[cham][p];
      z[nClust] = geo.Z(p + 1);
      cosA[nClust] = geo.GAngleCos(p + 1);
      sinA[nClust] = geo.GAngleSin(p + 1);
      lambdaX[nClust] = fLambdaX[cham][p];
      lambdaY[nClust] = fLambdaY[cham][p];
      chamber[nClust] = cham * kNPlane + p;
      nClust++;
    }
  }

  double par[5];
  FitCurvedTrack(nClust, u, w, z, cosA, sinA, fZ0, lambdaX, lambdaY,
                 par, 0);

  // Reinitialize lambda's based on the first-try fit
  double lx[kNCham][kNPlane], ly[kNCham][kNPlane];
  CalcLambda(fZ0, lx, ly, par);
  for (int i = 0; i < nClust; ++i) {
    int cham =chamber[i] / kNPlane;
    int pl = chamber[i] % kNPlane;
    lambdaX[i] = lx[cham][pl];
    lambdaY[i] = ly[cham][pl];
  }
  FitCurvedTrack(nClust, u, w, z, cosA, sinA, fZ0, lambdaX, lambdaY,
                 par, 0);
  double mom = par[4];
  SwimMIPP swim;
  swim.Swim(fZ0, par);
  for (int c = kDC1; c < kNCham; ++c) {
    double p[3];
    swim.GetPatZ(GMIPPGeo::Instance().Cham(c).ZGas(), p);
    dxdz[c] = p[0] / p[2];
  }

  // Loop through every cluster, exclude it from the fit and compute
  // residual
  for (int i = 0; i < nClust; ++i) {
    double saveW = w[i];
    w[i] = 0;
    FitCurvedTrack(nClust, u, w, z, cosA, sinA, fZ0, lambdaX, lambdaY,
           par, 0);
    w[i] = saveW;

    int c = chamber[i] / kNPlane;
    int p = chamber[i] % kNPlane;
    GChamGeo& geo = GMIPPGeo::Instance().Cham(c);

    double dz = geo.Z(p + 1) - fZ0;
    double x = par[0] + dz * par[2] + par[4] * lambdaX[i];
    double y = par[1] + dz * par[3] + par[4] * lambdaY[i];
    double ufit;
    geo.XYToU(p + 1, x, y, ufit);
    resid[c][p] = ufit - u[i];
  }

  return mom;
}

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

static void gsFCN(int& /*npar*/, double* /*gin*/, double& f,
                  double* par, int /*iflag*/)
{
  for (int i = 0; i < 12; ++i) {
    if (!isnormal(par[i]) && par[i] != 0) {
      f = 0;
      cout << "Abnormal parameters in TMinuit" << endl;
      return;
    }
  }
  cout << "Minimizing with " << par[7] << endl;

  // Setup geometry based on the new constants
  for (int i = 0; i < 6; ++i) {
    GMIPPGeo::Instance().Cham(i + kDC1).SetZGas(par[i]);
  }
  BFMagnet::JGG().SetCenterPos(par[6], par[7], par[8]);
  BFMagnet::Rosie().SetCenterPos(par[9], par[10], par[11]);
  
  // Compute the function we want to minimize
  f = gModule->ComputeF();
}

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

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

  cout << "Doing minimization with " << fTrackList.size() << " tracks"
       << endl;

  double lx[9][4], ly[9][4];
  CalcLambda(fZ0, lx, ly);
  for (int i = kDC1; i < 9; ++i) {
    cout << "Lambda X " << GChamGeo::Name(i) << ":";
    for (int j = 0; j < 4; ++j) cout << lx[i][j] << " ";
    cout << "\nLambda Y " << GChamGeo::Name(i) << ":";
    for (int j = 0; j < 4; ++j) cout << ly[i][j] << " ";
    cout << endl;
  }

  gModule = this;
  
  TMinuit minu(12);
  minu.SetFCN(gsFCN);

  char name[256];
  int ierr;
  double saveZ[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);
    saveZ[i] = geo.ZGas();
  }
  BFMagnet& jgg = BFMagnet::JGG();
  minu.mnparm(6, "JGGX", jgg.CenterPos()[0], 1, 0, 0, ierr);
  minu.mnparm(7, "JGGY", jgg.CenterPos()[1], 1, 0, 0, ierr);
  minu.mnparm(8, "JGGZ", jgg.CenterPos()[2], 1, 0, 0, ierr);

  BFMagnet& rosie = BFMagnet::Rosie();
  minu.mnparm(9, "RosieX", rosie.CenterPos()[0], 1, 0, 0, ierr);
  minu.mnparm(10, "RosieY", rosie.CenterPos()[1], 1, 0, 0, ierr);
  minu.mnparm(11, "RosieZ", rosie.CenterPos()[2], 1, 0, 0, ierr);

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

  minu.Migrad();

  MakeTree();

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

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

void AlignChamZ::MakeTree()
{
  double mom;
  double resid[kNCham][kNPlane];
  double  u[kNCham][kNPlane];
  double dxdz[kNCham];

  TTree* tree = new TTree("altree", "Post-alignment residuals");
  tree->Branch("mom", &mom, "mom/D");
  tree->Branch("resid", &resid, "resid[9][4]/D");
  tree->Branch("u", &u, "u[9][4]/D");
  tree->Branch("dxdz", &dxdz, "dxdz[9]/D");

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

    mom = Fit(ti, resid, dxdz);
    for (int c = 0; c < kNCham; ++c) {
      GChamGeo& geo = GMIPPGeo::Instance().Cham(c);
      for (int p = 0; p < kNPlane; ++p) {
    if (ti.fWeight[c][p] < 1e-3) u[c][p] = 1e6;
    else geo.WireToU(p + 1, ti.fWire[c][p], u[c][p]);
      }
    }
    tree->Fill();
  }

  tree->Write();
  delete tree;
}

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

void AlignChamZ::Update(const CfgConfig& c)
{
  c("Debug").Get(fDebug);
  c("Folder").Get(fFolder);
  c("LoadAlignment").Get(fLoadAlignment);
  c("MinNClust").Get(fMinNClust);
  fFixParam.clear();
  try { c("FixParam").Get(fFixParam); }
  catch (...) {
    int i;
    c("FixParam").Get(i);
    fFixParam.push_back(i);
  }
  c("UseCorr").Get(fUseCorr);

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

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

---