AlignTPCDrift.cxx

Go to the documentation of this file.

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

#include "Alignment/AlignTPCDrift.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/SwimMIPP.h"
#include "TPCRecoJP/TPCRHit.h"
#include "TrkRBase/TrkRUtil.h"

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

MODULE_DECL(AlignTPCDrift);

using namespace std;

AlignTPCDrift* AlignTPCDrift::gModule = 0;

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

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

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

AlignTPCDrift::~AlignTPCDrift()
{
}

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

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

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

  TrkInfo_t ti;
  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) {
    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();
      if (fabs(ti.fPar[4]) < 1e-3) ti.fPar[4] = 1e-3;

      ti.fSwim = new SwimMIPP;
      ti.fSwim->Swim(ti.fZ0, ti.fPar);

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

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

  return (nGood > 0);
}

//......................................................................
  
void AlignTPCDrift::NewRun(int run, int subrun)
{
  GMIPPGeo::Instance(run);
  config_bfield(run, subrun);

  fTrigOffset = 0.1;
  fXOffset = 0;
  fYOffset = GMIPPGeo::Instance().TPC().YOffset();
  fZOffset = 0;

}

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

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

  TMinuit* minu = new TMinuit(15);
  minu->SetMaxIterations(1000000);
  minu->SetFCN(AlignTPCDrift::MinFCN);
  char name[32];
  int nPar = 0;
  for (int i = 0; i < 6; ++i) {
    sprintf(name, "X par %d", i);
    minu->DefineParameter(nPar, name, fInitParam[nPar], 
              0.3 * fabs(fInitParam[nPar]) + 0.001, -5, 5);
    ++nPar;
  }
  for (int i = 0; i < 4; ++i) {
    sprintf(name, "Y par %d", i);
    minu->DefineParameter(nPar, name, fInitParam[nPar], 
              0.3 * fabs(fInitParam[nPar]) + 0.001, -5, 5);
    ++nPar;
  }
  for (int i = 0; i < 5; ++i) {
    sprintf(name, "Z par %d", i);
    minu->DefineParameter(nPar, name, fInitParam[nPar], 
              0.3 * fabs(fInitParam[nPar]) + 0.001, -5, 5);
    ++nPar;
  }

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

  SwimDrift::Instance().SetUseDriftMap(2);
  minu->Migrad();

  delete minu;
}

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

void AlignTPCDrift::MinFCN(int& /*npar*/, double* /*gin*/, double& f,
               double* par, int /*iflag*/)
{
  SwimDrift::Instance().SetFitDrift(&par[0], &par[6], &par[6+4]);

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

}

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

void AlignTPCDrift::Update(const CfgConfig& c)
{
  c("Debug").Get(fDebug);
  c("Folder").Get(fFolder);
  c("InitParam").Get(fInitParam);
  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 AlignTPCDrift::ComputeTrkChi2(TrkInfo_t& ti)
{

  GTPCGeo& geo = GMIPPGeo::Instance().TPC();
  SwimDrift& drift = SwimDrift::Instance();
  double xPad[3], xGlob[3], xCorr[3];

  double chi2 = 0;
  for (int i = 0; i < ti.fNTPCHit; ++i) {
    xPad[0] = ti.fPadX[i];
    xPad[1] = 0;
    xPad[2] = ti.fPadZ[i];

    geo.PPLNToCAVE(xPad, xGlob);

    drift.CorrectDistortion(xGlob, xCorr, ti.fDriftT[i]);

    chi2 += ti.fSwim->DsqrPointToTrack(xCorr[0], xCorr[1], xCorr[2]);
  }

  return chi2;
}

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

---