287 {
288
289
290
291 MsgStream log(
msgSvc(), name());
292 log << MSG::INFO << "in execute()" << endreq;
293 m_cout_all ++;
294 StatusCode sc=StatusCode::SUCCESS;
295
296 setFilterPassed(false);
297
298 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
299 if(!eventHeader){
300 log << MSG::ERROR << "EventHeader not found" << endreq;
301 return StatusCode::SUCCESS;
302 }
303 int run(eventHeader->runNumber());
304 int event(eventHeader->eventNumber());
305 if(event%1000==0) cout << "run: " << run << " event: " << event << endl;
306
307 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), "/Event/EvtRec/EvtRecEvent");
308 if(!evtRecEvent){
309 log << MSG::ERROR << "EvtRecEvent not found" << endreq;
310 return StatusCode::SUCCESS;
311 }
312 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
313 << evtRecEvent->totalCharged() << " , "
314 << evtRecEvent->totalNeutral() << " , "
315 << evtRecEvent->totalTracks() <<endreq;
316
317 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), "/Event/EvtRec/EvtRecTrackCol");
318 if(!evtRecTrkCol){
319 log << MSG::ERROR << "EvtRecTrackCol" << endreq;
320 return StatusCode::SUCCESS;
321 }
322
323 if(m_trigger_flag){
325 if (!trigData) {
326 log << MSG::FATAL << "Could not find Trigger Data for physics analysis" << endreq;
327 return StatusCode::FAILURE;
328 }
329
330 log << MSG::DEBUG << "Trigger conditions: " << endreq;
331 for(int i=0; i < 48; i++){
332 log << MSG::DEBUG << "i:" << i << " name:" << trigData->getTrigCondName(i) << " cond:" << trigData->getTrigCondition(i) << endreq;
333 }
334
335 int m_trig_tot(0), m_trig_which(-1);
336 if(m_eventrate){
337 for(int j=0; j<16; j++){
338 if(trigData->getTrigChannel(j)){
339 m_trig_tot ++;
340 m_trig_which = j+1;
341 }
342 }
343 if(m_trig_tot==1 && m_trig_which==m_chan_det)
m_cout_everat++;
344 return sc;
345 }
346 }
347
349 if(evtRecEvent->totalCharged()<3 || evtRecTrkCol->size()<3 || evtRecEvent->totalTracks()>99 || evtRecTrkCol->size()>99) return StatusCode::SUCCESS;
351
352
353 Vint iGood; iGood.clear();
354 int m_num[4]={0,0,0,0};
355 int nCharge = 0;
356 m_pion_matched = 0; m_lep_matched = 0;
357 HepLorentzVector m_lv_pionp, m_lv_pionm, m_lv_lepp, m_lv_lepm, m_lv_ele, m_lv_pos, m_lv_mum, m_lv_mup;
358
359 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
361 if(!(*itTrk)->isMdcKalTrackValid()) continue;
363
368 if(fabs(m_vz0) >= m_vz0cut) continue;
369 if(m_vr0 >= m_vr0cut) continue;
370 iGood.push_back(i);
371 nCharge += mdcTrk->
charge();
372 if(mdcTrk->
p()<1.0){
if((*itTrk)->isEmcShowerValid()) m_pion_matched ++; }
373 else{ if((*itTrk)->isEmcShowerValid()) m_lep_matched ++; }
374
378
379 m_lv_pionp = mdcTrk->
p4(
xmass[2]);
380 m_num[0] ++;
381 }
382 else{
384 m_lv_pos = mdcTrk->
p4(
xmass[0]);
386 m_lv_mup = mdcTrk->
p4(
xmass[1]);
387 m_num[2] ++;
388 }
389 }
390 else{
393 m_lv_pionm = mdcTrk->
p4(
xmass[2]);
394 m_num[1] ++;
395 }
396 else{
398 m_lv_ele = mdcTrk->
p4(
xmass[0]);
400 m_lv_mum = mdcTrk->
p4(
xmass[1]);
401 m_num[3] ++;
402 }
403 }
404 }
405
406 int nGood = iGood.size();
407 log << MSG::DEBUG << "With KalmanTrack, ngood, totcharge = " << nGood << " , " << nCharge << endreq;
408 if(nGood<3 || nGood>4) return sc;
410
411 m_ep_ratio = 0;
412 for(int i=0; i< evtRecEvent->totalTracks(); i++){
414 if(!(*itTrk)->isEmcShowerValid()) continue;
416 m_ep_ratio += emcTrk->
energy();
417 }
418
419 if(m_ep_ratio > m_distin_emuon){
420 m_lv_lepp = m_lv_pos;
421 m_lv_lepm = m_lv_ele;
422 }
423 else{
424 m_lv_lepp = m_lv_mup;
425 m_lv_lepm = m_lv_mum;
426 }
427
428 HepLorentzVector m_lv_lab(0.04,0,0,3.686);
429 if(nGood==4){
430 if(nCharge) return sc;
431 m_event_flag = 4;
432 }
433 else{
434 if(m_num[0]>1 || m_num[1]>1 || m_num[2]>1 || m_num[3]>1) return sc;
435 if(m_num[0]==0){
436 if(nCharge != -1) return StatusCode::SUCCESS;
437 m_lv_pionp = m_lv_lab - m_lv_pionm - m_lv_lepp - m_lv_lepm;
438 if(m_lv_pionp.vect().cosTheta()>m_cosThetaCut)return StatusCode::SUCCESS;
439 m_event_flag = 0;
440 }
441 if(m_num[1]==0){
442 if(nCharge != 1) return StatusCode::SUCCESS;
443 m_lv_pionm = m_lv_lab - m_lv_pionp - m_lv_lepp - m_lv_lepm;
444 if(m_lv_pionm.vect().cosTheta()>m_cosThetaCut)return StatusCode::SUCCESS;
445 m_event_flag = 1;
446 }
447 if(m_num[2]==0){
448 if(nCharge != -1) return StatusCode::SUCCESS;
449 m_lv_lepp = m_lv_lab - m_lv_pionp - m_lv_pionm - m_lv_lepm;
450 if(m_lv_lepp.vect().cosTheta()>m_cosThetaCut)return StatusCode::SUCCESS;
451 m_event_flag = 2;
452 }
453 if(m_num[3]==0){
454 if(nCharge != 1) return StatusCode::SUCCESS;
455 m_lv_lepm = m_lv_lab - m_lv_pionp - m_lv_pionm - m_lv_lepp;
456 if(m_lv_lepm.vect().cosTheta()>m_cosThetaCut)return StatusCode::SUCCESS;
457 m_event_flag = 3;
458 }
459 }
461
462
463
464 HepLorentzVector m_lv_recoil, m_lv_jpsi;
465 m_lv_recoil = m_lv_lab - m_lv_pionp - m_lv_pionm;
466 m_lv_jpsi = m_lv_lepp + m_lv_lepm;
467
468 m_mass_twopi = (m_lv_pionp + m_lv_pionm).m();
469 m_mass_recoil = m_lv_recoil.m();
470 m_mass_jpsi = m_lv_jpsi.m();
471
472
473 if( m_mass_recoil < 3.05 || m_mass_recoil > 3.15 ) return sc;
474 if( m_mass_jpsi < 3.0 || m_mass_jpsi > 3.2 ) return sc;
476
477 HepLorentzVector m_ttm(m_lv_jpsi + m_lv_pionp + m_lv_pionm);
478 if(m_ttm.m()>4 || m_ttm.m()<3) return sc;
479
480
481 m_pipi_dang = m_lv_pionp.vect().cosTheta(m_lv_pionm.vect());
482
483 m_mom_pionp = m_lv_pionp.vect().mag();
484 m_mom_pionm = m_lv_pionm.vect().mag();
485 m_mom_lepp = m_lv_lepp.vect().mag();
486 m_mom_lepm = m_lv_lepm.vect().mag();
487 m_trans_ratio_lepp = m_lv_lepp.vect().perp()/m_lv_lepp.vect().mag();
488 m_trans_ratio_lepm = m_lv_lepm.vect().perp()/m_lv_lepm.vect().mag();
489 m_trans_ratio_pionp = m_lv_pionp.vect().perp()/m_lv_pionp.vect().mag();
490 m_trans_ratio_pionm = m_lv_pionm.vect().perp()/m_lv_pionm.vect().mag();
491
492 Hep3Vector m_boost_jpsi(m_lv_recoil.boostVector());
493 HepLorentzVector m_lv_cms_lepp(boostOf(m_lv_lepp,-m_boost_jpsi));
494 HepLorentzVector m_lv_cms_lepm(boostOf(m_lv_lepm,-m_boost_jpsi));
495 m_cms_lepm = m_lv_cms_lepm.vect().mag();
496 m_cms_lepp = m_lv_cms_lepp.vect().mag();
497 log << MSG::DEBUG << "jpsi four momentum in cms " << m_lv_cms_lepp + m_lv_cms_lepm <<endreq;
498
499 m_inv_mass = m_ttm.m();
500 m_tot_e = m_ttm.e();
501 m_tot_px = m_ttm.px();
502 m_tot_py = m_ttm.py();
503 m_tot_pz = m_ttm.pz();
504 m_run = run;
505 m_event = event;
506 HepLorentzVector m_lv_book(0,0,0,0);
507 for(m_index=0; m_index<4; m_index++){
508 switch(m_index){
509 case 0: m_lv_book = m_lv_pionp;
510 break;
511 case 1: m_lv_book = m_lv_pionm;
512 break;
513 case 2: m_lv_book = m_lv_lepp;
514 break;
515 case 3: m_lv_book = m_lv_lepm;
516 break;
517 default: m_lv_book.setE(2008);
518 }
519 m_cos_theta[m_index] = m_lv_book.vect().cosTheta();
520 m_phi[m_index] = m_lv_book.vect().phi();
521 m_four_mom[m_index][0] = m_lv_book.e();
522 m_four_mom[m_index][1] = m_lv_book.px();
523 m_four_mom[m_index][2] = m_lv_book.py();
524 m_four_mom[m_index][3] = m_lv_book.pz();
525 }
526
527 if(m_subsample_flag) setFilterPassed(true);
528 else if(m_mass_recoil>3.08 && m_mass_recoil<3.12 && m_mass_jpsi>3.0 && m_mass_jpsi<3.2 && m_cms_lepp<1.7 && m_cms_lepp>1.4 && m_cms_lepm<1.7 && m_cms_lepm>1.4 && m_event_flag==4 && m_pipi_dang<m_pipi_dang_cut) setFilterPassed(true);
529
530
531
532 SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(), "/Event/MC/McParticleCol");
533 if(m_run<0){
534 int m_numParticle(0), m_true_pid(0);
535 if(!mcParticleCol){
536 log << MSG::ERROR << "Could not retrieve McParticelCol" << endreq;
537 return StatusCode::FAILURE;
538 }
539 else{
540 bool psipDecay(false);
541 int rootIndex(-1);
542 Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
543 for (; iter_mc != mcParticleCol->end(); iter_mc++){
544 if ((*iter_mc)->primaryParticle()) continue;
545 if (!(*iter_mc)->decayFromGenerator()) continue;
546
547 if ((*iter_mc)->particleProperty()==100443){
548 psipDecay = true;
549 rootIndex = (*iter_mc)->trackIndex();
550 }
551 if (!psipDecay) continue;
552 int mcidx = ((*iter_mc)->mother()).trackIndex() - rootIndex;
553 int pdgid = (*iter_mc)->particleProperty();
554 m_pdgid[m_numParticle] = pdgid;
555 m_motheridx[m_numParticle] = mcidx;
556 m_numParticle ++;
557
558
559 if((*iter_mc)->particleProperty() == 211) m_true_pionp = (*iter_mc)->initialFourMomentum().vect().mag();
560 if((*iter_mc)->particleProperty() == -211) m_true_pionm = (*iter_mc)->initialFourMomentum().vect().mag();
561 }
562 m_idxmc = m_numParticle;
563 }
564 }
565
566 m_tuple1->write();
567 m_tuple8->write();
568
569
570
571 Vint iGam; iGam.clear();
572 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
574 if(!(*itTrk)->isEmcShowerValid()) continue;
576 Hep3Vector emcpos(emcTrk->
x(), emcTrk->
y(), emcTrk->
z());
577
578 double dthe = 200.;
579 double dphi = 200.;
580 double dang = 200.;
581 for(int j = 0; j < evtRecEvent->totalCharged(); j++) {
583 if(!(*jtTrk)->isExtTrackValid()) continue;
587
588 double angd = extpos.angle(emcpos);
589 double thed = extpos.theta() - emcpos.theta();
590 double phid = extpos.deltaPhi(emcpos);
591 thed = fmod(thed+CLHEP::twopi+CLHEP::twopi+
pi, CLHEP::twopi) - CLHEP::pi;
592 phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+
pi, CLHEP::twopi) - CLHEP::pi;
593
594 if(fabs(thed) < fabs(dthe)) dthe = thed;
595 if(fabs(phid) < fabs(dphi)) dphi = phid;
596 if(angd < dang) dang = angd;
597 }
598 if(dang>=200) continue;
599 double eraw = emcTrk->
energy();
600 dthe = dthe * 180 / (CLHEP::pi);
601 dphi = dphi * 180 / (CLHEP::pi);
602 dang = dang * 180 / (CLHEP::pi);
603 m_dthe = dthe;
604 m_dphi = dphi;
605 m_dang = dang;
606 m_eraw = eraw;
607
608
609
610 iGam.push_back((*itTrk)->trackId());
611 }
612
613 m_nGam = iGam.size();
614 log << MSG::DEBUG << "num Good Photon " << m_nGam << " , " <<evtRecEvent->totalNeutral()<<endreq;
615 m_tuple2->write();
616
617
618
619
620 if(m_checkDedx) {
621 int m_dedx_cout(0);
622 for(int i = 0; i < nGood; i++) {
624 if(!(*itTrk)->isMdcDedxValid())continue;
627
628 m_ptrk = mdcTrk->
p();
629 m_chie = dedxTrk->
chiE();
630 m_chimu = dedxTrk->
chiMu();
631 m_chipi = dedxTrk->
chiPi();
632 m_chik = dedxTrk->
chiK();
633 m_chip = dedxTrk->
chiP();
636 m_probPH = dedxTrk->
probPH();
637 m_normPH = dedxTrk->
normPH();
638
639 m_tuple3->write();
640 }
641 }
642
643
644
645
646 if(m_checkTof) {
647 int m_endcap_cout(0), m_layer1_cout(0), m_layer2_cout(0);
648 for(int i = 0; i < nGood; i++) {
650 if(!(*itTrk)->isTofTrackValid()) continue;
651
653 SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
654
655 double ptrk = mdcTrk->
p();
656
657 for( SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin() ;iter_tof != tofTrkCol.end(); iter_tof++ ) {
659 status->
setStatus((*iter_tof)->status());
662 if( status->
layer()!=0 )
continue;
663 double path = (*iter_tof)->path();
664 double tof = (*iter_tof)->tof();
665 double ph = (*iter_tof)->ph();
666 double rhit = (*iter_tof)->zrhit();
667 double qual = 0.0 + (*iter_tof)->quality();
668 double cntr = 0.0 + (*iter_tof)->tofID();
669 double texp[5];
670 for(int j = 0; j < 5; j++) {
671 double gb =
xmass[j]/ptrk;
672 double beta = sqrt(1+gb*gb);
673 texp[j] = path*beta/
velc;
674 }
675 m_cntr_etof = cntr;
676 m_ptot_etof = ptrk;
677 m_path_etof = path;
678 m_ph_etof = ph;
679 m_rhit_etof = rhit;
680 m_qual_etof = qual;
681 m_tof_etof = tof;
682 m_te_etof = tof - texp[0];
683 m_tmu_etof = tof - texp[1];
684 m_tpi_etof = tof - texp[2];
685 m_tk_etof = tof - texp[3];
686 m_tp_etof = tof - texp[4];
687
688 m_tuple4->write();
689 }
690 else {
692 if(status->
layer()==1){
693 double path=(*iter_tof)->path();
694 double tof = (*iter_tof)->tof();
695 double ph = (*iter_tof)->ph();
696 double rhit = (*iter_tof)->zrhit();
697 double qual = 0.0 + (*iter_tof)->quality();
698 double cntr = 0.0 + (*iter_tof)->tofID();
699 double texp[5];
700 for(int j = 0; j < 5; j++) {
701 double gb =
xmass[j]/ptrk;
702 double beta = sqrt(1+gb*gb);
703 texp[j] = path*beta/
velc;
704 }
705
706 m_cntr_btof1 = cntr;
707 m_ptot_btof1 = ptrk;
708 m_path_btof1 = path;
709 m_ph_btof1 = ph;
710 m_zhit_btof1 = rhit;
711 m_qual_btof1 = qual;
712 m_tof_btof1 = tof;
713 m_te_btof1 = tof - texp[0];
714 m_tmu_btof1 = tof - texp[1];
715 m_tpi_btof1 = tof - texp[2];
716 m_tk_btof1 = tof - texp[3];
717 m_tp_btof1 = tof - texp[4];
718
719 m_tuple5->write();
720 }
721
722 if(status->
layer()==2){
723 double path=(*iter_tof)->path();
724 double tof = (*iter_tof)->tof();
725 double ph = (*iter_tof)->ph();
726 double rhit = (*iter_tof)->zrhit();
727 double qual = 0.0 + (*iter_tof)->quality();
728 double cntr = 0.0 + (*iter_tof)->tofID();
729 double texp[5];
730 for(int j = 0; j < 5; j++) {
731 double gb =
xmass[j]/ptrk;
732 double beta = sqrt(1+gb*gb);
733 texp[j] = path*beta/
velc;
734 }
735
736 m_cntr_btof2 = cntr;
737 m_ptot_btof2 = ptrk;
738 m_path_btof2 = path;
739 m_ph_btof2 = ph;
740 m_zhit_btof2 = rhit;
741 m_qual_btof2 = qual;
742 m_tof_btof2 = tof;
743 m_te_btof2 = tof - texp[0];
744 m_tmu_btof2 = tof - texp[1];
745 m_tpi_btof2 = tof - texp[2];
746 m_tk_btof2 = tof - texp[3];
747 m_tp_btof2 = tof - texp[4];
748
749 m_tuple6->write();
750 }
751 }
752
753 delete status;
754 }
755 }
756 }
757
758
759 return StatusCode::SUCCESS;
760}
EvtRecTrackCol::iterator EvtRecTrackIterator
static long m_cout_col(0)
static long m_cout_everat(0)
static long m_cout_charge(0)
static long m_cout_nGood(0)
static long m_cout_recoil(0)
static long m_cout_mom(0)
const Hep3Vector emcPosition() const
const int emcVolumeNumber() const
static void setPidType(PidType pidType)
const HepLorentzVector p4() const
unsigned int layer() const
void setStatus(unsigned int status)
_EXTERN_ std::string TrigData