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