8bad243b583bdeb688a8d8ba913eee148aa788db
[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 "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 /**************************************************************************/