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