]>
Commit | Line | Data |
---|---|---|
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 | 46 | ClassImp(AliITSNeuralTracker) |
47 | ||
48 | //-------------------------------------------------------------------------------------------- | |
49 | ||
cb729508 | 50 | AliITSNeuralTracker::AliITSNeuralTracker(): |
51 | fSectorNum(1), | |
52 | fSectorWidth(0), | |
53 | fPolarInterval(45.0), | |
54 | fCurvNum(0), | |
55 | fCurvCut(0), | |
56 | fStructureOK(kFALSE), | |
57 | fVX(0.0), | |
58 | fVY(0.0), | |
59 | fVZ(0.0), | |
60 | fActMinimum(0.5), | |
61 | fEdge1(0.3), | |
62 | fEdge2(0.7), | |
63 | fTemperature(1.0), | |
64 | fStabThreshold(0.001), | |
65 | fGain2CostRatio(1.0), | |
66 | fAlignExponent(1.0), | |
67 | fChains(0), | |
68 | fNeurons(0){ | |
b9d722bc | 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 | ||
b9d722bc | 80 | |
b9d722bc | 81 | fSectorWidth = TMath::Pi() * 2.0; |
b9d722bc | 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 | ||
b9d722bc | 94 | |
b9d722bc | 95 | |
b9d722bc | 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 | } | |
cb729508 | 105 | |
106 | AliITSNeuralTracker::AliITSNeuralTracker(const AliITSNeuralTracker &n):TObject(n), | |
107 | fSectorNum(n.fSectorNum), | |
108 | fSectorWidth(n.fSectorWidth), | |
109 | fPolarInterval(n.fPolarInterval), | |
110 | fCurvNum(n.fCurvNum), | |
111 | fCurvCut(n.fCurvCut), | |
112 | fStructureOK(n.fStructureOK), | |
113 | fVX(n.fVX), | |
114 | fVY(n.fVY), | |
115 | fVZ(n.fVZ), | |
116 | fActMinimum(n.fActMinimum), | |
117 | fEdge1(n.fEdge1), | |
118 | fEdge2(n.fEdge2), | |
119 | fTemperature(n.fTemperature), | |
120 | fStabThreshold(n.fStabThreshold), | |
121 | fGain2CostRatio(n.fGain2CostRatio), | |
122 | fAlignExponent(n.fAlignExponent), | |
123 | fChains(n.fChains), | |
124 | fNeurons(n.fNeurons){ | |
125 | //copy contructor | |
126 | } | |
127 | ||
128 | AliITSNeuralTracker& AliITSNeuralTracker::operator=(const AliITSNeuralTracker& t){ | |
129 | //assignment operator | |
130 | this->~AliITSNeuralTracker(); | |
131 | new(this) AliITSNeuralTracker(t); | |
132 | return *this; | |
133 | } | |
b9d722bc | 134 | // |
135 | //-------------------------------------------------------------------------------------------- | |
136 | // | |
137 | AliITSNeuralTracker::~AliITSNeuralTracker() | |
138 | { | |
139 | // DESTRUCTOR | |
140 | // | |
141 | // It Destroys all the dynamic arrays and | |
142 | // clears the TCollections and the points tree | |
143 | ||
b9d722bc | 144 | delete [] fCurvCut; |
145 | ||
146 | Int_t ilayer, itheta; | |
147 | if (fStructureOK) { | |
b9d722bc | 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 | } | |
b9d722bc | 154 | fNeurons->SetOwner(); |
155 | delete fNeurons; | |
156 | } | |
157 | ||
b9d722bc | 158 | } |
159 | // | |
160 | //-------------------------------------------------------------------------------------------- | |
161 | // | |
853a0f19 | 162 | void AliITSNeuralTracker::Display(TCanvas*& canv) const |
b9d722bc | 163 | { |
853a0f19 | 164 | // Displays the neural segments |
165 | ||
b9d722bc | 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; | |
ddcce466 | 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(); | |
b9d722bc | 177 | TLine *line = new TLine(x1, y1, x2, y2); |
178 | canv->cd(); | |
179 | line->Draw(); | |
180 | } | |
181 | canv->Update(); | |
182 | } | |
183 | // | |
184 | //-------------------------------------------------------------------------------------------- | |
185 | // | |
186 | void AliITSNeuralTracker::SetThetaCuts2D(Double_t *min, Double_t *max) | |
187 | { | |
853a0f19 | 188 | // Cut setter |
189 | ||
b9d722bc | 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 | // | |
208 | void AliITSNeuralTracker::SetThetaCuts3D(Double_t *min, Double_t *max) | |
209 | { | |
853a0f19 | 210 | // Cut setter |
211 | ||
b9d722bc | 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 | // | |
230 | void AliITSNeuralTracker::SetHelixMatchCuts(Double_t *min, Double_t *max) | |
231 | { | |
853a0f19 | 232 | // Cut setter |
233 | ||
b9d722bc | 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 | // | |
252 | void 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 | // | |
275 | void 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 | // | |
853a0f19 | 314 | Int_t AliITSNeuralTracker::ArrangePoints(TTree* ptstree) |
b9d722bc | 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... | |
853a0f19 | 323 | if ( !ptstree || (ptstree && !(Int_t)ptstree->GetEntries()) ) { |
b9d722bc | 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 | ||
853a0f19 | 338 | ptstree->SetBranchAddress("pos", &pos); |
339 | ptstree->SetBranchAddress("Points", &cursor); | |
340 | nentries = (Int_t)ptstree->GetEntries(); | |
b9d722bc | 341 | |
342 | for (ientry = 0; ientry < nentries; ientry++) { | |
853a0f19 | 343 | ptstree->GetEntry(ientry); |
b9d722bc | 344 | // creates the object |
345 | created = new AliITSNode(cursor, kTRUE); | |
346 | created->SetUser(-1); | |
ddcce466 | 347 | created->PosInTree() = pos; |
b9d722bc | 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 | // | |
367 | void 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 | // | |
387 | Bool_t AliITSNeuralTracker::PassCurvCut | |
853a0f19 | 388 | (AliITSNode *p1, AliITSNode *p2, Int_t curvindex, Double_t vx, Double_t vy, Double_t vz) |
b9d722bc | 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 | ||
853a0f19 | 395 | if (curvindex < 0 || curvindex >= fCurvNum) { |
396 | Error("PassCurvCut", "Curv index %d out of range", curvindex); | |
b9d722bc | 397 | return kFALSE; |
398 | } | |
399 | ||
400 | // Find the reference layer | |
401 | Int_t lay1 = p1->GetLayer(); | |
402 | Int_t lay2 = p2->GetLayer(); | |
853a0f19 | 403 | Int_t reflayer = (lay1 < lay2) ? lay1 : lay2; |
b9d722bc | 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); | |
853a0f19 | 422 | if (curv > fCurvCut[curvindex]) return kFALSE; |
b9d722bc | 423 | return kTRUE; |
424 | } | |
425 | else | |
426 | return kFALSE; | |
427 | */ | |
428 | // NEW VERSION | |
429 | if (den != 0.) { | |
b9f05b32 | 430 | curv = TMath::Abs(num / den); |
853a0f19 | 431 | if (curv > fCurvCut[curvindex]) return kFALSE; |
b9d722bc | 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; | |
853a0f19 | 438 | Double_t helmatch = 0.0; |
b9d722bc | 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; | |
b9f05b32 | 446 | helmatch = TMath::Abs(z1 / arc1 - z2 / arc2); |
853a0f19 | 447 | return (helmatch >= fHelixMatchCutMin[reflayer] && helmatch <= fHelixMatchCutMax[reflayer]); |
b9d722bc | 448 | // END NEW VERSION |
449 | } | |
450 | // | |
451 | //-------------------------------------------------------------------------------------------- | |
452 | // | |
453 | Int_t AliITSNeuralTracker::PassAllCuts | |
853a0f19 | 454 | (AliITSNode *p1, AliITSNode *p2, Int_t curvindex, Double_t vx, Double_t vy, Double_t vz) |
b9d722bc | 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 | ||
853a0f19 | 469 | if (curvindex < 0 || curvindex >= fCurvNum) return 6; |
b9d722bc | 470 | |
471 | // Find the reference layer | |
472 | Int_t lay1 = p1->GetLayer(); | |
473 | Int_t lay2 = p2->GetLayer(); | |
853a0f19 | 474 | Int_t reflayer = (lay1 < lay2) ? lay1 : lay2; |
b9d722bc | 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(); | |
853a0f19 | 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; | |
b9d722bc | 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.) { | |
b9f05b32 | 512 | curv = TMath::Abs(num / den); |
853a0f19 | 513 | if (curv > fCurvCut[curvindex]) return 3; |
b9d722bc | 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; | |
853a0f19 | 521 | Double_t helmatch = 0.0; |
b9d722bc | 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; | |
b9f05b32 | 529 | helmatch = TMath::Abs(z1 / arc1 - z2 / arc2); |
853a0f19 | 530 | if (helmatch < fHelixMatchCutMin[reflayer] || helmatch > fHelixMatchCutMax[reflayer]) return 5; |
b9d722bc | 531 | |
532 | return 0; | |
533 | } | |
534 | // | |
535 | //-------------------------------------------------------------------------------------------- | |
536 | // | |
537 | void 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; | |
853a0f19 | 546 | Double_t thetamin, thetamax; |
b9d722bc | 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; | |
ddcce466 | 557 | node1->Matches()->Clear(); |
853a0f19 | 558 | thetamin = node1->ThetaDeg() - fPolarInterval; |
559 | thetamax = node1->ThetaDeg() + fPolarInterval; | |
560 | imin = (Int_t)thetamin; | |
561 | imax = (Int_t)thetamax; | |
b9d722bc | 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: | |
ddcce466 | 571 | node1->Matches()->AddLast(node2); |
b9d722bc | 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 | // | |
604 | void AliITSNeuralTracker::PrintMatches(Bool_t stop) | |
605 | { | |
853a0f19 | 606 | // Prints the matches for each point |
607 | ||
b9d722bc | 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); | |
ddcce466 | 632 | nF = (Int_t)node1->Matches()->GetSize(); |
b9d722bc | 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 | } | |
ddcce466 | 646 | TObjArrayIter matches(node1->Matches()); |
853a0f19 | 647 | while ( (node2 = (AliITSNode*)matches.Next()) ) { |
b9d722bc | 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 | // | |
666 | void AliITSNeuralTracker::ResetNodes(Int_t isector) | |
667 | { | |
853a0f19 | 668 | // Clears some utility arrays for each node |
669 | ||
b9d722bc | 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; | |
ddcce466 | 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; | |
b9d722bc | 686 | } |
687 | } | |
688 | } | |
689 | } | |
690 | // | |
691 | //-------------------------------------------------------------------------------------------- | |
692 | // | |
693 | void AliITSNeuralTracker::NeuralTracking(const char* filesave, TCanvas*& display) | |
853a0f19 | 694 | { |
b9d722bc | 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. | |
853a0f19 | 701 | |
b9d722bc | 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 | // | |
853a0f19 | 759 | Int_t AliITSNeuralTracker::CreateNeurons(Int_t sectoridx, Int_t curvidx) |
760 | { | |
b9d722bc | 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) | |
853a0f19 | 770 | |
771 | ResetNodes(sectoridx); | |
b9d722bc | 772 | |
773 | if (fNeurons) delete fNeurons; | |
774 | fNeurons = new TObjArray; | |
775 | ||
776 | AliITSneuron *unit = 0; | |
777 | Int_t itheta, neurons = 0; | |
853a0f19 | 778 | TObjArray *lstsector = 0; |
b9d722bc | 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++) { | |
853a0f19 | 784 | lstsector = (TObjArray*)fPoints[0][itheta]->At(sectoridx); |
785 | TObjArrayIter lay0(lstsector); | |
b9d722bc | 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; | |
ddcce466 | 791 | TObjArrayIter lay1(p[0]->Matches()); |
b9d722bc | 792 | while ( (p[1] = (AliITSNode*)lay1.Next()) ) { |
793 | if (p[1]->GetUser() >= 0) continue; | |
853a0f19 | 794 | if (!PassCurvCut(p[0], p[1], curvidx, fVX, fVY, fVZ)) continue; |
b9d722bc | 795 | unit = new AliITSneuron; |
ddcce466 | 796 | unit->Inner() = p[0]; |
797 | unit->Outer() = p[1]; | |
798 | unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2; | |
799 | unit->Gain() = new TObjArray; | |
b9d722bc | 800 | fNeurons->AddLast(unit); |
ddcce466 | 801 | p[0]->InnerOf()->AddLast(unit); |
802 | p[1]->OuterOf()->AddLast(unit); | |
b9d722bc | 803 | neurons++; |
804 | vx[1] = p[0]->X(); | |
805 | vy[1] = p[0]->Y(); | |
806 | vz[1] = p[0]->Z(); | |
ddcce466 | 807 | TObjArrayIter lay2(p[1]->Matches()); |
b9d722bc | 808 | while ( (p[2] = (AliITSNode*)lay2.Next()) ) { |
809 | if (p[2]->GetUser() >= 0) continue; | |
853a0f19 | 810 | if (!PassCurvCut(p[1], p[2], curvidx, vx[1], vy[1], vz[1])) continue; |
b9d722bc | 811 | unit = new AliITSneuron; |
ddcce466 | 812 | unit->Inner() = p[1]; |
813 | unit->Outer() = p[2]; | |
814 | unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2; | |
815 | unit->Gain() = new TObjArray; | |
b9d722bc | 816 | fNeurons->AddLast(unit); |
ddcce466 | 817 | p[1]->InnerOf()->AddLast(unit); |
818 | p[2]->OuterOf()->AddLast(unit); | |
b9d722bc | 819 | neurons++; |
820 | vx[2] = p[1]->X(); | |
821 | vy[2] = p[1]->Y(); | |
822 | vz[2] = p[1]->Z(); | |
ddcce466 | 823 | TObjArrayIter lay3(p[2]->Matches()); |
b9d722bc | 824 | while ( (p[3] = (AliITSNode*)lay3.Next()) ) { |
825 | if (p[3]->GetUser() >= 0) continue; | |
853a0f19 | 826 | if (!PassCurvCut(p[2], p[3], curvidx, vx[2], vy[2], vz[2])) continue; |
b9d722bc | 827 | unit = new AliITSneuron; |
ddcce466 | 828 | unit->Inner() = p[2]; |
829 | unit->Outer() = p[3]; | |
830 | unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2; | |
831 | unit->Gain() = new TObjArray; | |
b9d722bc | 832 | fNeurons->AddLast(unit); |
ddcce466 | 833 | p[2]->InnerOf()->AddLast(unit); |
834 | p[3]->OuterOf()->AddLast(unit); | |
b9d722bc | 835 | neurons++; |
836 | vx[3] = p[2]->X(); | |
837 | vy[3] = p[2]->Y(); | |
838 | vz[3] = p[2]->Z(); | |
ddcce466 | 839 | TObjArrayIter lay4(p[3]->Matches()); |
b9d722bc | 840 | while ( (p[4] = (AliITSNode*)lay4.Next()) ) { |
841 | if (p[4]->GetUser() >= 0) continue; | |
853a0f19 | 842 | if (!PassCurvCut(p[3], p[4], curvidx, vx[3], vy[3], vz[3])) continue; |
b9d722bc | 843 | unit = new AliITSneuron; |
ddcce466 | 844 | unit->Inner() = p[3]; |
845 | unit->Outer() = p[4]; | |
846 | unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2; | |
847 | unit->Gain() = new TObjArray; | |
b9d722bc | 848 | fNeurons->AddLast(unit); |
ddcce466 | 849 | p[3]->InnerOf()->AddLast(unit); |
850 | p[4]->OuterOf()->AddLast(unit); | |
b9d722bc | 851 | neurons++; |
852 | vx[4] = p[3]->X(); | |
853 | vy[4] = p[3]->Y(); | |
854 | vz[4] = p[3]->Z(); | |
ddcce466 | 855 | TObjArrayIter lay5(p[4]->Matches()); |
b9d722bc | 856 | while ( (p[5] = (AliITSNode*)lay5.Next()) ) { |
857 | if (p[5]->GetUser() >= 0) continue; | |
853a0f19 | 858 | if (!PassCurvCut(p[4], p[5], curvidx, vx[4], vy[4], vz[4])) continue; |
b9d722bc | 859 | unit = new AliITSneuron; |
ddcce466 | 860 | unit->Inner() = p[4]; |
861 | unit->Outer() = p[5]; | |
862 | unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2; | |
863 | unit->Gain() = new TObjArray; | |
b9d722bc | 864 | fNeurons->AddLast(unit); |
ddcce466 | 865 | p[4]->InnerOf()->AddLast(unit); |
866 | p[5]->OuterOf()->AddLast(unit); | |
b9d722bc | 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++) { | |
853a0f19 | 880 | lstsector = (TObjArray*)fPoints[ilayer][itheta]->At(sectoridx); |
881 | TObjArrayIter inners(lstsector); | |
b9d722bc | 882 | while ( (inner = (AliITSNode*)inners.Next()) ) { |
883 | if (inner->GetUser() >= 0) continue; | |
ddcce466 | 884 | TObjArrayIter outers(inner->Matches()); |
b9d722bc | 885 | while ( (outer = (AliITSNode*)outers.Next()) ) { |
886 | if (outer->GetUser() >= 0) continue; | |
853a0f19 | 887 | if (!PassCurvCut(inner, outer, curvidx, fVX, fVY, fVZ)) continue; |
b9d722bc | 888 | unit = new AliITSneuron; |
ddcce466 | 889 | unit->Inner() = inner; |
890 | unit->Outer() = outer; | |
891 | unit->Activation() = gRandom->Rndm() * (fEdge1 - fEdge2) + fEdge2; | |
892 | unit->Gain() = new TObjArray; | |
b9d722bc | 893 | fNeurons->AddLast(unit); |
ddcce466 | 894 | inner->InnerOf()->AddLast(unit); |
895 | outer->OuterOf()->AddLast(unit); | |
b9d722bc | 896 | neurons++; |
897 | } // for (;;) | |
898 | } // for (;;) | |
899 | } // for (itheta...) | |
900 | } // for (ilayer...) | |
901 | */ | |
902 | ||
903 | fNeurons->SetOwner(); | |
904 | return neurons; | |
905 | } | |
906 | // | |
907 | // | |
908 | // | |
853a0f19 | 909 | Int_t AliITSNeuralTracker::LinkNeurons() const |
910 | { | |
b9d722bc | 911 | // Creates the neural synapses among all neurons |
912 | // which share a point. | |
853a0f19 | 913 | |
b9d722bc | 914 | Int_t total = 0; |
853a0f19 | 915 | TObjArrayIter iter(fNeurons), *testiter; |
b9d722bc | 916 | AliITSneuron *neuron = 0, *test = 0; |
917 | for (;;) { | |
918 | neuron = (AliITSneuron*)iter.Next(); | |
919 | if (!neuron) break; | |
920 | // gain contributors | |
ddcce466 | 921 | testiter = (TObjArrayIter*)neuron->Inner()->OuterOf()->MakeIterator(); |
b9d722bc | 922 | for (;;) { |
853a0f19 | 923 | test = (AliITSneuron*)testiter->Next(); |
b9d722bc | 924 | if (!test) break; |
925 | neuron->Add2Gain(test, fGain2CostRatio, fAlignExponent); | |
926 | total++; | |
927 | } | |
853a0f19 | 928 | delete testiter; |
ddcce466 | 929 | testiter = (TObjArrayIter*)neuron->Outer()->InnerOf()->MakeIterator(); |
b9d722bc | 930 | for (;;) { |
853a0f19 | 931 | test = (AliITSneuron*)testiter->Next(); |
b9d722bc | 932 | if (!test) break; |
933 | neuron->Add2Gain(test, fGain2CostRatio, fAlignExponent); | |
934 | total++; | |
935 | } | |
853a0f19 | 936 | delete testiter; |
b9d722bc | 937 | } |
938 | return total; | |
939 | } | |
940 | // | |
941 | // | |
942 | // | |
943 | Bool_t AliITSNeuralTracker::Update() | |
853a0f19 | 944 | { |
b9d722bc | 945 | // A complete updating cycle with the asynchronous method |
946 | // when the neural network is stable, kTRUE is returned. | |
853a0f19 | 947 | |
b9d722bc | 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 | // | |
964 | void AliITSNeuralTracker::CleanNetwork() | |
853a0f19 | 965 | { |
b9d722bc | 966 | // Removes all deactivated neurons, and resolves the competitions |
967 | // to the advantage of the most active unit | |
853a0f19 | 968 | |
b9d722bc | 969 | AliITSneuron *unit = 0; |
970 | TObjArrayIter neurons(fNeurons); | |
971 | while ( (unit = (AliITSneuron*)neurons.Next()) ) { | |
ddcce466 | 972 | if (unit->Activation() < fActMinimum) { |
973 | unit->Inner()->InnerOf()->Remove(unit); | |
974 | unit->Outer()->OuterOf()->Remove(unit); | |
b9d722bc | 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()) ) { | |
ddcce466 | 984 | nIn = (Int_t)unit->Inner()->InnerOf()->GetSize(); |
985 | nOut = (Int_t)unit->Outer()->OuterOf()->GetSize(); | |
b9d722bc | 986 | if (nIn < 2 && nOut < 2) continue; |
987 | removed = kFALSE; | |
988 | if (nIn > 1) { | |
ddcce466 | 989 | TObjArrayIter competing(unit->Inner()->InnerOf()); |
b9d722bc | 990 | while ( (enemy = (AliITSneuron*)competing.Next()) ) { |
ddcce466 | 991 | if (unit->Activation() > enemy->Activation()) { |
992 | enemy->Inner()->InnerOf()->Remove(enemy); | |
993 | enemy->Outer()->OuterOf()->Remove(enemy); | |
b9d722bc | 994 | delete fNeurons->Remove(enemy); |
995 | } | |
996 | else { | |
ddcce466 | 997 | unit->Inner()->InnerOf()->Remove(unit); |
998 | unit->Outer()->OuterOf()->Remove(unit); | |
b9d722bc | 999 | delete fNeurons->Remove(unit); |
1000 | removed = kTRUE; | |
1001 | break; | |
1002 | } | |
1003 | } | |
1004 | if (removed) continue; | |
1005 | } | |
1006 | if (nOut > 1) { | |
ddcce466 | 1007 | TObjArrayIter competing(unit->Outer()->OuterOf()); |
b9d722bc | 1008 | while ( (enemy = (AliITSneuron*)competing.Next()) ) { |
ddcce466 | 1009 | if (unit->Activation() > enemy->Activation()) { |
1010 | enemy->Inner()->InnerOf()->Remove(enemy); | |
1011 | enemy->Outer()->OuterOf()->Remove(enemy); | |
b9d722bc | 1012 | delete fNeurons->Remove(enemy); |
1013 | } | |
1014 | else { | |
ddcce466 | 1015 | unit->Inner()->InnerOf()->Remove(unit); |
1016 | unit->Outer()->OuterOf()->Remove(unit); | |
b9d722bc | 1017 | delete fNeurons->Remove(unit); |
1018 | removed = kTRUE; | |
1019 | break; | |
1020 | } | |
1021 | } | |
1022 | } | |
1023 | } | |
1024 | } | |
1025 | // | |
1026 | // | |
1027 | // | |
853a0f19 | 1028 | Bool_t AliITSNeuralTracker::CheckOccupation() const |
b9d722bc | 1029 | { |
853a0f19 | 1030 | // checks if a track recognized has a point for each layer |
b9d722bc | 1031 | Int_t i; |
1032 | for (i = 0; i < 6; i++) if (fPoint[i] < 0) return kFALSE; | |
1033 | return kTRUE; | |
1034 | } | |
1035 | // | |
1036 | // | |
1037 | // | |
853a0f19 | 1038 | Int_t AliITSNeuralTracker::Save(Int_t sectorid) |
1039 | { | |
b9d722bc | 1040 | // Creates chains of active neurons, in order to |
1041 | // find the tracks obtained as the neural network output. | |
853a0f19 | 1042 | |
b9d722bc | 1043 | // every chain is started from the neurons in the first 2 layers |
1044 | Int_t i, check, stored = 0; | |
853a0f19 | 1045 | Int_t a = sectorid; a = 0; |
1046 | Double_t testact = 0; | |
b9d722bc | 1047 | AliITSneuron *unit = 0, *cursor = 0, *fwd = 0; |
1048 | AliITSNode *node = 0; | |
853a0f19 | 1049 | TObjArrayIter iter(fNeurons), *fwditer; |
b9d722bc | 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; | |
ddcce466 | 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(); | |
b9d722bc | 1062 | removedUnits->AddLast(unit); |
ddcce466 | 1063 | removedPoints->AddAt(unit->Inner(), unit->Inner()->GetLayer()); |
1064 | removedPoints->AddAt(unit->Outer(), unit->Outer()->GetLayer()); | |
b9d722bc | 1065 | while (node) { |
853a0f19 | 1066 | testact = fActMinimum; |
ddcce466 | 1067 | fwditer = (TObjArrayIter*)node->InnerOf()->MakeIterator(); |
b9d722bc | 1068 | fwd = 0; |
1069 | for (;;) { | |
853a0f19 | 1070 | cursor = (AliITSneuron*)fwditer->Next(); |
b9d722bc | 1071 | if (!cursor) break; |
ddcce466 | 1072 | if (cursor->Used()) continue; |
1073 | if (cursor->Activation() >= testact) { | |
1074 | testact = cursor->Activation(); | |
b9d722bc | 1075 | fwd = cursor; |
1076 | } | |
1077 | } | |
1078 | if (!fwd) break; | |
1079 | removedUnits->AddLast(fwd); | |
ddcce466 | 1080 | node = fwd->Outer(); |
1081 | fPoint[node->GetLayer()] = node->PosInTree(); | |
b9d722bc | 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 | } | |
853a0f19 | 1093 | fwditer = (TObjArrayIter*)removedUnits->MakeIterator(); |
b9d722bc | 1094 | for (;;) { |
853a0f19 | 1095 | cursor = (AliITSneuron*)fwditer->Next(); |
b9d722bc | 1096 | if(!cursor) break; |
ddcce466 | 1097 | cursor->Used() = 1; |
b9d722bc | 1098 | } |
1099 | } | |
1100 | } | |
1101 | ||
1102 | return stored; | |
1103 | } | |
1104 | /* | |
853a0f19 | 1105 | Int_t AliITSNeuralTracker::Save(Int_t sectoridx) |
b9d722bc | 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 | |
853a0f19 | 1115 | Int_t allowmask[] = { 0x3F, 0x1F, 0x2F, 0x37, |
b9d722bc | 1116 | 0x3B, 0x3D, 0x3E }; |
1117 | // 00111011 00111101 00111110 | |
1118 | ||
853a0f19 | 1119 | Double_t test_fwd = 0., testback = 0.; |
b9d722bc | 1120 | Int_t ilayer, itheta; |
1121 | AliITSNode *node = 0; | |
1122 | AliITSneuron *fwd = 0, *back = 0; | |
853a0f19 | 1123 | TList *listsector = 0; |
b9d722bc | 1124 | |
1125 | cout << "A -" << fActMinimum << "-" << flush; | |
1126 | for (ilayer = 0; ilayer < 6; ilayer++) { | |
1127 | for (itheta = 0; itheta < 180; itheta++) { | |
853a0f19 | 1128 | listsector = (TList*)fPoints[ilayer][itheta]->At(sectoridx); |
1129 | TListIter iter(listsector); | |
b9d722bc | 1130 | while ( (node = (AliITSNode*)iter.Next()) ) { |
ddcce466 | 1131 | TListIter fwditer(node->InnerOf()); |
1132 | TListIter backiter(node->OuterOf()); | |
853a0f19 | 1133 | testfwd = testback = fActMinimum; |
1134 | while ( (fwd = (AliITSneuron*)fwditer.Next()) ) { | |
ddcce466 | 1135 | if (fwd->Activation() > testfwd) { |
1136 | testfwd = fwd->Activation(); | |
1137 | node->fNext = fwd->Outer(); | |
b9d722bc | 1138 | } |
1139 | } | |
853a0f19 | 1140 | while ( (back = (AliITSneuron*)backiter.Next()) ) { |
ddcce466 | 1141 | if (back->Activation() > testback) { |
1142 | testback = back->Activation(); | |
1143 | node->fPrev = back->Inner(); | |
b9d722bc | 1144 | } |
1145 | } | |
1146 | } | |
1147 | } | |
1148 | } | |
1149 | ||
1150 | cout << "B" << flush; | |
1151 | for (ilayer = 0; ilayer < 5; ilayer++) { | |
1152 | for (itheta = 0; itheta < 180; itheta++) { | |
853a0f19 | 1153 | listsector = (TList*)fPoints[ilayer][itheta]->At(sectoridx); |
1154 | TListIter iter(listsector); | |
b9d722bc | 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++) { | |
853a0f19 | 1169 | list[0] = (TList*)fPoints[0][itheta]->At(sectoridx); |
1170 | list[1] = (TList*)fPoints[1][itheta]->At(sectoridx); | |
b9d722bc | 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(); | |
853a0f19 | 1181 | for (im = 0; im < 7; im++) ok = ok || (mask == allowmask[im]); |
b9d722bc | 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 | // | |
1203 | Double_t AliITSNeuralTracker::Activate(AliITSneuron* &unit) | |
853a0f19 | 1204 | { |
b9d722bc | 1205 | // Computes the new neural activation and |
1206 | // returns the variation w.r.t the previous one | |
853a0f19 | 1207 | |
b9d722bc | 1208 | if (!unit) return 0.0; |
853a0f19 | 1209 | Double_t sumgain = 0.0, sumcost = 0.0, input, actOld, actNew; |
b9d722bc | 1210 | |
1211 | // sum gain contributions | |
ddcce466 | 1212 | TObjArrayIter *iter = (TObjArrayIter*)unit->Gain()->MakeIterator(); |
b9d722bc | 1213 | AliITSlink *link; |
1214 | for(;;) { | |
1215 | link = (AliITSlink*)iter->Next(); | |
1216 | if (!link) break; | |
ddcce466 | 1217 | sumgain += link->Linked()->Activation() * link->Weight(); |
b9d722bc | 1218 | } |
1219 | ||
1220 | // sum cost contributions | |
ddcce466 | 1221 | TObjArrayIter *testiter = (TObjArrayIter*)unit->Inner()->InnerOf()->MakeIterator(); |
b9d722bc | 1222 | AliITSneuron *linked = 0; |
1223 | for (;;) { | |
853a0f19 | 1224 | linked = (AliITSneuron*)testiter->Next(); |
b9d722bc | 1225 | if (!linked) break; |
1226 | if (linked == unit) continue; | |
ddcce466 | 1227 | sumcost += linked->Activation(); |
b9d722bc | 1228 | } |
853a0f19 | 1229 | delete testiter; |
ddcce466 | 1230 | testiter = (TObjArrayIter*)unit->Outer()->OuterOf()->MakeIterator(); |
b9d722bc | 1231 | for (;;) { |
853a0f19 | 1232 | linked = (AliITSneuron*)testiter->Next(); |
b9d722bc | 1233 | if (!linked) break; |
1234 | if (linked == unit) continue; | |
ddcce466 | 1235 | sumcost += linked->Activation(); |
b9d722bc | 1236 | } |
1237 | ||
853a0f19 | 1238 | //cout << "gain = " << sumgain << ", cost = " << sumcost << endl; |
b9d722bc | 1239 | |
853a0f19 | 1240 | input = (sumgain - sumcost) / fTemperature; |
ddcce466 | 1241 | actOld = unit->Activation(); |
b9d722bc | 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 | |
ddcce466 | 1252 | unit->Activation() = actNew; |
b9d722bc | 1253 | return TMath::Abs(actNew - actOld); |
1254 | } | |
1255 | // | |
1256 | // | |
1257 | // | |
853a0f19 | 1258 | void AliITSNeuralTracker::AliITSneuron::Add2Gain(AliITSneuron *n, Double_t multconst, Double_t exponent) |
b9d722bc | 1259 | { |
853a0f19 | 1260 | // Adds a neuron to the collection of sequenced ones of another neuron |
1261 | ||
b9d722bc | 1262 | AliITSlink *link = new AliITSlink; |
ddcce466 | 1263 | link->Linked() = n; |
b9d722bc | 1264 | Double_t weight = Weight(n); |
ddcce466 | 1265 | link->Weight() = multconst * TMath::Power(weight, exponent); |
b9d722bc | 1266 | fGain->AddLast(link); |
1267 | } | |
1268 | // | |
1269 | // | |
1270 | // | |
1271 | Double_t AliITSNeuralTracker::AliITSneuron::Weight(AliITSneuron *n) | |
1272 | { | |
853a0f19 | 1273 | // computes synaptic weight |
ddcce466 | 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()); | |
b9d722bc | 1276 | |
1277 | Double_t angle = me.Angle(it); | |
1278 | Double_t weight = 1.0 - sin(angle); | |
1279 | return weight; | |
1280 | } | |
ddcce466 | 1281 | // |
1282 | // | |
1283 | // | |
1284 | // | |
cb729508 | 1285 | |
1286 | ||
1287 | AliITSNeuralTracker::AliITSNode::AliITSNode(const AliITSNode &t): AliITSNeuralPoint((AliITSNeuralPoint&)t), | |
1288 | fPosInTree(t.fPosInTree), | |
1289 | fInnerOf(t.fInnerOf), | |
1290 | fOuterOf(t.fOuterOf), | |
1291 | fMatches(t.fMatches), | |
1292 | fNext(t.fNext), | |
1293 | fPrev(t.fPrev) | |
ddcce466 | 1294 | { |
cb729508 | 1295 | //copy constructor |
ddcce466 | 1296 | } |
cb729508 | 1297 | |
1298 | AliITSNeuralTracker::AliITSNode& AliITSNeuralTracker::AliITSNode::operator=(const AliITSNode& t) | |
ddcce466 | 1299 | { |
cb729508 | 1300 | // assignment operator |
1301 | this->~AliITSNode(); | |
1302 | new(this) AliITSNode(t); | |
1303 | return *this; | |
ddcce466 | 1304 | } |
cb729508 | 1305 | |
1306 | ||
ddcce466 | 1307 | AliITSNeuralTracker::AliITSneuron::AliITSneuron(const AliITSneuron &t) |
cb729508 | 1308 | : TObject((TObject&)t), |
1309 | fUsed(t.fUsed), | |
1310 | fActivation(t.fActivation), | |
1311 | fInner(t.fInner), | |
1312 | fOuter(t.fOuter), | |
1313 | fGain(t.fGain) | |
ddcce466 | 1314 | { |
cb729508 | 1315 | //copy constructor |
ddcce466 | 1316 | } |
cb729508 | 1317 | |
1318 | AliITSNeuralTracker::AliITSneuron& AliITSNeuralTracker::AliITSneuron::operator=(const AliITSneuron& t) | |
ddcce466 | 1319 | { |
cb729508 | 1320 | //assignment operator |
1321 | this->~AliITSneuron(); | |
1322 | new(this) AliITSneuron(t); | |
1323 | return *this; | |
ddcce466 | 1324 | } |
cb729508 | 1325 | |
1326 | ||
ddcce466 | 1327 | AliITSNeuralTracker::AliITSlink::AliITSlink(const AliITSlink &t) |
cb729508 | 1328 | : TObject((TObject&)t), |
1329 | fWeight(t.fWeight), | |
1330 | fLinked(t.fLinked) | |
ddcce466 | 1331 | { |
cb729508 | 1332 | //copy constructor |
ddcce466 | 1333 | } |
cb729508 | 1334 | |
1335 | AliITSNeuralTracker::AliITSlink& AliITSNeuralTracker::AliITSlink::operator=(const AliITSlink& t) | |
ddcce466 | 1336 | { |
cb729508 | 1337 | // assignment operator |
1338 | this->~AliITSlink(); | |
1339 | new(this) AliITSlink(t); | |
1340 | return *this; | |
ddcce466 | 1341 | } |
1342 | ||
1343 | ||
1344 |