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