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