Loading...
Searching...
No Matches
Ppjrhopi Class Reference
#include <Ppjrhopi.h>
Inheritance diagram for Ppjrhopi:
Public Member Functions | |
| Ppjrhopi (const std::string &name, ISvcLocator *pSvcLocator) | |
| StatusCode | initialize () |
| StatusCode | execute () |
| StatusCode | finalize () |
Detailed Description
Definition at line 8 of file Ppjrhopi.h.
Constructor & Destructor Documentation
◆ Ppjrhopi()
| Ppjrhopi::Ppjrhopi | ( | const std::string & | name, |
| ISvcLocator * | pSvcLocator ) |
Definition at line 40 of file Ppjrhopi.cxx.
41 : Algorithm( name, pSvcLocator ) {
42
43 // Declare the properties
44 declareProperty( "Vr0cut", m_vr0cut = 5.0 );
45 declareProperty( "Vz0cut", m_vz0cut = 20.0 );
46 declareProperty( "Vr1cut", m_vr1cut = 1.0 );
47 declareProperty( "Vz1cut", m_vz1cut = 5.0 );
48 declareProperty( "Vctcut", m_cthcut = 0.93 );
49 declareProperty( "EnergyThreshold", m_energyThreshold = 0.04 );
50 declareProperty( "GammaAngCut", m_gammaAngCut = 20.0 );
51 declareProperty( "Test4C", m_test4C = 1 );
52 declareProperty( "Test5C", m_test5C = 1 );
53 declareProperty( "CheckDedx", m_checkDedx = 1 );
54 declareProperty( "CheckTof", m_checkTof = 1 );
55}
Referenced by Ppjrhopi().
Member Function Documentation
◆ execute()
| StatusCode Ppjrhopi::execute | ( | ) |
Definition at line 279 of file Ppjrhopi.cxx.
279 {
280
281 // std::cout << "execute()" << std::endl;
282
283 MsgStream log( msgSvc(), name() );
284 log << MSG::INFO << "in execute()" << endmsg;
285
286 setFilterPassed( false );
287
288 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
290 int event = eventHeader->eventNumber();
292
293 Ncut0++;
294
295 // FOR 6.4.0 EventModel::Recon--->EventModel::EvtRec IN 6.4.1
296 SmartDataPtr<EvtRecEvent> evtRecEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
297 log << MSG::INFO << "get event tag OK" << endmsg;
298 log << MSG::DEBUG << "ncharg, nneu, tottks = " << evtRecEvent->totalCharged() << " , "
299 << evtRecEvent->totalNeutral() << " , " << evtRecEvent->totalTracks() << endmsg;
300
301 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
302 //
303 // check x0, y0, z0, r0
304 // suggest cut: |z0|<5 && r0<1
305 //
306 if ( evtRecEvent->totalNeutral() > 100 ) { return StatusCode::SUCCESS; }
307
308 Vint iGood, ipip, ipim;
309 iGood.clear();
310 ipip.clear();
311 ipim.clear();
312 Vp4 ppip, ppim;
313 ppip.clear();
314 ppim.clear();
315
316 Hep3Vector xorigin( 0, 0, 0 );
317
318 // if (m_reader.isRunNumberValid(runNo)) {
319 IVertexDbSvc* vtxsvc;
320 Gaudi::svcLocator()->service( "VertexDbSvc", vtxsvc ).ignore();
322 {
325 // HepVector dbv = m_reader.PrimaryVertex(runNo);
326 // HepVector vv = m_reader.SigmaPrimaryVertex(runNo);
327 xorigin.setX( dbv[0] );
328 xorigin.setY( dbv[1] );
329 xorigin.setZ( dbv[2] );
330 }
331
332 int nCharge = 0;
333 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
334 {
335 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
336 if ( !( *itTrk )->isMdcTrackValid() ) continue;
337 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
338 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
339
345 double xv = xorigin.x();
346 double yv = xorigin.y();
348 // 2009//4
349 double m_vx0 = x0;
350 double m_vy0 = y0;
351 double m_vz0 = z0 - xorigin.z();
352 double m_vr0 = Rxy;
353 double m_Vctc = z0 / sqrt( Rxy * Rxy + z0 * z0 );
355
356 // m_tuple1->write();
357 // 2009//4
358 if ( fabs( m_vz0 ) >= m_vz0cut ) continue;
359 if ( m_vr0 >= m_vr0cut ) continue;
360 // if(fabs(m_Vct)>=m_cthcut) continue;
361 iGood.push_back( ( *itTrk )->trackId() );
362 nCharge += mdcTrk->charge();
363 }
364
365 //
366 // Finish Good Charged Track Selection
367 //
368 int nGood = iGood.size();
369
370 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endmsg;
371 if ( ( nGood != 4 ) || ( nCharge != 0 ) ) { return StatusCode::SUCCESS; }
372
373 Ncut1++;
374
375 Vint iGam;
376 iGam.clear();
377 for ( int i = evtRecEvent->totalCharged(); i < evtRecEvent->totalTracks(); i++ )
378 {
379 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
380 if ( !( *itTrk )->isEmcShowerValid() ) continue;
381 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
383 // find the nearest charged track
384 double dthe = 200.;
385 double dphi = 200.;
386 double dang = 200.;
387 log << MSG::DEBUG << "liuf neu= " << i << endmsg;
388 for ( int j = 0; j < evtRecEvent->totalCharged(); j++ )
389 {
390 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + j;
391 if ( !( *jtTrk )->isExtTrackValid() ) continue;
392 RecExtTrack* extTrk = ( *jtTrk )->extTrack();
394 Hep3Vector extpos = extTrk->emcPosition();
395 log << MSG::DEBUG << "liuf charge= " << j << endmsg;
396 // double ctht = extpos.cosTheta(emcpos);
397 double angd = extpos.angle( emcpos );
398 double thed = extpos.theta() - emcpos.theta();
399 double phid = extpos.deltaPhi( emcpos );
400 thed = fmod( thed + CLHEP::twopi + CLHEP::twopi + pi, CLHEP::twopi ) - CLHEP::pi;
401 phid = fmod( phid + CLHEP::twopi + CLHEP::twopi + pi, CLHEP::twopi ) - CLHEP::pi;
402
403 // if(fabs(thed) < fabs(dthe)) dthe = thed;
404 // if(fabs(phid) < fabs(dphi)) dphi = phid;
405 if ( angd < dang )
406 {
407 dang = angd;
408 dthe = thed;
409 dphi = phid;
410 }
411 }
412 if ( dang >= 200 ) continue;
414 dthe = dthe * 180 / ( CLHEP::pi );
415 dphi = dphi * 180 / ( CLHEP::pi );
416 dang = dang * 180 / ( CLHEP::pi );
417 // 2009//4
418 double m_dthe = dthe;
419 double m_dphi = dphi;
420 double m_dang = dang;
421 double m_eraw = eraw;
422 // m_tuple2->write();
423 // 2009//4
424 log << MSG::DEBUG << "eraw dang= " << eraw << " , " << dang << "," << i << endmsg;
425 if ( eraw < m_energyThreshold ) continue;
426 if ( dang < m_gammaAngCut ) continue;
427 // if((fabs(dthe) < m_gammaThetaCut) && (fabs(dphi)<m_gammaPhiCut) ) continue;
428 //
429 // good photon cut will be set here
430 //
431 iGam.push_back( ( *itTrk )->trackId() );
432 }
433
434 //
435 // Finish Good Photon Selection
436 //
437 int nGam = iGam.size();
438
439 log << MSG::DEBUG << "num Good Photon " << nGam << " , " << evtRecEvent->totalNeutral()
440 << endmsg;
441 if ( nGam < 2 ) { return StatusCode::SUCCESS; }
442
443 Ncut2++;
444
445 //
446 // Assign 4-momentum to each photon
447 //
448
449 Vp4 pGam;
450 pGam.clear();
451 for ( int i = 0; i < nGam; i++ )
452 {
453 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGam[i];
454 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
458 HepLorentzVector ptrk;
461 ptrk.setPz( eraw * cos( the ) );
462 ptrk.setE( eraw );
463
464 // ptrk = ptrk.boost(-0.011,0,0);// boost to cms
465
466 pGam.push_back( ptrk );
467 }
468
469 for ( int i = 0; i < nGood; i++ )
470 { // for rhopi without PID
471 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
472 if ( !( *itTrk )->isMdcTrackValid() ) continue;
473 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
474 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
475 RecMdcKalTrack* mdcKalTrk = ( *itTrk )->mdcKalTrack();
478 {
479 ipip.push_back( iGood[i] );
480 HepLorentzVector ptrk;
481 ptrk.setPx( mdcKalTrk->px() );
482 ptrk.setPy( mdcKalTrk->py() );
483 ptrk.setPz( mdcKalTrk->pz() );
484 double p3 = ptrk.mag();
486 ppip.push_back( ptrk );
487 }
488 else
489 {
490 ipim.push_back( iGood[i] );
491 HepLorentzVector ptrk;
492 ptrk.setPx( mdcKalTrk->px() );
493 ptrk.setPy( mdcKalTrk->py() );
494 ptrk.setPz( mdcKalTrk->pz() );
495 double p3 = ptrk.mag();
497 ppim.push_back( ptrk );
498 }
499 } // without PID
500
501 int npip = ipip.size();
502 int npim = ipim.size();
503 if ( npip != 2 || npim != 2 ) return StatusCode::SUCCESS;
504 /*
505 log << MSG::DEBUG << "ngood track ID = " << ipip[0] << " , "
506 << ipim[0]<< " , " << ipip[1] << " , " << ipim[1] << endmsg;
507 */
508 Ncut3++;
509
510 //
511 // find the two pi from the primary vetex
512 // ipip[0] && ipim[0] from ppsi
513 // ipip[1] && ipim[1] from jpsi
514 // should change track ID
515 //
516 HepLorentzVector pTot0( 0.011 * 3.6862, 0, 0, 3.6862 );
517 HepLorentzVector pTrec1, pTrec2, pTrec3, pTrec4;
518 HepLorentzVector pTrecf;
519 double m_recjpsi1, m_recjpsi2, m_recjpsi3, m_recjpsi4, m_recppf;
520 double deljp1, deljp2, deljp3, deljp4;
521 pTrec1 = pTot0 - ppip[0] - ppim[0];
522 pTrec2 = pTot0 - ppip[0] - ppim[1];
523 pTrec3 = pTot0 - ppip[1] - ppim[0];
524 pTrec4 = pTot0 - ppip[1] - ppim[1];
525 m_recjpsi1 = pTrec1.m();
526 m_recjpsi2 = pTrec2.m();
527 m_recjpsi3 = pTrec3.m();
528 m_recjpsi4 = pTrec4.m();
529 deljp1 = fabs( m_recjpsi1 - 3.097 );
530 deljp2 = fabs( m_recjpsi2 - 3.097 );
531 deljp3 = fabs( m_recjpsi3 - 3.097 );
532 deljp4 = fabs( m_recjpsi4 - 3.097 );
533
534 int itmp, itmp1, itmp2;
535 HepLorentzVector ptmp, ptmp1, ptmp2;
536
537 pTrecf = pTrec1;
538 m_recppf = pTrec1.m();
539
540 if ( deljp2 < deljp1 && deljp2 < deljp3 && deljp2 < deljp4 )
541 {
542 itmp = ipim[1];
543 ipim[1] = ipim[0];
544 ipim[0] = itmp;
545
546 ptmp = ppim[1];
547 ppim[1] = ppim[0];
548 ppim[0] = ptmp;
549
550 pTrecf = pTrec2;
551 m_recppf = pTrec2.m();
552 }
553
554 if ( deljp3 < deljp1 && deljp3 < deljp2 && deljp3 < deljp4 )
555 {
556 itmp = ipip[1];
557 ipip[1] = ipip[0];
558 ipip[0] = itmp;
559
560 ptmp = ppip[1];
561 ppip[1] = ppip[0];
562 ppip[0] = ptmp;
563
564 pTrecf = pTrec3;
565 m_recppf = pTrec3.m();
566 }
567
568 if ( deljp4 < deljp1 && deljp4 < deljp2 && deljp4 < deljp3 )
569 {
570 itmp1 = ipip[1];
571 ipip[1] = ipip[0];
572 ipip[0] = itmp1;
573 itmp2 = ipim[1];
574 ipim[1] = ipim[0];
575 ipim[0] = itmp2;
576
577 ptmp1 = ppip[1];
578 ppip[1] = ppip[0];
579 ppip[0] = ptmp1;
580 ptmp2 = ppim[1];
581 ppim[1] = ppim[0];
582 ppim[0] = ptmp2;
583
584 pTrecf = pTrec4;
585 m_recppf = pTrec4.m();
586 }
587
588 if ( fabs( m_recppf - 3.097 ) > 0.2 ) return StatusCode::SUCCESS;
589
590 log << MSG::DEBUG << "ngood track ID after jpsi = " << ipip[0] << " , " << ipim[0] << " , "
591 << ipip[1] << " , " << ipim[1] << endmsg;
592 Ncut4++;
593
594 HepLorentzVector ppi2_no1 = ppip[0] + ppim[0];
595 HepLorentzVector ppi2_no2 = ppip[1] + ppim[1];
596 HepLorentzVector ppi2_no3 = ppip[0] + ppim[1];
597 HepLorentzVector ppi2_no4 = ppip[1] + ppim[0];
598 HepLorentzVector p4pi_no = ppi2_no1 + ppi2_no2;
599
600 double emcTg1 = 0.0;
601 double emcTg2 = 0.0;
602 double emcTg3 = 0.0;
603 double emcTg4 = 0.0;
604 double laypi1 = -1.0;
605 double laypi2 = -1.0;
606 double laypi3 = -1.0;
607 double laypi4 = -1.0;
608
609 EvtRecTrackIterator itTrkp1 = evtRecTrkCol->begin() + ipip[0];
610 RecMdcTrack* mdcTrkp1 = ( *itTrkp1 )->mdcTrack();
611 RecMdcKalTrack* mdcKalTrkp1 = ( *itTrkp1 )->mdcKalTrack();
612 RecEmcShower* emcTrkp1 = ( *itTrkp1 )->emcShower();
613 RecMucTrack* mucTrkp1 = ( *itTrkp1 )->mucTrack();
614
623 double m_p1pxy =
625
628 double m_laypip1 = laypi1;
629
630 EvtRecTrackIterator itTrkm1 = evtRecTrkCol->begin() + ipim[0];
631 RecMdcTrack* mdcTrkm1 = ( *itTrkm1 )->mdcTrack();
632 RecMdcKalTrack* mdcKalTrkm1 = ( *itTrkm1 )->mdcKalTrack();
633 RecEmcShower* emcTrkm1 = ( *itTrkm1 )->emcShower();
634 RecMucTrack* mucTrkm1 = ( *itTrkm1 )->mucTrack();
635
644 double m_m1pxy =
646
649 double m_laypim1 = laypi2;
650
651 EvtRecTrackIterator itTrkp2 = evtRecTrkCol->begin() + ipip[1];
652 RecMdcTrack* mdcTrkp2 = ( *itTrkp2 )->mdcTrack();
653 RecMdcKalTrack* mdcKalTrkp2 = ( *itTrkp2 )->mdcKalTrack();
654 RecEmcShower* emcTrkp2 = ( *itTrkp2 )->emcShower();
655 RecMucTrack* mucTrkp2 = ( *itTrkp2 )->mucTrack();
656
665 double m_p2pxy =
667
670 double m_laypip2 = laypi3;
671
672 EvtRecTrackIterator itTrkm2 = evtRecTrkCol->begin() + ipim[1];
673 RecMdcTrack* mdcTrkm2 = ( *itTrkm2 )->mdcTrack();
674 RecMdcKalTrack* mdcKalTrkm2 = ( *itTrkm2 )->mdcKalTrack();
675 RecEmcShower* emcTrkm2 = ( *itTrkm2 )->emcShower();
676 RecMucTrack* mucTrkm2 = ( *itTrkm2 )->mucTrack();
677
686 double m_m2pxy =
688
691 double m_laypim2 = laypi4;
692
693 double m_emcTp1 = emcTg1;
694 double m_emcTm1 = emcTg2;
695 double m_emcTp2 = emcTg3;
696 double m_emcTm2 = emcTg4;
697
698 if ( fabs( m_p1vz ) >= m_vz1cut ) return StatusCode::SUCCESS;
699 if ( m_p1vr >= m_vr1cut ) return StatusCode::SUCCESS;
700 if ( fabs( m_p1vct ) >= m_cthcut ) return StatusCode::SUCCESS;
701
702 if ( fabs( m_m1vz ) >= m_vz1cut ) return StatusCode::SUCCESS;
703 if ( m_m1vr >= m_vr1cut ) return StatusCode::SUCCESS;
704 if ( fabs( m_m1vct ) >= m_cthcut ) return StatusCode::SUCCESS;
705 Ncut5++;
706
707 HepLorentzVector p4muonp = ppip[1];
708 HepLorentzVector p4muonm = ppim[1];
709 HepLorentzVector p4uu = pTrecf;
710
711 // Lorentz transformation : boost and rotate
712 Hep3Vector p3jpsiUnit = ( p4uu.vect() ).unit();
713 double jBeta = p4uu.beta(); // just same as the P/E
714
715 // std::cout << jBeta << " " << p4uu.beta() << std::endl;
716
717 //
718 // Loop each gamma pair, check ppi0, pTot
719 // and other mass from MDC momentum
720 //
721
722 HepLorentzVector pTot;
723 double minpi0 = 999.0;
724 for ( int i = 0; i < nGam - 1; i++ )
725 {
726 for ( int j = i + 1; j < nGam; j++ )
727 {
728 HepLorentzVector p2g = pGam[i] + pGam[j];
729 pTot = ppip[0] + ppim[0] + ppip[1] + ppim[1];
730 pTot += p2g;
731 if ( fabs( p2g.m() - 0.135 ) < minpi0 )
732 {
733 minpi0 = fabs( p2g.m() - 0.135 );
734 // 2009//4
735 double m_m2gg = p2g.m();
736 // 2009//4
737 HepLorentzVector prho0_no = ppi2_no2;
738 HepLorentzVector prhop_no = ppip[1] + p2g;
739 HepLorentzVector prhom_no = ppim[1] + p2g;
740 HepLorentzVector prho0pi0 = ppi2_no2 + p2g;
741 HepLorentzVector frho1pi0 = ppi2_no1 + p2g;
742 HepLorentzVector frho2pi0 = ppi2_no3 + p2g;
743 HepLorentzVector frho3pi0 = ppi2_no4 + p2g;
744 HepLorentzVector prho0g1 = ppi2_no2 + pGam[i];
745 HepLorentzVector prho0g2 = ppi2_no2 + pGam[j];
746 HepLorentzVector frho1g1 = ppi2_no1 + pGam[i];
747 HepLorentzVector frho1g2 = ppi2_no1 + pGam[j];
748 HepLorentzVector frho2g1 = ppi2_no3 + pGam[i];
749 HepLorentzVector frho2g2 = ppi2_no3 + pGam[j];
750 HepLorentzVector frho3g1 = ppi2_no4 + pGam[i];
751 HepLorentzVector frho3g2 = ppi2_no4 + pGam[j];
752 HepLorentzVector p5pi_no = p4pi_no + p2g;
753
754 // 2009//4
755 double m_prho0_no = prho0_no.m();
756 double m_prhop_no = prhop_no.m();
757 double m_prhom_no = prhom_no.m();
758 double m_prho0pi0 = prho0pi0.m();
759 double m_frho1pi0 = frho1pi0.m();
760 double m_frho2pi0 = frho2pi0.m();
761 double m_frho3pi0 = frho3pi0.m();
762 double m_prho0g1 = prho0g1.m();
763 double m_prho0g2 = prho0g2.m();
764 double m_frho1g1 = frho1g1.m();
765 double m_frho1g2 = frho1g2.m();
766 double m_frho2g1 = frho2g1.m();
767 double m_frho2g2 = frho2g2.m();
768 double m_frho3g1 = frho3g1.m();
769 double m_frho3g2 = frho3g2.m();
770 double m_p4pi_no = p4pi_no.m();
771 double m_p5pi_no = p5pi_no.m();
772 double m_mdcpx1 = ppip[0].px();
773 double m_mdcpy1 = ppip[0].py();
774 double m_mdcpz1 = ppip[0].pz();
775 double m_mdcpe1 = ppip[0].e();
776 double m_mdcpx2 = ppim[0].px();
777 double m_mdcpy2 = ppim[0].py();
778 double m_mdcpz2 = ppim[0].pz();
779 double m_mdcpe2 = ppim[0].e();
780 double m_mdcpx3 = ppip[1].px();
781 double m_mdcpy3 = ppip[1].py();
782 double m_mdcpz3 = ppip[1].pz();
783 double m_mdcpe3 = ppip[1].e();
784 double m_mdcpx4 = ppim[1].px();
785 double m_mdcpy4 = ppim[1].py();
786 double m_mdcpz4 = ppim[1].pz();
787 double m_mdcpe4 = ppim[1].e();
788 double m_mdcpxg1 = pGam[i].px();
789 double m_mdcpyg1 = pGam[i].py();
790 double m_mdcpzg1 = pGam[i].pz();
791 double m_mdcpeg1 = pGam[i].e();
792 double m_mdcpxg2 = pGam[j].px();
793 double m_mdcpyg2 = pGam[j].py();
794 double m_mdcpzg2 = pGam[j].pz();
795 double m_mdcpeg2 = pGam[j].e();
796 double m_etot = pTot.e();
797 double m_mrecjp1 = m_recjpsi1;
798 double m_mrecjp2 = m_recjpsi2;
799 double m_mrecjp3 = m_recjpsi3;
800 double m_mrecjp4 = m_recjpsi4;
801 // m_tuple3 -> write();
802 // 2009//4
803 }
804 }
805 }
806 Ncut6++;
807
808 //
809 // Test vertex fit
810 //
811
812 HepPoint3D vx( 0., 0., 0. );
813 HepSymMatrix Evx( 3, 0 );
814 double bx = 1E+6;
815 double by = 1E+6;
816 double bz = 1E+6;
817 Evx[0][0] = bx * bx;
818 Evx[1][1] = by * by;
819 Evx[2][2] = bz * bz;
820
821 VertexParameter vxpar;
822 vxpar.setVx( vx );
823 vxpar.setEvx( Evx );
824
825 VertexFit* vtxfit = VertexFit::instance();
826 vtxfit->init();
827
828 // assume charged tracks to pi
829
830 RecMdcKalTrack* pipTrk1 = ( *( evtRecTrkCol->begin() + ipip[0] ) )->mdcKalTrack();
831 RecMdcKalTrack* pimTrk1 = ( *( evtRecTrkCol->begin() + ipim[0] ) )->mdcKalTrack();
832 RecMdcKalTrack* pipTrk2 = ( *( evtRecTrkCol->begin() + ipip[1] ) )->mdcKalTrack();
833 RecMdcKalTrack* pimTrk2 = ( *( evtRecTrkCol->begin() + ipim[1] ) )->mdcKalTrack();
834
835 WTrackParameter wvpipTrk1, wvpimTrk1, wvpipTrk2, wvpimTrk2;
840
841 vtxfit->AddTrack( 0, wvpipTrk1 );
842 vtxfit->AddTrack( 1, wvpimTrk1 );
843 vtxfit->AddVertex( 0, vxpar, 0, 1 );
845 vtxfit->Swim( 0 );
846
847 Ncut7++;
848
849 WTrackParameter wpip1 = vtxfit->wtrk( 0 );
850 WTrackParameter wpim1 = vtxfit->wtrk( 1 );
851
852 KinematicFit* kmfit = KinematicFit::instance();
853
854 //
855 // Apply Kinematic 4C fit
856 //
857 int igbf1 = -1;
858 int igbf2 = -1;
859 HepLorentzVector pTgam1( 0, 0, 0, 0 );
860 HepLorentzVector pTgam2( 0, 0, 0, 0 );
861
862 if ( m_test4C == 1 )
863 {
864 // double ecms = 3.097;
865 HepLorentzVector ecms( 0.011 * 3.6862, 0, 0, 3.6862 );
866
867 //
868 // kinematic fit to pi pi K K pi0
869 //
870 WTrackParameter wvkipTrk2, wvkimTrk2;
873 double chisq = 9999.;
874 int ig11 = -1;
875 int ig21 = -1;
876 double chikk = 9999.;
877 for ( int i = 0; i < nGam - 1; i++ )
878 {
879 RecEmcShower* g1Trk = ( *( evtRecTrkCol->begin() + iGam[i] ) )->emcShower();
880 for ( int j = i + 1; j < nGam; j++ )
881 {
882 RecEmcShower* g2Trk = ( *( evtRecTrkCol->begin() + iGam[j] ) )->emcShower();
883 kmfit->init();
884 kmfit->setDynamicerror( 1 );
885 kmfit->AddTrack( 0, wpip1 );
886 kmfit->AddTrack( 1, wpim1 );
887 kmfit->AddTrack( 2, wvkipTrk2 );
888 kmfit->AddTrack( 3, wvkimTrk2 );
889 kmfit->AddTrack( 4, 0.0, g1Trk );
890 kmfit->AddTrack( 5, 0.0, g2Trk );
894 }
895 }
896 Ncut8++;
897
898 //
899 // kinematic fit to pi pi pi pi pi0
900 //
901
902 chisq = 9999.;
903 int ig1 = -1;
904 int ig2 = -1;
905 for ( int i = 0; i < nGam - 1; i++ )
906 {
907 RecEmcShower* g1Trk = ( *( evtRecTrkCol->begin() + iGam[i] ) )->emcShower();
908 for ( int j = i + 1; j < nGam; j++ )
909 {
910 RecEmcShower* g2Trk = ( *( evtRecTrkCol->begin() + iGam[j] ) )->emcShower();
911 kmfit->init();
912 kmfit->setDynamicerror( 1 );
913 kmfit->AddTrack( 0, wpip1 );
914 kmfit->AddTrack( 1, wpim1 );
915 kmfit->AddTrack( 2, wvpipTrk2 );
916 kmfit->AddTrack( 3, wvpimTrk2 );
917 kmfit->AddTrack( 4, 0.0, g1Trk );
918 kmfit->AddTrack( 5, 0.0, g2Trk );
921 if ( oksq )
922 {
924 if ( chi2 < chisq )
925 {
926 chisq = chi2;
927 ig1 = iGam[i];
928 ig2 = iGam[j];
929 igbf1 = iGam[i];
930 igbf2 = iGam[j];
931 pTgam1 = pGam[i];
932 pTgam2 = pGam[j];
933 }
934 }
935 }
936 }
937 // log << MSG::DEBUG << "photon ID from 4c fit to 4pi+pi0 " << ig1 << " , "
938 // << ig2 << endmsg;
939 if ( chisq > 200 ) return StatusCode::SUCCESS;
940 Ncut9++;
941
942 // select charge track and nneu track
943 Vint jGood;
944 jGood.clear();
945 jGood.push_back( ipip[0] );
946 jGood.push_back( ipim[0] );
947 jGood.push_back( ipip[1] );
948 jGood.push_back( ipim[1] );
949
950 Vint jGgam;
951 jGgam.clear();
952 jGgam.push_back( igbf1 );
953 jGgam.push_back( igbf2 );
954
955 double chi1_pp = 9999.0;
956
957 RecEmcShower* g1Trk = ( *( evtRecTrkCol->begin() + ig1 ) )->emcShower();
958 RecEmcShower* g2Trk = ( *( evtRecTrkCol->begin() + ig2 ) )->emcShower();
959 kmfit->init();
960 kmfit->AddTrack( 0, wpip1 );
961 kmfit->AddTrack( 1, wpim1 );
962 kmfit->AddTrack( 2, wvpipTrk2 );
963 kmfit->AddTrack( 3, wvpimTrk2 );
964 kmfit->AddTrack( 4, 0.0, g1Trk );
965 kmfit->AddTrack( 5, 0.0, g2Trk );
968 if ( oksq )
969 {
970 chi1_pp = kmfit->chisq();
973 HepLorentzVector prhop = kmfit->pfit( 2 ) + ppi0;
974 HepLorentzVector prhom = kmfit->pfit( 3 ) + ppi0;
975 HepLorentzVector pjjj = prho0 + ppi0;
976
978 HepLorentzVector f2pi1g1 = p2pi1 + kmfit->pfit( 4 );
979 HepLorentzVector f2pi1g2 = p2pi1 + kmfit->pfit( 5 );
980 HepLorentzVector f2pi1pi0 = p2pi1 + ppi0;
981
982 HepLorentzVector t2pi2g1 = prho0 + kmfit->pfit( 4 );
983 HepLorentzVector t2pi2g2 = prho0 + kmfit->pfit( 5 );
984
986 HepLorentzVector f2pi3g1 = p2pi3 + kmfit->pfit( 4 );
987 HepLorentzVector f2pi3g2 = p2pi3 + kmfit->pfit( 5 );
988 HepLorentzVector f2pi3pi0 = p2pi3 + ppi0;
989
991 HepLorentzVector f2pi4g1 = p2pi4 + kmfit->pfit( 4 );
992 HepLorentzVector f2pi4g2 = p2pi4 + kmfit->pfit( 5 );
993 HepLorentzVector f2pi4pi0 = p2pi4 + ppi0;
994
995 HepLorentzVector p4pi = p2pi1 + prho0;
996 HepLorentzVector p4pig1 = p4pi + kmfit->pfit( 4 );
997 HepLorentzVector p4pig2 = p4pi + kmfit->pfit( 5 );
998 HepLorentzVector ppptot = p4pi + ppi0;
999
1000 // add
1001 HepLorentzVector be4cpi0 = pTgam1 + pTgam2;
1002 HepLorentzVector be4c_ppi1 = ppip[0] + ppim[0];
1003 HepLorentzVector be4c_ppi2 = ppip[1] + ppim[1];
1004 HepLorentzVector be4cjp = be4cpi0 + be4c_ppi2;
1005
1006 //**********************************//
1007 // final event selection //
1008 // for pion control sample //
1009 //**********************************//
1010 /*
1011 if(fabs(ppi0.m()-0.135)>0.02) return StatusCode::SUCCESS;
1012 if(fabs(m_recppf-3.097)>0.01) return StatusCode::SUCCESS;
1013 if(fabs(m_emcTp2+m_emcTm2)>2.6) return StatusCode::SUCCESS;
1014 if(chi1_pp>chikk) return StatusCode::SUCCESS;
1015
1016 */
1017 //**********************************//
1018
1019 /*
1020 m_run = eventHeader->runNumber();
1021 m_rec = eventHeader->eventNumber();
1022 m_nch = evtRecEvent->totalCharged();
1023 m_nneu = evtRecEvent->totalNeutral();
1024
1025 m_gdgam=nGam;
1026 m_recpp=m_recppf;
1027 m_ngch = jGood.size();
1028 m_nggneu=jGgam.size();
1029
1030 m_chi1=chi1_pp;
1031 m_bepi0=be4cpi0.m();
1032 m_be4cjpsi = be4cjp.m();
1033 //
1034 m_mpi0 = ppi0.m();
1035 m_mprho0=prho0.m();
1036 m_mprhop=prhop.m();
1037 m_mprhom=prhom.m();
1038 m_mpjjj =pjjj.m();
1039 m_mp2pi1 = p2pi1.m();
1040 m_mf2pi1g1=f2pi1g1.m();
1041 m_mf2pi1g2=f2pi1g2.m();
1042 m_mf2pi1pi0=f2pi1pi0.m();
1043 m_mt2pi2g1=t2pi2g1.m();
1044 m_mt2pi2g2=t2pi2g2.m();
1045 m_mp2pi3 = p2pi3.m();
1046 m_mf2pi3g1=f2pi3g1.m();
1047 m_mf2pi3g2=f2pi3g2.m();
1048 m_mf2pi3pi0=f2pi3pi0.m();
1049 m_mp2pi4 = p2pi4.m();
1050 m_mf2pi4g1=f2pi4g1.m();
1051 m_mf2pi4g2=f2pi4g2.m();
1052 m_mf2pi4pi0=f2pi4pi0.m();
1053 m_mp4pi = p4pi.m();
1054 m_mppptot=ppptot.m();
1055 m_mp4pig1=p4pig1.m();
1056 m_mp4pig2=p4pig2.m();
1057 m_mpx1=kmfit->pfit(0).px();
1058 m_mpy1=kmfit->pfit(0).py();
1059 m_mpz1=kmfit->pfit(0).pz();
1060 m_mpe1=kmfit->pfit(0).e();
1061 m_mpx2=kmfit->pfit(1).px();
1062 m_mpy2=kmfit->pfit(1).py();
1063 m_mpz2=kmfit->pfit(1).pz();
1064 m_mpe2=kmfit->pfit(1).e();
1065 m_mpx3=kmfit->pfit(2).px();
1066 m_mpy3=kmfit->pfit(2).py();
1067 m_mpz3=kmfit->pfit(2).pz();
1068 m_mpe3=kmfit->pfit(2).e();
1069 m_mpx4=kmfit->pfit(3).px();
1070 m_mpy4=kmfit->pfit(3).py();
1071 m_mpz4=kmfit->pfit(3).pz();
1072 m_mpe4=kmfit->pfit(3).e();
1073 m_mpxg1=kmfit->pfit(4).px();
1074 m_mpyg1=kmfit->pfit(4).py();
1075 m_mpzg1=kmfit->pfit(4).pz();
1076 m_mpeg1=kmfit->pfit(4).e();
1077 m_mpxg2=kmfit->pfit(5).px();
1078 m_mpyg2=kmfit->pfit(5).py();
1079 m_mpzg2=kmfit->pfit(5).pz();
1080 m_mpeg2=kmfit->pfit(5).e();
1081 m_good = nGood;
1082 m_chikk=chikk;
1083 m_gam = nGam;
1084 m_pip = npip;
1085 m_pim = npim;
1086
1087 //
1088 // fill information of dedx and tof
1089 //
1090
1091 for(int jk = 0; jk < 2; jk++) {
1092 m_ipipin[jk]=0;
1093 m_ipimin[jk]=0;
1094 }
1095 int ipidpip=0;
1096 int ipidpim=0;
1097
1098 ParticleID *pid = ParticleID::instance();
1099 for(int i = 0; i < npip; i++) {
1100 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + ipip[i];
1101 pid->init();
1102 pid->setMethod(pid->methodProbability());
1103 pid->setChiMinCut(4);
1104 pid->setRecTrack(*itTrk);
1105 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2()); // use PID
1106 sub-system pid->identify(pid->onlyPion() | pid->onlyKaon()); // seperater Pion/Kaon
1107 pid->calculate();
1108 if(!(pid->IsPidInfoValid())) continue;
1109 if(pid->probPion() > pid->probKaon()) {
1110 m_ipipin[ipidpip]=1;
1111 ipidpip++;
1112 }
1113 }
1114
1115 // ParticleID *pid = ParticleID::instance();
1116 for(int j = 0; j < npim; j++) {
1117 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + ipim[j];
1118 pid->init();
1119 pid->setMethod(pid->methodProbability());
1120 pid->setChiMinCut(4);
1121 pid->setRecTrack(*jtTrk);
1122 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2()); // use PID
1123 sub-system pid->identify(pid->onlyPion() | pid->onlyKaon()); // seperater Pion/Kaon
1124 pid->calculate();
1125 if(!(pid->IsPidInfoValid())) continue;
1126 if(pid->probPion() > pid->probKaon()) {
1127 m_ipimin[ipidpim]=1;
1128 ipidpim++;
1129 }
1130 }
1131
1132 m_pidpip=ipidpip;
1133 m_pidpim=ipidpim;
1134
1135 //
1136 // check dedx infomation
1137 //
1138
1139 for(int ii = 0; ii < 4; ii++) {
1140 // dedx
1141 m_ptrk[ii] = 9999.0;
1142 m_chie[ii] = 9999.0;
1143 m_chimu[ii] = 9999.0;
1144 m_chipi[ii] = 9999.0;
1145 m_chik[ii] = 9999.0;
1146 m_chip[ii] = 9999.0;
1147 m_ghit[ii] = 9999.0;
1148 m_thit[ii] = 9999.0;
1149 m_probPH[ii] = 9999.0;
1150 m_normPH[ii] = 9999.0;
1151
1152 //endtof
1153 m_cntr_etof[ii] = 9999.0;
1154 m_ptot_etof[ii] = 9999.0;
1155 m_ph_etof[ii] = 9999.0;
1156 m_rhit_etof[ii] = 9999.0;
1157 m_qual_etof[ii] = 9999.0;
1158 m_te_etof[ii] = 9999.0;
1159 m_tmu_etof[ii] = 9999.0;
1160 m_tpi_etof[ii] = 9999.0;
1161 m_tk_etof[ii] = 9999.0;
1162 m_tp_etof[ii] = 9999.0;
1163 m_ec_tof[ii] = 9999.0;
1164 m_ec_toff_e[ii] = 9999.0;
1165 m_ec_toff_mu[ii] = 9999.0;
1166 m_ec_toff_pi[ii] = 9999.0;
1167 m_ec_toff_k[ii] = 9999.0;
1168 m_ec_toff_p[ii] = 9999.0;
1169 m_ec_tsig_e[ii] = 9999.0;
1170 m_ec_tsig_mu[ii] = 9999.0;
1171 m_ec_tsig_pi[ii] = 9999.0;
1172 m_ec_tsig_k[ii] = 9999.0;
1173 m_ec_tsig_p[ii] = 9999.0;
1174
1175 // barrel tof
1176 m_cntr_btof1[ii] = 9999.0;
1177 m_ptot_btof1[ii] = 9999.0;
1178 m_ph_btof1[ii] = 9999.0;
1179 m_zhit_btof1[ii] = 9999.0;
1180 m_qual_btof1[ii] = 9999.0;
1181 m_te_btof1[ii] = 9999.0;
1182 m_tmu_btof1[ii] = 9999.0;
1183 m_tpi_btof1[ii] = 9999.0;
1184 m_tk_btof1[ii] = 9999.0;
1185 m_tp_btof1[ii] = 9999.0;
1186 m_b1_tof[ii] = 9999.0;
1187 m_b1_toff_e[ii] = 9999.0;
1188 m_b1_toff_mu[ii] = 9999.0;
1189 m_b1_toff_pi[ii] = 9999.0;
1190 m_b1_toff_k[ii] = 9999.0;
1191 m_b1_toff_p[ii] = 9999.0;
1192 m_b1_tsig_e[ii] = 9999.0;
1193 m_b1_tsig_mu[ii] = 9999.0;
1194 m_b1_tsig_pi[ii] = 9999.0;
1195 m_b1_tsig_k[ii] = 9999.0;
1196 m_b1_tsig_p[ii] = 9999.0;
1197 //pid
1198 m_dedx_pid[ii] = 9999.0;
1199 m_tof1_pid[ii] = 9999.0;
1200 m_tof2_pid[ii] = 9999.0;
1201 m_prob_pid[ii] = 9999.0;
1202 m_ptrk_pid[ii] = 9999.0;
1203 m_cost_pid[ii] = 9999.0;
1204 }
1205
1206
1207 int indx0=0;
1208 for(int i = 0; i < m_ngch; i++) {
1209 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + jGood[i];
1210 if(!(*itTrk)->isMdcTrackValid()) continue;
1211 if(!(*itTrk)->isMdcDedxValid())continue;
1212 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
1213 RecMdcDedx* dedxTrk = (*itTrk)->mdcDedx();
1214 m_ptrk[indx0] = mdcTrk->p();
1215
1216 m_chie[indx0] = dedxTrk->chiE();
1217 m_chimu[indx0] = dedxTrk->chiMu();
1218 m_chipi[indx0] = dedxTrk->chiPi();
1219 m_chik[indx0] = dedxTrk->chiK();
1220 m_chip[indx0] = dedxTrk->chiP();
1221 m_ghit[indx0] = dedxTrk->numGoodHits();
1222 m_thit[indx0] = dedxTrk->numTotalHits();
1223 m_probPH[indx0] = dedxTrk->probPH();
1224 m_normPH[indx0] = dedxTrk->normPH();
1225 indx0++;
1226 // m_tuple7->write();
1227 }
1228
1229
1230 //
1231 // check TOF infomation
1232 //
1233
1234
1235 int indx1=0;
1236 for(int i = 0; i < m_ngch; i++) {
1237 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + jGood[i];
1238 if(!(*itTrk)->isMdcTrackValid()) continue;
1239 if(!(*itTrk)->isTofTrackValid()) continue;
1240
1241 RecMdcTrack * mdcTrk = (*itTrk)->mdcTrack();
1242 SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
1243
1244 double ptrk = mdcTrk->p();
1245 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
1246 for(;iter_tof != tofTrkCol.end(); iter_tof++ ) {
1247 TofHitStatus *status = new TofHitStatus;
1248 status->setStatus((*iter_tof)->status());
1249 if(!(status->is_barrel())){//endcap
1250 if( !(status->is_counter()) ) continue; // ?
1251 if( status->layer()!=1 ) continue;//layer1
1252 double path=(*iter_tof)->path(); // ?
1253 double tof = (*iter_tof)->tof();
1254 double ph = (*iter_tof)->ph();
1255 double rhit = (*iter_tof)->zrhit();
1256 double qual = 0.0 + (*iter_tof)->quality();
1257 double cntr = 0.0 + (*iter_tof)->tofID();
1258 double texp[5];
1259 double tsig[5];
1260 for(int j = 0; j < 5; j++) {//0 e, 1 mu, 2 pi, 3 K, 4 p
1261 texp[j] = (*iter_tof)->texp(j);
1262 // tsig[j] = (*iter_tof)->sigma(j);
1263 // toffset[j] = (*iter_tof)->offset(j);
1264 }
1265 m_cntr_etof[indx1] = cntr;
1266 m_ptot_etof[indx1] = ptrk;
1267 m_ph_etof[indx1] = ph;
1268 m_rhit_etof[indx1] = rhit;
1269 m_qual_etof[indx1] = qual;
1270 m_te_etof[indx1] = tof - texp[0];
1271 m_tmu_etof[indx1] = tof - texp[1];
1272 m_tpi_etof[indx1] = tof - texp[2];
1273 m_tk_etof[indx1] = tof - texp[3];
1274 m_tp_etof[indx1] = tof - texp[4];
1275
1276 m_ec_tof[indx1] = tof;
1277
1278 m_ec_toff_e[indx1] = (*iter_tof)->toffset(0);
1279 m_ec_toff_mu[indx1] = (*iter_tof)->toffset(1);
1280 m_ec_toff_pi[indx1] = (*iter_tof)->toffset(2);
1281 m_ec_toff_k[indx1] = (*iter_tof)->toffset(3);
1282 m_ec_toff_p[indx1] = (*iter_tof)->toffset(4);
1283
1284 m_ec_tsig_e[indx1] = (*iter_tof)->sigma(0);
1285 m_ec_tsig_mu[indx1] = (*iter_tof)->sigma(1);
1286 m_ec_tsig_pi[indx1] = (*iter_tof)->sigma(2);
1287 m_ec_tsig_k[indx1] = (*iter_tof)->sigma(3);
1288 m_ec_tsig_p[indx1] = (*iter_tof)->sigma(4);
1289
1290 // m_tuple8->write();
1291 }
1292 else {//barrel
1293 if( !(status->is_cluster()) ) continue; // ?
1294 double path=(*iter_tof)->path(); // ?
1295 double tof = (*iter_tof)->tof();
1296 double ph = (*iter_tof)->ph();
1297 double rhit = (*iter_tof)->zrhit();
1298 double qual = 0.0 + (*iter_tof)->quality();
1299 double cntr = 0.0 + (*iter_tof)->tofID();
1300 double texp[5];
1301 for(int j = 0; j < 5; j++) {
1302 texp[j] = (*iter_tof)->texp(j);
1303 }
1304 m_cntr_btof1[indx1] = cntr;
1305 m_ptot_btof1[indx1] = ptrk;
1306 m_ph_btof1[indx1] = ph;
1307 m_zhit_btof1[indx1] = rhit;
1308 m_qual_btof1[indx1] = qual;
1309 m_te_btof1[indx1] = tof - texp[0];
1310 m_tmu_btof1[indx1] = tof - texp[1];
1311 m_tpi_btof1[indx1] = tof - texp[2];
1312 m_tk_btof1[indx1] = tof - texp[3];
1313 m_tp_btof1[indx1] = tof - texp[4];
1314
1315 m_b1_tof[indx1] = tof;
1316
1317 m_b1_toff_e[indx1] = (*iter_tof)->toffset(0);
1318 m_b1_toff_mu[indx1] = (*iter_tof)->toffset(1);
1319 m_b1_toff_pi[indx1] = (*iter_tof)->toffset(2);
1320 m_b1_toff_k[indx1] = (*iter_tof)->toffset(3);
1321 m_b1_toff_p[indx1] = (*iter_tof)->toffset(4);
1322
1323 m_b1_tsig_e[indx1] = (*iter_tof)->sigma(0);
1324 m_b1_tsig_mu[indx1] = (*iter_tof)->sigma(1);
1325 m_b1_tsig_pi[indx1] = (*iter_tof)->sigma(2);
1326 m_b1_tsig_k[indx1] = (*iter_tof)->sigma(3);
1327 m_b1_tsig_p[indx1] = (*iter_tof)->sigma(4);
1328
1329 // m_tuple9->write();
1330 }
1331 delete status;
1332 }
1333 indx1++;
1334 } // loop all charged track
1335
1336 //
1337 // Assign 4-momentum to each charged track
1338 //
1339 int indx2=0;
1340 // ParticleID *pid = ParticleID::instance();
1341 for(int i = 0; i < m_ngch; i++) {
1342 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + jGood[i];
1343 // if(pid) delete pid;
1344 pid->init();
1345 pid->setMethod(pid->methodProbability());
1346 // pid->setMethod(pid->methodLikelihood()); //for Likelihood Method
1347
1348 pid->setChiMinCut(4);
1349 pid->setRecTrack(*itTrk);
1350 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2()); // use PID
1351 sub-system pid->identify(pid->onlyPion() | pid->onlyKaon()); // seperater Pion/Kaon
1352 // pid->identify(pid->onlyPion());
1353 // pid->identify(pid->onlyKaon());
1354 pid->calculate();
1355 if(!(pid->IsPidInfoValid())) continue;
1356 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
1357
1358 m_dedx_pid[indx2] = pid->chiDedx(2);
1359 m_tof1_pid[indx2] = pid->chiTof1(2);
1360 m_tof2_pid[indx2] = pid->chiTof2(2);
1361 m_prob_pid[indx2] = pid->probPion();
1362
1363 // if(pid->probPion() < 0.001 || (pid->probPion() < pid->probKaon())) continue;
1364 // if(pid->probPion() < 0.001) continue;
1365 // if(pid->pdf(2)<pid->pdf(3)) continue; // for Likelihood Method(0=electron 1=muon
1366 2=pion 3=kaon 4=proton)
1367
1368 RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();//After ParticleID, use
1369 RecMdcKalTrack substitute RecMdcTrack RecMdcKalTrack::setPidType
1370 (RecMdcKalTrack::pion);//PID can set to electron, muon, pion, kaon and proton;The default
1371 setting is pion
1372
1373 m_ptrk_pid[indx2] = mdcKalTrk->p();
1374 m_cost_pid[indx2] = cos(mdcTrk->theta());
1375 }
1376
1377
1378 int iphoton = 0;
1379 for (int i=0; i<m_nggneu; i++)
1380 {
1381 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + jGgam[i];
1382 if (!(*itTrk)->isEmcShowerValid()) continue;
1383 RecEmcShower *emcTrk = (*itTrk)->emcShower();
1384 m_numHits[iphoton] = emcTrk->numHits();
1385 m_secondmoment[iphoton] = emcTrk->secondMoment();
1386 m_x[iphoton] = emcTrk->x();
1387 m_y[iphoton]= emcTrk->y();
1388 m_z[iphoton]= emcTrk->z();
1389 m_cosemc[iphoton] = cos(emcTrk->theta());
1390 m_phiemc[iphoton] = emcTrk->phi();
1391 m_energy[iphoton] = emcTrk->energy();
1392 m_eSeed[iphoton] = emcTrk->eSeed();
1393 m_e3x3[iphoton] = emcTrk->e3x3();
1394 m_e5x5[iphoton] = emcTrk->e5x5();
1395 m_lat[iphoton] = emcTrk->latMoment();
1396 m_a20[iphoton] = emcTrk->a20Moment();
1397 m_a42[iphoton] = emcTrk->a42Moment();
1398 iphoton++;
1399 }
1400 // m_tuple4->write();
1401 */
1402 Ncut10++;
1403 // log << MSG::DEBUG << "chisquare from 4c fit to 4pi+pi0 " << m_chi1 << endmsg;
1404 }
1405 }
1406
1407 //
1408 // Apply Kinematic 5C Fit
1409 //
1410
1411 // find the best combination over all possible pi+ pi- gamma gamma pair
1412
1413 setFilterPassed( true );
1414
1415 ( *( evtRecTrkCol->begin() + ipip[0] ) )->setPartId( 2 );
1416 ( *( evtRecTrkCol->begin() + ipim[0] ) )->setPartId( 2 );
1417 ( *( evtRecTrkCol->begin() + ipip[1] ) )->setPartId( 2 );
1418 ( *( evtRecTrkCol->begin() + ipim[1] ) )->setPartId( 2 );
1419
1420 return StatusCode::SUCCESS;
1421}
EvtRecTrackCol::iterator EvtRecTrackIterator
double sin(const BesAngle a)
double cos(const BesAngle a)
double energy() const
const Hep3Vector emcPosition() const
const int emcVolumeNumber() const
const double px() const
static void setPidType(PidType pidType)
const double pz() const
const double py() const
const double p() const
const int charge() const
const double theta() const
const int charge() const
int numLayers() const
virtual bool isVertexValid()=0
virtual double * SigmaPrimaryVertex()=0
virtual double * PrimaryVertex()=0
void setDynamicerror(const bool dynamicerror=1)
double chisq() const
void AddFourMomentum(int number, HepLorentzVector p4)
Definition KinematicFit.cxx:341
HepLorentzVector pfit(int n) const
const HepVector & getZHelix() const
HepVector & getZHelixK()
const HepSymMatrix & getZError() const
HepSymMatrix & getZErrorK()
void AddTrack(const int number, const double mass, const RecMdcTrack *trk)
Definition TrackPool.cxx:21
WTrackParameter wtrk(int n) const
void AddVertex(int number, VertexParameter vpar, std::vector< int > lis)
Definition VertexFit.cxx:85
void setEvx(const HepSymMatrix &eVx)
void setVx(const HepPoint3D &vx)
_EXTERN_ std::string EvtRecEvent
_EXTERN_ std::string EvtRecTrackCol
◆ finalize()
| StatusCode Ppjrhopi::finalize | ( | ) |
Definition at line 1424 of file Ppjrhopi.cxx.
1424 {
1436 MsgStream log( msgSvc(), name() );
1437 log << MSG::INFO << "in finalize()" << endmsg;
1438 return StatusCode::SUCCESS;
1439}
◆ initialize()
| StatusCode Ppjrhopi::initialize | ( | ) |
Definition at line 58 of file Ppjrhopi.cxx.
58 {
59 MsgStream log( msgSvc(), name() );
60
61 log << MSG::INFO << "in initialize()" << endmsg;
62
63 StatusCode status;
64
65 if ( m_test4C == 1 )
66 {
67 /*
68 NTuplePtr nt4(ntupleSvc(), "FILE1/fit4c");
69 if ( nt4 ) m_tuple4 = nt4;
70 else {
71 m_tuple4 = ntupleSvc()->book ("FILE1/fit4c", CLID_ColumnWiseTuple, "ks N-Tuple
72 example"); if ( m_tuple4 ) { status = m_tuple4->addItem ("run", m_run); status =
73 m_tuple4->addItem ("rec", m_rec); status = m_tuple4->addItem ("nch", m_nch);
74 status = m_tuple4->addItem ("nneu", m_nneu);
75 status = m_tuple4->addItem ("mgdgam", m_gdgam);
76 status = m_tuple4->addItem ("recpp", m_recpp);
77 status = m_tuple4->addItem ("chi2", m_chi1);
78 status = m_tuple4->addItem ("mpi0", m_mpi0);
79 status = m_tuple4->addItem ("mprho0", m_mprho0);
80 status = m_tuple4->addItem ("mprhop", m_mprhop);
81 status = m_tuple4->addItem ("mprhom", m_mprhom);
82 status = m_tuple4->addItem ("mpjjj", m_mpjjj);
83 status = m_tuple4->addItem ("mbepi0", m_bepi0);
84 status = m_tuple4->addItem ("mbe4cjpsi",m_be4cjpsi);
85 status = m_tuple4->addItem ("mp2pi1", m_mp2pi1);
86 status = m_tuple4->addItem ("mf2pi1g1", m_mf2pi1g1);
87 status = m_tuple4->addItem ("mf2pi1g2", m_mf2pi1g2);
88 status = m_tuple4->addItem ("mf2pi1pi0",m_mf2pi1pi0);
89 status = m_tuple4->addItem ("mt2pi2g1", m_mt2pi2g1);
90 status = m_tuple4->addItem ("mt2pi2g2", m_mt2pi2g2);
91 status = m_tuple4->addItem ("mp2pi3", m_mp2pi3);
92 status = m_tuple4->addItem ("mf2pi3g1", m_mf2pi3g1);
93 status = m_tuple4->addItem ("mf2pi3g2", m_mf2pi3g2);
94 status = m_tuple4->addItem ("mf2pi3pi0",m_mf2pi3pi0);
95 status = m_tuple4->addItem ("mp2pi4", m_mp2pi4);
96 status = m_tuple4->addItem ("mf2pi4g1", m_mf2pi4g1);
97 status = m_tuple4->addItem ("mf2pi4g2", m_mf2pi4g2);
98 status = m_tuple4->addItem ("mf2pi4pi0",m_mf2pi4pi0);
99 status = m_tuple4->addItem ("mp4pi", m_mp4pi);
100 status = m_tuple4->addItem ("mppptot", m_mppptot);
101 status = m_tuple4->addItem ("mp4pig1", m_mp4pig1);
102 status = m_tuple4->addItem ("mp4pig2", m_mp4pig2);
103 status = m_tuple4->addItem ("mpx1", m_mpx1);
104 status = m_tuple4->addItem ("mpy1", m_mpy1);
105 status = m_tuple4->addItem ("mpz1", m_mpz1);
106 status = m_tuple4->addItem ("mpe1", m_mpe1);
107 status = m_tuple4->addItem ("mpx2", m_mpx2);
108 status = m_tuple4->addItem ("mpy2", m_mpy2);
109 status = m_tuple4->addItem ("mpz2", m_mpz2);
110 status = m_tuple4->addItem ("mpe2", m_mpe2);
111 status = m_tuple4->addItem ("mpx3", m_mpx3);
112 status = m_tuple4->addItem ("mpy3", m_mpy3);
113 status = m_tuple4->addItem ("mpz3", m_mpz3);
114 status = m_tuple4->addItem ("mpe3", m_mpe3);
115 status = m_tuple4->addItem ("mpx4", m_mpx4);
116 status = m_tuple4->addItem ("mpy4", m_mpy4);
117 status = m_tuple4->addItem ("mpz4", m_mpz4);
118 status = m_tuple4->addItem ("mpe4", m_mpe4);
119 status = m_tuple4->addItem ("mpxg1", m_mpxg1);
120 status = m_tuple4->addItem ("mpyg1", m_mpyg1);
121 status = m_tuple4->addItem ("mpzg1", m_mpzg1);
122 status = m_tuple4->addItem ("mpeg1", m_mpeg1);
123 status = m_tuple4->addItem ("mpxg2", m_mpxg2);
124 status = m_tuple4->addItem ("mpyg2", m_mpyg2);
125 status = m_tuple4->addItem ("mpzg2", m_mpzg2);
126 status = m_tuple4->addItem ("mpeg2", m_mpeg2);
127 status = m_tuple4->addItem ("chikk", m_chikk);
128 status = m_tuple4->addItem ("p1vx", m_p1vx);
129 status = m_tuple4->addItem ("p1vy", m_p1vy);
130 status = m_tuple4->addItem ("p1vz", m_p1vz);
131 status = m_tuple4->addItem ("p1vr", m_p1vr);
132 status = m_tuple4->addItem ("p1vct", m_p1vct);
133 status = m_tuple4->addItem ("m1vx", m_m1vx);
134 status = m_tuple4->addItem ("m1vy", m_m1vy);
135 status = m_tuple4->addItem ("m1vz", m_m1vz);
136 status = m_tuple4->addItem ("m1vr", m_m1vr);
137 status = m_tuple4->addItem ("m1vct", m_m1vct);
138 status = m_tuple4->addItem ("p2vx", m_p2vx);
139 status = m_tuple4->addItem ("p2vy", m_p2vy);
140 status = m_tuple4->addItem ("p2vz", m_p2vz);
141 status = m_tuple4->addItem ("p2vr", m_p2vr);
142 status = m_tuple4->addItem ("p2vct", m_p2vct);
143 status = m_tuple4->addItem ("m2vx", m_m2vx);
144 status = m_tuple4->addItem ("m2vy", m_m2vy);
145 status = m_tuple4->addItem ("m2vz", m_m2vz);
146 status = m_tuple4->addItem ("m2vr", m_m2vr);
147 status = m_tuple4->addItem ("m2vct", m_m2vct);
148 status = m_tuple4->addItem ("mgood", m_good);
149 status = m_tuple4->addItem ("mgam", m_gam);
150 status = m_tuple4->addItem ("mpip", m_pip);
151 status = m_tuple4->addItem ("mpim", m_pim);
152 status = m_tuple4->addItem ("mp1ptot", m_p1ptot);
153 status = m_tuple4->addItem ("memcTp1", m_emcTp1);
154 status = m_tuple4->addItem ("mm1ptot", m_m1ptot);
155 status = m_tuple4->addItem ("memcTm1", m_emcTm1);
156 status = m_tuple4->addItem ("mp2ptot", m_p2ptot);
157 status = m_tuple4->addItem ("memcTp2", m_emcTp2);
158 status = m_tuple4->addItem ("mm2ptot", m_m2ptot);
159 status = m_tuple4->addItem ("memcTm2", m_emcTm2);
160 status = m_tuple4->addItem ("p1pxy", m_p1pxy);
161 status = m_tuple4->addItem ("m1pxy", m_m1pxy);
162 status = m_tuple4->addItem ("p2pxy", m_p2pxy);
163 status = m_tuple4->addItem ("m2pxy", m_m2pxy);
164
165 status = m_tuple4->addItem ("mpidpip", m_pidpip, 0, 10);
166 status = m_tuple4->addIndexedItem ("mipipin" , m_pidpip, m_ipipin);
167 status = m_tuple4->addItem ("mpidpim", m_pidpim, 0, 10);
168 status = m_tuple4->addIndexedItem ("mipimin" , m_pidpim, m_ipimin);
169
170 status = m_tuple4->addItem ("laypip1", m_laypip1);
171 status = m_tuple4->addItem ("laypim1", m_laypim1);
172 status = m_tuple4->addItem ("laypip2", m_laypip2);
173 status = m_tuple4->addItem ("laypim2", m_laypim2);
174 status = m_tuple4->addItem ("mangle", m_angle);
175 status = m_tuple4->addItem ("cosuubr",m_cosuubr );
176 status = m_tuple4->addItem ("cosmupbr", m_cosmupbr);
177 status = m_tuple4->addItem ("cosmumbr", m_cosmumbr);
178 status = m_tuple4->addItem ("phimupbr", m_phimupbr);
179 status = m_tuple4->addItem ("phimumbr", m_phimumbr);
180 status = m_tuple4->addItem ("ngch", m_ngch, 0, 10);
181 status = m_tuple4->addItem ("nggneu", m_nggneu,0, 10);
182 // modifiey by Zhu
183 // status = m_tuple4->addItem ("indx0", indx0, 0, 10);
184 status = m_tuple4->addIndexedItem ("mptrk" , m_ngch, m_ptrk);
185 status = m_tuple4->addIndexedItem ("chie", m_ngch, m_chie);
186 status = m_tuple4->addIndexedItem ("chimu", m_ngch,m_chimu);
187 status = m_tuple4->addIndexedItem ("chipi", m_ngch,m_chipi);
188 status = m_tuple4->addIndexedItem ("chik", m_ngch,m_chik);
189 status = m_tuple4->addIndexedItem ("chip", m_ngch,m_chip);
190 status = m_tuple4->addIndexedItem ("probPH", m_ngch,m_probPH);
191 status = m_tuple4->addIndexedItem ("normPH", m_ngch,m_normPH);
192 status = m_tuple4->addIndexedItem ("ghit", m_ngch,m_ghit);
193 status = m_tuple4->addIndexedItem ("thit", m_ngch,m_thit);
194
195 status = m_tuple4->addIndexedItem ("ptot_etof", m_ngch,m_ptot_etof);
196 status = m_tuple4->addIndexedItem ("cntr_etof", m_ngch,m_cntr_etof);
197 status = m_tuple4->addIndexedItem ("te_etof", m_ngch,m_te_etof);
198 status = m_tuple4->addIndexedItem ("tmu_etof", m_ngch,m_tmu_etof);
199 status = m_tuple4->addIndexedItem ("tpi_etof", m_ngch,m_tpi_etof);
200 status = m_tuple4->addIndexedItem ("tk_etof", m_ngch,m_tk_etof);
201 status = m_tuple4->addIndexedItem ("tp_etof", m_ngch,m_tp_etof);
202 status = m_tuple4->addIndexedItem ("ph_etof", m_ngch,m_ph_etof);
203 status = m_tuple4->addIndexedItem ("rhit_etof", m_ngch,m_rhit_etof);
204 status = m_tuple4->addIndexedItem ("qual_etof", m_ngch,m_qual_etof);
205 status = m_tuple4->addIndexedItem ("ec_toff_e", m_ngch,m_ec_toff_e);
206 status = m_tuple4->addIndexedItem ("ec_toff_mu",m_ngch,m_ec_toff_mu);
207 status = m_tuple4->addIndexedItem ("ec_toff_pi",m_ngch,m_ec_toff_pi);
208 status = m_tuple4->addIndexedItem ("ec_toff_k", m_ngch,m_ec_toff_k);
209 status = m_tuple4->addIndexedItem ("ec_toff_p", m_ngch,m_ec_toff_p);
210 status = m_tuple4->addIndexedItem ("ec_tsig_e", m_ngch,m_ec_tsig_e);
211 status = m_tuple4->addIndexedItem ("ec_tsig_mu",m_ngch,m_ec_tsig_mu);
212 status = m_tuple4->addIndexedItem ("ec_tsig_pi",m_ngch,m_ec_tsig_pi);
213 status = m_tuple4->addIndexedItem ("ec_tsig_k", m_ngch,m_ec_tsig_k);
214 status = m_tuple4->addIndexedItem ("ec_tsig_p", m_ngch,m_ec_tsig_p);
215 status = m_tuple4->addIndexedItem ("ec_tof", m_ngch,m_ec_tof);
216 status = m_tuple4->addIndexedItem ("ptot_btof1",m_ngch,m_ptot_btof1);
217 status = m_tuple4->addIndexedItem ("cntr_btof1",m_ngch,m_cntr_btof1);
218 status = m_tuple4->addIndexedItem ("te_btof1", m_ngch,m_te_btof1);
219 status = m_tuple4->addIndexedItem ("tmu_btof1", m_ngch,m_tmu_btof1);
220 status = m_tuple4->addIndexedItem ("tpi_btof1", m_ngch,m_tpi_btof1);
221 status = m_tuple4->addIndexedItem ("tk_btof1", m_ngch,m_tk_btof1);
222 status = m_tuple4->addIndexedItem ("tp_btof1", m_ngch,m_tp_btof1);
223 status = m_tuple4->addIndexedItem ("ph_btof1", m_ngch,m_ph_btof1);
224 status = m_tuple4->addIndexedItem ("zhit_btof1",m_ngch,m_zhit_btof1);
225 status = m_tuple4->addIndexedItem ("qual_btof1",m_ngch,m_qual_btof1);
226 status = m_tuple4->addIndexedItem ("b1_toff_e", m_ngch,m_b1_toff_e);
227 status = m_tuple4->addIndexedItem ("b1_toff_mu",m_ngch,m_b1_toff_mu);
228 status = m_tuple4->addIndexedItem ("b1_toff_pi",m_ngch,m_b1_toff_pi);
229 status = m_tuple4->addIndexedItem ("b1_toff_k", m_ngch,m_b1_toff_k);
230 status = m_tuple4->addIndexedItem ("b1_toff_p", m_ngch,m_b1_toff_p);
231 status = m_tuple4->addIndexedItem ("b1_tsig_e", m_ngch,m_b1_tsig_e);
232 status = m_tuple4->addIndexedItem ("b1_tsig_mu",m_ngch,m_b1_tsig_mu);
233 status = m_tuple4->addIndexedItem ("b1_tsig_pi",m_ngch,m_b1_tsig_pi);
234 status = m_tuple4->addIndexedItem ("b1_tsig_k", m_ngch,m_b1_tsig_k);
235 status = m_tuple4->addIndexedItem ("b1_tsig_p", m_ngch,m_b1_tsig_p);
236 status = m_tuple4->addIndexedItem ("b1_tof", m_ngch,m_b1_tof);
237
238 status = m_tuple4->addIndexedItem ("mdedx_pid", m_ngch,m_dedx_pid);
239 status = m_tuple4->addIndexedItem ("mtof1_pid", m_ngch,m_tof1_pid);
240 status = m_tuple4->addIndexedItem ("mtof2_pid", m_ngch,m_tof2_pid);
241 status = m_tuple4->addIndexedItem ("mprob_pid", m_ngch,m_prob_pid);
242 status = m_tuple4->addIndexedItem ("mptrk_pid", m_ngch,m_ptrk_pid);
243 status = m_tuple4->addIndexedItem ("mcost_pid", m_ngch,m_cost_pid);
244
245 status = m_tuple4->addIndexedItem ("numHits", m_nggneu,m_numHits); // Total
246 number of hits status = m_tuple4->addIndexedItem ("secondmoment",
247 m_nggneu,m_secondmoment); status = m_tuple4->addIndexedItem ("mx", m_nggneu,m_x); //
248 Shower coordinates and errors status = m_tuple4->addIndexedItem ("my", m_nggneu,m_y);
249 status = m_tuple4->addIndexedItem ("mz", m_nggneu,m_z);
250 status = m_tuple4->addIndexedItem ("cosemc", m_nggneu,m_cosemc); // Shower
251 Counter angles and errors status = m_tuple4->addIndexedItem ("phiemc",
252 m_nggneu,m_phiemc); status = m_tuple4->addIndexedItem ("energy", m_nggneu,m_energy); //
253 Total energy observed in Emc status = m_tuple4->addIndexedItem ("eseed",
254 m_nggneu,m_eSeed); status = m_tuple4->addIndexedItem ("me9", m_nggneu,m_e3x3);
255 status = m_tuple4->addIndexedItem ("me25", m_nggneu,m_e5x5);
256 status = m_tuple4->addIndexedItem ("mlat", m_nggneu,m_lat);
257 status = m_tuple4->addIndexedItem ("ma20", m_nggneu,m_a20);
258 status = m_tuple4->addIndexedItem ("ma42", m_nggneu,m_a42);
259
260
261 }
262 else {
263 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple4) << endmsg;
264 return StatusCode::FAILURE;
265 }
266 }
267 */
268 } // test 4C
269
270 //
271 //--------end of book--------
272 //
273
274 log << MSG::INFO << "successfully return from initialize()" << endmsg;
275 return StatusCode::SUCCESS;
276}
The documentation for this class was generated from the following files:
