]>
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 | ||
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 | |
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); | |
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 | ||
180 | void 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 | /**************************************************************************/ | |
257 | void 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 | ||
63845220 | 267 | /**************************************************************************/ |
268 | void 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 | /**************************************************************************/ | |
293 | Double_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 | ||
e7aa1ed8 | 310 | /**************************************************************************/ |
311 | void 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 | } |