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