30 #include <TApplication.h>
31 #include <TEveBrowser.h>
32 #include <TEveManager.h>
33 #include <TEveEventManager.h>
34 #include <TEveGeoNode.h>
35 #include <TEveGeoShape.h>
36 #include <TEveStraightLineSet.h>
37 #include <TEveTriangleSet.h>
38 #include <TDecompSVD.h>
41 #include <TGNumberEntry.h>
43 #include <TGeoManager.h>
44 #include <TGeoMatrix.h>
46 #include <TGeoSphere.h>
50 #include <TMatrixTSym.h>
51 #include <TMatrixDSymEigen.h>
57 #include "boost/scoped_ptr.hpp"
74 drawTrackMarkers_(
true),
82 drawCardinalRep_(
true),
90 squareRootFormalism_(
false),
100 if((!gApplication) || (gApplication && gApplication->TestBit(TApplication::kDefaultApplication))) {
101 std::cout <<
"In EventDisplay ctor: gApplication not found, creating..." << std::flush;
102 new TApplication(
"ROOT_application", 0, 0);
103 std::cout <<
"done!" << std::endl;
106 std::cout <<
"In EventDisplay ctor: gEve not found, creating..." << std::flush;
107 TEveManager::Create();
108 std::cout <<
"done!" << std::endl;
118 for(
size_t i = 0; i < opts.length(); ++i) {
153 for(
unsigned int i = 0; i <
events_.size(); i++) {
155 for(
unsigned int j = 0; j <
events_[i]->size(); j++) {
169 std::vector<Track*>* vec =
new std::vector<Track*>;
171 for(
unsigned int i = 0; i < tracks.size(); i++) {
172 vec->push_back(
new Track(*(tracks[i])));
181 std::vector<Track*>* vec =
new std::vector<Track*>;
183 for(
unsigned int i = 0; i < tracks.size(); i++) {
184 vec->push_back(
new Track(*(tracks[i])));
193 std::vector<Track*>* vec =
new std::vector<Track*>;
194 vec->push_back(
new Track(*tr));
207 if(
events_.size() == 0)
return;
226 bool resetCam =
true;
233 std::cout <<
"At event " <<
id << std::endl;
234 if (gEve->GetCurrentEvent()) {
235 gEve->GetCurrentEvent()->DestroyElements();
240 if (gEve->GetCurrentEvent()) {
241 gEve->GetCurrentEvent()->DestroyElements();
251 std::cout <<
"EventDisplay::open(); " <<
getNEvents() <<
" events loaded" << std::endl;
257 std::cout <<
"autoscaling changed the error, draw again." << std::endl;
266 gApplication->Run(kTRUE);
269 std::cout <<
"opened" << std::endl;
276 std::cout <<
"EventDisplay::drawEvent(" <<
id <<
")" << std::endl;
281 TGeoNode* top_node = gGeoManager->GetTopNode();
282 assert(top_node != NULL);
285 TObjArray* volumes = gGeoManager->GetListOfVolumes();
286 for(
int i = 0; i < volumes->GetEntriesFast(); i++) {
287 TGeoVolume* volume = dynamic_cast<TGeoVolume*>(volumes->At(i));
288 assert(volume != NULL);
289 volume->SetLineColor(12);
290 volume->SetTransparency(50);
293 TEveGeoTopNode* eve_top_node =
new TEveGeoTopNode(gGeoManager, top_node);
294 eve_top_node->IncDenyDestroy();
295 gEve->AddElement(eve_top_node);
299 for(
unsigned int i = 0; i <
events_.at(
id)->size(); i++) {
304 Track* track =
events_[id]->at(i);
305 if (! track->checkConsistency()){
306 std::cerr<<
"track is not consistent"<<std::endl;
311 boost::scoped_ptr<Track> refittedTrack(NULL);
314 std::cout <<
"Refit track:" << std::endl;
316 boost::scoped_ptr<AbsKalmanFitter> fitter;
320 fitter->setMultipleMeasurementHandling(
mmHandling_);
326 fitter->setMultipleMeasurementHandling(
mmHandling_);
327 static_cast<KalmanFitterRefTrack*>(fitter.get())->setDeltaChi2Ref(
dChi2Ref_);
331 fitter.reset(
new DAF(
false));
332 ( static_cast<KalmanFitter*>( (static_cast<DAF*>(fitter.get()))->getKalman() ) )->useSquareRootFormalism(
squareRootFormalism_);
335 fitter.reset(
new DAF());
336 ( static_cast<KalmanFitterRefTrack*>( (static_cast<DAF*>(fitter.get()))->getKalman() ) )->setDeltaChi2Ref(
dChi2Ref_);
340 fitter->setDebugLvl(std::max(0, (
int)
debugLvl_-1));
347 refittedTrack.reset(
new Track(*track));
348 refittedTrack->deleteFitterInfo();
351 refittedTrack->Print(
"C");
353 timeval startcputime, endcputime;
356 gettimeofday(&startcputime, NULL);
357 fitter->processTrack(refittedTrack.get(),
resort_);
358 gettimeofday(&endcputime, NULL);
361 std::cerr << e.
what();
362 std::cerr <<
"Exception, could not refit track" << std::endl;
366 int microseconds = 1000000*(endcputime.tv_sec - startcputime.tv_sec) + (endcputime.tv_usec - startcputime.tv_usec);
367 std::cout <<
"it took " << double(microseconds) / 1000 <<
" ms of CPU to fit the track\n";
369 if (! refittedTrack->checkConsistency()){
370 std::cerr<<
"refittedTrack is not consistent"<<std::endl;
374 track = refittedTrack.get();
383 rep = track->getCardinalRep();
384 std::cout <<
"Draw cardinal rep" << std::endl;
387 if (
repId_ >= track->getNumReps())
388 repId_ = track->getNumReps() - 1;
389 rep = track->getTrackRep(
repId_);
390 std::cout <<
"Draw rep" <<
repId_ << std::endl;
394 std::cout <<
"track " << i << std::endl;
397 track->getFitStatus(rep)->Print();
399 if (track->getFitStatus(rep)->isFitted()) {
401 std::cout <<
"fitted state: \n";
402 track->getFittedState().Print();
404 catch (Exception& e) {
405 std::cerr << e.what();
414 unsigned int numhits = track->getNumPointsWithMeasurement();
416 KalmanFitterInfo* fi;
417 KalmanFitterInfo* prevFi = 0;
418 const MeasuredStateOnPlane* fittedState(NULL);
419 const MeasuredStateOnPlane* prevFittedState(NULL);
421 for(
unsigned int j = 0; j < numhits; j++) {
425 TrackPoint* tp = track->getPointWithMeasurement(j);
426 if (! tp->hasRawMeasurements()) {
427 std::cerr<<
"trackPoint has no raw measurements"<<std::endl;
431 const AbsMeasurement* m = tp->getRawMeasurement();
432 int hit_coords_dim = m->getDim();
435 if (tp->getNumRawMeasurements() > 1) {
436 bool sameTypes(
true);
437 for (
unsigned int iM=1; iM<tp->getNumRawMeasurements(); ++iM) {
438 if (
typeid(*(tp->getRawMeasurement(iM))) !=
typeid(*m))
442 std::cerr<<
"cannot draw trackpoint containing multiple Measurements of differend types"<<std::endl;
450 if (! tp->hasFitterInfo(rep)) {
451 std::cerr<<
"trackPoint has no fitterInfo for rep"<<std::endl;
455 AbsFitterInfo* fitterInfo = tp->getFitterInfo(rep);
457 fi = dynamic_cast<KalmanFitterInfo*>(fitterInfo);
459 std::cerr<<
"can only display KalmanFitterInfo"<<std::endl;
462 if (! fi->hasPredictionsAndUpdates()) {
463 std::cerr<<
"KalmanFitterInfo does not have all predictions and updates"<<std::endl;
468 fittedState = &(fi->getFittedState(
true));
470 catch (Exception& e) {
471 std::cerr << e.what();
472 std::cerr<<
"can not get fitted state"<<std::endl;
475 prevFittedState = fittedState;
480 if (fittedState == NULL) {
481 if (fi->hasForwardUpdate()) {
482 fittedState = fi->getForwardUpdate();
484 else if (fi->hasBackwardUpdate()) {
485 fittedState = fi->getBackwardUpdate();
487 else if (fi->hasForwardPrediction()) {
488 fittedState = fi->getForwardPrediction();
490 else if (fi->hasBackwardPrediction()) {
491 fittedState = fi->getBackwardPrediction();
495 if (fittedState == NULL) {
496 std::cout <<
"cannot get any state from fitterInfo, continue.\n";
498 prevFittedState = fittedState;
502 TVector3 track_pos = fittedState->getPos();
503 double charge = fittedState->getCharge();
509 bool full_hit = (dynamic_cast<const FullMeasurement*>(m) != NULL);
510 bool planar_hit = (dynamic_cast<const PlanarMeasurement*>(m) != NULL);
511 bool planar_pixel_hit = planar_hit && hit_coords_dim == 2;
512 bool space_hit = (dynamic_cast<const SpacepointMeasurement*>(m) != NULL);
513 bool wire_hit = m && m->isLeftRightMeasurement();
514 bool wirepoint_hit = wire_hit && (dynamic_cast<const WirePointMeasurement*>(m) != NULL);
515 if (!full_hit && !planar_hit && !planar_pixel_hit && !space_hit && !wire_hit && !wirepoint_hit) {
516 std::cout <<
"Track " << i <<
", Hit " << j <<
": Unknown measurement type: skipping hit!" << std::endl;
522 unsigned int nMeas = fi->getNumMeasurements();
523 for (
unsigned int iMeas = 0; iMeas < nMeas; ++iMeas) {
525 if (iMeas > 0 && wire_hit)
528 const MeasurementOnPlane* mop = fi->getMeasurementOnPlane(iMeas);
529 const TVectorT<double>& hit_coords = mop->getState();
530 const TMatrixTSym<double>& hit_cov = mop->getCov();
535 TVector3 o = fittedState->getPlane()->getO();
536 TVector3 u = fittedState->getPlane()->getU();
537 TVector3 v = fittedState->getPlane()->getV();
541 double_t plane_size = 4;
542 TVector2 stripDir(1,0);
545 if(!planar_pixel_hit) {
546 if (dynamic_cast<RKTrackRep*>(rep) != NULL) {
547 const TMatrixD& H = mop->getHMatrix()->getMatrix();
548 stripDir.Set(H(0,3), H(0,4));
550 hit_u = hit_coords(0);
552 hit_u = hit_coords(0);
553 hit_v = hit_coords(1);
555 }
else if (wire_hit) {
556 hit_u = fabs(hit_coords(0));
557 hit_v = v*(track_pos-o);
559 hit_v = hit_coords(1);
563 if(plane_size < 4) plane_size = 4;
569 TVector3 move(0,0,0);
570 if (planar_hit) move = track_pos-o;
571 if (wire_hit) move = v*(v*(track_pos-o));
572 TEveBox* box =
boxCreator(o + move, u, v, plane_size, plane_size, 0.01);
574 box->SetMainColor(kCyan);
576 box->SetMainColor(kGray);
578 box->SetMainTransparency(50);
579 gEve->AddElement(box);
587 MeasuredStateOnPlane update ( *fi->getBackwardUpdate() );
588 update.extrapolateBy(-3.);
592 if (j > 0 && prevFi != NULL) {
596 makeLines(prevFittedState, fittedState, rep, charge > 0 ? kRed : kBlue, 1,
false,
drawErrors_, 0, 0);
601 if (j == numhits-1) {
602 MeasuredStateOnPlane update ( *fi->getForwardUpdate() );
603 update.extrapolateBy(3.);
611 if(
drawRefTrack_ && fi->hasReferenceState() && prevFi->hasReferenceState())
612 makeLines(prevFi->getReferenceState(), fi->getReferenceState(), rep, charge > 0 ? kRed + 2 : kBlue + 2, 2,
drawTrackMarkers_,
false, 3);
614 else if (j > 0 && prevFi == NULL) {
615 std::cout <<
"previous FitterInfo == NULL \n";
618 catch (Exception& e) {
619 std::cerr <<
"extrapolation failed, cannot draw track" << std::endl;
620 std::cerr << e.what();
627 TEveGeoShape* det_shape =
new TEveGeoShape(
"det_shape");
628 det_shape->IncDenyDestroy();
629 det_shape->SetShape(
new TGeoTube(std::max(0., (
double)(hit_u-0.0105/2.)), hit_u+0.0105/2., plane_size));
631 TVector3 norm = u.Cross(v);
632 TGeoRotation det_rot(
"det_rot", (u.Theta()*180)/TMath::Pi(), (u.Phi()*180)/TMath::Pi(),
633 (norm.Theta()*180)/TMath::Pi(), (norm.Phi()*180)/TMath::Pi(),
634 (v.Theta()*180)/TMath::Pi(), (v.Phi()*180)/TMath::Pi());
635 TVector3 move = v*(v*(track_pos-o));
636 TGeoCombiTrans det_trans(o(0) + move.X(),
640 det_shape->SetTransMatrix(det_trans);
641 det_shape->SetMainColor(kCyan);
642 det_shape->SetMainTransparency(25);
644 gEve->AddElement(det_shape);
656 StateOnPlane dummy(rep);
657 StateOnPlane dummy2(TVectorD(rep->getDim()), static_cast<const FullMeasurement*>(m)->constructPlane(dummy), rep);
658 MeasuredStateOnPlane sop = *(static_cast<const FullMeasurement*>(m)->constructMeasurementsOnPlane(dummy2)[0]);
661 MeasuredStateOnPlane prevSop(sop);
662 prevSop.extrapolateBy(-3);
663 makeLines(&sop, &prevSop, rep, kYellow, 1,
false,
true, 0, 0);
666 prevSop.extrapolateBy(3);
667 makeLines(&sop, &prevSop, rep, kYellow, 1,
false,
true, 0, 0);
671 if(!planar_pixel_hit) {
673 TVector3 stripDir3 = stripDir.X()*u + stripDir.Y()*v;
674 TVector3 stripDir3perp = stripDir.Y()*u - stripDir.X()*v;
675 TVector3 move = stripDir3perp*(stripDir3perp*(track_pos-o));
676 hit_box =
boxCreator((o + move + hit_u*stripDir3), stripDir3, stripDir3perp,
errorScale_*std::sqrt(hit_cov(0,0)), plane_size, 0.0105);
677 hit_box->SetMainColor(kYellow);
678 hit_box->SetMainTransparency(0);
679 gEve->AddElement(hit_box);
682 TMatrixDSymEigen eigen_values(hit_cov);
683 TEveGeoShape* cov_shape =
new TEveGeoShape(
"cov_shape");
684 cov_shape->IncDenyDestroy();
685 TVectorT<double> ev = eigen_values.GetEigenValues();
686 TMatrixT<double> eVec = eigen_values.GetEigenVectors();
687 double pseudo_res_0 =
errorScale_*std::sqrt(ev(0));
688 double pseudo_res_1 =
errorScale_*std::sqrt(ev(1));
693 double min_cov = std::min(pseudo_res_0,pseudo_res_1);
695 std::cout <<
"Track " << i <<
", Hit " << j <<
": Invalid covariance matrix (Eigenvalue < 1e-5), autoscaling not possible!" << std::endl;
697 if(min_cov < 0.049) {
698 double cor = 0.05 / min_cov;
699 std::cout <<
"Track " << i <<
", Hit " << j <<
": Pixel covariance too small, rescaling by " << cor;
710 cov_shape->SetShape(
new TGeoEltu(pseudo_res_0, pseudo_res_1, 0.0105));
711 TVector3 pix_pos = o + hit_u*u + hit_v*v;
712 TVector3 u_semiaxis = (pix_pos + eVec(0,0)*u + eVec(1,0)*v)-pix_pos;
713 TVector3 v_semiaxis = (pix_pos + eVec(0,1)*u + eVec(1,1)*v)-pix_pos;
714 TVector3 norm = u.Cross(v);
718 TGeoRotation det_rot(
"det_rot", (u_semiaxis.Theta()*180)/TMath::Pi(), (u_semiaxis.Phi()*180)/TMath::Pi(),
719 (v_semiaxis.Theta()*180)/TMath::Pi(), (v_semiaxis.Phi()*180)/TMath::Pi(),
720 (norm.Theta()*180)/TMath::Pi(), (norm.Phi()*180)/TMath::Pi());
721 TGeoCombiTrans det_trans(pix_pos(0),pix_pos(1),pix_pos(2), &det_rot);
722 cov_shape->SetTransMatrix(det_trans);
725 cov_shape->SetMainColor(kYellow);
726 cov_shape->SetMainTransparency(0);
727 gEve->AddElement(cov_shape);
736 TMatrixDSymEigen eigen_values(m->getRawHitCov());
737 TEveGeoShape* cov_shape =
new TEveGeoShape(
"cov_shape");
738 cov_shape->IncDenyDestroy();
739 cov_shape->SetShape(
new TGeoSphere(0.,1.));
740 TVectorT<double> ev = eigen_values.GetEigenValues();
741 TMatrixT<double> eVec = eigen_values.GetEigenVectors();
742 TVector3 eVec1(eVec(0,0),eVec(1,0),eVec(2,0));
743 TVector3 eVec2(eVec(0,1),eVec(1,1),eVec(2,1));
744 TVector3 eVec3(eVec(0,2),eVec(1,2),eVec(2,2));
745 const TVector3 norm = u.Cross(v);
748 static const double radDeg(180./TMath::Pi());
749 TGeoRotation det_rot(
"det_rot", eVec1.Theta()*radDeg, eVec1.Phi()*radDeg,
750 eVec2.Theta()*radDeg, eVec2.Phi()*radDeg,
751 eVec3.Theta()*radDeg, eVec3.Phi()*radDeg);
753 if (! det_rot.IsValid()){
755 if (fabs(eVec2*eVec3) > 1.e-10)
756 eVec3 = eVec1.Cross(eVec2);
758 det_rot.SetAngles(eVec1.Theta()*radDeg, eVec1.Phi()*radDeg,
759 eVec2.Theta()*radDeg, eVec2.Phi()*radDeg,
760 eVec3.Theta()*radDeg, eVec3.Phi()*radDeg);
764 double pseudo_res_0 =
errorScale_*std::sqrt(ev(0));
765 double pseudo_res_1 =
errorScale_*std::sqrt(ev(1));
766 double pseudo_res_2 =
errorScale_*std::sqrt(ev(2));
776 double min_cov = std::min(pseudo_res_0,std::min(pseudo_res_1,pseudo_res_2));
778 std::cout <<
"Track " << i <<
", Hit " << j <<
": Invalid covariance matrix (Eigenvalue < 1e-5), autoscaling not possible!" << std::endl;
780 if(min_cov <= 0.149) {
781 double cor = 0.15 / min_cov;
782 std::cout <<
"Track " << i <<
", Hit " << j <<
": Space hit covariance too small, rescaling by " << cor;
795 TGeoGenTrans det_trans(o(0),o(1),o(2),
798 pseudo_res_0, pseudo_res_1, pseudo_res_2,
800 cov_shape->SetTransMatrix(det_trans);
803 cov_shape->SetMainColor(kYellow);
804 cov_shape->SetMainTransparency(10);
805 gEve->AddElement(cov_shape);
811 TMatrixDSymEigen eigen_values(hit_cov);
812 TEveGeoShape* cov_shape =
new TEveGeoShape(
"cov_shape");
813 cov_shape->IncDenyDestroy();
814 TVectorT<double> ev = eigen_values.GetEigenValues();
815 TMatrixT<double> eVec = eigen_values.GetEigenVectors();
816 double pseudo_res_0 =
errorScale_*std::sqrt(ev(0));
817 double pseudo_res_1 =
errorScale_*std::sqrt(ev(1));
822 double min_cov = std::min(pseudo_res_0,pseudo_res_1);
824 std::cout <<
"Track " << i <<
", Hit " << j <<
": Invalid covariance matrix (Eigenvalue < 1e-5), autoscaling not possible!" << std::endl;
826 if(min_cov < 0.049) {
827 double cor = 0.05 / min_cov;
828 std::cout <<
"Track " << i <<
", Hit " << j <<
": Pixel covariance too small, rescaling by " << cor;
839 cov_shape->SetShape(
new TGeoEltu(pseudo_res_0, pseudo_res_1, 0.0105));
840 TVector3 pix_pos = o + hit_u*u + hit_v*v;
841 TVector3 u_semiaxis = (pix_pos + eVec(0,0)*u + eVec(1,0)*v)-pix_pos;
842 TVector3 v_semiaxis = (pix_pos + eVec(0,1)*u + eVec(1,1)*v)-pix_pos;
843 TVector3 norm = u.Cross(v);
847 static const double radDeg(180./TMath::Pi());
848 TGeoRotation det_rot(
"det_rot", u_semiaxis.Theta()*radDeg, u_semiaxis.Phi()*radDeg,
849 v_semiaxis.Theta()*radDeg, v_semiaxis.Phi()*radDeg,
850 norm.Theta()*radDeg, norm.Phi()*radDeg);
856 TGeoCombiTrans det_trans(pix_pos(0),pix_pos(1),pix_pos(2), &det_rot);
857 cov_shape->SetTransMatrix(det_trans);
860 cov_shape->SetMainColor(kYellow);
861 cov_shape->SetMainTransparency(0);
862 gEve->AddElement(cov_shape);
869 TEveGeoShape* cov_shape =
new TEveGeoShape(
"cov_shape");
870 cov_shape->IncDenyDestroy();
871 double pseudo_res_0 =
errorScale_*std::sqrt(hit_cov(0,0));
872 double pseudo_res_1 = plane_size;
873 if (wirepoint_hit) pseudo_res_1 =
errorScale_*std::sqrt(hit_cov(1,1));
877 if(pseudo_res_0 < 1e-5) {
878 std::cout <<
"Track " << i <<
", Hit " << j <<
": Invalid wire resolution (< 1e-5), autoscaling not possible!" << std::endl;
880 if(pseudo_res_0 < 0.0049) {
881 double cor = 0.005 / pseudo_res_0;
882 std::cout <<
"Track " << i <<
", Hit " << j <<
": Wire covariance too small, rescaling by " << cor;
889 if(wirepoint_hit && pseudo_res_1 < 1e-5) {
890 std::cout <<
"Track " << i <<
", Hit " << j <<
": Invalid wire resolution (< 1e-5), autoscaling not possible!" << std::endl;
892 if(pseudo_res_1 < 0.0049) {
893 double cor = 0.005 / pseudo_res_1;
894 std::cout <<
"Track " << i <<
", Hit " << j <<
": Wire covariance too small, rescaling by " << cor;
904 TVector3 move = v*(v*(track_pos-o));
905 hit_box =
boxCreator((o + move + hit_u*u), u, v,
errorScale_*std::sqrt(hit_cov(0,0)), pseudo_res_1, 0.0105);
906 hit_box->SetMainColor(kYellow);
907 hit_box->SetMainTransparency(0);
908 gEve->AddElement(hit_box);
910 hit_box =
boxCreator((o + move - hit_u*u), u, v,
errorScale_*std::sqrt(hit_cov(0,0)), pseudo_res_1, 0.0105);
911 hit_box->SetMainColor(kYellow);
912 hit_box->SetMainTransparency(0);
913 gEve->AddElement(hit_box);
923 prevFittedState = fittedState;
929 gEve->Redraw3D(resetCam);
938 TEveBox* box =
new TEveBox(
"detPlane_shape");
941 TVector3 norm = u.Cross(v);
946 vertices[0] = o(0) - u(0) - v(0) - norm(0);
947 vertices[1] = o(1) - u(1) - v(1) - norm(1);
948 vertices[2] = o(2) - u(2) - v(2) - norm(2);
950 vertices[3] = o(0) + u(0) - v(0) - norm(0);
951 vertices[4] = o(1) + u(1) - v(1) - norm(1);
952 vertices[5] = o(2) + u(2) - v(2) - norm(2);
954 vertices[6] = o(0) + u(0) - v(0) + norm(0);
955 vertices[7] = o(1) + u(1) - v(1) + norm(1);
956 vertices[8] = o(2) + u(2) - v(2) + norm(2);
958 vertices[9] = o(0) - u(0) - v(0) + norm(0);
959 vertices[10] = o(1) - u(1) - v(1) + norm(1);
960 vertices[11] = o(2) - u(2) - v(2) + norm(2);
962 vertices[12] = o(0) - u(0) + v(0) - norm(0);
963 vertices[13] = o(1) - u(1) + v(1) - norm(1);
964 vertices[14] = o(2) - u(2) + v(2) - norm(2);
966 vertices[15] = o(0) + u(0) + v(0) - norm(0);
967 vertices[16] = o(1) + u(1) + v(1) - norm(1);
968 vertices[17] = o(2) + u(2) + v(2) - norm(2);
970 vertices[18] = o(0) + u(0) + v(0) + norm(0);
971 vertices[19] = o(1) + u(1) + v(1) + norm(1);
972 vertices[20] = o(2) + u(2) + v(2) + norm(2);
974 vertices[21] = o(0) - u(0) + v(0) + norm(0);
975 vertices[22] = o(1) - u(1) + v(1) + norm(1);
976 vertices[23] = o(2) - u(2) + v(2) + norm(2);
979 for(
int k = 0; k < 24; k += 3) box->SetVertex((k/3), vertices[k], vertices[k+1], vertices[k+2]);
987 const Color_t& color,
const Style_t& style,
bool drawMarkers,
bool drawErrors,
double lineWidth,
int markerPos)
989 if (prevState == NULL || state == NULL) {
990 std::cerr <<
"prevState == NULL || state == NULL\n";
994 TVector3 pos, dir, oldPos, oldDir;
995 rep->getPosDir(*state, pos, dir);
996 rep->getPosDir(*prevState, oldPos, oldDir);
998 double distA = (pos-oldPos).Mag();
999 double distB = distA;
1000 if ((pos-oldPos)*oldDir < 0)
1002 if ((pos-oldPos)*dir < 0)
1004 TVector3 intermediate1 = oldPos + 0.3 * distA * oldDir;
1005 TVector3 intermediate2 = pos - 0.3 * distB * dir;
1006 TEveStraightLineSet* lineSet =
new TEveStraightLineSet;
1007 lineSet->AddLine(oldPos(0), oldPos(1), oldPos(2), intermediate1(0), intermediate1(1), intermediate1(2));
1008 lineSet->AddLine(intermediate1(0), intermediate1(1), intermediate1(2), intermediate2(0), intermediate2(1), intermediate2(2));
1009 lineSet->AddLine(intermediate2(0), intermediate2(1), intermediate2(2), pos(0), pos(1), pos(2));
1010 lineSet->SetLineColor(color);
1011 lineSet->SetLineStyle(style);
1012 lineSet->SetLineWidth(lineWidth);
1015 lineSet->AddMarker(oldPos(0), oldPos(1), oldPos(2));
1017 lineSet->AddMarker(pos(0), pos(1), pos(2));
1021 gEve->AddElement(lineSet);
1025 const MeasuredStateOnPlane* measuredState;
1027 measuredState = dynamic_cast<const MeasuredStateOnPlane*>(prevState);
1029 measuredState = dynamic_cast<const MeasuredStateOnPlane*>(state);
1031 if (measuredState != NULL) {
1035 if (markerPos == 0) {
1036 if (fabs(distA) < 1.) {
1037 distA < 0 ? distA = -1 : distA = 1;
1039 eval = 0.2 * distA * oldDir;
1042 if (fabs(distB) < 1.) {
1043 distB < 0 ? distB = -1 : distB = 1;
1045 eval = -0.2 * distB * dir;
1051 TVector3 position, direction;
1052 rep->getPosMomCov(*measuredState, position, direction, cov);
1055 TMatrixDSymEigen eigen_values(cov.GetSub(0,2, 0,2));
1056 TVectorT<double> ev = eigen_values.GetEigenValues();
1057 TMatrixT<double> eVec = eigen_values.GetEigenVectors();
1058 TVector3 eVec1, eVec2;
1060 static const double maxErr = 1000.;
1061 double ev0 = std::min(ev(0), maxErr);
1062 double ev1 = std::min(ev(1), maxErr);
1063 double ev2 = std::min(ev(2), maxErr);
1066 if (ev0 < ev1 && ev0 < ev2) {
1067 eVec1.SetXYZ(eVec(0,1),eVec(1,1),eVec(2,1));
1069 eVec2.SetXYZ(eVec(0,2),eVec(1,2),eVec(2,2));
1072 else if (ev1 < ev0 && ev1 < ev2) {
1073 eVec1.SetXYZ(eVec(0,0),eVec(1,0),eVec(2,0));
1075 eVec2.SetXYZ(eVec(0,2),eVec(1,2),eVec(2,2));
1079 eVec1.SetXYZ(eVec(0,0),eVec(1,0),eVec(2,0));
1081 eVec2.SetXYZ(eVec(0,1),eVec(1,1),eVec(2,1));
1085 if (eVec1.Cross(eVec2)*eval < 0)
1089 const TVector3 oldEVec1(eVec1);
1090 const TVector3 oldEVec2(eVec2);
1092 const int nEdges = 24;
1093 std::vector<TVector3> vertices;
1095 vertices.push_back(position);
1098 for (
int i=0; i<nEdges; ++i) {
1099 const double angle = 2*TMath::Pi()/nEdges * i;
1100 vertices.push_back(position + cos(angle)*eVec1 + sin(angle)*eVec2);
1105 DetPlane* newPlane =
new DetPlane(*(measuredState->getPlane()));
1106 newPlane->setO(position + eval);
1108 MeasuredStateOnPlane stateCopy(*measuredState);
1112 catch(Exception& e){
1113 std::cerr<<e.what();
1118 rep->getPosMomCov(stateCopy, position, direction, cov);
1121 TMatrixDSymEigen eigen_values2(cov.GetSub(0,2, 0,2));
1122 ev = eigen_values2.GetEigenValues();
1123 eVec = eigen_values2.GetEigenVectors();
1125 ev0 = std::min(ev(0), maxErr);
1126 ev1 = std::min(ev(1), maxErr);
1127 ev2 = std::min(ev(2), maxErr);
1130 if (ev0 < ev1 && ev0 < ev2) {
1131 eVec1.SetXYZ(eVec(0,1),eVec(1,1),eVec(2,1));
1133 eVec2.SetXYZ(eVec(0,2),eVec(1,2),eVec(2,2));
1136 else if (ev1 < ev0 && ev1 < ev2) {
1137 eVec1.SetXYZ(eVec(0,0),eVec(1,0),eVec(2,0));
1139 eVec2.SetXYZ(eVec(0,2),eVec(1,2),eVec(2,2));
1143 eVec1.SetXYZ(eVec(0,0),eVec(1,0),eVec(2,0));
1145 eVec2.SetXYZ(eVec(0,1),eVec(1,1),eVec(2,1));
1149 if (eVec1.Cross(eVec2)*eval < 0)
1153 if (oldEVec1*eVec1 < 0) {
1159 double angle0 = eVec1.Angle(oldEVec1);
1160 if (eVec1*(eval.Cross(oldEVec1)) < 0)
1162 for (
int i=0; i<nEdges; ++i) {
1163 const double angle = 2*TMath::Pi()/nEdges * i - angle0;
1164 vertices.push_back(position + cos(angle)*eVec1 + sin(angle)*eVec2);
1167 vertices.push_back(position);
1170 TEveTriangleSet* error_shape =
new TEveTriangleSet(vertices.size(), nEdges*2);
1171 for(
unsigned int k = 0; k < vertices.size(); ++k) {
1172 error_shape->SetVertex(k, vertices[k].X(), vertices[k].Y(), vertices[k].Z());
1175 assert(vertices.size() == 2*nEdges+2);
1178 for (
int i=0; i<nEdges; ++i) {
1180 error_shape->SetTriangle(iTri++, i+1, i+1+nEdges, (i+1)%nEdges+1);
1181 error_shape->SetTriangle(iTri++, (i+1)%nEdges+1, i+1+nEdges, (i+1)%nEdges+1+nEdges);
1187 error_shape->SetMainColor(color);
1188 error_shape->SetMainTransparency(25);
1189 gEve->AddElement(error_shape);
1197 TEveBrowser* browser = gEve->GetBrowser();
1198 browser->StartEmbedding(TRootBrowser::kLeft);
1200 TGMainFrame* frmMain =
new TGMainFrame(gClient->GetRoot(), 1000, 600);
1201 frmMain->SetWindowName(
"XX GUI");
1202 frmMain->SetCleanup(kDeepCleanup);
1205 TGTextButton* tb = 0;
1208 TGHorizontalFrame* hf =
new TGHorizontalFrame(frmMain); {
1210 lbl =
new TGLabel(hf,
"Go to event: ");
1212 guiEvent =
new TGNumberEntry(hf, 0, 9,999, TGNumberFormat::kNESInteger,
1213 TGNumberFormat::kNEANonNegative,
1214 TGNumberFormat::kNELLimitMinMax,
1217 guiEvent->Connect(
"ValueSet(Long_t)",
"genfit::EventDisplay", fh,
"guiGoto()");
1220 tb =
new TGTextButton(hf,
"Redraw Event");
1222 tb->Connect(
"Clicked()",
"genfit::EventDisplay", fh,
"guiGoto()");
1224 frmMain->AddFrame(hf);
1227 hf =
new TGHorizontalFrame(frmMain); {
1228 lbl =
new TGLabel(hf,
"\n Draw Options");
1231 frmMain->AddFrame(hf);
1233 hf =
new TGHorizontalFrame(frmMain); {
1237 guiDrawGeometry_->Connect(
"Toggled(Bool_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1239 frmMain->AddFrame(hf);
1241 hf =
new TGHorizontalFrame(frmMain); {
1245 guiDrawDetectors_->Connect(
"Toggled(Bool_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1247 frmMain->AddFrame(hf);
1249 hf =
new TGHorizontalFrame(frmMain); {
1253 guiDrawHits_->Connect(
"Toggled(Bool_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1255 frmMain->AddFrame(hf);
1259 hf =
new TGHorizontalFrame(frmMain); {
1263 guiDrawPlanes_->Connect(
"Toggled(Bool_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1265 frmMain->AddFrame(hf);
1267 hf =
new TGHorizontalFrame(frmMain); {
1271 guiDrawTrackMarkers_->Connect(
"Toggled(Bool_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1273 frmMain->AddFrame(hf);
1276 hf =
new TGHorizontalFrame(frmMain); {
1280 guiDrawTrack_->Connect(
"Toggled(Bool_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1282 frmMain->AddFrame(hf);
1284 hf =
new TGHorizontalFrame(frmMain); {
1288 guiDrawRefTrack_->Connect(
"Toggled(Bool_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1290 frmMain->AddFrame(hf);
1292 hf =
new TGHorizontalFrame(frmMain); {
1296 guiDrawErrors_->Connect(
"Toggled(Bool_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1298 frmMain->AddFrame(hf);
1300 hf =
new TGHorizontalFrame(frmMain); {
1304 guiDrawForward_->Connect(
"Toggled(Bool_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1306 frmMain->AddFrame(hf);
1308 hf =
new TGHorizontalFrame(frmMain); {
1312 guiDrawBackward_->Connect(
"Toggled(Bool_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1314 frmMain->AddFrame(hf);
1317 hf =
new TGHorizontalFrame(frmMain); {
1321 guiDrawAutoScale_->Connect(
"Toggled(Bool_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1323 frmMain->AddFrame(hf);
1325 hf =
new TGHorizontalFrame(frmMain); {
1329 guiDrawScaleMan_->Connect(
"Toggled(Bool_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1331 frmMain->AddFrame(hf);
1333 hf =
new TGHorizontalFrame(frmMain); {
1335 TGNumberFormat::kNEANonNegative,
1336 TGNumberFormat::kNELLimitMinMax,
1339 guiErrorScale_->Connect(
"ValueSet(Long_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1340 lbl =
new TGLabel(hf,
"Error scale");
1343 frmMain->AddFrame(hf);
1347 hf =
new TGHorizontalFrame(frmMain); {
1348 lbl =
new TGLabel(hf,
"\n TrackRep options");
1351 frmMain->AddFrame(hf);
1353 hf =
new TGHorizontalFrame(frmMain); {
1357 guiDrawCardinalRep_->Connect(
"Toggled(Bool_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1359 frmMain->AddFrame(hf);
1361 hf =
new TGHorizontalFrame(frmMain); {
1362 guiRepId_ =
new TGNumberEntry(hf,
repId_, 6,999, TGNumberFormat::kNESInteger,
1363 TGNumberFormat::kNEANonNegative,
1364 TGNumberFormat::kNELLimitMinMax,
1367 guiRepId_->Connect(
"ValueSet(Long_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1368 lbl =
new TGLabel(hf,
"Else draw rep with id");
1371 frmMain->AddFrame(hf);
1373 hf =
new TGHorizontalFrame(frmMain); {
1377 guiDrawAllTracks_->Connect(
"Toggled(Bool_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1379 frmMain->AddFrame(hf);
1381 hf =
new TGHorizontalFrame(frmMain); {
1383 TGNumberFormat::kNEANonNegative,
1384 TGNumberFormat::kNELLimitMinMax,
1387 guiTrackId_->Connect(
"ValueSet(Long_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1388 lbl =
new TGLabel(hf,
"Else draw track nr. ");
1391 frmMain->AddFrame(hf);
1395 frmMain->MapSubwindows();
1397 frmMain->MapWindow();
1399 browser->StopEmbedding();
1400 browser->SetTabTitle(
"Draw Control", 0);
1403 browser->StartEmbedding(TRootBrowser::kLeft);
1404 TGMainFrame* frmMain2 =
new TGMainFrame(gClient->GetRoot(), 1000, 600);
1405 frmMain2->SetWindowName(
"XX GUI");
1406 frmMain2->SetCleanup(kDeepCleanup);
1408 hf =
new TGHorizontalFrame(frmMain2); {
1410 lbl =
new TGLabel(hf,
"Go to event: ");
1412 guiEvent2 =
new TGNumberEntry(hf, 0, 9,999, TGNumberFormat::kNESInteger,
1413 TGNumberFormat::kNEANonNegative,
1414 TGNumberFormat::kNELLimitMinMax,
1417 guiEvent2->Connect(
"ValueSet(Long_t)",
"genfit::EventDisplay", fh,
"guiGoto2()");
1420 tb =
new TGTextButton(hf,
"Redraw Event");
1422 tb->Connect(
"Clicked()",
"genfit::EventDisplay", fh,
"guiGoto()");
1424 frmMain2->AddFrame(hf);
1426 hf =
new TGHorizontalFrame(frmMain2); {
1427 lbl =
new TGLabel(hf,
"\n Fitting options");
1430 frmMain2->AddFrame(hf);
1432 hf =
new TGHorizontalFrame(frmMain2); {
1433 guiRefit_ =
new TGCheckButton(hf,
"Refit");
1436 guiRefit_->Connect(
"Toggled(Bool_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1438 frmMain2->AddFrame(hf);
1440 hf =
new TGHorizontalFrame(frmMain2); {
1442 TGNumberFormat::kNEANonNegative,
1443 TGNumberFormat::kNELLimitMinMax,
1446 guiDebugLvl_->Connect(
"ValueSet(Long_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1447 lbl =
new TGLabel(hf,
"debug level");
1450 frmMain2->AddFrame(hf);
1452 hf =
new TGHorizontalFrame(frmMain2); {
1454 guiFitterId_->Connect(
"Clicked(Int_t)",
"genfit::EventDisplay", fh,
"guiSelectFitterId(int)");
1455 hf->AddFrame(
guiFitterId_,
new TGLayoutHints(kLHintsTop));
1456 TGRadioButton* fitterId_button =
new TGRadioButton(
guiFitterId_,
"Simple Kalman");
1458 new TGRadioButton(
guiFitterId_,
"DAF w/ simple Kalman");
1459 new TGRadioButton(
guiFitterId_,
"DAF w/ reference Kalman");
1460 fitterId_button->SetDown(
true,
false);
1463 frmMain2->AddFrame(hf);
1465 hf =
new TGHorizontalFrame(frmMain2); {
1466 guiMmHandling_ =
new TGButtonGroup(hf,
"Multiple measurement handling in Kalman:");
1467 guiMmHandling_->Connect(
"Clicked(Int_t)",
"genfit::EventDisplay", fh,
"guiSelectMmHandling(int)");
1469 TGRadioButton* mmHandling_button =
new TGRadioButton(
guiMmHandling_,
"weighted average");
1471 new TGRadioButton(
guiMmHandling_,
"weighted, closest to reference");
1472 new TGRadioButton(
guiMmHandling_,
"unweighted, closest to reference");
1473 new TGRadioButton(
guiMmHandling_,
"weighted, closest to prediction");
1474 new TGRadioButton(
guiMmHandling_,
"unweighted, closest to prediction");
1475 new TGRadioButton(
guiMmHandling_,
"weighted, closest to reference for WireMeasurements, weighted average else");
1476 new TGRadioButton(
guiMmHandling_,
"unweighted, closest to reference for WireMeasurements, unweighted average else");
1477 new TGRadioButton(
guiMmHandling_,
"weighted, closest to prediction for WireMeasurements, weighted average else");
1478 new TGRadioButton(
guiMmHandling_,
"unweighted, closest to prediction for WireMeasurements, unweighted average else");
1479 mmHandling_button->SetDown(
true,
false);
1482 frmMain2->AddFrame(hf);
1484 hf =
new TGHorizontalFrame(frmMain2); {
1490 frmMain2->AddFrame(hf);
1492 hf =
new TGHorizontalFrame(frmMain2); {
1493 guiDPVal_ =
new TGNumberEntry(hf,
dPVal_, 6,9999, TGNumberFormat::kNESReal,
1494 TGNumberFormat::kNEANonNegative,
1495 TGNumberFormat::kNELLimitMinMax,
1498 guiDPVal_->Connect(
"ValueSet(Long_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1499 lbl =
new TGLabel(hf,
"delta pVal (convergence criterium)");
1502 frmMain2->AddFrame(hf);
1504 hf =
new TGHorizontalFrame(frmMain2); {
1506 TGNumberFormat::kNEANonNegative,
1507 TGNumberFormat::kNELLimitMinMax,
1510 guiRelChi2_->Connect(
"ValueSet(Long_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1511 lbl =
new TGLabel(hf,
"rel chi^2 change (non-convergence criterium)");
1514 frmMain2->AddFrame(hf);
1516 hf =
new TGHorizontalFrame(frmMain2); {
1518 TGNumberFormat::kNEANonNegative,
1519 TGNumberFormat::kNELLimitMinMax,
1522 guiDChi2Ref_->Connect(
"ValueSet(Long_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1523 lbl =
new TGLabel(hf,
"min chi^2 change for re-calculating reference track (Ref Kalman)");
1526 frmMain2->AddFrame(hf);
1528 hf =
new TGHorizontalFrame(frmMain2); {
1530 TGNumberFormat::kNEANonNegative,
1531 TGNumberFormat::kNELLimitMinMax,
1534 guiNMinIter_->Connect(
"ValueSet(Long_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1535 lbl =
new TGLabel(hf,
"Minimum nr of iterations");
1538 frmMain2->AddFrame(hf);
1540 hf =
new TGHorizontalFrame(frmMain2); {
1542 TGNumberFormat::kNEANonNegative,
1543 TGNumberFormat::kNELLimitMinMax,
1546 guiNMaxIter_->Connect(
"ValueSet(Long_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1547 lbl =
new TGLabel(hf,
"Maximum nr of iterations");
1550 frmMain2->AddFrame(hf);
1552 hf =
new TGHorizontalFrame(frmMain2); {
1554 TGNumberFormat::kNEAAnyNumber,
1555 TGNumberFormat::kNELLimitMinMax,
1558 guiNMaxFailed_->Connect(
"ValueSet(Long_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1559 lbl =
new TGLabel(hf,
"Maximum nr of failed hits");
1562 frmMain2->AddFrame(hf);
1565 hf =
new TGHorizontalFrame(frmMain2); {
1566 guiResort_ =
new TGCheckButton(hf,
"Resort track");
1569 guiResort_->Connect(
"Toggled(Bool_t)",
"genfit::EventDisplay", fh,
"guiSetDrawParams()");
1571 frmMain2->AddFrame(hf);
1576 frmMain2->MapSubwindows();
1578 frmMain2->MapWindow();
1580 browser->StopEmbedding();
1581 browser->SetTabTitle(
"Refit Control", 0);
1586 Long_t n =
guiEvent->GetNumberEntry()->GetIntNumber();
1592 Long_t n =
guiEvent2->GetNumberEntry()->GetIntNumber();