1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: Alberto Pulvirenti. *
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. It is *
12 * provided "as is" without express or implied warranty. *
14 * track finder for ITS stand-alone with neural network algorithm *
15 * This class defines a Hopfield MFT neural network simulation *
16 * which reads all recpoints from an event and produces a tree with *
17 * the points belonging to recognized tracks *
18 * the TTree obtained as the output file will be saved in a .root file *
19 * and it must be read by groups of six entries at a time *
20 **************************************************************************/
28 #include <TObjArray.h>
36 #include "AliITSgeom.h"
37 #include "AliITSgeomMatrix.h"
38 #include "AliITSRecPoint.h"
39 #include "AliITSglobalRecPoint.h"
40 #include "AliITSneuralTrack.h"
41 #include "AliITSneuralTracker.h"
43 ClassImp(AliITSneuralTracker)
46 // ====================================================================================================
47 // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> AliITSneuralTracker <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
48 // ====================================================================================================
50 AliITSneuralTracker::AliITSneuralTracker()
53 // Default constructor (not for use)
57 for (i = 0; i < 6; i++) {
59 if (i < 5) fNeurons[i] = 0;
63 //---------------------------------------------------------------------------------------------------------
64 AliITSneuralTracker::AliITSneuralTracker (Int_t nsecs, Int_t ncurv, Double_t *curv, Double_t theta)
66 // Argumented constructor
67 // ----------------------
68 // gets the number of azymuthal sectors,
69 // the number of curvature cut steps and their cut values,
70 // and then the polar angle cut.
71 // Be careful not to put 'ncurv' greater than the dimension of the
72 // 'curv' array (--> segmentation fault)
75 fSectorWidth = 2.0 * TMath::Pi()/(Double_t)nsecs;
78 fCurvCut = new Double_t[ncurv];
80 for (i = 0; i < ncurv; i++) fCurvCut[i] = curv[i];
83 fThetaNum = (Int_t)(TMath::Pi() / theta) + 1;
84 if (fThetaNum < 1) fThetaNum = 1;
86 cout << "\nClass AliITSneuralTracker invoked with " << fSectorNum << " phi secs.";
87 cout << " and " << fThetaNum << " theta secs.\n" << endl;
90 for (k = 0; k < 6; k++) {
91 fPoints[k] = new TObjArray**[fSectorNum];
92 for (i = 0; i < fSectorNum; i++) {
93 fPoints[k][i] = new TObjArray*[fThetaNum];
94 for (j = 0; j < fThetaNum; j++) fPoints[k][i][j] = new TObjArray(0);
96 if (k < 6) fNeurons[k] = new TObjArray(0);
99 fTracks = new TObjArray(1);
101 //---------------------------------------------------------------------------------------------------------
102 AliITSneuralTracker::~AliITSneuralTracker()
106 // clears the TObjArrays and all heap objects
111 for (lay = 0; lay < 6; lay++) {
112 for (sec = 0; sec < fSectorNum; sec++) {
113 for (k = 0; k < fThetaNum; k++) {
114 fPoints[lay][sec][k]->Delete();
117 delete [] fPoints[lay];
119 fNeurons[lay]->Delete();
120 delete fNeurons[lay];
126 //---------------------------------------------------------------------------------------------------------
127 Int_t AliITSneuralTracker::ReadFile(const Text_t *fname, Int_t evnum)
131 // Reads a ROOT file in order to retrieve recpoints.
132 // the 'evnum' argument can have two meanings:
133 // if it is >= 0 it represents the event number to retrieve the correct gAlice
134 // if it is < 0, instead, this means that the open file contains
135 // a 'TreeP' tree, with all points which remained unused
136 // after the Kalman tracking procedure (produced via the ITSafterKalman.C macro)
137 // the return value is the number of collected points
138 // (then, if errors occur, the method returns 0)
140 TFile *file = new TFile(fname, "read");
142 Int_t i, nentries, total, sector, mesh;
145 tree = (TTree*)file->Get("TreeP0");
147 Error("", "Specifying evnum < 0, I expect to find a 'TreeP' tree in the file, but there isn't any");
150 AliITSglobalRecPoint *p = 0;
151 tree->SetBranchAddress("Points", &p);
152 nentries = (Int_t)tree->GetEntries();
154 for (i = 0; i < nentries; i++) {
156 AliITSglobalRecPoint *q = new AliITSglobalRecPoint(*p);
157 sector = (Int_t)(q->fPhi / fSectorWidth);
158 mesh = (Int_t)(q->fTheta / fThetaCut);
160 fPoints[q->fLayer][sector][mesh]->AddLast(q);
166 if (gAlice) { delete gAlice; gAlice = 0; }
167 gAlice = (AliRun*)file->Get("gAlice");
169 Error("", "gAlice is NULL, maybe wrong filename...");
173 Int_t nparticles = (Int_t)gAlice->GetEvent(evnum);
175 Error("", "No particles found");
179 AliITS *its = (AliITS*)gAlice->GetModule("ITS");
181 Error("", "No ITS found");
185 AliITSgeom *geom = (AliITSgeom*)its->GetITSgeom();
187 Error("", "AliITSgeom is NULL");
191 tree = gAlice->TreeR();
193 Error("", "TreeR() returned NULL");
197 nentries = (Int_t)tree->GetEntries();
199 Error("", "TreeR is empty");
203 // Reading objects initialization
204 Int_t k, entry, lay, lad, det;
205 Double_t l[3], g[3], ll[3][3], gg[3][3];
207 TObjArray *recsArray = its->RecPoints();
208 AliITSRecPoint *recp;
209 AliITSglobalRecPoint *pnt = 0;
212 for(entry = 0; entry < nentries; entry++) {
213 if (entry > geom->GetLastSSD()) continue;
214 its->ResetRecPoints();
215 tree->GetEvent(entry);
216 TObjArrayIter recs(recsArray);
217 AliITSgeomMatrix *gm = geom->GetGeomMatrix(entry);
218 geom->GetModuleId(entry, lay, lad, det);
220 while ((recp = (AliITSRecPoint*)recs.Next())) {
221 // conversion to global coordinate system
222 for (i = 0; i < 3; i++) {
224 for (k = 0; k < 3; k++) {
229 // local to global conversions of coords
233 gm->LtoGPosition(l, g);
234 // local to global conversions of sigmas
235 ll[0][0] = recp->fSigmaX2;
236 ll[2][2] = recp->fSigmaZ2;
237 gm->LtoGPositionError(ll, gg);
238 // instantiation of a new global rec-point
239 pnt = new AliITSglobalRecPoint(g[0], g[1], g[2], gg[0][0], gg[1][1], gg[2][2], lay);
240 // copy of owner track labels
241 for (i = 0; i < 3; i++) pnt->fLabel[i] = recp->fTracks[i];
242 sector = (Int_t)(pnt->fPhi / fSectorWidth);
243 mesh = (Int_t)(pnt->fTheta / fThetaCut);
245 fPoints[lay][sector][mesh]->AddLast(pnt);
252 //---------------------------------------------------------------------------------------------------------
253 void AliITSneuralTracker::Go(const char* filesave, Bool_t flag)
255 // Global working method
256 // ---------------------
257 // It's the method which calls all other ones,
258 // in order to analyze the whole ITS points ensemble
259 // It's the only method to use, besides the parameter setters,
260 // so, all other methods are defined as private
262 // the flag meaning is explained in the 'Initialize' method
265 for (i = 0; i < fSectorNum; i++) DoSector(i, flag);
266 cout << endl << "Saving results in file " << filesave << "..." << flush;
267 TFile *f = new TFile(filesave, "recreate");
269 cout << "done" << endl;
273 //=========================================================================================================
274 // * * * P R I V A T E S E C T I O N * * *
275 //=========================================================================================================
277 Int_t AliITSneuralTracker::Initialize(Int_t secnum, Int_t curvnum, Bool_t flag)
279 // Network initializer
280 // -------------------
281 // Fills the neuron arrays, following the cut criteria for the selected step
282 // (secnum = sector to analyze, curvnum = curvature cut step to use)
283 // It also sets the initial random activation.
284 // In order to avoid a large number of 'new' operations, all existing neurons
285 // are reset and initialized with the new values, and are created new unit only if
286 // it turns out to be necessary
287 // the 'flag' argument is used to decide if the lower edge in the intevral
288 // of the accepted curvatures is given by zero (kFALSE) or by the preceding used cut (kTRUE)
289 // (if this is the first step, anyway, the minimum is always zero)
291 Int_t l, m, k, neurons = 0;
293 for (l = 0; l < 5; l++) fNeurons[l]->Delete();
295 AliITSneuron *unit = 0;
296 Double_t abscurv, max = fCurvCut[curvnum], min = 0.0;
297 if (flag && curvnum > 0) min = fCurvCut[curvnum - 1];
298 AliITSglobalRecPoint *inner = 0, *outer = 0;
299 for (l = 0; l < 5; l++) {
300 for (m = 0; m < fThetaNum; m++) {
301 TObjArrayIter inners(fPoints[l][secnum][m]);
302 while ((inner = (AliITSglobalRecPoint*)inners.Next())) {
303 if (inner->fUsed > 0) continue; // points can't be recycled
304 for (k = m-1; k <= m+1; k++) {
305 if (k < 0 || k >= fThetaNum) continue; // to avoid seg faults
306 TObjArrayIter outers(fPoints[l+1][secnum][k]);
307 while ((outer = (AliITSglobalRecPoint*)outers.Next())) {
308 if (outer->fUsed > 0) continue; // points can't be recycled
309 if (inner->DTheta(outer) > fThetaCut) continue;
310 unit = new AliITSneuron;
311 unit->fInner = inner;
312 unit->fOuter = outer;
314 abscurv = TMath::Abs(unit->fCurv);
315 if (unit->fDiff > fDiff || abscurv < min || abscurv > max) {
319 unit->fActivation = gRandom->Rndm() * (fMax - fMin) + fMin;
320 fNeurons[l]->AddLast(unit);
322 } // end loop on candidates for outer point for neurons
324 } // end loop on candidates inner points for neuron
326 } // end loop on layers
329 //---------------------------------------------------------------------------------------------------------
330 Bool_t AliITSneuralTracker::Update()
334 // Performs an updating cycle, by summing all the
335 // gain and cost contribution for each neuron and
336 // updating its activation.
337 // An asynchronous method is followed, with the order
338 // according that the neurons have been created
339 // Time wasting is avoided by dividing the neurons
340 // into different arrays, depending on the layer of their inner point.
341 // The return value answers the question: "The network has stabilized?"
344 Double_t actOld, actNew, argNew, totDiff = 0.0;
345 Double_t sumInh, sumExc, units = 0.0;
346 AliITSneuron *me = 0, *it = 0;
347 for (l = 0; l < 5; l++) {
348 TObjArrayIter meIter(fNeurons[l]);
349 while ((me = (AliITSneuron*)meIter.Next())) {
351 TObjArrayIter inhIter(fNeurons[l]);
353 // inhibition (search only for neurons starting in the same layer)
355 while((it = (AliITSneuron*)inhIter.Next())) {
356 if (it->fOuter == me->fOuter || it->fInner == me->fInner)
357 sumInh += it->fActivation;
360 // excitation (search only for neurons
361 // which start in the previous or next layers)
363 for (j = l - 1; j <= l + 1; j += 2) {
364 if (j < 0 || j >= 5) continue;
365 TObjArrayIter itIter(fNeurons[j]);
366 while ((it = (AliITSneuron*)itIter.Next())) {
367 if (it->fInner == me->fOuter || it->fOuter == me->fInner && it->fCurv * me->fCurv > 0.0) {
368 sumExc += Weight(me, it) * it->fActivation;
371 } // end search for excitories
373 actOld = me->fActivation;
374 argNew = fRatio * sumExc - sumInh;
375 actNew = 1.0 / (1.0 + TMath::Exp(-argNew / fTemp));
376 me->fActivation = actNew;
377 // calculation the relative activation variation
378 // (for stabilization check)
379 totDiff += TMath::Abs((actNew - actOld) / actOld);
380 } // end loop over 'me' (updated neuron)
381 } // end loop over layers
384 return (totDiff < fVar);
386 //---------------------------------------------------------------------------------------------------------
387 Int_t AliITSneuralTracker::Save()
389 // Tracki saving method
390 // --------------------
391 // Stores the recognized tracks into the TObjArray 'fTracks'
392 // but doesn't save it into a file until the whole work has ended
397 // Before all, for each neuron is found the best active sequences
398 // among the incoming and outgoing other units
399 AliITSneuron *main, *fwd, *back, *start;
400 for (l = 0; l < 5; l++) {
401 TObjArrayIter mainIter(fNeurons[l]);
402 while((main = (AliITSneuron*)mainIter.Next())) {
404 TObjArrayIter fwdIter(fNeurons[l+1]);
406 while ((fwd = (AliITSneuron*)fwdIter.Next())) {
407 if (main->fOuter == fwd->fInner && fwd->fActivation > test) {
408 test = fwd->fActivation;
414 TObjArrayIter backIter(fNeurons[l-1]);
416 while ((back = (AliITSneuron*)backIter.Next())) {
417 if (main->fInner == back->fOuter && back->fActivation > test) {
418 test = back->fActivation;
426 // Then, are resolved the mismatches in the chains found above
427 // (when unit->next->prev != unit or unit->prev->next != unit)
428 for (l = 0; l < 5; l++) {
429 TObjArrayIter mainIter(fNeurons[l]);
430 while((main = (AliITSneuron*)mainIter.Next())) {
433 if (fwd != 0 && fwd->fPrev != main) main->fNext = 0;
434 if (back != 0 && back->fNext != main) main->fPrev = 0;
438 Int_t pos, neg, incr, decr, ntracks;
439 AliITSneuralTrack *trk = 0;
441 TObjArrayIter startIter(fNeurons[0]);
443 start = (AliITSneuron*)startIter.Next();
446 // with the chain above defined, a track can be followed
447 // like a linked list structure
448 for (main = start; main; main = main->fNext) pts++;
449 // the count will be > 5 only if the track contains 6 points
451 if (pts < 5) continue;
453 trk = new AliITSneuralTrack;
454 trk->CopyPoint(start->fInner);
455 for (main = start; main; main = main->fNext) {
456 //main->fInner->fUsed = kTRUE;
457 //main->fOuter->fUsed = kTRUE;
458 trk->CopyPoint(main->fOuter);
460 trk->Kinks(pos, neg, incr, decr);
461 if (pos != 0 && neg != 0) {
462 cout << "pos, neg kinks = " << pos << " " << neg << endl;
465 if (incr != 0 && decr != 0) {
466 cout << "increment, decrements in delta phi = " << incr << " " << decr << endl;
469 for (main = start; main; main = main->fNext) {
470 main->fInner->fUsed++;
471 main->fOuter->fUsed++;
473 fTracks->AddLast(trk);
475 ntracks = (Int_t)fTracks->GetEntriesFast();
478 //---------------------------------------------------------------------------------------------------------
479 void AliITSneuralTracker::DoSector(Int_t sect, Bool_t flag)
481 // Sector recognition
482 // ------------------
483 // This is a private method which works on a single sector
484 // just for an easier readability of the code.
485 // The 'flag' argument is needed to be given to the 'Initialize' method (see it)
487 Int_t curvnum, nTracks = 0, nUnits;
488 Bool_t isStable = kFALSE;
489 for (curvnum = 0; curvnum < fCurvNum; curvnum++) {
490 cout << "\rSector " << sect << ":" << ((curvnum+1)*100)/fCurvNum << "%" << flush;
491 nUnits = Initialize(sect, curvnum, flag);
492 if (!nUnits) continue;
498 cout << "\rTracks stored after sector " << sect << ": " << nTracks << endl;
500 //---------------------------------------------------------------------------------------------------------
501 Int_t AliITSneuralTracker::CountExits(AliITSneuron *n)
503 // Counter for neurons which enter in 'n' tail
505 Int_t count = 0, l = n->fOuter->fLayer;
506 if (l == 5) return 0;
508 AliITSneuron *test = 0;
509 TObjArrayIter iter(fNeurons[l]);
510 while ((test = (AliITSneuron*)iter.Next())) {
511 if (test->fInner == n->fOuter) count++;
515 //---------------------------------------------------------------------------------------------------------
516 Int_t AliITSneuralTracker::CountEnters(AliITSneuron *n)
518 // Counter for neurons which exit from 'n' head
520 Int_t count = 0, l = n->fInner->fLayer;
521 if (l == 0) return 0;
523 AliITSneuron *test = 0;
524 TObjArrayIter iter(fNeurons[l-1]);
525 while ((test = (AliITSneuron*)iter.Next())) {
526 if (test->fOuter == n->fInner) count++;
530 //---------------------------------------------------------------------------------------------------------
531 Bool_t AliITSneuralTracker::ResetNeuron(AliITSneuron *n)
533 // A method which resets the neuron's pointers to zero
535 if (!n) return kFALSE;
540 // the other datamembers don't need to be reinitialized
543 //---------------------------------------------------------------------------------------------------------
544 Bool_t AliITSneuralTracker::SetEdge(AliITSneuron *n, AliITSglobalRecPoint *p, const char what)
546 // Sets the indicated edge point to the neuron
547 // (if the arguments are suitable ones...)
549 if (!n || !p || (what != 'i' && what != 'I' && what != 'o' && what != 'O')) return kFALSE;
550 if (what == 'i' || what == 'I')
556 //---------------------------------------------------------------------------------------------------------
557 Bool_t AliITSneuralTracker::CalcParams(AliITSneuron *n)
559 // This method evaluates the helix parameters from the edged of a neuron
560 // (curvature, phase parameter, tangent of lambda angle)
561 // if the p[arameters assume too strange (and unphysical) values, the
562 // method return kFALSE, else it return kTRUE;
564 if (!n) return kFALSE;
566 Double_t sign(1.0), det, rc = 0.0;
567 Double_t tanl1, tanl2, l1, l2;
568 // changing the reference frame into the one centered in the vertex
569 Double_t x1 = n->fInner->fGX - fVPos[0];
570 Double_t y1 = n->fInner->fGY - fVPos[1];
571 Double_t z1 = n->fInner->fGZ - fVPos[2];
572 Double_t r1 = x1 * x1 + y1 * y1;
573 Double_t x2 = n->fOuter->fGX - fVPos[0];
574 Double_t y2 = n->fOuter->fGY - fVPos[1];
575 Double_t z2 = n->fOuter->fGZ - fVPos[2];
576 Double_t r2 = x2 * x2 + y2 * y2;
577 // initialization of the XY-plane data (curvature and centre)
578 // (will remain these values if we encounter errors)
579 n->fCurv = n->fCX = n->fCY = n->fTanL = n->fPhase = 0.0;
580 n->fLength = (n->fOuter->fGX - n->fInner->fGX) * (n->fOuter->fGX - n->fInner->fGX);
581 n->fLength += (n->fOuter->fGY - n->fInner->fGY) * (n->fOuter->fGY - n->fInner->fGY);
582 n->fLength += (n->fOuter->fGZ - n->fInner->fGZ) * (n->fOuter->fGZ - n->fInner->fGZ);
583 n->fLength = TMath::Sqrt(n->fLength);
584 // this determinant which is zero when the points are in a line
585 det = 2.0 * (x1 * y2 - x2 * y1);
587 return kFALSE; // it is difficult having perfectly aligned points --- it's more likely to be an error
589 n->fCX = (r1*y2 - r2*y1) / det;
590 n->fCY = (r2*x1 - r1*x2) / det;
591 rc = TMath::Sqrt(n->fCX * n->fCX + n->fCY * n->fCY);
592 r1 = TMath::Sqrt(r1);
593 r2 = TMath::Sqrt(r2);
594 sign = (n->fOuter->fPhi >= n->fInner->fPhi) ? 1.0 : -1.0;
596 n->fCurv = sign / rc;
597 n->fPhase = TMath::ATan2(-n->fCX, -n->fCY);
598 if (n->fPhase < 0.0) n->fPhase += 2.0 * TMath::Pi(); // angles in 0 --> 2pi
599 if (r1 <= 2.0 * rc && r2 <= 2.0 * rc) {
600 l1 = 2.0 * rc * TMath::ASin(r1 / (2.0 * rc));
601 l2 = 2.0 * rc * TMath::ASin(r2 / (2.0 * rc));
604 n->fDiff = TMath::Abs(tanl1 - tanl2);
605 n->fTanL = 0.5 * (tanl1 + tanl2);
610 return kFALSE; // it' even more difficult that the arguments above aren't acceptable for arcsine
618 //---------------------------------------------------------------------------------------------------------
619 Double_t AliITSneuralTracker::Angle(AliITSneuron *n, AliITSneuron *m)
621 // calculates the angle between two segments
622 // but doesn't check if they have a common point.
623 // The calculation is made by the ratio of their scalar product
624 // to the product of their magnitudes
626 Double_t x = (m->fOuter->fGX - m->fInner->fGX) * (n->fOuter->fGX - n->fInner->fGX);
627 Double_t y = (m->fOuter->fGY - m->fInner->fGY) * (n->fOuter->fGY - n->fInner->fGY);
628 Double_t z = (m->fOuter->fGZ - m->fInner->fGZ) * (n->fOuter->fGZ - n->fInner->fGZ);
630 Double_t cosine = (x + y + z) / (n->fLength * m->fLength);
632 if (cosine > 1.0 || cosine < -1.0) {
633 Warning("AliITSneuronV1::Angle", "Strange value of cosine");
637 return TMath::ACos(cosine);
639 //---------------------------------------------------------------------------------------------------------
640 Double_t AliITSneuralTracker::Weight(AliITSneuron *n, AliITSneuron *m)
642 // calculation of the excitoy weight
644 // Double_t v1 = TMath::Abs(fTanL - n->fTanL);
645 // Double_t v2 = TMath::Abs(TMath::Abs(fCurv) - TMath::Abs(n->fCurv));
646 // Double_t v3 = TMath::Abs(fPhase - n->fPhase);
647 // Double_t q = (v1 + v2 + v3) / 3.0;
648 // Double_t a = 1.0; // 0.02 / (fDiff + n->fDiff);
649 Double_t b = 1.0 - TMath::Sin(Angle(n, m));
650 return TMath::Power(b, fExp);
651 // return TMath::Exp(q / par1);
654 //---------------------------------------------------------------------------------------------------------
655 // ====================================================================================================
656 // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> AliITSneuron <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
657 // ====================================================================================================
659 AliITSneuron::AliITSneuron()
661 // Default constructor
662 // -------------------
663 // sets all parameters to default (but unuseful)
664 // (in particular all pointers are set to zero)
665 // There is no alternative constructor, because
666 // this class should not be used by itself, but
667 // it is only necessary to allow the neural tracker to work
668 // In fact, all methods which operate on neuron's datamembers
669 // are stored in the AliITSneuralTracker class anyway
684 //---------------------------------------------------------------------------------------------------------
685 AliITSneuron::~AliITSneuron()
689 // does nothing, because there is no need to affect parameters,
690 // while the AliITSglobalRecPoint pointers can't be deleted,
691 // and the AliITSneuron pointers belong to a TObjArray which will
692 // be cleared when necessary...