EVE/EveBase/AliEveCosmicRayFitter.cxx

   1 // $Id$
   2 // Main authors: Matevz Tadel & Alja Mrak-Tadel: 2006, 2007
   3
   4 /**************************************************************************
   5  * Copyright(c) 1998-2008, ALICE Experiment at CERN, all rights reserved. *
   6  * See http://aliceinfo.cern.ch/Offline/AliRoot/License.html for          *
   7  * full copyright notice.                                                 *
   8  **************************************************************************/
   9
  10 #include "AliEveCosmicRayFitter.h"
  11
  12 #include "AliLog.h"
  13
  14 #include "TEveTrack.h"
  15 #include "TEveTrackPropagator.h"
  16 #include "TEveVSDStructs.h"
  17 #include "TEveManager.h"
  18
  19 #include "TCanvas.h"
  20 #include "TPad.h"
  21 #include "TGraph.h"
  22 #include "TGraph2D.h"
  23 #include "TGraphErrors.h"
  24 #include "TLinearFitter.h"
  25 #include <Math/Vector3D.h>
  26
  27 //______________________________________________________________________
  28 // AliEveCosmicRayFitter
  29 //
  30 //  Same schema used in AliEveTrackFitter class.
  31 //
  32 //  AliEveCosmicRayFitter is a TEvePointSet interface to allow
  33 //  straight line fit. It creates a set of points by listening to
  34 //  selection signal of any TEvePointSet object. After selection, the
  35 //  list is feeded to TVirtualFitter class, which returns straight
  36 //  line parameters. The fit result is visualized with a TEveLine
  37 //  object.
  38 //
  39 //  Thanks to L. Moneta for the algorithm prepared to make a 3D
  40 //  straight line fit.
  41 //
  42 //  Author: A. De Caro
  43 //
  44
  45
  46 ClassImp(AliEveCosmicRayFitter)
  47
  48 AliEveCosmicRayFitter::AliEveCosmicRayFitter(const Text_t* name, Int_t n_points) :
  49     TEvePointSet    (name, n_points),
  50
  51     fConnected      (kFALSE),
  52     fSPMap          (),
  53     fTrackList      (0),
  54
  55     fGraphPicked1   (0),
  56     fGraphLinear1   (0),
  57     fGraphPicked2   (0),
  58     fGraphLinear2   (0)
  59 {
  60   // Constructor.
  61
  62   SetMarkerColor(3);
  63   SetOwnIds(kFALSE);
  64
  65   fTrackList = new TEveTrackList("Cosmic ray");
  66   fTrackList->SetTitle("muons");
  67   fTrackList->SetMainColor((Color_t)8);
  68   gEve->AddElement(fTrackList);
  69   UpdateItems();
  70
  71   fGraphPicked1 = new TGraph();
  72   fGraphPicked1->SetName("Selected points");
  73   fGraphPicked1->SetMarkerColor(4);
  74   fGraphPicked1->SetMarkerStyle(4);
  75   fGraphPicked1->SetMarkerSize(2);
  76
  77   fGraphLinear1 = new TGraphErrors();
  78   fGraphLinear1->SetName("Fitted points");
  79   fGraphLinear1->SetMarkerColor(2);
  80
  81   fGraphPicked2 = new TGraph();
  82   fGraphPicked2->SetName("Selected points");
  83   fGraphPicked2->SetMarkerColor(4);
  84   fGraphPicked2->SetMarkerStyle(4);
  85   fGraphPicked2->SetMarkerSize(2);
  86
  87   fGraphLinear2 = new TGraphErrors();
  88   fGraphLinear2->SetName("Fitted points");
  89   fGraphLinear2->SetMarkerColor(2);
  90 }
  91
  92 AliEveCosmicRayFitter::~AliEveCosmicRayFitter()
  93 {
  94   // Destructor.
  95
  96   if(fLineFitter1) delete fLineFitter1;
  97   if(fLineFitter2) delete fLineFitter2;
  98
  99   fTrackList->DecDenyDestroy();
 100   delete fTrackList;
 101 }
 102
 103 /**************************************************************************/
 104 void AliEveCosmicRayFitter::Start()
 105 {
 106   // Clear existing point selection and maintain connection to the
 107   // TEvePointSet signal.
 108
 109   Reset();
 110   if(fConnected == kFALSE)
 111   {
 112     TQObject::Connect("TEvePointSet", "PointSelected(Int_t)",
 113                       "AliEveCosmicRayFitter", this, "AddFitPoint(Int_t)");
 114     fConnected = kTRUE;
 115   }
 116 }
 117
 118 /**************************************************************************/
 119 void AliEveCosmicRayFitter::Stop()
 120 {
 121   // Stop selection of points.
 122
 123   if(fConnected)
 124   {
 125     TQObject::Disconnect("TEvePointSet", "AddFitPoint(Int_t)");
 126     fConnected = kFALSE;
 127   }
 128 }
 129
 130 /**************************************************************************/
 131
 132 void AliEveCosmicRayFitter::Reset(Int_t /*n*/, Int_t /*ids*/)
 133 {
 134   // Reset selection.
 135
 136   TEvePointSet::Reset();
 137   fSPMap.clear();
 138
 139   if(fLineFitter1) fLineFitter1->Clear();
 140   if(fLineFitter2) fLineFitter2->Clear();
 141 }
 142
 143 /**************************************************************************/
 144
 145 void AliEveCosmicRayFitter::AddFitPoint(Int_t n)
 146 {
 147   // Add/remove given point depending if exists in the fSPMap.
 148
 149   Float_t x, y, z;
 150   TEvePointSet* ps = dynamic_cast<TEvePointSet*>((TQObject*) gTQSender);
 151
 152   std::map<Point_t, Int_t>::iterator g = fSPMap.find(Point_t(ps, n));
 153   if(g != fSPMap.end())
 154   {
 155     Int_t idx = g->second;
 156     if(idx != fLastPoint)
 157     {
 158       GetPoint(fLastPoint, x, y, z);
 159       SetPoint(idx, x, y, z);
 160     }
 161     fSPMap.erase(g);
 162     fLastPoint--;
 163   }
 164   else
 165   {
 166     fSPMap[Point_t(ps, n)] = Size();
 167     ps->GetPoint(n, x, y, z);
 168     SetNextPoint(x, y, z);
 169   }
 170
 171   ComputeBBox();
 172   ElementChanged(kTRUE, kTRUE);
 173 }
 174
 175 /**************************************************************************/
 176
 177 void AliEveCosmicRayFitter::FitTrack()
 178 {
 179   // Fit selected points with TLinearFitter fitter.
 180
 181   static const TEveException eH("CosmicRayFitter::FitTrack ");
 182
 183   if(fLineFitter1) delete fLineFitter1;
 184   if(fLineFitter2) delete fLineFitter2;
 185   fLineFitter1 = new TLinearFitter(1, "hyp1");
 186   fLineFitter2 = new TLinearFitter(1, "hyp1");
 187
 188   AliDebug(2, Form(" Selected %3d points:\n", fLastPoint+1));
 189
 190   Double_t *posXpoint = new Double_t[fLastPoint+1];
 191   Double_t *posYpoint = new Double_t[fLastPoint+1];
 192   Double_t *posZpoint = new Double_t[fLastPoint+1];
 193   Double_t *errYpoint = new Double_t[fLastPoint+1];
 194   Double_t *errZpoint = new Double_t[fLastPoint+1];
 195   Double_t x, y, z;
 196   for (Int_t i=0; i<=fLastPoint; i++) {
 197     x=0., y=0., z=0.;
 198     GetPoint(i, x, y, z);
 199     posXpoint[i] = x;
 200     posYpoint[i] = y;
 201     posZpoint[i] = z;
 202     errYpoint[i] = 10.;//0.05*posYpoint[i];//arbitrary
 203     errZpoint[i] = 10.;//0.05*posZpoint[i];//arbitrary
 204
 205     AliDebug(2, Form("(%f, %f, %f) \n", posXpoint[i], posYpoint[i], posZpoint[i]));
 206   }
 207
 208   fLineFitter1->AssignData(fLastPoint+1, 1, posXpoint, posYpoint, errYpoint);
 209   fLineFitter2->AssignData(fLastPoint+1, 1, posXpoint, posZpoint, errZpoint);
 210   Int_t fitResult1 = fLineFitter1->Eval();
 211   Int_t fitResult2 = fLineFitter2->Eval();
 212
 213   if (fitResult1 || fitResult2)
 214     throw(eH + Form(" linear fit result is not ok (%1i, %1i)", fitResult1, fitResult2));
 215
 216
 217   delete [] posXpoint;
 218   delete [] posYpoint;
 219   delete [] posZpoint;
 220   delete [] errYpoint;
 221   delete [] errZpoint;
 222
 223   // make track
 224   TVectorD params1;
 225   TVectorD params2;
 226   fLineFitter1->GetParameters(params1);
 227   fLineFitter2->GetParameters(params2);
 228
 229   TEveRecTrack rc;
 230   GetPoint(0, x, y, z);
 231   y = params1(0) + x*params1(1);
 232   z = params2(0) + x*params2(1);
 233   rc.fV.Set(x, y, z);
 234   rc.fP.Set(1, params1(1), params2(1));
 235
 236   TEveTrack* track = new TEveTrack(&rc, fTrackList->GetPropagator());
 237   track->SetAttLineAttMarker(fTrackList);
 238
 239   TEveTrackPropagator* tp = fTrackList->GetPropagator();
 240
 241   tp->InitTrack(rc.fV, rc.fP, 1, 0);
 242   tp->GoToBounds(rc.fP);
 243
 244   rc.fP.Set(-1, -params1(1), -params2(1));
 245   tp->InitTrack(rc.fV, rc.fP, 1, 0);
 246   tp->GoToBounds(rc.fP);
 247   tp->FillPointSet(track);
 248   tp->ResetTrack();
 249
 250   fTrackList->AddElement(track);
 251 }
 252
 253 /**************************************************************************/
 254 void AliEveCosmicRayFitter::DestroyElements()
 255 {
 256   // Virtual method of base class Reve::RenderElement.
 257   // It preserves track list to have coomon track propagator attributes.
 258
 259   TEveElement::DestroyElements();
 260   gEve->AddElement(fTrackList, this);
 261   fTrackList->DestroyElements();
 262   UpdateItems();
 263 }
 264
 265 /**************************************************************************/
 266 void AliEveCosmicRayFitter::SumDistance3D(Int_t &, Double_t *, Double_t & sum, Double_t * par, Int_t )
 267 {
 268   //
 269   // Square sum of distances
 270   //
 271
 272   TGraph2D * gr =
 273     dynamic_cast<TGraph2D*>((TVirtualFitter::GetFitter())->GetObjectFit() );
 274   assert(gr != 0);
 275
 276   Double_t * x = gr->GetX();
 277   Double_t * y = gr->GetY();
 278   Double_t * z = gr->GetZ();
 279   Int_t npoints = gr->GetN();
 280   sum = 0;
 281   Double_t d = 0;
 282   for (Int_t i  = 0; i < npoints; ++i) {
 283     d = Distance3D(x[i],y[i],z[i],par);
 284     sum += d;
 285     // printf("%d (%f, %f, %f) dist:%d", i,  x[i], y[i], z[i], TMath::Sqrt(d));
 286   }
 287   // printf("Total sum %f\n", sum);
 288 }
 289
 290 /**************************************************************************/
 291 Double_t AliEveCosmicRayFitter::Distance3D(Double_t x, Double_t y, Double_t z, Double_t *p)
 292 {
 293   //
 294   // distance line point is D = | (xp-x0) cross  ux |
 295   // where ux is direction of line and x0 is a point in the line (like t = 0)
 296   //
 297
 298   using namespace ROOT::Math;
 299
 300   XYZVector xp(x,y,z);
 301   XYZVector x0(0., p[0], p[2]);
 302   XYZVector x1(1., p[0] + p[1], p[2] + p[3]);
 303   XYZVector u = (x1-x0).Unit();
 304   Double_t d2 = ((xp-x0).Cross(u)) .Mag2();
 305   return d2;
 306 }
 307
 308 /**************************************************************************/
 309 void AliEveCosmicRayFitter::DrawDebugGraph()
 310 {
 311   //
 312   // Draw graph of Line fit.
 313   //
 314
 315   static const TEveException eH("CosmicRayFitter::DrawDebugGraph ");
 316
 317   if(fLineFitter1 == 0 || fLineFitter2 == 0)
 318     throw(eH + "fitter not set.");
 319
 320   TVectorD params1;
 321   TVectorD params2;
 322   fLineFitter1->GetParameters(params1);
 323   fLineFitter2->GetParameters(params2);
 324
 325   // fill graphs
 326
 327   Int_t nR = fLastPoint+1;
 328   Double_t *x = new Double_t[nR];
 329   Double_t *y = new Double_t[nR];
 330   Double_t *z = new Double_t[nR];
 331   Float_t xx=0., yy=0., zz=0.;
 332   for (Int_t i=0; i<nR; i++) {
 333     GetPoint(i, xx, yy, zz);
 334     x[i] = (Double_t)xx;
 335     y[i] = (Double_t)yy;
 336     z[i] = (Double_t)zz;
 337     AliDebug(2, Form("%d  (%f, %f, %f)\n", i, x[i], y[i], z[i]));
 338   }
 339
 340   fGraphPicked1->Set(nR);
 341   fGraphLinear1->Set(nR);
 342   fGraphPicked2->Set(nR);
 343   fGraphLinear2->Set(nR);
 344   Double_t yCoor = 0.;
 345   Double_t zCoor = 0.;
 346   for (Int_t i=0; i<nR; i++)
 347   {
 348     yCoor = params1(0) + x[i]*params1(1);
 349     zCoor = params2(0) + x[i]*params2(1);
 350     fGraphPicked1->SetPoint(i, x[i], y[i]);
 351     fGraphLinear1->SetPoint(i, x[i], yCoor);
 352     fGraphLinear1->SetPointError(i, TMath::Abs(x[i])*0.001, TMath::Abs(yCoor)*0.05);
 353     fGraphPicked2->SetPoint(i, x[i], z[i]);
 354     fGraphLinear2->SetPoint(i, x[i], zCoor);
 355     fGraphLinear2->SetPointError(i, TMath::Abs(x[i])*0.001, TMath::Abs(zCoor)*0.05);
 356   }
 357   delete [] x;
 358   delete [] y;
 359   delete [] z;
 360
 361
 362   // draw graphs
 363
 364   TCanvas * canvas = 0;
 365   if (gPad) gPad->Clear();
 366   else if (gPad==0 || gPad->GetCanvas()->IsEditable() == kFALSE) {
 367     canvas = new TCanvas("canvas", "CosmicRayFitter", 800, 400);
 368     canvas->Clear();
 369   }
 370
 371   canvas->SetHighLightColor(2);
 372   canvas->Range(0,0,1,1);
 373   canvas->SetFillColor(0);
 374   canvas->SetBorderMode(0);
 375   canvas->SetBorderSize(2);
 376
 377   TPad *canvas_1 = new TPad("canvas_1", "canvas_1",0.01,0.01,0.49,0.99);
 378   canvas_1->Draw();
 379   canvas_1->cd();
 380   canvas_1->SetFillColor(0);
 381   canvas_1->SetBorderSize(2);
 382   canvas_1->SetFrameBorderMode(0);
 383   fGraphPicked1->Draw("AP");
 384   fGraphLinear1->Draw("SAME P");
 385   canvas_1->Modified();
 386   canvas->cd();
 387
 388   TPad *canvas_2 = new TPad("canvas_2", "canvas_2",0.51,0.01,0.99,0.99);
 389   canvas_2->Draw();
 390   canvas_2->cd();
 391   canvas_2->SetFillColor(0);
 392   canvas_2->SetBorderSize(2);
 393   canvas_2->SetFrameBorderMode(0);
 394   fGraphPicked2->Draw("AP");
 395   fGraphLinear2->Draw("SAME P");
 396   canvas_2->Modified();
 397   canvas->cd();
 398 }