]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSNeuralTracker.cxx
Neural ITS standalone tracking revisited
[u/mrichter/AliRoot.git] / ITS / AliITSNeuralTracker.cxx
CommitLineData
b9d722bc 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. *
12 * It is provided "as is" without express or implied warranty. *
13 * *
14 * AN ITS STAND-ALONE "NEURAL" TRACK FINDER *
15 * ---------------------------------------- *
16 * This class implements the Denby-Peterson algorithm for track finding *
17 * in the ITS stand-alone, by means of a neural network simulation. *
18 * Many parameters have to be set for the neural network to operate *
19 * correctly and with a good efficiency. *
20 * The neural tracker must be feeded with a TTree filled with objects *
21 * of the class "AliITSNeuralPoint", into a single branch called *
22 * "Points". *
23 **************************************************************************/
24
25//#define NEURAL_LINEAR
26#include <fstream>
27#include <Riostream.h>
28#include <stdlib.h>
29
30#include <TROOT.h>
31#include <TFile.h>
32#include <TTree.h>
33#include <TMath.h>
34#include <TLine.h>
35#include <TMarker.h>
36#include <TRandom.h>
37#include <TString.h>
38#include <TCanvas.h>
39#include <TVector3.h>
40#include <TParticle.h>
41#include <TObjArray.h>
42#include <TList.h>
43
44#include "AliITSNeuralPoint.h"
45#include "AliITSNeuralTracker.h"
46
47using namespace std;
48ClassImp(AliITSNeuralTracker)
49
50//--------------------------------------------------------------------------------------------
51
52AliITSNeuralTracker::AliITSNeuralTracker()
53{
54 // CONSTRUCTOR
55 //
56 // Initializes some data-members:
57 // - all pointers to NULL
58 // - theta cut to 180 deg. (= no theta cut)
59 // - initial choice of only 1 azimuthal sector
60 // - initial values for the neural network parameters.
61 //
62 // With these settings the neural tracker can't work
63 // because it has not any curvature cut set.
64
65 fCurvNum = 0;
66 fCurvCut = 0;
67
68 fSectorNum = 1;
69 fSectorWidth = TMath::Pi() * 2.0;
70 fPolarInterval = 45.0;
71
72 fStructureOK = kFALSE;
73
74 fVX = fVY = fVZ = 0.0;
75
76 Int_t ilayer, itheta;
77 for (ilayer = 0; ilayer < 6; ilayer++) {
78 for (itheta = 0; itheta < 180; itheta++) fPoints[ilayer][itheta] = 0;
79 fThetaCut2DMin[ilayer] = 0.0;
80 fThetaCut2DMax[ilayer] = TMath::Pi();
81 fThetaCut3DMin[ilayer] = 0.0;
82 fThetaCut3DMax[ilayer] = TMath::Pi();
83 fHelixMatchCutMin[ilayer] = 1.0;
84 fHelixMatchCutMax[ilayer] = 1.0;
85 }
86
87 fEdge1 = 0.3;
88 fEdge2= 0.7;
89
90 fTemperature = 1.0;
91 fStabThreshold = 0.001;
92 fGain2CostRatio = 1.0;
93 fAlignExponent = 1.0;
94 fActMinimum = 0.5;
95
96 fNeurons = 0;
97
98 fChains = new TTree("TreeC", "Sextines of points");
99 fChains->Branch("l0", &fPoint[0], "l0/I");
100 fChains->Branch("l1", &fPoint[1], "l1/I");
101 fChains->Branch("l2", &fPoint[2], "l2/I");
102 fChains->Branch("l3", &fPoint[3], "l3/I");
103 fChains->Branch("l4", &fPoint[4], "l4/I");
104 fChains->Branch("l5", &fPoint[5], "l5/I");
105}
106//
107//--------------------------------------------------------------------------------------------
108//
109AliITSNeuralTracker::~AliITSNeuralTracker()
110{
111 // DESTRUCTOR
112 //
113 // It Destroys all the dynamic arrays and
114 // clears the TCollections and the points tree
115
116 cout << "Starting destructor..." << endl;
117
118 delete [] fCurvCut;
119
120 Int_t ilayer, itheta;
121 if (fStructureOK) {
122 cout << "Deleting points..." << endl;
123 for (ilayer = 0; ilayer < 6; ilayer++) {
124 for (itheta = 0; itheta < 180; itheta++) {
125 fPoints[ilayer][itheta]->SetOwner();
126 delete fPoints[ilayer][itheta];
127 }
128 }
129 cout << "Deleting neurons..." << endl;
130 fNeurons->SetOwner();
131 delete fNeurons;
132 }
133
134 cout << "AliITSNeuralTracker destructed completely!" << endl;
135}
136//
137//--------------------------------------------------------------------------------------------
138//
139void AliITSNeuralTracker::Display(TCanvas*& canv)
140{
141 Double_t x1, y1, x2, y2;
142 canv->Clear();
143 TObjArrayIter iter(fNeurons);
144 for (;;) {
145 AliITSneuron *unit = (AliITSneuron*)iter.Next();
146 if (!unit) break;
147 if (unit->fActivation < fActMinimum) continue;
148 x1 = unit->fInner->X();
149 x2 = unit->fOuter->X();
150 y1 = unit->fInner->Y();
151 y2 = unit->fOuter->Y();
152 TLine *line = new TLine(x1, y1, x2, y2);
153 canv->cd();
154 line->Draw();
155 }
156 canv->Update();
157}
158//
159//--------------------------------------------------------------------------------------------
160//
161void AliITSNeuralTracker::SetThetaCuts2D(Double_t *min, Double_t *max)
162{
163 Int_t i;
164 Double_t temp;
165 for (i = 0; i < 5; i++) {
166 if (min[i] > max[i]) {
167 temp = min[i];
168 min[i] = max[i];
169 max[i] = temp;
170 }
171 fThetaCut2DMin[i] = min[i];
172 fThetaCut2DMax[i] = max[i];
173 }
174 for (i = 0; i < 5; i++) {
175 cout << "Theta 2D cut for layer " << i << " = " << fThetaCut2DMin[i] << " to " << fThetaCut2DMax[i] << endl;
176 }
177}
178//
179//--------------------------------------------------------------------------------------------
180//
181void AliITSNeuralTracker::SetThetaCuts3D(Double_t *min, Double_t *max)
182{
183 Int_t i;
184 Double_t temp;
185 for (i = 0; i < 5; i++) {
186 if (min[i] > max[i]) {
187 temp = min[i];
188 min[i] = max[i];
189 max[i] = temp;
190 }
191 fThetaCut3DMin[i] = min[i];
192 fThetaCut3DMax[i] = max[i];
193 }
194 for (i = 0; i < 5; i++) {
195 cout << "Theta 3D cut for layer " << i << " = " << fThetaCut3DMin[i] << " to " << fThetaCut3DMax[i] << endl;
196 }
197}
198//
199//--------------------------------------------------------------------------------------------
200//
201void AliITSNeuralTracker::SetHelixMatchCuts(Double_t *min, Double_t *max)
202{
203 Int_t i;
204 Double_t temp;
205 for (i = 0; i < 5; i++) {
206 if (min[i] > max[i]) {
207 temp = min[i];
208 min[i] = max[i];
209 max[i] = temp;
210 }
211 fHelixMatchCutMin[i] = min[i];
212 fHelixMatchCutMax[i] = max[i];
213 }
214 for (i = 0; i < 5; i++) {
215 cout << "Helix-match cut for layer " << i << " = " << fHelixMatchCutMin[i] << " to " << fHelixMatchCutMax[i] << endl;
216 }
217}
218//
219//--------------------------------------------------------------------------------------------
220//
221void AliITSNeuralTracker::SetCurvatureCuts(Int_t ncuts, Double_t *curv)
222{
223 // CURVATURE CUTS SETTER
224 //
225 // Requires an array of double values and its dimension
226 // After sorting it in increasing order, the array of curvature cuts
227 // is dinamically allocated, and filled with the sorted cuts array.
228 // A message is shown which lists all the curvature cuts.
229
230 Int_t i, *ind = new Int_t[ncuts];
231 TMath::Sort(ncuts, curv, ind, kFALSE);
232 fCurvCut = new Double_t[ncuts];
233 cout << "\n" << ncuts << " curvature cuts defined" << endl << "-----" << endl;
234 for (i = 0; i < ncuts; i++) {
235 fCurvCut[i] = curv[ind[i]];
236 cout << "Cut #" << i + 1 << " : " << fCurvCut[i] << endl;
237 }
238 cout << "-----" << endl;
239 fCurvNum = ncuts;
240}
241//
242//--------------------------------------------------------------------------------------------
243//
244void AliITSNeuralTracker::CreateArrayStructure(Int_t nsectors)
245{
246 // ARRAY CREATOR
247 //
248 // Organizes the array structure to store all points in.
249 //
250 // The array is organized into a "multi-level" TCollection:
251 // - 6 fPoints[] TObjArray containing a TObjArray for each layer
252 // - each TObject contains a TObjArray for each sector.
253
254 // sets the number of sectors and their width.
255 fSectorNum = nsectors;
256 fSectorWidth = TMath::Pi() * 2.0 / (Double_t)fSectorNum;
257
258 // creates the TCollection structure
259 Int_t ilayer, isector, itheta;
260 TObjArray *sector = 0;
261 for (ilayer = 0; ilayer < 6; ilayer++) {
262 for (itheta = 0; itheta < 180; itheta++) {
263 if (fPoints[ilayer][itheta]) {
264 fPoints[ilayer][itheta]->SetOwner();
265 delete fPoints[ilayer][itheta];
266 }
267 fPoints[ilayer][itheta] = new TObjArray(nsectors);
268 for (isector = 0; isector < nsectors; isector++) {
269 sector = new TObjArray;
270 sector->SetOwner();
271 fPoints[ilayer][itheta]->AddAt(sector, isector);
272 }
273 }
274 }
275
276 // Sets a checking flag to TRUE.
277 // This flag is checked before filling up the arrays with the points.
278 fStructureOK = kTRUE;
279}
280//
281//--------------------------------------------------------------------------------------------
282//
283Int_t AliITSNeuralTracker::ArrangePoints(TTree* pts_tree)
284{
285 // POINTS STORAGE INTO ARRAY
286 //
287 // Reads points from the tree and creates AliITSNode objects for each one,
288 // storing them into the array structure defined above.
289 // Returns the number of points collected (if successful) or 0 (otherwise)
290
291 // check: if the points tree is NULL or empty, there is nothing to do...
292 if ( !pts_tree || (pts_tree && !(Int_t)pts_tree->GetEntries()) ) {
293 Error("ArrangePoints", "Points tree is NULL or empty: no points to arrange");
294 return 0;
295 }
296
297 if (!fStructureOK) {
298 Error("ArrangePoints", "Structure NOT defined. Call CreateArrayStructure() first");
299 return 0;
300 }
301
302 Int_t isector, itheta, ientry, ilayer, nentries, pos;
303 TObjArray *sector = 0;
304 AliITSNode *created = 0;
305 AliITSNeuralPoint *cursor = 0;
306
307 pts_tree->SetBranchAddress("pos", &pos);
308 pts_tree->SetBranchAddress("Points", &cursor);
309 nentries = (Int_t)pts_tree->GetEntries();
310
311 for (ientry = 0; ientry < nentries; ientry++) {
312 pts_tree->GetEntry(ientry);
313 // creates the object
314 created = new AliITSNode(cursor, kTRUE);
315 created->SetUser(-1);
316 created->fPosInTree = pos;
317 // finds the sector in phi
318 isector = created->GetSector(fSectorWidth);
319 itheta = created->GetThetaCell();
320 ilayer = created->GetLayer();
321 if (ilayer < 0 || ilayer > 5) {
322 Error("ArrangePoints", "Layer value %d not allowed. Aborting...", ilayer);
323 return 0;
324 }
325 // selects the right TObjArray to store object into
326 sector = (TObjArray*)fPoints[ilayer][itheta]->At(isector);
327 sector->AddLast(created);
328 }
329
330 // returns the number of points processed
331 return ientry;
332}
333//
334//--------------------------------------------------------------------------------------------
335//
336void AliITSNeuralTracker::PrintPoints()
337{
338 // creates the TCollection structure
339 TObjArray *sector = 0;
340 Int_t ilayer, isector, itheta;
341 fstream file("test.log", ios::out);
342 for (ilayer = 0; ilayer < 6; ilayer++) {
343 for (isector = 0; isector < fSectorNum; isector++) {
344 for (itheta = 0; itheta < 180; itheta++) {
345 sector = (TObjArray*)fPoints[ilayer][itheta]->At(isector);
346 file << ilayer << " " << isector << " " << itheta;
347 file << " " << sector->GetSize() << " points" << endl;
348 }
349 }
350 }
351 file.close();
352}
353//
354//--------------------------------------------------------------------------------------------
355//
356Bool_t AliITSNeuralTracker::PassCurvCut
357(AliITSNode *p1, AliITSNode *p2, Int_t curv_index, Double_t vx, Double_t vy, Double_t vz)
358{
359 // CURVATURE CUT EVALUATOR
360 //
361 // Checks the passsed point pair w.r. to the current curvature cut
362 // Returns the result of the check.
363
364 if (curv_index < 0 || curv_index >= fCurvNum) {
365 Error("PassCurvCut", "Curv index %d out of range", curv_index);
366 return kFALSE;
367 }
368
369 // Find the reference layer
370 Int_t lay1 = p1->GetLayer();
371 Int_t lay2 = p2->GetLayer();
372 Int_t ref_layer = (lay1 < lay2) ? lay1 : lay2;
373
374 Double_t x1 = p1->X() - vx;
375 Double_t x2 = p2->X() - vx;
376 Double_t y1 = p1->Y() - vy;
377 Double_t y2 = p2->Y() - vy;
378 Double_t z1 = p1->Z() - vz;
379 Double_t z2 = p2->Z() - vz;
380 Double_t r1 = sqrt(x1*x1 + y1*y1);
381 Double_t r2 = sqrt(x2*x2 + y2*y2);
382
383 // calculation of curvature
384 Double_t dx = p1->X() - p2->X(), dy = p1->Y() - p2->Y();
385 Double_t num = 2 * (x1*y2 - x2*y1);
386 Double_t den = r1*r2*sqrt(dx*dx + dy*dy);
387 Double_t curv = 0.;
388 /* FOR OLD VERSION
389 if (den != 0.) {
390 curv = fabs(num / den);
391 if (curv > fCurvCut[curv_index]) return kFALSE;
392 return kTRUE;
393 }
394 else
395 return kFALSE;
396 */
397 // NEW VERSION
398 if (den != 0.) {
399 curv = fabs(num / den);
400 if (curv > fCurvCut[curv_index]) return kFALSE;
401 }
402 else
403 return kFALSE;
404 // calculation of helix matching
405 Double_t arc1 = 2.0 * r1 * curv;
406 Double_t arc2 = 2.0 * r2 * curv;
407 Double_t hel_match = 0.0;
408 if (arc1 > -1.0 && arc1 < 1.0) arc1 = asin(arc1);
409 else arc1 = ((arc1 > 0.0) ? 0.5 : 1.5) * TMath::Pi();
410 if (arc2 > -1.0 && arc2 < 1.0) arc2 = asin(arc2);
411 else arc2 = ((arc2 > 0.0) ? 0.5 : 1.5) * TMath::Pi();
412 arc1 /= 2.0 * curv;
413 arc2 /= 2.0 * curv;
414 if (arc1 == 0.0 || arc2 == 0.0) return kFALSE;
415 hel_match = fabs(z1 / arc1 - z2 / arc2);
416 return (hel_match >= fHelixMatchCutMin[ref_layer] && hel_match <= fHelixMatchCutMax[ref_layer]);
417 // END NEW VERSION
418}
419//
420//--------------------------------------------------------------------------------------------
421//
422Int_t AliITSNeuralTracker::PassAllCuts
423(AliITSNode *p1, AliITSNode *p2, Int_t curv_index, Double_t vx, Double_t vy, Double_t vz)
424{
425 // GLOBAL CUT EVALUATOR
426 //
427 // Checks all cuts for the passed point pair.
428 // Return values:
429 //
430 // 0 - All cuts passed
431 // 1 - theta 2D cut not passed
432 // 2 - theta 3D cut not passed
433 // 3 - curvature calculated but cut not passed
434 // 4 - curvature not calculated (division by zero)
435 // 5 - helix cut not passed
436 // 6 - curvature inxed out of range
437
438 if (curv_index < 0 || curv_index >= fCurvNum) return 6;
439
440 // Find the reference layer
441 Int_t lay1 = p1->GetLayer();
442 Int_t lay2 = p2->GetLayer();
443 Int_t ref_layer = (lay1 < lay2) ? lay1 : lay2;
444
445 // Swap points in order that r1 < r2
446 AliITSNode *temp = 0;
447 if (p2->GetLayer() < p1->GetLayer()) {
448 temp = p2;
449 p2 = p1;
450 p1 = temp;
451 }
452
453 // shift XY coords to the reference to the vertex position,
454 // for easier calculus of quantities.
455 Double_t x1 = p1->X() - vx;
456 Double_t x2 = p2->X() - vx;
457 Double_t y1 = p1->Y() - vy;
458 Double_t y2 = p2->Y() - vy;
459 Double_t z1 = p1->Z() - vz;
460 Double_t z2 = p2->Z() - vz;
461 Double_t r1 = sqrt(x1*x1 + y1*y1);
462 Double_t r2 = sqrt(x2*x2 + y2*y2);
463
464 // Check for theta cut
465 Double_t dtheta, dtheta3;
466 TVector3 v01(z1, r1, 0.0);
467 TVector3 v12(z2 - z1, r2 - r1, 0.0);
468 dtheta = v01.Angle(v12) * 180.0 / TMath::Pi();
469 TVector3 V01(x1, y1, z1);
470 TVector3 V12(x2 - x1, y2 - y1, z2 - z1);
471 dtheta3 = V01.Angle(V12) * 180.0 / TMath::Pi();
472 if (dtheta < fThetaCut2DMin[ref_layer] || dtheta > fThetaCut2DMax[ref_layer]) return 1;
473 if (dtheta3 < fThetaCut3DMin[ref_layer] || dtheta3 > fThetaCut3DMax[ref_layer]) return 2;
474
475 // calculation of curvature
476 Double_t dx = x1 - x2, dy = y1 - y2;
477 Double_t num = 2 * (x1*y2 - x2*y1);
478 Double_t den = r1*r2*sqrt(dx*dx + dy*dy);
479 Double_t curv = 0.;
480 if (den != 0.) {
481 curv = fabs(num / den);
482 if (curv > fCurvCut[curv_index]) return 3;
483 }
484 else
485 return 4;
486
487 // calculation of helix matching
488 Double_t arc1 = 2.0 * r1 * curv;
489 Double_t arc2 = 2.0 * r2 * curv;
490 Double_t hel_match = 0.0;
491 if (arc1 > -1.0 && arc1 < 1.0) arc1 = asin(arc1);
492 else arc1 = ((arc1 > 0.0) ? 0.5 : 1.5) * TMath::Pi();
493 if (arc2 > -1.0 && arc2 < 1.0) arc2 = asin(arc2);
494 else arc2 = ((arc2 > 0.0) ? 0.5 : 1.5) * TMath::Pi();
495 arc1 /= 2.0 * curv;
496 arc2 /= 2.0 * curv;
497 if (arc1 == 0.0 || arc2 == 0.0) return kFALSE;
498 hel_match = fabs(z1 / arc1 - z2 / arc2);
499 if (hel_match < fHelixMatchCutMin[ref_layer] || hel_match > fHelixMatchCutMax[ref_layer]) return 5;
500
501 return 0;
502}
503//
504//--------------------------------------------------------------------------------------------
505//
506void AliITSNeuralTracker::StoreAbsoluteMatches()
507{
508 // Stores in the 'fMatches' array of each node all the points in the
509 // adjacent layers which allow to create neurons accordin to the
510 // helix and theta cut, and also to the largest curvature cut
511
512 Int_t ilayer, isector, itheta1, itheta2, check;
513 TObjArray *list1 = 0, *list2 = 0;
514 AliITSNode *node1 = 0, *node2 = 0;
515 Double_t theta_min, theta_max;
516 Int_t imin, imax;
517
518 for (isector = 0; isector < fSectorNum; isector++) {
519 // sector is chosen once for both lists
520 for (ilayer = 0; ilayer < 5; ilayer++) {
521 for (itheta1 = 0; itheta1 < 180; itheta1++) {
522 list1 = (TObjArray*)fPoints[ilayer][itheta1]->At(isector);
523 TObjArrayIter iter1(list1);
524 while ( (node1 = (AliITSNode*)iter1.Next()) ) {
525 if (node1->GetUser() >= 0) continue;
526 node1->fMatches->Clear();
527 theta_min = node1->ThetaDeg() - fPolarInterval;
528 theta_max = node1->ThetaDeg() + fPolarInterval;
529 imin = (Int_t)theta_min;
530 imax = (Int_t)theta_max;
531 if (imin < 0) imin = 0;
532 if (imax > 179) imax = 179;
533 for (itheta2 = imin; itheta2 <= imax; itheta2++) {
534 list2 = (TObjArray*)fPoints[ilayer + 1][itheta2]->At(isector);
535 TObjArrayIter iter2(list2);
536 while ( (node2 = (AliITSNode*)iter2.Next()) ) {
537 check = PassAllCuts(node1, node2, fCurvNum - 1, fVX, fVY, fVZ);
538 switch (check) {
539 case 0:
540 node1->fMatches->AddLast(node2);
541 break;
542 case 1:
543 //Info("StoreAbsoluteMatches", "Layer %d: THETA 2D cut not passed", ilayer);
544 break;
545 case 2:
546 //Info("StoreAbsoluteMatches", "Layer %d: THETA 3D cut not passed", ilayer);
547 break;
548 case 3:
549 //Info("StoreAbsoluteMatches", "Layer %d: CURVATURE cut not passed", ilayer);
550 break;
551 case 4:
552 //Info("StoreAbsoluteMatches", "Layer %d: Division by ZERO in curvature evaluation", ilayer);
553 break;
554 case 5:
555 //Info("StoreAbsoluteMatches", "Layer %d: HELIX-MATCH cut not passed", ilayer);
556 break;
557 case 6:
558 //Error("PassAllCuts", "Layer %d: Curv index out of range", ilayer);
559 break;
560 default:
561 Warning("StoreAbsoluteMatches", "Layer %d: %d: unrecognized return value", ilayer, check);
562 }
563 } // while (node2...)
564 } // for (itheta2...)
565 } // while (node1...)
566 } // for (itheta...)
567 } // for (ilayer...)
568 } // for (isector...)
569}
570//
571//--------------------------------------------------------------------------------------------
572//
573void AliITSNeuralTracker::PrintMatches(Bool_t stop)
574{
575 TFile *ft = new TFile("its_findables_v1.root");
576 TTree *tt = (TTree*)ft->Get("Tracks");
577 Int_t it, nP, nU, lab, nF = (Int_t)tt->GetEntries();
578 tt->SetBranchAddress("nhitsP", &nP);
579 tt->SetBranchAddress("nhitsU", &nU);
580 tt->SetBranchAddress("label", &lab);
581 TString strP("|"), strU("|");
582 for (it = 0; it < nF; it++) {
583 tt->GetEntry(it);
584 if (nP >= 5) strP.Append(Form("%d|", lab));
585 if (nU >= 5) strU.Append(Form("%d|", lab));
586 }
587
588 TObjArray *sector = 0;
589 Int_t ilayer, isector, itheta, id[3];
590 AliITSNode *node1 = 0, *node2 = 0;
591
592 for (ilayer = 0; ilayer < 6; ilayer++) {
593 for (isector = 0; isector < fSectorNum; isector++) {
594 for (itheta = 0; itheta < 180; itheta++) {
595 sector = (TObjArray*)fPoints[ilayer][itheta]->At(isector);
596 TObjArrayIter points(sector);
597 while ( (node1 = (AliITSNode*)points.Next()) ) {
598 for (it = 0; it < 3; it++) id[it] = node1->GetLabel(it);
599 nF = (Int_t)node1->fMatches->GetSize();
600 cout << "Node layer: " << node1->GetLayer() << ", labels: ";
601 cout << id[0] << " " << id[1] << " " << id[2] << " --> ";
602 if (!nF) {
603 cout << "NO MatchES!!!" << endl;
604 continue;
605 }
606 cout << nF << " Matches" << endl;
607 for (it = 0; it < 3; it++) {
608 if (strP.Contains(Form("|%d|", id[it])))
609 cout << "Belongs to findable (originary) track " << id[it] << endl;
610 if (strU.Contains(Form("|%d|", id[it])))
611 cout << "Belongs to findable (post-Kalman) track " << id[it] << endl;
612 }
613 TObjArrayIter Matches(node1->fMatches);
614 while ( (node2 = (AliITSNode*)Matches.Next()) ) {
615 cout << "Match with " << node2;
616 Int_t *sh = node1->SharedID(node2);
617 for (Int_t k = 0; k < 3; k++)
618 if (sh[k] > -1) cout << " " << sh[k];
619 cout << endl;
620 }
621 if (stop) {
622 cout << "Press a key" << endl;
623 cin.get();
624 }
625 }
626 }
627 }
628 }
629}
630//
631//--------------------------------------------------------------------------------------------
632//
633void AliITSNeuralTracker::ResetNodes(Int_t isector)
634{
635 Int_t ilayer, itheta;
636 TObjArray *sector = 0;
637 AliITSNode *node = 0;
638 for (ilayer = 0; ilayer < 6; ilayer++) {
639 for (itheta = 0; itheta < 180; itheta++) {
640 sector = (TObjArray*)fPoints[ilayer][itheta]->At(isector);
641 TObjArrayIter iter(sector);
642 for (;;) {
643 node = (AliITSNode*)iter.Next();
644 if (!node) break;
645 node->fInnerOf->Clear();
646 node->fOuterOf->Clear();
647 delete node->fInnerOf;
648 delete node->fOuterOf;
649 node->fInnerOf = new TObjArray;
650 node->fOuterOf = new TObjArray;
651 }
652 }
653 }
654}
655//
656//--------------------------------------------------------------------------------------------
657//
658void AliITSNeuralTracker::NeuralTracking(const char* filesave, TCanvas*& display)
659// This is the public method that operates the tracking.
660// It works sector by sector, and at the end saves the found tracks.
661// Other methods are privare because they have no meaning id used alone,
662// and sometimes they could get segmentation faults due to uninitialized
663// datamembert they require to work on.
664// The argument is a file where the final results have to be stored.
665{
666 Bool_t isStable = kFALSE;
667 Int_t i, nUnits = 0, nLinks = 0, nTracks = 0, sectTracks = 0, totTracks = 0;
668
669 // tracking through sectors
670 cout << endl;
671 Int_t sector, curv;
672 for (sector = 0; sector < fSectorNum; sector++) {
673 cout << "\rSector " << sector << ": " << endl;
674 sectTracks = 0;
675 for(curv = 0; curv < fCurvNum; curv++) {
676 cout << "- curvature step " << curv + 1;
677 cout << " (" << fCurvCut[curv] << "): ";
678 // units creation
679 nUnits = CreateNeurons(sector, curv);
680 if (!nUnits) {
681 cout << "no neurons --> skipped" << endl;
682 continue;
683 }
684 cout << nUnits << " units, " << flush;
685 // units connection
686 nLinks = LinkNeurons();
687 if (!nLinks) {
688 cout << "no connections --> skipped" << endl;
689 continue;
690 }
691 cout << nLinks << " connections, " << flush;
692 // neural updating
693 for (i = 0;; i++) {
694 isStable = Update();
695 if (display) Display(display);
696 TObjArrayIter iter(fNeurons);
697 for (;;) {
698 AliITSneuron *n = (AliITSneuron*)iter.Next();
699 if (!n) break;
700 }
701 if (isStable) break;
702 }
703 cout << i << " cycles --> " << flush;
704 // tracks saving
705 CleanNetwork();
706 nTracks = Save(sector);
707 cout << nTracks << " tracks" << endl;
708 sectTracks += nTracks;
709 }
710 totTracks += sectTracks;
711 //cout << sectTracks << " tracks found (total = " << totTracks << ") " << flush;
712 }
713
714 cout << endl << totTracks << " tracks found!!!" << endl;
715 cout << endl << "Saving results in file " << filesave << "..." << flush;
716 TFile *f = new TFile(filesave, "recreate");
717 fChains->Write("TreeC");
718 f->Close();
719}
720//
721//--------------------------------------------------------------------------------------------
722//
723Int_t AliITSNeuralTracker::CreateNeurons(Int_t sector_idx, Int_t curv_idx)
724// Fills the neuron arrays, following the cut criteria for the selected step
725// (secnum = sector to analyze, curvnum = curvature cut step to use)
726// It also sets the initial random activation.
727// In order to avoid a large number of 'new' operations, all existing neurons
728// are reset and initialized with the new values, and are created new unit only if
729// it turns out to be necessary
730// the 'flag' argument is used to decide if the lower edge in the intevral
731// of the accepted curvatures is given by zero (kFALSE) or by the preceding used cut (kTRUE)
732// (if this is the first step, anyway, the minimum is always zero)
733{
734 ResetNodes(sector_idx);
735
736 if (fNeurons) delete fNeurons;
737 fNeurons = new TObjArray;
738
739 AliITSneuron *unit = 0;
740 Int_t itheta, neurons = 0;
741 TObjArray *lst_sector = 0;
742
743 // NEW VERSION
744 Double_t vx[6], vy[6], vz[6];
745 AliITSNode *p[6] = {0, 0, 0, 0, 0, 0};
746 for (itheta = 0; itheta < 180; itheta++) {
747 lst_sector = (TObjArray*)fPoints[0][itheta]->At(sector_idx);
748 TObjArrayIter lay0(lst_sector);
749 while ( (p[0] = (AliITSNode*)lay0.Next()) ) {
750 if (p[0]->GetUser() >= 0) continue;
751 vx[0] = fVX;
752 vy[0] = fVY;
753 vz[0] = fVZ;
754 TObjArrayIter lay1(p[0]->fMatches);
755 while ( (p[1] = (AliITSNode*)lay1.Next()) ) {
756 if (p[1]->GetUser() >= 0) continue;
757 if (!PassCurvCut(p[0], p[1], curv_idx, fVX, fVY, fVZ)) continue;
758 unit = new AliITSneuron;
759 unit->fInner = p[0];
760 unit->fOuter = p[1];
761 unit->fActivation = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
762 unit->fGain = new TObjArray;
763 fNeurons->AddLast(unit);
764 p[0]->fInnerOf->AddLast(unit);
765 p[1]->fOuterOf->AddLast(unit);
766 neurons++;
767 vx[1] = p[0]->X();
768 vy[1] = p[0]->Y();
769 vz[1] = p[0]->Z();
770 TObjArrayIter lay2(p[1]->fMatches);
771 while ( (p[2] = (AliITSNode*)lay2.Next()) ) {
772 if (p[2]->GetUser() >= 0) continue;
773 if (!PassCurvCut(p[1], p[2], curv_idx, vx[1], vy[1], vz[1])) continue;
774 unit = new AliITSneuron;
775 unit->fInner = p[1];
776 unit->fOuter = p[2];
777 unit->fActivation = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
778 unit->fGain = new TObjArray;
779 fNeurons->AddLast(unit);
780 p[1]->fInnerOf->AddLast(unit);
781 p[2]->fOuterOf->AddLast(unit);
782 neurons++;
783 vx[2] = p[1]->X();
784 vy[2] = p[1]->Y();
785 vz[2] = p[1]->Z();
786 TObjArrayIter lay3(p[2]->fMatches);
787 while ( (p[3] = (AliITSNode*)lay3.Next()) ) {
788 if (p[3]->GetUser() >= 0) continue;
789 if (!PassCurvCut(p[2], p[3], curv_idx, vx[2], vy[2], vz[2])) continue;
790 unit = new AliITSneuron;
791 unit->fInner = p[2];
792 unit->fOuter = p[3];
793 unit->fActivation = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
794 unit->fGain = new TObjArray;
795 fNeurons->AddLast(unit);
796 p[2]->fInnerOf->AddLast(unit);
797 p[3]->fOuterOf->AddLast(unit);
798 neurons++;
799 vx[3] = p[2]->X();
800 vy[3] = p[2]->Y();
801 vz[3] = p[2]->Z();
802 TObjArrayIter lay4(p[3]->fMatches);
803 while ( (p[4] = (AliITSNode*)lay4.Next()) ) {
804 if (p[4]->GetUser() >= 0) continue;
805 if (!PassCurvCut(p[3], p[4], curv_idx, vx[3], vy[3], vz[3])) continue;
806 unit = new AliITSneuron;
807 unit->fInner = p[3];
808 unit->fOuter = p[4];
809 unit->fActivation = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
810 unit->fGain = new TObjArray;
811 fNeurons->AddLast(unit);
812 p[3]->fInnerOf->AddLast(unit);
813 p[4]->fOuterOf->AddLast(unit);
814 neurons++;
815 vx[4] = p[3]->X();
816 vy[4] = p[3]->Y();
817 vz[4] = p[3]->Z();
818 TObjArrayIter lay5(p[4]->fMatches);
819 while ( (p[5] = (AliITSNode*)lay5.Next()) ) {
820 if (p[5]->GetUser() >= 0) continue;
821 if (!PassCurvCut(p[4], p[5], curv_idx, vx[4], vy[4], vz[4])) continue;
822 unit = new AliITSneuron;
823 unit->fInner = p[4];
824 unit->fOuter = p[5];
825 unit->fActivation = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
826 unit->fGain = new TObjArray;
827 fNeurons->AddLast(unit);
828 p[4]->fInnerOf->AddLast(unit);
829 p[5]->fOuterOf->AddLast(unit);
830 neurons++;
831 } // while (p[5])
832 } // while (p[4])
833 } // while (p[3])
834 } // while (p[2])
835 } // while (p[1])
836 } // while (p[0])
837 } // for (itheta...)
838 // END OF NEW VERSION
839
840 /* OLD VERSION
841 for (ilayer = 0; ilayer < 6; ilayer++) {
842 for (itheta = 0; itheta < 180; itheta++) {
843 lst_sector = (TObjArray*)fPoints[ilayer][itheta]->At(sector_idx);
844 TObjArrayIter inners(lst_sector);
845 while ( (inner = (AliITSNode*)inners.Next()) ) {
846 if (inner->GetUser() >= 0) continue;
847 TObjArrayIter outers(inner->fMatches);
848 while ( (outer = (AliITSNode*)outers.Next()) ) {
849 if (outer->GetUser() >= 0) continue;
850 if (!PassCurvCut(inner, outer, curv_idx, fVX, fVY, fVZ)) continue;
851 unit = new AliITSneuron;
852 unit->fInner = inner;
853 unit->fOuter = outer;
854 unit->fActivation = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2;
855 unit->fGain = new TObjArray;
856 fNeurons->AddLast(unit);
857 inner->fInnerOf->AddLast(unit);
858 outer->fOuterOf->AddLast(unit);
859 neurons++;
860 } // for (;;)
861 } // for (;;)
862 } // for (itheta...)
863 } // for (ilayer...)
864 */
865
866 fNeurons->SetOwner();
867 return neurons;
868}
869//
870//
871//
872Int_t AliITSNeuralTracker::LinkNeurons()
873// Creates the neural synapses among all neurons
874// which share a point.
875{
876 Int_t total = 0;
877 TObjArrayIter iter(fNeurons), *test_iter;
878 AliITSneuron *neuron = 0, *test = 0;
879 for (;;) {
880 neuron = (AliITSneuron*)iter.Next();
881 if (!neuron) break;
882 // gain contributors
883 test_iter = (TObjArrayIter*)neuron->fInner->fOuterOf->MakeIterator();
884 for (;;) {
885 test = (AliITSneuron*)test_iter->Next();
886 if (!test) break;
887 neuron->Add2Gain(test, fGain2CostRatio, fAlignExponent);
888 total++;
889 }
890 delete test_iter;
891 test_iter = (TObjArrayIter*)neuron->fOuter->fInnerOf->MakeIterator();
892 for (;;) {
893 test = (AliITSneuron*)test_iter->Next();
894 if (!test) break;
895 neuron->Add2Gain(test, fGain2CostRatio, fAlignExponent);
896 total++;
897 }
898 delete test_iter;
899 }
900 return total;
901}
902//
903//
904//
905Bool_t AliITSNeuralTracker::Update()
906// A complete updating cycle with the asynchronous method
907// when the neural network is stable, kTRUE is returned.
908{
909 Double_t actVar = 0.0, totDiff = 0.0;
910 TObjArrayIter iter(fNeurons);
911 AliITSneuron *unit;
912 for (;;) {
913 unit = (AliITSneuron*)iter.Next();
914 if (!unit) break;
915 actVar = Activate(unit);
916 // calculation the relative activation variation
917 totDiff += actVar;
918 }
919 totDiff /= fNeurons->GetSize();
920 return (totDiff < fStabThreshold);
921}
922//
923//
924//
925void AliITSNeuralTracker::CleanNetwork()
926// Removes all deactivated neurons, and resolves the competitions
927// to the advantage of the most active unit
928{
929 AliITSneuron *unit = 0;
930 TObjArrayIter neurons(fNeurons);
931 while ( (unit = (AliITSneuron*)neurons.Next()) ) {
932 if (unit->fActivation < fActMinimum) {
933 unit->fInner->fInnerOf->Remove(unit);
934 unit->fOuter->fOuterOf->Remove(unit);
935 delete fNeurons->Remove(unit);
936 }
937 }
938 return;
939 Bool_t removed;
940 Int_t nIn, nOut;
941 AliITSneuron *enemy = 0;
942 neurons.Reset();
943 while ( (unit = (AliITSneuron*)neurons.Next()) ) {
944 nIn = (Int_t)unit->fInner->fInnerOf->GetSize();
945 nOut = (Int_t)unit->fOuter->fOuterOf->GetSize();
946 if (nIn < 2 && nOut < 2) continue;
947 removed = kFALSE;
948 if (nIn > 1) {
949 TObjArrayIter competing(unit->fInner->fInnerOf);
950 while ( (enemy = (AliITSneuron*)competing.Next()) ) {
951 if (unit->fActivation > enemy->fActivation) {
952 enemy->fInner->fInnerOf->Remove(enemy);
953 enemy->fOuter->fOuterOf->Remove(enemy);
954 delete fNeurons->Remove(enemy);
955 }
956 else {
957 unit->fInner->fInnerOf->Remove(unit);
958 unit->fOuter->fOuterOf->Remove(unit);
959 delete fNeurons->Remove(unit);
960 removed = kTRUE;
961 break;
962 }
963 }
964 if (removed) continue;
965 }
966 if (nOut > 1) {
967 TObjArrayIter competing(unit->fOuter->fOuterOf);
968 while ( (enemy = (AliITSneuron*)competing.Next()) ) {
969 if (unit->fActivation > enemy->fActivation) {
970 enemy->fInner->fInnerOf->Remove(enemy);
971 enemy->fOuter->fOuterOf->Remove(enemy);
972 delete fNeurons->Remove(enemy);
973 }
974 else {
975 unit->fInner->fInnerOf->Remove(unit);
976 unit->fOuter->fOuterOf->Remove(unit);
977 delete fNeurons->Remove(unit);
978 removed = kTRUE;
979 break;
980 }
981 }
982 }
983 }
984}
985//
986//
987//
988Bool_t AliITSNeuralTracker::CheckOccupation()
989{
990 Int_t i;
991 for (i = 0; i < 6; i++) if (fPoint[i] < 0) return kFALSE;
992 return kTRUE;
993}
994//
995//
996//
997Int_t AliITSNeuralTracker::Save(Int_t sector_id)
998// Creates chains of active neurons, in order to
999// find the tracks obtained as the neural network output.
1000{
1001 // every chain is started from the neurons in the first 2 layers
1002 Int_t i, check, stored = 0;
1003 Double_t test_act = 0;
1004 AliITSneuron *unit = 0, *cursor = 0, *fwd = 0;
1005 AliITSNode *node = 0;
1006 TObjArrayIter iter(fNeurons), *fwd_iter;
1007 TObjArray *removedUnits = new TObjArray(0);
1008 TObjArray *removedPoints = new TObjArray(6);
1009 removedUnits->SetOwner(kFALSE);
1010 removedPoints->SetOwner(kFALSE);
1011 for (;;) {
1012 for (i = 0; i < 6; i++) fPoint[i] = -1;
1013 unit = (AliITSneuron*)iter.Next();
1014 if (!unit) break;
1015 if (unit->fInner->GetLayer() > 0) continue;
1016 fPoint[unit->fInner->GetLayer()] = unit->fInner->fPosInTree;
1017 fPoint[unit->fOuter->GetLayer()] = unit->fOuter->fPosInTree;
1018 node = unit->fOuter;
1019 removedUnits->AddLast(unit);
1020 removedPoints->AddAt(unit->fInner, unit->fInner->GetLayer());
1021 removedPoints->AddAt(unit->fOuter, unit->fOuter->GetLayer());
1022 while (node) {
1023 test_act = fActMinimum;
1024 fwd_iter = (TObjArrayIter*)node->fInnerOf->MakeIterator();
1025 fwd = 0;
1026 for (;;) {
1027 cursor = (AliITSneuron*)fwd_iter->Next();
1028 if (!cursor) break;
1029 if (cursor->fUsed) continue;
1030 if (cursor->fActivation >= test_act) {
1031 test_act = cursor->fActivation;
1032 fwd = cursor;
1033 }
1034 }
1035 if (!fwd) break;
1036 removedUnits->AddLast(fwd);
1037 node = fwd->fOuter;
1038 fPoint[node->GetLayer()] = node->fPosInTree;
1039 removedPoints->AddAt(node, node->GetLayer());
1040 }
1041 check = 0;
1042 for (i = 0; i < 6; i++) if (fPoint[i] >= 0) check++;
1043 if (check >= 6) {
1044 stored++;
1045 fChains->Fill();
1046 for (i = 0; i < 6; i++) {
1047 node = (AliITSNode*)removedPoints->At(i);
1048 node->SetUser(1);
1049 }
1050 fwd_iter = (TObjArrayIter*)removedUnits->MakeIterator();
1051 for (;;) {
1052 cursor = (AliITSneuron*)fwd_iter->Next();
1053 if(!cursor) break;
1054 cursor->fUsed = 1;
1055 }
1056 }
1057 }
1058
1059 return stored;
1060}
1061/*
1062Int_t AliITSNeuralTracker::Save(Int_t sector_idx)
1063// Creates chains of active neurons, in order to
1064// find the tracks obtained as the neural network output.
1065{
1066 // Reads the final map of neurons and removes all
1067 // units with too small activation
1068 cout << "saving..." << flush;
1069
1070 // every chain is started from the neurons in the first 2 layers
1071 // 00111111 00011111 00101111 00110111
1072 Int_t allow_mask[] = { 0x3F, 0x1F, 0x2F, 0x37,
1073 0x3B, 0x3D, 0x3E };
1074 // 00111011 00111101 00111110
1075
1076 Double_t test_fwd = 0., test_back = 0.;
1077 Int_t ilayer, itheta;
1078 AliITSNode *node = 0;
1079 AliITSneuron *fwd = 0, *back = 0;
1080 TList *list_sector = 0;
1081
1082 cout << "A -" << fActMinimum << "-" << flush;
1083 for (ilayer = 0; ilayer < 6; ilayer++) {
1084 for (itheta = 0; itheta < 180; itheta++) {
1085 list_sector = (TList*)fPoints[ilayer][itheta]->At(sector_idx);
1086 TListIter iter(list_sector);
1087 while ( (node = (AliITSNode*)iter.Next()) ) {
1088 TListIter fwd_iter(node->fInnerOf);
1089 TListIter back_iter(node->fOuterOf);
1090 test_fwd = test_back = fActMinimum;
1091 while ( (fwd = (AliITSneuron*)fwd_iter.Next()) ) {
1092 if (fwd->fActivation > test_fwd) {
1093 test_fwd = fwd->fActivation;
1094 node->fNext = fwd->fOuter;
1095 }
1096 }
1097 while ( (back = (AliITSneuron*)back_iter.Next()) ) {
1098 if (back->fActivation > test_back) {
1099 test_back = back->fActivation;
1100 node->fPrev = back->fInner;
1101 }
1102 }
1103 }
1104 }
1105 }
1106
1107 cout << "B" << flush;
1108 for (ilayer = 0; ilayer < 5; ilayer++) {
1109 for (itheta = 0; itheta < 180; itheta++) {
1110 list_sector = (TList*)fPoints[ilayer][itheta]->At(sector_idx);
1111 TListIter iter(list_sector);
1112 while ( (node = (AliITSNode*)iter.Next()) ) {
1113 if (node->fNext) {
1114 if (node->fNext->fPrev != node) node->fNext = 0;
1115 }
1116 }
1117 }
1118 }
1119
1120 cout << "C" << flush;
1121 Bool_t ok;
1122 Int_t stored = 0, l1, l2, i, mask, im;
1123 AliITSNode *temp = 0;
1124 TList *list[2];
1125 for (itheta = 0; itheta < 180; itheta++) {
1126 list[0] = (TList*)fPoints[0][itheta]->At(sector_idx);
1127 list[1] = (TList*)fPoints[1][itheta]->At(sector_idx);
1128 for (i = 0; i < 2; i++) {
1129 TListIter iter(list[i]);
1130 while ( (node = (AliITSNode*)iter.Next()) ) {
1131 //cout << endl << node << flush;
1132 AliITSneuralChain *chain = new AliITSneuralChain;
1133 for (temp = node; temp; temp = temp->fNext) {
1134 chain->Insert((AliITSNeuralPoint*)temp);
1135 }
1136 ok = kFALSE;
1137 mask = chain->OccupationMask();
1138 for (im = 0; im < 7; im++) ok = ok || (mask == allow_mask[im]);
1139 if (!ok) {
1140 delete chain;
1141 continue;
1142 }
1143 //cout << " lab " << flush;
1144 chain->AssignLabel();
1145 //cout << " add " << flush;
1146 fTracks->AddLast(chain);
1147 stored++;
1148 //cout << " " << stored << flush;
1149 }
1150 }
1151 }
1152
1153 cout << "D" << flush;
1154 return stored;
1155}
1156*/
1157//
1158//
1159//
1160Double_t AliITSNeuralTracker::Activate(AliITSneuron* &unit)
1161// Computes the new neural activation and
1162// returns the variation w.r.t the previous one
1163{
1164 if (!unit) return 0.0;
1165 Double_t sum_gain = 0.0, sum_cost = 0.0, input, actOld, actNew;
1166
1167 // sum gain contributions
1168 TObjArrayIter *iter = (TObjArrayIter*)unit->fGain->MakeIterator();
1169 AliITSlink *link;
1170 for(;;) {
1171 link = (AliITSlink*)iter->Next();
1172 if (!link) break;
1173 sum_gain += link->fLinked->fActivation * link->fWeight;
1174 }
1175
1176 // sum cost contributions
1177 TObjArrayIter *test_iter = (TObjArrayIter*)unit->fInner->fInnerOf->MakeIterator();
1178 AliITSneuron *linked = 0;
1179 for (;;) {
1180 linked = (AliITSneuron*)test_iter->Next();
1181 if (!linked) break;
1182 if (linked == unit) continue;
1183 sum_cost += linked->fActivation;
1184 }
1185 delete test_iter;
1186 test_iter = (TObjArrayIter*)unit->fOuter->fOuterOf->MakeIterator();
1187 for (;;) {
1188 linked = (AliITSneuron*)test_iter->Next();
1189 if (!linked) break;
1190 if (linked == unit) continue;
1191 sum_cost += linked->fActivation;
1192 }
1193
1194 //cout << "gain = " << sum_gain << ", cost = " << sum_cost << endl;
1195
1196 input = (sum_gain - sum_cost) / fTemperature;
1197 actOld = unit->fActivation;
1198#ifdef NEURAL_LINEAR
1199 if (input <= -2.0 * fTemperature)
1200 actNew = 0.0;
1201 else if (input >= 2.0 * fTemperature)
1202 actNew = 1.0;
1203 else
1204 actNew = input / (4.0 * fTemperature) + 0.5;
1205#else
1206 actNew = 1.0 / (1.0 + TMath::Exp(-input));
1207#endif
1208 unit->fActivation = actNew;
1209 return TMath::Abs(actNew - actOld);
1210}
1211//
1212//
1213//
1214void AliITSNeuralTracker::AliITSneuron::Add2Gain(AliITSneuron *n, Double_t mult_const, Double_t exponent)
1215{
1216 AliITSlink *link = new AliITSlink;
1217 link->fLinked = n;
1218 Double_t weight = Weight(n);
1219 link->fWeight = mult_const * TMath::Power(weight, exponent);
1220 fGain->AddLast(link);
1221}
1222//
1223//
1224//
1225Double_t AliITSNeuralTracker::AliITSneuron::Weight(AliITSneuron *n)
1226{
1227 TVector3 me(fOuter->X() - fInner->X(), fOuter->Y() - fInner->Y(), fOuter->Z() - fInner->Z());
1228 TVector3 it(n->fOuter->X() - n->fInner->X(), n->fOuter->Y() - n->fInner->Y(), n->fOuter->Z() - n->fInner->Z());
1229
1230 Double_t angle = me.Angle(it);
1231 Double_t weight = 1.0 - sin(angle);
1232 return weight;
1233}