DTagSetAlg Class Reference

#include <DTagSetAlg.h>

Inheritance diagram for DTagSetAlg:

Public Member Functions
	DTagSetAlg (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode	initialize ()
StatusCode	execute ()
StatusCode	finalize ()
unsigned int	IntToTag1 (int val1, int val2, int val3, int val4)
unsigned int	IntToTag2 (int val1, int val2, int val3, int val4, int val5, int val6)

Detailed Description

Definition at line 15 of file DTagSetAlg.h.

Constructor & Destructor Documentation

DTagSetAlg()

DTagSetAlg::DTagSetAlg	(	const std::string &	name,
		ISvcLocator *	pSvcLocator
	)

Definition at line 45 of file DTagSetAlg.cxx.

                                                                      :
        Algorithm(name, pSvcLocator) {
        }

Member Function Documentation

execute()

StatusCode DTagSetAlg::execute

(

)

Definition at line 153 of file DTagSetAlg.cxx.

                              {
 
        int nGoodCharged   = 0;  //1
        int nGoodChargedp = 0;   //2
        int nGoodChargedm = 0;   //3
        int nCharged = 0;        //4
        int nNeutrk  = 0;        //5
        int nTottrk = 0;         //6
        int totCharged = 0;      //7
        int npionp = 0;          //8
        int npionm = 0;          //9
        int nprotonp = 0;        //10
        int nprotonm = 0;        //11
        int nkaonp = 0;          //12
        int nkaonm = 0;          //13
        int nlambda = 0;         //14
        int nalambda= 0;         //15
        int nks = 0;             //16
        int ngamma = 0;          //17
        int neta = 0;            //18
        int npi0 = 0;            //19
        int nmuonp = 0;          //20
        int nmuonm = 0;          //21
        int nelectronp = 0;      //22
        int nelectronm = 0;      //23
 
        Hep3Vector xorigin(0,0,0);
        if(m_vtxsvc->isVertexValid()){
                double* dbv = m_vtxsvc->PrimaryVertex();
                xorigin.setX(dbv[0]);
                xorigin.setY(dbv[1]);
                xorigin.setZ(dbv[2]);
        }
 
        SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
        SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
 
        nCharged = evtRecEvent->totalCharged();
        nNeutrk = evtRecEvent->totalNeutral();
        nTottrk = evtRecEvent->totalTracks();
 
        // Good charged track selection
        vector<int> iGood;
        iGood.clear();
        for ( Int_t iCharge = 0; iCharge < evtRecEvent->totalCharged(); ++iCharge ) {
                EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + iCharge;
                if(!(*itTrk)->isMdcTrackValid()) continue;
                RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
                                double theta = mdcTrk->theta();
                HepVector a = mdcTrk->helix();
                HepSymMatrix Ea = mdcTrk->err();
                HepPoint3D point0(0.,0.,0.);   // the initial point for MDC reconstruction
                HepPoint3D IP(xorigin[0],xorigin[1],xorigin[2]);
                VFHelix helixip(point0,a,Ea);
                helixip.pivot(IP);
                HepVector vecipa = helixip.a();
                double  Rvxy0=fabs(vecipa[0]);  //the nearest distance to IP in xy plane
                double  Rvz0=vecipa[3];         //the nearest distance to IP in z direction
                double  Rvphi0=vecipa[1];
                if(fabs(Rvxy0) >= 1.0) continue;
                if(fabs(Rvz0) >= 10.0) continue;
                                if(fabs(cos(theta))>=0.93) continue;
                if(mdcTrk->charge() > 0) nGoodChargedp++;
                if(mdcTrk->charge() < 0) nGoodChargedm++;
                ++nGoodCharged;
                totCharged +=mdcTrk->charge();
                iGood.push_back(iCharge);
        }
        // std::cout <<"iEvt: "<<iEvt <<"  nGoodCharged: "<<nGoodCharged<<std::endl;
 
        // get photons
        for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
                EvtRecTrackIterator iTrk=evtRecTrkCol->begin() + i;
                if((*iTrk)->isEmcShowerValid()){
                        RecEmcShower* emcShower = (*iTrk)->emcShower();
                        if ( (emcShower->energy()>0.025) && (emcShower->time()>0) && (emcShower->time()<14) && ( ( (abs(cos(emcShower->theta()))<0.80) && (emcShower->energy()>0.025) ) || ( (abs(cos(emcShower->theta()))<0.92) && (abs(cos(emcShower->theta()))>0.86) && (emcShower->energy()>0.050) ) ) ){
                                ngamma++;
                        }
                }
        }
 
        // get tracks
        for (int i=0;i<iGood.size();i++){
                EvtRecTrackIterator iTrk = evtRecTrkCol->begin() + iGood[i];
                RecMdcKalTrack* mdcKalTrack = (*iTrk)->mdcKalTrack();
 
                // set up pid (but only make very rough pid cuts)
                ParticleID* pid = ParticleID::instance();
                pid->init();
                pid->setMethod(pid->methodProbability());
                pid->setRecTrack(*iTrk);
                pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2());
                pid->identify(pid->onlyPion() | pid->onlyKaon() | pid->onlyProton() | pid->onlyElectron() | pid->onlyMuon());
                pid->calculate();
                if(!(pid->IsPidInfoValid())) continue;
 
                if ((pid->probPion() > pid->probKaon()) && (pid->probPion() > pid->probProton())){
                        if (mdcKalTrack->charge() > 0) npionp++;
                        if (mdcKalTrack->charge() < 0) npionm++;
                }
 
                if ((pid->probKaon() > pid->probPion()) && (pid->probKaon() > pid->probProton())){
                        if (mdcKalTrack->charge() > 0) nkaonp++;
                        if (mdcKalTrack->charge() < 0) nkaonm++;
                }
 
                if ((pid->probProton() > pid->probPion()) && (pid->probProton() > pid->probKaon())){
                        if (mdcKalTrack->charge() > 0) nprotonp++;
                        if (mdcKalTrack->charge() < 0) nprotonm++;
                }
 
                if (pid->probElectron() > 1.0e-5){
                        if (mdcKalTrack->charge() > 0) nelectronp++;
                        if (mdcKalTrack->charge() < 0) nelectronm++;
                }
 
                if (pid->probMuon() > 1.0e-5){
                        if (mdcKalTrack->charge() > 0) nmuonp++;
                        if (mdcKalTrack->charge() < 0) nmuonm++;
                } 
        }
 
        // get etas, eta --> gamma gamma
        SmartDataPtr<EvtRecEtaToGGCol> evtRecEtaToGGCol(eventSvc(), EventModel::EvtRec::EvtRecEtaToGGCol);
        for (EvtRecEtaToGGCol::iterator iEta = evtRecEtaToGGCol->begin(); iEta != evtRecEtaToGGCol->end(); iEta++){
                if ((((*iEta)->chisq() < 2500) && ((*iEta)->unconMass() > 0.40) && ((*iEta)->unconMass() < 0.70))){
                        EvtRecTrack* lo = const_cast<EvtRecTrack*>((*iEta)->loEnGamma());
                        RecEmcShower* loShower = lo->emcShower();
                        if    ( (loShower->energy() > 0.025) && (loShower->time() > 0) && (loShower->time() < 14) && ( ( (abs(cos(loShower->theta())) < 0.80) && (loShower->energy() > 0.025) ) || ( (abs(cos(loShower->theta())) < 0.92) && (abs(cos(loShower->theta())) > 0.86) && (loShower->energy() > 0.050) ) ) ){
                                EvtRecTrack* hi = const_cast<EvtRecTrack*>((*iEta)->hiEnGamma());
                                RecEmcShower* hiShower = hi->emcShower();
                                if ( (hiShower->energy() > 0.025) && (hiShower->time() > 0) && (hiShower->time() < 14) && ( ( (abs(cos(hiShower->theta())) < 0.80) && (hiShower->energy() > 0.025) ) || ( (abs(cos(hiShower->theta())) < 0.92) && (abs(cos(hiShower->theta())) > 0.86) && (hiShower->energy() > 0.050) ) ) ){
                                        neta++;
                                }
                        }
                }
        }
 
        // get pi0s, pi0 --> gamma gamma
        SmartDataPtr<EvtRecPi0Col> evtRecPi0Col(eventSvc(), EventModel::EvtRec::EvtRecPi0Col);
        for (EvtRecPi0Col::iterator iPi0 = evtRecPi0Col->begin(); iPi0 != evtRecPi0Col->end(); iPi0++){
                if ((((*iPi0)->chisq() < 2500) && ((*iPi0)->unconMass() > 0.107) && ((*iPi0)->unconMass() < 0.163))){
                        EvtRecTrack* lo = const_cast<EvtRecTrack*>((*iPi0)->loEnGamma());
                        RecEmcShower* loShower = lo->emcShower();
                        if    ( (loShower->energy() > 0.025) && (loShower->time() > 0) && (loShower->time() < 14) && ( (( abs(cos(loShower->theta())) < 0.80) && (loShower->energy() > 0.025) ) || ( (abs(cos(loShower->theta())) < 0.92) && (abs(cos(loShower->theta())) > 0.86) && (loShower->energy() > 0.050) ) ) ){
                                EvtRecTrack* hi = const_cast<EvtRecTrack*>((*iPi0)->hiEnGamma());
                                RecEmcShower* hiShower = hi->emcShower();
                                if ( (hiShower->energy() > 0.025) && (hiShower->time() > 0) && (hiShower->time() < 14) && ( ( (abs(cos(hiShower->theta())) < 0.80) && (hiShower->energy() > 0.025) ) || ( (abs(cos(hiShower->theta())) < 0.92) && (abs(cos(hiShower->theta())) > 0.86) && (hiShower->energy() > 0.050) ) ) ){
                                        npi0++;
                                }
                        }
                }
        }
 
        // get kshorts
        // [[ vertexId is the pdgID ]]
        // Ks --> pi+ pi- and Lambda --> p+ pi- and ALambda --> p- pi+
        SmartDataPtr<EvtRecVeeVertexCol> evtRecVeeVertexCol(eventSvc(), EventModel::EvtRec::EvtRecVeeVertexCol);
        for (EvtRecVeeVertexCol::iterator iKs = evtRecVeeVertexCol->begin(); iKs != evtRecVeeVertexCol->end(); iKs++){
                if ((*iKs)->vertexId() == 310){
                        if ( ((*iKs)->mass() > 0.471) && ((*iKs)->mass() < 0.524) && ((*iKs)->chi2() < 100) ){
                                nks++;
                        }
                }
        }
 
        for (EvtRecVeeVertexCol::iterator iL = evtRecVeeVertexCol->begin(); iL != evtRecVeeVertexCol->end(); iL++){
                if ((*iL)->vertexId() == 3122){
                        if ( ((*iL)->mass() > 1.100) && ((*iL)->mass() < 1.130) && ((*iL)->chi2() < 100) ){
                                nlambda++;
                        }
                }
                if ((*iL)->vertexId() == -3122){
                        if ( ((*iL)->mass() > 1.100) && ((*iL)->mass() < 1.130) && ((*iL)->chi2() < 100) ){
                                nalambda++;
                        }
                }
        }
 
        unsigned int tagdata1 = nGoodCharged;
        unsigned int tagdata2 = IntToTag1(nNeutrk, nTottrk, ngamma, npi0);
        unsigned int tagdata3 = IntToTag2(npionp, npionm, nkaonp, nkaonm, nprotonp, nprotonm);
        unsigned int tagdata4 = IntToTag2(nlambda, nalambda, nelectronp, nelectronm, nmuonp, nmuonm);
        unsigned int tagdata5 = IntToTag2(nks, neta, nCharged, nGoodChargedp, nGoodChargedm, totCharged);
 
        m_tagFilterSvc->setTagData1(tagdata1);
        m_tagFilterSvc->setTagData2(tagdata2);
        m_tagFilterSvc->setTagData3(tagdata3);
        m_tagFilterSvc->setTagData4(tagdata4);
        m_tagFilterSvc->setTagData5(tagdata5);
 
        SmartDataPtr<EvtRecDTagCol> evtRecDTagCol(eventSvc(), EventModel::EvtRec::EvtRecDTagCol);
 
        std::set<UInt_t> Dsmode0, Dsmode1, D0mode0, D0mode1, Dmode0, Dmode1;
        Dsmode0.clear(); Dsmode1.clear(); D0mode0.clear(); D0mode1.clear(), Dmode0.clear(); Dmode1.clear();
        if ( evtRecDTagCol ) {
                EvtRecDTagCol::iterator iter_begin = evtRecDTagCol->begin();
                EvtRecDTagCol::iterator iter_end   = evtRecDTagCol->end();
                //cout << "Mode = ";
                for (EvtRecDTagCol::iterator iter = iter_begin; iter != iter_end; iter++) {
                        Int_t type = (*iter)->type();
                        Int_t mode = (*iter)->decayMode();
                        //cout << mode << "," << type << " " ;
                        if ( mode>=400 && mode<505 ) {
                                if (type==1) {
                                        Dsmode1.insert(mode);
                                } else {
                                        Dsmode0.insert(mode);
                                }
                        }
                        else if ( mode>=0 && mode<200 ) {
                                if (type==1) {
                                        D0mode1.insert(mode);
                                } else {
                                        D0mode0.insert(mode);
                                }
                        }
                        else if ( mode>=200 && mode<400 ) {
                                if (type==1) {
                                        Dmode1.insert(mode);
                                } else {
                                        Dmode0.insert(mode);
                                }
                        }
                }
                //cout << endl;
        }
 
        UInt_t Ds_modemap = 0, D0_modemap = 0, D_modemap = 0;
        int data0=3;
        std::set<UInt_t>::iterator it1, it0;
 
        if ( Dsmode1.size()>0 ) {
                std::cout << "Evt: " << iEvt << " Ds_1 : ";
                Ds_modemap |= (1<<31);
                for (it1=Dsmode1.begin(); it1!=Dsmode1.end(); ++it1) {
                        std::cout << ' ' << *it1;
                        Ds_modemap |= (1<<Ds_modeTag[*it1]);
                }
                std::cout << endl;
        } else {
                if ( Dsmode0.size()>0 ) {
                        std::cout << "Evt: " << iEvt << " Ds_0 : ";
                        for (it0=Dsmode0.begin(); it0!=Dsmode0.end(); ++it0) {
                                std::cout << ' ' << *it0;
                                Ds_modemap |= (1<<Ds_modeTag[*it0]);
                        }
                        std::cout << endl;
                }
        }
 
        if ( D0mode1.size()>0 ) {
                std::cout << "Evt: " << iEvt << " D0_1 : ";
                D0_modemap |= (1<<31);
                for (it1=D0mode1.begin(); it1!=D0mode1.end(); ++it1) {
                        std::cout << ' ' << *it1;
                        D0_modemap |= (1<<D0_modeTag[*it1]);
                }
                std::cout << endl;
        } else {
                if ( D0mode0.size()>0 ) {
                        std::cout << "Evt: " << iEvt << " D0_0 : ";
                        for (it0=D0mode0.begin(); it0!=D0mode0.end(); ++it0) {
                                std::cout << ' ' << *it0;
                                D0_modemap |= (1<<D0_modeTag[*it0]);
                        }
                        std::cout << endl;
                }
        }
 
        if ( Dmode1.size()>0 ) {
                std::cout << "Evt: " << iEvt << " D+_1 : ";
                D_modemap |= (1<<31);
                for (it1=Dmode1.begin(); it1!=Dmode1.end(); ++it1) {
                        std::cout << ' ' << *it1;
                        D_modemap |= (1<<D_modeTag[*it1]);
                }
                std::cout << endl;
        } else {
                if ( Dmode0.size()>0 ) {
                        std::cout << "Evt: " << iEvt << " D+_0 : ";
                        for (it0=Dmode0.begin(); it0!=Dmode0.end(); ++it0) {
                                std::cout << ' ' << *it0;
                                D_modemap |= (1<<D_modeTag[*it0]);
                        }
                        std::cout << endl;
                }
        }
 
        if ( Dsmode1.size()>0 ) {
                data0 = 0;
        } else if ( D0mode1.size()>0 || Dmode1.size()>0 ) {
                data0 = 1;
        } else if ( Dsmode0.size()>0 || D0mode0.size()>0 || Dmode0.size()>0 ) {
                data0 = 2;
        } 
        /*
             if (Ds_modemap>0) {
             cout << "Ds_mode = " << " ";
             for(int i = 0; i < 32; i ++) {
             int bit = Ds_modemap&(1<<(32-i-1));
             if(bit == 0) cout << 0;
             else cout << 1;
             if(i % 4 == 3)cout << ' ';
             }
             cout << " " << Ds_modemap << endl;
             }
             if (D0_modemap>0) {
             cout << "D0_mode = " << " ";
             for(int i = 0; i < 32; i ++) {
             int bit = D0_modemap&(1<<(32-i-1));
             if(bit == 0) cout << 0;
             else cout << 1;
             if(i % 4 == 3)cout << ' ';
             }
             cout << " " << D0_modemap << endl;
             }
             if (D_modemap>0) {
             cout << "Dp_mode = " << " ";
             for(int i = 0; i < 32; i ++) {
             int bit = D_modemap&(1<<(32-i-1));
             if(bit == 0) cout << 0;
             else cout << 1;
             if(i % 4 == 3)cout << ' ';
             }
             cout << " " << D_modemap << endl;
             }
             */
        m_tagFilterSvc->setTagData0(data0);
        m_tagFilterSvc->setTagData6(Ds_modemap);
        m_tagFilterSvc->setTagData7(D0_modemap);
        m_tagFilterSvc->setTagData8(D_modemap);
        //cout << iEvt << " Data0= " << tagdata1 << " 1= " << tagdata2 << " 2= " << tagdata3 << " 3= " << tagdata4 << " 4= " << tagdata5 << endl;
        //cout << "nNeutrk= " << nNeutrk << ", nTottrk= " << nTottrk << ", ngamma = "<< ngamma << ", npi0= " << npi0 << endl;
        //cout << "npionp=" << npionp << ", npionm=" << npionm << ", nkaonp= " << nkaonp << ", nkaonm= " <<  nkaonm << ", nprotonp= " <<  nprotonp << ", nprotonm= " <<  nprotonm << endl;
        //cout << "nlambda= " << nlambda << ", nalambda= " << nalambda << ", nelectronp= " <<  nelectronp << ", nelectronm= " <<  nelectronm << ", nmuonp= " <<  nmuonp << ", nmuonm= " <<  nmuonm << endl;
        //cout << "nks= " << nks << ", neta= " <<  neta << ", nCharged= " <<  nCharged << ", nGoodChargedp= " <<  nGoodChargedp << ", nGoodChargedm= " << nGoodChargedm << ", totCharged= " << totCharged << endl;
 
        iEvt++;
}

finalize()

StatusCode DTagSetAlg::finalize

(

)

Definition at line 494 of file DTagSetAlg.cxx.

                                {
        MsgStream log(msgSvc(), name());
        log << MSG::INFO << "In DTagSetAlg finalize()" << iEvt << "Events processed" << endmsg;  
        return StatusCode::SUCCESS;
}

initialize()

StatusCode DTagSetAlg::initialize

(

)

Definition at line 49 of file DTagSetAlg.cxx.

                                 {
        MsgStream log(msgSvc(), name());
 
        log << MSG::INFO << "in initialize()" << endmsg;
        iEvt=0; 
 
        Gaudi::svcLocator()->service("VertexDbSvc", m_vtxsvc);
        Gaudi::svcLocator()->service("TagFilterSvc", m_tagFilterSvc);
 
        Ds_modeTag[400] =  1;
        Ds_modeTag[401] =  0;
        Ds_modeTag[402] =  2;
        Ds_modeTag[403] = 13;
        Ds_modeTag[404] =  3;
        Ds_modeTag[405] =  4;
        Ds_modeTag[406] =  5;
        Ds_modeTag[407] = 14;
        Ds_modeTag[420] = 28;
        Ds_modeTag[421] =  6;
        Ds_modeTag[422] = 24;
        Ds_modeTag[423] = 25;
        Ds_modeTag[424] = 26;
        Ds_modeTag[425] = 27;
        Ds_modeTag[440] =  7;
        Ds_modeTag[441] =  8;
        Ds_modeTag[442] = 20;
        Ds_modeTag[450] = 21;
        Ds_modeTag[451] = 22;
        Ds_modeTag[452] = 23;
        Ds_modeTag[460] =  9;
        Ds_modeTag[461] = 15;
        Ds_modeTag[470] = 16;
        Ds_modeTag[471] = 17;
        Ds_modeTag[480] = 10;
        Ds_modeTag[481] = 11;
        Ds_modeTag[500] = 29;
        Ds_modeTag[501] = 18;
        Ds_modeTag[502] = 12;
        Ds_modeTag[503] = 19;
        Ds_modeTag[504] = 30;
 
        D0_modeTag[0]   =  0;
        D0_modeTag[1]   =  1;
        D0_modeTag[2]   =  2;
        D0_modeTag[3]   =  3;
        D0_modeTag[4]   =  4;
        D0_modeTag[5]   =  5;
        D0_modeTag[50]  =  6;
        D0_modeTag[51]  =  7;
        D0_modeTag[100] =  8;
        D0_modeTag[101] =  9;
        D0_modeTag[102] = 10;
        D0_modeTag[103] = 11;
        D0_modeTag[104] = 12;
        D0_modeTag[105] = 13;
        D0_modeTag[106] = 14;
        D0_modeTag[107] = 15;
        D0_modeTag[108] = 16;
        D0_modeTag[109] = 17;
        D0_modeTag[110] = 18;
        D0_modeTag[111] = 19;
        D0_modeTag[112] = 20;
        D0_modeTag[113] = 21;
        D0_modeTag[114] = 22;
        D0_modeTag[115] = 23;
        D0_modeTag[116] = 24;
        D0_modeTag[117] = 25;
        D0_modeTag[118] = 26;
        D0_modeTag[119] = 27;
        D0_modeTag[120] = 28;
        D0_modeTag[121] = 29;
        D0_modeTag[122] = 30;
        D0_modeTag[123] = 30;
        D0_modeTag[124] = 30;
        D0_modeTag[125] = 30;
        D0_modeTag[126] = 30;
 
        D_modeTag[200]  =  0;
        D_modeTag[201]  =  1;
        D_modeTag[202]  =  2;
        D_modeTag[203]  =  3;
        D_modeTag[204]  =  4;
        D_modeTag[205]  =  5;
        D_modeTag[206]  =  6;
        D_modeTag[207]  =  7;
        D_modeTag[208]  =  8;
        D_modeTag[209]  =  9;
        D_modeTag[210]  = 10;
        D_modeTag[211]  = 11;
        D_modeTag[212]  = 12;
        D_modeTag[213]  = 13;
        D_modeTag[214]  = 14;
        D_modeTag[215]  = 15;
        D_modeTag[216]  = 16;
        D_modeTag[217]  = 17;
        D_modeTag[218]  = 18;
        D_modeTag[219]  = 19;
        D_modeTag[220]  = 20;
        D_modeTag[221]  = 21;
        D_modeTag[222]  = 22;
 
        return StatusCode::SUCCESS;
}

IntToTag1()

unsigned int DTagSetAlg::IntToTag1	(	int	val1,
		int	val2,
		int	val3,
		int	val4
	)

Definition at line 500 of file DTagSetAlg.cxx.

{
        unsigned int res=0;
        res = (val4 +(val3<<8) +(val2<<16) +(val1<<24));
        return res;
}

Referenced by execute().

IntToTag2()

unsigned int DTagSetAlg::IntToTag2	(	int	val1,
		int	val2,
		int	val3,
		int	val4,
		int	val5,
		int	val6
	)

Definition at line 507 of file DTagSetAlg.cxx.

{
        unsigned int res=0;
        res = (val6 +(val5<<5) +(val4<<10) +(val3<<15) +(val2<<20) +(val1<<26));
        return res;
}

Referenced by execute().

---

The documentation for this class was generated from the following files:

• DTagSetAlg.h
• DTagSetAlg.cxx