error_ellipse.C

Go to the documentation of this file.

/// \file
/// \ingroup tutorial_eve_7
///
/// \macro_code
///
/// \author Alja Mrak-Tadel
 
#include <ROOT/REveElement.hxx>
#include <ROOT/REveScene.hxx>
#include <ROOT/REveViewer.hxx>
#include <ROOT/REveManager.hxx>
#include <ROOT/REvePointSet.hxx>
 
using namespace ROOT::Experimental;
 
void makeTracks(int N_Tracks, REveElement *trackHolder, float *pos)
{
   TRandom &r = *gRandom;
   auto prop = new REveTrackPropagator();
   prop->SetMagFieldObj(new REveMagFieldDuo(5, 3.5, -2.0));
   prop->SetMaxR(9);
   prop->SetMaxZ(15);
 
   double v = 0.1;
   double m = 5;
   for (int i = 0; i < N_Tracks; i++) {
      auto p = new TParticle();
      int pdg = 11 * (r.Integer(2) > 0 ? 1 : -1);
      p->SetPdgCode(pdg);
      p->SetProductionVertex(r.Uniform(-v, v) + pos[0], r.Uniform(-v, v) + pos[1], r.Uniform(-v, v) + pos[2], 1);
      p->SetMomentum(r.Uniform(-m, m), r.Uniform(-m, m), r.Uniform(-m, m) * r.Uniform(1, 3), 1);
      auto track = new REveTrack(p, 1, prop);
      track->MakeTrack();
      track->SetMainColor(kBlue);
      track->SetName(Form("Track_%d", i));
      trackHolder->AddElement(track);
   }
}
 
void makeProjected(REveElement *el, const char *pname, REveProjection::EPType_e t)
{
   auto eveMng = ROOT::Experimental::gEve;
   auto eventScene = eveMng->SpawnNewScene(Form("%s Event Data", pname), pname);
   auto mng = new REveProjectionManager();
   mng->SetProjection(t);
   mng->ImportElements(el, eventScene);
   auto view = eveMng->SpawnNewViewer(pname);
   view->SetCameraType(REveViewer::kCameraOrthoXOY);
   view->AddScene(eventScene);
}
 
void error_ellipse()
{
   auto eveMng = REveManager::Create();
   float pos[3] = {1.46589e-06, -1.30522e-05, -1.98267e-05};
 
   // symnetric matrix
   double a[16] = {1.46589e-01,
                   -1.30522e-02,
                   -1.98267e-02,
                   0,
                   -1.30522e-02,
                   4.22955e-02,
                   -5.86628e-03,
                   0,
                   -1.98267e-02,
                   -5.86628e-03,
                   2.12836e-01,
                   0,
                   0,
                   0,
                   0,
                   1};
 
   REveTrans t;
   t.SetFrom(a);
   TMatrixDSym xxx(3);
   for (int i = 0; i < 3; i++)
      for (int j = 0; j < 3; j++) {
         xxx(i, j) = t(i + 1, j + 1);
      }
 
   TMatrixDEigen eig(xxx);
   TVectorD xxxEig(eig.GetEigenValues());
   xxxEig = xxxEig.Sqrt();
 
   TMatrixD vecEig = eig.GetEigenVectors();
   REveVector v[3];
   int ei = 0;
   for (int i = 0; i < 3; ++i) {
      v[i].Set(vecEig(0, i), vecEig(1, i), vecEig(2, i));
      v[i] *= xxxEig(i);
   }
   REveElement *event = eveMng->GetEventScene();
 
   REveEllipsoid *ellipse = new REveEllipsoid("VertexError");
   ellipse->InitMainTrans();
   ellipse->SetMainColor(kGreen + 10);
   ellipse->SetLineWidth(2);
   ellipse->SetBaseVectors(v[0], v[1], v[2]);
   ellipse->Outline();
   event->AddElement(ellipse);
 
   // add TGeoSphere
   auto sph = new REveGeoShape("Sphere");
   sph->SetShape(new TGeoSphere(0.99f, 1.f));
   sph->SetMainColor(kGreen);
   sph->SetMainTransparency(80);
   sph->SetNSegments(80);
 
   float m0 = v[0].Mag();
   v[0].Normalize();
   float m1 = v[1].Mag();
   v[1].Normalize();
   float m2 = v[2].Mag();
   v[2].Normalize();
 
   sph->InitMainTrans();
   sph->RefMainTrans().SetBaseVec(1, v[0].fX, v[0].fY, v[0].fZ);
   sph->RefMainTrans().SetBaseVec(2, v[1].fX, v[1].fY, v[1].fZ);
   sph->RefMainTrans().SetBaseVec(3, v[2].fX, v[2].fY, v[2].fZ);
   sph->RefMainTrans().SetScale(m0, m1, m2);
   event->AddElement(sph);
 
   auto ps = new REvePointSet("Vertices");
   ps->SetMainColor(kYellow);
   ps->SetNextPoint(pos[0], pos[1], pos[2]);
   ps->SetMarkerStyle(4);
   ps->SetMarkerSize(4);
   float rng = 1;
   for (int i = 0; i < 6; ++i)
      ps->SetNextPoint(i * (rng / 3) - rng, 0, 0);
   event->AddElement(ps);
 
   auto trackHolder = new REveElement("Tracks");
   eveMng->GetEventScene()->AddElement(trackHolder);
   makeTracks(10, trackHolder, pos);
 
   makeProjected(eveMng->GetEventScene(), "RPhi", REveProjection::kPT_RPhi);
   makeProjected(eveMng->GetEventScene(), "RhoZ", REveProjection::kPT_RhoZ);
   eveMng->Show();
}