BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
utility.cxx
Go to the documentation of this file.
1#include "EvtRecEvent/EvtRecVeeVertex.h"
2#include "DTagAlg/LocalLambdaSelector.h"
3#include "DTagAlg/utility.h"
4
5HepLorentzVector utility::getp4(RecMdcKalTrack* mdcKalTrack, int pid){
6
7 HepVector zhelix;
8 double mass=0;
9
10 if(pid==0){
11 zhelix=mdcKalTrack->getZHelixE();
12 mass=0.000511;
13 }
14 else if(pid==1){
15 zhelix=mdcKalTrack->getZHelixMu();
16 mass= 0.105658;
17 }
18 else if(pid==2){
19 zhelix=mdcKalTrack->getZHelix();
20 mass=0.139570;
21 }
22 else if(pid==3){
23 zhelix=mdcKalTrack->getZHelixK();
24 mass= 0.493677;
25 }
26 else{
27 zhelix=mdcKalTrack->getZHelixP();
28 mass= 0.938272;
29 }
30
31 double dr(0),phi0(0),kappa(0),dz(0),tanl(0);
32 dr=zhelix[0];
33 phi0=zhelix[1];
34 kappa=zhelix[2];
35 dz=zhelix[3];
36 tanl=zhelix[4];
37
38 int charge=0;
39 if (kappa > 0.0000000001)
40 charge = 1;
41 else if (kappa < -0.0000000001)
42 charge = -1;
43
44 double pxy=0;
45 if(kappa!=0) pxy = 1.0/fabs(kappa);
46
47 double px = pxy * (-sin(phi0));
48 double py = pxy * cos(phi0);
49 double pz = pxy * tanl;
50
51 double e = sqrt( pxy*pxy + pz*pz + mass*mass );
52
53 return HepLorentzVector(px, py, pz, e);
54
55
56}
57
58
59HepLorentzVector utility::vfitref(string channel, vector<int> kaonid, vector<int> pionid, HepPoint3D vx, EvtRecTrackIterator charged_begin){
60 //use charged tracks only, except pions from Ks
61
62 HepLorentzVector pchange(0,0,0,0);
63
64 int nvalid= kaonid.size()+pionid.size();
65 if(nvalid<=1)
66 return pchange;
67
68 HepSymMatrix Evx(3, 0);
69 double bx = 1E+6; Evx[0][0] = bx*bx;
70 double by = 1E+6; Evx[1][1] = by*by;
71 double bz = 1E+6; Evx[2][2] = bz*bz;
72 VertexParameter vxpar; vxpar.setVx(vx); vxpar.setEvx(Evx);
73 double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
74
75
76 VertexFitRefine* vtxfit = VertexFitRefine::instance();
77 vtxfit->init();
78
79
80 HepLorentzVector pold(0,0,0,0);
81
82 for(int i=0; i<kaonid.size();++i){
83 EvtRecTrackIterator itTrk=charged_begin + kaonid[i];
84 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
85 pold+=getp4(mdcKalTrk, 3);
86 WTrackParameter wtrk(xmass[3],mdcKalTrk->getZHelixK(),mdcKalTrk->getZErrorK() );
87 vtxfit->AddTrack(i, mdcKalTrk, RecMdcKalTrack::kaon);
88 }
89
90 for(int i=0; i<pionid.size();++i){
91 EvtRecTrackIterator itTrk=charged_begin + pionid[i];
92 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
93 pold+=getp4(mdcKalTrk, 2);
94 WTrackParameter wtrk(xmass[2],mdcKalTrk->getZHelix(),mdcKalTrk->getZError() );
95 vtxfit->AddTrack(kaonid.size()+i, mdcKalTrk, RecMdcKalTrack::pion);
96 }
97
98 vector<int> trkIdxCol;
99 trkIdxCol.clear();
100 for (int i = 0; i < nvalid; i++)
101 trkIdxCol.push_back(i);
102 vtxfit->AddVertex(0, vxpar, trkIdxCol);
103
104 bool fitok = vtxfit->Fit();
105
106 if(!fitok){
107 return pchange;
108 }
109
110 HepLorentzVector pnew(0,0,0,0);
111
112 for(int i=0; i<nvalid;++i){
113 WTrackParameter wtrk = vtxfit->wtrk(i);
114 HepVector trk_val = HepVector(7,0);
115 trk_val = wtrk.w();
116 HepLorentzVector P_trk(trk_val[0],trk_val[1],trk_val[2],trk_val[3]);
117 pnew+=P_trk;
118 }
119
120 return (pnew-pold);
121
122}
123
124HepLorentzVector utility::vfitref(string channel, vector<int> kaonid, vector<int> pionid, vector<int> protonid, HepPoint3D vx, EvtRecTrackIterator charged_begin){
125 //use charged tracks only, except pions from Ks
126
127 HepLorentzVector pchange(0,0,0,0);
128
129 int nvalid= kaonid.size()+pionid.size()+protonid.size();
130 if(nvalid<=1)
131 return pchange;
132
133 HepSymMatrix Evx(3, 0);
134 double bx = 1E+6; Evx[0][0] = bx*bx;
135 double by = 1E+6; Evx[1][1] = by*by;
136 double bz = 1E+6; Evx[2][2] = bz*bz;
137 VertexParameter vxpar; vxpar.setVx(vx); vxpar.setEvx(Evx);
138 double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
139
140
141 VertexFitRefine* vtxfit = VertexFitRefine::instance();
142 vtxfit->init();
143
144
145 HepLorentzVector pold(0,0,0,0);
146
147 for(int i=0; i<kaonid.size();++i){
148 EvtRecTrackIterator itTrk=charged_begin + kaonid[i];
149 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
150 pold+=getp4(mdcKalTrk, 3);
151 WTrackParameter wtrk(xmass[3],mdcKalTrk->getZHelixK(),mdcKalTrk->getZErrorK() );
152 vtxfit->AddTrack(i, mdcKalTrk, RecMdcKalTrack::kaon);
153 }
154
155 for(int i=0; i<pionid.size();++i){
156 EvtRecTrackIterator itTrk=charged_begin + pionid[i];
157 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
158 pold+=getp4(mdcKalTrk, 2);
159 WTrackParameter wtrk(xmass[2],mdcKalTrk->getZHelix(),mdcKalTrk->getZError() );
160 vtxfit->AddTrack(kaonid.size()+i, mdcKalTrk, RecMdcKalTrack::pion);
161 }
162
163 for(int i=0; i<protonid.size();++i){
164 EvtRecTrackIterator itTrk=charged_begin + protonid[i];
165 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
166 pold+=getp4(mdcKalTrk, 4);
167 WTrackParameter wtrk(xmass[4],mdcKalTrk->getZHelixP(),mdcKalTrk->getZErrorP() );
168 vtxfit->AddTrack(kaonid.size()+pionid.size()+i, mdcKalTrk, RecMdcKalTrack::proton);
169 }
170
171 vector<int> trkIdxCol;
172 trkIdxCol.clear();
173 for (int i = 0; i < nvalid; i++)
174 trkIdxCol.push_back(i);
175 vtxfit->AddVertex(0, vxpar, trkIdxCol);
176
177 bool fitok = vtxfit->Fit();
178
179 if(!fitok){
180 return pchange;
181 }
182
183 HepLorentzVector pnew(0,0,0,0);
184
185 for(int i=0; i<nvalid;++i){
186 WTrackParameter wtrk = vtxfit->wtrk(i);
187 HepVector trk_val = HepVector(7,0);
188 trk_val = wtrk.w();
189 HepLorentzVector P_trk(trk_val[0],trk_val[1],trk_val[2],trk_val[3]);
190 pnew+=P_trk;
191 }
192
193 return (pnew-pold);
194
195}
196
197
198
199vector<double> utility::SecondaryVFitref(EvtRecVeeVertex* ks, IVertexDbSvc* vtxsvc){
200
201 // by default: chi2 = -999, length = -999, error = 999
202 double vfitchi2 = -999;
203 double vfitlength = -999;
204 double vfiterror = 999;
205
206 vector<double> results;
207 results.push_back(vfitchi2);
208 results.push_back(vfitlength);
209 results.push_back(vfiterror);
210
211 // --------------------------------------------------
212 // Do a secondary vertex fit and check the flight significance
213 // --------------------------------------------------
214
215 EvtRecTrack* veeTrack1 = ks->pairDaughters().first;
216 RecMdcKalTrack* veeKalTrack1 = veeTrack1->mdcKalTrack();
217 WTrackParameter veeInitialWTrack1 = WTrackParameter(0.13957018, veeKalTrack1->getZHelix(), veeKalTrack1->getZError());
218
219 EvtRecTrack* veeTrack2 = ks->pairDaughters().second;
220 RecMdcKalTrack* veeKalTrack2 = veeTrack2->mdcKalTrack();
221 WTrackParameter veeInitialWTrack2 = WTrackParameter(0.13957018, veeKalTrack2->getZHelix(), veeKalTrack2->getZError());
222
223 VertexParameter wideVertex;
224 HepPoint3D vWideVertex(0., 0., 0.);
225 HepSymMatrix evWideVertex(3, 0);
226
227 evWideVertex[0][0] = 1.0e12;
228 evWideVertex[1][1] = 1.0e12;
229 evWideVertex[2][2] = 1.0e12;
230
231 wideVertex.setVx(vWideVertex);
232 wideVertex.setEvx(evWideVertex);
233
234 // First, perform a vertex fit
235 VertexFitRefine* vtxfit = VertexFitRefine::instance();
236 vtxfit->init();
237
238 // add the daughters
239 vtxfit->AddTrack(0, veeKalTrack1, RecMdcKalTrack::pion);
240 vtxfit->AddTrack(1, veeKalTrack2, RecMdcKalTrack::pion);
241 vtxfit->AddVertex(0,wideVertex,0,1);
242
243 // do the fit
244 bool fitok = vtxfit->Fit();
245
246 // Now perform the secondary vertex fit
247 SecondVertexFit* svtxfit = SecondVertexFit::instance();
248 svtxfit->init();
249
250 // add the primary vertex
251 VertexParameter beamSpot;
252 HepPoint3D vBeamSpot(0., 0., 0.);
253 HepSymMatrix evBeamSpot(3, 0);
254
255 if(vtxsvc->isVertexValid()){
256 double* dbv = vtxsvc->PrimaryVertex();
257 double* vv = vtxsvc->SigmaPrimaryVertex();
258 for (unsigned int ivx = 0; ivx < 3; ivx++){
259 vBeamSpot[ivx] = dbv[ivx];
260 evBeamSpot[ivx][ivx] = vv[ivx] * vv[ivx];
261 }
262 }
263 else{
264 cout << "KSSELECTOR ERROR: Could not find a vertex from VertexDbSvc" << endl;
265 return results;
266 }
267
268 beamSpot.setVx(vBeamSpot);
269 beamSpot.setEvx(evBeamSpot);
270
271 VertexParameter primaryVertex(beamSpot);
272 svtxfit->setPrimaryVertex(primaryVertex);
273
274 // add the secondary vertex
275 svtxfit->setVpar(vtxfit->vpar(0));
276
277 // add the Ks or lambda
278 svtxfit->AddTrack(0,vtxfit->wVirtualTrack(0));
279
280 // do the second vertex fit
281 // if the fit fails, the default values will not be changed
282 if( !svtxfit->Fit() ) return results;
283
284 // save the new ks parameters
285 vfitlength = svtxfit->decayLength();
286 vfiterror = svtxfit->decayLengthError();
287 vfitchi2 = svtxfit->chisq();
288
289 results.clear();
290 results.push_back(vfitchi2);
291 results.push_back(vfitlength);
292 results.push_back(vfiterror);
293
294 return results;
295}
296
297vector<double> utility::SecondaryVFit_Lambdaref(EvtRecVeeVertex* lambda, IVertexDbSvc* vtxsvc){
298
299 // by default: chi2 = -999, length = -999, error = 999
300 double vfitchi2 = -999;
301 double vfitlength = -999;
302 double vfiterror = 999;
303
304 vector<double> results;
305 results.push_back(vfitchi2);
306 results.push_back(vfitlength);
307 results.push_back(vfiterror);
308
309 int index[2] = {0, 1};
310 if (lambda->vertexId() < 0){
311 index[0] = 1;
312 index[1] = 0;
313 }
314
315 EvtRecTrack* recTrk_proton = lambda->daughter(index[0]);
316 EvtRecTrack* recTrk_pion = lambda->daughter(index[1]);
317
318 // --------------------------------------------------
319 // Do a secondary vertex fit and check the flight significance
320 // --------------------------------------------------
321
322 RecMdcKalTrack* mdcKalTrk_proton = recTrk_proton->mdcKalTrack();
323 mdcKalTrk_proton->setPidType(RecMdcKalTrack::proton);
324 WTrackParameter WTrk_proton;
325 if(lambdaSelector.IfProtonPID()) WTrk_proton = WTrackParameter(0.9314940054, mdcKalTrk_proton->getZHelixP(), mdcKalTrk_proton->getZErrorP());
326 else WTrk_proton = WTrackParameter(0.9314940054, mdcKalTrk_proton->getZHelix() , mdcKalTrk_proton->getZError());
327
328 RecMdcKalTrack* mdcKalTrk_pion = recTrk_pion->mdcKalTrack();
329 mdcKalTrk_pion->setPidType(RecMdcKalTrack::pion);
330 WTrackParameter WTrk_pion = WTrackParameter(0.13957018, mdcKalTrk_pion->getZHelix(), mdcKalTrk_pion->getZError());
331
332 VertexParameter wideVertex;
333 HepPoint3D vWideVertex(0., 0., 0.);
334 HepSymMatrix evWideVertex(3, 0);
335
336 evWideVertex[0][0] = 1.0e12;
337 evWideVertex[1][1] = 1.0e12;
338 evWideVertex[2][2] = 1.0e12;
339
340 wideVertex.setVx(vWideVertex);
341 wideVertex.setEvx(evWideVertex);
342
343 // First, perform a vertex fit
344 VertexFitRefine* vtxfit = VertexFitRefine::instance();
345 vtxfit->init();
346
347 // add the daughters
348 vtxfit->AddTrack(0,mdcKalTrk_proton,RecMdcKalTrack::proton);
349 vtxfit->AddTrack(1,mdcKalTrk_pion,RecMdcKalTrack::pion);
350 vtxfit->AddVertex(0,wideVertex,0,1);
351
352 // do the fit
353 bool fitok = vtxfit->Fit();
354
355 // Now perform the secondary vertex fit
356 SecondVertexFit* svtxfit = SecondVertexFit::instance();
357 svtxfit->init();
358
359 // add the primary vertex
360 VertexParameter beamSpot;
361 HepPoint3D vBeamSpot(0., 0., 0.);
362 HepSymMatrix evBeamSpot(3, 0);
363
364 if(vtxsvc->isVertexValid()){
365 double* dbv = vtxsvc->PrimaryVertex();
366 double* vv = vtxsvc->SigmaPrimaryVertex();
367 for (unsigned int ivx = 0; ivx < 3; ivx++){
368 vBeamSpot[ivx] = dbv[ivx];
369 evBeamSpot[ivx][ivx] = vv[ivx] * vv[ivx];
370 }
371 }
372 else{
373 cout << "LambdaSELECTOR ERROR: Could not find a vertex from VertexDbSvc" << endl;
374 return results;
375 }
376
377 beamSpot.setVx(vBeamSpot);
378 beamSpot.setEvx(evBeamSpot);
379
380 VertexParameter primaryVertex(beamSpot);
381 svtxfit->setPrimaryVertex(primaryVertex);
382
383 // add the secondary vertex
384 svtxfit->setVpar(vtxfit->vpar(0));
385
386 // add the Ks or lambda
387 svtxfit->AddTrack(0,vtxfit->wVirtualTrack(0));
388
389 // do the second vertex fit
390 // if the fit fails, the default values will not be changed
391 if( !svtxfit->Fit() ) return results;
392
393 // save the new ks parameters
394 vfitlength = svtxfit->decayLength();
395 vfiterror = svtxfit->decayLengthError();
396 vfitchi2 = svtxfit->chisq();
397
398 results.clear();
399 results.push_back(vfitchi2);
400 results.push_back(vfitlength);
401 results.push_back(vfiterror);
402
403 return results;
404}
405
406HepLorentzVector utility::vfit(string channel, vector<int> kaonid, vector<int> pionid, HepPoint3D vx, EvtRecTrackIterator charged_begin){
407 //use charged tracks only, except pions from Ks
408
409 HepLorentzVector pchange(0,0,0,0);
410
411 int nvalid= kaonid.size()+pionid.size();
412 if(nvalid<=1)
413 return pchange;
414
415 HepSymMatrix Evx(3, 0);
416 double bx = 1E+6; Evx[0][0] = bx*bx;
417 double by = 1E+6; Evx[1][1] = by*by;
418 double bz = 1E+6; Evx[2][2] = bz*bz;
419 VertexParameter vxpar; vxpar.setVx(vx); vxpar.setEvx(Evx);
420 double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
421
422
423 VertexFit* vtxfit = VertexFit::instance();
424 vtxfit->init();
425
426
427 HepLorentzVector pold(0,0,0,0);
428
429 for(int i=0; i<kaonid.size();++i){
430 EvtRecTrackIterator itTrk=charged_begin + kaonid[i];
431 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
432 pold+=getp4(mdcKalTrk, 3);
433 WTrackParameter wtrk(xmass[3],mdcKalTrk->getZHelixK(),mdcKalTrk->getZErrorK() );
434 vtxfit->AddTrack(i, wtrk);
435 }
436
437 for(int i=0; i<pionid.size();++i){
438 EvtRecTrackIterator itTrk=charged_begin + pionid[i];
439 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
440 pold+=getp4(mdcKalTrk, 2);
441 WTrackParameter wtrk(xmass[2],mdcKalTrk->getZHelix(),mdcKalTrk->getZError() );
442 vtxfit->AddTrack(kaonid.size()+i, wtrk);
443 }
444
445 vector<int> trkIdxCol;
446 trkIdxCol.clear();
447 for (int i = 0; i < nvalid; i++)
448 trkIdxCol.push_back(i);
449 vtxfit->AddVertex(0, vxpar, trkIdxCol);
450 if(!vtxfit->Fit(0))
451 return pchange;
452
453 vtxfit->Swim(0);
454
455 HepLorentzVector pnew(0,0,0,0);
456
457 for(int i=0; i<nvalid;++i){
458 WTrackParameter wtrk = vtxfit->wtrk(i);
459 HepVector trk_val = HepVector(7,0);
460 trk_val = wtrk.w();
461 HepLorentzVector P_trk(trk_val[0],trk_val[1],trk_val[2],trk_val[3]);
462 pnew+=P_trk;
463 }
464
465 return (pnew-pold);
466
467}
468
469HepLorentzVector utility::vfit(string channel, vector<int> kaonid, vector<int> pionid, vector<int> protonid, HepPoint3D vx, EvtRecTrackIterator charged_begin){
470 //use charged tracks only, except pions from Ks
471
472 HepLorentzVector pchange(0,0,0,0);
473
474 int nvalid= kaonid.size()+pionid.size()+protonid.size();
475 if(nvalid<=1)
476 return pchange;
477
478 HepSymMatrix Evx(3, 0);
479 double bx = 1E+6; Evx[0][0] = bx*bx;
480 double by = 1E+6; Evx[1][1] = by*by;
481 double bz = 1E+6; Evx[2][2] = bz*bz;
482 VertexParameter vxpar; vxpar.setVx(vx); vxpar.setEvx(Evx);
483 double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
484
485
486 VertexFit* vtxfit = VertexFit::instance();
487 vtxfit->init();
488
489
490 HepLorentzVector pold(0,0,0,0);
491
492 for(int i=0; i<kaonid.size();++i){
493 EvtRecTrackIterator itTrk=charged_begin + kaonid[i];
494 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
495 pold+=getp4(mdcKalTrk, 3);
496 WTrackParameter wtrk(xmass[3],mdcKalTrk->getZHelixK(),mdcKalTrk->getZErrorK() );
497 vtxfit->AddTrack(i, wtrk);
498 }
499
500 for(int i=0; i<pionid.size();++i){
501 EvtRecTrackIterator itTrk=charged_begin + pionid[i];
502 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
503 pold+=getp4(mdcKalTrk, 2);
504 WTrackParameter wtrk(xmass[2],mdcKalTrk->getZHelix(),mdcKalTrk->getZError() );
505 vtxfit->AddTrack(kaonid.size()+i, wtrk);
506 }
507
508 for(int i=0; i<protonid.size();++i){
509 EvtRecTrackIterator itTrk=charged_begin + protonid[i];
510 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
511 pold+=getp4(mdcKalTrk, 4);
512 WTrackParameter wtrk(xmass[4],mdcKalTrk->getZHelixP(),mdcKalTrk->getZErrorP() );
513 vtxfit->AddTrack(kaonid.size()+pionid.size()+i, wtrk);
514 }
515
516 vector<int> trkIdxCol;
517 trkIdxCol.clear();
518 for (int i = 0; i < nvalid; i++)
519 trkIdxCol.push_back(i);
520 vtxfit->AddVertex(0, vxpar, trkIdxCol);
521 if(!vtxfit->Fit(0))
522 return pchange;
523
524 vtxfit->Swim(0);
525
526 HepLorentzVector pnew(0,0,0,0);
527
528 for(int i=0; i<nvalid;++i){
529 WTrackParameter wtrk = vtxfit->wtrk(i);
530 HepVector trk_val = HepVector(7,0);
531 trk_val = wtrk.w();
532 HepLorentzVector P_trk(trk_val[0],trk_val[1],trk_val[2],trk_val[3]);
533 pnew+=P_trk;
534 }
535
536 return (pnew-pold);
537
538}
539
540vector<double> utility::SecondaryVFit(EvtRecVeeVertex* ks, IVertexDbSvc* vtxsvc){
541
542 // by default: chi2 = -999, length = -999, error = 999
543 double vfitchi2 = -999;
544 double vfitlength = -999;
545 double vfiterror = 999;
546
547 vector<double> results;
548 results.push_back(vfitchi2);
549 results.push_back(vfitlength);
550 results.push_back(vfiterror);
551
552 // --------------------------------------------------
553 // Do a secondary vertex fit and check the flight significance
554 // --------------------------------------------------
555
556 EvtRecTrack* veeTrack1 = ks->pairDaughters().first;
557 RecMdcKalTrack* veeKalTrack1 = veeTrack1->mdcKalTrack();
558 WTrackParameter veeInitialWTrack1 = WTrackParameter(0.13957018, veeKalTrack1->getZHelix(), veeKalTrack1->getZError());
559
560 EvtRecTrack* veeTrack2 = ks->pairDaughters().second;
561 RecMdcKalTrack* veeKalTrack2 = veeTrack2->mdcKalTrack();
562 WTrackParameter veeInitialWTrack2 = WTrackParameter(0.13957018, veeKalTrack2->getZHelix(), veeKalTrack2->getZError());
563
564 VertexParameter wideVertex;
565 HepPoint3D vWideVertex(0., 0., 0.);
566 HepSymMatrix evWideVertex(3, 0);
567
568 evWideVertex[0][0] = 1.0e12;
569 evWideVertex[1][1] = 1.0e12;
570 evWideVertex[2][2] = 1.0e12;
571
572 wideVertex.setVx(vWideVertex);
573 wideVertex.setEvx(evWideVertex);
574
575 // First, perform a vertex fit
576 VertexFit* vtxfit = VertexFit::instance();
577 vtxfit->init();
578
579 // add the daughters
580 vtxfit->AddTrack(0,veeInitialWTrack1);
581 vtxfit->AddTrack(1,veeInitialWTrack2);
582 vtxfit->AddVertex(0,wideVertex,0,1);
583
584 // do the fit
585 vtxfit->Fit(0);
586 vtxfit->Swim(0);
587 vtxfit->BuildVirtualParticle(0);
588
589 // Now perform the secondary vertex fit
590 SecondVertexFit* svtxfit = SecondVertexFit::instance();
591 svtxfit->init();
592
593 // add the primary vertex
594 VertexParameter beamSpot;
595 HepPoint3D vBeamSpot(0., 0., 0.);
596 HepSymMatrix evBeamSpot(3, 0);
597
598 if(vtxsvc->isVertexValid()){
599 double* dbv = vtxsvc->PrimaryVertex();
600 double* vv = vtxsvc->SigmaPrimaryVertex();
601 for (unsigned int ivx = 0; ivx < 3; ivx++){
602 vBeamSpot[ivx] = dbv[ivx];
603 evBeamSpot[ivx][ivx] = vv[ivx] * vv[ivx];
604 }
605 }
606 else{
607 cout << "KSSELECTOR ERROR: Could not find a vertex from VertexDbSvc" << endl;
608 return results;
609 }
610
611 beamSpot.setVx(vBeamSpot);
612 beamSpot.setEvx(evBeamSpot);
613
614 VertexParameter primaryVertex(beamSpot);
615 svtxfit->setPrimaryVertex(primaryVertex);
616
617 // add the secondary vertex
618 svtxfit->setVpar(vtxfit->vpar(0));
619
620 // add the Ks or lambda
621 svtxfit->AddTrack(0,vtxfit->wVirtualTrack(0));
622
623 // do the second vertex fit
624 // if the fit fails, the default values will not be changed
625 if( !svtxfit->Fit() ) return results;
626
627 // save the new ks parameters
628 vfitlength = svtxfit->decayLength();
629 vfiterror = svtxfit->decayLengthError();
630 vfitchi2 = svtxfit->chisq();
631
632 results.clear();
633 results.push_back(vfitchi2);
634 results.push_back(vfitlength);
635 results.push_back(vfiterror);
636
637 return results;
638}
639
640vector<double> utility::SecondaryVFit_Lambda(EvtRecVeeVertex* lambda, IVertexDbSvc* vtxsvc){
641
642 // by default: chi2 = -999, length = -999, error = 999
643 double vfitchi2 = -999;
644 double vfitlength = -999;
645 double vfiterror = 999;
646
647 vector<double> results;
648 results.push_back(vfitchi2);
649 results.push_back(vfitlength);
650 results.push_back(vfiterror);
651
652 int index[2] = {0, 1};
653 if (lambda->vertexId() < 0){
654 index[0] = 1;
655 index[1] = 0;
656 }
657
658 EvtRecTrack* recTrk_proton = lambda->daughter(index[0]);
659 EvtRecTrack* recTrk_pion = lambda->daughter(index[1]);
660
661 // --------------------------------------------------
662 // Do a secondary vertex fit and check the flight significance
663 // --------------------------------------------------
664
665 RecMdcKalTrack* mdcKalTrk_proton = recTrk_proton->mdcKalTrack();
666 mdcKalTrk_proton->setPidType(RecMdcKalTrack::proton);
667 WTrackParameter WTrk_proton;
668 if(lambdaSelector.IfProtonPID()) WTrk_proton = WTrackParameter(0.9314940054, mdcKalTrk_proton->getZHelixP(), mdcKalTrk_proton->getZErrorP());
669 else WTrk_proton = WTrackParameter(0.9314940054, mdcKalTrk_proton->getZHelix() , mdcKalTrk_proton->getZError());
670
671 RecMdcKalTrack* mdcKalTrk_pion = recTrk_pion->mdcKalTrack();
672 mdcKalTrk_pion->setPidType(RecMdcKalTrack::pion);
673 WTrackParameter WTrk_pion = WTrackParameter(0.13957018, mdcKalTrk_pion->getZHelix(), mdcKalTrk_pion->getZError());
674
675 VertexParameter wideVertex;
676 HepPoint3D vWideVertex(0., 0., 0.);
677 HepSymMatrix evWideVertex(3, 0);
678
679 evWideVertex[0][0] = 1.0e12;
680 evWideVertex[1][1] = 1.0e12;
681 evWideVertex[2][2] = 1.0e12;
682
683 wideVertex.setVx(vWideVertex);
684 wideVertex.setEvx(evWideVertex);
685
686 // First, perform a vertex fit
687 VertexFit* vtxfit = VertexFit::instance();
688 vtxfit->init();
689
690 // add the daughters
691 vtxfit->AddTrack(0,WTrk_proton);
692 vtxfit->AddTrack(1,WTrk_pion);
693 vtxfit->AddVertex(0,wideVertex,0,1);
694
695 // do the fit
696 vtxfit->Fit(0);
697 vtxfit->Swim(0);
698 vtxfit->BuildVirtualParticle(0);
699
700 // Now perform the secondary vertex fit
701 SecondVertexFit* svtxfit = SecondVertexFit::instance();
702 svtxfit->init();
703
704 // add the primary vertex
705 VertexParameter beamSpot;
706 HepPoint3D vBeamSpot(0., 0., 0.);
707 HepSymMatrix evBeamSpot(3, 0);
708
709 if(vtxsvc->isVertexValid()){
710 double* dbv = vtxsvc->PrimaryVertex();
711 double* vv = vtxsvc->SigmaPrimaryVertex();
712 for (unsigned int ivx = 0; ivx < 3; ivx++){
713 vBeamSpot[ivx] = dbv[ivx];
714 evBeamSpot[ivx][ivx] = vv[ivx] * vv[ivx];
715 }
716 }
717 else{
718 cout << "LambdaSELECTOR ERROR: Could not find a vertex from VertexDbSvc" << endl;
719 return results;
720 }
721
722 beamSpot.setVx(vBeamSpot);
723 beamSpot.setEvx(evBeamSpot);
724
725 VertexParameter primaryVertex(beamSpot);
726 svtxfit->setPrimaryVertex(primaryVertex);
727
728 // add the secondary vertex
729 svtxfit->setVpar(vtxfit->vpar(0));
730
731 // add the Ks or lambda
732 svtxfit->AddTrack(0,vtxfit->wVirtualTrack(0));
733
734 // do the second vertex fit
735 // if the fit fails, the default values will not be changed
736 if( !svtxfit->Fit() ) return results;
737
738 // save the new ks parameters
739 vfitlength = svtxfit->decayLength();
740 vfiterror = svtxfit->decayLengthError();
741 vfitchi2 = svtxfit->chisq();
742
743 results.clear();
744 results.push_back(vfitchi2);
745 results.push_back(vfitlength);
746 results.push_back(vfiterror);
747
748 return results;
749}
750
751
752vector<double> utility::UpdatedKsIfo(EvtRecVeeVertex* ks, IVertexDbSvc* vtxsvc, bool m_useVFrefine){
753
754 // by default: px = -999, py = -999, pz = -999, E = -999, chisq = -999, ifok = -999.
755 double m_pxks = -999;
756 double m_pyks = -999;
757 double m_pzks = -999;
758 double m_Eks = -999;
759 double m_chisq = -999;
760 double m_ifok = -999;
761
762 vector<double> results;
763 results.clear();
764 results.push_back(m_pxks);
765 results.push_back(m_pyks);
766 results.push_back(m_pzks);
767 results.push_back(m_Eks);
768 results.push_back(m_chisq);
769 results.push_back(m_ifok);
770
771 EvtRecTrack* pion1Trk = ks->daughter(0);
772 RecMdcKalTrack* veeKalTrack1 = pion1Trk->mdcKalTrack();
773 veeKalTrack1->setPidType(RecMdcKalTrack::pion);
774 WTrackParameter veeInitialWTrack1 = WTrackParameter(0.13957018, veeKalTrack1->getZHelix(), veeKalTrack1->getZError());
775
776 EvtRecTrack* pion2Trk = ks->daughter(1);
777 RecMdcKalTrack* veeKalTrack2 = pion2Trk->mdcKalTrack();
778 veeKalTrack2->setPidType(RecMdcKalTrack::pion);
779 WTrackParameter veeInitialWTrack2 = WTrackParameter(0.13957018, veeKalTrack2->getZHelix(), veeKalTrack2->getZError());
780
781 //VertexParameter wideVertex;
782 HepPoint3D vWideVertex(0., 0., 0.);
783 HepSymMatrix evWideVertex(3, 0);
784
785 evWideVertex[0][0] = 1.0e12;
786 evWideVertex[1][1] = 1.0e12;
787 evWideVertex[2][2] = 1.0e12;
788
789 // add the primary vertex
790 HepPoint3D vBeamSpot(0., 0., 0.);
791 HepSymMatrix evBeamSpot(3, 0);
792
793 if(vtxsvc->isVertexValid()){
794 double* dbv = vtxsvc->PrimaryVertex();
795 double* vv = vtxsvc->SigmaPrimaryVertex();
796 for (unsigned int ivx = 0; ivx < 3; ivx++){
797 vBeamSpot[ivx] = dbv[ivx];
798 evBeamSpot[ivx][ivx] = vv[ivx] * vv[ivx];
799 }
800 }
801 else{
802 cout << "UpdateKs ERROR: Could not find a vertex from VertexDbSvc" << endl;
803 return results;
804 }
805
806 if(m_useVFrefine){
807
808 VertexParameter wideVertex;
809
810 wideVertex.setVx(vWideVertex);
811 wideVertex.setEvx(evWideVertex);
812
813 // First, perform a vertex fit refine
814 VertexFitRefine* vtxfitr = VertexFitRefine::instance();
815 vtxfitr->init();
816
817 vtxfitr->AddTrack(0, veeKalTrack1, RecMdcKalTrack::pion);
818 vtxfitr->AddTrack(1, veeKalTrack2, RecMdcKalTrack::pion);
819 vtxfitr->AddVertex(0, wideVertex, 0, 1);
820
821 bool fitok = vtxfitr->Fit();
822
823 if(fitok){
824 double chisq = vtxfitr->chisq(0);
825 HepLorentzVector pks = vtxfitr->pfit(0);
826
827 m_pxks = pks.px();
828 m_pyks = pks.py();
829 m_pzks = pks.pz();
830 m_Eks = pks.e();
831 m_chisq = chisq;
832 m_ifok = 1;
833
834 }
835
836 }else{
837 VertexParameter wideVertex;
838
839 wideVertex.setVx(vWideVertex);
840 wideVertex.setEvx(evWideVertex);
841
842 // First, perform a vertex fit
843 VertexFit* vtxfit = VertexFit::instance();
844 vtxfit->init();
845
846 // add the daughters
847 vtxfit->AddTrack(0,veeInitialWTrack1);
848 vtxfit->AddTrack(1,veeInitialWTrack2);
849 vtxfit->AddVertex(0,wideVertex,0,1);
850
851 // do the fit
852 if(vtxfit->Fit(0)){
853 vtxfit->Swim(0);
854 vtxfit->BuildVirtualParticle(0);
855
856 WTrackParameter wks = vtxfit->wVirtualTrack(0);
857 double chisq = vtxfit->chisq(0);
858 HepLorentzVector pks = wks.p();
859
860 m_pxks = pks.px();
861 m_pyks = pks.py();
862 m_pzks = pks.pz();
863 m_Eks = pks.e();
864 m_chisq = chisq;
865 m_ifok = 1;
866
867 }
868
869 }
870
871 results.clear();
872 results.push_back(m_pxks);
873 results.push_back(m_pyks);
874 results.push_back(m_pzks);
875 results.push_back(m_Eks);
876 results.push_back(m_chisq);
877 results.push_back(m_ifok);
878
879}
880
881
882vector<double> utility::UpdatedLambdaIfo(EvtRecVeeVertex* lambda, IVertexDbSvc* vtxsvc, bool m_useVFrefine){
883
884 // by default: px = -999, py = -999, pz = -999, E = -999, chisq = -999, ifok = -999.
885 double m_pxlambda = -999;
886 double m_pylambda = -999;
887 double m_pzlambda = -999;
888 double m_Elambda = -999;
889 double m_chisq = -999;
890 double m_ifok = -999;
891
892 vector<double> results;
893 results.clear();
894 results.push_back(m_pxlambda);
895 results.push_back(m_pylambda);
896 results.push_back(m_pzlambda);
897 results.push_back(m_Elambda);
898 results.push_back(m_chisq);
899 results.push_back(m_ifok);
900
901 int index[2] = {0, 1};
902 if(lambda->vertexId() < 0){
903 index[0] = 1;
904 index[1] = 0;
905 }
906
907 EvtRecTrack* protonTrk = lambda->daughter(index[0]);
908 RecMdcKalTrack* mdcKalTrk_proton = protonTrk->mdcKalTrack();
909 mdcKalTrk_proton->setPidType(RecMdcKalTrack::proton);
910 WTrackParameter WTrk_proton = WTrackParameter(0.9314940054, mdcKalTrk_proton->getZHelixP(), mdcKalTrk_proton->getZErrorP());
911
912 EvtRecTrack* pionTrk = lambda->daughter(index[1]);
913 RecMdcKalTrack* mdcKalTrk_pion = pionTrk->mdcKalTrack();
914 mdcKalTrk_pion->setPidType(RecMdcKalTrack::pion);
915 WTrackParameter WTrk_pion = WTrackParameter(0.13957018, mdcKalTrk_pion->getZHelix(), mdcKalTrk_pion->getZError());
916
917 //VertexParameter wideVertex;
918 HepPoint3D vWideVertex(0., 0., 0.);
919 HepSymMatrix evWideVertex(3, 0);
920
921 evWideVertex[0][0] = 1.0e12;
922 evWideVertex[1][1] = 1.0e12;
923 evWideVertex[2][2] = 1.0e12;
924
925 // add the primary vertex
926 HepPoint3D vBeamSpot(0., 0., 0.);
927 HepSymMatrix evBeamSpot(3, 0);
928
929 if(vtxsvc->isVertexValid()){
930 double* dbv = vtxsvc->PrimaryVertex();
931 double* vv = vtxsvc->SigmaPrimaryVertex();
932 for (unsigned int ivx = 0; ivx < 3; ivx++){
933 vBeamSpot[ivx] = dbv[ivx];
934 evBeamSpot[ivx][ivx] = vv[ivx] * vv[ivx];
935 }
936 }
937 else{
938 cout << "UpdateLambda ERROR: Could not find a vertex from VertexDbSvc" << endl;
939 return results;
940 }
941
942 if(m_useVFrefine){
943
944 VertexParameter wideVertex;
945
946 wideVertex.setVx(vWideVertex);
947 wideVertex.setEvx(evWideVertex);
948
949 // First, perform a vertex fit refine
950 VertexFitRefine* vtxfitr = VertexFitRefine::instance();
951 vtxfitr->init();
952
953 vtxfitr->AddTrack(0, mdcKalTrk_proton, RecMdcKalTrack::proton);
954 vtxfitr->AddTrack(1, mdcKalTrk_pion, RecMdcKalTrack::pion);
955 vtxfitr->AddVertex(0, wideVertex, 0, 1);
956
957 bool fitok = vtxfitr->Fit();
958
959 if(fitok){
960 double chisq = vtxfitr->chisq(0);
961 HepLorentzVector plambda = vtxfitr->pfit(0);
962
963 m_pxlambda = plambda.px();
964 m_pylambda = plambda.py();
965 m_pzlambda = plambda.pz();
966 m_Elambda = plambda.e();
967 m_chisq = chisq;
968 m_ifok = 1;
969
970 }
971
972 }else{
973 VertexParameter wideVertex;
974
975 wideVertex.setVx(vWideVertex);
976 wideVertex.setEvx(evWideVertex);
977
978 // First, perform a vertex fit
979 VertexFit* vtxfit = VertexFit::instance();
980 vtxfit->init();
981
982 // add the daughters
983 vtxfit->AddTrack(0,WTrk_proton);
984 vtxfit->AddTrack(1,WTrk_pion);
985 vtxfit->AddVertex(0,wideVertex,0,1);
986
987 // do the fit
988 if(vtxfit->Fit(0)){
989 vtxfit->Swim(0);
990 vtxfit->BuildVirtualParticle(0);
991
992 WTrackParameter wlambda = vtxfit->wVirtualTrack(0);
993 double chisq = vtxfit->chisq(0);
994 HepLorentzVector plambda = wlambda.p();
995
996 m_pxlambda = plambda.px();
997 m_pylambda = plambda.py();
998 m_pzlambda = plambda.pz();
999 m_Elambda = plambda.e();
1000 m_chisq = chisq;
1001 m_ifok = 1;
1002
1003 }
1004
1005 }
1006
1007 results.clear();
1008 results.push_back(m_pxlambda);
1009 results.push_back(m_pylambda);
1010 results.push_back(m_pzlambda);
1011 results.push_back(m_Elambda);
1012 results.push_back(m_chisq);
1013 results.push_back(m_ifok);
1014
1015}
1016
1017
1018
1019
sin
double sin(const BesAngle a)
Definition BesAngle.h:210
cos
double cos(const BesAngle a)
Definition BesAngle.h:213
mass
double mass
Definition CosmicGenerator.cxx:138
EvtRecTrackIterator
EvtRecTrackCol::iterator EvtRecTrackIterator
Definition EvtRecTrack.h:111
xmass
const double xmass[5]
Definition Gam4pikp.cxx:50
lambdaSelector
LocalLambdaSelector lambdaSelector
Definition LocalLambdaSelector.cxx:269
DstMdcKalTrack::setPidType
static void setPidType(PidType pidType)
Definition DstMdcKalTrack.h:59
EvtRecTrack::mdcKalTrack
RecMdcKalTrack * mdcKalTrack()
Definition EvtRecTrack.h:54
EvtRecVeeVertex::pairDaughters
std::pair< SmartRef< EvtRecTrack >, SmartRef< EvtRecTrack > > & pairDaughters()
Definition EvtRecVeeVertex.h:46
EvtRecVeeVertex::daughter
SmartRef< EvtRecTrack > & daughter(int i)
Definition EvtRecVeeVertex.h:89
EvtRecVeeVertex::vertexId
int vertexId() const
Definition EvtRecVeeVertex.h:35
IVertexDbSvc::isVertexValid
virtual bool isVertexValid()=0
IVertexDbSvc::SigmaPrimaryVertex
virtual double * SigmaPrimaryVertex()=0
IVertexDbSvc::PrimaryVertex
virtual double * PrimaryVertex()=0
RecMdcKalTrack::getZHelix
const HepVector & getZHelix() const
Definition RecMdcKalTrack.h:66
RecMdcKalTrack::getZHelixP
HepVector & getZHelixP()
Definition RecMdcKalTrack.h:98
RecMdcKalTrack::getZHelixK
HepVector & getZHelixK()
Definition RecMdcKalTrack.h:91
RecMdcKalTrack::getZHelixE
HepVector & getZHelixE()
Definition RecMdcKalTrack.h:77
RecMdcKalTrack::getZHelixMu
HepVector & getZHelixMu()
Definition RecMdcKalTrack.h:84
RecMdcKalTrack::getZError
const HepSymMatrix & getZError() const
Definition RecMdcKalTrack.h:67
RecMdcKalTrack::getZErrorP
HepSymMatrix & getZErrorP()
Definition RecMdcKalTrack.h:99
RecMdcKalTrack::getZErrorK
HepSymMatrix & getZErrorK()
Definition RecMdcKalTrack.h:92
SecondVertexFit::setPrimaryVertex
void setPrimaryVertex(const VertexParameter vpar)
Definition SecondVertexFit.h:37
SecondVertexFit::decayLength
double decayLength() const
Definition SecondVertexFit.h:46
SecondVertexFit::decayLengthError
double decayLengthError() const
Definition SecondVertexFit.h:47
SecondVertexFit::instance
static SecondVertexFit * instance()
Definition SecondVertexFit.cxx:45
SecondVertexFit::setVpar
void setVpar(const VertexParameter vpar)
Definition SecondVertexFit.h:35
SecondVertexFit::chisq
double chisq() const
Definition SecondVertexFit.h:45
TrackPool::AddTrack
void AddTrack(const int number, const double mass, const RecMdcTrack *trk)
Definition TrackPool.cxx:22
VertexFitRefine::wVirtualTrack
WTrackParameter wVirtualTrack(int n) const
Definition VertexFitRefine.cc:177
VertexFitRefine::pfit
HepLorentzVector pfit(int n) const
Definition VertexFitRefine.cc:127
VertexFitRefine::chisq
double chisq() const
Definition VertexFitRefine.cc:117
VertexFitRefine::AddTrack
void AddTrack(const int index, RecMdcKalTrack *p, const RecMdcKalTrack::PidType pid)
Definition VertexFitRefine.cxx:78
VertexFitRefine::wtrk
WTrackParameter wtrk(int n) const
Definition VertexFitRefine.cc:147
VertexFitRefine::instance
static VertexFitRefine * instance()
Definition VertexFitRefine.cxx:34
VertexFitRefine::AddVertex
void AddVertex(int number, VertexParameter vpar, std::vector< int > lis)
Definition VertexFitRefine.cc:17
VertexFitRefine::vpar
VertexParameter vpar(int n) const
Definition VertexFitRefine.cc:172
VertexFit::wtrk
WTrackParameter wtrk(int n) const
Definition VertexFit.h:79
VertexFit::chisq
double chisq() const
Definition VertexFit.h:66
VertexFit::wVirtualTrack
WTrackParameter wVirtualTrack(int n) const
Definition VertexFit.h:92
VertexFit::init
void init()
Definition VertexFit.cxx:29
VertexFit::AddVertex
void AddVertex(int number, VertexParameter vpar, std::vector< int > lis)
Definition VertexFit.cxx:89
VertexFit::instance
static VertexFit * instance()
Definition VertexFit.cxx:15
VertexFit::vpar
VertexParameter vpar(int n) const
Definition VertexFit.h:89
VertexFit::BuildVirtualParticle
void BuildVirtualParticle(int number)
Definition VertexFit.cxx:619
VertexFit::Swim
void Swim(int n)
Definition VertexFit.h:59
VertexFit::Fit
bool Fit()
Definition VertexFit.cxx:301
VertexParameter::setEvx
void setEvx(const HepSymMatrix &eVx)
Definition VertexParameter.h:29
VertexParameter::setVx
void setVx(const HepPoint3D &vx)
Definition VertexParameter.h:54
WTrackParameter::p
HepLorentzVector p() const
Definition WTrackParameter.h:69
WTrackParameter::w
HepVector w() const
Definition WTrackParameter.h:67
utility::vfitref
HepLorentzVector vfitref(string channel, vector< int > kaonid, vector< int > pionid, HepPoint3D vx, EvtRecTrackIterator charged_begin)
Definition utility.cxx:59
utility::SecondaryVFit_Lambdaref
vector< double > SecondaryVFit_Lambdaref(EvtRecVeeVertex *lambda, IVertexDbSvc *vtxsvc)
Definition utility.cxx:297
utility::SecondaryVFit_Lambda
vector< double > SecondaryVFit_Lambda(EvtRecVeeVertex *lambda, IVertexDbSvc *vtxsvc)
Definition utility.cxx:640
utility::SecondaryVFitref
vector< double > SecondaryVFitref(EvtRecVeeVertex *ks, IVertexDbSvc *vtxsvc)
Definition utility.cxx:199
utility::UpdatedKsIfo
vector< double > UpdatedKsIfo(EvtRecVeeVertex *ks, IVertexDbSvc *vtxsvc, bool m_useVFrefine)
Definition utility.cxx:752
utility::UpdatedLambdaIfo
vector< double > UpdatedLambdaIfo(EvtRecVeeVertex *lambda, IVertexDbSvc *vtxsvc, bool m_useVFrefine)
Definition utility.cxx:882
utility::vfit
HepLorentzVector vfit(string channel, vector< int > kaonid, vector< int > pionid, HepPoint3D vx, EvtRecTrackIterator charged_begin)
Definition utility.cxx:406
utility::getp4
HepLorentzVector getp4(RecMdcKalTrack *mdcKalTrack, int pid)
Definition utility.cxx:5
utility::SecondaryVFit
vector< double > SecondaryVFit(EvtRecVeeVertex *ks, IVertexDbSvc *vtxsvc)
Definition utility.cxx:540
charge
float charge
Definition prepare_data.cxx:28