#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 26 of file AlgFitTrackAtNu.cxx. References MSG. 00027 {
00028 MSG("AlgFitTrackAtNu",Msg::kVerbose) << " *** AlgFitTrackAtNu::AlgFitTrackAtNu() *** " << endl;
00029 }
|
|
|
Definition at line 31 of file AlgFitTrackAtNu.cxx. References MSG. 00032 {
00033 MSG("AlgFitTrackAtNu",Msg::kVerbose) << " *** AlgFitTrackAtNu::~AlgFitTrackAtNu() *** " << endl;
00034 }
|
|
||||||||||||||||||||||||||||||||||||
|
Definition at line 779 of file AlgFitTrackAtNu.cxx. References StandardMaterial::A(), C, and MSG. 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 }
|
|
||||||||||||||||
|
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 }
|
|
||||||||||||||||
|
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::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(), CandRecoHandle::GetVtxPlane(), 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(), and CandRecoHandle::SetVtxZ(). 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 }
|
|
|
Reimplemented from AlgBase. Definition at line 749 of file AlgFitTrackAtNu.cxx. 00750 {
00751
00752 }
|
|
|
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