AlgFitTrackAtNu Class Reference

#include <AlgFitTrackAtNu.h>

Inheritance diagram for AlgFitTrackAtNu:
AlgBase

List of all members.

Public Member Functions

 AlgFitTrackAtNu ()
 ~AlgFitTrackAtNu ()
void RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx)
void Trace (const char *c) const

Private Member Functions

Double_t GetTrackQP (TVector3 P, TVector3 dPdZ, TVector3 B)
Int_t FitTrack (Int_t N, TMatrixD *x, TMatrixD *y, TMatrixD *err, Int_t M, TMatrixD *p, TMatrixD *dp, TMatrixD *dpdx)

Private Attributes

TObjArray fTrkStrpList [500]
TObjArray fSliStrpList [500]

Detailed Description

Definition at line 11 of file AlgFitTrackAtNu.h.


Constructor & Destructor Documentation

AlgFitTrackAtNu::AlgFitTrackAtNu (  ) 

Definition at line 26 of file AlgFitTrackAtNu.cxx.

References Msg::kVerbose, and MSG.

00027 {
00028   MSG("AlgFitTrackAtNu",Msg::kVerbose) << " *** AlgFitTrackAtNu::AlgFitTrackAtNu() *** " << endl;
00029 }

AlgFitTrackAtNu::~AlgFitTrackAtNu (  ) 

Definition at line 31 of file AlgFitTrackAtNu.cxx.

References Msg::kVerbose, and MSG.

00032 {
00033   MSG("AlgFitTrackAtNu",Msg::kVerbose) << " *** AlgFitTrackAtNu::~AlgFitTrackAtNu() *** " << endl;
00034 }


Member Function Documentation

Int_t AlgFitTrackAtNu::FitTrack ( Int_t  N,
TMatrixD *  x,
TMatrixD *  y,
TMatrixD *  err,
Int_t  M,
TMatrixD *  p,
TMatrixD *  dp,
TMatrixD *  dpdx 
) [private]

Definition at line 779 of file AlgFitTrackAtNu.cxx.

References MinosMaterial::A(), C, Msg::kDebug, Munits::m, MSG, and n.

Referenced by RunAlg().

00780 {
00781   // This method performs M-th order polynomial fit
00782   // N measurements, M parameters
00783 
00784   Double_t u,v;
00785   Int_t i,n,m,q; 
00786 
00787   if(M>N){
00788     MSG("AlgFitTrackAtNu", Msg::kDebug) << " *** not enough measurements for fit *** " << endl;
00789     return -1;
00790   }
00791 
00792   // expanion matrix L(m,n) = F_m( X(n) )   
00793   TMatrixD L(M,N); // M down, N across  
00794   for(n=0;n<N;n++){
00795     for(m=0;m<M;m++){
00796       u=(*x)(n,0); v=0.0;
00797       switch(m){
00798         case 0 : v=1.0; break; 
00799         case 1 : v=1.0*u; break; 
00800         case 2 : v=1.0*u*u; break; 
00801         case 3 : v=1.0*u*u*u; break; 
00802         case 4 : v=1.0*u*u*u*u; break; 
00803         case 5 : v=1.0*u*u*u*u*u; break;
00804         default : break; 
00805       }
00806       L(m,n) = v;
00807     }
00808   }
00809 
00810   // measurement matrix
00811   TMatrixD Y(N,1); // N down, 1 across
00812   for(n=0;n<N;n++){
00813     Y(n,0) = (*y)(n,0);
00814   }
00815 
00816   // measurement error matrix 
00817   TMatrixD dY(N,N); // N down, N across
00818   for(n=0;n<N;n++){
00819     for(i=0;i<N;i++){
00820       if(i==n) dY(n,i)=((*err)(n,0))*((*err)(n,0)); else dY(n,i)=0.0;
00821     }
00822   }
00823 
00824   // parameter matrix A(m) = coefficient of F_m
00825   TMatrixD A(M,1); // M down, 1 across
00826 
00827   // parameter error matrix
00828   TMatrixD dA(M,M); // M down, M across
00829 
00830   // normal equations:
00831   // L(M,N) C(N,N) Y(N,0)  =  L(M,N) C(N,N) L_t(N,M) A(M,0)
00832   // C(N,N) = covariance matrix = inv(dY)
00833 
00834   // solutions:
00835   // A(M)  =  T(M,N) Y(N)  =  inv(L C L_t) (L C) Y(N)
00836   // dA(M,M) = T(M,N) dY(N,N) T_t(N,M)
00837 
00838   // Tranpose L
00839   TMatrixD Lt(N,M);
00840   for(n=0;n<N;n++){
00841     for(m=0;m<M;m++){
00842       Lt(n,m)=L(m,n);
00843     }
00844   }
00845 
00846   // Calculate A
00847   TMatrixD C(N,N); C=dY; C.Invert();
00848   TMatrixD LC(M,N); LC.Mult(L,C);
00849   TMatrixD LCL(M,M); LCL.Mult(LC,Lt); LCL.Invert();
00850   TMatrixD T(M,N); T.Mult(LCL,LC);
00851   A.Mult(T,Y);
00852 
00853   // Tranpose T
00854   TMatrixD Tt(N,M);
00855   for(n=0;n<N;n++){
00856     for(m=0;m<M;m++){
00857       Tt(n,m)=T(m,n);
00858     }
00859   }
00860 
00861   // Calculate dA
00862   TMatrixD dYTt(N,M); dYTt.Mult(dY,Tt);
00863   dA.Mult(T,dYTt);
00864 
00865   // return dA/dx
00866   for(n=0;n<N;n++){
00867     for(m=0;m<M;m++){
00868       (*dpdx)(n,m)=Tt(n,m);
00869     }
00870   }
00871 
00872   // return A and dA
00873   for(m=0;m<M;m++){
00874     (*p)(m,0)=A(m,0);
00875     for(q=0;q<M;q++){
00876       (*dp)(m,q)=dA(m,q);
00877     }
00878   }
00879 
00880   return 1;
00881 }

Double_t AlgFitTrackAtNu::GetTrackQP ( TVector3  P,
TVector3  dPdZ,
TVector3  B 
) [private]

Definition at line 754 of file AlgFitTrackAtNu.cxx.

Referenced by RunAlg().

00755 {
00756   // This method calculates Q/P from fit
00757   Double_t qp;
00758   Double_t dUdZ,d2UdZ2,dVdZ,d2VdZ2;
00759   Double_t Bu,Bv,Bz;
00760   TVector3 p,b,pxb,dpds,ndpds,pdnds;
00761 
00762   qp=0.0;
00763   dUdZ=P.X(); dVdZ=P.Y(); 
00764   d2UdZ2=dPdZ.X(); d2VdZ2=dPdZ.Y();
00765   Bu=B.X(); Bv=B.Y(); Bz=B.Z();
00766   p.SetXYZ(dUdZ,dVdZ,1.0);
00767   b.SetXYZ(Bu,Bv,Bz);
00768   if(p.Mag2()>0.0 && b.Mag2()>0.0){
00769     pxb=(1.0/p.Mag())*p.Cross(b);
00770     ndpds.SetXYZ(d2UdZ2,d2VdZ2,0.0);
00771     pdnds.SetXYZ(dUdZ,dVdZ,1.0); 
00772     dpds=(1.0/p.Mag2())*(ndpds-pdnds*((d2UdZ2*dUdZ+d2VdZ2*dVdZ)/(p.Mag2())));
00773     qp=(dpds.Dot(pxb))/(0.3*pxb.Mag2()); 
00774   }
00775 
00776   return qp;
00777 }

void AlgFitTrackAtNu::RunAlg ( AlgConfig ac,
CandHandle ch,
CandContext cx 
) [virtual]

Implements AlgBase.

Definition at line 36 of file AlgFitTrackAtNu.cxx.

References CandHandle::AddDaughterLink(), FitTrack(), fSliStrpList, fTrkStrpList, BField::GetBField(), CandContext::GetCandRecord(), CandRecoHandle::GetCandSlice(), CandStripHandle::GetCharge(), CandFitTrackHandle::GetChi2(), CandFitTrackAtNuHandle::GetChi2Lin(), CandContext::GetDataIn(), CandHandle::GetDaughterIterator(), VldContext::GetDetector(), CandRecoHandle::GetDirCosU(), CandRecoHandle::GetDirCosV(), CandRecoHandle::GetDirCosZ(), CandFitTrackHandle::GetEMCharge(), CandRecoHandle::GetEndDirCosU(), CandRecoHandle::GetEndDirCosV(), CandRecoHandle::GetEndDirCosZ(), CandRecoHandle::GetEndPlane(), CandRecoHandle::GetEndU(), CandRecoHandle::GetEndV(), CandRecoHandle::GetEndZ(), Registry::GetInt(), CandTrackHandle::GetMomentum(), CandFitTrackHandle::GetMomentumCurve(), CandFitTrackAtNuHandle::GetMomentumCurveErr(), CandFitTrackHandle::GetMomentumRange(), CandFitTrackHandle::GetPass(), CandStripHandle::GetPlaneView(), CandFitTrackAtNuHandle::GetQPcorr(), CandFitTrackAtNuHandle::GetQPerr(), CandFitTrackAtNuHandle::GetQPmean(), CandFitTrackAtNuHandle::GetQPplns(), CandFitTrackAtNuHandle::GetQPwidth(), CandTrackHandle::GetRange(), UgliGeomHandle::GetSteelPlnHandle(), VHS::GetStrip(), CandStripHandle::GetStrip(), CandRecoHandle::GetTimeSlope(), CandStripHandle::GetTPos(), GetTrackQP(), RecMinos::GetVldContext(), CandRecoHandle::GetVtxPlane(), CandFitTrackHandle::GetVtxQPError(), CandRecoHandle::GetVtxU(), CandRecoHandle::GetVtxV(), CandRecoHandle::GetVtxZ(), UgliSteelPlnHandle::GetZ0(), CandStripHandle::GetZPos(), PlaneView::kA, PlaneView::kB, Detector::kCalib, Msg::kDebug, Detector::kFar, Detector::kNear, BfldInterpMethod::kPlanar, PlaneView::kU, PlaneView::kV, PlaneView::kX, PlaneView::kY, MSG, CandFitTrackHandle::SetChi2(), CandFitTrackAtNuHandle::SetChi2Lin(), CandRecoHandle::SetDirCosU(), CandRecoHandle::SetDirCosV(), CandRecoHandle::SetDirCosZ(), CandFitTrackHandle::SetEMCharge(), CandRecoHandle::SetEndDirCosU(), CandRecoHandle::SetEndDirCosV(), CandRecoHandle::SetEndDirCosZ(), CandRecoHandle::SetEndPlane(), CandRecoHandle::SetEndT(), CandTrackHandle::SetEndTrace(), CandRecoHandle::SetEndU(), CandRecoHandle::SetEndV(), CandRecoHandle::SetEndZ(), CandFitTrackHandle::SetFinderTrack(), BField::SetInterpMethod(), CandTrackHandle::SetMomentum(), CandFitTrackHandle::SetMomentumCurve(), CandFitTrackAtNuHandle::SetMomentumCurveErr(), CandFitTrackHandle::SetMomentumRange(), CandFitTrackHandle::SetNIterate(), CandFitTrackHandle::SetPass(), CandFitTrackAtNuHandle::SetQPcorr(), CandFitTrackAtNuHandle::SetQPerr(), CandFitTrackAtNuHandle::SetQPmean(), CandFitTrackAtNuHandle::SetQPplns(), CandFitTrackAtNuHandle::SetQPwidth(), CandRecoHandle::SetTimeOffset(), CandRecoHandle::SetTimeSlope(), CandRecoHandle::SetVtxPlane(), CandFitTrackHandle::SetVtxQPError(), CandRecoHandle::SetVtxT(), CandTrackHandle::SetVtxTrace(), CandRecoHandle::SetVtxU(), CandRecoHandle::SetVtxV(), CandRecoHandle::SetVtxZ(), and MinosMaterial::Z().

00037 {
00038 
00039   MSG("AlgFitTrackAtNu", Msg::kDebug) << " *** AlgFitTrackAtNu::RunAlg(...) *** " << endl;
00040 
00041   CandFitTrackAtNuHandle& fittrack = dynamic_cast<CandFitTrackAtNuHandle&>(ch);
00042 
00043   // Unpack AlgConfig
00044   Int_t fBFieldMap=201;
00045   fBFieldMap = ac.GetInt("BFieldMap");
00046   MSG("AlgFitTrackAtNu", Msg::kDebug) << "   BFieldMap = " << fBFieldMap << endl;
00047 
00048   // Unpack CandContext
00049   const CandTrackHandle* track = dynamic_cast<const CandTrackHandle*>(cx.GetDataIn());
00050   MSG("AlgFitTrackAtNu", Msg::kDebug) << " *** FOUND THE TRACK *** " << endl;
00051 
00052   const CandSliceHandle* slice = dynamic_cast<const CandSliceHandle*>(track->GetCandSlice());
00053   MSG("AlgFitTrackAtNu", Msg::kDebug) << " *** FOUND THE SLICE *** " << endl;
00054   
00055   CandRecord* candrec = (CandRecord*)(cx.GetCandRecord());
00056   VldContext *vldc = (VldContext*)(candrec->GetVldContext());
00057 
00058   // Set up the BField
00059   MSG("AlgFitTrackAtNu", Msg::kDebug) << " *** SETTING UP MAGNETIC FIELD *** " << endl;
00060   BField* bf = new BField(*vldc,-1,fBFieldMap);
00061   bf->SetInterpMethod(BfldInterpMethod::kPlanar);
00062 
00063   // Set up the Geometry
00064   MSG("AlgFitTrackAtNu", Msg::kDebug) << " *** SETTING UP GEOMETRY *** " << endl;
00065   UgliGeomHandle ugh(*vldc);
00066 
00067   Int_t modtype,modnum,modpln,modstr;
00068   Int_t fEndPln,fEndStr;
00069   switch(vldc->GetDetector()){
00070     case Detector::kFar:
00071       modtype=1; modnum=2; modpln=248; modstr=192; break;
00072     case Detector::kNear:
00073       modtype=2; modnum=1; modpln=282; modstr=96; break;
00074     case Detector::kCalib:
00075       modtype=3; modnum=1; modpln=60; modstr=24; break;
00076     default:
00077       modtype=-1; modnum=0; modpln=0; break;
00078   }
00079   fEndPln = modnum*(modpln+1); fEndStr = modstr;
00080 
00081   MSG("AlgFitTrackAtNu", Msg::kDebug) << " modtype=" << modtype 
00082                                       << " modnum=" << modnum 
00083                                       << " modpln=" << modpln 
00084                                       << " modstr=" << modstr << endl;
00085 
00086   MSG("AlgFitTrackAtNu", Msg::kDebug) << " *** UNPACKING STRIPS *** " << endl;
00087 
00088   Int_t nplns=0;
00089   Int_t i,j,k;
00090   Double_t minz,maxz;
00091   Int_t pln,mod,vuw,shw;
00092   Int_t bstr=-1,estr=-1,bpln=-1,epln=-1;
00093   Int_t bplnu=-1,eplnu=-1,bplnv=-1,eplnv=-1;
00094   Double_t dir=0.0,range=0.0,ptrk=0.0;
00095 
00096   Int_t* moduplns = new Int_t[modnum];
00097   Int_t* modvplns = new Int_t[modnum];
00098   for(mod=0;mod<modnum;mod++){
00099     moduplns[mod]=0; modvplns[mod]=0;
00100   }
00101 
00102   TIter trackitr(track->GetDaughterIterator());
00103   while(CandStripHandle* strp = (CandStripHandle*)(trackitr())){
00104     vuw=-1; 
00105     pln=strp->GetPlane(); mod=(Int_t)((pln)/(modpln+1.0));
00106     if( strp->GetPlaneView()==PlaneView::kU
00107      || strp->GetPlaneView()==PlaneView::kX 
00108      || strp->GetPlaneView()==PlaneView::kA ){ vuw=0; }
00109     if( strp->GetPlaneView()==PlaneView::kV
00110      || strp->GetPlaneView()==PlaneView::kY
00111      || strp->GetPlaneView()==PlaneView::kB ){ vuw=1; }
00112     if(pln>0&&pln<500) fTrkStrpList[pln].Add(strp);  
00113     if(bpln<0||pln<bpln){ bpln=pln; minz=strp->GetZPos(); } 
00114     if(epln<0||pln>epln){ epln=pln; maxz=strp->GetZPos(); }
00115     if(vuw==0){ moduplns[mod]++; }  if(vuw==1){ modvplns[mod]++; }
00116     if(vuw==0){ if(bplnu<0||pln<bplnu) bplnu=pln; if(eplnu<0||pln>eplnu) eplnu=pln; }
00117     if(vuw==1){ if(bplnv<0||pln<bplnv) bplnv=pln; if(eplnv<0||pln>eplnv) eplnv=pln; }
00118     fittrack.AddDaughterLink(*strp);
00119   }
00120   if(bpln>0 && epln>0){
00121     nplns=1+epln-bpln;
00122   }
00123 
00124   TIter sliceitr(slice->GetDaughterIterator());
00125   while(CandStripHandle* strp = (CandStripHandle*)(sliceitr())){
00126     pln=strp->GetPlane(); if(pln>0&&pln<500) fSliStrpList[pln].Add(strp);
00127   }
00128 
00129   ptrk=track->GetMomentum();
00130   if(track->GetEndPlane()-track->GetVtxPlane()>=0) dir=1.0; 
00131   if(track->GetEndPlane()-track->GetVtxPlane()<0) dir=-1.0;
00132 
00133   for(i=bpln;i<epln+1;i++){
00134     if(track->GetRange(i)>range) range=track->GetRange(i);
00135   }
00136 
00137   Double_t* U = new Double_t[nplns];
00138   Double_t* V = new Double_t[nplns];
00139   Double_t* Z = new Double_t[nplns];
00140   Double_t* dU = new Double_t[nplns];
00141   Double_t* dV = new Double_t[nplns];
00142   Double_t* dE = new Double_t[nplns];
00143   Int_t* plnshw = new Int_t[nplns];
00144   Int_t* plnvuw = new Int_t[nplns];
00145   Int_t* plnpln = new Int_t[nplns];
00146   Int_t* plnpos = new Int_t[nplns];
00147 
00148   Int_t Nplns=0,Uplns=0,Vplns=0;
00149   Int_t Npos=0,Upos=0,Vpos=0;
00150   Double_t Sn,Sq,Sqt,Sqtt;
00151   Double_t u,v,z,s,du,dv,dz,dt;
00152   for(i=0;i<nplns;i++){
00153     plnvuw[i]=-1;
00154 
00155     pln=bpln+i; 
00156     mod=(Int_t)((pln)/(modpln+1.0));
00157     if(moduplns[mod]>2 && modvplns[mod]>2){ 
00158 
00159       vuw=-1; shw=0; bstr=-1; estr=-1;
00160       u=0.0; v=0.0; z=0.0; s=0.0;
00161       du=0.0; dv=0.0; dt=0.0;
00162 
00163       if(1+fTrkStrpList[pln].GetLast()>0){
00164         Sq=0.0; Sqt=0.0;
00165         for(j=0;j<1+fTrkStrpList[pln].GetLast();j++){
00166           CandStripHandle* strp = (CandStripHandle*)(fTrkStrpList[pln].At(j));
00167           z=strp->GetZPos();
00168           // range defined from END, s defined from BEGINNING
00169           if(dir>=0) s=range-track->GetRange(pln); else s=track->GetRange(pln);
00170           // range defined from BEGINNING, s defined from BEGINNING
00171           // if(dir>=0) s=track->GetRange(pln); else s=range-track->GetRange(pln);
00172           if( strp->GetPlaneView()==PlaneView::kU
00173            || strp->GetPlaneView()==PlaneView::kX 
00174            || strp->GetPlaneView()==PlaneView::kA ){ vuw=0; }
00175           if( strp->GetPlaneView()==PlaneView::kV
00176            || strp->GetPlaneView()==PlaneView::kY
00177            || strp->GetPlaneView()==PlaneView::kB ){ vuw=1; }
00178           if(bstr<0||strp->GetStrip()<bstr) bstr=strp->GetStrip();
00179           if(estr<0||strp->GetStrip()>estr) estr=strp->GetStrip();
00180           Sqt+=strp->GetTPos()*strp->GetCharge(); Sq+=strp->GetCharge();
00181         }
00182         if(Sq>0.0){
00183           if(vuw==0){ u=Sqt/Sq; Uplns++; } if(vuw==1){ v=Sqt/Sq; Vplns++; }
00184         }
00185       }
00186 
00187       if(1+fSliStrpList[pln].GetLast()>0){
00188         Sn=0.0; Sq=0.0; Sqt=0.0; Sqtt=0.0;
00189         for(j=0;j<1+fSliStrpList[pln].GetLast();j++){
00190           CandStripHandle* strp = (CandStripHandle*)(fSliStrpList[pln].At(j));
00191           if( bstr>-1 && estr>-1
00192            && strp->GetStrip()>bstr-3 && strp->GetStrip()<estr+3 ){
00193             Sqtt+=strp->GetTPos()*strp->GetTPos()*strp->GetCharge(); 
00194             Sqt+=strp->GetTPos()*strp->GetCharge(); 
00195             Sq+=strp->GetCharge();
00196           }
00197           if( strp->GetCharge()>3.0 ) Sn+=1.0;
00198         }
00199         if( Sq>0.0 ){
00200           if( (Sqtt/Sq)-(Sqt/Sq)*(Sqt/Sq)>0.0 ){
00201             dt=sqrt( (Sqtt/Sq)-(Sqt/Sq)*(Sqt/Sq) );  
00202           }   
00203           if(vuw==0){ du=dt+0.012; } if(vuw==1){ dv=dt+0.012; }
00204         }
00205         if( Sn>2.0 ){ shw=1; }
00206       }
00207 
00208       if( vuw>-1 && Nplns<nplns ){
00209         plnvuw[Nplns]=vuw; plnshw[Nplns]=shw; plnpln[Nplns]=pln; plnpos[Nplns]=-1;
00210         U[Nplns]=u; V[Nplns]=v; Z[Nplns]=z;
00211         dU[Nplns]=du; dV[Nplns]=dv; dE[Nplns]=ptrk*s/range;
00212         Nplns++;
00213       }
00214     }
00215 
00216   }
00217     
00218   Int_t im=-1,ip=-1;
00219   for(i=0;i<Nplns;i++){
00220     pln=plnpln[i];
00221     if( pln>bplnu && pln>bplnv && pln<eplnu && pln<eplnv ){
00222       plnpos[i]=Npos; Npos++;
00223       if(plnvuw[i]==0) Upos++; if(plnvuw[i]==1) Vpos++;
00224       if(im<0||i<im) im=i; if(ip<0||i>ip) ip=i;
00225     }
00226   }
00227 
00228   MSG("AlgFitTrackAtNu", Msg::kDebug) << " bpln=" << bpln << " epln=" << epln << endl;
00229     
00230   Double_t bvtxu,bvtxv,bvtxz;
00231   Double_t bdiru,bdirv,bdirz;
00232   Double_t evtxu,evtxv,evtxz;
00233   Double_t ediru,edirv,edirz;
00234 
00235   Double_t pcurv=0.0,pcurverr=0.0;
00236   Double_t qpmean=0.0,qperr=0.0,qpcorr=0.0,qpwidth=0.0;
00237   Double_t chisqpol=9999.9,chisqlin=9999.9;
00238   Double_t rmspol=9999.9,rmslin=9999.9;
00239   Int_t qpplns=0,emcharge=0,pass=0;
00240 
00241   if( Uplns>2 && Vplns>2 && Upos>0 && Vpos>0 ){
00242     MSG("AlgFitTrackAtNu", Msg::kDebug) << " *** FITTING TRACK *** " << endl;
00243 
00244    /*******************************
00245     * P O L Y N O M I A L   F I T *
00246     *******************************/
00247 
00248     MSG("AlgFitTrackAtNu", Msg::kDebug) << " *** POLYNOMIAL FIT *** " << endl;
00249 
00250     Int_t ipos,jpos;
00251     Int_t begi,endi;
00252     Double_t b,de,qp,tmpqp;
00253     Double_t x0,y0,z0;
00254     Double_t u0,u1,u2,v0,v1,v2;
00255     Double_t cu,cv,pu,pv;
00256     Int_t modbpln,modepln;
00257     Int_t utmp,vtmp,utot,vtot;
00258     Int_t Ntmp,Nlin;
00259     Int_t Upars,Vpars;
00260     Double_t dparam,dqpdau,dqpdav;
00261     Double_t dau,dav,dqpu,dqpv,dqp;
00262     Double_t bx,by,bz;
00263     Double_t Srms,Stot;
00264     TVector3 RXYZ,BXYZ;
00265     TVector3 B,P,dPdZ;
00266     Int_t flagU,flagV;
00267 
00268     TMatrixD* fB = new TMatrixD(Npos,1);
00269     TMatrixD* fQP = new TMatrixD(Npos,1);
00270     TMatrixD* fdE = new TMatrixD(Npos,1);
00271     TMatrixD* fdT = new TMatrixD(Npos,1);
00272     TMatrixD* fdQP = new TMatrixD(Npos,1);
00273     TMatrixD* fdQPdX = new TMatrixD(Nplns,Npos); 
00274     TMatrixD* fdQPdXt = new TMatrixD(Npos,Nplns);  
00275     TMatrixD* fdXdX = new TMatrixD(Nplns,Nplns); 
00276 
00277     for(i=0;i<Nplns;i++){
00278       for(j=0;j<Npos;j++){
00279         (*fdQPdX)(i,j)=0.0; (*fdQPdXt)(j,i)=0.0;
00280       }
00281       if(plnvuw[i]==0) (*fdXdX)(i,i)=dU[i]*dU[i]; if(plnvuw[i]==1) (*fdXdX)(i,i)=dV[i]*dV[i];
00282     }
00283 
00284 
00285     // CALCULATE Q/P
00286     for(i=im;i<=ip;i++){
00287 
00288       de=dE[i];
00289       qp=0.0; dqp=9999.9; dt=0.0; b=0.0;
00290       pln=plnpln[i]; ipos=plnpos[i];
00291       utot=0; vtot=0; begi=i; endi=i; 
00292       if(plnvuw[i]==0) utot=1; if(plnvuw[i]==1) vtot=1;
00293       mod=(Int_t)((pln)/(modpln+1.0));
00294       modbpln=mod*(modpln+1); modepln=(mod+1)*(modpln+1);
00295       PlexPlaneId planeid(vldc->GetDetector(),pln,1);
00296       UgliSteelPlnHandle planehandle=ugh.GetSteelPlnHandle(planeid);
00297       z0=planehandle.GetZ0();
00298 
00299       utmp=0; vtmp=0; k=0;
00300       while( i-k>0 && plnpln[i-k]>modbpln && plnpln[i-k]<modepln && (utmp<3||vtmp<3) ){
00301         k++; if(plnvuw[i-k]==0) utmp++; if(plnvuw[i-k]==1) vtmp++;
00302       }
00303       utot+=utmp; vtot+=vtmp; begi=i-k;
00304 
00305       utmp=0; vtmp=0; k=0;
00306       while( i+k<Nplns-1 && plnpln[i+k]>modbpln && plnpln[i+k]<modepln && (utmp<3||vtmp<3) ){
00307         k++; if(plnvuw[i+k]==0) utmp++; if(plnvuw[i+k]==1) vtmp++;
00308       }
00309       utot+=utmp; vtot+=vtmp; endi=i+k;
00310 
00311       // fit U view
00312       TMatrixD* ZUtmp = new TMatrixD(utot,1);
00313       TMatrixD* Utmp = new TMatrixD(utot,1);
00314       TMatrixD* dUtmp = new TMatrixD(utot,1);
00315       Int_t* JUtmp = new Int_t[utot];
00316       Ntmp=0; 
00317 
00318       for(j=begi;j<=endi;j++){
00319         if(Ntmp<utot && plnvuw[j]==0){
00320           u=U[j]; du=dU[j]; dz=Z[j]-z0;
00321           (*Utmp)(Ntmp,0)=u;
00322           (*dUtmp)(Ntmp,0)=du;
00323           (*ZUtmp)(Ntmp,0)=dz;
00324           JUtmp[Ntmp]=j; Ntmp++;
00325         }
00326       }
00327 
00328       Upars=3;
00329       TMatrixD* PUtmp = new TMatrixD(3,1);
00330       TMatrixD* dPUtmp = new TMatrixD(3,3);
00331       TMatrixD* dPdUtmp = new TMatrixD(Ntmp,3);
00332       flagU=this->FitTrack(Ntmp,ZUtmp,Utmp,dUtmp,Upars,PUtmp,dPUtmp,dPdUtmp);
00333       u0=(*PUtmp)(0,0); u1=(*PUtmp)(1,0); u2=(*PUtmp)(2,0); 
00334 
00335 
00336       // fit V view
00337       TMatrixD* ZVtmp = new TMatrixD(vtot,1);
00338       TMatrixD* Vtmp = new TMatrixD(vtot,1);
00339       TMatrixD* dVtmp = new TMatrixD(vtot,1);
00340       Int_t* JVtmp = new Int_t[vtot];
00341       Ntmp=0; 
00342 
00343       for(j=begi;j<=endi;j++){
00344         if(Ntmp<vtot && plnvuw[j]==1){
00345           v=V[j]; dv=dV[j]; dz=Z[j]-z0;
00346           (*Vtmp)(Ntmp,0)=v;
00347           (*dVtmp)(Ntmp,0)=dv;
00348           (*ZVtmp)(Ntmp,0)=dz;
00349           JVtmp[Ntmp]=j; Ntmp++;
00350         }
00351       }
00352 
00353       Vpars=3;
00354       TMatrixD* PVtmp = new TMatrixD(3,1);
00355       TMatrixD* dPVtmp = new TMatrixD(3,3);
00356       TMatrixD* dPdVtmp = new TMatrixD(Ntmp,3);
00357       flagV=this->FitTrack(Ntmp,ZVtmp,Vtmp,dVtmp,Vpars,PVtmp,dPVtmp,dPdVtmp);
00358       v0=(*PVtmp)(0,0); v1=(*PVtmp)(1,0); v2=(*PVtmp)(2,0); 
00359 
00360 
00361       // Calculate Q/P
00362       if( flagU>-1 && flagV>-1 ){
00363 
00364         if(plnvuw[i]==0) dt=u0-U[i]; if(plnvuw[i]==1) dt=v0-V[i];
00365         x0=0.7071*(u0-v0); y0=0.7071*(u0+v0);
00366         RXYZ.SetXYZ(x0,y0,z0); 
00367         BXYZ=bf->GetBField(RXYZ);
00368         bx=BXYZ.X(); by=BXYZ.Y(); bz=BXYZ.Z(); b=BXYZ.Mag();
00369         B.SetXYZ(0.7071*(by+bx),0.7071*(by-bx),bz);
00370         P.SetXYZ(u1,v1,1.0); dPdZ.SetXYZ(2.0*u2,2.0*v2,0.0);
00371         qp=this->GetTrackQP( P,dPdZ,B );
00372   
00373         dau=sqrt((*dPUtmp)(2,2)); dparam=0.01*dau;
00374         dPdZ.SetXYZ(2.0*(u2+dparam),2.0*v2,0.0);          
00375         tmpqp=this->GetTrackQP( P,dPdZ,B );
00376         dqpdau=(tmpqp-qp)/dparam; dqpu=dqpdau*dau;
00377         for(j=0;j<utot;j++){
00378           jpos=JUtmp[j];
00379           (*fdQPdX)(jpos,ipos)=(dqpdau)*((*dPdUtmp)(j,2));
00380           (*fdQPdXt)(ipos,jpos)=(*fdQPdX)(jpos,ipos);
00381         }
00382       
00383         dav=sqrt((*dPVtmp)(2,2)); dparam=0.01*dav;
00384         dPdZ.SetXYZ(2.0*u2,2.0*(v2+dparam),0.0);          
00385         tmpqp=this->GetTrackQP( P,dPdZ,B );
00386         dqpdav=(tmpqp-qp)/dparam; dqpv=dqpdav*dav;      
00387         for(j=0;j<vtot;j++){
00388           jpos=JVtmp[j];
00389           (*fdQPdX)(jpos,ipos)=(dqpdav)*((*dPdVtmp)(j,2));
00390           (*fdQPdXt)(ipos,jpos)=(*fdQPdX)(jpos,ipos);
00391         }
00392 
00393         dqp=sqrt(dqpu*dqpu+dqpv*dqpv); 
00394       }
00395 
00396       MSG("AlgFitTrackAtNu", Msg::kDebug) << " PLN=" << pln 
00397                                           << " QP=" << qp 
00398                                           << " dQP=" << dqp 
00399                                           << " dT=" << dt << endl;
00400 
00401       delete [] JUtmp;
00402       delete [] JVtmp;
00403 
00404       delete ZUtmp;
00405       delete Utmp;
00406       delete dUtmp;
00407       delete PUtmp;
00408       delete dPUtmp;
00409       delete dPdUtmp;
00410 
00411       delete ZVtmp;
00412       delete Vtmp;
00413       delete dVtmp;
00414       delete PVtmp;
00415       delete dPVtmp;
00416       delete dPdVtmp;
00417 
00418       (*fB)(ipos,0)=b; (*fdT)(ipos,0)=dt; (*fdE)(ipos,0)=de;
00419       (*fQP)(ipos,0)=qp; (*fdQP)(ipos,0)=dqp; 
00420     }
00421 
00422     
00423     // CALCULATE BEST FIT MOMENTUM
00424     MSG("AlgFitTrackAtNu", Msg::kDebug) << " *** calculating best fit momentum *** " << endl;
00425     Int_t myqpplns=0,myplns=0,mypass=0;
00426     Double_t myqpmean=0.0,myqperr=0.0,myqpcorr=0.0,myqpwidth=0.0;
00427     Double_t si,sj,xi,xj,dxi,dxj,dxij;
00428     Double_t Sw2dx2,Swx2,Swx,Sws,Sw;
00429     Double_t sigma,mean;
00430     Int_t itr;
00431 
00432     // best fit Q/P;
00433     TMatrixD* fQPdQP = new TMatrixD(Nplns,Npos);
00434     TMatrixD* fdQPdQP = new TMatrixD(Npos,Npos);
00435 
00436     fQPdQP->Mult(*fdXdX,*fdQPdX);
00437     fdQPdQP->Mult(*fdQPdXt,*fQPdQP);
00438 
00439     mean=0.0; sigma=4.0;
00440     for(itr=0;itr<2;itr++){
00441       myplns=0; mypass=0;
00442       Sw2dx2=0.0; Swx2=0.0; Swx=0.0; Sws=0.0; Sw=0.0; 
00443       for(i=0;i<Npos;i++){
00444         xi=(*fQP)(i,0); dxi=(*fdQP)(i,0); si=(*fdE)(i,0);
00445         if( dxi!=0.0 && xi>mean-2.0*sigma-0.5 && xi<mean+2.0*sigma+0.5 ){
00446           Swx2+=xi*xi/(dxi*dxi); Swx+=xi/(dxi*dxi); Sws+=xi*si/(dxi*dxi); 
00447           Sw+=1.0/(dxi*dxi); myplns++;
00448           for(j=0;j<Npos;j++){
00449             xj=(*fQP)(j,0); dxj=(*fdQP)(j,0); sj=(*fdE)(j,0);
00450             if( dxj!=0.0 && xj>mean-2.0*sigma-0.5 && xj<mean+2.0*sigma+0.5 ){
00451               dxij=(*fdQPdQP)(i,j); Sw2dx2+=dxij/(dxi*dxi*dxj*dxj);
00452             }
00453           }
00454         }
00455       }
00456       if(Sw>0.0){
00457         myqpmean=Swx/Sw; myqpcorr=Sws/Sw; myqpplns=myplns;
00458         if(Sw2dx2>=0.0) myqperr=sqrt(Sw2dx2)/Sw; else myqperr=-sqrt(-Sw2dx2)/Sw;
00459         if((Swx2/Sw)-(Swx/Sw)*(Swx/Sw)>0.0) myqpwidth=sqrt((Swx2/Sw)-(Swx/Sw)*(Swx/Sw)); else myqpwidth=0.0;
00460         mean=myqpmean; sigma=myqpwidth; mypass=1;
00461       }
00462     }
00463 
00464     qpmean=2.3*myqpmean; qperr=2.3*myqperr; qpwidth=2.3*myqpwidth;
00465     if(myqpmean>=0.0) qpcorr=2.3*myqpmean/(1.0+2.3*myqpcorr); else qpcorr=2.3*myqpmean/(1.0-2.3*myqpcorr);
00466     qpplns=myqpplns; pass=mypass;
00467 
00468     MSG("AlgFitTrackAtNu", Msg::kDebug) << " BEST FIT Q/P : " << endl
00469                                         << "   Q/P mean = " << qpmean << endl
00470                                         << "   Q/P err = " << qperr << endl
00471                                         << "   Q/P corr = " << qpcorr << endl 
00472                                         << "   Q/P width = " << qpwidth << endl
00473                                         << "   Q/P plns = " << qpplns << endl;
00474     
00475 
00476     pcurv=0.0; pcurverr=0.0;
00477     if( qpcorr!=0.0 ) pcurv=1.0/qpcorr;
00478 
00479 
00480     // calculate chisqpol
00481     Srms=0.0; Stot=0.0;
00482     for(i=0;i<Npos;i++){  
00483       dt=(*fdT)(i,0);  
00484       if(dt*dt>0.0){
00485         Srms+=dt*dt; Stot+=1.0;
00486       }  
00487     }
00488     if( Stot>0.0 ){
00489       rmspol=sqrt(Srms/Stot);
00490     }
00491 
00492     chisqpol=rmspol;
00493 
00494     MSG("AlgFitTrackAtNu", Msg::kDebug) << "  chisqpol = " << chisqpol << endl;
00495 
00496     delete fB;
00497     delete fQP;
00498     delete fdE;
00499     delete fdT;
00500     delete fdQP;
00501 
00502     delete fdQPdX;
00503     delete fdQPdXt;
00504     delete fdXdX;
00505     delete fQPdQP;
00506     delete fdQPdQP;
00507      
00508 
00509     /***********************
00510      * L I N E A R   F I T *
00511      ***********************/
00512   
00513     MSG("AlgFitTrackAtNu", Msg::kDebug) << " *** LINEAR FIT *** " << endl; 
00514 
00515     // fit U view
00516     TMatrixD* ZUlin = new TMatrixD(Uplns,1);
00517     TMatrixD* Ulin = new TMatrixD(Uplns,1);
00518     TMatrixD* dUlin = new TMatrixD(Uplns,1);
00519     Nlin=0; 
00520 
00521     for(i=0;i<Nplns;i++){
00522       if(Nlin<Uplns && plnvuw[i]==0){
00523         (*Ulin)(Nlin,0)=U[i];
00524         (*dUlin)(Nlin,0)=dU[i];
00525         (*ZUlin)(Nlin,0)=Z[i];
00526         Nlin++;
00527       }
00528     }
00529 
00530     Upars=2;
00531     TMatrixD* PUlin = new TMatrixD(2,1);
00532     TMatrixD* dPUlin = new TMatrixD(2,2);
00533     TMatrixD* dPdUlin = new TMatrixD(Nlin,2);
00534     flagU=this->FitTrack(Nlin,ZUlin,Ulin,dUlin,Upars,PUlin,dPUlin,dPdUlin);
00535     cu=(*PUlin)(0,0); pu=(*PUlin)(1,0); 
00536 
00537     // fit V view
00538     TMatrixD* ZVlin = new TMatrixD(Vplns,1);
00539     TMatrixD* Vlin = new TMatrixD(Vplns,1);
00540     TMatrixD* dVlin = new TMatrixD(Vplns,1);
00541     Nlin=0; 
00542 
00543     for(i=0;i<Nplns;i++){
00544       if(Nlin<Vplns && plnvuw[i]==1){
00545         (*Vlin)(Nlin,0)=V[i];
00546         (*dVlin)(Nlin,0)=dV[i];
00547         (*ZVlin)(Nlin,0)=Z[i];
00548         Nlin++;
00549       }
00550     }
00551 
00552     Vpars=2;
00553     TMatrixD* PVlin = new TMatrixD(2,1);
00554     TMatrixD* dPVlin = new TMatrixD(2,2);
00555     TMatrixD* dPdVlin = new TMatrixD(Nlin,2);
00556     flagV=this->FitTrack(Nlin,ZVlin,Vlin,dVlin,Vpars,PVlin,dPVlin,dPdVlin);
00557     cv=(*PVlin)(0,0); pv=(*PVlin)(1,0); 
00558   
00559     // calculate chisqlin
00560     rmslin = -9999.9;
00561     Srms=0.0; Stot=0.0;
00562     if( flagU>-1 && flagV>-1 ){
00563       for(i=0;i<Nplns;i++){
00564         dt=0.0;
00565         if(plnvuw[i]==0) dt=U[i]-(cu+pu*Z[i]); if(plnvuw[i]==1) dt=V[i]-(cv+pv*Z[i]);
00566         if(dt*dt>0.0){
00567           Srms+=dt*dt; Stot+=1.0;
00568         }
00569       }
00570       if(Stot>0.0){
00571         rmslin=sqrt(Srms/Stot);
00572       }
00573     }
00574 
00575     chisqlin=rmslin;
00576 
00577     MSG("AlgFitTrackAtNu", Msg::kDebug) << "  chisqlin=" << chisqlin << endl;
00578 
00579     delete ZUlin;
00580     delete Ulin;
00581     delete dUlin;
00582     delete PUlin;
00583     delete dPUlin;
00584     delete dPdUlin;
00585 
00586     delete ZVlin;
00587     delete Vlin;
00588     delete dVlin;
00589     delete PVlin;
00590     delete dPVlin;
00591     delete dPdVlin;
00592 
00593   }
00594 
00595 
00596   // set charge
00597   if(qpcorr>0) emcharge=1; if(qpcorr<0) emcharge=-1;
00598 
00599   MSG("AlgFitTrackAtNu", Msg::kDebug) << " *** Q=" << emcharge << " *** " << endl;
00600 
00601   delete [] moduplns;
00602   delete [] modvplns;
00603 
00604   delete [] U;
00605   delete [] V;
00606   delete [] Z;
00607   delete [] dU;
00608   delete [] dV;
00609   delete [] dE;
00610 
00611   delete [] plnshw;
00612   delete [] plnvuw;
00613   delete [] plnpln;
00614   delete [] plnpos;
00615 
00616 
00617   if(track->GetTimeSlope()>=0.0){
00618     bvtxu=track->GetVtxU();
00619     bvtxv=track->GetVtxV();
00620     bvtxz=track->GetVtxZ();
00621     bdiru=track->GetDirCosU();
00622     bdirv=track->GetDirCosV();
00623     bdirz=track->GetDirCosZ();
00624     evtxu=track->GetEndU();
00625     evtxv=track->GetEndV();
00626     evtxz=track->GetEndZ();
00627     ediru=track->GetEndDirCosU();
00628     edirv=track->GetEndDirCosV();
00629     edirz=track->GetEndDirCosZ();
00630   }
00631   else{
00632     bvtxu=track->GetEndU();
00633     bvtxv=track->GetEndV();
00634     bvtxz=track->GetEndZ();
00635     bdiru=-track->GetEndDirCosU();
00636     bdirv=-track->GetEndDirCosV();
00637     bdirz=-track->GetEndDirCosZ();
00638     evtxu=track->GetVtxU();
00639     evtxv=track->GetVtxV();
00640     evtxz=track->GetVtxZ();
00641     ediru=-track->GetDirCosU();
00642     edirv=-track->GetDirCosV();
00643     edirz=-track->GetDirCosZ();
00644   }
00645 
00646   fittrack.SetTimeSlope((dir)/(3.0e8));
00647   fittrack.SetTimeOffset(0.0);
00648   fittrack.SetMomentumRange(ptrk);
00649   fittrack.SetMomentum(0.0);
00650 
00651   if(dir>=0.0){
00652     fittrack.SetVtxU(bvtxu);
00653     fittrack.SetVtxV(bvtxv);
00654     fittrack.SetVtxZ(bvtxz);
00655     fittrack.SetVtxT(0.0);
00656     fittrack.SetVtxPlane(bpln);
00657     fittrack.SetVtxTrace(0.0);
00658     fittrack.SetDirCosU(bdiru);
00659     fittrack.SetDirCosV(bdirv);
00660     fittrack.SetDirCosZ(bdirz);  
00661     fittrack.SetEndU(evtxu);
00662     fittrack.SetEndV(evtxv);
00663     fittrack.SetEndZ(evtxz);
00664     fittrack.SetEndT(0.0);
00665     fittrack.SetEndPlane(epln);
00666     fittrack.SetEndTrace(0.0);
00667     fittrack.SetEndDirCosU(ediru);
00668     fittrack.SetEndDirCosV(edirv);
00669     fittrack.SetEndDirCosZ(edirz);
00670   }
00671   else{
00672     fittrack.SetVtxU(evtxu);
00673     fittrack.SetVtxV(evtxv);
00674     fittrack.SetVtxZ(evtxz);
00675     fittrack.SetVtxT(0.0);
00676     fittrack.SetVtxPlane(epln);
00677     fittrack.SetVtxTrace(0.0);
00678     fittrack.SetDirCosU(-ediru);
00679     fittrack.SetDirCosV(-edirv);
00680     fittrack.SetDirCosZ(-edirz);
00681     fittrack.SetEndU(bvtxu);
00682     fittrack.SetEndV(bvtxv);
00683     fittrack.SetEndZ(bvtxz);
00684     fittrack.SetEndT(0.0);
00685     fittrack.SetEndPlane(bpln);
00686     fittrack.SetEndTrace(0.0);
00687     fittrack.SetEndDirCosU(-bdiru);
00688     fittrack.SetEndDirCosV(-bdirv);
00689     fittrack.SetEndDirCosZ(-bdirz);
00690   }
00691 
00692   // charge
00693   if(dir>=0.0){
00694     fittrack.SetEMCharge(emcharge);  
00695     fittrack.SetQPmean(qpmean);
00696     fittrack.SetQPerr(qperr);
00697     fittrack.SetQPcorr(qpcorr);
00698     fittrack.SetQPwidth(qpwidth);
00699     fittrack.SetQPplns(qpplns);
00700     fittrack.SetMomentumCurve(emcharge*pcurv);
00701     fittrack.SetMomentumCurveErr(pcurverr);
00702     fittrack.SetVtxQPError(qperr);
00703   }
00704   else{
00705     fittrack.SetEMCharge(-emcharge);  
00706     fittrack.SetQPmean(-qpmean);
00707     fittrack.SetQPerr(qperr);
00708     fittrack.SetQPcorr(-qpcorr);
00709     fittrack.SetQPwidth(qpwidth);
00710     fittrack.SetQPplns(qpplns);
00711     fittrack.SetMomentumCurve((-emcharge)*(-pcurv));
00712     fittrack.SetMomentumCurveErr(pcurverr);
00713     fittrack.SetVtxQPError(qperr);
00714   }
00715 
00716   fittrack.SetChi2Lin(chisqlin);
00717   fittrack.SetChi2(chisqpol);
00718   fittrack.SetPass(pass);  
00719   fittrack.SetNIterate(1);
00720   fittrack.SetFinderTrack((CandTrackHandle*)(track));
00721   
00722 
00723   MSG("AlgFitTrackAtNu", Msg::kDebug) 
00724     << " *** CREATING NEW FITTED TRACK *** " << endl
00725     << "  EMcharge = " << fittrack.GetEMCharge() << endl
00726     << "  Chisq = " << fittrack.GetChi2() << endl
00727     << "  ChisqLin = " << fittrack.GetChi2Lin() << endl
00728     << "  MomentumRange = " << fittrack.GetMomentumRange() << endl
00729     << "  MomentumCurve = " << fittrack.GetMomentumCurve() << endl
00730     << "  MomentumCurveErr = " << fittrack.GetMomentumCurveErr() << endl 
00731     << "  VtxQPError = " << fittrack.GetVtxQPError() << endl
00732     << "  QPmean = " << fittrack.GetQPmean() << endl
00733     << "  QPerr = " << fittrack.GetQPerr() << endl
00734     << "  QPcorr = " << fittrack.GetQPcorr() << endl
00735     << "  QPwidth = " << fittrack.GetQPwidth() << endl
00736     << "  QPplns = " << fittrack.GetQPplns() << endl
00737     << "  Pass = " << fittrack.GetPass() << endl;
00738 
00739   // clear banks
00740   for(i=0;i<500;i++){
00741     fTrkStrpList[i].Clear(); fSliStrpList[i].Clear();
00742   }
00743   
00744   delete bf;
00745 
00746   return;
00747 }

void AlgFitTrackAtNu::Trace ( const char *  c  )  const [virtual]

Reimplemented from AlgBase.

Definition at line 749 of file AlgFitTrackAtNu.cxx.

00750 {
00751 
00752 }


Member Data Documentation

TObjArray AlgFitTrackAtNu::fSliStrpList[500] [private]

Definition at line 24 of file AlgFitTrackAtNu.h.

Referenced by RunAlg().

TObjArray AlgFitTrackAtNu::fTrkStrpList[500] [private]

Definition at line 23 of file AlgFitTrackAtNu.h.

Referenced by RunAlg().


The documentation for this class was generated from the following files:

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1