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