diff --git a/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.cxx b/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.cxx index 8aaacd7c28b4dce0e6870f973b3f324cbecfc97b..69b4a256e3d7eabe99bcd70dc4fc190c948790b2 100755 --- a/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.cxx +++ b/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.cxx @@ -70,11 +70,12 @@ namespace Tracker , m_residual_x_plane6(0) , m_residual_x_plane7(0) , m_residual_x_plane8(0) - , m_chi2(0) - , m_edm(0) - , m_ndf(0) + // , m_chi2(0) + // , m_edm(0) + // , m_ndf(0) , m_thistSvc("THistSvc", name) { + declareProperty("MaxChi2", m_maxchi2=20); } //----------------------------------------------------------------------- @@ -135,9 +136,9 @@ namespace Tracker m_residual_x_plane6=new TH1D("sp_residual_x_plane6","sp_residual_x_plane6",2000,-1,1); m_residual_x_plane7=new TH1D("sp_residual_x_plane7","sp_residual_x_plane7",2000,-1,1); m_residual_x_plane8=new TH1D("sp_residual_x_plane8","sp_residual_x_plane8",2000,-1,1); - m_chi2=new TH1D("chi2","chi2",100,0,100); - m_edm=new TH1D("edm","edm",100,0,0.01); - m_ndf=new TH1D("ndf","ndf",10,0,10); + // m_chi2=new TH1D("chi2","chi2",100,0,100); + // m_edm=new TH1D("edm","edm",100,0,0.01); + // m_ndf=new TH1D("ndf","ndf",10,0,10); CHECK(m_thistSvc->regTree("/TrackerSPFit/spfit",m_tree)); CHECK(m_thistSvc->regHist("/TrackerSPFit/sp/sp_x",m_hist_x)); CHECK(m_thistSvc->regHist("/TrackerSPFit/sp/sp_y",m_hist_y)); @@ -175,9 +176,9 @@ namespace Tracker CHECK(m_thistSvc->regHist("/TrackerSPFit/sp/sp_residual_x_plane6",m_residual_x_plane6)); CHECK(m_thistSvc->regHist("/TrackerSPFit/sp/sp_residual_x_plane7",m_residual_x_plane7)); CHECK(m_thistSvc->regHist("/TrackerSPFit/sp/sp_residual_x_plane8",m_residual_x_plane8)); - CHECK(m_thistSvc->regHist("/TrackerSPFit/sp/sp_chi2",m_chi2)); - CHECK(m_thistSvc->regHist("/TrackerSPFit/sp/sp_edm",m_edm)); - CHECK(m_thistSvc->regHist("/TrackerSPFit/sp/sp_ndf",m_ndf)); + // CHECK(m_thistSvc->regHist("/TrackerSPFit/sp/sp_chi2",m_chi2)); + // CHECK(m_thistSvc->regHist("/TrackerSPFit/sp/sp_edm",m_edm)); + // CHECK(m_thistSvc->regHist("/TrackerSPFit/sp/sp_ndf",m_ndf)); ATH_MSG_INFO( "TrackerSPFit::initialized" ); return StatusCode::SUCCESS; } @@ -198,9 +199,9 @@ namespace Tracker m_tree->Branch("sp_y_residual",&m_sp_y_residual); m_tree->Branch("sp_z_residual",&m_sp_z_residual); m_tree->Branch("trackId",&m_trackId); - m_tree->Branch("track_chi2",&m_track_chi2); - m_tree->Branch("track_edm",&m_track_edm); - m_tree->Branch("track_ndf",&m_track_ndf); + // m_tree->Branch("track_chi2",&m_track_chi2); + // m_tree->Branch("track_edm",&m_track_edm); + // m_tree->Branch("track_ndf",&m_track_ndf); // m_tree->Branch("track_p0",&m_track_p0); // m_tree->Branch("track_p1",&m_track_p1); // m_tree->Branch("track_p2",&m_track_p2); @@ -224,9 +225,9 @@ namespace Tracker m_sp_y_residual.clear(); m_sp_z_residual.clear(); m_trackId.clear(); - m_track_chi2.clear(); - m_track_edm.clear(); - m_track_ndf.clear(); + // m_track_chi2.clear(); + // m_track_edm.clear(); + // m_track_ndf.clear(); m_sp_station.clear(); m_sp_layer.clear(); m_sp_module.clear(); @@ -325,10 +326,9 @@ namespace Tracker } ATH_MSG_INFO( "TrackerSPFit number of spacepoints in each stations: "<<sp_sta0.size()<<" "<<sp_sta1.size()<<" "<<sp_sta2.size() ); - double maxchi2=10.; bool ifTrack=false; if(sp_sta0.size()>2) { - auto tracks_sta0=makeTrackSeg(sp_sta0,maxchi2); + auto tracks_sta0=makeTrackSeg(sp_sta0,m_maxchi2); ATH_MSG_INFO( "TrackerSPFit found "<<tracks_sta0.size()<<" track segments in station 0 " ); if(tracks_sta0.size()>0){ ifTrack=true; @@ -349,12 +349,13 @@ namespace Tracker m_sp_x_residual.push_back(pre.x()-pos.x()); m_sp_y_residual.push_back(pre.y()-pos.y()); m_sp_z_residual.push_back(pre.z()-pos.z()); + ATH_MSG_DEBUG(" fill the result "<<pre.x()<<" "<<pre.y()<<" "<<pre.z()<<" "<<pos.x()<<" "<<pos.y()<<" "<<pos.z()); m_sp_station.push_back(0); m_sp_layer.push_back(track.layer); m_sp_module.push_back(track.module); - m_track_chi2.push_back(track.chi2); - m_track_edm.push_back(track.edm); - m_track_ndf.push_back(track.ndf); + // m_track_chi2.push_back(track.chi2); + // m_track_edm.push_back(track.edm); + // m_track_ndf.push_back(track.ndf); double* par=const_cast<double*>(track.p); // m_track_p0.push_back(par[0]); // m_track_p1.push_back(par[1]); @@ -367,7 +368,7 @@ namespace Tracker } } if(sp_sta1.size()>2) { - auto tracks_sta1=makeTrackSeg(sp_sta1,maxchi2); + auto tracks_sta1=makeTrackSeg(sp_sta1,m_maxchi2); ATH_MSG_DEBUG( "TrackerSPFit found "<<tracks_sta1.size()<<" track segments in station 1 " ); if(tracks_sta1.size()>0){ ifTrack=true; @@ -391,9 +392,9 @@ namespace Tracker m_sp_station.push_back(1); m_sp_layer.push_back(track.layer); m_sp_module.push_back(track.module); - m_track_chi2.push_back(track.chi2); - m_track_edm.push_back(track.edm); - m_track_ndf.push_back(track.ndf); + // m_track_chi2.push_back(track.chi2); + // m_track_edm.push_back(track.edm); + // m_track_ndf.push_back(track.ndf); double* par=const_cast<double*>(track.p); // m_track_p0.push_back(par[0]); // m_track_p1.push_back(par[1]); @@ -406,7 +407,7 @@ namespace Tracker } } if(sp_sta2.size()>2) { - auto tracks_sta2=makeTrackSeg(sp_sta2,maxchi2); + auto tracks_sta2=makeTrackSeg(sp_sta2,m_maxchi2); ATH_MSG_DEBUG( "TrackerSPFit found "<<tracks_sta2.size()<<" track segments in station 2 " ); if(tracks_sta2.size()>0){ ifTrack=true; @@ -430,9 +431,9 @@ namespace Tracker m_sp_station.push_back(2); m_sp_layer.push_back(track.layer); m_sp_module.push_back(track.module); - m_track_chi2.push_back(track.chi2); - m_track_edm.push_back(track.edm); - m_track_ndf.push_back(track.ndf); + // m_track_chi2.push_back(track.chi2); + // m_track_edm.push_back(track.edm); + // m_track_ndf.push_back(track.ndf); double* par=const_cast<double*>(track.p); // m_track_p0.push_back(par[0]); // m_track_p1.push_back(par[1]); @@ -489,50 +490,94 @@ namespace Tracker std::vector<TrackerSPFit::SP_TSOS> TrackerSPFit::makeTrackSeg(SP_Seed sp0, SP_Seed sp1, SP_Seed sp2, double maxchi2) const{ std::vector<SP_TSOS> spt; spt.clear(); - TGraph2DErrors* gra= new TGraph2DErrors(); - gra->SetPoint(0,sp0.pos.x(),sp0.pos.y(),sp0.pos.z()); - gra->SetPoint(1,sp1.pos.x(),sp1.pos.y(),sp1.pos.z()); - gra->SetPoint(2,sp2.pos.x(),sp2.pos.y(),sp2.pos.z()); - gra->SetPointError(0,sqrt((sp0.cov)(0,0)),sqrt((sp0.cov)(1,1)),sqrt((sp0.cov)(2,2))); - gra->SetPointError(1,sqrt((sp1.cov)(0,0)),sqrt((sp1.cov)(1,1)),sqrt((sp1.cov)(2,2))); - gra->SetPointError(2,sqrt((sp2.cov)(0,0)),sqrt((sp2.cov)(1,1)),sqrt((sp2.cov)(2,2))); - - ROOT::Fit::Fitter fitter; - SumDistance2 sumdist(gra); - ROOT::Math::Functor fcn(sumdist,4); - double initParam[4]={1,1,1,1}; - fitter.SetFCN(fcn,initParam); - // for (int i = 0; i < 4; ++i) fitter.Config().ParSettings(i).SetStepSize(0.01); - if(fitter.FitFCN()){ - const ROOT::Fit::FitResult & result =fitter.Result(); - double chi2=result.MinFcnValue(); - double edm=result.Edm(); - int ndf=result.Ndf(); - ATH_MSG_DEBUG( "TrackerSPFit::makeTrackSeg(), track chi2 = "<<chi2<<" ; edm = "<<edm<<" ; ndf = "<<ndf ); - m_chi2->Fill(chi2); - m_edm->Fill(edm); - m_ndf->Fill(ndf); - - if(chi2<maxchi2){ - //result.Print(std::cout); - const double *fitParam=result.GetParams(); - Amg::Vector3D err0(sqrt((sp0.cov)(0,0)),sqrt((sp0.cov)(1,1)),sqrt((sp0.cov)(2,2))); - Amg::Vector3D err1(sqrt((sp1.cov)(0,0)),sqrt((sp1.cov)(1,1)),sqrt((sp1.cov)(2,2))); - Amg::Vector3D err2(sqrt((sp2.cov)(0,0)),sqrt((sp2.cov)(1,1)),sqrt((sp2.cov)(2,2))); - struct SP_TSOS tsos0{sp0.pos,predicted(sp0.pos.z(),fitParam),err0,sp0.layer,sp0.module,chi2,edm,ndf,fitParam}; - struct SP_TSOS tsos1{sp1.pos,predicted(sp1.pos.z(),fitParam),err1,sp1.layer,sp0.module,chi2,edm,ndf,fitParam}; - struct SP_TSOS tsos2{sp2.pos,predicted(sp2.pos.z(),fitParam),err2,sp2.layer,sp0.module,chi2,edm,ndf,fitParam}; + for(int isp=0;isp<3;isp++){ + Amg::Vector3D error; + + if(isp==0){ + error=Amg::Vector3D(sqrt((sp0.cov)(0,0)),sqrt((sp0.cov)(1,1)),sqrt((sp0.cov)(2,2))); + double resx= sp0.pos.x()-(sp1.pos.x()+sp2.pos.x())/2.; + double resy= sp0.pos.y()-(sp1.pos.y()-((sp1.pos.y()-sp2.pos.y())*(sp1.pos.z()-sp0.pos.z())/(sp1.pos.z()-sp2.pos.z()))); + Amg::Vector3D predictedpoint2(resx,resy,sp0.pos.z()); + + struct SP_TSOS tsos0{sp0.pos,predictedpoint2,error,sp0.layer,sp0.module}; + spt.push_back(tsos0); + } + else if(isp==1){ + error=Amg::Vector3D(sqrt((sp1.cov)(0,0)),sqrt((sp1.cov)(1,1)),sqrt((sp1.cov)(2,2))); + double resx= sp1.pos.x()-(sp2.pos.x()+sp0.pos.x())/2.; + double resy= sp1.pos.y()-(sp2.pos.y()-((sp2.pos.y()-sp0.pos.y())*(sp2.pos.z()-sp1.pos.z())/(sp2.pos.z()-sp0.pos.z()))); + Amg::Vector3D predictedpoint2(resx,resy,sp1.pos.z()); + + struct SP_TSOS tsos0{sp1.pos,predictedpoint2,error,sp1.layer,sp1.module}; + spt.push_back(tsos0); + } + else{ + error=Amg::Vector3D(sqrt((sp2.cov)(0,0)),sqrt((sp2.cov)(1,1)),sqrt((sp2.cov)(2,2))); + double resx= sp2.pos.x()-(sp1.pos.x()+sp0.pos.x())/2.; + double resy= sp2.pos.y()-(sp1.pos.y()-((sp1.pos.y()-sp0.pos.y())*(sp1.pos.z()-sp2.pos.z())/(sp1.pos.z()-sp0.pos.z()))); + Amg::Vector3D predictedpoint2(resx,resy,sp2.pos.z()); + + struct SP_TSOS tsos0{sp2.pos,predictedpoint2,error,sp2.layer,sp2.module}; spt.push_back(tsos0); - spt.push_back(tsos1); - spt.push_back(tsos2); } - } + } return spt; } + /* + std::vector<TrackerSPFit::SP_TSOS> TrackerSPFit::makeTrackSeg(SP_Seed sp0, SP_Seed sp1, SP_Seed sp2, double maxchi2) const{ + std::vector<SP_TSOS> spt; + spt.clear(); + TGraph2DErrors* gra= new TGraph2DErrors(); + gra->SetPoint(0,sp0.pos.x(),sp0.pos.y(),sp0.pos.z()); + gra->SetPoint(1,sp1.pos.x(),sp1.pos.y(),sp1.pos.z()); + gra->SetPoint(2,sp2.pos.x(),sp2.pos.y(),sp2.pos.z()); + gra->SetPointError(0,sqrt((sp0.cov)(0,0)),sqrt((sp0.cov)(1,1)),sqrt((sp0.cov)(2,2))); + gra->SetPointError(1,sqrt((sp1.cov)(0,0)),sqrt((sp1.cov)(1,1)),sqrt((sp1.cov)(2,2))); + gra->SetPointError(2,sqrt((sp2.cov)(0,0)),sqrt((sp2.cov)(1,1)),sqrt((sp2.cov)(2,2))); + + std::cout<<" spacepoints in the station"<<std::endl; + std::cout<<"sp 0 "<<sp0.pos.x()<<" "<<sp0.pos.y()<<" "<<sp0.pos.z()<<" "<<sqrt((sp0.cov)(0,0))<<" "<<sqrt((sp0.cov)(1,1))<<" "<<sqrt((sp0.cov)(2.2))<<std::endl; + std::cout<<"sp 1 "<<sp1.pos.x()<<" "<<sp1.pos.y()<<" "<<sp1.pos.z()<<" "<<sqrt((sp1.cov)(0,0))<<" "<<sqrt((sp1.cov)(1,1))<<" "<<sqrt((sp1.cov)(2.2))<<std::endl; + std::cout<<"sp 2 "<<sp2.pos.x()<<" "<<sp2.pos.y()<<" "<<sp2.pos.z()<<" "<<sqrt((sp2.cov)(0,0))<<" "<<sqrt((sp2.cov)(1,1))<<" "<<sqrt((sp2.cov)(2.2))<<std::endl; + ROOT::Fit::Fitter fitter; + SumDistance2 sumdist(gra); + ROOT::Math::Functor fcn(sumdist,4); + double initParam[4]={1,1,1,1}; + fitter.SetFCN(fcn,initParam); +// for (int i = 0; i < 4; ++i) fitter.Config().ParSettings(i).SetStepSize(0.01); +if(fitter.FitFCN()){ +const ROOT::Fit::FitResult & result =fitter.Result(); +double chi2=result.MinFcnValue(); +double edm=result.Edm(); +int ndf=result.Ndf(); +ATH_MSG_DEBUG( "TrackerSPFit::makeTrackSeg(), track chi2 = "<<chi2<<" ; edm = "<<edm<<" ; ndf = "<<ndf ); +m_chi2->Fill(chi2); +m_edm->Fill(edm); +m_ndf->Fill(ndf); + +if(chi2<maxchi2){ +result.Print(std::cout); +const double *fitParam=result.GetParams(); +std::cout<<" fit status: ook "<<chi2<<" "<<edm<<" "<<ndf<<" "<<fitParam[0]<<" "<<fitParam[1]<<" "<<fitParam[2]<<" "<<fitParam[3]<<std::endl; +Amg::Vector3D err0(sqrt((sp0.cov)(0,0)),sqrt((sp0.cov)(1,1)),sqrt((sp0.cov)(2,2))); +Amg::Vector3D err1(sqrt((sp1.cov)(0,0)),sqrt((sp1.cov)(1,1)),sqrt((sp1.cov)(2,2))); +Amg::Vector3D err2(sqrt((sp2.cov)(0,0)),sqrt((sp2.cov)(1,1)),sqrt((sp2.cov)(2,2))); +struct SP_TSOS tsos0{sp0.pos,predicted(sp0.pos.z(),fitParam),err0,sp0.layer,sp0.module,chi2,edm,ndf,fitParam}; +struct SP_TSOS tsos1{sp1.pos,predicted(sp1.pos.z(),fitParam),err1,sp1.layer,sp0.module,chi2,edm,ndf,fitParam}; +struct SP_TSOS tsos2{sp2.pos,predicted(sp2.pos.z(),fitParam),err2,sp2.layer,sp0.module,chi2,edm,ndf,fitParam}; +spt.push_back(tsos0); +spt.push_back(tsos1); +spt.push_back(tsos2); +} +} - Amg::Vector3D TrackerSPFit::predicted(double z, const double *p) const { - return Amg::Vector3D(p[0]+p[1]*z,p[2]+p[3]*z,z); - } +return spt; +} +*/ + +Amg::Vector3D TrackerSPFit::predicted(double z, const double *p) const { + return Amg::Vector3D(p[0]+p[1]*z,p[2]+p[3]*z,z); +} } diff --git a/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.h b/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.h index 6047dc8d39b77b7bbdcf374f9fd856107a4d400f..432befef2fc31d446aa2fd25bfbce8afc0c0d0c6 100755 --- a/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.h +++ b/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.h @@ -60,7 +60,8 @@ namespace Tracker virtual StatusCode finalize() override; struct SP_Seed{Amg::Vector3D pos; Amg::MatrixX cov; int layer; int module;}; - struct SP_TSOS{Amg::Vector3D pos; Amg::Vector3D predicted; Amg::Vector3D err; int layer;int module; double chi2; double edm; int ndf; const double *p;}; + //struct SP_TSOS{Amg::Vector3D pos; Amg::Vector3D predicted; Amg::Vector3D err; int layer;int module; double chi2; double edm; int ndf; const double *p;}; + struct SP_TSOS{Amg::Vector3D pos; Amg::Vector3D predicted; Amg::Vector3D err; int layer;int module; const double *p;}; std::vector<SP_TSOS> makeTrackSeg(std::vector<SP_Seed> sps, double maxchi2) const; std::vector<SP_TSOS> makeTrackSeg(SP_Seed sp0, SP_Seed sp1, SP_Seed sp2, double maxchi2) const; Amg::Vector3D predicted(double z, const double *p) const; @@ -161,9 +162,9 @@ TH1* m_residual_x_plane5; TH1* m_residual_x_plane6; TH1* m_residual_x_plane7; TH1* m_residual_x_plane8; -TH1* m_chi2; -TH1* m_edm; -TH1* m_ndf; +//TH1* m_chi2; +//TH1* m_edm; +//TH1* m_ndf; ServiceHandle<ITHistSvc> m_thistSvc; mutable int m_eventNumber; @@ -180,9 +181,9 @@ mutable std::vector<double> m_sp_x_residual; mutable std::vector<double> m_sp_y_residual; mutable std::vector<double> m_sp_z_residual; mutable std::vector<int> m_trackId; -mutable std::vector<double> m_track_chi2; -mutable std::vector<double> m_track_edm; -mutable std::vector<int> m_track_ndf; +//mutable std::vector<double> m_track_chi2; +//mutable std::vector<double> m_track_edm; +//mutable std::vector<int> m_track_ndf; // mutable std::vector<double> m_track_p0; // mutable std::vector<double> m_track_p1; // mutable std::vector<double> m_track_p2; @@ -190,6 +191,7 @@ mutable std::vector<int> m_track_ndf; mutable std::vector<int> m_sp_station; mutable std::vector<int> m_sp_layer; mutable std::vector<int> m_sp_module; +double m_maxchi2; };