AlignChamW0.cxx

Go to the documentation of this file.

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

#include "Alignment/AlignChamW0.h"

#include "Bfield/BFMagnet.h"
#include "Bfield/bfield.h"
#include "Config/CfgConfig.h"
#include "Config/CfgParam.h"
#include "EventDataModel/EDMEventHandle.h"
#include "EventDataModel/EDMDataBucket.h"
#include "EventDataModel/EDMEventHeader.h"
#include "Geometry/GTargetGeo.h"
#include "NumericalMethods/NMFit.h"
#include "RunInfo/RIRunConfig.h"
#include "TrkRBase/TrkCand.h"
#include "TrkRBase/WireCross.h"
#include "TrkRBase/WireClust.h"
#include "TrkRBase/TrkRUtil.h"
#include "TrkRBase/TrkChamGeo.h"

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

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

MODULE_DECL(AlignChamW0);
ClassImp(AlResidNt);

using namespace std;
using namespace TrkRUtil;

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

int AlignChamW0::fsCham = 0;

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

AlignChamW0::AlignChamW0(const char* version) :
  JobCModule("AlignChamW0"),
  fFolder("")
{
  this->SetCfgVersion(version);
  CheckInit();
}

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

AlignChamW0::~AlignChamW0()
{
  fTrackList.clear();
}

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

JobCResult AlignChamW0::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 (fFieldOn && 
    (fabs(trk->QoverP()) < fMaxMom || fabs(trk->QoverP()) > fMinMom)) {
      continue;
    }
    memset(&ti, 0, sizeof(TrkInfo_t));
    
    int nClu = 0;
    bool haveBC = false;
    for (int j = 0; j < trk->NClust(); ++j) {
      const WireClust* wc = trk->Clust(j);
      if (wc->NWire() > 1) continue;

      int cham = wc->Chamber();
      int pl = wc->Plane() - 1;

      if (cham > kBC3) nClu++;
      else {
    if (!fUseBC) continue;
    haveBC = true;
      }

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

    ti.fPoQ = trk->PQ();
    ti.fChi2 = trk->Chi2() / trk->NDF();
    fTrackList.push_back(ti);
    if (haveBC) fNTrackWithBC++;
  }
  return kPassed;
}

//----------------------------------------------------------------------
  
void AlignChamW0::NewRun(int run, int subrun)
{
  fNTrackWithBC = 0;
  fTrackList.clear();
  for (int i = 0; i < kNCham; ++i) {
    for (int j = 0; j < kNPlane; ++j) fNClust[i][j] = 0;
  }
  GMIPPGeo::Instance(run);
  if (fLoadAlignment) TrkChamGeo::Instance().Init(run);  
  config_bfield(run, subrun);
  if (BFMagnet::JGG().ByAtCenter() > 0.5 && 
      BFMagnet::Rosie().ByAtCenter() < -0.5) {
    fFieldOn = true;
  }
  else fFieldOn = false;

}

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

void AlignChamW0::EndRun(int run, int /*subrun*/)
{
  if ((int) fTrackList.size() < fMinTracks) {
    cout << "Don't have enough tracks (" << fTrackList.size() << ") "
     << " to do chamber alignment" << endl;
    return;
  }

  // Choose z0 to be just upstream of DC1, since a lot of tracks
  // will only have chambers 1-6, w/o BC's
  fZ0 = GMIPPGeo::Instance(run).Target().CenterPos()[2];
  TrkCand::Init(fZ0);

  double saveBy = BFMagnet::JGG().ByAtCenter();
  
  // First align w/o shifting BC positions, then allow to change 
  // BC positions
  if (!RoughAlign(false)) {
    cout << "AlignChamW0: alignment failed" << endl;
    return;
  }

  TrkChamGeo& db = TrkChamGeo::Instance();
  db.SetRunRange(run, run);

  // Save alignment constants with BC's fixed
  for (int c = kBC1; c < kNCham; ++c) db.SetCham(c);
  db.WriteToDB(TrkChamGeo::kBCFixed, 0770);

  // The other steps (adjusting JGG field scale and BC alignment)
  // are only possible if we have enough BC tracks
  if (fNTrackWithBC >= fMinTracksBC) {
    if (fAdjustBField  && fabs(saveBy) > 0.5) {
      // Find the new minimum for B-field, and if the field was changed,
      // then we make sure not to finish the loop
      AdjustBField();
      
      if (!RoughAlign(false)) {
    cout << "AlignChamW0: alignment failed after JGG scale was changed"
         << endl;
    return;
      }
      
      cout << "AlignChamW0 run " << run << " beam: " 
       << RIRunConfig::Instance(run).Momentum() << " By=" << saveBy
       << " now " << BFMagnet::JGG().ByAtCenter()
       << " ratio=" << BFMagnet::JGG().ByAtCenter() / saveBy
       << endl;
      
      // Save alignment constants having tweaked JGG field
      for (int c = kBC1; c < kNCham; ++c) db.SetCham(c);
      db.WriteToDB(TrkChamGeo::kPostJGGScale, 0770);
    }
    
    if (fUseBC) {
      if (!RoughAlign(true)) {
    cout << "AlignChamW0: couldn't do alignment with BC's" << endl;
      }
      else {
    for (int c = kBC1; c < kNCham; ++c) db.SetCham(c);
    db.WriteToDB(TrkChamGeo::kAllCham, 0777);
      }      
    }
  }
  else cout << "Number of tracks with BC: " << fNTrackWithBC << "/"
        << fMinTracksBC << endl;


  AnaResid(0);

}

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

void AlignChamW0::Update(const CfgConfig& c)
{
  c("LoadAlignment").Get(fLoadAlignment);
  c("MaxAlignAttempt").Get(fMaxAlignAttempt);
  c("MinNClust").Get(fMinNClust);
  c("UseBC").Get(fUseBC);
  c("Folder").Get(fFolder);
  c("OkShift").Get(fOkShift);
  c("MinTracks").Get(fMinTracks);
  c("MinTracksBC").Get(fMinTracksBC);
  c("AdjustBField").Get(fAdjustBField);
  c("MinMom").Get(fMinMom);
  c("MaxMom").Get(fMaxMom);

  assert(fMinMom > 0);
  assert(fMaxMom > 0);
  fMinMom = 1 / fMinMom;
  fMaxMom = 1 / fMaxMom;

  fIsInit = true;
}

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

double AlignChamW0::Fit(AlignChamW0::TrkInfo_t& trk, 
             float* poq /* = 0*/)
{
  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) {
      if (trk.fWeight[cham][p] < 1e-3) continue;
      geo.WireToU(p + 1, trk.fWire[cham][p], u[nClust]);
      w[nClust] = trk.fWeight[cham][p];
      z[nClust] = geo.Z(p + 1);
      cosA[nClust] = geo.GAngleCos(p + 1);
      sinA[nClust] = geo.GAngleSin(p + 1);
      lambdaX[nClust] = TrkCand::LambdaX(cham, p);
      lambdaY[nClust] = TrkCand::LambdaY(cham, p);
      chamber[nClust] = cham * kNPlane + p;
      nClust++;
    }
  }

  double par[5];

  // Make first-pass fit
  FitCurvedTrack(nClust, u, w, z, cosA, sinA, fZ0, lambdaX, lambdaY,
         par, 0);

  // Initialize lambda's based on the fit
  double lamX[kNCham][kNPlane], lamY[kNCham][kNPlane];
  CalcLambda(fZ0, lamX, lamY, par, chamber[0] / 4, chamber[nClust - 1] / 4);

  // Adjust lambda for each point based on the rough track parameters
  for (int i = 0; i < nClust; ++i) {
    lambdaX[i] = lamX[chamber[i] / kNPlane][chamber[i] % kNPlane];;
    lambdaY[i] = lamY[chamber[i] / kNPlane][chamber[i] % kNPlane];
  }

  // Return chi squared of the fit using better-evaluated lambda's
  double chi2 = FitCurvedTrack(nClust, u, w, z, cosA, sinA,
                   fZ0, lambdaX, lambdaY, par, 0);
  if (!poq) return chi2;
  
  chi2 /= (double) (nClust - 5);
  if (fabs(par[4]) > 1e-3) *poq = 1.0 / par[4];
  else *poq = 1e3;

  return chi2;
}

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

void AlignChamW0::AnaResid(int i)
{
  char name[32];
  sprintf(name, "resid%d", i);
  TTree* tree = new TTree(name, "Chamber Track Residual Tree");
  AlResidNt* nt = new AlResidNt;
  tree->Branch("fNt", "AlResidNt", &nt);

  for (unsigned int i = 0; i < fTrackList.size(); ++i) {
    Fit(fTrackList[i], nt);
    tree->Fill();
  }

  tree->Write();
  delete tree;
  delete nt;
}

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

void AlignChamW0::Fit(AlignChamW0::TrkInfo_t& trk, AlResidNt* nt,
               int minCham /* = 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];

  nt->nclust = 0;

  for (int cham = minCham; cham < kNCham; ++cham) {
    GChamGeo& geo = GMIPPGeo::Instance().Cham(cham);
    for (int p = 0; p < kNPlane; ++p) {
      if (trk.fWeight[cham][p] < 1e-3) {
    nt->u[cham][p] = 1e6;
      }
      else {
    geo.WireToU(p + 1, trk.fWire[cham][p], u[nt->nclust]);
    nt->u[cham][p] = u[nt->nclust];
    w[nt->nclust] = trk.fWeight[cham][p];
    z[nt->nclust] = geo.Z(p + 1);
    cosA[nt->nclust] = geo.GAngleCos(p + 1);
    sinA[nt->nclust] = geo.GAngleSin(p + 1);
    lambdaX[nt->nclust] = TrkCand::LambdaX(cham, p);
    lambdaY[nt->nclust] = TrkCand::LambdaY(cham, p);
    chamber[nt->nclust] = cham * kNPlane + p;
    nt->nclust++;
      }
    }
  }
   
  double par[5];
  
  // Make first-pass fit
  FitCurvedTrack(nt->nclust, u, w, z, cosA, sinA, 
         fZ0, lambdaX, lambdaY, par, 0);
  
  // Initialize lambda's based on the fit
  double lamX[kNCham][kNPlane], lamY[kNCham][kNPlane]; 
  CalcLambda(fZ0, lamX, lamY, par, 0, kPWC6);
  
  // Adjust lambda for each point based on the rough track parameters
  for (int i = 0; i < nt->nclust; ++i) {
    lambdaX[i] = lamX[chamber[i] / kNPlane][chamber[i] % kNPlane];
    lambdaY[i] = lamY[chamber[i] / kNPlane][chamber[i] % kNPlane];
  }
  
  // Return chi squared of the fit using better-evaluated lambda's
  nt->chi2 = FitCurvedTrack(nt->nclust, u, w, z, cosA, sinA,
                fZ0, lambdaX, lambdaY, par, 0);
  
  nt->poq = (fabs(par[4]) > 1e-3) ? (1 / par[4]) : 1e3;
  
  double x, y;
  for (int cham = 0; cham < kNCham; ++cham) {
    GChamGeo& geo = GMIPPGeo::Instance().Cham(cham);

    for (int p = 0; p < kNPlane; ++p) {
      x = par[0] + (geo.Z(p + 1) - fZ0) * par[2] + par[4] * lamX[cham][p];
      y = par[1] + (geo.Z(p + 1) - fZ0) * par[3] + par[4] * lamY[cham][p];
    
      nt->ufit[cham][p] = x * geo.GAngleCos(p + 1) - 
    y * geo.GAngleSin(p + 1);
      nt->uperp[cham][p] = x * geo.GAngleSin(p + 1) + 
    y * geo.GAngleCos(p + 1);
    }
  }

  for (int c = 0; c < kNCham; ++c) {
    for (int p = 0; p < kNPlane; ++p) {
      nt->poqrm[c][p] = 1e6;
      nt->chi2rm[c][p] = 1e6;
      nt->resid[c][p] = 1e6;
    }
  }
  
  double wCopy[kNCham * kNPlane];
  // Exclude every hit one by one and record residual
  for (int i = 0; i < nt->nclust; ++i) {
    int cham = chamber[i] / kNPlane;
    int p = chamber[i] % kNPlane;
    GChamGeo& geo = GMIPPGeo::Instance().Cham(cham);

    memcpy(wCopy, w, sizeof(w));
    wCopy[i] = 0;

    nt->chi2rm[cham][p] = FitCurvedTrack(nt->nclust, u, wCopy, z, cosA, sinA,
                     fZ0, lambdaX, lambdaY, par, 0);
    
    nt->poqrm[cham][p] = (fabs(par[4]) > 1e-3) ? (1 / par[4]) : 1e3;

    x = par[0] + (geo.Z(p + 1) - fZ0) * par[2] + par[4] * lamX[cham][p];
    y = par[1] + (geo.Z(p + 1) - fZ0) * par[3] + par[4] * lamY[cham][p];
    nt->resid[cham][p] = x * geo.GAngleCos(p + 1) - 
      y * geo.GAngleSin(p + 1) - u[i];
  }
}

//----------------------------------------------------------------------
/// Loop through beam tracks, and shift individual wire planes around
/// by 1/2 of the average unbiased residual.
///
bool AlignChamW0::RoughAlign(bool withBC)
{
  TH1F* resH[kNCham][kNPlane];
  char str[256];
  for (int c = 0; c < kNCham; ++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 sp. bins
      resH[c][p]->SetXTitle("Residual (wire spacing)");
    }
  }

  AlResidNt nt;
  TDirectory* dir = gDirectory;

  for (int at = 0; ; ++at) {
    // Fail this run if we can't quickly converge
    if (at >= fMaxAlignAttempt) return false;
    int ntrkBC = 0;

    cout << "Chamber align attempt " << at << " with "
     << fTrackList.size() << " tracks" << endl;
    
    if (withBC) sprintf(str, "attemptbc%d", at);
    else sprintf(str, "attempt%d", at);
    gDirectory->mkdir(str);
    gDirectory->cd(str);
    
    for (unsigned int i = 0; i < fTrackList.size(); ++i) {
      // Check that BC1 is present
      TrkInfo_t& trk = fTrackList[i];
      if (trk.fWeight[0][0] > 0 ||
      trk.fWeight[0][1] > 0 ||
      trk.fWeight[0][2] > 0 ||
      trk.fWeight[0][3] > 0) ntrkBC++;
      
      Fit(trk, &nt);
      for (int c = 0; c < kNCham; ++c) {
    for (int p = 0; p < kNPlane; ++p) {
      if (nt.resid[c][p] > 1e5) continue;
      resH[c][p]->Fill(nt.resid[c][p] / GChamGeo::Pitch(c));
    }
      }
    }

    if (withBC && (ntrkBC < fMinTracksBC)) {
      cout << "Can't do alignment with " << ntrkBC << " tracks with "
       << " BC info" << endl;
      return false;
    }

    // Now extract means of the histograms
    double maxRes = 0;
    int maxResP = -1;
    double res[kNCham][kNPlane];
    for (int c = 0; c < kNCham; ++c) {
      GChamGeo& geo = GMIPPGeo::Instance().Cham(c);
      for (int p = 0; p < kNPlane; ++p) {
    if (resH[c][p]->Integral() < 2) {
      res[c][p] = 0;
      continue;
    }
    
    res[c][p] = resH[c][p]->GetMean(); 
    resH[c][p]->Write();
    resH[c][p]->Reset();

    if (withBC || c > kBC3) {
      if (fabs(res[c][p]) > fabs(maxRes)) {
        maxRes = res[c][p];
        maxResP = c * kNPlane + p;
      }

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

    if (maxResP < 0) {
      cout << "Found no non-zero residuals?" << endl;
      return false;
    }

    cout << "Maximum residual in " << GChamGeo::Name(maxResP / 4)
     << " P" << maxResP % 4 + 1 << ": " << maxRes << endl;
  
    // Preen track list to throw out those
    if (fabs(maxRes) < 0.25 && PreenTrackList(res, withBC)) continue;

    if (fabs(maxRes) < fOkShift) break;
  }
  for (int c = 0; c < kNCham; ++c) {
    for (int p = 0; p < kNPlane; ++p) delete resH[c][p];
  }
  return true;
}

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

int AlignChamW0::PreenTrackList(double res[kNCham][kNPlane], bool withBC)
{
  vector<TrkInfo_t>::iterator itr(fTrackList.begin());

  AlResidNt nt;
  int removed = 0;

  while (itr != fTrackList.end()) {
    TrkInfo_t& trk = *itr;
    Fit(trk, &nt);
    
    int c = withBC ? kBC1 : kDC1;
    int p;

    for ( ; c < kNCham; ++c) {
      for (p = 0; p < kNPlane; ++p) {
    if (trk.fWeight[c][p] <= 0) continue;
    double r = fabs(nt.resid[c][p] / GChamGeo::Pitch(c) - res[c][p]);
    if (r > 1.5 + fabs(res[c][p])) break;
      }
      if (p != kNPlane) break;
    }
    if (c != kNCham) {
      itr = fTrackList.erase(itr);
      removed++;
    }
    else itr++;
  }
  return removed;
}

//----------------------------------------------------------------------
/// Calculates the sum of chi squared for current B-field, 
/// +/-0.5%, and +/-1%. Based on those 5 points, determines the 
/// new best setting for JGG field.
/// Returns true if the field was changed by more than 0.25%
///
bool AlignChamW0::AdjustBField()
{
  double bNom[5];
  bNom[2] = BFMagnet::JGG().ByAtCenter();
  if (bNom[2] < 0.5) return false;

  bNom[0] = bNom[2] * 0.990;
  bNom[1] = bNom[2] * 0.995;
  bNom[3] = bNom[2] * 1.005;
  bNom[4] = bNom[2] * 1.010;

  double sum[5];

  for (int i = 0; i < 5; ++i) {
    BFMagnet::JGG().SetByAtCenter(bNom[i]);
    TrkCand::Init(fZ0);

    sum[i] = 0;
    for (unsigned int j = 0; j < fTrackList.size(); ++j) {
      sum[i] += Fit(fTrackList[j]);
    }
  }

  cout << "Chi^2 sums: ";
  for (int i = 0; i < 5; ++i) cout << sum[i] << " ";
  cout << endl;

  double par[3];
  NMFit::FitPol(2, 5, bNom, sum, 0, par);

  double minPos;
  if (par[2] > 0) { // We are guaranteed a minimum
    minPos = -par[1] / (2.0 * par[2]);
    cout << "Minimum is at " << minPos << endl;
    if (minPos < bNom[0]) minPos = bNom[0];
    else if (minPos > bNom[4]) minPos = bNom[4];
  }
  else {
    cout << "AlignChamW0::AdjustBField(): "
     << "we have a problem parameters are: " 
     << par[0] << ", " << par[1] << ", " << par[2] << endl;
    double slope = 2 * par[2] * bNom[2] + par[1];
    if (slope > 0) minPos = bNom[3];
    else minPos = bNom[1];
  }

  BFMagnet::JGG().SetByAtCenter(minPos);
  return (fabs(minPos / bNom[2] - 1) > 2.5e-3);
}

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

---