AlignChamZ Class Reference

#include <AlignChamZ.h>

Inheritance diagram for AlignChamZ:

JobCModule CfgObserver List of all members.

Public Types

typedef AlignChamZ::track__ TrkInfo_t

Public Member Functions

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

Private Member Functions

double Fit (const TrkInfo_t &ti, double resid[kNCham][kNPlane], double dxdz[kNCham])
void MakeTree ()
double ComputeCorrF ()
double ComputeChi2F ()

Private Attributes

std::vector< TrkInfo_t > fTrackList
double fLambdaX [kNCham][kNPlane]
double fLambdaY [kNCham][kNPlane]
double fZ0
int fDebug
std::string fFolder
int fLoadAlignment
int fMinNClust
std::vector< int > fFixParam
int fUseCorr

Classes

struct  track__

Detailed Description

Author:
Andre Lebedev This module can be used to determine alignment in z for DC1-PWC6 using field-on data. In principle, this should also work for field-off data. Reconstructed TrkCand's are picked up, and all single-wire clusters are used, except for those belonging to BC's. At the end of the run, TMinuit is used to minimize either the sum of all chi squared (faster and I think more reliable) or the sum of the squares of correlation coefficients between residuals and dx/dz.

Definition at line 27 of file AlignChamZ.h.

Member Typedef Documentation

typedef struct AlignChamZ::track__ AlignChamZ::TrkInfo_t

Constructor & Destructor Documentation

AlignChamZ::AlignChamZ ( const char *  version  ) 

Definition at line 35 of file AlignChamZ.cxx.

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

00035                                           :
00036   JobCModule("AlignChamZ")
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 }

AlignChamZ::~AlignChamZ (  ) 

Definition at line 46 of file AlignChamZ.cxx.

00047 {
00048 }

Member Function Documentation

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

Reimplemented from JobCModule.

Definition at line 52 of file AlignChamZ.cxx.

References cham, fFolder, fMinNClust, fTrackList, EDMDataBucket::Get(), kBC3, JobCModule::kFailed, JobCModule::kPassed, p, and EDMEventHandle::Reco().

00053 {
00054   vector<const TrkCand*> vec(0);
00055 
00056   try { evt.Reco().Get(fFolder.c_str(), vec); }
00057   catch (...) { return kFailed; }
00058 
00059   TrkInfo_t ti;
00060 
00061   for (unsigned int i = 0; i < vec.size(); ++i) {
00062     const TrkCand* trk = vec[i];
00063 
00064     memset(&ti, 0, sizeof(TrkInfo_t));
00065 
00066     int nClu = 0;
00067     for (int j = 0; j < trk->NClust(); ++j) {
00068       const WireClust* wc = trk->Clust(j);
00069       int cham = wc->Chamber();
00070       int p = wc->Plane() - 1;
00071 
00072       if (cham > kBC3) nClu++;
00073 
00074       ti.fWire[cham][p] = wc->AvgWire();
00075       ti.fWeight[cham][p] = 1.0 / (wc->USig() * wc->USig());
00076     }
00077     if (nClu < fMinNClust) continue;
00078     fTrackList.push_back(ti);
00079   }
00080 
00081   return kPassed;
00082 }

double AlignChamZ::ComputeChi2F (  )  [private]

Definition at line 106 of file AlignChamZ.cxx.

References TrkRUtil::CalcLambda(), cham, f, TrkRUtil::FitCurvedTrack(), fLambdaX, fLambdaY, fTrackList, fZ0, geo, GMIPPGeo::Instance(), kNCham, kNPlane, nClust, p, and par.

Referenced by ComputeF().

00107 {
00108   double u[kNCham * kNPlane];
00109   double w[kNCham * kNPlane];
00110   double z[kNCham * kNPlane];
00111   double cosA[kNCham * kNPlane];
00112   double sinA[kNCham * kNPlane];
00113   double lambdaX[kNCham * kNPlane];
00114   double lambdaY[kNCham * kNPlane];
00115   int chamber[kNCham * kNPlane];
00116 
00117   double f = 0;
00118   for (unsigned int i = 0; i < fTrackList.size(); ++i) {
00119     TrkInfo_t& ti = fTrackList[i];
00120 
00121     int nClust = 0;
00122     for (int cham = 0; cham < kNCham; ++cham) {
00123       GChamGeo& geo = GMIPPGeo::Instance().Cham(cham);
00124       for (int p = 0; p < kNPlane; ++p) {
00125     if (ti.fWeight[cham][p] < 1e-3) continue;
00126     geo.WireToU(p + 1, ti.fWire[cham][p], u[nClust]);
00127     w[nClust] = ti.fWeight[cham][p];
00128     z[nClust] = geo.Z(p + 1);
00129     cosA[nClust] = geo.GAngleCos(p + 1);
00130     sinA[nClust] = geo.GAngleSin(p + 1);
00131     lambdaX[nClust] = fLambdaX[cham][p];
00132     lambdaY[nClust] = fLambdaY[cham][p];
00133     chamber[nClust] = cham * kNPlane + p;
00134     nClust++;
00135       }
00136     }
00137 
00138     double par[5];
00139     FitCurvedTrack(nClust, u, w, z, cosA, sinA, fZ0, lambdaX, lambdaY,
00140            par, 0);
00141     
00142     // Reinitialize lambda's based on the first-try fit
00143     double lx[kNCham][kNPlane], ly[kNCham][kNPlane];
00144     CalcLambda(fZ0, lx, ly, par);
00145     for (int i = 0; i < nClust; ++i) {
00146       int cham =chamber[i] / kNPlane;
00147       int pl = chamber[i] % kNPlane;
00148       lambdaX[i] = lx[cham][pl];
00149       lambdaY[i] = ly[cham][pl];
00150     }
00151     double chi2 = FitCurvedTrack(nClust, u, w, z, cosA, sinA, fZ0, 
00152                  lambdaX, lambdaY, par, 0);
00153     if (chi2 < 1e5) f += chi2;
00154     else {
00155       cout << "Failed fit?: chi2=" << chi2 << " (" 
00156        <<  par[0] << ", " << par[1]
00157        << "; " << par[2] << ", " << par[3] << "; " << par[4]
00158        << endl;
00159     }
00160   }
00161   return f;
00162 }

double AlignChamZ::ComputeCorrF (  )  [private]

Definition at line 166 of file AlignChamZ.cxx.

References dxdz, f, Fit(), fTrackList, kDC1, kNCham, kNPlane, p, x, and y.

Referenced by ComputeF().

00167 {
00168   double resid[kNCham][kNPlane];
00169   int n[kNCham][kNPlane];
00170   double x[kNCham][kNPlane];
00171   double y[kNCham][kNPlane];
00172   double xx[kNCham][kNPlane];
00173   double yy[kNCham][kNPlane];
00174   double xy[kNCham][kNPlane];
00175   memset(x, 0, sizeof(x));
00176   memset(n, 0, sizeof(n));
00177   memset(y, 0, sizeof(y));
00178   memset(xx, 0, sizeof(xx));
00179   memset(yy, 0, sizeof(yy));
00180   memset(xy, 0, sizeof(xy));
00181 
00182   double dxdz[kNCham];
00183 
00184   for (unsigned int i = 0; i < fTrackList.size(); ++i) {
00185     TrkInfo_t& ti = fTrackList[i];
00186     
00187     Fit(ti, resid, dxdz);
00188 
00189     for (int c = kDC1; c < kNCham; ++c) {
00190       for (int p = 0; p < kNPlane; ++p) {
00191     if (ti.fWeight[c][p] < 1) continue;
00192     n[c][p]++;
00193     x[c][p] += dxdz[c];
00194     xx[c][p] += dxdz[c] * dxdz[c];
00195     y[c][p] += resid[c][p];
00196     yy[c][p] += resid[c][p] * resid[c][p];
00197     xy[c][p] += dxdz[c] * resid[c][p];
00198       }
00199     }
00200   }
00201 
00202   double f = 0;
00203   for (int c = kDC1; c < kNCham; ++c) {
00204     for (int p = 0; p < kNPlane; ++p) {
00205       if (n[c][p] < 10) continue;
00206       
00207       // Add square of correlation coefficient for this plane
00208       double corr = (n[c][p] * xy[c][p] - x[c][p] * y[c][p]);
00209       corr *= corr;
00210       corr /= (n[c][p] * xx[c][p] - x[c][p] * x[c][p]) *
00211     (n[c][p] * yy[c][p] - y[c][p] * y[c][p]);
00212       assert(fabs(corr) < 1);
00213 
00214       //cout << "Correlation coeff for " << c << p << ": " << corr << endl;
00215       
00216       f += corr;
00217     }
00218   }
00219   return f;
00220 }

double AlignChamZ::ComputeF (  ) 

Definition at line 96 of file AlignChamZ.cxx.

References TrkRUtil::CalcLambda(), ComputeChi2F(), ComputeCorrF(), fLambdaX, fLambdaY, fUseCorr, and fZ0.

Referenced by gsFCN().

00097 {
00098   CalcLambda(fZ0, fLambdaX, fLambdaY);
00099   
00100   if (fUseCorr) return ComputeCorrF();
00101   else return ComputeChi2F();
00102 }

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

Reimplemented from JobCModule.

Definition at line 329 of file AlignChamZ.cxx.

References TrkRUtil::CalcLambda(), BFMagnet::CenterPos(), fFixParam, fTrackList, fZ0, geo, gModule, gsFCN(), GMIPPGeo::Instance(), BFMagnet::JGG(), kDC1, MakeTree(), minu, name, GChamGeo::Name(), par, and BFMagnet::Rosie().

00330 {
00331   if (fTrackList.size() < 100) {
00332     cout << "Can't do minimization with less than 100 tracks!" << endl;
00333     return;
00334   }
00335 
00336   cout << "Doing minimization with " << fTrackList.size() << " tracks"
00337        << endl;
00338 
00339   double lx[9][4], ly[9][4];
00340   CalcLambda(fZ0, lx, ly);
00341   for (int i = kDC1; i < 9; ++i) {
00342     cout << "Lambda X " << GChamGeo::Name(i) << ":";
00343     for (int j = 0; j < 4; ++j) cout << lx[i][j] << " ";
00344     cout << "\nLambda Y " << GChamGeo::Name(i) << ":";
00345     for (int j = 0; j < 4; ++j) cout << ly[i][j] << " ";
00346     cout << endl;
00347   }
00348 
00349   gModule = this;
00350   
00351   TMinuit minu(12);
00352   minu.SetFCN(gsFCN);
00353 
00354   char name[256];
00355   int ierr;
00356   double saveZ[6];
00357   for (int i = 0; i < 6; ++i) {
00358     GChamGeo& geo = GMIPPGeo::Instance().Cham(i + kDC1);
00359     sprintf(name, "%s_z", geo.Name());
00360     minu.mnparm(i, name, geo.ZGas(), 0.1, 
00361         geo.ZGas() - 10, geo.ZGas() + 10, ierr);
00362     saveZ[i] = geo.ZGas();
00363   }
00364   BFMagnet& jgg = BFMagnet::JGG();
00365   minu.mnparm(6, "JGGX", jgg.CenterPos()[0], 1, 0, 0, ierr);
00366   minu.mnparm(7, "JGGY", jgg.CenterPos()[1], 1, 0, 0, ierr);
00367   minu.mnparm(8, "JGGZ", jgg.CenterPos()[2], 1, 0, 0, ierr);
00368 
00369   BFMagnet& rosie = BFMagnet::Rosie();
00370   minu.mnparm(9, "RosieX", rosie.CenterPos()[0], 1, 0, 0, ierr);
00371   minu.mnparm(10, "RosieY", rosie.CenterPos()[1], 1, 0, 0, ierr);
00372   minu.mnparm(11, "RosieZ", rosie.CenterPos()[2], 1, 0, 0, ierr);
00373 
00374   for (unsigned int i = 0; i < fFixParam.size(); ++i) 
00375     minu.FixParameter(fFixParam[i]);
00376 
00377   minu.Migrad();
00378 
00379   MakeTree();
00380 
00381   for (int i = 0; i < 6; ++i) {
00382     double par, perr;
00383     minu.GetParameter(i, par, perr);
00384     cout << GChamGeo::Name(i + kDC1) << " z shift: " 
00385      << par - saveZ[i] << " err: " << perr << endl;
00386   }
00387 }

double AlignChamZ::Fit ( const TrkInfo_t &  ti,
double  resid[kNCham][kNPlane],
double  dxdz[kNCham] 
) [private]

Definition at line 224 of file AlignChamZ.cxx.

References TrkRUtil::CalcLambda(), GMIPPGeo::Cham(), cham, dxdz, TrkRUtil::FitCurvedTrack(), fLambdaX, fLambdaY, AlignChamZ::track__::fWeight, AlignChamZ::track__::fWire, fZ0, geo, SwimMIPP::GetPatZ(), GMIPPGeo::Instance(), kDC1, mom, nClust, p, par, SwimMIPP::Swim(), x, and y.

Referenced by ComputeCorrF(), and MakeTree().

00226 {
00227   int nClust = 0;
00228   double u[kNCham * kNPlane];
00229   double w[kNCham * kNPlane];
00230   double z[kNCham * kNPlane];
00231   double cosA[kNCham * kNPlane];
00232   double sinA[kNCham * kNPlane];
00233   double lambdaX[kNCham * kNPlane];
00234   double lambdaY[kNCham * kNPlane];
00235   short chamber[kNCham * kNPlane];
00236   
00237   for (int cham = 0; cham < kNCham; ++cham) {
00238     GChamGeo& geo = GMIPPGeo::Instance().Cham(cham);
00239     for (int p = 0; p < kNPlane; ++p) {
00240       resid[cham][p] = 0;
00241       if (ti.fWeight[cham][p] < 1e-3) continue;
00242       geo.WireToU(p + 1, ti.fWire[cham][p], u[nClust]);
00243       w[nClust] = ti.fWeight[cham][p];
00244       z[nClust] = geo.Z(p + 1);
00245       cosA[nClust] = geo.GAngleCos(p + 1);
00246       sinA[nClust] = geo.GAngleSin(p + 1);
00247       lambdaX[nClust] = fLambdaX[cham][p];
00248       lambdaY[nClust] = fLambdaY[cham][p];
00249       chamber[nClust] = cham * kNPlane + p;
00250       nClust++;
00251     }
00252   }
00253 
00254   double par[5];
00255   FitCurvedTrack(nClust, u, w, z, cosA, sinA, fZ0, lambdaX, lambdaY,
00256                  par, 0);
00257 
00258   // Reinitialize lambda's based on the first-try fit
00259   double lx[kNCham][kNPlane], ly[kNCham][kNPlane];
00260   CalcLambda(fZ0, lx, ly, par);
00261   for (int i = 0; i < nClust; ++i) {
00262     int cham =chamber[i] / kNPlane;
00263     int pl = chamber[i] % kNPlane;
00264     lambdaX[i] = lx[cham][pl];
00265     lambdaY[i] = ly[cham][pl];
00266   }
00267   FitCurvedTrack(nClust, u, w, z, cosA, sinA, fZ0, lambdaX, lambdaY,
00268                  par, 0);
00269   double mom = par[4];
00270   SwimMIPP swim;
00271   swim.Swim(fZ0, par);
00272   for (int c = kDC1; c < kNCham; ++c) {
00273     double p[3];
00274     swim.GetPatZ(GMIPPGeo::Instance().Cham(c).ZGas(), p);
00275     dxdz[c] = p[0] / p[2];
00276   }
00277 
00278   // Loop through every cluster, exclude it from the fit and compute
00279   // residual
00280   for (int i = 0; i < nClust; ++i) {
00281     double saveW = w[i];
00282     w[i] = 0;
00283     FitCurvedTrack(nClust, u, w, z, cosA, sinA, fZ0, lambdaX, lambdaY,
00284            par, 0);
00285     w[i] = saveW;
00286 
00287     int c = chamber[i] / kNPlane;
00288     int p = chamber[i] % kNPlane;
00289     GChamGeo& geo = GMIPPGeo::Instance().Cham(c);
00290 
00291     double dz = geo.Z(p + 1) - fZ0;
00292     double x = par[0] + dz * par[2] + par[4] * lambdaX[i];
00293     double y = par[1] + dz * par[3] + par[4] * lambdaY[i];
00294     double ufit;
00295     geo.XYToU(p + 1, x, y, ufit);
00296     resid[c][p] = ufit - u[i];
00297   }
00298 
00299   return mom;
00300 }

void AlignChamZ::MakeTree (  )  [private]

Definition at line 391 of file AlignChamZ.cxx.

References dxdz, Fit(), fTrackList, geo, GMIPPGeo::Instance(), kNCham, kNPlane, mom, p, and tree.

Referenced by EndRun().

00392 {
00393   double mom;
00394   double resid[kNCham][kNPlane];
00395   double  u[kNCham][kNPlane];
00396   double dxdz[kNCham];
00397 
00398   TTree* tree = new TTree("altree", "Post-alignment residuals");
00399   tree->Branch("mom", &mom, "mom/D");
00400   tree->Branch("resid", &resid, "resid[9][4]/D");
00401   tree->Branch("u", &u, "u[9][4]/D");
00402   tree->Branch("dxdz", &dxdz, "dxdz[9]/D");
00403 
00404   for (unsigned int i = 0; i < fTrackList.size(); ++i) {
00405     TrkInfo_t& ti = fTrackList[i];
00406 
00407     mom = Fit(ti, resid, dxdz);
00408     for (int c = 0; c < kNCham; ++c) {
00409       GChamGeo& geo = GMIPPGeo::Instance().Cham(c);
00410       for (int p = 0; p < kNPlane; ++p) {
00411     if (ti.fWeight[c][p] < 1e-3) u[c][p] = 1e6;
00412     else geo.WireToU(p + 1, ti.fWire[c][p], u[c][p]);
00413       }
00414     }
00415     tree->Fill();
00416   }
00417 
00418   tree->Write();
00419   delete tree;
00420 }

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

Reimplemented from JobCModule.

Definition at line 86 of file AlignChamZ.cxx.

References GMIPPGeo::Cham(), config_bfield(), fLoadAlignment, fZ0, TrkChamGeo::Init(), TrkChamGeo::Instance(), GMIPPGeo::Instance(), kDC1, and GChamGeo::ZGas().

00087 {
00088   GMIPPGeo::Instance(run);
00089   config_bfield(run, subrun);
00090   if (fLoadAlignment) TrkChamGeo::Instance().Init(run);
00091   fZ0 = GMIPPGeo::Instance().Cham(kDC1).ZGas() - 20;
00092 }

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

Implements CfgObserver.

Definition at line 424 of file AlignChamZ.cxx.

References fDebug, fFixParam, fFolder, JobCModule::fIsInit, fLoadAlignment, fMinNClust, and fUseCorr.

00425 {
00426   c("Debug").Get(fDebug);
00427   c("Folder").Get(fFolder);
00428   c("LoadAlignment").Get(fLoadAlignment);
00429   c("MinNClust").Get(fMinNClust);
00430   fFixParam.clear();
00431   try { c("FixParam").Get(fFixParam); }
00432   catch (...) {
00433     int i;
00434     c("FixParam").Get(i);
00435     fFixParam.push_back(i);
00436   }
00437   c("UseCorr").Get(fUseCorr);
00438 
00439   // Change the flag to signify that our configuration has been read
00440   fIsInit = true;
00441 }

Member Data Documentation

int AlignChamZ::fDebug [private]

Definition at line 59 of file AlignChamZ.h.

Referenced by Update().

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

Definition at line 63 of file AlignChamZ.h.

Referenced by EndRun(), and Update().

std::string AlignChamZ::fFolder [private]

Definition at line 60 of file AlignChamZ.h.

Referenced by Ana(), and Update().

double AlignChamZ::fLambdaX[kNCham][kNPlane] [private]

Definition at line 55 of file AlignChamZ.h.

Referenced by ComputeChi2F(), ComputeF(), and Fit().

double AlignChamZ::fLambdaY[kNCham][kNPlane] [private]

Definition at line 56 of file AlignChamZ.h.

Referenced by ComputeChi2F(), ComputeF(), and Fit().

int AlignChamZ::fLoadAlignment [private]

Definition at line 61 of file AlignChamZ.h.

Referenced by NewRun(), and Update().

int AlignChamZ::fMinNClust [private]

Definition at line 62 of file AlignChamZ.h.

Referenced by Ana(), and Update().

std::vector<TrkInfo_t> AlignChamZ::fTrackList [private]

Definition at line 54 of file AlignChamZ.h.

Referenced by Ana(), ComputeChi2F(), ComputeCorrF(), EndRun(), and MakeTree().

int AlignChamZ::fUseCorr [private]

Definition at line 64 of file AlignChamZ.h.

Referenced by ComputeF(), and Update().

double AlignChamZ::fZ0 [private]

Definition at line 57 of file AlignChamZ.h.

Referenced by ComputeChi2F(), ComputeF(), EndRun(), Fit(), and NewRun().

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