AlignTPCDrift2.cxx

Go to the documentation of this file.

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

#include "Alignment/AlignTPCDrift2.h"

#include "Bfield/bfield.h"
#include "Bfield/BFMagnet.h"
#include "Config/CfgConfig.h"
#include "Config/CfgParam.h"
#include "EventDataModel/EDMEventHandle.h"
#include "Geometry/GTPCGeo.h"
#include "RecoBase/RBKTrack.h"
#include "RecoBase/RBKVertex.h"
#include "RecoBase/RBViewHit.h"
#include "Swimmer/SwimDrift.h"
#include "Swimmer/SwimDriftMap.h"
#include "Swimmer/SwimMIPP.h"
#include "TPCRecoJP/TPCRHit.h"
#include "TrkRBase/TrkRUtil.h"

#include <TMinuit.h>
#include <cmath>

MODULE_DECL(AlignTPCDrift2);

using namespace std;

AlignTPCDrift2* AlignTPCDrift2::gModule = 0;

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

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

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

AlignTPCDrift2::~AlignTPCDrift2()
{
}

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

JobCResult AlignTPCDrift2::Ana(const EDMEventHandle& evt)
{
  vector<const RBKVertex*> vtx(0);

  evt.Cal();
  try { evt.Reco().Get(fFolder.c_str(), vtx); }
  catch (...) { return kFailed; }
  if (vtx.empty()) { return kFailed;}

  TrkInfo_t ti;
  ti.fSwim = 0;
  GTPCGeo& tpcgeo = GMIPPGeo::Instance().TPC();
  
  int nGood = 0;
  for (unsigned int i = 0; i < vtx.size(); ++i) {
    const RBKVertex* v = vtx[i];
    if (fRequireBeamTrk && !v->InTrack()) continue;
    int nTrk = v->NTracks() + (v->InTrack() ? 1 : 0);
    if (nTrk < fMinTrkVtx) continue;

    for (int j = 0; j < v->NTracks(); ++j) {
      const RBKTrack* t = v->OutTrack(j);

      if (t->NSpacePoints() < fMinTPCHits ||
      t->NSpacePoints() > fMaxTPCHits ||
      t->NViewHits() < fMinChamHits ||
      fabs(t->QoP()) > fMinP ||
      t->Time() < fMinTrkTime ||
      t->Time() > fMaxTrkTime ||
      t->GoF() < fMinGoF ||
      fabs(t->QoP()) < 5e-3) {
    if (fDebug) {
      cout << "Skipping track: nSpacePts=" << t->NSpacePoints()
           << " nViewHits=" << t->NViewHits()
           << " QoP=" << t->QoP()
           << " T=" << t->Time()
           << " GoF=" << t->GoF()
           << endl;
    }
    continue;
      }
      
      ti.fNTPCHit = t->NSpacePoints(); 
      if (ti.fNTPCHit > 128) ti.fNTPCHit = 128;
      for (int j = 0; j < ti.fNTPCHit; ++j) {
    const TPCRHit* hit = 
      & dynamic_cast<const TPCRHit&>(t->SpacePoint(j));
    assert (hit);
    ti.fPadX[j] = tpcgeo.PadColToX(hit->PadCol());
    ti.fPadZ[j] = tpcgeo.PadRowToZ(hit->PadRow());
    ti.fDriftT[j] = hit->TDC() * GTPCGeo::kUsecPerBucket + 
      tpcgeo.TOffset();
      }
      
      ti.fNChamHit = t->NViewHits();
      if (ti.fNChamHit > 36) ti.fNChamHit = 36;
      for (int j = 0; j < ti.fNChamHit; ++j) {
    const RBViewHit& hit = t->ViewHit(j);
    ti.fChamZ[j] = hit.Z();
    ti.fU[j] = hit.U();
    ti.fW[j] = 1 / (hit.USig() * hit.USig());
    ti.fCosA[j] = hit.CosTheta();
    ti.fSinA[j] = hit.SinTheta();
      }
      ti.fZ0 = t->VtxZ();
      ti.fPar[0] = t->VtxX();
      ti.fPar[1] = t->VtxY();
      ti.fPar[2] = t->VtxDXDZ();
      ti.fPar[3] = t->VtxDYDZ();
      ti.fPar[4] = t->QoP();

      double chi2 = ComputeTrkChi2(ti);
      if (chi2 > 1e6)  continue;

      if (fDebug) {
    printf("Using track: N=%d/%d q/p=%.1e vtx=(%.1f,%.1f;%.1e,%.1e) chi2=%.1f\n",
           t->NSpacePoints(), t->NViewHits(), t->QoP(),
           t->VtxX(), t->VtxY(), t->VtxDXDZ(), t->VtxDYDZ(),
           t->Chi2());
    printf("My fit: q/p=%.1e vtx=(%.1f,%.1f;%.1e,%.1e) chi2=%.1f\n",
           ti.fPar[4], ti.fPar[0], ti.fPar[1], ti.fPar[2], ti.fPar[3],
           chi2);
      }
      

      fTrkInfo.push_back(ti);
      ++nGood;
    }
  }

  if (fDebug) cout << "  picked " << nGood << " tracks" << endl;

  return (nGood > 0);
}

//......................................................................
  
void AlignTPCDrift2::NewRun(int run, int subrun)
{
  GMIPPGeo::Instance(run);
  config_bfield(run, subrun);
  SwimDriftMap::Instance().Init(run);
  SwimDrift::Instance().SetUseDriftMap(1);
}

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

void AlignTPCDrift2::EndRun(int /*run*/, int /*subrun*/)
{
  if (fTrkInfo.size() < 1) {
    cout << "Can't do minimization with " << fTrkInfo.size() << " tracks"
     << endl;
    return;
  }
  cout << "AlignTPCDrift2: doing minimization with " << fTrkInfo.size()
       << " tracks" << endl;
  gModule = this;

  TMinuit* minu = new TMinuit(7);
  minu->SetMaxIterations(1000000);
  minu->SetFCN(AlignTPCDrift2::MinFCN);
  char name[32];
  for (int i = 0; i < 3; ++i) {
    sprintf(name, "V%c Scale", 'x' + i);
    minu->DefineParameter(i, name, 1, 0.01, 0.1, 3);
  }
  for (int i = 0; i < 3; ++i) {
    sprintf(name, "V%c LinScale", 'x' + i);
    minu->DefineParameter(i + 3, name, 0, 0.01, -1, 1);
  }
  minu->DefineParameter(6, "Bx Bz rotation", 0.007, 1e-3, -2e-2, 2e-2);

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

  minu->Migrad();

  delete minu;
}

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

void AlignTPCDrift2::MinFCN(int& /*npar*/, double* /*gin*/, double& f,
               double* par, int /*iflag*/)
{
  SwimDriftMap::Instance().SetScale(par[0], par[1], par[2],
                    par[3], par[4], par[5]);
  SwimDriftMap::Instance().SetBxBzRotation(par[6]);

  f = 0;
  for (unsigned int i = 0; i < gModule->fTrkInfo.size(); ++i)
    f += gModule->ComputeTrkChi2(gModule->fTrkInfo[i]);

  cout << "par=(" << par[0] << ", " << par[1] << ", " << par[2]
       << "; " << par[3] << ", " << par[4] << ", " << par[5]
       << ", " << par[6] << ") f=" << f << endl;

}

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

void AlignTPCDrift2::Update(const CfgConfig& c)
{
  c("Debug").Get(fDebug);
  c("Folder").Get(fFolder);
  fFixParam.clear();
  try { c("FixParam").Get(fFixParam); }
  catch (...) { 
    int val;
    c("FixParam").Get(val);
    fFixParam.push_back(val);
  }
  c("MinTrkVtx").Get(fMinTrkVtx);
  c("RequireBeamTrk").Get(fRequireBeamTrk);
  c("MinTPCHits").Get(fMinTPCHits);
  c("MaxTPCHits").Get(fMaxTPCHits);
  c("MinChamHits").Get(fMinChamHits);
  c("MinP").Get(fMinP);
  fMinP = (fMinP > 0) ? 1 / fMinP : 1e6; // Convert to 1/p
  c("MinGoF").Get(fMinGoF);
  c("TPCSigmaX").Get(fTPCSigmaX);
  fTPCSigmaX = 1 / (fTPCSigmaX * fTPCSigmaX);
  c("TPCSigmaY").Get(fTPCSigmaY);
  fTPCSigmaX = 1 / (fTPCSigmaY * fTPCSigmaY);

  c("MinTrkTime").Get(fMinTrkTime);
  c("MaxTrkTime").Get(fMaxTrkTime);

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

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

double AlignTPCDrift2::ComputeTrkChi2(TrkInfo_t& ti)
{
  double u[512];
  double w[512];
  double z[512];
  double cosA[512];
  double sinA[512];
  double lamX[512];
  double lamY[512];

  if (!ti.fSwim) {
    ti.fSwim = new SwimMIPP;
    ti.fSwim->Swim(ti.fZ0, ti.fPar);
  }
  double p = 1 / ti.fPar[4];

  double padPos[3];
  double cavePos[3];
  double hitPos[3];
  SwimDrift& drift = SwimDrift::Instance();
  GTPCGeo& tpcgeo = GMIPPGeo::Instance().TPC();

  for (int i = 0; i < ti.fNTPCHit; ++i) {
    // Convert hit position from pad plane to cave
    padPos[0] = ti.fPadX[i];
    padPos[1] = 0;
    padPos[2] = ti.fPadZ[i];
    tpcgeo.PPLNToCAVE(padPos, cavePos);    

    // Drift hit to figure out its actual position
    drift.CorrectDistortion(cavePos, hitPos, ti.fDriftT[i]);

    z[i * 2 + 1] = z[i * 2] = hitPos[2];
    u[i * 2] = hitPos[0];
    u[i * 2 + 1] = hitPos[1];
    w[i * 2]  = fTPCSigmaX; // sigma's have been converted to weights
    w[i * 2 + 1]  = fTPCSigmaY;

    // Compute template coefficients
    if (!ti.fSwim->GetXY(hitPos[2], &hitPos[0], &hitPos[1])) return 1e6;
    lamX[i * 2 + 1] = lamX[i * 2] = 
      p * (hitPos[0] - ti.fPar[0] - ti.fPar[2] * (hitPos[2] - ti.fZ0));
    lamY[i * 2 + 1] = lamY[i * 2] = 
      p * (hitPos[1] - ti.fPar[1] - ti.fPar[3] * (hitPos[2] - ti.fZ0));

    cosA[i * 2] = 1;
    sinA[i * 2] = 0;
    cosA[i * 2 + 1] = 0;
    sinA[i * 2 + 1] = -1;
  }


  int i = ti.fNTPCHit * 2;
  for (int j = 0; j < ti.fNChamHit; ++i, ++j) {
    u[i]    = ti.fU[j];
    w[i]    = ti.fW[j];
    cosA[i] = ti.fCosA[j];
    sinA[i] = ti.fSinA[j];
    z[i]    = ti.fChamZ[j];
    if (!ti.fSwim->GetXY(z[i], &hitPos[0], &hitPos[1])) return 1e6;
    lamX[i] = 
      p * (hitPos[0] - ti.fPar[0] - ti.fPar[2] * (z[i] - ti.fZ0));
    lamY[i] = 
      p * (hitPos[1] - ti.fPar[1] - ti.fPar[3] * (z[i] - ti.fZ0));
  }

//   for (int j = 0; j < i; ++j) {
//     printf("z=%.1f u=%.1f cos=%.1f sin=%.1f lamx=%.1f lamy=%.1f\n", 
//         z[j], u[j], cosA[j], sinA[j], lamX[j], lamY[j]);
//   }

  double par[5];
  double chi2 =
    TrkRUtil::FitCurvedTrack(i, u, w, z, cosA, sinA, ti.fZ0, lamX, lamY, 
                 par, 0);

  // Recompute template coefficients
  SwimMIPP swim;
  swim.Swim(ti.fZ0, par);
  p = 1 / par[4];

  for (int j = 0; j < i; ++j) {
    if (!ti.fSwim->GetXY(z[i], &hitPos[0], &hitPos[1])) return 1e6;
    lamX[i] = 
      p * (hitPos[0] - ti.fPar[0] - ti.fPar[2] * (z[i] - ti.fZ0));
    lamY[i] = 
      p * (hitPos[1] - ti.fPar[1] - ti.fPar[3] * (z[i] - ti.fZ0));    
  }


  chi2 =
    TrkRUtil::FitCurvedTrack(i, u, w, z, cosA, sinA, ti.fZ0, lamX, lamY, 
                 par, 0);

  return chi2;
}

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

---