1f26c323 |
1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * |
3 | * * |
4 | * Author: Alberto Pulvirenti. * |
5 | * * |
6 | * Permission to use, copy, modify and distribute this software and its * |
7 | * documentation strictly for non-commercial purposes is hereby granted * |
8 | * without fee, provided that the above copyright notice appears in all * |
9 | * copies and that both the copyright notice and this permission notice * |
10 | * appear in the supporting documentation. The authors make no claims * |
11 | * about the suitability of this software for any purpose. It is * |
12 | * provided "as is" without express or implied warranty. * |
13 | * * |
14 | * track finder for ITS stand-alone with neural network algorithm * |
15 | * This class defines a Hopfield MFT neural network simulation * |
16 | * which reads all recpoints from an event and produces a tree with * |
17 | * the points belonging to recognized tracks * |
18 | * the TTree obtained as the output file will be saved in a .root file * |
19 | * and it must be read by groups of six entries at a time * |
20 | **************************************************************************/ |
21 | |
4ae5bbc4 |
22 | #include <Riostream.h> |
23 | #include <Riostream.h> |
1f26c323 |
24 | #include <stdlib.h> |
25 | |
26 | #include <TROOT.h> |
27 | #include <TFile.h> |
28 | #include <TObjArray.h> |
29 | #include <TTree.h> |
30 | #include <TMath.h> |
31 | #include <TRandom.h> |
32 | #include <TVector3.h> |
33 | |
34 | #include "AliRun.h" |
35 | #include "AliITS.h" |
36 | #include "AliITSgeom.h" |
37 | #include "AliITSgeomMatrix.h" |
38 | #include "AliITSRecPoint.h" |
39 | #include "AliITSglobalRecPoint.h" |
40 | #include "AliITSneuralTrack.h" |
41 | #include "AliITSneuralTracker.h" |
42 | |
43 | ClassImp(AliITSneuralTracker) |
44 | |
45 | |
46 | // ==================================================================================================== |
47 | // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> AliITSneuralTracker <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< |
48 | // ==================================================================================================== |
49 | |
50 | AliITSneuralTracker::AliITSneuralTracker() |
51 | { |
52 | // |
53 | // Default constructor (not for use) |
54 | // |
55 | Int_t i = 0; |
56 | fCurvCut = 0; |
57 | for (i = 0; i < 6; i++) { |
58 | fPoints[i] = 0; |
59 | if (i < 5) fNeurons[i] = 0; |
60 | } |
61 | fTracks = 0; |
62 | } |
63 | //--------------------------------------------------------------------------------------------------------- |
64 | AliITSneuralTracker::AliITSneuralTracker (Int_t nsecs, Int_t ncurv, Double_t *curv, Double_t theta) |
65 | { |
66 | // Argumented constructor |
67 | // ---------------------- |
68 | // gets the number of azymuthal sectors, |
69 | // the number of curvature cut steps and their cut values, |
70 | // and then the polar angle cut. |
71 | // Be careful not to put 'ncurv' greater than the dimension of the |
72 | // 'curv' array (--> segmentation fault) |
73 | |
74 | fSectorNum = nsecs; |
75 | fSectorWidth = 2.0 * TMath::Pi()/(Double_t)nsecs; |
76 | |
77 | fCurvNum = ncurv; |
78 | fCurvCut = new Double_t[ncurv]; |
79 | Int_t i; |
80 | for (i = 0; i < ncurv; i++) fCurvCut[i] = curv[i]; |
81 | |
82 | fThetaCut = theta; |
83 | fThetaNum = (Int_t)(TMath::Pi() / theta) + 1; |
84 | if (fThetaNum < 1) fThetaNum = 1; |
85 | |
86 | cout << "\nClass AliITSneuralTracker invoked with " << fSectorNum << " phi secs."; |
87 | cout << " and " << fThetaNum << " theta secs.\n" << endl; |
88 | |
89 | Int_t k, j; |
90 | for (k = 0; k < 6; k++) { |
91 | fPoints[k] = new TObjArray**[fSectorNum]; |
92 | for (i = 0; i < fSectorNum; i++) { |
93 | fPoints[k][i] = new TObjArray*[fThetaNum]; |
94 | for (j = 0; j < fThetaNum; j++) fPoints[k][i][j] = new TObjArray(0); |
95 | } |
96 | if (k < 6) fNeurons[k] = new TObjArray(0); |
97 | } |
98 | |
99 | fTracks = new TObjArray(1); |
100 | } |
101 | //--------------------------------------------------------------------------------------------------------- |
102 | AliITSneuralTracker::~AliITSneuralTracker() |
103 | { |
104 | // Destructor |
105 | // ---------- |
106 | // clears the TObjArrays and all heap objects |
107 | |
108 | Int_t lay, sec, k; |
109 | delete [] fCurvCut; |
110 | |
111 | for (lay = 0; lay < 6; lay++) { |
112 | for (sec = 0; sec < fSectorNum; sec++) { |
113 | for (k = 0; k < fThetaNum; k++) { |
114 | fPoints[lay][sec][k]->Delete(); |
115 | } |
116 | } |
117 | delete [] fPoints[lay]; |
118 | if (lay < 5) { |
119 | fNeurons[lay]->Delete(); |
120 | delete fNeurons[lay]; |
121 | } |
122 | } |
123 | |
124 | fTracks->Delete(); |
125 | } |
126 | //--------------------------------------------------------------------------------------------------------- |
127 | Int_t AliITSneuralTracker::ReadFile(const Text_t *fname, Int_t evnum) |
128 | { |
129 | // File reader |
130 | // ----------- |
131 | // Reads a ROOT file in order to retrieve recpoints. |
132 | // the 'evnum' argument can have two meanings: |
133 | // if it is >= 0 it represents the event number to retrieve the correct gAlice |
134 | // if it is < 0, instead, this means that the open file contains |
135 | // a 'TreeP' tree, with all points which remained unused |
136 | // after the Kalman tracking procedure (produced via the ITSafterKalman.C macro) |
137 | // the return value is the number of collected points |
138 | // (then, if errors occur, the method returns 0) |
139 | |
140 | TFile *file = new TFile(fname, "read"); |
141 | TTree *tree = 0; |
142 | Int_t i, nentries, total, sector, mesh; |
143 | |
144 | if (evnum < 0) { |
145 | tree = (TTree*)file->Get("TreeP0"); |
146 | if (!tree) { |
147 | Error("", "Specifying evnum < 0, I expect to find a 'TreeP' tree in the file, but there isn't any"); |
148 | return 0; |
149 | } |
150 | AliITSglobalRecPoint *p = 0; |
151 | tree->SetBranchAddress("Points", &p); |
152 | nentries = (Int_t)tree->GetEntries(); |
153 | total = 0; |
154 | for (i = 0; i < nentries; i++) { |
155 | tree->GetEntry(i); |
156 | AliITSglobalRecPoint *q = new AliITSglobalRecPoint(*p); |
157 | sector = (Int_t)(q->fPhi / fSectorWidth); |
158 | mesh = (Int_t)(q->fTheta / fThetaCut); |
159 | q->fUsed = 0; |
160 | fPoints[q->fLayer][sector][mesh]->AddLast(q); |
161 | total++; |
162 | } |
163 | return total; |
164 | } |
165 | |
166 | if (gAlice) { delete gAlice; gAlice = 0; } |
167 | gAlice = (AliRun*)file->Get("gAlice"); |
168 | if (!gAlice) { |
169 | Error("", "gAlice is NULL, maybe wrong filename..."); |
170 | return 0; |
171 | } |
172 | |
173 | Int_t nparticles = (Int_t)gAlice->GetEvent(evnum); |
174 | if (!nparticles) { |
175 | Error("", "No particles found"); |
176 | return 0; |
177 | } |
178 | |
179 | AliITS *its = (AliITS*)gAlice->GetModule("ITS"); |
180 | if (!its) { |
181 | Error("", "No ITS found"); |
182 | return 0; |
183 | } |
184 | |
185 | AliITSgeom *geom = (AliITSgeom*)its->GetITSgeom(); |
186 | if (!geom) { |
187 | Error("", "AliITSgeom is NULL"); |
188 | return 0; |
189 | } |
190 | |
191 | tree = gAlice->TreeR(); |
192 | if (!tree) { |
193 | Error("", "TreeR() returned NULL"); |
194 | return 0; |
195 | } |
196 | |
197 | nentries = (Int_t)tree->GetEntries(); |
198 | if (!nentries) { |
199 | Error("", "TreeR is empty"); |
200 | return 0; |
201 | } |
202 | |
203 | // Reading objects initialization |
204 | Int_t k, entry, lay, lad, det; |
205 | Double_t l[3], g[3], ll[3][3], gg[3][3]; |
206 | |
207 | TObjArray *recsArray = its->RecPoints(); |
208 | AliITSRecPoint *recp; |
209 | AliITSglobalRecPoint *pnt = 0; |
210 | |
211 | total = 0; |
212 | for(entry = 0; entry < nentries; entry++) { |
213 | if (entry > geom->GetLastSSD()) continue; |
214 | its->ResetRecPoints(); |
215 | tree->GetEvent(entry); |
216 | TObjArrayIter recs(recsArray); |
217 | AliITSgeomMatrix *gm = geom->GetGeomMatrix(entry); |
218 | geom->GetModuleId(entry, lay, lad, det); |
219 | lay--; |
220 | while ((recp = (AliITSRecPoint*)recs.Next())) { |
221 | // conversion to global coordinate system |
222 | for (i = 0; i < 3; i++) { |
223 | l[i] = g[i] = 0.0; |
224 | for (k = 0; k < 3; k++) { |
225 | ll[i][k] = 0.0; |
226 | gg[i][k] = 0.0; |
227 | } |
228 | } |
229 | // local to global conversions of coords |
230 | l[0] = recp->fX; |
231 | l[1] = 0.0; |
232 | l[2] = recp->fZ; |
233 | gm->LtoGPosition(l, g); |
234 | // local to global conversions of sigmas |
235 | ll[0][0] = recp->fSigmaX2; |
236 | ll[2][2] = recp->fSigmaZ2; |
237 | gm->LtoGPositionError(ll, gg); |
238 | // instantiation of a new global rec-point |
239 | pnt = new AliITSglobalRecPoint(g[0], g[1], g[2], gg[0][0], gg[1][1], gg[2][2], lay); |
240 | // copy of owner track labels |
241 | for (i = 0; i < 3; i++) pnt->fLabel[i] = recp->fTracks[i]; |
242 | sector = (Int_t)(pnt->fPhi / fSectorWidth); |
243 | mesh = (Int_t)(pnt->fTheta / fThetaCut); |
244 | pnt->fUsed = 0; |
245 | fPoints[lay][sector][mesh]->AddLast(pnt); |
246 | total++; |
247 | } |
248 | } |
249 | |
250 | return total; |
251 | } |
252 | //--------------------------------------------------------------------------------------------------------- |
253 | void AliITSneuralTracker::Go(const char* filesave, Bool_t flag) |
254 | { |
255 | // Global working method |
256 | // --------------------- |
257 | // It's the method which calls all other ones, |
258 | // in order to analyze the whole ITS points ensemble |
259 | // It's the only method to use, besides the parameter setters, |
260 | // so, all other methods are defined as private |
261 | // (see header file) |
262 | // the flag meaning is explained in the 'Initialize' method |
263 | |
264 | Int_t i; |
265 | for (i = 0; i < fSectorNum; i++) DoSector(i, flag); |
266 | cout << endl << "Saving results in file " << filesave << "..." << flush; |
267 | TFile *f = new TFile(filesave, "recreate"); |
268 | fTracks->Write(); |
269 | cout << "done" << endl; |
270 | f->Close(); |
271 | } |
272 | |
273 | //========================================================================================================= |
274 | // * * * P R I V A T E S E C T I O N * * * |
275 | //========================================================================================================= |
276 | |
277 | Int_t AliITSneuralTracker::Initialize(Int_t secnum, Int_t curvnum, Bool_t flag) |
278 | { |
279 | // Network initializer |
280 | // ------------------- |
281 | // Fills the neuron arrays, following the cut criteria for the selected step |
282 | // (secnum = sector to analyze, curvnum = curvature cut step to use) |
283 | // It also sets the initial random activation. |
284 | // In order to avoid a large number of 'new' operations, all existing neurons |
285 | // are reset and initialized with the new values, and are created new unit only if |
286 | // it turns out to be necessary |
287 | // the 'flag' argument is used to decide if the lower edge in the intevral |
288 | // of the accepted curvatures is given by zero (kFALSE) or by the preceding used cut (kTRUE) |
289 | // (if this is the first step, anyway, the minimum is always zero) |
290 | |
291 | Int_t l, m, k, neurons = 0; |
292 | |
293 | for (l = 0; l < 5; l++) fNeurons[l]->Delete(); |
294 | |
295 | AliITSneuron *unit = 0; |
296 | Double_t abscurv, max = fCurvCut[curvnum], min = 0.0; |
297 | if (flag && curvnum > 0) min = fCurvCut[curvnum - 1]; |
298 | AliITSglobalRecPoint *inner = 0, *outer = 0; |
299 | for (l = 0; l < 5; l++) { |
300 | for (m = 0; m < fThetaNum; m++) { |
301 | TObjArrayIter inners(fPoints[l][secnum][m]); |
302 | while ((inner = (AliITSglobalRecPoint*)inners.Next())) { |
303 | if (inner->fUsed > 0) continue; // points can't be recycled |
304 | for (k = m-1; k <= m+1; k++) { |
305 | if (k < 0 || k >= fThetaNum) continue; // to avoid seg faults |
306 | TObjArrayIter outers(fPoints[l+1][secnum][k]); |
307 | while ((outer = (AliITSglobalRecPoint*)outers.Next())) { |
308 | if (outer->fUsed > 0) continue; // points can't be recycled |
309 | if (inner->DTheta(outer) > fThetaCut) continue; |
310 | unit = new AliITSneuron; |
311 | unit->fInner = inner; |
312 | unit->fOuter = outer; |
313 | CalcParams(unit); |
314 | abscurv = TMath::Abs(unit->fCurv); |
315 | if (unit->fDiff > fDiff || abscurv < min || abscurv > max) { |
316 | delete unit; |
317 | continue; |
318 | } |
319 | unit->fActivation = gRandom->Rndm() * (fMax - fMin) + fMin; |
320 | fNeurons[l]->AddLast(unit); |
321 | neurons++; |
322 | } // end loop on candidates for outer point for neurons |
323 | } |
324 | } // end loop on candidates inner points for neuron |
325 | } |
326 | } // end loop on layers |
327 | return neurons; |
328 | } |
329 | //--------------------------------------------------------------------------------------------------------- |
330 | Bool_t AliITSneuralTracker::Update() |
331 | { |
332 | // Updating cycle |
333 | // -------------- |
334 | // Performs an updating cycle, by summing all the |
335 | // gain and cost contribution for each neuron and |
336 | // updating its activation. |
337 | // An asynchronous method is followed, with the order |
338 | // according that the neurons have been created |
339 | // Time wasting is avoided by dividing the neurons |
340 | // into different arrays, depending on the layer of their inner point. |
341 | // The return value answers the question: "The network has stabilized?" |
342 | |
343 | Int_t j, l; |
344 | Double_t actOld, actNew, argNew, totDiff = 0.0; |
345 | Double_t sumInh, sumExc, units = 0.0; |
346 | AliITSneuron *me = 0, *it = 0; |
347 | for (l = 0; l < 5; l++) { |
348 | TObjArrayIter meIter(fNeurons[l]); |
349 | while ((me = (AliITSneuron*)meIter.Next())) { |
350 | units++; |
351 | TObjArrayIter inhIter(fNeurons[l]); |
352 | |
353 | // inhibition (search only for neurons starting in the same layer) |
354 | sumInh = 0.0; |
355 | while((it = (AliITSneuron*)inhIter.Next())) { |
356 | if (it->fOuter == me->fOuter || it->fInner == me->fInner) |
357 | sumInh += it->fActivation; |
358 | } |
359 | |
360 | // excitation (search only for neurons |
361 | // which start in the previous or next layers) |
362 | sumExc = 0.0; |
363 | for (j = l - 1; j <= l + 1; j += 2) { |
364 | if (j < 0 || j >= 5) continue; |
365 | TObjArrayIter itIter(fNeurons[j]); |
366 | while ((it = (AliITSneuron*)itIter.Next())) { |
367 | if (it->fInner == me->fOuter || it->fOuter == me->fInner && it->fCurv * me->fCurv > 0.0) { |
368 | sumExc += Weight(me, it) * it->fActivation; |
369 | } |
370 | } |
371 | } // end search for excitories |
372 | |
373 | actOld = me->fActivation; |
374 | argNew = fRatio * sumExc - sumInh; |
375 | actNew = 1.0 / (1.0 + TMath::Exp(-argNew / fTemp)); |
376 | me->fActivation = actNew; |
377 | // calculation the relative activation variation |
378 | // (for stabilization check) |
379 | totDiff += TMath::Abs((actNew - actOld) / actOld); |
380 | } // end loop over 'me' (updated neuron) |
381 | } // end loop over layers |
382 | totDiff /= units; |
383 | |
384 | return (totDiff < fVar); |
385 | } |
386 | //--------------------------------------------------------------------------------------------------------- |
387 | Int_t AliITSneuralTracker::Save() |
388 | { |
389 | // Tracki saving method |
390 | // -------------------- |
391 | // Stores the recognized tracks into the TObjArray 'fTracks' |
392 | // but doesn't save it into a file until the whole work has ended |
393 | |
394 | Int_t l; |
395 | Double_t test = 0.5; |
396 | |
397 | // Before all, for each neuron is found the best active sequences |
398 | // among the incoming and outgoing other units |
399 | AliITSneuron *main, *fwd, *back, *start; |
400 | for (l = 0; l < 5; l++) { |
401 | TObjArrayIter mainIter(fNeurons[l]); |
402 | while((main = (AliITSneuron*)mainIter.Next())) { |
403 | if (l < 4) { |
404 | TObjArrayIter fwdIter(fNeurons[l+1]); |
405 | test = 0.5; |
406 | while ((fwd = (AliITSneuron*)fwdIter.Next())) { |
407 | if (main->fOuter == fwd->fInner && fwd->fActivation > test) { |
408 | test = fwd->fActivation; |
409 | main->fNext = fwd; |
410 | } |
411 | }; |
412 | } |
413 | if (l > 0) { |
414 | TObjArrayIter backIter(fNeurons[l-1]); |
415 | test = 0.5; |
416 | while ((back = (AliITSneuron*)backIter.Next())) { |
417 | if (main->fInner == back->fOuter && back->fActivation > test) { |
418 | test = back->fActivation; |
419 | main->fPrev = back; |
420 | } |
421 | }; |
422 | } |
423 | } |
424 | } |
425 | |
426 | // Then, are resolved the mismatches in the chains found above |
427 | // (when unit->next->prev != unit or unit->prev->next != unit) |
428 | for (l = 0; l < 5; l++) { |
429 | TObjArrayIter mainIter(fNeurons[l]); |
430 | while((main = (AliITSneuron*)mainIter.Next())) { |
431 | fwd = main->fNext; |
432 | back = main->fPrev; |
433 | if (fwd != 0 && fwd->fPrev != main) main->fNext = 0; |
434 | if (back != 0 && back->fNext != main) main->fPrev = 0; |
435 | } |
436 | } |
437 | |
438 | Int_t pos, neg, incr, decr, ntracks; |
439 | AliITSneuralTrack *trk = 0; |
440 | |
441 | TObjArrayIter startIter(fNeurons[0]); |
442 | for (;;) { |
443 | start = (AliITSneuron*)startIter.Next(); |
444 | if (!start) break; |
445 | Int_t pts = 0; |
446 | // with the chain above defined, a track can be followed |
447 | // like a linked list structure |
448 | for (main = start; main; main = main->fNext) pts++; |
449 | // the count will be > 5 only if the track contains 6 points |
450 | // (what we want) |
451 | if (pts < 5) continue; |
452 | // track storage |
453 | trk = new AliITSneuralTrack; |
454 | trk->CopyPoint(start->fInner); |
455 | for (main = start; main; main = main->fNext) { |
456 | //main->fInner->fUsed = kTRUE; |
457 | //main->fOuter->fUsed = kTRUE; |
458 | trk->CopyPoint(main->fOuter); |
459 | } |
460 | trk->Kinks(pos, neg, incr, decr); |
461 | if (pos != 0 && neg != 0) { |
462 | cout << "pos, neg kinks = " << pos << " " << neg << endl; |
463 | continue; |
464 | } |
465 | if (incr != 0 && decr != 0) { |
466 | cout << "increment, decrements in delta phi = " << incr << " " << decr << endl; |
467 | continue; |
468 | } |
469 | for (main = start; main; main = main->fNext) { |
470 | main->fInner->fUsed++; |
471 | main->fOuter->fUsed++; |
472 | } |
473 | fTracks->AddLast(trk); |
474 | } |
475 | ntracks = (Int_t)fTracks->GetEntriesFast(); |
476 | return ntracks; |
477 | } |
478 | //--------------------------------------------------------------------------------------------------------- |
479 | void AliITSneuralTracker::DoSector(Int_t sect, Bool_t flag) |
480 | { |
481 | // Sector recognition |
482 | // ------------------ |
483 | // This is a private method which works on a single sector |
484 | // just for an easier readability of the code. |
485 | // The 'flag' argument is needed to be given to the 'Initialize' method (see it) |
486 | |
487 | Int_t curvnum, nTracks = 0, nUnits; |
488 | Bool_t isStable = kFALSE; |
489 | for (curvnum = 0; curvnum < fCurvNum; curvnum++) { |
490 | cout << "\rSector " << sect << ":" << ((curvnum+1)*100)/fCurvNum << "%" << flush; |
491 | nUnits = Initialize(sect, curvnum, flag); |
492 | if (!nUnits) continue; |
493 | do { |
494 | isStable = Update(); |
495 | } while (!isStable); |
496 | nTracks = Save(); |
497 | } |
498 | cout << "\rTracks stored after sector " << sect << ": " << nTracks << endl; |
499 | } |
500 | //--------------------------------------------------------------------------------------------------------- |
501 | Int_t AliITSneuralTracker::CountExits(AliITSneuron *n) |
502 | { |
503 | // Counter for neurons which enter in 'n' tail |
504 | |
505 | Int_t count = 0, l = n->fOuter->fLayer; |
506 | if (l == 5) return 0; |
507 | |
508 | AliITSneuron *test = 0; |
509 | TObjArrayIter iter(fNeurons[l]); |
510 | while ((test = (AliITSneuron*)iter.Next())) { |
511 | if (test->fInner == n->fOuter) count++; |
512 | } |
513 | return count; |
514 | } |
515 | //--------------------------------------------------------------------------------------------------------- |
516 | Int_t AliITSneuralTracker::CountEnters(AliITSneuron *n) |
517 | { |
518 | // Counter for neurons which exit from 'n' head |
519 | |
520 | Int_t count = 0, l = n->fInner->fLayer; |
521 | if (l == 0) return 0; |
522 | |
523 | AliITSneuron *test = 0; |
524 | TObjArrayIter iter(fNeurons[l-1]); |
525 | while ((test = (AliITSneuron*)iter.Next())) { |
526 | if (test->fOuter == n->fInner) count++; |
527 | } |
528 | return count; |
529 | } |
530 | //--------------------------------------------------------------------------------------------------------- |
531 | Bool_t AliITSneuralTracker::ResetNeuron(AliITSneuron *n) |
532 | { |
533 | // A method which resets the neuron's pointers to zero |
534 | |
535 | if (!n) return kFALSE; |
536 | n->fNext = 0; |
537 | n->fPrev = 0; |
538 | n->fInner = 0; |
539 | n->fOuter = 0; |
540 | // the other datamembers don't need to be reinitialized |
541 | return kTRUE; |
542 | } |
543 | //--------------------------------------------------------------------------------------------------------- |
544 | Bool_t AliITSneuralTracker::SetEdge(AliITSneuron *n, AliITSglobalRecPoint *p, const char what) |
545 | { |
546 | // Sets the indicated edge point to the neuron |
547 | // (if the arguments are suitable ones...) |
548 | |
549 | if (!n || !p || (what != 'i' && what != 'I' && what != 'o' && what != 'O')) return kFALSE; |
550 | if (what == 'i' || what == 'I') |
551 | n->fInner = p; |
552 | else |
553 | n->fOuter = p; |
554 | return kTRUE; |
555 | } |
556 | //--------------------------------------------------------------------------------------------------------- |
557 | Bool_t AliITSneuralTracker::CalcParams(AliITSneuron *n) |
558 | { |
559 | // This method evaluates the helix parameters from the edged of a neuron |
560 | // (curvature, phase parameter, tangent of lambda angle) |
561 | // if the p[arameters assume too strange (and unphysical) values, the |
562 | // method return kFALSE, else it return kTRUE; |
563 | |
564 | if (!n) return kFALSE; |
565 | |
566 | Double_t sign(1.0), det, rc = 0.0; |
567 | Double_t tanl1, tanl2, l1, l2; |
568 | // changing the reference frame into the one centered in the vertex |
569 | Double_t x1 = n->fInner->fGX - fVPos[0]; |
570 | Double_t y1 = n->fInner->fGY - fVPos[1]; |
571 | Double_t z1 = n->fInner->fGZ - fVPos[2]; |
572 | Double_t r1 = x1 * x1 + y1 * y1; |
573 | Double_t x2 = n->fOuter->fGX - fVPos[0]; |
574 | Double_t y2 = n->fOuter->fGY - fVPos[1]; |
575 | Double_t z2 = n->fOuter->fGZ - fVPos[2]; |
576 | Double_t r2 = x2 * x2 + y2 * y2; |
577 | // initialization of the XY-plane data (curvature and centre) |
578 | // (will remain these values if we encounter errors) |
579 | n->fCurv = n->fCX = n->fCY = n->fTanL = n->fPhase = 0.0; |
580 | n->fLength = (n->fOuter->fGX - n->fInner->fGX) * (n->fOuter->fGX - n->fInner->fGX); |
581 | n->fLength += (n->fOuter->fGY - n->fInner->fGY) * (n->fOuter->fGY - n->fInner->fGY); |
582 | n->fLength += (n->fOuter->fGZ - n->fInner->fGZ) * (n->fOuter->fGZ - n->fInner->fGZ); |
583 | n->fLength = TMath::Sqrt(n->fLength); |
584 | // this determinant which is zero when the points are in a line |
585 | det = 2.0 * (x1 * y2 - x2 * y1); |
586 | if (det == 0.0) |
587 | return kFALSE; // it is difficult having perfectly aligned points --- it's more likely to be an error |
588 | else { |
589 | n->fCX = (r1*y2 - r2*y1) / det; |
590 | n->fCY = (r2*x1 - r1*x2) / det; |
591 | rc = TMath::Sqrt(n->fCX * n->fCX + n->fCY * n->fCY); |
592 | r1 = TMath::Sqrt(r1); |
593 | r2 = TMath::Sqrt(r2); |
594 | sign = (n->fOuter->fPhi >= n->fInner->fPhi) ? 1.0 : -1.0; |
595 | if (rc > 0.0) { |
596 | n->fCurv = sign / rc; |
597 | n->fPhase = TMath::ATan2(-n->fCX, -n->fCY); |
598 | if (n->fPhase < 0.0) n->fPhase += 2.0 * TMath::Pi(); // angles in 0 --> 2pi |
599 | if (r1 <= 2.0 * rc && r2 <= 2.0 * rc) { |
600 | l1 = 2.0 * rc * TMath::ASin(r1 / (2.0 * rc)); |
601 | l2 = 2.0 * rc * TMath::ASin(r2 / (2.0 * rc)); |
602 | tanl1 = z1 / l1; |
603 | tanl2 = z2 / l2; |
604 | n->fDiff = TMath::Abs(tanl1 - tanl2); |
605 | n->fTanL = 0.5 * (tanl1 + tanl2); |
606 | } |
607 | else { |
608 | n->fTanL = 0.0; |
609 | n->fDiff = 100.0; |
610 | return kFALSE; // it' even more difficult that the arguments above aren't acceptable for arcsine |
611 | } |
612 | } |
613 | else |
614 | return kFALSE; |
615 | } |
616 | return kTRUE; |
617 | } |
618 | //--------------------------------------------------------------------------------------------------------- |
619 | Double_t AliITSneuralTracker::Angle(AliITSneuron *n, AliITSneuron *m) |
620 | { |
621 | // calculates the angle between two segments |
622 | // but doesn't check if they have a common point. |
623 | // The calculation is made by the ratio of their scalar product |
624 | // to the product of their magnitudes |
625 | |
626 | Double_t x = (m->fOuter->fGX - m->fInner->fGX) * (n->fOuter->fGX - n->fInner->fGX); |
627 | Double_t y = (m->fOuter->fGY - m->fInner->fGY) * (n->fOuter->fGY - n->fInner->fGY); |
628 | Double_t z = (m->fOuter->fGZ - m->fInner->fGZ) * (n->fOuter->fGZ - n->fInner->fGZ); |
629 | |
630 | Double_t cosine = (x + y + z) / (n->fLength * m->fLength); |
631 | |
632 | if (cosine > 1.0 || cosine < -1.0) { |
633 | Warning("AliITSneuronV1::Angle", "Strange value of cosine"); |
634 | return 0.0; |
635 | } |
636 | |
637 | return TMath::ACos(cosine); |
638 | } |
639 | //--------------------------------------------------------------------------------------------------------- |
640 | Double_t AliITSneuralTracker::Weight(AliITSneuron *n, AliITSneuron *m) |
641 | { |
642 | // calculation of the excitoy weight |
643 | |
644 | // Double_t v1 = TMath::Abs(fTanL - n->fTanL); |
645 | // Double_t v2 = TMath::Abs(TMath::Abs(fCurv) - TMath::Abs(n->fCurv)); |
646 | // Double_t v3 = TMath::Abs(fPhase - n->fPhase); |
647 | // Double_t q = (v1 + v2 + v3) / 3.0; |
648 | // Double_t a = 1.0; // 0.02 / (fDiff + n->fDiff); |
649 | Double_t b = 1.0 - TMath::Sin(Angle(n, m)); |
650 | return TMath::Power(b, fExp); |
651 | // return TMath::Exp(q / par1); |
652 | |
653 | } |
654 | //--------------------------------------------------------------------------------------------------------- |
655 | // ==================================================================================================== |
656 | // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> AliITSneuron <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< |
657 | // ==================================================================================================== |
e73ba7c8 |
658 | ClassImp(AliITSneuron) |
1f26c323 |
659 | AliITSneuron::AliITSneuron() |
660 | { |
661 | // Default constructor |
662 | // ------------------- |
663 | // sets all parameters to default (but unuseful) |
664 | // (in particular all pointers are set to zero) |
665 | // There is no alternative constructor, because |
666 | // this class should not be used by itself, but |
667 | // it is only necessary to allow the neural tracker to work |
668 | // In fact, all methods which operate on neuron's datamembers |
669 | // are stored in the AliITSneuralTracker class anyway |
670 | |
671 | fActivation = 0.0; |
672 | |
673 | fCurv = 0.0; |
674 | fCX = 0.0; |
675 | fCY = 0.0; |
676 | fTanL = 0.0; |
677 | fPhase = 0.0; |
678 | fDiff = 0.0; |
679 | fLength = 0.0; |
680 | |
681 | fNext = fPrev = 0; |
682 | fInner = fOuter = 0; |
683 | } |
684 | //--------------------------------------------------------------------------------------------------------- |
685 | AliITSneuron::~AliITSneuron() |
686 | { |
687 | // Destructor |
688 | // ---------- |
689 | // does nothing, because there is no need to affect parameters, |
690 | // while the AliITSglobalRecPoint pointers can't be deleted, |
691 | // and the AliITSneuron pointers belong to a TObjArray which will |
692 | // be cleared when necessary... |
693 | } |
694 | |