00001 #include "DTagAlg/utility.h"
00002
00003 HepLorentzVector utility::getp4(RecMdcKalTrack* mdcKalTrack, int pid){
00004
00005 HepVector zhelix;
00006 double mass=0;
00007
00008 if(pid==0){
00009 zhelix=mdcKalTrack->getZHelixE();
00010 mass=0.000511;
00011 }
00012 else if(pid==1){
00013 zhelix=mdcKalTrack->getZHelixMu();
00014 mass= 0.105658;
00015 }
00016 else if(pid==2){
00017 zhelix=mdcKalTrack->getZHelix();
00018 mass=0.139570;
00019 }
00020 else if(pid==3){
00021 zhelix=mdcKalTrack->getZHelixK();
00022 mass= 0.493677;
00023 }
00024 else{
00025 zhelix=mdcKalTrack->getZHelixP();
00026 mass= 0.938272;
00027 }
00028
00029 double dr(0),phi0(0),kappa(0),dz(0),tanl(0);
00030 dr=zhelix[0];
00031 phi0=zhelix[1];
00032 kappa=zhelix[2];
00033 dz=zhelix[3];
00034 tanl=zhelix[4];
00035
00036 int charge=0;
00037 if (kappa > 0.0000000001)
00038 charge = 1;
00039 else if (kappa < -0.0000000001)
00040 charge = -1;
00041
00042 double pxy=0;
00043 if(kappa!=0) pxy = 1.0/fabs(kappa);
00044
00045 double px = pxy * (-sin(phi0));
00046 double py = pxy * cos(phi0);
00047 double pz = pxy * tanl;
00048
00049 double e = sqrt( pxy*pxy + pz*pz + mass*mass );
00050
00051 return HepLorentzVector(px, py, pz, e);
00052
00053
00054 }
00055
00056
00057 HepLorentzVector utility::vfit(string channel, vector<int> kaonid, vector<int> pionid, HepPoint3D vx, EvtRecTrackIterator charged_begin){
00058
00059
00060 HepLorentzVector pchange(0,0,0,0);
00061
00062 int nvalid= kaonid.size()+pionid.size();
00063 if(nvalid<=1)
00064 return pchange;
00065
00066 HepSymMatrix Evx(3, 0);
00067 double bx = 1E+6; Evx[0][0] = bx*bx;
00068 double by = 1E+6; Evx[1][1] = by*by;
00069 double bz = 1E+6; Evx[2][2] = bz*bz;
00070 VertexParameter vxpar; vxpar.setVx(vx); vxpar.setEvx(Evx);
00071 double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
00072
00073
00074 VertexFit* vtxfit = VertexFit::instance();
00075 vtxfit->init();
00076
00077
00078 HepLorentzVector pold(0,0,0,0);
00079
00080 for(int i=0; i<kaonid.size();++i){
00081 EvtRecTrackIterator itTrk=charged_begin + kaonid[i];
00082 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
00083 pold+=getp4(mdcKalTrk, 3);
00084 WTrackParameter wtrk(xmass[3],mdcKalTrk->getZHelixK(),mdcKalTrk->getZErrorK() );
00085 vtxfit->AddTrack(i, wtrk);
00086 }
00087
00088 for(int i=0; i<pionid.size();++i){
00089 EvtRecTrackIterator itTrk=charged_begin + pionid[i];
00090 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
00091 pold+=getp4(mdcKalTrk, 2);
00092 WTrackParameter wtrk(xmass[2],mdcKalTrk->getZHelix(),mdcKalTrk->getZError() );
00093 vtxfit->AddTrack(kaonid.size()+i, wtrk);
00094 }
00095
00096 vector<int> trkIdxCol;
00097 trkIdxCol.clear();
00098 for (int i = 0; i < nvalid; i++)
00099 trkIdxCol.push_back(i);
00100 vtxfit->AddVertex(0, vxpar, trkIdxCol);
00101 if(!vtxfit->Fit(0))
00102 return pchange;
00103
00104 vtxfit->Swim(0);
00105
00106 HepLorentzVector pnew(0,0,0,0);
00107
00108 for(int i=0; i<nvalid;++i){
00109 WTrackParameter wtrk = vtxfit->wtrk(i);
00110 HepVector trk_val = HepVector(7,0);
00111 trk_val = wtrk.w();
00112 HepLorentzVector P_trk(trk_val[0],trk_val[1],trk_val[2],trk_val[3]);
00113 pnew+=P_trk;
00114 }
00115
00116 return (pnew-pold);
00117
00118 }
00119
00120 vector<double> utility::SecondaryVFit(EvtRecVeeVertex* ks, IVertexDbSvc* vtxsvc){
00121
00122
00123 double vfitchi2 = -999;
00124 double vfitlength = -999;
00125 double vfiterror = 999;
00126
00127 vector<double> results;
00128 results.push_back(vfitchi2);
00129 results.push_back(vfitlength);
00130 results.push_back(vfiterror);
00131
00132
00133
00134
00135
00136 EvtRecTrack* veeTrack1 = ks->pairDaughters().first;
00137 RecMdcKalTrack* veeKalTrack1 = veeTrack1->mdcKalTrack();
00138 WTrackParameter veeInitialWTrack1 = WTrackParameter(0.13957018, veeKalTrack1->getZHelix(), veeKalTrack1->getZError());
00139
00140 EvtRecTrack* veeTrack2 = ks->pairDaughters().second;
00141 RecMdcKalTrack* veeKalTrack2 = veeTrack2->mdcKalTrack();
00142 WTrackParameter veeInitialWTrack2 = WTrackParameter(0.13957018, veeKalTrack2->getZHelix(), veeKalTrack2->getZError());
00143
00144 VertexParameter wideVertex;
00145 HepPoint3D vWideVertex(0., 0., 0.);
00146 HepSymMatrix evWideVertex(3, 0);
00147
00148 evWideVertex[0][0] = 1.0e12;
00149 evWideVertex[1][1] = 1.0e12;
00150 evWideVertex[2][2] = 1.0e12;
00151
00152 wideVertex.setVx(vWideVertex);
00153 wideVertex.setEvx(evWideVertex);
00154
00155
00156 VertexFit* vtxfit = VertexFit::instance();
00157 vtxfit->init();
00158
00159
00160 vtxfit->AddTrack(0,veeInitialWTrack1);
00161 vtxfit->AddTrack(1,veeInitialWTrack2);
00162 vtxfit->AddVertex(0,wideVertex,0,1);
00163
00164
00165 vtxfit->Fit(0);
00166 vtxfit->Swim(0);
00167 vtxfit->BuildVirtualParticle(0);
00168
00169
00170 SecondVertexFit* svtxfit = SecondVertexFit::instance();
00171 svtxfit->init();
00172
00173
00174 VertexParameter beamSpot;
00175 HepPoint3D vBeamSpot(0., 0., 0.);
00176 HepSymMatrix evBeamSpot(3, 0);
00177
00178 if(vtxsvc->isVertexValid()){
00179 double* dbv = vtxsvc->PrimaryVertex();
00180 double* vv = vtxsvc->SigmaPrimaryVertex();
00181 for (unsigned int ivx = 0; ivx < 3; ivx++){
00182 vBeamSpot[ivx] = dbv[ivx];
00183 evBeamSpot[ivx][ivx] = vv[ivx] * vv[ivx];
00184 }
00185 }
00186 else{
00187 cout << "KSSELECTOR ERROR: Could not find a vertex from VertexDbSvc" << endl;
00188 return results;
00189 }
00190
00191 beamSpot.setVx(vBeamSpot);
00192 beamSpot.setEvx(evBeamSpot);
00193
00194 VertexParameter primaryVertex(beamSpot);
00195 svtxfit->setPrimaryVertex(primaryVertex);
00196
00197
00198 svtxfit->setVpar(vtxfit->vpar(0));
00199
00200
00201 svtxfit->AddTrack(0,vtxfit->wVirtualTrack(0));
00202
00203
00204
00205 if( !svtxfit->Fit() ) return results;
00206
00207
00208 vfitlength = svtxfit->decayLength();
00209 vfiterror = svtxfit->decayLengthError();
00210 vfitchi2 = svtxfit->chisq();
00211
00212 results.clear();
00213 results.push_back(vfitchi2);
00214 results.push_back(vfitlength);
00215 results.push_back(vfiterror);
00216
00217 return results;
00218 }