CgemSegmentFitAlg Class Reference

#include <CgemSegmentFitAlg.h>

Inheritance diagram for CgemSegmentFitAlg:

Public Member Functions
	CgemSegmentFitAlg (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode	execute ()
StatusCode	initialize ()
StatusCode	finalize ()
StatusCode	exe_v1 ()
StatusCode	exe_v2 ()
	CgemSegmentFitAlg (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode	execute ()
StatusCode	initialize ()
StatusCode	finalize ()
StatusCode	exe_v1 ()
StatusCode	exe_v2 ()

Detailed Description

Definition at line 31 of file InstallArea/include/CgemSegmentFitAlg/CgemSegmentFitAlg/CgemSegmentFitAlg.h.

Constructor & Destructor Documentation

CgemSegmentFitAlg() [1/2]

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

Definition at line 64 of file CgemSegmentFitAlg.cxx.

                                                                                    :
    Algorithm(name, pSvcLocator)
{
    declareProperty("debug",          m_debug = 0);
    declareProperty("check",          m_checkIdx = 0);
    declareProperty("version",        m_ver = 2);
}

CgemSegmentFitAlg() [2/2]

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

Member Function Documentation

exe_v1() [1/2]

StatusCode CgemSegmentFitAlg::exe_v1

(

)

Definition at line 164 of file CgemSegmentFitAlg.cxx.

{
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "in excute()" << endreq;
 
    const double pi = 4.*atan(1);
    n_layer = 3;
 
    //   SmartDataPtr<CgemGeomSvc> cgemGeom(eventSvc(),  "CgemGeomSvc");
    //   if (!cgemGeom){
    //      log << MSG::WARNING << "Could not retrieve Cgem Geomtry list" << endreq;
    //      return StatusCode::FAILURE;
    //   }
    //   CgemGeomSvc::iterator itGeom = cgemGeom->begin();
 
    int vec_clusterid[3];
    m_seg_map.clear(); 
 
    SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
    if (!eventHeader) {
        log << MSG::FATAL << "Could not find Event Header" << endreq;
        return StatusCode::FAILURE;
    }
    cgem_evt = eventHeader->eventNumber();
    cgem_run = eventHeader->runNumber();
 
    //   if(cgem_evt<=221) return StatusCode::FAILURE;
    //   cout<<"event_ID = "<<cgem_evt<<endl;
    //if(cgem_evt%1000==0)cout<<"event_ID = "<<cgem_evt<<endl;
    //   if(cgem_evt > 62)return StatusCode::SUCCESS;
 
    SmartDataPtr<RecCgemSegmentCol> recCgemSegmentCol(eventSvc(),  "/Event/Recon/RecCgemSegmentCol");
    if (!recCgemSegmentCol){
        log << MSG::WARNING << "Could not retrieve Cgem segment list" << endreq;
        check_segNum++;
        //      return StatusCode::FAILURE;
        return StatusCode::SUCCESS;
    }
    RecCgemSegmentCol::iterator itSeg;
    /*
       itSeg = recCgemSegmentCol->begin();
       for(;itSeg != recCgemSegmentCol->end(); itSeg++){
       cout<<"event+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
       cout<<"event ID = "<<(*itSeg)->geteventid()<<endl;;
       cout<<"segment ID = "<<(*itSeg)->getsegmentid()<<endl;;
       cout<<"cluster ID_1 = "<<(*itSeg)->getclusterid_1()<<endl;;
       cout<<"cluster ID_2 = "<<(*itSeg)->getclusterid_2()<<endl;;
       cout<<"cluster ID_3 = "<<(*itSeg)->getclusterid_3()<<endl;;
       cout<<"event+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
       }
       */   
    SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(),  "/Event/Recon/RecCgemClusterCol");
    if (!recCgemClusterCol){
        log << MSG::WARNING << "Could not retrieve Cgem cluster list" << endreq;
        return StatusCode::FAILURE;
    }
 
    SmartDataPtr<Event::CgemMcHitCol> recCgemMcHitCol(eventSvc(), "/Event/MC/CgemMcHitCol");
    if (!recCgemMcHitCol)
    {
        log << MSG::WARNING << "Could not find event" << endreq;
        return StatusCode::FAILURE;
    }
 
    SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
    if (!mcParticleCol) {
        log << MSG::ERROR<< "Could not find McParticle" << endreq;
    } 
    else{
 
        Hep3Vector truth_p; HepPoint3D truth_position; double mc_charge(0.);
        McParticleCol::iterator iter_mc = mcParticleCol->begin();
        for (;iter_mc != mcParticleCol->end(); ++iter_mc ) {
            int pdg=(*iter_mc)->particleProperty();
            if(fabs(pdg)!=13) continue;
            if(truth_ntrack >=10)continue;
            if(pdg<0)mc_charge =  1.;
            else mc_charge =-1.; 
            double truth_x = (*iter_mc)->initialPosition().x();
            double truth_y = (*iter_mc)->initialPosition().y();
            double truth_z = (*iter_mc)->initialPosition().z();
 
            double truth_px = (*iter_mc)->initialFourMomentum().px();
            double truth_py = (*iter_mc)->initialFourMomentum().py();
            double truth_pz = (*iter_mc)->initialFourMomentum().pz();
            double truth_P  = sqrt(truth_px*truth_px+truth_py*truth_py+truth_pz*truth_pz);
            truth_CosTheta[truth_ntrack] = truth_pz/truth_P;
 
            truth_position = HepPoint3D(truth_x, truth_y, truth_z);
            truth_p = Hep3Vector(truth_px, truth_py, truth_pz);
 
            KalmanFit::Helix* tmp_helix = new KalmanFit::Helix(truth_position, truth_p, mc_charge);
            HepPoint3D tmp_pivot(0., 0., 0.);
            tmp_helix->pivot(tmp_pivot);
            truth_dr[truth_ntrack] = tmp_helix->dr();
            truth_dz[truth_ntrack] = tmp_helix->dz();
            truth_phi0[truth_ntrack] = tmp_helix->phi0();
            truth_kappa[truth_ntrack] = tmp_helix->kappa();
            truth_tanl[truth_ntrack] = tmp_helix->tanl();
            truth_charge[truth_ntrack] = mc_charge;
            if(m_debug){
                cout<<"Track Helix : "<<endl
                    <<"dr = "<<tmp_helix->dr()<<"\tphi0 = "<<tmp_helix->phi0()<<"\tkappa = "<<tmp_helix->kappa()
                    <<"\ndz = "<<tmp_helix->dz()<<"\ttanl = "<<tmp_helix->tanl()<<"\tpt = "<<1./tmp_helix->kappa()<<endl<<endl;
            }
 
            truth_ntrack++;
            delete tmp_helix;
        }
    }
    bool use_truth_hits = true;
    if(use_truth_hits){
        SmartDataPtr<Event::CgemMcHitCol> lv_CgemMcHitCol(eventSvc(), "/Event/MC/CgemMcHitCol");
        if (!lv_CgemMcHitCol)
        {
            log << MSG::WARNING << "Could not find event" << endreq;
            return StatusCode::FAILURE;
        }
        Event::CgemMcHitCol::iterator iter_truth = lv_CgemMcHitCol->begin();
        for( ; iter_truth != lv_CgemMcHitCol->end(); ++iter_truth){
            if(fabs((*iter_truth)->GetPDGCode())!=13) continue;
            double tmp_x = ((*iter_truth)->GetPositionXOfPrePoint()+(*iter_truth)->GetPositionXOfPostPoint())/2;
            double tmp_y = ((*iter_truth)->GetPositionYOfPrePoint()+(*iter_truth)->GetPositionYOfPostPoint())/2;
            double tmp_z = ((*iter_truth)->GetPositionZOfPrePoint()+(*iter_truth)->GetPositionZOfPostPoint())/2;
            double tmp_phi = atan2(tmp_y,tmp_x);
            while(tmp_phi < 0) tmp_phi += 2*pi;
            mc_phi[(*iter_truth)->GetLayerID()] = tmp_phi;
            double tmp_v = ((*iter_truth)->GetPositionZOfPostPoint()+(*iter_truth)->GetPositionZOfPrePoint())/2./10.*sin(CgemSegmentFun::a_stero[(*iter_truth)->GetLayerID()])*CgemSegmentFun::r_layer[(*iter_truth)->GetLayerID()]/CgemSegmentFun::rl[(*iter_truth)->GetLayerID()]+CgemSegmentFun::r_layer[(*iter_truth)->GetLayerID()]*tmp_phi*cos(CgemSegmentFun::a_stero[(*iter_truth)->GetLayerID()]);
            mc_v[(*iter_truth)->GetLayerID()] = tmp_v;
            CgemSegmentFun::vec_z[(*iter_truth)->GetLayerID()] = ((*iter_truth)->GetPositionZOfPostPoint()+(*iter_truth)->GetPositionZOfPrePoint())/2./10.;//cm
            CgemSegmentFun::vec_phi[(*iter_truth)->GetLayerID()] = tmp_phi;
            //   cout<<"truth_z = "<<tmp_z<<"\ttruth_v = "<<tmp_v<<"\ttruth_phi = "<<tmp_phi<<endl;
            //   cout<<"truth_z = "<<tmp_z<<"\ttruth_v = "<<tmp_v<<"\ttruth_phi = "<<tmp_phi<<endl;
        }
        double Ini_d0(0.), Ini_phi0(0.), Ini_kappa(0.), Ini_z0(0.), Ini_tanl(0.);
        CgemSegmentFun::GetHelixVarBeforeFit( -1, Ini_d0, Ini_phi0, Ini_kappa, Ini_z0, Ini_tanl);
    }
 
    itSeg = recCgemSegmentCol->begin();
    cgem_nsegment=0;
    for(;itSeg != recCgemSegmentCol->end(); ++itSeg){
        if(cgem_nsegment >= 1000){check_segOver++;continue;}
        cgem_segmentid[cgem_nsegment] = (*itSeg)->getsegmentid();
        vec_clusterid[0] = (*itSeg)->getclusterid_1();
        vec_clusterid[1] = (*itSeg)->getclusterid_2();
        vec_clusterid[2] = (*itSeg)->getclusterid_3();
        int IsTruth = (*itSeg)->getmatch();
 
        // initialize minuit
        Int_t ierflg;
        Double_t arglist[10];
        m_failed = 0;
 
        m_pMnTrk = new TMinuit(5);
        m_pMnTrk -> SetPrintLevel(-1);
        m_pMnTrk -> SetFCN(CgemSegmentFun::fcnTrk);
        m_pMnTrk -> SetErrorDef(1.0); // For Chisq
        // m_pMnTrk -> SetErrorDef(0.5); // For likelihood 
        m_pMnTrk -> mnparm(0, "dr"   , 0, 0.1 , -10, 10, ierflg);//Set starting values and step sizes for parameters
        m_pMnTrk -> mnparm(1, "phi0" , 0, 0.16, 0, 6.3, ierflg);
        m_pMnTrk -> mnparm(2, "kappa", 0, 0.4, -100, 100, ierflg);
        m_pMnTrk -> mnparm(3, "dz"   , 0, 0.5 , -50, 50, ierflg);
        m_pMnTrk -> mnparm(4, "tanL" , 0, 0.005 , -10, 10, ierflg);
        //m_pMnTrk -> mnparm(0, "dr"   , 0, 0.09 , -10, 10, ierflg);//Set starting values and step sizes for parameters
        //m_pMnTrk -> mnparm(1, "phi0" , 0, 0.016, 0, 6.3, ierflg);
        //m_pMnTrk -> mnparm(2, "kappa", 0, 0.01, -30, 30, ierflg);
        //m_pMnTrk -> mnparm(3, "dz"   , 0, 0.05 , -50, 50, ierflg);
        //m_pMnTrk -> mnparm(4, "tanL" , 0, 0.004 , -10, 10, ierflg);
        arglist[0] = 0;
        m_pMnTrk -> mnexcm("Set NOW", arglist, 0, ierflg);
 
 
        RecCgemClusterCol::iterator itTrk=recCgemClusterCol->begin();
        for(; itTrk != recCgemClusterCol->end(); ++itTrk){
            if((*itTrk)->getflag() == 0){
                n_stripX[(*itTrk)->getlayerid()] = (*itTrk)->getclusterflage()-(*itTrk)->getclusterflagb()+1;
                CgemSegmentFun::nStrip_X[(*itTrk)->getlayerid()] = (*itTrk)->getclusterflage()-(*itTrk)->getclusterflagb()+1;//cm
            }
            if((*itTrk)->getflag() == 1){
                n_stripV[(*itTrk)->getlayerid()] = (*itTrk)->getclusterflage()-(*itTrk)->getclusterflagb()+1;
                CgemSegmentFun::nStrip_V[(*itTrk)->getlayerid()] = (*itTrk)->getclusterflage()-(*itTrk)->getclusterflagb()+1;
            }
            if((*itTrk)->getflag() != 2) continue;
            if(!CheckClusterID((*itTrk)->getlayerid(), (*itTrk)->getclusterid(), vec_clusterid)) continue;
            //CgemSegmentFun::sheet_ID[(*itTrk)->getlayerid()] = (*itTrk)->getsheetid();
            double tmp_phi = (*itTrk)->getrecphi();
 
            CgemSegmentFun::vec_z[(*itTrk)->getlayerid()] = (*itTrk)->getRecZ()/10.;//cm
            CgemSegmentFun::vec_phi[(*itTrk)->getlayerid()] = (*itTrk)->getrecphi();
            //CgemSegmentFun::vec_v[(*itTrk)->getlayerid()] = (*itTrk)->getrecv()/10.;//cm
            rec_phi[(*itTrk)->getlayerid()] = (*itTrk)->getrecphi();
            rec_v[(*itTrk)->getlayerid()] = (*itTrk)->getRecZ()/10.*sin(CgemSegmentFun::a_stero[(*itTrk)->getlayerid()])*CgemSegmentFun::r_layer[(*itTrk)->getlayerid()]/CgemSegmentFun::rl[(*itTrk)->getlayerid()]+CgemSegmentFun::r_layer[(*itTrk)->getlayerid()]*tmp_phi*cos(CgemSegmentFun::a_stero[(*itTrk)->getlayerid()]);//cm
            CgemSegmentFun::vec_v[(*itTrk)->getlayerid()] = rec_v[(*itTrk)->getlayerid()];//cm
            //cout<<"vec_Z = "<<(*itTrk)->getRecZ()/10.<<"\tvec_phi = "<<(*itTrk)->getrecphi()<<endl;
            //cout<<"layerID = "<<(*itTrk)->getlayerid()<<"\tX_n_strip = "<<n_stripX[(*itTrk)->getlayerid()] <<"\tV_n_strip = "<<n_stripV[(*itTrk)->getlayerid()] <<endl;
            //cout<<"vec_v["<<(*itTrk)->getlayerid()<<"] = "<<(*itTrk)->getrecv()/10.<<"\t"<<rec_v[(*itTrk)->getlayerid()]<<endl;
            //cout<<"sheetID["<<(*itTrk)->getlayerid()<<"] = "<<(*itTrk)->getsheetid()<<endl;
        }
 
        double Ini_d0(0.), Ini_phi0(0.), Ini_kappa(0.), Ini_z0(0.), Ini_tanl(0.);
        double chisq(0.);
 
        int charge = GetChargeOfTrack();
        cgem_charge[cgem_nsegment]  = charge;
        CgemSegmentFun::GetHelixVarBeforeFit( charge, Ini_d0, Ini_phi0, Ini_kappa, Ini_z0, Ini_tanl);
 
        double ErrMatrix[5][5] ;
        double TrkPar[5] = {Ini_d0, Ini_phi0, Ini_kappa, Ini_z0, Ini_tanl};
        //double TrkPar[5] = {truth_dr[0], truth_phi0[0], truth_kappa[0], truth_dz[0], truth_tanl[0]};
        double TrkParErr[5] ;
        double helixErr[15];
        int fit = mnTrkFit(TrkPar, TrkParErr, &ErrMatrix[0][0], chisq);
        int fitAgain = 0; double kappa_err;
        if(!fit){
            //   cout<< "fit failed"<<"\tevent ID = "<<cgem_evt<<"\trun = "<<cgem_run << endl;
            m_failed = 1;
            double tmp_recphi[3], tmp_kappa[8];
            for(int ii = 0, mm = 0; ii < 2; ii++)
                for(int jj = 0; jj < 2; jj++)
                    for(int kk = 0; kk < 2; kk++, mm++){
                        tmp_recphi[0] = CgemSegmentFun::vec_phi[0] + pow(-1.,ii)*(wide_delta*n_stripX[0]+pitch*(n_stripX[0]-1))/CgemSegmentFun::r_layer[0];
                        tmp_recphi[1] = CgemSegmentFun::vec_phi[1] + pow(-1.,jj)*(wide_delta*n_stripX[1]+pitch*(n_stripX[1]-1))/CgemSegmentFun::r_layer[1];
                        tmp_recphi[2] = CgemSegmentFun::vec_phi[2] + pow(-1.,kk)*(wide_delta*n_stripX[2]+pitch*(n_stripX[2]-1))/CgemSegmentFun::r_layer[2];
                        charge = GetChargeOfTrack(tmp_recphi);
                        tmp_kappa[mm] = CgemSegmentFun::GetKappaAfterFit(charge, tmp_recphi);
                        //cout<<"tmp kappa["<<mm<<"] = "<<tmp_kappa[mm]<<endl;
                    }
            double kappa_s = GetMin(tmp_kappa);
            double kappa_l = GetMax(tmp_kappa);
            kappa_err = (kappa_s - kappa_l)/2.;
            double TrkPar[5] = {Ini_d0, Ini_phi0, (kappa_l+kappa_s)/2., Ini_z0, Ini_tanl};
            TrkPar[2] = kappa_l;
            fitAgain = mnTrkFit(TrkPar, TrkParErr, &ErrMatrix[0][0], chisq);
 
        }
        //cout<<"m_failed = "<<m_failed<<endl;
 
        //cout<<"charge = "<<charge<<"\t"<<fit<<"\t"<<fitAgain<<endl<<endl;
        //      cout<<
        //   "delta_dr = "<<Ini_d0-truth_dr[0]<<endl<<
        //   "delta_phi0 = "<<Ini_phi0-truth_phi0[0]<<endl<<
        //   "delta_kappa = "<<Ini_kappa-truth_kappa[0]<<endl<<
        //   "delta_dz = "<<Ini_z0-truth_dz[0]<<endl<<
        //   "delta_tanl = "<<Ini_tanl-truth_tanl[0]<<endl<<endl;
 
        //RecCgemSegment* recsegment = new RecCgemSegment();
        cgem_fitFlag[cgem_nsegment]  = fit;
        if(fit==1||fitAgain){
            if(fit){
                for (int ie = 0, k = 0 ; ie < 5 ; ie++){
                    for (int je = ie ; je < 5 ; je ++, k++) {
                        helixErr[k] = ErrMatrix[ie][je];
                    }
                }  
            }
 
            if(fitAgain){
                for (int ie = 0, k = 0 ; ie < 5 ; ie++){
                    for (int je = ie ; je < 5 ; je ++, k++) {
                        if(ie==2&&je<2)helixErr[k] = ErrMatrix[ie+2][je];         //k = 2,6
                        else if(je<2&&ie>2) helixErr[k] = ErrMatrix[ie-1][je];  //k = 3,4,7,8
                        else if(je>2&&ie>2) helixErr[k] = ErrMatrix[ie-1][je-1]; //k = 12,13,14
                        else if(je==2) helixErr[k] = 0.; //k = 9,10,11
                        else helixErr[k] = ErrMatrix[ie][je]; // k = 0,1,5
 
                        if(k==9)helixErr[k] = pow(kappa_err,2);
                        // cout<<"ErrMatrix["<<ie<<"]["<<je<<"] = "<<ErrMatrix[ie][je]<<"\t"<<helixErr[k]<<endl;
                    }
                }
                TrkParErr[0] = 0.8812/0.0677*TrkParErr[0];
                TrkParErr[1] = 0.266/0.004*TrkParErr[1];
                TrkParErr[2] = kappa_err;
                helixErr[0] = 0.8812/0.0677*helixErr[0];
                helixErr[5] = 0.266/0.004*helixErr[5];
            }
            //   cout<<"helixErr[0] = "<<helixErr[0]<<"\tTrkParErr[0] = "<<TrkParErr[0]<<"\n"<<
            //         "helixErr[5] = "<<helixErr[5]<<"\tTrkParErr[1] = "<<TrkParErr[1]<<"\n"<<
            //         "helixErr[9] = "<<helixErr[9]<<"\tTrkParErr[2] = "<<TrkParErr[2]<<"\n"<<
            //         "helixErr[12] = "<<helixErr[12]<<"\tTrkParErr[3] = "<<TrkParErr[3]<<"\n"<<
            //         "helixErr[14] = "<<helixErr[14]<<"\tTrkParErr[4] = "<<TrkParErr[4]<<"\n"<<endl<<endl;
            //recsegment->sethelix(TrkPar);
            //recsegment->sethelix_err(helixErr);
            //recsegment->setmatch(IsTruth);
            (*itSeg)->sethelix(TrkPar);
            (*itSeg)->sethelix_err(helixErr);
            (*itSeg)->setmatch(IsTruth);
            cgem_dr[cgem_nsegment]    = TrkPar[0];
            cgem_phi0[cgem_nsegment]  = TrkPar[1];
            cgem_kappa[cgem_nsegment] = TrkPar[2];
            cgem_dz[cgem_nsegment]    = TrkPar[3];
            cgem_tanl[cgem_nsegment]  = TrkPar[4];
            cgem_chisq[cgem_nsegment]  = chisq;
            if(m_debug){
                cout<< "fit success" << endl
                    <<"dr = "  <<TrkPar[0]<<" +- "<<TrkParErr[0]<<"\t\t"
                    <<"phi0 = "<<TrkPar[1]<<" +- "<<TrkParErr[1]<<"\t\t"
                    <<"kappa= "<<TrkPar[2]<<" +- "<<TrkParErr[2]<<"\n"
                    <<"dz = "  <<TrkPar[3]<<" +- "<<TrkParErr[3]<<"\t\t"
                    <<"tanl = "<<TrkPar[4]<<" +- "<<TrkParErr[4]<<"\n\n"
                    <<"kappa_err= "<<TrkParErr[2]<<endl;
                cout<<"chisq = "<<chisq<<endl;
            }
        }
        else{
            //cout<< "fit failed"<<"\tevent ID = "<<cgem_evt<<"\trun = "<<cgem_run << endl;
 
            double fail_helix[5], fail_helix_err[15] ;
            for(int i = 0; i <15; i++){
                if(i<5)fail_helix[i] = -999.;
                fail_helix_err[i] = -999.;
            }
            (*itSeg)->sethelix(fail_helix);
            (*itSeg)->sethelix_err(fail_helix_err);
            (*itSeg)->setmatch(IsTruth);
            cgem_dr[cgem_nsegment]    = -999.;
            cgem_phi0[cgem_nsegment]  = -999.;
            cgem_kappa[cgem_nsegment] = -999.;
            cgem_dz[cgem_nsegment]    = -999.;
            cgem_tanl[cgem_nsegment]  = -999.;
            cgem_chisq[cgem_nsegment]  = -999.;
        }
        //m_seg_map.insert(pair<int,RecCgemSegment*>(cgem_nsegment,recsegment));  
 
        cgem_nsegment++;
 
        //cout<<"***************************#######*******************************"<<endl<<endl;
        delete m_pMnTrk;
        //delete recsegment;
    }
    /*
       RecCgemSegmentCol* segmentcol = new RecCgemSegmentCol;
 
       IDataProviderSvc* evtSvc = NULL;
       Gaudi::svcLocator()->service("EventDataSvc",evtSvc);
       if (evtSvc) {
       log << MSG::INFO << "makeTds:event Svc has been found" << endreq;
       } else {
       log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
       return StatusCode::SUCCESS;
       }
 
       StatusCode segmentsc;
       IDataManagerSvc *dataManSvc;
       dataManSvc= dynamic_cast<IDataManagerSvc*>(evtSvc);
       DataObject *aSegmentCol;
       evtSvc->findObject("/Event/Recon/RecCgemSegmentCol",aSegmentCol);
       if(aSegmentCol != NULL) {
       dataManSvc->clearSubTree("/Event/Recon/RecCgemSegmentCol");
       evtSvc->unregisterObject("/Event/Recon/RecCgemSegmentCol");
       }
 
 
    //push back the segment from the map//
    for(m_seg_map_it = m_seg_map.begin();m_seg_map_it!=m_seg_map.end();++m_seg_map_it){
    RecCgemSegment* recsegment = (*m_seg_map_it).second;
    segmentcol->push_back(recsegment);
    cout<<">>>>>>> push recsegment, clusterid_1="<<recsegment->getclusterid_1()<<"kappa: "<<recsegment->gethelix(2)<<endl;
 
    }
 
    segmentsc = evtSvc->registerObject("/Event/Recon/RecCgemSegmentCol", segmentcol);
    if( segmentsc.isFailure() ) {
    log << MSG::FATAL << "Could not register RecCgemSegment" << endreq;
    return StatusCode::SUCCESS;
    }
    log << MSG::INFO << "RecCgemSegmentCol registered successfully!" <<endreq;
    */
 
 
 
    helix_ntuple->write();
    /*
       SmartDataPtr<RecCgemSegmentCol> recCgemSegmentCol2(eventSvc(),  "/Event/Recon/RecCgemSegmentCol");
       itSeg = recCgemSegmentCol2->begin();
       for(;itSeg != recCgemSegmentCol2->end(); itSeg++){
       cout<<"event+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
       cout<<"event ID = "<<(*itSeg)->geteventid()<<endl;;
       cout<<"segment ID = "<<(*itSeg)->getsegmentid()<<endl;;
       cout<<"cluster ID_1 = "<<(*itSeg)->getclusterid_1()<<endl;;
       cout<<"cluster ID_2 = "<<(*itSeg)->getclusterid_2()<<endl;;
       cout<<"cluster ID_3 = "<<(*itSeg)->getclusterid_3()<<endl;;
       cout<<"event+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
       }
       */
 
 
    return StatusCode::SUCCESS;
}

Referenced by execute().

exe_v1() [2/2]

StatusCode CgemSegmentFitAlg::exe_v1

(

)

exe_v2() [1/2]

StatusCode CgemSegmentFitAlg::exe_v2

(

)

Definition at line 553 of file CgemSegmentFitAlg.cxx.

{
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "in excute()" << endreq;
    
    SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
    if (!eventHeader) {
        log << MSG::FATAL << "Could not find Event Header" << endreq;
        return StatusCode::FAILURE;
    }
    m_run = eventHeader->runNumber();
    m_evt = eventHeader->eventNumber();
    getMcPar();
 
    SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(),  "/Event/Recon/RecCgemClusterCol");
    if (!recCgemClusterCol){
        log << MSG::WARNING << "Could not retrieve Cgem cluster list" << endreq;
        return StatusCode::FAILURE;
    }
    RecCgemClusterCol::iterator itCluBegin=recCgemClusterCol->begin();
 
    SmartDataPtr<RecCgemSegmentCol> recCgemSegmentCol(eventSvc(),  "/Event/Recon/RecCgemSegmentCol");
    if (!recCgemSegmentCol){
        log << MSG::WARNING << "Could not retrieve Cgem segment list" << endreq;
        check_segNum++;
        //      return StatusCode::FAILURE;
        return StatusCode::SUCCESS;
    }
    if(m_checkIdx) {
        cgem_run = m_run;
        cgem_evt = m_evt;
        cgem_nsegment=0;
    }
    int cluId[3];
    RecCgemSegmentCol::iterator itSeg = recCgemSegmentCol->begin();
    for( ;itSeg != recCgemSegmentCol->end(); ++itSeg)
    {
        cluId[0] = (*itSeg)->getclusterid_1();
        cluId[1] = (*itSeg)->getclusterid_2();
        cluId[2] = (*itSeg)->getclusterid_3();
        if(m_debug) 
        {
            cout<<"input clusters: "<<endl;
        }
        double vec_Q[3], vec_T[3];
        for(int i=0; i<3; i++)
        {
            RecCgemClusterCol::iterator itClu = itCluBegin+cluId[i];
            CgemSegmentFun::vec_phi[i] = (*itClu)->getrecphi();// in radian
            CgemSegmentFun::vec_z[i]   = (*itClu)->getRecZ()/10.;  // in cm
            vec_Q[i]=(*itClu)->getenergydeposit();
            vec_T[i]=100;
            if(m_debug) 
            {
                cout    <<"layer "<<i<<", cluster id = "<<cluId[i]<<", itCluId = "<<(*itClu)->getclusterid()
                    <<", phi="<<CgemSegmentFun::vec_phi[i]<<", z="<<CgemSegmentFun::vec_z[i]<<endl;
            }
        }
        //cout<<"charge = "<<GetChargeOfTrack()<<", "<<estiChargeOfSeg()<<endl;
        int charge = estiChargeOfSeg();
 
        // initial Helix parameters
        double dr_ini, phi0_ini, kappa_ini, dZ_ini, tanL_ini;
        CgemSegmentFun::GetHelixVarBeforeFit(charge, dr_ini, phi0_ini, kappa_ini, dZ_ini, tanL_ini);
        double dr_ini0 = dr_ini;
        double phi0_ini0 = phi0_ini; 
        double kappa_ini0 = kappa_ini;
            double dZ_ini0 = dZ_ini;
            double tanL_ini0 = tanL_ini;
 
        double helix[5] = {dr_ini, phi0_ini, kappa_ini, dZ_ini, tanL_ini};
        double helixErr[15];
        HepPoint3D pivot(0., 0., 0.);
        HepVector vec_trkpar(5);
        for(int i = 0; i < 5; i++){
            vec_trkpar[i] = helix[i];
        }
        //KalmanFit::Helix* cgem_helix = new KalmanFit::Helix(pivot, vec_trkpar);
        KalmanFit::Helix cgem_helix(pivot, vec_trkpar);
 
        /* initialize error matrix */
        int iView(0), mode(2);
        for(int i=0; i<3; i++)
        {
            double dphi = CgemSegmentFun::IntersectCylinder(CgemSegmentFun::rl[i], cgem_helix);
            HepPoint3D pos = cgem_helix.x(dphi);
            double phi_pos = pos.phi();
            Hep3Vector p3 = cgem_helix.momentum(dphi);
            double phi_p = p3.phi();
            double dPhi = phi_p-phi_pos;
            while(dPhi>CLHEP::pi)  dPhi-=CLHEP::twopi;
            while(dPhi<-CLHEP::pi) dPhi+=CLHEP::twopi;
            CgemSegmentFun::x_reso[i]= (myCgemCalibSvc->getSigma(i,0,mode,dPhi,vec_Q[i],vec_T[i]))/10.;// mm -> cm
            CgemSegmentFun::v_reso[i]= (myCgemCalibSvc->getSigma(i,1,mode,dPhi,vec_Q[i],vec_T[i]))/10.;// mm -> cm
        }
 
        // initialize minuit
        Int_t ierflg, istat;
        Double_t arglist[10];
        m_failed = 0;
 
        m_pMnTrk = new TMinuit(5);// num of par
        m_pMnTrk -> SetPrintLevel(-1);// no print out
        m_pMnTrk -> mnexcm("Set NOW",  arglist, 1, ierflg);// no warning
        m_pMnTrk -> SetFCN(CgemSegmentFun::fcnTrk);
        m_pMnTrk -> SetErrorDef(1.0); // For Chisq
        // m_pMnTrk -> SetErrorDef(0.5); // For likelihood 
        bool fitSucceed = false;
        bool fitUsable = false;
        int istat_last = 0;
        int nFitTried = 0; int iFitStored = 0;
        double chi2_min = 999;
        double chi2_saved = 999;
        double chi2_fixKapa_min = 990;
        while((nFitTried<3&&!fitSucceed)||(nFitTried<5&&!fitUsable)) {
            m_pMnTrk-> mnparm(0, "dr"   , dr_ini,    0.001 ,               dr_ini-10,     dr_ini+10,    ierflg);//Set starting values and step sizes for parameters
            m_pMnTrk-> mnparm(1, "phi0" , phi0_ini,  0.001,                phi0_ini-1.0,  phi0_ini+1.0, ierflg);
            m_pMnTrk-> mnparm(2, "kappa", kappa_ini, max(fabs(kappa_ini)*0.01,0.01), kappa_ini-5.0, kappa_ini+5.0,ierflg);
            if(nFitTried>=3) m_pMnTrk->FixParameter(2);
            m_pMnTrk -> mnparm(3, "dz"   , dZ_ini,    0.01 ,                dZ_ini-20.0,   dZ_ini+20.0,  ierflg);
            m_pMnTrk -> mnparm(4, "tanL" , tanL_ini,  0.001 ,               tanL_ini-1.0,  tanL_ini+1.0, ierflg);
            arglist[0] = 2000; //arglist[1] = 0.1;
            m_pMnTrk -> mnexcm("MIGRAD", arglist, 1, ierflg);
            nFitTried++;
            double trkpar[5], trkparErr[5], emat[25], chisq;
            Double_t fmin, edm, errdef;
            Int_t nvpar, nparx;
            m_pMnTrk -> mnstat(fmin, edm, errdef, nvpar, nparx, istat);
            for(int i=0; i<5; i++){
                m_pMnTrk -> GetParameter(i, trkpar[i], trkparErr[i]);
            }
            // error matrix
            m_pMnTrk -> mnemat(emat, 5);
            chisq = fmin;
            
            //HepVector aNew = CgemSegmentFun::CalculateHelix(vec_trkpar, pivot);
 
            if(m_debug) {
                cout<<"Fit ierflg="<<ierflg<<", istat="<<istat<<endl;
                cout.precision(4);
                cout<<setw(10)<<" "<<setw(10)<<"dr"<<setw(10)<<"phi0"<<setw(10)<<"kappa"<<setw(10)<<"dz"<<setw(10)<<"tanL"<<endl;
                cout<<setw(10)<<"MC truth"<<setw(10)<<m_helix_mc[0]<<setw(10)<<m_helix_mc[1]<<setw(10)<<m_helix_mc[2]<<setw(10)<<m_helix_mc[3]<<setw(10)<<m_helix_mc[4]<<endl;
                cout<<setw(10)<<"ini"<<setw(10)<<dr_ini<<setw(10)<<phi0_ini<<setw(10)<<kappa_ini<<setw(10)<<dZ_ini<<setw(10)<<tanL_ini<<endl;
                //cout<<setw(10)<<"calcu"<<setw(10)<<aNew[0]<<setw(10)<<aNew[1]<<setw(10)<<aNew[2]<<setw(10)<<aNew[3]<<setw(10)<<aNew[4]<<"chi2="<<chisq<<endl;
                cout<<setw(10)<<"fit"<<setw(10)<<trkpar[0]<<setw(10)<<trkpar[1]<<setw(10)<<trkpar[2]<<setw(10)<<trkpar[3]<<setw(10)<<trkpar[4]<<"chi2="<<chisq<<endl;
                cout<<setw(10)<<"fit err"<<setw(10)<<trkparErr[0]<<setw(10)<<trkparErr[1]<<setw(10)<<trkparErr[2]<<setw(10)<<trkparErr[3]<<setw(10)<<trkparErr[4]<<endl;
                cout.precision(6);
            }
 
            /*
             *-*  IERFLG is now (94.5) defined the same as ICONDN in MNCOMD
             *-*            = 0: command executed normally
             *-*              1: command is blank, ignored
             *-*              2: command line unreadable, ignored
             *-*              3: unknown command, ignored
             *-*              4: abnormal termination (e.g., MIGRAD not converged)
             *-*              9: reserved
             *-*             10: END command
             *-*             11: EXIT or STOP command
             *-*             12: RETURN command
             *
             *    ISTAT: a status integer indicating how good is the covariance
             *-*           matrix:  0= not calculated at all
             *-*                    1= approximation only, not accurate
             *-*                    2= full matrix, but forced positive-definite
             *-*                    3= full accurate covariance matrix
             */
            if( (0==ierflg) && (istat==3) ){
                fitSucceed=true;
            }
            else {
                //cout<<"TMinuit failed"<<endl;
                if(m_debug) {
                    cout<<"TMinuit tried "<<nFitTried<<" times!"<<endl;
                }
                dr_ini=trkpar[0];
                phi0_ini=trkpar[1];
                    //kappa_ini=trkpar[2]+pow(-1,nFitTried);
                    kappa_ini=trkpar[2]+0.5*pow(-1,nFitTried);
                    dZ_ini=trkpar[3]; 
                tanL_ini=trkpar[4];
            }
            if(nFitTried>=3) kappa_ini=0.5*pow(-1,nFitTried);
            
            if( (0==ierflg) && (istat>0) ){
                if(nFitTried<4) fitUsable = true;
                if((nFitTried<4&&istat>istat_last)||(nFitTried>3&&chisq<chi2_fixKapa_min)) {
                    int k=0;
                    chi2_fixKapa_min = chisq;
                    chi2_saved = chisq;
                    for(int ii=0; ii<5; ii++)
                    {
                        helix[ii]=trkpar[ii];
                        for(int jj=ii; jj<5; jj++)
                        {
                            helixErr[k]=emat[ii*5+jj];
                            k++;
                        }
                    }
                    istat_last = istat;
                    iFitStored = nFitTried;
 
                }
            }
        }
        //if(!fitUsable) cout<<"TMinuit failed after "<<nFitTried<<" attemps!"<<endl;
        if(m_debug) 
            cout<<"Fit results "<<iFitStored<<" are saved!"<<endl;
 
        if(m_checkIdx&&cgem_nsegment<1000) {
            cgem_segmentid[cgem_nsegment]=(*itSeg)->getsegmentid();
            cgem_dr_ini[cgem_nsegment]=dr_ini0;
            cgem_dz_ini[cgem_nsegment]=dZ_ini0;
            cgem_phi0_ini[cgem_nsegment]=phi0_ini0;
            cgem_kappa_ini[cgem_nsegment]=kappa_ini0;
            cgem_tanl_ini[cgem_nsegment]=tanL_ini0;
            cgem_dr[cgem_nsegment]=helix[0];
            cgem_dz[cgem_nsegment]=helix[3];
            cgem_phi0[cgem_nsegment]=helix[1];
            cgem_kappa[cgem_nsegment]=helix[2];
            cgem_tanl[cgem_nsegment]=helix[4];
            cgem_chisq[cgem_nsegment]=chi2_saved;
            cgem_fitFlag[cgem_nsegment]=iFitStored;
            cgem_charge[cgem_nsegment]=charge;
            cgem_nsegment++;
        }
        
        // fill out
        (*itSeg)->sethelix(helix);
        (*itSeg)->sethelix_err(helixErr);
 
        delete m_pMnTrk;
 
    }// loop CGEM segments
    //cout<<"end exe_v2()"<<endl;
    if(m_checkIdx) helix_ntuple->write();
 
    return StatusCode::SUCCESS;
 
}

Referenced by execute().

exe_v2() [2/2]

StatusCode CgemSegmentFitAlg::exe_v2

(

)

execute() [1/2]

StatusCode CgemSegmentFitAlg::execute

(

)

Definition at line 156 of file CgemSegmentFitAlg.cxx.

{
    if(m_ver==1) exe_v1();
    else if(m_ver==2) exe_v2();
 
    return StatusCode::SUCCESS;
}

execute() [2/2]

StatusCode CgemSegmentFitAlg::execute

(

)

finalize() [1/2]

StatusCode CgemSegmentFitAlg::finalize

(

)

Definition at line 1000 of file CgemSegmentFitAlg.cxx.

                                       {
 
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "in finalize()" << endreq;
    //cout<<"leave segment number by dZ dPhi = "<<check_segNum<<endl
    //   <<"leave segment number over 1000  = "<<check_segOver<<endl;
 
    //cout<<endl<<"******* CgemSegmentFitAlg *******"<<endl;
 
    return StatusCode::SUCCESS;
}

finalize() [2/2]

StatusCode CgemSegmentFitAlg::finalize

(

)

initialize() [1/2]

StatusCode CgemSegmentFitAlg::initialize

(

)

Definition at line 73 of file CgemSegmentFitAlg.cxx.

                                        {
 
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "in initialize()" << endreq;
    if(m_checkIdx>0) {
        /////////////////////////////// book ntuple //////////////////////////////////////
        INTupleSvc* ntupleSvc;
        StatusCode sc = service("NTupleSvc", ntupleSvc);
        NTuplePtr mt(ntupleSvc, "FILE169/CgemSegmentHelix");
        if ( mt ) helix_ntuple = mt;
        else{
            helix_ntuple = ntupleSvc->book("FILE169/CgemSegmentHelix", CLID_ColumnWiseTuple, "SegFit");
            StatusCode status;
            if(helix_ntuple){
                status = helix_ntuple->addItem("run", cgem_run);
                status = helix_ntuple->addItem("evt", cgem_evt);
                status = helix_ntuple->addItem("ntrack",truth_ntrack,0,10);
                status = helix_ntuple->addIndexedItem("truth_dr",        truth_ntrack ,  truth_dr);
                status = helix_ntuple->addIndexedItem("truth_dz",        truth_ntrack ,  truth_dz);
                status = helix_ntuple->addIndexedItem("truth_phi0",      truth_ntrack ,  truth_phi0);
                status = helix_ntuple->addIndexedItem("truth_kappa",         truth_ntrack ,  truth_kappa);
                status = helix_ntuple->addIndexedItem("truth_tanl",      truth_ntrack ,  truth_tanl);
                status = helix_ntuple->addIndexedItem("truth_charge",        truth_ntrack ,  truth_charge);
                status = helix_ntuple->addIndexedItem("truth_CosTheta",     truth_ntrack ,   truth_CosTheta);
                status = helix_ntuple->addItem("nsegment",cgem_nsegment,0,1000);
                status = helix_ntuple->addIndexedItem("segmentid",  cgem_nsegment ,   cgem_segmentid);
                status = helix_ntuple->addIndexedItem("dr_ini",      cgem_nsegment ,      cgem_dr_ini);
                status = helix_ntuple->addIndexedItem("dz_ini",      cgem_nsegment ,      cgem_dz_ini);
                status = helix_ntuple->addIndexedItem("phi0_ini",        cgem_nsegment ,      cgem_phi0_ini);
                status = helix_ntuple->addIndexedItem("kappa_ini",       cgem_nsegment ,      cgem_kappa_ini);
                status = helix_ntuple->addIndexedItem("tanl_ini",        cgem_nsegment ,      cgem_tanl_ini);
                status = helix_ntuple->addIndexedItem("dr",      cgem_nsegment ,      cgem_dr);
                status = helix_ntuple->addIndexedItem("dz",      cgem_nsegment ,      cgem_dz);
                status = helix_ntuple->addIndexedItem("phi0",        cgem_nsegment ,      cgem_phi0);
                status = helix_ntuple->addIndexedItem("kappa",       cgem_nsegment ,      cgem_kappa);
                status = helix_ntuple->addIndexedItem("tanl",        cgem_nsegment ,      cgem_tanl);
                status = helix_ntuple->addIndexedItem("chisq",       cgem_nsegment ,      cgem_chisq);
                status = helix_ntuple->addIndexedItem("fitFlag",    cgem_nsegment ,   cgem_fitFlag);
                status = helix_ntuple->addIndexedItem("charge",     cgem_nsegment ,   cgem_charge);
 
                status = helix_ntuple->addItem("n_layer",n_layer,0,3);
                status = helix_ntuple->addIndexedItem("n_stripX",  n_layer ,      n_stripX);
                status = helix_ntuple->addIndexedItem("n_stripV",  n_layer ,      n_stripV);
                status = helix_ntuple->addIndexedItem("mc_phi",    n_layer ,      mc_phi);
                status = helix_ntuple->addIndexedItem("mc_v",       n_layer ,     mc_v);
                status = helix_ntuple->addIndexedItem("rec_phi",   n_layer ,      rec_phi);
                status = helix_ntuple->addIndexedItem("rec_v",      n_layer ,     rec_v);
            }
        }
    }
    check_segNum = 0;
    check_segOver = 0;
    wide_delta = 570./2/10000; //cm
    pitch = (650.-570.)/2/10000; //cm
 
    ISvcLocator* svcLocator = Gaudi::svcLocator();
    ICgemGeomSvc* ISvc;
    StatusCode sc=svcLocator->service("CgemGeomSvc", ISvc);
    CgemGeomSvc* myCgemGeomSvc=dynamic_cast<CgemGeomSvc *>(ISvc);
    int nCgemLayer = myCgemGeomSvc->getNumberOfCgemLayer();
    for(int i=0; i<nCgemLayer; i++)
    {
        CgemGeoLayer* CgemLayer = myCgemGeomSvc->getCgemLayer(i);
        CgemSegmentFun::rl[i]=(CgemLayer->getMiddleROfGapD())/10.0;//cm
        CgemSegmentFun::a_stero[i]=(CgemLayer->getAngleOfStereo())/180.0*CLHEP::pi;// degree -> radians
        CgemGeoReadoutPlane* readoutPlane = myCgemGeomSvc->getReadoutPlane(i,0);
        CgemSegmentFun::r_X[i]=readoutPlane->getRX()/10;//cm
        CgemSegmentFun::r_V[i]=readoutPlane->getRV()/10;//cm
        if(m_debug) {
            cout<<"iLayer, r, ang_stereo, Rx, Rv = "<<i<<", "<<CgemSegmentFun::rl[i]<<", "<<CgemSegmentFun::a_stero[i]<<", "<<CgemSegmentFun::r_X[i]<<", "<<CgemSegmentFun::r_V[i]<<endl;
        }
    }
    
    sc = service ("CgemCalibFunSvc", myCgemCalibSvc);
    if ( sc.isFailure() ){
        cout<< name() << "Could not load MdcCalibFunSvc!" << endl;
        return sc;
    }
 
    return StatusCode::SUCCESS;
 
}

initialize() [2/2]

StatusCode CgemSegmentFitAlg::initialize

(

)

---

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

• InstallArea/include/CgemSegmentFitAlg/CgemSegmentFitAlg/CgemSegmentFitAlg.h
• Reconstruction/CgemSegmentFitAlg/CgemSegmentFitAlg-00-00-04/CgemSegmentFitAlg/CgemSegmentFitAlg.h
• CgemSegmentFitAlg.cxx