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