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