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