#include <ExtMdcTrack.h>
|
|
00029 :myMsgFlag(true),myInputTrk("Kal"),myParID(2),myHelixPar(NumHelixPar),myMdcErr(NdimMdcErr,0),myPhiTerm(0.){}
|
|
|
00031 {}
|
|
|
|
|
|
|
|
|
|
|
|
00203 {
00204 myHelixPar(1)*=10.;
00205 myHelixPar(3)*=0.001;
00206 myHelixPar(4)*=10.;
00207 myPivot*=10.;
00208 myMdcErr.fast(1,1)*=100.;
00209 myMdcErr.fast(2,1)*=10.;
00210 myMdcErr.fast(3,1)*=0.01;
00211 myMdcErr.fast(3,2)*=0.001;
00212 myMdcErr.fast(3,3)*=0.000001;
00213 myMdcErr.fast(4,1)*=100.;
00214 myMdcErr.fast(4,2)*=10.;
00215 myMdcErr.fast(4,3)*=0.01;
00216 myMdcErr.fast(4,4)*=100.;
00217 myMdcErr.fast(5,1)*=10.;
00218 myMdcErr.fast(5,3)*=0.001;
00219 myMdcErr.fast(5,4)*=10.;
00220 }
|
|
|
|
|
|
00252 {
00253 Ext_Helix aHelix(myPivot,myHelixPar,myMdcErr);
00254 // aHelix.pivot(aHelix.x(myPhiTerm));//Move the pivot to the PhiTerm position.
00255
00256 Hep3Vector aPoint = aHelix.x(myPhiTerm);
00257 double phi0 = aHelix.phi0();
00258 double kappaInv = 1.0/aHelix.kappa();
00259 double cosPhi0 = cos(phi0);
00260 double sinPhi0 = sin(phi0);
00261 // 0 -> myPhiTerm, 2009-5-29, by wangll
00262 // double px = aHelix.momentum(0).x();
00263 // double py = aHelix.momentum(0).y();
00264 // double pz = aHelix.momentum(0).z();
00265 double px = aHelix.momentum(myPhiTerm).x();
00266 double py = aHelix.momentum(myPhiTerm).y();
00267 double pz = aHelix.momentum(myPhiTerm).z();
00268
00269 //Calculate the Jacobian: d(xp)/d(helix).
00270 HepMatrix Helix2XpJcb(NdimExtErr,NdimMdcErr,0);
00271
00272 Helix2XpJcb(1,1) = cosPhi0;
00273 Helix2XpJcb(1,2) = myPivot.y()-aPoint.y();
00274 Helix2XpJcb(1,3) = (myPivot.x()-aPoint.x()+myHelixPar(1)*cosPhi0)*kappaInv;
00275 Helix2XpJcb(2,1) = sinPhi0;
00276 Helix2XpJcb(2,2) = aPoint.x()-myPivot.x();
00277 Helix2XpJcb(2,2) = (myPivot.y()-aPoint.y()+myHelixPar(1)*sinPhi0)*kappaInv;
00278 Helix2XpJcb(3,3) = (myPivot.z()-aPoint.z()+myHelixPar(4))*kappaInv;
00279 Helix2XpJcb(3,4) = 1.0;
00280 Helix2XpJcb(3,5) = (aPoint.z()-myPivot.z()-myHelixPar(4))/myHelixPar(5);
00281 Helix2XpJcb(4,2) = -py;
00282 Helix2XpJcb(4,3) = -px * kappaInv;
00283 Helix2XpJcb(5,2) = px;
00284 Helix2XpJcb(5,3) = -py * kappaInv;
00285 Helix2XpJcb(6,3) = -pz * kappaInv;
00286 Helix2XpJcb(6,5) = fabs(kappaInv);
00287
00288 //Calculate the Ext error matrix(6x6) in (x,p).
00289 // if(myMsgFlag) cout<<"Error matrix at new pivot: "<<endl;
00290 // aHelix.Ea();
00291 // cout<<"haha"<<endl;
00292 HepSymMatrix aErrorMatrix = aHelix.Ea().similarity(Helix2XpJcb);
00293 if(myMsgFlag) cout<<"PhiTerm Error matrix:"<<aErrorMatrix<<endl;
00294
00295 //Uints Conversion.
00296 for(int i=1;i<=6;i++)
00297 {
00298 for(int j=1;j<=i;j++)
00299 {
00300 if(i<=3 && j<=3) aErrorMatrix.fast(i,j)*=100.0;
00301 if(i>=4 && j>=4) aErrorMatrix.fast(i,j)*=1.0e+6;
00302 if(i>=4 && j<=3) aErrorMatrix.fast(i,j)*=10000.0;
00303 }
00304 }
00305
00306 return aErrorMatrix;
00307 }
|
|
|
|
|
|
00223 {
00224 return (*myMdcRecTrkIter);
00225 }
|
|
|
|
|
|
00243 {
00244 Ext_Helix aHelix(myPivot,myHelixPar,myMdcErr);
00245 Hep3Vector aMomentum = aHelix.momentum(myPhiTerm);//In GeV/c units.
00246 if(myMsgFlag) cout<<"PhiTerm Momentum: "<<aMomentum<<endl;
00247 aMomentum*=1000.0;//In MeV/c uints.
00248 return (aMomentum);
00249 }
|
|
|
|
|
|
00052 {
00053 if(myMsgFlag) cout<<"ExtMdcTrack::GetOneGoodTrk() begin"<<endl;
00054 if(myInputTrk=="Mdc")
00055 {
00056 myMdcRecTrkIter++;
00057 if(myMdcRecTrkIter!=myMdcTrackCol->end()) return true;
00058 else
00059 {
00060 if(myMsgFlag) cout<<"No more track!"<<endl;
00061 return false;
00062 }
00063 }
00064 else if(myInputTrk=="Kal") // use KalTrack as input
00065 {
00066 myMdcKalTrkIter++;
00067 if(myMdcKalTrkIter!=myMdcKalTrackCol->end()) return true;
00068 else
00069 {
00070 if(myMsgFlag) cout<<"No more track!"<<endl;
00071 return false;
00072 }
00073 }
00074 }
|
|
|
|
|
|
00357 {
00358 return ((myHelixPar[2]>0.)? 1.0:-1.0);
00359 }
|
|
|
|
|
|
00228 {
00229 // if(myInputTrk=="Mdc") {
00230 Ext_Helix aHelix(myPivot,myHelixPar,myMdcErr);
00231 Hep3Vector aPoint = aHelix.x(myPhiTerm);
00232 if(myMsgFlag) cout<<"PhiTerm Position: "<<aPoint<<endl;
00233 aPoint*=10.0;
00234 return (aPoint);
00235 /* }
00236 else if(myInputTrk=="Kal") {
00237 return (myLPosition);
00238 }
00239 */
00240 }
|
|
|
|
|
|
00352 {
00353 return (1.0/fabs(myHelixPar[2]));
00354 }
|
|
|
00050 {return myTrackID;}//Get RecMdcTrackl ID.
|
|
|
00050 {return myTrackID;}//Get RecMdcTrackl ID.
|
|
|
|
|
|
00310 {
00311 if(myInputTrk=="Mdc") {
00312
00313 /* //old version
00314 Helix aHelix(myPivot,myHelixPar,myMdcErr);
00315
00316 double phi0 = aHelix.phi0();
00317 HepPoint3D w = -(*this).GetParticleCharge() * aHelix.center();
00318 Hep3Vector center( w.x(),w.y(),w.z() );
00319 double centerX = center.x();
00320 double centerY = center.y();
00321 double phiCenter;
00322 if(fabs(centerX) > 1.0e-10)
00323 {
00324 phiCenter = atan2(centerY,centerX);
00325 if(phiCenter < 0.) phiCenter += 2.0*M_PI;
00326 }
00327 else
00328 {
00329 phiCenter = (centerY>0)? M_PI_2:3.0*M_PI_2;
00330 }
00331 double dPhi = fabs(phi0 + myPhiTerm - phiCenter);
00332 if(dPhi >= 2.0*M_PI) dPhi -= 2.0*M_PI;
00333 double tanLambda = aHelix.tanl();
00334 double cosLambdaInv = sqrt(tanLambda*tanLambda + 1.0);
00335 return (10*fabs(aHelix.radius() * dPhi * cosLambdaInv));//10* for cm -> mm
00336 */
00337
00338 Ext_Helix aHelix(myPivot,myHelixPar,myMdcErr);
00339 double dPhi = fabs(myPhiTerm);
00340 double tanLambda = aHelix.tanl();
00341 double cosLambdaInv = sqrt(tanLambda*tanLambda + 1.0);
00342 return (10*fabs(aHelix.radius() * dPhi * cosLambdaInv));//10* for cm -> mm
00343
00344 }
00345 else if(myInputTrk=="Kal")
00346 {
00347 return 10*myTrkLength; //10* for cm -> mm
00348 }
00349 }
|
|
|
00056 {return myTrkTof;};
|
|
|
00056 {return myTrkTof;};
|
|
|
|
|
|
00078 {
00079 if(myMsgFlag) cout<<"ExtMdcTrack::ReadTrk() begin"<<endl;
00080 myParID = pid;
00081 if(myInputTrk=="Mdc")
00082 {
00083 if(myMdcRecTrkIter!=myMdcTrackCol->end())
00084 {
00085 RecMdcTrack *aMdcTrack = *myMdcRecTrkIter;
00086 // if(!aMdcTrack->getStat())//If good track.Now status default 0
00087 {
00088 if(myMsgFlag) cout<<"Read a good track!"<<endl;
00089
00090 //get MdcTrk data
00091 myTrackID = aMdcTrack->trackId();
00092 myHelixPar = aMdcTrack->helix();
00093 myPivot = aMdcTrack->getPivot();
00094 myMdcErr = aMdcTrack->err();
00095 // double aPhiTerm = aMdcTrack->getFiTerm();
00096 myPhiTerm = aMdcTrack->getFiTerm();
00097 myTrkLength = GetTrackLength();
00098 double p = (GetMomentum()).mag();
00099 double beta = p/sqrt(mass[myParID]*mass[myParID]+p*p);
00100 myTrkTof = myTrkLength/(beta*299.792458);
00101
00102 /* //Start caculation of myPhiTerm.(when pivot is (0,0,0).)
00103 double r = DBL_MAX;
00104 double rMdc = 81.0;
00105 double dro = myHelixPar(1);
00106 if(myHelixPar(3)!=0.0)
00107 r = 333.564095/myHelixPar(3);
00108 if( fabs(2.0*r+dro)>rMdc && fabs(dro)<rMdc )
00109 {
00110 double aValue = ((dro+r)*(dro+r)+r*r-rMdc*rMdc)/2.0/fabs(r)/fabs(dro+r);
00111 if(r<0.0) myPhiTerm = acos(aValue);
00112 else myPhiTerm = -1.0*acos(aValue);
00113 }
00114 else
00115 {
00116 if(myMsgFlag) cout<<"Get a track which can't go out of MDC!"<<endl;
00117 myMdcRecTrkIter++;
00118 continue;
00119 }
00120 */
00121 // double aTrackLength=fabs(r*myPhiTerm)*sqrt(myHelixPar(5)*myHelixPar(5)+1);
00122 if(myMsgFlag)
00123 {
00124 cout<<"MDC track:"<<myTrackID<<endl;
00125 cout<<"Helix: "<<myHelixPar;
00126 cout<<"pivot: "<<myPivot<<"cm"<<endl;
00127 cout<<"Error: "<<myMdcErr<<endl;
00128
00129 /* Convert();
00130 cout<<"********after convert:"<<endl;
00131 cout<<"Helix: "<<myHelixPar;
00132 cout<<"pivot: "<<myPivot;
00133 cout<<"Error: "<<myMdcErr<<endl;
00134 */
00135 cout<<"Pt: "<<GetPt()<<"GeV/c"<<endl;
00136 cout<<"PhiTerm: "<<myPhiTerm<<endl;
00137 // cout<<"recPhiTerm: "<<aPhiTerm<<endl;
00138 cout<<"trackLength(mm): "<<GetTrackLength()<<endl;
00139 // cout<<"trackLengthTest(cm): "<<aTrackLength<<endl;
00140 }
00141
00142 return true;
00143 }
00144 // else
00145 // {
00146 // myMdcRecTrkIter++;
00147 // cout<<"Get a bad track!"<<endl;
00148 // }
00149 }
00150 else
00151 {
00152 if(myMsgFlag) cout<<"No more track!"<<endl;
00153 return false;
00154 }
00155 }
00156 else if(myInputTrk=="Kal") // use KalTrack as input
00157 {
00158 if(myMdcKalTrkIter!=myMdcKalTrackCol->end())
00159 {
00160 RecMdcKalTrack *aMdcKalTrack = *myMdcKalTrkIter;
00161 {
00162 if(myMsgFlag) cout<<"Get a good track!"<<endl;
00163
00164 //get MdcTrk data
00165 myTrackID = aMdcKalTrack->getTrackId();
00166 myHelixPar = aMdcKalTrack->getLHelix(myParID);
00167 //myPivot = HepPoint3D(0,0,0);
00168 myPivot = aMdcKalTrack->getLPivot(myParID);
00169 myMdcErr = aMdcKalTrack->getLError(myParID);
00170 myPhiTerm = aMdcKalTrack->getFiTerm(myParID);
00171 myPhiTerm = 0.0;
00172 myTrkLength = aMdcKalTrack->getPathSM(myParID);
00173 myLPosition = aMdcKalTrack->getLPoint(myParID);
00174 myTrkTof = aMdcKalTrack->getTof(myParID);
00175 //cout<<"myInitialTof = "<<myTrkTof<<endl;
00176
00177 if(myMsgFlag)
00178 {
00179 cout<<"Kal track:"<<myTrackID<<endl;
00180 cout<<"myParID:"<<myParID<<endl;
00181 cout<<"Helix: "<<myHelixPar;
00182 cout<<"pivot: "<<myPivot<<"cm"<<endl;
00183 cout<<"Error: "<<myMdcErr<<endl;
00184 cout<<"Pt: "<<GetPt()<<"GeV/c"<<endl;
00185 cout<<"PhiTerm: "<<myPhiTerm<<endl;
00186 // cout<<"trackLength(mm): "<<GetTrackLength()<<endl;
00187 cout<<"trackLength(cm): "<<myTrkLength<<endl;
00188 cout<<"myTrkTof: "<<myTrkTof<<endl;
00189 cout<<"Last point position: "<<myLPosition<<endl;
00190 }
00191 return true;
00192 }
00193 }
00194 else
00195 {
00196 if(myMsgFlag) cout<<"No more track!"<<endl;
00197 return false;
00198 }
00199 }
00200 }
|
|
|
|
|
|
00043 {
00044 myMdcKalTrackCol = aPointer;
00045 myMdcKalTrkIter = myMdcKalTrackCol->begin();
00046 myMdcKalTrkIter--;
00047 myInputTrk="Kal";
00048 return true;
00049 }
|
|
|
|
|
|
00034 {
00035 myMdcTrackCol = aPointer;
00036 myMdcRecTrkIter = myMdcTrackCol->begin();
00037 myMdcRecTrkIter--;
00038 myInputTrk="Mdc";
00039 return true;
00040 }
|
|
|
00044 {myMsgFlag=aFlag;};
|
|
|
00044 {myMsgFlag=aFlag;};
|
|
|
00047 { myParID=pid; return true;};
|
|
|
00047 { myParID=pid; return true;};
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1.3.9.1