]>
Commit | Line | Data |
---|---|---|
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 | ||
14 | #include "TEveTrack.h" | |
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 | |
46 | ClassImp(AliEveCosmicRayFitter) | |
47 | ||
48 | AliEveCosmicRayFitter::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 | ||
94 | AliEveCosmicRayFitter::~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 | /**************************************************************************/ | |
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 | ||
63845220 | 134 | void 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 | ||
147 | void 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); | |
153 | ||
e7aa1ed8 | 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 | } | |
63845220 | 166 | else |
e7aa1ed8 | 167 | { |
168 | fSPMap[Point_t(ps, n)] = Size(); | |
169 | ps->GetPoint(n, x, y, z); | |
63845220 | 170 | SetNextPoint(x, y, z); |
e7aa1ed8 | 171 | } |
63845220 | 172 | |
173 | ComputeBBox(); | |
e7aa1ed8 | 174 | ElementChanged(kTRUE, kTRUE); |
175 | } | |
176 | ||
177 | /**************************************************************************/ | |
178 | ||
179 | void AliEveCosmicRayFitter::FitTrack() | |
180 | { | |
181 | // Fit selected points with TLinearFitter fitter. | |
182 | ||
63845220 | 183 | static const TEveException eH("CosmicRayFitter::FitTrack "); |
e7aa1ed8 | 184 | |
3dc80ba2 | 185 | if(fLineFitter1) delete fLineFitter1; |
186 | if(fLineFitter2) delete fLineFitter2; | |
3dc80ba2 | 187 | fLineFitter1 = new TLinearFitter(1, "hyp1"); |
188 | fLineFitter2 = new TLinearFitter(1, "hyp1"); | |
e7aa1ed8 | 189 | |
63845220 | 190 | AliDebug(2, Form(" Selected %3d points:\n", fLastPoint+1)); |
e7aa1ed8 | 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]; | |
3dc80ba2 | 195 | Double_t *errYpoint = new Double_t[fLastPoint+1]; |
196 | Double_t *errZpoint = new Double_t[fLastPoint+1]; | |
e7aa1ed8 | 197 | Double_t x, y, z; |
63845220 | 198 | for (Int_t i=0; i<=fLastPoint; i++) { |
e7aa1ed8 | 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; | |
3dc80ba2 | 204 | errYpoint[i] = 10.;//0.05*posYpoint[i];//arbitrary |
205 | errZpoint[i] = 10.;//0.05*posZpoint[i];//arbitrary | |
e7aa1ed8 | 206 | |
63845220 | 207 | AliDebug(2, Form("(%f, %f, %f) \n", posXpoint[i], posYpoint[i], posZpoint[i])); |
e7aa1ed8 | 208 | } |
e7aa1ed8 | 209 | |
3dc80ba2 | 210 | fLineFitter1->AssignData(fLastPoint+1, 1, posXpoint, posYpoint, errYpoint); |
211 | fLineFitter2->AssignData(fLastPoint+1, 1, posXpoint, posZpoint, errZpoint); | |
3dc80ba2 | 212 | Int_t fitResult1 = fLineFitter1->Eval(); |
213 | Int_t fitResult2 = fLineFitter2->Eval(); | |
63845220 | 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; | |
3dc80ba2 | 224 | |
225 | // make track | |
226 | TVectorD params1; | |
227 | TVectorD params2; | |
228 | fLineFitter1->GetParameters(params1); | |
229 | fLineFitter2->GetParameters(params2); | |
230 | ||
63845220 | 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)); | |
3dc80ba2 | 237 | |
63845220 | 238 | TEveTrack* track = new TEveTrack(&rc, fTrackList->GetPropagator()); |
239 | track->SetAttLineAttMarker(fTrackList); | |
3dc80ba2 | 240 | |
63845220 | 241 | TEveTrackPropagator* tp = fTrackList->GetPropagator(); |
3dc80ba2 | 242 | |
68ca2fe7 | 243 | tp->InitTrack(rc.fV, 0); |
63845220 | 244 | tp->GoToBounds(rc.fP); |
3dc80ba2 | 245 | |
63845220 | 246 | rc.fP.Set(-1, -params1(1), -params2(1)); |
68ca2fe7 | 247 | tp->InitTrack(rc.fV, 0); |
63845220 | 248 | tp->GoToBounds(rc.fP); |
249 | tp->FillPointSet(track); | |
250 | tp->ResetTrack(); | |
e7aa1ed8 | 251 | |
63845220 | 252 | fTrackList->AddElement(track); |
e7aa1ed8 | 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(); | |
3dc80ba2 | 262 | gEve->AddElement(fTrackList, this); |
263 | fTrackList->DestroyElements(); | |
e7aa1ed8 | 264 | } |
265 | ||
63845220 | 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 | ||
e7aa1ed8 | 309 | /**************************************************************************/ |
310 | void AliEveCosmicRayFitter::DrawDebugGraph() | |
311 | { | |
312 | // | |
313 | // Draw graph of Line fit. | |
314 | // | |
315 | ||
316 | static const TEveException eH("CosmicRayFitter::DrawDebugGraph "); | |
317 | ||
3dc80ba2 | 318 | if(fLineFitter1 == 0 || fLineFitter2 == 0) |
e7aa1ed8 | 319 | throw(eH + "fitter not set."); |
320 | ||
3dc80ba2 | 321 | TVectorD params1; |
322 | TVectorD params2; | |
323 | fLineFitter1->GetParameters(params1); | |
324 | fLineFitter2->GetParameters(params2); | |
63845220 | 325 | |
326 | // fill graphs | |
327 | ||
328 | Int_t nR = fLastPoint+1; | |
e7aa1ed8 | 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; | |
63845220 | 338 | AliDebug(2, Form("%d (%f, %f, %f)\n", i, x[i], y[i], z[i])); |
e7aa1ed8 | 339 | } |
340 | ||
63845220 | 341 | fGraphPicked1->Set(nR); |
342 | fGraphLinear1->Set(nR); | |
343 | fGraphPicked2->Set(nR); | |
344 | fGraphLinear2->Set(nR); | |
e7aa1ed8 | 345 | Double_t yCoor = 0.; |
346 | Double_t zCoor = 0.; | |
347 | for (Int_t i=0; i<nR; i++) | |
348 | { | |
3dc80ba2 | 349 | yCoor = params1(0) + x[i]*params1(1); |
350 | zCoor = params2(0) + x[i]*params2(1); | |
e7aa1ed8 | 351 | fGraphPicked1->SetPoint(i, x[i], y[i]); |
352 | fGraphLinear1->SetPoint(i, x[i], yCoor); | |
3dc80ba2 | 353 | fGraphLinear1->SetPointError(i, TMath::Abs(x[i])*0.001, TMath::Abs(yCoor)*0.05); |
e7aa1ed8 | 354 | fGraphPicked2->SetPoint(i, x[i], z[i]); |
355 | fGraphLinear2->SetPoint(i, x[i], zCoor); | |
3dc80ba2 | 356 | fGraphLinear2->SetPointError(i, TMath::Abs(x[i])*0.001, TMath::Abs(zCoor)*0.05); |
e7aa1ed8 | 357 | } |
63845220 | 358 | delete [] x; |
359 | delete [] y; | |
360 | delete [] z; | |
361 | ||
362 | ||
363 | // draw graphs | |
364 | ||
3dc80ba2 | 365 | TCanvas * canvas = 0; |
e7aa1ed8 | 366 | if (gPad) gPad->Clear(); |
367 | else if (gPad==0 || gPad->GetCanvas()->IsEditable() == kFALSE) { | |
63845220 | 368 | canvas = new TCanvas("canvas", "CosmicRayFitter", 800, 400); |
e7aa1ed8 | 369 | canvas->Clear(); |
370 | } | |
63845220 | 371 | |
e7aa1ed8 | 372 | canvas->SetHighLightColor(2); |
373 | canvas->Range(0,0,1,1); | |
374 | canvas->SetFillColor(0); | |
375 | canvas->SetBorderMode(0); | |
376 | canvas->SetBorderSize(2); | |
377 | ||
63845220 | 378 | TPad *canvas_1 = new TPad("canvas_1", "canvas_1",0.01,0.01,0.49,0.99); |
e7aa1ed8 | 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"); | |
e7aa1ed8 | 386 | canvas_1->Modified(); |
387 | canvas->cd(); | |
388 | ||
63845220 | 389 | TPad *canvas_2 = new TPad("canvas_2", "canvas_2",0.51,0.01,0.99,0.99); |
e7aa1ed8 | 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"); | |
e7aa1ed8 | 397 | canvas_2->Modified(); |
398 | canvas->cd(); | |
e7aa1ed8 | 399 | } |