AlignTPCDrift2 Class Reference

#include <AlignTPCDrift2.h>

Inheritance diagram for AlignTPCDrift2:

JobCModule CfgObserver List of all members.

Public Member Functions

 AlignTPCDrift2 (const char *version)
 ~AlignTPCDrift2 ()
JobCResult Ana (const EDMEventHandle &evt)
void NewRun (int run, int subrun)
void EndRun (int run, int subrun)
void Update (const CfgConfig &c)

Static Public Member Functions

static void MinFCN (int &, double *, double &f, double *par, int)

Private Member Functions

double ComputeTrkChi2 (TrkInfo_t &trk)

Private Attributes

std::vector< TrkInfo_t > fTrkInfo
int fDebug
std::vector< int > fFixParam
std::string fFolder
int fMinTrkVtx
int fRequireBeamTrk
int fMinTPCHits
int fMaxTPCHits
int fMinChamHits
float fMinP
float fMinGoF
float fTPCSigmaX
float fTPCSigmaY
float fMinTrkTime
float fMaxTrkTime

Static Private Attributes

static AlignTPCDrift2 * gModule = 0

Classes

struct  TrkInfo_t

Detailed Description

Author:
lebedev@physics.harvard.edu Class which picks up global tracks and minimizes the sum of chi squared by scaling drift velocity vector components

Definition at line 21 of file AlignTPCDrift2.h.

Constructor & Destructor Documentation

AlignTPCDrift2::AlignTPCDrift2 ( const char *  version  ) 

Definition at line 35 of file AlignTPCDrift2.cxx.

References JobCModule::CheckInit(), and JobCModule::SetCfgVersion().

00035                                                   :
00036   JobCModule("AlignTPCDrift2")
00037 {
00038   // Register module to be updated with a given config version
00039   SetCfgVersion(version);
00040   // Make sure that configuration has been loaded
00041   CheckInit();
00042 }

AlignTPCDrift2::~AlignTPCDrift2 (  ) 

Definition at line 46 of file AlignTPCDrift2.cxx.

00047 {
00048 }

Member Function Documentation

JobCResult AlignTPCDrift2::Ana ( const EDMEventHandle &  evt  )  [virtual]

Reimplemented from JobCModule.

Definition at line 52 of file AlignTPCDrift2.cxx.

References EDMEventHandle::Cal(), ComputeTrkChi2(), AlignTPCDrift2::TrkInfo_t::fChamZ, AlignTPCDrift2::TrkInfo_t::fCosA, fDebug, AlignTPCDrift2::TrkInfo_t::fDriftT, fFolder, fMaxTPCHits, fMaxTrkTime, fMinChamHits, fMinGoF, fMinP, fMinTPCHits, fMinTrkTime, fMinTrkVtx, AlignTPCDrift2::TrkInfo_t::fNChamHit, AlignTPCDrift2::TrkInfo_t::fNTPCHit, AlignTPCDrift2::TrkInfo_t::fPadX, AlignTPCDrift2::TrkInfo_t::fPadZ, AlignTPCDrift2::TrkInfo_t::fPar, fRequireBeamTrk, AlignTPCDrift2::TrkInfo_t::fSinA, AlignTPCDrift2::TrkInfo_t::fSwim, fTrkInfo, AlignTPCDrift2::TrkInfo_t::fU, AlignTPCDrift2::TrkInfo_t::fW, AlignTPCDrift2::TrkInfo_t::fZ0, EDMDataBucket::Get(), GMIPPGeo::Instance(), JobCModule::kFailed, GTPCGeo::kUsecPerBucket, GTPCGeo::PadColToX(), GTPCGeo::PadRowToZ(), EDMEventHandle::Reco(), GTPCGeo::TOffset(), and GMIPPGeo::TPC().

00053 {
00054   vector<const RBKVertex*> vtx(0);
00055 
00056   evt.Cal();
00057   try { evt.Reco().Get(fFolder.c_str(), vtx); }
00058   catch (...) { return kFailed; }
00059   if (vtx.empty()) { return kFailed;}
00060 
00061   TrkInfo_t ti;
00062   ti.fSwim = 0;
00063   GTPCGeo& tpcgeo = GMIPPGeo::Instance().TPC();
00064   
00065   int nGood = 0;
00066   for (unsigned int i = 0; i < vtx.size(); ++i) {
00067     const RBKVertex* v = vtx[i];
00068     if (fRequireBeamTrk && !v->InTrack()) continue;
00069     int nTrk = v->NTracks() + (v->InTrack() ? 1 : 0);
00070     if (nTrk < fMinTrkVtx) continue;
00071 
00072     for (int j = 0; j < v->NTracks(); ++j) {
00073       const RBKTrack* t = v->OutTrack(j);
00074 
00075       if (t->NSpacePoints() < fMinTPCHits ||
00076       t->NSpacePoints() > fMaxTPCHits ||
00077       t->NViewHits() < fMinChamHits ||
00078       fabs(t->QoP()) > fMinP ||
00079       t->Time() < fMinTrkTime ||
00080       t->Time() > fMaxTrkTime ||
00081       t->GoF() < fMinGoF ||
00082       fabs(t->QoP()) < 5e-3) {
00083     if (fDebug) {
00084       cout << "Skipping track: nSpacePts=" << t->NSpacePoints()
00085            << " nViewHits=" << t->NViewHits()
00086            << " QoP=" << t->QoP()
00087            << " T=" << t->Time()
00088            << " GoF=" << t->GoF()
00089            << endl;
00090     }
00091     continue;
00092       }
00093       
00094       ti.fNTPCHit = t->NSpacePoints(); 
00095       if (ti.fNTPCHit > 128) ti.fNTPCHit = 128;
00096       for (int j = 0; j < ti.fNTPCHit; ++j) {
00097     const TPCRHit* hit = 
00098       & dynamic_cast<const TPCRHit&>(t->SpacePoint(j));
00099     assert (hit);
00100     ti.fPadX[j] = tpcgeo.PadColToX(hit->PadCol());
00101     ti.fPadZ[j] = tpcgeo.PadRowToZ(hit->PadRow());
00102     ti.fDriftT[j] = hit->TDC() * GTPCGeo::kUsecPerBucket + 
00103       tpcgeo.TOffset();
00104       }
00105       
00106       ti.fNChamHit = t->NViewHits();
00107       if (ti.fNChamHit > 36) ti.fNChamHit = 36;
00108       for (int j = 0; j < ti.fNChamHit; ++j) {
00109     const RBViewHit& hit = t->ViewHit(j);
00110     ti.fChamZ[j] = hit.Z();
00111     ti.fU[j] = hit.U();
00112     ti.fW[j] = 1 / (hit.USig() * hit.USig());
00113     ti.fCosA[j] = hit.CosTheta();
00114     ti.fSinA[j] = hit.SinTheta();
00115       }
00116       ti.fZ0 = t->VtxZ();
00117       ti.fPar[0] = t->VtxX();
00118       ti.fPar[1] = t->VtxY();
00119       ti.fPar[2] = t->VtxDXDZ();
00120       ti.fPar[3] = t->VtxDYDZ();
00121       ti.fPar[4] = t->QoP();
00122 
00123       double chi2 = ComputeTrkChi2(ti);
00124       if (chi2 > 1e6)  continue;
00125 
00126       if (fDebug) {
00127     printf("Using track: N=%d/%d q/p=%.1e vtx=(%.1f,%.1f;%.1e,%.1e) chi2=%.1f\n",
00128            t->NSpacePoints(), t->NViewHits(), t->QoP(),
00129            t->VtxX(), t->VtxY(), t->VtxDXDZ(), t->VtxDYDZ(),
00130            t->Chi2());
00131     printf("My fit: q/p=%.1e vtx=(%.1f,%.1f;%.1e,%.1e) chi2=%.1f\n",
00132            ti.fPar[4], ti.fPar[0], ti.fPar[1], ti.fPar[2], ti.fPar[3],
00133            chi2);
00134       }
00135       
00136 
00137       fTrkInfo.push_back(ti);
00138       ++nGood;
00139     }
00140   }
00141 
00142   if (fDebug) cout << "  picked " << nGood << " tracks" << endl;
00143 
00144   return (nGood > 0);
00145 }

double AlignTPCDrift2::ComputeTrkChi2 ( TrkInfo_t &  trk  )  [private]

Definition at line 246 of file AlignTPCDrift2.cxx.

References SwimDrift::CorrectDistortion(), AlignTPCDrift2::TrkInfo_t::fChamZ, AlignTPCDrift2::TrkInfo_t::fCosA, AlignTPCDrift2::TrkInfo_t::fDriftT, TrkRUtil::FitCurvedTrack(), AlignTPCDrift2::TrkInfo_t::fNChamHit, AlignTPCDrift2::TrkInfo_t::fNTPCHit, AlignTPCDrift2::TrkInfo_t::fPadX, AlignTPCDrift2::TrkInfo_t::fPadZ, AlignTPCDrift2::TrkInfo_t::fPar, AlignTPCDrift2::TrkInfo_t::fSinA, AlignTPCDrift2::TrkInfo_t::fSwim, fTPCSigmaX, fTPCSigmaY, AlignTPCDrift2::TrkInfo_t::fU, AlignTPCDrift2::TrkInfo_t::fW, AlignTPCDrift2::TrkInfo_t::fZ0, SwimMIPP::GetXY(), GMIPPGeo::Instance(), SwimDrift::Instance(), p, par, GTPCGeo::PPLNToCAVE(), SwimMIPP::Swim(), and GMIPPGeo::TPC().

Referenced by Ana(), and MinFCN().

00247 {
00248   double u[512];
00249   double w[512];
00250   double z[512];
00251   double cosA[512];
00252   double sinA[512];
00253   double lamX[512];
00254   double lamY[512];
00255 
00256   if (!ti.fSwim) {
00257     ti.fSwim = new SwimMIPP;
00258     ti.fSwim->Swim(ti.fZ0, ti.fPar);
00259   }
00260   double p = 1 / ti.fPar[4];
00261 
00262   double padPos[3];
00263   double cavePos[3];
00264   double hitPos[3];
00265   SwimDrift& drift = SwimDrift::Instance();
00266   GTPCGeo& tpcgeo = GMIPPGeo::Instance().TPC();
00267 
00268   for (int i = 0; i < ti.fNTPCHit; ++i) {
00269     // Convert hit position from pad plane to cave
00270     padPos[0] = ti.fPadX[i];
00271     padPos[1] = 0;
00272     padPos[2] = ti.fPadZ[i];
00273     tpcgeo.PPLNToCAVE(padPos, cavePos);    
00274 
00275     // Drift hit to figure out its actual position
00276     drift.CorrectDistortion(cavePos, hitPos, ti.fDriftT[i]);
00277 
00278     z[i * 2 + 1] = z[i * 2] = hitPos[2];
00279     u[i * 2] = hitPos[0];
00280     u[i * 2 + 1] = hitPos[1];
00281     w[i * 2]  = fTPCSigmaX; // sigma's have been converted to weights
00282     w[i * 2 + 1]  = fTPCSigmaY;
00283 
00284     // Compute template coefficients
00285     if (!ti.fSwim->GetXY(hitPos[2], &hitPos[0], &hitPos[1])) return 1e6;
00286     lamX[i * 2 + 1] = lamX[i * 2] = 
00287       p * (hitPos[0] - ti.fPar[0] - ti.fPar[2] * (hitPos[2] - ti.fZ0));
00288     lamY[i * 2 + 1] = lamY[i * 2] = 
00289       p * (hitPos[1] - ti.fPar[1] - ti.fPar[3] * (hitPos[2] - ti.fZ0));
00290 
00291     cosA[i * 2] = 1;
00292     sinA[i * 2] = 0;
00293     cosA[i * 2 + 1] = 0;
00294     sinA[i * 2 + 1] = -1;
00295   }
00296 
00297 
00298   int i = ti.fNTPCHit * 2;
00299   for (int j = 0; j < ti.fNChamHit; ++i, ++j) {
00300     u[i]    = ti.fU[j];
00301     w[i]    = ti.fW[j];
00302     cosA[i] = ti.fCosA[j];
00303     sinA[i] = ti.fSinA[j];
00304     z[i]    = ti.fChamZ[j];
00305     if (!ti.fSwim->GetXY(z[i], &hitPos[0], &hitPos[1])) return 1e6;
00306     lamX[i] = 
00307       p * (hitPos[0] - ti.fPar[0] - ti.fPar[2] * (z[i] - ti.fZ0));
00308     lamY[i] = 
00309       p * (hitPos[1] - ti.fPar[1] - ti.fPar[3] * (z[i] - ti.fZ0));
00310   }
00311 
00312 //   for (int j = 0; j < i; ++j) {
00313 //     printf("z=%.1f u=%.1f cos=%.1f sin=%.1f lamx=%.1f lamy=%.1f\n", 
00314 //         z[j], u[j], cosA[j], sinA[j], lamX[j], lamY[j]);
00315 //   }
00316 
00317   double par[5];
00318   double chi2 =
00319     TrkRUtil::FitCurvedTrack(i, u, w, z, cosA, sinA, ti.fZ0, lamX, lamY, 
00320                  par, 0);
00321 
00322   // Recompute template coefficients
00323   SwimMIPP swim;
00324   swim.Swim(ti.fZ0, par);
00325   p = 1 / par[4];
00326 
00327   for (int j = 0; j < i; ++j) {
00328     if (!ti.fSwim->GetXY(z[i], &hitPos[0], &hitPos[1])) return 1e6;
00329     lamX[i] = 
00330       p * (hitPos[0] - ti.fPar[0] - ti.fPar[2] * (z[i] - ti.fZ0));
00331     lamY[i] = 
00332       p * (hitPos[1] - ti.fPar[1] - ti.fPar[3] * (z[i] - ti.fZ0));    
00333   }
00334 
00335 
00336   chi2 =
00337     TrkRUtil::FitCurvedTrack(i, u, w, z, cosA, sinA, ti.fZ0, lamX, lamY, 
00338                  par, 0);
00339 
00340   return chi2;
00341 }

void AlignTPCDrift2::EndRun ( int  run,
int  subrun 
) [virtual]

Reimplemented from JobCModule.

Definition at line 159 of file AlignTPCDrift2.cxx.

References fFixParam, fTrkInfo, gModule, MinFCN(), minu, and name.

00160 {
00161   if (fTrkInfo.size() < 1) {
00162     cout << "Can't do minimization with " << fTrkInfo.size() << " tracks"
00163      << endl;
00164     return;
00165   }
00166   cout << "AlignTPCDrift2: doing minimization with " << fTrkInfo.size()
00167        << " tracks" << endl;
00168   gModule = this;
00169 
00170   TMinuit* minu = new TMinuit(7);
00171   minu->SetMaxIterations(1000000);
00172   minu->SetFCN(AlignTPCDrift2::MinFCN);
00173   char name[32];
00174   for (int i = 0; i < 3; ++i) {
00175     sprintf(name, "V%c Scale", 'x' + i);
00176     minu->DefineParameter(i, name, 1, 0.01, 0.1, 3);
00177   }
00178   for (int i = 0; i < 3; ++i) {
00179     sprintf(name, "V%c LinScale", 'x' + i);
00180     minu->DefineParameter(i + 3, name, 0, 0.01, -1, 1);
00181   }
00182   minu->DefineParameter(6, "Bx Bz rotation", 0.007, 1e-3, -2e-2, 2e-2);
00183 
00184   for (unsigned int i = 0; i < fFixParam.size(); ++i) 
00185     minu->FixParameter(fFixParam[i]);
00186 
00187   minu->Migrad();
00188 
00189   delete minu;
00190 }

void AlignTPCDrift2::MinFCN ( int &  ,
double *  ,
double &  f,
double *  par,
int   
) [static]

Definition at line 194 of file AlignTPCDrift2.cxx.

References ComputeTrkChi2(), fTrkInfo, gModule, SwimDriftMap::Instance(), SwimDriftMap::SetBxBzRotation(), and SwimDriftMap::SetScale().

Referenced by EndRun().

00196 {
00197   SwimDriftMap::Instance().SetScale(par[0], par[1], par[2],
00198                     par[3], par[4], par[5]);
00199   SwimDriftMap::Instance().SetBxBzRotation(par[6]);
00200 
00201   f = 0;
00202   for (unsigned int i = 0; i < gModule->fTrkInfo.size(); ++i)
00203     f += gModule->ComputeTrkChi2(gModule->fTrkInfo[i]);
00204 
00205   cout << "par=(" << par[0] << ", " << par[1] << ", " << par[2]
00206        << "; " << par[3] << ", " << par[4] << ", " << par[5]
00207        << ", " << par[6] << ") f=" << f << endl;
00208 
00209 }

void AlignTPCDrift2::NewRun ( int  run,
int  subrun 
) [virtual]

Reimplemented from JobCModule.

Definition at line 149 of file AlignTPCDrift2.cxx.

References config_bfield(), SwimDriftMap::Init(), SwimDrift::Instance(), SwimDriftMap::Instance(), GMIPPGeo::Instance(), and SwimDrift::SetUseDriftMap().

00150 {
00151   GMIPPGeo::Instance(run);
00152   config_bfield(run, subrun);
00153   SwimDriftMap::Instance().Init(run);
00154   SwimDrift::Instance().SetUseDriftMap(1);
00155 }

void AlignTPCDrift2::Update ( const CfgConfig &  c  )  [virtual]

Implements CfgObserver.

Definition at line 213 of file AlignTPCDrift2.cxx.

References fDebug, fFixParam, fFolder, JobCModule::fIsInit, fMaxTPCHits, fMaxTrkTime, fMinChamHits, fMinGoF, fMinP, fMinTPCHits, fMinTrkTime, fMinTrkVtx, fRequireBeamTrk, fTPCSigmaX, and fTPCSigmaY.

00214 {
00215   c("Debug").Get(fDebug);
00216   c("Folder").Get(fFolder);
00217   fFixParam.clear();
00218   try { c("FixParam").Get(fFixParam); }
00219   catch (...) { 
00220     int val;
00221     c("FixParam").Get(val);
00222     fFixParam.push_back(val);
00223   }
00224   c("MinTrkVtx").Get(fMinTrkVtx);
00225   c("RequireBeamTrk").Get(fRequireBeamTrk);
00226   c("MinTPCHits").Get(fMinTPCHits);
00227   c("MaxTPCHits").Get(fMaxTPCHits);
00228   c("MinChamHits").Get(fMinChamHits);
00229   c("MinP").Get(fMinP);
00230   fMinP = (fMinP > 0) ? 1 / fMinP : 1e6; // Convert to 1/p
00231   c("MinGoF").Get(fMinGoF);
00232   c("TPCSigmaX").Get(fTPCSigmaX);
00233   fTPCSigmaX = 1 / (fTPCSigmaX * fTPCSigmaX);
00234   c("TPCSigmaY").Get(fTPCSigmaY);
00235   fTPCSigmaX = 1 / (fTPCSigmaY * fTPCSigmaY);
00236 
00237   c("MinTrkTime").Get(fMinTrkTime);
00238   c("MaxTrkTime").Get(fMaxTrkTime);
00239 
00240   // Change the flag to signify that our configuration has been read
00241   fIsInit = true;
00242 }

Member Data Documentation

int AlignTPCDrift2::fDebug [private]

Definition at line 61 of file AlignTPCDrift2.h.

Referenced by Ana(), and Update().

std::vector<int> AlignTPCDrift2::fFixParam [private]

Definition at line 62 of file AlignTPCDrift2.h.

Referenced by EndRun(), and Update().

std::string AlignTPCDrift2::fFolder [private]

Definition at line 63 of file AlignTPCDrift2.h.

Referenced by Ana(), and Update().

int AlignTPCDrift2::fMaxTPCHits [private]

Definition at line 67 of file AlignTPCDrift2.h.

Referenced by Ana(), and Update().

float AlignTPCDrift2::fMaxTrkTime [private]

Definition at line 74 of file AlignTPCDrift2.h.

Referenced by Ana(), and Update().

int AlignTPCDrift2::fMinChamHits [private]

Definition at line 68 of file AlignTPCDrift2.h.

Referenced by Ana(), and Update().

float AlignTPCDrift2::fMinGoF [private]

Definition at line 70 of file AlignTPCDrift2.h.

Referenced by Ana(), and Update().

float AlignTPCDrift2::fMinP [private]

Definition at line 69 of file AlignTPCDrift2.h.

Referenced by Ana(), and Update().

int AlignTPCDrift2::fMinTPCHits [private]

Definition at line 66 of file AlignTPCDrift2.h.

Referenced by Ana(), and Update().

float AlignTPCDrift2::fMinTrkTime [private]

Definition at line 73 of file AlignTPCDrift2.h.

Referenced by Ana(), and Update().

int AlignTPCDrift2::fMinTrkVtx [private]

Definition at line 64 of file AlignTPCDrift2.h.

Referenced by Ana(), and Update().

int AlignTPCDrift2::fRequireBeamTrk [private]

Definition at line 65 of file AlignTPCDrift2.h.

Referenced by Ana(), and Update().

float AlignTPCDrift2::fTPCSigmaX [private]

Definition at line 71 of file AlignTPCDrift2.h.

Referenced by ComputeTrkChi2(), and Update().

float AlignTPCDrift2::fTPCSigmaY [private]

Definition at line 72 of file AlignTPCDrift2.h.

Referenced by ComputeTrkChi2(), and Update().

std::vector<TrkInfo_t> AlignTPCDrift2::fTrkInfo [private]

Definition at line 59 of file AlignTPCDrift2.h.

Referenced by Ana(), EndRun(), and MinFCN().

AlignTPCDrift2 * AlignTPCDrift2::gModule = 0 [static, private]

Definition at line 57 of file AlignTPCDrift2.h.

Referenced by EndRun(), and MinFCN().

The documentation for this class was generated from the following files:
Generated on Mon Nov 23 08:03:56 2009 for MIPP(E907) by  doxygen 1.4.7