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 "Riostream.h"
  11
  12 #include "AliEveCosmicRayFitter.h"
  13
  14 #include "TCanvas.h"
  15 #include "TPad.h"
  16 #include "TGraph.h"
  17 #include "TGraphErrors.h"
  18 #include "TGraph2D.h"
  19 #include "TGraph2DErrors.h"
  20
  21 #include "TVirtualFitter.h"
  22
  23 #include "TEveLine.h"
  24
  25 #include "AliLog.h"
  26
  27 #include <TEveManager.h>
  28 #include <Math/Vector3D.h>
  29
  30 //______________________________________________________________________
  31 // AliEveCosmicRayFitter
  32 //
  33 //  Same schema used in AliEveTrackFitter class.
  34 //
  35 //  AliEveCosmicRayFitter is a TEvePointSet interface to allow
  36 //  straight line fit. It creates a set of points by listening to
  37 //  selection signal of any TEvePointSet object. After selection, the
  38 //  list is feeded to TVirtualFitter class, which returns straight
  39 //  line parameters. The fit result is visualized with a TEveLine
  40 //  object.
  41 //
  42 //  Thanks to L. Moneta for the algorithm prepared to make a 3D
  43 //  straight line fit.
  44 //
  45 //  Author: A. De Caro
  46 //
  47
  48 Double_t Distance3D(Double_t x, Double_t y, Double_t z, Double_t *p);
  49 void SumDistance3D(Int_t &, Double_t *, Double_t & sum, Double_t * par, Int_t );
  50
  51 ClassImp(AliEveCosmicRayFitter)
  52
  53 AliEveCosmicRayFitter::AliEveCosmicRayFitter(const Text_t* name, Int_t n_points) :
  54     TEvePointSet   (name, n_points),
  55
  56     fLine3DFitter  (0),
  57
  58     fConnected     (kFALSE),
  59     fSPMap         (),
  60
  61     fGraphPicked1   (0),
  62     fGraphLinear1   (0),
  63     fGraphPicked2   (0),
  64     fGraphLinear2   (0),
  65     fGraphPicked3   (),
  66     fGraphLinear3   ()
  67
  68 {
  69   // Constructor.
  70
  71   fLine3DFitter = TVirtualFitter::Fitter(0, 4);
  72
  73   SetMarkerColor(3);
  74   SetOwnIds(kFALSE);
  75
  76   fGraphPicked1 = new TGraph("Y vs X");
  77   fGraphPicked1->SetName("Selected points");
  78   fGraphPicked1->SetMarkerColor(4);
  79   fGraphPicked1->SetMarkerStyle(4);
  80   fGraphPicked1->SetMarkerSize(2);
  81
  82   fGraphLinear1 = new TGraphErrors();
  83   fGraphLinear1->SetName("Fitted points");
  84   fGraphLinear1->SetMarkerColor(2);
  85
  86   fGraphPicked2 = new TGraph("Z vs X");
  87   fGraphPicked2->SetName("Selected points");
  88   fGraphPicked2->SetMarkerColor(4);
  89   fGraphPicked2->SetMarkerStyle(4);
  90   fGraphPicked2->SetMarkerSize(2);
  91
  92   fGraphLinear2 = new TGraphErrors();
  93   fGraphLinear2->SetName("Fitted points");
  94   fGraphLinear2->SetMarkerColor(2);
  95
  96   fGraphPicked3 = new TGraph2D("Z vs Y vs X");
  97   fGraphPicked3->SetName("Selected points");
  98   fGraphPicked3->SetMarkerColor(4);
  99   fGraphPicked3->SetMarkerStyle(4);
 100   fGraphPicked3->SetMarkerSize(2);
 101
 102   fGraphLinear3 = new TGraph2DErrors();
 103   fGraphLinear3->SetName("Fitted points");
 104   fGraphLinear3->SetMarkerColor(2);
 105
 106 }
 107
 108 AliEveCosmicRayFitter::~AliEveCosmicRayFitter()
 109 {
 110   // Destructor.
 111
 112   if(fLine3DFitter) delete fLine3DFitter;
 113
 114 }
 115
 116 /**************************************************************************/
 117 void AliEveCosmicRayFitter::Start()
 118 {
 119   // Clear existing point selection and maintain connection to the
 120   // TEvePointSet signal.
 121
 122   Reset();
 123   if(fConnected == kFALSE)
 124   {
 125     TQObject::Connect("TEvePointSet", "PointSelected(Int_t)",
 126                       "AliEveCosmicRayFitter", this, "AddFitPoint(Int_t)");
 127     fConnected = kTRUE;
 128   }
 129 }
 130
 131 /**************************************************************************/
 132 void AliEveCosmicRayFitter::Stop()
 133 {
 134   // Stop selection of points.
 135
 136   if(fConnected)
 137   {
 138     TQObject::Disconnect("TEvePointSet", "AddFitPoint(Int_t)");
 139     fConnected = kFALSE;
 140   }
 141 }
 142
 143 /**************************************************************************/
 144
 145 void AliEveCosmicRayFitter::Reset(Int_t n, Int_t ids)
 146 {
 147   // Reset selection.
 148
 149   if(fLine3DFitter) fLine3DFitter->Clear();
 150   TEvePointSet::Reset(n, ids);
 151   fSPMap.clear();
 152 }
 153
 154 /**************************************************************************/
 155
 156 void AliEveCosmicRayFitter::AddFitPoint(Int_t n)
 157 {
 158   // Add/remove given point depending if exists in the fSPMap.
 159
 160   Float_t x, y, z;
 161   TEvePointSet* ps = dynamic_cast<TEvePointSet*>((TQObject*) gTQSender);
 162
 163
 164   std::map<Point_t, Int_t>::iterator g = fSPMap.find(Point_t(ps, n));
 165   if(g != fSPMap.end())
 166   {
 167     Int_t idx = g->second;
 168     if(idx != fLastPoint)
 169     {
 170       GetPoint(fLastPoint, x, y, z);
 171       SetPoint(idx, x, y, z);
 172     }
 173     fSPMap.erase(g);
 174     fLastPoint--;
 175   }
 176   else
 177   {
 178     fSPMap[Point_t(ps, n)] = Size();
 179     ps->GetPoint(n, x, y, z);
 180     SetNextPoint(x, y, z);
 181     //  SetPointId(ps->GetPointId(n));
 182   }
 183   ResetBBox();
 184   ElementChanged(kTRUE, kTRUE);
 185 }
 186
 187 /**************************************************************************/
 188
 189 void AliEveCosmicRayFitter::FitTrack()
 190 {
 191   // Fit selected points with TLinearFitter fitter.
 192
 193   using namespace TMath;
 194
 195   if(fLine3DFitter) delete fLine3DFitter;
 196   fLine3DFitter = TVirtualFitter::Fitter(0, 4);
 197
 198   AliDebug(2, Form(" Selected point number %3i   %3i", fLastPoint+1, Size()));
 199
 200   Double_t *posXpoint = new Double_t[fLastPoint+1];
 201   Double_t *posYpoint = new Double_t[fLastPoint+1];
 202   Double_t *posZpoint = new Double_t[fLastPoint+1];
 203
 204   //Double_t *errXpoint = new Double_t[fLastPoint+1];
 205   //Double_t *errYpoint = new Double_t[fLastPoint+1];
 206   //Double_t *errZpoint = new Double_t[fLastPoint+1];
 207
 208   Double_t x, y, z;
 209
 210   for (Int_t i=0; i<=fLastPoint; i++) {
 211     x=0., y=0., z=0.;
 212     GetPoint(i, x, y, z);
 213     posXpoint[i] = x;
 214     posYpoint[i] = y;
 215     posZpoint[i] = z;
 216     //errXpoint[i] = 0.001*posXpoint[i];
 217     //errYpoint[i] = 0.05*posYpoint[i];
 218     //errZpoint[i] = 0.05*posZpoint[i];
 219
 220     fGraphPicked3->SetPoint(i,x,y,z);
 221     //fGraphPicked3->SetPointError(i,errXpoint[i],errYpoint[i],errZpoint[i]);
 222
 223     AliDebug(2, Form("  %f    %f    %f",
 224                      posXpoint[i], posYpoint[i], posZpoint[i]));
 225   }
 226
 227   fLine3DFitter->SetObjectFit(fGraphPicked3);
 228   fLine3DFitter->SetFCN(SumDistance3D);
 229   Double_t arglist[10];
 230   arglist[0] = 3;
 231   fLine3DFitter->ExecuteCommand("SET PRINT",arglist,1);
 232
 233   Double_t pStart[4] = {1,1,1,1};
 234   fLine3DFitter->SetParameter(0,"y0",pStart[0],0.01,0,0);
 235   fLine3DFitter->SetParameter(1,"Ay",pStart[1],0.01,0,0);
 236   fLine3DFitter->SetParameter(2,"z0",pStart[2],0.01,0,0);
 237   fLine3DFitter->SetParameter(3,"Az",pStart[3],0.01,0,0);
 238
 239   arglist[0] = 1000; // number of function calls
 240   arglist[1] = 0.001; // tolerance
 241   fLine3DFitter->ExecuteCommand("MIGRAD",arglist,2);
 242
 243   //if (minos) fLine3DFitter->ExecuteCommand("MINOS",arglist,0);
 244   Int_t nvpar,nparx;
 245   Double_t amin,edm, errdef;
 246   fLine3DFitter->GetStats(amin,edm,errdef,nvpar,nparx);
 247   fLine3DFitter->PrintResults(1,amin);
 248
 249   TEveLine* line3D = new TEveLine();
 250   line3D->SetLineColor(2);
 251   line3D->SetLineWidth(2);
 252   line3D->SetLineStyle(2);
 253
 254   Double_t parFit3D[4];
 255   for (int i = 0; i <4; ++i)
 256     parFit3D[i] = fLine3DFitter->GetParameter(i);
 257   Float_t xCoor = 16000.;
 258   Float_t yCoorP = parFit3D[0] + xCoor*parFit3D[1];
 259   Float_t zCoorP = parFit3D[2] + xCoor*parFit3D[3];
 260   Float_t yCoorM = parFit3D[0] - xCoor*parFit3D[1];
 261   Float_t zCoorM = parFit3D[2] - xCoor*parFit3D[3];
 262
 263   line3D->SetNextPoint( xCoor, yCoorP, zCoorP);
 264   line3D->SetNextPoint(-xCoor, yCoorM, zCoorM);
 265
 266   gEve->AddElement(line3D);
 267
 268   delete posXpoint;
 269   delete posYpoint;
 270   delete posZpoint;
 271   //delete errXpoint;
 272   //delete errYpoint;
 273   //delete errZpoint;
 274
 275 }
 276
 277 /**************************************************************************/
 278 void AliEveCosmicRayFitter::DestroyElements()
 279 {
 280   // Virtual method of base class Reve::RenderElement.
 281   // It preserves track list to have coomon track propagator attributes.
 282
 283   TEveElement::DestroyElements();
 284   UpdateItems();
 285 }
 286
 287 /**************************************************************************/
 288 void AliEveCosmicRayFitter::DrawDebugGraph()
 289 {
 290   //
 291   // Draw graph of Line fit.
 292   //
 293
 294   static const TEveException eH("CosmicRayFitter::DrawDebugGraph ");
 295
 296   if(fLine3DFitter == 0)
 297     throw(eH + "fitter not set.");
 298
 299   AliDebug(2,Form("           %3i   %3i\n", fLastPoint+1, Size()));
 300
 301   Int_t nR = fLastPoint+1;
 302   fGraphPicked1->Set(nR);
 303   fGraphLinear1->Set(nR);
 304   fGraphPicked2->Set(nR);
 305   fGraphLinear2->Set(nR);
 306   fGraphLinear3->Set(nR);
 307
 308   Double_t *x = new Double_t[nR];
 309   Double_t *y = new Double_t[nR];
 310   Double_t *z = new Double_t[nR];
 311   Float_t xx=0., yy=0., zz=0.;
 312   for (Int_t i=0; i<nR; i++) {
 313     GetPoint(i, xx, yy, zz);
 314     x[i] = (Double_t)xx;
 315     y[i] = (Double_t)yy;
 316     z[i] = (Double_t)zz;
 317     AliDebug(2,Form("  %f      %f      %f\n", x[i], y[i], z[i]));
 318   }
 319
 320   Double_t parFit3D[4];
 321    for (int i = 0; i <4; ++i)
 322       parFit3D[i] = fLine3DFitter->GetParameter(i);
 323
 324   Double_t yCoor = 0.;
 325   Double_t zCoor = 0.;
 326   for (Int_t i=0; i<nR; i++)
 327   {
 328     yCoor = parFit3D[0] + x[i]*parFit3D[1];
 329     zCoor = parFit3D[2] + x[i]*parFit3D[3];
 330     fGraphPicked1->SetPoint(i, x[i], y[i]);
 331     fGraphLinear1->SetPoint(i, x[i], yCoor);
 332     //fGraphLinear1->SetPointError(i, TMath::Abs(x[i])*0.001, TMath::Abs(yCoor)*0.05);
 333     fGraphPicked2->SetPoint(i, x[i], z[i]);
 334     fGraphLinear2->SetPoint(i, x[i], zCoor);
 335     //fGraphLinear2->SetPointError(i, TMath::Abs(x[i])*0.001, TMath::Abs(zCoor)*0.05);
 336
 337     fGraphLinear3->SetPoint(i, x[i], yCoor, zCoor);
 338     //fGraphLinear3->SetPointError(i, TMath::Abs(x[i])*0.001, TMath::Abs(yCoor)*0.05, TMath::Abs(zCoor)*0.05);
 339   }
 340
 341   TCanvas * canvas=0;
 342   if (gPad) gPad->Clear();
 343   else if (gPad==0 || gPad->GetCanvas()->IsEditable() == kFALSE) {
 344     canvas = new TCanvas("canvas", "CosmicRayFitter", 800, 800);
 345     canvas->Clear();
 346   }
 347   canvas->SetHighLightColor(2);
 348   canvas->Range(0,0,1,1);
 349   canvas->SetFillColor(0);
 350   canvas->SetBorderMode(0);
 351   canvas->SetBorderSize(2);
 352
 353   canvas->Divide(1,2);
 354
 355   TPad *canvas_1 = new TPad("canvas_1", "canvas_1",0.01,0.01,0.49,0.49);
 356   canvas_1->Draw();
 357   canvas_1->cd();
 358   canvas_1->SetFillColor(0);
 359   canvas_1->SetBorderSize(2);
 360   canvas_1->SetFrameBorderMode(0);
 361   fGraphPicked1->Draw("AP");
 362   fGraphLinear1->Draw("SAME P");
 363
 364   canvas_1->Modified();
 365   canvas->cd();
 366
 367   TPad *canvas_2 = new TPad("canvas_2", "canvas_2",0.51,0.01,0.99,0.49);
 368   canvas_2->Draw();
 369   canvas_2->cd();
 370   canvas_2->SetFillColor(0);
 371   canvas_2->SetBorderSize(2);
 372   canvas_2->SetFrameBorderMode(0);
 373   fGraphPicked2->Draw("AP");
 374   fGraphLinear2->Draw("SAME P");
 375
 376   canvas_2->Modified();
 377   canvas->cd();
 378
 379   TPad *canvas_3 = new TPad("canvas_2", "canvas_2",0.01,0.51,0.99,0.99);
 380   canvas_3->Draw();
 381   canvas_3->cd();
 382   canvas_3->SetFillColor(0);
 383   canvas_3->SetBorderSize(2);
 384   canvas_3->SetFrameBorderMode(0);
 385   fGraphPicked3->Draw("AP");
 386   fGraphLinear3->Draw("SAME P");
 387
 388   canvas_3->Modified();
 389   canvas->cd();
 390
 391   canvas->Modified();
 392   canvas->cd();
 393
 394   delete x;
 395   delete y;
 396   delete z;
 397
 398 }
 399
 400 /**************************************************************************/
 401 void SumDistance3D(Int_t &, Double_t *, Double_t & sum, Double_t * par, Int_t )
 402 {
 403   //
 404   // Square sum of distances
 405   //
 406
 407   TGraph2D * gr =
 408     dynamic_cast<TGraph2D*>((TVirtualFitter::GetFitter())->GetObjectFit() );
 409   assert(gr != 0);
 410
 411   Double_t * x = gr->GetX();
 412   Double_t * y = gr->GetY();
 413   Double_t * z = gr->GetZ();
 414   Int_t npoints = gr->GetN();
 415   sum = 0;
 416   Double_t d = 0;
 417   for (Int_t i  = 0; i < npoints; ++i) {
 418     d = Distance3D(x[i],y[i],z[i],par);
 419     sum += d;
 420 #ifdef DEBUG
 421     std::cout << " point " << i << "\t"
 422               << x[i] << "\t"
 423               << y[i] << "\t"
 424               << z[i] << "\t"
 425               << std::sqrt(d) << std::endl;
 426     std::cout << " Total sum2 = " << sum << std::endl;
 427 #endif
 428
 429   }
 430
 431 }
 432
 433 /**************************************************************************/
 434 Double_t Distance3D(Double_t x, Double_t y, Double_t z, Double_t *p)
 435 {
 436   //
 437   // distance line point is D = | (xp-x0) cross  ux |
 438   // where ux is direction of line and x0 is a point in the line (like t = 0)
 439   //
 440
 441   using namespace ROOT::Math;
 442
 443   XYZVector xp(x,y,z);
 444   XYZVector x0(0., p[0], p[2]);
 445   XYZVector x1(1., p[0] + p[1], p[2] + p[3]);
 446   XYZVector u = (x1-x0).Unit();
 447   Double_t d2 = ((xp-x0).Cross(u)) .Mag2();
 448   return d2;
 449 }
 450
 451 /**************************************************************************/