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