BOSS 6.6.4.p03
BESIII Offline Software System
Loading...
Searching...
No Matches
DTagTool.cxx
Go to the documentation of this file.
1#include "GaudiKernel/MsgStream.h"
2#include "GaudiKernel/AlgFactory.h"
3#include "GaudiKernel/SmartDataPtr.h"
4#include "GaudiKernel/IDataProviderSvc.h"
5#include "GaudiKernel/PropertyMgr.h"
6#include "GaudiKernel/Service.h"
7#include "GaudiKernel/ISvcLocator.h"
8#include "GaudiKernel/SvcFactory.h"
9
15#include "EventModel/Event.h"
16
17#include "McTruth/McParticle.h"
18
19#include "TMath.h"
20#include "GaudiKernel/INTupleSvc.h"
21#include "GaudiKernel/NTuple.h"
22#include "GaudiKernel/Bootstrap.h"
23#include "GaudiKernel/IHistogramSvc.h"
24
26#include "VertexFit/Helix.h"
27#include "VertexFit/VertexFit.h"
28
29#include "DTagTool/DTagTool.h"
30
31
32DTagTool::DTagTool() : m_evtSvc(0), m_iterbegin(0), m_iterend(0),
33 m_iterstag(0), m_iterdtag1(0), m_iterdtag2(0),
34 m_chargebegin(0),m_chargeend(0),m_showerbegin(0),m_showerend(0){
35
36 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), "/Event/EvtRec/EvtRecEvent");
37 if ( ! evtRecEvent ) {
38 cout << MSG::FATAL << "Could not find EvtRecEvent" << endl;
39 exit(1);
40 }
41
42 SmartDataPtr<EvtRecTrackCol> evtRecTrackCol( eventSvc(), "/Event/EvtRec/EvtRecTrackCol");
43 if ( ! evtRecTrackCol ) {
44 cout << MSG::FATAL << "Could not find EvtRecTrackCol" << endl;
45 exit(1);
46 }
47
48
49 StatusCode sc = Gaudi::svcLocator()->service("SimplePIDSvc", m_simplePIDSvc);
50 if ( sc.isFailure() )
51 {
52 //log << MSG::FATAL << "Could not load SimplePIDSvc!" << endreq;
53 exit(1);
54 }
55
56
57
58 m_chargebegin=evtRecTrackCol->begin();
59 m_chargeend=evtRecTrackCol->begin()+evtRecEvent->totalCharged();
60 m_showerbegin=evtRecTrackCol->begin()+evtRecEvent->totalCharged();
61 m_showerend=evtRecTrackCol->begin()+evtRecEvent->totalTracks();
62
63 /// Accessing DTagList
64 SmartDataPtr<EvtRecDTagCol> evtRecDTagCol(eventSvc(), EventModel::EvtRec::EvtRecDTagCol);
65 if ( ! evtRecDTagCol ) {
66 cout << "Could not find EvtRecDTagCol" << endl;
67 exit(1);
68 }
69
70 /// Accessing Ks list
71 SmartDataPtr<EvtRecVeeVertexCol> evtRecVeeVertexCol(eventSvc(), "/Event/EvtRec/EvtRecVeeVertexCol");
72 if ( ! evtRecVeeVertexCol ) {
73 cout<< "Could not find EvtRecVeeVertexCol" << endl;
74 exit(1);
75 }
76
77 /// Accessing pi0 list
78 SmartDataPtr<EvtRecPi0Col> recPi0Col(eventSvc(), "/Event/EvtRec/EvtRecPi0Col");
79 if ( ! recPi0Col ) {
80 cout<< "Could not find EvtRecPi0Col" << endl;
81 exit(1);
82 }
83
84
85
86 m_iterbegin=evtRecDTagCol->begin();
87 m_iterend=evtRecDTagCol->end();
88 m_pi0iterbegin=recPi0Col->begin();
89 m_pi0iterend=recPi0Col->end();
90 m_ksiterbegin=evtRecVeeVertexCol->begin();
91 m_ksiterend=evtRecVeeVertexCol->end();
92
93 if(evtRecDTagCol->size() > 0)
94 m_isdtaglistempty=false;
95 else
96 m_isdtaglistempty=true;
97
98 //set initial pid requirement
99 m_pid=true;
100
101
102 //fill d0, dp, ds modes seprately
103 int index=0;
104 for(EvtRecDTagCol::iterator iter=m_iterbegin; iter != m_iterend; iter++){
105
106 if( (int)( (*iter)->decayMode())< 200 )
107 m_d0modes.push_back( index );
108 else if( (int)( (*iter)->decayMode())< 400 )
109 m_dpmodes.push_back( index );
110 else if( (int)( (*iter)->decayMode())< 1000 )
111 m_dsmodes.push_back( index );
112
113 index++;
114 }
115
116}
117
119
120 m_d0modes.clear();
121 m_dpmodes.clear();
122 m_dsmodes.clear();
123
124}
125
127
128 m_d0modes.clear();
129 m_dpmodes.clear();
130 m_dsmodes.clear();
131
132}
133
134
136
137 vector<int> mode;
138 int index=0;
139 for(EvtRecDTagCol::iterator iter=m_iterbegin; iter != m_iterend; iter++){
140
141 if(m_pid){
142 if( (*iter)->decayMode() == decaymode && (*iter)->type() == 1 )
143 mode.push_back( index );
144 }
145 else{
146 if( (*iter)->decayMode() == decaymode )
147 mode.push_back( index );
148 }
149
150 index++;
151 }
152
153 return mode;
154}
155
156
157vector<int> DTagTool::mode(int decaymode){
158
159 vector<int> mode;
160 int index=0;
161 for(EvtRecDTagCol::iterator iter=m_iterbegin; iter != m_iterend; iter++){
162
163 if(m_pid){
164 if( (*iter)->decayMode() == decaymode && (*iter)->type() == 1 )
165 mode.push_back( index );
166 }
167 else{
168 if( (*iter)->decayMode() == decaymode )
169 mode.push_back( index );
170 }
171
172 index++;
173 }
174
175 return mode;
176}
177
178
179
180//search single tag with charm
182
183 return findSTag(static_cast<int>(mode), tagcharm);
184
185}//end of stag
186
187
188//search single tag without charm
190 return findSTag(static_cast<int>(mode));
191
192}//end of stag
193
194
195
196
197//use integer as argument
198bool DTagTool::findSTag(int mode, int tagcharm){
199
200 bool isStcand=false;
201 double de_min=1;
202
203 //loop over the dtag list
204 EvtRecDTagCol::iterator iter_dtag=m_iterbegin;
205 for ( ; iter_dtag != m_iterend; iter_dtag++){
206
207 if(m_pid){
208 if( (*iter_dtag)->type()!=1 || (*iter_dtag)->decayMode() != mode || (*iter_dtag)->charm() != tagcharm )
209 continue;
210 }
211 else{
212 if( (*iter_dtag)->decayMode() != mode || (*iter_dtag)->charm() != tagcharm )
213 continue;
214 }
215 if(fabs((*iter_dtag)->deltaE())<fabs(de_min)){
216 isStcand=true;
217 m_iterstag=iter_dtag;
218 de_min=(*iter_dtag)->deltaE();
219 }
220
221 } //end of looping over the entire DTag list
222
223 return isStcand;
224
225}//end of stag
226
227
228bool DTagTool::findSTag(int mode){
229
230 bool isStcand=false;
231 double de_min=1;
232
233 //loop over the dtag list
234 EvtRecDTagCol::iterator iter_dtag=m_iterbegin;
235 for ( ; iter_dtag != m_iterend; iter_dtag++){
236
237 if(m_pid){
238 if( (*iter_dtag)->type()!=1 || (*iter_dtag)->decayMode() != mode )
239 continue;
240 }
241 else{
242 if( (*iter_dtag)->decayMode() != mode )
243 continue;
244 }
245
246 if(fabs((*iter_dtag)->deltaE())<fabs(de_min)){
247 isStcand=true;
248 m_iterstag=iter_dtag;
249 de_min=(*iter_dtag)->deltaE();
250 }
251
252 } //end of looping over the entire DTag list
253
254 return isStcand;
255
256}//end of stag
257
258
259////////////////////////////////////////////
260////////////////////////////////////////////
261
262
263
264//generic double tagging searches
266
267 return findDTag(static_cast<int>(mode1), static_cast<int>(mode2), smass);
268
269}
270
271bool DTagTool::findDTag(EvtRecDTag::DecayMode mode1, int tagcharm1, EvtRecDTag::DecayMode mode2, int tagcharm2, string smass){
272
273 return findDTag(static_cast<int>(mode1), tagcharm1, static_cast<int>(mode2), tagcharm2, smass);
274
275} //end of findDtag
276
277
278//find all the double tags
280
281 return findADTag(static_cast<int>(mode1), static_cast<int>(mode2));
282
283}
284
285bool DTagTool::findADTag(EvtRecDTag::DecayMode mode1, int tagcharm1, EvtRecDTag::DecayMode mode2, int tagcharm2){
286
287 return findADTag(static_cast<int>(mode1), tagcharm1, static_cast<int>(mode2), tagcharm2);
288
289} //end of findDtag
290
291
292
293
294//use integer as argument
295bool DTagTool::findDTag(int mode1, int mode2, string smass){
296
297 int tagcharm1= (mode1<10 || mode1>=200)?+1:0;
298 int tagcharm2= (mode2<10 || mode2>=200)?-1:0;
299
300 if(tagcharm1*tagcharm2>0){
301 cout<<"double tagging warning! two modes can't have same nonzero charmness"<<endl;
302 return false;
303 }
304
305 //define D candidate mass
306 double mDcand=0;
307 if( mode1 < 200 && mode2 < 200)
308 mDcand = 1.8645;
309 else if ( mode1>=200 && mode1 < 400 && mode2>=200 && mode2 < 400)
310 mDcand = 1.8693;
311 else if ( mode1>=400 && mode1 < 1000 && mode2>=400 && mode2 < 1000)
312 mDcand = 1.9682;
313 else{
314 cout<<"double tag modes are not from same particle ! "<<endl;
315 return false;
316 }
317
318
319 vector<int> igood1, igood2;
320 igood1.clear(),igood2.clear();
321 int index=0;
322 EvtRecDTagCol::iterator iter_dtag=m_iterbegin;
323
324 //charge conjucation considered
325 for ( ; iter_dtag != m_iterend; iter_dtag++){
326 int iter_mode=(*iter_dtag)->decayMode();
327 int iter_charm=(*iter_dtag)->charm();
328 int iter_type=(*iter_dtag)->type();
329
330 if(m_pid){
331 if( iter_mode == mode1 && iter_charm == tagcharm1 && iter_type==1 )
332 igood1.push_back(index);
333 if( tagcharm1!=0 && iter_mode == mode1 && iter_charm == -tagcharm1 && iter_type==1 )
334 igood1.push_back(index);
335
336 if( iter_mode == mode2 && iter_charm == tagcharm2 && iter_type==1 )
337 igood2.push_back(index);
338 if( tagcharm2!=0 && iter_mode == mode2 && iter_charm == -tagcharm2 && iter_type==1 )
339 igood2.push_back(index);
340 }
341 else{
342 if( iter_mode == mode1 && iter_charm == tagcharm1 )
343 igood1.push_back(index);
344 if( tagcharm1!=0 && iter_mode == mode1 && iter_charm == -tagcharm1 )
345 igood1.push_back(index);
346
347 if( iter_mode == mode2 && iter_charm == tagcharm2 )
348 igood2.push_back(index);
349 if( tagcharm2!=0 && iter_mode == mode2 && iter_charm == -tagcharm2 )
350 igood2.push_back(index);
351 }
352
353 index++;
354 }
355
356 //look for the best pair of double-tagged event
357
358 bool isDtcand=false;
359 double deltaM=1.00;
360 int index_i=0, index_j=0;
361 EvtRecDTagCol::iterator iter_i, iter_j;
362
363 for(int i=0; i<igood1.size(); i++){
364
365 iter_i=m_iterbegin+igood1[i];
366 double mass_i=(*iter_i)->mBC();
367 if(smass=="inv")
368 mass_i=(*iter_i)->mass();
369 int charm_i=(*iter_i)->charm();
370 for(int j=0;j<igood2.size();j++){
371 iter_j=m_iterbegin+igood2[j];
372 double mass_j=(*iter_j)->mBC();
373 if(smass=="inv")
374 mass_j=(*iter_j)->mass();
375 int charm_j=(*iter_j)->charm();
376 if( charm_i*charm_j>0 || igood2[j] == igood1[i] ) continue;
377
378 if(shareTracks(iter_i,iter_j)) continue;
379
380 if( fabs(0.5*(mass_i+mass_j)-mDcand) < deltaM){
381 deltaM = fabs(0.5*(mass_i+mass_j)-mDcand);
382 index_i = i;
383 index_j = j;
384 isDtcand = true;
385 }
386
387 } //end of j loop
388 } //end of i loop
389
390 if(isDtcand){
391 m_iterdtag1=m_iterbegin+igood1[index_i];
392 m_iterdtag2=m_iterbegin+igood2[index_j];
393 }
394
395 return isDtcand;
396}
397
398
399bool DTagTool::findDTag(int mode1, int tagcharm1, int mode2, int tagcharm2, string smass){
400
401 if(tagcharm1*tagcharm2>0){
402 cout<<"double tagging warning! two modes can't have same nonzero charmness"<<endl;
403 return false;
404 }
405
406 //define D candidate mass
407 double mDcand=0;
408 if( mode1 < 200 && mode2 < 200)
409 mDcand = 1.8645;
410 else if (mode1>=200 && mode1 < 400 && mode2>=200 && mode2 < 400)
411 mDcand = 1.8693;
412 else if (mode1>=400 && mode1 < 1000 && mode2>=400 && mode2 < 1000)
413 mDcand = 1.9682;
414 else{
415 cout<<"double tag modes are not from same particle ! "<<endl;
416 return false;
417 }
418
419
420 vector<int> igood1, igood2;
421 igood1.clear(),igood2.clear();
422 int index=0;
423 EvtRecDTagCol::iterator iter_dtag=m_iterbegin;
424
425 for ( ; iter_dtag != m_iterend; iter_dtag++){
426 int iter_mode=(*iter_dtag)->decayMode();
427 int iter_charm=(*iter_dtag)->charm();
428 int iter_type=(*iter_dtag)->type();
429
430 if(m_pid){
431 if( iter_mode == mode1 && iter_charm == tagcharm1 && iter_type==1 )
432 igood1.push_back(index);
433
434 if( iter_mode == mode2 && iter_charm == tagcharm2 && iter_type==1)
435 igood2.push_back(index);
436 }
437 else{
438 if( iter_mode == mode1 && iter_charm == tagcharm1)
439 igood1.push_back(index);
440
441 if( iter_mode == mode2 && iter_charm == tagcharm2)
442 igood2.push_back(index);
443 }
444
445
446 index++;
447 }
448
449 //look for the best pair of double-tagged event
450
451 bool isDtcand=false;
452 double deltaM=1.00;
453 int index_i=0, index_j=0;
454 EvtRecDTagCol::iterator iter_i, iter_j;
455
456 for(int i=0; i<igood1.size(); i++){
457
458 iter_i=m_iterbegin+igood1[i];
459 double mass_i=(*iter_i)->mBC();
460 if(smass=="inv")
461 mass_i=(*iter_i)->mass();
462 int charm_i=(*iter_i)->charm();
463 for(int j=0;j<igood2.size();j++){
464 iter_j=m_iterbegin+igood2[j];
465 double mass_j=(*iter_j)->mBC();
466 if(smass=="inv")
467 mass_j=(*iter_j)->mass();
468 int charm_j=(*iter_j)->charm();
469 if( charm_i*charm_j>0 || igood2[j] == igood1[i] ) continue;
470
471 if(shareTracks(iter_i,iter_j)) continue;
472
473 if( fabs(0.5*(mass_i+mass_j)-mDcand) < deltaM){
474 deltaM = fabs(0.5*(mass_i+mass_j)-mDcand);
475 index_i = i;
476 index_j = j;
477 isDtcand = true;
478 }
479
480 } //end of j loop
481 } //end of i loop
482
483 if(isDtcand){
484 m_iterdtag1=m_iterbegin+igood1[index_i];
485 m_iterdtag2=m_iterbegin+igood2[index_j];
486 }
487
488 return isDtcand;
489
490} //end of findDtag
491
492bool DTagTool::findADTag(int mode1, int mode2){
493
494 int tagcharm1= (mode1<10 || mode1>=200)?+1:0;
495 int tagcharm2= (mode2<10 || mode2>=200)?-1:0;
496
497 if(tagcharm1*tagcharm2>0){
498 cout<<"double tagging warning! two modes can't have same nonzero charmness"<<endl;
499 return false;
500 }
501
502 //define D candidate mass
503 double mDcand=0;
504 if( mode1 < 200 && mode2 < 200)
505 mDcand = 1.8645;
506 else if ( mode1>=200 && mode1 < 400 && mode2>=200 && mode2 < 400)
507 mDcand = 1.8693;
508 else if ( mode1>=400 && mode1 < 1000 && mode2>=400 && mode2 < 1000)
509 mDcand = 1.9682;
510 else{
511 cout<<"double tag modes are not from same particle ! "<<endl;
512 return false;
513 }
514
515
516 vector<int> igood1, igood2;
517 igood1.clear(),igood2.clear();
518 int index=0;
519 EvtRecDTagCol::iterator iter_dtag=m_iterbegin;
520
521 //charge conjucation considered
522 for ( ; iter_dtag != m_iterend; iter_dtag++){
523 int iter_mode=(*iter_dtag)->decayMode();
524 int iter_charm=(*iter_dtag)->charm();
525 int iter_type=(*iter_dtag)->type();
526
527 if(m_pid){
528 if( iter_mode == mode1 && iter_charm == tagcharm1 && iter_type==1 )
529 igood1.push_back(index);
530 if( tagcharm1!=0 && iter_mode == mode1 && iter_charm == -tagcharm1 && iter_type==1 )
531 igood1.push_back(index);
532
533 if( iter_mode == mode2 && iter_charm == tagcharm2 && iter_type==1 )
534 igood2.push_back(index);
535 if( tagcharm2!=0 && iter_mode == mode2 && iter_charm == -tagcharm2 && iter_type==1 )
536 igood2.push_back(index);
537 }
538 else{
539 if( iter_mode == mode1 && iter_charm == tagcharm1 )
540 igood1.push_back(index);
541 if( tagcharm1!=0 && iter_mode == mode1 && iter_charm == -tagcharm1 )
542 igood1.push_back(index);
543
544 if( iter_mode == mode2 && iter_charm == tagcharm2 )
545 igood2.push_back(index);
546 if( tagcharm2!=0 && iter_mode == mode2 && iter_charm == -tagcharm2 )
547 igood2.push_back(index);
548 }
549
550 index++;
551 }
552
553 //look for the best pair of double-tagged event
554
555 bool isDtcand=false;
556 EvtRecDTagCol::iterator iter_i, iter_j;
557
558 for(int i=0; i<igood1.size(); i++){
559
560 iter_i=m_iterbegin+igood1[i];
561 int charm_i=(*iter_i)->charm();
562
563 for(int j=0;j<igood2.size();j++){
564 iter_j=m_iterbegin+igood2[j];
565 int charm_j=(*iter_j)->charm();
566 if( charm_i*charm_j>0 || igood2[j] == igood1[i] ) continue;
567 if(shareTracks(iter_i,iter_j)) continue;
568
569 m_viterdtag1.push_back(m_iterbegin+igood1[i]);
570 m_viterdtag2.push_back(m_iterbegin+igood2[j]);
571
572
573 } //end of j loop
574 } //end of i loop
575
576 if(m_viterdtag1.size()>0){
577 isDtcand=true;
578 }
579
580 return isDtcand;
581}
582
583bool DTagTool::findADTag(int mode1, int tagcharm1, int mode2, int tagcharm2){
584
585 if(tagcharm1*tagcharm2>0){
586 cout<<"double tagging warning! two modes can't have same nonzero charmness"<<endl;
587 return false;
588 }
589
590 //define D candidate mass
591 double mDcand=0;
592 if( mode1 < 200 && mode2 < 200)
593 mDcand = 1.8645;
594 else if (mode1>=200 && mode1 < 400 && mode2>=200 && mode2 < 400)
595 mDcand = 1.8693;
596 else if (mode1>=400 && mode1 < 1000 && mode2>=400 && mode2 < 1000)
597 mDcand = 1.9682;
598 else{
599 cout<<"double tag modes are not from same particle ! "<<endl;
600 return false;
601 }
602
603
604 vector<int> igood1, igood2;
605 igood1.clear(),igood2.clear();
606 int index=0;
607 EvtRecDTagCol::iterator iter_dtag=m_iterbegin;
608
609 for ( ; iter_dtag != m_iterend; iter_dtag++){
610 int iter_mode=(*iter_dtag)->decayMode();
611 int iter_charm=(*iter_dtag)->charm();
612 int iter_type=(*iter_dtag)->type();
613
614 if(m_pid){
615 if( iter_mode == mode1 && iter_charm == tagcharm1 && iter_type==1 )
616 igood1.push_back(index);
617
618 if( iter_mode == mode2 && iter_charm == tagcharm2 && iter_type==1)
619 igood2.push_back(index);
620 }
621 else{
622 if( iter_mode == mode1 && iter_charm == tagcharm1)
623 igood1.push_back(index);
624
625 if( iter_mode == mode2 && iter_charm == tagcharm2)
626 igood2.push_back(index);
627 }
628
629
630 index++;
631 }
632
633 //look for the best pair of double-tagged event
634
635 bool isDtcand=false;
636 double deltaM=1.00;
637 int index_i=0, index_j=0;
638 EvtRecDTagCol::iterator iter_i, iter_j;
639
640 for(int i=0; i<igood1.size(); i++){
641
642 iter_i=m_iterbegin+igood1[i];
643 int charm_i=(*iter_i)->charm();
644 for(int j=0;j<igood2.size();j++){
645 iter_j=m_iterbegin+igood2[j];
646 int charm_j=(*iter_j)->charm();
647 if( charm_i*charm_j>0 || igood2[j] == igood1[i] ) continue;
648
649 if(shareTracks(iter_i,iter_j)) continue;
650
651 m_viterdtag1.push_back(m_iterbegin+igood1[i]);
652 m_viterdtag2.push_back(m_iterbegin+igood2[j]);
653
654 } //end of j loop
655 } //end of i loop
656
657 if(m_viterdtag1.size()>0)
658 isDtcand=true;
659
660
661 return isDtcand;
662
663} //end of findADtag
664
665
666
667//////////////////////////////////
668//////////////////////////////////
669
670
671
672void DTagTool::operator<< ( EvtRecDTagCol::iterator iter){
673
674 cout<<" print mode:"<< (*iter)->decayMode()<<endl;
675 cout<<"beam energy is:"<< (*iter)->beamE()<<endl;
676 cout<<"mBC is:"<< (*iter)->mBC()<<endl;
677 cout<<"deltaE is:"<< (*iter)->deltaE()<<endl;
678 cout<<"inv mass is:"<< (*iter)->mass()<<endl;
679 cout<<"charge is:"<< (*iter)->charge()<<endl;
680 cout<<"charm is:"<< (*iter)->charm()<<endl;
681 cout<<"num of children is:"<< (*iter)->numOfChildren()<<endl;
682
683 cout<<"found "<< (*iter)->tracks().size()<<" same side tracks."<<endl;
684 cout<<"found "<< (*iter)->otherTracks().size()<<" other side tracks."<<endl;
685 cout<<"found "<< (*iter)->showers().size()<<" same side showers."<<endl;
686 cout<<"found "<< (*iter)->otherShowers().size()<<" other side showers."<<endl;
687
688}
689
690HepLorentzVector DTagTool::pi0p4(EvtRecPi0Col::iterator pi0Itr, bool isconstrain){
691
692
693 if(isconstrain){
694 HepLorentzVector p41= (*pi0Itr)->hiPfit();
695 HepLorentzVector p42= (*pi0Itr)->loPfit();
696 return (p41+p42);
697 }
698 else {
699 EvtRecTrack* trk1 = const_cast<EvtRecTrack*>((*pi0Itr)->hiEnGamma());
700 EvtRecTrack* trk2 = const_cast<EvtRecTrack*>((*pi0Itr)->loEnGamma());
701
702 RecEmcShower* emctrk1 = (trk1)->emcShower();
703 RecEmcShower* emctrk2 = (trk2)->emcShower();
704
705 HepLorentzVector p41=p4shower(emctrk1);
706 HepLorentzVector p42=p4shower(emctrk2);
707 return (p41+p42);
708 }
709
710}
711
712
713
714vector<int> DTagTool::pi0Id(EvtRecDTagCol::iterator iter_dtag, int numpi0){
715
716 SmartRefVector<EvtRecTrack> showers=(*iter_dtag)->showers();
717 if(showers.size()<2*numpi0){
718 cout<<"too many pi0 required"<<endl;
719 exit(1);
720 }
721
722 vector<int> canid;
723 canid.clear();
724
725 for(EvtRecPi0Col::iterator pi0Itr = m_pi0iterbegin;
726 pi0Itr < m_pi0iterend; pi0Itr++){
727
728 /// Access pi0 children
729 EvtRecTrack* heGammaTrk = const_cast<EvtRecTrack*>((*pi0Itr)->hiEnGamma());
730 EvtRecTrack* leGammaTrk = const_cast<EvtRecTrack*>((*pi0Itr)->loEnGamma());
731
732 int heGammaTrkId = heGammaTrk->trackId();
733 int leGammaTrkId = leGammaTrk->trackId();
734
735 for(int index=0; index<numpi0; ++index){
736 bool isheid=false;
737 bool isleid=false;
738
739 for(int i=index*2; i<index*2+2; ++i){
740
741 if(!isheid && heGammaTrkId == showers[i]->trackId())
742 isheid=true;
743 if(!isleid && leGammaTrkId == showers[i]->trackId())
744 isleid=true;
745 }
746
747 if(isheid && isleid)
748 canid.push_back( pi0Itr - m_pi0iterbegin);
749 }
750
751
752 if(canid.size()==numpi0){
753 return canid;
754 break;
755 }
756
757 } // End "pi0Itr" FOR Loop
758
759 return canid;
760
761}
762
763vector<int> DTagTool::ksId(EvtRecDTagCol::iterator iter_dtag, int numks){
764
765 SmartRefVector<EvtRecTrack> tracks=(*iter_dtag)->tracks();
766 if(tracks.size()<2*numks){
767 cout<<"too many kshort required"<<endl;
768 exit(1);
769 }
770 vector<int> canid;
771 canid.clear();
772
773 if(tracks.size()==0)
774 return canid;
775
776
777 for(EvtRecVeeVertexCol::iterator ksItr = m_ksiterbegin;
778 ksItr < m_ksiterend; ksItr++){
779
780 /// Needed to reject Lambda (and conversion?) combinations
781 if((*ksItr)->vertexId() != 310) continue;
782
783 EvtRecTrack* aKsChild1Trk = (*ksItr)->daughter(0);
784 EvtRecTrack* aKsChild2Trk = (*ksItr)->daughter(1);
785
786 int ksChild1TrkId = aKsChild1Trk->trackId();
787 int ksChild2TrkId = aKsChild2Trk->trackId();
788
789 for(int index=0; index<numks; ++index){
790 bool isc1id=false;
791 bool isc2id=false;
792
793 for(int i=index*2; i<index*2+2; ++i){
794 if(!isc1id && ksChild1TrkId == tracks[i]->trackId())
795 isc1id=true;
796 if(!isc2id && ksChild2TrkId == tracks[i]->trackId())
797 isc2id=true;
798 }
799
800 if(isc1id && isc2id)
801 canid.push_back( ksItr - m_ksiterbegin);
802 }
803
804 if(canid.size()==numks){
805 return canid;
806 break;
807 }
808 } // End "ksItr" FOR Loop
809
810 return canid;
811}
812
813HepLorentzVector DTagTool::p4shower(RecEmcShower* shower ){
814
815 double Eemc=shower->energy();
816 double phi=shower->phi();
817 double theta=shower->theta();
818 HepLorentzVector p4(Eemc*sin(theta)*cos(phi),Eemc*sin(theta)*sin(phi),Eemc*cos(theta),Eemc);
819 return p4;
820}
821
822HepLorentzVector DTagTool::p4(RecMdcKalTrack* mdcKalTrack, int pid){
823
824 HepVector zhelix;
825 double mass=0;
826
827 if(pid==0){
828 zhelix=mdcKalTrack->getZHelixE();
829 mass=0.000511;
830 }
831 else if(pid==1){
832 zhelix=mdcKalTrack->getZHelixMu();
833 mass= 0.105658;
834 }
835 else if(pid==2){
836 zhelix=mdcKalTrack->getZHelix();
837 mass=0.139570;
838 }
839 else if(pid==3){
840 zhelix=mdcKalTrack->getZHelixK();
841 mass= 0.493677;
842 }
843 else{
844 zhelix=mdcKalTrack->getZHelixP();
845 mass= 0.938272;
846 }
847
848 double dr(0),phi0(0),kappa(0),dz(0),tanl(0);
849 dr=zhelix[0];
850 phi0=zhelix[1];
851 kappa=zhelix[2];
852 dz=zhelix[3];
853 tanl=zhelix[4];
854
855 int charge=0;
856
857 if (kappa > 0.0000000001)
858 charge = 1;
859 else if (kappa < -0.0000000001)
860 charge = -1;
861
862 double pxy=0;
863 if(kappa!=0) pxy = 1.0/fabs(kappa);
864
865 double px = pxy * (-sin(phi0));
866 double py = pxy * cos(phi0);
867 double pz = pxy * tanl;
868
869 double e = sqrt( pxy*pxy + pz*pz + mass*mass );
870
871 return HepLorentzVector(px, py, pz, e);
872
873
874}
875
877
878 SmartRefVector<EvtRecTrack> pionid=(*m_iterbegin)->pionId();
879
880 for(int i=0; i < pionid.size() ;i++){
881 if( trk->trackId() == pionid[i]->trackId()){
882 return true;
883 break;
884 }
885 }
886
887 return false;
888}
889
890
892 SmartRefVector<EvtRecTrack> kaonid=(*m_iterbegin)->kaonId();
893
894 for(int i=0; i < kaonid.size() ;i++){
895 if( trk->trackId() == kaonid[i]->trackId()){
896 return true;
897 break;
898 }
899 }
900
901 return false;
902}
903
904
905
907
908 m_simplePIDSvc->preparePID(trk);
909 return ( m_simplePIDSvc->iselectron(true));
910
911 /*
912
913 double dedxchi=-99;
914 double Eemc=0;
915 double ptrk=-99;
916
917 if(trk->isMdcDedxValid()){
918 RecMdcDedx* dedxTrk=trk->mdcDedx();
919 dedxchi=dedxTrk->chiE();
920 }
921
922 if( trk->isMdcKalTrackValid() ){
923 RecMdcKalTrack *mdcKalTrk = trk->mdcKalTrack();
924 ptrk= mdcKalTrk->p();
925 }
926 if( trk->isEmcShowerValid()){
927 RecEmcShower *emcTrk = trk->emcShower();
928 Eemc=emcTrk->energy();
929 }
930
931 double eop = Eemc/ptrk;
932
933 if( fabs(eop)>0.8 && fabs(dedxchi)<5)
934 return true;
935 else
936 return false;
937 */
938}
939
941
942 double depth=-99;
943
944 double dedxchi=-99;
945 double Eemc=0;
946 double ptrk=-99;
947
948 if(trk->isMdcDedxValid()){
949 RecMdcDedx* dedxTrk=trk->mdcDedx();
950 dedxchi=dedxTrk->chiMu();
951 }
952
953 if( trk->isMdcKalTrackValid() ){
954 RecMdcKalTrack *mdcKalTrk = trk->mdcKalTrack();
955 ptrk= mdcKalTrk->p();
956 }
957 if( trk->isEmcShowerValid()){
958 RecEmcShower *emcTrk = trk->emcShower();
959 Eemc=emcTrk->energy();
960 }
961
962 double eop = Eemc/ptrk;
963
964 if( trk->isMucTrackValid() ){
965 RecMucTrack* mucTrk=trk->mucTrack();
966 depth=mucTrk->depth();
967 }
968
969 if( fabs(dedxchi)<5 && fabs(Eemc)>0.15 && fabs(Eemc)<0.3 && (depth>=80*ptrk-60 || depth>40))
970 return true;
971 return false;
972}
973
975
976 if( !trk->isMdcKalTrackValid()) {
977 return false;
978 }
979
980 Hep3Vector xorigin(0,0,0);
981 IVertexDbSvc* vtxsvc;
982 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
983 if(vtxsvc->isVertexValid()){
984 double* dbv = vtxsvc->PrimaryVertex();
985 double* vv = vtxsvc->SigmaPrimaryVertex();
986 xorigin.setX(dbv[0]);
987 xorigin.setY(dbv[1]);
988 xorigin.setZ(dbv[2]);
989 }
990
991
992 RecMdcKalTrack *mdcKalTrk = trk->mdcKalTrack();
994 HepVector a = mdcKalTrk->getZHelix();
995 HepSymMatrix Ea = mdcKalTrk->getZError();
996 HepPoint3D pivot(0.,0.,0.);
997 HepPoint3D IP(xorigin[0],xorigin[1],xorigin[2]);
998 VFHelix helixp(pivot,a,Ea);
999 helixp.pivot(IP);
1000 HepVector vec = helixp.a();
1001 double vrl = vec[0];
1002 double vzl = vec[3];
1003 double costheta = cos(mdcKalTrk->theta());
1004
1005 if(fabs(vrl)<1 && fabs(vzl)<10 && fabs(costheta)<0.93)
1006 return true;
1007 return false;
1008}
1009
1011 RecEmcShower *emcTrk = trk->emcShower();
1012 double eraw = emcTrk->energy();
1013 double time = emcTrk->time();
1014 double costh = cos(emcTrk->theta());
1015 if( (
1016 (fabs(costh)<0.80 && eraw>0.025)
1017 || (fabs(costh)>0.84 && eraw>0.05)
1018 ) && time>0 && time<14 )
1019 return true;
1020
1021 return false;
1022}
1023
1025
1026 //good track list
1027 vector<EvtRecTrackIterator> iGood;
1028 iGood.clear();
1029 for(EvtRecTrackIterator iter=m_chargebegin; iter!= m_chargeend; iter++){
1030 if(isGoodTrack(*iter))
1031 iGood.push_back(iter);
1032 }
1033
1034 if(iGood.size() != 2)
1035 return true;
1036
1037 //cosmic veto
1038 double time1=-99,time2=-99;
1039 for(vector<EvtRecTrackIterator>::size_type i=0;i<2;i++){
1040 if( (*iGood[i])->isTofTrackValid() ){
1041 SmartRefVector<RecTofTrack> tofTrkCol= (*iGood[i])->tofTrack();
1042 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
1043
1044 for(;iter_tof!=tofTrkCol.end();iter_tof++){
1045 TofHitStatus* status =new TofHitStatus;
1046 status->setStatus( (*iter_tof)->status() );
1047 if(status->is_cluster()){
1048 if(i==0)
1049 time1=(*iter_tof)->tof();
1050 else
1051 time2=(*iter_tof)->tof();
1052 }
1053 delete status;
1054 }
1055 }
1056 }
1057 if( time1!=-99 && time2!=-99 && fabs(time1-time2)>5)
1058 return false;
1059
1060 //rad bhabha veto
1061 if(isElectron( *iGood[0]) && isElectron(*iGood[1]))
1062 return false;
1063
1064 //rad dimu veto
1065 if(isMuon( *iGood[0]) && isMuon(*iGood[1]))
1066 return false;
1067
1068 return true;
1069}
1070
1071bool DTagTool::shareTracks(EvtRecDTagCol::iterator iter1,EvtRecDTagCol::iterator iter2){
1072
1073 SmartRefVector<EvtRecTrack> tracks1=(*iter1)->tracks();
1074 SmartRefVector<EvtRecTrack> showers1=(*iter1)->showers();
1075 SmartRefVector<EvtRecTrack> tracks2=(*iter2)->tracks();
1076 SmartRefVector<EvtRecTrack> showers2=(*iter2)->showers();
1077
1078 //charged tracks
1079 for(int i=0; i<tracks1.size(); i++){
1080 for(int j=0; j<tracks2.size(); j++){
1081 if(tracks1[i]==tracks2[j])
1082 return true;
1083 }
1084 }
1085
1086 //neutral showers
1087 for(int i=0; i<showers1.size(); i++){
1088 for(int j=0; j<showers2.size(); j++){
1089 if(showers1[i]==showers2[j])
1090 return true;
1091 }
1092 }
1093
1094 return false;
1095}
1096
1097IDataProviderSvc* DTagTool::eventSvc(){
1098
1099 if(m_evtSvc == 0){
1100
1101 StatusCode sc = Gaudi::svcLocator()->service ( "EventDataSvc", m_evtSvc, true);
1102 if( sc.isFailure() ) {
1103 assert(false);
1104 }
1105
1106 }
1107
1108 return m_evtSvc;
1109
1110}
double sin(const BesAngle a)
Definition: BesAngle.h:210
double cos(const BesAngle a)
Definition: BesAngle.h:213
double mass
Double_t time
EvtRecTrackCol::iterator EvtRecTrackIterator
Definition: EvtRecTrack.h:111
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
vector< int > pi0Id(EvtRecDTagCol::iterator iter, int numpi0=1)
Definition: DTagTool.cxx:714
bool isGoodTrack(EvtRecTrack *trk)
Definition: DTagTool.cxx:974
bool isPion(EvtRecTrack *trk)
Definition: DTagTool.cxx:876
void clear()
Definition: DTagTool.cxx:118
bool isElectron(EvtRecTrack *trk)
Definition: DTagTool.cxx:906
bool isKaon(EvtRecTrack *trk)
Definition: DTagTool.cxx:891
bool findDTag(EvtRecDTag::DecayMode mode1, EvtRecDTag::DecayMode mode2, string smass="mbc")
Definition: DTagTool.cxx:265
bool shareTracks(EvtRecDTagCol::iterator iter1, EvtRecDTagCol::iterator iter2)
Definition: DTagTool.cxx:1071
vector< int > mode(EvtRecDTag::DecayMode decaymode)
Definition: DTagTool.cxx:135
bool findSTag(EvtRecDTag::DecayMode mode, int tagcharm)
Definition: DTagTool.cxx:181
HepLorentzVector pi0p4(EvtRecPi0Col::iterator pi0Itr, bool isconstrain=true)
Definition: DTagTool.cxx:690
vector< int > ksId(EvtRecDTagCol::iterator iter, int numks=1)
Definition: DTagTool.cxx:763
void operator<<(EvtRecDTagCol::iterator iter)
Definition: DTagTool.cxx:672
IDataProviderSvc * eventSvc()
Definition: DTagTool.cxx:1097
bool isGoodShower(EvtRecTrack *trk)
Definition: DTagTool.cxx:1010
HepLorentzVector p4(RecMdcKalTrack *mdcKalTrack, int pid)
Definition: DTagTool.cxx:822
bool cosmicandleptonVeto()
Definition: DTagTool.cxx:1024
bool isMuon(EvtRecTrack *trk)
Definition: DTagTool.cxx:940
HepLorentzVector p4shower(RecEmcShower *shower)
Definition: DTagTool.cxx:813
bool findADTag(EvtRecDTag::DecayMode mode1, EvtRecDTag::DecayMode mode2)
Definition: DTagTool.cxx:279
DTagTool()
Definition: DTagTool.cxx:32
double theta() const
Definition: DstEmcShower.h:38
double phi() const
Definition: DstEmcShower.h:39
double time() const
Definition: DstEmcShower.h:50
double energy() const
Definition: DstEmcShower.h:45
double chiMu() const
Definition: DstMdcDedx.h:60
const double theta() const
static void setPidType(PidType pidType)
const double p() const
double depth() const
Definition: DstMucTrack.h:45
bool isMdcDedxValid()
Definition: EvtRecTrack.h:45
int trackId() const
Definition: EvtRecTrack.h:32
RecMucTrack * mucTrack()
Definition: EvtRecTrack.h:59
RecMdcDedx * mdcDedx()
Definition: EvtRecTrack.h:55
bool isMdcKalTrackValid()
Definition: EvtRecTrack.h:44
bool isMucTrackValid()
Definition: EvtRecTrack.h:48
RecEmcShower * emcShower()
Definition: EvtRecTrack.h:58
bool isEmcShowerValid()
Definition: EvtRecTrack.h:47
RecMdcKalTrack * mdcKalTrack()
Definition: EvtRecTrack.h:54
virtual void preparePID(EvtRecTrack *track)=0
virtual bool iselectron(bool eop=false)=0
virtual bool isVertexValid()=0
virtual double * SigmaPrimaryVertex()=0
virtual double * PrimaryVertex()=0
const HepVector & getZHelix() const
HepVector & getZHelixP()
HepVector & getZHelixK()
HepVector & getZHelixE()
HepVector & getZHelixMu()
const HepSymMatrix & getZError() const
bool is_cluster() const
Definition: TofHitStatus.h:25
void setStatus(unsigned int status)
const HepPoint3D & pivot(void) const
returns pivot position.
const HepVector & a(void) const
returns helix parameters.
_EXTERN_ std::string EvtRecDTagCol
Definition: EventModel.h:117
dble_vec_t vec[12]
Definition: ranlxd.c:372