AlignRotB.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// $Id: AlignRotB.cxx,v 1.1 2007/01/09 21:45:45 lebedev Exp $
//
// lebedev@physics.harvard.edu
////////////////////////////////////////////////////////////////////////

#include "Alignment/AlignRotB.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(AlignRotB);

using namespace std;
using namespace TrkRUtil;

static AlignRotB* gModule = 0;

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

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

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

AlignRotB::~AlignRotB()
{
}

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

JobCResult AlignRotB::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];
    if (!trk->SegBC()) continue; // We must have tracks with BC info

    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 (wc->NWire() > 1) continue;

      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 AlignRotB::NewRun(int run, int subrun)
{
  GMIPPGeo::Instance(run);
  config_bfield(run, subrun);
  TrkChamGeo::Instance().Init(run);
  fZ0 = GMIPPGeo::Instance().Cham(kBC3).ZGas();
}

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

double AlignRotB::ComputeF()
{
  CalcLambda(fZ0, fLambdaX, fLambdaY);
  
  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];
    }
    f += FitCurvedTrack(nClust, u, w, z, cosA, sinA, fZ0, lambdaX, lambdaY,
            par, 0);  
  }
  return f;
}

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

double AlignRotB::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, kDC1, kPWC6);
  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*/)
{
  // Setup field rotation
  Swimmer::SetFieldRotZ(par[0], par[1]);
  
  // Compute the function we want to minimize
  f = gModule->ComputeF();
}

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

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

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

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

  minu.DefineParameter(0, "JGG field z-rotation", 0, 1e-3, -0.02, 0.02);
  minu.DefineParameter(1, "Rosie field z-rotation", 0, 1e-3, -0.02, 0.02);

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

  minu.Migrad();

  //MakeTree();

  for (int i = 0; i < 2; ++i) {
    double par, perr;
    minu.GetParameter(i, par, perr);
    cout << ((i == 0) ? "JGG " : "Rosie ")
     << " field rotation: " << par*1e3 << " mrad, err: " << perr*1e3 << endl;
  }
}

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

void AlignRotB::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 AlignRotB::Update(const CfgConfig& c)
{
  c("Debug").Get(fDebug);
  c("Folder").Get(fFolder);
  c("MinNClust").Get(fMinNClust);
  fFixParam.clear();
  try { c("FixParam").Get(fFixParam); }
  catch (...) {
    int i;
    c("FixParam").Get(i);
    fFixParam.push_back(i);
  }

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

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

---