78942fdf6eb2b19e72ba7549f1d582e4a03137c4
[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 "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 }